geomAcceptEff/sim/Fun4All_HFG.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
0013 #include <phpythia8/PHPy8ParticleTrigger.h>
0014
0015 #include <decayfinder/DecayFinder.h>
0016 #include <hftrackefficiency/HFTrackEfficiency.h>
0017 #include <kfparticle_sphenix/KFParticle_sPHENIX.h>
0018
0019 #include <ffamodules/FlagHandler.h>
0020 #include <ffamodules/HeadReco.h>
0021 #include <ffamodules/SyncReco.h>
0022 #include <ffamodules/CDBInterface.h>
0023 #include <phool/PHRandomSeed.h>
0024 #include <phool/recoConsts.h>
0025
0026 #include <fun4all/Fun4AllRunNodeInputManager.h>
0027 #include <fun4all/Fun4AllDstOutputManager.h>
0028 #include <fun4all/Fun4AllOutputManager.h>
0029 #include <fun4all/Fun4AllServer.h>
0030
0031 #include <simqa_modules/QAG4SimulationTracking.h>
0032 #include <qautils/QAHistManagerDef.h>
0033
0034 #include "HF_selections.C"
0035
0036 #include <Trkr_Eval.C>
0037
0038
0039 R__LOAD_LIBRARY(libfun4all.so)
0040 R__LOAD_LIBRARY(libffamodules.so)
0041 R__LOAD_LIBRARY(libdecayfinder.so)
0042 R__LOAD_LIBRARY(libhftrackefficiency.so)
0043 R__LOAD_LIBRARY(libsimqa_modules.so)
0044
0045 int Fun4All_HFG(std::string processID = "0", std::string channel = "cascade")
0046 {
0047   int nEvents = 1e3;
0048   //int nEvents = 1e1;
0049   string date = "20260513";
0050   std::string outDir = "./" + channel + "_" + date + "/";
0051
0052   string makeDirectory = "mkdir -p " + outDir + "hfEff";
0053   system(makeDirectory.c_str());
0054
0055   //F4A setup
0056   Fun4AllServer *se = Fun4AllServer::instance();
0057   se->Verbosity(1);
0058
0059   PHRandomSeed::Verbosity(1);
0060   recoConsts *rc = recoConsts::instance();
0061
0062   //Generator setup
0063
0064   Input::PYTHIA8 = true;
0065   int particleID = 421;
0066   PYTHIA8::config_file[0] = "steeringCards/pythia8_MB_Detroit_Tony.cfg";
0067   if (channel == "Kshort2pipi")
0068   {
0069     run_pipi_reco = false;
0070     particleID = 310;
0071   }
0072   else if (channel == "Lambda2ppi")
0073   {
0074     run_ppi_reco = false;
0075     particleID = 3122;
0076   }
0077   else if (channel == "Lambdabar2ppi")
0078   {
0079     run_anti_ppi_reco = false;
0080     particleID = -3122;
0081   }
0082   else if (channel == "cascade")
0083   {
0084     run_cascade_reco = false;
0085     particleID = 3312;
0086   }
0087   else if (channel == "minBias")
0088   {
0089     std::cout << "Min bias simulations" << std::endl;
0090   }
0091   else
0092   {
0093     std::cout << "Your decay channel " << channel << " is not known" << std::endl;
0094     exit(1);
0095   }
0096   Input::BEAM_CONFIGURATION = Input::pp_COLLISION;
0097
0098   InputInit();
0099
0100   float abs_eta = 10;
0101
0102   if (channel != "minBias")
0103   {
0104
0105     int pidtrig = particleID;
0106
0107     auto trigger = new PHPy8ParticleTrigger();
0108     trigger->Verbosity(0);
0109     trigger->AddParticles(pidtrig);
0110     trigger->AddParticles(-1 * pidtrig);
0111     trigger->SetYHighLow(1.2, -1.2);
0112     trigger->SetStableParticleOnly(false);
0113     trigger->PrintConfig();
0114     // assign
0115     INPUTGENERATOR::Pythia8[0]->register_trigger(trigger);
0116     INPUTGENERATOR::Pythia8[0]->set_trigger_OR();
0117     Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8[0]);
0118
0119 /*
0120     PHPy8ParticleTrigger * p8_hf_signal_trigger = new PHPy8ParticleTrigger("thisTrigger");
0121     p8_hf_signal_trigger->SetPtLow(0.);
0122     p8_hf_signal_trigger->SetPtHigh(5.);
0123     p8_hf_signal_trigger->SetYHighLow(1, -1); // sample a rapidity range higher than the sPHENIX tracking pseudorapidity
0124     p8_hf_signal_trigger->SetStableParticleOnly(false); // process unstable particles that include quarks
0125     p8_hf_signal_trigger->SetParticleRadialDecayVertexHigh(4);
0126     p8_hf_signal_trigger->AddParticles(particleID);
0127     p8_hf_signal_trigger->AddParticles(-1*particleID);
0128     p8_hf_signal_trigger->PrintConfig();
0129     INPUTGENERATOR::Pythia8[0]->register_trigger(p8_hf_signal_trigger);
0130     INPUTGENERATOR::Pythia8[0]->set_trigger_OR();
0131 */
0132 /*
0133     vector<int> particleList;
0134     if (channel == "Kshort2pipi") particleList = {211, -211};
0135     else if (channel == "Lambda2ppi") particleList = {2212, -211};
0136     else if (channel == "cascade") particleList = {3122, -211};
0137     else particleList = {-2212, 211};
0138
0139     for (unsigned int i = 0; i < particleList.size(); ++i)
0140     {
0141       string trigger_name = "particle_trigger_" + to_string(i);
0142       PHPy8ParticleTrigger * p8_hf_signal_trigger = new PHPy8ParticleTrigger(trigger_name.c_str());
0143       p8_hf_signal_trigger->SetEtaHighLow(abs_eta, -1*abs_eta);
0144       p8_hf_signal_trigger->SetPtLow(0.15);
0145       p8_hf_signal_trigger->SetParentRadialDecayVertexHigh(4);
0146       p8_hf_signal_trigger->SetStableParticleOnly(false);
0147       p8_hf_signal_trigger->AddParents(particleID);
0148       p8_hf_signal_trigger->AddParticles(particleList[i]);
0149       p8_hf_signal_trigger->PrintConfig();
0150       INPUTGENERATOR::Pythia8[0]->register_trigger(p8_hf_signal_trigger);
0151     }
0152     INPUTGENERATOR::Pythia8[0]->set_trigger_AND();
0153 */
0154   }
0155
0156   Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8[0]);
0157
0158   InputRegister();
0159
0160   //CDB flags and such
0161
0162   Enable::CDB = true;
0163   rc->set_StringFlag("CDB_GLOBALTAG","MDC2");
0164   rc->set_uint64Flag("TIMESTAMP",1);
0165   rc->set_IntFlag("RUNNUMBER",1);
0166
0167   Enable::MVTX_APPLYMISALIGNMENT = true; //Check if this is default
0168   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0169
0170   SyncReco *sync = new SyncReco();
0171   se->registerSubsystem(sync);
0172
0173   HeadReco *head = new HeadReco();
0174   se->registerSubsystem(head);
0175
0176   FlagHandler *flag = new FlagHandler();
0177   se->registerSubsystem(flag);
0178
0179   //Simulation setup
0180   Enable::MBDFAKE = true;
0181   Enable::PIPE = true;
0182   Enable::PIPE_ABSORBER = true;
0183   Enable::MVTX = true;
0184   Enable::INTT = true;
0185   Enable::TPC = true;
0186   Enable::MICROMEGAS = true;
0187
0188   //Tracking setup
0189
0190   G4Init();
0191   MagnetInit();
0192   MagnetFieldInit();
0193
0194   G4Setup();
0195
0196   //Tagging stuff
0197   DecayFinder *myFinder = new DecayFinder("myFinder");
0198   myFinder->Verbosity(INT_MAX);
0199   if (channel == "Kshort2pipi") myFinder->setDecayDescriptor(pipi_decay_descriptor);
0200   else if (channel == "Lambda2ppi") myFinder->setDecayDescriptor(ppi_decay_descriptor);
0201   else if (channel == "cascade") myFinder->setDecayDescriptor(cascade_decay_descriptor);
0202   else myFinder->setDecayDescriptor(anti_ppi_decay_descriptor);
0203   myFinder->saveDST(1);
0204   myFinder->allowPi0(1);
0205   myFinder->allowPhotons(1);
0206   myFinder->triggerOnDecay(0);
0207   myFinder->setPTmin(0.); //Note: sPHENIX min pT is 0.2 GeV for tracking
0208   myFinder->setEtaRange(-1*abs_eta, abs_eta); //Note: sPHENIX acceptance is |eta| <= 1.1
0209   myFinder->useDecaySpecificEtaRange(false);
0210   if (channel != "minBias") se->registerSubsystem(myFinder);
0211
0212   Mbd_Reco();
0213   Mvtx_Cells();
0214   Intt_Cells();
0215   TPC_Cells();
0216   Micromegas_Cells();
0217
0218   TrackingInit();
0219
0220   Mvtx_Clustering();
0221   Intt_Clustering();
0222   TPC_Clustering();
0223   Micromegas_Clustering();
0224
0225   Tracking_Reco();
0226
0227   auto vtxfinder = new PHSimpleVertexFinder;
0228   vtxfinder->Verbosity(0);
0229   vtxfinder->setDcaCut(1);
0230   vtxfinder->setTrackPtCut(-99999.);
0231   vtxfinder->setBeamLineCut(1);
0232   vtxfinder->setTrackQualityCut(1000000000);
0233   vtxfinder->setRequireMVTX(false); //Set back to true?
0234   //vtxfinder->setNmvtxRequired(2);
0235   vtxfinder->setOutlierPairCut(1);
0236   se->registerSubsystem(vtxfinder);
0237
0238   Global_Reco();
0239
0240   build_truthreco_tables();
0241
0242 /*
0243  *
0244  * Thomas, I also added the lines below (ending before HFTrackEfficiency)
0245  *
0246  */
0247   string evalDir = outDir + "evaluatorFiles/" + channel + "/" + date + "/";
0248   makeDirectory = "mkdir -p " + evalDir;
0249   system(makeDirectory.c_str());
0250   string evalFile = evalDir + "evalFile_" + channel + "_" + processID + ".root";
0251
0252   SvtxEvaluator* eval;
0253   eval = new SvtxEvaluator("SVTXEVALUATOR", evalFile, "SvtxTrackMap",
0254                            G4MVTX::n_maps_layer,
0255                            G4INTT::n_intt_layer,
0256                            G4TPC::n_gas_layer,
0257                            G4MICROMEGAS::n_micromegas_layer);
0258   eval->do_cluster_eval(false);
0259   eval->do_g4hit_eval(false);
0260   eval->do_hit_eval(false);  // enable to see the hits that includes the chamber physics...
0261   eval->do_gpoint_eval(false);
0262   eval->do_vtx_eval_light(false);
0263   eval->do_eval_light(false);
0264   eval->do_info_eval(false);
0265   eval->do_vertex_eval(false);
0266   eval->do_g4cluster_eval(false);
0267   eval->do_track_eval(false);
0268   eval->do_gtrack_eval(true);
0269   eval->do_track_match(true);
0270   eval->set_use_initial_vertex(G4TRACKING::g4eval_use_initial_vertex);
0271   bool embed_scan = false;
0272   if (TRACKING::pp_mode) embed_scan = false;
0273   eval->scan_for_embedded(embed_scan);   // take all tracks if false - take only embedded tracks if true
0274   eval->scan_for_primaries(embed_scan);  // defaults to only thrown particles for ntp_gtrack
0275   eval->Verbosity(0);
0276   se->registerSubsystem(eval);
0277
0278   HFTrackEfficiency *myTrackEff = new HFTrackEfficiency("myTrackEff");
0279   myTrackEff->Verbosity(INT_MAX);
0280   myTrackEff->setDFNodeName("myFinder");
0281   myTrackEff->triggerOnDecay(0);
0282   myTrackEff->writeSelectedTrackMap(true);
0283   myTrackEff->writeOutputFile(true);
0284   std::string outputHFEffFile = outDir + "/hfEff/outputHFTrackEff_" + channel + "_" + processID + ".root";
0285   myTrackEff->setOutputFileName(outputHFEffFile);
0286   if (channel != "minBias" && false) se->registerSubsystem(myTrackEff);
0287
0288   output_dir = outDir;
0289
0290   if (run_pipi_reco) create_hf_directories(pipi_reconstruction_name, pipi_output_dir, pipi_output_reco_file, processID);
0291   if (run_ppi_reco) create_hf_directories(ppi_reconstruction_name, ppi_output_dir, ppi_output_reco_file, processID);
0292   if (run_anti_ppi_reco) create_hf_directories(anti_ppi_reconstruction_name, anti_ppi_output_dir, anti_ppi_output_reco_file, processID);
0293   if (run_cascade_reco) create_hf_directories(cascade_reconstruction_name, cascade_output_dir, cascade_output_reco_file, processID);
0294
0295   if (run_pipi_reco) reconstruct_pipi_mass();
0296   if (run_ppi_reco) reconstruct_ppi_mass();
0297   if (run_anti_ppi_reco) reconstruct_ppi_mass();
0298   if (run_cascade_reco) reconstruct_Lambdapi_mass();
0299 /*
0300   //Output file handling
0301   makeDirectory = "mkdir -p " + outDir + "DST";
0302   system(makeDirectory.c_str());
0303
0304   string FullOutFile = outDir + "/DST/" + channel + "_DST_" + processID + ".root";
0305   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
0306   out->StripNode("G4HIT_PIPE");
0307   out->StripNode("G4HIT_SVTXSUPPORT");
0308   //out->StripNode("PHG4INEVENT");
0309   //out->StripNode("Sync");
0310   out->StripNode("myFinder_DecayMap");
0311   out->StripNode("G4HIT_PIPE");
0312   out->StripNode("G4HIT_MVTX");
0313   out->StripNode("G4HIT_INTT");
0314   out->StripNode("G4HIT_TPC");
0315   out->StripNode("G4HIT_MICROMEGAS");
0316   out->StripNode("TRKR_HITSET");
0317   out->StripNode("TRKR_HITTRUTHASSOC");
0318   //out->StripNode("TRKR_CLUSTER");
0319   //out->StripNode("TRKR_CLUSTERHITASSOC");
0320   out->StripNode("TRKR_CLUSTERCROSSINGASSOC");
0321   out->StripNode("TRAINING_HITSET");
0322   out->StripNode("TRKR_TRUTHTRACKCONTAINER");
0323   out->StripNode("TRKR_TRUTHCLUSTERCONTAINER");
0324   out->StripNode("alignmentTransformationContainer");
0325   out->StripNode("alignmentTransformationContainerTransient");
0326   //out->StripNode("SiliconTrackSeedContainer");
0327   //out->StripNode("TpcTrackSeedContainer");
0328   //out->StripNode("SvtxTrackSeedContainer");
0329   out->StripNode("ActsTrajectories");
0330   //out->StripNode("SvtxTrackMap");
0331   out->StripNode("SvtxAlignmentStateMap");
0332   //out->SaveRunNode(0);
0333   se->registerOutputManager(out);
0334 */
0335
0336   se->run(nEvents);
0337
0338   se->End();
0339
0340   if (run_pipi_reco) end_kfparticle(pipi_output_reco_file, pipi_output_dir);
0341   if (run_ppi_reco) end_kfparticle(ppi_output_reco_file, ppi_output_dir);
0342   if (run_anti_ppi_reco) end_kfparticle(ppi_output_reco_file, ppi_output_dir);
0343   if (run_cascade_reco) end_kfparticle(cascade_output_reco_file, cascade_output_dir);
0344
0345   gSystem->Exit(0);
0346
0347   return 0;
0348 }
0349
0350 #endif