KFPCutEff/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 #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(std::string processID = "0", std::string channel = "cascade")
0044 {
0045   int nEvents = 1e3;
0046   //int nEvents = 10;
0047   std::string outDir = "./" + channel + "_20260513_DetroitMB_CR_2_mode_pTref_1p4/";
0048
0049   string makeDirectory = "mkdir -p " + outDir + "hfEff";
0050   system(makeDirectory.c_str());
0051
0052   //F4A setup
0053   Fun4AllServer *se = Fun4AllServer::instance();
0054   se->Verbosity(1);
0055
0056   PHRandomSeed::Verbosity(1);
0057   recoConsts *rc = recoConsts::instance();
0058
0059   //Generator setup
0060
0061   Input::PYTHIA8 = true;
0062   int particleID = 421;
0063   PYTHIA8::config_file[0] = "steeringCards/pythia8_MB_Detroit_Tony.cfg";
0064   if (channel == "cascade_minus")
0065   {
0066     particleID = 3312;
0067   }
0068   else if (channel == "cascade_plus")
0069   {
0070     particleID = -3312;
0071   }
0072   else if (channel == "minBias")
0073   {
0074     std::cout << "Min bias simulations" << std::endl;
0075   }
0076   else
0077   {
0078     std::cout << "Your decay channel " << channel << " is not known" << std::endl;
0079     exit(1);
0080   }
0081   Input::BEAM_CONFIGURATION = Input::pp_COLLISION;
0082
0083   InputInit();
0084
0085   float abs_eta = 1.1;
0086   float pt_cut = 0.0;
0087
0088   if (channel != "minBias")
0089   {
0090
0091     vector<int> particleList;
0092     int lambdaID;
0093     if (channel == "cascade_minus")
0094     {
0095       particleList = {2212,-211};
0096       lambdaID = 3122;
0097     }
0098     else if (channel == "cascade_plus")
0099     {
0100       particleList = {-2212,211};
0101       lambdaID = -3122;
0102     }
0103
0104     for (unsigned int i = 0; i < particleList.size(); ++i)
0105     {
0106       string trigger_name = "particle_trigger_" + to_string(i);
0107       PHPy8ParticleTrigger * p8_hf_signal_trigger = new PHPy8ParticleTrigger(trigger_name.c_str());
0108       p8_hf_signal_trigger->SetEtaHighLow(abs_eta, -1*abs_eta);
0109       p8_hf_signal_trigger->SetPtLow(pt_cut);
0110       p8_hf_signal_trigger->SetStableParticleOnly(false);
0111       p8_hf_signal_trigger->AddParents(lambdaID);
0112       p8_hf_signal_trigger->AddParticles(particleList[i]);
0113       p8_hf_signal_trigger->SetParentRadialDecayVertexHigh(3.5);
0114       //p8_hf_signal_trigger->AddParticles(particleID);
0115       //p8_hf_signal_trigger->AddParticles(-1*particleID);
0116       p8_hf_signal_trigger->PrintConfig();
0117       INPUTGENERATOR::Pythia8[0]->register_trigger(p8_hf_signal_trigger);
0118     }
0119     PHPy8ParticleTrigger * p8_hf_signal_trigger_bachelor = new PHPy8ParticleTrigger("particle_trigger_bachelor");
0120     p8_hf_signal_trigger_bachelor->SetEtaHighLow(abs_eta, -1*abs_eta);
0121     p8_hf_signal_trigger_bachelor->SetPtLow(pt_cut);
0122     p8_hf_signal_trigger_bachelor->SetParentRadialDecayVertexHigh(3.5);
0123     p8_hf_signal_trigger_bachelor->SetStableParticleOnly(false);
0124     p8_hf_signal_trigger_bachelor->AddParents(particleID);
0125     if (channel == "cascade_minus")
0126     {
0127       p8_hf_signal_trigger_bachelor->AddParticles(-211);
0128     }
0129     else if (channel == "cascade_plus")
0130     {
0131       p8_hf_signal_trigger_bachelor->AddParticles(211);
0132     }
0133     //p8_hf_signal_trigger_bachelor->AddParticles(particleID);
0134     //p8_hf_signal_trigger_bachelor->AddParticles(-1*particleID);
0135     p8_hf_signal_trigger_bachelor->PrintConfig();
0136     INPUTGENERATOR::Pythia8[0]->register_trigger(p8_hf_signal_trigger_bachelor);
0137
0138     INPUTGENERATOR::Pythia8[0]->set_trigger_AND();
0139     //INPUTGENERATOR::Pythia8[0]->set_trigger_OR();
0140
0141     Input::ApplysPHENIXBeamParameter(INPUTGENERATOR::Pythia8);
0142   }
0143
0144   InputRegister();
0145
0146   //CDB flags and such
0147
0148   Enable::CDB = true;
0149   rc->set_StringFlag("CDB_GLOBALTAG",CDB::global_tag);
0150   rc->set_uint64Flag("TIMESTAMP",1);
0151   rc->set_IntFlag("RUNNUMBER",1);
0152
0153   Enable::MVTX_APPLYMISALIGNMENT = true;
0154   ACTSGEOM::mvtx_applymisalignment = Enable::MVTX_APPLYMISALIGNMENT;
0155
0156   SyncReco *sync = new SyncReco();
0157   se->registerSubsystem(sync);
0158
0159   HeadReco *head = new HeadReco();
0160   se->registerSubsystem(head);
0161
0162   FlagHandler *flag = new FlagHandler();
0163   se->registerSubsystem(flag);
0164
0165   //Simulation setup
0166   Enable::MBDFAKE = true;
0167   Enable::PIPE = true;
0168   Enable::PIPE_ABSORBER = true;
0169   Enable::MVTX = true;
0170   Enable::INTT = true;
0171   Enable::TPC = true;
0172   Enable::MICROMEGAS = true;
0173
0174   //Tracking setup
0175
0176   G4Init();
0177   MagnetInit();
0178   MagnetFieldInit();
0179
0180   G4Setup();
0181
0182   //Tagging stuff
0183   DecayFinder *myFinder = new DecayFinder("myFinder");
0184   myFinder->Verbosity(INT_MAX);
0185   if (channel == "cascade_minus" || channel == "cascade_plus") myFinder->setDecayDescriptor("[Xi- -> {Lambda0 -> proton^+ pi^-} pi^-]cc");
0186   myFinder->saveDST(1);
0187   myFinder->allowPi0(1);
0188   myFinder->allowPhotons(1);
0189   myFinder->triggerOnDecay(1);
0190   myFinder->setPTmin(pt_cut); //Note: sPHENIX min pT is 0.2 GeV for tracking
0191   myFinder->setEtaRange(-1*abs_eta, abs_eta); //Note: sPHENIX acceptance is |eta| <= 1.1
0192   myFinder->useDecaySpecificEtaRange(false);
0193   if (channel != "minBias") se->registerSubsystem(myFinder);
0194
0195   Mbd_Reco();
0196   Mvtx_Cells();
0197   Intt_Cells();
0198   TPC_Cells();
0199   Micromegas_Cells();
0200
0201   TrackingInit();
0202
0203   Mvtx_Clustering();
0204   Intt_Clustering();
0205   TPC_Clustering();
0206   Micromegas_Clustering();
0207
0208   Tracking_Reco();
0209
0210   auto vtxfinder = new PHSimpleVertexFinder;
0211   vtxfinder->Verbosity(0);
0212   vtxfinder->setDcaCut(1);
0213   vtxfinder->setTrackPtCut(-99999.);
0214   vtxfinder->setBeamLineCut(1);
0215   vtxfinder->setTrackQualityCut(1000000000);
0216   vtxfinder->setRequireMVTX(false);
0217   //vtxfinder->setNmvtxRequired(2);
0218   vtxfinder->setOutlierPairCut(1);
0219   se->registerSubsystem(vtxfinder);
0220
0221   Global_Reco();
0222
0223   build_truthreco_tables();
0224
0225   output_dir = outDir;
0226   if (true) create_hf_directories(Lambdapi_reconstruction_name, Lambdapi_output_dir, Lambdapi_output_reco_file, processID);
0227
0228   HFTrackEfficiency *myTrackEff = new HFTrackEfficiency("myTrackEff");
0229   myTrackEff->Verbosity(INT_MAX);
0230   myTrackEff->setDFNodeName("myFinder");
0231   myTrackEff->triggerOnDecay(1);
0232   myTrackEff->writeSelectedTrackMap(true);
0233   myTrackEff->writeOutputFile(true);
0234   std::string outputHFEffFile = outDir + "/hfEff/outputHFTrackEff_" + channel + "_" + processID + ".root";
0235   myTrackEff->setOutputFileName(outputHFEffFile);
0236   if (channel != "minBias" && false) se->registerSubsystem(myTrackEff);
0237
0238   std::string dir_eval = "TrackingEval";
0239   if ((channel == "cascade_minus" || channel == "cascade_plus") && false)
0240   {
0241     system(("mkdir -p " + Lambdapi_output_dir + "/" + dir_eval).c_str());
0242     Tracking_Eval(Form("%s/%s/outputTrackingEval_%s_%s.root", Lambdapi_output_dir.c_str(), dir_eval.c_str(), channel.c_str(), processID.c_str()));
0243   }
0244
0245   //if (run_pipi_reco) create_hf_directories(pipi_reconstruction_name, pipi_output_dir, pipi_output_reco_file);
0246   //if (run_ppi_reco) create_hf_directories(ppi_reconstruction_name, ppi_output_dir, ppi_output_reco_file);
0247   //if (run_anti_ppi_reco) create_hf_directories(anti_ppi_reconstruction_name, anti_ppi_output_dir, anti_ppi_output_reco_file);
0248
0249
0250   //if (run_pipi_reco) reconstruct_pipi_mass();
0251   //if (run_ppi_reco) reconstruct_ppi_mass();
0252   //if (run_anti_ppi_reco) reconstruct_ppi_mass();
0253   if (run_Lambdapi_reco) reconstruct_Lambdapi_mass();
0254 /*
0255   //Output file handling
0256   makeDirectory = "mkdir -p " + outDir + "DST";
0257   system(makeDirectory.c_str());
0258
0259   string FullOutFile = outDir + "/DST/" + channel + "_DST_" + processID + ".root";
0260   Fun4AllDstOutputManager *out = new Fun4AllDstOutputManager("DSTOUT", FullOutFile);
0261   out->StripNode("G4HIT_PIPE");
0262   out->StripNode("G4HIT_SVTXSUPPORT");
0263   //out->StripNode("PHG4INEVENT");
0264   //out->StripNode("Sync");
0265   out->StripNode("myFinder_DecayMap");
0266   out->StripNode("G4HIT_PIPE");
0267   out->StripNode("G4HIT_MVTX");
0268   out->StripNode("G4HIT_INTT");
0269   out->StripNode("G4HIT_TPC");
0270   out->StripNode("G4HIT_MICROMEGAS");
0271   out->StripNode("TRKR_HITSET");
0272   out->StripNode("TRKR_HITTRUTHASSOC");
0273   //out->StripNode("TRKR_CLUSTER");
0274   //out->StripNode("TRKR_CLUSTERHITASSOC");
0275   out->StripNode("TRKR_CLUSTERCROSSINGASSOC");
0276   out->StripNode("TRAINING_HITSET");
0277   out->StripNode("TRKR_TRUTHTRACKCONTAINER");
0278   out->StripNode("TRKR_TRUTHCLUSTERCONTAINER");
0279   out->StripNode("alignmentTransformationContainer");
0280   out->StripNode("alignmentTransformationContainerTransient");
0281   //out->StripNode("SiliconTrackSeedContainer");
0282   //out->StripNode("TpcTrackSeedContainer");
0283   //out->StripNode("SvtxTrackSeedContainer");
0284   out->StripNode("ActsTrajectories");
0285   //out->StripNode("SvtxTrackMap");
0286   out->StripNode("SvtxAlignmentStateMap");
0287   //out->SaveRunNode(0);
0288   se->registerOutputManager(out);
0289 */
0290   se->run(nEvents);
0291
0292   se->End();
0293
0294   //if (run_pipi_reco) end_kfparticle(pipi_output_reco_file, pipi_output_dir);
0295   //if (run_ppi_reco) end_kfparticle(ppi_output_reco_file, ppi_output_dir);
0296   //if (run_anti_ppi_reco) end_kfparticle(ppi_output_reco_file, ppi_output_dir);
0297   if (run_Lambdapi_reco) end_kfparticle(Lambdapi_output_reco_file, Lambdapi_output_dir);
0298
0299   gSystem->Exit(0);
0300
0301   return 0;
0302 }
0303
0304 #endif