File indexing completed on 2025-08-06 08:18:21
0001 #include "TrkrNtuplizer.h"
0002
0003 #include <trackbase/ActsGeometry.h>
0004 #include <trackbase/ClusterErrorPara.h>
0005 #include <trackbase/TrackFitUtils.h>
0006 #include <trackbase/TrkrCluster.h>
0007 #include <trackbase/TrkrDefs.h>
0008 #include <trackbase/TrkrHit.h>
0009 #include <trackbase/TrkrHitSetContainer.h>
0010 #include <trackbase/InttDefs.h>
0011 #include <trackbase/MvtxDefs.h>
0012 #include <trackbase/TpcDefs.h>
0013 #include <trackbase/TrkrClusterContainer.h>
0014 #include <trackbase/TrkrClusterHitAssoc.h>
0015 #include <trackbase/TrkrClusterIterationMapv1.h>
0016 #include <trackbase/TrkrHitSet.h>
0017
0018 #include <trackbase_historic/TrackSeedContainer_v1.h>
0019 #include <trackbase_historic/ActsTransformations.h>
0020 #include <trackbase_historic/SvtxTrack.h>
0021 #include <trackbase_historic/SvtxTrackMap.h>
0022 #include <trackbase_historic/TrackAnalysisUtils.h>
0023 #include <trackbase_historic/TrackSeed.h>
0024
0025 #include <micromegas/MicromegasDefs.h>
0026
0027
0028 #include <globalvertex/GlobalVertex.h>
0029 #include <globalvertex/GlobalVertexMap.h>
0030 #include <globalvertex/SvtxVertex.h>
0031 #include <globalvertex/SvtxVertexMap.h>
0032
0033
0034 #include <g4detectors/PHG4TpcCylinderGeom.h>
0035 #include <g4detectors/PHG4TpcCylinderGeomContainer.h>
0036
0037 #include <trackermillepedealignment/HelicalFitter.h>
0038
0039 #include <ffamodules/CDBInterface.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 <cmath>
0054 #include <iomanip>
0055 #include <iostream>
0056 #include <iterator>
0057 #include <map>
0058 #include <memory> // for shared_ptr
0059 #include <set> // for _Rb_tree_cons...
0060 #include <utility>
0061 #include <vector>
0062
0063 enum n_event
0064 {
0065 evnev,
0066 evnseed,
0067 evnrun,
0068 evnseg,
0069 evnjob,
0070 evsize = evnjob + 1,
0071 };
0072
0073 enum n_info
0074 {
0075 infonocc11,
0076 infonocc116,
0077 infonocc21,
0078 infonocc216,
0079 infonocc31,
0080 infonocc316,
0081 infontrk,
0082 infonntpcseed,
0083 infonnsiseed,
0084 infonhitmvtx,
0085 infonhitintt,
0086 infonhittpot,
0087 infonhittpcall,
0088 infonhittpcin,
0089 infonhittpcmid,
0090 infonhittpcout,
0091 infonclusall,
0092 infonclustpc,
0093 infonclustpcpos,
0094 infonclustpcneg,
0095 infonclusintt,
0096 infonclusmvtx,
0097 infonclustpot,
0098 infosize = infonclustpot + 1
0099 };
0100
0101 enum n_vertex
0102 {
0103 vtxnvertexID,
0104 vtxnvx,
0105 vtxnvy,
0106 vtxnvz,
0107 vtxnntracks,
0108 vtxnchi2,
0109 vtxnndof,
0110 vtxsize = vtxnndof + 1
0111 };
0112
0113 enum n_hit
0114 {
0115 nhitID,
0116 nhite,
0117 nhitadc,
0118 nhitlayer,
0119 nhitphielem,
0120 nhitzelem,
0121 nhitcellID,
0122 nhitecell,
0123 nhitphibin,
0124 nhittbin,
0125 nhitphi,
0126 nhitr,
0127 nhitx,
0128 nhity,
0129 nhitz,
0130 hitsize = nhitz + 1
0131 };
0132
0133 enum n_seed
0134 {
0135 nseedtrackID,
0136 nseedniter,
0137 nseedpt,
0138 nseedptot,
0139 nseedeta,
0140 nseedphi,
0141 nseedsyxint,
0142 nseedsrzint,
0143 nseedsxyslope,
0144 nseedsrzslope,
0145 nseedX0,
0146 nseedY0,
0147 nseedZ0,
0148 nseedR0,
0149 nseedcharge,
0150 nseeddedx,
0151 nseedpidedx,
0152 nseedkdedx,
0153 nseedprdedx,
0154 nseedn1pix,
0155 nseednhits,
0156 seedsize = nseednhits + 1
0157 };
0158
0159 enum n_residual
0160 {
0161 nresalpha,
0162 nresbeta,
0163 nresphio,
0164 nresphi,
0165 nresz,
0166 ressize = nresz + 1
0167 };
0168
0169 enum n_track
0170 {
0171 ntrktrackID,
0172 ntrkcrossing,
0173 ntrkpx,
0174 ntrkpy,
0175 ntrkpz,
0176 ntrkpt,
0177 ntrketa,
0178 ntrkphi,
0179 ntrkdeltapt,
0180 ntrkdeltaeta,
0181 ntrkdeltaphi,
0182 ntrkcharge,
0183 ntrkquality,
0184 ntrkchisq,
0185 ntrkndf,
0186 ntrknhits,
0187 ntrknmaps,
0188 ntrknintt,
0189 ntrkntpc,
0190 ntrknmms,
0191 ntrkntpc1,
0192 ntrkntpc11,
0193 ntrkntpc2,
0194 ntrkntpc3,
0195 ntrkndedx,
0196 ntrknpidedx,
0197 ntrknkdedx,
0198 ntrknprdedx,
0199 ntrkvertexID,
0200 ntrkvx,
0201 ntrkvy,
0202 ntrkvz,
0203 ntrkdca2d,
0204 ntrkdca2dsigma,
0205 ntrkdca3dxy,
0206 ntrkdca3dxysigma,
0207 ntrkdca3dz,
0208 ntrkdca3dzsigma,
0209 ntrkpcax,
0210 ntrkpcay,
0211 ntrkpcaz,
0212 ntrkhlxpt,
0213 ntrkhlxeta,
0214 ntrkhlxphi,
0215 ntrkhlxX0,
0216 ntrkhlxY0,
0217 ntrkhlxZ0,
0218 ntrkhlxcharge,
0219 trksize = ntrkhlxcharge + 1
0220 };
0221
0222 enum n_cluster
0223 {
0224 nclulocx,
0225 nclulocy,
0226 nclux,
0227 ncluy,
0228 ncluz,
0229 nclur,
0230 ncluphi,
0231 nclueta,
0232 nclutheta,
0233 ncluphibin,
0234 nclutbin,
0235 ncluex,
0236 ncluey,
0237 ncluez,
0238 ncluephi,
0239 nclupez,
0240 nclupephi,
0241 nclue,
0242 ncluadc,
0243 nclumaxadc,
0244 ncluthick,
0245 ncluafac,
0246 nclubfac,
0247 ncludcal,
0248 nclulayer,
0249 ncluphielem,
0250 ncluzelem,
0251 nclusize,
0252 ncluphisize,
0253 ncluzsize,
0254 nclupedge,
0255 ncluredge,
0256 ncluovlp,
0257 nclutrackID,
0258 ncluniter,
0259 clusize = ncluniter + 1
0260 };
0261
0262 TrkrNtuplizer::TrkrNtuplizer(const std::string& , const std::string& filename, const std::string& trackmapname,
0263 unsigned int nlayers_maps,
0264 unsigned int nlayers_intt,
0265 unsigned int nlayers_tpc,
0266 unsigned int nlayers_mms)
0267 : SubsysReco("TrkrNtuplizer")
0268 , _nlayers_maps(nlayers_maps)
0269 , _nlayers_intt(nlayers_intt)
0270 , _nlayers_tpc(nlayers_tpc)
0271 , _nlayers_mms(nlayers_mms)
0272 , _filename(filename)
0273 , _trackmapname(trackmapname)
0274 {
0275 }
0276
0277 TrkrNtuplizer::~TrkrNtuplizer()
0278 {
0279 delete _timer;
0280 }
0281
0282 int TrkrNtuplizer::Init(PHCompositeNode* )
0283 {
0284 _ievent = 0;
0285
0286 _tfile = new TFile(_filename.c_str(), "RECREATE");
0287 _tfile->SetCompressionLevel(7);
0288 std::string str_vertex = {"vertexID:vx:vy:vz:ntracks:chi2:ndof"};
0289 std::string str_event = {"event:seed:run:seg:job"};
0290 std::string str_hit = {"hitID:e:adc:layer:phielem:zelem:cellID:ecell:phibin:tbin:phi:r:x:y:z"};
0291 std::string str_cluster = {"locx:locy:x:y:z:r:phi:eta:theta:phibin:tbin:ex:ey:ez:ephi:pez:pephi:e:adc:maxadc:thick:afac:bfac:dcal:layer:phielem:zelem:size:phisize:zsize:pedge:redge:ovlp:trackID:niter"};
0292 std::string str_seed = {"seedID:siter:spt:sptot:seta:sphi:syxint:srzint:sxyslope:srzslope:sX0:sY0:sdZ0:sR0:scharge:sdedx:spidedx:skdedx:sprdedx:sn1pix:snhits"};
0293 std::string str_residual = {"alpha:beta:resphio:resphi:resz"};
0294 std::string str_track = {"trackID:crossing:px:py:pz:pt:eta:phi:deltapt:deltaeta:deltaphi:charge:quality:chisq:ndf:nhits:nmaps:nintt:ntpc:nmms:ntpc1:ntpc11:ntpc2:ntpc3:dedx:pidedx:kdedx:prdedx:nlmaps:nlintt:nltpc:nlmms:layers:vertexID:vx:vy:vz:dca2d:dca2dsigma:dca3dxy:dca3dxysigma:dca3dz:dca3dzsigma:pcax:pcay:pcaz:hlxpt:hlxeta:hlxphi:hlxX0:hlxY0:hlxZ0:hlxcharge"};
0295 std::string str_info = {"occ11:occ116:occ21:occ216:occ31:occ316:ntrk:ntpcseed:nsiseed:nhitmvtx:nhitintt:nhittpot:nhittpcall:nhittpcin:nhittpcmid:nhittpcout:nclusall:nclustpc:nclustpcpos:nclustpcneg:nclusintt:nclusmaps:nclusmms"};
0296
0297 if (_do_info_eval)
0298 {
0299 std::string ntp_varlist_info = str_event + ":" + str_info;
0300 _ntp_info = new TNtuple("ntp_info", "event info", ntp_varlist_info.c_str());
0301 }
0302
0303 if (_do_vertex_eval)
0304 {
0305 std::string ntp_varlist_vtx = str_event + ":" + str_vertex + ":" + str_info;
0306 _ntp_vertex = new TNtuple("ntp_vertex", "vertex => max truth", ntp_varlist_vtx.c_str());
0307 }
0308
0309 if (_do_hit_eval)
0310 {
0311 std::string ntp_varlist_ev = str_event + ":" + str_hit + ":" + str_info;
0312 _ntp_hit = new TNtuple("ntp_hit", "svtxhit => max truth", ntp_varlist_ev.c_str());
0313 }
0314
0315 if (_do_cluster_eval)
0316 {
0317 std::string ntp_varlist_clu = str_event + ":" + str_cluster + ":" + str_info;
0318 _ntp_cluster = new TNtuple("ntp_cluster", "svtxcluster => max truth", ntp_varlist_clu.c_str());
0319 }
0320 if (_do_clus_trk_eval)
0321 {
0322 std::string ntp_varlist_clut = str_event + ":" + str_cluster + ":" + str_residual + ":" + str_seed + ":" + str_info;
0323 _ntp_clus_trk = new TNtuple("ntp_clus_trk", "cluster on track", ntp_varlist_clut.c_str());
0324 }
0325
0326 if (_do_track_eval)
0327 {
0328 std::string ntp_varlist_trk = str_event + ":" + str_track + ":" + str_info;
0329 _ntp_track = new TNtuple("ntp_track", "svtxtrack => max truth", ntp_varlist_trk.c_str());
0330 }
0331
0332 if (_do_tpcseed_eval)
0333 {
0334 std::string ntp_varlist_tsee = str_event + ":" + str_seed + ":" + str_info;
0335 _ntp_tpcseed = new TNtuple("ntp_tpcseed", "seeds from truth", ntp_varlist_tsee.c_str());
0336 }
0337 if (_do_siseed_eval)
0338 {
0339 std::string ntp_varlist_ssee = str_event + ":" + str_seed + ":" + str_info;
0340 _ntp_siseed = new TNtuple("ntp_siseed", "seeds from truth", ntp_varlist_ssee.c_str());
0341 }
0342
0343 std::string dedx_fitparams = CDBInterface::instance()->getUrl("TPC_DEDX_FITPARAM");
0344 TFile *filefit = TFile::Open(dedx_fitparams.c_str());
0345
0346 if (!filefit->IsOpen())
0347 {
0348 std::cout << "Error opening filefit!" << std::endl;
0349 return Fun4AllReturnCodes::ABORTEVENT;
0350 }
0351
0352 filefit->GetObject("f_piband", f_pion_plus);
0353 filefit->GetObject("f_Kband", f_kaon_plus);
0354 filefit->GetObject("f_pband", f_proton_plus);
0355 filefit->GetObject("f_piminus_band", f_pion_minus);
0356 filefit->GetObject("f_Kminus_band", f_kaon_minus);
0357 filefit->GetObject("f_pbar_band", f_proton_minus);
0358
0359 _timer = new PHTimer("_eval_timer");
0360 _timer->stop();
0361
0362
0363 dedxcorr[0][0][0] = 0.54;
0364 dedxcorr[0][0][1] = 0.90;
0365 dedxcorr[0][0][2] = 1.11;
0366 dedxcorr[0][1][0] = 0.72;
0367 dedxcorr[0][1][1] = 0.98;
0368 dedxcorr[0][1][2] = 0.99;
0369 dedxcorr[0][2][0] = 0.83;
0370 dedxcorr[0][2][1] = 0.81;
0371 dedxcorr[0][2][2] = 1.03;
0372 dedxcorr[0][3][0] = 0.82;
0373 dedxcorr[0][3][1] = 0.52;
0374 dedxcorr[0][3][2] = 1.02;
0375 dedxcorr[0][4][0] = 0.85;
0376 dedxcorr[0][4][1] = 0.78;
0377 dedxcorr[0][4][2] = 0.82;
0378 dedxcorr[0][5][0] = 0.85;
0379 dedxcorr[0][5][1] = 0.70;
0380 dedxcorr[0][5][2] = 0.92;
0381 dedxcorr[0][6][0] = 0.78;
0382 dedxcorr[0][6][1] = 0.94;
0383 dedxcorr[0][6][2] = 1.03;
0384 dedxcorr[0][7][0] = 0.64;
0385 dedxcorr[0][7][1] = 0.94;
0386 dedxcorr[0][7][2] = 1.19;
0387 dedxcorr[0][8][0] = 0.95;
0388 dedxcorr[0][8][1] = 0.47;
0389 dedxcorr[0][8][2] = 1.12;
0390 dedxcorr[0][9][0] = 0.74;
0391 dedxcorr[0][9][1] = 1.07;
0392 dedxcorr[0][9][2] = 0.85;
0393 dedxcorr[0][10][0] = 0.86;
0394 dedxcorr[0][10][1] = 0.74;
0395 dedxcorr[0][10][2] = 1.08;
0396 dedxcorr[0][11][0] = 0.60;
0397 dedxcorr[0][11][1] = 1.03;
0398 dedxcorr[0][11][2] = 0.89;
0399 dedxcorr[1][0][0] = 0.74;
0400 dedxcorr[1][0][1] = 0.76;
0401 dedxcorr[1][0][2] = 0.89;
0402 dedxcorr[1][1][0] = 0.61;
0403 dedxcorr[1][1][1] = 0.87;
0404 dedxcorr[1][1][2] = 1.10;
0405 dedxcorr[1][2][0] = 0.73;
0406 dedxcorr[1][2][1] = 0.86;
0407 dedxcorr[1][2][2] = 1.70;
0408 dedxcorr[1][3][0] = 0.40;
0409 dedxcorr[1][3][1] = 0.76;
0410 dedxcorr[1][3][2] = 0.99;
0411 dedxcorr[1][4][0] = 0.82;
0412 dedxcorr[1][4][1] = 0.71;
0413 dedxcorr[1][4][2] = 0.86;
0414 dedxcorr[1][5][0] = 0.65;
0415 dedxcorr[1][5][1] = 1.00;
0416 dedxcorr[1][5][2] = 0.79;
0417 dedxcorr[1][6][0] = 0.75;
0418 dedxcorr[1][6][1] = 0.84;
0419 dedxcorr[1][6][2] = 0.40;
0420 dedxcorr[1][7][0] = 0.95;
0421 dedxcorr[1][7][1] = 0.81;
0422 dedxcorr[1][7][2] = 0.85;
0423 dedxcorr[1][8][0] = 0.79;
0424 dedxcorr[1][8][1] = 0.67;
0425 dedxcorr[1][8][2] = 1.17;
0426 dedxcorr[1][9][0] = 0.73;
0427 dedxcorr[1][9][1] = 0.94;
0428 dedxcorr[1][9][2] = 0.95;
0429 dedxcorr[1][10][0] = 0.40;
0430 dedxcorr[1][10][1] = 0.83;
0431 dedxcorr[1][10][2] = 0.87;
0432 dedxcorr[1][11][0] = 0.80;
0433 dedxcorr[1][11][1] = 0.63;
0434 dedxcorr[1][11][2] = 0.68;
0435 dedxcorr[0][0][1] *= 0.88 ;
0436 dedxcorr[0][0][2] *= 0.92 ;
0437 dedxcorr[0][1][0] *= 0.90 ;
0438 dedxcorr[0][1][1] *= 0.92 ;
0439 dedxcorr[0][1][2] *= 0.94 ;
0440 dedxcorr[0][2][0] *= 0.80 ;
0441 dedxcorr[0][2][1] *= 0.81 ;
0442 dedxcorr[0][2][2] *= 0.88 ;
0443 dedxcorr[0][3][0] *= 0.76 ;
0444 dedxcorr[0][3][1] *= 0.76 ;
0445 dedxcorr[0][3][2] *= 0.77 ;
0446 dedxcorr[0][4][0] *= 0.87 ;
0447 dedxcorr[0][4][1] *= 0.88 ;
0448 dedxcorr[0][4][2] *= 0.85 ;
0449 dedxcorr[0][5][0] *= 0.86 ;
0450 dedxcorr[0][5][1] *= 0.87 ;
0451 dedxcorr[0][5][2] *= 0.92 ;
0452 dedxcorr[0][6][0] *= 0.89 ;
0453 dedxcorr[0][6][1] *= 0.91 ;
0454 dedxcorr[0][6][2] *= 0.91 ;
0455 dedxcorr[0][7][0] *= 0.90 ;
0456 dedxcorr[0][7][1] *= 0.90 ;
0457 dedxcorr[0][7][2] *= 0.94 ;
0458 dedxcorr[0][8][0] *= 0.87 ;
0459 dedxcorr[0][8][1] *= 0.90 ;
0460 dedxcorr[0][8][2] *= 0.91 ;
0461 dedxcorr[0][9][0] *= 0.92 ;
0462 dedxcorr[0][9][1] *= 0.92 ;
0463 dedxcorr[0][9][2] *= 0.94 ;
0464 dedxcorr[0][10][0] *= 0.91 ;
0465 dedxcorr[0][10][1] *= 0.92 ;
0466 dedxcorr[0][10][2] *= 0.93 ;
0467 dedxcorr[0][11][0] *= 0.87 ;
0468 dedxcorr[0][11][1] *= 0.88 ;
0469 dedxcorr[0][11][2] *= 0.91 ;
0470 dedxcorr[1][0][0] *= 0.91 ;
0471 dedxcorr[1][0][1] *= 0.92 ;
0472 dedxcorr[1][0][2] *= 0.92 ;
0473 dedxcorr[1][1][0] *= 0.92 ;
0474 dedxcorr[1][1][1] *= 0.89 ;
0475 dedxcorr[1][1][2] *= 0.91 ;
0476 dedxcorr[1][2][0] *= 0.84 ;
0477 dedxcorr[1][2][1] *= 0.87 ;
0478 dedxcorr[1][2][2] *= 1.05 ;
0479 dedxcorr[1][3][0] *= 0.83 ;
0480 dedxcorr[1][3][1] *= 0.89 ;
0481 dedxcorr[1][3][2] *= 0.90 ;
0482 dedxcorr[1][4][0] *= 0.90 ;
0483 dedxcorr[1][4][1] *= 0.91 ;
0484 dedxcorr[1][4][2] *= 0.91 ;
0485 dedxcorr[1][5][0] *= 0.92 ;
0486 dedxcorr[1][5][1] *= 0.93 ;
0487 dedxcorr[1][5][2] *= 0.93 ;
0488 dedxcorr[1][6][0] *= 0.87 ;
0489 dedxcorr[1][6][1] *= 0.85 ;
0490 dedxcorr[1][6][2] *= 0.99 ;
0491 dedxcorr[1][7][0] *= 0.92 ;
0492 dedxcorr[1][7][1] *= 0.89 ;
0493 dedxcorr[1][7][2] *= 0.89 ;
0494 dedxcorr[1][8][0] *= 0.88 ;
0495 dedxcorr[1][8][1] *= 0.89 ;
0496 dedxcorr[1][8][2] *= 0.95 ;
0497 dedxcorr[1][9][0] *= 0.89 ;
0498 dedxcorr[1][9][1] *= 0.90 ;
0499 dedxcorr[1][9][2] *= 0.91 ;
0500 dedxcorr[1][10][0] *= 0.86 ;
0501 dedxcorr[1][10][1] *= 0.92 ;
0502 dedxcorr[1][10][2] *= 0.95 ;
0503 dedxcorr[1][11][0] *= 0.85 ;
0504 dedxcorr[1][11][1] *= 0.85 ;
0505 dedxcorr[1][11][2] *= 0.88 ;
0506
0507
0508
0509 return Fun4AllReturnCodes::EVENT_OK;
0510 }
0511
0512 int TrkrNtuplizer::InitRun(PHCompositeNode* topNode)
0513 {
0514
0515 auto *geom =
0516 findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0517 if (!geom)
0518 {
0519 std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0520 return Fun4AllReturnCodes::ABORTRUN;
0521 }
0522 AdcClockPeriod = geom->GetFirstLayerCellGeom()->get_zstep();
0523 return Fun4AllReturnCodes::EVENT_OK;
0524 }
0525
0526 int TrkrNtuplizer::process_event(PHCompositeNode* topNode)
0527 {
0528 if ((Verbosity() > 1) && (_ievent % 100 == 0))
0529 {
0530 std::cout << "TrkrNtuplizer::process_event - Event = " << _ievent << std::endl;
0531 }
0532
0533 recoConsts* rc = recoConsts::instance();
0534 if (rc->FlagExist("RANDOMSEED"))
0535 {
0536 _iseed = rc->get_IntFlag("RANDOMSEED");
0537 m_fSeed = _iseed;
0538 }
0539 else
0540 {
0541 _iseed = 0;
0542 m_fSeed = std::numeric_limits<float>::quiet_NaN();
0543 }
0544 if (_trackmap == nullptr)
0545 {
0546 _trackmap = findNode::getClass<SvtxTrackMap>(topNode, _trackmapname);
0547 }
0548
0549 _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "CORRECTED_TRKR_CLUSTER");
0550 if (!_cluster_map)
0551 {
0552 _cluster_map = findNode::getClass<TrkrClusterContainer>(topNode, "TRKR_CLUSTER");
0553 }
0554
0555 _tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0556 if (!_tgeometry)
0557 {
0558 std::cout << "No Acts geometry on node tree. Can't continue."
0559 << std::endl;
0560 return Fun4AllReturnCodes::ABORTEVENT;
0561 }
0562
0563 _geom_container = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0564 {
0565 if (!_geom_container)
0566 {
0567 std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0568 return Fun4AllReturnCodes::ABORTEVENT;
0569 }
0570 }
0571
0572 if (Verbosity() > 1)
0573 {
0574 std::cout << "TrkrNtuplizer::process_event - Seed = " << _iseed << std::endl;
0575 }
0576
0577
0578
0579
0580 printInputInfo(topNode);
0581
0582
0583
0584
0585
0586 fillOutputNtuples(topNode);
0587
0588
0589
0590
0591
0592
0593
0594 ++_ievent;
0595 return Fun4AllReturnCodes::EVENT_OK;
0596 }
0597
0598 int TrkrNtuplizer::End(PHCompositeNode* )
0599 {
0600 _tfile->cd();
0601
0602 if (_ntp_info)
0603 {
0604 _ntp_info->Write();
0605 }
0606 if (_ntp_vertex)
0607 {
0608 _ntp_vertex->Write();
0609 }
0610 if (_ntp_hit)
0611 {
0612 _ntp_hit->Write();
0613 }
0614 if (_ntp_cluster)
0615 {
0616 _ntp_cluster->Write();
0617 }
0618 if (_ntp_clus_trk)
0619 {
0620 _ntp_clus_trk->Write();
0621 }
0622 if (_ntp_track)
0623 {
0624 _ntp_track->Write();
0625 }
0626 if (_ntp_tpcseed)
0627 {
0628 _ntp_tpcseed->Write();
0629 }
0630 if (_ntp_siseed)
0631 {
0632 _ntp_siseed->Write();
0633 }
0634
0635 _tfile->Close();
0636
0637 delete _tfile;
0638
0639 if (Verbosity() > 1)
0640 {
0641 std::cout << "========================= TrkrNtuplizer::End() ============================" << std::endl;
0642 std::cout << " " << _ievent << " events of output written to: " << _filename << std::endl;
0643 std::cout << "===========================================================================" << std::endl;
0644 }
0645
0646 return Fun4AllReturnCodes::EVENT_OK;
0647 }
0648
0649 void TrkrNtuplizer::printInputInfo(PHCompositeNode* topNode)
0650 {
0651 if (Verbosity() > 1)
0652 {
0653 std::cout << "TrkrNtuplizer::printInputInfo() entered" << std::endl;
0654 }
0655
0656 if (Verbosity() > 3)
0657 {
0658
0659 std::cout << std::endl;
0660 std::cout << PHWHERE << " INPUT FOR EVENT " << _ievent << std::endl;
0661
0662 std::cout << std::endl;
0663
0664 std::cout << "---SVTXCLUSTERS-------------" << std::endl;
0665
0666 if (_cluster_map != nullptr)
0667 {
0668 unsigned int icluster = 0;
0669 for (const auto& hitsetkey : _cluster_map->getHitSetKeys())
0670 {
0671 auto range = _cluster_map->getClusters(hitsetkey);
0672 for (auto iter = range.first; iter != range.second; ++iter)
0673 {
0674 TrkrDefs::cluskey cluster_key = iter->first;
0675 std::cout << icluster << " with key " << cluster_key << " of " << _cluster_map->size();
0676 std::cout << ": SvtxCluster: " << std::endl;
0677 iter->second->identify();
0678 ++icluster;
0679 }
0680 }
0681 }
0682
0683 std::cout << "---SVXTRACKS-------------" << std::endl;
0684
0685 if (_trackmap)
0686 {
0687 unsigned int itrack = 0;
0688 for (auto& iter : *_trackmap)
0689 {
0690 std::cout << itrack << " of " << _trackmap->size();
0691 SvtxTrack* track = iter.second;
0692 std::cout << " : SvtxTrack:" << std::endl;
0693 track->identify();
0694 std::cout << std::endl;
0695 ++itrack;
0696 }
0697 }
0698
0699 std::cout << "---SVXVERTEXES-------------" << std::endl;
0700 SvtxVertexMap* vertexmap = nullptr;
0701 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs");
0702
0703 if (vertexmap)
0704 {
0705 unsigned int ivertex = 0;
0706 for (SvtxVertexMap::Iter iter = vertexmap->begin();
0707 iter != vertexmap->end();
0708 ++iter)
0709 {
0710 std::cout << ivertex << " of " << vertexmap->size();
0711 SvtxVertex* vertex = iter->second;
0712 std::cout << " : SvtxVertex:" << std::endl;
0713 vertex->identify();
0714 std::cout << std::endl;
0715 }
0716 }
0717 }
0718
0719 return;
0720 }
0721
0722 void TrkrNtuplizer::printOutputInfo(PHCompositeNode* topNode)
0723 {
0724 if (Verbosity() > 1)
0725 {
0726 std::cout << "TrkrNtuplizer::printOutputInfo() entered" << std::endl;
0727 }
0728
0729
0730
0731
0732
0733 if (Verbosity() > 100)
0734 {
0735
0736 std::cout << std::endl;
0737 std::cout << PHWHERE << " NEW OUTPUT FOR EVENT " << _ievent << std::endl;
0738 std::cout << std::endl;
0739
0740 float vx = std::numeric_limits<float>::quiet_NaN();
0741 float vy = std::numeric_limits<float>::quiet_NaN();
0742 float vz = std::numeric_limits<float>::quiet_NaN();
0743
0744 SvtxVertexMap* vertexmap = nullptr;
0745
0746 vertexmap = findNode::getClass<SvtxVertexMap>(topNode, "SvtxVertexMapActs");
0747
0748 if (vertexmap)
0749 {
0750 if (!vertexmap->empty())
0751 {
0752 SvtxVertex* vertex = (vertexmap->begin()->second);
0753
0754 vx = vertex->get_x();
0755 vy = vertex->get_y();
0756 vz = vertex->get_z();
0757 }
0758 }
0759
0760 std::cout << "===Vertex Reconstruction=======================" << std::endl;
0761 std::cout << "vreco = (" << vx << "," << vy << "," << vz << ")" << std::endl;
0762 std::cout << std::endl;
0763
0764 std::cout << "===Tracking Summary============================" << std::endl;
0765
0766 TrkrHitSetContainer* hitsetmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0767
0768 unsigned int nclusters[100] = {0};
0769 unsigned int nhits[100] = {0};
0770
0771 _tgeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
0772
0773 if (!_tgeometry)
0774 {
0775 std::cout << PHWHERE << "No Acts geometry on node tree. Can't continue."
0776 << std::endl;
0777 }
0778
0779 if (hitsetmap)
0780 {
0781 TrkrHitSetContainer::ConstRange all_hitsets = hitsetmap->getHitSets();
0782 for (TrkrHitSetContainer::ConstIterator hitsetiter = all_hitsets.first;
0783 hitsetiter != all_hitsets.second;
0784 ++hitsetiter)
0785 {
0786
0787 unsigned int layer = TrkrDefs::getLayer(hitsetiter->first);
0788
0789
0790 TrkrHitSet::ConstRange hitrangei = hitsetiter->second->getHits();
0791 for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0792 hitr != hitrangei.second;
0793 ++hitr)
0794 {
0795 ++nhits[layer];
0796 }
0797 auto range = _cluster_map->getClusters(hitsetiter->first);
0798 for (auto clusIter = range.first; clusIter != range.second; ++clusIter)
0799 {
0800 const auto cluskey = clusIter->first;
0801 nclusters[TrkrDefs::getLayer(cluskey)]++;
0802 }
0803 }
0804 }
0805
0806 for (unsigned int ilayer = 0; ilayer < _nlayers_maps + _nlayers_intt + _nlayers_tpc; ++ilayer)
0807 {
0808 PHG4TpcCylinderGeom* GeoLayer = _geom_container->GetLayerCellGeom(ilayer);
0809
0810 std::cout << "layer " << ilayer
0811 << " => nHits = " << nhits[ilayer]
0812 << " => nClusters = " << nclusters[ilayer] << std::endl;
0813 if (ilayer >= _nlayers_maps + _nlayers_intt && ilayer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
0814 {
0815 std::cout << "layer " << ilayer
0816 << " => nphi = " << GeoLayer->get_phibins()
0817 << " => nz = " << GeoLayer->get_zbins()
0818 << " => ntot = " << GeoLayer->get_phibins() * GeoLayer->get_zbins()
0819 << std::endl;
0820 }
0821 }
0822
0823 std::cout << " => nTracks = ";
0824 if (_trackmap)
0825 {
0826 std::cout << _trackmap->size() << std::endl;
0827 }
0828 else
0829 {
0830 std::cout << 0 << std::endl;
0831 }
0832
0833 std::cout << std::endl;
0834
0835 }
0836
0837 return;
0838 }
0839
0840 void TrkrNtuplizer::fillOutputNtuples(PHCompositeNode* topNode)
0841 {
0842 if (Verbosity() > 1)
0843 {
0844 std::cout << "TrkrNtuplizer::fillOutputNtuples() entered" << std::endl;
0845 }
0846
0847 float fx_event[n_event::evsize]{(float) _ievent, (float) _iseed, (float) m_runnumber,(float) m_segment, (float) m_job};
0848 float fx_info[((int) (n_info::infosize))] = {0};
0849
0850 float nhit[100];
0851 for (float& i : nhit)
0852 {
0853 i = 0;
0854 }
0855
0856 if (_ntp_info){
0857 TrkrHitSetContainer* hitmap_in = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
0858 if (hitmap_in)
0859 {
0860 TrkrHitSetContainer::ConstRange all_hitsets = hitmap_in->getHitSets();
0861 for (TrkrHitSetContainer::ConstIterator hitsetiter = all_hitsets.first;
0862 hitsetiter != all_hitsets.second;
0863 ++hitsetiter)
0864 {
0865
0866 unsigned int layer = TrkrDefs::getLayer(hitsetiter->first);
0867
0868 TrkrHitSet::ConstRange hitrangei = hitsetiter->second->getHits();
0869 for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
0870 hitr != hitrangei.second;
0871 ++hitr)
0872 {
0873 nhit[layer]++;
0874
0875 if (layer < _nlayers_maps)
0876 {
0877 fx_info[n_info::infonhitmvtx]++;
0878 }
0879 if (layer >= _nlayers_maps && layer < _nlayers_maps + _nlayers_intt)
0880 {
0881 fx_info[n_info::infonhitintt]++;
0882 }
0883 if ((float) layer >= _nlayers_maps + _nlayers_intt &&
0884 (float) layer < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
0885 {
0886 fx_info[n_info::infonhittpcall]++;
0887 }
0888 if ((float) layer >= _nlayers_maps + _nlayers_intt + _nlayers_tpc)
0889 {
0890 fx_info[n_info::infonhittpot]++;
0891 }
0892 if ((float) layer == _nlayers_maps + _nlayers_intt)
0893 {
0894 fx_info[n_info::infonhittpcin]++;
0895 }
0896 if ((float) layer == _nlayers_maps + _nlayers_intt + _nlayers_tpc - 1)
0897 {
0898 fx_info[n_info::infonhittpcout]++;
0899 }
0900
0901 if ((float) layer == _nlayers_maps + _nlayers_intt + _nlayers_tpc / 2 - 1)
0902 {
0903 fx_info[n_info::infonhittpcmid]++;
0904 }
0905 }
0906 }
0907 }
0908
0909 _geom_container = findNode::getClass<PHG4TpcCylinderGeomContainer>(topNode, "CYLINDERCELLGEOM_SVTX");
0910 if (!_geom_container)
0911 {
0912 std::cout << PHWHERE << "ERROR: Can't find node CYLINDERCELLGEOM_SVTX" << std::endl;
0913 return;
0914 }
0915 PHG4TpcCylinderGeom* GeoLayer;
0916 int layer = _nlayers_maps + _nlayers_intt;
0917 GeoLayer = _geom_container->GetLayerCellGeom(layer);
0918 int nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0919 float nhits = nhit[layer];
0920 fx_info[n_info::infonocc11] = nhits / nbins;
0921
0922 layer = _nlayers_maps + _nlayers_intt + 15;
0923 GeoLayer = _geom_container->GetLayerCellGeom(layer);
0924 nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0925 nhits = nhit[layer];
0926 fx_info[n_info::infonocc116] = nhits / nbins;
0927
0928 layer = _nlayers_maps + _nlayers_intt + 16;
0929 GeoLayer = _geom_container->GetLayerCellGeom(layer);
0930 nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0931 nhits = nhit[layer];
0932 fx_info[n_info::infonocc21] = nhits / nbins;
0933 layer = _nlayers_maps + _nlayers_intt + 31;
0934 GeoLayer = _geom_container->GetLayerCellGeom(layer);
0935 nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0936 nhits = nhit[layer];
0937 fx_info[n_info::infonocc216] = nhits / nbins;
0938 layer = _nlayers_maps + _nlayers_intt + 32;
0939 GeoLayer = _geom_container->GetLayerCellGeom(layer);
0940 nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0941 nhits = nhit[layer];
0942 fx_info[n_info::infonocc31] = nhits / nbins;
0943 layer = _nlayers_maps + _nlayers_intt + 47;
0944 GeoLayer = _geom_container->GetLayerCellGeom(layer);
0945 nbins = GeoLayer->get_phibins() * GeoLayer->get_zbins();
0946 nhits = nhit[layer];
0947 fx_info[n_info::infonocc316] = nhits / nbins;
0948
0949 if (Verbosity() > 1)
0950 {
0951 std::cout << " occ11 = " << fx_info[n_info::infonocc11]
0952 << " occ116 = " << fx_info[n_info::infonocc116]
0953 << " occ21 = " << fx_info[n_info::infonocc21]
0954 << " occ216 = " << fx_info[n_info::infonocc216]
0955 << " occ31 = " << fx_info[n_info::infonocc31]
0956 << " occ316 = " << fx_info[n_info::infonocc316]
0957 << std::endl;
0958 }
0959
0960 if (_cluster_map)
0961 {
0962 fx_info[n_info::infonclusall] = _cluster_map->size();
0963
0964 for (const auto& hitsetkey : _cluster_map->getHitSetKeys())
0965 {
0966 auto range = _cluster_map->getClusters(hitsetkey);
0967 for (auto iter_cin = range.first; iter_cin != range.second; ++iter_cin)
0968 {
0969 TrkrDefs::cluskey cluster_key = iter_cin->first;
0970 unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
0971 unsigned int side_local = TpcDefs::getSide(cluster_key);
0972 if (_nlayers_maps > 0)
0973 {
0974 if (layer_local < _nlayers_maps)
0975 {
0976 fx_info[n_info::infonclusmvtx]++;
0977 }
0978 }
0979 if (_nlayers_intt > 0)
0980 {
0981 if ((layer_local >= _nlayers_maps) && (layer_local < (_nlayers_maps + _nlayers_intt)))
0982 {
0983 fx_info[n_info::infonclusintt]++;
0984 }
0985 }
0986 if (_nlayers_tpc > 0)
0987 {
0988 if (layer_local >= (_nlayers_maps + _nlayers_intt) && layer_local < (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
0989 {
0990 fx_info[n_info::infonclustpc]++;
0991 if(side_local==0){
0992 fx_info[n_info::infonclustpcneg]++;
0993 }
0994 if(side_local==1){
0995 fx_info[n_info::infonclustpcpos]++;
0996 }
0997 }
0998 }
0999 if (_nlayers_mms > 0)
1000 {
1001 if (layer_local >= (_nlayers_maps + _nlayers_intt + _nlayers_tpc))
1002 {
1003 fx_info[n_info::infonclustpot]++;
1004 }
1005 }
1006 }
1007 }
1008 }
1009 }
1010
1011
1012
1013 if (_ntp_info)
1014 {
1015 if (Verbosity() > 1)
1016 {
1017 std::cout << "Filling ntp_info " << std::endl;
1018 }
1019
1020 if (_trackmap)
1021 {
1022 fx_info[n_info::infontrk] = (float) _trackmap->size();
1023 }
1024 auto *siseedmap = findNode::getClass<TrackSeedContainer>(topNode, "SiliconTrackSeedContainer");
1025 auto *tpcseedmap = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
1026 if (tpcseedmap){
1027 fx_info[n_info::infonntpcseed] = tpcseedmap->size();
1028 }else{
1029 fx_info[n_info::infonntpcseed] = 0;
1030 }
1031 if (siseedmap){
1032 fx_info[n_info::infonnsiseed] = siseedmap->size();
1033 }else{
1034 fx_info[n_info::infonnsiseed] = 0;
1035 }
1036 if (Verbosity() > 0)
1037 {
1038 std::cout << "EVENTINFO SEED: " << m_fSeed << std::endl;
1039 std::cout << "EVENTINFO NHIT: " << std::setprecision(9) << fx_info[n_info::infonclusall] << std::endl;
1040 std::cout << "EVENTINFO CLUSTPC: " << fx_info[n_info::infonclustpc] << std::endl;
1041 std::cout << "EVENTINFO NTRKREC: " << fx_info[n_info::infontrk] << std::endl;
1042 }
1043
1044 float* info_data = new float[((int) (n_info::infosize)) + n_event::evsize];
1045 std::copy(fx_event, fx_event + n_event::evsize, info_data);
1046 std::copy(fx_info, fx_info + ((int) (n_info::infosize)), info_data + n_event::evsize);
1047 _ntp_info->Fill(info_data);
1048 delete[] info_data;
1049 }
1050
1051
1052
1053
1054 bool doit = true;
1055 if (_ntp_vertex && doit)
1056 {
1057 if (Verbosity() > 1)
1058 {
1059 std::cout << "Filling ntp_vertex " << std::endl;
1060 std::cout << "start vertex time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
1061 _timer->restart();
1062 }
1063 float fx_vertex[n_vertex::vtxsize];
1064 for (float& i : fx_vertex)
1065 {
1066 i = 0;
1067 }
1068
1069
1070
1071
1072
1073 float vx = std::numeric_limits<float>::quiet_NaN();
1074 float vy = std::numeric_limits<float>::quiet_NaN();
1075 float vz = std::numeric_limits<float>::quiet_NaN();
1076 float ntracks = std::numeric_limits<float>::quiet_NaN();
1077 fx_vertex[vtxnvx] = vx;
1078 fx_vertex[vtxnvy] = vy;
1079 fx_vertex[vtxnvz] = vz;
1080 fx_vertex[vtxnntracks] = ntracks;
1081 if (Verbosity() > 1)
1082 {
1083 std::cout << " adding vertex data " << std::endl;
1084 }
1085 float* vertex_data = new float[((int) (n_info::infosize)) + n_event::evsize + n_vertex::vtxsize];
1086 std::copy(fx_event, fx_event + n_event::evsize, vertex_data);
1087 std::copy(fx_vertex, fx_vertex + n_vertex::vtxsize, vertex_data + n_event::evsize);
1088 std::copy(fx_info, fx_info + ((int) (n_info::infosize)), vertex_data + n_event::evsize + n_vertex::vtxsize);
1089 _ntp_vertex->Fill(vertex_data);
1090 delete[] vertex_data;
1091 }
1092 if (Verbosity() > 1)
1093 {
1094 _timer->stop();
1095 std::cout << "vertex time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
1096 }
1097
1098
1099
1100
1101 if (_ntp_hit)
1102 {
1103 float* hit_data = new float[((int) (n_info::infosize)) + n_event::evsize + n_hit::hitsize];
1104 float fx_hit_0[((int) (n_info::infosize)) + n_event::evsize +n_hit::hitsize] = {0};
1105 float fx_hit[((int) (n_hit::hitsize))] = {0};
1106 auto *m_tGeometry = findNode::getClass<ActsGeometry>(topNode, "ActsGeometry");
1107
1108 if (Verbosity() >= 1)
1109 {
1110 std::cout << "Filling ntp_hit " << std::endl;
1111 _timer->restart();
1112 }
1113
1114 TrkrHitSetContainer* hitmap = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1115
1116 if (hitmap)
1117 {
1118 TrkrHitSetContainer::ConstRange all_hitsets = hitmap->getHitSets();
1119 for (TrkrHitSetContainer::ConstIterator iter = all_hitsets.first;
1120 iter != all_hitsets.second;
1121 ++iter)
1122 {
1123 TrkrDefs::hitsetkey hitset_key = iter->first;
1124 TrkrHitSet* hitset = iter->second;
1125
1126
1127 TrkrHitSet::ConstRange hitrangei = hitset->getHits();
1128 for (TrkrHitSet::ConstIterator hitr = hitrangei.first;
1129 hitr != hitrangei.second;
1130 ++hitr)
1131 {
1132 TrkrDefs::hitkey hit_key = hitr->first;
1133 TrkrHit* hit = hitr->second;
1134
1135 fx_hit[n_hit::nhitID] = hit_key;
1136 fx_hit[n_hit::nhite] = hit->getEnergy();
1137 fx_hit[n_hit::nhitadc] = hit->getAdc();
1138 unsigned int layer_local = TrkrDefs::getLayer(hitset_key);
1139 fx_hit[n_hit::nhitlayer] = (float) layer_local;
1140 fx_hit[n_hit::nhitphielem] = -666;
1141 fx_hit[n_hit::nhitzelem] = -666;
1142
1143 if (layer_local >= 3 && layer_local < 7)
1144 {
1145 fx_hit[n_hit::nhitphielem] = InttDefs::getLadderPhiId(hitset_key);
1146 fx_hit[n_hit::nhitzelem] = InttDefs::getLadderZId(hitset_key);
1147 }
1148 if (layer_local >= 7 && layer_local < 55)
1149 {
1150 fx_hit[n_hit::nhitphielem] = TpcDefs::getSectorId(hitset_key);
1151 fx_hit[n_hit::nhitzelem] = TpcDefs::getSide(hitset_key);
1152 }
1153
1154
1155
1156
1157
1158
1159
1160
1161
1162
1163
1164 fx_hit[n_hit::nhitphielem] = TpcDefs::getSectorId(hitset_key);
1165 fx_hit[n_hit::nhitzelem] = TpcDefs::getSide(hitset_key);
1166 fx_hit[n_hit::nhitcellID] = 0;
1167 fx_hit[n_hit::nhitecell] = hit->getAdc();
1168 fx_hit[n_hit::nhitphibin] = std::numeric_limits<float>::quiet_NaN();
1169 fx_hit[n_hit::nhittbin] = std::numeric_limits<float>::quiet_NaN();
1170 fx_hit[n_hit::nhitphi] = std::numeric_limits<float>::quiet_NaN();
1171 fx_hit[n_hit::nhitr] = std::numeric_limits<float>::quiet_NaN();
1172 fx_hit[n_hit::nhitx] = std::numeric_limits<float>::quiet_NaN();
1173 fx_hit[n_hit::nhity] = std::numeric_limits<float>::quiet_NaN();
1174 fx_hit[n_hit::nhitz] = std::numeric_limits<float>::quiet_NaN();
1175
1176 if (layer_local >= _nlayers_maps + _nlayers_intt && layer_local < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
1177 {
1178 PHG4TpcCylinderGeom* GeoLayer_local = _geom_container->GetLayerCellGeom(layer_local);
1179 double radius = GeoLayer_local->get_radius();
1180 fx_hit[n_hit::nhitphibin] = (float) TpcDefs::getPad(hit_key);
1181 fx_hit[n_hit::nhittbin] = (float) TpcDefs::getTBin(hit_key);
1182 fx_hit[n_hit::nhitphi] = GeoLayer_local->get_phicenter(fx_hit[n_hit::nhitphibin], fx_hit[n_hit::nhitzelem]);
1183 float phi = GeoLayer_local->get_phicenter(TpcDefs::getPad(hit_key), fx_hit[n_hit::nhitzelem]);
1184 float clockperiod = GeoLayer_local->get_zstep();
1185 auto glob = m_tGeometry->getGlobalPositionTpc(hitset_key, hit_key, phi, radius, clockperiod);
1186
1187 fx_hit[n_hit::nhitz] = glob.z();
1188 fx_hit[n_hit::nhitr] = radius;
1189 fx_hit[n_hit::nhitx] = glob.x();
1190 fx_hit[n_hit::nhity] = glob.y();
1191 }
1192
1193 std::copy(fx_event, fx_event + n_event::evsize, hit_data);
1194 std::copy(fx_hit, fx_hit + n_hit::hitsize, hit_data + n_event::evsize);
1195 std::copy(fx_info, fx_info + ((int) (n_info::infosize)), hit_data + n_event::evsize + n_hit::hitsize);
1196 _ntp_hit->Fill(hit_data);
1197
1198 std::copy(fx_hit_0, fx_hit_0 + ((int) (n_info::infosize)) + n_event::evsize + n_hit::hitsize, hit_data);
1199 }
1200 }
1201 delete[] hit_data;
1202 }
1203 if (Verbosity() >= 1)
1204 {
1205 _timer->stop();
1206 std::cout << "hit time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
1207 }
1208 }
1209
1210
1211
1212
1213
1214 if (Verbosity() >= 1)
1215 {
1216 std::cout << "check for ntp_cluster" << std::endl;
1217 _timer->restart();
1218 }
1219
1220 if (_ntp_cluster)
1221 {
1222 if (Verbosity() > 1)
1223 {
1224 std::cout << "Filling ntp_cluster (all of them) " << std::endl;
1225 }
1226
1227 TrkrHitSetContainer* hitsets = findNode::getClass<TrkrHitSetContainer>(topNode, "TRKR_HITSET");
1228
1229 ClusterErrorPara ClusErrPara;
1230
1231 if (Verbosity() > 1)
1232 {
1233 if (_cluster_map != nullptr)
1234 {
1235 std::cout << "got clustermap" << std::endl;
1236 }
1237 else
1238 {
1239 std::cout << "no clustermap" << std::endl;
1240 }
1241
1242 if (hitsets != nullptr)
1243 {
1244 std::cout << "got hitsets" << std::endl;
1245 }
1246 else
1247 {
1248 std::cout << "no hitsets" << std::endl;
1249 }
1250 }
1251
1252 if (_cluster_map && hitsets){
1253 float* cluster_data = new float[((int) (n_info::infosize)) + n_event::evsize + n_cluster::clusize];
1254 Float_t fx_cluster_0[((int) (n_info::infosize)) + n_event::evsize + n_cluster::clusize] = {0};
1255 for (const auto& hitsetkey : _cluster_map->getHitSetKeys())
1256 {
1257 auto range = _cluster_map->getClusters(hitsetkey);
1258 for (auto iter = range.first; iter != range.second; ++iter)
1259 {
1260 TrkrDefs::cluskey cluster_key = iter->first;
1261 Float_t fx_cluster[n_cluster::clusize];
1262 FillCluster(&fx_cluster[0], cluster_key);
1263
1264 std::copy(fx_event, fx_event + n_event::evsize, cluster_data);
1265 std::copy(fx_cluster, fx_cluster + n_cluster::clusize, cluster_data + n_event::evsize);
1266 std::copy(fx_info, fx_info + ((int) (n_info::infosize)), cluster_data + n_event::evsize + n_cluster::clusize);
1267 _ntp_cluster->Fill(cluster_data);
1268 std::copy(fx_cluster_0, fx_cluster_0 + ((int) (n_info::infosize)) + n_event::evsize + n_cluster::clusize, cluster_data);
1269 }
1270 }
1271 delete[] cluster_data;
1272 }
1273 }
1274
1275 if (Verbosity() >= 1)
1276 {
1277 _timer->stop();
1278 std::cout << "cluster time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
1279 }
1280 if (_ntp_clus_trk)
1281 {
1282
1283
1284 if (Verbosity() > 1)
1285 {
1286 std::cout << "Filling ntp_clus_trk " << std::endl;
1287 if (_cluster_map != nullptr)
1288 {
1289 std::cout << "got clustermap" << std::endl;
1290 }
1291 else
1292 {
1293 std::cout << "no clustermap" << std::endl;
1294 }
1295 }
1296
1297 TrackSeedContainer* _tpc_seeds = findNode::getClass<TrackSeedContainer>(topNode, "TpcTrackSeedContainer");
1298 if (!_tpc_seeds)
1299 {
1300 std::cout << PHWHERE << " ERROR: Can't find "
1301 << "TpcTrackSeedContainer" << std::endl;
1302 return;
1303 }
1304
1305 if (_trackmap)
1306 {
1307 int trackID = 0;
1308 float* clus_trk_data = new float[((int) (n_info::infosize)) + n_cluster::clusize + n_residual::ressize + n_seed::seedsize + n_event::evsize];
1309 float fx_clus_trk_0[((int) (n_info::infosize)) + n_cluster::clusize + n_residual::ressize + n_seed::seedsize + n_event::evsize] = {0};
1310 for (auto& iter : *_trackmap)
1311 {
1312 trackID++;
1313 SvtxTrack* track = iter.second;
1314 TrackSeed* tpcseed = track->get_tpc_seed();
1315 TrackSeed* siseed = track->get_silicon_seed();
1316 std::vector<Acts::Vector3> clusterPositions;
1317 std::vector<TrkrDefs::cluskey> clusterKeys;
1318
1319 TrackFitUtils::position_vector_t xypoints;
1320 TrackFitUtils::position_vector_t rzpoints;
1321
1322 clusterKeys.insert(clusterKeys.end(), tpcseed->begin_cluster_keys(),
1323 tpcseed->end_cluster_keys());
1324
1325 TrackFitUtils::getTrackletClusters(_tgeometry, _cluster_map,
1326 clusterPositions, clusterKeys);
1327 for (auto& pos : clusterPositions)
1328 {
1329 float clusr = sqrt((pos.x() * pos.x()) + (pos.y() * pos.y()));
1330 if (pos.y() < 0)
1331 {
1332 clusr *= -1;
1333 }
1334
1335
1336 if (std::fabs(clusr) > 80 || std::fabs(clusr) < 30)
1337 {
1338 continue;
1339 }
1340 rzpoints.emplace_back(pos.z(), clusr);
1341 xypoints.emplace_back(pos.x(), pos.y());
1342 }
1343
1344 std::vector<float> fitparams = TrackFitUtils::fitClusters(clusterPositions, clusterKeys);
1345 if (fitparams.empty())
1346 {
1347 std::cout << "fit failed bailing...." << std::endl;
1348 continue;
1349 }
1350
1351 float charge = std::numeric_limits<float>::quiet_NaN();
1352 if (tpcseed->get_qOverR() > 0)
1353 {
1354 charge = 1;
1355 }
1356 else
1357 {
1358 charge = -1;
1359 }
1360
1361
1362 float tpt = tpcseed->get_pt();
1363 float tptot = tpcseed->get_p();
1364 float teta = tpcseed->get_eta();
1365 float tphi = tpcseed->get_phi();
1366 auto xyparams = TrackFitUtils::line_fit(xypoints);
1367 auto rzparams = TrackFitUtils::line_fit(rzpoints);
1368 float xyint = std::get<1>(xyparams);
1369 float xyslope = std::get<0>(xyparams);
1370 float rzint = std::get<1>(rzparams);
1371 float rzslope = std::get<0>(rzparams);
1372 float R0 = abs(-1 * xyint) / std::sqrt((xyslope * xyslope) + 1);
1373 float tX0 = tpcseed->get_X0();
1374 float tY0 = tpcseed->get_Y0();
1375 float tZ0 = tpcseed->get_Z0();
1376
1377 float nhits_local = clusterPositions.size();
1378 if (Verbosity() > 1)
1379 {
1380 std::cout << " tpc: " << tpcseed->size_cluster_keys() << std::endl;
1381 if (siseed)
1382 {
1383 std::cout << " si " << siseed->size_cluster_keys() << std::endl;
1384 }
1385 std::cout << "done seedsize" << std::endl;
1386 }
1387
1388
1389
1390 float dedx = calc_dedx(tpcseed);
1391 float pidedx = 0;
1392 float kdedx = 0;
1393 float prdedx = 0;
1394 if(charge>0){
1395 pidedx = f_pion_plus->Eval(tptot);
1396 kdedx = f_kaon_plus->Eval(tptot);
1397 prdedx = f_proton_plus->Eval(tptot);
1398 }else{
1399 pidedx = f_pion_minus->Eval(tptot);
1400 kdedx = f_kaon_minus->Eval(tptot);
1401 prdedx = f_proton_plus->Eval(tptot);
1402 }
1403 float n1pix = get_n1pix(tpcseed);
1404 float fx_seed[n_seed::seedsize] = {(float) trackID, 0, tpt, tptot, teta, tphi, xyint, rzint, xyslope, rzslope, tX0, tY0, tZ0, R0, charge, dedx, pidedx, kdedx, prdedx, n1pix, nhits_local};
1405
1406 if (_ntp_tpcseed)
1407 {
1408 float* tpcseed_data = new float[((int) (n_info::infosize)) + n_cluster::clusize + n_residual::ressize + n_seed::seedsize + n_event::evsize];
1409 std::copy(fx_event, fx_event + n_event::evsize, tpcseed_data);
1410 std::copy(fx_seed, fx_seed + n_seed::seedsize, tpcseed_data + n_event::evsize);
1411 std::copy(fx_info, fx_info + ((int) (n_info::infosize)), tpcseed_data + n_event::evsize + n_seed::seedsize);
1412 _ntp_tpcseed->Fill(tpcseed_data);
1413 delete[] tpcseed_data;
1414 }
1415 for (unsigned int i = 0; i < clusterPositions.size(); i++)
1416 {
1417 TrkrDefs::cluskey cluster_key = clusterKeys.at(i);
1418 unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
1419 Acts::Vector3 position = clusterPositions[i];
1420 Acts::Vector3 pca = TrackFitUtils::get_helix_pca(fitparams, position);
1421
1422 float cluster_phi = atan2(position(1), position(0));
1423 float pca_phi = atan2(pca(1), pca(0));
1424 float dphi = cluster_phi - pca_phi;
1425 if (dphi > M_PI)
1426 {
1427 dphi = 2 * M_PI - dphi;
1428 }
1429 if (dphi < -M_PI)
1430 {
1431 dphi = 2 * M_PI + dphi;
1432 }
1433
1434 float dz = position(2) - pca(2);
1435
1436 float resr = sqrt((position(0)*position(0))+(position(1)*position(1))+(position(2)*position(2)));
1437 float seedR = std::abs(1.0 / tpcseed->get_qOverR());
1438 float alpha = (resr * resr) / (2 * resr * seedR);
1439 float beta = std::abs(std::atan(tpcseed->get_slope()));
1440 float fx_cluster[n_cluster::clusize];
1441
1442 if(layer_local>=7&&layer_local<55){
1443 PHG4TpcCylinderGeom* GeoLayer_local = _geom_container->GetLayerCellGeom(layer_local);
1444 float thick = GeoLayer_local->get_thickness();
1445
1446 float alphacorr = std::cos(alpha);
1447 if(alphacorr<0||alphacorr>4){
1448 alphacorr=4;
1449 }
1450 float betacorr = std::cos(beta);
1451 if(betacorr<0||betacorr>4){
1452 betacorr=4;
1453 }
1454 fx_cluster[ncluthick] = thick;
1455 fx_cluster[ncluafac] = alphacorr;
1456 fx_cluster[nclubfac] =betacorr;
1457 }
1458 else{
1459 fx_cluster[ncluthick] = 0;
1460 fx_cluster[ncluafac] = 0;
1461 fx_cluster[nclubfac] = 0;
1462 }
1463 float fx_res[n_residual::ressize] = {alpha,beta,dphi,dphi,dz};
1464
1465 FillCluster(&fx_cluster[0], cluster_key);
1466 std::copy(fx_event, fx_event + n_event::evsize, clus_trk_data);
1467 std::copy(fx_cluster, fx_cluster + n_cluster::clusize, clus_trk_data + n_event::evsize);
1468 std::copy(fx_res, fx_res + n_residual::ressize, clus_trk_data + n_event::evsize + n_cluster::clusize);
1469 std::copy(fx_seed, fx_seed + n_seed::seedsize, clus_trk_data + n_event::evsize + n_cluster::clusize + n_residual::ressize);
1470 std::copy(fx_info, fx_info + ((int) (n_info::infosize)), clus_trk_data + n_event::evsize + n_cluster::clusize + n_residual::ressize + n_seed::seedsize);
1471 _ntp_clus_trk->Fill(clus_trk_data);
1472 std::copy(fx_clus_trk_0, fx_clus_trk_0 + ((int) (n_info::infosize)) + n_cluster::clusize + n_residual::ressize + n_seed::seedsize + n_event::evsize,clus_trk_data);
1473
1474 }
1475 }
1476 delete[] clus_trk_data;
1477 }
1478 }
1479
1480
1481
1482
1483
1484 if (_ntp_track)
1485 {
1486 if (Verbosity() >= 1)
1487 {
1488 std::cout << "Filling ntp_track " << std::endl;
1489 _timer->restart();
1490 }
1491 GlobalVertexMap* vertexmap = findNode::getClass<GlobalVertexMap>(topNode, "GlobalVertexMap");
1492 if (_trackmap)
1493 {
1494 for (auto& iter : *_trackmap)
1495 {
1496 SvtxTrack* track = iter.second;
1497 float fx_track[n_track::trksize];
1498 FillTrack(&fx_track[0], track, vertexmap);
1499 float* track_data = new float[((int) (n_info::infosize)) + n_track::trksize + n_event::evsize];
1500 std::copy(fx_event, fx_event + n_event::evsize, track_data);
1501 std::copy(fx_track, fx_track + n_track::trksize, track_data + n_event::evsize);
1502 std::copy(fx_info, fx_info + ((int) (n_info::infosize)), track_data + n_event::evsize + n_track::trksize);
1503 _ntp_track->Fill(track_data);
1504 delete[] track_data;
1505 }
1506 }
1507 if (Verbosity() > 1)
1508 {
1509 _timer->stop();
1510 std::cout << "track time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
1511 }
1512 }
1513
1514 if (_ntp_siseed)
1515 {
1516 if (Verbosity() > 1)
1517 {
1518 std::cout << "Filling ntp_tpcseed " << std::endl;
1519 _timer->restart();
1520 }
1521
1522 if (Verbosity() > 1)
1523 {
1524 _timer->stop();
1525 std::cout << "tpcseed time: " << _timer->get_accumulated_time() / 1000. << " sec" << std::endl;
1526 }
1527 }
1528 return;
1529 }
1530
1531 void TrkrNtuplizer::FillTrack(float fX[50], SvtxTrack* track, GlobalVertexMap* vertexmap)
1532 {
1533 float trackID = track->get_id();
1534 TrackSeed* tpcseed = track->get_tpc_seed();
1535 TrackSeed* silseed = track->get_silicon_seed();
1536 short int crossing_int = track->get_crossing();
1537 fX[n_track::ntrktrackID] = trackID;
1538 fX[n_track::ntrkcrossing] = std::numeric_limits<float>::quiet_NaN();
1539 if (crossing_int == SHRT_MAX)
1540 {
1541 fX[n_track::ntrkcrossing] = std::numeric_limits<float>::quiet_NaN();
1542 }
1543 else
1544 {
1545 fX[n_track::ntrkcrossing] = (float) crossing_int;
1546 }
1547 fX[n_track::ntrkcharge] = track->get_charge();
1548 fX[n_track::ntrkquality] = track->get_quality();
1549 fX[n_track::ntrkchisq] = track->get_chisq();
1550 fX[n_track::ntrkndf] = track->get_ndf();
1551 fX[n_track::ntrknhits] = 0;
1552 if (tpcseed)
1553 {
1554 fX[n_track::ntrknhits] += tpcseed->size_cluster_keys();
1555 }
1556 if (silseed)
1557 {
1558 fX[n_track::ntrknhits] += silseed->size_cluster_keys();
1559 }
1560
1561 fX[n_track::ntrknmaps] = 0;
1562 fX[n_track::ntrknintt] = 0;
1563 fX[n_track::ntrknmms] = 0;
1564 fX[n_track::ntrkntpc] = 0;
1565 fX[n_track::ntrkntpc1] = 0;
1566 fX[n_track::ntrkntpc11] = 0;
1567 fX[n_track::ntrkntpc2] = 0;
1568 fX[n_track::ntrkntpc3] = 0;
1569 fX[n_track::ntrkndedx] = 0;
1570 fX[n_track::ntrknpidedx] = 0;
1571 fX[n_track::ntrknkdedx] = 0;
1572 fX[n_track::ntrknprdedx] = 0;
1573 if (tpcseed)
1574 {
1575 fX[n_track::ntrkndedx] = calc_dedx(tpcseed);
1576 float trptot = track->get_p();
1577 if(track->get_charge()>0){
1578 fX[n_track::ntrknpidedx] = f_pion_plus->Eval(trptot);
1579 fX[n_track::ntrknkdedx] = f_kaon_plus->Eval(trptot);
1580 fX[n_track::ntrknprdedx] = f_proton_plus->Eval(trptot);
1581 }else{
1582 fX[n_track::ntrknpidedx] = f_pion_minus->Eval(trptot);
1583 fX[n_track::ntrknkdedx] = f_kaon_minus->Eval(trptot);
1584 fX[n_track::ntrknprdedx] = f_proton_minus->Eval(trptot);
1585 }
1586
1587 for (SvtxTrack::ConstClusterKeyIter iter_local = tpcseed->begin_cluster_keys();
1588 iter_local != tpcseed->end_cluster_keys();
1589 ++iter_local)
1590 {
1591 TrkrDefs::cluskey cluster_key = *iter_local;
1592 unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
1593 switch (TrkrDefs::getTrkrId(cluster_key))
1594 {
1595 case TrkrDefs::TrkrId::mvtxId:
1596 fX[n_track::ntrknmaps]++;
1597 break;
1598 case TrkrDefs::TrkrId::inttId:
1599 fX[n_track::ntrknintt]++;
1600 break;
1601 case TrkrDefs::TrkrId::tpcId:
1602 fX[n_track::ntrkntpc]++;
1603 if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 8)
1604 {
1605 fX[n_track::ntrkntpc11]++;
1606 }
1607 if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 16)
1608 {
1609 fX[n_track::ntrkntpc1]++;
1610 }
1611 else if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 32)
1612 {
1613 fX[n_track::ntrkntpc2]++;
1614 }
1615 else if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 48)
1616 {
1617 fX[n_track::ntrkntpc3]++;
1618 }
1619 break;
1620 case TrkrDefs::TrkrId::micromegasId:
1621 fX[n_track::ntrknmms]++;
1622 break;
1623 default:
1624 break;
1625 }
1626 }
1627 }
1628 if (silseed)
1629 {
1630 for (SvtxTrack::ConstClusterKeyIter iter_local = silseed->begin_cluster_keys();
1631 iter_local != silseed->end_cluster_keys();
1632 ++iter_local)
1633 {
1634 TrkrDefs::cluskey cluster_key = *iter_local;
1635 unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
1636 switch (TrkrDefs::getTrkrId(cluster_key))
1637 {
1638 case TrkrDefs::TrkrId::mvtxId:
1639 fX[n_track::ntrknmaps]++;
1640 break;
1641 case TrkrDefs::TrkrId::inttId:
1642 fX[n_track::ntrknintt]++;
1643 break;
1644 case TrkrDefs::TrkrId::tpcId:
1645 fX[n_track::ntrkntpc]++;
1646 if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 8)
1647 {
1648 fX[n_track::ntrkntpc11]++;
1649 }
1650 if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 16)
1651 {
1652 fX[n_track::ntrkntpc1]++;
1653 }
1654 else if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 32)
1655 {
1656 fX[n_track::ntrkntpc2]++;
1657 }
1658 else if ((layer_local - (_nlayers_maps + _nlayers_intt)) < 48)
1659 {
1660 fX[n_track::ntrkntpc3]++;
1661 }
1662 break;
1663 case TrkrDefs::TrkrId::micromegasId:
1664 fX[n_track::ntrknmms]++;
1665 break;
1666 default:
1667 break;
1668 }
1669 }
1670 }
1671 fX[n_track::ntrkdca3dxy] = std::numeric_limits<float>::quiet_NaN();
1672 fX[n_track::ntrkdca3dz] = std::numeric_limits<float>::quiet_NaN();
1673 fX[n_track::ntrkdca3dxysigma] = std::numeric_limits<float>::quiet_NaN();
1674 fX[n_track::ntrkdca3dzsigma] = std::numeric_limits<float>::quiet_NaN();
1675 fX[n_track::ntrkdca2d] = std::numeric_limits<float>::quiet_NaN();
1676 fX[n_track::ntrkdca2dsigma] = std::numeric_limits<float>::quiet_NaN();
1677 fX[n_track::ntrkvx] = std::numeric_limits<float>::quiet_NaN();
1678 fX[n_track::ntrkvy] = std::numeric_limits<float>::quiet_NaN();
1679 fX[n_track::ntrkvz] = std::numeric_limits<float>::quiet_NaN();
1680
1681 int vertexID = track->get_vertex_id();
1682 fX[n_track::ntrkvertexID] = vertexID;
1683
1684 if (vertexID >= 0&& vertexmap !=nullptr)
1685 {
1686 if(vertexmap->size()>100000){
1687 std::cout << "too many vtx's" << std::endl;
1688 }
1689
1690
1691
1692
1693
1694
1695
1696
1697
1698
1699
1700
1701
1702
1703
1704
1705
1706
1707
1708
1709
1710
1711
1712 }
1713
1714 fX[n_track::ntrkcharge] = track->get_charge();
1715 float px = track->get_px();
1716 float py = track->get_py();
1717 float pz = track->get_pz();
1718 fX[n_track::ntrkpx] = px;
1719 fX[n_track::ntrkpy] = py;
1720 fX[n_track::ntrkpz] = pz;
1721 TVector3 v(px, py, pz);
1722 fX[n_track::ntrkpt] = v.Pt();
1723 fX[n_track::ntrketa] = v.Eta();
1724 fX[n_track::ntrkphi] = v.Phi();
1725 float CVxx = track->get_error(3, 3);
1726 float CVxy = track->get_error(3, 4);
1727 float CVxz = track->get_error(3, 5);
1728 float CVyy = track->get_error(4, 4);
1729 float CVyz = track->get_error(4, 5);
1730 float CVzz = track->get_error(5, 5);
1731 fX[n_track::ntrkdeltapt] = std::sqrt((CVxx * px * px + 2 * CVxy * px * py + CVyy * py * py) / (px * px + py * py));
1732 fX[n_track::ntrkdeltaeta] = 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)));
1733 fX[n_track::ntrkdeltaphi] = std::sqrt((CVyy * px * px - 2 * CVxy * px * py + CVxx * py * py) / ((px * px + py * py) * (px * px + py * py)));
1734
1735 fX[n_track::ntrkpcax] = track->get_x();
1736 fX[n_track::ntrkpcay] = track->get_y();
1737 fX[n_track::ntrkpcaz] = track->get_z();
1738 }
1739
1740 float TrkrNtuplizer::calc_dedx(TrackSeed *tpcseed){
1741
1742 std::vector<TrkrDefs::cluskey> clusterKeys;
1743 clusterKeys.insert(clusterKeys.end(), tpcseed->begin_cluster_keys(),
1744 tpcseed->end_cluster_keys());
1745
1746 std::vector<float> dedxlist;
1747 for (unsigned long cluster_key : clusterKeys){
1748 unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
1749 if(TrkrDefs::getTrkrId(cluster_key) != TrkrDefs::TrkrId::tpcId){
1750 continue;
1751 }
1752 TrkrCluster* cluster = _cluster_map->findCluster(cluster_key);
1753 int cside = TpcDefs::getSide(cluster_key);
1754 int csector = TpcDefs::getSectorId(cluster_key);
1755 int csegment = 0;
1756 if(layer_local >= 7+32){
1757 csegment = 2;
1758 }else if(layer_local >=7+16){
1759 csegment = 1;
1760 }
1761 float adc = cluster->getAdc();
1762 if(_do_dedx_calib==true){
1763 if(dedxcorr[cside][csector][csegment]!=0){
1764 adc/=dedxcorr[cside][csector][csegment];
1765 }
1766 }
1767 PHG4TpcCylinderGeom* GeoLayer_local = _geom_container->GetLayerCellGeom(layer_local);
1768 float thick = GeoLayer_local->get_thickness();
1769
1770 float r = GeoLayer_local->get_radius();
1771 float alpha = (r * r) / (2 * r * std::abs(1.0 / tpcseed->get_qOverR()));
1772 float beta = std::atan(tpcseed->get_slope());
1773 float alphacorr = std::cos(alpha);
1774 if(alphacorr<0||alphacorr>4){
1775 alphacorr=4;
1776 }
1777 float betacorr = std::cos(beta);
1778 if(betacorr<0||betacorr>4){
1779 betacorr=4;
1780 }
1781 adc/=thick;
1782 adc*=alphacorr;
1783 adc*=betacorr;
1784 dedxlist.push_back(adc);
1785 sort(dedxlist.begin(), dedxlist.end());
1786 }
1787 int trunc_min = 0;
1788 int trunc_max = (int)dedxlist.size()*0.7;
1789 float sumdedx = 0;
1790 int ndedx = 0;
1791 for(int j = trunc_min; j<=trunc_max;j++){
1792 sumdedx+=dedxlist.at(j);
1793 ndedx++;
1794 }
1795 sumdedx/=ndedx;
1796 return sumdedx;
1797 }
1798
1799 float TrkrNtuplizer::get_n1pix(TrackSeed *tpcseed){
1800
1801 std::vector<TrkrDefs::cluskey> clusterKeys;
1802 clusterKeys.insert(clusterKeys.end(), tpcseed->begin_cluster_keys(),
1803 tpcseed->end_cluster_keys());
1804
1805 float n1pix = 0;
1806 for (unsigned long cluster_key : clusterKeys){
1807 TrkrCluster* cluster = _cluster_map->findCluster(cluster_key);
1808 if(cluster->getPhiSize()==1){
1809 n1pix++;
1810 }
1811 }
1812 return n1pix;
1813 }
1814
1815 void TrkrNtuplizer::FillCluster(float fXcluster[n_cluster::clusize], TrkrDefs::cluskey cluster_key)
1816 {
1817 unsigned int layer_local = TrkrDefs::getLayer(cluster_key);
1818 TrkrCluster* cluster = _cluster_map->findCluster(cluster_key);
1819
1820 Acts::Vector3 cglob;
1821 cglob = _tgeometry->getGlobalPosition(cluster_key, cluster);
1822 float x = cglob(0);
1823 float y = cglob(1);
1824 float z = cglob(2);
1825 TVector3 pos(x, y, z);
1826 float r = pos.Perp();
1827 float phi = pos.Phi();
1828 auto para_errors = _ClusErrPara.get_clusterv5_modified_error(cluster, r, cluster_key);
1829 float phibin = std::numeric_limits<float>::quiet_NaN();
1830 float tbin = std::numeric_limits<float>::quiet_NaN();
1831 float locx = cluster->getLocalX();
1832 float locy = cluster->getLocalY();
1833 fXcluster[n_cluster::ncludcal] = 1;
1834 if (layer_local >= _nlayers_maps + _nlayers_intt && layer_local < _nlayers_maps + _nlayers_intt + _nlayers_tpc)
1835 {
1836 int side_tpc = TpcDefs::getSide(cluster_key);
1837 PHG4TpcCylinderGeom* GeoLayer_local = _geom_container->GetLayerCellGeom(layer_local);
1838 phibin = GeoLayer_local->get_pad_float(phi, side_tpc);
1839 tbin = GeoLayer_local->get_tbin_float(locy - 39.6);
1840 int cside = TpcDefs::getSide(cluster_key);
1841 int csector = TpcDefs::getSectorId(cluster_key);
1842 int csegment = 0;
1843 if(layer_local >= 7+32){
1844 csegment = 2;
1845 }else if(layer_local >=7+16){
1846 csegment = 1;
1847 }
1848 if(_do_dedx_calib==true){
1849 fXcluster[n_cluster::ncludcal] = dedxcorr[cside][csector][csegment];
1850 }else{
1851 fXcluster[n_cluster::ncludcal] = 1;
1852 }
1853 }
1854 else
1855 {
1856 phibin = locx;
1857 tbin = locy;
1858 }
1859 fXcluster[n_cluster::nclulocx] = locx;
1860 fXcluster[n_cluster::nclulocy] = locy;
1861 fXcluster[n_cluster::nclux] = cglob(0);
1862 fXcluster[n_cluster::ncluy] = cglob(1);
1863 fXcluster[n_cluster::ncluz] = cglob(2);
1864 fXcluster[n_cluster::nclur] = pos.Perp();
1865 fXcluster[n_cluster::ncluphi] = pos.Phi();
1866 fXcluster[n_cluster::nclueta] = pos.Eta();
1867 fXcluster[n_cluster::nclutheta] = pos.Theta();
1868
1869 fXcluster[n_cluster::ncluphibin] = phibin;
1870 fXcluster[n_cluster::nclutbin] = tbin;
1871
1872 fXcluster[n_cluster::ncluex] = std::numeric_limits<float>::quiet_NaN();
1873 fXcluster[n_cluster::ncluey] = std::numeric_limits<float>::quiet_NaN();
1874 fXcluster[n_cluster::ncluez] = sqrt(para_errors.second);
1875 fXcluster[n_cluster::ncluephi] = sqrt(para_errors.first);
1876 fXcluster[n_cluster::nclupez] = std::numeric_limits<float>::quiet_NaN();
1877 fXcluster[n_cluster::nclupephi] = std::numeric_limits<float>::quiet_NaN();
1878 fXcluster[n_cluster::nclue] = cluster->getAdc();
1879 fXcluster[n_cluster::ncluadc] = cluster->getAdc();
1880 fXcluster[n_cluster::nclumaxadc] = cluster->getMaxAdc();
1881 fXcluster[n_cluster::nclulayer] = layer_local;
1882
1883 if (layer_local < 3)
1884 {
1885 fXcluster[n_cluster::ncluphielem] = MvtxDefs::getStaveId(cluster_key);
1886 fXcluster[n_cluster::ncluzelem] = MvtxDefs::getChipId(cluster_key);
1887 }
1888 if (layer_local >= 3 && layer_local < 7)
1889 {
1890 fXcluster[n_cluster::ncluphielem] = InttDefs::getLadderPhiId(cluster_key);
1891 fXcluster[n_cluster::ncluzelem] = InttDefs::getLadderZId(cluster_key);
1892 }
1893 if (layer_local >= 7 && layer_local < 55)
1894 {
1895 fXcluster[n_cluster::ncluphielem] = TpcDefs::getSectorId(cluster_key);
1896 fXcluster[n_cluster::ncluzelem] = TpcDefs::getSide(cluster_key);
1897 }
1898
1899
1900
1901
1902
1903
1904
1905
1906
1907
1908 fXcluster[n_cluster::nclusize] = cluster->getSize();
1909 fXcluster[n_cluster::ncluphisize] = cluster->getPhiSize();
1910 fXcluster[n_cluster::ncluzsize] = cluster->getZSize();
1911 fXcluster[n_cluster::nclupedge] = cluster->getEdge();
1912 if (layer_local == 7 || layer_local == 22 ||
1913 layer_local == 23 || layer_local == 28 ||
1914 layer_local == 39 || layer_local == 54)
1915 {
1916 fXcluster[n_cluster::ncluredge] = 1;
1917 }
1918
1919 fXcluster[n_cluster::ncluovlp] = 3;
1920 fXcluster[n_cluster::nclutrackID] = std::numeric_limits<float>::quiet_NaN();
1921 fXcluster[n_cluster::ncluniter] = 0;
1922
1923 return;
1924 }
1925
1926 std::vector<TrkrDefs::cluskey> TrkrNtuplizer::get_track_ckeys(SvtxTrack* track)
1927 {
1928 std::vector<TrkrDefs::cluskey> cluster_keys;
1929 TrackSeed* tpcseed = track->get_tpc_seed();
1930 TrackSeed* silseed = track->get_silicon_seed();
1931 if (silseed)
1932 {
1933 for (auto iter = silseed->begin_cluster_keys();
1934 iter != silseed->end_cluster_keys();
1935 ++iter)
1936 {
1937 cluster_keys.push_back(*iter);
1938 }
1939 }
1940 if (tpcseed)
1941 {
1942 for (auto iter = tpcseed->begin_cluster_keys();
1943 iter != tpcseed->end_cluster_keys();
1944 ++iter)
1945 {
1946 cluster_keys.push_back(*iter);
1947 }
1948 }
1949
1950 return cluster_keys;
1951 }
1952
1953 SvtxTrack* TrkrNtuplizer::best_track_from(TrkrDefs::cluskey cluster_key)
1954 {
1955 std::map<TrkrDefs::cluskey, SvtxTrack*>::iterator find_iter =
1956 _cache_best_track_from_cluster.find(cluster_key);
1957 if (find_iter != _cache_best_track_from_cluster.end())
1958 {
1959 return find_iter->second;
1960 }
1961
1962 SvtxTrack* best_track = nullptr;
1963 float best_quality = FLT_MAX;
1964
1965 std::set<SvtxTrack*> tracks = all_tracks_from(cluster_key);
1966
1967 for (auto *candidate : tracks)
1968 {
1969 if (candidate->get_quality() < best_quality)
1970 {
1971 best_quality = candidate->get_quality();
1972 best_track = candidate;
1973 }
1974 }
1975
1976 _cache_best_track_from_cluster.insert(std::make_pair(cluster_key, best_track));
1977 return best_track;
1978 }
1979
1980 std::set<SvtxTrack*> TrkrNtuplizer::all_tracks_from(TrkrDefs::cluskey cluster_key)
1981 {
1982 std::set<SvtxTrack*> tracks;
1983
1984 if (_cache_track_from_cluster_exists == false)
1985 {
1986 create_cache_track_from_cluster();
1987 }
1988 std::map<TrkrDefs::cluskey, std::set<SvtxTrack*> >::iterator find_iter =
1989 _cache_all_tracks_from_cluster.find(cluster_key);
1990 if (find_iter != _cache_all_tracks_from_cluster.end())
1991 {
1992 return find_iter->second;
1993 }
1994
1995 return tracks;
1996
1997
1998
1999 for (auto& iter : *_trackmap)
2000 {
2001 SvtxTrack* track = iter.second;
2002 std::vector<TrkrDefs::cluskey> cluster_keys = get_track_ckeys(track);
2003
2004
2005 for (const auto& candidate : cluster_keys)
2006 {
2007
2008
2009
2010
2011
2012
2013
2014
2015
2016
2017 if (cluster_key == candidate)
2018 {
2019 tracks.insert(track);
2020 }
2021 }
2022 }
2023
2024 _cache_all_tracks_from_cluster.insert(std::make_pair(cluster_key, tracks));
2025
2026 return tracks;
2027 }
2028
2029 void TrkrNtuplizer::create_cache_track_from_cluster()
2030 {
2031 if (!_trackmap)
2032 {
2033 return;
2034 }
2035
2036
2037 for (auto& iter : *_trackmap)
2038 {
2039 SvtxTrack* track = iter.second;
2040 std::vector<TrkrDefs::cluskey> cluster_keys = get_track_ckeys(track);
2041
2042
2043 for (const auto& candidate_key : cluster_keys)
2044 {
2045
2046 std::map<TrkrDefs::cluskey, std::set<SvtxTrack*> >::iterator cliter =
2047 _cache_all_tracks_from_cluster.find(candidate_key);
2048 if (cliter != _cache_all_tracks_from_cluster.end())
2049 {
2050 cliter->second.insert(track);
2051 }
2052 else
2053 {
2054 std::set<SvtxTrack*> tracks;
2055 tracks.insert(track);
2056 _cache_all_tracks_from_cluster.insert(std::make_pair(candidate_key, tracks));
2057 }
2058 }
2059 }
2060 _cache_track_from_cluster_exists = true;
2061
2062 return;
2063 }
2064
2065 TMatrixF TrkrNtuplizer::calculateClusterError(TrkrCluster* c, float& clusphi)
2066 {
2067 TMatrixF localErr(3, 3);
2068 localErr[0][0] = 0.;
2069 localErr[0][1] = 0.;
2070 localErr[0][2] = 0.;
2071 localErr[1][0] = 0.;
2072 localErr[1][1] = c->getActsLocalError(0, 0);
2073 localErr[1][2] = c->getActsLocalError(0, 1);
2074 localErr[2][0] = 0.;
2075 localErr[2][1] = c->getActsLocalError(1, 0);
2076 localErr[2][2] = c->getActsLocalError(1, 1);
2077
2078 TMatrixF ROT(3, 3);
2079 ROT[0][0] = std::cos(clusphi);
2080 ROT[0][1] = -std::sin(clusphi);
2081 ROT[0][2] = 0.0;
2082 ROT[1][0] = std::sin(clusphi);
2083 ROT[1][1] = std::cos(clusphi);
2084 ROT[1][2] = 0.0;
2085 ROT[2][0] = 0.0;
2086 ROT[2][1] = 0.0;
2087 ROT[2][2] = 1.0;
2088 TMatrixF ROT_T(3, 3);
2089 ROT_T.Transpose(ROT);
2090
2091 TMatrixF err(3, 3);
2092 err = ROT * localErr * ROT_T;
2093 return err;
2094 }
