ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4HcalCellReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4HcalCellReco.cc
1 #include "PHG4HcalCellReco.h"
2 
3 #include "PHG4Cell.h" // for PHG4Cell
4 #include "PHG4CellContainer.h"
5 #include "PHG4CellDefs.h" // for genkey, keytype
6 #include "PHG4Cellv1.h"
7 
8 #include <phparameter/PHParameterInterface.h> // for PHParameterInterface
9 
10 #include <g4main/PHG4Hit.h>
12 #include <g4main/PHG4HitDefs.h> // for hit_idbits
13 
15 #include <fun4all/SubsysReco.h> // for SubsysReco
16 
17 #include <phool/PHCompositeNode.h>
18 #include <phool/PHIODataNode.h>
19 #include <phool/PHNode.h> // for PHNode
20 #include <phool/PHNodeIterator.h>
21 #include <phool/PHObject.h> // for PHObject
22 #include <phool/getClass.h>
23 #include <phool/phool.h> // for PHWHERE
24 
25 #include <TSystem.h>
26 
27 #include <array> // for array, array<>::value_...
28 #include <cmath>
29 #include <cstdlib>
30 #include <iostream>
31 #include <map> // for _Rb_tree_const_iterator
32 #include <sstream>
33 #include <utility> // for pair
34 
35 // for hcal dimension
36 #define ROWDIM 320
37 #define COLUMNDIM 27
38 
39 static std::array<std::array<PHG4Cell *, COLUMNDIM>, ROWDIM> slatarray = {{{nullptr}}};
40 
42  : SubsysReco(name)
43  , PHParameterInterface(name)
44 {
46 }
47 
49 {
50  PHNodeIterator iter(topNode);
51 
52  // Looking for the DST node
53  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
54  if (!dstNode)
55  {
56  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
57  exit(1);
58  }
59  PHCompositeNode *runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
60  PHCompositeNode *parNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "PAR"));
61 
62  std::string paramnodename = "G4CELLPARAM_" + detector;
63  std::string geonodename = "G4CELLGEO_" + detector;
64  hitnodename = "G4HIT_" + detector;
65  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
66  if (!g4hit)
67  {
68  std::cout << Name() << " Could not locate G4HIT node " << hitnodename << std::endl;
69  topNode->print();
70  gSystem->Exit(1);
71  exit(1);
72  }
73  cellnodename = "G4CELL_" + detector;
74  PHG4CellContainer *slats = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
75  if (!slats)
76  {
77  PHNodeIterator dstiter(dstNode);
78  PHCompositeNode *DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", detector));
79  if (!DetNode)
80  {
81  DetNode = new PHCompositeNode(detector);
82  dstNode->addNode(DetNode);
83  }
84  slats = new PHG4CellContainer();
85  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(slats, cellnodename, "PHObject");
86  DetNode->addNode(newNode);
87  }
89  // save this to the run wise tree to store on DST
90  PHNodeIterator runIter(runNode);
91  PHCompositeNode *RunDetNode = dynamic_cast<PHCompositeNode *>(runIter.findFirst("PHCompositeNode", detector));
92  if (!RunDetNode)
93  {
94  RunDetNode = new PHCompositeNode(detector);
95  runNode->addNode(RunDetNode);
96  }
97  SaveToNodeTree(RunDetNode, paramnodename);
98  // save this to the parNode for use
99  PHNodeIterator parIter(parNode);
100  PHCompositeNode *ParDetNode = dynamic_cast<PHCompositeNode *>(parIter.findFirst("PHCompositeNode", detector));
101  if (!ParDetNode)
102  {
103  ParDetNode = new PHCompositeNode(detector);
104  parNode->addNode(ParDetNode);
105  }
106  PutOnParNode(ParDetNode, geonodename);
107  tmin = get_double_param("tmin");
108  tmax = get_double_param("tmax");
110 }
111 
113 {
114  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
115  if (!g4hit)
116  {
117  std::cout << "Could not locate g4 hit node " << hitnodename << std::endl;
118  exit(1);
119  }
120  PHG4CellContainer *slats = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
121  if (!slats)
122  {
123  std::cout << "could not locate cell node " << cellnodename << std::endl;
124  exit(1);
125  }
127  {
128  int maxcolumn = 24;
129  int maxrow = 320;
130  if (detector == "HCALIN")
131  {
132  maxrow = 256;
133  }
134  for (int icolumn = 0 ; icolumn < maxcolumn; icolumn++)
135  {
136  for (int irow = 0; irow < maxrow; irow++)
137  {
139  PHG4Cell *cell = new PHG4Cellv1(key);
140  cell->add_edep(m_FixedEnergy);
141  cell->add_eion(m_FixedEnergy);
143  slats->AddCell(cell);
144  }
145  }
147  }
149  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
150  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
151  {
152  if (hiter->second->get_t(0) > tmax) continue;
153  if (hiter->second->get_t(1) < tmin) continue;
154  short icolumn = hiter->second->get_scint_id();
155  int introw = (hiter->second->get_hit_id() >> PHG4HitDefs::hit_idbits);
156  if (introw >= ROWDIM || introw < 0)
157  {
158  std::cout << __PRETTY_FUNCTION__ << " row " << introw
159  << " exceed array size: " << ROWDIM
160  << " adjust ROWDIM and recompile" << std::endl;
161  exit(1);
162  }
163  // after checking for size of introw so we do not run into
164  // overflow issues, put this into the short we want later
165  short irow = introw;
166  if (icolumn >= COLUMNDIM || icolumn < 0)
167  {
168  std::cout << __PRETTY_FUNCTION__ << " column: " << icolumn
169  << " exceed array size: " << COLUMNDIM
170  << " adjust COLUMNDIM and recompile" << std::endl;
171  exit(1);
172  }
173 
174  if (!slatarray[irow][icolumn])
175  {
176  // hcal has no layers so far, I do not want to make an expensive
177  // call to the g4hits to find that out use 0 as layer number
179  slatarray[irow][icolumn] = new PHG4Cellv1(key);
180  }
181  slatarray[irow][icolumn]->add_edep(hiter->second->get_edep());
182  slatarray[irow][icolumn]->add_eion(hiter->second->get_eion());
183  slatarray[irow][icolumn]->add_light_yield(hiter->second->get_light_yield());
184  slatarray[irow][icolumn]->add_edep(hiter->first, hiter->second->get_edep());
185  slatarray[irow][icolumn]->add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep());
186  } // end loop over g4hits
187  int nslathits = 0;
188  for (int irow = 0; irow < ROWDIM; irow++)
189  {
190  for (int icolumn = 0; icolumn < COLUMNDIM; icolumn++)
191  {
192  if (slatarray[irow][icolumn])
193  {
194  slats->AddCell(slatarray[irow][icolumn]);
195  slatarray[irow][icolumn] = nullptr;
196  nslathits++;
197  }
198  }
199  }
200  if (Verbosity() > 0)
201  {
202  std::cout << Name() << ": found " << nslathits << " slats with energy deposition" << std::endl;
203  }
204 
206  {
207  CheckEnergy(topNode);
208  }
210 }
211 
213 {
215 }
216 
218 {
219  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
220  PHG4CellContainer *slats = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
221  double sum_energy_g4hit = 0.;
222  double sum_energy_cells = 0.;
223  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
225  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
226  {
227  sum_energy_g4hit += hiter->second->get_edep();
228  }
229  PHG4CellContainer::ConstRange cell_begin_end = slats->getCells();
231  for (citer = cell_begin_end.first; citer != cell_begin_end.second; ++citer)
232  {
233  sum_energy_cells += citer->second->get_edep();
234  }
235  // the fractional eloss for particles traversing eta bins leads to minute rounding errors
236  if (fabs(sum_energy_cells - sum_energy_g4hit) / sum_energy_g4hit > 1e-6)
237  {
238  std::cout << "hint: timing cuts tmin/tmax will do this to you" << std::endl;
239  std::cout << "energy mismatch between cells: " << sum_energy_cells
240  << " and hits: " << sum_energy_g4hit
241  << " diff sum(cells) - sum(hits): " << sum_energy_cells - sum_energy_g4hit
242  << std::endl;
243  return -1;
244  }
245  else
246  {
247  if (Verbosity() > 0)
248  {
249  std::cout << Name() << ": total energy for this event: " << sum_energy_g4hit << " GeV" << std::endl;
250  }
251  }
252  return 0;
253 }
254 
256 {
257  set_default_double_param("tmax", 60.0);
258  set_default_double_param("tmin", -20.0); // collision has a timing spread around the triggered event. Accepting negative time too.
259  return;
260 }
261 
262 void PHG4HcalCellReco::set_timing_window(const double tmi, const double tma)
263 {
264  set_double_param("tmin", tmi);
265  set_double_param("tmax", tma);
266 }