resolutionplotterbias.C
Go to the documentation of this file.
1 
2 
3 
4 #include <iostream>
5 #include <cmath>
6 #include <iostream>
7 #include <vector>
8 #include <list>
9 #include <sstream>
10 #include <string>
11 #include <sstream>
12 #include <fstream>
13 #include <iomanip>
14 #include "TFile.h"
15 #include "TH2.h"
16 #include "TCanvas.h"
17 #include "TROOT.h"
18 #include "TStyle.h"
19 #include "TLatex.h"
20 #include "TGaxis.h"
21 #include "TLegend.h"
22 #include "TLegendEntry.h"
23 #include "Utilities/rootlogon.C"
24 using namespace std;
25 void NDSimulation()
26 {
27  TLatex* prelim = new TLatex(.9, .95, "NOvA Simulation");
28  prelim->SetTextColor(kGray+1);
29  prelim->SetNDC();
30  prelim->SetTextSize(2/30.);
31  prelim->SetTextAlign(32);
32  TGaxis::SetMaxDigits(3);
33  prelim->Draw();
34 }
35 
36 //****************************Updated Date and string path for saving files **********************************************************///
37  TString outDir2 = "/nova/app/users/acedeno/tag_releaseS17-09-24/NDAna/ncpi0_Systematics/Resolution/Correction";
38  struct FileUp{
39  std::string name;
40  };
41  struct DateUp{
42  std::string date;
43  };
44 
45  const int File =10;
46  const FileUp CurrentFile[File] = {
47  {"/AbsoluteResTruthTotEnWithCorrection"},
48  {"/BiasRecoTotEnWithCorrection"},
49  {"/LogDeltaEnVsPi0RecoWithCorrection"},
50  {"/AbsoluteResRecoTotEnWithCorrection"},
51  {"/BiasTrueTotEnWithCorrection"},
52  {"/LogDeltaEnVsPi0TrueWithCorrection"},
53 
54  };
55 
56  const int Date = 1;
57  const DateUp CurrentDate[Date] = {
58  {"_05-21-18_09_UncorrectedAngle"},
59  };
60 
61 
62 void resolutionplotterbias()
63 {
64  TFile *f = new TFile("resolutionstudybiasMay24-18-Png2_0_8_Cut_BeamAngle_uncorrected.root", "READ");
65  if( !(f) ) return;
66 
67 
68 
69 //////////////Get Spectrums////////////////////////////
70 ////ABSOLUTE RESOLUTION PLOTS with True/Reco///////
71  TH2D* h=(TH2D*)f->Get("selResolutionEn"); //kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
72  TH2D* hh=(TH2D*)f->Get("selSigResolutionEn"); //kIsSignal01L && kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
73  TH2D* hhh=(TH2D*)f->Get("selBkgResolutionEn"); //!kIsSignal01L && kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
74 ////BIAS PLOTS with True/Reco///////
75  TH2D* b=(TH2D*)f->Get("selResolutionEn1Bias"); //kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
76  TH2D* bb=(TH2D*)f->Get("selSigResolutionEn1Bias"); //kIsSignal01L && kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
77  TH2D* bbb=(TH2D*)f->Get("selBkgResolutionEn1Bias"); //!kIsSignal01L && kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
78 // ////ABSOLUTE RESOLUTION PLOTS with True/Reco///////
79 // TH2D* h=(TH2D*)f->Get("selResolutionEnn"); //kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
80 // TH2D* hh=(TH2D*)f->Get("selSigResolutionEnn"); //kIsSignal01L && kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
81 // TH2D* hhh=(TH2D*)f->Get("selBkgResolutionEnn"); //!kIsSignal01L && kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
82 // ////BIAS PLOTS with True/Reco///////
83 // TH2D* b=(TH2D*)f->Get("selResolutionEn2Bias"); //kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
84 // TH2D* bb=(TH2D*)f->Get("selSigResolutionEn2Bias"); //kIsSignal01L && kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
85 // TH2D* bbb=(TH2D*)f->Get("selBkgResolutionEn2Bias"); //!kIsSignal01L && kIsPreSel && kNpCVN1gammaID && kNpCVN2gammaID
86 ////////////////Scale the Spectrum/////////////
87 float datapot = 8.09e20/3.54e21;
88 
89 h->Scale(datapot);
90 hh->Scale(datapot);
91 hhh->Scale(datapot);
92 b->Scale(datapot);
93 bb->Scale(datapot);
94 bbb->Scale(datapot);
95 h->Sumw2();
96 hh->Sumw2();
97 hhh->Sumw2();
98 b->Sumw2();
99 bb->Sumw2();
100 bbb->Sumw2();
101 ///////Get Bins////////////////////////
102 ///////For RESOLUTION PLOTS/////////
103 
104 ///Total PreSelection/////
105  int nbins = h->GetNbinsX();
106  double xmin = h->GetXaxis()->GetXmin();
107  double xmax = h->GetXaxis()->GetXmax();
108 ////Signal//////
109 int nnbins = hh->GetNbinsX();
110 double xxmin = hh->GetXaxis()->GetXmin();
111 double xxmax = hh->GetXaxis()->GetXmax();
112 ////BACKGROUND/////
113 int nnnbins = hhh->GetNbinsX();
114 double xxxmin = hhh->GetXaxis()->GetXmin();
115 double xxxmax = hhh->GetXaxis()->GetXmax();
116 
117 ///////Get Bins////////////////////////
118 ///////For BIAS PLOTS/////////
119 ///Total PreSelection/////
120 int mbins = b->GetNbinsX();
121 double bxmin = b->GetXaxis()->GetXmin();
122 double bxmax = b->GetXaxis()->GetXmax();
123 ////Signal//////
124 int mmbins = bb->GetNbinsX();
125 double bbxmin = bb->GetXaxis()->GetXmin();
126 double bbxmax = bb->GetXaxis()->GetXmax();
127 ////BACKGROUND/////
128 int mmmbins = bbb->GetNbinsX();
129 double bbbxmin = bbb->GetXaxis()->GetXmin();
130 double bbbxmax = bbb->GetXaxis()->GetXmax();
131 
132 /////////Make 1D Histograms///////////
133 
134  TH1D *hrms = new TH1D("hrms" ,"", nbins ,xmin,xmax);
135  TH1D *hhrms = new TH1D("hhrms" ,"", nnbins ,xxmin,xxmax);
136  TH1D *hhhrms = new TH1D("hhhrms" ,"", nnnbins ,xxxmin,xxxmax);
137 
138  TH1D *biashmean = new TH1D("brms" ,"", mbins ,bxmin,bxmax);
139  TH1D *biashhmean = new TH1D("bbrms" ,"", mmbins ,bbxmin,bbxmax);
140  TH1D *biashhhmean = new TH1D("bbbrms" ,"", mmmbins ,bbbxmin,bbbxmax);
141 
142  ////Project a 2D histogram onto a 1D histogram RESOLUTION/Bias
143  ////Total PreSELECTIOn LOOOp////////////////////////////
144 for( int i = 1; i <= nbins; ++i){
145 
146  TH1D *projection = (TH1D*) h->ProjectionY("",i,i);
147  TH1D *projectionb = (TH1D*) b->ProjectionY("",i,i);
148  projection->Sumw2();
149  projectionb->Sumw2();
150  //////bias///////////
151  double mean = projectionb->GetMean(); //for bias vs reco
152  double meanerr = projectionb->GetMeanError(); //for bias vs reco
153  biashmean->SetBinContent(i,mean);
154  biashmean->SetBinError(i,meanerr);
155 
156  ////////resolution/////
157  double rms = projection->GetRMS(); ///for resolution vs true
158  double rmserr = projection->GetRMSError(); ///for resolution vs true
159  hrms->SetBinContent(i,rms);
160  hrms->SetBinError(i,rmserr);
161  }
162 
163 ////Project a 2D histogram onto a 1D histogram RESOLUTION/Bias
164 /// SIGNAL LOOOp////////////////////////////
165 for( int i = 1; i <= nnbins; ++i){
166 
167  TH1D *projection1 = (TH1D*) hh->ProjectionY("",i,i);
168  TH1D *projection1b = (TH1D*) bb->ProjectionY("",i,i);
169  projection1->Sumw2();
170  projection1b->Sumw2();
171  //////bias///////////
172  double hhmean = projection1b->GetMean(); //for bias vs reco
173  double hhmeanerr = projection1b->GetMeanError(); //for bias vs reco
174  biashhmean->SetBinContent(i,hhmean);
175  biashhmean->SetBinError(i,hhmeanerr);
176  ////////resolution/////
177  double rms1 = projection1->GetRMS(); ///for resolution vs true
178  double rmserr1 = projection1->GetRMSError(); //for resolution vs true
179 
180  hhrms->SetBinContent(i,rms1);
181  hhrms->SetBinError(i,rmserr1);
182 
183  }
184 
185 ////Project a 2D histogram onto a 1D histogram RESOLUTION/Bias
186 //// BACKGROUND LOOOp////////////////////////////
187 for( int i = 1; i <= nnnbins; ++i){
188 
189  TH1D *projection2 = (TH1D*) hhh->ProjectionY("",i,i);
190  TH1D *projection2b = (TH1D*) bbb->ProjectionY("",i,i);
191  projection2->Sumw2();
192  projection2b->Sumw2();
193  //////bias///////////
194  double hhhmean = projection2b->GetMean(); //for bias vs reco
195  double hhhmeanerr = projection2b->GetMeanError(); //for bias vs reco
196  biashhhmean->SetBinContent(i,hhhmean);
197  biashhhmean->SetBinError(i,hhhmeanerr);
198  ////////resolution/////
199  double rms2 = projection2->GetRMS(); ///for resolution vs true
200  double rmserr2 = projection2->GetRMSError(); //for resolution vs true
201 
202  hhhrms->SetBinContent(i,rms2);
203  hhhrms->SetBinError(i,rmserr2);
204 
205  }
206 
207 //////Plot Resolution Study///////
208 
209  TCanvas *c = new TCanvas ();
210  c->SetLeftMargin(0.15);
211  c->SetTickx(1);
212  c->SetTicky(1);
213  hrms->SetLineColor(46);
214  hhrms->SetLineColor(31);
215  hhhrms->SetLineColor(38);
216  hrms->Draw("e1");
217  hhrms->Draw("e1 same");
218  hhhrms->Draw("e1 same");
219  hrms->GetYaxis()->SetTitle("Absolute Resolution [Deg.]");
220  //hrms->GetXaxis()->SetTitle("True Total Energy [GeV]");
221  hrms->GetXaxis()->SetTitle("Reco #theta_{beam} [Deg.]");
222  TLegend * leg2 = new TLegend(0.22,0.67,0.32,0.82);
223  //leg2->SetTextAlign(40);
224  leg2->SetTextFont(42);
225  leg2->SetTextSize(0.05);
226  leg2->SetFillStyle(0);
227  leg2->AddEntry( hrms,"Alternate Selection","kFullDotSmall");
228  leg2->AddEntry(hhhrms,"Background","kFullDotSmall");
229  leg2->AddEntry( hhrms,"NC w/Pi0 > 0.3 E","kFullDotSmall");
230  leg2->Draw("same");
231  NDSimulation();
232  c->SaveAs(outDir2+CurrentFile[0].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
233  // c->SaveAs(outDir2+CurrentFile[3].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
234 
235 
236  //////Plot Bias Study///////
237 
238  TCanvas *c1 = new TCanvas ();
239  c1->SetLeftMargin(0.15);
240  c1->SetTickx(1);
241  c1->SetTicky(1);
242  biashmean->SetLineColor(46);
243  biashhmean->SetLineColor(31);
244  biashhhmean->SetLineColor(38);
245  biashmean->Draw("e1");
246  biashhmean->Draw("e1 same");
247  biashhhmean->Draw("e1 same");
248  biashmean->GetYaxis()->SetTitle("Bias [Deg.]");
249  biashmean->GetXaxis()->SetTitle("Reco #theta_{beam} [Deg.]");
250  // biashmean->GetXaxis()->SetTitle("True Total Energy [GeV]");
251  TLegend * leg3 = new TLegend(0.22,0.67,0.32,0.82);
252  //leg2->SetTextAlign(40);
253  leg3->SetTextFont(42);
254  leg3->SetTextSize(0.05);
255  leg3->SetFillStyle(0);
256  leg3->AddEntry( biashmean,"Alternate Selection","kFullDotSmall");
257  leg3->AddEntry(biashhhmean,"Background","kFullDotSmall");
258  leg3->AddEntry( biashhmean,"NC w/Pi0 > 0.3 E","kFullDotSmall");
259  leg3->Draw("same");
260  NDSimulation();
261 c1->SaveAs(outDir2+CurrentFile[1].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
262 //c1->SaveAs(outDir2+CurrentFile[4].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
263 
264 /// Plot 2 D ////
265 TCanvas *cs14 = new TCanvas("cs14","cs14",10,10,855,500);
266  cs14->cd();
267  gPad->SetLogz(1);
268  b->Draw("colz");
269  //b->GetXaxis()->SetTitle("#pi^{0} True Energy [GeV]");
270  b->GetXaxis()->SetTitle("Reco #theta_{beam} [Deg.] ");
271  b->GetYaxis()->SetTitle("(Reco-True) #Delta#theta_{beam} [Deg.]");
272  b->SetTitleOffset(0.9,"Y");
273 
274  NDSimulation();
275  cs14->SaveAs(outDir2+CurrentFile[2].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
276  //cs14->SaveAs(outDir2+CurrentFile[5].name.c_str()+CurrentDate[0].date.c_str()+".pdf");
277 
278 //f->Close();
279 
280 
281 }
282 
name
const XML_Char * name
Definition: expat.h:151
xml_schema::date
::xsd::cxx::tree::date< char, simple_type > date
Definition: Database.h:186
xmax
std::map< std::string, double > xmax
Definition: wrong_sign_uncertainty.C:79
File
const int File
Definition: resolutionplotterbias.C:45
resolutionplotterbias
void resolutionplotterbias()
Definition: resolutionplotterbias.C:62
make_mec_shifts_plots.xmin
float xmin
Definition: make_mec_shifts_plots.py:164
NDSimulation
void NDSimulation()
Definition: resolutionplotterbias.C:25
nbins
const int nbins
Definition: cellShifts.C:15
CurrentFile
const FileUp CurrentFile[File]
Definition: resolutionplotterbias.C:46
make_syst_table_plots.h
h
Definition: make_syst_table_plots.py:240
plot_validation_datamc.c
c
Definition: plot_validation_datamc.py:198
Date
const int Date
Definition: resolutionplotterbias.C:56
prelim
TLatex * prelim
Definition: Xsec_final.C:133
DateUp::date
std::string date
Definition: plot_uncertainty.C:58
outDir2
TString outDir2
Definition: resolutionplotterbias.C:37
b
const hit & b
Definition: hits.cxx:21
demo5.c1
c1
Definition: demo5.py:24
hhh
TH1F * hhh
Definition: AddMC.C:5
CurrentDate
const DateUp CurrentDate[Date]
Definition: resolutionplotterbias.C:57
hh
hh[ixs]
Definition: PlotSingle.C:6
string
enum BeamMode string
Generated on Wed Feb 24 2021 14:12:59 for NOvA by   doxygen 1.8.11