packages/TrackingDiagnostics/TrkrNtuplizer.cc

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