JetBkgdSub/macro/Fun4All_JetBkgd_Embed.C

0001 #pragma once
0002 #if ROOT_VERSION_CODE >= ROOT_VERSION(6,00,0)
0003 #include <fun4all/SubsysReco.h>
0004 #include <fun4all/Fun4AllServer.h>
0005 #include <fun4all/Fun4AllInputManager.h>
0006 #include <fun4all/Fun4AllSyncManager.h>
0007 #include <fun4all/Fun4AllDstInputManager.h>
0008 #include <fun4all/Fun4AllNoSyncDstInputManager.h>
0009
0010 #include <fun4all/Fun4AllDstOutputManager.h>
0011 #include <fun4all/Fun4AllOutputManager.h>
0012
0013 #include <phool/PHRandomSeed.h>
0014 #include <phool/recoConsts.h>
0015
0016 #include <centrality/CentralityReco.h>
0017
0018 #include <caloreco/CaloTowerStatus.h>
0019
0020 #include <jetbackground/RetowerCEMC.h>
0021 #include <jetbackground/FastJetAlgoSub.h>
0022
0023 #include <jetbase/JetReco.h>
0024 #include <jetbase/TowerJetInput.h>
0025 #include <jetbase/FastJetAlgo.h>
0026 #include <g4jets/TruthJetInput.h>
0027
0028 #include <caloembedding/caloTowerEmbed.h>
0029
0030 #include <jetbkgdsub/JetBkgdSub.h>
0031
0032 R__LOAD_LIBRARY(libfun4all.so)
0033 R__LOAD_LIBRARY(libg4jets.so)
0034 R__LOAD_LIBRARY(libjetbackground.so)
0035 R__LOAD_LIBRARY(libjetbase.so)
0036 R__LOAD_LIBRARY(libg4dst.so)
0037 R__LOAD_LIBRARY(libJetBkgdSub.so)
0038 R__LOAD_LIBRARY(libcentrality_io.so)
0039 R__LOAD_LIBRARY(libcentrality.so)
0040 R__LOAD_LIBRARY(libCaloEmbedding.so)
0041 R__LOAD_LIBRARY(libcalo_reco.so)
0042
0043
0044 #endif
0045
0046
0047 void Fun4All_JetBkgd_Embed(
0048   const char *filelistdata = "dst_data.list",
0049   const char *filelistsim = "dst_sim.list",
0050   const char * output = "output",
0051   const double jet_parameter = 0.4
0052 )
0053 {
0054
0055   std::string outputbase = output;
0056
0057   bool doCEMC = true;
0058   bool doIHCAL = true;
0059   bool doOHCAL = true;
0060
0061   //-----------------------------------
0062   // Fun4All server initialization
0063   //-----------------------------------
0064
0065   // create fun4all server
0066   Fun4AllServer *se = Fun4AllServer::instance();
0067   int verbosity = 1;
0068   se->Verbosity(verbosity);
0069   recoConsts *rc = recoConsts::instance();
0070
0071   rc->set_StringFlag("CDB_GLOBALTAG","ProdA_2023");
0072   rc->set_uint64Flag("TIMESTAMP",23745);
0073
0074
0075   CentralityReco *cent = new CentralityReco("CentralityReco");
0076   se->registerSubsystem(cent);
0077
0078
0079   /*
0080   // retower CEMC
0081   RetowerCEMC *rcemc = new RetowerCEMC();
0082   rcemc->set_towerinfo(true);
0083   rcemc->Verbosity(verbosity);
0084   se->registerSubsystem(rcemc);
0085   */
0086
0087   /*
0088   //-----------------------------------
0089   // Jet reco
0090   //-----------------------------------
0091   Enable::HIJETS_TRUTH=false;
0092   Enable::HIJETS_VERBOSITY=0;
0093   HIJetReco();
0094
0095   // tower jets
0096   // create jetreco and jettruth node names
0097   string rawname = "AntiKt_Tower_r0" + to_string(int(jet_parameter * 10));
0098   JetReco *towerjetreco = new JetReco();
0099   towerjetreco->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER));
0100   towerjetreco->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO));
0101   towerjetreco->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO));
0102   towerjetreco->add_algo(new FastJetAlgo(Jet::ANTIKT, jet_parameter), rawname);
0103   towerjetreco->set_algo_node("ANTIKT");
0104   towerjetreco->set_input_node("TOWER");
0105   towerjetreco->Verbosity(verbosity);
0106   se->registerSubsystem(towerjetreco);
0107
0108   */
0109
0110
0111   if(doCEMC){
0112     CaloTowerStatus *statusEMC = new CaloTowerStatus("CEMCSTATUS");
0113     statusEMC->set_detector_type(CaloTowerDefs::CEMC);
0114     statusEMC->set_time_cut(1);
0115     statusEMC->set_inputNodePrefix("TOWERINFO_CALIB_");
0116     se->registerSubsystem(statusEMC);
0117
0118     caloTowerEmbed *embedder_CEMC = new caloTowerEmbed("embedder_CEMC");
0119     embedder_CEMC->set_detector_type(CaloTowerDefs::CEMC);
0120     embedder_CEMC->set_removeBadTowers(true);
0121     embedder_CEMC->set_inputNodePrefix("TOWERINFO_CALIB_");
0122     se->registerSubsystem(embedder_CEMC);
0123   }
0124
0125   if(doIHCAL){
0126     CaloTowerStatus *statusIHCAL = new CaloTowerStatus("IHCALSTATUS");
0127     statusIHCAL->set_detector_type(CaloTowerDefs::HCALIN);
0128     statusIHCAL->set_time_cut(2);
0129     statusIHCAL->set_inputNodePrefix("TOWERINFO_CALIB_");
0130     se->registerSubsystem(statusIHCAL);
0131
0132     caloTowerEmbed *embedder_IHCAL = new caloTowerEmbed("embedder_IHCAL");
0133     embedder_IHCAL->set_detector_type(CaloTowerDefs::HCALIN);
0134     embedder_IHCAL->set_removeBadTowers(true);
0135     embedder_IHCAL->set_inputNodePrefix("TOWERINFO_CALIB_");
0136     se->registerSubsystem(embedder_IHCAL);
0137   }
0138
0139   if(doOHCAL){
0140     caloTowerEmbed *embedder_OHCAL = new caloTowerEmbed("embedder_OHCal");
0141     embedder_OHCAL->set_detector_type(CaloTowerDefs::HCALOUT);
0142     embedder_OHCAL->set_removeBadTowers(true);
0143     embedder_OHCAL->set_inputNodePrefix("TOWERINFO_CALIB_");
0144     se->registerSubsystem(embedder_OHCAL);
0145
0146     CaloTowerStatus *statusOHCAL = new CaloTowerStatus("OHCALSTATUS");
0147     statusOHCAL->set_detector_type(CaloTowerDefs::HCALOUT);
0148     statusOHCAL->set_time_cut(2);
0149     statusOHCAL->set_inputNodePrefix("TOWERINFO_CALIB_");
0150     se->registerSubsystem(statusOHCAL);
0151   }
0152
0153   // ==============
0154   // Jet Bkgd Sub
0155   // ==============
0156   double etamin = -1.1 + jet_parameter;
0157   double etamax = 1.1 - jet_parameter;
0158   //data
0159   JetBkgdSub *myJetTree = new JetBkgdSub(jet_parameter,outputbase + ".root");
0160   if(doCEMC) myJetTree->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO));
0161   if(doIHCAL) myJetTree->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO));
0162   if(doOHCAL) myJetTree->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO));
0163   myJetTree->doIterative(false);
0164   myJetTree->doAreaSub(true);
0165   myJetTree->doMultSub(false);
0166   myJetTree->doTruth(false);
0167   myJetTree->doData(true);
0168   myJetTree->doEmbed(false);
0169   myJetTree->doTowerECut(true);
0170   myJetTree->setTowerThreshold(0.05);
0171   myJetTree->setMinRecoPt(5.0); // only sets range for reco jets
0172   myJetTree->setEtaRange(etamin, etamax);
0173   myJetTree->setPtRange(10, 100); // only sets range for truth jets
0174   myJetTree->Verbosity(verbosity);
0175   se->registerSubsystem(myJetTree);
0176
0177
0178   //embed
0179   JetBkgdSub *myJetTreeEmbed = new JetBkgdSub(jet_parameter,outputbase + "_embed.root");
0180   if(doCEMC) myJetTreeEmbed->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO_EMBED));
0181   if(doIHCAL) myJetTreeEmbed->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO_EMBED));
0182   if(doOHCAL) myJetTreeEmbed->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO_EMBED));
0183   myJetTreeEmbed->doIterative(false);
0184   myJetTreeEmbed->doAreaSub(true);
0185   myJetTreeEmbed->doMultSub(false);
0186   myJetTreeEmbed->doTruth(false);
0187   myJetTreeEmbed->doData(true);
0188   myJetTreeEmbed->doEmbed(true);
0189   myJetTreeEmbed->doTowerECut(true);
0190   myJetTreeEmbed->setTowerThreshold(0.05);
0191   myJetTreeEmbed->setMinRecoPt(5.0); // only sets range for reco jets
0192   myJetTreeEmbed->setEtaRange(etamin, etamax);
0193   myJetTreeEmbed->setPtRange(10, 100); // only sets range for truth jets
0194   myJetTreeEmbed->Verbosity(verbosity);
0195   se->registerSubsystem(myJetTreeEmbed);
0196
0197   //sim
0198   JetBkgdSub *myJetTreeSim = new JetBkgdSub(jet_parameter,outputbase + "_sim.root");
0199   if(doCEMC) myJetTreeSim->add_input(new TowerJetInput(Jet::CEMC_TOWERINFO_SIM));
0200   if(doIHCAL) myJetTreeSim->add_input(new TowerJetInput(Jet::HCALIN_TOWERINFO_SIM));
0201   if(doOHCAL) myJetTreeSim->add_input(new TowerJetInput(Jet::HCALOUT_TOWERINFO_SIM));
0202   myJetTreeSim->doIterative(false);
0203   myJetTreeSim->doAreaSub(true);
0204   myJetTreeSim->doMultSub(false);
0205   myJetTreeSim->doTruth(false);
0206   myJetTreeSim->doData(false);
0207   myJetTreeSim->doEmbed(false);
0208   myJetTreeSim->doTowerECut(true);
0209   myJetTreeSim->setTowerThreshold(0.05);
0210   myJetTreeSim->setMinRecoPt(5.0); // only sets range for reco jets
0211   myJetTreeSim->setEtaRange(etamin, etamax);
0212   myJetTreeSim->setPtRange(10, 100); // only sets range for truth jets
0213   myJetTreeSim->Verbosity(verbosity);
0214   se->registerSubsystem(myJetTreeSim);
0215
0216   //-----------------------------------
0217   // Input managers
0218   //-----------------------------------
0219
0220   Fun4AllSyncManager *sync = se->getSyncManager();
0221   sync->MixRunsOk(true);
0222
0223   Fun4AllInputManager *inData = new Fun4AllDstInputManager("DSTData","DST");
0224   inData->AddListFile(filelistdata);
0225   se->registerInputManager(inData);
0226
0227   Fun4AllInputManager *inSim = new Fun4AllNoSyncDstInputManager("DSTSim","DST","TOPSim");
0228   inSim->AddListFile(filelistsim);
0229   se->registerInputManager(inSim);
0230   //-----------------------------------
0231   // Run the analysis
0232   //-----------------------------------
0233
0234   se->run(100);
0235   se->End();
0236
0237   gSystem->Exit(0);
0238   return 0;
0239
0240 }