Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
CalculateSigmaDca.cxx
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file CalculateSigmaDca.cxx
1 // ----------------------------------------------------------------------------
2 // 'CalculateSigmaDca.cxx'
3 // Derek Anderson
4 // 10.11.2023
5 //
6 // Quick script to fit sigma-dca/dca distributions
7 // to figure out pt-dependent cut from SCorrelatorJetTree
8 // track QA tuples.
9 
10 // c utilities
11 #include <string>
12 #include <vector>
13 #include <utility>
14 #include <iostream>
15 // root classes
16 #include <TH1.h>
17 #include <TH2.h>
18 #include <TF1.h>
19 #include <TFile.h>
20 #include <TMath.h>
21 #include <TError.h>
22 #include <TNtuple.h>
23 #include <TString.h>
24 #include <TCanvas.h>
25 #include <TObjArray.h>
26 #include <TPaveText.h>
27 #include <TGraphErrors.h>
28 
29 using namespace std;
30 
31 // global constants
32 static const size_t NDca(2);
33 static const size_t NPar(4);
34 static const size_t NObj(2);
35 static const size_t NVtx(4);
36 
37 
38 
39 void CalculateSigmaDca() {
40 
41  // lower verbosity
42  gErrorIgnoreLevel = kError;
43  cout << "\n Starting sigma-dca/dca calculation..." << endl;
44 
45  // options ------------------------------------------------------------------
46 
47  // io parameters
48  const string sInput("../SCorrelatorJetTree/output/condor/final_merge/correlatorJetTree.pp200py8jet10run8_forVtxCutTest_onlyPrimTrks.d22m10y2023.root");
49  const string sTuple("QA/Tracks/ntTrkQA");
50  const string sOutput("dcaSigmaCalc.pp200py8jet10run8_smallSample_onlyPrimTrks_withCubicFitAndPtRange.d22m10y2023.root");
51 
52  // histogram parameters
53  const array<string, NDca> sDcaVsPtAll = {"hDcaXYvsPtAll", "hDcaZvsPtAll"};
54  const array<string, NDca> sDcaVsPtSel = {"hDcaXYvsPtSel", "hDcaZvsPtSel"};
55  const array<string, NDca> sDcaName = {"hDcaXY", "hDcaZ"};
56  const array<string, NDca> sWidthName = {"hWidthDcaXY", "hWidthDcaZ"};
57  const array<string, NDca> sTitleWidth = {"#DeltaDCA_{xy} [cm]", "#DeltaDCA_{z} [cm]"};
58  const array<string, NDca> sTitleDca = {"DCA_{xy} [cm]", "DCA_{z} [cm]"};
59  const string sTitlePt("p_{T}^{trk} [GeV/c]");
60  const string sDcaXYvsZAll("hDcaXYvsZAll");
61  const string sDcaXYvsZSel("hDcaXYvsZSel");
62 
63  // fit parameters
64  const uint16_t cut(0);
65  const string sDcaFit("gaus(0)");
66  const string sWidthFit("[0]+[1]/x+[2]/(x*x)+[3]/(x*x*x)");
67  const string sDcaOpt("QNR");
68  const string sWidthOpt("R");
69  const string sParam("_2");
70  const pair<float, float> ptFitRange = {1., 15.};
71  const array<string, NDca> sDcaFitName = {"fFitDcaXY", "fFitDcaZ"};
72  const array<string, NDca> sWidthFitName = {"fFitWidthXY", "fFitWidthZ"};
73  const array<float, NPar> fWidthGuess = {1., 1., 1., 1.};
74 
75  // cut parameters
76  const bool doOtherCuts(true);
77  const float nMvtxMin(2.);
78  const float nInttMin(1.);
79  const float nTpcMin(24.);
80  const float qualMax(10.);
81  const float etaMax(1.1);
82  const float ptMin(0.1);
83  const float nCut(3.);
84  const string sAll("all tracks");
85  const string sSel("selected tracks");
86 
87  // plot options
88  const pair<float, float> ptPlotRange = {0., 25.};
89  const pair<float, float> dcaPlotRange = {-0.5, 0.5};
90  const pair<float, float> widthPlotRange = {0.001, 0.05};
91  const array<string, NDca> sDcaWidthPlot = {"cWidthDcaXY", "cWidthDcaZ"};
92  const array<string, NDca> sDcaVsPtPlot = {"cDcaXYvsPt", "cDcaZvsPt"};
93  const string sDcaXYvsZPlot("cDcaXYvsZ");
94 
95  // text parameters
96  const vector<string> sText = {
97  "#bf{#it{sPHENIX}} simulation [Run 8]",
98  "PYTHIA-8, JS 10 GeV sample",
99  "#bf{reco. tracks}"
100  };
101 
102  // open files & initialize tuple --------------------------------------------
103 
104  // open file
105  TFile* fOutput = new TFile(sOutput.data(), "recreate");
106  TFile* fInput = new TFile(sInput.data(), "read");
107  if (!fOutput || !fInput) {
108  cerr << "PANIC: couldn't open a file!\n"
109  << " fOutput = " << fOutput << ", " << fInput << "\n"
110  << endl;
111  return;
112  }
113  cout << " Opened files." << endl;
114 
115  // grab input tuple
116  TNtuple* ntInput = (TNtuple*) fInput -> Get(sTuple.data());
117  if (!ntInput) {
118  cerr << "PANIC: couldn't grab input tuple!\n" << endl;
119  return;
120  }
121  cout << " Grabbed input tuple." << endl;
122 
123  // declare leaves
124  float pt;
125  float eta;
126  float phi;
127  float energy;
128  float dcaxy;
129  float dcaz;
130  float deltapt;
131  float quality;
132  float nmvtxlayer;
133  float ninttlayer;
134  float ntpclayer;
135 
136  // set branch addresses
137  ntInput -> SetBranchAddress("pt", &pt);
138  ntInput -> SetBranchAddress("eta", &eta);
139  ntInput -> SetBranchAddress("phi", &phi);
140  ntInput -> SetBranchAddress("energy", &energy);
141  ntInput -> SetBranchAddress("dcaxy", &dcaxy);
142  ntInput -> SetBranchAddress("dcaz", &dcaz);
143  ntInput -> SetBranchAddress("deltapt", &deltapt);
144  ntInput -> SetBranchAddress("quality", &quality);
145  ntInput -> SetBranchAddress("nmvtxlayer",&nmvtxlayer);
146  ntInput -> SetBranchAddress("ninttlayer",&ninttlayer);
147  ntInput -> SetBranchAddress("ntpclayer", &ntpclayer);
148  cout << " Set input tuple branches." << endl;
149 
150  // initialize histograms ----------------------------------------------------
151 
152  // histogram binning
153  const tuple<size_t, pair<float, float>> dcaBins = {2000, make_pair(-5., 5.)};
154  const tuple<size_t, pair<float, float>> ptBins = {100, make_pair(0., 100.)};
155 
156  // declare histograms
157  array<TH1D*, NDca> arrDcaWidth;
158  array<TH2D*, NDca> arrDcaVsPtAll;
159  array<TH2D*, NDca> arrDcaVsPtSel;
160  for (size_t iDca = 0; iDca < NDca; iDca++) {
161  arrDcaWidth[iDca] = new TH1D(
162  sWidthName[iDca].data(),
163  "",
164  get<0>(ptBins),
165  get<1>(ptBins).first,
166  get<1>(ptBins).second
167  );
168  arrDcaVsPtAll[iDca] = new TH2D(
169  sDcaVsPtAll[iDca].data(),
170  sAll.data(),
171  get<0>(ptBins),
172  get<1>(ptBins).first,
173  get<1>(ptBins).second,
174  get<0>(dcaBins),
175  get<1>(dcaBins).first,
176  get<1>(dcaBins).second
177  );
178  arrDcaVsPtSel[iDca] = new TH2D(
179  sDcaVsPtSel[iDca].data(),
180  sSel.data(),
181  get<0>(ptBins),
182  get<1>(ptBins).first,
183  get<1>(ptBins).second,
184  get<0>(dcaBins),
185  get<1>(dcaBins).first,
186  get<1>(dcaBins).second
187  );
188  }
189 
190  TH2D* hDcaXYvsZAll = new TH2D(
191  sDcaXYvsZAll.data(),
192  sAll.data(),
193  get<0>(dcaBins),
194  get<1>(dcaBins).first,
195  get<1>(dcaBins).second,
196  get<0>(dcaBins),
197  get<1>(dcaBins).first,
198  get<1>(dcaBins).second
199  );
200  TH2D* hDcaXYvsZSel = new TH2D(
201  sDcaXYvsZSel.data(),
202  sSel.data(),
203  get<0>(dcaBins),
204  get<1>(dcaBins).first,
205  get<1>(dcaBins).second,
206  get<0>(dcaBins),
207  get<1>(dcaBins).first,
208  get<1>(dcaBins).second
209  );
210  cout << " Initialized output histograms." << endl;
211 
212  // 1st entry loop -------------------------------------------------------------
213 
214  // start entry loop
215  uint64_t nEntries = ntInput -> GetEntries();
216  cout << " Starting 1st entry loop: " << nEntries << " to process" << endl;
217 
218  uint64_t nBytes = 0;
219  for (uint64_t iEntry = 0; iEntry < nEntries; iEntry++) {
220 
221  // grab entry
222  const uint64_t bytes = ntInput -> GetEntry(iEntry);
223  if (bytes < 0) {
224  cerr << "WARNING: issue with entry " << iEntry << "! Aborting loop!" << endl;
225  break;
226  } else {
227  nBytes += bytes;
228  }
229 
230  // announce progress
231  const uint64_t iProg = iEntry + 1;
232  if (iProg == nEntries) {
233  cout << " Processing entry " << iProg << "/" << nEntries << "..." << endl;
234  } else {
235  cout << " Processing entry " << iProg << "/" << nEntries << "...\r" << flush;
236  }
237 
238  // apply cuts other than dca if needed
239  if (doOtherCuts) {
240  const bool isNumMvtxGood = (nmvtxlayer > nMvtxMin);
241  const bool isNumInttGood = (ninttlayer > nInttMin);
242  const bool isNumTpcGood = (ntpclayer > nTpcMin);
243  const bool isQualityGood = (quality < qualMax);
244  const bool isEtaGood = (abs(eta) < etaMax);
245  const bool isPtGood = (pt > ptMin);
246  const bool isTrackGood = (isNumMvtxGood && isNumInttGood && isNumTpcGood && isQualityGood && isEtaGood && isPtGood);
247  if (!isTrackGood) continue;
248  }
249 
250  // fill histograms
251  arrDcaVsPtAll[0] -> Fill(pt, dcaxy);
252  arrDcaVsPtAll[1] -> Fill(pt, dcaz);
253  hDcaXYvsZAll -> Fill(dcaz, dcaxy);
254 
255  } // end entry loop 1
256  cout << " Finished 1st entry loop!" << endl;
257 
258  // get & fit widths ---------------------------------------------------------
259 
260  array<TF1*, NDca> arrWidthFits;
261  for (size_t iDca = 0; iDca < NDca; iDca++) {
262 
263  // create result name
264  TString sResult(sDcaVsPtAll[iDca].data());
265  sResult.Append(sParam.data());
266 
267  // determine dca range
268  const uint64_t iStart = arrDcaVsPtAll[iDca] -> FindFirstBinAbove(0., 1);
269  const uint64_t iStop = arrDcaVsPtAll[iDca] -> FindLastBinAbove(0., 1);
270  const float start = arrDcaVsPtAll[iDca] -> GetXaxis() -> GetBinLowEdge(iStart);
271  const float stop = arrDcaVsPtAll[iDca] -> GetXaxis() -> GetBinLowEdge(iStop + 1);
272 
273  // declare fit and parameter pointers
274  TF1* fDcaFit = new TF1(sDcaFitName[iDca].data(), sDcaFit.data(), start, stop);
275  TObjArray* arrParams = NULL;
276 
277  // fit slices of dca vs pt
278  arrDcaVsPtAll[iDca] -> FitSlicesY(fDcaFit, iStart, iStop, cut, sDcaOpt.data(), arrParams);
279  arrDcaWidth[iDca] = (TH1D*) gDirectory -> Get(sResult.Data()) -> Clone();
280 
281  // declare width fit function
282  arrWidthFits[iDca] = new TF1(sWidthFitName[iDca].data(), sWidthFit.data(), start, stop);
283  for (size_t iPar = 0; iPar < NPar; iPar++) {
284  arrWidthFits[iDca] -> SetParameters(iPar, fWidthGuess[iPar]);
285  }
286 
287  // fit width of dca distributions
288  arrDcaWidth[iDca] -> Fit(arrWidthFits[iDca], sWidthOpt.data(), "", ptFitRange.first, ptFitRange.second);
289  arrDcaWidth[iDca] -> SetName(sWidthName[iDca].data());
290  } // end dca variable loop
291 
292  // 2nd entry loop -------------------------------------------------------------
293 
294  // start entry loop
295  cout << " Starting 2nd entry loop: applying a cut of " << nCut << " times the DCA width" << endl;
296 
297  nBytes = 0;
298  for (uint64_t iEntry = 0; iEntry < nEntries; iEntry++) {
299 
300  // grab entry
301  const uint64_t bytes = ntInput -> GetEntry(iEntry);
302  if (bytes < 0) {
303  cerr << "WARNING: issue with entry " << iEntry << "! Aborting loop!" << endl;
304  break;
305  } else {
306  nBytes += bytes;
307  }
308 
309  // announce progress
310  const uint64_t iProg = iEntry + 1;
311  if (iProg == nEntries) {
312  cout << " Processing entry " << iProg << "/" << nEntries << "..." << endl;
313  } else {
314  cout << " Processing entry " << iProg << "/" << nEntries << "...\r" << flush;
315  }
316 
317  // apply cuts other than dca if needed
318  if (doOtherCuts) {
319  const bool isNumMvtxGood = (nmvtxlayer > nMvtxMin);
320  const bool isNumInttGood = (ninttlayer > nInttMin);
321  const bool isNumTpcGood = (ntpclayer > nTpcMin);
322  const bool isQualityGood = (quality < qualMax);
323  const bool isEtaGood = (abs(eta) < etaMax);
324  const bool isPtGood = (pt > ptMin);
325  const bool isTrackGood = (isNumMvtxGood && isNumInttGood && isNumTpcGood && isQualityGood && isEtaGood && isPtGood);
326  if (!isTrackGood) continue;
327  }
328 
329  // calculate max dca
330  const double dcaWidthXY = arrWidthFits[0] -> Eval(pt);
331  const double dcaWidthZ = arrWidthFits[1] -> Eval(pt);
332  const double dcaMaxXY = nCut * dcaWidthXY;
333  const double dcaMaxZ = nCut * dcaWidthZ;
334 
335  // apply cuts
336  const bool isTrkInDcaCutXY = (abs(dcaxy) < dcaMaxXY);
337  const bool isTrkInDcaCutZ = (abs(dcaz) < dcaMaxZ);
338  const bool isTrkInDcaCut = (isTrkInDcaCutXY && isTrkInDcaCutZ);
339  if (!isTrkInDcaCut) continue;
340 
341  // fill histograms
342  arrDcaVsPtSel[0] -> Fill(pt, dcaxy);
343  arrDcaVsPtSel[1] -> Fill(pt, dcaz);
344  hDcaXYvsZSel -> Fill(dcaz, dcaxy);
345 
346  } // end entry loop 1
347  cout << " Finished 2nd entry loop!" << endl;
348 
349  // make cut lines -----------------------------------------------------------
350 
351  // line options
352  const uint16_t fColLin(2);
353  const uint16_t fStyLin(1);
354  const uint16_t fWidLin(2);
355 
356  array<TF1*, NDca> arrFitsNeg;
357  array<TF1*, NDca> arrFitsPos;
358  array<float, NPar> arrNegParams;
359  array<float, NPar> arrPosParams;
360  for (size_t iDca = 0; iDca < NDca; iDca++) {
361 
362  // make name
363  const TString sFitBase = arrWidthFits[iDca] -> GetName();
364  TString sFitNeg = sFitBase;
365  TString sFitPos = sFitBase;
366  sFitNeg.Append("_Neg");
367  sFitPos.Append("_Pos");
368 
369  // set parameters
370  for (size_t iPar = 0; iPar < NPar; iPar++) {
371  arrNegParams[iPar] = arrWidthFits[iDca] -> GetParameter(iPar);
372  arrPosParams[iPar] = arrWidthFits[iDca] -> GetParameter(iPar);
373  arrNegParams[iPar] *= (-1. * nCut);
374  arrPosParams[iPar] *= nCut;
375  }
376 
377  // make functions
378  arrFitsNeg[iDca] = (TF1*) arrWidthFits[iDca] -> Clone();
379  arrFitsPos[iDca] = (TF1*) arrWidthFits[iDca] -> Clone();
380  arrFitsNeg[iDca] -> SetName(sFitNeg.Data());
381  arrFitsPos[iDca] -> SetName(sFitPos.Data());
382  arrFitsNeg[iDca] -> SetLineColor(fColLin);
383  arrFitsPos[iDca] -> SetLineColor(fColLin);
384  arrFitsNeg[iDca] -> SetLineStyle(fStyLin);
385  arrFitsPos[iDca] -> SetLineStyle(fStyLin);
386  arrFitsNeg[iDca] -> SetLineWidth(fWidLin);
387  arrFitsPos[iDca] -> SetLineWidth(fWidLin);
388  for (size_t iPar = 0; iPar < NPar; iPar++) {
389  arrFitsNeg[iDca] -> SetParameter(iPar, arrNegParams[iPar]);
390  arrFitsPos[iDca] -> SetParameter(iPar, arrPosParams[iPar]);
391  }
392  } // end variable loops
393  cout << " Made cut lines." << endl;
394 
395  // make plots ---------------------------------------------------------------
396 
397  // hist options
398  const uint16_t fCol(1);
399  const uint16_t fMar(20);
400  const uint16_t fLin(1);
401  const uint16_t fWid(1);
402  const uint16_t fFil(0);
403  const uint16_t fTxt(42);
404  const uint16_t fCnt(1);
405  const float fOffX(1.1);
406  const float fOffY(1.3);
407  const float fTtl(0.04);
408  const float fLbl(0.03);
409 
410  // set dca-specific hist options
411  for (size_t iDca = 0; iDca < NDca; iDca++) {
412 
413  // figure out dca vs. pt z-axis ranges
414  const float minVsPtAll = arrDcaVsPtAll[iDca] -> GetMinimum(0.);
415  const float minVsPtSel = arrDcaVsPtSel[iDca] -> GetMinimum(0.);
416  const float maxVsPtAll = arrDcaVsPtAll[iDca] -> GetMaximum();
417  const float maxVsPtSel = arrDcaVsPtSel[iDca] -> GetMaximum();
418  const float minVsPtZ = TMath::Min(minVsPtAll, minVsPtSel);
419  const float maxVsPtZ = TMath::Max(maxVsPtAll, maxVsPtSel);
420 
421  // set 1d hist options
422  arrDcaWidth[iDca] -> SetMarkerColor(fCol);
423  arrDcaWidth[iDca] -> SetMarkerStyle(fMar);
424  arrDcaWidth[iDca] -> SetLineColor(fCol);
425  arrDcaWidth[iDca] -> SetLineStyle(fLin);
426  arrDcaWidth[iDca] -> SetLineWidth(fWid);
427  arrDcaWidth[iDca] -> SetFillColor(fCol);
428  arrDcaWidth[iDca] -> SetFillStyle(fFil);
429  arrDcaWidth[iDca] -> SetTitleFont(fTxt);
430  arrDcaWidth[iDca] -> SetTitle("");
431  arrDcaWidth[iDca] -> GetXaxis() -> SetRangeUser(ptPlotRange.first, ptPlotRange.second);
432  arrDcaWidth[iDca] -> GetXaxis() -> SetLabelSize(fLbl);
433  arrDcaWidth[iDca] -> GetXaxis() -> SetTitle(sTitlePt.data());
434  arrDcaWidth[iDca] -> GetXaxis() -> SetTitleFont(fTxt);
435  arrDcaWidth[iDca] -> GetXaxis() -> SetTitleSize(fTtl);
436  arrDcaWidth[iDca] -> GetXaxis() -> SetTitleOffset(fOffX);
437  arrDcaWidth[iDca] -> GetXaxis() -> CenterTitle(fCnt);
438  arrDcaWidth[iDca] -> GetYaxis() -> SetRangeUser(widthPlotRange.first, widthPlotRange.second);
439  arrDcaWidth[iDca] -> GetYaxis() -> SetLabelSize(fLbl);
440  arrDcaWidth[iDca] -> GetYaxis() -> SetTitle(sTitleWidth[iDca].data());
441  arrDcaWidth[iDca] -> GetYaxis() -> SetTitleFont(fTxt);
442  arrDcaWidth[iDca] -> GetYaxis() -> SetTitleSize(fTtl);
443  arrDcaWidth[iDca] -> GetYaxis() -> SetTitleOffset(fOffY);
444  arrDcaWidth[iDca] -> GetYaxis() -> CenterTitle(fCnt);
445 
446  // set dca vs. pt hist options
447  arrDcaVsPtAll[iDca] -> SetMarkerColor(fCol);
448  arrDcaVsPtAll[iDca] -> SetMarkerStyle(fMar);
449  arrDcaVsPtAll[iDca] -> SetLineColor(fCol);
450  arrDcaVsPtAll[iDca] -> SetLineStyle(fLin);
451  arrDcaVsPtAll[iDca] -> SetLineWidth(fWid);
452  arrDcaVsPtAll[iDca] -> SetFillColor(fCol);
453  arrDcaVsPtAll[iDca] -> SetFillStyle(fFil);
454  arrDcaVsPtAll[iDca] -> SetTitleFont(fTxt);
455  arrDcaVsPtAll[iDca] -> GetXaxis() -> SetRangeUser(ptPlotRange.first, ptPlotRange.second);
456  arrDcaVsPtAll[iDca] -> GetXaxis() -> SetLabelSize(fLbl);
457  arrDcaVsPtAll[iDca] -> GetXaxis() -> SetTitle(sTitlePt.data());
458  arrDcaVsPtAll[iDca] -> GetXaxis() -> SetTitleFont(fTxt);
459  arrDcaVsPtAll[iDca] -> GetXaxis() -> SetTitleSize(fTtl);
460  arrDcaVsPtAll[iDca] -> GetXaxis() -> SetTitleOffset(fOffX);
461  arrDcaVsPtAll[iDca] -> GetXaxis() -> CenterTitle(fCnt);
462  arrDcaVsPtAll[iDca] -> GetYaxis() -> SetRangeUser(dcaPlotRange.first, dcaPlotRange.second);
463  arrDcaVsPtAll[iDca] -> GetYaxis() -> SetLabelSize(fLbl);
464  arrDcaVsPtAll[iDca] -> GetYaxis() -> SetTitle(sTitleDca[iDca].data());
465  arrDcaVsPtAll[iDca] -> GetYaxis() -> SetTitleFont(fTxt);
466  arrDcaVsPtAll[iDca] -> GetYaxis() -> SetTitleSize(fTtl);
467  arrDcaVsPtAll[iDca] -> GetYaxis() -> SetTitleOffset(fOffY);
468  arrDcaVsPtAll[iDca] -> GetYaxis() -> CenterTitle(fCnt);
469  arrDcaVsPtAll[iDca] -> GetZaxis() -> SetRangeUser(minVsPtZ, maxVsPtZ);
470  arrDcaVsPtAll[iDca] -> GetZaxis() -> SetLabelSize(fLbl);
471  arrDcaVsPtSel[iDca] -> SetMarkerColor(fCol);
472  arrDcaVsPtSel[iDca] -> SetMarkerStyle(fMar);
473  arrDcaVsPtSel[iDca] -> SetLineColor(fCol);
474  arrDcaVsPtSel[iDca] -> SetLineStyle(fLin);
475  arrDcaVsPtSel[iDca] -> SetLineWidth(fWid);
476  arrDcaVsPtSel[iDca] -> SetFillColor(fCol);
477  arrDcaVsPtSel[iDca] -> SetFillStyle(fFil);
478  arrDcaVsPtSel[iDca] -> SetTitleFont(fTxt);
479  arrDcaVsPtSel[iDca] -> GetXaxis() -> SetRangeUser(ptPlotRange.first, ptPlotRange.second);
480  arrDcaVsPtSel[iDca] -> GetXaxis() -> SetLabelSize(fLbl);
481  arrDcaVsPtSel[iDca] -> GetXaxis() -> SetTitle(sTitlePt.data());
482  arrDcaVsPtSel[iDca] -> GetXaxis() -> SetTitleFont(fTxt);
483  arrDcaVsPtSel[iDca] -> GetXaxis() -> SetTitleSize(fTtl);
484  arrDcaVsPtSel[iDca] -> GetXaxis() -> SetTitleOffset(fOffX);
485  arrDcaVsPtSel[iDca] -> GetXaxis() -> CenterTitle(fCnt);
486  arrDcaVsPtSel[iDca] -> GetYaxis() -> SetRangeUser(dcaPlotRange.first, dcaPlotRange.second);
487  arrDcaVsPtSel[iDca] -> GetYaxis() -> SetLabelSize(fLbl);
488  arrDcaVsPtSel[iDca] -> GetYaxis() -> SetTitle(sTitleDca[iDca].data());
489  arrDcaVsPtSel[iDca] -> GetYaxis() -> SetTitleFont(fTxt);
490  arrDcaVsPtSel[iDca] -> GetYaxis() -> SetTitleSize(fTtl);
491  arrDcaVsPtSel[iDca] -> GetYaxis() -> SetTitleOffset(fOffY);
492  arrDcaVsPtSel[iDca] -> GetYaxis() -> CenterTitle(fCnt);
493  arrDcaVsPtSel[iDca] -> GetZaxis() -> SetRangeUser(minVsPtZ, maxVsPtZ);
494  arrDcaVsPtSel[iDca] -> GetZaxis() -> SetLabelSize(fLbl);
495  } // end variable loop
496 
497  // figure out dca xy vs. z z-axis ranges
498  const float minVsDcaAll = hDcaXYvsZAll -> GetMinimum(0.);
499  const float minVsDcaSel = hDcaXYvsZSel -> GetMinimum(0.);
500  const float maxVsDcaAll = hDcaXYvsZAll -> GetMaximum();
501  const float maxVsDcaSel = hDcaXYvsZSel -> GetMaximum();
502  const float minVsDcaZ = TMath::Min(minVsDcaAll, minVsDcaSel);
503  const float maxVsDcaZ = TMath::Max(maxVsDcaAll, maxVsDcaSel);
504 
505  // set dca xy vs. z hist options
506  hDcaXYvsZAll -> SetMarkerColor(fCol);
507  hDcaXYvsZAll -> SetMarkerStyle(fMar);
508  hDcaXYvsZAll -> SetLineColor(fCol);
509  hDcaXYvsZAll -> SetLineStyle(fLin);
510  hDcaXYvsZAll -> SetLineWidth(fWid);
511  hDcaXYvsZAll -> SetFillColor(fCol);
512  hDcaXYvsZAll -> SetFillStyle(fFil);
513  hDcaXYvsZAll -> SetTitleFont(fTxt);
514  hDcaXYvsZAll -> GetXaxis() -> SetRangeUser(dcaPlotRange.first, dcaPlotRange.second);
515  hDcaXYvsZAll -> GetXaxis() -> SetLabelSize(fLbl);
516  hDcaXYvsZAll -> GetXaxis() -> SetTitle(sTitleDca[1].data());
517  hDcaXYvsZAll -> GetXaxis() -> SetTitleFont(fTxt);
518  hDcaXYvsZAll -> GetXaxis() -> SetTitleSize(fTtl);
519  hDcaXYvsZAll -> GetXaxis() -> SetTitleOffset(fOffX);
520  hDcaXYvsZAll -> GetXaxis() -> CenterTitle(fCnt);
521  hDcaXYvsZAll -> GetYaxis() -> SetRangeUser(dcaPlotRange.first, dcaPlotRange.second);
522  hDcaXYvsZAll -> GetYaxis() -> SetLabelSize(fLbl);
523  hDcaXYvsZAll -> GetYaxis() -> SetTitle(sTitleDca[0].data());
524  hDcaXYvsZAll -> GetYaxis() -> SetTitleFont(fTxt);
525  hDcaXYvsZAll -> GetYaxis() -> SetTitleSize(fTtl);
526  hDcaXYvsZAll -> GetYaxis() -> SetTitleOffset(fOffY);
527  hDcaXYvsZAll -> GetYaxis() -> CenterTitle(fCnt);
528  hDcaXYvsZAll -> GetZaxis() -> SetRangeUser(minVsDcaZ, maxVsDcaZ);
529  hDcaXYvsZAll -> GetZaxis() -> SetLabelSize(fLbl);
530  hDcaXYvsZSel -> SetMarkerColor(fCol);
531  hDcaXYvsZSel -> SetMarkerStyle(fMar);
532  hDcaXYvsZSel -> SetLineColor(fCol);
533  hDcaXYvsZSel -> SetLineStyle(fLin);
534  hDcaXYvsZSel -> SetLineWidth(fWid);
535  hDcaXYvsZSel -> SetFillColor(fCol);
536  hDcaXYvsZSel -> SetFillStyle(fFil);
537  hDcaXYvsZSel -> SetTitleFont(fTxt);
538  hDcaXYvsZSel -> GetXaxis() -> SetRangeUser(dcaPlotRange.first, dcaPlotRange.second);
539  hDcaXYvsZSel -> GetXaxis() -> SetLabelSize(fLbl);
540  hDcaXYvsZSel -> GetXaxis() -> SetTitle(sTitleDca[1].data());
541  hDcaXYvsZSel -> GetXaxis() -> SetTitleFont(fTxt);
542  hDcaXYvsZSel -> GetXaxis() -> SetTitleSize(fTtl);
543  hDcaXYvsZSel -> GetXaxis() -> SetTitleOffset(fOffX);
544  hDcaXYvsZSel -> GetXaxis() -> CenterTitle(fCnt);
545  hDcaXYvsZSel -> GetYaxis() -> SetRangeUser(dcaPlotRange.first, dcaPlotRange.second);
546  hDcaXYvsZSel -> GetYaxis() -> SetLabelSize(fLbl);
547  hDcaXYvsZSel -> GetYaxis() -> SetTitle(sTitleDca[0].data());
548  hDcaXYvsZSel -> GetYaxis() -> SetTitleFont(fTxt);
549  hDcaXYvsZSel -> GetYaxis() -> SetTitleSize(fTtl);
550  hDcaXYvsZSel -> GetYaxis() -> SetTitleOffset(fOffY);
551  hDcaXYvsZSel -> GetYaxis() -> CenterTitle(fCnt);
552  hDcaXYvsZSel -> GetZaxis() -> SetRangeUser(minVsDcaZ, maxVsDcaZ);
553  hDcaXYvsZSel -> GetZaxis() -> SetLabelSize(fLbl);
554  cout << " Set histogram options." << endl;
555 
556  // text box options
557  const uint16_t fAln = 12;
558  const uint16_t fColTxt = 0;
559  const uint16_t fFilTxt = 0;
560  const uint16_t fLinTxt = 0;
561  const size_t nTxt = sText.size();
562  const float yTxtStart = 0.1;
563  const float yTxtStop = yTxtStart + (nTxt * 0.05);
564 
565  // text box dimensions
566  const array<float, NVtx> xyTxt = {0.1, yTxtStart, 0.3, yTxtStop};
567 
568  // create text box
569  TPaveText *ptInfo = new TPaveText(xyTxt[0], xyTxt[1], xyTxt[2], xyTxt[3], "NDC NB");
570  ptInfo -> SetFillColor(fColTxt);
571  ptInfo -> SetFillStyle(fFilTxt);
572  ptInfo -> SetLineColor(fColTxt);
573  ptInfo -> SetLineStyle(fLinTxt);
574  ptInfo -> SetTextFont(fTxt);
575  ptInfo -> SetTextAlign(fAln);
576  for (const string text : sText) {
577  ptInfo -> AddText(text.data());
578  }
579  cout << " Made text box." << endl;
580 
581  // plot options
582  const uint16_t fMode(0);
583  const uint16_t fBord(2);
584  const uint16_t fGrid(0);
585  const uint16_t fTick(1);
586  const uint16_t fLogZ(1);
587  const uint16_t fFrame(0);
588  const string sPadAll("pAll");
589  const string sPadSel("pSelect");
590 
591  // plot margins
592  const array<float, NVtx> xyMarginWidth = {0.02, 0.02, 0.15, 0.15};
593  const array<float, NVtx> xyMarginPad = {0.1, 0.15, 0.15, 0.15};
594 
595  // canvas & pad dimensions
596  const pair<uint32_t, uint32_t> widthDim = {750, 750};
597  const pair<uint32_t, uint32_t> allVsSelDim = {1500, 750};
598  const array<float, NVtx> xyAllPad = {0.0, 0.0, 0.5, 1.0};
599  const array<float, NVtx> xySelPad = {0.5, 0.0, 1.0, 1.0};
600 
601  // make dca-specific plots
602  array<TCanvas*, NDca> arrWidthPlots;
603  array<TCanvas*, NDca> arrDcaVsPtPlots;
604  array<TPad*, NDca> arrDcaVsPtAllPad;
605  array<TPad*, NDca> arrDcaVsPtSelPad;
606  for (size_t iDca = 0; iDca < NDca; iDca++) {
607 
608  // make width plots
609  arrWidthPlots[iDca] = new TCanvas(sDcaWidthPlot[iDca].data(), "", widthDim.first, widthDim.second);
610  arrWidthPlots[iDca] -> SetGrid(fGrid, fGrid);
611  arrWidthPlots[iDca] -> SetTicks(fTick, fTick);
612  arrWidthPlots[iDca] -> SetBorderMode(fMode);
613  arrWidthPlots[iDca] -> SetBorderSize(fBord);
614  arrWidthPlots[iDca] -> SetTopMargin(xyMarginWidth[0]);
615  arrWidthPlots[iDca] -> SetRightMargin(xyMarginWidth[1]);
616  arrWidthPlots[iDca] -> SetBottomMargin(xyMarginWidth[2]);
617  arrWidthPlots[iDca] -> SetLeftMargin(xyMarginWidth[3]);
618  arrWidthPlots[iDca] -> cd();
619  arrDcaWidth[iDca] -> Draw();
620  ptInfo -> Draw();
621  fOutput -> cd();
622  arrWidthPlots[iDca] -> Write();
623  arrWidthPlots[iDca] -> Close();
624 
625  // make dca vs. pt plots
626  arrDcaVsPtPlots[iDca] = new TCanvas(sDcaVsPtPlot[iDca].data(), "", allVsSelDim.first, allVsSelDim.second);
627  arrDcaVsPtAllPad[iDca] = new TPad(sPadAll.data(), "", xyAllPad[0], xyAllPad[1], xyAllPad[2], xyAllPad[3]);
628  arrDcaVsPtSelPad[iDca] = new TPad(sPadSel.data(), "", xySelPad[0], xySelPad[1], xySelPad[2], xySelPad[3]);
629  arrDcaVsPtPlots[iDca] -> SetGrid(fGrid, fGrid);
630  arrDcaVsPtPlots[iDca] -> SetTicks(fTick, fTick);
631  arrDcaVsPtPlots[iDca] -> SetBorderMode(fMode);
632  arrDcaVsPtPlots[iDca] -> SetBorderSize(fBord);
633  arrDcaVsPtAllPad[iDca] -> SetGrid(fGrid, fGrid);
634  arrDcaVsPtAllPad[iDca] -> SetTicks(fTick, fTick);
635  arrDcaVsPtAllPad[iDca] -> SetBorderMode(fMode);
636  arrDcaVsPtAllPad[iDca] -> SetBorderSize(fBord);
637  arrDcaVsPtAllPad[iDca] -> SetLogz(fLogZ);
638  arrDcaVsPtAllPad[iDca] -> SetTopMargin(xyMarginPad[0]);
639  arrDcaVsPtAllPad[iDca] -> SetRightMargin(xyMarginPad[1]);
640  arrDcaVsPtAllPad[iDca] -> SetBottomMargin(xyMarginPad[2]);
641  arrDcaVsPtAllPad[iDca] -> SetLeftMargin(xyMarginPad[3]);
642  arrDcaVsPtSelPad[iDca] -> SetGrid(fGrid, fGrid);
643  arrDcaVsPtSelPad[iDca] -> SetTicks(fTick, fTick);
644  arrDcaVsPtSelPad[iDca] -> SetBorderMode(fMode);
645  arrDcaVsPtSelPad[iDca] -> SetBorderSize(fBord);
646  arrDcaVsPtSelPad[iDca] -> SetLogz(fLogZ);
647  arrDcaVsPtSelPad[iDca] -> SetTopMargin(xyMarginPad[0]);
648  arrDcaVsPtSelPad[iDca] -> SetRightMargin(xyMarginPad[1]);
649  arrDcaVsPtSelPad[iDca] -> SetBottomMargin(xyMarginPad[2]);
650  arrDcaVsPtSelPad[iDca] -> SetLeftMargin(xyMarginPad[3]);
651  arrDcaVsPtPlots[iDca] -> cd();
652  arrDcaVsPtAllPad[iDca] -> Draw();
653  arrDcaVsPtSelPad[iDca] -> Draw();
654  arrDcaVsPtAllPad[iDca] -> cd();
655  arrDcaVsPtAll[iDca] -> Draw("colz");
656  arrFitsNeg[iDca] -> Draw("same");
657  arrFitsPos[iDca] -> Draw("same");
658  arrDcaVsPtSelPad[iDca] -> cd();
659  arrDcaVsPtSel[iDca] -> Draw("colz");
660  ptInfo -> Draw();
661  fOutput -> cd();
662  arrDcaVsPtPlots[iDca] -> Write();
663  arrDcaVsPtPlots[iDca] -> Close();
664  } // end variable loop
665 
666  // make dca xy vs. z plot
667  TCanvas* cDcaXYvsZ = new TCanvas(sDcaXYvsZPlot.data(), "", allVsSelDim.first, allVsSelDim.second);
668  TPad* pDcaXYvsZAll = new TPad(sPadAll.data(), "", xyAllPad[0], xyAllPad[1], xyAllPad[2], xyAllPad[3]);
669  TPad* pDcaXYvsZSel = new TPad(sPadSel.data(), "", xySelPad[0], xySelPad[1], xySelPad[2], xySelPad[3]);
670  cDcaXYvsZ -> SetGrid(fGrid, fGrid);
671  cDcaXYvsZ -> SetTicks(fTick, fTick);
672  cDcaXYvsZ -> SetBorderMode(fMode);
673  cDcaXYvsZ -> SetBorderSize(fBord);
674  pDcaXYvsZAll -> SetGrid(fGrid, fGrid);
675  pDcaXYvsZAll -> SetTicks(fTick, fTick);
676  pDcaXYvsZAll -> SetBorderMode(fMode);
677  pDcaXYvsZAll -> SetBorderSize(fBord);
678  pDcaXYvsZAll -> SetLogz(fLogZ);
679  pDcaXYvsZAll -> SetTopMargin(xyMarginPad[0]);
680  pDcaXYvsZAll -> SetRightMargin(xyMarginPad[1]);
681  pDcaXYvsZAll -> SetBottomMargin(xyMarginPad[2]);
682  pDcaXYvsZAll -> SetLeftMargin(xyMarginPad[3]);
683  pDcaXYvsZSel -> SetGrid(fGrid, fGrid);
684  pDcaXYvsZSel -> SetTicks(fTick, fTick);
685  pDcaXYvsZSel -> SetBorderMode(fMode);
686  pDcaXYvsZSel -> SetBorderSize(fBord);
687  pDcaXYvsZSel -> SetLogz(fLogZ);
688  pDcaXYvsZSel -> SetTopMargin(xyMarginPad[0]);
689  pDcaXYvsZSel -> SetRightMargin(xyMarginPad[1]);
690  pDcaXYvsZSel -> SetBottomMargin(xyMarginPad[2]);
691  pDcaXYvsZSel -> SetLeftMargin(xyMarginPad[3]);
692  cDcaXYvsZ -> cd();
693  pDcaXYvsZAll -> Draw();
694  pDcaXYvsZSel -> Draw();
695  pDcaXYvsZAll -> cd();
696  hDcaXYvsZAll -> Draw("colz");
697  pDcaXYvsZSel -> cd();
698  hDcaXYvsZSel -> Draw("colz");
699  ptInfo -> Draw();
700  fOutput -> cd();
701  cDcaXYvsZ -> Write();
702  cDcaXYvsZ -> Close();
703  cout << " Made plots." << endl;
704 
705  // save output & close ------------------------------------------------------
706 
707  // save histograms
708  fOutput -> cd();
709  for (size_t iDca = 0; iDca < NDca; iDca++) {
710  arrDcaWidth[iDca] -> Write();
711  arrWidthFits[iDca] -> Write();
712  arrDcaVsPtAll[iDca] -> Write();
713  arrDcaVsPtSel[iDca] -> Write();
714  arrFitsNeg[iDca] -> Write();
715  arrFitsPos[iDca] -> Write();
716  }
717  hDcaXYvsZAll -> Write();
718  hDcaXYvsZSel -> Write();
719 
720  // close flies
721  fOutput -> cd();
722  fOutput -> Close();
723  fInput -> cd();
724  fInput -> Close();
725  cout << " Finished with calculation!\n" << endl;
726 
727 }
728 
729 // end ------------------------------------------------------------------------