QA/Jet/UEvsEtaCentrality.cc

0001 ///===========================================================================
0002 /*!Migrated from analysis/JS-Jet/JetUESize-v2Psi2Centrality
0003  * Original Author: Micah Meskowitz
0004  * Author: Jinglin Liu
0005  * First Commit Date: 08.17.2025
0006  *
0007  *
0008  * Underlying event QA module
0009  */
0010 ///===========================================================================
0011
0012 #include "UEvsEtaCentrality.h"
0013
0014 #include <fun4all/Fun4AllReturnCodes.h>
0015 #include <fun4all/PHTFileServer.h>
0016
0017 #include <phool/PHCompositeNode.h>
0018
0019 #include <phool/getClass.h>
0020
0021 //#include <g4jets/JetMap.h>
0022 //#include <g4jets/Jetv1.h>
0023
0024 #include <centrality/CentralityInfo.h>
0025
0026 #include <calobase/RawTower.h>
0027 #include <calobase/RawTowerContainer.h>
0028 #include <calobase/RawTowerGeom.h>
0029 #include <calobase/RawTowerGeomContainer.h>
0030 #include <calobase/TowerInfoContainer.h>
0031 #include <calobase/TowerInfo.h>
0032
0033 #include <calotrigger/TriggerAnalyzer.h>
0034
0035 #include <jetbackground/TowerBackground.h>
0036
0037 #include <qautils/QAHistManagerDef.h>
0038
0039 #include <TTree.h>
0040 #include <TStyle.h>
0041 #include <TH1F.h>
0042 #include <TH2F.h>
0043
0044 //____________________________________________________________________________..
0045 UEvsEtaCentrality::UEvsEtaCentrality(const std::string& modulename)
0046   : SubsysReco(modulename)
0047   , m_moduleName(modulename)
0048 {
0049   if (m_config.debug && (Verbosity() > 1)){
0050     std::cout << "UEvsEtaCentrality::UEvsEtaCentrality(const std::string &moduleName) Calling ctor" << std::endl;
0051   }
0052 }
0053
0054 //____________________________________________________________________________..
0055 // Module constructor accepting a configuration
0056 UEvsEtaCentrality::UEvsEtaCentrality(const Config& config)
0057   : SubsysReco(config.moduleName)
0058   , m_config(config)
0059 {
0060   // print debug message
0061   if (m_config.debug && (Verbosity() > 1))
0062   {
0063     std::cout << "UEvsEtaCentrality::UEvsEtaCentrality(Config&) Calling ctor" << std::endl;
0064   }
0065 }  // end ctor(Config&)
0066
0067 //____________________________________________________________________________..
0068 UEvsEtaCentrality::~UEvsEtaCentrality()
0069 {
0070   // print debug message
0071   if (m_config.debug && (Verbosity() > 1)){
0072     std::cout << "UEvsEtaCentrality::~UEvsEtaCentrality() Calling dtor" << std::endl;
0073   }
0074 }
0075
0076 //____________________________________________________________________________..
0077 int UEvsEtaCentrality::Init(PHCompositeNode* /*topNode*/)
0078 {
0079   if (m_config.debug && (Verbosity() > 1)){
0080     std::cout << "UEvsEtaCentrality::Init(PHCompositeNode *topNode) Initializing" << std::endl;
0081   }
0082
0083   // initialize trigger analyzer
0084   delete m_analyzer;
0085   m_analyzer = new TriggerAnalyzer();
0086
0087   InitHistManager();
0088
0089   // make sure module name is lower case
0090   std::string smallModuleName = m_moduleName;
0091   std::transform(
0092       smallModuleName.begin(),
0093       smallModuleName.end(),
0094       smallModuleName.begin(),
0095       ::tolower);
0096
0097   // construct histogram names
0098   std::vector<std::string> vecHistNames = {
0099       "v2_cent",
0100       "psi2_cent",
0101       "ihcaleta_cent0_20",
0102       "ohcaleta_cent0_20",
0103       "emcaleta_cent0_20",
0104       "ihcaleta_cent20_50",
0105       "ohcaleta_cent20_50",
0106       "emcaleta_cent20_50",
0107       "ihcaleta_cent50_100",
0108       "ohcaleta_cent50_100",
0109       "emcaleta_cent50_100",
0110       "ihcaleta",
0111       "ohcaleta",
0112       "emcaleta",
0113   };
0114   for (auto &vecHistName : vecHistNames)
0115   {
0116     vecHistName.insert(0, "h_" + smallModuleName + "_");
0117     if (!m_config.histTag.empty())
0118     {
0119       vecHistName.append("_" + m_config.histTag);
0120     }
0121   }
0122
0123   // Initiate histograms
0124   hv2_cent = new TH2F(vecHistNames[0].data(),"",10,0,100,50,0,0.5);
0125   hPsi2_cent = new TH2F(vecHistNames[1].data(),"",10,0,100,50,-1.57,1.57);
0126   hUEiHcalEta_Cent0_20 = new TH2F(vecHistNames[2].data(), "", 24, -1.1,1.1, 48, 0, 0.25);
0127   hUEoHcalEta_Cent0_20 = new TH2F(vecHistNames[3].data(), "", 24, -1.1,1.1, 48, 0, 0.5);
0128   hUEemcalEta_Cent0_20 = new TH2F(vecHistNames[4].data(), "", 24, -1.1,1.1, 48, 0, 1.5);
0129
0130   hUEiHcalEta_Cent20_50 = new TH2F(vecHistNames[5].data(), "", 24, -1.1,1.1, 48, 0, 0.25);
0131   hUEoHcalEta_Cent20_50 = new TH2F(vecHistNames[6].data(), "", 24, -1.1,1.1, 48, 0, 0.5);
0132   hUEemcalEta_Cent20_50 = new TH2F(vecHistNames[7].data(), "", 24, -1.1,1.1, 48, 0, 1.5);
0133
0134   hUEiHcalEta_Cent50_100 = new TH2F(vecHistNames[8].data(), "", 24, -1.1,1.1, 48, 0, 0.25);
0135   hUEoHcalEta_Cent50_100 = new TH2F(vecHistNames[9].data(), "", 24, -1.1,1.1, 48, 0, 0.5);
0136   hUEemcalEta_Cent50_100 = new TH2F(vecHistNames[10].data(), "", 24, -1.1,1.1, 48, 0, 1.5);
0137
0138   hUEiHcalEta = new TH2F(vecHistNames[11].data(), "", 24, -1.1,1.1, 48, 0, 0.25);
0139   hUEoHcalEta = new TH2F(vecHistNames[12].data(), "", 24, -1.1,1.1, 48, 0, 0.5);
0140   hUEemcalEta = new TH2F(vecHistNames[13].data(), "", 24, -1.1,1.1, 48, 0, 1.5);
0141
0142   return Fun4AllReturnCodes::EVENT_OK;
0143 }
0144
0145 //____________________________________________________________________________..
0146 int UEvsEtaCentrality::InitRun(PHCompositeNode* /*topNode*/)
0147 {
0148   if (m_config.debug && (Verbosity() > 1))
0149   {
0150     std::cout << "UEvsEtaCentrality::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
0151   }
0152   return Fun4AllReturnCodes::EVENT_OK;
0153 }
0154
0155 //____________________________________________________________________________..
0156 int UEvsEtaCentrality::process_event(PHCompositeNode *topNode)
0157 {
0158
0159   if (m_config.debug && (Verbosity() > 1))
0160   {
0161     std::cout << "UEvsEtaCentrality::process_event(PHCompositeNode* topNode) Processing Event" << std::endl;
0162   }
0163
0164   // if needed, check if selected trigger fired
0165   if (m_config.doTrgSelect)
0166   {
0167     m_analyzer->decodeTriggers(topNode);
0168     bool hasTrigger = JetQADefs::DidTriggerFire(m_config.trgToSelect, m_analyzer);
0169     if (!hasTrigger)
0170     {
0171       return Fun4AllReturnCodes::EVENT_OK;
0172     }
0173   }
0174
0175   //centrality
0176   CentralityInfo* cent_node = findNode::getClass<CentralityInfo>(topNode, "CentralityInfo");
0177   if (!cent_node)
0178   {
0179     std::cout
0180       << "UEvsEtaCentrality::process_event - Error can not find centrality node "
0181       << std::endl;
0182     exit(-1);
0183   }
0184
0185   //calorimeter towers
0186   TowerInfoContainer *towersEM3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER");
0187   TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0188   TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0189   RawTowerGeomContainer *tower_geom = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0190   RawTowerGeomContainer *tower_geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0191   if(!towersEM3 || !towersIH3 || !towersOH3){
0192     std::cout
0193       <<"UEvsEtaCentrality::process_event - Error can not find raw tower node "
0194       << std::endl;
0195     exit(-1);
0196   }
0197
0198   if(!tower_geom || !tower_geomOH){
0199     std::cout
0200       <<"UEvsEtaCentrality::process_event - Error can not find raw tower geometry "
0201       << std::endl;
0202     exit(-1);
0203   }
0204
0205   //underlying event
0206   TowerBackground *background = findNode::getClass<TowerBackground>(topNode, "TowerInfoBackground_Sub2");
0207   if(!background){
0208     std::cout<<"Can't get background. Exiting"<<std::endl;
0209     return Fun4AllReturnCodes::EVENT_OK;
0210   }
0211
0212   float background_v2 = 0;
0213   float background_Psi2 = 0;
0214   float m_centrality = 0;
0215
0216   m_centrality =  cent_node->get_centile(CentralityInfo::PROP::mbd_NS);
0217   background_v2 = background->get_v2();
0218   background_Psi2 = background->get_Psi2();
0219
0220   // print centrality debug message
0221   if (m_config.debug && (Verbosity() > 1))
0222   {
0223     std::cout << "UEvsEtaCentrality::process_event, centrality = " << m_centrality << std::endl;
0224   }
0225
0226   hv2_cent->Fill(m_centrality, background_v2);
0227   hPsi2_cent->Fill(m_centrality, background_Psi2);
0228
0229   int i=0;
0230   for (i=0;i<=23;i++){
0231     float UEi = background->get_UE(1).at(i); //24 eta bins for each detector Inner Hcal
0232     float UEo = background->get_UE(2).at(i); //24 eta bins for each detector, Outer Hcal
0233     float UEe = background->get_UE(0).at(i); //24 eta bins for each detector, EMCal
0234
0235     double eta = tower_geom->get_etacenter(i); //eta comes from IHCal
0236
0237     // all centrality histograms
0238     hUEiHcalEta->Fill(eta, UEi);
0239     hUEoHcalEta->Fill(eta, UEo);
0240     hUEemcalEta->Fill(eta, UEe);
0241
0242     if (m_centrality > 0 && m_centrality <= 20){
0243       hUEiHcalEta_Cent0_20->Fill(eta, UEi);
0244       hUEoHcalEta_Cent0_20->Fill(eta, UEo);
0245       hUEemcalEta_Cent0_20->Fill(eta, UEe);
0246     }
0247     if (m_centrality > 20 && m_centrality <= 50){
0248       hUEiHcalEta_Cent20_50->Fill(eta, UEi);
0249       hUEoHcalEta_Cent20_50->Fill(eta, UEo);
0250       hUEemcalEta_Cent20_50->Fill(eta, UEe);
0251     }
0252     if (m_centrality > 50 && m_centrality <= 100){
0253       hUEiHcalEta_Cent50_100->Fill(eta, UEi);
0254       hUEoHcalEta_Cent50_100->Fill(eta, UEo);
0255       hUEemcalEta_Cent50_100->Fill(eta, UEe);
0256     }
0257   } // end of eta bin loop
0258
0259   // std::cout << "UEvsEtaCentrality::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
0260   return Fun4AllReturnCodes::EVENT_OK;
0261 } // end process_event
0262
0263 //____________________________________________________________________________..
0264 int UEvsEtaCentrality::ResetEvent(PHCompositeNode* /*topNode*/)
0265 {
0266   // std::cout << "UEvsEtaCentrality::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
0267
0268   // IMPORTANT!! CLEAR YOUR VECTORS AND RESET YOUR TREE VARIABLES HERE!!
0269
0270   return Fun4AllReturnCodes::EVENT_OK;
0271 }
0272
0273 //____________________________________________________________________________..
0274 int UEvsEtaCentrality::EndRun(const int runnumber)
0275 {
0276   if (m_config.debug && (Verbosity() > 1))
0277   {
0278     std::cout << "UEvsEtaCentrality::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
0279   }
0280
0281   return Fun4AllReturnCodes::EVENT_OK;
0282 }
0283
0284 //____________________________________________________________________________..
0285 int UEvsEtaCentrality::End(PHCompositeNode* /*topNode*/)
0286 {
0287
0288   m_manager->registerHisto(hv2_cent);
0289   m_manager->registerHisto(hPsi2_cent);
0290   m_manager->registerHisto(hUEiHcalEta_Cent0_20);
0291   m_manager->registerHisto(hUEoHcalEta_Cent0_20);
0292   m_manager->registerHisto(hUEemcalEta_Cent0_20);
0293   m_manager->registerHisto(hUEiHcalEta_Cent20_50);
0294   m_manager->registerHisto(hUEoHcalEta_Cent20_50);
0295   m_manager->registerHisto(hUEemcalEta_Cent20_50);
0296   m_manager->registerHisto(hUEiHcalEta_Cent50_100);
0297   m_manager->registerHisto(hUEoHcalEta_Cent50_100);
0298   m_manager->registerHisto(hUEemcalEta_Cent50_100);
0299   m_manager->registerHisto(hUEiHcalEta);
0300   m_manager->registerHisto(hUEoHcalEta);
0301   m_manager->registerHisto(hUEemcalEta);
0302
0303   if (m_config.debug && (Verbosity() > 1))
0304   {
0305     std::cout << "UEvsEtaCentrality::End(PHCompositeNode *topNode) This is the End..." << std::endl;
0306   }
0307
0308   return Fun4AllReturnCodes::EVENT_OK;
0309 }
0310
0311 //____________________________________________________________________________..
0312 int UEvsEtaCentrality::Reset(PHCompositeNode* /*topNode*/)
0313 {
0314   if (m_config.debug && (Verbosity() > 1))
0315   {
0316     std::cout << "UEvsEtaCentrality::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
0317   }
0318
0319   return Fun4AllReturnCodes::EVENT_OK;
0320 }
0321
0322 //____________________________________________________________________________..
0323 void UEvsEtaCentrality::Print(const std::string &what) const
0324 {
0325   std::cout << "UEvsEtaCentrality::Print(const std::string &what) const Printing info for " << what << std::endl;
0326 }
0327
0328 // Private methods =============================================================
0329
0330 //____________________________________________________________________________..
0331 // Initialize histogram manager
0332 void UEvsEtaCentrality::InitHistManager()
0333 {
0334
0335   // print debug message
0336   if (m_config.debug && (Verbosity() > 1))
0337   {
0338     std::cout << "UEvsEtaCentrality::InitHistManager() Initializing histogram manager" << std::endl;
0339   }
0340
0341   gStyle->SetOptTitle(0);
0342   m_manager = QAHistManagerDef::getHistoManager();
0343   if (!m_manager)
0344   {
0345     std::cerr << PHWHERE << ": PANIC! Couldn't grab histogram manager!" << std::endl;
0346     assert(m_manager);
0347   }
0348   return;
0349
0350 }  // end 'InitHistManager()'