Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawTowerDeadTowerInterp.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawTowerDeadTowerInterp.cc
2 
3 #include <calobase/TowerInfo.h>
4 #include <calobase/TowerInfoContainer.h>
5 #include <calobase/RawTowerDeadMap.h>
6 #include <calobase/RawTowerDefs.h>
7 #include <calobase/RawTowerGeom.h>
8 #include <calobase/RawTowerGeomContainer.h>
9 
10 #include <fun4all/Fun4AllBase.h>
12 #include <fun4all/SubsysReco.h>
13 
14 #include <phool/PHCompositeNode.h>
15 #include <phool/PHNode.h>
16 #include <phool/PHNodeIterator.h>
17 #include <phool/getClass.h>
18 
19 #include <cassert>
20 #include <cstdlib>
21 #include <exception>
22 #include <iostream>
23 #include <set> // for _Rb_tree_const_iterator
24 #include <stdexcept>
25 #include <string>
26 #include <utility>
27 #include <vector>
28 
30  : SubsysReco(name)
31 {
32 }
33 
35 {
36  PHNodeIterator iter(topNode);
37 
38  // Looking for the DST node
39  PHCompositeNode *dstNode;
40  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode",
41  "DST"));
42  if (!dstNode)
43  {
44  std::cout << Name() << "::" << m_detector << "::" << __PRETTY_FUNCTION__
45  << "DST Node missing, doing nothing." << std::endl;
46  exit(1);
47  }
48 
49  try
50  {
51  CreateNodes(topNode);
52  }
53  catch (std::exception &e)
54  {
55  std::cout << e.what() << std::endl;
57  }
59 }
60 
62 {
63  if (Verbosity())
64  {
65  std::cout << Name() << "::" << m_detector << "::"
66  << "process_event"
67  << "Process event entered" << std::endl;
68  }
69 
70  double recovery_energy = 0;
71  int recoverTower = 0;
72  int deadTowerCnt = 0;
73 
74  if (m_deadTowerMap)
75  {
76  const int eta_bins = m_geometry->get_etabins();
77  const int phi_bins = m_geometry->get_phibins();
78 
79  // assume cylindrical calorimeters for now. Need to add swithc for planary calorimeter
80  assert(eta_bins > 0);
81  assert(phi_bins > 0);
82 
84  {
85  ++deadTowerCnt;
86 
87  if (Verbosity() >= VERBOSITY_MORE)
88  {
89  std::cout << Name() << "::" << m_detector << "::"
90  << "process_event"
91  << " - processing tower " << key;
92  }
93 
94  // check deadmap validity
95  RawTowerGeom *towerGeom = m_geometry->get_tower_geometry(key);
96  if (towerGeom == nullptr)
97  {
98  std::cout << Name() << "::" << m_detector << "::"
99  << "process_event"
100  << " - invalid dead tower ID " << key << std::endl;
101 
102  exit(2);
103  }
104 
105  const int bineta = towerGeom->get_bineta();
106  const int binphi = towerGeom->get_binphi();
107 
108  if (Verbosity() >= VERBOSITY_MORE)
109  {
110  std::cout << " bin " << bineta << "-" << binphi;
111  std::cout << ". Add neighbors: ";
112  }
113 
114  assert(bineta >= 0);
115  assert(bineta <= eta_bins);
116  assert(binphi >= 0);
117  assert(binphi <= phi_bins);
118 
119  // eight neighbors
120  static const std::vector<std::pair<int, int>> neighborIndexs =
121  {{+1, 0}, {+1, +1}, {0, +1}, {-1, +1}, {-1, 0}, {-1, -1}, {0, -1}, {+1, -1}};
122 
123  int n_neighbor = 0;
124  double E_SumNeighbor = 0;
125  for (const std::pair<int, int> &neighborIndex : neighborIndexs)
126  {
127  int ieta = bineta + neighborIndex.first;
128  int iphi = binphi + neighborIndex.second;
129 
130  if (ieta >= eta_bins)
131  {
132  continue;
133  }
134  if (ieta < 0)
135  {
136  continue;
137  }
138 
139  if (iphi >= phi_bins)
140  {
141  iphi -= phi_bins;
142  }
143  if (iphi < 0)
144  {
145  iphi += phi_bins;
146  }
147 
148  assert(ieta >= 0);
149  assert(ieta <= eta_bins);
150  assert(iphi >= 0);
151  assert(iphi <= phi_bins);
152 
153  if (m_deadTowerMap->isDeadTower(ieta, iphi))
154  {
155  continue;
156  }
157 
158  ++n_neighbor;
159 
161  unsigned int towerkey = (ieta << 16U) + iphi;
162  TowerInfo *neighTower = m_calibTowers->get_tower_at_key(towerkey);
163  // RawTower *neighTower = m_calibTowers->getTower(ieta, iphi);
164  if (neighTower == nullptr)
165  {
166  continue;
167  }
168 
169  if (Verbosity() >= VERBOSITY_MORE)
170  {
171  std::cout << neighTower->get_energy() << " (" << ieta << "-" << iphi << "), ";
172  }
173 
174  E_SumNeighbor += neighTower->get_energy();
175 
176  } // for (const pair<int, int> &neighborIndex : neighborIndexs)
177 
178  if (n_neighbor > 0 and E_SumNeighbor != 0)
179  {
180 
181  unsigned int deadtowerkey = (bineta << 16U) + binphi;
182 
183  TowerInfo *deadTower = m_calibTowers->get_tower_at_key(deadtowerkey);
184 
185  // if (deadTower == nullptr)
186  // {
187  // deadTower = new TowerInfo();
188  // }
189  assert(deadTower);
190 
191  deadTower->set_energy(E_SumNeighbor / n_neighbor);
192  // m_calibTowers->AddTower(key, deadTower);
193 
194  recovery_energy += deadTower->get_energy();
195  ++recoverTower;
196 
197  if (Verbosity() >= VERBOSITY_MORE)
198  {
199  std::cout << " -> " << deadTower->get_energy() << " GeV @ etabin: " << bineta << " , phi bin: " << binphi ;
200  }
201 
202  } // if (n_neighbor>0)
203  else
204  {
205  if (Verbosity() >= VERBOSITY_MORE)
206  {
207  std::cout << "No neighbor towers found.";
208  }
209 
210  } // if (n_neighbor>0) -else
211 
212  if (Verbosity() >= VERBOSITY_MORE)
213  {
214  std::cout << std::endl;
215  }
216  }
217 
218  } // if (m_deadTowerMap)
219  else
220  {
221  static bool once = true;
222 
223  if (Verbosity() and once)
224  {
225  once = false;
226 
227  std::cout << Name() << "::" << m_detector << "::"
228  << "process_event"
229  << " - missing dead map node. Do nothing ..."
230  << std::endl;
231  }
232  }
233 
234  if (Verbosity())
235  {
236  std::cout << Name() << "::" << m_detector << "::"
237  << "process_event"
238  << "recovery_energy = " << recovery_energy
239  << " GeV from " << recoverTower << " towers out of total " << deadTowerCnt << " dead towers"
240  << std::endl;
241  }
243 }
244 
246 {
248 }
249 
251 {
252  PHNodeIterator iter(topNode);
253  PHCompositeNode *runNode = static_cast<PHCompositeNode *>(iter.findFirst(
254  "PHCompositeNode", "RUN"));
255  if (!runNode)
256  {
257  std::cout << Name() << "::" << m_detector << "::"
258  << "CreateNodes"
259  << "Run Node missing, doing nothing." << std::endl;
260  throw std::runtime_error(
261  "Failed to find Run node in RawTowerDeadTowerInterp::CreateNodes");
262  }
263 
264  const std::string deadMapName = "DEADMAP_" + m_detector;
265  m_deadTowerMap = findNode::getClass<RawTowerDeadMap>(topNode, deadMapName);
266  if (m_deadTowerMap)
267  {
268  std::cout << Name() << "::" << m_detector << "::"
269  << "CreateNodes"
270  << " use dead map: ";
272  }
273 
274  const std::string geometry_node = "TOWERGEOM_" + m_detector;
275  m_geometry = findNode::getClass<RawTowerGeomContainer>(topNode, geometry_node);
276  if (!m_geometry)
277  {
278  std::cout << Name() << "::" << m_detector << "::"
279  << "CreateNodes"
280  << " " << geometry_node << " Node missing, doing bail out!"
281  << std::endl;
282  throw std::runtime_error(
283  "Failed to find " + geometry_node + " node in RawTowerDeadTowerInterp::CreateNodes");
284  }
285 
286  if (Verbosity() >= 1)
287  {
288  m_geometry->identify();
289  }
290 
291  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst(
292  "PHCompositeNode", "DST"));
293  if (!dstNode)
294  {
295  std::cout << Name() << "::" << m_detector << "::"
296  << "CreateNodes"
297  << "DST Node missing, doing nothing." << std::endl;
298  throw std::runtime_error(
299  "Failed to find DST node in RawTowerDeadTowerInterp::CreateNodes");
300  }
301 
302  const std::string rawTowerNodeName = "TOWERINFO_" + _calib_tower_node_prefix + "_" + m_detector;
303  m_calibTowers = findNode::getClass<TowerInfoContainer>(dstNode, rawTowerNodeName);
304  if (!m_calibTowers)
305  {
306  std::cout << Name() << "::" << m_detector << "::" << __PRETTY_FUNCTION__
307  << " " << rawTowerNodeName << " Node missing, doing bail out!"
308  << std::endl;
309  throw std::runtime_error(
310  "Failed to find " + rawTowerNodeName + " node in RawTowerDeadTowerInterp::CreateNodes");
311  }
312 
313  return;
314 }