ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4CellNtuple.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4CellNtuple.cc
1 #include "G4CellNtuple.h"
2 
3 #include <g4detectors/PHG4Cell.h>
5 #include <g4detectors/PHG4CellDefs.h> // for get_phibin
8 
10 #include <fun4all/SubsysReco.h> // for SubsysReco
11 
12 #include <phool/getClass.h>
13 
14 #include <TFile.h>
15 #include <TH1.h>
16 #include <TNtuple.h>
17 
18 #include <cassert>
19 #include <cmath> // for isfinite
20 #include <iostream> // for operator<<
21 #include <sstream>
22 #include <utility> // for pair
23 
24 using namespace std;
25 
26 G4CellNtuple::G4CellNtuple(const std::string &name, const std::string &filename)
27  : SubsysReco(name)
28  , nblocks(0)
29  , hm(nullptr)
30  , _filename(filename)
31  , ntup(nullptr)
32  , outfile(nullptr)
33 {
34 }
35 
37 {
38  // delete ntup;
39  delete hm;
40 }
41 
43 {
44  hm = new Fun4AllHistoManager(Name());
45  outfile = new TFile(_filename.c_str(), "RECREATE");
46  ntup = new TNtuple("cellntup", "G4Cells", "detid:layer:phi:eta:edep");
47  // ntup->SetDirectory(0);
48  TH1 *h1 = new TH1F("edep1GeV", "edep 0-1GeV", 1000, 0, 1);
49  eloss.push_back(h1);
50  h1 = new TH1F("edep100GeV", "edep 0-100GeV", 1000, 0, 100);
51  eloss.push_back(h1);
52  return 0;
53 }
54 
56 {
57  ostringstream nodename;
58  ostringstream geonodename;
59  set<string>::const_iterator iter;
60  vector<TH1 *>::const_iterator eiter;
61  for (iter = _node_postfix.begin(); iter != _node_postfix.end(); ++iter)
62  {
63  int detid = (_detid.find(*iter))->second;
64  nodename.str("");
65  nodename << "G4CELL_" << *iter;
66  geonodename.str("");
67  geonodename << "CYLINDERCELLGEOM_" << *iter;
68  PHG4CylinderCellGeomContainer *cellgeos = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, geonodename.str());
69  if (!cellgeos)
70  {
71  cout << "no geometry node " << geonodename.str() << " for " << *iter << endl;
72  continue;
73  }
74  PHG4CellContainer *cells = findNode::getClass<PHG4CellContainer>(topNode, nodename.str());
75  if (cells)
76  {
77  int previouslayer = -1;
78  PHG4CylinderCellGeom *cellgeom = nullptr;
79  double esum = 0;
80  // double numcells = cells->size();
81  // ncells[i]->Fill(numcells);
82  // cout << "number of cells: " << cells->size() << endl;
83  PHG4CellContainer::ConstRange cell_range = cells->getCells();
84  for (PHG4CellContainer::ConstIterator cell_iter = cell_range.first; cell_iter != cell_range.second; cell_iter++)
85 
86  {
87  double edep = cell_iter->second->get_edep();
88  if (!isfinite(edep))
89  {
90  cout << "invalid edep: " << edep << endl;
91  }
92  esum += cell_iter->second->get_edep();
93  int phibin = ~0x0;
94  int etabin = ~0x0;
95  double phi = NAN;
96  double eta = NAN;
97  int layer = cell_iter->second->get_layer();
98  // to search the map fewer times, cache the geom object until the layer changes
99  if (layer != previouslayer)
100  {
101  cellgeom = cellgeos->GetLayerCellGeom(layer);
102  previouslayer = layer;
103  }
104  if (cell_iter->second->has_binning(PHG4CellDefs::etaphibinning))
105  {
106  phibin = PHG4CellDefs::EtaPhiBinning::get_phibin(cell_iter->second->get_cellid());
107  etabin = PHG4CellDefs::EtaPhiBinning::get_etabin(cell_iter->second->get_cellid());
108  phi = cellgeom->get_phicenter(phibin);
109  eta = cellgeom->get_etacenter(etabin);
110  }
111  else if (cell_iter->second->has_binning(PHG4CellDefs::sizebinning))
112  {
113  phibin = PHG4CellDefs::SizeBinning::get_phibin(cell_iter->second->get_cellid());
114  etabin = PHG4CellDefs::SizeBinning::get_zbin(cell_iter->second->get_cellid());
115  phi = cellgeom->get_phicenter(phibin);
116  eta = cellgeom->get_zcenter(etabin);
117  }
118  assert(cellgeom != nullptr);
119  ntup->Fill(detid,
120  cell_iter->second->get_layer(),
121  phi,
122  eta,
123  cell_iter->second->get_edep());
124  }
125  for (eiter = eloss.begin(); eiter != eloss.end(); ++eiter)
126  {
127  (*eiter)->Fill(esum);
128  }
129  }
130  }
131  return 0;
132 }
133 
135 {
136  outfile->cd();
137  ntup->Write();
138  outfile->Write();
139  outfile->Close();
140  delete outfile;
141  hm->dumpHistos(_filename, "UPDATE");
142  return 0;
143 }
144 
145 void G4CellNtuple::AddNode(const std::string &name, const int detid)
146 {
147  _node_postfix.insert(name);
148  _detid[name] = detid;
149  return;
150 }