ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
PHG4TpcPadBaselineShift.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file PHG4TpcPadBaselineShift.cc
2 
3 #include <tpc/TpcDefs.h>
4 
5 #include <trackbase/ActsSurfaceMaps.h> // for ActsSurfaceMaps
6 #include <trackbase/ActsTrackingGeometry.h> // for ActsTrackingG...
7 #include <trackbase/TrkrClusterContainer.h> // for TrkrClusterCo...
8 #include <trackbase/TrkrClusterHitAssoc.h> // for TrkrClusterHi...
9 #include <trackbase/TrkrHit.h> // for TrkrHit
10 
11 #include <trackbase/TrkrDefs.h> // for hitkey, getLayer
12 #include <trackbase/TrkrHitSet.h>
14 
16 #include <fun4all/SubsysReco.h> // for SubsysReco
17 
20 
21 #include <phool/PHCompositeNode.h>
22 #include <phool/PHNode.h> // for PHNode
23 #include <phool/PHNodeIterator.h>
24 #include <phool/getClass.h>
25 #include <phool/phool.h> // for PHWHERE
26 
27 #include <TF1.h>
28 #include <TFile.h>
29 #include <TTree.h>
30 
31 #include <cmath> // for sqrt, cos, sin
32 #include <iostream>
33 #include <map> // for _Rb_tree_cons...
34 #include <string>
35 #include <utility> // for pair
36 #include <vector>
37 
38 namespace
39 {
40  template <class T>
41  inline constexpr T square(const T &x)
42  {
43  return x * x;
44  }
45 } // namespace
46 
47 int findRBin(float R)
48 {
49  // Finding pad number from the center (bin) for hits
50  int binR = -1;
51  //Realistic binning
52  //double r_bins[r_bins_N+1] = {217.83,
53  // 311.05,317.92,323.31,329.27,334.63,340.59,345.95,351.91,357.27,363.23,368.59,374.55,379.91,385.87,391.23,397.19,402.49,
54  // 411.53,421.70,431.90,442.11,452.32,462.52,472.73,482.94,493.14,503.35,513.56,523.76,533.97,544.18,554.39,564.59,574.76,
55  // 583.67,594.59,605.57,616.54,627.51,638.48,649.45,660.42,671.39,682.36,693.33,704.30,715.27,726.24,737.21,748.18,759.11};
56  const int r_bins_N = 53;
57  double r_bins[r_bins_N + 1];
58  r_bins[0] = 30.3125;
59  double bin_width = 0.625;
60  for (int i = 1; i < r_bins_N; i++)
61  {
62  if (i == 16) bin_width = 0.9375;
63  if (i > 16) bin_width = 1.25;
64  if (i == 31) bin_width = 1.1562;
65  if (i > 31) bin_width = 1.0624;
66 
67  r_bins[i] = r_bins[i - 1] + bin_width;
68  }
69 
70  double R_min = 30;
71  while (R > R_min)
72  {
73  binR += 1;
74  R_min = r_bins[binR];
75  }
76  return binR;
77 }
78 
79 //____________________________________________________________________________..
81  : SubsysReco(name)
82 {
83  std::cout << "PHG4TpcPadBaselineShift::PHG4TpcPadBaselineShift(const std::string &name) Calling ctor" << std::endl;
84 }
85 
87 {
88  bool reject_it = false;
89 
90  // sector boundaries occur every 1/12 of the full phi bin range
91  int PhiBins = layergeom->get_phibins();
92  int PhiBinsSector = PhiBins / 12;
93 
94  double radius = layergeom->get_radius();
95  double PhiBinSize = 2.0 * radius * M_PI / (double) PhiBins;
96 
97  // sector starts where?
98  int sector_lo = sector * PhiBinsSector;
99  int sector_hi = sector_lo + PhiBinsSector - 1;
100 
101  int sector_fiducial_bins = (int) (SectorFiducialCut / PhiBinSize);
102 
103  if (phibin < sector_lo + sector_fiducial_bins || phibin > sector_hi - sector_fiducial_bins)
104  {
105  reject_it = true;
106  }
107 
108  return reject_it;
109 }
110 //____________________________________________________________________________..
112 {
113  std::cout << "PHG4TpcPadBaselineShift::~PHG4TpcPadBaselineShift() Calling dtor" << std::endl;
114 }
115 
116 //____________________________________________________________________________..
118 {
119  //outfile = new TFile(_filename.c_str(), "RECREATE");
120  _hit_z = 0;
121  _hit_r = 0;
122  _hit_phi = 0;
123  _hit_e = 0;
124  _hit_adc = 0;
125  _hit_adc_bls = 0;
126  _hit_layer = -1;
127  _hit_sector = -1;
128  if (_writeTree == 1)
129  {
130  outfile = new TFile(_filename.c_str(), "RECREATE");
131  _rawHits = new TTree("hTree", "tpc hit tree for base-line shift tests");
132  _rawHits->Branch("z", &_hit_z);
133  _rawHits->Branch("r", &_hit_r);
134  _rawHits->Branch("phi", &_hit_phi);
135  _rawHits->Branch("e", &_hit_e);
136  _rawHits->Branch("adc", &_hit_adc);
137  _rawHits->Branch("adc_BLS", &_hit_adc_bls);
138  _rawHits->Branch("hit_layer", &_hit_layer);
139  _rawHits->Branch("_hit_sector", &_hit_sector);
140  }
141  return 0;
142  //std::cout << "PHG4TpcPadBaselineShift::Init(PHCompositeNode *topNode) Initializing" << std::endl;
143  //return Fun4AllReturnCodes::EVENT_OK;
144 }
145 
146 //____________________________________________________________________________..
148 {
149  PHNodeIterator iter(topNode);
150 
151  // Looking for the DST node
152  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
153  if (!dstNode)
154  {
155  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
157  }
158 
159  std::cout << "PHG4TpcPadBaselineShift::InitRun(PHCompositeNode *topNode) Initializing for Run XXX" << std::endl;
161 }
162 
163 //____________________________________________________________________________..
165 {
166  // int print_layer = 18;
167 
168  if (Verbosity() > 1000)
169  {
170  std::cout << "PHG4TpcPadBaselineShift::Process_Event" << std::endl;
171  }
172 
173  PHNodeIterator iter(topNode);
174  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
175  if (!dstNode)
176  {
177  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
179  }
180 
181  // get node containing the digitized hits
182  m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
183  if (!m_hits)
184  {
185  std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
187  }
188 
189  // get node for clusters
190  m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
191  if (!m_clusterlist)
192  {
193  std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << std::endl;
195  }
196 
197  // get node for cluster hit associations
198  m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
199  if (!m_clusterhitassoc)
200  {
201  std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << std::endl;
203  }
204 
205  PHG4CylinderCellGeomContainer *geom_container =
206  findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
207  if (!geom_container)
208  {
209  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
211  }
212 
213  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
214  "ActsTrackingGeometry");
215  if (!m_tGeometry)
216  {
217  std::cout << PHWHERE
218  << "ActsTrackingGeometry not found on node tree. Exiting"
219  << std::endl;
221  }
222 
223  m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
224  "ActsSurfaceMaps");
225  if (!m_surfMaps)
226  {
227  std::cout << PHWHERE
228  << "ActsSurfaceMaps not found on node tree. Exiting"
229  << std::endl;
231  }
232 
233  // 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
234  // The TPC clustering is more complicated than for the silicon, because we have to deal with overlapping clusters
235 
236  // loop over the TPC HitSet objects
238  //const int num_hitsets = std::distance(hitsetrange.first,hitsetrange.second);
239 
240  //Loop over R positions & sectors
241  for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
242  hitsetitr != hitsetrange.second;
243  ++hitsetitr)
244  {
245  TrkrHitSet *hitset = hitsetitr->second;
246  unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
247  int side = TpcDefs::getSide(hitsetitr->first);
248  unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
249  PHG4CylinderCellGeom *layergeom = geom_container->GetLayerCellGeom(layer);
250 
251  _hit_sector = sector;
252  _hit_layer = layer;
253 
254  double radius = layergeom->get_radius(); //in cm
255 
256  _hit_r = radius;
257 
258  unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
259  unsigned short NPhiBinsSector = NPhiBins / 12;
260  unsigned short NZBins = (unsigned short) layergeom->get_zbins();
261  unsigned short NZBinsSide = NZBins / 2;
262  unsigned short NZBinsMin = 0;
263  unsigned short PhiOffset = NPhiBinsSector * sector;
264 
265  if (side == 0)
266  {
267  NZBinsMin = 0;
268  NZBinsMax = NZBins / 2 - 1;
269  }
270  else
271  {
272  NZBinsMin = NZBins / 2;
273  NZBinsMax = NZBins;
274  }
275  unsigned short ZOffset = NZBinsMin;
276  //Gives per pad ADC for particular R and sector in event
277  int perPadADC = 0;
278  unsigned short phibins = NPhiBinsSector;
279  unsigned short phioffset = PhiOffset;
280  unsigned short zbins = NZBinsSide;
281  unsigned short zoffset = ZOffset;
282  float sumADC = perPadADC; // 14 lines later just set to zero
283 
284  //phibins - number of pads in the sector
285  //The Sampa clock time is 18.8 MHz, so the sampling time is 53.2 ns = 1 Z bin.
286  //The Z bin range for one side of the TPC is 0-248.
287  //Check: 249 bins x 53.2 ns = 13.2 microseconds.
288  //The maximum drift time in the TPC is 13.2 microseconds.
289  //So, 53.2 ns / Z bin.
290 
291  TF1 *f1 = new TF1("f1", "[0]*exp(-(x-[1])/[2])", 0, 1000);
292  f1->SetParameter(0, 0.005);
293  f1->SetParameter(1, 0);
294  f1->SetParameter(2, 60); // in terms of 50nsec time bins
295 
296  // sumADC=0.; // this is set to zero 14 lines up (perPadADC is zero)
297  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
298 
299  std::vector<unsigned short> adcval(zbins, 0);
300  // std::multimap<unsigned short, ihit> all_hit_map;
301  // std::vector<ihit> hit_vect;
302  // Loop over phi & z
303  for (TrkrHitSet::ConstIterator hitr = hitrangei.first; hitr != hitrangei.second; ++hitr)
304  {
305  unsigned short phibin = TpcDefs::getPad(hitr->first) - phioffset;
306  unsigned short zbin = TpcDefs::getTBin(hitr->first) - zoffset;
307  float_t fadc = (hitr->second->getAdc()); // proper int rounding +0.5
308  unsigned short adc = 0;
309  if (fadc > 0)
310  {
311  adc = (unsigned short) fadc;
312  }
313  if (phibin >= phibins) continue;
314  if (zbin >= zbins) continue; // zbin is unsigned int, <0 cannot happen
315  adcval[zbin] = (unsigned short) adc;
316  sumADC += adc;
317  }
318  //Define ion-induced charge
319  sumADC /= phibins;
320  float ind_charge = -0.5 * sumADC * _CScale; //CScale is the coefficient related to the capacitance of the bottom layer of the bottom GEM
321  double pi = M_PI;
322 
323  for (TrkrHitSet::ConstIterator hitr = hitrangei.first; hitr != hitrangei.second; ++hitr)
324  {
325  unsigned short phibin = TpcDefs::getPad(hitr->first);
326  unsigned short zbin = TpcDefs::getTBin(hitr->first);
327  // Get the hitkey
329  TrkrHit *hit = nullptr;
330  hit = hitsetitr->second->getHit(hitkey);
331 
332  zbin = TpcDefs::getTBin(hitr->first);
333  phibin = TpcDefs::getPad(hitr->first);
334  double phi_center = layergeom->get_phicenter(phibin);
335  if (phi_center < 0) phi_center += 2 * pi;
336  _hit_phi = phi_center;
337  _hit_z = layergeom->get_zcenter(zbin);
338 
339  if (hit) _hit_e = hit->getEnergy();
340  _hit_adc = 0;
341  _hit_adc_bls = 0;
342 
343  float_t fadc = (hitr->second->getAdc()); // proper int rounding +0.5
344  _hit_adc = fadc;
345  _hit_adc_bls = _hit_adc + int(ind_charge);
346 
347  if (hit && _hit_adc > 0)
348  {
349  //Trkr hit has only one value m_adc which is energy and ADC at the same time
350  if (_hit_adc_bls > 0)
351  {
352  hit->setAdc(_hit_adc_bls);
353  }
354  else
355  {
356  hit->setAdc(0);
357  }
358  }
359  if (_writeTree == 1) _rawHits->Fill();
360  }
361 
362  //hitsetitr++;
363  }
364 
365  //pthread_attr_destroy(&attr);
366 
367  // wait for completion of all threads
368  //for( const auto& thread_pair:threads )
369  //{
370  // int rc2 = pthread_join(thread_pair.thread, nullptr);
371  // if (rc2)
372  // { std::cout << "Error:unable to join," << rc2 << std::endl; }
373  //}
374 
375  if (Verbosity() > 0)
376  std::cout << "TPC Clusterizer found " << m_clusterlist->size() << " Clusters " << std::endl;
377  std::cout << "PHG4TpcPadBaselineShift::process_event(PHCompositeNode *topNode) Processing Event" << std::endl;
379 }
380 
381 //____________________________________________________________________________..
382 //int PHG4TpcPadBaselineShift::ResetEvent(PHCompositeNode *topNode)
383 //{
384 // std::cout << "PHG4TpcPadBaselineShift::ResetEvent(PHCompositeNode *topNode) Resetting internal structures, prepare for next event" << std::endl;
385 // return Fun4AllReturnCodes::EVENT_OK;
386 //}
387 
388 //____________________________________________________________________________..
389 //int PHG4TpcPadBaselineShift::EndRun(const int runnumber)
390 //{
391 // std::cout << "PHG4TpcPadBaselineShift::EndRun(const int runnumber) Ending Run for Run " << runnumber << std::endl;
392 // return Fun4AllReturnCodes::EVENT_OK;
393 //}
394 
395 //____________________________________________________________________________..
397 {
398  if (_writeTree == 1)
399  {
400  outfile->cd();
401  outfile->Write();
402  outfile->Close();
403  std::cout << "PHG4TpcPadBaselineShift::End(PHCompositeNode *topNode) This is the End..." << std::endl;
404  }
405  //return 0;
407 }
408 
409 //____________________________________________________________________________..
410 //int PHG4TpcPadBaselineShift::Reset(PHCompositeNode *topNode)
411 //{
412 // std::cout << "PHG4TpcPadBaselineShift::Reset(PHCompositeNode *topNode) being Reset" << std::endl;
413 // return Fun4AllReturnCodes::EVENT_OK;
414 //}
415 
416 //____________________________________________________________________________..
417 //void PHG4TpcPadBaselineShift::Print(const std::string &what) const
418 //{
419 // std::cout << "PHG4TpcPadBaselineShift::Print(const std::string &what) const Printing info for " << what << std::endl;
420 //}
422 {
423  _CScale = CScale;
424  std::cout << "PHG4TpcPadBaselineShift::setFileName: Scale factor is set to:" << CScale << std::endl;
425 }
427 {
429  std::cout << "PHG4TpcPadBaselineShift::setFileName: Output file name for PHG4TpcPadBaselineShift is set to:" << filename << std::endl;
430 }
432 {
433  _writeTree = f_writeTree;
434  std::cout << "PHG4TpcPadBaselineShift::writeTree: True" << std::endl;
435 }