packages/tpc2019/TpcPrototypeUnpacker.cc

0001 /*
0002  * TpcPrototypeUnpacker.cc
0003  *
0004  *  Created on: Sep 19, 2018
0005  *      Author: jinhuang
0006  */
0007
0008 #include "TpcPrototypeUnpacker.h"
0009
0010 #include "ChanMap.h"
0011 #include "TpcPrototypeCluster.h"
0012 #include "TpcPrototypeDefs.h"
0013
0014 #include <g4tpc/PHG4TpcPadPlane.h>  // for PHG4TpcPadPlane
0015
0016 #include <tpc/TpcDefs.h>
0017
0018 #include <g4detectors/PHG4CylinderCellGeom.h>
0019 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
0020
0021 #include <fun4all/Fun4AllBase.h>
0022 #include <fun4all/Fun4AllHistoManager.h>
0023 #include <fun4all/Fun4AllReturnCodes.h>
0024 #include <fun4all/Fun4AllServer.h>
0025 #include <fun4all/PHTFileServer.h>
0026 #include <fun4all/SubsysReco.h>
0027
0028 #include <phfield/PHFieldConfig.h>  // for PHFie...
0029 #include <phfield/PHFieldConfigv2.h>
0030 #include <phfield/PHFieldUtility.h>
0031
0032 #include <tpc/TpcDefs.h>
0033 #include <tpc/TpcHit.h>
0034 #include <trackbase/TrkrClusterContainer.h>
0035 #include <trackbase/TrkrDefs.h>  // for hitkey, getLayer
0036 #include <trackbase/TrkrHit.h>   // for TrkrHit
0037 #include <trackbase/TrkrHitSet.h>
0038 #include <trackbase/TrkrHitSetContainer.h>
0039
0040 #include <phool/PHCompositeNode.h>
0041 #include <phool/PHIODataNode.h>    // for PHIOD...
0042 #include <phool/PHNode.h>          // for PHNode
0043 #include <phool/PHNodeIterator.h>  // for PHNod...
0044 #include <phool/PHObject.h>        // for PHObject
0045 #include <phool/getClass.h>
0046 #include <phool/phool.h>  // for PHWHERE
0047
0048 #include <Event/Event.h>
0049 #include <Event/EventTypes.h>
0050 #include <Event/packet.h>
0051
0052 #include <TAxis.h>  // for TAxis
0053 #include <TClonesArray.h>
0054 #include <TH1.h>            // for TH1D
0055 #include <TMatrixFfwd.h>    // for TMatrixF
0056 #include <TMatrixT.h>       // for TMatrixT
0057 #include <TMatrixTUtils.h>  // for TMatr...
0058 #include <TNamed.h>         // for TNamed
0059 #include <TTree.h>
0060
0061 #include <boost/bimap.hpp>
0062 #include <boost/bind.hpp>
0063 #include <boost/format.hpp>
0064 #include <boost/graph/adjacency_list.hpp>
0065 #include <boost/graph/connected_components.hpp>
0066
0067 #include <algorithm>
0068 #include <cassert>
0069 #include <cmath>
0070 #include <cstdint>  // for uint32_t
0071 #include <iostream>
0072 #include <iterator>  // for rever...
0073 #include <map>
0074 #include <memory>
0075 #include <sstream>
0076 #include <stdexcept>
0077
0078 class PHField;
0079
0080 using namespace std;
0081 using namespace TpcPrototypeDefs::FEEv2;
0082
0083 TpcPrototypeUnpacker::TpcPrototypeUnpacker(const std::string& outputfilename)
0084   : SubsysReco("TpcPrototypeUnpacker")
0085   , padplane(nullptr)
0086   , tpcCylinderCellGeom(nullptr)
0087   , hitsetcontainer(nullptr)
0088   , trkrclusters(nullptr)
0089   , m_outputFileName(outputfilename)
0090   , m_eventT(nullptr)
0091   , m_peventHeader(&m_eventHeader)
0092   , m_nClusters(-1)
0093   , m_IOClusters(nullptr)
0094   , m_chanT(nullptr)
0095   , m_pchanHeader(&m_chanHeader)
0096   , m_chanData(kSAMPLE_LENGTH, 0)
0097   , enableClustering(true)
0098   , m_clusteringZeroSuppression(50)
0099   , m_nPreSample(5)
0100   , m_nPostSample(5)
0101   , m_XRayLocationX(-1)
0102   , m_XRayLocationY(-1)
0103   , m_pdfMaker(nullptr)
0104 {
0105 }
0106
0107 TpcPrototypeUnpacker::~TpcPrototypeUnpacker()
0108 {
0109   if (m_IOClusters)
0110   {
0111     m_IOClusters->Clear();
0112     delete m_IOClusters;
0113   }
0114   if (m_pdfMaker)
0115   {
0116     delete m_pdfMaker;
0117   }
0118   if (padplane)
0119   {
0120     delete padplane;
0121   }
0122 }
0123
0124 int TpcPrototypeUnpacker::ResetEvent(PHCompositeNode* topNode)
0125 {
0126   m_eventHeader = EventHeader();
0127   m_padPlaneData.Reset();
0128   m_clusters.clear();
0129   m_chanHeader = ChannelHeader();
0130
0131   m_nClusters = -1;
0132   assert(m_IOClusters);
0133   m_IOClusters->Clear("C");
0134
0135   return Fun4AllReturnCodes::EVENT_OK;
0136 }
0137
0138 void TpcPrototypeUnpacker::ClusterData::Clear(Option_t*)
0139 {
0140   pad_radials.clear();
0141   pad_azimuths.clear();
0142   samples.clear();
0143
0144   //  for (auto& s : pad_radial_samples) s.second.clear();
0145   //  pad_radial_samples.clear();
0146   //
0147   //  for (auto& s : pad_azimuth_samples) s.second.clear();
0148   //  pad_azimuth_samples.clear();
0149
0150   while (pad_radial_samples.begin() != pad_radial_samples.end())
0151   {
0152     pad_radial_samples.begin()->second.clear();
0153     pad_radial_samples.erase(pad_radial_samples.begin());
0154   }
0155
0156   while (pad_azimuth_samples.begin() != pad_azimuth_samples.end())
0157   {
0158     pad_azimuth_samples.begin()->second.clear();
0159     pad_azimuth_samples.erase(pad_azimuth_samples.begin());
0160   }
0161
0162   while (pad_azimuth_peaks.begin() != pad_azimuth_peaks.end())
0163   {
0164     pad_azimuth_peaks.erase(pad_azimuth_peaks.begin());
0165   }
0166
0167   //  while (sum_samples.begin() != sum_samples.end())
0168   //  {
0169   //    sum_samples.erase(sum_samples.begin());
0170   //  }
0171   sum_samples.clear();
0172   sum_samples.shrink_to_fit();
0173 }
0174
0175 int TpcPrototypeUnpacker::Init(PHCompositeNode* topNode)
0176 {
0177   if (padplane)
0178     padplane->Init(topNode);
0179
0180   return Fun4AllReturnCodes::EVENT_OK;
0181 }
0182
0183 int TpcPrototypeUnpacker::InitRun(PHCompositeNode* topNode)
0184 {
0185   // Looking for the DST node
0186   PHNodeIterator iter(topNode);
0187   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0188   if (!dstNode)
0189   {
0190     cout << PHWHERE << "DST Node missing, doing nothing." << endl;
0191     return Fun4AllReturnCodes::ABORTRUN;
0192   }
0193
0194   if (padplane)
0195   {
0196     if (Verbosity() >= VERBOSITY_SOME)
0197     {
0198       cout << "TpcPrototypeUnpacker::InitRun - making pad plane"
0199            << endl;
0200     }
0201
0202     //setup the constant field
0203     const int field_ret = InitField(topNode);
0204     if (field_ret != Fun4AllReturnCodes::EVENT_OK)
0205     {
0206       cout << "TpcPrototypeUnpacker::InitRun- Error - Failed field init with status = " << field_ret << endl;
0207       return field_ret;
0208     }
0209
0210     string seggeonodename = "CYLINDERCELLGEOM_SVTX";  // + detector;
0211     tpcCylinderCellGeom = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, seggeonodename.c_str());
0212     if (!tpcCylinderCellGeom)
0213     {
0214       tpcCylinderCellGeom = new PHG4CylinderCellGeomContainer();
0215       PHNodeIterator iter(topNode);
0216       PHCompositeNode* runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
0217       PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(tpcCylinderCellGeom, seggeonodename.c_str(), "PHObject");
0218       runNode->addNode(newNode);
0219     }
0220
0221     padplane->InitRun(topNode);
0222     padplane->CreateReadoutGeometry(topNode, tpcCylinderCellGeom);
0223   }
0224
0225   // Create the TrkrHit node if required
0226
0227   hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0228   if (!hitsetcontainer)
0229   {
0230     PHNodeIterator dstiter(dstNode);
0231     PHCompositeNode* DetNode =
0232         dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0233     if (!DetNode)
0234     {
0235       DetNode = new PHCompositeNode("TRKR");
0236       dstNode->addNode(DetNode);
0237     }
0238
0239     hitsetcontainer = new TrkrHitSetContainer();
0240     PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(hitsetcontainer, "TRKR_HITSET", "PHObject");
0241     DetNode->addNode(newNode);
0242   }
0243
0244   // Create the Cluster node if required
0245   trkrclusters = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0246   if (!trkrclusters)
0247   {
0248
0249     PHNodeIterator dstiter(dstNode);
0250     PHCompositeNode* DetNode =
0251         dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0252     if (!DetNode)
0253     {
0254       DetNode = new PHCompositeNode("TRKR");
0255       dstNode->addNode(DetNode);
0256     }
0257
0258     trkrclusters = new TrkrClusterContainer();
0259     PHIODataNode<PHObject>* TrkrClusterContainerNode =
0260         new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
0261     DetNode->addNode(TrkrClusterContainerNode);
0262   }
0263
0264   if (Verbosity() >= VERBOSITY_SOME)
0265   {
0266     cout << "TpcPrototypeUnpacker::InitRun - Making PHTFileServer " << m_outputFileName
0267          << endl;
0268
0269     m_pdfMaker = new SampleFit_PowerLawDoubleExp_PDFMaker();
0270   }
0271   PHTFileServer::get().open(m_outputFileName, "RECREATE");
0272
0273   Fun4AllHistoManager* hm = getHistoManager();
0274   assert(hm);
0275
0276   TH1D* h = new TH1D("hNormalization",  //
0277                      "Normalization;Items;Summed quantity", 10, .5, 10.5);
0278   int i = 1;
0279   h->GetXaxis()->SetBinLabel(i++, "Event count");
0280   h->GetXaxis()->SetBinLabel(i++, "Collision count");
0281   h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit");
0282   h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit Edep");
0283   h->GetXaxis()->SetBinLabel(i++, "TPC Pad Hit");
0284   h->GetXaxis()->SetBinLabel(i++, "TPC Charge e");
0285   h->GetXaxis()->SetBinLabel(i++, "TPC Charge fC");
0286   h->GetXaxis()->LabelsOption("v");
0287   hm->registerHisto(h);
0288
0289   m_eventT = new TTree("eventT", "TPC FEE per-event Tree");
0290   assert(m_eventT);
0291   m_eventT->Branch("evthdr", &m_peventHeader);
0292   m_eventT->Branch("nClusters", &m_nClusters, "nClusters/I");
0293   m_IOClusters = new TClonesArray("TpcPrototypeUnpacker::ClusterData", 1000);
0294   m_eventT->Branch("Clusters", &m_IOClusters);
0295
0296   m_chanT = new TTree("chanT", "TPC FEE per-channel Tree");
0297   assert(m_chanT);
0298   m_chanT->Branch("event", &m_eventHeader.event, "event/I");
0299   m_chanT->Branch("chanhdr", &m_pchanHeader);
0300   m_chanT->Branch("adc", m_chanData.data(), str(boost::format("adc[%d]/i") % kSAMPLE_LENGTH).c_str());
0301
0302   //  for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
0303   //  {
0304   //    const PHG4CylinderCellGeom* layer_geom = seggeo->GetLayerCellGeom(layer);
0305
0306   //    const string histNameCellHit(boost::str(boost::format{"hCellHit_Layer%1%"} % layer));
0307   //    const string histNameCellCharge(boost::str(boost::format{"hCellCharge_Layer%1%"} % layer));
0308
0309   //  }
0310
0311   //  hm->registerHisto(new TH2D("hLayerCellHit",  //
0312   //                             "Number of ADC time-bin hit per channel;Layer ID;Hit number",
0313   //                             m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0314   //                             300, -.5, 299.5));
0315   //  hm->registerHisto(new TH2D("hLayerCellCharge",  //
0316   //                             "Charge integrated over drift window per channel;Layer ID;Charge [fC]",
0317   //                             m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0318   //                             1000, 0, 1e7 * eplus / (1e-15 * coulomb)));
0319   //
0320   //  hm->registerHisto(new TH2D("hLayerSumCellHit",  //
0321   //                             "Number of ADC time-bin hit integrated over channels per layer;Layer ID;Hit number",
0322   //                             m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0323   //                             10000, -.5, 99999.5));
0324   //  hm->registerHisto(new TH2D("hLayerSumCellCharge",  //
0325   //                             "Charge integrated over drift window and channel per layer;Layer ID;Charge [fC]",
0326   //                             m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
0327   //                             10000, 0, 1000 * 4e6 * eplus / (1e-15 * coulomb)));
0328
0329   return Fun4AllReturnCodes::EVENT_OK;
0330 }
0331
0332 int TpcPrototypeUnpacker::End(PHCompositeNode* topNode)
0333 {
0334   if (Verbosity() >= VERBOSITY_SOME)
0335   {
0336     cout << "TpcPrototypeUnpacker::End - write to " << m_outputFileName << endl;
0337   }
0338   PHTFileServer::get().cd(m_outputFileName);
0339
0340   Fun4AllHistoManager* hm = getHistoManager();
0341   assert(hm);
0342   for (unsigned int i = 0; i < hm->nHistos(); i++)
0343     hm->getHisto(i)->Write();
0344
0345   // help index files with TChain
0346   TTree* T_Index = new TTree("T_Index", "T_Index");
0347   assert(T_Index);
0348   T_Index->Write();
0349
0350   m_eventT->Write();
0351   m_chanT->Write();
0352
0353   if (m_pdfMaker)
0354   {
0355     delete m_pdfMaker;
0356     m_pdfMaker = nullptr;
0357   }
0358   return Fun4AllReturnCodes::EVENT_OK;
0359 }
0360
0361 int TpcPrototypeUnpacker::process_event(PHCompositeNode* topNode)
0362 {
0363   Fun4AllHistoManager* hm = getHistoManager();
0364   assert(hm);
0365   TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto("hNormalization"));
0366   assert(h_norm);
0367
0368   Event* event = findNode::getClass<Event>(topNode, "PRDF");
0369   if (event == nullptr)
0370   {
0371     if (Verbosity() >= VERBOSITY_SOME)
0372       cout << "TpcPrototypeUnpacker::Process_Event - Event not found" << endl;
0373     return Fun4AllReturnCodes::DISCARDEVENT;
0374   }
0375
0376   if (Verbosity() >= VERBOSITY_SOME)
0377     event->identify();
0378
0379   // search for data event
0380   if (event->getEvtType() == BEGRUNEVENT)
0381   {
0382     //    get_motor_loc(event);
0383
0384     return Fun4AllReturnCodes::EVENT_OK;
0385   }
0386   if (event->getEvtType() != DATAEVENT)
0387     return Fun4AllReturnCodes::DISCARDEVENT;
0388
0389   m_eventHeader.run = event->getRunNumber();
0390   m_eventHeader.event = event->getEvtSequence();
0391
0392   //  m_eventHeader.xray_x = m_XRayLocationX;
0393   //  m_eventHeader.xray_y = m_XRayLocationY;
0394
0395   if (m_pdfMaker)
0396   {
0397     m_pdfMaker->MakeSectionPage(str(boost::format("ADC signal fits for Run %1% and event %2%") % m_eventHeader.run % m_eventHeader.event));
0398   }
0399
0400   static const char* MAX_LENGTH = "MAX_SAMPLES";
0401   static const char* IS_PRESENT = "IS_PRESENT";
0402   static const char* BX_COUNTER = "BX";
0403
0404   unique_ptr<Packet> p(event->getPacket(kPACKET_ID));
0405   if (p == nullptr)
0406     return Fun4AllReturnCodes::DISCARDEVENT;
0407
0408   if (Verbosity() >= VERBOSITY_SOME) p->identify();
0409
0410   if (Verbosity() >= VERBOSITY_EVEN_MORE)
0411   {
0412     cout << "TpcPrototypeUnpacker::process_event - package dump" << endl;
0413     p->dump();
0414   }
0415
0416   int max_length[kN_FEES] = {-1};
0417   for (unsigned int i = 0; i < kN_FEES; i++)
0418   {
0419     try
0420     {
0421       max_length[i] = p->iValue(i, MAX_LENGTH);
0422     }
0423     catch (const std::out_of_range& e)
0424     {
0425       max_length[i] = -1;
0426       break;
0427     }
0428
0429     if (Verbosity() >= VERBOSITY_MORE)
0430     {
0431       cout << "TpcPrototypeUnpacker::process_event - max_length[" << i << "]=" << max_length[i] << endl;
0432     }
0433   }
0434
0435   m_eventHeader.bx_counter = 0;
0436   //  bool first_channel = true;
0437
0438   for (unsigned int fee = 0; fee < kN_FEES; ++fee)
0439   {
0440     try
0441     {
0442       if (!p->iValue(fee, IS_PRESENT))
0443       {
0444         if (Verbosity() >= VERBOSITY_MORE)
0445         {
0446           cout << "TpcPrototypeUnpacker::process_event - missing fee[" << fee << "]=" << endl;
0447         }
0448         continue;
0449       }
0450     }
0451     catch (const std::out_of_range& e)
0452     {
0453       cout << "TpcPrototypeUnpacker::process_event - out_of_range in p->iValue(" << fee << ", IS_PRESENT)" << endl;
0454
0455       break;
0456     }
0457
0458     if (Verbosity() >= VERBOSITY_MORE)
0459     {
0460       cout << "TpcPrototypeUnpacker::process_event - processing fee[" << fee << "]=" << endl;
0461     }
0462
0463     for (unsigned int channel = 0; channel < kN_CHANNELS; channel++)
0464     {
0465       m_chanHeader = ChannelHeader();
0466       m_chanHeader.fee_id = fee;
0467       m_chanHeader.size = p->iValue(fee, channel, "NR_SAMPLES");  // number of words until the next channel (header included). this is the real packet_length
0468
0469       if (Verbosity() >= VERBOSITY_MORE)
0470       {
0471         cout << "TpcPrototypeUnpacker::process_event - processing fee["
0472              << fee << "], chan[" << channel << "] NR_SAMPLES = "
0473              << p->iValue(fee, channel, "NR_SAMPLES") << endl;
0474       }
0475
0476       unsigned int real_t = 0;
0477       uint32_t old_bx_count = 0;
0478       uint16_t adc = 0;
0479       uint32_t bx_count = 0;
0480       m_chanData.resize(kSAMPLE_LENGTH, 0);
0481
0482       for (unsigned int t = 0; t < kSAMPLE_LENGTH; t++)
0483       {
0484         try
0485         {
0486           adc = p->iValue(fee, channel, t);
0487           bx_count = p->iValue(fee, channel, t, BX_COUNTER);
0488
0489           if (real_t == 0)
0490             m_chanHeader.bx_counter = p->iValue(fee, channel, t, BX_COUNTER);
0491         }
0492         catch (const std::out_of_range& e)
0493         {
0494           adc = 0;
0495           bx_count = 0;
0496           continue;
0497         }
0498
0499         if (bx_count >= old_bx_count + 128)
0500         {
0501           old_bx_count = bx_count;
0502           real_t = 0;
0503           m_chanData.resize(kSAMPLE_LENGTH, 0);
0504         }
0505
0506         assert(real_t < kSAMPLE_LENGTH);
0507         m_chanData[real_t] = adc;
0508         //          h_raw_wave->Fill(real_t, adc);
0509         //          h_raw->Fill(real_t, total_chan, adc);
0510         //          adc_data->push_back(adc);
0511         real_t++;
0512       }  //      for (unsigned int t = 3; t < kSAMPLE_LENGTH; t++)
0513
0514       //      m_chanHeader.packet_type = p->iValue(channel * kPACKET_LENGTH + 2) & 0xffff;  // that's the Elink packet type
0515       //      m_chanHeader.bx_counter = ((p->iValue(channel * kPACKET_LENGTH + 4) & 0xffff) << 4) | (p->iValue(channel * kPACKET_LENGTH + 5) & 0xffff);
0516       //      m_chanHeader.sampa_address = (p->iValue(channel * kPACKET_LENGTH + 3) >> 5) & 0xf;
0517       //      m_chanHeader.sampa_channel = p->iValue(channel * kPACKET_LENGTH + 3) & 0x1f;
0518       m_chanHeader.fee_channel = channel;
0519
0520       //            const pair<int, int> pad = SAMPAChan2PadXY(m_chanHeader.fee_channel);
0521
0522       m_chanHeader.pad_x = TpcR2Map.GetRpos(fee, channel);
0523       m_chanHeader.pad_y = TpcR2Map.Getphipos(fee, channel);
0524
0525       //
0526       if (Verbosity() >= VERBOSITY_MORE)
0527       {
0528         cout << "TpcPrototypeUnpacker::process_event - "
0529              << "m_chanHeader.m_size = " << int(m_chanHeader.size) << ", "
0530              << "m_chanHeader.m_packet_type = " << int(m_chanHeader.packet_type) << ", "
0531              << "m_chanHeader.m_bx_counter = " << int(m_chanHeader.bx_counter) << ", "
0532              << "m_chanHeader.m_sampa_address = " << int(m_chanHeader.sampa_address) << ", "
0533              << "m_chanHeader.m_sampa_channel = " << int(m_chanHeader.sampa_channel) << ", "
0534              << "m_chanHeader.m_fee_channel = " << int(m_chanHeader.fee_channel) << ": "
0535              << " ";
0536
0537         for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
0538         {
0539           cout << "data[" << sample << "] = " << int(m_chanData[sample]) << " ";
0540         }
0541
0542         cout << endl;
0543       }
0544
0545       // fill event data
0546       //      if (PadPlaneData::IsValidPad(m_chanHeader.pad_x, m_chanHeader.pad_y))
0547       //      {
0548       vector<int>& paddata = m_padPlaneData.getPad(m_chanHeader.pad_x, m_chanHeader.pad_y);
0549       //
0550       for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
0551       {
0552         paddata[sample] = int(m_chanData[sample]);
0553       }
0554       //
0555       auto pedestal_max = roughZeroSuppression(paddata);
0556       m_chanHeader.pedestal = pedestal_max.first;
0557       m_chanHeader.max = pedestal_max.second;
0558       //      }
0559       // output per-channel TTree
0560       m_chanT->Fill();
0561     }  //    for (unsigned int channel = 0; channel < kN_CHANNELS; channel++)
0562   }    //  for (unsigned int fee = 0; fee < kN_FEES; ++fee)
0563
0564   exportDSTHits();
0565
0566   if (enableClustering)
0567   {
0568     int ret = Clustering();
0569
0570     if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
0571   }
0572   h_norm->Fill("Event count", 1);
0573   m_eventT->Fill();
0574
0575   return Fun4AllReturnCodes::EVENT_OK;
0576 }
0577
0578 int TpcPrototypeUnpacker::exportDSTHits()
0579 {
0580   int nHits(0);
0581
0582   assert(hitsetcontainer);
0583
0584   for (unsigned int pad_radial = 0; pad_radial < kMaxPadX; ++pad_radial)
0585   {
0586     TrkrDefs::hitsetkey hitsetkey = TpcDefs::genHitSetKey(pad_radial, 0, 0);
0587     // Use existing hitset or add new one if needed
0588     TrkrHitSetContainer::Iterator hitsetit = hitsetcontainer->findOrAddHitSet(hitsetkey);
0589
0590     for (unsigned int pad_azimuth = 0; pad_azimuth < kMaxPadY; ++pad_azimuth)
0591     {
0592       for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
0593       {
0594         unsigned int adc = static_cast<unsigned int>(m_padPlaneData.getData()[pad_azimuth][pad_radial][sample]);
0595
0596         if (adc)
0597         {
0598           // generate the key for this hit, requires zbin and phibin
0599           TrkrDefs::hitkey hitkey = TpcDefs::genHitKey(pad_azimuth, sample);
0600           // See if this hit already exists
0601           TrkrHit* hit = nullptr;
0602           hit = hitsetit->second->getHit(hitkey);
0603           if (!hit)
0604           {
0605             // create a new one
0606             hit = new TpcHit();
0607             hitsetit->second->addHitSpecificKey(hitkey, hit);
0608           }
0609           else
0610           {
0611             cout << "TpcPrototypeUnpacker::exportDSTHits "
0612                  << " - Fatal error - hit with "
0613                  << " hitkey = " << hitkey << " already exist";
0614             exit(1);
0615           }
0616
0617           // Either way, add the energy to it  -- adc values will be added at digitization
0618           hit->setAdc(adc);
0619           ++nHits;
0620
0621         }  //        if (m_padPlaneData.m_data[pad_azimuth][pad_radial][sample])
0622
0623       }  //      for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
0624     }    //   for (unsigned int pad_azimuth = 0; pad_azimuth < kMaxPadY; ++pad_azimuth)
0625   }      //  for (unsigned int pad_radial = 0; pad_radial < kMaxPadX; ++pad_radial)
0626
0627   return nHits;
0628 }
0629
0630 int TpcPrototypeUnpacker::Clustering()
0631 {
0632   if (Verbosity())
0633     cout << __PRETTY_FUNCTION__ << " entry" << endl;
0634
0635   // find cluster
0636   m_padPlaneData.Clustering(m_clusteringZeroSuppression, Verbosity() >= VERBOSITY_SOME);
0637   const multimap<int, PadPlaneData::SampleID>& groups = m_padPlaneData.getGroups();
0638
0639   // export clusters
0640   assert(m_clusters.size() == 0);  //already cleared.
0641   for (const auto& iter : groups)
0642   {
0643     const int& i = iter.first;
0644     const PadPlaneData::SampleID& id = iter.second;
0645     m_clusters[i].pad_radials.insert(id.pad_radial);
0646     m_clusters[i].pad_azimuths.insert(id.pad_azimuth);
0647     m_clusters[i].samples.insert(id.sample);
0648   }
0649
0650   // process cluster
0651   for (auto& iter : m_clusters)
0652   {
0653     ClusterData& cluster = iter.second;
0654
0655     assert(cluster.pad_radials.size() > 0);
0656     assert(cluster.pad_azimuths.size() > 0);
0657     assert(cluster.samples.size() > 0);
0658
0659     //expand cluster by +/-1 in y
0660     if (*(cluster.pad_azimuths.begin()) - 1 >= 0)
0661       cluster.pad_azimuths.insert(*(cluster.pad_azimuths.begin()) - 1);
0662     if (*(cluster.pad_azimuths.rbegin()) + 1 < (int) kMaxPadY)
0663       cluster.pad_azimuths.insert(*(cluster.pad_azimuths.rbegin()) + 1);
0664
0665     cluster.min_sample = max(0, *cluster.samples.begin() - m_nPreSample);
0666     cluster.max_sample = min((int) (kSAMPLE_LENGTH) -1, *cluster.samples.rbegin() + m_nPostSample);
0667     const int n_sample = cluster.max_sample - cluster.min_sample + 1;
0668
0669     cluster.sum_samples.assign(n_sample, 0);
0670     for (int pad_x = *cluster.pad_radials.begin(); pad_x <= *cluster.pad_radials.rbegin(); ++pad_x)
0671     {
0672       cluster.pad_radial_samples[pad_x].assign(n_sample, 0);
0673     }
0674     for (int pad_y = *cluster.pad_azimuths.begin(); pad_y <= *cluster.pad_azimuths.rbegin(); ++pad_y)
0675     {
0676       cluster.pad_azimuth_samples[pad_y].assign(n_sample, 0);
0677     }
0678
0679     for (int pad_x = *cluster.pad_radials.begin(); pad_x <= *cluster.pad_radials.rbegin(); ++pad_x)
0680     {
0681       for (int pad_y = *cluster.pad_azimuths.begin(); pad_y <= *cluster.pad_azimuths.rbegin(); ++pad_y)
0682       {
0683         assert(m_padPlaneData.IsValidPad(pad_x, pad_y));
0684
0685         vector<int>& padsamples = m_padPlaneData.getPad(pad_x, pad_y);
0686
0687         for (int i = 0; i < n_sample; ++i)
0688         {
0689           int adc = padsamples.at(cluster.min_sample + i);
0690           cluster.sum_samples[i] += adc;
0691           cluster.pad_radial_samples[pad_x][i] += adc;
0692           cluster.pad_azimuth_samples[pad_y][i] += adc;
0693         }
0694
0695       }  //         for (int pad_y = *cluster.pad_azimuths.begin(); pad_y<=*cluster.pad_azimuths.rbegin() ;++pad_azimuth)
0696
0697     }  //    for (int pad_x = *cluster.pad_radials.begin(); pad_x<=*cluster.pad_radials.rbegin() ;++pad_radial)
0698
0699     if (m_pdfMaker)
0700     {
0701       m_pdfMaker->MakeSectionPage(str(boost::format("Event %1% Cluster %2%: sum all channel fit followed by fit of each pad component") % m_eventHeader.event % iter.first));
0702     }
0703
0704     // fit - overal cluster
0705     map<int, double> parameters_constraints;
0706     {
0707       double peak = NAN;
0708       double peak_sample = NAN;
0709       double pedstal = NAN;
0710       map<int, double> parameters_io;
0711       SampleFit_PowerLawDoubleExp(cluster.sum_samples, peak,
0712                                   peak_sample, pedstal, parameters_io, Verbosity());
0713
0714       parameters_constraints[1] = parameters_io[1];
0715       parameters_constraints[2] = parameters_io[2];
0716       parameters_constraints[3] = parameters_io[3];
0717       parameters_constraints[5] = parameters_io[5];
0718       parameters_constraints[6] = parameters_io[6];
0719
0720       cluster.peak = peak;
0721       cluster.peak_sample = peak_sample;
0722       cluster.pedstal = pedstal;
0723     }
0724
0725     // fit - X -> radial direction
0726     {
0727       //      double sum_peak = 0;
0728       //      double sum_peak_pad_radial = 0;
0729       //      for (int pad_x = *cluster.pad_radials.begin(); pad_x <= *cluster.pad_radials.rbegin(); ++pad_x)
0730       //      {
0731       //        double peak = NAN;
0732       //        double peak_sample = NAN;
0733       //        double pedstal = NAN;
0734       //        map<int, double> parameters_io(parameters_constraints);
0735       //
0736       //        SampleFit_PowerLawDoubleExp(cluster.pad_radial_samples[pad_x], peak,
0737       //                                    peak_sample, pedstal, parameters_io, Verbosity());
0738       //
0739       //        cluster.pad_radial_peaks[pad_x] = peak;
0740       //        sum_peak += peak;
0741       //        sum_peak_pad_radial += peak * pad_x;
0742       //      }
0743       //      cluster.avg_pad_radial = sum_peak_pad_radial / sum_peak;
0744       //      cluster.size_pad_radial = cluster.pad_radials.size();
0745
0746       assert(cluster.pad_radials.size() == 1);
0747       cluster.avg_pad_radial = *cluster.pad_radials.begin();
0748       cluster.size_pad_radial = cluster.pad_radials.size();
0749     }
0750
0751     // fit - Y -> azimuthal direction
0752     {
0753       double sum_peak = 0;
0754       double sum_peak_pad_azimuth = 0;
0755       for (int pad_y = *cluster.pad_azimuths.begin(); pad_y <= *cluster.pad_azimuths.rbegin(); ++pad_y)
0756       {
0757         double peak = NAN;
0758         double peak_sample = NAN;
0759         double pedstal = NAN;
0760         map<int, double> parameters_io(parameters_constraints);
0761
0762         SampleFit_PowerLawDoubleExp(cluster.pad_azimuth_samples[pad_y], peak,
0763                                     peak_sample, pedstal, parameters_io, Verbosity());
0764
0765         cluster.pad_azimuth_peaks[pad_y] = peak;
0766         sum_peak += peak;
0767         sum_peak_pad_azimuth += peak * pad_y;
0768       }
0769       cluster.avg_pad_azimuth = sum_peak_pad_azimuth / sum_peak;
0770       cluster.size_pad_azimuth = cluster.pad_azimuths.size();
0771
0772       cluster.min_pad_azimuth = *cluster.pad_azimuths.begin();
0773       cluster.max_pad_azimuth = *cluster.pad_azimuths.rbegin();
0774     }
0775   }  //   for (auto& iter : m_clusters)
0776
0777   // sort by energy
0778   map<double, int> cluster_energy;
0779   for (auto& iter : m_clusters)
0780   {
0781     //reverse energy sorting
0782     cluster_energy[-iter.second.peak] = iter.first;
0783   }
0784
0785   // save clusters
0786   m_nClusters = 0;
0787   assert(m_IOClusters);
0788   for (const auto& iter : cluster_energy)
0789   {
0790     ClusterData& cluster = m_clusters[iter.second];
0791
0792     //output to DST clusters
0793     cluster.clusterID = m_nClusters;  // sync cluster id from cluster container to m_nClusters::ClusterData
0794     int ret = exportDSTCluster(cluster, m_nClusters);
0795     if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
0796
0797     // super awkward ways of ROOT filling TClonesArray
0798     new ((*m_IOClusters)[m_nClusters++]) ClusterData(cluster);
0799   }
0800   return Fun4AllReturnCodes::EVENT_OK;
0801 }
0802
0803 int TpcPrototypeUnpacker::exportDSTCluster(ClusterData& cluster, const int iclus)
0804 {
0805   assert(tpcCylinderCellGeom);
0806   assert(trkrclusters);
0807
0808   assert(cluster.avg_pad_radial >= 0);
0809   assert(cluster.avg_pad_radial < (int) kMaxPadX);
0810   const uint8_t layer = static_cast<uint8_t>(cluster.avg_pad_radial);
0811   const PHG4CylinderCellGeom* layergeom = tpcCylinderCellGeom->GetLayerCellGeom(layer);
0812   const int NPhiBins = layergeom->get_phibins();
0813   const int NZBins = layergeom->get_zbins();
0814
0815   // create the cluster entry directly in the node tree
0816   const TrkrDefs::hitsetkey hit_set_key(TpcDefs::genHitSetKey(layer, 0, 0));
0817   const TrkrDefs::cluskey ckey = TpcDefs::genClusKey(hit_set_key, iclus);
0818
0819   // make sure this cluster location is available
0820   if (trkrclusters->findCluster(ckey))
0821   {
0822     cout << "TpcPrototypeUnpacker::exportDSTCluster - fatal error - cluster already exist which is not expected for prototype analysis. "
0823          << "iclus = " << iclus << ", ckey = " << ckey << ". Offending cluster: " << endl;
0824     trkrclusters->findCluster(ckey)->identify();
0825     exit(1);
0826   }
0827   TpcPrototypeCluster* clus = new TpcPrototypeCluster();
0828   assert(clus);
0829   clus->setClusKey(ckey);
0830   trkrclusters->addCluster(clus);
0831
0832   //  calculate geometry
0833   const double radius = layergeom->get_radius();  // returns center of layer
0834   if (cluster.avg_pad_azimuth < 0)
0835   {
0836     cout << __PRETTY_FUNCTION__ << " WARNING - cluster.avg_pad_azimuth = " << cluster.avg_pad_azimuth << " < 0"
0837          << ", cluster.avg_pad_radial = " << cluster.avg_pad_radial << endl;
0838
0839     return Fun4AllReturnCodes::ABORTEVENT;
0840   }
0841   if (cluster.avg_pad_azimuth >= NPhiBins)
0842   {
0843     cout << __PRETTY_FUNCTION__ << " WARNING - cluster.avg_pad_azimuth = " << cluster.avg_pad_azimuth << " > " << NPhiBins
0844          << ", cluster.avg_pad_radial = " << cluster.avg_pad_radial << endl;
0845
0846     return Fun4AllReturnCodes::ABORTEVENT;
0847   }
0848   const int lowery = floor(cluster.avg_pad_azimuth);
0849
0850   if (lowery < 0)
0851   {
0852     cout << __PRETTY_FUNCTION__ << " WARNING - cluster.avg_pad_azimuth = " << cluster.avg_pad_azimuth << " -> "
0853          << " lower = " << lowery << " < 0"
0854          << ", cluster.avg_pad_radial = " << cluster.avg_pad_radial << endl;
0855
0856     return Fun4AllReturnCodes::ABORTEVENT;
0857   }
0858
0859   const double clusphi = layergeom->get_phicenter(lowery)                                             //
0860                          + (layergeom->get_phicenter(lowery + 1) - layergeom->get_phicenter(lowery))  //
0861                                * (cluster.avg_pad_azimuth - lowery);
0862
0863   assert(cluster.min_sample >= 0);
0864   assert(cluster.min_sample + cluster.peak_sample < NZBins);
0865   const double clusz = layergeom->get_zcenter(cluster.min_sample)  //
0866                        + (layergeom->get_zcenter(cluster.min_sample + 1) - layergeom->get_zcenter(cluster.min_sample)) * cluster.peak_sample;
0867
0868   const double phi_size = cluster.size_pad_azimuth;                 // * radius * layergeom->get_phistep();
0869   const double z_size = (cluster.max_sample - cluster.min_sample);  // * layergeom->get_zstep();
0870
0871   static const double phi_err = 170e-4;
0872
0873   static const double z_err = 1000e-4;
0874
0875   // Fill in the cluster details
0876   //================
0877   clus->setAdc(cluster.peak);
0878   clus->setPosition(0, radius * cos(clusphi));
0879   clus->setPosition(1, radius * sin(clusphi));
0880   clus->setPosition(2, clusz);
0881   clus->setGlobal();
0882
0883   //update cluster positions
0884   cluster.avg_pos_x = clus->getPosition(0);
0885   cluster.avg_pos_y = clus->getPosition(1);
0886   cluster.avg_pos_z = clus->getPosition(2);
0887
0888   cluster.delta_z = (layergeom->get_zcenter(cluster.min_sample + 1) - layergeom->get_zcenter(cluster.min_sample));
0889   cluster.delta_azimuth_bin = (layergeom->get_phicenter(lowery + 1) - layergeom->get_phicenter(lowery));
0890
0891   // output prototype specifics
0892
0893   //  std::set<int> pad_radials;
0894   clus->setPadRadials(cluster.pad_radials);
0895   //  std::set<int> pad_azimuths;
0896   clus->setPadAzimuths(cluster.pad_azimuths);
0897   //  std::set<int> samples;
0898   clus->setSamples(cluster.samples);
0899   //
0900   //  std::map<int, std::vector<double>> pad_radial_samples;
0901   clus->setPadRadialSamples(cluster.pad_radial_samples);
0902   //  std::map<int, std::vector<double>> pad_azimuth_samples;
0903   clus->setPadAzimuthSamples(cluster.pad_azimuth_samples);
0904   //  std::vector<double> sum_samples;
0905   clus->setSumSamples(cluster.sum_samples);
0906   //
0907   //  int min_sample;
0908   clus->setMinSample(cluster.min_sample);
0909   //  int max_sample;
0910   clus->setMaxSample(cluster.max_sample);
0911   //  int min_pad_azimuth;
0912   clus->setMinPadAzimuth(cluster.min_pad_azimuth);
0913   //  int max_pad_azimuth;
0914   clus->setMaxPadAzimuth(cluster.max_pad_azimuth);
0915   //
0916   //  double peak;
0917   clus->setPeak(cluster.peak);
0918   //  double peak_sample;
0919   clus->setPeakSample(cluster.peak_sample);
0920   //  double pedstal;
0921   clus->setPedstal(cluster.pedstal);
0922   //
0923   //  std::map<int, double> pad_azimuth_peaks;
0924   clus->setPadAzimuthPeaks(cluster.pad_azimuth_peaks);
0925   //
0926   //  //! pad coordinate
0927   //  int avg_pad_radial;
0928   clus->setAvgPadRadial(cluster.avg_pad_radial);
0929   //  double avg_pad_azimuth;
0930   clus->setAvgPadAzimuth(cluster.avg_pad_azimuth);
0931   //
0932   //  //! cluster size in units of pad bins
0933   //  int size_pad_radial;
0934   clus->setSizePadRadial(cluster.size_pad_radial);
0935   //  int size_pad_azimuth;
0936   clus->setSizePadAzimuth(cluster.size_pad_azimuth);
0937   //
0938   //
0939   //  //! pad bin size
0940   //  //! phi size per pad in rad
0941   //  double delta_azimuth_bin;
0942   clus->setDeltaAzimuthBin(cluster.delta_azimuth_bin);
0943   //  //! z size per ADC sample bin
0944   //  double delta_z;
0945   clus->setDeltaZ(cluster.delta_z);
0946
0947   // update error matrix
0948
0949   TMatrixF DIM(3, 3);
0950   DIM[0][0] = 0.0;
0951   DIM[0][1] = 0.0;
0952   DIM[0][2] = 0.0;
0953   DIM[1][0] = 0.0;
0954   DIM[1][1] = phi_size * phi_size;  //cluster_v1 expects polar coordinates covariance
0955   DIM[1][2] = 0.0;
0956   DIM[2][0] = 0.0;
0957   DIM[2][1] = 0.0;
0958   DIM[2][2] = z_size * z_size;
0959
0960   TMatrixF ERR(3, 3);
0961   ERR[0][0] = 0.0;
0962   ERR[0][1] = 0.0;
0963   ERR[0][2] = 0.0;
0964   ERR[1][0] = 0.0;
0965   ERR[1][1] = phi_err * phi_err;  //cluster_v1 expects rad, arc, z as elementsof covariance
0966   ERR[1][2] = 0.0;
0967   ERR[2][0] = 0.0;
0968   ERR[2][1] = 0.0;
0969   ERR[2][2] = z_err * z_err;
0970
0971   TMatrixF ROT(3, 3);
0972   ROT[0][0] = cos(clusphi);
0973   ROT[0][1] = -sin(clusphi);
0974   ROT[0][2] = 0.0;
0975   ROT[1][0] = sin(clusphi);
0976   ROT[1][1] = cos(clusphi);
0977   ROT[1][2] = 0.0;
0978   ROT[2][0] = 0.0;
0979   ROT[2][1] = 0.0;
0980   ROT[2][2] = 1.0;
0981
0982   TMatrixF ROT_T(3, 3);
0983   ROT_T.Transpose(ROT);
0984
0985   TMatrixF COVAR_DIM(3, 3);
0986   COVAR_DIM = ROT * DIM * ROT_T;
0987
0988   clus->setSize(0, 0, COVAR_DIM[0][0]);
0989   clus->setSize(0, 1, COVAR_DIM[0][1]);
0990   clus->setSize(0, 2, COVAR_DIM[0][2]);
0991   clus->setSize(1, 0, COVAR_DIM[1][0]);
0992   clus->setSize(1, 1, COVAR_DIM[1][1]);
0993   clus->setSize(1, 2, COVAR_DIM[1][2]);
0994   clus->setSize(2, 0, COVAR_DIM[2][0]);
0995   clus->setSize(2, 1, COVAR_DIM[2][1]);
0996   clus->setSize(2, 2, COVAR_DIM[2][2]);
0997   //cout << " covar_dim[2][2] = " <<  COVAR_DIM[2][2] << endl;
0998
0999   TMatrixF COVAR_ERR(3, 3);
1000   COVAR_ERR = ROT * ERR * ROT_T;
1001
1002   clus->setError(0, 0, COVAR_ERR[0][0]);
1003   clus->setError(0, 1, COVAR_ERR[0][1]);
1004   clus->setError(0, 2, COVAR_ERR[0][2]);
1005   clus->setError(1, 0, COVAR_ERR[1][0]);
1006   clus->setError(1, 1, COVAR_ERR[1][1]);
1007   clus->setError(1, 2, COVAR_ERR[1][2]);
1008   clus->setError(2, 0, COVAR_ERR[2][0]);
1009   clus->setError(2, 1, COVAR_ERR[2][1]);
1010   clus->setError(2, 2, COVAR_ERR[2][2]);
1011
1012   // output prototype specifics
1013
1014   if (Verbosity() >= 2)
1015   {
1016     cout << __PRETTY_FUNCTION__ << "Dump clusters after TpcPrototypeClusterizer" << endl;
1017     clus->identify();
1018   }
1019
1020   return Fun4AllReturnCodes::EVENT_OK;
1021 }
1022
1023 TpcPrototypeUnpacker::PadPlaneData::
1024     PadPlaneData()
1025   : m_data(kMaxPadY, vector<vector<int>>(kMaxPadX, vector<int>(kSAMPLE_LENGTH, 0)))
1026 {
1027 }
1028
1029 void TpcPrototypeUnpacker::PadPlaneData::Reset()
1030 {
1031   for (auto& padrow : m_data)
1032   {
1033     for (auto& pad : padrow)
1034     {
1035       fill(pad.begin(), pad.end(), 0);
1036     }
1037   }
1038
1039   m_groups.clear();
1040 }
1041
1042 bool TpcPrototypeUnpacker::PadPlaneData::IsValidPad(const int pad_x, const int pad_y)
1043 {
1044   return (pad_x >= 0) and
1045          (pad_x < int(kMaxPadX)) and
1046          (pad_y >= 0) and
1047          (pad_y < int(kMaxPadY));
1048 }
1049
1050 vector<int>& TpcPrototypeUnpacker::PadPlaneData::getPad(const int pad_x, const int pad_y)
1051 {
1052   assert(pad_x >= 0);
1053   assert(pad_x < int(kMaxPadX));
1054   assert(pad_y >= 0);
1055   assert(pad_y < int(kMaxPadY));
1056
1057   return m_data[pad_y][pad_x];
1058 }
1059
1060 std::pair<int, int> TpcPrototypeUnpacker::roughZeroSuppression(std::vector<int>& data)
1061 {
1062   std::vector<int> sorted_data(data);
1063
1064   sort(sorted_data.begin(), sorted_data.end());
1065
1066   const int pedestal = sorted_data[sorted_data.size() / 2];
1067   const int max = sorted_data.back();
1068
1069   for (auto& d : data)
1070     d -= pedestal;
1071
1072   return make_pair(pedestal, max);
1073 }
1074
1075 bool operator<(const TpcPrototypeUnpacker::PadPlaneData::SampleID& s1, const TpcPrototypeUnpacker::PadPlaneData::SampleID& s2)
1076 {
1077   if (s1.pad_azimuth == s2.pad_azimuth)
1078   {
1079     if (s1.pad_radial == s2.pad_radial)
1080     {
1081       return s1.sample < s2.sample;
1082     }
1083     else
1084       return s1.pad_radial < s2.pad_radial;
1085   }
1086   else
1087     return s1.pad_azimuth < s2.pad_azimuth;
1088 }
1089
1090 //! 3-D Graph clustering based on PHMakeGroups()
1091 void TpcPrototypeUnpacker::PadPlaneData::Clustering(int zero_suppression, bool verbosity)
1092 {
1093   using namespace boost;
1094   typedef adjacency_list<vecS, vecS, undirectedS> Graph;
1095   typedef bimap<Graph::vertex_descriptor, SampleID> VertexList;
1096
1097   Graph G;
1098   VertexList vertex_list;
1099
1100   for (unsigned int pad_azimuth = 0; pad_azimuth < kMaxPadY; ++pad_azimuth)
1101   {
1102     for (unsigned int pad_radial = 0; pad_radial < kMaxPadX; ++pad_radial)
1103     {
1104       for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
1105       {
1106         if (m_data[pad_azimuth][pad_radial][sample] > zero_suppression)
1107         {
1108           SampleID id{(int) (pad_azimuth), (int) (pad_radial), (int) (sample)};
1109           Graph::vertex_descriptor v = boost::add_vertex(G);
1110           vertex_list.insert(VertexList::value_type(v, id));
1111
1112           add_edge(v, v, G);
1113         }
1114       }  //      for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
1115     }
1116   }  //   for (unsigned int pad_azimuth = 0; pad_azimuth < kMaxPadY; ++pad_azimuth)
1117
1118   // connect 2-D adjacent samples within each pad_x
1119   vector<SampleID> search_directions;
1120   search_directions.push_back(SampleID{0, 0, 1});
1121   //  search_directions.push_back(SampleID{0, 1, 0});
1122   search_directions.push_back(SampleID{1, 0, 0});
1123
1124   for (const auto& it : vertex_list.right)
1125   {
1126     const SampleID id = it.first;
1127     const Graph::vertex_descriptor v = it.second;
1128
1129     for (const SampleID& search_direction : search_directions)
1130     {
1131       //      const SampleID next_id = id + search_direction;
1132       SampleID next_id(id);
1133       next_id.adjust(search_direction);
1134
1135       auto next_it = vertex_list.right.find(next_id);
1136       if (next_it != vertex_list.right.end())
1137       {
1138         add_edge(v, next_it->second, G);
1139       }
1140     }
1141
1142   }  //  for (const auto & it : vertex_list)
1143
1144   // Find the connections between the vertices of the graph (vertices are the rawhits,
1145   // connections are made when they are adjacent to one another)
1146   std::vector<int> component(num_vertices(G));
1147   connected_components(G, &component[0]);
1148
1149   // Loop over the components(vertices) compiling a list of the unique
1150   // connections (ie clusters).
1151   set<int> comps;                // Number of unique components
1152   assert(m_groups.size() == 0);  // no overwrite
1153
1154   for (unsigned int i = 0; i < component.size(); i++)
1155   {
1156     comps.insert(component[i]);
1157     m_groups.insert(make_pair(component[i], vertex_list.left.find(vertex(i, G))->second));
1158   }
1159
1160   //debug prints
1161   if (verbosity)
1162     for (const int& comp : comps)
1163     {
1164       cout << "TpcPrototypeUnpacker::PadPlaneData::Clustering - find cluster " << comp << " containing ";
1165       const auto range = m_groups.equal_range(comp);
1166
1167       for (auto iter = range.first; iter != range.second; ++iter)
1168       {
1169         const SampleID& id = iter->second;
1170         cout << "adc[" << id.pad_azimuth << "][" << id.pad_radial << "][" << id.sample << "] = " << m_data[id.pad_azimuth][id.pad_radial][id.sample] << ", ";
1171       }
1172       cout << endl;
1173     }  //  for (const int& comp : comps)
1174 }
1175
1176 Fun4AllHistoManager*
1177 TpcPrototypeUnpacker::getHistoManager()
1178 {
1179   static string histname("TpcPrototypeUnpacker_HISTOS");
1180
1181   Fun4AllServer* se = Fun4AllServer::instance();
1182   Fun4AllHistoManager* hm = se->getHistoManager(histname);
1183
1184   if (not hm)
1185   {
1186     cout
1187         << "TpcPrototypeUnpacker::get_HistoManager - Making Fun4AllHistoManager "
1188         << histname << endl;
1189     hm = new Fun4AllHistoManager(histname);
1190     se->registerHistoManager(hm);
1191   }
1192
1193   assert(hm);
1194
1195   return hm;
1196 }
1197
1198 void TpcPrototypeUnpacker::registerPadPlane(PHG4TpcPadPlane* inpadplane)
1199 {
1200   cout << "Registering padplane " << endl;
1201   padplane = inpadplane;
1202   padplane->Detector("TPC");
1203   padplane->UpdateInternalParameters();
1204   if (Verbosity())
1205     cout << __PRETTY_FUNCTION__ << " padplane registered and parameters updated" << endl;
1206
1207   return;
1208 }
1209
1210 int TpcPrototypeUnpacker::InitField(PHCompositeNode* topNode)
1211 {
1212   if (Verbosity() >= 1) cout << "TpcPrototypeUnpacker::InitField - create magnetic field setup" << endl;
1213
1214   unique_ptr<PHFieldConfig> default_field_cfg(nullptr);
1215
1216   default_field_cfg.reset(new PHFieldConfigv2(0, 0, 0));
1217
1218   PHField* phfield = PHFieldUtility::GetFieldMapNode(default_field_cfg.get(), topNode, Verbosity() + 1);
1219   assert(phfield);
1220
1221   return Fun4AllReturnCodes::EVENT_OK;
1222 }