QA/SimulationModules/QAG4SimulationMvtx.cc

0001 #include "QAG4SimulationMvtx.h"
0002
0003 #include <qautils/QAHistManagerDef.h>
0004 #include <qautils/QAUtil.h>
0005
0006 #include <g4detectors/PHG4CylinderGeomContainer.h>
0007
0008 #include <g4main/PHG4Hit.h>
0009 #include <g4main/PHG4HitContainer.h>
0010
0011 #include <trackbase/ActsGeometry.h>
0012 #include <trackbase/ClusterErrorPara.h>
0013 #include <trackbase/MvtxDefs.h>
0014 #include <trackbase/TrkrCluster.h>
0015 #include <trackbase/TrkrClusterContainer.h>
0016 #include <trackbase/TrkrClusterHitAssoc.h>
0017 #include <trackbase/TrkrDefs.h>  // for getTrkrId, getHit...
0018 #include <trackbase/TrkrHitTruthAssoc.h>
0019
0020 #include <trackbase_historic/ActsTransformations.h>
0021
0022 #include <fun4all/Fun4AllHistoManager.h>
0023 #include <fun4all/Fun4AllReturnCodes.h>
0024 #include <fun4all/SubsysReco.h>  // for SubsysReco
0025
0026 #include <phool/getClass.h>
0027 #include <phool/phool.h>  // for PHWHERE
0028
0029 #include <TAxis.h>  // for TAxis
0030 #include <TH1.h>
0031
0032 #include <cassert>
0033 #include <cmath>  // for atan2
0034 #include <format>
0035 #include <iostream>  // for operator<<, basic...
0036 #include <iterator>  // for distance
0037 #include <map>       // for map
0038 #include <set>
0039 #include <string>
0040 #include <utility>  // for pair, make_pair
0041
0042 namespace
0043 {
0044
0045   //! range adaptor to be able to use range-based for loop
0046   template <class T>
0047   class range_adaptor
0048   {
0049    public:
0050     explicit range_adaptor(const T& range)
0051       : m_range(range)
0052     {
0053     }
0054     const typename T::first_type& begin() { return m_range.first; }
0055     const typename T::second_type& end() { return m_range.second; }
0056
0057    private:
0058     T m_range;
0059   };
0060
0061 }  // namespace
0062
0063 //________________________________________________________________________
0064 QAG4SimulationMvtx::QAG4SimulationMvtx(const std::string& name)
0065   : SubsysReco(name)
0066 {
0067 }
0068
0069 //________________________________________________________________________
0070 int QAG4SimulationMvtx::InitRun(PHCompositeNode* topNode)
0071 {
0072   // prevent multiple creations of histograms
0073   if (m_initialized)
0074   {
0075     return Fun4AllReturnCodes::EVENT_OK;
0076   }
0077
0078   m_initialized = true;
0079
0080   // find mvtx geometry
0081   auto* geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
0082   if (!geom_container)
0083   {
0084     std::cout << PHWHERE << " unable to find DST node CYLINDERGEOM_MVTX" << std::endl;
0085     return Fun4AllReturnCodes::ABORTRUN;
0086   }
0087
0088   // get layers from mvtx geometry
0089   const auto range = geom_container->get_begin_end();
0090   for (auto iter = range.first; iter != range.second; ++iter)
0091   {
0092     m_layers.insert(iter->first);
0093   }
0094
0095   // histogram manager
0096   auto* hm = QAHistManagerDef::getHistoManager();
0097   assert(hm);
0098
0099   // create histograms
0100   for (const auto& layer : m_layers)
0101   {
0102     if (Verbosity())
0103     {
0104       std::cout << PHWHERE << " adding layer " << layer << std::endl;
0105     }
0106     {
0107       // rphi residuals (cluster - truth)
0108       auto* h = new TH1F(std::format("{}drphi_{}", get_histo_prefix(), layer).c_str(), std::format("MVTX r#Delta#phi_{{cluster-truth}} layer_{}", layer).c_str(), 100, -2e-3, 2e-3);
0109       h->GetXaxis()->SetTitle("r#Delta#phi_{cluster-truth} (cm)");
0110       hm->registerHisto(h);
0111     }
0112
0113     {
0114       // rphi cluster errors
0115       auto* h = new TH1F(std::format("{}rphi_error_{}", get_histo_prefix(), layer).c_str(), std::format("MVTX r#Delta#phi error layer_{}", layer).c_str(), 100, 0, 2e-3);
0116       h->GetXaxis()->SetTitle("r#Delta#phi error (cm)");
0117       hm->registerHisto(h);
0118     }
0119
0120     {
0121       // phi pulls (cluster - truth)
0122       auto* h = new TH1F(std::format("{}phi_pulls_{}", get_histo_prefix(), layer).c_str(), std::format("MVTX #Delta#phi_{{cluster-truth}}/#sigma#phi layer_{}", layer).c_str(), 100, -3, 3);
0123       h->GetXaxis()->SetTitle("#Delta#phi_{cluster-truth}/#sigma#phi");
0124       hm->registerHisto(h);
0125     }
0126
0127     {
0128       // z residuals (cluster - truth)
0129       auto* h = new TH1F(std::format("{}dz_{}", get_histo_prefix(), layer).c_str(), std::format("MVTX #Deltaz_{{cluster-truth}} layer_{}", layer).c_str(), 100, -3e-3, 3e-3);
0130       h->GetXaxis()->SetTitle("#Deltaz_{cluster-truth} (cm)");
0131       hm->registerHisto(h);
0132     }
0133
0134     {
0135       // z cluster errors
0136       auto* h = new TH1F(std::format("{}z_error_{}", get_histo_prefix(), layer).c_str(), std::format("MVTX z error layer_{}", layer).c_str(), 100, 0, 3e-3);
0137       h->GetXaxis()->SetTitle("z error (cm)");
0138       hm->registerHisto(h);
0139     }
0140
0141     {
0142       // z pulls (cluster - truth)
0143       auto* h = new TH1F(std::format("{}z_pulls_{}", get_histo_prefix(), layer).c_str(), std::format("MVTX #Deltaz_{{cluster-truth}}/#sigmaz layer_{}", layer).c_str(), 100, -3, 3);
0144       h->GetXaxis()->SetTitle("#Deltaz_{cluster-truth}/#sigmaz");
0145       hm->registerHisto(h);
0146     }
0147
0148     {
0149       // total cluster size
0150       auto* h = new TH1F(std::format("{}clus_size_{}", get_histo_prefix(), layer).c_str(), std::format("MVTX cluster size layer_{}", layer).c_str(), 20, 0, 20);
0151       h->GetXaxis()->SetTitle("csize");
0152       hm->registerHisto(h);
0153     }
0154
0155     {
0156       // cluster size in phi
0157       auto* h = new TH1F(std::format("{}clus_size_phi_{}", get_histo_prefix(), layer).c_str(), std::format("MVTX cluster size (#phi) layer_{}", layer).c_str(), 10, 0, 10);
0158       h->GetXaxis()->SetTitle("csize_{#phi}");
0159       hm->registerHisto(h);
0160     }
0161
0162     {
0163       // cluster size in z
0164       auto* h = new TH1F(std::format("{}clus_size_z_{}", get_histo_prefix(), layer).c_str(), std::format("MVTX cluster size (z) layer_{}", layer).c_str(), 10, 0, 10);
0165       h->GetXaxis()->SetTitle("csize_{z}");
0166       hm->registerHisto(h);
0167     }
0168   }
0169
0170   return Fun4AllReturnCodes::EVENT_OK;
0171 }
0172
0173 //_____________________________________________________________________
0174 int QAG4SimulationMvtx::process_event(PHCompositeNode* topNode)
0175 {
0176   // load nodes
0177   auto res = load_nodes(topNode);
0178   if (res != Fun4AllReturnCodes::EVENT_OK)
0179   {
0180     return res;
0181   }
0182   // run evaluation
0183   evaluate_clusters();
0184   return Fun4AllReturnCodes::EVENT_OK;
0185 }
0186
0187 //________________________________________________________________________
0188 std::string QAG4SimulationMvtx::get_histo_prefix() const
0189 {
0190   return std::string("h_") + Name() + std::string("_");
0191 }
0192
0193 //________________________________________________________________________
0194 int QAG4SimulationMvtx::load_nodes(PHCompositeNode* topNode)
0195 {
0196   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0197   if (!m_tGeometry)
0198   {
0199     std::cout << PHWHERE << "No acts tracking geometry, exiting."
0200               << std::endl;
0201     return Fun4AllReturnCodes::ABORTEVENT;
0202   }
0203
0204   m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0205   if (!m_cluster_map)
0206   {
0207     std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTER" << std::endl;
0208     return Fun4AllReturnCodes::ABORTEVENT;
0209   }
0210
0211   m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0212   if (!m_cluster_hit_map)
0213   {
0214     std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTERHITASSOC" << std::endl;
0215     return Fun4AllReturnCodes::ABORTEVENT;
0216   }
0217
0218   m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0219   if (!m_hit_truth_map)
0220   {
0221     std::cout << PHWHERE << " unable to find DST node TRKR_HITTRUTHASSOC" << std::endl;
0222     return Fun4AllReturnCodes::ABORTEVENT;
0223   }
0224
0225   m_g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
0226   if (!m_g4hits_mvtx)
0227   {
0228     std::cout << PHWHERE << " unable to find DST node G4HIT_MVTX" << std::endl;
0229     return Fun4AllReturnCodes::ABORTEVENT;
0230   }
0231
0232   return Fun4AllReturnCodes::EVENT_OK;
0233 }
0234
0235 //________________________________________________________________________
0236 void QAG4SimulationMvtx::evaluate_clusters()
0237 {
0238   // histogram manager
0239   auto* hm = QAHistManagerDef::getHistoManager();
0240   assert(hm);
0241
0242   // load relevant histograms
0243   struct HistogramList
0244   {
0245     TH1* drphi = nullptr;
0246     TH1* rphi_error = nullptr;
0247     TH1* phi_pulls = nullptr;
0248
0249     TH1* dz = nullptr;
0250     TH1* z_error = nullptr;
0251     TH1* z_pulls = nullptr;
0252
0253     TH1* csize = nullptr;
0254     TH1* csize_phi = nullptr;
0255     TH1* csize_z = nullptr;
0256   };
0257
0258   using HistogramMap = std::map<int, HistogramList>;
0259   HistogramMap histograms;
0260
0261   for (const auto& layer : m_layers)
0262   {
0263     HistogramList h;
0264     h.drphi = dynamic_cast<TH1*>(hm->getHisto(std::format("{}drphi_{}", get_histo_prefix(), layer)));
0265     h.rphi_error = dynamic_cast<TH1*>(hm->getHisto(std::format("{}rphi_error_{}", get_histo_prefix(), layer)));
0266     h.phi_pulls = dynamic_cast<TH1*>(hm->getHisto(std::format("{}phi_pulls_{}", get_histo_prefix(), layer)));
0267
0268     h.dz = dynamic_cast<TH1*>(hm->getHisto(std::format("{}dz_{}", get_histo_prefix(), layer)));
0269     h.z_error = dynamic_cast<TH1*>(hm->getHisto(std::format("{}z_error_{}", get_histo_prefix(), layer)));
0270     h.z_pulls = dynamic_cast<TH1*>(hm->getHisto(std::format("{}z_pulls_{}", get_histo_prefix(), layer)));
0271
0272     h.csize = dynamic_cast<TH1*>(hm->getHisto(std::format("{}clus_size_{}", get_histo_prefix(), layer)));
0273     h.csize_phi = dynamic_cast<TH1*>(hm->getHisto(std::format("{}clus_size_phi_{}", get_histo_prefix(), layer)));
0274     h.csize_z = dynamic_cast<TH1*>(hm->getHisto(std::format("{}clus_size_z_{}", get_histo_prefix(), layer)));
0275
0276     histograms.insert(std::make_pair(layer, h));
0277   }
0278
0279   for (const auto& hitsetkey : m_cluster_map->getHitSetKeys(TrkrDefs::TrkrId::mvtxId))
0280   {
0281     auto range = m_cluster_map->getClusters(hitsetkey);
0282     for (auto clusterIter = range.first; clusterIter != range.second; ++clusterIter)
0283     {
0284       // get cluster key, detector id and check
0285       const auto& key = clusterIter->first;
0286       // get cluster
0287       const auto& cluster = clusterIter->second;
0288       const auto global = m_tGeometry->getGlobalPosition(key, cluster);
0289
0290       // get relevant cluster information
0291       const auto r_cluster = QAG4Util::get_r(global(0), global(1));
0292       const auto z_cluster = global(2);
0293       const auto phi_cluster = (float) std::atan2(global(1), global(0));
0294
0295       double const phi_error = cluster->getRPhiError() / r_cluster;
0296       double const z_error = cluster->getZError();
0297
0298       // find associated g4hits
0299       const auto g4hits = find_g4hits(key);
0300
0301       // get relevant truth information
0302       const auto x_truth = QAG4Util::interpolate<&PHG4Hit::get_x>(g4hits, r_cluster);
0303       const auto y_truth = QAG4Util::interpolate<&PHG4Hit::get_y>(g4hits, r_cluster);
0304       const auto z_truth = QAG4Util::interpolate<&PHG4Hit::get_z>(g4hits, r_cluster);
0305       const auto phi_truth = std::atan2(y_truth, x_truth);
0306
0307       const auto dphi = QAG4Util::delta_phi(phi_cluster, phi_truth);
0308       const auto dz = z_cluster - z_truth;
0309
0310       // get layer, get histograms
0311       const auto layer = TrkrDefs::getLayer(key);
0312       const auto hiter = histograms.find(layer);
0313       if (hiter == histograms.end())
0314       {
0315         continue;
0316       }
0317
0318       // fill phi residuals, errors and pulls
0319       auto fill = [](TH1* h, float value)
0320       { if( h ) { h->Fill( value ); } };
0321
0322       fill(hiter->second.drphi, r_cluster * dphi);
0323       fill(hiter->second.rphi_error, r_cluster * phi_error);
0324       fill(hiter->second.phi_pulls, dphi / phi_error);
0325
0326       // fill z residuals, errors and pulls
0327       fill(hiter->second.dz, dz);
0328       fill(hiter->second.z_error, z_error);
0329       fill(hiter->second.z_pulls, dz / z_error);
0330
0331       // cluster sizes
0332       // get associated hits
0333       const auto hit_range = m_cluster_hit_map->getHits(key);
0334       fill(hiter->second.csize, std::distance(hit_range.first, hit_range.second));
0335
0336       std::set<int> phibins;
0337       std::set<int> zbins;
0338       for (auto hit_iter = hit_range.first; hit_iter != hit_range.second; ++hit_iter)
0339       {
0340         const auto& hit_key = hit_iter->second;
0341         phibins.insert(MvtxDefs::getRow(hit_key));
0342         zbins.insert(MvtxDefs::getCol(hit_key));
0343       }
0344
0345       fill(hiter->second.csize_phi, phibins.size());
0346       fill(hiter->second.csize_z, zbins.size());
0347     }
0348   }
0349 }
0350 //_____________________________________________________________________
0351 QAG4SimulationMvtx::G4HitSet QAG4SimulationMvtx::find_g4hits(TrkrDefs::cluskey cluster_key) const
0352 {
0353   // find hitset associated to cluster
0354   G4HitSet out;
0355
0356   const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
0357
0358   // get detector id
0359   const auto trkrId = TrkrDefs::getTrkrId(hitset_key);
0360
0361   /*
0362    * for MVTX,
0363    * also get bare (== strobe 0) hitsetkey,
0364    * since this is the one recorded in the HitTruth association map
0365    */
0366   const auto bare_hitset_key =
0367       trkrId == TrkrDefs::TrkrId::mvtxId ? MvtxDefs::resetStrobe(hitset_key) : hitset_key;
0368
0369   // loop over hits associated to clusters
0370   const auto range = m_cluster_hit_map->getHits(cluster_key);
0371   for (const auto& [ckey, hit_key] : range_adaptor(range))
0372   {
0373     // store hits to g4hit associations
0374     TrkrHitTruthAssoc::MMap g4hit_map;
0375     m_hit_truth_map->getG4Hits(bare_hitset_key, hit_key, g4hit_map);
0376
0377     // find corresponding g4 hist
0378     for (auto& truth_iter : g4hit_map)
0379     {
0380       // g4hit key
0381       const auto g4hit_key = truth_iter.second.second;
0382
0383       // g4 hit
0384       PHG4Hit* g4hit = (TrkrDefs::getTrkrId(hitset_key) == TrkrDefs::mvtxId) ? m_g4hits_mvtx->findHit(g4hit_key) : nullptr;
0385
0386       // insert in set
0387       if (g4hit)
0388       {
0389         out.insert(g4hit);
0390       }
0391     }
0392   }
0393
0394   return out;
0395 }