File indexing completed on 2025-08-03 08:15:53
0001 #include <g4eval/SvtxEvalStack.h>
0002
0003 #include <phool/PHCompositeNode.h>
0004 #include <fun4all/Fun4AllReturnCodes.h>
0005 #include <phool/getClass.h>
0006
0007 #include <g4hough/SvtxVertexMap.h>
0008 #include <g4hough/SvtxVertex.h>
0009 #include <g4hough/SvtxTrackMap.h>
0010 #include <g4hough/SvtxTrack.h>
0011 #include <g4hough/SvtxClusterMap.h>
0012 #include <g4hough/SvtxCluster.h>
0013 #include <g4hough/SvtxHitMap.h>
0014 #include <g4hough/SvtxHit.h>
0015
0016 #include <g4main/PHG4Hit.h>
0017 #include <g4main/PHG4Particle.h>
0018 #include <g4main/PHG4VtxPoint.h>
0019 #include <g4main/PHG4TruthInfoContainer.h>
0020
0021 #include <g4detectors/PHG4Cell.h>
0022
0023 #include <TFile.h>
0024 #include <TNtuple.h>
0025
0026 #include <iostream>
0027 #include <set>
0028 #include <cmath>
0029 #include <cassert>
0030 #include "SvtxEvaluatorHaiwang.h"
0031
0032 using namespace std;
0033
0034 SvtxEvaluatorHaiwang::SvtxEvaluatorHaiwang(const string &name, const string &filename) :
0035 SubsysReco("SvtxEvaluatorHaiwang"),
0036 _ievent(0),
0037 _svtxevalstack(nullptr),
0038 _strict(false),
0039 _errors(0),
0040 _do_vertex_eval(true),
0041 _do_gpoint_eval(true),
0042 _do_g4hit_eval(true),
0043 _do_hit_eval(true),
0044 _do_cluster_eval(true),
0045 _do_gtrack_eval(true),
0046 _do_track_eval(true),
0047 _scan_for_embedded(false),
0048 _ntp_vertex(nullptr),
0049 _ntp_gpoint(nullptr),
0050 _ntp_g4hit(nullptr),
0051 _ntp_hit(nullptr),
0052 _ntp_cluster(nullptr),
0053 _ntp_gtrack(nullptr),
0054 _ntp_track(nullptr),
0055 _filename(filename),
0056 _tfile(nullptr) {
0057 verbosity = 0;
0058 }
0059
0060 int SvtxEvaluatorHaiwang::Init(PHCompositeNode *topNode) {
0061
0062 _ievent = 0;
0063
0064 _tfile = new TFile(_filename.c_str(), "RECREATE");
0065
0066 if (_do_vertex_eval) _ntp_vertex = new TNtuple("ntp_vertex","vertex => max truth",
0067 "event:vx:vy:vz:ntracks:"
0068 "gvx:gvy:gvz:gvt:gntracks:"
0069 "nfromtruth");
0070
0071 if (_do_gpoint_eval) _ntp_gpoint = new TNtuple("ntp_gpoint","g4point => best vertex",
0072 "event:gvx:gvy:gvz:gvt:gntracks:"
0073 "vx:vy:vz:ntracks:"
0074 "nfromtruth");
0075
0076 if (_do_g4hit_eval) _ntp_g4hit = new TNtuple("ntp_g4hit","g4hit => best svtxcluster",
0077 "event:g4hitID:gx:gy:gz:gt:gedep:"
0078 "glayer:gtrackID:gflavor:"
0079 "gpx:gpy:gpz:gvx:gvy:gvz:"
0080 "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
0081 "gembed:gprimary:nclusters:"
0082 "clusID:x:y:z:e:adc:layer:size:"
0083 "phisize:zsize:efromtruth");
0084
0085 if (_do_hit_eval) _ntp_hit = new TNtuple("ntp_hit","svtxhit => max truth",
0086 "event:hitID:e:adc:layer:"
0087 "cellID:ecell:"
0088 "g4hitID:gedep:gx:gy:gz:gt:"
0089 "gtrackID:gflavor:"
0090 "gpx:gpy:gpz:gvx:gvy:gvz:"
0091 "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
0092 "gembed:gprimary:efromtruth");
0093
0094 if (_do_cluster_eval) _ntp_cluster = new TNtuple("ntp_cluster","svtxcluster => max truth",
0095 "event:hitID:x:y:z:ex:ey:ez:ephi:"
0096 "e:adc:layer:size:phisize:"
0097 "zsize:trackID:g4hitID:gx:"
0098 "gy:gz:gt:gtrackID:gflavor:"
0099 "gpx:gpy:gpz:gvx:gvy:gvz:"
0100 "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
0101 "gembed:gprimary:efromtruth");
0102
0103 if (_do_gtrack_eval) _ntp_gtrack = new TNtuple("ntp_gtrack","g4particle => best svtxtrack",
0104 "event:gtrackID:gflavor:gnhits:"
0105 "gpx:gpy:gpz:"
0106 "gvx:gvy:gvz:gvt:"
0107 "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
0108 "gembed:gprimary:"
0109 "trackID:px:py:pz:charge:quality:chisq:ndf:nhits:layers:"
0110 "dca2d:dca2dsigma:pcax:pcay:pcaz:nfromtruth:layersfromtruth:"
0111 "nmaps:nintt:ntpc");
0112
0113 if (_do_track_eval) _ntp_track = new TNtuple("ntp_track","svtxtrack => max truth",
0114 "event:trackID:px:py:pz:charge:"
0115 "quality:chisq:ndf:nhits:layers:"
0116 "dca2d:dca2dsigma:pcax:pcay:pcaz:"
0117 "presdphi:presdeta:prese3x3:prese:"
0118 "cemcdphi:cemcdeta:cemce3x3:cemce:"
0119 "hcalindphi:hcalindeta:hcaline3x3:hcaline:"
0120 "hcaloutdphi:hcaloutdeta:hcaloute3x3:hcaloute:"
0121 "gtrackID:gflavor:gnhits:"
0122 "gpx:gpy:gpz:"
0123 "gvx:gvy:gvz:gvt:"
0124 "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
0125 "gembed:gprimary:nfromtruth:layersfromtruth:"
0126 "nmaps:nintt:ntpc");
0127
0128 return Fun4AllReturnCodes::EVENT_OK;
0129 }
0130
0131 int SvtxEvaluatorHaiwang::InitRun(PHCompositeNode *topNode) {
0132 return Fun4AllReturnCodes::EVENT_OK;
0133 }
0134
0135 int SvtxEvaluatorHaiwang::process_event(PHCompositeNode *topNode) {
0136
0137 if ((verbosity > 0)&&(_ievent%100==0)) {
0138 cout << "SvtxEvaluatorHaiwang::process_event - Event = " << _ievent << endl;
0139 }
0140
0141 if (!_svtxevalstack) {
0142 _svtxevalstack = new SvtxEvalStack(topNode);
0143 _svtxevalstack->set_strict(_strict);
0144 _svtxevalstack->set_verbosity(verbosity+1);
0145 } else {
0146 _svtxevalstack->next_event(topNode);
0147 }
0148
0149
0150
0151
0152
0153 printInputInfo(topNode);
0154
0155
0156
0157
0158
0159 fillOutputNtuples(topNode);
0160
0161
0162
0163
0164
0165 printOutputInfo(topNode);
0166
0167 ++_ievent;
0168 return Fun4AllReturnCodes::EVENT_OK;
0169 }
0170
0171 int SvtxEvaluatorHaiwang::End(PHCompositeNode *topNode) {
0172
0173 _tfile->cd();
0174
0175 if (_ntp_vertex) _ntp_vertex->Write();
0176 if (_ntp_gpoint) _ntp_gpoint->Write();
0177 if (_ntp_g4hit) _ntp_g4hit->Write();
0178 if (_ntp_hit) _ntp_hit->Write();
0179 if (_ntp_cluster) _ntp_cluster->Write();
0180 if (_ntp_gtrack) _ntp_gtrack->Write();
0181 if (_ntp_track) _ntp_track->Write();
0182
0183 _tfile->Close();
0184
0185 delete _tfile;
0186
0187 if (verbosity > 0) {
0188 cout << "========================= SvtxEvaluatorHaiwang::End() ============================" << endl;
0189 cout << " " << _ievent << " events of output written to: " << _filename << endl;
0190 cout << "===========================================================================" << endl;
0191 }
0192
0193 _errors += _svtxevalstack->get_errors();
0194
0195 if (verbosity > -1) {
0196 if ((_errors > 0)||(verbosity > 0)) {
0197 cout << "SvtxEvaluatorHaiwang::End() - Error Count: " << _errors << endl;
0198 }
0199 }
0200
0201 delete _svtxevalstack;
0202
0203 return Fun4AllReturnCodes::EVENT_OK;
0204 }
0205
0206 void SvtxEvaluatorHaiwang::printInputInfo(PHCompositeNode *topNode) {
0207
0208 if (verbosity > 1) cout << "SvtxEvaluatorHaiwang::printInputInfo() entered" << endl;
0209
0210 if (verbosity > 3) {
0211
0212
0213 cout << endl;
0214 cout << PHWHERE << " INPUT FOR EVENT " << _ievent << endl;
0215
0216 cout << endl;
0217 cout << "---PHG4HITS-------------" << endl;
0218 _svtxevalstack->get_truth_eval()->set_strict(_strict);
0219 std::set<PHG4Hit*> g4hits = _svtxevalstack->get_truth_eval()->all_truth_hits();
0220 unsigned int ig4hit = 0;
0221 for(std::set<PHG4Hit*>::iterator iter = g4hits.begin();
0222 iter != g4hits.end();
0223 ++iter) {
0224 PHG4Hit *g4hit = *iter;
0225 cout << ig4hit << " of " << g4hits.size();
0226 cout << ": PHG4Hit: " << endl;
0227 g4hit->identify();
0228 ++ig4hit;
0229 }
0230
0231 cout << "---SVTXCLUSTERS-------------" << endl;
0232 SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
0233 if (clustermap) {
0234 unsigned int icluster = 0;
0235 for (SvtxClusterMap::Iter iter = clustermap->begin();
0236 iter != clustermap->end();
0237 ++iter) {
0238 SvtxCluster* cluster = iter->second;
0239 cout << icluster << " of " << clustermap->size();
0240 cout << ": SvtxCluster: " << endl;
0241 cluster->identify();
0242 ++icluster;
0243 }
0244 }
0245
0246 cout << "---SVXTRACKS-------------" << endl;
0247 SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
0248 if (trackmap) {
0249 unsigned int itrack = 0;
0250 for (SvtxTrackMap::Iter iter = trackmap->begin();
0251 iter != trackmap->end();
0252 ++iter) {
0253 cout << itrack << " of " << trackmap->size();
0254 SvtxTrack *track = iter->second;
0255 cout << " : SvtxTrack:" << endl;
0256 track->identify();
0257 cout << endl;
0258 }
0259 }
0260
0261 cout << "---SVXVERTEXES-------------" << endl;
0262 SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
0263 if (vertexmap) {
0264 unsigned int ivertex = 0;
0265 for (SvtxVertexMap::Iter iter = vertexmap->begin();
0266 iter != vertexmap->end();
0267 ++iter) {
0268 cout << ivertex << " of " << vertexmap->size();
0269 SvtxVertex *vertex = iter->second;
0270 cout << " : SvtxVertex:" << endl;
0271 vertex->identify();
0272 cout << endl;
0273 }
0274 }
0275 }
0276
0277 return;
0278 }
0279
0280 void SvtxEvaluatorHaiwang::printOutputInfo(PHCompositeNode *topNode) {
0281
0282 if (verbosity > 1) cout << "SvtxEvaluatorHaiwang::printOutputInfo() entered" << endl;
0283
0284
0285
0286
0287
0288 if (verbosity > 0) {
0289
0290 SvtxTrackEval* trackeval = _svtxevalstack->get_track_eval();
0291 SvtxClusterEval* clustereval = _svtxevalstack->get_cluster_eval();
0292 SvtxTruthEval* trutheval = _svtxevalstack->get_truth_eval();
0293
0294
0295 cout << endl;
0296 cout << PHWHERE << " NEW OUTPUT FOR EVENT " << _ievent << endl;
0297 cout << endl;
0298
0299 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
0300
0301 PHG4VtxPoint *gvertex = truthinfo->GetPrimaryVtx( truthinfo->GetPrimaryVertexIndex() );
0302 float gvx = gvertex->get_x();
0303 float gvy = gvertex->get_y();
0304 float gvz = gvertex->get_z();
0305
0306 float vx = NAN;
0307 float vy = NAN;
0308 float vz = NAN;
0309
0310 SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
0311 if (vertexmap) {
0312 if (!vertexmap->empty()) {
0313 SvtxVertex* vertex = (vertexmap->begin()->second);
0314
0315 vx = vertex->get_x();
0316 vy = vertex->get_y();
0317 vz = vertex->get_z();
0318 }
0319 }
0320
0321 cout << "===Vertex Reconstruction=======================" << endl;
0322 cout << "vtrue = (" << gvx << "," << gvy << "," << gvz << ") => vreco = (" << vx << "," << vy << "," << vz << ")" << endl;
0323 cout << endl;
0324
0325 cout << "===Tracking Summary============================" << endl;
0326 unsigned int ng4hits[100] = {0};
0327 std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits();
0328 for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
0329 iter != g4hits.end();
0330 ++iter) {
0331 PHG4Hit *g4hit = *iter;
0332 ++ng4hits[g4hit->get_layer()];
0333 }
0334
0335 SvtxHitMap* hitmap = findNode::getClass<SvtxHitMap>(topNode,"SvtxHitMap");
0336 unsigned int nhits[100] = {0};
0337 if (hitmap) {
0338 for (SvtxHitMap::Iter iter = hitmap->begin();
0339 iter != hitmap->end();
0340 ++iter) {
0341 SvtxHit* hit = iter->second;
0342 ++nhits[hit->get_layer()];
0343 }
0344 }
0345
0346 SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
0347 unsigned int nclusters[100] = {0};
0348 if (clustermap) {
0349 for (SvtxClusterMap::Iter iter = clustermap->begin();
0350 iter != clustermap->end();
0351 ++iter) {
0352 SvtxCluster* cluster = iter->second;
0353 ++nclusters[cluster->get_layer()];
0354 }
0355 }
0356
0357 for (unsigned int ilayer = 0; ilayer < 100; ++ilayer) {
0358 cout << "layer " << ilayer << ": nG4hits = " << ng4hits[ilayer]
0359 << " => nHits = " << nhits[ilayer]
0360 << " => nClusters = " << nclusters[ilayer] << endl;
0361 }
0362
0363 SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
0364
0365 cout << "nGtracks = " << std::distance(truthinfo->GetPrimaryParticleRange().first,
0366 truthinfo->GetPrimaryParticleRange().second);
0367 cout << " => nTracks = ";
0368 if (trackmap) cout << trackmap->size() << endl;
0369 else cout << 0 << endl;
0370
0371
0372 if (verbosity > 1) {
0373
0374 for(std::set<PHG4Hit*>::iterator iter = g4hits.begin();
0375 iter != g4hits.end();
0376 ++iter) {
0377 PHG4Hit *g4hit = *iter;
0378
0379 cout << endl;
0380 cout << "===PHG4Hit===================================" << endl;
0381 cout << " PHG4Hit: "; g4hit->identify();
0382
0383 std::set<SvtxCluster*> clusters = clustereval->all_clusters_from(g4hit);
0384
0385 for (std::set<SvtxCluster*>::iterator jter = clusters.begin();
0386 jter != clusters.end();
0387 ++jter) {
0388 SvtxCluster *cluster = *jter;
0389 cout << "===Created-SvtxCluster================" << endl;
0390 cout << "SvtxCluster: "; cluster->identify();
0391 }
0392 }
0393
0394 PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
0395 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0396 iter != range.second;
0397 ++iter) {
0398
0399 PHG4Particle *particle = iter->second;
0400
0401
0402 cout << endl;
0403
0404 cout << "=== Gtrack ===================================================" << endl;
0405 cout << " PHG4Particle id = " << particle->get_track_id() << endl;
0406 particle->identify();
0407 cout << " ptrue = (";
0408 cout.width(5); cout << particle->get_px();
0409 cout << ",";
0410 cout.width(5); cout << particle->get_py();
0411 cout << ",";
0412 cout.width(5); cout << particle->get_pz();
0413 cout << ")" << endl;
0414
0415 cout << " vtrue = (";
0416 cout.width(5); cout << truthinfo->GetVtx(particle->get_vtx_id())->get_x();
0417 cout << ",";
0418 cout.width(5); cout << truthinfo->GetVtx(particle->get_vtx_id())->get_y();
0419 cout << ",";
0420 cout.width(5); cout << truthinfo->GetVtx(particle->get_vtx_id())->get_z();
0421 cout << ")" << endl;
0422
0423 cout << " pt = " << sqrt(pow(particle->get_px(),2)+pow(particle->get_py(),2)) << endl;
0424 cout << " phi = " << atan2(particle->get_py(),particle->get_px()) << endl;
0425 cout << " eta = " << asinh(particle->get_pz()/sqrt(pow(particle->get_px(),2)+pow(particle->get_py(),2))) << endl;
0426
0427 cout << " embed flag = " << truthinfo->isEmbeded(particle->get_track_id()) << endl;
0428
0429 cout << " ---Associated-PHG4Hits-----------------------------------------" << endl;
0430 std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(particle);
0431 for(std::set<PHG4Hit*>::iterator jter = g4hits.begin();
0432 jter != g4hits.end();
0433 ++jter) {
0434 PHG4Hit *g4hit = *jter;
0435
0436 float x = 0.5*(g4hit->get_x(0)+g4hit->get_x(1));
0437 float y = 0.5*(g4hit->get_y(0)+g4hit->get_y(1));
0438 float z = 0.5*(g4hit->get_z(0)+g4hit->get_z(1));
0439
0440 cout << " #" << g4hit->get_hit_id() << " xtrue = (";
0441 cout.width(5); cout << x;
0442 cout << ",";
0443 cout.width(5); cout << y;
0444 cout << ",";
0445 cout.width(5); cout << z;
0446 cout << ")";
0447
0448 std::set<SvtxCluster*> clusters = clustereval->all_clusters_from(g4hit);
0449 for (std::set<SvtxCluster*>::iterator kter = clusters.begin();
0450 kter != clusters.end();
0451 ++kter) {
0452
0453 SvtxCluster *cluster = *kter;
0454
0455 float x = cluster->get_x();
0456 float y = cluster->get_y();
0457 float z = cluster->get_z();
0458
0459 cout << " => #" << cluster->get_id();
0460 cout << " xreco = (";
0461 cout.width(5); cout << x;
0462 cout << ",";
0463 cout.width(5); cout << y;
0464 cout << ",";
0465 cout.width(5); cout << z;
0466 cout << ")";
0467 }
0468
0469 cout << endl;
0470 }
0471
0472 if (trackmap&&clustermap) {
0473
0474 std::set<SvtxTrack*> tracks = trackeval->all_tracks_from(particle);
0475 for (std::set<SvtxTrack*>::iterator jter = tracks.begin();
0476 jter != tracks.end();
0477 ++jter) {
0478
0479 SvtxTrack *track = *jter;
0480
0481 float px = track->get_px();
0482 float py = track->get_py();
0483 float pz = track->get_pz();
0484
0485 cout << "===Created-SvtxTrack==========================================" << endl;
0486 cout << " SvtxTrack id = " << track->get_id() << endl;
0487 cout << " preco = (";
0488 cout.width(5); cout << px;
0489 cout << ",";
0490 cout.width(5); cout << py;
0491 cout << ",";
0492 cout.width(5); cout << pz;
0493 cout << ")" << endl;
0494 cout << " quality = " << track->get_quality() << endl;
0495 cout << " nfromtruth = " << trackeval->get_nclusters_contribution(track,particle) << endl;
0496
0497 cout << " ---Associated-SvtxClusters-to-PHG4Hits-------------------------" << endl;
0498
0499 for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
0500 iter != track->end_clusters();
0501 ++iter) {
0502 unsigned int cluster_id = *iter;
0503 SvtxCluster* cluster = clustermap->get(cluster_id);
0504
0505 float x = cluster->get_x();
0506 float y = cluster->get_y();
0507 float z = cluster->get_z();
0508
0509 cout << " #" << cluster->get_id() << " xreco = (";
0510 cout.width(5); cout << x;
0511 cout << ",";
0512 cout.width(5); cout << y;
0513 cout << ",";
0514 cout.width(5); cout << z;
0515 cout << ") =>";
0516
0517 PHG4Hit* g4hit = clustereval->max_truth_hit_by_energy(cluster);
0518 if ((g4hit) && (g4hit->get_trkid() == particle->get_track_id())) {
0519
0520 x = 0.5*(g4hit->get_x(0)+g4hit->get_x(1));
0521 y = 0.5*(g4hit->get_y(0)+g4hit->get_y(1));
0522 z = 0.5*(g4hit->get_z(0)+g4hit->get_z(1));
0523
0524 cout << " #" << g4hit->get_hit_id()
0525 << " xtrue = (";
0526 cout.width(5); cout << x;
0527 cout << ",";
0528 cout.width(5); cout << y;
0529 cout << ",";
0530 cout.width(5); cout << z;
0531 cout << ") => Gtrack id = " << g4hit->get_trkid();
0532 } else {
0533 cout << " noise hit";
0534 }
0535 }
0536
0537 cout << endl;
0538 }
0539 }
0540 }
0541 }
0542
0543 cout << endl;
0544
0545 }
0546
0547 return;
0548 }
0549
0550 void SvtxEvaluatorHaiwang::fillOutputNtuples(PHCompositeNode *topNode) {
0551
0552 if (verbosity > 1) cout << "SvtxEvaluatorHaiwang::fillOutputNtuples() entered" << endl;
0553
0554 SvtxVertexEval* vertexeval = _svtxevalstack->get_vertex_eval();
0555 SvtxTrackEval* trackeval = _svtxevalstack->get_track_eval();
0556 SvtxClusterEval* clustereval = _svtxevalstack->get_cluster_eval();
0557 SvtxHitEval* hiteval = _svtxevalstack->get_hit_eval();
0558 SvtxTruthEval* trutheval = _svtxevalstack->get_truth_eval();
0559
0560
0561
0562
0563
0564 if (_ntp_vertex) {
0565
0566 SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
0567 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
0568 if (vertexmap && truthinfo) {
0569 for (SvtxVertexMap::Iter iter = vertexmap->begin();
0570 iter != vertexmap->end();
0571 ++iter) {
0572 SvtxVertex* vertex = iter->second;
0573 PHG4VtxPoint* point = vertexeval->max_truth_point_by_ntracks(vertex);
0574
0575 float vx = vertex->get_x();
0576 float vy = vertex->get_y();
0577 float vz = vertex->get_z();
0578 float ntracks = vertex->size_tracks();
0579
0580 float gvx = NAN;
0581 float gvy = NAN;
0582 float gvz = NAN;
0583 float gvt = NAN;
0584 float gntracks = truthinfo->GetNumPrimaryVertexParticles();
0585 float nfromtruth = NAN;
0586
0587 if (point) {
0588 gvx = point->get_x();
0589 gvy = point->get_y();
0590 gvz = point->get_z();
0591 gvt = point->get_t();
0592 gntracks = truthinfo->GetNumPrimaryVertexParticles();
0593 nfromtruth = vertexeval->get_ntracks_contribution(vertex,point);
0594 }
0595
0596 float vertex_data[11] = {(float) _ievent,
0597 vx,
0598 vy,
0599 vz,
0600 ntracks,
0601 gvx,
0602 gvy,
0603 gvz,
0604 gvt,
0605 gntracks,
0606 nfromtruth
0607 };
0608
0609
0610
0611
0612
0613
0614
0615
0616
0617
0618 _ntp_vertex->Fill(vertex_data);
0619 }
0620 }
0621 }
0622
0623
0624
0625
0626
0627 if (_ntp_gpoint) {
0628
0629 SvtxVertexMap* vertexmap = findNode::getClass<SvtxVertexMap>(topNode,"SvtxVertexMap");
0630 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
0631
0632 if (vertexmap && truthinfo) {
0633
0634 PHG4VtxPoint* point = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
0635
0636 if (point) {
0637
0638 SvtxVertex* vertex = vertexeval->best_vertex_from(point);
0639
0640 float gvx = point->get_x();
0641 float gvy = point->get_y();
0642 float gvz = point->get_z();
0643 float gvt = point->get_t();
0644 float gntracks = truthinfo->GetNumPrimaryVertexParticles();
0645 float vx = NAN;
0646 float vy = NAN;
0647 float vz = NAN;
0648 float ntracks = NAN;
0649 float nfromtruth = NAN;
0650
0651 if (vertex) {
0652 vx = vertex->get_x();
0653 vy = vertex->get_y();
0654 vz = vertex->get_z();
0655 ntracks = vertex->size_tracks();
0656 nfromtruth = vertexeval->get_ntracks_contribution(vertex,point);
0657 }
0658
0659 float gpoint_data[11] = {(float) _ievent,
0660 gvx,
0661 gvy,
0662 gvz,
0663 gvt,
0664 gntracks,
0665 vx,
0666 vy,
0667 vz,
0668 ntracks,
0669 nfromtruth
0670 };
0671
0672 _ntp_gpoint->Fill(gpoint_data);
0673 }
0674 }
0675 }
0676
0677
0678
0679
0680
0681 if (_ntp_g4hit) {
0682
0683 std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits();
0684 for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
0685 iter != g4hits.end();
0686 ++iter) {
0687
0688 PHG4Hit *g4hit = *iter;
0689 PHG4Particle *g4particle = trutheval->get_particle(g4hit);
0690
0691 float g4hitID = g4hit->get_hit_id();
0692 float gx = g4hit->get_avg_x();
0693 float gy = g4hit->get_avg_y();
0694 float gz = g4hit->get_avg_z();
0695 float gt = g4hit->get_avg_t();
0696 float gedep = g4hit->get_edep();
0697 float glayer = g4hit->get_layer();
0698
0699 float gtrackID = g4hit->get_trkid();
0700
0701 float gflavor = NAN;
0702 float gpx = NAN;
0703 float gpy = NAN;
0704 float gpz = NAN;
0705
0706 float gvx = NAN;
0707 float gvy = NAN;
0708 float gvz = NAN;
0709
0710 float gembed = NAN;
0711 float gprimary = NAN;
0712
0713 float gfpx = 0.;
0714 float gfpy = 0.;
0715 float gfpz = 0.;
0716 float gfx = 0.;
0717 float gfy = 0.;
0718 float gfz = 0.;
0719
0720 if (g4particle) {
0721
0722 if (_scan_for_embedded) {
0723 if (trutheval->get_embed(g4particle) <= 0) continue;
0724 }
0725
0726 gflavor = g4particle->get_pid();
0727 gpx = g4particle->get_px();
0728 gpy = g4particle->get_py();
0729 gpz = g4particle->get_pz();
0730
0731 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
0732
0733 if (vtx) {
0734 gvx = vtx->get_x();
0735 gvy = vtx->get_y();
0736 gvz = vtx->get_z();
0737 }
0738
0739 PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
0740 if (outerhit) {
0741 gfpx = outerhit->get_px(1);
0742 gfpy = outerhit->get_py(1);
0743 gfpz = outerhit->get_pz(1);
0744 gfx = outerhit->get_x(1);
0745 gfy = outerhit->get_y(1);
0746 gfz = outerhit->get_z(1);
0747 }
0748
0749 gembed = trutheval->get_embed(g4particle);
0750 gprimary = trutheval->is_primary(g4particle);
0751 }
0752
0753 std::set<SvtxCluster*> clusters = clustereval->all_clusters_from(g4hit);
0754 float nclusters = clusters.size();
0755
0756
0757 SvtxCluster* cluster = clustereval->best_cluster_from(g4hit);
0758
0759 float clusID = NAN;
0760 float x = NAN;
0761 float y = NAN;
0762 float z = NAN;
0763 float e = NAN;
0764 float adc = NAN;
0765 float layer = NAN;
0766 float size = NAN;
0767 float phisize = NAN;
0768 float zsize = NAN;
0769 float efromtruth = NAN;
0770
0771 if (cluster) {
0772 clusID = cluster->get_id();
0773 x = cluster->get_x();
0774 y = cluster->get_y();
0775 z = cluster->get_z();
0776 e = cluster->get_e();
0777 adc = cluster->get_adc();
0778 layer = cluster->get_layer();
0779 size = cluster->size_hits();
0780 phisize = cluster->get_phi_size();
0781 zsize = cluster->get_z_size();
0782 if (g4particle) {
0783 efromtruth = clustereval->get_energy_contribution(cluster,g4particle);
0784 }
0785 }
0786
0787 float g4hit_data[36] = {(float) _ievent,
0788 g4hitID,
0789 gx,
0790 gy,
0791 gz,
0792 gt,
0793 gedep,
0794 glayer,
0795 gtrackID,
0796 gflavor,
0797 gpx,
0798 gpy,
0799 gpz,
0800 gvx,
0801 gvy,
0802 gvz,
0803 gfpx,
0804 gfpy,
0805 gfpz,
0806 gfx,
0807 gfy,
0808 gfz,
0809 gembed,
0810 gprimary,
0811 nclusters,
0812 clusID,
0813 x,
0814 y,
0815 z,
0816 e,
0817 adc,
0818 layer,
0819 size,
0820 phisize,
0821 zsize,
0822 efromtruth
0823 };
0824
0825 _ntp_g4hit->Fill(g4hit_data);
0826 }
0827 }
0828
0829
0830
0831
0832
0833 if (_ntp_hit) {
0834
0835
0836 SvtxHitMap* hitmap = findNode::getClass<SvtxHitMap>(topNode,"SvtxHitMap");
0837 if (hitmap) {
0838
0839 for (SvtxHitMap::Iter iter = hitmap->begin();
0840 iter != hitmap->end();
0841 ++iter) {
0842
0843 SvtxHit* hit = iter->second;
0844 PHG4Hit* g4hit = hiteval->max_truth_hit_by_energy(hit);
0845 PHG4Cell* g4cell = hiteval->get_cell(hit);
0846 PHG4Particle* g4particle = trutheval->get_particle(g4hit);
0847
0848 float event = _ievent;
0849 float hitID = hit->get_id();
0850 float e = hit->get_e();
0851 float adc = hit->get_adc();
0852 float layer = hit->get_layer();
0853 float cellID = hit->get_cellid();
0854 float ecell = g4cell->get_edep();
0855
0856 float g4hitID = NAN;
0857 float gedep = NAN;
0858 float gx = NAN;
0859 float gy = NAN;
0860 float gz = NAN;
0861 float gt = NAN;
0862 float gtrackID = NAN;
0863 float gflavor = NAN;
0864 float gpx = NAN;
0865 float gpy = NAN;
0866 float gpz = NAN;
0867 float gvx = NAN;
0868 float gvy = NAN;
0869 float gvz = NAN;
0870 float gfpx = NAN;
0871 float gfpy = NAN;
0872 float gfpz = NAN;
0873 float gfx = NAN;
0874 float gfy = NAN;
0875 float gfz = NAN;
0876 float gembed = NAN;
0877 float gprimary = NAN;
0878
0879 float efromtruth = NAN;
0880
0881 if (g4hit) {
0882 g4hitID = g4hit->get_hit_id();
0883 gedep = g4hit->get_edep();
0884 gx = g4hit->get_avg_x();
0885 gy = g4hit->get_avg_y();
0886 gz = g4hit->get_avg_z();
0887 gt = g4hit->get_avg_t();
0888
0889 if (g4particle) {
0890
0891 if (_scan_for_embedded) {
0892 if (trutheval->get_embed(g4particle) <= 0) continue;
0893 }
0894
0895 gtrackID = g4particle->get_track_id();
0896 gflavor = g4particle->get_pid();
0897 gpx = g4particle->get_px();
0898 gpy = g4particle->get_py();
0899 gpz = g4particle->get_pz();
0900
0901 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
0902
0903 if (vtx) {
0904 gvx = vtx->get_x();
0905 gvy = vtx->get_y();
0906 gvz = vtx->get_z();
0907 }
0908
0909 PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
0910 if (outerhit) {
0911 gfpx = outerhit->get_px(1);
0912 gfpy = outerhit->get_py(1);
0913 gfpz = outerhit->get_pz(1);
0914 gfx = outerhit->get_x(1);
0915 gfy = outerhit->get_y(1);
0916 gfz = outerhit->get_z(1);
0917 }
0918 gembed = trutheval->get_embed(g4particle);
0919 gprimary = trutheval->is_primary(g4particle);
0920 }
0921 }
0922
0923 if (g4particle) {
0924 efromtruth = hiteval->get_energy_contribution(hit,g4particle);
0925 }
0926
0927 float hit_data[33] = {
0928 event,
0929 hitID,
0930 e,
0931 adc,
0932 layer,
0933 cellID,
0934 ecell,
0935 g4hitID,
0936 gedep,
0937 gx,
0938 gy,
0939 gz,
0940 gt,
0941 gtrackID,
0942 gflavor,
0943 gpx,
0944 gpy,
0945 gpz,
0946 gvx,
0947 gvy,
0948 gvz,
0949 gfpx,
0950 gfpy,
0951 gfpz,
0952 gfx,
0953 gfy,
0954 gfz,
0955 gembed,
0956 gprimary,
0957 efromtruth
0958 };
0959
0960 _ntp_hit->Fill(hit_data);
0961 }
0962 }
0963 }
0964
0965
0966
0967
0968
0969
0970
0971 if (_ntp_cluster && !_scan_for_embedded) {
0972
0973
0974 SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
0975 if (clustermap) {
0976
0977 for (SvtxClusterMap::Iter iter = clustermap->begin();
0978 iter != clustermap->end();
0979 ++iter) {
0980
0981 SvtxCluster* cluster = iter->second;
0982 SvtxTrack* track = trackeval->best_track_from(cluster);
0983
0984 PHG4Hit *g4hit = clustereval->max_truth_hit_by_energy(cluster);
0985 PHG4Particle *g4particle = trutheval->get_particle(g4hit);
0986
0987 float hitID = cluster->get_id();
0988 float x = cluster->get_x();
0989 float y = cluster->get_y();
0990 float z = cluster->get_z();
0991
0992 float ex = sqrt(cluster->get_error(0,0));
0993 float ey = sqrt(cluster->get_error(1,1));
0994 float ez = cluster->get_z_error();
0995
0996 float ephi = cluster->get_rphi_error();
0997
0998 float e = cluster->get_e();
0999 float adc = cluster->get_adc();
1000 float layer = cluster->get_layer();
1001 float size = cluster->size_hits();
1002 float phisize = cluster->get_phi_size();
1003 float zsize = cluster->get_z_size();
1004
1005 float trackID = NAN;
1006 if (track) trackID = track->get_id();
1007
1008 float g4hitID = NAN;
1009 float gx = NAN;
1010 float gy = NAN;
1011 float gz = NAN;
1012 float gt = NAN;
1013 float gtrackID = NAN;
1014 float gflavor = NAN;
1015 float gpx = NAN;
1016 float gpy = NAN;
1017 float gpz = NAN;
1018 float gvx = NAN;
1019 float gvy = NAN;
1020 float gvz = NAN;
1021 float gfpx = NAN;
1022 float gfpy = NAN;
1023 float gfpz = NAN;
1024 float gfx = NAN;
1025 float gfy = NAN;
1026 float gfz = NAN;
1027 float gembed = NAN;
1028 float gprimary = NAN;
1029
1030 float efromtruth = NAN;
1031
1032 if (g4hit) {
1033 g4hitID = g4hit->get_hit_id();
1034 gx = g4hit->get_avg_x();
1035 gy = g4hit->get_avg_y();
1036 gz = g4hit->get_avg_z();
1037 gt = g4hit->get_avg_t();
1038
1039 if (g4particle) {
1040
1041 gtrackID = g4particle->get_track_id();
1042 gflavor = g4particle->get_pid();
1043 gpx = g4particle->get_px();
1044 gpy = g4particle->get_py();
1045 gpz = g4particle->get_pz();
1046
1047 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1048 if (vtx) {
1049 gvx = vtx->get_x();
1050 gvy = vtx->get_y();
1051 gvz = vtx->get_z();
1052 }
1053
1054 PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
1055 if (outerhit) {
1056 gfpx = outerhit->get_px(1);
1057 gfpy = outerhit->get_py(1);
1058 gfpz = outerhit->get_pz(1);
1059 gfx = outerhit->get_x(1);
1060 gfy = outerhit->get_y(1);
1061 gfz = outerhit->get_z(1);
1062 }
1063
1064 gembed = trutheval->get_embed(g4particle);
1065 gprimary = trutheval->is_primary(g4particle);
1066 }
1067 }
1068
1069 if (g4particle){
1070 efromtruth = clustereval->get_energy_contribution(cluster,g4particle);
1071 }
1072
1073 float cluster_data[38] = {(float) _ievent,
1074 hitID,
1075 x,
1076 y,
1077 z,
1078 ex,
1079 ey,
1080 ez,
1081 ephi,
1082 e,
1083 adc,
1084 layer,
1085 size,
1086 phisize,
1087 zsize,
1088 trackID,
1089 g4hitID,
1090 gx,
1091 gy,
1092 gz,
1093 gt,
1094 gtrackID,
1095 gflavor,
1096 gpx,
1097 gpy,
1098 gpz,
1099 gvx,
1100 gvy,
1101 gvz,
1102 gfpx,
1103 gfpy,
1104 gfpz,
1105 gfx,
1106 gfy,
1107 gfz,
1108 gembed,
1109 gprimary,
1110 efromtruth};
1111
1112 _ntp_cluster->Fill(cluster_data);
1113 }
1114 }
1115
1116 } else if (_ntp_cluster && _scan_for_embedded) {
1117
1118
1119
1120
1121
1122
1123
1124
1125 SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
1126 SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
1127 if (trackmap) {
1128
1129 for (SvtxTrackMap::Iter iter = trackmap->begin();
1130 iter != trackmap->end();
1131 ++iter) {
1132
1133 SvtxTrack* track = iter->second;
1134 PHG4Particle* truth = trackeval->max_truth_particle_by_nclusters(track);
1135 if (truth) {
1136 if (trutheval->get_embed(truth) <= 0) continue;
1137 }
1138
1139 for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
1140 iter != track->end_clusters();
1141 ++iter) {
1142
1143 unsigned int cluster_id = *iter;
1144 SvtxCluster* cluster = clustermap->get(cluster_id);
1145
1146 PHG4Hit *g4hit = clustereval->max_truth_hit_by_energy(cluster);
1147 PHG4Particle *g4particle = trutheval->get_particle(g4hit);
1148
1149 float hitID = cluster->get_id();
1150 float x = cluster->get_x();
1151 float y = cluster->get_y();
1152 float z = cluster->get_z();
1153
1154 float ex = sqrt(cluster->get_error(0,0));
1155 float ey = sqrt(cluster->get_error(1,1));
1156 float ez = cluster->get_z_error();
1157
1158 float ephi = cluster->get_rphi_error();
1159
1160 float e = cluster->get_e();
1161 float adc = cluster->get_adc();
1162 float layer = cluster->get_layer();
1163 float size = cluster->size_hits();
1164 float phisize = cluster->get_phi_size();
1165 float zsize = cluster->get_z_size();
1166
1167 float trackID = NAN;
1168 if (track) trackID = track->get_id();
1169
1170 float g4hitID = NAN;
1171 float gx = NAN;
1172 float gy = NAN;
1173 float gz = NAN;
1174 float gt = NAN;
1175 float gtrackID = NAN;
1176 float gflavor = NAN;
1177 float gpx = NAN;
1178 float gpy = NAN;
1179 float gpz = NAN;
1180 float gvx = NAN;
1181 float gvy = NAN;
1182 float gvz = NAN;
1183 float gfpx = NAN;
1184 float gfpy = NAN;
1185 float gfpz = NAN;
1186 float gfx = NAN;
1187 float gfy = NAN;
1188 float gfz = NAN;
1189 float gembed = NAN;
1190 float gprimary = NAN;
1191
1192 float efromtruth = NAN;
1193
1194 if (g4hit) {
1195 g4hitID = g4hit->get_hit_id();
1196 gx = g4hit->get_avg_x();
1197 gy = g4hit->get_avg_y();
1198 gz = g4hit->get_avg_z();
1199 gt = g4hit->get_avg_t();
1200
1201 if (g4particle) {
1202
1203 gtrackID = g4particle->get_track_id();
1204 gflavor = g4particle->get_pid();
1205 gpx = g4particle->get_px();
1206 gpy = g4particle->get_py();
1207 gpz = g4particle->get_pz();
1208
1209 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1210 if (vtx) {
1211 gvx = vtx->get_x();
1212 gvy = vtx->get_y();
1213 gvz = vtx->get_z();
1214 }
1215
1216 PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
1217 if (outerhit) {
1218 gfpx = outerhit->get_px(1);
1219 gfpy = outerhit->get_py(1);
1220 gfpz = outerhit->get_pz(1);
1221 gfx = outerhit->get_x(1);
1222 gfy = outerhit->get_y(1);
1223 gfz = outerhit->get_z(1);
1224 }
1225
1226 gembed = trutheval->get_embed(g4particle);
1227 gprimary = trutheval->is_primary(g4particle);
1228 }
1229 }
1230
1231 if (g4particle){
1232 efromtruth = clustereval->get_energy_contribution(cluster,g4particle);
1233 }
1234
1235 float cluster_data[38] = {(float) _ievent,
1236 hitID,
1237 x,
1238 y,
1239 z,
1240 ex,
1241 ey,
1242 ez,
1243 ephi,
1244 e,
1245 adc,
1246 layer,
1247 size,
1248 phisize,
1249 zsize,
1250 trackID,
1251 g4hitID,
1252 gx,
1253 gy,
1254 gz,
1255 gt,
1256 gtrackID,
1257 gflavor,
1258 gpx,
1259 gpy,
1260 gpz,
1261 gvx,
1262 gvy,
1263 gvz,
1264 gfpx,
1265 gfpy,
1266 gfpz,
1267 gfx,
1268 gfy,
1269 gfz,
1270 gembed,
1271 gprimary,
1272 efromtruth};
1273
1274 _ntp_cluster->Fill(cluster_data);
1275 }
1276 }
1277 }
1278 }
1279
1280
1281
1282
1283
1284
1285
1286
1287
1288 if (_ntp_gtrack) {
1289
1290
1291 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode,"G4TruthInfo");
1292 SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
1293 if (truthinfo) {
1294
1295 PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
1296 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
1297 iter != range.second;
1298 ++iter) {
1299
1300 PHG4Particle* g4particle = iter->second;
1301
1302 if (_scan_for_embedded) {
1303 if (trutheval->get_embed(g4particle) <= 0) continue;
1304 }
1305
1306 float gtrackID = g4particle->get_track_id();
1307 float gflavor = g4particle->get_pid();
1308
1309 std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(g4particle);
1310 float ng4hits = g4hits.size();
1311 float gpx = g4particle->get_px();
1312 float gpy = g4particle->get_py();
1313 float gpz = g4particle->get_pz();
1314
1315 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1316 float gvx = vtx->get_x();
1317 float gvy = vtx->get_y();
1318 float gvz = vtx->get_z();
1319 float gvt = vtx->get_t();
1320
1321 float gfpx = 0.;
1322 float gfpy = 0.;
1323 float gfpz = 0.;
1324 float gfx = 0.;
1325 float gfy = 0.;
1326 float gfz = 0.;
1327
1328 PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
1329 if (outerhit) {
1330 gfpx = outerhit->get_px(1);
1331 gfpy = outerhit->get_py(1);
1332 gfpz = outerhit->get_pz(1);
1333 gfx = outerhit->get_x(1);
1334 gfy = outerhit->get_y(1);
1335 gfz = outerhit->get_z(1);
1336 }
1337
1338 float gembed = trutheval->get_embed(g4particle);
1339 float gprimary = trutheval->is_primary(g4particle);
1340
1341 SvtxTrack* track = trackeval->best_track_from(g4particle);
1342
1343 float trackID = NAN;
1344 float charge = NAN;
1345 float quality = NAN;
1346 float chisq = NAN;
1347 float ndf = NAN;
1348 float nhits = NAN;
1349 unsigned int layers = 0x0;
1350 float dca2d = NAN;
1351 float dca2dsigma = NAN;
1352 float px = NAN;
1353 float py = NAN;
1354 float pz = NAN;
1355 float pcax = NAN;
1356 float pcay = NAN;
1357 float pcaz = NAN;
1358
1359 float nfromtruth = NAN;
1360 float layersfromtruth = NAN;
1361
1362 float nmaps = 0;
1363 float nintt = 0;
1364 float ntpc = 0;
1365
1366 if (track) {
1367 trackID = track->get_id();
1368 charge = track->get_charge();
1369 quality = track->get_quality();
1370 chisq = track->get_chisq();
1371 ndf = track->get_ndf();
1372 nhits = track->size_clusters();
1373
1374 for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
1375 iter != track->end_clusters();
1376 ++iter) {
1377 unsigned int cluster_id = *iter;
1378 SvtxCluster* cluster = clustermap->get(cluster_id);
1379 unsigned int layer = cluster->get_layer();
1380 if (layer < 32) layers |= (0x1 << layer);
1381
1382
1383 if(layer >= 0 and layer <= 2) ++nmaps;
1384 if(layer >= 3 and layer <= 6) ++nintt;
1385 if(layer >= 7 and layer <= 66) ++ntpc;
1386 }
1387
1388 dca2d = track->get_dca2d();
1389 dca2dsigma = track->get_dca2d_error();
1390 px = track->get_px();
1391 py = track->get_py();
1392 pz = track->get_pz();
1393 pcax = track->get_x();
1394 pcay = track->get_y();
1395 pcaz = track->get_z();
1396
1397 nfromtruth = trackeval->get_nclusters_contribution(track,g4particle);
1398 layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track,g4particle);
1399 }
1400
1401 float gtrack_data[39] = {(float) _ievent,
1402 gtrackID,
1403 gflavor,
1404 ng4hits,
1405 gpx,
1406 gpy,
1407 gpz,
1408 gvx,
1409 gvy,
1410 gvz,
1411 gvt,
1412 gfpx,
1413 gfpy,
1414 gfpz,
1415 gfx,
1416 gfy,
1417 gfz,
1418 gembed,
1419 gprimary,
1420 trackID,
1421 px,
1422 py,
1423 pz,
1424 charge,
1425 quality,
1426 chisq,
1427 ndf,
1428 nhits,
1429 (float) layers,
1430 dca2d,
1431 dca2dsigma,
1432 pcax,
1433 pcay,
1434 pcaz,
1435 nfromtruth,
1436 layersfromtruth,
1437 nmaps,
1438 nintt,
1439 ntpc
1440 };
1441
1442
1443
1444
1445
1446
1447
1448
1449
1450 _ntp_gtrack->Fill(gtrack_data);
1451
1452 }
1453 }
1454 }
1455
1456
1457
1458
1459
1460
1461
1462 if (_ntp_track) {
1463
1464
1465
1466 SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode,"SvtxTrackMap");
1467 SvtxClusterMap* clustermap = findNode::getClass<SvtxClusterMap>(topNode,"SvtxClusterMap");
1468 if (trackmap) {
1469
1470 for (SvtxTrackMap::Iter iter = trackmap->begin();
1471 iter != trackmap->end();
1472 ++iter) {
1473
1474 SvtxTrack* track = iter->second;
1475
1476 float trackID = track->get_id();
1477 float charge = track->get_charge();
1478 float quality = track->get_quality();
1479 float chisq = track->get_chisq();
1480 float ndf = track->get_ndf();
1481 float nhits = track->size_clusters();
1482
1483 float nmaps = 0;
1484 float nintt = 0;
1485 float ntpc = 0;
1486
1487 unsigned int layers = 0x0;
1488 for (SvtxTrack::ConstClusterIter iter = track->begin_clusters();
1489 iter != track->end_clusters();
1490 ++iter) {
1491 unsigned int cluster_id = *iter;
1492 SvtxCluster* cluster = clustermap->get(cluster_id);
1493 unsigned int layer = cluster->get_layer();
1494 if (layer < 31) layers |= (0x1 << layer);
1495
1496 if(layer >= 0 and layer <= 2) ++nmaps;
1497 if(layer >= 3 and layer <= 6) ++nintt;
1498 if(layer >= 7 and layer <= 66) ++ntpc;
1499 }
1500
1501 float dca2d = track->get_dca2d();
1502 float dca2dsigma = track->get_dca2d_error();
1503 float px = track->get_px();
1504 float py = track->get_py();
1505 float pz = track->get_pz();
1506 float pcax = track->get_x();
1507 float pcay = track->get_y();
1508 float pcaz = track->get_z();
1509
1510 float presdphi = track->get_cal_dphi(SvtxTrack::PRES);
1511 float presdeta = track->get_cal_deta(SvtxTrack::PRES);
1512 float prese3x3 = track->get_cal_energy_3x3(SvtxTrack::PRES);
1513 float prese = track->get_cal_cluster_e(SvtxTrack::PRES);
1514
1515 float cemcdphi = track->get_cal_dphi(SvtxTrack::CEMC);
1516 float cemcdeta = track->get_cal_deta(SvtxTrack::CEMC);
1517 float cemce3x3 = track->get_cal_energy_3x3(SvtxTrack::CEMC);
1518 float cemce = track->get_cal_cluster_e(SvtxTrack::CEMC);
1519
1520 float hcalindphi = track->get_cal_dphi(SvtxTrack::HCALIN);
1521 float hcalindeta = track->get_cal_deta(SvtxTrack::HCALIN);
1522 float hcaline3x3 = track->get_cal_energy_3x3(SvtxTrack::HCALIN);
1523 float hcaline = track->get_cal_cluster_e(SvtxTrack::HCALIN);
1524
1525 float hcaloutdphi = track->get_cal_dphi(SvtxTrack::HCALOUT);
1526 float hcaloutdeta = track->get_cal_deta(SvtxTrack::HCALOUT);
1527 float hcaloute3x3 = track->get_cal_energy_3x3(SvtxTrack::HCALOUT);
1528 float hcaloute = track->get_cal_cluster_e(SvtxTrack::HCALOUT);
1529
1530 float gtrackID = NAN;
1531 float gflavor = NAN;
1532 float ng4hits = NAN;
1533 float gpx = NAN;
1534 float gpy = NAN;
1535 float gpz = NAN;
1536 float gvx = NAN;
1537 float gvy = NAN;
1538 float gvz = NAN;
1539 float gvt = NAN;
1540 float gfpx = NAN;
1541 float gfpy = NAN;
1542 float gfpz = NAN;
1543 float gfx = NAN;
1544 float gfy = NAN;
1545 float gfz = NAN;
1546 float gembed = NAN;
1547 float gprimary = NAN;
1548
1549 float nfromtruth = NAN;
1550 float layersfromtruth = NAN;
1551
1552 PHG4Particle* g4particle = trackeval->max_truth_particle_by_nclusters(track);
1553
1554 if (g4particle) {
1555
1556 if (_scan_for_embedded) {
1557 if (trutheval->get_embed(g4particle) <= 0) continue;
1558 }
1559
1560 gtrackID = g4particle->get_track_id();
1561 gflavor = g4particle->get_pid();
1562
1563 std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits(g4particle);
1564 ng4hits = g4hits.size();
1565 gpx = g4particle->get_px();
1566 gpy = g4particle->get_py();
1567 gpz = g4particle->get_pz();
1568
1569 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1570 gvx = vtx->get_x();
1571 gvy = vtx->get_y();
1572 gvz = vtx->get_z();
1573 gvt = vtx->get_t();
1574
1575 PHG4Hit* outerhit = trutheval->get_outermost_truth_hit(g4particle);
1576 if (outerhit) {
1577 gfpx = outerhit->get_px(1);
1578 gfpy = outerhit->get_py(1);
1579 gfpz = outerhit->get_pz(1);
1580 gfx = outerhit->get_x(1);
1581 gfy = outerhit->get_y(1);
1582 gfz = outerhit->get_z(1);
1583 }
1584 gembed = trutheval->get_embed(g4particle);
1585 gprimary = trutheval->is_primary(g4particle);
1586
1587 nfromtruth = trackeval->get_nclusters_contribution(track,g4particle);
1588 layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track,g4particle);
1589 }
1590
1591 float track_data[55] = {(float) _ievent,
1592 trackID,
1593 px,
1594 py,
1595 pz,
1596 charge,
1597 quality,
1598 chisq,
1599 ndf,
1600 nhits,
1601 (float) layers,
1602 dca2d,
1603 dca2dsigma,
1604 pcax,
1605 pcay,
1606 pcaz,
1607 presdphi,
1608 presdeta,
1609 prese3x3,
1610 prese,
1611 cemcdphi,
1612 cemcdeta,
1613 cemce3x3,
1614 cemce,
1615 hcalindphi,
1616 hcalindeta,
1617 hcaline3x3,
1618 hcaline,
1619 hcaloutdphi,
1620 hcaloutdeta,
1621 hcaloute3x3,
1622 hcaloute,
1623 gtrackID,
1624 gflavor,
1625 ng4hits,
1626 gpx,
1627 gpy,
1628 gpz,
1629 gvx,
1630 gvy,
1631 gvz,
1632 gvt,
1633 gfpx,
1634 gfpy,
1635 gfpz,
1636 gfx,
1637 gfy,
1638 gfz,
1639 gembed,
1640 gprimary,
1641 nfromtruth,
1642 layersfromtruth,
1643 nmaps,
1644 nintt,
1645 ntpc
1646 };
1647
1648
1649
1650
1651
1652
1653
1654
1655
1656
1657
1658
1659 _ntp_track->Fill(track_data);
1660 }
1661 }
1662 }
1663
1664 return;
1665 }