File indexing completed on 2025-08-06 08:16:11
0001 #include "TPCIntegratedCharge.h"
0002 #include "TPCDaqDefs.h"
0003
0004 #include <g4detectors/PHG4Cell.h>
0005 #include <g4detectors/PHG4CellContainer.h>
0006 #include <g4detectors/PHG4CylinderCellGeom.h>
0007 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0008 #include <g4main/PHG4Hit.h>
0009 #include <g4main/PHG4HitContainer.h>
0010
0011 #include <fun4all/Fun4AllHistoManager.h>
0012 #include <fun4all/Fun4AllReturnCodes.h>
0013 #include <fun4all/Fun4AllServer.h>
0014 #include <fun4all/PHTFileServer.h>
0015
0016 #include <phhepmc/PHHepMCGenEvent.h>
0017 #include <phhepmc/PHHepMCGenEventMap.h>
0018 #include <phool/PHCompositeNode.h>
0019 #include <phool/getClass.h>
0020
0021 #include <TFile.h>
0022 #include <TH1D.h>
0023 #include <TH2D.h>
0024 #include <TString.h>
0025 #include <TTree.h>
0026 #include <TVector3.h>
0027
0028 #include <CLHEP/Units/SystemOfUnits.h>
0029
0030 #include <boost/format.hpp>
0031
0032 #include <algorithm>
0033 #include <array>
0034 #include <cassert>
0035 #include <cmath>
0036 #include <iostream>
0037 #include <limits>
0038 #include <stdexcept>
0039
0040 using namespace std;
0041 using namespace CLHEP;
0042
0043 TPCIntegratedCharge::TPCIntegratedCharge(
0044 unsigned int minLayer,
0045 unsigned int m_maxLayer,
0046 const std::string& outputfilename)
0047 : SubsysReco("TPCIntegratedCharge")
0048 , m_outputFileName(outputfilename)
0049 , m_minLayer(minLayer)
0050 , m_maxLayer(m_maxLayer)
0051 {
0052 }
0053
0054 TPCIntegratedCharge::~TPCIntegratedCharge()
0055 {
0056 }
0057
0058 int TPCIntegratedCharge::Init(PHCompositeNode* topNode)
0059 {
0060 return Fun4AllReturnCodes::EVENT_OK;
0061 }
0062
0063 int TPCIntegratedCharge::End(PHCompositeNode* topNode)
0064 {
0065 if (Verbosity() >= VERBOSITY_SOME)
0066 cout << "TPCIntegratedCharge::End - write to " << m_outputFileName << endl;
0067 PHTFileServer::get().cd(m_outputFileName);
0068
0069 Fun4AllHistoManager* hm = getHistoManager();
0070 assert(hm);
0071 for (unsigned int i = 0; i < hm->nHistos(); i++)
0072 hm->getHisto(i)->Write();
0073
0074
0075 TTree* T_Index = new TTree("T_Index", "T_Index");
0076 assert(T_Index);
0077 T_Index->Write();
0078
0079 return Fun4AllReturnCodes::EVENT_OK;
0080 }
0081
0082 int TPCIntegratedCharge::InitRun(PHCompositeNode* topNode)
0083 {
0084
0085
0086
0087
0088
0089
0090
0091 if (Verbosity() >= VERBOSITY_SOME)
0092 cout << "TPCIntegratedCharge::get_HistoManager - Making PHTFileServer " << m_outputFileName
0093 << endl;
0094 PHTFileServer::get().open(m_outputFileName, "RECREATE");
0095
0096 Fun4AllHistoManager* hm = getHistoManager();
0097 assert(hm);
0098
0099 TH1D* h = new TH1D("hNormalization",
0100 "Normalization;Items;Summed quantity", 10, .5, 10.5);
0101 int i = 1;
0102 h->GetXaxis()->SetBinLabel(i++, "Event count");
0103 h->GetXaxis()->SetBinLabel(i++, "Collision count");
0104 h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit");
0105 h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit Edep");
0106 h->GetXaxis()->SetBinLabel(i++, "TPC Pad Hit");
0107 h->GetXaxis()->SetBinLabel(i++, "TPC Charge e");
0108 h->GetXaxis()->SetBinLabel(i++, "TPC Charge fC");
0109 h->GetXaxis()->LabelsOption("v");
0110 hm->registerHisto(h);
0111
0112
0113
0114
0115
0116
0117
0118
0119
0120
0121 hm->registerHisto(new TH2D("hLayerCellHit",
0122 "Number of ADC time-bin hit per channel;Layer ID;Hit number",
0123 m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0124 300, -.5, 299.5));
0125 hm->registerHisto(new TH2D("hLayerCellCharge",
0126 "Charge integrated over drift window per channel;Layer ID;Charge [fC]",
0127 m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0128 1000, 0, 1e7 * eplus / (1e-15 * coulomb)));
0129
0130 hm->registerHisto(new TH2D("hLayerSumCellHit",
0131 "Number of ADC time-bin hit integrated over channels per layer;Layer ID;Hit number",
0132 m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0133 10000, -.5, 99999.5));
0134 hm->registerHisto(new TH2D("hLayerSumCellCharge",
0135 "Charge integrated over drift window and channel per layer;Layer ID;Charge [fC]",
0136 m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0137 10000, 0, 1000 * 4e6 * eplus / (1e-15 * coulomb)));
0138
0139 return Fun4AllReturnCodes::EVENT_OK;
0140 }
0141
0142 int TPCIntegratedCharge::process_event(PHCompositeNode* topNode)
0143 {
0144 Fun4AllHistoManager* hm = getHistoManager();
0145 assert(hm);
0146 TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto("hNormalization"));
0147 assert(h_norm);
0148
0149 PHG4HitContainer* g4hit = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_SVTX");
0150 if (!g4hit)
0151 {
0152 cout << "TPCIntegratedCharge::process_event - Could not locate g4 hit node G4HIT_SVTX" << endl;
0153 exit(1);
0154 }
0155 PHG4CellContainer* cells = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_SVTX");
0156 if (!cells)
0157 {
0158 cout << "TPCIntegratedCharge::process_event - could not locate cell node "
0159 << "G4CELL_SVTX" << endl;
0160 exit(1);
0161 }
0162
0163 PHG4CylinderCellGeomContainer* seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0164 if (!seggeo)
0165 {
0166 cout << "TPCIntegratedCharge::process_event - could not locate geo node "
0167 << "CYLINDERCELLGEOM_SVTX" << endl;
0168 exit(1);
0169 }
0170 PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0171 if (!geneventmap)
0172 {
0173 static bool once = true;
0174 if (once)
0175 {
0176 once = false;
0177
0178 std::cout << "TPCIntegratedCharge::process_event - missing node PHHepMCGenEventMap. Skipping HepMC stat." << std::endl;
0179 }
0180 }
0181 else
0182 {
0183 h_norm->Fill("Collision count", geneventmap->size());
0184 }
0185
0186 for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0187 {
0188 PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(layer);
0189
0190 for (PHG4HitContainer::ConstIterator hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
0191 {
0192 const double edep = hiter->second->get_edep();
0193 h_norm->Fill("TPC G4Hit Edep", edep);
0194 h_norm->Fill("TPC G4Hit", 1);
0195 }
0196
0197 }
0198
0199
0200 unsigned int nZBins = 0;
0201 vector<array<vector<int>, 2> > layerChanCellHit(m_maxLayer + 1);
0202 vector<array<vector<double>, 2> > layerChanCellCharge(m_maxLayer + 1);
0203 for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0204 {
0205 PHG4CylinderCellGeom* layerGeom =
0206 seggeo->GetLayerCellGeom(layer);
0207 assert(layerGeom);
0208
0209
0210 const int nphibins = layerGeom->get_phibins();
0211 assert(nphibins > 0);
0212
0213 if (Verbosity() >= VERBOSITY_MORE)
0214 {
0215 cout << "TPCIntegratedCharge::process_event - init layer " << layer << " with "
0216 << "nphibins = " << nphibins
0217 << ", layerGeom->get_zbins() = " << layerGeom->get_zbins() << endl;
0218 }
0219
0220 if (nZBins <= 0)
0221 {
0222 nZBins = layerGeom->get_zbins();
0223 assert(nZBins > 0);
0224 }
0225 else
0226 {
0227 if ((int) nZBins != layerGeom->get_zbins())
0228 {
0229 cout << "TPCIntegratedCharge::process_event - Fatal Error - nZBin at layer " << layer << " is " << layerGeom->get_zbins()
0230 << ", which is different from previous layers of nZBin = " << nZBins << endl;
0231 exit(1);
0232 }
0233 }
0234
0235 layerChanCellHit[layer][0].resize(nphibins, 0);
0236 layerChanCellHit[layer][1].resize(nphibins, 0);
0237
0238 layerChanCellCharge[layer][0].resize(nphibins, 0);
0239 layerChanCellCharge[layer][1].resize(nphibins, 0);
0240 }
0241 assert(nZBins > 0);
0242
0243
0244 PHG4CellContainer::ConstRange cellrange = cells->getCells();
0245 for (PHG4CellContainer::ConstIterator celliter = cellrange.first;
0246 celliter != cellrange.second;
0247 ++celliter)
0248 {
0249 PHG4Cell* cell = celliter->second;
0250 const unsigned int layer = cell->get_layer();
0251 const unsigned int phibin = PHG4CellDefs::SizeBinning::get_phibin(cell->get_cellid());
0252 const unsigned int zbin = PHG4CellDefs::SizeBinning::get_zbin(cell->get_cellid());
0253
0254 if (Verbosity() >= VERBOSITY_MORE)
0255 {
0256 cout << "TPCIntegratedCharge::process_event - process cell: "
0257 << "layer = " << layer
0258 << ", get_phibin = " << phibin
0259 << ", get_zbin = " << zbin
0260 << ", get_edep = " << cell->get_edep()
0261 << ", get_eion = " << cell->get_eion() << ": "
0262 << "Class = " << cell->ClassName();
0263 cell->identify();
0264 }
0265
0266 if (layer < m_minLayer or layer > m_maxLayer) continue;
0267
0268 if (not(zbin >= 0))
0269 {
0270 cout << "TPCIntegratedCharge::process_event - invalid cell: "
0271 << "layer = " << layer
0272 << ", get_phibin = " << phibin
0273 << ", get_zbin = " << zbin
0274 << ", get_edep = " << cell->get_edep()
0275 << ", get_eion = " << cell->get_eion() << ": "
0276 << "Class = " << cell->ClassName();
0277 cell->identify();
0278 }
0279
0280 assert(zbin >= 0);
0281 assert(zbin < nZBins);
0282 const int side = (zbin < nZBins / 2) ? 0 : 1;
0283
0284 vector<int>& ChanCellHit = layerChanCellHit[layer][side];
0285 vector<double>& ChanCellCharge = layerChanCellCharge[layer][side];
0286
0287 assert(phibin >= 0);
0288 assert(phibin < ChanCellHit.size());
0289 assert(phibin < ChanCellHit.size());
0290
0291 ChanCellHit[phibin] += 1;
0292 ChanCellCharge[phibin] += cell->get_edep();
0293 }
0294
0295
0296 TH2* hLayerCellHit = dynamic_cast<TH2*>(hm->getHisto("hLayerCellHit"));
0297 assert(hLayerCellHit);
0298 TH2* hLayerCellCharge = dynamic_cast<TH2*>(hm->getHisto("hLayerCellCharge"));
0299 assert(hLayerCellCharge);
0300 TH2* hLayerSumCellHit = dynamic_cast<TH2*>(hm->getHisto("hLayerSumCellHit"));
0301 assert(hLayerSumCellHit);
0302 TH2* hLayerSumCellCharge = dynamic_cast<TH2*>(hm->getHisto("hLayerSumCellCharge"));
0303 assert(hLayerSumCellCharge);
0304
0305 for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0306 {
0307 for (unsigned int side = 0; side < 2; ++side)
0308 {
0309 int sumHit = 0;
0310 for (unsigned int phibin = 0; phibin < layerChanCellHit[layer][side].size(); ++phibin)
0311 {
0312 const int& hit = layerChanCellHit[layer][side][phibin];
0313 sumHit += hit;
0314
0315 hLayerCellHit->Fill(layer, hit);
0316 h_norm->Fill("TPC Pad Hit", hit);
0317
0318 if (Verbosity() >= VERBOSITY_MORE)
0319 {
0320 cout << "TPCIntegratedCharge::process_event - hit: "
0321 << "hit = " << hit
0322 << "at layer = " << layer
0323 << ", side = " << side
0324 << ", phibin = " << phibin
0325 << endl;
0326 }
0327 }
0328
0329 if (Verbosity() >= VERBOSITY_MORE)
0330 {
0331 cout << "TPCIntegratedCharge::process_event - hLayerSumCellHit->Fill(" << layer << ", " << sumHit << ")" << endl;
0332 }
0333
0334 hLayerSumCellHit->Fill(layer, sumHit);
0335
0336 double sum_charge_fC = 0;
0337 for (unsigned int phibin = 0; phibin < layerChanCellCharge[layer][side].size(); ++phibin)
0338 {
0339 const double& charge_e = layerChanCellCharge[layer][side][phibin];
0340 const double charge_fC = charge_e * eplus / (1e-15 * coulomb);
0341 sum_charge_fC += charge_fC;
0342
0343 hLayerCellCharge->Fill(layer, charge_fC);
0344
0345 h_norm->Fill("TPC Charge e", charge_e);
0346 h_norm->Fill("TPC Charge fC", charge_fC);
0347 }
0348 hLayerSumCellCharge->Fill(layer, sum_charge_fC);
0349 }
0350
0351 }
0352
0353 h_norm->Fill("Event count", 1);
0354
0355 return Fun4AllReturnCodes::EVENT_OK;
0356 }
0357
0358 Fun4AllHistoManager*
0359 TPCIntegratedCharge::getHistoManager()
0360 {
0361 static string histname("TPCIntegratedCharge_HISTOS");
0362
0363 Fun4AllServer* se = Fun4AllServer::instance();
0364 Fun4AllHistoManager* hm = se->getHistoManager(histname);
0365
0366 if (not hm)
0367 {
0368 cout
0369 << "TPCIntegratedCharge::get_HistoManager - Making Fun4AllHistoManager " << histname
0370 << endl;
0371 hm = new Fun4AllHistoManager(histname);
0372 se->registerHistoManager(hm);
0373 }
0374
0375 assert(hm);
0376
0377 return hm;
0378 }
