Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TPCMLDataInterface.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TPCMLDataInterface.cc
1 // $Id: $
2 
11 #include "TPCMLDataInterface.h"
12 
13 #include <tpc/TpcDefs.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 <g4eval/SvtxEvalStack.h>
27 #include <g4eval/SvtxEvaluator.h>
28 #include <g4eval/SvtxHitEval.h>
29 #include <g4eval/SvtxTruthEval.h>
30 
31 #include <trackbase/TrkrClusterContainer.h>
32 #include <trackbase/TrkrClusterHitAssoc.h>
33 #include <trackbase/TrkrClusterv1.h>
34 #include <trackbase/TrkrHit.h>
35 #include <trackbase/TrkrHitSet.h>
36 #include <trackbase/TrkrHitSetContainer.h>
37 
38 #include <fun4all/Fun4AllHistoManager.h>
39 #include <fun4all/Fun4AllReturnCodes.h>
40 #include <fun4all/Fun4AllServer.h>
41 #include <fun4all/PHTFileServer.h>
42 
43 #include <phhepmc/PHHepMCGenEvent.h>
44 #include <phhepmc/PHHepMCGenEventMap.h>
45 #include <phool/PHCompositeNode.h>
46 #include <phool/getClass.h>
47 
48 #include <TDatabasePDG.h>
49 #include <TFile.h>
50 #include <TH1D.h>
51 #include <TH2D.h>
52 #include <TString.h>
53 #include <TTree.h>
54 #include <TVector3.h>
55 
56 #include <HepMC/GenEvent.h>
57 
58 #include <CLHEP/Units/SystemOfUnits.h>
59 
60 #include <H5Cpp.h>
61 
62 #include <boost/format.hpp>
63 
64 #include <algorithm>
65 #include <array>
66 #include <cassert>
67 #include <cmath>
68 #include <iostream>
69 #include <limits>
70 #include <stdexcept>
71 
72 using namespace std;
73 using namespace CLHEP;
74 using namespace H5;
75 
76 TPCMLDataInterface::TPCMLDataInterface(
77  unsigned int minLayer,
78  unsigned int m_maxLayer,
79  const std::string& outputfilename)
80  : SubsysReco("TPCMLDataInterface")
81  , m_svtxevalstack(nullptr)
82  , m_strict(false)
83  , m_saveDataStreamFile(true)
84  , m_outputFileNameBase(outputfilename)
85  , m_h5File(nullptr)
86  , m_minLayer(minLayer)
87  , m_maxLayer(m_maxLayer)
88  , m_evtCounter(-1)
89  , m_vertexZAcceptanceCut(10)
90  , m_etaAcceptanceCut(1.1)
91  , m_momentumCut(.1)
92  , m_hDataSize(nullptr)
93  , m_hWavelet(nullptr)
94  , m_hNChEta(nullptr)
95  , m_hLayerWaveletSize(nullptr)
96  , m_hLayerHit(nullptr)
97  , m_hLayerZBinHit(nullptr)
98  , m_hLayerZBinADC(nullptr)
99  , m_hLayerDataSize(nullptr)
100  , m_hLayerSumHit(nullptr)
101  , m_hLayerSumDataSize(nullptr)
102 {
103 }
104 
105 TPCMLDataInterface::~TPCMLDataInterface()
106 {
107  if (m_h5File) delete m_h5File;
108 }
109 
110 int TPCMLDataInterface::Init(PHCompositeNode* topNode)
111 {
112  return Fun4AllReturnCodes::EVENT_OK;
113 }
114 
115 int TPCMLDataInterface::End(PHCompositeNode* topNode)
116 {
117  if (Verbosity() >= VERBOSITY_SOME)
118  cout << "TPCMLDataInterface::End - write to " << m_outputFileNameBase + ".root" << endl;
119  PHTFileServer::get().cd(m_outputFileNameBase + ".root");
120 
121  Fun4AllHistoManager* hm = getHistoManager();
122  assert(hm);
123  for (unsigned int i = 0; i < hm->nHistos(); i++)
124  hm->getHisto(i)->Write();
125 
126  // help index files with TChain
127  TTree* T_Index = new TTree("T_Index", "T_Index");
128  assert(T_Index);
129  T_Index->Write();
130 
131  if (m_h5File)
132  {
133  delete m_h5File;
134  m_h5File = nullptr;
135  }
136 
137  return Fun4AllReturnCodes::EVENT_OK;
138 }
139 
140 int TPCMLDataInterface::InitRun(PHCompositeNode* topNode)
141 {
142  PHG4CylinderCellGeomContainer* seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
143  if (!seggeo)
144  {
145  cout << "could not locate geo node "
146  << "CYLINDERCELLGEOM_SVTX" << endl;
147  exit(1);
148  }
149 
150  int nZBins = 0;
151  for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
152  {
153  PHG4CylinderCellGeom* layerGeom =
154  seggeo->GetLayerCellGeom(layer);
155  assert(layerGeom);
156 
157  if (nZBins <= 0)
158  {
159  nZBins = layerGeom->get_zbins();
160  assert(nZBins > 0);
161  }
162  else
163  {
164  if ((int) nZBins != layerGeom->get_zbins())
165  {
166  cout << "TPCMLDataInterface::InitRun - Fatal Error - nZBin at layer " << layer << " is " << layerGeom->get_zbins()
167  << ", which is different from previous layers of nZBin = " << nZBins << endl;
168  exit(1);
169  }
170  }
171  } // for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
172 
173  if (Verbosity() >= VERBOSITY_SOME)
174  cout << "TPCMLDataInterface::get_HistoManager - Making PHTFileServer " << m_outputFileNameBase + ".root"
175  << endl;
176  PHTFileServer::get().open(m_outputFileNameBase + ".root", "RECREATE");
177 
178  Fun4AllHistoManager* hm = getHistoManager();
179  assert(hm);
180 
181  TH1D* h = new TH1D("hNormalization", //
182  "Normalization;Items;Summed quantity", 10, .5, 10.5);
183  int i = 1;
184  h->GetXaxis()->SetBinLabel(i++, "Event count");
185  h->GetXaxis()->SetBinLabel(i++, "Collision count");
186  h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit");
187  h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit Edep");
188  h->GetXaxis()->SetBinLabel(i++, "TPC Hit");
189  h->GetXaxis()->SetBinLabel(i++, "TPC Wavelet");
190  h->GetXaxis()->SetBinLabel(i++, "TPC DataSize");
191 
192  h->GetXaxis()->LabelsOption("v");
193  hm->registerHisto(h);
194 
195  // for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
196  // {
197  // const PHG4CylinderCellGeom* layer_geom = seggeo->GetLayerCellGeom(layer);
198 
199  // const string histNameCellHit(boost::str(boost::format{"hCellHit_Layer%1%"} % layer));
200  // const string histNameCellCharge(boost::str(boost::format{"hCellCharge_Layer%1%"} % layer));
201 
202  // }
203 
204  hm->registerHisto(m_hDataSize =
205  new TH1D("hDataSize", //
206  "TPC Data Size per Event;Data size [Byte];Count",
207  10000, 0, 20e6));
208 
209  hm->registerHisto(m_hWavelet =
210  new TH1D("hWavelet", //
211  "TPC Recorded Wavelet per Event;Wavelet count;Count",
212  10000, 0, 4e6));
213 
214  hm->registerHisto(m_hNChEta =
215  new TH1D("hNChEta", //
216  "Charged particle #eta distribution;#eta;Count",
217  1000, -5, 5));
218 
219  hm->registerHisto(m_hLayerWaveletSize =
220  new TH2D("hLayerWaveletSize", //
221  "Number of Recorded ADC sample per Wavelet;Layer ID;ADC Sample Count per Wavelet",
222  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
223  nZBins, -.5, nZBins - .5));
224 
225  hm->registerHisto(m_hLayerHit =
226  new TH2D("hLayerHit", //
227  "Number of Recorded ADC sample per channel;Layer ID;ADC Sample Count",
228  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
229  nZBins, -.5, nZBins - .5));
230 
231  hm->registerHisto(m_hLayerDataSize =
232  new TH2D("hLayerDataSize", //
233  "Data size per channel;Layer ID;Data size [Byte]",
234  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
235  2 * nZBins, 0, 2 * nZBins));
236 
237  hm->registerHisto(m_hLayerSumHit =
238  new TH2D("hLayerSumHit", //
239  "Number of Recorded ADC sample per layer;Layer ID;ADC Sample Count",
240  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
241  10000, -.5, 99999.5));
242 
243  hm->registerHisto(m_hLayerSumDataSize =
244  new TH2D("hLayerSumDataSize", //
245  "Data size per trigger per layer;Layer ID;Data size [Byte]",
246  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
247  1000, 0, .5e6));
248 
249  hm->registerHisto(m_hLayerZBinHit =
250  new TH2D("hLayerZBinHit", //
251  "Number of Recorded ADC sample per Time Bin;z bin ID;Layer ID",
252  nZBins, -.5, nZBins - .5,
253  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5));
254 
255  hm->registerHisto(m_hLayerZBinADC =
256  new TH2D("hLayerZBinADC", //
257  "Sum ADC per Time Bin;z bin ID;Layer ID",
258  nZBins, -.5, nZBins - .5,
259  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5));
260 
261  // init HDF5 output
262 
263  if (Verbosity() >= VERBOSITY_SOME)
264  cout << "TPCMLDataInterface::get_HistoManager - Making H5File " << m_outputFileNameBase + ".h5"
265  << endl;
266  m_h5File = new H5File(m_outputFileNameBase + ".h5", H5F_ACC_TRUNC);
267 
268  return Fun4AllReturnCodes::EVENT_OK;
269 }
270 
271 int TPCMLDataInterface::process_event(PHCompositeNode* topNode)
272 {
273  m_evtCounter += 1;
274 
275  assert(m_h5File);
276 
277  Fun4AllHistoManager* hm = getHistoManager();
278  assert(hm);
279  TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto("hNormalization"));
280  assert(h_norm);
281  h_norm->Fill("Event count", 1);
282 
283  if (!m_svtxevalstack)
284  {
285  m_svtxevalstack = new SvtxEvalStack(topNode);
286  m_svtxevalstack->set_strict(m_strict);
287  m_svtxevalstack->set_verbosity(Verbosity());
288  m_svtxevalstack->set_use_initial_vertex(m_use_initial_vertex);
289  m_svtxevalstack->set_use_genfit_vertex(m_use_genfit_vertex);
290  m_svtxevalstack->next_event(topNode);
291  }
292  else
293  {
294  m_svtxevalstack->next_event(topNode);
295  }
296 
297  SvtxHitEval* hiteval = m_svtxevalstack->get_hit_eval();
298  SvtxTruthEval* trutheval = m_svtxevalstack->get_truth_eval();
299 
300  float gpx = 0;
301  float gpy = 0;
302  float gpz = 0;
303 
304  PHG4HitContainer* g4hit = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
305  cout << "TPCMLDataInterface::process_event - g4 hit node G4HIT_TPC" << endl;
306  if (!g4hit)
307  {
308  cout << "TPCMLDataInterface::process_event - Could not locate g4 hit node G4HIT_TPC" << endl;
309  return Fun4AllReturnCodes::ABORTRUN;
310  }
311 
312  // get node containing the digitized hits
313  TrkrHitSetContainer* hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
314  if (!hits)
315  {
316  cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << endl;
317  return Fun4AllReturnCodes::ABORTRUN;
318  }
319 
320  // PHG4CellContainer* cells = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_SVTX");
321  // if (!cells)
322  // {
323  // cout << "TPCMLDataInterface::process_event - could not locate cell node "
324  // << "G4CELL_SVTX" << endl;
325  // exit(1);
326  // }
327 
328  PHG4CylinderCellGeomContainer* seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
329  if (!seggeo)
330  {
331  cout << "TPCMLDataInterface::process_event - could not locate geo node "
332  << "CYLINDERCELLGEOM_SVTX" << endl;
333  exit(1);
334  }
335 
336  PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
337  if (!geneventmap)
338  {
339  static bool once = true;
340  if (once)
341  {
342  once = false;
343 
344  cout << "TPCMLDataInterface::process_event - - missing node PHHepMCGenEventMap. Skipping HepMC stat." << std::endl;
345  }
346  }
347  else
348  {
349  h_norm->Fill("Collision count", geneventmap->size());
350  }
351 
352  PHG4TruthInfoContainer* truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode,
353  "G4TruthInfo");
354  if (!truthInfoList)
355  {
356  cout << "TPCMLDataInterface::process_event - Fatal Error - "
357  << "unable to find DST node "
358  << "G4TruthInfo" << endl;
359  assert(truthInfoList);
360  }
361 
362  PHG4TruthInfoContainer::ConstRange primary_range =
363  truthInfoList->GetPrimaryParticleRange();
364 
365  for (PHG4TruthInfoContainer::ConstIterator particle_iter =
366  primary_range.first;
367  particle_iter != primary_range.second;
368  ++particle_iter)
369  {
370  const PHG4Particle* p = particle_iter->second;
371  assert(p);
372 
373  TParticlePDG* pdg_p = TDatabasePDG::Instance()->GetParticle(
374  p->get_pid());
375  assert(pdg_p);
376 
377  if (fabs(pdg_p->Charge()) > 0)
378  {
379  TVector3 pvec(p->get_px(), p->get_py(), p->get_pz());
380 
381  if (pvec.Perp2() > 0)
382  {
383  assert(m_hNChEta);
384  m_hNChEta->Fill(pvec.PseudoRapidity());
385  }
386  }
387 
388  } // if (_load_all_particle) else
389 
390  for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
391  {
392  PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(layer);
393 
394  for (PHG4HitContainer::ConstIterator hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
395  {
396  const double edep = hiter->second->get_edep();
397  h_norm->Fill("TPC G4Hit Edep", edep);
398  h_norm->Fill("TPC G4Hit", 1);
399  } // for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
400 
401  } // for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
402 
403  //H5 init
404  string h5GroupName = boost::str(boost::format("TimeFrame_%1%") % m_evtCounter);
405  unique_ptr<Group> h5Group(new Group(m_h5File->createGroup("/" + h5GroupName)));
406  h5Group->setComment(boost::str(boost::format("Collection of ADC data matrixes in Time Frame #%1%") % m_evtCounter));
407  if (Verbosity())
408  cout << "TPCMLDataInterface::process_event - save to H5 group " << h5GroupName << endl;
409  map<int, shared_ptr<DataSet>> layerH5DataSetMap;
410  map<int, shared_ptr<DataSet>> layerH5SignalBackgroundMap;
411  map<int, shared_ptr<DataSpace>> layerH5DataSpaceMap;
412  map<int, shared_ptr<DataSpace>> layerH5SignalBackgroundDataSpaceMap;
413 
414  //wavelet size stat.
415  vector<uint32_t> layerWaveletDataSize(m_maxLayer + 1, 0);
416  vector<hsize_t> layerSize(1);
417  layerSize[0] = static_cast<hsize_t>(m_maxLayer - m_minLayer + 1);
418  shared_ptr<DataSpace> H5DataSpaceLayerWaveletDataSize(new DataSpace(1, layerSize.data()));
419  shared_ptr<DataSet> H5DataSetLayerWaveletDataSize(new DataSet(h5Group->createDataSet(
420  "sPHENIXRawDataSizeBytePerLayer",
421  PredType::NATIVE_UINT32,
422  *(H5DataSpaceLayerWaveletDataSize))));
423 
424  // prepreare stat. storage
425  int nZBins = 0;
426  vector<array<vector<int>, 2>> layerChanHit(m_maxLayer + 1);
427  vector<array<vector<int>, 2>> layerChanDataSize(m_maxLayer + 1);
428  vector<vector<uint16_t>> layerDataBuffer(m_maxLayer + 1);
429  vector<vector<uint8_t>> layerSignalBackgroundBuffer(m_maxLayer + 1);
430  vector<vector<hsize_t>> layerDataBufferSize(m_maxLayer + 1, vector<hsize_t>({0, 0}));
431  vector<vector<hsize_t>> layerSignalBackgroundDataBufferSize(m_maxLayer + 1, vector<hsize_t>({0, 0}));
432 
433  int nWavelet = 0;
434  int sumDataSize = 0;
435  for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
436  {
437  PHG4CylinderCellGeom* layerGeom =
438  seggeo->GetLayerCellGeom(layer);
439  assert(layerGeom);
440 
441  // start with an empty vector of cells for each phibin
442  const int nphibins = layerGeom->get_phibins();
443  assert(nphibins > 0);
444 
445  if (Verbosity() >= VERBOSITY_MORE)
446  {
447  cout << "TPCMLDataInterface::process_event - init layer " << layer << " with "
448  << "nphibins = " << nphibins
449  << ", layerGeom->get_zbins() = " << layerGeom->get_zbins() << endl;
450  }
451 
452  if (nZBins <= 0)
453  {
454  nZBins = layerGeom->get_zbins();
455  assert(nZBins > 0);
456  }
457  else
458  {
459  if ((int) nZBins != layerGeom->get_zbins())
460  {
461  cout << "TPCMLDataInterface::process_event - Fatal Error - nZBin at layer " << layer << " is " << layerGeom->get_zbins()
462  << ", which is different from previous layers of nZBin = " << nZBins << endl;
463  exit(1);
464  }
465  }
466 
467  for (unsigned int side = 0; side < 2; ++side)
468  {
469  layerChanHit[layer][side].resize(nphibins, 0);
470 
471  layerChanDataSize[layer][side].resize(nphibins, 0);
472  } // for (unsigned int side = 0; side < 2; ++side)
473 
474  // buffer init
475  layerDataBufferSize[layer][0] = nphibins;
476  layerDataBufferSize[layer][1] = layerGeom->get_zbins();
477  layerSignalBackgroundDataBufferSize[layer][0] = nphibins;
478  layerSignalBackgroundDataBufferSize[layer][1] = layerGeom->get_zbins();
479  layerDataBuffer[layer].resize(layerDataBufferSize[layer][0] * layerDataBufferSize[layer][1], 0);
480  layerSignalBackgroundBuffer[layer].resize(layerSignalBackgroundDataBufferSize[layer][0] * layerSignalBackgroundDataBufferSize[layer][1], 0);
481 
482  static const vector<hsize_t> cdims({32, 32});
483  DSetCreatPropList ds_creatplist;
484  ds_creatplist.setChunk(2, cdims.data()); // then modify it for compression
485  ds_creatplist.setDeflate(6);
486 
487  layerH5DataSpaceMap[layer] = shared_ptr<DataSpace>(new DataSpace(2, layerDataBufferSize[layer].data()));
488  layerH5SignalBackgroundDataSpaceMap[layer] = shared_ptr<DataSpace>(new DataSpace(2, layerSignalBackgroundDataBufferSize[layer].data()));
489 
490  layerH5DataSetMap[layer] = shared_ptr<DataSet>(new DataSet(h5Group->createDataSet(
491  boost::str(boost::format("Data_Layer%1%") % (layer - m_minLayer)),
492  PredType::NATIVE_UINT16,
493  *(layerH5DataSpaceMap[layer]),
494  ds_creatplist)));
495 
496  layerH5SignalBackgroundMap[layer] = shared_ptr<DataSet>(new DataSet(h5Group->createDataSet(
497  boost::str(boost::format("Data_Layer_SignalBackground%1%") % (layer - m_minLayer)),
498  PredType::NATIVE_UINT8,
499  *(layerH5SignalBackgroundDataSpaceMap[layer]),
500  ds_creatplist)));
501 
502  } // for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
503 
504  assert(nZBins > 0);
505 
506  // count hits and make wavelets
507  int last_layer = -1;
508  int last_side = -1;
509  int last_phibin = -1;
510  int last_zbin = -1;
511  vector<unsigned int> last_wavelet;
512  int last_wavelet_hittime = -1;
513 
514  // for (SvtxHitMap::Iter iter = hits->begin(); iter != hits->end(); ++iter)
515  // {
516  // SvtxHit* hit = iter->second;
517  // loop over the TPC HitSet objects
518  TrkrHitSetContainer::ConstRange hitsetrange = hits->getHitSets(TrkrDefs::TrkrId::tpcId);
519  for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
520  hitsetitr != hitsetrange.second;
521  ++hitsetitr)
522  {
523  TrkrHitSet* hitset = hitsetitr->second;
524 
525  if (Verbosity() > 2) hitset->identify();
526 
527  // const int layer = hit->get_layer();
528  // we have a single hitset, get the info that identifies the module
529  const int layer = TrkrDefs::getLayer(hitsetitr->first);
530 
531  if (Verbosity() >= 2)
532  {
533  cout << "TPCMLDataInterface::process_event - TrkrHitSet @ layer " << layer << " with " << hitset->size() << " hits" << endl;
534  }
535 
536  if (layer < m_minLayer or layer > m_maxLayer) continue;
537 
538  // PHG4Cell* cell = cells->findCell(hit->get_cellid()); //not needed once geofixed
539 
540  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
541  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
542  hitr != hitrangei.second;
543  ++hitr)
544  {
545  // const int phibin = PHG4CellDefs::SizeBinning::get_phibin(cell->get_cellid()); //cell->get_binphi();
546  // const int zbin = PHG4CellDefs::SizeBinning::get_zbin(cell->get_cellid()); //cell->get_binz();
547  int phibin = TpcDefs::getPad(hitr->first);
548  int zbin = TpcDefs::getTBin(hitr->first);
549  const int side = (zbin < nZBins / 2) ? 0 : 1;
550 
551  if (Verbosity() >= 2)
552  {
553  cout << "TPCMLDataInterface::process_event - hit @ layer " << layer << " phibin = " << phibin << " zbin = " << zbin << " side = " << side << endl;
554  }
555 
556  assert(phibin >= 0);
557  assert(zbin >= 0);
558  assert(side >= 0);
559 
560  // new wavelet?
561  if (last_layer != layer or last_phibin != phibin or last_side != side or abs(last_zbin - zbin) != 1)
562  {
563  // save last wavelet
564  if (last_wavelet.size() > 0)
565  {
566  const int datasize = writeWavelet(last_layer, last_side, last_phibin, last_wavelet_hittime, last_wavelet);
567  assert(datasize > 0);
568 
569  nWavelet += 1;
570  sumDataSize += datasize;
571  layerChanDataSize[last_layer][last_side][last_phibin] += datasize;
572 
573  last_wavelet.clear();
574  last_zbin = -1;
575  }
576 
577  // z-R cut on digitized wavelet
578  PHG4CylinderCellGeom* layerGeom =
579  seggeo->GetLayerCellGeom(layer);
580  assert(layerGeom);
581 
582  const double z_abs = fabs(layerGeom->get_zcenter(zbin));
583  const double r = layerGeom->get_radius();
584  TVector3 acceptanceVec(r, 0, z_abs - m_vertexZAcceptanceCut);
585  const double eta = acceptanceVec.PseudoRapidity();
586 
587  if (eta > m_etaAcceptanceCut) continue;
588 
589  // make new wavelet
590  last_layer = layer;
591  last_side = side;
592  last_phibin = phibin;
593 
594  // time check
595  last_wavelet_hittime = (side == 0) ? (zbin) : (nZBins - 1 - zbin);
596  assert(last_wavelet_hittime >= 0);
597  assert(last_wavelet_hittime <= nZBins / 2);
598 
599  } // if (last_layer != layer or last_phibin != phibin)
600 
601  if (Verbosity() >= VERBOSITY_A_LOT)
602  {
603  cout << "TPCMLDataInterface::process_event - layer " << layer << " hit with "
604 
605  << "phibin = " << phibin
606  << ",zbin = " << zbin
607  << ",side = " << side
608  << ",last_wavelet.size() = " << last_wavelet.size()
609  << ",last_zbin = " << last_zbin
610  << endl;
611  }
612 
613  // more checks on signal continuity
614  if (last_wavelet.size() > 0)
615  {
616  // if (side == 0)
617  // {
618  assert(zbin - last_zbin == 1);
619  // }
620  // else
621  // {
622  // assert(last_zbin - zbin == 1);
623  // }
624  }
625 
626  // record adc
627  unsigned int adc = hitr->second->getAdc();
628  last_wavelet.push_back(adc);
629  last_zbin = zbin;
630 
631  // statistics
632  layerChanHit[layer][side][phibin] += 1;
633  assert(m_hLayerZBinHit);
634  m_hLayerZBinHit->Fill(zbin, layer, 1);
635  assert(m_hLayerZBinADC);
636  m_hLayerZBinADC->Fill(zbin, layer, adc);
637  assert(adc < (1 << 10));
638  assert((hsize_t) phibin < layerDataBufferSize[layer][0]);
639  assert((hsize_t) zbin < layerDataBufferSize[layer][1]);
640  const size_t hitindex(layerDataBufferSize[layer][1] * phibin + zbin);
641  assert((hsize_t) phibin < layerSignalBackgroundDataBufferSize[layer][0]);
642  assert((hsize_t) zbin < layerSignalBackgroundDataBufferSize[layer][1]);
643  const size_t hitindexSB(layerSignalBackgroundDataBufferSize[layer][1] * phibin + zbin);
644  assert(hitindex < layerDataBuffer[layer].size());
645  assert(hitindexSB < layerSignalBackgroundBuffer[layer].size());
646 
647  if (layerDataBuffer[layer][hitindex] != 0)
648  {
649  cout << "TPCMLDataInterface::process_event - WARNING - hit @ layer "
650  << layer << " phibin = " << phibin << " zbin = " << zbin << " side = " << side
651  << " overwriting previous hit with ADC = " << layerDataBuffer[layer][hitindex]
652  << endl;
653  }
654  layerDataBuffer[layer][hitindex] = adc;
655 
656  if (layerSignalBackgroundBuffer[layer][hitindexSB] != 0)
657  {
658  cout << "TPCMLDataInterface::process_event - WARNING - signal/background @ layer "
659  << layer << " phibin = " << phibin << " zbin = " << zbin << " side = " << side
660  << " overwriting previous hit with = " << layerSignalBackgroundBuffer[layer][hitindexSB]
661  << endl;
662  }
663 
664  // momentum cut
665 
666  bool signal_or_bgd = false;
667 
668  TrkrDefs::hitkey hit_key = hitr->first;
669  // PHG4Hit* g4hit = hiteval->max_truth_hit_by_energy(hit_key);
670  std::set<PHG4Hit*> g4hit_set = hiteval->all_truth_hits(hit_key);
671 
672  for (PHG4Hit* g4hit : g4hit_set)
673  {
674  PHG4Particle* g4particle = trutheval->get_particle(g4hit);
675 
676  if (g4particle != nullptr)
677  {
678  gpx = g4particle->get_px();
679  gpy = g4particle->get_py();
680  gpz = g4particle->get_pz();
681  }
682 
683  if (sqrt(gpx * gpx + gpy * gpy + gpz * gpz) > m_momentumCut) signal_or_bgd = true;
684 
685  } // for (auto g4hit: g4hit_set)
686 
687  layerSignalBackgroundBuffer[layer][hitindexSB] = signal_or_bgd;
688 
689  } // for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
690 
691  } // for(SvtxHitMap::Iter iter = hits->begin(); iter != hits->end(); ++iter) {
692 
693  // save last wavelet
694  if (last_wavelet.size() > 0)
695  {
696  const int datasize = writeWavelet(last_layer, last_side, last_phibin, last_wavelet_hittime, last_wavelet);
697  assert(datasize > 0);
698 
699  nWavelet += 1;
700  sumDataSize += datasize;
701  layerChanDataSize[last_layer][last_side][last_phibin] += datasize;
702  }
703 
704  // statistics
705  for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
706  {
707  for (unsigned int side = 0; side < 2; ++side)
708  {
709  int sumHit = 0;
710  for (const int& hit : layerChanHit[layer][side])
711  {
712  sumHit += hit;
713 
714  assert(m_hLayerHit);
715  m_hLayerHit->Fill(layer, hit);
716  h_norm->Fill("TPC Hit", hit);
717 
718  if (Verbosity() >= VERBOSITY_MORE)
719  {
720  cout << "TPCMLDataInterface::process_event - layerChanHit: "
721  << "hit = " << hit
722  << "at layer = " << layer
723  << ", side = " << side
724  << endl;
725  }
726  }
727 
728  if (Verbosity() >= VERBOSITY_MORE)
729  {
730  cout << "TPCMLDataInterface::process_event - hLayerSumCellHit->Fill(" << layer << ", " << sumHit << ")" << endl;
731  }
732  assert(m_hLayerSumHit);
733  m_hLayerSumHit->Fill(layer, sumHit);
734 
735  double sumData = 0;
736  for (const int& data : layerChanDataSize[layer][side])
737  {
738  sumData += data;
739 
740  assert(m_hLayerDataSize);
741  m_hLayerDataSize->Fill(layer, data);
742  }
743 
744  assert(m_hLayerSumDataSize);
745  m_hLayerSumDataSize->Fill(layer, sumData);
746  layerWaveletDataSize[(layer - m_minLayer)] += sumData;
747  } // for (unsigned int side = 0; side < 2; ++side)
748 
749  // store in H5
750 
751  assert(layerH5DataSetMap[layer]);
752  assert(layerH5SignalBackgroundMap[layer]);
753  layerH5DataSetMap[layer]->write(layerDataBuffer[layer].data(), PredType::NATIVE_UINT16);
754  layerH5SignalBackgroundMap[layer]->write(layerSignalBackgroundBuffer[layer].data(), PredType::NATIVE_UINT8);
755  H5DataSetLayerWaveletDataSize->write(layerWaveletDataSize.data(), PredType::NATIVE_UINT32);
756 
757  } // for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
758 
759  assert(m_hWavelet);
760  m_hWavelet->Fill(nWavelet);
761  h_norm->Fill("TPC Wavelet", nWavelet);
762  assert(m_hDataSize);
763  m_hDataSize->Fill(sumDataSize);
764  h_norm->Fill("TPC DataSize", sumDataSize);
765 
766  return Fun4AllReturnCodes::EVENT_OK;
767 }
768 
769 int TPCMLDataInterface::writeWavelet(int layer, int side, int phibin, int hittime, const vector<unsigned int>& wavelet)
770 {
771  static const int headersize = 2; // 2-byte header per wavelet
772 
773  //data in byte aligned and padding
774  const int datasizebit = wavelet.size() * 10;
775  int datasizebyte = datasizebit / 8;
776  if (datasizebyte * 8 < datasizebit) datasizebyte += 1;
777 
778  assert(m_hLayerWaveletSize);
779  m_hLayerWaveletSize->Fill(layer, wavelet.size());
780 
781  return headersize + datasizebyte;
782 }
783 
784 Fun4AllHistoManager*
785 TPCMLDataInterface::getHistoManager()
786 {
787  static string histname("TPCMLDataInterface_HISTOS");
788 
789  Fun4AllServer* se = Fun4AllServer::instance();
790  Fun4AllHistoManager* hm = se->getHistoManager(histname);
791 
792  if (not hm)
793  {
794  cout
795  << "TPCMLDataInterface::get_HistoManager - Making Fun4AllHistoManager " << histname
796  << endl;
797  hm = new Fun4AllHistoManager(histname);
798  se->registerHistoManager(hm);
799  }
800 
801  assert(hm);
802 
803  return hm;
804 }