ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CaloAna.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CaloAna.cc
1 #include "CaloAna.h"
2 
3 // G4Hits includes
4 #include <g4main/PHG4Hit.h>
6 
7 // G4Cells includes
8 #include <g4detectors/PHG4Cell.h>
10 
11 // Tower includes
12 #include <calobase/RawTower.h>
13 #include <calobase/RawTowerContainer.h>
14 #include <calobase/RawTowerGeom.h>
15 #include <calobase/RawTowerGeomContainer.h>
16 
17 // Cluster includes
18 #include <calobase/RawCluster.h>
19 #include <calobase/RawClusterContainer.h>
20 
23 
24 #include <phool/getClass.h>
25 
26 #include <TFile.h>
27 #include <TNtuple.h>
28 
29 #include <cassert>
30 #include <sstream>
31 #include <string>
32 
33 using namespace std;
34 
35 CaloAna::CaloAna(const std::string& name, const std::string& filename)
36  : SubsysReco(name)
37  , detector("HCALIN")
38  , outfilename(filename)
39  , hm(nullptr)
40  , outfile(nullptr)
41  , g4hitntuple(nullptr)
42  , g4cellntuple(nullptr)
43  , towerntuple(nullptr)
44  , clusterntuple(nullptr)
45 {
46 }
47 
49 {
50  delete hm;
51  delete g4hitntuple;
52  delete g4cellntuple;
53  delete towerntuple;
54  delete clusterntuple;
55 }
56 
58 {
59  hm = new Fun4AllHistoManager(Name());
60  // create and register your histos (all types) here
61  // TH1 *h1 = new TH1F("h1",....)
62  // hm->registerHisto(h1);
63  outfile = new TFile(outfilename.c_str(), "RECREATE");
64  g4hitntuple = new TNtuple("hitntup", "G4Hits", "x0:y0:z0:x1:y1:z1:edep");
65  g4cellntuple = new TNtuple("cellntup", "G4Cells", "phi:eta:edep");
66  towerntuple = new TNtuple("towerntup", "Towers", "phi:eta:energy");
67  clusterntuple = new TNtuple("clusterntup", "Clusters", "phi:z:energy:towers");
68  return 0;
69 }
70 
72 {
73  // For the calorimeters we have the following node name convention
74  // where detector is the calorimeter name (CEMC, HCALIN, HCALOUT)
75  // this is the order in which they are reconstructed
76  // G4HIT_<detector>: G4 Hits
77  // G4CELL_<detector>: Cells (combined hits inside a cell - scintillator, eta/phi bin)
78  // TOWER_SIM_<detector>: simulated tower (before pedestal and noise)
79  // TOWER_RAW_<detector>: Raw Tower (adc/tdc values - from sims or real data)
80  // TOWER_CALIB_<detector>: Calibrated towers
81  // CLUSTER_<detector>: clusters
82  process_g4hits(topNode);
83  process_g4cells(topNode);
84  process_towers(topNode);
85  process_clusters(topNode);
87 }
88 
90 {
91  ostringstream nodename;
92 
93  // loop over the G4Hits
94  nodename.str("");
95  nodename << "G4HIT_" << detector;
96  PHG4HitContainer* hits = findNode::getClass<PHG4HitContainer>(topNode, nodename.str().c_str());
97  if (hits)
98  {
99  // this returns an iterator to the beginning and the end of our G4Hits
100  PHG4HitContainer::ConstRange hit_range = hits->getHits();
101  for (PHG4HitContainer::ConstIterator hit_iter = hit_range.first; hit_iter != hit_range.second; hit_iter++)
102 
103  {
104  // the pointer to the G4Hit is hit_iter->second
105  g4hitntuple->Fill(hit_iter->second->get_x(0),
106  hit_iter->second->get_y(0),
107  hit_iter->second->get_z(0),
108  hit_iter->second->get_x(1),
109  hit_iter->second->get_y(1),
110  hit_iter->second->get_z(1),
111  hit_iter->second->get_edep());
112  }
113  }
115 }
116 
118 {
119  ostringstream nodename;
120 
121  // loop over the G4Hits
122  nodename.str("");
123  nodename << "G4CELL_" << detector;
124  PHG4CellContainer* cells = findNode::getClass<PHG4CellContainer>(topNode, nodename.str());
125  if (cells)
126  {
127  PHG4CellContainer::ConstRange cell_range = cells->getCells();
128  int phibin = -999;
129  int etabin = -999;
130  for (PHG4CellContainer::ConstIterator cell_iter = cell_range.first; cell_iter != cell_range.second; cell_iter++)
131  {
132  if (cell_iter->second->has_binning(PHG4CellDefs::scintillatorslatbinning))
133  {
134  phibin = PHG4CellDefs::ScintillatorSlatBinning::get_row(cell_iter->second->get_cellid());
135  etabin = PHG4CellDefs::ScintillatorSlatBinning::get_column(cell_iter->second->get_cellid());
136  }
137  else if (cell_iter->second->has_binning(PHG4CellDefs::sizebinning))
138  {
139  phibin = PHG4CellDefs::SizeBinning::get_phibin(cell_iter->second->get_cellid());
140  etabin = PHG4CellDefs::SizeBinning::get_zbin(cell_iter->second->get_cellid());
141  }
142  else if (cell_iter->second->has_binning(PHG4CellDefs::spacalbinning))
143  {
144  phibin = PHG4CellDefs::SpacalBinning::get_phibin(cell_iter->second->get_cellid());
145  etabin = PHG4CellDefs::SpacalBinning::get_etabin(cell_iter->second->get_cellid());
146  }
147  else
148  {
149  cout << "unknown cell binning, implement 0x" << hex << PHG4CellDefs::get_binning(cell_iter->second->get_cellid()) << dec << endl;
150  }
151  g4cellntuple->Fill(
152  phibin,
153  etabin,
154  cell_iter->second->get_edep());
155  }
156  }
158 }
159 
161 {
162  ostringstream nodename;
163  ostringstream geonodename;
164 
165  // loop over the G4Hits
166  nodename.str("");
167  nodename << "TOWER_CALIB_" << detector;
168  geonodename.str("");
169  geonodename << "TOWERGEOM_" << detector;
170  RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, geonodename.str().c_str());
171  if (!towergeom)
172  {
174  }
175  RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(topNode, nodename.str().c_str());
176  if (towers)
177  {
178  // again pair of iterators to begin and end of tower map
179  RawTowerContainer::ConstRange tower_range = towers->getTowers();
180  for (RawTowerContainer::ConstIterator tower_iter = tower_range.first; tower_iter != tower_range.second; tower_iter++)
181 
182  {
183  int phibin = tower_iter->second->get_binphi();
184  int etabin = tower_iter->second->get_bineta();
185  double phi = towergeom->get_phicenter(phibin);
186  double eta = towergeom->get_etacenter(etabin);
187  towerntuple->Fill(phi,
188  eta,
189  tower_iter->second->get_energy());
190  }
191  }
193 }
194 
196 {
197  ostringstream nodename;
198 
199  // loop over the G4Hits
200  nodename.str("");
201  nodename << "CLUSTER_" << detector;
202  RawClusterContainer* clusters = findNode::getClass<RawClusterContainer>(topNode, nodename.str().c_str());
203  if (clusters)
204  {
205  RawClusterContainer::ConstRange cluster_range = clusters->getClusters();
206  for (RawClusterContainer::ConstIterator cluster_iter = cluster_range.first; cluster_iter != cluster_range.second; cluster_iter++)
207  {
208  clusterntuple->Fill(cluster_iter->second->get_phi(),
209  cluster_iter->second->get_z(),
210  cluster_iter->second->get_energy(),
211  cluster_iter->second->getNTowers());
212  }
213  }
215 }
216 
218 {
219  outfile->cd();
220  g4hitntuple->Write();
221  g4cellntuple->Write();
222  towerntuple->Write();
223  clusterntuple->Write();
224  outfile->Write();
225  outfile->Close();
226  delete outfile;
227  hm->dumpHistos(outfilename, "UPDATE");
228  return 0;
229 }