File indexing completed on 2025-08-06 08:17:36
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "ClusterIso.h"
0011
0012 #include <calobase/RawCluster.h>
0013 #include <calobase/RawClusterContainer.h>
0014 #include <calobase/RawClusterUtility.h>
0015 #include <calobase/RawTowerGeom.h>
0016 #include <calobase/RawTowerGeomContainer.h>
0017 #include <calobase/TowerInfo.h>
0018 #include <calobase/TowerInfoContainer.h>
0019
0020 #include <globalvertex/GlobalVertex.h>
0021 #include <globalvertex/GlobalVertexMap.h>
0022
0023 #include <fun4all/Fun4AllBase.h> // for Fun4AllBase::VERBOSITY_MORE
0024 #include <fun4all/SubsysReco.h>
0025
0026 #include <phool/getClass.h>
0027
0028 #include <CLHEP/Vector/ThreeVector.h>
0029
0030 #include <iostream>
0031 #include <map>
0032 #include <utility>
0033
0034
0035
0036
0037
0038
0039
0040 double ClusterIso::getTowerEta(RawTowerGeom *tower_geom, double vx, double vy, double vz)
0041 {
0042 float r;
0043 if (vx == 0 && vy == 0 && vz == 0)
0044 {
0045 r = tower_geom->get_eta();
0046 }
0047 else
0048 {
0049 double radius = sqrt((tower_geom->get_center_x() - vx) * (tower_geom->get_center_x() - vx) + (tower_geom->get_center_y() - vy) * (tower_geom->get_center_y() - vy));
0050 double theta = atan2(radius, tower_geom->get_center_z() - vz);
0051 r = -log(tan(theta / 2.));
0052 }
0053 return r;
0054 }
0055
0056
0057
0058
0059
0060 ClusterIso::ClusterIso(const std::string &kname, float eTCut = 0.0, int coneSize = 3, bool do_subtracted = true, bool do_unsubtracted = true)
0061 : SubsysReco(kname)
0062 , m_do_subtracted(do_subtracted)
0063 , m_do_unsubtracted(do_unsubtracted)
0064 , m_use_towerinfo(true)
0065 {
0066 if (Verbosity() >= VERBOSITY_SOME)
0067 {
0068 std::cout << Name() << "::ClusterIso constructed" << '\n';
0069 }
0070 if (coneSize == 0 && Verbosity() >= VERBOSITY_QUIET)
0071 {
0072 std::cout << "WARNING in " << Name() << "ClusterIso:: cone size is zero" << '\n';
0073 }
0074 m_vx = m_vy = m_vz = 0;
0075 setConeSize(coneSize);
0076 seteTCut(eTCut);
0077 if (Verbosity() >= VERBOSITY_EVEN_MORE)
0078 {
0079 std::cout << Name() << "::ClusterIso::m_coneSize is:" << m_coneSize << '\n';
0080 std::cout << Name() << "::ClusterIso::m_eTCut is:" << m_eTCut << '\n';
0081 }
0082 if (!do_subtracted && !do_unsubtracted && Verbosity() >= VERBOSITY_QUIET)
0083 {
0084 std::cout << "WARNING in " << Name() << "ClusterIso:: all processes turned off doing nothing" << '\n';
0085 }
0086 }
0087
0088 int ClusterIso::Init(PHCompositeNode * )
0089 {
0090 return 0;
0091 }
0092
0093
0094
0095
0096 void ClusterIso::seteTCut(float eTCut)
0097 {
0098 this->m_eTCut = eTCut;
0099 }
0100
0101
0102
0103
0104 void ClusterIso::setConeSize(int coneSize)
0105 {
0106 this->m_coneSize = coneSize / 10.0;
0107 }
0108
0109
0110
0111
0112 float ClusterIso::geteTCut()
0113 {
0114 return m_eTCut;
0115 }
0116
0117
0118
0119
0120 int ClusterIso::getConeSize()
0121 {
0122 return (int) m_coneSize * 10;
0123 }
0124
0125
0126
0127
0128 CLHEP::Hep3Vector ClusterIso::getVertex()
0129 {
0130 return CLHEP::Hep3Vector(m_vx, m_vy, m_vz);
0131 }
0132
0133
0134
0135
0136
0137
0138
0139 int ClusterIso::process_event(PHCompositeNode *topNode)
0140 {
0141 if (Verbosity() >= VERBOSITY_MORE)
0142 {
0143 std::cout << Name() << "::ClusterIso::process_event" << '\n';
0144 }
0145
0146
0147
0148
0149
0150
0151
0152
0153 std::string RawCemcClusterNodeName = "CLUSTER_CEMC";
0154 if (m_use_towerinfo)
0155 {
0156 RawCemcClusterNodeName = m_cluster_node_name;
0157 }
0158
0159 if (m_do_subtracted)
0160 {
0161 {
0162 if (Verbosity() >= VERBOSITY_EVEN_MORE)
0163 {
0164 std::cout << Name() << "::ClusterIso starting subtracted calculation" << '\n';
0165 }
0166
0167 TowerInfoContainer *towersEM3old = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
0168 if (towersEM3old == nullptr)
0169 {
0170 m_do_subtracted = false;
0171 if (Verbosity() >= VERBOSITY_SOME)
0172 {
0173 std::cout << "In " << Name() << "::ClusterIso WARNING substracted towers do not exist subtracted isolation cannot be preformed \n";
0174 }
0175 }
0176 if (Verbosity() >= VERBOSITY_MORE)
0177 {
0178 std::cout << Name() << "::ClusterIso::process_event: " << towersEM3old->size() << " TOWERINFO_CALIB_CEMC_RETOWER_SUB1 towers" << '\n';
0179 }
0180
0181
0182 TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
0183 if (Verbosity() >= VERBOSITY_MORE)
0184 {
0185 std::cout << Name() << "::ClusterIso::process_event: " << towersIH3->size() << " TOWERINFO_CALIB_HCALIN_SUB1 towers" << '\n';
0186 }
0187
0188
0189 TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
0190 if (Verbosity() >= VERBOSITY_MORE)
0191 {
0192 std::cout << Name() << "::ClusterIso::process_event: " << towersOH3->size() << " TOWERINFO_CALIB_HCALOUT_SUB1 towers" << std::endl;
0193 }
0194
0195
0196 RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0197 RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0198 RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0199
0200 {
0201 RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, RawCemcClusterNodeName.c_str());
0202 RawClusterContainer::ConstRange begin_end = clusters->getClusters();
0203 RawClusterContainer::ConstIterator rtiter;
0204 if (Verbosity() >= VERBOSITY_SOME)
0205 {
0206 std::cout << Name() << "::ClusterIso sees " << clusters->size() << " clusters " << '\n';
0207 }
0208
0209
0210 GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0211 m_vx = m_vy = m_vz = 0;
0212 if (vertexmap && !vertexmap->empty())
0213 {
0214 GlobalVertex *vertex = (vertexmap->begin()->second);
0215 m_vx = vertex->get_x();
0216 m_vy = vertex->get_y();
0217 m_vz = vertex->get_z();
0218 if (Verbosity() >= VERBOSITY_SOME)
0219 {
0220 std::cout << Name() << "::ClusterIso Event Vertex Calculated at x:" << m_vx << " y:" << m_vy << " z:" << m_vz << '\n';
0221 }
0222 }
0223
0224 for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0225 {
0226 RawCluster *cluster = rtiter->second;
0227
0228 CLHEP::Hep3Vector vertex(m_vx, m_vy, m_vz);
0229 CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
0230 double cluster_energy = E_vec_cluster.mag();
0231 double cluster_eta = E_vec_cluster.pseudoRapidity();
0232 double cluster_phi = E_vec_cluster.phi();
0233 double et = cluster_energy / cosh(cluster_eta);
0234 double isoEt = 0;
0235
0236 if (et < m_eTCut)
0237 {
0238 continue;
0239 }
0240
0241
0242 {
0243 unsigned int ntowers = towersEM3old->size();
0244 for (unsigned int channel = 0; channel < ntowers; channel++)
0245 {
0246 TowerInfo *tower = towersEM3old->get_tower_at_channel(channel);
0247 if (!IsAcceptableTower(tower))
0248 {
0249 continue;
0250 }
0251 unsigned int towerkey = towersEM3old->encode_key(channel);
0252 int ieta = towersEM3old->getTowerEtaBin(towerkey);
0253 int iphi = towersEM3old->getTowerPhiBin(towerkey);
0254
0255 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0256
0257 RawTowerGeom *tower_geom = geomEM->get_tower_geometry(key);
0258 double this_phi = tower_geom->get_phi();
0259 double this_eta = tower_geom->get_eta();
0260 if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0261 {
0262 isoEt += tower->get_energy() / cosh(this_eta);
0263 }
0264 }
0265 }
0266
0267
0268 {
0269 unsigned int ntowers = towersIH3->size();
0270 for (unsigned int channel = 0; channel < ntowers; channel++)
0271 {
0272 TowerInfo *tower = towersIH3->get_tower_at_channel(channel);
0273 if (!IsAcceptableTower(tower))
0274 {
0275 continue;
0276 }
0277 unsigned int towerkey = towersIH3->encode_key(channel);
0278 int ieta = towersIH3->getTowerEtaBin(towerkey);
0279 int iphi = towersIH3->getTowerPhiBin(towerkey);
0280 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0281 RawTowerGeom *tower_geom = geomIH->get_tower_geometry(key);
0282 double this_phi = tower_geom->get_phi();
0283 double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
0284 if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0285 {
0286 isoEt += tower->get_energy() / cosh(this_eta);
0287 }
0288 }
0289 }
0290
0291
0292 {
0293 unsigned int ntowers = towersOH3->size();
0294 for (unsigned int channel = 0; channel < ntowers; channel++)
0295 {
0296 TowerInfo *tower = towersOH3->get_tower_at_channel(channel);
0297 if (!IsAcceptableTower(tower))
0298 {
0299 continue;
0300 }
0301 unsigned int towerkey = towersOH3->encode_key(channel);
0302 int ieta = towersOH3->getTowerEtaBin(towerkey);
0303 int iphi = towersOH3->getTowerPhiBin(towerkey);
0304 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0305 RawTowerGeom *tower_geom = geomOH->get_tower_geometry(key);
0306 double this_phi = tower_geom->get_phi();
0307 double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
0308 if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0309 {
0310 isoEt += tower->get_energy() / cosh(this_eta);
0311 }
0312 }
0313 }
0314
0315 isoEt -= et;
0316 if (Verbosity() >= VERBOSITY_EVEN_MORE)
0317 {
0318 std::cout << Name() << "::ClusterIso iso_et for ";
0319 cluster->identify();
0320 std::cout << "=" << isoEt << '\n';
0321 }
0322 cluster->set_et_iso(isoEt, (int) 10 * m_coneSize, true, true);
0323 }
0324 }
0325 }
0326 }
0327 if (m_do_unsubtracted)
0328 {
0329
0330
0331
0332 if (Verbosity() >= VERBOSITY_EVEN_MORE)
0333 {
0334 std::cout << Name() << "::ClusterIso starting unsubtracted calculation" << '\n';
0335 }
0336 {
0337
0338 TowerInfoContainer *towersEM3old = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0339 if (Verbosity() >= VERBOSITY_MORE)
0340 {
0341 std::cout << "ClusterIso::process_event: " << towersEM3old->size() << " TOWERINFO_CALIB_CEMC towers" << '\n';
0342 }
0343
0344
0345 TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0346 if (Verbosity() >= VERBOSITY_MORE)
0347 {
0348 std::cout << "ClusterIso::process_event: " << towersIH3->size() << " TOWERINFO_CALIB_HCALIN towers" << '\n';
0349 }
0350
0351
0352 TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0353 if (Verbosity() >= VERBOSITY_MORE)
0354 {
0355 std::cout << "ClusterIso::process_event: " << towersOH3->size() << " TOWERINFO_CALIB_HCALOUT towers" << std::endl;
0356 }
0357
0358
0359 RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0360 RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0361 RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0362
0363 {
0364 RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, RawCemcClusterNodeName.c_str());
0365 RawClusterContainer::ConstRange begin_end = clusters->getClusters();
0366 RawClusterContainer::ConstIterator rtiter;
0367 if (Verbosity() >= VERBOSITY_SOME)
0368 {
0369 std::cout << " ClusterIso sees " << clusters->size() << " clusters " << '\n';
0370 }
0371
0372
0373 GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0374 m_vx = m_vy = m_vz = 0;
0375 if (vertexmap && !vertexmap->empty())
0376 {
0377 GlobalVertex *vertex = (vertexmap->begin()->second);
0378 m_vx = vertex->get_x();
0379 m_vy = vertex->get_y();
0380 m_vz = vertex->get_z();
0381 if (Verbosity() >= VERBOSITY_SOME)
0382 {
0383 std::cout << Name() << "ClusterIso Event Vertex Calculated at x:" << m_vx << " y:" << m_vy << " z:" << m_vz << '\n';
0384 }
0385 }
0386
0387 for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0388 {
0389 RawCluster *cluster = rtiter->second;
0390
0391 CLHEP::Hep3Vector vertex(m_vx, m_vy, m_vz);
0392 CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
0393 double cluster_energy = E_vec_cluster.mag();
0394 double cluster_eta = E_vec_cluster.pseudoRapidity();
0395 double cluster_phi = E_vec_cluster.phi();
0396 double et = cluster_energy / cosh(cluster_eta);
0397 double isoEt = 0;
0398 if (Verbosity() >= VERBOSITY_MAX)
0399 {
0400 std::cout << Name() << "::ClusterIso processing";
0401 cluster->identify();
0402 std::cout << '\n';
0403 }
0404 if (et < m_eTCut)
0405 {
0406 if (Verbosity() >= VERBOSITY_MAX)
0407 {
0408 std::cout << "\t does not pass eT cut" << '\n';
0409 }
0410 continue;
0411 }
0412
0413
0414 {
0415 unsigned int ntowers = towersEM3old->size();
0416 for (unsigned int channel = 0; channel < ntowers; channel++)
0417 {
0418 TowerInfo *tower = towersEM3old->get_tower_at_channel(channel);
0419 if (!IsAcceptableTower(tower))
0420 {
0421 continue;
0422 }
0423 unsigned int towerkey = towersEM3old->encode_key(channel);
0424 int ieta = towersEM3old->getTowerEtaBin(towerkey);
0425 int iphi = towersEM3old->getTowerPhiBin(towerkey);
0426 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
0427 RawTowerGeom *tower_geom = geomEM->get_tower_geometry(key);
0428 double this_phi = tower_geom->get_phi();
0429 double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
0430 if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0431 {
0432 isoEt += tower->get_energy() / cosh(this_eta);
0433 }
0434 }
0435 }
0436 if (Verbosity() >= VERBOSITY_MAX)
0437 {
0438 std::cout << "\t after EMCal isoEt:" << isoEt << '\n';
0439 }
0440
0441 {
0442 unsigned int ntowers = towersIH3->size();
0443 for (unsigned int channel = 0; channel < ntowers; channel++)
0444 {
0445 TowerInfo *tower = towersIH3->get_tower_at_channel(channel);
0446 if (!IsAcceptableTower(tower))
0447 {
0448 continue;
0449 }
0450 unsigned int towerkey = towersIH3->encode_key(channel);
0451 int ieta = towersIH3->getTowerEtaBin(towerkey);
0452 int iphi = towersIH3->getTowerPhiBin(towerkey);
0453 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0454 RawTowerGeom *tower_geom = geomIH->get_tower_geometry(key);
0455 double this_phi = tower_geom->get_phi();
0456 double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
0457 if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0458 {
0459 isoEt += tower->get_energy() / cosh(this_eta);
0460 }
0461 }
0462 }
0463 if (Verbosity() >= VERBOSITY_MAX)
0464 {
0465 std::cout << "\t after innerHCal isoEt:" << isoEt << '\n';
0466 }
0467
0468 {
0469 unsigned int ntowers = towersOH3->size();
0470 for (unsigned int channel = 0; channel < ntowers; channel++)
0471 {
0472 TowerInfo *tower = towersOH3->get_tower_at_channel(channel);
0473 if (!IsAcceptableTower(tower))
0474 {
0475 continue;
0476 }
0477 unsigned int towerkey = towersOH3->encode_key(channel);
0478 int ieta = towersOH3->getTowerEtaBin(towerkey);
0479 int iphi = towersOH3->getTowerPhiBin(towerkey);
0480 const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0481 RawTowerGeom *tower_geom = geomOH->get_tower_geometry(key);
0482 double this_phi = tower_geom->get_phi();
0483 double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
0484 if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0485 {
0486 isoEt += tower->get_energy() / cosh(this_eta);
0487 }
0488 }
0489 }
0490 if (Verbosity() >= VERBOSITY_MAX)
0491 {
0492 std::cout << "\t after outerHCal isoEt:" << isoEt << '\n';
0493 }
0494 isoEt -= et;
0495 if (Verbosity() >= VERBOSITY_EVEN_MORE)
0496 {
0497 std::cout << Name() << "::ClusterIso iso_et for ";
0498 cluster->identify();
0499 std::cout << "=" << isoEt << '\n';
0500 }
0501 cluster->set_et_iso(isoEt, (int) 10 * m_coneSize, false, true);
0502 }
0503 }
0504 }
0505 }
0506 return 0;
0507 }
0508
0509 int ClusterIso::End(PHCompositeNode * )
0510 {
0511 return 0;
0512 }
0513
0514 bool ClusterIso::IsAcceptableTower(TowerInfo *tower)
0515 {
0516 if (!tower->get_isGood())
0517 {
0518 return false;
0519 }
0520 return true;
0521 }