packages/ClusterIso/ClusterIso.cc

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 <cmath>
0031 #include <iostream>
0032 #include <map>
0033 #include <utility>
0034
0035 namespace
0036 {
0037   /** \Brief Function to find delta R between 2 objects
0038    *
0039    * Takes the eta and phi of each object and returns the difference
0040    * of the etas and phis added in quadrature. Used to find towers
0041    * inside a cone of delta R around a cluster.
0042    */
0043   float deltaR(float eta1, float eta2, float phi1, float phi2)
0044   {
0045     float deta = eta1 - eta2;
0046     float dphi = phi1 - phi2;
0047     if (dphi > M_PI)
0048     {
0049       dphi -= 2 * M_PI;  // corrects to keep range -pi to pi
0050     }
0051     if (dphi < -1 * M_PI)
0052     {
0053       dphi += 2 * M_PI;  // corrects to keep range -pi to pi
0054     }
0055     return std::sqrt(deta * deta + dphi * dphi);
0056   }
0057 }  // namespace
0058
0059 /** \Brief Function to get correct tower eta
0060  *
0061  * Each calorimeter tower's eta is calculated using the vertex (0,0,0)
0062  * which is incorrect in many collisions. This function
0063  * uses geometry to find a given tower's eta using the correct vertex.
0064  */
0065 double ClusterIso::getTowerEta(RawTowerGeom *tower_geom, double vx, double vy, double vz)
0066 {
0067   float r;
0068   if (vx == 0 && vy == 0 && vz == 0)
0069   {
0070     r = tower_geom->get_eta();
0071   }
0072   else
0073   {
0074     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));
0075     double theta = atan2(radius, tower_geom->get_center_z() - vz);
0076     r = -log(tan(theta / 2.));
0077   }
0078   return r;
0079 }
0080
0081 /**
0082  * Contructor takes the argument of the class name, the minimum eT of the clusters which defaults to 0,
0083  * and the isolation cone size which defaults to 0.3.
0084  */
0085 ClusterIso::ClusterIso(const std::string &kname, float eTCut = 0.0, int coneSize = 3, bool do_subtracted = true, bool do_unsubtracted = true)
0086   : SubsysReco(kname)
0087   , m_do_subtracted(do_subtracted)
0088   , m_do_unsubtracted(do_unsubtracted)
0089   , m_use_towerinfo(true)
0090 {
0091   if (Verbosity() >= VERBOSITY_SOME)
0092   {
0093     std::cout << Name() << "::ClusterIso constructed" << '\n';
0094   }
0095   if (coneSize == 0 && Verbosity() >= VERBOSITY_QUIET)
0096   {
0097     std::cout << "WARNING in " << Name() << "ClusterIso:: cone size is zero" << '\n';
0098   }
0099   m_vx = m_vy = m_vz = 0;
0100   setConeSize(coneSize);
0101   seteTCut(eTCut);
0102   if (Verbosity() >= VERBOSITY_EVEN_MORE)
0103   {
0104     std::cout << Name() << "::ClusterIso::m_coneSize is:" << m_coneSize << '\n';
0105     std::cout << Name() << "::ClusterIso::m_eTCut is:" << m_eTCut << '\n';
0106   }
0107   if (!do_subtracted && !do_unsubtracted && Verbosity() >= VERBOSITY_QUIET)
0108   {
0109     std::cout << "WARNING in " << Name() << "ClusterIso:: all processes turned off doing nothing" << '\n';
0110   }
0111 }
0112
0113 int ClusterIso::Init(PHCompositeNode * /*topNode*/)
0114 {
0115   return 0;
0116 }
0117
0118 /**
0119  * Set the minimum transverse energy required for a cluster to have its isolation calculated
0120  */
0121 void ClusterIso::seteTCut(float eTCut)
0122 {
0123   this->m_eTCut = eTCut;
0124 }
0125
0126 /**
0127  * Set the size of isolation cone as integer multiple of 0.1, (i.e. 3 will use an R=0.3 cone)
0128  */
0129 void ClusterIso::setConeSize(int coneSize)
0130 {
0131   this->m_coneSize = coneSize / 10.0;
0132 }
0133
0134 /**
0135  * Returns the minimum transverse energy required for a cluster to have its isolation calculated
0136  */
0137 float ClusterIso::geteTCut() const
0138 {
0139   return m_eTCut;
0140 }
0141
0142 /**
0143  * Returns size of isolation cone as integer multiple of 0.1 (i.e. 3 is an R=0.3 cone)
0144  */
0145 int ClusterIso::getConeSize() const
0146 {
0147   return (int) m_coneSize * 10;
0148 }
0149
0150 /**
0151  * Must be called to set the new vertex for the cluster
0152  */
0153 CLHEP::Hep3Vector ClusterIso::getVertex() const
0154 {
0155   return CLHEP::Hep3Vector(m_vx, m_vy, m_vz);
0156 }
0157
0158 /** \Brief Calculates isolation energy for all electromagnetic calorimeter clusters over the specified eT cut.
0159  *
0160  * For each cluster in the EMCal this iterates through all of the towers in each calorimeter,
0161  * if the towers are within the isolation cone their energy is added to the sum of isolation energy.
0162  * Finally subtract the cluster energy from the sum
0163  */
0164 int ClusterIso::process_event(PHCompositeNode *topNode)
0165 {
0166   if (Verbosity() >= VERBOSITY_MORE)
0167   {
0168     std::cout << Name() << "::ClusterIso::process_event" << '\n';
0169   }
0170   /**
0171    * If there event is embedded in Au+Au or another larger background we want to
0172    * get isolation energy from the towers with a subtracted background. This first section
0173    * looks at those towers instead of the original objects which include the background.
0174    * NOTE: that during the background event subtraction the EMCal towers are grouped
0175    * together so we have to use the inner HCal geometry.
0176    */
0177
0178   std::string RawCemcClusterNodeName = "CLUSTER_CEMC";
0179   if (m_use_towerinfo)
0180   {
0181     RawCemcClusterNodeName = m_cluster_node_name;
0182   }
0183
0184   if (m_do_subtracted)
0185   {
0186     {
0187       if (Verbosity() >= VERBOSITY_EVEN_MORE)
0188       {
0189         std::cout << Name() << "::ClusterIso starting subtracted calculation" << '\n';
0190       }
0191       // get EMCal towers
0192       TowerInfoContainer *towersEM3old = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC_RETOWER_SUB1");
0193       if (towersEM3old == nullptr)
0194       {
0195         m_do_subtracted = false;
0196         if (Verbosity() >= VERBOSITY_SOME)
0197         {
0198           std::cout << "In " << Name() << "::ClusterIso WARNING substracted towers do not exist subtracted isolation cannot be preformed \n";
0199         }
0200       }
0201       if (Verbosity() >= VERBOSITY_MORE)
0202       {
0203         std::cout << Name() << "::ClusterIso::process_event: " << towersEM3old->size() << " TOWERINFO_CALIB_CEMC_RETOWER_SUB1 towers" << '\n';
0204       }
0205
0206       // get InnerHCal towers
0207       TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN_SUB1");
0208       if (Verbosity() >= VERBOSITY_MORE)
0209       {
0210         std::cout << Name() << "::ClusterIso::process_event: " << towersIH3->size() << " TOWERINFO_CALIB_HCALIN_SUB1 towers" << '\n';
0211       }
0212
0213       // get outerHCal towers
0214       TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT_SUB1");
0215       if (Verbosity() >= VERBOSITY_MORE)
0216       {
0217         std::cout << Name() << "::ClusterIso::process_event: " << towersOH3->size() << " TOWERINFO_CALIB_HCALOUT_SUB1 towers" << std::endl;
0218       }
0219
0220       // get geometry of calorimeter towers
0221       RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0222       RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0223       RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0224
0225       {
0226         RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, RawCemcClusterNodeName);
0227         RawClusterContainer::ConstRange begin_end = clusters->getClusters();
0228         RawClusterContainer::ConstIterator rtiter;
0229         if (Verbosity() >= VERBOSITY_SOME)
0230         {
0231           std::cout << Name() << "::ClusterIso sees " << clusters->size() << " clusters " << '\n';
0232         }
0233
0234         // vertexmap is used to get correct collision vertex
0235         GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0236         m_vx = m_vy = m_vz = 0;
0237         if (vertexmap && !vertexmap->empty())
0238         {
0239           GlobalVertex *vertex = (vertexmap->begin()->second);
0240           m_vx = vertex->get_x();
0241           m_vy = vertex->get_y();
0242           m_vz = vertex->get_z();
0243           if (Verbosity() >= VERBOSITY_SOME)
0244           {
0245             std::cout << Name() << "::ClusterIso Event Vertex Calculated at x:" << m_vx << " y:" << m_vy << " z:" << m_vz << '\n';
0246           }
0247         }
0248
0249         for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0250         {
0251           RawCluster *cluster = rtiter->second;
0252
0253           CLHEP::Hep3Vector vertex(m_vx, m_vy, m_vz);
0254           CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
0255           double cluster_energy = E_vec_cluster.mag();
0256           double cluster_eta = E_vec_cluster.pseudoRapidity();
0257           double cluster_phi = E_vec_cluster.phi();
0258           double et = cluster_energy / cosh(cluster_eta);
0259           double isoEt = 0;
0260
0261           if (et < m_eTCut)
0262           {
0263             continue;
0264           }  // skip if cluster is under eT cut
0265
0266           // calculate EMCal tower contribution to isolation energy
0267           {
0268             unsigned int ntowers = towersEM3old->size();
0269             for (unsigned int channel = 0; channel < ntowers; channel++)
0270             {
0271               TowerInfo *tower = towersEM3old->get_tower_at_channel(channel);
0272               if (!IsAcceptableTower(tower))
0273               {
0274                 continue;
0275               }
0276               if (tower->get_energy() < m_minTowerEnergy)
0277               {
0278                 continue;
0279               }
0280               unsigned int towerkey = towersEM3old->encode_key(channel);
0281               int ieta = towersEM3old->getTowerEtaBin(towerkey);
0282               int iphi = towersEM3old->getTowerPhiBin(towerkey);
0283
0284               const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0285
0286               RawTowerGeom *tower_geom = geomEM->get_tower_geometry(key);
0287               double this_phi = tower_geom->get_phi();
0288               double this_eta = tower_geom->get_eta();
0289               if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0290               {
0291                 isoEt += tower->get_energy() / cosh(this_eta);  // if tower is in cone, add energy
0292               }
0293             }
0294           }
0295
0296           // calculate Inner HCal tower contribution to isolation energy
0297           {
0298             unsigned int ntowers = towersIH3->size();
0299             for (unsigned int channel = 0; channel < ntowers; channel++)
0300             {
0301               TowerInfo *tower = towersIH3->get_tower_at_channel(channel);
0302               if (!IsAcceptableTower(tower))
0303               {
0304                 continue;
0305               }
0306               if (tower->get_energy() < m_minTowerEnergy)
0307               {
0308                 continue;
0309               }
0310               unsigned int towerkey = towersIH3->encode_key(channel);
0311               int ieta = towersIH3->getTowerEtaBin(towerkey);
0312               int iphi = towersIH3->getTowerPhiBin(towerkey);
0313               const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0314               RawTowerGeom *tower_geom = geomIH->get_tower_geometry(key);
0315               double this_phi = tower_geom->get_phi();
0316               double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
0317               if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0318               {
0319                 isoEt += tower->get_energy() / cosh(this_eta);  // if tower is in cone, add energy
0320               }
0321             }
0322           }
0323
0324           // calculate Outer HCal tower contribution to isolation energy
0325           {
0326             unsigned int ntowers = towersOH3->size();
0327             for (unsigned int channel = 0; channel < ntowers; channel++)
0328             {
0329               TowerInfo *tower = towersOH3->get_tower_at_channel(channel);
0330               if (!IsAcceptableTower(tower))
0331               {
0332                 continue;
0333               }
0334               if (tower->get_energy() < m_minTowerEnergy)
0335               {
0336                 continue;
0337               }
0338               unsigned int towerkey = towersOH3->encode_key(channel);
0339               int ieta = towersOH3->getTowerEtaBin(towerkey);
0340               int iphi = towersOH3->getTowerPhiBin(towerkey);
0341               const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0342               RawTowerGeom *tower_geom = geomOH->get_tower_geometry(key);
0343               double this_phi = tower_geom->get_phi();
0344               double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
0345               if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0346               {
0347                 isoEt += tower->get_energy() / cosh(this_eta);  // if tower is in cone, add energy
0348               }
0349             }
0350           }
0351
0352           isoEt -= et;  // Subtract cluster eT from isoET
0353           if (Verbosity() >= VERBOSITY_EVEN_MORE)
0354           {
0355             std::cout << Name() << "::ClusterIso iso_et for ";
0356             cluster->identify();
0357             std::cout << "=" << isoEt << '\n';
0358           }
0359           cluster->set_et_iso(isoEt, 10 * m_coneSize, true, true);
0360         }
0361       }
0362     }
0363   }
0364   if (m_do_unsubtracted)
0365   {
0366     /**
0367      * This second section repeats the isolation calculation without any background subtraction
0368      */
0369     if (Verbosity() >= VERBOSITY_EVEN_MORE)
0370     {
0371       std::cout << Name() << "::ClusterIso starting unsubtracted calculation" << '\n';
0372     }
0373     {
0374       // get EMCal towers
0375       TowerInfoContainer *towersEM3old = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_CEMC");
0376       if (Verbosity() >= VERBOSITY_MORE)
0377       {
0378         std::cout << "ClusterIso::process_event: " << towersEM3old->size() << " TOWERINFO_CALIB_CEMC towers" << '\n';
0379       }
0380
0381       // get InnerHCal towers
0382       TowerInfoContainer *towersIH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALIN");
0383       if (Verbosity() >= VERBOSITY_MORE)
0384       {
0385         std::cout << "ClusterIso::process_event: " << towersIH3->size() << " TOWERINFO_CALIB_HCALIN towers" << '\n';
0386       }
0387
0388       // get outerHCal towers
0389       TowerInfoContainer *towersOH3 = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_HCALOUT");
0390       if (Verbosity() >= VERBOSITY_MORE)
0391       {
0392         std::cout << "ClusterIso::process_event: " << towersOH3->size() << " TOWERINFO_CALIB_HCALOUT towers" << std::endl;
0393       }
0394
0395       // get geometry of calorimeter towers
0396       RawTowerGeomContainer *geomEM = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_CEMC");
0397       RawTowerGeomContainer *geomIH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALIN");
0398       RawTowerGeomContainer *geomOH = findNode::getClass<RawTowerGeomContainer>(topNode, "TOWERGEOM_HCALOUT");
0399
0400       {
0401         RawClusterContainer *clusters = findNode::getClass<RawClusterContainer>(topNode, RawCemcClusterNodeName);
0402         RawClusterContainer::ConstRange begin_end = clusters->getClusters();
0403         RawClusterContainer::ConstIterator rtiter;
0404         if (Verbosity() >= VERBOSITY_SOME)
0405         {
0406           std::cout << " ClusterIso sees " << clusters->size() << " clusters " << '\n';
0407         }
0408
0409         // vertexmap is used to get correct collision vertex
0410         GlobalVertexMap *vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
0411         m_vx = m_vy = m_vz = 0;
0412         if (vertexmap && !vertexmap->empty())
0413         {
0414           GlobalVertex *vertex = (vertexmap->begin()->second);
0415           m_vx = vertex->get_x();
0416           m_vy = vertex->get_y();
0417           m_vz = vertex->get_z();
0418           if (Verbosity() >= VERBOSITY_SOME)
0419           {
0420             std::cout << Name() << "ClusterIso Event Vertex Calculated at x:" << m_vx << " y:" << m_vy << " z:" << m_vz << '\n';
0421           }
0422         }
0423
0424         for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
0425         {
0426           RawCluster *cluster = rtiter->second;
0427
0428           CLHEP::Hep3Vector vertex(m_vx, m_vy, m_vz);
0429           CLHEP::Hep3Vector E_vec_cluster = RawClusterUtility::GetEVec(*cluster, vertex);
0430           double cluster_energy = E_vec_cluster.mag();
0431           double cluster_eta = E_vec_cluster.pseudoRapidity();
0432           double cluster_phi = E_vec_cluster.phi();
0433           double et = cluster_energy / cosh(cluster_eta);
0434           double isoEt = 0;
0435           if (Verbosity() >= VERBOSITY_MAX)
0436           {
0437             std::cout << Name() << "::ClusterIso processing";
0438             cluster->identify();
0439             std::cout << '\n';
0440           }
0441           if (et < m_eTCut)
0442           {
0443             if (Verbosity() >= VERBOSITY_MAX)
0444             {
0445               std::cout << "\t does not pass eT cut" << '\n';
0446             }
0447             continue;
0448           }  // skip if cluster is below eT cut
0449
0450           // calculate EMCal tower contribution to isolation energy
0451           {
0452             unsigned int ntowers = towersEM3old->size();
0453             for (unsigned int channel = 0; channel < ntowers; channel++)
0454             {
0455               TowerInfo *tower = towersEM3old->get_tower_at_channel(channel);
0456               if (!IsAcceptableTower(tower))
0457               {
0458                 continue;
0459               }
0460               if (tower->get_energy() < m_minTowerEnergy)
0461               {
0462                 continue;
0463               }
0464               unsigned int towerkey = towersEM3old->encode_key(channel);
0465               int ieta = towersEM3old->getTowerEtaBin(towerkey);
0466               int iphi = towersEM3old->getTowerPhiBin(towerkey);
0467               const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::CEMC, ieta, iphi);
0468               RawTowerGeom *tower_geom = geomEM->get_tower_geometry(key);
0469               double this_phi = tower_geom->get_phi();
0470               double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
0471               if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0472               {
0473                 isoEt += tower->get_energy() / cosh(this_eta);  // if tower is in cone, add energy
0474               }
0475             }
0476           }
0477           if (Verbosity() >= VERBOSITY_MAX)
0478           {
0479             std::cout << "\t after EMCal isoEt:" << isoEt << '\n';
0480           }
0481           // calculate Inner HCal tower contribution to isolation energy
0482           {
0483             unsigned int ntowers = towersIH3->size();
0484             for (unsigned int channel = 0; channel < ntowers; channel++)
0485             {
0486               TowerInfo *tower = towersIH3->get_tower_at_channel(channel);
0487               if (!IsAcceptableTower(tower))
0488               {
0489                 continue;
0490               }
0491               if (tower->get_energy() < m_minTowerEnergy)
0492               {
0493                 continue;
0494               }
0495               unsigned int towerkey = towersIH3->encode_key(channel);
0496               int ieta = towersIH3->getTowerEtaBin(towerkey);
0497               int iphi = towersIH3->getTowerPhiBin(towerkey);
0498               const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALIN, ieta, iphi);
0499               RawTowerGeom *tower_geom = geomIH->get_tower_geometry(key);
0500               double this_phi = tower_geom->get_phi();
0501               double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
0502               if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0503               {
0504                 isoEt += tower->get_energy() / cosh(this_eta);  // if tower is in cone, add energy
0505               }
0506             }
0507           }
0508           if (Verbosity() >= VERBOSITY_MAX)
0509           {
0510             std::cout << "\t after innerHCal isoEt:" << isoEt << '\n';
0511           }
0512           // calculate Outer HCal tower contribution to isolation energy
0513           {
0514             unsigned int ntowers = towersOH3->size();
0515             for (unsigned int channel = 0; channel < ntowers; channel++)
0516             {
0517               TowerInfo *tower = towersOH3->get_tower_at_channel(channel);
0518               if (!IsAcceptableTower(tower))
0519               {
0520                 continue;
0521               }
0522               if (tower->get_energy() < m_minTowerEnergy)
0523               {
0524                 continue;
0525               }
0526               unsigned int towerkey = towersOH3->encode_key(channel);
0527               int ieta = towersOH3->getTowerEtaBin(towerkey);
0528               int iphi = towersOH3->getTowerPhiBin(towerkey);
0529               const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(RawTowerDefs::CalorimeterId::HCALOUT, ieta, iphi);
0530               RawTowerGeom *tower_geom = geomOH->get_tower_geometry(key);
0531               double this_phi = tower_geom->get_phi();
0532               double this_eta = getTowerEta(tower_geom, m_vx, m_vy, m_vz);
0533               if (deltaR(cluster_eta, this_eta, cluster_phi, this_phi) < m_coneSize)
0534               {
0535                 isoEt += tower->get_energy() / cosh(this_eta);  // if tower is in cone, add energy
0536               }
0537             }
0538           }
0539           if (Verbosity() >= VERBOSITY_MAX)
0540           {
0541             std::cout << "\t after outerHCal isoEt:" << isoEt << '\n';
0542           }
0543           isoEt -= et;  // Subtract cluster eT from isoET
0544           if (Verbosity() >= VERBOSITY_EVEN_MORE)
0545           {
0546             std::cout << Name() << "::ClusterIso iso_et for ";
0547             cluster->identify();
0548             std::cout << "=" << isoEt << '\n';
0549           }
0550           cluster->set_et_iso(isoEt, 10 * m_coneSize, false, true);
0551         }
0552       }
0553     }
0554   }
0555   return 0;
0556 }
0557
0558 int ClusterIso::End(PHCompositeNode * /*topNode*/)
0559 {
0560   return 0;
0561 }
0562
0563 bool ClusterIso::IsAcceptableTower(TowerInfo *tower)
0564 {
0565   if (!tower->get_isGood())
0566   {
0567     return false;
0568   }
0569   return true;
0570 }