QA/Jet/StructureinJets.cc

0001 #include "StructureinJets.h"
0002
0003 #include <calotrigger/TriggerAnalyzer.h>
0004
0005 #include <centrality/CentralityInfo.h>
0006
0007 #include <jetbase/Jet.h>
0008 #include <jetbase/JetContainer.h>
0009 #include <jetbase/JetInput.h>
0010
0011 #include <qautils/QAHistManagerDef.h>
0012
0013 #include <trackbase_historic/SvtxTrack.h>
0014 #include <trackbase_historic/SvtxTrackMap.h>
0015 #include <trackbase_historic/TrackSeed.h>
0016
0017 #include <fun4all/Fun4AllHistoManager.h>
0018 #include <fun4all/Fun4AllReturnCodes.h>
0019 #include <fun4all/PHTFileServer.h>
0020
0021 #include <phool/PHCompositeNode.h>
0022 #include <phool/getClass.h>
0023
0024 #include <TH1.h>
0025 #include <TH2.h>
0026 #include <TH3.h>
0027 #include <TVector3.h>
0028 #include <TStyle.h>//for gStyle
0029
0030 #include <algorithm>
0031 #include <cassert>
0032 #include <cmath>
0033 #include <cstdlib>
0034 #include <iostream>
0035 #include <map>
0036 #include <utility>
0037
0038 //____________________________________________________________________________..
0039 StructureinJets::StructureinJets(const std::string& moduleName, const std::string& recoJetName, const std::string& trkNodeName, const std::string& histTag, const std::string& outputfilename)
0040   : SubsysReco(moduleName)
0041   , m_moduleName(moduleName)
0042   , m_recoJetName(recoJetName)
0043   , m_trkNodeName(trkNodeName)
0044   , m_histTag(histTag)
0045   , m_outputFileName(outputfilename)
0046 {
0047   if(Verbosity() > 1 )
0048   {
0049     std::cout << "StructureinJets::StructureinJets(const std::string& x 5) Calling ctor" << std::endl;
0050   }
0051 }
0052
0053 //____________________________________________________________________________..
0054 StructureinJets::~StructureinJets()
0055 {
0056   if (Verbosity() > 1)
0057   {
0058     std::cout << "StructureinJets::~StructureinJets() Calling dtor" << std::endl;
0059   }
0060 }
0061
0062 //____________________________________________________________________________..
0063 int StructureinJets::Init(PHCompositeNode* /*topNode*/)
0064 {
0065   if (Verbosity() > 0)
0066   {
0067     std::cout << "StructureinJets::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0068   }
0069
0070   if (m_writeToOutputFileFlag)
0071   {
0072     PHTFileServer::open(m_outputFileName, "RECREATE");
0073   }
0074
0075   delete m_analyzer;
0076   m_analyzer = new TriggerAnalyzer();
0077   m_manager = QAHistManagerDef::getHistoManager();
0078
0079   // make sure module name is lower case
0080   std::string smallModuleName = m_moduleName;
0081   std::transform(
0082       smallModuleName.begin(),
0083       smallModuleName.end(),
0084       smallModuleName.begin(),
0085       ::tolower);
0086
0087   // construct histogram names
0088   std::vector<std::string> vecHistNames = {
0089       "sumtrkvsjetptvscent",
0090       "sumtrkvsjetpt",
0091       "sumtrkoverjtpt",
0092       "sumtrkpt"};
0093   for (auto& vecHistName : vecHistNames)
0094   {
0095     vecHistName.insert(0, "h_" + smallModuleName + "_");
0096     if (!m_histTag.empty())
0097     {
0098       vecHistName.append("_" + m_histTag);
0099     }
0100   }
0101
0102   m_hSumTrkVsJetPtVsCent = new TH3F(vecHistNames[0].data(), "", 100, 0, 100, 200, 0, 100, 10, 0, 100);
0103   m_hSumTrkVsJetPtVsCent->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0104   m_hSumTrkVsJetPtVsCent->GetYaxis()->SetTitle("Sum track p_{T} [GeV]");
0105   m_hSumTrkVsJetPt = new TH2F(vecHistNames[1].data(), "", 100, 0, 100, 200, 0, 100);
0106   m_hSumTrkVsJetPt->GetXaxis()->SetTitle("Jet p_{T} [GeV]");
0107   m_hSumTrkVsJetPt->GetYaxis()->SetTitle("Sum track p_{T} [GeV]");
0108   m_hSumTrkOverJetPt = new TH1F(vecHistNames[2].data(), "", 50, 0, 5);
0109   m_hSumTrkOverJetPt->GetXaxis()->SetTitle("Sum track p_{T} / Jet p_{T}");
0110   m_hSumTrkOverJetPt->GetYaxis()->SetTitle("Counts");
0111   m_hSumTrkPt = new TH1F(vecHistNames[3].data(), "", 200, 0, 100);
0112   m_hSumTrkPt->GetXaxis()->SetTitle("Sum track p_{T}");
0113   m_hSumTrkPt->GetYaxis()->SetTitle("Counts");
0114
0115   if (m_isAAFlag)
0116   {
0117     m_manager->registerHisto(m_hSumTrkVsJetPtVsCent);
0118   }
0119   m_manager->registerHisto(m_hSumTrkVsJetPt);
0120   m_manager->registerHisto(m_hSumTrkOverJetPt);
0121   m_manager->registerHisto(m_hSumTrkPt);
0122   return Fun4AllReturnCodes::EVENT_OK;
0123 }
0124
0125 //____________________________________________________________________________..
0126 int StructureinJets::InitRun(PHCompositeNode* /*topNode*/)
0127 {
0128   if (Verbosity() > 0)
0129   {
0130     std::cout << "StructureinJets::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0131   }
0132   return Fun4AllReturnCodes::EVENT_OK;
0133 }
0134
0135 //____________________________________________________________________________..
0136 int StructureinJets::process_event(PHCompositeNode* topNode)
0137 {
0138
0139   if (Verbosity() > 1)
0140   {
0141     std::cout << "StructureinJets::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0142   }
0143
0144   // if needed, check if selected trigger fired
0145   if (m_doTrgSelect)
0146   {
0147     m_analyzer->decodeTriggers(topNode);
0148     bool hasTrigger = JetQADefs::DidTriggerFire(m_trgToSelect, m_analyzer);
0149     if (!hasTrigger)
0150     {
0151       return Fun4AllReturnCodes::EVENT_OK;
0152     }
0153   }
0154
0155   // interface to reco jets
0156   JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName);
0157   if (!jets)
0158   {
0159     std::cout
0160         << "StructureInJets::process_event - Error can not find DST Reco JetContainer node "
0161         << m_recoJetName << std::endl;
0162     return Fun4AllReturnCodes::EVENT_OK;
0163   }
0164
0165   // get reco tracks
0166   SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, m_trkNodeName);
0167   if (!trackmap)
0168   {
0169     // if specified node name not found, try looking for this node;
0170     // and if that's not found, exit
0171     trackmap = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0172     if (!trackmap)
0173     {
0174       std::cout
0175           << "StructureinJets::process_event - Error can not find DST trackmap node SvtxTrackMap" << std::endl;
0176       return Fun4AllReturnCodes::EVENT_OK;
0177     }
0178   }
0179
0180   // get event centrality
0181   int cent = -1;
0182   if (m_isAAFlag)
0183   {
0184     CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0185     if (!cent_node)
0186     {
0187       std::cout
0188           << "StructureinJets::process_event - Error can not find centrality node "
0189           << std::endl;
0190       return Fun4AllReturnCodes::EVENT_OK;
0191     }
0192     cent = cent_node->get_centile(CentralityInfo::PROP::mbd_NS);
0193   }
0194
0195   // Loop through jets
0196   for (auto *jet : *jets)
0197   {
0198     if (!jet)
0199     {
0200       if (Verbosity() > 2)
0201       {
0202         std::cout << "WARNING!!! Jet not found" << std::endl;
0203       }
0204       continue;
0205     }
0206
0207     // remove noise
0208     if (jet->get_pt() < m_ptJetRange.first)
0209     {
0210       continue;
0211     }
0212
0213     // apply jet kinematic cuts
0214     const bool inJetEtaCut = (jet->get_eta() >= m_etaJetRange.first) && (jet->get_eta() <= m_etaJetRange.second);
0215     const bool inJetPtCut = (jet->get_pt() >= m_ptJetRange.first) && (jet->get_pt() <= m_ptJetRange.second);
0216     if (!inJetEtaCut || !inJetPtCut)
0217     {
0218       continue;
0219     }
0220
0221     // sum up tracks in jet
0222     TVector3 sumtrk(0, 0, 0);
0223     for (auto& iter : *trackmap)
0224     {
0225       SvtxTrack* track = iter.second;
0226       float quality = track->get_quality();
0227       auto *silicon_seed = track->get_silicon_seed();
0228       auto *tpc_seed = track->get_tpc_seed();
0229
0230       // get no. of clusters in silicon seed
0231       int nsiliconclusts = 0;
0232       if (silicon_seed)
0233       {
0234         nsiliconclusts = silicon_seed->size_cluster_keys();
0235       }
0236
0237       // get no. of clusters in tpc seed
0238       int ntpcclusts = 0;
0239       if (tpc_seed)
0240       {
0241         ntpcclusts = tpc_seed->size_cluster_keys();
0242       }
0243
0244       // do some basic quality selection on tracks
0245       const bool inTrkPtCut = (track->get_pt() >= m_trk_pt_cut);
0246       const bool inTrkQualCut = (quality <= m_trk_qual_cut);
0247       const bool inTrkNSilCut = (nsiliconclusts >= m_trk_nsil_cut);
0248       const bool inTrkNTPCCut = (ntpcclusts >= m_trk_ntpc_cut);
0249       if (!inTrkPtCut || !inTrkQualCut || !inTrkNSilCut || !inTrkNTPCCut)
0250       {
0251         continue;
0252       }
0253
0254       // Calculate delta eta, delta phi, and dR
0255       TVector3 v(track->get_px(), track->get_py(), track->get_pz());
0256       double dEta = v.Eta() - jet->get_eta();
0257       double dPhi = v.Phi() - jet->get_phi();
0258       while (dPhi > M_PI)
0259       {
0260         dPhi -= 2 * M_PI;
0261       }
0262       while (dPhi < -M_PI)
0263       {
0264         dPhi += 2 * M_PI;
0265       }
0266       double dR = sqrt((dEta * dEta) + (dPhi * dPhi));
0267
0268       // Check if track is within jet radius
0269       if (dR < m_jetRadius)
0270       {
0271         sumtrk += v;
0272       }
0273     }
0274
0275     // Fill histogram for the current jet
0276     assert(m_hSumTrkVsJetPtVsCent);
0277     assert(m_hSumTrkVsJetPt);
0278     assert(m_hSumTrkOverJetPt);
0279     assert(m_hSumTrkPt);
0280
0281     // Fill TH3 histogram for Au+Au collisions
0282     if (m_isAAFlag)
0283     {
0284       m_hSumTrkVsJetPtVsCent->Fill(jet->get_pt(), sumtrk.Perp(), cent);
0285     }
0286
0287     // Always fill others
0288     m_hSumTrkVsJetPt->Fill(jet->get_pt(), sumtrk.Perp());
0289     m_hSumTrkOverJetPt->Fill(sumtrk.Perp()/jet->get_pt());
0290     m_hSumTrkPt->Fill(sumtrk.Perp());
0291   }
0292   return Fun4AllReturnCodes::EVENT_OK;
0293 }
0294
0295 //____________________________________________________________________________..
0296 int StructureinJets::ResetEvent(PHCompositeNode* /*topNode*/)
0297 {
0298   if (Verbosity() > 1)
0299   {
0300     std::cout << "StructureinJets::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0301   }
0302   return Fun4AllReturnCodes::EVENT_OK;
0303 }
0304
0305 //____________________________________________________________________________..
0306 int StructureinJets::EndRun(const int runnumber)
0307 {
0308   if (Verbosity() > 0)
0309   {
0310   std::cout << "StructureinJets::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0311   }
0312   return Fun4AllReturnCodes::EVENT_OK;
0313 }
0314
0315 //____________________________________________________________________________..
0316 int StructureinJets::End(PHCompositeNode* /*topNode*/)
0317 {
0318   // if flag is true, write to output file
0319   // otherwise rely on histogram manager
0320   if (m_writeToOutputFileFlag)
0321   {
0322     if (Verbosity() > 1)
0323     {
0324       std::cout << "StructureinJets::End - Output to " << m_outputFileName << std::endl;
0325     }
0326     PHTFileServer::cd(m_outputFileName);
0327   }
0328   else
0329   {
0330     if (Verbosity() > 1)
0331     {
0332       std::cout << "StructureinJets::End - Output to histogram manager" << std::endl;
0333     }
0334   }
0335
0336   if (m_isAAFlag)
0337   {
0338     TH2* h_proj;
0339     for (int i = 0; i < m_hSumTrkVsJetPtVsCent->GetNbinsZ(); i++)
0340     {
0341       // construct histogram name for projection
0342       std::string name = m_hSumTrkVsJetPtVsCent->GetName();
0343       name.append("_centbin" + std::to_string(i + 1));
0344
0345       m_hSumTrkVsJetPtVsCent->GetZaxis()->SetRange(i + 1, i + 1);
0346       h_proj = (TH2*) m_hSumTrkVsJetPtVsCent->Project3D("yx");
0347       h_proj->SetName(name.data());
0348       if (m_writeToOutputFileFlag)
0349       {
0350         h_proj->Write();
0351       }
0352       else
0353       {
0354         m_manager->registerHisto(h_proj);
0355       }
0356     }
0357   }
0358   else
0359   {
0360     if (m_writeToOutputFileFlag)
0361     {
0362       m_hSumTrkVsJetPt->Write();  // if pp, do not project onto centrality bins
0363     }
0364   }
0365   if (Verbosity() > 0)
0366   {
0367     std::cout << "StructureinJets::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0368   }
0369   return Fun4AllReturnCodes::EVENT_OK;
0370 }
0371
0372 //____________________________________________________________________________..
0373 int StructureinJets::Reset(PHCompositeNode* /*topNode*/)
0374 {
0375   if (Verbosity() > 0)
0376   {
0377     std::cout << "StructureinJets::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0378   }
0379   return Fun4AllReturnCodes::EVENT_OK;
0380 }
0381
0382 //____________________________________________________________________________..
0383 void StructureinJets::Print(const std::string& what) const
0384 {
0385   std::cout << "StructureinJets::Print(const std::string &what) const Printing info for " << what << std::endl;
0386 }