QA/Jet/DijetQA.cc

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