QA/SimulationModules/QAG4SimulationIntt.cc

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