ATLAS Offline Software
Loading...
Searching...
No Matches
MonitoringFile_IDEnhancedPrimaryVertex.cxx
Go to the documentation of this file.
1/*
2 Copyright (C) 2002-2017 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 kgs_z_vs_ntrk->Fit("pol2", "Q", "", minFitRange, maxFitRange);
224 kgs_z_vs_ntrk->GetFunction("pol2")->SetLineColor(kRed);
225 kgs_z_ntrk_fit = kgs_z_vs_ntrk->GetFunction("pol2");
226 kgs_z_ntrk_fit_er = kgs_z_ntrk_fit->GetParErrors();
227 } else if (fitResKFactorMethod == 2) {
228 //Fit with a pol1
229 kgs_z_vs_ntrk->Fit("pol1", "Q", "", minFitRange, maxFitRange);
230 kgs_z_vs_ntrk->GetFunction("pol1")->SetLineColor(kRed);
231 kgs_z_ntrk_fit = kgs_z_vs_ntrk->GetFunction("pol1");
232 kgs_z_ntrk_fit_er = kgs_z_ntrk_fit->GetParErrors();
233 } else if (fitResKFactorMethod == 3) {
234 TF1* kgsFitFcn = new TF1("kgsFitFcn", scaleFactorFitFcn, minFitRange, maxFitRange, 3);
235 kgsFitFcn->SetParameter(0, minFitRange);
236 kgsFitFcn->SetParameter(1, 1.0);
237 kgsFitFcn->SetParameter(2, 1.0);
238 for (int ifit = 0; ifit < 1; ifit++) //initial estimation of best parameters
239 kgs_z_vs_ntrk->Fit(kgsFitFcn, "Q");
240 kgs_z_vs_ntrk->Fit(kgsFitFcn, "Q"); //perform actual fit
241 kgs_z_ntrk_fit = kgsFitFcn;
242 kgs_z_ntrk_fit_er = kgsFitFcn->GetParErrors();
243/* cout << "ScaleFactor fit for " << coordinate << " vs " << versus << " (method " << fitResKFactorMethod << ")= "
244 << kgsFitFcn->GetParameter(0) << " +/- " << kgsFitFcn->GetParError(0) << " "
245 << kgsFitFcn->GetParameter(1) << " +/- " << kgsFitFcn->GetParError(1) << " "
246 << kgsFitFcn->GetParameter(2) << " +/- " << kgsFitFcn->GetParError(2) << endl; */
247 } else if (fitResKFactorMethod == 4) {
248 //constant fit
249 kgs_z_vs_ntrk->Fit("pol0", "Q", "", minFitRange, maxFitRange);
250 kgs_z_vs_ntrk->GetFunction("pol0")->SetLineColor(kRed);
251 kgs_z_ntrk_fit = kgs_z_vs_ntrk->GetFunction("pol0");
252 kgs_z_ntrk_fit_er = kgs_z_ntrk_fit->GetParErrors();
253/* cout << "ScaleFactor fit for " << coordinate << " vs " << versus << " (method " << fitResKFactorMethod << ")= "
254 << kgs_z_ntrk_fit->GetParameter(0) << " +/- " << kgs_z_ntrk_fit->GetParError(0) << endl; */
255 }
256
257// plotting the fit error of the unconstrained primary vertex and correcting them
258 int nbins_z_err_ntrk = h_Vrt_err_vs_Something->GetNbinsX();
259
260 std::vector<float> av_err_z;
261 std::vector<float> av_err_z_er;
262// std::vector<float> av_err_tag_z;
263// std::vector<float> av_err_tag_z_er;
264 std::vector<float> err_bins_z_nt;
265 std::vector<float> err_bins_z_nt_er;
266 std::vector<float> res_z;
267 std::vector<float> res_z_er;
268// std::vector<float> res_tag_z;
269// std::vector<float> res_tag_z_er;
270
271 for (int bin_count = 1; bin_count <= nbins_z_err_ntrk; ++bin_count) {
272 err_bins_z_nt.push_back(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count));
273 err_bins_z_nt_er.push_back(h_Vrt_err_vs_Something->GetXaxis()->GetBinWidth(bin_count) / 2.);
274
275 TH1D* profileY = h_Vrt_err_vs_Something->ProjectionY("projectionErrors", bin_count, bin_count, "e");
276// TH1D * profileYTag(0);
277// if (h_Vrt_Tag_err_vs_Something)
278// profileYTag = h_Vrt_Tag_err_vs_Something->ProjectionY("projectionErrorsTag",bin_count,bin_count,"e");
279
280 float mean = profileY->GetMean();
281 float mean_error = profileY->GetMeanError();
282// float meanTag(0);
283// float mean_errorTag(0);
284// if (profileYTag) {
285// meanTag = profileYTag->GetMean();
286// mean_errorTag = profileYTag->GetMeanError();
287// }
288 delete profileY;
289// delete profileYTag;
290
291 av_err_z.push_back(mean);
292 av_err_z_er.push_back(mean_error);
293// av_err_tag_z.push_back(meanTag);
294// av_err_tag_z_er.push_back(mean_errorTag);
295
296 //estimating the approximate k-factor and the error value
297 double pr_er = 0.0;
298 float val(0.);
299 if (fitResKFactorMethod == 1) {
300 pr_er = error_func(bin_count, kgs_z_ntrk_fit_er);
301 } else if (fitResKFactorMethod == 2) {
302 val = h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count);
303 pr_er = TMath::Power(kgs_z_ntrk_fit_er[1] * val, 2) + TMath::Power(kgs_z_ntrk_fit_er[0], 2);
304 pr_er = TMath::Sqrt(pr_er);
305// cout << "val = " << val << ", pr_er = " << pr_er << ", p0er = " << kgs_z_ntrk_fit_er[0] << ", p1er = "<<
306// kgs_z_ntrk_fit_er[1] << endl;
307 } else if (fitResKFactorMethod == 3) {
308 val = h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count);
309 //approximately the error on the plateau
310 pr_er = kgs_z_ntrk_fit_er[2];
311 } else if (fitResKFactorMethod == 4) {
312 pr_er = kgs_z_ntrk_fit_er[0];
313 }
314
315 res_z.push_back(mean * kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count)));
316 res_z_er.push_back(TMath::Sqrt(TMath::Power(mean_error *
317 kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(
318 bin_count)),
319 2) + TMath::Power(pr_er * mean, 2)));
320// res_tag_z.push_back(meanTag * kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count)));
321// res_tag_z_er.push_back(TMath::Sqrt(TMath::Power(mean_errorTag *
322// kgs_z_ntrk_fit->Eval(h_Vrt_err_vs_Something->GetXaxis()->GetBinCenter(bin_count)),2) + TMath::Power( pr_er * mean
323// ,2)));
324 }
325 TGraphErrors* res_z_vs_ntrk =
326 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]));
327 res_z_vs_ntrk->GetYaxis()->SetTitle(coordinate + " Vertex Resolution [mm]");
328 res_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel);
329 res_z_vs_ntrk->SetTitle(coordinate + " Vertex Resolution");
330 res_z_vs_ntrk->SetName("resolution_" + coordinate + "_" + versus);
331
332// TGraphErrors * res_tag_z_vs_ntrk = new TGraphErrors(err_bins_z_nt.size(),
333// &(err_bins_z_nt[0]),&(res_tag_z[0]),&(err_bins_z_nt_er[0]), &(res_tag_z_er[0]) );
334// res_tag_z_vs_ntrk->GetYaxis()->SetTitle(coordinate+" Tagged Vertex Resolution [mm]");
335// res_tag_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel);
336// res_tag_z_vs_ntrk->SetTitle(coordinate+" Tagged Vertex Resolution");
337// res_tag_z_vs_ntrk->SetName("resolution_tag_"+coordinate+"_"+versus);
338
339// TGraphErrors * mean_err_z_vs_ntrk = new TGraphErrors(err_bins_z_nt.size(),
340// &(err_bins_z_nt[0]),&(av_err_z[0]),&(err_bins_z_nt_er[0]), &(av_err_z_er[0]) );
341// mean_err_z_vs_ntrk->GetYaxis()->SetTitle(coordinate+" mean vertex error [mm]");
342// mean_err_z_vs_ntrk->GetXaxis()->SetTitle(xAxisLabel);
343// mean_err_z_vs_ntrk->SetTitle(coordinate+" Mean Vertex Error");
344// mean_err_z_vs_ntrk->SetName("resolution_"+coordinate+"_"+versus+"_Uncorrected"); //not corrected with k-factor
345
347 // Writing out
349 if (versus == "Ntrk") res_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
350 else res_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
351 res_z_vs_ntrk->GetYaxis()->SetRangeUser(0.0025, 1.);
352 res_z_vs_ntrk->Write("", TObject::kOverwrite);
353 delete res_z_vs_ntrk;
354
355 if (versus == "Ntrk") krms_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
356 else krms_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
357 krms_z_vs_ntrk->GetYaxis()->SetRangeUser(0.5, 1.3);
358 krms_z_vs_ntrk->Write("", TObject::kOverwrite);
359 delete krms_z_vs_ntrk;
360
361 if (versus == "Ntrk") kgs_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 100.);
362 else kgs_z_vs_ntrk->GetXaxis()->SetRangeUser(0., 20.);
363 kgs_z_vs_ntrk->GetYaxis()->SetRangeUser(0.5, 1.3);
364 kgs_z_vs_ntrk->Write("", TObject::kOverwrite);
365 delete kgs_z_vs_ntrk;
366
367// nTrksPerVertex->GetYaxis()->SetTitle("Entries");
368// nTrksPerVertex->Draw();
369// nTrksPerVertex->Write("", TObject::kOverwrite);
370
371 return;
372 }
373
374// Reconstruction efficiency
375 void plotEfficiency() {
376// cout << "Creating and writing histos for efficiency" << endl;
377
378 // get histos but return if any histo is not there in input
379 TH1F* h_Vrt_split_tag_ntrk = (TH1F*) gDirectory->Get("Vrt_split_tag_ntrk");
380
381 if (h_Vrt_split_tag_ntrk == 0) return;
382
383 TH1F* h_Vrt_split_probe_ntrk = (TH1F*) gDirectory->Get("Vrt_split_probe_ntrk");
384 if (h_Vrt_split_probe_ntrk == 0) return;
385
386 TH1F* h_Vrt_split_matched_tag_ntrk = (TH1F*) gDirectory->Get("Vrt_split_matched_tag_ntrk");
387 if (h_Vrt_split_matched_tag_ntrk == 0) return;
388
389 TH1F* h_Vrt_split_matched_probe_ntrk = (TH1F*) gDirectory->Get("Vrt_split_matched_probe_ntrk");
390 if (h_Vrt_split_matched_probe_ntrk == 0) return;
391
392 TH1F* h_Vrt_split_dist_tag = (TH1F*) gDirectory->Get("Vrt_split_dist_tag");
393 if (h_Vrt_split_dist_tag == 0) return;
394
395 TH1F* h_Vrt_split_dist_probe = (TH1F*) gDirectory->Get("Vrt_split_dist_probe");
396 if (h_Vrt_split_dist_probe == 0) return;
397
398 // Use BayesDivide routing of TGraphAsymmErrors
399 TGraphAsymmErrors* g_Vrt_rec_eff_m1_split_vs_ntrk = new TGraphAsymmErrors();
400
401 g_Vrt_rec_eff_m1_split_vs_ntrk->BayesDivide(h_Vrt_split_probe_ntrk, h_Vrt_split_tag_ntrk);
402 g_Vrt_rec_eff_m1_split_vs_ntrk->SetName("g_RecEff_M1");
403
404 TGraphAsymmErrors* g_Vrt_sel_eff_m1_split_vs_ntrk = new TGraphAsymmErrors();
405 g_Vrt_sel_eff_m1_split_vs_ntrk->BayesDivide(h_Vrt_split_matched_probe_ntrk, h_Vrt_split_matched_tag_ntrk);
406 g_Vrt_sel_eff_m1_split_vs_ntrk->SetName("g_SelEff_M1");
407
408 // formatting and writing out
409 g_Vrt_rec_eff_m1_split_vs_ntrk->GetHistogram()->GetXaxis()->SetTitle("Number of tracks");
410 g_Vrt_rec_eff_m1_split_vs_ntrk->GetHistogram()->GetYaxis()->SetTitle("Reconstruction efficiency");
411 g_Vrt_rec_eff_m1_split_vs_ntrk->SetMarkerStyle(20);
412 g_Vrt_rec_eff_m1_split_vs_ntrk->Write("", TObject::kOverwrite);
413 delete g_Vrt_rec_eff_m1_split_vs_ntrk;
414
415 g_Vrt_sel_eff_m1_split_vs_ntrk->GetHistogram()->GetXaxis()->SetTitle("Number of tracks");
416 g_Vrt_sel_eff_m1_split_vs_ntrk->GetHistogram()->GetYaxis()->SetTitle("Selection Efficiency");
417 g_Vrt_sel_eff_m1_split_vs_ntrk->SetMarkerStyle(20);
418 g_Vrt_sel_eff_m1_split_vs_ntrk->Write("", TObject::kOverwrite);
419 delete g_Vrt_sel_eff_m1_split_vs_ntrk;
420
421 return;
422 }
423
424 std::vector<float> stableGaussianFit(TH1* histo) {
425 std::vector<float> results;
426 if (histo) {
427 double sigmas = 2.;
428
429 histo->Fit("gaus", "Q0", "", -2, 2);
430 TF1* func = histo->GetFunction("gaus");
431 double actualSigma = func->GetParameter(2);
432 double actualSigmaErr = func->GetParError(2);
433
434 for (int u = 0; u < 10; u++) {
435 histo->Fit("gaus", "Q0", "", -sigmas * actualSigma, sigmas * actualSigma);
436 func = histo->GetFunction("gaus");
437 actualSigma = func->GetParameter(2);
438 actualSigmaErr = func->GetParError(2);
439 }//end of the fitting loop
440
441 results.push_back(actualSigma);
442 results.push_back(actualSigmaErr);
443 } else {
444 results.push_back(0.);
445 results.push_back(0.);
446 }//end of protection against an empty histogram
447
448 return results;
449 }//end of stableGaussian Fit function
450
451 double error_func(float x, const Double_t* par) {
452//calculating the square of the propagated error on the fit values
453 return(TMath::Power(par[0], 2) + x * TMath::Power(par[1], 2) + TMath::Power(x * par[2], 2));
454 }
455
456 double scaleFactorFitFcn(double* x, double* par) {
457 // Truncated gaus-smeared step function
458 // par 0: mean of turn-on (and truncation point)
459 // par 1: slope of turn-on
460 // par 2: plateau
461 return (x[0] >= par[0]) * TMath::Erf(par[1] * x[0] - par[0]) * par[2];
462 }
463}
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:424
dqutils::error_func
double error_func(float x, const Double_t *par)
Definition MonitoringFile_IDEnhancedPrimaryVertex.cxx:451
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:456
Generated on for ATLAS Offline Software by doxygen 1.14.0