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