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