NueCCIncTemplateFit.h
Go to the documentation of this file.
1 #pragma once
2 
3 #include "CAFAna/Core/SpectrumLoader.h"
4 #include "CAFAna/Core/Spectrum.h"
5 #include "CAFAna/Core/ISyst.h"
6 #include "CAFAna/Core/SystShifts.h"
7 #include "CAFAna/Cuts/SpillCuts.h"
8 #include "CAFAna/Cuts/TruthCuts.h"
9 #include "CAFAna/Systs/XSecSystLists.h"
10 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
11 #include "3FlavorAna/Vars/NueVars.h"
12 #include "3FlavorAna/Vars/NueVarsExtra.h"
13 #include "CAFAna/Vars/GenieWeights.h"
14 #include "CAFAna/Vars/PPFXWeights.h"
15 #include "CAFAna/XSec/GenieMultiverseSyst.h"
16 #include "CAFAna/XSec/FluxMultiverseSyst.h"
17 #include "CAFAna/XSec/Flux.h"
18 #include "CAFAna/XSec/TargetCount.h"
19 #include "CAFAna/Core/ReweightableSpectrum.h"
20 #include "CAFAna/Unfold/UnfoldIterative.h"
21 #include "CAFAna/Analysis/Plots.h"
22 
23 //ROOT Includes
24 #include "Utilities/rootlogon.C"
25 #include "TH1.h"
26 #include "TH2.h"
27 #include "THStack.h"
28 #include "TCanvas.h"
29 #include "TLegend.h"
30 #include "TStyle.h"
31 #include "TFile.h"
32 #include <vector>
33 #include <math.h>
34 #include "TLine.h"
35 #include "TFile.h"
36 #include "TTree.h"
37 #include "TCanvas.h"
38 #include "TF1.h"
39 #include "TH1F.h"
40 #include "TH2F.h"
41 #include "TF2.h"
42 #include "TH2.h"
43 #include "TCutG.h"
44 #include "TMath.h"
45 #include "TCanvas.h"
46 #include "TStyle.h"
47 #include "TRandom.h"
48 #include "TGraph.h"
49 #include "TColor.h"
50 #include "TF3.h"
51 #include "TH3.h"
52 #include "TH3F.h"
53 #include "TMinuit.h"
54 #include "TLegend.h"
55 #include "THStack.h"
56 #include "TMatrix.h"
57 #include "TGraphErrors.h"
58 #include "TGraph.h"
59 #include "TDecompLU.h"
60 #include "TDecompSVD.h"
61 #include "TRandom.h"
62 
63 //C++
64 #include <iostream>
65 #include <fstream>
66 #include <iomanip>
67 #include <vector>
68 #include <cmath>
69 #include <math.h>
70 #include <string>
71 #include <stdint.h>
72 #include <sstream>
73 
74 using namespace ana;
75 
76 const int pidBins = 21;
77 const bool shouldRebin = false;
78 double covMatrix[pidBins][pidBins];
79 double stat_covMatrix[pidBins][pidBins];
80 const float s_b_ratio = 0.40; //ratio of signal/bkgd (in signal region)
81  // that needs to be exceeded
82  // to perform fit in analysis bin
83 const float signalRegion = 0.85;
84 
85 float fDA[pidBins];
86 float fSigLike[pidBins];
87 float fNumuLike[pidBins];
88 float fNC[pidBins];
89 float fOther[pidBins];
90 
91 const std::vector<int> kNueCCColorDef = {
92  632, //kRed, Total MC
93  600, //kBlue Signal
94  616, //kMagenta nuebar CC
95  416-3, //kGreen - 3 NumuCC
96  416+3, //kGreen + 3 Anti NumuCC
97  800, //kOrange NC
98  400-3, //kYellow - 3 Other
99  416-3, //kGreen - 3 Total Background
100 };
101 
102 
103 std::string getCaption2D(TH3F* sample, int xbin, int ybin)
104 {
105  std::stringstream low_X;
106  low_X << std::fixed << std::setprecision(2) <<
107  sample->GetXaxis()->GetBinLowEdge(xbin);
108  std::stringstream hi_X;
109  hi_X << std::fixed << std::setprecision(2) <<
110  sample->GetXaxis()->GetBinUpEdge(xbin);
111 
112  std::stringstream low_Y;
113  low_Y << std::fixed << std::setprecision(2) <<
114  sample->GetYaxis()->GetBinLowEdge(ybin);
115  std::stringstream hi_Y;
116  hi_Y << std::fixed << std::setprecision(2) <<
117  sample->GetYaxis()->GetBinUpEdge(ybin);
118 
119  std::string caption = low_X.str() + " #leq cos #theta_{e} < "+
120  hi_X.str() + " " + low_Y.str() + " #leq E_{e} < " + hi_Y.str();
121  return caption;
122 }
123 
124 void PrintCaption(TString str)
125 {
126 
127  TLatex* prelim = new TLatex(0.23, .80, str);
128  prelim->SetTextColor(kBlack);
129  prelim->SetNDC();
130  prelim->SetTextSize(0.05);
131  prelim->Draw();
132 }
133 
134 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){
135 
136  double chisquare = 0;
137  float data = 0;
138  float signal = 0;
139  float numu = 0;
140  float nc = 0;
141  float other = 0;
142 
143  // Covariance Matrix Inversion First
144  TMatrixF Vij_stat = TMatrix(pidBins,pidBins);
145  TMatrixF Vij = TMatrix(pidBins,pidBins);
146  for(int i =0; i < pidBins; i++){
147  for(int j = 0; j < pidBins; j++){
148  Vij_stat[i][j] = stat_covMatrix[i][j];
149  Vij[i][j] = covMatrix[i][j];
150  }
151  }
152 
153  //Inversion
154  TMatrixD H3 = Vij_stat; //+ Vij;
155  TDecompSVD svd(H3);
156  TMatrixD VijInv = svd.Invert();
157 
158 
159  float delta_i = 0;
160  float delta_j = 0;
161 
162  TMatrixD FitVector = TMatrixD(1,pidBins);
163  TMatrixD FitVector2 = TMatrixD(pidBins,1);
164 
165  //Calculate Chi Sq
166  for(int i =0; i < pidBins; i++){
167  data = fDA[i];
168  signal = fSigLike[i];
169  numu = fNumuLike[i];
170  nc = fNC[i];
171  other = fOther[i];
172 
173  if(data < 30){
174  FitVector[0][i] = 0;
175  FitVector2[i][0] = 0;
176  continue;
177  }
178 
179  FitVector[0][i] = ( data - (par[1]*signal + par[0]*numu +
180  par[2]*nc + other));
181 
182  FitVector2[i][0] = ( data -(par[1]*signal + par[0]*numu +
183  par[2]*nc + other));
184  }
185 
186  TMatrixD Chi = FitVector * VijInv;
187  TMatrixD Chi2 = Chi * FitVector2;
188  float secondterm = Chi2(0,0);
189 
190  chisquare = secondterm;
191  f = chisquare;
192 
193 }
194 
195 void fcn2Var(Int_t &npar, Double_t *gin,
196  Double_t &f, Double_t *par, Int_t iflag){
197 
198  double chisquare = 0;
199  float data = 0;
200  float signal = 0;
201  float numu = 0;
202  float nc = 0;
203  float other = 0;
204 
205  // Covariance Matrix Inversion First
206  TMatrixF Vij_stat = TMatrix(pidBins,pidBins);
207  TMatrixF Vij = TMatrix(pidBins,pidBins);
208  for(int i =0; i < pidBins; i++){
209  for(int j = 0; j < pidBins; j++){
210  Vij_stat[i][j] = stat_covMatrix[i][j];
211  Vij[i][j] = covMatrix[i][j];
212  }
213  }
214 
215  //Inversion
216  TMatrixD H3 = Vij_stat;// + Vij;
217  TDecompSVD svd(H3);
218  TMatrixD VijInv = svd.Invert();
219 
220  TMatrixD FitVector = TMatrixD(1,pidBins);
221  TMatrixD FitVector2 = TMatrixD(pidBins,1);
222 
223  //Calculate Chi Sq
224  for(int i =0; i < pidBins; i++){
225  data = fDA[i];
226  signal = fSigLike[i];
227  numu = fNumuLike[i];
228  nc = fNC[i];
229  other = fOther[i];
230 
231  if(data < 30){
232  FitVector[0][i] = 0;
233  FitVector2[i][0] = 0;
234  continue;
235  }
236 
237  FitVector[0][i] = ( data - (par[1]*signal + par[0]*numu +
238  par[0]*nc + other));
239 
240  FitVector2[i][0] = ( data -(par[1]*signal + par[0]*numu +
241  par[0]*nc + other));
242  }
243 
244  TMatrixD Chi = FitVector * VijInv;
245  TMatrixD Chi2 = Chi * FitVector2;
246  float secondterm = Chi2(0,0);
247 
248  f = secondterm;
249 }
250 
251 void MakeWeightAndErrorPlot(TH2F* hWeights, TH2F* hErr,
252  std::string pidName, std::string dataName,
253  std::string outnameWei, std::string outnameErr,
254  std::string ytitle, std::string xtitle,
255  std::string title,std::string ztitle_wei,
256  std::string ztitle_err, float x_lo = 0.75,
257  float x_hi = 1.0, float y_lo = 1.0,
258  float y_hi = 10.0)
259 {
260  hWeights->GetXaxis()->SetTitle(xtitle.c_str());
261  hWeights->GetYaxis()->SetTitle(ytitle.c_str());
262  hWeights->SetTitle(title.c_str());
263  hWeights->GetXaxis()->SetRangeUser(x_lo,x_hi);
264  hWeights->GetYaxis()->SetRangeUser(y_lo,y_hi);
265  hWeights->GetZaxis()->SetTitle(ztitle_wei.c_str());
266  hWeights->GetXaxis()->CenterTitle();
267  hWeights->GetYaxis()->CenterTitle();
268  hWeights->GetZaxis()->CenterTitle();
269 
270  hErr->GetXaxis()->SetTitle(xtitle.c_str());
271  hErr->GetYaxis()->SetTitle(ytitle.c_str());
272  hErr->SetTitle(title.c_str());
273  hErr->GetXaxis()->SetRangeUser(x_lo,x_hi);
274  hErr->GetYaxis()->SetRangeUser(y_lo,y_hi);
275  hErr->GetZaxis()->SetTitle(ztitle_err.c_str());
276  hErr->GetXaxis()->CenterTitle();
277  hErr->GetYaxis()->CenterTitle();
278  hErr->GetZaxis()->CenterTitle();
279 
280  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
281  hWeights->Draw("COLZ");
282  char out[50];
283  sprintf(out, "%s%s_%s_%s.png","Plots/",pidName.c_str(),
284  dataName.c_str(),outnameWei.c_str());
285  c1->SaveAs(out);
286  delete c1;
287 
288  TCanvas *c2 = new TCanvas(ana::UniqueName().c_str());
289  hErr->Draw("COLZ");
290  sprintf(out, "%s%s_%s_%s.png","Plots/",pidName.c_str(),
291  dataName.c_str(),outnameErr.c_str());
292  c2->SaveAs(out);
293  delete c2;
294 
295 }
296 
297 
298 
299 void PlotFitResults(TH1F* hData,
300  std::vector<TH1F*>hMCnom, std::vector<TH1F*>hMCwei,
301  TGraphAsymmErrors* g, std::string pidName,
302  std::string dataName, std::string caption,
303  int xbin,int ybin)
304 {
305 
306  TH1F* hMC_nom_signal = (TH1F*)hMCnom[1]->Clone(ana::UniqueName().c_str());
307  hMC_nom_signal->Add(hMCnom[2],1);
308  hMC_nom_signal->Add(hMCnom[7],1);
309  TH1F* hMC_nom_bkgd = (TH1F*)hMCnom[3]->Clone(ana::UniqueName().c_str());
310  hMC_nom_bkgd->Add(hMCnom[4],1);
311  hMC_nom_bkgd->Add(hMCnom[5],1);
312  hMC_nom_bkgd->Add(hMCnom[6],1);
313 
314  TH1F* hMC_wgt_signal = (TH1F*)hMCwei[1]->Clone(ana::UniqueName().c_str());
315  hMC_wgt_signal->Add(hMCwei[2],1);
316  hMC_wgt_signal->Add(hMCwei[7],1);
317  TH1F* hMC_wgt_bkgd = (TH1F*)hMCwei[3]->Clone(ana::UniqueName().c_str());
318  hMC_wgt_bkgd->Add(hMCwei[2],1);
319  hMC_wgt_bkgd->Add(hMCwei[4],1);
320  hMC_wgt_bkgd->Add(hMCwei[5],1);
321  hMC_wgt_bkgd->Add(hMCwei[6],1);
322 
323  hMC_nom_signal->SetLineColor(kBlue+3);
324  hMC_nom_bkgd->SetLineColor(kGreen+3);
325 
326  TH1F* hMCWeight = (TH1F*)hMCwei[6]->Clone(ana::UniqueName().c_str());
327  hMCWeight->Add(hMCwei[1]);
328  hMCWeight->Add(hMCwei[2]);
329  hMCWeight->Add(hMCwei[3]);
330  hMCWeight->Add(hMCwei[4]);
331  hMCWeight->Add(hMCwei[5]);
332  hMCWeight->Add(hMCwei[7]);
333  hMCWeight->SetLineColor(kRed);
334 
335  hMCnom[0]->SetLineColor(kMagenta-3);
336  hData->SetMarkerStyle(20);
337  hData->SetTitle(caption.c_str());
338  hMCnom[0]->SetTitle(caption.c_str());
339 
340  if(hData->GetMaximum() > hMCnom[0]->GetMaximum())
341  hMCnom[0]->GetYaxis()->SetRangeUser(0,hData->GetMaximum()*1.3);
342  else(hMCnom[0]->GetYaxis()->SetRangeUser(0,hMCnom[0]->GetMaximum()*1.3));
343 
344  TCanvas* cCompare = new TCanvas(ana::UniqueName().c_str(),"cCompare");
345  cCompare->SetLeftMargin(0.10);
346  cCompare->SetRightMargin(0.10);
347  TPad *pad1 = new TPad(ana::UniqueName().c_str(), "pad1",0,0,1,1);
348  TPad *pad2 = new TPad(ana::UniqueName().c_str(), "pad2",0,0,1,1);
349  pad1->SetFillStyle(0);
350  pad2->SetFillStyle(0);
351  pad1->SetBottomMargin(0.30);
352  pad2->SetTopMargin(1-0.30);
353  pad2->SetGridx();
354  pad2->SetGridy();
355  cCompare->cd();
356  pad1->Draw();
357  pad1->cd();
358  hMCWeight->SetTitle("");
359  hMCnom[0]->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} POT");
360  hMCnom[0]->GetXaxis()->SetLabelColor(kWhite);
361  hMCnom[0]->GetYaxis()->SetLabelSize(0.03);
362  hMCnom[0]->GetYaxis()->SetTitleSize(0.035);
363  hMCnom[0]->SetTitle("");
364  hMCnom[0]->Draw("hist");
365  PrintCaption(caption.c_str());
366  g->Draw("e2 same");
367  hMCnom[0]->Draw("hist same");
368  hMCWeight->Draw("hist same");
369  hData->Draw("same");
370  hMC_wgt_signal->Draw("hist same");
371  hMC_nom_signal->Draw("hist same");
372  hMC_wgt_bkgd->Draw("hist same");
373  hMC_nom_bkgd->Draw("hist same");
374  hData->Draw("same");
375  TLegend* leg = ana::AutoPlaceLegend(0.3,0.3);
376  leg->AddEntry(hData, "Fake Data", "p");
377  leg->AddEntry(hMCWeight, "Scaled Total MC", "l");
378  leg->AddEntry(hMC_wgt_signal, "Scaled Signal", "l");
379  leg->AddEntry(hMC_wgt_bkgd,"Scaled Background", "l");
380  leg->AddEntry(hMCnom[0], "Total MC", "l");
381  leg->AddEntry(hMC_nom_signal, "Signal", "l");
382  leg->AddEntry(hMC_nom_bkgd,"Background", "l");
383  leg->Draw();
384  Simulation();
385  gPad->RedrawAxis();
386  cCompare->cd();
387  pad2->Draw();
388  pad2->cd();
389  TH1F* hratio = (TH1F*)hMCWeight->Clone(ana::UniqueName().c_str());
390  TH1F* hratio2 = (TH1F*)hMCnom[0]->Clone(ana::UniqueName().c_str());
391  hratio->GetYaxis()->SetTitle("MC/Data");
392  hratio->GetXaxis()->SetTitle("Electron Prong CVN Score");
393  hratio->GetXaxis()->CenterTitle();
394  hratio->SetMarkerColor(kBlack);
395  hratio->SetLineColor(kBlack);
396  hratio->SetMarkerStyle(20);
397  hratio->Divide(hData);
398  hratio2->Divide(hData);
399  hratio->GetYaxis()->SetLabelSize(0.03);
400  hratio->GetYaxis()->SetTitleSize(0.035);
401  hratio->GetYaxis()->SetRangeUser(0.4,1.5);
402  hratio->Draw("same");
403  hratio2->Draw("hist same");
404  //hratio2->Draw("hist same");
405  pad2->Draw("same");
406  gPad->RedrawAxis();
407  cCompare->cd();
408  cCompare->Update();
409  char outname[50];
410  sprintf(outname, "%s%s_%s_postfit_ratio_bin_%i_%i.png", "Plots/",
411  pidName.c_str(),
412  dataName.c_str(),xbin,ybin);
413  cCompare->SaveAs(outname);
414  cCompare->Close();
415 }
416 
417 TGraphAsymmErrors* PlotSystErrorBand(TH1F*& nom,
418  std::vector<TH1F*>& ups,
419  std::vector<TH1F*>& dns,
420  int col, int errCol,
421  float headroom, bool newaxis)
422 {
423  if(col == -1){
424  col = kRed;
425  errCol = kRed - 9;
426  }
427  else if(errCol == -1) errCol = col-9; // hopefully a lighter version
428 
429  nom->SetLineColor(col);
430  nom->GetXaxis()->CenterTitle();
431  nom->GetYaxis()->CenterTitle();
432 
433  double yMax = nom->GetBinContent(nom->GetMaximumBin());
434 
435  TGraphAsymmErrors* g = new TGraphAsymmErrors;
436 
437  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
438  const double y = nom->GetBinContent(binIdx);
439  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
440 
441  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
442 
443  double errUp = 0;
444  double errDn = 0;
445 
446  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
447  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
448  double lo = y-dns[systIdx]->GetBinContent(binIdx);
449 
450 
451  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
452 
453  errUp += hi*hi;
454  errDn += lo*lo;
455 
456  // TODO: what happens if they're both high or both low?
457  } // end for systIdx
458 
459  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
460  } // end for i
461 
462  g->SetFillColor(errCol);
463  return g;
464 }
465 
466 
467 void CalculateCovarianceMatrix(TH1F* nom,std::vector<TH1F*> hSysts,
468  int nPseudo,std::string pidName,
469  int xbin,int ybin)
470 {
471  std::vector<TH1F*> vPseudoExp;
472 
473 
474  TH2F* covariance_hist =
475  new TH2F(ana::UniqueName().c_str(),
476  "Covariance Matrix;PID Bin;PID Bin;Covariance",
477  nom->GetXaxis()->GetNbins(),
478  0,nom->GetXaxis()->GetNbins(),
479  nom->GetXaxis()->GetNbins(),
480  0,nom->GetXaxis()->GetNbins());
481 
482 
483  TRandom *r0 = new TRandom();
484  for(int i = 0; i < nPseudo; i++){
485  TH1F* hClone = (TH1F*)nom->Clone(ana::UniqueName().c_str());
486 
487  for(uint j = 0; j < hSysts.size(); j++){
488  double wei = r0->Gaus(0,1);
489 
490  /*std::cout << nom->GetXaxis()->GetNbins() << "\t"
491  << hSysts[j]->GetXaxis()->GetNbins() <<"\t"
492  << nom->GetXaxis()->GetBinLowEdge(1) << "\t"
493  << hSysts[j]->GetBinLowEdge(1) << "\t"
494  << nom->GetXaxis()->GetBinLowEdge(21) << "\t"
495  << hSysts[j]->GetBinLowEdge(21) << std::endl;*/
496  hClone->Add(hSysts[j],wei);
497  }
498 
499  vPseudoExp.push_back(hClone);
500  }
501 
502  //Mij = < (Val_i)*(Val_j)> - <Val_i> <Val_j>
503  //covMatrix has to be filled
504  for(int j = 1; j <= hSysts[0]->GetXaxis()->GetNbins(); j++){
505  for(int k = 1; k <= hSysts[0]->GetXaxis()->GetNbins(); k++){
506  double multiple = 0;
507  double x_j = 0;
508  double x_k = 0;
509  double counter = 0;
510  for(int i = 0; i < nPseudo; i++){
511  double val_j = vPseudoExp[i]->GetBinContent(j);
512  double val_k = vPseudoExp[i]->GetBinContent(k);
513  multiple += val_j * val_k;
514  x_j += val_j;
515  x_k += val_k;
516  counter++;
517  }
518  double covariance = multiple/counter - ((x_j/counter)*(x_k/counter));
519  covMatrix[j-1][k-1] = covariance;
520  covariance_hist->SetBinContent(j,k,covariance);
521  }
522  }
523 
524  char out[50];
525  sprintf(out, "Plots/covariance_%s_%i_%i.png", pidName.c_str(),xbin,ybin);
526  TCanvas *c1 = new TCanvas("c1","c1");
527  c1->SetLogz(1);
528  c1->SetLeftMargin(0.15);
529  c1->SetRightMargin(0.15);
530  covariance_hist->GetXaxis()->CenterTitle();
531  covariance_hist->GetYaxis()->CenterTitle();
532  covariance_hist->GetZaxis()->CenterTitle();
533  covariance_hist->Draw("COLZ");
534  c1->SaveAs(out);
535  delete c1;
536  return;
537 }
538 
539 //ToDo: Look through this function and make sure it is calculating
540 //things properly
541 void CalculateTotalCovariance(TH1F* nominal,std::vector<TH1F*> hSysts,
542  std::string pidName,
543  int xbin,int ybin)
544 {
545  TH2F* covariance_hist =
546  new TH2F(ana::UniqueName().c_str(),
547  "Covariance Matrix;Systematic Number;Systematic Number;Covariance",
548  (int)hSysts.size(),
549  0,(int)hSysts.size(),
550  (int)hSysts.size(),
551  0,(int)hSysts.size());
552 
553  TH2F* correlation_hist =
554  new TH2F(ana::UniqueName().c_str(),
555  "Correlation Matrix;;;Correlation",
556  (int)hSysts.size(),
557  0,(int)hSysts.size(),
558  (int)hSysts.size(),
559  0,(int)hSysts.size());
560 
561 
562  std::vector<double> vAverages;
563  for(int j = 0; j < (int)hSysts.size(); j++){
564  double avg = 0;
565  for(int i = 1; i <= nominal->GetXaxis()->GetNbins(); i++){
566  avg += hSysts[j]->GetBinContent(i);
567  }
568  avg /= nominal->GetXaxis()->GetNbins();
569  vAverages.push_back(avg);
570  }
571 
572  const char *labels[6] = {"Cal.Shp.","Ckv","Calib","Light","XSec","PPFX"};
573 
574 
575  //Sample Covariance = 1/N Sum( (x_i - <x>) (y_i - <y>))
576  for(int j = 0; j < (int) hSysts.size(); j++){
577  for(int k = 0; k < (int)hSysts.size(); k++){
578  double summed_term = 0;
579  double nterms = 0;
580  for(int i = 1; i <= nominal->GetXaxis()->GetNbins(); i++){
581  summed_term += (hSysts[j]->GetBinContent(i) - vAverages[j]) *
582  (hSysts[k]->GetBinContent(i) - vAverages[k]);
583  nterms++;
584  }
585  double covariance = 1./nterms * summed_term;
586  covariance_hist->SetBinContent(j+1,k+1,covariance);
587  }
588  }
589 
590  for(int j = 1; j <= covariance_hist->GetXaxis()->GetNbins(); j++){
591  for(int k = 1; k <= covariance_hist->GetXaxis()->GetNbins(); k++){
592  double covariance = covariance_hist->GetBinContent(j,k);
593  double correlation =
594  covariance/(sqrt(covariance_hist->GetBinContent(j,j)) *
595  sqrt(covariance_hist->GetBinContent(k,k)));
596  correlation_hist->Fill(labels[j-1],labels[k-1],correlation);
597  }
598  }
599 
600 
601  char out[50];
602  sprintf(out, "Plots/covariance_systs_%s_%i_%i.png", pidName.c_str(),xbin,ybin);
603  TCanvas *c1 = new TCanvas("c1","c1");
604  c1->SetLeftMargin(0.15);
605  c1->SetRightMargin(0.15);
606  covariance_hist->GetXaxis()->CenterTitle();
607  covariance_hist->GetYaxis()->CenterTitle();
608  covariance_hist->GetZaxis()->CenterTitle();
609  covariance_hist->Draw("COLZ");
610  c1->SaveAs(out);
611  delete c1;
612 
613 
614  correlation_hist->LabelsDeflate("X");
615  correlation_hist->LabelsDeflate("Y");
616  correlation_hist->LabelsOption("v");
617  sprintf(out, "Plots/correlation_systs_%s_%i_%i.png",
618  pidName.c_str(),xbin,ybin);
619  TCanvas *c2 = new TCanvas("c2","c2");
620  c2->SetLeftMargin(0.15);
621  c2->SetRightMargin(0.15);
622  correlation_hist->GetXaxis()->CenterTitle();
623  correlation_hist->GetYaxis()->CenterTitle();
624  correlation_hist->GetZaxis()->CenterTitle();
625  correlation_hist->Draw("COLZ text");
626  c2->SaveAs(out);
627  delete c2;
628  return;
629 
630 }
631 /*
632 void CalculateTotalCovariance(TH1F* nominal,std::vector<TH1F*> hSysts,
633  std::string pidName,
634  int xbin,int ybin,std::string checkForAllSysts)
635 {
636  //This function is designed to try to
637  //identify the relationship between different systematic uncertainties
638  //by calculate the covariance between the systematics after integrating
639  //over all of the PID bins
640 
641  if(checkForAllSysts != "tot") return;
642 
643  const char *labels[6] = {"Cal.Shp.","Ckv","Calib","Light","XSec","PPFX"};
644 
645  std::vector<TH1F*> vPseudoExp_CalibShape;
646  std::vector<TH1F*> vPseudoExp_Ckv;
647  std::vector<TH1F*> vPseudoExp_Calib;
648  std::vector<TH1F*> vPseudoExp_Light;
649  std::vector<TH1F*> vPseudoExp_XSec;
650  std::vector<TH1F*> vPseudoExp_PPFX;
651 
652 
653  TH2F* covariance_hist =
654  new TH2F(ana::UniqueName().c_str(),
655  "Covariance Matrix;PID Bin;PID Bin;Covariance",
656  nom->GetXaxis()->GetNbins(),
657  0,nom->GetXaxis()->GetNbins(),
658  nom->GetXaxis()->GetNbins(),
659  0,nom->GetXaxis()->GetNbins());
660 
661 
662  TRandom *r0 = new TRandom();
663  for(int i = 0; i < nPseudo; i++){
664  TH1F* hClone_CalibShape = (TH1F*)nom->Clone(ana::UniqueName().c_str());
665  TH1F* hClone_Ckv = (TH1F*)nom->Clone(ana::UniqueName().c_str());
666  TH1F* hClone_Calib = (TH1F*)nom->Clone(ana::UniqueName().c_str());
667  TH1F* hClone_Light = (TH1F*)nom->Clone(ana::UniqueName().c_str());
668  TH1F* hClone_XSec = (TH1F*)nom->Clone(ana::UniqueName().c_str());
669  TH1F* hClone_PPFX = (TH1F*)nom->Clone(ana::UniqueName().c_str());
670 
671  double wei_calibshape = r0->Gaus(0,1);
672  double wei_ckv = r0->Gaus(0,1);
673  double wei_calib = r0->Gaus(0,1);
674  double wei_light = r0->Gaus(0,1);
675  double wei_xsec = r0->Gaus(0,1);
676  double wei_ppfx = r0->Gaus(0,1);
677 
678  hClone_CalibShape->Add(hSysts[0],wei_calibshape);
679  hClone_Ckv->Add(hSysts[0],wei_ckv);
680  hClone_Calib->Add(hSysts[0],wei_calib);
681  hClone_Light->Add(hSysts[0],wei_light);
682  hClone_XSec->Add(hSysts[0],wei_xsec);
683  hClone_PPFX->Add(hSysts[0],wei_ppfx);
684 
685  vPseudoExp_CalibShape.push_back(hClone_CalibShape);
686  vPseudoExp_Ckv.push_back(hClone_Ckv);
687  vPseudoExp_Calib.push_back(hClone_Calib);
688  vPseudoExp_Light.push_back(hClone_Light);
689  vPseudoExp_XSec.push_back(hClone_XSec);
690  vPseudoExp_PPFX.push_back(hClone_PPFX);
691  };
692 
693  //Mij = < (Val_i)*(Val_j)> - <Val_i> <Val_j>
694  //covMatrix has to be filled
695  for(int j = 1; j <= nom->GetXaxis()->GetNbins(); j++){
696  for(int k = 1; k <= nom->GetXaxis()->GetNbins(); k++){
697  double multiple = 0;
698  double x_j = 0;
699  double x_k = 0;
700  double counter = 0;
701  for(int i = 0; i < nPseudo; i++){
702  double val_j = vPseudoExp[i]->GetBinContent(j);
703  double val_k = vPseudoExp[i]->GetBinContent(k);
704  multiple += val_j * val_k;
705  x_j += val_j;
706  x_k += val_k;
707  counter++;
708  }
709  double covariance = multiple/counter - ((x_j/counter)*(x_k/counter));
710  covMatrix[j-1][k-1] = covariance;
711  covariance_hist->SetBinContent(j,k,covariance);
712  }
713  }
714 
715 
716 
717 
718 }
719 */
720 
721 std::vector<TH3F*> BinByBinTemplateFit(TH3F* data,
722  std::vector<TH3F*> templates,
723  std::vector<TH3F*> systs_hists,
724  std::vector<TH3F*> systs_hists_up,
725  std::vector<TH3F*> systs_hists_down,
726  std::vector<TH3F*> analysis_templates,
727  std::string pidName, std::string varName,
728  std::string dataName)
729 {
730  std::vector<TH3F*> badentry;
731  if(systs_hists_up.size() != systs_hists_down.size()){ std::cerr << "Unequal up and down shifts. Exiting" << std::endl;
732  }
733 
734  std::vector<TH3F*> hResults;
735  for(uint i = 0; i < analysis_templates.size(); i++){
736  hResults.push_back((TH3F*)analysis_templates[i]->
737  Clone(ana::UniqueName().c_str()));
738  }
739 
740 
741  TH2F* hSigWeights = (TH2F*)data->Project3D("yx");
742  hSigWeights->SetName(ana::UniqueName().c_str());
743  TH2F* hSigError = (TH2F*)data->Project3D("yx");
744  hSigError->SetName(ana::UniqueName().c_str());
745  TH2F* hNumuWeights = (TH2F*)data->Project3D("yx");
746  hNumuWeights->SetName(ana::UniqueName().c_str());
747  TH2F* hNumuError = (TH2F*)data->Project3D("yx");
748  hNumuError->SetName(ana::UniqueName().c_str());
749  TH2F* hNCWeights = (TH2F*)data->Project3D("yx");
750  hNCWeights->SetName(ana::UniqueName().c_str());
751  TH2F* hNCError = (TH2F*)data->Project3D("yx");
752  hNCError->SetName(ana::UniqueName().c_str());
753 
754  for(int xbin = 1; xbin <= data->GetXaxis()->GetNbins(); xbin++){
755  for(int ybin = 1; ybin <= data->GetYaxis()->GetNbins(); ybin++){
756 
757  std::string caption = getCaption2D(templates[0],xbin,ybin);
758 
759 
760  //Variables to hold fit results
761  std::pair<double,double> nue_wei = std::make_pair (0,0);
762  std::pair<double,double> numu_wei = std::make_pair (0,0);
763  std::pair<double,double> nc_wei = std::make_pair (0,0);
764  Double_t covariance[3][3];
765 
766  //Project Data and MC to PID axis for template fit
767  TH1F* hDataZ = (TH1F*)data->ProjectionZ(ana::UniqueName().c_str(),
768  xbin,xbin,ybin,ybin);
769 
770  std::vector<TH1F*> hMCZ;
771  for(uint numChns = 0; numChns < templates.size(); numChns++){
772  hMCZ.push_back((TH1F*)templates[numChns]->
773  ProjectionZ(ana::UniqueName().c_str(),
774  xbin,xbin,ybin,ybin));
775  hMCZ[numChns]->SetLineColor(kNueCCColorDef[numChns]);
776  }
777 
778  //Make Sure Fit should be performed in this bin (xbin,ybin)
779  int pid_bin = hDataZ->GetXaxis()->FindBin(signalRegion);
780  float signal_like = hMCZ[1]->Integral(pid_bin,-1) +
781  hMCZ[2]->Integral(pid_bin,-1) + hMCZ[7]->Integral(pid_bin,-1);
782  float bkgd_amount = hMCZ[0]->Integral(pid_bin,-1) - signal_like;
783 
784  if(signal_like > 100 &&
785  signal_like/bkgd_amount > s_b_ratio){
786 
787  //Calculate Uncertainties due to each systematic sample
788  std::vector<TH1F*> hSystsZ;
789 
790  //For One-Sided Systematics
791  for(int systNum = 0; systNum < (int) systs_hists.size(); systNum++){
792  TH1F* hHolder =
793  (TH1F*)systs_hists[systNum]->ProjectionZ(ana::UniqueName().c_str(),
794  xbin,xbin,ybin,ybin);
795  //Loop Over Bins, Calculate the absolute distance of systematic from
796  //Nominal MC
797  for(int bin = 1; bin <= hHolder->GetXaxis()->GetNbins(); bin++){
798  double err = hHolder->GetBinContent(bin);
799  double cv = hMCZ[0]->GetBinContent(bin);
800  hHolder->SetBinContent(bin,fabs(err-cv));
801  }
802  hSystsZ.push_back(hHolder);
803  }
804 
805  //For Two-Side Systematics
806  //Let's take a conservative approach, so furthest value away from
807  //nominal is taken as the uncertainty for that bin. Should deal with
808  //highly assymmetric systematics better.
809  for(int systNum = 0; systNum < (int)systs_hists_up.size(); systNum++){
810  TH1F* hHolder = (TH1F*)systs_hists_up[systNum]->
811  ProjectionZ(ana::UniqueName().c_str(),xbin,xbin,ybin,ybin);
812  TH1F* hHolder_dwn = (TH1F*)systs_hists_down[systNum]->
813  ProjectionZ(ana::UniqueName().c_str(),xbin,xbin,ybin,ybin);
814 
815  //Loop over bins, convert two-sided uncertainty into a single value
816  //which is the largest distance of each shift from nominal mc
817  for(int bin = 1; bin <= hHolder->GetXaxis()->GetNbins(); bin++){
818  double hi = hHolder->GetBinContent(bin);
819  double lo = hHolder_dwn->GetBinContent(bin);
820  double cv = hMCZ[0]->GetBinContent(bin);
821  double value_hi = fabs(hi - cv);
822  double value_lo = fabs(lo - cv);
823  double value = 0.5*(fabs(cv-hi) + fabs(cv-lo));
824  (value_hi >= value_lo)? value = value_hi : value = value_lo;
825  //std::cout << bin << "\t" << value << std::endl;
826  hHolder->SetBinContent(bin,value);
827  }
828  hSystsZ.push_back(hHolder);
829  }
830 
831  CalculateTotalCovariance(hMCZ[0],hSystsZ,pidName,xbin,ybin);
832  //Calculate Systematic Covariance Matrix
833  //Throw nPseudo experiments
834  int nPseudo = 1000;
835  CalculateCovarianceMatrix(hMCZ[0],hSystsZ,nPseudo,pidName,xbin,ybin);
836 
837  //CheckHere
838 
839  //Calculate Statistical Covariance Matrix
840  //Diagonal with Data statistics for each pid bin
841  for(int jj = 1; jj <= hSystsZ[0]->GetXaxis()->GetNbins(); jj++){
842  for(int kk = 1; kk <= hSystsZ[0]->GetXaxis()->GetNbins(); kk++){
843  if(jj==kk){
844  stat_covMatrix[jj-1][kk-1] =
845  hDataZ->GetBinError(jj)*hDataZ->GetBinError(jj)
846  + hMCZ[0]->GetBinError(jj)*hMCZ[0]->GetBinError(jj);
847  }
848  else{
849  stat_covMatrix[jj-1][kk-1] = 0.0;
850  }
851  }
852  }
853 
854 
855  //Arrays are most(?) efficienct way to feed values to Minuit
856  //Fill Arrays for MinuitFit
857  for(int x = 1; x <= hDataZ->GetXaxis()->GetNbins(); x++){
858  float data = hDataZ->GetBinContent(x);
859  float signal = hMCZ[1]->GetBinContent(x);
860  float nuebar = hMCZ[2]->GetBinContent(x);
861  float numu = hMCZ[3]->GetBinContent(x);
862  float numubar = hMCZ[4]->GetBinContent(x);
863  float nc = hMCZ[5]->GetBinContent(x);
864  float other = hMCZ[6]->GetBinContent(x);
865  float nuenonfiducial = hMCZ[7]->GetBinContent(x);
866 
867  std::cout << "Data: " << data << "\tMC: " <<
868  signal + nuebar + numu+numubar+nc+other+nuenonfiducial << std::endl;
869 
870  fDA[x-1] = data;
871  fSigLike[x-1] = signal + nuebar + nuenonfiducial;
872  fNumuLike[x-1] = numu + numubar;
873  fNC[x-1] = nc;
874  fOther[x-1] = other;
875 
876  //std::cout << fDA[x-1] << "\t" << fSigLike[x-1] << "\t"
877  //<< fNumuLike[x-1] << "\t" << fNC[x-1] << "\t" << fOther[x-1]
878  // << std::endl;
879 
880  }
881 
882 
883  //Now perform the fit
884  //Fit 1: Get Good Starting Fit Values, Quickly
885  //Three parameter fit
886  TMinuit *gMinuit = new TMinuit(3);
887  gMinuit->SetPrintLevel(-1);
888  gMinuit->SetFCN(fcn);
889  Double_t arglist[4];
890  Int_t ierflg = 0;
891  arglist[0] = 1;
892  gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
893 
894  double vstart[3] = {0.8,1.2,0.9};
895  double verror[3] = {0,0,0};
896  double vstep[3] = {0.01,0.01,0.01};
897 
898  //Define Fit Parameters
899  gMinuit->mnparm(0, "Numu-Like Scaling",
900  vstart[0], vstep[0], 0, 0, ierflg);
901  gMinuit->mnparm(1, "Nue-Like Scaling",
902  vstart[1], vstep[1], 0, 0, ierflg);
903  gMinuit->mnparm(2, "NC Scaling", vstart[2], vstep[2], 0, 0, ierflg);
904 
905  arglist[0] = 0; //number of iterations 0 == no limit
906  arglist[1] = 1; //1 == standard minimization strategy
907 
908  //SIMPLEX should provide a good starting point for further fits
909  gMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
910 
911  //Fit 2: Improved Fit + Better Eror Estimation
912  arglist[1] = 2; //2 == try to improve minimum (slower)
913  gMinuit->mnexcm("SET STR", arglist,1,ierflg);
914 
915  gMinuit->mnexcm("SIMPLEX", arglist, 3, ierflg);
916  gMinuit->mnexcm("MIGRAD", arglist, 5, ierflg);
917  gMinuit->mnexcm("SEEk", arglist, 4, ierflg);
918  gMinuit->mnexcm("MINOS", arglist, 4, ierflg);
919 
920  //Grab Final Results
921  double fParNewStart;
922  double fParNewErr;
923  gMinuit->GetParameter(0,fParNewStart,fParNewErr);
924  vstart[0] = fParNewStart;
925  verror[0] = fParNewErr;
926  gMinuit->GetParameter(1,fParNewStart,fParNewErr);
927  vstart[1] = fParNewStart;
928  verror[1] = fParNewErr;
929  gMinuit->GetParameter(2,fParNewStart,fParNewErr);
930  vstart[2] = fParNewStart;
931  verror[2] = fParNewErr;
932 
933  //Before Ending Fits, Make Sure There Was No Issue With the Fit
934  bool troubleFitting = false;
935  Double_t fmin = 0,fedm = 0,ferrdef = 0;
936  Int_t npari =0, nparx = 0, istat = 0;
937  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
938 
939  Double_t corr_mat_test[3][3];
940  gMinuit->mnemat(*corr_mat_test,3);
941 
942  //Turn this test error matix into a correlation matrix
943  double corr_01 = corr_mat_test[0][1]/(corr_mat_test[0][0] *
944  corr_mat_test[1][1]);
945  double corr_02= corr_mat_test[0][2]/(corr_mat_test[0][0] *
946  corr_mat_test[2][2]);
947  double corr_12= corr_mat_test[1][2]/(corr_mat_test[1][1] *
948  corr_mat_test[2][2]);
949  if(istat == 2 || vstart[0] < 0 || vstart[1] < 0 || vstart[2] < 0 ||
950  fabs(corr_01) > 0.96 || fabs(corr_02) > 0.96 || fabs(corr_12) > 0.96)
951  troubleFitting=true;
952 
953  if(troubleFitting){
954  //There Was trouble with the original fits,
955  //most likely this was due to extremely high correlation (approx 1)
956  //between two of the fitting parameters
957  //Switch to a fit that adjusts all backgrounds
958  TMinuit *jMinuit = new TMinuit(2);
959  jMinuit->SetPrintLevel(-1);
960  jMinuit->SetFCN(fcn2Var);
961  arglist[0] = 1;
962  jMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
963 
964  jMinuit->mnparm(0, "Background", vstart[0],
965  vstep[0], 0, 0, ierflg);
966  jMinuit->mnparm(1, "Nue CC Scaling", vstart[1],
967  vstep[1], 0, 0, ierflg);
968 
969  arglist[0] = 0; //number of iterations 0 == no limit
970  arglist[1] = 1;
971 
972  jMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
973  jMinuit->mnexcm("MIGRAD", arglist, 3, ierflg);
974  jMinuit->mnexcm("SEEk", arglist, 3, ierflg);
975  jMinuit->mnexcm("MINOS", arglist, 3, ierflg);
976 
977 
978  jMinuit->GetParameter(0,fParNewStart,fParNewErr);
979  vstart[0] = fParNewStart;
980  verror[0] = fParNewErr;
981  jMinuit->GetParameter(0,fParNewStart,fParNewErr);
982  vstart[2] = fParNewStart;
983  verror[2] = fParNewErr;
984  jMinuit->GetParameter(1,fParNewStart,fParNewErr);
985  vstart[1] = fParNewStart;
986  verror[1] = fParNewErr;
987  //Get the covariance matrix
988  Double_t emat[2][2];
989  jMinuit->mnemat(*emat,2);
990 
991  nue_wei = std::make_pair (vstart[1],verror[1]);
992  numu_wei = std::make_pair (vstart[0],verror[0]);
993  nc_wei = std::make_pair (vstart[2],verror[2]);
994  covariance[0][0] = emat[0][0];
995  covariance[0][1] = emat[0][1];
996  covariance[0][2] = 0;
997  covariance[1][0] = emat[1][0];
998  covariance[1][1] = emat[1][1];
999  covariance[1][2] = 0;
1000  covariance[2][0] = emat[0][0];
1001  covariance[2][1] = emat[0][1];
1002  covariance[2][2] = emat[0][0];
1003  }
1004  else{
1005  gMinuit->GetParameter(0,fParNewStart,fParNewErr);
1006  vstart[0] = fParNewStart;
1007  verror[0] = fParNewErr;
1008  gMinuit->GetParameter(2,fParNewStart,fParNewErr);
1009  vstart[2] = fParNewStart;
1010  verror[2] = fParNewErr;
1011  gMinuit->GetParameter(1,fParNewStart,fParNewErr);
1012  vstart[1] = fParNewStart;
1013  verror[1] = fParNewErr;
1014 
1015  //Get the covariance matrix
1016  Double_t emat[3][3];
1017  gMinuit->mnemat(*emat,3);
1018 
1019  nue_wei = std::make_pair (vstart[1],verror[1]);
1020  numu_wei = std::make_pair (vstart[0],verror[0]);
1021  nc_wei = std::make_pair (vstart[2],verror[2]);
1022  covariance[0][0] = emat[0][0];
1023  covariance[0][1] = emat[0][1];
1024  covariance[0][2] = emat[0][2];
1025  covariance[1][0] = emat[1][0];
1026  covariance[1][1] = emat[1][1];
1027  covariance[1][2] = emat[1][2];
1028  covariance[2][0] = emat[2][0];
1029  covariance[2][1] = emat[2][1];
1030  covariance[2][2] = emat[2][2];
1031  }
1032  }
1033  else{
1034  std::cout << "******** Skipping Fit For Bin "
1035  << caption << "\t\n\n\n\n" << std::endl;
1036  std::cout << "******** Skipping Fit For Bin "
1037  << caption << "\t" << signal_like << "\t"
1038  << signal_like/bkgd_amount << "\n\n\n\n" << std::endl;
1039 
1040  nue_wei = std::make_pair (1,0);
1041  numu_wei = std::make_pair (1,0);
1042  nc_wei = std::make_pair (1,0);
1043  covariance[0][0] = 0;
1044  covariance[0][1] = 0;
1045  covariance[0][2] = 0;
1046  covariance[1][0] = 0;
1047  covariance[1][1] = 0;
1048  covariance[1][2] = 0;
1049  covariance[2][0] = 0;
1050  covariance[2][1] = 0;
1051  covariance[2][2] = 0;
1052  }
1053 
1054  std::cout << "*********Fit for bin " << caption << " : " << std::endl;
1055  std::cout << "NumuCC Scaling Factor= " << numu_wei.first << " +- " <<
1056  numu_wei.second << std::endl;
1057  std::cout << "NueCC Scaling Factor= " << nue_wei.first << " +- " <<
1058  nue_wei.second << std::endl;
1059  std::cout << "NC Scaling Factor= " << nc_wei.first << " +- " <<
1060  nc_wei.second << std::endl;
1061 
1062  //Push Weights and Uncertainties to 2D Histograms
1063  hSigWeights->SetBinContent(xbin,ybin,nue_wei.first);
1064  hSigError->SetBinContent(xbin,ybin,nue_wei.second);
1065  hNumuWeights->SetBinContent(xbin,ybin,numu_wei.first);
1066  hNumuError->SetBinContent(xbin,ybin,numu_wei.second);
1067  hNCWeights->SetBinContent(xbin,ybin,nc_wei.first);
1068  hNCError->SetBinContent(xbin,ybin,nc_wei.second);
1069 
1070 
1071  //Projections For Systematic Error Band Drawing
1072  std::vector<TH1F*> hSysts;
1073  for(int z = 0; z < (int)systs_hists.size(); z++){
1074  hSysts.push_back((TH1F*)systs_hists[z]->
1075  ProjectionZ(ana::UniqueName().c_str(),
1076  xbin,xbin,ybin,ybin));
1077  }
1078  std::vector<TH1F*> hSystsUp;
1079  for(int z = 0; z < (int)systs_hists_up.size();z++){
1080  hSystsUp.push_back((TH1F*)systs_hists_up[z]->
1081  ProjectionZ(ana::UniqueName().c_str(),
1082  xbin,xbin,
1083  ybin,ybin));
1084  }
1085  std::vector<TH1F*> hSystsDown;
1086  for(int z = 0; z < (int)systs_hists_down.size();z++){
1087  hSystsDown.push_back((TH1F*)systs_hists_down[z]->
1088  ProjectionZ(ana::UniqueName().c_str(),
1089  xbin,xbin,
1090  ybin,ybin));
1091  }
1092 
1093  std::vector<TH1F*> shifts_up;
1094  std::vector<TH1F*> shifts_down;
1095 
1096  for(int z = 0; z < (int)hSysts.size();z++){
1097  shifts_up.push_back(hSysts[z]);
1098  shifts_down.push_back((TH1F*)hMCZ[0]->Clone(ana::UniqueName().c_str()));
1099  }
1100 
1101  for(int z = 0; z < (int)hSystsUp.size(); z++){
1102  shifts_up.push_back(hSystsUp[z]);
1103  shifts_down.push_back(hSystsDown[z]);
1104  }
1105 
1106  //Get Systematic Error Band
1107  TGraphAsymmErrors* g2 = PlotSystErrorBand(hMCZ[0],
1108  shifts_up,
1109  shifts_down,
1110  -1,-1,1.3,false);
1111 
1112  //Make Copies of Nominal MC template predictions
1113  std::vector<TH1F*> hMCnom;
1114  for(int z = 0; z < (int)hMCZ.size(); z++){
1115  hMCnom.push_back((TH1F*)hMCZ[z]->Clone(ana::UniqueName().c_str()));
1116  }
1117 
1118  //Scale MC templates by fitted weights
1119  hMCZ[1]->Scale(nue_wei.first); //Signal
1120  hMCZ[2]->Scale(nue_wei.first); //Nuebar
1121  hMCZ[3]->Scale(numu_wei.first); //Numu
1122  hMCZ[4]->Scale(numu_wei.first); //Numubar
1123  hMCZ[5]->Scale(nc_wei.first); //NC
1124  hMCZ[6]->Scale(1); //Other
1125  hMCZ[7]->Scale(nue_wei.first); //Non-fiducial Nue
1126 
1127 
1128  //Plot Template Before after Fit
1129  PlotFitResults(hDataZ,hMCnom,hMCZ,g2, pidName,dataName,caption,xbin,ybin);
1130 
1131  //Apply Fit Results to Analysis Histogram
1132  for(int numChn = 0; numChn < (int)hResults.size(); numChn++){
1133  for(int zbin = 1; zbin <= hResults[0]->GetZaxis()->GetNbins();zbin++){
1134  float nsig = hResults[1]->GetBinContent(xbin,ybin,zbin);
1135  float nnuebar = hResults[2]->GetBinContent(xbin,ybin,zbin);
1136  float nnumu = hResults[3]->GetBinContent(xbin,ybin,zbin);
1137  float nnumubar = hResults[4]->GetBinContent(xbin,ybin,zbin);
1138  float nnc = hResults[5]->GetBinContent(xbin,ybin,zbin);
1139  float nother = hResults[6]->GetBinContent(xbin,ybin,zbin);
1140  float nnf = hResults[7]->GetBinContent(xbin,ybin,zbin);
1141  if(numChn == 0){
1142  //Add all Channels to total
1143  float binCont = nsig * nue_wei.first;
1144  binCont += nnuebar * nue_wei.first;
1145  binCont += nnumu * numu_wei.first;
1146  binCont += nnumubar * numu_wei.first;
1147  binCont += nnc * nc_wei.first;
1148  binCont += nother;
1149  binCont += nnf * nue_wei.first;
1150  hResults[0]->SetBinContent(xbin,ybin,zbin,binCont);
1151 
1152  float binErr2 =
1153  pow(covariance[0][0],2)*pow((nnumu+nnumubar),2) +
1154  pow(covariance[1][1],2)* pow((nsig+nnuebar+nnf),2) +
1155  pow(covariance[2][2],2)* pow((nnc),2) +
1156  2 * covariance[0][1] * (nnumu+nnumubar)*(nsig+nnuebar+nnf) +
1157  2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) +
1158  2 * covariance[1][2] * (nsig + nnuebar + nnf) * (nnc) +
1159  nsig * pow(nue_wei.first,2) + nnuebar * pow(nue_wei.first,2) +
1160  nnf * pow(nue_wei.first,2) + nnumu * pow(numu_wei.first,2) +
1161  nnumubar * pow(numu_wei.first,2) + nnc * pow(nc_wei.first,2) +
1162  nother;
1163 
1164  if(covariance[0][1] == 0 && covariance[0][2] == 0 &&
1165  covariance[1][2] == 0 && nue_wei.second == 0 &&
1166  numu_wei.second == 0 && nc_wei.second == 0){
1167  binErr2 = nsig + nnuebar + nnf + nnumu + nnumubar +
1168  nnc + nother;
1169  }
1170 
1171 
1172  //if(binErr2 < 0) binErr2 = fabs(binErr2);
1173 
1174  hResults[0]->SetBinError(xbin,ybin,zbin,sqrt(binErr2));
1175 
1176  std::cout << zbin << "\t" << binCont << "\t" <<
1177  sqrt(binErr2) << std::endl;
1178  if(std::isnan(sqrt((binErr2)))){
1179  std::cout << pow(covariance[0][0],2)*pow((nnumu+nnumubar),2) <<
1180  "\t" << pow(covariance[1][1],2)* pow((nsig+nnuebar+nnf),2) <<
1181  "\t" << pow(covariance[2][2],2)* pow((nnc),2) << "\t" <<
1182  2 * covariance[0][1] * (nnumu+nnumubar)*(nsig+nnuebar+nnf) <<
1183  "\t" << 2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) <<
1184  "\t" << 2 * covariance[1][2] * (nsig + nnuebar + nnf) * (nnc) <<
1185  "\t" << nsig * pow(nue_wei.first,2) << "\t" <<
1186  nnuebar * pow(nue_wei.first,2) << "\t" <<
1187  nnf * pow(nue_wei.first,2) << "\t" <<
1188  nnumu * pow(numu_wei.first,2) << "\t" <<
1189  nnumubar * pow(numu_wei.first,2) << "\t" <<
1190  nnc * pow(nc_wei.first,2) << "\t" << nother << std::endl;
1191 
1192  }
1193  }
1194  if(numChn == 1){
1195  //Signal Estimate
1196  float binCont = nsig * nue_wei.first;
1197 
1198  //Just uncertainty on Signal estimate
1199  //Signal + Background (taken into account in TMinuit internally
1200  //when calculating the internal error matrix)
1201  //Template Fit Uncertainties
1202  //float binErr = sqrt( pow(nue_wei.second,2)*pow(nsig,2) +
1203  // nsig * pow(nue_wei.first,2));
1204 
1205  //Need to add wrong sign estimate uncertainty
1206  //Lets just add this in quadrature to the
1207  //signal estimate + background estimate uncertainties
1208  //We add uncertainties on the signal template estimate and
1209  //statistical uncertainties on the WS component
1210  float binErr = sqrt( pow(nue_wei.second,2)*pow(nsig,2) +
1211  nsig * pow(nue_wei.first,2) +
1212  pow(nue_wei.second,2)*pow(nnuebar,2) +
1213  nnuebar * pow(nue_wei.first,2));
1214 
1215  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1216  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1217  }
1218  if(numChn == 2){
1219  float binCont = nnuebar * nue_wei.first;
1220  float binErr = sqrt( pow(nue_wei.second,2)*pow(nnuebar,2) +
1221  nnuebar * pow(nue_wei.first,2));
1222  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1223  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1224  }
1225  if(numChn == 3){
1226  float binCont = nnumu * numu_wei.first;
1227  float binErr = sqrt( pow(numu_wei.second,2)*pow(nnumu,2) +
1228  nnumu * pow(numu_wei.first,2));
1229  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1230  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1231  }
1232  if(numChn == 4){
1233  float binCont = nnumubar * numu_wei.first;
1234  float binErr = sqrt( pow(numu_wei.second,2)*pow(nnumubar,2) +
1235  nnumubar * pow(numu_wei.first,2));
1236  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1237  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1238  }
1239  if(numChn == 5){
1240  float binCont = nnc * nc_wei.first;
1241  float binErr = sqrt( pow(nc_wei.second,2)*pow(nnc,2) +
1242  nnc * pow(nc_wei.first,2));
1243  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1244  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1245  }
1246  if(numChn == 7){
1247  float binCont = nnf * nue_wei.first;
1248  float binErr = sqrt( pow(nue_wei.second,2)*pow(nnf,2) +
1249  nnf * pow(nue_wei.first,2));
1250  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1251  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1252  }
1253  }//loop over zbins
1254  } // loop over interaction channels
1255  }//loop over ybins
1256  }//loop over xbins
1257 
1258  //Draw Weights and Relative Uncertainty
1259  MakeWeightAndErrorPlot(hSigWeights,hSigError, pidName,dataName,
1260  "nue_weights", "nue_errors",
1261  "Reconstructed Electron Energy, E_{e} (GeV)",
1262  "Reconstructed cos #theta_{e}",
1263  "Signal-Like Template Fit Results",
1264  "Normalization", "Uncertainty");
1265  MakeWeightAndErrorPlot(hNumuWeights,hNumuError, pidName,dataName,
1266  "numu_weights", "numu_errors",
1267  "Reconstructed Electron Energy, E_{e} (GeV)",
1268  "Reconstructed cos #theta_{e}",
1269  "Numu-Like Template Fit Results",
1270  "Normalization", "Uncertainty");
1271  MakeWeightAndErrorPlot(hNCWeights,hNCError, pidName,dataName,
1272  "nc_weights", "nc_errors",
1273  "Reconstructed Electron Energy, E_{e} (GeV)",
1274  "Reconstructed cos #theta_{e}",
1275  "NC Template Fit Results",
1276  "Normalization", "Uncertainty");
1277 
1278  return hResults;
1279 }
1280 
1281 //Get Fitted Template Results Full
1282 //Necessary to calculate Chi^{2} Value correctly w/ all degrees of freedom
1283 std::vector<TH3F*> BinByBinTemplateFit_TemplateResults(TH3F* data,
1284  std::vector<TH3F*> templates,
1285  std::vector<TH3F*> systs_hists,
1286  std::vector<TH3F*> systs_hists_up,
1287  std::vector<TH3F*> systs_hists_down,
1288  std::vector<TH3F*> analysis_templates,
1289  std::string pidName, std::string varName,
1290  std::string dataName)
1291 {
1292  std::vector<TH3F*> badentry;
1293  if(systs_hists_up.size() != systs_hists_down.size()){
1294  std::cerr << "Unequal up and down shifts. Exiting" << std::endl;
1295  }
1296 
1297  std::vector<TH3F*> hResults;
1298  for(uint i = 0; i < templates.size(); i++){
1299  hResults.push_back((TH3F*)templates[i]->
1300  Clone(ana::UniqueName().c_str()));
1301  }
1302 
1303  TH2F* hSigWeights = (TH2F*)data->Project3D("yx");
1304  hSigWeights->SetName(ana::UniqueName().c_str());
1305  TH2F* hSigError = (TH2F*)data->Project3D("yx");
1306  hSigError->SetName(ana::UniqueName().c_str());
1307  TH2F* hNumuWeights = (TH2F*)data->Project3D("yx");
1308  hNumuWeights->SetName(ana::UniqueName().c_str());
1309  TH2F* hNumuError = (TH2F*)data->Project3D("yx");
1310  hNumuError->SetName(ana::UniqueName().c_str());
1311  TH2F* hNCWeights = (TH2F*)data->Project3D("yx");
1312  hNCWeights->SetName(ana::UniqueName().c_str());
1313  TH2F* hNCError = (TH2F*)data->Project3D("yx");
1314  hNCError->SetName(ana::UniqueName().c_str());
1315 
1316  for(int xbin = 1; xbin <= data->GetXaxis()->GetNbins(); xbin++){
1317  for(int ybin = 1; ybin <= data->GetYaxis()->GetNbins(); ybin++){
1318 
1319  std::string caption = getCaption2D(templates[0],xbin,ybin);
1320 
1321 
1322  //Variables to hold fit results
1323  std::pair<double,double> nue_wei = std::make_pair (0,0);
1324  std::pair<double,double> numu_wei = std::make_pair (0,0);
1325  std::pair<double,double> nc_wei = std::make_pair (0,0);
1326  Double_t covariance[3][3];
1327 
1328  //Project Data and MC to PID axis for template fit
1329  TH1F* hDataZ = (TH1F*)data->ProjectionZ(ana::UniqueName().c_str(),
1330  xbin,xbin,ybin,ybin);
1331 
1332  std::vector<TH1F*> hMCZ;
1333  for(uint numChns = 0; numChns < templates.size(); numChns++){
1334  hMCZ.push_back((TH1F*)templates[numChns]->
1335  ProjectionZ(ana::UniqueName().c_str(),
1336  xbin,xbin,ybin,ybin));
1337  hMCZ[numChns]->SetLineColor(kNueCCColorDef[numChns]);
1338  }
1339 
1340  //Make Sure Fit should be performed in this bin (xbin,ybin)
1341  int pid_bin = hDataZ->GetXaxis()->FindBin(signalRegion);
1342  float signal_like = hMCZ[1]->Integral(pid_bin,-1) +
1343  hMCZ[2]->Integral(pid_bin,-1) + hMCZ[7]->Integral(pid_bin,-1);
1344  float bkgd_amount = hMCZ[0]->Integral(pid_bin,-1) - signal_like;
1345 
1346  if(signal_like > 100 &&
1347  signal_like/bkgd_amount > s_b_ratio){
1348 
1349  //Calculate Uncertainties due to each systematic sample
1350  std::vector<TH1F*> hSystsZ;
1351 
1352  //For One-Sided Systematics
1353  for(int systNum = 0; systNum < (int) systs_hists.size(); systNum++){
1354  TH1F* hHolder =
1355  (TH1F*)systs_hists[systNum]->ProjectionZ(ana::UniqueName().c_str(),
1356  xbin,xbin,ybin,ybin);
1357  //Loop Over Bins, Calculate the absolute distance of systematic from
1358  //Nominal MC
1359  for(int bin = 1; bin <= hHolder->GetXaxis()->GetNbins(); bin++){
1360  double err = hHolder->GetBinContent(bin);
1361  double cv = hMCZ[0]->GetBinContent(bin);
1362  hHolder->SetBinContent(bin,fabs(err-cv));
1363  }
1364  hSystsZ.push_back(hHolder);
1365  }
1366 
1367  /*
1368  //For Two-Side Systematics
1369  for(int systNum = 0; systNum < (int)systs_hists_up.size(); systNum++){
1370  TH1F* hHolder = (TH1F*)systs_hists_up[systNum]->
1371  ProjectionZ(ana::UniqueName().c_str(),xbin,xbin,ybin,ybin);
1372  TH1F* hHolder_dwn = (TH1F*)systs_hists_down[systNum]->
1373  ProjectionZ(ana::UniqueName().c_str(),xbin,xbin,ybin,ybin);
1374 
1375  //Loop over bins, convert two-sided uncertainty into a single value
1376  //which is the average distance of each shift from nominal mc
1377  for(int bin = 1; bin <= hHolder->GetXaxis()->GetNbins(); bin++){
1378  double hi = hHolder->GetBinContent(bin);
1379  double lo = hHolder_dwn->GetBinContent(bin);
1380  double cv = hMCZ[0]->GetBinContent(bin);
1381  double value = 0.5*(fabs(cv-hi) + fabs(cv-lo));
1382  hHolder->SetBinContent(bin,value);
1383  }
1384  hSystsZ.push_back(hHolder);
1385  }
1386  */
1387 
1388  //For Two-Side Systematics
1389  //Let's take a conservative approach, so furthest value away from
1390  //nominal is taken as the uncertainty for that bin. Should deal with
1391  //highly assymmetric systematics better.
1392  for(int systNum = 0; systNum < (int)systs_hists_up.size(); systNum++){
1393  TH1F* hHolder = (TH1F*)systs_hists_up[systNum]->
1394  ProjectionZ(ana::UniqueName().c_str(),xbin,xbin,ybin,ybin);
1395  TH1F* hHolder_dwn = (TH1F*)systs_hists_down[systNum]->
1396  ProjectionZ(ana::UniqueName().c_str(),xbin,xbin,ybin,ybin);
1397 
1398  //Loop over bins, convert two-sided uncertainty into a single value
1399  //which is the largest distance of each shift from nominal mc
1400  for(int bin = 1; bin <= hHolder->GetXaxis()->GetNbins(); bin++){
1401  double hi = hHolder->GetBinContent(bin);
1402  double lo = hHolder_dwn->GetBinContent(bin);
1403  double cv = hMCZ[0]->GetBinContent(bin);
1404  double value_hi = fabs(hi - cv);
1405  double value_lo = fabs(lo - cv);
1406  double value = 0.5*(fabs(cv-hi) + fabs(cv-lo));
1407  (value_hi >= value_lo)? value = value_hi : value = value_lo;
1408  hHolder->SetBinContent(bin,value);
1409  }
1410  hSystsZ.push_back(hHolder);
1411  }
1412 
1413 
1414  CalculateTotalCovariance(hMCZ[0],hSystsZ,pidName,xbin,ybin);
1415  //Calculate Systematic Covariance Matrix
1416  //Throw nPseudo experiments
1417  int nPseudo = 1000;
1418  CalculateCovarianceMatrix(hMCZ[0],hSystsZ,nPseudo,pidName,xbin,ybin);
1419 
1420  //Calculate Statistical Covariance Matrix
1421  //Diagonal with Data statistics for each pid bin
1422  for(int jj = 1; jj <= hSystsZ[0]->GetXaxis()->GetNbins(); jj++){
1423  for(int kk = 1; kk <= hSystsZ[0]->GetXaxis()->GetNbins(); kk++){
1424  if(jj==kk){
1425  stat_covMatrix[jj-1][kk-1] =
1426  hDataZ->GetBinError(jj)*hDataZ->GetBinError(jj)
1427  +hMCZ[0]->GetBinError(jj)*hMCZ[0]->GetBinError(jj);
1428  }
1429  else{
1430  stat_covMatrix[jj-1][kk-1] = 0.0;
1431  }
1432  }
1433  }
1434 
1435 
1436  //Arrays are most(?) efficienct way to feed values to Minuit
1437  //Fill Arrays for MinuitFit
1438  for(int x = 1; x <= hDataZ->GetXaxis()->GetNbins(); x++){
1439  float data = hDataZ->GetBinContent(x);
1440  float signal = hMCZ[1]->GetBinContent(x);
1441  float nuebar = hMCZ[2]->GetBinContent(x);
1442  float numu = hMCZ[3]->GetBinContent(x);
1443  float numubar = hMCZ[4]->GetBinContent(x);
1444  float nc = hMCZ[5]->GetBinContent(x);
1445  float other = hMCZ[6]->GetBinContent(x);
1446  float nuenonfiducial = hMCZ[7]->GetBinContent(x);
1447 
1448  fDA[x-1] = data;
1449  fSigLike[x-1] = signal + nuebar + nuenonfiducial;
1450  fNumuLike[x-1] = numu + numubar;
1451  fNC[x-1] = nc;
1452  fOther[x-1] = other;
1453  }
1454 
1455 
1456  //Now perform the fit
1457  //Fit 1: Get Good Starting Fit Values, Quickly
1458  //Three parameter fit
1459  TMinuit *gMinuit = new TMinuit(3);
1460  gMinuit->SetPrintLevel(-1);
1461  gMinuit->SetFCN(fcn);
1462  Double_t arglist[4];
1463  Int_t ierflg = 0;
1464  arglist[0] = 1;
1465  gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1466 
1467  double vstart[3] = {0.8,1.2,0.9};
1468  double verror[3] = {0,0,0};
1469  double vstep[3] = {0.01,0.01,0.01};
1470 
1471  //Define Fit Parameters
1472  gMinuit->mnparm(0, "Numu-Like Scaling",
1473  vstart[0], vstep[0], 0, 0, ierflg);
1474  gMinuit->mnparm(1, "Nue-Like Scaling",
1475  vstart[1], vstep[1], 0, 0, ierflg);
1476  gMinuit->mnparm(2, "NC Scaling", vstart[2], vstep[2], 0, 0, ierflg);
1477 
1478  arglist[0] = 0; //number of iterations 0 == no limit
1479  arglist[1] = 1; //1 == standard minimization strategy
1480 
1481  //SIMPLEX should provide a good starting point for further fits
1482  gMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
1483 
1484  //Fit 2: Improved Fit + Better Eror Estimation
1485  arglist[1] = 2; //2 == try to improve minimum (slower)
1486  gMinuit->mnexcm("SET STR", arglist,1,ierflg);
1487 
1488  gMinuit->mnexcm("SIMPLEX", arglist, 3, ierflg);
1489  gMinuit->mnexcm("MIGRAD", arglist, 5, ierflg);
1490  gMinuit->mnexcm("SEEk", arglist, 4, ierflg);
1491  gMinuit->mnexcm("MINOS", arglist, 4, ierflg);
1492 
1493  //Grab Final Results
1494  double fParNewStart;
1495  double fParNewErr;
1496  gMinuit->GetParameter(0,fParNewStart,fParNewErr);
1497  vstart[0] = fParNewStart;
1498  verror[0] = fParNewErr;
1499  gMinuit->GetParameter(1,fParNewStart,fParNewErr);
1500  vstart[1] = fParNewStart;
1501  verror[1] = fParNewErr;
1502  gMinuit->GetParameter(2,fParNewStart,fParNewErr);
1503  vstart[2] = fParNewStart;
1504  verror[2] = fParNewErr;
1505 
1506  //Before Ending Fits, Make Sure There Was No Issue With the Fit
1507  bool troubleFitting = false;
1508  Double_t fmin = 0,fedm = 0,ferrdef = 0;
1509  Int_t npari =0, nparx = 0, istat = 0;
1510  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1511 
1512  Double_t corr_mat_test[3][3];
1513  gMinuit->mnemat(*corr_mat_test,3);
1514 
1515  //Turn this test error matix into a correlation matrix
1516  double corr_01 = corr_mat_test[0][1]/(corr_mat_test[0][0] *
1517  corr_mat_test[1][1]);
1518  double corr_02= corr_mat_test[0][2]/(corr_mat_test[0][0] *
1519  corr_mat_test[2][2]);
1520  double corr_12= corr_mat_test[1][2]/(corr_mat_test[1][1] *
1521  corr_mat_test[2][2]);
1522  if(istat == 2 || vstart[0] < 0 || vstart[1] < 0 || vstart[2] < 0 ||
1523  fabs(corr_01) > 0.96 || fabs(corr_02) > 0.96 || fabs(corr_12) > 0.96)
1524  troubleFitting=true;
1525 
1526  if(troubleFitting){
1527  //There Was trouble with the original fits,
1528  //most likely this was due to extremely high correlation (approx 1)
1529  //between two of the fitting parameters
1530  //Switch to a fit that adjusts all backgrounds
1531  TMinuit *jMinuit = new TMinuit(2);
1532  jMinuit->SetPrintLevel(-1);
1533  jMinuit->SetFCN(fcn2Var);
1534  arglist[0] = 1;
1535  jMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1536 
1537  jMinuit->mnparm(0, "Background", vstart[0],
1538  vstep[0], 0, 0, ierflg);
1539  jMinuit->mnparm(1, "Nue CC Scaling", vstart[1],
1540  vstep[1], 0, 0, ierflg);
1541 
1542  arglist[0] = 0; //number of iterations 0 == no limit
1543  arglist[1] = 1;
1544 
1545  jMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
1546  jMinuit->mnexcm("MIGRAD", arglist, 3, ierflg);
1547  jMinuit->mnexcm("SEEk", arglist, 3, ierflg);
1548  jMinuit->mnexcm("MINOS", arglist, 3, ierflg);
1549 
1550 
1551  jMinuit->GetParameter(0,fParNewStart,fParNewErr);
1552  vstart[0] = fParNewStart;
1553  verror[0] = fParNewErr;
1554  jMinuit->GetParameter(0,fParNewStart,fParNewErr);
1555  vstart[2] = fParNewStart;
1556  verror[2] = fParNewErr;
1557  jMinuit->GetParameter(1,fParNewStart,fParNewErr);
1558  vstart[1] = fParNewStart;
1559  verror[1] = fParNewErr;
1560  //Get the covariance matrix
1561  Double_t emat[2][2];
1562  jMinuit->mnemat(*emat,2);
1563 
1564  nue_wei = std::make_pair (vstart[1],verror[1]);
1565  numu_wei = std::make_pair (vstart[0],verror[0]);
1566  nc_wei = std::make_pair (vstart[2],verror[2]);
1567  covariance[0][0] = emat[0][0];
1568  covariance[0][1] = emat[0][1];
1569  covariance[0][2] = 0;
1570  covariance[1][0] = emat[1][0];
1571  covariance[1][1] = emat[1][1];
1572  covariance[1][2] = 0;
1573  covariance[2][0] = emat[0][0];
1574  covariance[2][1] = emat[0][1];
1575  covariance[2][2] = emat[0][0];
1576  }
1577  else{
1578  gMinuit->GetParameter(0,fParNewStart,fParNewErr);
1579  vstart[0] = fParNewStart;
1580  verror[0] = fParNewErr;
1581  gMinuit->GetParameter(2,fParNewStart,fParNewErr);
1582  vstart[2] = fParNewStart;
1583  verror[2] = fParNewErr;
1584  gMinuit->GetParameter(1,fParNewStart,fParNewErr);
1585  vstart[1] = fParNewStart;
1586  verror[1] = fParNewErr;
1587 
1588  //Get the covariance matrix
1589  Double_t emat[3][3];
1590  gMinuit->mnemat(*emat,3);
1591 
1592  nue_wei = std::make_pair (vstart[1],verror[1]);
1593  numu_wei = std::make_pair (vstart[0],verror[0]);
1594  nc_wei = std::make_pair (vstart[2],verror[2]);
1595  covariance[0][0] = emat[0][0];
1596  covariance[0][1] = emat[0][1];
1597  covariance[0][2] = emat[0][2];
1598  covariance[1][0] = emat[1][0];
1599  covariance[1][1] = emat[1][1];
1600  covariance[1][2] = emat[1][2];
1601  covariance[2][0] = emat[2][0];
1602  covariance[2][1] = emat[2][1];
1603  covariance[2][2] = emat[2][2];
1604  }
1605  }
1606  else{
1607  std::cout << "******** Skipping Fit For Bin "
1608  << caption << "\t" << signal_like << "\t"
1609  << signal_like/bkgd_amount << "\n\n\n\n" << std::endl;
1610 
1611  nue_wei = std::make_pair (1,0);
1612  numu_wei = std::make_pair (1,0);
1613  nc_wei = std::make_pair (1,0);
1614  covariance[0][0] = 0;
1615  covariance[0][1] = 0;
1616  covariance[0][2] = 0;
1617  covariance[1][0] = 0;
1618  covariance[1][1] = 0;
1619  covariance[1][2] = 0;
1620  covariance[2][0] = 0;
1621  covariance[2][1] = 0;
1622  covariance[2][2] = 0;
1623  }
1624 
1625  std::cout << "*********Fit for bin " << caption << " : " << std::endl;
1626  std::cout << "NumuCC Scaling Factor= " << numu_wei.first << " +- " <<
1627  numu_wei.second << std::endl;
1628  std::cout << "NueCC Scaling Factor= " << nue_wei.first << " +- " <<
1629  nue_wei.second << std::endl;
1630  std::cout << "NC Scaling Factor= " << nc_wei.first << " +- " <<
1631  nc_wei.second << std::endl;
1632 
1633  //Push Weights and Uncertainties to 2D Histograms
1634  hSigWeights->SetBinContent(xbin,ybin,nue_wei.first);
1635  hSigError->SetBinContent(xbin,ybin,nue_wei.second);
1636  hNumuWeights->SetBinContent(xbin,ybin,numu_wei.first);
1637  hNumuError->SetBinContent(xbin,ybin,numu_wei.second);
1638  hNCWeights->SetBinContent(xbin,ybin,nc_wei.first);
1639  hNCError->SetBinContent(xbin,ybin,nc_wei.second);
1640 
1641 
1642  //Projections For Systematic Error Band Drawing
1643  std::vector<TH1F*> hSysts;
1644  for(int z = 0; z < (int)systs_hists.size(); z++){
1645  hSysts.push_back((TH1F*)systs_hists[z]->
1646  ProjectionZ(ana::UniqueName().c_str(),
1647  xbin,xbin,ybin,ybin));
1648  }
1649  std::vector<TH1F*> hSystsUp;
1650  for(int z = 0; z < (int)systs_hists_up.size();z++){
1651  hSystsUp.push_back((TH1F*)systs_hists_up[z]->
1652  ProjectionZ(ana::UniqueName().c_str(),
1653  xbin,xbin,
1654  ybin,ybin));
1655  }
1656  std::vector<TH1F*> hSystsDown;
1657  for(int z = 0; z < (int)systs_hists_down.size();z++){
1658  hSystsDown.push_back((TH1F*)systs_hists_down[z]->
1659  ProjectionZ(ana::UniqueName().c_str(),
1660  xbin,xbin,
1661  ybin,ybin));
1662  }
1663 
1664  std::vector<TH1F*> shifts_up;
1665  std::vector<TH1F*> shifts_down;
1666 
1667  for(int z = 0; z < (int)hSysts.size();z++){
1668  shifts_up.push_back(hSysts[z]);
1669  shifts_down.push_back((TH1F*)hMCZ[0]->Clone(ana::UniqueName().c_str()));
1670  }
1671 
1672  for(int z = 0; z < (int)hSystsUp.size(); z++){
1673  shifts_up.push_back(hSystsUp[z]);
1674  shifts_down.push_back(hSystsDown[z]);
1675  }
1676 
1677  //Get Systematic Error Band
1678  TGraphAsymmErrors* g2 = PlotSystErrorBand(hMCZ[0],
1679  shifts_up,
1680  shifts_down,
1681  -1,-1,1.3,false);
1682 
1683  //Make Copies of Nominal MC template predictions
1684  std::vector<TH1F*> hMCnom;
1685  for(int z = 0; z < (int)hMCZ.size(); z++){
1686  hMCnom.push_back((TH1F*)hMCZ[z]->Clone(ana::UniqueName().c_str()));
1687  }
1688 
1689  //Scale MC templates by fitted weights
1690  hMCZ[1]->Scale(nue_wei.first); //Signal
1691  hMCZ[2]->Scale(nue_wei.first); //Nuebar
1692  hMCZ[3]->Scale(numu_wei.first); //Numu
1693  hMCZ[4]->Scale(numu_wei.first); //Numubar
1694  hMCZ[5]->Scale(nc_wei.first); //NC
1695  hMCZ[6]->Scale(1); //Other
1696  hMCZ[7]->Scale(nue_wei.first); //Non-fiducial Nue
1697 
1698 
1699  //Plot Template Before after Fit
1700  PlotFitResults(hDataZ,hMCnom,hMCZ,g2, pidName,dataName,caption,xbin,ybin);
1701 
1702  //Apply Fit Results to Analysis Histogram
1703  for(int numChn = 0; numChn < (int)hResults.size(); numChn++){
1704  for(int zbin = 1; zbin <= hResults[0]->GetZaxis()->GetNbins();zbin++){
1705  float nsig = hResults[1]->GetBinContent(xbin,ybin,zbin);
1706  float nnuebar = hResults[2]->GetBinContent(xbin,ybin,zbin);
1707  float nnumu = hResults[3]->GetBinContent(xbin,ybin,zbin);
1708  float nnumubar = hResults[4]->GetBinContent(xbin,ybin,zbin);
1709  float nnc = hResults[5]->GetBinContent(xbin,ybin,zbin);
1710  float nother = hResults[6]->GetBinContent(xbin,ybin,zbin);
1711  float nnf = hResults[7]->GetBinContent(xbin,ybin,zbin);
1712  if(numChn == 0){
1713  float binCont = nsig * nue_wei.first;
1714  binCont += nnuebar * nue_wei.first;
1715  binCont += nnumu * numu_wei.first;
1716  binCont += nnumubar * numu_wei.first;
1717  binCont += nnc * nc_wei.first;
1718  binCont += nother;
1719  binCont += nnf * nue_wei.first;
1720  hResults[0]->SetBinContent(xbin,ybin,zbin,binCont);
1721 
1722  float binErr2 =
1723  pow(covariance[0][0],2)*pow((nnumu+nnumubar),2) +
1724  pow(covariance[1][1],2)* pow((nsig+nnuebar+nnf),2) +
1725  pow(covariance[2][2],2)* pow((nnc),2) +
1726  2 * covariance[0][1] * (nnumu+nnumubar)*(nsig+nnuebar+nnf) +
1727  2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) +
1728  2 * covariance[1][2] * (nsig + nnuebar + nnf) * (nnc) +
1729  nsig * pow(nue_wei.first,2) + nnuebar * pow(nue_wei.first,2) +
1730  nnf * pow(nue_wei.first,2) + nnumu * pow(numu_wei.first,2) +
1731  nnumubar * pow(numu_wei.first,2) + nnc * pow(nc_wei.first,2) +
1732  nother;
1733 
1734  if(covariance[0][1] == 0 && covariance[0][2] == 0 &&
1735  covariance[1][2] == 0 && nue_wei.second == 0 &&
1736  numu_wei.second == 0 && nc_wei.second == 0){
1737  binErr2 = nsig + nnuebar + nnf + nnumu + nnumubar +
1738  nnc + nother;
1739  }
1740 
1741 
1742  //if(binErr2 < 0) binErr2 = fabs(binErr2);
1743 
1744  hResults[0]->SetBinError(xbin,ybin,zbin,sqrt(binErr2));
1745 
1746  std::cout << zbin << "\t" << binCont << "\t" <<
1747  sqrt(binErr2) << std::endl;
1748  if(std::isnan(sqrt((binErr2)))){
1749  std::cout << pow(covariance[0][0],2)*pow((nnumu+nnumubar),2) <<
1750  "\t" << pow(covariance[1][1],2)* pow((nsig+nnuebar+nnf),2) <<
1751  "\t" << pow(covariance[2][2],2)* pow((nnc),2) << "\t" <<
1752  2 * covariance[0][1] * (nnumu+nnumubar)*(nsig+nnuebar+nnf) <<
1753  "\t" << 2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) <<
1754  "\t" << 2 * covariance[1][2] * (nsig + nnuebar + nnf) * (nnc) <<
1755  "\t" << nsig * pow(nue_wei.first,2) << "\t" <<
1756  nnuebar * pow(nue_wei.first,2) << "\t" <<
1757  nnf * pow(nue_wei.first,2) << "\t" <<
1758  nnumu * pow(numu_wei.first,2) << "\t" <<
1759  nnumubar * pow(numu_wei.first,2) << "\t" <<
1760  nnc * pow(nc_wei.first,2) << "\t" << nother << std::endl;
1761  }
1762  }
1763  if(numChn == 1){
1764  float binCont = nsig * nue_wei.first;
1765  float binErr = sqrt( pow(nue_wei.second,2)*pow(nsig,2) +
1766  nsig * pow(nue_wei.first,2));
1767  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1768  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1769  }
1770  if(numChn == 2){
1771  float binCont = nnuebar * nue_wei.first;
1772  float binErr = sqrt( pow(nue_wei.second,2)*pow(nnuebar,2) +
1773  nnuebar * pow(nue_wei.first,2));
1774  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1775  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1776  }
1777  if(numChn == 3){
1778  float binCont = nnumu * numu_wei.first;
1779  float binErr = sqrt( pow(numu_wei.second,2)*pow(nnumu,2) +
1780  nnumu * pow(numu_wei.first,2));
1781  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1782  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1783  }
1784  if(numChn == 4){
1785  float binCont = nnumubar * numu_wei.first;
1786  float binErr = sqrt( pow(numu_wei.second,2)*pow(nnumubar,2) +
1787  nnumubar * pow(numu_wei.first,2));
1788  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1789  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1790  }
1791  if(numChn == 5){
1792  float binCont = nnc * nc_wei.first;
1793  float binErr = sqrt( pow(nc_wei.second,2)*pow(nnc,2) +
1794  nnc * pow(nc_wei.first,2));
1795  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1796  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1797  }
1798  if(numChn == 7){
1799  float binCont = nnf * nue_wei.first;
1800  float binErr = sqrt( pow(nue_wei.second,2)*pow(nnf,2) +
1801  nnf * pow(nue_wei.first,2));
1802  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1803  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1804  }
1805  }//loop over zbins
1806  } // loop over interaction channels
1807  }//loop over ybins
1808  }//loop over xbins
1809 
1810  //Draw Weights and Relative Uncertainty
1811  MakeWeightAndErrorPlot(hSigWeights,hSigError, pidName,dataName,
1812  "nue_weights", "nue_errors",
1813  "Reconstructed Electron Energy, E_{e} (GeV)",
1814  "Reconstructed cos #theta_{e}",
1815  "Signal-Like Template Fit Results",
1816  "Normalization", "Uncertainty");
1817  MakeWeightAndErrorPlot(hNumuWeights,hNumuError, pidName,dataName,
1818  "numu_weights", "numu_errors",
1819  "Reconstructed Electron Energy, E_{e} (GeV)",
1820  "Reconstructed cos #theta_{e}",
1821  "Numu-Like Template Fit Results",
1822  "Normalization", "Uncertainty");
1823  MakeWeightAndErrorPlot(hNCWeights,hNCError, pidName,dataName,
1824  "nc_weights", "nc_errors",
1825  "Reconstructed Electron Energy, E_{e} (GeV)",
1826  "Reconstructed cos #theta_{e}",
1827  "NC Template Fit Results",
1828  "Normalization", "Uncertainty");
1829 
1830  return hResults;
1831 }
1832 
1833 
1834 
hadd_many_files.counter
int counter
Definition: hadd_many_files.py:18
Simulation
void Simulation()
Definition: tools.h:16
fNumuLike
float fNumuLike[pidBins]
Definition: NueCCIncTemplateFit.h:87
stan::math::fmin
fvar< T > fmin(const fvar< T > &x1, const fvar< T > &x2)
Definition: fmin.hpp:14
lo
TSpline3 lo("lo", xlo, ylo, 12,"0")
kRed
enum BeamMode kRed
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
stan::math::fabs
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
fcn2Var
void fcn2Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: NueCCIncTemplateFit.h:195
BinByBinTemplateFit
std::vector< TH3F * > BinByBinTemplateFit(TH3F *data, std::vector< TH3F * > templates, std::vector< TH3F * > systs_hists, std::vector< TH3F * > systs_hists_up, std::vector< TH3F * > systs_hists_down, std::vector< TH3F * > analysis_templates, std::string pidName, std::string varName, std::string dataName)
Definition: NueCCIncTemplateFit.h:721
GetYaxis
correl_xv GetYaxis() -> SetDecimals()
MakeWeightAndErrorPlot
void MakeWeightAndErrorPlot(TH2F *hWeights, TH2F *hErr, std::string pidName, std::string dataName, std::string outnameWei, std::string outnameErr, std::string ytitle, std::string xtitle, std::string title, std::string ztitle_wei, std::string ztitle_err, float x_lo=0.75, float x_hi=1.0, float y_lo=1.0, float y_hi=10.0)
Definition: NueCCIncTemplateFit.h:251
BinByBinTemplateFit_TemplateResults
std::vector< TH3F * > BinByBinTemplateFit_TemplateResults(TH3F *data, std::vector< TH3F * > templates, std::vector< TH3F * > systs_hists, std::vector< TH3F * > systs_hists_up, std::vector< TH3F * > systs_hists_down, std::vector< TH3F * > analysis_templates, std::string pidName, std::string varName, std::string dataName)
Definition: NueCCIncTemplateFit.h:1283
make_syst_table_plots.xtitle
string xtitle
Definition: make_syst_table_plots.py:204
fSigLike
float fSigLike[pidBins]
Definition: NueCCIncTemplateFit.h:86
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
plot_validation_datamc.Clone
def Clone(obj)
Definition: plot_validation_datamc.py:15
cet::pow
constexpr T pow(T x)
Definition: pow.h:75
CalculateCovarianceMatrix
void CalculateCovarianceMatrix(TH1F *nom, std::vector< TH1F * > hSysts, int nPseudo, std::string pidName, int xbin, int ybin)
Definition: NueCCIncTemplateFit.h:467
par
Int_t par
Definition: SimpleIterate.C:24
covariance
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
Definition: chains_test.cpp:464
om::cerr
OStream cerr
Definition: OStream.cxx:7
make_pair
std::pair< Spectrum *, CheatDecomp * > make_pair(SpectrumLoaderBase &loader_data, SpectrumLoaderBase &loader_mc, HistAxis *axis, Cut *cut, const SystShifts &shift, const Var &wei)
Definition: DataMCLoad.C:336
PandAna.Demos.pi0_spectra.labels
labels
Definition: pi0_spectra.py:178
fcn
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: NueCCIncTemplateFit.h:134
npar
static const unsigned int npar
Definition: CherenkovCalc.cxx:46
s_b_ratio
const float s_b_ratio
Definition: NueCCIncTemplateFit.h:80
submit_syst.y
y
Definition: submit_syst.py:249
std::isnan
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
fDA
float fDA[pidBins]
Definition: NueCCIncTemplateFit.h:85
ana::pidBins
const Binning pidBins
Definition: NumuCC_CPiBin.h:5
kNueCCColorDef
const std::vector< int > kNueCCColorDef
Definition: NueCCIncTemplateFit.h:91
data
const XML_Char const XML_Char * data
Definition: expat.h:268
std::swap
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
Definition: HLTGlobalStatus.h:199
demo5.c2
c2
Definition: demo5.py:33
stat_covMatrix
double stat_covMatrix[pidBins][pidBins]
Definition: NueCCIncTemplateFit.h:79
make_true_q0q3_plots.caption
caption
Definition: make_true_q0q3_plots.py:404
hi
TSpline3 hi("hi", xhi, yhi, 18,"0")
submit_syst.x
x
Definition: submit_syst.py:249
value
const XML_Char int const XML_Char * value
Definition: expat.h:331
col
Int_t col[ntarg]
Definition: Style.C:29
cv
const std::string cv[Ncv]
Definition: Calculate_covariances.C:39
fNC
float fNC[pidBins]
Definition: NueCCIncTemplateFit.h:88
shouldRebin
const bool shouldRebin
Definition: NueCCIncTemplateFit.h:77
calib::j
const double j
Definition: BetheBloch.cxx:29
TMatrixD
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
getCaption2D
std::string getCaption2D(TH3F *sample, int xbin, int ybin)
Definition: NueCCIncTemplateFit.h:103
bin
float bin[41]
Definition: plottest35.C:14
submit_syst.str
str
Definition: submit_syst.py:150
prelim
TLatex * prelim
Definition: Xsec_final.C:133
PlotSystErrorBand
TGraphAsymmErrors * PlotSystErrorBand(TH1F *&nom, std::vector< TH1F * > &ups, std::vector< TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)
Definition: NueCCIncTemplateFit.h:417
test.z
z
Definition: test.py:28
covMatrix
double covMatrix[pidBins][pidBins]
Definition: NueCCIncTemplateFit.h:78
signalRegion
const float signalRegion
Definition: NueCCIncTemplateFit.h:83
om::cout
OStream cout
Definition: OStream.cxx:6
make_mec_shifts_plots.GetMaximum
def GetMaximum(hist)
Definition: make_mec_shifts_plots.py:20
nc
enum BeamMode nc
PrintCaption
void PrintCaption(TString str)
Definition: NueCCIncTemplateFit.h:124
demo5.c1
c1
Definition: demo5.py:24
dataName
std::string dataName
Definition: make_starplots.C:22
ana::AutoPlaceLegend
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:717
PlotFitResults
void PlotFitResults(TH1F *hData, std::vector< TH1F * >hMCnom, std::vector< TH1F * >hMCwei, TGraphAsymmErrors *g, std::string pidName, std::string dataName, std::string caption, int xbin, int ybin)
Definition: NueCCIncTemplateFit.h:299
pad2
TPad * pad2
Definition: analysis.C:13
fOther
float fOther[pidBins]
Definition: NueCCIncTemplateFit.h:89
CalculateTotalCovariance
void CalculateTotalCovariance(TH1F *nominal, std::vector< TH1F * > hSysts, std::string pidName, int xbin, int ybin)
Definition: NueCCIncTemplateFit.h:541
kGreen
enum BeamMode kGreen
kBlue
enum BeamMode kBlue
w
Float_t w
Definition: plot.C:20
ana::UniqueName
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
pad1
TPad * pad1
Definition: analysis.C:13
string
enum BeamMode string
uint
unsigned int uint
Definition: DAQChannelMapBaseOLD.h:22
Generated on Sun Mar 7 2021 14:13:07 for NOvA by   doxygen 1.8.11