Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-06 08:18:02 

0001 // One-stop header
0002 // Must include first to avoid conflict with "ClassDef" in Rtypes.h
0003 #include <torch/script.h>
0004 
0005 #include "TpcClusterizer.h"
0006 
0007 #include "LaserEventInfo.h"
0008 
0009 #include "TrainingHits.h"
0010 #include "TrainingHitsContainer.h"
0011 
0012 #include <trackbase/ClusHitsVerbosev1.h>
0013 #include <trackbase/TpcDefs.h>
0014 #include <trackbase/TrkrClusterContainerv4.h>
0015 #include <trackbase/TrkrClusterHitAssocv3.h>
0016 #include <trackbase/TrkrClusterv3.h>
0017 #include <trackbase/TrkrClusterv4.h>
0018 #include <trackbase/TrkrClusterv5.h>
0019 #include <trackbase/TrkrDefs.h>  // for hitkey, getLayer
0020 #include <trackbase/TrkrHit.h>
0021 #include <trackbase/TrkrHitSet.h>
0022 #include <trackbase/TrkrHitSetContainer.h>
0023 #include <trackbase/alignmentTransformationContainer.h>
0024 
0025 #include <trackbase/RawHit.h>
0026 #include <trackbase/RawHitSet.h>
0027 #include <trackbase/RawHitSetContainer.h>
0028 #include <trackbase/RawHitSet.h>
0029 
0030 #include <fun4all/Fun4AllReturnCodes.h>
0031 #include <fun4all/SubsysReco.h>  // for SubsysReco
0032 
0033 #include <g4detectors/PHG4TpcCylinderGeom.h>
0034 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0035 
0036 #include <Acts/Definitions/Units.hpp>
0037 #include <Acts/Surfaces/Surface.hpp>
0038 
0039 #include <phool/PHCompositeNode.h>
0040 #include <phool/PHIODataNode.h>  // for PHIODataNode
0041 #include <phool/PHNode.h>        // for PHNode
0042 #include <phool/PHNodeIterator.h>
0043 #include <phool/PHObject.h>  // for PHObject
0044 #include <phool/getClass.h>
0045 #include <phool/phool.h>  // for PHWHERE
0046 
0047 #include <TMatrixFfwd.h>    // for TMatrixF
0048 #include <TMatrixT.h>       // for TMatrixT, ope...
0049 #include <TMatrixTUtils.h>  // for TMatrixTRow
0050 
0051 #include <TFile.h>
0052 
0053 #include <array>
0054 #include <cmath>  // for sqrt, cos, sin
0055 #include <iostream>
0056 #include <limits>
0057 #include <map>  // for _Rb_tree_cons...
0058 #include <string>
0059 #include <utility>  // for pair
0060 #include <vector>
0061 // Terra incognita....
0062 #include <pthread.h>
0063 
0064 namespace
0065 {
0066   template <class T>
0067   inline constexpr T square(const T &x)
0068   {
0069     return x * x;
0070   }
0071 
0072   using assoc = std::pair<TrkrDefs::cluskey, TrkrDefs::hitkey>;
0073 
0074   struct ihit
0075   {
0076     unsigned short iphi = 0;
0077     unsigned short it = 0;
0078     unsigned short adc = 0;
0079     unsigned short edge = 0;
0080   };
0081 
0082   using vec_dVerbose = std::vector<std::vector<std::pair<int, int>>>;
0083 
0084   // Neural network parameters and modules
0085   bool gen_hits = false;
0086   bool use_nn = false;
0087   const int nd = 5;
0088   torch::jit::script::Module module_pos;
0089 
0090   struct thread_data
0091   {
0092     PHG4TpcCylinderGeom *layergeom = nullptr;
0093     TrkrHitSet *hitset = nullptr;
0094     RawHitSet *rawhitset = nullptr;
0095     ActsGeometry *tGeometry = nullptr;
0096     unsigned int layer = 0;
0097     int side = 0;
0098     unsigned int sector = 0;
0099     float radius = 0;
0100     float drift_velocity = 0;
0101     unsigned short pads_per_sector = 0;
0102     float phistep = 0;
0103     float pedestal = 0;
0104     float seed_threshold = 0;
0105     float edge_threshold = 0;
0106     float min_err_squared = 0;
0107     float min_clus_size = 0;
0108     float min_adc_sum = 0;
0109     bool do_assoc = true;
0110     bool do_wedge_emulation = true;
0111     bool do_singles = true;
0112     bool do_split = true;
0113     int FixedWindow = 0;
0114     unsigned short phibins = 0;
0115     unsigned short phioffset = 0;
0116     unsigned short tbins = 0;
0117     unsigned short toffset = 0;
0118     unsigned short maxHalfSizeT = 0;
0119     unsigned short maxHalfSizePhi = 0;
0120     double m_tdriftmax = 0;
0121     double sampa_tbias = 0;
0122     std::vector<assoc> association_vector;
0123     std::vector<TrkrCluster *> cluster_vector;
0124     std::vector<TrainingHits *> v_hits;
0125     int verbosity = 0;
0126     bool fillClusHitsVerbose = false;
0127     vec_dVerbose phivec_ClusHitsVerbose;  // only fill if fillClusHitsVerbose
0128     vec_dVerbose zvec_ClusHitsVerbose;    // only fill if fillClusHitsVerbose
0129   };
0130 
0131   pthread_mutex_t mythreadlock;
0132 
0133   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)
0134   {
0135     using hit_iterator = std::multimap<unsigned short, ihit>::iterator;
0136     std::pair<hit_iterator, hit_iterator> iterpair = all_hit_map.equal_range(adc);
0137     hit_iterator it = iterpair.first;
0138     for (; it != iterpair.second; ++it)
0139     {
0140       if (it->second.iphi == phibin && it->second.it == tbin)
0141       {
0142         all_hit_map.erase(it);
0143         break;
0144       }
0145     }
0146     if (edge)
0147     {
0148       adcval[phibin][tbin] = USHRT_MAX;
0149     }
0150     else
0151     {
0152       adcval[phibin][tbin] = 0;
0153     }
0154   }
0155 
0156   void remove_hits(std::vector<ihit> &ihit_list, std::multimap<unsigned short, ihit> &all_hit_map, std::vector<std::vector<unsigned short>> &adcval)
0157   {
0158     for (auto &iter : ihit_list)
0159     {
0160       unsigned short adc = iter.adc;
0161       unsigned short phibin = iter.iphi;
0162       unsigned short tbin = iter.it;
0163       unsigned short edge = iter.edge;
0164       remove_hit(adc, phibin, tbin, edge, all_hit_map, adcval);
0165     }
0166   }
0167 
0168   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)
0169   {
0170     const int FitRangeT = (int) my_data.maxHalfSizeT;
0171     const int NTBinsMax = (int) my_data.tbins;
0172     const int FixedWindow = (int) my_data.FixedWindow;
0173     tup = 0;
0174     tdown = 0;
0175     if (FixedWindow != 0)
0176     {
0177       tup = FixedWindow;
0178       tdown = FixedWindow;
0179       if (tbin + tup >= NTBinsMax)
0180       {
0181         tup = NTBinsMax - tbin - 1;
0182         edge++;
0183       }
0184       if ((tbin - tdown) <= 0)
0185       {
0186         tdown = tbin;
0187         edge++;
0188       }
0189       return;
0190     }
0191     for (int it = 0; it < FitRangeT; it++)
0192     {
0193       int ct = tbin + it;
0194 
0195       if (ct <= 0 || ct >= NTBinsMax)
0196       {
0197         // tup = it;
0198         edge++;
0199         break;  // truncate edge
0200       }
0201 
0202       if (adcval[phibin][ct] <= 0)
0203       {
0204         break;
0205       }
0206       if (adcval[phibin][ct] == USHRT_MAX)
0207       {
0208         touch++;
0209         break;
0210       }
0211       if (my_data.do_split)
0212       {
0213         // check local minima and break at minimum.
0214         if (ct < NTBinsMax - 4)
0215         {  // make sure we stay clear from the edge
0216           if (adcval[phibin][ct] + adcval[phibin][ct + 1] <
0217               adcval[phibin][ct + 2] + adcval[phibin][ct + 3])
0218           {  // rising again
0219             tup = it + 1;
0220             touch++;
0221             break;
0222           }
0223         }
0224       }
0225       tup = it;
0226     }
0227     for (int it = 0; it < FitRangeT; it++)
0228     {
0229       int ct = tbin - it;
0230       if (ct <= 0 || ct >= NTBinsMax)
0231       {
0232         //      tdown = it;
0233         edge++;
0234         break;  // truncate edge
0235       }
0236       if (adcval[phibin][ct] <= 0)
0237       {
0238         break;
0239       }
0240       if (adcval[phibin][ct] == USHRT_MAX)
0241       {
0242         touch++;
0243         break;
0244       }
0245       if (my_data.do_split)
0246       {  // check local minima and break at minimum.
0247         if (ct > 4)
0248         {  // make sure we stay clear from the edge
0249           if (adcval[phibin][ct] + adcval[phibin][ct - 1] <
0250               adcval[phibin][ct - 2] + adcval[phibin][ct - 3])
0251           {  // rising again
0252             tdown = it + 1;
0253             touch++;
0254             break;
0255           }
0256         }
0257       }
0258       tdown = it;
0259     }
0260     return;
0261   }
0262 
0263   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)
0264   {
0265     int FitRangePHI = (int) my_data.maxHalfSizePhi;
0266     int NPhiBinsMax = (int) my_data.phibins;
0267     const int FixedWindow = (int) my_data.FixedWindow;
0268     phidown = 0;
0269     phiup = 0;
0270     if (FixedWindow != 0)
0271     {
0272       phiup = FixedWindow;
0273       phidown = FixedWindow;
0274       if (phibin + phiup >= NPhiBinsMax)
0275       {
0276         phiup = NPhiBinsMax - phibin - 1;
0277         edge++;
0278       }
0279       if (phibin - phidown <= 0)
0280       {
0281         phidown = phibin;
0282         edge++;
0283       }
0284       return;
0285     }
0286     for (int iphi = 0; iphi < FitRangePHI; iphi++)
0287     {
0288       int cphi = phibin + iphi;
0289       if (cphi < 0 || cphi >= NPhiBinsMax)
0290       {
0291         // phiup = iphi;
0292         edge++;
0293         break;  // truncate edge
0294       }
0295 
0296       // break when below minimum
0297       if (adcval[cphi][tbin] <= 0)
0298       {
0299         // phiup = iphi;
0300         break;
0301       }
0302       if (adcval[cphi][tbin] == USHRT_MAX)
0303       {
0304         touch++;
0305         break;
0306       }
0307       if (my_data.do_split)
0308       {  // check local minima and break at minimum.
0309         if (cphi < NPhiBinsMax - 4)
0310         {  // make sure we stay clear from the edge
0311           if (adcval[cphi][tbin] + adcval[cphi + 1][tbin] <
0312               adcval[cphi + 2][tbin] + adcval[cphi + 3][tbin])
0313           {  // rising again
0314             phiup = iphi + 1;
0315             touch++;
0316             break;
0317           }
0318         }
0319       }
0320       phiup = iphi;
0321     }
0322 
0323     for (int iphi = 0; iphi < FitRangePHI; iphi++)
0324     {
0325       int cphi = phibin - iphi;
0326       if (cphi < 0 || cphi >= NPhiBinsMax)
0327       {
0328         // phidown = iphi;
0329         edge++;
0330         break;  // truncate edge
0331       }
0332 
0333       if (adcval[cphi][tbin] <= 0)
0334       {
0335         // phidown = iphi;
0336         break;
0337       }
0338       if (adcval[cphi][tbin] == USHRT_MAX)
0339       {
0340         touch++;
0341         break;
0342       }
0343       if (my_data.do_split)
0344       {  // check local minima and break at minimum.
0345         if (cphi > 4)
0346         {  // make sure we stay clear from the edge
0347           if (adcval[cphi][tbin] + adcval[cphi - 1][tbin] <
0348               adcval[cphi - 2][tbin] + adcval[cphi - 3][tbin])
0349           {  // rising again
0350             phidown = iphi + 1;
0351             touch++;
0352             break;
0353           }
0354         }
0355       }
0356       phidown = iphi;
0357     }
0358     return;
0359   }
0360 
0361   int is_hit_isolated(int iphi, int it, int NPhiBinsMax, int NTBinsMax, const std::vector<std::vector<unsigned short>> &adcval)
0362   {
0363     // check isolated hits
0364     //  const int NPhiBinsMax = (int) my_data.phibins;
0365     // const int NTBinsMax = (int) my_data.tbins;
0366 
0367     int isosum = 0;
0368     int isophimin = iphi - 1;
0369     if (isophimin < 0)
0370     {
0371       isophimin = 0;
0372     }
0373     int isophimax = iphi + 1;
0374     if (!(isophimax < NPhiBinsMax))
0375     {
0376       isophimax = NPhiBinsMax - 1;
0377     }
0378     int isotmin = it - 1;
0379     if (isotmin < 0)
0380     {
0381       isotmin = 0;
0382     }
0383     int isotmax = it + 1;
0384     if (!(isotmax < NTBinsMax))
0385     {
0386       isotmax = NTBinsMax - 1;
0387     }
0388     for (int isophi = isophimin; isophi <= isophimax; isophi++)
0389     {
0390       for (int isot = isotmin; isot <= isotmax; isot++)
0391       {
0392         if (isophi == iphi && isot == it)
0393         {
0394           continue;
0395         }
0396         /*
0397         std::cout <<" isominphi: " << isophimin
0398                   << "phi: " << isophi
0399                   <<" isomaxphi: " << isophimax
0400                   <<" maxphi: " << NPhiBinsMax
0401                   <<" isomint: " << isotmin
0402                   << " t: " << isot
0403                   <<" isomaxt: " << isotmax
0404                   << " maxt: " << NTBinsMax
0405                   << std::endl;
0406         */
0407         isosum += adcval[isophi][isot];
0408         /*
0409         std::cout << "adc " << adcval[isophi][isot]
0410                   << " isosum " << isosum
0411                   << " done"
0412                   << std::endl;
0413         */
0414       }
0415     }
0416     int isiso = 0;
0417     if (isosum == 0)
0418     {
0419       isiso = 1;
0420     }
0421     return isiso;
0422   }
0423 
0424   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)
0425   {
0426     // search along phi at the peak in t
0427     //    const int NPhiBinsMax = (int) my_data.phibins;
0428     // const int NTBinsMax = (int) my_data.tbins;
0429     int tup = 0;
0430     int tdown = 0;
0431     find_t_range(phibin, tbin, my_data, adcval, tdown, tup, touch, edge);
0432     // now we have the t extent of the cluster, go find the phi edges
0433     for (int it = tbin - tdown; it <= (tbin + tup); it++)
0434     {
0435       int phiup = 0;
0436       int phidown = 0;
0437       find_phi_range(phibin, it, my_data, adcval, phidown, phiup, touch, edge);
0438       for (int iphi = (phibin - phidown); iphi <= (phibin + phiup); iphi++)
0439       {
0440         if (adcval[iphi][it] > 0 && adcval[iphi][it] != USHRT_MAX)
0441         {
0442           if (my_data.do_singles)
0443           {
0444             if (is_hit_isolated(iphi, it, (int) my_data.phibins, (int) my_data.tbins, adcval))
0445             {
0446               continue;
0447             }
0448           }
0449           ihit hit;
0450           hit.iphi = iphi;
0451           hit.it = it;
0452 
0453           hit.adc = adcval[iphi][it];
0454           if (touch > 0)
0455           {
0456             if ((iphi == (phibin - phidown)) ||
0457                 (iphi == (phibin + phiup)))
0458             {
0459               hit.edge = 1;
0460             }
0461           }
0462           ihit_list.push_back(hit);
0463         }
0464       }
0465     }
0466     return;
0467   }
0468 
0469   void calc_cluster_parameter(const int iphi_center, const int it_center,
0470                               const std::vector<ihit> &ihit_list, thread_data &my_data, int ntouch, int nedge)
0471   {
0472     //
0473     // get z range from layer geometry
0474     /* these are used for rescaling the drift velocity */
0475     // const double z_min = -105.5;
0476     // const double z_max = 105.5;
0477     //  std::cout << "calc clus" << std::endl;
0478     //  loop over the hits in this cluster
0479     double t_sum = 0.0;
0480     // double phi_sum = 0.0;
0481     double adc_sum = 0.0;
0482     double t2_sum = 0.0;
0483     // double phi2_sum = 0.0;
0484 
0485     double iphi_sum = 0.0;
0486     double iphi2_sum = 0.0;
0487 
0488     double radius = my_data.layergeom->get_radius();  // returns center of layer
0489 
0490     int phibinhi = -1;
0491     int phibinlo = 666666;
0492     int tbinhi = -1;
0493     int tbinlo = 666666;
0494     int clus_size = ihit_list.size();
0495     int max_adc = 0;
0496     if (clus_size <= my_data.min_clus_size)
0497     {
0498       return;
0499     }
0500 
0501     // training information
0502     TrainingHits *training_hits = nullptr;
0503     if (gen_hits)
0504     {
0505       training_hits = new TrainingHits;
0506       assert(training_hits);
0507       training_hits->radius = radius;
0508       training_hits->phi = my_data.layergeom->get_phicenter(iphi_center + my_data.phioffset, my_data.side);
0509       double center_t = my_data.layergeom->get_zcenter(it_center + my_data.toffset) + my_data.sampa_tbias;
0510       training_hits->z = (my_data.m_tdriftmax - center_t) * my_data.tGeometry->get_drift_velocity();
0511       if (my_data.side == 0)
0512       {
0513         training_hits->z = -training_hits->z;
0514       }
0515       training_hits->phistep = my_data.layergeom->get_phistep();
0516       training_hits->zstep = my_data.layergeom->get_zstep() * my_data.tGeometry->get_drift_velocity();
0517       training_hits->layer = my_data.layer;
0518       training_hits->ntouch = ntouch;
0519       training_hits->nedge = nedge;
0520       training_hits->v_adc.fill(0);
0521     }
0522 
0523     //      std::cout << "process list" << std::endl;
0524     std::vector<TrkrDefs::hitkey> hitkeyvec;
0525 
0526     // keep track of the hit locations in a given cluster
0527     std::map<int, unsigned int> m_phi{};
0528     std::map<int, unsigned int> m_z{};
0529 
0530     for (auto iter : ihit_list)
0531     {
0532       int iphi = iter.iphi + my_data.phioffset;
0533       int it = iter.it + my_data.toffset;
0534       double adc = iter.adc;
0535 
0536       if (adc <= 0)
0537       {
0538         continue;
0539       }
0540 
0541       if (adc > max_adc)
0542       {
0543         max_adc = adc;
0544       }
0545 
0546       if (iphi > phibinhi)
0547       {
0548         phibinhi = iphi;
0549       }
0550 
0551       if (iphi < phibinlo)
0552       {
0553         phibinlo = iphi;
0554       }
0555 
0556       if (it > tbinhi)
0557       {
0558         tbinhi = it;
0559       }
0560 
0561       if (it < tbinlo)
0562       {
0563         tbinlo = it;
0564       }
0565 
0566       // if(it==it_center){ yg_sum += adc; }
0567       // update phi sums
0568       //    double phi_center = my_data.layergeom->get_phicenter(iphi);
0569 
0570       // phi_sum += phi_center * adc;
0571       // phi2_sum += square(phi_center)*adc;
0572       //    std::cout << "phi_center: " << phi_center << " adc: " << adc <<std::endl;
0573       iphi_sum += iphi * adc;
0574       iphi2_sum += square(iphi) * adc;
0575 
0576       // update t sums
0577       double t = my_data.layergeom->get_zcenter(it);
0578       t_sum += t * adc;
0579       t2_sum += square(t) * adc;
0580 
0581       adc_sum += adc;
0582 
0583       if (my_data.fillClusHitsVerbose)
0584       {
0585         auto pnew = m_phi.try_emplace(iphi, adc);
0586         if (!pnew.second)
0587         {
0588           pnew.first->second += adc;
0589         }
0590 
0591         pnew = m_z.try_emplace(it, adc);
0592         if (!pnew.second)
0593         {
0594           pnew.first->second += adc;
0595         }
0596       }
0597 
0598       // capture the hitkeys for all adc values above a certain threshold
0599       TrkrDefs::hitkey hitkey = TpcDefs::genHitKey(iphi, it);
0600       // if(adc>5)
0601       hitkeyvec.push_back(hitkey);
0602 
0603       // training adc
0604       if (gen_hits && training_hits)
0605       {
0606         int iphi_diff = iter.iphi - iphi_center;
0607         int it_diff = iter.it - it_center;
0608         if (std::abs(iphi_diff) <= nd && std::abs(it_diff) <= nd)
0609         {
0610           training_hits->v_adc[(iphi_diff + nd) * (2 * nd + 1) + (it_diff + nd)] = adc;
0611         }
0612       }
0613     }
0614     //      std::cout << "done process list" << std::endl;
0615     if (adc_sum < my_data.min_adc_sum)
0616     {
0617       hitkeyvec.clear();
0618       return;  // skip obvious noise "clusters"
0619     }
0620 
0621     // This is the global position
0622     double clusiphi = iphi_sum / adc_sum;
0623     double clusphi = my_data.layergeom->get_phi(clusiphi, my_data.side);
0624 
0625     float clusx = radius * cos(clusphi);
0626     float clusy = radius * sin(clusphi);
0627     double clust = t_sum / adc_sum;
0628     // needed for surface identification
0629     double zdriftlength = clust * my_data.tGeometry->get_drift_velocity();
0630     // convert z drift length to z position in the TPC
0631     double clusz = my_data.m_tdriftmax * my_data.tGeometry->get_drift_velocity() - zdriftlength;
0632     if (my_data.side == 0)
0633     {
0634       clusz = -clusz;
0635     }
0636     // std::cout << " side " << my_data.side << " clusz " << clusz << " clust " << clust << " driftmax " << my_data.m_tdriftmax << std::endl;
0637     const double phi_cov = (iphi2_sum / adc_sum - square(clusiphi)) * pow(my_data.layergeom->get_phistep(), 2);
0638     const double t_cov = t2_sum / adc_sum - square(clust);
0639 
0640     // Get the surface key to find the surface from the
0641     TrkrDefs::hitsetkey tpcHitSetKey = TpcDefs::genHitSetKey(my_data.layer, my_data.sector, my_data.side);
0642     Acts::Vector3 global(clusx, clusy, clusz);
0643     TrkrDefs::subsurfkey subsurfkey = 0;
0644 
0645     Surface surface = my_data.tGeometry->get_tpc_surface_from_coords(
0646         tpcHitSetKey,
0647         global,
0648         subsurfkey);
0649 
0650     if (!surface)
0651     {
0652       /// If the surface can't be found, we can't track with it. So
0653       /// just return and don't add the cluster to the container
0654       hitkeyvec.clear();
0655       return;
0656     }
0657 
0658     // Estimate the errors
0659     // Blow up error on single pixel clusters by a factor 3 to compensate for threshold effects
0660     const double phi_err_square = (phibinhi == phibinlo) ? 9*(square(radius * my_data.layergeom->get_phistep()) / 12) : square(radius) * phi_cov / (adc_sum * 0.14);
0661   
0662   const double t_err_square = (tbinhi == tbinlo) ? 9*(square(my_data.layergeom->get_zstep()) / 12) : t_cov / (adc_sum * 0.14);
0663   
0664     char tsize = tbinhi - tbinlo + 1;
0665     char phisize = phibinhi - phibinlo + 1;
0666     // std::cout << "phisize: "  << (int) phisize << " phibinhi " << phibinhi << " phibinlo " << phibinlo << std::endl;
0667     // 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:
0668     // 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
0669     //    - N is the number of electrons that drift to the readout plane
0670     // We have to convert (sum of adc units for all bins in the cluster) to number of ionization electrons N
0671     // 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
0672     // 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
0673     // 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
0674 
0675     // SAMPA shaping bias correction
0676     clust = clust + my_data.sampa_tbias;
0677 
0678     /// convert to Acts units
0679     global *= Acts::UnitConstants::cm;
0680     // std::cout << "transform" << std::endl;
0681     Acts::Vector3 local = surface->transform(my_data.tGeometry->geometry().getGeoContext()).inverse() * global;
0682     local /= Acts::UnitConstants::cm;
0683     // std::cout << "done transform" << std::endl;
0684     //  we need the cluster key and all associated hit keys (note: the cluster key includes the hitset key)
0685 
0686     TrkrCluster *clus_base = nullptr;
0687     bool b_made_cluster{false};
0688 
0689     // std::cout << "ver5" << std::endl;
0690     //  std::cout << "clus num" << my_data.cluster_vector.size() << " X " << local(0) << " Y " << clust << std::endl;
0691     if (sqrt(phi_err_square) > my_data.min_err_squared)
0692     {
0693       auto clus = new TrkrClusterv5;
0694       // auto clus = std::make_unique<TrkrClusterv3>();
0695       clus_base = clus;
0696       clus->setAdc(adc_sum);
0697       clus->setMaxAdc(max_adc);
0698       clus->setEdge(nedge);
0699       clus->setPhiSize(phisize);
0700       clus->setZSize(tsize);
0701       clus->setSubSurfKey(subsurfkey);
0702       clus->setOverlap(ntouch);
0703       clus->setLocalX(local(0));
0704       clus->setLocalY(clust);
0705       clus->setPhiError(sqrt(phi_err_square));
0706       clus->setZError(sqrt(t_err_square * pow(my_data.tGeometry->get_drift_velocity(), 2)));
0707       my_data.cluster_vector.push_back(clus);
0708       b_made_cluster = true;
0709     }
0710 
0711     if (use_nn && clus_base && training_hits)
0712     {
0713       try
0714       {
0715         // Create a vector of inputs
0716         std::vector<torch::jit::IValue> inputs;
0717         inputs.emplace_back(torch::stack({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),
0718                                           torch::full({1, 2 * nd + 1, 2 * nd + 1}, std::clamp((training_hits->layer - 7) / 16, 0, 2), torch::kFloat32),
0719                                           torch::full({1, 2 * nd + 1, 2 * nd + 1}, training_hits->z / radius, torch::kFloat32)},
0720                                          1));
0721 
0722         // Execute the model and turn its output into a tensor
0723         at::Tensor ten_pos = module_pos.forward(inputs).toTensor();
0724         float nn_phi = training_hits->phi + std::clamp(ten_pos[0][0][0].item<float>(), -(float) nd, (float) nd) * training_hits->phistep;
0725         float nn_z = training_hits->z + std::clamp(ten_pos[0][1][0].item<float>(), -(float) nd, (float) nd) * training_hits->zstep;
0726         float nn_x = radius * std::cos(nn_phi);
0727         float nn_y = radius * std::sin(nn_phi);
0728         Acts::Vector3 nn_global(nn_x, nn_y, nn_z);
0729         nn_global *= Acts::UnitConstants::cm;
0730         Acts::Vector3 nn_local = surface->transform(my_data.tGeometry->geometry().geoContext).inverse() * nn_global;
0731         nn_local /= Acts::UnitConstants::cm;
0732         float nn_t = my_data.m_tdriftmax - std::fabs(nn_z) / my_data.tGeometry->get_drift_velocity();
0733         clus_base->setLocalX(nn_local(0));
0734         clus_base->setLocalY(nn_t);
0735       }
0736       catch (const c10::Error &e)
0737       {
0738         std::cout << PHWHERE << "Error: Failed to execute NN modules" << std::endl;
0739       }
0740     }  // use_nn
0741 
0742     if (my_data.fillClusHitsVerbose && b_made_cluster)
0743     {
0744       // push the data back to
0745       my_data.phivec_ClusHitsVerbose.push_back(std::vector<std::pair<int, int>>{});
0746       my_data.zvec_ClusHitsVerbose.push_back(std::vector<std::pair<int, int>>{});
0747 
0748       auto &vphi = my_data.phivec_ClusHitsVerbose.back();
0749       auto &vz = my_data.zvec_ClusHitsVerbose.back();
0750 
0751       for (auto &entry : m_phi)
0752       {
0753         vphi.push_back({entry.first, entry.second});
0754       }
0755       for (auto &entry : m_z)
0756       {
0757         vz.push_back({entry.first, entry.second});
0758       }
0759     }
0760 
0761     // std::cout << "end clus out" << std::endl;
0762     //       if(my_data.do_assoc && my_data.clusterhitassoc){
0763     if (my_data.do_assoc)
0764     {
0765       // get cluster index in vector. It is used to store associations, and build relevant cluster keys when filling the containers
0766       uint32_t index = my_data.cluster_vector.size() - 1;
0767       for (unsigned int &i : hitkeyvec)
0768       {
0769         my_data.association_vector.emplace_back(index, i);
0770       }
0771       if (gen_hits && training_hits)
0772       {
0773         training_hits->cluskey = TrkrDefs::genClusKey(tpcHitSetKey, index);
0774       }
0775     }
0776     hitkeyvec.clear();
0777     if (gen_hits && training_hits)
0778     {
0779       my_data.v_hits.emplace_back(training_hits);
0780     }
0781     //      std::cout << "done calc" << std::endl;
0782   }
0783 
0784   void ProcessSectorData(thread_data *my_data)
0785   {
0786     const auto &pedestal = my_data->pedestal;
0787     const auto &phibins = my_data->phibins;
0788     const auto &phioffset = my_data->phioffset;
0789     const auto &tbins = my_data->tbins;
0790     const auto &toffset = my_data->toffset;
0791     const auto &layer = my_data->layer;
0792     //    int nhits = 0;
0793     // for convenience, create a 2D vector to store adc values in and initialize to zero
0794     std::vector<std::vector<unsigned short>> adcval(phibins, std::vector<unsigned short>(tbins, 0));
0795     std::multimap<unsigned short, ihit> all_hit_map;
0796     std::vector<ihit> hit_vect;
0797 
0798     int tbinmax = tbins;
0799     int tbinmin = 0;
0800     if (my_data->do_wedge_emulation)
0801     {
0802       if (layer >= 7 && layer < 22)
0803       {
0804         int etacut = (tbins / 2.) - ((50 + (layer - 7)) / 105.5) * (tbins / 2.);
0805         tbinmin = etacut;
0806         tbinmax -= etacut;
0807       }
0808       if (layer >= 22 && layer <= 48)
0809       {
0810         int etacut = (tbins / 2.) - ((65 + ((40.5 / 26) * (layer - 22))) / 105.5) * (tbins / 2.);
0811         tbinmin = etacut;
0812         tbinmax -= etacut;
0813       }
0814     }
0815 
0816     if (my_data->hitset != nullptr)
0817     {
0818       TrkrHitSet *hitset = my_data->hitset;
0819       TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0820 
0821       for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0822            hitr != hitrangei.second;
0823            ++hitr)
0824       {
0825         if (TpcDefs::getPad(hitr->first) - phioffset < 0)
0826         {
0827           // std::cout << "WARNING phibin out of range: " << TpcDefs::getPad(hitr->first) - phioffset << " | " << phibins << std::endl;
0828           continue;
0829         }
0830         if (TpcDefs::getTBin(hitr->first) - toffset < 0)
0831         {
0832           // std::cout << "WARNING tbin out of range: " << TpcDefs::getTBin(hitr->first) - toffset  << " | " << tbins <<std::endl;
0833         }
0834         unsigned short phibin = TpcDefs::getPad(hitr->first) - phioffset;
0835         unsigned short tbin = TpcDefs::getTBin(hitr->first) - toffset;
0836         unsigned short tbinorg = TpcDefs::getTBin(hitr->first);
0837         if (phibin >= phibins)
0838         {
0839           // std::cout << "WARNING phibin out of range: " << phibin << " | " << phibins << std::endl;
0840           continue;
0841         }
0842         if (tbin >= tbins)
0843         {
0844           // std::cout << "WARNING z bin out of range: " << tbin << " | " << tbins << std::endl;
0845           continue;
0846         }
0847         if (tbinorg > tbinmax || tbinorg < tbinmin)
0848         {
0849           continue;
0850         }
0851         float_t fadc = (hitr->second->getAdc()) - pedestal;  // proper int rounding +0.5
0852         unsigned short adc = 0;
0853         if (fadc > 0)
0854         {
0855           adc = (unsigned short) fadc;
0856         }
0857         if (phibin >= phibins)
0858         {
0859           continue;
0860         }
0861         if (tbin >= tbins)
0862         {
0863           continue;  // tbin is unsigned int, <0 cannot happen
0864         }
0865 
0866         if (adc > 0)
0867         {
0868           if (adc > (my_data->seed_threshold))
0869           {
0870             ihit thisHit;
0871 
0872             thisHit.iphi = phibin;
0873             thisHit.it = tbin;
0874             thisHit.adc = adc;
0875             thisHit.edge = 0;
0876             all_hit_map.insert(std::make_pair(adc, thisHit));
0877           }
0878           if (adc > my_data->edge_threshold)
0879           {
0880             adcval[phibin][tbin] = (unsigned short) adc;
0881           }
0882         }
0883       }
0884     }
0885     else if (my_data->rawhitset != nullptr)
0886     {
0887       RawHitSet *hitset = my_data->rawhitset;
0888       /*std::cout << "Layer: " << my_data->layer
0889                 << "Side: " << my_data->side
0890                 << "Sector: " << my_data->sector
0891                 << " nhits:  " << hitset.size()
0892                 << std::endl;
0893       */
0894       for (int nphi = 0; nphi < phibins; nphi++)
0895       {
0896         if (hitset->size(nphi) == 0)
0897         {
0898           continue;
0899         }
0900 
0901         int pindex = 0;
0902         for (unsigned int nt = 0; nt < hitset->size(nphi); nt++)
0903         {
0904           unsigned short val = (*(hitset->getHits(nphi)))[nt];
0905 
0906           if (val == 0)
0907           {
0908             pindex++;
0909           }
0910           else
0911           {
0912             if (nt == 0)
0913             {
0914               if (val > 5)
0915               {
0916                 ihit thisHit;
0917                 thisHit.iphi = nphi;
0918                 thisHit.it = pindex;
0919                 thisHit.adc = val;
0920                 thisHit.edge = 0;
0921                 all_hit_map.insert(std::make_pair(val, thisHit));
0922               }
0923               adcval[nphi][pindex++] = val;
0924             }
0925             else
0926             {
0927               if (((*(hitset->getHits(nphi)))[nt - 1] == 0) && ((*(hitset->getHits(nphi)))[nt + 1] == 0))
0928               {  // found zero count
0929                 pindex += val;
0930               }
0931               else
0932               {
0933                 if (val > 5)
0934                 {
0935                   ihit thisHit;
0936                   thisHit.iphi = nphi;
0937                   thisHit.it = pindex;
0938                   thisHit.adc = val;
0939                   thisHit.edge = 0;
0940                   all_hit_map.insert(std::make_pair(val, thisHit));
0941                 }
0942                 adcval[nphi][pindex++] = val;
0943               }
0944             }
0945           }
0946         }
0947       }
0948     }
0949     /*
0950     if (my_data->do_singles)
0951     {
0952       for (auto ahit : all_hit_map)
0953       {
0954         ihit hiHit = ahit.second;
0955         int iphi = hiHit.iphi;
0956         int it = hiHit.it;
0957         unsigned short edge = hiHit.edge;
0958         double adc = hiHit.adc;
0959         if (it > 0 && it < tbins)
0960         {
0961           if (adcval[iphi][it - 1] == 0 &&
0962               adcval[iphi][it + 1] == 0)
0963           {
0964             remove_hit(adc, iphi, it, edge, all_hit_map, adcval);
0965           }
0966         }
0967       }
0968     }
0969     */
0970     // std::cout << "done filling " << std::endl;
0971     while (all_hit_map.size() > 0)
0972     {
0973       // std::cout << "all hit map size: " << all_hit_map.size() << std::endl;
0974       auto iter = all_hit_map.rbegin();
0975       if (iter == all_hit_map.rend())
0976       {
0977         break;
0978       }
0979       ihit hiHit = iter->second;
0980       int iphi = hiHit.iphi;
0981       int it = hiHit.it;
0982       unsigned short edge = hiHit.edge;
0983       double adc = hiHit.adc;
0984       if (my_data->do_singles)
0985       {
0986         if (is_hit_isolated(iphi, it, (int) my_data->phibins, (int) my_data->tbins, adcval))
0987         {
0988           remove_hit(adc, iphi, it, edge, all_hit_map, adcval);
0989           continue;
0990         }
0991       }
0992 
0993       // put all hits in the all_hit_map (sorted by adc)
0994       // start with highest adc hit
0995       //  -> cluster around it and get vector of hits
0996       std::vector<ihit> ihit_list;
0997       int ntouch = 0;
0998       int nedge = 0;
0999       get_cluster(iphi, it, *my_data, adcval, ihit_list, ntouch, nedge);
1000 
1001       if (my_data->FixedWindow > 0)
1002       {
1003         // get cluster dimensions and check if they hit the window boundaries
1004         int wphibinhi = -1;
1005         int wphibinlo = 666666;
1006         int wtbinhi = -1;
1007         int wtbinlo = 666666;
1008         for (auto witer : ihit_list)
1009         {
1010           int wiphi = witer.iphi + my_data->phioffset;
1011           int wit = witer.it + my_data->toffset;
1012           double wadc = witer.adc;
1013           if (wadc <= 0)
1014           {
1015             continue;
1016           }
1017           if (wiphi > wphibinhi)
1018           {
1019             wphibinhi = wiphi;
1020           }
1021           if (wiphi < wphibinlo)
1022           {
1023             wphibinlo = wiphi;
1024           }
1025           if (wit > wtbinhi)
1026           {
1027             wtbinhi = wit;
1028           }
1029           if (wit < wtbinlo)
1030           {
1031             wtbinlo = wit;
1032           }
1033         }
1034         char wtsize = wtbinhi - wtbinlo + 1;
1035         char wphisize = wphibinhi - wphibinlo + 1;
1036 
1037         // check if we have a super big cluster and switch from fixed window to step down
1038         if (ihit_list.size() > (0.5 * pow((2 * my_data->FixedWindow + 1), 2)) ||
1039             wphisize >= (2 * my_data->FixedWindow + 1) ||
1040             wtsize >= (2 * my_data->FixedWindow + 1))
1041         {
1042           int window_cache = my_data->FixedWindow;
1043           // std::cout << " fixed size before " << ihit_list.size() << " | " << pow((2*my_data->FixedWindow+1),2) << std::endl;
1044           my_data->FixedWindow = 0;
1045           // reset hit list and try again without fixed window
1046           ihit_list.clear();
1047           get_cluster(iphi, it, *my_data, adcval, ihit_list, ntouch, nedge);
1048           // std::cout << " stepdown size after " << ihit_list.size() << std::endl;
1049           my_data->FixedWindow = window_cache;
1050         }
1051       }
1052       if (ihit_list.size() <= 1)
1053       {
1054         remove_hits(ihit_list, all_hit_map, adcval);
1055         ihit_list.clear();
1056         remove_hit(adc, iphi, it, edge, all_hit_map, adcval);
1057       }
1058       // -> calculate cluster parameters
1059       // -> add hits to truth association
1060       // remove hits from all_hit_map
1061       // repeat untill all_hit_map empty
1062       calc_cluster_parameter(iphi, it, ihit_list, *my_data, ntouch, nedge);
1063       remove_hits(ihit_list, all_hit_map, adcval);
1064       ihit_list.clear();
1065     }
1066     /*    if( my_data->rawhitset!=nullptr){
1067       RawHitSetv1 *hitset = my_data->rawhitset;
1068       std::cout << "Layer: " << my_data->layer
1069                 << " Side: " << my_data->side
1070                 << " Sector: " << my_data->sector
1071                 << " nhits:  " << hitset->size()
1072                 << " nhits coutn :  " << nhits
1073                 << " nclus: " << my_data->cluster_vector.size()
1074                 << std::endl;
1075     }
1076     */
1077     // pthread_exit(nullptr);
1078   }
1079   void *ProcessSector(void *threadarg)
1080   {
1081     auto my_data = static_cast<thread_data *>(threadarg);
1082     ProcessSectorData(my_data);
1083     pthread_exit(nullptr);
1084   }
1085 }  // namespace
1086 
1087 TpcClusterizer::TpcClusterizer(const std::string &name)
1088   : SubsysReco(name)
1089   , m_training(nullptr)
1090 {
1091 }
1092 
1093 bool TpcClusterizer::is_in_sector_boundary(int phibin, int sector, PHG4TpcCylinderGeom *layergeom) const
1094 {
1095   bool reject_it = false;
1096 
1097   // sector boundaries occur every 1/12 of the full phi bin range
1098   int PhiBins = layergeom->get_phibins();
1099   int PhiBinsSector = PhiBins / 12;
1100 
1101   double radius = layergeom->get_radius();
1102   double PhiBinSize = 2.0 * radius * M_PI / (double) PhiBins;
1103 
1104   // sector starts where?
1105   int sector_lo = sector * PhiBinsSector;
1106   int sector_hi = sector_lo + PhiBinsSector - 1;
1107 
1108   int sector_fiducial_bins = (int) (SectorFiducialCut / PhiBinSize);
1109 
1110   if (phibin < sector_lo + sector_fiducial_bins || phibin > sector_hi - sector_fiducial_bins)
1111   {
1112     reject_it = true;
1113     /*
1114     int layer = layergeom->get_layer();
1115     std::cout << " local maximum is in sector fiducial boundary: layer " << layer << " radius " << radius << " sector " << sector
1116     << " PhiBins " << PhiBins << " sector_fiducial_bins " << sector_fiducial_bins
1117     << " PhiBinSize " << PhiBinSize << " phibin " << phibin << " sector_lo " << sector_lo << " sector_hi " << sector_hi << std::endl;
1118     */
1119   }
1120 
1121   return reject_it;
1122 }
1123 
1124 int TpcClusterizer::InitRun(PHCompositeNode *topNode)
1125 {
1126   PHNodeIterator iter(topNode);
1127 
1128   // Looking for the DST node
1129   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
1130   if (!dstNode)
1131   {
1132     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
1133     return Fun4AllReturnCodes::ABORTRUN;
1134   }
1135 
1136   // Create the Cluster node if required
1137   auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
1138   if (!trkrclusters)
1139   {
1140     PHNodeIterator dstiter(dstNode);
1141     PHCompositeNode *DetNode =
1142         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
1143     if (!DetNode)
1144     {
1145       DetNode = new PHCompositeNode("TRKR");
1146       dstNode->addNode(DetNode);
1147     }
1148 
1149     trkrclusters = new TrkrClusterContainerv4;
1150     PHIODataNode<PHObject> *TrkrClusterContainerNode =
1151         new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
1152     DetNode->addNode(TrkrClusterContainerNode);
1153   }
1154 
1155   auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
1156   if (!clusterhitassoc)
1157   {
1158     PHNodeIterator dstiter(dstNode);
1159     PHCompositeNode *DetNode =
1160         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
1161     if (!DetNode)
1162     {
1163       DetNode = new PHCompositeNode("TRKR");
1164       dstNode->addNode(DetNode);
1165     }
1166 
1167     clusterhitassoc = new TrkrClusterHitAssocv3;
1168     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(clusterhitassoc, "TRKR_CLUSTERHITASSOC", "PHObject");
1169     DetNode->addNode(newNode);
1170   }
1171 
1172   auto training_container = findNode::getClass<TrainingHitsContainer>(dstNode, "TRAINING_HITSET");
1173   if (!training_container)
1174   {
1175     PHNodeIterator dstiter(dstNode);
1176     PHCompositeNode *DetNode =
1177         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
1178     if (!DetNode)
1179     {
1180       DetNode = new PHCompositeNode("TRKR");
1181       dstNode->addNode(DetNode);
1182     }
1183 
1184     training_container = new TrainingHitsContainer;
1185     PHIODataNode<PHObject> *TrainingHitsContainerNode =
1186         new PHIODataNode<PHObject>(training_container, "TRAINING_HITSET", "PHObject");
1187     DetNode->addNode(TrainingHitsContainerNode);
1188   }
1189 
1190   gen_hits = _store_hits || _use_nn;
1191   use_nn = _use_nn;
1192   if (use_nn)
1193   {
1194     const char *offline_main = std::getenv("OFFLINE_MAIN");
1195     assert(offline_main);
1196     std::string net_model = std::string(offline_main) + "/share/tpc/net_model.pt";
1197     try
1198     {
1199       // Deserialize the ScriptModule from a file using torch::jit::load()
1200       module_pos = torch::jit::load(net_model);
1201       std::cout << PHWHERE << "Load NN module: " << net_model << std::endl;
1202     }
1203     catch (const c10::Error &e)
1204     {
1205       std::cout << PHWHERE << "Error: Cannot load module " << net_model << std::endl;
1206       exit(1);
1207     }
1208   }
1209   else
1210   {
1211     std::cout << PHWHERE << "Use traditional clustering" << std::endl;
1212   }
1213 
1214   if (record_ClusHitsVerbose)
1215   {
1216     // get the node
1217     mClusHitsVerbose = findNode::getClass<ClusHitsVerbosev1>(topNode, "Trkr_SvtxClusHitsVerbose");
1218     if (!mClusHitsVerbose)
1219     {
1220       PHNodeIterator dstiter(dstNode);
1221       auto DetNode = dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
1222       if (!DetNode)
1223       {
1224         DetNode = new PHCompositeNode("TRKR");
1225         dstNode->addNode(DetNode);
1226       }
1227       mClusHitsVerbose = new ClusHitsVerbosev1();
1228       auto newNode = new PHIODataNode<PHObject>(mClusHitsVerbose, "Trkr_SvtxClusHitsVerbose", "PHObject");
1229       DetNode->addNode(newNode);
1230     }
1231   }
1232   auto geom =
1233       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
1234   if (!geom)
1235   {
1236     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
1237     return Fun4AllReturnCodes::ABORTRUN;
1238   }
1239   
1240   AdcClockPeriod = geom->GetFirstLayerCellGeom()->get_zstep();
1241 
1242   return Fun4AllReturnCodes::EVENT_OK;
1243 }
1244 
1245 int TpcClusterizer::process_event(PHCompositeNode *topNode)
1246 {
1247   // The TPC is the only subsystem that clusters in global coordinates. For consistency,
1248   // we must use the construction transforms to get the local coordinates.
1249   // Set the flag to use ideal transforms for the duration of this process_event, for thread safety
1250   alignmentTransformationContainer::use_alignment = false;
1251 
1252   //  int print_layer = 18;
1253 
1254   if (Verbosity() > 1000)
1255   {
1256     std::cout << "TpcClusterizer::Process_Event" << std::endl;
1257   }
1258 
1259   PHNodeIterator iter(topNode);
1260   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
1261   if (!dstNode)
1262   {
1263     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
1264     return Fun4AllReturnCodes::ABORTRUN;
1265   }
1266   if (!do_read_raw)
1267   {
1268     // get node containing the digitized hits
1269     m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1270     if (!m_hits)
1271     {
1272       std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
1273       return Fun4AllReturnCodes::ABORTRUN;
1274     }
1275   }
1276   else
1277   {
1278     // get node containing the digitized hits
1279     m_rawhits = findNode::getClass<RawHitSetContainer>(topNode, "TRKR_RAWHITSET");
1280     if (!m_rawhits)
1281     {
1282       std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
1283       return Fun4AllReturnCodes::ABORTRUN;
1284     }
1285   }
1286 
1287   // get laser event info, if exists and event has laser and rejection is on, don't bother with clustering
1288   LaserEventInfo *laserInfo = findNode::getClass<LaserEventInfo>(topNode, "LaserEventInfo");
1289   if (m_rejectEvent && laserInfo && laserInfo->isLaserEvent())
1290   {
1291     return Fun4AllReturnCodes::EVENT_OK;
1292   }
1293 
1294   // get node for clusters
1295   m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1296   if (!m_clusterlist)
1297   {
1298     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << std::endl;
1299     return Fun4AllReturnCodes::ABORTRUN;
1300   }
1301 
1302   // get node for cluster hit associations
1303   m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
1304   if (!m_clusterhitassoc)
1305   {
1306     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << std::endl;
1307     return Fun4AllReturnCodes::ABORTRUN;
1308   }
1309 
1310   // get node for training hits
1311   m_training = findNode::getClass<TrainingHitsContainer>(topNode, "TRAINING_HITSET");
1312   if (!m_training)
1313   {
1314     std::cout << PHWHERE << " ERROR: Can't find TRAINING_HITSET." << std::endl;
1315     return Fun4AllReturnCodes::ABORTRUN;
1316   }
1317 
1318   PHG4TpcCylinderGeomContainer *geom_container =
1319       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
1320   if (!geom_container)
1321   {
1322     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
1323     return Fun4AllReturnCodes::ABORTRUN;
1324   }
1325 
1326   m_tGeometry = findNode::getClass<ActsGeometry>(topNode,
1327                                                  "ActsGeometry");
1328   if (!m_tGeometry)
1329   {
1330     std::cout << PHWHERE
1331               << "ActsGeometry not found on node tree. Exiting"
1332               << std::endl;
1333     return Fun4AllReturnCodes::ABORTRUN;
1334   }
1335 
1336   // 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
1337   // The TPC clustering is more complicated than for the silicon, because we have to deal with overlapping clusters
1338 
1339   TrkrHitSetContainer::ConstRange hitsetrange;
1340   RawHitSetContainer::ConstRange rawhitsetrange;
1341   int num_hitsets = 0;
1342 
1343   if (!do_read_raw)
1344   {
1345     hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
1346     num_hitsets = std::distance(hitsetrange.first, hitsetrange.second);
1347   }
1348   else
1349   {
1350     rawhitsetrange = m_rawhits->getHitSets(TrkrDefs::TrkrId::tpcId);
1351     num_hitsets = std::distance(rawhitsetrange.first, rawhitsetrange.second);
1352   }
1353 
1354   // create structure to store given thread and associated data
1355   struct thread_pair_t
1356   {
1357     pthread_t thread{};
1358     thread_data data;
1359   };
1360 
1361   // create vector of thread pairs and reserve the right size upfront to avoid reallocation
1362   std::vector<thread_pair_t> threads;
1363   threads.reserve(num_hitsets);
1364 
1365   pthread_attr_t attr;
1366   pthread_attr_init(&attr);
1367   pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
1368 
1369   if (pthread_mutex_init(&mythreadlock, nullptr) != 0)
1370   {
1371     std::cout << std::endl
1372               << " mutex init failed" << std::endl;
1373     return 1;
1374   }
1375 //  int count = 0;
1376 
1377   if (!do_read_raw)
1378   {
1379     for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
1380          hitsetitr != hitsetrange.second;
1381          ++hitsetitr)
1382     {
1383       // if(count>0)continue;
1384       TrkrHitSet *hitset = hitsetitr->second;
1385       unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
1386       int side = TpcDefs::getSide(hitsetitr->first);
1387       unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
1388       PHG4TpcCylinderGeom *layergeom = geom_container->GetLayerCellGeom(layer);
1389 
1390       // instanciate new thread pair, at the end of thread vector
1391       thread_pair_t &thread_pair = threads.emplace_back();
1392       if (mClusHitsVerbose)
1393       {
1394         thread_pair.data.fillClusHitsVerbose = true;
1395       };
1396 
1397       thread_pair.data.layergeom = layergeom;
1398       thread_pair.data.hitset = hitset;
1399       thread_pair.data.rawhitset = nullptr;
1400       thread_pair.data.layer = layer;
1401       thread_pair.data.pedestal = pedestal;
1402       thread_pair.data.seed_threshold = seed_threshold;
1403       thread_pair.data.edge_threshold = edge_threshold;
1404       thread_pair.data.sector = sector;
1405       thread_pair.data.side = side;
1406       thread_pair.data.do_assoc = do_hit_assoc;
1407       thread_pair.data.do_wedge_emulation = do_wedge_emulation;
1408       thread_pair.data.do_singles = do_singles;
1409       thread_pair.data.tGeometry = m_tGeometry;
1410       thread_pair.data.maxHalfSizeT = MaxClusterHalfSizeT;
1411       thread_pair.data.maxHalfSizePhi = MaxClusterHalfSizePhi;
1412       thread_pair.data.sampa_tbias = m_sampa_tbias;
1413       thread_pair.data.verbosity = Verbosity();
1414       thread_pair.data.do_split = do_split;
1415       thread_pair.data.FixedWindow = do_fixed_window;
1416       thread_pair.data.min_err_squared = min_err_squared;
1417       thread_pair.data.min_clus_size = min_clus_size;
1418       thread_pair.data.min_adc_sum = min_adc_sum;
1419       unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
1420       unsigned short NPhiBinsSector = NPhiBins / 12;
1421       unsigned short NTBins = 0;
1422       if (is_reco)
1423       {
1424         NTBins = NZBinsSide;
1425       }
1426       else
1427       {
1428         NTBins = (unsigned short) layergeom->get_zbins();
1429       }
1430       unsigned short NTBinsSide = NTBins;
1431       unsigned short NTBinsMin = 0;
1432       unsigned short PhiOffset = NPhiBinsSector * sector;
1433       unsigned short TOffset = NTBinsMin;
1434 
1435       m_tdriftmax = AdcClockPeriod * NZBinsSide;
1436       thread_pair.data.m_tdriftmax = m_tdriftmax;
1437 
1438       thread_pair.data.phibins = NPhiBinsSector;
1439       thread_pair.data.phioffset = PhiOffset;
1440       thread_pair.data.tbins = NTBinsSide;
1441       thread_pair.data.toffset = TOffset;
1442 
1443       thread_pair.data.radius = layergeom->get_radius();
1444       thread_pair.data.drift_velocity = m_tGeometry->get_drift_velocity();
1445       thread_pair.data.pads_per_sector = 0;
1446       thread_pair.data.phistep = 0;
1447       int rc;
1448       rc = pthread_create(&thread_pair.thread, &attr, ProcessSector, (void *) &thread_pair.data);
1449 
1450       if (rc)
1451       {
1452         std::cout << "Error:unable to create thread," << rc << std::endl;
1453       }
1454       if (do_sequential)
1455       {
1456         int rc2 = pthread_join(thread_pair.thread, nullptr);
1457         if (rc2)
1458         {
1459           std::cout << "Error:unable to join," << rc2 << std::endl;
1460         }
1461 
1462         // get the hitsetkey from thread data
1463         const auto &data(thread_pair.data);
1464         const auto hitsetkey = TpcDefs::genHitSetKey(data.layer, data.sector, data.side);
1465 
1466         // copy clusters to map
1467         for (uint32_t index = 0; index < data.cluster_vector.size(); ++index)
1468         {
1469           // generate cluster key
1470           const auto ckey = TrkrDefs::genClusKey(hitsetkey, index);
1471 
1472           // get cluster
1473           auto cluster = data.cluster_vector[index];
1474 
1475           // insert in map
1476           m_clusterlist->addClusterSpecifyKey(ckey, cluster);
1477 
1478           if (mClusHitsVerbose)
1479           {
1480             for (auto &hit : data.phivec_ClusHitsVerbose[index])
1481             {
1482               mClusHitsVerbose->addPhiHit(hit.first, hit.second);
1483             }
1484             for (auto &hit : data.zvec_ClusHitsVerbose[index])
1485             {
1486               mClusHitsVerbose->addZHit(hit.first, hit.second);
1487             }
1488             mClusHitsVerbose->push_hits(ckey);
1489           }
1490         }
1491 
1492         // copy hit associations to map
1493         for (const auto &[index, hkey] : thread_pair.data.association_vector)
1494         {
1495           // generate cluster key
1496           const auto ckey = TrkrDefs::genClusKey(hitsetkey, index);
1497 
1498           // add to association table
1499           m_clusterhitassoc->addAssoc(ckey, hkey);
1500         }
1501       }
1502 //      count++;
1503     }
1504   }
1505   else
1506   {
1507     for (RawHitSetContainer::ConstIterator hitsetitr = rawhitsetrange.first;
1508          hitsetitr != rawhitsetrange.second;
1509          ++hitsetitr)
1510     {
1511       //    if(count>0)continue;
1512       //    const auto hitsetid = hitsetitr->first;
1513       //    std::cout << " starting thread # " << count << std::endl;
1514 
1515       RawHitSet *hitset = hitsetitr->second;
1516       unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
1517       int side = TpcDefs::getSide(hitsetitr->first);
1518       unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
1519       PHG4TpcCylinderGeom *layergeom = geom_container->GetLayerCellGeom(layer);
1520 
1521       // instanciate new thread pair, at the end of thread vector
1522       thread_pair_t &thread_pair = threads.emplace_back();
1523 
1524       thread_pair.data.layergeom = layergeom;
1525       thread_pair.data.hitset = nullptr;
1526       thread_pair.data.rawhitset = hitset;
1527       thread_pair.data.layer = layer;
1528       thread_pair.data.pedestal = pedestal;
1529       thread_pair.data.sector = sector;
1530       thread_pair.data.side = side;
1531       thread_pair.data.do_assoc = do_hit_assoc;
1532       thread_pair.data.do_wedge_emulation = do_wedge_emulation;
1533       thread_pair.data.tGeometry = m_tGeometry;
1534       thread_pair.data.maxHalfSizeT = MaxClusterHalfSizeT;
1535       thread_pair.data.maxHalfSizePhi = MaxClusterHalfSizePhi;
1536       thread_pair.data.sampa_tbias = m_sampa_tbias;
1537       thread_pair.data.verbosity = Verbosity();
1538 
1539       unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
1540       unsigned short NPhiBinsSector = NPhiBins / 12;
1541       unsigned short NTBins = (unsigned short) layergeom->get_zbins();
1542       unsigned short NTBinsSide = NTBins;
1543       unsigned short NTBinsMin = 0;
1544       unsigned short PhiOffset = NPhiBinsSector * sector;
1545       unsigned short TOffset = NTBinsMin;
1546 
1547       m_tdriftmax = AdcClockPeriod * NZBinsSide;
1548       thread_pair.data.m_tdriftmax = m_tdriftmax;
1549 
1550       thread_pair.data.phibins = NPhiBinsSector;
1551       thread_pair.data.phioffset = PhiOffset;
1552       thread_pair.data.tbins = NTBinsSide;
1553       thread_pair.data.toffset = TOffset;
1554 
1555       /*
1556       PHG4TpcCylinderGeom *testlayergeom = geom_container->GetLayerCellGeom(32);
1557       for( float iphi = 1408; iphi < 1408+ 128;iphi+=0.1){
1558         double clusiphi = iphi;
1559         double clusphi = testlayergeom->get_phi(clusiphi);
1560         double radius = layergeom->get_radius();
1561         float clusx = radius * cos(clusphi);
1562         float clusy = radius * sin(clusphi);
1563         float clusz  = -37.524;
1564 
1565         TrkrDefs::hitsetkey tpcHitSetKey = TpcDefs::genHitSetKey( 32,11, 0 );
1566         Acts::Vector3 global(clusx, clusy, clusz);
1567         TrkrDefs::subsurfkey subsurfkey = 0;
1568 
1569         Surface surface = m_tGeometry->get_tpc_surface_from_coords(
1570                                                                    tpcHitSetKey,
1571                                                                    global,
1572                                                                    subsurfkey);
1573         std::cout << " iphi: " << iphi << " clusphi: " << clusphi << " surfkey " << subsurfkey << std::endl;
1574         //  std::cout << "surfkey" << subsurfkey << std::endl;
1575       }
1576       continue;
1577       */
1578       int rc = 0;
1579       //      if(layer==32)
1580       rc = pthread_create(&thread_pair.thread, &attr, ProcessSector, (void *) &thread_pair.data);
1581       //      else
1582       // continue;
1583 
1584       if (rc)
1585       {
1586         std::cout << "Error:unable to create thread," << rc << std::endl;
1587       }
1588 
1589       if (do_sequential)
1590       {
1591         int rc2 = pthread_join(thread_pair.thread, nullptr);
1592         if (rc2)
1593         {
1594           std::cout << "Error:unable to join," << rc2 << std::endl;
1595         }
1596 
1597         // get the hitsetkey from thread data
1598         const auto &data(thread_pair.data);
1599         const auto hitsetkey = TpcDefs::genHitSetKey(data.layer, data.sector, data.side);
1600 
1601         // copy clusters to map
1602         for (uint32_t index = 0; index < data.cluster_vector.size(); ++index)
1603         {
1604           // generate cluster key
1605           const auto ckey = TrkrDefs::genClusKey(hitsetkey, index);
1606 
1607           // get cluster
1608           auto cluster = data.cluster_vector[index];
1609 
1610           // insert in map
1611           m_clusterlist->addClusterSpecifyKey(ckey, cluster);
1612         }
1613 
1614         // copy hit associations to map
1615         for (const auto &[index, hkey] : thread_pair.data.association_vector)
1616         {
1617           // generate cluster key
1618           const auto ckey = TrkrDefs::genClusKey(hitsetkey, index);
1619 
1620           // add to association table
1621           m_clusterhitassoc->addAssoc(ckey, hkey);
1622         }
1623       }
1624 //      count++;
1625     }
1626   }
1627 
1628   pthread_attr_destroy(&attr);
1629 //  count = 0;
1630   // wait for completion of all threads
1631   if (!do_sequential)
1632   {
1633     for (const auto &thread_pair : threads)
1634     {
1635       int rc2 = pthread_join(thread_pair.thread, nullptr);
1636       if (rc2)
1637       {
1638         std::cout << "Error:unable to join," << rc2 << std::endl;
1639       }
1640 
1641       // get the hitsetkey from thread data
1642       const auto &data(thread_pair.data);
1643       const auto hitsetkey = TpcDefs::genHitSetKey(data.layer, data.sector, data.side);
1644 
1645       // copy clusters to map
1646       for (uint32_t index = 0; index < data.cluster_vector.size(); ++index)
1647       {
1648         // generate cluster key
1649         const auto ckey = TrkrDefs::genClusKey(hitsetkey, index);
1650 
1651         // get cluster
1652         auto cluster = data.cluster_vector[index];
1653 
1654         // insert in map
1655         // std::cout << "X: " << cluster->getLocalX() << "Y: " << cluster->getLocalY() << std::endl;
1656         m_clusterlist->addClusterSpecifyKey(ckey, cluster);
1657 
1658         if (mClusHitsVerbose)
1659         {
1660           for (auto &hit : data.phivec_ClusHitsVerbose[index])
1661           {
1662             mClusHitsVerbose->addPhiHit(hit.first, (float) hit.second);
1663           }
1664           for (auto &hit : data.zvec_ClusHitsVerbose[index])
1665           {
1666             mClusHitsVerbose->addZHit(hit.first, (float) hit.second);
1667           }
1668           mClusHitsVerbose->push_hits(ckey);
1669         }
1670       }
1671 
1672       // copy hit associations to map
1673       for (const auto &[index, hkey] : thread_pair.data.association_vector)
1674       {
1675         // generate cluster key
1676         const auto ckey = TrkrDefs::genClusKey(hitsetkey, index);
1677 
1678         // add to association table
1679         m_clusterhitassoc->addAssoc(ckey, hkey);
1680       }
1681 
1682       for (auto v_hit : thread_pair.data.v_hits)
1683       {
1684         if (_store_hits)
1685         {
1686           m_training->v_hits.emplace_back(*v_hit);
1687         }
1688         delete v_hit;
1689       }
1690     }
1691   }
1692 
1693   // set the flag to use alignment transformations, needed by the rest of reconstruction
1694   alignmentTransformationContainer::use_alignment = true;
1695 
1696   if (Verbosity() > 0)
1697   {
1698     std::cout << "TPC Clusterizer found " << m_clusterlist->size() << " Clusters " << std::endl;
1699   }
1700   return Fun4AllReturnCodes::EVENT_OK;
1701 }
1702 
1703 int TpcClusterizer::End(PHCompositeNode * /*topNode*/)
1704 {
1705   return Fun4AllReturnCodes::EVENT_OK;
1706 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!