NuebarCCIncTemplateFit.h
Go to the documentation of this file.
1 #pragma once
2 
4 #include "CAFAna/Core/Spectrum.h"
5 #include "CAFAna/Core/ISyst.h"
16 #include "CAFAna/XSec/Flux.h"
20 #include "CAFAna/Analysis/Plots.h"
21 
22 //ROOT Includes
23 #include "Utilities/rootlogon.C"
24 #include "TH1.h"
25 #include "TH2.h"
26 #include "THStack.h"
27 #include "TCanvas.h"
28 #include "TLegend.h"
29 #include "TStyle.h"
30 #include "TFile.h"
31 #include <vector>
32 #include <math.h>
33 #include "TLine.h"
34 #include "TFile.h"
35 #include "TTree.h"
36 #include "TCanvas.h"
37 #include "TF1.h"
38 #include "TH1F.h"
39 #include "TH2F.h"
40 #include "TF2.h"
41 #include "TH2.h"
42 #include "TCutG.h"
43 #include "TMath.h"
44 #include "TCanvas.h"
45 #include "TStyle.h"
46 #include "TRandom.h"
47 #include "TGraph.h"
48 #include "TColor.h"
49 #include "TF3.h"
50 #include "TH3.h"
51 #include "TH3F.h"
52 #include "TMinuit.h"
53 #include "TLegend.h"
54 #include "THStack.h"
55 #include "TMatrix.h"
56 #include "TGraphErrors.h"
57 #include "TGraph.h"
58 #include "TDecompLU.h"
59 #include "TDecompSVD.h"
60 #include "TRandom.h"
61 
62 //C++
63 #include <iostream>
64 #include <fstream>
65 #include <iomanip>
66 #include <vector>
67 #include <cmath>
68 #include <math.h>
69 #include <string>
70 #include <stdint.h>
71 #include <sstream>
72 
73 using namespace ana;
74 
75 const int pidBins = 21;
76 const bool shouldRebin = false;
79 const float s_b_ratio = 0.5; //0.40; //ratio of signal/bkgd (in signal region)
80  // that needs to be exceeded
81  // to perform fit in analysis bin
82 const float signalRegion = 0.85;
83 
84 float fDA[pidBins];
87 float fNC[pidBins];
88 float fOther[pidBins];
89 
90 const std::vector<int> kNueCCColorDef = {
91  632, //kRed, Total MC
92  600, //kBlue Signal
93  616, //kMagenta nuebar CC
94  416-3, //kGreen - 3 NumuCC
95  416+3, //kGreen + 3 Anti NumuCC
96  800, //kOrange NC
97  400-3, //kYellow - 3 Other
98  416-3, //kGreen - 3 Total Background
99 };
100 
101 
102 std::string getCaption2D(TH3F* sample, int xbin, int ybin)
103 {
104  std::stringstream low_X;
105  low_X << std::fixed << std::setprecision(2) <<
106  sample->GetXaxis()->GetBinLowEdge(xbin);
107  std::stringstream hi_X;
108  hi_X << std::fixed << std::setprecision(2) <<
109  sample->GetXaxis()->GetBinUpEdge(xbin);
110 
111  std::stringstream low_Y;
112  low_Y << std::fixed << std::setprecision(2) <<
113  sample->GetYaxis()->GetBinLowEdge(ybin);
114  std::stringstream hi_Y;
115  hi_Y << std::fixed << std::setprecision(2) <<
116  sample->GetYaxis()->GetBinUpEdge(ybin);
117 
118  std::string caption = low_X.str() + " #leq cos #theta_{e} < "+
119  hi_X.str() + " " + low_Y.str() + " #leq E_{e} < " + hi_Y.str();
120  return caption;
121 }
122 
123 void PrintCaption(TString str)
124 {
125 
126  TLatex* prelim = new TLatex(0.23, .80, str);
127  prelim->SetTextColor(kBlack);
128  prelim->SetNDC();
129  prelim->SetTextSize(0.05);
130  prelim->Draw();
131 }
132 
133 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){
134 
135  double chisquare = 0;
136  float data = 0;
137  float signal = 0;
138  float numu = 0;
139  float nc = 0;
140  float other = 0;
141 
142  // Covariance Matrix Inversion First
143  TMatrixF Vij_stat = TMatrix(pidBins,pidBins);
144  TMatrixF Vij = TMatrix(pidBins,pidBins);
145  for(int i =0; i < pidBins; i++){
146  for(int j = 0; j < pidBins; j++){
147  Vij_stat[i][j] = stat_covMatrix[i][j];
148  Vij[i][j] = covMatrix[i][j];
149  }
150  }
151 
152  //Inversion
153  TMatrixD H3 = Vij_stat + Vij;
154  TDecompSVD svd(H3);
155  TMatrixD VijInv = svd.Invert();
156 
157 
158  // float delta_i = 0;
159  // float delta_j = 0;
160 
161  TMatrixD FitVector = TMatrixD(1,pidBins);
162  TMatrixD FitVector2 = TMatrixD(pidBins,1);
163 
164  //Calculate Chi Sq
165  for(int i =0; i < pidBins; i++){
166  data = fDA[i];
167  signal = fSigLike[i];
168  numu = fNumuLike[i];
169  nc = fNC[i];
170  other = fOther[i];
171 
172  if(data < 30){
173  FitVector[0][i] = 0;
174  FitVector2[i][0] = 0;
175  continue;
176  }
177 
178  FitVector[0][i] = ( data - (par[1]*signal + par[0]*numu +
179  par[2]*nc + other));
180 
181  FitVector2[i][0] = ( data -(par[1]*signal + par[0]*numu +
182  par[2]*nc + other));
183  }
184 
185  TMatrixD Chi = FitVector * VijInv;
186  TMatrixD Chi2 = Chi * FitVector2;
187  float secondterm = Chi2(0,0);
188 
189  chisquare = secondterm;
190  f = chisquare;
191 
192 }
193 
194 void fcn2Var(Int_t &npar, Double_t *gin,
195  Double_t &f, Double_t *par, Int_t iflag){
196 
197  // double chisquare = 0;
198  float data = 0;
199  float signal = 0;
200  float numu = 0;
201  float nc = 0;
202  float other = 0;
203 
204  // Covariance Matrix Inversion First
205  TMatrixF Vij_stat = TMatrix(pidBins,pidBins);
206  TMatrixF Vij = TMatrix(pidBins,pidBins);
207  for(int i =0; i < pidBins; i++){
208  for(int j = 0; j < pidBins; j++){
209  Vij_stat[i][j] = stat_covMatrix[i][j];
210  Vij[i][j] = covMatrix[i][j];
211  }
212  }
213 
214  //Inversion
215  TMatrixD H3 = Vij_stat + Vij;
216  TDecompSVD svd(H3);
217  TMatrixD VijInv = svd.Invert();
218 
219  TMatrixD FitVector = TMatrixD(1,pidBins);
220  TMatrixD FitVector2 = TMatrixD(pidBins,1);
221 
222  //Calculate Chi Sq
223  for(int i =0; i < pidBins; i++){
224  data = fDA[i];
225  signal = fSigLike[i];
226  numu = fNumuLike[i];
227  nc = fNC[i];
228  other = fOther[i];
229 
230  if(data < 30){
231  FitVector[0][i] = 0;
232  FitVector2[i][0] = 0;
233  continue;
234 
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 
252 
253 
254 
255 void PlotFitResults(TH1F* hData,
256  std::vector<TH1F*>hMCnom, std::vector<TH1F*>hMCwei,
257  TGraphAsymmErrors* g, std::string pidName,
259  int xbin,int ybin)
260 {
261 
262  TH1F* hMC_nom_signal = (TH1F*)hMCnom[1]->Clone(ana::UniqueName().c_str());
263  hMC_nom_signal->Add(hMCnom[2],1);
264  hMC_nom_signal->Add(hMCnom[7],1);
265  TH1F* hMC_nom_bkgd = (TH1F*)hMCnom[3]->Clone(ana::UniqueName().c_str());
266  hMC_nom_bkgd->Add(hMCnom[4],1);
267  hMC_nom_bkgd->Add(hMCnom[5],1);
268  hMC_nom_bkgd->Add(hMCnom[6],1);
269 
270  TH1F* hMC_wgt_signal = (TH1F*)hMCwei[1]->Clone(ana::UniqueName().c_str());
271  hMC_wgt_signal->Add(hMCwei[2],1);
272  hMC_wgt_signal->Add(hMCwei[7],1);
273  TH1F* hMC_wgt_bkgd = (TH1F*)hMCwei[3]->Clone(ana::UniqueName().c_str());
274  hMC_wgt_bkgd->Add(hMCwei[2],1);
275  hMC_wgt_bkgd->Add(hMCwei[4],1);
276  hMC_wgt_bkgd->Add(hMCwei[5],1);
277  hMC_wgt_bkgd->Add(hMCwei[6],1);
278 
279  hMC_nom_signal->SetLineColor(kBlue+3);
280  hMC_nom_bkgd->SetLineColor(kGreen+3);
281 
282  TH1F* hMCWeight = (TH1F*)hMCwei[6]->Clone(ana::UniqueName().c_str());
283  hMCWeight->Add(hMCwei[1]);
284  hMCWeight->Add(hMCwei[2]);
285  hMCWeight->Add(hMCwei[3]);
286  hMCWeight->Add(hMCwei[4]);
287  hMCWeight->Add(hMCwei[5]);
288  hMCWeight->Add(hMCwei[7]);
289  hMCWeight->SetLineColor(kRed);
290 
291  hMCnom[0]->SetLineColor(kMagenta-3);
292  hData->SetMarkerStyle(20);
293  hData->SetTitle(caption.c_str());
294  hMCnom[0]->SetTitle(caption.c_str());
295 
296  if(hData->GetMaximum() > hMCnom[0]->GetMaximum())
297  hMCnom[0]->GetYaxis()->SetRangeUser(0,hData->GetMaximum()*1.3);
298  else(hMCnom[0]->GetYaxis()->SetRangeUser(0,hMCnom[0]->GetMaximum()*1.3));
299 
300  TCanvas* cCompare = new TCanvas(ana::UniqueName().c_str(),"cCompare");
301  cCompare->SetLeftMargin(0.10);
302  cCompare->SetRightMargin(0.10);
303  TPad *pad1 = new TPad(ana::UniqueName().c_str(), "pad1",0,0,1,1);
304  TPad *pad2 = new TPad(ana::UniqueName().c_str(), "pad2",0,0,1,1);
305  pad1->SetFillStyle(0);
306  pad2->SetFillStyle(0);
307  pad1->SetBottomMargin(0.30);
308  pad2->SetTopMargin(1-0.30);
309  pad2->SetGridx();
310  pad2->SetGridy();
311  cCompare->cd();
312  pad1->Draw();
313  pad1->cd();
314  hMCWeight->SetTitle("");
315  hMCnom[0]->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} POT");
316  hMCnom[0]->GetXaxis()->SetLabelColor(kWhite);
317  hMCnom[0]->GetYaxis()->SetLabelSize(0.03);
318  hMCnom[0]->GetYaxis()->SetTitleSize(0.035);
319  hMCnom[0]->SetTitle("");
320  hMCnom[0]->Draw("hist");
321  PrintCaption(caption.c_str());
322  g->Draw("e2 same");
323  hMCnom[0]->Draw("hist same");
324  hMCWeight->Draw("hist same");
325  hData->Draw("same");
326  hMC_wgt_signal->Draw("hist same");
327  hMC_nom_signal->Draw("hist same");
328  hMC_wgt_bkgd->Draw("hist same");
329  hMC_nom_bkgd->Draw("hist same");
330  hData->Draw("same");
331  TLegend* leg = ana::AutoPlaceLegend(0.3,0.3);
332  leg->AddEntry(hData, "Fake Data", "p");
333  leg->AddEntry(hMCWeight, "Scaled Total MC", "l");
334  leg->AddEntry(hMC_wgt_signal, "Scaled Signal", "l");
335  leg->AddEntry(hMC_wgt_bkgd,"Scaled Background", "l");
336  leg->AddEntry(hMCnom[0], "Total MC", "l");
337  leg->AddEntry(hMC_nom_signal, "Signal", "l");
338  leg->AddEntry(hMC_nom_bkgd,"Background", "l");
339  leg->Draw();
340  Simulation();
341  gPad->RedrawAxis();
342  cCompare->cd();
343  pad2->Draw();
344  pad2->cd();
345  TH1F* hratio = (TH1F*)hMCWeight->Clone(ana::UniqueName().c_str());
346  TH1F* hratio2 = (TH1F*)hMCnom[0]->Clone(ana::UniqueName().c_str());
347  hratio->GetYaxis()->SetTitle("MC/Data");
348  hratio->GetXaxis()->SetTitle("Electron Prong CVN Score");
349  hratio->GetXaxis()->CenterTitle();
350  hratio->SetMarkerColor(kBlack);
351  hratio->SetLineColor(kBlack);
352  hratio->SetMarkerStyle(20);
353  hratio->Divide(hData);
354  hratio2->Divide(hData);
355  hratio->GetYaxis()->SetLabelSize(0.03);
356  hratio->GetYaxis()->SetTitleSize(0.035);
357  hratio->GetYaxis()->SetRangeUser(0.4,1.5);
358  hratio->Draw("same");
359  hratio2->Draw("hist same");
360  //hratio2->Draw("hist same");
361  pad2->Draw("same");
362  gPad->RedrawAxis();
363  cCompare->cd();
364  cCompare->Update();
365  char outname[50];
366  sprintf(outname, "%s%s_%s_postfit_ratio_bin_%i_%i.png", "Plots/",
367  pidName.c_str(),
368  dataName.c_str(),xbin,ybin);
369  cCompare->SaveAs(outname);
370  cCompare->Close();
371 }
372 
373 TGraphAsymmErrors* PlotSystErrorBand(TH1F*& nom,
374  std::vector<TH1F*>& ups,
375  std::vector<TH1F*>& dns,
376  int col, int errCol,
377  float headroom, bool newaxis)
378 {
379  if(col == -1){
380  col = kRed;
381  errCol = kRed - 9;
382  }
383  else if(errCol == -1) errCol = col-9; // hopefully a lighter version
384 
385  nom->SetLineColor(col);
386  nom->GetXaxis()->CenterTitle();
387  nom->GetYaxis()->CenterTitle();
388 
389  // double yMax = nom->GetBinContent(nom->GetMaximumBin());
390 
391  TGraphAsymmErrors* g = new TGraphAsymmErrors;
392 
393  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
394  const double y = nom->GetBinContent(binIdx);
395  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
396 
397  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
398 
399  double errUp = 0;
400  double errDn = 0;
401 
402  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
403  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
404  double lo = y-dns[systIdx]->GetBinContent(binIdx);
405 
406 
407  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
408  if(lo < 0 && hi > 0){ // if both are high, take largest deviation
409  double max = std::max(abs(lo), abs(hi));
410  lo = 0;
411  hi = max;
412  }
413  if(lo > 0 && hi < 0) { // if both are low, take largest deviation
414  double max = std::max(abs(lo), abs(hi));
415  lo = max;
416  hi = 0;
417  }
418 
419  errUp += hi*hi;
420  errDn += lo*lo;
421 
422  // TODO: what happens if they're both high or both low?
423  } // end for systIdx
424 
425  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
426  } // end for i
427 
428  g->SetFillColor(errCol);
429  return g;
430 }
431 
432 
433 void CalculateCovarianceMatrix(TH1F* nom,std::vector<TH1F*> hSysts,
434  int nPseudo,std::string pidName,
435  int xbin,int ybin)
436 {
437  std::vector<TH1F*> vPseudoExp;
438 
439  TH2F* covariance_hist =
440  new TH2F(ana::UniqueName().c_str(),
441  "Covariance Matrix;PID Bin;PID Bin;Covariance",
442  nom->GetXaxis()->GetNbins(),
443  0,nom->GetXaxis()->GetNbins(),
444  nom->GetXaxis()->GetNbins(),
445  0,nom->GetXaxis()->GetNbins());
446 
447 
448  TRandom *r0 = new TRandom();
449  for(int i = 0; i < nPseudo; i++){
450  TH1F* hClone = (TH1F*)nom->Clone(ana::UniqueName().c_str());
451 
452  for(uint j = 0; j < hSysts.size(); j++){
453  double wei = r0->Gaus(0,1);
454  //std::cout << wei << std::endl;
455  hClone->Add(hSysts[j],wei);
456  }
457 
458  vPseudoExp.push_back(hClone);
459  }
460 
461  //Mij = < (Val_i)*(Val_j)> - <Val_i> <Val_j>
462  //covMatrix has to be filled
463  for(int j = 1; j <= hSysts[0]->GetXaxis()->GetNbins(); j++){
464  for(int k = 1; k <= hSysts[0]->GetXaxis()->GetNbins(); k++){
465  double multiple = 0;
466  double x_j = 0;
467  double x_k = 0;
468  double counter = 0;
469  for(int i = 0; i < nPseudo; i++){
470  double val_j = vPseudoExp[i]->GetBinContent(j);
471  double val_k = vPseudoExp[i]->GetBinContent(k);
472  multiple += val_j * val_k;
473  x_j += val_j;
474  x_k += val_k;
475  counter++;
476  }
477  double covariance = multiple/counter - ((x_j/counter)*(x_k/counter));
478  covMatrix[j-1][k-1] = covariance;
479  covariance_hist->SetBinContent(j,k,covariance);
480  }
481  }
482 
483  char out[50];
484  sprintf(out, "Plots/covariance_%s_%i_%i.png", pidName.c_str(),xbin,ybin);
485  TCanvas *c1 = new TCanvas("c1","c1");
486  c1->SetLeftMargin(0.15);
487  c1->SetRightMargin(0.15);
488  covariance_hist->GetXaxis()->CenterTitle();
489  covariance_hist->GetYaxis()->CenterTitle();
490  covariance_hist->GetZaxis()->CenterTitle();
491  covariance_hist->Draw("COLZ");
492  c1->SaveAs(out);
493  delete c1;
494  return;
495 }
496 
497 void CalculateTotalCovariance(TH1F* nom,std::vector<TH1F*> hSysts,
498  std::string pidName,
499  int xbin,int ybin)
500 {
501  std::vector<TH1F*> vPseudoExp;
502 
503  std::cout << "Hey Hey Right Here " << hSysts.size() << std::endl;
504 
505  TH2F* covariance_hist =
506  new TH2F(ana::UniqueName().c_str(),
507  "Covariance Matrix;Systematic Number;Systematic Number;Covariance",
508  (int)hSysts.size(),
509  0,(int)hSysts.size(),
510  (int)hSysts.size(),
511  0,(int)hSysts.size());
512 
513  TH2F* correlation_hist =
514  new TH2F(ana::UniqueName().c_str(),
515  "Correlation Matrix;;;Correlation",
516  (int)hSysts.size(),
517  0,(int)hSysts.size(),
518  (int)hSysts.size(),
519  0,(int)hSysts.size());
520 
521 
522  std::vector<double> vAverages;
523  for(int j = 0; j < (int)hSysts.size(); j++){
524  double avg = 0;
525  for(int i = 0; i < nom->GetXaxis()->GetNbins(); i++){
526  avg += hSysts[j]->GetBinContent(i);
527  }
528  avg /= nom->GetXaxis()->GetNbins();
529  vAverages.push_back(avg);
530  }
531 
532  for(int i = 0; i < (int)vAverages.size(); i++){
533  std::cout << vAverages[i] << std::endl;
534  }
535 
536  const char *labels[6] = {"Cal.Shp.","Ckv","Calib","Light","XSec","PPFX"};
537 
538 
539  //Sample Covariance = 1/N Sum( (x_i - <x>) (y_i - <y>))
540  for(int j = 0; j < (int) hSysts.size(); j++){
541  for(int k = 0; k < (int)hSysts.size(); k++){
542  double summed_term = 0;
543  double nterms = 0;
544  for(int i = 0; i < nom->GetXaxis()->GetNbins(); i++){
545  summed_term += (hSysts[j]->GetBinContent(i) - vAverages[j]) *
546  (hSysts[k]->GetBinContent(i) - vAverages[k]);
547  nterms++;
548  }
549  double covariance = 1./nterms * summed_term;
550  covariance_hist->SetBinContent(j+1,k+1,covariance);
551  std::cout << j+1 << " " << k+1 << " " << covariance << std::endl;
552  }
553  }
554 
555  for(int j = 1; j <= covariance_hist->GetXaxis()->GetNbins(); j++){
556  for(int k = 1; k <= covariance_hist->GetXaxis()->GetNbins(); k++){
557  double covariance = covariance_hist->GetBinContent(j,k);
558  double correlation =
559  covariance/(sqrt(covariance_hist->GetBinContent(j,j)) *
560  sqrt(covariance_hist->GetBinContent(k,k)));
561  correlation_hist->Fill(labels[j-1],labels[k-1],correlation);
562  }
563  }
564 
565 
566  char out[50];
567  sprintf(out, "Plots/covariance_systs_%s_%i_%i.png", pidName.c_str(),xbin,ybin);
568  TCanvas *c1 = new TCanvas("c1","c1");
569  c1->SetLeftMargin(0.15);
570  c1->SetRightMargin(0.15);
571  covariance_hist->GetXaxis()->CenterTitle();
572  covariance_hist->GetYaxis()->CenterTitle();
573  covariance_hist->GetZaxis()->CenterTitle();
574  covariance_hist->Draw("COLZ");
575  c1->SaveAs(out);
576  delete c1;
577 
578 
579  correlation_hist->LabelsDeflate("X");
580  correlation_hist->LabelsDeflate("Y");
581  correlation_hist->LabelsOption("v");
582  sprintf(out, "Plots/correlation_systs_%s_%i_%i.png",
583  pidName.c_str(),xbin,ybin);
584  TCanvas *c2 = new TCanvas("c2","c2");
585  c2->SetLeftMargin(0.15);
586  c2->SetRightMargin(0.15);
587  correlation_hist->GetXaxis()->CenterTitle();
588  correlation_hist->GetYaxis()->CenterTitle();
589  correlation_hist->GetZaxis()->CenterTitle();
590  correlation_hist->Draw("COLZ text");
591  c2->SaveAs(out);
592  delete c2;
593  return;
594 
595 }
596 
597 std::vector<TH3F*> BinByBinTemplateFit(TH3F* data,
598  std::vector<TH3F*> templates,
599  std::vector<TH3F*> systs_hists,
600  std::vector<TH3F*> systs_hists_up,
601  std::vector<TH3F*> systs_hists_down,
602  std::vector<TH3F*> analysis_templates,
603  std::string pidName, std::string varName,
605 {
606  std::vector<TH3F*> badentry;
607  if(systs_hists_up.size() != systs_hists_down.size()){
608  std::cerr << "Unequal up and down shifts. Exiting" << std::endl;
609  }
610 
611  std::vector<TH3F*> hResults;
612  for(uint i = 0; i < analysis_templates.size(); i++){
613  hResults.push_back((TH3F*)analysis_templates[i]->
614  Clone(ana::UniqueName().c_str()));
615  }
616 
617 
618  for(int xbin = 0; xbin <= data->GetXaxis()->GetNbins(); xbin++){
619  for(int ybin = 0; ybin <= data->GetYaxis()->GetNbins(); ybin++){
620 
621  std::string caption = getCaption2D(templates[0],xbin,ybin);
622 
623 
624  //Variables to hold fit results
625  std::pair<double,double> nue_wei = std::make_pair (0,0);
626  std::pair<double,double> numu_wei = std::make_pair (0,0);
627  std::pair<double,double> nc_wei = std::make_pair (0,0);
628  Double_t covariance[3][3];
629 
630  //Project Data and MC to PID axis for template fit
631  TH1F* hDataZ = (TH1F*)data->ProjectionZ(ana::UniqueName().c_str(),
632  xbin,xbin,ybin,ybin);
633 
634  std::vector<TH1F*> hMCZ;
635  for(uint numChns = 0; numChns < templates.size(); numChns++){
636  hMCZ.push_back((TH1F*)templates[numChns]->
637  ProjectionZ(ana::UniqueName().c_str(),
638  xbin,xbin,ybin,ybin));
639  hMCZ[numChns]->SetLineColor(kNueCCColorDef[numChns]);
640  }
641 
642  //Make Sure Fit should be performed in this bin (xbin,ybin)
643  int pid_bin = hDataZ->GetXaxis()->FindBin(signalRegion);
644  float signal_like = hMCZ[1]->Integral(pid_bin,-1) +
645  hMCZ[2]->Integral(pid_bin,-1) + hMCZ[7]->Integral(pid_bin,-1);
646  float bkgd_amount = hMCZ[0]->Integral(pid_bin,-1) - signal_like;
647 
648  if(signal_like > 100 &&
649  signal_like/bkgd_amount > s_b_ratio){
650 
651  //Calculate Uncertainties due to each systematic sample
652  std::vector<TH1F*> hSystsZ;
653 
654  //For One-Sided Systematics
655  for(int systNum = 0; systNum < (int) systs_hists.size(); systNum++){
656  TH1F* hHolder =
657  (TH1F*)systs_hists[systNum]->ProjectionZ(ana::UniqueName().c_str(),
658  xbin,xbin,ybin,ybin);
659  //Loop Over Bins, Calculate the absolute distance of systematic from
660  //Nominal MC
661  for(int bin = 1; bin <= hHolder->GetXaxis()->GetNbins(); bin++){
662  double err = hHolder->GetBinContent(bin);
663  double cv = hMCZ[0]->GetBinContent(bin);
664  hHolder->SetBinContent(bin,fabs(err-cv));
665  }
666  hSystsZ.push_back(hHolder);
667  }
668 
669  //For Two-Side Systematics
670  // Lets take a conservative approach, so furthest value away from
671  // nominal is taken as the uncertainty for that bin. Should deal with
672  // highly asymmetric systematics better
673  for(int systNum = 0; systNum < (int)systs_hists_up.size(); systNum++){
674  TH1F* hHolder = (TH1F*)systs_hists_up[systNum]->
675  ProjectionZ(ana::UniqueName().c_str(),xbin,xbin,ybin,ybin);
676  TH1F* hHolder_dwn = (TH1F*)systs_hists_down[systNum]->
677  ProjectionZ(ana::UniqueName().c_str(),xbin,xbin,ybin,ybin);
678 
679  //Loop over bins, convert two-sided uncertainty into a single value
680  //which is the largest distance of each shift from nominal mc
681  for(int bin = 1; bin <= hHolder->GetXaxis()->GetNbins(); bin++){
682  double hi = hHolder->GetBinContent(bin);
683  double lo = hHolder_dwn->GetBinContent(bin);
684  double cv = hMCZ[0]->GetBinContent(bin);
685  double val_hi = fabs(hi - cv);
686  double val_lo = fabs(lo - cv);
687  // double value = 0.5*(fabs(cv-hi) + fabs(cv-lo));
688  double value = (val_hi > val_lo)? val_hi : val_lo;
689  hHolder->SetBinContent(bin,value);
690  }
691  hSystsZ.push_back(hHolder);
692  }
693 
694  CalculateTotalCovariance(hMCZ[0],hSystsZ,pidName,xbin,ybin);
695  //Calculate Systematic Covariance Matrix
696  //Throw nPseudo experiments
697  int nPseudo = 1000;
698  CalculateCovarianceMatrix(hMCZ[0],hSystsZ,nPseudo,pidName,xbin,ybin);
699 
700  //Calculate Statistical Covariance Matrix
701  //Diagonal with Data statistics for each pid bin
702  for(int jj = 1; jj <= hSystsZ[0]->GetXaxis()->GetNbins(); jj++){
703  for(int kk = 1; kk <= hSystsZ[0]->GetXaxis()->GetNbins(); kk++){
704  if(jj==kk){
705  stat_covMatrix[jj-1][kk-1] =
706  hDataZ->GetBinError(jj)*hDataZ->GetBinError(jj)
707  +hMCZ[0]->GetBinError(jj)*hMCZ[0]->GetBinError(jj);
708  }
709  else{
710  stat_covMatrix[jj-1][kk-1] = 0.0;
711  }
712  }
713  }
714 
715 
716  //Arrays are most(?) efficienct way to feed values to Minuit
717  //Fill Arrays for MinuitFit
718  for(int x = 1; x <= hDataZ->GetXaxis()->GetNbins(); x++){
719  float data = hDataZ->GetBinContent(x);
720  // float signal = hMCZ[1]->GetBinContent(x);
721  float nue = hMCZ[1]->GetBinContent(x);
722  float nuebar = hMCZ[2]->GetBinContent(x);
723  float numu = hMCZ[3]->GetBinContent(x);
724  float numubar = hMCZ[4]->GetBinContent(x);
725  float nc = hMCZ[5]->GetBinContent(x);
726  float other = hMCZ[6]->GetBinContent(x);
727  float nuenonfiducial = hMCZ[7]->GetBinContent(x);
728 
729  fDA[x-1] = data;
730  fSigLike[x-1] = nue + nuebar + nuenonfiducial;
731  fNumuLike[x-1] = numu + numubar;
732  fNC[x-1] = nc;
733  fOther[x-1] = other;
734  }
735 
736 
737  //Now perform the fit
738  //Fit 1: Get Good Starting Fit Values, Quickly
739  //Three parameter fit
740  TMinuit *gMinuit = new TMinuit(3);
741  gMinuit->SetPrintLevel(-1);
742  gMinuit->SetFCN(fcn);
743  Double_t arglist[4];
744  Int_t ierflg = 0;
745  arglist[0] = 1;
746  gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
747 
748  double vstart[3] = {0.8,1.2,0.9};
749  double verror[3] = {0,0,0};
750  double vstep[3] = {0.01,0.01,0.01};
751 
752  //Define Fit Parameters
753  gMinuit->mnparm(0, "Numu-Like Scaling",
754  vstart[0], vstep[0], 0, 0, ierflg);
755  gMinuit->mnparm(1, "Nue-Like Scaling",
756  vstart[1], vstep[1], 0, 0, ierflg);
757  gMinuit->mnparm(2, "NC Scaling", vstart[2], vstep[2], 0, 0, ierflg);
758 
759  arglist[0] = 0; //number of iterations 0 == no limit
760  arglist[1] = 1; //1 == standard minimization strategy
761 
762  //SIMPLEX should provide a good starting point for further fits
763  gMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
764 
765  //Fit 2: Improved Fit + Better Eror Estimation
766  arglist[1] = 2; //2 == try to improve minimum (slower)
767  gMinuit->mnexcm("SET STR", arglist,1,ierflg);
768 
769  gMinuit->mnexcm("SIMPLEX", arglist, 3, ierflg);
770  gMinuit->mnexcm("MIGRAD", arglist, 5, ierflg);
771  gMinuit->mnexcm("SEEk", arglist, 4, ierflg);
772  gMinuit->mnexcm("MINOS", arglist, 4, ierflg);
773 
774  //Grab Final Results
775  double fParNewStart;
776  double fParNewErr;
777  gMinuit->GetParameter(0,fParNewStart,fParNewErr);
778  vstart[0] = fParNewStart;
779  verror[0] = fParNewErr;
780  gMinuit->GetParameter(1,fParNewStart,fParNewErr);
781  vstart[1] = fParNewStart;
782  verror[1] = fParNewErr;
783  gMinuit->GetParameter(2,fParNewStart,fParNewErr);
784  vstart[2] = fParNewStart;
785  verror[2] = fParNewErr;
786 
787  //Before Ending Fits, Make Sure There Was No Issue With the Fit
788  bool troubleFitting = false;
789  Double_t fmin = 0,fedm = 0,ferrdef = 0;
790  Int_t npari =0, nparx = 0, istat = 0;
791  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
792 
793  Double_t corr_mat_test[3][3];
794  gMinuit->mnemat(*corr_mat_test,3);
795 
796  //Turn this test error matix into a correlation matrix
797  double corr_01 = corr_mat_test[0][1]/(corr_mat_test[0][0] *
798  corr_mat_test[1][1]);
799  double corr_02= corr_mat_test[0][2]/(corr_mat_test[0][0] *
800  corr_mat_test[2][2]);
801  double corr_12= corr_mat_test[1][2]/(corr_mat_test[1][1] *
802  corr_mat_test[2][2]);
803  if(istat == 2 || vstart[0] < 0 || vstart[1] < 0 || vstart[2] < 0 ||
804  fabs(corr_01) > 0.96 || fabs(corr_02) > 0.96 || fabs(corr_12) > 0.96)
805  troubleFitting=true;
806 
807  if(troubleFitting){
808  //There Was trouble with the original fits,
809  //most likely this was due to extremely high correlation (approx 1)
810  //between two of the fitting parameters
811  //Switch to a fit that adjusts all backgrounds
812  TMinuit *jMinuit = new TMinuit(2);
813  jMinuit->SetPrintLevel(-1);
814  jMinuit->SetFCN(fcn2Var);
815  arglist[0] = 1;
816  jMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
817 
818  jMinuit->mnparm(0, "Background", vstart[0],
819  vstep[0], 0, 0, ierflg);
820  jMinuit->mnparm(1, "Nue CC Scaling", vstart[1],
821  vstep[1], 0, 0, ierflg);
822 
823  arglist[0] = 0; //number of iterations 0 == no limit
824  arglist[1] = 1;
825 
826  jMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
827  jMinuit->mnexcm("MIGRAD", arglist, 3, ierflg);
828  jMinuit->mnexcm("SEEk", arglist, 3, ierflg);
829  jMinuit->mnexcm("MINOS", arglist, 3, ierflg);
830 
831 
832  jMinuit->GetParameter(0,fParNewStart,fParNewErr);
833  vstart[0] = fParNewStart;
834  verror[0] = fParNewErr;
835  jMinuit->GetParameter(0,fParNewStart,fParNewErr);
836  vstart[2] = fParNewStart;
837  verror[2] = fParNewErr;
838  jMinuit->GetParameter(1,fParNewStart,fParNewErr);
839  vstart[1] = fParNewStart;
840  verror[1] = fParNewErr;
841  //Get the covariance matrix
842  Double_t emat[2][2];
843  jMinuit->mnemat(*emat,2);
844 
845  nue_wei = std::make_pair (vstart[1],verror[1]);
846  numu_wei = std::make_pair (vstart[0],verror[0]);
847  nc_wei = std::make_pair (vstart[2],verror[2]);
848  covariance[0][0] = emat[0][0];
849  covariance[0][1] = emat[0][1];
850  covariance[0][2] = 0;
851  covariance[1][0] = emat[1][0];
852  covariance[1][1] = emat[1][1];
853  covariance[1][2] = 0;
854  covariance[2][0] = emat[0][0];
855  covariance[2][1] = emat[0][1];
856  covariance[2][2] = emat[0][0];
857  }
858  else{
859  gMinuit->GetParameter(0,fParNewStart,fParNewErr);
860  vstart[0] = fParNewStart;
861  verror[0] = fParNewErr;
862  gMinuit->GetParameter(2,fParNewStart,fParNewErr);
863  vstart[2] = fParNewStart;
864  verror[2] = fParNewErr;
865  gMinuit->GetParameter(1,fParNewStart,fParNewErr);
866  vstart[1] = fParNewStart;
867  verror[1] = fParNewErr;
868 
869  //Get the covariance matrix
870  Double_t emat[3][3];
871  gMinuit->mnemat(*emat,3);
872 
873  nue_wei = std::make_pair (vstart[1],verror[1]);
874  numu_wei = std::make_pair (vstart[0],verror[0]);
875  nc_wei = std::make_pair (vstart[2],verror[2]);
876  covariance[0][0] = emat[0][0];
877  covariance[0][1] = emat[0][1];
878  covariance[0][2] = emat[0][2];
879  covariance[1][0] = emat[1][0];
880  covariance[1][1] = emat[1][1];
881  covariance[1][2] = emat[1][2];
882  covariance[2][0] = emat[2][0];
883  covariance[2][1] = emat[2][1];
884  covariance[2][2] = emat[2][2];
885  }
886  }
887  else{
888  std::cout << "******** Skipping Fit For Bin "
889  << caption << "\t\n\n\n\n" << std::endl;
890  nue_wei = std::make_pair (1,0);
891  numu_wei = std::make_pair (1,0);
892  nc_wei = std::make_pair (1,0);
893  covariance[0][0] = 0;
894  covariance[0][1] = 0;
895  covariance[0][2] = 0;
896  covariance[1][0] = 0;
897  covariance[1][1] = 0;
898  covariance[1][2] = 0;
899  covariance[2][0] = 0;
900  covariance[2][1] = 0;
901  covariance[2][2] = 0;
902  }
903 
904  std::cout << "*********Fit for bin " << caption << " : " << std::endl;
905  std::cout << "NumuCC Scaling Factor= " << numu_wei.first << " +- " <<
906  numu_wei.second << std::endl;
907  std::cout << "NueCC Scaling Factor= " << nue_wei.first << " +- " <<
908  nue_wei.second << std::endl;
909  std::cout << "NC Scaling Factor= " << nc_wei.first << " +- " <<
910  nc_wei.second << std::endl;
911 
912 
913  //Projections For Systematic Error Band Drawing
914  std::vector<TH1F*> hSysts;
915  for(int z = 0; z < (int)systs_hists.size(); z++){
916  hSysts.push_back((TH1F*)systs_hists[z]->
917  ProjectionZ(ana::UniqueName().c_str(),
918  xbin,xbin,ybin,ybin));
919  }
920  std::vector<TH1F*> hSystsUp;
921  for(int z = 0; z < (int)systs_hists_up.size();z++){
922  hSystsUp.push_back((TH1F*)systs_hists_up[z]->
923  ProjectionZ(ana::UniqueName().c_str(),
924  xbin,xbin,
925  ybin,ybin));
926  }
927  std::vector<TH1F*> hSystsDown;
928  for(int z = 0; z < (int)systs_hists_down.size();z++){
929  hSystsDown.push_back((TH1F*)systs_hists_down[z]->
930  ProjectionZ(ana::UniqueName().c_str(),
931  xbin,xbin,
932  ybin,ybin));
933  }
934 
935  std::vector<TH1F*> shifts_up;
936  std::vector<TH1F*> shifts_down;
937 
938  for(int z = 0; z < (int)hSysts.size();z++){
939  shifts_up.push_back(hSysts[z]);
940  shifts_down.push_back((TH1F*)hMCZ[0]->Clone(ana::UniqueName().c_str()));
941  }
942 
943  for(int z = 0; z < (int)hSystsUp.size(); z++){
944  shifts_up.push_back(hSystsUp[z]);
945  shifts_down.push_back(hSystsDown[z]);
946  }
947 
948  //Get Systematic Error Band
949  TGraphAsymmErrors* g2 = PlotSystErrorBand(hMCZ[0],
950  shifts_up,
951  shifts_down,
952  -1,-1,1.3,false);
953 
954  //Make Copies of Nominal MC template predictions
955  std::vector<TH1F*> hMCnom;
956  for(int z = 0; z < (int)hMCZ.size(); z++){
957  hMCnom.push_back((TH1F*)hMCZ[z]->Clone(ana::UniqueName().c_str()));
958  }
959 
960  //Scale MC templates by fitted weights
961  hMCZ[1]->Scale(nue_wei.first); //Nue
962  hMCZ[2]->Scale(nue_wei.first); //Nuebar
963  hMCZ[3]->Scale(numu_wei.first); //Numu
964  hMCZ[4]->Scale(numu_wei.first); //Numubar
965  hMCZ[5]->Scale(nc_wei.first); //NC
966  hMCZ[6]->Scale(1); //Other
967  hMCZ[7]->Scale(nue_wei.first); //Non-fiducial Nuebar
968 
969 
970  //Plot Template Before after Fit
971  PlotFitResults(hDataZ,hMCnom,hMCZ,g2, pidName,dataName,caption,xbin,ybin);
972 
973  //Apply Fit Results to Analysis Histogram
974  for(int numChn = 0; numChn < (int)hResults.size(); numChn++){
975  for(int zbin = 0; zbin < hResults[0]->GetZaxis()->GetNbins();zbin++){
976  float nnue = hResults[1]->GetBinContent(xbin,ybin,zbin);
977  float nnuebar = hResults[2]->GetBinContent(xbin,ybin,zbin);
978  float nnumu = hResults[3]->GetBinContent(xbin,ybin,zbin);
979  float nnumubar = hResults[4]->GetBinContent(xbin,ybin,zbin);
980  float nnc = hResults[5]->GetBinContent(xbin,ybin,zbin);
981  float nother = hResults[6]->GetBinContent(xbin,ybin,zbin);
982  float nnf = hResults[7]->GetBinContent(xbin,ybin,zbin);
983  if(numChn == 0){
984  float binCont = nnue * nue_wei.first;
985  binCont += nnuebar * nue_wei.first;
986  binCont += nnumu * numu_wei.first;
987  binCont += nnumubar * numu_wei.first;
988  binCont += nnc * nc_wei.first;
989  binCont += nother;
990  binCont += nnf * nue_wei.first;
991  hResults[0]->SetBinContent(xbin,ybin,zbin,binCont);
992 
993  float binErr2 =
994  pow(covariance[0][0],2)*pow((nnumu+nnumubar),2) +
995  pow(covariance[1][1],2)* pow((nnue+nnuebar+nnf),2) +
996  pow(covariance[2][2],2)* pow((nnc),2) +
997  2 * covariance[0][1] * (nnumu+nnumubar)*(nnue+nnuebar+nnf) +
998  2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) +
999  2 * covariance[1][2] * (nnue + nnuebar + nnf) * (nnc) +
1000  nnue * pow(nue_wei.first,2) + nnuebar * pow(nue_wei.first,2) +
1001  nnf * pow(nue_wei.first,2) + nnumu * pow(numu_wei.first,2) +
1002  nnumubar * pow(numu_wei.first,2) + nnc * pow(nc_wei.first,2) +
1003  nother;
1004 
1005  if(covariance[0][1] == 0 && covariance[0][2] == 0 &&
1006  covariance[1][2] == 0 && nue_wei.second == 0 &&
1007  numu_wei.second == 0 && nc_wei.second == 0){
1008  binErr2 = nnue + nnuebar + nnf + nnumu + nnumubar +
1009  nnc + nother;
1010  }
1011 
1012 
1013  //if(binErr2 < 0) binErr2 = fabs(binErr2);
1014 
1015  hResults[0]->SetBinError(xbin,ybin,zbin,sqrt(binErr2));
1016 
1017  std::cout << zbin << "\t" << binCont << "\t" <<
1018  sqrt(binErr2) << std::endl;
1019  if(std::isnan(sqrt((binErr2))) || binErr2 > 1e6){
1020  std::cout << pow(covariance[0][0],2)*pow((nnumu+nnumubar),2) <<
1021  "\t" << pow(covariance[1][1],2)* pow((nnue+nnuebar+nnf),2) <<
1022  "\t" << pow(covariance[2][2],2)* pow((nnc),2) << "\t" <<
1023  2 * covariance[0][1] * (nnumu+nnumubar)*(nnue+nnuebar+nnf) <<
1024  "\t" << 2 * covariance[0][2] * (nnumu+nnumubar)*(nnc) <<
1025  "\t" << 2 * covariance[1][2] * (nnue + nnuebar + nnf) * (nnc) <<
1026  "\t" << nnue * pow(nue_wei.first,2) << "\t" <<
1027  nnuebar * pow(nue_wei.first,2) << "\t" <<
1028  nnf * pow(nue_wei.first,2) << "\t" <<
1029  nnumu * pow(numu_wei.first,2) << "\t" <<
1030  nnumubar * pow(numu_wei.first,2) << "\t" <<
1031  nnc * pow(nc_wei.first,2) << "\t" << nother << std::endl;
1032  }
1033  }
1034  if(numChn == 1){
1035  float binCont = nnue * nue_wei.first;
1036  float binErr = sqrt( pow(nue_wei.second,2)*pow(nnue,2) +
1037  nnue * pow(nue_wei.first,2));
1038  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1039  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1040  }
1041  if(numChn == 2){
1042  float binCont = nnuebar * nue_wei.first;
1043  float binErr = sqrt( pow(nue_wei.second,2)*pow(nnuebar,2) +
1044  nnuebar * pow(nue_wei.first,2) +
1045  pow(nue_wei.second,2)*pow(nnue,2) + // Add uncertainty in ws to signal component
1046  nnue * pow(nue_wei.first,2)); //
1047  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1048  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1049  }
1050  if(numChn == 3){
1051  float binCont = nnumu * numu_wei.first;
1052  float binErr = sqrt( pow(numu_wei.second,2)*pow(nnumu,2) +
1053  nnumu * pow(numu_wei.first,2));
1054  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1055  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1056  }
1057  if(numChn == 4){
1058  float binCont = nnumubar * numu_wei.first;
1059  float binErr = sqrt( pow(numu_wei.second,2)*pow(nnumubar,2) +
1060  nnumubar * pow(numu_wei.first,2));
1061  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1062  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1063  }
1064  if(numChn == 5){
1065  float binCont = nnc * nc_wei.first;
1066  float binErr = sqrt( pow(nc_wei.second,2)*pow(nnc,2) +
1067  nnc * pow(nc_wei.first,2));
1068  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1069  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1070  }
1071  if(numChn == 7){
1072  float binCont = nnf * nue_wei.first;
1073  float binErr = sqrt( pow(nue_wei.second,2)*pow(nnf,2) +
1074  nnf * pow(nue_wei.first,2));
1075  hResults[numChn]->SetBinContent(xbin,ybin,zbin,binCont);
1076  hResults[numChn]->SetBinError(xbin,ybin,zbin,binErr);
1077  }
1078  }//loop over zbins
1079  } // loop over interaction channels
1080  }//loop over ybins
1081  }//loop over xbins
1082 
1083  return hResults;
1084 }
1085 
void Simulation()
Definition: tools.h:16
T max(const caf::Proxy< T > &a, T b)
fvar< T > fmin(const fvar< T > &x1, const fvar< T > &x2)
Definition: fmin.hpp:14
TSpline3 lo("lo", xlo, ylo, 12,"0")
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const bool shouldRebin
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
correl_xv GetYaxis() -> SetDecimals()
T sqrt(T number)
Definition: d0nt_math.hpp:156
double stat_covMatrix[pidBins][pidBins]
constexpr T pow(T x)
Definition: pow.h:75
unsigned int nnue
Definition: runWimpSim.h:90
Int_t par
Definition: SimpleIterate.C:24
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
OStream cerr
Definition: OStream.cxx:7
void abs(TH1 *hist)
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
float fSigLike[pidBins]
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const Binning pidBins
Definition: NumuCC_CPiBin.h:5
double covMatrix[pidBins][pidBins]
const XML_Char const XML_Char * data
Definition: expat.h:268
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
c2
Definition: demo5.py:33
TSpline3 hi("hi", xhi, yhi, 18,"0")
const XML_Char int const XML_Char * value
Definition: expat.h:331
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
std::string getCaption2D(TH3F *sample, int xbin, int ybin)
Int_t col[ntarg]
Definition: Style.C:29
const std::string cv[Ncv]
float fOther[pidBins]
void CalculateCovarianceMatrix(TH1F *nom, std::vector< TH1F * > hSysts, int nPseudo, std::string pidName, int xbin, int ybin)
const double j
Definition: BetheBloch.cxx:29
const float s_b_ratio
const std::vector< int > kNueCCColorDef
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
const float signalRegion
float bin[41]
Definition: plottest35.C:14
TLatex * prelim
Definition: Xsec_final.C:133
z
Definition: test.py:28
OStream cout
Definition: OStream.cxx:6
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)
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)
void PrintCaption(TString str)
float fDA[pidBins]
float fNC[pidBins]
enum BeamMode nc
void fcn2Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
c1
Definition: demo5.py:24
std::string dataName
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
TPad * pad2
Definition: analysis.C:13
float fNumuLike[pidBins]
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode kGreen
enum BeamMode kBlue
Float_t w
Definition: plot.C:20
void CalculateTotalCovariance(TH1F *nom, std::vector< TH1F * > hSysts, std::string pidName, int xbin, int ybin)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
TGraphAsymmErrors * PlotSystErrorBand(TH1F *&nom, std::vector< TH1F * > &ups, std::vector< TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)
TPad * pad1
Definition: analysis.C:13
enum BeamMode string
unsigned int uint