NueCCIncGlobalFitter.cxx
Go to the documentation of this file.
1 #include "NueCCIncGlobalFitter.h"
2 #include <functional>
3 #include <string>
4 #include <algorithm>
5 #include <map>
7 
8 namespace nueccinc_test{
9 
10  const std::string covarianceDir = "Covariance/";
11  const bool printVerbose = true;
13 
15  (TH1F* da, std::vector<TH1F*> templates, std::vector<TH3F*> analysis,
16  TH1F* cv, std::vector<TH1F*> systs_up,
17  std::vector<TH1F*> systs_dwn, std::vector<TH1F*> systs_oneside,
18  std::vector<TH1F*> multiverse)
19  : hData1D(da),
20  hTemplates1D(templates),
21  hAnalysis3D(analysis),
22  hCV1D(cv),
23  hSystsUp1D(systs_up),
24  hSystsDown1D(systs_dwn),
25  hSystsOneSided(systs_oneside),
26  hSystsMultiverse(multiverse)
27  {
28  hData1DReduced = new TH1F(ana::UniqueName().c_str(),
29  ";Reduced Template Bins;",
30  hData1D->GetXaxis()->GetNbins(),
31  hData1D->GetXaxis()->GetXmin(),
32  hData1D->GetXaxis()->GetXmax());
33  hCV1DReduced = new TH1F(ana::UniqueName().c_str(),
34  ";Reduced Template Bins;",
35  hData1D->GetXaxis()->GetNbins(),
36  hData1D->GetXaxis()->GetXmin(),
37  hData1D->GetXaxis()->GetXmax());
38  for(uint i = 0; i < hSystsUp1D.size();i++)
39  hSystsUp1DReduced.push_back(new TH1F(ana::UniqueName().c_str(),
40  ";Reduced Template Bins;",
41  hData1D->GetXaxis()->GetNbins(),
42  hData1D->GetXaxis()->GetXmin(),
43  hData1D->GetXaxis()->GetXmax()));
44  for(uint i = 0; i < hSystsDown1D.size();i++)
45  hSystsDown1DReduced.push_back(new TH1F(ana::UniqueName().c_str(),
46  ";Reduced Template Bins;",
47  hData1D->GetXaxis()->GetNbins(),
48  hData1D->GetXaxis()->GetXmin(),
49  hData1D->GetXaxis()->GetXmax()));
50  for(uint i = 0; i < hSystsOneSided.size();i++)
51  hSystsOneSidedReduced.push_back(new TH1F(ana::UniqueName().c_str(),
52  ";Reduced Template Bins;",
53  hData1D->GetXaxis()->GetNbins(),
54  hData1D->GetXaxis()->GetXmin(),
55  hData1D->GetXaxis()->GetXmax()));
56  for(uint i = 0; i < hSystsMultiverse.size();i++)
57  hSystsMultiverseReduced.push_back(new TH1F(ana::UniqueName().c_str(),
58  ";Reduced Template Bins;",
59  hData1D->GetXaxis()->GetNbins(),
60  hData1D->GetXaxis()->GetXmin(),
61  hData1D->GetXaxis()->GetXmax()));
62 
63 
64  parNumMin = new TH2F(ana::UniqueName().c_str(),
65  ";;;",
68  parNumMax = new TH2F(ana::UniqueName().c_str(),
69  ";;;",
72 
73  reducedBinsMin = new TH2F(ana::UniqueName().c_str(),
74  ";;;",
77 
78  numu_normalization = new TH2F(ana::UniqueName().c_str(),
79  ";;;",
82 
83  numu_uncertainty = new TH2F(ana::UniqueName().c_str(),
84  ";;;",
87 
88  nue_normalization = new TH2F(ana::UniqueName().c_str(),
89  ";;;",
92 
93  nue_uncertainty = new TH2F(ana::UniqueName().c_str(),
94  ";;;",
97 
98  nc_normalization = new TH2F(ana::UniqueName().c_str(),
99  ";;;",
102 
103  nc_uncertainty = new TH2F(ana::UniqueName().c_str(),
104  ";;;",
107 
108 
109 
110  for(uint i = 1; i <= TotalXBins; i++){
111  for(uint j = 1; j <= TotalYBins; j++){
112  parNumMin->SetBinContent(i,j,-1);
113  parNumMax->SetBinContent(i,j,-1);
114  reducedBinsMin->SetBinContent(i,j,-1);
115  numu_normalization->SetBinContent(i,j,1);
116  numu_uncertainty->SetBinContent(i,j,0);
117  nue_normalization->SetBinContent(i,j,1);
118  nue_uncertainty->SetBinContent(i,j,0);
119  nc_normalization->SetBinContent(i,j,1);
120  nc_uncertainty->SetBinContent(i,j,0);
121  }
122  }
123 
124  chisquared = 0;
125  }
126 
128  //////////// //////////// //////////// //////////// ////////////
129  //Private Functions/////////////////////////////////////////////////
130  //////////// //////////// //////////// //////////// ////////////
131  //Check if the fit can be performed
132  //Are Up and Down Syst Vectors the same size?
133  //Are there an equal number of bins in all histograms?
135  {
136  bool result = true;
137  int nbins = hData1D->GetXaxis()->GetNbins();
138 
139  if( (hSystsUp1D.size() != hSystsDown1D.size())) result = false;
140  for(uint i = 0; i < hTemplates1D.size(); i++){
141  if(nbins != hTemplates1D[i]->GetXaxis()->GetNbins()) result = false;
142  }
143  for(uint i = 0; i < hSystsUp1D.size(); i++){
144  if(nbins != hSystsUp1D[i]->GetXaxis()->GetNbins()) result = false;
145  }
146  for(uint i = 0; i < hSystsDown1D.size(); i++){
147  if(nbins != hSystsDown1D[i]->GetXaxis()->GetNbins()) result = false;
148  }
149  for(uint i = 0; i < hSystsOneSided.size(); i++){
150  if(nbins != hSystsOneSided[i]->GetXaxis()->GetNbins()) result = false;
151  }
152  for(uint i = 0; i < hSystsMultiverse.size(); i++){
153  if(nbins != hSystsMultiverse[i]->GetXaxis()->GetNbins()) result = false;
154  }
155  if(nbins != hCV1D->GetXaxis()->GetNbins()) result = false;
156 
157  return result;
158  }
159 
160  //Fill Vectors that hold the values in all histograms
162  {
163  for(int i = 0; i < numPIDBins; i++){
164  vDataValues.push_back((double)hData1D->GetBinContent(i+1));
165  vMCValues.push_back((double)hTemplates1D[0]->GetBinContent(i+1));
166  vSignalLike_Values.push_back
167  ((double)hTemplates1D[1]->GetBinContent(i+1)
168  +(double)hTemplates1D[2]->GetBinContent(i+1)
169  +(double)hTemplates1D[7]->GetBinContent(i+1));
170  vNumuLike_Values.push_back
171  ((double)hTemplates1D[3]->GetBinContent(i+1)
172  +(double)hTemplates1D[4]->GetBinContent(i+1));
173  vNC_Values.push_back((double)hTemplates1D[5]->GetBinContent(i+1));
174  vOther_Values.push_back((double)hTemplates1D[6]->GetBinContent(i+1));
175 
176  if(printVerbose)
177  std::cout << i << "\t" << vDataValues[i] << "\t" <<
180 
181  }//loop over template bins
182 
183  return;
184  }
185 
186  //Calculate the one sigma shifts used to produce the covariance matrix
187  //for some of the systematic shifts (Onesided and Up/Down Systs)
189  {
190  std::vector<TH1F*> hResults; // vector to hold the one sigma shifts
191 
192  //One sided systematics
193  for(int i = 0; i < (int)hSystsOneSided.size(); i++){
194  TH1F* hHolder =
195  (TH1F*)hSystsOneSided[i]->Clone(ana::UniqueName().c_str());
196 
197  hHolder->Add(hCV1D,-1);
198  hResults.push_back(hHolder);
199  }
200 
201  //Two Sided Systematics
202  for(int i = 0; i < (int)hSystsUp1D.size(); i++){
203  TH1F* hHolderUp =
204  (TH1F*)hSystsUp1D[i]->Clone(ana::UniqueName().c_str());
205  TH1F* hHolderDwn =
206  (TH1F*)hSystsDown1D[i]->Clone(ana::UniqueName().c_str());
207 
208  hHolderUp->Add(hCV1D,-1);
209  hHolderDwn->Add(hCV1D,-1);
210 
211  for(int bin = 1; bin <= hHolderUp->GetXaxis()->GetNbins();bin++){
212  if(fabs(hHolderDwn->GetBinContent(bin)) >
213  fabs(hHolderUp->GetBinContent(bin)))
214  hHolderUp->SetBinContent(bin,hHolderDwn->GetBinContent(bin));
215  }
216  hResults.push_back((TH1F*)hHolderUp->Clone(ana::UniqueName().c_str()));
217  }
218 
219  return hResults;
220  }
221 
222  //Calculate Full Covariance Matrix, without the reduced phasespace used to fit
224  {
225  //Turn systematics samples that are not generated from a multiverse into
226  // a "one-sigma" shift that we can generate a covariance matrix from
227  std::vector<TH1F*> hSystsShifts =
229 
230  //Now create objects to hold the covariance and correlation matrices
231  std::map<std::string, TH2F*> mapCovariances;
232  std::map<std::string, TH2F*> mapCorrelations;
233 
234  //Hard-coded within the Inc. NueCC analysis scripts is the order of
235  //systematic samples
236  std::vector<std::string> vSystNames =
237  { "ckv", "calib_shape", "light", "calib", "genie", "ppfx"};
238 
239  //Use these to grab the correct histograms to plot some
240  //diagnostic plots
241  hCV1D->GetYaxis()->SetRangeUser(0,hCV1D->GetMaximum()*1.5);
242  std::string binTitle = "Full Template Space";
243  int singlesided_counter = 0;
244  int doublesided_counter = 0;
245  for(uint systNum = 0; systNum < hSystsShifts.size(); systNum++){
246  //Draw the systematic shift around the CV within template space
247  if( (vSystNames[systNum] == "ckv" ||
248  vSystNames[systNum] == "calib_shape") && makePlots){
249  hSystsOneSided[singlesided_counter]->SetLineColorAlpha(kRed-9,0.5);
250  hSystsOneSided[singlesided_counter]->SetMarkerColorAlpha(kRed-9,0.5);
251  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
252  hCV1D->GetXaxis()->SetTitle("Template Bins");
253  hCV1D->GetXaxis()->CenterTitle();
254  hCV1D->GetYaxis()->CenterTitle();
255  hCV1D->SetTitle(binTitle.c_str());
256  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
257  c1->SetLeftMargin(0.15);
258  c1->SetRightMargin(0.15);
259  hCV1D->Draw();
260  hSystsOneSided[singlesided_counter]->Draw("same");
261  TLegend* leg = ana::AutoPlaceLegend(0.3,0.4);
262  leg->AddEntry(hSystsOneSided[singlesided_counter], "Sigma Shift", "fl");
263  leg->AddEntry(hCV1D, "Central Value", "fl");
264  leg->Draw();
265  c1->SaveAs(Form("%s%s_%s_%s_only.png",covarianceDir.c_str(),
266  "standard_template", vSystNames[systNum].c_str(),
267  "full_template"));
268  c1->Close();
269  delete leg;
270  delete c1;
271  singlesided_counter++;
272  }
273  else if(makePlots){
274  hSystsUp1D[doublesided_counter]->SetLineColorAlpha(kRed-9,0.5);
275  hSystsDown1D[doublesided_counter]->SetLineColorAlpha(kBlue-9,0.5);
276  hSystsUp1D[doublesided_counter]->SetMarkerColorAlpha(kRed-9,0.5);
277  hSystsDown1D[doublesided_counter]->SetMarkerColorAlpha(kBlue-9,0.5);
278  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
279  hCV1D->GetXaxis()->SetTitle("Template Bins");
280  hCV1D->GetXaxis()->CenterTitle();
281  hCV1D->GetYaxis()->CenterTitle();
282  hCV1D->SetTitle(binTitle.c_str());
283  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
284  c1->SetLeftMargin(0.15);
285  c1->SetBottomMargin(0.15);
286  hCV1D->Draw();
287  hSystsUp1D[doublesided_counter]->Draw("same");
288  hSystsDown1D[doublesided_counter]->Draw("same");
289  TLegend* leg = ana::AutoPlaceLegend(0.3,0.4);
290  leg->AddEntry(hSystsUp1D[doublesided_counter], "Sigma Up", "fl");
291  leg->AddEntry(hSystsDown1D[doublesided_counter], "Sigma Down", "fl");
292  leg->AddEntry(hCV1D, "Central Value", "fl");
293  leg->Draw();
294  c1->SaveAs(Form("%s%s_%s_%s_only.png",covarianceDir.c_str(),
295  "standard_template", vSystNames[systNum].c_str(),
296  "full_template"));
297  c1->Close();
298  delete leg;
299  delete c1;
300  doublesided_counter++;
301  }
302 
303  //Vector to hold the one-sigma shifted sample
304  std::vector<TH1F*> vPseudoExp;
305  TH1F* hClone = (TH1F*)hCV1D->Clone(ana::UniqueName().c_str());
306  double wei = 1;
307  hClone->Add(hSystsShifts[systNum],wei);
308  vPseudoExp.push_back(hClone);
309 
310 
311  //Draw the one-sigma sample
312  if(makePlots){
313  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
314  hCV1D->Draw();
315  for(uint i = 0; i < vPseudoExp.size(); i++){
316  vPseudoExp[i]->SetLineColorAlpha(kBlue-9,0.50);
317  vPseudoExp[i]->Draw("same");
318  }
319  hCV1D->Draw("hist same");
320  c1->SaveAs(Form("%s%s_%s_%s_only.png","Covariance/",
321  "standard_template_pseudo_exp",
322  vSystNames[systNum].c_str(),"full_template"));
323  c1->Close();
324  delete c1;
325  }
326 
327 
328  TH2F* hCovHolder =
329  new TH2F(ana::UniqueName().c_str(),
330  ";Template Bins;Template Bins; Covariance (Events^{2})",
333  TH2F* hCorHolder =
334  new TH2F(ana::UniqueName().c_str(),
335  ";Template Bins;Template Bins; Correlation",
338  hCovHolder->SetTitle(binTitle.c_str());
339  hCorHolder->SetTitle(binTitle.c_str());
340 
341  //Calculate Covariance Matrix For the Single Systematic
342  for(uint iu = 0; iu < vPseudoExp.size(); iu++){
343  for(int i = 0; i < numPIDBins; i++){
344  double xi=vPseudoExp[iu]->GetBinContent(i+1);
345  double ximean=hCV1D->GetBinContent(i+1);
346  for(int j = i; j < numPIDBins; j++){
347  double xj=vPseudoExp[iu]->GetBinContent(j+1);
348  double xjmean=hCV1D->GetBinContent(j+1);
349  double current_value = hCovHolder->GetBinContent(i+1,j+1);
350  double this_value = (xi-ximean)*(xj-xjmean);
351  hCovHolder->SetBinContent(i+1,j+1,current_value+this_value);
352  }
353  }
354  }
355 
356  //Only calculate with the 1 sigma error band
357  const double N=(double)2;
358  hCovHolder->Scale(1.0/(N-1.0));
359  for(int i=0; i<numPIDBins; i++){
360  for(int j=i; j<numPIDBins; j++){
361  hCovHolder->SetBinContent(j+1,i+1,
362  hCovHolder->GetBinContent(i+1,j+1));
363  }
364  }
365 
366  for(int i=0; i<numPIDBins; i++){
367  for(int j=0; j<numPIDBins; j++){
368  double covariance = hCovHolder->GetBinContent(i+1,j+1);
369  double variance_i = sqrt(hCovHolder->GetBinContent(i+1,i+1));
370  double variance_j = sqrt(hCovHolder->GetBinContent(j+1,j+1));
371  double correlation = covariance/(variance_i*variance_j);
372  if(isnan(correlation)) correlation = 0;
373  hCorHolder->SetBinContent(i+1,j+1,
374  correlation);
375  }
376  }
377 
378 
379  mapCovariances.insert(std::make_pair(vSystNames[systNum],
380  hCovHolder));
381  mapCorrelations.insert(std::make_pair(vSystNames[systNum],
382  hCorHolder));
383  }//loop over systematic samples
384 
385  //Now to calculate the covariance matrices from the multiverse samples
386  //These are calculated from shape only GENIE and PPFX Multiverses
387  {
388  TH2F* hCovHolder =
389  new TH2F(ana::UniqueName().c_str(),
390  ";Template Bins;Template Bins; Covariance (Events^{2})",
393  TH2F* hCorHolder =
394  new TH2F(ana::UniqueName().c_str(),
395  ";Template Bins;Template Bins;Correlation",
398 
399  TH2F* hCovHolderPPFX =
400  new TH2F(ana::UniqueName().c_str(),
401  ";Template Bins;Template Bins; Covariance (Events^{2})",
404  TH2F* hCorHolderPPFX =
405  new TH2F(ana::UniqueName().c_str(),
406  ";Template Bins;Template Bins;Correlation",
409 
410  hCovHolder->SetTitle(binTitle.c_str());
411  hCorHolder->SetTitle(binTitle.c_str());
412  hCovHolderPPFX->SetTitle(binTitle.c_str());
413  hCorHolderPPFX->SetTitle(binTitle.c_str());
414 
415  //Get Ready to calculate the covariance and correlation matrices
416  std::vector<TH1D*> hMultiverseHistsGENIE;
417  std::vector<TH1D*> hMultiverseHistsPPFX;
418 
419  //Performing area normalization
420  //GENIE Universes
421  for(uint i = 0; i < hSystsMultiverse.size()- 100; i++){
422  TH1D* hHolder=
423  (TH1D*)hSystsMultiverse[i]->Clone(ana::UniqueName().c_str());
424  float ymax = hHolder->Integral();
425  float cvmax = hCV1D->Integral();
426  hHolder->Scale(cvmax/ymax);
427  hMultiverseHistsGENIE.push_back(hHolder);
428  }
429 
430  //PPFX Universes
431  for(uint i = hSystsMultiverse.size()-100;
432  i < hSystsMultiverse.size(); i++){
433  TH1D* hHolder=
434  (TH1D*)hSystsMultiverse[i]->Clone(ana::UniqueName().c_str());
435  float ymax = hHolder->Integral();
436  float cvmax = hCV1D->Integral();
437  hHolder->Scale(cvmax/ymax);
438  hMultiverseHistsPPFX.push_back(hHolder);
439  }
440 
441  //Covariance and Correlation Calculations
442  //Use an external Class from trunk/CAFAna/XSec/MultiverseCorrelation.h
443  ana::MultiverseCorrelation* multi_correlation_calc_genie =
444  new ana::MultiverseCorrelation(hMultiverseHistsGENIE);
445 
446  TH2F* hCov =
447  (TH2F*)multi_correlation_calc_genie->GetCovarianceMatrix();
448  TH2F* hCor =
449  (TH2F*)multi_correlation_calc_genie->GetCorrelationMatrix();
450 
451  ana::MultiverseCorrelation* multi_correlation_calc_ppfx =
452  new ana::MultiverseCorrelation(hMultiverseHistsPPFX);
453 
454  TH2F* hCovPPFX =
455  (TH2F*)multi_correlation_calc_ppfx->GetCovarianceMatrix();
456  TH2F* hCorPPFX =
457  (TH2F*)multi_correlation_calc_ppfx->GetCorrelationMatrix();
458 
459  for(int i = 1; i <= hCov->GetXaxis()->GetNbins(); i++){
460  for(int j = 1; j <= hCov->GetXaxis()->GetNbins(); j++){
461  hCovHolder->SetBinContent(i,j,hCov->GetBinContent(i,j));
462  hCorHolder->SetBinContent(i,j,hCor->GetBinContent(i,j));
463  hCovHolderPPFX->SetBinContent(i,j,hCovPPFX->GetBinContent(i,j));
464  hCorHolderPPFX->SetBinContent(i,j,hCorPPFX->GetBinContent(i,j));
465  }
466  }
467 
468  //Make Plots
469  if(makePlots){
470  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
471  hCV1D->GetYaxis()->SetTitle("Template Bins");
472  hCV1D->GetYaxis()->CenterTitle();
473  hCV1D->GetXaxis()->CenterTitle();
474  hCV1D->SetTitle(binTitle.c_str());
475 
476  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
477  hCV1D->Draw();
478  for(uint i = 0; i < hMultiverseHistsPPFX.size(); i++){
479  hMultiverseHistsPPFX[i]->SetLineColorAlpha(kBlue-9,0.50);
480  hMultiverseHistsPPFX[i]->Draw("hist same");
481  }
482  hCV1D->Draw("same");
483  c1->SaveAs(Form("%s%s_%s_%s_only.png","Covariance/",
484  "standard_template_pseudo_exp",
485  "ppfx","full_template"));
486  c1->Close();
487  }
488 
489  if(makePlots){
490  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
491  hCV1D->GetYaxis()->SetTitle("Template Bins");
492  hCV1D->GetYaxis()->CenterTitle();
493  hCV1D->GetXaxis()->CenterTitle();
494  hCV1D->SetTitle(binTitle.c_str());
495 
496  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
497  hCV1D->Draw();
498  for(uint i = 0; i < hMultiverseHistsGENIE.size(); i++){
499  hMultiverseHistsGENIE[i]->SetLineColorAlpha(kBlue-9,0.50);
500  hMultiverseHistsGENIE[i]->Draw("hist same");
501  }
502  hCV1D->Draw("same");
503  c1->SaveAs(Form("%s%s_%s_%s_only.png","Covariance/",
504  "standard_template_pseudo_exp",
505  "genie","full_template"));
506  c1->Close();
507  }
508 
509  mapCovariances.insert(std::make_pair("genie",
510  hCovHolder));
511  mapCorrelations.insert(std::make_pair("genie",
512  hCorHolder));
513  mapCovariances.insert(std::make_pair("ppfx",
514  hCovHolderPPFX));
515  mapCorrelations.insert(std::make_pair("ppfx",
516  hCorHolderPPFX));
517  delete hCorPPFX;
518  delete hCovPPFX;
519  delete multi_correlation_calc_ppfx;
520  for(auto ele :hMultiverseHistsPPFX)
521  delete ele;
522  delete hCor;
523  delete hCov;
524  for(auto ele :hMultiverseHistsGENIE)
525  delete ele;
526  delete multi_correlation_calc_genie;
527 
528  }//multiverse sample calculation
529 
530 
531  //Now start adding up each of the individual covariance matrices
532  TH2F* hCovariance =
533  (TH2F*)mapCovariances["genie"]->Clone(ana::UniqueName().c_str());
534 
535 
536 
537 
538  std::string output = covarianceDir + "covariance_addition";
539  for(auto pair: mapCovariances){
540  if(pair.first == "genie")continue;
541  output+="_";
542  output+=pair.first;
543  hCovariance->Add(pair.second);
544  if(makePlots){
545  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
546  c3->SetLeftMargin(0.1);
547  c3->SetRightMargin(0.15);
548  c3->SetBottomMargin(0.1);
549  hCovariance->SetTitle(binTitle.c_str());
550  hCovariance->GetZaxis()->SetMaxDigits(3);
551  hCovariance->GetZaxis()->SetTitle("Covariance (Events^{2})");
552  hCovariance->Draw("COLZ");
553  c3->SaveAs(Form("%s_%s.png",output.c_str(),"full_template"));
554  c3->Close();
555  delete c3;
556  }
557  }
558 
559  TH2F* hCorrelation =
560  (TH2F*)hCovariance->Clone(ana::UniqueName().c_str());
561 
562  for(int i=0; i<numPIDBins; i++){
563  for(int j=0; j<numPIDBins; j++){
564  double covariance = hCovariance->GetBinContent(i+1,j+1);
565  double variance_i = sqrt(hCovariance->GetBinContent(i+1,i+1));
566  double variance_j = sqrt(hCovariance->GetBinContent(j+1,j+1));
567  double correlation = covariance/(variance_i*variance_j);
568  if(isnan(correlation)) correlation = 0;
569  hCorrelation->SetBinContent(i+1,j+1,
570  correlation);
571  }
572  }
573 
574  if(makePlots){
575  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
576  hCorrelation->GetXaxis()->CenterTitle();
577  hCorrelation->GetYaxis()->CenterTitle();
578  hCorrelation->GetZaxis()->CenterTitle();
579  hCorrelation->SetTitle(binTitle.c_str());
580  hCorrelation->GetZaxis()->SetTitle("Correlation");
581  hCorrelation->GetZaxis()->SetRangeUser(-1,1);
582  c3->SetLeftMargin(0.1);
583  c3->SetRightMargin(0.1);
584  c3->SetBottomMargin(0.1);
585  hCorrelation->GetYaxis()->SetMaxDigits(3);
586  hCorrelation->Draw("COLZ");
587  c3->SaveAs(Form("%s%s_%s.png",covarianceDir.c_str(),
588  "correlation_all","full_template"));
589  c3->Close();
590  delete c3;
591  }
592 
593  //Move the covariance matrix to the holder so that it can be used
594  for(uint i = 0; i < numPIDBins; i++){
595  std::vector<float> hCovCol;
596  for(uint j = 0; j < numPIDBins; j++){
597  hCovCol.push_back(hCovariance->GetBinContent(i+1,j+1));
598  }
599  vCovarianceMatrixFull.push_back(hCovCol);
600  }
601 
602  //Clean up memory a bit
603  for(uint i = 0; i < hSystsShifts.size(); i++){
604  delete hSystsShifts[i];
605  }
606  hSystsShifts.clear();
607  }
608 
609  //Fill the diagonal statistical covariance matrix from the full template
610  //space
612  {
613  for(int i = 0; i < numPIDBins; i++){
614  std::vector<float> StatCovCol;
615  for(int j = 0; j < numPIDBins; j++){
616  if(i == j){
617  StatCovCol.push_back( hCV1D->GetBinError(i+1)*
618  hCV1D->GetBinError(i+1));
619  }
620  else
621  StatCovCol.push_back(0);
622  }
623  vStatCovarianceMatrixFull.push_back(StatCovCol);
624  }
625  }
626 
627  //Get Values from the Systematic or Statistical Covariance Matrix
629  {
630  return vCovarianceMatrixFull[binx][biny];
631  }
632 
634  {
635  return vStatCovarianceMatrixFull[binx][biny];
636  }
637 
638  //After the Covariance Matrices have been filled
639  //Reduce the phase space of the fit to that used within the analysis itself
640  //And remove any additional areas of phase space that will
641  //cause the covariance matrix to be non-invertable.
643  {
644 
645  //Count the number of bins (or columns of the covariance matrix)
646  //that will be used within the analysis
647  int iCols = 0;
648 
649  //WARNING: Very hard coded!
650  for(int xcol = 0; xcol < numPIDBins; xcol++){
651  const int nmodnynz = xcol%TotalYZBins;
652  const int ix = xcol/TotalYZBins;
653  const int iy = nmodnynz/TotalZBins;
654  if((vCovarianceMatrixFull[xcol][xcol] != 0 &&
655  vStatCovarianceMatrixFull[xcol][xcol] != 0) &&
656  ( ((ix+1 == 5) && (iy+1 < 5) && (iy+1 > 2)) ||
657  ((ix+1 == 6) && (iy+1 < 7) && (iy+1 > 2)) ||
658  ((ix+1 == 7) && (iy+1 < 9) && (iy+1 > 2)) ||
659  ((ix+1 == 8) && (iy+1 < 12) && (iy+1 > 2)) ))
660  iCols++;
661  }
662 
663  //Now apply the Reduced Phase Space to the all input histograms.
664  hData1DReduced->SetBins(iCols,0,iCols);
665  hCV1DReduced->SetBins(iCols,0,iCols);
666  for(uint i = 0; i < hSystsUp1D.size();i++)
667  hSystsUp1DReduced[i]->SetBins(iCols,0,iCols);
668  for(uint i = 0; i < hSystsDown1D.size();i++)
669  hSystsDown1DReduced[i]->SetBins(iCols,0,iCols);
670  for(uint i = 0; i < hSystsOneSided.size();i++)
671  hSystsOneSidedReduced[i]->SetBins(iCols,0,iCols);
672  for(uint i = 0; i < hSystsMultiverse.size();i++)
673  hSystsMultiverseReduced[i]->SetBins(iCols,0,iCols);
674 
675 
676  //Fill the reduced phase space histograms
677  int binNum = 0;
678  for(int xcol = 0; xcol < numPIDBins; xcol++){
679  const int nmodnynz = xcol%TotalYZBins;
680  const int ix = xcol/TotalYZBins;
681  const int iy = nmodnynz/TotalZBins;
682  if((vCovarianceMatrixFull[xcol][xcol] != 0 &&
683  vStatCovarianceMatrixFull[xcol][xcol] != 0) &&
684  ( ((ix+1 == 5) && (iy+1 < 5) && (iy+1 > 2)) ||
685  ((ix+1 == 6) && (iy+1 < 7) && (iy+1 > 2)) ||
686  ((ix+1 == 7) && (iy+1 < 9) && (iy+1 > 2)) ||
687  ((ix+1 == 8) && (iy+1 < 12) && (iy+1 > 2)) ))
688  {
689  binNum++;
690  hData1DReduced->SetBinContent(binNum,
691  hData1D->GetBinContent(xcol+1));
692  hCV1DReduced->SetBinContent(binNum,
693  hCV1D->GetBinContent(xcol+1));
694  for(uint i = 0; i < hSystsUp1D.size();i++)
696  hSystsUp1D[i]->
697  GetBinContent(xcol+1));
698  for(uint i = 0; i < hSystsDown1D.size();i++)
700  hSystsDown1D[i]->
701  GetBinContent(xcol+1));
702  for(uint i = 0; i < hSystsOneSided.size();i++)
704  hSystsOneSided[i]->
705  GetBinContent(xcol+1));
706  for(uint i = 0; i < hSystsMultiverse.size();i++)
709  GetBinContent(xcol+1));
710 
711  vDataValuesReduced.push_back((double)hData1D->GetBinContent(xcol+1));
712  vMCValuesReduced.push_back((double)hTemplates1D[0]->GetBinContent(xcol+1));
713  vSignalLike_ValuesReduced.push_back
714  ((double)hTemplates1D[1]->GetBinContent(xcol+1)
715  +(double)hTemplates1D[2]->GetBinContent(xcol+1)
716  +(double)hTemplates1D[7]->GetBinContent(xcol+1));
717  vNumuLike_ValuesReduced.push_back
718  ((double)hTemplates1D[3]->GetBinContent(xcol+1)
719  +(double)hTemplates1D[4]->GetBinContent(xcol+1));
720  vNC_ValuesReduced.push_back
721  ((double)hTemplates1D[5]->GetBinContent(xcol+1));
722  vOther_ValuesReduced.push_back
723  ((double)hTemplates1D[6]->GetBinContent(xcol+1));
724 
725  int i = binNum - 1;
726  if(printVerbose)
727  std::cout << i << "\t" << vDataValuesReduced[i] << "\t" <<
731 
732  }
733  }
734 
735  return;
736  }
737 
738  //Calculate the one sigma shifts used to produce the reduced covariance matrix
739  //for some of the systematic shifts (Onesided and Up/Down Systs)
741  {
742 
743  std::vector<TH1F*> hResults; // vector to hold the one sigma shifts
744 
745  //One sided systematics
746  for(int i = 0; i < (int)hSystsOneSidedReduced.size(); i++){
747  TH1F* hHolder =
748  (TH1F*)hSystsOneSidedReduced[i]->Clone(ana::UniqueName().c_str());
749 
750  hHolder->Add(hCV1DReduced,-1);
751  hResults.push_back(hHolder);
752  }
753 
754  //Two Sided Systematics
755  for(int i = 0; i < (int)hSystsUp1DReduced.size(); i++){
756  TH1F* hHolderUp =
757  (TH1F*)hSystsUp1DReduced[i]->Clone(ana::UniqueName().c_str());
758  TH1F* hHolderDwn =
759  (TH1F*)hSystsDown1DReduced[i]->Clone(ana::UniqueName().c_str());
760 
761  hHolderUp->Add(hCV1DReduced,-1);
762  hHolderDwn->Add(hCV1DReduced,-1);
763 
764  for(int bin = 1; bin <= hHolderUp->GetXaxis()->GetNbins();bin++){
765  if(fabs(hHolderDwn->GetBinContent(bin)) >
766  fabs(hHolderUp->GetBinContent(bin)))
767  hHolderUp->SetBinContent(bin,hHolderDwn->GetBinContent(bin));
768  }
769  hResults.push_back((TH1F*)hHolderUp->Clone(ana::UniqueName().c_str()));
770  }
771 
772  return hResults;
773  }
774 
775  //Calculate Full Covariance Matrix, without the reduced phasespace used to fit
777  {
778  //Turn systematics samples that are not generated from a multiverse into
779  // a "one-sigma" shift that we can generate a covariance matrix from
780  std::vector<TH1F*> hSystsShifts =
782 
783  //Now create objects to hold the covariance and correlation matrices
784  std::map<std::string, TH2F*> mapCovariances;
785  std::map<std::string, TH2F*> mapCorrelations;
786 
787  //Hard-coded within the Inc. NueCC analysis scripts is the order of
788  //systematic samples
789  std::vector<std::string> vSystNames =
790  { "ckv", "calib_shape", "light", "calib", "genie", "ppfx"};
791 
792  //Use these to grab the correct histograms to plot some
793  //diagnostic plots
794 
795  hCV1DReduced->GetYaxis()->SetRangeUser(0,
796  hCV1DReduced->GetMaximum()*1.5);
797 
798  std::string binTitle = "Reduced Template Space";
799  int singlesided_counter = 0;
800  int doublesided_counter = 0;
801  for(uint systNum = 0; systNum < hSystsShifts.size(); systNum++){
802  //Draw the systematic shift around the CV within template space
803  if( (vSystNames[systNum] == "ckv" ||
804  vSystNames[systNum] == "calib_shape") && makePlots){
805  hSystsOneSidedReduced[singlesided_counter]->
806  SetLineColorAlpha(kRed-9,0.5);
807  hSystsOneSidedReduced[singlesided_counter]->
808  SetMarkerColorAlpha(kRed-9,0.5);
809  hCV1DReduced->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
810  hCV1DReduced->GetXaxis()->SetTitle("Template Bins");
811  hCV1DReduced->GetXaxis()->CenterTitle();
812  hCV1DReduced->GetYaxis()->CenterTitle();
813  hCV1DReduced->SetTitle(binTitle.c_str());
814  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
815  c1->SetLeftMargin(0.15);
816  c1->SetRightMargin(0.15);
817  hCV1DReduced->Draw();
818  hSystsOneSidedReduced[singlesided_counter]->Draw("same");
819  TLegend* leg = ana::AutoPlaceLegend(0.3,0.4);
820  leg->AddEntry(hSystsOneSidedReduced[singlesided_counter],
821  "Sigma Shift", "fl");
822  leg->AddEntry(hCV1DReduced, "Central Value", "fl");
823  leg->Draw();
824  c1->SaveAs(Form("%s%s_%s_%s_only.png",covarianceDir.c_str(),
825  "standard_template", vSystNames[systNum].c_str(),
826  "reduced_template"));
827  c1->Close();
828  delete leg;
829  delete c1;
830  singlesided_counter++;
831  }
832  else if(makePlots){
833  hSystsUp1DReduced[doublesided_counter]->SetLineColorAlpha(kRed-9,0.5);
834  hSystsDown1DReduced[doublesided_counter]->SetLineColorAlpha(kBlue-9,0.5);
835  hSystsUp1DReduced[doublesided_counter]->SetMarkerColorAlpha(kRed-9,0.5);
836  hSystsDown1DReduced[doublesided_counter]->SetMarkerColorAlpha(kBlue-9,0.5);
837  hCV1DReduced->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
838  hCV1DReduced->GetXaxis()->SetTitle("Template Bins");
839  hCV1DReduced->GetXaxis()->CenterTitle();
840  hCV1DReduced->GetYaxis()->CenterTitle();
841  hCV1DReduced->SetTitle(binTitle.c_str());
842  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
843  c1->SetLeftMargin(0.15);
844  c1->SetBottomMargin(0.15);
845  hCV1DReduced->Draw();
846  hSystsUp1DReduced[doublesided_counter]->Draw("same");
847  hSystsDown1DReduced[doublesided_counter]->Draw("same");
848  TLegend* leg = ana::AutoPlaceLegend(0.3,0.4);
849  leg->AddEntry(hSystsUp1DReduced[doublesided_counter], "Sigma Up", "fl");
850  leg->AddEntry(hSystsDown1DReduced[doublesided_counter],
851  "Sigma Down", "fl");
852  leg->AddEntry(hCV1DReduced, "Central Value", "fl");
853  leg->Draw();
854  c1->SaveAs(Form("%s%s_%s_%s_only.png",covarianceDir.c_str(),
855  "standard_template", vSystNames[systNum].c_str(),
856  "reduced_template"));
857  c1->Close();
858  delete leg;
859  delete c1;
860  doublesided_counter++;
861  }
862 
863  //Vector to hold the one-sigma shifted sample
864  std::vector<TH1F*> vPseudoExp;
865  TH1F* hClone = (TH1F*)hCV1DReduced->Clone(ana::UniqueName().c_str());
866  double wei = 1;
867  hClone->Add(hSystsShifts[systNum],wei);
868  vPseudoExp.push_back(hClone);
869 
870 
871  //Draw the one-sigma sample
872  if(makePlots){
873  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
874  hCV1DReduced->Draw();
875  for(uint i = 0; i < vPseudoExp.size(); i++){
876  vPseudoExp[i]->SetLineColorAlpha(kBlue-9,0.50);
877  vPseudoExp[i]->Draw("same");
878  }
879  hCV1DReduced->Draw("hist same");
880  c1->SaveAs(Form("%s%s_%s_%s_only.png","Covariance/",
881  "standard_template_pseudo_exp",
882  vSystNames[systNum].c_str(),"reduced_template"));
883  c1->Close();
884  delete c1;
885  }
886 
887 
888  TH2F* hCovHolder =
889  new TH2F(ana::UniqueName().c_str(),
890  ";Template Bins;Template Bins; Covariance (Events^{2})",
891  hData1DReduced->GetXaxis()->GetNbins(),0,
892  hData1DReduced->GetXaxis()->GetNbins(),
893  hData1DReduced->GetXaxis()->GetNbins(),0,
894  hData1DReduced->GetXaxis()->GetNbins());
895  TH2F* hCorHolder =
896  new TH2F(ana::UniqueName().c_str(),
897  ";Template Bins;Template Bins; Correlation",
898  hData1DReduced->GetXaxis()->GetNbins(),0,
899  hData1DReduced->GetXaxis()->GetNbins(),
900  hData1DReduced->GetXaxis()->GetNbins(),0,
901  hData1DReduced->GetXaxis()->GetNbins());
902  hCovHolder->SetTitle(binTitle.c_str());
903  hCorHolder->SetTitle(binTitle.c_str());
904 
905  //Calculate Covariance Matrix For the Single Systematic
906  for(uint iu = 0; iu < vPseudoExp.size(); iu++){
907  for(int i = 0; i < hCV1DReduced->GetXaxis()->GetNbins(); i++){
908  double xi=vPseudoExp[iu]->GetBinContent(i+1);
909  double ximean=hCV1DReduced->GetBinContent(i+1);
910  for(int j = i; j < hCV1DReduced->GetXaxis()->GetNbins() ; j++){
911  double xj=vPseudoExp[iu]->GetBinContent(j+1);
912  double xjmean=hCV1DReduced->GetBinContent(j+1);
913  double current_value = hCovHolder->GetBinContent(i+1,j+1);
914  double this_value = (xi-ximean)*(xj-xjmean);
915  hCovHolder->SetBinContent(i+1,j+1,current_value+this_value);
916  }
917  }
918  }
919 
920  //Only calculate with the 1 sigma error band
921  const double N=(double)2;
922  hCovHolder->Scale(1.0/(N-1.0));
923  for(int i=0; i< hCV1DReduced->GetXaxis()->GetNbins(); i++){
924  for(int j=i; j< hCV1DReduced->GetXaxis()->GetNbins(); j++){
925  hCovHolder->SetBinContent(j+1,i+1,
926  hCovHolder->GetBinContent(i+1,j+1));
927  }
928  }
929 
930  for(int i=0; i< hCV1DReduced->GetXaxis()->GetNbins(); i++){
931  for(int j=0; j< hCV1DReduced->GetXaxis()->GetNbins(); j++){
932  double covariance = hCovHolder->GetBinContent(i+1,j+1);
933  double variance_i = sqrt(hCovHolder->GetBinContent(i+1,i+1));
934  double variance_j = sqrt(hCovHolder->GetBinContent(j+1,j+1));
935  double correlation = covariance/(variance_i*variance_j);
936  if(isnan(correlation)) correlation = 0;
937  hCorHolder->SetBinContent(i+1,j+1,
938  correlation);
939  }
940  }
941 
942 
943  mapCovariances.insert(std::make_pair(vSystNames[systNum],
944  hCovHolder));
945  mapCorrelations.insert(std::make_pair(vSystNames[systNum],
946  hCorHolder));
947  }//loop over systematic samples
948 
949  //Now to calculate the covariance matrices from the multiverse samples
950  //These are calculated from shape only GENIE and PPFX Multiverses
951  {
952  TH2F* hCovHolder =
953  new TH2F(ana::UniqueName().c_str(),
954  ";Template Bins;Template Bins; Covariance (Events^{2})",
955  hData1DReduced->GetXaxis()->GetNbins(),0,
956  hData1DReduced->GetXaxis()->GetNbins(),
957  hData1DReduced->GetXaxis()->GetNbins(),0,
958  hData1DReduced->GetXaxis()->GetNbins());
959  TH2F* hCorHolder =
960  new TH2F(ana::UniqueName().c_str(),
961  ";Template Bins;Template Bins;Correlation",
962  hData1DReduced->GetXaxis()->GetNbins(),0,
963  hData1DReduced->GetXaxis()->GetNbins(),
964  hData1DReduced->GetXaxis()->GetNbins(),0,
965  hData1DReduced->GetXaxis()->GetNbins());
966 
967  TH2F* hCovHolderPPFX =
968  new TH2F(ana::UniqueName().c_str(),
969  ";Template Bins;Template Bins; Covariance (Events^{2})",
970  hData1DReduced->GetXaxis()->GetNbins(),0,
971  hData1DReduced->GetXaxis()->GetNbins(),
972  hData1DReduced->GetXaxis()->GetNbins(),0,
973  hData1DReduced->GetXaxis()->GetNbins());
974 
975  TH2F* hCorHolderPPFX =
976  new TH2F(ana::UniqueName().c_str(),
977  ";Template Bins;Template Bins;Correlation",
978  hData1DReduced->GetXaxis()->GetNbins(),0,
979  hData1DReduced->GetXaxis()->GetNbins(),
980  hData1DReduced->GetXaxis()->GetNbins(),0,
981  hData1DReduced->GetXaxis()->GetNbins());
982 
983 
984  hCovHolder->SetTitle(binTitle.c_str());
985  hCorHolder->SetTitle(binTitle.c_str());
986  hCovHolderPPFX->SetTitle(binTitle.c_str());
987  hCorHolderPPFX->SetTitle(binTitle.c_str());
988 
989  //Get Ready to calculate the covariance and correlation matrices
990  std::vector<TH1D*> hMultiverseHistsGENIE;
991  std::vector<TH1D*> hMultiverseHistsPPFX;
992 
993  //Performing area normalization
994  //GENIE Universes
995  for(uint i = 0; i < hSystsMultiverse.size()- 100; i++){
996  TH1D* hHolder=
997  (TH1D*)hSystsMultiverseReduced[i]->Clone(ana::UniqueName().c_str());
998  float ymax = hHolder->Integral();
999  float cvmax = hCV1DReduced->Integral();
1000  hHolder->Scale(cvmax/ymax);
1001  hMultiverseHistsGENIE.push_back(hHolder);
1002  }
1003 
1004  //PPFX Universes
1005  for(uint i = hSystsMultiverse.size()-100;
1006  i < hSystsMultiverse.size(); i++){
1007  TH1D* hHolder=
1008  (TH1D*)hSystsMultiverseReduced[i]->Clone(ana::UniqueName().c_str());
1009  float ymax = hHolder->Integral();
1010  float cvmax = hCV1DReduced->Integral();
1011  hHolder->Scale(cvmax/ymax);
1012  hMultiverseHistsPPFX.push_back(hHolder);
1013  }
1014 
1015  //Covariance and Correlation Calculations
1016  //Use an external Class from trunk/CAFAna/XSec/MultiverseCorrelation.h
1017  ana::MultiverseCorrelation* multi_correlation_calc_genie =
1018  new ana::MultiverseCorrelation(hMultiverseHistsGENIE);
1019 
1020  TH2F* hCov =
1021  (TH2F*)multi_correlation_calc_genie->GetCovarianceMatrix();
1022  TH2F* hCor =
1023  (TH2F*)multi_correlation_calc_genie->GetCorrelationMatrix();
1024 
1025  ana::MultiverseCorrelation* multi_correlation_calc_ppfx =
1026  new ana::MultiverseCorrelation(hMultiverseHistsPPFX);
1027 
1028  TH2F* hCovPPFX =
1029  (TH2F*)multi_correlation_calc_ppfx->GetCovarianceMatrix();
1030  TH2F* hCorPPFX =
1031  (TH2F*)multi_correlation_calc_ppfx->GetCorrelationMatrix();
1032 
1033  for(int i = 1; i <= hCov->GetXaxis()->GetNbins(); i++){
1034  for(int j = 1; j <= hCov->GetXaxis()->GetNbins(); j++){
1035  hCovHolder->SetBinContent(i,j,hCov->GetBinContent(i,j));
1036  hCorHolder->SetBinContent(i,j,hCor->GetBinContent(i,j));
1037  hCovHolderPPFX->SetBinContent(i,j,hCovPPFX->GetBinContent(i,j));
1038  hCorHolderPPFX->SetBinContent(i,j,hCorPPFX->GetBinContent(i,j));
1039  }
1040  }
1041 
1042  //Make Plots
1043  if(makePlots){
1044  hCV1DReduced->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
1045  hCV1DReduced->GetYaxis()->SetTitle("Template Bins");
1046  hCV1DReduced->GetYaxis()->CenterTitle();
1047  hCV1DReduced->GetXaxis()->CenterTitle();
1048  hCV1DReduced->SetTitle(binTitle.c_str());
1049 
1050  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
1051  hCV1DReduced->Draw();
1052  for(uint i = 0; i < hMultiverseHistsPPFX.size(); i++){
1053  hMultiverseHistsPPFX[i]->SetLineColorAlpha(kBlue-9,0.50);
1054  hMultiverseHistsPPFX[i]->Draw("hist same");
1055  }
1056  hCV1DReduced->Draw("same");
1057  c1->SaveAs(Form("%s%s_%s_%s_only.png","Covariance/",
1058  "standard_template_pseudo_exp",
1059  "ppfx","reduced_template"));
1060  c1->Close();
1061  }
1062 
1063  if(makePlots){
1064  hCV1DReduced->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
1065  hCV1DReduced->GetYaxis()->SetTitle("Template Bins");
1066  hCV1DReduced->GetYaxis()->CenterTitle();
1067  hCV1DReduced->GetXaxis()->CenterTitle();
1068  hCV1DReduced->SetTitle(binTitle.c_str());
1069 
1070  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
1071  hCV1DReduced->Draw();
1072  for(uint i = 0; i < hMultiverseHistsGENIE.size(); i++){
1073  hMultiverseHistsGENIE[i]->SetLineColorAlpha(kBlue-9,0.50);
1074  hMultiverseHistsGENIE[i]->Draw("hist same");
1075  }
1076  hCV1DReduced->Draw("same");
1077  c1->SaveAs(Form("%s%s_%s_%s_only.png","Covariance/",
1078  "standard_template_pseudo_exp",
1079  "genie","reduced_template"));
1080  c1->Close();
1081  }
1082 
1083  mapCovariances.insert(std::make_pair("genie",
1084  hCovHolder));
1085  mapCorrelations.insert(std::make_pair("genie",
1086  hCorHolder));
1087  mapCovariances.insert(std::make_pair("ppfx",
1088  hCovHolderPPFX));
1089  mapCorrelations.insert(std::make_pair("ppfx",
1090  hCorHolderPPFX));
1091  delete hCorPPFX;
1092  delete hCovPPFX;
1093  delete multi_correlation_calc_ppfx;
1094  for(auto ele :hMultiverseHistsPPFX)
1095  delete ele;
1096  delete hCor;
1097  delete hCov;
1098  for(auto ele :hMultiverseHistsGENIE)
1099  delete ele;
1100  delete multi_correlation_calc_genie;
1101 
1102  }//multiverse sample calculation
1103 
1104  //Now Plot All Covariance Matrices Seperately
1105  if(makePlots){
1106  for(auto pair: mapCovariances){
1107  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
1108  pair.second->GetXaxis()->CenterTitle();
1109  pair.second->GetYaxis()->CenterTitle();
1110  pair.second->GetZaxis()->CenterTitle();
1111  pair.second->GetZaxis()->SetTitle("Covariance (Events^{2})");
1112  pair.second->SetTitle(binTitle.c_str());
1113  c3->SetLeftMargin(0.1);
1114  c3->SetRightMargin(0.15);
1115  c3->SetBottomMargin(0.1);
1116  pair.second->GetZaxis()->SetMaxDigits(3);
1117  pair.second->Draw("COLZ");
1118  c3->SaveAs(Form("%s%s_%s_%s_only.png","Covariance/",
1119  "covariance_single",
1120  pair.first.c_str(),"reduced_template"));
1121  c3->Close();
1122  delete c3;
1123  }
1124  }//if makePlots
1125 
1126  if(makePlots){
1127  for(auto pair: mapCorrelations){
1128 
1129  for(uint i = 1; i < pair.second->GetXaxis()->GetNbins(); i++){
1130  if(pair.second->GetBinContent(i,i) != 1)
1131  std::cout << "We have a problem here" << std::endl;
1132  }
1133 
1134  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
1135  pair.second->GetXaxis()->CenterTitle();
1136  pair.second->GetYaxis()->CenterTitle();
1137  pair.second->GetZaxis()->CenterTitle();
1138  pair.second->SetTitle(binTitle.c_str());
1139  pair.second->GetZaxis()->SetTitle("Correlation");
1140  pair.second->GetZaxis()->SetRangeUser(-1,1);
1141  c3->SetLeftMargin(0.1);
1142  c3->SetRightMargin(0.15);
1143  c3->SetBottomMargin(0.1);
1144  //pair.second->GetZaxis()->SetMaxDigits(3);
1145  pair.second->Draw("COLZ");
1146  c3->SaveAs(Form("%s%s_%s_%s_only.png",covarianceDir.c_str(),
1147  "correlation_single",
1148  pair.first.c_str(),"reduced_template"));
1149  c3->Close();
1150  delete c3;
1151  }
1152  }//if makePlots
1153 
1154  //Now start adding up each of the individual covariance matrices
1155  TH2F* hCovariance =
1156  new TH2F(ana::UniqueName().c_str(),
1157  ";Template Bins;Template Bins; Covariance (Events^{2})",
1158  hData1DReduced->GetXaxis()->GetNbins(),0,
1159  hData1DReduced->GetXaxis()->GetNbins(),
1160  hData1DReduced->GetXaxis()->GetNbins(),0,
1161  hData1DReduced->GetXaxis()->GetNbins());
1162 
1163 
1164 
1165 
1166  std::string output = covarianceDir + "covariance_addition";
1167  for(auto pair: mapCovariances){
1168  //if(pair.first == "genie")continue;
1169  output+="_";
1170  output+=pair.first;
1171  hCovariance->Add(pair.second);
1172  if(makePlots){
1173  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
1174  c3->SetLeftMargin(0.1);
1175  c3->SetRightMargin(0.15);
1176  c3->SetBottomMargin(0.1);
1177  hCovariance->SetTitle(binTitle.c_str());
1178  hCovariance->GetZaxis()->SetMaxDigits(3);
1179  hCovariance->GetZaxis()->SetTitle("Covariance (Events^{2})");
1180  hCovariance->Draw("COLZ");
1181  c3->SaveAs(Form("%s_%s.png",output.c_str(),"reduced_template"));
1182  c3->Close();
1183  delete c3;
1184  }
1185  }
1186 
1187  TH2F* hCorrelation =
1188  (TH2F*)hCovariance->Clone(ana::UniqueName().c_str());
1189 
1190  std::cout << hCorrelation->GetXaxis()->GetNbins() << "\t"
1191  << hCorrelation->GetYaxis()->GetNbins() << std::endl;
1192 
1193  for(int i=0; i<hCovariance->GetXaxis()->GetNbins(); i++){
1194  for(int j=0; j<hCovariance->GetXaxis()->GetNbins(); j++){
1195  double covariance = hCovariance->GetBinContent(i+1,j+1);
1196  double variance_i = sqrt(hCovariance->GetBinContent(i+1,i+1));
1197  double variance_j = sqrt(hCovariance->GetBinContent(j+1,j+1));
1198  double correlation = covariance/(variance_i*variance_j);
1199  if(isnan(correlation)) correlation = 0;
1200  if(i == j && correlation != 1)
1201  std::cout << "We have a problem here " << i << "\t" <<
1202  correlation << std::endl;
1203  if(i == j) correlation = 1;
1204  hCorrelation->SetBinContent(i+1,j+1,
1205  correlation);
1206  }
1207  }
1208 
1209  if(makePlots){
1210  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
1211  hCorrelation->GetXaxis()->CenterTitle();
1212  hCorrelation->GetYaxis()->CenterTitle();
1213  hCorrelation->GetZaxis()->CenterTitle();
1214  hCorrelation->SetTitle(binTitle.c_str());
1215  hCorrelation->GetZaxis()->SetTitle("Correlation");
1216  hCorrelation->GetZaxis()->SetRangeUser(-1,1);
1217  c3->SetLeftMargin(0.1);
1218  c3->SetRightMargin(0.15);
1219  c3->SetBottomMargin(0.1);
1220  //hCorrelation->GetYaxis()->SetMaxDigits(3);
1221  hCorrelation->Draw("COLZ");
1222  c3->SaveAs(Form("%s%s_%s.png",covarianceDir.c_str(),
1223  "correlation_all","reduced_template"));
1224  c3->Close();
1225  delete c3;
1226  }
1227 
1228  //Move the covariance matrix to the holder so that it can be used
1229  for(int i = 0; i < hCovariance->GetXaxis()->GetNbins(); i++){
1230  std::vector<float> hCovCol;
1231  for(int j = 0; j < hCovariance->GetXaxis()->GetNbins(); j++){
1232  hCovCol.push_back(hCovariance->GetBinContent(i+1,j+1));
1233  }
1234  vCovarianceMatrixReduced.push_back(hCovCol);
1235  }
1236 
1237  //Clean up memory a bit
1238  for(uint i = 0; i < hSystsShifts.size(); i++){
1239  delete hSystsShifts[i];
1240  }
1241  hSystsShifts.clear();
1242  }
1243 
1244  //Fill the diagonal statistical covariance matrix from the full template
1245  //space
1247  {
1248  for(int i = 0; i < hData1DReduced->GetXaxis()->GetNbins(); i++){
1249  std::vector<float> StatCovCol;
1250  for(int j = 0; j < hData1DReduced->GetXaxis()->GetNbins(); j++){
1251  if(i == j){
1252  StatCovCol.push_back(hCV1DReduced->GetBinError(i+1)*
1253  hCV1DReduced->GetBinError(i+1) +
1254  hData1DReduced->GetBinError(i+1)*
1255  hData1DReduced->GetBinError(i+1));
1256  }
1257  else
1258  StatCovCol.push_back(0);
1259  }
1260  vStatCovarianceMatrixReduced.push_back(StatCovCol);
1261  }
1262  }
1263 
1264  //Get Values from the Systematic or Statistical Covariance Matrix
1266  {
1267  return vCovarianceMatrixReduced[binx][biny];
1268  }
1269 
1271  {
1272  return vStatCovarianceMatrixReduced[binx][biny];
1273  }
1274 
1275  //Invert Syst + Stat Covariance Matrix
1277  {
1278  std::cout << "Inverting" << std::endl;
1279  TMatrixD statMatrixReduced= TMatrixD(vDataValuesReduced.size(),
1280  vDataValuesReduced.size());
1281  TMatrixD covMatrixReduced= TMatrixD(vDataValuesReduced.size(),
1282  vDataValuesReduced.size());
1283 
1284  for(int i = 0; i < statMatrixReduced.GetNcols(); i++){
1285  for(int j = 0; j < statMatrixReduced.GetNcols(); j++){
1286  statMatrixReduced[i][j] = vStatCovarianceMatrixReduced[i][j];
1287  covMatrixReduced[i][j] = vCovarianceMatrixReduced[i][j];
1288  }
1289  }
1290 
1291  TMatrixD H = statMatrixReduced + covMatrixReduced;
1292  TDecompSVD svd(H);
1293  TMatrixD VijInv = svd.Invert();
1294 
1295 
1296  for(int i = 0; i < VijInv.GetNcols(); i++){
1297  std::vector<float> vInvCovCol;
1298  for(int j = 0; j < VijInv.GetNcols(); j++){
1299  vInvCovCol.push_back(VijInv[i][j]);
1300  }
1301  vInvertedCovarianceMatrixReduced.push_back(vInvCovCol);
1302  }
1303  std::cout << "Done" << std::endl;
1304  }
1305 
1306  //Initiallizing Fit Parameters
1308  {
1309 
1310  int numBin = 0;
1311  int parNumber = 0;
1312  for(int xcol = 0; xcol < numPIDBins; xcol++){
1313  const int nmodnynz = xcol%TotalYZBins;
1314  const int ix = xcol/TotalYZBins;
1315  const int iy = nmodnynz/TotalZBins;
1316 
1317  if((vCovarianceMatrixFull[xcol][xcol] != 0 &&
1318  vStatCovarianceMatrixFull[xcol][xcol] != 0) &&
1319  ( ((ix+1 == 5) && (iy+1 < 5) && (iy+1 > 2)) ||
1320  ((ix+1 == 6) && (iy+1 < 7) && (iy+1 > 2)) ||
1321  ((ix+1 == 7) && (iy+1 < 9) && (iy+1 > 2)) ||
1322  ((ix+1 == 8) && (iy+1 < 12) && (iy+1 > 2)) ))
1323  {
1324  numBin++;
1325  if(parNumMin->GetBinContent(ix+1,iy+1) == -1){
1326  /*if( (ix+1 ==5) && (iy+1 == 3)){
1327  parNumMin->SetBinContent(ix+1,iy+1,parNumber);
1328  parNumMax->SetBinContent(ix+1,iy+1,parNumber+1);
1329  reducedBinsMin->SetBinContent(ix+1,iy+1,numBin);
1330  if(printVerbose)
1331  std::cout << xcol << "\t" << numBin <<
1332  "\t" << ix+1 << "\t" << iy+1 << "\t" << parNumber <<
1333  "\t" << parNumber +1 << std::endl;
1334  parNumber += 2;
1335  }*/
1336  if( (ix+1 ==5) && (iy+1 == 4)){
1337  parNumMin->SetBinContent(ix+1,iy+1,parNumber);
1338  parNumMax->SetBinContent(ix+1,iy+1,parNumber+1);
1339  reducedBinsMin->SetBinContent(ix+1,iy+1,numBin);
1340  if(printVerbose)
1341  std::cout << xcol << "\t" << numBin <<
1342  "\t" << ix+1 << "\t" << iy+1 << "\t" << parNumber <<
1343  "\t" << parNumber +1 << std::endl;
1344  parNumber += 2;
1345  }
1346  else if( (ix+1 ==6) && (iy+1 == 6)){
1347  parNumMin->SetBinContent(ix+1,iy+1,parNumber);
1348  parNumMax->SetBinContent(ix+1,iy+1,parNumber+1);
1349  reducedBinsMin->SetBinContent(ix+1,iy+1,numBin);
1350  if(printVerbose)
1351  std::cout << xcol << "\t" << numBin <<
1352  "\t" << ix+1 << "\t" << iy+1 << "\t" << parNumber <<
1353  "\t" << parNumber +1 << std::endl;
1354  parNumber += 2;
1355  }
1356  else if( (ix+1 ==7) && (iy+1 == 7)){
1357  parNumMin->SetBinContent(ix+1,iy+1,parNumber);
1358  parNumMax->SetBinContent(ix+1,iy+1,parNumber+1);
1359  reducedBinsMin->SetBinContent(ix+1,iy+1,numBin);
1360  if(printVerbose)
1361  std::cout << xcol << "\t" << numBin <<
1362  "\t" << ix+1 << "\t" << iy+1 << "\t" << parNumber <<
1363  "\t" << parNumber +1 << std::endl;
1364  parNumber += 2;
1365  }
1366  else if( (ix+1 ==7) && (iy+1 == 8)){
1367  parNumMin->SetBinContent(ix+1,iy+1,parNumber);
1368  parNumMax->SetBinContent(ix+1,iy+1,parNumber+1);
1369  reducedBinsMin->SetBinContent(ix+1,iy+1,numBin);
1370  if(printVerbose)
1371  std::cout << xcol << "\t" << numBin <<
1372  "\t" << ix+1 << "\t" << iy+1 << "\t" << parNumber <<
1373  "\t" << parNumber +1 << std::endl;
1374  parNumber += 2;
1375  }
1376  else if( (ix+1 ==8) && (iy+1 == 10)){
1377  parNumMin->SetBinContent(ix+1,iy+1,parNumber);
1378  parNumMax->SetBinContent(ix+1,iy+1,parNumber+1);
1379  reducedBinsMin->SetBinContent(ix+1,iy+1,numBin);
1380  if(printVerbose)
1381  std::cout << xcol << "\t" << numBin <<
1382  "\t" << ix+1 << "\t" << iy+1 << "\t" << parNumber <<
1383  "\t" << parNumber +1 << std::endl;
1384  parNumber += 2;
1385  }
1386  else if( (ix+1 ==8) && (iy+1 == 11)){
1387  parNumMin->SetBinContent(ix+1,iy+1,parNumber);
1388  parNumMax->SetBinContent(ix+1,iy+1,parNumber+1);
1389  reducedBinsMin->SetBinContent(ix+1,iy+1,numBin);
1390  if(printVerbose)
1391  std::cout << xcol << "\t" << numBin <<
1392  "\t" << ix+1 << "\t" << iy+1 << "\t" << parNumber <<
1393  "\t" << parNumber +1 << std::endl;
1394  parNumber += 2;
1395  }
1396  else{
1397  parNumMin->SetBinContent(ix+1,iy+1,parNumber);
1398  parNumMax->SetBinContent(ix+1,iy+1,parNumber+2);
1399  reducedBinsMin->SetBinContent(ix+1,iy+1,numBin);
1400 
1401  if(printVerbose)
1402  std::cout << xcol << "\t" << numBin <<
1403  "\t" << ix+1 << "\t" << iy+1 << "\t" << parNumber <<
1404  "\t" << parNumber +2 << std::endl;
1405  parNumber += 3;
1406  }
1407  }//change parNumMin the first time
1408  }//if in phase space
1409  }//loop over full bins
1410  }
1411 
1412  std::pair<float,float> NueCCIncGlobalFitter::searchFitInfo
1413  (float value,std::string fitSearch)
1414  {
1415 
1416  int binX = -1;
1417  int binY = -1;
1418  for(int i = 1; i <= TotalXBins; i++){
1419  for(int j = 1; j <= TotalYBins; j++){
1420  if(fitSearch == "parNum")
1421  {
1422  if(value <= parNumMax->GetBinContent(i,j) &&
1423  value >= parNumMin->GetBinContent(i,j))
1424  {
1425  return std::make_pair(parNumMin->GetBinContent(i,j),
1426  parNumMax->GetBinContent(i,j));
1427  }
1428  }
1429  if(fitSearch == "binNum")
1430  {
1431  if((value >= reducedBinsMin->GetBinContent(i,j)) &&
1432  (reducedBinsMin->GetBinContent(i,j) > 0))
1433  {
1434  binX = i;
1435  binY = j;
1436  }
1437  }
1438  if(fitSearch == "analysisBin")
1439  {
1440  if(value <= parNumMax->GetBinContent(i,j) &&
1441  value >= parNumMin->GetBinContent(i,j))
1442  {
1443  return std::make_pair(i,j);
1444  }
1445  }
1446  }//loop over y bins
1447  }//loop over x bins
1448 
1449  if((fitSearch == "binNum") && (binX > 0) && (binY > 0))
1450  {
1451  return std::make_pair(parNumMin->GetBinContent(binX,binY),
1452  parNumMax->GetBinContent(binX,binY));
1453  }
1454 
1455 
1456  return std::make_pair(-1.f,-1.f);
1457  }
1458 
1459  //Fitting Function
1460  void NueCCIncGlobalFitter::fcn(Int_t &npar, Double_t *gin,
1461  Double_t &f, Double_t *par,
1462  Int_t iflag)
1463  {
1464  std::vector<double> parameters;
1465  for(int i = 0; i < numberOfFitParametersFCN; i++){
1466  parameters.push_back(par[i]);
1467  }
1468  f = Fitter_obj->myFunction(parameters);
1469  }
1470 
1471 
1472  double NueCCIncGlobalFitter::myFunction(std::vector<double>par)
1473  {
1474  float data,signal,numu,nc,other = 0;
1475 
1478 
1479  for(int i = 0; i < VijInv.GetNcols(); i++){
1480  for(int j = 0; j < VijInv.GetNcols(); j++){
1481  VijInv[i][j] = vInvertedCovarianceMatrixReduced[i][j];
1482  }
1483  }
1484 
1485  TMatrixD FitVectorRow = TMatrixD(1,vInvertedCovarianceMatrixReduced.size());
1486  TMatrixD FitVectorCol = TMatrixD(vInvertedCovarianceMatrixReduced.size(),1);
1487  //Calculate Chi Squared
1488 
1489  for(int i = 0; i < (int)vDataValuesReduced.size(); i++){
1490 
1491  std::pair<float,float> minMaxParNum =
1493  "binNum");
1494 
1495  data = Fitter_obj->vDataValuesReduced[i];
1496  signal = Fitter_obj->vSignalLike_ValuesReduced[i];
1497  numu = Fitter_obj->vNumuLike_ValuesReduced[i];
1498  nc = Fitter_obj->vNC_ValuesReduced[i];
1499  other = Fitter_obj->vOther_ValuesReduced[i];
1500 
1501  if( (data < 30) || (signal + numu + nc + other < 30)){
1502  FitVectorRow[0][i] = 0;
1503  FitVectorCol[i][0] = 0;
1504  continue;
1505  }
1506 
1507  //Background + Signal Normalization Parameters Only
1508  if((minMaxParNum.second - minMaxParNum.first) == 1){
1509  FitVectorRow[0][i] = ( data - (par[minMaxParNum.second]*signal +
1510  par[minMaxParNum.first]*numu +
1511  par[minMaxParNum.first]*nc + other));
1512 
1513  FitVectorCol[i][0] = ( data -(par[minMaxParNum.second]*signal +
1514  par[minMaxParNum.first]*numu +
1515  par[minMaxParNum.first]*nc + other));
1516  }
1517  //NumuCC, Signal, NC Normalization Parameters
1518  else if((minMaxParNum.second - minMaxParNum.first) == 2){
1519  FitVectorRow[0][i] = ( data - (par[minMaxParNum.first+1]*signal +
1520  par[minMaxParNum.first]*numu +
1521  par[minMaxParNum.second]*nc + other));
1522 
1523  FitVectorCol[i][0] = ( data -(par[minMaxParNum.first+1]*signal +
1524  par[minMaxParNum.first]*numu +
1525  par[minMaxParNum.second]*nc + other));
1526  }
1527  else{
1528  FitVectorRow[0][i] = 0;
1529  FitVectorCol[i][0] = 0;
1530  }
1531 
1532  }
1533 
1534  TMatrixD Chi = FitVectorRow * VijInv;
1535  TMatrix Chi2 = Chi * FitVectorCol;
1536  return Chi2(0,0);
1537  }
1538 
1539  void NueCCIncGlobalFitter::FillParameterCovariance(std::vector<std::vector<double>> cov)
1540  {
1541  for(int i = 0; i < numberOfFitParametersFCN; i++){
1542  std::vector<double> covHolder;
1543  std::vector<double> corHolder;
1544  for(int j = 0; j < numberOfFitParametersFCN; j++){
1545  covHolder.push_back(cov[i][j]);
1546  corHolder.push_back(cov[i][j]/(pow(cov[i][i],0.5)*
1547  pow(cov[j][j],0.5)));
1548  }
1549  parCovariance.push_back(covHolder);
1550  parCorrelation.push_back(corHolder);
1551  }
1552  }
1553 
1555  {
1556  for(int binx = 1; binx <= hAnalysis3D[0]->GetXaxis()->GetNbins(); binx++){
1557  for(int biny = 1; biny <= hAnalysis3D[0]->GetYaxis()->GetNbins(); biny++){
1558 
1559  std::pair<double,double>
1560  nue_wei(nue_normalization->GetBinContent(binx,biny),
1561  nue_uncertainty->GetBinContent(binx,biny));
1562  std::pair<double,double>
1563  numu_wei(numu_normalization->GetBinContent(binx,biny),
1564  numu_uncertainty->GetBinContent(binx,biny));
1565  std::pair<double,double>
1566  nc_wei(nc_normalization->GetBinContent(binx,biny),
1567  nc_uncertainty->GetBinContent(binx,biny));
1568 
1569  //Hard coded a bit
1570  //gives weights from the fit, to adjacent areas just outside of the
1571  //fit phase space, used to stabilize unfolding results
1572  //if(binx == 5 && biny > 4 && nue_wei.second == 0 )
1573  //nue_wei.first = nue_normalization->GetBinContent(5,4);
1574  //else if(binx <= 5 && nue_wei.second == 0 )
1575  //nue_wei.first = nue_normalization->GetBinContent(5,3);
1576  /*
1577  if(binx == 8 && biny > 11 && nue_wei.second == 0 )
1578  nue_wei.first = nue_normalization->GetBinContent(8,11);
1579  if(binx == 7 && biny > 8)
1580  nue_wei.first = nue_normalization->GetBinContent(7,8);
1581  if(binx == 6 && biny > 6)
1582  nue_wei.first = nue_normalization->GetBinContent(6,6);
1583 
1584  if(binx == 8 && biny < 4 && nue_wei.second == 0 )
1585  nue_wei.first = nue_normalization->GetBinContent(8,4);
1586  if(binx == 7 && biny < 3)
1587  nue_wei.first = nue_normalization->GetBinContent(7,3);
1588  if(binx == 6 && biny < 3)
1589  nue_wei.first = nue_normalization->GetBinContent(6,3);
1590  */
1591 
1592 
1593  for(int zbin = 1; zbin <= hAnalysis3D[0]->GetZaxis()->GetNbins();zbin++)
1594  {
1595  float nsig = hAnalysis3D[1]->GetBinContent(binx,biny,zbin);
1596  float nnuebar = hAnalysis3D[2]->GetBinContent(binx,biny,zbin);
1597  float nnumu = hAnalysis3D[3]->GetBinContent(binx,biny,zbin);
1598  float nnumubar = hAnalysis3D[4]->GetBinContent(binx,biny,zbin);
1599  float nnc = hAnalysis3D[5]->GetBinContent(binx,biny,zbin);
1600  float nother = hAnalysis3D[6]->GetBinContent(binx,biny,zbin);
1601  float nnf = hAnalysis3D[7]->GetBinContent(binx,biny,zbin);
1602 
1603 
1604  float err_nue =
1605  (nue_wei.first * nsig) *
1606  sqrt( pow(nue_wei.second/nue_wei.first,2) +
1607  pow( sqrt(nsig + nnuebar + nnf)/(nsig+nnuebar+nnf),2));
1608 
1609  float err_nuebar =
1610  (nue_wei.first * nnuebar) *
1611  sqrt( pow(nue_wei.second/nue_wei.first,2) +
1612  pow( sqrt(nsig + nnuebar + nnf)/(nsig+nnuebar+nnf),2));
1613 
1614  float err_nnf =
1615  (nue_wei.first * nnf) *
1616  sqrt( pow(nue_wei.second/nue_wei.first,2) +
1617  pow( sqrt(nsig + nnuebar + nnf)/(nsig+nnuebar+nnf),2));
1618 
1619  float err_numu =
1620  (numu_wei.first * nnumu) *
1621  sqrt( pow(numu_wei.second/numu_wei.first,2) +
1622  pow( sqrt(nnumu + nnumubar)/(nnumu + nnumubar),2));
1623 
1624  float err_numubar =
1625  (numu_wei.first * nnumubar) *
1626  sqrt( pow(numu_wei.second/numu_wei.first,2) +
1627  pow( sqrt(nnumu + nnumubar)/(nnumu + nnumubar),2));
1628 
1629  float err_nc =
1630  (nc_wei.first * nnc) *
1631  sqrt( pow(nc_wei.second/nc_wei.first,2) +
1632  pow( sqrt(nnc)/(nnc),2));
1633 
1634 
1635  float err_total =
1636  sqrt( pow(err_nue,2) + pow(err_numu,2) + pow(err_nc,2)
1637  + nother);
1638 
1639  if(isnan(err_total) ||
1640  nue_wei.second == 0 ||
1641  numu_wei.second == 0 || nc_wei.second == 0) continue;
1642  float tot_nue = nue_wei.first * nsig;
1643  float tot_nuebar = nue_wei.first * nnuebar;
1644  float tot_numu = numu_wei.first * nnumu;
1645  float tot_numubar = numu_wei.first * nnumubar;
1646  float tot_nc = nc_wei.first * nnc;
1647  float tot_nnf = nue_wei.first * nnf;
1648 
1649  hAnalysis3D[0]->SetBinContent(binx,biny,zbin,tot_nue+
1650  tot_nuebar +
1651  tot_numu +
1652  tot_numubar +
1653  tot_nc + tot_nnf + nother);
1654  hAnalysis3D[1]->SetBinContent(binx,biny,zbin,tot_nue);
1655  hAnalysis3D[2]->SetBinContent(binx,biny,zbin,tot_nuebar);
1656  hAnalysis3D[3]->SetBinContent(binx,biny,zbin,tot_numu);
1657  hAnalysis3D[4]->SetBinContent(binx,biny,zbin,tot_numubar);
1658  hAnalysis3D[5]->SetBinContent(binx,biny,zbin,tot_nc);
1659  hAnalysis3D[7]->SetBinContent(binx,biny,zbin,tot_nnf);
1660 
1661 
1662  hAnalysis3D[0]->SetBinError(binx,biny,zbin,err_total);
1663  hAnalysis3D[1]->SetBinError(binx,biny,zbin,err_nue);
1664  hAnalysis3D[2]->SetBinError(binx,biny,zbin,err_nuebar);
1665  hAnalysis3D[3]->SetBinError(binx,biny,zbin,err_numu);
1666  hAnalysis3D[4]->SetBinError(binx,biny,zbin,err_numubar);
1667  hAnalysis3D[5]->SetBinError(binx,biny,zbin,err_nc);
1668  hAnalysis3D[7]->SetBinError(binx,biny,zbin,err_nnf);
1669 
1670 
1671  std::cout << binx << "\t" << biny << "\t" <<
1672  zbin << "\t" << err_total << "\t" <<
1673  err_nue/tot_nue << "\t" << err_numu/tot_numu <<
1674  "\t" << err_nc/tot_nc << std::endl;
1675 
1676 
1677  }//loop over zbins
1678 
1679 
1680  }//loop over ybins
1681  }//loop over xbins
1682  }
1683 
1685  {
1686  std::vector<TH1F*> hResults;
1687  hResults.push_back((TH1F*)hData1D->Clone(ana::UniqueName().c_str()));
1688  for(uint i = 0; i < hTemplates1D.size(); i++){
1689  hResults.push_back((TH1F*)
1690  hTemplates1D[i]->Clone(ana::UniqueName().c_str()));
1691  }
1692  return hResults;
1693  }
1694 
1696  {
1697 
1698  TH1F* hDataTemplate =
1699  new TH1F(ana::UniqueName().c_str(),
1700  ";Template Bins; Events/8.09#times10^{20} POT",
1701  vDataValuesReduced.size(),
1702  0,vDataValuesReduced.size());
1703  TH1F* hSignalLikeTemplate =
1704  new TH1F(ana::UniqueName().c_str(),
1705  ";Template Bins; Events/8.09#times10^{20} POT",
1706  vDataValuesReduced.size(),
1707  0,vDataValuesReduced.size());
1708  TH1F* hNumuLikeTemplate =
1709  new TH1F(ana::UniqueName().c_str(),
1710  ";Template Bins; Events/8.09#times10^{20} POT",
1711  vDataValuesReduced.size(),
1712  0,vDataValuesReduced.size());
1713  TH1F* hNCTemplate =
1714  new TH1F(ana::UniqueName().c_str(),
1715  ";Template Bins; Events/8.09#times10^{20} POT",
1716  vDataValuesReduced.size(),
1717  0,vDataValuesReduced.size());
1718  TH1F* hOtherTemplate =
1719  new TH1F(ana::UniqueName().c_str(),
1720  ";Template Bins; Events/8.09#times10^{20} POT",
1721  vDataValuesReduced.size(),
1722  0,vDataValuesReduced.size());
1723 
1724  for(int i = 0; i < (int)vDataValuesReduced.size(); i++){
1725  hDataTemplate->SetBinContent(i+1,vDataValuesReduced[i]);
1726  hSignalLikeTemplate->SetBinContent(i+1,vSignalLike_ValuesReduced[i]);
1727  hNumuLikeTemplate->SetBinContent(i+1,vNumuLike_ValuesReduced[i]);
1728  hNCTemplate->SetBinContent(i+1,vNC_ValuesReduced[i]);
1729  hOtherTemplate->SetBinContent(i+1,vOther_ValuesReduced[i]);
1730  }
1731 
1732  return {hDataTemplate,hSignalLikeTemplate,hNumuLikeTemplate,hNCTemplate,
1733  hOtherTemplate};
1734  }
1735 
1737  {
1738 
1739  TH1F* hDataTemplate =
1740  new TH1F(ana::UniqueName().c_str(),
1741  ";Template Bins; Events/8.09#times10^{20} POT",
1742  vDataValuesReduced.size(),
1743  0,vDataValuesReduced.size());
1744  TH1F* hSignalLikeTemplate =
1745  new TH1F(ana::UniqueName().c_str(),
1746  ";Template Bins; Events/8.09#times10^{20} POT",
1747  vDataValuesReduced.size(),
1748  0,vDataValuesReduced.size());
1749  TH1F* hNumuLikeTemplate =
1750  new TH1F(ana::UniqueName().c_str(),
1751  ";Template Bins; Events/8.09#times10^{20} POT",
1752  vDataValuesReduced.size(),
1753  0,vDataValuesReduced.size());
1754  TH1F* hNCTemplate =
1755  new TH1F(ana::UniqueName().c_str(),
1756  ";Template Bins; Events/8.09#times10^{20} POT",
1757  vDataValuesReduced.size(),
1758  0,vDataValuesReduced.size());
1759  TH1F* hOtherTemplate =
1760  new TH1F(ana::UniqueName().c_str(),
1761  ";Template Bins; Events/8.09#times10^{20} POT",
1762  vDataValuesReduced.size(),
1763  0,vDataValuesReduced.size());
1764 
1765  for(int i = 0; i < (int)vDataValuesReduced.size(); i++){
1766  std::pair<float,float> reducedBinParNums =
1768  "binNum");
1769  std::pair<float,float> analysisBinNums =
1770  NueCCIncGlobalFitter::searchFitInfo(reducedBinParNums.first,
1771  "analysisBin");
1772  float signal_weight =
1773  nue_normalization->GetBinContent(analysisBinNums.first,
1774  analysisBinNums.second);
1775  float numu_weight =
1776  numu_normalization->GetBinContent(analysisBinNums.first,
1777  analysisBinNums.second);
1778  float nc_weight =
1779  nc_normalization->GetBinContent(analysisBinNums.first,
1780  analysisBinNums.second);
1781  hDataTemplate->SetBinContent(i+1,vDataValuesReduced[i]);
1782  hSignalLikeTemplate->SetBinContent(i+1,vSignalLike_ValuesReduced[i]*
1783  signal_weight);
1784  hNumuLikeTemplate->SetBinContent(i+1,vNumuLike_ValuesReduced[i]*
1785  numu_weight);
1786  hNCTemplate->SetBinContent(i+1,vNC_ValuesReduced[i]*
1787  nc_weight);
1788  hOtherTemplate->SetBinContent(i+1,vOther_ValuesReduced[i]);
1789  }
1790 
1791  return {hDataTemplate,hSignalLikeTemplate,hNumuLikeTemplate,hNCTemplate,
1792  hOtherTemplate};
1793  }
1794 
1796  {
1797  std::vector<TH3F*> hResults;
1798  for(uint i = 0; i < hAnalysis3D.size(); i++){
1799  hResults.push_back((TH3F*)
1800  hAnalysis3D[i]->Clone(ana::UniqueName().c_str()));
1801  }
1802  return hResults;
1803  }
1804 
1805  //////////// //////////// //////////// //////////// ////////////
1806  //Public Functions//////////////////////////////////////////////////
1807  //////////// //////////// //////////// //////////// ////////////
1809  {
1810  Fitter_obj = this;
1811  }
1812 
1813  //The main function of this class
1814  //Fills Everything needed for the fit
1815  //Performs the fit
1816  //Returns the results of the fit
1817  std::vector<TH3F*> NueCCIncGlobalFitter::doFullFit()
1818  {
1819  std::vector<TH3F*> badentry;
1820  bool CheckSysts = NueCCIncGlobalFitter::CheckSystsSize();
1821  if(!CheckSysts){
1822  std::cerr << "Unequal up and down shifts. Exiting" << std::endl;
1823  return badentry;
1824  }
1825 
1835 
1836  //Prepare the starting seeds
1837  std::vector<float> nue_starting_points;
1838  std::vector<float> numu_starting_points;
1839  std::vector<float> nc_starting_points;
1840 
1841  for(float fillVal = 0.61; fillVal < 1.41; fillVal += 0.30){
1842  nue_starting_points.push_back(fillVal);
1843  numu_starting_points.push_back(fillVal);
1844  nc_starting_points.push_back(fillVal);
1845  }
1846 
1847  //Start Preparing the Fitting Parameters
1848  int numberPar = -1;
1849  for(int xbin = 1; xbin <= TotalXBins; xbin++){
1850  for(int ybin = 1;ybin <= TotalYBins; ybin++){
1851  if(parNumMax->GetBinContent(xbin,ybin) > numberPar){
1852  numberPar = parNumMax->GetBinContent(xbin,ybin);
1853  }
1854  }
1855  }
1856 
1857  if(printVerbose)
1858  std::cout << "Number of Parameters: " << numberPar+1 << std::endl;
1859 
1860  const int parametersForFit = numberPar+1; //Have to count 0 as a parameter;
1861  numberOfFitParametersFCN = parametersForFit;
1862 
1863  //Hold Fit Results
1864  std::map<std::string, float> mFitChiSquared;
1865  std::map<std::string, float> mFitNueStart;
1866  std::map<std::string, float> mFitNumuStart;
1867  std::map<std::string, float> mFitNCStart;
1868 
1869  int fitNum = -1;
1870 
1871  bool beQuick = false;
1872 
1873  for(float nuestart : nue_starting_points){
1874  if(beQuick) continue;
1875  for(float numustart : numu_starting_points){
1876  if(beQuick)continue;
1877  for(float ncstart : nc_starting_points){
1878  fitNum++;
1879  TMinuit *gMinuit = new TMinuit(parametersForFit);
1880  gMinuit->SetPrintLevel(-1);
1881  gMinuit->SetFCN(fcn);
1882 
1883  Double_t arglist[4];
1884  Int_t ierflg = 0;
1885  arglist[0] = 1;
1886  gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1887 
1888  //Define Fit Parameters
1889  for(int parNum = 0; parNum < parametersForFit; parNum++)
1890  {
1891  std::pair<float,float> minMaxParNum =
1892  NueCCIncGlobalFitter::searchFitInfo((float)parNum,
1893  "parNum");
1894 
1895  if(minMaxParNum.first < 0 || minMaxParNum.second < 0) break;
1896 
1897  if(minMaxParNum.second - minMaxParNum.first == 2){
1898  if(minMaxParNum.second - parNum == 2)
1899  gMinuit->mnparm(parNum,Form("%i Numu Scaling",parNum),
1900  numustart,0.15,0,2,ierflg);
1901  else if(minMaxParNum.second - parNum == 1)
1902  gMinuit->mnparm(parNum,Form("%i Nue Scaling",parNum),
1903  nuestart,0.15,0,0,ierflg);
1904  else
1905  gMinuit->mnparm(parNum,Form("%i NC Scaling",parNum),
1906  ncstart,0.15,0,0,ierflg);
1907  }
1908  else if(minMaxParNum.second - minMaxParNum.first == 1){
1909  if(minMaxParNum.second - parNum == 1)
1910  gMinuit->mnparm(parNum,Form("%i Bkg Scaling",parNum),
1911  numustart,0.15,0,2,ierflg);
1912  else
1913  gMinuit->mnparm(parNum,Form("%i Nue Scaling",parNum),
1914  nuestart,0.15,0,0,ierflg);
1915  }
1916 
1917  }//loop over parameters
1918 
1919  std::cout << "Number of parameters: " << gMinuit->GetNumPars()
1920  << std::endl;
1921 
1922  arglist[0] = 0; //number of iterations 0 == no limit
1923  arglist[1] = 1; //1 == standard minimization strategy
1924  gMinuit->Command("SIMPLEX");
1925 
1926  Double_t fmin = 0,fedm = 0,ferrdef = 0;
1927  Int_t npari =0, nparx = 0, istat = 0;
1928  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1929  std::cout << "SIMPLEX " << istat << std::endl;
1930 
1931  gMinuit->Command("HESSE");
1932  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1933  std::cout << "HESSE " << istat << std::endl;
1934 
1935  gMinuit->Command("MIGRAD");
1936  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1937  std::cout << "MIGRAD " << istat << std::endl;
1938 
1939  if(istat == 3){
1940  gMinuit->Command("MINOS");
1941  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1942  std::cout << "MINOS " << istat << std::endl;
1943  }
1944  else{
1945  std::cout << "Fit is messed up" << std::endl;
1946  continue;
1947  }
1948 
1949 
1950  //Grab Final Results
1951  bool troubleFitting = false;
1952  //double vstart[parametersForFit];
1953  //double verror[parametersForFit];
1954  double fParStart;
1955  double fParErr;
1956 
1957  for(int parNum = 0; parNum < parametersForFit; parNum++){
1958  gMinuit->GetParameter(parNum,fParStart,fParErr);
1959  //vstart[parNum] = fParStart;
1960  //verror[parNum] = fParErr;
1961  std::cout << parNum << "\t" << fParStart << "\t"
1962  << fParErr << std::endl;
1963  if(fParStart > 2 || fParStart < 0) troubleFitting = true;
1964  }
1965 
1966  //Before Ending Fits, Make Sure There Was No Issue With the Fit
1967  Double_t corr_mat_test[parametersForFit][parametersForFit];
1968  gMinuit->mnemat(*corr_mat_test,parametersForFit);
1969 
1970  //Turn this test error matix into a correlation matrix
1971  if(istat !=3)std::cout << "Errors in fit " << istat << std::endl;
1972  if(istat != 3) troubleFitting = true;
1973 
1974  for(int i = 0; i < parametersForFit; i++){
1975  for(int j = 0; j < parametersForFit; j++){
1976  if(i == j) continue;
1977  if( corr_mat_test[i][j]/(pow(corr_mat_test[i][i],0.5)*
1978  pow(corr_mat_test[j][j],0.5))
1979  > 0.99)
1980  troubleFitting = true;
1981  std::cout << i << "\t" << j << "\t"
1982  << corr_mat_test[i][j]/(pow(corr_mat_test[i][i],0.5)*
1983  pow(corr_mat_test[j][j],0.5))
1984  << std::endl;
1985  }
1986  }
1987 
1988  if(!troubleFitting){
1989  std::cout << "**********Booyah!*************" << std::endl;
1990  std::cout << "Chi Squared " << fmin << std::endl;
1991  std::string fitName = "Fit_"+std::to_string(fitNum);
1992  mFitChiSquared.insert(std::make_pair(fitName,fmin));
1993  mFitNueStart.insert(std::make_pair(fitName,nuestart));
1994  mFitNumuStart.insert(std::make_pair(fitName,numustart));
1995  mFitNCStart.insert(std::make_pair(fitName,ncstart));
1996  beQuick = true;
1997  if(beQuick)break;
1998  }
1999 
2000  }//loop ncstart
2001  }//loop numustart
2002  }//loop nuestart
2003 
2004 
2005  //Find the best fits (lowest ChiSquared)
2006  std::string bestFitName = "NULL";
2007  float bestChiSquared = 100000;
2008  for(auto element : mFitChiSquared){
2009  if(floor(element.second*1e3) < bestChiSquared*1e3){
2010  bestFitName = element.first;
2011  bestChiSquared = element.second;
2012  }
2013  }
2014 
2015  if(bestFitName == "NULL"){
2016  std::cout << "No successful fits" << std::endl;
2017  return badentry;
2018  }
2019 
2020  std::cout << "***********************************************" << std::endl;
2021  std::cout << "***********************************************" << std::endl;
2022  std::cout << "***********************************************" << std::endl;
2023  std::cout << "***********************************************" << std::endl;
2024  std::cout << "***********************************************" << std::endl;
2025  std::cout << "***********************************************" << std::endl;
2026  std::cout << "Moving on to final fit" << std::endl;
2027  std::cout << "***********************************************" << std::endl;
2028  std::cout << "***********************************************" << std::endl;
2029  std::cout << "***********************************************" << std::endl;
2030  std::cout << "***********************************************" << std::endl;
2031  std::cout << "***********************************************" << std::endl;
2032  std::cout << "***********************************************" << std::endl;
2033 
2034  //Final Fit
2035  //This performs the fit one last time and
2036  //saves all necessary information from the fit
2037  TMinuit *jMinuit = new TMinuit(parametersForFit);
2038  jMinuit->SetPrintLevel(-1);
2039  jMinuit->SetFCN(fcn);
2040 
2041  Double_t arglist[4];
2042  Int_t ierflg = 0;
2043  arglist[0] = 1;
2044  jMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
2045 
2046  //Define Fit Parameters
2047  for(int parNum = 0; parNum < parametersForFit; parNum++)
2048  {
2049  std::pair<float,float> minMaxParNum =
2050  NueCCIncGlobalFitter::searchFitInfo((float)parNum,
2051  "parNum");
2052 
2053  if(minMaxParNum.first < 0 || minMaxParNum.second < 0) break;
2054 
2055  if(minMaxParNum.second - minMaxParNum.first == 2){
2056  if(minMaxParNum.second - parNum == 2)
2057  jMinuit->mnparm(parNum,Form("%i Numu Scaling",parNum),
2058  mFitNumuStart[bestFitName],0.15,0,0,ierflg);
2059  else if(minMaxParNum.second - parNum == 1)
2060  jMinuit->mnparm(parNum,Form("%i Nue Scaling",parNum),
2061  mFitNueStart[bestFitName],0.15,0,0,ierflg);
2062  else
2063  jMinuit->mnparm(parNum,Form("%i NC Scaling",parNum),
2064  mFitNCStart[bestFitName],0.15,0,0,ierflg);
2065  }
2066  else if(minMaxParNum.second - minMaxParNum.first == 1){
2067  if(minMaxParNum.second - parNum == 1)
2068  jMinuit->mnparm(parNum,Form("%i Bkg Scaling",parNum),
2069  mFitNumuStart[bestFitName],0.15,0,0,ierflg);
2070  else
2071  jMinuit->mnparm(parNum,Form("%i Nue Scaling",parNum),
2072  mFitNueStart[bestFitName],0.15,0,0,ierflg);
2073  }
2074 
2075  }//loop over parameters
2076 
2077  std::cout << "Number of parameters: " << gMinuit->GetNumPars()
2078  << std::endl;
2079 
2080  arglist[0] = 0; //number of iterations 0 == no limit
2081  arglist[1] = 1; //1 == standard minimization strategy
2082  jMinuit->Command("SIMPLEX");
2083 
2084  Double_t fmin = 0,fedm = 0,ferrdef = 0;
2085  Int_t npari =0, nparx = 0, istat = 0;
2086  jMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
2087  std::cout << "SIMPLEX " << istat << std::endl;
2088 
2089  jMinuit->Command("HESSE");
2090  jMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
2091  std::cout << "HESSE " << istat << std::endl;
2092 
2093  jMinuit->Command("MIGRAD");
2094  jMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
2095  std::cout << "MIGRAD " << istat << std::endl;
2096 
2097  if(istat == 3){
2098  jMinuit->Command("MINOS");
2099  jMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
2100  std::cout << "MINOS " << istat << std::endl;
2101  }
2102  else{
2103  std::cout << "Fit is messed up" << std::endl;
2104  return badentry;
2105  }
2106 
2107  //Grab the Final Fit Results
2108  bool troubleFitting = false;
2109  double vstart[parametersForFit];
2110  double verror[parametersForFit];
2111  double fParValue;
2112  double fParError;
2113 
2114  for(int parNum = 0; parNum < parametersForFit; parNum++){
2115  jMinuit->GetParameter(parNum,fParValue,fParError);
2116  vstart[parNum] = fParValue;
2117  verror[parNum] = fParError;
2118  if(fParValue > 2 || fParError < 0) troubleFitting = true;
2119  }
2120 
2121  Double_t corr_mat_test[parametersForFit][parametersForFit];
2122  jMinuit->mnemat(*corr_mat_test,parametersForFit);
2123 
2124  std::vector<std::vector<double>> covariance_matrix_pars;
2125  for(int i = 0; i < parametersForFit; i++){
2126  std::vector<double> covHolder;
2127  for(int j =0; j < parametersForFit; j++){
2128  covHolder.push_back(corr_mat_test[i][j]);
2129  }
2130  covariance_matrix_pars.push_back(covHolder);
2131  }
2132 
2133  NueCCIncGlobalFitter::FillParameterCovariance(covariance_matrix_pars);
2134 
2135 
2136  //Test Correlations Between the Normalization Coefficiencts,
2137  //too high of correlations lead to issues in uncertainty estimates
2138  for(uint i = 0; i < numberParameters; i++){
2139  for(uint j = 0; j < numberParameters; j++){
2140  if(i == j) continue;
2141  if( corr_mat_test[i][j]/(pow(corr_mat_test[i][i],0.5)*
2142  pow(corr_mat_test[j][j],0.5))
2143  > 0.99)
2144  troubleFitting = true;
2145  }
2146  }
2147 
2148  //Propagate Normalizations and Uncertainties to the Histograms that
2149  //Hold the results
2150  if(!troubleFitting){
2151  int parNum = 0;
2152  while(parNum < parametersForFit){
2153  std::pair<float,float> analysisBin =
2155  "analysisBin");
2156  std::pair<float,float> minMaxParNum =
2158  "parNum");
2159 
2160  if(minMaxParNum.first < 0 || minMaxParNum.second < 0) continue;
2161 
2162  if(minMaxParNum.second - minMaxParNum.first == 2){
2163  numu_normalization->SetBinContent(analysisBin.first,
2164  analysisBin.second,
2165  vstart[parNum]);
2166  numu_uncertainty->SetBinContent(analysisBin.first,
2167  analysisBin.second,
2168  verror[parNum]);
2169  nue_normalization->SetBinContent(analysisBin.first,
2170  analysisBin.second,
2171  vstart[parNum+1]);
2172  nue_uncertainty->SetBinContent(analysisBin.first,
2173  analysisBin.second,
2174  verror[parNum+1]);
2175  nc_normalization->SetBinContent(analysisBin.first,
2176  analysisBin.second,
2177  vstart[parNum+2]);
2178  nc_uncertainty->SetBinContent(analysisBin.first,
2179  analysisBin.second,
2180  verror[parNum+2]);
2181  parNum += 3;
2182  }
2183  else if(minMaxParNum.second - minMaxParNum.first == 1){
2184  numu_normalization->SetBinContent(analysisBin.first,
2185  analysisBin.second,
2186  vstart[parNum]);
2187  numu_uncertainty->SetBinContent(analysisBin.first,
2188  analysisBin.second,
2189  verror[parNum]);
2190  nue_normalization->SetBinContent(analysisBin.first,
2191  analysisBin.second,
2192  vstart[parNum+1]);
2193  nue_uncertainty->SetBinContent(analysisBin.first,
2194  analysisBin.second,
2195  verror[parNum+1]);
2196  nc_normalization->SetBinContent(analysisBin.first,
2197  analysisBin.second,
2198  vstart[parNum]);
2199  nc_uncertainty->SetBinContent(analysisBin.first,
2200  analysisBin.second,
2201  verror[parNum]);
2202  parNum += 2;
2203  }
2204  };
2205  }
2206  else{
2207  std::cout << "There was an issue with the fit" << std::endl;
2208  return badentry;
2209  }
2210 
2211 
2213  //Now that we have gone over the fit and propagated uncertainty
2214  //Return the final results
2215  std::vector<TH3F*> hResultsPostFit =
2217 
2218  return hResultsPostFit;
2219  }
2220 
2221  std::vector<TH1F*> NueCCIncGlobalFitter::getFullTemplates(bool isReduced)
2222  {
2223  //Order of this is {Data,Signal,NC,Numu,Other};
2224  if(isReduced){
2225  std::vector<TH1F*> hResultsHolder =
2227 
2228  TH1F* hTotalMC =
2229  (TH1F*)hResultsHolder[1]->Clone(ana::UniqueName().c_str());
2230  hTotalMC->Add(hResultsHolder[2]);
2231  hTotalMC->Add(hResultsHolder[3]);
2232  hTotalMC->Add(hResultsHolder[4]);
2233 
2234  return {hResultsHolder[0], hTotalMC, hResultsHolder[1],
2235  hResultsHolder[2], hResultsHolder[3],hResultsHolder[4]};
2236  }
2237  else{
2238  std::vector<TH1F*> hResultsHolder =
2240 
2241  TH1F* hSignalTemplate =
2242  (TH1F*)hResultsHolder[2]->Clone(ana::UniqueName().c_str());
2243  hSignalTemplate->Add(hResultsHolder[3],1);
2244  hSignalTemplate->Add(hResultsHolder[8],1);
2245 
2246  TH1F* hNumuTemplate =
2247  (TH1F*)hResultsHolder[4]->Clone(ana::UniqueName().c_str());
2248  hSignalTemplate->Add(hResultsHolder[5],1);
2249 
2250  TH1F* hNCTemplate =
2251  (TH1F*)hResultsHolder[6]->Clone(ana::UniqueName().c_str());
2252 
2253  TH1F* hOtherTemplate =
2254  (TH1F*)hResultsHolder[7]->Clone(ana::UniqueName().c_str());
2255 
2256 
2257  TH1F* hTotalMC =
2258  (TH1F*)hSignalTemplate->Clone(ana::UniqueName().c_str());
2259  hTotalMC->Add(hNumuTemplate);
2260  hTotalMC->Add(hNCTemplate);
2261  hTotalMC->Add(hOtherTemplate);
2262 
2263  return { hResultsHolder[0],hTotalMC, hSignalTemplate, hNCTemplate,
2264  hNumuTemplate,hOtherTemplate};
2265  }
2266  }
2267 
2269  {
2270  //Order of this is {Data,Signal,NC,Numu,Other};
2271  std::vector<TH1F*> hResultsHolder =
2273 
2274  TH1F* hTotalMC = (TH1F*)hResultsHolder[1]->Clone(ana::UniqueName().c_str());
2275  hTotalMC->Add(hResultsHolder[2]);
2276  hTotalMC->Add(hResultsHolder[3]);
2277  hTotalMC->Add(hResultsHolder[4]);
2278 
2279  return {hResultsHolder[0], hTotalMC, hResultsHolder[1],
2280  hResultsHolder[2], hResultsHolder[3],hResultsHolder[4]};
2281 
2282  }
2283 
2285  {
2286  TH2F* hResult = new TH2F(ana::UniqueName().c_str(),
2287  ";Template Bins;Template Bins",
2288  (int)vCovarianceMatrixFull.size(),
2289  0,(int)vCovarianceMatrixFull.size(),
2290  (int)vCovarianceMatrixFull.size(),
2291  0,(int)vCovarianceMatrixFull.size());
2292 
2293  for(uint i = 0; i < vCovarianceMatrixFull.size(); i++){
2294  for(uint j = 0; j < vCovarianceMatrixFull.size(); j++){
2295  hResult->SetBinContent(i+1,j+1,vCovarianceMatrixFull[i][j]);
2296  }
2297  }
2298  return hResult;
2299  }
2300 
2302  {
2303  TH2F* hResult = new TH2F(ana::UniqueName().c_str(),
2304  ";Template Bins;Template Bins",
2305  (int)vCovarianceMatrixReduced.size(),
2306  0,(int)vCovarianceMatrixReduced.size(),
2307  (int)vCovarianceMatrixReduced.size(),
2308  0,(int)vCovarianceMatrixReduced.size());
2309 
2310  for(uint i = 0; i < vCovarianceMatrixReduced.size(); i++){
2311  for(uint j = 0; j < vCovarianceMatrixReduced.size(); j++){
2312  hResult->SetBinContent(i+1,j+1,vCovarianceMatrixReduced[i][j]);
2313  }
2314  }
2315  return hResult;
2316  }
2317 
2319  {
2320  TH2F* hResult = new TH2F(ana::UniqueName().c_str(),
2321  ";Template Bins;Template Bins",
2322  (int)vCovarianceMatrixFull.size(),
2323  0,(int)vCovarianceMatrixFull.size(),
2324  (int)vCovarianceMatrixFull.size(),
2325  0,(int)vCovarianceMatrixFull.size());
2326 
2327  for(uint i = 0; i < vCovarianceMatrixFull.size(); i++){
2328  for(uint j = 0; j < vCovarianceMatrixFull.size(); j++){
2329  hResult->SetBinContent(i+1,j+1,
2331  }
2332  }
2333  return hResult;
2334  }
2335 
2337  {
2338  TH2F* hResult = new TH2F(ana::UniqueName().c_str(),
2339  ";Template Bins;Template Bins",
2340  (int)vCovarianceMatrixReduced.size(),
2341  0,(int)vCovarianceMatrixReduced.size(),
2342  (int)vCovarianceMatrixReduced.size(),
2343  0,(int)vCovarianceMatrixReduced.size());
2344 
2345  for(uint i = 0; i < vCovarianceMatrixReduced.size(); i++){
2346  for(uint j = 0; j < vCovarianceMatrixReduced.size(); j++){
2347  hResult->SetBinContent(i+1,j+1,
2349  }
2350  }
2351  return hResult;
2352  }
2353 
2354 
2356  {
2357  TH2F* hResult = new TH2F(ana::UniqueName().c_str(),
2358  ";Parameter Number;Parameter Number;",
2359  (int)parCovariance.size(),
2360  0, (int)parCovariance.size(),
2361  (int)parCovariance.size(),
2362  0, (int)parCovariance.size());
2363 
2364  for(uint i = 0; i < parCovariance.size(); i++){
2365  for(uint j = 0; j < parCovariance.size(); j++){
2366  hResult->SetBinContent(i+1,j+1,parCovariance[i][j]);
2367  }
2368  }
2369 
2370  return hResult;
2371  }
2372 
2373 
2375  {
2376  TH2F* hResult = new TH2F(ana::UniqueName().c_str(),
2377  ";Parameter Number;Parameter Number;",
2378  (int)parCorrelation.size(),
2379  0, (int)parCorrelation.size(),
2380  (int)parCorrelation.size(),
2381  0, (int)parCorrelation.size());
2382 
2383  for(uint i = 0; i < parCorrelation.size(); i++){
2384  for(uint j = 0; j < parCorrelation.size(); j++){
2385  hResult->SetBinContent(i+1,j+1,parCorrelation[i][j]);
2386  }
2387  }
2388 
2389  return hResult;
2390  }
2391 
2392  std::vector<std::pair<TH2F*,TH2F*>>
2394  {
2395 
2396  std::vector<std::pair<TH2F*,TH2F*>> vResult;
2397 
2398  TH3F* hHolder3D = (TH3F*) hAnalysis3D[0]->Clone(ana::UniqueName().c_str());
2399  TH2F* hHolder2D = (TH2F*)hHolder3D->Project3D("yx");
2400  hHolder2D->SetName(ana::UniqueName().c_str());
2401 
2402  TH2F* hNueWei2D = (TH2F*)hHolder2D->Clone(ana::UniqueName().c_str());
2403  TH2F* hNueUnc2D = (TH2F*)hHolder2D->Clone(ana::UniqueName().c_str());
2404  TH2F* hNumuWei2D = (TH2F*)hHolder2D->Clone(ana::UniqueName().c_str());
2405  TH2F* hNumuUnc2D = (TH2F*)hHolder2D->Clone(ana::UniqueName().c_str());
2406  TH2F* hNCWei2D = (TH2F*)hHolder2D->Clone(ana::UniqueName().c_str());
2407  TH2F* hNCUnc2D = (TH2F*)hHolder2D->Clone(ana::UniqueName().c_str());
2408 
2409  for(int i=1; i<=hHolder2D->GetXaxis()->GetNbins();i++){
2410  for(int j=1;j<=hHolder2D->GetYaxis()->GetNbins();j++){
2411  float nue_norm = nue_normalization->GetBinContent(i,j);
2412  float nue_err = nue_uncertainty->GetBinContent(i,j);
2413  float numu_norm = numu_normalization->GetBinContent(i,j);
2414  float numu_err = numu_uncertainty->GetBinContent(i,j);
2415  float nc_norm = nc_normalization->GetBinContent(i,j);
2416  float nc_err = nc_uncertainty->GetBinContent(i,j);
2417  hNueWei2D->SetBinContent(i,j,nue_norm);
2418  hNumuWei2D->SetBinContent(i,j,numu_norm);
2419  hNCWei2D->SetBinContent(i,j,nc_norm);
2420  hNueUnc2D->SetBinContent(i,j,nue_err);
2421  hNumuUnc2D->SetBinContent(i,j,numu_err);
2422  hNCUnc2D->SetBinContent(i,j,nc_err);
2423  }
2424  }
2425 
2426 
2427  vResult.push_back(std::make_pair(hNueWei2D,hNueUnc2D));
2428  vResult.push_back(std::make_pair(hNumuWei2D,hNumuUnc2D));
2429  vResult.push_back(std::make_pair(hNCWei2D,hNCUnc2D));
2430 
2431  return vResult;
2432  }
2433 
2434 }//end namespace nueccinc_test
std::vector< TH1F * > GetTemplatesReducedPostFit()
void FillCovarianceMatrixReduced(bool makePlots=false)
ofstream output
fvar< T > fmin(const fvar< T > &x1, const fvar< T > &x2)
Definition: fmin.hpp:14
double GetStatisticMatrixReducedValue(int binx, int biny)
enum BeamMode kRed
std::vector< std::pair< TH2F *, TH2F * > > GetTemplateWeightsAndErrors()
std::vector< std::vector< double > > parCorrelation
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
Calculate bin to bin correlation matrix from the universe technique.
std::vector< std::vector< float > > vCovarianceMatrixFull
TH2D * GetCovarianceMatrix()
Get covariance matrix.
void makePlots()
Definition: makePlots.C:17
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
Int_t par
Definition: SimpleIterate.C:24
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
OStream cerr
Definition: OStream.cxx:7
void FillParameterCovariance(std::vector< std::vector< double >> cov)
std::vector< std::vector< float > > vInvertedCovarianceMatrixReduced
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
std::vector< std::vector< double > > parCovariance
double myFunction(std::vector< double > par)
std::vector< std::vector< float > > vStatCovarianceMatrixFull
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
const XML_Char const XML_Char * data
Definition: expat.h:268
Double_t ymax
Definition: plot.C:25
static void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
const int nbins
Definition: cellShifts.C:15
std::vector< std::vector< float > > vStatCovarianceMatrixReduced
const XML_Char int const XML_Char * value
Definition: expat.h:331
TH2D * GetCorrelationMatrix()
Get correlation matrix.
correl_xv GetXaxis() -> SetDecimals()
const std::string cv[Ncv]
std::pair< float, float > searchFitInfo(float parNum, std::string fitSearch)
double GetCovarianceMatrixFullValue(int binx, int biny)
void FillCovarianceMatrixFull(bool makePlots=false)
std::vector< TH1F * > CalculateOneSigmaShiftReduced()
const double j
Definition: BetheBloch.cxx:29
const std::string covarianceDir
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
NueCCIncGlobalFitter(TH1F *da, std::vector< TH1F * > templates, std::vector< TH3F * > analysis, TH1F *cv, std::vector< TH1F * > systs_up, std::vector< TH1F * > systs_down, std::vector< TH1F * > systs_oneside, std::vector< TH1F * > multiverse)
float bin[41]
Definition: plottest35.C:14
OStream cout
Definition: OStream.cxx:6
cout<< "--"<< endl;for(Int_t iP=1;iP<=hyz->GetNbinsX();iP++){for(Int_t iC=1;iC<=hyz->GetNbinsY();iC++){if(hyv->GetBinContent(iP, iC)>-999){goal_hyv-> SetBinContent(iP, iC,-(dy[iP-1][iC-1]))
const int numberParameters
NueCCIncGlobalFitter * Fitter_obj
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
enum BeamMode nc
std::vector< TH1F * > CalculateOneSigmaShiftFull()
c1
Definition: demo5.py:24
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
std::vector< TH1F * > getFullTemplates(bool isReduced)
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
std::vector< std::vector< float > > vCovarianceMatrixReduced
double GetStatisticMatrixFullValue(int binx, int biny)
enum BeamMode kBlue
double GetCovarianceMatrixReducedValue(int binx, int biny)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
enum BeamMode string
unsigned int uint