58 struct interpolation_data_t
60 using list = std::vector<interpolation_data_t>;
65 double px()
const {
return momentum.x(); }
66 double py()
const {
return momentum.y(); }
67 double pz()
const {
return momentum.z(); }
75 template <
double (
interpolation_data_t::*accessor)() const>
76 double interpolate(
const interpolation_data_t::list&
hits,
double y_extrap)
87 for (
const auto& hit : hits)
89 const double x = (hit.*accessor)();
90 const double w = hit.weight;
94 const double y = hit.y();
103 if (!
valid)
return NAN;
105 const auto alpha = (sw * swyx - swy * swx);
106 const auto beta = (swy2 * swx - swy * swyx);
130 auto geom_container = findNode::getClass<PHG4CylinderGeomContainer>(topNode,
"CYLINDERGEOM_MICROMEGAS_FULL");
133 std::cout <<
PHWHERE <<
" unable to find DST node CYLINDERGEOM_MICROMEGAS_FULL" << std::endl;
138 const auto range = geom_container->get_begin_end();
139 for (
auto iter = range.first; iter != range.second; ++iter)
158 const auto segmentation_type = layergeom->get_segmentation_type();
164 Form(
"micromegas ADC distribution layer_%i",
layer),
166 h->GetXaxis()->SetTitle(
"ADC");
167 hm->registerHisto(
h);
172 const double max_residual = is_segmentation_phi ? 0.04 : 0.08;
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);
183 const double max_error = is_segmentation_phi ? 0.04 : 0.08;
185 Form(
"micromegas %s error layer_%i",
186 is_segmentation_phi ?
"r#Delta#phi" :
"#Deltaz",
layer),
188 h->GetXaxis()->SetTitle(Form(
"%s error (cm)", is_segmentation_phi ?
"r#Delta#phi" :
"#Deltaz"));
189 hm->registerHisto(
h);
195 Form(
"micromegas %s layer_%i",
196 is_segmentation_phi ?
"#Delta#phi/#sigma#phi" :
"#Deltaz/#sigmaz",
layer),
198 h->GetXaxis()->SetTitle(is_segmentation_phi ?
"#Delta#phi/#sigma#phi" :
"#Deltaz/#sigmaz");
199 hm->registerHisto(
h);
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);
229 return std::string(
"h_") +
Name() + std::string(
"_");
235 m_surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,
"ActsSurfaceMaps");
238 std::cout <<
PHWHERE <<
"Error: can't find Acts surface maps"
243 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
246 std::cout <<
PHWHERE <<
"No acts tracking geometry, exiting."
251 m_micromegas_geonode = findNode::getClass<PHG4CylinderGeomContainer>(topNode,
"CYLINDERGEOM_MICROMEGAS_FULL");
254 std::cout <<
PHWHERE <<
" ERROR: Can't find CYLINDERGEOM_MICROMEGAS_FULL." << std::endl;
258 m_hitsets = findNode::getClass<TrkrHitSetContainer>(topNode,
"TRKR_HITSET");
261 std::cout <<
PHWHERE <<
" ERROR: Can't find TrkrHitSetContainer." << std::endl;
265 m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
268 std::cout <<
PHWHERE <<
" unable to find DST node TRKR_CLUSTER" << std::endl;
272 m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
275 std::cout <<
PHWHERE <<
" unable to find DST node TRKR_CLUSTERHITASSOC" << std::endl;
279 m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
282 std::cout <<
PHWHERE <<
" unable to find DST node TRKR_HITTRUTHASSOC" << std::endl;
289 std::cout <<
PHWHERE <<
" unable to find DST node G4HIT_MVTX" << std::endl;
309 using HistogramMap = std::map<int, HistogramList>;
310 HistogramMap histograms;
316 histograms.insert(std::make_pair(
layer, h));
321 for (
auto hitsetiter = hitsetrange.first; hitsetiter != hitsetrange.second; ++hitsetiter)
324 const auto hitsetkey = hitsetiter->first;
328 const auto hiter = histograms.find(
layer);
329 if (hiter == histograms.end())
continue;
336 for (
auto hit_it = hitrange.first; hit_it != hitrange.second; ++hit_it)
339 auto fill = [](TH1*
h,
float value) {
if( h ) h->Fill(
value ); };
340 fill(hiter->second.adc, hit_it->second->getAdc());
355 TH1* residual =
nullptr;
356 TH1* residual_error =
nullptr;
357 TH1* pulls =
nullptr;
359 TH1* csize =
nullptr;
362 using HistogramMap = std::map<int, HistogramList>;
363 HistogramMap histograms;
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)));
373 histograms.insert(std::make_pair(
layer, h));
380 for (
auto hitsetiter = hitsetrange.first; hitsetiter != hitsetrange.second; ++hitsetiter)
383 const auto hitsetkey = hitsetiter->first;
393 if (!layergeom)
continue;
396 const auto hiter = histograms.find(
layer);
397 if (hiter == histograms.end())
continue;
401 for (
auto clusterIter = range.first; clusterIter != range.second; ++clusterIter)
404 const auto& key = clusterIter->first;
407 const auto& cluster = clusterIter->second;
416 const auto rphi_error = cluster->getRPhiError();
417 const auto z_error = cluster->getZError();
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);
427 interpolation_data_t::list
hits;
428 for (
const auto& g4hit : g4hits)
430 const auto weight = g4hit->get_edep();
431 for (
int i = 0; i < 2; ++i)
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);
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);
441 hits.push_back({.position = g4hit_local, .momentum = momentum_local, .weight =
weight});
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));
453 auto fill = [](TH1*
h,
float value) {
if( h ) h->Fill(
value ); };
454 switch (segmentation_type)
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);
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);
478 fill(hiter->second.csize, std::distance(hit_range.first, hit_range.second));
491 for (
auto iter = range.first; iter != range.second; ++iter)
494 const auto& hit_key = iter->second;
501 for (
auto truth_iter = g4hit_map.begin(); truth_iter != g4hit_map.end(); ++truth_iter)
504 const auto g4hit_key = truth_iter->second.second;
510 if (g4hit) out.insert(g4hit);