ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
QAG4SimulationMicromegas.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file QAG4SimulationMicromegas.cc
2 
3 #include "QAHistManagerDef.h"
4 
5 #include <g4detectors/PHG4CylinderGeom.h> // for PHG4CylinderGeom
7 
8 #include <g4main/PHG4Hit.h>
10 
11 #include <micromegas/MicromegasDefs.h> // for SegmentationType
12 
14 
17 #include <trackbase/TrkrCluster.h>
20 #include <trackbase/TrkrDefs.h> // for getTrkrId, getHit...
21 #include <trackbase/TrkrHit.h>
22 #include <trackbase/TrkrHitSet.h>
25 
28 #include <fun4all/SubsysReco.h> // for SubsysReco
29 
31 #include <phool/getClass.h>
32 #include <phool/phool.h> // for PHWHERE
33 
34 #include <TAxis.h> // for TAxis
35 #include <TH1.h>
36 #include <TString.h> // for Form
37 #include <TVector3.h>
38 
39 #include <cassert>
40 #include <cmath> // for cos, sin, NAN
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 //_____________________________________________________________________
48 namespace
49 {
51  template <class T>
52  inline constexpr T square(T x)
53  {
54  return x * x;
55  }
56 
58  struct interpolation_data_t
59  {
60  using list = std::vector<interpolation_data_t>;
61  double x() const { return position.x(); }
62  double y() const { return position.y(); }
63  double z() const { return position.z(); }
64 
65  double px() const { return momentum.x(); }
66  double py() const { return momentum.y(); }
67  double pz() const { return momentum.z(); }
68 
69  TVector3 position;
70  TVector3 momentum;
71  double weight = 1;
72  };
73 
75  template <double (interpolation_data_t::*accessor)() const>
76  double interpolate(const interpolation_data_t::list& hits, double y_extrap)
77  {
78  // calculate all terms needed for the interpolation
79  // need to use double everywhere here due to numerical divergences
80  double sw = 0;
81  double swy = 0;
82  double swy2 = 0;
83  double swx = 0;
84  double swyx = 0;
85 
86  bool valid(false);
87  for (const auto& hit : hits)
88  {
89  const double x = (hit.*accessor)();
90  const double w = hit.weight;
91  if (w <= 0) continue;
92 
93  valid = true;
94  const double y = hit.y();
95 
96  sw += w;
97  swy += w * y;
98  swy2 += w * square(y);
99  swx += w * x;
100  swyx += w * x * y;
101  }
102 
103  if (!valid) return NAN;
104 
105  const auto alpha = (sw * swyx - swy * swx);
106  const auto beta = (swy2 * swx - swy * swyx);
107  const auto denom = (sw * swy2 - square(swy));
108 
109  return (alpha * y_extrap + beta) / denom;
110  }
111 
112 } // namespace
113 
114 //________________________________________________________________________
116  : SubsysReco(name)
117 {
118 }
119 
120 //________________________________________________________________________
122 {
123  // prevent multiple creations of histograms
124  if (m_initialized)
126  else
127  m_initialized = true;
128 
129  // find mvtx geometry
130  auto geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
131  if (!geom_container)
132  {
133  std::cout << PHWHERE << " unable to find DST node CYLINDERGEOM_MICROMEGAS_FULL" << std::endl;
135  }
136 
137  // get layers from mvtx geometry
138  const auto range = geom_container->get_begin_end();
139  for (auto iter = range.first; iter != range.second; ++iter)
140  {
141  m_layers.insert(iter->first);
142  }
143 
144  // histogram manager
146  assert(hm);
147 
148  // create histograms
149  for (const auto& layer : m_layers)
150  {
151  if (Verbosity()) std::cout << PHWHERE << " adding layer " << layer << std::endl;
152 
153  // get layer geometry
154  const auto layergeom = dynamic_cast<CylinderGeomMicromegas*>(geom_container->GetLayerGeom(layer));
155  assert(layergeom);
156 
157  // get segmentation type
158  const auto segmentation_type = layergeom->get_segmentation_type();
159  const bool is_segmentation_phi = (segmentation_type == MicromegasDefs::SegmentationType::SEGMENTATION_PHI);
160 
161  {
162  // ADC distributions
163  auto h = new TH1F(Form("%sadc_%i", get_histo_prefix().c_str(), layer),
164  Form("micromegas ADC distribution layer_%i", layer),
165  1024, 0, 1024);
166  h->GetXaxis()->SetTitle("ADC");
167  hm->registerHisto(h);
168  }
169 
170  {
171  // residuals (cluster - truth)
172  const double max_residual = is_segmentation_phi ? 0.04 : 0.08;
173  auto h = new TH1F(Form("%sresidual_%i", get_histo_prefix().c_str(), layer),
174  Form("micromegas %s_{cluster-truth} layer_%i",
175  is_segmentation_phi ? "r#Delta#phi" : "#Deltaz", layer),
176  100, -max_residual, max_residual);
177  h->GetXaxis()->SetTitle(Form("%s_{cluster-truth} (cm)", is_segmentation_phi ? "r#Delta#phi" : "#Deltaz"));
178  hm->registerHisto(h);
179  }
180 
181  {
182  // cluster errors
183  const double max_error = is_segmentation_phi ? 0.04 : 0.08;
184  auto h = new TH1F(Form("%sresidual_error_%i", get_histo_prefix().c_str(), layer),
185  Form("micromegas %s error layer_%i",
186  is_segmentation_phi ? "r#Delta#phi" : "#Deltaz", layer),
187  100, 0, max_error);
188  h->GetXaxis()->SetTitle(Form("%s error (cm)", is_segmentation_phi ? "r#Delta#phi" : "#Deltaz"));
189  hm->registerHisto(h);
190  }
191 
192  {
193  // pulls (cluster - truth)
194  auto h = new TH1F(Form("%scluster_pulls_%i", get_histo_prefix().c_str(), layer),
195  Form("micromegas %s layer_%i",
196  is_segmentation_phi ? "#Delta#phi/#sigma#phi" : "#Deltaz/#sigmaz", layer),
197  100, -5, 5);
198  h->GetXaxis()->SetTitle(is_segmentation_phi ? "#Delta#phi/#sigma#phi" : "#Deltaz/#sigmaz");
199  hm->registerHisto(h);
200  }
201 
202  {
203  // cluster size
204  auto h = new TH1F(Form("%sclus_size_%i", get_histo_prefix().c_str(), layer), Form("micromegas cluster size layer_%i", layer), 20, 0, 20);
205  h->GetXaxis()->SetTitle("csize");
206  hm->registerHisto(h);
207  }
208  }
209 
211 }
212 
213 //_____________________________________________________________________
215 {
216  // load nodes
217  auto res = load_nodes(topNode);
218  if (res != Fun4AllReturnCodes::EVENT_OK) return res;
219 
220  // run evaluation
221  evaluate_hits();
224 }
225 
226 //________________________________________________________________________
228 {
229  return std::string("h_") + Name() + std::string("_");
230 }
231 
232 //________________________________________________________________________
234 {
235  m_surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode, "ActsSurfaceMaps");
236  if (!m_surfmaps)
237  {
238  std::cout << PHWHERE << "Error: can't find Acts surface maps"
239  << std::endl;
241  }
242 
243  m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode, "ActsTrackingGeometry");
244  if (!m_tGeometry)
245  {
246  std::cout << PHWHERE << "No acts tracking geometry, exiting."
247  << std::endl;
249  }
250 
251  m_micromegas_geonode = findNode::getClass<PHG4CylinderGeomContainer>(topNode, "CYLINDERGEOM_MICROMEGAS_FULL");
253  {
254  std::cout << PHWHERE << " ERROR: Can't find CYLINDERGEOM_MICROMEGAS_FULL." << std::endl;
256  }
257 
258  m_hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
259  if (!m_hitsets)
260  {
261  std::cout << PHWHERE << " ERROR: Can't find TrkrHitSetContainer." << std::endl;
263  }
264 
265  m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
266  if (!m_cluster_map)
267  {
268  std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTER" << std::endl;
270  }
271 
272  m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
273  if (!m_cluster_hit_map)
274  {
275  std::cout << PHWHERE << " unable to find DST node TRKR_CLUSTERHITASSOC" << std::endl;
277  }
278 
279  m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode, "TRKR_HITTRUTHASSOC");
280  if (!m_hit_truth_map)
281  {
282  std::cout << PHWHERE << " unable to find DST node TRKR_HITTRUTHASSOC" << std::endl;
284  }
285 
286  m_g4hits_micromegas = findNode::getClass<PHG4HitContainer>(topNode, "G4HIT_MICROMEGAS");
287  if (!m_g4hits_micromegas)
288  {
289  std::cout << PHWHERE << " unable to find DST node G4HIT_MVTX" << std::endl;
291  }
292 
294 }
295 
296 //________________________________________________________________________
298 {
299  // histogram manager
301  assert(hm);
302 
303  // load relevant histograms
304  struct HistogramList
305  {
306  TH1* adc = nullptr;
307  };
308 
309  using HistogramMap = std::map<int, HistogramList>;
310  HistogramMap histograms;
311 
312  for (const auto& layer : m_layers)
313  {
314  HistogramList h;
315  h.adc = dynamic_cast<TH1*>(hm->getHisto(Form("%sadc_%i", get_histo_prefix().c_str(), layer)));
316  histograms.insert(std::make_pair(layer, h));
317  }
318 
319  // loop over hitsets
320  const auto hitsetrange = m_hitsets->getHitSets(TrkrDefs::TrkrId::micromegasId);
321  for (auto hitsetiter = hitsetrange.first; hitsetiter != hitsetrange.second; ++hitsetiter)
322  {
323  // get hitsetkey and layer
324  const auto hitsetkey = hitsetiter->first;
325  const auto layer = TrkrDefs::getLayer(hitsetkey);
326 
327  // get relevant histograms
328  const auto hiter = histograms.find(layer);
329  if (hiter == histograms.end()) continue;
330 
331  // get all of the hits from this hitset
332  TrkrHitSet* hitset = hitsetiter->second;
333  TrkrHitSet::ConstRange hitrange = hitset->getHits();
334 
335  // loop over hits
336  for (auto hit_it = hitrange.first; hit_it != hitrange.second; ++hit_it)
337  {
338  // store ADC
339  auto fill = [](TH1* h, float value) { if( h ) h->Fill( value ); };
340  fill(hiter->second.adc, hit_it->second->getAdc());
341  }
342  }
343 }
344 
345 //________________________________________________________________________
347 {
348  // histogram manager
350  assert(hm);
351 
352  // load relevant histograms
353  struct HistogramList
354  {
355  TH1* residual = nullptr;
356  TH1* residual_error = nullptr;
357  TH1* pulls = nullptr;
358 
359  TH1* csize = nullptr;
360  };
361 
362  using HistogramMap = std::map<int, HistogramList>;
363  HistogramMap histograms;
364 
365  for (const auto& layer : m_layers)
366  {
367  HistogramList h;
368  h.residual = dynamic_cast<TH1*>(hm->getHisto(Form("%sresidual_%i", get_histo_prefix().c_str(), layer)));
369  h.residual_error = dynamic_cast<TH1*>(hm->getHisto(Form("%sresidual_error_%i", get_histo_prefix().c_str(), layer)));
370  h.pulls = dynamic_cast<TH1*>(hm->getHisto(Form("%scluster_pulls_%i", get_histo_prefix().c_str(), layer)));
371  h.csize = dynamic_cast<TH1*>(hm->getHisto(Form("%sclus_size_%i", get_histo_prefix().c_str(), layer)));
372 
373  histograms.insert(std::make_pair(layer, h));
374  }
375 
376  ActsTransformations transformer;
377 
378  // loop over hitsets
379  const auto hitsetrange = m_hitsets->getHitSets(TrkrDefs::TrkrId::micromegasId);
380  for (auto hitsetiter = hitsetrange.first; hitsetiter != hitsetrange.second; ++hitsetiter)
381  {
382  // get hitsetkey, layer and tileid
383  const auto hitsetkey = hitsetiter->first;
384 
385  // get layer
386  const auto layer = TrkrDefs::getLayer(hitsetkey);
387 
388  // get tileid
389  const auto tileid = MicromegasDefs::getTileId(hitsetkey);
390 
391  // load geometry
392  const auto layergeom = dynamic_cast<CylinderGeomMicromegas*>(m_micromegas_geonode->GetLayerGeom(layer));
393  if (!layergeom) continue;
394 
395  // get relevant histograms
396  const auto hiter = histograms.find(layer);
397  if (hiter == histograms.end()) continue;
398 
399  // get associated clusters
400  const auto range = m_cluster_map->getClusters(hitsetkey);
401  for (auto clusterIter = range.first; clusterIter != range.second; ++clusterIter)
402  {
403  // get cluster key
404  const auto& key = clusterIter->first;
405 
406  // get cluster
407  const auto& cluster = clusterIter->second;
408 
409  const auto global = transformer.getGlobalPosition(cluster, m_surfmaps,
410  m_tGeometry);
411 
412  // get segmentation type
413  const auto segmentation_type = MicromegasDefs::getSegmentationType(key);
414 
415  // get relevant cluster information
416  const auto rphi_error = cluster->getRPhiError();
417  const auto z_error = cluster->getZError();
418 
419  // convert cluster position to local tile coordinates
420  const TVector3 cluster_world(global(0), global(1), global(2));
421  const TVector3 cluster_local = layergeom->get_local_from_world_coords(tileid, cluster_world);
422 
423  // find associated g4hits
424  const auto g4hits = find_g4hits(key);
425 
426  // convert hits to list of interpolation_data_t
427  interpolation_data_t::list hits;
428  for (const auto& g4hit : g4hits)
429  {
430  const auto weight = g4hit->get_edep();
431  for (int i = 0; i < 2; ++i)
432  {
433  // convert position to local
434  TVector3 g4hit_world(g4hit->get_x(i), g4hit->get_y(i), g4hit->get_z(i));
435  TVector3 g4hit_local = layergeom->get_local_from_world_coords(tileid, g4hit_world);
436 
437  // convert momentum to local
438  TVector3 momentum_world(g4hit->get_px(i), g4hit->get_py(i), g4hit->get_pz(i));
439  TVector3 momentum_local = layergeom->get_local_from_world_vect(tileid, momentum_world);
440 
441  hits.push_back({.position = g4hit_local, .momentum = momentum_local, .weight = weight});
442  }
443  }
444 
445  // do position interpolation
446  const auto y_extrap = cluster_local.y();
447  const TVector3 interpolation_local(
448  interpolate<&interpolation_data_t::x>(hits, y_extrap),
449  interpolate<&interpolation_data_t::y>(hits, y_extrap),
450  interpolate<&interpolation_data_t::z>(hits, y_extrap));
451 
452  // fill phi residuals, errors and pulls
453  auto fill = [](TH1* h, float value) { if( h ) h->Fill( value ); };
454  switch (segmentation_type)
455  {
457  {
458  const auto drphi = cluster_local.x() - interpolation_local.x();
459  fill(hiter->second.residual, drphi);
460  fill(hiter->second.residual_error, rphi_error);
461  fill(hiter->second.pulls, drphi / rphi_error);
462  break;
463  }
464 
466  {
467  const auto dz = cluster_local.z() - interpolation_local.z();
468  fill(hiter->second.residual, dz);
469  fill(hiter->second.residual_error, z_error);
470  fill(hiter->second.pulls, dz / z_error);
471  break;
472  }
473  }
474 
475  // cluster size
476  // get associated hits
477  const auto hit_range = m_cluster_hit_map->getHits(key);
478  fill(hiter->second.csize, std::distance(hit_range.first, hit_range.second));
479  }
480  }
481 }
482 //_____________________________________________________________________
484 {
485  // find hitset associated to cluster
486  G4HitSet out;
487  const auto hitset_key = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
488 
489  // loop over hits associated to clusters
490  const auto range = m_cluster_hit_map->getHits(cluster_key);
491  for (auto iter = range.first; iter != range.second; ++iter)
492  {
493  // hit key
494  const auto& hit_key = iter->second;
495 
496  // store hits to g4hit associations
497  TrkrHitTruthAssoc::MMap g4hit_map;
498  m_hit_truth_map->getG4Hits(hitset_key, hit_key, g4hit_map);
499 
500  // find corresponding g4 hist
501  for (auto truth_iter = g4hit_map.begin(); truth_iter != g4hit_map.end(); ++truth_iter)
502  {
503  // g4hit key
504  const auto g4hit_key = truth_iter->second.second;
505 
506  // g4 hit
507  PHG4Hit* g4hit = (TrkrDefs::getTrkrId(hitset_key) == TrkrDefs::micromegasId) ? m_g4hits_micromegas->findHit(g4hit_key) : nullptr;
508 
509  // insert in set
510  if (g4hit) out.insert(g4hit);
511  }
512  }
513 
514  return out;
515 }