// ----------------------------------------------------------------------------- void runRecoFOPI_batch_matching_Felix(int runNumber, TString jobname, TString outpath="$RECOOUT", bool forceDriftVel = false, unsigned int deepPersist = 0, //persist Digis, clusters unsigned int smoothing = 0, unsigned int gaincorr = 1, unsigned int nEvents = 0, unsigned int nEvStart = 0, TString decodedInPath="/nfs/hicran/data/tpc/fopi/2011/decoded/", //TString fopiRecoPath="/nfs/nas/data/panda/tpc/SIM/fcusanno/FopiFiles/", TString fopiRecoPath="/nfs/mds/data/tpc/fopi/2011/fopi_reco_sverre/rev176/", TString SGE_IDENTIFIER="" ) { // ======================================================================== // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug) // Load basic libraries in rootlogon gROOT->LoadMacro("$VMCWORKDIR/gconfig/rootlogon.C"); rootlogon(); Int_t iVerbose = 1; TStopwatch timer; timer.Start(); unsigned int nFopiEvents = 0; gROOT->SetBatch(); //to ensure running in batch TString basedir = gSystem->Getenv("VMCWORKDIR"); FairRunAna* fRun = new FairRunAna(); if(smoothing) jobname += "_smoothed.reco.root"; else jobname += ".reco.root"; TString outFile = outpath+"/"; outFile += jobname; TString zCorrFile = outpath+"/analysis/zShifts.root"; //std::cout<SetOutputFile(outFile); std::stringstream ss; ss<GetRuntimeDb(); FairParAsciiFileIo* parInput1 = new FairParAsciiFileIo(); TString parFile = "/tpc/FOPI/par/tpc.run"+runNr+".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(); //root file parameters TString copy = jobname; TString parfile = outFile.ReplaceAll("reco.root","param.root"); FairParRootFileIo* parOut=new FairParRootFileIo(kTRUE); parOut->open(parfile.Data()); rtdb->setOutput(parOut); //TpcAlignmentManager::init(par->getAlignmentFile()); TString alignFile = basedir; alignFile += "/tpc/FOPI/par/alignmentAfter_4.txt"; TpcAlignmentManager::init(alignFile); TpcAlignmentManager::getInstance(); TpcDigiMapper::getInstance()->init(par); TpcDigiMapper::getInstance()->forceManualDriftVel(forceDriftVel); double vdrift=TpcDigiMapper::getInstance()->getGas()->VDrift(); //################################################################################ //################################################################################ // 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 = -61.; //FIXED from 65->61 on 30.7.2014 F.B Double_t targetpos = taroff+41.; // possible materials: copper (default), lead, carbon TString targetmat = TString("carbon"); TString alignmentFileName = par->getAlignmentFile(); TString geoFile; //IF we run parallel, each job has to write his OWN geometry file to avoid conflicts : if(SGE_IDENTIFIER.Length()>0) geoFile.Form("/tmp/FopiGeom_s339%5.1f_%s_%s.root", targetpos,targetmat.Data(),SGE_IDENTIFIER.Data()); else geoFile.Form("$FOPI2ROOT/fopigeometry/FopiGeom_s339%5.1f_%s.root",targetpos,targetmat.Data()); TString cmd; cmd.Form(".! root -b -q '$FOPI2ROOT/fopigeometry/FopiGeom.C(%f,\"%s\",\"%s\",\"%s\",kFALSE)'",targetpos,targetmat.Data(),alignmentFileName.Data(),geoFile.Data()); std::cout<ProcessLine(cmd); std::cout<<"################################################################################"<GetGeoFile() == NULL) { std::cout<<"\n trying to set geometry file... try #"<SetGeomFile(geoFile); ntries++; } // 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); //const field: //PndConstField *fMagField=new PndConstField(); //fMagField->SetField(0., 0., 6.16 ); // values are in kG //fMagField->SetFieldRegion(-500, 500,-500, 500, -2000, 2000); // values are in cm fRun->SetField(fMagField); //initialize GF field GFFieldManager::getInstance()->init(new PndFieldAdaptor(fRun->GetField())); GFMaterialEffects::getInstance()->init(new GFTGeoMaterialInterface()); //GFMaterialEffects::getInstance()->setNoEffects(); //parse input files: TString pythonScript=basedir+"/macro/tpc/FOPI/parseFiles.py"; TPython::LoadMacro(pythonScript); PyFileParser fp; if(decodedInPath.Last('/') != decodedInPath.Length()-1) decodedInPath+="/"; fp.setFiles(TString(decodedInPath+"*"+runNr+"*").Data()); int nFiles = (int)fp.getNum(); std::cout<<"FOUND "<AddFile((char*)fp.getFile(i)); std::cout<<"Add file "<<(char*)fp.getFile(i)<<" to decoded chain"<1) Fatal("MACRO:","More than 1 decoded file found in input path"); if(fopiRecoPath.Last('/') != fopiRecoPath.Length()-1) fopiRecoPath+="/"; fp.setFiles(TString(fopiRecoPath+"FOPI_*"+runNr+"*").Data()); nFiles = (int)fp.getNum(); TChain* recoTestChain = new TChain("FOPITree"); for(int i=0; iAddFile((char*)fp.getFile(i)); std::cout<<"Add file "<<(char*)fp.getFile(i)<<" to reco chain"<1) Fatal("MACRO:","More than 1 reco file found in input path"); if(nEvents==0) nEvents=decodedTestChain->GetEntries(); std::cout<<"Found "<GetEntries()<<" events in TPC input data file"<GetEntries()<<" events in fopi input data file"<SetPersistence(kFALSE); read->SetDatafile(decodedTestChain->GetFile()->GetName()); read->SetSampleBranchName("TpcSample"); read->SetStartEvent(nEvStart); //read->SetCutSmallPad(); //read->SetMinSamples(1000); fRun->AddTask(read); FopiDataReaderTask* readFopi = new FopiDataReaderTask(); readFopi->SetPersistence(false); readFopi->SetCdcHitPersistence(true); readFopi->SetCdcTrackPersistence(true); readFopi->SetRpcTrackPersistence(true); //readFopi->SetBarrelOn(); //not (re)implemented yet, but it is only a double-check //readFopi->SetRunNr(runNr); readFopi->SetInputFile(recoTestChain->GetFile()->GetName()); //readFopi->SetNevent(nFopiEvents); readFopi->SetDisableHeliTracks(true); readFopi->SetBarrelOn(); fRun->AddTask(readFopi); //CORRECT FOR CDC SHIFTS CdcCorrectionTask* corTask=new CdcCorrectionTask(); corTask->SetCorrection(0.18); //default is 0.2 //corTask->SetCorrectionFile("testShifts.root"); corTask->SetCorrectionFile(zCorrFile); fRun->AddTask(corTask); //Conversion from fopi2root RpcTrack objects to GF strip hits FopiRpcInitTask* rpcInit = new FopiRpcInitTask(); rpcInit->SetPersistence(true); 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(true); 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(DP); tpsa->SetSamplePersistence(DP); tpsa->SetSampleBranchName("TpcSample"); // Input of PSA //tpsa->UseNoiseMaxSampleCut(decodedTestChain->GetFile()->GetName(),7.); fRun->AddTask(tpsa); if(gaincorr==1){ TString gainfile = basedir; TpcDigiAmpCorrectionTask* tpcCalib = new TpcDigiAmpCorrectionTask(); tpcCalib->SetVerbose(kFALSE); tpcCalib->SetPersistence(DP); fRun->AddTask(tpcCalib); } //bool SimpleClustering = true; TpcClusterFinderTask* tpcCF = new TpcClusterFinderTask(); tpcCF->SetPersistence(DP); // keep Digis refs in clusters tpcCF->SetDigiPersistence(DP); //tpcCF->SetMode(2); // individual timeslice tpcCF->timeslice(6); //in samples //tpcCF->SetSingleDigiClusterAmpCut(10); //tpcCF->SetClusterAmpCut(5); // cut on mean digi amplitude tpcCF->SetErrorPars(-1,10.); tpcCF->SetSimpleClustering(); // use TpcClusterFinderSimple fRun->AddTask(tpcCF); //find TpcRiemannTracks in the TPC alone TpcRiemannTrackingTask* tpcSPR = new TpcRiemannTrackingTask(); tpcSPR->SetClusterBranchName("TpcCluster"); tpcSPR->SetSortingParameters(true,3,42.78); tpcSPR->SetPersistence(DP); //tpcSPR->SetVerbose(1); fRun->AddTask(tpcSPR); //build GFTracks from TpcRiemannTracks TpcTrackInitTask* trackInit=new TpcTrackInitTask(); trackInit->SetPersistence(DP); trackInit->SetAllRecoHitPersistence(true); trackInit->SetTrackOutBranchName("TpcTrackPreFit"); trackInit->SetClusterBranchName("TpcCluster"); trackInit->SetRecoHitOutBranchName("TpcSPHit"); TVector3 IP(0,0,42.78); trackInit->SetInteractionPoint(IP); trackInit->SetPDG(211); trackInit->SetSmoothing(smoothing); trackInit->SetWeightedPlaneConstruction(true); trackInit->SetCutCov(true); fRun->AddTask(trackInit); TpcAlignmentTask* align = new TpcAlignmentTask(); align->SetRecoHitBranchName("TpcSPHit"); fRun->AddTask(align); //TPC standalone fitting. I guess we don't really need to do that... KalmanTask* kalman = new KalmanTask(); //kalman->SetUseDAF(); kalman->SetPersistence(true); //kalman->SetTpcClusterBranchName("TpcCluster"); kalman->SetTpcClusterBranchName("TpcSPHit"); kalman->SetTrackBranchName("TpcTrackPreFit"); kalman->SetOutBranchName("TpcTrackPostFit"); kalman->SetNumIterations(1); // number of fitting iterations (back and forth) fRun->AddTask(kalman); //build GFTracks for CDC hits only. NO BARREL/RPC FopiTrackInitTask* trackCdcInit=new FopiTrackInitTask(); trackCdcInit->SetPersistence(false); //affects prefit tracks trackCdcInit->SetAllRecoHitPersistence(true); //trackCdcInit->SetTrackPersistence(true); trackCdcInit->SetHitBranchName("CdcHit"); trackCdcInit->SetTrackBranchName("CdcTrack"); trackCdcInit->SetOutBranchName("CdcTrackPreTempFit"); //still mising BAR/RPC trackCdcInit->SetRecoHitOutBranchName("CdcSpacepoint"); trackCdcInit->SetTrackPersistence(DP); //trackCdcInit->SetPiKaonThreshold(0.32); //trackCdcInit->SetMinCDCHits(5); trackCdcInit->SetSmoothing(smoothing); //trackCdcInit->SetUsePseudoSPHits(); trackCdcInit->SetWeightedPlaneConstruction(true); trackCdcInit->SetCutCov(true); //in cm: trackCdcInit->SetConstantPointErrors(3.5e-2, 3.5e-2, 10.); fRun->AddTask(trackCdcInit); //refit with added BAR / RPC hits KalmanTask* KalmanCdc1 = new KalmanTask(); //KalmanCdc1->SetUseDAF(); KalmanCdc1->SetPersistence(false); KalmanCdc1->SetTpcClusterBranchName("TpcSPHit"); KalmanCdc1->SetFopiCDCBranchName("CdcSpacepoint"); KalmanCdc1->SetTrackBranchName("CdcTrackPreTempFit"); KalmanCdc1->SetOutBranchName("CdcTrackPostTempFit"); KalmanCdc1->SetNumIterations(1); // number of fitting iterations fRun->AddTask(KalmanCdc1); //match TPC and CDC tracks (no GFTracks!) TpcCdcMatchingTask* tpcCdcMatch = new TpcCdcMatchingTask(); tpcCdcMatch->SetPersistence(kTRUE); tpcCdcMatch->SetTpcTrackBranchName("TpcTrackPostFit"); tpcCdcMatch->SetCdcTrackBranchName("CdcTrack"); tpcCdcMatch->SetCdcHitBranchName("CdcHit"); tpcCdcMatch->SetWriteFailedMatches(); fRun->AddTask(tpcCdcMatch); //merge CDC and TPC FopiTrackMergingTask* merge = new FopiTrackMergingTask(); merge->SetPersistence(); merge->SetTrackPersistence(false); //prefit tracks merge->SetTpcTrackBranchName("TpcTrackPostFit"); merge->AddTupleBranchName("TpcCdcMatchingPair"); //merge->SetTrackOutBranchName("FopiTrackPreFit"); //TEST merge->SetTrackOutBranchName("FopiTrackPreTempFit"); merge->SetCdcGFTrackBranchName("CdcTrackPreTempFit"); fRun->AddTask(merge); //obtain fits for TPC-CDC tracks. Needed for RPC/BAR matching decision KalmanTask* KalmanComb = new KalmanTask(); KalmanComb->SetUseDAF(); KalmanComb->SetPersistence(true); KalmanComb->SetTpcClusterBranchName("TpcSPHit"); KalmanComb->SetFopiCDCBranchName("CdcSpacepoint"); KalmanComb->SetTrackBranchName("FopiTrackPreTempFit"); KalmanComb->SetOutBranchName("FopiTrackPostTempFit"); //KalmanComb->SetNumIterations(1); // number of fitting iterations fRun->AddTask(KalmanComb); //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(true); addBarRpc->SetSmoothing(smoothing); addBarRpc->SetTrackBranchName("FopiTrackPostTempFit"); 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.); fRun->AddTask(addBarRpc); //refit with added BAR / RPC hits KalmanTask* KalmanCdc = new KalmanTask(); KalmanCdc->SetUseDAF(); KalmanCdc->SetPersistence(kTRUE); KalmanCdc->SetTpcClusterBranchName("TpcSPHit"); KalmanCdc->SetFopiCDCBranchName("CdcSpacepoint"); KalmanCdc->SetTrackBranchName("CdcTrackPreFit"); KalmanCdc->SetOutBranchName("CdcTrackPostFit"); //KalmanCdc->SetNumIterations(2); // number of fitting iterations fRun->AddTask(KalmanCdc); Chi2CleanupTask* cleanupCDC = new Chi2CleanupTask(); cleanupCDC->SetTrackBranchName("CdcTrackPostFit"); //cleanupCDC->AddDiffCutWindow(-1.e6,-10.); cleanupCDC->AddDiffCutWindow(10.,1.e6); cleanupCDC->SetHitWindow(0.93,1.); //fRun->AddTask(cleanupCDC); 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(kTRUE); 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); Chi2CleanupTask* cleanupCOMB = new Chi2CleanupTask(); cleanupCOMB->SetTrackBranchName("FopiTrackPostFit"); //cleanupCOMB->AddDiffCutWindow(-1.e6,-10.); cleanupCOMB->AddDiffCutWindow(10.,1.e6); cleanupCOMB->SetHitWindow(0.95,1.); //fRun->AddTask(cleanupCOMB); // dE/dx task TpcdEdxTask* dEdx = new TpcdEdxTask(); dEdx->SetTrackBranchName("FopiTrackPostFit"); dEdx->SetClusterBranchName("TpcSPHit"); //dEdx->SetIdealdEdx(); dEdx->SetPersistence(); //dEdx->SetVerbose(1); dEdx->SetUseFirstDigiPos(); dEdx->SetIgnoreEdgeDigis(5.8,14.5); dEdx->SetDXgrid(0.5); fRun->AddTask(dEdx); //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(decodedTestChain->GetFile()->GetName()); fRun->AddTask(qa); // ----- Intialise and run -------------------------------------------- fRun->Init(); std::cout<<"\n****** gGeoManager->Print() ***************** \n"<Print(); std::cout<<"\n*********************************************\n\n"<saveOutput(); rtdb->print(); std::cout<<"start:"<Run(nEvStart,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: "<