ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawClusterDeadAreaMask.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawClusterDeadAreaMask.cc
2 
3 #include <calobase/RawCluster.h>
4 #include <calobase/RawClusterContainer.h>
5 #include <calobase/RawTower.h>
6 #include <calobase/RawTowerContainer.h>
7 #include <calobase/RawTowerDeadMap.h>
8 #include <calobase/RawTowerGeomContainer.h>
9 
10 #include <fun4all/Fun4AllBase.h>
12 #include <fun4all/SubsysReco.h>
13 
14 #include <phool/getClass.h>
15 #include <phool/PHCompositeNode.h>
16 #include <phool/PHNode.h>
17 #include <phool/PHNodeIterator.h>
18 
19 #include <cassert>
20 #include <cmath>
21 #include <fstream>
22 #include <iostream>
23 #include <map>
24 #include <sstream>
25 #include <stdexcept>
26 #include <string>
27 #include <utility> // for pair
28 #include <vector>
29 
30 using namespace std;
31 
33  : SubsysReco(name)
34  , m_detector("NONE")
35  , m_deadTowerMaskHalfWidth(1.6)
36  , m_rawClusters(nullptr)
37  , m_deadMap(nullptr)
38  , m_calibTowers(nullptr)
39  , m_geometry(nullptr)
40 {
41 }
42 
44 {
45  CreateNodeTree(topNode);
46 
48 }
49 
51 {
52  if (Verbosity() >= VERBOSITY_SOME)
53  {
54  cout << Name() << "::" << m_detector << "::process_event - Entry" << endl;
55  }
56  int nMasked = 0;
57 
58  const int eta_bins = m_geometry->get_etabins();
59  const int phi_bins = m_geometry->get_phibins();
60  assert(eta_bins > 0);
61  assert(phi_bins > 0);
62 
63  //loop over the clusters
66 
67  for (iter = begin_end.first; iter != begin_end.second;)
68  {
69  // RawClusterDefs::keytype key = iter->first;
70  const RawCluster *cluster = iter->second;
71 
72  std::vector<float> toweretas;
73  std::vector<float> towerphis;
74  std::vector<float> towerenergies;
75 
76  //loop over the towers in the cluster
77  RawCluster::TowerConstRange towers = cluster->get_towers();
79 
80  for (toweriter = towers.first;
81  toweriter != towers.second;
82  ++toweriter)
83  {
84  RawTower *tower = m_calibTowers->getTower(toweriter->first);
85  assert(tower);
86 
87  int towereta = tower->get_bineta();
88  int towerphi = tower->get_binphi();
89  double towerenergy = tower->get_energy();
90 
91  //put the etabin, phibin, and energy into the corresponding vectors
92  toweretas.push_back(towereta);
93  towerphis.push_back(towerphi);
94  towerenergies.push_back(towerenergy);
95  }
96 
97  int ntowers = toweretas.size();
98  assert(ntowers >= 1);
99 
100  //loop over the towers to determine the energy
101  //weighted eta and phi position of the cluster
102 
103  float etamult = 0;
104  float etasum = 0;
105  float phimult = 0;
106  float phisum = 0;
107 
108  for (int j = 0; j < ntowers; j++)
109  {
110  float energymult = towerenergies.at(j) * toweretas.at(j);
111  etamult += energymult;
112  etasum += towerenergies.at(j);
113 
114  int phibin = towerphis.at(j);
115 
116  if (phibin - towerphis.at(0) < -phi_bins / 2.0)
117  phibin += phi_bins;
118  else if (phibin - towerphis.at(0) > +phi_bins / 2.0)
119  phibin -= phi_bins;
120  assert(abs(phibin - towerphis.at(0)) <= phi_bins / 2.0);
121 
122  energymult = towerenergies.at(j) * phibin;
123  phimult += energymult;
124  phisum += towerenergies.at(j);
125  }
126 
127  float avgphi = phimult / phisum;
128  float avgeta = etamult / etasum;
129 
130  if (Verbosity() > VERBOSITY_MORE)
131  {
132  cout << Name() << "::" << m_detector << "::process_event - "
133  << "process cluster at average location " << avgeta << "," << avgphi << " : ";
134  cluster->identify();
135  }
136 
137  // rejection if close to a dead tower
138  bool rejecCluster = false;
139 
140  for (int search_eta = ceil(avgeta - m_deadTowerMaskHalfWidth); search_eta <= floor(avgeta + m_deadTowerMaskHalfWidth); ++search_eta)
141  {
142  for (int search_phi = ceil(avgphi - m_deadTowerMaskHalfWidth); search_phi <= floor(avgphi + m_deadTowerMaskHalfWidth); ++search_phi)
143  {
144  int ieta = search_eta;
145  int iphi = search_phi;
146 
147  if (ieta >= eta_bins) continue;
148  if (ieta < 0) continue;
149 
150  if (iphi >= phi_bins) iphi -= phi_bins;
151  if (iphi < 0) iphi += phi_bins;
152 
153  const bool isDead = m_deadMap->isDeadTower(ieta, iphi);
154 
155  // dead tower found in cluster
156  if (Verbosity() > VERBOSITY_MORE)
157  {
158  cout << "\t"
159  << "tower " << ieta
160  << "," << iphi
161  << (isDead ? ": is dead." : "OK")
162  << endl;
163  }
164 
165  if (isDead)
166  {
167  rejecCluster = true;
168  break;
169  }
170  }
171 
172  if (rejecCluster) break;
173  }
174 
175  //container operation
176  if (rejecCluster)
177  {
178  if (Verbosity() > VERBOSITY_MORE)
179  {
180  cout << Name() << "::" << m_detector << "::process_event - "
181  << "reject cluster " << cluster->get_id() << endl;
182  cluster->identify();
183  }
184 
185  ++nMasked;
186  m_rawClusters->getClustersMap().erase(iter++);
187  }
188  else
189  {
190  if (Verbosity() > VERBOSITY_MORE)
191  {
192  cout << Name() << "::" << m_detector << "::process_event - "
193  << "keep cluster " << cluster->get_id() << endl;
194  }
195  ++iter;
196  }
197 
198  } // for (iter = begin_end.first; iter != begin_end.second;)
199 
200  if (Verbosity() >= VERBOSITY_SOME)
201  {
202  cout << Name() << "::" << m_detector << "::process_event - masked "
203  << nMasked << " clusters. Final cluster containers has "
204  << m_rawClusters->size() << " clusters"
205  << endl;
206  }
208 }
209 
211 {
212  PHNodeIterator iter(topNode);
213 
214  //Get the DST Node
215  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
216 
217  //Check that it is there
218  if (!dstNode)
219  {
220  std::cerr << "DST Node missing, quitting" << std::endl;
221  throw std::runtime_error("failed to find DST node in RawClusterDeadAreaMask::CreateNodeTree");
222  }
223 
224  m_rawClusters = findNode::getClass<RawClusterContainer>(topNode, "CLUSTER_" + m_detector);
225  if (!m_rawClusters)
226  {
227  cout << Name() << "::" << m_detector << "::"
228  << "CreateNodeTree "
229  "No "
230  << m_detector << " Cluster Container found while in RawClusterDeadAreaMask, can't proceed!!!" << std::endl;
231  topNode->print();
232  throw std::runtime_error("failed to find CLUSTER node in RawClusterDeadAreaMask::CreateNodeTree");
233  }
234 
235  m_calibTowers = findNode::getClass<RawTowerContainer>(topNode, "TOWER_CALIB_" + m_detector);
236  if (!m_calibTowers)
237  {
238  cout << Name() << "::" << m_detector << "::"
239  << "CreateNodeTree "
240  << "No calibrated " << m_detector << " tower info found while in RawClusterDeadAreaMask, can't proceed!!!" << std::endl;
241  throw std::runtime_error("failed to find TOWER_CALIB node in RawClusterDeadAreaMask::CreateNodeTree");
242  }
243 
244  string towergeomnodename = "TOWERGEOM_" + m_detector;
245  m_geometry = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnodename);
246  if (!m_geometry)
247  {
248  cout << Name() << "::" << m_detector << "::"
249  << "CreateNodeTree"
250  << ": Could not find node " << towergeomnodename << endl;
251  throw std::runtime_error("failed to find TOWERGEOM node in RawClusterDeadAreaMask::CreateNodeTree");
252  }
253 
254  const string deadMapName = "DEADMAP_" + m_detector;
255  m_deadMap = findNode::getClass<RawTowerDeadMap>(topNode, deadMapName);
256  if (m_deadMap)
257  {
258  cout << Name() << "::" << m_detector << "::"
259  << "CreateNodeTree"
260  << " use dead map: ";
261  m_deadMap->identify();
262  }
263 }
264 
266 {
268 }