packages/intt/InttClusterizer.cc

0001 #include "InttClusterizer.h"
0002 #include "CylinderGeomIntt.h"
0003
0004 #include <trackbase/InttDefs.h>
0005 #include <trackbase/TrkrClusterContainerv4.h>
0006 #include <trackbase/TrkrClusterCrossingAssocv1.h>
0007 #include <trackbase/TrkrClusterHitAssocv3.h>
0008 #include <trackbase/TrkrClusterv5.h>
0009 #include <trackbase/TrkrHit.h>
0010 #include <trackbase/TrkrHitSet.h>
0011 #include <trackbase/TrkrHitSetContainer.h>
0012
0013 #include <trackbase/ClusHitsVerbosev1.h>
0014 #include <trackbase/RawHit.h>
0015 #include <trackbase/RawHitSet.h>
0016 #include <trackbase/RawHitSetContainer.h>
0017
0018 #include <g4detectors/PHG4CylinderGeom.h>
0019 #include <g4detectors/PHG4CylinderGeomContainer.h>
0020
0021 #include <fun4all/Fun4AllReturnCodes.h>
0022 #include <fun4all/SubsysReco.h>
0023
0024 #include <phool/PHCompositeNode.h>
0025 #include <phool/PHIODataNode.h>
0026 #include <phool/PHNode.h>
0027 #include <phool/PHNodeIterator.h>
0028 #include <phool/PHObject.h>  // for PHObject
0029 #include <phool/getClass.h>
0030 #include <phool/phool.h>
0031
0032 #include <boost/graph/adjacency_list.hpp>
0033 #include <boost/graph/connected_components.hpp>
0034
0035 #include <array>
0036 #include <cmath>
0037 #include <iostream>
0038 #include <limits>
0039 #include <memory>  // for unique_ptr, make_...
0040 #include <set>
0041 #include <vector>  // for vector
0042
0043 namespace
0044 {
0045
0046   /// convenience square method
0047   template <class T>
0048   inline constexpr T square(const T& x)
0049   {
0050     return x * x;
0051   }
0052 }  // namespace
0053
0054 bool InttClusterizer::ladder_are_adjacent(const std::pair<TrkrDefs::hitkey, TrkrHit*>& lhs, const std::pair<TrkrDefs::hitkey, TrkrHit*>& rhs, const int layer)
0055 {
0056   if (get_z_clustering(layer))
0057   {
0058     if (std::abs(InttDefs::getCol(lhs.first) - InttDefs::getCol(rhs.first)) <= 1)
0059     {
0060       if (std::abs(InttDefs::getRow(lhs.first) - InttDefs::getRow(rhs.first)) <= 1)
0061       {
0062         return true;
0063       }
0064     }
0065   }
0066   else if (std::abs(InttDefs::getCol(lhs.first) - InttDefs::getCol(rhs.first)) == 0)
0067   {
0068     if (std::abs(InttDefs::getRow(lhs.first) - InttDefs::getRow(rhs.first)) <= 1)
0069     {
0070       return true;
0071     }
0072   }
0073
0074   return false;
0075 }
0076
0077 bool InttClusterizer::ladder_are_adjacent(RawHit* lhs, RawHit* rhs, const int layer)
0078 {
0079   if (get_z_clustering(layer))
0080   {
0081     if (std::abs((int) (lhs->getPhiBin()) - (int) (rhs->getPhiBin())) <= 1)  // col
0082     {
0083       if (std::abs((int) (lhs->getTBin()) - (int) (rhs->getTBin())) <= 1)  // Row
0084       {
0085         return true;
0086       }
0087     }
0088   }
0089   else if (std::abs((int) (lhs->getPhiBin()) - (int)(rhs->getPhiBin())) <= 1)
0090   {
0091     if (std::abs((int) (lhs->getTBin()) - (int) (rhs->getTBin())) == 0)
0092     {
0093       return true;
0094     }
0095   }
0096
0097   return false;
0098 }
0099
0100 InttClusterizer::InttClusterizer(const std::string& name,
0101                                  unsigned int /*min_layer*/,
0102                                  unsigned int /*max_layer*/)
0103   : SubsysReco(name)
0104 {
0105 }
0106
0107 int InttClusterizer::InitRun(PHCompositeNode* topNode)
0108 {
0109   /*
0110   // get node containing the digitized hits
0111   _hits = findNode::getClass<SvtxHitMap>(topNode, "SvtxHitMap");
0112   if (!_hits)
0113   {
0114     std::cout << PHWHERE << "ERROR: Can't find node SvtxHitMap" << std::endl;
0115     return Fun4AllReturnCodes::ABORTRUN;
0116   }
0117   */
0118
0119   //-----------------
0120   // Add Cluster Node
0121   //-----------------
0122
0123   PHNodeIterator iter(topNode);
0124
0125   // Looking for the DST node
0126   PHCompositeNode* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0127   if (!dstNode)
0128   {
0129     std::cout << PHWHERE << "DST Node missing, doing nothing." << std::endl;
0130     return Fun4AllReturnCodes::ABORTRUN;
0131   }
0132   PHNodeIterator iter_dst(dstNode);
0133
0134   // Create the Cluster and association nodes if required
0135   auto trkrclusters = findNode::getClass<TrkrClusterContainer>(dstNode, "TRKR_CLUSTER");
0136   if (!trkrclusters)
0137   {
0138     PHNodeIterator dstiter(dstNode);
0139     PHCompositeNode* DetNode =
0140         dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0141     if (!DetNode)
0142     {
0143       DetNode = new PHCompositeNode("TRKR");
0144       dstNode->addNode(DetNode);
0145     }
0146
0147     trkrclusters = new TrkrClusterContainerv4;
0148     PHIODataNode<PHObject>* TrkrClusterContainerNode =
0149         new PHIODataNode<PHObject>(trkrclusters, "TRKR_CLUSTER", "PHObject");
0150     DetNode->addNode(TrkrClusterContainerNode);
0151   }
0152
0153   auto clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0154   if (!clusterhitassoc)
0155   {
0156     PHNodeIterator dstiter(dstNode);
0157     PHCompositeNode* DetNode =
0158         dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0159     if (!DetNode)
0160     {
0161       DetNode = new PHCompositeNode("TRKR");
0162       dstNode->addNode(DetNode);
0163     }
0164
0165     clusterhitassoc = new TrkrClusterHitAssocv3;
0166     PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(clusterhitassoc, "TRKR_CLUSTERHITASSOC", "PHObject");
0167     DetNode->addNode(newNode);
0168   }
0169
0170   // Add the multimap holding the INTT cluster-crossing associations
0171   {
0172     PHNodeIterator dstiter(dstNode);
0173     PHCompositeNode* DetNode =
0174         dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0175     if (!DetNode)
0176     {
0177       DetNode = new PHCompositeNode("TRKR");
0178       dstNode->addNode(DetNode);
0179     }
0180
0181     auto clustercrossingassoc = new TrkrClusterCrossingAssocv1;
0182     PHIODataNode<PHObject>* newNode = new PHIODataNode<PHObject>(clustercrossingassoc, "TRKR_CLUSTERCROSSINGASSOC", "PHObject");
0183     DetNode->addNode(newNode);
0184   }
0185
0186   //---------------------
0187   // Calculate Thresholds
0188   //---------------------
0189
0190   CalculateLadderThresholds(topNode);
0191
0192   //----------------
0193   // Report Settings
0194   //----------------
0195
0196   if (Verbosity() > 0)
0197   {
0198     std::cout << "====================== InttClusterizer::InitRun() =====================" << std::endl;
0199     std::cout << " Fraction of expected thickness MIP energy = " << _fraction_of_mip << std::endl;
0200     for (auto& iter_tmp : _thresholds_by_layer)
0201     {
0202       std::cout << " Cluster Threshold in Layer #" << iter_tmp.first << " = " << 1.0e6 * iter_tmp.second << " keV" << std::endl;
0203     }
0204     for (auto& iter_tmp : _make_z_clustering)
0205     {
0206       std::cout << " Z-dimension Clustering in Layer #" << iter_tmp.first << " = " << std::boolalpha << iter_tmp.second << std::noboolalpha << std::endl;
0207     }
0208     for (auto& _make_e_weight : _make_e_weights)
0209     {
0210       std::cout << " Energy weighting clusters in Layer #" << _make_e_weight.first << " = " << std::boolalpha << _make_e_weight.second << std::noboolalpha << std::endl;
0211     }
0212     std::cout << "===========================================================================" << std::endl;
0213   }
0214
0215   if (record_ClusHitsVerbose)
0216   {
0217     // get the node
0218     mClusHitsVerbose = findNode::getClass<ClusHitsVerbosev1>(topNode, "Trkr_SvtxClusHitsVerbose");
0219     if (!mClusHitsVerbose)
0220     {
0221       PHNodeIterator dstiter(dstNode);
0222       auto DetNode = dynamic_cast<PHCompositeNode*>(dstiter.findFirst("PHCompositeNode", "TRKR"));
0223       if (!DetNode)
0224       {
0225         DetNode = new PHCompositeNode("TRKR");
0226         dstNode->addNode(DetNode);
0227       }
0228       mClusHitsVerbose = new ClusHitsVerbosev1();
0229       auto newNode = new PHIODataNode<PHObject>(mClusHitsVerbose, "Trkr_SvtxClusHitsVerbose", "PHObject");
0230       DetNode->addNode(newNode);
0231     }
0232   }
0233
0234   return Fun4AllReturnCodes::EVENT_OK;
0235 }
0236
0237 int InttClusterizer::process_event(PHCompositeNode* topNode)
0238 {
0239   // get node containing the digitized hits
0240   if (!do_read_raw)
0241   {
0242     m_hits = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0243     if (!m_hits)
0244     {
0245       std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
0246       return Fun4AllReturnCodes::ABORTRUN;
0247     }
0248   }
0249   else
0250   {
0251     // get node containing the digitized hits
0252     m_rawhits = findNode::getClass<RawHitSetContainer>(topNode, "TRKR_RAWHITSET");
0253     if (!m_rawhits)
0254     {
0255       std::cout << PHWHERE << "ERROR: Can't find node TRKR_HITSET" << std::endl;
0256       return Fun4AllReturnCodes::ABORTRUN;
0257     }
0258   }
0259
0260   // get node for clusters
0261   m_clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0262   if (!m_clusterlist)
0263   {
0264     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTER." << std::endl;
0265     return Fun4AllReturnCodes::ABORTRUN;
0266   }
0267
0268   // get node for cluster hit associations
0269   m_clusterhitassoc = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0270   if (!m_clusterhitassoc)
0271   {
0272     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERHITASSOC" << std::endl;
0273     return Fun4AllReturnCodes::ABORTRUN;
0274   }
0275
0276   // get node for cluster-crossing associations
0277   m_clustercrossingassoc = findNode::getClass<TrkrClusterCrossingAssoc>(topNode, "TRKR_CLUSTERCROSSINGASSOC");
0278   if (!m_clustercrossingassoc)
0279   {
0280     std::cout << PHWHERE << " ERROR: Can't find TRKR_CLUSTERCROSINGASSOC" << std::endl;
0281     return Fun4AllReturnCodes::ABORTRUN;
0282   }
0283
0284   if (!do_read_raw)
0285   {
0286     ClusterLadderCells(topNode);
0287   }
0288   else
0289   {
0290     ClusterLadderCellsRaw(topNode);
0291   }
0292   PrintClusters(topNode);
0293
0294   return Fun4AllReturnCodes::EVENT_OK;
0295 }
0296
0297 void InttClusterizer::CalculateLadderThresholds(PHCompositeNode* topNode)
0298 {
0299   /*
0300   PHG4CellContainer* cells = findNode::getClass<PHG4CellContainer>(topNode, "G4CELL_INTT");
0301   if (!cells) return;
0302   */
0303
0304   PHG4CylinderGeomContainer* geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
0305   if (!geom_container)
0306   {
0307     return;
0308   }
0309
0310   PHG4CylinderGeomContainer::ConstRange layerrange = geom_container->get_begin_end();
0311   for (PHG4CylinderGeomContainer::ConstIterator layeriter = layerrange.first;
0312        layeriter != layerrange.second;
0313        ++layeriter)
0314   {
0315     // index mapping
0316     int layer = layeriter->second->get_layer();
0317
0318     // cluster threshold
0319     float thickness = (layeriter->second)->get_thickness();
0320     float threshold = _fraction_of_mip * 0.003876 * thickness;
0321     _thresholds_by_layer.insert(std::make_pair(layer, threshold));
0322
0323     // fill in a default z_clustering value if not present
0324     if (_make_z_clustering.find(layer) == _make_z_clustering.end())
0325     {
0326       _make_z_clustering.insert(std::make_pair(layer, false));
0327     }
0328
0329     if (_make_e_weights.find(layer) == _make_e_weights.end())
0330     {
0331       _make_e_weights.insert(std::make_pair(layer, false));
0332     }
0333   }
0334
0335   return;
0336 }
0337
0338 void InttClusterizer::ClusterLadderCells(PHCompositeNode* topNode)
0339 {
0340   if (Verbosity() > 0)
0341   {
0342     std::cout << "Entering InttClusterizer::ClusterLadderCells " << std::endl;
0343   }
0344
0345   //----------
0346   // Get Nodes
0347   //----------
0348
0349   // get the geometry node
0350   PHG4CylinderGeomContainer* geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
0351   if (!geom_container)
0352   {
0353     return;
0354   }
0355
0356   //-----------
0357   // Clustering
0358   //-----------
0359
0360   // loop over the InttHitSet objects
0361   TrkrHitSetContainer::ConstRange hitsetrange =
0362       m_hits->getHitSets(TrkrDefs::TrkrId::inttId);
0363   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0364        hitsetitr != hitsetrange.second;
0365        ++hitsetitr)
0366   {
0367     // Each hitset contains only hits that are clusterizable - i.e. belong to a single sensor
0368     TrkrHitSet* hitset = hitsetitr->second;
0369
0370     if (Verbosity() > 1)
0371     {
0372       std::cout << "InttClusterizer found hitsetkey " << hitsetitr->first << std::endl;
0373     }
0374     if (Verbosity() > 2)
0375     {
0376       hitset->identify();
0377     }
0378
0379     // we have a single hitset, get the info that identifies the sensor
0380     int layer = TrkrDefs::getLayer(hitsetitr->first);
0381     int ladder_z_index = InttDefs::getLadderZId(hitsetitr->first);
0382     int type = (ladder_z_index == 0 || ladder_z_index == 2) ? 0 : 1; // ladder ID 0 and 2 are type-A (1.6 cm), ladder ID 1 and 3 are type-B (2.0 cm)
0383
0384     // we will need the geometry object for this layer to get the global position
0385     CylinderGeomIntt* geom = dynamic_cast<CylinderGeomIntt*>(geom_container->GetLayerGeom(layer));
0386     float pitch = geom->get_strip_y_spacing();
0387     float length = geom->get_strip_z_spacing(type);
0388
0389     // fill a vector of hits to make things easier - gets every hit in the hitset
0390     std::vector<std::pair<TrkrDefs::hitkey, TrkrHit*>> hitvec;
0391     TrkrHitSet::ConstRange hitrangei = hitset->getHits();
0392     for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0393          hitr != hitrangei.second;
0394          ++hitr)
0395     {
0396       hitvec.emplace_back(hitr->first, hitr->second);
0397     }
0398     if (Verbosity() > 2)
0399     {
0400       std::cout << "hitvec.size(): " << hitvec.size() << std::endl;
0401     }
0402
0403     using Graph = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS>;
0404     Graph G;
0405
0406     // Find adjacent strips
0407     for (unsigned int i = 0; i < hitvec.size(); i++)
0408     {
0409       for (unsigned int j = i + 1; j < hitvec.size(); j++)
0410       {
0411         if (ladder_are_adjacent(hitvec[i], hitvec[j], layer))
0412         {
0413           add_edge(i, j, G);
0414         }
0415       }
0416
0417       add_edge(i, i, G);
0418     }
0419
0420     // Find the connections between the vertices of the graph (vertices are the rawhits,
0421     // connections are made when they are adjacent to one another)
0422     std::vector<int> component(num_vertices(G));
0423
0424     // this is the actual clustering, performed by boost
0425     connected_components(G, &component[0]);
0426
0427     // Loop over the components(hit cells) compiling a list of the
0428     // unique connected groups (ie. clusters).
0429     std::set<int> cluster_ids;  // unique components
0430
0431     std::multimap<int, std::pair<TrkrDefs::hitkey, TrkrHit*>> clusters;
0432     for (unsigned int i = 0; i < component.size(); i++)
0433     {
0434       cluster_ids.insert(component[i]);                          // one entry per unique cluster id
0435       clusters.insert(std::make_pair(component[i], hitvec[i]));  // multiple entries per unique cluster id
0436     }
0437
0438     // loop over the cluster ID's and make the clusters from the connected hits
0439     for (int clusid : cluster_ids)
0440     {
0441       // std::cout << " intt clustering: add cluster number " << clusid << std::endl;
0442
0443       // make the cluster directly in the node tree
0444       TrkrDefs::cluskey ckey = TrkrDefs::genClusKey(hitset->getHitSetKey(), clusid);
0445
0446       if (Verbosity() > 2)
0447       {
0448         std::cout << "Filling cluster with key " << ckey << std::endl;
0449       }
0450
0451       // get the bunch crossing number from the hitsetkey
0452       short int crossing = InttDefs::getTimeBucketId(hitset->getHitSetKey());
0453
0454       // Add clusterkey/bunch crossing to mmap
0455       m_clustercrossingassoc->addAssoc(ckey, crossing);
0456
0457       // determine the size of the cluster in phi and z, useful for track fitting the cluster
0458       std::set<int> phibins;
0459       std::set<int> zbins;
0460
0461       // determine the cluster position...
0462       double xlocalsum = 0.0;
0463       double ylocalsum = 0.0;
0464       double zlocalsum = 0.0;
0465       unsigned int clus_adc = 0.0;
0466       unsigned int clus_maxadc = 0.0;
0467       unsigned nhits = 0;
0468       // std::cout << PHWHERE << " ckey " << ckey << ":" << std::endl;
0469
0470       // get all hits for this cluster ID only
0471       const auto clusrange = clusters.equal_range(clusid);
0472       for (auto mapiter = clusrange.first; mapiter != clusrange.second; ++mapiter)
0473       {
0474         // mapiter->second.first  is the hit key
0475         // std::cout << " adding hitkey " << mapiter->second.first << std::endl;
0476         int col = InttDefs::getCol((mapiter->second).first);
0477         int row = InttDefs::getRow((mapiter->second).first);
0478         zbins.insert(col);
0479         phibins.insert(row);
0480
0481         // mapiter->second.second is the hit
0482         unsigned int hit_adc = (mapiter->second).second->getAdc();
0483
0484         // now get the positions from the geometry
0485         double local_hit_location[3] = {0., 0., 0.};
0486         geom->find_strip_center_localcoords(ladder_z_index,
0487                                             row, col,
0488                                             local_hit_location);
0489
0490         if (_make_e_weights[layer])
0491         {
0492           xlocalsum += local_hit_location[0] * (double) hit_adc;
0493           ylocalsum += local_hit_location[1] * (double) hit_adc;
0494           zlocalsum += local_hit_location[2] * (double) hit_adc;
0495         }
0496         else
0497         {
0498           xlocalsum += local_hit_location[0];
0499           ylocalsum += local_hit_location[1];
0500           zlocalsum += local_hit_location[2];
0501         }
0502         if (hit_adc > clus_maxadc)
0503         {
0504           clus_maxadc = hit_adc;
0505         }
0506         clus_adc += hit_adc;
0507         ++nhits;
0508
0509         // add this cluster-hit association to the association map of (clusterkey,hitkey)
0510         m_clusterhitassoc->addAssoc(ckey, mapiter->second.first);
0511
0512         if (Verbosity() > 2)
0513         {
0514           std::cout << "     nhits = " << nhits << std::endl;
0515         }
0516         if (Verbosity() > 2)
0517         {
0518           std::cout << "  From  geometry object: hit x " << local_hit_location[0] << " hit y " << local_hit_location[1] << " hit z " << local_hit_location[2] << std::endl;
0519           std::cout << "     nhits " << nhits << " clusx  = " << xlocalsum / nhits << " clusy " << ylocalsum / nhits << " clusz " << zlocalsum / nhits << " hit_adc " << hit_adc << std::endl;
0520         }
0521       }
0522
0523       static const float invsqrt12 = 1. / sqrt(12);
0524
0525       // scale factors (phi direction)
0526       /*
0527         they corresponds to clusters of size 1 and 2 in phi
0528         other clusters, which are very few and pathological, get a scale factor of 1
0529         These scale factors are applied to produce cluster pulls with width unity
0530       */
0531
0532       float phierror = pitch * invsqrt12;
0533
0534       static constexpr std::array<double, 3> scalefactors_phi = {{0.85, 0.4, 0.33}};
0535       if (phibins.size() == 1 && layer < 5)
0536       {
0537         phierror *= scalefactors_phi[0];
0538       }
0539       else if (phibins.size() == 2 && layer < 5)
0540       {
0541         phierror *= scalefactors_phi[1];
0542       }
0543       else if (phibins.size() == 2 && layer > 4)
0544       {
0545         phierror *= scalefactors_phi[2];
0546       }
0547       // z error.
0548       const float zerror = zbins.size() * length * invsqrt12;
0549
0550       double cluslocaly = std::numeric_limits<double>::quiet_NaN();
0551       double cluslocalz = std::numeric_limits<double>::quiet_NaN();
0552
0553       if (_make_e_weights[layer])
0554       {
0555         cluslocaly = ylocalsum / (double) clus_adc;
0556         cluslocalz = zlocalsum / (double) clus_adc;
0557       }
0558       else
0559       {
0560         cluslocaly = ylocalsum / nhits;
0561         cluslocalz = zlocalsum / nhits;
0562       }
0563
0564       auto clus = std::make_unique<TrkrClusterv5>();
0565       clus->setAdc(clus_adc);
0566       clus->setMaxAdc(clus_maxadc);
0567       clus->setLocalX(cluslocaly);
0568       clus->setLocalY(cluslocalz);
0569       clus->setPhiError(phierror);
0570       clus->setZError(zerror);
0571       clus->setPhiSize(phibins.size());
0572       clus->setZSize(zbins.size());
0573       // All silicon surfaces have a 1-1 map to hitsetkey.
0574       // So set subsurface key to 0
0575       clus->setSubSurfKey(0);
0576
0577       if (Verbosity() > 2)
0578       {
0579         clus->identify();
0580       }
0581
0582       m_clusterlist->addClusterSpecifyKey(ckey, clus.release());
0583
0584     }  // end loop over cluster ID's
0585   }    // end loop over hitsets
0586
0587   if (Verbosity() > 2)
0588   {
0589     // check that the associations were written correctly
0590     std::cout << "After InttClusterizer, cluster-hit associations are:" << std::endl;
0591     m_clusterhitassoc->identify();
0592   }
0593
0594   if (Verbosity() > 0)
0595   {
0596     std::cout << " Cluster-crossing associations are:" << std::endl;
0597     m_clustercrossingassoc->identify();
0598   }
0599
0600   return;
0601 }
0602 void InttClusterizer::ClusterLadderCellsRaw(PHCompositeNode* topNode)
0603 {
0604   if (Verbosity() > 0)
0605   {
0606     std::cout << "Entering InttClusterizer::ClusterLadderCells " << std::endl;
0607   }
0608
0609   //----------
0610   // Get Nodes
0611   //----------
0612
0613   // get the geometry node
0614   PHG4CylinderGeomContainer* geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
0615   if (!geom_container)
0616   {
0617     return;
0618   }
0619
0620   //-----------
0621   // Clustering
0622   //-----------
0623
0624   // loop over the InttHitSet objects
0625   RawHitSetContainer::ConstRange hitsetrange =
0626       m_rawhits->getHitSets(TrkrDefs::TrkrId::inttId);
0627   for (RawHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0628        hitsetitr != hitsetrange.second;
0629        ++hitsetitr)
0630   {
0631     // Each hitset contains only hits that are clusterizable - i.e. belong to a single sensor
0632     RawHitSet* hitset = hitsetitr->second;
0633
0634     if (Verbosity() > 1)
0635     {
0636       std::cout << "InttClusterizer found hitsetkey " << hitsetitr->first << std::endl;
0637     }
0638     if (Verbosity() > 2)
0639     {
0640       hitset->identify();
0641     }
0642
0643     // we have a single hitset, get the info that identifies the sensor
0644     int layer = TrkrDefs::getLayer(hitsetitr->first);
0645     int ladder_z_index = InttDefs::getLadderZId(hitsetitr->first);
0646     int type = (ladder_z_index == 0 || ladder_z_index == 2) ? 0 : 1; // ladder ID 0 and 2 are type-A (1.6 cm), ladder ID 1 and 3 are type-B (2.0 cm)
0647
0648     // we will need the geometry object for this layer to get the global position
0649     CylinderGeomIntt* geom = dynamic_cast<CylinderGeomIntt*>(geom_container->GetLayerGeom(layer));
0650     float pitch = geom->get_strip_y_spacing();
0651     float length = geom->get_strip_z_spacing(type);
0652
0653     // fill a vector of hits to make things easier - gets every hit in the hitset
0654     std::vector<RawHit*> hitvec;
0655     // int sector = InttDefs::getLadderPhiId(hitsetitr->first);
0656     // int side = InttDefs::getLadderZId(hitsetitr->first);
0657
0658     RawHitSet::ConstRange hitrangei = hitset->getHits();
0659     for (RawHitSet::ConstIterator hitr = hitrangei.first;
0660          hitr != hitrangei.second;
0661          ++hitr)
0662     {
0663       // unsigned short iphi = (*hitr)->getPhiBin();
0664       // unsigned short it = (*hitr)->getTBin();
0665       //    std::cout << " intt layer " << layer << " sector: " << sector << " side " << side << " col: " << iphi << " row " << it << std::endl;
0666       hitvec.push_back((*hitr));
0667     }
0668     if (Verbosity() > 2)
0669     {
0670       std::cout << "hitvec.size(): " << hitvec.size() << std::endl;
0671     }
0672
0673     using Graph = boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS>;
0674     Graph G;
0675
0676     // Find adjacent strips
0677     for (unsigned int i = 0; i < hitvec.size(); i++)
0678     {
0679       for (unsigned int j = i + 1; j < hitvec.size(); j++)
0680       {
0681         if (ladder_are_adjacent(hitvec[i], hitvec[j], layer))
0682         {
0683           add_edge(i, j, G);
0684         }
0685       }
0686
0687       add_edge(i, i, G);
0688     }
0689
0690     // Find the connections between the vertices of the graph (vertices are the rawhits,
0691     // connections are made when they are adjacent to one another)
0692     std::vector<int> component(num_vertices(G));
0693
0694     // this is the actual clustering, performed by boost
0695     connected_components(G, &component[0]);
0696
0697     // Loop over the components(hit cells) compiling a list of the
0698     // unique connected groups (ie. clusters).
0699     std::set<int> cluster_ids;  // unique components
0700
0701     std::multimap<int, RawHit*> clusters;
0702     for (unsigned int i = 0; i < component.size(); i++)
0703     {
0704       cluster_ids.insert(component[i]);                          // one entry per unique cluster id
0705       clusters.insert(std::make_pair(component[i], hitvec[i]));  // multiple entries per unique cluster id
0706     }
0707
0708     // loop over the cluster ID's and make the clusters from the connected hits
0709     for (int clusid : cluster_ids)
0710     {
0711       // std::cout << " intt clustering: add cluster number " << clusid << std::endl;
0712       // make the cluster directly in the node tree
0713       TrkrDefs::cluskey ckey = TrkrDefs::genClusKey(hitset->getHitSetKey(), clusid);
0714
0715       if (Verbosity() > 2)
0716       {
0717         std::cout << "Filling cluster with key " << ckey << std::endl;
0718       }
0719
0720       // get the bunch crossing number from the hitsetkey
0721       short int crossing = InttDefs::getTimeBucketId(hitset->getHitSetKey());
0722
0723       // Add clusterkey/bunch crossing to mmap
0724       m_clustercrossingassoc->addAssoc(ckey, crossing);
0725
0726       // determine the size of the cluster in phi and z, useful for track fitting the cluster
0727       std::set<int> phibins;
0728       std::set<int> zbins;
0729
0730       // determine the cluster position...
0731       double xlocalsum = 0.0;
0732       double ylocalsum = 0.0;
0733       double zlocalsum = 0.0;
0734       unsigned int clus_adc = 0.0;
0735       unsigned nhits = 0;
0736
0737       // std::cout << PHWHERE << " ckey " << ckey << ":" << std::endl;
0738       std::map<int, unsigned int> m_phi, m_z;  // hold data for
0739
0740       // get all hits for this cluster ID only
0741       const auto clusrange = clusters.equal_range(clusid);
0742       for (auto mapiter = clusrange.first; mapiter != clusrange.second; ++mapiter)
0743       {
0744         // mapiter->second.first  is the hit key
0745         // std::cout << " adding hitkey " << mapiter->second.first << std::endl;
0746         const auto energy = (mapiter->second)->getAdc();
0747         int col = (mapiter->second)->getPhiBin();
0748         int row = (mapiter->second)->getTBin();
0749         //      std::cout << " found Tbin(row) " << row << " Phibin(col) " << col << std::endl;
0750         zbins.insert(col);
0751         phibins.insert(row);
0752
0753         if (mClusHitsVerbose)
0754         {
0755           auto pnew = m_phi.try_emplace(row, energy);
0756           if (!pnew.second)
0757           {
0758             pnew.first->second += energy;
0759           }
0760
0761           pnew = m_z.try_emplace(col, energy);
0762           if (!pnew.second)
0763           {
0764             pnew.first->second += energy;
0765           }
0766         }
0767
0768         // mapiter->second.second is the hit
0769         unsigned int hit_adc = (mapiter->second)->getAdc();
0770
0771         // now get the positions from the geometry
0772         double local_hit_location[3] = {0., 0., 0.};
0773         geom->find_strip_center_localcoords(ladder_z_index,
0774                                             row, col,
0775                                             local_hit_location);
0776
0777         if (_make_e_weights[layer])
0778         {
0779           xlocalsum += local_hit_location[0] * (double) hit_adc;
0780           ylocalsum += local_hit_location[1] * (double) hit_adc;
0781           zlocalsum += local_hit_location[2] * (double) hit_adc;
0782         }
0783         else
0784         {
0785           xlocalsum += local_hit_location[0];
0786           ylocalsum += local_hit_location[1];
0787           zlocalsum += local_hit_location[2];
0788         }
0789
0790         clus_adc += hit_adc;
0791         ++nhits;
0792
0793         // add this cluster-hit association to the association map of (clusterkey,hitkey)
0794         //      m_clusterhitassoc->addAssoc(ckey, mapiter->second.first);
0795
0796         if (Verbosity() > 2)
0797         {
0798           std::cout << "     nhits = " << nhits << std::endl;
0799         }
0800         if (Verbosity() > 2)
0801         {
0802           std::cout << "  From  geometry object: hit x " << local_hit_location[0] << " hit y " << local_hit_location[1] << " hit z " << local_hit_location[2] << std::endl;
0803           std::cout << "     nhits " << nhits << " clusx  = " << xlocalsum / nhits << " clusy " << ylocalsum / nhits << " clusz " << zlocalsum / nhits << " hit_adc " << hit_adc << std::endl;
0804         }
0805       }
0806
0807       if (mClusHitsVerbose)
0808       {
0809         if (Verbosity() > 10)
0810         {
0811           for (auto const& hit : m_phi)
0812           {
0813             std::cout << " m_phi(" << hit.first << " : " << hit.second << ") " << std::endl;
0814           }
0815         }
0816         for (auto& hit : m_phi)
0817         {
0818           mClusHitsVerbose->addPhiHit(hit.first, (float) hit.second);
0819         }
0820         for (auto& hit : m_z)
0821         {
0822           mClusHitsVerbose->addZHit(hit.first, (float) hit.second);
0823         }
0824         mClusHitsVerbose->push_hits(ckey);
0825       }
0826
0827       static const float invsqrt12 = 1. / sqrt(12);
0828
0829       // scale factors (phi direction)
0830       /*
0831         they corresponds to clusters of size 1 and 2 in phi
0832         other clusters, which are very few and pathological, get a scale factor of 1
0833         These scale factors are applied to produce cluster pulls with width unity
0834       */
0835
0836       float phierror = pitch * invsqrt12;
0837
0838       static constexpr std::array<double, 3> scalefactors_phi = {{0.85, 0.4, 0.33}};
0839       if (phibins.size() == 1 && layer < 5)
0840       {
0841         phierror *= scalefactors_phi[0];
0842       }
0843       else if (phibins.size() == 2 && layer < 5)
0844       {
0845         phierror *= scalefactors_phi[1];
0846       }
0847       else if (phibins.size() == 2 && layer > 4)
0848       {
0849         phierror *= scalefactors_phi[2];
0850       }
0851       // z error.
0852       const float zerror = zbins.size() * length * invsqrt12;
0853
0854       double cluslocaly = std::numeric_limits<double>::quiet_NaN();
0855       double cluslocalz = std::numeric_limits<double>::quiet_NaN();
0856
0857       if (_make_e_weights[layer])
0858       {
0859         cluslocaly = ylocalsum / (double) clus_adc;
0860         cluslocalz = zlocalsum / (double) clus_adc;
0861       }
0862       else
0863       {
0864         cluslocaly = ylocalsum / nhits;
0865         cluslocalz = zlocalsum / nhits;
0866       }
0867
0868       auto clus = std::make_unique<TrkrClusterv5>();
0869       clus->setAdc(clus_adc);
0870       clus->setLocalX(cluslocaly);
0871       clus->setLocalY(cluslocalz);
0872       clus->setPhiError(phierror);
0873       clus->setZError(zerror);
0874       clus->setPhiSize(phibins.size());
0875       clus->setZSize(zbins.size());
0876       // All silicon surfaces have a 1-1 map to hitsetkey.
0877       // So set subsurface key to 0
0878       clus->setSubSurfKey(0);
0879
0880       if (Verbosity() > 2)
0881       {
0882         clus->identify();
0883       }
0884
0885       m_clusterlist->addClusterSpecifyKey(ckey, clus.release());
0886
0887     }  // end loop over cluster ID's
0888   }    // end loop over hitsets
0889
0890   if (Verbosity() > 2)
0891   {
0892     // check that the associations were written correctly
0893     std::cout << "After InttClusterizer, cluster-hit associations are:" << std::endl;
0894     m_clusterhitassoc->identify();
0895   }
0896
0897   if (Verbosity() > 0)
0898   {
0899     std::cout << " Cluster-crossing associations are:" << std::endl;
0900     m_clustercrossingassoc->identify();
0901   }
0902
0903   return;
0904 }
0905
0906 void InttClusterizer::PrintClusters(PHCompositeNode* topNode)
0907 {
0908   if (Verbosity() > 1)
0909   {
0910     TrkrClusterContainer* clusterlist = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0911     if (!clusterlist)
0912     {
0913       return;
0914     }
0915
0916     std::cout << "================= After InttClusterizer::process_event() ====================" << std::endl;
0917
0918     std::cout << " There are " << clusterlist->size() << " clusters recorded: " << std::endl;
0919
0920     clusterlist->identify();
0921
0922     std::cout << "===========================================================================" << std::endl;
0923   }
0924
0925   return;
0926 }