void runReco(TString digifile, int field, TString gas, int gain,bool reclust, unsigned int nEvents = 0, unsigned int nEvStart = 0, unsigned int outnum=0, int fADCSens=600, int fC=400) { gRandom->SetSeed(1234); // ======================================================================== // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug) Int_t iVerbose = 0; //gROOT->ProcessLine(".x $VMCWORKDIR/gconfig/rootlogon.C"); // Input file // Input file TString inDigiFile = digifile; TString inSimFile = inDigiFile; inSimFile.ReplaceAll("raw.root", "mc.root"); // Parameter file TString parFile = inDigiFile; parFile.ReplaceAll("raw.root", "param.root"); // Output file TString outFile = inDigiFile; outFile.ReplaceAll(".raw.root", Form("_%i_%i.reco.root", fADCSens, fC)); //if (outnum>0) // outFile.ReplaceAll(".root",Form("_%i.root", outnum)); // Number of events to process // Int_t nEvents = 0; // ---- Load libraries ------------------------------------------------- TString sysFile = gSystem->Getenv("VMCWORKDIR"); // ------------------------------------------------------------------------ // In general, the following parts need not be touched // ======================================================================== // ----- Timer -------------------------------------------------------- TStopwatch timer; timer.Start(); // ------------------------------------------------------------------------ TDatabasePDG::Instance()->AddParticle("deuteron","deuteron", 1.875612793, kTRUE, 0, 3,"nucleus", 1000010020, -1, 45); // ----- Digitization run ------------------------------------------- FairRunAna *fRun= new FairRunAna(); fRun->SetInputFile(inDigiFile); cout<<"Set Input File: "<AddFriend(inSimFile); cout<<"Set Friend MC File: "<SetOutputFile(outFile); cout<<"Set Output File: "<Getenv("VMCWORKDIR"); tpcDigiFile += Form("/tpc/FOPI/par/tpc.%i%s%iMC_long.par",field,gas.Data(),gain); // parInput1->open(tpcDigiFile.Data(),"in"); FairRuntimeDb* rtdb = fRun->GetRuntimeDb(); FairParRootFileIo* parInput1 = new FairParRootFileIo(); parInput1->open(parFile.Data()); FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo(); parIo1->open(tpcDigiFile.Data(),"in"); rtdb->setFirstInput(parInput1); rtdb->setSecondInput(parIo1); //PndGeoHandling* geoH = PndGeoHandling::Instance(); TString geoFile="$VMCWORKDIR/geometry/tpc_prototype_long_ArCO2.root"; fRun->SetGeomFile(geoFile); // ------- RECO procedure ------------------------------------------------ //TpcAlignmentManager::init(basedir+"/tpc/parfiles/dummyAlignment.Prototype.txt"); //TpcAlignmentManager * align = TpcAlignmentManager::getInstance(); TpcEventCounter* evCount = new TpcEventCounter(); evCount->SetnEvts(nEvents); //evCount->SetStep(1); fRun->AddTask(evCount); TpcClusterFinderTask* tpcCF = new TpcClusterFinderTask(); if (reclust) tpcCF->SetClusterBranchName("TpcPreCluster"); else tpcCF->SetClusterBranchName("TpcCluster"); tpcCF->SetDigiPersistence(); // keep reference to digis in clusters tpcCF->SetPersistence(); // keep Clusters tpcCF->timeslice(8); //in samples tpcCF->SetThreshold(0); tpcCF->SetSingleDigiClusterAmpCut(0.); tpcCF->SetClusterAmpCut(0.); // cut on mean digi amplitude tpcCF->SetErrorPars(fADCSens, fC); tpcCF->SetSimpleClustering(); // use TpcClusterFinderSimple tpcCF->SetVerbose(0); fRun->AddTask(tpcCF); //the actual calculator: TpcMCResCalcSample* mc3dSample = new TpcMCResCalcSample(); mc3dSample->addBranchName("TpcPoint"); mc3dSample->addBranchName("TpcSample"); mc3dSample->skipSecondaries(); //the task: TpcRefTrackResidualTask* resSample = new TpcRefTrackResidualTask(); resSample->SetOutBranch("MC_Sample_Residuals"); resSample->SetResCalculator(mc3dSample); resSample->SetPersistence(); //fRun->AddTask(resSample); //the actual calculator: TpcMCResCalcDigi* mc3dDigi = new TpcMCResCalcDigi(); mc3dDigi->addBranchName("TpcPoint"); mc3dDigi->addBranchName("TpcDigi"); mc3dDigi->skipSecondaries(); //the task: TpcRefTrackResidualTask* resDigi = new TpcRefTrackResidualTask(); resDigi->SetOutBranch("MC_Digi_Residuals"); resDigi->SetResCalculator(mc3dDigi); resDigi->SetPersistence(); //fRun->AddTask(resDigi); //the actual calculator: if (reclust) { cout<<"Setting up Rieman Task"<SetClusterBranchName("TpcPreCluster"); tpcSPR->SetRiemannTrackBranchName("PreRiemannTrack"); tpcSPR->SetRiemannHitBranchName("PreRiemannHit"); //tpcSPR->SetPersistence(); tpcSPR->SetMergeTracks(); //tpcSPR->SetTrkFinderParameters(10,4,3,2);//prox,helix,minpoints,zstretch tpcSPR->SetTrkMergerParameters(40,10,10,10); //prox,dip,helix,plane //tpcSPR->SetRiemannScale(10); //tpcSPR->SetSortingParameters(true,1,0); //tpcSPR->SetVerbose(1); fRun->AddTask(tpcSPR); cout<<"Setting up TrackInitTask"<SetPersistence(); trackInit->SetRiemannBranchName("PreRiemannTrack"); trackInit->SetTrackOutBranchName("TpcPreTrackPreFit"); trackInit->SetClusterBranchName("TpcPreCluster"); trackInit->SetRecoHitOutBranchName("TpcPreSPHit"); trackInit->SetVerbose(1); trackInit->SetPDG(211); trackInit->SetSmoothing(true); fRun->AddTask(trackInit); // cout<<"Setting up AlignmentTask"<SetRecoHitBranchName("TpcPreSPHit"); // fRun->AddTask(align); cout<<"Setting up KalmanTask 1"<SetTpcClusterBranchName("TpcPreSPHit"); kalman1->SetTrackBranchName("TpcPreTrackPreFit"); kalman1->SetOutBranchName("TrackPrePostFit"); //kalman1->SetPersistence(); kalman1->SetNumIterations(3); // number of fitting iterations (back and forth) fRun->AddTask(kalman1); cout<<"Setting up Reclusterizer"<SetTrackBranchName("TrackPrePostFit"); tpcReCl->SetClusterBranchName("TpcPreCluster"); tpcReCl->SetRecoHitBranchName("TpcPreSPHit"); tpcReCl->SetReClusterBranchName("TpcCluster"); tpcReCl->SetStepsize(1); //tpcReCl->SetVerbose(); tpcReCl->SetPersistence(); tpcReCl->SetUseChamberEdge(); tpcReCl->SetUseAllDigis(); tpcReCl->SetRDigiCutoff(2); //tpcReCl->SetUseFirstDigiPos(); fRun->AddTask(tpcReCl); // //the actual calculator: // TpcMCResCalcCluster* mc3dCluster2 = new TpcMCResCalcCluster(); // mc3dCluster2->addBranchName("TpcPoint"); // mc3dCluster2->addBranchName("TpcCluster"); // mc3dCluster2->skipSecondaries(); // //the task: // TpcRefTrackResidualTask* resCluster2 = new TpcRefTrackResidualTask(); // resCluster2->SetOutBranch("MC_Cluster_Residuals"); // resCluster2->SetResCalculator(mc3dCluster); // resCluster2->SetPersistence(); // fRun->AddTask(resCluster2); } TpcMCResCalcCluster* mc3dCluster = new TpcMCResCalcCluster(); mc3dCluster->addBranchName("TpcPoint"); mc3dCluster->addBranchName("TpcCluster"); mc3dCluster->skipSecondaries(); //the task: TpcRefTrackResidualTask* resCluster = new TpcRefTrackResidualTask(); resCluster->SetOutBranch("MC_Cluster_Residuals"); resCluster->SetResCalculator(mc3dCluster); resCluster->SetPersistence(); fRun->AddTask(resCluster); cout<<"Setting up Rieman Task"<SetClusterBranchName("TpcCluster"); tpcSPR2->SetPersistence(); tpcSPR2->SetSortingParameters(true,1,0); //tpcSPR->SetVerbose(1); fRun->AddTask(tpcSPR2); TpcTrackInitTask* trackInit2=new TpcTrackInitTask(); trackInit2->SetPersistence(); trackInit2->SetTrackOutBranchName("TpcTrackPreFit"); trackInit2->SetClusterBranchName("TpcCluster"); trackInit2->SetRecoHitOutBranchName("TpcSPHit"); trackInit2->SetVerbose(1); trackInit2->SetPDG(211); trackInit2->SetSmoothing(true); fRun->AddTask(trackInit2); cout<<"Setting up KalmanTask 2"<SetTpcClusterBranchName("TpcSPHit"); kalman2->SetTrackBranchName("TpcTrackPreFit"); kalman2->SetOutBranchName("TrackPostFit"); kalman2->SetPersistence(); kalman2->SetNumIterations(3); // number of fitting iterations (back and forth) fRun->AddTask(kalman2); cout<<"Setting up ResidualTask"<SetPersistence(); Res->SetNumberOfTrackReps(1); // set to 2 if you use GeaneTrackrep (tpcSPR->useGeane();) Res->SetUnbiased(); fRun->AddTask(Res); // ----- Intialise and run -------------------------------------------- fRun->Init(); // std::cout<<"Number of events to process:"<GetRuntimeDb()->getContainer("TpcDigiPar"); TpcDigiMapper::getInstance()->init(tpcpar); TpcAlignmentManager::init(tpcpar->getAlignmentFile()); std::cout<<"post init"<Run(nEvStart, nEvents); rtdb->saveOutput(); rtdb->print(); // ------------------------------------------------------------------------ // ----- Finish ------------------------------------------------------- timer.Stop(); Double_t rtime = timer.RealTime(); Double_t ctime = timer.CpuTime(); cout << endl << endl; cout << "Macro finished succesfully." << 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; // ------------------------------------------------------------------------ }