Xdiff_gwt.C
Go to the documentation of this file.
1 {
2  //extrat differential cross section using BDT
3  double dataPOT=0.260292;//E21
4  //double T=3.24;//(E31);
5  double T=2.68;//(E31);
6  double pot_flux=1.;//(E10);
7 
8  double Ebins[11]={0.};
9  for( int i=0; i<11; ++i )
10  Ebins[i] = i*0.5;
11 
12  //signal
13  TFile *fgenie=new TFile("nue_Xtrue_v280.root");
14  TH1F *hgenie=(TH1F*)fgenie->Get("Xnue");
15  double Xtrue[10]={0.};
16  for( int ibin=1; ibin<hgenie->GetNbinsX()+1; ++ibin )
17  Xtrue[ibin-1] = hgenie->GetBinContent(ibin);
18 
19  //genie-uncertainty on true Xnue
20  TFile *fgwt=new TFile("Etrue_gwt_total_ccnue.root");
21  TH1F *gwt_1p=(TH1F*)fgwt->Get("ccnue_1p_total");
22  TH1F *gwt_2p=(TH1F*)fgwt->Get("ccnue_2p_total");
23  TH1F *gwt_1m=(TH1F*)fgwt->Get("ccnue_1m_total");
24  TH1F *gwt_2m=(TH1F*)fgwt->Get("ccnue_2m_total");
25 
26  vector<double> Xtrue_1p;
27  vector<double> Xtrue_2p;
28  vector<double> Xtrue_1m;
29  vector<double> Xtrue_2m;
30 
31  for( int ibin=1; ibin<gwt_1p->GetNbinsX()+1; ++ibin ){
32  double binctr=gwt_1p->GetBinCenter(ibin);
33  if( binctr<1. || binctr>3. ) continue;
34  Xtrue_1p.push_back((1.+gwt_1p->GetBinContent(ibin))*hgenie->GetBinContent(ibin));
35  Xtrue_2p.push_back((1.+gwt_2p->GetBinContent(ibin))*hgenie->GetBinContent(ibin));
36  Xtrue_1m.push_back((1.-gwt_1m->GetBinContent(ibin))*hgenie->GetBinContent(ibin));
37  Xtrue_2m.push_back((1.-gwt_2m->GetBinContent(ibin))*hgenie->GetBinContent(ibin));
38  //cout<<binctr<<" "<<gwt_1p->GetBinContent(ibin)<<" "<<gwt_2p->GetBinContent(ibin)<<" "<<gwt_1m->GetBinContent(ibin)<<" "<<gwt_2m->GetBinContent(ibin)<<" "<<hgenie->GetBinContent(ibin)<<endl;
39  //cout<<binctr<<" "<<Xtrue_1p[ibin-1]<<" "<<Xtrue_1m[ibin-1]<<" "<<Xtrue_2p[ibin-1]<<" "<<Xtrue_2m[ibin-1]<<endl;
40  }
41 
42  //cout<<"no. of items: "<<Xtrue_1p.size()<<endl;
43  //for( int i=0; i<Xtrue_1p.size(); ++i ){
44  //cout<<i<<" "<<Xtrue_1p[i]<<" "<<Xtrue_1m[i]<<" "<<Xtrue_2p[i]<<" "<<Xtrue_2m[i]<<endl;
45  //}
46 
47  //data - background
48  TFile *fEreco=new TFile("Ereco.root");
49  double Nsub[10]={0.};//background-subtracted data
50  double Nsub_err[10]={0.};
51  double Nsub_dataerr[10]={0.};
52  double Nsub_bkgerr[10]={0.};
53  double Nbkg[10]={0.};
54  double Nbkg_syst[10]={0.};//relative bkg normalization uncertainty
55  double Nsig[10]={0.};
56 
57  TH1F *hd=(TH1F*)fEreco->Get("Ereco_data");
58  TH1F *hbkg=(TH1F*)fEreco->Get("Ereco_bkg");
59  TH1F *hsig=(TH1F*)fEreco->Get("Ereco_sig");
60 
61  double flux_syst[10]={0.,0.,0.101,0.114,0.107,0.121,0,0,0,0};
62  double ehadron_syst[10]={0.,0.,0.099,0.069,0.054,0.015,0,0,0,0};
63  double seleff_syst[10]={0.,0.,0.071,0.071,0.071,0.071,0,0,0,0};
64 
65  //double bkg_sf[10]={0,0,0.892,0.816,0.756,0.834,0,0,0,0};
66  //double bkg_sf_err[10]={0,0,0.028,0.030,0.045,0.149,0,0,0,0};
67  double bkg_sf[10]={0,0,0.95,0.95,0.95,0.95,0,0,0,0};
68  double bkg_sf_err[10]={0,0,0.21,0.21,0.21,0.21,0,0,0,0};
69  double bkg_sf_err_rel[10]={0.};
70  for( int i=0; i<10; ++i ){
71  if( bkg_sf[i]>0. ) bkg_sf_err_rel[i]=bkg_sf_err[i]/bkg_sf[i];
72  }
73 
74  for(int ibin=1; ibin<hd->GetNbinsX()+1; ++ibin ){
75  Nsig[ibin-1] = hsig->GetBinContent(ibin);
76 
77  double Ndata=hd->GetBinContent(ibin);
78  double Nbackground=(hbkg->GetBinContent(ibin))*bkg_sf[ibin-1];//scaling bkg
79 
80  Nsub[ibin-1] = Ndata - Nbackground;
81  Nsub_err[ibin-1] = sqrt((hd->GetBinError(ibin))**2 + (hbkg->GetBinError(ibin))**2);
82  //stat err from data
83  Nsub_dataerr[ibin-1] = (hd->GetBinError(ibin))/Nsub[ibin-1];
84  //stat err from MC
85  Nsub_bkgerr[ibin-1] = (hbkg->GetBinError(ibin))/Nsub[ibin-1];
86 
87  Nbkg[ibin-1] = hbkg->GetBinContent(ibin);
88  if( Nsub[ibin-1]>0. )
89  Nbkg_syst[ibin-1] = (Ndata - Nbackground*(1.-bkg_sf_err_rel[ibin-1]))/Nsub[ibin-1] - 1.;
90  //cout<<Ndata<<" "<<Nbackground<<" "<<bkg_sf_err_rel[ibin-1]<<endl;
91  cout<<"bkg norm uncertainty: "<<ibin<<" "<<Nbkg_syst[ibin-1]<<endl;
92  }
93 
94  //unfolding scale
95  double unfold_scale[10]={0.};
96  TFile *funfold=new TFile("BinByBinUnfold.root");
97  TH1F *hunfold=(TH1F*)funfold->Get("ratio");
98  for( int ibin=1; ibin<hunfold->GetNbinsX()+1; ++ibin ){
99  //unfold_scale[ibin-1]=hunfold->GetBinContent(ibin);
100  unfold_scale[ibin-1]=1.0;
101  }
102 
103  //acc and id efficiency
104  double ideff[10]={0.};
105  TFile *feff=new TFile("SelEff_nueCC.root");
106  TH1F *hid=(TH1F*)feff->Get("toteff");
107  for( int ibin=1; ibin<hid->GetNbinsX()+1; ++ibin ){
108  double binctr=hid->GetBinCenter(ibin);
109  if( binctr>5. ) continue;
110 
111  ideff[ibin-1] = hid->GetBinContent(ibin);
112 
113  if( binctr>1. && binctr<3. )
114  cout<<binctr<<" eff "<<hid->GetBinContent(ibin)<<" "<<hid->GetBinError(ibin)<<endl;
115  }
116 
117  //flux
118  TFile *fflux=new TFile("newflux.root");
119  double flux[10]={0.};
120  double flux_weight[10]={0,0,0.903,0.906,0.915,0.932,0,0,0,0};
121  TH1F *hflux=(TH1F*)fflux->Get("nue_05_fid");
122  for( int ibin=1; ibin<hflux->GetNbinsX()+1; ++ibin ){
123  flux[ibin-1]=hflux->GetBinContent(ibin)/7.84*flux_weight[ibin-1];//7.84=2.8m*2.8m
124  }
125 
126  //systematic error
127  double total_syst[10]={0.};//sqrt(0.2*0.2+0.15*0.15+0.05*0.05+0.1*0.1);
128 
129  double Xdata[10]={0.};
130  double Xbins[10]={0.};
131  double Xbins_err[10]={0.};
132  double Xdata_err[10]={0.};
133  double Xdata_allerr[10]={0.};
134 
135  for( int i=0; i<10; ++i ){
136  double binctr=0.25+i*0.5;
137  if( binctr<1. || binctr>3. ) continue;
138  Xbins[i] = binctr;
139  double Xsec=(Nsub[i]/(flux[i]*dataPOT*T*ideff[i]))*unfold_scale[i];
140  Xdata[i] = Xsec;
141  Xdata_err[i] = Xsec*(Nsub_err[i]/Nsub[i]);
142 
143  //total_syst[i] = sqrt(0.2*0.2+0.05*0.05+0.05*0.05+Nbkg_syst[i]**2);
144  total_syst[i] = sqrt(flux_syst[i]**2+ehadron_syst[i]**2+seleff_syst[i]**2+0.04*0.04+0.03*0.03+0.02*0.02+Nbkg_syst[i]**2);
145  Xdata_allerr[i] = Xsec*sqrt(total_syst[i]**2+(Nsub_err[i]/Nsub[i])**2);
146  //cout<<Ebins[i]<<" - "<<Ebins[i+1]<<": "<<flux[i]<<" "<<Xtrue[i]<<" "<<Xsec<<" +/- "<<Xsec*(Nsub_err[i]/Nsub[i])<<" +/- "<<Xsec*total_syst[i]<<endl;
147 
148  printf("%7.1f $-$ %7.1f & %7.1f & %7.1f $\\\pm$ %7.1f $(\\\delta_{\\\\text{stat}})$ $\\\pm$ %7.1f $(\\\delta_{\\\\text{syst}})$ & %7.1f \\\\\\\hline \n",
149  Ebins[i],Ebins[i+1], flux[i], Xsec, Xdata_err[i], Xsec*total_syst[i], Xtrue[i]);
150 
151  //cout<<"stat: "<<Nsub_dataerr[i]<<" "<<Nsub_bkgerr[i]<<endl;
152  //cout<<"*** "<<Nbkg_syst[i]<<" "<<total_syst[i]<<endl;
153 
154  cout<<"uncertainty: "<<Xdata_err[i]/Xsec<<" "<<total_syst[i]<<endl;
155  cout<<binctr<<" "<<Xsec*10.<<" "<<Xdata_allerr[i]*10.<<" "<<Xdata_err[i]*10.<<endl;
156  cout<<"bkg subtracted data "<<Nsub[i]<<endl;
157  cout<<"syst err "<<total_syst[i]*Xsec<<" "<<Xtrue[i]*10.<<endl;
158  }
159 
160  gStyle->SetOptStat(0);
161 
162  TCanvas *c1=new TCanvas("Xdiff_bdt_gwt","",700,600);
163  c1->cd();
164 
165  TH2F *h2f=new TH2F("h2f","NOvA, 2.6 #times 10^{20} POT",10,1.001,2.999,10,0,5.49);
166  h2f->SetXTitle("E_{#nu_{e}} (GeV)");
167  h2f->SetYTitle("#nu_{e} CC #sigma (10^{-38} cm^{2}/nucleon)");
168  h2f->Draw();
169 
170 
171  TBox *bin1_2sigma=new TBox(1.,Xtrue_2m[0],1.5,Xtrue_2p[0]);
172  TBox *bin2_2sigma=new TBox(1.5,Xtrue_2m[1],2.0,Xtrue_2p[1]);
173  TBox *bin3_2sigma=new TBox(2.0,Xtrue_2m[2],2.5,Xtrue_2p[2]);
174  TBox *bin4_2sigma=new TBox(2.5,Xtrue_2m[3],3.0,Xtrue_2p[3]);
175  bin1_2sigma->SetLineColor(3);
176  bin1_2sigma->SetFillColor(3);
177  bin1_2sigma->Draw("same");
178  bin2_2sigma->SetLineColor(3);
179  bin2_2sigma->SetFillColor(3);
180  bin2_2sigma->Draw("same");
181  bin3_2sigma->SetLineColor(3);
182  bin3_2sigma->SetFillColor(3);
183  bin3_2sigma->Draw("same");
184  bin4_2sigma->SetLineColor(3);
185  bin4_2sigma->SetFillColor(3);
186  bin4_2sigma->Draw("same");
187 
188  TBox *bin1_1sigma=new TBox(1.,Xtrue_1m[0],1.5,Xtrue_1p[0]);
189  TBox *bin2_1sigma=new TBox(1.5,Xtrue_1m[1],2.0,Xtrue_1p[1]);
190  TBox *bin3_1sigma=new TBox(2.0,Xtrue_1m[2],2.5,Xtrue_1p[2]);
191  TBox *bin4_1sigma=new TBox(2.5,Xtrue_1m[3],3.0,Xtrue_1p[3]);
192  bin1_1sigma->SetLineColor(5);
193  bin1_1sigma->SetFillColor(5);
194  bin1_1sigma->Draw("same");
195  bin2_1sigma->SetLineColor(5);
196  bin2_1sigma->SetFillColor(5);
197  bin2_1sigma->Draw("same");
198  bin3_1sigma->SetLineColor(5);
199  bin3_1sigma->SetFillColor(5);
200  bin3_1sigma->Draw("same");
201  bin4_1sigma->SetLineColor(5);
202  bin4_1sigma->SetFillColor(5);
203  bin4_1sigma->Draw("same");
204 
205  cout<<"genie syst: "<<Xtrue_1p[0]<<" "<<Xtrue_1p[1]<<" "<<Xtrue_1p[2]<<" "<<Xtrue_1p[3]<<endl;
206 
207 
208  hgenie->SetLineWidth(3);
209  hgenie->SetLineColor(4);
210  hgenie->SetLineStyle(2);
211  hgenie->Draw("HIST same");
212 
213  GRdata_sys = new TGraphErrors(10,Xbins,Xdata,Xbins_err,Xdata_allerr);
214  GRdata_sys->SetLineWidth(2);
215  GRdata_sys->SetLineColor(2);
216  GRdata_sys->Draw("P z same");
217 
218  gStyle->SetEndErrorSize(4);
219 
220  GRdata = new TGraphErrors(10,Xbins,Xdata,Xbins_err,Xdata_err);
221  GRdata->SetMarkerStyle(8);
222  GRdata->SetMarkerSize(1.0);
223  GRdata->SetLineWidth(2);
224  GRdata->Draw("P same");
225 
226  TLegend *leg = new TLegend(0.6,0.6,0.89,0.89);
227  leg->SetLineColor(0);
228  leg->SetFillColor(0);
229  leg->AddEntry(GRdata,"data","p");
230  leg->AddEntry(GRdata,"stat. uncert.","E");
231  leg->AddEntry(GRdata_sys,"stat. + syst. uncert.","E");
232  leg->AddEntry(hgenie,"GENIE","l");
233  leg->AddEntry(bin1_1sigma,"GENIE #pm 1#sigma","f");
234  leg->AddEntry(bin1_2sigma,"GENIE #pm 2#sigma","f");
235  leg->Draw("same");
236 
237 }
TH1F * hid
Definition: Xdiff_gwt.C:106
TH1F * gwt_1m
Definition: Xdiff_gwt.C:23
double total_syst[10]
Definition: Xdiff_gwt.C:127
TH1F * gwt_2p
Definition: Xdiff_gwt.C:22
T sqrt(T number)
Definition: d0nt_math.hpp:156
double Nbkg_syst[10]
Definition: Xdiff_gwt.C:54
TBox * bin4_1sigma
Definition: Xdiff_gwt.C:191
double Nsub[10]
Definition: Xdiff_gwt.C:49
TH1F * gwt_2m
Definition: Xdiff_gwt.C:24
TBox * bin1_2sigma
Definition: Xdiff_gwt.C:171
double bkg_sf_err_rel[10]
Definition: Xdiff_gwt.C:69
TH1F * hsig
Definition: Xdiff_gwt.C:59
double flux_syst[10]
Definition: Xdiff_gwt.C:61
TBox * bin2_2sigma
Definition: Xdiff_gwt.C:172
TFile * funfold
Definition: Xdiff_gwt.C:96
double Xbins_err[10]
Definition: Xdiff_gwt.C:131
double flux_weight[10]
Definition: Xdiff_gwt.C:120
double bkg_sf_err[10]
Definition: Xdiff_gwt.C:68
double seleff_syst[10]
Definition: Xdiff_gwt.C:63
double pot_flux
Definition: Xdiff_gwt.C:6
TBox * bin1_1sigma
Definition: Xdiff_gwt.C:188
TLegend * leg
Definition: Xdiff_gwt.C:226
TH1F * hflux
Definition: Xdiff_gwt.C:121
double Ebins[11]
Definition: Xdiff_gwt.C:8
TH1F * gwt_1p
Definition: Xdiff_gwt.C:21
double Xdata[10]
Definition: Xdiff_gwt.C:129
double Nbkg[10]
Definition: Xdiff_gwt.C:53
printf("%d Experimental points found\n", nlines)
GRdata_sys
Definition: Xdiff_gwt.C:213
TBox * bin2_1sigma
Definition: Xdiff_gwt.C:189
double unfold_scale[10]
Definition: Xdiff_gwt.C:95
TFile * fEreco
Definition: Xdiff_gwt.C:48
double Xdata_err[10]
Definition: Xdiff_gwt.C:132
TFile * fgwt
Definition: Xdiff_gwt.C:20
TBox * bin3_2sigma
Definition: Xdiff_gwt.C:173
double Nsub_err[10]
Definition: Xdiff_gwt.C:50
vector< double > Xtrue_1p
Definition: Xdiff_gwt.C:26
OStream cout
Definition: OStream.cxx:6
TH2F * h2f
Definition: Xdiff_gwt.C:165
double flux[10]
Definition: Xdiff_gwt.C:119
double bkg_sf[10]
Definition: Xdiff_gwt.C:67
TBox * bin4_2sigma
Definition: Xdiff_gwt.C:174
TFile * fgenie
Definition: Xsec_final.C:56
TFile * feff
Definition: Xdiff_gwt.C:105
TFile * fflux
Definition: Xdiff_gwt.C:118
GRdata
Definition: Xdiff_gwt.C:220
vector< double > Xtrue_2p
Definition: Xdiff_gwt.C:27
double Nsub_bkgerr[10]
Definition: Xdiff_gwt.C:52
double Xbins[10]
Definition: Xdiff_gwt.C:130
TCanvas * c1
Definition: Xdiff_gwt.C:162
TH1F * hunfold
Definition: Xdiff_gwt.C:97
double Xdata_allerr[10]
Definition: Xdiff_gwt.C:133
vector< double > Xtrue_2m
Definition: Xdiff_gwt.C:29
double T
Definition: Xdiff_gwt.C:5
double Nsub_dataerr[10]
Definition: Xdiff_gwt.C:51
vector< double > Xtrue_1m
Definition: Xdiff_gwt.C:28
TH1F * hd
Definition: Xdiff_gwt.C:57
double ehadron_syst[10]
Definition: Xdiff_gwt.C:62
double Nsig[10]
Definition: Xdiff_gwt.C:55
TBox * bin3_1sigma
Definition: Xdiff_gwt.C:190
double ideff[10]
Definition: Xdiff_gwt.C:104
TH1F * hbkg
Definition: Xdiff_gwt.C:58