macros/common/HIJetReco.C

0001 #ifndef MACRO_HIJETRECO_C
0002 #define MACRO_HIJETRECO_C
0003
0004 #include <GlobalVariables.C>
0005
0006 #include <g4jets/TruthJetInput.h>
0007
0008 #include <globalvertex/GlobalVertex.h>
0009
0010 #include <jetbackground/CopyAndSubtractJets.h>
0011 #include <jetbackground/DetermineEventRho.h>
0012 #include <jetbackground/DetermineTowerBackground.h>
0013 #include <jetbackground/DetermineTowerRho.h>
0014 #include <jetbackground/FastJetAlgoSub.h>
0015 #include <jetbackground/RetowerCEMC.h>
0016 #include <jetbackground/SubtractTowers.h>
0017 #include <jetbackground/SubtractTowersCS.h>
0018 #include <jetbackground/TowerRho.h>
0019
0020 #include <jetbase/FastJetOptions.h>
0021 #include <jetbase/JetReco.h>
0022 #include <jetbase/TowerJetInput.h>
0023 #include <jetbase/TrackJetInput.h>
0024
0025 #include <particleflowreco/ParticleFlowJetInput.h>
0026
0027 #include <eventplaneinfo/EventPlaneReco.h>
0028 #include <eventplaneinfo/Eventplaneinfo.h>
0029
0030 #include <fun4all/Fun4AllServer.h>
0031
0032 R__LOAD_LIBRARY(libg4jets.so)
0033 R__LOAD_LIBRARY(libglobalvertex.so)
0034 R__LOAD_LIBRARY(libjetbackground.so)
0035 R__LOAD_LIBRARY(libjetbase.so)
0036 R__LOAD_LIBRARY(libparticleflow.so)
0037 R__LOAD_LIBRARY(libeventplaneinfo.so)
0038
0039 // ----------------------------------------------------------------------------
0040 //! General options for background-subtracted jet reconstruction
0041 // ----------------------------------------------------------------------------
0042 namespace Enable
0043 {
0044   bool HIJETS = false;        ///< do HI jet reconstruction
0045   int HIJETS_VERBOSITY = 0;   ///< verbosity
0046   bool HIJETS_MC = false;     ///< is simulation
0047   bool HIJETS_TRUTH = false;  ///< make truth jets
0048   bool HIJETS_TOWER = true;   ///< make tower jets
0049   bool HIJETS_TRACK = false;  ///< make track jets
0050   bool HIJETS_PFLOW = false;  ///< make particle flow jets
0051 }  // namespace Enable
0052
0053 // ----------------------------------------------------------------------------
0054 //! Options specific to background-subtracted jet reconstruction
0055 // ----------------------------------------------------------------------------
0056 namespace HIJETS
0057 {
0058   // do_flow = 0 --noflow
0059   // do_flow = 1 --psi2 derived from calo
0060   // do_flow = 2 --psi2 derived from HIJING
0061   // do_flow = 3 --psi2 derived from sEPD
0062   int do_flow = 0;
0063
0064   ///! do constituent subtraction
0065   bool do_CS = false;
0066
0067   ///! turn on/off functionality only relevant for
0068   ///! nucleon collisions
0069   bool is_pp = true;
0070
0071   ///! sets prefix of nodes to use as tower jet
0072   ///! input
0073   std::string tower_prefix = "TOWERINFO_CALIB";
0074
0075   ///! turn on/off vertex specification
0076   bool do_vertex_type = true;
0077
0078   ///! if specifying vertex, set vertex
0079   ///! to use to this one
0080   GlobalVertex::VTXTYPE vertex_type = GlobalVertex::MBD;
0081
0082   ///! Base fastjet options to use. Note that the
0083   ///! resolution parameter will be overwritten
0084   ///! to R = 0.2, 0.3, 0.4, and 0.5
0085   FastJetOptions fj_opts({Jet::ANTIKT, JET_R, 0.4, VERBOSITY, static_cast<float>(Enable::HIJETS_VERBOSITY)});
0086
0087   ///! sets jet node name
0088   std::string jet_node = "ANTIKT";
0089
0090   ///! sets prefix of nodes to store jets
0091   std::string algo_prefix = "AntiKt";
0092
0093   ///! sets the embedding g4 flag used to determine which particles are classified as 'truth' in MakeHITruthJets
0094   ///! if you don't know what this is, you probably don't need to change it
0095   ///! 0 = all particles (use for all g4particles !not useful for HIJING+Pythia samples)
0096   ///! 1 = pythia/herwig particles only (use for pythia/herwig jets in HIJING+Pythia samples)
0097   ///! 2 = pythia particles from the HIJING+Pythia samples (use for pythia jets in HIJING+Pythia samples)
0098   ///! negative values are typically background particles (HIJING particles)
0099   int embedding_flag = 1;
0100
0101   ///! enumerates reconstructed resolution
0102   ///! parameters
0103   enum Res
0104   {
0105     R02 = 0,
0106     R03 = 1,
0107     R04 = 2,
0108     R05 = 3
0109   };
0110
0111   // --------------------------------------------------------------------------
0112   //! Helper method to generate releveant FastJet algorithms
0113   // --------------------------------------------------------------------------
0114   FastJetAlgoSub *GetFJAlgo(const float reso)
0115   {
0116     // grab current options & update
0117     // reso parameter
0118     FastJetOptions opts = fj_opts;
0119     opts.jet_R = reso;
0120
0121     // create new algorithm
0122     return new FastJetAlgoSub(opts);
0123
0124   }  // end 'GetFJAlgo()'
0125 }  // namespace HIJETS
0126
0127 // ----------------------------------------------------------------------------
0128 //! Make jets out of appropriate truth particles
0129 // ----------------------------------------------------------------------------
0130 void MakeHITruthJets()
0131 {
0132   // set verbosity
0133   int verbosity = std::max(Enable::VERBOSITY, Enable::HIJETS_VERBOSITY);
0134
0135   //---------------
0136   // Fun4All server
0137   //---------------
0138   Fun4AllServer *se = Fun4AllServer::instance();
0139
0140   // if making track jets, make truth jets out of only charged particles
0141   if (Enable::HIJETS_TRACK)
0142   {
0143     // configure truth jet input for charged particles
0144     TruthJetInput *ctji = new TruthJetInput(Jet::SRC::CHARGED_PARTICLE);
0145     ctji->add_embedding_flag(HIJETS::embedding_flag);  // changes depending on signal vs. embedded
0146
0147     // book jet reconstruction on chargedparticles
0148     JetReco *chargedtruthjetreco = new JetReco();
0149     chargedtruthjetreco->add_input(ctji);
0150     chargedtruthjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_ChargedTruth_r02");
0151     chargedtruthjetreco->add_algo(HIJETS::GetFJAlgo(0.3), HIJETS::algo_prefix + "_ChargedTruth_r03");
0152     chargedtruthjetreco->add_algo(HIJETS::GetFJAlgo(0.4), HIJETS::algo_prefix + "_ChargedTruth_r04");
0153     chargedtruthjetreco->add_algo(HIJETS::GetFJAlgo(0.5), HIJETS::algo_prefix + "_ChargedTruth_r05");
0154     chargedtruthjetreco->set_algo_node(HIJETS::jet_node);
0155     chargedtruthjetreco->set_input_node("TRUTH");
0156     chargedtruthjetreco->Verbosity(verbosity);
0157     se->registerSubsystem(chargedtruthjetreco);
0158   }
0159
0160   // if making tower or pflow jets, make truth jets out of all particles
0161   if (Enable::HIJETS_TOWER || Enable::HIJETS_PFLOW)
0162   {
0163     // configure truth jet input for all particles
0164     TruthJetInput *tji = new TruthJetInput(Jet::PARTICLE);
0165     tji->add_embedding_flag(HIJETS::embedding_flag);  // changes depending on signal vs. embedded
0166
0167     // book jet reconstruction on all particles
0168     JetReco *truthjetreco = new JetReco();
0169     truthjetreco->add_input(tji);
0170     truthjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_Truth_r02");
0171     truthjetreco->add_algo(HIJETS::GetFJAlgo(0.3), HIJETS::algo_prefix + "_Truth_r03");
0172     truthjetreco->add_algo(HIJETS::GetFJAlgo(0.4), HIJETS::algo_prefix + "_Truth_r04");
0173     truthjetreco->add_algo(HIJETS::GetFJAlgo(0.5), HIJETS::algo_prefix + "_Truth_r05");
0174     truthjetreco->set_algo_node(HIJETS::jet_node);
0175     truthjetreco->set_input_node("TRUTH");
0176     truthjetreco->Verbosity(verbosity);
0177     se->registerSubsystem(truthjetreco);
0178   }
0179
0180   // exit back to HIJetReco()
0181   return;
0182 }
0183
0184 // ----------------------------------------------------------------------------
0185 //! Make jets out of subtracted towers
0186 // ----------------------------------------------------------------------------
0187 void MakeHITowerJets()
0188 {
0189   int verbosity = std::max(Enable::VERBOSITY, Enable::HIJETS_VERBOSITY);
0190
0191   //---------------
0192   // Fun4All server
0193   //---------------
0194   Fun4AllServer *se = Fun4AllServer::instance();
0195
0196   if (HIJETS::do_flow == 3)
0197   {
0198     EventPlaneReco *epreco = new EventPlaneReco();
0199     epreco->set_sepd_epreco(true);
0200     se->registerSubsystem(epreco);
0201   }
0202
0203   RetowerCEMC *rcemc = new RetowerCEMC();
0204   rcemc->Verbosity(verbosity);
0205   rcemc->set_towerinfo(true);
0206   rcemc->set_frac_cut(0.5);  // fraction of retower that must be masked to mask the full retower
0207   rcemc->set_towerNodePrefix(HIJETS::tower_prefix);
0208   se->registerSubsystem(rcemc);
0209
0210   JetReco *towerjetreco = new JetReco();
0211   TowerJetInput *incemc = new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER, HIJETS::tower_prefix);
0212   TowerJetInput *inihcal = new TowerJetInput(Jet::HCALIN_TOWERINFO, HIJETS::tower_prefix);
0213   TowerJetInput *inohcal = new TowerJetInput(Jet::HCALOUT_TOWERINFO, HIJETS::tower_prefix);
0214   if (HIJETS::do_vertex_type)
0215   {
0216     incemc->set_GlobalVertexType(HIJETS::vertex_type);
0217     inihcal->set_GlobalVertexType(HIJETS::vertex_type);
0218     inohcal->set_GlobalVertexType(HIJETS::vertex_type);
0219   }
0220   towerjetreco->add_input(incemc);
0221   towerjetreco->add_input(inihcal);
0222   towerjetreco->add_input(inohcal);
0223   towerjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_TowerInfo_HIRecoSeedsRaw_r02");
0224   towerjetreco->set_algo_node(HIJETS::jet_node);
0225   towerjetreco->set_input_node("TOWER");
0226   towerjetreco->Verbosity(verbosity);
0227   se->registerSubsystem(towerjetreco);
0228
0229   DetermineTowerBackground *dtb = new DetermineTowerBackground();
0230   dtb->SetBackgroundOutputName("TowerInfoBackground_Sub1");
0231   dtb->SetFlow(HIJETS::do_flow);
0232   dtb->SetSeedType(0);
0233   dtb->SetSeedJetD(3);
0234   // dtb->set_towerinfo(true);
0235   dtb->Verbosity(verbosity);
0236   dtb->set_towerNodePrefix(HIJETS::tower_prefix);
0237   se->registerSubsystem(dtb);
0238
0239   CopyAndSubtractJets *casj = new CopyAndSubtractJets();
0240   casj->SetFlowModulation(HIJETS::do_flow);
0241   casj->Verbosity(verbosity);
0242   casj->set_towerinfo(true);
0243   casj->set_towerNodePrefix(HIJETS::tower_prefix);
0244   se->registerSubsystem(casj);
0245
0246   DetermineTowerBackground *dtb2 = new DetermineTowerBackground();
0247   dtb2->SetBackgroundOutputName("TowerInfoBackground_Sub2");
0248   dtb2->SetFlow(HIJETS::do_flow);
0249   dtb2->SetSeedType(1);
0250   dtb2->SetSeedJetPt(7);
0251   dtb2->Verbosity(verbosity);
0252   // dtb2->set_towerinfo(true);
0253   dtb2->set_towerNodePrefix(HIJETS::tower_prefix);
0254   se->registerSubsystem(dtb2);
0255
0256   SubtractTowers *st = new SubtractTowers();
0257   st->SetFlowModulation(HIJETS::do_flow);
0258   st->Verbosity(verbosity);
0259   st->set_towerinfo(true);
0260   st->set_towerNodePrefix(HIJETS::tower_prefix);
0261   se->registerSubsystem(st);
0262
0263   towerjetreco = new JetReco();
0264   incemc = new TowerJetInput(Jet::CEMC_TOWERINFO_SUB1, HIJETS::tower_prefix);
0265   inihcal = new TowerJetInput(Jet::HCALIN_TOWERINFO_SUB1, HIJETS::tower_prefix);
0266   inohcal = new TowerJetInput(Jet::HCALOUT_TOWERINFO_SUB1, HIJETS::tower_prefix);
0267   if (HIJETS::do_vertex_type)
0268   {
0269     incemc->set_GlobalVertexType(HIJETS::vertex_type);
0270     inihcal->set_GlobalVertexType(HIJETS::vertex_type);
0271     inohcal->set_GlobalVertexType(HIJETS::vertex_type);
0272   }
0273   towerjetreco->add_input(incemc);
0274   towerjetreco->add_input(inihcal);
0275   towerjetreco->add_input(inohcal);
0276   towerjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_Tower_r02_Sub1");
0277   towerjetreco->add_algo(HIJETS::GetFJAlgo(0.3), HIJETS::algo_prefix + "_Tower_r03_Sub1");
0278   towerjetreco->add_algo(HIJETS::GetFJAlgo(0.4), HIJETS::algo_prefix + "_Tower_r04_Sub1");
0279   towerjetreco->add_algo(HIJETS::GetFJAlgo(0.5), HIJETS::algo_prefix + "_Tower_r05_Sub1");
0280   towerjetreco->set_algo_node(HIJETS::jet_node);
0281   towerjetreco->set_input_node("TOWER");
0282   towerjetreco->Verbosity(verbosity);
0283   se->registerSubsystem(towerjetreco);
0284
0285   return;
0286 }
0287
0288 // ----------------------------------------------------------------------------
0289 //! Make jets out of tracks with background subtraction
0290 // ----------------------------------------------------------------------------
0291 void MakeHITrackJets()
0292 {
0293   // set verbosity
0294   int verbosity = std::max(Enable::VERBOSITY, Enable::HIJETS_VERBOSITY);
0295
0296   //---------------
0297   // Fun4All server
0298   //---------------
0299   Fun4AllServer *se = Fun4AllServer::instance();
0300
0301   // emit warning: background sub will be added later
0302   std::cerr << "WARNING: Background subtraction for track jets is still in development!\n"
0303             << "  If you want to do jet reco without background subtraction, please\n"
0304             << "  use NoBkgdSubJetReco()"
0305             << std::endl;
0306
0307   // book jet reconstruction routines on tracks
0308   JetReco *trackjetreco = new JetReco();
0309   trackjetreco->add_input(new TrackJetInput(Jet::SRC::TRACK));
0310   trackjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_Track_r02");
0311   trackjetreco->add_algo(HIJETS::GetFJAlgo(0.3), HIJETS::algo_prefix + "_Track_r03");
0312   trackjetreco->add_algo(HIJETS::GetFJAlgo(0.4), HIJETS::algo_prefix + "_Track_r04");
0313   trackjetreco->add_algo(HIJETS::GetFJAlgo(0.5), HIJETS::algo_prefix + "_Track_r05");
0314   trackjetreco->set_algo_node(HIJETS::jet_node);
0315   trackjetreco->set_input_node("TRACK");
0316   trackjetreco->Verbosity(verbosity);
0317   se->registerSubsystem(trackjetreco);
0318
0319   // exit back to HIJetReco()
0320   return;
0321 }
0322
0323 // ----------------------------------------------------------------------------
0324 //! Make jets out of particle-flow elements with background subtraction
0325 // ----------------------------------------------------------------------------
0326 void MakeHIPFlowJets()
0327 {
0328   // set verbosity
0329   int verbosity = std::max(Enable::VERBOSITY, Enable::HIJETS_VERBOSITY);
0330
0331   //---------------
0332   // Fun4All server
0333   //---------------
0334   Fun4AllServer *se = Fun4AllServer::instance();
0335
0336   // emit warning: background sub will be added later
0337   std::cerr << "WARNING: Background subtraction for particle-flow jets is still in development!\n"
0338             << "  If you want to do jet reco without background subtraction, please\n"
0339             << "  use NoBkgdSubJetReco.C"
0340             << std::endl;
0341
0342   // book jet reconstruction routines on pflow elements
0343   JetReco *pflowjetreco = new JetReco();
0344   pflowjetreco->add_input(new ParticleFlowJetInput());
0345   pflowjetreco->add_algo(HIJETS::GetFJAlgo(0.2), HIJETS::algo_prefix + "_ParticleFlow_r02");
0346   pflowjetreco->add_algo(HIJETS::GetFJAlgo(0.3), HIJETS::algo_prefix + "_ParticleFlow_r03");
0347   pflowjetreco->add_algo(HIJETS::GetFJAlgo(0.4), HIJETS::algo_prefix + "_ParticleFlow_r04");
0348   pflowjetreco->add_algo(HIJETS::GetFJAlgo(0.5), HIJETS::algo_prefix + "_ParticleFlow_r05");
0349   pflowjetreco->set_algo_node(HIJETS::jet_node);
0350   pflowjetreco->set_input_node("ELEMENT");
0351   pflowjetreco->Verbosity(verbosity);
0352   se->registerSubsystem(pflowjetreco);
0353
0354   // exit back to HIJetReco()
0355   return;
0356 }
0357
0358 // ----------------------------------------------------------------------------
0359 //! Run background-subtracted jet reconstruction
0360 // ----------------------------------------------------------------------------
0361 void HIJetReco()
0362 {
0363   // if simulation, make appropriate truth jets
0364   if (Enable::HIJETS_MC && Enable::HIJETS_TRUTH) MakeHITruthJets();
0365
0366   // run approriate jet reconstruction routines
0367   if (Enable::HIJETS_TOWER) MakeHITowerJets();
0368   if (Enable::HIJETS_TRACK) MakeHITrackJets();
0369   if (Enable::HIJETS_PFLOW) MakeHIPFlowJets();
0370 }
0371
0372 // ----------------------------------------------------------------------------
0373 //! Determine rho from tower input to jet reco (necessary for jet QA)
0374 // ----------------------------------------------------------------------------
0375 void DoRhoCalculation()
0376 {
0377   // set verbosity
0378   //  int verbosity = std::max(Enable::VERBOSITY, Enable::HIJETS_VERBOSITY);
0379
0380   //---------------
0381   // Fun4All server
0382   //---------------
0383   Fun4AllServer *se = Fun4AllServer::instance();
0384
0385   // run rho calculations w/ towers
0386   if (Enable::HIJETS_TOWER)
0387   {
0388     DetermineTowerRho *towRhoCalc = new DetermineTowerRho();
0389     towRhoCalc->add_method(TowerRho::Method::AREA);
0390     towRhoCalc->add_method(TowerRho::Method::MULT);
0391     TowerJetInput *rho_incemc = new TowerJetInput(Jet::CEMC_TOWERINFO_RETOWER, HIJETS::tower_prefix);
0392     TowerJetInput *rho_inihcal = new TowerJetInput(Jet::HCALIN_TOWERINFO, HIJETS::tower_prefix);
0393     TowerJetInput *rho_inohcal = new TowerJetInput(Jet::HCALOUT_TOWERINFO, HIJETS::tower_prefix);
0394     if (HIJETS::do_vertex_type)
0395     {
0396       rho_incemc->set_GlobalVertexType(HIJETS::vertex_type);
0397       rho_inihcal->set_GlobalVertexType(HIJETS::vertex_type);
0398       rho_inohcal->set_GlobalVertexType(HIJETS::vertex_type);
0399     }
0400     towRhoCalc->add_tower_input(rho_incemc);
0401     towRhoCalc->add_tower_input(rho_inihcal);
0402     towRhoCalc->add_tower_input(rho_inohcal);
0403     se->registerSubsystem(towRhoCalc);
0404   }
0405
0406   // run rho calculations w/ tracks
0407   if (Enable::HIJETS_TRACK)
0408   {
0409     DetermineEventRho *trkRhoCalc = new DetermineEventRho();
0410     trkRhoCalc->add_method(EventRho::Method::AREA, "EventRho_AREA");
0411     trkRhoCalc->add_method(EventRho::Method::MULT, "EventRho_MULT");
0412     trkRhoCalc->add_input(new TrackJetInput(Jet::SRC::TRACK));
0413     se->registerSubsystem(trkRhoCalc);
0414   }
0415
0416   // exit back to main macro
0417   return;
0418 }
0419
0420 #endif