File indexing completed on 2025-12-16 09:21:48
0001 #include "SvtxEvaluator.h"
0002
0003 #include "SvtxClusterEval.h"
0004 #include "SvtxEvalStack.h"
0005 #include "SvtxHitEval.h"
0006 #include "SvtxTrackEval.h"
0007 #include "SvtxTruthEval.h"
0008 #include "SvtxVertexEval.h"
0009
0010 #include <trackbase/ActsGeometry.h>
0011 #include <trackbase/ClusterErrorPara.h>
0012 #include <trackbase/TpcDefs.h>
0013 #include <trackbase/TrkrCluster.h>
0014 #include <trackbase/TrkrClusterContainer.h>
0015 #include <trackbase/TrkrClusterHitAssoc.h>
0016 #include <trackbase/TrkrClusterIterationMapv1.h>
0017 #include <trackbase/TrkrDefs.h>
0018 #include <trackbase/TrkrHit.h>
0019 #include <trackbase/TrkrHitSet.h>
0020 #include <trackbase/TrkrHitSetContainer.h>
0021
0022 #include <trackbase_historic/ActsTransformations.h>
0023 #include <trackbase_historic/SvtxTrack.h>
0024 #include <trackbase_historic/SvtxTrackMap.h>
0025 #include <trackbase_historic/TrackAnalysisUtils.h>
0026 #include <trackbase_historic/TrackSeed.h>
0027
0028 #include <globalvertex/GlobalVertex.h>
0029 #include <globalvertex/GlobalVertexMap.h>
0030 #include <globalvertex/SvtxVertex.h>
0031 #include <globalvertex/SvtxVertexMap.h>
0032
0033 #include <g4main/PHG4Hit.h>
0034 #include <g4main/PHG4Particle.h>
0035 #include <g4main/PHG4TruthInfoContainer.h>
0036 #include <g4main/PHG4VtxPoint.h>
0037
0038 #include <g4detectors/PHG4TpcGeom.h>
0039 #include <g4detectors/PHG4TpcGeomContainer.h>
0040
0041 #include <fun4all/Fun4AllReturnCodes.h>
0042 #include <fun4all/SubsysReco.h>
0043
0044 #include <phool/PHTimer.h>
0045 #include <phool/getClass.h>
0046 #include <phool/phool.h>
0047 #include <phool/recoConsts.h>
0048
0049 #include <TFile.h>
0050 #include <TNtuple.h>
0051 #include <TVector3.h>
0052
0053 #include <algorithm>
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 = std::numeric_limits<float>::quiet_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);
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 = std::numeric_limits<float>::quiet_NaN();
0493 float vy = std::numeric_limits<float>::quiet_NaN();
0494 float vz = std::numeric_limits<float>::quiet_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 PHG4TpcGeomContainer* geom_container =
0583 findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0584 {
0585 if (!geom_container)
0586 {
0587 std::cout << PHWHERE << "ERROR: Can't find node TPCGEOMCONTAINER" << std::endl;
0588 }
0589 return;
0590 }
0591
0592 for (unsigned int ilayer = 0; ilayer < _nlayers_maps + _nlayers_intt + _nlayers_tpc; ++ilayer)
0593 {
0594 PHG4TpcGeom* 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);
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 PHG4TpcGeomContainer* geom_container =
0900 findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
0901 if (!geom_container)
0902 {
0903 std::cout << PHWHERE << "ERROR: Can't find node TPCGEOMCONTAINER" << std::endl;
0904 return;
0905 }
0906 PHG4TpcGeom* 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);
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 = std::numeric_limits<float>::quiet_NaN();
1118 float geta = std::numeric_limits<float>::quiet_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 && std::abs(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 = std::numeric_limits<float>::quiet_NaN();
1197 float gvy = std::numeric_limits<float>::quiet_NaN();
1198 float gvz = std::numeric_limits<float>::quiet_NaN();
1199 float gvt = std::numeric_limits<float>::quiet_NaN();
1200 float gembed = std::numeric_limits<float>::quiet_NaN();
1201 float gntracks = truthinfo->GetNumPrimaryVertexParticles();
1202 float gntracksmaps = std::numeric_limits<float>::quiet_NaN();
1203 float gnembed = std::numeric_limits<float>::quiet_NaN();
1204 float nfromtruth = std::numeric_limits<float>::quiet_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 && std::abs(gvt) < 2000. && std::abs(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 = std::numeric_limits<float>::quiet_NaN();
1261 float vy = std::numeric_limits<float>::quiet_NaN();
1262 float vz = std::numeric_limits<float>::quiet_NaN();
1263 float ntracks = std::numeric_limits<float>::quiet_NaN();
1264
1265 float gvx = std::numeric_limits<float>::quiet_NaN();
1266 float gvy = std::numeric_limits<float>::quiet_NaN();
1267 float gvz = std::numeric_limits<float>::quiet_NaN();
1268 float gvt = std::numeric_limits<float>::quiet_NaN();
1269 float gembed = iter->first;
1270 float gntracks = std::numeric_limits<float>::quiet_NaN();
1271 float gntracksmaps = std::numeric_limits<float>::quiet_NaN();
1272 float gnembed = std::numeric_limits<float>::quiet_NaN();
1273 float nfromtruth = std::numeric_limits<float>::quiet_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 && std::abs(gvt) < 2000 && std::abs(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 = std::numeric_limits<float>::quiet_NaN();
1372 float vy = std::numeric_limits<float>::quiet_NaN();
1373 float vz = std::numeric_limits<float>::quiet_NaN();
1374 float ntracks = std::numeric_limits<float>::quiet_NaN();
1375 float nfromtruth = std::numeric_limits<float>::quiet_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 = std::abs(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 = std::numeric_limits<float>::quiet_NaN();
1447 float gpx = std::numeric_limits<float>::quiet_NaN();
1448 float gpy = std::numeric_limits<float>::quiet_NaN();
1449 float gpz = std::numeric_limits<float>::quiet_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 = std::numeric_limits<float>::quiet_NaN();
1454 float gvy = std::numeric_limits<float>::quiet_NaN();
1455 float gvz = std::numeric_limits<float>::quiet_NaN();
1456
1457 float gembed = std::numeric_limits<float>::quiet_NaN();
1458 float gprimary = std::numeric_limits<float>::quiet_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 = std::numeric_limits<float>::quiet_NaN();
1517 float x = std::numeric_limits<float>::quiet_NaN();
1518 float y = std::numeric_limits<float>::quiet_NaN();
1519 float z = std::numeric_limits<float>::quiet_NaN();
1520 float eta = std::numeric_limits<float>::quiet_NaN();
1521 float phi = std::numeric_limits<float>::quiet_NaN();
1522 float e = std::numeric_limits<float>::quiet_NaN();
1523 float adc = std::numeric_limits<float>::quiet_NaN();
1524 float local_layer = std::numeric_limits<float>::quiet_NaN();
1525 float size = std::numeric_limits<float>::quiet_NaN();
1526 float efromtruth = std::numeric_limits<float>::quiet_NaN();
1527 float dphitru = std::numeric_limits<float>::quiet_NaN();
1528 float detatru = std::numeric_limits<float>::quiet_NaN();
1529 float dztru = std::numeric_limits<float>::quiet_NaN();
1530 float drtru = std::numeric_limits<float>::quiet_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 PHG4TpcGeomContainer* local_geom_container =
1651 findNode::getClass<PHG4TpcGeomContainer>(topNode, "TPCGEOMCONTAINER");
1652 if (!local_geom_container)
1653 {
1654 std::cout << PHWHERE << "ERROR: Can't find node TPCGEOMCONTAINER" << 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 = std::numeric_limits<float>::quiet_NaN();
1689 float zbin = std::numeric_limits<float>::quiet_NaN();
1690 float phi = std::numeric_limits<float>::quiet_NaN();
1691 float r = std::numeric_limits<float>::quiet_NaN();
1692 float x = std::numeric_limits<float>::quiet_NaN();
1693 float y = std::numeric_limits<float>::quiet_NaN();
1694 float z = std::numeric_limits<float>::quiet_NaN();
1695
1696 if (local_layer >= _nlayers_maps + _nlayers_intt && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
1697 {
1698 PHG4TpcGeom* 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 * std::cos(phi);
1705 y = r * std::sin(phi);
1706 }
1707
1708 float g4hitID = std::numeric_limits<float>::quiet_NaN();
1709 float gedep = std::numeric_limits<float>::quiet_NaN();
1710 float gx = std::numeric_limits<float>::quiet_NaN();
1711 float gy = std::numeric_limits<float>::quiet_NaN();
1712 float gz = std::numeric_limits<float>::quiet_NaN();
1713 float gt = std::numeric_limits<float>::quiet_NaN();
1714 float gtrackID = std::numeric_limits<float>::quiet_NaN();
1715 float gflavor = std::numeric_limits<float>::quiet_NaN();
1716 float gpx = std::numeric_limits<float>::quiet_NaN();
1717 float gpy = std::numeric_limits<float>::quiet_NaN();
1718 float gpz = std::numeric_limits<float>::quiet_NaN();
1719 float gvx = std::numeric_limits<float>::quiet_NaN();
1720 float gvy = std::numeric_limits<float>::quiet_NaN();
1721 float gvz = std::numeric_limits<float>::quiet_NaN();
1722 float gvt = std::numeric_limits<float>::quiet_NaN();
1723 float gfpx = std::numeric_limits<float>::quiet_NaN();
1724 float gfpy = std::numeric_limits<float>::quiet_NaN();
1725 float gfpz = std::numeric_limits<float>::quiet_NaN();
1726 float gfx = std::numeric_limits<float>::quiet_NaN();
1727 float gfy = std::numeric_limits<float>::quiet_NaN();
1728 float gfz = std::numeric_limits<float>::quiet_NaN();
1729 float gembed = std::numeric_limits<float>::quiet_NaN();
1730 float gprimary = std::numeric_limits<float>::quiet_NaN();
1731
1732 float efromtruth = std::numeric_limits<float>::quiet_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 phibin,
1805 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 = std::numeric_limits<float>::quiet_NaN();
1947 float pedge = std::numeric_limits<float>::quiet_NaN();
1948 float ovlp = std::numeric_limits<float>::quiet_NaN();
1949
1950 auto para_errors = ClusterErrorPara::get_clusterv5_modified_error(cluster, r, cluster_key);
1951 phisize = cluster->getPhiSize();
1952 zsize = cluster->getZSize();
1953
1954 ez = sqrt(para_errors.second);
1955 ephi = sqrt(para_errors.first);
1956 maxadc = cluster->getMaxAdc();
1957 pedge = cluster->getEdge();
1958 ovlp = cluster->getOverlap();
1959
1960 if (hitsetlayer == 7 || hitsetlayer == 22 || hitsetlayer == 23 || hitsetlayer == 38 || hitsetlayer == 39)
1961 {
1962 redge = 1;
1963 }
1964
1965 float e = cluster->getAdc();
1966 float adc = cluster->getAdc();
1967 float local_layer = (float) TrkrDefs::getLayer(cluster_key);
1968 float sector = TpcDefs::getSectorId(cluster_key);
1969 float side = TpcDefs::getSide(cluster_key);
1970
1971
1972 TrkrDefs::hitsetkey local_hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
1973 int hitsetlayer2 = TrkrDefs::getLayer(local_hitsetkey);
1974 if (hitsetlayer != local_layer)
1975 {
1976 std::cout << "WARNING hitset layer " << hitsetlayer << "| " << hitsetlayer2 << " layer " << local_layer << std::endl;
1977 }
1978
1979
1980
1981
1982 float sumadc = 0;
1983 TrkrHitSetContainer::Iterator hitset = hitsets->findOrAddHitSet(local_hitsetkey);
1984 std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
1985 hitrange = clusterhitmap->getHits(cluster_key);
1986 for (std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
1987 clushititer = hitrange.first;
1988 clushititer != hitrange.second; ++clushititer)
1989 {
1990 TrkrHit* hit = hitset->second->getHit(clushititer->second);
1991 if (!hit)
1992 {
1993 continue;
1994 }
1995
1996 ++size;
1997 sumadc += (hit->getAdc() - 70);
1998 maxadc = std::max<float>(hit->getAdc() - 70, maxadc);
1999 }
2000 e = sumadc;
2001
2002 float trackID = std::numeric_limits<float>::quiet_NaN();
2003 if (track != nullptr)
2004 {
2005 trackID = track->get_id();
2006 }
2007
2008 float g4hitID = std::numeric_limits<float>::quiet_NaN();
2009 float gx = std::numeric_limits<float>::quiet_NaN();
2010 float gy = std::numeric_limits<float>::quiet_NaN();
2011 float gz = std::numeric_limits<float>::quiet_NaN();
2012 float gr = std::numeric_limits<float>::quiet_NaN();
2013 float gphi = std::numeric_limits<float>::quiet_NaN();
2014
2015 float geta = std::numeric_limits<float>::quiet_NaN();
2016 float gt = std::numeric_limits<float>::quiet_NaN();
2017 float gtrackID = std::numeric_limits<float>::quiet_NaN();
2018 float gflavor = std::numeric_limits<float>::quiet_NaN();
2019 float gpx = std::numeric_limits<float>::quiet_NaN();
2020 float gpy = std::numeric_limits<float>::quiet_NaN();
2021 float gpz = std::numeric_limits<float>::quiet_NaN();
2022 float gvx = std::numeric_limits<float>::quiet_NaN();
2023 float gvy = std::numeric_limits<float>::quiet_NaN();
2024 float gvz = std::numeric_limits<float>::quiet_NaN();
2025 float gvt = std::numeric_limits<float>::quiet_NaN();
2026 float gfpx = std::numeric_limits<float>::quiet_NaN();
2027 float gfpy = std::numeric_limits<float>::quiet_NaN();
2028 float gfpz = std::numeric_limits<float>::quiet_NaN();
2029 float gfx = std::numeric_limits<float>::quiet_NaN();
2030 float gfy = std::numeric_limits<float>::quiet_NaN();
2031 float gfz = std::numeric_limits<float>::quiet_NaN();
2032 float gembed = std::numeric_limits<float>::quiet_NaN();
2033 float gprimary = std::numeric_limits<float>::quiet_NaN();
2034
2035 float efromtruth = std::numeric_limits<float>::quiet_NaN();
2036
2037 if (Verbosity() > 1)
2038 {
2039 std::cout << PHWHERE << " **** reco: layer " << local_layer << std::endl;
2040 std::cout << " reco cluster key " << cluster_key << " r " << r << " x " << x << " y " << y << " z " << z << " phi " << phi << " adc " << adc << std::endl;
2041 }
2042 float nparticles = std::numeric_limits<float>::quiet_NaN();
2043
2044
2045 const auto [truth_ckey, truth_cluster] = clustereval->max_truth_cluster_by_energy(cluster_key);
2046 if (truth_cluster)
2047 {
2048 if (Verbosity() > 1)
2049 {
2050 std::cout << "Found matching truth cluster with key " << truth_ckey << " for reco cluster key " << cluster_key << " in layer " << local_layer << std::endl;
2051 }
2052
2053 g4hitID = 0;
2054 gx = truth_cluster->getX();
2055 gy = truth_cluster->getY();
2056 gz = truth_cluster->getZ();
2057 efromtruth = truth_cluster->getError(0, 0);
2058
2059 TVector3 gpos(gx, gy, gz);
2060 gr = gpos.Perp();
2061 gphi = gpos.Phi();
2062 geta = gpos.Eta();
2063
2064 if (g4particle)
2065 {
2066 gtrackID = g4particle->get_track_id();
2067 gflavor = g4particle->get_pid();
2068 gpx = g4particle->get_px();
2069 gpy = g4particle->get_py();
2070 gpz = g4particle->get_pz();
2071
2072 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
2073 if (vtx)
2074 {
2075 gvx = vtx->get_x();
2076 gvy = vtx->get_y();
2077 gvz = vtx->get_z();
2078 gvt = vtx->get_t();
2079 }
2080
2081 PHG4Hit* outerhit = nullptr;
2082 if (_do_eval_light == false)
2083 {
2084 outerhit = trutheval->get_outermost_truth_hit(g4particle);
2085 }
2086 if (outerhit)
2087 {
2088 gfpx = outerhit->get_px(1);
2089 gfpy = outerhit->get_py(1);
2090 gfpz = outerhit->get_pz(1);
2091 gfx = outerhit->get_x(1);
2092 gfy = outerhit->get_y(1);
2093 gfz = outerhit->get_z(1);
2094 }
2095
2096 gembed = trutheval->get_embed(g4particle);
2097 gprimary = trutheval->is_primary(g4particle);
2098 }
2099
2100 if (Verbosity() > 1)
2101 {
2102 std::cout << " truth cluster key " << truth_ckey << " gr " << gr << " gx " << gx << " gy " << gy << " gz " << gz << " gphi " << gphi << " efromtruth " << efromtruth << std::endl;
2103 }
2104 }
2105
2106 if (g4particle)
2107 {
2108 }
2109 nparticles = clustereval->all_truth_particles(cluster_key).size();
2110
2111 float cluster_data[] = {(float) _ievent,
2112 (float) _iseed,
2113 hitID,
2114 x,
2115 y,
2116 z,
2117 r,
2118 phi,
2119 eta,
2120 theta,
2121 ex,
2122 ey,
2123 ez,
2124 ephi,
2125 pez,
2126 pephi,
2127 e,
2128 adc,
2129 maxadc,
2130 local_layer,
2131 sector,
2132 side,
2133 size,
2134 phisize,
2135 zsize,
2136 pedge,
2137 redge,
2138 ovlp,
2139 trackID,
2140 niter,
2141 g4hitID,
2142 gx,
2143 gy,
2144 gz,
2145 gr,
2146 gphi,
2147 geta,
2148 gt,
2149 gtrackID,
2150 gflavor,
2151 gpx,
2152 gpy,
2153 gpz,
2154 gvx,
2155 gvy,
2156 gvz,
2157 gvt,
2158 gfpx,
2159 gfpy,
2160 gfpz,
2161 gfx,
2162 gfy,
2163 gfz,
2164 gembed,
2165 gprimary,
2166 efromtruth,
2167 nparticles,
2168 nhit_tpc_all,
2169 nhit_tpc_in,
2170 nhit_tpc_mid,
2171 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
2172
2173 _ntp_cluster->Fill(cluster_data);
2174 }
2175 }
2176 }
2177 }
2178 else if (_ntp_cluster && _scan_for_embedded)
2179 {
2180 if (Verbosity() > 1)
2181 {
2182 std::cout << "Filling ntp_cluster (embedded only) " << std::endl;
2183 }
2184
2185
2186
2187
2188
2189
2190 SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname);
2191
2192 TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
2193 if (!clustermap)
2194 {
2195 clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
2196 }
2197 ClusterErrorPara ClusErrPara;
2198
2199 TrkrClusterHitAssoc* clusterhitmap = findNode::getClass<TrkrClusterHitAssoc>(topNode, "TRKR_CLUSTERHITASSOC");
2200 TrkrHitSetContainer* hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
2201 TrkrClusterIterationMapv1* _iteration_map = findNode::getClass<TrkrClusterIterationMapv1>(topNode, "CLUSTER_ITERATION_MAP");
2202
2203 if (trackmap != nullptr && clustermap != nullptr && clusterhitmap != nullptr && hitsets != nullptr)
2204 {
2205 for (auto& iter : *trackmap)
2206 {
2207 SvtxTrack* track = iter.second;
2208
2209 PHG4Particle* truth = trackeval->max_truth_particle_by_nclusters(track);
2210 if (truth)
2211 {
2212 if (trutheval->get_embed(truth) <= 0)
2213 {
2214 continue;
2215 }
2216 }
2217
2218 std::vector<TrkrDefs::cluskey> clusters;
2219 auto* siseed = track->get_silicon_seed();
2220 if (siseed)
2221 {
2222 for (auto local_iter = siseed->begin_cluster_keys();
2223 local_iter != siseed->end_cluster_keys();
2224 ++local_iter)
2225 {
2226 TrkrDefs::cluskey cluster_key = *local_iter;
2227 clusters.push_back(cluster_key);
2228 }
2229 }
2230 auto* tpcseed = track->get_tpc_seed();
2231 if (tpcseed)
2232 {
2233 for (auto local_iter = tpcseed->begin_cluster_keys();
2234 local_iter != tpcseed->end_cluster_keys();
2235 ++local_iter)
2236 {
2237 TrkrDefs::cluskey cluster_key = *local_iter;
2238 clusters.push_back(cluster_key);
2239 }
2240 }
2241
2242
2243 for (unsigned long cluster_key : clusters)
2244 {
2245 TrkrCluster* cluster = clustermap->findCluster(cluster_key);
2246 if (!cluster)
2247 {
2248 continue;
2249 }
2250
2251 PHG4Hit* g4hit = clustereval->max_truth_hit_by_energy(cluster_key);
2252 PHG4Particle* g4particle = trutheval->get_particle(g4hit);
2253
2254 float niter = 0;
2255 if (_iteration_map != nullptr)
2256 {
2257 niter = _iteration_map->getIteration(cluster_key);
2258 }
2259 float hitID = (float) TrkrDefs::getClusIndex(cluster_key);
2260 Acts::Vector3 glob = tgeometry->getGlobalPosition(cluster_key, cluster);
2261 float x = glob(0);
2262 float y = glob(1);
2263 float z = glob(2);
2264 TVector3 pos(x, y, z);
2265 float r = pos.Perp();
2266 float phi = pos.Phi();
2267 float eta = pos.Eta();
2268 float theta = pos.Theta();
2269 float ex = 0;
2270 float ey = 0;
2271 float ez = 0;
2272 float ephi = 0;
2273 float pez = 0;
2274 float pephi = 0;
2275 float size = 0;
2276 float phisize = 0;
2277 float zsize = 0;
2278 float maxadc = -999;
2279 float redge = std::numeric_limits<float>::quiet_NaN();
2280 float pedge = std::numeric_limits<float>::quiet_NaN();
2281 float ovlp = std::numeric_limits<float>::quiet_NaN();
2282
2283 auto para_errors = ClusterErrorPara::get_clusterv5_modified_error(cluster, r, cluster_key);
2284 phisize = cluster->getPhiSize();
2285 zsize = cluster->getZSize();
2286
2287 ez = sqrt(para_errors.second);
2288 ephi = sqrt(para_errors.first);
2289 maxadc = cluster->getMaxAdc();
2290 pedge = cluster->getEdge();
2291 ovlp = cluster->getOverlap();
2292
2293 float e = cluster->getAdc();
2294 float adc = cluster->getAdc();
2295 float local_layer = (float) TrkrDefs::getLayer(cluster_key);
2296 float sector = TpcDefs::getSectorId(cluster_key);
2297 float side = TpcDefs::getSide(cluster_key);
2298
2299 if (local_layer == 7 || local_layer == 22 || local_layer == 23 || local_layer == 38 || local_layer == 39)
2300 {
2301 redge = 1;
2302 }
2303
2304
2305 TrkrDefs::hitsetkey hitsetkey = TrkrDefs::getHitSetKeyFromClusKey(cluster_key);
2306 TrkrHitSetContainer::Iterator hitset = hitsets->findOrAddHitSet(hitsetkey);
2307 std::pair<std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator, std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator>
2308 hitrange = clusterhitmap->getHits(cluster_key);
2309 for (std::multimap<TrkrDefs::cluskey, TrkrDefs::hitkey>::const_iterator
2310 clushititer = hitrange.first;
2311 clushititer != hitrange.second; ++clushititer)
2312 {
2313 TrkrHit* hit = hitset->second->getHit(clushititer->second);
2314 ++size;
2315 maxadc = std::max<float>(hit->getAdc(), maxadc);
2316 }
2317
2318 float trackID = std::numeric_limits<float>::quiet_NaN();
2319 trackID = track->get_id();
2320
2321 float g4hitID = std::numeric_limits<float>::quiet_NaN();
2322 float gx = std::numeric_limits<float>::quiet_NaN();
2323 float gy = std::numeric_limits<float>::quiet_NaN();
2324 float gz = std::numeric_limits<float>::quiet_NaN();
2325 float gr = std::numeric_limits<float>::quiet_NaN();
2326 float gphi = std::numeric_limits<float>::quiet_NaN();
2327
2328 float geta = std::numeric_limits<float>::quiet_NaN();
2329 float gt = std::numeric_limits<float>::quiet_NaN();
2330 float gtrackID = std::numeric_limits<float>::quiet_NaN();
2331 float gflavor = std::numeric_limits<float>::quiet_NaN();
2332 float gpx = std::numeric_limits<float>::quiet_NaN();
2333 float gpy = std::numeric_limits<float>::quiet_NaN();
2334 float gpz = std::numeric_limits<float>::quiet_NaN();
2335 float gvx = std::numeric_limits<float>::quiet_NaN();
2336 float gvy = std::numeric_limits<float>::quiet_NaN();
2337 float gvz = std::numeric_limits<float>::quiet_NaN();
2338 float gvt = std::numeric_limits<float>::quiet_NaN();
2339 float gfpx = std::numeric_limits<float>::quiet_NaN();
2340 float gfpy = std::numeric_limits<float>::quiet_NaN();
2341 float gfpz = std::numeric_limits<float>::quiet_NaN();
2342 float gfx = std::numeric_limits<float>::quiet_NaN();
2343 float gfy = std::numeric_limits<float>::quiet_NaN();
2344 float gfz = std::numeric_limits<float>::quiet_NaN();
2345 float gembed = std::numeric_limits<float>::quiet_NaN();
2346 float gprimary = std::numeric_limits<float>::quiet_NaN();
2347
2348 float efromtruth = std::numeric_limits<float>::quiet_NaN();
2349
2350
2351 const auto [truth_ckey, truth_cluster] = clustereval->max_truth_cluster_by_energy(cluster_key);
2352 if (truth_cluster)
2353 {
2354 if (Verbosity() > 1)
2355 {
2356 std::cout << " Found matching truth cluster with key " << truth_ckey << " for reco cluster key " << cluster_key << " in layer " << local_layer << std::endl;
2357 }
2358
2359 g4hitID = 0;
2360 gx = truth_cluster->getX();
2361 gy = truth_cluster->getY();
2362 gz = truth_cluster->getZ();
2363 efromtruth = truth_cluster->getError(0, 0);
2364
2365 TVector3 gpos(gx, gy, gz);
2366 gr = gpos.Perp();
2367 gphi = gpos.Phi();
2368 geta = gpos.Eta();
2369
2370 if (g4particle)
2371 {
2372 gtrackID = g4particle->get_track_id();
2373 gflavor = g4particle->get_pid();
2374 gpx = g4particle->get_px();
2375 gpy = g4particle->get_py();
2376 gpz = g4particle->get_pz();
2377
2378 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
2379 if (vtx)
2380 {
2381 gvx = vtx->get_x();
2382 gvy = vtx->get_y();
2383 gvz = vtx->get_z();
2384 gvt = vtx->get_t();
2385 }
2386 PHG4Hit* outerhit = nullptr;
2387 if (_do_eval_light == false)
2388 {
2389 outerhit = trutheval->get_outermost_truth_hit(g4particle);
2390 }
2391 if (outerhit)
2392 {
2393 gfpx = outerhit->get_px(1);
2394 gfpy = outerhit->get_py(1);
2395 gfpz = outerhit->get_pz(1);
2396 gfx = outerhit->get_x(1);
2397 gfy = outerhit->get_y(1);
2398 gfz = outerhit->get_z(1);
2399 }
2400
2401 gembed = trutheval->get_embed(g4particle);
2402 gprimary = trutheval->is_primary(g4particle);
2403 }
2404 }
2405
2406 float nparticles = clustereval->all_truth_particles(cluster_key).size();
2407
2408 float cluster_data[] = {(float) _ievent,
2409 (float) _iseed,
2410 hitID,
2411 x,
2412 y,
2413 z,
2414 r,
2415 phi,
2416 eta,
2417 theta,
2418 ex,
2419 ey,
2420 ez,
2421 ephi,
2422 pez,
2423 pephi,
2424 e,
2425 adc,
2426 maxadc,
2427 local_layer,
2428 sector,
2429 side,
2430 size,
2431 phisize,
2432 zsize,
2433 pedge,
2434 redge,
2435 ovlp,
2436 trackID,
2437 niter,
2438 g4hitID,
2439 gx,
2440 gy,
2441 gz,
2442 gr,
2443 gphi,
2444 geta,
2445 gt,
2446 gtrackID,
2447 gflavor,
2448 gpx,
2449 gpy,
2450 gpz,
2451 gvx,
2452 gvy,
2453 gvz,
2454 gvt,
2455 gfpx,
2456 gfpy,
2457 gfpz,
2458 gfx,
2459 gfy,
2460 gfz,
2461 gembed,
2462 gprimary,
2463 efromtruth,
2464 nparticles,
2465 nhit_tpc_all,
2466 nhit_tpc_in,
2467 nhit_tpc_mid,
2468 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
2469
2470 _ntp_cluster->Fill(cluster_data);
2471 }
2472 }
2473 }
2474 }
2475 if (Verbosity() >= 1)
2476 {
2477 _timer->stop();
2478 std::cout << "cluster time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
2479 }
2480
2481
2482
2483
2484 if (Verbosity() > 1)
2485 {
2486 std::cout << "check for ntp_g4cluster" << std::endl;
2487 _timer->restart();
2488 }
2489
2490 if (_ntp_g4cluster)
2491 {
2492 if (Verbosity() >= 1)
2493 {
2494 std::cout << "Filling ntp_g4cluster " << std::endl;
2495 }
2496 ClusterErrorPara ClusErrPara;
2497 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
2498 PHG4TruthInfoContainer::ConstRange range = truthinfo->GetParticleRange();
2499 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
2500 iter != range.second;
2501 ++iter)
2502 {
2503 PHG4Particle* g4particle = iter->second;
2504
2505 if (_scan_for_embedded)
2506 {
2507 if (trutheval->get_embed(g4particle) <= 0)
2508 {
2509 continue;
2510 }
2511 }
2512
2513 float gtrackID = g4particle->get_track_id();
2514 float gflavor = g4particle->get_pid();
2515 float gembed = trutheval->get_embed(g4particle);
2516 float gprimary = trutheval->is_primary(g4particle);
2517
2518 if (Verbosity() > 1)
2519 {
2520 std::cout << PHWHERE << " PHG4Particle ID " << gtrackID << " gflavor " << gflavor << " gprimary " << gprimary << std::endl;
2521 }
2522
2523
2524 std::map<TrkrDefs::cluskey, std::shared_ptr<TrkrCluster> > truth_clusters = trutheval->all_truth_clusters(g4particle);
2525
2526
2527 for (const auto& [ckey, gclus] : truth_clusters)
2528 {
2529 unsigned int local_layer = TrkrDefs::getLayer(ckey);
2530
2531 float gx = gclus->getX();
2532 float gy = gclus->getY();
2533 float gz = gclus->getZ();
2534 float gt = std::numeric_limits<float>::quiet_NaN();
2535 float gedep = gclus->getError(0, 0);
2536 float gadc = (float) gclus->getAdc();
2537
2538 TVector3 gpos(gx, gy, gz);
2539 float gr = std::sqrt(gx * gx + gy * gy);
2540 float gphi = gpos.Phi();
2541 float geta = gpos.Eta();
2542
2543 if (Verbosity() > 1)
2544 {
2545 std::cout << PHWHERE << " **** truth: layer " << local_layer << std::endl;
2546 std::cout << " truth cluster key " << ckey << " gr " << gr << " gx " << gx << " gy " << gy << " gz " << gz
2547 << " gphi " << gphi << " gedep " << gedep << " gadc " << gadc << std::endl;
2548 }
2549
2550 float gphisize = gclus->getSize(1, 1);
2551 float gzsize = gclus->getSize(2, 2);
2552
2553
2554
2555 float x = std::numeric_limits<float>::quiet_NaN();
2556 float y = std::numeric_limits<float>::quiet_NaN();
2557 float z = std::numeric_limits<float>::quiet_NaN();
2558 float r = std::numeric_limits<float>::quiet_NaN();
2559 float phi = std::numeric_limits<float>::quiet_NaN();
2560 float eta = std::numeric_limits<float>::quiet_NaN();
2561 float ex = std::numeric_limits<float>::quiet_NaN();
2562 float ey = std::numeric_limits<float>::quiet_NaN();
2563 float ez = std::numeric_limits<float>::quiet_NaN();
2564 float ephi = std::numeric_limits<float>::quiet_NaN();
2565 float adc = std::numeric_limits<float>::quiet_NaN();
2566
2567 float nreco = 0;
2568 float phisize = 0;
2569 float zsize = 0;
2570 const auto [reco_ckey, reco_cluster] = clustereval->reco_cluster_from_truth_cluster(ckey, gclus);
2571 if (reco_cluster)
2572 {
2573 nreco = 1;
2574
2575 Acts::Vector3 glob;
2576 glob = tgeometry->getGlobalPosition(reco_ckey, reco_cluster);
2577 x = glob(0);
2578 y = glob(1);
2579 z = glob(2);
2580
2581 TVector3 pos(x, y, z);
2582 r = std::sqrt(x * x + y * y);
2583 phi = pos.Phi();
2584 eta = pos.Eta();
2585
2586 auto para_errors = ClusterErrorPara::get_clusterv5_modified_error(reco_cluster, r, ckey);
2587
2588 phisize = reco_cluster->getPhiSize();
2589 zsize = reco_cluster->getZSize();
2590 ez = sqrt(para_errors.second);
2591 ephi = sqrt(para_errors.first);
2592
2593 adc = reco_cluster->getAdc();
2594
2595 if (Verbosity() > 1)
2596 {
2597 std::cout << " reco cluster key " << reco_ckey << " r " << r << " x " << x << " y " << y << " z " << z << " phi " << phi << " adc " << adc << std::endl;
2598 }
2599 }
2600 if (nreco == 0 && Verbosity() > 1)
2601 {
2602 if (Verbosity() > 1)
2603 {
2604 std::cout << " ----------- Failed to find matching reco cluster " << std::endl;
2605 }
2606 }
2607
2608
2609
2610 float g4cluster_data[] = {(float) _ievent,
2611 (float) local_layer,
2612 gx,
2613 gy,
2614 gz,
2615 gt,
2616 gedep,
2617 gr,
2618 gphi,
2619 geta,
2620 gtrackID,
2621 gflavor,
2622 gembed,
2623 gprimary,
2624 gphisize,
2625 gzsize,
2626 gadc,
2627 nreco,
2628 x,
2629 y,
2630 z,
2631 r,
2632 phi,
2633 eta,
2634 ex,
2635 ey,
2636 ez,
2637 ephi,
2638 adc,
2639 phisize,
2640 zsize};
2641 _ntp_g4cluster->Fill(g4cluster_data);
2642 }
2643 }
2644 }
2645
2646
2647
2648
2649
2650
2651
2652
2653
2654 if (_ntp_gtrack)
2655 {
2656 if (Verbosity() > 1)
2657 {
2658 std::cout << "Filling ntp_gtrack " << std::endl;
2659 _timer->restart();
2660 }
2661 ClusterErrorPara ClusErrPara;
2662 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
2663 GlobalVertexMap* gvertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
2664 TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
2665 if (!clustermap)
2666 {
2667 clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
2668 }
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 (std::abs(gflavor) == 211 || std::abs(gflavor) == 321 || std::abs(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 auto g4clustermap = trutheval->all_truth_clusters(g4particle);
2716 std::set<TrkrDefs::cluskey> g4clusters;
2717 for(const auto& [key, cluster]: g4clustermap)
2718 {
2719 g4clusters.insert(key);
2720 }
2721
2722 float ng4hits = g4clusters.size();
2723 unsigned int ngmaps = 0;
2724 unsigned int ngmms = 0;
2725 unsigned int ngintt = 0;
2726 unsigned int ngintt1 = 0;
2727 unsigned int ngintt2 = 0;
2728 unsigned int ngintt3 = 0;
2729 unsigned int ngintt4 = 0;
2730 unsigned int ngintt5 = 0;
2731 unsigned int ngintt6 = 0;
2732 unsigned int ngintt7 = 0;
2733 unsigned int ngintt8 = 0;
2734 unsigned int ngtpc = 0;
2735 unsigned int nglmaps = 0;
2736 unsigned int nglintt = 0;
2737 unsigned int ngltpc = 0;
2738 unsigned int nglmms = 0;
2739
2740 std::vector<int> lmaps(_nlayers_maps + 1, 0);
2741 std::vector<int> lintt(_nlayers_intt + 1, 0);
2742 std::vector<int> ltpc(_nlayers_tpc + 1, 0);
2743 std::vector<int> lmms(_nlayers_mms + 1, 0);
2744
2745 for (const TrkrDefs::cluskey g4cluster : g4clusters)
2746 {
2747 unsigned int local_layer = TrkrDefs::getLayer(g4cluster);
2748
2749 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
2750 {
2751 lmaps[local_layer] = 1;
2752 ngmaps++;
2753 }
2754 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
2755 {
2756 lmms[local_layer - _nlayers_tpc - _nlayers_intt - _nlayers_maps] = 1;
2757 ngmms++;
2758 }
2759 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
2760 {
2761 lintt[local_layer - _nlayers_maps] = 1;
2762 ngintt++;
2763 }
2764
2765 if (_nlayers_intt > 0 && local_layer == _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
2766 {
2767 ngintt1++;
2768 }
2769
2770 if (_nlayers_intt > 1 && local_layer == _nlayers_maps + 1 && local_layer < _nlayers_maps + _nlayers_intt)
2771 {
2772 ngintt2++;
2773 }
2774
2775 if (_nlayers_intt > 2 && local_layer == _nlayers_maps + 2 && local_layer < _nlayers_maps + _nlayers_intt)
2776 {
2777 ngintt3++;
2778 }
2779
2780 if (_nlayers_intt > 3 && local_layer == _nlayers_maps + 3 && local_layer < _nlayers_maps + _nlayers_intt)
2781 {
2782 ngintt4++;
2783 }
2784
2785 if (_nlayers_intt > 4 && local_layer == _nlayers_maps + 4 && local_layer < _nlayers_maps + _nlayers_intt)
2786 {
2787 ngintt5++;
2788 }
2789
2790 if (_nlayers_intt > 5 && local_layer == _nlayers_maps + 5 && local_layer < _nlayers_maps + _nlayers_intt)
2791 {
2792 ngintt6++;
2793 }
2794
2795 if (_nlayers_intt > 6 && local_layer == _nlayers_maps + 6 && local_layer < _nlayers_maps + _nlayers_intt)
2796 {
2797 ngintt7++;
2798 }
2799
2800 if (_nlayers_intt > 7 && local_layer == _nlayers_maps + 7 && local_layer < _nlayers_maps + _nlayers_intt)
2801 {
2802 ngintt8++;
2803 }
2804 if (_nlayers_tpc > 0 && local_layer >= _nlayers_maps + _nlayers_intt && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
2805 {
2806 ltpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
2807 ngtpc++;
2808 }
2809 }
2810 if (_nlayers_maps > 0)
2811 {
2812 for (unsigned int i = 0; i < _nlayers_maps; i++)
2813 {
2814 nglmaps += lmaps[i];
2815 }
2816 }
2817 if (_nlayers_intt > 0)
2818 {
2819 for (unsigned int i = 0; i < _nlayers_intt; i++)
2820 {
2821 nglintt += lintt[i];
2822 }
2823 }
2824 if (_nlayers_tpc > 0)
2825 {
2826 for (unsigned int i = 0; i < _nlayers_tpc; i++)
2827 {
2828 ngltpc += ltpc[i];
2829 }
2830 }
2831 if (_nlayers_mms > 0)
2832 {
2833 for (unsigned int i = 0; i < _nlayers_mms; i++)
2834 {
2835 nglmms += lmms[i];
2836 }
2837 }
2838
2839 float gpx = g4particle->get_px();
2840 float gpy = g4particle->get_py();
2841 float gpz = g4particle->get_pz();
2842 float gpt = std::numeric_limits<float>::quiet_NaN();
2843 float geta = std::numeric_limits<float>::quiet_NaN();
2844 float gphi = std::numeric_limits<float>::quiet_NaN();
2845 if (gpx != 0 && gpy != 0)
2846 {
2847 TVector3 gv(gpx, gpy, gpz);
2848 gpt = gv.Pt();
2849 geta = gv.Eta();
2850 gphi = gv.Phi();
2851 }
2852 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
2853 float gvx = vtx->get_x();
2854 float gvy = vtx->get_y();
2855 float gvz = vtx->get_z();
2856 float gvt = vtx->get_t();
2857
2858 float gfpx = 0.;
2859 float gfpy = 0.;
2860 float gfpz = 0.;
2861 float gfx = 0.;
2862 float gfy = 0.;
2863 float gfz = 0.;
2864
2865 PHG4Hit* outerhit = nullptr;
2866 if (_do_eval_light == false)
2867 {
2868 outerhit = trutheval->get_outermost_truth_hit(g4particle);
2869 }
2870
2871 if (outerhit)
2872 {
2873 gfpx = outerhit->get_px(1);
2874 gfpy = outerhit->get_py(1);
2875 gfpz = outerhit->get_pz(1);
2876 gfx = outerhit->get_x(1);
2877 gfy = outerhit->get_y(1);
2878 gfz = outerhit->get_z(1);
2879 }
2880
2881 float gembed = trutheval->get_embed(g4particle);
2882 float gprimary = trutheval->is_primary(g4particle);
2883
2884 float trackID = std::numeric_limits<float>::quiet_NaN();
2885 float charge = std::numeric_limits<float>::quiet_NaN();
2886 float quality = std::numeric_limits<float>::quiet_NaN();
2887 float chisq = std::numeric_limits<float>::quiet_NaN();
2888 float ndf = std::numeric_limits<float>::quiet_NaN();
2889 float local_nhits = 0;
2890 float nmaps = 0;
2891 float nintt = 0;
2892 float ntpc = 0;
2893 float nmms = 0;
2894 float ntpc1 = 0;
2895 float ntpc11 = 0;
2896 float ntpc2 = 0;
2897 float ntpc3 = 0;
2898 float nlintt = 0;
2899 float nlmaps = 0;
2900 float nltpc = 0;
2901 float nlmms = 0;
2902 unsigned int layers = 0x0;
2903 int vertexID = -1;
2904 float vx = std::numeric_limits<float>::quiet_NaN();
2905 float vy = std::numeric_limits<float>::quiet_NaN();
2906 float vz = std::numeric_limits<float>::quiet_NaN();
2907 float dca2d = std::numeric_limits<float>::quiet_NaN();
2908 float dca2dsigma = std::numeric_limits<float>::quiet_NaN();
2909 float dca3dxy = std::numeric_limits<float>::quiet_NaN();
2910 float dca3dxysigma = std::numeric_limits<float>::quiet_NaN();
2911 float dca3dz = std::numeric_limits<float>::quiet_NaN();
2912 float dca3dzsigma = std::numeric_limits<float>::quiet_NaN();
2913 float px = std::numeric_limits<float>::quiet_NaN();
2914 float py = std::numeric_limits<float>::quiet_NaN();
2915 float pz = std::numeric_limits<float>::quiet_NaN();
2916 float pt = std::numeric_limits<float>::quiet_NaN();
2917 float eta = std::numeric_limits<float>::quiet_NaN();
2918 float phi = std::numeric_limits<float>::quiet_NaN();
2919 float deltapt = std::numeric_limits<float>::quiet_NaN();
2920 float deltaeta = std::numeric_limits<float>::quiet_NaN();
2921 float deltaphi = std::numeric_limits<float>::quiet_NaN();
2922 float pcax = std::numeric_limits<float>::quiet_NaN();
2923 float pcay = std::numeric_limits<float>::quiet_NaN();
2924 float pcaz = std::numeric_limits<float>::quiet_NaN();
2925
2926 float nfromtruth = std::numeric_limits<float>::quiet_NaN();
2927 float nwrong = std::numeric_limits<float>::quiet_NaN();
2928 float ntrumaps = std::numeric_limits<float>::quiet_NaN();
2929 float nwrongmaps = std::numeric_limits<float>::quiet_NaN();
2930 float ntruintt = std::numeric_limits<float>::quiet_NaN();
2931 float nwrongintt = std::numeric_limits<float>::quiet_NaN();
2932 float ntrumms = std::numeric_limits<float>::quiet_NaN();
2933 float nwrongmms = std::numeric_limits<float>::quiet_NaN();
2934 float ntrutpc = std::numeric_limits<float>::quiet_NaN();
2935 float nwrongtpc = std::numeric_limits<float>::quiet_NaN();
2936 float ntrutpc1 = std::numeric_limits<float>::quiet_NaN();
2937 float nwrongtpc1 = std::numeric_limits<float>::quiet_NaN();
2938 float ntrutpc11 = std::numeric_limits<float>::quiet_NaN();
2939 float nwrongtpc11 = std::numeric_limits<float>::quiet_NaN();
2940 float ntrutpc2 = std::numeric_limits<float>::quiet_NaN();
2941 float nwrongtpc2 = std::numeric_limits<float>::quiet_NaN();
2942 float ntrutpc3 = std::numeric_limits<float>::quiet_NaN();
2943 float nwrongtpc3 = std::numeric_limits<float>::quiet_NaN();
2944 float layersfromtruth = std::numeric_limits<float>::quiet_NaN();
2945 float npedge = 0;
2946 float nredge = 0;
2947 float nbig = 0;
2948 float novlp = 0;
2949 float merr = 0;
2950 float msize = 0;
2951 float siqr = std::numeric_limits<float>::quiet_NaN();
2952 float siphi = std::numeric_limits<float>::quiet_NaN();
2953 float sithe = std::numeric_limits<float>::quiet_NaN();
2954 float six0 = std::numeric_limits<float>::quiet_NaN();
2955 float siy0 = std::numeric_limits<float>::quiet_NaN();
2956 float tpqr = std::numeric_limits<float>::quiet_NaN();
2957 float tpphi = std::numeric_limits<float>::quiet_NaN();
2958 float tpthe = std::numeric_limits<float>::quiet_NaN();
2959 float tpx0 = std::numeric_limits<float>::quiet_NaN();
2960 float tpy0 = std::numeric_limits<float>::quiet_NaN();
2961
2962 if (_do_track_match)
2963 {
2964 SvtxTrack* track = trackeval->best_track_from(g4particle);
2965
2966 if (track)
2967 {
2968 trackID = track->get_id();
2969 charge = track->get_charge();
2970 quality = track->get_quality();
2971 chisq = track->get_chisq();
2972 ndf = track->get_ndf();
2973 TrackSeed* silseed = track->get_silicon_seed();
2974 TrackSeed* tpcseed = track->get_tpc_seed();
2975 if (tpcseed)
2976 {
2977 local_nhits += tpcseed->size_cluster_keys();
2978 }
2979 if (silseed)
2980 {
2981 local_nhits += silseed->size_cluster_keys();
2982 }
2983
2984 std::vector<int> maps(_nlayers_maps + 1, 0);
2985 std::vector<int> intt(_nlayers_intt + 1, 0);
2986 std::vector<int> tpc(_nlayers_tpc + 1, 0);
2987 std::vector<int> mms(_nlayers_mms + 1, 0);
2988
2989 if (tpcseed)
2990 {
2991 tpqr = tpcseed->get_qOverR();
2992 tpphi = tpcseed->get_phi();
2993 tpthe = tpcseed->get_theta();
2994 tpx0 = tpcseed->get_X0();
2995 tpy0 = tpcseed->get_Y0();
2996 for (TrackSeed::ConstClusterKeyIter local_iter = tpcseed->begin_cluster_keys();
2997 local_iter != tpcseed->end_cluster_keys();
2998 ++local_iter)
2999 {
3000 TrkrDefs::cluskey cluster_key = *local_iter;
3001 TrkrCluster* cluster = clustermap->findCluster(cluster_key);
3002 unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
3003 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
3004 {
3005 maps[local_layer] = 1;
3006 nmaps++;
3007 }
3008 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
3009 {
3010 intt[local_layer - _nlayers_maps] = 1;
3011 nintt++;
3012 }
3013 if (_nlayers_tpc > 0 &&
3014 local_layer >= (_nlayers_maps + _nlayers_intt) &&
3015 local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
3016 {
3017 tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
3018 ntpc++;
3019
3020 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 8)
3021 {
3022 ntpc11++;
3023 }
3024
3025 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 16)
3026 {
3027
3028 ntpc1++;
3029 }
3030 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 32)
3031 {
3032
3033 ntpc2++;
3034 }
3035 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 48)
3036 {
3037
3038 ntpc3++;
3039 }
3040 }
3041
3042 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
3043 {
3044 mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3045 nmms++;
3046 }
3047
3048 Acts::Vector3 glob;
3049 glob = tgeometry->getGlobalPosition(cluster_key, cluster);
3050 float x = glob(0);
3051 float y = glob(1);
3052 float z = glob(2);
3053
3054 TVector3 pos(x, y, z);
3055 float r = std::sqrt(x * x + y * y);
3056 phi = pos.Phi();
3057 eta = pos.Eta();
3058 float gphisize = 0;
3059
3060 float gphierr = 0;
3061
3062 float govlp = 0;
3063 float gedge = 0;
3064 if (cluster != nullptr)
3065 {
3066 gphisize = cluster->getPhiSize();
3067
3068 auto para_errors = ClusterErrorPara::get_clusterv5_modified_error(cluster, r, cluster_key);
3069
3070 gphierr = sqrt(para_errors.first);
3071 govlp = cluster->getOverlap();
3072 gedge = cluster->getEdge();
3073
3074 if (gedge > 0)
3075 {
3076 npedge++;
3077 }
3078 if (gphisize >= 4)
3079 {
3080 nbig++;
3081 }
3082 if (govlp >= 2)
3083 {
3084 novlp++;
3085 }
3086 merr += gphierr;
3087 msize += gphisize;
3088 }
3089 if (local_layer == 7 || local_layer == 22 || local_layer == 23 || local_layer == 38 || local_layer == 39)
3090 {
3091 nredge++;
3092 }
3093 }
3094 }
3095
3096 if (silseed)
3097 {
3098 siqr = silseed->get_qOverR();
3099 siphi = silseed->get_phi();
3100 sithe = silseed->get_theta();
3101 six0 = silseed->get_X0();
3102 siy0 = silseed->get_Y0();
3103 for (TrackSeed::ConstClusterKeyIter local_iter = silseed->begin_cluster_keys();
3104 local_iter != silseed->end_cluster_keys();
3105 ++local_iter)
3106 {
3107 TrkrDefs::cluskey cluster_key = *local_iter;
3108
3109 unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
3110 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
3111 {
3112 maps[local_layer] = 1;
3113 nmaps++;
3114 }
3115 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
3116 {
3117 intt[local_layer - _nlayers_maps] = 1;
3118 nintt++;
3119 }
3120 if (_nlayers_tpc > 0 &&
3121 local_layer >= (_nlayers_maps + _nlayers_intt) &&
3122 local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
3123 {
3124 tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
3125 ntpc++;
3126
3127 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 8)
3128 {
3129 ntpc11++;
3130 }
3131
3132 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 16)
3133 {
3134
3135 ntpc1++;
3136 }
3137 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 32)
3138 {
3139
3140 ntpc2++;
3141 }
3142 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 48)
3143 {
3144
3145 ntpc3++;
3146 }
3147 }
3148
3149 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
3150 {
3151 mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3152 nmms++;
3153 }
3154 }
3155 }
3156
3157 if (_nlayers_maps > 0)
3158 {
3159 for (unsigned int i = 0; i < _nlayers_maps; i++)
3160 {
3161 nlmaps += maps[i];
3162 }
3163 }
3164 if (_nlayers_intt > 0)
3165 {
3166 for (unsigned int i = 0; i < _nlayers_intt; i++)
3167 {
3168 nlintt += intt[i];
3169 }
3170 }
3171 if (_nlayers_tpc > 0)
3172 {
3173 for (unsigned int i = 0; i < _nlayers_tpc; i++)
3174 {
3175 nltpc += tpc[i];
3176 }
3177 }
3178 if (_nlayers_mms > 0)
3179 {
3180 for (unsigned int i = 0; i < _nlayers_mms; i++)
3181 {
3182 nlmms += mms[i];
3183 }
3184 }
3185
3186 layers = nlmaps + nlintt + nltpc + nlmms;
3187
3188
3189
3190
3191
3192
3193
3194
3195
3196
3197
3198 vertexID = track->get_vertex_id();
3199
3200 GlobalVertex* vertex = gvertexmap->get(vertexID);
3201 if (vertex)
3202 {
3203 vx = vertex->get_x();
3204 vy = vertex->get_y();
3205 vz = vertex->get_z();
3206 Acts::Vector3 vert(vx, vy, vz);
3207 auto dcapair = TrackAnalysisUtils::get_dca(track, vert);
3208 dca3dxy = dcapair.first.first;
3209 dca3dxysigma = dcapair.first.second;
3210 dca3dz = dcapair.second.first;
3211 dca3dzsigma = dcapair.second.second;
3212 }
3213
3214 px = track->get_px();
3215 py = track->get_py();
3216 pz = track->get_pz();
3217 TVector3 v(px, py, pz);
3218 pt = v.Pt();
3219 eta = v.Eta();
3220 phi = v.Phi();
3221 float CVxx = track->get_error(3, 3);
3222 float CVxy = track->get_error(3, 4);
3223 float CVxz = track->get_error(3, 5);
3224 float CVyy = track->get_error(4, 4);
3225 float CVyz = track->get_error(4, 5);
3226 float CVzz = track->get_error(5, 5);
3227 deltapt = std::sqrt((CVxx * px * px + 2 * CVxy * px * py + CVyy * py * py) / (px * px + py * py));
3228 deltaeta = std::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)));
3229 deltaphi = std::sqrt((CVyy * px * px - 2 * CVxy * px * py + CVxx * py * py) / ((px * px + py * py) * (px * px + py * py)));
3230 pcax = track->get_x();
3231 pcay = track->get_y();
3232 pcaz = track->get_z();
3233
3234 nfromtruth = trackeval->get_nclusters_contribution(track, g4particle);
3235 nwrong = trackeval->get_nwrongclusters_contribution(track, g4particle);
3236
3237 if (_nlayers_maps == 0)
3238 {
3239 ntrumaps = 0;
3240 }
3241 else
3242 {
3243 auto pair = trackeval->get_layer_range_contribution(track, g4particle, 0, _nlayers_maps);
3244 ntrumaps = pair.first;
3245 nwrongmaps = pair.second;
3246 }
3247 if (_nlayers_intt == 0)
3248 {
3249 ntruintt = 0;
3250 }
3251 else
3252 {
3253 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps, _nlayers_maps + _nlayers_intt);
3254 ntruintt = pair.first;
3255 nwrongintt = pair.second;
3256 }
3257 if (_nlayers_mms == 0)
3258 {
3259 ntrumms = 0;
3260 }
3261 else
3262 {
3263 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + _nlayers_tpc, _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms);
3264 ntrumms = pair.first;
3265 nwrongmms = pair.second;
3266 }
3267 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
3268 ntrutpc = pair.first;
3269 nwrongtpc = pair.second;
3270 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 16);
3271 ntrutpc1 = pair.first;
3272 nwrongtpc1 = pair.second;
3273 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 8);
3274 ntrutpc11 = pair.first;
3275 nwrongtpc11 = pair.second;
3276 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + 16, _nlayers_maps + _nlayers_intt + 32);
3277 ntrutpc2 = pair.first;
3278 nwrongtpc2 = pair.second;
3279 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + 32, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
3280 ntrutpc3 = pair.first;
3281 nwrongtpc3 = pair.first;
3282 layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track, g4particle);
3283 }
3284 }
3285
3286 float gtrack_data[] = {(float) _ievent, m_fSeed,
3287 gntracks,
3288 gnchghad,
3289 gtrackID,
3290 gflavor,
3291 ng4hits,
3292 (float) ngmaps,
3293 (float) ngintt,
3294 (float) ngmms,
3295 (float) ngintt1,
3296 (float) ngintt2,
3297 (float) ngintt3,
3298 (float) ngintt4,
3299 (float) ngintt5,
3300 (float) ngintt6,
3301 (float) ngintt7,
3302 (float) ngintt8,
3303 (float) ngtpc,
3304 (float) nglmaps,
3305 (float) nglintt,
3306 (float) ngltpc,
3307 (float) nglmms,
3308 gpx,
3309 gpy,
3310 gpz,
3311 gpt,
3312 geta,
3313 gphi,
3314 gvx,
3315 gvy,
3316 gvz,
3317 gvt,
3318 gfpx,
3319 gfpy,
3320 gfpz,
3321 gfx,
3322 gfy,
3323 gfz,
3324 gembed,
3325 gprimary,
3326 trackID,
3327 px,
3328 py,
3329 pz,
3330 pt,
3331 eta,
3332 phi,
3333 deltapt,
3334 deltaeta,
3335 deltaphi,
3336 siqr,
3337 siphi,
3338 sithe,
3339 six0,
3340 siy0,
3341 tpqr,
3342 tpphi,
3343 tpthe,
3344 tpx0,
3345 tpy0,
3346 charge,
3347 quality,
3348 chisq,
3349 ndf,
3350 local_nhits,
3351 (float) layers,
3352 nmaps,
3353 nintt,
3354 ntpc,
3355 nmms,
3356 ntpc1,
3357 ntpc11,
3358 ntpc2,
3359 ntpc3,
3360 nlmaps,
3361 nlintt,
3362 nltpc,
3363 nlmms,
3364 (float) vertexID,
3365 vx,
3366 vy,
3367 vz,
3368 dca2d,
3369 dca2dsigma,
3370 dca3dxy,
3371 dca3dxysigma,
3372 dca3dz,
3373 dca3dzsigma,
3374 pcax,
3375 pcay,
3376 pcaz,
3377 nfromtruth,
3378 nwrong,
3379 ntrumaps,
3380 nwrongmaps,
3381 ntruintt,
3382 nwrongintt,
3383 ntrutpc,
3384 nwrongtpc,
3385 ntrumms,
3386 nwrongmms,
3387 ntrutpc1,
3388 nwrongtpc1,
3389 ntrutpc11,
3390 nwrongtpc11,
3391 ntrutpc2,
3392 nwrongtpc2,
3393 ntrutpc3,
3394 nwrongtpc3,
3395 layersfromtruth,
3396 npedge,
3397 nredge,
3398 nbig,
3399 novlp,
3400 merr,
3401 msize,
3402 nhit_tpc_all,
3403 nhit_tpc_in,
3404 nhit_tpc_mid,
3405 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
3406
3407
3408
3409
3410
3411
3412
3413
3414
3415 _ntp_gtrack->Fill(gtrack_data);
3416 }
3417 }
3418 if (Verbosity() > 1)
3419 {
3420 _timer->stop();
3421 std::cout << "gtrack time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
3422 }
3423 }
3424
3425
3426
3427
3428
3429 if (_ntp_track)
3430 {
3431 if (Verbosity() > 1)
3432 {
3433 std::cout << "Filling ntp_track " << std::endl;
3434 _timer->restart();
3435 }
3436
3437
3438 SvtxTrackMap* trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname);
3439
3440 GlobalVertexMap* gvertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
3441 TrkrClusterContainer* clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
3442 if (!clustermap)
3443 {
3444 clustermap = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
3445 }
3446 ClusterErrorPara ClusErrPara;
3447 if (trackmap)
3448 {
3449 for (auto& iter : *trackmap)
3450 {
3451 SvtxTrack* track = iter.second;
3452 float trackID = track->get_id();
3453
3454 TrackSeed* tpcseed = track->get_tpc_seed();
3455 TrackSeed* silseed = track->get_silicon_seed();
3456 short int crossing_int = track->get_crossing();
3457 float crossing;
3458 if (crossing_int == SHRT_MAX)
3459 {
3460 crossing = std::numeric_limits<float>::quiet_NaN();
3461 }
3462 else
3463 {
3464 crossing = (float) crossing_int;
3465 }
3466 float charge = track->get_charge();
3467 float quality = track->get_quality();
3468 float chisq = track->get_chisq();
3469 float ndf = track->get_ndf();
3470 float local_nhits = 0;
3471 if (tpcseed)
3472 {
3473 local_nhits += tpcseed->size_cluster_keys();
3474 }
3475 if (silseed)
3476 {
3477 local_nhits += silseed->size_cluster_keys();
3478 }
3479 unsigned int layers = 0x0;
3480
3481 std::vector<int> maps(_nlayers_maps + 1, 0);
3482 std::vector<int> intt(_nlayers_intt + 1, 0);
3483 std::vector<int> tpc(_nlayers_tpc + 1, 0);
3484 std::vector<int> mms(_nlayers_mms + 1, 0);
3485
3486 float nmaps = 0;
3487 float nintt = 0;
3488 float nmms = 0;
3489 float ntpc = 0;
3490 float ntpc1 = 0;
3491 float ntpc11 = 0;
3492 float ntpc2 = 0;
3493 float ntpc3 = 0;
3494 float nlmaps = 0;
3495 float nlintt = 0;
3496 float nltpc = 0;
3497 float nlmms = 0;
3498 float npedge = 0;
3499 float nredge = 0;
3500 float nbig = 0;
3501 float novlp = 0;
3502 float merr = 0;
3503 float msize = 0;
3504 float siqr = std::numeric_limits<float>::quiet_NaN();
3505 float siphi = std::numeric_limits<float>::quiet_NaN();
3506 float sithe = std::numeric_limits<float>::quiet_NaN();
3507 float six0 = std::numeric_limits<float>::quiet_NaN();
3508 float siy0 = std::numeric_limits<float>::quiet_NaN();
3509 float tpqr = std::numeric_limits<float>::quiet_NaN();
3510 float tpphi = std::numeric_limits<float>::quiet_NaN();
3511 float tpthe = std::numeric_limits<float>::quiet_NaN();
3512 float tpx0 = std::numeric_limits<float>::quiet_NaN();
3513 float tpy0 = std::numeric_limits<float>::quiet_NaN();
3514
3515 if (tpcseed)
3516 {
3517 tpqr = tpcseed->get_qOverR();
3518 tpphi = tpcseed->get_phi();
3519 tpthe = tpcseed->get_theta();
3520 tpx0 = tpcseed->get_X0();
3521 tpy0 = tpcseed->get_Y0();
3522 for (SvtxTrack::ConstClusterKeyIter local_iter = tpcseed->begin_cluster_keys();
3523 local_iter != tpcseed->end_cluster_keys();
3524 ++local_iter)
3525 {
3526 TrkrDefs::cluskey cluster_key = *local_iter;
3527 TrkrCluster* cluster = clustermap->findCluster(cluster_key);
3528 unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
3529
3530 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
3531 {
3532 maps[local_layer] = 1;
3533 nmaps++;
3534 }
3535 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
3536 {
3537 intt[local_layer - _nlayers_maps] = 1;
3538 nintt++;
3539 }
3540 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
3541 {
3542 mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3543 nmms++;
3544 }
3545 if (_nlayers_tpc > 0 && local_layer >= (_nlayers_maps + _nlayers_intt) && local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
3546 {
3547 tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
3548 ntpc++;
3549 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 8)
3550 {
3551 ntpc11++;
3552 }
3553
3554 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 16)
3555 {
3556 ntpc1++;
3557 }
3558 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 32)
3559 {
3560 ntpc2++;
3561 }
3562 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 48)
3563 {
3564 ntpc3++;
3565 }
3566 }
3567
3568 Acts::Vector3 glob;
3569 glob = tgeometry->getGlobalPosition(cluster_key, cluster);
3570 float x = glob(0);
3571 float y = glob(1);
3572 float z = glob(2);
3573
3574 TVector3 pos(x, y, z);
3575 float r = std::sqrt(x * x + y * y);
3576
3577
3578 float rphisize = 0;
3579
3580 float rphierr = 0;
3581
3582 float rovlp = 0;
3583 float pedge = 0;
3584 if (cluster != nullptr)
3585 {
3586 rphisize = cluster->getPhiSize();
3587
3588 auto para_errors = ClusterErrorPara::get_clusterv5_modified_error(cluster, r, cluster_key);
3589
3590 rphierr = sqrt(para_errors.first);
3591 rovlp = cluster->getOverlap();
3592 pedge = cluster->getEdge();
3593
3594 if (pedge > 0)
3595 {
3596 npedge++;
3597 }
3598 if (rphisize >= 4)
3599 {
3600 nbig++;
3601 }
3602 if (rovlp >= 2)
3603 {
3604 novlp++;
3605 }
3606 merr += rphierr;
3607 msize += rphisize;
3608 }
3609 if (local_layer == 7 || local_layer == 22 || local_layer == 23 || local_layer == 38 || local_layer == 39)
3610 {
3611 nredge++;
3612 }
3613 if (Verbosity() > 2)
3614 {
3615 std::cout << " lay: " << local_layer
3616 << " pedge " << pedge
3617 << " | " << npedge
3618 << " nredge " << nredge
3619 << " rphisize " << rphisize
3620 << " | " << nbig
3621 << " rovlp " << rovlp
3622 << " | " << novlp
3623 << std::endl;
3624 }
3625 }
3626 }
3627
3628 if (silseed)
3629 {
3630 siqr = silseed->get_qOverR();
3631 siphi = silseed->get_phi();
3632 sithe = silseed->get_theta();
3633 six0 = silseed->get_X0();
3634 siy0 = silseed->get_Y0();
3635 for (SvtxTrack::ConstClusterKeyIter local_iter = silseed->begin_cluster_keys();
3636 local_iter != silseed->end_cluster_keys();
3637 ++local_iter)
3638 {
3639 TrkrDefs::cluskey cluster_key = *local_iter;
3640
3641 unsigned int local_layer = TrkrDefs::getLayer(cluster_key);
3642
3643 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
3644 {
3645 maps[local_layer] = 1;
3646 nmaps++;
3647 }
3648 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
3649 {
3650 intt[local_layer - _nlayers_maps] = 1;
3651 nintt++;
3652 }
3653 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
3654 {
3655 mms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3656 nmms++;
3657 }
3658 if (_nlayers_tpc > 0 && local_layer >= (_nlayers_maps + _nlayers_intt) && local_layer < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
3659 {
3660 tpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
3661 ntpc++;
3662 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 8)
3663 {
3664 ntpc11++;
3665 }
3666
3667 if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 16)
3668 {
3669 ntpc1++;
3670 }
3671 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 32)
3672 {
3673 ntpc2++;
3674 }
3675 else if ((local_layer - (_nlayers_maps + _nlayers_intt)) < 48)
3676 {
3677 ntpc3++;
3678 }
3679 }
3680 }
3681 }
3682
3683 if (_nlayers_maps > 0)
3684 {
3685 for (unsigned int i = 0; i < _nlayers_maps; i++)
3686 {
3687 nlmaps += maps[i];
3688 }
3689 }
3690 if (_nlayers_intt > 0)
3691 {
3692 for (unsigned int i = 0; i < _nlayers_intt; i++)
3693 {
3694 nlintt += intt[i];
3695 }
3696 }
3697 if (_nlayers_tpc > 0)
3698 {
3699 for (unsigned int i = 0; i < _nlayers_tpc; i++)
3700 {
3701 nltpc += tpc[i];
3702 }
3703 }
3704 if (_nlayers_mms > 0)
3705 {
3706 for (unsigned int i = 0; i < _nlayers_mms; i++)
3707 {
3708 nlmms += mms[i];
3709 }
3710 }
3711 layers = nlmaps + nlintt + nltpc + nlmms;
3712
3713 float dca3dxy = std::numeric_limits<float>::quiet_NaN();
3714 float dca3dz = std::numeric_limits<float>::quiet_NaN();
3715 float dca3dxysigma = std::numeric_limits<float>::quiet_NaN();
3716 float dca3dzsigma = std::numeric_limits<float>::quiet_NaN();
3717 float dca2d = std::numeric_limits<float>::quiet_NaN();
3718 float dca2dsigma = std::numeric_limits<float>::quiet_NaN();
3719
3720
3721 int vertexID = track->get_vertex_id();
3722 GlobalVertex* vertex = gvertexmap->get(vertexID);
3723 float vx = std::numeric_limits<float>::quiet_NaN();
3724 float vy = std::numeric_limits<float>::quiet_NaN();
3725 float vz = std::numeric_limits<float>::quiet_NaN();
3726 if (vertex)
3727 {
3728 vx = vertex->get_x();
3729 vy = vertex->get_y();
3730 vz = vertex->get_z();
3731 Acts::Vector3 vert(vx, vy, vz);
3732
3733 auto dcapair = TrackAnalysisUtils::get_dca(track, vert);
3734 dca3dxy = dcapair.first.first;
3735 dca3dxysigma = dcapair.first.second;
3736 dca3dz = dcapair.second.first;
3737 dca3dzsigma = dcapair.second.second;
3738 }
3739 float px = track->get_px();
3740 float py = track->get_py();
3741 float pz = track->get_pz();
3742 TVector3 v(px, py, pz);
3743 float pt = v.Pt();
3744 float eta = v.Eta();
3745 float phi = v.Phi();
3746 float CVxx = track->get_error(3, 3);
3747 float CVxy = track->get_error(3, 4);
3748 float CVxz = track->get_error(3, 5);
3749 float CVyy = track->get_error(4, 4);
3750 float CVyz = track->get_error(4, 5);
3751 float CVzz = track->get_error(5, 5);
3752 float deltapt = std::sqrt((CVxx * px * px + 2 * CVxy * px * py + CVyy * py * py) / (px * px + py * py));
3753 float deltaeta = std::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)));
3754 float deltaphi = std::sqrt((CVyy * px * px - 2 * CVxy * px * py + CVxx * py * py) / ((px * px + py * py) * (px * px + py * py)));
3755 float pcax = track->get_x();
3756 float pcay = track->get_y();
3757 float pcaz = track->get_z();
3758
3759 float gtrackID = std::numeric_limits<float>::quiet_NaN();
3760 float gflavor = std::numeric_limits<float>::quiet_NaN();
3761 float ng4hits = std::numeric_limits<float>::quiet_NaN();
3762 unsigned int ngmaps = 0;
3763 unsigned int ngintt = 0;
3764 unsigned int ngmms = 0;
3765 unsigned int ngtpc = 0;
3766 unsigned int nglmaps = 0;
3767 unsigned int nglintt = 0;
3768 unsigned int nglmms = 0;
3769 unsigned int ngltpc = 0;
3770 float gpx = std::numeric_limits<float>::quiet_NaN();
3771 float gpy = std::numeric_limits<float>::quiet_NaN();
3772 float gpt = std::numeric_limits<float>::quiet_NaN();
3773 float geta = std::numeric_limits<float>::quiet_NaN();
3774 float gphi = std::numeric_limits<float>::quiet_NaN();
3775 float gpz = std::numeric_limits<float>::quiet_NaN();
3776 float gvx = std::numeric_limits<float>::quiet_NaN();
3777 float gvy = std::numeric_limits<float>::quiet_NaN();
3778 float gvz = std::numeric_limits<float>::quiet_NaN();
3779 float gvt = std::numeric_limits<float>::quiet_NaN();
3780 float gfpx = std::numeric_limits<float>::quiet_NaN();
3781 float gfpy = std::numeric_limits<float>::quiet_NaN();
3782 float gfpz = std::numeric_limits<float>::quiet_NaN();
3783 float gfx = std::numeric_limits<float>::quiet_NaN();
3784 float gfy = std::numeric_limits<float>::quiet_NaN();
3785 float gfz = std::numeric_limits<float>::quiet_NaN();
3786 float gembed = std::numeric_limits<float>::quiet_NaN();
3787 float gprimary = std::numeric_limits<float>::quiet_NaN();
3788
3789 int ispure = 0;
3790 float nfromtruth = std::numeric_limits<float>::quiet_NaN();
3791 float nwrong = std::numeric_limits<float>::quiet_NaN();
3792 float ntrumaps = std::numeric_limits<float>::quiet_NaN();
3793 float nwrongmaps = std::numeric_limits<float>::quiet_NaN();
3794 float ntruintt = std::numeric_limits<float>::quiet_NaN();
3795 float nwrongintt = std::numeric_limits<float>::quiet_NaN();
3796 float ntrumms = std::numeric_limits<float>::quiet_NaN();
3797 float nwrongmms = std::numeric_limits<float>::quiet_NaN();
3798 float ntrutpc = std::numeric_limits<float>::quiet_NaN();
3799 float nwrongtpc = std::numeric_limits<float>::quiet_NaN();
3800 float ntrutpc1 = std::numeric_limits<float>::quiet_NaN();
3801 float nwrongtpc1 = std::numeric_limits<float>::quiet_NaN();
3802 float ntrutpc11 = std::numeric_limits<float>::quiet_NaN();
3803 float nwrongtpc11 = std::numeric_limits<float>::quiet_NaN();
3804 float ntrutpc2 = std::numeric_limits<float>::quiet_NaN();
3805 float nwrongtpc2 = std::numeric_limits<float>::quiet_NaN();
3806 float ntrutpc3 = std::numeric_limits<float>::quiet_NaN();
3807 float nwrongtpc3 = std::numeric_limits<float>::quiet_NaN();
3808 float layersfromtruth = std::numeric_limits<float>::quiet_NaN();
3809
3810 if (_do_track_match)
3811 {
3812 PHG4Particle* g4particle = trackeval->max_truth_particle_by_nclusters(track);
3813 if (g4particle)
3814 {
3815 if (_scan_for_embedded)
3816 {
3817 if (trutheval->get_embed(g4particle) <= 0)
3818 {
3819 continue;
3820 }
3821 }
3822 SvtxTrack* truthrecotrk = trackeval->best_track_from(g4particle);
3823 if (truthrecotrk)
3824 {
3825 if (truthrecotrk->get_id() == track->get_id())
3826 {
3827 ispure = 1;
3828 }
3829 }
3830 gtrackID = g4particle->get_track_id();
3831 gflavor = g4particle->get_pid();
3832
3833 std::set<TrkrDefs::cluskey> g4clusters = clustereval->all_clusters_from(g4particle);
3834 ng4hits = g4clusters.size();
3835 gpx = g4particle->get_px();
3836 gpy = g4particle->get_py();
3837 gpz = g4particle->get_pz();
3838
3839 std::vector<int> lmaps(_nlayers_maps + 1, 0);
3840 std::vector<int> lintt(_nlayers_intt + 1, 0);
3841 std::vector<int> ltpc(_nlayers_tpc + 1, 0);
3842 std::vector<int> lmms(_nlayers_mms + 1, 0);
3843
3844 for (const TrkrDefs::cluskey g4cluster : g4clusters)
3845 {
3846 unsigned int local_layer = TrkrDefs::getLayer(g4cluster);
3847 if (_nlayers_maps > 0 && local_layer < _nlayers_maps)
3848 {
3849 lmaps[local_layer] = 1;
3850 ngmaps++;
3851 }
3852
3853 if (_nlayers_intt > 0 && local_layer >= _nlayers_maps && local_layer < _nlayers_maps + _nlayers_intt)
3854 {
3855 lintt[local_layer - _nlayers_maps] = 1;
3856 ngintt++;
3857 }
3858
3859 if (_nlayers_tpc > 0 && local_layer >= _nlayers_maps + _nlayers_intt && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
3860 {
3861 ltpc[local_layer - (_nlayers_maps + _nlayers_intt)] = 1;
3862 ngtpc++;
3863 }
3864
3865 if (_nlayers_mms > 0 && local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc && local_layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
3866 {
3867 lmms[local_layer - (_nlayers_maps + _nlayers_intt + _nlayers_tpc)] = 1;
3868 ngmms++;
3869 }
3870 }
3871 if (_nlayers_maps > 0)
3872 {
3873 for (unsigned int i = 0; i < _nlayers_maps; i++)
3874 {
3875 nglmaps += lmaps[i];
3876 }
3877 }
3878 if (_nlayers_intt > 0)
3879 {
3880 for (unsigned int i = 0; i < _nlayers_intt; i++)
3881 {
3882 nglintt += lintt[i];
3883 }
3884 }
3885 if (_nlayers_tpc > 0)
3886 {
3887 for (unsigned int i = 0; i < _nlayers_tpc; i++)
3888 {
3889 ngltpc += ltpc[i];
3890 }
3891 }
3892 if (_nlayers_mms > 0)
3893 {
3894 for (unsigned int i = 0; i < _nlayers_mms; i++)
3895 {
3896 nglmms += lmms[i];
3897 }
3898 }
3899
3900 TVector3 gv(gpx, gpy, gpz);
3901 gpt = gv.Pt();
3902 geta = gv.Eta();
3903 gphi = gv.Phi();
3904 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
3905 gvx = vtx->get_x();
3906 gvy = vtx->get_y();
3907 gvz = vtx->get_z();
3908 gvt = vtx->get_t();
3909
3910 PHG4Hit* outerhit = nullptr;
3911 if (_do_eval_light == false)
3912 {
3913 outerhit = trutheval->get_outermost_truth_hit(g4particle);
3914 }
3915 if (outerhit)
3916 {
3917 gfpx = outerhit->get_px(1);
3918 gfpy = outerhit->get_py(1);
3919 gfpz = outerhit->get_pz(1);
3920 gfx = outerhit->get_x(1);
3921 gfy = outerhit->get_y(1);
3922 gfz = outerhit->get_z(1);
3923 }
3924 gembed = trutheval->get_embed(g4particle);
3925 gprimary = trutheval->is_primary(g4particle);
3926
3927 nfromtruth = trackeval->get_nclusters_contribution(track, g4particle);
3928 nwrong = trackeval->get_nwrongclusters_contribution(track, g4particle);
3929 if (_nlayers_maps == 0)
3930 {
3931 ntrumaps = 0;
3932 }
3933 else
3934 {
3935 auto pair = trackeval->get_layer_range_contribution(track, g4particle, 0, _nlayers_maps);
3936 ntrumaps = pair.first;
3937 nwrongmaps = pair.second;
3938 }
3939 if (_nlayers_intt == 0)
3940 {
3941 ntruintt = 0;
3942 }
3943 else
3944 {
3945 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps, _nlayers_maps + _nlayers_intt);
3946 ntruintt = pair.first;
3947 nwrongintt = pair.second;
3948 }
3949 if (_nlayers_mms == 0)
3950 {
3951 ntrumms = 0;
3952 }
3953 else
3954 {
3955 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + _nlayers_tpc, _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms);
3956 ntrumms = pair.first;
3957 nwrongmms = pair.second;
3958 }
3959 auto pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
3960 ntrutpc = pair.first;
3961 nwrongtpc = pair.second;
3962 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 16);
3963 ntrutpc1 = pair.first;
3964 nwrongtpc1 = pair.second;
3965 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt, _nlayers_maps + _nlayers_intt + 8);
3966 ntrutpc11 = pair.first;
3967 nwrongtpc11 = pair.second;
3968 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + 16, _nlayers_maps + _nlayers_intt + 32);
3969 ntrutpc2 = pair.first;
3970 nwrongtpc2 = pair.second;
3971 pair = trackeval->get_layer_range_contribution(track, g4particle, _nlayers_maps + _nlayers_intt + 32, _nlayers_maps + _nlayers_intt + _nlayers_tpc);
3972 ntrutpc3 = pair.first;
3973 nwrongtpc3 = pair.second;
3974 layersfromtruth = trackeval->get_nclusters_contribution_by_layer(track, g4particle);
3975 }
3976 }
3977 if (Verbosity() > 2)
3978 {
3979 std::cout << " npedge " << npedge
3980 << " nredge " << nredge
3981 << " nbig " << nbig
3982 << " novlp " << novlp
3983 << std::endl;
3984 }
3985 float track_data[] = {(float) _ievent, m_fSeed,
3986 trackID,
3987 crossing,
3988 px,
3989 py,
3990 pz,
3991 pt,
3992 eta,
3993 phi,
3994 deltapt,
3995 deltaeta,
3996 deltaphi,
3997 siqr,
3998 siphi,
3999 sithe,
4000 six0,
4001 siy0,
4002 tpqr,
4003 tpphi,
4004 tpthe,
4005 tpx0,
4006 tpy0,
4007 charge,
4008 quality,
4009 chisq,
4010 ndf,
4011 local_nhits, nmaps, nintt, ntpc, nmms,
4012 ntpc1, ntpc11, ntpc2, ntpc3,
4013 nlmaps, nlintt, nltpc, nlmms,
4014 (float) layers,
4015 (float) vertexID,
4016 vx,
4017 vy,
4018 vz,
4019 dca2d,
4020 dca2dsigma,
4021 dca3dxy,
4022 dca3dxysigma,
4023 dca3dz,
4024 dca3dzsigma,
4025 pcax,
4026 pcay,
4027 pcaz,
4028 gtrackID,
4029 (float) ispure,
4030 gflavor,
4031 ng4hits,
4032 (float) ngmaps,
4033 (float) ngintt,
4034 (float) ngtpc,
4035 (float) ngmms,
4036 (float) nglmaps,
4037 (float) nglintt,
4038 (float) ngltpc,
4039 (float) nglmms,
4040 gpx,
4041 gpy,
4042 gpz,
4043 gpt,
4044 geta,
4045 gphi,
4046 gvx,
4047 gvy,
4048 gvz,
4049 gvt,
4050 gfpx,
4051 gfpy,
4052 gfpz,
4053 gfx,
4054 gfy,
4055 gfz,
4056 gembed,
4057 gprimary,
4058 nfromtruth,
4059 nwrong,
4060 ntrumaps,
4061 nwrongmaps,
4062 ntruintt,
4063 nwrongintt,
4064 ntrutpc,
4065 nwrongtpc,
4066 ntrumms,
4067 nwrongmms,
4068 ntrutpc1,
4069 nwrongtpc1,
4070 ntrutpc11,
4071 nwrongtpc11,
4072 ntrutpc2,
4073 nwrongtpc2,
4074 ntrutpc3,
4075 nwrongtpc3,
4076 layersfromtruth,
4077 npedge,
4078 nredge,
4079 nbig,
4080 novlp,
4081 merr,
4082 msize,
4083 nhit_tpc_all,
4084 nhit_tpc_in,
4085 nhit_tpc_mid,
4086 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
4087
4088 if (Verbosity() >= 1)
4089 {
4090 std::cout << "ievent " << _ievent
4091 << " trackID " << trackID
4092 << " nhits " << local_nhits
4093 << " px " << px
4094 << " py " << py
4095 << " pz " << pz
4096 << " gembed " << gembed
4097 << " gprimary " << gprimary
4098 << std::endl;
4099 }
4100
4101 _ntp_track->Fill(track_data);
4102 }
4103 }
4104 if (Verbosity() > 1)
4105 {
4106 _timer->stop();
4107 std::cout << "track time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
4108 }
4109 }
4110
4111
4112
4113
4114
4115 if (_ntp_gseed)
4116 {
4117 if (Verbosity() > 1)
4118 {
4119 std::cout << "Filling ntp_gseed " << std::endl;
4120 _timer->restart();
4121 }
4122
4123 PHG4TruthInfoContainer* truthinfo = findNode::getClass<PHG4TruthInfoContainer>(topNode, "G4TruthInfo");
4124
4125 float gx = std::numeric_limits<float>::quiet_NaN();
4126 float gy = std::numeric_limits<float>::quiet_NaN();
4127 float gz = std::numeric_limits<float>::quiet_NaN();
4128 float gr = std::numeric_limits<float>::quiet_NaN();
4129 float geta = std::numeric_limits<float>::quiet_NaN();
4130 float gphi = std::numeric_limits<float>::quiet_NaN();
4131 float glayer = std::numeric_limits<float>::quiet_NaN();
4132 float gpx = std::numeric_limits<float>::quiet_NaN();
4133 float gpy = std::numeric_limits<float>::quiet_NaN();
4134 float gpz = std::numeric_limits<float>::quiet_NaN();
4135 float gtpt = std::numeric_limits<float>::quiet_NaN();
4136 float gtphi = std::numeric_limits<float>::quiet_NaN();
4137 float gteta = std::numeric_limits<float>::quiet_NaN();
4138 float gvx = std::numeric_limits<float>::quiet_NaN();
4139 float gvy = std::numeric_limits<float>::quiet_NaN();
4140 float gvz = std::numeric_limits<float>::quiet_NaN();
4141 float gembed = std::numeric_limits<float>::quiet_NaN();
4142 float gprimary = std::numeric_limits<float>::quiet_NaN();
4143 float gflav = std::numeric_limits<float>::quiet_NaN();
4144 float dphiprev = std::numeric_limits<float>::quiet_NaN();
4145 float detaprev = std::numeric_limits<float>::quiet_NaN();
4146
4147 float* xval = new float[_nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms];
4148 float* yval = new float[_nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms];
4149 float* zval = new float[_nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms];
4150 if (truthinfo)
4151 {
4152 int ntrk = 0;
4153 PHG4TruthInfoContainer::ConstRange range = truthinfo->GetPrimaryParticleRange();
4154 for (PHG4TruthInfoContainer::ConstIterator iter = range.first;
4155 iter != range.second;
4156 ++iter)
4157 {
4158 ntrk++;
4159 PHG4Particle* g4particle = iter->second;
4160 for (unsigned int i = 0; i < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms; i++)
4161 {
4162 xval[i] = 0;
4163 yval[i] = 0;
4164 zval[i] = 0;
4165 }
4166 std::set<PHG4Hit*> truth_hits = trutheval->all_truth_hits(g4particle);
4167 for (auto* g4hit : truth_hits)
4168 {
4169 unsigned int local_layer = g4hit->get_layer();
4170
4171 if (local_layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms)
4172 {
4173
4174 continue;
4175 }
4176 xval[local_layer] = g4hit->get_avg_x();
4177 yval[local_layer] = g4hit->get_avg_y();
4178 zval[local_layer] = g4hit->get_avg_z();
4179 }
4180
4181 for (unsigned int i = 0; i < _nlayers_maps + _nlayers_intt + _nlayers_tpc + _nlayers_mms; i++)
4182 {
4183 gx = xval[i];
4184 gy = yval[i];
4185 gz = zval[i];
4186 if (gx == 0 && gy == 0)
4187 {
4188 continue;
4189 }
4190
4191 TVector3 vg4(gx, gy, gz);
4192 glayer = i;
4193 gr = vg4.Perp();
4194 geta = vg4.Eta();
4195 gphi = vg4.Phi();
4196 gpx = g4particle->get_px();
4197 gpy = g4particle->get_py();
4198 gpz = g4particle->get_pz();
4199 TVector3 vg4p(gpx, gpy, gpz);
4200
4201 gtpt = vg4p.Perp();
4202 gtphi = vg4p.Phi();
4203 gteta = vg4p.Eta();
4204
4205 PHG4VtxPoint* vtx = trutheval->get_vertex(g4particle);
4206
4207 if (vtx)
4208 {
4209 gvx = vtx->get_x();
4210 gvy = vtx->get_y();
4211 gvz = vtx->get_z();
4212 }
4213
4214 gembed = trutheval->get_embed(g4particle);
4215 gprimary = trutheval->is_primary(g4particle);
4216 gflav = g4particle->get_pid();
4217 if (i >= 1)
4218 {
4219 if (xval[i - 1] != 0 && yval[i - 1] != 0)
4220 {
4221 TVector3 vg4prev(xval[i - 1], yval[i - 1], zval[i - 1]);
4222 dphiprev = vg4.DeltaPhi(vg4prev);
4223 detaprev = geta - vg4prev.Eta();
4224 }
4225 }
4226
4227 float ntrk_f = ntrk;
4228 float _ievent_f = _ievent;
4229 float gseed_data[] = {_ievent_f,
4230 ntrk_f,
4231 gx,
4232 gy,
4233 gz,
4234 gr,
4235 geta,
4236 gphi,
4237 glayer,
4238 gpx,
4239 gpy,
4240 gpz,
4241 gtpt,
4242 gtphi,
4243 gteta,
4244 gvx,
4245 gvy,
4246 gvz,
4247 gembed,
4248 gprimary,
4249 gflav,
4250 dphiprev,
4251 detaprev,
4252 nhit_tpc_all,
4253 nhit_tpc_in,
4254 nhit_tpc_mid,
4255 nhit_tpc_out, nclus_all, nclus_tpc, nclus_intt, nclus_maps, nclus_mms};
4256
4257 _ntp_gseed->Fill(gseed_data);
4258 }
4259 delete[] xval;
4260 delete[] yval;
4261 delete[] zval;
4262 }
4263 }
4264
4265 if (Verbosity() > 1)
4266 {
4267 _timer->stop();
4268 std::cout << "g4hit time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
4269 }
4270 }
4271 return;
4272 }
4273
4274 TMatrixF SvtxEvaluator::calculateClusterError(TrkrCluster* c, float& clusphi)
4275 {
4276 TMatrixF localErr(3, 3);
4277 localErr[0][0] = 0.;
4278 localErr[0][1] = 0.;
4279 localErr[0][2] = 0.;
4280 localErr[1][0] = 0.;
4281 localErr[1][1] = c->getActsLocalError(0, 0);
4282 localErr[1][2] = c->getActsLocalError(0, 1);
4283 localErr[2][0] = 0.;
4284 localErr[2][1] = c->getActsLocalError(1, 0);
4285 localErr[2][2] = c->getActsLocalError(1, 1);
4286
4287 TMatrixF ROT(3, 3);
4288 ROT[0][0] = std::cos(clusphi);
4289 ROT[0][1] = -std::sin(clusphi);
4290 ROT[0][2] = 0.0;
4291 ROT[1][0] = std::sin(clusphi);
4292 ROT[1][1] = std::cos(clusphi);
4293 ROT[1][2] = 0.0;
4294 ROT[2][0] = 0.0;
4295 ROT[2][1] = 0.0;
4296 ROT[2][2] = 1.0;
4297 TMatrixF ROT_T(3, 3);
4298 ROT_T.Transpose(ROT);
4299
4300 TMatrixF err(3, 3);
4301 err = ROT * localErr * ROT_T;
4302 return err;
4303 }
