50 , m_truthContainer(nullptr)
74 std::cout <<
"QAG4SimulationTpc::InitRun - Fatal Error - "
75 <<
"unable to find node G4TruthInfo" << std::endl;
81 findNode::getClass<PHG4CylinderCellGeomContainer>(topNode,
"CYLINDERCELLGEOM_SVTX");
85 std::cout <<
PHWHERE <<
" unable to find DST node CYLINDERCELLGEOM_SVTX" << std::endl;
90 std::vector<int> region_layer_low = {7, 23, 39};
92 std::vector<int> region_layer_high = {22, 38, 54};
97 for (
auto iter = range.first; iter != range.second; ++iter)
101 for (
int region = 0; region < 3; ++region)
103 if (iter->first >= region_layer_low[region] && iter->first <= region_layer_high[region])
116 auto h =
new TH1F(Form(
"%sefficiency_0",
get_histo_prefix().c_str()), Form(
"TPC_truth_clusters"), 48, 7, 54);
117 h->GetXaxis()->SetTitle(
"sPHENIX layer");
118 hm->registerHisto(
h);
122 auto h =
new TH1F(Form(
"%sefficiency_1",
get_histo_prefix().c_str()), Form(
"TPC_matched_clusters"), 48, 7, 54);
123 h->GetXaxis()->SetTitle(
"sPHENIX layer");
124 hm->registerHisto(
h);
128 for (
int region = 0; region < 3; ++region)
130 if (
Verbosity()) std::cout <<
PHWHERE <<
" adding region " << region <<
" with layers " << region_layer_low[region] <<
" to " << region_layer_high[region] << std::endl;
133 auto h =
new TH1F(Form(
"%sdrphi_%i",
get_histo_prefix().c_str(), region), Form(
"TPC r#Delta#phi_{cluster-truth} region_%i", region), 100, -0.079, 0.075);
134 h->GetXaxis()->SetTitle(
"r#Delta#phi_{cluster-truth} (cm)");
135 hm->registerHisto(
h);
140 auto h =
new TH1F(Form(
"%srphi_error_%i",
get_histo_prefix().c_str(), region), Form(
"TPC r#Delta#phi error region_%i", region), 100, 0, 0.075);
141 h->GetXaxis()->SetTitle(
"r#Delta#phi error (cm)");
142 hm->registerHisto(
h);
147 auto h =
new TH1F(Form(
"%sphi_pulls_%i",
get_histo_prefix().c_str(), region), Form(
"TPC #Delta#phi_{cluster-truth}/#sigma#phi region_%i", region), 100, -3, 3);
148 h->GetXaxis()->SetTitle(
"#Delta#phi_{cluster-truth}/#sigma#phi (cm)");
149 hm->registerHisto(
h);
154 auto h =
new TH1F(Form(
"%sdz_%i",
get_histo_prefix().c_str(), region), Form(
"TPC #Deltaz_{cluster-truth} region_%i", region), 100, -0.19, 0.19);
155 h->GetXaxis()->SetTitle(
"#Delta#z_{cluster-truth} (cm)");
156 hm->registerHisto(
h);
161 auto h =
new TH1F(Form(
"%sz_error_%i",
get_histo_prefix().c_str(), region), Form(
"TPC z error region_%i", region), 100, 0, 0.18);
162 h->GetXaxis()->SetTitle(
"z error (cm)");
163 hm->registerHisto(
h);
168 auto h =
new TH1F(Form(
"%sz_pulls_%i",
get_histo_prefix().c_str(), region), Form(
"TPC #Deltaz_{cluster-truth}/#sigmaz region_%i", region), 100, -3, 3);
169 h->GetXaxis()->SetTitle(
"#Delta#z_{cluster-truth}/#sigmaz (cm)");
170 hm->registerHisto(
h);
175 auto h =
new TH1F(Form(
"%sclus_size_%i",
get_histo_prefix().c_str(), region), Form(
"TPC cluster size region_%i", region), 30, 0, 30);
176 h->GetXaxis()->SetTitle(
"csize");
177 hm->registerHisto(
h);
182 auto h =
new TH1F(Form(
"%sclus_size_phi_%i",
get_histo_prefix().c_str(), region), Form(
"TPC cluster size (#phi) region_%i", region), 10, 0, 10);
183 h->GetXaxis()->SetTitle(
"csize_{#phi}");
184 hm->registerHisto(
h);
189 auto h =
new TH1F(Form(
"%sclus_size_z_%i",
get_histo_prefix().c_str(), region), Form(
"TPC cluster size (z) region_%i", region), 12, 0, 12);
190 h->GetXaxis()->SetTitle(
"csize_{z}");
191 hm->registerHisto(
h);
216 return std::string(
"h_") +
Name() + std::string(
"_");
222 m_cluster_map = findNode::getClass<TrkrClusterContainer>(topNode,
"TRKR_CLUSTER");
225 std::cout <<
PHWHERE <<
" unable to find DST node TRKR_CLUSTER" << std::endl;
229 m_surfmaps = findNode::getClass<ActsSurfaceMaps>(topNode,
"ActsSurfaceMaps");
232 std::cout <<
PHWHERE <<
"Error: can't find Acts surface maps"
237 m_tGeometry = findNode::getClass<ActsTrackingGeometry>(topNode,
"ActsTrackingGeometry");
240 std::cout <<
PHWHERE <<
"No acts tracking geometry, exiting."
245 m_cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode,
"TRKR_CLUSTERHITASSOC");
248 std::cout <<
PHWHERE <<
" unable to find DST node TRKR_CLUSTERHITASSOC" << std::endl;
252 m_hit_truth_map = findNode::getClass<TrkrHitTruthAssoc>(topNode,
"TRKR_HITTRUTHASSOC");
255 std::cout <<
PHWHERE <<
" unable to find DST node TRKR_HITTRUTHASSOC" << std::endl;
259 m_g4hits_tpc = findNode::getClass<PHG4HitContainer>(topNode,
"G4HIT_TPC");
262 std::cout <<
PHWHERE <<
" unable to find DST node G4HIT_TPC" << std::endl;
282 TH1* h_eff0 =
dynamic_cast<TH1*
>(hm->getHisto(Form(
"%sefficiency_0",
get_histo_prefix().c_str())));
284 TH1* h_eff1 =
dynamic_cast<TH1*
>(hm->getHisto(Form(
"%sefficiency_1",
get_histo_prefix().c_str())));
290 TH1* drphi =
nullptr;
291 TH1* rphi_error =
nullptr;
292 TH1* phi_pulls =
nullptr;
295 TH1* z_error =
nullptr;
296 TH1* z_pulls =
nullptr;
298 TH1* csize =
nullptr;
299 TH1* csize_phi =
nullptr;
300 TH1* csize_z =
nullptr;
303 using HistogramMap = std::map<int, HistogramList>;
304 HistogramMap histograms;
306 for (
int region = 0; region < 3; ++region)
309 h.drphi =
dynamic_cast<TH1*
>(hm->getHisto(Form(
"%sdrphi_%i",
get_histo_prefix().c_str(), region)));
310 h.rphi_error =
dynamic_cast<TH1*
>(hm->getHisto(Form(
"%srphi_error_%i",
get_histo_prefix().c_str(), region)));
311 h.phi_pulls =
dynamic_cast<TH1*
>(hm->getHisto(Form(
"%sphi_pulls_%i",
get_histo_prefix().c_str(), region)));
313 h.dz =
dynamic_cast<TH1*
>(hm->getHisto(Form(
"%sdz_%i",
get_histo_prefix().c_str(), region)));
314 h.z_error =
dynamic_cast<TH1*
>(hm->getHisto(Form(
"%sz_error_%i",
get_histo_prefix().c_str(), region)));
315 h.z_pulls =
dynamic_cast<TH1*
>(hm->getHisto(Form(
"%sz_pulls_%i",
get_histo_prefix().c_str(), region)));
317 h.csize =
dynamic_cast<TH1*
>(hm->getHisto(Form(
"%sclus_size_%i",
get_histo_prefix().c_str(), region)));
318 h.csize_phi =
dynamic_cast<TH1*
>(hm->getHisto(Form(
"%sclus_size_phi_%i",
get_histo_prefix().c_str(), region)));
319 h.csize_z =
dynamic_cast<TH1*
>(hm->getHisto(Form(
"%sclus_size_z_%i",
get_histo_prefix().c_str(), region)));
321 histograms.insert(std::make_pair(region, h));
327 std::cout <<
PHWHERE <<
" get all truth clusters for primary particles " << std::endl;
335 iter != range.second;
341 float gflavor = g4particle->
get_pid();
342 float gembed = trutheval->
get_embed(g4particle);
343 float gprimary = trutheval->
is_primary(g4particle);
346 std::cout <<
PHWHERE <<
" PHG4Particle ID " << gtrackID <<
" gembed " << gembed <<
" gflavor " << gflavor <<
" gprimary " << gprimary << std::endl;
349 std::map<unsigned int, std::shared_ptr<TrkrCluster> > truth_clusters = trutheval->
all_truth_clusters(g4particle);
350 for (
auto it = truth_clusters.begin();
it != truth_clusters.end(); ++
it)
352 unsigned int layer =
it->first;
353 std::shared_ptr<TrkrCluster> gclus =
it->second;
354 const auto gkey = gclus->getClusKey();
357 if (layer < 7)
continue;
359 float gx = gclus->getX();
360 float gy = gclus->getY();
361 float gz = gclus->getZ();
362 float gedep = gclus->getError(0, 0);
363 float ng4hits = gclus->getAdc();
366 const auto gphi = std::atan2(gclus->getY(), gclus->getX());
370 std::cout <<
" gkey " << gkey <<
" detID " << detID <<
" tpcId " <<
TrkrDefs::tpcId <<
" layer " << layer <<
" truth clus " << gkey <<
" ng4hits " << ng4hits <<
" gr " << gr <<
" gx " << gx <<
" gy " << gy <<
" gz " << gz
371 <<
" gphi " << gphi <<
" gedep " << gedep << std::endl;
390 const auto z_cluster = global(2);
391 const auto phi_cluster = (float) std::atan2(global(1), global(0));
392 const auto phi_error = rclus->
getRPhiError() / r_cluster;
396 const auto dz = z_cluster - gz;
400 int region =
it->second;
404 std::cout <<
" Found match in layer " << layer <<
" region " << region <<
" for gtrackID " << gtrackID << std::endl;
405 std::cout <<
" x " << rclus->
getX() <<
" y " << rclus->
getY() <<
" z " << rclus->
getZ() << std::endl;
406 std::cout <<
" gx " << gclus->getX() <<
" gy " << gclus->getY() <<
" gz " << gclus->getZ() << std::endl;
407 std::cout <<
" drphi " << r_cluster * dphi <<
" rphi_error " << r_cluster * phi_error <<
" dz " <<
dz <<
" z_error " << z_error << std::endl;
410 const auto hiter = histograms.find(region);
411 if (hiter == histograms.end())
continue;
414 auto fill = [](TH1*
h,
float value) {
if( h ) h->Fill(
value ); };
415 fill(hiter->second.drphi, r_cluster * dphi);
416 fill(hiter->second.rphi_error, r_cluster * phi_error);
417 fill(hiter->second.phi_pulls, dphi / phi_error);
420 fill(hiter->second.dz,
dz);
421 fill(hiter->second.z_error, z_error);
422 fill(hiter->second.z_pulls,
dz / z_error);
427 fill(hiter->second.csize, std::distance(hit_range.first, hit_range.second));
431 for (
auto hit_iter = hit_range.first; hit_iter != hit_range.second; ++hit_iter)
433 const auto& hit_key = hit_iter->second;
438 fill(hiter->second.csize_phi, phibins.size());
439 fill(hiter->second.csize_z, zbins.size());
454 for (
auto iter = range.first; iter != range.second; ++iter)
457 const auto& hit_key = iter->second;
464 for (
auto truth_iter = g4hit_map.begin(); truth_iter != g4hit_map.end(); ++truth_iter)
467 const auto g4hit_key = truth_iter->second.second;
473 if (g4hit) out.insert(g4hit);