File indexing completed on 2025-08-05 08:17:59
0001 #include "SvtxEvaluator.h"
0002
0003 #include "SvtxEvalStack.h"
0004
0005 #include "SvtxClusterEval.h"
0006 #include "SvtxHitEval.h"
0007 #include "SvtxTrackEval.h"
0008 #include "SvtxTruthEval.h"
0009 #include "SvtxVertexEval.h"
0010
0011 #include <trackbase/ActsGeometry.h>
0012 #include <trackbase/ClusterErrorPara.h>
0013 #include <trackbase/TpcDefs.h>
0014 #include <trackbase/TrkrCluster.h>
0015 #include <trackbase/TrkrClusterContainer.h>
0016 #include <trackbase/TrkrClusterHitAssoc.h>
0017 #include <trackbase/TrkrClusterIterationMapv1.h>
0018 #include <trackbase/TrkrDefs.h>
0019 #include <trackbase/TrkrHit.h>
0020 #include <trackbase/TrkrHitSet.h>
0021 #include <trackbase/TrkrHitSetContainer.h>
0022
0023 #include <trackbase_historic/ActsTransformations.h>
0024 #include <trackbase_historic/SvtxTrack.h>
0025 #include <trackbase_historic/SvtxTrackMap.h>
0026 #include <trackbase_historic/TrackSeed.h>
0027 #include <trackbase_historic/TrackAnalysisUtils.h>
0028
0029 #include <globalvertex/SvtxVertex.h>
0030 #include <globalvertex/SvtxVertexMap.h>
0031 #include <globalvertex/GlobalVertex.h>
0032 #include <globalvertex/GlobalVertexMap.h>
0033
0034 #include <g4main/PHG4Hit.h>
0035 #include <g4main/PHG4Particle.h>
0036 #include <g4main/PHG4TruthInfoContainer.h>
0037 #include <g4main/PHG4VtxPoint.h>
0038
0039 #include <g4detectors/PHG4TpcCylinderGeom.h>
0040 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0041
0042 #include <fun4all/Fun4AllReturnCodes.h>
0043 #include <fun4all/SubsysReco.h>
0044
0045 #include <phool/PHTimer.h>
0046 #include <phool/getClass.h>
0047 #include <phool/phool.h>
0048 #include <phool/recoConsts.h>
0049
0050 #include <TFile.h>
0051 #include <TNtuple.h>
0052 #include <TVector3.h>
0053
0054 #include <cmath>
0055 #include <iomanip>
0056 #include <iostream>
0057 #include <iterator>
0058 #include <map>
0059 #include <memory> // for shared_ptr
0060 #include <set> // for _Rb_tree_cons...
0061 #include <utility>
0062 #include <vector>
0063
0064 SvtxEvaluator::SvtxEvaluator(const std::string& , const std::string& filename, const std::string& trackmapname,
0065 unsigned int nlayers_maps,
0066 unsigned int nlayers_intt,
0067 unsigned int nlayers_tpc,
0068 unsigned int nlayers_mms)
0069 : SubsysReco("SvtxEvaluator")
0070 , _nlayers_maps(nlayers_maps)
0071 , _nlayers_intt(nlayers_intt)
0072 , _nlayers_tpc(nlayers_tpc)
0073 , _nlayers_mms(nlayers_mms)
0074 , _filename(filename)
0075 , _trackmapname(trackmapname)
0076 {
0077 }
0078
0079 SvtxEvaluator::~SvtxEvaluator()
0080 {
0081 delete _timer;
0082 }
0083
0084 int SvtxEvaluator::Init(PHCompositeNode* )
0085 {
0086 _ievent = 0;
0087
0088 _tfile = new TFile(_filename.c_str(), "RECREATE");
0089 _tfile->SetCompressionLevel(0);
0090 if (_do_info_eval)
0091 {
0092 _ntp_info = new TNtuple("ntp_info", "event info",
0093 "event:seed:"
0094 "occ11:occ116:occ21:occ216:occ31:occ316:"
0095 "gntrkall:gntrkprim:ntrk:"
0096 "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
0097 }
0098
0099 if (_do_vertex_eval)
0100 {
0101 _ntp_vertex = new TNtuple("ntp_vertex", "vertex => max truth",
0102 "event:seed:vertexID:vx:vy:vz:ntracks:chi2:ndof:"
0103 "gvx:gvy:gvz:gvt:gembed:gntracks:gntracksmaps:"
0104 "gnembed:nfromtruth:"
0105 "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
0106 }
0107
0108 if (_do_gpoint_eval)
0109 {
0110 _ntp_gpoint = new TNtuple("ntp_gpoint", "g4point => best vertex",
0111 "event:seed:gvx:gvy:gvz:gvt:gntracks:gembed:"
0112 "vx:vy:vz:ntracks:"
0113 "nfromtruth:"
0114 "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
0115 }
0116
0117 if (_do_g4hit_eval)
0118 {
0119 _ntp_g4hit = new TNtuple("ntp_g4hit", "g4hit => best svtxcluster",
0120 "event:seed:g4hitID:gx:gy:gz:gt:gpl:gedep:geta:gphi:"
0121 "gdphi:gdz:"
0122 "glayer:gtrackID:gflavor:"
0123 "gpx:gpy:gpz:"
0124 "gvx:gvy:gvz:"
0125 "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
0126 "gembed:gprimary:nclusters:"
0127 "clusID:x:y:z:eta:phi:e:adc:layer:size:"
0128 "efromtruth:dphitru:detatru:dztru:drtru:"
0129 "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
0130 }
0131
0132 if (_do_hit_eval)
0133 {
0134 _ntp_hit = new TNtuple("ntp_hit", "svtxhit => max truth",
0135 "event:seed:hitID:e:adc:layer:phielem:zelem:"
0136 "cellID:ecell:phibin:zbin:phi:x:y:z:"
0137 "g4hitID:gedep:gx:gy:gz:gt:"
0138 "gtrackID:gflavor:"
0139 "gpx:gpy:gpz:gvx:gvy:gvz:gvt:"
0140 "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
0141 "gembed:gprimary:efromtruth:"
0142 "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
0143 }
0144
0145 if (_do_cluster_eval)
0146 {
0147 _ntp_cluster = new TNtuple("ntp_cluster", "svtxcluster => max truth",
0148 "event:seed:hitID:x:y:z:r:phi:eta:theta:ex:ey:ez:ephi:pez:pephi:"
0149 "e:adc:maxadc:layer:phielem:zelem:size:phisize:zsize:"
0150 "pedge:redge:ovlp:"
0151 "trackID:niter:g4hitID:gx:"
0152 "gy:gz:gr:gphi:geta:gt:gtrackID:gflavor:"
0153 "gpx:gpy:gpz:gvx:gvy:gvz:gvt:"
0154 "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
0155 "gembed:gprimary:efromtruth:nparticles:"
0156 "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
0157 }
0158
0159 if (_do_g4cluster_eval)
0160 {
0161 _ntp_g4cluster = new TNtuple("ntp_g4cluster", "g4cluster => max truth",
0162 "event:layer:gx:gy:gz:gt:gedep:gr:gphi:geta:gtrackID:gflavor:gembed:gprimary:gphisize:gzsize:gadc:nreco:x:y:z:r:phi:eta:ex:ey:ez:ephi:adc:phisize:zsize");
0163 }
0164
0165 if (_do_gtrack_eval)
0166 {
0167 _ntp_gtrack = new TNtuple("ntp_gtrack", "g4particle => best svtxtrack",
0168 "event:seed:gntracks:gnchghad:gtrackID:gflavor:gnhits:gnmaps:gnintt:gnmms:"
0169 "gnintt1:gnintt2:gnintt3:gnintt4:"
0170 "gnintt5:gnintt6:gnintt7:gnintt8:"
0171 "gntpc:gnlmaps:gnlintt:gnltpc:gnlmms:"
0172 "gpx:gpy:gpz:gpt:geta:gphi:"
0173 "gvx:gvy:gvz:gvt:"
0174 "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
0175 "gembed:gprimary:"
0176 "trackID:px:py:pz:pt:eta:phi:deltapt:deltaeta:deltaphi:"
0177 "siqr:siphi:sithe:six0:siy0:tpqr:tpphi:tpthe:tpx0:tpy0:"
0178 "charge:quality:chisq:ndf:nhits:layers:nmaps:nintt:ntpc:nmms:ntpc1:ntpc11:ntpc2:ntpc3:nlmaps:nlintt:nltpc:nlmms:"
0179 "vertexID:vx:vy:vz:dca2d:dca2dsigma:dca3dxy:dca3dxysigma:dca3dz:dca3dzsigma:pcax:pcay:pcaz:nfromtruth:nwrong:ntrumaps:nwrongmaps:ntruintt:nwrongintt:ntrutpc:nwrongtpc:ntrumms:nwrongmms:ntrutpc1:nwrongtpc1:ntrutpc11:nwrongtpc11:ntrutpc2:nwrongtpc2:ntrutpc3:nwrongtpc3:layersfromtruth:"
0180 "npedge:nredge:nbig:novlp:merr:msize:"
0181 "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
0182 }
0183
0184 if (_do_track_eval)
0185 {
0186 _ntp_track = new TNtuple("ntp_track", "svtxtrack => max truth",
0187 "event:seed:trackID:crossing:px:py:pz:pt:eta:phi:deltapt:deltaeta:deltaphi:"
0188 "siqr:siphi:sithe:six0:siy0:tpqr:tpphi:tpthe:tpx0:tpy0:"
0189 "charge:quality:chisq:ndf:nhits:nmaps:nintt:ntpc:nmms:ntpc1:ntpc11:ntpc2:ntpc3:nlmaps:nlintt:nltpc:nlmms:layers:"
0190 "vertexID:vx:vy:vz:dca2d:dca2dsigma:dca3dxy:dca3dxysigma:dca3dz:dca3dzsigma:pcax:pcay:pcaz:"
0191 "gtrackID:singlematch:gflavor:gnhits:gnmaps:gnintt:gntpc:gnmms:gnlmaps:gnlintt:gnltpc:gnlmms:"
0192 "gpx:gpy:gpz:gpt:geta:gphi:"
0193 "gvx:gvy:gvz:gvt:"
0194 "gfpx:gfpy:gfpz:gfx:gfy:gfz:"
0195 "gembed:gprimary:nfromtruth:nwrong:ntrumaps:nwrongmaps:ntruintt:nwrongintt:"
0196 "ntrutpc:nwrongtpc:ntrumms:nwrongmms:ntrutpc1:nwrongtpc1:ntrutpc11:nwrongtpc11:ntrutpc2:nwrongtpc2:ntrutpc3:nwrongtpc3:layersfromtruth:"
0197 "npedge:nredge:nbig:novlp:merr:msize:"
0198 "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
0199 }
0200
0201 if (_do_gseed_eval)
0202 {
0203 _ntp_gseed = new TNtuple("ntp_gseed", "seeds from truth",
0204 "event:seed:ntrk:gx:gy:gz:gr:geta:gphi:"
0205 "glayer:"
0206 "gpx:gpy:gpz:gtpt:gtphi:gteta:"
0207 "gvx:gvy:gvz:"
0208 "gembed:gprimary:gflav:"
0209 "dphiprev:detaprev:"
0210 "nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclusintt:nclusmaps:nclusmms");
0211 }
0212
0213 _timer = new PHTimer("_eval_timer");
0214 _timer->stop();
0215
0216 return Fun4AllReturnCodes::EVENT_OK;
0217 }
0218
0219 int SvtxEvaluator::InitRun(PHCompositeNode* )
0220 {
0221 return Fun4AllReturnCodes::EVENT_OK;
0222 }
0223
0224 int SvtxEvaluator::process_event(PHCompositeNode* topNode)
0225 {
0226 if ((Verbosity() > 1) && (_ievent % 100 == 0))
0227 {
0228 std::cout << "SvtxEvaluator::process_event - Event = " << _ievent << std::endl;
0229 }
0230
0231 recoConsts* rc = recoConsts::instance();
0232 if (rc->FlagExist("RANDOMSEED"))
0233 {
0234 _iseed = rc->get_IntFlag("RANDOMSEED");
0235 m_fSeed = _iseed;
0236 }
0237 else
0238 {
0239 _iseed = 0;
0240 m_fSeed = NAN;
0241 }
0242
0243 if (Verbosity() > 1)
0244 {
0245 std::cout << "SvtxEvaluator::process_event - Seed = " << _iseed << std::endl;
0246 }
0247
0248 if (!_svtxevalstack)
0249 {
0250 _svtxevalstack = new SvtxEvalStack(topNode);
0251 _svtxevalstack->set_strict(_strict);
0252 _svtxevalstack->set_verbosity(Verbosity());
0253 _svtxevalstack->set_use_initial_vertex(_use_initial_vertex);
0254 _svtxevalstack->set_use_genfit_vertex(_use_genfit_vertex);
0255 _svtxevalstack->next_event(topNode);
0256 }
0257 else
0258 {
0259 _svtxevalstack->next_event(topNode);
0260 }
0261
0262
0263
0264
0265
0266 printInputInfo(topNode);
0267
0268
0269
0270
0271
0272 fillOutputNtuples(topNode);
0273
0274
0275
0276
0277
0278
0279
0280 ++_ievent;
0281 return Fun4AllReturnCodes::EVENT_OK;
0282 }
0283
0284 int SvtxEvaluator::End(PHCompositeNode* )
0285 {
0286 _tfile->cd();
0287
0288 if (_ntp_info)
0289 {
0290 _ntp_info->Write();
0291 }
0292 if (_ntp_vertex)
0293 {
0294 _ntp_vertex->Write();
0295 }
0296 if (_ntp_gpoint)
0297 {
0298 _ntp_gpoint->Write();
0299 }
0300 if (_ntp_g4hit)
0301 {
0302 _ntp_g4hit->Write();
0303 }
0304 if (_ntp_hit)
0305 {
0306 _ntp_hit->Write();
0307 }
0308 if (_ntp_cluster)
0309 {
0310 _ntp_cluster->Write();
0311 }
0312 if (_ntp_g4cluster)
0313 {
0314 _ntp_g4cluster->Write();
0315 }
0316 if (_ntp_gtrack)
0317 {
0318 _ntp_gtrack->Write();
0319 }
0320 if (_ntp_track)
0321 {
0322 _ntp_track->Write();
0323 }
0324 if (_ntp_gseed)
0325 {
0326 _ntp_gseed->Write();
0327 }
0328
0329 _tfile->Close();
0330
0331 delete _tfile;
0332
0333 if (Verbosity() > 1)
0334 {
0335 std::cout << "========================= SvtxEvaluator::End() ============================" << std::endl;
0336 std::cout << " " << _ievent << " events of output written to: " << _filename << std::endl;
0337 std::cout << "===========================================================================" << std::endl;
0338 }
0339
0340 if (_svtxevalstack)
0341 {
0342 _errors += _svtxevalstack->get_errors();
0343
0344 if (Verbosity() > 1)
0345 {
0346 if ((_errors > 0) || (Verbosity() > 1))
0347 {
0348 std::cout << "SvtxEvaluator::End() - Error Count: " << _errors << std::endl;
0349 }
0350 }
0351
0352 delete _svtxevalstack;
0353 }
0354
0355 return Fun4AllReturnCodes::EVENT_OK;
0356 }
0357
0358 void SvtxEvaluator::printInputInfo(PHCompositeNode* topNode)
0359 {
0360 if (Verbosity() > 1)
0361 {
0362 std::cout << "SvtxEvaluator::printInputInfo() entered" << std::endl;
0363 }
0364
0365 if (Verbosity() > 3)
0366 {
0367
0368 std::cout << std::endl;
0369 std::cout << PHWHERE << " INPUT FOR EVENT " << _ievent << std::endl;
0370
0371 std::cout << std::endl;
0372 std::cout << "---PHG4HITS-------------" << std::endl;
0373 _svtxevalstack->get_truth_eval()->set_strict(_strict);
0374 std::set<PHG4Hit*> g4hits = _svtxevalstack->get_truth_eval()->all_truth_hits();
0375 unsigned int ig4hit = 0;
0376 for (std::set<PHG4Hit*>::iterator iter = g4hits.begin();
0377 iter != g4hits.end();
0378 ++iter)
0379 {
0380 PHG4Hit* g4hit = *iter;
0381 std::cout << ig4hit << " of " << g4hits.size();
0382 std::cout << ": PHG4Hit: " << std::endl;
0383 g4hit->identify();
0384 ++ig4hit;
0385 }
0386
0387 std::cout << "---SVTXCLUSTERS-------------" << std::endl;
0388 TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
0389 if (!clustermap)
0390 {
0391 clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0392 }
0393
0394 if (clustermap != nullptr)
0395 {
0396 unsigned int icluster = 0;
0397 for (const auto& hitsetkey : clustermap->getHitSetKeys())
0398 {
0399 auto range = clustermap->getClusters(hitsetkey);
0400 for (auto iter = range.first; iter != range.second; ++iter)
0401 {
0402 TrkrDefs::cluskey cluster_key = iter->first;
0403 std::cout << icluster << " with key " << cluster_key << " of " << clustermap->size();
0404 std::cout << ": SvtxCluster: " << std::endl;
0405 iter->second->identify();
0406 ++icluster;
0407 }
0408 }
0409 }
0410
0411 std::cout << "---SVXTRACKS-------------" << std::endl;
0412 SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname.c_str());
0413 if (trackmap)
0414 {
0415 unsigned int itrack = 0;
0416 for (SvtxTrackMap::Iter iter = trackmap->begin();
0417 iter != trackmap->end();
0418 ++iter)
0419 {
0420 std::cout << itrack << " of " << trackmap->size();
0421 SvtxTrack* track = iter->second;
0422 std::cout << " : SvtxTrack:" << std::endl;
0423 track->identify();
0424 std::cout << std::endl;
0425 ++itrack;
0426 }
0427 }
0428
0429 std::cout << "---SVXVERTEXES-------------" << std::endl;
0430 SvtxVertexMap* vertexmap = nullptr;
0431 if (_use_initial_vertex)
0432 {
0433 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0434 }
0435 else if (_use_genfit_vertex)
0436 {
0437 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapRefit");
0438 }
0439 else
0440 {
0441 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs");
0442 }
0443
0444 if (vertexmap)
0445 {
0446 unsigned int ivertex = 0;
0447 for (SvtxVertexMap::Iter iter = vertexmap->begin();
0448 iter != vertexmap->end();
0449 ++iter)
0450 {
0451 std::cout << ivertex << " of " << vertexmap->size();
0452 SvtxVertex* vertex = iter->second;
0453 std::cout << " : SvtxVertex:" << std::endl;
0454 vertex->identify();
0455 std::cout << std::endl;
0456 }
0457 }
0458 }
0459
0460 return;
0461 }
0462
0463 void SvtxEvaluator::printOutputInfo(PHCompositeNode* topNode)
0464 {
0465 if (Verbosity() > 1)
0466 {
0467 std::cout << "SvtxEvaluator::printOutputInfo() entered" << std::endl;
0468 }
0469
0470
0471
0472
0473
0474 if (Verbosity() > 100)
0475 {
0476 SvtxTrackEval* trackeval = _svtxevalstack->get_track_eval();
0477 SvtxClusterEval* clustereval = _svtxevalstack->get_cluster_eval();
0478 SvtxTruthEval* trutheval = _svtxevalstack->get_truth_eval();
0479
0480
0481 std::cout << std::endl;
0482 std::cout << PHWHERE << " NEW OUTPUT FOR EVENT " << _ievent << std::endl;
0483 std::cout << std::endl;
0484
0485 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
0486
0487 PHG4VtxPoint* gvertex = truthinfo->GetPrimaryVtx(truthinfo->GetPrimaryVertexIndex());
0488 float gvx = gvertex->get_x();
0489 float gvy = gvertex->get_y();
0490 float gvz = gvertex->get_z();
0491
0492 float vx = NAN;
0493 float vy = NAN;
0494 float vz = NAN;
0495
0496 SvtxVertexMap* vertexmap = nullptr;
0497 if (_use_initial_vertex)
0498 {
0499 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
0500 }
0501 else if (_use_genfit_vertex)
0502 {
0503 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapRefit");
0504 }
0505 else
0506 {
0507 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs");
0508 }
0509
0510 if (vertexmap)
0511 {
0512 if (!vertexmap->empty())
0513 {
0514 SvtxVertex* vertex = (vertexmap->begin()->second);
0515
0516 vx = vertex->get_x();
0517 vy = vertex->get_y();
0518 vz = vertex->get_z();
0519 }
0520 }
0521
0522 std::cout << "===Vertex Reconstruction=======================" << std::endl;
0523 std::cout << "vtrue = (" << gvx << "," << gvy << "," << gvz << ") => vreco = (" << vx << "," << vy << "," << vz << ")" << std::endl;
0524 std::cout << std::endl;
0525
0526 std::cout << "===Tracking Summary============================" << std::endl;
0527 unsigned int ng4hits[100] = {0};
0528 std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits();
0529 for (auto g4hit : g4hits)
0530 {
0531 ++ng4hits[g4hit->get_layer()];
0532 }
0533
0534 TrkrHitSetContainer* hitsetmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0535
0536 TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
0537 if (!clustermap)
0538 {
0539 clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0540 }
0541
0542 unsigned int nclusters[100] = {0};
0543 unsigned int nhits[100] = {0};
0544
0545 ActsGeometry* tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0546
0547 if (!tgeometry)
0548 {
0549 std::cout << PHWHERE << "No Acts geometry on node tree. Can't continue."
0550 << std::endl;
0551 }
0552
0553 if (hitsetmap)
0554 {
0555 TrkrHitSetContainer::ConstRange all_hitsets = hitsetmap->getHitSets();
0556 for (TrkrHitSetContainer::ConstIterator hitsetiter = all_hitsets.first;
0557 hitsetiter != all_hitsets.second;
0558 ++hitsetiter)
0559 {
0560
0561 unsigned int layer = TrkrDefs::getLayer(hitsetiter->first);
0562
0563
0564 TrkrHitSet::ConstRange hitrangei = hitsetiter->second->getHits();
0565 for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0566 hitr != hitrangei.second;
0567 ++hitr)
0568 {
0569 ++nhits[layer];
0570 }
0571 auto range = clustermap->getClusters(hitsetiter->first);
0572 for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0573 {
0574 const auto cluskey = clusIter->first;
0575
0576 unsigned int local_layer = TrkrDefs::getLayer(cluskey);
0577 nclusters[local_layer]++;
0578 }
0579 }
0580 }
0581
0582 PHG4TpcCylinderGeomContainer* geom_container =
0583 findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0584 {
0585 if (!geom_container)
0586 {
0587 std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0588 }
0589 return;
0590 }
0591
0592 for (unsigned int ilayer = 0; ilayer < _nlayers_maps + _nlayers_intt + _nlayers_tpc; ++ilayer)
0593 {
0594 PHG4TpcCylinderGeom* GeoLayer = geom_container->GetLayerCellGeom(ilayer);
0595
0596 std::cout << "layer " << ilayer << ": nG4hits = " << ng4hits[ilayer]
0597 << " => nHits = " << nhits[ilayer]
0598 << " => nClusters = " << nclusters[ilayer] << std::endl;
0599 if (ilayer >= _nlayers_maps + _nlayers_intt && ilayer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
0600 {
0601 std::cout << "layer " << ilayer
0602 << " => nphi = " << GeoLayer->get_phibins()
0603 << " => nz = " << GeoLayer->get_zbins()
0604 << " => ntot = " << GeoLayer->get_phibins() * GeoLayer->get_zbins()
0605 << std::endl;
0606 }
0607 }
0608
0609 SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname.c_str());
0610
0611 std::cout << "nGtracks = " << std::distance(truthinfo->GetPrimaryParticleRange().first, truthinfo->GetPrimaryParticleRange().second);
0612 std::cout << " => nTracks = ";
0613 if (trackmap)
0614 {
0615 std::cout << trackmap->size() << std::endl;
0616 }
0617 else
0618 {
0619 std::cout << 0 << std::endl;
0620 }
0621
0622
0623 if (Verbosity() > 1)
0624 {
0625 for (auto g4hit : g4hits)
0626 {
0627 std::cout << std::endl;
0628 std::cout << "===PHG4Hit===================================" << std::endl;
0629 std::cout << " PHG4Hit: ";
0630 g4hit->identify();
0631
0632 std::set<TrkrDefs::cluskey> clusters = clustereval->all_clusters_from(g4hit);
0633
0634 for (unsigned long cluster_key : clusters)
0635 {
0636 std::cout << "===Created-SvtxCluster================" << std::endl;
0637 std::cout << "SvtxCluster: ";
0638 TrkrCluster* cluster = clustermap->findCluster(cluster_key);
0639 cluster->identify();
0640 }
0641 }
0642
0643 PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
0644 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
0645 iter != range.second;
0646 ++iter)
0647 {
0648 PHG4Particle* particle = iter->second;
0649
0650
0651 std::cout << std::endl;
0652
0653 std::cout << "=== Gtrack ===================================================" << std::endl;
0654 std::cout << " PHG4Particle id = " << particle->get_track_id() << std::endl;
0655 particle->identify();
0656 std::cout << " ptrue = (";
0657 std::cout.width(5);
0658 std::cout << particle->get_px();
0659 std::cout << ",";
0660 std::cout.width(5);
0661 std::cout << particle->get_py();
0662 std::cout << ",";
0663 std::cout.width(5);
0664 std::cout << particle->get_pz();
0665 std::cout << ")" << std::endl;
0666
0667 std::cout << " vtrue = (";
0668 std::cout.width(5);
0669 std::cout << truthinfo->GetVtx(particle->get_vtx_id())->get_x();
0670 std::cout << ",";
0671 std::cout.width(5);
0672 std::cout << truthinfo->GetVtx(particle->get_vtx_id())->get_y();
0673 std::cout << ",";
0674 std::cout.width(5);
0675 std::cout << truthinfo->GetVtx(particle->get_vtx_id())->get_z();
0676 std::cout << ")" << std::endl;
0677
0678 std::cout << " pt = " << sqrt(pow(particle->get_px(), 2) + pow(particle->get_py(), 2)) << std::endl;
0679 std::cout << " phi = " << atan2(particle->get_py(), particle->get_px()) << std::endl;
0680 std::cout << " eta = " << asinh(particle->get_pz() / sqrt(pow(particle->get_px(), 2) + pow(particle->get_py(), 2))) << std::endl;
0681
0682 std::cout << " embed flag = " << truthinfo->isEmbeded(particle->get_track_id()) << std::endl;
0683
0684 std::cout << " ---Associated-PHG4Hits-----------------------------------------" << std::endl;
0685 std::set<PHG4Hit*> local_g4hits = trutheval->all_truth_hits(particle);
0686 for (auto g4hit : local_g4hits)
0687 {
0688 float x = 0.5 * (g4hit->get_x(0) + g4hit->get_x(1));
0689 float y = 0.5 * (g4hit->get_y(0) + g4hit->get_y(1));
0690 float z = 0.5 * (g4hit->get_z(0) + g4hit->get_z(1));
0691
0692 std::cout << " #" << g4hit->get_hit_id() << " xtrue = (";
0693 std::cout.width(5);
0694 std::cout << x;
0695 std::cout << ",";
0696 std::cout.width(5);
0697 std::cout << y;
0698 std::cout << ",";
0699 std::cout.width(5);
0700 std::cout << z;
0701 std::cout << ")";
0702
0703 std::set<TrkrDefs::cluskey> clusters = clustereval->all_clusters_from(g4hit);
0704 for (unsigned long cluster_key : clusters)
0705 {
0706 TrkrCluster* cluster = clustermap->findCluster(cluster_key);
0707
0708 Acts::Vector3 glob;
0709 glob = tgeometry->getGlobalPosition(cluster_key, cluster);
0710
0711 float local_x = glob(0);
0712 float local_y = glob(1);
0713 float local_z = glob(2);
0714
0715 std::cout << " => #" << cluster_key;
0716 std::cout << " xreco = (";
0717 std::cout.width(5);
0718 std::cout << local_x;
0719 std::cout << ",";
0720 std::cout.width(5);
0721 std::cout << local_y;
0722 std::cout << ",";
0723 std::cout.width(5);
0724 std::cout << local_z;
0725 std::cout << ")";
0726 }
0727
0728 std::cout << std::endl;
0729 }
0730
0731 if (trackmap && clustermap)
0732 {
0733 std::set<SvtxTrack*> tracks = trackeval->all_tracks_from(particle);
0734 for (auto track : tracks)
0735 {
0736 float px = track->get_px();
0737 float py = track->get_py();
0738 float pz = track->get_pz();
0739
0740 std::cout << "===Created-SvtxTrack==========================================" << std::endl;
0741 std::cout << " SvtxTrack id = " << track->get_id() << std::endl;
0742 std::cout << " preco = (";
0743 std::cout.width(5);
0744 std::cout << px;
0745 std::cout << ",";
0746 std::cout.width(5);
0747 std::cout << py;
0748 std::cout << ",";
0749 std::cout.width(5);
0750 std::cout << pz;
0751 std::cout << ")" << std::endl;
0752 std::cout << " quality = " << track->get_quality() << std::endl;
0753 std::cout << " nfromtruth = " << trackeval->get_nclusters_contribution(track, particle) << std::endl;
0754
0755 std::cout << " ---Associated-SvtxClusters-to-PHG4Hits-------------------------" << std::endl;
0756
0757 for (SvtxTrack::ConstClusterKeyIter local_iter = track->begin_cluster_keys();
0758 local_iter != track->end_cluster_keys();
0759 ++local_iter)
0760 {
0761 TrkrDefs::cluskey cluster_key = *local_iter;
0762 TrkrCluster* cluster = clustermap->findCluster(cluster_key);
0763
0764 Acts::Vector3 glob;
0765 glob = tgeometry->getGlobalPosition(cluster_key, cluster);
0766
0767 float x = glob(0);
0768 float y = glob(1);
0769 float z = glob(2);
0770
0771 std::cout << " #" << cluster_key << " xreco = (";
0772 std::cout.width(5);
0773 std::cout << x;
0774 std::cout << ",";
0775 std::cout.width(5);
0776 std::cout << y;
0777 std::cout << ",";
0778 std::cout.width(5);
0779 std::cout << z;
0780 std::cout << ") =>";
0781
0782 PHG4Hit* g4hit = clustereval->max_truth_hit_by_energy(cluster_key);
0783 if ((g4hit) && (g4hit->get_trkid() == particle->get_track_id()))
0784 {
0785 x = 0.5 * (g4hit->get_x(0) + g4hit->get_x(1));
0786 y = 0.5 * (g4hit->get_y(0) + g4hit->get_y(1));
0787 z = 0.5 * (g4hit->get_z(0) + g4hit->get_z(1));
0788
0789 std::cout << " #" << g4hit->get_hit_id()
0790 << " xtrue = (";
0791 std::cout.width(5);
0792 std::cout << x;
0793 std::cout << ",";
0794 std::cout.width(5);
0795 std::cout << y;
0796 std::cout << ",";
0797 std::cout.width(5);
0798 std::cout << z;
0799 std::cout << ") => Gtrack id = " << g4hit->get_trkid();
0800 }
0801 else
0802 {
0803 std::cout << " noise hit";
0804 }
0805 }
0806
0807 std::cout << std::endl;
0808 }
0809 }
0810 }
0811 }
0812
0813 std::cout << std::endl;
0814
0815 }
0816
0817 return;
0818 }
0819
0820 void SvtxEvaluator::fillOutputNtuples(PHCompositeNode* topNode)
0821 {
0822 if (Verbosity() > 1)
0823 {
0824 std::cout << "SvtxEvaluator::fillOutputNtuples() entered" << std::endl;
0825 }
0826
0827 ActsGeometry* tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0828 if (!tgeometry)
0829 {
0830 std::cout << "No Acts geometry on node tree. Can't continue."
0831 << std::endl;
0832 return;
0833 }
0834
0835 SvtxVertexEval* vertexeval = _svtxevalstack->get_vertex_eval();
0836 SvtxClusterEval* clustereval = _svtxevalstack->get_cluster_eval();
0837 SvtxTrackEval* trackeval = _svtxevalstack->get_track_eval();
0838 SvtxHitEval* hiteval = _svtxevalstack->get_hit_eval();
0839 SvtxTruthEval* trutheval = _svtxevalstack->get_truth_eval();
0840
0841 float nhit_tpc_all = 0;
0842 float nhit_tpc_in = 0;
0843 float nhit_tpc_mid = 0;
0844 float nhit_tpc_out = 0;
0845 float nclus_all = 0;
0846 float nclus_tpc = 0;
0847 float nclus_intt = 0;
0848 float nclus_maps = 0;
0849 float nclus_mms = 0;
0850 float nhit[100];
0851 for (float& i : nhit)
0852 {
0853 i = 0;
0854 }
0855 float occ11 = 0;
0856 float occ116 = 0;
0857 float occ21 = 0;
0858 float occ216 = 0;
0859 float occ31 = 0;
0860 float occ316 = 0;
0861
0862 TrkrHitSetContainer* hitmap_in = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0863 if (hitmap_in)
0864 {
0865 TrkrHitSetContainer::ConstRange all_hitsets = hitmap_in->getHitSets();
0866 for (TrkrHitSetContainer::ConstIterator hitsetiter = all_hitsets.first;
0867 hitsetiter != all_hitsets.second;
0868 ++hitsetiter)
0869 {
0870
0871 unsigned int layer = TrkrDefs::getLayer(hitsetiter->first);
0872 if (layer >= _nlayers_maps + _nlayers_intt && layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
0873 {
0874
0875 TrkrHitSet::ConstRange hitrangei = hitsetiter->second->getHits();
0876 for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0877 hitr != hitrangei.second;
0878 ++hitr)
0879 {
0880 nhit[layer]++;
0881 nhit_tpc_all++;
0882 if ((float) layer == _nlayers_maps + _nlayers_intt)
0883 {
0884 nhit_tpc_in++;
0885 }
0886 if ((float) layer == _nlayers_maps + _nlayers_intt + _nlayers_tpc - 1)
0887 {
0888 nhit_tpc_out++;
0889 }
0890 if ((float) layer == _nlayers_maps + _nlayers_intt + _nlayers_tpc / 2 - 1)
0891 {
0892 nhit_tpc_mid++;
0893 }
0894 }
0895 }
0896 }
0897 }
0898
0899 PHG4TpcCylinderGeomContainer* geom_container =
0900 findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0901 if (!geom_container)
0902 {
0903 std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0904 return;
0905 }
0906 PHG4TpcCylinderGeom* GeoLayer;
0907 int layer = _nlayers_maps + _nlayers_intt;
0908 GeoLayer = geom_container->GetLayerCellGeom(layer);
0909 int nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0910 float nhits = nhit[layer];
0911 occ11 = nhits / nbins;
0912 if (Verbosity() > 1)
0913 {
0914 std::cout << " occ11 = " << occ11
0915 << " nbins = " << nbins
0916 << " nhits = " << nhits
0917 << " layer = " << layer
0918 << std::endl;
0919 }
0920 layer = _nlayers_maps + _nlayers_intt + 15;
0921 GeoLayer = geom_container->GetLayerCellGeom(layer);
0922 nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0923 nhits = nhit[layer];
0924 occ116 = nhits / nbins;
0925
0926 layer = _nlayers_maps + _nlayers_intt + 16;
0927 GeoLayer = geom_container->GetLayerCellGeom(layer);
0928 nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0929 nhits = nhit[layer];
0930 occ21 = nhits / nbins;
0931 layer = _nlayers_maps + _nlayers_intt + 31;
0932 GeoLayer = geom_container->GetLayerCellGeom(layer);
0933 nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0934 nhits = nhit[layer];
0935 occ216 = nhits / nbins;
0936 layer = _nlayers_maps + _nlayers_intt + 32;
0937 GeoLayer = geom_container->GetLayerCellGeom(layer);
0938 nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0939 nhits = nhit[layer];
0940 occ31 = nhits / nbins;
0941 layer = _nlayers_maps + _nlayers_intt + 47;
0942 GeoLayer = geom_container->GetLayerCellGeom(layer);
0943 nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0944 nhits = nhit[layer];
0945 occ316 = nhits / nbins;
0946
0947 if (Verbosity() > 1)
0948 {
0949 std::cout << " occ11 = " << occ11
0950 << " occ116 = " << occ116
0951 << " occ21 = " << occ21
0952 << " occ216 = " << occ216
0953 << " occ31 = " << occ31
0954 << " occ316 = " << occ316
0955 << std::endl;
0956 }
0957 TrkrClusterContainer* clustermap_in = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
0958 if (!clustermap_in)
0959 {
0960 clustermap_in = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0961 }
0962
0963 if (clustermap_in)
0964 {
0965 nclus_all = clustermap_in->size();
0966
0967 for (const auto& hitsetkey : clustermap_in->getHitSetKeys())
0968 {
0969 auto range = clustermap_in->getClusters(hitsetkey);
0970 for (auto iter_cin = range.first; iter_cin != range.second; ++iter_cin)
0971 {
0972 TrkrDefs::cluskey cluster_key = iter_cin->first;
0973 unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
0974 if (_nlayers_maps > 0)
0975 {
0976 if (local_layer < _nlayers_maps)
0977 {
0978 nclus_maps++;
0979 }
0980 }
0981 if (_nlayers_intt > 0)
0982 {
0983 if (local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
0984 {
0985 nclus_intt++;
0986 }
0987 }
0988 if (_nlayers_tpc > 0)
0989 {
0990 if (local_layer >= _nlayers_maps + _nlayers_intt && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
0991 {
0992 nclus_tpc++;
0993 }
0994 }
0995 if (_nlayers_mms > 0)
0996 {
0997 if (local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
0998 {
0999 nclus_mms++;
1000 }
1001 }
1002 }
1003 }
1004 }
1005
1006
1007
1008 if (_ntp_info)
1009 {
1010 if (Verbosity() > 1)
1011 {
1012 std::cout << "Filling ntp_info " << std::endl;
1013 }
1014 float ntrk = 0;
1015 SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname.c_str());
1016 if (trackmap)
1017 {
1018 ntrk = (float) trackmap->size();
1019 }
1020 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1021 int nprim = truthinfo->GetNumPrimaryVertexParticles();
1022 if (Verbosity() > 0)
1023 {
1024 std::cout << "EVENTINFO SEED: " << m_fSeed << std::endl;
1025 std::cout << "EVENTINFO NHIT: " << std::setprecision(9) << nhit_tpc_all << std::endl;
1026 std::cout << "EVENTINFO NTRKGEN: " << nprim << std::endl;
1027 std::cout << "EVENTINFO NTRKREC: " << ntrk << std::endl;
1028 }
1029 float info_data[] = {(float) _ievent, m_fSeed,
1030 occ11, occ116, occ21, occ216, occ31, occ316,
1031 (float) nprim,
1032 0,
1033 ntrk,
1034 nhit_tpc_all,
1035 nhit_tpc_in,
1036 nhit_tpc_mid,
1037 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1038
1039 _ntp_info->Fill(info_data);
1040 }
1041
1042
1043
1044
1045 bool doit = true;
1046 if (_ntp_vertex && doit)
1047 {
1048 if (Verbosity() > 1)
1049 {
1050 std::cout << "Filling ntp_vertex " << std::endl;
1051 std::cout << "start vertex time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
1052 _timer->restart();
1053 }
1054
1055 GlobalVertexMap* gvertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
1056 SvtxVertexMap* vertexmap = nullptr;
1057 if (_use_initial_vertex)
1058 {
1059 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMap");
1060 }
1061 else if (_use_genfit_vertex)
1062 {
1063 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapRefit");
1064 }
1065 else
1066 {
1067 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs");
1068 }
1069
1070 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1071
1072 if (gvertexmap && vertexmap && truthinfo)
1073 {
1074 const auto prange = truthinfo->GetPrimaryParticleRange();
1075 std::map<int, unsigned int> embedvtxid_particle_count;
1076 std::map<int, unsigned int> embedvtxid_maps_particle_count;
1077 std::map<int, unsigned int> vertex_particle_count;
1078
1079 if (_do_vtx_eval_light == false)
1080 {
1081 for (auto iter = prange.first; iter != prange.second; ++iter)
1082 {
1083 const int point_id = iter->second->get_vtx_id();
1084 int gembed = truthinfo->isEmbededVtx(iter->second->get_vtx_id());
1085 ++vertex_particle_count[point_id];
1086 ++embedvtxid_particle_count[gembed];
1087 PHG4Particle* g4particle = iter->second;
1088
1089 if (_scan_for_embedded && gembed <= 0)
1090 {
1091 continue;
1092 }
1093
1094 std::set<TrkrDefs::cluskey> g4clusters = clustereval->all_clusters_from(g4particle);
1095 unsigned int nglmaps = 0;
1096
1097 std::vector<int> lmaps(_nlayers_maps + 1,0);
1098 for (const TrkrDefs::cluskey g4cluster : g4clusters)
1099 {
1100 unsigned int local_layer = TrkrDefs::getLayer(g4cluster);
1101
1102 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
1103 {
1104 lmaps[local_layer] = 1;
1105 }
1106 }
1107 if (_nlayers_maps > 0)
1108 {
1109 for (unsigned int i = 0; i < _nlayers_maps; i++)
1110 {
1111 nglmaps += lmaps[i];
1112 }
1113 }
1114 float gpx = g4particle->get_px();
1115 float gpy = g4particle->get_py();
1116 float gpz = g4particle->get_pz();
1117 float gpt = NAN;
1118 float geta = NAN;
1119
1120 if (gpx != 0 && gpy != 0)
1121 {
1122 TVector3 gv(gpx, gpy, gpz);
1123 gpt = gv.Pt();
1124 geta = gv.Eta();
1125
1126 }
1127
1128 if (nglmaps == 3 && fabs(geta) < 1.0 && gpt > 0.5)
1129 {
1130 ++embedvtxid_maps_particle_count[gembed];
1131 }
1132 }
1133 }
1134
1135 auto vrange = truthinfo->GetPrimaryVtxRange();
1136 std::map<int, bool> embedvtxid_found;
1137 std::map<int, int> embedvtxid_vertex_id;
1138 std::map<int, PHG4VtxPoint*> embedvtxid_vertex;
1139 for (auto iter = vrange.first; iter != vrange.second; ++iter)
1140 {
1141 const int point_id = iter->first;
1142 int gembed = truthinfo->isEmbededVtx(point_id);
1143
1144 if (_scan_for_embedded && gembed <= 0)
1145 {
1146 continue;
1147 }
1148
1149 auto search = embedvtxid_found.find(gembed);
1150 if (search != embedvtxid_found.end())
1151 {
1152 embedvtxid_vertex_id[gembed] = point_id;
1153 embedvtxid_vertex[gembed] = iter->second;
1154 }
1155 else
1156 {
1157 if (vertex_particle_count[embedvtxid_vertex_id[gembed]] < vertex_particle_count[point_id])
1158 {
1159 embedvtxid_vertex_id[gembed] = point_id;
1160 embedvtxid_vertex[gembed] = iter->second;
1161 }
1162 }
1163 embedvtxid_found[gembed] = false;
1164 }
1165
1166 unsigned int ngembed = 0;
1167 for (auto& iter : embedvtxid_found)
1168 {
1169 if (iter.first >= 0)
1170 {
1171 continue;
1172 }
1173 ++ngembed;
1174 }
1175
1176 for (auto& iter : *gvertexmap)
1177 {
1178 GlobalVertex* gvertex = iter.second;
1179 auto svtxviter = gvertex->find_vertexes(GlobalVertex::SVTX);
1180 if(svtxviter == gvertex->end_vertexes())
1181 {
1182 continue;
1183 }
1184
1185 GlobalVertex::VertexVector vertices = svtxviter->second;
1186 for(auto& vertex : vertices)
1187 {
1188 PHG4VtxPoint* point = vertexeval->max_truth_point_by_ntracks(vertex);
1189 int vertexID = vertex->get_id();
1190 float vx = vertex->get_x();
1191 float vy = vertex->get_y();
1192 float vz = vertex->get_z();
1193 float ntracks = vertex->size_tracks();
1194 float chi2 = vertex->get_chisq();
1195 float ndof = vertex->get_ndof();
1196 float gvx = NAN;
1197 float gvy = NAN;
1198 float gvz = NAN;
1199 float gvt = NAN;
1200 float gembed = NAN;
1201 float gntracks = truthinfo->GetNumPrimaryVertexParticles();
1202 float gntracksmaps = NAN;
1203 float gnembed = NAN;
1204 float nfromtruth = NAN;
1205 if (point)
1206 {
1207 const int point_id = point->get_id();
1208 gvx = point->get_x();
1209 gvy = point->get_y();
1210 gvz = point->get_z();
1211 gvt = point->get_t();
1212 gembed = truthinfo->isEmbededVtx(point_id);
1213 gntracks = embedvtxid_particle_count[(int) gembed];
1214 if (embedvtxid_maps_particle_count[(int) gembed] > 0 && fabs(gvt) < 2000. && fabs(gvz) < 13.0)
1215 {
1216 gntracksmaps = embedvtxid_maps_particle_count[(int) gembed];
1217 }
1218 gnembed = (float) ngembed;
1219 nfromtruth = vertexeval->get_ntracks_contribution(vertex, point);
1220 embedvtxid_found[(int) gembed] = true;
1221 }
1222
1223 float vertex_data[] = {(float) _ievent, m_fSeed,
1224 (float) vertexID,
1225 vx,
1226 vy,
1227 vz,
1228 ntracks,
1229 chi2,
1230 ndof,
1231 gvx,
1232 gvy,
1233 gvz,
1234 gvt,
1235 gembed,
1236 gntracks,
1237 gntracksmaps,
1238 gnembed,
1239 nfromtruth,
1240 nhit_tpc_all,
1241 nhit_tpc_in,
1242 nhit_tpc_mid,
1243 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1244
1245 _ntp_vertex->Fill(vertex_data);
1246 }
1247 }
1248
1249 if (!_scan_for_embedded)
1250 {
1251 for (std::map<int, bool>::iterator iter = embedvtxid_found.begin();
1252 iter != embedvtxid_found.end();
1253 ++iter)
1254 {
1255 if (embedvtxid_found[iter->first])
1256 {
1257 continue;
1258 }
1259
1260 float vx = NAN;
1261 float vy = NAN;
1262 float vz = NAN;
1263 float ntracks = NAN;
1264
1265 float gvx = NAN;
1266 float gvy = NAN;
1267 float gvz = NAN;
1268 float gvt = NAN;
1269 float gembed = iter->first;
1270 float gntracks = NAN;
1271 float gntracksmaps = NAN;
1272 float gnembed = NAN;
1273 float nfromtruth = NAN;
1274
1275 PHG4VtxPoint* point = embedvtxid_vertex[gembed];
1276
1277 if (point)
1278 {
1279 const int point_id = point->get_id();
1280 gvx = point->get_x();
1281 gvy = point->get_y();
1282 gvz = point->get_z();
1283 gvt = point->get_t();
1284 gembed = truthinfo->isEmbededVtx(point_id);
1285 gntracks = embedvtxid_particle_count[(int) gembed];
1286 if (embedvtxid_maps_particle_count[(int) gembed] > 0 && fabs(gvt) < 2000 && fabs(gvz) < 13.0)
1287 {
1288 gntracksmaps = embedvtxid_maps_particle_count[(int) gembed];
1289 }
1290 gnembed = (float) ngembed;
1291
1292 }
1293
1294 if (Verbosity() > 1)
1295 {
1296 std::cout << " adding vertex data " << std::endl;
1297 }
1298
1299 float vertex_data[] = {(float) _ievent, m_fSeed,
1300 vx,
1301 vy,
1302 vz,
1303 ntracks,
1304 gvx,
1305 gvy,
1306 gvz,
1307 gvt,
1308 gembed,
1309 gntracks,
1310 gntracksmaps,
1311 gnembed,
1312 nfromtruth,
1313 nhit_tpc_all,
1314 nhit_tpc_in,
1315 nhit_tpc_mid,
1316 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1317
1318 _ntp_vertex->Fill(vertex_data);
1319 }
1320 }
1321 }
1322 if (Verbosity() > 1)
1323 {
1324 _timer->stop();
1325 std::cout << "vertex time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
1326 }
1327 }
1328
1329
1330
1331
1332
1333 if (_ntp_gpoint)
1334 {
1335 if (Verbosity() > 1)
1336 {
1337 std::cout << "Filling ntp_gpoint " << std::endl;
1338 _timer->restart();
1339 }
1340 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
1341
1342 if (truthinfo)
1343 {
1344 auto vrange = truthinfo->GetPrimaryVtxRange();
1345 const auto prange = truthinfo->GetPrimaryParticleRange();
1346
1347 std::map<int, unsigned int> vertex_particle_count;
1348 for (auto iter = prange.first; iter != prange.second; ++iter)
1349 {
1350 ++vertex_particle_count[iter->second->get_vtx_id()];
1351 }
1352
1353 for (auto iter = vrange.first; iter != vrange.second; ++iter)
1354 {
1355 const int point_id = iter->first;
1356 PHG4VtxPoint* point = iter->second;
1357
1358
1359
1360 if (point)
1361 {
1362 const Vertex* vertex = vertexeval->best_vertex_from(point);
1363
1364 float gvx = point->get_x();
1365 float gvy = point->get_y();
1366 float gvz = point->get_z();
1367 float gvt = point->get_t();
1368 float gntracks = vertex_particle_count[point_id];
1369
1370 float gembed = truthinfo->isEmbededVtx(point_id);
1371 float vx = NAN;
1372 float vy = NAN;
1373 float vz = NAN;
1374 float ntracks = NAN;
1375 float nfromtruth = NAN;
1376
1377 if (vertex)
1378 {
1379 vx = vertex->get_x();
1380 vy = vertex->get_y();
1381 vz = vertex->get_z();
1382 ntracks = vertex->size_tracks();
1383 nfromtruth = vertexeval->get_ntracks_contribution(vertex, point);
1384 }
1385
1386 float gpoint_data[] = {(float) _ievent, m_fSeed,
1387 gvx,
1388 gvy,
1389 gvz,
1390 gvt,
1391 gntracks,
1392 gembed,
1393 vx,
1394 vy,
1395 vz,
1396 ntracks,
1397 nfromtruth,
1398 nhit_tpc_all,
1399 nhit_tpc_in,
1400 nhit_tpc_mid,
1401 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1402
1403 _ntp_gpoint->Fill(gpoint_data);
1404 }
1405 }
1406 }
1407 if (Verbosity() > 1)
1408 {
1409 _timer->stop();
1410 std::cout << "gpoint time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
1411 }
1412 }
1413
1414
1415
1416
1417
1418 if (_ntp_g4hit)
1419 {
1420 if (Verbosity() > 1)
1421 {
1422 std::cout << "Filling ntp_g4hit " << std::endl;
1423 _timer->restart();
1424 }
1425 std::set<PHG4Hit*> g4hits = trutheval->all_truth_hits();
1426 for (auto g4hit : g4hits)
1427 {
1428 PHG4Particle* g4particle = trutheval->get_particle(g4hit);
1429
1430 float g4hitID = g4hit->get_hit_id();
1431 float gx = g4hit->get_avg_x();
1432 float gy = g4hit->get_avg_y();
1433 float gz = g4hit->get_avg_z();
1434 TVector3 vg4(gx, gy, gz);
1435 float gt = g4hit->get_avg_t();
1436 float gpl = g4hit->get_path_length();
1437 TVector3 vin(g4hit->get_x(0), g4hit->get_y(0), g4hit->get_z(0));
1438 TVector3 vout(g4hit->get_x(1), g4hit->get_y(1), g4hit->get_z(1));
1439 float gdphi = vin.DeltaPhi(vout);
1440 float gdz = fabs(g4hit->get_z(1) - g4hit->get_z(0));
1441 float gedep = g4hit->get_edep();
1442 float glayer = g4hit->get_layer();
1443
1444 float gtrackID = g4hit->get_trkid();
1445
1446 float gflavor = NAN;
1447 float gpx = NAN;
1448 float gpy = NAN;
1449 float gpz = NAN;
1450 TVector3 vec(g4hit->get_avg_x(), g4hit->get_avg_y(), g4hit->get_avg_z());
1451 float geta = vec.Eta();
1452 float gphi = vec.Phi();
1453 float gvx = NAN;
1454 float gvy = NAN;
1455 float gvz = NAN;
1456
1457 float gembed = NAN;
1458 float gprimary = NAN;
1459
1460 float gfpx = 0.;
1461 float gfpy = 0.;
1462 float gfpz = 0.;
1463 float gfx = 0.;
1464 float gfy = 0.;
1465 float gfz = 0.;
1466
1467 if (g4particle)
1468 {
1469 if (_scan_for_embedded)
1470 {
1471 if (trutheval->get_embed(g4particle) <= 0)
1472 {
1473 continue;
1474 }
1475 }
1476
1477 gflavor = g4particle->get_pid();
1478 gpx = g4particle->get_px();
1479 gpy = g4particle->get_py();
1480 gpz = g4particle->get_pz();
1481
1482 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1483
1484 if (vtx)
1485 {
1486 gvx = vtx->get_x();
1487 gvy = vtx->get_y();
1488 gvz = vtx->get_z();
1489 }
1490 PHG4Hit* outerhit = nullptr;
1491 if (_do_eval_light == false)
1492 {
1493 outerhit = trutheval->get_outermost_truth_hit(g4particle);
1494 }
1495
1496 if (outerhit)
1497 {
1498 gfpx = outerhit->get_px(1);
1499 gfpy = outerhit->get_py(1);
1500 gfpz = outerhit->get_pz(1);
1501 gfx = outerhit->get_x(1);
1502 gfy = outerhit->get_y(1);
1503 gfz = outerhit->get_z(1);
1504 }
1505
1506 gembed = trutheval->get_embed(g4particle);
1507 gprimary = trutheval->is_primary(g4particle);
1508 }
1509
1510 std::set<TrkrDefs::cluskey> clusters = clustereval->all_clusters_from(g4hit);
1511 float nclusters = clusters.size();
1512
1513
1514 TrkrDefs::cluskey cluster_key = clustereval->best_cluster_from(g4hit);
1515
1516 float clusID = NAN;
1517 float x = NAN;
1518 float y = NAN;
1519 float z = NAN;
1520 float eta = NAN;
1521 float phi = NAN;
1522 float e = NAN;
1523 float adc = NAN;
1524 float local_layer = NAN;
1525 float size = NAN;
1526 float efromtruth = NAN;
1527 float dphitru = NAN;
1528 float detatru = NAN;
1529 float dztru = NAN;
1530 float drtru = NAN;
1531
1532 TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
1533 if (!clustermap)
1534 {
1535 clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1536 }
1537
1538 TrkrCluster* cluster = clustermap->findCluster(cluster_key);
1539
1540 if (cluster)
1541 {
1542 clusID = cluster_key;
1543
1544 Acts::Vector3 global;
1545 global = tgeometry->getGlobalPosition(cluster_key, cluster);
1546 x = global(0);
1547 y = global(1);
1548 z = global(2);
1549 TVector3 vec2(x, y, z);
1550 eta = vec2.Eta();
1551 phi = vec2.Phi();
1552 e = cluster->getAdc();
1553 adc = cluster->getAdc();
1554 local_layer = (float) TrkrDefs::getLayer(cluster_key);
1555 size = 0.0;
1556
1557
1558 TrkrClusterHitAssoc* cluster_hit_map = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
1559 std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
1560 hitrange = cluster_hit_map->getHits(cluster_key);
1561 for (std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
1562 clushititer = hitrange.first;
1563 clushititer != hitrange.second; ++clushititer)
1564 {
1565 ++size;
1566 }
1567
1568 dphitru = vec2.DeltaPhi(vg4);
1569 detatru = eta - geta;
1570 dztru = z - gz;
1571 drtru = vec2.DeltaR(vg4);
1572 if (g4particle)
1573 {
1574 efromtruth = clustereval->get_energy_contribution(cluster_key, g4particle);
1575 }
1576 }
1577
1578 float g4hit_data[] = {(float) _ievent, m_fSeed,
1579 g4hitID,
1580 gx,
1581 gy,
1582 gz,
1583 gt,
1584 gpl,
1585 gedep,
1586 geta,
1587 gphi,
1588 gdphi,
1589 gdz,
1590 glayer,
1591 gtrackID,
1592 gflavor,
1593 gpx,
1594 gpy,
1595 gpz,
1596 gvx,
1597 gvy,
1598 gvz,
1599 gfpx,
1600 gfpy,
1601 gfpz,
1602 gfx,
1603 gfy,
1604 gfz,
1605 gembed,
1606 gprimary,
1607 nclusters,
1608 clusID,
1609 x,
1610 y,
1611 z,
1612 eta,
1613 phi,
1614 e,
1615 adc,
1616 local_layer,
1617 size,
1618 efromtruth,
1619 dphitru,
1620 detatru,
1621 dztru,
1622 drtru,
1623 nhit_tpc_all,
1624 nhit_tpc_in,
1625 nhit_tpc_mid,
1626 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1627
1628 _ntp_g4hit->Fill(g4hit_data);
1629 }
1630 if (Verbosity() > 1)
1631 {
1632 _timer->stop();
1633 std::cout << "g4hit time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
1634 }
1635 }
1636
1637
1638
1639
1640
1641 if (_ntp_hit)
1642 {
1643 if (Verbosity() >= 1)
1644 {
1645 std::cout << "Filling ntp_hit " << std::endl;
1646 _timer->restart();
1647 }
1648
1649 TrkrHitSetContainer* hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1650 PHG4TpcCylinderGeomContainer* local_geom_container =
1651 findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
1652 if (!local_geom_container)
1653 {
1654 std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
1655 return;
1656 }
1657
1658 if (hitmap)
1659 {
1660 TrkrHitSetContainer::ConstRange all_hitsets = hitmap->getHitSets();
1661 for (TrkrHitSetContainer::ConstIterator iter = all_hitsets.first;
1662 iter != all_hitsets.second;
1663 ++iter)
1664 {
1665 TrkrDefs::hitsetkey hitset_key = iter->first;
1666 TrkrHitSet* hitset = iter->second;
1667
1668
1669 TrkrHitSet::ConstRange hitrangei = hitset->getHits();
1670 for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
1671 hitr != hitrangei.second;
1672 ++hitr)
1673 {
1674 TrkrDefs::hitkey hit_key = hitr->first;
1675 TrkrHit* hit = hitr->second;
1676 PHG4Hit* g4hit = hiteval->max_truth_hit_by_energy(hit_key);
1677 PHG4Particle* g4particle = trutheval->get_particle(g4hit);
1678 float event = _ievent;
1679 float hitID = hit_key;
1680 float e = hit->getEnergy();
1681 float adc = hit->getAdc();
1682 float local_layer = TrkrDefs::getLayer(hitset_key);
1683 float sector = TpcDefs::getSectorId(hitset_key);
1684 float side = TpcDefs::getSide(hitset_key);
1685 float cellID = 0;
1686 float ecell = hit->getAdc();
1687
1688 float phibin = NAN;
1689 float zbin = NAN;
1690 float phi = NAN;
1691 float r = NAN;
1692 float x = NAN;
1693 float y = NAN;
1694 float z = NAN;
1695
1696 if (local_layer >= _nlayers_maps + _nlayers_intt && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
1697 {
1698 PHG4TpcCylinderGeom* local_GeoLayer = local_geom_container->GetLayerCellGeom(local_layer);
1699 phibin = (float) TpcDefs::getPad(hit_key);
1700 zbin = (float) TpcDefs::getTBin(hit_key);
1701 phi = local_GeoLayer->get_phicenter(phibin, side);
1702 z = local_GeoLayer->get_zcenter(zbin);
1703 r = local_GeoLayer->get_radius();
1704 x = r*cos(phi);
1705 y = r*sin(phi);
1706 }
1707
1708 float g4hitID = NAN;
1709 float gedep = NAN;
1710 float gx = NAN;
1711 float gy = NAN;
1712 float gz = NAN;
1713 float gt = NAN;
1714 float gtrackID = NAN;
1715 float gflavor = NAN;
1716 float gpx = NAN;
1717 float gpy = NAN;
1718 float gpz = NAN;
1719 float gvx = NAN;
1720 float gvy = NAN;
1721 float gvz = NAN;
1722 float gvt = NAN;
1723 float gfpx = NAN;
1724 float gfpy = NAN;
1725 float gfpz = NAN;
1726 float gfx = NAN;
1727 float gfy = NAN;
1728 float gfz = NAN;
1729 float gembed = NAN;
1730 float gprimary = NAN;
1731
1732 float efromtruth = NAN;
1733
1734 if (g4hit)
1735 {
1736 g4hitID = g4hit->get_hit_id();
1737 gedep = g4hit->get_edep();
1738 gx = g4hit->get_avg_x();
1739 gy = g4hit->get_avg_y();
1740 gz = g4hit->get_avg_z();
1741 gt = g4hit->get_avg_t();
1742
1743 if (g4particle)
1744 {
1745 if (_scan_for_embedded)
1746 {
1747 if (trutheval->get_embed(g4particle) <= 0)
1748 {
1749 continue;
1750 }
1751 }
1752
1753 gtrackID = g4particle->get_track_id();
1754 gflavor = g4particle->get_pid();
1755 gpx = g4particle->get_px();
1756 gpy = g4particle->get_py();
1757 gpz = g4particle->get_pz();
1758
1759 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
1760
1761 if (vtx)
1762 {
1763 gvx = vtx->get_x();
1764 gvy = vtx->get_y();
1765 gvz = vtx->get_z();
1766 gvt = vtx->get_t();
1767 }
1768
1769 PHG4Hit* outerhit = nullptr;
1770 if (_do_eval_light == false)
1771 {
1772 outerhit = trutheval->get_outermost_truth_hit(g4particle);
1773 }
1774 if (outerhit)
1775 {
1776 gfpx = outerhit->get_px(1);
1777 gfpy = outerhit->get_py(1);
1778 gfpz = outerhit->get_pz(1);
1779 gfx = outerhit->get_x(1);
1780 gfy = outerhit->get_y(1);
1781 gfz = outerhit->get_z(1);
1782 }
1783 gembed = trutheval->get_embed(g4particle);
1784 gprimary = trutheval->is_primary(g4particle);
1785 }
1786 }
1787
1788 if (g4particle)
1789 {
1790 efromtruth = hiteval->get_energy_contribution(hit_key, g4particle);
1791 }
1792
1793 float hit_data[] = {
1794 event,
1795 (float) _iseed,
1796 hitID,
1797 e,
1798 adc,
1799 local_layer,
1800 sector,
1801 side,
1802 cellID,
1803 ecell,
1804 (float) phibin,
1805 (float) zbin,
1806 phi,
1807 x,
1808 y,
1809 z,
1810 g4hitID,
1811 gedep,
1812 gx,
1813 gy,
1814 gz,
1815 gt,
1816 gtrackID,
1817 gflavor,
1818 gpx,
1819 gpy,
1820 gpz,
1821 gvx,
1822 gvy,
1823 gvz,
1824 gvt,
1825 gfpx,
1826 gfpy,
1827 gfpz,
1828 gfx,
1829 gfy,
1830 gfz,
1831 gembed,
1832 gprimary,
1833 efromtruth,
1834 nhit_tpc_all,
1835 nhit_tpc_in,
1836 nhit_tpc_mid,
1837 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
1838
1839 _ntp_hit->Fill(hit_data);
1840 }
1841 }
1842 }
1843 if (Verbosity() >= 1)
1844 {
1845 _timer->stop();
1846 std::cout << "hit time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
1847 }
1848 }
1849
1850
1851
1852
1853
1854 if (Verbosity() >= 1)
1855 {
1856 std::cout << "check for ntp_cluster" << std::endl;
1857 _timer->restart();
1858 }
1859
1860 if (_ntp_cluster && !_scan_for_embedded)
1861 {
1862 if (Verbosity() > 1)
1863 {
1864 std::cout << "Filling ntp_cluster (all of them) " << std::endl;
1865 }
1866
1867 TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
1868 if (!clustermap)
1869 {
1870 clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
1871 }
1872
1873 TrkrClusterHitAssoc* clusterhitmap = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
1874 TrkrHitSetContainer* hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1875 TrkrClusterIterationMapv1* _iteration_map = findNode::getClass<TrkrClusterIterationMapv1>(topNode, "CLUSTER_ITERATION_MAP");
1876 ClusterErrorPara ClusErrPara;
1877
1878 if (Verbosity() > 1)
1879 {
1880 if (clustermap != nullptr)
1881 {
1882 std::cout << "got clustermap" << std::endl;
1883 }
1884 else
1885 {
1886 std::cout << "no clustermap" << std::endl;
1887 }
1888 if (clusterhitmap != nullptr)
1889 {
1890 std::cout << "got clusterhitmap" << std::endl;
1891 }
1892 else
1893 {
1894 std::cout << "no clusterhitmap" << std::endl;
1895 }
1896 if (hitsets != nullptr)
1897 {
1898 std::cout << "got hitsets" << std::endl;
1899 }
1900 else
1901 {
1902 std::cout << "no hitsets" << std::endl;
1903 }
1904 }
1905
1906 if (clustermap && clusterhitmap && hitsets)
1907 {
1908 for (const auto& hitsetkey : clustermap->getHitSetKeys())
1909 {
1910 int hitsetlayer = TrkrDefs::getLayer(hitsetkey);
1911 auto range = clustermap->getClusters(hitsetkey);
1912 for (auto iter = range.first; iter != range.second; ++iter)
1913 {
1914 TrkrDefs::cluskey cluster_key = iter->first;
1915 TrkrCluster* cluster = clustermap->findCluster(cluster_key);
1916 SvtxTrack* track = trackeval->best_track_from(cluster_key);
1917 PHG4Particle* g4particle = clustereval->max_truth_particle_by_cluster_energy(cluster_key);
1918 float niter = 0;
1919 if (_iteration_map != nullptr)
1920 {
1921 niter = _iteration_map->getIteration(cluster_key);
1922 }
1923 float hitID = (float) cluster_key;
1924
1925 Acts::Vector3 cglob;
1926 cglob = tgeometry->getGlobalPosition(cluster_key, cluster);
1927 float x = cglob(0);
1928 float y = cglob(1);
1929 float z = cglob(2);
1930 TVector3 pos(x, y, z);
1931 float r = pos.Perp();
1932 float phi = pos.Phi();
1933 float eta = pos.Eta();
1934 float theta = pos.Theta();
1935
1936 float ex = 0;
1937 float ey = 0;
1938 float ez = 0;
1939 float ephi = 0;
1940 float pez = 0;
1941 float pephi = 0;
1942 float size = 0;
1943 float phisize = 0;
1944 float zsize = 0;
1945 float maxadc = -999;
1946 float redge = NAN;
1947 float pedge = NAN;
1948 float ovlp = NAN;
1949
1950
1951 auto para_errors = ClusErrPara.get_clusterv5_modified_error(cluster, r, cluster_key);
1952 phisize = cluster->getPhiSize();
1953 zsize = cluster->getZSize();
1954
1955 ez = sqrt(para_errors.second);
1956 ephi = sqrt(para_errors.first);
1957 maxadc = cluster->getMaxAdc();
1958 pedge = cluster->getEdge();
1959 ovlp = cluster->getOverlap();
1960
1961 if(hitsetlayer==7||hitsetlayer==22||hitsetlayer==23||hitsetlayer==38||hitsetlayer==39) redge = 1;
1962
1963 float e = cluster->getAdc();
1964 float adc = cluster->getAdc();
1965 float local_layer = (float) TrkrDefs::getLayer(cluster_key);
1966 float sector = TpcDefs::getSectorId(cluster_key);
1967 float side = TpcDefs::getSide(cluster_key);
1968
1969
1970 TrkrDefs::hitsetkey local_hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
1971 int hitsetlayer2 = TrkrDefs::getLayer(local_hitsetkey);
1972 if (hitsetlayer != local_layer)
1973 {
1974 std::cout << "WARNING hitset layer " << hitsetlayer << "| " << hitsetlayer2 << " layer " << local_layer << std::endl;
1975 }
1976
1977
1978
1979
1980 float sumadc = 0;
1981 TrkrHitSetContainer::Iterator hitset = hitsets->findOrAddHitSet(local_hitsetkey);
1982 std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
1983 hitrange = clusterhitmap->getHits(cluster_key);
1984 for (std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
1985 clushititer = hitrange.first;
1986 clushititer != hitrange.second; ++clushititer)
1987 {
1988 TrkrHit* hit = hitset->second->getHit(clushititer->second);
1989 if (!hit)
1990 {
1991 continue;
1992 }
1993
1994 ++size;
1995 sumadc += (hit->getAdc() - 70);
1996 if ((hit->getAdc() - 70) > maxadc)
1997 {
1998 maxadc = (hit->getAdc() - 70);
1999 }
2000 }
2001 e = sumadc;
2002
2003 float trackID = NAN;
2004 if (track != nullptr)
2005 {
2006 trackID = track->get_id();
2007 }
2008
2009 float g4hitID = NAN;
2010 float gx = NAN;
2011 float gy = NAN;
2012 float gz = NAN;
2013 float gr = NAN;
2014 float gphi = NAN;
2015
2016 float geta = NAN;
2017 float gt = NAN;
2018 float gtrackID = NAN;
2019 float gflavor = NAN;
2020 float gpx = NAN;
2021 float gpy = NAN;
2022 float gpz = NAN;
2023 float gvx = NAN;
2024 float gvy = NAN;
2025 float gvz = NAN;
2026 float gvt = NAN;
2027 float gfpx = NAN;
2028 float gfpy = NAN;
2029 float gfpz = NAN;
2030 float gfx = NAN;
2031 float gfy = NAN;
2032 float gfz = NAN;
2033 float gembed = NAN;
2034 float gprimary = NAN;
2035
2036 float efromtruth = NAN;
2037
2038 if (Verbosity() > 1)
2039 {
2040 std::cout << PHWHERE << " **** reco: layer " << local_layer << std::endl;
2041 std::cout << " reco cluster key " << cluster_key << " r " << r << " x " << x << " y " << y << " z " << z << " phi " << phi << " adc " << adc << std::endl;
2042 }
2043 float nparticles = NAN;
2044
2045
2046 const auto [truth_ckey, truth_cluster] = clustereval->max_truth_cluster_by_energy(cluster_key);
2047 if (truth_cluster)
2048 {
2049 if (Verbosity() > 1)
2050 {
2051 std::cout << "Found matching truth cluster with key " << truth_ckey << " for reco cluster key " << cluster_key << " in layer " << local_layer << std::endl;
2052 }
2053
2054 g4hitID = 0;
2055 gx = truth_cluster->getX();
2056 gy = truth_cluster->getY();
2057 gz = truth_cluster->getZ();
2058 efromtruth = truth_cluster->getError(0, 0);
2059
2060 TVector3 gpos(gx, gy, gz);
2061 gr = gpos.Perp();
2062 gphi = gpos.Phi();
2063 geta = gpos.Eta();
2064
2065 if (g4particle)
2066 {
2067 gtrackID = g4particle->get_track_id();
2068 gflavor = g4particle->get_pid();
2069 gpx = g4particle->get_px();
2070 gpy = g4particle->get_py();
2071 gpz = g4particle->get_pz();
2072
2073 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
2074 if (vtx)
2075 {
2076 gvx = vtx->get_x();
2077 gvy = vtx->get_y();
2078 gvz = vtx->get_z();
2079 gvt = vtx->get_t();
2080 }
2081
2082 PHG4Hit* outerhit = nullptr;
2083 if (_do_eval_light == false)
2084 {
2085 outerhit = trutheval->get_outermost_truth_hit(g4particle);
2086 }
2087 if (outerhit)
2088 {
2089 gfpx = outerhit->get_px(1);
2090 gfpy = outerhit->get_py(1);
2091 gfpz = outerhit->get_pz(1);
2092 gfx = outerhit->get_x(1);
2093 gfy = outerhit->get_y(1);
2094 gfz = outerhit->get_z(1);
2095 }
2096
2097 gembed = trutheval->get_embed(g4particle);
2098 gprimary = trutheval->is_primary(g4particle);
2099 }
2100
2101 if (Verbosity() > 1)
2102 {
2103 std::cout << " truth cluster key " << truth_ckey << " gr " << gr << " gx " << gx << " gy " << gy << " gz " << gz << " gphi " << gphi << " efromtruth " << efromtruth << std::endl;
2104 }
2105 }
2106
2107 if (g4particle)
2108 {
2109 }
2110 nparticles = clustereval->all_truth_particles(cluster_key).size();
2111
2112 float cluster_data[] = {(float) _ievent,
2113 (float) _iseed,
2114 hitID,
2115 x,
2116 y,
2117 z,
2118 r,
2119 phi,
2120 eta,
2121 theta,
2122 ex,
2123 ey,
2124 ez,
2125 ephi,
2126 pez,
2127 pephi,
2128 e,
2129 adc,
2130 maxadc,
2131 local_layer,
2132 sector,
2133 side,
2134 size,
2135 phisize,
2136 zsize,
2137 pedge,
2138 redge,
2139 ovlp,
2140 trackID,
2141 niter,
2142 g4hitID,
2143 gx,
2144 gy,
2145 gz,
2146 gr,
2147 gphi,
2148 geta,
2149 gt,
2150 gtrackID,
2151 gflavor,
2152 gpx,
2153 gpy,
2154 gpz,
2155 gvx,
2156 gvy,
2157 gvz,
2158 gvt,
2159 gfpx,
2160 gfpy,
2161 gfpz,
2162 gfx,
2163 gfy,
2164 gfz,
2165 gembed,
2166 gprimary,
2167 efromtruth,
2168 nparticles,
2169 nhit_tpc_all,
2170 nhit_tpc_in,
2171 nhit_tpc_mid,
2172 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
2173
2174 _ntp_cluster->Fill(cluster_data);
2175 }
2176 }
2177 }
2178 }
2179 else if (_ntp_cluster && _scan_for_embedded)
2180 {
2181 if (Verbosity() > 1)
2182 {
2183 std::cout << "Filling ntp_cluster (embedded only) " << std::endl;
2184 }
2185
2186
2187
2188
2189
2190
2191 SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname.c_str());
2192
2193 TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
2194 if (!clustermap)
2195 {
2196 clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
2197 }
2198 ClusterErrorPara ClusErrPara;
2199
2200 TrkrClusterHitAssoc* clusterhitmap = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
2201 TrkrHitSetContainer* hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
2202 TrkrClusterIterationMapv1* _iteration_map = findNode::getClass<TrkrClusterIterationMapv1>(topNode, "CLUSTER_ITERATION_MAP");
2203
2204 if (trackmap != nullptr && clustermap != nullptr && clusterhitmap != nullptr && hitsets != nullptr)
2205 {
2206 for (auto& iter : *trackmap)
2207 {
2208 SvtxTrack* track = iter.second;
2209
2210 PHG4Particle* truth = trackeval->max_truth_particle_by_nclusters(track);
2211 if (truth)
2212 {
2213 if (trutheval->get_embed(truth) <= 0)
2214 {
2215 continue;
2216 }
2217 }
2218
2219 std::vector<TrkrDefs::cluskey> clusters;
2220 auto siseed = track->get_silicon_seed();
2221 if (siseed)
2222 {
2223 for (auto local_iter = siseed->begin_cluster_keys();
2224 local_iter != siseed->end_cluster_keys();
2225 ++local_iter)
2226 {
2227 TrkrDefs::cluskey cluster_key = *local_iter;
2228 clusters.push_back(cluster_key);
2229 }
2230 }
2231 auto tpcseed = track->get_tpc_seed();
2232 if (tpcseed)
2233 {
2234 for (auto local_iter = tpcseed->begin_cluster_keys();
2235 local_iter != tpcseed->end_cluster_keys();
2236 ++local_iter)
2237 {
2238 TrkrDefs::cluskey cluster_key = *local_iter;
2239 clusters.push_back(cluster_key);
2240 }
2241 }
2242
2243
2244 for (unsigned long cluster_key : clusters)
2245 {
2246 TrkrCluster* cluster = clustermap->findCluster(cluster_key);
2247 if (!cluster)
2248 {
2249 continue;
2250 }
2251
2252 PHG4Hit* g4hit = clustereval->max_truth_hit_by_energy(cluster_key);
2253 PHG4Particle* g4particle = trutheval->get_particle(g4hit);
2254
2255 float niter = 0;
2256 if (_iteration_map != nullptr)
2257 {
2258 niter = _iteration_map->getIteration(cluster_key);
2259 }
2260 float hitID = (float) TrkrDefs::getClusIndex(cluster_key);
2261 Acts::Vector3 glob = tgeometry->getGlobalPosition(cluster_key, cluster);
2262 float x = glob(0);
2263 float y = glob(1);
2264 float z = glob(2);
2265 TVector3 pos(x, y, z);
2266 float r = pos.Perp();
2267 float phi = pos.Phi();
2268 float eta = pos.Eta();
2269 float theta = pos.Theta();
2270 float ex = 0;
2271 float ey = 0;
2272 float ez = 0;
2273 float ephi = 0;
2274 float pez = 0;
2275 float pephi = 0;
2276 float size = 0;
2277 float phisize = 0;
2278 float zsize = 0;
2279 float maxadc = -999;
2280 float redge = NAN;
2281 float pedge = NAN;
2282 float ovlp = NAN;
2283
2284 auto para_errors = ClusErrPara.get_clusterv5_modified_error(cluster, r, cluster_key);
2285 phisize = cluster->getPhiSize();
2286 zsize = cluster->getZSize();
2287
2288 ez = sqrt(para_errors.second);
2289 ephi = sqrt(para_errors.first);
2290 maxadc = cluster->getMaxAdc();
2291 pedge = cluster->getEdge();
2292 ovlp = cluster->getOverlap();
2293
2294 float e = cluster->getAdc();
2295 float adc = cluster->getAdc();
2296 float local_layer = (float) TrkrDefs::getLayer(cluster_key);
2297 float sector = TpcDefs::getSectorId(cluster_key);
2298 float side = TpcDefs::getSide(cluster_key);
2299
2300 if(local_layer==7||local_layer==22||local_layer==23||local_layer==38||local_layer==39) redge = 1;
2301
2302
2303 TrkrDefs::hitsetkey hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
2304 TrkrHitSetContainer::Iterator hitset = hitsets->findOrAddHitSet(hitsetkey);
2305 std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
2306 hitrange = clusterhitmap->getHits(cluster_key);
2307 for (std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
2308 clushititer = hitrange.first;
2309 clushititer != hitrange.second; ++clushititer)
2310 {
2311 TrkrHit* hit = hitset->second->getHit(clushititer->second);
2312 ++size;
2313 if (hit->getAdc() > maxadc)
2314 {
2315 maxadc = hit->getAdc();
2316 }
2317 }
2318
2319 float trackID = NAN;
2320 trackID = track->get_id();
2321
2322 float g4hitID = NAN;
2323 float gx = NAN;
2324 float gy = NAN;
2325 float gz = NAN;
2326 float gr = NAN;
2327 float gphi = NAN;
2328
2329 float geta = NAN;
2330 float gt = NAN;
2331 float gtrackID = NAN;
2332 float gflavor = NAN;
2333 float gpx = NAN;
2334 float gpy = NAN;
2335 float gpz = NAN;
2336 float gvx = NAN;
2337 float gvy = NAN;
2338 float gvz = NAN;
2339 float gvt = NAN;
2340 float gfpx = NAN;
2341 float gfpy = NAN;
2342 float gfpz = NAN;
2343 float gfx = NAN;
2344 float gfy = NAN;
2345 float gfz = NAN;
2346 float gembed = NAN;
2347 float gprimary = NAN;
2348
2349 float efromtruth = NAN;
2350
2351
2352 const auto [truth_ckey, truth_cluster] = clustereval->max_truth_cluster_by_energy(cluster_key);
2353 if (truth_cluster)
2354 {
2355 if (Verbosity() > 1)
2356 {
2357 std::cout << " Found matching truth cluster with key " << truth_ckey << " for reco cluster key " << cluster_key << " in layer " << local_layer << std::endl;
2358 }
2359
2360 g4hitID = 0;
2361 gx = truth_cluster->getX();
2362 gy = truth_cluster->getY();
2363 gz = truth_cluster->getZ();
2364 efromtruth = truth_cluster->getError(0, 0);
2365
2366 TVector3 gpos(gx, gy, gz);
2367 gr = gpos.Perp();
2368 gphi = gpos.Phi();
2369 geta = gpos.Eta();
2370
2371 if (g4particle)
2372 {
2373 gtrackID = g4particle->get_track_id();
2374 gflavor = g4particle->get_pid();
2375 gpx = g4particle->get_px();
2376 gpy = g4particle->get_py();
2377 gpz = g4particle->get_pz();
2378
2379 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
2380 if (vtx)
2381 {
2382 gvx = vtx->get_x();
2383 gvy = vtx->get_y();
2384 gvz = vtx->get_z();
2385 gvt = vtx->get_t();
2386 }
2387 PHG4Hit* outerhit = nullptr;
2388 if (_do_eval_light == false)
2389 {
2390 outerhit = trutheval->get_outermost_truth_hit(g4particle);
2391 }
2392 if (outerhit)
2393 {
2394 gfpx = outerhit->get_px(1);
2395 gfpy = outerhit->get_py(1);
2396 gfpz = outerhit->get_pz(1);
2397 gfx = outerhit->get_x(1);
2398 gfy = outerhit->get_y(1);
2399 gfz = outerhit->get_z(1);
2400 }
2401
2402 gembed = trutheval->get_embed(g4particle);
2403 gprimary = trutheval->is_primary(g4particle);
2404 }
2405 }
2406
2407 float nparticles = clustereval->all_truth_particles(cluster_key).size();
2408
2409 float cluster_data[] = {(float) _ievent,
2410 (float) _iseed,
2411 hitID,
2412 x,
2413 y,
2414 z,
2415 r,
2416 phi,
2417 eta,
2418 theta,
2419 ex,
2420 ey,
2421 ez,
2422 ephi,
2423 pez,
2424 pephi,
2425 e,
2426 adc,
2427 maxadc,
2428 local_layer,
2429 sector,
2430 side,
2431 size,
2432 phisize,
2433 zsize,
2434 pedge,
2435 redge,
2436 ovlp,
2437 trackID,
2438 niter,
2439 g4hitID,
2440 gx,
2441 gy,
2442 gz,
2443 gr,
2444 gphi,
2445 geta,
2446 gt,
2447 gtrackID,
2448 gflavor,
2449 gpx,
2450 gpy,
2451 gpz,
2452 gvx,
2453 gvy,
2454 gvz,
2455 gvt,
2456 gfpx,
2457 gfpy,
2458 gfpz,
2459 gfx,
2460 gfy,
2461 gfz,
2462 gembed,
2463 gprimary,
2464 efromtruth,
2465 nparticles,
2466 nhit_tpc_all,
2467 nhit_tpc_in,
2468 nhit_tpc_mid,
2469 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
2470
2471 _ntp_cluster->Fill(cluster_data);
2472 }
2473 }
2474 }
2475 }
2476 if (Verbosity() >= 1)
2477 {
2478 _timer->stop();
2479 std::cout << "cluster time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
2480 }
2481
2482
2483
2484
2485 if (Verbosity() > 1)
2486 {
2487 std::cout << "check for ntp_g4cluster" << std::endl;
2488 _timer->restart();
2489 }
2490
2491 if (_ntp_g4cluster)
2492 {
2493 if (Verbosity() >= 1)
2494 {
2495 std::cout << "Filling ntp_g4cluster " << std::endl;
2496 }
2497 ClusterErrorPara ClusErrPara;
2498 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
2499 PHG4TruthInfoContainer::ConstRange range = truthinfo->GetParticleRange();
2500 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
2501 iter != range.second;
2502 ++iter)
2503 {
2504 PHG4Particle* g4particle = iter->second;
2505
2506 if (_scan_for_embedded)
2507 {
2508 if (trutheval->get_embed(g4particle) <= 0)
2509 {
2510 continue;
2511 }
2512 }
2513
2514 float gtrackID = g4particle->get_track_id();
2515 float gflavor = g4particle->get_pid();
2516 float gembed = trutheval->get_embed(g4particle);
2517 float gprimary = trutheval->is_primary(g4particle);
2518
2519 if (Verbosity() > 1)
2520 {
2521 std::cout << PHWHERE << " PHG4Particle ID " << gtrackID << " gflavor " << gflavor << " gprimary " << gprimary << std::endl;
2522 }
2523
2524
2525 std::map<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster> > truth_clusters = trutheval->all_truth_clusters(g4particle);
2526
2527
2528 for (const auto& [ckey, gclus] : truth_clusters)
2529 {
2530 unsigned int local_layer = TrkrDefs::getLayer(ckey);
2531
2532 float gx = gclus->getX();
2533 float gy = gclus->getY();
2534 float gz = gclus->getZ();
2535 float gt = NAN;
2536 float gedep = gclus->getError(0, 0);
2537 float gadc = (float) gclus->getAdc();
2538
2539 TVector3 gpos(gx, gy, gz);
2540 float gr = sqrt(gx * gx + gy * gy);
2541 float gphi = gpos.Phi();
2542 float geta = gpos.Eta();
2543
2544 if (Verbosity() > 1)
2545 {
2546 std::cout << PHWHERE << " **** truth: layer " << local_layer << std::endl;
2547 std::cout << " truth cluster key " << ckey << " gr " << gr << " gx " << gx << " gy " << gy << " gz " << gz
2548 << " gphi " << gphi << " gedep " << gedep << " gadc " << gadc << std::endl;
2549 }
2550
2551 float gphisize = gclus->getSize(1, 1);
2552 float gzsize = gclus->getSize(2, 2);
2553
2554
2555
2556 float x = NAN;
2557 float y = NAN;
2558 float z = NAN;
2559 float r = NAN;
2560 float phi = NAN;
2561 float eta = NAN;
2562 float ex = NAN;
2563 float ey = NAN;
2564 float ez = NAN;
2565 float ephi = NAN;
2566 float adc = NAN;
2567
2568 float nreco = 0;
2569 float phisize = 0;
2570 float zsize = 0;
2571 const auto [reco_ckey, reco_cluster] = clustereval->reco_cluster_from_truth_cluster(ckey, gclus);
2572 if (reco_cluster)
2573 {
2574 nreco = 1;
2575
2576 Acts::Vector3 glob;
2577 glob = tgeometry->getGlobalPosition(reco_ckey, reco_cluster);
2578 x = glob(0);
2579 y = glob(1);
2580 z = glob(2);
2581
2582 TVector3 pos(x, y, z);
2583 r = sqrt(x * x + y * y);
2584 phi = pos.Phi();
2585 eta = pos.Eta();
2586
2587 auto para_errors = ClusErrPara.get_clusterv5_modified_error(reco_cluster, r, ckey);
2588
2589 phisize = reco_cluster->getPhiSize();
2590 zsize = reco_cluster->getZSize();
2591 ez = sqrt(para_errors.second);
2592 ephi = sqrt(para_errors.first);
2593
2594
2595 adc = reco_cluster->getAdc();
2596
2597 if (Verbosity() > 1)
2598 {
2599 std::cout << " reco cluster key " << reco_ckey << " r " << r << " x " << x << " y " << y << " z " << z << " phi " << phi << " adc " << adc << std::endl;
2600 }
2601 }
2602 if (nreco == 0 && Verbosity() > 1)
2603 {
2604 if (Verbosity() > 1)
2605 {
2606 std::cout << " ----------- Failed to find matching reco cluster " << std::endl;
2607 }
2608 }
2609
2610
2611
2612 float g4cluster_data[] = {(float) _ievent,
2613 (float) local_layer,
2614 gx,
2615 gy,
2616 gz,
2617 gt,
2618 gedep,
2619 gr,
2620 gphi,
2621 geta,
2622 gtrackID,
2623 gflavor,
2624 gembed,
2625 gprimary,
2626 gphisize,
2627 gzsize,
2628 gadc,
2629 nreco,
2630 x,
2631 y,
2632 z,
2633 r,
2634 phi,
2635 eta,
2636 ex,
2637 ey,
2638 ez,
2639 ephi,
2640 adc,
2641 phisize,
2642 zsize};
2643 _ntp_g4cluster->Fill(g4cluster_data);
2644 }
2645 }
2646 }
2647
2648
2649
2650
2651
2652
2653
2654
2655
2656 if (_ntp_gtrack)
2657 {
2658 if (Verbosity() > 1)
2659 {
2660 std::cout << "Filling ntp_gtrack " << std::endl;
2661 _timer->restart();
2662 }
2663 ClusterErrorPara ClusErrPara;
2664 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
2665 GlobalVertexMap* gvertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
2666 TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
2667 if(!clustermap)
2668 clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
2669
2670 if (truthinfo)
2671 {
2672 PHG4TruthInfoContainer::ConstRange range = truthinfo->GetParticleRange();
2673 if (_scan_for_primaries)
2674 {
2675 range = truthinfo->GetPrimaryParticleRange();
2676 }
2677
2678 Float_t gntracks = (Float_t) truthinfo->GetNumPrimaryVertexParticles();
2679 Float_t gnchghad = 0;
2680 for(auto iter = range.first; iter!= range.second; ++iter)
2681 {
2682 PHG4Particle* g4particle = iter->second;
2683
2684 if (_scan_for_embedded)
2685 {
2686 if (trutheval->get_embed(g4particle) <= 0)
2687 {
2688 continue;
2689 }
2690 }
2691
2692 float gflavor = g4particle->get_pid();
2693 if(fabs(gflavor)==211 || fabs(gflavor)==321 || fabs(gflavor)==2212)
2694 {
2695 gnchghad++;
2696 }
2697 }
2698
2699 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
2700 iter != range.second;
2701 ++iter)
2702 {
2703 PHG4Particle* g4particle = iter->second;
2704
2705 if (_scan_for_embedded)
2706 {
2707 if (trutheval->get_embed(g4particle) <= 0)
2708 {
2709 continue;
2710 }
2711 }
2712
2713 float gtrackID = g4particle->get_track_id();
2714 float gflavor = g4particle->get_pid();
2715
2716 std::set<TrkrDefs::cluskey> g4clusters = clustereval->all_clusters_from(g4particle);
2717
2718 float ng4hits = g4clusters.size();
2719 unsigned int ngmaps = 0;
2720 unsigned int ngmms = 0;
2721 unsigned int ngintt = 0;
2722 unsigned int ngintt1 = 0;
2723 unsigned int ngintt2 = 0;
2724 unsigned int ngintt3 = 0;
2725 unsigned int ngintt4 = 0;
2726 unsigned int ngintt5 = 0;
2727 unsigned int ngintt6 = 0;
2728 unsigned int ngintt7 = 0;
2729 unsigned int ngintt8 = 0;
2730 unsigned int ngtpc = 0;
2731 unsigned int nglmaps = 0;
2732 unsigned int nglintt = 0;
2733 unsigned int ngltpc = 0;
2734 unsigned int nglmms = 0;
2735
2736 std::vector<int> lmaps(_nlayers_maps + 1,0);
2737 std::vector<int> lintt(_nlayers_intt + 1,0);
2738 std::vector<int> ltpc(_nlayers_tpc + 1,0);
2739 std::vector<int> lmms(_nlayers_mms + 1,0);
2740
2741 for (const TrkrDefs::cluskey g4cluster : g4clusters)
2742 {
2743 unsigned int local_layer = TrkrDefs::getLayer(g4cluster);
2744
2745 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
2746 {
2747 lmaps[local_layer] = 1;
2748 ngmaps++;
2749 }
2750 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
2751 {
2752 lmms[local_layer - _nlayers_tpc - _nlayers_intt - _nlayers_maps] = 1;
2753 ngmms++;
2754 }
2755 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
2756 {
2757 lintt[local_layer - _nlayers_maps] = 1;
2758 ngintt++;
2759 }
2760
2761 if (_nlayers_intt > 0 && local_layer == _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
2762 {
2763 ngintt1++;
2764 }
2765
2766 if (_nlayers_intt > 1 && local_layer == _nlayers_maps + 1 && local_layer < _nlayers_maps + _nlayers_intt)
2767 {
2768 ngintt2++;
2769 }
2770
2771 if (_nlayers_intt > 2 && local_layer == _nlayers_maps + 2 && local_layer < _nlayers_maps + _nlayers_intt)
2772 {
2773 ngintt3++;
2774 }
2775
2776 if (_nlayers_intt > 3 && local_layer == _nlayers_maps + 3 && local_layer < _nlayers_maps + _nlayers_intt)
2777 {
2778 ngintt4++;
2779 }
2780
2781 if (_nlayers_intt > 4 && local_layer == _nlayers_maps + 4 && local_layer < _nlayers_maps + _nlayers_intt)
2782 {
2783 ngintt5++;
2784 }
2785
2786 if (_nlayers_intt > 5 && local_layer == _nlayers_maps + 5 && local_layer < _nlayers_maps + _nlayers_intt)
2787 {
2788 ngintt6++;
2789 }
2790
2791 if (_nlayers_intt > 6 && local_layer == _nlayers_maps + 6 && local_layer < _nlayers_maps + _nlayers_intt)
2792 {
2793 ngintt7++;
2794 }
2795
2796 if (_nlayers_intt > 7 && local_layer == _nlayers_maps + 7 && local_layer < _nlayers_maps + _nlayers_intt)
2797 {
2798 ngintt8++;
2799 }
2800 if (_nlayers_tpc > 0 && local_layer >= _nlayers_maps + _nlayers_intt && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
2801 {
2802 ltpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
2803 ngtpc++;
2804 }
2805 }
2806 if (_nlayers_maps > 0)
2807 {
2808 for (unsigned int i = 0; i < _nlayers_maps; i++)
2809 {
2810 nglmaps += lmaps[i];
2811 }
2812 }
2813 if (_nlayers_intt > 0)
2814 {
2815 for (unsigned int i = 0; i < _nlayers_intt; i++)
2816 {
2817 nglintt += lintt[i];
2818 }
2819 }
2820 if (_nlayers_tpc > 0)
2821 {
2822 for (unsigned int i = 0; i < _nlayers_tpc; i++)
2823 {
2824 ngltpc += ltpc[i];
2825 }
2826 }
2827 if (_nlayers_mms > 0)
2828 {
2829 for (unsigned int i = 0; i < _nlayers_mms; i++)
2830 {
2831 nglmms += lmms[i];
2832 }
2833 }
2834
2835 float gpx = g4particle->get_px();
2836 float gpy = g4particle->get_py();
2837 float gpz = g4particle->get_pz();
2838 float gpt = NAN;
2839 float geta = NAN;
2840 float gphi = NAN;
2841 if (gpx != 0 && gpy != 0)
2842 {
2843 TVector3 gv(gpx, gpy, gpz);
2844 gpt = gv.Pt();
2845 geta = gv.Eta();
2846 gphi = gv.Phi();
2847 }
2848 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
2849 float gvx = vtx->get_x();
2850 float gvy = vtx->get_y();
2851 float gvz = vtx->get_z();
2852 float gvt = vtx->get_t();
2853
2854 float gfpx = 0.;
2855 float gfpy = 0.;
2856 float gfpz = 0.;
2857 float gfx = 0.;
2858 float gfy = 0.;
2859 float gfz = 0.;
2860
2861 PHG4Hit* outerhit = nullptr;
2862 if (_do_eval_light == false)
2863 {
2864 outerhit = trutheval->get_outermost_truth_hit(g4particle);
2865 }
2866
2867 if (outerhit)
2868 {
2869 gfpx = outerhit->get_px(1);
2870 gfpy = outerhit->get_py(1);
2871 gfpz = outerhit->get_pz(1);
2872 gfx = outerhit->get_x(1);
2873 gfy = outerhit->get_y(1);
2874 gfz = outerhit->get_z(1);
2875 }
2876
2877 float gembed = trutheval->get_embed(g4particle);
2878 float gprimary = trutheval->is_primary(g4particle);
2879
2880 float trackID = NAN;
2881 float charge = NAN;
2882 float quality = NAN;
2883 float chisq = NAN;
2884 float ndf = NAN;
2885 float local_nhits = 0;
2886 float nmaps = 0;
2887 float nintt = 0;
2888 float ntpc = 0;
2889 float nmms = 0;
2890 float ntpc1 = 0;
2891 float ntpc11 = 0;
2892 float ntpc2 = 0;
2893 float ntpc3 = 0;
2894 float nlintt = 0;
2895 float nlmaps = 0;
2896 float nltpc = 0;
2897 float nlmms = 0;
2898 unsigned int layers = 0x0;
2899 int vertexID = -1;
2900 float vx = NAN;
2901 float vy = NAN;
2902 float vz = NAN;
2903 float dca2d = NAN;
2904 float dca2dsigma = NAN;
2905 float dca3dxy = NAN;
2906 float dca3dxysigma = NAN;
2907 float dca3dz = NAN;
2908 float dca3dzsigma = NAN;
2909 float px = NAN;
2910 float py = NAN;
2911 float pz = NAN;
2912 float pt = NAN;
2913 float eta = NAN;
2914 float phi = NAN;
2915 float deltapt = NAN;
2916 float deltaeta = NAN;
2917 float deltaphi = NAN;
2918 float pcax = NAN;
2919 float pcay = NAN;
2920 float pcaz = NAN;
2921
2922 float nfromtruth = NAN;
2923 float nwrong = NAN;
2924 float ntrumaps = NAN;
2925 float nwrongmaps = NAN;
2926 float ntruintt = NAN;
2927 float nwrongintt = NAN;
2928 float ntrumms = NAN;
2929 float nwrongmms = NAN;
2930 float ntrutpc = NAN;
2931 float nwrongtpc = NAN;
2932 float ntrutpc1 = NAN;
2933 float nwrongtpc1 = NAN;
2934 float ntrutpc11 = NAN;
2935 float nwrongtpc11 =NAN;
2936 float ntrutpc2 = NAN;
2937 float nwrongtpc2 = NAN;
2938 float ntrutpc3 = NAN;
2939 float nwrongtpc3 = NAN;
2940 float layersfromtruth = NAN;
2941 float npedge = 0;
2942 float nredge = 0;
2943 float nbig = 0;
2944 float novlp = 0;
2945 float merr = 0;
2946 float msize = 0;
2947 float siqr = NAN;
2948 float siphi = NAN;
2949 float sithe = NAN;
2950 float six0 = NAN;
2951 float siy0 = NAN;
2952 float tpqr = NAN;
2953 float tpphi = NAN;
2954 float tpthe = NAN;
2955 float tpx0 = NAN;
2956 float tpy0 = NAN;
2957
2958 if (_do_track_match)
2959 {
2960 SvtxTrack* track = trackeval->best_track_from(g4particle);
2961
2962 if (track)
2963 {
2964 trackID = track->get_id();
2965 charge = track->get_charge();
2966 quality = track->get_quality();
2967 chisq = track->get_chisq();
2968 ndf = track->get_ndf();
2969 TrackSeed* silseed = track->get_silicon_seed();
2970 TrackSeed* tpcseed = track->get_tpc_seed();
2971 if (tpcseed)
2972 {
2973 local_nhits += tpcseed->size_cluster_keys();
2974 }
2975 if (silseed)
2976 {
2977 local_nhits += silseed->size_cluster_keys();
2978 }
2979
2980 std::vector<int> maps(_nlayers_maps+1, 0);
2981 std::vector<int> intt(_nlayers_intt+1, 0);
2982 std::vector<int> tpc(_nlayers_tpc+1, 0);
2983 std::vector<int> mms(_nlayers_mms+1, 0);
2984
2985
2986 if (tpcseed)
2987 {
2988 tpqr = tpcseed->get_qOverR();
2989 tpphi = tpcseed->get_phi();
2990 tpthe = tpcseed->get_theta();
2991 tpx0 = tpcseed->get_X0();
2992 tpy0 = tpcseed->get_Y0();
2993 for (TrackSeed::ConstClusterKeyIter local_iter = tpcseed->begin_cluster_keys();
2994 local_iter != tpcseed->end_cluster_keys();
2995 ++local_iter)
2996 {
2997 TrkrDefs::cluskey cluster_key = *local_iter;
2998 TrkrCluster* cluster = clustermap->findCluster(cluster_key);
2999 unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
3000 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
3001 {
3002 maps[local_layer] = 1;
3003 nmaps++;
3004 }
3005 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
3006 {
3007 intt[local_layer - _nlayers_maps] = 1;
3008 nintt++;
3009 }
3010 if (_nlayers_tpc > 0 &&
3011 local_layer >= (_nlayers_maps + _nlayers_intt) &&
3012 local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
3013 {
3014 tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
3015 ntpc++;
3016
3017 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 8)
3018 {
3019 ntpc11++;
3020 }
3021
3022 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 16)
3023 {
3024
3025 ntpc1++;
3026 }
3027 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 32)
3028 {
3029
3030 ntpc2++;
3031 }
3032 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 48)
3033 {
3034
3035 ntpc3++;
3036 }
3037 }
3038
3039 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
3040 {
3041 mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3042 nmms++;
3043 }
3044
3045 Acts::Vector3 glob;
3046 glob = tgeometry->getGlobalPosition(cluster_key, cluster);
3047 float x = glob(0);
3048 float y = glob(1);
3049 float z = glob(2);
3050
3051 TVector3 pos(x, y, z);
3052 float r = sqrt(x*x+y*y);
3053 phi = pos.Phi();
3054 eta = pos.Eta();
3055 float gphisize = 0;
3056
3057 float gphierr = 0;
3058
3059 float govlp = 0 ;
3060 float gedge = 0;
3061 if(cluster!=nullptr){
3062
3063 gphisize = cluster->getPhiSize();
3064
3065 auto para_errors = ClusErrPara.get_clusterv5_modified_error(cluster,r,cluster_key);
3066
3067 gphierr = sqrt(para_errors.first);
3068 govlp = cluster->getOverlap();
3069 gedge = cluster->getEdge();
3070
3071 if(gedge>0) npedge++;
3072 if(gphisize>=4) nbig++;
3073 if(govlp>=2) novlp++;
3074 merr+=gphierr;
3075 msize+=gphisize;
3076 }
3077 if(local_layer==7||local_layer==22||local_layer==23||local_layer==38||local_layer==39) nredge++;
3078 }
3079 }
3080
3081 if (silseed)
3082 {
3083 siqr = silseed->get_qOverR();
3084 siphi = silseed->get_phi();
3085 sithe = silseed->get_theta();
3086 six0 = silseed->get_X0();
3087 siy0 = silseed->get_Y0();
3088 for (TrackSeed::ConstClusterKeyIter local_iter = silseed->begin_cluster_keys();
3089 local_iter != silseed->end_cluster_keys();
3090 ++local_iter)
3091 {
3092 TrkrDefs::cluskey cluster_key = *local_iter;
3093
3094 unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
3095 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
3096 {
3097 maps[local_layer] = 1;
3098 nmaps++;
3099 }
3100 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
3101 {
3102 intt[local_layer - _nlayers_maps] = 1;
3103 nintt++;
3104 }
3105 if (_nlayers_tpc > 0 &&
3106 local_layer >= (_nlayers_maps + _nlayers_intt) &&
3107 local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
3108 {
3109 tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
3110 ntpc++;
3111
3112 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 8)
3113 {
3114 ntpc11++;
3115 }
3116
3117 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 16)
3118 {
3119
3120 ntpc1++;
3121 }
3122 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 32)
3123 {
3124
3125 ntpc2++;
3126 }
3127 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 48)
3128 {
3129
3130 ntpc3++;
3131 }
3132 }
3133
3134 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
3135 {
3136 mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3137 nmms++;
3138 }
3139 }
3140 }
3141
3142 if (_nlayers_maps > 0)
3143 {
3144 for (unsigned int i = 0; i < _nlayers_maps; i++)
3145 {
3146 nlmaps += maps[i];
3147 }
3148 }
3149 if (_nlayers_intt > 0)
3150 {
3151 for (unsigned int i = 0; i < _nlayers_intt; i++)
3152 {
3153 nlintt += intt[i];
3154 }
3155 }
3156 if (_nlayers_tpc > 0)
3157 {
3158 for (unsigned int i = 0; i < _nlayers_tpc; i++)
3159 {
3160 nltpc += tpc[i];
3161 }
3162 }
3163 if (_nlayers_mms > 0)
3164 {
3165 for (unsigned int i = 0; i < _nlayers_mms; i++)
3166 {
3167 nlmms += mms[i];
3168 }
3169 }
3170
3171 layers = nlmaps + nlintt + nltpc + nlmms;
3172
3173
3174
3175
3176
3177
3178
3179
3180
3181
3182
3183 vertexID = track->get_vertex_id();
3184
3185 GlobalVertex* vertex = gvertexmap->get(vertexID);
3186 if (vertex)
3187 {
3188 vx = vertex->get_x();
3189 vy = vertex->get_y();
3190 vz = vertex->get_z();
3191 Acts::Vector3 vert(vx,vy,vz);
3192 auto dcapair = TrackAnalysisUtils::get_dca(track, vert);
3193 dca3dxy = dcapair.first.first;
3194 dca3dxysigma = dcapair.first.second;
3195 dca3dz = dcapair.second.first;
3196 dca3dzsigma = dcapair.second.second;
3197 }
3198
3199 px = track->get_px();
3200 py = track->get_py();
3201 pz = track->get_pz();
3202 TVector3 v(px, py, pz);
3203 pt = v.Pt();
3204 eta = v.Eta();
3205 phi = v.Phi();
3206 float CVxx = track->get_error(3, 3);
3207 float CVxy = track->get_error(3, 4);
3208 float CVxz = track->get_error(3, 5);
3209 float CVyy = track->get_error(4, 4);
3210 float CVyz = track->get_error(4, 5);
3211 float CVzz = track->get_error(5, 5);
3212 deltapt = sqrt((CVxx * px * px + 2 * CVxy * px * py + CVyy * py * py) / (px * px + py * py));
3213 deltaeta = sqrt((CVzz * (px * px + py * py) * (px * px + py * py) + pz * (-2 * (CVxz * px + CVyz * py) * (px * px + py * py) + CVxx * px * px * pz + CVyy * py * py * pz + 2 * CVxy * px * py * pz)) / ((px * px + py * py) * (px * px + py * py) * (px * px + py * py + pz * pz)));
3214 deltaphi = sqrt((CVyy * px * px - 2 * CVxy * px * py + CVxx * py * py) / ((px * px + py * py) * (px * px + py * py)));
3215 pcax = track->get_x();
3216 pcay = track->get_y();
3217 pcaz = track->get_z();
3218
3219 nfromtruth = trackeval->get_nclusters_contribution(track, g4particle);
3220 nwrong = trackeval->get_nwrongclusters_contribution(track, g4particle);
3221
3222 if (_nlayers_maps == 0)
3223 {
3224 ntrumaps = 0;
3225 }
3226 else
3227 {
3228 auto pair = trackeval->get_layer_range_contribution(track, g4particle, 0, _nlayers_maps);
3229 ntrumaps = pair.first;
3230 nwrongmaps = pair.second;
3231 }
3232 if (_nlayers_intt == 0)
3233 {
3234 ntruintt = 0;
3235 }
3236 else
3237 {
3238 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps, _nlayers_maps + _nlayers_intt);
3239 ntruintt = pair.first;
3240 nwrongintt = pair.second;
3241 }
3242 if (_nlayers_mms == 0)
3243 {
3244 ntrumms = 0;
3245 }
3246 else
3247 {
3248 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + _nlayers_tpc, _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms);
3249 ntrumms = pair.first;
3250 nwrongmms = pair.second;
3251 }
3252 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
3253 ntrutpc = pair.first;
3254 nwrongtpc = pair.second;
3255 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 16);
3256 ntrutpc1 = pair.first;
3257 nwrongtpc1 = pair.second;
3258 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 8);
3259 ntrutpc11 = pair.first;
3260 nwrongtpc11 = pair.second;
3261 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + 16, _nlayers_maps + _nlayers_intt + 32);
3262 ntrutpc2 = pair.first;
3263 nwrongtpc2 = pair.second;
3264 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + 32, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
3265 ntrutpc3 = pair.first;
3266 nwrongtpc3 = pair.first;
3267 layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track, g4particle);
3268 }
3269 }
3270
3271 float gtrack_data[] = {(float) _ievent, m_fSeed,
3272 gntracks,
3273 gnchghad,
3274 gtrackID,
3275 gflavor,
3276 ng4hits,
3277 (float) ngmaps,
3278 (float) ngintt,
3279 (float) ngmms,
3280 (float) ngintt1,
3281 (float) ngintt2,
3282 (float) ngintt3,
3283 (float) ngintt4,
3284 (float) ngintt5,
3285 (float) ngintt6,
3286 (float) ngintt7,
3287 (float) ngintt8,
3288 (float) ngtpc,
3289 (float) nglmaps,
3290 (float) nglintt,
3291 (float) ngltpc,
3292 (float) nglmms,
3293 gpx,
3294 gpy,
3295 gpz,
3296 gpt,
3297 geta,
3298 gphi,
3299 gvx,
3300 gvy,
3301 gvz,
3302 gvt,
3303 gfpx,
3304 gfpy,
3305 gfpz,
3306 gfx,
3307 gfy,
3308 gfz,
3309 gembed,
3310 gprimary,
3311 trackID,
3312 px,
3313 py,
3314 pz,
3315 pt,
3316 eta,
3317 phi,
3318 deltapt,
3319 deltaeta,
3320 deltaphi,
3321 siqr,
3322 siphi,
3323 sithe,
3324 six0,
3325 siy0,
3326 tpqr,
3327 tpphi,
3328 tpthe,
3329 tpx0,
3330 tpy0,
3331 charge,
3332 quality,
3333 chisq,
3334 ndf,
3335 local_nhits,
3336 (float) layers,
3337 nmaps,
3338 nintt,
3339 ntpc,
3340 nmms,
3341 ntpc1,
3342 ntpc11,
3343 ntpc2,
3344 ntpc3,
3345 nlmaps,
3346 nlintt,
3347 nltpc,
3348 nlmms,
3349 (float) vertexID,
3350 vx,
3351 vy,
3352 vz,
3353 dca2d,
3354 dca2dsigma,
3355 dca3dxy,
3356 dca3dxysigma,
3357 dca3dz,
3358 dca3dzsigma,
3359 pcax,
3360 pcay,
3361 pcaz,
3362 nfromtruth,
3363 nwrong,
3364 ntrumaps,
3365 nwrongmaps,
3366 ntruintt,
3367 nwrongintt,
3368 ntrutpc,
3369 nwrongtpc,
3370 ntrumms,
3371 nwrongmms,
3372 ntrutpc1,
3373 nwrongtpc1,
3374 ntrutpc11,
3375 nwrongtpc11,
3376 ntrutpc2,
3377 nwrongtpc2,
3378 ntrutpc3,
3379 nwrongtpc3,
3380 layersfromtruth,
3381 npedge,
3382 nredge,
3383 nbig,
3384 novlp,
3385 merr,
3386 msize,
3387 nhit_tpc_all,
3388 nhit_tpc_in,
3389 nhit_tpc_mid,
3390 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
3391
3392
3393
3394
3395
3396
3397
3398
3399
3400 _ntp_gtrack->Fill(gtrack_data);
3401 }
3402 }
3403 if (Verbosity() > 1)
3404 {
3405 _timer->stop();
3406 std::cout << "gtrack time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
3407 }
3408 }
3409
3410
3411
3412
3413
3414 if (_ntp_track)
3415 {
3416 if (Verbosity() > 1)
3417 {
3418 std::cout << "Filling ntp_track " << std::endl;
3419 _timer->restart();
3420 }
3421
3422
3423 SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname.c_str());
3424
3425 GlobalVertexMap* gvertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
3426 TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
3427 if(!clustermap)
3428 clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
3429 ClusterErrorPara ClusErrPara;
3430 if (trackmap)
3431 {
3432 for (auto& iter : *trackmap)
3433 {
3434 SvtxTrack* track = iter.second;
3435 float trackID = track->get_id();
3436
3437 TrackSeed* tpcseed = track->get_tpc_seed();
3438 TrackSeed* silseed = track->get_silicon_seed();
3439 short int crossing_int = track->get_crossing();
3440 float crossing;
3441 if (crossing_int == SHRT_MAX)
3442 {
3443 crossing = NAN;
3444 }
3445 else
3446 {
3447 crossing = (float) crossing_int;
3448 }
3449 float charge = track->get_charge();
3450 float quality = track->get_quality();
3451 float chisq = track->get_chisq();
3452 float ndf = track->get_ndf();
3453 float local_nhits = 0;
3454 if (tpcseed)
3455 {
3456 local_nhits += tpcseed->size_cluster_keys();
3457 }
3458 if (silseed)
3459 {
3460 local_nhits += silseed->size_cluster_keys();
3461 }
3462 unsigned int layers = 0x0;
3463
3464 std::vector<int> maps(_nlayers_maps+1,0);
3465 std::vector<int> intt(_nlayers_intt+1,0);
3466 std::vector<int> tpc(_nlayers_tpc+1,0);
3467 std::vector<int> mms(_nlayers_mms+1,0);
3468
3469 float nmaps = 0;
3470 float nintt = 0;
3471 float nmms = 0;
3472 float ntpc = 0;
3473 float ntpc1 = 0;
3474 float ntpc11 = 0;
3475 float ntpc2 = 0;
3476 float ntpc3 = 0;
3477 float nlmaps = 0;
3478 float nlintt = 0;
3479 float nltpc = 0;
3480 float nlmms = 0;
3481 float npedge = 0;
3482 float nredge = 0;
3483 float nbig = 0;
3484 float novlp = 0;
3485 float merr = 0;
3486 float msize = 0;
3487 float siqr = NAN;
3488 float siphi = NAN;
3489 float sithe = NAN;
3490 float six0 = NAN;
3491 float siy0 = NAN;
3492 float tpqr = NAN;
3493 float tpphi = NAN;
3494 float tpthe = NAN;
3495 float tpx0 = NAN;
3496 float tpy0 = NAN;
3497
3498 if (tpcseed)
3499 {
3500 tpqr = tpcseed->get_qOverR();
3501 tpphi = tpcseed->get_phi();
3502 tpthe = tpcseed->get_theta();
3503 tpx0 = tpcseed->get_X0();
3504 tpy0 = tpcseed->get_Y0();
3505 for (SvtxTrack::ConstClusterKeyIter local_iter = tpcseed->begin_cluster_keys();
3506 local_iter != tpcseed->end_cluster_keys();
3507 ++local_iter){
3508 TrkrDefs::cluskey cluster_key = *local_iter;
3509 TrkrCluster* cluster = clustermap->findCluster(cluster_key);
3510 unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
3511
3512 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
3513 {
3514 maps[local_layer] = 1;
3515 nmaps++;
3516 }
3517 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
3518 {
3519 intt[local_layer - _nlayers_maps] = 1;
3520 nintt++;
3521 }
3522 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
3523 {
3524 mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3525 nmms++;
3526 }
3527 if (_nlayers_tpc > 0 && local_layer >= (_nlayers_maps + _nlayers_intt) && local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
3528 {
3529 tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
3530 ntpc++;
3531 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 8)
3532 {
3533 ntpc11++;
3534 }
3535
3536 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 16)
3537 {
3538 ntpc1++;
3539 }
3540 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 32)
3541 {
3542 ntpc2++;
3543 }
3544 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 48)
3545 {
3546 ntpc3++;
3547 }
3548 }
3549
3550 Acts::Vector3 glob;
3551 glob = tgeometry->getGlobalPosition(cluster_key, cluster);
3552 float x = glob(0);
3553 float y = glob(1);
3554 float z = glob(2);
3555
3556 TVector3 pos(x, y, z);
3557 float r = sqrt(x*x+y*y);
3558
3559
3560 float rphisize = 0;
3561
3562 float rphierr = 0;
3563
3564 float rovlp = 0 ;
3565 float pedge = 0;
3566 if(cluster!=nullptr){
3567
3568 rphisize = cluster->getPhiSize();
3569
3570 auto para_errors = ClusErrPara.get_clusterv5_modified_error(cluster,r,cluster_key);
3571
3572 rphierr = sqrt(para_errors.first);
3573 rovlp = cluster->getOverlap();
3574 pedge = cluster->getEdge();
3575
3576 if(pedge>0) npedge++;
3577 if(rphisize>=4) nbig++;
3578 if(rovlp>=2) novlp++;
3579 merr+=rphierr;
3580 msize+=rphisize;
3581 }
3582 if(local_layer==7||local_layer==22||local_layer==23||local_layer==38||local_layer==39) nredge++;
3583 if(Verbosity() > 2)
3584 {
3585 std::cout << " lay: " << local_layer
3586 << " pedge " << pedge
3587 << " | " << npedge
3588 << " nredge " << nredge
3589 << " rphisize " << rphisize
3590 << " | " << nbig
3591 << " rovlp " << rovlp
3592 << " | " << novlp
3593 << std::endl;
3594 }
3595 }
3596 }
3597
3598 if (silseed)
3599 {
3600 siqr = silseed->get_qOverR();
3601 siphi = silseed->get_phi();
3602 sithe = silseed->get_theta();
3603 six0 = silseed->get_X0();
3604 siy0 = silseed->get_Y0();
3605 for (SvtxTrack::ConstClusterKeyIter local_iter = silseed->begin_cluster_keys();
3606 local_iter != silseed->end_cluster_keys();
3607 ++local_iter)
3608 {
3609 TrkrDefs::cluskey cluster_key = *local_iter;
3610
3611 unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
3612
3613 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
3614 {
3615 maps[local_layer] = 1;
3616 nmaps++;
3617 }
3618 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
3619 {
3620 intt[local_layer - _nlayers_maps] = 1;
3621 nintt++;
3622 }
3623 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
3624 {
3625 mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3626 nmms++;
3627 }
3628 if (_nlayers_tpc > 0 && local_layer >= (_nlayers_maps + _nlayers_intt) && local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
3629 {
3630 tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
3631 ntpc++;
3632 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 8)
3633 {
3634 ntpc11++;
3635 }
3636
3637 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 16)
3638 {
3639 ntpc1++;
3640 }
3641 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 32)
3642 {
3643 ntpc2++;
3644 }
3645 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 48)
3646 {
3647 ntpc3++;
3648 }
3649 }
3650 }
3651 }
3652
3653 if (_nlayers_maps > 0)
3654 {
3655 for (unsigned int i = 0; i < _nlayers_maps; i++)
3656 {
3657 nlmaps += maps[i];
3658 }
3659 }
3660 if (_nlayers_intt > 0)
3661 {
3662 for (unsigned int i = 0; i < _nlayers_intt; i++)
3663 {
3664 nlintt += intt[i];
3665 }
3666 }
3667 if (_nlayers_tpc > 0)
3668 {
3669 for (unsigned int i = 0; i < _nlayers_tpc; i++)
3670 {
3671 nltpc += tpc[i];
3672 }
3673 }
3674 if (_nlayers_mms > 0)
3675 {
3676 for (unsigned int i = 0; i < _nlayers_mms; i++)
3677 {
3678 nlmms += mms[i];
3679 }
3680 }
3681 layers = nlmaps + nlintt + nltpc + nlmms;
3682
3683 float dca3dxy = NAN, dca3dz = NAN,
3684 dca3dxysigma = NAN, dca3dzsigma = NAN;
3685 float dca2d = NAN, dca2dsigma = NAN;
3686
3687
3688 int vertexID = track->get_vertex_id();
3689 GlobalVertex* vertex = gvertexmap->get(vertexID);
3690 float vx = NAN;
3691 float vy = NAN;
3692 float vz = NAN;
3693 if (vertex)
3694 {
3695 vx = vertex->get_x();
3696 vy = vertex->get_y();
3697 vz = vertex->get_z();
3698 Acts::Vector3 vert(vx,vy,vz);
3699
3700 auto dcapair = TrackAnalysisUtils::get_dca(track, vert);
3701 dca3dxy = dcapair.first.first;
3702 dca3dxysigma = dcapair.first.second;
3703 dca3dz = dcapair.second.first;
3704 dca3dzsigma = dcapair.second.second;
3705 }
3706 float px = track->get_px();
3707 float py = track->get_py();
3708 float pz = track->get_pz();
3709 TVector3 v(px, py, pz);
3710 float pt = v.Pt();
3711 float eta = v.Eta();
3712 float phi = v.Phi();
3713 float CVxx = track->get_error(3, 3);
3714 float CVxy = track->get_error(3, 4);
3715 float CVxz = track->get_error(3, 5);
3716 float CVyy = track->get_error(4, 4);
3717 float CVyz = track->get_error(4, 5);
3718 float CVzz = track->get_error(5, 5);
3719 float deltapt = sqrt((CVxx * px * px + 2 * CVxy * px * py + CVyy * py * py) / (px * px + py * py));
3720 float deltaeta = sqrt((CVzz * (px * px + py * py) * (px * px + py * py) + pz * (-2 * (CVxz * px + CVyz * py) * (px * px + py * py) + CVxx * px * px * pz + CVyy * py * py * pz + 2 * CVxy * px * py * pz)) / ((px * px + py * py) * (px * px + py * py) * (px * px + py * py + pz * pz)));
3721 float deltaphi = sqrt((CVyy * px * px - 2 * CVxy * px * py + CVxx * py * py) / ((px * px + py * py) * (px * px + py * py)));
3722 float pcax = track->get_x();
3723 float pcay = track->get_y();
3724 float pcaz = track->get_z();
3725
3726 float gtrackID = NAN;
3727 float gflavor = NAN;
3728 float ng4hits = NAN;
3729 unsigned int ngmaps = 0;
3730 unsigned int ngintt = 0;
3731 unsigned int ngmms = 0;
3732 unsigned int ngtpc = 0;
3733 unsigned int nglmaps = 0;
3734 unsigned int nglintt = 0;
3735 unsigned int nglmms = 0;
3736 unsigned int ngltpc = 0;
3737 float gpx = NAN;
3738 float gpy = NAN;
3739 float gpt = NAN;
3740 float geta = NAN;
3741 float gphi = NAN;
3742 float gpz = NAN;
3743 float gvx = NAN;
3744 float gvy = NAN;
3745 float gvz = NAN;
3746 float gvt = NAN;
3747 float gfpx = NAN;
3748 float gfpy = NAN;
3749 float gfpz = NAN;
3750 float gfx = NAN;
3751 float gfy = NAN;
3752 float gfz = NAN;
3753 float gembed = NAN;
3754 float gprimary = NAN;
3755
3756 int ispure = 0;
3757 float nfromtruth = NAN;
3758 float nwrong = NAN;
3759 float ntrumaps = NAN;
3760 float nwrongmaps = NAN;
3761 float ntruintt = NAN;
3762 float nwrongintt = NAN;
3763 float ntrumms = NAN;
3764 float nwrongmms = NAN;
3765 float ntrutpc = NAN;
3766 float nwrongtpc = NAN;
3767 float ntrutpc1 = NAN;
3768 float nwrongtpc1 = NAN;
3769 float ntrutpc11 = NAN;
3770 float nwrongtpc11 = NAN;
3771 float ntrutpc2 = NAN;
3772 float nwrongtpc2 = NAN;
3773 float ntrutpc3 = NAN;
3774 float nwrongtpc3 = NAN;
3775 float layersfromtruth = NAN;
3776
3777 if (_do_track_match)
3778 {
3779 PHG4Particle* g4particle = trackeval->max_truth_particle_by_nclusters(track);
3780 if (g4particle)
3781 {
3782 if (_scan_for_embedded)
3783 {
3784 if (trutheval->get_embed(g4particle) <= 0)
3785 {
3786 continue;
3787 }
3788 }
3789 SvtxTrack* truthrecotrk = trackeval->best_track_from(g4particle);
3790 if(truthrecotrk)
3791 {
3792 if(truthrecotrk->get_id() == track->get_id())
3793 {
3794 ispure = 1;
3795 }
3796 }
3797 gtrackID = g4particle->get_track_id();
3798 gflavor = g4particle->get_pid();
3799
3800 std::set<TrkrDefs::cluskey> g4clusters = clustereval->all_clusters_from(g4particle);
3801 ng4hits = g4clusters.size();
3802 gpx = g4particle->get_px();
3803 gpy = g4particle->get_py();
3804 gpz = g4particle->get_pz();
3805
3806 std::vector<int> lmaps(_nlayers_maps + 1,0);
3807 std::vector<int> lintt(_nlayers_intt + 1,0);
3808 std::vector<int> ltpc(_nlayers_tpc + 1,0);
3809 std::vector<int> lmms(_nlayers_mms + 1,0);
3810
3811 for (const TrkrDefs::cluskey g4cluster : g4clusters)
3812 {
3813 unsigned int local_layer = TrkrDefs::getLayer(g4cluster);
3814 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
3815 {
3816 lmaps[local_layer] = 1;
3817 ngmaps++;
3818 }
3819
3820 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
3821 {
3822 lintt[local_layer - _nlayers_maps] = 1;
3823 ngintt++;
3824 }
3825
3826 if (_nlayers_tpc > 0 && local_layer >= _nlayers_maps + _nlayers_intt && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
3827 {
3828 ltpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
3829 ngtpc++;
3830 }
3831
3832 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
3833 {
3834 lmms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3835 ngmms++;
3836 }
3837 }
3838 if (_nlayers_maps > 0)
3839 {
3840 for (unsigned int i = 0; i < _nlayers_maps; i++)
3841 {
3842 nglmaps += lmaps[i];
3843 }
3844 }
3845 if (_nlayers_intt > 0)
3846 {
3847 for (unsigned int i = 0; i < _nlayers_intt; i++)
3848 {
3849 nglintt += lintt[i];
3850 }
3851 }
3852 if (_nlayers_tpc > 0)
3853 {
3854 for (unsigned int i = 0; i < _nlayers_tpc; i++)
3855 {
3856 ngltpc += ltpc[i];
3857 }
3858 }
3859 if (_nlayers_mms > 0)
3860 {
3861 for (unsigned int i = 0; i < _nlayers_mms; i++)
3862 {
3863 nglmms += lmms[i];
3864 }
3865 }
3866
3867 TVector3 gv(gpx, gpy, gpz);
3868 gpt = gv.Pt();
3869 geta = gv.Eta();
3870 gphi = gv.Phi();
3871 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
3872 gvx = vtx->get_x();
3873 gvy = vtx->get_y();
3874 gvz = vtx->get_z();
3875 gvt = vtx->get_t();
3876
3877 PHG4Hit* outerhit = nullptr;
3878 if (_do_eval_light == false)
3879 {
3880 outerhit = trutheval->get_outermost_truth_hit(g4particle);
3881 }
3882 if (outerhit)
3883 {
3884 gfpx = outerhit->get_px(1);
3885 gfpy = outerhit->get_py(1);
3886 gfpz = outerhit->get_pz(1);
3887 gfx = outerhit->get_x(1);
3888 gfy = outerhit->get_y(1);
3889 gfz = outerhit->get_z(1);
3890 }
3891 gembed = trutheval->get_embed(g4particle);
3892 gprimary = trutheval->is_primary(g4particle);
3893
3894 nfromtruth = trackeval->get_nclusters_contribution(track, g4particle);
3895 nwrong = trackeval->get_nwrongclusters_contribution(track, g4particle);
3896 if (_nlayers_maps == 0)
3897 {
3898 ntrumaps = 0;
3899 }
3900 else
3901 {
3902 auto pair = trackeval->get_layer_range_contribution(track, g4particle, 0, _nlayers_maps);
3903 ntrumaps = pair.first;
3904 nwrongmaps = pair.second;
3905 }
3906 if (_nlayers_intt == 0)
3907 {
3908 ntruintt = 0;
3909 }
3910 else
3911 {
3912 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps, _nlayers_maps + _nlayers_intt);
3913 ntruintt = pair.first;
3914 nwrongintt = pair.second;
3915 }
3916 if (_nlayers_mms == 0)
3917 {
3918 ntrumms = 0;
3919 }
3920 else
3921 {
3922 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + _nlayers_tpc, _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms);
3923 ntrumms = pair.first;
3924 nwrongmms = pair.second;
3925 }
3926 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
3927 ntrutpc = pair.first;
3928 nwrongtpc = pair.second;
3929 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 16);
3930 ntrutpc1 = pair.first;
3931 nwrongtpc1 = pair.second;
3932 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 8);
3933 ntrutpc11 = pair.first;
3934 nwrongtpc11 = pair.second;
3935 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + 16, _nlayers_maps + _nlayers_intt + 32);
3936 ntrutpc2 = pair.first;
3937 nwrongtpc2 = pair.second;
3938 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + 32, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
3939 ntrutpc3 = pair.first;
3940 nwrongtpc3 = pair.second;
3941 layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track, g4particle);
3942 }
3943 }
3944 if(Verbosity() > 2)
3945 {
3946 std::cout << " npedge " << npedge
3947 << " nredge " << nredge
3948 << " nbig " << nbig
3949 << " novlp "<< novlp
3950 << std::endl;
3951 }
3952 float track_data[] = {(float) _ievent, m_fSeed,
3953 trackID,
3954 crossing,
3955 px,
3956 py,
3957 pz,
3958 pt,
3959 eta,
3960 phi,
3961 deltapt,
3962 deltaeta,
3963 deltaphi,
3964 siqr,
3965 siphi,
3966 sithe,
3967 six0,
3968 siy0,
3969 tpqr,
3970 tpphi,
3971 tpthe,
3972 tpx0,
3973 tpy0,
3974 charge,
3975 quality,
3976 chisq,
3977 ndf,
3978 local_nhits, nmaps, nintt, ntpc, nmms,
3979 ntpc1, ntpc11, ntpc2, ntpc3,
3980 nlmaps, nlintt, nltpc, nlmms,
3981 (float) layers,
3982 (float) vertexID,
3983 vx,
3984 vy,
3985 vz,
3986 dca2d,
3987 dca2dsigma,
3988 dca3dxy,
3989 dca3dxysigma,
3990 dca3dz,
3991 dca3dzsigma,
3992 pcax,
3993 pcay,
3994 pcaz,
3995 gtrackID,
3996 (float)ispure,
3997 gflavor,
3998 ng4hits,
3999 (float) ngmaps,
4000 (float) ngintt,
4001 (float) ngtpc,
4002 (float) ngmms,
4003 (float) nglmaps,
4004 (float) nglintt,
4005 (float) ngltpc,
4006 (float) nglmms,
4007 gpx,
4008 gpy,
4009 gpz,
4010 gpt,
4011 geta,
4012 gphi,
4013 gvx,
4014 gvy,
4015 gvz,
4016 gvt,
4017 gfpx,
4018 gfpy,
4019 gfpz,
4020 gfx,
4021 gfy,
4022 gfz,
4023 gembed,
4024 gprimary,
4025 nfromtruth,
4026 nwrong,
4027 ntrumaps,
4028 nwrongmaps,
4029 ntruintt,
4030 nwrongintt,
4031 ntrutpc,
4032 nwrongtpc,
4033 ntrumms,
4034 nwrongmms,
4035 ntrutpc1,
4036 nwrongtpc1,
4037 ntrutpc11,
4038 nwrongtpc11,
4039 ntrutpc2,
4040 nwrongtpc2,
4041 ntrutpc3,
4042 nwrongtpc3,
4043 layersfromtruth,
4044 npedge,
4045 nredge,
4046 nbig,
4047 novlp,
4048 merr,
4049 msize,
4050 nhit_tpc_all,
4051 nhit_tpc_in,
4052 nhit_tpc_mid,
4053 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
4054
4055 if (Verbosity() >= 1)
4056 {
4057 std::cout << "ievent " << _ievent
4058 << " trackID " << trackID
4059 << " nhits " << local_nhits
4060 << " px " << px
4061 << " py " << py
4062 << " pz " << pz
4063 << " gembed " << gembed
4064 << " gprimary " << gprimary
4065 << std::endl;
4066 }
4067
4068 _ntp_track->Fill(track_data);
4069 }
4070 }
4071 if (Verbosity() > 1)
4072 {
4073 _timer->stop();
4074 std::cout << "track time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
4075 }
4076 }
4077
4078
4079
4080
4081
4082 if (_ntp_gseed)
4083 {
4084 if (Verbosity() > 1)
4085 {
4086 std::cout << "Filling ntp_gseed " << std::endl;
4087 _timer->restart();
4088 }
4089
4090 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
4091
4092 float gx = NAN;
4093 float gy = NAN;
4094 float gz = NAN;
4095 float gr = NAN;
4096 float geta = NAN;
4097 float gphi = NAN;
4098 float glayer = NAN;
4099 float gpx = NAN;
4100 float gpy = NAN;
4101 float gpz = NAN;
4102 float gtpt = NAN;
4103 float gtphi = NAN;
4104 float gteta = NAN;
4105 float gvx = NAN;
4106 float gvy = NAN;
4107 float gvz = NAN;
4108 float gembed = NAN;
4109 float gprimary = NAN;
4110 float gflav = NAN;
4111 float dphiprev = NAN;
4112 float detaprev = NAN;
4113
4114 float *xval = new float[_nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms];
4115 float *yval = new float[_nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms];
4116 float *zval = new float[_nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms];
4117 if (truthinfo)
4118 {
4119 int ntrk = 0;
4120 PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
4121 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
4122 iter != range.second;
4123 ++iter)
4124 {
4125 ntrk++;
4126 PHG4Particle* g4particle = iter->second;
4127 for (unsigned int i = 0; i < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms; i++)
4128 {
4129 xval[i] = 0;
4130 yval[i] = 0;
4131 zval[i] = 0;
4132 }
4133 std::set<PHG4Hit*> truth_hits = trutheval->all_truth_hits(g4particle);
4134 for (auto g4hit : truth_hits)
4135 {
4136 unsigned int local_layer = g4hit->get_layer();
4137
4138 if (local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
4139 {
4140
4141 continue;
4142 }
4143 xval[local_layer] = g4hit->get_avg_x();
4144 yval[local_layer] = g4hit->get_avg_y();
4145 zval[local_layer] = g4hit->get_avg_z();
4146 }
4147
4148 for (unsigned int i = 0; i < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms; i++)
4149 {
4150 gx = xval[i];
4151 gy = yval[i];
4152 gz = zval[i];
4153 if (gx == 0 && gy == 0)
4154 {
4155 continue;
4156 }
4157
4158 TVector3 vg4(gx, gy, gz);
4159 glayer = i;
4160 gr = vg4.Perp();
4161 geta = vg4.Eta();
4162 gphi = vg4.Phi();
4163 gpx = g4particle->get_px();
4164 gpy = g4particle->get_py();
4165 gpz = g4particle->get_pz();
4166 TVector3 vg4p(gpx, gpy, gpz);
4167
4168 gtpt = vg4p.Perp();
4169 gtphi = vg4p.Phi();
4170 gteta = vg4p.Eta();
4171
4172 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
4173
4174 if (vtx)
4175 {
4176 gvx = vtx->get_x();
4177 gvy = vtx->get_y();
4178 gvz = vtx->get_z();
4179 }
4180
4181 gembed = trutheval->get_embed(g4particle);
4182 gprimary = trutheval->is_primary(g4particle);
4183 gflav = g4particle->get_pid();
4184 if (i >= 1)
4185 {
4186 if (xval[i - 1] != 0 && yval[i - 1] != 0)
4187 {
4188 TVector3 vg4prev(xval[i - 1], yval[i - 1], zval[i - 1]);
4189 dphiprev = vg4.DeltaPhi(vg4prev);
4190 detaprev = geta - vg4prev.Eta();
4191 }
4192 }
4193
4194 float ntrk_f = ntrk;
4195 float _ievent_f = _ievent;
4196 float gseed_data[] = {_ievent_f,
4197 ntrk_f,
4198 gx,
4199 gy,
4200 gz,
4201 gr,
4202 geta,
4203 gphi,
4204 glayer,
4205 gpx,
4206 gpy,
4207 gpz,
4208 gtpt,
4209 gtphi,
4210 gteta,
4211 gvx,
4212 gvy,
4213 gvz,
4214 gembed,
4215 gprimary,
4216 gflav,
4217 dphiprev,
4218 detaprev,
4219 nhit_tpc_all,
4220 nhit_tpc_in,
4221 nhit_tpc_mid,
4222 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
4223
4224 _ntp_gseed->Fill(gseed_data);
4225 }
4226 delete [] xval;
4227 delete [] yval;
4228 delete [] zval;
4229 }
4230 }
4231
4232 if (Verbosity() > 1)
4233 {
4234 _timer->stop();
4235 std::cout << "g4hit time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
4236 }
4237 }
4238 return;
4239 }
4240
4241 TMatrixF SvtxEvaluator::calculateClusterError(TrkrCluster* c, float& clusphi)
4242 {
4243 TMatrixF localErr(3, 3);
4244 localErr[0][0] = 0.;
4245 localErr[0][1] = 0.;
4246 localErr[0][2] = 0.;
4247 localErr[1][0] = 0.;
4248 localErr[1][1] = c->getActsLocalError(0, 0);
4249 localErr[1][2] = c->getActsLocalError(0, 1);
4250 localErr[2][0] = 0.;
4251 localErr[2][1] = c->getActsLocalError(1, 0);
4252 localErr[2][2] = c->getActsLocalError(1, 1);
4253
4254 TMatrixF ROT(3, 3);
4255 ROT[0][0] = cos(clusphi);
4256 ROT[0][1] = -sin(clusphi);
4257 ROT[0][2] = 0.0;
4258 ROT[1][0] = sin(clusphi);
4259 ROT[1][1] = cos(clusphi);
4260 ROT[1][2] = 0.0;
4261 ROT[2][0] = 0.0;
4262 ROT[2][1] = 0.0;
4263 ROT[2][2] = 1.0;
4264 TMatrixF ROT_T(3, 3);
4265 ROT_T.Transpose(ROT);
4266
4267 TMatrixF err(3, 3);
4268 err = ROT * localErr * ROT_T;
4269 return err;
4270 }
