packages/tpc/TpcSimpleClusterizer.cc

0001 #include "TpcSimpleClusterizer.h"
0002
0003 #include <trackbase/TpcDefs.h>
0004
0005 #include <trackbase/TrkrClusterContainerv4.h>
0006 #include <trackbase/TrkrClusterHitAssocv3.h>
0007 #include <trackbase/TrkrClusterv3.h>
0008 #include <trackbase/TrkrDefs.h>  // for hitkey, getLayer
0009 #include <trackbase/TrkrHitSet.h>
0010 #include <trackbase/TrkrHitSetContainer.h>
0011 #include <trackbase/TrkrHitv2.h>
0012
0013 #include <fun4all/Fun4AllReturnCodes.h>
0014 #include <fun4all/SubsysReco.h>  // for SubsysReco
0015
0016 #include <g4detectors/PHG4TpcCylinderGeom.h>
0017 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0018
0019 #include <Acts/Definitions/Units.hpp>
0020 #include <Acts/Surfaces/Surface.hpp>
0021
0022 #include <phool/PHCompositeNode.h>
0023 #include <phool/PHIODataNode.h>  // for PHIODataNode
0024 #include <phool/PHNode.h>        // for PHNode
0025 #include <phool/PHNodeIterator.h>
0026 #include <phool/PHObject.h>  // for PHObject
0027 #include <phool/getClass.h>
0028 #include <phool/phool.h>  // for PHWHERE
0029
0030 #include <TMatrixFfwd.h>    // for TMatrixF
0031 #include <TMatrixT.h>       // for TMatrixT, ope...
0032 #include <TMatrixTUtils.h>  // for TMatrixTRow
0033
0034 #include <TFile.h>
0035
0036 #include <array>
0037 #include <cmath>  // for sqrt, cos, sin
0038 #include <iostream>
0039 #include <map>  // for _Rb_tree_cons...
0040 #include <string>
0041 #include <utility>  // for pair
0042 #include <vector>
0043 // Terra incognita....
0044 #include <pthread.h>
0045
0046 namespace
0047 {
0048   template <class T>
0049   inline constexpr T square(const T &x)
0050   {
0051     return x * x;
0052   }
0053
0054   using iphiz = std::pair<unsigned short, unsigned short>;
0055   using ihit = std::pair<unsigned short, iphiz>;
0056   using assoc = std::pair<uint32_t, TrkrDefs::hitkey>;
0057
0058   struct thread_data
0059   {
0060     PHG4TpcCylinderGeom *layergeom = nullptr;
0061     TrkrHitSet *hitset = nullptr;
0062     ActsGeometry *tGeometry = nullptr;
0063     unsigned int layer = 0;
0064     int side = 0;
0065     unsigned int sector = 0;
0066     float pedestal = 0;
0067     bool do_assoc = true;
0068     unsigned short phibins = 0;
0069     unsigned short phioffset = 0;
0070     unsigned short zbins = 0;
0071     unsigned short zoffset = 0;
0072     double par0_neg = 0;
0073     double par0_pos = 0;
0074     std::vector<assoc> association_vector;
0075     std::vector<TrkrCluster *> cluster_vector;
0076   };
0077
0078   pthread_mutex_t mythreadlock;
0079
0080   void remove_hit(double adc, int phibin, int zbin, std::multimap<unsigned short, ihit> &all_hit_map, std::vector<std::vector<unsigned short>> &adcval)
0081   {
0082     using hit_iterator = std::multimap<unsigned short, ihit>::iterator;
0083     std::pair<hit_iterator, hit_iterator> iterpair = all_hit_map.equal_range(adc);
0084     hit_iterator it = iterpair.first;
0085     for (; it != iterpair.second; ++it)
0086     {
0087       if (it->second.second.first == phibin && it->second.second.second == zbin)
0088       {
0089         all_hit_map.erase(it);
0090         break;
0091       }
0092     }
0093     adcval[phibin][zbin] = 0;
0094   }
0095
0096   void remove_hits(std::vector<ihit> &ihit_list, std::multimap<unsigned short, ihit> &all_hit_map, std::vector<std::vector<unsigned short>> &adcval)
0097   {
0098     for (auto &iter : ihit_list)
0099     {
0100       unsigned short adc = iter.first;
0101       unsigned short phibin = iter.second.first;
0102       unsigned short zbin = iter.second.second;
0103       remove_hit(adc, phibin, zbin, all_hit_map, adcval);
0104     }
0105   }
0106
0107   void get_cluster(int phibin, int zbin, const std::vector<std::vector<unsigned short>> &adcval, std::vector<ihit> &ihit_list)
0108   {
0109     // on hit = one cluster
0110     const int &iphi = phibin;
0111     const int &iz = zbin;
0112     iphiz iCoord(std::make_pair(iphi, iz));
0113     ihit thisHit(adcval[iphi][iz], iCoord);
0114     ihit_list.push_back(thisHit);
0115   }
0116
0117   void calc_cluster_parameter(std::vector<ihit> &ihit_list, thread_data &my_data)
0118   {
0119     // loop over the hits in this cluster
0120     double z_sum = 0.0;
0121     double phi_sum = 0.0;
0122     double adc_sum = 0.0;
0123     double z2_sum = 0.0;
0124     double phi2_sum = 0.0;
0125
0126     double radius = my_data.layergeom->get_radius();  // returns center of layer
0127
0128     int phibinhi = -1;
0129     int phibinlo = 666666;
0130     int zbinhi = -1;
0131     int zbinlo = 666666;
0132
0133     std::vector<TrkrDefs::hitkey> hitkeyvec;
0134     for (auto &iter : ihit_list)
0135     {
0136       double adc = iter.first;
0137
0138       if (adc <= 0)
0139       {
0140         continue;
0141       }
0142
0143       int iphi = iter.second.first + my_data.phioffset;
0144       int iz = iter.second.second + my_data.zoffset;
0145       if (iphi > phibinhi)
0146       {
0147         phibinhi = iphi;
0148       }
0149       if (iphi < phibinlo)
0150       {
0151         phibinlo = iphi;
0152       }
0153       if (iz > zbinhi)
0154       {
0155         zbinhi = iz;
0156       }
0157       if (iz < zbinlo)
0158       {
0159         zbinlo = iz;
0160       }
0161
0162       // update phi sums
0163       double phi_center = my_data.layergeom->get_phicenter(iphi, my_data.side);
0164       phi_sum += phi_center * adc;
0165       phi2_sum += square(phi_center) * adc;
0166
0167       // update z sums
0168       double z = my_data.layergeom->get_zcenter(iz);
0169       z_sum += z * adc;
0170       z2_sum += square(z) * adc;
0171
0172       adc_sum += adc;
0173
0174       // capture the hitkeys for all adc values above a certain threshold
0175       TrkrDefs::hitkey hitkey = TpcDefs::genHitKey(iphi, iz);
0176       // if(adc>5)
0177       hitkeyvec.push_back(hitkey);
0178     }
0179     // if (adc_sum < 10) return;  // skip obvious noise "clusters"
0180
0181     // This is the global position
0182     double clusphi = phi_sum / adc_sum;
0183     double clusz = z_sum / adc_sum;
0184     const double clusx = radius * std::cos(clusphi);
0185     const double clusy = radius * std::sin(clusphi);
0186
0187     const double phi_cov = phi2_sum / adc_sum - square(clusphi);
0188     const double z_cov = z2_sum / adc_sum - square(clusz);
0189
0190     // create the cluster entry directly in the node tree
0191     // Estimate the errors
0192     const double phi_err_square = (phibinhi == phibinlo) ? square(radius * my_data.layergeom->get_phistep()) / 12 : square(radius) * phi_cov / (adc_sum * 0.14);
0193
0194     const double z_err_square = (zbinhi == zbinlo) ? square(my_data.layergeom->get_zstep()) / 12 : z_cov / (adc_sum * 0.14);
0195
0196     // 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:
0197     // 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
0198     //    - N is the number of electrons that drift to the readout plane
0199     // We have to convert (sum of adc units for all bins in the cluster) to number of ionization electrons N
0200     // 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
0201     // To get equivalent charge per Z bin, so that summing ADC input voltage over all Z bins returns total input charge, divide voltages by 2.4 for 80 ns SAMPA
0202     // Equivalent charge per Z bin is then  (ADU x 2200 mV / 1024) / 2.4 x (1/20) fC/mV x (1/1.6e-04) electrons/fC x (1/2000) = ADU x 0.14
0203
0204     // cluster z correction
0205     clusz -= (clusz < 0) ? my_data.par0_neg : my_data.par0_pos;
0206
0207     // create cluster and fill
0208     auto clus = new TrkrClusterv3;
0209     clus->setAdc(adc_sum);
0210
0211     /// Get the surface key to find the surface from the map
0212     TrkrDefs::hitsetkey tpcHitSetKey = TpcDefs::genHitSetKey(my_data.layer, my_data.sector, my_data.side);
0213
0214     Acts::Vector3 global(clusx, clusy, clusz);
0215
0216     TrkrDefs::subsurfkey subsurfkey;
0217     Surface surface = my_data.tGeometry->get_tpc_surface_from_coords(
0218         tpcHitSetKey,
0219         global,
0220         subsurfkey);
0221
0222     if (!surface)
0223     {
0224       /// If the surface can't be found, we can't track with it. So
0225       /// just return and don't add the cluster to the container
0226       return;
0227     }
0228
0229     clus->setSubSurfKey(subsurfkey);
0230
0231     Acts::Vector3 center = surface->center(my_data.tGeometry->geometry().getGeoContext()) / Acts::UnitConstants::cm;
0232
0233     /// no conversion needed, only used in acts
0234     const Acts::Vector3 normal = surface->normal(my_data.tGeometry->geometry().getGeoContext(), Acts::Vector3(1,1,1), Acts::Vector3(1,1,1));
0235     const double clusRadius = std::sqrt(square(clusx) + square(clusy));
0236     const double rClusPhi = clusRadius * clusphi;
0237     const double surfRadius = sqrt(center(0) * center(0) + center(1) * center(1));
0238     const double surfPhiCenter = atan2(center[1], center[0]);
0239     const double surfRphiCenter = surfPhiCenter * surfRadius;
0240     const double surfZCenter = center[2];
0241     auto local = surface->globalToLocal(my_data.tGeometry->geometry().getGeoContext(),
0242                                         global * Acts::UnitConstants::cm,
0243                                         normal);
0244     Acts::Vector2 localPos;
0245
0246     /// Prefer Acts transformation since we build the TPC surfaces manually
0247     if (local.ok())
0248     {
0249       localPos = local.value() / Acts::UnitConstants::cm;
0250     }
0251     else
0252     {
0253       /// otherwise take the manual calculation
0254       localPos(0) = rClusPhi - surfRphiCenter;
0255       localPos(1) = clusz - surfZCenter;
0256     }
0257
0258     clus->setLocalX(localPos(0));
0259     clus->setLocalY(localPos(1));
0260     clus->setActsLocalError(0, 0, phi_err_square);
0261     clus->setActsLocalError(1, 0, 0);
0262     clus->setActsLocalError(0, 1, 0);
0263     clus->setActsLocalError(1, 1, z_err_square);
0264
0265     my_data.cluster_vector.push_back(clus);
0266
0267     // Add the hit associations to the TrkrClusterHitAssoc node
0268     // we need the cluster key and all associated hit keys (note: the cluster key includes the hitset key)
0269
0270     if (my_data.do_assoc)
0271     {
0272       // get cluster index in vector. It is used to store associations, and build relevant cluster keys when filling the containers
0273       uint32_t index = my_data.cluster_vector.size() - 1;
0274       for (unsigned int &i : hitkeyvec)
0275       {
0276         my_data.association_vector.emplace_back(index, i);
0277       }
0278     }
0279   }
0280
0281   void *ProcessSector(void *threadarg)
0282   {
0283     auto my_data = (struct thread_data *) threadarg;
0284
0285     const auto &pedestal = my_data->pedestal;
0286     const auto &phibins = my_data->phibins;
0287     const auto &phioffset = my_data->phioffset;
0288     const auto &zbins = my_data->zbins;
0289     const auto &zoffset = my_data->zoffset;
0290
0291     TrkrHitSet *hitset = my_data->hitset;
0292     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0293
0294     // for convenience, create a 2D vector to store adc values in and initialize to zero
0295     std::vector<std::vector<unsigned short>> adcval(phibins, std::vector<unsigned short>(zbins, 0));
0296     std::multimap<unsigned short, ihit> all_hit_map;
0297     std::vector<ihit> hit_vect;
0298
0299     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0300          hitr != hitrangei.second;
0301          ++hitr)
0302     {
0303       unsigned short phibin = TpcDefs::getPad(hitr->first) - phioffset;
0304       unsigned short zbin = TpcDefs::getTBin(hitr->first) - zoffset;
0305
0306       float_t fadc = (hitr->second->getAdc()) - pedestal;  // proper int rounding +0.5
0307       // std::cout << " layer: " << my_data->layer  << " phibin " << phibin << " zbin " << zbin << " fadc " << hitr->second->getAdc() << " pedestal " << pedestal << " fadc " << std::endl
0308
0309       unsigned short adc = 0;
0310       if (fadc > 0)
0311       {
0312         adc = (unsigned short) fadc;
0313       }
0314
0315       //     if(phibin < 0) continue; // phibin is unsigned int, <0 cannot happen
0316       if (phibin >= phibins)
0317       {
0318         continue;
0319       }
0320       //     if(zbin   < 0) continue;
0321       if (zbin >= zbins)
0322       {
0323         continue;  // zbin is unsigned int, <0 cannot happen
0324       }
0325
0326       if (adc > 0)
0327       {
0328         iphiz iCoord(std::make_pair(phibin, zbin));
0329         ihit thisHit(adc, iCoord);
0330         if (adc > 5)
0331         {
0332           all_hit_map.insert(std::make_pair(adc, thisHit));
0333         }
0334         // adcval[phibin][zbin] = (unsigned short) adc;
0335         adcval[phibin][zbin] = (unsigned short) adc;
0336       }
0337     }
0338
0339     while (all_hit_map.size() > 0)
0340     {
0341       auto iter = all_hit_map.rbegin();
0342       if (iter == all_hit_map.rend())
0343       {
0344         break;
0345       }
0346       ihit hiHit = iter->second;
0347       int iphi = hiHit.second.first;
0348       int iz = hiHit.second.second;
0349
0350       // put all hits in the all_hit_map (sorted by adc)
0351       // start with highest adc hit
0352       //  -> cluster around it and get vector of hits
0353       std::vector<ihit> ihit_list;
0354       get_cluster(iphi, iz, adcval, ihit_list);
0355
0356       // -> calculate cluster parameters
0357       // -> add hits to truth association
0358       // remove hits from all_hit_map
0359       // repeat untill all_hit_map empty
0360       calc_cluster_parameter(ihit_list, *my_data);
0361       remove_hits(ihit_list, all_hit_map, adcval);
0362     }
0363     pthread_exit(nullptr);
0364   }
0365 }  // namespace
0366
0367 TpcSimpleClusterizer::TpcSimpleClusterizer(const std::string &name)
0368   : SubsysReco(name)
0369 {
0370 }
0371
0372 bool TpcSimpleClusterizer::is_in_sector_boundary(int phibin, int sector, PHG4TpcCylinderGeom *layergeom) const
0373 {
0374   bool reject_it = false;
0375
0376   // sector boundaries occur every 1/12 of the full phi bin range
0377   int PhiBins = layergeom->get_phibins();
0378   int PhiBinsSector = PhiBins / 12;
0379
0380   double radius = layergeom->get_radius();
0381   double PhiBinSize = 2.0 * radius * M_PI / (double) PhiBins;
0382
0383   // sector starts where?
0384   int sector_lo = sector * PhiBinsSector;
0385   int sector_hi = sector_lo + PhiBinsSector - 1;
0386
0387   int sector_fiducial_bins = (int) (SectorFiducialCut / PhiBinSize);
0388
0389   if (phibin < sector_lo + sector_fiducial_bins || phibin > sector_hi - sector_fiducial_bins)
0390   {
0391     reject_it = true;
0392     /*
0393     int layer = layergeom->get_layer();
0394     std::cout << " local maximum is in sector fiducial boundary: layer " << layer << " radius " << radius << " sector " << sector
0395     << " PhiBins " << PhiBins << " sector_fiducial_bins " << sector_fiducial_bins
0396     << " PhiBinSize " << PhiBinSize << " phibin " << phibin << " sector_lo " << sector_lo << " sector_hi " << sector_hi << std::endl;
0397     */
0398   }
0399
0400   return reject_it;
0401 }
0402
0403 int TpcSimpleClusterizer::InitRun(PHCompositeNode *topNode)
0404 {
0405   PHNodeIterator iter(topNode);
0406
0407   // Looking for the DST node
0408   PHCompositeNode *dstNode = dynamic_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0409   if (!dstNode)
0410   {
0411     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0412     return Fun4AllReturnCodes::ABORTRUN;
0413   }
0414
0415   // Create the Cluster node if required
0416   auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
0417   if (!trkrclusters)
0418   {
0419     PHNodeIterator dstiter(dstNode);
0420     PHCompositeNode *DetNode =
0421         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0422     if (!DetNode)
0423     {
0424       DetNode = new PHCompositeNode("TRKR");
0425       dstNode->addNode(DetNode);
0426     }
0427
0428     trkrclusters = new TrkrClusterContainerv4;
0429     PHIODataNode<PHObject> *TrkrClusterContainerNode =
0430         new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
0431     DetNode->addNode(TrkrClusterContainerNode);
0432   }
0433
0434   auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0435   if (!clusterhitassoc)
0436   {
0437     PHNodeIterator dstiter(dstNode);
0438     PHCompositeNode *DetNode =
0439         dynamic_cast<PHCompositeNode *>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0440     if (!DetNode)
0441     {
0442       DetNode = new PHCompositeNode("TRKR");
0443       dstNode->addNode(DetNode);
0444     }
0445
0446     clusterhitassoc = new TrkrClusterHitAssocv3;
0447     PHIODataNode<PHObject> *newNode = new PHIODataNode<PHObject>(clusterhitassoc, "TRKR_CLUSTERHITASSOC", "PHObject");
0448     DetNode->addNode(newNode);
0449   }
0450
0451   return Fun4AllReturnCodes::EVENT_OK;
0452 }
0453
0454 int TpcSimpleClusterizer::process_event(PHCompositeNode *topNode)
0455 {
0456   //  int print_layer = 18;
0457
0458   if (Verbosity() > 1000)
0459   {
0460     std::cout << "TpcSimpleClusterizer::Process_Event" << std::endl;
0461   }
0462
0463   PHNodeIterator iter(topNode);
0464   PHCompositeNode *dstNode = static_cast<PHCompositeNode *>(iter.findFirst("PHCompositeNode", "DST"));
0465   if (!dstNode)
0466   {
0467     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0468     return Fun4AllReturnCodes::ABORTRUN;
0469   }
0470
0471   // get node containing the digitized hits
0472   m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0473   if (!m_hits)
0474   {
0475     std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
0476     return Fun4AllReturnCodes::ABORTRUN;
0477   }
0478
0479   // get node for clusters
0480   m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0481   if (!m_clusterlist)
0482   {
0483     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << std::endl;
0484     return Fun4AllReturnCodes::ABORTRUN;
0485   }
0486
0487   // get node for cluster hit associations
0488   m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0489   if (!m_clusterhitassoc)
0490   {
0491     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << std::endl;
0492     return Fun4AllReturnCodes::ABORTRUN;
0493   }
0494
0495   PHG4TpcCylinderGeomContainer *geom_container =
0496       findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0497   if (!geom_container)
0498   {
0499     std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0500     return Fun4AllReturnCodes::ABORTRUN;
0501   }
0502
0503   m_tGeometry = findNode::getClass<ActsGeometry>(topNode,
0504                                                  "ActsGeometry");
0505   if (!m_tGeometry)
0506   {
0507     std::cout << PHWHERE
0508               << "ActsGeometry not found on node tree. Exiting"
0509               << std::endl;
0510     return Fun4AllReturnCodes::ABORTRUN;
0511   }
0512
0513   // 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
0514   // The TPC clustering is more complicated than for the silicon, because we have to deal with overlapping clusters
0515
0516   // loop over the TPC HitSet objects
0517   TrkrHitSetContainer::ConstRange hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
0518   const int num_hitsets = std::distance(hitsetrange.first, hitsetrange.second);
0519
0520   // create structure to store given thread and associated data
0521   struct thread_pair_t
0522   {
0523     pthread_t thread{};
0524     thread_data data;
0525   };
0526
0527   // create vector of thread pairs and reserve the right size upfront to avoid reallocation
0528   std::vector<thread_pair_t> threads;
0529   threads.reserve(num_hitsets);
0530
0531   pthread_attr_t attr;
0532   pthread_attr_init(&attr);
0533   pthread_attr_setdetachstate(&attr, PTHREAD_CREATE_JOINABLE);
0534
0535   if (pthread_mutex_init(&mythreadlock, nullptr) != 0)
0536   {
0537     std::cout << std::endl << " mutex init failed" << std::endl;
0538     return 1;
0539   }
0540
0541   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0542        hitsetitr != hitsetrange.second;
0543        ++hitsetitr)
0544   {
0545     TrkrHitSet *hitset = hitsetitr->second;
0546     unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
0547     int side = TpcDefs::getSide(hitsetitr->first);
0548     unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
0549     PHG4TpcCylinderGeom *layergeom = geom_container->GetLayerCellGeom(layer);
0550
0551     // instanciate new thread pair, at the end of thread vector
0552     thread_pair_t &thread_pair = threads.emplace_back();
0553
0554     thread_pair.data.layergeom = layergeom;
0555     thread_pair.data.hitset = hitset;
0556     thread_pair.data.layer = layer;
0557     thread_pair.data.pedestal = pedestal;
0558     thread_pair.data.sector = sector;
0559     thread_pair.data.side = side;
0560     thread_pair.data.do_assoc = do_hit_assoc;
0561     thread_pair.data.tGeometry = m_tGeometry;
0562     thread_pair.data.par0_neg = par0_neg;
0563     thread_pair.data.par0_pos = par0_pos;
0564
0565     unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
0566     unsigned short NPhiBinsSector = NPhiBins / 12;
0567     unsigned short NZBins = (unsigned short) layergeom->get_zbins();
0568     unsigned short NZBinsSide = NZBins / 2;
0569     unsigned short NZBinsMin = 0;
0570     unsigned short PhiOffset = NPhiBinsSector * sector;
0571
0572     if (side == 0)
0573     {
0574       NZBinsMin = 0;
0575     }
0576     else
0577     {
0578       NZBinsMin = NZBins / 2;
0579     }
0580
0581     unsigned short ZOffset = NZBinsMin;
0582
0583     thread_pair.data.phibins = NPhiBinsSector;
0584     thread_pair.data.phioffset = PhiOffset;
0585     thread_pair.data.zbins = NZBinsSide;
0586     thread_pair.data.zoffset = ZOffset;
0587
0588     int rc = pthread_create(&thread_pair.thread, &attr, ProcessSector, (void *) &thread_pair.data);
0589     if (rc)
0590     {
0591       std::cout << "Error:unable to create thread," << rc << std::endl;
0592     }
0593   }
0594
0595   pthread_attr_destroy(&attr);
0596
0597   // wait for completion of all threads
0598   for (const auto &thread_pair : threads)
0599   {
0600     int rc2 = pthread_join(thread_pair.thread, nullptr);
0601     if (rc2)
0602     {
0603       std::cout << "Error:unable to join," << rc2 << std::endl;
0604     }
0605
0606     // get the hitsetkey from thread data
0607     const auto &data(thread_pair.data);
0608     const auto hitsetkey = TpcDefs::genHitSetKey(data.layer, data.sector, data.side);
0609
0610     // copy clusters to map
0611     for (uint32_t index = 0; index < data.cluster_vector.size(); ++index)
0612     {
0613       // generate cluster key
0614       const auto ckey = TrkrDefs::genClusKey(hitsetkey, index);
0615
0616       // get cluster
0617       auto cluster = data.cluster_vector[index];
0618
0619       // insert in map
0620       m_clusterlist->addClusterSpecifyKey(ckey, cluster);
0621     }
0622
0623     // copy hit associations to map
0624     for (const auto &[index, hkey] : thread_pair.data.association_vector)
0625     {
0626       // generate cluster key
0627       const auto ckey = TrkrDefs::genClusKey(hitsetkey, index);
0628
0629       // add to association table
0630       m_clusterhitassoc->addAssoc(ckey, hkey);
0631     }
0632   }
0633
0634   if (Verbosity() > 0)
0635   {
0636     std::cout << "TPC Clusterizer found " << m_clusterlist->size() << " Clusters " << std::endl;
0637   }
0638   return Fun4AllReturnCodes::EVENT_OK;
0639 }
0640
0641 int TpcSimpleClusterizer::End(PHCompositeNode * /*topNode*/)
0642 {
0643   return Fun4AllReturnCodes::EVENT_OK;
0644 }