Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:20:23 

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