Back to home page

sPhenix code displayed by LXR

 
 

    


File indexing completed on 2025-08-06 08:17:36

0001 /*!
0002  * \file ClusterIso.cc
0003  * \brief
0004  * \author Francesco Vassalli <Francesco.Vassalli@colorado.edu>
0005  * \author Chase Smith <chsm5267@colorado.edu>
0006  * \version $Revision: 2 $
0007  * \date $Date: July 17th, 2018 $
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 /** \Brief Function to get correct tower eta
0035  *
0036  * Each calorimeter tower's eta is calculated using the vertex (0,0,0)
0037  * which is incorrect in many collisions. This function
0038  * uses geometry to find a given tower's eta using the correct vertex.
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  * Contructor takes the argument of the class name, the minimum eT of the clusters which defaults to 0,
0058  * and the isolation cone size which defaults to 0.3.
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 * /*topNode*/)
0089 {
0090   return 0;
0091 }
0092 
0093 /**
0094  * Set the minimum transverse energy required for a cluster to have its isolation calculated
0095  */
0096 void ClusterIso::seteTCut(float eTCut)
0097 {
0098   this->m_eTCut = eTCut;
0099 }
0100 
0101 /**
0102  * Set the size of isolation cone as integer multiple of 0.1, (i.e. 3 will use an R=0.3 cone)
0103  */
0104 void ClusterIso::setConeSize(int coneSize)
0105 {
0106   this->m_coneSize = coneSize / 10.0;
0107 }
0108 
0109 /**
0110  * Returns the minimum transverse energy required for a cluster to have its isolation calculated
0111  */
0112 /*const*/ float ClusterIso::geteTCut()
0113 {
0114   return m_eTCut;
0115 }
0116 
0117 /**
0118  * Returns size of isolation cone as integer multiple of 0.1 (i.e. 3 is an R=0.3 cone)
0119  */
0120 /*const*/ int ClusterIso::getConeSize()
0121 {
0122   return (int) m_coneSize * 10;
0123 }
0124 
0125 /**
0126  * Must be called to set the new vertex for the cluster
0127  */
0128 /*const*/ CLHEP::Hep3Vector ClusterIso::getVertex()
0129 {
0130   return CLHEP::Hep3Vector(m_vx, m_vy, m_vz);
0131 }
0132 
0133 /** \Brief Calculates isolation energy for all electromagnetic calorimeter clusters over the specified eT cut.
0134  *
0135  * For each cluster in the EMCal this iterates through all of the towers in each calorimeter,
0136  * if the towers are within the isolation cone their energy is added to the sum of isolation energy.
0137  * Finally subtract the cluster energy from the sum
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    * If there event is embedded in Au+Au or another larger background we want to
0147    * get isolation energy from the towers with a subtracted background. This first section
0148    * looks at those towers instead of the original objects which include the background.
0149    * NOTE: that during the background event subtraction the EMCal towers are grouped
0150    * together so we have to use the inner HCal geometry.
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       // get EMCal towers
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       // get InnerHCal towers
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       // get outerHCal towers
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       // get geometry of calorimeter towers
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         // vertexmap is used to get correct collision vertex
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           }  // skip if cluster is under eT cut
0240 
0241           // calculate EMCal tower contribution to isolation energy
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);  // if tower is in cone, add energy
0263               }
0264             }
0265           }
0266 
0267           // calculate Inner HCal tower contribution to isolation energy
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);  // if tower is in cone, add energy
0287               }
0288             }
0289           }
0290 
0291           // calculate Outer HCal tower contribution to isolation energy
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);  // if tower is in cone, add energy
0311               }
0312             }
0313           }
0314 
0315           isoEt -= et;  // Subtract cluster eT from isoET
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      * This second section repeats the isolation calculation without any background subtraction
0331      */
0332     if (Verbosity() >= VERBOSITY_EVEN_MORE)
0333     {
0334       std::cout << Name() << "::ClusterIso starting unsubtracted calculation" << '\n';
0335     }
0336     {
0337       // get EMCal towers
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       // get InnerHCal towers
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       // get outerHCal towers
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       // get geometry of calorimeter towers
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         // vertexmap is used to get correct collision vertex
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           }  // skip if cluster is below eT cut
0412 
0413           // calculate EMCal tower contribution to isolation energy
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);  // if tower is in cone, add energy
0433               }
0434             }
0435           }
0436           if (Verbosity() >= VERBOSITY_MAX)
0437           {
0438             std::cout << "\t after EMCal isoEt:" << isoEt << '\n';
0439           }
0440           // calculate Inner HCal tower contribution to isolation energy
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);  // if tower is in cone, add energy
0460               }
0461             }
0462           }
0463           if (Verbosity() >= VERBOSITY_MAX)
0464           {
0465             std::cout << "\t after innerHCal isoEt:" << isoEt << '\n';
0466           }
0467           // calculate Outer HCal tower contribution to isolation energy
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);  // if tower is in cone, add energy
0487               }
0488             }
0489           }
0490           if (Verbosity() >= VERBOSITY_MAX)
0491           {
0492             std::cout << "\t after outerHCal isoEt:" << isoEt << '\n';
0493           }
0494           isoEt -= et;  // Subtract cluster eT from isoET
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 * /*topNode*/)
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 }