myjetanalysis/src/MyJetAnalysis.cc

0001 #include "MyJetAnalysis.h"
0002
0003 #include <trackbase_historic/SvtxTrackMap.h>
0004
0005 #include <jetbase/JetMap.h>
0006
0007 #include <fun4all/Fun4AllReturnCodes.h>
0008 #include <fun4all/Fun4AllServer.h>
0009 #include <fun4all/PHTFileServer.h>
0010
0011 #include <phool/PHCompositeNode.h>
0012 #include <phool/getClass.h>
0013
0014 #include <TFile.h>
0015 #include <TH1F.h>
0016 #include <TH2F.h>
0017 #include <TString.h>
0018 #include <TTree.h>
0019 #include <TVector3.h>
0020
0021 #include <algorithm>
0022 #include <cassert>
0023 #include <cmath>
0024 #include <iostream>
0025 #include <limits>
0026 #include <stdexcept>
0027
0028 MyJetAnalysis::MyJetAnalysis(const std::string& recojetname, const std::string& truthjetname, const std::string& outputfilename)
0029   : SubsysReco("MyJetAnalysis_" + recojetname + "_" + truthjetname)
0030   , m_recoJetName(recojetname)
0031   , m_truthJetName(truthjetname)
0032   , m_outputFileName(outputfilename)
0033   , m_etaRange(-1, 1)
0034   , m_ptRange(5, 100)
0035   , m_trackJetMatchingRadius(.7)
0036 {
0037   m_trackdR.fill(std::numeric_limits<float>::signaling_NaN());
0038   m_trackpT.fill(std::numeric_limits<float>::signaling_NaN());
0039 }
0040
0041 int MyJetAnalysis::Init(PHCompositeNode* topNode)
0042 {
0043   if (Verbosity() >= MyJetAnalysis::VERBOSITY_SOME)
0044     std::cout << "MyJetAnalysis::Init - Outoput to " << m_outputFileName << std::endl;
0045
0046   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0047
0048   // Histograms
0049   m_hInclusiveE = new TH1F(
0050       "hInclusive_E",  //
0051       TString(m_recoJetName) + " inclusive jet E;Total jet energy (GeV)", 100, 0, 100);
0052
0053   m_hInclusiveEta =
0054       new TH1F(
0055           "hInclusive_eta",  //
0056           TString(m_recoJetName) + " inclusive jet #eta;#eta;Jet energy density", 50, -1, 1);
0057   m_hInclusivePhi =
0058       new TH1F(
0059           "hInclusive_phi",  //
0060           TString(m_recoJetName) + " inclusive jet #phi;#phi;Jet energy density", 50, -M_PI, M_PI);
0061
0062   //Trees
0063   m_T = new TTree("T", "MyJetAnalysis Tree");
0064
0065   //      int m_event;
0066   m_T->Branch("m_event", &m_event, "event/I");
0067   //      int m_id;
0068   m_T->Branch("id", &m_id, "id/I");
0069   //      int m_nComponent;
0070   m_T->Branch("nComponent", &m_nComponent, "nComponent/I");
0071   //      float m_eta;
0072   m_T->Branch("eta", &m_eta, "eta/F");
0073   //      float m_phi;
0074   m_T->Branch("phi", &m_phi, "phi/F");
0075   //      float m_e;
0076   m_T->Branch("e", &m_e, "e/F");
0077   //      float m_pt;
0078   m_T->Branch("pt", &m_pt, "pt/F");
0079   //
0080   //      int m_truthID;
0081   m_T->Branch("truthID", &m_truthID, "truthID/I");
0082   //      int m_truthNComponent;
0083   m_T->Branch("truthNComponent", &m_truthNComponent, "truthNComponent/I");
0084   //      float m_truthEta;
0085   m_T->Branch("truthEta", &m_truthEta, "truthEta/F");
0086   //      float m_truthPhi;
0087   m_T->Branch("truthPhi", &m_truthPhi, "truthPhi/F");
0088   //      float m_truthE;
0089   m_T->Branch("truthE", &m_truthE, "truthE/F");
0090   //      float m_truthPt;
0091   m_T->Branch("truthPt", &m_truthPt, "truthPt/F");
0092   //
0093   //      //! number of matched tracks
0094   //      int m_nMatchedTrack;
0095   m_T->Branch("nMatchedTrack", &m_nMatchedTrack, "nMatchedTrack/I");
0096   //      std::array<float, kMaxMatchedTrack> m_trackdR;
0097   m_T->Branch("id", m_trackdR.data(), "trackdR[nMatchedTrack]/F");
0098   //      std::array<float, kMaxMatchedTrack> m_trackpT;
0099   m_T->Branch("id", m_trackpT.data(), "trackpT[nMatchedTrack]/F");
0100
0101   return Fun4AllReturnCodes::EVENT_OK;
0102 }
0103
0104 int MyJetAnalysis::End(PHCompositeNode* topNode)
0105 {
0106   std::cout << "MyJetAnalysis::End - Output to " << m_outputFileName << std::endl;
0107   PHTFileServer::get().cd(m_outputFileName);
0108
0109   m_hInclusiveE->Write();
0110   m_hInclusiveEta->Write();
0111   m_hInclusivePhi->Write();
0112   m_T->Write();
0113
0114   return Fun4AllReturnCodes::EVENT_OK;
0115 }
0116
0117 int MyJetAnalysis::InitRun(PHCompositeNode* topNode)
0118 {
0119   m_jetEvalStack = std::shared_ptr<JetEvalStack>(new JetEvalStack(topNode, m_recoJetName, m_truthJetName));
0120   m_jetEvalStack->get_stvx_eval_stack()->set_use_initial_vertex(initial_vertex);
0121   return Fun4AllReturnCodes::EVENT_OK;
0122 }
0123
0124 int MyJetAnalysis::process_event(PHCompositeNode* topNode)
0125 {
0126   if (Verbosity() >= MyJetAnalysis::VERBOSITY_SOME)
0127     std::cout << "MyJetAnalysis::process_event() entered" << std::endl;
0128
0129   m_jetEvalStack->next_event(topNode);
0130   JetRecoEval* recoeval = m_jetEvalStack->get_reco_eval();
0131   ++m_event;
0132
0133   // interface to jets
0134   JetMap* jets = findNode::getClass<JetMap>(topNode, m_recoJetName);
0135   if (!jets)
0136   {
0137     std::cout
0138         << "MyJetAnalysis::process_event - Error can not find DST JetMap node "
0139         << m_recoJetName << std::endl;
0140     exit(-1);
0141   }
0142
0143   // interface to tracks
0144   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0145   if (!trackmap)
0146   {
0147     trackmap = findNode::getClass<SvtxTrackMap>(topNode, "TrackMap");
0148     if (!trackmap)
0149     {
0150       std::cout
0151           << "MyJetAnalysis::process_event - Error can not find DST trackmap node SvtxTrackMap" << std::endl;
0152       exit(-1);
0153     }
0154   }
0155   for (JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter)
0156   {
0157     Jet* jet = iter->second;
0158     assert(jet);
0159
0160     bool eta_cut = (jet->get_eta() >= m_etaRange.first) and (jet->get_eta() <= m_etaRange.second);
0161     bool pt_cut = (jet->get_pt() >= m_ptRange.first) and (jet->get_pt() <= m_ptRange.second);
0162     if ((not eta_cut) or (not pt_cut))
0163     {
0164       if (Verbosity() >= MyJetAnalysis::VERBOSITY_MORE)
0165       {
0166         std::cout << "MyJetAnalysis::process_event() - jet failed acceptance cut: ";
0167         std::cout << "eta cut: " << eta_cut << ", ptcut: " << pt_cut << std::endl;
0168         std::cout << "jet eta: " << jet->get_eta() << ", jet pt: " << jet->get_pt() << std::endl;
0169         jet->identify();
0170       }
0171       continue;
0172     }
0173
0174     // fill histograms
0175     assert(m_hInclusiveE);
0176     m_hInclusiveE->Fill(jet->get_e());
0177     assert(m_hInclusiveEta);
0178     m_hInclusiveEta->Fill(jet->get_eta());
0179     assert(m_hInclusivePhi);
0180     m_hInclusivePhi->Fill(jet->get_phi());
0181
0182     // fill trees - jet spectrum
0183     Jet* truthjet = recoeval->max_truth_jet_by_energy(jet);
0184
0185     m_id = jet->get_id();
0186     m_nComponent = jet->size_comp();
0187     m_eta = jet->get_eta();
0188     m_phi = jet->get_phi();
0189     m_e = jet->get_e();
0190     m_pt = jet->get_pt();
0191
0192     m_truthID = -1;
0193     m_truthNComponent = -1;
0194     m_truthEta = NAN;
0195     m_truthPhi = NAN;
0196     m_truthE = NAN;
0197     m_truthPt = NAN;
0198
0199     if (truthjet)
0200     {
0201       m_truthID = truthjet->get_id();
0202       m_truthNComponent = truthjet->size_comp();
0203       m_truthEta = truthjet->get_eta();
0204       m_truthPhi = truthjet->get_phi();
0205       m_truthE = truthjet->get_e();
0206       m_truthPt = truthjet->get_pt();
0207     }
0208
0209     // fill trees - jet track matching
0210     m_nMatchedTrack = 0;
0211
0212     for (SvtxTrackMap::Iter iter = trackmap->begin();
0213          iter != trackmap->end();
0214          ++iter)
0215     {
0216       SvtxTrack* track = iter->second;
0217
0218       TVector3 v(track->get_px(), track->get_py(), track->get_pz());
0219       const double dEta = v.Eta() - m_eta;
0220       const double dPhi = v.Phi() - m_phi;
0221       const double dR = sqrt(dEta * dEta + dPhi * dPhi);
0222
0223       if (dR < m_trackJetMatchingRadius)
0224       {
0225         //matched track to jet
0226
0227         assert(m_nMatchedTrack < kMaxMatchedTrack);
0228
0229         m_trackdR[m_nMatchedTrack] = dR;
0230         m_trackpT[m_nMatchedTrack] = v.Perp();
0231
0232         ++m_nMatchedTrack;
0233       }
0234
0235       if (m_nMatchedTrack >= kMaxMatchedTrack)
0236       {
0237         std::cout << "MyJetAnalysis::process_event() - reached max track that matching a jet. Quit iterating tracks" << std::endl;
0238         break;
0239       }
0240
0241     }  //    for (SvtxTrackMap::Iter iter = trackmap->begin();
0242
0243     m_T->Fill();
0244   }  //   for (JetMap::Iter iter = jets->begin(); iter != jets->end(); ++iter)
0245
0246   return Fun4AllReturnCodes::EVENT_OK;
0247 }