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