packages/CaloReco/RawClusterDeadHotMask.cc

0001 #include "RawClusterDeadHotMask.h"
0002
0003 #include <calobase/RawCluster.h>
0004 #include <calobase/RawClusterContainer.h>
0005 #include <calobase/RawTower.h>
0006 #include <calobase/RawTowerContainer.h>
0007 #include <calobase/TowerInfo.h>
0008 #include <calobase/TowerInfoContainer.h>
0009 #include <calobase/TowerInfoDefs.h>
0010 #include <calobase/RawTowerDeadMap.h>
0011 #include <calobase/RawTowerGeomContainer.h>
0012 #include <calobase/RawTowerGeom.h>
0013
0014 #include <fun4all/Fun4AllBase.h>
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016 #include <fun4all/SubsysReco.h>
0017
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHNode.h>
0020 #include <phool/PHNodeIterator.h>
0021 #include <phool/getClass.h>
0022
0023 #include <cassert>
0024 #include <cmath>
0025 #include <fstream>
0026 #include <iostream>
0027 #include <map>
0028 #include <sstream>
0029 #include <stdexcept>
0030 #include <string>
0031 #include <utility>  // for pair
0032 #include <vector>
0033
0034 RawClusterDeadHotMask::RawClusterDeadHotMask(const std::string &name)
0035   : SubsysReco(name)
0036   , m_detector("NONE")
0037   , m_deadTowerMaskHalfWidth(1.6)
0038   , m_rawClusters(nullptr)
0039   , m_towerinfoClusters(nullptr)
0040   , m_deadMap(nullptr)
0041   , m_calibTowers(nullptr)
0042   , m_calibTowerInfos(nullptr)
0043   , m_geometry(nullptr)
0044 {
0045 }
0046
0047 int RawClusterDeadHotMask::InitRun(PHCompositeNode *topNode)
0048 {
0049   CreateNodeTree(topNode);
0050
0051   return Fun4AllReturnCodes::EVENT_OK;
0052 }
0053
0054 int RawClusterDeadHotMask::process_event(PHCompositeNode * /*topNode*/)
0055 {
0056   if (Verbosity() >= VERBOSITY_SOME)
0057   {
0058     std::cout << Name() << "::" << m_detector << "::process_event - Entry" << std::endl;
0059   }
0060   int nMasked = 0;
0061
0062   const int eta_bins = m_geometry->get_etabins();
0063   const int phi_bins = m_geometry->get_phibins();
0064   assert(eta_bins > 0);
0065   assert(phi_bins > 0);
0066
0067   // loop over the clusters
0068   RawClusterContainer::ConstRange begin_end;
0069   if(m_UseTowerInfo)
0070   {
0071     begin_end = m_towerinfoClusters->getClusters();
0072   }
0073   else
0074   {
0075     begin_end = m_rawClusters->getClusters();
0076   }
0077   RawClusterContainer::ConstIterator iter;
0078
0079   for (iter = begin_end.first; iter != begin_end.second;)
0080   {
0081     //    RawClusterDefs::keytype key = iter->first;
0082     const RawCluster *cluster = iter->second;
0083
0084     std::vector<float> toweretas;
0085     std::vector<float> towerphis;
0086     std::vector<float> towerenergies;
0087
0088     // loop over the towers in the cluster
0089     RawCluster::TowerConstRange towers = cluster->get_towers();
0090     RawCluster::TowerConstIterator toweriter;
0091
0092     for (toweriter = towers.first; toweriter != towers.second; ++toweriter)
0093     {
0094       if(m_UseTowerInfo)
0095       {
0096         int towereta = m_geometry->get_tower_geometry(toweriter->first)->get_bineta();
0097         int towerphi = m_geometry->get_tower_geometry(toweriter->first)->get_binphi();
0098         unsigned int key = TowerInfoDefs::encode_emcal(towereta, towerphi);
0099         TowerInfo *towerinfo = m_calibTowerInfos->get_tower_at_key(key);
0100
0101         assert(towerinfo);
0102
0103         double towerenergy = towerinfo->get_energy();
0104
0105         toweretas.push_back(towereta);
0106         towerphis.push_back(towerphi);
0107         towerenergies.push_back(towerenergy);
0108       }
0109       else
0110       {
0111         RawTower *tower = m_calibTowers->getTower(toweriter->first);
0112
0113         assert(tower);
0114
0115         int towereta = tower->get_bineta();
0116         int towerphi = tower->get_binphi();
0117         double towerenergy = tower->get_energy();
0118
0119         // put the etabin, phibin, and energy into the corresponding vectors
0120         toweretas.push_back(towereta);
0121         towerphis.push_back(towerphi);
0122         towerenergies.push_back(towerenergy);
0123       }
0124     }
0125
0126     int ntowers = toweretas.size();
0127     assert(ntowers >= 1);
0128
0129     // loop over the towers to determine the energy
0130     // weighted eta and phi position of the cluster
0131
0132     float etamult = 0;
0133     float etasum = 0;
0134     float phimult = 0;
0135     float phisum = 0;
0136
0137     for (int j = 0; j < ntowers; j++)
0138     {
0139       float energymult = towerenergies.at(j) * toweretas.at(j);
0140       etamult += energymult;
0141       etasum += towerenergies.at(j);
0142
0143       int phibin = towerphis.at(j);
0144
0145       if (phibin - towerphis.at(0) < -phi_bins / 2.0)
0146       {
0147         phibin += phi_bins;
0148       }
0149       else if (phibin - towerphis.at(0) > +phi_bins / 2.0)
0150       {
0151         phibin -= phi_bins;
0152       }
0153       assert(std::abs(phibin - towerphis.at(0)) <= phi_bins / 2.0);
0154
0155       energymult = towerenergies.at(j) * phibin;
0156       phimult += energymult;
0157       phisum += towerenergies.at(j);
0158     }
0159
0160     float avgphi = phimult / phisum;
0161     float avgeta = etamult / etasum;
0162
0163     if (Verbosity() > VERBOSITY_MORE)
0164     {
0165       std::cout << Name() << "::" << m_detector << "::process_event - process cluster at average location " << avgeta << "," << avgphi << " : ";
0166       cluster->identify();
0167     }
0168
0169     // rejection if close to a dead tower
0170     bool rejecCluster = false;
0171
0172     for (int search_eta = ceil(avgeta - m_deadTowerMaskHalfWidth); search_eta <= floor(avgeta + m_deadTowerMaskHalfWidth); ++search_eta)
0173     {
0174       for (int search_phi = ceil(avgphi - m_deadTowerMaskHalfWidth); search_phi <= floor(avgphi + m_deadTowerMaskHalfWidth); ++search_phi)
0175       {
0176         int ieta = search_eta;
0177         int iphi = search_phi;
0178
0179         if (ieta >= eta_bins)
0180         {
0181           continue;
0182         }
0183         if (ieta < 0)
0184         {
0185           continue;
0186         }
0187
0188         if (iphi >= phi_bins)
0189         {
0190           iphi -= phi_bins;
0191         }
0192         if (iphi < 0)
0193         {
0194           iphi += phi_bins;
0195         }
0196
0197         const bool isDead = m_deadMap->isDeadTower(ieta, iphi);
0198
0199         // dead tower found in cluster
0200         if (Verbosity() > VERBOSITY_MORE)
0201         {
0202           std::cout << "\t"
0203                     << "tower " << ieta
0204                     << "," << iphi
0205                     << (isDead ? ": is dead." : "OK")
0206                     << std::endl;
0207         }
0208
0209         if (isDead)
0210         {
0211           rejecCluster = true;
0212           break;
0213         }
0214       }
0215
0216       if (rejecCluster)
0217       {
0218         break;
0219       }
0220     }
0221
0222     // container operation
0223     if (rejecCluster)
0224     {
0225       if (Verbosity() > VERBOSITY_MORE)
0226       {
0227         std::cout << Name() << "::" << m_detector << "::process_event - " << "reject cluster " << cluster->get_id() << std::endl;
0228         cluster->identify();
0229       }
0230
0231       ++nMasked;
0232       if(m_UseTowerInfo)
0233       {
0234         m_towerinfoClusters->getClustersMap().erase(iter++);
0235       }
0236       else
0237       {
0238         m_rawClusters->getClustersMap().erase(iter++);
0239       }
0240     }
0241     else
0242     {
0243       if (Verbosity() > VERBOSITY_MORE)
0244       {
0245         std::cout << Name() << "::" << m_detector << "::process_event - " << "keep cluster " << cluster->get_id() << std::endl;
0246       }
0247       ++iter;
0248     }
0249
0250   }  //  for (iter = begin_end.first; iter != begin_end.second;)
0251
0252   if (Verbosity() >= VERBOSITY_SOME)
0253   {
0254     unsigned int clus_size;
0255     if(m_UseTowerInfo)
0256     {
0257       clus_size = m_towerinfoClusters->size();
0258     }
0259     else
0260     {
0261       clus_size = m_rawClusters->size();
0262     }
0263     std::cout << Name() << "::" << m_detector << "::process_event - masked "
0264               << nMasked << " clusters. Final cluster containers has "
0265               << clus_size << " clusters"
0266               << std::endl;
0267   }
0268   return Fun4AllReturnCodes::EVENT_OK;
0269 }
0270
0271 void RawClusterDeadHotMask::CreateNodeTree(PHCompositeNode *topNode)
0272 {
0273   PHNodeIterator iter(topNode);
0274
0275   // Get the DST Node
0276   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0277
0278   // Check that it is there
0279   if (!dstNode)
0280   {
0281     std::cout << "DST Node missing, quitting" << std::endl;
0282     throw std::runtime_error("failed to find DST node in RawClusterDeadHotMask::CreateNodeTree");
0283   }
0284
0285   m_rawClusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_" + m_detector);
0286   if (!m_rawClusters && m_UseTowerInfo == false)
0287   {
0288     std::cout << Name() << "::" << m_detector << "::"
0289               << "CreateNodeTree "
0290                  "No "
0291               << m_detector << " Cluster Container found while in RawClusterDeadHotMask, can't proceed!!!" << std::endl;
0292     topNode->print();
0293     throw std::runtime_error("failed to find CLUSTER node in RawClusterDeadHotMask::CreateNodeTree");
0294   }
0295
0296   m_towerinfoClusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTERINFO_" + m_detector);
0297   if (!m_towerinfoClusters && m_UseTowerInfo == true)
0298   {
0299     std::cout << Name() << "::" << m_detector << "::"
0300               << "CreateNodeTree "
0301                  "No "
0302               << m_detector << " Cluster Container found while in RawClusterDeadHotMask, can't proceed!!!" << std::endl;
0303     topNode->print();
0304     throw std::runtime_error("failed to find CLUSTERINFO node in RawClusterDeadHotMask::CreateNodeTree");
0305   }
0306
0307   m_calibTowers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_" + m_detector);
0308   if (!m_calibTowers  && m_UseTowerInfo == false)
0309   {
0310     std::cout << Name() << "::" << m_detector << "::"
0311               << "CreateNodeTree "
0312               << "No calibrated " << m_detector << " tower info found while in RawClusterDeadHotMask, can't proceed!!!" << std::endl;
0313     throw std::runtime_error("failed to find TOWER_CALIB node in RawClusterDeadHotMask::CreateNodeTree");
0314   }
0315
0316   m_calibTowerInfos = findNode::getClass<TowerInfoContainer>(topNode, "TOWERINFO_CALIB_" + m_detector);
0317   if (!m_calibTowerInfos && m_UseTowerInfo == true)
0318   {
0319     std::cout << Name() << "::" << m_detector << "::"
0320               << "CreateNodeTree "
0321               << "No calibrated " << m_detector << " tower info found while in RawClusterDeadHotMask, can't proceed!!!" << std::endl;
0322     throw std::runtime_error("failed to find TOWERINFO_CALIB node in RawClusterDeadHotMask::CreateNodeTree");
0323   }
0324
0325   std::string towergeomnodename = "TOWERGEOM_" + m_detector;
0326   m_geometry = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
0327   if (!m_geometry)
0328   {
0329     std::cout << Name() << "::" << m_detector << "::"
0330               << "CreateNodeTree"
0331               << ": Could not find node " << towergeomnodename << std::endl;
0332     throw std::runtime_error("failed to find TOWERGEOM node in RawClusterDeadHotMask::CreateNodeTree");
0333   }
0334
0335   const std::string deadMapName = "DEADMAP_" + m_detector;
0336   m_deadMap = findNode::getClass<RawTowerDeadMap>(topNode, deadMapName);
0337   if (m_deadMap)
0338   {
0339     std::cout << Name() << "::" << m_detector << "::"
0340               << "CreateNodeTree"
0341               << " use dead map: ";
0342     m_deadMap->identify();
0343   }
0344 }
0345
0346 int RawClusterDeadHotMask::End(PHCompositeNode * /*topNode*/)
0347 {
0348   return Fun4AllReturnCodes::EVENT_OK;
0349 }