ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4ForwardCalCellReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4ForwardCalCellReco.cc
2 
3 #include <g4detectors/PHG4CylinderCell.h> // for PHG4CylinderCell
6 
7 #include <g4main/PHG4Hit.h>
9 
11 #include <fun4all/SubsysReco.h> // for SubsysReco
12 
13 #include <phool/PHCompositeNode.h>
14 #include <phool/PHIODataNode.h>
15 #include <phool/PHNode.h> // for PHNode
16 #include <phool/PHNodeIterator.h>
17 #include <phool/PHObject.h> // for PHObject
18 #include <phool/getClass.h>
19 #include <phool/phool.h> // for PHWHERE
20 
21 #include <cmath>
22 #include <cstdlib>
23 #include <cstring> // for memset
24 #include <iostream>
25 #include <sstream>
26 
27 using namespace std;
28 
30  : SubsysReco(name)
31  , chkenergyconservation(0)
32  , tmin_default(0.0)
33  , // ns
34  tmax_default(60.0)
35  , // ns
36  tmin_max()
37 {
38  memset(nbins, 0, sizeof(nbins));
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", "DST"));
48  if (!dstNode)
49  {
50  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
51  exit(1);
52  }
53  hitnodename = "G4HIT_" + detector;
54  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
55  if (!g4hit)
56  {
57  cout << "Could not locate g4 hit node " << hitnodename << endl;
58  exit(1);
59  }
60  cellnodename = "G4CELL_" + detector;
61  PHG4CylinderCellContainer *cells = findNode::getClass<PHG4CylinderCellContainer>(topNode, cellnodename);
62  if (!cells)
63  {
64  PHNodeIterator dstiter(dstNode);
65  PHCompositeNode *DetNode =
66  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode",
67  detector));
68  if (!DetNode)
69  {
70  DetNode = new PHCompositeNode(detector);
71  dstNode->addNode(DetNode);
72  }
73  cells = new PHG4CylinderCellContainer();
74  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(cells, cellnodename.c_str(), "PHObject");
75  DetNode->addNode(newNode);
76  }
77 
79 }
80 
82 {
83  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
84  if (!g4hit)
85  {
86  cout << "Could not locate g4 hit node " << hitnodename << endl;
87  exit(1);
88  }
89  PHG4CylinderCellContainer *cells = findNode::getClass<PHG4CylinderCellContainer>(topNode, cellnodename);
90  if (!cells)
91  {
92  cout << "could not locate cell node " << cellnodename << endl;
93  exit(1);
94  }
95 
97  pair<PHG4HitContainer::LayerIter, PHG4HitContainer::LayerIter> layer_begin_end = g4hit->getLayers();
98  for (layer = layer_begin_end.first; layer != layer_begin_end.second; ++layer)
99  {
101  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(*layer);
102  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
103  {
104  // checking ADC timing integration window cut
105  if (hiter->second->get_t(0) > tmax_default) continue;
106  if (hiter->second->get_t(1) < tmin_default) continue;
107 
108  // only hits that deposited energy (or geantinos)
109  if (hiter->second->get_edep() <= 0)
110  continue;
111 
112  unsigned int key = (hiter->second->get_index_j() << 16) + hiter->second->get_index_k();
113  if (celllist.find(key) == celllist.end())
114  {
115  celllist[key] = new PHG4CylinderCellv3();
116  celllist[key]->set_layer(*layer);
117  celllist[key]->set_j_index(hiter->second->get_index_j());
118  celllist[key]->set_k_index(hiter->second->get_index_k());
119  celllist[key]->set_l_index(hiter->second->get_index_l());
120  }
121 
122  celllist[key]->add_edep(hiter->first, hiter->second->get_edep(), hiter->second->get_light_yield());
123  celllist[key]->add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep());
124  }
125 
126  int numcells = 0;
127  for (map<unsigned int, PHG4CylinderCell *>::const_iterator mapiter = celllist.begin(); mapiter != celllist.end(); ++mapiter)
128  {
129  cells->AddCylinderCellSpecifyKey(mapiter->first, mapiter->second);
130  numcells++;
131  }
132  celllist.clear();
133  if (Verbosity() > 0)
134  {
135  cout << Name() << ": found " << numcells << " eta/slat cells with energy deposition" << endl;
136  }
137  }
138 
140  {
141  CheckEnergy(topNode);
142  }
144 }
145 
147 {
149 }
150 
152 {
153  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
154  PHG4CylinderCellContainer *cells = findNode::getClass<PHG4CylinderCellContainer>(topNode, cellnodename);
155  double sum_energy_g4hit = 0.;
156  double sum_energy_cells = 0.;
157  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits();
159  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
160  {
161  sum_energy_g4hit += hiter->second->get_edep();
162  }
163  PHG4CylinderCellContainer::ConstRange cell_begin_end = cells->getCylinderCells();
165  for (citer = cell_begin_end.first; citer != cell_begin_end.second; ++citer)
166  {
167  sum_energy_cells += citer->second->get_edep();
168  }
169  // the fractional eloss for particles traversing eta bins leads to minute rounding errors
170  if (fabs(sum_energy_cells - sum_energy_g4hit) / sum_energy_g4hit > 1e-6)
171  {
172  cout << "energy mismatch between cells: " << sum_energy_cells
173  << " and hits: " << sum_energy_g4hit
174  << " diff sum(cells) - sum(hits): " << sum_energy_cells - sum_energy_g4hit
175  << endl;
176  return -1;
177  }
178  else
179  {
180  if (Verbosity() > 0)
181  {
182  cout << Name() << ": total energy for this event: " << sum_energy_g4hit << " GeV" << endl;
183  }
184  }
185  return 0;
186 }