void runReco_mb(unsigned int run,
		TString outpath,
		TString outFileName = "NULL",
		bool clcorr=false,
		bool reclust=false,
		unsigned int nEvents = 0,
		unsigned int nEvStart = 0,
		unsigned int bfield=1,
		unsigned int smoothing = 1,
		int outnum=-1,
		unsigned int gaincorr=1,
		double fac = 0.0,
		double stepoff= 1.0,
		double fC=14, //this is just a linear scaling of the ols x,y,z independent error
		double fG=-1, //-1 std error; -2 unscaled shape cov; <0 scaled shape cov; >0 covOnPlane
		Int_t scaler=99, //this controls the scaling of the cov based errors 55=covOnPlane; 56=covOnPlane scaled
		bool fopi=false,
		bool doPRF=true
)
{
	// ========================================================================
	// Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug)
	Int_t iVerbose = 0;
	Int_t persistence=0;


	TString zCorrFile = "datafiles/zShifts.root";
	TStopwatch timer;
	timer.Start();
	TString filename="";
	TString fopiFilename="";
	TString outfilename="";
	if (outFileName != "NULL")
	{
		outfilename += outFileName;
		outfilename += ".root";

		//TString filename="CosmicData_e12/";
		filename+= gSystem->Getenv("datapath");
		if(fopi)
		{
			filename+="/run_";
			filename+=run;
			filename+=".lmd_decoded.root";

			fopiFilename=gSystem->Getenv("fopipath");
			fopiFilename+="/FOPI_";
			fopiFilename+=run;
			fopiFilename+="_t0_0";
			fopiFilename+=".root";
		}
		else
		{
			//filename+="/runC_";
			filename+="/run_";
			filename+=run;
			filename+=".lmd_decoded.root";
		}
	}
	else
	{
		filename+="decoded/";
		filename+="runC_";
		filename+=run;
		filename+=".lmd_decoded.root";
	}

	std::cout<<"using files:\n"<<filename<<"\n"<<fopiFilename<<"\n";
	TString basedir = gSystem->Getenv("VMCWORKDIR");
	FairRunAna* fRun = new FairRunAna();

	TString jobdir = outpath;
	std::string jobname(filename.Data());
	int last = jobname.rfind("/");
	if(last>0)
		jobname = jobname.substr(last+1,jobname.size()+1);

	//create output file -----------------------------------------------
	TString outName(jobname);
	if(outName.Contains("repaired."))
		outName.ReplaceAll(".lmd_decoded_repaired","");
	if(outName.Contains("decoded."))
		outName.ReplaceAll(".lmd_decoded","");
	if(smoothing)
		outName.ReplaceAll(".root", "_smoothed.root");
	if(clcorr)
		outName.ReplaceAll(".root","_clcorr.root");
	if (reclust)
		outName.ReplaceAll(".root","_recl.root");
	if(outnum>-1)
		outName.ReplaceAll(".root",Form("part_%03i.root",outnum));


	TString outFile = outpath+"/";
	if (outFileName != "NULL")
	{
		std::cout<<outFile<<outfilename<<"\n";
		outFile += outfilename;
	}
	else
	{
		outFile += outName;
	}

	outFile.ReplaceAll(".root",".reco.root");
	std::cout<<"OUTFILE:"<<outFile<<"\n";
	fRun->SetOutputFile(outFile);
	//DONE create output file -----------------------------------------------

	FairRuntimeDb* rtdb = fRun->GetRuntimeDb();
	FairParAsciiFileIo* parInput1 = new FairParAsciiFileIo();
	TString parFile = basedir;
	parFile+="/tpc/FOPI/par/tpc.run";
	parFile+=run;
	parFile +=".par";
	std::cout<<parFile<<std::endl;
	std::fstream fin;
	fin.open(parFile.Data(),std::ios::in);
	if( !fin.is_open() )
	{
		std::cout<<"parfile doesnt exist"<<std::endl;
		return;
	}
	fin.close();
	parInput1->open(parFile.Data(),"in");
	rtdb->setFirstInput(parInput1);

	FairParRootFileIo* parOut=new FairParRootFileIo(kTRUE);
	TString parfileOut(outFile);
	parfileOut.ReplaceAll("reco.root","param.root");
	if (outnum<=0)
	{
		parOut->open(parfileOut.Data());
		rtdb->setOutput(parOut);
	}
	//TString parOutFile=outFile;
	//parOutFile.ReplaceAll("reco.root","param.root");
	//FairParRootFileIo* parOut=new FairParRootFileIo(kTRUE);
	//parOut->open(parOutFile.Data());
	//rtdb->setOutput(parOut);

	TpcDigiPar* par = (TpcDigiPar*) rtdb->getContainer("TpcDigiPar");
	par->setInputVersion(fRun->GetRunId(),1);
	par->setChanged(kTRUE);
	//rtdb->Print();
	par->init();


	Double_t taroff = -61;
	Double_t targetpos = taroff+41.;
	// possible materials: copper (default), lead, carbon
	TString targetmat = TString("carbon");
	TString alignmentFileName = par->getAlignmentFile();
	TString geoFile;
	TString realData;
	realData=(fopi) ? "kTRUE":"kFALSE";

	geoFile.Form("$FOPI2ROOT/fopigeometry/FopiGeom_s339%5.1f.root",targetpos);
	TFile* geotmp =new TFile(geoFile.Data(),"READ");
	if(geotmp->IsOpen())
	{
		std::cout<<"geofile exists already\n";
		geotmp->Close();
	}
	else
	{
		TString cmd;
		cmd.Form(".! root -b -q '$FOPI2ROOT/fopigeometry/FopiGeom.C(%f,\"%s\",\"%s\",\"%s\",%s)'",targetpos,targetmat.Data(),alignmentFileName.Data(),geoFile.Data(),realData.Data());
		gROOT->ProcessLine(cmd);
	}
	// add geometry file to fRun
	fRun->SetGeomFile(geoFile);

	if(bfield==1)
	{
		FOPIField *fMagField = new FOPIField(0.616);
		fMagField->SetTargetOffset(taroff);
		fRun->SetField(fMagField);
	}
	else
	{
		std::cout<<"running with 0 Magnetic Field"<<std::endl;
		PndConstField *fMagField=new PndConstField();
		fMagField->SetField(0., 0. , 0. ); // values are in kG
		fMagField->SetFieldRegion(-50, 50,-50, 50, -2000, 2000); // values are in cm
		fRun->SetField(fMagField);
	}


	//initialize GF field
	GFFieldManager::getInstance()->init(new PndFieldAdaptor(fRun->GetField()));
	GFMaterialEffects::getInstance()->init(new GFTGeoMaterialInterface());



	//--------------------SET UP TASKS ------------------------------
	TpcEventCounter* evCount =  new TpcEventCounter();
	if (nEvents!=nEvStart)
		evCount->SetnEvts(nEvents-nEvStart);
	fRun->AddTask(evCount);


	TpcDataReaderTask* read = new TpcDataReaderTask();
	read->SetPersistence(0);
	read->SetDatafile(filename);
	read->SetSampleBranchName("TpcSample");
	read->SetStartEvent(nEvStart);
	//read->SetCutSmallPad();
	//read->SetMinSamples(1000);
	fRun->AddTask(read);

	if (fopi)
	{


		FopiDataReaderTask* readFopi = new FopiDataReaderTask();
		readFopi->SetPersistence(persistence);
		readFopi->SetCdcHitPersistence(persistence);
		readFopi->SetCdcTrackPersistence(persistence);
		readFopi->SetRpcTrackPersistence(persistence);
		readFopi->SetBarrelOn();
		//not (re)implemented yet, but it is only a double-check
		//readFopi->SetRunNr(runNr);
		readFopi->SetInputFile(fopiFilename);
		//readFopi->SetNevent(nFopiEvents);
		readFopi->SetDisableHeliTracks(true);
		fRun->AddTask(readFopi);

		TString pythonScript2=basedir+"/macro/tpc/FOPI/getCdcPars.py";
		TPython::LoadMacro(pythonScript2);
		PyGetCdcPars cdcParLoader;
		TString bla = basedir+"/tpc/FOPI/par/cdc_par_AE.par";
		cdcParLoader.loadInputFile(bla.Data());

		//CORRECT FOR CDC SHIFTS
		double cdcT0=cdcParLoader.getT0(int(run));
		std::cout<<"cdcT0:"<<cdcT0<<std::endl;
		double cdcVDrift=cdcParLoader.getVDrift(int(run));
		std::cout<<"cdcVDrift:"<<cdcVDrift<<std::endl;
		double cdcLorentzAngle=cdcParLoader.getLorentzAngle(int(run));
		std::cout<<"cdcLorentzAngle:"<<cdcLorentzAngle<<std::endl;

		CdcCorrectionTask* corTask=new CdcCorrectionTask();
		corTask->SetCdcVdr(cdcVDrift);
		corTask->SetCdcLorentz(cdcLorentzAngle);
		corTask->SetCdcT0(cdcT0);
		//corTask->SetEnableFlightTimeCorr(true);
		corTask->SetCorrection(0.18); //default is 0.2
		//corTask->SetCorrectionFile("testShifts.root");
		corTask->SetCorrectionFile(zCorrFile);
		corTask->SetCdcT0Corr(0);
		corTask->SetVerbose(0);
		fRun->AddTask(corTask);


		//Conversion from fopi2root RpcTrack objects to GF strip hits
		FopiRpcInitTask* rpcInit = new FopiRpcInitTask();
		rpcInit->SetPersistence(persistence);
		rpcInit->SetUseConstantOffset(true);
		rpcInit->SetConstantErrorZ(2.5);
		rpcInit->SetConstantErrorPhi(0.2);
		rpcInit->SetOutBranchName("RpcPixHit");
		fRun->AddTask(rpcInit);

		//Conversion from fopi2root BarTrack objects to GF strip hits
		FopiBarInitTask* barInit = new FopiBarInitTask();
		barInit->SetPersistence(persistence);
		barInit->SetConstantErrorZ(8.);
		barInit->SetConstantErrorPhi(3./TMath::Sqrt(12.));
		barInit->SetOutBranchName("BarPixHit");
		//barInit->SetConstantErrorPhi(100.); //basically turning phi info off
		fRun->AddTask(barInit);
	}

	TpcPSATask* tpsa = new  TpcPSATask();
	tpsa->SetPersistence(0);
	tpsa->SetSamplePersistence(0);
	tpsa->SetSampleBranchName("TpcSample"); // Input of PSA
	fRun->AddTask(tpsa);


	if (gaincorr==1)
	{
		TpcDigiAmpCorrectionTask* tpcCalib = new TpcDigiAmpCorrectionTask();
		fRun->AddTask(tpcCalib);
	}

	TpcClusterFinderTask* tpcCF = new TpcClusterFinderTask();
	tpcCF->SetPersistence(persistence);
	if (clcorr)
		tpcCF->SetClusterBranchName("TpcCluster_raw");
	else if (reclust)
		tpcCF->SetClusterBranchName("TpcPreCluster");
	else
	{
		tpcCF->SetClusterBranchName("TpcCluster");
		tpcCF->SetPersistence(persistence); // keep Clusters
		tpcCF->SetDigiPersistence(persistence); // keep Digis copys in clusters
	}
	tpcCF->timeslice(6); //in samples
	tpcCF->SetErrorPars(-1,22);
	tpcCF->SetSingleDigiClusterAmpCut(0.);
	tpcCF->SetClusterAmpCut(0.);
	tpcCF->SetVerbose(iVerbose);
	tpcCF->SetSimpleClustering(); // use TpcClusterFinderSimple
	tpcCF->SetDigiPersistence(0); // keep reference to digis in clusters
	tpcCF->SetPersistence(persistence);
	tpcCF->SetUsePRF(false);
	//tpcCF->SetUsePRF(true); actually this makes no sense. prf correction also need track information
	fRun->AddTask(tpcCF);


	if (clcorr)
	{
		std::cout<<"Setting up cluster correction task"<<std::endl;
		TpcClusterCorrectionTask* tpcCC = new TpcClusterCorrectionTask();
		tpcCC->SetVerbose(0);
		tpcCC->SetPersistence(persistence);
		tpcCC->SetClusterBranchName("TpcCluster_raw");
		if (reclust)
			tpcCC->SetClusterOutBranchName("TpcPreCluster");
		else
			tpcCC->SetClusterOutBranchName("TpcCluster");
		//tpcCC->SetDevMap("datafiles/full_field_-148.93_-312.89_4s3_4.7__merged7.txt");
		tpcCC->SetDevMap(""); //if no devmap name for correction is given. it is read from the parfile
		fRun->AddTask(tpcCC);
	}

	std::cout<<"Setting up Rieman Task"<<std::endl;
	//find TpcRiemannTracks in the TPC alone
	TpcRiemannTrackingTask* tpcSPR1 = new TpcRiemannTrackingTask();
	if (reclust)
		tpcSPR1->SetClusterBranchName("TpcPreCluster");
	else
		tpcSPR1->SetClusterBranchName("TpcCluster");
	tpcSPR1->SetRiemannTrackBranchName("TpcRiemannTrack");
	tpcSPR1->SetRiemannHitBranchName("TpcRiemannHit");
	tpcSPR1->SetPersistence(persistence);
	tpcSPR1->SetSortingParameters(true,1,0);
	//tpcSPR1->SetTrkFinderParameters(2.0,1.6,3,0.4);
	tpcSPR1->SetTrkFinderParameters(2.6,1.7,3,1.45); //proxcut,helixcut,min hits,z-cut
	fRun->AddTask(tpcSPR1);


	std::cout<<"Setting up TrackInitTask"<<std::endl;
	//build GFTracks from TpcRiemannTracks
	TpcTrackInitTask* trackInit=new TpcTrackInitTask();
	trackInit->SetPersistence(persistence);
	trackInit->SetRiemannBranchName("TpcRiemannTrack");
	if (reclust)
	{
		trackInit->SetTrackOutBranchName("TpcPreTrackPreFit");
		trackInit->SetClusterBranchName("TpcPreCluster");
		trackInit->SetRecoHitOutBranchName("TpcPreSPHit");
	}
	else
	{
		trackInit->SetTrackOutBranchName("TpcTrackPreFit");
		trackInit->SetClusterBranchName("TpcCluster");
		trackInit->SetRecoHitOutBranchName("TpcSPHit");
	}
	if (fopi)
	{
		TVector3 IP(0,0,42.78);
		trackInit->SetInteractionPoint(IP);
		trackInit->SetPDG(211);
	}
	else
	{
		trackInit->SetCosmicSorting();
		trackInit->SetPDG(13);
	}
	trackInit->SetVerbose(iVerbose);
	//trackInit->SetWeightedPlaneConstruction(weightedPlaneConstruction);
	//trackInit->SetCutCov(cutCov);
	//trackInit->SetSimulation();
	trackInit->SetDipCut(0.05);
	trackInit->SetSmoothing(smoothing);
	trackInit->SetAlign(kFALSE);
	if(fG>0)
		trackInit->SetUseFullCov();
	fRun->AddTask(trackInit);

	if (reclust)
	{

		std::cout<<"Setting up KalmanTask 1"<<std::endl;
		KalmanTask* kalman1 =new KalmanTask();
		kalman1->SetTpcClusterBranchName("TpcPreSPHit");
		kalman1->SetTrackBranchName("TpcPreTrackPreFit");
		kalman1->SetOutBranchName("TpcTrackPrePostFit");
		kalman1->SetPersistence(persistence);
		kalman1->SetNumIterations(3); // number of fitting iterations (back and forth)
		fRun->AddTask(kalman1);

		std::cout<<"Setting up Reclusterizer"<<std::endl;
		TpcReClusterizerTask* tpcReCl = new TpcReClusterizerTask();
		tpcReCl->SetTrackBranchName("TpcTrackPrePostFit");
		tpcReCl->SetClusterBranchName("TpcPreCluster");
		tpcReCl->SetRecoHitBranchName("TpcPreSPHit");
		tpcReCl->SetReClusterBranchName("TpcCluster");
		tpcReCl->SetErrorPars(fG,fC);
		tpcReCl->SetAlign(kFALSE);
		//tpcReCl->SetStepSize(1);
		//tpcReCl->SetSkipUnFitted();
		tpcReCl->SetMinClusterSize(0);
		tpcReCl->SetUseFirstDigiPos();
		tpcReCl->SetVarStepSqrt();
		tpcReCl->SetStepOff_sqrt(stepoff);
		tpcReCl->SetFac(fac);
		tpcReCl->SetStepOff_0_sqrt(0);
		tpcReCl->SetVerbose(0);
		//tpcReCl->SetTimeInfo();
		tpcReCl->SetPersistence(persistence);
		tpcReCl->SetAllRecoHitPersistence(persistence);
		tpcReCl->SetDoPRF(doPRF);

		//tpcReCl->SetUseChamberEdge();
		if(!fopi)
			tpcReCl->SetUseCosmics();
		//tpcReCl->SetCutCov();
		//tpcReCl->SetWeightedPlaneConstruction(true);
		//tpcReCl->SetUseAllDigis();
		if(clcorr)
			tpcReCl->SetClusterCorr("");
		tpcReCl->SetTrackInit();

		if (fG>0)
		{
			tpcReCl->SetUseFullCov();

			if ( scaler>=55 )
			{
				tpcReCl->SetUseOnPlaneCov();
			}

			if(scaler==63)
			{
				//use the parametrization of the error directly inside the sphitoutfiles_e12/bat/bat_rev28544_target_fitCovOnPlane_cut_thirdIteration.pdf
				tpcReCl->SetUseSPHitParam(true);
			}
		}
		fRun->AddTask(tpcReCl);
	}
	//else
	//{
	TpcAlignmentTask* align = new TpcAlignmentTask();
	align->SetRecoHitBranchName("TpcSPHit");
	fRun->AddTask(align);
	//}

	KalmanTask* kalman =new KalmanTask();
	kalman->SetTpcClusterBranchName("TpcSPHit");
	kalman->SetTrackBranchName("TpcTrackPreFit");
	kalman->SetOutBranchName("TpcTrackPostFit");
	kalman->SetPersistence(persistence);
	kalman->SetVerbose(iVerbose);
	kalman->SetNumIterations(3); // number of fitting iterations (back and forth)
	kalman->SetUseDAF();
	fRun->AddTask(kalman);


	if(fopi)
	{
		//build GFTracks for CDC hits only
		FopiTrackInitTask* trackCdcInit=new FopiTrackInitTask();
		trackCdcInit->SetPersistence(persistence);  //affects prefit tracks
		trackCdcInit->SetAllRecoHitPersistence(persistence);
		trackCdcInit->SetExtrapToZ(true);
		//trackCdcInit->SetTrackPersistence(true);
		trackCdcInit->SetHitBranchName("CdcHit");
		trackCdcInit->SetTrackBranchName("CdcTrack");
		trackCdcInit->SetOutBranchName("CdcTrackPreTempFit");
		trackCdcInit->SetRecoHitOutBranchName("CdcSpacepoint");
		trackCdcInit->SetTrackPersistence(persistence);
		//trackCdcInit->SetPiKaonThreshold(0.32);
		//trackCdcInit->SetMinCDCHits(5);
		trackCdcInit->SetSmoothing(smoothing);
		//trackCdcInit->SetUsePseudoSPHits();
		trackCdcInit->SetWeightedPlaneConstruction(true);
		trackCdcInit->SetVerbose(iVerbose);
		//trackCdcInit->SetCutCov(true);
		//in cm:
		trackCdcInit->SetConstantPointErrors(3.5e-2, 3.5e-2, 10.);
		fRun->AddTask(trackCdcInit);


		CollectParticlesTask* tpcOnlyColl=new CollectParticlesTask();
		tpcOnlyColl->SetVerbose(0);
		tpcOnlyColl->SetPersistence(kTRUE);
		tpcOnlyColl->SetTpcHitBranchName("TpcSPHit");
		tpcOnlyColl->SetCdcHitBranchName("CdcSpacepoint");
		tpcOnlyColl->SetTrackBranchName("TpcTrackPostFit");
		tpcOnlyColl->SetOutBranchName("TpcOnlyParticles");
		//fRun->AddTask(tpcOnlyColl);

		//track fitting for cdc
		KalmanTask* KalmanCdc = new KalmanTask();
		KalmanCdc->SetPersistence(persistence);
		KalmanCdc->SetTrackBranchName("CdcTrackPreTempFit");
		KalmanCdc->SetFopiCDCBranchName("CdcSpacepoint");
		KalmanCdc->SetOutBranchName("CdcTrackPostTempFit");
		KalmanCdc->SetNumIterations(3); // number of fitting iterations
		KalmanCdc->SetUseDAF();
		KalmanCdc->SetVerbose(iVerbose);
		fRun->AddTask(KalmanCdc);

		//match TPC and CDC tracks (no GFTracks!)
		TpcCdcMatchingTask* tpcCdcMatch = new TpcCdcMatchingTask();
		tpcCdcMatch->SetPersistence(persistence);
		tpcCdcMatch->SetTpcTrackBranchName("TpcTrackPostFit");
		tpcCdcMatch->SetCdcTrackBranchName("CdcTrack");
		tpcCdcMatch->SetGFCdcTrackBranchName("CdcTrackPostTempFit");
		tpcCdcMatch->SetCdcHitBranchName("CdcHit");
		tpcCdcMatch->SetWriteFailedMatches();
		tpcCdcMatch->SetVerbose(iVerbose);
		fRun->AddTask(tpcCdcMatch);

		FopiTrackMergingTask* merge = new FopiTrackMergingTask();
		merge->SetPersistence(persistence);
		merge->SetTrackPersistence(persistence); //prefit tracks
		merge->SetTpcTrackBranchName("TpcTrackPostFit");
		merge->AddTupleBranchName("TpcCdcMatchingPair");
		merge->SetCdcGFTrackBranchName("CdcTrackPreTempFit");
		merge->SetTrackOutBranchName("FopiTrackPreTempFit");
		merge->SetVerbose(iVerbose);
		fRun->AddTask(merge);


		KalmanTask* KalmanComb = new KalmanTask();
		KalmanComb->SetUseDAF();
		KalmanComb->SetPersistence(persistence);
		KalmanComb->SetTpcClusterBranchName("TpcSPHit");
		KalmanComb->SetFopiCDCBranchName("CdcSpacepoint");
		KalmanComb->SetTrackBranchName("FopiTrackPreTempFit");
		KalmanComb->SetOutBranchName("FopiTrackPostTempFit");
		KalmanComb->SetNumIterations(3); // number of fitting iterations
		KalmanComb->SetVerbose(iVerbose);
		//fRun->AddTask(KalmanComb);

		CollectParticlesTask* tpccdc1Coll=new CollectParticlesTask();
		tpccdc1Coll->SetVerbose(0);
		tpccdc1Coll->SetPersistence(kTRUE);
		tpccdc1Coll->SetTpcHitBranchName("TpcSPHit");
		tpccdc1Coll->SetCdcHitBranchName("CdcSpacepoint");
		tpccdc1Coll->SetTrackBranchName("FopiTrackPostTempFit");
		tpccdc1Coll->SetOutBranchName("TpcCdc1Particles");
		//fRun->AddTask(tpccdc1Coll);

		//RPC / Barrel merging goes here. RPC residual check is performed
		//creates both TPC-CDC + BAR/RPC and CDC + BAR/RPC
		FopiRpcBarMergingTask* addBarRpc = new FopiRpcBarMergingTask();
		addBarRpc->SetPersistence(persistence);
		addBarRpc->SetExtrapToZ(true);
		addBarRpc->SetSmoothing(smoothing);
		addBarRpc->SetTrackBranchName("FopiTrackPreTempFit");
		addBarRpc->SetCdcGFTrackBranchName("CdcTrackPostTempFit");
		addBarRpc->SetTrackOutBranchName("FopiTrackPreFit");
		addBarRpc->SetCdcTrackOutBranchName("CdcTrackPreFit");
		addBarRpc->SetDissolveFOPIRpcMatching(false); //DISSOLVE FOPI's CDC-RPC matching if residual too large
		//addBarRpc->SetExtrapResMeanZ(2.82e-1);
		//addBarRpc->SetExtrapResSigZ(1.40);
		//addBarRpc->SetExtrapWindow(2.);
		addBarRpc->SetVerbose(iVerbose);
		fRun->AddTask(addBarRpc);

		//refit Cdc with added BAR / RPC hits
		KalmanTask* KalmanCombCdc = new KalmanTask();
		KalmanCombCdc->SetUseDAF();
		KalmanCombCdc->SetPersistence(persistence);
		//KalmanCombCdc->SetTpcClusterBranchName("TpcSPHit");
		KalmanCombCdc->SetFopiCDCBranchName("CdcSpacepoint");
		KalmanCombCdc->SetTrackBranchName("CdcTrackPreFit");
		KalmanCombCdc->SetOutBranchName("CdcTrackPostFit");
		KalmanCombCdc->SetNumIterations(3); // number of fitting iterations
		KalmanCombCdc->SetVerbose(iVerbose);
		fRun->AddTask(KalmanCombCdc);

		//refit all with added BAR / RPC hits
		KalmanTask* KalmanComb2 = new KalmanTask();
		KalmanComb2->SetUseDAF();
		KalmanComb2->SetPersistence(persistence);
		KalmanComb2->SetTpcClusterBranchName("TpcSPHit");
		KalmanComb2->SetFopiCDCBranchName("CdcSpacepoint");
		KalmanComb2->SetTrackBranchName("FopiTrackPreFit");
		KalmanComb2->SetOutBranchName("FopiTrackPostFit");
		KalmanComb2->SetNumIterations(3); // number of fitting iterations
		KalmanComb2->SetVerbose(iVerbose);
		fRun->AddTask(KalmanComb2);

		// AddTrackRepTask* addReps = new AddTrackRepTask();
		// addReps->SetTrackBranchName("FopiTrackPreFit");
		// addReps->AddHypothesis("proton");
		// addReps->AddHypothesis("kaon");
		// //  fRun->AddTask(addReps);


		// //refit with added BAR / RPC hits
		// KalmanTask* KalmanComb2 = new KalmanTask();
		// KalmanComb2->SetUseDAF();
		// KalmanComb2->SetPersistence(kFALSE);
		// KalmanComb2->SetTpcClusterBranchName("TpcSPHit");
		// KalmanComb2->SetFopiCDCBranchName("CdcSpacepoint");
		// KalmanComb2->SetTrackBranchName("FopiTrackPreFit");
		// KalmanComb2->SetOutBranchName("FopiTrackPostFit");
		// //KalmanComb2->SetNumIterations(2); // number of fitting iterations
		// KalmanComb2->AddIgnorePDG(-2212);  //ignore antiprotons in the fit
		// //  fRun->AddTask(KalmanComb2);


		// dE/dx task
		TpcdEdxTask* dEdx = new TpcdEdxTask();
		dEdx->SetTrackBranchName("FopiTrackPostFit");
		dEdx->SetOutBranchName("FopidEdx");
		dEdx->SetClusterBranchName("TpcSPHit");
		//dEdx->SetIdealdEdx();E
		dEdx->SetPersistence(persistence);
		dEdx->SetVerbose(0);
		dEdx->SetUseFirstDigiPos();
		dEdx->SetIgnoreEdgeDigis(5.8,14.5);
		dEdx->SetDXgrid(0.5);
		dEdx->SetVerbose(iVerbose);
		fRun->AddTask(dEdx);

		CollectParticlesTask* coll=new CollectParticlesTask();
		coll->SetVerbose(0);
		coll->SetPersistence(kTRUE);
		coll->SetTrackBranchName("FopiTrackPostFit");
		coll->SetdEdxBranchName("FopidEdx");
		coll->SetSkipUnfitted();
		coll->SetdEdxOnly();
		fRun->AddTask(coll);

		CollectParticlesTask* tpcColl=new CollectParticlesTask();
		tpcColl->SetVerbose(0);
		tpcColl->SetPersistence(kTRUE);
		tpcColl->SetTpcHitBranchName("TpcSPHit");
		tpcColl->SetTrackBranchName("TpcTrackPostFit");
		tpcColl->SetOutBranchName("TpcParticles");
		//fRun->AddTask(tpcColl);

		CollectParticlesTask* cdcColl=new CollectParticlesTask();
		cdcColl->SetVerbose(0);
		cdcColl->SetPersistence(kTRUE);
		cdcColl->SetTpcHitBranchName("TpcSPHit");
		cdcColl->SetCdcHitBranchName("CdcSpacepoint");
		cdcColl->SetTrackBranchName("CdcTrackPostTempFit");
		cdcColl->SetOutBranchName("CdcParticles");
		//fRun->AddTask(cdcColl);

		if(1)
		{
			//the actual calculator:
			TpcCdc2DMatchedResCalc* cdcRefRes = new TpcCdc2DMatchedResCalc();
			cdcRefRes->addBranchName("FopiTrackTuples");
			cdcRefRes->addBranchName("TpcSPHit");
			cdcRefRes->setTpcTrackBranchName("TpcTrackPostFit");
			cdcRefRes->setCdcTrackBranchName("CdcTrack");
			//the task:
			TpcRefTrackResidualTask* res = new TpcRefTrackResidualTask();
			res->SetResCalculator(cdcRefRes);
			res->SetOutBranch("TpcCdcRefResiduals");
			res->SetPersistence();
			//fRun->AddTask(res);

			TpcQATask* qa = new TpcQATask();
			qa->SetTpcCdcTrackBranchName("FopiTrackPostFit");
			qa->SetTpcDecodedFile(filename);
			//fRun->AddTask(qa);
		}
	}

	if( (fG > 0) && (scaler==63) )
	{
		std::cout<<"Setting up ClusterErrorRecalcTask"<<std::endl;
		TpcClusterErrorRecalcTask* clErrRec=new TpcClusterErrorRecalcTask();
		fRun->AddTask(clErrRec);
	}

	if (! fopi)
	{
		std::cout<<"Setting up CosmicsTask"<<std::endl;
		TpcCosmicsTask* Cos = new TpcCosmicsTask();
		Cos->SetClusterBranchName("TpcCluster");
		Cos->SetPersistence(0);
		Cos->SetNumberOfTrackReps(1); // set to 2 if you use GeaneTrackrep (tpcSPR->useGeane();)
		Cos->SetDatafile(filename);
		if (smoothing==1)
			Cos->SetUnbiased();
		//Cos->SetVerbose();
		fRun->AddTask(Cos);

		std::cout<<"Setting up Track Cleaner"<<std::endl;
		TpcTrackCleaner* tpcTrCl = new TpcTrackCleaner();
		tpcTrCl->SetCosmics();
		//tpcTrCl->SetThetaCut(35,145);
		tpcTrCl->SetClMeanCut(200);
		tpcTrCl->SetTrackLenCut(3);
		//tpcTrCl->SetPhiCut(70,110);
		tpcTrCl->SetVerbose(iVerbose);
		tpcTrCl->SetKeepDelTracks();
		fRun->AddTask(tpcTrCl);
	}

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


	std::cout<<"************start ini***************"<<std::endl;
	fRun->Init();
	//Initialize te Digimapper:
	TpcDigiPar* tpcpar = FairRun::Instance()->GetRuntimeDb()->getContainer("TpcDigiPar");
	TpcDigiMapper::getInstance()->init(tpcpar);
	//TpcDigiMapper::getInstance()->forceManualDriftVel(kTRUE);
	TpcAlignmentManager::init(tpcpar->getAlignmentFile());
	//TpcAlignmentManager::getInstance();
	std::cout<<"********************************************"<<std::endl;
	std::cout<<"pads:"<<TpcDigiMapper::getInstance()->getPadPlane()->GetNPads()<<std::endl;
	//rtdb->print();
	std::cout<<"start:"<<nEvStart<<" end:"<<nEvents<<std::endl;
	if (outnum<=0)
		rtdb->saveOutput();

	switch (scaler)
	{
	case 0:
		TpcClusterErrorScaler::getInstance()->init("./tpc/FOPI/par/tpc.errorScale.recl.physNO1.86.par",0);
		break;
	case 10:
		TpcClusterErrorScaler::getInstance()->init("./tpc/FOPI/par/tpc.errorScale.recl.phys.86.par",0);
		break;
	case 1:
		TpcClusterErrorScaler::getInstance()->init("./tpc/FOPI/par/tpc.errorScale.recl.physNum.86.par",0);
		break;
	case 15:
		TpcClusterErrorScaler::getInstance()->init("./tpc/FOPI/par/tpc.errorScale.recl.physNumProj.86.par",0);
		break;
	case 4:
		TpcClusterErrorScaler::getInstance()->init("./tpc/FOPI/par/tpc.errorScale.recl.physMeanUV.86.par",4);
		break;
	case 44:
		TpcClusterErrorScaler::getInstance()->init("./tpc/FOPI/par/tpc.errorScale.recl.physMeanUVcut.86.par",4);
		break;
	case 20:
		TpcClusterErrorScaler::getInstance()->init("./tpc/FOPI/par/tpc.errorScale.recl.exp.86.par",20);
		break;
	case 98:
		TpcClusterErrorScaler::getInstance()->init("./tpc/FOPI/par/tpc.errorScale.recl.Ax1_1000.par",0);
		break;
	case 55:
		break;
	case 61: //parametrization in XYZ and cut by plane
		TpcClusterErrorScaler::getInstance()->init("./tpc/FOPI/par/tpc.errorParam.86.data.par",61);
		break;
	case 63:
		//all these parameter sets were obtaind from cosmics data with field:86 gas:ar gain:85
		//TpcClusterErrorScaler::getInstance()->init("./outfiles_e12/bat/bat_rev28521_data_fitCovOnPlane_dataCuts_firstIteration_param.txt",61);

		//TpcClusterErrorScaler::getInstance()->init("./outfiles_e12/bat/bat_rev28546_data_fitCovOnPlane_dataCuts_secondIteration_param.txt",61);
		//TpcClusterErrorScaler::getInstance()->init("./outfiles_e12/bat/bat_rev28528_data_fitCovOnPlane_dataCuts_secondIteration_param.txt",61);

		//TpcClusterErrorScaler::getInstance()->init("./outfiles_e12/bat/bat_rev28578_data_fitCovOnPlane_dataCuts_firstIteration_param.txt",61);
		//TpcClusterErrorScaler::getInstance()->init("./outfiles_e12/bat/bat_rev28585_data_fitCovOnPlane_dataCuts_secondIteration_param.txt",61);
		//TpcClusterErrorScaler::getInstance()->init("./outfiles_e12/bat/bat_rev28591_data_fitCovOnPlane_dataCuts_thirdIteration_param.txt",61);
		//TpcClusterErrorScaler::getInstance()->init("./outfiles_e12/bat/bat_rev28604_data_fitCovOnPlane_dataCuts_fourthIteration_param.txt",61);

		//new iterations with z-field-correction of cluster
		//TpcClusterErrorScaler::getInstance()->init("./outfiles_e12/bat/bat_rev28640_data_fitCovOnPlane_dataCuts_firstIteration_param.txt",61);
		TpcClusterErrorScaler::getInstance()->init("./outfiles_e12/bat/bat_rev28645_data_fitCovOnPlane_dataCuts_secondIteration_param.txt",61);

		break;
	case 64:
		TpcClusterErrorScaler::getInstance()->init("",61);
		break;
	case 56:
		TpcClusterErrorScaler::getInstance()->init("./outfiles_e12/error_lookup.txt",56);
		break;
	case 99:
		TpcClusterErrorScaler::getInstance()->init("",-1);
		break;
	}
	RooMsgService::instance().setGlobalKillBelow(RooFit::ERROR);
	RooMsgService::instance().setSilentMode(true);

	fRun->Run(nEvStart,nEvents);

	if (outnum<=0)
		rtdb->saveOutput();

	rtdb->print();
	timer.Stop();
	Double_t rtime = timer.RealTime();
	Double_t ctime = timer.CpuTime();
	printf("RealTime=%f seconds, CpuTime=%f seconds\n",rtime,ctime);

	std::cout<<"OutputFile: "<<outFile<<std::endl;
}