ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
HcalRawTowerBuilder.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file HcalRawTowerBuilder.cc
1 #include "HcalRawTowerBuilder.h"
2 
3 #include <calobase/RawTower.h> // for RawTower
4 #include <calobase/RawTowerContainer.h>
5 #include <calobase/RawTowerDefs.h> // for convert_name_...
6 #include <calobase/RawTowerGeom.h> // for RawTowerGeom
7 #include <calobase/RawTowerGeomContainer.h> // for RawTowerGeomC...
8 #include <calobase/RawTowerGeomContainer_Cylinderv1.h>
9 #include <calobase/RawTowerGeomv1.h>
10 #include <calobase/RawTowerv1.h>
11 
12 #include <g4detectors/PHG4Cell.h>
16 
17 #include <phparameter/PHParameterInterface.h> // for PHParameterIn...
18 #include <phparameter/PHParameters.h>
19 
20 #include <g4main/PHG4Utils.h>
21 
23 #include <fun4all/Fun4AllServer.h>
24 #include <fun4all/SubsysReco.h> // for SubsysReco
25 
26 #include <pdbcalbase/PdbParameterMapContainer.h>
27 
28 #include <phool/PHCompositeNode.h>
29 #include <phool/PHIODataNode.h>
30 #include <phool/PHNode.h> // for PHNode
31 #include <phool/PHNodeIterator.h>
32 #include <phool/PHObject.h> // for PHObject
33 #include <phool/getClass.h>
34 #include <phool/phool.h> // for PHWHERE
35 
36 #include <TSystem.h>
37 
38 #include <cmath> // for fabs, NAN, cos
39 #include <exception> // for exception
40 #include <filesystem>
41 #include <fstream>
42 #include <iostream>
43 #include <map>
44 #include <stdexcept>
45 #include <utility> // for make_pair, pair
46 
48  : SubsysReco(name)
49  , PHParameterInterface(name)
50 {
52 
53 }
54 
56 {
57  if (m_Detector.empty())
58  {
59  std::cout << PHWHERE
60  << " Detector name not set, use HcalRawTowerBuilder::Detector(string) to set, exiting"
61  << std::endl;
62  gSystem->Exit(1);
63  exit(1);
64  }
65  PHNodeIterator iter(topNode);
66 
67  // Looking for the DST node
68  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
69  if (!dstNode)
70  {
71  std::cout << PHWHERE << "DST Node missing, exiting" << std::endl;
72  gSystem->Exit(1);
73  exit(1);
74  }
75  PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
76  std::string paramnodename = "TOWERPARAM_" + m_Detector;
77 
78  try
79  {
80  CreateNodes(topNode);
81  }
82  catch (std::exception &e)
83  {
84  std::cout << e.what() << std::endl;
85  gSystem->Exit(1);
86  exit(1);
87  }
88 
89  // order first default,
90  // then parameter from g4detector on node tree
91  ReadParamsFromNodeTree(topNode);
92  // then macro setting
94  PHNodeIterator runIter(runNode);
95  PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", m_Detector));
96  if (!RunDetNode)
97  {
98  RunDetNode = new PHCompositeNode(m_Detector);
99  runNode->addNode(RunDetNode);
100  }
101  SaveToNodeTree(RunDetNode, paramnodename);
102  PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
103  std::string geonodename = "TOWERGEO_" + m_Detector;
104 
105  PHNodeIterator parIter(parNode);
106  PHCompositeNode *ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", m_Detector));
107  if (!ParDetNode)
108  {
109  ParDetNode = new PHCompositeNode(m_Detector);
110  parNode->addNode(ParDetNode);
111  }
112  PutOnParNode(ParDetNode, geonodename);
113  m_TowerEnergySrc = get_int_param("tower_energy_source");
114  m_Emin = get_double_param("emin");
115  m_NcellToTower = get_int_param("n_scinti_plates_per_tower");
116  if (!m_TowerDecalFactors.empty())
117  {
119  }
120  if (Verbosity() >= 1)
121  {
122  std::cout << "HcalRawTowerBuilder::InitRun :";
124  {
125  std::cout << "save Geant4 energy deposition in towers" << std::endl;
126  }
127  else if (m_TowerEnergySrc == kLightYield)
128  {
129  std::cout << "save light yield in towers" << std::endl;
130  }
132  {
133  std::cout << "save ionization energy in towers" << std::endl;
134  }
135  else
136  {
137  std::cout << "unknown energy source" << std::endl;
138  }
139  }
140  m_TowerGeomNodeName = "TOWERGEOM_" + m_Detector;
141  m_RawTowerGeom = findNode::getClass<RawTowerGeomContainer>(topNode, m_TowerGeomNodeName);
142  if (!m_RawTowerGeom)
143  {
146  RunDetNode->addNode(newNode);
147  }
150  m_RawTowerGeom->set_radius(innerrad);
151  m_RawTowerGeom->set_thickness(thickness);
154  double geom_ref_radius = innerrad + thickness / 2.;
155  double phistart = 0;
156  for (int i = 0; i < get_int_param(PHG4HcalDefs::n_towers); i++)
157  {
158  double phiend = phistart + 2. * M_PI / get_int_param(PHG4HcalDefs::n_towers);
159  std::pair<double, double> range = std::make_pair(phistart, phiend);
160  phistart = phiend;
161  m_RawTowerGeom->set_phibounds(i, range);
162  }
163  //double etalowbound = -1.1;
164  double etalowbound = -get_double_param("scinti_eta_coverage_neg");
165  for (int i = 0; i < get_int_param("etabins"); i++)
166  {
167  //double etahibound = etalowbound + 2.2 / get_int_param("etabins");
168  double etahibound = etalowbound +
169  (get_double_param("scinti_eta_coverage_neg")+get_double_param("scinti_eta_coverage_pos")) / get_int_param("etabins");
170  std::pair<double, double> range = std::make_pair(etalowbound, etahibound);
171  m_RawTowerGeom->set_etabounds(i, range);
172  etalowbound = etahibound;
173  }
174  for (int iphi = 0; iphi < m_RawTowerGeom->get_phibins(); iphi++)
175  {
176  for (int ieta = 0; ieta < m_RawTowerGeom->get_etabins(); ieta++)
177  {
179 
180  const double x(geom_ref_radius * cos(m_RawTowerGeom->get_phicenter(iphi)));
181  const double y(geom_ref_radius * sin(m_RawTowerGeom->get_phicenter(iphi)));
182  const double z(geom_ref_radius / tan(PHG4Utils::get_theta(m_RawTowerGeom->get_etacenter(ieta))));
183 
185  if (tg)
186  {
187  if (Verbosity() > 0)
188  {
189  std::cout << "HcalRawTowerBuilder::InitRun - Tower geometry " << key << " already exists" << std::endl;
190  }
191 
192  if (fabs(tg->get_center_x() - x) > 1e-4)
193  {
194  std::cout << "HcalRawTowerBuilder::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing x = " << tg->get_center_x() << " and expected x = " << x
195  << std::endl;
196 
198  }
199  if (fabs(tg->get_center_y() - y) > 1e-4)
200  {
201  std::cout << "HcalRawTowerBuilder::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing y = " << tg->get_center_y() << " and expected y = " << y
202  << std::endl;
204  }
205  if (fabs(tg->get_center_z() - z) > 1e-4)
206  {
207  std::cout << "HcalRawTowerBuilder::InitRun - Fatal Error - duplicated Tower geometry " << key << " with existing z= " << tg->get_center_z() << " and expected z = " << z
208  << std::endl;
210  }
211  }
212  else
213  {
214  if (Verbosity() > 0)
215  {
216  std::cout << "HcalRawTowerBuilder::InitRun - building tower geometry " << key << "" << std::endl;
217  }
218 
219  tg = new RawTowerGeomv1(key);
220 
221  tg->set_center_x(x);
222  tg->set_center_y(y);
223  tg->set_center_z(z);
225  }
226  }
227  }
228  if (Verbosity() > 0)
229  {
231  }
232 
233  int m = m_DecalArray[0].size();
234  int n = m_DecalArray.size();
235 
236  for(int i = 0; i < n; i++){
237  for(int j= 0; j < m; j++){
238  m_DecalArray[i][j] = 1.;
239  }}
240 
241  if (!m_DeCalibrationFileName.empty())
242  {
244  {
245  std::ifstream decalibrate_tower;
246  decalibrate_tower.open(m_DeCalibrationFileName, std::ifstream::in);
247  if (decalibrate_tower.is_open())
248  {
249  while (!decalibrate_tower.eof())
250  {
251  int etabin = -1;
252  int phibin = -1;
253  for(int i = 0; i < n ; i++)
254  {
255  for(int j = 0; j < m; j++)
256  {
257  decalibrate_tower >> etabin >> phibin >> m_DecalArray[i][j];
258  if (!std::isfinite(m_DecalArray[i][j]))
259  {
260  std::cout << "Calibration constant at etabin " << etabin
261  << ", phibin " << phibin << " in " << m_DeCalibrationFileName
262  << " is not finite: " << m_DecalArray[i][j] << std::endl;
263  gSystem->Exit(1);
264  exit(1);
265  }
266 
267  }
268  }
269  }
270  decalibrate_tower.close();
271  }
272  }
273  }
274 
276 }
277 
279 {
280  /* decalibration occurs if user supplies a non empty decalMap.txt
281  file, otherwise code will proceed with no de-calibration (as is)
282 */
283 
284  double cell_weight = 0.0;
285  if (Verbosity() > 3)
286  {
287  std::cout << PHWHERE << "Process event entered" << std::endl;
288  }
289 
290  // get cells
291  std::string cellnodename = "G4CELL_" + m_Detector;
292  PHG4CellContainer *slats = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
293  if (!slats)
294  {
295  std::cout << PHWHERE << " Node " << cellnodename
296  << " missing, quitting" << std::endl;
297  gSystem->Exit(1);
298  exit(1);
299  }
300 
301  // loop over all slats in an event
303  PHG4CellContainer::ConstRange cell_range = slats->getCells();
304  for (cell_iter = cell_range.first; cell_iter != cell_range.second;
305  ++cell_iter)
306  {
307  PHG4Cell *cell = cell_iter->second;
308 
310 
312  if (!tower)
313  {
314  tower = new RawTowerv1();
315  tower->set_energy(0.0);
317  }
318 
320  {
321  cell_weight = cell->get_edep();
322  }
323  else if (m_TowerEnergySrc == kLightYield)
324  {
325  cell_weight = cell->get_light_yield();
326  }
328  {
329  cell_weight = cell->get_eion();
330  }
331  else
332  {
333  std::cout << Name() << ": unknown tower energy source "
334  << m_TowerEnergySrc << std::endl;
335  gSystem->Exit(1);
336  exit(1);
337  }
338 
341  tower->add_ecell(cell->get_cellid(), cell_weight);
342 
344  for (PHG4Cell::ShowerEdepConstIterator shower_iter = range.first;
345  shower_iter != range.second;
346  ++shower_iter)
347  {
348  tower->add_eshower(shower_iter->first, shower_iter->second);
349  }
350 
351  tower->set_energy(tower->get_energy() + cell_weight);
352  }
353 
354  double towerE = 0;
356  {
357  double cellE = slats->getTotalEdep();
358  towerE = m_Towers->getTotalEdep();
359  if (fabs(cellE - towerE) / cellE > 1e-5)
360  {
361  std::cout << "towerE: " << towerE << ", cellE: " << cellE << ", delta: "
362  << cellE - towerE << std::endl;
363  }
364  }
365 
366  if (Verbosity())
367  {
368  towerE = m_Towers->getTotalEdep();
369  }
370 
372  if (Verbosity())
373  {
374  std::cout << "Energy lost by dropping towers with less than " << m_Emin
375  << " energy, lost energy: " << towerE - m_Towers->getTotalEdep()
376  << std::endl;
377  m_Towers->identify();
380  for (iter = begin_end.first; iter != begin_end.second; ++iter)
381  {
382  iter->second->identify();
383  }
384  }
385 
387 }
388 
390 {
391  PHNodeIterator iter(topNode);
392  PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
393  if (!runNode)
394  {
395  std::cout << PHWHERE << "Run Node missing, exiting." << std::endl;
396  gSystem->Exit(1);
397  exit(1);
398  }
399  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
400  if (!dstNode)
401  {
402  std::cout << PHWHERE << "DST Node missing, exiting." << std::endl;
403  gSystem->Exit(1);
404  exit(1);
405  }
406 
407  PHNodeIterator dstiter(dstNode);
408  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", m_Detector));
409  if (!DetNode)
410  {
411  DetNode = new PHCompositeNode(m_Detector);
412  dstNode->addNode(DetNode);
413  }
414 
415  // Create the tower nodes on the tree
416  if (m_SimTowerNodePrefix.empty())
417  {
418  // no prefix, consistent with older convension
419  m_TowerNodeName = "TOWER_" + m_Detector;
420  }
421  else
422  {
423  m_TowerNodeName = "TOWER_" + m_SimTowerNodePrefix + "_" + m_Detector;
424  }
425  m_Towers = findNode::getClass<RawTowerContainer>(DetNode, m_TowerNodeName);
426  if (!m_Towers)
427  {
429 
431  DetNode->addNode(towerNode);
432  }
433  return;
434 }
435 
436 short HcalRawTowerBuilder::get_tower_row(const short cellrow) const
437 {
438  short twrrow = cellrow / m_NcellToTower;
439  return twrrow;
440 }
441 
443 {
445  set_default_int_param("tower_energy_source", kLightYield);
447  set_default_int_param("etabins", 24);
448 
449  set_default_double_param("emin", 1.e-6);
452 
453  set_default_double_param("scinti_eta_coverage_neg", 1.1);
454  set_default_double_param("scinti_eta_coverage_pos", 1.1);
455 
456 }
457 
459 {
460  PHParameters *pars = new PHParameters("temp");
461  // we need the number of scintillator plates per tower
462  std::string geonodename = "G4GEOPARAM_" + m_Detector;
463  PdbParameterMapContainer *saveparams = findNode::getClass<PdbParameterMapContainer>(topNode, geonodename);
464  if (!saveparams)
465  {
466  std::cout << "could not find " << geonodename << std::endl;
468  se->Print("NODETREE");
469  return;
470  }
471  pars->FillFrom(saveparams, 0);
476 
477  int nTiles = 2 * pars->get_int_param(PHG4HcalDefs::n_scinti_tiles);
479  if (nTiles <= 0)
480  {
482  set_double_param("scinti_eta_coverage_neg",pars->get_double_param("scinti_eta_coverage_neg"));
483  set_double_param("scinti_eta_coverage_pos",pars->get_double_param("scinti_eta_coverage_pos"));
484  }
485  set_int_param("etabins", nTiles);
486  m_DecalArray.resize(nTiles, std::vector<double> (nPhislices));
487 
488  delete pars;
489  return;
490 }
491 
492 void HcalRawTowerBuilder::set_cell_decal_factor(const int etabin, const int phibin, const double d)
493 {
494  m_DecalArray.at(etabin).at(phibin) = d;
495 }
496 
498 {
499  for (auto iter = m_TowerDecalFactors.begin(); iter != m_TowerDecalFactors.end(); ++iter)
500  {
501  set_tower_decal_factor_real(iter->first.first,iter->first.second,iter->second);
502  }
503 }
504 
505 void HcalRawTowerBuilder::set_tower_decal_factor(const int etabin, const int phibin, const double d)
506 {
507  // since we do not have the number of scintillators per tower at this point
508  // the decal values are cached in m_TowerDecalFactors to be set during the InitRun
509  std::pair<int, int> etaphi = std::make_pair(etabin,phibin);
510  m_TowerDecalFactors[etaphi] = d;
511 }
512 
513 void HcalRawTowerBuilder::set_tower_decal_factor_real(const int etabin, const int phibin, const double d)
514 {
515  for (int i=0; i<m_NcellToTower; i++)
516  {
517  int istart = phibin*m_NcellToTower + i;
518  m_DecalArray.at(etabin).at(istart) = d;
519  }
520 }