ATLAS Offline Software
Loading...
Searching...
No Matches
MonitoringFile_IDEnhancedPrimaryVertex.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2025 CERN for the benefit of the ATLAS collaboration
3 */
4
5// **********************************************************************
6// $Id: MonitoringFile_IDEnhancedPrimaryVertex.cxx,v 1.1 2009-04-03 08:48:00 ponyisi Exp $
7// **********************************************************************
8
9#include "DataQualityUtils/MonitoringFile.h"
10
11// #include <cmath>
12#include <vector>
13
14#include <TGraph.h>
15#include <TGraphErrors.h>
16#include <TGraphAsymmErrors.h>
17#include <TCanvas.h>
18#include <TF1.h>
19#include <TFile.h>
20#include <TH1.h>
21#include <TH2.h>
22#include <TKey.h>
23#include <TProfile.h>
24#include <TMath.h>
25#include <TString.h>
26
27namespace dqutils {
28// define helper methods here rather then in this central MonitoringFile.h "beast"
29 void plotResolution(const TString& coordinate, const TString& versus);
30 void plotEfficiency();
31 double error_func(float x, const Double_t* par);
32 double scaleFactorFitFcn(double* x, double* par);
33 std::vector<float> stableGaussianFit(TH1* histo);
34
35 void MonitoringFile::pv_PrimaryVertexMonitoring_calcResoAndEfficiency(const std::string& inFilename,
36 bool /* isIncremental */) {
37// std::cout << "\n";
38// std::cout << "Running Inner-Detector primary vertexing monitoring analysis\n";
39// std::cout << "\nWarning messages from fitting and histogram updating follow:\n\n";
40
41 TFile* f = TFile::Open(inFilename.c_str(), "UPDATE");
42
43 if (f == 0 || !f->IsOpen()) {
44 //std::cerr << "MonitoringFile::PrimaryVertexMonitoring_calcResoAndEfficiency(): "
45// << "Input file not opened\n";
46 delete f;
47 return;
48 }
49 if (f->GetSize() < 1000.) {
50// std::cerr << "MonitoringFile::PrimaryVertexMonitoring_calcResoAndEfficiency(): "
51// << "MonitoringFile empty\n";
52 delete f;
53 return;
54 }
55
56 bool dirExists = false;
57 TString run_dir;
58 TIter next_run(f->GetListOfKeys());
59 TKey* key_run(0);
60 while ((key_run = dynamic_cast<TKey*>(next_run())) != 0) {
61 TObject* obj_run = key_run->ReadObj();
62 TDirectory* tdir_run = dynamic_cast<TDirectory*>(obj_run);
63 if (tdir_run != 0) {
64 std::string tdir_run_name(tdir_run->GetName());
65 if (tdir_run_name.find("run") != std::string::npos) {
66 run_dir = tdir_run_name;
67
68 dirExists = f->GetDirectory(run_dir + "/InDetGlobal/PrimaryVertex");
69 if (dirExists) {
70// std::cout << run_dir << "/InDetGlobal/PrimaryVertex exists. Creating and saving plots:" << std::endl;
71 f->cd(run_dir + "/InDetGlobal/PrimaryVertex");
72 }
73 }
74 } else {
75 delete obj_run;
76 }
77 }
78
79 if (!dirExists) {
80 //std::cerr << "Either enhanced vertex info is not there or monitoring file was produced outside T0 context.
81 // Trying dir: InDetGlobal/PrimaryVertex\n";
82 dirExists = f->GetDirectory("InDetGlobal/PrimaryVertex");
83 if (dirExists) f->cd("InDetGlobal/PrimaryVertex");
84 }
85
86 if (dirExists) {
87 // create histos
88 plotResolution("X", "Ntrk");
89 plotResolution("Y", "Ntrk");
90 plotResolution("Z", "Ntrk");
91
92 plotResolution("X", "SumPt2");
93 plotResolution("Y", "SumPt2");
94 plotResolution("Z", "SumPt2");
95
96 plotEfficiency();
97 }
98
99 f->Close();
100 delete f;
101// std::cout << "\n";
102// std::cout << "Finish Inner-Detector primary vertex monitoring analysis"<<std::endl;
103 return;
104 }
105
106 void
107 plotResolution(const TString& coordinate = "Z", const TString& versus = "Ntrk") {
108// cout << "Creating and writing histos for resolution " << coordinate << " versus " << versus << endl;
109
110 TH2F* h_Vrt_pullVsSomething_split(0);
111 TH2F* h_Vrt_err_vs_Something(0);
112// TH2F* h_Vrt_Tag_err_vs_Something(0);
113 TString xAxisLabel("");
114
115 if (versus == "Ntrk") {
116 h_Vrt_pullVsSomething_split = (TH2F*) gDirectory->Get("Vrt_" + coordinate + "pullVsNtrkAverage_split");
117 h_Vrt_err_vs_Something = (TH2F*) gDirectory->Get("Vrt_" + coordinate + "err_vs_ntrk");
118// h_Vrt_Tag_err_vs_Something = (TH2F*)gDirectory->Get("Vrt_Tag_"+coordinate+"err_vs_ntrk");
119 xAxisLabel = "Number of fitted tracks";
120 } else if (versus == "SumPt2") {
121 h_Vrt_pullVsSomething_split = (TH2F*) gDirectory->Get("Vrt_" + coordinate + "pullVsPt2Average_split");
122 h_Vrt_err_vs_Something = (TH2F*) gDirectory->Get("Vrt_" + coordinate + "err_vs_pt2");
123// h_Vrt_Tag_err_vs_Something = (TH2F*)gDirectory->Get("Vrt_Tag_"+coordinate+"err_vs_pt2");
124 xAxisLabel = "#sqrt{#sum p_{T}^{2}} [GeV]";
125 } else return;
126
127 //if (h_Vrt_pullVsSomething_split == 0) std::cout << "h_Vrt_pullVsSomething_split has zero pointer!" << std::endl;
128 //if (h_Vrt_err_vs_Something == 0) std::cout << "h_Vrt_err_vs_Something has zero pointer!" << std::endl;
129 if (h_Vrt_pullVsSomething_split == 0 or h_Vrt_err_vs_Something == 0) return;
130
131 int n_bins = h_Vrt_pullVsSomething_split->GetNbinsX();
132 std::vector<float> rms_z;
133 std::vector<float> rms_z_er;
134 std::vector<float> sigma_z;
135 std::vector<float> sigma_z_er;
136 std::vector<float> bins_z_nt;
137 std::vector<float> bins_z_nt_er;
138
139// root starts counting the bins at 1, i.e. bin 1 holds NTrk = 0. or sqrt(sumpt2) = 0. - 0.25. GeV
140// std::cout << "n bins: " << n_bins << "\tTotal entries: " << h_Vrt_pullVsSomething_split->GetEntries() << std::endl;
141// TH1D * nTrksPerVertex = h_Vrt_pullVsSomething_split->ProjectionX("projectionNTrks_"+coordinate+"_"+versus);
142
143// TH1D * profileZFull = h_Vrt_pullVsSomething_split->ProjectionY("projectionPullsFull_"+coordinate+"_"+versus);
144
145 Int_t startBin = 0;
146 TH1D* profileZ = 0;
147 const Int_t minEntriesForKFactorBin = 1000;
148 for (int bin_count = 1; bin_count < n_bins + 1; bin_count++) {
149 //Fixed binning
150// TH1D *profileZ = h_Vrt_pullVsSomething_split->ProjectionY("projectionPulls", bin_count, bin_count,"e");
151// Double_t binCenter = h_Vrt_pullVsSomething_split->GetXaxis()->GetBinCenter(bin_count);
152// Double_t binWidth = h_Vrt_pullVsSomething_split->GetXaxis()->GetBinWidth(bin_count)/2.;
153
154 //Variable binning
155 TH1D* profileZTmp = h_Vrt_pullVsSomething_split->ProjectionY("projectionPulls", bin_count, bin_count, "e");
156 //cout << "Bin: " << bin_count << ", Entries: " << profileZTmp->GetEntries() << endl;
157 if (profileZ == 0) {
158 startBin = bin_count;
159 profileZ = (TH1D*) profileZTmp->Clone("projectionPulls_Integrated");
160 //cout << "StartBin = " << startBin << endl;
161 } else {
162 profileZ->Add(profileZTmp);
163 }
164 delete profileZTmp;
165 profileZTmp = 0;
166 if ((profileZ->GetEntries() < minEntriesForKFactorBin) && (bin_count < n_bins)) //not enough entries, not last bin
167 continue;
168
169 Double_t lowEdge = h_Vrt_pullVsSomething_split->GetXaxis()->GetBinLowEdge(startBin);
170 Double_t highEdge = h_Vrt_pullVsSomething_split->GetXaxis()->GetBinLowEdge(bin_count) +
171 h_Vrt_pullVsSomething_split->GetXaxis()->GetBinWidth(bin_count);
172 Double_t binCenter = (lowEdge + highEdge) / 2;
173 Double_t binWidth = (highEdge - lowEdge) / 2; //half of the bin width
174 //cout << "Bin done: " << binCenter << " +/- " << binWidth << ", Entries: " << profileZ->GetEntries() << endl;
175 // variable binning end
176
177 bins_z_nt.push_back(binCenter);
178 bins_z_nt_er.push_back(binWidth); // dummy error of binwidth for now
179
180 rms_z.push_back(profileZ->GetRMS());
181 rms_z_er.push_back(profileZ->GetRMSError());
182
183 //making a gaussian fit if there is anough entries
184 if (profileZ->GetEntries() > 100.) {
185 std::vector<float> fit_res = stableGaussianFit(profileZ);
186 sigma_z.push_back(fit_res[0]);
187 sigma_z_er.push_back(fit_res[1]);
188 } else {
189 sigma_z.push_back(0.);
190 sigma_z_er.push_back(0.);
191 }//end of good number of bins selection
192
193 delete profileZ; // must keep this to delete the projection from memory (next one has same name!)
194 profileZ = 0;
195 }//end of loop over all the ntrk bins
196
197 TGraphErrors* krms_z_vs_ntrk = new TGraphErrors(
198 bins_z_nt.size(), &(bins_z_nt[0]), &(rms_z[0]), &(bins_z_nt_er[0]), &(rms_z_er[0]));
199 krms_z_vs_ntrk->GetYaxis()->SetTitle(coordinate + " scale factor from RMS");
200 krms_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel);
201 krms_z_vs_ntrk->SetTitle("scaleFactor" + coordinate + "_RMS");
202 krms_z_vs_ntrk->SetName("scaleFactor" + coordinate + "_" + versus + "_RMS");
203
204 TGraphErrors* kgs_z_vs_ntrk = new TGraphErrors(
205 bins_z_nt.size(), &(bins_z_nt[0]), &(sigma_z[0]), &(bins_z_nt_er[0]), &(sigma_z_er[0]));
206 kgs_z_vs_ntrk->GetYaxis()->SetTitle(coordinate + " scale factor from gauss fit");
207 kgs_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel);
208 kgs_z_vs_ntrk->SetTitle("scaleFactor" + coordinate + "_Fit");
209 kgs_z_vs_ntrk->SetName("scaleFactor_" + coordinate + "_" + versus + "_Fit");
210
211// approximating the graph with 2nd order polynomial.
212 float maxFitRange(100.);
213 float minFitRange(2.);
214 if (versus == "SumPt2") {
215 minFitRange = 0.5;
216 maxFitRange = 20.;
217 }
218 TF1* kgs_z_ntrk_fit;
219 const Double_t* kgs_z_ntrk_fit_er;
220 int fitResKFactorMethod = 2; // set by hand for now
221 if (fitResKFactorMethod == 1) {
222 //Fit with a pol2
223 //coverity[DEADCODE]
224 kgs_z_vs_ntrk->Fit("pol2", "Q", "", minFitRange, maxFitRange);
225 kgs_z_vs_ntrk->GetFunction("pol2")->SetLineColor(kRed);
226 kgs_z_ntrk_fit = kgs_z_vs_ntrk->GetFunction("pol2");
227 kgs_z_ntrk_fit_er = kgs_z_ntrk_fit->GetParErrors();
228 } else if (fitResKFactorMethod == 2) {
229 //Fit with a pol1
230 kgs_z_vs_ntrk->Fit("pol1", "Q", "", minFitRange, maxFitRange);
231 kgs_z_vs_ntrk->GetFunction("pol1")->SetLineColor(kRed);
232 kgs_z_ntrk_fit = kgs_z_vs_ntrk->GetFunction("pol1");
233 kgs_z_ntrk_fit_er = kgs_z_ntrk_fit->GetParErrors();
234 //coverity[DEADCODE]
235 } else if (fitResKFactorMethod == 3) {
236 TF1* kgsFitFcn = new TF1("kgsFitFcn", scaleFactorFitFcn, minFitRange, maxFitRange, 3);
237 kgsFitFcn->SetParameter(0, minFitRange);
238 kgsFitFcn->SetParameter(1, 1.0);
239 kgsFitFcn->SetParameter(2, 1.0);
240 for (int ifit = 0; ifit < 1; ifit++) //initial estimation of best parameters
241 kgs_z_vs_ntrk->Fit(kgsFitFcn, "Q");
242 kgs_z_vs_ntrk->Fit(kgsFitFcn, "Q"); //perform actual fit
243 kgs_z_ntrk_fit = kgsFitFcn;
244 kgs_z_ntrk_fit_er = kgsFitFcn->GetParErrors();
245/* cout << "ScaleFactor fit for " << coordinate << " vs " << versus << " (method " << fitResKFactorMethod << ")= "
246 << kgsFitFcn->GetParameter(0) << " +/- " << kgsFitFcn->GetParError(0) << " "
247 << kgsFitFcn->GetParameter(1) << " +/- " << kgsFitFcn->GetParError(1) << " "
248 << kgsFitFcn->GetParameter(2) << " +/- " << kgsFitFcn->GetParError(2) << endl; */
249 } else if (fitResKFactorMethod == 4) {
250 //constant fit
251 kgs_z_vs_ntrk->Fit("pol0", "Q", "", minFitRange, maxFitRange);
252 kgs_z_vs_ntrk->GetFunction("pol0")->SetLineColor(kRed);
253 kgs_z_ntrk_fit = kgs_z_vs_ntrk->GetFunction("pol0");
254 kgs_z_ntrk_fit_er = kgs_z_ntrk_fit->GetParErrors();
255/* cout << "ScaleFactor fit for " << coordinate << " vs " << versus << " (method " << fitResKFactorMethod << ")= "
256 << kgs_z_ntrk_fit->GetParameter(0) << " +/- " << kgs_z_ntrk_fit->GetParError(0) << endl; */
257 }
258
259// plotting the fit error of the unconstrained primary vertex and correcting them
260 int nbins_z_err_ntrk = h_Vrt_err_vs_Something->GetNbinsX();
261
262 std::vector<float> av_err_z;
263 std::vector<float> av_err_z_er;
264// std::vector<float> av_err_tag_z;
265// std::vector<float> av_err_tag_z_er;
266 std::vector<float> err_bins_z_nt;
267 std::vector<float> err_bins_z_nt_er;
268 std::vector<float> res_z;
269 std::vector<float> res_z_er;
270// std::vector<float> res_tag_z;
271// std::vector<float> res_tag_z_er;
272
273 for (int bin_count = 1; bin_count <= nbins_z_err_ntrk; ++bin_count) {
274 err_bins_z_nt.push_back(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count));
275 err_bins_z_nt_er.push_back(h_Vrt_err_vs_Something->GetXaxis()->GetBinWidth(bin_count) / 2.);
276
277 TH1D* profileY = h_Vrt_err_vs_Something->ProjectionY("projectionErrors", bin_count, bin_count, "e");
278// TH1D * profileYTag(0);
279// if (h_Vrt_Tag_err_vs_Something)
280// profileYTag = h_Vrt_Tag_err_vs_Something->ProjectionY("projectionErrorsTag",bin_count,bin_count,"e");
281
282 float mean = profileY->GetMean();
283 float mean_error = profileY->GetMeanError();
284// float meanTag(0);
285// float mean_errorTag(0);
286// if (profileYTag) {
287// meanTag = profileYTag->GetMean();
288// mean_errorTag = profileYTag->GetMeanError();
289// }
290 delete profileY;
291// delete profileYTag;
292
293 av_err_z.push_back(mean);
294 av_err_z_er.push_back(mean_error);
295// av_err_tag_z.push_back(meanTag);
296// av_err_tag_z_er.push_back(mean_errorTag);
297
298 //estimating the approximate k-factor and the error value
299 double pr_er = 0.0;
300 float val(0.);
301 if (fitResKFactorMethod == 1) {
302 //coverity[DEADCODE]
303 pr_er = error_func(bin_count, kgs_z_ntrk_fit_er);
304 } else if (fitResKFactorMethod == 2) {
305 val = h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count);
306 pr_er = TMath::Power(kgs_z_ntrk_fit_er[1] * val, 2) + TMath::Power(kgs_z_ntrk_fit_er[0], 2);
307 pr_er = TMath::Sqrt(pr_er);
308// cout << "val = " << val << ", pr_er = " << pr_er << ", p0er = " << kgs_z_ntrk_fit_er[0] << ", p1er = "<<
309// kgs_z_ntrk_fit_er[1] << endl;
310 //coverity[DEADCODE]
311 } else if (fitResKFactorMethod == 3) {
312 val = h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count);
313 //approximately the error on the plateau
314 pr_er = kgs_z_ntrk_fit_er[2];
315 } else if (fitResKFactorMethod == 4) {
316 pr_er = kgs_z_ntrk_fit_er[0];
317 }
318
319 res_z.push_back(mean * kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count)));
320 res_z_er.push_back(TMath::Sqrt(TMath::Power(mean_error *
321 kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(
322 bin_count)),
323 2) + TMath::Power(pr_er * mean, 2)));
324// res_tag_z.push_back(meanTag * kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count)));
325// res_tag_z_er.push_back(TMath::Sqrt(TMath::Power(mean_errorTag *
326// kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count)),2) + TMath::Power( pr_er * mean
327// ,2)));
328 }
329 TGraphErrors* res_z_vs_ntrk =
330 new TGraphErrors(err_bins_z_nt.size(), &(err_bins_z_nt[0]), &(res_z[0]), &(err_bins_z_nt_er[0]), &(res_z_er[0]));
331 res_z_vs_ntrk->GetYaxis()->SetTitle(coordinate + " Vertex Resolution [mm]");
332 res_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel);
333 res_z_vs_ntrk->SetTitle(coordinate + " Vertex Resolution");
334 res_z_vs_ntrk->SetName("resolution_" + coordinate + "_" + versus);
335
336// TGraphErrors * res_tag_z_vs_ntrk = new TGraphErrors(err_bins_z_nt.size(),
337// &(err_bins_z_nt[0]),&(res_tag_z[0]),&(err_bins_z_nt_er[0]), &(res_tag_z_er[0]) );
338// res_tag_z_vs_ntrk->GetYaxis()->SetTitle(coordinate+" Tagged Vertex Resolution [mm]");
339// res_tag_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel);
340// res_tag_z_vs_ntrk->SetTitle(coordinate+" Tagged Vertex Resolution");
341// res_tag_z_vs_ntrk->SetName("resolution_tag_"+coordinate+"_"+versus);
342
343// TGraphErrors * mean_err_z_vs_ntrk = new TGraphErrors(err_bins_z_nt.size(),
344// &(err_bins_z_nt[0]),&(av_err_z[0]),&(err_bins_z_nt_er[0]), &(av_err_z_er[0]) );
345// mean_err_z_vs_ntrk->GetYaxis()->SetTitle(coordinate+" mean vertex error [mm]");
346// mean_err_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel);
347// mean_err_z_vs_ntrk->SetTitle(coordinate+" Mean Vertex Error");
348// mean_err_z_vs_ntrk->SetName("resolution_"+coordinate+"_"+versus+"_Uncorrected"); //not corrected with k-factor
349
351 // Writing out
353 if (versus == "Ntrk") res_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
354 else res_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
355 res_z_vs_ntrk->GetYaxis()->SetRangeUser(0.0025, 1.);
356 res_z_vs_ntrk->Write("", TObject::kOverwrite);
357 delete res_z_vs_ntrk;
358
359 if (versus == "Ntrk") krms_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
360 else krms_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
361 krms_z_vs_ntrk->GetYaxis()->SetRangeUser(0.5, 1.3);
362 krms_z_vs_ntrk->Write("", TObject::kOverwrite);
363 delete krms_z_vs_ntrk;
364
365 if (versus == "Ntrk") kgs_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
366 else kgs_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
367 kgs_z_vs_ntrk->GetYaxis()->SetRangeUser(0.5, 1.3);
368 kgs_z_vs_ntrk->Write("", TObject::kOverwrite);
369 delete kgs_z_vs_ntrk;
370
371// nTrksPerVertex->GetYaxis()->SetTitle("Entries");
372// nTrksPerVertex->Draw();
373// nTrksPerVertex->Write("", TObject::kOverwrite);
374
375 return;
376 }
377
378// Reconstruction efficiency
379 void plotEfficiency() {
380// cout << "Creating and writing histos for efficiency" << endl;
381
382 // get histos but return if any histo is not there in input
383 TH1F* h_Vrt_split_tag_ntrk = (TH1F*) gDirectory->Get("Vrt_split_tag_ntrk");
384
385 if (h_Vrt_split_tag_ntrk == 0) return;
386
387 TH1F* h_Vrt_split_probe_ntrk = (TH1F*) gDirectory->Get("Vrt_split_probe_ntrk");
388 if (h_Vrt_split_probe_ntrk == 0) return;
389
390 TH1F* h_Vrt_split_matched_tag_ntrk = (TH1F*) gDirectory->Get("Vrt_split_matched_tag_ntrk");
391 if (h_Vrt_split_matched_tag_ntrk == 0) return;
392
393 TH1F* h_Vrt_split_matched_probe_ntrk = (TH1F*) gDirectory->Get("Vrt_split_matched_probe_ntrk");
394 if (h_Vrt_split_matched_probe_ntrk == 0) return;
395
396 TH1F* h_Vrt_split_dist_tag = (TH1F*) gDirectory->Get("Vrt_split_dist_tag");
397 if (h_Vrt_split_dist_tag == 0) return;
398
399 TH1F* h_Vrt_split_dist_probe = (TH1F*) gDirectory->Get("Vrt_split_dist_probe");
400 if (h_Vrt_split_dist_probe == 0) return;
401
402 // Use BayesDivide routing of TGraphAsymmErrors
403 TGraphAsymmErrors* g_Vrt_rec_eff_m1_split_vs_ntrk = new TGraphAsymmErrors();
404
405 g_Vrt_rec_eff_m1_split_vs_ntrk->BayesDivide(h_Vrt_split_probe_ntrk, h_Vrt_split_tag_ntrk);
406 g_Vrt_rec_eff_m1_split_vs_ntrk->SetName("g_RecEff_M1");
407
408 TGraphAsymmErrors* g_Vrt_sel_eff_m1_split_vs_ntrk = new TGraphAsymmErrors();
409 g_Vrt_sel_eff_m1_split_vs_ntrk->BayesDivide(h_Vrt_split_matched_probe_ntrk, h_Vrt_split_matched_tag_ntrk);
410 g_Vrt_sel_eff_m1_split_vs_ntrk->SetName("g_SelEff_M1");
411
412 // formatting and writing out
413 g_Vrt_rec_eff_m1_split_vs_ntrk->GetHistogram()->GetXaxis()->SetTitle("Number of tracks");
414 g_Vrt_rec_eff_m1_split_vs_ntrk->GetHistogram()->GetYaxis()->SetTitle("Reconstruction efficiency");
415 g_Vrt_rec_eff_m1_split_vs_ntrk->SetMarkerStyle(20);
416 g_Vrt_rec_eff_m1_split_vs_ntrk->Write("", TObject::kOverwrite);
417 delete g_Vrt_rec_eff_m1_split_vs_ntrk;
418
419 g_Vrt_sel_eff_m1_split_vs_ntrk->GetHistogram()->GetXaxis()->SetTitle("Number of tracks");
420 g_Vrt_sel_eff_m1_split_vs_ntrk->GetHistogram()->GetYaxis()->SetTitle("Selection Efficiency");
421 g_Vrt_sel_eff_m1_split_vs_ntrk->SetMarkerStyle(20);
422 g_Vrt_sel_eff_m1_split_vs_ntrk->Write("", TObject::kOverwrite);
423 delete g_Vrt_sel_eff_m1_split_vs_ntrk;
424
425 return;
426 }
427
428 std::vector<float> stableGaussianFit(TH1* histo) {
429 std::vector<float> results;
430 if (histo) {
431 double sigmas = 2.;
432
433 histo->Fit("gaus", "Q0", "", -2, 2);
434 TF1* func = histo->GetFunction("gaus");
435 double actualSigma = func->GetParameter(2);
436 double actualSigmaErr = func->GetParError(2);
437
438 for (int u = 0; u < 10; u++) {
439 histo->Fit("gaus", "Q0", "", -sigmas * actualSigma, sigmas * actualSigma);
440 func = histo->GetFunction("gaus");
441 actualSigma = func->GetParameter(2);
442 actualSigmaErr = func->GetParError(2);
443 }//end of the fitting loop
444
445 results.push_back(actualSigma);
446 results.push_back(actualSigmaErr);
447 } else {
448 results.push_back(0.);
449 results.push_back(0.);
450 }//end of protection against an empty histogram
451
452 return results;
453 }//end of stableGaussian Fit function
454
455 double error_func(float x, const Double_t* par) {
456//calculating the square of the propagated error on the fit values
457 return(TMath::Power(par[0], 2) + x * TMath::Power(par[1], 2) + TMath::Power(x * par[2], 2));
458 }
459
460 double scaleFactorFitFcn(double* x, double* par) {
461 // Truncated gaus-smeared step function
462 // par 0: mean of turn-on (and truncation point)
463 // par 1: slope of turn-on
464 // par 2: plateau
465 return (x[0] >= par[0]) * TMath::Erf(par[1] * x[0] - par[0]) * par[2];
466 }
467}
x
#define x
dqutils::MonitoringFile::pv_PrimaryVertexMonitoring_calcResoAndEfficiency
static void pv_PrimaryVertexMonitoring_calcResoAndEfficiency(const std::string &inFilename, bool isIncremental=false)
Definition MonitoringFile_IDEnhancedPrimaryVertex.cxx:35
mean
void mean(std::vector< double > &bins, std::vector< double > &values, const std::vector< std::string > &files, const std::string &histname, const std::string &tplotname, const std::string &label="")
Definition dependence.cxx:254
binWidth
void binWidth(TH1 *h)
Definition listroot.cxx:80
dqutils::stableGaussianFit
std::vector< float > stableGaussianFit(TH1 *histo)
Definition MonitoringFile_IDEnhancedPrimaryVertex.cxx:428
dqutils::error_func
double error_func(float x, const Double_t *par)
Definition MonitoringFile_IDEnhancedPrimaryVertex.cxx:455
dqutils::plotResolution
void plotResolution(const TString &coordinate, const TString &versus)
Definition MonitoringFile_IDEnhancedPrimaryVertex.cxx:107
dqutils::scaleFactorFitFcn
double scaleFactorFitFcn(double *x, double *par)
Definition MonitoringFile_IDEnhancedPrimaryVertex.cxx:460
Generated on for ATLAS Offline Software by doxygen 1.14.0