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