void runAna_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(inSimFile); cout<<"Set Input File: "<AddFriend(inRecoFile); cout<<"Set Friend MC File: "<SetOutputFile(outFile); cout<<"Set Output File: "<GetRuntimeDb(); FairParRootFileIo* parInput1 = new FairParRootFileIo(); parInput1->open(parFile.Data()); // FairParAsciiFileIo* parIo1 = new FairParAsciiFileIo(); // parIo1->open(allDigiFile.Data(),"in"); rtdb->setFirstInput(parInput1); // rtdb->printParamContexts(); // 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(100); fRun->AddTask(evCount); TpcMCResCalcCluster* mc3dCluster = new TpcMCResCalcCluster(); mc3dCluster->addBranchName("TpcPoint"); mc3dCluster->addBranchName("TpcCluster"); mc3dCluster->skipSecondaries(); //the task: TpcRefTrackResidualTask* resCluster = new TpcRefTrackResidualTask(); resCluster->SetOutBranch("MC_Cluster_Residuals"); resCluster->SetResCalculator(mc3dCluster); resCluster->SetPersistence(); // fRun->AddTask(resCluster); TpcMCResCalcSPHit* mc3dSPHit = new TpcMCResCalcSPHit(); mc3dSPHit->addBranchName("TpcPoint"); mc3dSPHit->addBranchName("TpcSPHit"); mc3dSPHit->addBranchName("TpcCluster"); mc3dSPHit->skipSecondaries(); //the task: TpcRefTrackResidualTask* resSPHit = new TpcRefTrackResidualTask(); resSPHit->SetOutBranch("MC_SPHit_Residuals"); resSPHit->SetResCalculator(mc3dSPHit); resSPHit->SetPersistence(); fRun->AddTask(resSPHit); MCTrackResCalc* mcTrackResCalc = new MCTrackResCalc(); mcTrackResCalc->addBranchName("TpcTrackPostFit"); mcTrackResCalc->addBranchName("MCTrack"); mcTrackResCalc->setExtrapolateToMcStartPos(); // mcTrackResCalc->addBranchName("TpcPoint"); 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(); 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; // ------------------------------------------------------------------------ }