KFParticle_sPHENIX/hf_generator/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
0018 #include <ffamodules/FlagHandler.h>
0019 #include <ffamodules/HeadReco.h>
0020 #include <ffamodules/SyncReco.h>
0021 #include <ffamodules/CDBInterface.h>
0022 #include <phool/PHRandomSeed.h>
0023 #include <phool/recoConsts.h>
0024
0025 #include <fun4all/Fun4AllDstOutputManager.h>
0026 #include <fun4all/Fun4AllOutputManager.h>
0027 #include <fun4all/Fun4AllServer.h>
0028
0029 #include <simqa_modules/QAG4SimulationTracking.h>
0030 #include <qautils/QAHistManagerDef.h>
0031
0032 R__LOAD_LIBRARY(libfun4all.so)
0033 R__LOAD_LIBRARY(libffamodules.so)
0034 R__LOAD_LIBRARY(libdecayfinder.so)
0035 R__LOAD_LIBRARY(libhftrackefficiency.so)
0036 R__LOAD_LIBRARY(libsimqa_modules.so)
0037
0038 int Fun4All_HFG(std::string processID = "0")
0039 {
0040   int nEvents = 2e2;
0041   std::string channel = "Ds2KKpi";
0042   std::string outDir = "/sphenix/tg/tg01/hf/cdean/selected_hf_simulations/" + channel;
0043   //F4A setup
0044
0045   Fun4AllServer *se = Fun4AllServer::instance();
0046   se->Verbosity(1);
0047
0048   PHRandomSeed::Verbosity(1);
0049   recoConsts *rc = recoConsts::instance();
0050
0051   //Generator setup
0052
0053   Input::PYTHIA8 = true;
0054   int particleID = 421;
0055
0056   if (channel == "bs2jpsiks0")
0057   {
0058     PYTHIA8::config_file = "steeringCards/pythia8_bs2jpsiks0.cfg";
0059     EVTGENDECAYER::DecayFile = "decFiles/Bs2JpsiKS0.DEC";
0060     particleID = 531;
0061   }
0062   else if (channel == "b2DX")
0063   {
0064     PYTHIA8::config_file = "steeringCards/pythia8_b2DX.cfg";
0065     EVTGENDECAYER::DecayFile = "decFiles/b2DX.DEC";
0066     particleID = 5;
0067   }
0068   else if (channel == "D2Kpi")
0069   {
0070     PYTHIA8::config_file = "steeringCards/pythia8_D2Kpi.cfg";
0071     EVTGENDECAYER::DecayFile = "decFiles/D2Kpi.DEC";
0072   }
0073   else if ((channel == "Dp2KKpi") || (channel == "Ds2KKpi"))
0074   {
0075     PYTHIA8::config_file = "steeringCards/pythia8_D2KKpi.cfg";
0076     EVTGENDECAYER::DecayFile = "decFiles/D2KKpi.DEC";
0077     if (channel == "Dp2KKpi") particleID = 411;
0078     else particleID = 431;
0079   }
0080   else if (channel == "Kshort2pipi")
0081   {
0082     PYTHIA8::config_file = "steeringCards/pythia8_K2pipi.cfg";
0083     EVTGENDECAYER::DecayFile = "decFiles/K2pipi.DEC";
0084     particleID = 310;
0085   }
0086   else if (channel == "minBias")
0087   {
0088     std::cout << "Min bias simulations" << std::endl;
0089   }
0090   else
0091   {
0092     std::cout << "Your decay channel " << channel << " is not known" << std::endl;
0093     exit(1);
0094   }
0095
0096   Input::BEAM_CONFIGURATION = Input::pp_COLLISION;
0097
0098   InputInit();
0099
0100   PHPy8ParticleTrigger * p8_hf_signal_trigger = new PHPy8ParticleTrigger();
0101   p8_hf_signal_trigger->SetPtLow(0.);
0102   p8_hf_signal_trigger->SetEtaHighLow(1.3, -1.3); // sample a rapidity range higher than the sPHENIX tracking pseudorapidity
0103   p8_hf_signal_trigger->SetStableParticleOnly(false); // process unstable particles that include quarks
0104   p8_hf_signal_trigger->PrintConfig();
0105
0106   if (channel != "minBias")
0107   {
0108     p8_hf_signal_trigger->AddParticles(particleID);
0109     p8_hf_signal_trigger->AddParticles(-1*particleID);
0110     INPUTGENERATOR::Pythia8->register_trigger(p8_hf_signal_trigger);
0111     INPUTGENERATOR::Pythia8->set_trigger_OR();
0112
0113     Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8);
0114   }
0115
0116   InputRegister();
0117
0118   //CDB flags and such
0119
0120   rc->set_IntFlag("RUNNUMBER",1);
0121
0122   SyncReco *sync = new SyncReco();
0123   se->registerSubsystem(sync);
0124
0125   HeadReco *head = new HeadReco();
0126   se->registerSubsystem(head);
0127
0128   FlagHandler *flag = new FlagHandler();
0129   se->registerSubsystem(flag);
0130
0131   DecayFinder *myFinder = new DecayFinder("myFinder");
0132   myFinder->Verbosity(INT_MAX);
0133   if (channel == "bs2jpsiks0") myFinder->setDecayDescriptor("[B_s0 -> {J/psi -> e^+ e^-} {K_S0 -> pi^+ pi^-}]cc");
0134   else if  (channel == "Kshort2pipi") myFinder->setDecayDescriptor("K_S0 -> pi^- pi^+");
0135   else if  (channel == "Dp2KKpi") myFinder->setDecayDescriptor("[D+ -> {phi -> K^+ K^-} pi^+]cc");
0136   else if  (channel == "Ds2KKpi") myFinder->setDecayDescriptor("[D_s+ -> {phi -> K^+ K^-} pi^+]cc");
0137   else myFinder->setDecayDescriptor("[D0 -> K^- pi^+]cc");
0138   myFinder->saveDST(1);
0139   myFinder->allowPi0(1);
0140   myFinder->allowPhotons(1);
0141   myFinder->triggerOnDecay(1);
0142   myFinder->setPTmin(0.16); //Note: sPHENIX min pT is 0.2 GeV for tracking
0143   myFinder->setEtaRange(-1.2, 1.2); //Note: sPHENIX acceptance is |eta| <= 1.1
0144   myFinder->useDecaySpecificEtaRange(false);
0145   if (channel != "minBias") se->registerSubsystem(myFinder);
0146
0147   //Simulation setup
0148
0149   Enable::MBDFAKE = true;
0150   Enable::PIPE = true;
0151   Enable::PIPE_ABSORBER = true;
0152   Enable::MVTX = true;
0153   Enable::INTT = true;
0154   Enable::TPC = true;
0155   Enable::MICROMEGAS = true;
0156
0157   //Tracking setup
0158
0159   Enable::CDB = true;
0160   rc->set_StringFlag("CDB_GLOBALTAG",CDB::global_tag);
0161   rc->set_uint64Flag("TIMESTAMP",CDB::timestamp);
0162
0163   MagnetInit();
0164   MagnetFieldInit();
0165
0166   G4Init();
0167   G4Setup();
0168
0169   Mbd_Reco();
0170   Mvtx_Cells();
0171   Intt_Cells();
0172   TPC_Cells();
0173   Micromegas_Cells();
0174
0175   TrackingInit();
0176
0177   Mvtx_Clustering();
0178   Intt_Clustering();
0179   TPC_Clustering();
0180   Micromegas_Clustering();
0181
0182   Tracking_Reco();
0183   Global_Reco();
0184
0185   build_truthreco_tables();
0186
0187   //HF stuff
0188
0189   HFTrackEfficiency *myTrackEff = new HFTrackEfficiency("myTrackEff");
0190   myTrackEff->Verbosity(INT_MAX);
0191   //myTrackEff->setTruthRecoMatchingPercentage(100.);
0192   myTrackEff->setDFNodeName("myFinder");
0193   myTrackEff->triggerOnDecay(1);
0194   myTrackEff->writeSelectedTrackMap(true);
0195   myTrackEff->writeOutputFile(true);
0196   std::string outputHFEffFile = outDir + "/hfEff/outputHFTrackEff_" + channel + "_" + processID + ".root";
0197   myTrackEff->setOutputFileName(outputHFEffFile);
0198   if (channel != "minBias") se->registerSubsystem(myTrackEff);
0199
0200   //Output file handling
0201
0202   string FullOutFile = outDir + "/" + channel + "_DST_" + processID + ".root";
0203   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
0204   out->StripNode("G4HIT_PIPE");
0205   out->StripNode("G4HIT_SVTXSUPPORT");
0206   out->StripNode("PHG4INEVENT");
0207   out->StripNode("Sync");
0208   out->StripNode("myFinder_DecayMap");
0209   out->StripNode("G4HIT_PIPE");
0210   out->StripNode("G4HIT_MVTX");
0211   out->StripNode("G4HIT_INTT");
0212   out->StripNode("G4HIT_TPC");
0213   out->StripNode("G4HIT_MICROMEGAS");
0214   out->StripNode("TRKR_HITSET");
0215   out->StripNode("TRKR_HITTRUTHASSOC");
0216   //out->StripNode("TRKR_CLUSTER");
0217   out->StripNode("TRKR_CLUSTERHITASSOC");
0218   out->StripNode("TRKR_CLUSTERCROSSINGASSOC");
0219   out->StripNode("TRAINING_HITSET");
0220   out->StripNode("TRKR_TRUTHTRACKCONTAINER");
0221   out->StripNode("TRKR_TRUTHCLUSTERCONTAINER");
0222   out->StripNode("alignmentTransformationContainer");
0223   out->StripNode("alignmentTransformationContainerTransient");
0224   //out->StripNode("SiliconTrackSeedContainer");
0225   //out->StripNode("TpcTrackSeedContainer");
0226   //out->StripNode("SvtxTrackSeedContainer");
0227   out->StripNode("ActsTrajectories");
0228   //out->StripNode("SvtxTrackMap");
0229   out->StripNode("SvtxAlignmentStateMap");
0230   out->SaveRunNode(0);
0231   se->registerOutputManager(out);
0232
0233   se->run(nEvents);
0234
0235   se->End();
0236   gSystem->Exit(0);
0237
0238   return 0;
0239 }
0240
0241 #endif