ATLAS Offline Software
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 
27 namespace 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 }
python.CaloRecoConfig.f
f
Definition: CaloRecoConfig.py:127
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
TH2F
Definition: rootspy.cxx:420
TH1D
Definition: rootspy.cxx:342
NSWL1::coordinate
float coordinate(const Vertex &v)
Definition: GeoUtils.h:56
x
#define x
Trk::u
@ u
Enums for curvilinear frames.
Definition: ParamDefs.h:83
binWidth
void binWidth(TH1 *h)
Definition: listroot.cxx:80
ParseInputs.gDirectory
gDirectory
Definition: Final2012/ParseInputs.py:133
dqutils::scaleFactorFitFcn
double scaleFactorFitFcn(double *x, double *par)
Definition: MonitoringFile_IDEnhancedPrimaryVertex.cxx:456
dqutils::plotResolution
void plotResolution(const TString &coordinate, const TString &versus)
Definition: MonitoringFile_IDEnhancedPrimaryVertex.cxx:107
dqutils
Definition: CoolMdt.h:76
createCoolChannelIdFile.par
par
Definition: createCoolChannelIdFile.py:29
dqutils::stableGaussianFit
std::vector< float > stableGaussianFit(TH1 *histo)
Definition: MonitoringFile_IDEnhancedPrimaryVertex.cxx:424
python.ami.results
def results
Definition: ami.py:386
MonitoringFile.h
TH1F
Definition: rootspy.cxx:320
Pythia8_RapidityOrderMPI.val
val
Definition: Pythia8_RapidityOrderMPI.py:14
TH1
Definition: rootspy.cxx:268
dqutils::MonitoringFile::pv_PrimaryVertexMonitoring_calcResoAndEfficiency
static void pv_PrimaryVertexMonitoring_calcResoAndEfficiency(const std::string &inFilename, bool isIncremental=false)
Definition: MonitoringFile_IDEnhancedPrimaryVertex.cxx:35
dqutils::plotEfficiency
void plotEfficiency()
Definition: MonitoringFile_IDEnhancedPrimaryVertex.cxx:375
plotBeamSpotCompare.histo
histo
Definition: plotBeamSpotCompare.py:415
dqutils::error_func
double error_func(float x, const Double_t *par)
Definition: MonitoringFile_IDEnhancedPrimaryVertex.cxx:451
Generated on Tue Jul 2 2024 21:22:02 for ATLAS Offline Software by   doxygen 1.8.18