TracksInJets/src/StructureinJets.cc

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