// ----------------------------------------------------------------------------- // ----------Reconstruction Macro for ALICE GEM TPC IROC Test in PS------------- // ----------------------------------------------------------------------------- // ----------Philipp Gadow E18/E12 Bachelor student July 2013------------------- // ----------------------------------------------------------------------------- //Description: // This macro starts the reconstruction process of particle tracks from the samples obtained in the PS testrun in Winter 2012 of the GEM IROC prototype. // The input data is in the decodedInPath, the output data will be saved in the outpath // The Reconstruction macro calls different tasks: // Data Reader Task: Reads the samples from input file // Pulse Shape Analysis Task: Processes samples into digis // Gain Correction Task: Uses the gain map of the analysis of Jens Wiechula, is set in the parameter file (par folder) // Clustering Task: Use ClusterFinderSimple (and look in Preliminary Cluster class!), perhaps experiment with carefull cuts on small amplitudes? // Pattern Recognition Tasks: Use either Riemann PR or Hough PR, also set proper PR in TrackInit-Task! // Track Fitting: Fitting of Tracks, not much to change // dEdx-Task: The collection region around the track can be varied void runRecoALICE(int runNumber=681, unsigned int nEvents = 0, unsigned int nEvStart = 0, unsigned int dataNumber = 99, unsigned int smoothing = 1, unsigned int gaincorr = 1, unsigned int usehough = 1, TString decodedInPath="/nfs/mds/data/tpc/alice/ps_test_beam/convertedDigis/", TString outpath="/nfs/mds/data/tpc/alice/ps_test_beam/RecoOut", bool forceDriftVel = false ) { // ======================================================================== // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug) Int_t iVerbose = 0; //Stopwatch for performance checks TStopwatch timer; timer.Start(); //-------------------- SET UP ROOT --------------------------------------- gROOT->SetBatch(); //to ensure running in batch TString basedir = gSystem->Getenv("VMCWORKDIR"); gROOT->ProcessLine(".x $VMCWORKDIR/gconfig/rootlogon.C"); FairRunAna* fRun = new FairRunAna(); TString jobname = ""; jobname += "merge_"; jobname += runNumber; jobname += "_"; jobname += dataNumber; if(gaincorr) jobname += "_corrected"; if(usehough) jobname += "_hough"; if(smoothing) jobname += "_smoothed"; jobname += ".reco.root"; TString outFile = outpath+"/"; outFile += jobname; fRun->SetOutputFile(outFile); std::stringstream ss; ss<GetRuntimeDb(); FairParAsciiFileIo* parInput1 = new FairParAsciiFileIo(); TString parFile = "/tpc/ALICE/par/tpc.iroc.par"; parFile = basedir+parFile; std::cout << "parFile: " << parFile << std::endl; parInput1->open(parFile.Data(),"in"); rtdb->setFirstInput(parInput1); TpcDigiPar* par = (TpcDigiPar*) rtdb->getContainer("TpcDigiPar"); par->setInputVersion(fRun->GetRunId(),1); par->setChanged(kTRUE); par->init(); //--- add geometry file to fRun --- TString geoFile = basedir + "/tpc/ALICE/par/ALICEGeom_test.root"; fRun->SetGeomFile(geoFile); PndConstField *fMagField=new PndConstField(); fMagField->SetField(0., 0. , 0. ); // values are in kG fMagField->SetFieldRegion(-500, 500,-500, 500, -2000, 2000); // values are in cm fRun->SetField(fMagField); //--- load input file --- TString inputfile = ""; inputfile += decodedInPath; inputfile += "merge_"; inputfile += runNumber; inputfile += ".processed.root"; TChain* decodedTestChain = new TChain("TpcTree"); decodedTestChain->AddFile(inputfile); if(nEvents==0) nEvents=decodedTestChain->GetEntries(); std::cout<<"Found "<GetEntries()<<" events in TPC input data file"<SetPersistence(kTRUE); read->SetDatafile(inputfile); read->SetSampleBranchName("TpcSample"); read->SetStartEvent(nEvStart); //read->SetCutSmallPad(); //read->SetMinSamples(1000); fRun->AddTask(read); //################################################################ // Pulse Shape Analysis Task //################################################################ TpcPSATaskALICE* tpsa = new TpcPSATaskALICE(); tpsa->SetPersistence(kTRUE); // tpsa->SetTrivial(kTRUE); tpsa->SetTimeCut(0,90); tpsa->SetPadRowCut(1); // 0 no cut, 1 remove bad rows at beginning and end, 2 remove rows at spacer grids tpsa->SetSampleBranchName("TpcSample"); // Input of PSA tpsa->SetVerbose(kFALSE); fRun->AddTask(tpsa); //################################################################ // Gain Correction Task //################################################################ if(gaincorr==1){ // TString gainfile = basedir + "/tpc/ALICE/par/GainMapAliceGEMTPCIROC_PsTestBeam.csv"; // already configured in parfiles TpcDigiAmpCorrectionTask* tpcCalib = new TpcDigiAmpCorrectionTask(); tpcCalib->SetDigiBranchName("TpcDigi"); // tpcCalib->SetExternalInput(GainFile); tpcCalib->SetVerbose(kFALSE); tpcCalib->SetPersistence(kTRUE); fRun->AddTask(tpcCalib); //mal mit und mal ohne } //################################################################ // Clustering Task //################################################################ //bool SimpleClustering = true; TpcClusterFinderTask* tpcCF = new TpcClusterFinderTask(); tpcCF->SetPersistence(kTRUE); // keep Digis refs in clusters //tpcCF->SetMode(2); // individual timeslice tpcCF->timeslice(21); //in samples tpcCF->SetSingleDigiClusterAmpCut(0); //was 20, then 4 tpcCF->SetClusterAmpCut(0); // cut on mean digi amplitude, was 9.1, then 4 tpcCF->SetErrorPars(600,400); tpcCF->SetSimpleClustering(); // use TpcClusterFinderSimple fRun->AddTask(tpcCF); //################################################################ //Pattern Recognition Tasks //################################################################ // use either Hough PR or Riemann PR, never both!!! // Hough Transform Pattern Recognition //-------------------------------------- configure task TpcSLPatternRecoTaskAlice* tpcSLPR = new TpcSLPatternRecoTaskAlice(); tpcSLPR->SetPersistence(true); tpcSLPR->SetClusterBranchName("TpcCluster"); // tpcSLPR->SetDigiBranchName("TpcDigi"); tpcSLPR->SetTrackBranchName("TpcTrackPreFit"); double parMins[4] = {1.56,0.,1.55,0.}; //TMath::Pi()/2-0.05 TMath::Pi()/2-0.08 {1.54,0.,1.5,0.} double parMaxs[4] = {1.595,3.,1.59,5.}; //TMath::Pi()/2+0.05 TMath::Pi()/2+0.08 {1.595,4.,1.6,20.} tpcSLPR->SetParameterSpace(parMins, parMaxs); tpcSLPR->SetChamberGeometry(85.225, 131.725, -21.4,21.4,0,11); //added by pgadow for alice iroc test, tpcSLPR->SetDepth(5); // MBall 8 tpcSLPR->SetThresh(12); // MBall 15, either 15 (3d) or 30 (2d) tpcSLPR->SetMinCandHits(8); // MBall 15 tpcSLPR->SetClusterAmpCut(0.); // MBall 20 tpcSLPR->SetCutTracksParallelZ(0); //MBall 5 tpcSLPR->SetNClusterLimit(500); //500 or 30 tpcSLPR->SetAbsMomentum(1); //set absolute momentum MBall 1000 tpcSLPR->SetSortMode(1); // tpcSLPR->CutSmallPad(); // tpcSLPR->CutBigPad(); // tpcSLPR->SetSmoothing(); // tpcSLPR->SetDebugOutput(); tpcSLPR->SetVerbose(kFALSE); // tpcSLPR->SetStoreHistograms(outpath+"/HoughHisto12.root"); tpcSLPR->SetWriteHoughTracks(kTRUE); tpcSLPR->SetHoughTrackBranchName("HoughTracks"); // fRun->AddTask(tpcSLPR); /// uncomment if you want to use Hough PR //-------------------------------------- // Riemann Pattern Recognition //-------------------------------------- configure task TpcRiemannTrackingTask* tpcSPR = new TpcRiemannTrackingTask(); tpcSPR->SetClusterBranchName("TpcCluster"); tpcSPR->SetPersistence(kTRUE); tpcSPR->SetSortingParameters(true,0,0); // -1: no sorting, 0: sort Clusters by X, 1: Y, 2: Z, 3: R, 4: distance to interaction point, 5: Phi, -5: -Phi tpcSPR->SetTrkFinderParameters(1.5,0.4,6,1.0); //proximity cut, helix cut, minimum number for points to fit, stretch proximity cut in z direction) tested parameters for 3 d: (4.7,1.1,3,1.2) for 2d: (1.5,0.6,6,1.0) for 20d: (1.,5,0.4,1.) tpcSPR->SetMaxRMS(0.4); //tested parameter for 3d: 0.5 for 2d:0.6 or 20d: 0.4 tpcSPR->SetRiemannScale(80); tpcSPR->SetMultistepParameters(false); tpcSPR->SetMergeTracks(true); tpcSPR->SetMergeCurlers(false); tpcSPR->SetTrkMergerParameters(5,0.8,2.0,1.0);// double TTproxcut, proximity cut in 3Ddouble || TTdipcut, // cut on difference of dip angles of tracklets || double TThelixcut, // distance of the two helices || double TTplanecut);// cut on rms of distances of the riemann hits to intersection of the plane with the sphere of a combined fit // for 3d: (10,0.8,0.5,0.3) for 2d: (5,0.8,2.0,1.0) // tpcSPR->SetVerbose(1); fRun->AddTask(tpcSPR); /// uncomment if you want to use Riemann PR //-------------------------------------- //build GFTracks from TpcRiemannTracks/HoughTracks TpcTrackInitTaskALICE* trackInit=new TpcTrackInitTaskALICE(); trackInit->SetPersistence(); trackInit->SetAllRecoHitPersistence(); trackInit->SetTrackOutBranchName("TpcTrackPreFit"); trackInit->SetClusterBranchName("TpcCluster"); trackInit->SetRecoHitOutBranchName("TpcSPHit"); trackInit->SetPDG(211); // trackInit->useGeane(); // uses RKTrackrep and GeaneTrackrep trackInit->SetSmoothing(smoothing); // trackInit->SetHoughTrackBranchName("HoughTracks"); /// uncomment if you want to use Hough PR // trackInit->SetHough(); /// uncomment if you want to use Hough PR fRun->AddTask(trackInit); // TpcAlignmentTask* align = new TpcAlignmentTask(); // align->SetRecoHitBranchName("TpcSPHit"); // fRun->AddTask(align); //################################################################ //Track Fitting - Kalman Filter Task //################################################################ KalmanTask* kalman = new KalmanTask(); kalman->SetPersistence(kTRUE); // kalman->SetTpcClusterBranchName("TpcCluster"); kalman->SetTpcClusterBranchName("TpcSPHit"); kalman->SetTrackBranchName("TpcTrackPreFit"); kalman->SetOutBranchName("TrackPostFit"); kalman->SetNumIterations(2); // number of fitting iterations (back and forth) // kalman->SetVerbose(); fRun->AddTask(kalman); //################################################################ //dEdx Task //################################################################ // dE/dx task Alice TpcdEdxTaskAlice* dEdx = new TpcdEdxTaskAlice(); dEdx->SetTrackBranchName("TrackPostFit"); dEdx->SetClusterBranchName("TpcSPHit"); // dEdx->SetClusterBranchName("TpcCluster"); //dEdx->SetIdealdEdx(); dEdx->SetPersistence(); //dEdx->SetVerbose(1); dEdx->SetSearchwidth(1); //searchwidth in cm, defines a square in which is searched around the track fRun->AddTask(dEdx); //################################################################ //Residual Task //################################################################ TpcResidualTask* Res = new TpcResidualTask(); Res->SetPersistence(kTRUE); Res->SetNumberOfTrackReps(1); // set to 2 if you use GeaneTrackrep (tpcSPR->useGeane();) if (smoothing==1) Res->SetUnbiased(); fRun->AddTask(Res); // ----- Intialise and run -------------------------------------------- fRun->Init(); // Initialize the Digimapper: TpcDigiPar* tpcpar = FairRun::Instance()->GetRuntimeDb()->getContainer("TpcDigiPar"); // TpcDigiMapper::getInstance()->init(par); TpcDigiMapper::getInstance()->init(tpcpar); TpcDigiMapper::getInstance()->forceManualDriftVel(forceDriftVel); // TpcAlignmentManager::getInstance(); // TpcAlignmentManager::init(tpcpar->getAlignmentFile()); std::cout<<"start:"<Run(nEvStart,nEvStart+nEvents); 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: "<