QA/SimulationModules/QAG4SimulationTpc.cc

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