void emc_reco(const char* digiFile, const char* simFile, const char* outFile, Int_t events=0, const char* parFile="simparams_hc7g5kG4.root", Double_t dtau=1000., Double_t dt=300., Bool_t runOnlineClus = kTRUE, Bool_t runDistrClus = kTRUE, Bool_t useTimebasedSim = kTRUE) {
	// Loads files with hits and digitised hits and performs reconstruction for EMC

// PARAMETERS:
	// dtau: value used by gapfinder for timebunch cutting. The value specified here will be used if SetAutoDetermineTimecuts is set to kFALSE
	// dt: value used by cluster task to determine which digis are temporal neighbours. The value specified here will be used if SetAutoDetermineTimecuts is set to kFALSE
	// Set runOnlineClus to kFALSE to use default clustering algorithms (will use online version if set to kTRUE)
	// Set runDistrClus to kTRUE to enable Distributed Cluster Finding (runOnlineClus must be set to kTRUE first)
	// >>Set useTimebasedSim to kTRUE when using timebased sim, kFALSE for eventbased sim (DON'T FORGET TO SET THIS)

	// Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
	Int_t iVerbose = 0; 

	// Input file (MC events)
	//
	// Number of events to process
	Int_t nEvents=events;  // if 0 all the events will be processed
	if (nEvents !=0) std::cout << "will process " << nEvents << " events" << std::endl;
	else std::cout << "all events will be processed" << std::endl;

	// Digitisation file (ascii)
	TString parFile_Ascii = "emc.par";

	// Loading libraries
	// If the macro gives error messages in loading libraries, please check the path of the libs and put it by hands

	//  gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C");
	// gROOT->LoadMacro("$VMCWORKDIR/gconfig/basiclibs.C");
	//rootlogon();
	//basiclibs();


	// -----   Timer   -----------------------------------------
	TStopwatch timer;
	timer.Start();
	// ---------------------------------------------------------

	// -----   Reconstruction run   ----------------------------
	FairRunAna *fRun= new FairRunAna();
	fRun->SetInputFile(digiFile);
	fRun->AddFriend(simFile);
	fRun->SetOutputFile(outFile);

	fRun->SetUseFairLinks(kTRUE);

	// -----  Parameter database   -----------------------------
	FairRuntimeDb* rtdb = fRun->GetRuntimeDb();

	FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo();
	parIo1->open(parFile_Ascii.Data(),"in");
	rtdb->setFirstInput(parIo1);

	FairParRootFileIo* parInput1 = new FairParRootFileIo();
	parInput1->open(parFile);
	rtdb->setSecondInput(parInput1);

	// ---------------------------------------------------------

	Bool_t storeclusters = kTRUE;
	Bool_t storebumps = kTRUE;
//	Bool_t runOnlineClus = kTRUE; // set to kTRUE to enable the online cluster finding algorithms, or to kFALSE to use the default one
//	Bool_t runDistrClus = kTRUE; // set to kTRUE to enable Distributed Cluster Finding (runOnlineClus must be set to kTRUE first)

//	Bool_t useTimebasedSim = kTRUE; // set to kTRUE when using timebased sim, kFALSE for eventbased sim (DON'T FORGET TO SET THIS)

	if (useTimebasedSim) fRun->RunWithTimeStamps();

	if (runOnlineClus) {
		if (runDistrClus) { // Use Distributed Clustering
			PndEmcDistributedClustering* clusterTask = new PndEmcDistributedClustering(iVerbose, storeclusters);
			clusterTask->SetPositionMethod(1); // set (pre)cluster position and radius method: 
//														0 = default method from the parfile (by default, logarithmic energy weighing, "lilo"), 
//														1 = simplified method for preclusters only; default method for clusters (RECOMMENDED), 
//														2 = simplified method for preclusters only; position of digi with highest energy for clusters
			clusterTask->SetNeighbourMethod(0); // Set cluster neighbour method: 
//					0 = default method (for position, logarithmic weighing; for radius, distance of digi furthest from position) -> use this CIRCLE for neighbour relations),  
//					1 = simplified method (for position, xpos=(xmax-xmin)/2 ypos=(ymax-ymin)/2; for radius, r=max(ysize, xsize) -> use this CIRCLE for neighbour relations), 
//					2 = simplified method (for position, xpos=(xmax+xmin)/2 ypos=(ymax+ymin)/2; for radius, rx=(xmax-xmin)/2, ry=(ymax-ymin)/2 -> use this RECTANGULAR BOX for neighbour relations). 
//					3 = simplified method (for position, xpos=(xmax+xmin)/2 ypos=(ymax+ymin)/2; for radius, r=max(ysize, xsize) -> use this SQUARE BOX for neighbour relations). 
		}
		else { // Use Online Clustering
			PndEmcMakeClusterOnline* clusterTask = new PndEmcMakeClusterOnline(iVerbose, storeclusters);
		}
	}
	else { // Use Default PandaRoot Clustering
		PndEmcMakeCluster* clusterTask = new PndEmcMakeCluster(iVerbose, storeclusters);	
	}

//--- Set some common cluster finding parameters ---//
	clusterTask->StoreClusterBaseDigis(kTRUE);
	clusterTask->SetTimebunchCutTime(dtau); // change timebunch cutting time (minimal time between "events")
	//clusterTask->SetClusterMinimumEnergy(0.040); // change minimum cluster energy here (set to 0.03 GeV by default)
	//clusterTask->EnableRemovalOfLowEnergyClusters(kFALSE); // enable/disable removal of low energy (i.e. below the minimum energy) clusters (enabled by default). Slows down cluster finding task, but speeds up future processing, such as analysis.
	if (runOnlineClus) clusterTask->SetClusterActiveTime(dt); // change threshold for considering digis to belong to the same cluster (only for online cluster finding tasks)
	if (useTimebasedSim) clusterTask->SetAutoDetermineTimecuts(kTRUE); // automatically determine time cuts. Will override values set by SetTimebunchCutTime and SetClusterActiveTime
	fRun->AddTask(clusterTask);	
//----   ----//

	if (useTimebasedSim) {
		// Sorting (only for timebased sim)
		PndEmcClusterSorterTask* sorterTask = new PndEmcClusterSorterTask(10000, 1, "EmcClusterTemp", "EmcClusterSorted", "Emc");
		//sorterTask->SetVerbose(iVerbose); //												^					^
		fRun->AddTask(sorterTask);//														|					|
//																						input branch		output branch
	}

	//PndEmcMakeBump* bumperTask = new PndEmcMakeBump(iVerbose, storebumps); // bump splitting
	//fRun->AddTask(bumperTask);

	// -----   Intialise and run   -----------------------------

	gRandom->SetSeed();
	cout << "fRun->Init()" << endl;
	fRun->Init();

	cout << "fRun->Run()" << endl;
	fRun->Run(0,nEvents);

	// ---------------------------------------------------------


	// -----   Finish   ----------------------------------------
	timer.Stop();
	Double_t rtime = timer.RealTime();
	Double_t ctime = timer.CpuTime();
	cout << endl << endl;
	cout << "Macro finished successfully." << endl;
	cout << "Output file is "    << outFile << endl;
	cout << "Parameter file is " << parFile << endl;
	cout << "Real time " << rtime << " s, CPU time " << ctime << " s" << endl;
	cout << endl;
	// ---------------------------------------------------------
  exit(0);
}