QA/SimulationModules/QAG4SimulationTpc.cc

0001 #include "QAG4SimulationTpc.h"
0002
0003 #include <qautils/QAHistManagerDef.h>
0004 #include <qautils/QAUtil.h>
0005
0006 #include <g4detectors/PHG4TpcGeomContainer.h>
0007
0008 #include <g4main/PHG4HitContainer.h>
0009 #include <g4main/PHG4Particle.h>
0010 #include <g4main/PHG4TruthInfoContainer.h>
0011
0012 #include <trackbase_historic/ActsTransformations.h>
0013
0014 #include <trackbase/ActsGeometry.h>
0015 #include <trackbase/TpcDefs.h>
0016 #include <trackbase/TrackFitUtils.h>
0017 #include <trackbase/TrkrCluster.h>
0018 #include <trackbase/TrkrClusterContainer.h>
0019 #include <trackbase/TrkrClusterHitAssoc.h>
0020 #include <trackbase/TrkrDefs.h>  // for getTrkrId
0021 #include <trackbase/TrkrHitTruthAssoc.h>
0022
0023 #include <g4eval/SvtxClusterEval.h>  // for SvtxClusterEval
0024 #include <g4eval/SvtxEvalStack.h>
0025 #include <g4eval/SvtxTruthEval.h>  // for SvtxTruthEval
0026
0027 #include <fun4all/Fun4AllHistoManager.h>
0028 #include <fun4all/Fun4AllReturnCodes.h>
0029 #include <fun4all/SubsysReco.h>  // for SubsysReco
0030
0031 #include <phool/getClass.h>
0032 #include <phool/phool.h>  // for PHWHERE
0033
0034 #include <TAxis.h>  // for TAxis
0035 #include <TH1.h>
0036
0037 #include <cassert>
0038 #include <cmath>  // for atan2
0039 #include <format>
0040 #include <iostream>  // for operator<<, basic...
0041 #include <iterator>  // for distance
0042 #include <map>       // for map
0043 #include <set>
0044 #include <string>
0045 #include <utility>  // for pair, make_pair
0046 #include <vector>   // for vector
0047
0048 //________________________________________________________________________
0049 QAG4SimulationTpc::QAG4SimulationTpc(const std::string& name)
0050   : SubsysReco(name)
0051   , m_truthContainer(nullptr)
0052 {
0053 }
0054
0055 //________________________________________________________________________
0056 int QAG4SimulationTpc::InitRun(PHCompositeNode* topNode)
0057 {
0058   // prevent multiple creations of histograms
0059   if (m_initialized)
0060   {
0061     return Fun4AllReturnCodes::EVENT_OK;
0062   }
0063
0064   m_initialized = true;
0065
0066   if (!m_svtxEvalStack)
0067   {
0068     m_svtxEvalStack.reset(new SvtxEvalStack(topNode));
0069     m_svtxEvalStack->set_strict(true);
0070     m_svtxEvalStack->set_verbosity(Verbosity());
0071   }
0072
0073   m_truthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode,
0074                                                                 "G4TruthInfo");
0075   if (!m_truthContainer)
0076   {
0077     std::cout << "QAG4SimulationTpc::InitRun - Fatal Error - "
0078               << "unable to find node G4TruthInfo" << std::endl;
0079     assert(m_truthContainer);
0080   }
0081
0082   // find tpc geometry
0083   PHG4TpcGeomContainer* geom_container =
0084       findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0085   // auto geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "TPCGEOMCONTAINER");
0086   if (!geom_container)
0087   {
0088     std::cout << PHWHERE << " unable to find DST node TPCGEOMCONTAINER" << std::endl;
0089     return Fun4AllReturnCodes::ABORTRUN;
0090   }
0091
0092   // TPC has 3 regions, inner, mid and outer
0093   std::vector<int> region_layer_low = {7, 23, 39};
0094   ;
0095   std::vector<int> region_layer_high = {22, 38, 54};
0096
0097   // get layers from tpc geometry
0098   // make a layer to region multimap
0099   const auto range = geom_container->get_begin_end();
0100   for (auto iter = range.first; iter != range.second; ++iter)
0101   {
0102     m_layers.insert(iter->first);
0103
0104     for (int region = 0; region < 3; ++region)
0105     {
0106       if (iter->first >= region_layer_low[region] && iter->first <= region_layer_high[region])
0107       {
0108         m_layer_region_map.insert(std::make_pair(iter->first, region));
0109       }
0110     }
0111   }
0112
0113   // histogram manager
0114   auto* hm = QAHistManagerDef::getHistoManager();
0115   assert(hm);
0116
0117   // create histograms
0118
0119   // truth clusters
0120   {
0121     auto* h = new TH1F(std::format("{}efficiency_0", get_histo_prefix()).c_str(), "TPC_truth_clusters", 48, 7, 54);
0122     h->GetXaxis()->SetTitle("sPHENIX layer");
0123     hm->registerHisto(h);
0124   }
0125   // matched clusters
0126   {
0127     auto* h = new TH1F(std::format("{}efficiency_1", get_histo_prefix()).c_str(), "TPC_matched_clusters", 48, 7, 54);
0128     h->GetXaxis()->SetTitle("sPHENIX layer");
0129     hm->registerHisto(h);
0130   }
0131
0132   // cluster parameters
0133   for (int region = 0; region < 3; ++region)
0134   {
0135     if (Verbosity())
0136     {
0137       std::cout << PHWHERE << " adding region " << region << " with layers " << region_layer_low[region] << " to " << region_layer_high[region] << std::endl;
0138     }
0139     {
0140       // rphi residuals (cluster - truth)
0141       auto* h = new TH1F(std::format("{}drphi_{}", get_histo_prefix(), region).c_str(), std::format("TPC r#Delta#phi_{{cluster-truth}} region_{}", region).c_str(), 100, -0.079, 0.075);
0142       h->GetXaxis()->SetTitle("r#Delta#phi_{cluster-truth} (cm)");
0143       hm->registerHisto(h);
0144     }
0145
0146     {
0147       // rphi cluster errors
0148       auto* h = new TH1F(std::format("{}rphi_error_{}", get_histo_prefix(), region).c_str(), std::format("TPC r#Delta#phi error region_{}", region).c_str(), 100, 0, 0.075);
0149       h->GetXaxis()->SetTitle("r#Delta#phi error (cm)");
0150       hm->registerHisto(h);
0151     }
0152
0153     {
0154       // phi pulls (cluster - truth)
0155       auto* h = new TH1F(std::format("{}phi_pulls_{}", get_histo_prefix(), region).c_str(), std::format("TPC #Delta#phi_{{cluster-truth}}/#sigma#phi region_{}", region).c_str(), 100, -3, 3);
0156       h->GetXaxis()->SetTitle("#Delta#phi_{cluster-truth}/#sigma#phi (cm)");
0157       hm->registerHisto(h);
0158     }
0159
0160     {
0161       // z residuals (cluster - truth)
0162       auto* h = new TH1F(std::format("{}dz_{}", get_histo_prefix(), region).c_str(), std::format("TPC #Deltaz_{{cluster-truth}} region_{}", region).c_str(), 100, -0.19, 0.19);
0163       h->GetXaxis()->SetTitle("#Delta#z_{cluster-truth} (cm)");
0164       hm->registerHisto(h);
0165     }
0166
0167     {
0168       // z cluster errors
0169       auto* h = new TH1F(std::format("{}z_error_{}", get_histo_prefix(), region).c_str(), std::format("TPC z error region_{}", region).c_str(), 100, 0, 0.18);
0170       h->GetXaxis()->SetTitle("z error (cm)");
0171       hm->registerHisto(h);
0172     }
0173
0174     {
0175       // z pulls (cluster - truth)
0176       auto* h = new TH1F(std::format("{}z_pulls_{}", get_histo_prefix(), region).c_str(), std::format("TPC #Deltaz_{{cluster-truth}}/#sigmaz region_{}", region).c_str(), 100, -3, 3);
0177       h->GetXaxis()->SetTitle("#Delta#z_{cluster-truth}/#sigmaz (cm)");
0178       hm->registerHisto(h);
0179     }
0180
0181     {
0182       // total cluster size
0183       auto* h = new TH1F(std::format("{}clus_size_{}", get_histo_prefix(), region).c_str(), std::format("TPC cluster size region_{}", region).c_str(), 30, 0, 30);
0184       h->GetXaxis()->SetTitle("csize");
0185       hm->registerHisto(h);
0186     }
0187
0188     {
0189       // cluster size in phi
0190       auto* h = new TH1F(std::format("{}clus_size_phi_{}", get_histo_prefix(), region).c_str(), std::format("TPC cluster size (#phi) region_{}", region).c_str(), 10, 0, 10);
0191       h->GetXaxis()->SetTitle("csize_{#phi}");
0192       hm->registerHisto(h);
0193     }
0194
0195     {
0196       // cluster size in z
0197       auto* h = new TH1F(std::format("{}clus_size_z_{}", get_histo_prefix(), region).c_str(), std::format("TPC cluster size (z) region_{}", region).c_str(), 12, 0, 12);
0198       h->GetXaxis()->SetTitle("csize_{z}");
0199       hm->registerHisto(h);
0200     }
0201   }
0202
0203   return Fun4AllReturnCodes::EVENT_OK;
0204 }
0205
0206 //_____________________________________________________________________
0207 int QAG4SimulationTpc::process_event(PHCompositeNode* topNode)
0208 {
0209   // load nodes
0210   auto res = load_nodes(topNode);
0211   if (res != Fun4AllReturnCodes::EVENT_OK)
0212   {
0213     return res;
0214   }
0215
0216   if (m_svtxEvalStack)
0217   {
0218     m_svtxEvalStack->next_event(topNode);
0219   }
0220
0221   // run evaluation
0222   evaluate_clusters();
0223   return Fun4AllReturnCodes::EVENT_OK;
0224 }
0225
0226 //________________________________________________________________________
0227 std::string QAG4SimulationTpc::get_histo_prefix() const
0228 {
0229   return std::string("h_") + Name() + std::string("_");
0230 }
0231
0232 //________________________________________________________________________
0233 int QAG4SimulationTpc::load_nodes(PHCompositeNode* topNode)
0234 {
0235   m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0236   if (!m_cluster_map)
0237   {
0238     std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTER" << std::endl;
0239     return Fun4AllReturnCodes::ABORTEVENT;
0240   }
0241
0242   m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0243   if (!m_tGeometry)
0244   {
0245     std::cout << PHWHERE << "No acts tracking geometry, exiting."
0246               << std::endl;
0247     return Fun4AllReturnCodes::ABORTEVENT;
0248   }
0249
0250   m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
0251   if (!m_cluster_hit_map)
0252   {
0253     std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTERHITASSOC" << std::endl;
0254     return Fun4AllReturnCodes::ABORTEVENT;
0255   }
0256
0257   m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
0258   if (!m_hit_truth_map)
0259   {
0260     std::cout << PHWHERE << " unable to find DST node TRKR_HITTRUTHASSOC" << std::endl;
0261     return Fun4AllReturnCodes::ABORTEVENT;
0262   }
0263
0264   m_g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
0265   if (!m_g4hits_tpc)
0266   {
0267     std::cout << PHWHERE << " unable to find DST node G4HIT_TPC" << std::endl;
0268     return Fun4AllReturnCodes::ABORTEVENT;
0269   }
0270
0271   return Fun4AllReturnCodes::EVENT_OK;
0272 }
0273
0274 //________________________________________________________________________
0275 void QAG4SimulationTpc::evaluate_clusters()
0276 {
0277   SvtxTruthEval* trutheval = m_svtxEvalStack->get_truth_eval();
0278   assert(trutheval);
0279   SvtxClusterEval* clustereval = m_svtxEvalStack->get_cluster_eval();
0280   assert(clustereval);
0281
0282   // histogram manager
0283   auto* hm = QAHistManagerDef::getHistoManager();
0284   assert(hm);
0285
0286   // get histograms for cluster efficiency
0287   TH1* h_eff0 = dynamic_cast<TH1*>(hm->getHisto(std::format("{}efficiency_0", get_histo_prefix())));
0288   assert(h_eff0);
0289   TH1* h_eff1 = dynamic_cast<TH1*>(hm->getHisto(std::format("{}efficiency_1", get_histo_prefix())));
0290   assert(h_eff1);
0291
0292   // get histograms for cluster parameters vs truth
0293   struct HistogramList
0294   {
0295     TH1* drphi = nullptr;
0296     TH1* rphi_error = nullptr;
0297     TH1* phi_pulls = nullptr;
0298
0299     TH1* dz = nullptr;
0300     TH1* z_error = nullptr;
0301     TH1* z_pulls = nullptr;
0302
0303     TH1* csize = nullptr;
0304     TH1* csize_phi = nullptr;
0305     TH1* csize_z = nullptr;
0306   };
0307
0308   using HistogramMap = std::map<int, HistogramList>;
0309   HistogramMap histograms;
0310
0311   for (int region = 0; region < 3; ++region)
0312   {
0313     HistogramList h;
0314     h.drphi = dynamic_cast<TH1*>(hm->getHisto(std::format("{}drphi_{}", get_histo_prefix(), region)));
0315     h.rphi_error = dynamic_cast<TH1*>(hm->getHisto(std::format("{}rphi_error_{}", get_histo_prefix(), region)));
0316     h.phi_pulls = dynamic_cast<TH1*>(hm->getHisto(std::format("{}phi_pulls_{}", get_histo_prefix(), region)));
0317
0318     h.dz = dynamic_cast<TH1*>(hm->getHisto(std::format("{}dz_{}", get_histo_prefix(), region)));
0319     h.z_error = dynamic_cast<TH1*>(hm->getHisto(std::format("{}z_error_{}", get_histo_prefix(), region)));
0320     h.z_pulls = dynamic_cast<TH1*>(hm->getHisto(std::format("{}z_pulls_{}", get_histo_prefix(), region)));
0321
0322     h.csize = dynamic_cast<TH1*>(hm->getHisto(std::format("{}clus_size_{}", get_histo_prefix(), region)));
0323     h.csize_phi = dynamic_cast<TH1*>(hm->getHisto(std::format("{}clus_size_phi_{}", get_histo_prefix(), region)));
0324     h.csize_z = dynamic_cast<TH1*>(hm->getHisto(std::format("{}clus_size_z_{}", get_histo_prefix(), region)));
0325
0326     histograms.insert(std::make_pair(region, h));
0327   }
0328
0329   // Get all truth clusters
0330   //===============
0331   if (Verbosity() > 0)
0332   {
0333     std::cout << PHWHERE << " get all truth clusters for primary particles " << std::endl;
0334   }
0335
0336   // PHG4TruthInfoContainer::ConstRange range = m_truthContainer->GetParticleRange();  // all truth cluters
0337   PHG4TruthInfoContainer::ConstRange const range = m_truthContainer->GetPrimaryParticleRange();  // only from primary particles
0338
0339   for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0340        iter != range.second;
0341        ++iter)
0342   {
0343     PHG4Particle* g4particle = iter->second;
0344
0345     float const gtrackID = g4particle->get_track_id();
0346     float const gflavor = g4particle->get_pid();
0347     float const gembed = trutheval->get_embed(g4particle);
0348     float const gprimary = trutheval->is_primary(g4particle);
0349
0350     if (Verbosity() > 0)
0351     {
0352       std::cout << PHWHERE << " PHG4Particle ID " << gtrackID << " gembed " << gembed << " gflavor " << gflavor << " gprimary " << gprimary << std::endl;
0353     }
0354
0355     // Get the truth clusters from this particle
0356     const auto truth_clusters = trutheval->all_truth_clusters(g4particle);
0357
0358     // get circle fit parameters first
0359     TrackFitUtils::position_vector_t xy_pts;
0360     TrackFitUtils::position_vector_t rz_pts;
0361
0362     for (const auto& [gkey, gclus] : truth_clusters)
0363     {
0364       const auto layer = TrkrDefs::getLayer(gkey);
0365       if (layer < 7)
0366       {
0367         continue;
0368       }
0369
0370       float const gx = gclus->getX();
0371       float const gy = gclus->getY();
0372       float const gz = gclus->getZ();
0373
0374       xy_pts.emplace_back(gx, gy);
0375       rz_pts.emplace_back(std::sqrt(gx * gx + gy * gy), gz);
0376     }
0377
0378     // fit a circle through x,y coordinates
0379     const auto [R, X0, Y0] = TrackFitUtils::circle_fit_by_taubin(xy_pts);
0380
0381     // skip chain entirely if fit fails
0382     if (std::isnan(R))
0383     {
0384       continue;
0385     }
0386
0387     // process residuals and pulls
0388     for (const auto& [gkey, gclus] : truth_clusters)
0389     {
0390       const auto layer = TrkrDefs::getLayer(gkey);
0391       const auto detID = TrkrDefs::getTrkrId(gkey);
0392       // if (detID != TrkrDefs::tpcId) continue;
0393       if (layer < 7)
0394       {
0395         continue;
0396       }
0397
0398       float const gx = gclus->getX();
0399       float const gy = gclus->getY();
0400       float const gz = gclus->getZ();
0401       float const gedep = gclus->getError(0, 0);
0402       float const ng4hits = gclus->getAdc();
0403
0404       const auto gr = QAG4Util::get_r(gclus->getX(), gclus->getY());
0405       const auto gphi = std::atan2(gclus->getY(), gclus->getX());
0406
0407       if (Verbosity() > 0)
0408       {
0409         std::cout << "     gkey " << gkey << " detID " << detID << " tpcId " << TrkrDefs::tpcId << " layer " << layer << "  truth clus " << gkey << " ng4hits " << ng4hits << " gr " << gr << " gx " << gx << " gy " << gy << " gz " << gz
0410                   << " gphi " << gphi << " gedep " << gedep << std::endl;
0411       }
0412
0413       // fill the truth cluster histo
0414       h_eff0->Fill(layer);
0415
0416       // find matching reco cluster histo
0417       const auto [rkey, rclus] = clustereval->reco_cluster_from_truth_cluster(gkey, gclus);
0418       if (rclus)
0419       {
0420         // fill the matched cluster histo
0421         h_eff1->Fill(layer);
0422
0423         const auto global = m_tGeometry->getGlobalPosition(rkey, rclus);
0424
0425         // get relevant cluster information
0426         const auto r_cluster = QAG4Util::get_r(global(0), global(1));
0427         const auto z_cluster = global(2);
0428         const auto phi_cluster = (float) std::atan2(global(1), global(0));
0429
0430         double const phi_error = rclus->getRPhiError() / r_cluster;
0431         double const z_error = rclus->getZError();
0432
0433         const auto dphi = QAG4Util::delta_phi(phi_cluster, gphi);
0434         const auto dz = z_cluster - gz;
0435
0436         // get region from layer, fill histograms
0437         const auto it = m_layer_region_map.find(layer);
0438         int const region = it->second;
0439
0440         if (Verbosity() > 0)
0441         {
0442           std::cout << "   Found match in layer " << layer << " region " << region << " for gtrackID " << gtrackID << std::endl;
0443           std::cout << "      x " << rclus->getX() << " y " << rclus->getY() << " z " << rclus->getZ() << std::endl;
0444           std::cout << "     gx " << gclus->getX() << " gy " << gclus->getY() << " gz " << gclus->getZ() << std::endl;
0445           std::cout << "     drphi " << r_cluster * dphi << " rphi_error " << r_cluster * phi_error << " dz " << dz << " z_error " << z_error << std::endl;
0446         }
0447
0448         const auto hiter = histograms.find(region);
0449         if (hiter == histograms.end())
0450         {
0451           continue;
0452         }
0453
0454         // fill phi residuals, errors and pulls
0455         auto fill = [](TH1* h, float value)
0456         { if( h ) { h->Fill( value );
0457 } };
0458         fill(hiter->second.drphi, r_cluster * dphi);
0459         fill(hiter->second.rphi_error, r_cluster * phi_error);
0460         fill(hiter->second.phi_pulls, dphi / phi_error);
0461
0462         // fill z residuals, errors and pulls
0463         fill(hiter->second.dz, dz);
0464         fill(hiter->second.z_error, z_error);
0465         fill(hiter->second.z_pulls, dz / z_error);
0466
0467         // cluster sizes
0468         // get associated hits
0469         const auto hit_range = m_cluster_hit_map->getHits(rkey);
0470         fill(hiter->second.csize, std::distance(hit_range.first, hit_range.second));
0471
0472         std::set<int> phibins;
0473         std::set<int> zbins;
0474         for (auto hit_iter = hit_range.first; hit_iter != hit_range.second; ++hit_iter)
0475         {
0476           const auto& hit_key = hit_iter->second;
0477           phibins.insert(TpcDefs::getPad(hit_key));
0478           zbins.insert(TpcDefs::getTBin(hit_key));
0479         }
0480
0481         fill(hiter->second.csize_phi, phibins.size());
0482         fill(hiter->second.csize_z, zbins.size());
0483       }
0484     }
0485   }
0486 }
0487
0488 //_____________________________________________________________________
0489 QAG4SimulationTpc::G4HitSet QAG4SimulationTpc::find_g4hits(TrkrDefs::cluskey cluster_key) const
0490 {
0491   // find hitset associated to cluster
0492   G4HitSet out;
0493   const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
0494
0495   // loop over hits associated to clusters
0496   const auto range = m_cluster_hit_map->getHits(cluster_key);
0497   for (auto iter = range.first; iter != range.second; ++iter)
0498   {
0499     // hit key
0500     const auto& hit_key = iter->second;
0501
0502     // store hits to g4hit associations
0503     TrkrHitTruthAssoc::MMap g4hit_map;
0504     m_hit_truth_map->getG4Hits(hitset_key, hit_key, g4hit_map);
0505
0506     // find corresponding g4 hist
0507     for (auto& truth_iter : g4hit_map)
0508     {
0509       // g4hit key
0510       const auto g4hit_key = truth_iter.second.second;
0511
0512       // g4 hit
0513       PHG4Hit* g4hit = (TrkrDefs::getTrkrId(hitset_key) == TrkrDefs::tpcId) ? m_g4hits_tpc->findHit(g4hit_key) : nullptr;
0514
0515       // insert in set
0516       if (g4hit)
0517       {
0518         out.insert(g4hit);
0519       }
0520     }
0521   }
0522
0523   return out;
0524 }