13 #include <g4vertex/GlobalVertex.h>
14 #include <g4vertex/GlobalVertexMap.h>
16 #include <calobase/RawCluster.h>
17 #include <calobase/RawClusterContainer.h>
18 #include <calobase/RawClusterUtility.h>
19 #include <calobase/RawTower.h>
20 #include <calobase/RawTowerContainer.h>
21 #include <calobase/RawTowerGeom.h>
22 #include <calobase/RawTowerGeomContainer.h>
34 #include <CLHEP/Vector/ThreeVector.h>
57 , _truth_trace_embed_flags()
58 , _truth_e_threshold(0.0)
60 _reco_e_threshold(0.0)
62 _caloevalstack(nullptr)
64 , _do_gpoint_eval(
false)
65 , _do_gshower_eval(
true)
66 , _do_tower_eval(
false)
67 , _do_cluster_eval(
true)
68 , _ntp_gpoint(nullptr)
69 , _ntp_gshower(nullptr)
71 , _tower_debug(nullptr)
72 , _ntp_cluster(nullptr)
89 "event:gparticleID:gflavor:gnhits:"
90 "geta:gphi:ge:gpt:gvx:gvy:gvz:gembed:gedep:"
91 "clusterID:ntowers:eta:x:y:z:phi:e:efromtruth");
96 _ntp_tower =
new TNtuple(
"ntp_tower",
"tower => max truth primary",
97 "event:towerID:ieta:iphi:eta:phi:e:x:y:z:"
98 "gparticleID:gflavor:gnhits:"
99 "geta:gphi:ge:gpt:gvx:gvy:gvz:"
104 _tower_debug =
new TTree(
"tower_debug",
"tower => max truth primary");
119 _ntp_cluster =
new TNtuple(
"ntp_cluster",
"cluster => max truth primary",
120 "event:clusterID:ntowers:eta_detector:eta:x:y:z:phi:e:ecore:"
121 "gparticleID:gflavor:gnhits:"
122 "geta:gphi:ge:gpt:gvx:gvy:gvz:"
196 cout <<
"========================= " <<
Name() <<
"::End() ============================" << endl;
197 cout <<
" " <<
_ievent <<
" events of output written to: " <<
_filename << endl;
198 cout <<
"===========================================================================" << endl;
208 if (
Verbosity() > 2) cout <<
"CaloEvaluator::printInputInfo() entered" << endl;
222 cerr <<
PHWHERE <<
" ERROR: Can't find G4TruthInfo" << endl;
226 cout <<
Name() <<
": PHG4TruthInfoContainer contents: " << endl;
230 truthiter != truthrange.second;
235 cout << truthiter->first <<
" => pid: " << particle->
get_pid()
236 <<
" pt: " << sqrt(pow(particle->
get_px(), 2) + pow(particle->
get_py(), 2)) << endl;
245 if (
Verbosity() > 2) cout <<
"CaloEvaluator::printOutputInfo() entered" << endl;
265 cerr <<
PHWHERE <<
" ERROR: Can't find G4TruthInfo" << endl;
270 GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode,
"GlobalVertexMap");
273 float gvx = gvertex->
get_x();
274 float gvy = gvertex->get_y();
275 float gvz = gvertex->get_z();
282 if (!vertexmap->
empty())
286 vx = vertex->
get_x();
287 vy = vertex->
get_y();
288 vz = vertex->
get_z();
292 cout <<
"vtrue = (" << gvx <<
"," << gvy <<
"," << gvz <<
") => vreco = (" << vx <<
"," << vy <<
"," << vz <<
")" << endl;
296 iter != range.second;
303 cout <<
"===Primary PHG4Particle=========================================" << endl;
304 cout <<
" particle id = " << primary->
get_track_id() << endl;
305 cout <<
" flavor = " << primary->
get_pid() << endl;
306 cout <<
" (px,py,pz,e) = (";
311 float ge = primary->
get_e();
326 float gpt = sqrt(gpx * gpx + gpy * gpy);
328 if (gpt != 0.0) geta = asinh(gpz / gpt);
329 float gphi = atan2(gpy, gpx);
331 cout <<
"(eta,phi,e,pt) = (";
350 cout <<
" vtrue = (";
361 cout <<
" embed = " << trutheval->
get_embed(primary) << endl;
365 for (std::set<RawCluster*>::iterator clusiter = clusters.begin();
366 clusiter != clusters.end();
372 float x = cluster->
get_x();
373 float y = cluster->
get_y();
374 float z = cluster->
get_z();
380 cout <<
" => #" << cluster->
get_id() <<
" (x,y,z,phi,e) = (";
395 cout <<
"), ntowers = " << ntowers <<
", efromtruth = " << efromtruth << endl;
406 if (
Verbosity() > 2) cout <<
"CaloEvaluator::fillOutputNtuples() entered" << endl;
422 cerr <<
PHWHERE <<
" ERROR: Can't find G4TruthInfo" << endl;
427 GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode,
"GlobalVertexMap");
430 float gvx = gvertex->
get_x();
431 float gvy = gvertex->get_y();
432 float gvz = gvertex->get_z();
439 if (!vertexmap->
empty())
443 vx = vertex->
get_x();
444 vy = vertex->
get_y();
445 vz = vertex->
get_z();
449 float gpoint_data[7] = {(float)
_ievent,
466 if (
Verbosity() > 1) cout <<
Name() <<
" CaloEvaluator::filling gshower ntuple..." << endl;
468 GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode,
"GlobalVertexMap");
473 cerr <<
PHWHERE <<
" ERROR: Can't find G4TruthInfo" << endl;
479 iter != range.second;
504 float ge = primary->
get_e();
506 float gpt = sqrt(gpx * gpx + gpy * gpy);
508 if (gpt != 0.0) geta = asinh(gpz / gpt);
509 float gphi = atan2(gpy, gpx);
512 float gvx = vtx->
get_x();
513 float gvy = vtx->
get_y();
514 float gvz = vtx->
get_z();
516 float gembed = trutheval->
get_embed(primary);
521 float clusterID = NAN;
534 clusterID = cluster->
get_id();
536 x = cluster->
get_x();
537 y = cluster->
get_y();
538 z = cluster->
get_z();
547 if (!vertexmap->
empty())
559 float shower_data[] = {(float)
_ievent,
592 if (
Verbosity() > 1) cout <<
"CaloEvaluator::filling tower ntuple..." << endl;
594 string towernode =
"TOWER_CALIB_" +
_caloname;
595 RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(topNode, towernode.c_str());
598 cerr <<
PHWHERE <<
" ERROR: Can't find " << towernode << endl;
602 string towergeomnode =
"TOWERGEOM_" +
_caloname;
603 RawTowerGeomContainer* towergeom = findNode::getClass<RawTowerGeomContainer>(topNode, towergeomnode.c_str());
606 cerr <<
PHWHERE <<
" ERROR: Can't find " << towergeomnode << endl;
612 for (rtiter = begin_end.first; rtiter != begin_end.second; ++rtiter)
621 cerr <<
PHWHERE <<
" ERROR: Can't find tower geometry for this tower hit: ";
650 float gparticleID = NAN;
684 ge = primary->
get_e();
686 gpt = sqrt(gpx * gpx + gpy * gpy);
687 if (gpt != 0.0) geta = asinh(gpz / gpt);
688 gphi = atan2(gpy, gpx);
705 float tower_data[] = {(float)
_ievent,
740 if (
Verbosity() > 1) cout <<
"CaloEvaluator::filling gcluster ntuple..." << endl;
742 GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode,
"GlobalVertexMap");
744 string clusternode =
"CLUSTER_" +
_caloname;
748 cerr <<
PHWHERE <<
" ERROR: Can't find " << clusternode << endl;
752 string towernode =
"TOWER_CALIB_" +
_caloname;
753 RawTowerContainer* towers = findNode::getClass<RawTowerContainer>(topNode, towernode.c_str());
755 cerr <<
PHWHERE <<
" ERROR: Can't find towers " << towernode << endl;
771 float clusterID = cluster->
get_id();
773 float x = cluster->
get_x();
774 float y = cluster->
get_y();
775 float z = cluster->
get_z();
785 if (!vertexmap->
empty())
795 float gparticleID = NAN;
830 ge = primary->
get_e();
832 gpt = sqrt(gpx * gpx + gpy * gpy);
833 if (gpt != 0.0) geta = asinh(gpz / gpt);
834 gphi = atan2(gpy, gpx);
852 float cluster_data[] = {(float)
_ievent,
885 for (toweriter = towersConstRange.first; toweriter != towersConstRange.second; ++toweriter) {