ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4SimulationMvtx.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4SimulationMvtx.cc
1 #include "QAG4SimulationMvtx.h"
2 #include "QAG4Util.h"
3 #include "QAHistManagerDef.h"
4 
6 
7 #include <g4main/PHG4Hit.h>
9 
10 #include <mvtx/MvtxDefs.h>
11 
14 #include <trackbase/TrkrCluster.h>
17 #include <trackbase/TrkrDefs.h> // for getTrkrId, getHit...
20 
22 
25 #include <fun4all/SubsysReco.h> // for SubsysReco
26 
27 #include <phool/getClass.h>
28 #include <phool/phool.h> // for PHWHERE
29 
30 #include <TAxis.h> // for TAxis
31 #include <TH1.h>
32 #include <TString.h> // for Form
33 
34 #include <cassert>
35 #include <cmath> // for atan2
36 #include <iostream> // for operator<<, basic...
37 #include <iterator> // for distance
38 #include <map> // for map
39 #include <utility> // for pair, make_pair
40 
41 //________________________________________________________________________
43  : SubsysReco(name)
44 {
45 }
46 
47 //________________________________________________________________________
49 {
50  // prevent multiple creations of histograms
51  if (m_initialized)
53  else
54  m_initialized = true;
55 
56  // find mvtx geometry
57  auto geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MVTX");
58  if (!geom_container)
59  {
60  std::cout << PHWHERE << " unable to find DST node CYLINDERGEOM_MVTX" << std::endl;
62  }
63 
64  // get layers from mvtx geometry
65  const auto range = geom_container->get_begin_end();
66  for (auto iter = range.first; iter != range.second; ++iter)
67  {
68  m_layers.insert(iter->first);
69  }
70 
71  // histogram manager
73  assert(hm);
74 
75  // create histograms
76  for (const auto& layer : m_layers)
77  {
78  if (Verbosity()) std::cout << PHWHERE << " adding layer " << layer << std::endl;
79  {
80  // rphi residuals (cluster - truth)
81  auto h = new TH1F(Form("%sdrphi_%i", get_histo_prefix().c_str(), layer), Form("MVTX r#Delta#phi_{cluster-truth} layer_%i", layer), 100, -2e-3, 2e-3);
82  h->GetXaxis()->SetTitle("r#Delta#phi_{cluster-truth} (cm)");
83  hm->registerHisto(h);
84  }
85 
86  {
87  // rphi cluster errors
88  auto h = new TH1F(Form("%srphi_error_%i", get_histo_prefix().c_str(), layer), Form("MVTX r#Delta#phi error layer_%i", layer), 100, 0, 2e-3);
89  h->GetXaxis()->SetTitle("r#Delta#phi error (cm)");
90  hm->registerHisto(h);
91  }
92 
93  {
94  // phi pulls (cluster - truth)
95  auto h = new TH1F(Form("%sphi_pulls_%i", get_histo_prefix().c_str(), layer), Form("MVTX #Delta#phi_{cluster-truth}/#sigma#phi layer_%i", layer), 100, -3, 3);
96  h->GetXaxis()->SetTitle("#Delta#phi_{cluster-truth}/#sigma#phi");
97  hm->registerHisto(h);
98  }
99 
100  {
101  // z residuals (cluster - truth)
102  auto h = new TH1F(Form("%sdz_%i", get_histo_prefix().c_str(), layer), Form("MVTX #Deltaz_{cluster-truth} layer_%i", layer), 100, -3e-3, 3e-3);
103  h->GetXaxis()->SetTitle("#Deltaz_{cluster-truth} (cm)");
104  hm->registerHisto(h);
105  }
106 
107  {
108  // z cluster errors
109  auto h = new TH1F(Form("%sz_error_%i", get_histo_prefix().c_str(), layer), Form("MVTX z error layer_%i", layer), 100, 0, 3e-3);
110  h->GetXaxis()->SetTitle("z error (cm)");
111  hm->registerHisto(h);
112  }
113 
114  {
115  // z pulls (cluster - truth)
116  auto h = new TH1F(Form("%sz_pulls_%i", get_histo_prefix().c_str(), layer), Form("MVTX #Deltaz_{cluster-truth}/#sigmaz layer_%i", layer), 100, -3, 3);
117  h->GetXaxis()->SetTitle("#Deltaz_{cluster-truth}/#sigmaz");
118  hm->registerHisto(h);
119  }
120 
121  {
122  // total cluster size
123  auto h = new TH1F(Form("%sclus_size_%i", get_histo_prefix().c_str(), layer), Form("MVTX cluster size layer_%i", layer), 20, 0, 20);
124  h->GetXaxis()->SetTitle("csize");
125  hm->registerHisto(h);
126  }
127 
128  {
129  // cluster size in phi
130  auto h = new TH1F(Form("%sclus_size_phi_%i", get_histo_prefix().c_str(), layer), Form("MVTX cluster size (#phi) layer_%i", layer), 10, 0, 10);
131  h->GetXaxis()->SetTitle("csize_{#phi}");
132  hm->registerHisto(h);
133  }
134 
135  {
136  // cluster size in z
137  auto h = new TH1F(Form("%sclus_size_z_%i", get_histo_prefix().c_str(), layer), Form("MVTX cluster size (z) layer_%i", layer), 10, 0, 10);
138  h->GetXaxis()->SetTitle("csize_{z}");
139  hm->registerHisto(h);
140  }
141  }
142 
144 }
145 
146 //_____________________________________________________________________
148 {
149  // load nodes
150  auto res = load_nodes(topNode);
151  if (res != Fun4AllReturnCodes::EVENT_OK) return res;
152  // run evaluation
155 }
156 
157 //________________________________________________________________________
159 {
160  return std::string("h_") + Name() + std::string("_");
161 }
162 
163 //________________________________________________________________________
165 {
166  m_hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
167  if (!m_hitsets)
168  {
169  std::cout << PHWHERE << " ERROR: Can't find TrkrHitSetContainer." << std::endl;
171  }
172 
173  m_surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
174  if (!m_surfmaps)
175  {
176  std::cout << PHWHERE << "Error: can't find Acts surface maps"
177  << std::endl;
179  }
180 
181  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
182  if (!m_tGeometry)
183  {
184  std::cout << PHWHERE << "No acts tracking geometry, exiting."
185  << std::endl;
187  }
188 
189  m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
190  if (!m_cluster_map)
191  {
192  std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTER" << std::endl;
194  }
195 
196  m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
197  if (!m_cluster_hit_map)
198  {
199  std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTERHITASSOC" << std::endl;
201  }
202 
203  m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
204  if (!m_hit_truth_map)
205  {
206  std::cout << PHWHERE << " unable to find DST node TRKR_HITTRUTHASSOC" << std::endl;
208  }
209 
210  m_g4hits_mvtx = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MVTX");
211  if (!m_g4hits_mvtx)
212  {
213  std::cout << PHWHERE << " unable to find DST node G4HIT_MVTX" << std::endl;
215  }
216 
218 }
219 
220 //________________________________________________________________________
222 {
223  // histogram manager
225  assert(hm);
226 
227  // load relevant histograms
228  struct HistogramList
229  {
230  TH1* drphi = nullptr;
231  TH1* rphi_error = nullptr;
232  TH1* phi_pulls = nullptr;
233 
234  TH1* dz = nullptr;
235  TH1* z_error = nullptr;
236  TH1* z_pulls = nullptr;
237 
238  TH1* csize = nullptr;
239  TH1* csize_phi = nullptr;
240  TH1* csize_z = nullptr;
241  };
242 
243  using HistogramMap = std::map<int, HistogramList>;
244  HistogramMap histograms;
245 
246  for (const auto& layer : m_layers)
247  {
248  HistogramList h;
249  h.drphi = dynamic_cast<TH1*>(hm->getHisto(Form("%sdrphi_%i", get_histo_prefix().c_str(), layer)));
250  h.rphi_error = dynamic_cast<TH1*>(hm->getHisto(Form("%srphi_error_%i", get_histo_prefix().c_str(), layer)));
251  h.phi_pulls = dynamic_cast<TH1*>(hm->getHisto(Form("%sphi_pulls_%i", get_histo_prefix().c_str(), layer)));
252 
253  h.dz = dynamic_cast<TH1*>(hm->getHisto(Form("%sdz_%i", get_histo_prefix().c_str(), layer)));
254  h.z_error = dynamic_cast<TH1*>(hm->getHisto(Form("%sz_error_%i", get_histo_prefix().c_str(), layer)));
255  h.z_pulls = dynamic_cast<TH1*>(hm->getHisto(Form("%sz_pulls_%i", get_histo_prefix().c_str(), layer)));
256 
257  h.csize = dynamic_cast<TH1*>(hm->getHisto(Form("%sclus_size_%i", get_histo_prefix().c_str(), layer)));
258  h.csize_phi = dynamic_cast<TH1*>(hm->getHisto(Form("%sclus_size_phi_%i", get_histo_prefix().c_str(), layer)));
259  h.csize_z = dynamic_cast<TH1*>(hm->getHisto(Form("%sclus_size_z_%i", get_histo_prefix().c_str(), layer)));
260 
261  histograms.insert(std::make_pair(layer, h));
262  }
263 
264  ActsTransformations transformer;
265  auto hitsetrange = m_hitsets->getHitSets(TrkrDefs::TrkrId::mvtxId);
266  for (auto hitsetitr = hitsetrange.first;
267  hitsetitr != hitsetrange.second;
268  ++hitsetitr)
269  {
270  auto range = m_cluster_map->getClusters(hitsetitr->first);
271  for (auto clusterIter = range.first; clusterIter != range.second; ++clusterIter)
272  {
273  // get cluster key, detector id and check
274  const auto& key = clusterIter->first;
275  // get cluster
276  const auto& cluster = clusterIter->second;
277 
278  const auto global = transformer.getGlobalPosition(cluster, m_surfmaps,
279  m_tGeometry);
280 
281  // get relevant cluster information
282  const auto r_cluster = QAG4Util::get_r(global(0), global(1));
283  const auto z_cluster = global(2);
284  const auto phi_cluster = (float) std::atan2(global(1), global(0));
285  const auto phi_error = cluster->getRPhiError() / r_cluster;
286  const auto z_error = cluster->getZError();
287 
288  // find associated g4hits
289  const auto g4hits = find_g4hits(key);
290 
291  // get relevant truth information
292  const auto x_truth = QAG4Util::interpolate<&PHG4Hit::get_x>(g4hits, r_cluster);
293  const auto y_truth = QAG4Util::interpolate<&PHG4Hit::get_y>(g4hits, r_cluster);
294  const auto z_truth = QAG4Util::interpolate<&PHG4Hit::get_z>(g4hits, r_cluster);
295  const auto phi_truth = std::atan2(y_truth, x_truth);
296 
297  const auto dphi = QAG4Util::delta_phi(phi_cluster, phi_truth);
298  const auto dz = z_cluster - z_truth;
299 
300  // get layer, get histograms
301  const auto layer = TrkrDefs::getLayer(key);
302  const auto hiter = histograms.find(layer);
303  if (hiter == histograms.end()) continue;
304 
305  // fill phi residuals, errors and pulls
306  auto fill = [](TH1* h, float value) { if( h ) h->Fill( value ); };
307  fill(hiter->second.drphi, r_cluster * dphi);
308  fill(hiter->second.rphi_error, r_cluster * phi_error);
309  fill(hiter->second.phi_pulls, dphi / phi_error);
310 
311  // fill z residuals, errors and pulls
312  fill(hiter->second.dz, dz);
313  fill(hiter->second.z_error, z_error);
314  fill(hiter->second.z_pulls, dz / z_error);
315 
316  // cluster sizes
317  // get associated hits
318  const auto hit_range = m_cluster_hit_map->getHits(key);
319  fill(hiter->second.csize, std::distance(hit_range.first, hit_range.second));
320 
321  std::set<int> phibins;
322  std::set<int> zbins;
323  for (auto hit_iter = hit_range.first; hit_iter != hit_range.second; ++hit_iter)
324  {
325  const auto& hit_key = hit_iter->second;
326  phibins.insert(MvtxDefs::getRow(hit_key));
327  zbins.insert(MvtxDefs::getCol(hit_key));
328  }
329 
330  fill(hiter->second.csize_phi, phibins.size());
331  fill(hiter->second.csize_z, zbins.size());
332  }
333  }
334 }
335 //_____________________________________________________________________
337 {
338  // find hitset associated to cluster
339  G4HitSet out;
340  const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
341 
342  // loop over hits associated to clusters
343  const auto range = m_cluster_hit_map->getHits(cluster_key);
344  for (auto iter = range.first; iter != range.second; ++iter)
345  {
346  // hit key
347  const auto& hit_key = iter->second;
348 
349  // store hits to g4hit associations
350  TrkrHitTruthAssoc::MMap g4hit_map;
351  m_hit_truth_map->getG4Hits(hitset_key, hit_key, g4hit_map);
352 
353  // find corresponding g4 hist
354  for (auto truth_iter = g4hit_map.begin(); truth_iter != g4hit_map.end(); ++truth_iter)
355  {
356  // g4hit key
357  const auto g4hit_key = truth_iter->second.second;
358 
359  // g4 hit
360  PHG4Hit* g4hit = (TrkrDefs::getTrkrId(hitset_key) == TrkrDefs::mvtxId) ? m_g4hits_mvtx->findHit(g4hit_key) : nullptr;
361 
362  // insert in set
363  if (g4hit) out.insert(g4hit);
364  }
365  }
366 
367  return out;
368 }