Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4SimulationIntt.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4SimulationIntt.cc
1 #include "QAG4SimulationIntt.h"
2 #include <qautils/QAUtil.h>
3 #include <qautils/QAHistManagerDef.h>
4 
6 
7 #include <g4main/PHG4Hit.h>
9 
11 
12 #include <trackbase/ActsGeometry.h>
14 #include <trackbase/InttDefs.h>
15 #include <trackbase/TrkrCluster.h>
18 #include <trackbase/TrkrDefs.h> // for getTrkrId, getHit...
20 
23 #include <fun4all/SubsysReco.h> // for SubsysReco
24 
25 #include <phool/getClass.h>
26 #include <phool/phool.h> // for PHWHERE
27 
28 #include <TH1.h>
29 #include <TString.h> // for Form
30 
31 #include <cassert>
32 #include <cmath>
33 #include <iostream> // for operator<<, basic...
34 #include <iterator> // for distance
35 #include <map> // for map
36 #include <utility> // for pair, make_pair
37 
38 //________________________________________________________________________
40  : SubsysReco(name)
41 {
42 }
43 
44 //________________________________________________________________________
46 {
47  // prevent multiple creations of histograms
48  if (m_initialized)
50  else
51  m_initialized = true;
52 
53  // find intt geometry
54  auto geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_INTT");
55  if (!geom_container)
56  {
57  std::cout << PHWHERE << " unable to find DST node CYLINDERGEOM_INTT" << std::endl;
59  }
60 
61  // get layers from intt geometry
62  const auto range = geom_container->get_begin_end();
63  for (auto iter = range.first; iter != range.second; ++iter)
64  {
65  m_layers.insert(iter->first);
66  }
67 
68  // histogram manager
70  assert(hm);
71 
72  // create histograms
73  for (const auto& layer : m_layers)
74  {
75  if (Verbosity()) std::cout << PHWHERE << " adding layer " << layer << std::endl;
76  {
77  // rphi residuals (cluster - truth)
78  auto h = new TH1F(Form("%sdrphi_%i", get_histo_prefix().c_str(), layer), Form("INTT r#Delta#phi_{cluster-truth} layer_%i", layer), 100, -1e-2, 1e-2);
79  h->GetXaxis()->SetTitle("r#Delta#phi_{cluster-truth} (cm)");
80  hm->registerHisto(h);
81  }
82 
83  {
84  // rphi cluster errors
85  auto h = new TH1F(Form("%srphi_error_%i", get_histo_prefix().c_str(), layer), Form("INTT r#Delta#phi error layer_%i", layer), 100, 0, 1e-2);
86  h->GetXaxis()->SetTitle("r#Delta#phi error (cm)");
87  hm->registerHisto(h);
88  }
89 
90  {
91  // phi pulls (cluster - truth)
92  auto h = new TH1F(Form("%sphi_pulls_%i", get_histo_prefix().c_str(), layer), Form("INTT #Delta#phi_{cluster-truth}/#sigma#phi layer_%i", layer), 100, -3, 3);
93  h->GetXaxis()->SetTitle("#Delta#phi_{cluster-truth}/#sigma#phi");
94  hm->registerHisto(h);
95  }
96 
97  {
98  // z residuals (cluster - truth)
99  auto h = new TH1F(Form("%sdz_%i", get_histo_prefix().c_str(), layer), Form("INTT #Deltaz_{cluster-truth} layer_%i", layer), 100, -2.5, 2.5);
100  h->GetXaxis()->SetTitle("#Deltaz_{cluster-truth} (cm)");
101  hm->registerHisto(h);
102  }
103 
104  {
105  // z cluster errors
106  auto h = new TH1F(Form("%sz_error_%i", get_histo_prefix().c_str(), layer), Form("INTT z error layer_%i", layer), 100, 0, 2.5);
107  h->GetXaxis()->SetTitle("z error (cm)");
108  hm->registerHisto(h);
109  }
110 
111  {
112  // z pulls (cluster - truth)
113  auto h = new TH1F(Form("%sz_pulls_%i", get_histo_prefix().c_str(), layer), Form("INTT #Deltaz_{cluster-truth}/#sigmaz layer_%i", layer), 100, -3, 3);
114  h->GetXaxis()->SetTitle("#Deltaz_{cluster-truth}/#sigmaz");
115  hm->registerHisto(h);
116  }
117 
118  {
119  // total cluster size
120  auto h = new TH1F(Form("%sclus_size_%i", get_histo_prefix().c_str(), layer), Form("INTT cluster size layer_%i", layer), 10, 0, 10);
121  h->GetXaxis()->SetTitle("csize");
122  hm->registerHisto(h);
123  }
124 
125  {
126  // cluster size in phi
127  auto h = new TH1F(Form("%sclus_size_phi_%i", get_histo_prefix().c_str(), layer), Form("INTT cluster size (#phi) layer_%i", layer), 10, 0, 10);
128  h->GetXaxis()->SetTitle("csize_{#phi}");
129  hm->registerHisto(h);
130  }
131 
132  {
133  // cluster size in z
134  auto h = new TH1F(Form("%sclus_size_z_%i", get_histo_prefix().c_str(), layer), Form("INTT cluster size (z) layer_%i", layer), 10, 0, 10);
135  h->GetXaxis()->SetTitle("csize_{z}");
136  hm->registerHisto(h);
137  }
138  }
139 
141 }
142 
143 //_____________________________________________________________________
145 {
146  // load nodes
147  auto res = load_nodes(topNode);
148  if (res != Fun4AllReturnCodes::EVENT_OK) return res;
149 
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_intt = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_INTT");
194  if (!m_g4hits_intt)
195  {
196  std::cout << PHWHERE << " unable to find DST node G4HIT_INTT" << 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 
257  const auto global = m_tGeometry->getGlobalPosition(key, cluster);
258 
259  // get relevant cluster information
260  const auto r_cluster = QAG4Util::get_r(global(0), global(1));
261  const auto z_cluster = global(2);
262  const auto phi_cluster = (float) std::atan2(global(1), global(0));
263 
264  double phi_error = 0;
265  double z_error = 0;
266 
267  phi_error = cluster->getRPhiError() / r_cluster;
268  z_error = cluster->getZError();
269 
270 
271  // find associated g4hits
272  const auto g4hits = find_g4hits(key);
273 
274  // get relevant truth information
275  const auto x_truth = QAG4Util::interpolate<&PHG4Hit::get_x>(g4hits, r_cluster);
276  const auto y_truth = QAG4Util::interpolate<&PHG4Hit::get_y>(g4hits, r_cluster);
277  const auto z_truth = QAG4Util::interpolate<&PHG4Hit::get_z>(g4hits, r_cluster);
278  const auto phi_truth = std::atan2(y_truth, x_truth);
279 
280  const auto dphi = QAG4Util::delta_phi(phi_cluster, phi_truth);
281  const auto dz = z_cluster - z_truth;
282 
283  // get layer, get histograms
284  const auto layer = TrkrDefs::getLayer(key);
285  const auto hiter = histograms.find(layer);
286  if (hiter == histograms.end()) continue;
287 
288  // fill histograms
289  auto fill = [](TH1* h, float value)
290  { if( h ) h->Fill( value ); };
291  fill(hiter->second.drphi, r_cluster * dphi);
292  fill(hiter->second.rphi_error, r_cluster * phi_error);
293  fill(hiter->second.phi_pulls, dphi / phi_error);
294 
295  fill(hiter->second.dz, dz);
296  fill(hiter->second.z_error, z_error);
297  fill(hiter->second.z_pulls, dz / z_error);
298 
299  // cluster sizes
300  // get associated hits
301  const auto hit_range = m_cluster_hit_map->getHits(key);
302  fill(hiter->second.csize, std::distance(hit_range.first, hit_range.second));
303 
304  std::set<int> phibins;
305  std::set<int> zbins;
306  for (auto hit_iter = hit_range.first; hit_iter != hit_range.second; ++hit_iter)
307  {
308  const auto& hit_key = hit_iter->second;
309  phibins.insert(InttDefs::getRow(hit_key));
310  zbins.insert(InttDefs::getCol(hit_key));
311  }
312 
313  fill(hiter->second.csize_phi, phibins.size());
314  fill(hiter->second.csize_z, zbins.size());
315  }
316  }
317 }
318 
319 //_____________________________________________________________________
321 {
322  // find hitset associated to cluster
323  G4HitSet out;
324  const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
325 
326  // loop over hits associated to clusters
327  const auto range = m_cluster_hit_map->getHits(cluster_key);
328  for (auto iter = range.first; iter != range.second; ++iter)
329  {
330  // hit key
331  const auto& hit_key = iter->second;
332 
333  // store hits to g4hit associations
334  TrkrHitTruthAssoc::MMap g4hit_map;
335  m_hit_truth_map->getG4Hits(hitset_key, hit_key, g4hit_map);
336 
337  // find corresponding g4 hist
338  for (auto truth_iter = g4hit_map.begin(); truth_iter != g4hit_map.end(); ++truth_iter)
339  {
340  // g4hit key
341  const auto g4hit_key = truth_iter->second.second;
342 
343  // g4 hit
344  PHG4Hit* g4hit = (TrkrDefs::getTrkrId(hitset_key) == TrkrDefs::inttId) ? m_g4hits_intt->findHit(g4hit_key) : nullptr;
345 
346  // insert in set
347  if (g4hit) out.insert(g4hit);
348  }
349  }
350 
351  return out;
352 }