g4simulation/g4eval/g4evalfn.cc

0001 #include <g4tracking/EmbRecoMatchContainer.h>
0002
0003 #include <trackbase/ActsGeometry.h>
0004 #include <trackbase/TrkrCluster.h>
0005 #include <trackbase/TrkrClusterContainer.h>
0006 #include <trackbase_historic/SvtxTrackMap.h>
0007 #include "TrkrClusLoc.h"
0008 #include "TrkrClusterIsMatcher.h"
0009 #include "g4evalfn.h"
0010
0011 #include <algorithm>
0012 #include <cmath>
0013 #include <set>
0014
0015 namespace g4evalfn
0016 {
0017
0018   int trklayer_det(int layer)
0019   {
0020     if (layer < 3)
0021     {
0022       return DET::MVTX;
0023     }
0024     if (layer < 7)
0025     {
0026       return DET::INTT;
0027     }
0028     if (layer < 55)
0029     {
0030       return DET::TPC;
0031     }
0032     return DET::TPOT;
0033   }
0034
0035   int trklayer_det(TrkrDefs::hitsetkey key)
0036   {
0037     return trklayer_det(TrkrDefs::getLayer(key));
0038   }
0039
0040   int trklayer_det(TrkrDefs::cluskey key)
0041   {
0042     return trklayer_det(TrkrDefs::getLayer(key));
0043   }
0044
0045   std::vector<int> unmatchedSvtxTrkIds(EmbRecoMatchContainer* matches, SvtxTrackMap* m_SvtxTrackMap)
0046   {
0047     std::set<unsigned int> ids_matched{};
0048     for (auto trkid : matches->ids_RecoMatched())
0049     {
0050       ids_matched.insert(trkid);
0051     }
0052
0053     std::set<int> ids_unmatched;
0054     for (auto& reco : *m_SvtxTrackMap)
0055     {
0056       auto trkid = reco.first;
0057       if (ids_matched.count(trkid) == 0)
0058       {
0059         ids_unmatched.insert(trkid);
0060       }
0061     }
0062     std::vector<int> ids_vec;
0063     ids_vec.reserve(ids_unmatched.size());
0064     for (auto id : ids_unmatched)
0065     {
0066       ids_vec.push_back((int) id);
0067     }
0068     std::sort(ids_vec.begin(), ids_vec.end());
0069     return ids_vec;
0070   }
0071
0072   TrkrClusLoc clusloc_PHG4(TrkrClusterIsMatcher* ismatcher, TrkrDefs::cluskey key)
0073   {
0074     auto cluster = ismatcher->m_TruthClusters->findCluster(key);
0075     Eigen::Vector3d gloc = ismatcher->m_ActsGeometry->getGlobalPosition(key, cluster);
0076     return {TrkrDefs::getLayer(key), gloc,
0077             cluster->getPosition(0), cluster->getPhiSize(), cluster->getPosition(1), cluster->getZSize(), key};
0078   }
0079
0080   TrkrClusLoc clusloc_SVTX(TrkrClusterIsMatcher* ismatcher, TrkrDefs::cluskey key)
0081   {
0082     auto cluster = ismatcher->m_RecoClusters->findCluster(key);
0083     Eigen::Vector3d gloc = ismatcher->m_ActsGeometry->getGlobalPosition(key, cluster);
0084     return {TrkrDefs::getLayer(key), gloc,
0085             cluster->getPosition(0), cluster->getPhiSize(), cluster->getPosition(1), cluster->getZSize(), key};
0086   }
0087
0088   float calc_match_statistic(TrkrClusterIsMatcher* ismatcher, TrkrDefs::cluskey key_A, TrkrDefs::cluskey key_B)
0089   {
0090     int layer = TrkrDefs::getLayer(key_A);
0091     auto det_0123 = trklayer_det(layer);
0092     auto clus_T = ismatcher->m_TruthClusters->findCluster(key_A);
0093     auto clus_R = ismatcher->m_RecoClusters->findCluster(key_B);
0094     auto dphi = abs_dphi(clus_T->getPosition(0), clus_R->getPosition(0));
0095
0096     float stat = 0.;
0097     switch (det_0123)
0098     {
0099     case 0:
0100     {  // MVTX
0101       if (ismatcher->single_pixel_phi_MVTX)
0102       {
0103         stat += dphi / ismatcher->tol_pitch_phi[layer];
0104       }
0105       else
0106       {
0107         stat += dphi / (ismatcher->tol_pitch_phi[layer] * std::max(clus_T->getPhiSize(), clus_R->getPhiSize()));
0108       }
0109       const float delta_z = std::fabs(clus_T->getPosition(1) - clus_R->getPosition(1));
0110       if (ismatcher->single_pixel_z_MVTX)
0111       {
0112         stat += delta_z / ismatcher->tol_pitch_z_MVTX;
0113       }
0114       else
0115       {
0116         stat += delta_z / (ismatcher->tol_pitch_z_MVTX * std::max(clus_T->getZSize(), clus_R->getZSize()));
0117       }
0118     }
0119     break;
0120     case 1:
0121     {
0122       if (ismatcher->single_pixel_phi_INTT)
0123       {
0124         stat += dphi / ismatcher->tol_pitch_phi[layer];
0125       }
0126       else
0127       {
0128         stat += dphi / (ismatcher->tol_pitch_phi[layer] * std::max(clus_T->getPhiSize(), clus_R->getPhiSize()));
0129       }
0130     }
0131     break;
0132     case 2:
0133     {  // TPC
0134       if (ismatcher->single_bin_phi_TPC)
0135       {
0136         stat += dphi / ismatcher->tol_pitch_phi[layer];
0137       }
0138       else
0139       {
0140         stat += dphi / (ismatcher->tol_pitch_phi[layer] * std::max(clus_T->getPhiSize(), clus_R->getPhiSize()));
0141       }
0142       const float delta_t = std::fabs(clus_T->getPosition(1) - clus_R->getPosition(1));
0143       if (ismatcher->single_bin_t_TPC)
0144       {
0145         stat += delta_t / ismatcher->tol_step_t_TPC;
0146       }
0147       else
0148       {
0149         stat += delta_t / (ismatcher->tol_step_t_TPC * std::max(clus_T->getZSize(), clus_R->getZSize()));
0150       }
0151     }
0152     break;
0153     // case 3:  // TPOT // empty case triggers clang-tidy warning
0154     //   break;
0155     default:
0156       break;
0157     }
0158     return stat;
0159   }
0160
0161 }  // namespace g4evalfn