ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
TpcClusterizer.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file TpcClusterizer.cc
1 #include "TpcClusterizer.h"
2 
3 #include "TpcDefs.h"
4 
8 #include <trackbase/TrkrDefs.h> // for hitkey, getLayer
9 #include <trackbase/TrkrHitv2.h>
10 #include <trackbase/TrkrHitSet.h>
12 
14 #include <fun4all/SubsysReco.h> // for SubsysReco
15 
18 
19 #include <Acts/Utilities/Units.hpp>
21 
22 #include <phool/PHCompositeNode.h>
23 #include <phool/PHIODataNode.h> // for PHIODataNode
24 #include <phool/PHNode.h> // for PHNode
25 #include <phool/PHNodeIterator.h>
26 #include <phool/PHObject.h> // for PHObject
27 #include <phool/getClass.h>
28 #include <phool/phool.h> // for PHWHERE
29 
30 #include <TMatrixFfwd.h> // for TMatrixF
31 #include <TMatrixT.h> // for TMatrixT, ope...
32 #include <TMatrixTUtils.h> // for TMatrixTRow
33 
34 #include <TFile.h>
35 
36 #include <cmath> // for sqrt, cos, sin
37 #include <iostream>
38 #include <map> // for _Rb_tree_cons...
39 #include <string>
40 #include <utility> // for pair
41 #include <array>
42 #include <vector>
43 // Terra incognita....
44 #include <pthread.h>
45 
46 namespace
47 {
48  template<class T> inline constexpr T square( const T& x ) { return x*x; }
49 
50  typedef std::pair<unsigned short, unsigned short> iphiz;
51  typedef std::pair<unsigned short, iphiz> ihit;
52 
53  struct thread_data
54  {
55  PHG4CylinderCellGeom *layergeom = nullptr;
56  TrkrHitSet *hitset = nullptr;
57  ActsSurfaceMaps *surfmaps = nullptr;
59  unsigned int layer = 0;
60  int side = 0;
61  unsigned int sector = 0;
62  float pedestal = 0;
63  bool do_assoc = true;
64  unsigned short phibins = 0;
65  unsigned short phioffset = 0;
66  unsigned short zbins = 0;
67  unsigned short zoffset = 0;
68  unsigned short maxHalfSizeZ = 0;
69  unsigned short maxHalfSizePhi = 0;
70  double m_drift_velocity_scale = 1.0;
71  double par0_neg = 0;
72  double par0_pos = 0;
73  std::map<TrkrDefs::cluskey, TrkrCluster *> *clusterlist = nullptr;
74  std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey> *clusterhitassoc = nullptr;
75  };
76 
77  pthread_mutex_t mythreadlock;
78 
79  void remove_hit(double adc, int phibin, int zbin, std::multimap<unsigned short, ihit> &all_hit_map, std::vector<std::vector<unsigned short>> &adcval)
80  {
81  typedef std::multimap<unsigned short, ihit>::iterator hit_iterator;
82  std::pair<hit_iterator, hit_iterator> iterpair = all_hit_map.equal_range(adc);
83  hit_iterator it = iterpair.first;
84  for (; it != iterpair.second; ++it) {
85  if (it->second.second.first == phibin && it->second.second.second == zbin) {
86  all_hit_map.erase(it);
87  break;
88  }
89  }
90  adcval[phibin][zbin] = 0;
91  }
92 
93  void remove_hits(std::vector<ihit> &ihit_list, std::multimap<unsigned short, ihit> &all_hit_map,std::vector<std::vector<unsigned short>> &adcval)
94  {
95  for(auto iter = ihit_list.begin(); iter != ihit_list.end();++iter){
96  unsigned short adc = iter->first;
97  unsigned short phibin = iter->second.first;
98  unsigned short zbin = iter->second.second;
99  remove_hit(adc,phibin,zbin,all_hit_map,adcval);
100  }
101  }
102 
103  void find_z_range(int phibin, int zbin, const thread_data& my_data, const std::vector<std::vector<unsigned short>> &adcval, int& zdown, int& zup){
104 
105  const int FitRangeZ = (int) my_data.maxHalfSizeZ;
106  const int NZBinsMax = (int) my_data.zbins;
107  zup = 0;
108  zdown = 0;
109  for(int iz=0; iz< FitRangeZ; iz++){
110  int cz = zbin + iz;
111 
112  if(cz <= 0 || cz >= NZBinsMax){
113  // zup = iz;
114  break; // truncate edge
115  }
116 
117  if(adcval[phibin][cz] <= 0) {
118  break;
119  }
120  //check local minima and break at minimum.
121  if(cz<NZBinsMax-4){//make sure we stay clear from the edge
122  if(adcval[phibin][cz]+adcval[phibin][cz+1] <
123  adcval[phibin][cz+2]+adcval[phibin][cz+3]){//rising again
124  zup = iz+1;
125  break;
126  }
127  }
128  zup = iz;
129  }
130  for(int iz=0; iz< FitRangeZ; iz++){
131  int cz = zbin - iz;
132  if(cz <= 0 || cz >= NZBinsMax){
133  // zdown = iz;
134  break; // truncate edge
135  }
136  if(adcval[phibin][cz] <= 0) {
137  break;
138  }
139 
140  if(cz>4){//make sure we stay clear from the edge
141  if(adcval[phibin][cz]+adcval[phibin][cz-1] <
142  adcval[phibin][cz-2]+adcval[phibin][cz-3]){//rising again
143  zdown = iz+1;
144  break;
145  }
146  }
147  zdown = iz;
148  }
149  }
150 
151  void find_phi_range(int phibin, int zbin, const thread_data& my_data, const std::vector<std::vector<unsigned short>> &adcval, int& phidown, int& phiup)
152  {
153 
154  int FitRangePHI = (int) my_data.maxHalfSizePhi;
155  int NPhiBinsMax = (int) my_data.phibins;
156  phidown = 0;
157  phiup = 0;
158  for(int iphi=0; iphi< FitRangePHI; iphi++){
159  int cphi = phibin + iphi;
160  if(cphi < 0 || cphi >= NPhiBinsMax){
161  // phiup = iphi;
162  break; // truncate edge
163  }
164 
165  //break when below minimum
166  if(adcval[cphi][zbin] <= 0) {
167  // phiup = iphi;
168  break;
169  }
170  //check local minima and break at minimum.
171  if(cphi<NPhiBinsMax-4){//make sure we stay clear from the edge
172  if(adcval[cphi][zbin]+adcval[cphi+1][zbin] <
173  adcval[cphi+2][zbin]+adcval[cphi+3][zbin]){//rising again
174  phiup = iphi+1;
175  break;
176  }
177  }
178  phiup = iphi;
179  }
180 
181  for(int iphi=0; iphi< FitRangePHI; iphi++){
182  int cphi = phibin - iphi;
183  if(cphi < 0 || cphi >= NPhiBinsMax){
184  // phidown = iphi;
185  break; // truncate edge
186  }
187 
188  if(adcval[cphi][zbin] <= 0) {
189  //phidown = iphi;
190  break;
191  }
192 
193  if(cphi>4){//make sure we stay clear from the edge
194  if(adcval[cphi][zbin]+adcval[cphi-1][zbin] <
195  adcval[cphi-2][zbin]+adcval[cphi-3][zbin]){//rising again
196  phidown = iphi+1;
197  break;
198  }
199  }
200  phidown = iphi;
201  }
202  }
203 
204  void get_cluster(int phibin, int zbin, const thread_data& my_data, const std::vector<std::vector<unsigned short>> &adcval, std::vector<ihit> &ihit_list)
205  {
206  // search along phi at the peak in z
207 
208  int zup =0;
209  int zdown =0;
210  find_z_range(phibin, zbin, my_data, adcval, zdown, zup);
211  //now we have the z extent of the cluster, go find the phi edges
212 
213  for(int iz=zbin - zdown ; iz<= zbin + zup; iz++){
214  int phiup = 0;
215  int phidown = 0;
216  find_phi_range(phibin, iz, my_data, adcval, phidown, phiup);
217  for (int iphi = phibin - phidown; iphi <= (phibin + phiup); iphi++){
218  iphiz iCoord(std::make_pair(iphi,iz));
219  ihit thisHit(adcval[iphi][iz],iCoord);
220  ihit_list.push_back(thisHit);
221  }
222  }
223  }
224 
225  Surface get_tpc_surface_from_coords(TrkrDefs::hitsetkey hitsetkey,
227  ActsSurfaceMaps *surfMaps,
230  {
231 
232  unsigned int layer = TrkrDefs::getLayer(hitsetkey);
233  std::map<unsigned int, std::vector<Surface>>::iterator mapIter;
234  mapIter = surfMaps->tpcSurfaceMap.find(layer);
235 
236  if(mapIter == surfMaps->tpcSurfaceMap.end())
237  {
238  std::cout << PHWHERE
239  << "Error: hitsetkey not found in clusterSurfaceMap, hitsetkey = "
240  << hitsetkey << std::endl;
241  return nullptr;
242  }
243 
244  double world_phi = atan2(world[1], world[0]);
245  double world_z = world[2];
246 
247  std::vector<Surface> surf_vec = mapIter->second;
248  unsigned int surf_index = 999;
249 
250  // Predict which surface index this phi and z will correspond to
251  // assumes that the vector elements are ordered positive z, -pi to pi, then negative z, -pi to pi
252  double fraction = (world_phi + M_PI) / (2.0 * M_PI);
253  double rounded_nsurf = round( (double) (surf_vec.size()/2) * fraction - 0.5);
254  unsigned int nsurf = (unsigned int) rounded_nsurf;
255  if(world_z < 0)
256  nsurf += surf_vec.size()/2;
257 
258  double surfStepPhi = tGeometry->tpcSurfStepPhi;
259  double surfStepZ = tGeometry->tpcSurfStepZ;
260 
261  Surface this_surf = surf_vec[nsurf];
262  auto vec3d = this_surf->center(tGeometry->geoContext);
263  std::vector<double> surf_center = {vec3d(0) / 10.0, vec3d(1) / 10.0, vec3d(2) / 10.0}; // convert from mm to cm
264  double surf_z = surf_center[2];
265  double surf_phi = atan2(surf_center[1], surf_center[0]);
266 
267  if( (world_phi > surf_phi - surfStepPhi / 2.0 && world_phi < surf_phi + surfStepPhi / 2.0 ) &&
268  (world_z > surf_z - surfStepZ / 2.0 && world_z < surf_z + surfStepZ / 2.0) )
269  {
270  surf_index = nsurf;
271  subsurfkey = nsurf;
272  }
273  else
274  {
275  std::cout << PHWHERE
276  << "Error: TPC surface index not defined, skipping cluster!"
277  << std::endl;
278  std::cout << " coordinates: " << world[0] << " " << world[1] << " " << world[2]
279  << " radius " << sqrt(world[0]*world[0]+world[1]*world[1]) << std::endl;
280  std::cout << " world_phi " << world_phi << " world_z " << world_z << std::endl;
281  std::cout << " surf coords: " << surf_center[0] << " " << surf_center[1] << " " << surf_center[2] << std::endl;
282  std::cout << " surf_phi " << surf_phi << " surf_z " << surf_z << std::endl;
283  std::cout << " number of surfaces " << surf_vec.size() << std::endl;
284  return nullptr;
285  }
286 
287  return surf_vec[surf_index];
288 
289  }
290 
291  void calc_cluster_parameter(const std::vector<ihit> &ihit_list,int iclus, const thread_data& my_data )
292  {
293 
294  // get z range from layer geometry
295  /* these are used for rescaling the drift velocity */
296  const double z_min = -105.5;
297  const double z_max = 105.5;
298 
299  // loop over the hits in this cluster
300  double z_sum = 0.0;
301  double phi_sum = 0.0;
302  double adc_sum = 0.0;
303  double z2_sum = 0.0;
304  double phi2_sum = 0.0;
305 
306  double radius = my_data.layergeom->get_radius(); // returns center of layer
307 
308  int phibinhi = -1;
309  int phibinlo = 666666;
310  int zbinhi = -1;
311  int zbinlo = 666666;
312  int clus_size = ihit_list.size();
313 
314  if(clus_size == 1) return;
315 
316  std::vector<TrkrDefs::hitkey> hitkeyvec;
317  for(auto iter = ihit_list.begin(); iter != ihit_list.end();++iter){
318  double adc = iter->first;
319 
320  if (adc <= 0) continue;
321 
322  int iphi = iter->second.first + my_data.phioffset;
323  int iz = iter->second.second + my_data.zoffset;
324  if(iphi > phibinhi) phibinhi = iphi;
325  if(iphi < phibinlo) phibinlo = iphi;
326  if(iz > zbinhi) zbinhi = iz;
327  if(iz < zbinlo) zbinlo = iz;
328 
329  // update phi sums
330  double phi_center = my_data.layergeom->get_phicenter(iphi);
331  phi_sum += phi_center * adc;
332  phi2_sum += square(phi_center)*adc;
333 
334  // update z sums
335  double z = my_data.layergeom->get_zcenter(iz);
336 
337  // apply drift velocity scale
338  /* this formula ensures that z remains unchanged when located on one of the readout plane, at either z_max or z_min */
339  z =
340  z*my_data.m_drift_velocity_scale +
341  (z>0 ? z_max:z_min)*(1.-my_data.m_drift_velocity_scale);
342 
343  z_sum += z*adc;
344  z2_sum += square(z)*adc;
345 
346  adc_sum += adc;
347 
348  // capture the hitkeys for all adc values above a certain threshold
350  // if(adc>5)
351  hitkeyvec.push_back(hitkey);
352  }
353  if (adc_sum < 10) return; // skip obvious noise "clusters"
354 
355  // This is the global position
356  double clusphi = phi_sum / adc_sum;
357  double clusz = z_sum / adc_sum;
358 
359  const double phi_cov = phi2_sum/adc_sum - square(clusphi);
360  const double z_cov = z2_sum/adc_sum - square(clusz);
361 
362  // create the cluster entry directly in the node tree
363 
364  TrkrDefs::cluskey ckey = TpcDefs::genClusKey( my_data.hitset->getHitSetKey(), iclus );
365 
366  auto clus = std::make_unique<TrkrClusterv3>();
367  clus->setClusKey(ckey);
368 
369  // Estimate the errors
370  const double phi_err_square = (phibinhi == phibinlo) ?
371  square(radius*my_data.layergeom->get_phistep())/12:
372  square(radius)*phi_cov/(adc_sum*0.14);
373 
374  const double z_err_square = (zbinhi == zbinlo) ?
375  square(my_data.layergeom->get_zstep())/12:
376  z_cov/(adc_sum*0.14);
377 
378  // phi_cov = (weighted mean of dphi^2) - (weighted mean of dphi)^2, which is essentially the weighted mean of dphi^2. The error is then:
379  // e_phi = sigma_dphi/sqrt(N) = sqrt( sigma_dphi^2 / N ) -- where N is the number of samples of the distribution with standard deviation sigma_dphi
380  // - N is the number of electrons that drift to the readout plane
381  // We have to convert (sum of adc units for all bins in the cluster) to number of ionization electrons N
382  // Conversion gain is 20 mV/fC - relates total charge collected on pad to PEAK voltage out of ADC. The GEM gain is assumed to be 2000
383  // To get equivalent charge per Z bin, so that summing ADC input voltage over all Z bins returns total input charge, divide voltages by 2.4 for 80 ns SAMPA
384  // Equivalent charge per Z bin is then (ADU x 2200 mV / 1024) / 2.4 x (1/20) fC/mV x (1/1.6e-04) electrons/fC x (1/2000) = ADU x 0.14
385  clusz -= (clusz<0) ? my_data.par0_neg:my_data.par0_pos;
386 
387  // Fill in the cluster details
388  //================
389  clus->setAdc(adc_sum);
390 
391  float clusx = radius * cos(clusphi);
392  float clusy = radius * sin(clusphi);
393 
394  TMatrixF ERR(3, 3);
395  ERR[0][0] = 0.0;
396  ERR[0][1] = 0.0;
397  ERR[0][2] = 0.0;
398  ERR[1][0] = 0.0;
399  ERR[1][1] = phi_err_square; //cluster_v1 expects rad, arc, z as elementsof covariance
400  ERR[1][2] = 0.0;
401  ERR[2][0] = 0.0;
402  ERR[2][1] = 0.0;
403  ERR[2][2] = z_err_square;
404 
406  TrkrDefs::hitsetkey tpcHitSetKey = TpcDefs::genHitSetKey( my_data.layer, my_data.sector, my_data.side );
407 
408  Acts::Vector3D global(clusx, clusy, clusz);
409 
411  Surface surface = get_tpc_surface_from_coords(tpcHitSetKey,
412  global,
413  my_data.surfmaps,
414  my_data.tGeometry,
415  subsurfkey);
416 
417  if(!surface)
418  {
421  return;
422  }
423 
424  clus->setSubSurfKey(subsurfkey);
425 
426  Acts::Vector3D center = surface->center(my_data.tGeometry->geoContext)/Acts::UnitConstants::cm;
427 
429  Acts::Vector3D normal = surface->normal(my_data.tGeometry->geoContext);
430  double clusRadius = sqrt(clusx * clusx + clusy * clusy);
431  double rClusPhi = clusRadius * clusphi;
432  double surfRadius = sqrt(center(0)*center(0) + center(1)*center(1));
433  double surfPhiCenter = atan2(center[1], center[0]);
434  double surfRphiCenter = surfPhiCenter * surfRadius;
435  double surfZCenter = center[2];
436 
437  auto local = surface->globalToLocal(my_data.tGeometry->geoContext,
438  global * Acts::UnitConstants::cm,
439  normal);
440  Acts::Vector2D localPos;
441 
443  if(local.ok())
444  {
445  localPos = local.value() / Acts::UnitConstants::cm;
446  }
447  else
448  {
450  localPos(0) = rClusPhi - surfRphiCenter;
451  localPos(1) = clusz - surfZCenter;
452  }
453 
454  clus->setLocalX(localPos(0));
455  clus->setLocalY(localPos(1));
456  clus->setActsLocalError(0,0, ERR[1][1]);
457  clus->setActsLocalError(1,0, ERR[2][1]);
458  clus->setActsLocalError(0,1, ERR[1][2]);
459  clus->setActsLocalError(1,1, ERR[2][2]);
460 
461  // Add the hit associations to the TrkrClusterHitAssoc node
462  // we need the cluster key and all associated hit keys (note: the cluster key includes the hitset key)
463 
464  if(my_data.clusterlist)
465  {
466  const auto [iter, inserted] = my_data.clusterlist->insert(std::make_pair(ckey, clus.get()));
467 
468  /*
469  * if cluster was properly inserted in the map, one should release the unique_ptr,
470  * to make sure that the cluster is not deleted when going out-of-scope
471  */
472  if( inserted ) clus.release();
473  else {
474  // print error message. Duplicated cluster keys should not happen
475  std::cout
476  << PHWHERE
477  << "Error: duplicated cluster key: " << ckey << " - new cluster not inserted in map"
478  << std::endl;
479  }
480 
481  }
482 
483  if(my_data.do_assoc && my_data.clusterhitassoc){
484  for (unsigned int i = 0; i < hitkeyvec.size(); i++){
485  my_data.clusterhitassoc->insert(std::make_pair(ckey, hitkeyvec[i]));
486  }
487  }
488  }
489 
490  void *ProcessSector(void *threadarg) {
491 
492  auto my_data = static_cast<thread_data*>(threadarg);
493  const auto& pedestal = my_data->pedestal;
494  const auto& phibins = my_data->phibins;
495  const auto& phioffset = my_data->phioffset;
496  const auto& zbins = my_data->zbins ;
497  const auto& zoffset = my_data->zoffset ;
498 
499  TrkrHitSet *hitset = my_data->hitset;
500  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
501 
502  // for convenience, create a 2D vector to store adc values in and initialize to zero
503  std::vector<std::vector<unsigned short>> adcval(phibins, std::vector<unsigned short>(zbins, 0));
504  std::multimap<unsigned short, ihit> all_hit_map;
505  std::vector<ihit> hit_vect;
506 
507  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
508  hitr != hitrangei.second;
509  ++hitr){
510  unsigned short phibin = TpcDefs::getPad(hitr->first) - phioffset;
511  unsigned short zbin = TpcDefs::getTBin(hitr->first) - zoffset;
512 
513  float_t fadc = (hitr->second->getAdc()) - pedestal; // proper int rounding +0.5
514  //std::cout << " layer: " << my_data->layer << " phibin " << phibin << " zbin " << zbin << " fadc " << hitr->second->getAdc() << " pedestal " << pedestal << " fadc " << std::endl
515 
516  unsigned short adc = 0;
517  if(fadc>0)
518  adc = (unsigned short) fadc;
519 
520  // if(phibin < 0) continue; // phibin is unsigned int, <0 cannot happen
521  if(phibin >= phibins) continue;
522  // if(zbin < 0) continue;
523  if(zbin >= zbins) continue; // zbin is unsigned int, <0 cannot happen
524 
525  if(adc>0){
526  iphiz iCoord(std::make_pair(phibin,zbin));
527  ihit thisHit(adc,iCoord);
528  if(adc>5){
529  all_hit_map.insert(std::make_pair(adc, thisHit));
530  }
531  //adcval[phibin][zbin] = (unsigned short) adc;
532  adcval[phibin][zbin] = (unsigned short) adc;
533  }
534  }
535 
536  int nclus = 0;
537  while(all_hit_map.size()>0){
538 
539  auto iter = all_hit_map.rbegin();
540  if(iter == all_hit_map.rend()){
541  break;
542  }
543  ihit hiHit = iter->second;
544  int iphi = hiHit.second.first;
545  int iz = hiHit.second.second;
546 
547  //put all hits in the all_hit_map (sorted by adc)
548  //start with highest adc hit
549  // -> cluster around it and get vector of hits
550  std::vector<ihit> ihit_list;
551  get_cluster(iphi, iz, *my_data, adcval, ihit_list);
552  nclus++;
553 
554  // -> calculate cluster parameters
555  // -> add hits to truth association
556  // remove hits from all_hit_map
557  // repeat untill all_hit_map empty
558  calc_cluster_parameter(ihit_list,nclus++, *my_data );
559  remove_hits(ihit_list,all_hit_map, adcval);
560  }
561  pthread_exit(nullptr);
562  }
563 }
564 
566  : SubsysReco(name)
567 {}
568 
569 bool TpcClusterizer::is_in_sector_boundary(int phibin, int sector, PHG4CylinderCellGeom *layergeom) const
570 {
571  bool reject_it = false;
572 
573  // sector boundaries occur every 1/12 of the full phi bin range
574  int PhiBins = layergeom->get_phibins();
575  int PhiBinsSector = PhiBins/12;
576 
577  double radius = layergeom->get_radius();
578  double PhiBinSize = 2.0* radius * M_PI / (double) PhiBins;
579 
580  // sector starts where?
581  int sector_lo = sector * PhiBinsSector;
582  int sector_hi = sector_lo + PhiBinsSector - 1;
583 
584  int sector_fiducial_bins = (int) (SectorFiducialCut / PhiBinSize);
585 
586  if(phibin < sector_lo + sector_fiducial_bins || phibin > sector_hi - sector_fiducial_bins)
587  {
588  reject_it = true;
589  /*
590  int layer = layergeom->get_layer();
591  std::cout << " local maximum is in sector fiducial boundary: layer " << layer << " radius " << radius << " sector " << sector
592  << " PhiBins " << PhiBins << " sector_fiducial_bins " << sector_fiducial_bins
593  << " PhiBinSize " << PhiBinSize << " phibin " << phibin << " sector_lo " << sector_lo << " sector_hi " << sector_hi << std::endl;
594  */
595  }
596 
597  return reject_it;
598 }
599 
600 
602 {
603  PHNodeIterator iter(topNode);
604 
605  // Looking for the DST node
606  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
607  if (!dstNode)
608  {
609  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
611  }
612 
613  // Create the Cluster node if required
614  auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
615  if (!trkrclusters)
616  {
617  PHNodeIterator dstiter(dstNode);
618  PHCompositeNode *DetNode =
619  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
620  if (!DetNode)
621  {
622  DetNode = new PHCompositeNode("TRKR");
623  dstNode->addNode(DetNode);
624  }
625 
626  trkrclusters = new TrkrClusterContainerv3;
627  PHIODataNode<PHObject> *TrkrClusterContainerNode =
628  new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
629  DetNode->addNode(TrkrClusterContainerNode);
630  }
631 
632  auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
633  if (!clusterhitassoc)
634  {
635  PHNodeIterator dstiter(dstNode);
636  PHCompositeNode *DetNode =
637  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
638  if (!DetNode)
639  {
640  DetNode = new PHCompositeNode("TRKR");
641  dstNode->addNode(DetNode);
642  }
643 
644  clusterhitassoc = new TrkrClusterHitAssocv3;
645  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(clusterhitassoc, "TRKR_CLUSTERHITASSOC", "PHObject");
646  DetNode->addNode(newNode);
647  }
648 
650 }
651 
653 {
654  // int print_layer = 18;
655 
656  if (Verbosity() > 1000)
657  std::cout << "TpcClusterizer::Process_Event" << std::endl;
658 
659  PHNodeIterator iter(topNode);
660  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
661  if (!dstNode)
662  {
663  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
665  }
666 
667  // get node containing the digitized hits
668  m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
669  if (!m_hits)
670  {
671  std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
673  }
674 
675  // get node for clusters
676  m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
677  if (!m_clusterlist)
678  {
679  std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << std::endl;
681  }
682 
683  // get node for cluster hit associations
684  m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
685  if (!m_clusterhitassoc)
686  {
687  std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << std::endl;
689  }
690 
691  PHG4CylinderCellGeomContainer *geom_container =
692  findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
693  if (!geom_container)
694  {
695  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
697  }
698 
699  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
700  "ActsTrackingGeometry");
701  if(!m_tGeometry)
702  {
703  std::cout << PHWHERE
704  << "ActsTrackingGeometry not found on node tree. Exiting"
705  << std::endl;
707  }
708 
709  m_surfMaps = findNode::getClass<ActsSurfaceMaps>(topNode,
710  "ActsSurfaceMaps");
711  if(!m_surfMaps)
712  {
713  std::cout << PHWHERE
714  << "ActsSurfaceMaps not found on node tree. Exiting"
715  << std::endl;
717  }
718 
719  // 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
720  // The TPC clustering is more complicated than for the silicon, because we have to deal with overlapping clusters
721 
722  // loop over the TPC HitSet objects
724  const int num_hitsets = std::distance(hitsetrange.first,hitsetrange.second);
725 
726  // create structure to store given thread and associated data
727  struct thread_pair_t
728  {
729  pthread_t thread;
730  thread_data data;
731  };
732 
733  // create vector of thread pairs and reserve the right size upfront to avoid reallocation
734  std::vector<thread_pair_t> threads;
735  threads.reserve( num_hitsets );
736 
737  pthread_attr_t attr;
738  pthread_attr_init(&attr);
739  pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
740 
741  if (pthread_mutex_init(&mythreadlock, nullptr) != 0)
742  {
743  printf("\n mutex init failed\n");
744  return 1;
745  }
746 
747  for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
748  hitsetitr != hitsetrange.second;
749  ++hitsetitr)
750  {
751  const auto hitsetid = hitsetitr->first;
752  TrkrHitSet *hitset = hitsetitr->second;
753  unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
754  int side = TpcDefs::getSide(hitsetitr->first);
755  unsigned int sector= TpcDefs::getSectorId(hitsetitr->first);
756  PHG4CylinderCellGeom *layergeom = geom_container->GetLayerCellGeom(layer);
757 
758  // instanciate new thread pair, at the end of thread vector
759  thread_pair_t& thread_pair = threads.emplace_back();
760 
761  thread_pair.data.layergeom = layergeom;
762  thread_pair.data.hitset = hitset;
763  thread_pair.data.layer = layer;
764  thread_pair.data.pedestal = pedestal;
765  thread_pair.data.sector = sector;
766  thread_pair.data.side = side;
767  thread_pair.data.do_assoc = do_hit_assoc;
768  thread_pair.data.clusterlist = m_clusterlist->getClusterMap(hitsetid);
769  thread_pair.data.clusterhitassoc = m_clusterhitassoc->getClusterMap(hitsetid);
770  thread_pair.data.tGeometry = m_tGeometry;
771  thread_pair.data.surfmaps = m_surfMaps;
772  thread_pair.data.maxHalfSizeZ = MaxClusterHalfSizeZ;
773  thread_pair.data.maxHalfSizePhi = MaxClusterHalfSizePhi;
774  thread_pair.data.m_drift_velocity_scale = m_drift_velocity_scale;
775  thread_pair.data.par0_neg = par0_neg;
776  thread_pair.data.par0_pos = par0_pos;
777 
778  unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
779  unsigned short NPhiBinsSector = NPhiBins/12;
780  unsigned short NZBins = (unsigned short)layergeom->get_zbins();
781  unsigned short NZBinsSide = NZBins/2;
782  unsigned short NZBinsMin = 0;
783  unsigned short PhiOffset = NPhiBinsSector * sector;
784 
785  if (side == 0){
786  NZBinsMin = 0;
787  }
788  else{
789  NZBinsMin = NZBins / 2;
790  }
791 
792  unsigned short ZOffset = NZBinsMin;
793 
794  thread_pair.data.phibins = NPhiBinsSector;
795  thread_pair.data.phioffset = PhiOffset;
796  thread_pair.data.zbins = NZBinsSide;
797  thread_pair.data.zoffset = ZOffset ;
798 
799  int rc = pthread_create(&thread_pair.thread, &attr, ProcessSector, (void *)&thread_pair.data);
800  if (rc) {
801  std::cout << "Error:unable to create thread," << rc << std::endl;
802  }
803  }
804 
805  pthread_attr_destroy(&attr);
806 
807  // wait for completion of all threads
808  for( const auto& thread_pair:threads )
809  {
810  int rc2 = pthread_join(thread_pair.thread, nullptr);
811  if (rc2)
812  { std::cout << "Error:unable to join," << rc2 << std::endl; }
813  }
814 
815  if (Verbosity() > 0)
816  std::cout << "TPC Clusterizer found " << m_clusterlist->size() << " Clusters " << std::endl;
818 }
819 
821 {
823 }