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