Analysis Software
Documentation for sPHENIX simulation 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 // One-stop header
2 // Must include first to avoid conflict with "ClassDef" in Rtypes.h
3 #pragma GCC diagnostic push
4 #pragma GCC diagnostic ignored "-Wdeprecated-declarations"
5 #pragma GCC diagnostic ignored "-Wunused-parameter"
6 #pragma GCC diagnostic ignored "-Wshadow"
7 #include <torch/script.h>
8 #pragma GCC diagnostic pop
9 
10 #include "TpcClusterizer.h"
11 
12 #include "TrainingHitsContainer.h"
13 #include "TrainingHits.h"
14 
15 #include <trackbase/ClusHitsVerbosev1.h>
16 #include <trackbase/TpcDefs.h>
17 #include <trackbase/TrkrClusterContainerv4.h>
18 #include <trackbase/TrkrClusterv3.h>
19 #include <trackbase/TrkrClusterv4.h>
20 #include <trackbase/TrkrClusterv5.h>
21 #include <trackbase/TrkrClusterHitAssocv3.h>
22 #include <trackbase/TrkrDefs.h> // for hitkey, getLayer
23 #include <trackbase/TrkrHit.h>
24 #include <trackbase/TrkrHitSet.h>
25 #include <trackbase/TrkrHitSetContainer.h>
26 #include <trackbase/alignmentTransformationContainer.h>
27 
28 #include <trackbase/RawHit.h>
29 #include <trackbase/RawHitSet.h>
30 #include <trackbase/RawHitSetv1.h>
31 #include <trackbase/RawHitSetContainer.h>
32 
33 #include <fun4all/Fun4AllReturnCodes.h>
34 #include <fun4all/SubsysReco.h> // for SubsysReco
35 
36 #include <g4detectors/PHG4TpcCylinderGeom.h>
37 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
38 
39 #include <Acts/Definitions/Units.hpp>
40 #include <Acts/Surfaces/Surface.hpp>
41 
42 #include <phool/PHCompositeNode.h>
43 #include <phool/PHIODataNode.h> // for PHIODataNode
44 #include <phool/PHNode.h> // for PHNode
45 #include <phool/PHNodeIterator.h>
46 #include <phool/PHObject.h> // for PHObject
47 #include <phool/getClass.h>
48 #include <phool/phool.h> // for PHWHERE
49 
50 #include <TMatrixFfwd.h> // for TMatrixF
51 #include <TMatrixT.h> // for TMatrixT, ope...
52 #include <TMatrixTUtils.h> // for TMatrixTRow
53 
54 #include <TFile.h>
55 
56 #include <cmath> // for sqrt, cos, sin
57 #include <iostream>
58 #include <map> // for _Rb_tree_cons...
59 #include <string>
60 #include <utility> // for pair
61 #include <array>
62 #include <vector>
63 #include <limits>
64 // Terra incognita....
65 #include <pthread.h>
66 
67 namespace
68 {
69  template<class T> inline constexpr T square( const T& x ) { return x*x; }
70 
71  using assoc = std::pair<TrkrDefs::cluskey, TrkrDefs::hitkey>;
72 
73  struct ihit
74  {
75  unsigned short iphi = 0;
76  unsigned short it = 0;
77  unsigned short adc = 0;
78  unsigned short edge = 0;
79  };
80 
81  using vec_dVerbose = std::vector<std::vector<std::pair<int,int>>>;
82 
83  // Neural network parameters and modules
84  bool gen_hits = false;
85  bool use_nn = false;
86  const int nd = 5;
87  torch::jit::script::Module module_pos;
88 
89  struct thread_data
90  {
91  PHG4TpcCylinderGeom *layergeom = nullptr;
92  TrkrHitSet *hitset = nullptr;
93  RawHitSetv1 *rawhitset = nullptr;
94  ActsGeometry *tGeometry = nullptr;
95  unsigned int layer = 0;
96  int side = 0;
97  unsigned int sector = 0;
98  float radius = 0;
99  float drift_velocity = 0;
100  unsigned short pads_per_sector = 0;
101  float phistep = 0;
102  float pedestal = 0;
103  float seed_threshold = 0;
104  float edge_threshold = 0;
105  float min_err_squared = 0;
106  float min_clus_size = 0;
107  float min_adc_sum = 0;
108  bool do_assoc = true;
109  bool do_wedge_emulation = true;
110  bool do_singles = true;
111  bool do_split = true;
112  unsigned short phibins = 0;
113  unsigned short phioffset = 0;
114  unsigned short tbins = 0;
115  unsigned short toffset = 0;
116  unsigned short maxHalfSizeT = 0;
117  unsigned short maxHalfSizePhi = 0;
118  double m_tdriftmax = 0;
119  double sampa_tbias = 0;
120  std::vector<assoc> association_vector;
121  std::vector<TrkrCluster*> cluster_vector;
122  std::vector<TrainingHits*> v_hits;
123  int verbosity = 0;
124  bool fillClusHitsVerbose = false;
125  vec_dVerbose phivec_ClusHitsVerbose ; // only fill if fillClusHitsVerbose
126  vec_dVerbose zvec_ClusHitsVerbose ;// only fill if fillClusHitsVerbose
127  };
128 
129  pthread_mutex_t mythreadlock;
130 
131  void remove_hit(double adc, int phibin, int tbin, int edge, std::multimap<unsigned short, ihit> &all_hit_map, std::vector<std::vector<unsigned short>> &adcval)
132  {
133  typedef std::multimap<unsigned short, ihit>::iterator hit_iterator;
134  std::pair<hit_iterator, hit_iterator> iterpair = all_hit_map.equal_range(adc);
135  hit_iterator it = iterpair.first;
136  for (; it != iterpair.second; ++it) {
137  if (it->second.iphi == phibin && it->second.it == tbin) {
138  all_hit_map.erase(it);
139  break;
140  }
141  }
142  if(edge)
143  adcval[phibin][tbin] = USHRT_MAX;
144  else
145  adcval[phibin][tbin] = 0;
146  }
147 
148  void remove_hits(std::vector<ihit> &ihit_list, std::multimap<unsigned short, ihit> &all_hit_map,std::vector<std::vector<unsigned short>> &adcval)
149  {
150  for(auto iter = ihit_list.begin(); iter != ihit_list.end();++iter){
151  unsigned short adc = iter->adc;
152  unsigned short phibin = iter->iphi;
153  unsigned short tbin = iter->it;
154  unsigned short edge = iter->edge;
155  remove_hit(adc,phibin,tbin,edge,all_hit_map,adcval);
156  }
157  }
158 
159  void find_t_range(int phibin, int tbin, const thread_data& my_data, const std::vector<std::vector<unsigned short>> &adcval, int& tdown, int& tup, int &touch, int &edge){
160 
161  const int FitRangeT= (int) my_data.maxHalfSizeT;
162  const int NTBinsMax = (int) my_data.tbins;
163  tup = 0;
164  tdown = 0;
165  for(int it=0; it< FitRangeT; it++){
166  int ct = tbin + it;
167 
168  if(ct <= 0 || ct >= NTBinsMax){
169  // tup = it;
170  edge++;
171  break; // truncate edge
172  }
173 
174  if(adcval[phibin][ct] <= 0) {
175  break;
176  }
177  if(adcval[phibin][ct] == USHRT_MAX) {
178  touch++;
179  break;
180  }
181  if(my_data.do_split){
182  //check local minima and break at minimum.
183  if(ct<NTBinsMax-4){//make sure we stay clear from the edge
184  if(adcval[phibin][ct]+adcval[phibin][ct+1] <
185  adcval[phibin][ct+2]+adcval[phibin][ct+3]){//rising again
186  tup = it+1;
187  touch++;
188  break;
189  }
190  }
191  }
192  tup = it;
193  }
194  for(int it=0; it< FitRangeT; it++){
195  int ct = tbin - it;
196  if(ct <= 0 || ct >= NTBinsMax){
197  // tdown = it;
198  edge++;
199  break; // truncate edge
200  }
201  if(adcval[phibin][ct] <= 0) {
202  break;
203  }
204  if(adcval[phibin][ct] == USHRT_MAX) {
205  touch++;
206  break;
207  }
208  if(my_data.do_split){//check local minima and break at minimum.
209  if(ct>4){//make sure we stay clear from the edge
210  if(adcval[phibin][ct]+adcval[phibin][ct-1] <
211  adcval[phibin][ct-2]+adcval[phibin][ct-3]){//rising again
212  tdown = it+1;
213  touch++;
214  break;
215  }
216  }
217  }
218  tdown = it;
219  }
220  return;
221  }
222 
223  void find_phi_range(int phibin, int tbin, const thread_data& my_data, const std::vector<std::vector<unsigned short>> &adcval, int& phidown, int& phiup, int &touch, int &edge)
224  {
225 
226  int FitRangePHI = (int) my_data.maxHalfSizePhi;
227  int NPhiBinsMax = (int) my_data.phibins;
228  phidown = 0;
229  phiup = 0;
230  for(int iphi=0; iphi< FitRangePHI; iphi++){
231  int cphi = phibin + iphi;
232  if(cphi < 0 || cphi >= NPhiBinsMax){
233  // phiup = iphi;
234  edge++;
235  break; // truncate edge
236  }
237 
238  //break when below minimum
239  if(adcval[cphi][tbin] <= 0) {
240  // phiup = iphi;
241  break;
242  }
243  if(adcval[cphi][tbin] == USHRT_MAX) {
244  touch++;
245  break;
246  }
247  if(my_data.do_split){//check local minima and break at minimum.
248  if(cphi<NPhiBinsMax-4){//make sure we stay clear from the edge
249  if(adcval[cphi][tbin]+adcval[cphi+1][tbin] <
250  adcval[cphi+2][tbin]+adcval[cphi+3][tbin]){//rising again
251  phiup = iphi+1;
252  touch++;
253  break;
254  }
255  }
256  }
257  phiup = iphi;
258  }
259 
260  for(int iphi=0; iphi< FitRangePHI; iphi++){
261  int cphi = phibin - iphi;
262  if(cphi < 0 || cphi >= NPhiBinsMax){
263  // phidown = iphi;
264  edge++;
265  break; // truncate edge
266  }
267 
268  if(adcval[cphi][tbin] <= 0) {
269  //phidown = iphi;
270  break;
271  }
272  if(adcval[cphi][tbin] == USHRT_MAX) {
273  touch++;
274  break;
275  }
276  if(my_data.do_split){//check local minima and break at minimum.
277  if(cphi>4){//make sure we stay clear from the edge
278  if(adcval[cphi][tbin]+adcval[cphi-1][tbin] <
279  adcval[cphi-2][tbin]+adcval[cphi-3][tbin]){//rising again
280  phidown = iphi+1;
281  touch++;
282  break;
283  }
284  }
285  }
286  phidown = iphi;
287  }
288  return;
289  }
290 
291  void get_cluster(int phibin, int tbin, const thread_data& my_data, const std::vector<std::vector<unsigned short>> &adcval, std::vector<ihit> &ihit_list, int &touch, int &edge)
292  {
293  // search along phi at the peak in t
294 
295  int tup =0;
296  int tdown =0;
297  find_t_range(phibin, tbin, my_data, adcval, tdown, tup, touch, edge);
298  //now we have the t extent of the cluster, go find the phi edges
299 
300  for(int it=tbin - tdown ; it<= tbin + tup; it++){
301  int phiup = 0;
302  int phidown = 0;
303  find_phi_range(phibin, it, my_data, adcval, phidown, phiup, touch, edge);
304  for (int iphi = phibin - phidown; iphi <= (phibin + phiup); iphi++){
305  if(adcval[iphi][it]>0 && adcval[iphi][it]!=USHRT_MAX){
306  ihit hit;
307  hit.iphi = iphi;
308  hit.it = it;
309  hit.adc = adcval[iphi][it];
310  if(touch>0){
311  if((iphi == (phibin - phidown))||
312  (iphi == (phibin + phiup))){
313  hit.edge = 1;
314  }
315  }
316  ihit_list.push_back(hit);
317  }
318  }
319  }
320  return;
321  }
322 
323  void calc_cluster_parameter(const int iphi_center, const int it_center,
324  const std::vector<ihit> &ihit_list, thread_data& my_data, int ntouch, int nedge )
325  {
326  //
327  // get z range from layer geometry
328  /* these are used for rescaling the drift velocity */
329  //const double z_min = -105.5;
330  //const double z_max = 105.5;
331  // std::cout << "calc clus" << std::endl;
332  // loop over the hits in this cluster
333  double t_sum = 0.0;
334  //double phi_sum = 0.0;
335  double adc_sum = 0.0;
336  double t2_sum = 0.0;
337  // double phi2_sum = 0.0;
338 
339  double iphi_sum = 0.0;
340  double iphi2_sum = 0.0;
341 
342  double radius = my_data.layergeom->get_radius(); // returns center of layer
343 
344  int phibinhi = -1;
345  int phibinlo = 666666;
346  int tbinhi = -1;
347  int tbinlo = 666666;
348  int clus_size = ihit_list.size();
349  int max_adc = 0;
350  if(clus_size <= my_data.min_clus_size){
351  return;
352  }
353 
354  // training information
355  TrainingHits *training_hits = nullptr;
356  if(gen_hits)
357  {
358  training_hits = new TrainingHits;
359  assert(training_hits);
360  training_hits->radius = radius;
361  training_hits->phi = my_data.layergeom->get_phicenter(iphi_center+my_data.phioffset);
362  double center_t = my_data.layergeom->get_zcenter(it_center+my_data.toffset) + my_data.sampa_tbias;
363  training_hits->z = (my_data.m_tdriftmax - center_t) * my_data.tGeometry->get_drift_velocity();
364  if(my_data.side == 0)
365  training_hits->z = -training_hits->z;
366  training_hits->phistep = my_data.layergeom->get_phistep();
367  training_hits->zstep = my_data.layergeom->get_zstep() * my_data.tGeometry->get_drift_velocity();
368  training_hits->layer = my_data.layer;
369  training_hits->ntouch = ntouch;
370  training_hits->nedge = nedge;
371  training_hits->v_adc.fill(0);
372  }
373 
374  // std::cout << "process list" << std::endl;
375  std::vector<TrkrDefs::hitkey> hitkeyvec;
376 
377  // keep track of the hit locations in a given cluster
378  std::map<int,unsigned int> m_phi {};
379  std::map<int,unsigned int> m_z {};
380 
381  for(auto iter = ihit_list.begin(); iter != ihit_list.end();++iter){
382  double adc = iter->adc;
383 
384  if (adc <= 0) continue;
385  if(adc > max_adc)
386  max_adc = adc;
387  int iphi = iter->iphi + my_data.phioffset;
388  int it = iter->it + my_data.toffset;
389  if(iphi > phibinhi) phibinhi = iphi;
390  if(iphi < phibinlo) phibinlo = iphi;
391  if(it > tbinhi) tbinhi = it;
392  if(it < tbinlo) tbinlo = it;
393 
394  // update phi sums
395  // double phi_center = my_data.layergeom->get_phicenter(iphi);
396 
397  //phi_sum += phi_center * adc;
398  //phi2_sum += square(phi_center)*adc;
399  // std::cout << "phi_center: " << phi_center << " adc: " << adc <<std::endl;
400  iphi_sum += iphi * adc;
401  iphi2_sum += square(iphi)*adc;
402 
403  // update t sums
404  double t = my_data.layergeom->get_zcenter(it);
405  t_sum += t*adc;
406  t2_sum += square(t)*adc;
407 
408  adc_sum += adc;
409 
410  if (my_data.fillClusHitsVerbose) {
411  auto pnew = m_phi.try_emplace(iphi,adc);
412  if (!pnew.second) pnew.first->second += adc;
413 
414  pnew = m_z.try_emplace(it,adc);
415  if (!pnew.second) pnew.first->second += adc;
416  }
417 
418  // capture the hitkeys for all adc values above a certain threshold
419  TrkrDefs::hitkey hitkey = TpcDefs::genHitKey(iphi, it);
420  // if(adc>5)
421  hitkeyvec.push_back(hitkey);
422 
423  // training adc
424  if(gen_hits && training_hits)
425  {
426  int iphi_diff = iter->iphi - iphi_center;
427  int it_diff = iter->it - it_center;
428  if( std::abs(iphi_diff) <= nd && std::abs(it_diff) <= nd)
429  training_hits->v_adc[(iphi_diff+nd)*(2*nd+1)+(it_diff+nd)] = adc;
430  }
431  }
432  // std::cout << "done process list" << std::endl;
433  if (adc_sum < my_data.min_adc_sum){
434  hitkeyvec.clear();
435  return; // skip obvious noise "clusters"
436  }
437  // This is the global position
438  double clusiphi = iphi_sum / adc_sum;
439  double clusphi = my_data.layergeom->get_phi(clusiphi);
440 
441  float clusx = radius * cos(clusphi);
442  float clusy = radius * sin(clusphi);
443  double clust = t_sum / adc_sum;
444  // needed for surface identification
445  double zdriftlength = clust * my_data.tGeometry->get_drift_velocity();
446  // convert z drift length to z position in the TPC
447  double clusz = my_data.m_tdriftmax * my_data.tGeometry->get_drift_velocity() - zdriftlength;
448  if(my_data.side == 0)
449  clusz = -clusz;
450 
451  const double phi_cov = (iphi2_sum/adc_sum - square(clusiphi))* pow(my_data.layergeom->get_phistep(),2);
452  const double t_cov = t2_sum/adc_sum - square(clust);
453 
454  // Get the surface key to find the surface from the
455  TrkrDefs::hitsetkey tpcHitSetKey = TpcDefs::genHitSetKey( my_data.layer, my_data.sector, my_data.side );
456  Acts::Vector3 global(clusx, clusy, clusz);
457  TrkrDefs::subsurfkey subsurfkey = 0;
458 
459  Surface surface = my_data.tGeometry->get_tpc_surface_from_coords(
460  tpcHitSetKey,
461  global,
462  subsurfkey);
463 
464  if(!surface)
465  {
468  hitkeyvec.clear();
469  return;
470  }
471 
472  // Estimate the errors
473  const double phi_err_square = (phibinhi == phibinlo) ?
474  square(radius*my_data.layergeom->get_phistep())/12:
475  square(radius)*phi_cov/(adc_sum*0.14);
476 
477  const double t_err_square = (tbinhi == tbinlo) ?
478  square(my_data.layergeom->get_zstep())/12:
479  t_cov/(adc_sum*0.14);
480 
481  char tsize = tbinhi - tbinlo + 1;
482  char phisize = phibinhi - phibinlo + 1;
483  // 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:
484  // 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
485  // - N is the number of electrons that drift to the readout plane
486  // We have to convert (sum of adc units for all bins in the cluster) to number of ionization electrons N
487  // 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
488  // To get equivalent charge per T bin, so that summing ADC input voltage over all T bins returns total input charge, divide voltages by 2.4 for 80 ns SAMPA
489  // Equivalent charge per T 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
490 
491  // SAMPA shaping bias correction
492  clust = clust + my_data.sampa_tbias;
493 
495  global *= Acts::UnitConstants::cm;
496  //std::cout << "transform" << std::endl;
497  Acts::Vector3 local = surface->transform(my_data.tGeometry->geometry().getGeoContext()).inverse() * global;
498  local /= Acts::UnitConstants::cm;
499  //std::cout << "done transform" << std::endl;
500  // we need the cluster key and all associated hit keys (note: the cluster key includes the hitset key)
501 
502  TrkrCluster *clus_base = nullptr;
503  bool b_made_cluster { false };
504 
505 
506  //std::cout << "ver5" << std::endl;
507  // std::cout << "clus num" << my_data.cluster_vector.size() << " X " << local(0) << " Y " << clust << std::endl;
508  if(sqrt(phi_err_square) > my_data.min_err_squared){
509  auto clus = new TrkrClusterv5;
510  //auto clus = std::make_unique<TrkrClusterv3>();
511  clus_base = clus;
512  clus->setAdc(adc_sum);
513  clus->setMaxAdc(max_adc);
514  clus->setEdge(nedge);
515  clus->setPhiSize(phisize);
516  clus->setZSize(tsize);
517  clus->setSubSurfKey(subsurfkey);
518  clus->setOverlap(ntouch);
519  clus->setLocalX(local(0));
520  clus->setLocalY(clust);
521  clus->setPhiError(sqrt(phi_err_square));
522  clus->setZError(sqrt(t_err_square * pow(my_data.tGeometry->get_drift_velocity(),2)));
523  my_data.cluster_vector.push_back(clus);
524  b_made_cluster = true;
525 
526  }
527 
528  if(use_nn && clus_base && training_hits)
529  {
530  try
531  {
532  // Create a vector of inputs
533  std::vector<torch::jit::IValue> inputs;
534  inputs.emplace_back(torch::stack({
535  torch::from_blob(std::vector<float>(training_hits->v_adc.begin(), training_hits->v_adc.end()).data(), {1, 2*nd+1, 2*nd+1}, torch::kFloat32),
536  torch::full({1, 2*nd+1, 2*nd+1}, std::clamp((training_hits->layer - 7) / 16, 0, 2), torch::kFloat32),
537  torch::full({1, 2*nd+1, 2*nd+1}, training_hits->z / radius, torch::kFloat32)
538  }, 1));
539 
540  // Execute the model and turn its output into a tensor
541  at::Tensor ten_pos = module_pos.forward(inputs).toTensor();
542  float nn_phi = training_hits->phi + std::clamp(ten_pos[0][0][0].item<float>(), -(float)nd, (float)nd) * training_hits->phistep;
543  float nn_z = training_hits->z + std::clamp(ten_pos[0][1][0].item<float>(), -(float)nd, (float)nd) * training_hits->zstep;
544  float nn_x = radius * cos(nn_phi);
545  float nn_y = radius * sin(nn_phi);
546  Acts::Vector3 nn_global(nn_x, nn_y, nn_z);
547  nn_global *= Acts::UnitConstants::cm;
548  Acts::Vector3 nn_local = surface->transform(my_data.tGeometry->geometry().geoContext).inverse() * nn_global;
549  nn_local /= Acts::UnitConstants::cm;
550  float nn_t = my_data.m_tdriftmax - fabs(nn_z) / my_data.tGeometry->get_drift_velocity();
551  clus_base->setLocalX(nn_local(0));
552  clus_base->setLocalY(nn_t);
553  }
554  catch(const c10::Error &e)
555  {
556  std::cout << PHWHERE << "Error: Failed to execute NN modules" << std::endl;
557  }
558  } // use_nn
559 
560  if (my_data.fillClusHitsVerbose && b_made_cluster) {
561  // push the data back to
562  my_data.phivec_ClusHitsVerbose .push_back( std::vector<std::pair<int,int>> {} );
563  my_data.zvec_ClusHitsVerbose .push_back( std::vector<std::pair<int,int>> {} );
564 
565  auto& vphi = my_data.phivec_ClusHitsVerbose .back();
566  auto& vz = my_data.zvec_ClusHitsVerbose .back();
567 
568  for (auto& entry : m_phi ) vphi.push_back({entry.first, entry.second});
569  for (auto& entry : m_z ) vz .push_back({entry.first, entry.second});
570  }
571 
572  //std::cout << "end clus out" << std::endl;
573  // if(my_data.do_assoc && my_data.clusterhitassoc){
574  if(my_data.do_assoc)
575  {
576  // get cluster index in vector. It is used to store associations, and build relevant cluster keys when filling the containers
577  uint32_t index = my_data.cluster_vector.size()-1;
578  for (unsigned int i = 0; i < hitkeyvec.size(); i++){
579  my_data.association_vector.emplace_back(index, hitkeyvec[i]);
580  }
581  if(gen_hits && training_hits)
582  training_hits->cluskey = TrkrDefs::genClusKey(tpcHitSetKey, index);
583  }
584  hitkeyvec.clear();
585  if(gen_hits && training_hits)
586  my_data.v_hits.emplace_back(training_hits);
587  // std::cout << "done calc" << std::endl;
588  }
589 
590  void ProcessSectorData(thread_data* my_data) {
591 
592  const auto& pedestal = my_data->pedestal;
593  const auto& phibins = my_data->phibins;
594  const auto& phioffset = my_data->phioffset;
595  const auto& tbins = my_data->tbins ;
596  const auto& toffset = my_data->toffset ;
597  const auto& layer = my_data->layer ;
598  // int nhits = 0;
599  // for convenience, create a 2D vector to store adc values in and initialize to zero
600  std::vector<std::vector<unsigned short>> adcval(phibins, std::vector<unsigned short>(tbins, 0));
601  std::multimap<unsigned short, ihit> all_hit_map;
602  std::vector<ihit> hit_vect;
603 
604  int tbinmax = tbins;
605  int tbinmin = 0;
606  if(my_data->do_wedge_emulation){
607  if(layer>=7 && layer <22){
608  int etacut = (tbins / 2.) - ((50 + (layer - 7)) / 105.5) * (tbins / 2.);
609  tbinmin = etacut;
610  tbinmax -= etacut;
611  }
612  if(layer>=22 && layer <=48){
613  int etacut = (tbins / 2.) - ((65 + ((40.5 / 26) * (layer - 22))) / 105.5) * (tbins / 2.);
614  tbinmin = etacut;
615  tbinmax -= etacut;
616  }
617  }
618 
619  if( my_data->hitset!=nullptr){
620  TrkrHitSet *hitset = my_data->hitset;
621  TrkrHitSet::ConstRange hitrangei = hitset->getHits();
622 
623  for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
624  hitr != hitrangei.second;
625  ++hitr){
626 
627  if( TpcDefs::getPad(hitr->first) - phioffset < 0 ){
628  //std::cout << "WARNING phibin out of range: " << TpcDefs::getPad(hitr->first) - phioffset << " | " << phibins << std::endl;
629  continue;
630  }
631  if( TpcDefs::getTBin(hitr->first) - toffset < 0 ){
632  //std::cout << "WARNING tbin out of range: " << TpcDefs::getTBin(hitr->first) - toffset << " | " << tbins <<std::endl;
633  }
634  unsigned short phibin = TpcDefs::getPad(hitr->first) - phioffset;
635  unsigned short tbin = TpcDefs::getTBin(hitr->first) - toffset;
636  unsigned short tbinorg = TpcDefs::getTBin(hitr->first);
637  if(phibin>=phibins){
638  //std::cout << "WARNING phibin out of range: " << phibin << " | " << phibins << std::endl;
639  continue;
640  }
641  if(tbin>=tbins){
642  //std::cout << "WARNING z bin out of range: " << tbin << " | " << tbins << std::endl;
643  continue;
644  }
645  if(tbinorg>tbinmax||tbinorg<tbinmin)
646  continue;
647  float_t fadc = (hitr->second->getAdc()) - pedestal; // proper int rounding +0.5
648  unsigned short adc = 0;
649  if(fadc>0) adc = (unsigned short) fadc;
650  if(phibin >= phibins) continue;
651  if(tbin >= tbins) continue; // tbin is unsigned int, <0 cannot happen
652 
653  if(adc>0){
654  if(adc>(my_data->seed_threshold)){
655  ihit thisHit;
656 
657  thisHit.iphi = phibin;
658  thisHit.it = tbin;
659  thisHit.adc = adc;
660  thisHit.edge = 0;
661  all_hit_map.insert(std::make_pair(adc, thisHit));
662  }
663  if(adc>my_data->edge_threshold){
664  adcval[phibin][tbin] = (unsigned short) adc;
665  }
666  }
667  }
668  }else if( my_data->rawhitset!=nullptr){
669  RawHitSetv1 *hitset = my_data->rawhitset;
670  /*std::cout << "Layer: " << my_data->layer
671  << "Side: " << my_data->side
672  << "Sector: " << my_data->sector
673  << " nhits: " << hitset.size()
674  << std::endl;
675  */
676  for(int nphi= 0; nphi < phibins;nphi++){
677  // nhits += hitset->m_tpchits[nphi].size();
678  if(hitset->m_tpchits[nphi].size()==0) continue;
679 
680  int pindex = 0;
681  for(unsigned int nt = 0;nt<hitset->m_tpchits[nphi].size();nt++){
682  unsigned short val = hitset->m_tpchits[nphi][nt];
683 
684  if(val==0)
685  pindex++;
686  else{
687  if(nt==0){
688  if(val>5){
689  ihit thisHit;
690  thisHit.iphi = nphi;
691  thisHit.it = pindex;
692  thisHit.adc = val;
693  thisHit.edge = 0;
694  all_hit_map.insert(std::make_pair(val, thisHit));
695  }
696  adcval[nphi][pindex++]=val;
697  }else{
698  if((hitset->m_tpchits[nphi][nt-1]==0)&&(hitset->m_tpchits[nphi][nt+1]==0))//found zero count
699  pindex+=val;
700  else{
701  if(val>5){
702  ihit thisHit;
703  thisHit.iphi = nphi;
704  thisHit.it = pindex;
705  thisHit.adc = val;
706  thisHit.edge = 0;
707  all_hit_map.insert(std::make_pair(val, thisHit));
708  }
709  adcval[nphi][pindex++]=val;
710  }
711  }
712  }
713  }
714  }
715  }
716 
717  if(my_data->do_singles){
718  for(auto ahit:all_hit_map){
719  ihit hiHit = ahit.second;
720  int iphi = hiHit.iphi;
721  int it = hiHit.it;
722  unsigned short edge = hiHit.edge;
723  double adc = hiHit.adc;
724  if(it>0&&it<tbins){
725  if(adcval[iphi][it-1]==0&&
726  adcval[iphi][it+1]==0){
727  remove_hit(adc, iphi, it, edge, all_hit_map, adcval);
728 
729  }
730  }
731  }
732  }
733 
734  // std::cout << "done filling " << std::endl;
735  while(all_hit_map.size()>0){
736  //std::cout << "all hit map size: " << all_hit_map.size() << std::endl;
737  auto iter = all_hit_map.rbegin();
738  if(iter == all_hit_map.rend()){
739  break;
740  }
741  ihit hiHit = iter->second;
742  int iphi = hiHit.iphi;
743  int it = hiHit.it;
744  //put all hits in the all_hit_map (sorted by adc)
745  //start with highest adc hit
746  // -> cluster around it and get vector of hits
747  std::vector<ihit> ihit_list;
748  int ntouch = 0;
749  int nedge =0;
750  get_cluster(iphi, it, *my_data, adcval, ihit_list, ntouch, nedge );
751 
752  // -> calculate cluster parameters
753  // -> add hits to truth association
754  // remove hits from all_hit_map
755  // repeat untill all_hit_map empty
756  calc_cluster_parameter(iphi, it, ihit_list, *my_data, ntouch, nedge );
757  remove_hits(ihit_list,all_hit_map, adcval);
758  ihit_list.clear();
759  }
760  /* if( my_data->rawhitset!=nullptr){
761  RawHitSetv1 *hitset = my_data->rawhitset;
762  std::cout << "Layer: " << my_data->layer
763  << " Side: " << my_data->side
764  << " Sector: " << my_data->sector
765  << " nhits: " << hitset->size()
766  << " nhits coutn : " << nhits
767  << " nclus: " << my_data->cluster_vector.size()
768  << std::endl;
769  }
770  */
771  // pthread_exit(nullptr);
772  }
773  void *ProcessSector(void *threadarg) {
774 
775  auto my_data = static_cast<thread_data*>(threadarg);
776  ProcessSectorData(my_data);
777  pthread_exit(nullptr);
778  }
779 }
780 
781 TpcClusterizer::TpcClusterizer(const std::string &name)
782  : SubsysReco(name),
783  m_training(nullptr)
784 {}
785 
786 bool TpcClusterizer::is_in_sector_boundary(int phibin, int sector, PHG4TpcCylinderGeom *layergeom) const
787 {
788  bool reject_it = false;
789 
790  // sector boundaries occur every 1/12 of the full phi bin range
791  int PhiBins = layergeom->get_phibins();
792  int PhiBinsSector = PhiBins/12;
793 
794  double radius = layergeom->get_radius();
795  double PhiBinSize = 2.0* radius * M_PI / (double) PhiBins;
796 
797  // sector starts where?
798  int sector_lo = sector * PhiBinsSector;
799  int sector_hi = sector_lo + PhiBinsSector - 1;
800 
801  int sector_fiducial_bins = (int) (SectorFiducialCut / PhiBinSize);
802 
803  if(phibin < sector_lo + sector_fiducial_bins || phibin > sector_hi - sector_fiducial_bins)
804  {
805  reject_it = true;
806  /*
807  int layer = layergeom->get_layer();
808  std::cout << " local maximum is in sector fiducial boundary: layer " << layer << " radius " << radius << " sector " << sector
809  << " PhiBins " << PhiBins << " sector_fiducial_bins " << sector_fiducial_bins
810  << " PhiBinSize " << PhiBinSize << " phibin " << phibin << " sector_lo " << sector_lo << " sector_hi " << sector_hi << std::endl;
811  */
812  }
813 
814  return reject_it;
815 }
816 
817 
818 int TpcClusterizer::InitRun(PHCompositeNode *topNode)
819 {
820 
821  PHNodeIterator iter(topNode);
822 
823  // Looking for the DST node
824  PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
825  if (!dstNode)
826  {
827  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
828  return Fun4AllReturnCodes::ABORTRUN;
829  }
830 
831  // Create the Cluster node if required
832  auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
833  if (!trkrclusters)
834  {
835  PHNodeIterator dstiter(dstNode);
836  PHCompositeNode *DetNode =
837  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
838  if (!DetNode)
839  {
840  DetNode = new PHCompositeNode("TRKR");
841  dstNode->addNode(DetNode);
842  }
843 
844  trkrclusters = new TrkrClusterContainerv4;
845  PHIODataNode<PHObject> *TrkrClusterContainerNode =
846  new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
847  DetNode->addNode(TrkrClusterContainerNode);
848  }
849 
850  auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
851  if (!clusterhitassoc)
852  {
853  PHNodeIterator dstiter(dstNode);
854  PHCompositeNode *DetNode =
855  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
856  if (!DetNode)
857  {
858  DetNode = new PHCompositeNode("TRKR");
859  dstNode->addNode(DetNode);
860  }
861 
862  clusterhitassoc = new TrkrClusterHitAssocv3;
863  PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(clusterhitassoc, "TRKR_CLUSTERHITASSOC", "PHObject");
864  DetNode->addNode(newNode);
865  }
866 
867  auto training_container = findNode::getClass<TrainingHitsContainer>(dstNode, "TRAINING_HITSET");
868  if (!training_container)
869  {
870  PHNodeIterator dstiter(dstNode);
871  PHCompositeNode *DetNode =
872  dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
873  if (!DetNode)
874  {
875  DetNode = new PHCompositeNode("TRKR");
876  dstNode->addNode(DetNode);
877  }
878 
879  training_container = new TrainingHitsContainer;
880  PHIODataNode<PHObject> *TrainingHitsContainerNode =
881  new PHIODataNode<PHObject>(training_container, "TRAINING_HITSET", "PHObject");
882  DetNode->addNode(TrainingHitsContainerNode);
883  }
884 
885  gen_hits = _store_hits || _use_nn;
886  use_nn = _use_nn;
887  if(use_nn)
888  {
889  const char *offline_main = std::getenv("OFFLINE_MAIN");
890  assert(offline_main);
891  std::string net_model = std::string(offline_main) + "/share/tpc/net_model.pt";
892  try
893  {
894  // Deserialize the ScriptModule from a file using torch::jit::load()
895  module_pos = torch::jit::load(net_model);
896  std::cout << PHWHERE << "Load NN module: " << net_model << std::endl;
897  }
898  catch(const c10::Error &e)
899  {
900  std::cout << PHWHERE << "Error: Cannot load module " << net_model << std::endl;
901  exit(1);
902  }
903  }
904  else
905  {
906  std::cout << PHWHERE << "Use traditional clustering" << std::endl;
907  }
908 
909  if (record_ClusHitsVerbose) {
910  // get the node
911  mClusHitsVerbose = findNode::getClass<ClusHitsVerbosev1>(topNode, "Trkr_SvtxClusHitsVerbose");
912  if (!mClusHitsVerbose)
913  {
914  PHNodeIterator dstiter(dstNode);
915  auto DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
916  if (!DetNode)
917  {
918  DetNode = new PHCompositeNode("TRKR");
919  dstNode->addNode(DetNode);
920  }
921  mClusHitsVerbose = new ClusHitsVerbosev1();
922  auto newNode = new PHIODataNode<PHObject>(mClusHitsVerbose, "Trkr_SvtxClusHitsVerbose", "PHObject");
923  DetNode->addNode(newNode);
924  }
925  }
926 
927  return Fun4AllReturnCodes::EVENT_OK;
928 }
929 
930 int TpcClusterizer::process_event(PHCompositeNode *topNode)
931 {
932  // The TPC is the only subsystem that clusters in global coordinates. For consistency,
933  // we must use the construction transforms to get the local coordinates.
934  // Set the flag to use ideal transforms for the duration of this process_event, for thread safety
935  alignmentTransformationContainer::use_alignment = false;
936 
937  // int print_layer = 18;
938 
939  if (Verbosity() > 1000)
940  std::cout << "TpcClusterizer::Process_Event" << std::endl;
941 
942  PHNodeIterator iter(topNode);
943  PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
944  if (!dstNode)
945  {
946  std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
947  return Fun4AllReturnCodes::ABORTRUN;
948  }
949  if(!do_read_raw){
950  // get node containing the digitized hits
951  m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
952  if (!m_hits)
953  {
954  std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
955  return Fun4AllReturnCodes::ABORTRUN;
956  }
957  }else{
958  // get node containing the digitized hits
959  m_rawhits = findNode::getClass<RawHitSetContainer>(topNode, "TRKR_RAWHITSET");
960  if (!m_rawhits)
961  {
962  std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
963  return Fun4AllReturnCodes::ABORTRUN;
964  }
965  }
966 
967  // get node for clusters
968  m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
969  if (!m_clusterlist)
970  {
971  std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << std::endl;
972  return Fun4AllReturnCodes::ABORTRUN;
973  }
974 
975  // get node for cluster hit associations
976  m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
977  if (!m_clusterhitassoc)
978  {
979  std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << std::endl;
980  return Fun4AllReturnCodes::ABORTRUN;
981  }
982 
983  // get node for training hits
984  m_training = findNode::getClass<TrainingHitsContainer>(topNode, "TRAINING_HITSET");
985  if (!m_training)
986  {
987  std::cout << PHWHERE << " ERROR: Can't find TRAINING_HITSET." << std::endl;
988  return Fun4AllReturnCodes::ABORTRUN;
989  }
990 
991  PHG4TpcCylinderGeomContainer *geom_container =
992  findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
993  if (!geom_container)
994  {
995  std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
996  return Fun4AllReturnCodes::ABORTRUN;
997  }
998 
999  m_tGeometry = findNode::getClass<ActsGeometry>(topNode,
1000  "ActsGeometry");
1001  if(!m_tGeometry)
1002  {
1003  std::cout << PHWHERE
1004  << "ActsGeometry not found on node tree. Exiting"
1005  << std::endl;
1006  return Fun4AllReturnCodes::ABORTRUN;
1007  }
1008 
1009  // 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
1010  // The TPC clustering is more complicated than for the silicon, because we have to deal with overlapping clusters
1011 
1012  TrkrHitSetContainer::ConstRange hitsetrange;
1013  RawHitSetContainer::ConstRange rawhitsetrange;
1014  int num_hitsets = 0;
1015 
1016  if(!do_read_raw){
1017  hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
1018  num_hitsets = std::distance(hitsetrange.first,hitsetrange.second);
1019  }else{
1020  rawhitsetrange = m_rawhits->getHitSets(TrkrDefs::TrkrId::tpcId);
1021  num_hitsets = std::distance(rawhitsetrange.first,rawhitsetrange.second);
1022  }
1023 
1024  // create structure to store given thread and associated data
1025  struct thread_pair_t
1026  {
1027  pthread_t thread;
1028  thread_data data;
1029  };
1030 
1031  // create vector of thread pairs and reserve the right size upfront to avoid reallocation
1032  std::vector<thread_pair_t> threads;
1033  threads.reserve( num_hitsets );
1034 
1035  pthread_attr_t attr;
1036  pthread_attr_init(&attr);
1037  pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
1038 
1039  if (pthread_mutex_init(&mythreadlock, nullptr) != 0)
1040  {
1041  printf("\n mutex init failed\n");
1042  return 1;
1043  }
1044  int count = 0;
1045 
1046  if(!do_read_raw){
1047  for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
1048  hitsetitr != hitsetrange.second;
1049  ++hitsetitr)
1050  {
1051  //if(count>0)continue;
1052  TrkrHitSet *hitset = hitsetitr->second;
1053  unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
1054  int side = TpcDefs::getSide(hitsetitr->first);
1055  unsigned int sector= TpcDefs::getSectorId(hitsetitr->first);
1056  PHG4TpcCylinderGeom *layergeom = geom_container->GetLayerCellGeom(layer);
1057 
1058  // instanciate new thread pair, at the end of thread vector
1059  thread_pair_t& thread_pair = threads.emplace_back();
1060  if (mClusHitsVerbose) { thread_pair.data.fillClusHitsVerbose = true; };
1061 
1062  thread_pair.data.layergeom = layergeom;
1063  thread_pair.data.hitset = hitset;
1064  thread_pair.data.rawhitset = nullptr;
1065  thread_pair.data.layer = layer;
1066  thread_pair.data.pedestal = pedestal;
1067  thread_pair.data.seed_threshold = seed_threshold;
1068  thread_pair.data.edge_threshold = edge_threshold;
1069  thread_pair.data.sector = sector;
1070  thread_pair.data.side = side;
1071  thread_pair.data.do_assoc = do_hit_assoc;
1072  thread_pair.data.do_wedge_emulation = do_wedge_emulation;
1073  thread_pair.data.do_singles = do_singles;
1074  thread_pair.data.tGeometry = m_tGeometry;
1075  thread_pair.data.maxHalfSizeT = MaxClusterHalfSizeT;
1076  thread_pair.data.maxHalfSizePhi = MaxClusterHalfSizePhi;
1077  thread_pair.data.sampa_tbias = m_sampa_tbias;
1078  thread_pair.data.verbosity = Verbosity();
1079  thread_pair.data.do_split = do_split;
1080  thread_pair.data.min_err_squared = min_err_squared;
1081  thread_pair.data.min_clus_size = min_clus_size;
1082  thread_pair.data.min_adc_sum = min_adc_sum;
1083  unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
1084  unsigned short NPhiBinsSector = NPhiBins/12;
1085  unsigned short NTBins = (unsigned short)layergeom->get_zbins();
1086  unsigned short NTBinsSide = NTBins;
1087  unsigned short NTBinsMin = 0;
1088  unsigned short PhiOffset = NPhiBinsSector * sector;
1089  unsigned short TOffset = NTBinsMin;
1090 
1091  m_tdriftmax = AdcClockPeriod * NTBins / 2.0;
1092  thread_pair.data.m_tdriftmax = m_tdriftmax;
1093 
1094  thread_pair.data.phibins = NPhiBinsSector;
1095  thread_pair.data.phioffset = PhiOffset;
1096  thread_pair.data.tbins = NTBinsSide;
1097  thread_pair.data.toffset = TOffset ;
1098 
1099  thread_pair.data.radius = layergeom->get_radius();
1100  thread_pair.data.drift_velocity = m_tGeometry->get_drift_velocity();
1101  thread_pair.data.pads_per_sector = 0;
1102  thread_pair.data.phistep = 0;
1103  int rc;
1104  rc = pthread_create(&thread_pair.thread, &attr, ProcessSector, (void *)&thread_pair.data);
1105 
1106  if (rc) {
1107  std::cout << "Error:unable to create thread," << rc << std::endl;
1108  }
1109  if(do_sequential){
1110  int rc2 = pthread_join(thread_pair.thread, nullptr);
1111  if (rc2)
1112  { std::cout << "Error:unable to join," << rc2 << std::endl; }
1113 
1114  // get the hitsetkey from thread data
1115  const auto& data( thread_pair.data );
1116  const auto hitsetkey = TpcDefs::genHitSetKey( data.layer, data.sector, data.side );
1117 
1118  // copy clusters to map
1119  for( uint32_t index = 0; index < data.cluster_vector.size(); ++index )
1120  {
1121  // generate cluster key
1122  const auto ckey = TrkrDefs::genClusKey( hitsetkey, index );
1123 
1124  // get cluster
1125  auto cluster = data.cluster_vector[index];
1126 
1127  // insert in map
1128  m_clusterlist->addClusterSpecifyKey(ckey, cluster);
1129 
1130  if (mClusHitsVerbose) {
1131  for (auto& hit : data.phivec_ClusHitsVerbose[index]) {
1132  mClusHitsVerbose->addPhiHit (hit.first, hit.second);
1133  }
1134  for (auto& hit : data.zvec_ClusHitsVerbose[index]) {
1135  mClusHitsVerbose->addZHit (hit.first, hit.second);
1136  }
1137  mClusHitsVerbose->push_hits(ckey);
1138  }
1139  }
1140 
1141  // copy hit associations to map
1142  for( const auto& [index,hkey]:thread_pair.data.association_vector)
1143  {
1144  // generate cluster key
1145  const auto ckey = TrkrDefs::genClusKey( hitsetkey, index );
1146 
1147  // add to association table
1148  m_clusterhitassoc->addAssoc(ckey,hkey);
1149  }
1150  }
1151  count++;
1152  }
1153  }else{
1154 
1155  for (RawHitSetContainer::ConstIterator hitsetitr = rawhitsetrange.first;
1156  hitsetitr != rawhitsetrange.second;
1157  ++hitsetitr){
1158  // if(count>0)continue;
1159  // const auto hitsetid = hitsetitr->first;
1160  // std::cout << " starting thread # " << count << std::endl;
1161 
1162  RawHitSet *hitset = hitsetitr->second;
1163  unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
1164  int side = TpcDefs::getSide(hitsetitr->first);
1165  unsigned int sector= TpcDefs::getSectorId(hitsetitr->first);
1166  PHG4TpcCylinderGeom *layergeom = geom_container->GetLayerCellGeom(layer);
1167 
1168  // instanciate new thread pair, at the end of thread vector
1169  thread_pair_t& thread_pair = threads.emplace_back();
1170 
1171  thread_pair.data.layergeom = layergeom;
1172  thread_pair.data.hitset = nullptr;
1173  thread_pair.data.rawhitset = dynamic_cast<RawHitSetv1 *>(hitset);
1174  thread_pair.data.layer = layer;
1175  thread_pair.data.pedestal = pedestal;
1176  thread_pair.data.sector = sector;
1177  thread_pair.data.side = side;
1178  thread_pair.data.do_assoc = do_hit_assoc;
1179  thread_pair.data.do_wedge_emulation = do_wedge_emulation;
1180  thread_pair.data.tGeometry = m_tGeometry;
1181  thread_pair.data.maxHalfSizeT = MaxClusterHalfSizeT;
1182  thread_pair.data.maxHalfSizePhi = MaxClusterHalfSizePhi;
1183  thread_pair.data.sampa_tbias = m_sampa_tbias;
1184  thread_pair.data.verbosity = Verbosity();
1185 
1186  unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
1187  unsigned short NPhiBinsSector = NPhiBins/12;
1188  unsigned short NTBins = (unsigned short)layergeom->get_zbins();
1189  unsigned short NTBinsSide = NTBins;
1190  unsigned short NTBinsMin = 0;
1191  unsigned short PhiOffset = NPhiBinsSector * sector;
1192  unsigned short TOffset = NTBinsMin;
1193 
1194  m_tdriftmax = AdcClockPeriod * NTBins / 2.0;
1195  thread_pair.data.m_tdriftmax = m_tdriftmax;
1196 
1197  thread_pair.data.phibins = NPhiBinsSector;
1198  thread_pair.data.phioffset = PhiOffset;
1199  thread_pair.data.tbins = NTBinsSide;
1200  thread_pair.data.toffset = TOffset ;
1201 
1202  /*
1203  PHG4TpcCylinderGeom *testlayergeom = geom_container->GetLayerCellGeom(32);
1204  for( float iphi = 1408; iphi < 1408+ 128;iphi+=0.1){
1205  double clusiphi = iphi;
1206  double clusphi = testlayergeom->get_phi(clusiphi);
1207  double radius = layergeom->get_radius();
1208  float clusx = radius * cos(clusphi);
1209  float clusy = radius * sin(clusphi);
1210  float clusz = -37.524;
1211 
1212  TrkrDefs::hitsetkey tpcHitSetKey = TpcDefs::genHitSetKey( 32,11, 0 );
1213  Acts::Vector3 global(clusx, clusy, clusz);
1214  TrkrDefs::subsurfkey subsurfkey = 0;
1215 
1216  Surface surface = m_tGeometry->get_tpc_surface_from_coords(
1217  tpcHitSetKey,
1218  global,
1219  subsurfkey);
1220  std::cout << " iphi: " << iphi << " clusphi: " << clusphi << " surfkey " << subsurfkey << std::endl;
1221  // std::cout << "surfkey" << subsurfkey << std::endl;
1222  }
1223  continue;
1224  */
1225  int rc = 0;
1226  // if(layer==32)
1227  rc = pthread_create(&thread_pair.thread, &attr, ProcessSector, (void *)&thread_pair.data);
1228  // else
1229  //continue;
1230 
1231  if (rc) {
1232  std::cout << "Error:unable to create thread," << rc << std::endl;
1233  }
1234 
1235  if(do_sequential){
1236  int rc2 = pthread_join(thread_pair.thread, nullptr);
1237  if (rc2)
1238  { std::cout << "Error:unable to join," << rc2 << std::endl; }
1239 
1240  // get the hitsetkey from thread data
1241  const auto& data( thread_pair.data );
1242  const auto hitsetkey = TpcDefs::genHitSetKey( data.layer, data.sector, data.side );
1243 
1244  // copy clusters to map
1245  for( uint32_t index = 0; index < data.cluster_vector.size(); ++index )
1246  {
1247  // generate cluster key
1248  const auto ckey = TrkrDefs::genClusKey( hitsetkey, index );
1249 
1250  // get cluster
1251  auto cluster = data.cluster_vector[index];
1252 
1253  // insert in map
1254  m_clusterlist->addClusterSpecifyKey(ckey, cluster);
1255  }
1256 
1257  // copy hit associations to map
1258  for( const auto& [index,hkey]:thread_pair.data.association_vector)
1259  {
1260  // generate cluster key
1261  const auto ckey = TrkrDefs::genClusKey( hitsetkey, index );
1262 
1263  // add to association table
1264  m_clusterhitassoc->addAssoc(ckey,hkey);
1265  }
1266  }
1267  count++;
1268  }
1269  }
1270 
1271 
1272  pthread_attr_destroy(&attr);
1273  count =0;
1274  // wait for completion of all threads
1275  if(!do_sequential){
1276  for( const auto& thread_pair:threads )
1277  {
1278  int rc2 = pthread_join(thread_pair.thread, nullptr);
1279  if (rc2)
1280  { std::cout << "Error:unable to join," << rc2 << std::endl; }
1281 
1282  // get the hitsetkey from thread data
1283  const auto& data( thread_pair.data );
1284  const auto hitsetkey = TpcDefs::genHitSetKey( data.layer, data.sector, data.side );
1285 
1286  // copy clusters to map
1287  for( uint32_t index = 0; index < data.cluster_vector.size(); ++index )
1288  {
1289  // generate cluster key
1290  const auto ckey = TrkrDefs::genClusKey( hitsetkey, index );
1291 
1292  // get cluster
1293  auto cluster = data.cluster_vector[index];
1294 
1295  // insert in map
1296  //std::cout << "X: " << cluster->getLocalX() << "Y: " << cluster->getLocalY() << std::endl;
1297  m_clusterlist->addClusterSpecifyKey(ckey, cluster);
1298 
1299  if (mClusHitsVerbose) {
1300  for (auto& hit : data.phivec_ClusHitsVerbose[index]) {
1301  mClusHitsVerbose->addPhiHit (hit.first, (float)hit.second);
1302  }
1303  for (auto& hit : data.zvec_ClusHitsVerbose[index]) {
1304  mClusHitsVerbose->addZHit (hit.first, (float)hit.second);
1305  }
1306  mClusHitsVerbose->push_hits(ckey);
1307  }
1308 
1309  }
1310 
1311  // copy hit associations to map
1312  for( const auto& [index,hkey]:thread_pair.data.association_vector)
1313  {
1314  // generate cluster key
1315  const auto ckey = TrkrDefs::genClusKey( hitsetkey, index );
1316 
1317  // add to association table
1318  m_clusterhitassoc->addAssoc(ckey,hkey);
1319  }
1320 
1321  for(auto hiter = thread_pair.data.v_hits.begin(); hiter != thread_pair.data.v_hits.end(); hiter++)
1322  {
1323  if(_store_hits)
1324  m_training->v_hits.emplace_back(**hiter);
1325  delete *hiter;
1326  }
1327  }
1328  }
1329 
1330  // set the flag to use alignment transformations, needed by the rest of reconstruction
1331  alignmentTransformationContainer::use_alignment = true;
1332 
1333  if (Verbosity() > 0)
1334  std::cout << "TPC Clusterizer found " << m_clusterlist->size() << " Clusters " << std::endl;
1335  return Fun4AllReturnCodes::EVENT_OK;
1336 }
1337 
1338 int TpcClusterizer::End(PHCompositeNode */*topNode*/)
1339 {
1340  return Fun4AllReturnCodes::EVENT_OK;
1341 }