Back to home page

sPhenix code displayed by LXR
 
 

    

File indexing completed on 2025-12-16 09:21:43 

0001 /*!
0002  * \file DSTEmulation.cc
0003  * \author Hugo Pereira Da Costa <hugo.pereira-da-costa@cea.fr>
0004  */
0005 
0006 #include "DSTEmulator.h"
0007 
0008 #include "DSTCompressor.h"
0009 #include "TrackEvaluationContainerv1.h"
0010 
0011 #include <g4main/PHG4Hit.h>
0012 #include <g4main/PHG4HitContainer.h>
0013 #include <g4main/PHG4Particle.h>
0014 #include <g4main/PHG4TruthInfoContainer.h>
0015 
0016 #include <micromegas/MicromegasDefs.h>
0017 
0018 #include <trackbase/InttDefs.h>
0019 #include <trackbase/MvtxDefs.h>
0020 #include <trackbase/TpcDefs.h>
0021 #include <trackbase/TrkrClusterContainer.h>
0022 #include <trackbase/TrkrClusterHitAssoc.h>
0023 #include <trackbase/TrkrClusterv2.h>
0024 #include <trackbase/TrkrDefs.h>
0025 #include <trackbase/TrkrHit.h>
0026 #include <trackbase/TrkrHitSet.h>
0027 #include <trackbase/TrkrHitSetContainer.h>
0028 #include <trackbase/TrkrHitTruthAssoc.h>
0029 #include <trackbase_historic/SvtxTrack.h>
0030 #include <trackbase_historic/SvtxTrackMap.h>
0031 
0032 #include <fun4all/Fun4AllReturnCodes.h>
0033 
0034 #include <phool/PHCompositeNode.h>
0035 #include <phool/PHNodeIterator.h>
0036 #include <phool/PHRandomSeed.h>
0037 #include <phool/getClass.h>
0038 #include <phool/recoConsts.h>
0039 
0040 #include <Acts/Definitions/Units.hpp>
0041 #include <Acts/Surfaces/Surface.hpp>
0042 
0043 #include <TFile.h>
0044 #include <TNtuple.h>
0045 
0046 #include <algorithm>
0047 #include <bitset>
0048 #include <cassert>
0049 #include <iostream>
0050 #include <numeric>
0051 
0052 //_____________________________________________________________________
0053 namespace
0054 {
0055 
0056   //! range adaptor to be able to use range-based for loop
0057   template <class T>
0058   class range_adaptor
0059   {
0060    public:
0061     explicit range_adaptor(const T& range)
0062       : m_range(range)
0063     {
0064     }
0065     const typename T::first_type& begin() { return m_range.first; }
0066     const typename T::second_type& end() { return m_range.second; }
0067 
0068    private:
0069     T m_range;
0070   };
0071 
0072   //! square
0073   template <class T>
0074   constexpr T square(T x)
0075   {
0076     return x * x;
0077   }
0078 
0079   //! radius
0080   template <class T>
0081   constexpr T get_r(T x, T y)
0082   {
0083     return std::sqrt(square(x) + square(y));
0084   }
0085 
0086   //! pt
0087   template <class T>
0088   T get_pt(T px, T py)
0089   {
0090     return std::sqrt(square(px) + square(py));
0091   }
0092 
0093   //! p
0094   template <class T>
0095   T get_p(T px, T py, T pz)
0096   {
0097     return std::sqrt(square(px) + square(py) + square(pz));
0098   }
0099 
0100   //! eta
0101   template <class T>
0102   T get_eta(T p, T pz)
0103   {
0104     return std::log((p + pz) / (p - pz)) / 2;
0105   }
0106 
0107   /*  //! radius
0108   float get_r( PHG4Hit* hit, int i )
0109   {  return get_r( hit->get_x(i), hit->get_y(i) ); }
0110   */
0111   /*
0112   //! calculate the average of member function called on all members in collection
0113   template< float (PHG4Hit::*accessor)(int) const>
0114   float interpolate( std::set<PHG4Hit*> hits, float rextrap )
0115   {
0116     // calculate all terms needed for the interpolation
0117     // need to use double everywhere here due to numerical divergences
0118     double sw = 0;
0119     double swr = 0;
0120     double swr2 = 0;
0121     double swx = 0;
0122     double swrx = 0;
0123 
0124     bool valid( false );
0125     for( const auto& hit:hits )
0126     {
0127 
0128       const double x0 = (hit->*accessor)(0);
0129       const double x1 = (hit->*accessor)(1);
0130       if( std::isnan( x0 ) || std::isnan( x1 ) ) continue;
0131 
0132       const double w = hit->get_edep();
0133       if( w < 0 ) continue;
0134 
0135       valid = true;
0136       const double r0 = get_r( hit, 0 );
0137       const double r1 = get_r( hit, 1 );
0138 
0139       sw += w*2;
0140       swr += w*(r0 + r1);
0141       swr2 += w*(square(r0) + square(r1));
0142       swx += w*(x0 + x1);
0143       swrx += w*(r0*x0 + r1*x1);
0144     }
0145 
0146     if( !valid ) return std::numeric_limits<float>::quiet_NaN();
0147 
0148     const auto alpha = (sw*swrx - swr*swx);
0149     const auto beta = (swr2*swx - swr*swrx);
0150     const auto denom = (sw*swr2 - square(swr));
0151 
0152     return ( alpha*rextrap + beta )/denom;
0153   }
0154   */
0155   /*
0156   //! true if a track is a primary
0157   inline int is_primary( PHG4Particle* particle )
0158   { return particle->get_parent_id() == 0; }
0159   */
0160   //! get mask from track clusters
0161   int64_t get_mask(SvtxTrack* track)
0162   {
0163     return std::accumulate(track->begin_cluster_keys(), track->end_cluster_keys(), int64_t(0),
0164                            [](int64_t value, const TrkrDefs::cluskey& key)
0165                            {
0166                              return TrkrDefs::getLayer(key) < 64 ? value | (1ULL << TrkrDefs::getLayer(key)) : value;  // NOLINT(hicpp-signed-bitwise)
0167                            });
0168   }
0169 
0170   //! return number of clusters of a given type
0171   template <int type>
0172   int get_clusters(SvtxTrack* track)
0173   {
0174     return std::count_if(track->begin_cluster_keys(), track->end_cluster_keys(),
0175                          [](const TrkrDefs::cluskey& key)
0176                          { return TrkrDefs::getTrkrId(key) == type; });
0177   }
0178 
0179   //! create track struct from struct from svx track
0180   TrackEvaluationContainerv1::TrackStruct create_track(SvtxTrack* track)
0181   {
0182     TrackEvaluationContainerv1::TrackStruct trackStruct;
0183 
0184     trackStruct.charge = track->get_charge();
0185     trackStruct.nclusters = track->size_cluster_keys();
0186     trackStruct.mask = get_mask(track);
0187     trackStruct.nclusters_mvtx = get_clusters<TrkrDefs::mvtxId>(track);
0188     trackStruct.nclusters_intt = get_clusters<TrkrDefs::inttId>(track);
0189     trackStruct.nclusters_tpc = get_clusters<TrkrDefs::tpcId>(track);
0190     trackStruct.nclusters_micromegas = get_clusters<TrkrDefs::micromegasId>(track);
0191 
0192     trackStruct.chisquare = track->get_chisq();
0193     trackStruct.ndf = track->get_ndf();
0194 
0195     trackStruct.x = track->get_x();
0196     trackStruct.y = track->get_y();
0197     trackStruct.z = track->get_z();
0198     trackStruct.r = get_r(trackStruct.x, trackStruct.y);
0199     trackStruct.phi = std::atan2(trackStruct.y, trackStruct.x);
0200 
0201     trackStruct.px = track->get_px();
0202     trackStruct.py = track->get_py();
0203     trackStruct.pz = track->get_pz();
0204     trackStruct.pt = get_pt(trackStruct.px, trackStruct.py);
0205     trackStruct.p = get_p(trackStruct.px, trackStruct.py, trackStruct.pz);
0206     trackStruct.eta = get_eta(trackStruct.p, trackStruct.pz);
0207 
0208     return trackStruct;
0209   }
0210 #ifdef JUNK
0211   //! create cluster struct from svx cluster
0212   TrackEvaluationContainerv1::ClusterStruct create_cluster(TrkrDefs::cluskey key, TrkrCluster* cluster)
0213   {
0214     TrackEvaluationContainerv1::ClusterStruct cluster_struct;
0215     cluster_struct.layer = TrkrDefs::getLayer(key);
0216     cluster_struct.x = cluster->getX();
0217     cluster_struct.y = cluster->getY();
0218     cluster_struct.z = cluster->getZ();
0219     cluster_struct.r = get_r(cluster_struct.x, cluster_struct.y);
0220     cluster_struct.phi = std::atan2(cluster_struct.y, cluster_struct.x);
0221     cluster_struct.phi_error = cluster->getPhiError();
0222     cluster_struct.z_error = cluster->getZError();
0223     std::cout << " (x|y|z|r|l) "
0224               << cluster_struct.x << " | "
0225               << cluster_struct.y << " | "
0226               << cluster_struct.z << " | "
0227               << cluster_struct.r << " | "
0228               << cluster_struct.layer << " | "
0229               << std::endl;
0230     std::cout << " (xl|yl) "
0231               << cluster->getLocalX() << " | "
0232               << cluster->getLocalY()
0233               << std::endl;
0234     return cluster_struct;
0235   }
0236 
0237   //! add track information
0238   void add_trk_information(TrackEvaluationContainerv1::ClusterStruct& cluster, SvtxTrackState* state)
0239   {
0240     // need to extrapolate the track state to the right cluster radius to get proper residuals
0241     const auto trk_r = get_r(state->get_x(), state->get_y());
0242     const auto dr = cluster.r - trk_r;
0243     const auto trk_drdt = (state->get_x() * state->get_px() + state->get_y() * state->get_py()) / trk_r;
0244     const auto trk_dxdr = state->get_px() / trk_drdt;
0245     const auto trk_dydr = state->get_py() / trk_drdt;
0246     const auto trk_dzdr = state->get_pz() / trk_drdt;
0247 
0248     // store state position
0249     cluster.trk_x = state->get_x() + dr * trk_dxdr;
0250     cluster.trk_y = state->get_y() + dr * trk_dydr;
0251     cluster.trk_z = state->get_z() + dr * trk_dzdr;
0252     cluster.trk_r = get_r(cluster.trk_x, cluster.trk_y);
0253     cluster.trk_phi = std::atan2(cluster.trk_y, cluster.trk_x);
0254 
0255     /* store local momentum information */
0256     cluster.trk_px = state->get_px();
0257     cluster.trk_py = state->get_py();
0258     cluster.trk_pz = state->get_pz();
0259 
0260     /*
0261     store state angles in (r,phi) and (r,z) plans
0262     they are needed to study space charge distortions
0263     */
0264     const auto cosphi(std::cos(cluster.trk_phi));
0265     const auto sinphi(std::sin(cluster.trk_phi));
0266     const auto trk_pphi = -state->get_px() * sinphi + state->get_py() * cosphi;
0267     const auto trk_pr = state->get_px() * cosphi + state->get_py() * sinphi;
0268     const auto trk_pz = state->get_pz();
0269     cluster.trk_alpha = std::atan2(trk_pphi, trk_pr);
0270     cluster.trk_beta = std::atan2(trk_pz, trk_pr);
0271     cluster.trk_phi_error = state->get_phi_error();
0272     cluster.trk_z_error = state->get_z_error();
0273   }
0274 
0275   //! number of hits associated to cluster
0276   void add_cluster_size(TrackEvaluationContainerv1::ClusterStruct& cluster, TrkrDefs::cluskey clus_key, TrkrClusterHitAssoc* cluster_hit_map)
0277   {
0278     if (!cluster_hit_map) return;
0279     const auto range = cluster_hit_map->getHits(clus_key);
0280 
0281     // store full size
0282     cluster.size = std::distance(range.first, range.second);
0283 
0284     const auto detId = TrkrDefs::getTrkrId(clus_key);
0285     if (detId == TrkrDefs::micromegasId)
0286     {
0287       // for micromegas the directional cluster size depends on segmentation type
0288       auto segmentation_type = MicromegasDefs::getSegmentationType(clus_key);
0289       if (segmentation_type == MicromegasDefs::SegmentationType::SEGMENTATION_Z)
0290         cluster.z_size = cluster.size;
0291       else
0292         cluster.phi_size = cluster.size;
0293     }
0294     else
0295     {
0296       // for other detectors, one must loop over the constituting hits
0297       std::set<int> phibins;
0298       std::set<int> zbins;
0299       for (const auto& [first, hit_key] : range_adaptor(range))
0300       {
0301         switch (detId)
0302         {
0303         default:
0304           break;
0305         case TrkrDefs::mvtxId:
0306         {
0307           phibins.insert(MvtxDefs::getRow(hit_key));
0308           zbins.insert(MvtxDefs::getCol(hit_key));
0309           break;
0310         }
0311         case TrkrDefs::inttId:
0312         {
0313           phibins.insert(InttDefs::getRow(hit_key));
0314           zbins.insert(InttDefs::getCol(hit_key));
0315           break;
0316         }
0317         case TrkrDefs::tpcId:
0318         {
0319           phibins.insert(TpcDefs::getPad(hit_key));
0320           zbins.insert(TpcDefs::getTBin(hit_key));
0321           break;
0322         }
0323         }
0324       }
0325       cluster.phi_size = phibins.size();
0326       cluster.z_size = zbins.size();
0327     }
0328   }
0329 
0330   //! hit energy for a given cluster
0331   void add_cluster_energy(TrackEvaluationContainerv1::ClusterStruct& cluster, TrkrDefs::cluskey clus_key,
0332                           TrkrClusterHitAssoc* cluster_hit_map,
0333                           TrkrHitSetContainer* hitsetcontainer)
0334   {
0335     // check container
0336     if (!(cluster_hit_map && hitsetcontainer)) return;
0337 
0338     // for now this is only filled for micromegas
0339     const auto detId = TrkrDefs::getTrkrId(clus_key);
0340     if (detId != TrkrDefs::micromegasId) return;
0341 
0342     const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(clus_key);
0343     const auto hitset = hitsetcontainer->findHitSet(hitset_key);
0344     if (!hitset) return;
0345 
0346     const auto range = cluster_hit_map->getHits(clus_key);
0347     cluster.energy_max = 0;
0348     cluster.energy_sum = 0;
0349 
0350     for (const auto& pair : range_adaptor(range))
0351     {
0352       const auto hit = hitset->getHit(pair.second);
0353       if (hit)
0354       {
0355         const auto energy = hit->getEnergy();
0356         cluster.energy_sum += energy;
0357         if (energy > cluster.energy_max) cluster.energy_max = energy;
0358       }
0359     }
0360   }
0361 
0362   // add truth information
0363   void add_truth_information(TrackEvaluationContainerv1::ClusterStruct& cluster, std::set<PHG4Hit*> hits)
0364   {
0365     // need to extrapolate g4hits position to the cluster r
0366     /* this is done by using a linear extrapolation, with a straight line going through all provided G4Hits in and out positions */
0367     const auto rextrap = cluster.r;
0368     cluster.truth_size = hits.size();
0369     std::cout << "inter"
0370               << "\n";
0371 
0372     cluster.truth_x = interpolate<&PHG4Hit::get_x>(hits, rextrap);
0373     cluster.truth_y = interpolate<&PHG4Hit::get_y>(hits, rextrap);
0374     cluster.truth_z = interpolate<&PHG4Hit::get_z>(hits, rextrap);
0375     cluster.truth_r = get_r(cluster.truth_x, cluster.truth_y);
0376     cluster.truth_phi = std::atan2(cluster.truth_y, cluster.truth_x);
0377 
0378     /* add truth momentum information */
0379     cluster.truth_px = interpolate<&PHG4Hit::get_px>(hits, rextrap);
0380     cluster.truth_py = interpolate<&PHG4Hit::get_py>(hits, rextrap);
0381     cluster.truth_pz = interpolate<&PHG4Hit::get_pz>(hits, rextrap);
0382 
0383     std::cout << "inter2"
0384               << "\n";
0385 
0386     /*
0387     store state angles in (r,phi) and (r,z) plans
0388     they are needed to study space charge distortions
0389     */
0390     const auto cosphi(std::cos(cluster.truth_phi));
0391     const auto sinphi(std::sin(cluster.truth_phi));
0392     const auto truth_pphi = -cluster.truth_px * sinphi + cluster.truth_py * cosphi;
0393     const auto truth_pr = cluster.truth_px * cosphi + cluster.truth_py * sinphi;
0394 
0395     cluster.truth_alpha = std::atan2(truth_pphi, truth_pr);
0396     cluster.truth_beta = std::atan2(cluster.truth_pz, truth_pr);
0397     if (std::isnan(cluster.truth_alpha) || std::isnan(cluster.truth_beta))
0398     {
0399       // recalculate
0400       double truth_alpha = 0;
0401       double truth_beta = 0;
0402       double sum_w = 0;
0403       for (const auto& hit : hits)
0404       {
0405         const auto px = hit->get_x(1) - hit->get_x(0);
0406         const auto py = hit->get_y(1) - hit->get_y(0);
0407         const auto pz = hit->get_z(1) - hit->get_z(0);
0408         const auto pphi = -px * sinphi + py * cosphi;
0409         const auto pr = px * cosphi + py * sinphi;
0410 
0411         const auto w = hit->get_edep();
0412         if (w < 0) continue;
0413 
0414         sum_w += w;
0415         truth_alpha += w * std::atan2(pphi, pr);
0416         truth_beta += w * std::atan2(pz, pr);
0417       }
0418       truth_alpha /= sum_w;
0419       truth_beta /= sum_w;
0420       cluster.truth_alpha = truth_alpha;
0421       cluster.truth_beta = truth_beta;
0422     }
0423   }
0424 
0425   // ad}d truth information
0426   void add_truth_information(TrackEvaluationContainerv1::TrackStruct& track, PHG4Particle* particle)
0427   {
0428     if (particle)
0429     {
0430       track.is_primary = is_primary(particle);
0431       track.pid = particle->get_pid();
0432       track.truth_px = particle->get_px();
0433       track.truth_py = particle->get_py();
0434       track.truth_pz = particle->get_pz();
0435       track.truth_pt = get_pt(track.truth_px, track.truth_py);
0436       track.truth_p = get_p(track.truth_px, track.truth_py, track.truth_pz);
0437       track.truth_eta = get_eta(track.truth_p, track.truth_pz);
0438     }
0439   }
0440 #endif
0441 }  // namespace
0442 
0443 //_____________________________________________________________________
0444 DSTEmulator::DSTEmulator(const std::string& name, const std::string& filename, int inBits,
0445                          int inSabotage, bool compress)
0446   : SubsysReco(name)
0447   , _filename(filename)
0448   , _tfile(nullptr)
0449   , nBits(inBits)
0450   , sabotage(inSabotage)
0451   , apply_compression(compress)
0452 {
0453   rnd.SetSeed(PHRandomSeed());
0454 }
0455 
0456 //_____________________________________________________________________
0457 int DSTEmulator::Init(PHCompositeNode* topNode)
0458 {
0459   delete _tfile;
0460 
0461   _tfile = new TFile(_filename.c_str(), "RECREATE");
0462 
0463   _dst_data = new TNtuple("dst_data", "dst data",
0464                           "event:seed:"
0465                           "c3x:c3y:c3z:c3p:t3x:t3y:t3z:"
0466                           "c2x:c2y:c2r:c2l:t2x:t2y:t2r:t2l:"
0467                           "d2x:d2y:dr:"
0468                           "cmp_d2x:cmp_d2y:"
0469                           "pt:eta:phi:charge");
0470 
0471   // find DST node
0472   PHNodeIterator iter(topNode);
0473   auto* dstNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "DST"));
0474   if (!dstNode)
0475   {
0476     std::cout << "DSTEmulator::Init - DST Node missing" << std::endl;
0477     return Fun4AllReturnCodes::ABORTEVENT;
0478   }
0479 
0480   // get EVAL node
0481   iter = PHNodeIterator(dstNode);
0482   auto* evalNode = dynamic_cast<PHCompositeNode*>(iter.findFirst("PHCompositeNode", "EVAL"));
0483   if (!evalNode)
0484   {
0485     // create
0486     std::cout << "DSTEmulator::Init - EVAL node missing - creating" << std::endl;
0487     evalNode = new PHCompositeNode("EVAL");
0488     dstNode->addNode(evalNode);
0489   }
0490 
0491   auto* newNode = new PHIODataNode<PHObject>(new TrackEvaluationContainerv1, "TrackEvaluationContainer", "PHObject");
0492   evalNode->addNode(newNode);
0493 
0494   // m_compressor = new DSTCompressor(4.08407e-02,
0495   //                                  7.46530e-02,
0496   //                                  5.14381e-05,
0497   //                                  2.06291e-01,
0498   // with new distortion setting
0499   // m_compressor = new DSTCompressor(-2.70072e-02,
0500   //                                  2.49574e-02,
0501   //                                  1.12803e-03,
0502   //                                  5.91965e-02,
0503   // with mininum layer 7
0504   m_compressor = new DSTCompressor(6.96257e-04,
0505                                    3.16806e-02,
0506                                    7.32860e-05,
0507                                    5.93230e-02,
0508                                    nBits,
0509                                    nBits);
0510   return Fun4AllReturnCodes::EVENT_OK;
0511 }
0512 
0513 //_____________________________________________________________________
0514 int DSTEmulator::InitRun(PHCompositeNode* /*unused*/)
0515 {
0516   return Fun4AllReturnCodes::EVENT_OK;
0517 }
0518 
0519 //_____________________________________________________________________
0520 int DSTEmulator::process_event(PHCompositeNode* topNode)
0521 {
0522   // load nodes
0523   auto res = load_nodes(topNode);
0524   if (res != Fun4AllReturnCodes::EVENT_OK)
0525   {
0526     return res;
0527   }
0528 
0529   m_tGeometry = findNode::getClass<ActsGeometry>(topNode,
0530                                                  "ActsGeometry");
0531   if (!m_tGeometry)
0532   {
0533     std::cout << PHWHERE
0534               << "ActsTrackingGeometry not found on node tree. Exiting"
0535               << std::endl;
0536     return Fun4AllReturnCodes::ABORTRUN;
0537   }
0538 
0539   // cleanup output container
0540   if (m_container)
0541   {
0542     m_container->Reset();
0543   }
0544 
0545   evaluate_tracks();
0546 
0547   // clear maps
0548   m_g4hit_map.clear();
0549   return Fun4AllReturnCodes::EVENT_OK;
0550 }
0551 
0552 //_____________________________________________________________________
0553 int DSTEmulator::End(PHCompositeNode* /*unused*/)
0554 {
0555   _tfile->cd();
0556   _dst_data->Write();
0557   _tfile->Close();
0558   delete _tfile;
0559   return Fun4AllReturnCodes::EVENT_OK;
0560 }
0561 
0562 Acts::Vector3 DSTEmulator::getGlobalPosition(TrkrDefs::cluskey key, TrkrCluster* cluster) const
0563 {
0564   // get global position from Acts transform
0565   Acts::Vector3 globalpos;
0566   globalpos = m_tGeometry->getGlobalPosition(key, cluster);
0567 
0568   // check if TPC distortion correction are in place and apply
0569   // if( m_dcc ) { globalpos = m_distortionCorrection.get_corrected_position( globalpos, m_dcc ); }
0570 
0571   return globalpos;
0572 }
0573 
0574 //_____________________________________________________________________
0575 int DSTEmulator::load_nodes(PHCompositeNode* topNode)
0576 {
0577   // get necessary nodes
0578   m_track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMapTPCOnly");
0579   if (m_track_map)
0580   {
0581     std::cout << " DSTEmulator: Using TPC Only Track Map node " << std::endl;
0582   }
0583   else
0584   {
0585     std::cout << " DSTEmulator: TPC Only Track Map node not found, using default" << std::endl;
0586     m_track_map = findNode::getClass<SvtxTrackMap>(topNode, "SvtxTrackMap");
0587   }
0588   // cluster map
0589 
0590   m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
0591   if (m_cluster_map)
0592   {
0593     if (m_cluster_map->size() > 0)
0594     {
0595       std::cout << " DSTEmulator: Using CORRECTED_TRKR_CLUSTER node " << std::endl;
0596     }
0597     else
0598     {
0599       std::cout << " DSTEmulator: CORRECTED_TRKR_CLUSTER node not found, using TRKR_CLUSTER" << std::endl;
0600       m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0601     }
0602   }
0603   else
0604   {
0605     std::cout << " DSTEmulator: CORRECTED_TRKR_CLUSTER node not found at all, using TRKR_CLUSTER" << std::endl;
0606     m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0607   }
0608 
0609   // cluster hit association map
0610   m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0611 
0612   // cluster hit association map
0613   m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0614 
0615   // local container
0616   m_container = findNode::getClass<TrackEvaluationContainerv1>(topNode, "TrackEvaluationContainer");
0617 
0618   // hitset container
0619   m_hitsetcontainer = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0620 
0621   // g4hits
0622   m_g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
0623   m_g4hits_intt = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
0624   m_g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
0625   m_g4hits_micromegas = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
0626 
0627   // g4 truth info
0628   m_g4truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0629 
0630   return Fun4AllReturnCodes::EVENT_OK;
0631 }
0632 
0633 //_____________________________________________________________________
0634 void DSTEmulator::evaluate_tracks()
0635 {
0636   if (!(m_track_map && m_cluster_map && m_container))
0637   {
0638     return;
0639   }
0640 
0641   int iseed = 0;
0642   recoConsts* rc = recoConsts::instance();
0643   if (rc->FlagExist("RANDOMSEED"))
0644   {
0645     iseed = rc->get_IntFlag("RANDOMSEED");
0646   }
0647 
0648   // clear array
0649   m_container->clearTracks();
0650 
0651   for (const auto& trackpair : *m_track_map)
0652   {
0653     auto* const track = trackpair.second;
0654     auto track_struct = create_track(track);
0655 
0656     // truth information
0657     const auto [id, contributors] = get_max_contributor(track);
0658     track_struct.contributors = contributors;
0659 
0660     // get particle
0661     auto* particle = m_g4truthinfo->GetParticle(id);
0662     track_struct.embed = get_embed(particle);
0663     //    add_truth_information(track_struct, particle);
0664 
0665     // running iterator over track states, used to match a given cluster to a track state
0666     auto state_iter = track->begin_states();
0667 
0668     // loop over clusters
0669     for (auto key_iter = track->begin_cluster_keys(); key_iter != track->end_cluster_keys(); ++key_iter)
0670     {
0671       const auto& cluster_key = *key_iter;
0672       auto* cluster = m_cluster_map->findCluster(cluster_key);
0673       TrkrDefs::hitsetkey hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
0674 
0675       if (!cluster)
0676       {
0677         std::cout << "DSTEmulator::evaluate_tracks - unable to find cluster for key " << cluster_key << std::endl;
0678         continue;
0679       }
0680 
0681       if (TrkrDefs::getTrkrId(cluster_key) != TrkrDefs::tpcId)
0682       {
0683         continue;
0684       }
0685 
0686       // create new cluster struct
0687       // auto cluster_struct = create_cluster( cluster_key, cluster );
0688 
0689       // find track state that is the closest to cluster
0690       /* this assumes that both clusters and states are sorted along r */
0691       //     const auto radius( cluster_struct.r );
0692       const Acts::Vector3 globalpos_d = getGlobalPosition(cluster_key, cluster);
0693       float radius = get_r(globalpos_d[0], globalpos_d[1]);
0694       float clu_phi = std::atan2(globalpos_d[0], globalpos_d[1]);
0695       std::cout << "radius " << radius << std::endl;
0696       float dr_min = -1;
0697       for (auto iter = state_iter; iter != track->end_states(); ++iter)
0698       {
0699         const auto dr = std::abs(radius - get_r(iter->second->get_x(), iter->second->get_y()));
0700         if (dr_min < 0 || dr < dr_min)
0701         {
0702           state_iter = iter;
0703           dr_min = dr;
0704         }
0705         else
0706         {
0707           break;
0708         }
0709       }
0710       // Got cluster and got state
0711       /*
0712        */
0713 
0714       std::cout << "NEW (xg|yg) "
0715                 << globalpos_d[0] << " | "
0716                 << globalpos_d[1]
0717                 << std::endl;
0718       std::cout << "NEW (xl|yl) "
0719                 << cluster->getLocalX() << " | "
0720                 << cluster->getLocalY()
0721                 << std::endl;
0722 
0723       // state to local
0724       //  store track state in cluster struct
0725       //       add_trk_information( cluster_struct, state_iter->second );
0726       //    void add_trk_information( TrackEvaluationContainerv1::ClusterStruct& cluster, SvtxTrackState* state )
0727       //  need to extrapolate the track state to the right cluster radius to get proper residuals
0728       const auto trk_r = get_r(state_iter->second->get_x(), state_iter->second->get_y());
0729       std::cout << " trk_r  " << trk_r << std::endl;
0730       const auto dr = get_r(globalpos_d[0], globalpos_d[1]) - trk_r;
0731       std::cout << " dr  " << dr << std::endl;
0732       const auto trk_drdt = (state_iter->second->get_x() * state_iter->second->get_px() + state_iter->second->get_y() * state_iter->second->get_py()) / trk_r;
0733       std::cout << " trk_drdt  " << trk_drdt << std::endl;
0734       const auto trk_dxdr = state_iter->second->get_px() / trk_drdt;
0735       std::cout << " trk_dxdr " << trk_dxdr << std::endl;
0736       const auto trk_dydr = state_iter->second->get_py() / trk_drdt;
0737       std::cout << " trk_dydr  " << trk_dydr << std::endl;
0738       const auto trk_dzdr = state_iter->second->get_pz() / trk_drdt;
0739       std::cout << " trk_dzdr  " << trk_dzdr << std::endl;
0740 
0741       // store state position
0742       /*    */
0743       float trk_x = state_iter->second->get_x() + dr * trk_dxdr;
0744       float trk_y = state_iter->second->get_y() + dr * trk_dydr;
0745       float trk_z = state_iter->second->get_z() + dr * trk_dzdr;
0746       std::cout << "trk_x " << state_iter->second->get_x() << "trk_y" << state_iter->second->get_y() << "trk_z " << state_iter->second->get_z() << std::endl;
0747       //      float trk_r = get_r( trk_x, trk_y );
0748       // cluster.trk_phi = std::atan2( cluster.trk_y, cluster.trk_x );
0749 
0750       /* store local momentum information
0751          cluster.trk_px = state->get_px();
0752          cluster.trk_py = state->get_py();
0753          cluster.trk_pz = state->get_pz();
0754       */
0755       auto layer = TrkrDefs::getLayer(cluster_key);
0756       /*
0757       std::cout << " track 3D (x|y|z|r|l) "
0758                 << trk_x << " | "
0759                 << trk_y << " | "
0760                 << trk_z << " | "
0761                 << trk_r << " | "
0762                 << layer << " | "
0763                 << std::endl;
0764 
0765        std::cout << " cluster 3D (x|y|z|r|l) "
0766               << cluster->getX() << " | "
0767               << cluster->getY() << " | "
0768               << cluster->getZ() << " | "
0769               << get_r(cluster->getX() ,cluster->getY() ) << " | "
0770               << layer << " | "
0771               << std::endl;
0772       */
0773       // Get Surface
0774       Acts::Vector3 global(trk_x, trk_y, trk_z);
0775       //    TrkrDefs::subsurfkey subsurfkey = cluster->getSubSurfKey();
0776 
0777       //    std::cout << " subsurfkey: " << subsurfkey << std::endl;
0778       auto mapIter = m_tGeometry->maps().m_tpcSurfaceMap.find(layer);
0779 
0780       if (mapIter == m_tGeometry->maps().m_tpcSurfaceMap.end())
0781       {
0782         std::cout << PHWHERE
0783                   << "Error: hitsetkey not found in clusterSurfaceMap, layer = " << trk_r  // layer
0784                   << " hitsetkey = "
0785                   << hitsetkey << std::endl;
0786         continue;  // nullptr;
0787       }
0788       std::cout << " g0: " << global[0] << " g1: " << global[1] << " g2:" << global[2] << std::endl;
0789       // double global_phi = std::atan2(global[1], global[0]);
0790       // double global_z = global[2];
0791 
0792       // Predict which surface index this phi and z will correspond to
0793       // assumes that the vector elements are ordered positive z, -pi to pi, then negative z, -pi to pi
0794       std::vector<Surface> surf_vec = mapIter->second;
0795 
0796       Acts::Vector3 world(globalpos_d[0], globalpos_d[1], globalpos_d[2]);
0797       double world_phi = std::atan2(world[1], world[0]);
0798       double world_z = world[2];
0799 
0800       double fraction = (world_phi + M_PI) / (2.0 * M_PI);
0801       double rounded_nsurf = round((double) (surf_vec.size() / 2) * fraction - 0.5);  // NOLINT(bugprone-integer-division)
0802       unsigned int nsurf = (unsigned int) rounded_nsurf;
0803       if (world_z < 0)
0804       {
0805         nsurf += surf_vec.size() / 2;
0806       }
0807 
0808       const Surface& surface = surf_vec[nsurf];
0809 
0810       Acts::Vector3 center = surface->center(m_tGeometry->geometry().getGeoContext()) / Acts::UnitConstants::cm;
0811 
0812       // no conversion needed, only used in acts
0813       //    Acts::Vector3 normal = surface->normal(m_tGeometry->getGeoContext());
0814       double TrkRadius = std::sqrt(trk_x * trk_x + trk_y * trk_y);
0815       double rTrkPhi = TrkRadius * std::atan2(trk_y, trk_x);  // trkphi;
0816       double surfRadius = std::sqrt(center(0) * center(0) + center(1) * center(1));
0817       double surfPhiCenter = std::atan2(center[1], center[0]);
0818       double surfRphiCenter = surfPhiCenter * surfRadius;
0819       double surfZCenter = center[2];
0820 
0821       double trklocX = 0;
0822       double trklocY = 0;
0823       float delta_r = 0;
0824 
0825       delta_r = TrkRadius - surfRadius;
0826       trklocX = rTrkPhi - surfRphiCenter;
0827       trklocY = trk_z - surfZCenter;
0828       /*
0829       std::cout << " clus 2D: (xl|yl) "
0830                 << cluster->getLocalX() << " | "
0831                 << cluster->getLocalY()
0832                 << std::endl;
0833 
0834       std::cout << " trk 2D: (xl|yl) "
0835                 << trklocX << " | "
0836                 << trklocY
0837                 << std::endl;
0838       */
0839       float delta_x = trklocX - cluster->getLocalX();
0840       float delta_y = trklocY - cluster->getLocalY();
0841 
0842       float comp_dx = compress_dx(delta_x);
0843       float comp_dy = compress_dy(delta_y);
0844 
0845       // Sabotage the tracking by multiplying the residuals with a random number
0846       if (sabotage == -1)
0847       {
0848         comp_dx = rnd.Uniform(-1, 1) * delta_x;
0849         comp_dy = rnd.Uniform(-1, 1) * delta_y;
0850       }
0851       else if (sabotage == -2)
0852       {
0853         comp_dx = rnd.Uniform(-10, 10);
0854         comp_dy = rnd.Uniform(-10, 10);
0855       }
0856 
0857       /* "dst_data", "dst data","event:seed:"
0858                             "c3x:c3y:c3z:t3x:t3y:t3z:"
0859                             "c2x:c2y:c2r:c2l:t2x:t2y:t2r:t2l:"
0860                             "d2x:d2y:"
0861                             "cmp_d2x:cmp_d2y:"
0862                             "pt:eta:phi"
0863       */
0864       float data[] = {1,
0865                       (float) iseed,
0866                       (float) globalpos_d[0],
0867                       (float) globalpos_d[1],
0868                       (float) globalpos_d[2],
0869                       clu_phi,
0870                       trk_x,
0871                       trk_y,
0872                       trk_z,
0873                       cluster->getLocalX(),
0874                       cluster->getLocalY(),
0875                       get_r((float) globalpos_d[0], (float) globalpos_d[1]),
0876                       (float) layer,
0877                       (float) trklocX,
0878                       (float) trklocY,
0879                       get_r(trk_x, trk_y),
0880                       (float) layer,
0881                       delta_x,
0882                       delta_y,
0883                       delta_r,
0884                       comp_dx,
0885                       comp_dy,
0886                       track_struct.pt,
0887                       track_struct.eta,
0888                       track_struct.phi,
0889                       (float) track_struct.charge};
0890       _dst_data->Fill(data);
0891       std::cout << "filled"
0892                 << "\n";
0893 
0894       if (apply_compression)
0895       {
0896         cluster->setLocalX(trklocX - comp_dx);
0897         cluster->setLocalY(trklocY - comp_dy);
0898       }
0899 
0900 #ifdef JUNK
0901       /*
0902       store state angles in (r,phi) and (r,z) plans
0903       they are needed to study space charge distortions
0904       */
0905       /*
0906     const auto cosphi( std::cos( cluster.trk_phi ) );
0907     const auto sinphi( std::sin( cluster.trk_phi ) );
0908     const auto trk_pphi = -state->get_px()*sinphi + state->get_py()*cosphi;
0909     const auto trk_pr = state->get_px()*cosphi + state->get_py()*sinphi;
0910     const auto trk_pz = state->get_pz();
0911     cluster.trk_alpha = std::atan2( trk_pphi, trk_pr );
0912     cluster.trk_beta = std::atan2( trk_pz, trk_pr );
0913     cluster.trk_phi_error = state->get_phi_error();
0914     cluster.trk_z_error = state->get_z_error();
0915       */
0916 #endif
0917     }
0918   }
0919 
0920   std::cout << "DSTEmulator::evaluate_tracks - tracks: " << m_container->tracks().size() << std::endl;
0921 }
0922 
0923 //_____________________________________________________________________
0924 float DSTEmulator::compress_dx(float in_val)
0925 {
0926   unsigned short key = m_compressor->compressPhi(in_val);
0927   float out_val = m_compressor->decompressPhi(key);
0928   return out_val;
0929 }
0930 
0931 //_____________________________________________________________________
0932 float DSTEmulator::compress_dy(float in_val)
0933 {
0934   unsigned short key = m_compressor->compressZ(in_val);
0935   float out_val = m_compressor->decompressZ(key);
0936   return out_val;
0937 }
0938 
0939 //_____________________________________________________________________
0940 DSTEmulator::G4HitSet DSTEmulator::find_g4hits(TrkrDefs::cluskey cluster_key) const
0941 {
0942   // check maps
0943   if (!(m_cluster_hit_map && m_hit_truth_map))
0944   {
0945     return G4HitSet();
0946   }
0947 
0948   // check if in map
0949   auto map_iter = m_g4hit_map.lower_bound(cluster_key);
0950   if (map_iter != m_g4hit_map.end() && cluster_key == map_iter->first)
0951   {
0952     return map_iter->second;
0953   }
0954 
0955   // find hitset associated to cluster
0956   G4HitSet out;
0957   const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
0958   for (const auto& [first, hit_key] : range_adaptor(m_cluster_hit_map->getHits(cluster_key)))
0959   {
0960     // store hits to g4hit associations
0961     TrkrHitTruthAssoc::MMap g4hit_map;
0962     m_hit_truth_map->getG4Hits(hitset_key, hit_key, g4hit_map);
0963 
0964     // find corresponding g4 hist
0965     for (auto& truth_iter : g4hit_map)
0966     {
0967       const auto g4hit_key = truth_iter.second.second;
0968       PHG4Hit* g4hit = nullptr;
0969 
0970       switch (TrkrDefs::getTrkrId(hitset_key))
0971       {
0972       case TrkrDefs::mvtxId:
0973         if (m_g4hits_mvtx)
0974         {
0975           g4hit = m_g4hits_mvtx->findHit(g4hit_key);
0976         }
0977         break;
0978 
0979       case TrkrDefs::inttId:
0980         if (m_g4hits_intt)
0981         {
0982           g4hit = m_g4hits_intt->findHit(g4hit_key);
0983         }
0984         break;
0985 
0986       case TrkrDefs::tpcId:
0987         if (m_g4hits_tpc)
0988         {
0989           g4hit = m_g4hits_tpc->findHit(g4hit_key);
0990         }
0991         break;
0992 
0993       case TrkrDefs::micromegasId:
0994         if (m_g4hits_micromegas)
0995         {
0996           g4hit = m_g4hits_micromegas->findHit(g4hit_key);
0997         }
0998         break;
0999 
1000       default:
1001         break;
1002       }
1003 
1004       if (g4hit)
1005       {
1006         out.insert(g4hit);
1007       }
1008       else
1009       {
1010         std::cout << "DSTEmulator::find_g4hits - g4hit not found " << g4hit_key << std::endl;
1011       }
1012     }
1013   }
1014 
1015   // insert in map and return
1016   return m_g4hit_map.insert(map_iter, std::make_pair(cluster_key, std::move(out)))->second;
1017 }
1018 
1019 //_____________________________________________________________________
1020 std::pair<int, int> DSTEmulator::get_max_contributor(SvtxTrack* track) const
1021 {
1022   if (!(m_track_map && m_cluster_map && m_g4truthinfo))
1023   {
1024     return {0, 0};
1025   }
1026 
1027   // maps MC track id and number of matching g4hits
1028   using IdMap = std::map<int, int>;
1029   IdMap contributor_map;
1030 
1031   // loop over clusters
1032   for (auto key_iter = track->begin_cluster_keys(); key_iter != track->end_cluster_keys(); ++key_iter)
1033   {
1034     const auto& cluster_key = *key_iter;
1035     for (const auto& hit : find_g4hits(cluster_key))
1036     {
1037       const int trkid = hit->get_trkid();
1038       auto iter = contributor_map.lower_bound(trkid);
1039       if (iter == contributor_map.end() || iter->first != trkid)
1040       {
1041         contributor_map.insert(iter, std::make_pair(trkid, 1));
1042       }
1043       else
1044       {
1045         ++iter->second;
1046       }
1047     }
1048   }
1049 
1050   if (contributor_map.empty())
1051   {
1052     return {0, 0};
1053   }
1054 
1055   return *std::max_element(
1056       contributor_map.cbegin(), contributor_map.cend(),
1057       [](const IdMap::value_type& first, const IdMap::value_type& second)
1058       { return first.second < second.second; });
1059 }
1060 
1061 //_____________________________________________________________________
1062 int DSTEmulator::get_embed(PHG4Particle* particle) const
1063 {
1064   return (m_g4truthinfo && particle) ? m_g4truthinfo->isEmbeded(particle->get_primary_id()) : 0;
1065 }
[ 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!