void runRecoFOPI_GiBUU(int runNumber, 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("GiBUU_FopiMC_set%i_run%i.mc.root",jobset,runNumber); }else{ inSimFile=Form("GiBUU_FopiMC_%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 = khSlash*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); // ----- 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); PndGeoHandling* geoH = PndGeoHandling::Instance(); Double_t taroff = -65; Double_t targetpos = taroff+41.; // 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); // ------- 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(); // keep Clusters tpcCF->timeslice(6); //in samples tpcCF->SetThreshold(1); tpcCF->SetSingleDigiClusterAmpCut(4.); tpcCF->SetClusterAmpCut(2.); // 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(); //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->SetTpcClusterBranchName("TpcSPHit"); kalman->SetTrackBranchName("TpcTrackPreFit"); kalman->SetOutBranchName("TpcTrackPostFit"); kalman->SetPersistence(); kalman->SetNumIterations(3); // number of fitting iterations (back and forth) fRun->AddTask(kalman); TpcVertexingTask* vtxTask = new TpcVertexingTask(); vtxTask->SetTrackBranchName("TpcTrackPostFit"); vtxTask->SetPersistence(); vtxTask->SetBeamspot(TVector3(0,0,0),2); fRun->AddTask(vtxTask); TpcSimpleVertexingTask* vtxTaskSimple=new TpcSimpleVertexingTask(); vtxTaskSimple->SetTrackBranchName("TpcTrackPostFit"); vtxTaskSimple->SetPersistence(); fRun->AddTask(vtxTaskSimple); FopiLambdaRaveVtxTask* raveTask = new FopiLambdaRaveVtxTask(); raveTask->SetTrackBranchName("TpcTrackPostFit"); raveTask->SetVtxMethod("kalman-smoothing:0"); raveTask->SetUseMCids(); fRun->AddTask(raveTask); FopiLambdaAnaTask3* lTask=new FopiLambdaAnaTask3(); lTask->SetTrackBranchName("TpcTrackPostFit"); lTask->SetPersistence(); lTask->SetMcPid(); fRun->AddTask(lTask); // ----- Intialise and run -------------------------------------------- std::cerr<<"1"<getContainer("TpcDigiPar"); // par->setInputVersion(fRun->GetRunId(),1); // par->setChanged(kTRUE); // par->init(); fRun->Init(); TpcDigiPar* par = FairRun::Instance()->GetRuntimeDb()->getContainer("TpcDigiPar"); 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; // ------------------------------------------------------------------------ }