ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
RawTowerBuilder.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file RawTowerBuilder.cc
1 #include "RawTowerBuilder.h"
2 
3 #include <calobase/RawTower.h> // for RawTower
4 #include <calobase/RawTowerv1.h>
5 #include <calobase/RawTowerContainer.h>
6 #include <calobase/RawTowerDefs.h> // for encode_towerid
7 #include <calobase/RawTowerGeomContainer.h> // for RawTowerGeomC...
8 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
9 #include <calobase/RawTowerGeom.h> // for RawTowerGeom
10 #include <calobase/RawTowerGeomv1.h>
11 
14 
15 #include <g4detectors/PHG4Cell.h>
18 
19 #include <g4main/PHG4Utils.h>
20 
22 #include <fun4all/SubsysReco.h> // for SubsysReco
23 
24 #include <phool/PHCompositeNode.h>
25 #include <phool/PHIODataNode.h>
26 #include <phool/PHNode.h> // for PHNode
27 #include <phool/PHNodeIterator.h>
28 #include <phool/PHObject.h> // for PHObject
29 #include <phool/getClass.h>
30 #include <phool/phool.h> // for PHWHERE
31 
32 #include <boost/io/ios_state.hpp>
33 
34 #include <cmath> // for fabs, tan, atan2
35 #include <cstdlib> // for exit
36 #include <exception> // for exception
37 #include <iostream>
38 #include <map>
39 #include <stdexcept>
40 #include <utility> // for pair, make_pair
41 
42 using namespace std;
43 
45  : SubsysReco(name)
46  , m_TowerContainer(nullptr)
47  , m_RawTowerGeomContainer(nullptr)
48  , m_Detector("NONE")
49  , m_TowerEnergySrcEnum(kLightYield)
50  , m_CellBinning(PHG4CellDefs::undefined)
51  , m_ChkEnergyConservationFlag(0)
52  , m_NumLayers(-1)
53  , m_NumPhiBins(-1)
54  , m_NumEtaBins(-1)
55  , m_Emin(1e-6)
56  , m_EtaMin(NAN)
57  , m_PhiMin(NAN)
58  , m_EtaStep(NAN)
59  , m_PhiStep(NAN)
60 {
61 }
62 
64 {
65  std::string geonodename = "CYLINDERCELLGEOM_" + m_Detector;
67  PHG4CylinderCellGeomContainer>(topNode, geonodename);
68  if (!cellgeos)
69  {
70  cout << PHWHERE << " " << geonodename
71  << " Node missing, doing nothing." << std::endl;
72  throw std::runtime_error(
73  "Failed to find " + geonodename + " node in RawTowerBuilder::CreateNodes");
74  }
75 
76  // fill the number of layers in the calorimeter
77  m_NumLayers = cellgeos->get_NLayers();
78 
79  PHNodeIterator iter(topNode);
80 
81  // Looking for the DST node
82  PHCompositeNode *dstNode;
83  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
84  if (!dstNode)
85  {
86  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
87  exit(1);
88  }
89 
90  try
91  {
92  CreateNodes(topNode);
93  }
94  catch (std::exception &e)
95  {
96  std::cout << e.what() << std::endl;
97  //exit(1);
98  }
99 
100  if (Verbosity() >= 1)
101  {
102  cout << "RawTowerBuilder::InitRun :";
104  {
105  cout << "save Geant4 energy deposition as the weight of the cells"
106  << endl;
107  }
108  else if (m_TowerEnergySrcEnum == kLightYield)
109  {
110  cout << "save light yield as the weight of the cells" << endl;
111  }
112  }
113 
115 }
116 
118 {
119  if (Verbosity())
120  {
121  std::cout << PHWHERE << "Process event entered" << std::endl;
122  }
123 
124  // get cells
125  std::string cellnodename = "G4CELL_" + m_Detector;
126  PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
127  if (!cells)
128  {
129  cout << PHWHERE << " " << cellnodename
130  << " Node missing, doing nothing." << std::endl;
132  }
133 
134  // loop over all cells in an event
136  PHG4CellContainer::ConstRange cell_range = cells->getCells();
137  for (cell_iter = cell_range.first; cell_iter != cell_range.second; ++cell_iter)
138  {
139  PHG4Cell *cell = cell_iter->second;
140 
141  if (Verbosity() > 2)
142  {
143  std::cout << PHWHERE << " print out the cell:" << std::endl;
144  cell->identify();
145  }
146 
147  // add the energy to the corresponding tower
148  RawTower *tower = nullptr;
149  short int firstpar;
150  short int secondpar;
152  {
155  }
157  {
160  }
161  else
162  {
163  boost::io::ios_flags_saver ifs(cout);
164  cout << "unknown cell binning, implement 0x" << hex << PHG4CellDefs::get_binning(cell->get_cellid()) << endl;
165  exit(1);
166  }
167  tower = m_TowerContainer->getTower(firstpar, secondpar);
168  if (!tower)
169  {
170  tower = new RawTowerv1();
171  tower->set_energy(0);
172  m_TowerContainer->AddTower(firstpar, secondpar, tower);
173  }
174  float cell_weight = 0;
176  {
177  cell_weight = cell->get_edep();
178  }
179  else if (m_TowerEnergySrcEnum == kLightYield)
180  {
181  cell_weight = cell->get_light_yield();
182  }
183 
184  tower->add_ecell(cell->get_cellid(), cell_weight);
185 
187  for (PHG4Cell::ShowerEdepConstIterator shower_iter = range.first;
188  shower_iter != range.second;
189  ++shower_iter)
190  {
191  tower->add_eshower(shower_iter->first, shower_iter->second);
192  }
193 
194  tower->set_energy(tower->get_energy() + cell_weight);
195 
196  if (Verbosity() > 2)
197  {
198  m_RawTowerGeomContainer = findNode::getClass<RawTowerGeomContainer>(topNode, m_TowerGeomNodeName);
199  tower->identify();
200  }
201  }
202  double towerE = 0;
204  {
205  double cellE = cells->getTotalEdep();
206  towerE = m_TowerContainer->getTotalEdep();
207  if (fabs(cellE - towerE) / cellE > 1e-5)
208  {
209  cout << "towerE: " << towerE << ", cellE: " << cellE << ", delta: "
210  << cellE - towerE << endl;
211  }
212  }
213  if (Verbosity())
214  {
215  towerE = m_TowerContainer->getTotalEdep();
216  }
217 
219  if (Verbosity())
220  {
221  cout << "Energy lost by dropping towers with less than " << m_Emin
222  << " GeV energy, lost energy: " << towerE - m_TowerContainer->getTotalEdep()
223  << endl;
227  for (iter = begin_end.first; iter != begin_end.second; ++iter)
228  {
229  iter->second->identify();
230  }
231  }
232 
234 }
235 
237 {
238  PHNodeIterator iter(topNode);
239  PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
240  if (!runNode)
241  {
242  cout << PHWHERE << "Run Node missing, doing nothing." << std::endl;
243  throw std::runtime_error("Failed to find Run node in RawTowerBuilder::CreateNodes");
244  }
245 
246  PHNodeIterator runIter(runNode);
247  PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", m_Detector));
248  if (!RunDetNode)
249  {
250  RunDetNode = new PHCompositeNode(m_Detector);
251  runNode->addNode(RunDetNode);
252  }
253 
255 
256  // get the cell geometry and build up the tower geometry object
257  std::string geonodename = "CYLINDERCELLGEOM_" + m_Detector;
258  PHG4CylinderCellGeomContainer *cellgeos = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, geonodename);
259  if (!cellgeos)
260  {
261  cout << PHWHERE << " " << geonodename
262  << " Node missing, doing nothing." << std::endl;
263  throw std::runtime_error(
264  "Failed to find " + geonodename + " node in RawTowerBuilder::CreateNodes");
265  }
266  m_TowerGeomNodeName = "TOWERGEOM_" + m_Detector;
267  m_RawTowerGeomContainer = findNode::getClass<RawTowerGeomContainer>(topNode, m_TowerGeomNodeName);
269  {
272  RunDetNode->addNode(newNode);
273  }
274  // fill the number of layers in the calorimeter
275  m_NumLayers = cellgeos->get_NLayers();
276 
277  // Create the tower nodes on the tree
280  cellgeos->get_begin_end();
281  int ifirst = 1;
282  int first_layer = -1;
283  PHG4CylinderCellGeom *first_cellgeo = nullptr;
284  double inner_radius = 0;
285  double thickness = 0;
286  for (miter = begin_end.first; miter != begin_end.second; ++miter)
287  {
288  PHG4CylinderCellGeom *cellgeo = miter->second;
289 
290  if (Verbosity())
291  {
292  cellgeo->identify();
293  }
294  thickness += cellgeo->get_thickness();
295  if (ifirst)
296  {
297  first_cellgeo = miter->second;
298  m_CellBinning = cellgeo->get_binning();
299  m_NumPhiBins = cellgeo->get_phibins();
300  m_PhiMin = cellgeo->get_phimin();
301  m_PhiStep = cellgeo->get_phistep();
303  {
304  m_NumEtaBins = cellgeo->get_etabins();
305  m_EtaMin = cellgeo->get_etamin();
306  m_EtaStep = cellgeo->get_etastep();
307  }
309  {
310  m_NumEtaBins = cellgeo->get_zbins(); // bin eta in the same number of z bins
311  }
313  {
314  // use eta definiton for each row of towers
315  m_NumEtaBins = cellgeo->get_etabins();
316  }
317  else
318  {
319  cout << "RawTowerBuilder::CreateNodes::" << Name()
320  << " - Fatal Error - unsupported cell binning method "
321  << m_CellBinning << endl;
322  }
323  inner_radius = cellgeo->get_radius();
324  first_layer = cellgeo->get_layer();
325  ifirst = 0;
326  }
327  else
328  {
329  if (m_CellBinning != cellgeo->get_binning())
330  {
331  cout << "inconsistent binning method from " << m_CellBinning
332  << " layer " << cellgeo->get_layer() << ": "
333  << cellgeo->get_binning() << endl;
334  exit(1);
335  }
336  if (inner_radius > cellgeo->get_radius())
337  {
338  cout << "radius of layer " << cellgeo->get_layer() << " is "
339  << cellgeo->get_radius() << " which smaller than radius "
340  << inner_radius << " of first layer in list " << first_layer
341  << endl;
342  }
343  if (m_NumPhiBins != cellgeo->get_phibins())
344  {
345  cout << "mixing number of phibins, fisrt layer: " << m_NumPhiBins
346  << " layer " << cellgeo->get_layer() << ": "
347  << cellgeo->get_phibins() << endl;
348  exit(1);
349  }
350  if (m_PhiMin != cellgeo->get_phimin())
351  {
352  cout << "mixing number of phimin, fisrt layer: " << m_PhiMin
353  << " layer " << cellgeo->get_layer() << ": "
354  << cellgeo->get_phimin() << endl;
355  exit(1);
356  }
357  if (m_PhiStep != cellgeo->get_phistep())
358  {
359  cout << "mixing phisteps first layer: " << m_PhiStep << " layer "
360  << cellgeo->get_layer() << ": " << cellgeo->get_phistep()
361  << " diff: " << m_PhiStep - cellgeo->get_phistep() << endl;
362  exit(1);
363  }
365  {
366  if (m_NumEtaBins != cellgeo->get_etabins())
367  {
368  cout << "mixing number of EtaBins , first layer: "
369  << m_NumEtaBins << " layer " << cellgeo->get_layer() << ": "
370  << cellgeo->get_etabins() << endl;
371  exit(1);
372  }
373  if (fabs(m_EtaMin - cellgeo->get_etamin()) > 1e-9)
374  {
375  cout << "mixing etamin, fisrt layer: " << m_EtaMin << " layer "
376  << cellgeo->get_layer() << ": " << cellgeo->get_etamin()
377  << " diff: " << m_EtaMin - cellgeo->get_etamin() << endl;
378  exit(1);
379  }
380  if (fabs(m_EtaStep - cellgeo->get_etastep()) > 1e-9)
381  {
382  cout << "mixing eta steps first layer: " << m_EtaStep
383  << " layer " << cellgeo->get_layer() << ": "
384  << cellgeo->get_etastep() << " diff: "
385  << m_EtaStep - cellgeo->get_etastep() << endl;
386  exit(1);
387  }
388  }
389 
391  {
392  if (m_NumEtaBins != cellgeo->get_zbins())
393  {
394  cout << "mixing number of z bins , first layer: " << m_NumEtaBins
395  << " layer " << cellgeo->get_layer() << ": "
396  << cellgeo->get_zbins() << endl;
397  exit(1);
398  }
399  }
400  }
401  }
402  m_RawTowerGeomContainer->set_radius(inner_radius);
405  // m_RawTowerGeomContainer->set_phistep(m_PhiStep);
406  // m_RawTowerGeomContainer->set_phimin(m_PhiMin);
408 
409  if (!first_cellgeo)
410  {
411  cout << "RawTowerBuilder::CreateNodes - ERROR - can not find first layer of cells "
412  << endl;
413 
414  exit(1);
415  }
416 
417  for (int ibin = 0; ibin < first_cellgeo->get_phibins(); ibin++)
418  {
419  const pair<double, double> range = first_cellgeo->get_phibounds(ibin);
420 
422  }
423 
425  {
426  const double r = inner_radius;
427 
428  for (int ibin = 0; ibin < first_cellgeo->get_etabins(); ibin++)
429  {
430  const pair<double, double> range = first_cellgeo->get_etabounds(ibin);
431 
433  }
434 
435  // setup location of all towers
436  for (int iphi = 0; iphi < m_RawTowerGeomContainer->get_phibins(); iphi++)
437  {
438  for (int ieta = 0; ieta < m_RawTowerGeomContainer->get_etabins(); ieta++)
439  {
440  const RawTowerDefs::keytype key =
441  RawTowerDefs::encode_towerid(caloid, ieta, iphi);
442 
443  const double x(r * cos(m_RawTowerGeomContainer->get_phicenter(iphi)));
444  const double y(r * sin(m_RawTowerGeomContainer->get_phicenter(iphi)));
445  const double z(r / tan(PHG4Utils::get_theta(m_RawTowerGeomContainer->get_etacenter(ieta))));
446 
448  if (tg)
449  {
450  if (Verbosity() > 0)
451  {
452  cout << "RawTowerBuilder::CreateNodes - Tower geometry " << key << " already exists" << endl;
453  }
454 
455  if (fabs(tg->get_center_x() - x) > 1e-4)
456  {
457  cout << "RawTowerBuilder::CreateNodes - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
458  << endl;
459 
460  exit(1);
461  }
462  if (fabs(tg->get_center_y() - y) > 1e-4)
463  {
464  cout << "RawTowerBuilder::CreateNodes - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
465  << endl;
466  exit(1);
467  }
468  if (fabs(tg->get_center_z() - z) > 1e-4)
469  {
470  cout << "RawTowerBuilder::CreateNodes - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
471  << endl;
472  exit(1);
473  }
474  }
475  else
476  {
477  if (Verbosity() > 0)
478  {
479  cout << "RawTowerBuilder::CreateNodes - building tower geometry " << key << "" << endl;
480  }
481 
482  tg = new RawTowerGeomv1(key);
483 
484  tg->set_center_x(x);
485  tg->set_center_y(y);
486  tg->set_center_z(z);
488  }
489  }
490  }
491  }
493  {
494  const double r = inner_radius;
495 
496  for (int ibin = 0; ibin < first_cellgeo->get_zbins(); ibin++)
497  {
498  const pair<double, double> z_range = first_cellgeo->get_zbounds(ibin);
499  // const double r = first_cellgeo->get_radius();
500  const double eta1 = -log(tan(atan2(r, z_range.first) / 2.));
501  const double eta2 = -log(tan(atan2(r, z_range.second) / 2.));
502  m_RawTowerGeomContainer->set_etabounds(ibin, make_pair(eta1, eta2));
503  }
504 
505  // setup location of all towers
506  for (int iphi = 0; iphi < m_RawTowerGeomContainer->get_phibins(); iphi++)
507  {
508  for (int ibin = 0; ibin < first_cellgeo->get_zbins(); ibin++)
509  {
510  const RawTowerDefs::keytype key = RawTowerDefs::encode_towerid(caloid, ibin, iphi);
511 
512  const double x(r * cos(m_RawTowerGeomContainer->get_phicenter(iphi)));
513  const double y(r * sin(m_RawTowerGeomContainer->get_phicenter(iphi)));
514  const double z(first_cellgeo->get_zcenter(ibin));
515 
517  if (tg)
518  {
519  if (Verbosity() > 0)
520  {
521  cout << "RawTowerBuilder::CreateNodes - Tower geometry " << key << " already exists" << endl;
522  }
523 
524  if (fabs(tg->get_center_x() - x) > 1e-4)
525  {
526  cout << "RawTowerBuilder::CreateNodes - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
527  << endl;
528 
529  exit(1);
530  }
531  if (fabs(tg->get_center_y() - y) > 1e-4)
532  {
533  cout << "RawTowerBuilder::CreateNodes - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
534  << endl;
535  exit(1);
536  }
537  if (fabs(tg->get_center_z() - z) > 1e-4)
538  {
539  cout << "RawTowerBuilder::CreateNodes - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
540  << endl;
541  exit(1);
542  }
543  }
544  else
545  {
546  if (Verbosity() > 0)
547  {
548  cout << "RawTowerBuilder::CreateNodes - building tower geometry " << key << "" << endl;
549  }
550 
551  tg = new RawTowerGeomv1(key);
552 
553  tg->set_center_x(x);
554  tg->set_center_y(y);
555  tg->set_center_z(z);
557  }
558  }
559  }
560  }
561  else
562  {
563  cout << "RawTowerBuilder::CreateNodes - ERROR - unsupported cell geometry "
564  << m_CellBinning << endl;
565  exit(1);
566  }
567 
568  if (Verbosity() >= 1)
569  {
571  }
572 
573  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
574  if (!dstNode)
575  {
576  cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
577  throw std::runtime_error(
578  "Failed to find DST node in RawTowerBuilder::CreateNodes");
579  }
580 
581  PHNodeIterator dstiter(dstNode);
582  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", m_Detector));
583  if (!DetNode)
584  {
585  DetNode = new PHCompositeNode(m_Detector);
586  dstNode->addNode(DetNode);
587  }
588 
589  // Create the tower nodes on the tree
590  if (m_SimTowerNodePrefix.empty())
591  {
592  // no prefix, consistent with older convention
593  m_TowerNodeName = "TOWER_" + m_Detector;
594  }
595  else
596  {
597  m_TowerNodeName = "TOWER_" + m_SimTowerNodePrefix + "_" + m_Detector;
598  }
599  m_TowerContainer = findNode::getClass<RawTowerContainer>(DetNode, m_TowerNodeName.c_str());
600  if (m_TowerContainer == nullptr)
601  {
602  m_TowerContainer = new RawTowerContainer(caloid);
603 
605  DetNode->addNode(towerNode);
606  }
607 
608  return;
609 }