Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-08-05 08:18:18 

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/PHG4TpcCylinderGeom.h>
0012 
0013 #include <g4tracking/TrkrTruthTrack.h>
0014 #include <g4tracking/TrkrTruthTrackContainer.h>
0015 
0016 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0017 
0018 #include <boost/format.hpp>
0019 
0020 #include <algorithm>
0021 #include <cmath>  // for sqrt, cos, sin
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   inline 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     PHG4TpcCylinderGeom* 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, v_it;                                    // FIXME
0113     std::map<int, unsigned int> m_iphi, m_it, m_iphiCut, m_itCut;  // FIXME
0114     for (auto iter = ihit_list.first; iter != ihit_list.second; ++iter)
0115     {
0116       unsigned int adc = iter->second->getAdc();
0117       if (adc <= 0)
0118       {
0119         continue;
0120       }
0121       int iphi = TpcDefs::getPad(iter->first);  //- phioffset;
0122       int it = TpcDefs::getTBin(iter->first);   // - toffset;
0123       if (iphi < phioffset || iphi > iphi_max)
0124       {
0125         std::cout << "WARNING phibin out of range: " << iphi << " | " << phibins << std::endl;
0126         continue;
0127       }
0128       if (it < toffset || it > it_max)
0129       {
0130         std::cout << "WARNING z bin out of range: " << it << " | " << tbins << std::endl;
0131         continue;
0132       }
0133       if (adc < threshold)
0134       {
0135         if (mClusHitsVerbose)
0136         {
0137           auto pnew = m_iphiCut.try_emplace(iphi, adc);
0138           if (!pnew.second)
0139           {
0140             pnew.first->second += adc;
0141           }
0142 
0143           pnew = m_itCut.try_emplace(it, adc);
0144           if (!pnew.second)
0145           {
0146             pnew.first->second += adc;
0147           }
0148         }
0149         continue;
0150       }
0151 
0152       v_iphi.insert(iphi);
0153       v_it.insert(it);
0154 
0155       // fill m_iphi
0156       auto pnew = m_iphi.try_emplace(iphi, adc);
0157       if (!pnew.second)
0158       {
0159         pnew.first->second += adc;
0160       }
0161       pnew = m_it.try_emplace(it, adc);
0162       if (!pnew.second)
0163       {
0164         pnew.first->second += adc;
0165       }
0166 
0167       if (iphi > phibinhi)
0168       {
0169         phibinhi = iphi;
0170       }
0171       if (iphi < phibinlo)
0172       {
0173         phibinlo = iphi;
0174       }
0175       if (it > tbinhi)
0176       {
0177         tbinhi = it;
0178       }
0179       if (it < tbinlo)
0180       {
0181         tbinlo = it;
0182       }
0183 
0184       iphi_sum += iphi * adc;
0185       // phi2_sum += square(phi_center)*adc;
0186 
0187       // update t sums
0188       double t = layergeom->get_zcenter(it);
0189       t_sum += t * adc;
0190       // t2_sum += square(t)*adc;
0191 
0192       adc_sum += adc;
0193     }
0194     if (mClusHitsVerbose)
0195     {
0196       if (verbosity > 10)
0197       {
0198         for (auto& hit : m_iphi)
0199         {
0200           std::cout << " m_phi(" << hit.first << " : " << hit.second << ") ";
0201         }
0202       }
0203       /* std::cout << " MELON 0 m_iphi "; */
0204       /* for (auto& hit : m_iphi)   std::cout  << hit.first << "|" << hit.second << " "; */
0205       /* std::cout << std::endl; */
0206       /* std::cout << " MELON 1 m_it   "; */
0207       /* for (auto& hit : m_it) std::cout    << hit.first << "|" << hit.second << " "; */
0208       /* std::cout << std::endl; */
0209       for (auto& hit : m_iphi)
0210       {
0211         mClusHitsVerbose->addPhiHit(hit.first, hit.second);
0212       }
0213       for (auto& hit : m_it)
0214       {
0215         mClusHitsVerbose->addZHit(hit.first, hit.second);
0216       }
0217       for (auto& hit : m_iphiCut)
0218       {
0219         mClusHitsVerbose->addPhiCutHit(hit.first, hit.second);
0220       }
0221       for (auto& hit : m_itCut)
0222       {
0223         mClusHitsVerbose->addZCutHit(hit.first, hit.second);
0224       }
0225     }
0226 
0227     // This is the global position
0228     double clusiphi = iphi_sum / adc_sum;
0229     double clusphi = layergeom->get_phi(clusiphi, side);
0230     /* double clusphi = phi_sum / adc_sum; */
0231     float clusx = radius * cos(clusphi);
0232     float clusy = radius * sin(clusphi);
0233     double clust = t_sum / adc_sum;
0234     double zdriftlength = clust * m_tGeometry->get_drift_velocity();
0235     // convert z drift length to z position in the TPC
0236     double clusz = m_tdriftmax * m_tGeometry->get_drift_velocity() - zdriftlength;
0237     if (side == 0)
0238     {
0239       clusz = -clusz;
0240     }
0241 
0242     char tsize = tbinhi - tbinlo + 1;
0243     char phisize = phibinhi - phibinlo + 1;
0244 
0245     if (tsize < 0 && verbosity > 1)
0246     {
0247       std::cout << " FIXME z4 tsize: " << ((int) tsize) << " " << tbinlo << " to " << tbinhi << std::endl;
0248     }
0249 
0250     // -------------------------------------------------
0251     // -------------------------------------------------
0252     // debug here: FIXME
0253     //
0254     /* std::cout << " FIXME phisize " << ((int) phisize) << std::endl; */
0255     // FIXME
0256     if (false)
0257     {  // Printing for debugging
0258       if ((int) phisize > 10 || (int) tsize > 8)
0259       {
0260         int _size_phi = ((int) phisize);
0261         int _nbins_phi = v_iphi.size();
0262         int _delta_phi = abs(_size_phi - _nbins_phi);
0263         int _size_z = ((int) tsize);
0264         int _nbins_z = v_it.size();
0265         int _delta_z = abs(_size_z - _nbins_z);
0266 
0267         //        TString fmt;
0268         std::cout << " x|" << _delta_phi << "|" << _delta_z
0269                   << "| new node FIXME A1  layer(" << layer
0270                   << ") (nset:size) phi("
0271                   << _nbins_phi << ":" << _size_phi << ") z("
0272                   << _nbins_z << ":" << _size_z << ") "
0273                   << "trkId(" << track->getTrackid() << ") trkpT(" << track->getPt() << ")"
0274                   << std::endl;
0275         if (phisize > 10)
0276         {
0277           std::cout << "  iphi-from-(";
0278           int _prev = -1;
0279           double tempsum = 0.;
0280           for (auto _ : v_iphi)
0281           {
0282             if (_prev == -1)
0283             {
0284               std::cout << _ << "): ";
0285             }
0286             else
0287             {
0288               int _diff = ((int) _ - _prev - 1);
0289               if (_diff != 0)
0290               {
0291                 std::cout << ">" << _diff << ">";
0292               }
0293             }
0294             /* std::cout << std::setprecision(2) << std::fixed; */
0295             double _rat = (float) m_iphi[_] / (float) adc_sum;
0296             std::cout << boost::str(boost::format("%.2f") % _rat) << " ";
0297             tempsum += _rat;
0298             _prev = _;
0299           }
0300           if (tempsum < 0.999)
0301           {
0302             std::cout << " Z3 sumphirat: " << tempsum;
0303           }
0304           std::cout << std::endl;
0305         }
0306         if (tsize > 8)
0307         {
0308           int _prev = -1;
0309           double tempsum = 0.;
0310           std::cout << "  iz-from-(";
0311           for (auto _ : v_it)
0312           {
0313             if (_prev == -1)
0314             {
0315               std::cout << _ << "): ";
0316             }
0317             else
0318             {
0319               int _diff = ((int) _ - _prev - 1);
0320               if (_diff != 0)
0321               {
0322                 std::cout << ">" << _diff << ">";
0323               }
0324             }
0325             /* std::cout << std::setprecision(2) << std::fixed; */
0326             double _rat = (float) m_it[_] / (float) adc_sum;
0327             std::cout << boost::str(boost::format("%.2f") % _rat) << " ";
0328             tempsum += _rat;
0329             _prev = _;
0330           }
0331           if (tempsum < 0.999)
0332           {
0333             std::cout << " Z3 sumzrat: " << tempsum;
0334           }
0335         }
0336         std::cout << std::endl;
0337       }
0338     }  // end debug printing
0339 
0340     // get the global vector3 to then get the surface local phi and z
0341     Acts::Vector3 global(clusx, clusy, clusz);
0342     TrkrDefs::subsurfkey subsurfkey = 0;
0343 
0344     Surface surface = m_tGeometry->get_tpc_surface_from_coords(
0345         hitsetkey, global, subsurfkey);
0346 
0347     if (!surface)
0348     {
0349       if (verbosity)
0350       {
0351         std::cout << "Can't find the surface! with hitsetkey " << ((int) hitsetkey) << std::endl;
0352       }
0353       /* if (mClusHitsVerbose) mClusHitsVerbose->Reset(); */
0354       continue;
0355     }
0356 
0357     // SAMPA shaping bias correction
0358     clust = clust + m_sampa_tbias;
0359 
0360     global *= Acts::UnitConstants::cm;
0361 
0362     Acts::Vector3 local = surface->transform(m_tGeometry->geometry().getGeoContext()).inverse() * global;
0363     local /= Acts::UnitConstants::cm;
0364 
0365     auto cluster = new TrkrClusterv4;  //
0366     cluster->setAdc(adc_sum);
0367     /* cluster->setOverlap(ntouch); */
0368     /* cluster->setEdge(nedge); */
0369     cluster->setPhiSize(phisize);
0370     cluster->setZSize(tsize);
0371     cluster->setSubSurfKey(subsurfkey);
0372     cluster->setLocalX(local(0));
0373     cluster->setLocalY(clust);
0374 
0375     // add the clusterkey
0376     auto empl = hitsetkey_cnt.try_emplace(hitsetkey, 0);
0377     if (!empl.second)
0378     {
0379       empl.first->second += 1;
0380     }
0381     TrkrDefs::cluskey cluskey = TrkrDefs::genClusKey(hitsetkey, hitsetkey_cnt[hitsetkey]);
0382     m_clusterlist->addClusterSpecifyKey(cluskey, cluster);
0383 
0384     track->addCluster(cluskey);
0385     if (mClusHitsVerbose)
0386     {
0387       mClusHitsVerbose->push_hits(cluskey);
0388       /* std::cout << " FIXME z1 ClusHitsVerbose.size: " << mClusHitsVerbose->getMap().size() << std::endl; */
0389     }
0390   }
0391   m_hits->Reset();
0392 }
0393 
0394 /* void TpcClusterBuilder::set_current_track(TrkrTruthTrack* track) { */
0395 /* current_track = track; */
0396 /* ++ n_tracks; */
0397 /* } */
0398 
0399 void TpcClusterBuilder::addhitset(
0400     TrkrDefs::hitsetkey hitsetkey,
0401     TrkrDefs::hitkey hitkey,
0402     float neffelectrons)
0403 {
0404   // copy of code in PHG4TpcPadPlaneReadout::MapToPadPlane, with a switch
0405   // to ignore non embedded tracks
0406   if (!b_collect_hits)
0407   {
0408     return;
0409   }
0410 
0411   // Add the hitset to the current embedded track
0412   // Code from PHG4TpcPadPlaneReadout::MapToPadPlane (around lines {}.cc::386-401)
0413   TrkrHitSetContainer::Iterator hitsetit = m_hits->findOrAddHitSet(hitsetkey);
0414   // See if this hit already exists
0415   TrkrHit* hit = nullptr;
0416   hit = hitsetit->second->getHit(hitkey);
0417   if (!hit)
0418   {
0419     // create a new one
0420     hit = new TrkrHitv2();
0421     hitsetit->second->addHitSpecificKey(hitkey, hit);
0422   }
0423   // Either way, add the energy to it  -- adc values will be added at digitization
0424   hit->addEnergy(neffelectrons);
0425 }
0426 
0427 void TpcClusterBuilder::clear_hitsetkey_cnt()
0428 {
0429   hitsetkey_cnt.clear();
0430 }
0431 
0432 void TpcClusterBuilder::print(
0433     TrkrTruthTrackContainer* truth_tracks, int nclusprint)
0434 {
0435   std::cout << " ------------- content of TrkrTruthTrackContainer ---------- " << std::endl;
0436   auto& tmap = truth_tracks->getMap();
0437   std::cout << " Number of tracks:  xyz db : " << tmap.size() << std::endl;
0438   for (auto& _pair : tmap)
0439   {
0440     auto& track = _pair.second;
0441     std::cout << boost::str(boost::format("id(%2i) phi:eta:pt(") % ((int) track->getTrackid()))
0442               << boost::str(boost::format("%5.2f:%5.2f:%5.2f") % track->getPhi() % track->getPseudoRapidity() % track->getPt())
0443               << ") nclusters(" << track->getClusters().size() << ") ";
0444     if (verbosity <= 10)
0445     {
0446       std::cout << std::endl;
0447     }
0448     else
0449     {
0450       int nclus = 0;
0451       for (auto cluskey : track->getClusters())
0452       {
0453         std::cout << " "
0454                   << ((int) TrkrDefs::getHitSetKeyFromClusKey(cluskey)) << ":index(" << ((int) TrkrDefs::getClusIndex(cluskey)) << ")";
0455         ++nclus;
0456         if (nclusprint > 0 && nclus >= nclusprint)
0457         {
0458           std::cout << " ... ";
0459           break;
0460         }
0461       }
0462     }
0463   }
0464   std::cout << " ----- end of tracks in TrkrrTruthTrackContainer ------ " << std::endl;
0465 }
0466 
0467 void TpcClusterBuilder::print_file(
0468     TrkrTruthTrackContainer* truth_tracks, const std::string& ofile_name)
0469 {
0470   std::ofstream fout;
0471   fout.open(ofile_name.c_str());
0472   fout << " ------------- content of TrkrTruthTrackContainer ---------- " << std::endl;
0473   auto& tmap = truth_tracks->getMap();
0474   fout << " Number of tracks: " << tmap.size() << std::endl;
0475   for (auto& _pair : tmap)
0476   {
0477     auto& track = _pair.second;
0478     fout << " id( " << track->getTrackid() << ")  phi:eta:pt(" << track->getPhi() << ":" << track->getPseudoRapidity() << ":" << track->getPt() << ") nclusters("
0479          << track->getClusters().size() << ") ";
0480 //    int nclus = 0;
0481     for (auto cluskey : track->getClusters())
0482     {
0483       auto C = m_clusterlist->findCluster(cluskey);
0484       fout << " "
0485            << ((int) TrkrDefs::getHitSetKeyFromClusKey(cluskey)) << ":" << ((int) TrkrDefs::getClusIndex(cluskey)) << "->adc:X:phisize:Y:zsize("
0486            << C->getAdc() << ":"
0487            << C->getLocalX() << ":"
0488            << C->getPhiSize() << ":"
0489            << C->getLocalY() << ":"
0490            << C->getZSize() << ") ";
0491 //      ++nclus;
0492     }
0493     fout << std::endl;
0494   }
0495   fout << " ----- end of tracks in TrkrrTruthTrackContainer ------ " << std::endl;
0496   fout.close();
0497 }
0498 
0499 void TpcClusterBuilder::set_input_nodes(
0500     TrkrClusterContainer* _truth_cluster_container, ActsGeometry* ActsGeometry, PHG4TpcCylinderGeomContainer* _geom_container, ClusHitsVerbosev1* _clushitsverbose)
0501 {
0502   m_clusterlist = _truth_cluster_container;
0503   m_tGeometry = ActsGeometry;
0504   geom_container = _geom_container;
0505   mClusHitsVerbose = _clushitsverbose;
0506   m_needs_input_nodes = false;
0507 };
[ 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!