Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TpcClusterBuilder.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TpcClusterBuilder.cc
1 #include "TpcClusterBuilder.h"
2 /* #include <trackbase/TrkrClusterv3.h> */
3 #include <trackbase/ClusHitsVerbosev1.h>
4 #include <trackbase/TrkrClusterv4.h>
5 #include <trackbase/TpcDefs.h>
6 #include <g4detectors/PHG4TpcCylinderGeom.h>
7 #include <g4tracking/TrkrTruthTrack.h>
8 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
9 #include <algorithm>
10 
11 #include <trackbase/TrkrHitSet.h>
12 #include <trackbase/TrkrHitSetContainerv1.h>
13 #include <trackbase/TrkrHitv2.h> // for TrkrHit
14 #include <cmath> // for sqrt, cos, sin
15 #include <map> // for _Rb_tree_cons...
16 #include <TString.h>
17 
18 #include <trackbase/TrkrClusterContainer.h>
19 #include <g4tracking/TrkrTruthTrackContainer.h>
20 
21 #include <set>
22 
23 #include <iostream>
24 #include <ios>
25 
26 /* class TpcClusterBuilder; */
27 
28 using std::cout, std::endl, std::string, std::ofstream, std::ostream;
29 
30 double TpcClusterBuilder::square(double v) { return v*v; }
31 double TpcClusterBuilder::square(float v) { return v*v; }
32 
33 void TpcClusterBuilder::set_verbosity(int verbosity_level) { verbosity = verbosity_level; }
34 
35 void TpcClusterBuilder::cluster_hits(TrkrTruthTrack* track) {
36  // now build the TrkrCluster
37  TrkrHitSetContainer::ConstRange hitsetrange;
38  hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
39 
40  //-----------------------------------------------------
41  // from TpcClusterizer.cc ~line: 837
42  //-----------------------------------------------------
43 
44  for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
45  hitsetitr != hitsetrange.second;
46  ++hitsetitr)
47  {
48  //if (count>0)continue;
49  TrkrDefs::hitsetkey hitsetkey = hitsetitr->first;
50  TrkrHitSet *hitset = hitsetitr->second;
51  unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
52  int side = TpcDefs::getSide(hitsetitr->first);
53  unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
54  PHG4TpcCylinderGeom *layergeom = geom_container->GetLayerCellGeom(layer);
55 
56  //get the maximum and minimum phi and time
57  unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
58  unsigned short NPhiBinsSector = NPhiBins/12;
59  unsigned short NTBins = (unsigned short)layergeom->get_zbins();
60  unsigned short NTBinsSide = NTBins;
61  unsigned short NTBinsMin = 0;
62  unsigned short PhiOffset = NPhiBinsSector * sector;
63  unsigned short TOffset = NTBinsMin;
64 
65  double m_tdriftmax = AdcClockPeriod * NTBins / 2.0;
66 
67  unsigned short phibins = NPhiBinsSector;
68  unsigned short phioffset = PhiOffset;
69  unsigned short tbins = NTBinsSide;
70  unsigned short toffset = TOffset ;
71 
72  // loop over the hits in this cluster
73  double t_sum = 0.0;
74  double adc_sum = 0.0;
75 
76  /* double phi_sum = 0.0; */
77  double iphi_sum = 0.0;
78  // double t2_sum = 0.0;
79  // double phi2_sum = 0.0;
80 
81  double radius = layergeom->get_radius(); // returns center of layer
82 
83  int phibinhi = -1;
84  int phibinlo = INT_MAX;
85  int tbinhi = -1;
86  int tbinlo = INT_MAX;
87 
88  auto ihit_list = hitset->getHits();
89 
90  const int iphi_max = phioffset+phibins;
91  const int it_max = toffset+tbins;
92 
93  double sum_adc { 0. }; // energy = 4 * adc at this point
94  //
95  // accumulate energy from all hits that are not out of range
96  for(auto iter = ihit_list.first; iter != ihit_list.second; ++iter){
97  int iphi = TpcDefs::getPad(iter->first) ;//- phioffset;
98  int it = TpcDefs::getTBin(iter->first);// - toffset;
99  if (iphi < phioffset || iphi > iphi_max || it < toffset || it > it_max) continue;
100  sum_adc += iter->second->getAdc();
101  }
102  const double threshold = sum_adc * m_pixel_thresholdrat;
103 
104  // FIXME -- see why the hits are so scattered
105  std::set<int> v_iphi, v_it; // FIXME
106  std::map<int,unsigned int> m_iphi, m_it, m_iphiCut, m_itCut; // FIXME
107  for(auto iter = ihit_list.first; iter != ihit_list.second; ++iter)
108  {
109  unsigned int adc = iter->second->getAdc();
110  if (adc <= 0) continue;
111  int iphi = TpcDefs::getPad(iter->first) ;//- phioffset;
112  int it = TpcDefs::getTBin(iter->first);// - toffset;
113  if (iphi < phioffset || iphi > iphi_max) {
114  std::cout << "WARNING phibin out of range: " << iphi << " | " << phibins << std::endl;
115  continue;
116  }
117  if (it < toffset || it > it_max) {
118  std::cout << "WARNING z bin out of range: " << it << " | " << tbins << std::endl;
119  continue;
120  }
121  if (adc<threshold) {
122  if (mClusHitsVerbose) {
123  auto pnew = m_iphiCut.try_emplace(iphi,adc);
124  if (!pnew.second) pnew.first->second += adc;
125 
126  pnew = m_itCut.try_emplace(it,adc);
127  if (!pnew.second) pnew.first->second += adc;
128  }
129  continue;
130  }
131 
132  v_iphi.insert(iphi);
133  v_it.insert(it);
134 
135  // fill m_iphi
136  auto pnew = m_iphi.try_emplace(iphi,adc);
137  if (!pnew.second) pnew.first->second += adc;
138  pnew = m_it.try_emplace(it,adc);
139  if (!pnew.second) pnew.first->second += adc;
140 
141  if (iphi > phibinhi) phibinhi = iphi;
142  if (iphi < phibinlo) phibinlo = iphi;
143  if (it > tbinhi) tbinhi = it;
144  if (it < tbinlo) tbinlo = it;
145 
146  iphi_sum += iphi * adc;
147  // phi2_sum += square(phi_center)*adc;
148 
149  // update t sums
150  double t = layergeom->get_zcenter(it);
151  t_sum += t*adc;
152  // t2_sum += square(t)*adc;
153 
154  adc_sum += adc;
155  }
156  if (mClusHitsVerbose) {
157  if (verbosity>10) {
158  for (auto& hit : m_iphi) {
159  cout << " m_phi(" << hit.first <<" : " << hit.second<<") ";
160  }
161  }
162  /* cout << " MELON 0 m_iphi "; */
163  /* for (auto& hit : m_iphi) cout << hit.first << "|" << hit.second << " "; */
164  /* cout << endl; */
165  /* cout << " MELON 1 m_it "; */
166  /* for (auto& hit : m_it) cout << hit.first << "|" << hit.second << " "; */
167  /* cout << endl; */
168  for (auto& hit : m_iphi) mClusHitsVerbose->addPhiHit (hit.first, hit.second);
169  for (auto& hit : m_it) mClusHitsVerbose->addZHit (hit.first, hit.second);
170  for (auto& hit : m_iphiCut) mClusHitsVerbose->addPhiCutHit (hit.first, hit.second);
171  for (auto& hit : m_itCut) mClusHitsVerbose->addZCutHit (hit.first, hit.second);
172  }
173 
174  // This is the global position
175  double clusiphi = iphi_sum / adc_sum;
176  double clusphi = layergeom->get_phi(clusiphi);
177  /* double clusphi = phi_sum / adc_sum; */
178  float clusx = radius * cos(clusphi);
179  float clusy = radius * sin(clusphi);
180  double clust = t_sum / adc_sum;
181  double zdriftlength = clust * m_tGeometry->get_drift_velocity();
182  // convert z drift length to z position in the TPC
183  double clusz = m_tdriftmax * m_tGeometry->get_drift_velocity() - zdriftlength;
184  if (side == 0) clusz = -clusz;
185 
186  char tsize = tbinhi - tbinlo + 1;
187  char phisize = phibinhi - phibinlo + 1;
188 
189  if (tsize < 0 && verbosity > 1) std::cout << " FIXME z4 tsize: " << ((int)tsize) << " " << tbinlo << " to " << tbinhi << std::endl;
190 
191  // -------------------------------------------------
192  // -------------------------------------------------
193  // debug here: FIXME
194  //
195  /* cout << " FIXME phisize " << ((int) phisize) << endl; */
196  //FIXME
197  if (false) { // Printing for debugging
198  if ((int)phisize > 10 || (int)tsize > 8) {
199  int _size_phi = ((int)phisize);
200  int _nbins_phi = v_iphi.size();
201  int _delta_phi = abs(_size_phi-_nbins_phi);
202  int _size_z = ((int)tsize);
203  int _nbins_z = v_it.size();
204  int _delta_z = abs(_size_z - _nbins_z);
205 
206  TString fmt;
207  cout << " x|"<<_delta_phi<<"|"<<_delta_z
208  <<"| new node FIXME A1 layer("<< layer
209  <<") (nset:size) phi("
210  << _nbins_phi<<":"<<_size_phi << ") z("
211  <<_nbins_z<<":"<<_size_z<<") "
212  <<"trkId("<<track->getTrackid()<<") trkpT("<<track->getPt()<<")"
213  << endl;
214  if (phisize>10) {
215  cout << " iphi-from-(";
216  int _prev = -1;
217  double tempsum = 0.;
218  for (auto _ : v_iphi) {
219  if (_prev == -1) {
220  cout << _<<"): ";
221  } else {
222  int _diff = ((int)_-_prev-1);
223  if (_diff != 0) cout<<">"<<_diff<<">";
224  }
225  /* cout << std::setprecision(2) << std::fixed; */
226  double _rat = (float)m_iphi[_] / (float)adc_sum;
227  fmt.Form("%.2f",_rat);
228  cout << fmt <<" ";
229  tempsum += _rat;
230  _prev = _;
231  }
232  if (tempsum < 0.999) cout << " Z3 sumphirat: " << tempsum;
233  cout << endl;
234  }
235  if (tsize>8) {
236  int _prev = -1;
237  double tempsum = 0.;
238  cout << " iz-from-(";
239  for (auto _ : v_it) {
240  if (_prev == -1) {
241  cout << _<<"): ";
242  } else {
243  int _diff = ((int)_-_prev-1);
244  if (_diff != 0) cout<<">"<<_diff<<">";
245  }
246  /* cout << std::setprecision(2) << std::fixed; */
247  double _rat = (float)m_it[_] / (float)adc_sum;
248  fmt.Form("%.2f",_rat);
249  cout << fmt <<" ";
250  tempsum += _rat;
251  _prev = _;
252  }
253  if (tempsum < 0.999) cout << " Z3 sumzrat: " << tempsum;
254  }
255  cout << endl;
256  }
257  } // end debug printing
258 
259  // get the global vector3 to then get the surface local phi and z
260  Acts::Vector3 global(clusx, clusy, clusz);
261  TrkrDefs::subsurfkey subsurfkey = 0;
262 
263  Surface surface = m_tGeometry->get_tpc_surface_from_coords(
264  hitsetkey, global, subsurfkey);
265 
266  if (!surface) {
267  if (verbosity) std::cout << "Can't find the surface! with hitsetkey " << ((int)hitsetkey) << std::endl;
268  /* if (mClusHitsVerbose) mClusHitsVerbose->Reset(); */
269  continue;
270  }
271 
272  // SAMPA shaping bias correction
273  clust = clust + m_sampa_tbias;
274 
275  global *= Acts::UnitConstants::cm;
276 
277  Acts::Vector3 local=surface->transform(m_tGeometry->geometry().getGeoContext()).inverse() * global;
278  local /= Acts::UnitConstants::cm;
279 
280  auto cluster = new TrkrClusterv4; //
281  cluster->setAdc(adc_sum);
282  /* cluster->setOverlap(ntouch); */
283  /* cluster->setEdge(nedge); */
284  cluster->setPhiSize(phisize);
285  cluster->setZSize(tsize);
286  cluster->setSubSurfKey(subsurfkey);
287  cluster->setLocalX(local(0));
288  cluster->setLocalY(clust);
289 
290  // add the clusterkey
291  auto empl = hitsetkey_cnt.try_emplace(hitsetkey,0);
292  if (!empl.second) empl.first->second += 1;
293  TrkrDefs::cluskey cluskey = TrkrDefs::genClusKey(hitsetkey, hitsetkey_cnt[hitsetkey]);
294  m_clusterlist->addClusterSpecifyKey(cluskey, cluster);
295 
296  track->addCluster(cluskey);
297  if (mClusHitsVerbose) {
298  mClusHitsVerbose->push_hits(cluskey);
299  /* cout << " FIXME z1 ClusHitsVerbose.size: " << mClusHitsVerbose->getMap().size() << endl; */
300  }
301  }
302  m_hits->Reset();
303 }
304 
305 /* void TpcClusterBuilder::set_current_track(TrkrTruthTrack* track) { */
306  /* current_track = track; */
307  /* ++ n_tracks; */
308 /* } */
309 
310 void TpcClusterBuilder::addhitset(
311  TrkrDefs::hitsetkey hitsetkey,
312  TrkrDefs::hitkey hitkey,
313  float neffelectrons)
314 {
315  // copy of code in PHG4TpcPadPlaneReadout::MapToPadPlane, with a switch
316  // to ignore non embedded tracks
317  if (!b_collect_hits) return;
318 
319  // Add the hitset to the current embedded track
320  // Code from PHG4TpcPadPlaneReadout::MapToPadPlane (around lines {}.cc::386-401)
321  TrkrHitSetContainer::Iterator hitsetit = m_hits->findOrAddHitSet(hitsetkey);
322  // See if this hit already exists
323  TrkrHit *hit = nullptr;
324  hit = hitsetit->second->getHit(hitkey);
325  if (!hit)
326  {
327  // create a new one
328  hit = new TrkrHitv2();
329  hitsetit->second->addHitSpecificKey(hitkey, hit);
330  }
331  // Either way, add the energy to it -- adc values will be added at digitization
332  hit->addEnergy(neffelectrons);
333 }
334 
335 void TpcClusterBuilder::clear_hitsetkey_cnt() {
336  hitsetkey_cnt.clear();
337 }
338 
339 void TpcClusterBuilder::print(
340  TrkrTruthTrackContainer* truth_tracks, int nclusprint) {
341  cout << " ------------- content of TrkrTruthTrackContainer ---------- " << endl;
342  auto& tmap = truth_tracks->getMap();
343  cout << " Number of tracks: xyz db : " << tmap.size() << endl;
344  for (auto& _pair : tmap) {
345  auto& track = _pair.second;
346 
347  printf("id(%2i) phi:eta:pt(", (int)track->getTrackid());
348  cout << "phi:eta:pt(";
349  printf("%5.2f:%5.2f:%5.2f", track->getPhi(), track->getPseudoRapidity(), track->getPt());
350  /* Form("%5.2:%5.2:%5.2", track->getPhi(), track->getPseudoRapidity(), track->getPt()) */
351  //<<track->getPhi()<<":"<<track->getPseudoRapidity()<<":"<<track->getPt()
352  cout << ") nclusters(" << track->getClusters().size() <<") ";
353  if (verbosity <= 10) { cout << endl; }
354  else {
355  int nclus = 0;
356  for (auto cluskey : track->getClusters()) {
357  cout << " "
358  << ((int) TrkrDefs::getHitSetKeyFromClusKey(cluskey)) <<":index(" <<
359  ((int) TrkrDefs::getClusIndex(cluskey)) << ")";
360  ++nclus;
361  if (nclusprint > 0 && nclus >= nclusprint) {
362  cout << " ... ";
363  break;
364  }
365  }
366  }
367  }
368  cout << " ----- end of tracks in TrkrrTruthTrackContainer ------ " << endl;
369 }
370 
371 void TpcClusterBuilder::print_file(
372  TrkrTruthTrackContainer* truth_tracks, string ofile_name)
373 {
374  ofstream fout;
375  fout.open(ofile_name.c_str());
376  fout << " ------------- content of TrkrTruthTrackContainer ---------- " << endl;
377  auto& tmap = truth_tracks->getMap();
378  fout << " Number of tracks: " << tmap.size() << endl;
379  for (auto& _pair : tmap) {
380  auto &track = _pair.second;
381  fout << " id( " << track->getTrackid() << ") phi:eta:pt("<<
382  track->getPhi()<<":"<<track->getPseudoRapidity()<<":"<<track->getPt() << ") nclusters("
383  << track->getClusters().size() <<") ";
384  int nclus = 0;
385  for (auto cluskey : track->getClusters()) {
386  auto C = m_clusterlist->findCluster(cluskey);
387  fout << " "
388  << ((int) TrkrDefs::getHitSetKeyFromClusKey(cluskey)) <<":" <<
389  ((int) TrkrDefs::getClusIndex(cluskey)) << "->adc:X:phisize:Y:zsize("
390  << C->getAdc() <<":"
391  << C->getLocalX() <<":"
392  << C->getPhiSize() <<":"
393  << C->getLocalY() <<":"
394  << C->getZSize() <<") ";
395  ++nclus;
396  }
397  fout << endl;
398  }
399  fout << " ----- end of tracks in TrkrrTruthTrackContainer ------ " << endl;
400  fout.close();
401 }
402 
403 void TpcClusterBuilder::set_input_nodes(
404  TrkrClusterContainer* _truth_cluster_container
405  , ActsGeometry* _ActsGeometry
406  , PHG4TpcCylinderGeomContainer* _geom_container
407  , ClusHitsVerbosev1* _clushitsverbose
408 ) {
409  m_clusterlist = _truth_cluster_container;
410  m_tGeometry = _ActsGeometry;
411  geom_container = _geom_container;
412  mClusHitsVerbose = _clushitsverbose;
413  needs_input_nodes = false;
414 };
415