File indexing completed on 2025-08-06 08:16:10
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011 #include "TPCDataStreamEmulator.h"
0012
0013 #include "TPCDaqDefs.h"
0014
0015 #include <g4detectors/PHG4Cell.h>
0016 #include <g4detectors/PHG4CellContainer.h>
0017 #include <g4detectors/PHG4CylinderCellGeom.h>
0018 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0019 #include <trackbase_historic/SvtxHit.h>
0020 #include <trackbase_historic/SvtxHitMap.h>
0021 #include <g4main/PHG4Hit.h>
0022 #include <g4main/PHG4HitContainer.h>
0023 #include <g4main/PHG4Particle.h>
0024 #include <g4main/PHG4TruthInfoContainer.h>
0025
0026 #include <fun4all/Fun4AllHistoManager.h>
0027 #include <fun4all/Fun4AllReturnCodes.h>
0028 #include <fun4all/Fun4AllServer.h>
0029 #include <fun4all/PHTFileServer.h>
0030
0031 #include <phhepmc/PHHepMCGenEvent.h>
0032 #include <phhepmc/PHHepMCGenEventMap.h>
0033 #include <phool/PHCompositeNode.h>
0034 #include <phool/getClass.h>
0035
0036 #include <TDatabasePDG.h>
0037 #include <TFile.h>
0038 #include <TH1D.h>
0039 #include <TH2D.h>
0040 #include <TString.h>
0041 #include <TTree.h>
0042 #include <TVector3.h>
0043
0044 #include <HepMC/GenEvent.h>
0045
0046 #include <CLHEP/Units/SystemOfUnits.h>
0047
0048 #include <boost/format.hpp>
0049
0050 #include <algorithm>
0051 #include <array>
0052 #include <cassert>
0053 #include <cmath>
0054 #include <iostream>
0055 #include <limits>
0056 #include <stdexcept>
0057
0058 using namespace std;
0059 using namespace CLHEP;
0060
0061 TPCDataStreamEmulator::TPCDataStreamEmulator(
0062 unsigned int minLayer,
0063 unsigned int m_maxLayer,
0064 const std::string& outputfilename)
0065 : SubsysReco("TPCDataStreamEmulator")
0066 , m_saveDataStreamFile(true)
0067 , m_outputFileNameBase(outputfilename)
0068 , m_minLayer(minLayer)
0069 , m_maxLayer(m_maxLayer)
0070 , m_evtCounter(-1)
0071 , m_vertexZAcceptanceCut(10)
0072 , m_etaAcceptanceCut(1.1)
0073 , m_hDataSize(nullptr)
0074 , m_hWavelet(nullptr)
0075 , m_hNChEta(nullptr)
0076 , m_hLayerWaveletSize(nullptr)
0077 , m_hLayerHit(nullptr)
0078 , m_hLayerZBinHit(nullptr)
0079 , m_hLayerZBinADC(nullptr)
0080 , m_hLayerDataSize(nullptr)
0081 , m_hLayerSumHit(nullptr)
0082 , m_hLayerSumDataSize(nullptr)
0083 {
0084 }
0085
0086 TPCDataStreamEmulator::~TPCDataStreamEmulator()
0087 {
0088 }
0089
0090 int TPCDataStreamEmulator::Init(PHCompositeNode* topNode)
0091 {
0092 return Fun4AllReturnCodes::EVENT_OK;
0093 }
0094
0095 int TPCDataStreamEmulator::End(PHCompositeNode* topNode)
0096 {
0097 if (Verbosity() >= VERBOSITY_SOME)
0098 cout << "TPCDataStreamEmulator::End - write to " << m_outputFileNameBase + ".root" << endl;
0099 PHTFileServer::get().cd(m_outputFileNameBase + ".root");
0100
0101 Fun4AllHistoManager* hm = getHistoManager();
0102 assert(hm);
0103 for (unsigned int i = 0; i < hm->nHistos(); i++)
0104 hm->getHisto(i)->Write();
0105
0106
0107 TTree* T_Index = new TTree("T_Index", "T_Index");
0108 assert(T_Index);
0109 T_Index->Write();
0110
0111 return Fun4AllReturnCodes::EVENT_OK;
0112 }
0113
0114 int TPCDataStreamEmulator::InitRun(PHCompositeNode* topNode)
0115 {
0116 PHG4CylinderCellGeomContainer* seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0117 if (!seggeo)
0118 {
0119 cout << "could not locate geo node "
0120 << "CYLINDERCELLGEOM_SVTX" << endl;
0121 exit(1);
0122 }
0123
0124 int nZBins = 0;
0125 for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0126 {
0127 PHG4CylinderCellGeom* layerGeom =
0128 seggeo->GetLayerCellGeom(layer);
0129 assert(layerGeom);
0130
0131 if (nZBins <= 0)
0132 {
0133 nZBins = layerGeom->get_zbins();
0134 assert(nZBins > 0);
0135 }
0136 else
0137 {
0138 if ((int) nZBins != layerGeom->get_zbins())
0139 {
0140 cout << "TPCDataStreamEmulator::InitRun - Fatal Error - nZBin at layer " << layer << " is " << layerGeom->get_zbins()
0141 << ", which is different from previous layers of nZBin = " << nZBins << endl;
0142 exit(1);
0143 }
0144 }
0145 }
0146
0147 if (Verbosity() >= VERBOSITY_SOME)
0148 cout << "TPCDataStreamEmulator::get_HistoManager - Making PHTFileServer " << m_outputFileNameBase + ".root"
0149 << endl;
0150 PHTFileServer::get().open(m_outputFileNameBase + ".root", "RECREATE");
0151
0152 Fun4AllHistoManager* hm = getHistoManager();
0153 assert(hm);
0154
0155 TH1D* h = new TH1D("hNormalization",
0156 "Normalization;Items;Summed quantity", 10, .5, 10.5);
0157 int i = 1;
0158 h->GetXaxis()->SetBinLabel(i++, "Event count");
0159 h->GetXaxis()->SetBinLabel(i++, "Collision count");
0160 h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit");
0161 h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit Edep");
0162 h->GetXaxis()->SetBinLabel(i++, "TPC Hit");
0163 h->GetXaxis()->SetBinLabel(i++, "TPC Wavelet");
0164 h->GetXaxis()->SetBinLabel(i++, "TPC DataSize");
0165
0166 h->GetXaxis()->LabelsOption("v");
0167 hm->registerHisto(h);
0168
0169
0170
0171
0172
0173
0174
0175
0176
0177
0178 hm->registerHisto(m_hDataSize =
0179 new TH1D("hDataSize",
0180 "TPC Data Size per Event;Data size [Byte];Count",
0181 10000, 0, 20e6));
0182
0183 hm->registerHisto(m_hWavelet =
0184 new TH1D("hWavelet",
0185 "TPC Recorded Wavelet per Event;Wavelet count;Count",
0186 10000, 0, 4e6));
0187
0188 hm->registerHisto(m_hNChEta =
0189 new TH1D("hNChEta",
0190 "Charged particle #eta distribution;#eta;Count",
0191 1000, -5, 5));
0192
0193 hm->registerHisto(m_hLayerWaveletSize =
0194 new TH2D("hLayerWaveletSize",
0195 "Number of Recorded ADC sample per Wavelet;Layer ID;ADC Sample Count per Wavelet",
0196 m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0197 nZBins, -.5, nZBins - .5));
0198
0199 hm->registerHisto(m_hLayerHit =
0200 new TH2D("hLayerHit",
0201 "Number of Recorded ADC sample per channel;Layer ID;ADC Sample Count",
0202 m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0203 nZBins, -.5, nZBins - .5));
0204
0205 hm->registerHisto(m_hLayerDataSize =
0206 new TH2D("hLayerDataSize",
0207 "Data size per channel;Layer ID;Data size [Byte]",
0208 m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0209 2 * nZBins, 0, 2 * nZBins));
0210
0211 hm->registerHisto(m_hLayerSumHit =
0212 new TH2D("hLayerSumHit",
0213 "Number of Recorded ADC sample per layer;Layer ID;ADC Sample Count",
0214 m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0215 10000, -.5, 99999.5));
0216
0217 hm->registerHisto(m_hLayerSumDataSize =
0218 new TH2D("hLayerSumDataSize",
0219 "Data size per trigger per layer;Layer ID;Data size [Byte]",
0220 m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0221 1000, 0, .5e6));
0222
0223 hm->registerHisto(m_hLayerZBinHit =
0224 new TH2D("hLayerZBinHit",
0225 "Number of Recorded ADC sample per Time Bin;z bin ID;Layer ID",
0226 nZBins, -.5, nZBins - .5,
0227 m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5));
0228
0229 hm->registerHisto(m_hLayerZBinADC =
0230 new TH2D("hLayerZBinADC",
0231 "Sum ADC per Time Bin;z bin ID;Layer ID",
0232 nZBins, -.5, nZBins - .5,
0233 m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5));
0234 return Fun4AllReturnCodes::EVENT_OK;
0235 }
0236
0237 int TPCDataStreamEmulator::process_event(PHCompositeNode* topNode)
0238 {
0239 m_evtCounter += 1;
0240
0241 Fun4AllHistoManager* hm = getHistoManager();
0242 assert(hm);
0243 TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto("hNormalization"));
0244 assert(h_norm);
0245 h_norm->Fill("Event count", 1);
0246
0247 PHG4HitContainer* g4hit = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_SVTX");
0248 if (!g4hit)
0249 {
0250 cout << "TPCDataStreamEmulator::process_event - Could not locate g4 hit node G4HIT_SVTX" << endl;
0251 return Fun4AllReturnCodes::ABORTRUN;
0252 }
0253
0254 SvtxHitMap* hits = findNode::getClass<SvtxHitMap>(topNode, "SvtxHitMap");
0255 if (!hits)
0256 {
0257 cout << "PCDataStreamEmulator::process_event - ERROR: Can't find node SvtxHitMap" << endl;
0258 return Fun4AllReturnCodes::ABORTRUN;
0259 }
0260
0261 PHG4CellContainer* cells = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_SVTX");
0262 if (!cells)
0263 {
0264 cout << "TPCDataStreamEmulator::process_event - could not locate cell node "
0265 << "G4CELL_SVTX" << endl;
0266 exit(1);
0267 }
0268
0269 PHG4CylinderCellGeomContainer* seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0270 if (!seggeo)
0271 {
0272 cout << "TPCDataStreamEmulator::process_event - could not locate geo node "
0273 << "CYLINDERCELLGEOM_SVTX" << endl;
0274 exit(1);
0275 }
0276
0277 PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0278 if (!geneventmap)
0279 {
0280 static bool once = true;
0281 if (once)
0282 {
0283 once = false;
0284
0285 cout << "TPCDataStreamEmulator::process_event - - missing node PHHepMCGenEventMap. Skipping HepMC stat." << std::endl;
0286 }
0287 }
0288 else
0289 {
0290 h_norm->Fill("Collision count", geneventmap->size());
0291 }
0292
0293 PHG4TruthInfoContainer* truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode,
0294 "G4TruthInfo");
0295 if (!truthInfoList)
0296 {
0297 cout << "TPCDataStreamEmulator::process_event - Fatal Error - "
0298 << "unable to find DST node "
0299 << "G4TruthInfo" << endl;
0300 exit(1);
0301 }
0302
0303 PHG4TruthInfoContainer::ConstRange primary_range =
0304 truthInfoList->GetPrimaryParticleRange();
0305
0306 for (PHG4TruthInfoContainer::ConstIterator particle_iter =
0307 primary_range.first;
0308 particle_iter != primary_range.second;
0309 ++particle_iter)
0310 {
0311 const PHG4Particle* p = particle_iter->second;
0312 assert(p);
0313
0314 TParticlePDG* pdg_p = TDatabasePDG::Instance()->GetParticle(
0315 p->get_pid());
0316 assert(pdg_p);
0317
0318 if (fabs(pdg_p->Charge()) > 0)
0319 {
0320 TVector3 pvec(p->get_px(), p->get_py(), p->get_pz());
0321
0322 if (pvec.Perp2()>0)
0323 {
0324 assert(m_hNChEta);
0325 m_hNChEta->Fill(pvec.PseudoRapidity());
0326 }
0327 }
0328
0329 }
0330
0331 for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0332 {
0333 PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(layer);
0334
0335 for (PHG4HitContainer::ConstIterator hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
0336 {
0337 const double edep = hiter->second->get_edep();
0338 h_norm->Fill("TPC G4Hit Edep", edep);
0339 h_norm->Fill("TPC G4Hit", 1);
0340 }
0341
0342 }
0343
0344
0345 int nZBins = 0;
0346 vector<array<vector<int>, 2> > layerChanHit(m_maxLayer + 1);
0347 vector<array<vector<int>, 2> > layerChanDataSize(m_maxLayer + 1);
0348 int nWavelet = 0;
0349 int sumDataSize = 0;
0350 for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0351 {
0352 PHG4CylinderCellGeom* layerGeom =
0353 seggeo->GetLayerCellGeom(layer);
0354 assert(layerGeom);
0355
0356
0357 const int nphibins = layerGeom->get_phibins();
0358 assert(nphibins > 0);
0359
0360 if (Verbosity() >= VERBOSITY_MORE)
0361 {
0362 cout << "TPCDataStreamEmulator::process_event - init layer " << layer << " with "
0363 << "nphibins = " << nphibins
0364 << ", layerGeom->get_zbins() = " << layerGeom->get_zbins() << endl;
0365 }
0366
0367 if (nZBins <= 0)
0368 {
0369 nZBins = layerGeom->get_zbins();
0370 assert(nZBins > 0);
0371 }
0372 else
0373 {
0374 if ((int) nZBins != layerGeom->get_zbins())
0375 {
0376 cout << "TPCDataStreamEmulator::process_event - Fatal Error - nZBin at layer " << layer << " is " << layerGeom->get_zbins()
0377 << ", which is different from previous layers of nZBin = " << nZBins << endl;
0378 exit(1);
0379 }
0380 }
0381
0382 for (unsigned int side = 0; side < 2; ++side)
0383 {
0384 layerChanHit[layer][side].resize(nphibins, 0);
0385
0386 layerChanDataSize[layer][side].resize(nphibins, 0);
0387 }
0388
0389 }
0390
0391 assert(nZBins > 0);
0392
0393
0394 int last_layer = -1;
0395 int last_side = -1;
0396 int last_phibin = -1;
0397 int last_zbin = -1;
0398 vector<unsigned int> last_wavelet;
0399 int last_wavelet_hittime = -1;
0400
0401 for (SvtxHitMap::Iter iter = hits->begin(); iter != hits->end(); ++iter)
0402 {
0403 SvtxHit* hit = iter->second;
0404
0405 const int layer = hit->get_layer();
0406
0407 if (layer < m_minLayer or layer > m_maxLayer) continue;
0408
0409 PHG4Cell* cell = cells->findCell(hit->get_cellid());
0410 const int phibin = PHG4CellDefs::SizeBinning::get_phibin(cell->get_cellid());
0411 const int zbin = PHG4CellDefs::SizeBinning::get_zbin(cell->get_cellid());
0412 const int side = (zbin < nZBins / 2) ? 0 : 1;
0413
0414
0415 if (last_layer != layer or last_phibin != phibin or last_side != side or abs(last_zbin - zbin) != 1)
0416 {
0417
0418 if (last_wavelet.size() > 0)
0419 {
0420 const int datasize = writeWavelet(last_layer, last_side, last_phibin, last_wavelet_hittime, last_wavelet);
0421 assert(datasize > 0);
0422
0423 nWavelet += 1;
0424 sumDataSize += datasize;
0425 layerChanDataSize[last_layer][last_side][last_phibin] += datasize;
0426
0427 last_wavelet.clear();
0428 last_zbin = -1;
0429 }
0430
0431
0432 PHG4CylinderCellGeom* layerGeom =
0433 seggeo->GetLayerCellGeom(layer);
0434 assert(layerGeom);
0435 const double z_abs = fabs(layerGeom->get_zcenter(zbin));
0436 const double r = layerGeom->get_radius();
0437 TVector3 acceptanceVec(r, 0, z_abs - m_vertexZAcceptanceCut);
0438 const double eta = acceptanceVec.PseudoRapidity();
0439
0440 if (eta > m_etaAcceptanceCut) continue;
0441
0442
0443 last_layer = layer;
0444 last_side = side;
0445 last_phibin = phibin;
0446
0447
0448 last_wavelet_hittime = (side == 0) ? (zbin) : (nZBins - 1 - zbin);
0449 assert(last_wavelet_hittime >= 0);
0450 assert(last_wavelet_hittime <= nZBins / 2);
0451 }
0452
0453 if (Verbosity() >= VERBOSITY_A_LOT)
0454 {
0455 cout << "TPCDataStreamEmulator::process_event - layer " << layer << " hit with "
0456
0457 << "phibin = " << phibin
0458 << ",zbin = " << zbin
0459 << ",side = " << side
0460 << ",last_wavelet.size() = " << last_wavelet.size()
0461 << ",last_zbin = " << last_zbin
0462 << endl;
0463 }
0464
0465
0466 if (last_wavelet.size() > 0)
0467 {
0468 if (side == 0)
0469 {
0470 assert(zbin - last_zbin == 1);
0471 }
0472 else
0473 {
0474 assert(last_zbin - zbin == 1);
0475 }
0476 }
0477
0478
0479 unsigned int adc = hit->get_adc();
0480 last_wavelet.push_back(adc);
0481 last_zbin = zbin;
0482
0483
0484 layerChanHit[layer][side][phibin] += 1;
0485 assert(m_hLayerZBinHit);
0486 m_hLayerZBinHit->Fill(zbin, layer, 1);
0487 assert(m_hLayerZBinADC);
0488 m_hLayerZBinADC->Fill(zbin, layer, adc);
0489
0490 }
0491
0492
0493 if (last_wavelet.size() > 0)
0494 {
0495 const int datasize = writeWavelet(last_layer, last_side, last_phibin, last_wavelet_hittime, last_wavelet);
0496 assert(datasize > 0);
0497
0498 nWavelet += 1;
0499 sumDataSize += datasize;
0500 layerChanDataSize[last_layer][last_side][last_phibin] += datasize;
0501 }
0502
0503
0504 for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0505 {
0506 for (unsigned int side = 0; side < 2; ++side)
0507 {
0508 int sumHit = 0;
0509 for (const int& hit : layerChanHit[layer][side])
0510 {
0511 sumHit += hit;
0512
0513 assert(m_hLayerHit);
0514 m_hLayerHit->Fill(layer, hit);
0515 h_norm->Fill("TPC Hit", hit);
0516
0517 if (Verbosity() >= VERBOSITY_MORE)
0518 {
0519 cout << "TPCDataStreamEmulator::process_event - hit: "
0520 << "hit = " << hit
0521 << "at layer = " << layer
0522 << ", side = " << side
0523 << endl;
0524 }
0525 }
0526
0527 if (Verbosity() >= VERBOSITY_MORE)
0528 {
0529 cout << "TPCDataStreamEmulator::process_event - hLayerSumCellHit->Fill(" << layer << ", " << sumHit << ")" << endl;
0530 }
0531 assert(m_hLayerSumHit);
0532 m_hLayerSumHit->Fill(layer, sumHit);
0533
0534 double sumData = 0;
0535 for (const int& data : layerChanDataSize[layer][side])
0536 {
0537 sumData += data;
0538
0539 assert(m_hLayerDataSize);
0540 m_hLayerDataSize->Fill(layer, data);
0541 }
0542 assert(m_hLayerSumDataSize);
0543 m_hLayerSumDataSize->Fill(layer, sumData);
0544 }
0545
0546 }
0547
0548 assert(m_hWavelet);
0549 m_hWavelet->Fill(nWavelet);
0550 h_norm->Fill("TPC Wavelet", nWavelet);
0551 assert(m_hDataSize);
0552 m_hDataSize->Fill(sumDataSize);
0553 h_norm->Fill("TPC DataSize", sumDataSize);
0554
0555 return Fun4AllReturnCodes::EVENT_OK;
0556 }
0557
0558 int TPCDataStreamEmulator::writeWavelet(int layer, int side, int phibin, int hittime, const vector<unsigned int>& wavelet)
0559 {
0560 static const int headersize = 2;
0561
0562
0563 const int datasizebit = wavelet.size() * 10;
0564 int datasizebyte = datasizebit / 8;
0565 if (datasizebyte * 8 < datasizebit) datasizebyte += 1;
0566
0567 assert(m_hLayerWaveletSize);
0568 m_hLayerWaveletSize->Fill(layer, wavelet.size());
0569
0570 return headersize + datasizebyte;
0571 }
0572
0573 Fun4AllHistoManager*
0574 TPCDataStreamEmulator::getHistoManager()
0575 {
0576 static string histname("TPCDataStreamEmulator_HISTOS");
0577
0578 Fun4AllServer* se = Fun4AllServer::instance();
0579 Fun4AllHistoManager* hm = se->getHistoManager(histname);
0580
0581 if (not hm)
0582 {
0583 cout
0584 << "TPCDataStreamEmulator::get_HistoManager - Making Fun4AllHistoManager " << histname
0585 << endl;
0586 hm = new Fun4AllHistoManager(histname);
0587 se->registerHistoManager(hm);
0588 }
0589
0590 assert(hm);
0591
0592 return hm;
0593 }
