Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-17 09:22:13 

0001 #include "TpcClusterBuilder.h"
0002 
0003 #include <trackbase/ClusHitsVerbosev1.h>
0004 #include <trackbase/TpcDefs.h>
0005 #include <trackbase/TrkrClusterContainer.h>
0006 #include <trackbase/TrkrClusterv4.h>
0007 #include <trackbase/TrkrHitSet.h>
0008 #include <trackbase/TrkrHitSetContainerv1.h>
0009 #include <trackbase/TrkrHitv2.h>  // for TrkrHit
0010 
0011 #include <g4detectors/PHG4TpcGeom.h>
0012 #include <g4detectors/PHG4TpcGeomContainer.h>
0013 
0014 #include <g4tracking/TrkrTruthTrack.h>
0015 #include <g4tracking/TrkrTruthTrackContainer.h>
0016 
0017 
0018 
0019 #include <algorithm>
0020 #include <cmath>  // for sqrt, cos, sin
0021 #include <format>
0022 #include <ios>
0023 #include <iostream>
0024 #include <limits>
0025 #include <map>  // for _Rb_tree_cons...
0026 #include <set>
0027 
0028 namespace
0029 {
0030   template <class T>
0031   constexpr T square(const T& x)
0032   {
0033     return x * x;
0034   }
0035 }  // namespace
0036 
0037 void TpcClusterBuilder::cluster_hits(TrkrTruthTrack* track)
0038 {
0039   // now build the TrkrCluster
0040   TrkrHitSetContainer::ConstRange hitsetrange;
0041   hitsetrange = m_hits->getHitSets(TrkrDefs::TrkrId::tpcId);
0042 
0043   //-----------------------------------------------------
0044   // from TpcClusterizer.cc ~line: 837
0045   //-----------------------------------------------------
0046 
0047   for (TrkrHitSetContainer::ConstIterator hitsetitr = hitsetrange.first;
0048        hitsetitr != hitsetrange.second;
0049        ++hitsetitr)
0050   {
0051     // if (count>0)continue;
0052     TrkrDefs::hitsetkey hitsetkey = hitsetitr->first;
0053     TrkrHitSet* hitset = hitsetitr->second;
0054     unsigned int layer = TrkrDefs::getLayer(hitsetitr->first);
0055     int side = TpcDefs::getSide(hitsetitr->first);
0056     unsigned int sector = TpcDefs::getSectorId(hitsetitr->first);
0057     PHG4TpcGeom* layergeom = geom_container->GetLayerCellGeom(layer);
0058 
0059     // get the maximum and minimum phi and time
0060     unsigned short NPhiBins = (unsigned short) layergeom->get_phibins();
0061     unsigned short NPhiBinsSector = NPhiBins / 12;
0062     unsigned short NTBins = (unsigned short) layergeom->get_zbins();
0063     unsigned short NTBinsSide = NTBins;
0064     unsigned short NTBinsMin = 0;
0065     unsigned short PhiOffset = NPhiBinsSector * sector;
0066     unsigned short TOffset = NTBinsMin;
0067     AdcClockPeriod = layergeom->get_zstep();
0068     double m_tdriftmax = AdcClockPeriod * NTBins / 2.0;
0069 
0070     unsigned short phibins = NPhiBinsSector;
0071     unsigned short phioffset = PhiOffset;
0072     unsigned short tbins = NTBinsSide;
0073     unsigned short toffset = TOffset;
0074 
0075     // loop over the hits in this cluster
0076     double t_sum = 0.0;
0077     double adc_sum = 0.0;
0078 
0079     /* double phi_sum = 0.0; */
0080     double iphi_sum = 0.0;
0081     // double t2_sum = 0.0;
0082     // double phi2_sum = 0.0;
0083 
0084     double radius = layergeom->get_radius();  // returns center of layer
0085 
0086     int phibinhi = -1;
0087     int phibinlo = std::numeric_limits<int>::max();
0088     int tbinhi = -1;
0089     int tbinlo = std::numeric_limits<int>::max();
0090 
0091     auto ihit_list = hitset->getHits();
0092 
0093     const int iphi_max = phioffset + phibins;
0094     const int it_max = toffset + tbins;
0095 
0096     double sum_adc{0.};  // energy = 4 * adc at this point
0097                          //
0098     // accumulate energy from all hits that are not out of range
0099     for (auto iter = ihit_list.first; iter != ihit_list.second; ++iter)
0100     {
0101       int iphi = TpcDefs::getPad(iter->first);  //- phioffset;
0102       int it = TpcDefs::getTBin(iter->first);   // - toffset;
0103       if (iphi < phioffset || iphi > iphi_max || it < toffset || it > it_max)
0104       {
0105         continue;
0106       }
0107       sum_adc += iter->second->getAdc();
0108     }
0109     const double threshold = sum_adc * m_pixel_thresholdrat;
0110 
0111     // FIXME -- see why the hits are so scattered
0112     std::set<int> v_iphi;
0113     std::set<int> v_it;  // FIXME
0114     std::map<int, unsigned int> m_iphi;
0115     std::map<int, unsigned int> m_it;
0116     std::map<int, unsigned int> m_iphiCut;
0117     std::map<int, unsigned int> m_itCut;  // FIXME
0118     for (auto iter = ihit_list.first; iter != ihit_list.second; ++iter)
0119     {
0120       unsigned int adc = iter->second->getAdc();
0121       if (adc <= 0)
0122       {
0123         continue;
0124       }
0125       int iphi = TpcDefs::getPad(iter->first);  //- phioffset;
0126       int it = TpcDefs::getTBin(iter->first);   // - toffset;
0127       if (iphi < phioffset || iphi > iphi_max)
0128       {
0129         std::cout << "WARNING phibin out of range: " << iphi << " | " << phibins << std::endl;
0130         continue;
0131       }
0132       if (it < toffset || it > it_max)
0133       {
0134         std::cout << "WARNING z bin out of range: " << it << " | " << tbins << std::endl;
0135         continue;
0136       }
0137       if (adc < threshold)
0138       {
0139         if (mClusHitsVerbose)
0140         {
0141           auto pnew = m_iphiCut.try_emplace(iphi, adc);
0142           if (!pnew.second)
0143           {
0144             pnew.first->second += adc;
0145           }
0146 
0147           pnew = m_itCut.try_emplace(it, adc);
0148           if (!pnew.second)
0149           {
0150             pnew.first->second += adc;
0151           }
0152         }
0153         continue;
0154       }
0155 
0156       v_iphi.insert(iphi);
0157       v_it.insert(it);
0158 
0159       // fill m_iphi
0160       auto pnew = m_iphi.try_emplace(iphi, adc);
0161       if (!pnew.second)
0162       {
0163         pnew.first->second += adc;
0164       }
0165       pnew = m_it.try_emplace(it, adc);
0166       if (!pnew.second)
0167       {
0168         pnew.first->second += adc;
0169       }
0170 
0171       phibinhi = std::max(iphi, phibinhi);
0172       phibinlo = std::min(iphi, phibinlo);
0173       tbinhi = std::max(it, tbinhi);
0174       tbinlo = std::min(it, tbinlo);
0175 
0176       iphi_sum += iphi * adc;
0177       // phi2_sum += square(phi_center)*adc;
0178 
0179       // update t sums
0180       double t = layergeom->get_zcenter(it);
0181       t_sum += t * adc;
0182       // t2_sum += square(t)*adc;
0183 
0184       adc_sum += adc;
0185     }
0186     if (mClusHitsVerbose)
0187     {
0188       if (verbosity > 10)
0189       {
0190         for (auto& hit : m_iphi)
0191         {
0192           std::cout << " m_phi(" << hit.first << " : " << hit.second << ") ";
0193         }
0194       }
0195       /* std::cout << " MELON 0 m_iphi "; */
0196       /* for (auto& hit : m_iphi)   std::cout  << hit.first << "|" << hit.second << " "; */
0197       /* std::cout << std::endl; */
0198       /* std::cout << " MELON 1 m_it   "; */
0199       /* for (auto& hit : m_it) std::cout    << hit.first << "|" << hit.second << " "; */
0200       /* std::cout << std::endl; */
0201       for (auto& hit : m_iphi)
0202       {
0203         mClusHitsVerbose->addPhiHit(hit.first, hit.second);
0204       }
0205       for (auto& hit : m_it)
0206       {
0207         mClusHitsVerbose->addZHit(hit.first, hit.second);
0208       }
0209       for (auto& hit : m_iphiCut)
0210       {
0211         mClusHitsVerbose->addPhiCutHit(hit.first, hit.second);
0212       }
0213       for (auto& hit : m_itCut)
0214       {
0215         mClusHitsVerbose->addZCutHit(hit.first, hit.second);
0216       }
0217     }
0218 
0219     // This is the global position
0220     double clusiphi = iphi_sum / adc_sum;
0221     double clusphi = layergeom->get_phi(clusiphi, side);
0222     /* double clusphi = phi_sum / adc_sum; */
0223     float clusx = radius * cos(clusphi);
0224     float clusy = radius * sin(clusphi);
0225     double clust = t_sum / adc_sum;
0226     double zdriftlength = clust * m_tGeometry->get_drift_velocity();
0227     // convert z drift length to z position in the TPC
0228     double clusz = m_tdriftmax * m_tGeometry->get_drift_velocity() - zdriftlength;
0229     if (side == 0)
0230     {
0231       clusz = -clusz;
0232     }
0233 
0234     char tsize = tbinhi - tbinlo + 1;
0235     char phisize = phibinhi - phibinlo + 1;
0236 
0237     if (tsize < 0 && verbosity > 1)
0238     {
0239       std::cout << " FIXME z4 tsize: " << ((int) tsize) << " " << tbinlo << " to " << tbinhi << std::endl;
0240     }
0241 
0242     // -------------------------------------------------
0243     // -------------------------------------------------
0244     // debug here: FIXME
0245     //
0246     /* std::cout << " FIXME phisize " << ((int) phisize) << std::endl; */
0247     // FIXME
0248     if (false)
0249     {  // Printing for debugging
0250       if ((int) phisize > 10 || (int) tsize > 8)
0251       {
0252         int _size_phi = ((int) phisize);
0253         int _nbins_phi = v_iphi.size();
0254         int _delta_phi = abs(_size_phi - _nbins_phi);
0255         int _size_z = ((int) tsize);
0256         int _nbins_z = v_it.size();
0257         int _delta_z = abs(_size_z - _nbins_z);
0258 
0259         //        TString fmt;
0260         std::cout << " x|" << _delta_phi << "|" << _delta_z
0261                   << "| new node FIXME A1  layer(" << layer
0262                   << ") (nset:size) phi("
0263                   << _nbins_phi << ":" << _size_phi << ") z("
0264                   << _nbins_z << ":" << _size_z << ") "
0265                   << "trkId(" << track->getTrackid() << ") trkpT(" << track->getPt() << ")"
0266                   << std::endl;
0267         if (phisize > 10)
0268         {
0269           std::cout << "  iphi-from-(";
0270           int _prev = -1;
0271           double tempsum = 0.;
0272           for (auto _ : v_iphi)
0273           {
0274             if (_prev == -1)
0275             {
0276               std::cout << _ << "): ";
0277             }
0278             else
0279             {
0280               int _diff = ((int) _ - _prev - 1);
0281               if (_diff != 0)
0282               {
0283                 std::cout << ">" << _diff << ">";
0284               }
0285             }
0286             /* std::cout << std::setprecision(2) << std::fixed; */
0287             double _rat = (float) m_iphi[_] / (float) adc_sum;
0288             std::cout << std::format("{:.2f}", _rat) << " ";
0289             tempsum += _rat;
0290             _prev = _;
0291           }
0292           if (tempsum < 0.999)
0293           {
0294             std::cout << " Z3 sumphirat: " << tempsum;
0295           }
0296           std::cout << std::endl;
0297         }
0298         if (tsize > 8)
0299         {
0300           int _prev = -1;
0301           double tempsum = 0.;
0302           std::cout << "  iz-from-(";
0303           for (auto _ : v_it)
0304           {
0305             if (_prev == -1)
0306             {
0307               std::cout << _ << "): ";
0308             }
0309             else
0310             {
0311               int _diff = ((int) _ - _prev - 1);
0312               if (_diff != 0)
0313               {
0314                 std::cout << ">" << _diff << ">";
0315               }
0316             }
0317             /* std::cout << std::setprecision(2) << std::fixed; */
0318             double _rat = (float) m_it[_] / (float) adc_sum;
0319             std::cout << std::format("{:.2f}", _rat) << " ";
0320             tempsum += _rat;
0321             _prev = _;
0322           }
0323           if (tempsum < 0.999)
0324           {
0325             std::cout << " Z3 sumzrat: " << tempsum;
0326           }
0327         }
0328         std::cout << std::endl;
0329       }
0330     }  // end debug printing
0331 
0332     // get the global vector3 to then get the surface local phi and z
0333     Acts::Vector3 global(clusx, clusy, clusz);
0334     TrkrDefs::subsurfkey subsurfkey = 0;
0335 
0336     Surface surface = m_tGeometry->get_tpc_surface_from_coords(
0337         hitsetkey, global, subsurfkey);
0338 
0339     if (!surface)
0340     {
0341       if (verbosity)
0342       {
0343         std::cout << "Can't find the surface! with hitsetkey " << ((int) hitsetkey) << std::endl;
0344       }
0345       /* if (mClusHitsVerbose) mClusHitsVerbose->Reset(); */
0346       continue;
0347     }
0348 
0349     // SAMPA shaping bias correction
0350     clust = clust + m_sampa_tbias;
0351 
0352     global *= Acts::UnitConstants::cm;
0353 
0354     Acts::Vector3 local = surface->transform(m_tGeometry->geometry().getGeoContext()).inverse() * global;
0355     local /= Acts::UnitConstants::cm;
0356 
0357     auto* cluster = new TrkrClusterv4;  //
0358     cluster->setAdc(adc_sum);
0359     /* cluster->setOverlap(ntouch); */
0360     /* cluster->setEdge(nedge); */
0361     cluster->setPhiSize(phisize);
0362     cluster->setZSize(tsize);
0363     cluster->setSubSurfKey(subsurfkey);
0364     cluster->setLocalX(local(0));
0365     cluster->setLocalY(clust);
0366 
0367     // add the clusterkey
0368     auto empl = hitsetkey_cnt.try_emplace(hitsetkey, 0);
0369     if (!empl.second)
0370     {
0371       empl.first->second += 1;
0372     }
0373     TrkrDefs::cluskey cluskey = TrkrDefs::genClusKey(hitsetkey, hitsetkey_cnt[hitsetkey]);
0374     m_clusterlist->addClusterSpecifyKey(cluskey, cluster);
0375 
0376     track->addCluster(cluskey);
0377     if (mClusHitsVerbose)
0378     {
0379       mClusHitsVerbose->push_hits(cluskey);
0380       /* std::cout << " FIXME z1 ClusHitsVerbose.size: " << mClusHitsVerbose->getMap().size() << std::endl; */
0381     }
0382   }
0383   m_hits->Reset();
0384 }
0385 
0386 /* void TpcClusterBuilder::set_current_track(TrkrTruthTrack* track) { */
0387 /* current_track = track; */
0388 /* ++ n_tracks; */
0389 /* } */
0390 
0391 void TpcClusterBuilder::addhitset(
0392     TrkrDefs::hitsetkey hitsetkey,
0393     TrkrDefs::hitkey hitkey,
0394     float neffelectrons)
0395 {
0396   // copy of code in PHG4TpcPadPlaneReadout::MapToPadPlane, with a switch
0397   // to ignore non embedded tracks
0398   if (!b_collect_hits)
0399   {
0400     return;
0401   }
0402 
0403   // Add the hitset to the current embedded track
0404   // Code from PHG4TpcPadPlaneReadout::MapToPadPlane (around lines {}.cc::386-401)
0405   TrkrHitSetContainer::Iterator hitsetit = m_hits->findOrAddHitSet(hitsetkey);
0406   // See if this hit already exists
0407   TrkrHit* hit = nullptr;
0408   hit = hitsetit->second->getHit(hitkey);
0409   if (!hit)
0410   {
0411     // create a new one
0412     hit = new TrkrHitv2();
0413     hitsetit->second->addHitSpecificKey(hitkey, hit);
0414   }
0415   // Either way, add the energy to it  -- adc values will be added at digitization
0416   hit->addEnergy(neffelectrons);
0417 }
0418 
0419 void TpcClusterBuilder::clear_hitsetkey_cnt()
0420 {
0421   hitsetkey_cnt.clear();
0422 }
0423 
0424 void TpcClusterBuilder::print(TrkrTruthTrackContainer* truth_tracks, int nclusprint) const
0425 {
0426   std::cout << " ------------- content of TrkrTruthTrackContainer ---------- " << std::endl;
0427   auto& tmap = truth_tracks->getMap();
0428   std::cout << " Number of tracks:  xyz db : " << tmap.size() << std::endl;
0429   for (auto& _pair : tmap)
0430   {
0431     auto& track = _pair.second;
0432     std::cout << std::format("id({:2d}) phi:eta:pt(", static_cast<int>(track->getTrackid()))
0433           << std::format("{:5.2f}:{:5.2f}:{:5.2f}", track->getPhi(), track->getPseudoRapidity(), track->getPt())
0434               << ") nclusters(" << track->getClusters().size() << ") ";
0435     if (verbosity <= 10)
0436     {
0437       std::cout << std::endl;
0438     }
0439     else
0440     {
0441       int nclus = 0;
0442       for (auto cluskey : track->getClusters())
0443       {
0444         std::cout << " "
0445                   << ((int) TrkrDefs::getHitSetKeyFromClusKey(cluskey)) << ":index(" << ((int) TrkrDefs::getClusIndex(cluskey)) << ")";
0446         ++nclus;
0447         if (nclusprint > 0 && nclus >= nclusprint)
0448         {
0449           std::cout << " ... ";
0450           break;
0451         }
0452       }
0453     }
0454   }
0455   std::cout << " ----- end of tracks in TrkrrTruthTrackContainer ------ " << std::endl;
0456 }
0457 
0458 void TpcClusterBuilder::print_file(
0459     TrkrTruthTrackContainer* truth_tracks, const std::string& ofile_name)
0460 {
0461   std::ofstream fout;
0462   fout.open(ofile_name.c_str());
0463   fout << " ------------- content of TrkrTruthTrackContainer ---------- " << std::endl;
0464   auto& tmap = truth_tracks->getMap();
0465   fout << " Number of tracks: " << tmap.size() << std::endl;
0466   for (auto& _pair : tmap)
0467   {
0468     auto& track = _pair.second;
0469     fout << " id( " << track->getTrackid() << ")  phi:eta:pt(" << track->getPhi() << ":" << track->getPseudoRapidity() << ":" << track->getPt() << ") nclusters("
0470          << track->getClusters().size() << ") ";
0471     //    int nclus = 0;
0472     for (auto cluskey : track->getClusters())
0473     {
0474       auto* C = m_clusterlist->findCluster(cluskey);
0475       fout << " "
0476            << ((int) TrkrDefs::getHitSetKeyFromClusKey(cluskey)) << ":" << ((int) TrkrDefs::getClusIndex(cluskey)) << "->adc:X:phisize:Y:zsize("
0477            << C->getAdc() << ":"
0478            << C->getLocalX() << ":"
0479            << C->getPhiSize() << ":"
0480            << C->getLocalY() << ":"
0481            << C->getZSize() << ") ";
0482       //      ++nclus;
0483     }
0484     fout << std::endl;
0485   }
0486   fout << " ----- end of tracks in TrkrrTruthTrackContainer ------ " << std::endl;
0487   fout.close();
0488 }
0489 
0490 void TpcClusterBuilder::set_input_nodes(
0491     TrkrClusterContainer* _truth_cluster_container, ActsGeometry* ActsGeometry, PHG4TpcGeomContainer* _geom_container, ClusHitsVerbosev1* _clushitsverbose)
0492 {
0493   m_clusterlist = _truth_cluster_container;
0494   m_tGeometry = ActsGeometry;
0495   geom_container = _geom_container;
0496   mClusHitsVerbose = _clushitsverbose;
0497   m_needs_input_nodes = false;
0498 };
[ 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!