Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4BlockCellReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4BlockCellReco.cc
1 #include "PHG4BlockCellReco.h"
2 #include "PHG4BlockCellGeom.h"
3 #include "PHG4BlockCellGeomContainer.h"
4 #include "PHG4BlockGeom.h"
5 #include "PHG4BlockGeomContainer.h"
6 #include "PHG4Cell.h" // for PHG4Cell, PHG...
7 #include "PHG4CellContainer.h"
8 #include "PHG4Cellv1.h"
9 
10 #include "PHG4CellDefs.h"
11 
12 #include <phparameter/PHParameterContainerInterface.h> // for PHParameterCo...
13 #include <phparameter/PHParametersContainer.h>
14 
15 #include <g4main/PHG4Hit.h>
16 #include <g4main/PHG4HitContainer.h>
17 #include <g4main/PHG4Utils.h>
18 
19 #include <fun4all/Fun4AllReturnCodes.h>
20 #include <fun4all/SubsysReco.h> // for SubsysReco
21 
22 #include <phool/PHCompositeNode.h>
23 #include <phool/PHIODataNode.h>
24 #include <phool/PHNode.h> // for PHNode
25 #include <phool/PHNodeIterator.h>
26 #include <phool/PHObject.h> // for PHObject
27 #include <phool/getClass.h>
28 #include <phool/phool.h> // for PHWHERE
29 
30 #include <cmath>
31 #include <cstdlib>
32 #include <iostream>
33 #include <sstream>
34 #include <vector> // for vector
35 
36 using namespace std;
37 
38 static vector<PHG4Cell *> cellptarray;
39 
40 PHG4BlockCellReco::PHG4BlockCellReco(const string &name)
41  : SubsysReco(name)
42  , PHParameterContainerInterface(name)
43  , sum_energy_g4hit(0)
44  , chkenergyconservation(0)
45 {
46  SetDefaultParameters();
47 }
48 
49 int PHG4BlockCellReco::ResetEvent(PHCompositeNode * /*topNode*/)
50 {
51  sum_energy_g4hit = 0.;
52  return Fun4AllReturnCodes::EVENT_OK;
53 }
54 
55 int PHG4BlockCellReco::InitRun(PHCompositeNode *topNode)
56 {
57  PHNodeIterator iter(topNode);
58 
59  // Looking for the DST node
60  PHCompositeNode *dstNode;
61  dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
62  if (!dstNode)
63  {
64  cout << Name() << " DST Node missing, doing nothing." << std::endl;
65  exit(1);
66  }
67  PHNodeIterator dstiter(dstNode);
68  PHCompositeNode *DetNode =
69  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode",
70  detector));
71  if (!DetNode)
72  {
73  DetNode = new PHCompositeNode(detector);
74  dstNode->addNode(DetNode);
75  }
76 
77  PHCompositeNode *runNode;
78  runNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "RUN"));
79  if (!runNode)
80  {
81  cout << Name() << "DST Node missing, doing nothing." << endl;
82  exit(1);
83  }
84  PHNodeIterator runiter(runNode);
85  PHCompositeNode *RunDetNode =
86  dynamic_cast<PHCompositeNode *>(runiter.findFirst("PHCompositeNode",
87  detector));
88  if (!RunDetNode)
89  {
90  RunDetNode = new PHCompositeNode(detector);
91  runNode->addNode(RunDetNode);
92  }
93 
94  hitnodename = "G4HIT_" + detector;
95  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
96  if (!g4hit)
97  {
98  cout << Name() << " Could not locate g4 hit node " << hitnodename << endl;
99  exit(1);
100  }
101 
102  cellnodename = "G4CELL_" + detector;
103  PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
104  if (!cells)
105  {
106  cells = new PHG4CellContainer();
107  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(cells, cellnodename, "PHObject");
108  DetNode->addNode(newNode);
109  }
110 
111  geonodename = "BLOCKGEOM_" + detector;
112  PHG4BlockGeomContainer *geo = findNode::getClass<PHG4BlockGeomContainer>(topNode, geonodename);
113  if (!geo)
114  {
115  cout << Name() << " Could not locate geometry node " << geonodename << endl;
116  exit(1);
117  }
118 
119  if (Verbosity() > 0)
120  {
121  geo->identify();
122  }
123 
124  seggeonodename = "BLOCKCELLGEOM_" + detector;
125  PHG4BlockCellGeomContainer *seggeo = findNode::getClass<PHG4BlockCellGeomContainer>(topNode, seggeonodename);
126  if (!seggeo)
127  {
128  seggeo = new PHG4BlockCellGeomContainer();
129  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(seggeo, seggeonodename, "PHObject");
130  runNode->addNode(newNode);
131  }
132  GetParamsContainerModify()->set_name(detector);
133 
134  UpdateParametersWithMacro();
135 
136  map<int, PHG4BlockGeom *>::const_iterator miter;
137  pair<map<int, PHG4BlockGeom *>::const_iterator, map<int, PHG4BlockGeom *>::const_iterator> begin_end = geo->get_begin_end();
138  map<int, std::pair<double, double> >::iterator sizeiter;
139  for (miter = begin_end.first; miter != begin_end.second; ++miter)
140  {
141  PHG4BlockGeom *layergeom = miter->second;
142  int layer = layergeom->get_layer();
143  if (!ExistDetid(layer))
144  {
145  cout << Name() << ": No parameters for detid/layer " << layer
146  << ", hits from this detid/layer will not be accumulated into cells" << endl;
147  continue;
148  }
149  implemented_detid.insert(layer);
150  double radius = sqrt(pow(layergeom->get_center_x(), 2) + pow(layergeom->get_center_y(), 2));
151  double width = layergeom->get_size_x();
152  double length_in_z = layergeom->get_size_z();
153  double zmin = layergeom->get_center_z() - length_in_z / 2.;
154  double zmax = zmin + length_in_z;
155  set_size(layer, get_double_param(layer, "deltaeta"), get_double_param(layer, "deltax"), PHG4CellDefs::etaphibinning);
156  tmin_max.insert(std::make_pair(layer, std::make_pair(get_double_param(layer, "tmin"), get_double_param(layer, "tmax"))));
157  sizeiter = cell_size.find(layer);
158  if (sizeiter == cell_size.end())
159  {
160  cout << Name() << "no cell sizes for layer " << layer << endl;
161  exit(1);
162  }
163 
164  // create geo object and fill with variables common to all binning methods
165  PHG4BlockCellGeom *layerseggeo = new PHG4BlockCellGeom();
166  layerseggeo->set_layer(layergeom->get_layer());
167  // layerseggeo->set_radius(layergeom->get_radius());
168  // layerseggeo->set_thickness(layergeom->get_thickness());
169 
170  if (binning[layer] == PHG4CellDefs::etaphibinning)
171  {
172  // calculate eta at radius+ thickness (outer radius)
173  // length via eta coverage is calculated using the outer radius
174  double etamin = PHG4Utils::get_eta(radius + 0.5 * layergeom->get_size_y(), zmin);
175  double etamax = PHG4Utils::get_eta(radius + 0.5 * layergeom->get_size_y(), zmax);
176  zmin_max[layer] = make_pair(etamin, etamax);
177  double etastepsize = (sizeiter->second).first;
178  double d_etabins;
179  double fract = modf((etamax - etamin) / etastepsize, &d_etabins);
180  if (fract != 0)
181  {
182  d_etabins++;
183  }
184  etastepsize = (etamax - etamin) / d_etabins;
185  (sizeiter->second).first = etastepsize;
186  int etabins = d_etabins;
187  double etahi = etamin + etastepsize;
188  for (int i = 0; i < etabins; i++)
189  {
190  if (etahi > (etamax + 1.e-6)) // etahi is a tiny bit larger due to numerical uncertainties
191  {
192  cout << "etahi: " << etahi << ", etamax: " << etamax << endl;
193  }
194  etahi += etastepsize;
195  }
196 
197  double xmin = -layergeom->get_width() / 2.;
198  //double xmax = -xmin;
199  double xstepsize = (sizeiter->second).second;
200  double d_xbins;
201  fract = modf(width / xstepsize, &d_xbins);
202  if (fract != 0)
203  {
204  d_xbins++;
205  }
206 
207  xstepsize = width / d_xbins;
208  (sizeiter->second).second = xstepsize;
209  int xbins = d_xbins;
210  // double xlow = xmin;
211  // double xhi = xlow + xstepsize;
212 
213  // for (int i = 0; i < xbins; i++)
214  // {
215  // if (xhi > xmax)
216  // {
217  // cout << "xhi: " << xhi << ", xmax: " << xmax << endl;
218  // }
219  // xlow = xhi;
220  // xhi += xstepsize;
221  // }
222 
223  pair<int, int> x_z_bin = make_pair(xbins, etabins);
224  n_x_z_bins[layer] = x_z_bin;
225  layerseggeo->set_binning(PHG4CellDefs::etaphibinning);
226  layerseggeo->set_etabins(etabins);
227  layerseggeo->set_etamin(etamin);
228  layerseggeo->set_etastep(etastepsize);
229  layerseggeo->set_xmin(xmin);
230  layerseggeo->set_xbins(xbins);
231  layerseggeo->set_xstep(xstepsize);
232  xstep[layer] = xstepsize;
233  etastep[layer] = etastepsize;
234  }
235 
236  // add geo object filled by different binning methods
237  seggeo->AddLayerCellGeom(layerseggeo);
238  if (Verbosity() > 1)
239  {
240  layerseggeo->identify();
241  }
242  }
243  string nodename = "G4CELLPARAM_" + GetParamsContainer()->Name();
244  SaveToNodeTree(RunDetNode, nodename);
245  return Fun4AllReturnCodes::EVENT_OK;
246 }
247 
248 int PHG4BlockCellReco::process_event(PHCompositeNode *topNode)
249 {
250  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename);
251  if (!g4hit)
252  {
253  cout << "Could not locate g4 hit node " << hitnodename << endl;
254  exit(1);
255  }
256  PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
257  if (!cells)
258  {
259  cout << "could not locate cell node " << cellnodename << endl;
260  exit(1);
261  }
262 
263  PHG4BlockCellGeomContainer *seggeo = findNode::getClass<PHG4BlockCellGeomContainer>(topNode, seggeonodename);
264  if (!seggeo)
265  {
266  cout << "could not locate geo node " << seggeonodename << endl;
267  exit(1);
268  }
269 
270  PHG4HitContainer::LayerIter layer;
271  pair<PHG4HitContainer::LayerIter, PHG4HitContainer::LayerIter> layer_begin_end = g4hit->getLayers();
272  // cout << "number of layers: " << g4hit->num_layers() << endl;
273  // cout << "number of hits: " << g4hit->size() << endl;
274  // for (layer = layer_begin_end.first; layer != layer_begin_end.second; layer++)
275  // {
276  // cout << "layer number: " << *layer << endl;
277  // }
278 
279  for (layer = layer_begin_end.first; layer != layer_begin_end.second; layer++)
280  {
281  // only handle layers/detector ids which have parameters set
282  if (implemented_detid.find(*layer) == implemented_detid.end())
283  {
284  continue;
285  }
286  PHG4HitContainer::ConstIterator hiter;
287  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(*layer);
288  PHG4BlockCellGeom *geo = seggeo->GetLayerCellGeom(*layer);
289  int nxbins = n_x_z_bins[*layer].first;
290  int nzbins = n_x_z_bins[*layer].second;
291  unsigned int nbins = nxbins * nzbins;
292 
293  if (cellptarray.size() < nbins)
294  {
295  cellptarray.resize(nbins, nullptr);
296  }
297 
298  // ------- eta/x binning ------------------------------------------------------------------------
299  if (binning[*layer] == PHG4CellDefs::etaphibinning)
300  {
301  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
302  {
303  // checking ADC timing integration window cut
304  if (hiter->second->get_t(0) > tmin_max[*layer].second) continue;
305  if (hiter->second->get_t(1) < tmin_max[*layer].first) continue;
306 
307  pair<double, double> etax[2];
308  double xbin[2];
309  double etabin[2];
310  for (int i = 0; i < 2; i++)
311  {
312  etax[i] = PHG4Utils::get_etaphi(hiter->second->get_x(i), hiter->second->get_y(i), hiter->second->get_z(i));
313  etabin[i] = geo->get_etabin(etax[i].first);
314  xbin[i] = geo->get_xbin(etax[i].second);
315  }
316 
317  // check bin range
318  if (xbin[0] < 0 || xbin[0] >= nxbins || xbin[1] < 0 || xbin[1] >= nxbins)
319  {
320  continue;
321  }
322  if (etabin[0] < 0 || etabin[0] >= nzbins || etabin[1] < 0 || etabin[1] >= nzbins)
323  {
324  continue;
325  }
326 
327  if (etabin[0] < 0)
328  {
329  if (Verbosity() > 0)
330  {
331  hiter->second->identify();
332  }
333  continue;
334  }
335 
336  sum_energy_g4hit += hiter->second->get_edep();
337  int intxbin = xbin[0];
338  int intetabin = etabin[0];
339  int intxbinout = xbin[1];
340  int intetabinout = etabin[1];
341 
342  // Determine all fired cells
343 
344  double ax = (etax[0]).second; // x
345  double ay = (etax[0]).first; // eta
346  double bx = (etax[1]).second;
347  double by = (etax[1]).first;
348  if (intxbin > intxbinout)
349  {
350  int tmp = intxbin;
351  intxbin = intxbinout;
352  intxbinout = tmp;
353  }
354 
355  if (intetabin > intetabinout)
356  {
357  int tmp = intetabin;
358  intetabin = intetabinout;
359  intetabinout = tmp;
360  }
361 
362  double trklen = sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
363  // if entry and exit hit are the same (seems to happen rarely), trklen = 0
364  // which leads to a 0/0 and an NaN in edep later on
365  // this code does for particles in the same cell a trklen/trklen (vdedx[ii]/trklen)
366  // so setting this to any non zero number will do just fine
367  // I just pick -1 here to flag those strange hits in case I want t oanalyze them
368  // later on
369  if (trklen == 0)
370  {
371  trklen = -1.;
372  }
373  vector<int> vx;
374  vector<int> veta;
375  vector<double> vdedx;
376 
377  if (intxbin == intxbinout && intetabin == intetabinout) // single cell fired
378  {
379  if (Verbosity() > 0) cout << "SINGLE CELL FIRED: " << intxbin << " " << intetabin << endl;
380  vx.push_back(intxbin);
381  veta.push_back(intetabin);
382  vdedx.push_back(trklen);
383  }
384  else
385  {
386  for (int ibp = intxbin; ibp <= intxbinout; ibp++)
387  {
388  for (int ibz = intetabin; ibz <= intetabinout; ibz++)
389  {
390  double cx = geo->get_xcenter(ibp) - geo->get_xstep() / 2.;
391  double dx = geo->get_xcenter(ibp) + geo->get_xstep() / 2.;
392  double cy = geo->get_etacenter(ibz) - geo->get_etastep() / 2.;
393  double dy = geo->get_etacenter(ibz) + geo->get_etastep() / 2.;
394  //cout << "##### line: " << ax << " " << ay << " " << bx << " " << by << endl;
395  //cout << "####### cell: " << cx << " " << cy << " " << dx << " " << dy << endl;
396  pair<bool, double> intersect = PHG4Utils::line_and_rectangle_intersect(ax, ay, bx, by, cx, cy, dx, dy);
397  if (intersect.first)
398  {
399  if (Verbosity() > 0) cout << "CELL FIRED: " << ibp << " " << ibz << " " << intersect.second << endl;
400  vx.push_back(ibp);
401  veta.push_back(ibz);
402  vdedx.push_back(intersect.second);
403  }
404  }
405  }
406  }
407  if (Verbosity() > 0) cout << "NUMBER OF FIRED CELLS = " << vx.size() << endl;
408 
409  double tmpsum = 0.;
410  for (unsigned int ii = 0; ii < vx.size(); ii++)
411  {
412  tmpsum += vdedx[ii];
413  vdedx[ii] = vdedx[ii] / trklen;
414  if (Verbosity() > 0) cout << " CELL " << ii << " dE/dX = " << vdedx[ii] << endl;
415  }
416  if (Verbosity() > 0) cout << " TOTAL TRACK LENGTH = " << tmpsum << " " << trklen << endl;
417 
418  for (unsigned int i1 = 0; i1 < vx.size(); i1++) // loop over all fired cells
419  {
420  int ixbin = vx[i1];
421  int ietabin = veta[i1];
422  int ibin = ixbin * nzbins + ietabin;
423  if (!cellptarray[ibin])
424  {
425  PHG4CellDefs::keytype key = PHG4CellDefs::EtaXsizeBinning::genkey(*layer, ixbin, ietabin);
426  cellptarray[ibin] = new PHG4Cellv1(key);
427  }
428  cellptarray[ibin]->add_edep(hiter->first, hiter->second->get_edep() * vdedx[i1]);
429  cellptarray[ibin]->add_edep(hiter->second->get_edep() * vdedx[i1]);
430  if (hiter->second->has_property(PHG4Hit::prop_light_yield))
431  {
432  cellptarray[ibin]->add_light_yield(hiter->second->get_light_yield() * vdedx[i1]);
433  }
434  cellptarray[ibin]->add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep() * vdedx[i1]);
435 
436  // just a sanity check - we don't want to mess up by having Nan's or Infs in our energy deposition
437  if (!isfinite(hiter->second->get_edep() * vdedx[i1]))
438  {
439  cout << PHWHERE << " invalid energy dep " << hiter->second->get_edep()
440  << " or path length: " << vdedx[i1] << endl;
441  }
442  }
443 
444  vx.clear();
445  veta.clear();
446  } // end loop over g4hits
447 
448  int numcells = 0;
449  for (int ix = 0; ix < nxbins; ix++)
450  {
451  for (int iz = 0; iz < nzbins; iz++)
452  {
453  int ibin = ix * nzbins + iz;
454 
455  if (cellptarray[ibin])
456  {
457  cells->AddCell(cellptarray[ibin]);
458  numcells++;
459  if (Verbosity() > 1)
460  {
461  cout << "Adding cell in bin x: " << ix
462  << " x: " << geo->get_xcenter(ix) * 180. / M_PI
463  << ", eta bin: " << iz
464  << ", eta: " << geo->get_etacenter(iz)
465  << ", energy dep: " << cellptarray[ibin]->get_edep()
466  << endl;
467  }
468 
469  cellptarray[ibin] = nullptr;
470  }
471  }
472  }
473 
474  if (Verbosity() > 0)
475  {
476  cout << Name() << ": found " << numcells << " eta/x cells with energy deposition" << endl;
477  }
478  }
479  }
480 
481  if (chkenergyconservation)
482  {
483  CheckEnergy(topNode);
484  }
485 
486  return Fun4AllReturnCodes::EVENT_OK;
487 }
488 
489 void PHG4BlockCellReco::etaxsize(const int detid, const double deltaeta, const double deltax)
490 {
491  set_double_param(detid, "deltaeta", deltaeta);
492  set_double_param(detid, "deltax", deltax);
493  return;
494 }
495 
496 void PHG4BlockCellReco::set_timing_window(const int detid, const double tmin, const double tmax)
497 {
498  set_double_param(detid, "tmin", tmin);
499  set_double_param(detid, "tmax", tmax);
500  return;
501 }
502 
503 void PHG4BlockCellReco::set_size(const int i, const double sizeA, const double sizeB, const int what)
504 {
505  if (binning.find(i) != binning.end())
506  {
507  cout << "size for layer " << i << " already set" << endl;
508  return;
509  }
510 
511  binning[i] = what;
512  cell_size[i] = std::make_pair(sizeA, sizeB);
513 
514  return;
515 }
516 
517 //---------------------------------------------------------------
518 
519 int PHG4BlockCellReco::CheckEnergy(PHCompositeNode *topNode)
520 {
521  PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
522  double sum_energy_cells = 0.;
523  double sum_energy_stored_hits = 0.;
524  double sum_energy_stored_showers = 0.;
525 
526  PHG4CellContainer::ConstRange cell_begin_end = cells->getCells();
527  PHG4CellContainer::ConstIterator citer;
528  for (citer = cell_begin_end.first; citer != cell_begin_end.second; ++citer)
529  {
530  sum_energy_cells += citer->second->get_edep();
531  PHG4Cell::EdepConstRange cellrange = citer->second->get_g4hits();
532  for (PHG4Cell::EdepConstIterator iter = cellrange.first; iter != cellrange.second; ++iter)
533  {
534  sum_energy_stored_hits += iter->second;
535  }
536  PHG4Cell::ShowerEdepConstRange shwrrange = citer->second->get_g4showers();
537  for (PHG4Cell::ShowerEdepConstIterator iter = shwrrange.first; iter != shwrrange.second; ++iter)
538  {
539  sum_energy_stored_showers += iter->second;
540  }
541  }
542 
543  // the fractional eloss for particles traversing eta bins leads to minute rounding errors
544  if (sum_energy_stored_hits > 0)
545  {
546  if (fabs(sum_energy_cells - sum_energy_stored_hits) / sum_energy_cells > 1e-6)
547  {
548  cout << "energy mismatch between cell energy " << sum_energy_cells
549  << " and stored hit energies " << sum_energy_stored_hits
550  << endl;
551  }
552  }
553  if (sum_energy_stored_showers > 0)
554  {
555  if (fabs(sum_energy_cells - sum_energy_stored_showers) / sum_energy_cells > 1e-6)
556  {
557  cout << "energy mismatch between cell energy " << sum_energy_cells
558  << " and stored shower energies " << sum_energy_stored_showers
559  << endl;
560  }
561  }
562 
563  if (fabs(sum_energy_cells - sum_energy_g4hit) / sum_energy_g4hit > 1e-6)
564  {
565  cout << "energy mismatch between cells: " << sum_energy_cells
566  << " and hits: " << sum_energy_g4hit
567  << " diff sum(cells) - sum(hits): " << sum_energy_cells - sum_energy_g4hit
568  << endl;
569  return -1;
570  }
571 
572  else
573  {
574  if (Verbosity() > 0)
575  {
576  cout << Name() << ": sum hit energy: " << sum_energy_g4hit << " GeV" << endl;
577  cout << Name() << ": sum cell energy: " << sum_energy_cells << " GeV" << endl;
578  cout << Name() << ": sum shower energy: " << sum_energy_stored_showers << " GeV" << endl;
579  cout << Name() << ": sum stored hit energy: " << sum_energy_stored_hits << " GeV" << endl;
580  }
581  }
582 
583  return 0;
584 }
585 
586 void PHG4BlockCellReco::SetDefaultParameters()
587 {
588  set_default_double_param("deltaeta", NAN);
589  set_default_double_param("deltax", NAN);
590  set_default_double_param("tmax", 60.0);
591  set_default_double_param("tmin", 0.0);
592  return;
593 }