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