ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4SimulationTpc.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4SimulationTpc.cc
1 #include "QAG4SimulationTpc.h"
2 #include "QAG4Util.h"
3 #include "QAHistManagerDef.h"
4 
6 
7 #include <g4main/PHG4Particle.h>
9 
11 
12 #include <tpc/TpcDefs.h>
13 
15 
18 #include <trackbase/TrkrCluster.h>
21 #include <trackbase/TrkrDefs.h> // for getTrkrId
23 
24 #include <g4eval/SvtxClusterEval.h> // for SvtxClusterEval
25 #include <g4eval/SvtxEvalStack.h>
26 #include <g4eval/SvtxTruthEval.h> // for SvtxTruthEval
27 
30 #include <fun4all/SubsysReco.h> // for SubsysReco
31 
32 #include <phool/getClass.h>
33 #include <phool/phool.h> // for PHWHERE
34 
35 #include <TAxis.h> // for TAxis
36 #include <TH1.h>
37 #include <TString.h> // for Form
38 
39 #include <cassert>
40 #include <cmath> // for atan2
41 #include <iostream> // for operator<<, basic...
42 #include <iterator> // for distance
43 #include <map> // for map
44 #include <utility> // for pair, make_pair
45 #include <vector> // for vector
46 
47 //________________________________________________________________________
49  : SubsysReco(name)
50  , m_truthContainer(nullptr)
51 {
52 }
53 
54 //________________________________________________________________________
56 {
57  // prevent multiple creations of histograms
58  if (m_initialized)
60  else
61  m_initialized = true;
62 
63  if (!m_svtxEvalStack)
64  {
65  m_svtxEvalStack.reset(new SvtxEvalStack(topNode));
66  m_svtxEvalStack->set_strict(true);
67  m_svtxEvalStack->set_verbosity(Verbosity());
68  }
69 
70  m_truthContainer = findNode::getClass<PHG4TruthInfoContainer>(topNode,
71  "G4TruthInfo");
72  if (!m_truthContainer)
73  {
74  std::cout << "QAG4SimulationTpc::InitRun - Fatal Error - "
75  << "unable to find node G4TruthInfo" << std::endl;
76  assert(m_truthContainer);
77  }
78 
79  // find tpc geometry
80  PHG4CylinderCellGeomContainer* geom_container =
81  findNode::getClass<PHG4CylinderCellGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
82  //auto geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
83  if (!geom_container)
84  {
85  std::cout << PHWHERE << " unable to find DST node CYLINDERCELLGEOM_SVTX" << std::endl;
87  }
88 
89  // TPC has 3 regions, inner, mid and outer
90  std::vector<int> region_layer_low = {7, 23, 39};
91  ;
92  std::vector<int> region_layer_high = {22, 38, 54};
93 
94  // get layers from tpc geometry
95  // make a layer to region multimap
96  const auto range = geom_container->get_begin_end();
97  for (auto iter = range.first; iter != range.second; ++iter)
98  {
99  m_layers.insert(iter->first);
100 
101  for (int region = 0; region < 3; ++region)
102  {
103  if (iter->first >= region_layer_low[region] && iter->first <= region_layer_high[region])
104  m_layer_region_map.insert(std::make_pair(iter->first, region));
105  }
106  }
107 
108  // histogram manager
110  assert(hm);
111 
112  // create histograms
113 
114  // truth clusters
115  {
116  auto h = new TH1F(Form("%sefficiency_0", get_histo_prefix().c_str()), Form("TPC_truth_clusters"), 48, 7, 54);
117  h->GetXaxis()->SetTitle("sPHENIX layer");
118  hm->registerHisto(h);
119  }
120  // matched clusters
121  {
122  auto h = new TH1F(Form("%sefficiency_1", get_histo_prefix().c_str()), Form("TPC_matched_clusters"), 48, 7, 54);
123  h->GetXaxis()->SetTitle("sPHENIX layer");
124  hm->registerHisto(h);
125  }
126 
127  // cluster parameters
128  for (int region = 0; region < 3; ++region)
129  {
130  if (Verbosity()) std::cout << PHWHERE << " adding region " << region << " with layers " << region_layer_low[region] << " to " << region_layer_high[region] << std::endl;
131  {
132  // rphi residuals (cluster - truth)
133  auto h = new TH1F(Form("%sdrphi_%i", get_histo_prefix().c_str(), region), Form("TPC r#Delta#phi_{cluster-truth} region_%i", region), 100, -0.079, 0.075);
134  h->GetXaxis()->SetTitle("r#Delta#phi_{cluster-truth} (cm)");
135  hm->registerHisto(h);
136  }
137 
138  {
139  // rphi cluster errors
140  auto h = new TH1F(Form("%srphi_error_%i", get_histo_prefix().c_str(), region), Form("TPC r#Delta#phi error region_%i", region), 100, 0, 0.075);
141  h->GetXaxis()->SetTitle("r#Delta#phi error (cm)");
142  hm->registerHisto(h);
143  }
144 
145  {
146  // phi pulls (cluster - truth)
147  auto h = new TH1F(Form("%sphi_pulls_%i", get_histo_prefix().c_str(), region), Form("TPC #Delta#phi_{cluster-truth}/#sigma#phi region_%i", region), 100, -3, 3);
148  h->GetXaxis()->SetTitle("#Delta#phi_{cluster-truth}/#sigma#phi (cm)");
149  hm->registerHisto(h);
150  }
151 
152  {
153  // z residuals (cluster - truth)
154  auto h = new TH1F(Form("%sdz_%i", get_histo_prefix().c_str(), region), Form("TPC #Deltaz_{cluster-truth} region_%i", region), 100, -0.19, 0.19);
155  h->GetXaxis()->SetTitle("#Delta#z_{cluster-truth} (cm)");
156  hm->registerHisto(h);
157  }
158 
159  {
160  // z cluster errors
161  auto h = new TH1F(Form("%sz_error_%i", get_histo_prefix().c_str(), region), Form("TPC z error region_%i", region), 100, 0, 0.18);
162  h->GetXaxis()->SetTitle("z error (cm)");
163  hm->registerHisto(h);
164  }
165 
166  {
167  // z pulls (cluster - truth)
168  auto h = new TH1F(Form("%sz_pulls_%i", get_histo_prefix().c_str(), region), Form("TPC #Deltaz_{cluster-truth}/#sigmaz region_%i", region), 100, -3, 3);
169  h->GetXaxis()->SetTitle("#Delta#z_{cluster-truth}/#sigmaz (cm)");
170  hm->registerHisto(h);
171  }
172 
173  {
174  // total cluster size
175  auto h = new TH1F(Form("%sclus_size_%i", get_histo_prefix().c_str(), region), Form("TPC cluster size region_%i", region), 30, 0, 30);
176  h->GetXaxis()->SetTitle("csize");
177  hm->registerHisto(h);
178  }
179 
180  {
181  // cluster size in phi
182  auto h = new TH1F(Form("%sclus_size_phi_%i", get_histo_prefix().c_str(), region), Form("TPC cluster size (#phi) region_%i", region), 10, 0, 10);
183  h->GetXaxis()->SetTitle("csize_{#phi}");
184  hm->registerHisto(h);
185  }
186 
187  {
188  // cluster size in z
189  auto h = new TH1F(Form("%sclus_size_z_%i", get_histo_prefix().c_str(), region), Form("TPC cluster size (z) region_%i", region), 12, 0, 12);
190  h->GetXaxis()->SetTitle("csize_{z}");
191  hm->registerHisto(h);
192  }
193  }
194 
196 }
197 
198 //_____________________________________________________________________
200 {
201  // load nodes
202  auto res = load_nodes(topNode);
203  if (res != Fun4AllReturnCodes::EVENT_OK) return res;
204 
205  if (m_svtxEvalStack)
206  m_svtxEvalStack->next_event(topNode);
207 
208  // run evaluation
211 }
212 
213 //________________________________________________________________________
215 {
216  return std::string("h_") + Name() + std::string("_");
217 }
218 
219 //________________________________________________________________________
221 {
222  m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
223  if (!m_cluster_map)
224  {
225  std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTER" << std::endl;
227  }
228 
229  m_surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
230  if (!m_surfmaps)
231  {
232  std::cout << PHWHERE << "Error: can't find Acts surface maps"
233  << std::endl;
235  }
236 
237  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
238  if (!m_tGeometry)
239  {
240  std::cout << PHWHERE << "No acts tracking geometry, exiting."
241  << std::endl;
243  }
244 
245  m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
246  if (!m_cluster_hit_map)
247  {
248  std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTERHITASSOC" << std::endl;
250  }
251 
252  m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
253  if (!m_hit_truth_map)
254  {
255  std::cout << PHWHERE << " unable to find DST node TRKR_HITTRUTHASSOC" << std::endl;
257  }
258 
259  m_g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_TPC");
260  if (!m_g4hits_tpc)
261  {
262  std::cout << PHWHERE << " unable to find DST node G4HIT_TPC" << std::endl;
264  }
265 
267 }
268 
269 //________________________________________________________________________
271 {
272  SvtxTruthEval* trutheval = m_svtxEvalStack->get_truth_eval();
273  assert(trutheval);
274  SvtxClusterEval* clustereval = m_svtxEvalStack->get_cluster_eval();
275  assert(clustereval);
276 
277  // histogram manager
279  assert(hm);
280 
281  // get histograms for cluster efficiency
282  TH1* h_eff0 = dynamic_cast<TH1*>(hm->getHisto(Form("%sefficiency_0", get_histo_prefix().c_str())));
283  assert(h_eff0);
284  TH1* h_eff1 = dynamic_cast<TH1*>(hm->getHisto(Form("%sefficiency_1", get_histo_prefix().c_str())));
285  assert(h_eff1);
286 
287  // get histograms for cluster parameters vs truth
288  struct HistogramList
289  {
290  TH1* drphi = nullptr;
291  TH1* rphi_error = nullptr;
292  TH1* phi_pulls = nullptr;
293 
294  TH1* dz = nullptr;
295  TH1* z_error = nullptr;
296  TH1* z_pulls = nullptr;
297 
298  TH1* csize = nullptr;
299  TH1* csize_phi = nullptr;
300  TH1* csize_z = nullptr;
301  };
302 
303  using HistogramMap = std::map<int, HistogramList>;
304  HistogramMap histograms;
305 
306  for (int region = 0; region < 3; ++region)
307  {
308  HistogramList h;
309  h.drphi = dynamic_cast<TH1*>(hm->getHisto(Form("%sdrphi_%i", get_histo_prefix().c_str(), region)));
310  h.rphi_error = dynamic_cast<TH1*>(hm->getHisto(Form("%srphi_error_%i", get_histo_prefix().c_str(), region)));
311  h.phi_pulls = dynamic_cast<TH1*>(hm->getHisto(Form("%sphi_pulls_%i", get_histo_prefix().c_str(), region)));
312 
313  h.dz = dynamic_cast<TH1*>(hm->getHisto(Form("%sdz_%i", get_histo_prefix().c_str(), region)));
314  h.z_error = dynamic_cast<TH1*>(hm->getHisto(Form("%sz_error_%i", get_histo_prefix().c_str(), region)));
315  h.z_pulls = dynamic_cast<TH1*>(hm->getHisto(Form("%sz_pulls_%i", get_histo_prefix().c_str(), region)));
316 
317  h.csize = dynamic_cast<TH1*>(hm->getHisto(Form("%sclus_size_%i", get_histo_prefix().c_str(), region)));
318  h.csize_phi = dynamic_cast<TH1*>(hm->getHisto(Form("%sclus_size_phi_%i", get_histo_prefix().c_str(), region)));
319  h.csize_z = dynamic_cast<TH1*>(hm->getHisto(Form("%sclus_size_z_%i", get_histo_prefix().c_str(), region)));
320 
321  histograms.insert(std::make_pair(region, h));
322  }
323 
324  // Get all truth clusters
325  //===============
326  if (Verbosity() > 0)
327  std::cout << PHWHERE << " get all truth clusters for primary particles " << std::endl;
328 
329  //PHG4TruthInfoContainer::ConstRange range = m_truthContainer->GetParticleRange(); // all truth cluters
330  PHG4TruthInfoContainer::ConstRange range = m_truthContainer->GetPrimaryParticleRange(); // only from primary particles
331 
332  ActsTransformations transformer;
333 
334  for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
335  iter != range.second;
336  ++iter)
337  {
338  PHG4Particle* g4particle = iter->second;
339 
340  float gtrackID = g4particle->get_track_id();
341  float gflavor = g4particle->get_pid();
342  float gembed = trutheval->get_embed(g4particle);
343  float gprimary = trutheval->is_primary(g4particle);
344 
345  if (Verbosity() > 0)
346  std::cout << PHWHERE << " PHG4Particle ID " << gtrackID << " gembed " << gembed << " gflavor " << gflavor << " gprimary " << gprimary << std::endl;
347 
348  // Get the truth clusters from this particle
349  std::map<unsigned int, std::shared_ptr<TrkrCluster> > truth_clusters = trutheval->all_truth_clusters(g4particle);
350  for (auto it = truth_clusters.begin(); it != truth_clusters.end(); ++it)
351  {
352  unsigned int layer = it->first;
353  std::shared_ptr<TrkrCluster> gclus = it->second;
354  const auto gkey = gclus->getClusKey();
355  const auto detID = TrkrDefs::getTrkrId(gkey);
356  //if (detID != TrkrDefs::tpcId) continue;
357  if (layer < 7) continue;
358 
359  float gx = gclus->getX();
360  float gy = gclus->getY();
361  float gz = gclus->getZ();
362  float gedep = gclus->getError(0, 0);
363  float ng4hits = gclus->getAdc();
364 
365  const auto gr = QAG4Util::get_r(gclus->getX(), gclus->getY());
366  const auto gphi = std::atan2(gclus->getY(), gclus->getX());
367 
368  if (Verbosity() > 0)
369  {
370  std::cout << " gkey " << gkey << " detID " << detID << " tpcId " << TrkrDefs::tpcId << " layer " << layer << " truth clus " << gkey << " ng4hits " << ng4hits << " gr " << gr << " gx " << gx << " gy " << gy << " gz " << gz
371  << " gphi " << gphi << " gedep " << gedep << std::endl;
372  }
373 
374  // fill the truth cluster histo
375  h_eff0->Fill(layer);
376 
377  // find matching reco cluster histo
378  TrkrCluster* rclus = clustereval->reco_cluster_from_truth_cluster(gclus);
379  if (rclus)
380  {
381  // fill the matched cluster histo
382  h_eff1->Fill(layer);
383 
384  const auto global = transformer.getGlobalPosition(rclus, m_surfmaps,
385  m_tGeometry);
386 
387  // get relevant cluster information
388  const auto rkey = rclus->getClusKey();
389  const auto r_cluster = QAG4Util::get_r(global(0), global(1));
390  const auto z_cluster = global(2);
391  const auto phi_cluster = (float) std::atan2(global(1), global(0));
392  const auto phi_error = rclus->getRPhiError() / r_cluster;
393  const auto z_error = rclus->getZError();
394 
395  const auto dphi = QAG4Util::delta_phi(phi_cluster, gphi);
396  const auto dz = z_cluster - gz;
397 
398  // get region from layer, fill histograms
399  const auto it = m_layer_region_map.find(layer);
400  int region = it->second;
401 
402  if (Verbosity() > 0)
403  {
404  std::cout << " Found match in layer " << layer << " region " << region << " for gtrackID " << gtrackID << std::endl;
405  std::cout << " x " << rclus->getX() << " y " << rclus->getY() << " z " << rclus->getZ() << std::endl;
406  std::cout << " gx " << gclus->getX() << " gy " << gclus->getY() << " gz " << gclus->getZ() << std::endl;
407  std::cout << " drphi " << r_cluster * dphi << " rphi_error " << r_cluster * phi_error << " dz " << dz << " z_error " << z_error << std::endl;
408  }
409 
410  const auto hiter = histograms.find(region);
411  if (hiter == histograms.end()) continue;
412 
413  // fill phi residuals, errors and pulls
414  auto fill = [](TH1* h, float value) { if( h ) h->Fill( value ); };
415  fill(hiter->second.drphi, r_cluster * dphi);
416  fill(hiter->second.rphi_error, r_cluster * phi_error);
417  fill(hiter->second.phi_pulls, dphi / phi_error);
418 
419  // fill z residuals, errors and pulls
420  fill(hiter->second.dz, dz);
421  fill(hiter->second.z_error, z_error);
422  fill(hiter->second.z_pulls, dz / z_error);
423 
424  // cluster sizes
425  // get associated hits
426  const auto hit_range = m_cluster_hit_map->getHits(rkey);
427  fill(hiter->second.csize, std::distance(hit_range.first, hit_range.second));
428 
429  std::set<int> phibins;
430  std::set<int> zbins;
431  for (auto hit_iter = hit_range.first; hit_iter != hit_range.second; ++hit_iter)
432  {
433  const auto& hit_key = hit_iter->second;
434  phibins.insert(TpcDefs::getPad(hit_key));
435  zbins.insert(TpcDefs::getTBin(hit_key));
436  }
437 
438  fill(hiter->second.csize_phi, phibins.size());
439  fill(hiter->second.csize_z, zbins.size());
440  }
441  }
442  }
443 }
444 
445 //_____________________________________________________________________
447 {
448  // find hitset associated to cluster
449  G4HitSet out;
450  const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
451 
452  // loop over hits associated to clusters
453  const auto range = m_cluster_hit_map->getHits(cluster_key);
454  for (auto iter = range.first; iter != range.second; ++iter)
455  {
456  // hit key
457  const auto& hit_key = iter->second;
458 
459  // store hits to g4hit associations
460  TrkrHitTruthAssoc::MMap g4hit_map;
461  m_hit_truth_map->getG4Hits(hitset_key, hit_key, g4hit_map);
462 
463  // find corresponding g4 hist
464  for (auto truth_iter = g4hit_map.begin(); truth_iter != g4hit_map.end(); ++truth_iter)
465  {
466  // g4hit key
467  const auto g4hit_key = truth_iter->second.second;
468 
469  // g4 hit
470  PHG4Hit* g4hit = (TrkrDefs::getTrkrId(hitset_key) == TrkrDefs::tpcId) ? m_g4hits_tpc->findHit(g4hit_key) : nullptr;
471 
472  // insert in set
473  if (g4hit) out.insert(g4hit);
474  }
475  }
476 
477  return out;
478 }