SCorrelatorQAMaker/macros/QuickWeirdTuplePlotter.cxx

0001 // ----------------------------------------------------------------------------
0002 // 'QuickWeirdTuplePlotter.cxx
0003 // Derek Anderson
0004 // 10.26.2023
0005 //
0006 // Use to quick plot some leaves from the weird tracks ntuple.
0007 // ----------------------------------------------------------------------------
0008
0009 #include <string>
0010 #include <vector>
0011 #include <utility>
0012 #include <iostream>
0013 #include <TH2.h>
0014 #include <TFile.h>
0015 #include <TError.h>
0016 #include <TNtuple.h>
0017 #include <TDirectory.h>
0018 #include <ROOT/RDataFrame.hxx>
0019 #include <ROOT/RDF/HistoModels.hxx>
0020
0021 // make common namespaces implicit
0022 using namespace std;
0023 using namespace ROOT;
0024
0025 // set up aliases
0026 using TH2Def = ROOT::RDF::TH2DModel;
0027
0028
0029
0030 // quickly plot variables vs. delta-r froma tuple -----------------------------
0031
0032 void QuickWeirdTuplePlotter() {
0033
0034   // lower verbosity
0035   gErrorIgnoreLevel = kError;
0036   cout << "\n  Beginning weird track tuple plotting..." << endl;
0037
0038   // options, histogram binning, and histogram definitions --------------------
0039
0040   // i/o parameters
0041   const string sOutput  = "weirdTracksVsDr.pp200py8jet10run8.d27m10y2023.root";
0042   const string sInput   = "output/condor/final_merge/correlatorJetTree.pp200py8jet10run8_roughCutsWithWeirdTrackTuple.d27m10y2023.root";
0043   const string sInTuple = "ntWeirdTracks";
0044
0045   // cuts to apply, drawing options & delta-r leaf
0046   const string sCut("");
0047   const string sDeltaR("deltartrack");
0048
0049   // deltra-r binning
0050   const size_t         nDrBins(75);
0051   const vector<double> drBinEdges = {
0052     1e-05,
0053     1.16591e-05,
0054     1.35936e-05,
0055     1.58489e-05,
0056     1.84785e-05,
0057     2.15443e-05,
0058     2.51189e-05,
0059     2.92864e-05,
0060     3.41455e-05,
0061     3.98107e-05,
0062     4.64159e-05,
0063     5.4117e-05,
0064     6.30957e-05,
0065     7.35642e-05,
0066     8.57696e-05,
0067     0.0001,
0068     0.000116591,
0069     0.000135936,
0070     0.000158489,
0071     0.000184785,
0072     0.000215443,
0073     0.000251189,
0074     0.000292864,
0075     0.000341455,
0076     0.000398107,
0077     0.000464159,
0078     0.00054117,
0079     0.000630957,
0080     0.000735642,
0081     0.000857696,
0082     0.001,
0083     0.00116591,
0084     0.00135936,
0085     0.00158489,
0086     0.00184785,
0087     0.00215443,
0088     0.00251189,
0089     0.00292864,
0090     0.00341455,
0091     0.00398107,
0092     0.00464159,
0093     0.0054117,
0094     0.00630957,
0095     0.00735642,
0096     0.00857696,
0097     0.01,
0098     0.0116591,
0099     0.0135936,
0100     0.0158489,
0101     0.0184785,
0102     0.0215443,
0103     0.0251189,
0104     0.0292864,
0105     0.0341455,
0106     0.0398107,
0107     0.0464159,
0108     0.054117,
0109     0.0630957,
0110     0.0735642,
0111     0.0857696,
0112     0.1,
0113     0.116591,
0114     0.135936,
0115     0.158489,
0116     0.184785,
0117     0.215443,
0118     0.251189,
0119     0.292864,
0120     0.341455,
0121     0.398107,
0122     0.464159,
0123     0.54117,
0124     0.630957,
0125     0.735642,
0126     0.857696,
0127     1
0128   };  // end delta-r bin edges
0129
0130   // binning accessor
0131   enum Var {VNum, VEne, VEta, VPhi, VDca, VVtx, VQual, VDPt, VFrac};
0132
0133   // other variable binning
0134   //   <0> = no. of bins
0135   //   <1> = 1st bin lower edge
0136   //   <2> = last bin upper edge
0137   vector<tuple<size_t, double, double>> vecBins = {
0138     make_tuple(100,    0.,   100.),
0139     make_tuple(200,    0.,   100.),
0140     make_tuple(400,   -2.,   2.),
0141     make_tuple(360,   -3.15, 3.15),
0142     make_tuple(10000, -10.,  10.),
0143     make_tuple(10000, -10.,  10.),
0144     make_tuple(100,    0.,   10.),
0145     make_tuple(5000,   0.,   5.),
0146     make_tuple(5000,   0.,   5.)
0147   };
0148
0149   // histogram accessor
0150   enum Hist {HPt, HPFrac, HEta, HPhi, HEne, HDcaXY, HDcaZ, HVtxX, HVtxY, HVtxZ, HQual, HDPt, HNLMvtx, HNLIntt, HNLTpc, HNCMvtx, HNCIntt, HNCTpc, HNKey, HNSame, HNFrac};
0151
0152   // histogram titles
0153   const vector<string> vecTitles = {
0154     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};p_{T}^{trk} [GeV]",
0155     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};p_{T,i}^{trk} / p_{T,j}^{trk}",
0156     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};#eta^{trk}",
0157     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};#varphi^{trk}",
0158     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};E_{trk}",
0159     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};DCA_{xy} [cm]",
0160     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};DCZ_{z} [cm]",
0161     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};v_{x} [cm]",
0162     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};v_{y} [cm]",
0163     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};v_{z} [cm]",
0164     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};#chi^{2}/ndf",
0165     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};#deltap_{T}^{trk}/p_{T}^{trk}",
0166     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};N_{layer}^{mvtx}",
0167     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};N_{layer}^{intt}",
0168     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};N_{layer}^{tpc}",
0169     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};N_{clust}^{mvtx}",
0170     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};N_{clust}^{intt}",
0171     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};N_{clust}^{tpc}",
0172     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};N_{clust}^{tot}",
0173     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};N_{clust}^{same}",
0174     ";#DeltaR = #sqrt{(#Delta#varphi)^{2} + (#Delta#eta)^{2})};N_{clust}^{same} / N_{clust}^{tot}"
0175   };
0176
0177   // leaves to draw vs. delta-r & histogram definitions
0178   //  <0> = leaf to draw
0179   //  <1> = histogram name
0180   //  <2> = axis titles
0181   //  <3> = y-axis binning
0182   const vector<tuple<string, string, string, tuple<size_t, double, double>>> vecHistDef = {
0183     make_tuple("pt_a",          "hPtVsDr",         vecTitles[Hist::HPt],     vecBins[Var::VEne]),
0184     make_tuple("ptFrac",        "hPtFracVsDr",     vecTitles[Hist::HPFrac],  vecBins[Var::VFrac]),
0185     make_tuple("eta_a",         "hEtaVsDr",        vecTitles[Hist::HEta],    vecBins[Var::VEta]),
0186     make_tuple("phi_a",         "hPhiVsDr",        vecTitles[Hist::HPhi],    vecBins[Var::VPhi]),
0187     make_tuple("ene_a",         "hEneVsDr",        vecTitles[Hist::HEne],    vecBins[Var::VEne]),
0188     make_tuple("dcaxy_a",       "hDcaXYVsDr",      vecTitles[Hist::HDcaXY],  vecBins[Var::VDca]),
0189     make_tuple("dcaz_a",        "hDcaZVsDr",       vecTitles[Hist::HDcaZ],   vecBins[Var::VDca]),
0190     make_tuple("vtxx_a",        "hVtxXVsDr",       vecTitles[Hist::HVtxX],   vecBins[Var::VVtx]),
0191     make_tuple("vtxy_a",        "hVtxYVsDr",       vecTitles[Hist::HVtxY],   vecBins[Var::VVtx]),
0192     make_tuple("vtxz_a",        "hVtxZVsDr",       vecTitles[Hist::HVtxZ],   vecBins[Var::VVtx]),
0193     make_tuple("quality_a",     "hQualityVsDr",    vecTitles[Hist::HQual],   vecBins[Var::VQual]),
0194     make_tuple("deltapt_a",     "hDeltaPtVsDr",    vecTitles[Hist::HDPt],    vecBins[Var::VDPt]),
0195     make_tuple("nmvtxlayers_a", "hNMvtxLayerVsDr", vecTitles[Hist::HNLMvtx], vecBins[Var::VNum]),
0196     make_tuple("ninttlayers_a", "hNInttLayerVsDr", vecTitles[Hist::HNLIntt], vecBins[Var::VNum]),
0197     make_tuple("ntpclayers_a",  "hNTpcLayerVsDr",  vecTitles[Hist::HNLTpc],  vecBins[Var::VNum]),
0198     make_tuple("nmvtxclusts_a", "hNMvtxClustVsDr", vecTitles[Hist::HNCMvtx], vecBins[Var::VNum]),
0199     make_tuple("ninttclusts_a", "hNInttClustVsDr", vecTitles[Hist::HNCIntt], vecBins[Var::VNum]),
0200     make_tuple("ntpcclusts_a",  "hNTpcClustVsDr",  vecTitles[Hist::HNCTpc],  vecBins[Var::VNum]),
0201     make_tuple("nclustkey_a",   "hNClustKeyVsDr",  vecTitles[Hist::HNKey],   vecBins[Var::VNum]),
0202     make_tuple("nsameclustkey", "hNSameKeyVsDr",   vecTitles[Hist::HNSame],  vecBins[Var::VNum]),
0203     make_tuple("nSameFrac",     "hNSameFracVsDr",  vecTitles[Hist::HNFrac],  vecBins[Var::VFrac])
0204   };
0205   cout << "    Defined output histograms." << endl;
0206
0207   // create output histograms & TNtuple arguments
0208   //   <0> = what to draw on y-axis
0209   //   <1> = what to draw on x-axis
0210   //   <2> = pointer to histogram
0211   //vector<tuple<string, TH2D*>> vecArgAndHist;
0212   vector<tuple<string, string, TH2Def>> vecArgAndHist;
0213   for (const auto histDef : vecHistDef) {
0214
0215     // make histogram
0216     auto   bins = get<3>(histDef);
0217     TH2Def hist = TH2Def(get<1>(histDef).data(), get<2>(histDef).data(), nDrBins, drBinEdges.data(), get<0>(bins), get<1>(bins), get<2>(bins));
0218
0219     // add argument and hist to vector
0220     vecArgAndHist.push_back(make_tuple(get<0>(histDef), sDeltaR, hist));
0221   }
0222   cout << "    Created output histograms and TNtuple arguments." << endl;
0223
0224   // helper lambdas -----------------------------------------------------------
0225
0226   auto getFrac = [](const float a, const float b) {
0227     return a / b;
0228   };
0229
0230   // open files and make frame ------------------------------------------------
0231
0232   // open output file
0233   TFile* fOutput = new TFile(sOutput.data(), "recreate");
0234   if (!fOutput) {
0235     cerr << "PANIC: couldn't open a file!\n"
0236          << "       fOutput = " << fOutput <<  "\n"
0237          << endl;
0238     return;
0239   }
0240   cout << "    Opened output file." << endl;
0241
0242   // create data frame
0243   RDataFrame frame(sInTuple.data(), sInput.data());
0244
0245   // make sure file isn't empty
0246   auto pairs = frame.Count();
0247   if (pairs == 0) {
0248     cerr << "Error: No track pairs found!" << endl;
0249     return;
0250   }
0251   cout << "    Set up data frame." << endl;
0252
0253   // run any analysis and make histograms -------------------------------------
0254
0255   // define additional columns
0256   auto analysis = frame.Define("ptFrac",    getFrac, {"pt_a", "pt_b"})
0257                        .Define("nSameFrac", getFrac, {"nsameclustkey", "nclustkey_a"});
0258   cout << "    Defined additional columns.\n"
0259        << "    Beginning draw loop..."
0260        << endl;
0261
0262   // draw histograms & save to output
0263   for (auto argAndHist : vecArgAndHist) {
0264
0265     // draw histogram
0266     auto hist = analysis.Histo2D(get<2>(argAndHist), get<1>(argAndHist), get<0>(argAndHist));
0267     cout << "      Drawing '" << get<0>(argAndHist) << "'..." << endl;
0268
0269     // save to output file
0270     fOutput -> cd();
0271     hist    -> Write();
0272   }
0273   cout << "    Finished drawing and saved all histograms." << endl;
0274
0275   // close output and exit ----------------------------------------------------
0276
0277   // close files
0278   fOutput -> cd();
0279   fOutput -> Close();
0280   cout << "    Closed files." << endl;
0281
0282   // announce end and exit
0283   cout << "  Finished weird track tuple plotting!\n" << endl;
0284   return;
0285
0286 }
0287
0288 // end ------------------------------------------------------------------------