void runRecoFOPI_AlignmentTest(int runNumber, bool field, unsigned int smoothing = 0, unsigned int seed=0, int jobset=0) { // ======================================================================== // Verbosity level (0=quiet, 1=event level, 2=track level, 3=debug) Int_t iVerbose = 0; gRandom->SetSeed(seed); gROOT->ProcessLine(".x $VMCWORKDIR/gconfig/rootlogon.C"); // Input file TString sysFile = gSystem->Getenv("VMCWORKDIR"); // Parameter file TString inSimFile; if(jobset!=0){ inSimFile=Form("AlignmentTest_set%i_run%i.mc.root",jobset,runNumber); }else{ inSimFile=Form("AlignmentTest_%i.mc.root",runNumber); } cout<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); // const Double_t kAu2Gev=0.9314943228; // const Double_t khSlash = 1.0545726663e-27; // const Double_t kErg2Gev = 1/1.6021773349e-3; // const Double_t khShGev = khSla\sh*kErg2Gev; // const Double_t kYear2Sec = 3600*24*365.25; // TDatabasePDG::Instance()->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE, // khShGev/(12.33*kYear2Sec),3,"Ion",1000010030); // FairLogger::GetLogger()->SetLogScreenLevel("DEBUG2"); // ----- 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: "<GetRuntimeDb(); FairParRootFileIo* parInput1 = new FairParRootFileIo(); parInput1->open(parFile.Data()); parInput1->print(); //FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo(); // parIo1->open(allDigiFile.Data(),"in"); rtdb->setFirstInput(parInput1); // rtdb->setSecondInput(parIo1); //rtdb->readAll(); // rtdb->initContainers(fRun->GetRunId()); // r->print(); //parInput1->getRun(); TpcDigiPar* par = (TpcDigiPar*) rtdb->getContainer("TpcDigiPar"); // // par->setInputVersion(fRun->GetRunId(),1); // par->setChanged(kTRUE); // std::cout<<"################################################################################"<init(); // std::cout<<"################################################################################"<SetField(0., 0. , 0. ); // values are in kG fMagField->SetFieldRegion(-50, 50,-50, 50, -2000, 2000); // values are in cm fRun->SetField(fMagField); }else{ FOPIField *fMagField = new FOPIField(0.616); fMagField->SetTargetOffset(taroff); fRun->SetField(fMagField); } // ------- RECO procedure ------------------------------------------------ TpcEventCounter* evCount = new TpcEventCounter(); // evCount->SetnEvts(nEvents); evCount->SetStep(2); fRun->AddTask(evCount); TpcClusterFinderTask* tpcCF = new TpcClusterFinderTask(); // tpcCF->SetDigiPersistence(); // keep reference to digis in clusters // tpcCF->SetPersistence(kFALSE); // keep Clusters tpcCF->SetPersistence(kTRUE); // keep Clusters tpcCF->timeslice(6); //in samples tpcCF->SetThreshold(1); tpcCF->SetSingleDigiClusterAmpCut(4.); tpcCF->SetClusterAmpCut(2.); // cut on mean digi amplitude tpcCF->SetErrorPars(-1.,22.); tpcCF->SetSimpleClustering(); // use TpcClusterFinderSimple fRun->AddTask(tpcCF); //find TpcRiemannTracks in the TPC alone TpcRiemannTrackingTask* tpcSPR = new TpcRiemannTrackingTask(); tpcSPR->SetClusterBranchName("TpcCluster"); tpcSPR->SetPersistence(kFALSE); //tpcSPR->SetSortingParameters(true,4,0); //tpcSPR->SetVerbose(1); fRun->AddTask(tpcSPR); //build GFTracks from TpcRiemannTracks TpcTrackInitTask* trackInit=new TpcTrackInitTask(); trackInit->SetClusterBranchName("TpcCluster"); trackInit->SetTrackOutBranchName("TpcTrackPreFit"); trackInit->SetRecoHitOutBranchName("TpcSPHit"); trackInit->SetPersistence(); trackInit->SetAllRecoHitPersistence(); trackInit->SetMCPid(); // use ideal particle identification trackInit->SetVerbose(0); //trackInit->SetPDG(211); //trackInit->useGeane(); // uses RKTrackrep and GeaneTrackrep if (smoothing==1) trackInit->SetSmoothing(true); fRun->AddTask(trackInit); TpcAlignmentTask* align = new TpcAlignmentTask(); align->SetRecoHitBranchName("TpcSPHit"); fRun->AddTask(align); KalmanTask* kalman =new KalmanTask(); kalman->SetUseDAF(); kalman->SetTpcClusterBranchName("TpcSPHit"); kalman->SetTrackBranchName("TpcTrackPreFit"); kalman->SetOutBranchName("TpcTrackPostFit"); kalman->SetPersistence(); kalman->SetNumIterations(3); // number of fitting iterations (back and forth) fRun->AddTask(kalman); TpcShieldToCdcTrackTask *CdcTrackTask = new TpcShieldToCdcTrackTask(); CdcTrackTask->SetTrackOutBranchName("CdcTrackPostFit"); CdcTrackTask->SetPersistence(true); fRun->AddTask(CdcTrackTask); //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("CdcTrackPostFit"); fRun->AddTask(merge); // MCTrackResCalc* mcTrackResCalc = new MCTrackResCalc(); // mcTrackResCalc->addBranchName("TpcTrackPostFit"); // mcTrackResCalc->addBranchName("MCTrack"); // mcTrackResCalc->setDebug(); // mcTrackResCalc->setExtrapolateToMcStartPos(); // TpcMcTrackResidualTask* resTask = new TpcMcTrackResidualTask(); // resTask->SetResCalculator(mcTrackResCalc); // resTask->SetPersistence(); // fRun->AddTask(resTask); // ----- Intialise and run -------------------------------------------- std::cerr<<"1"<getContainer("TpcDigiPar"); // par->setInputVersion(fRun->GetRunId(),1); // par->setChanged(kTRUE); // par->init(); GFFieldManager::getInstance()->init(new PndFieldAdaptor(fRun->GetField())); GFMaterialEffects::getInstance()->init(new GFTGeoMaterialInterface()); fRun->Init(); TpcDigiPar* par = (TpcDigiPar*)FairRun::Instance()->GetRuntimeDb()->getContainer("TpcDigiPar"); // par->printParams(); TpcAlignmentManager::init(par->getAlignmentFile()); TpcDigiMapper::getInstance()->init(par); //TpcDigiMapper::getInstance()->forceManualDriftVel(forceDriftVel); double vdrift=TpcDigiMapper::getInstance()->getGas()->VDrift(); // std::cout<<"Number of events to process:"<init(par); // TpcAlignmentManager::init(par->getAlignmentFile()); TDatabasePDG* pdgDB = TDatabasePDG::Instance(); const Double_t kAu2Gev=0.9314943228; const Double_t khSlash = 1.0545726663e-27; const Double_t kErg2Gev = 1/1.6021773349e-3; const Double_t khShGev = khSlash*kErg2Gev; const Double_t kYear2Sec = 3600*24*365.25; // Ions // if ( !pdgDB->GetParticle(1000010020) ) pdgDB->AddParticle("Deuteron","Deuteron",2*kAu2Gev+8.071e-3,kTRUE, 0,3,"Ion",1000010020); if ( !pdgDB->GetParticle(1000010030) ) pdgDB->AddParticle("Triton","Triton",3*kAu2Gev+14.931e-3,kFALSE, khShGev/(12.33*kYear2Sec),3,"Ion",1000010030); if ( !pdgDB->GetParticle(1000020040) ) pdgDB->AddParticle("Alpha","Alpha",4*kAu2Gev+2.424e-3,kTRUE, khShGev/(12.33*kYear2Sec),6,"Ion",1000020040); if ( !pdgDB->GetParticle(1000020030) ) pdgDB->AddParticle("HE3","HE3",3*kAu2Gev+14.931e-3,kFALSE, 0,6,"Ion",1000020030); std::cout<<"post init"<Run(0,0); // 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; // ------------------------------------------------------------------------ }