absCal.cxx
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////
2 // //
3 // Absolute calibratrion parameters //
4 // //
5 // Created by: Luke Vinton < July/2015 //
6 // Modified by: Diana Patricia Mendez > July 2015 ... //
7 // //
8 // Output: cMCdataND.pdf - Corrected response/cm //
9 // = (PECorr/cm) = MEU //
10 // cMCND_truth.pdf - Simulated dE/dx (MeV/cm) //
11 // cMCdataND_cal.pdf - Calibrated dE/dx (MeV/cm) //
12 // //
13 // Calibrated dE/dx is multiplied by the conversion factor //
14 // //
15 ////////////////////////////////////////////////////////////////
16 
17 #include <iostream>
18 #include <vector>
19 
20 #include <TFile.h>
21 #include <TPad.h>
22 #include <TH1.h>
23 #include <TH2.h>
24 #include <TH2F.h>
25 #include <TH2D.h>
26 #include <TProfile2D.h>
27 #include <TH1D.h>
28 #include <TH1F.h>
29 #include <TH3D.h>
30 #include <TCanvas.h>
31 #include <TVector.h>
32 #include <TColor.h>
33 #include <TLegend.h>
34 #include <TF1.h>
35 #include <TLatex.h>
36 #include <TMath.h>
37 #include <TStyle.h>
38 #include <TROOT.h>
39 #include <TPad.h>
40 #include <TFile.h>
41 #include <TTree.h>
42 #include <TSystem.h>
43 #include <TColor.h>
44 #include "TChain.h"
45 
46 #include "Utilities/rootlogon.C"
47 
48 // canvas size
49 #define c1_x 700
50 #define c1_y 500
51 
52 #define c2_x 1000
53 #define c2_y 500
54 
55 
56 
57 
58 double Landau(TH1D* h, double& err);
59 double Gaus(TH1D* h, double& err, bool isTruth);
60 double median(TH1D* h);
61 std::vector<double> quantiles( TH1D* h );
62 
63 double errAoverB(double& err, double A, double B, double Aerr, double Berr);
64 double GetHistMax(TH1* h0, TH1* h1, TH1* h2, TH1* h3);
65 
66 
67 int main(){
68 
69  //This constants refer to the MEU value for mc and data
70  //This should be changed to get the right calibrated plot
71  double const mccal = 0.0451381;
72  double const datacal = 0.0442841;
73 
74  gROOT->ForceStyle();
75 
76  //Using TChain avoids using merge and then having an only input file
77  //which sometimes do not work due to size and nunber or files to merge
78  //specially for FD
79  std::cout << "Making MC chain ..." << std::endl;
80  TChain* tNDMC = new TChain("muondedxana/fTree");
81  tNDMC->Add("mc/*.root");
82 
83  double PECorr_MC, path_MC, trueE_MC, x_MC;
84  double planes_MC, xplanes_MC, yplanes_MC;
85  int cont_MC, nhits_MC,hitId_MC;
86 
87  //Declaring branches saves time compared to the Draw function
88  std::cout << "Setting MC branches ..." <<std::endl;
89  tNDMC->SetBranchAddress("PECorr", &PECorr_MC);
90  tNDMC->SetBranchAddress("path" , &path_MC);
91  tNDMC->SetBranchAddress("trueE" , &trueE_MC);
92  tNDMC->SetBranchAddress("x" , &x_MC);
93 
94 
95  std::cout << "Declaring MC histograms ..." <<std::endl;
96  TH1D* hMCND_PECorr = new TH1D("hMCND_PECorr", "hMCND_PECorr", 1000, 0, 500);
97  TH1D* hMCND_trueE = new TH1D("hMCND_trueE", "hMCND_trueE", 1000, 0, 50);
98  TH1D* hMCND_cal = new TH1D("hMCND_cal", "hMCND_cal", 60, 0, 6);
99 
100  std::cout << "Filling MC histos ..." <<std::endl;
101  unsigned int nEntries_MC = tNDMC->GetEntries();
102  for(unsigned int i=0; i<nEntries_MC; i++){
103  tNDMC->GetEntry(i);
104 
105  if( x_MC>100 && x_MC<200){
106  if(PECorr_MC>0){
107  hMCND_PECorr ->Fill(PECorr_MC/path_MC);
108  hMCND_cal ->Fill(mccal*PECorr_MC/path_MC);
109  }
110  if(trueE_MC>0){
111  hMCND_trueE ->Fill(trueE_MC/path_MC);
112  }
113  }
114 
115  }
116 
117 
118 
119  std::cout << "Making data chain ..." <<std::endl;
120  TChain* tNDdata = new TChain("muondedxana/fTree");
121  tNDdata->Add("data/*.root");
122 
123  double PECorr_data, path_data, trueE_data, x_data;
124  double planes_data, xplanes_data, yplanes_data;
125  int cont_data, nhits_data,hitId_data;
126 
127  std::cout << "Setting data branches ..." <<std::endl;
128  tNDdata->SetBranchAddress("PECorr", &PECorr_data);
129  tNDdata->SetBranchAddress("path" , &path_data);
130  tNDdata->SetBranchAddress("trueE" , &trueE_data);
131  tNDdata->SetBranchAddress("x" , &x_data);
132 
133  std::cout << "Declaring data histos ..." <<std::endl;
134 
135  TH2D* hdataND_PECorrX = new TH2D("hdataND_PECorrX", "hdataND_PECorrX", 160, -800, 800, 30, 0, 60);
136  TH1D* hdataND_PECorr = new TH1D("hdataND_PECorr","hdataND_PECorr", 1000, 0, 500);
137  TH1D* hdataND_cal = new TH1D("hdataND_cal","hdataND_cal", 60, 0, 6);
138 
139  std::cout << "Filling data histos ..." <<std::endl;
140 
141  unsigned int nEntries_data = tNDdata->GetEntries();
142  for(unsigned int i=0; i<nEntries_data; i++){
143  tNDdata->GetEntry(i);
144 
145  // if(hitId_data!=2) continue;
146  // if(!PassActivityTrigger(planes_data, xplanes_data, yplanes_data, nhits_data, cont_data)) continue;
147 
148  if(x_data>100 && x_data<200)
149  if(PECorr_data>0){
150  hdataND_PECorrX ->Fill(x_data,PECorr_data/path_data);
151  hdataND_PECorr ->Fill(PECorr_data/path_data);
152  hdataND_cal ->Fill(datacal*PECorr_data/path_data);
153 
154  }
155 
156  }
157 
158 
159 
160 
161  //landau
162  TH1D* hdataND_PECorr_landau = (TH1D*)hdataND_PECorr->Clone("hdataND_PECorr_landau");
163  TH1D* hMCND_PECorr_landau = (TH1D*)hMCND_PECorr ->Clone("hMCND_PECorr_landau");
164  TH1D* hMCND_trueE_landau = (TH1D*)hMCND_trueE ->Clone("hMCND_trueE_landau");
165  //gaus
166  TH1D* hdataND_PECorr_gaus = (TH1D*)hdataND_PECorr->Clone("hdataND_PECorr_gaus");
167  TH1D* hMCND_PECorr_gaus = (TH1D*)hMCND_PECorr ->Clone("hMCND_PECorr_gaus");
168  TH1D* hMCND_trueE_gaus = (TH1D*)hMCND_trueE ->Clone("hMCND_trueE_gaus");
169 
170  //set data lines to black:
171  hdataND_PECorr -> SetLineColor(1);
172  hdataND_PECorr_landau -> SetLineColor(1);
173  hdataND_PECorr_gaus -> SetLineColor(1);
174 
175  //set histo range
176  hdataND_PECorr -> GetXaxis() -> SetRangeUser(0,100);
177  hMCND_PECorr -> GetXaxis() -> SetRangeUser(0,100);
178  hMCND_trueE -> GetXaxis() -> SetRangeUser(0,6);
179  hdataND_cal -> GetXaxis() -> SetRangeUser(0,6);
180  hMCND_cal -> GetXaxis() -> SetRangeUser(0,6);
181 
182 
183  double mpvErr_dataND, mpvErr_MCND, mpvErr_trueE_MCND;
184  double mpv_dataND = Landau(hdataND_PECorr_landau, mpvErr_dataND);
185  double mpv_MCND = Landau(hMCND_PECorr_landau, mpvErr_MCND);
186  double mpv_trueE_MCND = Landau(hMCND_trueE_landau, mpvErr_trueE_MCND);
187 
188  double gausErr_dataND, gausErr_MCND, gausErr_trueE_MCND;
189  double gausMean_dataND = Gaus(hdataND_PECorr_gaus, gausErr_dataND, false);
190  double gausMean_MCND = Gaus(hMCND_PECorr_gaus, gausErr_MCND, false);
191  double gausMean_trueE_MCND = Gaus(hMCND_trueE_gaus, gausErr_trueE_MCND, true);
192 
193  double mean_dataND = hdataND_PECorr->GetMean(1);
194  double mean_MCND = hMCND_PECorr->GetMean(1);
195  double mean_trueE_MCND = hMCND_trueE->GetMean(1);
196  double meanErr_dataND = hdataND_PECorr->GetRMS()/(pow(hdataND_PECorr->GetEntries(), 0.5));
197  double meanErr_MCND = hMCND_PECorr->GetRMS()/(pow(hMCND_PECorr->GetEntries(), 0.5));
198  double meanErr_trueE_MCND = hMCND_trueE->GetRMS()/(pow(hMCND_trueE->GetEntries(), 0.5));
199 
200  double Median_dataND;
201  double pointFive = 0.5;
202  hdataND_PECorr->GetQuantiles(1,&Median_dataND,&pointFive);
203  double Median_MCND;
204  hMCND_PECorr->GetQuantiles(1,&Median_MCND,&pointFive);
205  double Median_trueE_MCND;
206  hMCND_trueE->GetQuantiles(1,&Median_trueE_MCND,&pointFive);
207 
208 
209  // print MEU values from all techniques along with errors
210  std::cout<< " ----------- Printing MEU values: --------" <<std::endl;
211  std::cout<< "*** Mean: *** " <<std::endl;
212  std::cout<< "dataND: " << mean_dataND << " +/- " << meanErr_dataND << " (stat.)" <<std::endl;
213  std::cout<< "MCND: " << mean_MCND << " +/- " << meanErr_MCND << " (stat.)" <<std::endl;
214  std::cout<< "trueE_MCND: " << mean_trueE_MCND << " +/- " << meanErr_trueE_MCND << " (stat.)" <<std::endl;
215 
216  std::cout<< "*** Landau: *** " <<std::endl;
217  std::cout<< "dataND: " << mpv_dataND << " +/- " << mpvErr_dataND << " (stat.)" <<std::endl;
218  std::cout<< "MCND: " << mpv_MCND << " +/- " << mpvErr_MCND << " (stat.)" <<std::endl;
219  std::cout<< "trueE_MCND: " << mpv_trueE_MCND << " +/- " << mpvErr_trueE_MCND << " (stat.)" <<std::endl;
220 
221  std::cout<< "*** Gaus: *** " <<std::endl;
222  std::cout<< "dataND: " << gausMean_dataND << " +/- " << gausErr_dataND << " (stat.)" <<std::endl;
223  std::cout<< "MCND: " << gausMean_MCND << " +/- " << gausErr_MCND << " (stat.)" <<std::endl;
224  std::cout<< "trueE_MCND: " << gausMean_trueE_MCND << " +/- " << gausErr_trueE_MCND << " (stat.)" <<std::endl;
225 
226 
227 
228  //abs. scales with diff. techniques
229  //landau:
230  double abs_landau_dataND = mpv_trueE_MCND/mpv_dataND;
231  double abs_landau_MCND = mpv_trueE_MCND/mpv_MCND;
232  double absErr_landau_dataND, absErr_landau_MCND;
233  errAoverB(absErr_landau_dataND, mpv_trueE_MCND, mpv_dataND, mpvErr_trueE_MCND, mpvErr_dataND);
234  errAoverB(absErr_landau_MCND, mpv_trueE_MCND, mpv_MCND, mpvErr_trueE_MCND, mpvErr_MCND);
235  //gaus:
236  double abs_gaus_dataND = gausMean_trueE_MCND/gausMean_dataND;
237  double abs_gaus_MCND = gausMean_trueE_MCND/gausMean_MCND;
238  double absErr_gaus_dataND, absErr_gaus_MCND;
239  errAoverB(absErr_gaus_dataND, gausMean_trueE_MCND, gausMean_dataND, gausErr_trueE_MCND, gausErr_dataND);
240  errAoverB(absErr_gaus_MCND, gausMean_trueE_MCND, gausMean_MCND, gausErr_trueE_MCND, gausErr_MCND);
241  //mean:
242  double abs_mean_dataND = mean_trueE_MCND/mean_dataND;
243  double abs_mean_MCND = mean_trueE_MCND/mean_MCND;
244  double absErr_mean_dataND, absErr_mean_MCND;
245  errAoverB(absErr_mean_dataND, mean_trueE_MCND, mean_dataND, meanErr_trueE_MCND, meanErr_dataND);
246  errAoverB(absErr_mean_MCND, mean_trueE_MCND, mean_MCND, meanErr_trueE_MCND, meanErr_MCND);
247  //median:
248  // Need to find technique to calculate error on the median...
249  double abs_median_dataND = Median_trueE_MCND/Median_dataND;
250  double abs_median_MCND = Median_trueE_MCND/Median_MCND;
251 
252  // Print Absolute scale values along with statistical errors
253  std::cout<< "------------- Printing Abs. scale values: ----------- " <<std::endl;
254  std::cout<< "*** Mean: *** " <<std::endl;
255  std::cout<< "dataND: " << abs_mean_dataND << " +/- " << absErr_mean_dataND << " (stat.)" <<std::endl;
256  std::cout<< "MCND: " << abs_mean_MCND << " +/- " << absErr_mean_MCND << " (stat.)" <<std::endl;
257  std::cout<< "*** Median: *** " <<std::endl;
258  std::cout<< "dataND: " << abs_median_dataND << " +/- " << absErr_mean_dataND << " (stat.)" <<std::endl;
259  std::cout<< "MCND: " << abs_median_MCND << " +/- " << absErr_mean_MCND << " (stat.)" <<std::endl;
260 
261  std::cout<< "*** Landau: *** " <<std::endl;
262  std::cout<< "dataND: " << abs_landau_dataND << " +/- " << absErr_landau_dataND << " (stat.)" <<std::endl;
263  std::cout<< "MCND: " << abs_landau_MCND << " +/- " << absErr_landau_MCND << " (stat.)" <<std::endl;
264  std::cout<< "*** Gaus: *** " <<std::endl;
265  std::cout<< "dataND: " << abs_gaus_dataND << " +/- " << absErr_gaus_dataND << " (stat.)" <<std::endl;
266  std::cout<< "MCND: " << abs_gaus_MCND << " +/- " << absErr_gaus_MCND << " (stat.)" <<std::endl;
267 
268 
269 
270 
271  std::ofstream values("absCal.txt");
272 
273  //printing values in created file
274  values<< " ----------- Printing MEU values: --------" <<std::endl;
275  values<< "*** Mean: *** " <<std::endl;
276  values<< "dataND: " << mean_dataND << " +/- " << meanErr_dataND << " (stat.)" <<std::endl;
277  values<< "MCND: " << mean_MCND << " +/- " << meanErr_MCND << " (stat.)" <<std::endl;
278  values<< "trueE_MCND: " << mean_trueE_MCND << " +/- " << meanErr_trueE_MCND << " (stat.)" <<std::endl;
279 
280  values<< "*** Landau: *** " <<std::endl;
281  values<< "dataND: " << mpv_dataND << " +/- " << mpvErr_dataND << " (stat.)" <<std::endl;
282  values<< "MCND: " << mpv_MCND << " +/- " << mpvErr_MCND << " (stat.)" <<std::endl;
283  values<< "trueE_MCND: " << mpv_trueE_MCND << " +/- " << mpvErr_trueE_MCND << " (stat.)" <<std::endl;
284 
285  values<< "*** Gaus: *** " <<std::endl;
286  values<< "dataND: " << gausMean_dataND << " +/- " << gausErr_dataND << " (stat.)" <<std::endl;
287  values<< "MCND: " << gausMean_MCND << " +/- " << gausErr_MCND << " (stat.)" <<std::endl;
288  values<< "trueE_MCND: " << gausMean_trueE_MCND << " +/- " << gausErr_trueE_MCND << " (stat.)" <<std::endl;
289 
290  values<< "------------- Printing Abs. scale values: ----------- " <<std::endl;
291  values<< "*** Mean: *** " <<std::endl;
292  values<< "dataND: " << abs_mean_dataND << " +/- " << absErr_mean_dataND << " (stat.)" <<std::endl;
293  values<< "MCND: " << abs_mean_MCND << " +/- " << absErr_mean_MCND << " (stat.)" <<std::endl;
294  values<< "*** Median: *** " <<std::endl;
295  values<< "dataND: " << abs_median_dataND << " +/- " << absErr_mean_dataND << " (stat.)" <<std::endl;
296  values<< "MCND: " << abs_median_MCND << " +/- " << absErr_mean_MCND << " (stat.)" <<std::endl;
297 
298  values<< "*** Landau: *** " <<std::endl;
299  values<< "dataND: " << abs_landau_dataND << " +/- " << absErr_landau_dataND << " (stat.)" <<std::endl;
300  values<< "MCND: " << abs_landau_MCND << " +/- " << absErr_landau_MCND << " (stat.)" <<std::endl;
301  values<< "*** Gaus: *** " <<std::endl;
302  values<< "dataND: " << abs_gaus_dataND << " +/- " << absErr_gaus_dataND << " (stat.)" <<std::endl;
303  values<< "MCND: " << abs_gaus_MCND << " +/- " << absErr_gaus_MCND << " (stat.)" <<std::endl;
304 
305  values.close();
306 
307 
308  //////////// plotting ///////////
309 
310  TCanvas *cDataND = new TCanvas("cDataND","cDataND",c1_x,c1_y);
311  cDataND->cd();
312  hdataND_PECorr->SetTitle("");
313  hdataND_PECorr->GetXaxis()->SetTitle("Calibrated detector response/cm");
314  hdataND_PECorr->GetYaxis()->SetTitle("Hits in track window");
315  CenterTitles(hdataND_PECorr);
316  hdataND_PECorr->Draw();
317  Preliminary();
318 
319  TCanvas *cDataND_landau = new TCanvas("cDataND_landau","cDataND_landau", c1_x, c1_y);
320  cDataND_landau->cd();
321  hdataND_PECorr_landau->GetXaxis()->SetTitle("Calibrated detector response/cm");
322  hdataND_PECorr_landau->SetTitle("");
323  CenterTitles(hdataND_PECorr_landau);
324  hdataND_PECorr_landau->Draw();
325  Preliminary();
326 
327 
328  TCanvas *cDataND_gaus = new TCanvas("cDataND_gaus","cDataND_gaus",c1_x ,c1_y);
329  cDataND_gaus->cd();
330  hdataND_PECorr_gaus->GetXaxis()->SetTitle("Calibrated detector response/cm");
331  hdataND_PECorr_gaus->SetTitle("");
332  CenterTitles(hdataND_PECorr_gaus);
333  hdataND_PECorr_gaus->Draw();
334  Preliminary();
335 
336 
337  TCanvas *cMCND = new TCanvas("cMCND","cMCND", c1_x, c1_y);
338  cMCND -> cd();
339  hMCND_PECorr->SetTitle("");
340  hMCND_PECorr->GetXaxis()->SetTitle("dE/dx (MeV/cm)");
341  hMCND_PECorr->GetYaxis()->SetTitle("Hits in track window");
342  hMCND_PECorr->SetLineColor(2);
343  CenterTitles(hMCND_PECorr);
344  hMCND_PECorr->Draw();
345  Simulation();
346 
347  cMCND->Print("cMCND.pdf");
348 
349 
350  TCanvas *cMCdataND = new TCanvas("cMCdataND","cMCdataND", c1_x, c1_y);
351  cMCdataND -> cd();
352  hMCND_PECorr->SetTitle("");
353  hMCND_PECorr->GetXaxis()->SetTitle("Corrected response/cm (PECorr/cm)");
354  hMCND_PECorr->GetYaxis()->SetTitle("Hits in track window");
355 
356  hMCND_PECorr->SetLineColor(2);
357  CenterTitles(hMCND_PECorr);
358  hMCND_PECorr->Draw();
359  hdataND_PECorr->Scale(hMCND_PECorr->Integral() / hdataND_PECorr->Integral());
360  // Call below is commented out because I cannot find a definition for
361  // SetTH1Max anywhere. BJR
362  //SetTH1Max(hMCND_PECorr, hdataND_PECorr);
363  hdataND_PECorr->SetTitle("");
364  CenterTitles(hdataND_PECorr);
365  hdataND_PECorr->Draw("same");
366  Preliminary();
367 
368  TLegend lMCdataND(0.6,0.55,0.85,0.7);
369  lMCdataND.AddEntry(hdataND_PECorr,"FD Data","l");
370  lMCdataND.AddEntry(hMCND_PECorr, "FD MC","l");
371  lMCdataND.SetBorderSize(0);
372  lMCdataND.SetTextFont(42);
373  lMCdataND.DrawClone();
374 
375  cMCdataND->Print("cMCdataND.pdf");
376 
377 
378  TCanvas *cMCND_landau = new TCanvas("cMCND_landau","cMCND_landau", c1_x, c1_y);
379  cMCND_landau -> cd();
380  hMCND_PECorr_landau->SetTitle("");
381  hMCND_PECorr_landau->GetXaxis()->SetTitle("Calibrated detector response/cm");
382  hMCND_PECorr_landau->SetLineColor(2);
383  CenterTitles(hMCND_PECorr_landau);
384  hMCND_PECorr_landau->Draw();
385  Simulation();
386 
387 
388  TCanvas *cMCND_gaus = new TCanvas("cMCND_gaus","cMCND_gaus",c1_x,c1_y);
389  cMCND_gaus -> cd();
390  hMCND_PECorr_gaus->SetTitle("");
391  hMCND_PECorr_gaus->GetXaxis()->SetTitle("Calibrated detector response/cm");
392  hMCND_PECorr_gaus->SetLineColor(2);
393  CenterTitles(hMCND_PECorr_gaus);
394  hMCND_PECorr_gaus->Draw();
395  Simulation();
396 
397 
398  TCanvas *cMCND_truth = new TCanvas("cMCND_truth","cMCND_truth",c1_x,c1_y);
399  cMCND_truth -> cd();
400  gStyle->SetOptStat(0);
401  hMCND_trueE->SetTitle("");
402  hMCND_trueE->GetXaxis()->SetTitle("Simulated dE/dx (MeV/cm)");
403  hMCND_trueE->GetYaxis()->SetTitle("Hits in track window");
404  hMCND_trueE->SetLineColor(2);
405  CenterTitles(hMCND_trueE);
406  hMCND_trueE->Draw();
407  Simulation();
408 
409  TLegend lTruthND(0.5,0.55,0.8,0.65);
410  lTruthND.AddEntry(hMCND_trueE,"FD MC","l");
411  lTruthND.SetBorderSize(0);
412  lTruthND.SetTextFont(42);
413  lTruthND.DrawClone();
414 
415  cMCND_truth->Print("cMCND_truth.pdf");
416 
417 
418  TCanvas *cMCdataND_cal = new TCanvas("cMCdataND_cal","cMCdataND_cal", c1_x, c1_y);
419  cMCdataND_cal -> cd();
420  hMCND_cal->SetTitle("");
421  hMCND_cal->GetXaxis()->SetTitle("Calibrated dE/dx (MeV/cm)");
422  hMCND_cal->GetYaxis()->SetTitle("Hits in track window");
423 
424  hMCND_cal->SetLineColor(2);
425  CenterTitles(hMCND_cal);
426  hMCND_cal->SetTitle("");
427  hMCND_cal->Draw();
428  hdataND_cal->Scale(hMCND_cal->Integral() / hdataND_cal->Integral());
429  // Call below is commented out because I cannot find a definition for
430  // SetTH1Max anywhere. BJR
431  //SetTH1Max(hMCND_cal, hdataND_cal);
432  CenterTitles(hdataND_cal);
433  hdataND_cal->Draw("same");
434  Preliminary();
435 
436  TLegend lMCdataND_cal(0.5,0.55,0.8,0.7);
437  lMCdataND_cal.AddEntry(hdataND_cal,"FD Data","l");
438  lMCdataND_cal.AddEntry(hMCND_cal, "FD MC","l");
439  lMCdataND_cal.SetBorderSize(0);
440  lMCdataND_cal.SetTextFont(42);
441  lMCdataND_cal.DrawClone();
442 
443  cMCdataND_cal->Print("cMCdataND_cal.pdf");
444 
445  TCanvas *cMCND_truth_landau = new TCanvas("cMCND_truth_landau","cMCND_truth_landau", c1_x, c1_y);
446  cMCND_truth_landau -> cd();
447  hMCND_trueE_landau->SetTitle("");
448  hMCND_trueE_landau->GetXaxis()->SetTitle("MeV/cm");
449  hMCND_trueE_landau->SetLineColor(2);
450  CenterTitles(hMCND_trueE_landau);
451  hMCND_trueE_landau->Draw();
452  Simulation();
453 
454 
455  TCanvas *cMCND_truth_gaus = new TCanvas("cMCND_truth_gaus","cMCND_truth_gaus", c1_x, c1_y);
456  cMCND_truth_gaus -> cd();
457 
458  hMCND_trueE_gaus->SetTitle("");
459  hMCND_trueE_gaus->GetXaxis()->SetTitle("MeV/cm");
460  hMCND_trueE_gaus->SetLineColor(2);
461  CenterTitles(hMCND_trueE_gaus);
462  hMCND_trueE_gaus->Draw();
463  Simulation();
464 
465 
466 }
467 
468 
469 
470 
471 //landau
472 double Landau(TH1D* h, double& err){
473  int bin1 = h -> FindFirstBinAbove(h -> GetMaximum()/2);
474  int bin2 = h -> FindLastBinAbove(h -> GetMaximum()/2);
475  double xLow = h -> GetBinCenter(bin1);
476  double xHigh = h -> GetBinCenter(bin2);
477 
478  //fit with landau
479  TF1 *f1 = new TF1("f1","landau",xLow,xHigh);
480  h -> Fit(f1,"","",xLow,xHigh);
481 
482  double mpv = f1->GetParameter(1);
483  err = f1->GetParError(1);
484  std::cout<< "err: " << err <<std::endl;
485  return mpv;
486 }
487 
488 
489 double Gaus(TH1D* h, double& err, bool isTruth){
490 
491  Double_t parameters[3];
492  parameters[0] = h->GetBinContent(h->GetMaximumBin());
493  parameters[1] = h->GetBinCenter(h->GetMaximumBin());
494  parameters[2] = h->GetRMS();
495 
496  // only care about finding peak position, use gaussian truncated about mode
497  //Double_t bound_size = (isTruth) ? 0.2 : 0.5; // narrower for truth, *not* MC reco... don't mess this up.
498 
499  Double_t bound_size;
500  if(isTruth) bound_size = 0.2;
501  else bound_size = 0.5;
502 
503  std::cout<<" ****** bound size: " << bound_size <<std::endl;
504 
505  //Double_t lowerBound = max(0,parameters[1] - bound_size*h->GetRMS());
506  Double_t lowerBound;
507  if(parameters[1] - bound_size*h->GetRMS() > 0)
508  lowerBound = parameters[1] - bound_size*h->GetRMS();
509  else lowerBound = 0;
510 
511  Double_t upperBound = parameters[1] + bound_size*h->GetRMS();
512 
513  TF1* fitfunction = new TF1("fitfunction","[0]*TMath::Gaus(x,[1],[2])",lowerBound,upperBound);
514  fitfunction->SetParameters(parameters);
515 
516  // fitting
517  h->Fit(fitfunction,"RQ");
518 
519  // retrieve values
520  err = fitfunction->GetParError(1);
521  return fitfunction->GetParameter(1);
522 }
523 
524 double median(TH1D* h){
525  return quantiles(h)[2];
526 }
527 
528 std::vector<double> quantiles( TH1D* h ){
529  //=================================
530  double q[5];
531  double probs[5] = {0.025, 0.16, 0.5, 1 - 0.16, 0.975 };
532  h->GetQuantiles( 5, q, probs );
533  std::vector<double> r(5);
534  for (int i=0; i<5; i++)
535  {
536 
537  r[i] = q[i];
538  }
539  return r;
540 }
541 
542 double errAoverB(double& err, double A, double B, double Aerr, double Berr){
543 
544  double C = A/B;
545  //std::cout<< "C: " << C <<std::endl;
546 
547  err = std::sqrt( std::pow((Aerr/B),2) + std::pow((Berr*A)/(B*B),2) );
548  //double err = std::sqrt(std::pow(C,2)) * std::sqrt( std::pow((Berr/B),2) + std::pow((Aerr/A),2));
549 
550  //std::cout<< "err: " << err <<std::endl;
551  return err;
552 
553 }
554 
555 
556 double GetHistMax(TH1* h0, TH1* h1, TH1* h2, TH1* h3){
557  double max[4];
558  max[0] = h0->GetMaximum();
559  max[1] = h1->GetMaximum();
560  max[2] = h2->GetMaximum();
561  max[3] = h3->GetMaximum();
562 
563  double hiy = 0;
564  for(int i=1;i<4;i++){
565  if(max[i]>hiy)
566  hiy = max[i];
567  }
568  return (hiy * 1.1);
569 }
570 
571 int absCal(){
572  return main();
573 }
574 
575 
void Simulation()
Definition: tools.h:16
std::vector< double > quantiles(TH1D *h)
Definition: absCal.cxx:528
double errAoverB(double &err, double A, double B, double Aerr, double Berr)
Definition: absCal.cxx:542
int main()
Definition: absCal.cxx:67
TH1F * h3
Definition: berger.C:36
fVtxDx Fit("f")
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
const double C
double Landau(TH1D *h, double &err)
Definition: absCal.cxx:472
void CenterTitles(TH1 *h)
correl_xv GetXaxis() -> SetDecimals()
hmean SetLineColor(4)
Float_t f1
double Gaus(TH1D *h, double &err, bool isTruth)
Definition: absCal.cxx:489
#define c1_x
Definition: absCal.cxx:49
double median(TH1D *h)
Definition: absCal.cxx:524
#define c1_y
Definition: absCal.cxx:50
void Preliminary()
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
double GetHistMax(TH1 *h0, TH1 *h1, TH1 *h2, TH1 *h3)
Definition: absCal.cxx:556
TH1F * h1
static const double A
Definition: Units.h:82
TRandom3 r(0)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
int absCal()
Definition: absCal.cxx:571
c cd(1)