Back to home page

sPhenix code displayed by LXR

 
 

    


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   // print what is coming into the code
0151   //-----------------------------------
0152   
0153   printInputInfo(topNode);
0154   
0155   //---------------------------
0156   // fill the Evaluator NTuples
0157   //---------------------------
0158 
0159   fillOutputNtuples(topNode);
0160   
0161   //--------------------------------------------------
0162   // Print out the ancestry information for this event
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     // event information
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   // print out some useful stuff for debugging
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     // event information
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     // cluster wise information
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     // track-wise information
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   } // if verbosity
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   // fill the Vertex NTuple
0562   //-----------------------
0563 
0564   if (_ntp_vertex) {
0565     //cout << "Filling ntp_vertex " << endl;
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     cout << "vertex: " 
0611          << " ievent " << vertex_data[0]
0612          << " vx " << vertex_data[1]
0613          << " vy " << vertex_data[2]
0614          << " vz " << vertex_data[3]
0615          << endl;
0616     */
0617 
0618     _ntp_vertex->Fill(vertex_data);      
0619       }
0620     }
0621   }
0622   
0623   //-----------------------
0624   // fill the gpoint NTuple
0625   //-----------------------
0626   
0627   if (_ntp_gpoint) {
0628     //cout << "Filling ntp_gpoint " << endl;    
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   // fill the G4hit NTuple
0679   //---------------------
0680 
0681   if (_ntp_g4hit) {
0682     //cout << "Filling ntp_g4hit " << endl;
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       } //       if (g4particle)
0752       
0753       std::set<SvtxCluster*> clusters = clustereval->all_clusters_from(g4hit);  
0754       float nclusters = clusters.size();
0755 
0756       // best cluster reco'd
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   // fill the Hit NTuple
0831   //--------------------
0832 
0833   if (_ntp_hit) {
0834     //cout << "Filling ntp_hit " << endl;
0835     // need things off of the DST...
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       } //   if (g4particle){
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   // fill the Cluster NTuple
0967   //------------------------
0968 
0969   //cout << "check for ntp_cluster" << endl;
0970 
0971   if (_ntp_cluster && !_scan_for_embedded) {
0972     //cout << "Filling ntp_cluster 1 " << endl;
0973     // need things off of the DST...
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       }      //   if (g4particle){
1067     } //  if (g4hit) {
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     //cout << "Filling ntp_cluster 2 " << endl;
1119 
1120     // if only scanning embedded signals, loop over all the tracks from
1121     // embedded particles and report all of their clusters, including those
1122     // from other sources (noise hits on the embedded track)
1123     
1124     // need things off of the DST...
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         }      //   if (g4particle){
1229       } //  if (g4hit) {
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   // fill the Gtrack NTuple
1282   //------------------------
1283 
1284   // need things off of the DST...
1285 
1286   //cout << "check for ntp_gtrack" << endl;
1287 
1288   if (_ntp_gtrack) {
1289     //cout << "Filling ntp_gtrack " << endl;
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     cout << " ievent " << _ievent
1444          << " gtrackID " << gtrackID
1445          << " gflavor " << gflavor
1446          << " ng4hits " << ng4hits
1447          << endl;
1448     */
1449 
1450     _ntp_gtrack->Fill(gtrack_data);
1451 
1452       }      
1453     }
1454   }
1455   
1456   //------------------------
1457   // fill the Track NTuple
1458   //------------------------
1459 
1460 
1461 
1462   if (_ntp_track) {
1463     //cout << "Filling ntp_track " << endl;
1464 
1465     // need things off of the DST...
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     cout << "ievent " << _ievent
1650          << " trackID " << trackID
1651          << " nhits " << nhits
1652          << " px " << px
1653          << " py " << py
1654          << " pz " << pz
1655          << " gembed " << gembed
1656          << " gprimary " << gprimary 
1657          << endl;
1658     */
1659     _ntp_track->Fill(track_data);
1660       }
1661     }
1662   }
1663   
1664   return;
1665 }