Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:21:18 

0001 //////////////////////////////////////////////////////////////////
0002 //                                  //
0003 //      Dijet QA module for JetQA               //
0004 //      QA module that looks at dijet Ajj, xj       //
0005 //                              //
0006 //                              //
0007 //                              //
0008 //      Author:Skadi                    //
0009 //      First commit:       13 Nov 24       //
0010 //      Most recent update: 14 Dec 25       //
0011 //      version:        v2.0            //
0012 //////////////////////////////////////////////////////////////////
0013 #include "DijetQA.h"
0014 
0015 #include <calotrigger/TriggerAnalyzer.h>
0016 
0017 #include <globalvertex/GlobalVertex.h>
0018 #include <globalvertex/GlobalVertexMap.h>
0019 
0020 #include <jetbase/Jet.h>
0021 #include <jetbase/JetContainer.h>
0022 
0023 #include <qautils/QAHistManagerDef.h>
0024 
0025 #include <fun4all/Fun4AllHistoManager.h>
0026 #include <fun4all/Fun4AllReturnCodes.h>
0027 #include <fun4all/SubsysReco.h>
0028 
0029 #include <phool/PHCompositeNode.h>
0030 #include <phool/getClass.h>
0031 #include <phool/phool.h>
0032 
0033 #include <TH1.h>
0034 #include <TH2.h>
0035 #include <TStyle.h>
0036 
0037 #include <cassert>
0038 #include <format>
0039 #include <string>
0040 #include <utility>
0041 #include <vector>
0042 
0043 //____________________________________________________________________________..
0044 DijetQA::DijetQA(const std::string& name, const std::string& recojetname)
0045   : SubsysReco(name)
0046   , m_moduleName(name)
0047   , m_recoJetName(recojetname)
0048 {
0049   if (Verbosity() > 1)
0050   {
0051     std::cout << "DijetQA::DijetQA(const std::string &name) Calling ctor" << std::endl;
0052   }
0053 }
0054 
0055 //____________________________________________________________________________..
0056 int DijetQA::Init(PHCompositeNode* /*topNode*/)
0057 {
0058   //  std::cout << "DijetQA::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0059   delete m_analyzer;  // make cppcheck happy
0060   m_analyzer = new TriggerAnalyzer();
0061 
0062   gStyle->SetOptTitle(0);
0063   m_manager = QAHistManagerDef::getHistoManager();  // get the histogram anager
0064 
0065   if (!m_manager)
0066   {
0067     std::cout << PHWHERE << ": PANIC: couldn't grab histogram manager!" << std::endl;
0068     assert(m_manager);
0069   }
0070   std::string smallModuleName = m_moduleName;  // make sure name is lowercase
0071   std::transform(
0072       smallModuleName.begin(),
0073       smallModuleName.end(),
0074       smallModuleName.begin(),
0075       ::tolower);
0076   h_Ajj = new TH1F(std::format("h_{}_Ajj", smallModuleName).c_str(), std::format("A_{{jj}} for identified jet pairs for {}; A_{{jj}}; N_{{pairs}}", m_recoJetName).c_str(), 100, -0.005, 0.995);
0077   h_xj = new TH1F(std::format("h_{}_xj", smallModuleName).c_str(), std::format("x_{{j}} for identified jet pairs for {}; x_{{j}}; N_{{pairs}}", m_recoJetName).c_str(), 100, -0.005, 0.995);
0078   h_pt = new TH1F(std::format("h_{}_pt", smallModuleName).c_str(), std::format("p_{{T}} for leading jets in identified pairs for {}; p_{{T}} [GeV/c]; N_{{jet}}", m_recoJetName).c_str(), 70, -0.5, 69.5);
0079   h_dphi = new TH1F(std::format("h_{}_dphi", smallModuleName).c_str(), std::format("#Delta #varphi for identified jet pairs for {}; #Delta #varphi; N_{{pairs}}", m_recoJetName).c_str(), 64, -M_PI, M_PI);
0080   h_Ajj_pt = new TH2F(std::format("h_{}_Ajj_pt", smallModuleName).c_str(), std::format("A_{{jj}} as a function of leading jet $p_{{T}}$ for {}; p_{{T}}^{{leading}} [GeV/c]; A_{{jj}}; N_{{pairs}}", m_recoJetName).c_str(), 70, -0.5, 69.5, 100, -0.005, 0.995);
0081   h_xj_pt = new TH2F(std::format("h_{}_xj_pt", smallModuleName).c_str(), std::format("x_{{j}} as a function of leading jet $p_{{T}}$ for {}; p_{{T}}^{{leading}} [GeV]; x_{{j}}; N_{{pairs}}", m_recoJetName).c_str(), 70, -0.5, 69.5, 100, -0.005, 0.995);
0082   h_dphi_pt = new TH2F(std::format("h_{}_dphi_pt", smallModuleName).c_str(), std::format("|#Delta #varphi| of dijet pair as a function of leading jet p_{{T}} for {}; p_{{T}}^{{leading}} [GeV/c]; |#Delta #varphi|; N_{{pairs}}", m_recoJetName).c_str(), 70, -0.5, 69.5, 64, 0, M_PI);
0083   h_dphi_Ajj = new TH2F(std::format("h_{}_dphi_Ajj", smallModuleName).c_str(), std::format("A_{{jj}} of dijet pair as a function of |#Delta #varphi| for {}; |#Delta #varphi|; A_{{jj}}; N_{{pairs}}", m_recoJetName).c_str(), 64, 0, M_PI, 100, -0.005, 0.995);
0084   h_Ajj_l = new TH1F(std::format("h_{}_Ajj_l", smallModuleName).c_str(), std::format("A_{{jj}} for event leading jet pairs for {}; A_{{jj}}; N_{{pairs}}", m_recoJetName).c_str(), 100, -0.005, 0.995);
0085   h_xj_l = new TH1F(std::format("h_{}_xj_l", smallModuleName).c_str(), std::format("x_{{j}} for event leading jet pairs for {}; x_{{j}}; N_{{pairs}}", m_recoJetName).c_str(), 100, -0.005, 0.995);
0086   h_pt_l = new TH1F(std::format("h_{}_pt_l", smallModuleName).c_str(), std::format("p_{{T}} for leading jets in event leading pair for {}; p_{{T}} [GeV/c]; N_{{jet}}", m_recoJetName).c_str(), 70, -0.5, 69.5);
0087   h_dphi_l = new TH1F(std::format("h_{}_dphi_l", smallModuleName).c_str(), std::format("#Delta #varphi for leading jet pairs for {}; #Delta #varphi; N_{{pairs}}", m_recoJetName).c_str(), 64, -M_PI, M_PI);
0088   h_Ajj_pt_l = new TH2F(std::format("h_{}_Ajj_pt_l", smallModuleName).c_str(), std::format("A_{{jj}} of event leading dijet pair as a function of leading jet p_{{T}} for {}; p_{{T}}^{{leading}} [GeV/c]; A_{{jj}}; N_{{pairs}}", m_recoJetName).c_str(), 70, -0.5, 69.5, 100, -0.005, 0.995);
0089   h_xj_pt_l = new TH2F(std::format("h_{}_xj_pt_l", smallModuleName).c_str(), std::format("x_{{j}} of event leading dijet pair as a function of leading jet p_{{T}} for {}; p_{{T}}^{{leading}} [GeV/c]; x_{{j}}; N_{{pairs}}", m_recoJetName).c_str(), 70, -0.5, 69.5, 100, -0.005, 0.995);
0090   h_dphi_pt_l = new TH2F(std::format("h_{}_dphi_pt_l", smallModuleName).c_str(), std::format("|#Delta #varphi| of event leading dijet pair as a function of leading jet p_{{T}} for {}; p_{{T}}^{{leading}} [GeV/c]; |#Delta #varphi|; N_{{pairs}}", m_recoJetName).c_str(), 70, -0.5, 69.5, 64, 0, M_PI);
0091   h_dphi_Ajj_l = new TH2F(std::format("h_{}_dphi_Ajj_l", smallModuleName).c_str(), std::format("A_{{jj}} of event leading dijet pair as a function of |#Delta #varphi| for {}; |#Delta #varphi|^{{leading}}; A_{{jj}}; N_{{pairs}}", m_recoJetName).c_str(), 64, 0, M_PI, 100, -0.005, 0.995);
0092 
0093   return Fun4AllReturnCodes::EVENT_OK;
0094 }
0095 
0096 //____________________________________________________________________________..
0097 int DijetQA::process_event(PHCompositeNode* topNode)
0098 {
0099   if (Verbosity() > 1)
0100   {
0101     std::cout << "DijetQA::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0102   }
0103   if (m_doTrgSelect)
0104   {
0105     m_analyzer->decodeTriggers(topNode);
0106     bool hasTrigger = JetQADefs::DidTriggerFire(m_trgToSelect, m_analyzer);
0107     if (!hasTrigger)
0108     {
0109       return Fun4AllReturnCodes::EVENT_OK;
0110     }
0111   }
0112   /* //removing the centrality for now, will add back in with a conditional flag for pp
0113     CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0114     if (!cent_node)
0115     {
0116       if (Verbosity() > 1)
0117       {
0118         std::cout << "No centrality info found" << std::endl;
0119       }
0120       m_centrality = 1.;
0121       m_impactparam = 0.;
0122     }
0123     else
0124     {
0125       m_centrality = cent_node->get_centile(CentralityInfo::PROP::bimp);
0126       m_impactparam = cent_node->get_quantity(CentralityInfo::PROP::bimp);
0127     }*/
0128   GlobalVertex* vtx = nullptr;
0129   GlobalVertexMap* vtxmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0130   if (!vtxmap || vtxmap->empty())
0131   {
0132     if (!vtxmap)
0133     {
0134       std::cout << "DijetQA::process_event - Error can not find vtxmap node " << "GlobalVertexMap" << std::endl;
0135     }
0136     if (Verbosity() > 1)
0137     {
0138       if (vtxmap->empty())
0139       {
0140         std::cout << "No vertex map found, assuming the vertex has z=0" << std::endl;
0141       }
0142     }
0143     m_zvtx = 0;
0144   }
0145   else
0146   {
0147     vtx = vtxmap->begin()->second;
0148     m_zvtx = vtx->get_z();
0149   }
0150   JetContainer* jets = findNode::getClass<JetContainer>(topNode, m_recoJetName);
0151   if (!jets)
0152   {
0153     std::cout << "DijetQA::process_event - Error can not find jets node " << m_recoJetName << std::endl;
0154     if (Verbosity() > 1)
0155     {
0156       std::cout << "No Jet container found" << std::endl;
0157     }
0158     return Fun4AllReturnCodes::EVENT_OK;
0159   }
0160 
0161   FindPairs(jets);
0162 
0163   return Fun4AllReturnCodes::EVENT_OK;
0164 }
0165 //____________________________________________________________________________..
0166 void DijetQA::FindPairs(JetContainer* jets)
0167 {
0168   // find all pairs that are potenital dijets and measure the kinematics
0169   if (Verbosity() > 1)
0170   {
0171     std::cout << "JetKinematicCheck::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0172   }
0173   Jet* jet_leading = nullptr;
0174   float pt_leading = 0;
0175   float pt1 = 0;
0176   float pt2 = 0;
0177   m_nJet = jets->size();
0178   if (Verbosity() > 2)
0179   {
0180     std::cout << "number of jets is" << m_nJet << std::endl;
0181   }
0182   std::vector<std::pair<Jet*, Jet*> > jet_pairs;
0183   bool set_leading = false;
0184   for (auto* j1 : *jets)
0185   {
0186     // assert(j1);
0187     Jet* jet_pair1 = nullptr;
0188     Jet* jet_pair2 = nullptr;
0189     if (j1->get_pt() < m_ptLeadRange.first || j1->get_eta() < m_etaRange.first || j1->get_pt() > m_ptLeadRange.second || j1->get_eta() > m_etaRange.second)
0190     {
0191       continue;  // cut on 1 GeV jets
0192     }
0193     if (j1->get_pt() > pt_leading)
0194     {
0195       set_leading = true;
0196       pt_leading = j1->get_pt();
0197       jet_leading = j1;
0198     }
0199     for (auto* j2 : (*jets))
0200     {
0201       if (j2 < j1)
0202       {
0203         continue;
0204       }
0205       if (/*j2 == j1 ||*/ j2->get_pt() < m_ptSubRange.first || j2->get_eta() < m_etaRange.first || j2->get_pt() > m_ptSubRange.second || j2->get_eta() > m_etaRange.second)
0206       {
0207         continue;
0208       }
0209       float dphil = j2->get_phi() - j1->get_phi();
0210       if (dphil > M_PI)
0211       {
0212         dphil = 2 * M_PI - dphil;
0213       }
0214       if (dphil < -M_PI)
0215       {
0216         dphil = -2 * M_PI - dphil;
0217       }
0218 
0219       if (std::abs(dphil) > M_PI - DeltaPhi)
0220       {
0221         if (j2->get_pt() > j1->get_pt())
0222         {
0223           jet_pair1 = j2;
0224           jet_pair2 = j1;
0225         }
0226         else
0227         {
0228           jet_pair1 = j1;
0229           jet_pair2 = j2;
0230         }
0231         jet_pairs.emplace_back(jet_pair1, jet_pair2);
0232         if (set_leading && jet_pair1 && jet_pair2)
0233         {
0234           pt1 = jet_pair1->get_pt();
0235           pt2 = jet_pair2->get_pt();
0236           m_Ajj = (pt1 - pt2) / (pt1 + pt2);
0237           m_xj = pt2 / pt1;
0238           m_ptl = pt1;
0239           m_ptsl = pt2;
0240           m_phil = jet_pair1->get_phi();
0241           m_phisl = jet_pair2->get_phi();
0242           m_dphil = dphil;
0243           m_etal = jet_pair1->get_eta();
0244           m_etasl = jet_pair2->get_eta();
0245           m_deltaeta = m_etal - m_etasl;
0246           h_Ajj_l->Fill(m_Ajj);
0247           h_xj_l->Fill(m_xj);
0248           h_pt_l->Fill(m_ptl);
0249           h_dphi_l->Fill(m_dphil);
0250           h_Ajj_pt_l->Fill(m_ptl, m_Ajj);
0251           h_xj_pt_l->Fill(m_ptl, m_xj);
0252           h_dphi_pt_l->Fill(m_ptl, std::abs(m_dphil));
0253           h_dphi_Ajj_l->Fill(std::abs(m_dphil), m_Ajj);
0254           //    m_T->Fill();
0255         }
0256         set_leading = false;
0257       }
0258     }
0259   }
0260   if (Verbosity() > 2)
0261   {
0262     std::cout << "Finished search for pairs" << std::endl;
0263   }
0264   m_nJetPair = jet_pairs.size();
0265   float Ajj = 0.;
0266   float xj = 0.;
0267   if (!jet_pairs.empty())
0268   {
0269     for (auto js : jet_pairs)
0270     {
0271       Jet* jet_pair1 = js.first;
0272       Jet* jet_pair2 = js.second;
0273       if (!jet_pair1 || !jet_pair2)
0274       {
0275         continue;
0276       }
0277       if (jet_pair1 && Verbosity() > 2)
0278       {
0279         std::cout << "jetpair 1 object has a pt of " << jet_pair1->get_pt() << std::endl;
0280       }
0281       pt1 = jet_pair1->get_pt();
0282       pt2 = jet_pair2->get_pt();
0283       if (pt1 < pt2)
0284       {
0285         auto* j = jet_pair1;
0286         jet_pair1 = jet_pair2;
0287         jet_pair2 = j;
0288         pt1 = jet_pair1->get_pt();
0289         pt2 = jet_pair2->get_pt();
0290       }
0291       float dphi = jet_pair1->get_phi() - jet_pair2->get_phi();
0292       if (dphi > M_PI)
0293       {
0294         dphi = 2 * M_PI - dphi;
0295       }
0296       if (dphi < -M_PI)
0297       {
0298         dphi = -2 * M_PI - dphi;
0299       }
0300       Ajj = (pt1 - pt2) / (pt1 + pt2);
0301       xj = pt2 / pt1;
0302       h_Ajj->Fill(Ajj);
0303       h_xj->Fill(xj);
0304       h_pt->Fill(pt1);
0305       h_dphi->Fill(dphi);
0306       h_Ajj_pt->Fill(pt1, Ajj);
0307       h_xj_pt->Fill(pt1, xj);
0308       h_dphi_pt->Fill(pt1, std::abs(dphi));
0309       h_dphi_Ajj->Fill(std::abs(dphi), Ajj);
0310       if (Verbosity() > 2)
0311       {
0312         std::cout << "highest pt jet is " << jet_leading->get_pt() << " and highest pt in a pair is " << jet_pair1->get_pt() << std::endl;
0313       }
0314     }
0315   }
0316   else
0317   {
0318     if (Verbosity() > 2)
0319     {
0320       std::cout << "Did not find a pair of jets" << std::endl;
0321     }
0322   }
0323   return;
0324 }
0325 
0326 //____________________________________________________________________________..
0327 int DijetQA::End(PHCompositeNode* /*topNode*/)
0328 {
0329   //  std::cout << "DijetQA::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0330   /*h_Ajj->SetStats(0);
0331   h_xj->SetStats(0);
0332   h_dphi->SetStats(0);
0333   h_Ajj_pt->SetStats(0);
0334   h_xj_pt->SetStats(0);
0335   h_dphi_pt->SetStats(0);
0336   h_dphi_Ajj->SetStats(0);
0337   h_Ajj_l->SetStats(0);
0338   h_xj_l->SetStats(0);
0339   h_dphi_l->SetStats(0);
0340   h_Ajj_pt_l->SetStats(0);
0341   h_xj_pt_l->SetStats(0);
0342   h_dphi_pt_l->SetStats(0);
0343   h_dphi_Ajj_l->SetStats(0);
0344   TLegend* l1=new TLegend(0.7, 0.9, 0.9, 1.0); //TLegend for the Ajj plots
0345   l1->SetFillStyle(0);
0346   l1->SetBorderSize(0);
0347   l1->SetTextSize(0.06f);
0348   l1->AddEntry((TObject*) nullptr, std::format("A_{{jj}} dijet pairs with pt_{{l}} #geq {}", m_ptLeadRange.first),"");
0349   h_Ajj->GetListOfFunctions()->Add(l1);*/
0350 
0351   m_manager->registerHisto(h_Ajj);
0352   m_manager->registerHisto(h_xj);
0353   // m_manager->registerHisto(h_pt); //this is turned off but can be turned on to diagnose issues
0354   m_manager->registerHisto(h_dphi);
0355   m_manager->registerHisto(h_Ajj_pt);
0356   m_manager->registerHisto(h_xj_pt);
0357   m_manager->registerHisto(h_dphi_pt);
0358   m_manager->registerHisto(h_dphi_Ajj);
0359   m_manager->registerHisto(h_Ajj_l);
0360   m_manager->registerHisto(h_xj_l);
0361   // m_manager->registerHisto(h_pt_l); //this is turned off but can be turned on to diagnose issues
0362   m_manager->registerHisto(h_dphi_l);
0363   m_manager->registerHisto(h_Ajj_pt_l);
0364   m_manager->registerHisto(h_xj_pt_l);
0365   m_manager->registerHisto(h_dphi_pt_l);
0366   m_manager->registerHisto(h_dphi_Ajj_l);
0367   return Fun4AllReturnCodes::EVENT_OK;
0368 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!