Analysis Software
Documentation for sPHENIX simulation software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
SDeltaPtCutStudy.ana.h
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file SDeltaPtCutStudy.ana.h
1 // ----------------------------------------------------------------------------
2 // 'SDeltaPtCutStudy.ana.h'
3 // Derek Anderson
4 // 07.07.2023
5 //
6 // Reads in the 'ntp_track' Ntuple
7 // generated by the SVtxEvaluator
8 // class and studies how deltapt/pt
9 // varies with quality cuts.
10 // ----------------------------------------------------------------------------
11 
12 #pragma once
13 
14 using namespace std;
15 
16 
17 
18 // analysis methods -----------------------------------------------------------
19 
20 void SDeltaPtCutStudy::ApplyFlatDeltaPtCuts() {
21 
22  // announce start of track loop
23  cout << " First loop over reco. tracks:" << endl;
24 
25  // 1st track loop
26  uint64_t nBytesTrk = 0;
27  for (uint64_t iTrk = 0; iTrk < nTrks; iTrk++) {
28 
29  // grab entry
30  const uint64_t bytesTrk = ntTrack -> GetEntry(iTrk);
31  if (bytesTrk < 0.) {
32  cerr << "WARNING: something wrong with track #" << iTrk << "! Aborting loop!" << endl;
33  break;
34  }
35  nBytesTrk += bytesTrk;
36 
37  // announce progress
38  const uint64_t iProgTrk = iTrk + 1;
39  if (iProgTrk == nTrks) {
40  cout << " Processing track " << iProgTrk << "/" << nTrks << "..." << endl;
41  } else {
42  cout << " Processing track " << iProgTrk << "/" << nTrks << "...\r" << flush;
43  }
44 
45  // do calculations
46  const double ptFrac = trk_pt / trk_gpt;
47  const double ptDelta = trk_deltapt / trk_pt;
48 
49  // apply trk cuts
50  const bool isInZVtxCut = (abs(trk_vz) < vzTrkMax);
51  const bool isInInttCut = (trk_nintt >= nInttTrkMin);
52  const bool isInMVtxCut = (trk_nlmaps > nMVtxTrkMin);
53  const bool isInTpcCut = (trk_ntpc > nTpcTrkMin);
54  const bool isInPtCut = (trk_pt > ptTrkMin);
55  const bool isInQualCut = (trk_quality < qualTrkMax);
56  const bool isGoodTrk = (isInZVtxCut && isInInttCut && isInMVtxCut && isInTpcCut && isInPtCut && isInQualCut);
57  if (!isGoodTrk) continue;
58 
59  // fill histograms
60  hPtDelta -> Fill(ptDelta);
61  hPtTrack -> Fill(trk_pt);
62  hPtFrac -> Fill(ptFrac);
63  hPtTrkTru -> Fill(trk_gpt);
64  hPtDeltaVsFrac -> Fill(ptFrac, ptDelta);
65  hPtDeltaVsTrue -> Fill(trk_gpt, ptDelta);
66  hPtDeltaVsTrack -> Fill(trk_pt, ptDelta);
67  hPtTrueVsTrack -> Fill(trk_pt, trk_gpt);
68 
69  // apply delta-pt cuts
70  const bool isNormalTrk = ((ptFrac > normRange[0]) && (ptFrac < normRange[1]));
71  for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
72  const bool isInDeltaPtCut = (ptDelta < ptDeltaMax[iCut]);
73  if (isInDeltaPtCut) {
74 
75  // fill histograms
76  hPtDeltaCut[iCut] -> Fill(ptDelta);
77  hPtTrackCut[iCut] -> Fill(trk_pt);
78  hPtFracCut[iCut] -> Fill(ptFrac);
79  hPtTrkTruCut[iCut] -> Fill(trk_gpt);
80  hPtDeltaVsFracCut[iCut] -> Fill(ptFrac, ptDelta);
81  hPtDeltaVsTrueCut[iCut] -> Fill(trk_gpt, ptDelta);
82  hPtDeltaVsTrackCut[iCut] -> Fill(trk_pt, ptDelta);
83  hPtTrueVsTrackCut[iCut] -> Fill(trk_pt, trk_gpt);
84 
85  // increment counters
86  if (isNormalTrk) {
87  ++nNormCut[iCut];
88  } else {
89  ++nWeirdCut[iCut];
90  }
91  }
92  } // end delta-pt cut
93  } // end 1st track loop
94 
95  cout << " First loop over reco. tracks finished!" << endl;
96  return;
97 
98 } // end 'ApplyFlatDeltaPtCuts()'
99 
100 
101 
102 void SDeltaPtCutStudy::ApplyPtDependentDeltaPtCuts() {
103 
104  // announce start of track loop
105  cout << " Second loop over reco. tracks:" << endl;
106 
107  // 2nd track loop
108  uint64_t nBytesTrk = 0;
109  for (uint64_t iTrk = 0; iTrk < nTrks; iTrk++) {
110 
111  // grab entry
112  const uint64_t bytesTrk = ntTrack -> GetEntry(iTrk);
113  if (bytesTrk < 0.) {
114  cerr << "WARNING: something wrong with track #" << iTrk << "! Aborting loop!" << endl;
115  break;
116  }
117  nBytesTrk += bytesTrk;
118 
119  // announce progress
120  const uint64_t iProgTrk = iTrk + 1;
121  if (iProgTrk == nTrks) {
122  cout << " Processing track " << iProgTrk << "/" << nTrks << "..." << endl;
123  } else {
124  cout << " Processing track " << iProgTrk << "/" << nTrks << "...\r" << flush;
125  }
126 
127  // do calculations
128  const double ptFrac = trk_pt / trk_gpt;
129  const double ptDelta = trk_deltapt / trk_pt;
130 
131  // apply trk cuts
132  const bool isInZVtxCut = (abs(trk_vz) < vzTrkMax);
133  const bool isInInttCut = (trk_nintt >= nInttTrkMin);
134  const bool isInMVtxCut = (trk_nlmaps > nMVtxTrkMin);
135  const bool isInTpcCut = (trk_ntpc > nTpcTrkMin);
136  const bool isInPtCut = (trk_pt > ptTrkMin);
137  const bool isInQualCut = (trk_quality < qualTrkMax);
138  const bool isGoodTrk = (isInZVtxCut && isInInttCut && isInMVtxCut && isInTpcCut && isInPtCut && isInQualCut);
139  if (!isGoodTrk) continue;
140 
141  // apply delta-pt cuts
142  const bool isNormalTrk = ((ptFrac > normRange[0]) && (ptFrac < normRange[1]));
143  for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
144 
145  // get bounds
146  const float ptDeltaMin = fMuLoProj[iSig] -> Eval(trk_pt);
147  const float ptDeltaMax = fMuHiProj[iSig] -> Eval(trk_pt);
148 
149  const bool isInDeltaPtSigma = ((ptDelta >= ptDeltaMin) && (ptDelta <= ptDeltaMax));
150  if (isInDeltaPtSigma) {
151 
152  // fill histograms
153  hPtDeltaSig[iSig] -> Fill(ptDelta);
154  hPtTrackSig[iSig] -> Fill(trk_pt);
155  hPtFracSig[iSig] -> Fill(ptFrac);
156  hPtTrkTruSig[iSig] -> Fill(trk_gpt);
157  hPtDeltaVsFracSig[iSig] -> Fill(ptFrac, ptDelta);
158  hPtDeltaVsTrueSig[iSig] -> Fill(trk_gpt, ptDelta);
159  hPtDeltaVsTrackSig[iSig] -> Fill(trk_pt, ptDelta);
160  hPtTrueVsTrackSig[iSig] -> Fill(trk_pt, trk_gpt);
161 
162  // increment counters
163  if (isNormalTrk) {
164  ++nNormSig[iSig];
165  } else {
166  ++nWeirdSig[iSig];
167  }
168  }
169  } // end delta-pt cut
170  } // end 1st track loop
171 
172  cout << " Second loop over reco. tracks finished!" << endl;
173  return;
174 
175 } // end 'ApplyPtDependentDeltaptCuts()'
176 
177 
178 
179 void SDeltaPtCutStudy::FillTruthHistograms() {
180 
181  // announce start of truth loop
182  cout << " Loop over particles:" << endl;
183 
184  // truth loop
185  uint64_t nBytesTru = 0;
186  for (uint64_t iTru = 0; iTru < nTrus; iTru++) {
187 
188  // grab entry
189  const uint64_t bytesTru = ntTruth -> GetEntry(iTru);
190  if (bytesTru < 0.) {
191  cerr << "WARNING: something wrong with particle #" << iTru << "! Aborting loop!" << endl;
192  break;
193  }
194  nBytesTru += bytesTru;
195 
196  // announce progress
197  const uint64_t iProgTru = iTru + 1;
198  if (iProgTru == nTrus) {
199  cout << " Processing particle " << iProgTru << "/" << nTrus << "..." << endl;
200  } else {
201  cout << " Processing particle " << iProgTru << "/" << nTrus << "...\r" << flush;
202  }
203 
204  // fill truth histogram
205  const bool isPrimary = (tru_gprimary == 1);
206  if (isPrimary) {
207  hPtTruth -> Fill(tru_gpt);
208  }
209  } // end track loop
210  cout << " Loop over particles finished!" << endl;
211 
212 } // end 'FillTruthHistograms()'
213 
214 
215 
216 void SDeltaPtCutStudy::CreateSigmaGraphs() {
217 
218  // for graph names
219  const TString sMuHiBase = "MeanPlusSigma";
220  const TString sMuLoBase = "MeanMinusSigma";
221  const TString sSigBase = "ProjectionSigma";
222  const TString sMuBase = "ProjectionMean";
223 
224  // projection fit names
225  vector<TString> sFitProj(nProj);
226  for (size_t iProj = 0; iProj < nProj; iProj++) {
227  sFitProj[iProj] = "f";
228  sFitProj[iProj].Append(sPtProjBase.Data());
229  sFitProj[iProj].Append(sProjSuffix[iProj].Data());
230  }
231 
232  // project slices of delta-pt and get sigmas
233  const uint32_t fWidFit = 2;
234  const uint32_t fLinFit = 1;
235  for (size_t iProj = 0; iProj < nProj; iProj++) {
236 
237  // do projection
238  const uint32_t iBinProj = hPtDeltaVsTrack -> GetXaxis() -> FindBin(ptProj[iProj]);
239  hPtDeltaProj[iProj] = hPtDeltaVsTrack -> ProjectionY(sPtProj[iProj], iBinProj, iBinProj, "");
240 
241  // get initial values for fit
242  const float ampGuess = hPtDeltaProj[iProj] -> GetMaximum();
243  const float muGuess = hPtDeltaProj[iProj] -> GetMean();
244  const float sigGuess = hPtDeltaProj[iProj] -> GetRMS();
245 
246  // fit with gaussian
247  fPtDeltaProj[iProj] = new TF1(sFitProj[iProj].Data(), "gaus", deltaFitRange[0], deltaFitRange[1]);
248  fPtDeltaProj[iProj] -> SetLineColor(fColFit[iProj]);
249  fPtDeltaProj[iProj] -> SetLineStyle(fLinFit);
250  fPtDeltaProj[iProj] -> SetLineWidth(fWidFit);
251  fPtDeltaProj[iProj] -> SetParameter(0, ampGuess);
252  fPtDeltaProj[iProj] -> SetParameter(1, muGuess);
253  fPtDeltaProj[iProj] -> SetParameter(2, sigGuess);
254  hPtDeltaProj[iProj] -> Fit(sFitProj[iProj].Data(), "R");
255 
256  // add values to arrays
257  muProj[iProj] = fPtDeltaProj[iProj] -> GetParameter(1);
258  sigProj[iProj] = fPtDeltaProj[iProj] -> GetParameter(2);
259  for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
260  muHiProj[iSig][iProj] = muProj[iProj] + (ptDeltaSig[iSig] * sigProj[iProj]);
261  muLoProj[iSig][iProj] = muProj[iProj] - (ptDeltaSig[iSig] * sigProj[iProj]);
262  }
263  } // end projection loop
264  cout << " Obtained delta-pt projections, fits, and sigmas." << endl;
265 
266  // sigma graph names
267  TString sMuProj("gr");
268  TString sSigProj("gr");
269  sMuProj.Append(sMuBase.Data());
270  sSigProj.Append(sSigBase.Data());
271 
272  vector<TString> sGrMuHiProj(nSigCuts);
273  vector<TString> sGrMuLoProj(nSigCuts);
274  vector<TString> sFnMuHiProj(nSigCuts);
275  vector<TString> sFnMuLoProj(nSigCuts);
276  for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
277  sGrMuHiProj[iSig] = "gr";
278  sGrMuLoProj[iSig] = "gr";
279  sFnMuHiProj[iSig] = "f";
280  sFnMuLoProj[iSig] = "f";
281  sGrMuHiProj[iSig].Append(sMuHiBase.Data());
282  sGrMuLoProj[iSig].Append(sMuLoBase.Data());
283  sFnMuHiProj[iSig].Append(sMuHiBase.Data());
284  sFnMuLoProj[iSig].Append(sMuLoBase.Data());
285  sGrMuHiProj[iSig].Append(sSigSuffix[iSig].Data());
286  sGrMuLoProj[iSig].Append(sSigSuffix[iSig].Data());
287  sFnMuHiProj[iSig].Append(sSigSuffix[iSig].Data());
288  sFnMuLoProj[iSig].Append(sSigSuffix[iSig].Data());
289  }
290 
291  // turn std::vectors into TVectors
292  TVectorD tvecPtProj(ptProj.size(), ptProj.data());
293  TVectorD tvecMuProj(muProj.size(), muProj.data());
294  TVectorD tvecSigProj(sigProj.size(), sigProj.data());
295 
296  vector<TVectorD> tvecMuHiProj;
297  vector<TVectorD> tvecMuLoProj;
298  for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
299  tvecMuHiProj.push_back(TVectorD(muHiProj[iSig].size(), muHiProj[iSig].data()));
300  tvecMuLoProj.push_back(TVectorD(muLoProj[iSig].size(), muLoProj[iSig].data()));
301  }
302 
303  // construct sigma graphs
304  grMuProj = new TGraph(tvecPtProj, tvecMuProj);
305  grSigProj = new TGraph(tvecPtProj, tvecSigProj);
306  grMuProj -> SetName(sMuProj);
307  grSigProj -> SetName(sSigProj);
308 
309  // fit sigma graphs
310  for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
311 
312  // create graphs
313  grMuHiProj[iSig] = new TGraph(tvecPtProj, tvecMuHiProj[iSig]);
314  grMuLoProj[iSig] = new TGraph(tvecPtProj, tvecMuLoProj[iSig]);
315  grMuHiProj[iSig] -> SetName(sGrMuHiProj[iSig].Data());
316  grMuLoProj[iSig] -> SetName(sGrMuLoProj[iSig].Data());
317 
318  // create fit functions
319  fMuHiProj[iSig] = new TF1(sFnMuHiProj[iSig].Data(), "pol2", rPtRange[0], rPtRange[1]);
320  fMuLoProj[iSig] = new TF1(sFnMuLoProj[iSig].Data(), "pol2", rPtRange[0], rPtRange[1]);
321  fMuHiProj[iSig] -> SetLineColor(fColSigFit[iSig]);
322  fMuLoProj[iSig] -> SetLineColor(fColSigFit[iSig]);
323  fMuHiProj[iSig] -> SetLineStyle(fLinFit);
324  fMuLoProj[iSig] -> SetLineStyle(fLinFit);
325  fMuHiProj[iSig] -> SetLineWidth(fWidFit);
326  fMuLoProj[iSig] -> SetLineWidth(fWidFit);
327  fMuHiProj[iSig] -> SetParameter(0, sigHiGuess[0]);
328  fMuLoProj[iSig] -> SetParameter(0, sigLoGuess[0]);
329  fMuHiProj[iSig] -> SetParameter(1, sigHiGuess[1]);
330  fMuLoProj[iSig] -> SetParameter(1, sigLoGuess[1]);
331  fMuHiProj[iSig] -> SetParameter(2, sigHiGuess[2]);
332  fMuLoProj[iSig] -> SetParameter(2, sigLoGuess[2]);
333 
334  // do fitting
335  grMuHiProj[iSig] -> Fit(sFnMuHiProj[iSig].Data(), "", "", ptFitRange[0], ptFitRange[1]);
336  grMuLoProj[iSig] -> Fit(sFnMuLoProj[iSig].Data(), "", "", ptFitRange[0], ptFitRange[1]);
337  }
338 
339  cout << " Created and fit sigma graphs." << endl;
340  return;
341 
342 } // end 'CreateSigmaGraphs()'
343 
344 
345 
346 void SDeltaPtCutStudy::CalculateRejectionFactors() {
347 
348  // for graph names
349  const TString sRejCutBase = "Reject_flatDPtCut";
350  const TString sRejSigBase = "Reject_sigmaCut";
351 
352  // calculate flat delta-pt rejection factors
353  for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
354  rejCut[iCut] = (double) nNormCut[iCut] / (double) nWeirdCut[iCut];
355  }
356  cout << " Calculated flat delta-pt rejection factors." << endl;
357 
358  // calculate pt-dependent delta-pt rejection factors
359  for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
360  rejSig[iSig] = (double) nNormSig[iSig] / (double) nWeirdSig[iSig];
361  }
362  cout << " Calculated pt-depdendent delta-pt rejection factors\n"
363  << " Rejection factors:\n"
364  << " Flat delta-pt cuts"
365  << endl;
366 
367  // announce flat delta-pt rejection factors
368  for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
369  cout << " n(Norm, Weird) = (" << nNormCut[iCut] << ", " << nWeirdCut[iCut] << "), rejection = " << rejCut[iCut] << endl;
370  }
371 
372  // announce pt-dependent delta-pt rejection factors
373  cout << " Pt-dependent delta-pt cuts" << endl;
374  for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
375  cout << " n(Norm, Weird) = (" << nNormSig[iSig] << ", " << nWeirdSig[iSig] << "), rejection = " << rejSig[iSig] << endl;
376  }
377 
378  // graph names
379  TString sRejCut("gr");
380  TString sRejSig("gr");
381  sRejCut.Append(sRejCutBase.Data());
382  sRejSig.Append(sRejSigBase.Data());
383 
384  // convert vectors to TVectors
385  TVectorD tvecPtDeltaMax(ptDeltaMax.size(), ptDeltaMax.data());
386  TVectorD tvecPtDeltaSig(ptDeltaSig.size(), ptDeltaSig.data());
387  TVectorD tvecRejCut(rejCut.size(), rejCut.data());
388  TVectorD tvecRejSig(rejSig.size(), rejSig.data());
389 
390  // make rejection graphs
391  grRejCut = new TGraph(tvecPtDeltaMax, tvecRejCut);
392  grRejSig = new TGraph(tvecPtDeltaSig, tvecRejSig);
393  grRejCut -> SetName(sRejCut.Data());
394  grRejSig -> SetName(sRejSig.Data());
395 
396  cout << " Made rejection factor graph." << endl;
397  return;
398 
399 } // end 'CalculateRejectionFactors()'
400 
401 
402 
403 void SDeltaPtCutStudy::CalculateEfficiencies() {
404 
405  // for histogram names
406  const TString sEffBase = "Efficiency";
407 
408  // rebin histograms if needed
409  if (doEffRebin) {
410  hPtTruth -> Rebin(nEffRebin);
411  hPtTrkTru -> Rebin(nEffRebin);
412  for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
413  hPtTrkTruCut[iCut] -> Rebin(nEffRebin);
414  }
415  for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
416  hPtTrkTruSig[iSig] -> Rebin(nEffRebin);
417  }
418  cout << " Rebinned efficiency histograms." << endl;
419  }
420 
421  TString sEff("h");
422  sEff.Append(sEffBase.Data());
423 
424  // create flat delta-pt cut efficiency names
425  vector<TString> sEffCut(nDPtCuts);
426  for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
427  sEffCut[iCut] = "h";
428  sEffCut[iCut].Append(sEffBase.Data());
429  sEffCut[iCut].Append(sDPtSuffix[iCut].Data());
430  }
431 
432  // create pt-dependent delta-pt cut efficiency names
433  vector<TString> sEffSig(nSigCuts);
434  for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
435  sEffSig[iSig] = "h";
436  sEffSig[iSig].Append(sEffBase.Data());
437  sEffSig[iSig].Append(sSigSuffix[iSig].Data());
438  }
439 
440  hEff = (TH1D*) hPtTruth -> Clone();
441  hEff -> SetName(sEff.Data());
442  hEff -> Reset("ICES");
443  hEff -> Divide(hPtTrkTru, hPtTruth, 1., 1.);
444 
445  // calculate flat delta-pt cut efficiencies
446  for (size_t iCut = 0; iCut < nDPtCuts; iCut++) {
447  hEffCut[iCut] = (TH1D*) hPtTruth -> Clone();
448  hEffCut[iCut] -> SetName(sEffCut[iCut].Data());
449  hEffCut[iCut] -> Reset("ICES");
450  hEffCut[iCut] -> Divide(hPtTrkTruCut[iCut], hPtTruth, 1., 1.);
451  }
452 
453  // calculate pt-dependent delta-pt cut efficiencies
454  for (size_t iSig = 0; iSig < nSigCuts; iSig++) {
455  hEffSig[iSig] = (TH1D*) hPtTruth -> Clone();
456  hEffSig[iSig] -> SetName(sEffSig[iSig].Data());
457  hEffSig[iSig] -> Reset("ICES");
458  hEffSig[iSig] -> Divide(hPtTrkTruSig[iSig], hPtTruth, 1., 1.);
459  }
460 
461  cout << " Calculated efficiencies." << endl;
462  return;
463 
464 } // end 'CalculateEfficiencies()'
465 
466 // end ------------------------------------------------------------------------