g4simulation/g4eval/JetEvaluator.cc

0001 #include "JetEvaluator.h"
0002
0003 #include "JetEvalStack.h"
0004 #include "JetRecoEval.h"
0005
0006 #include <jetbase/Jet.h>
0007 #include <jetbase/JetContainer.h>
0008
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010 #include <fun4all/SubsysReco.h>
0011
0012 #include <phool/getClass.h>
0013 #include <phool/phool.h>
0014
0015 #include <TFile.h>
0016 #include <TNtuple.h>
0017
0018 #include <cmath>
0019 #include <cstdlib>
0020 #include <iostream>
0021 #include <map>
0022 #include <utility>
0023
0024 JetEvaluator::JetEvaluator(const std::string &name,
0025                            const std::string &recojetname,
0026                            const std::string &truthjetname,
0027                            const std::string &filename)
0028   : SubsysReco(name)
0029   , _recojetname(recojetname)
0030   , _truthjetname(truthjetname)
0031   , _filename(filename)
0032 { }
0033
0034 int JetEvaluator::Init(PHCompositeNode * /*topNode*/)
0035 {
0036   _ievent = 0;
0037
0038   _tfile = new TFile(_filename.c_str(), "RECREATE");
0039
0040   if (_do_recojet_eval)
0041   {
0042     _ntp_recojet = new TNtuple("ntp_recojet", "reco jet => max truth jet",
0043                                "event:id:ncomp:eta:phi:e:pt:"
0044                                "gid:gncomp:geta:gphi:ge:gpt:"
0045                                "efromtruth");
0046   }
0047
0048   if (_do_truthjet_eval)
0049   {
0050     _ntp_truthjet = new TNtuple("ntp_truthjet", "truth jet => best reco jet",
0051                                 "event:gid:gncomp:geta:gphi:ge:gpt:"
0052                                 "id:ncomp:eta:phi:e:pt:"
0053                                 "efromtruth");
0054   }
0055
0056   return Fun4AllReturnCodes::EVENT_OK;
0057 }
0058
0059 int JetEvaluator::process_event(PHCompositeNode *topNode)
0060 {
0061   if (!_jetevalstack)
0062   {
0063     _jetevalstack = new JetEvalStack(topNode, _recojetname, _truthjetname);
0064     _jetevalstack->set_strict(_strict);
0065     _jetevalstack->set_verbosity(Verbosity() + 1);
0066   }
0067   else
0068   {
0069     _jetevalstack->next_event(topNode);
0070   }
0071
0072   //-----------------------------------
0073   // print what is coming into the code
0074   //-----------------------------------
0075
0076   printInputInfo(topNode);
0077
0078   //---------------------------
0079   // fill the Evaluator NTuples
0080   //---------------------------
0081
0082   fillOutputNtuples(topNode);
0083
0084   //--------------------------------------------------
0085   // Print out the ancestry information for this event
0086   //--------------------------------------------------
0087
0088   printOutputInfo(topNode);
0089
0090   ++_ievent;
0091
0092   return Fun4AllReturnCodes::EVENT_OK;
0093 }
0094
0095 int JetEvaluator::End(PHCompositeNode * /*topNode*/)
0096 {
0097   _tfile->cd();
0098
0099   if (_do_recojet_eval)
0100   {
0101     _ntp_recojet->Write();
0102   }
0103   if (_do_truthjet_eval)
0104   {
0105     _ntp_truthjet->Write();
0106   }
0107
0108   _tfile->Close();
0109
0110   delete _tfile;
0111
0112   if (Verbosity() > 0)
0113   {
0114     std::cout << "========================== JetEvaluator::End() ============================" << std::endl;
0115     std::cout << " " << _ievent << " events of output written to: " << _filename << std::endl;
0116     std::cout << "===========================================================================" << std::endl;
0117   }
0118
0119   delete _jetevalstack;
0120
0121   return Fun4AllReturnCodes::EVENT_OK;
0122 }
0123
0124 void JetEvaluator::printInputInfo(PHCompositeNode * /*topNode*/)
0125 {
0126   // to be implemented later if needed
0127   return;
0128 }
0129
0130 void JetEvaluator::printOutputInfo(PHCompositeNode * /*topNode*/)
0131 {
0132   // to be implemented later if needed
0133   return;
0134 }
0135
0136 void JetEvaluator::fillOutputNtuples(PHCompositeNode *topNode)
0137 {
0138   if (Verbosity() > 2)
0139   {
0140     std::cout << "JetEvaluator::fillOutputNtuples() entered" << std::endl;
0141   }
0142
0143
0144   JetRecoEval *recoeval = _jetevalstack->get_reco_eval();
0145   // JetTruthEval* trutheval = _jetevalstack->get_truth_eval();
0146
0147   //-------------------------
0148   // fill the reco jet ntuple
0149   //-------------------------
0150
0151   if (_do_recojet_eval)
0152   {
0153     if (Verbosity() > 1)
0154     {
0155       std::cout << "JetEvaluator::filling recojet ntuple..." << std::endl;
0156     }
0157
0158     JetContainer *recojets = findNode::getClass<JetContainer>(topNode, _recojetname);
0159     if (!recojets)
0160     {
0161       std::cout << PHWHERE << " ERROR: Can't find " << _recojetname << std::endl;
0162       exit(-1);
0163     }
0164
0165     // for every recojet
0166     for (auto recojet : *recojets)
0167     {
0168       /* Jet *recojet = iter.second; */
0169       Jet *truthjet = recoeval->max_truth_jet_by_energy(recojet);
0170
0171       float id = recojet->get_id();
0172       float ncomp = recojet->size_comp();
0173       float eta = recojet->get_eta();
0174       float phi = recojet->get_phi();
0175       float e = recojet->get_e();
0176       float pt = recojet->get_pt();
0177
0178       float gid = NAN;
0179       float gncomp = NAN;
0180       float geta = NAN;
0181       float gphi = NAN;
0182       float ge = NAN;
0183       float gpt = NAN;
0184       float efromtruth = NAN;
0185
0186       if (truthjet)
0187       {
0188         gid = truthjet->get_id();
0189         gncomp = truthjet->size_comp();
0190         geta = truthjet->get_eta();
0191         gphi = truthjet->get_phi();
0192         ge = truthjet->get_e();
0193         gpt = truthjet->get_pt();
0194         efromtruth = recoeval->get_energy_contribution(recojet, truthjet);
0195       }
0196
0197       float recojet_data[14] = {(float) _ievent,
0198                                 id,
0199                                 ncomp,
0200                                 eta,
0201                                 phi,
0202                                 e,
0203                                 pt,
0204                                 gid,
0205                                 gncomp,
0206                                 geta,
0207                                 gphi,
0208                                 ge,
0209                                 gpt,
0210                                 efromtruth};
0211
0212       _ntp_recojet->Fill(recojet_data);
0213     }
0214   }
0215
0216   //-------------------------
0217   // fill the truth jet ntuple
0218   //-------------------------
0219
0220   if (_do_truthjet_eval)
0221   {
0222     if (Verbosity() > 1)
0223     {
0224       std::cout << "JetEvaluator::filling truthjet ntuple..." << std::endl;
0225     }
0226
0227     JetContainer *truthjets = findNode::getClass<JetContainer>(topNode, _truthjetname);
0228     if (!truthjets)
0229     {
0230       std::cout << PHWHERE << " ERROR: Can't find " << _truthjetname << std::endl;
0231       exit(-1);
0232     }
0233
0234     // for every truthjet
0235     for (auto truthjet : *truthjets)
0236     {
0237       /* Jet *truthjet = iter.second; */
0238       Jet *recojet = recoeval->best_jet_from(truthjet);
0239
0240       float gid = truthjet->get_id();
0241       float gncomp = truthjet->size_comp();
0242       float geta = truthjet->get_eta();
0243       float gphi = truthjet->get_phi();
0244       float ge = truthjet->get_e();
0245       float gpt = truthjet->get_pt();
0246
0247       float id = NAN;
0248       float ncomp = NAN;
0249       float eta = NAN;
0250       float phi = NAN;
0251       float e = NAN;
0252       float pt = NAN;
0253       float efromtruth = NAN;
0254
0255       if (recojet)
0256       {
0257         id = recojet->get_id();
0258         ncomp = recojet->size_comp();
0259         eta = recojet->get_eta();
0260         phi = recojet->get_phi();
0261         e = recojet->get_e();
0262         pt = recojet->get_pt();
0263         efromtruth = recoeval->get_energy_contribution(recojet, truthjet);
0264       }
0265
0266       float truthjet_data[14] = {(float) _ievent,
0267                                  gid,
0268                                  gncomp,
0269                                  geta,
0270                                  gphi,
0271                                  ge,
0272                                  gpt,
0273                                  id,
0274                                  ncomp,
0275                                  eta,
0276                                  phi,
0277                                  e,
0278                                  pt,
0279                                  efromtruth};
0280
0281       _ntp_truthjet->Fill(truthjet_data);
0282     }
0283   }
0284
0285   return;
0286 }