packages/CaloReco/RawTowerDeadTowerInterp.cc

0001 #include "RawTowerDeadTowerInterp.h"
0002
0003 #include <calobase/TowerInfo.h>
0004 #include <calobase/TowerInfoContainer.h>
0005 #include <calobase/RawTowerDeadMap.h>
0006 #include <calobase/RawTowerDefs.h>
0007 #include <calobase/RawTowerGeom.h>
0008 #include <calobase/RawTowerGeomContainer.h>
0009
0010 #include <fun4all/Fun4AllBase.h>
0011 #include <fun4all/Fun4AllReturnCodes.h>
0012 #include <fun4all/SubsysReco.h>
0013
0014 #include <phool/PHCompositeNode.h>
0015 #include <phool/PHNode.h>
0016 #include <phool/PHNodeIterator.h>
0017 #include <phool/getClass.h>
0018
0019 #include <cassert>
0020 #include <cstdlib>
0021 #include <exception>
0022 #include <iostream>
0023 #include <set>  // for _Rb_tree_const_iterator
0024 #include <stdexcept>
0025 #include <string>
0026 #include <utility>
0027 #include <vector>
0028
0029 RawTowerDeadTowerInterp::RawTowerDeadTowerInterp(const std::string &name)
0030   : SubsysReco(name)
0031 {
0032 }
0033
0034 int RawTowerDeadTowerInterp::InitRun(PHCompositeNode *topNode)
0035 {
0036   PHNodeIterator iter(topNode);
0037
0038   // Looking for the DST node
0039   PHCompositeNode *dstNode;
0040   dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode",
0041                                                            "DST"));
0042   if (!dstNode)
0043   {
0044     std::cout << Name() << "::" << m_detector << "::" << __PRETTY_FUNCTION__
0045               << "DST Node missing, doing nothing." << std::endl;
0046     exit(1);
0047   }
0048
0049   try
0050   {
0051     CreateNodes(topNode);
0052   }
0053   catch (std::exception &e)
0054   {
0055     std::cout << e.what() << std::endl;
0056     return Fun4AllReturnCodes::ABORTRUN;
0057   }
0058   return Fun4AllReturnCodes::EVENT_OK;
0059 }
0060
0061 int RawTowerDeadTowerInterp::process_event(PHCompositeNode * /*topNode*/)
0062 {
0063   if (Verbosity())
0064   {
0065     std::cout << Name() << "::" << m_detector << "::"
0066               << "process_event"
0067               << "Process event entered" << std::endl;
0068   }
0069
0070   double recovery_energy = 0;
0071   int recoverTower = 0;
0072   int deadTowerCnt = 0;
0073
0074   if (m_deadTowerMap)
0075   {
0076     const int eta_bins = m_geometry->get_etabins();
0077     const int phi_bins = m_geometry->get_phibins();
0078
0079     // assume cylindrical calorimeters for now. Need to add swithc for planary calorimeter
0080     assert(eta_bins > 0);
0081     assert(phi_bins > 0);
0082
0083     for (const RawTowerDefs::keytype &key : m_deadTowerMap->getDeadTowers())
0084     {
0085       ++deadTowerCnt;
0086
0087       if (Verbosity() >= VERBOSITY_MORE)
0088       {
0089         std::cout << Name() << "::" << m_detector << "::"
0090                   << "process_event"
0091                   << " - processing tower " << key;
0092       }
0093
0094       // check deadmap validity
0095       RawTowerGeom *towerGeom = m_geometry->get_tower_geometry(key);
0096       if (towerGeom == nullptr)
0097       {
0098         std::cout << Name() << "::" << m_detector << "::"
0099                   << "process_event"
0100                   << " - invalid dead tower ID " << key << std::endl;
0101
0102         exit(2);
0103       }
0104
0105       const int bineta = towerGeom->get_bineta();
0106       const int binphi = towerGeom->get_binphi();
0107
0108       if (Verbosity() >= VERBOSITY_MORE)
0109       {
0110         std::cout << " bin " << bineta << "-" << binphi;
0111         std::cout << ". Add neighbors: ";
0112       }
0113
0114       assert(bineta >= 0);
0115       assert(bineta <= eta_bins);
0116       assert(binphi >= 0);
0117       assert(binphi <= phi_bins);
0118
0119       // eight neighbors
0120       static const std::vector<std::pair<int, int>> neighborIndexs =
0121           {{+1, 0}, {+1, +1}, {0, +1}, {-1, +1}, {-1, 0}, {-1, -1}, {0, -1}, {+1, -1}};
0122
0123       int n_neighbor = 0;
0124       double E_SumNeighbor = 0;
0125       for (const std::pair<int, int> &neighborIndex : neighborIndexs)
0126       {
0127         int ieta = bineta + neighborIndex.first;
0128         int iphi = binphi + neighborIndex.second;
0129
0130         if (ieta >= eta_bins)
0131         {
0132           continue;
0133         }
0134         if (ieta < 0)
0135         {
0136           continue;
0137         }
0138
0139         if (iphi >= phi_bins)
0140         {
0141           iphi -= phi_bins;
0142         }
0143         if (iphi < 0)
0144         {
0145           iphi += phi_bins;
0146         }
0147
0148         assert(ieta >= 0);
0149         assert(ieta <= eta_bins);
0150         assert(iphi >= 0);
0151         assert(iphi <= phi_bins);
0152
0153         if (m_deadTowerMap->isDeadTower(ieta, iphi))
0154         {
0155           continue;
0156         }
0157
0158         ++n_neighbor;
0159
0160         assert(m_calibTowers);
0161     unsigned int towerkey = ((unsigned int) (ieta) << 16U) + iphi;
0162     TowerInfo *neighTower = m_calibTowers->get_tower_at_key(towerkey);
0163         // RawTower *neighTower = m_calibTowers->getTower(ieta, iphi);
0164         if (neighTower == nullptr)
0165     {
0166       continue;
0167     }
0168
0169         if (Verbosity() >= VERBOSITY_MORE)
0170         {
0171           std::cout << neighTower->get_energy() << " (" << ieta << "-" << iphi << "), ";
0172         }
0173
0174         E_SumNeighbor += neighTower->get_energy();
0175
0176       }  // for (const pair<int, int> &neighborIndex : neighborIndexs)
0177
0178       if (n_neighbor > 0 and E_SumNeighbor != 0)
0179       {
0180
0181     unsigned int deadtowerkey = ((unsigned int) (bineta) << 16U) + binphi;
0182
0183         TowerInfo *deadTower = m_calibTowers->get_tower_at_key(deadtowerkey);
0184
0185         // if (deadTower == nullptr)
0186         // {
0187         //   deadTower = new TowerInfo();
0188         // }
0189         assert(deadTower);
0190
0191         deadTower->set_energy(E_SumNeighbor / n_neighbor);
0192         // m_calibTowers->AddTower(key, deadTower);
0193
0194         recovery_energy += deadTower->get_energy();
0195         ++recoverTower;
0196
0197         if (Verbosity() >= VERBOSITY_MORE)
0198         {
0199           std::cout << " -> " << deadTower->get_energy() << " GeV @ etabin: " << bineta << " , phi bin: " << binphi ;
0200         }
0201
0202       }  // if (n_neighbor>0)
0203       else
0204       {
0205         if (Verbosity() >= VERBOSITY_MORE)
0206         {
0207           std::cout << "No neighbor towers found.";
0208         }
0209
0210       }  // if (n_neighbor>0) -else
0211
0212       if (Verbosity() >= VERBOSITY_MORE)
0213       {
0214         std::cout << std::endl;
0215       }
0216     }
0217
0218   }  // if (m_deadTowerMap)
0219   else
0220   {
0221     static bool once = true;
0222
0223     if (Verbosity() and once)
0224     {
0225       once = false;
0226
0227       std::cout << Name() << "::" << m_detector << "::"
0228                 << "process_event"
0229                 << " - missing dead map node. Do nothing ..."
0230                 << std::endl;
0231     }
0232   }
0233
0234   if (Verbosity())
0235   {
0236     std::cout << Name() << "::" << m_detector << "::"
0237               << "process_event"
0238               << "recovery_energy = " << recovery_energy
0239               << " GeV from " << recoverTower << " towers out of total " << deadTowerCnt << " dead towers"
0240                << std::endl;
0241   }
0242   return Fun4AllReturnCodes::EVENT_OK;
0243 }
0244
0245 void RawTowerDeadTowerInterp::CreateNodes(PHCompositeNode *topNode)
0246 {
0247   PHNodeIterator iter(topNode);
0248   PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
0249   if (!runNode)
0250   {
0251     std::cout << Name() << "::" << m_detector << "::"
0252               << "CreateNodes"
0253               << "Run Node missing, doing nothing." << std::endl;
0254     throw std::runtime_error(
0255         "Failed to find Run node in RawTowerDeadTowerInterp::CreateNodes");
0256   }
0257
0258   const std::string deadMapName = "DEADMAP_" + m_detector;
0259   m_deadTowerMap = findNode::getClass<RawTowerDeadMap>(topNode, deadMapName);
0260   if (m_deadTowerMap)
0261   {
0262     std::cout << Name() << "::" << m_detector << "::"
0263               << "CreateNodes"
0264               << " use dead map: ";
0265     m_deadTowerMap->identify();
0266   }
0267
0268   const std::string geometry_node = "TOWERGEOM_" + m_detector;
0269   m_geometry = findNode::getClass<RawTowerGeomContainer>(topNode, geometry_node);
0270   if (!m_geometry)
0271   {
0272     std::cout << Name() << "::" << m_detector << "::"
0273               << "CreateNodes"
0274               << " " << geometry_node << " Node missing, doing bail out!"
0275               << std::endl;
0276     throw std::runtime_error(
0277         "Failed to find " + geometry_node + " node in RawTowerDeadTowerInterp::CreateNodes");
0278   }
0279
0280   if (Verbosity() >= 1)
0281   {
0282     m_geometry->identify();
0283   }
0284
0285   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst(
0286       "PHCompositeNode", "DST"));
0287   if (!dstNode)
0288   {
0289     std::cout << Name() << "::" << m_detector << "::"
0290               << "CreateNodes"
0291               << "DST Node missing, doing nothing." << std::endl;
0292     throw std::runtime_error(
0293         "Failed to find DST node in RawTowerDeadTowerInterp::CreateNodes");
0294   }
0295
0296   const std::string rawTowerNodeName = "TOWERINFO_" + _calib_tower_node_prefix + "_" + m_detector;
0297   m_calibTowers = findNode::getClass<TowerInfoContainer>(dstNode, rawTowerNodeName);
0298   if (!m_calibTowers)
0299   {
0300     std::cout << Name() << "::" << m_detector << "::" << __PRETTY_FUNCTION__
0301               << " " << rawTowerNodeName << " Node missing, doing bail out!"
0302               << std::endl;
0303     throw std::runtime_error(
0304         "Failed to find " + rawTowerNodeName + " node in RawTowerDeadTowerInterp::CreateNodes");
0305   }
0306
0307   return;
0308 }