NueCCIncCrossSectionFunctions.h
Go to the documentation of this file.
1 #pragma once
2 
4 #include "RooUnfoldResponse.h"
5 #include "RooUnfoldBayes.h"
6 
7 void PrintChiSq(TString str)
8 {
9  TLatex* prelim = new TLatex(.22, .85, str);
10  prelim->SetTextColor(kBlack);
11  prelim->SetNDC();
12  prelim->SetTextSize(0.05);
13  prelim->Draw();
14 }
15 
16 void Plot2D(TH2F* hist,std::string outName, bool setLog = false)//char* outName)
17 {
18  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str(),"c1");
19  c1->SetLeftMargin(0.20);
20  c1->SetRightMargin(0.20);
21  if(setLog)c1->SetLogz();
22  hist->GetZaxis()->SetTitleOffset(1.10);
23  hist->GetXaxis()->CenterTitle();
24  hist->GetYaxis()->CenterTitle();
25  hist->GetZaxis()->CenterTitle();
26  hist->Draw("COLZ");
27  char out[50];
28  sprintf(out, "Plots/%s", outName.c_str());
29  c1->SaveAs(out);
30  return;
31 }
32 /*
33 std::vector<TH3F*> ApplyFitResults(std::vector<TH3F*> mc,
34  std::vector<TH2F*> wei)
35 {
36  std::vector<TH3F*> results;
37 
38  for(uint i = 0; i < mc.size(); i++){
39  results.push_back((TH3F*)mc[i]->Clone(ana::UniqueName().c_str()));
40  }
41 
42  if(wei.size() != 4){
43  std::cout << "Something failed" << std::endl;
44  return results;
45  }
46 
47  for(int xbin = 0; xbin <= mc[0]->GetXaxis()->GetNbins(); xbin++){
48  for(int ybin = 0; ybin <= mc[0]->GetYaxis()->GetNbins(); ybin++){
49 
50  float nuewei = wei[0]->GetBinContent(xbin,ybin);
51  float nueerr = wei[0]->GetBinError(xbin,ybin);
52  float numuwei = wei[1]->GetBinContent(xbin,ybin);
53  float numuerr = wei[1]->GetBinError(xbin,ybin);
54  float ncwei = wei[2]->GetBinContent(xbin,ybin);
55  float ncerr = wei[2]->GetBinError(xbin,ybin);
56  float correl = wei[3]->GetBinContent(xbin,ybin);
57 
58  for(int zbin = 0; zbin <= mc[0]->GetZaxis()->GetNbins(); zbin++){
59 
60  float n_sig = results[1]->GetBinContent(xbin,ybin,zbin);
61  float n_nuebar = results[2]->GetBinContent(xbin,ybin,zbin);
62  float n_numu = results[3]->GetBinContent(xbin,ybin,zbin);
63  float n_numubar = results[4]->GetBinContent(xbin,ybin,zbin);
64  float n_nc = results[5]->GetBinContent(xbin,ybin,zbin);
65  float n_other = results[6]->GetBinContent(xbin,ybin,zbin);
66  float n_nuenonfid = results[7]->GetBinContent(xbin,ybin,zbin);
67 
68  float err_sig = results[1]->GetBinError(xbin,ybin,zbin);
69  float err_nuebar = results[2]->GetBinError(xbin,ybin,zbin);
70  float err_numu = results[3]->GetBinError(xbin,ybin,zbin);
71  float err_numubar = results[4]->GetBinError(xbin,ybin,zbin);
72  float err_nc = results[5]->GetBinError(xbin,ybin,zbin);
73  float err_other = results[6]->GetBinError(xbin,ybin,zbin);
74  float err_nuenonfid = results[7]->GetBinError(xbin,ybin,zbin);
75 
76 
77 
78  //Propagate Only Systematic Uncertainties
79  float val = nuewei*(n_sig+n_nuebar+n_nuenonfid)+
80  numuwei*(n_numu+n_numubar)+ncwei*(n_nc)+n_other;
81  float err_da = sqrt( pow(n_sig + n_nuebar + n_nuenonfid,2)*
82  pow(nueerr,2) +
83  pow(n_nc+n_numu+n_numubar,2)*pow(ncerr,2) +
84  2*(n_sig+n_nuebar+n_nuenonfid)*
85  (n_nc+n_numu+n_numubar)*correl*nueerr*ncerr);
86  results[0]->SetBinContent(xbin,ybin,zbin,val);
87  results[0]->SetBinError(xbin,ybin,zbin,err_da);
88 
89  float err_background = sqrt(pow(n_nuebar+n_nuenonfid,2)*pow(nueerr,2) +
90  pow(n_nc+n_numu+n_numubar,2)*pow(ncerr,2) +
91  2*(n_nuebar+n_nuenonfid)*
92  (n_nc+n_numu+n_numubar)*
93  correl*nueerr*ncerr);
94 
95 
96  //err = sqrt(pow(nueerr,2)*pow(n_sig,2) + n_sig*pow(nuewei,2));
97  //float err = sqrt(pow(err_da,2)*pow(err_background,2));
98  float err = sqrt(pow(nueerr,2)*pow(n_sig,2));
99  results[1]->SetBinContent(xbin,ybin,zbin,nuewei*n_sig);
100  results[1]->SetBinError(xbin,ybin,zbin,err);
101 
102  //err = sqrt(pow(nueerr,2)*pow(n_nuebar,2) + n_nuebar*pow(nuewei,2));
103  err = sqrt(pow(nueerr,2)*pow(n_nuebar,2));
104  results[2]->SetBinContent(xbin,ybin,zbin,nuewei*n_nuebar);
105  results[2]->SetBinError(xbin,ybin,zbin,err);
106 
107  //err = sqrt(pow(numuerr,2)*pow(n_numu,2) + n_numu*pow(numuwei,2));
108  err = sqrt(pow(numuerr,2)*pow(n_numu,2));
109  results[3]->SetBinContent(xbin,ybin,zbin,numuwei*n_numu);
110  results[3]->SetBinError(xbin,ybin,zbin,err);
111 
112  //err = sqrt(pow(numuerr,2)*pow(n_numubar,2)
113  // + n_numubar*pow(numuwei,2));
114  err = sqrt(pow(numuerr,2)*pow(n_numubar,2));
115  results[4]->SetBinContent(xbin,ybin,zbin,numuwei*n_numubar);
116  results[4]->SetBinError(xbin,ybin,zbin,err);
117 
118  //err = sqrt(pow(ncerr,2)*pow(n_nc,2) + n_nc*pow(ncwei,2));
119  err = sqrt(pow(ncerr,2)*pow(n_nc,2));
120  results[5]->SetBinContent(xbin,ybin,zbin,ncwei*n_nc);
121  results[5]->SetBinError(xbin,ybin,zbin,err);
122 
123  //err = sqrt(pow(nueerr,2)*pow(n_nuenonfid,2) +
124  // n_nuenonfid*pow(nuewei,2));
125  err = sqrt(pow(nueerr,2)*pow(n_nuenonfid,2));
126  results[7]->SetBinContent(xbin,ybin,zbin,nuewei*n_nuenonfid);
127  results[7]->SetBinError(xbin,ybin,zbin,err);
128  }
129  }
130  }
131 
132 
133  //results[0]->Add(mc[0],-1);
134 
135  //for(uint i = 1; i < results.size(); i++)
136  //results[0]->Add(results[i],1);
137 
138  return results;
139 }
140 */
141 
142 TH3F* GetEfficiency(TH3F* num, TH3F* denom)
143 {
144  TH3F* eff = (TH3F*)num->Clone("eff");
145  eff->Divide(denom);
146 
147  return eff;
148 }
149 
150 TH1F* DoUnfolding(TH1F* hMeasured, TH2F* hResponse, int iter)
151 {
152  /*
153  //Pretty Colors
154  //-------------------------------------------------------------------------
155  const Int_t NRGBs = 9;
156  const Int_t NCont = 255;
157 
158  Double_t stops[NRGBs] = { 0.0000, 0.1250, 0.2500, 0.3750, 0.5000, 0.6250, 0.7500, 0.8750, 1.0000};
159  Double_t red[9] = { 251./255., 235./255., 223./255., 214./255., 196./255., 188./255., 157./255., 102./255., 37./255.};
160  Double_t green[9] = { 251./255., 185./255., 132./255., 91./255., 67./255., 37./255., 25./255., 29./255., 37./255.};
161  Double_t blue[9] = { 251./255., 187./255., 137./255., 98./255., 66./255., 45./255., 33./255., 32./255., 37./255.};
162  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
163  gStyle->SetNumberContours(NCont);
164 */
165 
166 
167  RooUnfoldResponse ResponseMatrix(0,0,hResponse, "ResponseMatrix");
168  std::cout << "Using Overflow? : " <<
169  ResponseMatrix.UseOverflowStatus() << std::endl;
170  if(!ResponseMatrix.UseOverflowStatus()){
171  ResponseMatrix.UseOverflow(1);
172  std::cout << "Using Overflow, now? : " <<
173  ResponseMatrix.UseOverflowStatus() << std::endl;
174  }
175 
176  RooUnfoldBayes unfold(&ResponseMatrix, hMeasured, iter);
177  TH1F* hunfolded = (TH1F*)unfold.Hreco();
178  //TMatrixD covMatrix = unfold.Ereco();
179 
180  //Draw Response Matrix
181  TCanvas *c1 = new TCanvas("c1","c1");
182  hResponse->SetMarkerColor(kBlack);
183  hResponse->Draw("COLZ TEXT");
184  gStyle->SetPaintTextFormat("4.f");
185  Simulation();
186  c1->SaveAs("Plots/response_matrix.png");
187  c1->Close();
188  delete c1;
189 
190  /*
191  //Draw covarivance matrix
192  TCanvas *c2 = new TCanvas("c2","c2");
193  covMatrix.Draw("COLZ TEXT");
194  gStyle->SetPaintTextFormat("4.f");
195  Simulation();
196  c2->SaveAs("Plots/covariance_matrix.png");
197  c2->Close();
198  delete c2;
199  */
200 
201  //Draw Unfolded Histogram
202  TCanvas *c3 = new TCanvas("c3","c3");
203  hunfolded->SetTitle("");
204  hunfolded->Draw();
205  Simulation();
206  c3->SaveAs("Plots/unfolded_signal_estimate.png");
207  c3->Close();
208  delete c3;
209 
210  return hunfolded;
211 }
212 
213 void CalculateXSec(TH3F* hNumerator,TH1F* hFlux, TH3F* hEff,
214  std::string dataName,std::string xsecResult,
215  double numTargets)
216 {
217  TFile* outf = new TFile("CrossSection_MCResults.root","update");
218 
219  if(xsecResult == "energy"){
220  TH3F* hXSec = (TH3F*)hNumerator->Clone("hXSec");
221 
222  //3D Efficiency Correction
223  hXSec->Divide(hEff);
224 
225  TH1F* hXSec1D = (TH1F*)hXSec->ProjectionZ("hXSec1D");
226  TH1F* hTrueVar = (TH1F*)hXSec1D->Clone("clone");
227  hFlux->Scale(1e-4);//nu/cm^2 from nu/m^2
228  hXSec1D->Divide(hFlux);
229  hXSec1D->Scale(1./numTargets);
230  hXSec1D->SetTitle("");
231  hXSec1D->GetXaxis()->SetTitle("Neutrino Energy, E_{#nu} (GeV)");
232  hXSec1D->GetYaxis()->SetTitle("#sigma (cm^{2}/nucleon)");
233  hXSec1D->GetXaxis()->CenterTitle();
234  hXSec1D->GetYaxis()->CenterTitle();
235 
236  char name[50];
237  sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
238  "xsec");
239  hXSec1D->SetName(name);
240  hXSec1D->Write();
241  sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
242  "truevar_eff_corrected");
243  hTrueVar->SetName(name);
244  hTrueVar->Write();
245  sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
246  "efficiency");
247  hEff->SetName(name);
248  hEff->Write();
249 
250  outf->Close();
251  }
252  if(xsecResult == "electron"){
253  TH3F* hXSec = (TH3F*)hNumerator->Clone("hXSec");
254  //3D Efficiency Correction
255  hXSec->Divide(hEff);
256 
257  TH1F* hXSec1D = (TH1F*)hXSec->ProjectionY("hXSec1D");
258  TH1F* hTrueVar = (TH1F*)hXSec1D->Clone("clone");
259  hFlux->Scale(1e-4);//nu/cm^2 from nu/m^2
260  double fluxtot = hFlux->Integral();
261  hXSec1D->Scale(1./fluxtot);
262  hXSec1D->Scale(1./numTargets);
263  hXSec1D->Scale(1,"width");
264  hXSec1D->SetTitle("");
265  hXSec1D->GetXaxis()->SetTitle("Electron Energy, E_{e} (GeV)");
266  hXSec1D->GetYaxis()->SetTitle("#sigma (cm^{2}/nucleon GeV)");
267  hXSec1D->GetXaxis()->CenterTitle();
268  hXSec1D->GetYaxis()->CenterTitle();
269 
270  char name[50];
271  sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
272  "xsec");
273  hXSec1D->SetName(name);
274 
275  hXSec1D->Write();
276  sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
277  "truevar_eff_corrected");
278  hTrueVar->SetName(name);
279  hTrueVar->Write();
280  sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
281  "efficiency");
282  hEff->SetName(name);
283  hEff->Write();
284 
285  outf->Close();
286  }
287  if(xsecResult == "cos"){
288  TH3F* hXSec = (TH3F*)hNumerator->Clone("hXSec");
289 
290  //3D Efficiency Correction
291  hXSec->Divide(hEff);
292 
293  TH1F* hXSec1D = (TH1F*)hXSec->ProjectionX("hXSec1D");
294  TH1F* hTrueVar = (TH1F*)hXSec1D->Clone("clone");
295 
296  hFlux->Scale(1e-4);//nu/cm^2 from nu/m^2
297  double fluxtot = hFlux->Integral();
298  hXSec1D->Scale(1./fluxtot);
299  hXSec1D->Scale(1./numTargets);
300  hXSec1D->Scale(1,"width");
301  hXSec1D->SetTitle("");
302  hXSec1D->GetXaxis()->SetTitle("cos #theta_{e}");
303  hXSec1D->GetYaxis()->SetTitle("#sigma (cm^{2}/nucleon cos #theta)");
304  hXSec1D->GetXaxis()->CenterTitle();
305  hXSec1D->GetYaxis()->CenterTitle();
306 
307  char name[50];
308  sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
309  "xsec");
310  hXSec1D->SetName(name);
311 
312  hXSec1D->Write();
313  sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
314  "truevar_eff_corrected");
315  hTrueVar->SetName(name);
316  hTrueVar->Write();
317  sprintf(name, "%s_%s_%s_%s", "mc", dataName.c_str(), xsecResult.c_str(),
318  "efficiency");
319  hEff->SetName(name);
320  hEff->Write();
321 
322  outf->Close();
323  }
324 }
void Simulation()
Definition: tools.h:16
const XML_Char * name
Definition: expat.h:151
void CalculateXSec(TH3F *hNumerator, TH1F *hFlux, TH3F *hEff, std::string dataName, std::string xsecResult, double numTargets)
void Plot2D(TH2F *hist, std::string outName, bool setLog=false)
TH1F * DoUnfolding(TH1F *hMeasured, TH2F *hResponse, int iter)
TFile * outf
Definition: testXsec.C:51
TLatex * prelim
Definition: Xsec_final.C:133
OStream cout
Definition: OStream.cxx:6
int num
Definition: f2_nu.C:119
c1
Definition: demo5.py:24
std::string dataName
void PrintChiSq(TString str)
TH3F * GetEfficiency(TH3F *num, TH3F *denom)
Float_t e
Definition: plot.C:35
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
enum BeamMode string