MakePlots.C
Go to the documentation of this file.
1 //Script that does background estimation in pi^{0} kinetic energy bins
2 //It reads the spectra generated from the getSpectra_ForFitting.C macro
3 
5 #include "CAFAna/Core/Spectrum.h"
6 #include "CAFAna/Core/ISyst.h"
10 
11 #include "3FlavorAna/Cuts/NueCutsSecondAna.h"
14 
21 
22 
27 
28 #include "Utilities/rootlogon.C"
29 
30 #include "TH1.h"
31 #include "TH2.h"
32 #include "THStack.h"
33 #include "TCanvas.h"
34 #include "TLegend.h"
35 #include "TStyle.h"
36 #include "TFile.h"
37 #include <vector>
38 #include <math.h>
39 #include "TLine.h"
40 
41 #include <iostream>
42 #include <fstream>
43 #include <iomanip>
44 #include <vector>
45 #include <cmath>
46 #include <math.h>
47 #include <string>
48 #include <stdint.h>
49 #include "TFile.h"
50 #include "TTree.h"
51 #include "TCanvas.h"
52 #include "TF1.h"
53 #include "TH1F.h"
54 #include "TH2F.h"
55 #include <sstream>
56 #include <math.h>
57 #include <string>
58 #include <vector>
59 #include "TF2.h"
60 #include "TH2.h"
61 #include "TCutG.h"
62 #include "TMath.h"
63 #include "TCanvas.h"
64 #include "TStyle.h"
65 #include "TRandom.h"
66 #include "TGraph.h"
67 #include "TColor.h"
68 #include "TF3.h"
69 #include "TH3.h"
70 #include "TH3F.h"
71 #include "TMinuit.h"
72 #include "TLegend.h"
73 #include "THStack.h"
74 #include "TMatrix.h"
75 #include "TGraphErrors.h"
76 #include "TGraph.h"
77 //#include "TDecompLU.h"
78 #include "TDecompSVD.h"
79 
80 using namespace ana;
81 
82 const int xbins = 20;
83 const double pot = 8.09e20; //data pot
84 
85 float fDA[xbins];
86 float fSig[xbins];
87 float fCCBkg[xbins];
88 float fNCBkg[xbins];
89 
90 float eDA[xbins];
91 float eSig[xbins];
92 float eCCBkg[xbins];
93 float eNCBkg[xbins];
94 
97 
98 const std::vector<int> kNueCCColorDef = {
99  632, //kRed, Total MC
100  600, //kBlue Signal
101  616, //kMagenta nuebar CC
102  416-3, //kGrenn - 3 NumuCC
103 };
104 
105 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){
106 
107  double chisquare = 0;
108  float data = 0;
109  float signal = 0;
110  float ccbkg = 0;
111  float ncbkg = 0;
112 
113 
114  float data_j = 0;
115  float signal_j = 0;
116  float ccbkg_j = 0;
117  float ncbkg_j = 0;
118 
119  float error = 0;
120 
121  float delta = 0;
122  float firstterm = 0;
123  int dof_my = 10;
124 
125  // Covariance Matrix Inversion First
126  TMatrixF Vij_stat = TMatrix(xbins,xbins);
127  TMatrixF Vij = TMatrix(xbins,xbins);
128  for(int i =0; i < xbins; i++){
129  for(int j = 0; j < xbins; j++){
130  Vij_stat[i][j] = stat_covMatrix[i][j];
131  Vij[i][j] = covMatrix[i][j];
132  }
133  }
134 
135  //inversion
136  TMatrixD H3 = Vij_stat + Vij;
137  TDecompSVD svd(H3);
138  TMatrixD VijInv = svd.Invert();
139 
140  float delta_i = 0;
141  float delta_j = 0;
142 
143  TMatrixD FitVector = TMatrixD(1,xbins);
144  TMatrixD FitVector2 = TMatrixD(xbins,1);
145 
146  //Calculate Chi Sq
147  for(int i =0; i < xbins; i++){
148  data = fDA[i];
149  signal = fSig[i];
150  ccbkg = fCCBkg[i];
151  ncbkg = fNCBkg[i];
152 
153  if(data < 20){
154  FitVector[0][i] = 0;
155  FitVector2[i][0] = 0;
156  continue;
157 
158  }
159 
160  FitVector[0][i] = ( data - (par[0]*signal + par[1]*ccbkg + par[2]*ncbkg));
161  FitVector2[i][0] = ( data -(par[0]*signal+ par[1]*ccbkg + par[2]*ncbkg));
162  }
163 
164  TMatrixD Chi = FitVector * VijInv;
165  TMatrixD Chi2 = Chi * FitVector2;
166  float secondterm = Chi2(0,0);
167 
168  chisquare = secondterm;
169  f = chisquare;
170 }
171 
172 void fcn2Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag){
173 
174  double chisquare = 0;
175  float data = 0;
176  float signal = 0;
177  float ccbkg = 0;
178  float ncbkg = 0;
179 
180 
181  float data_j = 0;
182  float signal_j = 0;
183  float ccbkg_j = 0;
184  float ncbkg_j = 0;
185 
186 
187  float error = 0;
188 
189  float delta = 0;
190  float firstterm = 0;
191 
192  // Covariance Matrix Inversion First
193  TMatrixF Vij_stat = TMatrix(xbins,xbins);
194  TMatrixF Vij = TMatrix(xbins,xbins);
195  for(int i =0; i < xbins; i++){
196  for(int j = 0; j < xbins; j++){
197  Vij_stat[i][j] = stat_covMatrix[i][j];
198  Vij[i][j] = covMatrix[i][j];
199  }
200  }
201 
202 
203  //Inversion
204  TMatrixD H3 = Vij_stat + Vij;
205  TDecompSVD svd(H3);
206  TMatrixD VijInv = svd.Invert();
207 
208 
209  float delta_i = 0;
210  float delta_j = 0;
211 
212  TMatrixD FitVector = TMatrixD(1,xbins);
213  TMatrixD FitVector2 = TMatrixD(xbins,1);
214 
215  //Calculate Chi Sq
216  for(int i =0; i < xbins; i++){
217  data = fDA[i];
218  signal = fSig[i];
219  ccbkg = fCCBkg[i];
220  ncbkg = fNCBkg[i];
221 
222  if(data < 20){
223  FitVector[0][i] = 0;
224  FitVector2[i][0] = 0;
225  continue;
226 
227  }
228 
229  FitVector[0][i] = ( data - (par[0]*signal + par[1]*(ccbkg+ncbkg)));
230  FitVector2[i][0] = ( data -(par[0]*signal+ par[1]*(ccbkg+ncbkg)));
231  //dof_my++;
232  }
233 
234  TMatrixD Chi = FitVector * VijInv;
235  TMatrixD Chi2 = Chi * FitVector2;
236  float secondterm = Chi2(0,0);
237 
238 
239  chisquare = secondterm;
240 
241  f = chisquare;
242 
243 }
244 
245 
246 void PrintChiSq(TString str)
247 {
248 
249  TLatex* prelim = new TLatex(.55, .80, str);
250  prelim->SetTextColor(kBlack);
251  prelim->SetNDC();
252  prelim->SetTextSize(0.05);
253  //prelim->SetTextAlign(32);
254  prelim->Draw();
255 }
256 void PrintChiSq2(TString str)
257 {
258 
259  TLatex* prelim = new TLatex(.20, .75, str);
260  prelim->SetTextColor(kBlack);
261  prelim->SetNDC();
262  prelim->SetTextSize(0.05);
263  //prelim->SetTextAlign(32);
264  prelim->Draw();
265 }
266 
267 TGraphAsymmErrors* PlotWithSystErrorBand_mine(TH1F*& nom,
268  std::vector<TH1F*>& ups,
269  std::vector<TH1F*>& dns,
270  int col, int errCol,
271  float headroom, bool newaxis)
272 {
273  if(col == -1){
274  col = kRed;
275  errCol = kRed - 9;
276  }
277  else if(errCol == -1) errCol = col-9; // hopefully a lighter version
278 
279  nom->SetLineColor(col);
280  nom->GetXaxis()->CenterTitle();
281  nom->GetYaxis()->CenterTitle();
282  //if(newaxis) nom->Draw("hist ]["); // Set the axes up
283 
284  double yMax = nom->GetBinContent(nom->GetMaximumBin());
285 
286  TGraphAsymmErrors* g = new TGraphAsymmErrors;
287 
288  for(int binIdx = 0; binIdx < nom->GetNbinsX()+2; ++binIdx){
289  const double y = nom->GetBinContent(binIdx);
290  g->SetPoint(binIdx, nom->GetXaxis()->GetBinCenter(binIdx), y);
291 
292  const double w = nom->GetXaxis()->GetBinWidth(binIdx);
293 
294  double errUp = 0;
295  double errDn = 0;
296 
297  for(unsigned int systIdx = 0; systIdx < ups.size(); ++systIdx){
298  double hi = ups[systIdx]->GetBinContent(binIdx)-y;
299  double lo = y-dns[systIdx]->GetBinContent(binIdx);
300 
301 
302  if(lo <= 0 && hi <= 0) std::swap(lo, hi);
303 
304  errUp += hi*hi;
305  errDn += lo*lo;
306 
307  // TODO: what happens if they're both high or both low?
308  } // end for systIdx
309 
310  g->SetPointError(binIdx, w/2, w/2, sqrt(errDn), sqrt(errUp));
311  } // end for i
312 
313  g->SetFillColor(errCol);
314 
315  for(TH1* up: ups) delete up;
316  for(TH1* dn: dns) delete dn;
317  return g;
318 }
319 
320 void MakePlots(TH2F* data,
321  std::vector<TH2F*> mc,
322  TH1F* sig_weights, TH1F* cc_weights,
323  TH1F* nc_weights,
324  std::vector<TH2F*>shifts_up,
325  std::vector<TH2F*>shifts_down,
326  std::string dataname,
327  std::string weightname)
328 {
329  //PID Axis Plots
330  for(int x = 1; x <= data->GetYaxis()->GetNbins(); x++){
331 
332  std::stringstream low;
333  low<<std::fixed<<std::setprecision(2)<<
334  data->GetYaxis()->GetBinLowEdge(x);
335  std::stringstream high;
336  high<<std::fixed<<std::setprecision(2)<<
337  data->GetYaxis()->GetBinUpEdge(x);
338  std::string caption = low.str() + " #leq KE_{#pi^{0}} < " + high.str();
339 
340  TH1F* dataproj = (TH1F*)data->ProjectionX("dataproj",x,x);
341  dataproj->Rebin(2);
342  std::vector<TH1F*> shift_up_proj;
343  std::vector<TH1F*> shift_down_proj;
344  std::vector<TH1F*> shift_up_proj_ratio;
345  std::vector<TH1F*> shift_down_proj_ratio;
346 
347 
348 
349  for(uint j = 0; j < shifts_up.size(); j++){
350  char name_up[50];
351  char name_down[50];
352  sprintf(name_up, "%s_%s_%i", "proj", "up", j);
353  sprintf(name_down, "%s_%s_%i", "proj", "down", j);
354  shift_up_proj.push_back((TH1F*)shifts_up[j]->ProjectionX(name_up,x,x));
355  shift_down_proj.push_back((TH1F*)shifts_down[j]->ProjectionX(name_down,x,x));
356  sprintf(name_up, "%s_%s_%i", "proj", "up_ratio", j);
357  sprintf(name_down, "%s_%s_%i", "proj", "down_ratio", j);
358  shift_up_proj_ratio.push_back((TH1F*)shifts_up[j]->ProjectionX(name_up,x,x));
359  shift_down_proj_ratio.push_back((TH1F*)shifts_down[j]->ProjectionX(name_down,x,x));
360 
361  shift_up_proj[j]->Rebin(2);
362  shift_down_proj[j]->Rebin(2);
363  shift_up_proj_ratio[j]->Rebin(2);
364  shift_down_proj_ratio[j]->Rebin(2);
365  }
366 
367  std::vector<TH1F*> hMC;
368  std::vector<TH1F*> hMC_wgt;
369  for(uint i = 0; i < mc.size(); i++){
370  char name[50];
371  char name_wgt[50];
372 
373  sprintf(name, "%s_%i", "nominal", i);
374  sprintf(name_wgt, "%s_%i", "weight", i);
375 
376  hMC.push_back((TH1F*)mc[i]->ProjectionX(name, x,x));
377  hMC_wgt.push_back((TH1F*)mc[i]->ProjectionX(name_wgt, x,x));
378  hMC[i]->Rebin(2);
379  hMC_wgt[i]->Rebin(2);
380  }
381 
382  hMC_wgt[1]->Scale(sig_weights->GetBinContent(x));
383  hMC_wgt[6]->Scale(cc_weights->GetBinContent(x));
384  hMC_wgt[7]->Scale(nc_weights->GetBinContent(x));
385 
386  TH1F *hTotalMC_wgt = (TH1F*)hMC_wgt[1]->Clone("hTotalMC_wgt");
387  hTotalMC_wgt->Add(hMC_wgt[6]);
388  hTotalMC_wgt->Add(hMC_wgt[7]);
389  hTotalMC_wgt->SetLineColor(kBlue);
390 
391 
392  TH1F* hBkgNom = (TH1F*)hMC[6]->Clone("hBkgNom");
393  hBkgNom->Add(hMC[7]);
394  TH1F* hBkgScl = (TH1F*)hMC_wgt[6]->Clone("hBkgScl");
395  hBkgScl->Add(hMC_wgt[7]);
396 
397  TH1F* hSigNom = (TH1F*)hMC[1]->Clone("hSigNom");
398  TH1F* hSigScl = (TH1F*)hMC_wgt[1]->Clone("hSigScl");
399 
400  hSigNom->SetLineColor(30); // sig before scaling
401  hBkgNom->SetLineColor(kOrange-3); // bkg before scaling
402 
403  hSigScl->SetLineColor(3); // sig after scaling
404  hBkgScl->SetLineColor(kOrange+8); // bkg after scaling
405 
406 
407  TH1F* hDataRatio = (TH1F*) hMC[0]->Clone("data_ratio");
408  TH1F* hDataRatio_Wgt = (TH1F*)hTotalMC_wgt->Clone("data_ratio_wgt");
409 
410  hDataRatio_Wgt->SetTitle("");
411  hDataRatio->GetYaxis()->SetRangeUser(0.1,1.9);
412 
413  hDataRatio->Divide(dataproj);
414  hDataRatio_Wgt->Divide(dataproj);
415 
416  TH1F* hSignalRatio = (TH1F*)hSigScl->Clone("hSignalRatio");
417  hSignalRatio->Divide(hSigNom);
418 
419  TH1F* hBkgRatio = (TH1F*)hBkgScl->Clone("hBkgRatio");
420  hBkgRatio->Divide(hBkgNom);
421 
422  for(uint ibin = 0; ibin < shift_up_proj_ratio.size(); ibin++){
423  shift_up_proj_ratio[ibin]->Divide(dataproj);
424  shift_down_proj_ratio[ibin]->Divide(dataproj);
425  }
426 
427 
428  float chisquared = 0;
429 
430  for(int i = 1; i <= dataproj->GetXaxis()->GetNbins(); i++){
431  float d_i = dataproj->GetBinContent(i);
432  // float m_i = hTotalMC_wgt->GetBinContent(i);
433  float m_i = hTotalMC_wgt->GetBinContent(i);
434  float error = dataproj->GetBinError(i)*dataproj->GetBinError(i) +
435  hTotalMC_wgt->GetBinError(i)*hTotalMC_wgt->GetBinError(i);
436 
437  if(error == 0) continue;
438 
439  chisquared += (d_i - m_i) * (d_i - m_i)/error;
440 
441  }
442 
443  std::cout << "******ChiSq in bin :*****" << x << "\t" << chisquared << "\t" << "******And ChiSq/DoF is******* : " << chisquared/(dataproj->GetXaxis()->GetNbins() -3) << std::endl;
444 
445 
446  std::stringstream low_a;
447  low_a<<std::fixed<<std::setprecision(2)<<
448  chisquared/(dataproj->GetXaxis()->GetNbins() -3) ;
449  std::stringstream high_a;
450  high_a<<std::fixed<<std::setprecision(2)<<
451  dataproj->GetXaxis()->GetNbins() - 3;
452  //std::string caption = low.str() + " #leq E_{#nu} < " + high.str();
453 
454  std::string caption_a = "#chi^{2}/dof = " +
455  low_a.str() + " (" +
456  high_a.str() + ") ";
457 
458  hDataRatio->SetLineColor(kRed);
459  hDataRatio_Wgt->SetLineColor(kBlue);
460 
461  dataproj->GetYaxis()->SetTitle("Events/8.09#times10^{20}POT");
462  dataproj->SetMarkerStyle(20);
463  //
464  //TLegend *leg = new TLegend(0.3,0.4,0.89,0.89);
465  TLegend *leg = new TLegend(0.4,0.65,0.89,0.89);
466 
467  leg->SetNColumns(2);
468 
469  leg->AddEntry(hMC[0], "Nominal MC", "PL");
470  leg->AddEntry(hTotalMC_wgt, "Scaled MC", "PL");
471  leg->AddEntry(hSigNom, "Nominal Sig", "PL");
472  leg->AddEntry(hBkgNom, "Nominal Bkg", "PL");
473  leg->AddEntry(hSigScl, "Scaled Sig", "PL");
474  leg->AddEntry(hBkgScl, "Scaled Bkg", "PL");
475  leg->AddEntry(dataproj, "Fake Data", "PL");
476 
477 
478  TCanvas *cPID1 = new TCanvas("cPID1","cPID1",700,700);
479  gStyle->SetPadBorderMode(0);
480  gStyle->SetFrameBorderMode(0);
481  Float_t small = 1e-5;
482  cPID1->Divide(1,2,small,small);
483  cPID1->cd(1);
484  gPad->SetBottomMargin(small);
485  CenterTitles(dataproj);
486 
487 
488  // dataproj->GetYaxis()->SetRangeUser(0,2000);
489 
490  dataproj->GetXaxis()->SetRangeUser(-1,0.92);
491  //dataproj->GetYaxis()->SetRangeUser(0,1.3*dataproj->GetMaximum());
492  dataproj->SetTitle(caption.c_str());
493  dataproj->SetTitleOffset(0.08);
494  dataproj->GetXaxis()->SetTitle("");
495  dataproj->Draw();
496  hMC[0]->Draw("hist same");
497  hTotalMC_wgt->Draw("hist same");
498  hSigNom->Draw("hist same");
499  hBkgNom->Draw("hist same");
500  hSigScl->Draw("hist same");
501  hBkgScl->Draw("hist same");
502  TGraphAsymmErrors* asymm =
503  PlotWithSystErrorBand_mine(hMC[0],shift_up_proj,
504  shift_down_proj,
505  -1,-1,1.3,false);
506 
507  asymm->Draw("e2 same");
508  dataproj->Draw("same");
509  hMC[0]->Draw("hist same");
510  hTotalMC_wgt->Draw("hist same");
511  hSigNom->Draw("hist same");
512  hBkgNom->Draw("hist same");
513  hSigScl->Draw("hist same");
514  hBkgScl->Draw("hist same");
515  gPad->RedrawAxis();
516  leg->Draw();
517  // PrintChiSq(caption_a);
518  Simulation();
519  cPID1->cd(2);
520  gPad->SetTopMargin(small);
521  hDataRatio->GetYaxis()->SetTitle("MC/Fake Data");
522  hDataRatio->GetXaxis()->SetRangeUser(-1,0.92);
523  hDataRatio->GetYaxis()->SetRangeUser(0.1,1.9);
524  hDataRatio->GetYaxis()->CenterTitle();
525  hDataRatio->GetXaxis()->CenterTitle();
526  hDataRatio->Draw("hist");
527  TGraphAsymmErrors* asymm2 =
528  PlotWithSystErrorBand_mine(hDataRatio,shift_up_proj_ratio,
529  shift_down_proj_ratio,
530  -1,-1,1.3,false);
531  asymm2->Draw("e2 same");
532  hDataRatio->Draw("hist same");
533  hDataRatio_Wgt->Draw("hist same");
534  gPad->RedrawAxis();
535  char outname[50];
536  sprintf(outname, "%s%s_bin_%i.pdf", "PlotsFinal/", "bdt", x);
537  cPID1->SaveAs(outname);
538 
539 
540  TLegend *legA = new TLegend(0.4,0.6,0.89,0.89);
541  legA->AddEntry(hSignalRatio, "Sclaed Sig/Nom. Sig", "PL");
542  legA->AddEntry(hBkgRatio, "Sclaed Bkg/Nom. Bkg", "PL");
543 
544 
545  }//loop over analysis bins
546 
547  //Now For the Analysis Variable
548  TH1F* data_y = (TH1F*)data->ProjectionY("data_y");
549 
550  data_y->SetMarkerStyle(20);
551  data_y->GetXaxis()->SetRangeUser(0.2,1.5);
552 
553  std::vector<TH1F*> mc_y;
554  std::vector<TH1F*> mc_y_wgt;
555  for(uint i = 0; i < mc.size(); i++){
556  char name[50];
557  char name_wgt[50];
558  sprintf(name, "%s_%i", "yproj", i);
559  sprintf(name_wgt, "%s_%i", "yproj_wgt", i);
560  mc_y.push_back((TH1F*)mc[i]->ProjectionY(name));
561  mc_y_wgt.push_back((TH1F*)mc[i]->ProjectionY(name_wgt));
562  }
563 
564  std::vector<TH1F*> shift_up_y;
565  std::vector<TH1F*> shift_down_y;
566  std::vector<TH1F*> shift_up_y_ratio;
567  std::vector<TH1F*> shift_down_y_ratio;
568  for(uint i = 0; i < shifts_up.size(); i++){
569  char name[50];
570  char name_wgt[50];
571  sprintf(name, "%s_%i", "syst", i);
572  sprintf(name_wgt, "%s_%i", "syst_down", i);
573 
574  shift_up_y.push_back((TH1F*)shifts_up[i]->ProjectionY(name));
575  shift_down_y.push_back((TH1F*)shifts_up[i]->ProjectionY(name_wgt));
576 
577  sprintf(name, "%s_%i_ratio", "syst", i);
578  sprintf(name_wgt, "%s_%i_ratio", "syst_down", i);
579 
580  shift_up_y_ratio.push_back((TH1F*)shifts_up[i]->ProjectionY(name));
581  shift_down_y_ratio.push_back((TH1F*)shifts_up[i]->ProjectionY(name_wgt));
582 
583  }
584 
585  mc_y_wgt[1]->Multiply(sig_weights);
586  mc_y_wgt[6]->Multiply(cc_weights);
587  mc_y_wgt[7]->Multiply(nc_weights);
588 
589  TH1F* hMCWeightedTot = (TH1F*)mc_y_wgt[1]->Clone("hMCWeightedTotal");
590  hMCWeightedTot->Add(mc_y_wgt[6]);
591  hMCWeightedTot->Add(mc_y_wgt[7]);
592  hMCWeightedTot->SetLineColor(kBlue);
593 
594  TH1F* hDataRatio = (TH1F*)mc_y[0]->Clone("hDataRatio");
595  TH1F* hDataRatio_Wgt = (TH1F*)hMCWeightedTot->Clone("hDataRatio_Wgt");
596 
597  hDataRatio->SetLineColor(kRed);
598  hDataRatio_Wgt->SetLineColor(kBlue);
599  hDataRatio->GetYaxis()->SetRangeUser(0.5,1.5);
600 
601  hDataRatio->Divide(data_y);
602  hDataRatio_Wgt->Divide(data_y);
603 
604  for(uint i = 0; i < shifts_up.size(); i++){
605  shift_up_y_ratio[i]->Divide(data_y);
606  shift_down_y_ratio[i]->Divide(data_y);
607  }
608 
609  TLegend *leg = new TLegend(0.45,0.5,0.80,0.80);
610  leg->AddEntry(data_y, "Fake Data", "PL");
611  leg->AddEntry(mc_y[0], "Nominal MC", "PL");
612  leg->AddEntry(hMCWeightedTot, "Scaled MC", "PL");
613 
614  float chisquared = 0;
615 
616  for(int i = data_y->GetXaxis()->FindBin(0.21);
617  i <= data_y->GetXaxis()->FindBin(1.45); i++){
618  float d_i = data_y->GetBinContent(i);
619  float m_i = hMCWeightedTot->GetBinContent(i);
620  float error = data_y->GetBinError(i)*data_y->GetBinError(i) +
621  hMCWeightedTot->GetBinError(i)*hMCWeightedTot->GetBinError(i);
622 
623  if(error == 0) continue;
624 
625  chisquared += (d_i - m_i) * (d_i - m_i)/d_i;// error;
626  }
627 
628  std::stringstream low_a;
629  low_a<<std::fixed<<std::setprecision(2)<<
630  chisquared/(4) ;
631  std::stringstream high_a;
632  high_a<<std::fixed<<std::setprecision(2)<<
633  4;
634  //std::string caption = low.str() + " #leq E_{#nu} < " + high.str();
635 
636  std::string caption_a = "#chi^{2}/dof = " +
637  low_a.str() + " (" +
638  high_a.str() + ") ";
639 
640  std::cout << caption_a << std::endl;
641 
642  TCanvas *cPID2 = new TCanvas("cPID2","cPID2");
643  cPID2->SetBottomMargin(0.15);
644  gStyle->SetPadBorderMode(0);
645  gStyle->SetFrameBorderMode(0);
646  Float_t small = 1e-5;
647  cPID2->Divide(1,2,small,small);
648  cPID2->cd(1);
649  // data_y->GetYaxis()->SetRangeUser(0,1.4*data_y->GetMaximum());
650  data_y->GetYaxis()->SetTitle("Events/8.09#times10^{20}POT");
651  data_y->GetYaxis()->SetRangeUser(0,50000);
652  data_y->GetXaxis()->SetRangeUser(0.2,1.5);
653  data_y->GetYaxis()->CenterTitle();
654  data_y->GetXaxis()->SetTitle("");
655  data_y->GetXaxis()->CenterTitle();
656  data_y->Draw();
657  TGraphAsymmErrors* asymm3 =
658  PlotWithSystErrorBand_mine(mc_y[0],shift_up_y,shift_down_y,
659  -1,-1,1.3,false);
660  asymm3->Draw("e2 same");
661  data_y->Draw("same");
662  mc_y[0]->Draw("hist same");
663  hMCWeightedTot->Draw("hist same");
664  gPad->RedrawAxis();
665  Simulation();
666  leg->Draw();
667  gPad->RedrawAxis();
668  cPID2->cd(2);
669  gPad->SetTopMargin(small);
670  hDataRatio->GetYaxis()->SetTitle("Fake Data/MC");
671  hDataRatio->GetYaxis()->SetRangeUser(0.4,1.6);
672  hDataRatio->GetYaxis()->CenterTitle();
673  hDataRatio->GetXaxis()->CenterTitle();
674  hDataRatio->Draw("hist");
675  TGraphAsymmErrors* asymm4 =
676  PlotWithSystErrorBand_mine(hDataRatio,shift_up_y_ratio,
677  shift_down_y_ratio,
678  -1,-1,1.3,false);
679  asymm4->Draw("e2 same");
680  hDataRatio->Draw("hist same");
681  hDataRatio_Wgt->Draw("hist same");
682  gPad->RedrawAxis();
683  cPID2->SaveAs("PlotsFinal/RecoKE.pdf");
684 
685 }
686 
687 void CompareToTrueSignal(TH2F* data, TH2F* data_sig, std::vector<TH2F*> mc,
688  TH1F* sig_weights, TH1F* cc_weights,
689  TH1F* nc_weights,std::string dataname,
690  std::string weightname)
691 {
692  TH1F* data_y_signal = (TH1F*)data_sig->ProjectionY("data_y_signal");
693 
694  std::vector<TH1F*> hMC_y;
695  std::vector<TH1F*> hMC_y_wgt;
696  for(uint i = 0; i < mc.size(); i++){
697  char name[50];
698  sprintf(name, "mc_%i", i);
699  hMC_y.push_back((TH1F*)mc[i]->ProjectionY(name));
700  char name_wgt[50];
701  sprintf(name_wgt, "mc_wgt_%i", i);
702  hMC_y_wgt.push_back((TH1F*)mc[i]->ProjectionY(name_wgt));
703 
704  }
705 
706  hMC_y_wgt[1]->Multiply(sig_weights);
707  hMC_y_wgt[6]->Multiply(cc_weights);
708  hMC_y_wgt[7]->Multiply(nc_weights);
709 
710  for(int i = 1; i < hMC_y_wgt[1]->GetXaxis()->GetNbins(); i++){
711 
712  float fSignal = hMC_y[1]->GetBinContent(i);
713  float fSignalError = hMC_y[1]->GetBinError(i);
714  float fTotal = hMC_y[0]->GetBinContent(i) - fSignal;
715  float fWeight = sig_weights->GetBinContent(i);
716  float fError = sig_weights->GetBinError(i);
717 
718  float total_error = (fWeight * fSignal)* pow(pow((fSignalError/fSignal),2) + pow(fError/fWeight,2),0.5);
719 
720  hMC_y_wgt[1]->SetBinError(i,total_error);
721  }
722 
723 
724  //Data True Signal-Like
725  data_y_signal->SetLineColor(kOrange+7);
726  data_y_signal->SetMarkerColor(kOrange+7);
727 
728  //MC Signal Template
729  hMC_y[1]->SetLineColor(kRed);
730  hMC_y[1]->SetMarkerColor(kRed);
731 
732  //MC Signal Template
733  hMC_y_wgt[1]->SetLineColor(kBlue);
734  hMC_y_wgt[1]->SetMarkerColor(kBlue);
735  // hMC_y_wgt[1]->SetFillColor(kBlue);
736 
737  TLegend *leg = new TLegend(0.40,0.50,0.80,0.90);
738  leg->AddEntry(data_y_signal, "Fake Data True Signal", "PL");
739  leg->AddEntry(hMC_y[1], "Nominal Sig estimate", "PL");
740  leg->AddEntry(hMC_y_wgt[1], "Scaled Sig estimate", "PL");
741 
742  float chisquared = 0;
743  float max_value = 0;
744 
745  for(int i = data_y_signal->GetXaxis()->FindBin(0.21);
746  i <= data_y_signal->GetXaxis()->FindBin(1.45); i++){
747  float d_i = data_y_signal->GetBinContent(i);
748  float m_i = hMC_y_wgt[1]->GetBinContent(i);
749  float e_d_i = data_y_signal->GetBinError(i);
750  float e_m_i = hMC_y_wgt[1]->GetBinError(i);
751  if(d_i > max_value) max_value = d_i;
752  if(d_i == 0) continue;
753 
754  chisquared += (d_i - m_i) * (d_i - m_i)/d_i;// error;
755  }
756 
757  std::stringstream low_a;
758  low_a<<std::fixed<<std::setprecision(2)<<
759  chisquared/(4) ;
760  std::stringstream high_a;
761  high_a<<std::fixed<<std::setprecision(2)<<
762  4;
763 
764  std::string caption_a = "#chi^{2}/dof = " +
765  low_a.str() + " (" +
766  high_a.str() + ") ";
767 
768  std::cout << caption_a << std::endl;
769 
770  TCanvas *c3 = new TCanvas("c3","c3");
771  c3->SetBottomMargin(0.15);
772  gStyle->SetPadBorderMode(0);
773  gStyle->SetFrameBorderMode(0);
774  Float_t small = 1e-5;
775  c3->Divide(1,2,small,small);
776  c3->cd(1);
777  gPad->SetBottomMargin(small);
778  data_y_signal->GetXaxis()->SetRangeUser(0.2,1.5);
779  // data_y_signal->GetYaxis()->SetRangeUser(0, 50000);
780  //data_y_signal->GetYaxis()->SetRangeUser(0, 2.0*max_value);
781  data_y_signal->GetYaxis()->SetTitle("Signal events/8.09#times10^{20}POT");
782  data_y_signal->GetYaxis()->CenterTitle();
783  data_y_signal->SetMarkerStyle(20);
784  data_y_signal->SetMarkerColor(kBlack);
785  data_y_signal->SetLineColor(kBlack);
786  data_y_signal->Draw();
787  hMC_y_wgt[1]->Draw("hist same");
788  hMC_y[1]->Draw("hist same");
789  data_y_signal->Draw("same");
790  gPad->RedrawAxis();
791  Simulation();
792  leg->Draw();
793  // PrintChiSq2(caption_a);
794  c3->cd(2);
795  gPad->SetTopMargin(small);
796  TH1F* hRatio_mc = (TH1F*)data_y_signal->Clone("hRatio_mc");
797  TH1F* hRatio_wgt = (TH1F*)data_y_signal->Clone("hRatio_wgt");
798  hRatio_mc->Divide(hMC_y[1]);
799  hRatio_wgt->Divide(hMC_y_wgt[1]);
800  hRatio_mc->SetLineColor(kRed);
801  hRatio_wgt->SetLineColor(kBlue);
802  hRatio_mc->GetYaxis()->SetRangeUser(0.8,1.2);
803  hRatio_mc->GetYaxis()->SetTitle("Data/MC");
804  hRatio_mc->GetYaxis()->CenterTitle();
805  hRatio_mc->GetXaxis()->CenterTitle();
806  hRatio_mc->Draw("hist");
807  hRatio_wgt->Draw("hist same");
808  char outname[50];
809  sprintf(outname,"%s%s_%s_%s.pdf","PlotsFinal/", "signal_prediction",
810  dataname.c_str(), weightname.c_str());
811  c3->SaveAs("PlotsFinal/SignalPrediction.pdf");
812  // c3->Close();
813  //delete c3;
814 }
815 
816 void MakePlots()
817 {
818 
819 
820  std::string outname_flux =
821  "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/PPFX_universes_bdte_energy_new_bdt_bins.root";
822  TFile *flux_file = new TFile(outname_flux.c_str(), "read");
823 
824  std::string outname_genie =
825  "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/GENIE_universes_bdte_energy_new_bdt_bins.root";
826  TFile *genie_file = new TFile(outname_genie.c_str(), "read");
827 
828 
829  std::string outname_systs =
830  "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/Systs_samples_bdte_energy_new_bdt_bins.root";
831  TFile *systs_file = new TFile(outname_systs.c_str(), "read");
832 
833  std::string outname_nominal =
834  "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/Nominal_MC_bdte_energy_new_bdt_bins.root";
835  TFile *nominal_file = new TFile(outname_nominal.c_str(),"read");
836 
837  std::string outname_data =
838  "/nova/ana/users/kalra/tmva/SystOptimisation/BkgEstimation_WithSyst/forenergy/Data_bdte_energy_new_bdt_bins.root";
839  TFile *data_file = new TFile(outname_data.c_str(),"read");
840 
841 
842 
843  ////////////////////////////////////////////////////////////////////////
844  //////////// PPFX Multiverse
845  /*
846  std::vector<std::unique_ptr<ana::Spectrum>> ppfx_spects;
847  for(int i = 0; i <= 100; i++){
848  if( i < 100){
849  char name_ppfx[50];
850  sprintf(name_ppfx, "%s_%i", "flux_universe", i);
851  ppfx_spects.push_back(Spectrum::LoadFrom(flux_file->GetDirectory(name_ppfx)));
852  }
853  if( i == 100)
854  ppfx_spects.push_back(Spectrum::LoadFrom(nominal_file->GetDirectory("nominal_mc")));
855  }
856 
857  FluxMultiverseSyst fluxmulti(100,ppfx_spects,true);
858 
859  Spectrum* sPPFX_Nominal = fluxmulti.CV();
860  Spectrum* sPPFX_Upper = fluxmulti.UpperSigma();
861  Spectrum* sPPFX_Lower = fluxmulti.LowerSigma();
862 
863  TH2F* hPPFX_Nominal = (TH2F*)sPPFX_Nominal->ToTH2(pot);
864  TH2F* hPPFX_Upper = (TH2F*)sPPFX_Upper->ToTH2(pot);
865  TH2F* hPPFX_Lower = (TH2F*)sPPFX_Lower->ToTH2(pot);
866 
867  TH1F* hPPFX_Upper_1d = (TH1F*)sPPFX_Upper->ToTH1(pot);
868  TH1F* hPPFX_Lower_1d = (TH1F*)sPPFX_Lower->ToTH1(pot);
869 
870  std::vector<float> histvalue_up;
871  std::vector<float> histvalue_low;
872  for(int i = 1; i < hPPFX_Upper_1d->GetXaxis()->GetNbins(); i++){
873  histvalue_up.push_back(hPPFX_Upper_1d->GetBinContent(i));
874  histvalue_low.push_back(hPPFX_Lower_1d->GetBinContent(i));
875  }
876 
877  double counter = 0;
878  for(int i = 1; i <= hPPFX_Nominal->GetXaxis()->GetNbins(); i++){
879  for(int j = 1; j <= hPPFX_Nominal->GetYaxis()->GetNbins(); j++){
880  hPPFX_Upper->SetBinContent(i,j,histvalue_up[counter]);
881  hPPFX_Lower->SetBinContent(i,j,histvalue_low[counter]);
882  counter++;
883  }
884  }
885  */
886  ////////////////////////////////////////////////////////////////////////
887  //////////// Genie Multiverse
888  std::vector<std::unique_ptr<ana::Spectrum>>genie_spects;
889  for(int i = 0; i < 100; i++){
890  char name[50];
891  sprintf(name, "%s_%il", "genie_universe",i);
892  genie_spects.push_back(Spectrum::LoadFrom(genie_file->GetDirectory(name)));
893  }
894 
895  GenieMultiverseSpectra geniemulti =
896  GenieMultiverseSpectra(100,genie_spects,true);
897 
898  const Spectrum* sGenie_Nominal = geniemulti.Nominal();
899  const Spectrum* sGenie_Upper = geniemulti.UpperSigma();
900  const Spectrum* sGenie_Lower = geniemulti.LowerSigma();
901 
902  TH2F* hGenie_Nominal = (TH2F*)sGenie_Nominal->ToTH2(pot);
903  TH2F* hGenie_Upper = (TH2F*)sGenie_Nominal->ToTH2(pot);
904  TH2F* hGenie_Lower = (TH2F*)sGenie_Nominal->ToTH2(pot);
905 
906  TH1F* hGenie_Upper_1d = (TH1F*)sGenie_Upper->ToTH1(pot);
907  TH1F* hGenie_Lower_1d = (TH1F*)sGenie_Lower->ToTH1(pot);
908 
909  std::vector<float> histvalue_up_genie;
910  std::vector<float> histvalue_low_genie;
911  for(int i = 1; i < hGenie_Upper_1d->GetXaxis()->GetNbins(); i++){
912  histvalue_up_genie.push_back(hGenie_Upper_1d->GetBinContent(i));
913  histvalue_low_genie.push_back(hGenie_Lower_1d->GetBinContent(i));
914  }
915 
916  double counter = 0;
917  for(int i = 1; i <= hGenie_Nominal->GetXaxis()->GetNbins(); i++){
918  for(int j = 1; j <= hGenie_Nominal->GetYaxis()->GetNbins(); j++){
919  hGenie_Upper->SetBinContent(i,j,histvalue_up_genie[counter]);
920  hGenie_Lower->SetBinContent(i,j,histvalue_low_genie[counter]);
921  counter++;
922  }
923  }
924 
925  ////////////////////////////////////////////////////////////////////////
926  //////////// Other Systematics
927  std::vector<ana::Spectrum> systs_spects;
928  std::vector<ana::Spectrum> systs_spects_up;
929  std::vector<ana::Spectrum> systs_spects_down;
930 
931  systs_spects_up.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("upcal_systs")));
932  systs_spects_down.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("dwncal_systs")));
933  systs_spects_up.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("lightup_systs")));
934  systs_spects_down.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("lightdwn_systs")));
935 
936  systs_spects.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("calibshape_systs")));
937  systs_spects.push_back(*Spectrum::LoadFrom(systs_file->GetDirectory("cherenkov_systs")));
938 
939  std::vector<TH2F*> systs_hists;
940  std::vector<TH2F*> systs_hists_up;
941  std::vector<TH2F*> systs_hists_down;
942 
943  for(uint i = 0; i < systs_spects.size(); i++){
944  systs_hists.push_back((TH2F*)systs_spects[i].ToTH2(pot));
945  }
946  for(uint i = 0; i < systs_spects_up.size(); i++){
947  systs_hists_up.push_back((TH2F*)systs_spects_up[i].ToTH2(pot));
948  systs_hists_down.push_back((TH2F*)systs_spects_down[i].ToTH2(pot));
949  }
950 
951  systs_hists_up.push_back(hGenie_Upper);
952  systs_hists_down.push_back(hGenie_Lower);
953 
954  ////////////////////////////////////////////////////////////////////////
955  //////////// Nominal MC
956  std::vector<ana::Spectrum> nominal_spects;
957  for(int i = 0; i < knumchns; i++){
958  char name[50];
959  sprintf(name, "%s_%s", "template", chns[i].name.c_str());
960  nominal_spects.push_back(*Spectrum::LoadFrom(nominal_file->GetDirectory(name)));
961  }
962 
963  std::vector<TH2F*> nominal_hists;
964  for(uint i = 0; i < nominal_spects.size(); i++){
965  nominal_hists.push_back((TH2F*)nominal_spects[i].ToTH2(pot));
966  }
967 
968  ////////////////////////////////////////////////////////////////////////
969  //////////// Data
970 
971  std::string dataname="data";
972  std::string weightname = "nominal";
973  // std::string weightname = "maccres";
974 
975  char grabData[50];
976  sprintf(grabData, "%s_TotMC", dataname.c_str()); // only for nominal one
977  //sprintf(grabData, "%s_%s_TotMC", dataname.c_str(),weightname.c_str()); // fake data with sighigh and other
978 
979  Spectrum data_spect = *Spectrum::LoadFrom(data_file->GetDirectory(grabData));
980  TH2F* data_hist = (TH2F*)data_spect.ToTH2(pot);
981 
982  char grabData_signal[50];
983  sprintf(grabData_signal, "%s_Sig", dataname.c_str() );
984  // sprintf(grabData_signal, "%s_%s_Sig", dataname.c_str(),weightname.c_str());
985 
986  Spectrum data_signal_spect = *Spectrum::LoadFrom(data_file->GetDirectory(grabData_signal));
987  TH2F* data_hist_signal = (TH2F*)data_signal_spect.ToTH2(pot);
988 
989 
990  ////////////////////////////////////////////////////////////////////////
991  //////////// Create Histograms to hold weights
992  TH1F* sig_weights = (TH1F*)data_hist->ProjectionY("sig_weights");
993  TH1F* cc_weights = (TH1F*)data_hist->ProjectionY("cc_weights");
994  TH1F* nc_weights = (TH1F*)data_hist->ProjectionY("nc_weights");
995 
996 
997  //////////////////////////////////////////////////////////////////////////
998  ///////////////
999  //Start Doing the Projections and Fitting
1000  ///////////////
1001  //////////////////////////////////////////////////////////////////////////
1002  for(int i = 1; i <= data_hist->GetYaxis()->GetNbins(); i++){
1003 
1004  //for(int i = data_hist->GetYaxis()->FindBin(0.0); i <= data_hist->GetYaxis()->FindBin(0.99); i++){
1005  std::stringstream low;
1006  low<<std::fixed<<std::setprecision(2)<<
1007  data_hist->GetYaxis()->GetBinLowEdge(i);
1008  std::stringstream high;
1009  high<<std::fixed<<std::setprecision(2)<<
1010  data_hist->GetYaxis()->GetBinUpEdge(i);
1011  std::string caption = low.str() + " #leq KE_{#pi^{0}} < " + high.str();
1012 
1013 
1014  /////////////////////////////////////////////////////////////////
1015  //Project Data and MC templates
1016  TH1F* hData1D = (TH1F*)data_hist->ProjectionX("hData1D",i,i);
1017  hData1D->Rebin(2);
1018 
1019  std::vector<TH1D*>hMC;
1020  for(int j = 0; j < (int)nominal_hists.size(); j++){
1021  char name[50];
1022  sprintf(name, "mc_bin_%i_sample_%i", i,j);
1023  hMC.push_back((TH1D*)nominal_hists[j]->ProjectionX(name,i,i));
1024  hMC[j]->Rebin(2);
1025  }
1026 
1027 
1028 
1029  /////////////////////////////////////////////////////////////////
1030  //Calculate Errors due to Each Systematic Sample
1031  std::vector<TH1D*>hSysts;
1032  for(int j = 0; j < (int)systs_hists.size(); j++){
1033  char name[50];
1034  sprintf(name, "syst_bin_%i_sample_%i", i,j);
1035  TH1D* Systs_holder = systs_hists[j]->ProjectionX(name,i,i);
1036  for(int bin = 1; bin <= Systs_holder->GetXaxis()->GetNbins(); bin++){
1037  double err = Systs_holder->GetBinContent(bin);
1038  double cv = hMC[0]->GetBinContent(bin);
1039  Systs_holder->SetBinContent(bin,fabs(err-cv));
1040  }
1041  hSysts.push_back(Systs_holder);
1042  hSysts.back()->Rebin(2);
1043  }
1044 
1045 
1046  for(int j = 0; j < (int) systs_hists_up.size(); j++){
1047  char name[50];
1048  sprintf(name, "syst_up_bin_%i,smaple_%i", i, j);
1049  char name2[50];
1050  sprintf(name2, "syst_down_bin_%i,smaple_%i", i, j);
1051  TH1D* Systs_holder_up = systs_hists_up[j]->ProjectionX(name,i,i);
1052  TH1D* Systs_holder_down = systs_hists_down[j]->ProjectionX(name,i,i);
1053 
1054  for(int bin =1; bin <=Systs_holder_up->GetXaxis()->GetNbins(); bin++){
1055  double hi = Systs_holder_up->GetBinContent(bin);
1056  double lo = Systs_holder_down->GetBinContent(bin);
1057  double cv = hMC[0]->GetBinContent(bin);
1058  double value = 0.5 * (fabs(cv - hi) + fabs(cv - lo));
1059  Systs_holder_up->SetBinContent(bin, value);
1060  }
1061 
1062  hSysts.push_back(Systs_holder_up);
1063  hSysts.back()->Rebin(2);
1064  }
1065 
1066  /////////////////////////////////////////////////////////////////
1067  //Calculate Systematic Covariance Matrix Using hSysts and hMC[0]
1068  for(int j = 1; j <= hSysts[0]->GetXaxis()->GetNbins(); j++){
1069  for(int k = 1; k <= hSysts[0]->GetXaxis()->GetNbins(); k++){
1070 
1071  double multiple = 0;
1072  double x_j = 0;
1073  double x_k = 0;
1074  double counter = 0;
1075  for(int iSyst = 0; iSyst < (int)hSysts.size(); iSyst++){
1076  double val_j = hMC[0]->GetBinContent(j) +
1077  hSysts[iSyst]->GetBinContent(j);
1078  double val_k = hMC[0]->GetBinContent(k) +
1079  hSysts[iSyst]->GetBinContent(k);
1080  multiple += val_j * val_k;
1081  x_j += val_j;
1082  x_k += val_k;
1083  counter++;
1084  }
1085 
1086  double covariance = multiple/counter - ((x_j/counter)*(x_k/counter));
1087 
1088  covMatrix[j-1][k-1] = covariance;
1089 
1090  }
1091  }
1092 
1093 
1094  /////////////////////////////////////////////////////////////////
1095  //Create Matrix for Statistical Uncertainty
1096  for(int j = 1; j <= hSysts[0]->GetXaxis()->GetNbins(); j++){
1097  for(int k = 1; k <=hSysts[0]->GetXaxis()->GetNbins(); k++){
1098  if(j==k){
1099  stat_covMatrix[j-1][k-1] =
1100  hData1D->GetBinError(j)*hData1D->GetBinError(j) +
1101  hMC[0]->GetBinError(k)*hMC[0]->GetBinError(k);
1102 
1103  }
1104  else{
1105  stat_covMatrix[j-1][k-1] = 0.0;
1106  }
1107  }
1108  }
1109 
1110 
1111 
1112  /////////////////////////////////////////////////////////////////
1113  //Fill arrays for minuit fit
1114  for(int x = 1; x <= hData1D->GetXaxis()->GetNbins(); x++)
1115  {
1116  float data = hData1D->GetBinContent(x);
1117  float total = hMC[0]->GetBinContent(x);
1118  float signal = hMC[1]->GetBinContent(x);
1119  float ccbkg = hMC[6]->GetBinContent(x);
1120  float ncbkg = hMC[7]->GetBinContent(x);
1121 
1122  fDA[x-1] = data;
1123  fSig[x-1] = signal;
1124  fCCBkg[x-1] = ccbkg;
1125  fNCBkg[x-1] = ncbkg;
1126 
1127  }
1128 
1129  TH2F* covariance_hist = new TH2F("covariance_hist",";NC#pi^{0}ID;NC#pi^{0}ID",
1130  xbins, 0 ,xbins, xbins, 0, xbins);
1131  TH2F* stat_covariance_hist = new TH2F("stat_covariance_hist",
1132  ";NC#pi^{0}ID;NC#pi^{0}ID",
1133  xbins, 0 ,xbins, xbins, 0, xbins);
1134 
1135  TH2F* correlation_hist = new TH2F("correlation_hist",
1136  ";NC#pi^{0}ID;NC#pi^{0}ID",
1137  xbins, 0 ,xbins, xbins, 0, xbins);
1138 
1139  covariance_hist->GetXaxis()->CenterTitle();
1140  covariance_hist->GetYaxis()->CenterTitle();
1141  covariance_hist->GetZaxis()->CenterTitle();
1142 
1143  stat_covariance_hist->GetXaxis()->CenterTitle();
1144  stat_covariance_hist->GetYaxis()->CenterTitle();
1145  stat_covariance_hist->GetZaxis()->CenterTitle();
1146 
1147  correlation_hist->GetXaxis()->CenterTitle();
1148  correlation_hist->GetYaxis()->CenterTitle();
1149  correlation_hist->GetZaxis()->CenterTitle();
1150 
1151 
1152  //Covariance Matrix Inversion First
1153  TMatrixF Covij = TMatrix(xbins,xbins);
1154  for(int i =0; i < xbins; i++){
1155  for(int j = 0; j < xbins; j++){
1156  Covij[i][j] = covMatrix[i][j];
1157  covariance_hist->SetBinContent(i+1,j+1, covMatrix[i][j]);
1158  stat_covariance_hist->SetBinContent(i+1,j+1, stat_covMatrix[i][j]);
1159 
1160  correlation_hist->SetBinContent(i+1,j+1,
1161  covMatrix[i][j]/(sqrt(covMatrix[i][i]) *
1162  sqrt(covMatrix[j][j])));
1163  }
1164  }
1165 
1166  stat_covariance_hist->Add(covariance_hist);
1167 
1168  TCanvas *c9 = new TCanvas("c9","c9");
1169  gStyle->SetTitleAlign(32);
1170  c9->SetRightMargin(0.15);
1171  covariance_hist->SetTitle(caption.c_str());
1172  covariance_hist->Draw("COLZ");
1173  Simulation();
1174  char cov_name[50];
1175  sprintf(cov_name, "%s%s_bin_%i.pdf", "PlotsFinal/", "error_matrix", i);
1176  c9->SaveAs(cov_name);
1177 
1178 
1179  TCanvas *cCor = new TCanvas("cCor","cCor");
1180  gStyle->SetTitleAlign(32);
1181  cCor->SetRightMargin(0.15);
1182  correlation_hist->GetZaxis()->SetTitle("Correlation");
1183  correlation_hist->SetTitle(caption.c_str());
1184  correlation_hist->Draw("COLZ");
1185  Simulation();
1186  char cor_name[50];
1187  sprintf(cor_name, "%s%s_bin_%i.pdf", "PlotsFinal/", "CorMatrix", i);
1188  cCor->SaveAs(cor_name);
1189 
1190 
1191  TMinuit *gMinuit = new TMinuit(3);
1192  gMinuit->SetPrintLevel(-1);
1193  gMinuit->SetFCN(fcn);
1194  Double_t arglist[4];
1195  Int_t ierflg = 0;
1196  arglist[0] = 1;
1197  gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1198 
1199  // double vstart[3] = {1.2,0.8,0.9};
1200  double vstart[3] = {1.0,1.0,1.0};
1201  double verror[3] = {0,0,0};
1202  double vstep[3] = {0.1,0.1,0.1};
1203 
1204  gMinuit->mnparm(0, "Signal scaling", vstart[0], vstep[0], 0, 0, ierflg);
1205  gMinuit->mnparm(1, "CC-Bkg Scaling", vstart[1], vstep[1], 0, 0, ierflg);
1206  gMinuit->mnparm(2, "NC-Bkg Scaling", vstart[2], vstep[2], 0, 0, ierflg);
1207 
1208  arglist[0] = 0; //number of iterations 0 == no limit
1209  arglist[1] = 0;
1210 
1211  gMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
1212 
1213  double fParNewStart;
1214  double fParamNewErr;
1215  for(int ll = 0; ll < 3; ll++){
1216  gMinuit->GetParameter(ll,fParNewStart,fParamNewErr);
1217  vstart[ll] = fParNewStart;
1218  }
1219 
1220 
1221  gMinuit->mnexcm("MIGRAD", arglist, 3, ierflg);
1222  gMinuit->mnexcm("SEEk", arglist, 3, ierflg);
1223  for(int ll = 0; ll < 3; ll++){
1224  gMinuit->GetParameter(ll,fParNewStart,fParamNewErr);
1225  vstart[ll] = fParNewStart;
1226  //vstep[ll] = 0.05;
1227  }
1228 
1229 
1230  TMinuit *hMinuit = new TMinuit(3);
1231  // hMinuit->SetPrintLevel(-1);
1232  hMinuit->SetFCN(fcn);
1233  arglist[0] = 1;
1234  hMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1235 
1236  hMinuit->mnparm(0, "Signal Scaling", vstart[0], vstep[0], 0, 0, ierflg);
1237  hMinuit->mnparm(1, "CC-Bkg Scaling", vstart[1], vstep[1], 0, 0, ierflg);
1238  hMinuit->mnparm(2, "NC-Bkg Scaling", vstart[2], vstep[2], 0, 0, ierflg);
1239 
1240 
1241  arglist[0] = 2; //1 == standard minimization strategy
1242  //2 == try to improve minimium (slower)
1243 
1244 
1245  hMinuit->mnexcm("SET STR",arglist,1,ierflg);
1246  arglist[0] = 0;
1247  hMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
1248  hMinuit->mnexcm("MIGRAD", arglist, 5, ierflg);
1249  hMinuit->mnexcm("SEEk", arglist, 4, ierflg);
1250  hMinuit->mnexcm("MINOS", arglist, 4, ierflg);
1251 
1252  for(int ll = 0; ll < 3; ll++){
1253  hMinuit->GetParameter(ll,fParNewStart,fParamNewErr);
1254  vstart[ll] = fParNewStart;
1255  verror[ll] = fParamNewErr;
1256  }
1257 
1258  bool troubleFitting = false;
1259 
1260  Double_t fmin = 0;
1261  Double_t fedm = 0;
1262  Double_t errdef = 0;
1263  Int_t npari = 0;
1264  Int_t nparx = 0;
1265  Int_t istat = 0;
1266  hMinuit->mnstat(fmin,fedm,errdef,npari,nparx,istat);
1267  if(istat == 2)
1268  troubleFitting = true;
1269 
1270  if(!troubleFitting){
1271  hMinuit->Command("SHOw FCNvalue");
1272  }
1273 
1274  if(troubleFitting){
1275  //Fix Background Weights together
1276  TMinuit *jMinuit = new TMinuit(2);
1277  jMinuit->SetPrintLevel(-1);
1278  jMinuit->SetFCN(fcn2Var);
1279  //Double_t arglist[4];
1280  //Int_t ierflg = 0;
1281  arglist[0] = 1;
1282  jMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1283 
1284  jMinuit->mnparm(0, "Signal Scaling", vstart[0], vstep[0], 0, 0, ierflg);
1285  jMinuit->mnparm(1, "TotBkg Scaling", vstart[1], vstep[1], 0, 0, ierflg);
1286 
1287  arglist[0] = 0; //number of iterations 0 == no limit
1288  arglist[1] = 0;
1289 
1290  jMinuit->mnexcm("SIMPLEX", arglist, 2, ierflg);
1291  jMinuit->mnexcm("MIGRAD", arglist, 3, ierflg);
1292  jMinuit->mnexcm("SEEk", arglist, 3, ierflg);
1293 
1294  jMinuit->GetParameter(0,fParNewStart,fParamNewErr);
1295  vstart[0] = fParNewStart;
1296  verror[0] = fParamNewErr;
1297  // vstart[2] = fParNewStart;
1298  //verror[2] = fParamNewErr;
1299  jMinuit->GetParameter(1,fParNewStart,fParamNewErr);
1300  vstart[1] = fParNewStart;
1301  verror[1] = fParamNewErr;
1302  vstart[2] = fParNewStart;
1303  verror[2] = fParamNewErr;
1304 
1305  jMinuit->Command("SHOw FCNvalue");
1306  }
1307 
1308 
1309  // int stat = gMinuit->GetStatus();
1310  //std::cout << "Error status = " << stat << std::endl;
1311 
1312  double fParamVal;
1313  double fParamErr;
1314  hMinuit->GetParameter(0,fParamVal,fParamErr);
1315  double fParamVal_2;
1316  double fParamErr_2;
1317  hMinuit->GetParameter(1,fParamVal_2,fParamErr_2);
1318  double fParamVal_3;
1319  double fParamErr_3;
1320  hMinuit->GetParameter(2,fParamVal_3,fParamErr_3);
1321 
1322  std::cout << "*********Fit for bin " << i << ": " << std::endl;
1323  std::cout << "Signal Scaling Factor= " << vstart[0] << " +- " <<
1324  verror[0] << std::endl;
1325  sig_weights->SetBinContent(i, vstart[0]);
1326  sig_weights->SetBinError(i, verror[0]);
1327 
1328  std::cout << "CC-Bkg Scaling Factor= " << vstart[1] << " +- " <<
1329  verror[1] << std::endl;
1330  cc_weights->SetBinContent(i, vstart[1]);
1331  cc_weights->SetBinError(i, verror[1]);
1332  std::cout << "NC-Bkg Factor= " << vstart[2] << " +- " <<
1333  verror[2] << std::endl;
1334  nc_weights->SetBinContent(i, vstart[2]);
1335  nc_weights->SetBinError(i, verror[2]);
1336 
1337  delete gMinuit;
1338  delete hMinuit;
1339 
1340  } //end of loop over analysis bins
1341 
1342 
1343 
1344 
1345  MakePlots(data_hist,nominal_hists,sig_weights,cc_weights,
1346  nc_weights, systs_hists_up,systs_hists_down,dataname,weightname);
1347  CompareToTrueSignal(data_hist,data_hist_signal, nominal_hists,
1348  sig_weights, cc_weights,
1349  nc_weights,dataname,weightname);
1350 
1351 
1352  TCanvas *cWeights = new TCanvas("cWeights","cWeights");
1353  sig_weights->GetYaxis()->SetTitle("Signal Weight");
1354  sig_weights->GetXaxis()->SetTitleOffset(0.75);
1355  sig_weights->GetXaxis()->SetRangeUser(0.2,1.5);
1356  sig_weights->GetYaxis()->SetRangeUser(0,1.5);
1357  sig_weights->GetXaxis()->CenterTitle();
1358  sig_weights->GetYaxis()->CenterTitle();
1359  sig_weights->Draw("e1");
1360  Simulation();
1361 
1362  char out2[50];
1363  cWeights->SaveAs("PlotsFinal/SigWts.pdf");
1364 
1365 
1366  ///////////////////////////////////////////////////////////////////////////
1367  //CC-bkg Weights
1368  TCanvas *cWeights3 = new TCanvas("cWeights3","cWeights3");
1369  cc_weights->GetYaxis()->SetTitle("CCBkg Weight");
1370  cc_weights->GetXaxis()->SetRangeUser(0.2,1.5);
1371  cc_weights->GetYaxis()->SetRangeUser(0,1.5);
1372  cc_weights->GetXaxis()->SetTitleOffset(0.75);
1373  cc_weights->GetXaxis()->CenterTitle();
1374  cc_weights->GetYaxis()->CenterTitle();
1375  cc_weights->Draw("e1");
1376  Simulation();
1377 
1378  cWeights3->SaveAs("PlotsFinal/CCBkgWts.pdf");
1379 
1380 
1381  ///////////////////////////////////////////////////////////////////////////
1382  //NC-bkg Weights
1383  TCanvas *cWeights5 = new TCanvas("cWeights5","cWeights5");
1384  nc_weights->GetYaxis()->SetTitle("NC Weight");
1385  nc_weights->GetXaxis()->SetTitleOffset(0.75);
1386  nc_weights->GetXaxis()->CenterTitle();
1387  nc_weights->GetYaxis()->CenterTitle();
1388  nc_weights->GetXaxis()->SetRangeUser(0.2,1.5);
1389  nc_weights->GetYaxis()->SetRangeUser(0,1.5);
1390  nc_weights->Draw("e1");
1391  Simulation();
1392  cWeights5->SaveAs("PlotsFinal/NCWtd.pdf");
1393 
1394 
1395 
1396  // data_file->Close();
1397  // nominal_file->Close();
1398  //flux_file->Close();
1399  //genie_file->Close();
1400  //systs_file->Close();
1401 
1402 }
1403 
void Simulation()
Definition: tools.h:16
const std::vector< int > kNueCCColorDef
Definition: MakePlots.C:98
const XML_Char * name
Definition: expat.h:151
enum BeamMode kOrange
fvar< T > fmin(const fvar< T > &x1, const fvar< T > &x2)
Definition: fmin.hpp:14
TSpline3 lo("lo", xlo, ylo, 12,"0")
float covMatrix[xbins][xbins]
Definition: MakePlots.C:95
enum BeamMode kRed
TH2 * ToTH2(double exposure, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Spectrum must be 2D to obtain TH2.
Definition: Spectrum.cxx:165
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Eigen::ArrayXd ProjectionX(const Eigen::MatrixXd &mat)
Helper for Unweighted.
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
const SelDef chns[knumchns]
double delta
Definition: runWimpSim.h:98
float fDA[xbins]
Definition: MakePlots.C:85
T sqrt(T number)
Definition: d0nt_math.hpp:156
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: MakePlots.C:105
constexpr T pow(T x)
Definition: pow.h:75
Int_t par
Definition: SimpleIterate.C:24
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
void PrintChiSq2(TString str)
Definition: MakePlots.C:256
float fCCBkg[xbins]
Definition: MakePlots.C:87
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
const int knumchns
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
float stat_covMatrix[xbins][xbins]
Definition: MakePlots.C:96
void fcn2Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Definition: MakePlots.C:172
const XML_Char const XML_Char * data
Definition: expat.h:268
void swap(art::HLTGlobalStatus &lhs, art::HLTGlobalStatus &rhs)
TSpline3 hi("hi", xhi, yhi, 18,"0")
const XML_Char int const XML_Char * value
Definition: expat.h:331
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
const int xbins
Definition: MakePlots.C:82
void CompareToTrueSignal(TH2F *data, TH2F *data_sig, std::vector< TH2F * > mc, TH1F *sig_weights, TH1F *cc_weights, TH1F *nc_weights, std::string dataname, std::string weightname)
Definition: MakePlots.C:687
Int_t col[ntarg]
Definition: Style.C:29
const std::string cv[Ncv]
float eSig[xbins]
Definition: MakePlots.C:91
const double j
Definition: BetheBloch.cxx:29
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
float bin[41]
Definition: plottest35.C:14
TLatex * prelim
Definition: Xsec_final.C:133
OStream cout
Definition: OStream.cxx:6
float fNCBkg[xbins]
Definition: MakePlots.C:88
Eigen::ArrayXd ProjectionY(const Eigen::MatrixXd &mat)
Helper for WeightingVariable.
const Spectrum * UpperSigma(BandOptions opt=kBandFromNominal) const
void PrintChiSq(TString str)
Definition: MakePlots.C:246
const Spectrum * LowerSigma(BandOptions opt=kBandFromNominal) const
float eNCBkg[xbins]
Definition: MakePlots.C:93
void MakePlots(TH2F *data, std::vector< TH2F * > mc, TH1F *sig_weights, TH1F *cc_weights, TH1F *nc_weights, std::vector< TH2F * >shifts_up, std::vector< TH2F * >shifts_down, std::string dataname, std::string weightname)
Definition: MakePlots.C:320
const double pot
Definition: MakePlots.C:83
float eCCBkg[xbins]
Definition: MakePlots.C:92
Float_t e
Definition: plot.C:35
float fSig[xbins]
Definition: MakePlots.C:86
enum BeamMode kBlue
Float_t w
Definition: plot.C:20
TGraphAsymmErrors * PlotWithSystErrorBand_mine(TH1F *&nom, std::vector< TH1F * > &ups, std::vector< TH1F * > &dns, int col, int errCol, float headroom, bool newaxis)
Definition: MakePlots.C:267
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
Definition: UtilsExt.cxx:115
float eDA[xbins]
Definition: MakePlots.C:90
const Spectrum * Nominal() const
return the nominal universe
enum BeamMode string
unsigned int uint