macros/simulation/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 #include <Trkr_Eval.C>
0013
0014 #include <phpythia8/PHPy8ParticleTrigger.h>
0015
0016 #include <decayfinder/DecayFinder.h>
0017 #include <hftrackefficiency/HFTrackEfficiency.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 <kfparticle_sphenix/KFParticle_sPHENIX.h>
0035 #include <trackingdiagnostics/TrkrNtuplizer.h>
0036
0037 #include "../data/HF_selections_QM25.C"
0038
0039 //#include <dstarreco/DstarReco.h>
0040
0041 //R__LOAD_LIBRARY(libDstarReco.so)
0042 R__LOAD_LIBRARY(libfun4all.so)
0043 R__LOAD_LIBRARY(libffamodules.so)
0044 R__LOAD_LIBRARY(libdecayfinder.so)
0045 R__LOAD_LIBRARY(libhftrackefficiency.so)
0046 R__LOAD_LIBRARY(libsimqa_modules.so)
0047
0048 using namespace HeavyFlavorReco;
0049
0050 int Fun4All_HFG(std::string processID = "0",
0051                 int nEvents = 10,
0052                 std::string channel = "Kshort2pipi",
0053                 std::string outDir = "./")
0054 {
0055   //F4A setup
0056
0057   Fun4AllServer *se = Fun4AllServer::instance();
0058   se->Verbosity(1);
0059
0060   PHRandomSeed::Verbosity(1);
0061   recoConsts *rc = recoConsts::instance();
0062
0063   //Generator setup
0064
0065   Input::PYTHIA8 = true;
0066   int particleID = 421;
0067
0068   if (channel == "bs2jpsiks0")
0069   {
0070     PYTHIA8::config_file = "steeringCards/pythia8_bs2jpsiks0.cfg";
0071     EVTGENDECAYER::DecayFile = "decFiles/Bs2JpsiKS0.DEC";
0072     particleID = 531;
0073   }
0074   else if (channel == "b2DX")
0075   {
0076     PYTHIA8::config_file = "steeringCards/pythia8_b2DX.cfg";
0077     EVTGENDECAYER::DecayFile = "decFiles/b2DX.DEC";
0078     particleID = 5;
0079   }
0080   else if (channel == "D2Kpi")
0081   {
0082     PYTHIA8::config_file = "steeringCards/pythia8_D2Kpi.cfg";
0083     EVTGENDECAYER::DecayFile = "decFiles/D2Kpi.DEC";
0084   }
0085   else if ((channel == "Dp2KKpi") || (channel == "Ds2KKpi"))
0086   {
0087     PYTHIA8::config_file = "steeringCards/pythia8_D2KKpi.cfg";
0088     EVTGENDECAYER::DecayFile = "decFiles/D2KKpi.DEC";
0089     if (channel == "Dp2KKpi") particleID = 411;
0090     else particleID = 431;
0091   }
0092   else if (channel == "Dstar2D0pi")
0093   {
0094     PYTHIA8::config_file = "steeringCards/pythia8_Dstar2D0pi.cfg";
0095     EVTGENDECAYER::DecayFile = "decFiles/Dstar2D0pi.DEC";
0096     particleID = 413;
0097   }
0098   else if (channel == "Kshort2pipi")
0099   {
0100     PYTHIA8::config_file = "steeringCards/pythia8_K2pipi_Detroit.cfg";
0101     EVTGENDECAYER::DecayFile = "decFiles/K2pipi.DEC";
0102     particleID = 310;
0103   }
0104   else if (channel == "phi2KK")
0105   {
0106     PYTHIA8::config_file = "steeringCards/pythia8_phi2KK_Detroit.cfg";
0107     EVTGENDECAYER::DecayFile = "decFiles/phi2KK.DEC";
0108     particleID = 333;
0109   }
0110   else if (channel == "lambda2ppi")
0111   {
0112     PYTHIA8::config_file = "steeringCards/pythia8_lambda2ppi_Detroit.cfg";
0113     EVTGENDECAYER::DecayFile = "decFiles/lambda2ppi.DEC";
0114     particleID = 3122;
0115   }
0116   else if (channel == "minBias")
0117   {
0118     std::cout << "Min bias simulations" << std::endl;
0119   }
0120   else
0121   {
0122     std::cout << "Your decay channel " << channel << " is not known" << std::endl;
0123     exit(1);
0124   }
0125
0126   Input::BEAM_CONFIGURATION = Input::pp_COLLISION;
0127
0128   InputInit();
0129
0130   PHPy8ParticleTrigger * p8_hf_signal_trigger = new PHPy8ParticleTrigger();
0131   p8_hf_signal_trigger->SetPtLow(0.);
0132   p8_hf_signal_trigger->SetEtaHighLow(1.3, -1.3); // sample a rapidity range higher than the sPHENIX tracking pseudorapidity
0133   p8_hf_signal_trigger->SetStableParticleOnly(false); // process unstable particles that include quarks
0134   p8_hf_signal_trigger->PrintConfig();
0135
0136   if (channel != "minBias")
0137   {
0138     p8_hf_signal_trigger->AddParticles(particleID);
0139     p8_hf_signal_trigger->AddParticles(-1*particleID);
0140     INPUTGENERATOR::Pythia8->register_trigger(p8_hf_signal_trigger);
0141     INPUTGENERATOR::Pythia8->set_trigger_OR();
0142
0143     Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8);
0144   }
0145
0146   InputRegister();
0147
0148   //CDB flags and such
0149
0150   Enable::CDB = true;
0151   rc->set_StringFlag("CDB_GLOBALTAG","ProdA_2024");
0152   rc->set_uint64Flag("TIMESTAMP",1);
0153   rc->set_IntFlag("RUNNUMBER",1);
0154
0155   Enable::MVTX_APPLYMISALIGNMENT = true;
0156   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0157
0158   SyncReco *sync = new SyncReco();
0159   se->registerSubsystem(sync);
0160
0161   HeadReco *head = new HeadReco();
0162   se->registerSubsystem(head);
0163
0164   FlagHandler *flag = new FlagHandler();
0165   se->registerSubsystem(flag);
0166
0167   DecayFinder *myFinder = new DecayFinder("myFinder");
0168   myFinder->Verbosity(INT_MAX);
0169   if (channel == "bs2jpsiks0") myFinder->setDecayDescriptor("[B_s0 -> {J/psi -> e^+ e^-} {K_S0 -> pi^+ pi^-}]cc");
0170   else if  (channel == "Kshort2pipi") myFinder->setDecayDescriptor("K_S0 -> pi^- pi^+");
0171   else if  (channel == "Dp2KKpi") myFinder->setDecayDescriptor("[D+ -> {phi -> K^+ K^-} pi^+]cc");
0172   else if  (channel == "Ds2KKpi") myFinder->setDecayDescriptor("[D_s+ -> {phi -> K^+ K^-} pi^+]cc");
0173   else if  (channel == "Dstar2D0pi") myFinder->setDecayDescriptor("[D*+ -> {D0 -> K^- pi^+} pi^+]cc");
0174   else if  (channel == "phi2KK") myFinder->setDecayDescriptor("phi -> K^+ K^-");
0175   else if  (channel == "lambda2ppi") myFinder->setDecayDescriptor("[Lambda0 -> proton^+ pi^-]cc");
0176   else myFinder->setDecayDescriptor("[D0 -> K^- pi^+]cc");
0177   myFinder->saveDST(1);
0178   myFinder->allowPi0(1);
0179   myFinder->allowPhotons(1);
0180   myFinder->triggerOnDecay(1);
0181   myFinder->setPTmin(0.16); //Note: sPHENIX min pT is 0.2 GeV for tracking
0182   myFinder->setEtaRange(-1.2, 1.2); //Note: sPHENIX acceptance is |eta| <= 1.1
0183   myFinder->useDecaySpecificEtaRange(false);
0184   if (channel != "minBias") se->registerSubsystem(myFinder);
0185
0186   //Simulation setup
0187 /*
0188   std::string geofile = CDBInterface::instance()->getUrl("Tracking_Geometry");
0189
0190   Fun4AllRunNodeInputManager *ingeo = new Fun4AllRunNodeInputManager("GeoIn");
0191   ingeo->AddFile(geofile);
0192   se->registerInputManager(ingeo);
0193 */
0194   Enable::MBDFAKE = true;
0195   Enable::PIPE = true;
0196   Enable::PIPE_ABSORBER = true;
0197   Enable::MVTX = true;
0198   Enable::INTT = true;
0199   Enable::TPC = true;
0200   Enable::MICROMEGAS = true;
0201
0202   //Tracking setup
0203
0204   G4Init();
0205 /*
0206   PHG4Reco *g4Reco = new PHG4Reco();
0207   g4Reco->set_rapidity_coverage(1.8);  // according to drawings
0208   WorldInit(g4Reco);
0209 */
0210   MagnetInit();
0211   MagnetFieldInit();
0212
0213   G4Setup();
0214 /*
0215   PHG4TruthSubsystem *truth = new PHG4TruthSubsystem();
0216   g4Reco->registerSubsystem(truth);
0217
0218   se->registerSubsystem(g4Reco);
0219 */
0220   Mbd_Reco();
0221   Mvtx_Cells();
0222   Intt_Cells();
0223   TPC_Cells();
0224   Micromegas_Cells();
0225
0226   TrackingInit();
0227
0228   Mvtx_Clustering();
0229   Intt_Clustering();
0230   TPC_Clustering();
0231   Micromegas_Clustering();
0232
0233   Tracking_Reco();
0234
0235   string eval_dir = outDir + "/trutheval/";
0236   string eval_output_reco_file = eval_dir + processing_folder + "trutheval_output" + trailer;
0237
0238   std::string makeDirectory = "mkdir -p " + eval_dir + processing_folder;
0239   system(makeDirectory.c_str());
0240
0241   Tracking_Eval(eval_output_reco_file);
0242
0243   auto vtxfinder = new PHSimpleVertexFinder;
0244   vtxfinder->Verbosity(0);
0245   vtxfinder->setDcaCut(0.5);
0246   vtxfinder->setTrackPtCut(-99999.);
0247   vtxfinder->setBeamLineCut(1);
0248   vtxfinder->setTrackQualityCut(1000000000);
0249   vtxfinder->setNmvtxRequired(3);
0250   vtxfinder->setOutlierPairCut(0.1);
0251   se->registerSubsystem(vtxfinder);
0252
0253   Global_Reco();
0254
0255   build_truthreco_tables();
0256
0257   //HF stuff
0258
0259   output_dir = outDir; //Top dir of where the output nTuples will be written
0260   trailer = "_" + processID + ".root";
0261
0262   if (run_pipi_reco) create_hf_directories(pipi_reconstruction_name, pipi_output_dir, pipi_output_reco_file);
0263   if (run_ppi_reco) create_hf_directories(ppi_reconstruction_name, ppi_output_dir, ppi_output_reco_file);
0264   if (run_KK_reco) create_hf_directories(KK_reconstruction_name, KK_output_dir, KK_output_reco_file);
0265
0266   if (run_pipi_reco) reconstruct_pipi_mass();
0267   if (run_ppi_reco) reconstruct_ppi_mass();
0268   if (run_KK_reco) reconstruct_KK_mass();
0269
0270   string nTuplizer_dir = output_dir + "/nTP/";
0271   string NTP_output_reco_file = nTuplizer_dir + processing_folder + "nTuplizer_tracks" + trailer;
0272
0273   std::string makeDirectory = "mkdir -p " + nTuplizer_dir + processing_folder;
0274   system(makeDirectory.c_str());
0275
0276   TrkrNtuplizer* eval = new TrkrNtuplizer("TrkrNtuplizer", NTP_output_reco_file, "SvtxTrackMap", 3, 4, 48);
0277   int do_default = 0;
0278   eval->Verbosity(0);
0279   eval->do_hit_eval(false);
0280   eval->do_cluster_eval(false);
0281   eval->do_clus_trk_eval(false);
0282   eval->do_vertex_eval(false);
0283   eval->do_track_eval(true);
0284   eval->do_tpcseed_eval(false);
0285   eval->do_siseed_eval(false);
0286   se->registerSubsystem(eval);
0287
0288   HFTrackEfficiency *myTrackEff = new HFTrackEfficiency("myTrackEff");
0289   myTrackEff->Verbosity(INT_MAX);
0290   //myTrackEff->setTruthRecoMatchingPercentage(100.);
0291   myTrackEff->setDFNodeName("myFinder");
0292   myTrackEff->triggerOnDecay(1);
0293   myTrackEff->writeSelectedTrackMap(true);
0294   myTrackEff->writeOutputFile(true);
0295   std::string outputHFEffFile = outDir + "/hfEff/outputHFTrackEff_" + channel + "_" + processID + ".root";
0296   myTrackEff->setOutputFileName(outputHFEffFile);
0297   if (channel != "minBias") se->registerSubsystem(myTrackEff);
0298 /*
0299   if (channel == "Dstar2D0pi")
0300   {
0301     DstarReco *myDstarReco = new DstarReco();
0302     std::string outputDstarfFile = outDir + "/standalone/outputStandaloneReco_" + channel + "_" + processID + ".root";
0303     myDstarReco->setOutputName(outputDstarfFile);
0304     se->registerSubsystem(myDstarReco);
0305   }
0306
0307   if (run_pipi_reco || run_Kpi_reco || run_Kpipi_reco || run_KKpi_reco || run_pKpi_reco || run_Lambdapi_reco) init_kfp_dependencies();
0308   if (run_pipi_reco) reconstruct_pipi_mass();
0309 */
0310   //Output file handling
0311
0312   string FullOutFile = outDir + "/" + channel + "_DST_" + processID + ".root";
0313   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
0314   //out->StripNode("G4HIT_PIPE");
0315   //out->StripNode("G4HIT_SVTXSUPPORT");
0316   //out->StripNode("PHG4INEVENT");
0317   //out->StripNode("Sync");
0318   out->StripNode("myFinder_DecayMap");
0319   //out->StripNode("G4HIT_PIPE");
0320   //out->StripNode("G4HIT_MVTX");
0321   //out->StripNode("G4HIT_INTT");
0322   //out->StripNode("G4HIT_TPC");
0323   //out->StripNode("G4HIT_MICROMEGAS");
0324   //out->StripNode("TRKR_HITSET");
0325   //out->StripNode("TRKR_HITTRUTHASSOC");
0326   //out->StripNode("TRKR_CLUSTER");
0327   //out->StripNode("TRKR_CLUSTERHITASSOC");
0328   out->StripNode("TRKR_CLUSTERCROSSINGASSOC");
0329   out->StripNode("TRAINING_HITSET");
0330   out->StripNode("TRKR_TRUTHTRACKCONTAINER");
0331   out->StripNode("TRKR_TRUTHCLUSTERCONTAINER");
0332   out->StripNode("alignmentTransformationContainer");
0333   out->StripNode("alignmentTransformationContainerTransient");
0334   //out->StripNode("SiliconTrackSeedContainer");
0335   //out->StripNode("TpcTrackSeedContainer");
0336   //out->StripNode("SvtxTrackSeedContainer");
0337   out->StripNode("ActsTrajectories");
0338   //out->StripNode("SvtxTrackMap");
0339   out->StripNode("SvtxAlignmentStateMap");
0340   //out->SaveRunNode(0);
0341   se->registerOutputManager(out);
0342
0343   se->run(nEvents);
0344
0345   se->End();
0346
0347   end_kfparticle(NTP_output_reco_file, nTuplizer_dir);
0348
0349   if (run_pipi_reco) end_kfparticle(pipi_output_reco_file, pipi_output_dir);
0350   if (run_ppi_reco) end_kfparticle(ppi_output_reco_file, pipi_output_dir);
0351   if (run_KK_reco) end_kfparticle(KK_output_reco_file, KK_output_dir);
0352
0353   gSystem->Exit(0);
0354
0355   return 0;
0356 }
0357
0358 #endif