Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TPCIntegratedCharge.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TPCIntegratedCharge.cc
1 #include "TPCIntegratedCharge.h"
2 #include "TPCDaqDefs.h"
3 
4 #include <g4detectors/PHG4Cell.h>
5 #include <g4detectors/PHG4CellContainer.h>
6 #include <g4detectors/PHG4CylinderCellGeom.h>
7 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
8 #include <g4main/PHG4Hit.h>
9 #include <g4main/PHG4HitContainer.h>
10 
11 #include <fun4all/Fun4AllHistoManager.h>
12 #include <fun4all/Fun4AllReturnCodes.h>
13 #include <fun4all/Fun4AllServer.h>
14 #include <fun4all/PHTFileServer.h>
15 
16 #include <phhepmc/PHHepMCGenEvent.h>
17 #include <phhepmc/PHHepMCGenEventMap.h>
18 #include <phool/PHCompositeNode.h>
19 #include <phool/getClass.h>
20 
21 #include <TFile.h>
22 #include <TH1D.h>
23 #include <TH2D.h>
24 #include <TString.h>
25 #include <TTree.h>
26 #include <TVector3.h>
27 
28 #include <CLHEP/Units/SystemOfUnits.h>
29 
30 #include <boost/format.hpp>
31 
32 #include <algorithm>
33 #include <array>
34 #include <cassert>
35 #include <cmath>
36 #include <iostream>
37 #include <limits>
38 #include <stdexcept>
39 
40 using namespace std;
41 using namespace CLHEP;
42 
43 TPCIntegratedCharge::TPCIntegratedCharge(
44  unsigned int minLayer,
45  unsigned int m_maxLayer,
46  const std::string& outputfilename)
47  : SubsysReco("TPCIntegratedCharge")
48  , m_outputFileName(outputfilename)
49  , m_minLayer(minLayer)
50  , m_maxLayer(m_maxLayer)
51 {
52 }
53 
54 TPCIntegratedCharge::~TPCIntegratedCharge()
55 {
56 }
57 
58 int TPCIntegratedCharge::Init(PHCompositeNode* topNode)
59 {
60  return Fun4AllReturnCodes::EVENT_OK;
61 }
62 
63 int TPCIntegratedCharge::End(PHCompositeNode* topNode)
64 {
65  if (Verbosity() >= VERBOSITY_SOME)
66  cout << "TPCIntegratedCharge::End - write to " << m_outputFileName << endl;
67  PHTFileServer::get().cd(m_outputFileName);
68 
69  Fun4AllHistoManager* hm = getHistoManager();
70  assert(hm);
71  for (unsigned int i = 0; i < hm->nHistos(); i++)
72  hm->getHisto(i)->Write();
73 
74  // help index files with TChain
75  TTree* T_Index = new TTree("T_Index", "T_Index");
76  assert(T_Index);
77  T_Index->Write();
78 
79  return Fun4AllReturnCodes::EVENT_OK;
80 }
81 
82 int TPCIntegratedCharge::InitRun(PHCompositeNode* topNode)
83 {
84  // PHG4CylinderCellGeomContainer* seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
85  // if (!seggeo)
86  // {
87  // cout << "could not locate geo node " << "CYLINDERCELLGEOM_SVTX" << endl;
88  // exit(1);
89  // }
90 
91  if (Verbosity() >= VERBOSITY_SOME)
92  cout << "TPCIntegratedCharge::get_HistoManager - Making PHTFileServer " << m_outputFileName
93  << endl;
94  PHTFileServer::get().open(m_outputFileName, "RECREATE");
95 
96  Fun4AllHistoManager* hm = getHistoManager();
97  assert(hm);
98 
99  TH1D* h = new TH1D("hNormalization", //
100  "Normalization;Items;Summed quantity", 10, .5, 10.5);
101  int i = 1;
102  h->GetXaxis()->SetBinLabel(i++, "Event count");
103  h->GetXaxis()->SetBinLabel(i++, "Collision count");
104  h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit");
105  h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit Edep");
106  h->GetXaxis()->SetBinLabel(i++, "TPC Pad Hit");
107  h->GetXaxis()->SetBinLabel(i++, "TPC Charge e");
108  h->GetXaxis()->SetBinLabel(i++, "TPC Charge fC");
109  h->GetXaxis()->LabelsOption("v");
110  hm->registerHisto(h);
111 
112  // for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
113  // {
114  // const PHG4CylinderCellGeom* layer_geom = seggeo->GetLayerCellGeom(layer);
115 
116  // const string histNameCellHit(boost::str(boost::format{"hCellHit_Layer%1%"} % layer));
117  // const string histNameCellCharge(boost::str(boost::format{"hCellCharge_Layer%1%"} % layer));
118 
119  // }
120 
121  hm->registerHisto(new TH2D("hLayerCellHit", //
122  "Number of ADC time-bin hit per channel;Layer ID;Hit number",
123  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
124  300, -.5, 299.5));
125  hm->registerHisto(new TH2D("hLayerCellCharge", //
126  "Charge integrated over drift window per channel;Layer ID;Charge [fC]",
127  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
128  1000, 0, 1e7 * eplus / (1e-15 * coulomb)));
129 
130  hm->registerHisto(new TH2D("hLayerSumCellHit", //
131  "Number of ADC time-bin hit integrated over channels per layer;Layer ID;Hit number",
132  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
133  10000, -.5, 99999.5));
134  hm->registerHisto(new TH2D("hLayerSumCellCharge", //
135  "Charge integrated over drift window and channel per layer;Layer ID;Charge [fC]",
136  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
137  10000, 0, 1000 * 4e6 * eplus / (1e-15 * coulomb)));
138 
139  return Fun4AllReturnCodes::EVENT_OK;
140 }
141 
142 int TPCIntegratedCharge::process_event(PHCompositeNode* topNode)
143 {
144  Fun4AllHistoManager* hm = getHistoManager();
145  assert(hm);
146  TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto("hNormalization"));
147  assert(h_norm);
148 
149  PHG4HitContainer* g4hit = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_SVTX");
150  if (!g4hit)
151  {
152  cout << "TPCIntegratedCharge::process_event - Could not locate g4 hit node G4HIT_SVTX" << endl;
153  exit(1);
154  }
155  PHG4CellContainer* cells = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_SVTX");
156  if (!cells)
157  {
158  cout << "TPCIntegratedCharge::process_event - could not locate cell node "
159  << "G4CELL_SVTX" << endl;
160  exit(1);
161  }
162 
163  PHG4CylinderCellGeomContainer* seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
164  if (!seggeo)
165  {
166  cout << "TPCIntegratedCharge::process_event - could not locate geo node "
167  << "CYLINDERCELLGEOM_SVTX" << endl;
168  exit(1);
169  }
170  PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
171  if (!geneventmap)
172  {
173  static bool once = true;
174  if (once)
175  {
176  once = false;
177 
178  std::cout << "TPCIntegratedCharge::process_event - missing node PHHepMCGenEventMap. Skipping HepMC stat." << std::endl;
179  }
180  }
181  else
182  {
183  h_norm->Fill("Collision count", geneventmap->size());
184  }
185 
186  for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
187  {
188  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(layer);
189 
190  for (PHG4HitContainer::ConstIterator hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
191  {
192  const double edep = hiter->second->get_edep();
193  h_norm->Fill("TPC G4Hit Edep", edep);
194  h_norm->Fill("TPC G4Hit", 1);
195  } // for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
196 
197  } // for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
198 
199  // prepreare charge stat.
200  unsigned int nZBins = 0;
201  vector<array<vector<int>, 2> > layerChanCellHit(m_maxLayer + 1);
202  vector<array<vector<double>, 2> > layerChanCellCharge(m_maxLayer + 1);
203  for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
204  {
205  PHG4CylinderCellGeom* layerGeom =
206  seggeo->GetLayerCellGeom(layer);
207  assert(layerGeom);
208 
209  // start with an empty vector of cells for each phibin
210  const int nphibins = layerGeom->get_phibins();
211  assert(nphibins > 0);
212 
213  if (Verbosity() >= VERBOSITY_MORE)
214  {
215  cout << "TPCIntegratedCharge::process_event - init layer " << layer << " with "
216  << "nphibins = " << nphibins
217  << ", layerGeom->get_zbins() = " << layerGeom->get_zbins() << endl;
218  }
219 
220  if (nZBins <= 0)
221  {
222  nZBins = layerGeom->get_zbins();
223  assert(nZBins > 0);
224  }
225  else
226  {
227  if ((int) nZBins != layerGeom->get_zbins())
228  {
229  cout << "TPCIntegratedCharge::process_event - Fatal Error - nZBin at layer " << layer << " is " << layerGeom->get_zbins()
230  << ", which is different from previous layers of nZBin = " << nZBins << endl;
231  exit(1);
232  }
233  }
234 
235  layerChanCellHit[layer][0].resize(nphibins, 0);
236  layerChanCellHit[layer][1].resize(nphibins, 0);
237 
238  layerChanCellCharge[layer][0].resize(nphibins, 0);
239  layerChanCellCharge[layer][1].resize(nphibins, 0);
240  }
241  assert(nZBins > 0);
242 
243  // count all cells
244  PHG4CellContainer::ConstRange cellrange = cells->getCells();
245  for (PHG4CellContainer::ConstIterator celliter = cellrange.first;
246  celliter != cellrange.second;
247  ++celliter)
248  {
249  PHG4Cell* cell = celliter->second;
250  const unsigned int layer = cell->get_layer();
251  const unsigned int phibin = PHG4CellDefs::SizeBinning::get_phibin(cell->get_cellid());
252  const unsigned int zbin = PHG4CellDefs::SizeBinning::get_zbin(cell->get_cellid());
253 
254  if (Verbosity() >= VERBOSITY_MORE)
255  {
256  cout << "TPCIntegratedCharge::process_event - process cell: "
257  << "layer = " << layer
258  << ", get_phibin = " << phibin
259  << ", get_zbin = " << zbin
260  << ", get_edep = " << cell->get_edep()
261  << ", get_eion = " << cell->get_eion() << ": "
262  << "Class = " << cell->ClassName();
263  cell->identify();
264  }
265 
266  if (layer < m_minLayer or layer > m_maxLayer) continue;
267 
268  if (not(zbin >= 0))
269  {
270  cout << "TPCIntegratedCharge::process_event - invalid cell: "
271  << "layer = " << layer
272  << ", get_phibin = " << phibin
273  << ", get_zbin = " << zbin
274  << ", get_edep = " << cell->get_edep()
275  << ", get_eion = " << cell->get_eion() << ": "
276  << "Class = " << cell->ClassName();
277  cell->identify();
278  }
279 
280  assert(zbin >= 0);
281  assert(zbin < nZBins);
282  const int side = (zbin < nZBins / 2) ? 0 : 1;
283 
284  vector<int>& ChanCellHit = layerChanCellHit[layer][side];
285  vector<double>& ChanCellCharge = layerChanCellCharge[layer][side];
286 
287  assert(phibin >= 0);
288  assert(phibin < ChanCellHit.size());
289  assert(phibin < ChanCellHit.size());
290 
291  ChanCellHit[phibin] += 1;
292  ChanCellCharge[phibin] += cell->get_edep();
293  }
294 
295  // fill histograms
296  TH2* hLayerCellHit = dynamic_cast<TH2*>(hm->getHisto("hLayerCellHit"));
297  assert(hLayerCellHit);
298  TH2* hLayerCellCharge = dynamic_cast<TH2*>(hm->getHisto("hLayerCellCharge"));
299  assert(hLayerCellCharge);
300  TH2* hLayerSumCellHit = dynamic_cast<TH2*>(hm->getHisto("hLayerSumCellHit"));
301  assert(hLayerSumCellHit);
302  TH2* hLayerSumCellCharge = dynamic_cast<TH2*>(hm->getHisto("hLayerSumCellCharge"));
303  assert(hLayerSumCellCharge);
304 
305  for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
306  {
307  for (unsigned int side = 0; side < 2; ++side)
308  {
309  int sumHit = 0;
310  for (unsigned int phibin = 0; phibin < layerChanCellHit[layer][side].size(); ++phibin)
311  {
312  const int& hit = layerChanCellHit[layer][side][phibin];
313  sumHit += hit;
314 
315  hLayerCellHit->Fill(layer, hit);
316  h_norm->Fill("TPC Pad Hit", hit);
317 
318  if (Verbosity() >= VERBOSITY_MORE)
319  {
320  cout << "TPCIntegratedCharge::process_event - hit: "
321  << "hit = " << hit
322  << "at layer = " << layer
323  << ", side = " << side
324  << ", phibin = " << phibin
325  << endl;
326  }
327  }
328 
329  if (Verbosity() >= VERBOSITY_MORE)
330  {
331  cout << "TPCIntegratedCharge::process_event - hLayerSumCellHit->Fill(" << layer << ", " << sumHit << ")" << endl;
332  }
333 
334  hLayerSumCellHit->Fill(layer, sumHit);
335 
336  double sum_charge_fC = 0;
337  for (unsigned int phibin = 0; phibin < layerChanCellCharge[layer][side].size(); ++phibin)
338  {
339  const double& charge_e = layerChanCellCharge[layer][side][phibin];
340  const double charge_fC = charge_e * eplus / (1e-15 * coulomb);
341  sum_charge_fC += charge_fC;
342 
343  hLayerCellCharge->Fill(layer, charge_fC);
344 
345  h_norm->Fill("TPC Charge e", charge_e);
346  h_norm->Fill("TPC Charge fC", charge_fC);
347  }
348  hLayerSumCellCharge->Fill(layer, sum_charge_fC);
349  } // for (unsigned int side = 0; side < 2; ++side)
350 
351  } // for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
352 
353  h_norm->Fill("Event count", 1);
354 
355  return Fun4AllReturnCodes::EVENT_OK;
356 }
357 
358 Fun4AllHistoManager*
359 TPCIntegratedCharge::getHistoManager()
360 {
361  static string histname("TPCIntegratedCharge_HISTOS");
362 
363  Fun4AllServer* se = Fun4AllServer::instance();
364  Fun4AllHistoManager* hm = se->getHistoManager(histname);
365 
366  if (not hm)
367  {
368  cout
369  << "TPCIntegratedCharge::get_HistoManager - Making Fun4AllHistoManager " << histname
370  << endl;
371  hm = new Fun4AllHistoManager(histname);
372  se->registerHistoManager(hm);
373  }
374 
375  assert(hm);
376 
377  return hm;
378 }