Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TpcPrototypeUnpacker.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TpcPrototypeUnpacker.cc
1 /*
2  * TpcPrototypeUnpacker.cc
3  *
4  * Created on: Sep 19, 2018
5  * Author: jinhuang
6  */
7 
8 #include "TpcPrototypeUnpacker.h"
9 
10 #include "ChanMap.h"
11 #include "TpcPrototypeCluster.h"
12 #include "TpcPrototypeDefs.h"
13 
14 #include <g4tpc/PHG4TpcPadPlane.h> // for PHG4TpcPadPlane
15 
16 #include <tpc/TpcDefs.h>
17 
18 #include <g4detectors/PHG4CylinderCellGeom.h>
19 #include <g4detectors/PHG4CylinderCellGeomContainer.h>
20 
21 #include <fun4all/Fun4AllBase.h>
22 #include <fun4all/Fun4AllHistoManager.h>
23 #include <fun4all/Fun4AllReturnCodes.h>
24 #include <fun4all/Fun4AllServer.h>
25 #include <fun4all/PHTFileServer.h>
26 #include <fun4all/SubsysReco.h>
27 
28 #include <phfield/PHFieldConfig.h> // for PHFie...
29 #include <phfield/PHFieldConfigv2.h>
30 #include <phfield/PHFieldUtility.h>
31 
32 #include <tpc/TpcDefs.h>
33 #include <tpc/TpcHit.h>
34 #include <trackbase/TrkrClusterContainer.h>
35 #include <trackbase/TrkrDefs.h> // for hitkey, getLayer
36 #include <trackbase/TrkrHit.h> // for TrkrHit
37 #include <trackbase/TrkrHitSet.h>
38 #include <trackbase/TrkrHitSetContainer.h>
39 
40 #include <phool/PHCompositeNode.h>
41 #include <phool/PHIODataNode.h> // for PHIOD...
42 #include <phool/PHNode.h> // for PHNode
43 #include <phool/PHNodeIterator.h> // for PHNod...
44 #include <phool/PHObject.h> // for PHObject
45 #include <phool/getClass.h>
46 #include <phool/phool.h> // for PHWHERE
47 
48 #include <Event/Event.h>
49 #include <Event/EventTypes.h>
50 #include <Event/packet.h>
51 
52 #include <TAxis.h> // for TAxis
53 #include <TClonesArray.h>
54 #include <TH1.h> // for TH1D
55 #include <TMatrixFfwd.h> // for TMatrixF
56 #include <TMatrixT.h> // for TMatrixT
57 #include <TMatrixTUtils.h> // for TMatr...
58 #include <TNamed.h> // for TNamed
59 #include <TTree.h>
60 
61 #include <boost/bimap.hpp>
62 #include <boost/bind.hpp>
63 #include <boost/format.hpp>
64 #include <boost/graph/adjacency_list.hpp>
65 #include <boost/graph/connected_components.hpp>
66 
67 #include <algorithm>
68 #include <cassert>
69 #include <cmath>
70 #include <cstdint> // for uint32_t
71 #include <iostream>
72 #include <iterator> // for rever...
73 #include <map>
74 #include <memory>
75 #include <sstream>
76 #include <stdexcept>
77 
78 class PHField;
79 
80 using namespace std;
81 using namespace TpcPrototypeDefs::FEEv2;
82 
83 TpcPrototypeUnpacker::TpcPrototypeUnpacker(const std::string& outputfilename)
84  : SubsysReco("TpcPrototypeUnpacker")
85  , padplane(nullptr)
86  , tpcCylinderCellGeom(nullptr)
87  , hitsetcontainer(nullptr)
88  , trkrclusters(nullptr)
89  , m_outputFileName(outputfilename)
90  , m_eventT(nullptr)
91  , m_peventHeader(&m_eventHeader)
92  , m_nClusters(-1)
93  , m_IOClusters(nullptr)
94  , m_chanT(nullptr)
95  , m_pchanHeader(&m_chanHeader)
96  , m_chanData(kSAMPLE_LENGTH, 0)
97  , enableClustering(true)
98  , m_clusteringZeroSuppression(50)
99  , m_nPreSample(5)
100  , m_nPostSample(5)
101  , m_XRayLocationX(-1)
102  , m_XRayLocationY(-1)
103  , m_pdfMaker(nullptr)
104 {
105 }
106 
107 TpcPrototypeUnpacker::~TpcPrototypeUnpacker()
108 {
109  if (m_IOClusters)
110  {
111  m_IOClusters->Clear();
112  delete m_IOClusters;
113  }
114  if (m_pdfMaker)
115  {
116  delete m_pdfMaker;
117  }
118  if (padplane)
119  {
120  delete padplane;
121  }
122 }
123 
124 int TpcPrototypeUnpacker::ResetEvent(PHCompositeNode* topNode)
125 {
126  m_eventHeader = EventHeader();
127  m_padPlaneData.Reset();
128  m_clusters.clear();
129  m_chanHeader = ChannelHeader();
130 
131  m_nClusters = -1;
132  assert(m_IOClusters);
133  m_IOClusters->Clear("C");
134 
135  return Fun4AllReturnCodes::EVENT_OK;
136 }
137 
138 void TpcPrototypeUnpacker::ClusterData::Clear(Option_t*)
139 {
140  pad_radials.clear();
141  pad_azimuths.clear();
142  samples.clear();
143 
144  // for (auto& s : pad_radial_samples) s.second.clear();
145  // pad_radial_samples.clear();
146  //
147  // for (auto& s : pad_azimuth_samples) s.second.clear();
148  // pad_azimuth_samples.clear();
149 
150  while (pad_radial_samples.begin() != pad_radial_samples.end())
151  {
152  pad_radial_samples.begin()->second.clear();
153  pad_radial_samples.erase(pad_radial_samples.begin());
154  }
155 
156  while (pad_azimuth_samples.begin() != pad_azimuth_samples.end())
157  {
158  pad_azimuth_samples.begin()->second.clear();
159  pad_azimuth_samples.erase(pad_azimuth_samples.begin());
160  }
161 
162  while (pad_azimuth_peaks.begin() != pad_azimuth_peaks.end())
163  {
164  pad_azimuth_peaks.erase(pad_azimuth_peaks.begin());
165  }
166 
167  // while (sum_samples.begin() != sum_samples.end())
168  // {
169  // sum_samples.erase(sum_samples.begin());
170  // }
171  sum_samples.clear();
172  sum_samples.shrink_to_fit();
173 }
174 
175 int TpcPrototypeUnpacker::Init(PHCompositeNode* topNode)
176 {
177  if (padplane)
178  padplane->Init(topNode);
179 
180  return Fun4AllReturnCodes::EVENT_OK;
181 }
182 
183 int TpcPrototypeUnpacker::InitRun(PHCompositeNode* topNode)
184 {
185  // Looking for the DST node
186  PHNodeIterator iter(topNode);
187  PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
188  if (!dstNode)
189  {
190  cout << PHWHERE << "DST Node missing, doing nothing." << endl;
191  return Fun4AllReturnCodes::ABORTRUN;
192  }
193 
194  if (padplane)
195  {
196  if (Verbosity() >= VERBOSITY_SOME)
197  {
198  cout << "TpcPrototypeUnpacker::InitRun - making pad plane"
199  << endl;
200  }
201 
202  //setup the constant field
203  const int field_ret = InitField(topNode);
204  if (field_ret != Fun4AllReturnCodes::EVENT_OK)
205  {
206  cout << "TpcPrototypeUnpacker::InitRun- Error - Failed field init with status = " << field_ret << endl;
207  return field_ret;
208  }
209 
210  string seggeonodename = "CYLINDERCELLGEOM_SVTX"; // + detector;
211  tpcCylinderCellGeom = findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, seggeonodename.c_str());
212  if (!tpcCylinderCellGeom)
213  {
214  tpcCylinderCellGeom = new PHG4CylinderCellGeomContainer();
215  PHNodeIterator iter(topNode);
216  PHCompositeNode* runNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "RUN"));
217  PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(tpcCylinderCellGeom, seggeonodename.c_str(), "PHObject");
218  runNode->addNode(newNode);
219  }
220 
221  padplane->InitRun(topNode);
222  padplane->CreateReadoutGeometry(topNode, tpcCylinderCellGeom);
223  }
224 
225  // Create the TrkrHit node if required
226 
227  hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
228  if (!hitsetcontainer)
229  {
230  PHNodeIterator dstiter(dstNode);
231  PHCompositeNode* DetNode =
232  dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
233  if (!DetNode)
234  {
235  DetNode = new PHCompositeNode("TRKR");
236  dstNode->addNode(DetNode);
237  }
238 
239  hitsetcontainer = new TrkrHitSetContainer();
240  PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(hitsetcontainer, "TRKR_HITSET", "PHObject");
241  DetNode->addNode(newNode);
242  }
243 
244  // Create the Cluster node if required
245  trkrclusters = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
246  if (!trkrclusters)
247  {
248 
249  PHNodeIterator dstiter(dstNode);
250  PHCompositeNode* DetNode =
251  dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
252  if (!DetNode)
253  {
254  DetNode = new PHCompositeNode("TRKR");
255  dstNode->addNode(DetNode);
256  }
257 
258  trkrclusters = new TrkrClusterContainer();
259  PHIODataNode<PHObject>* TrkrClusterContainerNode =
260  new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
261  DetNode->addNode(TrkrClusterContainerNode);
262  }
263 
264  if (Verbosity() >= VERBOSITY_SOME)
265  {
266  cout << "TpcPrototypeUnpacker::InitRun - Making PHTFileServer " << m_outputFileName
267  << endl;
268 
269  m_pdfMaker = new SampleFit_PowerLawDoubleExp_PDFMaker();
270  }
271  PHTFileServer::get().open(m_outputFileName, "RECREATE");
272 
273  Fun4AllHistoManager* hm = getHistoManager();
274  assert(hm);
275 
276  TH1D* h = new TH1D("hNormalization", //
277  "Normalization;Items;Summed quantity", 10, .5, 10.5);
278  int i = 1;
279  h->GetXaxis()->SetBinLabel(i++, "Event count");
280  h->GetXaxis()->SetBinLabel(i++, "Collision count");
281  h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit");
282  h->GetXaxis()->SetBinLabel(i++, "TPC G4Hit Edep");
283  h->GetXaxis()->SetBinLabel(i++, "TPC Pad Hit");
284  h->GetXaxis()->SetBinLabel(i++, "TPC Charge e");
285  h->GetXaxis()->SetBinLabel(i++, "TPC Charge fC");
286  h->GetXaxis()->LabelsOption("v");
287  hm->registerHisto(h);
288 
289  m_eventT = new TTree("eventT", "TPC FEE per-event Tree");
290  assert(m_eventT);
291  m_eventT->Branch("evthdr", &m_peventHeader);
292  m_eventT->Branch("nClusters", &m_nClusters, "nClusters/I");
293  m_IOClusters = new TClonesArray("TpcPrototypeUnpacker::ClusterData", 1000);
294  m_eventT->Branch("Clusters", &m_IOClusters);
295 
296  m_chanT = new TTree("chanT", "TPC FEE per-channel Tree");
297  assert(m_chanT);
298  m_chanT->Branch("event", &m_eventHeader.event, "event/I");
299  m_chanT->Branch("chanhdr", &m_pchanHeader);
300  m_chanT->Branch("adc", m_chanData.data(), str(boost::format("adc[%d]/i") % kSAMPLE_LENGTH).c_str());
301 
302  // for (unsigned int layer = m_minLayer; layer <= m_maxLayer; ++layer)
303  // {
304  // const PHG4CylinderCellGeom* layer_geom = seggeo->GetLayerCellGeom(layer);
305 
306  // const string histNameCellHit(boost::str(boost::format{"hCellHit_Layer%1%"} % layer));
307  // const string histNameCellCharge(boost::str(boost::format{"hCellCharge_Layer%1%"} % layer));
308 
309  // }
310 
311  // hm->registerHisto(new TH2D("hLayerCellHit", //
312  // "Number of ADC time-bin hit per channel;Layer ID;Hit number",
313  // m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
314  // 300, -.5, 299.5));
315  // hm->registerHisto(new TH2D("hLayerCellCharge", //
316  // "Charge integrated over drift window per channel;Layer ID;Charge [fC]",
317  // m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
318  // 1000, 0, 1e7 * eplus / (1e-15 * coulomb)));
319  //
320  // hm->registerHisto(new TH2D("hLayerSumCellHit", //
321  // "Number of ADC time-bin hit integrated over channels per layer;Layer ID;Hit number",
322  // m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
323  // 10000, -.5, 99999.5));
324  // hm->registerHisto(new TH2D("hLayerSumCellCharge", //
325  // "Charge integrated over drift window and channel per layer;Layer ID;Charge [fC]",
326  // m_maxLayer - m_minLayer + 1, m_minLayer - .5, m_maxLayer + .5,
327  // 10000, 0, 1000 * 4e6 * eplus / (1e-15 * coulomb)));
328 
329  return Fun4AllReturnCodes::EVENT_OK;
330 }
331 
332 int TpcPrototypeUnpacker::End(PHCompositeNode* topNode)
333 {
334  if (Verbosity() >= VERBOSITY_SOME)
335  {
336  cout << "TpcPrototypeUnpacker::End - write to " << m_outputFileName << endl;
337  }
338  PHTFileServer::get().cd(m_outputFileName);
339 
340  Fun4AllHistoManager* hm = getHistoManager();
341  assert(hm);
342  for (unsigned int i = 0; i < hm->nHistos(); i++)
343  hm->getHisto(i)->Write();
344 
345  // help index files with TChain
346  TTree* T_Index = new TTree("T_Index", "T_Index");
347  assert(T_Index);
348  T_Index->Write();
349 
350  m_eventT->Write();
351  m_chanT->Write();
352 
353  if (m_pdfMaker)
354  {
355  delete m_pdfMaker;
356  m_pdfMaker = nullptr;
357  }
358  return Fun4AllReturnCodes::EVENT_OK;
359 }
360 
361 int TpcPrototypeUnpacker::process_event(PHCompositeNode* topNode)
362 {
363  Fun4AllHistoManager* hm = getHistoManager();
364  assert(hm);
365  TH1D* h_norm = dynamic_cast<TH1D*>(hm->getHisto("hNormalization"));
366  assert(h_norm);
367 
368  Event* event = findNode::getClass<Event>(topNode, "PRDF");
369  if (event == nullptr)
370  {
371  if (Verbosity() >= VERBOSITY_SOME)
372  cout << "TpcPrototypeUnpacker::Process_Event - Event not found" << endl;
373  return Fun4AllReturnCodes::DISCARDEVENT;
374  }
375 
376  if (Verbosity() >= VERBOSITY_SOME)
377  event->identify();
378 
379  // search for data event
380  if (event->getEvtType() == BEGRUNEVENT)
381  {
382  // get_motor_loc(event);
383 
384  return Fun4AllReturnCodes::EVENT_OK;
385  }
386  if (event->getEvtType() != DATAEVENT)
387  return Fun4AllReturnCodes::DISCARDEVENT;
388 
389  m_eventHeader.run = event->getRunNumber();
390  m_eventHeader.event = event->getEvtSequence();
391 
392  // m_eventHeader.xray_x = m_XRayLocationX;
393  // m_eventHeader.xray_y = m_XRayLocationY;
394 
395  if (m_pdfMaker)
396  {
397  m_pdfMaker->MakeSectionPage(str(boost::format("ADC signal fits for Run %1% and event %2%") % m_eventHeader.run % m_eventHeader.event));
398  }
399 
400  static const char* MAX_LENGTH = "MAX_SAMPLES";
401  static const char* IS_PRESENT = "IS_PRESENT";
402  static const char* BX_COUNTER = "BX";
403 
404  unique_ptr<Packet> p(event->getPacket(kPACKET_ID));
405  if (p == nullptr)
406  return Fun4AllReturnCodes::DISCARDEVENT;
407 
408  if (Verbosity() >= VERBOSITY_SOME) p->identify();
409 
410  if (Verbosity() >= VERBOSITY_EVEN_MORE)
411  {
412  cout << "TpcPrototypeUnpacker::process_event - package dump" << endl;
413  p->dump();
414  }
415 
416  int max_length[kN_FEES] = {-1};
417  for (unsigned int i = 0; i < kN_FEES; i++)
418  {
419  try
420  {
421  max_length[i] = p->iValue(i, MAX_LENGTH);
422  }
423  catch (const std::out_of_range& e)
424  {
425  max_length[i] = -1;
426  break;
427  }
428 
429  if (Verbosity() >= VERBOSITY_MORE)
430  {
431  cout << "TpcPrototypeUnpacker::process_event - max_length[" << i << "]=" << max_length[i] << endl;
432  }
433  }
434 
435  m_eventHeader.bx_counter = 0;
436  // bool first_channel = true;
437 
438  for (unsigned int fee = 0; fee < kN_FEES; ++fee)
439  {
440  try
441  {
442  if (!p->iValue(fee, IS_PRESENT))
443  {
444  if (Verbosity() >= VERBOSITY_MORE)
445  {
446  cout << "TpcPrototypeUnpacker::process_event - missing fee[" << fee << "]=" << endl;
447  }
448  continue;
449  }
450  }
451  catch (const std::out_of_range& e)
452  {
453  cout << "TpcPrototypeUnpacker::process_event - out_of_range in p->iValue(" << fee << ", IS_PRESENT)" << endl;
454 
455  break;
456  }
457 
458  if (Verbosity() >= VERBOSITY_MORE)
459  {
460  cout << "TpcPrototypeUnpacker::process_event - processing fee[" << fee << "]=" << endl;
461  }
462 
463  for (unsigned int channel = 0; channel < kN_CHANNELS; channel++)
464  {
465  m_chanHeader = ChannelHeader();
466  m_chanHeader.fee_id = fee;
467  m_chanHeader.size = p->iValue(fee, channel, "NR_SAMPLES"); // number of words until the next channel (header included). this is the real packet_length
468 
469  if (Verbosity() >= VERBOSITY_MORE)
470  {
471  cout << "TpcPrototypeUnpacker::process_event - processing fee["
472  << fee << "], chan[" << channel << "] NR_SAMPLES = "
473  << p->iValue(fee, channel, "NR_SAMPLES") << endl;
474  }
475 
476  unsigned int real_t = 0;
477  uint32_t old_bx_count = 0;
478  uint16_t adc = 0;
479  uint32_t bx_count = 0;
480  m_chanData.resize(kSAMPLE_LENGTH, 0);
481 
482  for (unsigned int t = 0; t < kSAMPLE_LENGTH; t++)
483  {
484  try
485  {
486  adc = p->iValue(fee, channel, t);
487  bx_count = p->iValue(fee, channel, t, BX_COUNTER);
488 
489  if (real_t == 0)
490  m_chanHeader.bx_counter = p->iValue(fee, channel, t, BX_COUNTER);
491  }
492  catch (const std::out_of_range& e)
493  {
494  adc = 0;
495  bx_count = 0;
496  continue;
497  }
498 
499  if (bx_count >= old_bx_count + 128)
500  {
501  old_bx_count = bx_count;
502  real_t = 0;
503  m_chanData.resize(kSAMPLE_LENGTH, 0);
504  }
505 
506  assert(real_t < kSAMPLE_LENGTH);
507  m_chanData[real_t] = adc;
508  // h_raw_wave->Fill(real_t, adc);
509  // h_raw->Fill(real_t, total_chan, adc);
510  // adc_data->push_back(adc);
511  real_t++;
512  } // for (unsigned int t = 3; t < kSAMPLE_LENGTH; t++)
513 
514  // m_chanHeader.packet_type = p->iValue(channel * kPACKET_LENGTH + 2) & 0xffff; // that's the Elink packet type
515  // m_chanHeader.bx_counter = ((p->iValue(channel * kPACKET_LENGTH + 4) & 0xffff) << 4) | (p->iValue(channel * kPACKET_LENGTH + 5) & 0xffff);
516  // m_chanHeader.sampa_address = (p->iValue(channel * kPACKET_LENGTH + 3) >> 5) & 0xf;
517  // m_chanHeader.sampa_channel = p->iValue(channel * kPACKET_LENGTH + 3) & 0x1f;
518  m_chanHeader.fee_channel = channel;
519 
520  // const pair<int, int> pad = SAMPAChan2PadXY(m_chanHeader.fee_channel);
521 
522  m_chanHeader.pad_x = TpcR2Map.GetRpos(fee, channel);
523  m_chanHeader.pad_y = TpcR2Map.Getphipos(fee, channel);
524 
525  //
526  if (Verbosity() >= VERBOSITY_MORE)
527  {
528  cout << "TpcPrototypeUnpacker::process_event - "
529  << "m_chanHeader.m_size = " << int(m_chanHeader.size) << ", "
530  << "m_chanHeader.m_packet_type = " << int(m_chanHeader.packet_type) << ", "
531  << "m_chanHeader.m_bx_counter = " << int(m_chanHeader.bx_counter) << ", "
532  << "m_chanHeader.m_sampa_address = " << int(m_chanHeader.sampa_address) << ", "
533  << "m_chanHeader.m_sampa_channel = " << int(m_chanHeader.sampa_channel) << ", "
534  << "m_chanHeader.m_fee_channel = " << int(m_chanHeader.fee_channel) << ": "
535  << " ";
536 
537  for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
538  {
539  cout << "data[" << sample << "] = " << int(m_chanData[sample]) << " ";
540  }
541 
542  cout << endl;
543  }
544 
545  // fill event data
546  // if (PadPlaneData::IsValidPad(m_chanHeader.pad_x, m_chanHeader.pad_y))
547  // {
548  vector<int>& paddata = m_padPlaneData.getPad(m_chanHeader.pad_x, m_chanHeader.pad_y);
549  //
550  for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
551  {
552  paddata[sample] = int(m_chanData[sample]);
553  }
554  //
555  auto pedestal_max = roughZeroSuppression(paddata);
556  m_chanHeader.pedestal = pedestal_max.first;
557  m_chanHeader.max = pedestal_max.second;
558  // }
559  // output per-channel TTree
560  m_chanT->Fill();
561  } // for (unsigned int channel = 0; channel < kN_CHANNELS; channel++)
562  } // for (unsigned int fee = 0; fee < kN_FEES; ++fee)
563 
564  exportDSTHits();
565 
566  if (enableClustering)
567  {
568  int ret = Clustering();
569 
570  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
571  }
572  h_norm->Fill("Event count", 1);
573  m_eventT->Fill();
574 
575  return Fun4AllReturnCodes::EVENT_OK;
576 }
577 
578 int TpcPrototypeUnpacker::exportDSTHits()
579 {
580  int nHits(0);
581 
582  assert(hitsetcontainer);
583 
584  for (unsigned int pad_radial = 0; pad_radial < kMaxPadX; ++pad_radial)
585  {
586  TrkrDefs::hitsetkey hitsetkey = TpcDefs::genHitSetKey(pad_radial, 0, 0);
587  // Use existing hitset or add new one if needed
588  TrkrHitSetContainer::Iterator hitsetit = hitsetcontainer->findOrAddHitSet(hitsetkey);
589 
590  for (unsigned int pad_azimuth = 0; pad_azimuth < kMaxPadY; ++pad_azimuth)
591  {
592  for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
593  {
594  unsigned int adc = static_cast<unsigned int>(m_padPlaneData.getData()[pad_azimuth][pad_radial][sample]);
595 
596  if (adc)
597  {
598  // generate the key for this hit, requires zbin and phibin
599  TrkrDefs::hitkey hitkey = TpcDefs::genHitKey(pad_azimuth, sample);
600  // See if this hit already exists
601  TrkrHit* hit = nullptr;
602  hit = hitsetit->second->getHit(hitkey);
603  if (!hit)
604  {
605  // create a new one
606  hit = new TpcHit();
607  hitsetit->second->addHitSpecificKey(hitkey, hit);
608  }
609  else
610  {
611  cout << "TpcPrototypeUnpacker::exportDSTHits "
612  << " - Fatal error - hit with "
613  << " hitkey = " << hitkey << " already exist";
614  exit(1);
615  }
616 
617  // Either way, add the energy to it -- adc values will be added at digitization
618  hit->setAdc(adc);
619  ++nHits;
620 
621  } // if (m_padPlaneData.m_data[pad_azimuth][pad_radial][sample])
622 
623  } // for (unsigned int sample = 0; sample < kSAMPLE_LENGTH; sample++)
624  } // for (unsigned int pad_azimuth = 0; pad_azimuth < kMaxPadY; ++pad_azimuth)
625  } // for (unsigned int pad_radial = 0; pad_radial < kMaxPadX; ++pad_radial)
626 
627  return nHits;
628 }
629 
630 int TpcPrototypeUnpacker::Clustering()
631 {
632  if (Verbosity())
633  cout << __PRETTY_FUNCTION__ << " entry" << endl;
634 
635  // find cluster
636  m_padPlaneData.Clustering(m_clusteringZeroSuppression, Verbosity() >= VERBOSITY_SOME);
637  const multimap<int, PadPlaneData::SampleID>& groups = m_padPlaneData.getGroups();
638 
639  // export clusters
640  assert(m_clusters.size() == 0); //already cleared.
641  for (const auto& iter : groups)
642  {
643  const int& i = iter.first;
644  const PadPlaneData::SampleID& id = iter.second;
645  m_clusters[i].pad_radials.insert(id.pad_radial);
646  m_clusters[i].pad_azimuths.insert(id.pad_azimuth);
647  m_clusters[i].samples.insert(id.sample);
648  }
649 
650  // process cluster
651  for (auto& iter : m_clusters)
652  {
653  ClusterData& cluster = iter.second;
654 
655  assert(cluster.pad_radials.size() > 0);
656  assert(cluster.pad_azimuths.size() > 0);
657  assert(cluster.samples.size() > 0);
658 
659  //expand cluster by +/-1 in y
660  if (*(cluster.pad_azimuths.begin()) - 1 >= 0)
661  cluster.pad_azimuths.insert(*(cluster.pad_azimuths.begin()) - 1);
662  if (*(cluster.pad_azimuths.rbegin()) + 1 < (int) kMaxPadY)
663  cluster.pad_azimuths.insert(*(cluster.pad_azimuths.rbegin()) + 1);
664 
665  cluster.min_sample = max(0, *cluster.samples.begin() - m_nPreSample);
666  cluster.max_sample = min((int) (kSAMPLE_LENGTH) -1, *cluster.samples.rbegin() + m_nPostSample);
667  const int n_sample = cluster.max_sample - cluster.min_sample + 1;
668 
669  cluster.sum_samples.assign(n_sample, 0);
670  for (int pad_x = *cluster.pad_radials.begin(); pad_x <= *cluster.pad_radials.rbegin(); ++pad_x)
671  {
672  cluster.pad_radial_samples[pad_x].assign(n_sample, 0);
673  }
674  for (int pad_y = *cluster.pad_azimuths.begin(); pad_y <= *cluster.pad_azimuths.rbegin(); ++pad_y)
675  {
676  cluster.pad_azimuth_samples[pad_y].assign(n_sample, 0);
677  }
678 
679  for (int pad_x = *cluster.pad_radials.begin(); pad_x <= *cluster.pad_radials.rbegin(); ++pad_x)
680  {
681  for (int pad_y = *cluster.pad_azimuths.begin(); pad_y <= *cluster.pad_azimuths.rbegin(); ++pad_y)
682  {
683  assert(m_padPlaneData.IsValidPad(pad_x, pad_y));
684 
685  vector<int>& padsamples = m_padPlaneData.getPad(pad_x, pad_y);
686 
687  for (int i = 0; i < n_sample; ++i)
688  {
689  int adc = padsamples.at(cluster.min_sample + i);
690  cluster.sum_samples[i] += adc;
691  cluster.pad_radial_samples[pad_x][i] += adc;
692  cluster.pad_azimuth_samples[pad_y][i] += adc;
693  }
694 
695  } // for (int pad_y = *cluster.pad_azimuths.begin(); pad_y<=*cluster.pad_azimuths.rbegin() ;++pad_azimuth)
696 
697  } // for (int pad_x = *cluster.pad_radials.begin(); pad_x<=*cluster.pad_radials.rbegin() ;++pad_radial)
698 
699  if (m_pdfMaker)
700  {
701  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));
702  }
703 
704  // fit - overal cluster
705  map<int, double> parameters_constraints;
706  {
707  double peak = NAN;
708  double peak_sample = NAN;
709  double pedstal = NAN;
710  map<int, double> parameters_io;
711  SampleFit_PowerLawDoubleExp(cluster.sum_samples, peak,
712  peak_sample, pedstal, parameters_io, Verbosity());
713 
714  parameters_constraints[1] = parameters_io[1];
715  parameters_constraints[2] = parameters_io[2];
716  parameters_constraints[3] = parameters_io[3];
717  parameters_constraints[5] = parameters_io[5];
718  parameters_constraints[6] = parameters_io[6];
719 
720  cluster.peak = peak;
721  cluster.peak_sample = peak_sample;
722  cluster.pedstal = pedstal;
723  }
724 
725  // fit - X -> radial direction
726  {
727  // double sum_peak = 0;
728  // double sum_peak_pad_radial = 0;
729  // for (int pad_x = *cluster.pad_radials.begin(); pad_x <= *cluster.pad_radials.rbegin(); ++pad_x)
730  // {
731  // double peak = NAN;
732  // double peak_sample = NAN;
733  // double pedstal = NAN;
734  // map<int, double> parameters_io(parameters_constraints);
735  //
736  // SampleFit_PowerLawDoubleExp(cluster.pad_radial_samples[pad_x], peak,
737  // peak_sample, pedstal, parameters_io, Verbosity());
738  //
739  // cluster.pad_radial_peaks[pad_x] = peak;
740  // sum_peak += peak;
741  // sum_peak_pad_radial += peak * pad_x;
742  // }
743  // cluster.avg_pad_radial = sum_peak_pad_radial / sum_peak;
744  // cluster.size_pad_radial = cluster.pad_radials.size();
745 
746  assert(cluster.pad_radials.size() == 1);
747  cluster.avg_pad_radial = *cluster.pad_radials.begin();
748  cluster.size_pad_radial = cluster.pad_radials.size();
749  }
750 
751  // fit - Y -> azimuthal direction
752  {
753  double sum_peak = 0;
754  double sum_peak_pad_azimuth = 0;
755  for (int pad_y = *cluster.pad_azimuths.begin(); pad_y <= *cluster.pad_azimuths.rbegin(); ++pad_y)
756  {
757  double peak = NAN;
758  double peak_sample = NAN;
759  double pedstal = NAN;
760  map<int, double> parameters_io(parameters_constraints);
761 
762  SampleFit_PowerLawDoubleExp(cluster.pad_azimuth_samples[pad_y], peak,
763  peak_sample, pedstal, parameters_io, Verbosity());
764 
765  cluster.pad_azimuth_peaks[pad_y] = peak;
766  sum_peak += peak;
767  sum_peak_pad_azimuth += peak * pad_y;
768  }
769  cluster.avg_pad_azimuth = sum_peak_pad_azimuth / sum_peak;
770  cluster.size_pad_azimuth = cluster.pad_azimuths.size();
771 
772  cluster.min_pad_azimuth = *cluster.pad_azimuths.begin();
773  cluster.max_pad_azimuth = *cluster.pad_azimuths.rbegin();
774  }
775  } // for (auto& iter : m_clusters)
776 
777  // sort by energy
778  map<double, int> cluster_energy;
779  for (auto& iter : m_clusters)
780  {
781  //reverse energy sorting
782  cluster_energy[-iter.second.peak] = iter.first;
783  }
784 
785  // save clusters
786  m_nClusters = 0;
787  assert(m_IOClusters);
788  for (const auto& iter : cluster_energy)
789  {
790  ClusterData& cluster = m_clusters[iter.second];
791 
792  //output to DST clusters
793  cluster.clusterID = m_nClusters; // sync cluster id from cluster container to m_nClusters::ClusterData
794  int ret = exportDSTCluster(cluster, m_nClusters);
795  if (ret != Fun4AllReturnCodes::EVENT_OK) return ret;
796 
797  // super awkward ways of ROOT filling TClonesArray
798  new ((*m_IOClusters)[m_nClusters++]) ClusterData(cluster);
799  }
800  return Fun4AllReturnCodes::EVENT_OK;
801 }
802 
803 int TpcPrototypeUnpacker::exportDSTCluster(ClusterData& cluster, const int iclus)
804 {
805  assert(tpcCylinderCellGeom);
806  assert(trkrclusters);
807 
808  assert(cluster.avg_pad_radial >= 0);
809  assert(cluster.avg_pad_radial < (int) kMaxPadX);
810  const uint8_t layer = static_cast<uint8_t>(cluster.avg_pad_radial);
811  const PHG4CylinderCellGeom* layergeom = tpcCylinderCellGeom->GetLayerCellGeom(layer);
812  const int NPhiBins = layergeom->get_phibins();
813  const int NZBins = layergeom->get_zbins();
814 
815  // create the cluster entry directly in the node tree
816  const TrkrDefs::hitsetkey hit_set_key(TpcDefs::genHitSetKey(layer, 0, 0));
817  const TrkrDefs::cluskey ckey = TpcDefs::genClusKey(hit_set_key, iclus);
818 
819  // make sure this cluster location is available
820  if (trkrclusters->findCluster(ckey))
821  {
822  cout << "TpcPrototypeUnpacker::exportDSTCluster - fatal error - cluster already exist which is not expected for prototype analysis. "
823  << "iclus = " << iclus << ", ckey = " << ckey << ". Offending cluster: " << endl;
824  trkrclusters->findCluster(ckey)->identify();
825  exit(1);
826  }
827  TpcPrototypeCluster* clus = new TpcPrototypeCluster();
828  assert(clus);
829  clus->setClusKey(ckey);
830  trkrclusters->addCluster(clus);
831 
832  // calculate geometry
833  const double radius = layergeom->get_radius(); // returns center of layer
834  if (cluster.avg_pad_azimuth < 0)
835  {
836  cout << __PRETTY_FUNCTION__ << " WARNING - cluster.avg_pad_azimuth = " << cluster.avg_pad_azimuth << " < 0"
837  << ", cluster.avg_pad_radial = " << cluster.avg_pad_radial << endl;
838 
839  return Fun4AllReturnCodes::ABORTEVENT;
840  }
841  if (cluster.avg_pad_azimuth >= NPhiBins)
842  {
843  cout << __PRETTY_FUNCTION__ << " WARNING - cluster.avg_pad_azimuth = " << cluster.avg_pad_azimuth << " > " << NPhiBins
844  << ", cluster.avg_pad_radial = " << cluster.avg_pad_radial << endl;
845 
846  return Fun4AllReturnCodes::ABORTEVENT;
847  }
848  const int lowery = floor(cluster.avg_pad_azimuth);
849 
850  if (lowery < 0)
851  {
852  cout << __PRETTY_FUNCTION__ << " WARNING - cluster.avg_pad_azimuth = " << cluster.avg_pad_azimuth << " -> "
853  << " lower = " << lowery << " < 0"
854  << ", cluster.avg_pad_radial = " << cluster.avg_pad_radial << endl;
855 
856  return Fun4AllReturnCodes::ABORTEVENT;
857  }
858 
859  const double clusphi = layergeom->get_phicenter(lowery) //
860  + (layergeom->get_phicenter(lowery + 1) - layergeom->get_phicenter(lowery)) //
861  * (cluster.avg_pad_azimuth - lowery);
862 
863  assert(cluster.min_sample >= 0);
864  assert(cluster.min_sample + cluster.peak_sample < NZBins);
865  const double clusz = layergeom->get_zcenter(cluster.min_sample) //
866  + (layergeom->get_zcenter(cluster.min_sample + 1) - layergeom->get_zcenter(cluster.min_sample)) * cluster.peak_sample;
867 
868  const double phi_size = cluster.size_pad_azimuth; // * radius * layergeom->get_phistep();
869  const double z_size = (cluster.max_sample - cluster.min_sample); // * layergeom->get_zstep();
870 
871  static const double phi_err = 170e-4;
872 
873  static const double z_err = 1000e-4;
874 
875  // Fill in the cluster details
876  //================
877  clus->setAdc(cluster.peak);
878  clus->setPosition(0, radius * cos(clusphi));
879  clus->setPosition(1, radius * sin(clusphi));
880  clus->setPosition(2, clusz);
881  clus->setGlobal();
882 
883  //update cluster positions
884  cluster.avg_pos_x = clus->getPosition(0);
885  cluster.avg_pos_y = clus->getPosition(1);
886  cluster.avg_pos_z = clus->getPosition(2);
887 
888  cluster.delta_z = (layergeom->get_zcenter(cluster.min_sample + 1) - layergeom->get_zcenter(cluster.min_sample));
889  cluster.delta_azimuth_bin = (layergeom->get_phicenter(lowery + 1) - layergeom->get_phicenter(lowery));
890 
891  // output prototype specifics
892 
893  // std::set<int> pad_radials;
894  clus->setPadRadials(cluster.pad_radials);
895  // std::set<int> pad_azimuths;
896  clus->setPadAzimuths(cluster.pad_azimuths);
897  // std::set<int> samples;
898  clus->setSamples(cluster.samples);
899  //
900  // std::map<int, std::vector<double>> pad_radial_samples;
901  clus->setPadRadialSamples(cluster.pad_radial_samples);
902  // std::map<int, std::vector<double>> pad_azimuth_samples;
903  clus->setPadAzimuthSamples(cluster.pad_azimuth_samples);
904  // std::vector<double> sum_samples;
905  clus->setSumSamples(cluster.sum_samples);
906  //
907  // int min_sample;
908  clus->setMinSample(cluster.min_sample);
909  // int max_sample;
910  clus->setMaxSample(cluster.max_sample);
911  // int min_pad_azimuth;
912  clus->setMinPadAzimuth(cluster.min_pad_azimuth);
913  // int max_pad_azimuth;
914  clus->setMaxPadAzimuth(cluster.max_pad_azimuth);
915  //
916  // double peak;
917  clus->setPeak(cluster.peak);
918  // double peak_sample;
919  clus->setPeakSample(cluster.peak_sample);
920  // double pedstal;
921  clus->setPedstal(cluster.pedstal);
922  //
923  // std::map<int, double> pad_azimuth_peaks;
924  clus->setPadAzimuthPeaks(cluster.pad_azimuth_peaks);
925  //
926  // //! pad coordinate
927  // int avg_pad_radial;
928  clus->setAvgPadRadial(cluster.avg_pad_radial);
929  // double avg_pad_azimuth;
930  clus->setAvgPadAzimuth(cluster.avg_pad_azimuth);
931  //
932  // //! cluster size in units of pad bins
933  // int size_pad_radial;
934  clus->setSizePadRadial(cluster.size_pad_radial);
935  // int size_pad_azimuth;
936  clus->setSizePadAzimuth(cluster.size_pad_azimuth);
937  //
938  //
939  // //! pad bin size
940  // //! phi size per pad in rad
941  // double delta_azimuth_bin;
942  clus->setDeltaAzimuthBin(cluster.delta_azimuth_bin);
943  // //! z size per ADC sample bin
944  // double delta_z;
945  clus->setDeltaZ(cluster.delta_z);
946 
947  // update error matrix
948 
949  TMatrixF DIM(3, 3);
950  DIM[0][0] = 0.0;
951  DIM[0][1] = 0.0;
952  DIM[0][2] = 0.0;
953  DIM[1][0] = 0.0;
954  DIM[1][1] = phi_size * phi_size; //cluster_v1 expects polar coordinates covariance
955  DIM[1][2] = 0.0;
956  DIM[2][0] = 0.0;
957  DIM[2][1] = 0.0;
958  DIM[2][2] = z_size * z_size;
959 
960  TMatrixF ERR(3, 3);
961  ERR[0][0] = 0.0;
962  ERR[0][1] = 0.0;
963  ERR[0][2] = 0.0;
964  ERR[1][0] = 0.0;
965  ERR[1][1] = phi_err * phi_err; //cluster_v1 expects rad, arc, z as elementsof covariance
966  ERR[1][2] = 0.0;
967  ERR[2][0] = 0.0;
968  ERR[2][1] = 0.0;
969  ERR[2][2] = z_err * z_err;
970 
971  TMatrixF ROT(3, 3);
972  ROT[0][0] = cos(clusphi);
973  ROT[0][1] = -sin(clusphi);
974  ROT[0][2] = 0.0;
975  ROT[1][0] = sin(clusphi);
976  ROT[1][1] = cos(clusphi);
977  ROT[1][2] = 0.0;
978  ROT[2][0] = 0.0;
979  ROT[2][1] = 0.0;
980  ROT[2][2] = 1.0;
981 
982  TMatrixF ROT_T(3, 3);
983  ROT_T.Transpose(ROT);
984 
985  TMatrixF COVAR_DIM(3, 3);
986  COVAR_DIM = ROT * DIM * ROT_T;
987 
988  clus->setSize(0, 0, COVAR_DIM[0][0]);
989  clus->setSize(0, 1, COVAR_DIM[0][1]);
990  clus->setSize(0, 2, COVAR_DIM[0][2]);
991  clus->setSize(1, 0, COVAR_DIM[1][0]);
992  clus->setSize(1, 1, COVAR_DIM[1][1]);
993  clus->setSize(1, 2, COVAR_DIM[1][2]);
994  clus->setSize(2, 0, COVAR_DIM[2][0]);
995  clus->setSize(2, 1, COVAR_DIM[2][1]);
996  clus->setSize(2, 2, COVAR_DIM[2][2]);
997  //cout << " covar_dim[2][2] = " << COVAR_DIM[2][2] << endl;
998 
999  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 
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 }