STrackMatchingComparator/src/STrackMatcherComparator.cc

0001 // ----------------------------------------------------------------------------
0002 // 'STrackMatcherComparator.cc'
0003 // Derek Anderson
0004 // 01.30.2024
0005 //
0006 // Small module to produce plots from the output of the
0007 // SvtxEvaluator and FillMatchedClusTree modules.
0008 // ----------------------------------------------------------------------------
0009
0010 // c++ utilities
0011 #include <cmath>
0012 #include <array>
0013 #include <string>
0014 #include <vector>
0015 #include <cassert>
0016 #include <utility>
0017 #include <iostream>
0018 // root libraries
0019 #include <TH1.h>
0020 #include <TH2.h>
0021 #include <TPad.h>
0022 #include <TFile.h>
0023 #include <TLine.h>
0024 #include <TError.h>
0025 #include <TString.h>
0026 #include <TLegend.h>
0027 #include <TCanvas.h>
0028 #include <TPaveText.h>
0029 #include <TDirectory.h>
0030 // analysis utilities
0031 #include "STrackMatcherComparator.h"
0032 #include "STrackMatcherComparatorConfig.h"
0033 #include "STrackMatcherComparatorHistDef.h"
0034 #include "STrackMatcherComparatorHistContent.h"
0035
0036 // make common namespaces implicit
0037 using namespace std;
0038
0039
0040
0041 // ctor/dtor ------------------------------------------------------------------
0042
0043 STrackMatcherComparator::STrackMatcherComparator(optional<STrackMatcherComparatorConfig> config) {
0044
0045   if (config.has_value()) {
0046     m_config = config.value();
0047   }
0048
0049   // make sure vectors are empty
0050   m_vecHistDirs.clear();
0051   m_vecRatioDirs.clear();
0052   m_vecPlotDirs.clear();
0053   m_vecTreeHists1D.clear();
0054   m_vecTupleHists1D.clear();
0055   m_vecOldHists1D.clear();
0056   m_vecTreeHists2D.clear();
0057   m_vecTupleHists2D.clear();
0058   m_vecOldHists2D.clear();
0059
0060 }  // end ctor
0061
0062
0063
0064 STrackMatcherComparator::~STrackMatcherComparator() {
0065
0066   /* nothing to do */
0067
0068 }  // end dtor
0069
0070
0071
0072 // public methods -------------------------------------------------------------
0073
0074 void STrackMatcherComparator::Init() {
0075
0076   // announce start
0077   cout << "\n  Beginning track matcher comparison...\n"
0078        << "    Initialization:"
0079        << endl;
0080
0081   // run initialization routines
0082   OpenOutput();
0083   OpenInput();
0084   InitHists();
0085   return;
0086
0087 }  // end 'Init()'
0088
0089
0090
0091 void STrackMatcherComparator::Analyze() {
0092
0093   // announce analyzing
0094   cout << "    Analysis:" << endl;
0095
0096   // run analysis routines
0097   GetNewTreeHists();
0098   GetNewTupleHists();
0099   GetOldTupleHists();
0100   return;
0101
0102 }  // end 'Analysis()'
0103
0104
0105
0106 void STrackMatcherComparator::End() {
0107
0108   // announce end
0109   cout << "    End:" << endl;
0110
0111   // run ending routines
0112   MakeRatiosAndPlots(m_vecTreeHists1D,  m_vecTreeHists2D,  Src::NewTree,  "VsNewTree");
0113   MakeRatiosAndPlots(m_vecTupleHists1D, m_vecTupleHists2D, Src::NewTuple, "VsNewTuple");
0114   SaveHistograms();
0115   CloseInput();
0116   CloseOutput();
0117
0118   // exit
0119   cout << "  Finished track matcher comparison!\n" << endl;
0120   return;
0121 }
0122
0123
0124
0125 // internal methods -----------------------------------------------------------
0126
0127 void STrackMatcherComparator::OpenOutput() {
0128
0129   // open output files
0130   m_outFile = new TFile(m_config.outFileName.data(), "recreate");
0131   if (!m_outFile) {
0132     cerr << "PANIC: couldn't open output file!" << endl;
0133     assert(m_outFile);
0134   }
0135   cout << "      Opened output file." << endl;
0136
0137   // create output directories
0138   m_vecHistDirs.resize(m_const.nDirHist);
0139   m_vecRatioDirs.resize(m_const.nDirRatio);
0140   m_vecPlotDirs.resize(m_const.nDirPlot);
0141
0142   m_vecHistDirs[Src::NewTree]   = (TDirectory*) m_outFile -> mkdir(m_config.histDirNames.at(Src::NewTree).data());
0143   m_vecHistDirs[Src::NewTuple]  = (TDirectory*) m_outFile -> mkdir(m_config.histDirNames.at(Src::NewTuple).data());
0144   m_vecHistDirs[Src::OldTuple]  = (TDirectory*) m_outFile -> mkdir(m_config.histDirNames.at(Src::OldTuple).data());
0145   m_vecRatioDirs[Src::NewTree]  = (TDirectory*) m_outFile -> mkdir(m_config.ratioDirNames.at(Src::NewTree).data());
0146   m_vecRatioDirs[Src::NewTuple] = (TDirectory*) m_outFile -> mkdir(m_config.ratioDirNames.at(Src::NewTuple).data());
0147   m_vecPlotDirs[Src::NewTree]   = (TDirectory*) m_outFile -> mkdir(m_config.plotDirNames.at(Src::NewTree).data());
0148   m_vecPlotDirs[Src::NewTuple]  = (TDirectory*) m_outFile -> mkdir(m_config.plotDirNames.at(Src::NewTuple).data());
0149
0150   // announce end and return
0151   cout << "      Created output directories." << endl;
0152   return;
0153
0154 }  // end 'OpenOutput()'
0155
0156
0157
0158 void STrackMatcherComparator::OpenInput() {
0159
0160   // open true files
0161   m_treeInFileTrue  = new TFile(m_config.newInFileName.data(), "read");
0162   m_tupleInFileTrue = new TFile(m_config.newInFileName.data(), "read");
0163   m_oldInFileTrue   = new TFile(m_config.oldInFileName.data(), "read");
0164   if (!m_treeInFileTrue || !m_tupleInFileTrue || !m_oldInFileTrue) {
0165     cerr << "PANIC: couldn't open an input file!\n"
0166          << "       m_treeInFileTrue  = " << m_treeInFileTrue  << "\n"
0167          << "       m_tupleInFileTrue = " << m_tupleInFileTrue << "\n"
0168          << "       m_oldInFileTrue   = " << m_oldInFileTrue
0169          << endl;
0170     assert(m_treeInFileTrue && m_tupleInFileTrue && m_oldInFileTrue);
0171   }
0172   cout << "      Opened truth input files." << endl;
0173
0174   // open reco files
0175   m_treeInFileReco  = new TFile(m_config.newInFileName.data(), "read");
0176   m_tupleInFileReco = new TFile(m_config.newInFileName.data(), "read");
0177   m_oldInFileReco   = new TFile(m_config.oldInFileName.data(), "read");
0178   if (!m_treeInFileReco || !m_tupleInFileReco || !m_oldInFileReco) {
0179     cerr << "PANIC: couldn't open an input file!\n"
0180          << "       m_treeInFileReco  = " << m_treeInFileReco  << "\n"
0181          << "       m_tupleInFileReco = " << m_tupleInFileReco << "\n"
0182          << "       m_oldInFileReco   = " << m_oldInFileReco
0183          << endl;
0184     assert(m_treeInFileReco && m_tupleInFileReco && m_oldInFileReco);
0185   }
0186   cout << "      Opened reco input files." << endl;
0187
0188   // grab true input trees/tuples
0189   m_tTreeTrue   = (TTree*)   m_treeInFileTrue  -> Get(m_config.newTreeTrueName.data());
0190   m_ntTupleTrue = (TNtuple*) m_tupleInFileTrue -> Get(m_config.newTupleTrueName.data());
0191   m_ntOldTrue   = (TNtuple*) m_oldInFileTrue   -> Get(m_config.oldTupleTrueName.data());
0192   if (!m_tTreeTrue || !m_ntTupleTrue || !m_ntOldTrue) {
0193     cerr << "PANIC: couldn't grab an input truth tree/tuple!\n"
0194          << "       m_tTreeTrue   = " << m_tTreeTrue   << "\n"
0195          << "       m_ntTupleTrue = " << m_ntTupleTrue << "\n"
0196          << "       m_ntOldTrue   = " << m_ntOldTrue
0197          << endl;
0198     assert(m_tTreeTrue && m_ntTupleTrue && m_ntOldTrue);
0199   }
0200   cout << "      Grabbed input truth trees/tuples." << endl;
0201
0202   // grab reco input trees/tuples
0203   m_tTreeReco   = (TTree*)   m_treeInFileReco  -> Get(m_config.newTreeRecoName.data());
0204   m_ntTupleReco = (TNtuple*) m_tupleInFileReco -> Get(m_config.newTupleRecoName.data());
0205   m_ntOldReco   = (TNtuple*) m_oldInFileReco   -> Get(m_config.oldTupleRecoName.data());
0206   if (!m_tTreeReco || !m_ntTupleReco || !m_ntOldReco) {
0207     cerr << "PANIC: couldn't grab an input truth tree/tuple!\n"
0208          << "       m_tTreeReco   = " << m_tTreeReco   << "\n"
0209          << "       m_ntTupleReco = " << m_ntTupleReco << "\n"
0210          << "       m_ntOldReco   = " << m_ntOldReco
0211          << endl;
0212     assert(m_tTreeReco && m_ntTupleReco && m_ntOldReco);
0213   }
0214
0215   // announce end and return
0216   cout << "      Grabbed input truth trees/tuples." << endl;
0217   return;
0218
0219 }  // end 'OpenInput()'
0220
0221
0222
0223 void STrackMatcherComparator::InitHists() {
0224
0225   // announce start of routine
0226   const size_t nHistRows   = m_hist.vecNameBase.size();
0227   const size_t nHistTypes  = m_hist.vecNameBase[0].size();
0228   const size_t nVsHistMods = m_hist.vecVsModifiers.size();
0229   cout << "      Initializing histograms." << endl;
0230
0231   // output histogram base binning definitions
0232   vector<tuple<uint32_t, pair<float, float>>> vecBaseHistBins = m_hist.GetVecHistBins();
0233   vector<tuple<uint32_t, pair<float, float>>> vecVsHistBins   = m_hist.GetVecVsHistBins();
0234
0235   // turn on errors
0236   TH1::SetDefaultSumw2(true);
0237   TH2::SetDefaultSumw2(true);
0238
0239   // set up 1d output histograms
0240   m_vecTreeHists1D.resize(nHistRows);
0241   m_vecTupleHists1D.resize(nHistRows);
0242   m_vecOldHists1D.resize(nHistRows);
0243   for (size_t iHistRow = 0; iHistRow < nHistRows; iHistRow++) {
0244     for (size_t iHistType = 0; iHistType < nHistTypes; iHistType++) {
0245
0246       // create names
0247       const string sTreeHist1D  = m_hist.vecNameBase[iHistRow][iHistType] + "_" + m_config.newTreeLabel;
0248       const string sTupleHist1D = m_hist.vecNameBase[iHistRow][iHistType] + "_" + m_config.newTupleLabel;
0249       const string sOldHist1D   = m_hist.vecNameBase[iHistRow][iHistType] + "_" + m_config.oldTupleLabel;
0250
0251       // create output histograms
0252       m_vecTreeHists1D[iHistRow].push_back(
0253         new TH1D(
0254           sTreeHist1D.data(),
0255           "",
0256           get<0>(vecBaseHistBins[iHistRow]),
0257           get<1>(vecBaseHistBins[iHistRow]).first,
0258           get<1>(vecBaseHistBins[iHistRow]).second
0259         )
0260       );
0261       m_vecTupleHists1D[iHistRow].push_back(
0262         new TH1D(
0263           sTupleHist1D.data(),
0264           "",
0265           get<0>(vecBaseHistBins[iHistRow]),
0266           get<1>(vecBaseHistBins[iHistRow]).first,
0267           get<1>(vecBaseHistBins[iHistRow]).second
0268         )
0269       );
0270       m_vecOldHists1D[iHistRow].push_back(
0271         new TH1D(
0272           sOldHist1D.data(),
0273           "",
0274           get<0>(vecBaseHistBins[iHistRow]),
0275           get<1>(vecBaseHistBins[iHistRow]).first,
0276           get<1>(vecBaseHistBins[iHistRow]).second
0277         )
0278       );
0279
0280       // set axis titles
0281       m_vecTreeHists1D[iHistRow].back()  -> GetXaxis() -> SetTitle(m_hist.vecBaseAxisVars.at(iHistRow).data());
0282       m_vecTupleHists1D[iHistRow].back() -> GetXaxis() -> SetTitle(m_hist.vecBaseAxisVars.at(iHistRow).data());
0283       m_vecOldHists1D[iHistRow].back()   -> GetXaxis() -> SetTitle(m_hist.vecBaseAxisVars.at(iHistRow).data());
0284       m_vecTreeHists1D[iHistRow].back()  -> GetYaxis() -> SetTitle(m_config.count.data());
0285       m_vecTupleHists1D[iHistRow].back() -> GetYaxis() -> SetTitle(m_config.count.data());
0286       m_vecOldHists1D[iHistRow].back()   -> GetYaxis() -> SetTitle(m_config.count.data());
0287     }  // end type loop
0288   }  // end row loop
0289   cout << "      Initialized 1D histograms." << endl;
0290
0291   // set up 2d output histograms
0292   m_vecTreeHists2D.resize(nHistRows);
0293   m_vecTupleHists2D.resize(nHistRows);
0294   m_vecOldHists2D.resize(nHistRows);
0295   for (size_t iHistRow = 0; iHistRow < nHistRows; iHistRow++) {
0296
0297     // reserve space for all histograms in row
0298     m_vecTreeHists2D[iHistRow].resize(nHistTypes);
0299     m_vecTupleHists2D[iHistRow].resize(nHistTypes);
0300     m_vecOldHists2D[iHistRow].resize(nHistTypes);
0301     for (size_t iHistType = 0; iHistType < nHistTypes; iHistType++) {
0302       for (size_t iVsHistMod = 0; iVsHistMod < nVsHistMods; iVsHistMod++) {
0303
0304         // create names
0305         const string sTreeHist2D  = m_hist.vecNameBase[iHistRow][iHistType] + m_hist.vecVsModifiers[iVsHistMod] + "_" + m_config.newTreeLabel;
0306         const string sTupleHist2D = m_hist.vecNameBase[iHistRow][iHistType] + m_hist.vecVsModifiers[iVsHistMod] + "_" + m_config.newTupleLabel;
0307         const string sOldHist2D   = m_hist.vecNameBase[iHistRow][iHistType] + m_hist.vecVsModifiers[iVsHistMod] + "_" + m_config.oldTupleLabel;
0308
0309         // create output histograms
0310         m_vecTreeHists2D[iHistRow][iHistType].push_back(
0311           new TH2D(
0312             sTreeHist2D.data(),
0313             "",
0314             get<0>(vecVsHistBins[iVsHistMod]),
0315             get<1>(vecVsHistBins[iVsHistMod]).first,
0316             get<1>(vecVsHistBins[iVsHistMod]).second,
0317             get<0>(vecBaseHistBins[iHistRow]),
0318             get<1>(vecBaseHistBins[iHistRow]).first,
0319             get<1>(vecBaseHistBins[iHistRow]).second
0320           )
0321         );
0322         m_vecTupleHists2D[iHistRow][iHistType].push_back(
0323           new TH2D(
0324             sTupleHist2D.data(),
0325             "",
0326             get<0>(vecVsHistBins[iVsHistMod]),
0327             get<1>(vecVsHistBins[iVsHistMod]).first,
0328             get<1>(vecVsHistBins[iVsHistMod]).second,
0329             get<0>(vecBaseHistBins[iHistRow]),
0330             get<1>(vecBaseHistBins[iHistRow]).first,
0331             get<1>(vecBaseHistBins[iHistRow]).second
0332           )
0333         );
0334         m_vecOldHists2D[iHistRow][iHistType].push_back(
0335           new TH2D(
0336             sOldHist2D.data(),
0337             "",
0338             get<0>(vecVsHistBins[iVsHistMod]),
0339             get<1>(vecVsHistBins[iVsHistMod]).first,
0340             get<1>(vecVsHistBins[iVsHistMod]).second,
0341             get<0>(vecBaseHistBins[iHistRow]),
0342             get<1>(vecBaseHistBins[iHistRow]).first,
0343             get<1>(vecBaseHistBins[iHistRow]).second
0344           )
0345         );
0346
0347         // set axis titles
0348         m_vecTreeHists2D[iHistRow][iHistType].back()  -> GetXaxis() -> SetTitle(m_hist.vecVsAxisVars.at(iVsHistMod).data());
0349         m_vecTupleHists2D[iHistRow][iHistType].back() -> GetXaxis() -> SetTitle(m_hist.vecVsAxisVars.at(iVsHistMod).data());
0350         m_vecOldHists2D[iHistRow][iHistType].back()   -> GetXaxis() -> SetTitle(m_hist.vecVsAxisVars.at(iVsHistMod).data());
0351         m_vecTreeHists2D[iHistRow][iHistType].back()  -> GetYaxis() -> SetTitle(m_hist.vecBaseAxisVars.at(iHistRow).data());
0352         m_vecTupleHists2D[iHistRow][iHistType].back() -> GetYaxis() -> SetTitle(m_hist.vecBaseAxisVars.at(iHistRow).data());
0353         m_vecOldHists2D[iHistRow][iHistType].back()   -> GetYaxis() -> SetTitle(m_hist.vecBaseAxisVars.at(iHistRow).data());
0354         m_vecTreeHists2D[iHistRow][iHistType].back()  -> GetZaxis() -> SetTitle(m_config.count.data());
0355         m_vecTupleHists2D[iHistRow][iHistType].back() -> GetZaxis() -> SetTitle(m_config.count.data());
0356         m_vecOldHists2D[iHistRow][iHistType].back()   -> GetZaxis() -> SetTitle(m_config.count.data());
0357       }  // end modifier loop
0358     }  // end type loop
0359   }  // end row loop
0360
0361   // announce end and exit
0362   cout << "      Initialized 2D output histograms." << endl;
0363   return;
0364
0365 }  // end 'InitHists()'
0366
0367
0368
0369 void STrackMatcherComparator::GetNewTreeHists() {
0370
0371   // announce start of routine
0372   cout << "      Grabbing new matcher tree histograms:" << endl;
0373
0374   // declare input leaves (ignore cluster branches for now)
0375   int   tru_event;
0376   int   tru_nphg4_part;
0377   float tru_centrality;
0378   int   tru_ntrackmatches;
0379   int   tru_nphg4;
0380   int   tru_nsvtx;
0381   int   tru_trackid;
0382   bool  tru_is_G4track;
0383   bool  tru_is_Svtrack;
0384   bool  tru_is_matched;
0385   float tru_trkpt;
0386   float tru_trketa;
0387   float tru_trkphi;
0388   int   tru_nclus;
0389   int   tru_nclustpc;
0390   int   tru_nclusmvtx;
0391   int   tru_nclusintt;
0392   float tru_matchrat;
0393   float tru_matchrat_intt;
0394   float tru_matchrat_mvtx;
0395   float tru_matchrat_tpc;
0396
0397   int   rec_event;
0398   int   rec_nphg4_part;
0399   float rec_centrality;
0400   int   rec_ntrackmatches;
0401   int   rec_nphg4;
0402   int   rec_nsvtx;
0403   int   rec_trackid;
0404   bool  rec_is_G4track;
0405   bool  rec_is_Svtrack;
0406   bool  rec_is_matched;
0407   float rec_trkpt;
0408   float rec_trketa;
0409   float rec_trkphi;
0410   int   rec_nclus;
0411   int   rec_nclustpc;
0412   int   rec_nclusmvtx;
0413   int   rec_nclusintt;
0414   float rec_matchrat;
0415   float rec_matchrat_intt;
0416   float rec_matchrat_mvtx;
0417   float rec_matchrat_tpc;
0418
0419   // set branch addresses (ignore cluster branches for now)
0420   m_tTreeTrue -> SetBranchAddress("event",         &tru_event);
0421   m_tTreeTrue -> SetBranchAddress("nphg4_part",    &tru_nphg4_part);
0422   m_tTreeTrue -> SetBranchAddress("centrality",    &tru_centrality);
0423   m_tTreeTrue -> SetBranchAddress("ntrackmatches", &tru_ntrackmatches);
0424   m_tTreeTrue -> SetBranchAddress("nphg4",         &tru_nphg4);
0425   m_tTreeTrue -> SetBranchAddress("nsvtx",         &tru_nsvtx);
0426   m_tTreeTrue -> SetBranchAddress("trackid",       &tru_trackid);
0427   m_tTreeTrue -> SetBranchAddress("is_G4track",    &tru_is_G4track);
0428   m_tTreeTrue -> SetBranchAddress("is_Svtrack",    &tru_is_Svtrack);
0429   m_tTreeTrue -> SetBranchAddress("is_matched",    &tru_is_matched);
0430   m_tTreeTrue -> SetBranchAddress("trkpt",         &tru_trkpt);
0431   m_tTreeTrue -> SetBranchAddress("trketa",        &tru_trketa);
0432   m_tTreeTrue -> SetBranchAddress("trkphi",        &tru_trkphi);
0433   m_tTreeTrue -> SetBranchAddress("nclus",         &tru_nclus);
0434   m_tTreeTrue -> SetBranchAddress("nclustpc",      &tru_nclustpc);
0435   m_tTreeTrue -> SetBranchAddress("nclusmvtx",     &tru_nclusmvtx);
0436   m_tTreeTrue -> SetBranchAddress("nclusintt",     &tru_nclusintt);
0437   m_tTreeTrue -> SetBranchAddress("matchrat",      &tru_matchrat);
0438   m_tTreeTrue -> SetBranchAddress("matchrat_intt", &tru_matchrat_intt);
0439   m_tTreeTrue -> SetBranchAddress("matchrat_mvtx", &tru_matchrat_mvtx);
0440   m_tTreeTrue -> SetBranchAddress("matchrat_tpc",  &tru_matchrat_tpc);
0441
0442   m_tTreeReco -> SetBranchAddress("event",         &rec_event);
0443   m_tTreeReco -> SetBranchAddress("nphg4_part",    &rec_nphg4_part);
0444   m_tTreeReco -> SetBranchAddress("centrality",    &rec_centrality);
0445   m_tTreeReco -> SetBranchAddress("ntrackmatches", &rec_ntrackmatches);
0446   m_tTreeReco -> SetBranchAddress("nphg4",         &rec_nphg4);
0447   m_tTreeReco -> SetBranchAddress("nsvtx",         &rec_nsvtx);
0448   m_tTreeReco -> SetBranchAddress("trackid",       &rec_trackid);
0449   m_tTreeReco -> SetBranchAddress("is_G4track",    &rec_is_G4track);
0450   m_tTreeReco -> SetBranchAddress("is_Svtrack",    &rec_is_Svtrack);
0451   m_tTreeReco -> SetBranchAddress("is_matched",    &rec_is_matched);
0452   m_tTreeReco -> SetBranchAddress("trkpt",         &rec_trkpt);
0453   m_tTreeReco -> SetBranchAddress("trketa",        &rec_trketa);
0454   m_tTreeReco -> SetBranchAddress("trkphi",        &rec_trkphi);
0455   m_tTreeReco -> SetBranchAddress("nclus",         &rec_nclus);
0456   m_tTreeReco -> SetBranchAddress("nclustpc",      &rec_nclustpc);
0457   m_tTreeReco -> SetBranchAddress("nclusmvtx",     &rec_nclusmvtx);
0458   m_tTreeReco -> SetBranchAddress("nclusintt",     &rec_nclusintt);
0459   m_tTreeReco -> SetBranchAddress("matchrat",      &rec_matchrat);
0460   m_tTreeReco -> SetBranchAddress("matchrat_intt", &rec_matchrat_intt);
0461   m_tTreeReco -> SetBranchAddress("matchrat_mvtx", &rec_matchrat_mvtx);
0462   m_tTreeReco -> SetBranchAddress("matchrat_tpc",  &rec_matchrat_tpc);
0463   cout << "        Set input branches." << endl;
0464
0465   // grab no. of entries
0466   const int64_t nTrueEntries = m_tTreeTrue -> GetEntries();
0467   const int64_t nRecoEntries = m_tTreeReco -> GetEntries();
0468   cout << "        Beginning truth particle loop: " << nTrueEntries << " to process" << endl;
0469
0470   // loop over truth particles
0471   int64_t nTrueBytes = 0;
0472   for (int64_t iTrueEntry = 0; iTrueEntry < nTrueEntries; iTrueEntry++) {
0473
0474     // grab truth particle entry
0475     const int64_t trueBytes = m_tTreeTrue -> GetEntry(iTrueEntry);
0476     if (trueBytes < 0) {
0477       cerr << "PANIC: issue with entry " << iTrueEntry << "! Aborting loop!\n" << endl;
0478       break;
0479     } else {
0480       nTrueBytes += trueBytes;
0481     }
0482
0483     const int64_t iTrueProg = iTrueEntry + 1;
0484     if (iTrueProg == nTrueEntries) {
0485       cout << "          Processing entry " << iTrueProg << "/" << nTrueEntries << "..." << endl;
0486     } else {
0487       cout << "          Processing entry " << iTrueProg << "/" << nTrueEntries << "...\r" << flush;
0488     }
0489
0490     // select truth particles
0491     //   - FIXME add cuts
0492     if (!tru_is_G4track) continue;
0493
0494     // bundle histogram values to fill
0495     const STrackMatcherComparatorHistContent content {
0496       (double) tru_nclus,
0497       (double) tru_nclusintt,
0498       (double) tru_nclusmvtx,
0499       (double) tru_nclustpc,
0500       1.,
0501       1.,
0502       1.,
0503       1.,
0504       tru_trkphi,
0505       tru_trketa,
0506       tru_trkpt,
0507       1.,
0508       1.,
0509       1.,
0510       1.,
0511       1.,
0512       1.,
0513       1.,
0514       1.
0515     };
0516
0517     // fill truth histograms
0518     FillHistogram1D(content, Type::Truth, m_vecTreeHists1D);
0519     FillHistogram2D(content, Type::Truth, Comp::VsTruthPt, tru_trkpt,    m_vecTreeHists2D);
0520     FillHistogram2D(content, Type::Truth, Comp::VsNumTpc,  tru_nclustpc, m_vecTreeHists2D);
0521   }  // end truth particle loop
0522
0523   // announce next entry loop
0524   cout << "        Finished truth particle loop.\n"
0525        << "        Beginning reconstructed track loop: " << nRecoEntries << " to process"
0526        << endl;
0527
0528   // loop over reco tracks
0529   //   - TODO identify matched truth particles
0530   int64_t nRecoBytes = 0;
0531   for (int64_t iRecoEntry = 0; iRecoEntry < nRecoEntries; iRecoEntry++) {
0532
0533     // grab reco track entry
0534     const int64_t recoBytes = m_tTreeReco -> GetEntry(iRecoEntry);
0535     if (recoBytes < 0) {
0536       cerr << "PANIC: issue with entry " << iRecoEntry << "! Aborting loop!\n" << endl;
0537       break;
0538     } else {
0539       nRecoBytes += recoBytes;
0540     }
0541
0542     const int64_t iRecoProg = iRecoEntry + 1;
0543     if (iRecoProg == nRecoEntries) {
0544       cout << "          Processing entry " << iRecoProg << "/" << nRecoEntries << "..." << endl;
0545     } else {
0546       cout << "          Processing entry " << iRecoProg << "/" << nRecoEntries << "...\r" << flush;
0547     }
0548
0549     // select only tracks matched to truth particle
0550     //   - FIXME add cuts
0551     if (!rec_is_matched || !rec_is_Svtrack) continue;
0552
0553     // bundle histogram values to fill
0554     //   - FIXME actually calculate pt fraction
0555     //   - FIXME include quality, errors, and resolutions
0556     const STrackMatcherComparatorHistContent content {
0557       (double) rec_nclus,
0558       (double) rec_nclusintt,
0559       (double) rec_nclusmvtx,
0560       (double) rec_nclustpc,
0561       rec_matchrat,
0562       rec_matchrat_intt,
0563       rec_matchrat_mvtx,
0564       rec_matchrat_tpc,
0565       rec_trkphi,
0566       rec_trketa,
0567       rec_trkpt,
0568       1.,
0569       1.,
0570       1.,
0571       1.,
0572       1.,
0573       1.,
0574       1.,
0575       1.
0576     };
0577
0578     // fill all matched reco histograms
0579     //   - FIXME use actual truth pt & ntpc of matched particle
0580     FillHistogram1D(content, Type::Track, m_vecTreeHists1D);
0581     FillHistogram2D(content, Type::Track, Comp::VsTruthPt, rec_trkpt,    m_vecTreeHists2D);
0582     FillHistogram2D(content, Type::Track, Comp::VsTruthPt, rec_nclustpc, m_vecTreeHists2D);
0583
0584     // fill weird and normal matched reco 1D histograms
0585     //   - FIXME actually check if is a weird track
0586     const bool isWeirdTrack = true;
0587     if (isWeirdTrack) {
0588       FillHistogram1D(content, Type::Weird, m_vecTreeHists1D);
0589       FillHistogram2D(content, Type::Weird, Comp::VsTruthPt, rec_trkpt,    m_vecTreeHists2D);
0590       FillHistogram2D(content, Type::Weird, Comp::VsNumTpc,  rec_nclustpc, m_vecTreeHists2D);
0591     } else {
0592       FillHistogram1D(content, Type::Normal, m_vecTreeHists1D);
0593       FillHistogram2D(content, Type::Normal, Comp::VsTruthPt, rec_trkpt,    m_vecTreeHists2D);
0594       FillHistogram2D(content, Type::Normal, Comp::VsNumTpc,  rec_nclustpc, m_vecTreeHists2D);
0595     }
0596   }  // end reco track loop
0597
0598   // announce method end and return
0599   cout << "        Finished reconstructed track loop.\n"
0600        << "      Finished getting histograms from new matcher cluster tree."
0601        << endl;
0602   return;
0603
0604 }  // end 'GetNewTreeHists()'
0605
0606
0607
0608 void STrackMatcherComparator::GetNewTupleHists() {
0609
0610   // announce start of method
0611   cout << "      Grabbing new matcher tuple histograms:" << endl;
0612
0613   // declare input leaves
0614   float tru_evtid;
0615   float tru_trkid;
0616   float tru_pt;
0617   float tru_eta;
0618   float tru_phi;
0619   float tru_nmvtxclust_trkmatcher;
0620   float tru_ninttclust_trkmatcher;
0621   float tru_ntpclust_trkmatcher;
0622   float tru_nmvtxclust_manual;
0623   float tru_ninttclust_manual;
0624   float tru_ntpcclust_manual;
0625   float tru_nmvtxlayer_;
0626   float tru_ninttlayer;
0627   float tru_ntpclayer;
0628   float tru_deltapt;
0629   float tru_quality;
0630   float tru_dcaxy;
0631   float tru_dcaz;
0632   float tru_sigmadcaxy;
0633   float tru_sigmadacz;
0634   float tru_vx;
0635   float tru_vy;
0636   float tru_vz;
0637   float tru_gtrkid;
0638   float tru_gpt;
0639   float tru_geta;
0640   float tru_gphi;
0641   float tru_gnmvtxclust_trkmatcher;
0642   float tru_gninttclust_trkmatcher;
0643   float tru_gntpclust_trkmatcher;
0644   float tru_gnmvtxclust_manual;
0645   float tru_gninttclust_manual;
0646   float tru_gntpcclust_manual;
0647   float tru_gnmvtxlayer;
0648   float tru_gninttlayer;
0649   float tru_gntpclayer;
0650   float tru_gdeltapt;
0651   float tru_gquality;
0652   float tru_gdcaxy;
0653   float tru_gdcaz;
0654   float tru_gsigmadcaxy;
0655   float tru_gsigmadacz;
0656   float tru_gvx;
0657   float tru_gvy;
0658   float tru_gvz;
0659   float tru_fracnmvtxmatched;
0660   float tru_fracninttmatched;
0661   float tru_fracntpcmatched;
0662
0663   float rec_evtid;
0664   float rec_trkid;
0665   float rec_pt;
0666   float rec_eta;
0667   float rec_phi;
0668   float rec_nmvtxclust_trkmatcher;
0669   float rec_ninttclust_trkmatcher;
0670   float rec_ntpclust_trkmatcher;
0671   float rec_nmvtxclust_manual;
0672   float rec_ninttclust_manual;
0673   float rec_ntpcclust_manual;
0674   float rec_nmvtxlayer_;
0675   float rec_ninttlayer;
0676   float rec_ntpclayer;
0677   float rec_deltapt;
0678   float rec_quality;
0679   float rec_dcaxy;
0680   float rec_dcaz;
0681   float rec_sigmadcaxy;
0682   float rec_sigmadacz;
0683   float rec_vx;
0684   float rec_vy;
0685   float rec_vz;
0686   float rec_gtrkid;
0687   float rec_gpt;
0688   float rec_geta;
0689   float rec_gphi;
0690   float rec_gnmvtxclust_trkmatcher;
0691   float rec_gninttclust_trkmatcher;
0692   float rec_gntpclust_trkmatcher;
0693   float rec_gnmvtxclust_manual;
0694   float rec_gninttclust_manual;
0695   float rec_gntpcclust_manual;
0696   float rec_gnmvtxlayer;
0697   float rec_gninttlayer;
0698   float rec_gntpclayer;
0699   float rec_gdeltapt;
0700   float rec_gquality;
0701   float rec_gdcaxy;
0702   float rec_gdcaz;
0703   float rec_gsigmadcaxy;
0704   float rec_gsigmadacz;
0705   float rec_gvx;
0706   float rec_gvy;
0707   float rec_gvz;
0708   float rec_fracnmvtxmatched;
0709   float rec_fracninttmatched;
0710   float rec_fracntpcmatched;
0711
0712   // set branch addresses
0713   m_ntTupleTrue -> SetBranchAddress("evtid",                 &tru_evtid);
0714   m_ntTupleTrue -> SetBranchAddress("trkid",                 &tru_trkid);
0715   m_ntTupleTrue -> SetBranchAddress("pt",                    &tru_pt);
0716   m_ntTupleTrue -> SetBranchAddress("eta",                   &tru_eta);
0717   m_ntTupleTrue -> SetBranchAddress("phi",                   &tru_phi);
0718   m_ntTupleTrue -> SetBranchAddress("nmvtxclust_trkmatcher", &tru_nmvtxclust_trkmatcher);
0719   m_ntTupleTrue -> SetBranchAddress("ninttclust_trkmatcher", &tru_ninttclust_trkmatcher);
0720   m_ntTupleTrue -> SetBranchAddress("ntpclust_trkmatcher",   &tru_ntpclust_trkmatcher);
0721   m_ntTupleTrue -> SetBranchAddress("nmvtxclust_manual",     &tru_nmvtxclust_manual);
0722   m_ntTupleTrue -> SetBranchAddress("ninttclust_manual",     &tru_ninttclust_manual);
0723   m_ntTupleTrue -> SetBranchAddress("ntpcclust_manual",      &tru_ntpcclust_manual);
0724   m_ntTupleTrue -> SetBranchAddress("nmvtxlayer_",           &tru_nmvtxlayer_);
0725   m_ntTupleTrue -> SetBranchAddress("ninttlayer",            &tru_ninttlayer);
0726   m_ntTupleTrue -> SetBranchAddress("ntpclayer",             &tru_ntpclayer);
0727   m_ntTupleTrue -> SetBranchAddress("deltapt",               &tru_deltapt);
0728   m_ntTupleTrue -> SetBranchAddress("quality",               &tru_quality);
0729   m_ntTupleTrue -> SetBranchAddress("dcaxy",                 &tru_dcaxy);
0730   m_ntTupleTrue -> SetBranchAddress("dcaz",                  &tru_dcaz);
0731   m_ntTupleTrue -> SetBranchAddress("sigmadcaxy",            &tru_sigmadcaxy);
0732   m_ntTupleTrue -> SetBranchAddress("sigmadacz",             &tru_sigmadacz);
0733   m_ntTupleTrue -> SetBranchAddress("vx",                    &tru_vx);
0734   m_ntTupleTrue -> SetBranchAddress("vy",                    &tru_vy);
0735   m_ntTupleTrue -> SetBranchAddress("vz",                    &tru_vz);
0736   m_ntTupleTrue -> SetBranchAddress("gtrkid",                &tru_gtrkid);
0737   m_ntTupleTrue -> SetBranchAddress("gpt",                   &tru_gpt);
0738   m_ntTupleTrue -> SetBranchAddress("geta",                  &tru_geta);
0739   m_ntTupleTrue -> SetBranchAddress("gphi",                  &tru_gphi);
0740   m_ntTupleTrue -> SetBranchAddress("gnmvtxclust_trkmatcher",&tru_gnmvtxclust_trkmatcher);
0741   m_ntTupleTrue -> SetBranchAddress("gninttclust_trkmatcher",&tru_gninttclust_trkmatcher);
0742   m_ntTupleTrue -> SetBranchAddress("gntpclust_trkmatcher",  &tru_gntpclust_trkmatcher);
0743   m_ntTupleTrue -> SetBranchAddress("gnmvtxclust_manual",    &tru_gnmvtxclust_manual);
0744   m_ntTupleTrue -> SetBranchAddress("gninttclust_manual",    &tru_gninttclust_manual);
0745   m_ntTupleTrue -> SetBranchAddress("gntpcclust_manual",     &tru_gntpcclust_manual);
0746   m_ntTupleTrue -> SetBranchAddress("gnmvtxlayer",           &tru_gnmvtxlayer);
0747   m_ntTupleTrue -> SetBranchAddress("gninttlayer",           &tru_gninttlayer);
0748   m_ntTupleTrue -> SetBranchAddress("gntpclayer",            &tru_gntpclayer);
0749   m_ntTupleTrue -> SetBranchAddress("gdeltapt",              &tru_gdeltapt);
0750   m_ntTupleTrue -> SetBranchAddress("gquality",              &tru_gquality);
0751   m_ntTupleTrue -> SetBranchAddress("gdcaxy",                &tru_gdcaxy);
0752   m_ntTupleTrue -> SetBranchAddress("gdcaz",                 &tru_gdcaz);
0753   m_ntTupleTrue -> SetBranchAddress("gsigmadcaxy",           &tru_gsigmadcaxy);
0754   m_ntTupleTrue -> SetBranchAddress("gsigmadacz",            &tru_gsigmadacz);
0755   m_ntTupleTrue -> SetBranchAddress("gvx",                   &tru_gvx);
0756   m_ntTupleTrue -> SetBranchAddress("gvy",                   &tru_gvy);
0757   m_ntTupleTrue -> SetBranchAddress("gvz",                   &tru_gvz);
0758   m_ntTupleTrue -> SetBranchAddress("fracnmvtxmatched",      &tru_fracnmvtxmatched);
0759   m_ntTupleTrue -> SetBranchAddress("fracninttmatched",      &tru_fracninttmatched);
0760   m_ntTupleTrue -> SetBranchAddress("fracntpcmatched",       &tru_fracntpcmatched);
0761
0762   m_ntTupleReco -> SetBranchAddress("evtid",                  &rec_evtid);
0763   m_ntTupleReco -> SetBranchAddress("trkid",                  &rec_trkid);
0764   m_ntTupleReco -> SetBranchAddress("pt",                     &rec_pt);
0765   m_ntTupleReco -> SetBranchAddress("eta",                    &rec_eta);
0766   m_ntTupleReco -> SetBranchAddress("phi",                    &rec_phi);
0767   m_ntTupleReco -> SetBranchAddress("nmvtxclust_trkmatcher",  &rec_nmvtxclust_trkmatcher);
0768   m_ntTupleReco -> SetBranchAddress("ninttclust_trkmatcher",  &rec_ninttclust_trkmatcher);
0769   m_ntTupleReco -> SetBranchAddress("ntpclust_trkmatcher",    &rec_ntpclust_trkmatcher);
0770   m_ntTupleReco -> SetBranchAddress("nmvtxclust_manual",      &rec_nmvtxclust_manual);
0771   m_ntTupleReco -> SetBranchAddress("ninttclust_manual",      &rec_ninttclust_manual);
0772   m_ntTupleReco -> SetBranchAddress("ntpcclust_manual",       &rec_ntpcclust_manual);
0773   m_ntTupleReco -> SetBranchAddress("nmvtxlayer_",            &rec_nmvtxlayer_);
0774   m_ntTupleReco -> SetBranchAddress("ninttlayer",             &rec_ninttlayer);
0775   m_ntTupleReco -> SetBranchAddress("ntpclayer",              &rec_ntpclayer);
0776   m_ntTupleReco -> SetBranchAddress("deltapt",                &rec_deltapt);
0777   m_ntTupleReco -> SetBranchAddress("quality",                &rec_quality);
0778   m_ntTupleReco -> SetBranchAddress("dcaxy",                  &rec_dcaxy);
0779   m_ntTupleReco -> SetBranchAddress("dcaz",                   &rec_dcaz);
0780   m_ntTupleReco -> SetBranchAddress("sigmadcaxy",             &rec_sigmadcaxy);
0781   m_ntTupleReco -> SetBranchAddress("sigmadacz",              &rec_sigmadacz);
0782   m_ntTupleReco -> SetBranchAddress("vx",                     &rec_vx);
0783   m_ntTupleReco -> SetBranchAddress("vy",                     &rec_vy);
0784   m_ntTupleReco -> SetBranchAddress("vz",                     &rec_vz);
0785   m_ntTupleReco -> SetBranchAddress("gtrkid",                 &rec_gtrkid);
0786   m_ntTupleReco -> SetBranchAddress("gpt",                    &rec_gpt);
0787   m_ntTupleReco -> SetBranchAddress("geta",                   &rec_geta);
0788   m_ntTupleReco -> SetBranchAddress("gphi",                   &rec_gphi);
0789   m_ntTupleReco -> SetBranchAddress("gnmvtxclust_trkmatcher", &rec_gnmvtxclust_trkmatcher);
0790   m_ntTupleReco -> SetBranchAddress("gninttclust_trkmatcher", &rec_gninttclust_trkmatcher);
0791   m_ntTupleReco -> SetBranchAddress("gntpclust_trkmatcher",   &rec_gntpclust_trkmatcher);
0792   m_ntTupleReco -> SetBranchAddress("gnmvtxclust_manual",     &rec_gnmvtxclust_manual);
0793   m_ntTupleReco -> SetBranchAddress("gninttclust_manual",     &rec_gninttclust_manual);
0794   m_ntTupleReco -> SetBranchAddress("gntpcclust_manual",      &rec_gntpcclust_manual);
0795   m_ntTupleReco -> SetBranchAddress("gnmvtxlayer",            &rec_gnmvtxlayer);
0796   m_ntTupleReco -> SetBranchAddress("gninttlayer",            &rec_gninttlayer);
0797   m_ntTupleReco -> SetBranchAddress("gntpclayer",             &rec_gntpclayer);
0798   m_ntTupleReco -> SetBranchAddress("gdeltapt",               &rec_gdeltapt);
0799   m_ntTupleReco -> SetBranchAddress("gquality",               &rec_gquality);
0800   m_ntTupleReco -> SetBranchAddress("gdcaxy",                 &rec_gdcaxy);
0801   m_ntTupleReco -> SetBranchAddress("gdcaz",                  &rec_gdcaz);
0802   m_ntTupleReco -> SetBranchAddress("gsigmadcaxy",            &rec_gsigmadcaxy);
0803   m_ntTupleReco -> SetBranchAddress("gsigmadacz",             &rec_gsigmadacz);
0804   m_ntTupleReco -> SetBranchAddress("gvx",                    &rec_gvx);
0805   m_ntTupleReco -> SetBranchAddress("gvy",                    &rec_gvy);
0806   m_ntTupleReco -> SetBranchAddress("gvz",                    &rec_gvz);
0807   m_ntTupleReco -> SetBranchAddress("fracnmvtxmatched",       &rec_fracnmvtxmatched);
0808   m_ntTupleReco -> SetBranchAddress("fracninttmatched",       &rec_fracninttmatched);
0809   m_ntTupleReco -> SetBranchAddress("fracntpcmatched",        &rec_fracntpcmatched);
0810   cout << "        Set input branches." << endl;
0811
0812   // grab no. of entries
0813   const int64_t nTrueEntries = m_ntTupleTrue -> GetEntries();
0814   const int64_t nRecoEntries = m_ntTupleReco -> GetEntries();
0815   cout << "        Beginning truth particle loop: " << nTrueEntries << " to process" << endl;
0816
0817   // loop over truth particles
0818   int64_t nTrueBytes = 0;
0819   for (int64_t iTrueEntry = 0; iTrueEntry < nTrueEntries; iTrueEntry++) {
0820
0821     // grab truth particle entry
0822     const int64_t trueBytes = m_ntTupleTrue -> GetEntry(iTrueEntry);
0823     if (trueBytes < 0) {
0824       cerr << "PANIC: issue with entry " << iTrueEntry << "! Aborting loop!\n" << endl;
0825       break;
0826     } else {
0827       nTrueBytes += trueBytes;
0828     }
0829
0830     const int64_t iTrueProg = iTrueEntry + 1;
0831     if (iTrueProg == nTrueEntries) {
0832       cout << "          Processing entry " << iTrueProg << "/" << nTrueEntries << "..." << endl;
0833     } else {
0834       cout << "          Processing entry " << iTrueProg << "/" << nTrueEntries << "...\r" << flush;
0835     }
0836
0837     // check if near sector
0838     bool isNearSector = false;
0839     if (m_config.doPhiCut) {
0840       isNearSector = IsNearSector(tru_gphi);
0841     }
0842
0843     // apply cuts
0844     const bool isInZVtxCut = ((tru_gvz >= m_config.zVtxRange.first) && (tru_gvz <= m_config.zVtxRange.second));
0845     if (m_config.doZVtxCut && !isInZVtxCut) continue;
0846     if (m_config.doPhiCut  && isNearSector) continue;
0847
0848     // run calculations
0849     const float tru_gnclust = tru_gnmvtxclust_trkmatcher + tru_gninttclust_trkmatcher + tru_gntpclust_trkmatcher;
0850
0851     // bundle histogram values to fill
0852     STrackMatcherComparatorHistContent content {
0853       tru_gnclust,
0854       tru_gninttclust_trkmatcher,
0855       tru_gnmvtxclust_trkmatcher,
0856       tru_gntpclust_trkmatcher,
0857       1.,
0858       1.,
0859       1.,
0860       1.,
0861       tru_gphi,
0862       tru_geta,
0863       tru_gpt,
0864       1.,
0865       1.,
0866       1.,
0867       1.,
0868       1.,
0869       1.,
0870       1.,
0871       1.
0872     };
0873
0874     // fill truth histograms
0875     FillHistogram1D(content, Type::Truth, m_vecTupleHists1D);
0876     FillHistogram2D(content, Type::Truth, Comp::VsTruthPt, tru_gpt,                 m_vecTupleHists2D);
0877     FillHistogram2D(content, Type::Truth, Comp::VsNumTpc,  tru_ntpclust_trkmatcher, m_vecTupleHists2D);
0878   }  // end truth particle loop
0879
0880   // announce next entry loop
0881   cout << "        Finished truth particle loop.\n"
0882        << "        Beginning reconstructed track loop: " << nRecoEntries << " to process"
0883        << endl;
0884
0885   // loop over reco tracks
0886   int64_t nRecoBytes = 0;
0887   for (int64_t iRecoEntry = 0; iRecoEntry < nRecoEntries; iRecoEntry++) {
0888
0889     // grab reco track entry
0890     const int64_t recoBytes = m_ntTupleReco -> GetEntry(iRecoEntry);
0891     if (recoBytes < 0) {
0892       cerr << "PANIC: issue with entry " << iRecoEntry << "! Aborting loop!\n" << endl;
0893       break;
0894     } else {
0895       nRecoBytes += recoBytes;
0896     }
0897
0898     const int64_t iRecoProg = iRecoEntry + 1;
0899     if (iRecoProg == nRecoEntries) {
0900       cout << "          Processing entry " << iRecoProg << "/" << nRecoEntries << "..." << endl;
0901     } else {
0902       cout << "          Processing entry " << iRecoProg << "/" << nRecoEntries << "...\r" << flush;
0903     }
0904
0905     // check if near sector
0906     bool isNearSector = false;
0907     if (m_config.doPhiCut) {
0908       isNearSector = IsNearSector(rec_phi);
0909     }
0910
0911     // apply cuts
0912     const bool isInZVtxCut = ((rec_vz >= m_config.zVtxRange.first) && (rec_vz <= m_config.zVtxRange.second));
0913     if (m_config.doZVtxCut && !isInZVtxCut) continue;
0914     if (m_config.doPhiCut  && isNearSector) continue;
0915
0916     // run calculations
0917     //   - FIXME add other errors
0918     const double rec_nclus  = rec_ninttclust_trkmatcher + rec_nmvtxclust_trkmatcher + rec_ntpclust_trkmatcher;
0919     const double rec_gnclus = rec_gninttclust_trkmatcher + rec_gnmvtxclust_trkmatcher + rec_gntpcclust_manual;
0920     const double rec_rnclus = rec_nclus / rec_gnclus;
0921     const double rec_rintt  = rec_ninttclust_trkmatcher / rec_gninttclust_trkmatcher;
0922     const double rec_rmaps  = rec_nmvtxclust_trkmatcher / rec_gnmvtxclust_trkmatcher;
0923     const double rec_rtpc   = rec_ntpclust_trkmatcher / rec_gntpcclust_manual;
0924     const double rec_ptfrac = rec_pt / rec_gpt;
0925     const double rec_ptper  = rec_deltapt / rec_pt;
0926     const double rec_etaper = 1.;
0927     const double rec_phiper = 1.;
0928     const double rec_ptres  = abs(rec_pt - rec_gpt) / rec_gpt;
0929     const double rec_etares = abs(rec_eta - rec_geta) / rec_geta;
0930     const double rec_phires = abs(rec_phi - rec_gphi) / rec_gphi;
0931
0932     // bundle histogram values to fill
0933     const STrackMatcherComparatorHistContent content {
0934       rec_nclus,
0935       rec_ninttclust_trkmatcher,
0936       rec_nmvtxclust_trkmatcher,
0937       rec_ntpclust_trkmatcher,
0938       rec_rnclus,
0939       rec_rintt,
0940       rec_rmaps,
0941       rec_rtpc,
0942       rec_phi,
0943       rec_eta,
0944       rec_pt,
0945       rec_ptfrac,
0946       rec_quality,
0947       rec_ptper,
0948       rec_etaper,
0949       rec_phiper,
0950       rec_ptres,
0951       rec_etares,
0952       rec_phires
0953     };
0954
0955     // fill all matched reco histograms
0956     FillHistogram1D(content, Type::Track, m_vecTupleHists1D);
0957     FillHistogram2D(content, Type::Track, Comp::VsTruthPt, rec_gpt,                 m_vecTupleHists2D);
0958     FillHistogram2D(content, Type::Track, Comp::VsNumTpc,  rec_ntpclust_trkmatcher, m_vecTupleHists2D);
0959
0960     // fill weird and normal matched reco 1D histograms
0961     const bool isNormalTrack = ((rec_ptfrac >= m_config.oddPtFrac.first) && (rec_ptfrac <= m_config.oddPtFrac.second));
0962     if (isNormalTrack) {
0963       FillHistogram1D(content, Type::Normal, m_vecTupleHists1D);
0964       FillHistogram2D(content, Type::Normal, Comp::VsTruthPt, rec_gpt,                 m_vecTupleHists2D);
0965       FillHistogram2D(content, Type::Normal, Comp::VsNumTpc,  rec_ntpclust_trkmatcher, m_vecTupleHists2D);
0966     } else {
0967       FillHistogram1D(content, Type::Weird, m_vecTupleHists1D);
0968       FillHistogram2D(content, Type::Weird, Comp::VsTruthPt, rec_gpt,                 m_vecTupleHists2D);
0969       FillHistogram2D(content, Type::Weird, Comp::VsNumTpc,  rec_ntpclust_trkmatcher, m_vecTupleHists2D);
0970     }
0971   }  // end reco track loop
0972
0973   // announce method end and return
0974   cout << "        Finished reconstructed track loop.\n"
0975        << "      Finished getting histograms from new matcher track tuple."
0976        << endl;
0977   return;
0978
0979 }  // end 'GetNewTupleHists()'
0980
0981
0982
0983 void STrackMatcherComparator::GetOldTupleHists() {
0984
0985   // announce start of method
0986   cout << "      Grabbing old matcher tuple histograms:" << endl;
0987
0988   // declare input leaves
0989   float tru_event;
0990   float tru_seed;
0991   float tru_gntracks;
0992   float tru_gnchghad;
0993   float tru_gtrackID;
0994   float tru_gflavor;
0995   float tru_gnhits;
0996   float tru_gnmaps;
0997   float tru_gnintt;
0998   float tru_gnmms;
0999   float tru_gnintt1;
1000   float tru_gnintt2;
1001   float tru_gnintt3;
1002   float tru_gnintt4;
1003   float tru_gnintt5;
1004   float tru_gnintt6;
1005   float tru_gnintt7;
1006   float tru_gnintt8;
1007   float tru_gntpc;
1008   float tru_gnlmaps;
1009   float tru_gnlintt;
1010   float tru_gnltpc;
1011   float tru_gnlmms;
1012   float tru_gpx;
1013   float tru_gpy;
1014   float tru_gpz;
1015   float tru_gpt;
1016   float tru_geta;
1017   float tru_gphi;
1018   float tru_gvx;
1019   float tru_gvy;
1020   float tru_gvz;
1021   float tru_gvt;
1022   float tru_gfpx;
1023   float tru_gfpy;
1024   float tru_gfpz;
1025   float tru_gfx;
1026   float tru_gfy;
1027   float tru_gfz;
1028   float tru_gembed;
1029   float tru_gprimary;
1030   float tru_trackID;
1031   float tru_px;
1032   float tru_py;
1033   float tru_pz;
1034   float tru_pt;
1035   float tru_eta;
1036   float tru_phi;
1037   float tru_deltapt;
1038   float tru_deltaeta;
1039   float tru_deltaphi;
1040   float tru_siqr;
1041   float tru_siphi;
1042   float tru_sithe;
1043   float tru_six0;
1044   float tru_siy0;
1045   float tru_tpqr;
1046   float tru_tpphi;
1047   float tru_tpthe;
1048   float tru_tpx0;
1049   float tru_tpy0;
1050   float tru_charge;
1051   float tru_quality;
1052   float tru_chisq;
1053   float tru_ndf;
1054   float tru_nhits;
1055   float tru_layers;
1056   float tru_nmaps;
1057   float tru_nintt;
1058   float tru_ntpc;
1059   float tru_nmms;
1060   float tru_ntpc1;
1061   float tru_ntpc11;
1062   float tru_ntpc2;
1063   float tru_ntpc3;
1064   float tru_nlmaps;
1065   float tru_nlintt;
1066   float tru_nltpc;
1067   float tru_nlmms;
1068   float tru_vertexID;
1069   float tru_vx;
1070   float tru_vy;
1071   float tru_vz;
1072   float tru_dca2d;
1073   float tru_dca2dsigma;
1074   float tru_dca3dxy;
1075   float tru_dca3dxysigma;
1076   float tru_dca3dz;
1077   float tru_dca3dzsigma;
1078   float tru_pcax;
1079   float tru_pcay;
1080   float tru_pcaz;
1081   float tru_nfromtruth;
1082   float tru_nwrong;
1083   float tru_ntrumaps;
1084   float tru_nwrongmaps;
1085   float tru_ntruintt;
1086   float tru_nwrongintt;
1087   float tru_ntrutpc;
1088   float tru_nwrongtpc;
1089   float tru_ntrumms;
1090   float tru_nwrongmms;
1091   float tru_ntrutpc1;
1092   float tru_nwrongtpc1;
1093   float tru_ntrutpc11;
1094   float tru_nwrongtpc11;
1095   float tru_ntrutpc2;
1096   float tru_nwrongtpc2;
1097   float tru_ntrutpc3;
1098   float tru_nwrongtpc3;
1099   float tru_layersfromtruth;
1100   float tru_npedge;
1101   float tru_nredge;
1102   float tru_nbig;
1103   float tru_novlp;
1104   float tru_merr;
1105   float tru_msize;
1106   float tru_nhittpcall;
1107   float tru_nhittpcin;
1108   float tru_nhittpcmid;
1109   float tru_nhittpcout;
1110   float tru_nclusall;
1111   float tru_nclustpc;
1112   float tru_nclusintt;
1113   float tru_nclusmaps;
1114   float tru_nclusmms;
1115
1116   float rec_event;
1117   float rec_seed;
1118   float rec_trackID;
1119   float rec_crossing;
1120   float rec_px;
1121   float rec_py;
1122   float rec_pz;
1123   float rec_pt;
1124   float rec_eta;
1125   float rec_phi;
1126   float rec_deltapt;
1127   float rec_deltaeta;
1128   float rec_deltaphi;
1129   float rec_siqr;
1130   float rec_siphi;
1131   float rec_sithe;
1132   float rec_six0;
1133   float rec_siy0;
1134   float rec_tpqr;
1135   float rec_tpphi;
1136   float rec_tpthe;
1137   float rec_tpx0;
1138   float rec_tpy0;
1139   float rec_charge;
1140   float rec_quality;
1141   float rec_chisq;
1142   float rec_ndf;
1143   float rec_nhits;
1144   float rec_nmaps;
1145   float rec_nintt;
1146   float rec_ntpc;
1147   float rec_nmms;
1148   float rec_ntpc1;
1149   float rec_ntpc11;
1150   float rec_ntpc2;
1151   float rec_ntpc3;
1152   float rec_nlmaps;
1153   float rec_nlintt;
1154   float rec_nltpc;
1155   float rec_nlmms;
1156   float rec_layers;
1157   float rec_vertexID;
1158   float rec_vx;
1159   float rec_vy;
1160   float rec_vz;
1161   float rec_dca2d;
1162   float rec_dca2dsigma;
1163   float rec_dca3dxy;
1164   float rec_dca3dxysigma;
1165   float rec_dca3dz;
1166   float rec_dca3dzsigma;
1167   float rec_pcax;
1168   float rec_pcay;
1169   float rec_pcaz;
1170   float rec_gtrackID;
1171   float rec_gflavor;
1172   float rec_gnhits;
1173   float rec_gnmaps;
1174   float rec_gnintt;
1175   float rec_gntpc;
1176   float rec_gnmms;
1177   float rec_gnlmaps;
1178   float rec_gnlintt;
1179   float rec_gnltpc;
1180   float rec_gnlmms;
1181   float rec_gpx;
1182   float rec_gpy;
1183   float rec_gpz;
1184   float rec_gpt;
1185   float rec_geta;
1186   float rec_gphi;
1187   float rec_gvx;
1188   float rec_gvy;
1189   float rec_gvz;
1190   float rec_gvt;
1191   float rec_gfpx;
1192   float rec_gfpy;
1193   float rec_gfpz;
1194   float rec_gfx;
1195   float rec_gfy;
1196   float rec_gfz;
1197   float rec_gembed;
1198   float rec_gprimary;
1199   float rec_nfromtruth;
1200   float rec_nwrong;
1201   float rec_ntrumaps;
1202   float rec_nwrongmaps;
1203   float rec_ntruintt;
1204   float rec_nwrongintt;
1205   float rec_ntrutpc;
1206   float rec_nwrongtpc;
1207   float rec_ntrumms;
1208   float rec_nwrongmms;
1209   float rec_ntrutpc1;
1210   float rec_nwrongtpc1;
1211   float rec_ntrutpc11;
1212   float rec_nwrongtpc11;
1213   float rec_ntrutpc2;
1214   float rec_nwrongtpc2;
1215   float rec_ntrutpc3;
1216   float rec_nwrongtpc3;
1217   float rec_layersfromtruth;
1218   float rec_npedge;
1219   float rec_nredge;
1220   float rec_nbig;
1221   float rec_novlp;
1222   float rec_merr;
1223   float rec_msize;
1224   float rec_nhittpcall;
1225   float rec_nhittpcin;
1226   float rec_nhittpcmid;
1227   float rec_nhittpcout;
1228   float rec_nclusall;
1229   float rec_nclustpc;
1230   float rec_nclusintt;
1231   float rec_nclusmaps;
1232   float rec_nclusmms;
1233
1234   // Set branch addresses.
1235   m_ntOldTrue -> SetBranchAddress("event",           &tru_event);
1236   m_ntOldTrue -> SetBranchAddress("seed",            &tru_seed);
1237   m_ntOldTrue -> SetBranchAddress("gntracks",        &tru_gntracks);
1238   m_ntOldTrue -> SetBranchAddress("gnchghad",        &tru_gnchghad);
1239   m_ntOldTrue -> SetBranchAddress("gtrackID",        &tru_gtrackID);
1240   m_ntOldTrue -> SetBranchAddress("gflavor",         &tru_gflavor);
1241   m_ntOldTrue -> SetBranchAddress("gnhits",          &tru_gnhits);
1242   m_ntOldTrue -> SetBranchAddress("gnmaps",          &tru_gnmaps);
1243   m_ntOldTrue -> SetBranchAddress("gnintt",          &tru_gnintt);
1244   m_ntOldTrue -> SetBranchAddress("gnmms",           &tru_gnmms);
1245   m_ntOldTrue -> SetBranchAddress("gnintt1",         &tru_gnintt1);
1246   m_ntOldTrue -> SetBranchAddress("gnintt2",         &tru_gnintt2);
1247   m_ntOldTrue -> SetBranchAddress("gnintt3",         &tru_gnintt3);
1248   m_ntOldTrue -> SetBranchAddress("gnintt4",         &tru_gnintt4);
1249   m_ntOldTrue -> SetBranchAddress("gnintt5",         &tru_gnintt5);
1250   m_ntOldTrue -> SetBranchAddress("gnintt6",         &tru_gnintt6);
1251   m_ntOldTrue -> SetBranchAddress("gnintt7",         &tru_gnintt7);
1252   m_ntOldTrue -> SetBranchAddress("gnintt8",         &tru_gnintt8);
1253   m_ntOldTrue -> SetBranchAddress("gntpc",           &tru_gntpc);
1254   m_ntOldTrue -> SetBranchAddress("gnlmaps",         &tru_gnlmaps);
1255   m_ntOldTrue -> SetBranchAddress("gnlintt",         &tru_gnlintt);
1256   m_ntOldTrue -> SetBranchAddress("gnltpc",          &tru_gnltpc);
1257   m_ntOldTrue -> SetBranchAddress("gnlmms",          &tru_gnlmms);
1258   m_ntOldTrue -> SetBranchAddress("gpx",             &tru_gpx);
1259   m_ntOldTrue -> SetBranchAddress("gpy",             &tru_gpy);
1260   m_ntOldTrue -> SetBranchAddress("gpz",             &tru_gpz);
1261   m_ntOldTrue -> SetBranchAddress("gpt",             &tru_gpt);
1262   m_ntOldTrue -> SetBranchAddress("geta",            &tru_geta);
1263   m_ntOldTrue -> SetBranchAddress("gphi",            &tru_gphi);
1264   m_ntOldTrue -> SetBranchAddress("gvx",             &tru_gvx);
1265   m_ntOldTrue -> SetBranchAddress("gvy",             &tru_gvy);
1266   m_ntOldTrue -> SetBranchAddress("gvz",             &tru_gvz);
1267   m_ntOldTrue -> SetBranchAddress("gvt",             &tru_gvt);
1268   m_ntOldTrue -> SetBranchAddress("gfpx",            &tru_gfpx);
1269   m_ntOldTrue -> SetBranchAddress("gfpy",            &tru_gfpy);
1270   m_ntOldTrue -> SetBranchAddress("gfpz",            &tru_gfpz);
1271   m_ntOldTrue -> SetBranchAddress("gfx",             &tru_gfx);
1272   m_ntOldTrue -> SetBranchAddress("gfy",             &tru_gfy);
1273   m_ntOldTrue -> SetBranchAddress("gfz",             &tru_gfz);
1274   m_ntOldTrue -> SetBranchAddress("gembed",          &tru_gembed);
1275   m_ntOldTrue -> SetBranchAddress("gprimary",        &tru_gprimary);
1276   m_ntOldTrue -> SetBranchAddress("trackID",         &tru_trackID);
1277   m_ntOldTrue -> SetBranchAddress("px",              &tru_px);
1278   m_ntOldTrue -> SetBranchAddress("py",              &tru_py);
1279   m_ntOldTrue -> SetBranchAddress("pz",              &tru_pz);
1280   m_ntOldTrue -> SetBranchAddress("pt",              &tru_pt);
1281   m_ntOldTrue -> SetBranchAddress("eta",             &tru_eta);
1282   m_ntOldTrue -> SetBranchAddress("phi",             &tru_phi);
1283   m_ntOldTrue -> SetBranchAddress("deltapt",         &tru_deltapt);
1284   m_ntOldTrue -> SetBranchAddress("deltaeta",        &tru_deltaeta);
1285   m_ntOldTrue -> SetBranchAddress("deltaphi",        &tru_deltaphi);
1286   m_ntOldTrue -> SetBranchAddress("siqr",            &tru_siqr);
1287   m_ntOldTrue -> SetBranchAddress("siphi",           &tru_siphi);
1288   m_ntOldTrue -> SetBranchAddress("sithe",           &tru_sithe);
1289   m_ntOldTrue -> SetBranchAddress("six0",            &tru_six0);
1290   m_ntOldTrue -> SetBranchAddress("siy0",            &tru_siy0);
1291   m_ntOldTrue -> SetBranchAddress("tpqr",            &tru_tpqr);
1292   m_ntOldTrue -> SetBranchAddress("tpphi",           &tru_tpphi);
1293   m_ntOldTrue -> SetBranchAddress("tpthe",           &tru_tpthe);
1294   m_ntOldTrue -> SetBranchAddress("tpx0",            &tru_tpx0);
1295   m_ntOldTrue -> SetBranchAddress("tpy0",            &tru_tpy0);
1296   m_ntOldTrue -> SetBranchAddress("charge",          &tru_charge);
1297   m_ntOldTrue -> SetBranchAddress("quality",         &tru_quality);
1298   m_ntOldTrue -> SetBranchAddress("chisq",           &tru_chisq);
1299   m_ntOldTrue -> SetBranchAddress("ndf",             &tru_ndf);
1300   m_ntOldTrue -> SetBranchAddress("nhits",           &tru_nhits);
1301   m_ntOldTrue -> SetBranchAddress("layers",          &tru_layers);
1302   m_ntOldTrue -> SetBranchAddress("nmaps",           &tru_nmaps);
1303   m_ntOldTrue -> SetBranchAddress("nintt",           &tru_nintt);
1304   m_ntOldTrue -> SetBranchAddress("ntpc",            &tru_ntpc);
1305   m_ntOldTrue -> SetBranchAddress("nmms",            &tru_nmms);
1306   m_ntOldTrue -> SetBranchAddress("ntpc1",           &tru_ntpc1);
1307   m_ntOldTrue -> SetBranchAddress("ntpc11",          &tru_ntpc11);
1308   m_ntOldTrue -> SetBranchAddress("ntpc2",           &tru_ntpc2);
1309   m_ntOldTrue -> SetBranchAddress("ntpc3",           &tru_ntpc3);
1310   m_ntOldTrue -> SetBranchAddress("nlmaps",          &tru_nlmaps);
1311   m_ntOldTrue -> SetBranchAddress("nlintt",          &tru_nlintt);
1312   m_ntOldTrue -> SetBranchAddress("nltpc",           &tru_nltpc);
1313   m_ntOldTrue -> SetBranchAddress("nlmms",           &tru_nlmms);
1314   m_ntOldTrue -> SetBranchAddress("vertexID",        &tru_vertexID);
1315   m_ntOldTrue -> SetBranchAddress("vx",              &tru_vx);
1316   m_ntOldTrue -> SetBranchAddress("vy",              &tru_vy);
1317   m_ntOldTrue -> SetBranchAddress("vz",              &tru_vz);
1318   m_ntOldTrue -> SetBranchAddress("dca2d",           &tru_dca2d);
1319   m_ntOldTrue -> SetBranchAddress("dca2dsigma",      &tru_dca2dsigma);
1320   m_ntOldTrue -> SetBranchAddress("dca3dxy",         &tru_dca3dxy);
1321   m_ntOldTrue -> SetBranchAddress("dca3dxysigma",    &tru_dca3dxysigma);
1322   m_ntOldTrue -> SetBranchAddress("dca3dz",          &tru_dca3dz);
1323   m_ntOldTrue -> SetBranchAddress("dca3dzsigma",     &tru_dca3dzsigma);
1324   m_ntOldTrue -> SetBranchAddress("pcax",            &tru_pcax);
1325   m_ntOldTrue -> SetBranchAddress("pcay",            &tru_pcay);
1326   m_ntOldTrue -> SetBranchAddress("pcaz",            &tru_pcaz);
1327   m_ntOldTrue -> SetBranchAddress("nfromtruth",      &tru_nfromtruth);
1328   m_ntOldTrue -> SetBranchAddress("nwrong",          &tru_nwrong);
1329   m_ntOldTrue -> SetBranchAddress("ntrumaps",        &tru_ntrumaps);
1330   m_ntOldTrue -> SetBranchAddress("nwrongmaps",      &tru_nwrongmaps);
1331   m_ntOldTrue -> SetBranchAddress("ntruintt",        &tru_ntruintt);
1332   m_ntOldTrue -> SetBranchAddress("nwrongintt",      &tru_nwrongintt);
1333   m_ntOldTrue -> SetBranchAddress("ntrutpc",         &tru_ntrutpc);
1334   m_ntOldTrue -> SetBranchAddress("nwrongtpc",       &tru_nwrongtpc);
1335   m_ntOldTrue -> SetBranchAddress("ntrumms",         &tru_ntrumms);
1336   m_ntOldTrue -> SetBranchAddress("nwrongmms",       &tru_nwrongmms);
1337   m_ntOldTrue -> SetBranchAddress("ntrutpc1",        &tru_ntrutpc1);
1338   m_ntOldTrue -> SetBranchAddress("nwrongtpc1",      &tru_nwrongtpc1);
1339   m_ntOldTrue -> SetBranchAddress("ntrutpc11",       &tru_ntrutpc11);
1340   m_ntOldTrue -> SetBranchAddress("nwrongtpc11",     &tru_nwrongtpc11);
1341   m_ntOldTrue -> SetBranchAddress("ntrutpc2",        &tru_ntrutpc2);
1342   m_ntOldTrue -> SetBranchAddress("nwrongtpc2",      &tru_nwrongtpc2);
1343   m_ntOldTrue -> SetBranchAddress("ntrutpc3",        &tru_ntrutpc3);
1344   m_ntOldTrue -> SetBranchAddress("nwrongtpc3",      &tru_nwrongtpc3);
1345   m_ntOldTrue -> SetBranchAddress("layersfromtruth", &tru_layersfromtruth);
1346   m_ntOldTrue -> SetBranchAddress("npedge",          &tru_npedge);
1347   m_ntOldTrue -> SetBranchAddress("nredge",          &tru_nredge);
1348   m_ntOldTrue -> SetBranchAddress("nbig",            &tru_nbig);
1349   m_ntOldTrue -> SetBranchAddress("novlp",           &tru_novlp);
1350   m_ntOldTrue -> SetBranchAddress("merr",            &tru_merr);
1351   m_ntOldTrue -> SetBranchAddress("msize",           &tru_msize);
1352   m_ntOldTrue -> SetBranchAddress("nhittpcall",      &tru_nhittpcall);
1353   m_ntOldTrue -> SetBranchAddress("nhittpcin",       &tru_nhittpcin);
1354   m_ntOldTrue -> SetBranchAddress("nhittpcmid",      &tru_nhittpcmid);
1355   m_ntOldTrue -> SetBranchAddress("nhittpcout",      &tru_nhittpcout);
1356   m_ntOldTrue -> SetBranchAddress("nclusall",        &tru_nclusall);
1357   m_ntOldTrue -> SetBranchAddress("nclustpc",        &tru_nclustpc);
1358   m_ntOldTrue -> SetBranchAddress("nclusintt",       &tru_nclusintt);
1359   m_ntOldTrue -> SetBranchAddress("nclusmaps",       &tru_nclusmaps);
1360   m_ntOldTrue -> SetBranchAddress("nclusmms",        &tru_nclusmms);
1361
1362   m_ntOldReco -> SetBranchAddress("event",           &rec_event);
1363   m_ntOldReco -> SetBranchAddress("seed",            &rec_seed);
1364   m_ntOldReco -> SetBranchAddress("trackID",         &rec_trackID);
1365   m_ntOldReco -> SetBranchAddress("crossing",        &rec_crossing);
1366   m_ntOldReco -> SetBranchAddress("px",              &rec_px);
1367   m_ntOldReco -> SetBranchAddress("py",              &rec_py);
1368   m_ntOldReco -> SetBranchAddress("pz",              &rec_pz);
1369   m_ntOldReco -> SetBranchAddress("pt",              &rec_pt);
1370   m_ntOldReco -> SetBranchAddress("eta",             &rec_eta);
1371   m_ntOldReco -> SetBranchAddress("phi",             &rec_phi);
1372   m_ntOldReco -> SetBranchAddress("deltapt",         &rec_deltapt);
1373   m_ntOldReco -> SetBranchAddress("deltaeta",        &rec_deltaeta);
1374   m_ntOldReco -> SetBranchAddress("deltaphi",        &rec_deltaphi);
1375   m_ntOldReco -> SetBranchAddress("siqr",            &rec_siqr);
1376   m_ntOldReco -> SetBranchAddress("siphi",           &rec_siphi);
1377   m_ntOldReco -> SetBranchAddress("sithe",           &rec_sithe);
1378   m_ntOldReco -> SetBranchAddress("six0",            &rec_six0);
1379   m_ntOldReco -> SetBranchAddress("siy0",            &rec_siy0);
1380   m_ntOldReco -> SetBranchAddress("tpqr",            &rec_tpqr);
1381   m_ntOldReco -> SetBranchAddress("tpphi",           &rec_tpphi);
1382   m_ntOldReco -> SetBranchAddress("tpthe",           &rec_tpthe);
1383   m_ntOldReco -> SetBranchAddress("tpx0",            &rec_tpx0);
1384   m_ntOldReco -> SetBranchAddress("tpy0",            &rec_tpy0);
1385   m_ntOldReco -> SetBranchAddress("charge",          &rec_charge);
1386   m_ntOldReco -> SetBranchAddress("quality",         &rec_quality);
1387   m_ntOldReco -> SetBranchAddress("chisq",           &rec_chisq);
1388   m_ntOldReco -> SetBranchAddress("ndf",             &rec_ndf);
1389   m_ntOldReco -> SetBranchAddress("nhits",           &rec_nhits);
1390   m_ntOldReco -> SetBranchAddress("nmaps",           &rec_nmaps);
1391   m_ntOldReco -> SetBranchAddress("nintt",           &rec_nintt);
1392   m_ntOldReco -> SetBranchAddress("ntpc",            &rec_ntpc);
1393   m_ntOldReco -> SetBranchAddress("nmms",            &rec_nmms);
1394   m_ntOldReco -> SetBranchAddress("ntpc1",           &rec_ntpc1);
1395   m_ntOldReco -> SetBranchAddress("ntpc11",          &rec_ntpc11);
1396   m_ntOldReco -> SetBranchAddress("ntpc2",           &rec_ntpc2);
1397   m_ntOldReco -> SetBranchAddress("ntpc3",           &rec_ntpc3);
1398   m_ntOldReco -> SetBranchAddress("nlmaps",          &rec_nlmaps);
1399   m_ntOldReco -> SetBranchAddress("nlintt",          &rec_nlintt);
1400   m_ntOldReco -> SetBranchAddress("nltpc",           &rec_nltpc);
1401   m_ntOldReco -> SetBranchAddress("nlmms",           &rec_nlmms);
1402   m_ntOldReco -> SetBranchAddress("layers",          &rec_layers);
1403   m_ntOldReco -> SetBranchAddress("vertexID",        &rec_vertexID);
1404   m_ntOldReco -> SetBranchAddress("vx",              &rec_vx);
1405   m_ntOldReco -> SetBranchAddress("vy",              &rec_vy);
1406   m_ntOldReco -> SetBranchAddress("vz",              &rec_vz);
1407   m_ntOldReco -> SetBranchAddress("dca2d",           &rec_dca2d);
1408   m_ntOldReco -> SetBranchAddress("dca2dsigma",      &rec_dca2dsigma);
1409   m_ntOldReco -> SetBranchAddress("dca3dxy",         &rec_dca3dxy);
1410   m_ntOldReco -> SetBranchAddress("dca3dxysigma",    &rec_dca3dxysigma);
1411   m_ntOldReco -> SetBranchAddress("dca3dz",          &rec_dca3dz);
1412   m_ntOldReco -> SetBranchAddress("dca3dzsigma",     &rec_dca3dzsigma);
1413   m_ntOldReco -> SetBranchAddress("pcax",            &rec_pcax);
1414   m_ntOldReco -> SetBranchAddress("pcay",            &rec_pcay);
1415   m_ntOldReco -> SetBranchAddress("pcaz",            &rec_pcaz);
1416   m_ntOldReco -> SetBranchAddress("gtrackID",        &rec_gtrackID);
1417   m_ntOldReco -> SetBranchAddress("gflavor",         &rec_gflavor);
1418   m_ntOldReco -> SetBranchAddress("gnhits",          &rec_gnhits);
1419   m_ntOldReco -> SetBranchAddress("gnmaps",          &rec_gnmaps);
1420   m_ntOldReco -> SetBranchAddress("gnintt",          &rec_gnintt);
1421   m_ntOldReco -> SetBranchAddress("gntpc",           &rec_gntpc);
1422   m_ntOldReco -> SetBranchAddress("gnmms",           &rec_gnmms);
1423   m_ntOldReco -> SetBranchAddress("gnlmaps",         &rec_gnlmaps);
1424   m_ntOldReco -> SetBranchAddress("gnlintt",         &rec_gnlintt);
1425   m_ntOldReco -> SetBranchAddress("gnltpc",          &rec_gnltpc);
1426   m_ntOldReco -> SetBranchAddress("gnlmms",          &rec_gnlmms);
1427   m_ntOldReco -> SetBranchAddress("gpx",             &rec_gpx);
1428   m_ntOldReco -> SetBranchAddress("gpy",             &rec_gpy);
1429   m_ntOldReco -> SetBranchAddress("gpz",             &rec_gpz);
1430   m_ntOldReco -> SetBranchAddress("gpt",             &rec_gpt);
1431   m_ntOldReco -> SetBranchAddress("geta",            &rec_geta);
1432   m_ntOldReco -> SetBranchAddress("gphi",            &rec_gphi);
1433   m_ntOldReco -> SetBranchAddress("gvx",             &rec_gvx);
1434   m_ntOldReco -> SetBranchAddress("gvy",             &rec_gvy);
1435   m_ntOldReco -> SetBranchAddress("gvz",             &rec_gvz);
1436   m_ntOldReco -> SetBranchAddress("gvt",             &rec_gvt);
1437   m_ntOldReco -> SetBranchAddress("gfpx",            &rec_gfpx);
1438   m_ntOldReco -> SetBranchAddress("gfpy",            &rec_gfpy);
1439   m_ntOldReco -> SetBranchAddress("gfpz",            &rec_gfpz);
1440   m_ntOldReco -> SetBranchAddress("gfx",             &rec_gfx);
1441   m_ntOldReco -> SetBranchAddress("gfy",             &rec_gfy);
1442   m_ntOldReco -> SetBranchAddress("gfz",             &rec_gfz);
1443   m_ntOldReco -> SetBranchAddress("gembed",          &rec_gembed);
1444   m_ntOldReco -> SetBranchAddress("gprimary",        &rec_gprimary);
1445   m_ntOldReco -> SetBranchAddress("nfromtruth",      &rec_nfromtruth);
1446   m_ntOldReco -> SetBranchAddress("nwrong",          &rec_nwrong);
1447   m_ntOldReco -> SetBranchAddress("ntrumaps",        &rec_ntrumaps);
1448   m_ntOldReco -> SetBranchAddress("nwrongmaps",      &rec_nwrongmaps);
1449   m_ntOldReco -> SetBranchAddress("ntruintt",        &rec_ntruintt);
1450   m_ntOldReco -> SetBranchAddress("nwrongintt",      &rec_nwrongintt);
1451   m_ntOldReco -> SetBranchAddress("ntrutpc",         &rec_ntrutpc);
1452   m_ntOldReco -> SetBranchAddress("nwrongtpc",       &rec_nwrongtpc);
1453   m_ntOldReco -> SetBranchAddress("ntrumms",         &rec_ntrumms);
1454   m_ntOldReco -> SetBranchAddress("nwrongmms",       &rec_nwrongmms);
1455   m_ntOldReco -> SetBranchAddress("ntrutpc1",        &rec_ntrutpc1);
1456   m_ntOldReco -> SetBranchAddress("nwrongtpc1",      &rec_nwrongtpc1);
1457   m_ntOldReco -> SetBranchAddress("ntrutpc11",       &rec_ntrutpc11);
1458   m_ntOldReco -> SetBranchAddress("nwrongtpc11",     &rec_nwrongtpc11);
1459   m_ntOldReco -> SetBranchAddress("ntrutpc2",        &rec_ntrutpc2);
1460   m_ntOldReco -> SetBranchAddress("nwrongtpc2",      &rec_nwrongtpc2);
1461   m_ntOldReco -> SetBranchAddress("ntrutpc3",        &rec_ntrutpc3);
1462   m_ntOldReco -> SetBranchAddress("nwrongtpc3",      &rec_nwrongtpc3);
1463   m_ntOldReco -> SetBranchAddress("layersfromtruth", &rec_layersfromtruth);
1464   m_ntOldReco -> SetBranchAddress("npedge",          &rec_npedge);
1465   m_ntOldReco -> SetBranchAddress("nredge",          &rec_nredge);
1466   m_ntOldReco -> SetBranchAddress("nbig",            &rec_nbig);
1467   m_ntOldReco -> SetBranchAddress("novlp",           &rec_novlp);
1468   m_ntOldReco -> SetBranchAddress("merr",            &rec_merr);
1469   m_ntOldReco -> SetBranchAddress("msize",           &rec_msize);
1470   m_ntOldReco -> SetBranchAddress("nhittpcall",      &rec_nhittpcall);
1471   m_ntOldReco -> SetBranchAddress("nhittpcin",       &rec_nhittpcin);
1472   m_ntOldReco -> SetBranchAddress("nhittpcmid",      &rec_nhittpcmid);
1473   m_ntOldReco -> SetBranchAddress("nhittpcout",      &rec_nhittpcout);
1474   m_ntOldReco -> SetBranchAddress("nclusall",        &rec_nclusall);
1475   m_ntOldReco -> SetBranchAddress("nclustpc",        &rec_nclustpc);
1476   m_ntOldReco -> SetBranchAddress("nclusintt",       &rec_nclusintt);
1477   m_ntOldReco -> SetBranchAddress("nclusmaps",       &rec_nclusmaps);
1478   m_ntOldReco -> SetBranchAddress("nclusmms",        &rec_nclusmms);
1479   cout << "        Set input branches." << endl;
1480
1481   // grab no. of entries
1482   const int64_t nTrueEntries = m_ntOldTrue -> GetEntries();
1483   const int64_t nRecoEntries = m_ntOldReco -> GetEntries();
1484   cout << "        Beginning truth particle loop: " << nTrueEntries << " to process" << endl;
1485
1486   // loop over truth particles
1487   int64_t nTrueBytes = 0;
1488   for (int64_t iTrueEntry = 0; iTrueEntry < nTrueEntries; iTrueEntry++) {
1489
1490     // grab truth particle entry
1491     const int64_t trueBytes = m_ntOldTrue -> GetEntry(iTrueEntry);
1492     if (trueBytes < 0) {
1493       cerr << "PANIC: issue with entry " << iTrueEntry << "! Aborting loop!\n" << endl;
1494       break;
1495     } else {
1496       nTrueBytes += trueBytes;
1497     }
1498
1499     const int64_t iTrueProg = iTrueEntry + 1;
1500     if (iTrueProg == nTrueEntries) {
1501       cout << "          Processing entry " << iTrueProg << "/" << nTrueEntries << "..." << endl;
1502     } else {
1503       cout << "          Processing entry " << iTrueProg << "/" << nTrueEntries << "...\r" << flush;
1504     }
1505
1506     // skip nan's
1507     if (isnan(tru_trackID)) continue;
1508
1509     // check if near sector
1510     bool isNearSector = false;
1511     if (m_config.doPhiCut) {
1512       isNearSector = IsNearSector(tru_gphi);
1513     }
1514
1515     // apply cuts
1516     const bool isPrimary   = (tru_gprimary == 1);
1517     const bool isInZVtxCut = ((tru_gvz > m_config.zVtxRange.first) && (tru_gvz < m_config.zVtxRange.second));
1518     if (m_config.useOnlyPrimTrks && !isPrimary)   continue;
1519     if (m_config.doZVtxCut       && !isInZVtxCut) continue;
1520     if (m_config.doPhiCut        && isNearSector) continue;
1521
1522     // run calculations
1523     const double tru_gntot = tru_gnintt + tru_gnmaps + tru_gntpc;
1524
1525     // bundle histogram values to fill
1526     STrackMatcherComparatorHistContent content {
1527       tru_gntot,
1528       tru_gnintt,
1529       tru_gnmaps,
1530       tru_gntpc,
1531       1.,
1532       1.,
1533       1.,
1534       1.,
1535       tru_gphi,
1536       tru_geta,
1537       tru_gpt,
1538       1.,
1539       1.,
1540       1.,
1541       1.,
1542       1.,
1543       1.,
1544       1.,
1545       1.
1546     };
1547
1548     // fill truth histograms
1549     FillHistogram1D(content, Type::Truth, m_vecOldHists1D);
1550     FillHistogram2D(content, Type::Truth, Comp::VsTruthPt, tru_gpt,   m_vecOldHists2D);
1551     FillHistogram2D(content, Type::Truth, Comp::VsNumTpc,  tru_gntpc, m_vecOldHists2D);
1552   }  // end truth particle loop
1553
1554   // announce next entry loop
1555   cout << "        Finished truth particle loop.\n"
1556        << "        Beginning reconstructed track loop: " << nRecoEntries << " to process"
1557        << endl;
1558
1559   // loop over reco tracks
1560   int64_t nRecoBytes = 0;
1561   for (int64_t iRecoEntry = 0; iRecoEntry < nRecoEntries; iRecoEntry++) {
1562
1563     // grab reco track entry
1564     const int64_t recoBytes = m_ntOldReco -> GetEntry(iRecoEntry);
1565     if (recoBytes < 0) {
1566       cerr << "PANIC: issue with entry " << iRecoEntry << "! Aborting loop!\n" << endl;
1567       break;
1568     } else {
1569       nRecoBytes += recoBytes;
1570     }
1571
1572     const int64_t iRecoProg = iRecoEntry + 1;
1573     if (iRecoProg == nRecoEntries) {
1574       cout << "          Processing entry " << iRecoProg << "/" << nRecoEntries << "..." << endl;
1575     } else {
1576       cout << "          Processing entry " << iRecoProg << "/" << nRecoEntries << "...\r" << flush;
1577     }
1578
1579     // skip nan's
1580     if (isnan(rec_gpt)) continue;
1581
1582     // check if near sector
1583     bool isNearSector = false;
1584     if (m_config.doPhiCut) {
1585       isNearSector = IsNearSector(rec_phi);
1586     }
1587
1588     // apply cuts
1589     const bool isPrimary   = (rec_gprimary == 1);
1590     const bool isInZVtxCut = ((rec_vz > m_config.zVtxRange.first) && (rec_vz < m_config.zVtxRange.second));
1591     if (m_config.useOnlyPrimTrks && !isPrimary)   continue;
1592     if (m_config.doZVtxCut       && !isInZVtxCut) continue;
1593     if (m_config.doPhiCut        && isNearSector) continue;
1594
1595     // run calculations
1596     const double rec_ntot   = rec_nintt + rec_nmaps + rec_ntpc;
1597     const double rec_gntot  = rec_gnintt + rec_gnmaps + rec_gntpc;
1598     const double rec_rtot   = rec_ntot / rec_gntot;
1599     const double rec_rintt  = rec_nintt / rec_gnintt;
1600     const double rec_rmaps  = rec_nmaps / rec_gnmaps;
1601     const double rec_rtpc   = rec_ntpc / rec_gntpc;
1602     const double rec_ptfrac = rec_pt / rec_gpt;
1603     const double rec_ptper  = rec_deltapt / rec_pt;
1604     const double rec_etaper = rec_deltaeta / rec_eta;
1605     const double rec_phiper = rec_deltaphi / rec_phi;
1606     const double rec_ptres  = abs(rec_pt - rec_gpt) / rec_gpt;
1607     const double rec_etares = abs(rec_eta - rec_geta) / rec_geta;
1608     const double rec_phires = abs(rec_phi - rec_gphi) / rec_gphi;
1609
1610     // bundle histogram values to fill
1611     STrackMatcherComparatorHistContent content {
1612       rec_ntot,
1613       rec_nintt,
1614       rec_nmaps,
1615       rec_ntpc,
1616       rec_rtot,
1617       rec_rintt,
1618       rec_rmaps,
1619       rec_rtpc,
1620       rec_phi,
1621       rec_eta,
1622       rec_pt,
1623       rec_ptfrac,
1624       rec_quality,
1625       rec_ptper,
1626       rec_etaper,
1627       rec_phiper,
1628       rec_ptres,
1629       rec_etares,
1630       rec_phires
1631     };
1632
1633     // fill all matched reco histograms
1634     FillHistogram1D(content, Type::Track, m_vecOldHists1D);
1635     FillHistogram2D(content, Type::Track, Comp::VsTruthPt, rec_gpt,   m_vecOldHists2D);
1636     FillHistogram2D(content, Type::Track, Comp::VsNumTpc,  rec_gntpc, m_vecOldHists2D);
1637
1638     // fill weird and normal matched reco 1D histograms
1639     const bool isNormalTrack = ((rec_ptfrac >= m_config.oddPtFrac.first) && (rec_ptfrac <= m_config.oddPtFrac.second));
1640     if (isNormalTrack) {
1641       FillHistogram1D(content, Type::Normal, m_vecOldHists1D);
1642       FillHistogram2D(content, Type::Normal, Comp::VsTruthPt, rec_gpt,   m_vecOldHists2D);
1643       FillHistogram2D(content, Type::Normal, Comp::VsNumTpc,  rec_gntpc, m_vecOldHists2D);
1644     } else {
1645       FillHistogram1D(content, Type::Weird, m_vecOldHists1D);
1646       FillHistogram2D(content, Type::Weird, Comp::VsTruthPt, rec_gpt,   m_vecOldHists2D);
1647       FillHistogram2D(content, Type::Weird, Comp::VsNumTpc,  rec_gntpc, m_vecOldHists2D);
1648     }
1649   }  // end reco track loop
1650
1651   // announce method end and return
1652   cout << "        Finished reconstructed track loop.\n"
1653        << "      Finished getting histograms from old matcher track tuple."
1654        << endl;
1655   return;
1656
1657 }  // end 'GetOldTupleHists()'
1658
1659
1660
1661 void STrackMatcherComparator::MakeRatiosAndPlots(
1662   const vector<vector<TH1D*>> vecNewHists1D,
1663   const vector<vector<vector<TH2D*>>> vecNewHists2D,
1664   const int iDir,
1665   const string sLabel
1666 ) {
1667
1668   // announce start of routine
1669   const size_t nHistRows1D  = m_vecOldHists1D.size();
1670   const size_t nHistRows2D  = m_vecOldHists2D.size();
1671   const size_t nHistTypes1D = m_vecOldHists1D.front().size();
1672   const size_t nHistTypes2D = m_vecOldHists2D.front().size();
1673   const size_t nHist2D      = m_vecOldHists2D.front().front().size();
1674   cout << "      Making ratios and plots:" << endl;
1675
1676   // create 1d canvas names
1677   vector<vector<string>> vecCanvasNames1D(nHistRows1D);
1678   for (size_t iHistRow1D = 0; iHistRow1D < nHistRows1D; iHistRow1D++) {
1679     vecCanvasNames1D[iHistRow1D].resize(nHistTypes1D);
1680     for (size_t iHistType1D = 0; iHistType1D < nHistTypes1D; iHistType1D++) {
1681       vecCanvasNames1D[iHistRow1D][iHistType1D] = m_hist.vecNameBase[iHistRow1D][iHistType1D];
1682       vecCanvasNames1D[iHistRow1D][iHistType1D].replace(0, 1, "c");
1683     }
1684   }
1685
1686   // create 2d canvas names
1687   vector<vector<vector<string>>> vecCanvasNames2D(nHistRows2D);
1688   for (size_t iHistRow2D = 0; iHistRow2D < nHistRows2D; iHistRow2D++) {
1689     vecCanvasNames2D[iHistRow2D].resize(nHistRows2D);
1690     for (size_t iHistType2D = 0; iHistType2D < nHistTypes2D; iHistType2D++) {
1691       vecCanvasNames2D[iHistRow2D][iHistType2D].resize(nHist2D);
1692       for (size_t iHist2D = 0; iHist2D < nHist2D; iHist2D++) {
1693         vecCanvasNames2D[iHistRow2D][iHistType2D][iHist2D] = m_hist.vecNameBase[iHistRow2D][iHistType2D];
1694         vecCanvasNames2D[iHistRow2D][iHistType2D][iHist2D] += m_hist.vecVsModifiers[iHist2D];
1695         vecCanvasNames2D[iHistRow2D][iHistType2D][iHist2D].replace(0, 1, "c");
1696       }
1697     }
1698   }
1699   cout << "        Set up canvas names." << endl;
1700
1701   // normalize histograms and set axes' ranges as needed ----------------------
1702
1703   // normalize by integral if needed
1704   if (m_config.doIntNorm) {
1705     for (auto oldHistRow1D : m_vecOldHists1D) {
1706       for (auto hOldHist1D : oldHistRow1D) {
1707         const double intOld1D = hOldHist1D -> Integral();
1708         if (intOld1D > 0.) {
1709           hOldHist1D -> Scale(1. / intOld1D);
1710         }
1711       }
1712     }
1713     for (auto newHistRow1D : vecNewHists1D) {
1714       for (auto hNewHist1D : newHistRow1D) {
1715         const double intNew1D = hNewHist1D -> Integral();
1716         if (intNew1D > 0.) {
1717           hNewHist1D -> Scale(1. / intNew1D);
1718         }
1719       }
1720     }
1721     for (auto oldHistRow2D : m_vecOldHists2D) {
1722       for (auto oldHistTypes2D : oldHistRow2D) {
1723         for (auto hOldHist2D : oldHistTypes2D) {
1724           const double intOld2D = hOldHist2D -> Integral();
1725           if (intOld2D > 0.) {
1726             hOldHist2D -> Scale(1. / intOld2D);
1727           }
1728         }
1729       }
1730     }
1731     for (auto newHistRow2D : vecNewHists2D) {
1732       for (auto newHistTypes2D : newHistRow2D) {
1733         for (auto hNewHist2D : newHistTypes2D) {
1734           const double intNew2D = hNewHist2D -> Integral();
1735           if (intNew2D > 0.) {
1736             hNewHist2D -> Scale(1. / intNew2D);
1737           }
1738         }
1739       }
1740     }
1741     cout << "        Normalized histograms." << endl;
1742   }
1743
1744   // set z-axis ranges if needed
1745   if (m_config.matchVertScales) {
1746     for (size_t iHistRow2D = 0; iHistRow2D < nHistRows2D; iHistRow2D++) {
1747       for (size_t iHistType2D = 0; iHistType2D < nHistTypes2D; iHistType2D++) {
1748         for (size_t iHist2D = 0; iHist2D < nHist2D; iHist2D++) {
1749           const float oldMin = m_vecOldHists2D[iHistRow2D][iHistType2D][iHist2D] -> GetMinimum(0.);
1750           const float oldMax = m_vecOldHists2D[iHistRow2D][iHistType2D][iHist2D] -> GetMaximum();
1751           const float newMin = vecNewHists2D[iHistRow2D][iHistType2D][iHist2D]   -> GetMinimum(0.);
1752           const float newMax = vecNewHists2D[iHistRow2D][iHistType2D][iHist2D]   -> GetMaximum();
1753           const float setMin = min(oldMin, newMin);
1754           const float setMax = max(oldMax, newMax);
1755           m_vecOldHists2D[iHistRow2D][iHistType2D][iHist2D] -> GetZaxis() -> SetRangeUser(setMin, setMax);
1756           vecNewHists2D[iHistRow2D][iHistType2D][iHist2D]   -> GetZaxis() -> SetRangeUser(setMin, setMax);
1757         }
1758       }
1759     }
1760     cout << "        Adjusted z-axis scales to match." << endl;
1761   }
1762
1763   // pick relevant count
1764   string countUse("");
1765   if (m_config.doIntNorm) {
1766     countUse = m_config.norm;
1767   } else {
1768     countUse = m_config.count;
1769   }
1770
1771   // calculate ratios ---------------------------------------------------------
1772
1773   // calculate 1d ratios
1774   vector<vector<TH1D*>> vecRatios1D(nHistRows1D);
1775   for (size_t iHistRow1D = 0; iHistRow1D < nHistRows1D; iHistRow1D++) {
1776     vecRatios1D[iHistRow1D].resize(nHistTypes1D);
1777     for (size_t iHist1D = 0; iHist1D < nHistTypes1D; iHist1D++) {
1778
1779       // make histogram name
1780       TString sRatioName = m_vecOldHists1D[iHistRow1D][iHist1D] -> GetName();
1781       sRatioName.Append("_");
1782       sRatioName.Append(sLabel.data());
1783
1784       vecRatios1D[iHistRow1D][iHist1D] = (TH1D*) m_vecOldHists1D[iHistRow1D][iHist1D] -> Clone();
1785       vecRatios1D[iHistRow1D][iHist1D] -> SetName(sRatioName.Data());
1786       vecRatios1D[iHistRow1D][iHist1D] -> Reset("ICES");
1787       vecRatios1D[iHistRow1D][iHist1D] -> Divide(m_vecOldHists1D[iHistRow1D][iHist1D], vecNewHists1D[iHistRow1D][iHist1D], 1., 1.);
1788     }
1789   }
1790
1791   // calculate 2d ratios
1792   vector<vector<vector<TH2D*>>> vecRatios2D(nHistRows2D);
1793   for (size_t iHistRow2D = 0; iHistRow2D < nHistRows2D; iHistRow2D++) {
1794     vecRatios2D[iHistRow2D].resize(nHistTypes2D);
1795     for (size_t iHistType2D = 0; iHistType2D < nHistTypes2D; iHistType2D++) {
1796       vecRatios2D[iHistRow2D][iHistType2D].resize(nHist2D);
1797       for (size_t iHist2D = 0; iHist2D < nHist2D; iHist2D++) {
1798
1799         // make histogram name
1800         TString sRatioName = m_vecOldHists2D[iHistRow2D][iHistType2D][iHist2D] -> GetName();
1801         sRatioName.Append("_");
1802         sRatioName.Append(sLabel.data());
1803
1804         vecRatios2D[iHistRow2D][iHistType2D][iHist2D] = (TH2D*) m_vecOldHists2D[iHistRow2D][iHistType2D][iHist2D] -> Clone();
1805         vecRatios2D[iHistRow2D][iHistType2D][iHist2D] -> SetName(sRatioName.Data());
1806         vecRatios2D[iHistRow2D][iHistType2D][iHist2D] -> Reset("ICES");
1807         vecRatios2D[iHistRow2D][iHistType2D][iHist2D] -> Divide(m_vecOldHists2D[iHistRow2D][iHistType2D][iHist2D], vecNewHists2D[iHistRow2D][iHistType2D][iHist2D], 1., 1.);
1808       }
1809     }
1810   }
1811   cout << "        Calculated ratios." << endl;
1812
1813   // set histogram styles -----------------------------------------------------
1814
1815   // style options
1816   const uint32_t fTxt(42);
1817   const uint32_t fAln(12);
1818   const uint32_t fCnt(1);
1819   const float    fLabH[m_const.nAxes]   = {0.04,  0.04,  0.03};
1820   const float    fLabR1[m_const.nAxes]  = {0.074, 0.074, 0.056};
1821   const float    fLabR2[m_const.nAxes]  = {0.04,  0.04,  0.03};
1822   const float    fTitH[m_const.nAxes]   = {0.04,  0.04,  0.04};
1823   const float    fTitR1[m_const.nAxes]  = {0.074, 0.074, 0.056};
1824   const float    fTitR2[m_const.nAxes]  = {0.04,  0.04,  0.04};
1825   const float    fOffTH[m_const.nAxes]  = {1.0,   1.3,   1.2};
1826   const float    fOffTR1[m_const.nAxes] = {0.8,   0.8,   1.0};
1827   const float    fOffTR2[m_const.nAxes] = {1.0,   1.3,   1.2};
1828   const float    fOffLH[m_const.nAxes]  = {0.005, 0.005, 0.000};
1829   const float    fOffLR1[m_const.nAxes] = {0.005, 0.005, 0.000};
1830   const float    fOffLR2[m_const.nAxes] = {0.005, 0.005, 0.000};
1831
1832   // set old histogram styles
1833   for (auto oldHistRow1D : m_vecOldHists1D) {
1834     for (auto hOldHist1D : oldHistRow1D) {
1835       hOldHist1D -> SetMarkerColor(m_config.fCol.first);
1836       hOldHist1D -> SetMarkerStyle(m_config.fMar.first);
1837       hOldHist1D -> SetFillColor(m_config.fCol.first);
1838       hOldHist1D -> SetFillStyle(m_config.fFil.first);
1839       hOldHist1D -> SetLineColor(m_config.fCol.first);
1840       hOldHist1D -> SetLineStyle(m_config.fLin.first);
1841       hOldHist1D -> SetLineWidth(m_config.fWid.first);
1842       hOldHist1D -> SetTitle("");
1843       hOldHist1D -> SetTitleFont(fTxt);
1844       hOldHist1D -> GetXaxis() -> SetTitleFont(fTxt);
1845       hOldHist1D -> GetXaxis() -> SetTitleSize(fTitH[0]);
1846       hOldHist1D -> GetXaxis() -> SetTitleOffset(fOffTH[0]);
1847       hOldHist1D -> GetXaxis() -> SetLabelFont(fTxt);
1848       hOldHist1D -> GetXaxis() -> SetLabelSize(fLabH[0]);
1849       hOldHist1D -> GetXaxis() -> SetLabelOffset(fOffLH[0]);
1850       hOldHist1D -> GetXaxis() -> CenterTitle(fCnt);
1851       hOldHist1D -> GetYaxis() -> SetTitle(countUse.data());
1852       hOldHist1D -> GetYaxis() -> SetTitleFont(fTxt);
1853       hOldHist1D -> GetYaxis() -> SetTitleSize(fTitH[1]);
1854       hOldHist1D -> GetYaxis() -> SetTitleOffset(fOffTH[1]);
1855       hOldHist1D -> GetYaxis() -> SetLabelFont(fTxt);
1856       hOldHist1D -> GetYaxis() -> SetLabelSize(fLabH[1]);
1857       hOldHist1D -> GetYaxis() -> SetLabelOffset(fOffLH[1]);
1858       hOldHist1D -> GetYaxis() -> CenterTitle(fCnt);
1859     }
1860   }
1861   for (auto oldHistRow2D : m_vecOldHists2D) {
1862     for (auto oldHistTypes2D : oldHistRow2D) {
1863       for (auto hOldHist2D : oldHistTypes2D) {
1864         hOldHist2D -> SetMarkerColor(m_config.fCol.first);
1865         hOldHist2D -> SetMarkerStyle(m_config.fMar.first);
1866         hOldHist2D -> SetFillColor(m_config.fCol.first);
1867         hOldHist2D -> SetFillStyle(m_config.fFil.first);
1868         hOldHist2D -> SetLineColor(m_config.fCol.first);
1869         hOldHist2D -> SetLineStyle(m_config.fLin.first);
1870         hOldHist2D -> SetLineWidth(m_config.fWid.first);
1871         hOldHist2D -> SetTitle(m_config.legOld.data());
1872         hOldHist2D -> SetTitleFont(fTxt);
1873         hOldHist2D -> GetXaxis() -> SetTitleFont(fTxt);
1874         hOldHist2D -> GetXaxis() -> SetTitleSize(fTitH[0]);
1875         hOldHist2D -> GetXaxis() -> SetTitleOffset(fOffTH[0]);
1876         hOldHist2D -> GetXaxis() -> SetLabelFont(fTxt);
1877         hOldHist2D -> GetXaxis() -> SetLabelSize(fLabH[0]);
1878         hOldHist2D -> GetXaxis() -> SetLabelOffset(fOffLH[0]);
1879         hOldHist2D -> GetXaxis() -> CenterTitle(fCnt);
1880         hOldHist2D -> GetYaxis() -> SetTitleFont(fTxt);
1881         hOldHist2D -> GetYaxis() -> SetTitleSize(fTitH[1]);
1882         hOldHist2D -> GetYaxis() -> SetTitleOffset(fOffTH[1]);
1883         hOldHist2D -> GetYaxis() -> SetLabelFont(fTxt);
1884         hOldHist2D -> GetYaxis() -> SetLabelSize(fLabH[1]);
1885         hOldHist2D -> GetYaxis() -> SetLabelOffset(fOffLH[1]);
1886         hOldHist2D -> GetYaxis() -> CenterTitle(fCnt);
1887         hOldHist2D -> GetZaxis() -> SetTitle(countUse.data());
1888         hOldHist2D -> GetZaxis() -> SetTitleFont(fTxt);
1889         hOldHist2D -> GetZaxis() -> SetTitleSize(fTitH[2]);
1890         hOldHist2D -> GetZaxis() -> SetTitleOffset(fOffTH[2]);
1891         hOldHist2D -> GetZaxis() -> SetLabelFont(fTxt);
1892         hOldHist2D -> GetZaxis() -> SetLabelSize(fLabH[2]);
1893         hOldHist2D -> GetZaxis() -> SetLabelOffset(fOffLH[2]);
1894         hOldHist2D -> GetZaxis() -> CenterTitle(fCnt);
1895       }
1896     }
1897   }
1898
1899   // set new histogram styles
1900   for (auto newHistRow1D : vecNewHists1D) {
1901     for (auto hNewHist1D : newHistRow1D) {
1902       hNewHist1D -> SetMarkerColor(m_config.fCol.second);
1903       hNewHist1D -> SetMarkerStyle(m_config.fMar.second);
1904       hNewHist1D -> SetFillColor(m_config.fCol.second);
1905       hNewHist1D -> SetFillStyle(m_config.fFil.second);
1906       hNewHist1D -> SetLineColor(m_config.fCol.second);
1907       hNewHist1D -> SetLineStyle(m_config.fLin.second);
1908       hNewHist1D -> SetLineWidth(m_config.fWid.second);
1909       hNewHist1D -> SetTitle("");
1910       hNewHist1D -> SetTitleFont(fTxt);
1911       hNewHist1D -> GetXaxis() -> SetTitleFont(fTxt);
1912       hNewHist1D -> GetXaxis() -> SetTitleSize(fTitH[0]);
1913       hNewHist1D -> GetXaxis() -> SetTitleOffset(fOffTH[0]);
1914       hNewHist1D -> GetXaxis() -> SetLabelFont(fTxt);
1915       hNewHist1D -> GetXaxis() -> SetLabelSize(fLabH[0]);
1916       hNewHist1D -> GetXaxis() -> SetLabelOffset(fOffLH[0]);
1917       hNewHist1D -> GetXaxis() -> CenterTitle(fCnt);
1918       hNewHist1D -> GetYaxis() -> SetTitle(countUse.data());
1919       hNewHist1D -> GetYaxis() -> SetTitleFont(fTxt);
1920       hNewHist1D -> GetYaxis() -> SetTitleSize(fTitH[1]);
1921       hNewHist1D -> GetYaxis() -> SetTitleOffset(fOffTH[1]);
1922       hNewHist1D -> GetYaxis() -> SetLabelFont(fTxt);
1923       hNewHist1D -> GetYaxis() -> SetLabelSize(fLabH[1]);
1924       hNewHist1D -> GetYaxis() -> SetLabelOffset(fOffLH[1]);
1925       hNewHist1D -> GetYaxis() -> CenterTitle(fCnt);
1926     }
1927   }
1928   for (auto newHistRow2D : vecNewHists2D) {
1929     for (auto newHistTypes2D : newHistRow2D) {
1930       for (auto hNewHist2D : newHistTypes2D) {
1931         hNewHist2D -> SetMarkerColor(m_config.fCol.first);
1932         hNewHist2D -> SetMarkerStyle(m_config.fMar.first);
1933         hNewHist2D -> SetFillColor(m_config.fCol.first);
1934         hNewHist2D -> SetFillStyle(m_config.fFil.first);
1935         hNewHist2D -> SetLineColor(m_config.fCol.first);
1936         hNewHist2D -> SetLineStyle(m_config.fLin.first);
1937         hNewHist2D -> SetLineWidth(m_config.fWid.first);
1938         hNewHist2D -> SetTitle(m_config.legNew.data());
1939         hNewHist2D -> SetTitleFont(fTxt);
1940         hNewHist2D -> GetXaxis() -> SetTitleFont(fTxt);
1941         hNewHist2D -> GetXaxis() -> SetTitleSize(fTitH[0]);
1942         hNewHist2D -> GetXaxis() -> SetTitleOffset(fOffTH[0]);
1943         hNewHist2D -> GetXaxis() -> SetLabelFont(fTxt);
1944         hNewHist2D -> GetXaxis() -> SetLabelSize(fLabH[0]);
1945         hNewHist2D -> GetXaxis() -> SetLabelOffset(fOffLH[0]);
1946         hNewHist2D -> GetXaxis() -> CenterTitle(fCnt);
1947         hNewHist2D -> GetYaxis() -> SetTitleFont(fTxt);
1948         hNewHist2D -> GetYaxis() -> SetTitleSize(fTitH[1]);
1949         hNewHist2D -> GetYaxis() -> SetTitleOffset(fOffTH[1]);
1950         hNewHist2D -> GetYaxis() -> SetLabelFont(fTxt);
1951         hNewHist2D -> GetYaxis() -> SetLabelSize(fLabH[1]);
1952         hNewHist2D -> GetYaxis() -> SetLabelOffset(fOffLH[1]);
1953         hNewHist2D -> GetYaxis() -> CenterTitle(fCnt);
1954         hNewHist2D -> GetZaxis() -> SetTitle(countUse.data());
1955         hNewHist2D -> GetZaxis() -> SetTitleFont(fTxt);
1956         hNewHist2D -> GetZaxis() -> SetTitleSize(fTitH[2]);
1957         hNewHist2D -> GetZaxis() -> SetTitleOffset(fOffTH[2]);
1958         hNewHist2D -> GetZaxis() -> SetLabelFont(fTxt);
1959         hNewHist2D -> GetZaxis() -> SetLabelSize(fLabH[2]);
1960         hNewHist2D -> GetZaxis() -> SetLabelOffset(fOffLH[2]);
1961         hNewHist2D -> GetZaxis() -> CenterTitle(fCnt);
1962       }
1963     }
1964   }
1965
1966   // set ratio styles
1967   for (auto ratioRow1D : vecRatios1D) {
1968     for (auto hRatio1D : ratioRow1D) {
1969       hRatio1D -> SetMarkerColor(m_config.fCol.first);
1970       hRatio1D -> SetMarkerStyle(m_config.fMar.first);
1971       hRatio1D -> SetFillColor(m_config.fCol.first);
1972       hRatio1D -> SetFillStyle(m_config.fFil.first);
1973       hRatio1D -> SetLineColor(m_config.fCol.first);
1974       hRatio1D -> SetLineStyle(m_config.fLin.first);
1975       hRatio1D -> SetLineWidth(m_config.fWid.first);
1976       hRatio1D -> SetTitle("");
1977       hRatio1D -> SetTitleFont(fTxt);
1978       hRatio1D -> GetXaxis() -> SetTitleFont(fTxt);
1979       hRatio1D -> GetXaxis() -> SetTitleSize(fTitR1[0]);
1980       hRatio1D -> GetXaxis() -> SetTitleOffset(fOffTR1[0]);
1981       hRatio1D -> GetXaxis() -> SetLabelFont(fTxt);
1982       hRatio1D -> GetXaxis() -> SetLabelSize(fLabR1[0]);
1983       hRatio1D -> GetXaxis() -> SetLabelOffset(fOffLR1[0]);
1984       hRatio1D -> GetXaxis() -> CenterTitle(fCnt);
1985       hRatio1D -> GetYaxis() -> SetTitle(m_config.ratio.data());
1986       hRatio1D -> GetYaxis() -> SetTitleFont(fTxt);
1987       hRatio1D -> GetYaxis() -> SetTitleSize(fTitR1[1]);
1988       hRatio1D -> GetYaxis() -> SetTitleOffset(fOffTR1[1]);
1989       hRatio1D -> GetYaxis() -> SetLabelFont(fTxt);
1990       hRatio1D -> GetYaxis() -> SetLabelSize(fLabR1[1]);
1991       hRatio1D -> GetYaxis() -> SetLabelOffset(fOffLR1[1]);
1992       hRatio1D -> GetYaxis() -> CenterTitle(fCnt);
1993     }
1994   }
1995   for (auto ratioRow2D : vecRatios2D) {
1996     for (auto ratioTypes2D : ratioRow2D) {
1997       for (auto hRatio2D : ratioTypes2D) {
1998         hRatio2D -> SetMarkerColor(m_config.fCol.first);
1999         hRatio2D -> SetMarkerStyle(m_config.fMar.first);
2000         hRatio2D -> SetFillColor(m_config.fCol.first);
2001         hRatio2D -> SetFillStyle(m_config.fFil.first);
2002         hRatio2D -> SetLineColor(m_config.fCol.first);
2003         hRatio2D -> SetLineStyle(m_config.fLin.first);
2004         hRatio2D -> SetLineWidth(m_config.fWid.first);
2005         hRatio2D -> SetTitle("");
2006         hRatio2D -> SetTitleFont(fTxt);
2007         hRatio2D -> GetXaxis() -> SetTitleFont(fTxt);
2008         hRatio2D -> GetXaxis() -> SetTitleSize(fTitR2[0]);
2009         hRatio2D -> GetXaxis() -> SetTitleOffset(fOffTR2[0]);
2010         hRatio2D -> GetXaxis() -> SetLabelFont(fTxt);
2011         hRatio2D -> GetXaxis() -> SetLabelSize(fLabR2[0]);
2012         hRatio2D -> GetXaxis() -> SetLabelOffset(fOffLR2[0]);
2013         hRatio2D -> GetXaxis() -> CenterTitle(fCnt);
2014         hRatio2D -> GetYaxis() -> SetTitleFont(fTxt);
2015         hRatio2D -> GetYaxis() -> SetTitleSize(fTitR2[1]);
2016         hRatio2D -> GetYaxis() -> SetTitleOffset(fOffTR2[1]);
2017         hRatio2D -> GetYaxis() -> SetLabelFont(fTxt);
2018         hRatio2D -> GetYaxis() -> SetLabelSize(fLabR2[1]);
2019         hRatio2D -> GetYaxis() -> SetLabelOffset(fOffLR2[1]);
2020         hRatio2D -> GetYaxis() -> CenterTitle(fCnt);
2021         hRatio2D -> GetZaxis() -> SetTitle(m_config.ratio.data());
2022         hRatio2D -> GetZaxis() -> SetTitleFont(fTxt);
2023         hRatio2D -> GetZaxis() -> SetTitleSize(fTitR2[2]);
2024         hRatio2D -> GetZaxis() -> SetTitleOffset(fOffTR2[2]);
2025         hRatio2D -> GetZaxis() -> SetLabelFont(fTxt);
2026         hRatio2D -> GetZaxis() -> SetLabelSize(fLabR2[2]);
2027         hRatio2D -> GetZaxis() -> SetLabelOffset(fOffLR2[2]);
2028         hRatio2D -> GetZaxis() -> CenterTitle(fCnt);
2029       }
2030     }
2031   }
2032   cout << "        Set styles." << endl;
2033
2034   // make legends and text boxes ----------------------------------------------
2035
2036   // make legend
2037   const uint32_t fColLe(0);
2038   const uint32_t fFilLe(0);
2039   const uint32_t fLinLe(0);
2040   const float    fLegXY[m_const.nVtx] = {0.1, 0.1, 0.3, 0.2};
2041
2042   TLegend *leg = new TLegend(fLegXY[0], fLegXY[1], fLegXY[2], fLegXY[3], m_config.header.data());
2043   leg -> SetFillColor(fColLe);
2044   leg -> SetFillStyle(fFilLe);
2045   leg -> SetLineColor(fColLe);
2046   leg -> SetLineStyle(fLinLe);
2047   leg -> SetTextFont(fTxt);
2048   leg -> SetTextAlign(fAln);
2049   leg -> AddEntry(m_vecOldHists1D.front().front(), m_config.legOld.data(), "pf");
2050   leg -> AddEntry(vecNewHists1D.front().front(),   m_config.legNew.data(), "pf");
2051   cout << "        Made legend." << endl;
2052
2053   // make text
2054   const uint32_t fColTx(0);
2055   const uint32_t fFilTx(0);
2056   const uint32_t fLinTx(0);
2057   const float    fTxtXY[m_const.nVtx] = {0.3, 0.1, 0.5, 0.25};
2058
2059   TPaveText *txt = new TPaveText(fTxtXY[0], fTxtXY[1], fTxtXY[2], fTxtXY[3], "NDC NB");
2060   txt -> SetFillColor(fColTx);
2061   txt -> SetFillStyle(fFilTx);
2062   txt -> SetLineColor(fColTx);
2063   txt -> SetLineStyle(fLinTx);
2064   txt -> SetTextFont(fTxt);
2065   txt -> SetTextAlign(fAln);
2066   for (string txtLine : m_config.info) {
2067     txt -> AddText(txtLine.data());
2068   }
2069   cout << "        Made text." << endl;
2070
2071   // make plots ---------------------------------------------------------------
2072
2073   // canvas parameters
2074   const uint32_t width1D(750);
2075   const uint32_t width2D(1500);
2076   const uint32_t width2DR(2250);
2077   const uint32_t height(750);
2078   const uint32_t heightR1(900);
2079   const uint32_t heightR2(500);
2080   const uint32_t fMode(0);
2081   const uint32_t fBord(2);
2082   const uint32_t fGrid(0);
2083   const uint32_t fTick(1);
2084   const uint32_t fLogX(0);
2085   const uint32_t fLogY(1);
2086   const uint32_t fLogZ(1);
2087   const uint32_t fFrame(0);
2088   const string   sOldPadName("pOld");
2089   const string   sNewPadName("pNew");
2090   const string   sHistPadName("pHist");
2091   const string   sRatPadName("pRatio");
2092   const float    fMargin1D[m_const.nSide]  = {0.02,  0.02, 0.15,  0.15};
2093   const float    fMargin1DH[m_const.nSide] = {0.02,  0.02, 0.005, 0.15};
2094   const float    fMargin1DR[m_const.nSide] = {0.005, 0.02, 0.2,   0.15};
2095   const float    fMargin2D[m_const.nSide]  = {0.10, 0.15, 0.15, 0.15};
2096   const float    xyOldPad[m_const.nVtx]    = {0.0,  0.0,  0.5,  1.0};
2097   const float    xyOldPadR[m_const.nVtx]   = {0.0,  0.0,  0.33, 1.0};
2098   const float    xyNewPad[m_const.nVtx]    = {0.5,  0.0,  1.0,  1.0};
2099   const float    xyNewPadR[m_const.nVtx]   = {0.33, 0.0,  0.66, 1.0};
2100   const float    xyHistPadR[m_const.nVtx]  = {0.0,  0.33, 1.0,  1.0};
2101   const float    xyRatPadR1[m_const.nVtx]  = {0.0,  0.0,  1.0,  0.33};
2102   const float    xyRatPadR2[m_const.nVtx]  = {0.66, 0.0,  1.0,  1.0};
2103
2104   // make 1d plot without ratio
2105   for (size_t iHistRow1D = 0; iHistRow1D < nHistRows1D; iHistRow1D++) {
2106     for (size_t iHist1D = 0; iHist1D < nHistTypes1D; iHist1D++) {
2107
2108       // make new name
2109       const string sName = vecCanvasNames1D[iHistRow1D][iHist1D] + "_" + sLabel;
2110
2111
2112       // construct canvas
2113       TCanvas* cPlot1D = new TCanvas(sName.data(), "", width1D, height);
2114       cPlot1D -> SetGrid(fGrid, fGrid);
2115       cPlot1D -> SetTicks(fTick, fTick);
2116       cPlot1D -> SetBorderMode(fMode);
2117       cPlot1D -> SetBorderSize(fBord);
2118       cPlot1D -> SetFrameBorderMode(fFrame);
2119       cPlot1D -> SetTopMargin(fMargin1D[0]);
2120       cPlot1D -> SetRightMargin(fMargin1D[1]);
2121       cPlot1D -> SetBottomMargin(fMargin1D[2]);
2122       cPlot1D -> SetLeftMargin(fMargin1D[3]);
2123       cPlot1D -> SetLogx(fLogX);
2124       cPlot1D -> SetLogy(fLogY);
2125       cPlot1D -> cd();
2126
2127       // draw old vs. new histograms
2128       m_vecOldHists1D[iHistRow1D][iHist1D] -> Draw();
2129       vecNewHists1D[iHistRow1D][iHist1D]   -> Draw("same");
2130
2131       // draw text and save
2132       leg                    -> Draw();
2133       txt                    -> Draw();
2134       m_vecPlotDirs.at(iDir) -> cd();
2135       cPlot1D                -> Write();
2136       cPlot1D                -> Close();
2137     }
2138   }
2139
2140   // make 1d plot with ratio
2141   for (size_t iHistRow1D = 0; iHistRow1D < nHistRows1D; iHistRow1D++) {
2142     for (size_t iHist1D = 0; iHist1D < nHistTypes1D; iHist1D++) {
2143
2144       // make new name
2145       const string sNameWithRatio = vecCanvasNames1D[iHistRow1D][iHist1D] + "WithRatio_" + sLabel;
2146
2147       // construct canvas
2148       TCanvas* cPlot1D = new TCanvas(sNameWithRatio.data(), "", width1D, heightR1);
2149       TPad*    pPadH1D = new TPad(sHistPadName.data(), "", xyHistPadR[0], xyHistPadR[1], xyHistPadR[2], xyHistPadR[3]);
2150       TPad*    pPadR1D = new TPad(sRatPadName.data(),  "", xyRatPadR1[0], xyRatPadR1[1], xyRatPadR1[2], xyRatPadR1[3]);
2151       cPlot1D -> SetGrid(fGrid, fGrid);
2152       cPlot1D -> SetTicks(fTick, fTick);
2153       cPlot1D -> SetBorderMode(fMode);
2154       cPlot1D -> SetBorderSize(fBord);
2155       pPadH1D -> SetGrid(fGrid, fGrid);
2156       pPadH1D -> SetTicks(fTick, fTick);
2157       pPadH1D -> SetBorderMode(fMode);
2158       pPadH1D -> SetBorderSize(fBord);
2159       pPadH1D -> SetFrameBorderMode(fFrame);
2160       pPadH1D -> SetTopMargin(fMargin1DH[0]);
2161       pPadH1D -> SetRightMargin(fMargin1DH[1]);
2162       pPadH1D -> SetBottomMargin(fMargin1DH[2]);
2163       pPadH1D -> SetLeftMargin(fMargin1DH[3]);
2164       pPadH1D -> SetLogx(fLogX);
2165       pPadH1D -> SetLogy(fLogY);
2166       pPadR1D -> SetGrid(fGrid, fGrid);
2167       pPadR1D -> SetTicks(fTick, fTick);
2168       pPadR1D -> SetBorderMode(fMode);
2169       pPadR1D -> SetBorderSize(fBord);
2170       pPadR1D -> SetFrameBorderMode(fFrame);
2171       pPadR1D -> SetTopMargin(fMargin1DR[0]);
2172       pPadR1D -> SetRightMargin(fMargin1DR[1]);
2173       pPadR1D -> SetBottomMargin(fMargin1DR[2]);
2174       pPadR1D -> SetLeftMargin(fMargin1DR[3]);
2175       pPadR1D -> SetLogx(fLogX);
2176       pPadR1D -> SetLogy(fLogY);
2177       cPlot1D -> cd();
2178       pPadH1D -> Draw();
2179       pPadR1D -> Draw();
2180
2181       // draw old vs. new histograms
2182       pPadH1D                              -> cd();
2183       m_vecOldHists1D[iHistRow1D][iHist1D] -> Draw();
2184       vecNewHists1D[iHistRow1D][iHist1D]   -> Draw("same");
2185       leg                                  -> Draw();
2186       txt                                  -> Draw();
2187
2188       // draw ratio
2189       pPadR1D                          -> cd();
2190       vecRatios1D[iHistRow1D][iHist1D] -> Draw();
2191
2192       // draw text and save
2193       m_vecPlotDirs.at(iDir) -> cd();
2194       cPlot1D                -> Write();
2195       cPlot1D                -> Close();
2196     }
2197   }
2198
2199   // make 2d plot without ratio
2200   for (size_t iHistRow2D = 0; iHistRow2D < nHistRows2D; iHistRow2D++) {
2201     for (size_t iHistType2D = 0; iHistType2D < nHistTypes2D; iHistType2D++) {
2202       for (size_t iHist2D = 0; iHist2D < nHist2D; iHist2D++) {
2203
2204         // construct canvas
2205         TCanvas* cPlot2D = new TCanvas(vecCanvasNames2D[iHistRow2D][iHistType2D][iHist2D].data(), "", width2D, height);
2206         TPad*    pOld    = new TPad(sOldPadName.data(), "", xyOldPad[0], xyOldPad[1], xyOldPad[2], xyOldPad[3]);
2207         TPad*    pNew    = new TPad(sNewPadName.data(), "", xyNewPad[0], xyNewPad[1], xyNewPad[2], xyNewPad[3]);
2208         cPlot2D -> SetGrid(fGrid, fGrid);
2209         cPlot2D -> SetTicks(fTick, fTick);
2210         cPlot2D -> SetBorderMode(fMode);
2211         cPlot2D -> SetBorderSize(fBord);
2212         pOld    -> SetGrid(fGrid, fGrid);
2213         pOld    -> SetTicks(fTick, fTick);
2214         pOld    -> SetBorderMode(fMode);
2215         pOld    -> SetBorderSize(fBord);
2216         pOld    -> SetFrameBorderMode(fFrame);
2217         pOld    -> SetTopMargin(fMargin2D[0]);
2218         pOld    -> SetRightMargin(fMargin2D[1]);
2219         pOld    -> SetBottomMargin(fMargin2D[2]);
2220         pOld    -> SetLeftMargin(fMargin2D[3]);
2221         pOld    -> SetLogx(fLogX);
2222         pOld    -> SetLogy(fLogY);
2223         pOld    -> SetLogz(fLogZ);
2224         pNew    -> SetGrid(fGrid, fGrid);
2225         pNew    -> SetTicks(fTick, fTick);
2226         pNew    -> SetBorderMode(fMode);
2227         pNew    -> SetBorderSize(fBord);
2228         pNew    -> SetFrameBorderMode(fFrame);
2229         pNew    -> SetTopMargin(fMargin2D[0]);
2230         pNew    -> SetRightMargin(fMargin2D[1]);
2231         pNew    -> SetBottomMargin(fMargin2D[2]);
2232         pNew    -> SetLeftMargin(fMargin2D[3]);
2233         pNew    -> SetLogx(fLogX);
2234         pNew    -> SetLogy(fLogY);
2235         pNew    -> SetLogz(fLogZ);
2236         cPlot2D -> cd();
2237         pOld    -> Draw();
2238         pNew    -> Draw();
2239
2240         // draw old vs. new histograms
2241         pOld                                              -> cd();
2242         m_vecOldHists2D[iHistRow2D][iHistType2D][iHist2D] -> Draw("colz");
2243         pNew                                              -> cd();
2244         vecNewHists2D[iHistRow2D][iHistType2D][iHist2D]   -> Draw("colz");
2245
2246         // draw text and save
2247         pNew                   -> cd();
2248         txt                    -> Draw();
2249         m_vecPlotDirs.at(iDir) -> cd();
2250         cPlot2D                -> Write();
2251         cPlot2D                -> Close();
2252       }
2253     }
2254   }
2255
2256   // make 2d plot with ratio
2257   for (size_t iHistRow2D = 0; iHistRow2D < nHistRows2D; iHistRow2D++) {
2258     for (size_t iHistType2D = 0; iHistType2D < nHistTypes2D; iHistType2D++) {
2259       for (size_t iHist2D = 0; iHist2D < nHist2D; iHist2D++) {
2260
2261         // make new name
2262         const string sNameWithRatio = vecCanvasNames2D[iHistRow2D][iHistType2D][iHist2D] + "_" + sLabel;
2263
2264         // construct canvas
2265         TCanvas* cPlot2D = new TCanvas(sNameWithRatio.data(), "", width2DR, heightR2);
2266         TPad*    pOld    = new TPad(sOldPadName.data(), "", xyOldPadR[0],  xyOldPadR[1],  xyOldPadR[2],  xyOldPadR[3]);
2267         TPad*    pNew    = new TPad(sNewPadName.data(), "", xyNewPadR[0],  xyNewPadR[1],  xyNewPadR[2],  xyNewPadR[3]);
2268         TPad*    pRat    = new TPad(sRatPadName.data(), "", xyRatPadR2[0], xyRatPadR2[1], xyRatPadR2[2], xyRatPadR2[3]);
2269         cPlot2D -> SetGrid(fGrid, fGrid);
2270         cPlot2D -> SetTicks(fTick, fTick);
2271         cPlot2D -> SetBorderMode(fMode);
2272         cPlot2D -> SetBorderSize(fBord);
2273         pOld    -> SetGrid(fGrid, fGrid);
2274         pOld    -> SetTicks(fTick, fTick);
2275         pOld    -> SetBorderMode(fMode);
2276         pOld    -> SetBorderSize(fBord);
2277         pOld    -> SetFrameBorderMode(fFrame);
2278         pOld    -> SetTopMargin(fMargin2D[0]);
2279         pOld    -> SetRightMargin(fMargin2D[1]);
2280         pOld    -> SetBottomMargin(fMargin2D[2]);
2281         pOld    -> SetLeftMargin(fMargin2D[3]);
2282         pOld    -> SetLogx(fLogX);
2283         pOld    -> SetLogy(fLogY);
2284         pOld    -> SetLogz(fLogZ);
2285         pNew    -> SetGrid(fGrid, fGrid);
2286         pNew    -> SetTicks(fTick, fTick);
2287         pNew    -> SetBorderMode(fMode);
2288         pNew    -> SetBorderSize(fBord);
2289         pNew    -> SetFrameBorderMode(fFrame);
2290         pNew    -> SetTopMargin(fMargin2D[0]);
2291         pNew    -> SetRightMargin(fMargin2D[1]);
2292         pNew    -> SetBottomMargin(fMargin2D[2]);
2293         pNew    -> SetLeftMargin(fMargin2D[3]);
2294         pNew    -> SetLogx(fLogX);
2295         pNew    -> SetLogy(fLogY);
2296         pNew    -> SetLogz(fLogZ);
2297         pRat    -> SetGrid(fGrid, fGrid);
2298         pRat    -> SetTicks(fTick, fTick);
2299         pRat    -> SetBorderMode(fMode);
2300         pRat    -> SetBorderSize(fBord);
2301         pRat    -> SetFrameBorderMode(fFrame);
2302         pRat    -> SetTopMargin(fMargin2D[0]);
2303         pRat    -> SetRightMargin(fMargin2D[1]);
2304         pRat    -> SetBottomMargin(fMargin2D[2]);
2305         pRat    -> SetLeftMargin(fMargin2D[3]);
2306         pRat    -> SetLogx(fLogX);
2307         pRat    -> SetLogy(fLogY);
2308         pRat    -> SetLogz(fLogZ);
2309         cPlot2D -> cd();
2310         pOld    -> Draw();
2311         pNew    -> Draw();
2312         pRat    -> Draw();
2313
2314         // draw old vs. new vs. ratio histograms
2315         pOld                                              -> cd();
2316         m_vecOldHists2D[iHistRow2D][iHistType2D][iHist2D] -> Draw("colz");
2317         txt                                               -> Draw();
2318         pNew                                              -> cd();
2319         vecNewHists2D[iHistRow2D][iHistType2D][iHist2D]   -> Draw("colz");
2320         pRat                                              -> cd();
2321         vecRatios2D[iHistRow2D][iHistType2D][iHist2D]     -> Draw("colz");
2322
2323         // draw text and save
2324         m_vecPlotDirs.at(iDir) -> cd();
2325         cPlot2D                -> Write();
2326         cPlot2D                -> Close();
2327       }
2328     }
2329   }
2330   cout << "        Made 2D plot with ratio." << endl;
2331
2332   // save ratio histograms
2333   m_vecRatioDirs.at(iDir) -> cd();
2334   for (auto histRow1D : vecRatios1D) {
2335     for (auto hRatio1D : histRow1D) {
2336       hRatio1D -> Write();
2337     }
2338   }
2339   for (auto histRow2D : vecRatios2D) {
2340     for (auto histTypes2D : histRow2D) {
2341       for (auto hRatio2D : histTypes2D) {
2342         hRatio2D -> Write();
2343       }
2344     }
2345   }
2346
2347   // announce end and return
2348   cout << "       Saved ratio histograms.\n"
2349        << "      Finished making ratios and plots."
2350        << endl;
2351   return;
2352
2353 }  // end 'MakeRatiosAndPlots(vector<TH1D*>, vector<TH2D*>, int, string)'
2354
2355
2356
2357 void STrackMatcherComparator::SaveHistograms() {
2358
2359   // save 1d histograms
2360   for (size_t iHistRow = 0; iHistRow < m_vecTreeHists1D.size(); iHistRow++) {
2361     for (size_t iHistType = 0; iHistType < m_vecTreeHists1D[iHistRow].size(); iHistType++) {
2362       m_vecHistDirs[Src::NewTree]            -> cd();
2363       m_vecTreeHists1D[iHistRow][iHistType]  -> Write();
2364       m_vecHistDirs[Src::NewTuple]           -> cd();
2365       m_vecTupleHists1D[iHistRow][iHistType] -> Write();
2366       m_vecHistDirs[Src::OldTuple]           -> cd();
2367       m_vecOldHists1D[iHistRow][iHistType]   -> Write();
2368     }  // end type loop
2369   }  // end row loop
2370   cout << "      Saved 1d histograms." << endl;
2371
2372   // save 2d histograms
2373   for (size_t iHistRow = 0; iHistRow < m_vecTreeHists2D.size(); iHistRow++) {
2374     for (size_t iHistType = 0; iHistType < m_vecTreeHists2D[iHistRow].size(); iHistType++) {
2375       for (size_t iVsHistMod = 0; iVsHistMod < m_vecTreeHists2D[iHistRow][iHistType].size(); iVsHistMod++) {
2376         m_vecHistDirs[Src::NewTree]                        -> cd();
2377         m_vecTreeHists2D[iHistRow][iHistType][iVsHistMod]  -> Write();
2378         m_vecHistDirs[Src::NewTuple]                       -> cd();
2379         m_vecTupleHists2D[iHistRow][iHistType][iVsHistMod] -> Write();
2380         m_vecHistDirs[Src::OldTuple]                       -> cd();
2381         m_vecOldHists2D[iHistRow][iHistType][iVsHistMod]   -> Write();
2382       }  // end modifier loop
2383     }  // end type loop
2384   }  // end row loop
2385
2386   // announce end and return
2387   cout << "      Saved 2d histograms." << endl;
2388   return;
2389
2390 }  // end 'SaveHistograms()'
2391
2392
2393
2394 void STrackMatcherComparator::CloseInput() {
2395
2396   m_treeInFileTrue  -> cd();
2397   m_treeInFileTrue  -> Close();
2398   m_treeInFileReco  -> cd();
2399   m_treeInFileReco  -> Close();
2400   m_tupleInFileTrue -> cd();
2401   m_tupleInFileTrue -> Close();
2402   m_tupleInFileReco -> cd();
2403   m_tupleInFileReco -> Close();
2404   m_oldInFileTrue   -> cd();
2405   m_oldInFileTrue   -> Close();
2406   m_oldInFileReco   -> cd();
2407   m_oldInFileReco   -> Close();
2408
2409   // announce end and return
2410   cout << "      Closed input files." << endl;
2411   return;
2412
2413 }  // end 'CloseInput()'
2414
2415
2416
2417 void STrackMatcherComparator::CloseOutput() {
2418
2419   m_outFile -> cd();
2420   m_outFile -> Close();
2421
2422   // announce end and return
2423   cout << "      Closed output file." << endl;
2424   return;
2425
2426 }  // end 'CloseOutput()'
2427
2428
2429
2430 // helper functions -----------------------------------------------------------
2431
2432 void STrackMatcherComparator::FillHistogram1D(
2433   const STrackMatcherComparatorHistContent& content,
2434   const Type type,
2435   const vector<vector<TH1D*>> vecHist1D
2436 ) {
2437
2438   vecHist1D.at(Var::NTot).at(type)   -> Fill(content.nTot);
2439   vecHist1D.at(Var::NIntt).at(type)  -> Fill(content.nIntt);
2440   vecHist1D.at(Var::NMvtx).at(type)  -> Fill(content.nMvtx);
2441   vecHist1D.at(Var::NTpc).at(type)   -> Fill(content.nTpc);
2442   vecHist1D.at(Var::RTot).at(type)   -> Fill(content.rTot);
2443   vecHist1D.at(Var::RIntt).at(type)  -> Fill(content.rIntt);
2444   vecHist1D.at(Var::RMvtx).at(type)  -> Fill(content.rMvtx);
2445   vecHist1D.at(Var::RTpc).at(type)   -> Fill(content.rTpc);
2446   vecHist1D.at(Var::Phi).at(type)    -> Fill(content.phi);
2447   vecHist1D.at(Var::Eta).at(type)    -> Fill(content.eta);
2448   vecHist1D.at(Var::Pt).at(type)     -> Fill(content.pt);
2449   vecHist1D.at(Var::Frac).at(type)   -> Fill(content.ptFrac);
2450   vecHist1D.at(Var::Qual).at(type)   -> Fill(content.quality);
2451   vecHist1D.at(Var::PtErr).at(type)  -> Fill(content.ptErr);
2452   vecHist1D.at(Var::EtaErr).at(type) -> Fill(content.etaErr);
2453   vecHist1D.at(Var::PhiErr).at(type) -> Fill(content.phiErr);
2454   vecHist1D.at(Var::PtRes).at(type)  -> Fill(content.ptRes);
2455   vecHist1D.at(Var::EtaRes).at(type) -> Fill(content.etaRes);
2456   vecHist1D.at(Var::PhiRes).at(type) -> Fill(content.phiRes);
2457   return;
2458
2459 }  // end 'FillHistogram1D(STrackMatcherComparatorHistContent&, int, vector<vector<TH1D*>>)'
2460
2461
2462
2463
2464 void STrackMatcherComparator::FillHistogram2D(
2465   const STrackMatcherComparatorHistContent& content,
2466   const Type type,
2467   const Comp comparison,
2468   const double value,
2469   const vector<vector<vector<TH2D*>>> vecHist2D
2470
2471 ) {
2472
2473   vecHist2D.at(Var::NTot).at(type).at(comparison)   -> Fill(value, content.nTot);
2474   vecHist2D.at(Var::NIntt).at(type).at(comparison)  -> Fill(value, content.nIntt);
2475   vecHist2D.at(Var::NMvtx).at(type).at(comparison)  -> Fill(value, content.nMvtx);
2476   vecHist2D.at(Var::NTpc).at(type).at(comparison)   -> Fill(value, content.nTpc);
2477   vecHist2D.at(Var::RTot).at(type).at(comparison)   -> Fill(value, content.rTot);
2478   vecHist2D.at(Var::RIntt).at(type).at(comparison)  -> Fill(value, content.rIntt);
2479   vecHist2D.at(Var::RMvtx).at(type).at(comparison)  -> Fill(value, content.rMvtx);
2480   vecHist2D.at(Var::RTpc).at(type).at(comparison)   -> Fill(value, content.rTpc);
2481   vecHist2D.at(Var::Phi).at(type).at(comparison)    -> Fill(value, content.phi);
2482   vecHist2D.at(Var::Eta).at(type).at(comparison)    -> Fill(value, content.eta);
2483   vecHist2D.at(Var::Pt).at(type).at(comparison)     -> Fill(value, content.pt);
2484   vecHist2D.at(Var::Frac).at(type).at(comparison)   -> Fill(value, content.ptFrac);
2485   vecHist2D.at(Var::Qual).at(type).at(comparison)   -> Fill(value, content.quality);
2486   vecHist2D.at(Var::PtErr).at(type).at(comparison)  -> Fill(value, content.ptErr);
2487   vecHist2D.at(Var::EtaErr).at(type).at(comparison) -> Fill(value, content.etaErr);
2488   vecHist2D.at(Var::PhiErr).at(type).at(comparison) -> Fill(value, content.phiErr);
2489   vecHist2D.at(Var::PtRes).at(type).at(comparison)  -> Fill(value, content.ptRes);
2490   vecHist2D.at(Var::EtaRes).at(type).at(comparison) -> Fill(value, content.etaRes);
2491   vecHist2D.at(Var::PhiRes).at(type).at(comparison) -> Fill(value, content.phiRes);
2492   return;
2493
2494 }  // end 'FillVsHistogram2D(STrackMatcherComparator&, int, int, double, vector<vector<vector<TH2D*>>>)'
2495
2496
2497
2498 bool STrackMatcherComparator::IsNearSector(const float phi) {
2499
2500   bool isNearSector = false;
2501   for (size_t iSector = 0; iSector < m_const.nSectors; iSector++) {
2502     const float cutVal = m_config.sigCutVal * m_config.phiSectors[iSector].second;
2503     const float minPhi = m_config.phiSectors[iSector].first - cutVal;
2504     const float maxPhi = m_config.phiSectors[iSector].first + cutVal;
2505     const bool  isNear = ((phi >= minPhi) && (phi <= maxPhi));
2506     if (isNear) {
2507       isNearSector = true;
2508       break;
2509     }
2510   }  // end sector loop
2511   return isNearSector;
2512
2513 }  // end 'IsNearSector(float)'
2514
2515 // end ------------------------------------------------------------------------