Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TpcRawDataDecoder.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TpcRawDataDecoder.cc
1 
2 #include "TpcRawDataDecoder.h"
3 
4 #include <trackbase/TpcDefs.h>
5 #include <trackbase/TrkrDefs.h>
6 #include <trackbase/TrkrHit.h>
7 #include <trackbase/TrkrHitSet.h>
9 
12 #include <trackbase/TrkrHitSetv1.h>
13 #include <trackbase/TrkrHitv2.h>
14 
15 #include <phool/PHCompositeNode.h>
16 #include <phool/PHIODataNode.h> // for PHIODataNode
17 #include <phool/PHNodeIterator.h> // for PHNodeIterator
18 #include <phool/PHObject.h> // for PHObject
19 #include <phool/getClass.h>
20 #include <phool/phool.h>
21 
26 
27 #include <cdbobjects/CDBTTree.h>
29 
30 #include <Event/Event.h>
31 #include <Event/EventTypes.h>
32 #include <Event/packet.h>
33 
35 
36 #include <TFile.h>
37 #include <TH2.h>
38 #include <TH3.h>
39 #include <TNtuple.h>
40 #include <limits.h>
41 
42 #include <string>
43 
44 //____________________________________________________________________________..
46  : SubsysReco(name)
47  , hm(nullptr)
48  , _filename("./outputfile.root")
49  {
50  std::cout << "TpcRawDataDecoder::TpcRawDataDecoder(const std::string &name)" << std::endl;
51  starting_BCO = -1;
52  rollover_value = 0;
53  current_BCOBIN = 0;
54  //M.setMapNames("AutoPad-R1-RevA.sch.ChannelMapping.csv", "AutoPad-R2-RevA-Pads.sch.ChannelMapping.csv", "AutoPad-R3-RevA.sch.ChannelMapping.csv");
55 
56  // Open a file, save the ntuple and close the file
57  //TFile in_file("/sphenix/user/shulga/Work/Files/pedestal-10616-outfile.root");
58  //in_file.GetObject("h_Alive",h_Alive);
59  //float chan_id,fee_id,module_id,pedMean,pedStdi, sec_id;
60  //float* row_content;
61 //
62  // if( Verbosity() )std::cout << "chan_id\t fee_id\t module_id\t pedMean\t pedStdi\t sec_id\n";
63  // for (int irow=0;irow<h_Alive->GetEntries();++irow)
64  // {
65  // h_Alive->GetEntry(irow);
66  // row_content = h_Alive->GetArgs();
67  // chan_id = row_content[0];
68  // fee_id = row_content[1];
69  // module_id = row_content[2];
70  // pedMean = row_content[3];
71  // pedStdi = row_content[4];
72  // sec_id = row_content[5];
73  // if( Verbosity() )
74  // {
75  // std::cout
76  // << chan_id << "\t"
77  // << fee_id << "\t"
78  // << module_id << "\t"
79  // << pedMean << "\t"
80  // << pedStdi << "\t"
81  // << sec_id << "\t"
82  // << std::endl;
83  // }
84 //
85  // struct ped_tpc_map x
86  // {
87  // };
88 //
89  // x.CHN_ID = chan_id ;
90  // x.FEE_ID = fee_id ;
91  // x.MOD_ID = module_id;
92  // x.PedMean = pedMean ;
93  // x.PedStdi = pedStdi ;
94  // x.SEC_ID = sec_id ;
95 //
96  // unsigned int key = 256 * (fee_id) + chan_id;
97  // tmap[key] = x;
98  // }
99 }
100 
101 //____________________________________________________________________________..
103 {
104  delete hm;
105  std::cout << "TpcRawDataDecoder::~TpcRawDataDecoder() Calling dtor" << std::endl;
106 }
107 
108 //____________________________________________________________________________..
110 {
111  std::cout << "TpcRawDataDecoder::Init(PHCompositeNode *topNode) Initializing" << std::endl;
112 
114  std::string calibdir = m_cdb->getUrl("TPC_FEE_CHANNEL_MAP");
115  //calibdir = "/sphenix/user/shulga/Work/TpcPadPlane_phi_coresoftware/macros/CDBTest/TPCChannelMap.root";
116  if (calibdir[0] == '/')
117  {
118  // use generic CDBTree to load
119  m_cdbttree = new CDBTTree(calibdir.c_str());
121  }
122  else
123  {
124  std::cout << "TpcRawDataDecoder::::InitRun No calibration file found" << std::endl;
125  exit(1);
126  }
127 
128  //m_cdbttree = new CDBTTree( "/sphenix/user/shulga/Work/TpcPadPlane_phi_coresoftware/macros/CDBTest/test.root" );
129 
130  if(m_Debug==1){
131  hm = new Fun4AllHistoManager("HITHIST");
132 
133  //_h_hit_XYT = new TH3F("_h_hit_XYT" ,"_h_hit_XYT;X, [mm];Y, [mm]; T [BCO]", 400, -800, 800, 400, -800, 800, 500, 128000000000,128050000000);
134  _h_hit_XY = new TH2F("_h_hit_XY" ,"_h_hit_XY;X, [mm];Y, [mm]", 400, -800, 800, 400, -800, 800);
135  _h_hit_XY_ADCcut = new TH2F("_h_hit_XY_ADCcut" ,"_h_hit_XY_ADCcut;X, [mm];Y, [mm]", 400, -800, 800, 400, -800, 800);
136  //_h_hit_PT_ADCcut = new TH3F("_h_hit_PT_ADCcut" ,"_h_hit_PT_ADCcut;Pad number;time [50 ns];BCO;", 400, -0.5, 399.5, 400, -0.5, 399.5, 500, 128000000000,128050000000);
137 
138  //hm->registerHisto(_h_hit_XYT );
139  //hm->registerHisto(_h_hit_PT_ADCcut);
142  }
143 
145 }
146 
147 //____________________________________________________________________________..
149 {
150 
151  // get dst node
152  PHNodeIterator iter(topNode);
153  auto dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
154  if (!dstNode)
155  {
156  std::cout << "TpcRawDataDecoder::InitRun - DST Node missing, doing nothing." << std::endl;
157  exit(1);
158  }
159 
160  // create hitset container if needed
161  auto hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
162  if (!hitsetcontainer)
163  {
164  std::cout << "TpcRawDataDecoder::InitRun - creating TRKR_HITSET." << std::endl;
165 
166  // find or create TRKR node
167  PHNodeIterator dstiter(dstNode);
168  auto trkrnode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
169  if (!trkrnode)
170  {
171  trkrnode = new PHCompositeNode("TRKR");
172  dstNode->addNode(trkrnode);
173  }
174 
175  // create container and add to the tree
176  hitsetcontainer = new TrkrHitSetContainerv1;
177  auto newNode = new PHIODataNode<PHObject>(hitsetcontainer, "TRKR_HITSET", "PHObject");
178  trkrnode->addNode(newNode);
179  }
180  topNode->print();
181 
182  // we reset the BCO for the new run
183  starting_BCO = -1;
184  rollover_value = 0;
185  current_BCOBIN = 0;
186 
187  //m_hits = new TrkrHitSetContainerv1();
188 
190 }
191 
192 //____________________________________________________________________________..
194 {
195  _ievent++;
196  //if(_ievent<1270 || _ievent>1270+200) return Fun4AllReturnCodes::DISCARDEVENT;
197  // load relevant nodes
198  // Get the TrkrHitSetContainer node
199  auto trkrhitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
200  assert(trkrhitsetcontainer);
201 
202  PHG4TpcCylinderGeomContainer *geom_container =
203  findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
204  if (!geom_container)
205  {
206  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
208  }
209 
210 
211  Event *_event = findNode::getClass<Event>(topNode, "PRDF");
212  assert( _event );
213 
214  if (_event == nullptr)
215  {
216  std::cout << "TpcRawDataDecoder::Process_Event - Event not found" << std::endl;
217  return -1;
218  }
219  if (_event->getEvtType() >= 8)
220  {
222  }
223  rollover_value = 0;
224  // check all possible TPC packets that we need to analyze
225  for(int ep=0;ep<2;ep++){
226  for (int sector = 0; sector<24; sector++)
227  {
228  const int packet = 4000 + sector*10 + ep;
229 
230  // Reading packet
231  Packet *p = _event->getPacket(packet);
232 
233  // Figure out which side
234  int side = 0;
235  if(sector>11) side=1;
236  if(m_Debug==1){
237  char buff[100];
238  snprintf(buff, sizeof(buff), "./outputfile_%i.root",sector );
239  if(p) _filename = buff;
240  }
241  if (p)
242  {
243  std::cout << "My Own TpcRawDataDecoder:: Event getPacket: "<< packet << "| Sector"<< sector << "| EndPoint "<< ep << std::endl;
244  }else{
245  continue;
246  }
247 
248  uint64_t triggerBCO = 0;
249  int n_tagger = p->lValue(0, "N_TAGGER");
250  for (int t = 0; t < n_tagger; t++){
251  // const auto tagger_type = static_cast<uint16_t>(p->lValue(t, "TAGGER_TYPE"));
252  // const auto is_endat = static_cast<uint8_t>(p->lValue(t, "IS_ENDAT"));
253  const auto is_lvl1 = static_cast<uint8_t>(p->lValue(t, "IS_LEVEL1_TRIGGER"));
254  const auto bco = static_cast<uint64_t>(p->lValue(t, "BCO"));
255  // const auto lvl1_count = static_cast<uint32_t>(p->lValue(t, "LEVEL1_COUNT"));
256  // const auto endat_count = static_cast<uint32_t>(p->lValue(t, "ENDAT_COUNT"));
257  // const auto last_bco = static_cast<uint64_t>(p->lValue(t, "LAST_BCO"));
258  // const auto modebits = static_cast<uint8_t>(p->lValue(t, "MODEBITS"));
259 
260  // only printout the is_lvl1 triggers
261  // if( Verbosity() )
262  if( is_lvl1 )
263  {
264  triggerBCO = bco;
265  std::cout << " is_lvl1: " << (bool) (is_lvl1)
266  << " bco: " << bco
267  << "_ievent: " << _ievent
268  << std::endl;
269  }
270  }
271  //if (triggerBCO != 128330850912) return Fun4AllReturnCodes::DISCARDEVENT;
272  int nr_of_waveforms = p->iValue(0, "NR_WF");
273 
274  // for (auto &l : m_hitset){
275  // l = new TrkrHitSetv1();
276 
277  int wf;
278  //find waveform with earliest BCO in packet
279  //take earliest BCO as starting BCO - works only for very long waveform setting
280  //current waveforms are ~360 timebins long, eariest BCO == starting BCO should hold for these conditions
281  //Check if BCO of previous event is within 360 time bins or 180 BCO. If yes we we are looking at a truncated event and we use the BCO of the previous one as startingBCO
282 
283  int earliest_BCO = INT_MAX;
284 
285  for (wf = 0; wf < nr_of_waveforms; wf++){
286  int current_BCO = p->iValue(wf, "BCO");
287  if(current_BCO < earliest_BCO) earliest_BCO = current_BCO;
288  // std::cout << " earliest BCO: " << earliest_BCO << " current BCO " << current_BCO << std::endl;
289  }
290  // if((earliest_BCO - starting_BCO > 0) && (earliest_BCO - starting_BCO < 180)){
291  // starting_BCO = starting_BCO;
292  //}else{
293  starting_BCO = earliest_BCO;
294  //}
295  if( Verbosity() )
296  std::cout << " _ievent: " << _ievent << " earliest BCO: " << earliest_BCO << " ep: " << ep << " sector " << sector << " trigBCO: " << triggerBCO << " diff " << triggerBCO - earliest_BCO << std::endl;
297 
298  for (wf = 0; wf < nr_of_waveforms; wf++){
299  int current_BCO = p->iValue(wf, "BCO") + rollover_value;
300 
301  /*
302  std::cout << " BCO: " << p->iValue(wf, "BCO") << std::endl;
303  std::cout << " rollover: " << rollover_value << std::endl;
304  std::cout << " current BCO: " << current_BCO << std::endl;
305  std::cout << " starting BCO: " << starting_BCO << std::endl;
306  std::cout << " earliest BCO: " << earliest_BCO << std::endl;
307  std::cout << " current BCOBIN: " << current_BCOBIN << std::endl;
308  */
309  if (starting_BCO < 0)
310  {
311  starting_BCO = current_BCO;
312  }
313 
314  if (current_BCO < starting_BCO) // we have a rollover
315  {
316  rollover_value = 0x100000;
317  current_BCO = p->iValue(wf, "BCO") + rollover_value;
318  starting_BCO = current_BCO;
319  current_BCOBIN++;
320  }
321  // std::cout << " nu current BCO: " << current_BCO << std::endl;
322  //std::cout << " nu starting BCO: " << starting_BCO << std::endl;
323  int sampa_nr = p->iValue(wf, "SAMPAADDRESS");
324  int channel = p->iValue(wf, "CHANNEL");
325 
326  int fee = p->iValue(wf, "FEE");
327  int samples = p->iValue( wf, "SAMPLES" );
328  // clockwise FEE mapping
329  //int FEE_map[26]={5, 6, 1, 3, 2, 12, 10, 11, 9, 8, 7, 1, 2, 4, 8, 7, 6, 5, 4, 3, 1, 3, 2, 4, 6, 5};
330  int FEE_R[26]={2, 2, 1, 1, 1, 3, 3, 3, 3, 3, 3, 2, 2, 1, 2, 2, 1, 1, 2, 2, 3, 3, 3, 3, 3, 3};
331  // conter clockwise FEE mapping (From Takao)
332  //int FEE_map[26]={3, 2, 5, 3, 4, 0, 2, 1, 3, 4, 5, 7, 6, 2, 0, 1, 0, 1, 4, 5, 11, 9, 10, 8, 6, 7};
333  //int pads_per_sector[3] = {96, 128, 192};
334  int FEE_map[26] = { 4, 5, 0, 2, 1, 11, 9, 10, 8, 7, 6, 0, 1, 3, 7, 6, 5, 4, 3, 2, 0, 2, 1, 3, 5, 4};
335  // setting the mapp of the FEE
336  int feeM = FEE_map[fee];
337  if(FEE_R[fee]==2) feeM += 6;
338  if(FEE_R[fee]==3) feeM += 14;
339  unsigned int key = 256 * (feeM) + channel;
340  //int layer = M.getLayer(feeM, channel);
341 
342  std::string varname = "layer";// + std::to_string(key);
343  int layer = m_cdbttree->GetIntValue(key,varname);
344  // antenna pads will be in 0 layer
345  if(layer==0)continue;
346 
347  PHG4TpcCylinderGeom *layergeom = geom_container->GetLayerCellGeom(layer);
348  //varname = "fee" + std::to_string(key);
349  //int feeM_cdbt = m_cdbttree->GetIntValue(key,varname);
350 
351  //varname = "channel" + std::to_string(key);
352  //int feeM_chn_cdbt = m_cdbttree->GetIntValue(key,varname);
353 
354  varname = "R";// + std::to_string(key);
355  double R = m_cdbttree->GetDoubleValue(key,varname);
356 
357  varname = "phi";// + std::to_string(key);
358  double phi = pow(-1,side)*m_cdbttree->GetDoubleValue(key,varname) + (sector - side*12 )* M_PI / 6 ;;
359 
360  //varname = "padN" + std::to_string(key);
361  //int padn_cdbt = m_cdbttree->GetIntValue(key,varname);
362 
363  //double R = M.getR(feeM, channel);
364  //double phi = M.getPhi(feeM, channel) + (sector - side*12 )* M_PI / 6 ;
365  //int pad = M.getPad(feeM, channel);
366  //varname = "pad";// + std::to_string(key);
367  //int pad = m_cdbttree->GetIntValue(key,varname);
368 
369  //std::cout<<"pad "<< padn_cdbt <<" "<< pad << " delta =" << padn_cdbt - pad
370  //<<"\n \t phi "<< phi_cdbt <<" "<< phi << " delta =" << phi_cdbt - phi
371  //<<"\n \t R "<< R_cdbt <<" "<< R << " delta =" << R_cdbt - R
372  //<<"\n \t FEE chn"<< feeM_chn_cdbt <<" "<< channel << " delta =" << feeM_chn_cdbt - channel
373  //<<"\n \t FEE "<< feeM_cdbt <<" "<< feeM << " delta =" << feeM_cdbt - feeM
374  //<<"\n \t layer "<< layer_cdbt <<" "<< layer << " delta =" << layer_cdbt - layer
375  //<< std::endl;
376  int mc_sectors[12] = { 5, 4, 3, 2, 1, 0, 11, 10, 9, 8, 7, 6};
377  //int mc_sectors[12] = {5, 6, 7, 8, 9, 10, 11, 0, 1, 2, 3, 4};
378 
379  float pedestal = 72.4;//round(tmap[key].PedMean);
380  TrkrDefs::hitsetkey tpcHitSetKey = TpcDefs::genHitSetKey(layer, (mc_sectors[sector - side*12] ), side);
381  TrkrHitSetContainer::Iterator hitsetit = trkrhitsetcontainer->findOrAddHitSet(tpcHitSetKey);
382 
383 
384  if( Verbosity()>1 ){
385  int sampaAddress = p->iValue(wf, "SAMPAADDRESS");
386  int sampaChannel = p->iValue(wf, "SAMPACHANNEL");
387  int checksum = p->iValue(wf, "CHECKSUM");
388  int checksumError = p->iValue(wf, "CHECKSUMERROR");
389  std::cout << "TpcRawDataDecoder::Process_Event Samples "<< samples
390  <<" Chn:"<< channel
391  <<" layer: " << layer
392  << " sampa: "<< sampa_nr
393  << " fee: "<< fee
394  << " Mapped fee: "<< feeM
395  << " sampaAddress: "<< sampaAddress
396  << " sampaChannel: "<< sampaChannel
397  << " checksum: "<< checksum
398  << " checksumError: "<< checksumError
399  << " hitsetkey "<< tpcHitSetKey
400  << " R = " << R
401  << " phi = " << phi
402  << std::endl;
403 
404  }
405  pedestal = 0.0;
406  for (int s = 0; s < 5; s++)
407  {
408  int adc = p->iValue(wf,s);
409  pedestal += adc;
410  }
411  pedestal /=5;
412  for (int s = 0; s < samples; s++)
413  {
414 
415  int t = s; // + 2 * (current_BCO - starting_BCO);
416  int adc = p->iValue(wf,s);
417  if(m_Debug==1){
418  _h_hit_XY->Fill(R*cos(phi),R*sin(phi),float(adc)-pedestal);
419  }
420  // if(adc-pedestal<4) continue;
421  // generate hit key
422  if(float(adc)-pedestal>2){
423  unsigned int phibin = layergeom->find_phibin(phi);
424  //if((int)phibin > pads_per_sector[FEE_R[fee]-1]*12) std::cout << "TpcRawDataDecoder:: phibin is out of range > "<< pads_per_sector[FEE_R[fee]-1]*12 << std::endl;
425  TrkrDefs::hitkey hitkey = TpcDefs::genHitKey( phibin, (unsigned int) t);
426  //double phi_center = layergeom->get_phicenter(phibin);
427  //if(phi_center<0) phi_center += 2*M_PI;
428  //int phibin_mc = layergeom->find_phibin(phi);
429  // find existing hit, or create
430  auto hit = hitsetit->second->getHit(hitkey);
431 
432  // create hit, assign adc and insert in hitset
433  if (!hit)
434  {
435 
436  // create a new one
437  hit = new TrkrHitv2();
438  hit->setAdc(float(adc)-pedestal);
439  //std::cout<< "ADC = " << adc << " Pedestal = "<< pedestal << "delta = "<< adc-pedestal << std::endl;
440  //if(adc - pedestal > 40) std::cout<< adc - pedestal << "| ";
441  //std::cout<< t << "| ";
442 
443  hitsetit->second->addHitSpecificKey(hitkey, hit);
444  }
445 
446 
447 
448  if(m_Debug==1){
449  //if(layer==16 ) _h_hit_PT_ADCcut->Fill(pad, s, triggerBCO,float(adc));
450  if(adc - pedestal > 15){
451  //std::cout << "pad = " << pad << "phibin = " << phibin << " phibin_mc = " << phibin_mc << " sector = " << sector << "mc_sectors = " << mc_sectors[sector - side*12] << "phibin/pads_per_sector = " << phibin/pads_per_sector[FEE_R[fee]-1] << " phi_center = " << phi_center << " phi = " << phi << std::endl;
452  _h_hit_XY_ADCcut->Fill(R*cos(phi),R*sin(phi),float(adc)-pedestal);
453  //_h_hit_XYT->Fill(R*cos(phi),R*sin(phi), triggerBCO,float(adc)-pedestal);
454 
455  }
456  }
457  }
458  //if (s==samples-1) std::cout << std::endl;
459  }
460 
461  }
462  // }//auto l:
463 
464  }// End of run over packets
465  }//End of ep loop
466 
467  // we skip the mapping to real pads at first. We just say
468  // that we get 16 rows (segment R2) with 128 pads
469  // so each FEE fills 2 rows. Not right, but one step at a time.
470  std::cout << "TpcRawDataDecoder:: done" << std::endl;
472 }
473 
474 //____________________________________________________________________________..
475 //int TpcRawDataDecoder::ResetEvent(PHCompositeNode * /*topNode*/)
476 //{
477 // std::cout << "TpcRawDataDecoder::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
478 // return Fun4AllReturnCodes::EVENT_OK;
479 //}
480 
481 //____________________________________________________________________________..
482 //int TpcRawDataDecoder::EndRun(const int runnumber)
483 //{
484 // std::cout << "TpcRawDataDecoder::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
485 // return Fun4AllReturnCodes::EVENT_OK;
486 //}
487 
488 //____________________________________________________________________________..
490 {
491  std::cout << "TpcRawDataDecoder::End(PHCompositeNode *topNode) This is the End..." << std::endl;
492  if(m_Debug==1)hm->dumpHistos(_filename, "RECREATE");
493 
495 }
496 
497 //____________________________________________________________________________..
498 //int TpcRawDataDecoder::Reset(PHCompositeNode * /*topNode*/)
499 //{
500 // std::cout << "TpcRawDataDecoder::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
501 // return Fun4AllReturnCodes::EVENT_OK;
502 //}
503 
504 //____________________________________________________________________________..
505 //void TpcRawDataDecoder::Print(const std::string &what) const
506 //{
507 // std::cout << "TpcRawDataDecoder::Print(const std::string &what) const Printing info for " << what << std::endl;
508 //}
509 
510 //____________________________________________________________________________..
511 //void TpcRawDataDecoder::setHistoFileName(const std::string &what)
512 //{
513 // _filename = what;
514 // std::cout << "TpcRawDataDecoder::setHistoFileName(const std::string &what) Histogram File Name is " << what << std::endl;
515 //}