macros/common/HIJetReco.C

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