ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawTowerCombiner.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawTowerCombiner.cc
1 #include "RawTowerCombiner.h"
2 
3 #include <calobase/RawTower.h>
4 #include <calobase/RawTowerContainer.h>
5 #include <calobase/RawTowerDefs.h>
6 #include <calobase/RawTowerGeomContainer.h>
7 #include <calobase/RawTowerGeomv1.h>
8 #include <calobase/RawTowerv1.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 
20 
21 #include <cassert>
22 #include <cmath>
23 #include <cstdlib>
24 #include <exception>
25 #include <iostream>
26 #include <map>
27 #include <stdexcept>
28 #include <utility>
29 #include <vector>
30 
31 using namespace std;
32 
34  : SubsysReco(name)
35  , //
36  /*std::string*/ _tower_node_prefix("SIM")
37  ,
38  /*unsigned int*/ _n_combine_eta(2)
39  ,
40  /*unsigned int*/ _n_combine_phi(2)
41  ,
42  /*RawTowerContainer**/ _towers(nullptr)
43  ,
44  /*std::string*/ detector("NONE")
45 {
46 }
47 
49 {
50  if (_n_combine_eta == 0)
51  {
52  cout << __PRETTY_FUNCTION__ << " Fatal error _n_combine_eta==0" << endl;
53 
54  exit(1243);
55  }
56 
57  if (_n_combine_phi == 0)
58  {
59  cout << __PRETTY_FUNCTION__ << " Fatal error _n_combine_phi==0" << endl;
60 
61  exit(1243);
62  }
63 
64  PHNodeIterator iter(topNode);
65 
66  // Looking for the DST node
67  PHCompositeNode *dstNode;
68  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode",
69  "DST"));
70  if (!dstNode)
71  {
72  std::cout << __PRETTY_FUNCTION__ << "DST Node missing, doing nothing."
73  << std::endl;
74  exit(1);
75  }
76 
77  try
78  {
79  CreateNodes(topNode);
80  }
81  catch (std::exception &e)
82  {
83  std::cout << e.what() << std::endl;
84  //exit(1);
85  }
86 
88 }
89 
91 {
92  assert(_towers);
93 
94  const double input_e_sum = _towers->getTotalEdep();
95  const double input_n_tower = _towers->size();
96 
97  if (Verbosity())
98  {
99  std::cout << __PRETTY_FUNCTION__ << "Process event entered" << std::endl;
100  }
101 
102  typedef pair<int, int> tower_id_t;
103  typedef map<tower_id_t, RawTower *> new_tower_map_t;
104  new_tower_map_t new_tower_map;
105 
107  for (RawTowerContainer::ConstIterator it = all_towers.first;
108  it != all_towers.second; ++it)
109  {
110  const RawTower *input_tower = it->second;
111  assert(input_tower);
112 
113  const int intput_eta = input_tower->get_bineta();
114  const int intput_phi = input_tower->get_binphi();
115  const int output_eta = get_output_bin_eta(intput_eta);
116  const int output_phi = get_output_bin_phi(intput_phi);
117 
118  RawTower *output_tower = nullptr;
119  new_tower_map_t::iterator it_new = new_tower_map.find(
120  make_pair(output_eta, output_phi));
121 
122  if (it_new == new_tower_map.end())
123  {
124  output_tower = new RawTowerv1(*input_tower);
125  assert(output_tower);
126  new_tower_map[make_pair(output_eta, output_phi)] = output_tower;
127 
128  if (Verbosity() >= VERBOSITY_MORE)
129  {
130  std::cout << __PRETTY_FUNCTION__ << "::" << detector << "::"
131  << " new output tower (prior to tower ID assignment): ";
132  output_tower->identify();
133  }
134  }
135  else
136  {
137  output_tower = it_new->second;
138 
139  output_tower->set_energy(
140  output_tower->get_energy() + input_tower->get_energy());
141 
142  output_tower->set_time(
143  (abs(output_tower->get_energy()) * output_tower->get_time() + abs(input_tower->get_energy()) * input_tower->get_time()) //
144  / (abs(output_tower->get_energy()) + abs(input_tower->get_energy()) + 1e-9) //avoid devide 0
145  );
146 
147  RawTower::CellConstRange cell_range = input_tower->get_g4cells();
148 
149  for (RawTower::CellConstIterator cell_iter = cell_range.first;
150  cell_iter != cell_range.second; ++cell_iter)
151  {
152  output_tower->add_ecell(cell_iter->first, cell_iter->second);
153  }
154 
155  RawTower::ShowerConstRange shower_range =
156  input_tower->get_g4showers();
157 
158  for (RawTower::ShowerConstIterator shower_iter = shower_range.first;
159  shower_iter != shower_range.second; ++shower_iter)
160  {
161  output_tower->add_eshower(shower_iter->first,
162  shower_iter->second);
163  }
164 
165  if (Verbosity() >= VERBOSITY_MORE)
166  {
167  std::cout << __PRETTY_FUNCTION__ << "::" << detector << "::"
168  << " merget into output tower (prior to tower ID assignment) : ";
169  output_tower->identify();
170  }
171  }
172  }
173 
174  // replace content in tower container
175  _towers->Reset();
176 
177  for (new_tower_map_t::iterator it = new_tower_map.begin();
178  it != new_tower_map.end(); ++it)
179  {
180  const int eta = it->first.first;
181  const int phi = it->first.second;
182  RawTower *tower = it->second;
183 
184  _towers->AddTower(eta, phi, tower);
185  }
186 
187  if (Verbosity())
188  {
189  std::cout << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
190  << "input sum energy = " << input_e_sum << " from " << input_n_tower << " towers, merged sum energy = "
191  << _towers->getTotalEdep() << " from " << _towers->size() << " towers" << std::endl;
192  assert(input_e_sum == _towers->getTotalEdep());
193  }
194 
196 }
197 
199 {
201 }
202 
204 {
205  PHNodeIterator iter(topNode);
206  PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst(
207  "PHCompositeNode", "RUN"));
208  if (!runNode)
209  {
210  std::cerr << __PRETTY_FUNCTION__ << "Run Node missing, doing nothing."
211  << std::endl;
212  throw std::runtime_error(
213  "Failed to find Run node in RawTowerCombiner::CreateNodes");
214  }
215 
216  const RawTowerDefs::CalorimeterId caloid =
218 
219  const string iTowerGeomNodeName = "TOWERGEOM_" + detector;
221  RawTowerGeomContainer>(topNode, iTowerGeomNodeName);
222  if (!towergeom)
223  {
224  std::cerr << __PRETTY_FUNCTION__ << " - " << iTowerGeomNodeName
225  << " Node missing, doing nothing." << std::endl;
226  throw std::runtime_error(
227  "Failed to find input tower geometry node in RawTowerCombiner::CreateNodes");
228  }
229  assert(towergeom);
230 
231  const double r = towergeom->get_radius();
232 
233  const int new_phibins = ceil(
234  double(towergeom->get_phibins()) / double(_n_combine_phi));
235  const int new_etabins = ceil(
236  double(towergeom->get_etabins()) / double(_n_combine_eta));
237 
238  typedef std::pair<double, double> bound_t;
239  typedef std::vector<bound_t> bound_map_t;
240 
241  bound_map_t eta_bound_map;
242  bound_map_t phi_bound_map;
243 
244  for (int ibin = 0; ibin < new_phibins; ibin++)
245  {
246  const int first_bin = ibin * _n_combine_phi;
247  assert(first_bin >= 0 && first_bin < towergeom->get_phibins());
248 
249  int last_bin = (ibin + 1) * _n_combine_phi - 1;
250  if (last_bin >= towergeom->get_phibins())
251  last_bin = towergeom->get_phibins();
252 
253  const pair<double, double> range1 = towergeom->get_phibounds(
254  first_bin);
255  const pair<double, double> range2 = towergeom->get_phibounds(
256  last_bin);
257 
258  phi_bound_map.push_back(make_pair(range1.first, range2.second));
259  }
260 
261  for (int ibin = 0; ibin < new_etabins; ibin++)
262  {
263  const int first_bin = ibin * _n_combine_eta;
264  assert(first_bin >= 0 && first_bin < towergeom->get_etabins());
265 
266  int last_bin = (ibin + 1) * _n_combine_eta - 1;
267  if (last_bin >= towergeom->get_etabins())
268  last_bin = towergeom->get_etabins();
269 
270  const pair<double, double> range1 = towergeom->get_etabounds(
271  first_bin);
272  const pair<double, double> range2 = towergeom->get_etabounds(
273  last_bin);
274 
275  eta_bound_map.push_back(make_pair(range1.first, range2.second));
276  }
277 
278  // now update the tower geometry object with the new tower structure.
279  towergeom->Reset();
280 
281  towergeom->set_phibins(new_phibins);
282  towergeom->set_etabins(new_etabins);
283 
284  for (int ibin = 0; ibin < new_phibins; ibin++)
285  {
286  towergeom->set_phibounds(ibin, phi_bound_map[ibin]);
287  }
288  for (int ibin = 0; ibin < new_etabins; ibin++)
289  {
290  towergeom->set_etabounds(ibin, eta_bound_map[ibin]);
291  }
292 
293  // setup location of all towers
294  for (int iphi = 0; iphi < towergeom->get_phibins(); iphi++)
295  for (int ieta = 0; ieta < towergeom->get_etabins(); ieta++)
296  {
297  RawTowerGeomv1 *tg = new RawTowerGeomv1(
298  RawTowerDefs::encode_towerid(caloid, ieta, iphi));
299 
300  CLHEP::Hep3Vector tower_pos;
301  tower_pos.setRhoPhiEta(r, towergeom->get_phicenter(iphi), towergeom->get_etacenter(ieta));
302 
303  tg->set_center_x(tower_pos.x());
304  tg->set_center_y(tower_pos.y());
305  tg->set_center_z(tower_pos.z());
306 
307  towergeom->add_tower_geometry(tg);
308  }
309  if (Verbosity() >= VERBOSITY_SOME)
310  {
311  towergeom->identify();
312  }
313 
314  const string input_TowerNodeName = "TOWER_" + _tower_node_prefix + "_" + detector;
315  _towers = findNode::getClass<RawTowerContainer>(topNode,
316  input_TowerNodeName);
317  if (!_towers)
318  {
319  std::cerr << Name() << "::" << detector << "::" << __PRETTY_FUNCTION__
320  << " " << input_TowerNodeName << " Node missing, doing bail out!"
321  << std::endl;
322  throw std::runtime_error(
323  "Failed to find " + input_TowerNodeName + " node in RawTowerCombiner::CreateNodes");
324  }
325 
326  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst(
327  "PHCompositeNode", "DST"));
328  if (!dstNode)
329  {
330  std::cerr << __PRETTY_FUNCTION__ << "DST Node missing, doing nothing."
331  << std::endl;
332  throw std::runtime_error(
333  "Failed to find DST node in RawTowerCombiner::CreateNodes");
334  }
335 
336  return;
337 }