MC-Calibrations/macro/Fun4All_CaloSimAna.C

0001 #ifndef MACRO_FUN4ALL_JESMCTREEGEN_C
0002 #define MACRO_FUN4ALL_JESMCTREEGEN_C
0003
0004 #include <ffamodules/CDBInterface.h>
0005 #include <ffamodules/FlagHandler.h>
0006 #include <ffamodules/HeadReco.h>
0007 #include <ffamodules/SyncReco.h>
0008
0009 #include <fun4all/Fun4AllInputManager.h>
0010 #include <fun4all/Fun4AllOutputManager.h>
0011 #include <fun4all/Fun4AllDstInputManager.h>
0012 #include <fun4all/Fun4AllDstOutputManager.h>
0013 #include <fun4all/Fun4AllServer.h>
0014 #include <fun4all/SubsysReco.h>
0015
0016 #include <phool/recoConsts.h>
0017 #include <phool/PHRandomSeed.h>
0018 #include <phool/recoConsts.h>
0019 #include <G4_Input.C>
0020 #include <GlobalVariables.C>
0021
0022 #include <caloreco/RawClusterCNNClassifier.h>
0023 #include <clusteriso/ClusterIso.h>
0024
0025 #include <jetbase/FastJetAlgo.h>
0026 #include <jetbase/JetReco.h>
0027 #include <jetbase/TowerJetInput.h>
0028 #include <g4jets/TruthJetInput.h>
0029
0030 #include <jetbackground/CopyAndSubtractJets.h>
0031 #include <jetbackground/DetermineTowerBackground.h>
0032 #include <jetbackground/FastJetAlgoSub.h>
0033 #include <jetbackground/RetowerCEMC.h>
0034 #include <jetbackground/SubtractTowers.h>
0035 #include <jetbackground/SubtractTowersCS.h>
0036
0037 #include <calotrigger/CaloTriggerEmulator.h>
0038
0039 #include <JESMCTreeGen/JESMCTreeGen.h>
0040
0041 R__LOAD_LIBRARY(libfun4all.so)
0042 R__LOAD_LIBRARY(libfun4allraw.so)
0043 R__LOAD_LIBRARY(libclusteriso.so)
0044 R__LOAD_LIBRARY(libcalo_reco.so)
0045 R__LOAD_LIBRARY(libffamodules.so)
0046 R__LOAD_LIBRARY(libffarawobjects.so)
0047 R__LOAD_LIBRARY(libJESMCTreeGen.so)
0048 R__LOAD_LIBRARY(libjetbase.so)
0049 R__LOAD_LIBRARY(libg4jets.so)
0050 R__LOAD_LIBRARY(libjetbackground.so)
0051
0052 void Fun4All_CaloSimAna(
0053     const int nEvents = 0,
0054     const string &inputFile0 = "example_list/g4hits.list",
0055     const string &inputFile1 = "example_list/dst_calo_waveform.list",
0056     const string &inputFile2 = "example_list/dst_global.list",
0057     const string &inputFile3 = "example_list/dst_truth_jet.list",
0058     const string &outputFile = "output_sim.root")
0059 {
0060   // ********** Setting ********** //
0061   int verbosity = 0;
0062   bool doZVertexReadout = true;
0063   bool doTowerInfoReadout = false;
0064     bool doTowerInfoReadout_RetowerEMCal = false;
0065   bool doClusterReadout = false;
0066     bool doClusterIsoReadout = false;
0067     bool doClusterCNNReadout = false;
0068   bool doRecoJetReadout = true;
0069     std::vector<float> doRecoJet_radius = {0.4};
0070     bool doSeedJetRawReadout = false;
0071     bool doSeedJetSubReadout = false;
0072   bool doTruthJetReadout = true;
0073     std::vector<float> doTruthJet_radius = {0.4};
0074   bool doTruthParticleReadout = false;
0075
0076   // ********** Initialization for Fun4All ********** //
0077   Fun4AllServer *se = Fun4AllServer::instance();
0078   se->Verbosity(verbosity);
0079   recoConsts *rc = recoConsts::instance();
0080
0081   Enable::CDB = true;
0082   rc->set_StringFlag("CDB_GLOBALTAG", "ProdA_2024");
0083   rc->set_uint64Flag("TIMESTAMP", 47289);
0084   CDBInterface::instance()->Verbosity(1);
0085
0086   FlagHandler *flag = new FlagHandler();
0087   se->registerSubsystem(flag);
0088
0089   // ********** Register subsystems ********** //
0090   // Cluster
0091   ClusterIso *cliso4;
0092   if (doClusterIsoReadout) {
0093     cliso4 = new ClusterIso("ClusterIso4",1,4,false,true);
0094     cliso4->Verbosity(0);
0095     cliso4->set_use_towerinfo(true);
0096     se->registerSubsystem(cliso4);
0097   }
0098
0099   RawClusterCNNClassifier *cnn;
0100   if (doClusterCNNReadout) {
0101     cnn = new RawClusterCNNClassifier("CNN_single");
0102     cnn->set_modelPath("/sphenix/u/shuhang98/core_patch/coresoftware/offline/packages/CaloReco/functional_model_single.onnx");
0103     cnn->set_inputNodeName("CLUSTERINFO_CEMC");
0104     cnn->set_outputNodeName("CLUSTERINFO_CEMC_CNN");
0105     se->registerSubsystem(cnn);
0106   }
0107
0108   // JetReco
0109   std::string jetreco_input_prefix = "TOWERINFO_CALIB";
0110   RetowerCEMC *_retowerCEMC;
0111   JetReco *_jetReco_seed;
0112   DetermineTowerBackground *_determineTowerBackground;
0113   CopyAndSubtractJets *_copyAndSubtractJets;
0114   DetermineTowerBackground *_determineTowerBackground2;
0115   SubtractTowers *_subtractTowers;
0116   JetReco *_jetReco;
0117   if (doRecoJetReadout) {
0118     _retowerCEMC = new RetowerCEMC();
0119     _retowerCEMC->Verbosity(verbosity);
0120     _retowerCEMC->set_towerinfo(true);
0121     _retowerCEMC->set_frac_cut(0.5); //fraction of retower that must be masked to mask the full retower
0122     _retowerCEMC->set_towerNodePrefix(jetreco_input_prefix);
0123     se->registerSubsystem(_retowerCEMC);
0124
0125     _jetReco_seed = new JetReco();
0126     _jetReco_seed->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER,jetreco_input_prefix));
0127     _jetReco_seed->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO,jetreco_input_prefix));
0128     _jetReco_seed->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO,jetreco_input_prefix));
0129     _jetReco_seed->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.2), "AntiKt_TowerInfo_HIRecoSeedsRaw_r02");
0130     _jetReco_seed->set_algo_node("ANTIKT");
0131     _jetReco_seed->set_input_node("TOWER");
0132     _jetReco_seed->Verbosity(verbosity);
0133     se->registerSubsystem(_jetReco_seed);
0134
0135     _determineTowerBackground = new DetermineTowerBackground();
0136     _determineTowerBackground->SetBackgroundOutputName("TowerInfoBackground_Sub1");
0137     _determineTowerBackground->SetFlow(false);
0138     _determineTowerBackground->SetSeedType(0);
0139     _determineTowerBackground->SetSeedJetD(3);
0140     _determineTowerBackground->set_towerinfo(true);
0141     _determineTowerBackground->Verbosity(verbosity);
0142     _determineTowerBackground->set_towerNodePrefix(jetreco_input_prefix);
0143     se->registerSubsystem(_determineTowerBackground);
0144
0145     _copyAndSubtractJets = new CopyAndSubtractJets();
0146     _copyAndSubtractJets->SetFlowModulation(false);
0147     _copyAndSubtractJets->Verbosity(verbosity);
0148     _copyAndSubtractJets->set_towerinfo(true);
0149     _copyAndSubtractJets->set_towerNodePrefix(jetreco_input_prefix);
0150     se->registerSubsystem(_copyAndSubtractJets);
0151
0152     _determineTowerBackground2 = new DetermineTowerBackground();
0153     _determineTowerBackground2->SetBackgroundOutputName("TowerInfoBackground_Sub2");
0154     _determineTowerBackground2->SetFlow(false);
0155     _determineTowerBackground2->SetSeedType(1);
0156     _determineTowerBackground2->SetSeedJetPt(7);
0157     _determineTowerBackground2->Verbosity(verbosity);
0158     _determineTowerBackground2->set_towerinfo(true);
0159     _determineTowerBackground2->set_towerNodePrefix(jetreco_input_prefix);
0160     se->registerSubsystem(_determineTowerBackground2);
0161
0162     _subtractTowers = new SubtractTowers();
0163     _subtractTowers->SetFlowModulation(false);
0164     _subtractTowers->Verbosity(verbosity);
0165     _subtractTowers->set_towerinfo(true);
0166     _subtractTowers->set_towerNodePrefix(jetreco_input_prefix);
0167     se->registerSubsystem(_subtractTowers);
0168
0169     _jetReco = new JetReco();
0170     _jetReco->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO_SUB1, jetreco_input_prefix));
0171     _jetReco->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO_SUB1, jetreco_input_prefix));
0172     _jetReco->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO_SUB1, jetreco_input_prefix));
0173     _jetReco->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.4), "AntiKt_subtracted_r04");
0174     _jetReco->set_algo_node("ANTIKT");
0175     _jetReco->set_input_node("TOWER");
0176     _jetReco->Verbosity(verbosity);
0177     se->registerSubsystem(_jetReco);
0178   }
0179
0180   // ********** Register input managers ********** //
0181   Fun4AllInputManager *in0 = new Fun4AllDstInputManager("in0");
0182   in0->AddListFile(inputFile0,1);
0183   se->registerInputManager(in0);
0184
0185   Fun4AllInputManager *in1 = new Fun4AllDstInputManager("in1");
0186   in1->AddListFile(inputFile1,1);
0187   se->registerInputManager(in1);
0188
0189   Fun4AllInputManager *in2 = new Fun4AllDstInputManager("in2");
0190   in2->AddListFile(inputFile2,1);
0191   se->registerInputManager(in2);
0192
0193   Fun4AllInputManager *in3 = new Fun4AllDstInputManager("in3");
0194   in3->AddListFile(inputFile3,1);
0195   se->registerInputManager(in3);
0196
0197
0198   // ********** Register output managers ********** //
0199   JESMCTreeGen *_JESMCTreeGen = new JESMCTreeGen("JESMCTreeGen",outputFile);
0200   _JESMCTreeGen->SetVerbosity(verbosity);
0201   _JESMCTreeGen->DoZVertexReadout(doZVertexReadout);
0202   _JESMCTreeGen->DoTowerInfoReadout(doTowerInfoReadout);
0203   _JESMCTreeGen->DoTowerInfoReadout_RetowerEMCal(doTowerInfoReadout_RetowerEMCal);
0204   _JESMCTreeGen->DoClusterReadout(doClusterReadout);
0205   _JESMCTreeGen->DoClusterIsoReadout(doClusterIsoReadout);
0206   _JESMCTreeGen->DoClusterCNNReadout(doClusterCNNReadout);
0207   _JESMCTreeGen->DoRecoJetReadout(doRecoJetReadout);
0208   _JESMCTreeGen->DoRecoJet_radius(doRecoJet_radius);
0209   _JESMCTreeGen->DoSeedJetRawReadout(doSeedJetRawReadout);
0210   _JESMCTreeGen->DoSeedJetSubReadout(doSeedJetSubReadout);
0211   _JESMCTreeGen->DoTruthJetReadout(doTruthJetReadout);
0212   _JESMCTreeGen->DoTruthJet_radius(doTruthJet_radius);
0213   _JESMCTreeGen->DoTruthParticleReadout(doTruthParticleReadout);
0214   se->registerSubsystem(_JESMCTreeGen);
0215
0216   se->run(nEvents);
0217   CDBInterface::instance()->Print();  // print used DB files
0218   se->End();
0219   cout << "JOB COMPLETE." <<endl;
0220   se->PrintTimer();
0221   gSystem->Exit(0);
0222 }
0223 #endif