TracksInJets/src/TracksInJets.cc

0001
0002
0003 #include "TracksInJets.h"
0004
0005 #include <trackbase_historic/SvtxTrack_v4.h>
0006 #include <trackbase_historic/SvtxTrackMap.h>
0007
0008 #include <trackbase_historic/SvtxTrackMap_v2.h>
0009
0010 #include <globalvertex/SvtxVertexMap.h>
0011 #include <globalvertex/GlobalVertex.h>
0012 #include <globalvertex/GlobalVertexMap.h>
0013 #include <globalvertex/GlobalVertexMapv1.h>
0014
0015
0016 #include <g4eval/SvtxEvalStack.h>
0017 #include <sstream>
0018
0019 #include <jetbase/JetMap.h>
0020 #include <jetbase/JetContainer.h>
0021 #include <jetbase/Jetv2.h>
0022 #include <jetbase/Jetv1.h>
0023 #include <jetbase/Jet.h>
0024 #include <jetbase/JetAlgo.h>
0025 #include <jetbase/FastJetAlgo.h>
0026 #include <jetbase/JetInput.h>
0027 #include <jetbase/TowerJetInput.h>
0028 #include <jetbase/JetMap.h>
0029 #include <jetbase/JetMapv1.h>
0030 #include <jetbase/JetContainer.h>
0031 #include <jetbase/JetContainerv1.h>
0032 #include <centrality/CentralityInfo.h>
0033
0034 #include <fun4all/Fun4AllBase.h>
0035 #include <fun4all/Fun4AllReturnCodes.h>
0036 #include <fun4all/PHTFileServer.h>
0037
0038 #include <phool/PHCompositeNode.h>
0039 #include <phool/getClass.h>
0040
0041
0042 #include "TH1.h"
0043 #include "TH2.h"
0044 #include "TH2F.h"
0045 #include <TTree.h>
0046 #include <TFile.h>
0047 #include <TH3F.h>
0048 #include <TH2F.h>
0049 #include <TString.h>
0050 #include <TVector3.h>
0051 #include <TTree.h>
0052 #include <cassert>
0053
0054 //____________________________________________________________________________..
0055 TracksInJets::TracksInJets(const std::string& recojetname, const std::string& outputfilename):
0056   SubsysReco("TracksInJets_" + recojetname)
0057   , m_recoJetName(recojetname)
0058   , m_trk_pt_cut(2)
0059   , m_jetRadius(0.4)
0060   , m_outputFileName(outputfilename)
0061 {
0062   std::cout << "TracksInJets::TracksInJets(const std::string &name) Calling ctor" << std::endl;
0063 }
0064
0065 //____________________________________________________________________________..
0066 TracksInJets::~TracksInJets()
0067 {
0068   std::cout << "TracksInJets::~TracksInJets() Calling dtor" << std::endl;
0069 }
0070
0071 //____________________________________________________________________________..
0072 int TracksInJets::Init(PHCompositeNode *topNode)
0073 {
0074   std::cout << "TracksInJets::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0075   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0076   m_h_track_vs_calo_pt = new TH3F("m_h_track_vs_calo_pt","",100,0,100,500,0,100,10,0,100);
0077   m_h_track_vs_calo_pt->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0078   m_h_track_vs_calo_pt->GetYaxis()->SetTitle("Sum track p_{T} [GeV]");
0079   m_h_track_pt = new TH2F("m_h_track_pt", "", 100, 0, 100, 100, 0, 100);
0080   m_h_track_pt->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0081   m_h_track_pt->GetYaxis()->SetTitle("Sum track p_{T} [GeV]");
0082   return Fun4AllReturnCodes::EVENT_OK;
0083 }
0084
0085 //____________________________________________________________________________..
0086 int TracksInJets::InitRun(PHCompositeNode *topNode)
0087 {
0088   std::cout << "TracksInJets::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0089   return Fun4AllReturnCodes::EVENT_OK;
0090 }
0091
0092 //____________________________________________________________________________..
0093 int TracksInJets::process_event(PHCompositeNode *topNode)
0094 {
0095   //std::cout << "TracksInJets::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0096  // interface to reco jets
0097   JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName);
0098   if (!jets)
0099     {
0100       std::cout
0101     << "MyJetAnalysis::process_event - Error can not find DST Reco JetContainer node "
0102     << m_recoJetName << std::endl;
0103       exit(-1);
0104     }
0105   //get reco tracks
0106   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0107   if (!trackmap)
0108     {
0109       trackmap = findNode::getClass<SvtxTrackMap>(topNode, "TrackMap");
0110       if (!trackmap)
0111     {
0112       std::cout
0113         << "TracksInJets::process_event - Error can not find DST trackmap node SvtxTrackMap" << std::endl;
0114       exit(-1);
0115     }
0116     }
0117
0118   //get event centrality
0119   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0120   if (!cent_node)
0121     {
0122       std::cout
0123         << "TracksInJets::process_event - Error can not find centrality node "
0124         << std::endl;
0125       exit(-1);
0126     }
0127   int cent = cent_node->get_centile(CentralityInfo::PROP::mbd_NS);
0128
0129
0130   // Loop through jets
0131   for(auto jet : *jets)
0132     {
0133       if(!jet)
0134     {
0135       std::cout << "WARNING!!! Jet not found" << std::endl;
0136       continue;
0137     }
0138       //sum up tracks in jet
0139       TVector3 sumtrk(0, 0, 0);
0140       for (SvtxTrackMap::Iter iter = trackmap->begin();
0141        iter != trackmap->end();
0142        ++iter)
0143     {
0144       SvtxTrack* track = iter->second;
0145       float quality = track->get_quality();
0146       auto silicon_seed = track->get_silicon_seed();
0147       int nmvtxhits = 0;
0148       if (silicon_seed)
0149         {
0150           nmvtxhits = silicon_seed->size_cluster_keys();
0151         }
0152       if(track->get_pt() < m_trk_pt_cut || quality > 6 || nmvtxhits < 4) //do some basic quality selections on tracks
0153         {
0154           continue;
0155         }
0156
0157       // Calculate delta eta, delta phi, and dR
0158       TVector3 v(track->get_px(), track->get_py(), track->get_pz());
0159       double dEta = v.Eta() - jet->get_eta();
0160       double dPhi = v.Phi() - jet->get_phi();
0161       while (dPhi > M_PI) dPhi -= 2 * M_PI;
0162       while (dPhi < -M_PI) dPhi += 2 * M_PI;
0163       double dR = sqrt(dEta * dEta + dPhi * dPhi);
0164
0165       // Check if track is within jet radius
0166       if (dR < m_jetRadius) {
0167     sumtrk += v;
0168       }
0169     }
0170
0171     // Fill histogram for the current jet
0172     assert(m_h_track_vs_calo_pt);
0173     assert(m_h_track_pt);
0174     // Fill TH3 histogram for Au+Au collisions
0175     if (isAA) {
0176       m_h_track_vs_calo_pt->Fill(jet->get_pt(), sumtrk.Perp(), cent);
0177     }
0178
0179     // Fill TH2 histogram for pp collisions
0180     else {
0181       m_h_track_pt->Fill(jet->get_pt(), sumtrk.Perp());
0182     }
0183     // Reset sumtrk for the next jet
0184     sumtrk.SetXYZ(0, 0, 0);
0185   }
0186   return Fun4AllReturnCodes::EVENT_OK;
0187 }
0188
0189 //____________________________________________________________________________..
0190 int TracksInJets::ResetEvent(PHCompositeNode *topNode)
0191 {
0192   //std::cout << "TracksInJets::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0193   return Fun4AllReturnCodes::EVENT_OK;
0194 }
0195
0196 //____________________________________________________________________________..
0197 int TracksInJets::EndRun(const int runnumber)
0198 {
0199   std::cout << "TracksInJets::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0200   return Fun4AllReturnCodes::EVENT_OK;
0201 }
0202
0203 //____________________________________________________________________________..
0204 int TracksInJets::End(PHCompositeNode *topNode)
0205 {
0206   std::cout << "TracksInJets::End - Output to " << m_outputFileName << std::endl;
0207   PHTFileServer::get().cd(m_outputFileName);
0208
0209   if (isAA) {
0210     TH2 *h_proj;
0211     for(int i = 0; i < m_h_track_vs_calo_pt->GetNbinsZ(); i++)
0212       {
0213     m_h_track_vs_calo_pt->GetZaxis()->SetRange(i+1,i+1);
0214     h_proj = (TH2*) m_h_track_vs_calo_pt->Project3D("yx");
0215     h_proj->Write(Form("h_track_vs_calo_%1.0f",m_h_track_vs_calo_pt->GetZaxis()->GetBinLowEdge(i+1)));
0216       }
0217   }
0218   else {
0219     m_h_track_pt->Write();  //if pp, do not project onto centrality bins
0220   }
0221   std::cout << "TracksInJets::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0222   return Fun4AllReturnCodes::EVENT_OK;
0223 }
0224
0225 //____________________________________________________________________________..
0226 int TracksInJets::Reset(PHCompositeNode *topNode)
0227 {
0228   std::cout << "TracksInJets::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0229   return Fun4AllReturnCodes::EVENT_OK;
0230 }
0231
0232 //____________________________________________________________________________..
0233 void TracksInJets::Print(const std::string &what) const
0234 {
0235   std::cout << "TracksInJets::Print(const std::string &what) const Printing info for " << what << std::endl;
0236 }