{ // ======================================================================== // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug) Int_t iVerbose = 0; // Input file TString inDigiFile = "evt_digi_tpc.root"; TString inSimFile = "evt_points_tpc.root"; // Parameter file TString parFile = "evt_params_tpc.root"; // Output file TString outFile = "evt_reco_tpc.root"; // Number of events to process Int_t nEvents = 0; // ---- Load libraries ------------------------------------------------- gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C"); rootlogon(); TString sysFile = gSystem->Getenv("VMCWORKDIR"); // ------------------------------------------------------------------------ // In general, the following parts need not be touched // ======================================================================== // ----- Timer -------------------------------------------------------- TStopwatch timer; timer.Start(); // ------------------------------------------------------------------------ // ----- Digitization run ------------------------------------------- FairRunAna *fRun= new FairRunAna(); fRun->SetInputFile(inDigiFile); fRun->AddFriend(inSimFile); fRun->SetOutputFile(outFile); FairGeane *Geane = new FairGeane(); fRun->AddTask(Geane); // ------------------------------------------------------------------------ // ----- Parameter database -------------------------------------------- TString allDigiFile = sysFile+"/macro/params/all.par"; FairRuntimeDb* rtdb = fRun->GetRuntimeDb(); FairParRootFileIo* parInput1 = new FairParRootFileIo(); parInput1->open(parFile.Data()); FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo(); parIo1->open(allDigiFile.Data(),"in"); rtdb->setFirstInput(parInput1); rtdb->setSecondInput(parIo1); PndGeoHandling* geoH = PndGeoHandling::Instance(); // ------------------------------------------------------------------------ // ------- RECO procedure ------------------------------------------------ //correct for unfortunate shift in TPC digi //PndTpcRoughAlignmentTask* align = new PndTpcRoughAlignmentTask(); //align->SetShift(TVector3(0.,0.,-3.71357e-01)); //old PSA //align->SetShift(TVector3(0.,0.,5.6E-2)); //new PSA //fRun->AddTask(align); //find track candidates in the TPC alone PndTpcRiemannTrackingTask* tpcSPR = new PndTpcRiemannTrackingTask(); tpcSPR->SetPersistence(); tpcSPR->SetSortingParameters( true, // false: sort only according to _sorting (see next argument); true: use internal sorting when adding hits to trackcands 3, // -1: no sorting, 0: sort Clusters by X, 1: Y, 2: Z, 3: R, 4: distance to origin 0.); // z-position of interaction point (for sorting 4) tpcSPR->SetTrkFinderParameters( 1.9, // proximity cut in 3D [cm] 0.4, // helix cut [cm] 5); // minimum hits for helix-fit tpcSPR->SetMergeTracks(); tpcSPR->SetTrkMergerParameters( 2.5, // proximity cut [cm] 0.1, // dip cut [rad] 0.6, // helix cut [cm] 0.025);// plane cut (RMS) //tpcSPR->SetRiemannScale(); // sets riemannscale for the prototype; tpcSPR->useGeane(); // use RKTrackrep and GeaneTrackrep tpcSPR->SetSmoothing(true); tpcSPR->SetMCPid(kFALSE); // use ideal particle identification fRun->AddTask(tpcSPR); KalmanTask* kalman =new KalmanTask(); kalman->SetPersistence(); kalman->SetNumIterations(3); // number of fitting iterations (back and forth) fRun->AddTask(kalman); // creates TrackPostFit branch //correlate fitted track with MVD pixels and strips PndTpcMVDCorrelatorTask* corr = new PndTpcMVDCorrelatorTask(); corr->SetMatchDistance(200.); //mutliple of MVD hit sigma (which 100 -> roughly 20 mu) corr->SetMinMVDHits(1); corr->SetOutTrackBranchName("TrackPreFitMVD"); corr->SetPersistence(true); fRun->AddTask(corr); //fit after MVD corr KalmanTask* kalman2 =new KalmanTask(); kalman2->SetPersistence(); kalman2->SetNumIterations(3); // number of fitting iterations (back and forth) kalman2->SetTrackBranchName("TrackPreFitMVD"); kalman2->SetOutBranchName("TrackPostFitMVD"); fRun->AddTask(kalman2); PndTpcGEMCorrelatorTask* corrG = new PndTpcGEMCorrelatorTask(); corrG->SetMatchDistance(100.); //mutliple of GEM hit sigma corrG->SetMinGEMHits(2); corrG->SetTrackBranchName("TrackPostFitMVD"); corrG->SetOutTrackBranchName("TrackPreFitGEM"); corrG->SetPersistence(true); fRun->AddTask(corrG); //final fit KalmanTask* kalman3 =new KalmanTask(); kalman3->SetPersistence(); kalman3->SetNumIterations(3); // number of fitting iterations (back and forth) kalman3->SetTrackBranchName("TrackPreFitGEM"); kalman3->SetOutBranchName("TrackPostFitComplete"); fRun->AddTask(kalman3); PndGFTrackToPndTrackConvertorTask* converter =new PndGFTrackToPndTrackConvertorTask(); converter->SetTrackInBranchName("TrackPostFitComplete"); converter->SetTrackOutBranchName("PndTrackPostFitComplete"); fRun->AddTask(converter); PndMCTrackAssociator* trackMC = new PndMCTrackAssociator(); trackMC->SetTrackInBranchName("PndTrackPostFitComplete"); trackMC->SetTrackOutBranchName("TrackPostFitCompleteID"); fRun->AddTask(trackMC); // ----- Intialise and run -------------------------------------------- PndEmcMapper::Init(6); fRun->Init(); fRun->Run(0, nEvents); rtdb->saveOutput(); rtdb->print(); corr->WriteHistograms("MVDRes.root"); corrG->WriteHistograms("GEMRes.root"); // ------------------------------------------------------------------------ // ----- 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; // ------------------------------------------------------------------------ }