Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TPCDataStreamEmulator.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TPCDataStreamEmulator.cc
1 // $Id: $
2 
11 #include "TPCDataStreamEmulator.h"
12 
13 #include "TPCDaqDefs.h"
14 
15 #include <g4detectors/PHG4Cell.h>
16 #include <g4detectors/PHG4CellContainer.h>
17 #include <g4detectors/PHG4CylinderCellGeom.h>
18 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
19 #include <trackbase_historic/SvtxHit.h>
20 #include <trackbase_historic/SvtxHitMap.h>
21 #include <g4main/PHG4Hit.h>
22 #include <g4main/PHG4HitContainer.h>
23 #include <g4main/PHG4Particle.h>
24 #include <g4main/PHG4TruthInfoContainer.h>
25 
26 #include <fun4all/Fun4AllHistoManager.h>
27 #include <fun4all/Fun4AllReturnCodes.h>
28 #include <fun4all/Fun4AllServer.h>
29 #include <fun4all/PHTFileServer.h>
30 
31 #include <phhepmc/PHHepMCGenEvent.h>
32 #include <phhepmc/PHHepMCGenEventMap.h>
33 #include <phool/PHCompositeNode.h>
34 #include <phool/getClass.h>
35 
36 #include <TDatabasePDG.h>
37 #include <TFile.h>
38 #include <TH1D.h>
39 #include <TH2D.h>
40 #include <TString.h>
41 #include <TTree.h>
42 #include <TVector3.h>
43 
44 #include <HepMC/GenEvent.h>
45 
46 #include <CLHEP/Units/SystemOfUnits.h>
47 
48 #include <boost/format.hpp>
49 
50 #include <algorithm>
51 #include <array>
52 #include <cassert>
53 #include <cmath>
54 #include <iostream>
55 #include <limits>
56 #include <stdexcept>
57 
58 using namespace std;
59 using namespace CLHEP;
60 
61 TPCDataStreamEmulator::TPCDataStreamEmulator(
62  unsigned int minLayer,
63  unsigned int m_maxLayer,
64  const std::string& outputfilename)
65  : SubsysReco("TPCDataStreamEmulator")
66  , m_saveDataStreamFile(true)
67  , m_outputFileNameBase(outputfilename)
68  , m_minLayer(minLayer)
69  , m_maxLayer(m_maxLayer)
70  , m_evtCounter(-1)
71  , m_vertexZAcceptanceCut(10)
72  , m_etaAcceptanceCut(1.1)
73  , m_hDataSize(nullptr)
74  , m_hWavelet(nullptr)
75  , m_hNChEta(nullptr)
76  , m_hLayerWaveletSize(nullptr)
77  , m_hLayerHit(nullptr)
78  , m_hLayerZBinHit(nullptr)
79  , m_hLayerZBinADC(nullptr)
80  , m_hLayerDataSize(nullptr)
81  , m_hLayerSumHit(nullptr)
82  , m_hLayerSumDataSize(nullptr)
83 {
84 }
85 
86 TPCDataStreamEmulator::~TPCDataStreamEmulator()
87 {
88 }
89 
90 int TPCDataStreamEmulator::Init(PHCompositeNode* topNode)
91 {
92  return Fun4AllReturnCodes::EVENT_OK;
93 }
94 
95 int TPCDataStreamEmulator::End(PHCompositeNode* topNode)
96 {
97  if (Verbosity() >= VERBOSITY_SOME)
98  cout << "TPCDataStreamEmulator::End - write to " << m_outputFileNameBase + ".root" << endl;
99  PHTFileServer::get().cd(m_outputFileNameBase + ".root");
100 
101  Fun4AllHistoManager* hm = getHistoManager();
102  assert(hm);
103  for (unsigned int i = 0; i < hm->nHistos(); i++)
104  hm->getHisto(i)->Write();
105 
106  // help index files with TChain
107  TTree* T_Index = new TTree("T_Index", "T_Index");
108  assert(T_Index);
109  T_Index->Write();
110 
111  return Fun4AllReturnCodes::EVENT_OK;
112 }
113 
114 int TPCDataStreamEmulator::InitRun(PHCompositeNode* topNode)
115 {
116  PHG4CylinderCellGeomContainer* seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
117  if (!seggeo)
118  {
119  cout << "could not locate geo node "
120  << "CYLINDERCELLGEOM_SVTX" << endl;
121  exit(1);
122  }
123 
124  int nZBins = 0;
125  for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
126  {
127  PHG4CylinderCellGeom* layerGeom =
128  seggeo->GetLayerCellGeom(layer);
129  assert(layerGeom);
130 
131  if (nZBins <= 0)
132  {
133  nZBins = layerGeom->get_zbins();
134  assert(nZBins > 0);
135  }
136  else
137  {
138  if ((int) nZBins != layerGeom->get_zbins())
139  {
140  cout << "TPCDataStreamEmulator::InitRun - Fatal Error - nZBin at layer " << layer << " is " << layerGeom->get_zbins()
141  << ", which is different from previous layers of nZBin = " << nZBins << endl;
142  exit(1);
143  }
144  }
145  } // for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
146 
147  if (Verbosity() >= VERBOSITY_SOME)
148  cout << "TPCDataStreamEmulator::get_HistoManager - Making PHTFileServer " << m_outputFileNameBase + ".root"
149  << endl;
150  PHTFileServer::get().open(m_outputFileNameBase + ".root", "RECREATE");
151 
152  Fun4AllHistoManager* hm = getHistoManager();
153  assert(hm);
154 
155  TH1D* h = new TH1D("hNormalization", //
156  "Normalization;Items;Summed quantity", 10, .5, 10.5);
157  int i = 1;
158  h->GetXaxis()->SetBinLabel(i++, "Event count");
159  h->GetXaxis()->SetBinLabel(i++, "Collision count");
160  h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit");
161  h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit Edep");
162  h->GetXaxis()->SetBinLabel(i++, "TPC Hit");
163  h->GetXaxis()->SetBinLabel(i++, "TPC Wavelet");
164  h->GetXaxis()->SetBinLabel(i++, "TPC DataSize");
165 
166  h->GetXaxis()->LabelsOption("v");
167  hm->registerHisto(h);
168 
169  // for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
170  // {
171  // const PHG4CylinderCellGeom* layer_geom = seggeo->GetLayerCellGeom(layer);
172 
173  // const string histNameCellHit(boost::str(boost::format{"hCellHit_Layer%1%"} % layer));
174  // const string histNameCellCharge(boost::str(boost::format{"hCellCharge_Layer%1%"} % layer));
175 
176  // }
177 
178  hm->registerHisto(m_hDataSize =
179  new TH1D("hDataSize", //
180  "TPC Data Size per Event;Data size [Byte];Count",
181  10000, 0, 20e6));
182 
183  hm->registerHisto(m_hWavelet =
184  new TH1D("hWavelet", //
185  "TPC Recorded Wavelet per Event;Wavelet count;Count",
186  10000, 0, 4e6));
187 
188  hm->registerHisto(m_hNChEta =
189  new TH1D("hNChEta", //
190  "Charged particle #eta distribution;#eta;Count",
191  1000, -5, 5));
192 
193  hm->registerHisto(m_hLayerWaveletSize =
194  new TH2D("hLayerWaveletSize", //
195  "Number of Recorded ADC sample per Wavelet;Layer ID;ADC Sample Count per Wavelet",
196  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
197  nZBins, -.5, nZBins - .5));
198 
199  hm->registerHisto(m_hLayerHit =
200  new TH2D("hLayerHit", //
201  "Number of Recorded ADC sample per channel;Layer ID;ADC Sample Count",
202  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
203  nZBins, -.5, nZBins - .5));
204 
205  hm->registerHisto(m_hLayerDataSize =
206  new TH2D("hLayerDataSize", //
207  "Data size per channel;Layer ID;Data size [Byte]",
208  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
209  2 * nZBins, 0, 2 * nZBins));
210 
211  hm->registerHisto(m_hLayerSumHit =
212  new TH2D("hLayerSumHit", //
213  "Number of Recorded ADC sample per layer;Layer ID;ADC Sample Count",
214  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
215  10000, -.5, 99999.5));
216 
217  hm->registerHisto(m_hLayerSumDataSize =
218  new TH2D("hLayerSumDataSize", //
219  "Data size per trigger per layer;Layer ID;Data size [Byte]",
220  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
221  1000, 0, .5e6));
222 
223  hm->registerHisto(m_hLayerZBinHit =
224  new TH2D("hLayerZBinHit", //
225  "Number of Recorded ADC sample per Time Bin;z bin ID;Layer ID",
226  nZBins, -.5, nZBins - .5,
227  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5));
228 
229  hm->registerHisto(m_hLayerZBinADC =
230  new TH2D("hLayerZBinADC", //
231  "Sum ADC per Time Bin;z bin ID;Layer ID",
232  nZBins, -.5, nZBins - .5,
233  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5));
234  return Fun4AllReturnCodes::EVENT_OK;
235 }
236 
237 int TPCDataStreamEmulator::process_event(PHCompositeNode* topNode)
238 {
239  m_evtCounter += 1;
240 
241  Fun4AllHistoManager* hm = getHistoManager();
242  assert(hm);
243  TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto("hNormalization"));
244  assert(h_norm);
245  h_norm->Fill("Event count", 1);
246 
247  PHG4HitContainer* g4hit = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_SVTX");
248  if (!g4hit)
249  {
250  cout << "TPCDataStreamEmulator::process_event - Could not locate g4 hit node G4HIT_SVTX" << endl;
251  return Fun4AllReturnCodes::ABORTRUN;
252  }
253 
254  SvtxHitMap* hits = findNode::getClass<SvtxHitMap>(topNode, "SvtxHitMap");
255  if (!hits)
256  {
257  cout << "PCDataStreamEmulator::process_event - ERROR: Can't find node SvtxHitMap" << endl;
258  return Fun4AllReturnCodes::ABORTRUN;
259  }
260 
261  PHG4CellContainer* cells = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_SVTX");
262  if (!cells)
263  {
264  cout << "TPCDataStreamEmulator::process_event - could not locate cell node "
265  << "G4CELL_SVTX" << endl;
266  exit(1);
267  }
268 
269  PHG4CylinderCellGeomContainer* seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
270  if (!seggeo)
271  {
272  cout << "TPCDataStreamEmulator::process_event - could not locate geo node "
273  << "CYLINDERCELLGEOM_SVTX" << endl;
274  exit(1);
275  }
276 
277  PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
278  if (!geneventmap)
279  {
280  static bool once = true;
281  if (once)
282  {
283  once = false;
284 
285  cout << "TPCDataStreamEmulator::process_event - - missing node PHHepMCGenEventMap. Skipping HepMC stat." << std::endl;
286  }
287  }
288  else
289  {
290  h_norm->Fill("Collision count", geneventmap->size());
291  }
292 
293  PHG4TruthInfoContainer* truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode,
294  "G4TruthInfo");
295  if (!truthInfoList)
296  {
297  cout << "TPCDataStreamEmulator::process_event - Fatal Error - "
298  << "unable to find DST node "
299  << "G4TruthInfo" << endl;
300  exit(1);
301  }
302 
303  PHG4TruthInfoContainer::ConstRange primary_range =
304  truthInfoList->GetPrimaryParticleRange();
305 
306  for (PHG4TruthInfoContainer::ConstIterator particle_iter =
307  primary_range.first;
308  particle_iter != primary_range.second;
309  ++particle_iter)
310  {
311  const PHG4Particle* p = particle_iter->second;
312  assert(p);
313 
314  TParticlePDG* pdg_p = TDatabasePDG::Instance()->GetParticle(
315  p->get_pid());
316  assert(pdg_p);
317 
318  if (fabs(pdg_p->Charge()) > 0)
319  {
320  TVector3 pvec(p->get_px(), p->get_py(), p->get_pz());
321 
322  if (pvec.Perp2()>0)
323  {
324  assert(m_hNChEta);
325  m_hNChEta->Fill(pvec.PseudoRapidity());
326  }
327  }
328 
329  } // if (_load_all_particle) else
330 
331  for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
332  {
333  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(layer);
334 
335  for (PHG4HitContainer::ConstIterator hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
336  {
337  const double edep = hiter->second->get_edep();
338  h_norm->Fill("TPC G4Hit Edep", edep);
339  h_norm->Fill("TPC G4Hit", 1);
340  } // for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
341 
342  } // for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
343 
344  // prepreare stat. storage
345  int nZBins = 0;
346  vector<array<vector<int>, 2> > layerChanHit(m_maxLayer + 1);
347  vector<array<vector<int>, 2> > layerChanDataSize(m_maxLayer + 1);
348  int nWavelet = 0;
349  int sumDataSize = 0;
350  for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
351  {
352  PHG4CylinderCellGeom* layerGeom =
353  seggeo->GetLayerCellGeom(layer);
354  assert(layerGeom);
355 
356  // start with an empty vector of cells for each phibin
357  const int nphibins = layerGeom->get_phibins();
358  assert(nphibins > 0);
359 
360  if (Verbosity() >= VERBOSITY_MORE)
361  {
362  cout << "TPCDataStreamEmulator::process_event - init layer " << layer << " with "
363  << "nphibins = " << nphibins
364  << ", layerGeom->get_zbins() = " << layerGeom->get_zbins() << endl;
365  }
366 
367  if (nZBins <= 0)
368  {
369  nZBins = layerGeom->get_zbins();
370  assert(nZBins > 0);
371  }
372  else
373  {
374  if ((int) nZBins != layerGeom->get_zbins())
375  {
376  cout << "TPCDataStreamEmulator::process_event - Fatal Error - nZBin at layer " << layer << " is " << layerGeom->get_zbins()
377  << ", which is different from previous layers of nZBin = " << nZBins << endl;
378  exit(1);
379  }
380  }
381 
382  for (unsigned int side = 0; side < 2; ++side)
383  {
384  layerChanHit[layer][side].resize(nphibins, 0);
385 
386  layerChanDataSize[layer][side].resize(nphibins, 0);
387  } // for (unsigned int side = 0; side < 2; ++side)
388 
389  } // for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
390 
391  assert(nZBins > 0);
392 
393  // count hits and make wavelets
394  int last_layer = -1;
395  int last_side = -1;
396  int last_phibin = -1;
397  int last_zbin = -1;
398  vector<unsigned int> last_wavelet;
399  int last_wavelet_hittime = -1;
400 
401  for (SvtxHitMap::Iter iter = hits->begin(); iter != hits->end(); ++iter)
402  {
403  SvtxHit* hit = iter->second;
404 
405  const int layer = hit->get_layer();
406 
407  if (layer < m_minLayer or layer > m_maxLayer) continue;
408 
409  PHG4Cell* cell = cells->findCell(hit->get_cellid()); //not needed once geofixed
410  const int phibin = PHG4CellDefs::SizeBinning::get_phibin(cell->get_cellid()); //cell->get_binphi();
411  const int zbin = PHG4CellDefs::SizeBinning::get_zbin(cell->get_cellid()); //cell->get_binz();
412  const int side = (zbin < nZBins / 2) ? 0 : 1;
413 
414  // new wavelet?
415  if (last_layer != layer or last_phibin != phibin or last_side != side or abs(last_zbin - zbin) != 1)
416  {
417  // save last wavelet
418  if (last_wavelet.size() > 0)
419  {
420  const int datasize = writeWavelet(last_layer, last_side, last_phibin, last_wavelet_hittime, last_wavelet);
421  assert(datasize > 0);
422 
423  nWavelet += 1;
424  sumDataSize += datasize;
425  layerChanDataSize[last_layer][last_side][last_phibin] += datasize;
426 
427  last_wavelet.clear();
428  last_zbin = -1;
429  }
430 
431  // z-R cut on digitized wavelet
432  PHG4CylinderCellGeom* layerGeom =
433  seggeo->GetLayerCellGeom(layer);
434  assert(layerGeom);
435  const double z_abs = fabs(layerGeom->get_zcenter(zbin));
436  const double r = layerGeom->get_radius();
437  TVector3 acceptanceVec(r, 0, z_abs - m_vertexZAcceptanceCut);
438  const double eta = acceptanceVec.PseudoRapidity();
439 
440  if (eta > m_etaAcceptanceCut) continue;
441 
442  // make new wavelet
443  last_layer = layer;
444  last_side = side;
445  last_phibin = phibin;
446 
447  // time check
448  last_wavelet_hittime = (side == 0) ? (zbin) : (nZBins - 1 - zbin);
449  assert(last_wavelet_hittime >= 0);
450  assert(last_wavelet_hittime <= nZBins / 2);
451  } // if (last_layer != layer or last_phibin != phibin)
452 
453  if (Verbosity() >= VERBOSITY_A_LOT)
454  {
455  cout << "TPCDataStreamEmulator::process_event - layer " << layer << " hit with "
456 
457  << "phibin = " << phibin
458  << ",zbin = " << zbin
459  << ",side = " << side
460  << ",last_wavelet.size() = " << last_wavelet.size()
461  << ",last_zbin = " << last_zbin
462  << endl;
463  }
464 
465  // more checks on signal continuity
466  if (last_wavelet.size() > 0)
467  {
468  if (side == 0)
469  {
470  assert(zbin - last_zbin == 1);
471  }
472  else
473  {
474  assert(last_zbin - zbin == 1);
475  }
476  }
477 
478  // record adc
479  unsigned int adc = hit->get_adc();
480  last_wavelet.push_back(adc);
481  last_zbin = zbin;
482 
483  // statistics
484  layerChanHit[layer][side][phibin] += 1;
485  assert(m_hLayerZBinHit);
486  m_hLayerZBinHit->Fill(zbin, layer, 1);
487  assert(m_hLayerZBinADC);
488  m_hLayerZBinADC->Fill(zbin, layer, adc);
489 
490  } // for(SvtxHitMap::Iter iter = hits->begin(); iter != hits->end(); ++iter) {
491 
492  // save last wavelet
493  if (last_wavelet.size() > 0)
494  {
495  const int datasize = writeWavelet(last_layer, last_side, last_phibin, last_wavelet_hittime, last_wavelet);
496  assert(datasize > 0);
497 
498  nWavelet += 1;
499  sumDataSize += datasize;
500  layerChanDataSize[last_layer][last_side][last_phibin] += datasize;
501  }
502 
503  // statistics
504  for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
505  {
506  for (unsigned int side = 0; side < 2; ++side)
507  {
508  int sumHit = 0;
509  for (const int& hit : layerChanHit[layer][side])
510  {
511  sumHit += hit;
512 
513  assert(m_hLayerHit);
514  m_hLayerHit->Fill(layer, hit);
515  h_norm->Fill("TPC Hit", hit);
516 
517  if (Verbosity() >= VERBOSITY_MORE)
518  {
519  cout << "TPCDataStreamEmulator::process_event - hit: "
520  << "hit = " << hit
521  << "at layer = " << layer
522  << ", side = " << side
523  << endl;
524  }
525  }
526 
527  if (Verbosity() >= VERBOSITY_MORE)
528  {
529  cout << "TPCDataStreamEmulator::process_event - hLayerSumCellHit->Fill(" << layer << ", " << sumHit << ")" << endl;
530  }
531  assert(m_hLayerSumHit);
532  m_hLayerSumHit->Fill(layer, sumHit);
533 
534  double sumData = 0;
535  for (const int& data : layerChanDataSize[layer][side])
536  {
537  sumData += data;
538 
539  assert(m_hLayerDataSize);
540  m_hLayerDataSize->Fill(layer, data);
541  }
542  assert(m_hLayerSumDataSize);
543  m_hLayerSumDataSize->Fill(layer, sumData);
544  } // for (unsigned int side = 0; side < 2; ++side)
545 
546  } // for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
547 
548  assert(m_hWavelet);
549  m_hWavelet->Fill(nWavelet);
550  h_norm->Fill("TPC Wavelet", nWavelet);
551  assert(m_hDataSize);
552  m_hDataSize->Fill(sumDataSize);
553  h_norm->Fill("TPC DataSize", sumDataSize);
554 
555  return Fun4AllReturnCodes::EVENT_OK;
556 }
557 
558 int TPCDataStreamEmulator::writeWavelet(int layer, int side, int phibin, int hittime, const vector<unsigned int>& wavelet)
559 {
560  static const int headersize = 2; // 2-byte header per wavelet
561 
562  //data in byte aligned and padding
563  const int datasizebit = wavelet.size() * 10;
564  int datasizebyte = datasizebit / 8;
565  if (datasizebyte * 8 < datasizebit) datasizebyte += 1;
566 
567  assert(m_hLayerWaveletSize);
568  m_hLayerWaveletSize->Fill(layer, wavelet.size());
569 
570  return headersize + datasizebyte;
571 }
572 
573 Fun4AllHistoManager*
574 TPCDataStreamEmulator::getHistoManager()
575 {
576  static string histname("TPCDataStreamEmulator_HISTOS");
577 
578  Fun4AllServer* se = Fun4AllServer::instance();
579  Fun4AllHistoManager* hm = se->getHistoManager(histname);
580 
581  if (not hm)
582  {
583  cout
584  << "TPCDataStreamEmulator::get_HistoManager - Making Fun4AllHistoManager " << histname
585  << endl;
586  hm = new Fun4AllHistoManager(histname);
587  se->registerHistoManager(hm);
588  }
589 
590  assert(hm);
591 
592  return hm;
593 }