26 #include <TProfile2D.h> 58 double Landau(TH1D*
h,
double& err);
59 double Gaus(TH1D*
h,
double& err,
bool isTruth);
63 double errAoverB(
double& err,
double A,
double B,
double Aerr,
double Berr);
71 double const mccal = 0.0451381;
72 double const datacal = 0.0442841;
80 TChain* tNDMC =
new TChain(
"muondedxana/fTree");
81 tNDMC->Add(
"mc/*.root");
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;
89 tNDMC->SetBranchAddress(
"PECorr", &PECorr_MC);
90 tNDMC->SetBranchAddress(
"path" , &path_MC);
91 tNDMC->SetBranchAddress(
"trueE" , &trueE_MC);
92 tNDMC->SetBranchAddress(
"x" , &x_MC);
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);
101 unsigned int nEntries_MC = tNDMC->GetEntries();
102 for(
unsigned int i=0;
i<nEntries_MC;
i++){
105 if( x_MC>100 && x_MC<200){
107 hMCND_PECorr ->Fill(PECorr_MC/path_MC);
108 hMCND_cal ->Fill(mccal*PECorr_MC/path_MC);
111 hMCND_trueE ->Fill(trueE_MC/path_MC);
120 TChain* tNDdata =
new TChain(
"muondedxana/fTree");
121 tNDdata->Add(
"data/*.root");
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;
128 tNDdata->SetBranchAddress(
"PECorr", &PECorr_data);
129 tNDdata->SetBranchAddress(
"path" , &path_data);
130 tNDdata->SetBranchAddress(
"trueE" , &trueE_data);
131 tNDdata->SetBranchAddress(
"x" , &x_data);
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);
141 unsigned int nEntries_data = tNDdata->GetEntries();
142 for(
unsigned int i=0;
i<nEntries_data;
i++){
143 tNDdata->GetEntry(
i);
148 if(x_data>100 && x_data<200)
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);
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");
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");
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);
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);
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);
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));
200 double Median_dataND;
201 double pointFive = 0.5;
202 hdataND_PECorr->GetQuantiles(1,&Median_dataND,&pointFive);
204 hMCND_PECorr->GetQuantiles(1,&Median_MCND,&pointFive);
205 double Median_trueE_MCND;
206 hMCND_trueE->GetQuantiles(1,&Median_trueE_MCND,&pointFive);
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;
217 std::cout<<
"dataND: " << mpv_dataND <<
" +/- " << mpvErr_dataND <<
" (stat.)" <<
std::endl;
219 std::cout<<
"trueE_MCND: " << mpv_trueE_MCND <<
" +/- " << mpvErr_trueE_MCND <<
" (stat.)" <<
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;
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);
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);
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);
249 double abs_median_dataND = Median_trueE_MCND/Median_dataND;
250 double abs_median_MCND = Median_trueE_MCND/Median_MCND;
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;
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;
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;
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;
271 std::ofstream
values(
"absCal.txt");
274 values<<
" ----------- Printing MEU values: --------" <<
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;
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;
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;
290 values<<
"------------- Printing Abs. scale values: ----------- " <<
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;
295 values<<
"dataND: " << abs_median_dataND <<
" +/- " << absErr_mean_dataND <<
" (stat.)" <<
std::endl;
296 values<<
"MCND: " << abs_median_MCND <<
" +/- " << absErr_mean_MCND <<
" (stat.)" <<
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;
302 values<<
"dataND: " << abs_gaus_dataND <<
" +/- " << absErr_gaus_dataND <<
" (stat.)" <<
std::endl;
303 values<<
"MCND: " << abs_gaus_MCND <<
" +/- " << absErr_gaus_MCND <<
" (stat.)" <<
std::endl;
310 TCanvas *cDataND =
new TCanvas(
"cDataND",
"cDataND",
c1_x,
c1_y);
312 hdataND_PECorr->SetTitle(
"");
313 hdataND_PECorr->GetXaxis()->SetTitle(
"Calibrated detector response/cm");
314 hdataND_PECorr->GetYaxis()->SetTitle(
"Hits in track window");
316 hdataND_PECorr->Draw();
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(
"");
324 hdataND_PECorr_landau->Draw();
328 TCanvas *cDataND_gaus =
new TCanvas(
"cDataND_gaus",
"cDataND_gaus",
c1_x ,
c1_y);
330 hdataND_PECorr_gaus->GetXaxis()->SetTitle(
"Calibrated detector response/cm");
331 hdataND_PECorr_gaus->SetTitle(
"");
333 hdataND_PECorr_gaus->Draw();
337 TCanvas *cMCND =
new TCanvas(
"cMCND",
"cMCND",
c1_x,
c1_y);
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);
344 hMCND_PECorr->Draw();
347 cMCND->Print(
"cMCND.pdf");
350 TCanvas *cMCdataND =
new TCanvas(
"cMCdataND",
"cMCdataND",
c1_x,
c1_y);
352 hMCND_PECorr->SetTitle(
"");
353 hMCND_PECorr->GetXaxis()->SetTitle(
"Corrected response/cm (PECorr/cm)");
354 hMCND_PECorr->GetYaxis()->SetTitle(
"Hits in track window");
356 hMCND_PECorr->SetLineColor(2);
358 hMCND_PECorr->Draw();
359 hdataND_PECorr->Scale(hMCND_PECorr->Integral() / hdataND_PECorr->Integral());
363 hdataND_PECorr->SetTitle(
"");
365 hdataND_PECorr->Draw(
"same");
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();
375 cMCdataND->Print(
"cMCdataND.pdf");
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);
384 hMCND_PECorr_landau->Draw();
388 TCanvas *cMCND_gaus =
new TCanvas(
"cMCND_gaus",
"cMCND_gaus",
c1_x,
c1_y);
390 hMCND_PECorr_gaus->SetTitle(
"");
391 hMCND_PECorr_gaus->GetXaxis()->SetTitle(
"Calibrated detector response/cm");
392 hMCND_PECorr_gaus->SetLineColor(2);
394 hMCND_PECorr_gaus->Draw();
398 TCanvas *cMCND_truth =
new TCanvas(
"cMCND_truth",
"cMCND_truth",
c1_x,
c1_y);
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);
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();
415 cMCND_truth->Print(
"cMCND_truth.pdf");
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");
424 hMCND_cal->SetLineColor(2);
426 hMCND_cal->SetTitle(
"");
428 hdataND_cal->Scale(hMCND_cal->Integral() / hdataND_cal->Integral());
433 hdataND_cal->Draw(
"same");
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();
443 cMCdataND_cal->Print(
"cMCdataND_cal.pdf");
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);
451 hMCND_trueE_landau->Draw();
455 TCanvas *cMCND_truth_gaus =
new TCanvas(
"cMCND_truth_gaus",
"cMCND_truth_gaus",
c1_x,
c1_y);
456 cMCND_truth_gaus ->
cd();
458 hMCND_trueE_gaus->SetTitle(
"");
459 hMCND_trueE_gaus->GetXaxis()->SetTitle(
"MeV/cm");
460 hMCND_trueE_gaus->SetLineColor(2);
462 hMCND_trueE_gaus->Draw();
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);
479 TF1 *
f1 =
new TF1(
"f1",
"landau",xLow,xHigh);
480 h ->
Fit(f1,
"",
"",xLow,xHigh);
482 double mpv = f1->GetParameter(1);
483 err = f1->GetParError(1);
489 double Gaus(TH1D*
h,
double& err,
bool isTruth){
491 Double_t parameters[3];
492 parameters[0] = h->GetBinContent(h->GetMaximumBin());
493 parameters[1] = h->GetBinCenter(h->GetMaximumBin());
494 parameters[2] = h->GetRMS();
500 if(isTruth) bound_size = 0.2;
501 else bound_size = 0.5;
507 if(parameters[1] - bound_size*h->GetRMS() > 0)
508 lowerBound = parameters[1] - bound_size*h->GetRMS();
511 Double_t upperBound = parameters[1] + bound_size*h->GetRMS();
513 TF1* fitfunction =
new TF1(
"fitfunction",
"[0]*TMath::Gaus(x,[1],[2])",lowerBound,upperBound);
514 fitfunction->SetParameters(parameters);
517 h->Fit(fitfunction,
"RQ");
520 err = fitfunction->GetParError(1);
521 return fitfunction->GetParameter(1);
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++)
542 double errAoverB(
double& err,
double A,
double B,
double Aerr,
double Berr){
558 max[0] = h0->GetMaximum();
559 max[1] = h1->GetMaximum();
560 max[2] = h2->GetMaximum();
561 max[3] = h3->GetMaximum();
564 for(
int i=1;
i<4;
i++){
std::vector< double > quantiles(TH1D *h)
double errAoverB(double &err, double A, double B, double Aerr, double Berr)
double Landau(TH1D *h, double &err)
void CenterTitles(TH1 *h)
correl_xv GetXaxis() -> SetDecimals()
double Gaus(TH1D *h, double &err, bool isTruth)
double GetHistMax(TH1 *h0, TH1 *h1, TH1 *h2, TH1 *h3)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)