ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4CylinderCellReco.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4CylinderCellReco.cc
1 #include "PHG4CylinderCellReco.h"
2 #include "PHG4CylinderCellGeom.h"
4 #include "PHG4CylinderGeom.h"
6 
7 #include "PHG4Cell.h" // for PHG4Cell, PHG...
8 #include "PHG4CellContainer.h"
9 #include "PHG4CellDefs.h"
10 #include "PHG4Cellv1.h"
11 
12 #include <phparameter/PHParameterContainerInterface.h> // for PHParameterCo...
13 #include <phparameter/PHParametersContainer.h>
14 
15 #include <g4main/PHG4Hit.h>
17 #include <g4main/PHG4Utils.h>
18 
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 <cstring> // for memset
33 #include <iostream>
34 #include <iterator> // for reverse_iterator
35 #include <sstream>
36 #include <vector> // for vector
37 
38 using namespace std;
39 
41  : SubsysReco(name)
43  , chkenergyconservation(0)
44  , sum_energy_before_cuts(0.)
45  , sum_energy_g4hit(0.)
46 {
48  memset(nbins, 0, sizeof(nbins));
49 }
50 
52 {
54  sum_energy_g4hit = 0.;
56 }
57 
59 {
60  PHNodeIterator nodeiter(topNode);
61 
62  // Looking for the DST node
63  PHCompositeNode *dstNode;
64  dstNode = dynamic_cast<PHCompositeNode *>(nodeiter.findFirst("PHCompositeNode", "DST"));
65  if (!dstNode)
66  {
67  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
68  exit(1);
69  }
70  PHCompositeNode *runNode;
71  runNode = dynamic_cast<PHCompositeNode *>(nodeiter.findFirst("PHCompositeNode", "RUN"));
72  if (!runNode)
73  {
74  cout << Name() << "DST Node missing, doing nothing." << endl;
75  exit(1);
76  }
77  PHNodeIterator runiter(runNode);
78  PHCompositeNode *RunDetNode =
79  dynamic_cast<PHCompositeNode *>(runiter.findFirst("PHCompositeNode",
80  detector));
81  if (!RunDetNode)
82  {
83  RunDetNode = new PHCompositeNode(detector);
84  runNode->addNode(RunDetNode);
85  }
86 
87  hitnodename = "G4HIT_" + detector;
88  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
89  if (!g4hit)
90  {
91  cout << "Could not locate g4 hit node " << hitnodename << endl;
92  exit(1);
93  }
94  cellnodename = "G4CELL_" + outdetector;
95  PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
96  if (!cells)
97  {
98  PHNodeIterator dstiter(dstNode);
99  PHCompositeNode *DetNode =
100  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode",
101  detector));
102  if (!DetNode)
103  {
104  DetNode = new PHCompositeNode(detector);
105  dstNode->addNode(DetNode);
106  }
107  cells = new PHG4CellContainer();
108  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(cells, cellnodename.c_str(), "PHObject");
109  DetNode->addNode(newNode);
110  }
111 
112  geonodename = "CYLINDERGEOM_" + detector;
113  PHG4CylinderGeomContainer *geo = findNode::getClass<PHG4CylinderGeomContainer>(topNode, geonodename.c_str());
114  if (!geo)
115  {
116  cout << "Could not locate geometry node " << geonodename << endl;
117  exit(1);
118  }
119  if (Verbosity() > 0)
120  {
121  geo->identify();
122  }
123  seggeonodename = "CYLINDERCELLGEOM_" + outdetector;
124  PHG4CylinderCellGeomContainer *seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, seggeonodename.c_str());
125  if (!seggeo)
126  {
127  seggeo = new PHG4CylinderCellGeomContainer();
128  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(seggeo, seggeonodename.c_str(), "PHObject");
129  runNode->addNode(newNode);
130  }
131 
133 
135 
136  map<int, PHG4CylinderGeom *>::const_iterator miter;
137  pair<map<int, PHG4CylinderGeom *>::const_iterator, map<int, PHG4CylinderGeom *>::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  PHG4CylinderGeom *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  set_size(layer, get_double_param(layer, "size_long"), get_double_param(layer, "size_perp"));
151  tmin_max.insert(std::make_pair(layer, std::make_pair(get_double_param(layer, "tmin"), get_double_param(layer, "tmax"))));
152  double circumference = layergeom->get_radius() * 2 * M_PI;
153  double length_in_z = layergeom->get_zmax() - layergeom->get_zmin();
154  sizeiter = cell_size.find(layer);
155  if (sizeiter == cell_size.end())
156  {
157  cout << "no cell sizes for layer " << layer << endl;
158  exit(1);
159  }
160  // create geo object and fill with variables common to all binning methods
161  PHG4CylinderCellGeom *layerseggeo = new PHG4CylinderCellGeom();
162  layerseggeo->set_layer(layergeom->get_layer());
163  layerseggeo->set_radius(layergeom->get_radius());
164  layerseggeo->set_thickness(layergeom->get_thickness());
165  if (binning[layer] == PHG4CellDefs::etaphibinning)
166  {
167  // calculate eta at radius+ thickness (outer radius)
168  // length via eta coverage is calculated using the outer radius
169  double etamin = PHG4Utils::get_eta(layergeom->get_radius() + layergeom->get_thickness(), layergeom->get_zmin());
170  double etamax = PHG4Utils::get_eta(layergeom->get_radius() + layergeom->get_thickness(), layergeom->get_zmax());
171  zmin_max[layer] = make_pair(etamin, etamax);
172  double etastepsize = (sizeiter->second).first;
173  double d_etabins;
174 // if the eta cell size is larger than the eta range, make one bin
175  if (etastepsize > etamax - etamin)
176  {
177  d_etabins = 1;
178  }
179  else
180  {
181  // it is unlikely that the eta range is a multiple of the eta cell size
182  // then fract is 0, if not - add 1 bin which makes the
183  // cells a tiny bit smaller but makes them fit
184  double fract = modf((etamax - etamin) / etastepsize, &d_etabins);
185  if (fract != 0)
186  {
187  d_etabins++;
188  }
189  }
190  etastepsize = (etamax - etamin) / d_etabins;
191  (sizeiter->second).first = etastepsize;
192  int etabins = d_etabins;
193  double etalow = etamin;
194  double etahi = etalow + etastepsize;
195  for (int i = 0; i < etabins; i++)
196  {
197  if (etahi > (etamax + 1.e-6)) // etahi is a tiny bit larger due to numerical uncertainties
198  {
199  cout << "etahi: " << etahi << ", etamax: " << etamax << endl;
200  }
201  etahi += etastepsize;
202  }
203 
204  double phimin = -M_PI;
205  double phimax = M_PI;
206  double phistepsize = (sizeiter->second).second;
207  double d_phibins;
208  if (phistepsize >= phimax-phimin)
209  {
210  d_phibins = 1;
211  }
212  else
213  {
214  // it is unlikely that the phi range is a multiple of the phi cell size
215  // then fract is 0, if not - add 1 bin which makes the
216  // cells a tiny bit smaller but makes them fit
217  double fract = modf((phimax - phimin) / phistepsize, &d_phibins);
218  if (fract != 0)
219  {
220  d_phibins++;
221  }
222  }
223  phistepsize = (phimax - phimin) / d_phibins;
224  (sizeiter->second).second = phistepsize;
225  int phibins = d_phibins;
226  double philow = phimin;
227  double phihi = philow + phistepsize;
228  for (int i = 0; i < phibins; i++)
229  {
230  if (phihi > phimax)
231  {
232  cout << "phihi: " << phihi << ", phimax: " << phimax << endl;
233  }
234  phihi += phistepsize;
235  }
236  pair<int, int> phi_z_bin = make_pair(phibins, etabins);
237  n_phi_z_bins[layer] = phi_z_bin;
239  layerseggeo->set_etabins(etabins);
240  layerseggeo->set_etamin(etamin);
241  layerseggeo->set_etastep(etastepsize);
242  layerseggeo->set_phimin(phimin);
243  layerseggeo->set_phibins(phibins);
244  layerseggeo->set_phistep(phistepsize);
245  phistep[layer] = phistepsize;
246  etastep[layer] = etastepsize;
247  }
248  else if (binning[layer] == PHG4CellDefs::sizebinning)
249  {
250  zmin_max[layer] = make_pair(layergeom->get_zmin(), layergeom->get_zmax());
251  double size_z = (sizeiter->second).second;
252  double size_r = (sizeiter->second).first;
253  double bins_r;
254  // if the size is larger than circumference, make it one bin
255  if (size_r >= circumference)
256  {
257  bins_r = 1;
258  }
259  else
260  {
261  // unlikely but if the circumference is a multiple of the cell size
262  // use result of division, if not - add 1 bin which makes the
263  // cells a tiny bit smaller but makes them fit
264  double fract = modf(circumference / size_r, &bins_r);
265  if (fract != 0)
266  {
267  bins_r++;
268  }
269  }
270  nbins[0] = bins_r;
271  size_r = circumference / bins_r;
272  (sizeiter->second).first = size_r;
273  double phistepsize = 2 * M_PI / bins_r;
274  double phimin = -M_PI;
275  double phimax = phimin + phistepsize;
276  phistep[layer] = phistepsize;
277  for (int i = 0; i < nbins[0]; i++)
278  {
279  if (phimax > (M_PI + 1e-9))
280  {
281  cout << "phimax: " << phimax << ", M_PI: " << M_PI
282  << "phimax-M_PI: " << phimax - M_PI << endl;
283  }
284  phimax += phistepsize;
285  }
286  // if the size is larger than length, make it one bin
287  if (size_z >= length_in_z)
288  {
289  bins_r = 1;
290  }
291  else
292  {
293  // unlikely but if the length is a multiple of the cell size
294  // use result of division, if not - add 1 bin which makes the
295  // cells a tiny bit smaller but makes them fit
296  double fract = modf(length_in_z / size_z, &bins_r);
297  if (fract != 0)
298  {
299  bins_r++;
300  }
301  }
302  nbins[1] = bins_r;
303  pair<int, int> phi_z_bin = make_pair(nbins[0], nbins[1]);
304  n_phi_z_bins[layer] = phi_z_bin;
305  // update our map with the new sizes
306  size_z = length_in_z / bins_r;
307  (sizeiter->second).second = size_z;
308  double zlow = layergeom->get_zmin();
309  double zhigh = zlow + size_z;
310  ;
311  for (int i = 0; i < nbins[1]; i++)
312  {
313  if (zhigh > (layergeom->get_zmax() + 1e-9))
314  {
315  cout << "zhigh: " << zhigh << ", zmax "
316  << layergeom->get_zmax()
317  << ", zhigh-zmax: " << zhigh - layergeom->get_zmax()
318  << endl;
319  }
320  zhigh += size_z;
321  }
323  layerseggeo->set_zbins(nbins[1]);
324  layerseggeo->set_zmin(layergeom->get_zmin());
325  layerseggeo->set_zstep(size_z);
326  layerseggeo->set_phibins(nbins[0]);
327  layerseggeo->set_phistep(phistepsize);
328  }
329  // add geo object filled by different binning methods
330  seggeo->AddLayerCellGeom(layerseggeo);
331  if (Verbosity() > 1)
332  {
333  layerseggeo->identify();
334  }
335  }
336 
337  // print out settings
338  if (Verbosity() > 0)
339  {
340  cout << "===================== PHG4CylinderCellReco::InitRun() =====================" << endl;
341  cout << " " << outdetector << " Segmentation Description: " << endl;
342  for (std::map<int, int>::iterator iter = binning.begin();
343  iter != binning.end(); ++iter)
344  {
345  int layer = iter->first;
346 
347  if (binning[layer] == PHG4CellDefs::etaphibinning)
348  {
349  // phi & eta bin is usually used to make projective towers
350  // so just print the first layer
351  cout << " Layer #" << binning.begin()->first << "-" << binning.rbegin()->first << endl;
352  cout << " Nbins (phi,eta): (" << n_phi_z_bins[layer].first << ", " << n_phi_z_bins[layer].second << ")" << endl;
353  cout << " Cell Size (phi,eta): (" << cell_size[layer].first << " rad, " << cell_size[layer].second << " units)" << endl;
354  break;
355  }
356  else if (binning[layer] == PHG4CellDefs::sizebinning)
357  {
358  cout << " Layer #" << layer << endl;
359  cout << " Nbins (phi,z): (" << n_phi_z_bins[layer].first << ", " << n_phi_z_bins[layer].second << ")" << endl;
360  cout << " Cell Size (phi,z): (" << cell_size[layer].first << " cm, " << cell_size[layer].second << " cm)" << endl;
361  }
362  }
363  cout << "===========================================================================" << endl;
364  }
365  string nodename = "G4CELLPARAM_" + GetParamsContainer()->Name();
366  SaveToNodeTree(RunDetNode, nodename);
368 }
369 
371 {
372  PHG4HitContainer *g4hit = findNode::getClass<PHG4HitContainer>(topNode, hitnodename.c_str());
373  if (!g4hit)
374  {
375  cout << "Could not locate g4 hit node " << hitnodename << endl;
376  exit(1);
377  }
378  PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
379  if (!cells)
380  {
381  cout << "could not locate cell node " << cellnodename << endl;
382  exit(1);
383  }
384 
385  PHG4CylinderCellGeomContainer *seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, seggeonodename.c_str());
386  if (!seggeo)
387  {
388  cout << "could not locate geo node " << seggeonodename << endl;
389  exit(1);
390  }
391 
392  map<int, std::pair<double, double> >::iterator sizeiter;
394  pair<PHG4HitContainer::LayerIter, PHG4HitContainer::LayerIter> layer_begin_end = g4hit->getLayers();
395  // cout << "number of layers: " << g4hit->num_layers() << endl;
396  // cout << "number of hits: " << g4hit->size() << endl;
397  // for (layer = layer_begin_end.first; layer != layer_begin_end.second; layer++)
398  // {
399  // cout << "layer number: " << *layer << endl;
400  // }
401  for (layer = layer_begin_end.first; layer != layer_begin_end.second; layer++)
402  {
403  // only handle layers/detector ids which have parameters set
404  if (implemented_detid.find(*layer) == implemented_detid.end())
405  {
406  continue;
407  }
409  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(*layer);
410  PHG4CylinderCellGeom *geo = seggeo->GetLayerCellGeom(*layer);
411  int nphibins = n_phi_z_bins[*layer].first;
412  int nzbins = n_phi_z_bins[*layer].second;
413 
414  // ------- eta/phi binning ------------------------------------------------------------------------
415  if (binning[*layer] == PHG4CellDefs::etaphibinning)
416  {
417  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
418  {
419  sum_energy_before_cuts += hiter->second->get_edep();
420  // checking ADC timing integration window cut
421  if (hiter->second->get_t(0) > tmin_max[*layer].second) continue;
422  if (hiter->second->get_t(1) < tmin_max[*layer].first) continue;
423 
424  pair<double, double> etaphi[2];
425  double phibin[2];
426  double etabin[2];
427  for (int i = 0; i < 2; i++)
428  {
429  etaphi[i] = PHG4Utils::get_etaphi(hiter->second->get_x(i), hiter->second->get_y(i), hiter->second->get_z(i));
430  etabin[i] = geo->get_etabin(etaphi[i].first);
431  phibin[i] = geo->get_phibin(etaphi[i].second);
432  }
433  // check bin range
434  if (phibin[0] < 0 || phibin[0] >= nphibins || phibin[1] < 0 || phibin[1] >= nphibins)
435  {
436  continue;
437  }
438  if (etabin[0] < 0 || etabin[0] >= nzbins || etabin[1] < 0 || etabin[1] >= nzbins)
439  {
440  continue;
441  }
442 
443  if (etabin[0] < 0)
444  {
445  if (Verbosity() > 0)
446  {
447  hiter->second->identify();
448  }
449  continue;
450  }
451  sum_energy_g4hit += hiter->second->get_edep();
452  int intphibin = phibin[0];
453  int intetabin = etabin[0];
454  int intphibinout = phibin[1];
455  int intetabinout = etabin[1];
456 
457  // Determine all fired cells
458 
459  double ax = (etaphi[0]).second; // phi
460  double ay = (etaphi[0]).first; // eta
461  double bx = (etaphi[1]).second;
462  double by = (etaphi[1]).first;
463  if (intphibin > intphibinout)
464  {
465  int tmp = intphibin;
466  intphibin = intphibinout;
467  intphibinout = tmp;
468  }
469  if (intetabin > intetabinout)
470  {
471  int tmp = intetabin;
472  intetabin = intetabinout;
473  intetabinout = tmp;
474  }
475 
476  double trklen = sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
477  // if entry and exit hit are the same (seems to happen rarely), trklen = 0
478  // which leads to a 0/0 and an NaN in edep later on
479  // this code does for particles in the same cell a trklen/trklen (vdedx[ii]/trklen)
480  // so setting this to any non zero number will do just fine
481  // I just pick -1 here to flag those strange hits in case I want t oanalyze them
482  // later on
483  if (trklen == 0)
484  {
485  trklen = -1.;
486  }
487  vector<int> vphi;
488  vector<int> veta;
489  vector<double> vdedx;
490 
491  if (intphibin == intphibinout && intetabin == intetabinout) // single cell fired
492  {
493  if (Verbosity() > 0) cout << "SINGLE CELL FIRED: " << intphibin << " " << intetabin << endl;
494  vphi.push_back(intphibin);
495  veta.push_back(intetabin);
496  vdedx.push_back(trklen);
497  }
498  else
499  {
500  double phistep_half = geo->get_phistep() / 2.;
501  double etastep_half = geo->get_etastep() / 2.;
502  for (int ibp = intphibin; ibp <= intphibinout; ibp++)
503  {
504  double cx = geo->get_phicenter(ibp) - phistep_half;
505  double dx = geo->get_phicenter(ibp) + phistep_half;
506  for (int ibz = intetabin; ibz <= intetabinout; ibz++)
507  {
508  double cy = geo->get_etacenter(ibz) - etastep_half;
509  double dy = geo->get_etacenter(ibz) + etastep_half;
510 
511  //cout << "##### line: " << ax << " " << ay << " " << bx << " " << by << endl;
512  //cout << "####### cell: " << cx << " " << cy << " " << dx << " " << dy << endl;
513  pair<bool, double> intersect = PHG4Utils::line_and_rectangle_intersect(ax, ay, bx, by, cx, cy, dx, dy);
514  if (intersect.first)
515  {
516  if (Verbosity() > 0) cout << "CELL FIRED: " << ibp << " " << ibz << " " << intersect.second << endl;
517  vphi.push_back(ibp);
518  veta.push_back(ibz);
519  vdedx.push_back(intersect.second);
520  }
521  }
522  }
523  }
524  if (Verbosity() > 0) cout << "NUMBER OF FIRED CELLS = " << vphi.size() << endl;
525 
526  double tmpsum = 0.;
527  for (unsigned int ii = 0; ii < vphi.size(); ii++)
528  {
529  tmpsum += vdedx[ii];
530  vdedx[ii] = vdedx[ii] / trklen;
531  if (Verbosity() > 0) cout << " CELL " << ii << " dE/dX = " << vdedx[ii] << endl;
532  }
533  if (Verbosity() > 0) cout << " TOTAL TRACK LENGTH = " << tmpsum << " " << trklen << endl;
534 
535  for (unsigned int i1 = 0; i1 < vphi.size(); i1++) // loop over all fired cells
536  {
537  int iphibin = vphi[i1];
538  int ietabin = veta[i1];
539 
540  // This is the key for cellptmap
541  // It is constructed using the phi and z (or eta) bin index values
542  // It will be unique for a given phi and z (or eta) bin combination
543  unsigned long long tmp = iphibin;
544  unsigned long long key = tmp << 32;
545  key += ietabin;
546  if (Verbosity() > 1)
547  {
548  cout << " iphibin " << iphibin << " ietabin " << ietabin << " key 0x" << hex << key << dec << endl;
549  }
550  PHG4Cell *cell = nullptr;
551  it = cellptmap.find(key);
552  if (it != cellptmap.end())
553  {
554  cell = it->second;
555  }
556  else
557  {
558  PHG4CellDefs::keytype cellkey = PHG4CellDefs::EtaPhiBinning::genkey(*layer, ietabin, iphibin);
559  cell = new PHG4Cellv1(cellkey);
560  cellptmap[key] = cell;
561  }
562  if (!isfinite(hiter->second->get_edep() * vdedx[i1]))
563  {
564  cout << "hit 0x" << hex << hiter->first << dec << " not finite, edep: "
565  << hiter->second->get_edep() << " weight " << vdedx[i1] << endl;
566  }
567  cell->add_edep(hiter->first, hiter->second->get_edep() * vdedx[i1]); // add hit with edep to g4hit list
568  cell->add_edep(hiter->second->get_edep() * vdedx[i1]); // add edep to cell
569  if (hiter->second->has_property(PHG4Hit::prop_light_yield))
570  {
571  cell->add_light_yield(hiter->second->get_light_yield() * vdedx[i1]);
572  }
573  cell->add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep() * vdedx[i1]);
574  // just a sanity check - we don't want to mess up by having Nan's or Infs in our energy deposition
575  if (!isfinite(hiter->second->get_edep() * vdedx[i1]))
576  {
577  cout << PHWHERE << " invalid energy dep " << hiter->second->get_edep()
578  << " or path length: " << vdedx[i1] << endl;
579  }
580  }
581  vphi.clear();
582  veta.clear();
583 
584  } // end loop over g4hits
585 
586  int numcells = 0;
587 
588  for (it = cellptmap.begin(); it != cellptmap.end(); ++it)
589  {
590  cells->AddCell(it->second);
591  numcells++;
592  if (Verbosity() > 1)
593  {
594  cout << "Adding cell in bin phi: " << PHG4CellDefs::EtaPhiBinning::get_phibin(it->second->get_cellid())
595  << " phi: " << geo->get_phicenter(PHG4CellDefs::EtaPhiBinning::get_phibin(it->second->get_cellid())) * 180. / M_PI
596  << ", eta bin: " << PHG4CellDefs::EtaPhiBinning::get_etabin(it->second->get_cellid())
597  << ", eta: " << geo->get_etacenter(PHG4CellDefs::EtaPhiBinning::get_etabin(it->second->get_cellid()))
598  << ", energy dep: " << it->second->get_edep()
599  << endl;
600  }
601  }
602 
603  if (Verbosity() > 0)
604  {
605  cout << Name() << ": found " << numcells << " eta/phi cells with energy deposition" << endl;
606  }
607  }
608 
609  else // ------ size binning ---------------------------------------------------------------
610  {
611  sizeiter = cell_size.find(*layer);
612  if (sizeiter == cell_size.end())
613  {
614  cout << "logical screwup!!! no sizes for layer " << *layer << endl;
615  exit(1);
616  }
617  double zstepsize = (sizeiter->second).second;
618  double phistepsize = phistep[*layer];
619 
620  for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
621  {
622  sum_energy_before_cuts += hiter->second->get_edep();
623  // checking ADC timing integration window cut
624  if (hiter->second->get_t(0) > tmin_max[*layer].second) continue;
625  if (hiter->second->get_t(1) < tmin_max[*layer].first) continue;
626 
627  double xinout[2];
628  double yinout[2];
629  double px[2];
630  double py[2];
631  double phi[2];
632  double z[2];
633  double phibin[2];
634  double zbin[2];
635  if (Verbosity() > 0) cout << "--------- new hit in layer # " << *layer << endl;
636 
637  for (int i = 0; i < 2; i++)
638  {
639  xinout[i] = hiter->second->get_x(i);
640  yinout[i] = hiter->second->get_y(i);
641  px[i] = hiter->second->get_px(i);
642  py[i] = hiter->second->get_py(i);
643  phi[i] = atan2(hiter->second->get_y(i), hiter->second->get_x(i));
644  z[i] = hiter->second->get_z(i);
645  phibin[i] = geo->get_phibin(phi[i]);
646  zbin[i] = geo->get_zbin(hiter->second->get_z(i));
647 
648  if (Verbosity() > 0) cout << " " << i << " phibin: " << phibin[i] << ", phi: " << phi[i] << ", stepsize: " << phistepsize << endl;
649  if (Verbosity() > 0) cout << " " << i << " zbin: " << zbin[i] << ", z = " << hiter->second->get_z(i) << ", stepsize: " << zstepsize << " offset: " << zmin_max[*layer].first << endl;
650  }
651  // check bin range
652  if (phibin[0] < 0 || phibin[0] >= nphibins || phibin[1] < 0 || phibin[1] >= nphibins)
653  {
654  continue;
655  }
656  if (zbin[0] < 0 || zbin[0] >= nzbins || zbin[1] < 0 || zbin[1] >= nzbins)
657  {
658  continue;
659  }
660 
661  if (zbin[0] < 0)
662  {
663  hiter->second->identify();
664  continue;
665  }
666  sum_energy_g4hit += hiter->second->get_edep();
667 
668  int intphibin = phibin[0];
669  int intzbin = zbin[0];
670  int intphibinout = phibin[1];
671  int intzbinout = zbin[1];
672 
673  if (Verbosity() > 0)
674  {
675  cout << " phi bin range: " << intphibin << " to " << intphibinout << " phi: " << phi[0] << " to " << phi[1] << endl;
676  cout << " Z bin range: " << intzbin << " to " << intzbinout << " Z: " << z[0] << " to " << z[1] << endl;
677  cout << " phi difference: " << (phi[1] - phi[0]) * 1000. << " milliradians." << endl;
678  cout << " phi difference: " << 2.5 * (phi[1] - phi[0]) * 10000. << " microns." << endl;
679  cout << " path length = " << sqrt((xinout[1] - xinout[0]) * (xinout[1] - xinout[0]) + (yinout[1] - yinout[0]) * (yinout[1] - yinout[0])) << endl;
680  cout << " px = " << px[0] << " " << px[1] << endl;
681  cout << " py = " << py[0] << " " << py[1] << endl;
682  cout << " x = " << xinout[0] << " " << xinout[1] << endl;
683  cout << " y = " << yinout[0] << " " << yinout[1] << endl;
684  }
685 
686  // Determine all fired cells
687 
688  double ax = phi[0];
689  double ay = z[0];
690  double bx = phi[1];
691  double by = z[1];
692  if (intphibin > intphibinout)
693  {
694  int tmp = intphibin;
695  intphibin = intphibinout;
696  intphibinout = tmp;
697  }
698  if (intzbin > intzbinout)
699  {
700  int tmp = intzbin;
701  intzbin = intzbinout;
702  intzbinout = tmp;
703  }
704 
705  double trklen = sqrt((ax - bx) * (ax - bx) + (ay - by) * (ay - by));
706  // if entry and exit hit are the same (seems to happen rarely), trklen = 0
707  // which leads to a 0/0 and an NaN in edep later on
708  // this code does for particles in the same cell a trklen/trklen (vdedx[ii]/trklen)
709  // so setting this to any non zero number will do just fine
710  // I just pick -1 here to flag those strange hits in case I want to analyze them
711  // later on
712  if (trklen == 0)
713  {
714  trklen = -1.;
715  }
716  vector<int> vphi;
717  vector<int> vz;
718  vector<double> vdedx;
719 
720  if (intphibin == intphibinout && intzbin == intzbinout) // single cell fired
721  {
722  if (Verbosity() > 0)
723  {
724  cout << "SINGLE CELL FIRED: " << intphibin << " " << intzbin << endl;
725  }
726  vphi.push_back(intphibin);
727  vz.push_back(intzbin);
728  vdedx.push_back(trklen);
729  }
730  else
731  {
732  double phistep_half = geo->get_phistep() / 2.;
733  double zstep_half = geo->get_zstep() / 2.;
734  for (int ibp = intphibin; ibp <= intphibinout; ibp++)
735  {
736  double cx = geo->get_phicenter(ibp) - phistep_half;
737  double dx = geo->get_phicenter(ibp) + phistep_half;
738  for (int ibz = intzbin; ibz <= intzbinout; ibz++)
739  {
740  double cy = geo->get_zcenter(ibz) - zstep_half;
741  double dy = geo->get_zcenter(ibz) + zstep_half;
742  //cout << "##### line: " << ax << " " << ay << " " << bx << " " << by << endl;
743  //cout << "####### cell: " << cx << " " << cy << " " << dx << " " << dy << endl;
744  pair<bool, double> intersect = PHG4Utils::line_and_rectangle_intersect(ax, ay, bx, by, cx, cy, dx, dy);
745  if (intersect.first)
746  {
747  if (Verbosity() > 0) cout << "CELL FIRED: " << ibp << " " << ibz << " " << intersect.second << endl;
748  vphi.push_back(ibp);
749  vz.push_back(ibz);
750  vdedx.push_back(intersect.second);
751  }
752  }
753  }
754  }
755  if (Verbosity() > 0) cout << "NUMBER OF FIRED CELLS = " << vz.size() << endl;
756 
757  double tmpsum = 0.;
758  for (unsigned int ii = 0; ii < vz.size(); ii++)
759  {
760  tmpsum += vdedx[ii];
761  vdedx[ii] = vdedx[ii] / trklen;
762  if (Verbosity() > 0) cout << " CELL " << ii << " dE/dX = " << vdedx[ii] << endl;
763  }
764  if (Verbosity() > 0) cout << " TOTAL TRACK LENGTH = " << tmpsum << " " << trklen << endl;
765 
766  for (unsigned int i1 = 0; i1 < vphi.size(); i1++) // loop over all fired cells
767  {
768  int iphibin = vphi[i1];
769  int izbin = vz[i1];
770 
771  unsigned long long tmp = iphibin;
772  unsigned long long key = tmp << 32;
773  key += izbin;
774  if (Verbosity() > 1)
775  {
776  cout << " iphibin " << iphibin << " izbin " << izbin << " key 0x" << hex << key << dec << endl;
777  }
778  // check to see if there is already an entry for this cell
779  PHG4Cell *cell = nullptr;
780  it = cellptmap.find(key);
781 
782  if (it != cellptmap.end())
783  {
784  cell = it->second;
785  if (Verbosity() > 1)
786  {
787  cout << " add energy to existing cell for key = " << cellptmap.find(key)->first << endl;
788  }
789 
790  if (Verbosity() > 1 && hiter->second->has_property(PHG4Hit::prop_light_yield) && std::isnan(hiter->second->get_light_yield() * vdedx[i1]))
791  {
792  cout << " NAN lighy yield with vdedx[i1] = " << vdedx[i1]
793  << " and hiter->second->get_light_yield() = " << hiter->second->get_light_yield() << endl;
794  }
795  }
796  else
797  {
798  if (Verbosity() > 1)
799  {
800  cout << " did not find a previous entry for key = 0x" << hex << key << dec << " create a new one" << endl;
801  }
802  PHG4CellDefs::keytype cellkey = PHG4CellDefs::SizeBinning::genkey(*layer, izbin, iphibin);
803  cell = new PHG4Cellv1(cellkey);
804  cellptmap[key] = cell;
805  }
806  if (!isfinite(hiter->second->get_edep() * vdedx[i1]))
807  {
808  cout << "hit 0x" << hex << hiter->first << dec << " not finite, edep: "
809  << hiter->second->get_edep() << " weight " << vdedx[i1] << endl;
810  }
811  cell->add_edep(hiter->first, hiter->second->get_edep() * vdedx[i1]);
812  cell->add_edep(hiter->second->get_edep() * vdedx[i1]); // add edep to cell
813  if (hiter->second->has_property(PHG4Hit::prop_light_yield))
814  {
815  cell->add_light_yield(hiter->second->get_light_yield() * vdedx[i1]);
816  if (Verbosity() > 1 && !std::isfinite(hiter->second->get_light_yield() * vdedx[i1]))
817  {
818  cout << " NAN lighy yield with vdedx[i1] = " << vdedx[i1]
819  << " and hiter->second->get_light_yield() = " << hiter->second->get_light_yield() << endl;
820  }
821  }
822  cell->add_shower_edep(hiter->second->get_shower_id(), hiter->second->get_edep() * vdedx[i1]);
823  }
824  vphi.clear();
825  vz.clear();
826 
827  } // end loop over hits
828 
829  int numcells = 0;
830 
831  for (it = cellptmap.begin(); it != cellptmap.end(); ++it)
832  {
833  cells->AddCell(it->second);
834  numcells++;
835  if (Verbosity() > 1)
836  {
837  cout << "Adding cell for key " << hex << it->first << dec << " in bin phi: " << PHG4CellDefs::SizeBinning::get_phibin(it->second->get_cellid())
838  << " phi: " << geo->get_phicenter(PHG4CellDefs::SizeBinning::get_phibin(it->second->get_cellid())) * 180. / M_PI
839  << ", z bin: " << PHG4CellDefs::SizeBinning::get_zbin(it->second->get_cellid())
840  << ", z: " << geo->get_zcenter(PHG4CellDefs::SizeBinning::get_zbin(it->second->get_cellid()))
841  << ", energy dep: " << it->second->get_edep()
842  << endl;
843  }
844  }
845 
846  if (Verbosity() > 0)
847  {
848  cout << "found " << numcells << " z/phi cells with energy deposition" << endl;
849  }
850  }
851 
852  //==========================================================
853  // now reset the cell map before moving on to the next layer
854  if (Verbosity() > 1)
855  {
856  cout << "cellptmap for layer " << *layer << " has final length " << cellptmap.size();
857  }
858  while (cellptmap.begin() != cellptmap.end())
859  {
860  // Assumes that memmory is freed by the cylinder cell container when it is destroyed
861  cellptmap.erase(cellptmap.begin());
862  }
863  if (Verbosity() > 1)
864  {
865  cout << " reset it to " << cellptmap.size() << endl;
866  }
867  }
869  {
870  CheckEnergy(topNode);
871  }
872 
874 }
875 
876 void PHG4CylinderCellReco::cellsize(const int detid, const double sr, const double sz)
877 {
878  if (binning.find(detid) != binning.end())
879  {
880  cout << "size for layer " << detid << " already set" << endl;
881  return;
882  }
884  set_double_param(detid, "size_long", sz);
885  set_double_param(detid, "size_perp", sr);
886 }
887 
888 void PHG4CylinderCellReco::etaphisize(const int detid, const double deltaeta, const double deltaphi)
889 {
890  if (binning.find(detid) != binning.end())
891  {
892  cout << "size for layer " << detid << " already set" << endl;
893  return;
894  }
896  set_double_param(detid, "size_long", deltaeta);
897  set_double_param(detid, "size_perp", deltaphi);
898  return;
899 }
900 
901 void PHG4CylinderCellReco::set_size(const int i, const double sizeA, const double sizeB)
902 {
903  cell_size[i] = std::make_pair(sizeA, sizeB);
904  return;
905 }
906 
907 void PHG4CylinderCellReco::set_timing_window(const int detid, const double tmin, const double tmax)
908 {
909  set_double_param(detid, "tmin", tmin);
910  set_double_param(detid, "tmax", tmax);
911  return;
912 }
913 
915 {
916  PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, cellnodename);
917  double sum_energy_cells = 0.;
918  double sum_energy_stored_hits = 0.;
919  double sum_energy_stored_showers = 0.;
920  PHG4CellContainer::ConstRange cell_begin_end = cells->getCells();
922  for (citer = cell_begin_end.first; citer != cell_begin_end.second; ++citer)
923  {
924  sum_energy_cells += citer->second->get_edep();
925  PHG4Cell::EdepConstRange cellrange = citer->second->get_g4hits();
926  for (PHG4Cell::EdepConstIterator iter = cellrange.first; iter != cellrange.second; ++iter)
927  {
928  sum_energy_stored_hits += iter->second;
929  }
930  PHG4Cell::ShowerEdepConstRange shwrrange = citer->second->get_g4showers();
931  for (PHG4Cell::ShowerEdepConstIterator iter = shwrrange.first; iter != shwrrange.second; ++iter)
932  {
933  sum_energy_stored_showers += iter->second;
934  }
935  }
936  // the fractional eloss for particles traversing eta bins leads to minute rounding errors
937  if (sum_energy_stored_hits > 0)
938  {
939  if (fabs(sum_energy_cells - sum_energy_stored_hits) / sum_energy_cells > 1e-6)
940  {
941  cout << "energy mismatch between cell energy " << sum_energy_cells
942  << " and stored hit energies " << sum_energy_stored_hits
943  << endl;
944  }
945  }
946  if (sum_energy_stored_showers > 0)
947  {
948  if (fabs(sum_energy_cells - sum_energy_stored_showers) / sum_energy_cells > 1e-6)
949  {
950  cout << "energy mismatch between cell energy " << sum_energy_cells
951  << " and stored shower energies " << sum_energy_stored_showers
952  << endl;
953  }
954  }
955  if (fabs(sum_energy_cells - sum_energy_g4hit) / sum_energy_g4hit > 1e-6)
956  {
957  cout << "energy mismatch between cells: " << sum_energy_cells
958  << " and hits: " << sum_energy_g4hit
959  << " diff sum(cells) - sum(hits): " << sum_energy_cells - sum_energy_g4hit
960  << " cut val " << fabs(sum_energy_cells - sum_energy_g4hit) / sum_energy_g4hit
961  << endl;
962  cout << Name() << ":total energy for this event: " << sum_energy_g4hit << " GeV" << endl;
963  cout << Name() << ": sum cell energy: " << sum_energy_cells << " GeV" << endl;
964  cout << Name() << ": sum shower energy: " << sum_energy_stored_showers << " GeV" << endl;
965  cout << Name() << ": sum stored hit energy: " << sum_energy_stored_hits << " GeV" << endl;
966  cout << Name() << ": hit energy before cuts: " << sum_energy_before_cuts << " GeV" << endl;
967  return -1;
968  }
969  else
970  {
971  if (Verbosity() > 0)
972  {
973  cout << Name() << ":total energy for this event: " << sum_energy_g4hit << " GeV" << endl;
974  cout << Name() << ": sum cell energy: " << sum_energy_cells << " GeV" << endl;
975  cout << Name() << ": sum shower energy: " << sum_energy_stored_showers << " GeV" << endl;
976  cout << Name() << ": sum stored hit energy: " << sum_energy_stored_hits << " GeV" << endl;
977  cout << Name() << ": hit energy before cuts: " << sum_energy_before_cuts << " GeV" << endl;
978  }
979  }
980  return 0;
981 }
982 
983 void PHG4CylinderCellReco::Detector(const std::string &d)
984 {
985  detector = d;
986  // only set the outdetector nodename if it wasn't set already
987  // in case the order in the macro is outdetector(); detector();
988  if (outdetector.size() == 0)
989  {
990  OutputDetector(d);
991  }
992  return;
993 }
994 
996 {
997  set_default_double_param("size_long", NAN);
998  set_default_double_param("size_perp", NAN);
999  set_default_double_param("tmax", 60.0);
1000  set_default_double_param("tmin", -20.0); // collision has a timing spread around the triggered event. Accepting negative time too.
1001  return;
1002 }