DAQ/TPCMLDataInterface/TPCMLDataInterface.cc

0001 // $Id: $
0002
0003 /*!
0004  * \file TPCMLDataInterface.cpp
0005  * \brief
0006  * \author Jin Huang <jhuang@bnl.gov>
0007  * \version $Revision:   $
0008  * \date $Date: $
0009  */
0010
0011 #include "TPCMLDataInterface.h"
0012
0013 #include <tpc/TpcDefs.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 <g4eval/SvtxEvalStack.h>
0027 #include <g4eval/SvtxEvaluator.h>
0028 #include <g4eval/SvtxHitEval.h>
0029 #include <g4eval/SvtxTruthEval.h>
0030
0031 #include <trackbase/TrkrClusterContainer.h>
0032 #include <trackbase/TrkrClusterHitAssoc.h>
0033 #include <trackbase/TrkrClusterv1.h>
0034 #include <trackbase/TrkrHit.h>
0035 #include <trackbase/TrkrHitSet.h>
0036 #include <trackbase/TrkrHitSetContainer.h>
0037
0038 #include <fun4all/Fun4AllHistoManager.h>
0039 #include <fun4all/Fun4AllReturnCodes.h>
0040 #include <fun4all/Fun4AllServer.h>
0041 #include <fun4all/PHTFileServer.h>
0042
0043 #include <phhepmc/PHHepMCGenEvent.h>
0044 #include <phhepmc/PHHepMCGenEventMap.h>
0045 #include <phool/PHCompositeNode.h>
0046 #include <phool/getClass.h>
0047
0048 #include <TDatabasePDG.h>
0049 #include <TFile.h>
0050 #include <TH1D.h>
0051 #include <TH2D.h>
0052 #include <TString.h>
0053 #include <TTree.h>
0054 #include <TVector3.h>
0055
0056 #include <HepMC/GenEvent.h>
0057
0058 #include <CLHEP/Units/SystemOfUnits.h>
0059
0060 #include <H5Cpp.h>
0061
0062 #include <boost/format.hpp>
0063
0064 #include <algorithm>
0065 #include <array>
0066 #include <cassert>
0067 #include <cmath>
0068 #include <iostream>
0069 #include <limits>
0070 #include <stdexcept>
0071
0072 using namespace std;
0073 using namespace CLHEP;
0074 using namespace H5;
0075
0076 TPCMLDataInterface::TPCMLDataInterface(
0077     unsigned int minLayer,
0078     unsigned int m_maxLayer,
0079     const std::string& outputfilename)
0080   : SubsysReco("TPCMLDataInterface")
0081   , m_svtxevalstack(nullptr)
0082   , m_strict(false)
0083   , m_saveDataStreamFile(true)
0084   , m_outputFileNameBase(outputfilename)
0085   , m_h5File(nullptr)
0086   , m_minLayer(minLayer)
0087   , m_maxLayer(m_maxLayer)
0088   , m_evtCounter(-1)
0089   , m_vertexZAcceptanceCut(10)
0090   , m_etaAcceptanceCut(1.1)
0091   , m_momentumCut(.1)
0092   , m_hDataSize(nullptr)
0093   , m_hWavelet(nullptr)
0094   , m_hNChEta(nullptr)
0095   , m_hLayerWaveletSize(nullptr)
0096   , m_hLayerHit(nullptr)
0097   , m_hLayerZBinHit(nullptr)
0098   , m_hLayerZBinADC(nullptr)
0099   , m_hLayerDataSize(nullptr)
0100   , m_hLayerSumHit(nullptr)
0101   , m_hLayerSumDataSize(nullptr)
0102 {
0103 }
0104
0105 TPCMLDataInterface::~TPCMLDataInterface()
0106 {
0107   if (m_h5File) delete m_h5File;
0108 }
0109
0110 int TPCMLDataInterface::Init(PHCompositeNode* topNode)
0111 {
0112   return Fun4AllReturnCodes::EVENT_OK;
0113 }
0114
0115 int TPCMLDataInterface::End(PHCompositeNode* topNode)
0116 {
0117   if (Verbosity() >= VERBOSITY_SOME)
0118     cout << "TPCMLDataInterface::End - write to " << m_outputFileNameBase + ".root" << endl;
0119   PHTFileServer::get().cd(m_outputFileNameBase + ".root");
0120
0121   Fun4AllHistoManager* hm = getHistoManager();
0122   assert(hm);
0123   for (unsigned int i = 0; i < hm->nHistos(); i++)
0124     hm->getHisto(i)->Write();
0125
0126   // help index files with TChain
0127   TTree* T_Index = new TTree("T_Index", "T_Index");
0128   assert(T_Index);
0129   T_Index->Write();
0130
0131   if (m_h5File)
0132   {
0133     delete m_h5File;
0134     m_h5File = nullptr;
0135   }
0136
0137   return Fun4AllReturnCodes::EVENT_OK;
0138 }
0139
0140 int TPCMLDataInterface::InitRun(PHCompositeNode* topNode)
0141 {
0142   PHG4CylinderCellGeomContainer* seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0143   if (!seggeo)
0144   {
0145     cout << "could not locate geo node "
0146          << "CYLINDERCELLGEOM_SVTX" << endl;
0147     exit(1);
0148   }
0149
0150   int nZBins = 0;
0151   for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0152   {
0153     PHG4CylinderCellGeom* layerGeom =
0154         seggeo->GetLayerCellGeom(layer);
0155     assert(layerGeom);
0156
0157     if (nZBins <= 0)
0158     {
0159       nZBins = layerGeom->get_zbins();
0160       assert(nZBins > 0);
0161     }
0162     else
0163     {
0164       if ((int) nZBins != layerGeom->get_zbins())
0165       {
0166         cout << "TPCMLDataInterface::InitRun - Fatal Error - nZBin at layer " << layer << " is " << layerGeom->get_zbins()
0167              << ", which is different from previous layers of nZBin = " << nZBins << endl;
0168         exit(1);
0169       }
0170     }
0171   }  //   for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0172
0173   if (Verbosity() >= VERBOSITY_SOME)
0174     cout << "TPCMLDataInterface::get_HistoManager - Making PHTFileServer " << m_outputFileNameBase + ".root"
0175          << endl;
0176   PHTFileServer::get().open(m_outputFileNameBase + ".root", "RECREATE");
0177
0178   Fun4AllHistoManager* hm = getHistoManager();
0179   assert(hm);
0180
0181   TH1D* h = new TH1D("hNormalization",  //
0182                      "Normalization;Items;Summed quantity", 10, .5, 10.5);
0183   int i = 1;
0184   h->GetXaxis()->SetBinLabel(i++, "Event count");
0185   h->GetXaxis()->SetBinLabel(i++, "Collision count");
0186   h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit");
0187   h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit Edep");
0188   h->GetXaxis()->SetBinLabel(i++, "TPC Hit");
0189   h->GetXaxis()->SetBinLabel(i++, "TPC Wavelet");
0190   h->GetXaxis()->SetBinLabel(i++, "TPC DataSize");
0191
0192   h->GetXaxis()->LabelsOption("v");
0193   hm->registerHisto(h);
0194
0195   //  for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0196   //  {
0197   //    const PHG4CylinderCellGeom* layer_geom = seggeo->GetLayerCellGeom(layer);
0198
0199   //    const string histNameCellHit(boost::str(boost::format{"hCellHit_Layer%1%"} % layer));
0200   //    const string histNameCellCharge(boost::str(boost::format{"hCellCharge_Layer%1%"} % layer));
0201
0202   //  }
0203
0204   hm->registerHisto(m_hDataSize =
0205                         new TH1D("hDataSize",  //
0206                                  "TPC Data Size per Event;Data size [Byte];Count",
0207                                  10000, 0, 20e6));
0208
0209   hm->registerHisto(m_hWavelet =
0210                         new TH1D("hWavelet",  //
0211                                  "TPC Recorded Wavelet per Event;Wavelet count;Count",
0212                                  10000, 0, 4e6));
0213
0214   hm->registerHisto(m_hNChEta =
0215                         new TH1D("hNChEta",  //
0216                                  "Charged particle #eta distribution;#eta;Count",
0217                                  1000, -5, 5));
0218
0219   hm->registerHisto(m_hLayerWaveletSize =
0220                         new TH2D("hLayerWaveletSize",  //
0221                                  "Number of Recorded ADC sample per Wavelet;Layer ID;ADC Sample Count per Wavelet",
0222                                  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0223                                  nZBins, -.5, nZBins - .5));
0224
0225   hm->registerHisto(m_hLayerHit =
0226                         new TH2D("hLayerHit",  //
0227                                  "Number of Recorded ADC sample per channel;Layer ID;ADC Sample Count",
0228                                  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0229                                  nZBins, -.5, nZBins - .5));
0230
0231   hm->registerHisto(m_hLayerDataSize =
0232                         new TH2D("hLayerDataSize",  //
0233                                  "Data size per channel;Layer ID;Data size [Byte]",
0234                                  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0235                                  2 * nZBins, 0, 2 * nZBins));
0236
0237   hm->registerHisto(m_hLayerSumHit =
0238                         new TH2D("hLayerSumHit",  //
0239                                  "Number of Recorded ADC sample per layer;Layer ID;ADC Sample Count",
0240                                  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0241                                  10000, -.5, 99999.5));
0242
0243   hm->registerHisto(m_hLayerSumDataSize =
0244                         new TH2D("hLayerSumDataSize",  //
0245                                  "Data size per trigger per layer;Layer ID;Data size [Byte]",
0246                                  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0247                                  1000, 0, .5e6));
0248
0249   hm->registerHisto(m_hLayerZBinHit =
0250                         new TH2D("hLayerZBinHit",  //
0251                                  "Number of Recorded ADC sample per Time Bin;z bin ID;Layer ID",
0252                                  nZBins, -.5, nZBins - .5,
0253                                  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5));
0254
0255   hm->registerHisto(m_hLayerZBinADC =
0256                         new TH2D("hLayerZBinADC",  //
0257                                  "Sum ADC per Time Bin;z bin ID;Layer ID",
0258                                  nZBins, -.5, nZBins - .5,
0259                                  m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5));
0260
0261   // init HDF5 output
0262
0263   if (Verbosity() >= VERBOSITY_SOME)
0264     cout << "TPCMLDataInterface::get_HistoManager - Making H5File " << m_outputFileNameBase + ".h5"
0265          << endl;
0266   m_h5File = new H5File(m_outputFileNameBase + ".h5", H5F_ACC_TRUNC);
0267
0268   return Fun4AllReturnCodes::EVENT_OK;
0269 }
0270
0271 int TPCMLDataInterface::process_event(PHCompositeNode* topNode)
0272 {
0273   m_evtCounter += 1;
0274
0275   assert(m_h5File);
0276
0277   Fun4AllHistoManager* hm = getHistoManager();
0278   assert(hm);
0279   TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto("hNormalization"));
0280   assert(h_norm);
0281   h_norm->Fill("Event count", 1);
0282
0283   if (!m_svtxevalstack)
0284   {
0285     m_svtxevalstack = new SvtxEvalStack(topNode);
0286     m_svtxevalstack->set_strict(m_strict);
0287     m_svtxevalstack->set_verbosity(Verbosity());
0288     m_svtxevalstack->set_use_initial_vertex(m_use_initial_vertex);
0289     m_svtxevalstack->set_use_genfit_vertex(m_use_genfit_vertex);
0290     m_svtxevalstack->next_event(topNode);
0291   }
0292   else
0293   {
0294     m_svtxevalstack->next_event(topNode);
0295   }
0296
0297   SvtxHitEval* hiteval = m_svtxevalstack->get_hit_eval();
0298   SvtxTruthEval* trutheval = m_svtxevalstack->get_truth_eval();
0299
0300   float gpx = 0;
0301   float gpy = 0;
0302   float gpz = 0;
0303
0304   PHG4HitContainer* g4hit = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
0305   cout << "TPCMLDataInterface::process_event - g4 hit node G4HIT_TPC" << endl;
0306   if (!g4hit)
0307   {
0308     cout << "TPCMLDataInterface::process_event - Could not locate g4 hit node G4HIT_TPC" << endl;
0309     return Fun4AllReturnCodes::ABORTRUN;
0310   }
0311
0312   // get node containing the digitized hits
0313   TrkrHitSetContainer* hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0314   if (!hits)
0315   {
0316     cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << endl;
0317     return Fun4AllReturnCodes::ABORTRUN;
0318   }
0319
0320   //  PHG4CellContainer* cells = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_SVTX");
0321   //  if (!cells)
0322   //  {
0323   //    cout << "TPCMLDataInterface::process_event - could not locate cell node "
0324   //         << "G4CELL_SVTX" << endl;
0325   //    exit(1);
0326   //  }
0327
0328   PHG4CylinderCellGeomContainer* seggeo = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0329   if (!seggeo)
0330   {
0331     cout << "TPCMLDataInterface::process_event - could not locate geo node "
0332          << "CYLINDERCELLGEOM_SVTX" << endl;
0333     exit(1);
0334   }
0335
0336   PHHepMCGenEventMap* geneventmap = findNode::getClass<PHHepMCGenEventMap>(topNode, "PHHepMCGenEventMap");
0337   if (!geneventmap)
0338   {
0339     static bool once = true;
0340     if (once)
0341     {
0342       once = false;
0343
0344       cout << "TPCMLDataInterface::process_event - - missing node PHHepMCGenEventMap. Skipping HepMC stat." << std::endl;
0345     }
0346   }
0347   else
0348   {
0349     h_norm->Fill("Collision count", geneventmap->size());
0350   }
0351
0352   PHG4TruthInfoContainer* truthInfoList = findNode::getClass<PHG4TruthInfoContainer>(topNode,
0353                                                                                      "G4TruthInfo");
0354   if (!truthInfoList)
0355   {
0356     cout << "TPCMLDataInterface::process_event - Fatal Error - "
0357          << "unable to find DST node "
0358          << "G4TruthInfo" << endl;
0359     assert(truthInfoList);
0360   }
0361
0362   PHG4TruthInfoContainer::ConstRange primary_range =
0363       truthInfoList->GetPrimaryParticleRange();
0364
0365   for (PHG4TruthInfoContainer::ConstIterator particle_iter =
0366            primary_range.first;
0367        particle_iter != primary_range.second;
0368        ++particle_iter)
0369   {
0370     const PHG4Particle* p = particle_iter->second;
0371     assert(p);
0372
0373     TParticlePDG* pdg_p = TDatabasePDG::Instance()->GetParticle(
0374         p->get_pid());
0375     assert(pdg_p);
0376
0377     if (fabs(pdg_p->Charge()) > 0)
0378     {
0379       TVector3 pvec(p->get_px(), p->get_py(), p->get_pz());
0380
0381       if (pvec.Perp2() > 0)
0382       {
0383         assert(m_hNChEta);
0384         m_hNChEta->Fill(pvec.PseudoRapidity());
0385       }
0386     }
0387
0388   }  //          if (_load_all_particle) else
0389
0390   for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0391   {
0392     PHG4HitContainer::ConstRange hit_begin_end = g4hit->getHits(layer);
0393
0394     for (PHG4HitContainer::ConstIterator hiter = hit_begin_end.first; hiter != hit_begin_end.second; ++hiter)
0395     {
0396       const double edep = hiter->second->get_edep();
0397       h_norm->Fill("TPC G4Hit Edep", edep);
0398       h_norm->Fill("TPC G4Hit", 1);
0399     }  //     for (hiter = hit_begin_end.first; hiter != hit_begin_end.second; hiter++)
0400
0401   }  //   for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0402
0403   //H5 init
0404   string h5GroupName = boost::str(boost::format("TimeFrame_%1%") % m_evtCounter);
0405   unique_ptr<Group> h5Group(new Group(m_h5File->createGroup("/" + h5GroupName)));
0406   h5Group->setComment(boost::str(boost::format("Collection of ADC data matrixes in Time Frame #%1%") % m_evtCounter));
0407   if (Verbosity())
0408     cout << "TPCMLDataInterface::process_event - save to H5 group " << h5GroupName << endl;
0409   map<int, shared_ptr<DataSet>> layerH5DataSetMap;
0410   map<int, shared_ptr<DataSet>> layerH5SignalBackgroundMap;
0411   map<int, shared_ptr<DataSpace>> layerH5DataSpaceMap;
0412   map<int, shared_ptr<DataSpace>> layerH5SignalBackgroundDataSpaceMap;
0413
0414   //wavelet size stat.
0415   vector<uint32_t> layerWaveletDataSize(m_maxLayer + 1, 0);
0416   vector<hsize_t> layerSize(1);
0417   layerSize[0] = static_cast<hsize_t>(m_maxLayer - m_minLayer + 1);
0418   shared_ptr<DataSpace> H5DataSpaceLayerWaveletDataSize(new DataSpace(1, layerSize.data()));
0419   shared_ptr<DataSet> H5DataSetLayerWaveletDataSize(new DataSet(h5Group->createDataSet(
0420       "sPHENIXRawDataSizeBytePerLayer",
0421       PredType::NATIVE_UINT32,
0422       *(H5DataSpaceLayerWaveletDataSize))));
0423
0424   // prepreare stat. storage
0425   int nZBins = 0;
0426   vector<array<vector<int>, 2>> layerChanHit(m_maxLayer + 1);
0427   vector<array<vector<int>, 2>> layerChanDataSize(m_maxLayer + 1);
0428   vector<vector<uint16_t>> layerDataBuffer(m_maxLayer + 1);
0429   vector<vector<uint8_t>> layerSignalBackgroundBuffer(m_maxLayer + 1);
0430   vector<vector<hsize_t>> layerDataBufferSize(m_maxLayer + 1, vector<hsize_t>({0, 0}));
0431   vector<vector<hsize_t>> layerSignalBackgroundDataBufferSize(m_maxLayer + 1, vector<hsize_t>({0, 0}));
0432
0433   int nWavelet = 0;
0434   int sumDataSize = 0;
0435   for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0436   {
0437     PHG4CylinderCellGeom* layerGeom =
0438         seggeo->GetLayerCellGeom(layer);
0439     assert(layerGeom);
0440
0441     // start with an empty vector of cells for each phibin
0442     const int nphibins = layerGeom->get_phibins();
0443     assert(nphibins > 0);
0444
0445     if (Verbosity() >= VERBOSITY_MORE)
0446     {
0447       cout << "TPCMLDataInterface::process_event - init layer " << layer << " with "
0448            << "nphibins = " << nphibins
0449            << ", layerGeom->get_zbins() = " << layerGeom->get_zbins() << endl;
0450     }
0451
0452     if (nZBins <= 0)
0453     {
0454       nZBins = layerGeom->get_zbins();
0455       assert(nZBins > 0);
0456     }
0457     else
0458     {
0459       if ((int) nZBins != layerGeom->get_zbins())
0460       {
0461         cout << "TPCMLDataInterface::process_event - Fatal Error - nZBin at layer " << layer << " is " << layerGeom->get_zbins()
0462              << ", which is different from previous layers of nZBin = " << nZBins << endl;
0463         exit(1);
0464       }
0465     }
0466
0467     for (unsigned int side = 0; side < 2; ++side)
0468     {
0469       layerChanHit[layer][side].resize(nphibins, 0);
0470
0471       layerChanDataSize[layer][side].resize(nphibins, 0);
0472     }  //     for (unsigned int side = 0; side < 2; ++side)
0473
0474     // buffer init
0475     layerDataBufferSize[layer][0] = nphibins;
0476     layerDataBufferSize[layer][1] = layerGeom->get_zbins();
0477     layerSignalBackgroundDataBufferSize[layer][0] = nphibins;
0478     layerSignalBackgroundDataBufferSize[layer][1] = layerGeom->get_zbins();
0479     layerDataBuffer[layer].resize(layerDataBufferSize[layer][0] * layerDataBufferSize[layer][1], 0);
0480     layerSignalBackgroundBuffer[layer].resize(layerSignalBackgroundDataBufferSize[layer][0] * layerSignalBackgroundDataBufferSize[layer][1], 0);
0481
0482     static const vector<hsize_t> cdims({32, 32});
0483     DSetCreatPropList ds_creatplist;
0484     ds_creatplist.setChunk(2, cdims.data());  // then modify it for compression
0485     ds_creatplist.setDeflate(6);
0486
0487     layerH5DataSpaceMap[layer] = shared_ptr<DataSpace>(new DataSpace(2, layerDataBufferSize[layer].data()));
0488     layerH5SignalBackgroundDataSpaceMap[layer] = shared_ptr<DataSpace>(new DataSpace(2, layerSignalBackgroundDataBufferSize[layer].data()));
0489
0490     layerH5DataSetMap[layer] = shared_ptr<DataSet>(new DataSet(h5Group->createDataSet(
0491         boost::str(boost::format("Data_Layer%1%") % (layer - m_minLayer)),
0492         PredType::NATIVE_UINT16,
0493         *(layerH5DataSpaceMap[layer]),
0494         ds_creatplist)));
0495
0496     layerH5SignalBackgroundMap[layer] = shared_ptr<DataSet>(new DataSet(h5Group->createDataSet(
0497         boost::str(boost::format("Data_Layer_SignalBackground%1%") % (layer - m_minLayer)),
0498         PredType::NATIVE_UINT8,
0499         *(layerH5SignalBackgroundDataSpaceMap[layer]),
0500         ds_creatplist)));
0501
0502   }  //   for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0503
0504   assert(nZBins > 0);
0505
0506   // count hits and make wavelets
0507   int last_layer = -1;
0508   int last_side = -1;
0509   int last_phibin = -1;
0510   int last_zbin = -1;
0511   vector<unsigned int> last_wavelet;
0512   int last_wavelet_hittime = -1;
0513
0514   //  for (SvtxHitMap::Iter iter = hits->begin(); iter != hits->end(); ++iter)
0515   //  {
0516   //    SvtxHit* hit = iter->second;
0517   // loop over the TPC HitSet objects
0518   TrkrHitSetContainer::ConstRange hitsetrange = hits->getHitSets(TrkrDefs::TrkrId::tpcId);
0519   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0520        hitsetitr != hitsetrange.second;
0521        ++hitsetitr)
0522   {
0523     TrkrHitSet* hitset = hitsetitr->second;
0524
0525     if (Verbosity() > 2) hitset->identify();
0526
0527     //    const int layer = hit->get_layer();
0528     // we have a single hitset, get the info that identifies the module
0529     const int layer = TrkrDefs::getLayer(hitsetitr->first);
0530
0531     if (Verbosity() >= 2)
0532     {
0533       cout << "TPCMLDataInterface::process_event - TrkrHitSet @ layer " << layer << " with " << hitset->size() << " hits" << endl;
0534     }
0535
0536     if (layer < m_minLayer or layer > m_maxLayer) continue;
0537
0538     //    PHG4Cell* cell = cells->findCell(hit->get_cellid());                           //not needed once geofixed
0539
0540     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0541     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0542          hitr != hitrangei.second;
0543          ++hitr)
0544     {
0545       //    const int phibin = PHG4CellDefs::SizeBinning::get_phibin(cell->get_cellid());  //cell->get_binphi();
0546       //    const int zbin = PHG4CellDefs::SizeBinning::get_zbin(cell->get_cellid());      //cell->get_binz();
0547       int phibin = TpcDefs::getPad(hitr->first);
0548       int zbin = TpcDefs::getTBin(hitr->first);
0549       const int side = (zbin < nZBins / 2) ? 0 : 1;
0550
0551       if (Verbosity() >= 2)
0552       {
0553         cout << "TPCMLDataInterface::process_event - hit @ layer " << layer << " phibin = " << phibin << " zbin = " << zbin << " side = " << side << endl;
0554       }
0555
0556       assert(phibin >= 0);
0557       assert(zbin >= 0);
0558       assert(side >= 0);
0559
0560       // new wavelet?
0561       if (last_layer != layer or last_phibin != phibin or last_side != side or abs(last_zbin - zbin) != 1)
0562       {
0563         // save last wavelet
0564         if (last_wavelet.size() > 0)
0565         {
0566           const int datasize = writeWavelet(last_layer, last_side, last_phibin, last_wavelet_hittime, last_wavelet);
0567           assert(datasize > 0);
0568
0569           nWavelet += 1;
0570           sumDataSize += datasize;
0571           layerChanDataSize[last_layer][last_side][last_phibin] += datasize;
0572
0573           last_wavelet.clear();
0574           last_zbin = -1;
0575         }
0576
0577         // z-R cut on digitized wavelet
0578         PHG4CylinderCellGeom* layerGeom =
0579             seggeo->GetLayerCellGeom(layer);
0580         assert(layerGeom);
0581
0582         const double z_abs = fabs(layerGeom->get_zcenter(zbin));
0583         const double r = layerGeom->get_radius();
0584         TVector3 acceptanceVec(r, 0, z_abs - m_vertexZAcceptanceCut);
0585         const double eta = acceptanceVec.PseudoRapidity();
0586
0587         if (eta > m_etaAcceptanceCut) continue;
0588
0589         // make new wavelet
0590         last_layer = layer;
0591         last_side = side;
0592         last_phibin = phibin;
0593
0594         // time check
0595         last_wavelet_hittime = (side == 0) ? (zbin) : (nZBins - 1 - zbin);
0596         assert(last_wavelet_hittime >= 0);
0597         assert(last_wavelet_hittime <= nZBins / 2);
0598
0599       }  //     if (last_layer != layer or last_phibin != phibin)
0600
0601       if (Verbosity() >= VERBOSITY_A_LOT)
0602       {
0603         cout << "TPCMLDataInterface::process_event -  layer " << layer << " hit with "
0604
0605              << "phibin = " << phibin
0606              << ",zbin = " << zbin
0607              << ",side = " << side
0608              << ",last_wavelet.size() = " << last_wavelet.size()
0609              << ",last_zbin = " << last_zbin
0610              << endl;
0611       }
0612
0613       // more checks on signal continuity
0614       if (last_wavelet.size() > 0)
0615       {
0616         //        if (side == 0)
0617         //        {
0618         assert(zbin - last_zbin == 1);
0619         //        }
0620         //        else
0621         //        {
0622         //          assert(last_zbin - zbin == 1);
0623         //        }
0624       }
0625
0626       // record adc
0627       unsigned int adc = hitr->second->getAdc();
0628       last_wavelet.push_back(adc);
0629       last_zbin = zbin;
0630
0631       // statistics
0632       layerChanHit[layer][side][phibin] += 1;
0633       assert(m_hLayerZBinHit);
0634       m_hLayerZBinHit->Fill(zbin, layer, 1);
0635       assert(m_hLayerZBinADC);
0636       m_hLayerZBinADC->Fill(zbin, layer, adc);
0637       assert(adc < (1 << 10));
0638       assert((hsize_t) phibin < layerDataBufferSize[layer][0]);
0639       assert((hsize_t) zbin < layerDataBufferSize[layer][1]);
0640       const size_t hitindex(layerDataBufferSize[layer][1] * phibin + zbin);
0641       assert((hsize_t) phibin < layerSignalBackgroundDataBufferSize[layer][0]);
0642       assert((hsize_t) zbin < layerSignalBackgroundDataBufferSize[layer][1]);
0643       const size_t hitindexSB(layerSignalBackgroundDataBufferSize[layer][1] * phibin + zbin);
0644       assert(hitindex < layerDataBuffer[layer].size());
0645       assert(hitindexSB < layerSignalBackgroundBuffer[layer].size());
0646
0647       if (layerDataBuffer[layer][hitindex] != 0)
0648       {
0649         cout << "TPCMLDataInterface::process_event - WARNING - hit @ layer "
0650              << layer << " phibin = " << phibin << " zbin = " << zbin << " side = " << side
0651              << " overwriting previous hit with ADC = " << layerDataBuffer[layer][hitindex]
0652              << endl;
0653       }
0654       layerDataBuffer[layer][hitindex] = adc;
0655
0656       if (layerSignalBackgroundBuffer[layer][hitindexSB] != 0)
0657       {
0658         cout << "TPCMLDataInterface::process_event - WARNING - signal/background @ layer "
0659              << layer << " phibin = " << phibin << " zbin = " << zbin << " side = " << side
0660              << " overwriting previous hit with = " << layerSignalBackgroundBuffer[layer][hitindexSB]
0661              << endl;
0662       }
0663
0664       // momentum cut
0665
0666       bool signal_or_bgd = false;
0667
0668       TrkrDefs::hitkey hit_key = hitr->first;
0669       //        PHG4Hit* g4hit = hiteval->max_truth_hit_by_energy(hit_key);
0670       std::set<PHG4Hit*> g4hit_set = hiteval->all_truth_hits(hit_key);
0671
0672       for (PHG4Hit* g4hit : g4hit_set)
0673       {
0674         PHG4Particle* g4particle = trutheval->get_particle(g4hit);
0675
0676         if (g4particle != nullptr)
0677         {
0678           gpx = g4particle->get_px();
0679           gpy = g4particle->get_py();
0680           gpz = g4particle->get_pz();
0681         }
0682
0683         if (sqrt(gpx * gpx + gpy * gpy + gpz * gpz) > m_momentumCut) signal_or_bgd = true;
0684
0685       }  //        for (auto g4hit: g4hit_set)
0686
0687       layerSignalBackgroundBuffer[layer][hitindexSB] = signal_or_bgd;
0688
0689     }  //     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0690
0691   }  //   for(SvtxHitMap::Iter iter = hits->begin(); iter != hits->end(); ++iter) {
0692
0693   // save last wavelet
0694   if (last_wavelet.size() > 0)
0695   {
0696     const int datasize = writeWavelet(last_layer, last_side, last_phibin, last_wavelet_hittime, last_wavelet);
0697     assert(datasize > 0);
0698
0699     nWavelet += 1;
0700     sumDataSize += datasize;
0701     layerChanDataSize[last_layer][last_side][last_phibin] += datasize;
0702   }
0703
0704   // statistics
0705   for (int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0706   {
0707     for (unsigned int side = 0; side < 2; ++side)
0708     {
0709       int sumHit = 0;
0710       for (const int& hit : layerChanHit[layer][side])
0711       {
0712         sumHit += hit;
0713
0714         assert(m_hLayerHit);
0715         m_hLayerHit->Fill(layer, hit);
0716         h_norm->Fill("TPC Hit", hit);
0717
0718         if (Verbosity() >= VERBOSITY_MORE)
0719         {
0720           cout << "TPCMLDataInterface::process_event - layerChanHit: "
0721                << "hit = " << hit
0722                << "at layer = " << layer
0723                << ", side = " << side
0724                << endl;
0725         }
0726       }
0727
0728       if (Verbosity() >= VERBOSITY_MORE)
0729       {
0730         cout << "TPCMLDataInterface::process_event - hLayerSumCellHit->Fill(" << layer << ", " << sumHit << ")" << endl;
0731       }
0732       assert(m_hLayerSumHit);
0733       m_hLayerSumHit->Fill(layer, sumHit);
0734
0735       double sumData = 0;
0736       for (const int& data : layerChanDataSize[layer][side])
0737       {
0738         sumData += data;
0739
0740         assert(m_hLayerDataSize);
0741         m_hLayerDataSize->Fill(layer, data);
0742       }
0743
0744       assert(m_hLayerSumDataSize);
0745       m_hLayerSumDataSize->Fill(layer, sumData);
0746       layerWaveletDataSize[(layer - m_minLayer)] += sumData;
0747     }  //    for (unsigned int side = 0; side < 2; ++side)
0748
0749     // store in H5
0750
0751     assert(layerH5DataSetMap[layer]);
0752     assert(layerH5SignalBackgroundMap[layer]);
0753     layerH5DataSetMap[layer]->write(layerDataBuffer[layer].data(), PredType::NATIVE_UINT16);
0754     layerH5SignalBackgroundMap[layer]->write(layerSignalBackgroundBuffer[layer].data(), PredType::NATIVE_UINT8);
0755     H5DataSetLayerWaveletDataSize->write(layerWaveletDataSize.data(), PredType::NATIVE_UINT32);
0756
0757   }  //  for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0758
0759   assert(m_hWavelet);
0760   m_hWavelet->Fill(nWavelet);
0761   h_norm->Fill("TPC Wavelet", nWavelet);
0762   assert(m_hDataSize);
0763   m_hDataSize->Fill(sumDataSize);
0764   h_norm->Fill("TPC DataSize", sumDataSize);
0765
0766   return Fun4AllReturnCodes::EVENT_OK;
0767 }
0768
0769 int TPCMLDataInterface::writeWavelet(int layer, int side, int phibin, int hittime, const vector<unsigned int>& wavelet)
0770 {
0771   static const int headersize = 2;  // 2-byte header per wavelet
0772
0773   //data in byte aligned and padding
0774   const int datasizebit = wavelet.size() * 10;
0775   int datasizebyte = datasizebit / 8;
0776   if (datasizebyte * 8 < datasizebit) datasizebyte += 1;
0777
0778   assert(m_hLayerWaveletSize);
0779   m_hLayerWaveletSize->Fill(layer, wavelet.size());
0780
0781   return headersize + datasizebyte;
0782 }
0783
0784 Fun4AllHistoManager*
0785 TPCMLDataInterface::getHistoManager()
0786 {
0787   static string histname("TPCMLDataInterface_HISTOS");
0788
0789   Fun4AllServer* se = Fun4AllServer::instance();
0790   Fun4AllHistoManager* hm = se->getHistoManager(histname);
0791
0792   if (not hm)
0793   {
0794     cout
0795         << "TPCMLDataInterface::get_HistoManager - Making Fun4AllHistoManager " << histname
0796         << endl;
0797     hm = new Fun4AllHistoManager(histname);
0798     se->registerHistoManager(hm);
0799   }
0800
0801   assert(hm);
0802
0803   return hm;
0804 }