packages/CaloReco/RawClusterBuilderGraph.cc

0001 #include "RawClusterBuilderGraph.h"
0002
0003 #include "PHMakeGroups.h"
0004
0005 #include <calobase/RawCluster.h>
0006 #include <calobase/RawClusterContainer.h>
0007 #include <calobase/RawClusterDefs.h>
0008 #include <calobase/RawClusterv1.h>
0009 #include <calobase/RawTower.h>
0010 #include <calobase/RawTowerContainer.h>
0011 #include <calobase/RawTowerDefs.h>
0012 #include <calobase/RawTowerGeom.h>
0013 #include <calobase/RawTowerGeomContainer.h>
0014
0015 #include <fun4all/Fun4AllReturnCodes.h>
0016 #include <fun4all/SubsysReco.h>
0017
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/PHIODataNode.h>
0020 #include <phool/PHNode.h>
0021 #include <phool/PHNodeIterator.h>
0022 #include <phool/PHObject.h>
0023 #include <phool/getClass.h>
0024 #include <phool/phool.h>
0025
0026 #include <algorithm>                         // for sort
0027 #include <cassert>
0028 #include <cmath>
0029 #include <exception>
0030 #include <iostream>
0031 #include <map>
0032 #include <stdexcept>
0033 #include <utility>
0034 #include <vector>
0035
0036 // this is just a helper class which enables us to handle rollovers
0037 // when checking for adjacent towers, it requires one bit of
0038 // information (the total number of phibins) which
0039 // is not in the tower class
0040 // NOLINTNEXTLINE(hicpp-special-member-functions)
0041 class twrs
0042 {
0043  public:
0044   explicit twrs(RawTower * /*rt*/);
0045   virtual ~twrs() = default;
0046   bool is_adjacent(const twrs & /*tower*/) const;
0047   void set_id(const int i)
0048   {
0049     id = i;
0050   }
0051   int get_id() const
0052   {
0053     return id;
0054   }
0055   void set_maxphibin(const int i)
0056   {
0057     maxphibin = i;
0058   }
0059   int get_maxphibin() const
0060   {
0061     return maxphibin;
0062   }
0063   int get_bineta() const
0064   {
0065     return bineta;
0066   }
0067   int get_binphi() const
0068   {
0069     return binphi;
0070   }
0071
0072  private:
0073   int bineta;
0074   int binphi;
0075   int maxphibin;
0076   RawTowerDefs::keytype id;
0077 };
0078
0079 twrs::twrs(RawTower *rt)
0080   : bineta(rt->get_bineta()), binphi(rt->get_binphi()), maxphibin(-10)
0081   , id(-1)
0082 {
0083 }
0084
0085 bool twrs::is_adjacent(const twrs &tower) const
0086 {
0087   if (bineta - 1 <= tower.get_bineta() && tower.get_bineta() <= bineta + 1)
0088   {
0089     if (binphi - 1 <= tower.get_binphi() && tower.get_binphi() <= binphi + 1)
0090 // NOLINTNEXTLINE(bugprone-branch-clone)
0091     {
0092       return true;
0093     }
0094     // cluster through the phi-wraparound
0095     if (((tower.get_binphi() == maxphibin - 1) && (binphi == 0)) ||
0096              ((tower.get_binphi() == 0) && (binphi == maxphibin - 1)))
0097     {
0098       return true;
0099     }
0100   }
0101
0102   return false;
0103 }
0104
0105 static bool operator<(const twrs &a, const twrs &b)// NOLINT(misc-use-anonymous-namespace)
0106 {
0107   if (a.get_bineta() != b.get_bineta())
0108   {
0109     return a.get_bineta() < b.get_bineta();
0110   }
0111   return a.get_binphi() < b.get_binphi();
0112 }
0113
0114 RawClusterBuilderGraph::RawClusterBuilderGraph(const std::string &name)
0115   : SubsysReco(name)
0116 {
0117 }
0118
0119 int RawClusterBuilderGraph::InitRun(PHCompositeNode *topNode)
0120 {
0121   try
0122   {
0123     CreateNodes(topNode);
0124   }
0125   catch (std::exception &e)
0126   {
0127     std::cout << PHWHERE << ": " << e.what() << std::endl;
0128     throw;
0129   }
0130
0131   return Fun4AllReturnCodes::EVENT_OK;
0132 }
0133
0134 int RawClusterBuilderGraph::process_event(PHCompositeNode *topNode)
0135 {
0136   std::string towernodename = "TOWER_CALIB_" + detector;
0137   // Grab the towers
0138   RawTowerContainer *towers = findNode::getClass<RawTowerContainer>(topNode, towernodename);
0139   if (!towers)
0140   {
0141     std::cout << PHWHERE << ": Could not find node " << towernodename << std::endl;
0142     return Fun4AllReturnCodes::DISCARDEVENT;
0143   }
0144   std::string towergeomnodename = "TOWERGEOM_" + detector;
0145   RawTowerGeomContainer *towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
0146   if (!towergeom)
0147   {
0148     std::cout << PHWHERE << ": Could not find node " << towergeomnodename << std::endl;
0149     return Fun4AllReturnCodes::ABORTEVENT;
0150   }
0151   // make the list of towers above threshold
0152   std::vector<twrs> towerVector;
0153   RawTowerContainer::ConstRange begin_end = towers->getTowers();
0154   RawTowerContainer::ConstIterator itr = begin_end.first;
0155   for (; itr != begin_end.second; ++itr)
0156   {
0157     RawTower *tower = itr->second;
0158     RawTowerDefs::keytype towerid = itr->first;
0159     if (tower->get_energy() > _min_tower_e)
0160     {
0161       twrs twr(tower);
0162       twr.set_maxphibin(towergeom->get_phibins());
0163       twr.set_id(towerid);
0164       towerVector.push_back(twr);
0165     }
0166   }
0167
0168   // cluster the towers
0169   std::multimap<int, twrs> clusteredTowers;
0170   PHMakeGroups(towerVector, clusteredTowers);
0171
0172   RawCluster *cluster = nullptr;
0173   int last_id = -1;
0174   std::multimap<int, twrs>::iterator ctitr = clusteredTowers.begin();
0175   std::multimap<int, twrs>::iterator lastct = clusteredTowers.end();
0176   for (; ctitr != lastct; ++ctitr)
0177   {
0178     int clusterid = ctitr->first;
0179
0180     if (last_id != clusterid)
0181     {
0182       // new cluster
0183       cluster = new RawClusterv1();
0184       _clusters->AddCluster(cluster);
0185
0186       last_id = clusterid;
0187     }
0188     assert(cluster);
0189
0190     const twrs &tmptower = ctitr->second;
0191     RawTower *rawtower = towers->getTower(tmptower.get_id());
0192
0193     const double e = rawtower->get_energy();
0194     cluster->addTower(rawtower->get_id(), e);
0195   }
0196
0197   for (const auto &cluster_pair : _clusters->getClustersMap())
0198   {
0199     RawClusterDefs::keytype clusterid = cluster_pair.first;
0200     RawCluster *clusterA = cluster_pair.second;
0201
0202     assert(clusterA);
0203     assert(clusterA->get_id() == clusterid);
0204
0205     double sum_x(0);
0206     double sum_y(0);
0207     double sum_z(0);
0208     double sum_e(0);
0209
0210     for (const auto tower_pair : clusterA->get_towermap())
0211     {
0212       const RawTower *rawtower = towers->getTower(tower_pair.first);
0213       const RawTowerGeom *rawtowergeom = towergeom->get_tower_geometry(tower_pair.first);
0214
0215       assert(rawtower);
0216       assert(rawtowergeom);
0217       const double e = rawtower->get_energy();
0218
0219       sum_e += e;
0220
0221       if (e > 0)
0222       {
0223         sum_x += e * rawtowergeom->get_center_x();
0224         sum_y += e * rawtowergeom->get_center_y();
0225         sum_z += e * rawtowergeom->get_center_z();
0226       }
0227     }  //     for (const auto tower_pair : clusterA->get_towermap())
0228
0229     clusterA->set_energy(sum_e);
0230
0231     if (sum_e > 0)
0232     {
0233       sum_x /= sum_e;
0234       sum_y /= sum_e;
0235       sum_z /= sum_e;
0236
0237       clusterA->set_r(sqrt((sum_y * sum_y) + (sum_x * sum_x)));
0238       clusterA->set_phi(atan2(sum_y, sum_x));
0239       clusterA->set_z(sum_z);
0240     }
0241
0242     if (Verbosity() > 1)
0243     {
0244       std::cout << "RawClusterBuilderGraph constucted ";
0245       clusterA->identify();
0246     }
0247   }  //  for (const auto & cluster_pair : _clusters->getClustersMap())
0248
0249   if (chkenergyconservation)
0250   {
0251     double ecluster = _clusters->getTotalEdep();
0252     double etower = towers->getTotalEdep();
0253     if (ecluster > 0)
0254     {
0255       if (fabs(etower - ecluster) / ecluster > 1e-9)
0256       {
0257         std::cout << "energy conservation violation: ETower: " << etower
0258                   << " ECluster: " << ecluster
0259                   << " diff: " << etower - ecluster << std::endl;
0260       }
0261     }
0262     else
0263     {
0264       if (etower != 0)
0265       {
0266         std::cout << "energy conservation violation: ETower: " << etower
0267                   << " ECluster: " << ecluster << std::endl;
0268       }
0269     }
0270   }
0271   return Fun4AllReturnCodes::EVENT_OK;
0272 }
0273
0274 int RawClusterBuilderGraph::End(PHCompositeNode * /*topNode*/)
0275 {
0276   return Fun4AllReturnCodes::EVENT_OK;
0277 }
0278
0279 void RawClusterBuilderGraph::CreateNodes(PHCompositeNode *topNode)
0280 {
0281   PHNodeIterator iter(topNode);
0282
0283   // Grab the cEMC node
0284   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0285   if (!dstNode)
0286   {
0287     std::cerr << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0288     throw std::runtime_error("Failed to find DST node in EmcRawTowerBuilder::CreateNodes");
0289   }
0290
0291   PHNodeIterator dstiter(dstNode);
0292   PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", detector));
0293   if (!DetNode)
0294   {
0295     DetNode = new PHCompositeNode(detector);
0296     dstNode->addNode(DetNode);
0297   }
0298
0299   _clusters = new RawClusterContainer();
0300   ClusterNodeName = "CLUSTER_" + detector;
0301   PHIODataNode<PHObject> *clusterNode = new PHIODataNode<PHObject>(_clusters, ClusterNodeName, "PHObject");
0302   DetNode->addNode(clusterNode);
0303 }