geometric_acceptance/simulation/Fun4All_HFG_MB.C

0001 #ifndef MACRO_FUN4ALLG4SPHENIX_C
0002 #define MACRO_FUN4ALLG4SPHENIX_C
0003
0004 #include "G4_Input.C"
0005 #include <G4_Global.C>
0006 #include <G4Setup_sPHENIX.C>
0007
0008 #include <Trkr_RecoInit.C>
0009 #include <Trkr_Clustering.C>
0010 #include <Trkr_TruthTables.C>
0011 #include <Trkr_Reco.C>
0012 #include <Trkr_Eval.C>
0013
0014 #include <phpythia8/PHPy8ParticleTrigger.h>
0015
0016 #include <decayfinder/DecayFinder.h>
0017 #include <hftrackefficiency/HFTrackEfficiency.h>
0018 #include <kfparticle_sphenix/KFParticle_sPHENIX.h>
0019
0020 #include <ffamodules/FlagHandler.h>
0021 #include <ffamodules/HeadReco.h>
0022 #include <ffamodules/SyncReco.h>
0023 #include <ffamodules/CDBInterface.h>
0024 #include <phool/PHRandomSeed.h>
0025 #include <phool/recoConsts.h>
0026
0027 #include <fun4all/Fun4AllRunNodeInputManager.h>
0028 #include <fun4all/Fun4AllDstOutputManager.h>
0029 #include <fun4all/Fun4AllOutputManager.h>
0030 #include <fun4all/Fun4AllServer.h>
0031
0032 #include <simqa_modules/QAG4SimulationTracking.h>
0033 #include <qautils/QAHistManagerDef.h>
0034
0035 #include "HF_selections.C"
0036
0037 R__LOAD_LIBRARY(libfun4all.so)
0038 R__LOAD_LIBRARY(libffamodules.so)
0039 R__LOAD_LIBRARY(libdecayfinder.so)
0040 R__LOAD_LIBRARY(libhftrackefficiency.so)
0041 R__LOAD_LIBRARY(libsimqa_modules.so)
0042
0043 int Fun4All_HFG_MB(std::string processID = "000000", std::string channel = "lambdaKshort")
0044 {
0045   int nEvents = 1000;
0046
0047   int ndigits = 6;
0048   //std::string infile_base_reco = "/sphenix/tg/tg01/hf/mjpeters/lambdaKshortMB/G4Hits_pileup_sample/DST_TRKR_G4HIT_pythia8_NONE-0000000029-";
0049   //std::string infile_base_truth = "/sphenix/tg/tg01/hf/mjpeters/lambdaKshortMB/G4Hits_pileup_sample/DST_TRUTH_G4HIT_pythia8_NONE-0000000029-";
0050   std::string infile_base_reco = "DST_TRKR_G4HIT_pythia8_Detroit_200kHz-0000000029-";
0051   std::string infile_base_truth = "DST_TRUTH_G4HIT_pythia8_Detroit_200kHz-0000000029-";
0052   std::string infile_base = "/sphenix/lustre01/sphnxpro/mdc2/js_pp200_signal/g4hits/run0029/detroit/G4Hits_pythia8_Detroit-0000000029-";
0053   //std::string process_str = std::string(ndigits - std::to_string(process).length(),'0') + std::to_string(process);
0054   std::string infile = infile_base+processID+".root";
0055   std::string infile_reco = infile_base_reco+processID+".root";
0056   std::string infile_truth = infile_base_truth+processID+".root";
0057
0058   //std::string outDir = "/sphenix/tg/tg01/hf/mjpeters/lambdaKshortMB/" + channel + "_20260422_DetroitMB_CR_2_mode_pTref_1p4/";
0059   //std::string outDir = "./";
0060
0061   std::string outDir = "/sphenix/tg/tg01/hf/mjpeters/LightFlavorProduction/closureTestSample/";
0062
0063   string makeDirectory = "mkdir -p " + outDir + "hfEff";
0064   system(makeDirectory.c_str());
0065   makeDirectory = "mkdir -p " + outDir + "evaluator";
0066   system(makeDirectory.c_str());
0067
0068   //F4A setup
0069   Fun4AllServer *se = Fun4AllServer::instance();
0070   se->Verbosity(1);
0071
0072   PHRandomSeed::Verbosity(1);
0073   recoConsts *rc = recoConsts::instance();
0074
0075   Enable::VERBOSITY = 0;
0076   Input::VERBOSITY = 0;
0077
0078   //Generator setup
0079
0080   //Input::PILEUPRATE = 3e6;
0081
0082
0083   Input::READHITS = false;
0084   //INPUTREADHITS::filename[0] = infile;
0085   INPUTREADHITS::filename[0] = infile_reco;
0086   INPUTREADHITS::filename[1] = infile_truth;
0087
0088   std::cout << "readhits ssetup" << std::endl;
0089   Input::SIMPLE = true;
0090
0091   Input::PYTHIA8 = true;
0092   int particleID = 421;
0093   PYTHIA8::config_file[0] = "steeringCards/pythia8_MB_Detroit_Tony.cfg";
0094   if (channel == "Kshort2pipi")
0095   {
0096     run_pipi_reco = true;
0097     particleID = 310;
0098   }
0099   else if (channel == "Lambda2ppi")
0100   {
0101     run_ppi_reco = true;
0102     particleID = 3122;
0103   }
0104   else if (channel == "Lambdabar2ppi")
0105   {
0106     run_anti_ppi_reco = true;
0107     particleID = -3122;
0108   }
0109   else if (channel == "cascade")
0110   {
0111     run_cascade_reco = true;
0112     particleID = 3312;
0113   }
0114   else if (channel == "minBias")
0115   {
0116     std::cout << "Min bias simulations" << std::endl;
0117   }
0118   else if (channel == "lambdaKshort")
0119   {
0120     run_pipi_reco = true;
0121     run_ppi_reco = true;
0122     run_anti_ppi_reco = true;
0123   }
0124   else
0125   {
0126     std::cout << "Your decay channel " << channel << " is not known" << std::endl;
0127     exit(1);
0128   }
0129   Input::BEAM_CONFIGURATION = Input::pp_COLLISION;
0130
0131   InputInit();
0132
0133   std::cout << "input init" << std::endl;
0134
0135   float abs_eta = 1.1;
0136
0137
0138 /*
0139   if (channel != "minBias" && channel != "lambdaKshort")
0140   {
0141
0142     int pidtrig = particleID;
0143
0144     auto trigger = new PHPy8ParticleTrigger();
0145     trigger->Verbosity(0);
0146     trigger->AddParticles(pidtrig);
0147     trigger->AddParticles(-1 * pidtrig);
0148     trigger->SetYHighLow(1.2, -1.2);
0149     trigger->SetStableParticleOnly(false);
0150     trigger->PrintConfig();
0151     // assign
0152     INPUTGENERATOR::Pythia8[0]->register_trigger(trigger);
0153     INPUTGENERATOR::Pythia8[0]->set_trigger_OR();
0154     Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8[0]);
0155 */
0156 /*
0157     PHPy8ParticleTrigger * p8_hf_signal_trigger = new PHPy8ParticleTrigger("thisTrigger");
0158     p8_hf_signal_trigger->SetPtLow(0.);
0159     p8_hf_signal_trigger->SetPtHigh(5.);
0160     p8_hf_signal_trigger->SetYHighLow(1, -1); // sample a rapidity range higher than the sPHENIX tracking pseudorapidity
0161     p8_hf_signal_trigger->SetStableParticleOnly(false); // process unstable particles that include quarks
0162     p8_hf_signal_trigger->SetParticleRadialDecayVertexHigh(4);
0163     p8_hf_signal_trigger->AddParticles(particleID);
0164     p8_hf_signal_trigger->AddParticles(-1*particleID);
0165     p8_hf_signal_trigger->PrintConfig();
0166     INPUTGENERATOR::Pythia8[0]->register_trigger(p8_hf_signal_trigger);
0167     INPUTGENERATOR::Pythia8[0]->set_trigger_OR();
0168
0169     vector<int> particleList;
0170     if (channel == "Kshort2pipi") particleList = {211, -211};
0171     else if (channel == "Lambda2ppi") particleList = {2212, -211};
0172     else if (channel == "cascade") particleList = {3122, -211};
0173     else particleList = {-2212, 211};
0174
0175     for (unsigned int i = 0; i < particleList.size(); ++i)
0176     {
0177       string trigger_name = "particle_trigger_" + to_string(i);
0178       PHPy8ParticleTrigger * p8_hf_signal_trigger = new PHPy8ParticleTrigger(trigger_name.c_str());
0179       p8_hf_signal_trigger->SetEtaHighLow(abs_eta, -1*abs_eta);
0180       p8_hf_signal_trigger->SetPtLow(0.15);
0181       //p8_hf_signal_trigger->SetParentRadialDecayVertexHigh(4);
0182       p8_hf_signal_trigger->SetStableParticleOnly(false);
0183       p8_hf_signal_trigger->AddParents(particleID);
0184       p8_hf_signal_trigger->AddParticles(particleList[i]);
0185       p8_hf_signal_trigger->PrintConfig();
0186       INPUTGENERATOR::Pythia8[0]->register_trigger(p8_hf_signal_trigger);
0187     }
0188     INPUTGENERATOR::Pythia8[0]->set_trigger_AND();
0189   }
0190 */
0191   if(channel == "lambdaKshort")
0192   {
0193     run_pipi_reco = true;
0194     run_ppi_reco = true;
0195     run_anti_ppi_reco = true;
0196 /*
0197     PHPy8ParticleTrigger* p8_lktrigger = new PHPy8ParticleTrigger("lktrigger");
0198     p8_lktrigger->SetEtaHighLow(abs_eta,-1*abs_eta);
0199     p8_lktrigger->SetPtLow(0.);
0200     //p8_lktrigger->SetParentRadialDecayVertexHigh(4);
0201     p8_lktrigger->SetStableParticleOnly(false);
0202
0203     std::vector<int> mothers = {310,3122,-3122};
0204
0205     p8_lktrigger->AddParticles(mothers);
0206     p8_lktrigger->PrintConfig();
0207     INPUTGENERATOR::Pythia8[0]->register_trigger(p8_lktrigger);
0208 */
0209   }
0210
0211   Input::BEAM_CONFIGURATION = Input::pp_COLLISION;
0212   Enable::MVTX_APPLYMISALIGNMENT = true;
0213
0214   Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8[0]);
0215
0216   if(Input::SIMPLE)
0217   {
0218     INPUTGENERATOR::SimpleEventGenerator[0]->add_particles(310,10);
0219     INPUTGENERATOR::SimpleEventGenerator[0]->add_particles(3122,5);
0220     INPUTGENERATOR::SimpleEventGenerator[0]->add_particles(-3122,5);
0221     INPUTGENERATOR::SimpleEventGenerator[0]->set_reuse_existing_vertex(true);
0222     INPUTGENERATOR::SimpleEventGenerator[0]->set_eta_range(-1.,1.);
0223     INPUTGENERATOR::SimpleEventGenerator[0]->set_phi_range(-M_PI,M_PI);
0224     INPUTGENERATOR::SimpleEventGenerator[0]->set_pt_range(0.,5.);
0225   }
0226
0227   InputRegister();
0228
0229   std::cout << "input registered" << std::endl;
0230
0231   //CDB flags and such
0232
0233   Enable::CDB = true;
0234   //rc->set_StringFlag("CDB_GLOBALTAG","MDC2");
0235   rc->set_StringFlag("CDB_GLOBALTAG",CDB::global_tag);
0236   //rc->set_StringFlag("CDB_GLOBALTAG","ProdA_2024");
0237   rc->set_uint64Flag("TIMESTAMP",CDB::timestamp);
0238   rc->set_IntFlag("RUNNUMBER",29);
0239
0240   Enable::MVTX_APPLYMISALIGNMENT = true;
0241   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0242 /*
0243   SyncReco *sync = new SyncReco();
0244   se->registerSubsystem(sync);
0245
0246   HeadReco *head = new HeadReco();
0247   se->registerSubsystem(head);
0248 */
0249   FlagHandler *flag = new FlagHandler();
0250   se->registerSubsystem(flag);
0251
0252   //Simulation setup
0253   Enable::MBDFAKE = true;
0254   Enable::PIPE = true;
0255   Enable::PIPE_ABSORBER = true;
0256   Enable::MVTX = true;
0257   Enable::INTT = true;
0258   Enable::TPC = true;
0259   Enable::MICROMEGAS = true;
0260
0261   //Tracking setup
0262   TRACKING::streaming_mode = true;
0263
0264   InputManagers();
0265
0266   G4Init();
0267   MagnetInit();
0268   MagnetFieldInit();
0269
0270   G4Setup();
0271
0272   std::cout << "g4setup" << std::endl;
0273 /*
0274   //Tagging stuff
0275   DecayFinder *myFinder = new DecayFinder("myFinder");
0276   myFinder->Verbosity(INT_MAX);
0277   if (channel == "Kshort2pipi") myFinder->setDecayDescriptor(pipi_decay_descriptor);
0278   else if (channel == "Lambda2ppi") myFinder->setDecayDescriptor(ppi_decay_descriptor);
0279   else if (channel == "cascade") myFinder->setDecayDescriptor(cascade_decay_descriptor);
0280   else myFinder->setDecayDescriptor(anti_ppi_decay_descriptor);
0281   myFinder->saveDST(1);
0282   myFinder->allowPi0(1);
0283   myFinder->allowPhotons(1);
0284   myFinder->triggerOnDecay(1);
0285   myFinder->setPTmin(0.); //Note: sPHENIX min pT is 0.2 GeV for tracking
0286   myFinder->setEtaRange(-1*abs_eta, abs_eta); //Note: sPHENIX acceptance is |eta| <= 1.1
0287   myFinder->useDecaySpecificEtaRange(false);
0288   if (channel != "minBias") se->registerSubsystem(myFinder);
0289 */
0290   Mbd_Reco();
0291   Mvtx_Cells();
0292   Intt_Cells();
0293   TPC_Cells();
0294   Micromegas_Cells();
0295
0296   TrackingInit();
0297
0298   Mvtx_Clustering();
0299   Intt_Clustering();
0300   TPC_Clustering();
0301   Micromegas_Clustering();
0302
0303   Tracking_Reco();
0304
0305   auto vtxfinder = new PHSimpleVertexFinder;
0306   vtxfinder->Verbosity(0);
0307   vtxfinder->setDcaCut(1);
0308   vtxfinder->setTrackPtCut(-99999.);
0309   vtxfinder->setBeamLineCut(1);
0310   vtxfinder->setTrackQualityCut(1000000000);
0311   vtxfinder->setRequireMVTX(false);
0312   //vtxfinder->setNmvtxRequired(2);
0313   vtxfinder->setOutlierPairCut(1);
0314   se->registerSubsystem(vtxfinder);
0315
0316   Global_Reco();
0317
0318   build_truthreco_tables();
0319
0320   Tracking_Eval(outDir+"evaluator/outputeval_"+processID+".root");
0321 /*
0322   HFTrackEfficiency *myTrackEff = new HFTrackEfficiency("myTrackEff");
0323   myTrackEff->Verbosity(INT_MAX);
0324   myTrackEff->setDFNodeName("myFinder");
0325   myTrackEff->triggerOnDecay(1);
0326   myTrackEff->writeSelectedTrackMap(true);
0327   myTrackEff->writeOutputFile(true);
0328   std::string outputHFEffFile = outDir + "/hfEff/outputHFTrackEff_" + channel + "_" + processID + ".root";
0329   myTrackEff->setOutputFileName(outputHFEffFile);
0330   if (channel != "minBias") se->registerSubsystem(myTrackEff);
0331 */
0332   output_dir = outDir;
0333
0334   if (run_pipi_reco) create_hf_directories(pipi_reconstruction_name, pipi_output_dir, pipi_output_reco_file, processID);
0335   if (run_ppi_reco || run_anti_ppi_reco) create_hf_directories(ppi_reconstruction_name, ppi_output_dir, ppi_output_reco_file, processID);
0336   //if (run_anti_ppi_reco) create_hf_directories(anti_ppi_reconstruction_name, anti_ppi_output_dir, anti_ppi_output_reco_file, processID);
0337   if (run_cascade_reco) create_hf_directories(cascade_reconstruction_name, cascade_output_dir, cascade_output_reco_file, processID);
0338
0339   if (run_pipi_reco) reconstruct_pipi_mass();
0340   if (run_ppi_reco || run_anti_ppi_reco) reconstruct_ppi_mass();
0341   //if (run_anti_ppi_reco) reconstruct_ppi_mass();
0342   if (run_cascade_reco) reconstruct_Lambdapi_mass();
0343 /*
0344   //Output file handling
0345   makeDirectory = "mkdir -p " + outDir + "DST";
0346   system(makeDirectory.c_str());
0347
0348   string FullOutFile = outDir + "/DST/" + channel + "_DST_" + processID + ".root";
0349   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
0350   out->StripNode("G4HIT_PIPE");
0351   out->StripNode("G4HIT_SVTXSUPPORT");
0352   //out->StripNode("PHG4INEVENT");
0353   //out->StripNode("Sync");
0354   out->StripNode("myFinder_DecayMap");
0355   out->StripNode("G4HIT_PIPE");
0356   out->StripNode("G4HIT_MVTX");
0357   out->StripNode("G4HIT_INTT");
0358   out->StripNode("G4HIT_TPC");
0359   out->StripNode("G4HIT_MICROMEGAS");
0360   out->StripNode("TRKR_HITSET");
0361   out->StripNode("TRKR_HITTRUTHASSOC");
0362   //out->StripNode("TRKR_CLUSTER");
0363   //out->StripNode("TRKR_CLUSTERHITASSOC");
0364   out->StripNode("TRKR_CLUSTERCROSSINGASSOC");
0365   out->StripNode("TRAINING_HITSET");
0366   out->StripNode("TRKR_TRUTHTRACKCONTAINER");
0367   out->StripNode("TRKR_TRUTHCLUSTERCONTAINER");
0368   out->StripNode("alignmentTransformationContainer");
0369   out->StripNode("alignmentTransformationContainerTransient");
0370   //out->StripNode("SiliconTrackSeedContainer");
0371   //out->StripNode("TpcTrackSeedContainer");
0372   //out->StripNode("SvtxTrackSeedContainer");
0373   out->StripNode("ActsTrajectories");
0374   //out->StripNode("SvtxTrackMap");
0375   out->StripNode("SvtxAlignmentStateMap");
0376   //out->SaveRunNode(0);
0377   se->registerOutputManager(out);
0378 */
0379   se->run(nEvents);
0380
0381   se->End();
0382
0383   if (run_pipi_reco) end_kfparticle(pipi_output_reco_file, pipi_output_dir);
0384   if (run_ppi_reco || run_anti_ppi_reco) end_kfparticle(ppi_output_reco_file, ppi_output_dir);
0385   //if (run_anti_ppi_reco) end_kfparticle(ppi_output_reco_file, ppi_output_dir);
0386   if (run_cascade_reco) end_kfparticle(cascade_output_reco_file, cascade_output_dir);
0387
0388   gSystem->Exit(0);
0389
0390   return 0;
0391 }
0392
0393 #endif