Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TpcRawWriter.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TpcRawWriter.cc
1 #include "TpcRawWriter.h"
2 
3 #include <trackbase/TpcDefs.h>
4 #include <trackbase/MvtxDefs.h>
5 #include <trackbase/InttDefs.h>
6 #include <micromegas/MicromegasDefs.h>
7 
8 #include <trackbase/TrkrClusterContainerv4.h>
9 #include <trackbase/TrkrCluster.h>
10 #include <trackbase/TrkrClusterHitAssocv3.h>
11 #include <trackbase/TrkrDefs.h> // for hitkey, getLayer
12 #include <trackbase/TrkrHit.h>
13 #include <trackbase/RawHitv1.h>
14 #include <trackbase/RawHitTpc.h>
15 #include <trackbase/TrkrHitSet.h>
16 #include <trackbase/RawHitSet.h>
17 #include <trackbase/RawHitSetv1.h>
18 #include <trackbase/TrkrHitSetContainer.h>
19 #include <trackbase/RawHitSetContainer.h>
20 #include <trackbase/RawHitSetContainerv1.h>
21 #include <trackbase/ActsGeometry.h>
22 #include <fun4all/Fun4AllReturnCodes.h>
23 #include <fun4all/SubsysReco.h> // for SubsysReco
24 
25 #include <g4detectors/PHG4TpcCylinderGeom.h>
26 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
27 
28 #include <Acts/Definitions/Units.hpp>
29 #include <Acts/Surfaces/Surface.hpp>
30 
31 #include <phool/PHCompositeNode.h>
32 #include <phool/PHIODataNode.h> // for PHIODataNode
33 #include <phool/PHNode.h> // for PHNode
34 #include <phool/PHNodeIterator.h>
35 #include <phool/PHObject.h> // for PHObject
36 #include <phool/getClass.h>
37 #include <phool/phool.h> // for PHWHERE
38 
39 #include <TMatrixFfwd.h> // for TMatrixF
40 #include <TMatrixT.h> // for TMatrixT, ope...
41 #include <TMatrixTUtils.h> // for TMatrixTRow
42 
43 #include <TFile.h>
44 
45 #include <cmath> // for sqrt, cos, sin
46 #include <iostream>
47 #include <map> // for _Rb_tree_cons...
48 #include <string>
49 #include <utility> // for pair
50 #include <array>
51 #include <vector>
52 #include <limits>
53 // Terra incognita....
54 #include <pthread.h>
55 
56 
57 
58 TpcRawWriter::TpcRawWriter(const std::string &name)
59  : SubsysReco(name)
60 { std::cout << PHWHERE << "Construct TpcRawWriter" << std::endl; }
61 
62 
63 int TpcRawWriter::InitRun(PHCompositeNode *topNode)
64 {
65  if(topNode)
66  std::cout << PHWHERE << "Init TpcRawWriter" << std::endl;
67  PHNodeIterator iter(topNode);
68 
69  // Looking for the DST node
70  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
71  if (!dstNode)
72  {
73  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
74  return Fun4AllReturnCodes::ABORTRUN;
75  }
76 
77  // Create the Cluster node if required
78  auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
79  if (!trkrclusters)
80  {
81  PHNodeIterator dstiter(dstNode);
82  PHCompositeNode *DetNode =
83  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
84  if (!DetNode)
85  {
86  DetNode = new PHCompositeNode("TRKR");
87  dstNode->addNode(DetNode);
88  }
89 
90  trkrclusters = new TrkrClusterContainerv4;
91  PHIODataNode<PHObject> *TrkrClusterContainerNode =
92  new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
93  DetNode->addNode(TrkrClusterContainerNode);
94  }
95 
96  auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
97  if (!clusterhitassoc)
98  {
99  PHNodeIterator dstiter(dstNode);
100  PHCompositeNode *DetNode =
101  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
102  if (!DetNode)
103  {
104  DetNode = new PHCompositeNode("TRKR");
105  dstNode->addNode(DetNode);
106  }
107 
108  clusterhitassoc = new TrkrClusterHitAssocv3;
109  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(clusterhitassoc, "TRKR_CLUSTERHITASSOC", "PHObject");
110  DetNode->addNode(newNode);
111  }
112 
113  // new containers
114  m_rawhits = findNode::getClass<RawHitSetContainerv1>(topNode, "TRKR_RAWHITSET");
115  if (!m_rawhits)
116  {
117  PHNodeIterator dstiter(dstNode);
118  auto DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
119  if (!DetNode)
120  {
121  DetNode = new PHCompositeNode("TRKR");
122  dstNode->addNode(DetNode);
123  }
124 
125  m_rawhits = new RawHitSetContainerv1;
126  auto newNode = new PHIODataNode<PHObject>(m_rawhits, "TRKR_RAWHITSET", "PHObject");
127  DetNode->addNode(newNode);
128  }
129 
130  return Fun4AllReturnCodes::EVENT_OK;
131 }
132 
133 int TpcRawWriter::process_event(PHCompositeNode *topNode)
134 {
135  // int print_layer = 18;
136 
137  // if (Verbosity() > 1000)
138  if(topNode)
139  std::cout << "TpcRawWriter::Process_Event" << std::endl;
140 
141  PHNodeIterator iter(topNode);
142  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
143  if (!dstNode)
144  {
145  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
146  return Fun4AllReturnCodes::ABORTRUN;
147  }
148 
149  // get node containing the digitized hits
150  m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
151  if (!m_hits)
152  {
153  std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
154  return Fun4AllReturnCodes::ABORTRUN;
155  }
156 
157  // new containers
158  m_rawhits = findNode::getClass<RawHitSetContainerv1>(topNode, "TRKR_RAWHITSET");
159  if (!m_rawhits)
160  {
161  std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
162  return Fun4AllReturnCodes::ABORTRUN;
163  }
164 
165  // get node for clusters
166  m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
167  if (!m_clusterlist)
168  {
169  std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << std::endl;
170  return Fun4AllReturnCodes::ABORTRUN;
171  }
172 
173  // get node for cluster hit associations
174  m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
175  if (!m_clusterhitassoc)
176  {
177  std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << std::endl;
178  return Fun4AllReturnCodes::ABORTRUN;
179  }
180 
181  PHG4TpcCylinderGeomContainer *geom_container =
182  findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
183  if (!geom_container)
184  {
185  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
186  return Fun4AllReturnCodes::ABORTRUN;
187  }
188 
189  m_tGeometry = findNode::getClass<ActsGeometry>(topNode,
190  "ActsGeometry");
191  if(!m_tGeometry)
192  {
193  std::cout << PHWHERE
194  << "ActsTrackingGeometry not found on node tree. Exiting"
195  << std::endl;
196  return Fun4AllReturnCodes::ABORTRUN;
197  }
198 
199 
200  // The hits are stored in hitsets, where each hitset contains all hits in a given TPC readout (layer, sector, side), so clusters are confined to a hitset
201  // The TPC clustering is more complicated than for the silicon, because we have to deal with overlapping clusters
202 
203  // loop over the MVTX HitSet objects
204  std::cout << "processing mvtx" << std::endl;
205  TrkrHitSetContainer::ConstRange mvtxhitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::mvtxId);
206  // const int num_hitsets = std::distance(hitsetrange.first,hitsetrange.second);
207 
208  int count = 0;
209  for (TrkrHitSetContainer::ConstIterator hitsetitr = mvtxhitsetrange.first;
210  hitsetitr != mvtxhitsetrange.second;
211  ++hitsetitr)
212  {
213  TrkrHitSet *hitset = hitsetitr->second;
214 
215  m_rawhits->findOrAddHitSet(hitsetitr->first);
216  RawHitSetv1 *rhitset = dynamic_cast<RawHitSetv1 *>(m_rawhits->findHitSet(hitsetitr->first));
217  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
218 
219  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
220  hitr != hitrangei.second;
221  ++hitr){
222  unsigned short adc = (unsigned short)(hitr->second->getAdc());
223  if(adc>0){
224  RawHitv1 *rHit = new RawHitv1;
225  unsigned short iphi = MvtxDefs::getCol(hitr->first);
226  unsigned short it = MvtxDefs::getRow(hitr->first);
227  rHit->setPhiBin(iphi);
228  rHit->setTBin(it);
229  rHit->setAdc(adc);
230  rhitset->addHit(rHit);
231  }
232  }
233 
234  count++;
235  }
236  std::cout << "processing intt" << std::endl;
237  // loop over the INTT HitSet objects
238  TrkrHitSetContainer::ConstRange intt_hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::inttId);
239  // const int num_hitsets = std::distance(intt_hitsetrange.first,intt_hitsetrange.second);
240 
241  for (TrkrHitSetContainer::ConstIterator hitsetitr = intt_hitsetrange.first;
242  hitsetitr != intt_hitsetrange.second;
243  ++hitsetitr)
244  {
245  TrkrHitSet *hitset = hitsetitr->second;
246 
247  m_rawhits->findOrAddHitSet(hitsetitr->first);
248  RawHitSet *rhitset = m_rawhits->findHitSet(hitsetitr->first);
249  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
250 
251  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
252  hitr != hitrangei.second;
253  ++hitr){
254  //int sector = InttDefs::getLadderPhiId(hitsetitr->first);
255  //int side = InttDefs::getLadderZId(hitsetitr->first);
256  //int layer = TrkrDefs::getLayer(hitsetitr->first);
257  unsigned short adc = (unsigned short)(hitr->second->getAdc());
258  if(adc>0){
259  RawHitv1 *rHit = new RawHitv1;
260  unsigned short iphi = InttDefs::getCol(hitr->first);
261  unsigned short it = InttDefs::getRow(hitr->first);
262  //std::cout << " layer " << layer << " sector: " << sector << " side " << side << " col: " << iphi << " row " << it << std::endl;
263  rHit->setPhiBin(iphi);
264  rHit->setTBin(it);
265  rHit->setAdc(adc);
266  rhitset->addHit(rHit);
267  }
268  }
269 
270  count++;
271  }
272  std::cout << "processing mm" << std::endl;
273  // loop over the micromega HitSet objects
274  TrkrHitSetContainer::ConstRange mm_hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::micromegasId);
275  // const int num_hitsets = std::distance(mm_hitsetrange.first,mm_hitsetrange.second);
276 
277  for (TrkrHitSetContainer::ConstIterator hitsetitr = mm_hitsetrange.first;
278  hitsetitr != mm_hitsetrange.second;
279  ++hitsetitr)
280  {
281  TrkrHitSet *hitset = hitsetitr->second;
282 
283  m_rawhits->findOrAddHitSet(hitsetitr->first);
284  RawHitSet *rhitset = m_rawhits->findHitSet(hitsetitr->first);
285  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
286 
287  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
288  hitr != hitrangei.second;
289  ++hitr){
290 
291  unsigned short adc = (unsigned short)(hitr->second->getAdc());
292  if(adc>0){
293  RawHitv1 *rHit = new RawHitv1;
294  // TrkrDefs::hitkey tmp = (hitr->first >> MicromegasDefs::kBitShiftStrip);
295  unsigned short iphi = MicromegasDefs::getStrip( hitr->first );
296  unsigned short it = 0;
297  rHit->setPhiBin(iphi);
298  rHit->setTBin(it);
299  rHit->setAdc(adc);
300  rhitset->addHit(rHit);
301  }
302  }
303 
304  count++;
305  }
306  std::cout << "processing tpc" << std::endl;
307  // loop over the TPC HitSet objects
308  TrkrHitSetContainer::ConstRange tpc_hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
309  // const int num_hitsets = std::distance(hitsetrange.first,hitsetrange.second);
310 
311  int ncheck = 0;
312  int allbad = 0;
313 
314  for (TrkrHitSetContainer::ConstIterator hitsetitr = tpc_hitsetrange.first;
315  hitsetitr != tpc_hitsetrange.second;
316  ++hitsetitr)
317  //TrkrHitSetContainer::ConstIterator hitsetitr = tpc_hitsetrange.first;
318  {
319  TrkrHitSet *hitset = hitsetitr->second;
320  unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
321  // if(layer!=7) continue;
322  // int side = TpcDefs::getSide(hitsetitr->first);
323  unsigned int sector= TpcDefs::getSectorId(hitsetitr->first);
324  PHG4TpcCylinderGeom *layergeom = geom_container->GetLayerCellGeom(layer);
325  unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
326  unsigned short NPhiBinsSector = NPhiBins/12;
327  unsigned short NZBins = (unsigned short)layergeom->get_zbins();
328  //unsigned short NZBinsSide = NZBins/2;
329  unsigned short NZBinsSide = NZBins;
330  unsigned short NZBinsMin = 0;
331  unsigned short PhiOffset = NPhiBinsSector * sector;
332  unsigned short ZOffset = NZBinsMin;
333  std::vector<int> nhits_thispad;
334  std::vector<std::vector<unsigned short>> adcval;
335  adcval.resize(NPhiBins);
336  for(int nphi= 0; nphi < NPhiBins;nphi++){
337  for(int nz= 0; nz < NZBins;nz++){
338  adcval[nphi].push_back(0);
339  nhits_thispad.push_back(0);
340  }
341  }
342  //std::vector<std::vector<unsigned short>> outadc;
343  //outadc.resize(NPhiBins);
344 
345  // RawHitSetContainerv1::Iterator rawhitset_iter =
346  m_rawhits->findOrAddHitSet(hitsetitr->first);
347  RawHitSetv1 *rhitset = dynamic_cast<RawHitSetv1 *>(m_rawhits->findHitSet(hitsetitr->first));
348  rhitset->setTpcPhiBins(NPhiBins);
349 
350  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
351  int zbinmax = 498;
352  int zbinmin = 0;
353  if(layer>=7 && layer <22){
354  int etacut = 249 - ((50+(layer-7))/105.5)*249;
355  zbinmin = etacut;
356  zbinmax -= etacut;
357  }
358  if(layer>=22 && layer <=48){
359  int etacut = 249 - ((65+((40.5/26)*(layer-22)))/105.5)*249;
360  zbinmin = etacut;
361  zbinmax -= etacut;
362  }
363 
364  // std::cout << " layer: " << layer << " zbin limit " << zbinmin << " | " << zbinmax <<std::endl;
365  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
366  hitr != hitrangei.second;
367  ++hitr){
368  if( TpcDefs::getPad(hitr->first) - PhiOffset < 0 ){
369  continue;
370  }
371  if( TpcDefs::getTBin(hitr->first) - ZOffset < 0 ){
372  //std::cout << "WARNING zbin out of range: " << TpcDefs::getTBin(hitr->first) - zoffset << " | " << zbins <<std::endl;
373  continue;
374  }
375  unsigned short phibin = TpcDefs::getPad(hitr->first) - PhiOffset;
376  unsigned short zbin = TpcDefs::getTBin(hitr->first) - ZOffset;
377  unsigned short zbinorg = TpcDefs::getTBin(hitr->first);
378  if(phibin>=NPhiBinsSector){
379  //std::cout << "WARNING phibin out of range: " << phibin << " | " << phibins << std::endl;
380  continue;
381  }
382  if(zbin>=NZBinsSide){
383  //std::cout << "WARNING z bin out of range: " << zbin << " | " << zbins << std::endl;
384  continue;
385  }
386  if(zbinorg>zbinmax||zbinorg<zbinmin)
387  continue;
388  float_t fadc = (hitr->second->getAdc()) - pedestal; // proper int rounding +0.5
389  unsigned short adc = 0;
390  if(fadc>0) adc = (unsigned short) fadc;
391 
392  if(adc>0){
393  nhits_thispad[phibin]++;
394  unsigned short it = TpcDefs::getTBin(hitr->first);
395  adcval[phibin][it] = adc;
396  //if(layer==20&&sector==4&&phibin==77)
397  //std::cout << "input sector: " << sector << " lay: " << layer << " # " << phibin << "|"<< it <<" adc: " << (int)adcval[phibin][it] << " pack: " << std::endl;
398  /*RawHitTpc *rHit = new RawHitTpc;
399  unsigned short it = TpcDefs::getTBin(hitr->first);
400  rHit->setTBin(it);
401  unsigned short iadc = adc;
402  rHit->setAdc(iadc);
403  rhitset->addTpcHit(phibin,rHit);
404  */
405  }
406  }
407  // assert(adcval);
408  for(int nphi= 0; nphi < NPhiBins;nphi++){
409  for(int nz= 0; nz < NZBins;nz++){
410  /*if(nphi>=1&&nz<NZBins-1)//Remove isolated single timebin hits
411  if(layer==29&&sector==0&&nphi==0){
412  if(nz>209&&nz<215){
413  std::cout << "single check sector: " << sector << " lay: " << layer << " # " << nphi << "|"<< nz <<" adc: " << (int)adcval[nphi][nz] << " pack: " << std::endl;
414  }
415  }
416  */
417  if(adcval[nphi][nz]>0)
418  if((adcval[nphi][nz-1]==0)&&(adcval[nphi][nz+1]==0)){
419  // std::cout << "single hit sector: " << sector << " lay: " << layer << " # " << nphi << "|"<< nz <<" adc: " << (int)adcval[nphi][nz] << " pack: " << std::endl;
420  adcval[nphi][nz] = 0;
421  }
422  }
423  }
424  for(int nphi= 0; nphi < NPhiBins;nphi++){
425  int zero_count = 0;
426  // int outpos = 0;
427  if(nhits_thispad[nphi]==0)continue;
428  for(int nz= 0; nz < NZBins;nz++){
429 
430  uint8_t thisadc = 0;
431  if(nphi>=1&&nz<NZBins-1)//Remove isolated single timebin hits
432  if(adcval[nphi][nz]>0)
433  if((adcval[nphi][nz-1]==0)&&(adcval[nphi][nz+1]==0))
434  adcval[nphi][nz] = 0;
435 
436  if(adcval[nphi][nz]>255){//take care of overflows, we want to store 8 bit ADC values
437  thisadc=255;
438  // adcval[nphi][nz]=255;
439  }
440  else{
441  thisadc = adcval[nphi][nz];
442  }
443  //std::cout << "nz: " << nz << " adc " << (int)thisadc << std::endl;
444  if(thisadc>0){//non zero adc, fill adc
445  zero_count = 0;
446  rhitset->m_tpchits[nphi].push_back(thisadc);
447  if(layer==222&&sector==6&&nphi==10)std::cout << "0#"<<nphi<<"nz= " << nz << " filling " << (int)thisadc << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
448  }
449  else{
450  zero_count++;
451  // if(nz==0){
452  // rhitset->m_tpchits[nphi].push_back(0);
453  //std::cout << " 1filling " << 0 << "|" << rhitset->m_tpchits[nphi].back() << " | " << rhitset->m_tpchits[nphi][outpos++] << std::endl;
454  // if(layer==222&&sector==6&&nphi==10)std::cout << "1#"<<nphi<<"nz= " << nz << " filling " << 0 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
455  // }else
456  /*
457  if(adcval[nphi][nz-1]>0){ //trailing edge 0
458  rhitset->m_tpchits[nphi].push_back(0);
459  if(layer==222&&sector==6&&nphi==10)std::cout << "2#"<<nphi<<"nz= " << nz << " filling " << 0 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
460  //std::cout << " 2filling " << 0 << "|" << rhitset->m_tpchits[nphi].back() << " | " << rhitset->m_tpchits[nphi][outpos++] << std::endl;
461  }
462  */
463  }
464  if(zero_count>1){
465  if(nz<NZBins-2){
466  if(adcval[nphi][nz+2]>0&&adcval[nphi][nz+1]==0){//fill zero count, end of zero series
467  rhitset->m_tpchits[nphi].push_back(zero_count-1);
468  zero_count = 0;
469  if(layer==222&&sector==6&&nphi==10)
470  std::cout << "3#"<<nphi<<"nz= " << nz << " filling " << zero_count-1 << " at " << rhitset->m_tpchits[nphi].size() << " nz+2 " << (nz+2) << " adcval[nz+2]" << (int) adcval[nphi][nz+2] << std::endl;
471  }
472  if(adcval[nphi][nz+1]>0){
473  rhitset->m_tpchits[nphi].push_back(0);//leading edge zero
474  if(layer==222&&sector==6&&nphi==10)std::cout << "4#"<<nphi<<"nz= " << nz << " filling " << 0 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
475 
476  //std::cout << " 4filling " << 0 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
477  }
478  }
479 
480  }else{
481  if(zero_count==1){
482  rhitset->m_tpchits[nphi].push_back(0);
483  if(layer==222&&sector==6&&nphi==10)std::cout << "5#"<<nphi<<"nz= " << nz << " filling " << 0 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
484  //std::cout << " 5filling " << 0 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
485  }
486  }
487  if(zero_count == 254){
488  zero_count = 0;
489  rhitset->m_tpchits[nphi].push_back(254-2);
490  rhitset->m_tpchits[nphi].push_back(0);
491  if(layer==222&&sector==6&&nphi==10)std::cout << "6#"<<nphi<<"nz= " << nz << " filling " << 254-1 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
492  //std::cout << "6 filling " << 254-1 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
493  // rhitset->m_tpchits[nphi].push_back(0);
494  // if(layer==222&&sector==6&&nphi==10)std::cout << "7#"<<nphi<<"nz= " << nz << " filling " << 0 << " at " << rhitset->m_tpchits[nphi].size() << std::endl;
495  //std::cout << " 7filling " << 0 << "|" << ((int)rhitset->m_tpchits[nphi].back()) << " | " << ((int)rhitset->m_tpchits[nphi][outpos++]) << std::endl;
496  }
497  /*
498  if(layer%15==0){
499  std::cout << "nz #" << nz << " adc " << adcval[nphi][nz] << " n zero " << zero_count << " last (" << rhitset->m_tpchits[nphi].size() << ") " << ((int)rhitset->m_tpchits[nphi].back()) << std::endl;
500  if(rhitset->m_tpchits[nphi].size()>=3)
501  std::cout << " [ " << rhitset->m_tpchits[nphi].size()-1 << ": "<< ((int)rhitset->m_tpchits[nphi][rhitset->m_tpchits[nphi].size()-1]) << " ] "
502  << " [ " << rhitset->m_tpchits[nphi].size()-2 << ": "<< ((int)rhitset->m_tpchits[nphi][rhitset->m_tpchits[nphi].size()-2]) << " ] "
503  << " [ " << rhitset->m_tpchits[nphi].size()-3<< ": "<< ((int)rhitset->m_tpchits[nphi][rhitset->m_tpchits[nphi].size()-3]) << " ] "
504  << std::endl;
505  }
506  */
507  }
508  if(zero_count>0&&rhitset->m_tpchits[nphi].back()!=0)rhitset->m_tpchits[nphi].push_back(0); //pad zero at the end
509  // std::cout << "sector: " << sector << " lay: " << layer << " nphi: " << nphi << "|" << NPhiBins<< " nzfilled: " << rhitset->m_tpchits[nphi].size() << std::endl;
510  }
511  count++;
512 
513  // std::cout << "unpacking..." << std::endl;
514  // unpack and check for correctness
515  std::vector<std::vector<uint8_t>> outval;
516  outval.resize(NPhiBins);
517  for(int nphi= 0; nphi < NPhiBins;nphi++){
518  outval[nphi].resize(NZBins,0);
519  for(int nz= 0; nz < NZBins;nz++){
520  if(outval[nphi][nz]!=0)
521  std::cout << "WARNING!" << std::endl;
522  }
523  }
524  // now we have a clean output array
525  for(int nphi= 0; nphi < NPhiBins;nphi++){
526  if(rhitset->m_tpchits[nphi].size()==0) continue;
527 
528  int pindex = 0;
529  for(unsigned int nzo = 0;nzo<rhitset->m_tpchits[nphi].size();nzo++){
530  uint8_t val = rhitset->m_tpchits[nphi][nzo];
531 
532  if(val==0)
533  pindex++;
534  else{
535  if(nzo==0){
536  outval[nphi][pindex++]=val;
537  }else{
538  if((rhitset->m_tpchits[nphi][nzo-1]==0)&&(rhitset->m_tpchits[nphi][nzo+1]==0))//found zero count
539  pindex+=val;
540  else{
541  outval[nphi][pindex++]=val;
542  }
543  }
544  }
545  }
546  }
547  {
548  // std::cout << "checking..." << std::endl;
549  for(int nphi= 0; nphi < NPhiBins;nphi++){
550  // std::cout << "nphi: " << nphi << "|" << NPhiBins<< std::endl;
551  // std::cout << " adcval size " << adcval[nphi].size() << std::endl;
552  // std::cout << " outval size " << outval[nphi].size() << std::endl;
553  int bad = 0;
554  for(int nz= 0; nz < NZBins;nz++){
555  ncheck++;
556  if(adcval[nphi][nz]!=outval[nphi][nz]){
557  if(outval[nphi][nz]!=255){
558  bad++;
559  allbad++;
560  }
561  }
562  }
563 
564  if(bad>0){
565  std::cout << "sector: " << sector << " lay: " << layer << " # " << nphi << " bad:" << bad << " packsize: "<< (int)rhitset->m_tpchits[nphi].size() << std::endl;
566 
567  for(int nz= 0; nz < NZBins;nz++){
568  std::cout << "sector: " << sector << " lay: " << layer << " # " << nphi << "|"<< nz << " bad:" << bad <<" org: " << (int)adcval[nphi][nz] << " pack: "<< (int)outval[nphi][nz] << std::endl;
569  }
570  for(int nz= 0; nz < (int)rhitset->m_tpchits[nphi].size() ;nz++){
571  std::cout << "lay: " << layer << " # " << nphi << "|"<< nz << " bad:" << bad <<" pack: " << (int) rhitset->m_tpchits[nphi][nz] << std::endl;
572  }
573 
574  }
575 
576  }
577 
578  }
579 
580  //#ifdef WORK // if(layer==7)
581 
582 
583  }
584  std::cout << " number of hits checked " << ncheck << " bad: " << allbad << std::endl;
585 
586  return Fun4AllReturnCodes::EVENT_OK;
587 }
588 
589 int TpcRawWriter::End(PHCompositeNode */*topNode*/)
590 {
591  return Fun4AllReturnCodes::EVENT_OK;
592 }