macros/JetProduction/Fun4All_JetSkimmedProductionYear2.C

0001 #ifndef FUN4ALL_JETSKIMMEDPRODUCTIONYEAR2_C
0002 #define FUN4ALL_JETSKIMMEDPRODUCTIONYEAR2_C
0003
0004 #include <QA.C>
0005 #include <Calo_Calib.C>
0006 #include <HIJetReco.C>
0007 #include <Jet_QA.C>
0008 #include <G4_Global.C>
0009 #include <G4_Centrality.C>
0010
0011 #include <mbd/MbdReco.h>
0012
0013 #include <globalvertex/GlobalVertexReco.h>
0014
0015 #include <ffamodules/CDBInterface.h>
0016 #include <ffamodules/FlagHandler.h>
0017 #include <ffamodules/HeadReco.h>
0018 #include <ffamodules/SyncReco.h>
0019
0020 #include <fun4allraw/Fun4AllPrdfInputManager.h>
0021
0022 #include <fun4all/Fun4AllDstInputManager.h>
0023 #include <fun4all/Fun4AllDstOutputManager.h>
0024 #include <fun4all/Fun4AllInputManager.h>
0025 #include <fun4all/Fun4AllRunNodeInputManager.h>
0026 #include <fun4all/Fun4AllServer.h>
0027 #include <fun4all/Fun4AllUtils.h>
0028 #include <fun4all/SubsysReco.h>
0029
0030 #include <phool/recoConsts.h>
0031
0032 #include <centrality/CentralityReco.h>
0033 #include <calotrigger/MinimumBiasClassifier.h>
0034
0035 #include <calovalid/CaloValid.h>
0036
0037 #include <jetdstskimmer/JetDSTSkimmer.h>
0038
0039
0040 R__LOAD_LIBRARY(libfun4all.so)
0041 R__LOAD_LIBRARY(libfun4allraw.so)
0042 R__LOAD_LIBRARY(libcalo_reco.so)
0043 R__LOAD_LIBRARY(libcalotrigger.so)
0044 R__LOAD_LIBRARY(libcentrality.so)
0045 R__LOAD_LIBRARY(libffamodules.so)
0046 R__LOAD_LIBRARY(libmbd.so)
0047 R__LOAD_LIBRARY(libglobalvertex.so)
0048 R__LOAD_LIBRARY(libcalovalid.so)
0049 R__LOAD_LIBRARY(libJetDSTSkimmer.so)
0050
0051 void Fun4All_JetSkimmedProductionYear2(int nEvents=1000,
0052                         const std::string &fname = "DST_CALOFITTING_run2pp_ana509_2024p022_v001-00047289-00000.root",
0053                         const std::string& outfile_low= "DST_JETCALO_run2pp_ana509_2024p022_v001-00047289-00000.root",
0054                         const std::string& outfile_high= "DST_Jet_run2pp_ana509_2024p022_v001-00047289-00000.root",
0055                         const std::string& outfile_hist= "HIST_JETQA_run2pp_ana509_2024p022_v001-00047289-00000.root",
0056                         const std::string& dbtag= "ProdA_2024"
0057   )
0058 {
0059
0060   Fun4AllServer *se = Fun4AllServer::instance();
0061   se->Verbosity(1);
0062
0063   recoConsts *rc = recoConsts::instance();
0064
0065   std::pair<int, int> runseg = Fun4AllUtils::GetRunSegment(fname);
0066   int runnumber = runseg.first;
0067
0068   // conditions DB flags and timestamp
0069   rc->set_StringFlag("CDB_GLOBALTAG", dbtag);
0070   rc->set_uint64Flag("TIMESTAMP", runnumber);
0071   CDBInterface::instance()->Verbosity(1);
0072
0073   FlagHandler *flag = new FlagHandler();
0074   se->registerSubsystem(flag);
0075
0076   // MBD/BBC Reconstruction
0077   MbdReco *mbdreco = new MbdReco();
0078   se->registerSubsystem(mbdreco);
0079
0080   // Official vertex storage
0081   GlobalVertexReco *gvertex = new GlobalVertexReco();
0082   se->registerSubsystem(gvertex);
0083
0084   /////////////////////////////////////////////////////
0085   // Set status of CALO towers, Calibrate towers,  Cluster
0086   std::cout << "Processing Calo Calibration" << std::endl;
0087   Process_Calo_Calib();
0088
0089   ///////////////////////////////////
0090   // Validation maybe no need to run for the skimmed production?
0091   /*
0092   CaloValid *ca = new CaloValid("CaloValid");
0093   ca->set_timing_cut_width(200);
0094   se->registerSubsystem(ca);
0095   */
0096
0097   // Jet reco related flags
0098   Enable::QA = true;
0099
0100   // turn on pp mode
0101   HIJETS::is_pp = true;
0102
0103   // qa options
0104   JetQA::HasTracks = false;
0105   JetQA::DoInclusive = true;
0106   JetQA::DoTriggered = true;
0107   JetQA::RestrictPtToTrig = false;
0108   JetQA::RestrictEtaByR = true;
0109
0110   if (!HIJETS::is_pp)
0111   {
0112     Centrality();
0113   }
0114
0115   GlobalVertex::VTXTYPE vertex_type = GlobalVertex::MBD;
0116
0117   std::string tower_prefix = "TOWERINFO_CALIB";
0118
0119   RetowerCEMC *rcemc = new RetowerCEMC();
0120   rcemc->Verbosity(0);
0121   rcemc->set_towerinfo(true);
0122   rcemc->set_frac_cut(0.5);  // fraction of retower that must be masked to mask the full retower
0123   rcemc->set_towerNodePrefix(tower_prefix);
0124   se->registerSubsystem(rcemc);
0125
0126
0127   // do unsubtracted jet reconstruction
0128   TowerJetInput *incemc = new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER,tower_prefix);
0129   TowerJetInput *inihcal = new TowerJetInput(Jet::HCALIN_TOWERINFO, tower_prefix);
0130   TowerJetInput *inohcal = new TowerJetInput(Jet::HCALOUT_TOWERINFO,tower_prefix);
0131   incemc->set_GlobalVertexType(vertex_type);
0132   inihcal->set_GlobalVertexType(vertex_type);
0133   inohcal->set_GlobalVertexType(vertex_type);
0134
0135   JetReco *_jetRecoUnsub = new JetReco();
0136   _jetRecoUnsub->add_input(incemc);
0137   _jetRecoUnsub->add_input(inihcal);
0138   _jetRecoUnsub->add_input(inohcal);
0139   _jetRecoUnsub->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.2), "AntiKt_unsubtracted_r02");
0140   _jetRecoUnsub->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.3), "AntiKt_unsubtracted_r03");
0141   _jetRecoUnsub->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.4), "AntiKt_unsubtracted_r04");
0142   _jetRecoUnsub->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.5), "AntiKt_unsubtracted_r05");
0143   _jetRecoUnsub->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.6), "AntiKt_unsubtracted_r06");
0144   _jetRecoUnsub->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.7), "AntiKt_unsubtracted_r07");
0145   _jetRecoUnsub->add_algo(new FastJetAlgoSub(Jet::ANTIKT, 0.8), "AntiKt_unsubtracted_r08");
0146
0147
0148   _jetRecoUnsub->set_algo_node("ANTIKT");
0149   _jetRecoUnsub->set_input_node("TOWER");
0150   se->registerSubsystem(_jetRecoUnsub);
0151
0152   JetDSTSkimmer *jetDSTSkimmer = new JetDSTSkimmer();
0153   std::map<std::string, float> jetNodePts;
0154   jetNodePts["AntiKt_unsubtracted_r02"] = 6.6;
0155   jetNodePts["AntiKt_unsubtracted_r03"] = 7.3;
0156   jetNodePts["AntiKt_unsubtracted_r04"] = 7.6;
0157   jetNodePts["AntiKt_unsubtracted_r05"] = 7.9;
0158   jetNodePts["AntiKt_unsubtracted_r06"] = 8.0;
0159   jetNodePts["AntiKt_unsubtracted_r07"] = 11.0;
0160   jetNodePts["AntiKt_unsubtracted_r08"] = 12.0;
0161   jetDSTSkimmer->SetJetNodeThresholds(jetNodePts);
0162   std::map<std::string, float> clusterNodePts;
0163   clusterNodePts["CLUSTERINFO_CEMC"] = 5;
0164   jetDSTSkimmer->SetClusterNodeThresholds(clusterNodePts);
0165   se->registerSubsystem(jetDSTSkimmer);
0166
0167   // register modules necessary for QA
0168   if (Enable::QA)
0169   {
0170     DoRhoCalculation();
0171     Jet_QA();
0172   }
0173
0174   Fun4AllInputManager *In = new Fun4AllDstInputManager("in");
0175   In->AddFile(fname);
0176   se->registerInputManager(In);
0177
0178   Fun4AllDstOutputManager *outlower = new Fun4AllDstOutputManager("DSTOUTLOW", outfile_low);
0179   outlower->AddNode("Sync");
0180   outlower->AddNode("EventHeader");
0181   //outlower->AddNode("PacketsKeep");
0182   outlower->AddNode("14001");
0183   outlower->AddNode("1001");
0184   outlower->AddNode("1002");
0185   outlower->AddNode("TOWERS_HCALIN");
0186   outlower->AddNode("TOWERS_HCALOUT");
0187   outlower->AddNode("TOWERS_CEMC");
0188   outlower->AddNode("TOWERS_SEPD");
0189   outlower->AddNode("TOWERS_ZDC");
0190   //outlower->AddNode("MbdOut");
0191   //outlower->AddNode("MbdPmtContainer");
0192   //outlower->AddNode("MBDPackets");
0193   outlower->AddNode("TriggerRunInfo");
0194   se->registerOutputManager(outlower);
0195
0196   Fun4AllDstOutputManager *outhigher = new Fun4AllDstOutputManager("DSTOUTHIGH", outfile_high);
0197   outhigher->StripNode("1001");
0198   outhigher->StripNode("1002");
0199   outhigher->StripNode("TOWERINFO_CALIB_HCALIN");
0200   outhigher->StripNode("TOWERS_HCALIN");
0201   outhigher->StripNode("TOWERINFO_CALIB_HCALOUT");
0202   outhigher->StripNode("TOWERS_HCALOUT");
0203   outhigher->StripNode("TOWERINFO_CALIB_CEMC");
0204   outhigher->StripNode("TOWERS_CEMC");
0205   outhigher->StripNode("TOWERS_SEPD");
0206   outhigher->StripNode("TOWERINFO_CALIB_ZDC");
0207   outhigher->StripNode("TOWERS_ZDC");
0208   se->registerOutputManager(outhigher);
0209
0210   se->run(nEvents);
0211   se->End();
0212
0213   if (Enable::QA)
0214   {
0215     QAHistManagerDef::saveQARootFile(outfile_hist);
0216   }
0217
0218   CDBInterface::instance()->Print();  // print used DB files
0219   se->PrintTimer();
0220   delete se;
0221   std::cout << "All done!" << std::endl;
0222   gSystem->Exit(0);
0223 }
0224
0225 #endif