// runRecoFOPI_batch_matching("/home/pbuehler/physics/projects/fopi/tpc/data/tpc/run_3280.lmd_decoded.root","/home/pbuehler/physics/projects/fopi/tpc/data/fopi/FOPI_","/home/pbuehler/physics/projects/fopi/tpc/results/",false,0,10,0) //FOR E18 machines: //runRecoFOPI_batch_matching("/nfs/hicran/data/tpc/fopi/2011/decoded/run_3280.lmd_decoded.root","/nfs/nas/data/panda/tpc/SIM/fcusanno/FopiFiles/FOPI_","/nfs/nas/data/panda/tpc/SIM/[username]/[FolderName]/",false,0,1,1000,0,0,0) #include <"GFException"> // ----------------------------------------------------------------------------- void runRecoFOPI_batch_matching(TString filename, TString fopiFilePrefix, TString outpath="$RECOOUT", bool forceDriftVel = false, unsigned int smoothing = 0, unsigned int gaincorr = 1, unsigned int nEvents = 0, unsigned int nEvStart = 0, unsigned int outnum = 0) { // ======================================================================== // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug) Int_t iVerbose = 1; TStopwatch timer; timer.Start(); unsigned int nFopiEvents = 0; // Load basic libraries in rootlogon gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C"); rootlogon(); gROOT->SetBatch(); //to ensure running in batch TString basedir = gSystem->Getenv("VMCWORKDIR"); FairRunAna* fRun = new FairRunAna(); //FairRunSim* fSim = new FairRunSim(); TString jobdir = outpath; std::string jobname(filename.Data()); int last = jobname.rfind("/"); if(last>0) jobname = jobname.substr(last+1,jobname.size()+1); TString outName(jobname); if(outName.Contains("repaired.")) outName.ReplaceAll(".lmd_decoded_repaired.root",".reco.root"); if(outName.Contains("decoded.")) outName.ReplaceAll(".lmd_decoded.root",".reco.root"); if(smoothing) outName.ReplaceAll("reco.root", "smoothed_reco.root"); if(outnum>0) outName.ReplaceAll(".root",Form("%i.root",outnum)); TString outFile = outpath+"/"; outFile += outName; TString fopiSuffix = outName; if(outnum==0){ fopiSuffix.ReplaceAll(".reco.root",".root"); } else{ fopiSuffix.ReplaceAll(Form(".reco%i.root",outnum),".root"); } std::string suffixname(fopiSuffix.Data()); int suffl = suffixname.rfind("_"); if(suffl>0) suffixname = suffixname.substr(suffl+1,suffixname.size()+1); TString fopiFilename(suffixname); TString runName = fopiFilename; fopiFilename = fopiFilePrefix + fopiFilename; // std::cout<<"FOPI input file: "<Setenv("PROUTFILENAME", PROutFile.Data()); gROOT->ProcessLine(".! rm $PROUTFILENAME"); gSystem->Unsetenv("PROUTFILENAME"); } fRun->SetOutputFile(outFile); //--- Parameter database --- FairRuntimeDb* rtdb = fRun->GetRuntimeDb(); FairParAsciiFileIo* parInput1 = new FairParAsciiFileIo(); TString parFile = "/tpc/FOPI/par/tpc.run"+runName+".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(); // build a run specific geometry // run the FopiGeom.C in a separate root session // // taroff: distance between target and center of magnet // targetpos: distance between real target position and nominal FOPI target // position // In experiment S339 the distance between the nominal FOPI target position // and the center of the magnet is 41 cm. Double_t taroff = -65; Double_t targetpos = taroff+41.; // possible materials: copper (default), lead, carbon TString targetmat = TString("carbon"); TString alignmentFileName = par->getAlignmentFile(); TString geoFile; geoFile.Form("$FOPI2ROOT/fopigeometry/FopiGeom_s339%5.1f.root",targetpos); TString cmd; cmd.Form(".! root -b -q '$FOPI2ROOT/fopigeometry/FopiGeom.C(%f,\"%s\",\"%s\",\"%s\",kFALSE)'",targetpos,targetmat.Data(),alignmentFileName.Data(),geoFile.Data()); gROOT->ProcessLine(cmd); // add geometry file to fRun fRun->SetGeomFile(geoFile); // FOPI magnetic field //For pion beam data, FOPIField(0.616) has to be used according to CDC/RPC analysis! // if the constructor is called with no argument, 0.6T is used FOPIField *fMagField = new FOPIField(0.616); fMagField->SetTargetOffset(taroff); fRun->SetField(fMagField); //extract number of entries in external data tree TFile testFile(filename); unsigned int numEvents = ((TTree*)testFile.Get("tpcEvent"))->GetEntries(); if(nEvents==0 || nEvents>numEvents) nEvents=numEvents; TFile testFile2(fopiFilename); unsigned int numFopiEvents = ((TTree*)testFile2.Get("FOPITree"))->GetEntries(); if(nFopiEvents==0 || nFopiEvents>numFopiEvents) {nFopiEvents=numFopiEvents;} std::cout<<"Found "<AddTask(Geane); TpcDataReaderTask* read = new TpcDataReaderTask(); read->SetPersistence(kFALSE); read->SetDatafile(filename); read->SetSampleBranchName("TpcSample"); read->SetStartEvent(nEvStart); //read->SetCutSmallPad(); //read->SetMinSamples(1000); fRun->AddTask(read); FopiDataReaderTask* readFopi = new FopiDataReaderTask(); readFopi->SetPersistence(kTRUE); //readFopi->SetBarrelOn(); //not (re)implemented yet, but it is only a double-check //readFopi->SetRunNr(runNr); readFopi->SetInputFile(fopiFilename); //readFopi->SetNevent(nFopiEvents); fRun->AddTask(readFopi); TpcPSATask* tpsa = new TpcPSATask(); tpsa->SetPersistence(kFALSE); tpsa->SetSampleBranchName("TpcSample"); // Input of PSA fRun->AddTask(tpsa); if(gaincorr==1){ TString gainfile = basedir; TpcDigiAmpCorrectionTask* tpcCalib = new TpcDigiAmpCorrectionTask(); tpcCalib->SetVerbose(kFALSE); tpcCalib->SetPersistence(kFALSE); fRun->AddTask(tpcCalib); } //bool SimpleClustering = true; TpcClusterFinderTask* tpcCF = new TpcClusterFinderTask(); tpcCF->SetPersistence(kFALSE); // keep Digis refs in clusters //tpcCF->SetMode(2); // individual timeslice tpcCF->timeslice(6); //in samples tpcCF->SetSingleDigiClusterAmpCut(20); tpcCF->SetClusterAmpCut(9.1); // cut on mean digi amplitude tpcCF->SetErrorPars(600,400); tpcCF->SetSimpleClustering(); // use TpcClusterFinderSimple fRun->AddTask(tpcCF); //find TpcRiemannTracks in the TPC alone TpcRiemannTrackingTask* tpcSPR = new TpcRiemannTrackingTask(); tpcSPR->SetClusterBranchName("TpcCluster"); tpcSPR->SetPersistence(kFALSE); //tpcSPR->SetVerbose(1); fRun->AddTask(tpcSPR); //build GFTracks from TpcRiemannTracks TpcTrackInitTask* trackInit=new TpcTrackInitTask(); trackInit->SetPersistence(kFALSE); trackInit->SetOutBranchName("TpcTrackPreFit"); //trackInit->SetVerbose(1); //trackInit->SetMCPid(); // use ideal particle identification trackInit->SetClusterBranchName("TpcCluster"); trackInit->SetPDG(211); // trackInit->useGeane(); // uses RKTrackrep and GeaneTrackrep trackInit->SetSmoothing(true); fRun->AddTask(trackInit); TpcSLPatternRecoTask* tpcSLPR = new TpcSLPatternRecoTask(); tpcSLPR->SetPersistence(kFALSE); tpcSLPR->SetStoreHistograms(PROutFile); tpcSLPR->SetClusterAmpCut(20.); tpcSLPR->SetCutTracksParallelZ(5); //tpcSLPR->SetXSorting(true); double parMins[4] = {-TMath::Pi(),0.,-TMath::Pi(),0.}; double parMaxs[4] = {TMath::Pi(),10.,TMath::Pi(),20.}; tpcSLPR->SetParameterSpace(parMins, parMaxs); tpcSLPR->SetDepth(8); tpcSLPR->SetThresh(15); tpcSLPR->SetMinCandHits(15); //tpcSLPR->SetClusterBranchName("TpcCluster_cut"); tpcSLPR->SetAbsMomentum(1000); //fRun->AddTask(tpcSLPR); KalmanTask* kalman = new KalmanTask(); kalman->SetPersistence(kTRUE); kalman->SetTpcClusterBranchName("TpcCluster"); kalman->SetTrackBranchName("TpcTrackPreFit"); kalman->SetOutBranchName("TpcTrackPostFit"); kalman->SetNumIterations(3); // number of fitting iterations (back and forth) fRun->AddTask(kalman); //build GFTracks for CDC hits only FopiTrackInitTask* trackCdcInit=new FopiTrackInitTask(); trackCdcInit->SetPersistence(kFALSE); //trackCdcInit->SetHitBranchName("CdcHit"); //trackCdcInit->SetTrackBranchName("CdcTrack"); //trackCdcInit->SetVerbose(1); //trackCdcInit->SetPiKaonThreshold(0.32); //trackCdcInit->SetMinCDCHits(5); //trackCdcInit->SetSmoothing(true); fRun->AddTask(trackCdcInit); KalmanTask* KalmanCdc = new KalmanTask(); KalmanCdc->SetPersistence(kTRUE); KalmanCdc->SetTrackBranchName("CdcTrackPreFit"); KalmanCdc->SetOutBranchName("CdcTrackPostFit"); KalmanCdc->SetFopiCDCBranchName("CdcHit"); KalmanCdc->SetNumIterations(3); // number of fitting iterations fRun->AddTask(KalmanCdc); //match TPC and CDC/RPC tracks TpcCdcMatchingTask* tpcCdcMatch = new TpcCdcMatchingTask(); //tpcCdcMatch->SetPersistence(kTRUE); //tpcCdcMatch->SetTpcTrackBranchName("TpcTrackPostFit"); //tpcCdcMatch->SetTpcClusterBranchName("TpcCluster"); //tpcCdcMatch->SetCdcTrackBranchName("CdcTrack"); //tpcCdcMatch->SetCdcHitBranchName("CdcHit"); //tpcCdcMatch->SetOutBranchName("TpcCdcPreFit"); tpcCdcMatch->SetMaxMatDistance(.5); //tpcCdcMatch->SetMaxMatPhi(acos(0.)); //tpcCdcMatch->SetMinHitsperLength(0.0); //tpcCdcMatch->SetMinClus(30); fRun->AddTask(tpcCdcMatch); KalmanTask* KalmanComb = new KalmanTask(); KalmanComb->SetPersistence(kTRUE); KalmanComb->SetTpcClusterBranchName("TpcCluster"); KalmanComb->SetFopiCDCBranchName("CdcHit"); KalmanComb->SetTrackBranchName("TpcCdcPreFit"); KalmanComb->SetOutBranchName("TpcCdcPostFit"); KalmanComb->SetNumIterations(3); // number of fitting iterations fRun->AddTask(KalmanComb); // build GFTracks for Helitron/Plawa hits only HeliPlawaTrackInitTask* trackHeliPlawaInit = new HeliPlawaTrackInitTask(); trackHeliPlawaInit->SetPersistence(kFALSE); trackHeliPlawaInit->SetVerbose(kFALSE); trackHeliPlawaInit->setWithPlawaHit(1); trackHeliPlawaInit->SetPosErr(5.,5.,0.1); trackHeliPlawaInit->SetMomErr(0.2,0.2,0.2); //trackHeliPlawaInit->SetHHitRes(0.1,1.0,0.1); trackHeliPlawaInit->SetdrMax(100.0); trackHeliPlawaInit->SetResScale(1.0); trackHeliPlawaInit->SetPHitRes(10.0,10.0,1.0); //fRun->AddTask(trackHeliPlawaInit); HeliKalmanTask* KalmanHeliPlawa = new HeliKalmanTask(); KalmanHeliPlawa->SetPersistence(kFALSE); KalmanHeliPlawa->SetVerbose(kFALSE); KalmanHeliPlawa->SetTrackBranchName("HeliPlawaTrackPreFit"); KalmanHeliPlawa->SetOutBranchName("HeliPlawaTrackPostFit"); KalmanHeliPlawa->SetNumIterations(3); // number of fitting iterations //fRun->AddTask(KalmanHeliPlawa); // match TPC and Helitron tracks TpcHelMatchingTask* tpcHelMatch = new TpcHelMatchingTask(); tpcHelMatch->SetPersistence(kFALSE); tpcHelMatch->SetVerbose(kFALSE); tpcHelMatch->SetTpcTrackFitBranchName(TString("TpcTrackPostFit")); // tpcHelMatch->SetTpcTrackFitBranchName(TString("TpcCdcPostFit")); tpcHelMatch->setNHeliHits(100); tpcHelMatch->setWithPlawaHit(0); tpcHelMatch->SetPosErr(0.1,0.1,0.1); tpcHelMatch->SetMomErr(0.2,0.2,0.2); // tpcHelMatch->SetHHitRes(0.1,1.0,0.1); tpcHelMatch->SetdrMax(100.0); tpcHelMatch->SetResScale(0.9); tpcHelMatch->SetminHPmatch(0.5); tpcHelMatch->SetPHitRes(5.,5.,5.); fRun->AddTask(tpcHelMatch); // fit matched TPC and Helitron tracks TpcHelTrackFittingTask* tpcHelTrackFit = new TpcHelTrackFittingTask(); tpcHelTrackFit->SetPersistence(kFALSE); tpcHelTrackFit->SetVerbose(kFALSE); tpcHelTrackFit->SetDefRep(0); fRun->AddTask(tpcHelTrackFit); // display tpc and Helitron/Plawa tracks together /* DO NOT USE IT FOR PRODUCTION! FOLLOWING TASK CAN BE ACTIVATED FOR FEW EVENTS. IN THIS CASE, REMEMBER TO REMOVE COMMENTS TO THE OPENING AND CLOSING OF THE RELATED PDF FILE*/ /* TpcHelDisplayTask* tpcHelDisp = new TpcHelDisplayTask(); tpcHelDisp->SetVerbose(kTRUE); tpcHelDisp->SetOutName(outpath); tpcHelDisp->SetDisplayType(1); tpcHelDisp->SetMaxChi2(1500.); // tt,hf,ht, v, t, c, h, p,th tpcHelDisp->SetToShow( 1, 1, 1, 1, 1, 1, 1, 1, 1); tpcHelDisp->SetToShow( 1, 0, 1, 1, 0, 1, 0, 1, 1); fRun->AddTask(tpcHelDisp); */ /* // Philipps tasks ----------------------------------------------------------- // match TPC and Helitron tracks TpcHelMatchingTaskPhilipp* tpcHelMatch = new TpcHelMatchingTaskPhilipp(); tpcHelMatch->SetPersistence(kTRUE); tpcHelMatch->SetVerbose(kTRUE); tpcHelMatch->SetTpcTrackFitBranchName(TString("TpcTrackPostFit")); tpcHelMatch->SetPosErr(0.1,0.1,0.1); tpcHelMatch->SetMomErr(0.2,0.2,0.2); //tpcHelMatch->SetHHitRes(0.1,1.0,0.1); tpcHelMatch->SetResScale(1.0); tpcHelMatch->SetPHitRes(10.,10.,1.); fRun->AddTask(tpcHelMatch); // fit matched TPC and Helitron tracks TpcHelTrackFittingTaskPhilipp* tpcHelTrackFit = new TpcHelTrackFittingTaskPhilipp(); tpcHelTrackFit->SetPersistence(kTRUE); tpcHelTrackFit->SetVerbose(kTRUE); // fRun->AddTask(tpcHelTrackFit); // display tpc and Helitron/Plawa tracks together TpcHelDisplayTaskPhilipp* tpcHelDisp = new TpcHelDisplayTaskPhilipp(); tpcHelDisp->SetVerbose(kTRUE); tpcHelDisp->SetOutPath(outpath); tpcHelDisp->SetDisplayType(1); // tt,hf,ht, v, t, c, h, p,th tpcHelDisp->SetToShow( 1, 1, 1, 1, 1, 1, 1, 1, 1); tpcHelDisp->SetToShow( 1, 0, 1, 1, 0, 1, 0, 1, 1); // fRun->AddTask(tpcHelDisp); // Philipps tasks ----------------------------------------------------------- */ // open pdf-file for TpcHelDisplayTask // TO BE USED WITH TpcHelDisplayTask TASK // TCanvas * c1 = new TCanvas(); //c1->Print(outpath+"/TpcHel_output.pdf["); // dE/dx task TpcdEdxTask* dEdx = new TpcdEdxTask(); dEdx->SetTrackBranchName("TpcTrackPostFit"); dEdx->SetClusterBranchName("TpcCluster"); //dEdx->SetIdealdEdx(); dEdx->SetPersistence(); dEdx->SetDXgrid(0.5); fRun->AddTask(dEdx); // TpcResidualTask* Res = new TpcResidualTask(); // Res->SetPersistence(); // //Res->SetNumberOfTrackReps(1); // set to 2 if you use GeaneTrackrep (tpcSPR->useGeane();) // //SLres->SetClusterBranchName("TpcCluster_cut"); // 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,nEvents); // close pdf-file // TO BE USED WITH TpcHelDisplayTask TASK //c1->Print(outpath+"/TpcHel_output.pdf]"); 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: "<