absCal_Diblocks.cxx
Go to the documentation of this file.
1 ///////////////////////////////////////////////////////////////////
2 // //
3 // Absolute calibratrion parameters //
4 // FD Calibration divided by diblocks //
5 // Planes per diblick = 64, diblocks = 14 //
6 // //
7 // The calibration group is starting to study the //
8 // posibility of calibrating by regions //
9 // This macro is a demo about how to compute the energy //
10 // scale factors for each FD diblock; is just a modfified //
11 // version of absCal.cxx with a simple array to store the //
12 // values for the different regions //
13 // //
14 // Created by: Diana Patricia Mendez, May/17/2016 //
15 // //
16 // Output: .txt file //
17 // Corrected detector response (PECorr/cm) = MEU //
18 // Truth energy deposited - Simulated dE/dx (MeV/cm) //
19 // Energy scale factors //
20 // Everything with statistical errors //
21 ///////////////////////////////////////////////////////////////////
22 
23 #include <iostream>
24 #include <vector>
25 #include <fstream>
26 
27 #include <TFile.h>
28 #include <TPad.h>
29 #include <TH1.h>
30 #include <TH2.h>
31 #include <TH2F.h>
32 #include <TH2D.h>
33 #include <TProfile2D.h>
34 #include <TH1D.h>
35 #include <TH1F.h>
36 #include <TH3D.h>
37 #include <TCanvas.h>
38 #include <TVector.h>
39 #include <TColor.h>
40 #include <TLegend.h>
41 #include <TF1.h>
42 #include <TLatex.h>
43 #include <TMath.h>
44 #include <TStyle.h>
45 #include <TROOT.h>
46 #include <TPad.h>
47 #include <TFile.h>
48 #include <TTree.h>
49 #include <TSystem.h>
50 #include <TColor.h>
51 #include "TChain.h"
52 
53 #include "Utilities/rootlogon.C"
54 
55 // canvas size
56 #define c1_x 700
57 #define c1_y 500
58 
59 #define c2_x 1000
60 #define c2_y 500
61 
62 
63 
64 
65 double Landau(TH1D* h, double& err);
66 double Gaus(TH1D* h, double& err, bool isTruth);
67 double median(TH1D* h);
68 std::vector<double> quantiles( TH1D* h );
69 
70 double errAoverB(double& err, double A, double B, double Aerr, double Berr);
71 double GetHistMax(TH1* h0, TH1* h1, TH1* h2, TH1* h3);
72 
73 
74 int main(){
75 
76  //main output file name with mean values
77 
78  TString detector = "FD";
79  TString period = "period2";
80  int const diblocks = 14;
81 
82  //This constants refer to the MEU value for mc and data
83  //This should be changed to get the right calibrated plot
84  double mccal[diblocks];
85  double datacal[diblocks];
86 
87 
88  //Using TChain avoids using "merge" and then having an only input file
89  //which sometimes does not work due to size and number of files
90  std::cout << "Making MC chain ..." << std::endl;
91  TChain* tMC = new TChain("muondedxana/fTree");
92  tMC->Add("/nova/ana/users/dmendez/Calibration/Files/FDmc/prod_pcliststop_R16-01-27-prod2calib.e_fd_cry_all_nova_v08_ideal_batch1_v1_gain100/*.root");
93 
94  double PECorr_MC, path_MC, trueE_MC, x_MC;
95  double planes_MC, xplanes_MC, yplanes_MC;
96  int cont_MC, nhits_MC,hitId_MC, db_MC;
97 
98  //Declaring branches saves time compared to the Draw function
99  std::cout << "Setting MC branches ..." <<std::endl;
100  tMC->SetBranchAddress("PECorr", &PECorr_MC);
101  tMC->SetBranchAddress("path" , &path_MC);
102  tMC->SetBranchAddress("trueE" , &trueE_MC);
103  tMC->SetBranchAddress("x" , &x_MC);
104  tMC->SetBranchAddress("diblock",&db_MC);
105 
106 
107  std::cout << "Declaring MC histograms ..." <<std::endl;
108  // adding an extra histogram for calibration over the whole detector
109  // histogram 0 is for all the detector
110  TH1D *hMC_PECorr[diblocks+1];
111  TH1D *hMC_trueE[diblocks+1];
112  TH1D *hMC_cal[diblocks+1];
113 
114  char name[20];
115  char title[20];
116 
117  for (unsigned int i=0;i<diblocks+1;i++) {
118  sprintf(title,"DiBlock_%d",i);
119 
120  sprintf(name,"MC_PECorr_%d",i);
121  hMC_PECorr[i] = new TH1D(name,title,1000,0,500);
122 
123  sprintf(name,"MC_trueE_%d",i);
124  hMC_trueE[i] = new TH1D(name,title,1000,0,50);
125 
126  sprintf(name,"MC_cal_%d",i);
127  hMC_cal[i] = new TH1D(name,title,60,0,6);
128  }
129 
130 
131  std::cout << "Filling MC histos ..." <<std::endl;
132  unsigned int nEntries_MC = tMC->GetEntries();
133  for(unsigned int i=0; i<nEntries_MC; i++){
134  tMC->GetEntry(i);
135 
136  if( x_MC>100 && x_MC<200){
137  if(PECorr_MC>0){
138  hMC_PECorr[0] ->Fill(PECorr_MC/path_MC);
139  // hMC_cal[0] ->Fill(mccal*PECorr_MC/path_MC);
140 
141  for(unsigned int ii=0; ii<diblocks; ii++){
142  if(db_MC==(ii+1))
143  hMC_PECorr[db_MC] ->Fill(PECorr_MC/path_MC);
144  // hMC_cal[db_MC] ->Fill(mccal*PECorr_MC/path_MC);
145  }
146  }
147 
148  if(trueE_MC>0){
149  hMC_trueE[0] ->Fill(trueE_MC/path_MC);
150  for(unsigned int ii=0; ii<diblocks; ii++){
151  if(db_MC==(ii+1))
152  hMC_trueE[db_MC] ->Fill(trueE_MC/path_MC);
153  }
154  }
155 
156  }
157 
158  }
159 
160 
161 
162  std::cout << "Making data chain ..." <<std::endl;
163  TChain* tdata = new TChain("muondedxana/fTree");
164  tdata->Add("/nova/ana/users/dmendez/Calibration/Files/FDdata/prod_pcliststop_R16-01-27-prod2calib.a_fd_cosmic_period2_quarter/*.root");
165 
166  double PECorr_data, path_data, trueE_data, x_data;
167  double planes_data, xplanes_data, yplanes_data;
168  int cont_data, nhits_data,hitId_data, db_data;
169 
170  std::cout << "Setting data branches ..." <<std::endl;
171  tdata->SetBranchAddress("PECorr", &PECorr_data);
172  tdata->SetBranchAddress("path" , &path_data);
173  tdata->SetBranchAddress("trueE" , &trueE_data);
174  tdata->SetBranchAddress("x" , &x_data);
175  tdata->SetBranchAddress("diblock",&db_data);
176 
177  std::cout << "Declaring data histos ..." <<std::endl;
178  TH1D *hdata_PECorr[diblocks+1];//adding an extra histogram for calibration over the whole detector
179  TH1D *hdata_cal[diblocks+1];
180 
181 
182  for (unsigned int i=0;i<diblocks+1;i++) {
183  sprintf(title,"DiBlock_%d",i);
184 
185  sprintf(name,"data_PECorr_%d",i);
186  hdata_PECorr[i] = new TH1D(name,title,1000,0,500);
187 
188  sprintf(name,"data_cal_%d",i);
189  hdata_cal[i] = new TH1D(name,title,60,0,6);
190  }
191 
192 
193  std::cout << "Filling data histos ..." <<std::endl;
194 
195  unsigned int nEntries_data = tdata->GetEntries();
196  for(unsigned int i=0; i<nEntries_data; i++){
197  tdata->GetEntry(i);
198 
199  if(x_data>100 && x_data<200)
200  if(PECorr_data>0){
201  hdata_PECorr[0] ->Fill(PECorr_data/path_data);
202  // hdata_cal[0] ->Fill(datacal*PECorr_data/path_data);
203 
204  for(unsigned int ii=0; ii<diblocks; ii++){
205  if(db_data==(ii+1))
206  hdata_PECorr[db_data] ->Fill(PECorr_data/path_data);
207  // hdata_cal[db_data] ->Fill(datacal*PECorr_data/path_data);
208  }
209  }
210 
211  }
212 
213 
214 
215 
216 
217  //mean of distributions with different fits
218  //landau:
219  double mpv_data[diblocks+1], mpv_MC[diblocks+1], mpv_trueE_MC[diblocks+1];
220  double mpvErr_data[diblocks+1], mpvErr_MC[diblocks+1], mpvErr_trueE_MC[diblocks+1];
221  //gaus:
222  double gausMean_data[diblocks+1], gausMean_MC[diblocks+1], gausMean_trueE_MC[diblocks+1];
223  double gausErr_data[diblocks+1], gausErr_MC[diblocks+1], gausErr_trueE_MC[diblocks+1];
224  //mean:
225  double mean_data[diblocks+1], mean_MC[diblocks+1], mean_trueE_MC[diblocks+1];
226  double meanErr_data[diblocks+1], meanErr_MC[diblocks+1], meanErr_trueE_MC[diblocks+1];
227  //median:
228  double Median_data[diblocks+1], Median_MC[diblocks+1], Median_trueE_MC[diblocks+1];
229 
230  //absolute scales with different techniques
231  //landau:
232  double abs_landau_data[diblocks+1], abs_landau_MC[diblocks+1];
233  double absErr_landau_data[diblocks+1], absErr_landau_MC[diblocks+1];
234  //gaus:
235  double abs_gaus_data[diblocks+1], abs_gaus_MC[diblocks+1];
236  double absErr_gaus_data[diblocks+1], absErr_gaus_MC[diblocks+1];
237  //mean:
238  double abs_mean_data[diblocks+1], abs_mean_MC[diblocks+1];
239  double absErr_mean_data[diblocks+1], absErr_mean_MC[diblocks+1];
240  //median:
241  // Need to find technique to calculate error on the median...
242  double abs_median_data[diblocks+1], abs_median_MC[diblocks+1];
243 
244 
245  for(unsigned int i=0; i<diblocks+1; i++){//loop over histos to get means and MEUs
246 
247  //landau
248  TH1D* hdata_PECorr_landau = (TH1D*)hdata_PECorr[i]->Clone("hdata_PECorr_landau");
249  TH1D* hMC_PECorr_landau = (TH1D*)hMC_PECorr[i] ->Clone("hMC_PECorr_landau");
250  TH1D* hMC_trueE_landau = (TH1D*)hMC_trueE[i] ->Clone("hMC_trueE_landau");
251  //gaus
252  TH1D* hdata_PECorr_gaus = (TH1D*)hdata_PECorr[i]->Clone("hdata_PECorr_gaus");
253  TH1D* hMC_PECorr_gaus = (TH1D*)hMC_PECorr[i] ->Clone("hMC_PECorr_gaus");
254  TH1D* hMC_trueE_gaus = (TH1D*)hMC_trueE[i] ->Clone("hMC_trueE_gaus");
255 
256  //set data lines to black:
257  hdata_PECorr[i] -> SetLineColor(1);
258  hdata_PECorr_landau -> SetLineColor(1);
259  hdata_PECorr_gaus -> SetLineColor(1);
260 
261  //set histo range
262  hdata_PECorr[i] -> GetXaxis() -> SetRangeUser(0,100);
263  hMC_PECorr[i] -> GetXaxis() -> SetRangeUser(0,100);
264  hMC_trueE[i] -> GetXaxis() -> SetRangeUser(0,6);
265  hdata_cal[i] -> GetXaxis() -> SetRangeUser(0,6);
266  hMC_cal[i] -> GetXaxis() -> SetRangeUser(0,6);
267 
268 
269  mpv_data[i] = Landau(hdata_PECorr_landau, mpvErr_data[i]);
270  mpv_MC[i] = Landau(hMC_PECorr_landau, mpvErr_MC[i]);
271  mpv_trueE_MC[i] = Landau(hMC_trueE_landau, mpvErr_trueE_MC[i]);
272 
273  gausMean_data[i] = Gaus(hdata_PECorr_gaus, gausErr_data[i], false);
274  gausMean_MC[i] = Gaus(hMC_PECorr_gaus, gausErr_MC[i], false);
275  gausMean_trueE_MC[i] = Gaus(hMC_trueE_gaus, gausErr_trueE_MC[i], true);
276 
277  mean_data[i] = hdata_PECorr[i]->GetMean(1);
278  mean_MC[i] = hMC_PECorr[i]->GetMean(1);
279  mean_trueE_MC[i] = hMC_trueE[i]->GetMean(1);
280  meanErr_data[i] = hdata_PECorr[i]->GetRMS()/(pow(hdata_PECorr[i]->GetEntries(), 0.5));
281  meanErr_MC[i] = hMC_PECorr[i]->GetRMS()/(pow(hMC_PECorr[i]->GetEntries(), 0.5));
282  meanErr_trueE_MC[i] = hMC_trueE[i]->GetRMS()/(pow(hMC_trueE[i]->GetEntries(), 0.5));
283 
284 
285  double pointFive = 0.5;
286  hdata_PECorr[i]->GetQuantiles(1,&Median_data[i],&pointFive);
287  hMC_PECorr[i] ->GetQuantiles(1,&Median_MC[i],&pointFive);
288  hMC_trueE[i] ->GetQuantiles(1,&Median_trueE_MC[i],&pointFive);
289 
290 
291  //absolute scales with different techniques
292  //landau:
293  abs_landau_data[i] = mpv_trueE_MC[i]/mpv_data[i];
294  abs_landau_MC[i] = mpv_trueE_MC[i]/mpv_MC[i];
295  errAoverB(absErr_landau_data[i], mpv_trueE_MC[i], mpv_data[i], mpvErr_trueE_MC[i], mpvErr_data[i]);
296  errAoverB(absErr_landau_MC[i], mpv_trueE_MC[i], mpv_MC[i], mpvErr_trueE_MC[i], mpvErr_MC[i]);
297  //gaus:
298  abs_gaus_data[i] = gausMean_trueE_MC[i]/gausMean_data[i];
299  abs_gaus_MC[i] = gausMean_trueE_MC[i]/gausMean_MC[i];
300  errAoverB(absErr_gaus_data[i], gausMean_trueE_MC[i], gausMean_data[i], gausErr_trueE_MC[i], gausErr_data[i]);
301  errAoverB(absErr_gaus_MC[i], gausMean_trueE_MC[i], gausMean_MC[i], gausErr_trueE_MC[i], gausErr_MC[i]);
302  //mean:
303  abs_mean_data[i] = mean_trueE_MC[i]/mean_data[i];
304  abs_mean_MC[i] = mean_trueE_MC[i]/mean_MC[i];
305  errAoverB(absErr_mean_data[i], mean_trueE_MC[i], mean_data[i], meanErr_trueE_MC[i], meanErr_data[i]);
306  errAoverB(absErr_mean_MC[i], mean_trueE_MC[i], mean_MC[i], meanErr_trueE_MC[i], meanErr_MC[i]);
307  //median:
308  // Need to find technique to calculate error on the median...
309  abs_median_data[i] = Median_trueE_MC[i]/Median_data[i];
310  abs_median_MC[i] = Median_trueE_MC[i]/Median_MC[i];
311 
312 
313  }// end of loop over histos to get means and MEUs
314 
315 
316  ////////////////////////////////////////////////////////////
317  // Printing results on screen //
318  //////////////////////////////////////////////////////////
319  // print MEU values from all techniques along with errors
320  std::cout<< " ----------- Printing MEU values (+/- stat.): --------" <<std::endl;
321 
322  std::cout<< "*** Mean: *** " <<std::endl;
323  std::cout<< " Diblock |" << " data |" << " MC |" << " trueE_MC " <<std::endl;
324  for(unsigned int i=0; i<diblocks+1; i++){
325  std::cout<< i << " " << mean_data[i] << " +/- " << meanErr_data[i] << " | " << mean_MC[i] << " +/- " << meanErr_MC[i] << " | " << mean_trueE_MC[i] << " +/- " << meanErr_trueE_MC[i] <<std::endl;
326  }
327 
328  std::cout<< "*** Landau: *** " <<std::endl;
329  std::cout<< " Diblock |" << " data |" << " MC |" << " trueE_MC " <<std::endl;
330  for(unsigned int i=0; i<diblocks+1; i++){
331  std::cout<< i << " " << mpv_data[i] << " +/- " << mpvErr_data[i] << " | " << mpv_MC[i] << " +/- " << mpvErr_MC[i] << " | " << mpv_trueE_MC[i] << " +/- " << mpvErr_trueE_MC[i] <<std::endl;
332  }
333 
334  std::cout<< "*** Gaus: *** " <<std::endl;
335  std::cout<< " Diblock |" << " data |" << " MC |" << " trueE_MC " <<std::endl;
336  for(unsigned int i=0; i<diblocks+1; i++){
337  std::cout<< i << " " << gausMean_data[i] << " +/- " << gausErr_data[i] << " | " << gausMean_MC[i] << " +/- " << gausErr_MC[i] << " | " << gausMean_trueE_MC[i] << " +/- " << gausErr_trueE_MC[i] <<std::endl;
338  }
339 
340  // Print Absolute scale values with statistical errors
341  std::cout<< "------------- Printing Abs. scale values (+/- stat.): ----------- " <<std::endl;
342 
343  std::cout<< "*** Mean: *** " << std::endl;
344  std::cout<< " Diblock |" << " data |" << " MC " <<std::endl;
345  for(unsigned int i=0; i<diblocks+1; i++){
346  std::cout<< i << " " << abs_mean_data[i] << " +/- " << absErr_mean_data[i] << " | " << abs_mean_MC[i] << " +/- " << absErr_mean_MC[i] << std::endl;
347  }
348  std::cout<< "*** Median: *** " << std::endl;
349  for(unsigned int i=0; i<diblocks+1; i++){
350  std::cout<< i << " " << abs_median_data[i] << " +/- " << absErr_mean_data[i] << " | " << abs_median_MC[i] << " +/- " << absErr_mean_MC[i] << std::endl;
351  }
352  std::cout<< "*** Landau: *** " << std::endl;
353  for(unsigned int i=0; i<diblocks+1; i++){
354  std::cout<< i << " " << abs_landau_data[i] << " +/- " << absErr_landau_data[i] << " | " << abs_landau_MC[i] << " +/- " << absErr_landau_MC[i] << std::endl;
355  }
356  std::cout<< "*** Gaus: *** " <<std::endl;
357  for(unsigned int i=0; i<diblocks+1; i++){
358  std::cout<< i << " " << abs_gaus_data[i] << " +/- " << absErr_gaus_data[i] << " | " << abs_gaus_MC[i] << " +/- " << absErr_gaus_MC[i] << std::endl;
359  }
360 
361 
362  ////////////////////////////////////////////////////////////
363  // Printing results in a .txt file //
364  //////////////////////////////////////////////////////////
365  TString filename = "absCal_DiBlocks_" + detector + period + ".txt";
366  std::ofstream mean(filename);
367  mean << "Includes statistical errors: val1, error1, val2, error1, ..." << std::endl;
368  mean << "Diblock |" << " data |" << " MC |" << " trueE_MC " << " MEU data |" << " MEU MC " <<std::endl;
369  for(unsigned int i=0; i<diblocks+1; i++){
370  mean<< i << " " << mean_data[i] << " " << meanErr_data[i] << " " << mean_MC[i] << " " << meanErr_MC[i] << " " << mean_trueE_MC[i] << " " << meanErr_trueE_MC[i] << " " << abs_mean_data[i] << " " << absErr_mean_data[i] << " " << abs_mean_MC[i] << " " << absErr_mean_MC[i] << std::endl;
371  }
372 
373  mean.close();
374 
375 
376 }
377 
378 
379 
380 
381 //landau
382 double Landau(TH1D* h, double& err){
383  int bin1 = h -> FindFirstBinAbove(h -> GetMaximum()/2);
384  int bin2 = h -> FindLastBinAbove(h -> GetMaximum()/2);
385  double xLow = h -> GetBinCenter(bin1);
386  double xHigh = h -> GetBinCenter(bin2);
387 
388  //fit with landau
389  TF1 *f1 = new TF1("f1","landau",xLow,xHigh);
390  h -> Fit(f1,"","",xLow,xHigh);
391 
392  double mpv = f1->GetParameter(1);
393  err = f1->GetParError(1);
394  std::cout<< "err: " << err <<std::endl;
395  return mpv;
396 }
397 
398 
399 double Gaus(TH1D* h, double& err, bool isTruth){
400 
401  Double_t parameters[3];
402  parameters[0] = h->GetBinContent(h->GetMaximumBin());
403  parameters[1] = h->GetBinCenter(h->GetMaximumBin());
404  parameters[2] = h->GetRMS();
405 
406  // only care about finding peak position, use gaussian truncated about mode
407  //Double_t bound_size = (isTruth) ? 0.2 : 0.5; // narrower for truth, *not* MC reco... don't mess this up.
408 
409  Double_t bound_size;
410  if(isTruth) bound_size = 0.2;
411  else bound_size = 0.5;
412 
413  std::cout<<" ****** bound size: " << bound_size <<std::endl;
414 
415  //Double_t lowerBound = max(0,parameters[1] - bound_size*h->GetRMS());
416  Double_t lowerBound;
417  if(parameters[1] - bound_size*h->GetRMS() > 0)
418  lowerBound = parameters[1] - bound_size*h->GetRMS();
419  else lowerBound = 0;
420 
421  Double_t upperBound = parameters[1] + bound_size*h->GetRMS();
422 
423  TF1* fitfunction = new TF1("fitfunction","[0]*TMath::Gaus(x,[1],[2])",lowerBound,upperBound);
424  fitfunction->SetParameters(parameters);
425 
426  // fitting
427  h->Fit(fitfunction,"RQ");
428 
429  // retrieve values
430  err = fitfunction->GetParError(1);
431  return fitfunction->GetParameter(1);
432 }
433 
434 double median(TH1D* h){
435  return quantiles(h)[2];
436 }
437 
438 std::vector<double> quantiles( TH1D* h ){
439  //=================================
440  double q[5];
441  double probs[5] = {0.025, 0.16, 0.5, 1 - 0.16, 0.975 };
442  h->GetQuantiles( 5, q, probs );
443  std::vector<double> r(5);
444  for (int i=0; i<5; i++)
445  {
446 
447  r[i] = q[i];
448  }
449  return r;
450 }
451 
452 double errAoverB(double& err, double A, double B, double Aerr, double Berr){
453 
454  double C = A/B;
455  //std::cout<< "C: " << C <<std::endl;
456 
457  err = std::sqrt( std::pow((Aerr/B),2) + std::pow((Berr*A)/(B*B),2) );
458  //double err = std::sqrt(std::pow(C,2)) * std::sqrt( std::pow((Berr/B),2) + std::pow((Aerr/A),2));
459 
460  //std::cout<< "err: " << err <<std::endl;
461  return err;
462 
463 }
464 
465 
467  return main();
468 }
469 
470 
const XML_Char * name
Definition: expat.h:151
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
string filename
Definition: shutoffs.py:106
double Gaus(TH1D *h, double &err, bool isTruth)
cout<< t1-> GetEntries()
Definition: plottest35.C:29
const double C
int main()
std::vector< double > quantiles(TH1D *h)
correl_xv GetXaxis() -> SetDecimals()
double median(TH1D *h)
hmean SetLineColor(4)
Float_t f1
double errAoverB(double &err, double A, double B, double Aerr, double Berr)
OStream cout
Definition: OStream.cxx:6
TH1F * h2
Definition: plot.C:45
TH1F * h1
static const double A
Definition: Units.h:82
int absCal_Diblocks()
double Landau(TH1D *h, double &err)
TRandom3 r(0)
double GetHistMax(TH1 *h0, TH1 *h1, TH1 *h2, TH1 *h3)
Definition: absCal.cxx:556