File indexing completed on 2026-05-23 08:13:54
0001
0002
0003
0004
0005 #include "EMCalShowerShapes.h"
0006
0007 #include <calobase/RawCluster.h>
0008 #include <calobase/RawClusterContainer.h>
0009 #include <calobase/RawClusterUtility.h>
0010 #include <calobase/RawTowerDefs.h>
0011 #include <calobase/RawTowerGeom.h>
0012 #include <calobase/RawTowerGeomContainer.h>
0013 #include <calobase/TowerInfo.h>
0014 #include <calobase/TowerInfoContainer.h>
0015 #include <calobase/TowerInfoDefs.h>
0016
0017 #include <calotrigger/TriggerAnalyzer.h>
0018
0019 #include <fun4all/Fun4AllHistoManager.h>
0020 #include <fun4all/Fun4AllReturnCodes.h>
0021
0022 #include <globalvertex/MbdVertex.h>
0023 #include <globalvertex/MbdVertexMap.h>
0024
0025 #include <phool/PHCompositeNode.h>
0026 #include <phool/getClass.h>
0027
0028 #include <qautils/QAHistManagerDef.h>
0029
0030 #include <TH1.h>
0031 #include <TH2.h>
0032 #include <TStyle.h>
0033 #include <TSystem.h>
0034
0035 #include <algorithm>
0036 #include <cassert>
0037 #include <cmath>
0038 #include <iostream>
0039 #include <set>
0040 #include <vector>
0041
0042 namespace
0043 {
0044 void shift_tower_index(int& ieta, int& iphi, int etadiv, int phidiv)
0045 {
0046 while (iphi < 0)
0047 {
0048 iphi += phidiv;
0049 }
0050 while (iphi >= phidiv)
0051 {
0052 iphi -= phidiv;
0053 }
0054 if (ieta < 0 || ieta >= etadiv)
0055 {
0056 ieta = -1;
0057 }
0058 }
0059 }
0060
0061 EMCalShowerShapes::EMCalShowerShapes(const std::string &modulename, const std::string &inputnode, const std::string &histtag)
0062 : SubsysReco(modulename)
0063 , m_modulename(modulename)
0064 , m_inputnode(inputnode)
0065 , m_histtag(histtag)
0066 , m_trgToSelect(JetQADefs::GL1::MBDNSPhoton1)
0067 , m_doTrgSelect(false)
0068 {
0069 }
0070
0071 EMCalShowerShapes::~EMCalShowerShapes()
0072 {
0073 delete m_analyzer;
0074 }
0075
0076 int EMCalShowerShapes::Init(PHCompositeNode* )
0077 {
0078 delete m_analyzer;
0079 m_analyzer = new TriggerAnalyzer();
0080
0081 m_manager = QAHistManagerDef::getHistoManager();
0082 if (!m_manager)
0083 {
0084 std::cerr << PHWHERE << "PANIC: couldn't grab histogram manager!" << std::endl;
0085 gSystem->Exit(1);
0086 }
0087
0088 std::string smallModuleName = m_modulename;
0089 std::transform(smallModuleName.begin(), smallModuleName.end(), smallModuleName.begin(), ::tolower);
0090
0091 std::vector<std::string> vecHistNames = {
0092 "cluster_et",
0093 "e11_to_e33",
0094 "e33_to_e55",
0095 "e55_to_e77",
0096 "e32_to_e35",
0097 "weta",
0098 "wphi",
0099 "weta_cogx",
0100 "wphi_cogx",
0101 "detamax",
0102 "dphimax",
0103 "mean_time",
0104 "iso04_emcal",
0105 "weta_vs_et",
0106 "wphi_vs_et"};
0107
0108 for (auto &histName : vecHistNames)
0109 {
0110 histName.insert(0, "h_" + smallModuleName + "_");
0111 if (!m_histtag.empty())
0112 {
0113 histName.append("_" + m_histtag);
0114 }
0115 }
0116
0117 h_cluster_et = new TH1F(vecHistNames[0].data(), "", 120, 0, 30);
0118 h_cluster_et->GetXaxis()->SetTitle("E_{T} [GeV]");
0119
0120 h_e11oe33 = new TH1F(vecHistNames[1].data(), "", 26, -0.02, 1.02);
0121 h_e11oe33->GetXaxis()->SetTitle("e11/e33");
0122
0123 h_e33oe55 = new TH1F(vecHistNames[2].data(), "", 26, -0.02, 1.02);
0124 h_e33oe55->GetXaxis()->SetTitle("e33/e55");
0125
0126 h_e55oe77 = new TH1F(vecHistNames[3].data(), "", 26, -0.02, 1.02);
0127 h_e55oe77->GetXaxis()->SetTitle("e55/e77");
0128
0129 h_e32oe35 = new TH1F(vecHistNames[4].data(), "", 26, -0.02, 1.02);
0130 h_e32oe35->GetXaxis()->SetTitle("e32/e35");
0131
0132 h_weta = new TH1F(vecHistNames[5].data(), "", 120, 0, 2);
0133 h_weta->GetXaxis()->SetTitle("w#eta");
0134
0135 h_wphi = new TH1F(vecHistNames[6].data(), "", 120, 0, 2);
0136 h_wphi->GetXaxis()->SetTitle("w#phi");
0137
0138 h_weta_cogx = new TH1F(vecHistNames[7].data(), "", 50, 0, 2);
0139 h_weta_cogx->GetXaxis()->SetTitle("w#eta_cogx");
0140
0141 h_wphi_cogx = new TH1F(vecHistNames[8].data(), "", 50, 0, 2);
0142 h_wphi_cogx->GetXaxis()->SetTitle("w#phi_cogx");
0143
0144 h_detamax = new TH1F(vecHistNames[9].data(), "", 10, -0.5, 9.5);
0145 h_detamax->GetXaxis()->SetTitle("detamax");
0146
0147 h_dphimax = new TH1F(vecHistNames[10].data(), "", 20, -0.5, 19.5);
0148 h_dphimax->GetXaxis()->SetTitle("dphimax");
0149
0150 h_mean_time = new TH1F(vecHistNames[11].data(), "", 200, -20, 20);
0151 h_mean_time->GetXaxis()->SetTitle("cluster mean time");
0152
0153 h_iso04_emcal = new TH1F(vecHistNames[12].data(), "", 200, -10, 40);
0154 h_iso04_emcal->GetXaxis()->SetTitle("iso_{0.4}^{EMCal} [GeV]");
0155
0156 h_weta_vs_et = new TH2F(vecHistNames[13].data(), "", 120, 0, 30, 120, 0, 6);
0157 h_weta_vs_et->GetXaxis()->SetTitle("E_{T} [GeV]");
0158 h_weta_vs_et->GetYaxis()->SetTitle("w#eta");
0159
0160 h_wphi_vs_et = new TH2F(vecHistNames[14].data(), "", 120, 0, 30, 120, 0, 6);
0161 h_wphi_vs_et->GetXaxis()->SetTitle("E_{T} [GeV]");
0162 h_wphi_vs_et->GetYaxis()->SetTitle("w#phi");
0163
0164
0165 m_manager->registerHisto(h_cluster_et);
0166 m_manager->registerHisto(h_e11oe33);
0167 m_manager->registerHisto(h_e33oe55);
0168 m_manager->registerHisto(h_e55oe77);
0169 m_manager->registerHisto(h_e32oe35);
0170 m_manager->registerHisto(h_weta);
0171 m_manager->registerHisto(h_wphi);
0172 m_manager->registerHisto(h_weta_cogx);
0173 m_manager->registerHisto(h_wphi_cogx);
0174 m_manager->registerHisto(h_detamax);
0175 m_manager->registerHisto(h_dphimax);
0176 m_manager->registerHisto(h_mean_time);
0177 m_manager->registerHisto(h_iso04_emcal);
0178 m_manager->registerHisto(h_weta_vs_et);
0179 m_manager->registerHisto(h_wphi_vs_et);
0180
0181 return Fun4AllReturnCodes::EVENT_OK;
0182 }
0183
0184 int EMCalShowerShapes::InitRun(PHCompositeNode* topNode)
0185 {
0186 if (!LoadEMCalNodes(topNode))
0187 {
0188 return Fun4AllReturnCodes::ABORTRUN;
0189 }
0190 return Fun4AllReturnCodes::EVENT_OK;
0191 }
0192
0193 bool EMCalShowerShapes::LoadEMCalNodes(PHCompositeNode *topNode)
0194 {
0195 m_emc_tower_container = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0196 m_geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0197
0198 const bool have_nodes = (m_emc_tower_container && m_geomEM);
0199 if (!have_nodes && !m_reportedMissingCaloNodes)
0200 {
0201 std::cout << PHWHERE << "EMCalShowerShapes::LoadEMCalNodes - missing TOWERINFO_CALIB_CEMC or TOWERGEOM_CEMC" << std::endl;
0202 m_reportedMissingCaloNodes = true;
0203 }
0204 return have_nodes;
0205 }
0206
0207 float EMCalShowerShapes::GetVertexZ(PHCompositeNode *topNode) const
0208 {
0209 MbdVertexMap* vertexmap = findNode::getClass<MbdVertexMap>(topNode, "MbdVertexMap");
0210 if (!vertexmap || vertexmap->empty())
0211 {
0212 return 0.0F;
0213 }
0214
0215 MbdVertex* vtx = vertexmap->begin()->second;
0216 if (!vtx)
0217 {
0218 return 0.0F;
0219 }
0220
0221 return vtx->get_z();
0222 }
0223
0224 int EMCalShowerShapes::process_event(PHCompositeNode *topNode)
0225 {
0226 RawClusterContainer* clusterContainer = findNode::getClass<RawClusterContainer>(topNode, m_inputnode);
0227 if (!clusterContainer)
0228 {
0229 if (!m_reportedMissingClusterNode)
0230 {
0231 std::cout << PHWHERE << "EMCalShowerShapes::process_event - missing node " << m_inputnode << std::endl;
0232 m_reportedMissingClusterNode = true;
0233 }
0234 return Fun4AllReturnCodes::ABORTRUN;
0235 }
0236
0237 if (!LoadEMCalNodes(topNode))
0238 {
0239 return Fun4AllReturnCodes::ABORTRUN;
0240 }
0241
0242 if (m_doTrgSelect)
0243 {
0244 m_analyzer->decodeTriggers(topNode);
0245 if (!JetQADefs::DidTriggerFire(m_trgToSelect, m_analyzer))
0246 {
0247 return Fun4AllReturnCodes::EVENT_OK;
0248 }
0249 }
0250
0251 const float vertex_z = GetVertexZ(topNode);
0252 if (m_doMbdZvtxCut && std::abs(vertex_z) >= m_mbdZvtxMax)
0253 {
0254 return Fun4AllReturnCodes::EVENT_OK;
0255 }
0256
0257 const CLHEP::Hep3Vector vertex_vec(0, 0, vertex_z);
0258
0259 RawClusterContainer::ConstRange clusters = clusterContainer->getClusters();
0260 for (auto iter = clusters.first; iter != clusters.second; ++iter)
0261 {
0262 RawCluster* cluster = iter->second;
0263 if (!cluster)
0264 {
0265 continue;
0266 }
0267
0268 const float eta = RawClusterUtility::GetPseudorapidity(*cluster, vertex_vec);
0269 if (m_doClusterEtaCut && std::abs(eta) >= m_clusterEtaMax)
0270 {
0271 continue;
0272 }
0273
0274 const float phi = RawClusterUtility::GetAzimuthAngle(*cluster, vertex_vec);
0275 const float et = cluster->get_energy() / std::cosh(eta);
0276 if (m_doClusterETCut && et < m_clusterETMin)
0277 {
0278 continue;
0279 }
0280
0281 ShowerShapeData data;
0282 if (!CalculateShowerShapes(cluster, eta, phi, et, vertex_z, data))
0283 {
0284 continue;
0285 }
0286
0287 h_cluster_et->Fill(et);
0288 if (data.e33 > 0) {
0289 h_e11oe33->Fill(data.e11 / data.e33);
0290 }
0291 if (data.e55 > 0) {
0292 h_e33oe55->Fill(data.e33 / data.e55);
0293 }
0294 if (data.e77 > 0) {
0295 h_e55oe77->Fill(data.e55 / data.e77);
0296 }
0297 if (data.e35 > 0) {
0298 h_e32oe35->Fill(data.e32 / data.e35);
0299 }
0300 h_weta->Fill(data.weta);
0301 h_wphi->Fill(data.wphi);
0302 h_weta_cogx->Fill(data.weta_cogx);
0303 h_wphi_cogx->Fill(data.wphi_cogx);
0304 h_detamax->Fill(data.detamax);
0305 h_dphimax->Fill(data.dphimax);
0306 h_mean_time->Fill(data.mean_time);
0307 h_iso04_emcal->Fill(data.iso04_emcal);
0308 h_weta_vs_et->Fill(et, data.weta);
0309 h_wphi_vs_et->Fill(et, data.wphi);
0310
0311 }
0312
0313 return Fun4AllReturnCodes::EVENT_OK;
0314 }
0315
0316 bool EMCalShowerShapes::CalculateShowerShapes(RawCluster* cluster, float cluster_eta, float cluster_phi, float cluster_et, float vertex_z, ShowerShapeData& data) const
0317 {
0318 std::vector<float> showershape = cluster->get_shower_shapes(m_shape_min_tower_E);
0319 if (showershape.empty())
0320 {
0321 return false;
0322 }
0323
0324 const std::pair<int, int> leadtowerindex = cluster->get_lead_tower();
0325 const int lead_ieta = leadtowerindex.first;
0326 const int lead_iphi = leadtowerindex.second;
0327
0328 const float avg_eta = showershape[4] + 0.5F;
0329 const float avg_phi = showershape[5] + 0.5F;
0330 const int maxieta = std::floor(avg_eta);
0331 const int maxiphi = std::floor(avg_phi);
0332
0333 int detamax = 0;
0334 int dphimax = 0;
0335 float clusteravgtime = 0.0F;
0336 float cluster_total_e = 0.0F;
0337 const RawCluster::TowerMap& tower_map = cluster->get_towermap();
0338 std::set<unsigned int> towers_in_cluster;
0339 for (auto tower_iter : tower_map)
0340 {
0341 RawTowerDefs::keytype tower_key = tower_iter.first;
0342 const int ieta = RawTowerDefs::decode_index1(tower_key);
0343 const int iphi = RawTowerDefs::decode_index2(tower_key);
0344
0345 const unsigned int towerinfokey = TowerInfoDefs::encode_emcal(ieta, iphi);
0346 towers_in_cluster.insert(towerinfokey);
0347 TowerInfo* towerinfo = m_emc_tower_container->get_tower_at_key(towerinfokey);
0348 if (towerinfo)
0349 {
0350 clusteravgtime += towerinfo->get_time() * towerinfo->get_energy();
0351 cluster_total_e += towerinfo->get_energy();
0352 }
0353
0354 int totalphibins = 256;
0355 auto dphiwrap = [totalphibins](int towerphi, int maxiphi_arg)
0356 {
0357 int idphi = towerphi - maxiphi_arg;
0358 if (idphi > totalphibins / 2)
0359 {
0360 idphi -= totalphibins;
0361 }
0362 if (idphi < -totalphibins / 2)
0363 {
0364 idphi += totalphibins;
0365 }
0366 return idphi;
0367 };
0368
0369 const int deta = ieta - lead_ieta;
0370 const int dphi_val = dphiwrap(iphi, lead_iphi);
0371 detamax = std::max(std::abs(deta), detamax);
0372 dphimax = std::max(std::abs(dphi_val), dphimax);
0373 }
0374
0375 if (cluster_total_e > 0)
0376 {
0377 clusteravgtime /= cluster_total_e;
0378 }
0379 else
0380 {
0381 std::cout << "cluster_total_e is 0(this should not happen!!!), setting clusteravgtime to NaN" << std::endl;
0382 clusteravgtime = std::numeric_limits<float>::quiet_NaN();
0383 }
0384
0385 float E77[7][7] = {{0.0F}};
0386 int E77_ownership[7][7] = {{0}};
0387
0388 for (int ieta = maxieta - 3; ieta < maxieta + 4; ++ieta)
0389 {
0390 for (int iphi = maxiphi - 3; iphi < maxiphi + 4; ++iphi)
0391 {
0392 if (ieta < 0 || ieta > 95)
0393 {
0394 E77[ieta - maxieta + 3][iphi - maxiphi + 3] = 0.0F;
0395 E77_ownership[ieta - maxieta + 3][iphi - maxiphi + 3] = 0;
0396 continue;
0397 }
0398
0399 int temp_ieta = ieta;
0400 int temp_iphi = iphi;
0401 shift_tower_index(temp_ieta, temp_iphi, 96, 256);
0402 if (temp_ieta < 0)
0403 {
0404 continue;
0405 }
0406
0407 const unsigned int towerinfokey = TowerInfoDefs::encode_emcal(temp_ieta, temp_iphi);
0408
0409 if (towers_in_cluster.contains(towerinfokey))
0410 {
0411 E77_ownership[ieta - maxieta + 3][iphi - maxiphi + 3] = 1;
0412 }
0413
0414 TowerInfo* towerinfo = m_emc_tower_container->get_tower_at_key(towerinfokey);
0415 if (towerinfo && towerinfo->get_isGood())
0416 {
0417 const float energy = towerinfo->get_energy();
0418 if (energy > m_shape_min_tower_E)
0419 {
0420 E77[ieta - maxieta + 3][iphi - maxiphi + 3] = energy;
0421 }
0422 }
0423 }
0424 }
0425
0426 float e11 = E77[3][3];
0427 float e32 = 0.0F;
0428 float e33 = 0.0F;
0429 float e35 = 0.0F;
0430 float e55 = 0.0F;
0431 float e77 = 0.0F;
0432 float weta = 0.0F;
0433 float wphi = 0.0F;
0434 float weta_cogx = 0.0F;
0435 float wphi_cogx = 0.0F;
0436 float Eetaphi = 0.0F;
0437
0438 const float shift_eta = avg_eta - std::floor(avg_eta) - 0.5F;
0439 const float shift_phi = avg_phi - std::floor(avg_phi) - 0.5F;
0440 const float cog_eta = 3 + shift_eta;
0441 const float cog_phi = 3 + shift_phi;
0442 const int signphi = (avg_phi - std::floor(avg_phi)) > 0.5 ? 1 : -1;
0443
0444 for (int i = 0; i < 7; ++i)
0445 {
0446 for (int j = 0; j < 7; ++j)
0447 {
0448 const int di = std::abs(i - 3);
0449 const int dj = std::abs(j - 3);
0450 const float di_float = i - cog_eta;
0451 const float dj_float = j - cog_phi;
0452
0453 if (E77_ownership[i][j] == 1)
0454 {
0455 weta += E77[i][j] * di * di;
0456 wphi += E77[i][j] * dj * dj;
0457 Eetaphi += E77[i][j];
0458 if (i != 3 || j != 3)
0459 {
0460 weta_cogx += E77[i][j] * di_float * di_float;
0461 wphi_cogx += E77[i][j] * dj_float * dj_float;
0462 }
0463 }
0464
0465 e77 += E77[i][j];
0466 if (di <= 1 && (dj == 0 || j == (3 + signphi)))
0467 {
0468 e32 += E77[i][j];
0469 }
0470 if (di <= 1 && dj <= 1)
0471 {
0472 e33 += E77[i][j];
0473 }
0474 if (di <= 1 && dj <= 2)
0475 {
0476 e35 += E77[i][j];
0477 }
0478 if (di <= 2 && dj <= 2)
0479 {
0480 e55 += E77[i][j];
0481 }
0482 }
0483 }
0484
0485 if (Eetaphi > 0)
0486 {
0487 weta /= Eetaphi;
0488 wphi /= Eetaphi;
0489 weta_cogx /= Eetaphi;
0490 wphi_cogx /= Eetaphi;
0491 }
0492
0493
0494
0495
0496
0497
0498
0499
0500 data.e11 = e11;
0501 data.e33 = e33;
0502 data.e32 = e32;
0503 data.e35 = e35;
0504 data.e55 = e55;
0505 data.e77 = e77;
0506 data.weta = weta;
0507 data.wphi = wphi;
0508 data.weta_cogx = weta_cogx;
0509 data.wphi_cogx = wphi_cogx;
0510 data.detamax = detamax;
0511 data.dphimax = dphimax;
0512 data.mean_time = clusteravgtime;
0513 data.iso04_emcal = CalculateLayerET(cluster_eta, cluster_phi, 0.4F, m_emc_tower_container, m_geomEM, vertex_z) - cluster_et;
0514
0515 return true;
0516 }
0517
0518 double EMCalShowerShapes::GetTowerEta(RawTowerGeom* tower_geom, double vx, double vy, double vz) const
0519 {
0520 if (!tower_geom)
0521 {
0522 return -9999;
0523 }
0524 if (vx == 0 && vy == 0 && vz == 0)
0525 {
0526 return tower_geom->get_eta();
0527 }
0528
0529 const double radius = std::sqrt((tower_geom->get_center_x() - vx) * (tower_geom->get_center_x() - vx) +
0530 (tower_geom->get_center_y() - vy) * (tower_geom->get_center_y() - vy));
0531 const double theta = std::atan2(radius, tower_geom->get_center_z() - vz);
0532 return -std::log(std::tan(theta / 2.));
0533 }
0534
0535 double EMCalShowerShapes::DeltaR(double eta1, double phi1, double eta2, double phi2) const
0536 {
0537 double dphi = phi1 - phi2;
0538 while (dphi > M_PI)
0539 {
0540 dphi -= 2 * M_PI;
0541 }
0542 while (dphi <= -M_PI)
0543 {
0544 dphi += 2 * M_PI;
0545 }
0546 return std::sqrt(std::pow(eta1 - eta2, 2) + std::pow(dphi, 2));
0547 }
0548
0549 float EMCalShowerShapes::CalculateLayerET(float seed_eta, float seed_phi, float radius, TowerInfoContainer* towerContainer, RawTowerGeomContainer* geomContainer, float vertex_z) const
0550 {
0551 if (!towerContainer || !geomContainer)
0552 {
0553 return std::numeric_limits<float>::quiet_NaN();
0554 }
0555
0556 float layer_et = 0.0F;
0557 const unsigned int ntowers = towerContainer->size();
0558 for (unsigned int channel = 0; channel < ntowers; ++channel)
0559 {
0560 TowerInfo* tower = towerContainer->get_tower_at_channel(channel);
0561 if (!tower || !tower->get_isGood())
0562 {
0563 continue;
0564 }
0565
0566 const unsigned int towerkey = towerContainer->encode_key(channel);
0567 const int ieta = towerContainer->getTowerEtaBin(towerkey);
0568 const int iphi = towerContainer->getTowerPhiBin(towerkey);
0569
0570 const RawTowerDefs::keytype geom_key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
0571 RawTowerGeom* tower_geom = geomContainer->get_tower_geometry(geom_key);
0572 if (!tower_geom)
0573 {
0574 continue;
0575 }
0576
0577 const double tower_eta = GetTowerEta(tower_geom, 0, 0, vertex_z);
0578 const double tower_phi = tower_geom->get_phi();
0579 if (DeltaR(seed_eta, seed_phi, tower_eta, tower_phi) >= radius)
0580 {
0581 continue;
0582 }
0583
0584 const float energy = tower->get_energy();
0585 if (energy <= m_shape_min_tower_E)
0586 {
0587 continue;
0588 }
0589
0590 layer_et += energy / std::cosh(tower_eta);
0591 }
0592
0593 return layer_et;
0594 }
0595
0596
0597
0598
0599
0600
0601
0602
0603
0604
0605
0606
0607
0608
0609
0610
0611
0612
0613
0614
0615
0616 void EMCalShowerShapes::Print(const std::string &what) const
0617 {
0618 std::cout << "EMCalShowerShapes::Print(" << what << ")" << std::endl;
0619 }
