ForwardCalo/Flux_Return_Study/JetEnergies.C

0001
0002 #include "JetEnergies.h"
0003
0004 #include <cassert>
0005
0006 #include <g4eval/JetEvalStack.h>
0007 #include <g4eval/JetRecoEval.h>
0008
0009 #include <fun4all/Fun4AllReturnCodes.h>
0010 #include <phool/getClass.h>
0011 #include <fun4all/SubsysReco.h>
0012 #include <phool/PHCompositeNode.h>
0013 #include <g4jets/JetMap.h>
0014 #include <g4jets/Jet.h>
0015
0016 #include <g4main/PHG4HitContainer.h>
0017
0018 #include <TNtuple.h>
0019 #include <TFile.h>
0020
0021 #include <iostream>
0022 #include <cmath>
0023
0024 using namespace std;
0025
0026 JetEnergies::JetEnergies(const string &name,
0027                const string &recojetname,
0028                const string &truthjetname,
0029                const string &filename)
0030   : SubsysReco(name),
0031     _recojetname(recojetname),
0032     _truthjetname(truthjetname),
0033     _ievent(0),
0034     _jetevalstack(NULL),
0035     _strict(false),
0036     _errors(0),
0037     _do_recojet_eval(true),
0038     _do_truthjet_eval(true),
0039     _ntp_recojet(NULL),
0040     _ntp_truthjet(NULL),
0041     _filename(filename),
0042     _tfile(NULL),
0043     _FluxReturn_plus_hit_container(NULL),
0044     _FluxReturn_minus_hit_container(NULL),
0045     _BH_1_hit_container(NULL),
0046     _BH_Forward_hit_container(NULL),
0047     _BH_Negative_hit_container(NULL)
0048 {
0049   verbosity = 0;
0050 }
0051
0052 int JetEnergies::Init(PHCompositeNode *topNode) {
0053
0054   _ievent = 0;
0055
0056   _tfile = new TFile(_filename.c_str(), "RECREATE");
0057
0058   if (_do_recojet_eval) _ntp_recojet = new TNtuple("ntp_recojet","reco jet => max truth jet",
0059                            "event:id:ncomp:eta:phi:e:pt:"
0060                            "gid:gncomp:geta:gphi:ge:gpt:"
0061                            "efromtruth");
0062
0063   if (_do_truthjet_eval) _ntp_truthjet = new TNtuple("ntp_truthjet","truth jet => best reco jet",
0064                              "event:gid:gncomp:geta:gphi:ge:gpt:"
0065                              "id:ncomp:eta:phi:e:pt:"
0066                              "efromtruth:e_FR_plus:e_FR_minus:e_BH1:e_BH_plus:e_BH_minus");
0067
0068
0069   return Fun4AllReturnCodes::EVENT_OK;
0070 }
0071
0072 //Added 'InitRun(...) since the hit container for the Flux Return was not being intialized in 'Init(...){...}'
0073 int JetEnergies::InitRun(PHCompositeNode *topNode)
0074 {
0075   //cout << "In InitRun" << endl;
0076   //Here I add the container for the Plug Door hits
0077   _FluxReturn_plus_hit_container = findNode::getClass<PHG4HitContainer>(topNode,      "G4HIT_FLUXRET_ETA_PLUS");
0078   _FluxReturn_minus_hit_container = findNode::getClass<PHG4HitContainer>(topNode,      "G4HIT_FLUXRET_ETA_MINUS");
0079
0080   _BH_1_hit_container = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_BH_1");
0081   _BH_Forward_hit_container = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_BH_FORWARD_PLUS");
0082   _BH_Negative_hit_container = findNode::getClass<PHG4HitContainer>(topNode,"G4HIT_BH_FORWARD_NEG");
0083
0084
0085   //Used for debugging
0086   /*
0087   if( _FluxReturn_hit_container )
0088     {
0089       cout << "Correct name" << endl;
0090     }
0091   */
0092
0093   return Fun4AllReturnCodes::EVENT_OK;
0094
0095 }
0096
0097 int JetEnergies::process_event(PHCompositeNode *topNode) {
0098
0099   //cout << "In process Event" << endl;
0100
0101   if (!_jetevalstack) {
0102     _jetevalstack = new JetEvalStack(topNode,_recojetname,_truthjetname);
0103     _jetevalstack->set_strict(_strict);
0104     _jetevalstack->set_verbosity(verbosity+1);
0105   } else {
0106     _jetevalstack->next_event(topNode);
0107   }
0108
0109   //-----------------------------------
0110   // print what is coming into the code
0111   //-----------------------------------
0112
0113   printInputInfo(topNode);
0114
0115   //---------------------------
0116   // fill the Evaluator NTuples
0117   //---------------------------
0118
0119   fillOutputNtuples(topNode);
0120
0121   //--------------------------------------------------
0122   // Print out the ancestry information for this event
0123   //--------------------------------------------------
0124
0125   printOutputInfo(topNode);
0126
0127   ++_ievent;
0128
0129   return Fun4AllReturnCodes::EVENT_OK;
0130 }
0131
0132 int JetEnergies::End(PHCompositeNode *topNode) {
0133
0134   _tfile->cd();
0135
0136   if (_do_recojet_eval) _ntp_recojet->Write();
0137   if (_do_truthjet_eval) _ntp_truthjet->Write();
0138
0139   _tfile->Close();
0140
0141   delete _tfile;
0142
0143   if (verbosity > 0) {
0144     cout << "========================== JetEnergies::End() ============================" << endl;
0145     cout << " " << _ievent << " events of output written to: " << _filename << endl;
0146     cout << "===========================================================================" << endl;
0147   }
0148
0149   delete _jetevalstack;
0150
0151   return Fun4AllReturnCodes::EVENT_OK;
0152 }
0153
0154 void JetEnergies::printInputInfo(PHCompositeNode *topNode) {
0155   // to be implemented later if needed
0156   return;
0157 }
0158
0159 void JetEnergies::printOutputInfo(PHCompositeNode *topNode) {
0160   // to be implemented later if needed
0161   return;
0162 }
0163
0164 void JetEnergies::fillOutputNtuples(PHCompositeNode *topNode) {
0165
0166   if (verbosity > 2) cout << "JetEnergies::fillOutputNtuples() entered" << endl;
0167
0168   float e_FR_plus = GetTotalEnergy( _FluxReturn_plus_hit_container );
0169   float e_FR_minus = GetTotalEnergy( _FluxReturn_minus_hit_container );
0170   float e_BH1 = GetTotalEnergy( _BH_1_hit_container );
0171   float e_BH_plus = GetTotalEnergy( _BH_Forward_hit_container );
0172   float e_BH_minus = GetTotalEnergy( _BH_Negative_hit_container );
0173
0174   JetRecoEval*   recoeval = _jetevalstack->get_reco_eval();
0175   //JetTruthEval* trutheval = _jetevalstack->get_truth_eval();
0176
0177   //-------------------------
0178   // fill the reco jet ntuple
0179   //-------------------------
0180
0181   if (_do_recojet_eval) {
0182     if (verbosity > 1) cout << "JetEnergies::filling recojet ntuple..." << endl;
0183
0184     JetMap* recojets = findNode::getClass<JetMap>(topNode,_recojetname.c_str());
0185     if (!recojets) {
0186       cerr << PHWHERE << " ERROR: Can't find " << _recojetname << endl;
0187       exit(-1);
0188     }
0189
0190     // for every recojet
0191     for (JetMap::Iter iter = recojets->begin();
0192      iter != recojets->end();
0193      ++iter) {
0194       Jet* recojet = iter->second;
0195       Jet* truthjet = recoeval->max_truth_jet_by_energy(recojet);
0196
0197       float id    = recojet->get_id();
0198       float ncomp = recojet->size_comp();
0199       float eta   = recojet->get_eta();
0200       float phi   = recojet->get_phi();
0201       float e     = recojet->get_e();
0202       float pt    = recojet->get_pt();
0203
0204       float gid        = NAN;
0205       float gncomp     = NAN;
0206       float geta       = NAN;
0207       float gphi       = NAN;
0208       float ge         = NAN;
0209       float gpt        = NAN;
0210       float efromtruth = NAN;
0211
0212       if (truthjet) {
0213     gid    = truthjet->get_id();
0214     gncomp = truthjet->size_comp();
0215     geta   = truthjet->get_eta();
0216     gphi   = truthjet->get_phi();
0217     ge     = truthjet->get_e();
0218     gpt    = truthjet->get_pt();
0219     efromtruth = recoeval->get_energy_contribution(recojet,truthjet);
0220       }
0221
0222       float recojet_data[14] = {(float) _ievent,
0223                 id,
0224                 ncomp,
0225                 eta,
0226                 phi,
0227                 e,
0228                 pt,
0229                 gid,
0230                 gncomp,
0231                 geta,
0232                 gphi,
0233                 ge,
0234                 gpt,
0235                 efromtruth
0236       };
0237
0238       _ntp_recojet->Fill(recojet_data);
0239     }
0240   }
0241
0242   //-------------------------
0243   // fill the truth jet ntuple
0244   //-------------------------
0245
0246   if (_do_truthjet_eval) {
0247     if (verbosity > 1) cout << "JetEnergies::filling truthjet ntuple..." << endl;
0248
0249     JetMap* truthjets = findNode::getClass<JetMap>(topNode,_truthjetname.c_str());
0250     if (!truthjets) {
0251       cerr << PHWHERE << " ERROR: Can't find " << _truthjetname << endl;
0252       exit(-1);
0253     }
0254
0255     // for every truthjet
0256     for (JetMap::Iter iter = truthjets->begin();
0257      iter != truthjets->end();
0258      ++iter) {
0259       Jet* truthjet = iter->second;
0260       Jet* recojet = recoeval->best_jet_from(truthjet);
0261
0262       float gid    = truthjet->get_id();
0263       float gncomp = truthjet->size_comp();
0264       float geta   = truthjet->get_eta();
0265       float gphi   = truthjet->get_phi();
0266       float ge     = truthjet->get_e();
0267       float gpt    = truthjet->get_pt();
0268
0269       float id         = NAN;
0270       float ncomp      = NAN;
0271       float eta        = NAN;
0272       float phi        = NAN;
0273       float e          = NAN;
0274       float pt         = NAN;
0275       float efromtruth = NAN;
0276
0277       if (recojet) {
0278     id         = recojet->get_id();
0279     ncomp      = recojet->size_comp();
0280     eta        = recojet->get_eta();
0281     phi        = recojet->get_phi();
0282     e          = recojet->get_e();
0283     pt         = recojet->get_pt();
0284     efromtruth = recoeval->get_energy_contribution(recojet,truthjet);
0285       }
0286
0287       float truthjet_data[19] = {(float) _ievent,
0288                  gid,
0289                  gncomp,
0290                  geta,
0291                  gphi,
0292                  ge,
0293                  gpt,
0294                  id,
0295                  ncomp,
0296                  eta,
0297                  phi,
0298                  e,
0299                  pt,
0300                  efromtruth,
0301                  e_FR_plus,
0302                  e_FR_minus,
0303                  e_BH1,
0304                  e_BH_plus,
0305                  e_BH_minus
0306       };
0307
0308       _ntp_truthjet->Fill(truthjet_data);
0309     }
0310   }
0311
0312   return;
0313 }
0314
0315
0316 //Function that iterates over all the hits for a given object and returns the total energy deposited in that object requires that the object is set as active in the G4_Setup file otherwise may not work
0317 float JetEnergies::GetTotalEnergy(PHG4HitContainer *hit_object)
0318 {
0319   //cout << "Executing Modified part of code" << endl;
0320   float e_object = 0;//Holds Energy deposited in object of interest
0321   //Check if hit container exists mostly used for debugging
0322   if( hit_object )
0323     {
0324       //cout << "Execute my code" << endl;
0325       PHG4HitContainer::ConstRange object_hit_range = hit_object->getHits();
0326       //For loop that will store the total energy deposited in the object in 'e_object'
0327       //cout << "Starting for loop" << endl;
0328       for (PHG4HitContainer::ConstIterator hit_iter = object_hit_range.first; hit_iter != object_hit_range.second; hit_iter++)
0329     {
0330       PHG4Hit *this_hit = hit_iter->second;
0331       assert(this_hit);
0332       e_object += this_hit->get_edep();
0333     }
0334
0335       //cout << "End of for loop" << endl;
0336       //cout << "Value of e_FR: " << e_FR << endl;
0337
0338     }
0339   //cout << "Executing Rest of original code" << endl;
0340   return e_object;
0341
0342 }