NueCCIncTemplateFitter.cxx
Go to the documentation of this file.
1 
3 #include <functional>
4 #include <string>
5 #include <algorithm>
7 
8 namespace ana{
9  namespace nueccinc
10  {
11  //Defines where the signal region begins in the template being
12  //used for the fit
13  const float signal_region = 0.8; //0.8 for CVN,
14  //0.9 for LID,
15  // 0.2 for ElectronID
16  const float sliding_signal_region[13] = {0.2,0.2,0.15,0.15,0.15,
17  0.25,0.25,0.25,0.25,0.25,0.25,0.25,
18  0.25};
19 
20  const bool useSlidingSignalRegion = true;
21  const bool printVerbose = false;
22 
23  //Where to Save the individual covariance/correlation matrices
24  const std::string covarianceDir = "Covariance/";
25 
26 
28  (TH3F* da,std::vector<TH3F*> templates,std::vector<TH3F*> analysis,
29  TH3F* cv,std::vector<TH3F*> systs_up,
30  std::vector<TH3F*> systs_down,std::vector<TH3F*> systs_oneside,
31  std::vector<TH3F*> multiverse,int numbins)
32  :hData3D(da),
33  hTemplates3D(templates),
34  hAnalysis3D(analysis),
35  hCV3D(cv),
36  hSystsUp3D(systs_up),
37  hSystsDown3D(systs_down),
38  hSystsOneSided(systs_oneside),
39  hSystsMultiverse(multiverse),
40  NumberTemplateBins(numbins)
41  {
42  hSignalWeights2D = (TH2F*)hCV3D->Project3D("yx");
43  hSignalWeights2D->SetName("hSignalWeights2D");
44  hNumuCCWeights2D = (TH2F*)hCV3D->Project3D("yx");
45  hNumuCCWeights2D->SetName("hNumuCCWeights2D");
46  hNCWeights2D = (TH2F*)hCV3D->Project3D("yx");
47  hNCWeights2D->SetName("hNCWeights2D");
48  hSignalError2D = (TH2F*)hCV3D->Project3D("yx");
49  hSignalError2D->SetName("hSignalError2D");
50  hNumuCCError2D = (TH2F*)hCV3D->Project3D("yx");
51  hNumuCCError2D->SetName("hNumuCCError2D");
52  hNCError2D = (TH2F*)hCV3D->Project3D("yx");
53  hNCError2D->SetName("hNCError2D");
54  hChiSquaredResults2D = (TH2F*)hCV3D->Project3D("yx");
55  hChiSquaredResults2D->SetName("hChiSquaredResults");
56 
57  nue_wei = std::make_pair(0,0);
58  numu_wei = std::make_pair(0,0);
59  nc_wei = std::make_pair(0,0);
60  chisquared = 0;
61  }
62 
64 
65  //////////// //////////// //////////// //////////// ////////////
66  //Protected Functions///////////////////////////////////////////////
67  //////////// //////////// //////////// //////////// ////////////
68 
70  {
71 
72  for(int i = 0; i < NumberTemplateBins; i++){
73  DataValues[i] = ((double)hData3D->GetBinContent(binx,biny,i+1));
75  (double)hTemplates3D[1]->GetBinContent(binx,
76  biny,i+1)
77  + (double)hTemplates3D[2]->GetBinContent(binx,
78  biny,
79  i+1)
80  + (double)hTemplates3D[7]->GetBinContent(binx,
81  biny,
82  i+1);
83  NumuLike_Values[i] =
84  (double)hTemplates3D[3]->GetBinContent(binx,
85  biny,i+1)
86  + (double)hTemplates3D[4]->GetBinContent(binx,
87  biny,
88  i+1);
89  NC_Values[i] =
90  (double)hTemplates3D[5]->GetBinContent(binx,biny,i+1);
91  Other_Values[i] =
92  (double)hTemplates3D[6]->GetBinContent(binx,biny,i+1);
93  }//loop over template bins
94  }
95 
97  {
98  return ((hSystsUp3D.size() == hSystsDown3D.size()) &&
99  (hData3D->GetXaxis()->GetNbins() ==
100  hCV3D->GetXaxis()->GetNbins())
101  &&
102  (hData3D->GetYaxis()->GetNbins() ==
103  hCV3D->GetYaxis()->GetNbins())
104  &&
105  (hData3D->GetZaxis()->GetNbins() ==
106  hCV3D->GetZaxis()->GetNbins()));
107  }
108 
110  {
111  std::vector<TH3F*> hResults;
112  for(unsigned int i = 0; i < hAnalysis3D.size(); i++)
113  hResults.push_back((TH3F*)
114  hAnalysis3D[i]->Clone
115  (ana::UniqueName().c_str()));
116  return hResults;
117  }
118 
120  int biny)
121  {
122  std::stringstream low_X;
123  low_X << std::fixed << std::setprecision(2) <<
124  hCV3D->GetXaxis()->GetBinLowEdge(binx);
125  std::stringstream hi_X;
126  hi_X << std::fixed << std::setprecision(2) <<
127  hCV3D->GetXaxis()->GetBinUpEdge(binx);
128 
129  std::stringstream low_Y;
130  low_Y << std::fixed << std::setprecision(2) <<
131  hCV3D->GetYaxis()->GetBinLowEdge(biny);
132  std::stringstream hi_Y;
133  hi_Y << std::fixed << std::setprecision(2) <<
134  hCV3D->GetYaxis()->GetBinUpEdge(biny);
135 
136  std::string caption = low_X.str() + " #leq cos #theta_{e} < "+
137  hi_X.str() + " " + low_Y.str() + " #leq E_{e} < " + hi_Y.str();
138  return caption;
139  }
140 
142  {
143  //Subtract one to compare with vector entries
144  int signal_region_bin = 0;
145  if(useSlidingSignalRegion)
146  signal_region_bin =
147  hCV3D->GetZaxis()->FindBin(sliding_signal_region[biny])-1;
148  else
149  signal_region_bin =
150  hCV3D->GetZaxis()->FindBin(signal_region) - 1;
151  std::cout << signal_region_bin << std::endl;
152 
153  float signal_amount = 0;
154  float bkgd_amount = 0;
155  for(int bin = signal_region_bin; bin < NumberTemplateBins; bin++)
156  {
157  signal_amount += SignalLike_Values[bin];
158  bkgd_amount += (NC_Values[bin]+NumuLike_Values[bin]+Other_Values[bin]);
159  }
160 
161  return ( ((signal_amount > 100) && (signal_amount/bkgd_amount > 0.25)) ||
162  ((signal_amount/bkgd_amount > 0.2) && (signal_amount > 300)));
163  }
164 
166  int biny)
167  {
168  //Function to calculate shift from one-sided and two sided systematics
169  std::vector<TH1F*> hResults;
170 
171 
172  TH1F* hCV =
173  (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
174  binx,binx,biny,biny);
175 
176  //One Sided Systematics
177  for(int i = 0; i < (int)hSystsOneSided.size();i++){
178  TH1F* hHolder =
179  (TH1F*)hSystsOneSided[i]->ProjectionZ(ana::UniqueName().c_str(),
180  binx,binx,biny,biny);
181  hHolder->Add(hCV,-1);
182  hResults.push_back(hHolder);
183  }
184 
185  //Two Sided Systematics
186  for(int i = 0; i < (int)hSystsUp3D.size(); i++){
187  TH1F* hHolderUp =
188  (TH1F*)hSystsUp3D[i]->ProjectionZ(ana::UniqueName().c_str(),
189  binx,binx,biny,biny);
190  TH1F* hHolderDwn =
191  (TH1F*)hSystsDown3D[i]->ProjectionZ(ana::UniqueName().c_str(),
192  binx,binx,biny,biny);
193 
194  hHolderUp->Add(hCV,-1);
195  hHolderDwn->Add(hCV,-1);
196 
197  for(int bin = 0; bin <= hHolderUp->GetXaxis()->GetNbins();bin++){
198  if(fabs(hHolderDwn->GetBinContent(bin)) >
199  fabs(hHolderUp->GetBinContent(bin)))
200  hHolderUp->SetBinContent(bin,hHolderDwn->GetBinContent(bin));
201  }
202  hResults.push_back((TH1F*)hHolderUp->Clone(ana::UniqueName().c_str()));
203  }
204 
205  return hResults;
206  }
207 
209  bool makePlots = false)
210  {
211  //Make Sure this bin isn't excluded by the phase space cuts
212  //If it is, don't worrky about calculating anything
213  //Just checking one of the inner bins
214  if(hCV3D->Integral(binx,binx,biny,biny,1,hCV3D->GetZaxis()->GetNbins()) ==
215  0) return;
216 
217  std::string binTitle =
219 
220 
221  //Systematic Samples
222  //Turn Systamtic samples that are not generated from a multiverse
223  //into a "one sigma" shift that we can randomly draw from in a moment
224  std::vector<TH1F*> hSystsShifts =
226 
227  std::map<std::string,TH2F*> mapCovariances;
228  std::map<std::string,TH2F*> mapCorrelations;
229  //Order of the systematics is Ckv, Calibration Shape, Light Levels, Cal
230  std::vector<std::string> vSystNames = {"ckv", "calib_shape", "light",
231  "calib", //"horn_current",
232  //"spot_size", "shiftx",
233  //"shifty", "horn1x", "horn1y",
234  //"horn2x", "horn2y", "targetz",
235  //"mag_field", "hornwater",
236  "multiverse"};
237 
238  //TRandom *r0 = new TRandom();
239 
240  int singlesided_counter = 0;
241  int doublesided_counter = 0;
242  for(uint systNum = 0; systNum < hSystsShifts.size(); systNum++){
243  //First draw systematic around the cv for the template
244  if((vSystNames[systNum] == "ckv" ||
245  vSystNames[systNum] == "calib_shape") && makePlots){
246  hSystsOneSided[singlesided_counter]->SetLineColorAlpha(kRed-9,0.5);
247  TH1F* hSyst =
248  (TH1F*)hSystsOneSided[singlesided_counter]->
249  ProjectionZ(ana::UniqueName().c_str(),
250  binx,binx,biny,biny);
251  TH1F* hCV1D =
252  (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
253  binx,binx,biny,biny);
254  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
255  hCV1D->GetXaxis()->SetTitle("Template Bins");
256  hCV1D->SetTitle(binTitle.c_str());
257  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
258  c1->SetLeftMargin(0.15);
259  c1->SetBottomMargin(0.15);
260  hCV1D->GetYaxis()->SetRangeUser(0,hCV1D->GetMaximum()*1.5);
261  hCV1D->Draw("hist");
262  hSyst->Draw("hist same");
263  TLegend* leg = ana::AutoPlaceLegend(0.3,0.4);
264  leg->AddEntry(hSyst, "Sigma Shift", "fl");
265  leg->AddEntry(hCV1D, "Central Value", "fl");
266  leg->Draw();
267  c1->SaveAs(Form("%s%s_%s_%i_%i_only.png",covarianceDir.c_str(),
268  "standard_template", vSystNames[systNum].c_str(),
269  binx,biny));
270  c1->Close();
271  delete leg;
272  delete c1;
273  singlesided_counter++;
274  }
275  else if(makePlots){
276  hSystsUp3D[doublesided_counter]->SetLineColorAlpha(kRed-9,0.5);
277  hSystsDown3D[doublesided_counter]->SetLineColorAlpha(kBlue-9,0.5);
278  TH1F* hSystsUp =
279  (TH1F*)hSystsUp3D[doublesided_counter]->
280  ProjectionZ(ana::UniqueName().c_str(),
281  binx,binx,biny,biny);
282  TH1F* hSystsDown =
283  (TH1F*)hSystsDown3D[doublesided_counter]->
284  ProjectionZ(ana::UniqueName().c_str(),
285  binx,binx,biny,biny);
286  TH1F* hCV1D =
287  (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
288  binx,binx,biny,biny);
289  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
290  hCV1D->GetXaxis()->SetTitle("Template Bins");
291  hCV1D->SetTitle(binTitle.c_str());
292  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
293  c1->SetLeftMargin(0.15);
294  c1->SetBottomMargin(0.15);
295  hCV1D->GetYaxis()->SetRangeUser(0,hCV1D->GetMaximum()*1.5);
296  hCV1D->Draw("hist");
297  hSystsUp->Draw("hist same");
298  hSystsDown->Draw("hist same");
299  TLegend* leg = ana::AutoPlaceLegend(0.3,0.4);
300  leg->AddEntry(hSystsUp, "Sigma Up", "fl");
301  leg->AddEntry(hSystsDown, "Sigma Down", "fl");
302  leg->AddEntry(hCV1D, "Central Value", "fl");
303  leg->Draw();
304  c1->SaveAs(Form("%s%s_%s_%i_%i_only.png",covarianceDir.c_str(),
305  "standard_template", vSystNames[systNum].c_str(),
306  binx,biny));
307  c1->Close();
308  delete leg;
309  delete c1;
310  doublesided_counter++;
311  }
312 
313  //Randomly draw from the one sigma shift samples
314  //Draw 1000 from each sample
315  std::vector<TH1F*> vPseudoExp;
316  for(uint i = 0; i < 1; i++){
317  TH1F* hClone = (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
318  binx,binx,biny,biny);
319  double wei = 1;
320  hClone->Add(hSystsShifts[systNum],wei);
321  vPseudoExp.push_back(hClone);
322  }//loop over pseudo experiments
323 
324 
325  TH1F* hCV1D = (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
326  binx,binx,biny,biny);
327  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
328  hCV1D->GetXaxis()->SetTitle("Template Bins");
329  hCV1D->GetYaxis()->CenterTitle();
330  hCV1D->GetXaxis()->CenterTitle();
331  hCV1D->SetTitle(binTitle.c_str());
332 
333 
334  //Draw the Pseudo experiments
335  if(makePlots){
336  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
337  hCV1D->GetYaxis()->SetRangeUser(0,hCV1D->GetMaximum()*1.5);
338  hCV1D->Draw("hist");
339  for(uint i = 0; i < vPseudoExp.size(); i++){
340  vPseudoExp[i]->SetLineColorAlpha(kBlue-9,0.50);
341  vPseudoExp[i]->Draw("hist same");
342  }
343  hCV1D->Draw("hist same");
344  c1->SaveAs(Form("%s%s_%s_%i_%i_only.png","Covariance/",
345  "standard_template_pseudo_exp",
346  vSystNames[systNum].c_str(),binx,biny));
347  c1->Close();
348  }
349 
350  TH2F* hCovHolder =
351  new TH2F(ana::UniqueName().c_str(),
352  ";Template Bins;Template Bins; Covariance (Events^{2})",
355  TH2F* hCorHolder =
356  new TH2F(ana::UniqueName().c_str(),
357  ";Template Bins;Template Bins; Correlation",
360  hCovHolder->SetTitle(binTitle.c_str());
361  hCorHolder->SetTitle(binTitle.c_str());
362 
363  //Calculate Covariance Matrix For the Single Systematic
364  for(uint iu = 0; iu < vPseudoExp.size(); iu++){
365  for(int i = 0; i < numPIDBins; i++){
366  double xi=vPseudoExp[iu]->GetBinContent(i+1);
367  double ximean=hCV1D->GetBinContent(i+1);
368  for(int j = i; j < numPIDBins; j++){
369  double xj=vPseudoExp[iu]->GetBinContent(j+1);
370  double xjmean=hCV1D->GetBinContent(j+1);
371  double current_value = hCovHolder->GetBinContent(i+1,j+1);
372  double this_value = (xi-ximean)*(xj-xjmean);
373  hCovHolder->SetBinContent(i+1,j+1,current_value+this_value);
374  }
375  }
376 
377  }
378 
379  //Only calculate with the 1 sigma error band
380  const double N=(double)2;
381  hCovHolder->Scale(1.0/(N-1.0));
382  for(int i=0; i<numPIDBins; i++){
383  for(int j=i; j<numPIDBins; j++){
384  hCovHolder->SetBinContent(j+1,i+1,
385  hCovHolder->GetBinContent(i+1,j+1));
386  }
387  }
388 
389  for(int i=0; i<numPIDBins; i++){
390  for(int j=0; j<numPIDBins; j++){
391  double covariance = hCovHolder->GetBinContent(i+1,j+1);
392  double variance_i = sqrt(hCovHolder->GetBinContent(i+1,i+1));
393  double variance_j = sqrt(hCovHolder->GetBinContent(j+1,j+1));
394  double correlation = covariance/(variance_i*variance_j);
395  if(isnan(correlation)) correlation = 0;
396  hCorHolder->SetBinContent(i+1,j+1,
397  correlation);
398  }
399  }
400 
401 
402  mapCovariances.insert(std::make_pair(vSystNames[systNum],
403  hCovHolder));
404  mapCorrelations.insert(std::make_pair(vSystNames[systNum],
405  hCorHolder));
406  }
407 
408 
409  //Now For the Multiverse
410  //GENIE and Shape Only PPFX
411  {
412 
413  TH2F* hCovHolder =
414  new TH2F(ana::UniqueName().c_str(),
415  ";Template Bins;Template Bins; Covariance (Events^{2})",
418  TH2F* hCorHolder =
419  new TH2F(ana::UniqueName().c_str(),
420  ";Template Bins;Template Bins;Correlation",
423 
424  hCovHolder->SetTitle(binTitle.c_str());
425  hCorHolder->SetTitle(binTitle.c_str());
426 
427  TH2F* hCovHolderPPFX =
428  new TH2F(ana::UniqueName().c_str(),
429  ";Template Bins;Template Bins; Covariance (Events^{2})",
432  TH2F* hCorHolderPPFX =
433  new TH2F(ana::UniqueName().c_str(),
434  ";Template Bins;Template Bins;Correlation",
437 
438  hCovHolderPPFX->SetTitle(binTitle.c_str());
439  hCorHolderPPFX->SetTitle(binTitle.c_str());
440 
441  TH1D* hCV_Tester = (TH1D*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
442  binx,binx,biny,biny);
443 
444  std::vector<TH1D*> hMultiverseHistsGENIE;
445  std::vector<TH1D*> hMultiverseHistsPPFX;
446 
447  //Let's try to area normalize these
448  for(uint i = 0; i < hSystsMultiverse.size()- 100; i++){
449  TH1D* hHolder=
450  (TH1D*)hSystsMultiverse[i]->ProjectionZ
451  (ana::UniqueName().c_str(),binx,binx,biny,biny);
452  float ymax = hHolder->Integral();
453  float cvmax = hCV_Tester->Integral();
454  hHolder->Scale(cvmax/ymax);
455  hMultiverseHistsGENIE.push_back(hHolder);
456  }
457 
458  //PPFX Universes
459  for(uint i = hSystsMultiverse.size()-100;
460  i < hSystsMultiverse.size(); i++){
461  TH1D* hHolder=
462  (TH1D*)hSystsMultiverse[i]->ProjectionZ
463  (ana::UniqueName().c_str(),binx,binx,biny,biny);
464  float ymax = hHolder->Integral();
465  float cvmax = hCV_Tester->Integral();
466  hHolder->Scale(cvmax/ymax);
467  hMultiverseHistsPPFX.push_back(hHolder);
468  }
469 
470  ana::MultiverseCorrelation* multi_correlation_calc_genie =
471  new ana::MultiverseCorrelation(hMultiverseHistsGENIE);
472 
473  TH2F* hCov =
474  (TH2F*)multi_correlation_calc_genie->GetCovarianceMatrix();
475  TH2F* hCor =
476  (TH2F*)multi_correlation_calc_genie->GetCorrelationMatrix();
477 
478  ana::MultiverseCorrelation* multi_correlation_calc_ppfx =
479  new ana::MultiverseCorrelation(hMultiverseHistsPPFX);
480 
481  TH2F* hCovPPFX =
482  (TH2F*)multi_correlation_calc_ppfx->GetCovarianceMatrix();
483  TH2F* hCorPPFX =
484  (TH2F*)multi_correlation_calc_ppfx->GetCorrelationMatrix();
485 
486  for(int i = 1; i <= hCov->GetXaxis()->GetNbins(); i++){
487  for(int j = 1; j <= hCov->GetXaxis()->GetNbins(); j++){
488  hCovHolder->SetBinContent(i,j,hCov->GetBinContent(i,j));
489  hCorHolder->SetBinContent(i,j,hCor->GetBinContent(i,j));
490  hCovHolderPPFX->SetBinContent(i,j,hCovPPFX->GetBinContent(i,j));
491  hCorHolderPPFX->SetBinContent(i,j,hCorPPFX->GetBinContent(i,j));
492  }
493  }
494 
495 
496  if(makePlots){
497  TH1F* hCV1D = (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
498  binx,binx,biny,biny);
499  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
500  hCV1D->GetYaxis()->SetTitle("Template Bins");
501  hCV1D->GetYaxis()->CenterTitle();
502  hCV1D->GetXaxis()->CenterTitle();
503  hCV1D->SetTitle(binTitle.c_str());
504 
505  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
506  hCV1D->Draw("hist");
507  for(uint i = 0; i < hMultiverseHistsPPFX.size(); i++){
508  hMultiverseHistsPPFX[i]->SetLineColorAlpha(kBlue-9,0.50);
509  hMultiverseHistsPPFX[i]->Draw("hist same");
510  }
511  hCV1D->Draw("hist same");
512  c1->SaveAs(Form("%s%s_%s_%i_%i_only.png","Covariance/",
513  "standard_template_pseudo_exp",
514  "ppfx",binx,biny));
515  c1->Close();
516  }
517 
518  if(makePlots){
519  TH1F* hCV1D = (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
520  binx,binx,biny,biny);
521  hCV1D->GetYaxis()->SetTitle("Events/8.09 #times 10^{20} (POT)");
522  hCV1D->GetXaxis()->SetTitle("Template Bins");
523  hCV1D->SetTitle(binTitle.c_str());
524  hCV1D->GetYaxis()->CenterTitle();
525  hCV1D->GetXaxis()->CenterTitle();
526  TCanvas *c1 = new TCanvas(ana::UniqueName().c_str());
527  hCV1D->Draw("hist");
528  for(uint i = 0; i < hMultiverseHistsGENIE.size()-1;i++){
529 
530  hMultiverseHistsGENIE[i]->SetLineColorAlpha(kBlue-9,0.50);
531  hMultiverseHistsGENIE[i]->Draw("hist same");
532  }
533  hCV1D->Draw("hist same");
534  c1->SaveAs(Form("%s%s_%s_%i_%i_only.png","Covariance/",
535  "standard_template_pseudo_exp",
536  "genie",binx,biny));
537  c1->Close();
538  }
539 
540 
541  mapCovariances.insert(std::make_pair("genie",
542  hCovHolder));
543  mapCorrelations.insert(std::make_pair("genie",
544  hCorHolder));
545  mapCovariances.insert(std::make_pair("ppfx",
546  hCovHolderPPFX));
547  mapCorrelations.insert(std::make_pair("ppfx",
548  hCorHolderPPFX));
549  delete hCorPPFX;
550  delete hCovPPFX;
551  delete multi_correlation_calc_ppfx;
552  for(auto ele :hMultiverseHistsPPFX)
553  delete ele;
554  delete hCor;
555  delete hCov;
556  for(auto ele :hMultiverseHistsGENIE)
557  delete ele;
558  delete multi_correlation_calc_genie;
559 
560  }
561 
562  //Now Plot All Covariance Matrices Seperately
563  if(makePlots){
564  for(auto pair: mapCovariances){
565  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
566  pair.second->GetXaxis()->CenterTitle();
567  pair.second->GetYaxis()->CenterTitle();
568  pair.second->GetZaxis()->CenterTitle();
569  pair.second->GetZaxis()->SetTitle("Covariance (Events^{2})");
570  pair.second->SetTitle(binTitle.c_str());
571  c3->SetLeftMargin(0.1);
572  c3->SetRightMargin(0.15);
573  c3->SetBottomMargin(0.1);
574  pair.second->GetZaxis()->SetMaxDigits(3);
575  pair.second->Draw("COLZ");
576  c3->SaveAs(Form("%s%s_%s_%i_%i_only.png","Covariance/",
577  "covariance_single",
578  pair.first.c_str(),binx,biny));
579  c3->Close();
580  delete c3;
581  }
582  }//if makePlots
583 
584  if(makePlots){
585  for(auto pair: mapCorrelations){
586  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
587  pair.second->GetXaxis()->CenterTitle();
588  pair.second->GetYaxis()->CenterTitle();
589  pair.second->GetZaxis()->CenterTitle();
590  pair.second->SetTitle(binTitle.c_str());
591  pair.second->GetZaxis()->SetTitle("Correlation");
592  pair.second->GetZaxis()->SetRangeUser(-1,1);
593  c3->SetLeftMargin(0.1);
594  c3->SetRightMargin(0.15);
595  c3->SetBottomMargin(0.1);
596  pair.second->GetZaxis()->SetMaxDigits(3);
597  pair.second->Draw("COLZ");
598  c3->SaveAs(Form("%s%s_%s_%i_%i_only.png",covarianceDir.c_str(),
599  "correlation_single",
600  pair.first.c_str(),binx,biny));
601  c3->Close();
602  delete c3;
603  }
604  }//if makePlots
605 
606  //Now start adding up each of the individual covariance matrices
607  TH2F* hCovariance =
608  (TH2F*)mapCovariances["genie"]->Clone(ana::UniqueName().c_str());
609 
610 
611 
612 
613  std::string output = covarianceDir + "covariance_addition";
614  for(auto pair: mapCovariances){
615  if(pair.first == "genie")continue;
616  output+="_";
617  output+=pair.first;
618  hCovariance->Add(pair.second);
619  if(makePlots){
620  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
621  c3->SetLeftMargin(0.1);
622  c3->SetRightMargin(0.15);
623  c3->SetBottomMargin(0.1);
624  hCovariance->SetTitle(binTitle.c_str());
625  hCovariance->GetZaxis()->SetMaxDigits(3);
626  hCovariance->GetZaxis()->SetTitle("Covariance (Events^{2})");
627  hCovariance->Draw("COLZ");
628  c3->SaveAs(Form("%s_%i_%i.png",output.c_str(),binx,biny));
629  c3->Close();
630  delete c3;
631  }
632  }
633 
634  TH2F* hCorrelation =
635  (TH2F*)hCovariance->Clone(ana::UniqueName().c_str());
636 
637  for(int i=0; i<numPIDBins; i++){
638  for(int j=0; j<numPIDBins; j++){
639  double covariance = hCovariance->GetBinContent(i+1,j+1);
640  double variance_i = sqrt(hCovariance->GetBinContent(i+1,i+1));
641  double variance_j = sqrt(hCovariance->GetBinContent(j+1,j+1));
642  double correlation = covariance/(variance_i*variance_j);
643  if(isnan(correlation)) correlation = 0;
644  hCorrelation->SetBinContent(i+1,j+1,
645  correlation);
646  }
647  }
648 
649  if(makePlots){
650  TCanvas* c3 = new TCanvas(ana::UniqueName().c_str());
651  hCorrelation->GetXaxis()->CenterTitle();
652  hCorrelation->GetYaxis()->CenterTitle();
653  hCorrelation->GetZaxis()->CenterTitle();
654  hCorrelation->SetTitle(binTitle.c_str());
655  hCorrelation->GetZaxis()->SetTitle("Correlation");
656  hCorrelation->GetZaxis()->SetRangeUser(-1,1);
657  c3->SetLeftMargin(0.1);
658  c3->SetRightMargin(0.1);
659  c3->SetBottomMargin(0.1);
660  hCorrelation->GetYaxis()->SetMaxDigits(3);
661  hCorrelation->Draw("COLZ");
662  c3->SaveAs(Form("%s%s_%i_%i.png",covarianceDir.c_str(),
663  "correlation_all",binx,biny));
664  c3->Close();
665  delete c3;
666  }
667 
668 
669  //Move the covariance matrix to the holder so that it can be used
670  //in the fit
671  for(int i = 0; i < NumberTemplateBins; i++){
672  for(int j = 0; j < NumberTemplateBins; j++){
673  CovarianceMatrix[i][j] = hCovariance->GetBinContent(i+1,j+1);
674  if((i == j) && hCovariance->GetBinContent(i+1,j+1) < 1)
675  CovarianceMatrix[i][j] = 10;
676  }
677  }
678 
679  for(uint i = 0; i < hSystsShifts.size(); i++){
680  delete hSystsShifts[i];
681  }
682  hSystsShifts.clear();
683  }
684 
686  int biny)
687  {
688  for(int i = 1; i <= NumberTemplateBins; i++){
689  for(int j = 1; j <= NumberTemplateBins; j++){
690  if(i == j){
691  StatCovarianceMatrix[i-1][j-1] =
692  pow(hData3D->GetBinError(binx,biny,i),2) +
693  pow(hCV3D->GetBinError(binx,biny,i),2);
694  }
695  else
696  StatCovarianceMatrix[i-1][j-1] = 0;
697  }
698  }
699  }
700 
702  {
703  return CovarianceMatrix[binx][biny];
704  }
706  int biny)
707  {
708  return StatCovarianceMatrix[binx][biny];
709  }
710 
711  void NueCCIncTemplateFitter::fcn3Var(Int_t &npar, Double_t *gin,
712  Double_t &f, Double_t *par,
713  Int_t iflag)
714  {
715  f = Fitter_obj->myFunction3Vars(par[0],par[1],par[2]);
716  }
717 
718 
720  double nue_par,
721  double nc_par)
722  {
723  float data = 0;
724  float signal = 0;
725  float numu = 0;
726  float nc = 0;
727  float other = 0;
728  // Covariance Matrix Inversion First
729  TMatrixF Vij_stat = TMatrix(NumberTemplateBins,NumberTemplateBins);
730  TMatrixF Vij = TMatrix(NumberTemplateBins,NumberTemplateBins);
731  for(int i =0; i < NumberTemplateBins; i++){
732  for(int j = 0; j < NumberTemplateBins; j++){
733  Vij_stat[i][j] =
735  Vij[i][j] =
737  }
738  }
739  //Inversion
740  TMatrixD H3 = Vij_stat + Vij;
741  TDecompSVD svd(H3);
742  TMatrixD VijInv = svd.Invert();
743 
744  TMatrixD FitVector = TMatrixD(1,NumberTemplateBins);
745  TMatrixD FitVector2 = TMatrixD(NumberTemplateBins,1);
746 
747  //Calculate Chi Sq
748  for(int i =0; i < NumberTemplateBins; i++){
749  data = Fitter_obj->DataValues[i];
750  signal = Fitter_obj->SignalLike_Values[i];
751  numu = Fitter_obj->NumuLike_Values[i];
752  nc = Fitter_obj->NC_Values[i];
753  other = Fitter_obj->Other_Values[i];
754 
755  if( (data < 30) || (signal + numu + nc + other < 30)){
756  FitVector[0][i] = 0;
757  FitVector2[i][0] = 0;
758  continue;
759  }
760 
761  FitVector[0][i] = ( data - (nue_par*signal + numu_par*numu +
762  nc_par*nc + other));
763 
764  FitVector2[i][0] = ( data -(nue_par*signal + numu_par*numu +
765  nc_par*nc + other));
766  }
767 
768  TMatrixD Chi = FitVector * VijInv;
769  TMatrixD Chi2 = Chi * FitVector2;
770  float chisquare = Chi2(0,0);
771  return chisquare;
772  }
773 
774  void NueCCIncTemplateFitter::fcn2Var(Int_t &npar, Double_t *gin,
775  Double_t &f, Double_t *par,
776  Int_t iflag)
777  {
778  f = Fitter_obj->myFunction2Vars(par[0],par[1]);
779  }
780 
782  double nue_par)
783  {
784  float data = 0;
785  float signal = 0;
786  float numu = 0;
787  float nc = 0;
788  float other = 0;
789  // Covariance Matrix Inversion First
790  TMatrixF Vij_stat = TMatrix(numPIDBins,numPIDBins);
791  TMatrixF Vij = TMatrix(numPIDBins,numPIDBins);
792  for(int i =0; i < numPIDBins; i++){
793  for(int j = 0; j < NumberTemplateBins; j++){
794  Vij_stat[i][j] =
796  Vij[i][j] =
798  }
799  }
800  //Inversion
801  TMatrixD H3 = Vij_stat + Vij;
802  TDecompSVD svd(H3);
803  TMatrixD VijInv = svd.Invert();
804 
805  TMatrixD FitVector = TMatrixD(1,NumberTemplateBins);
806  TMatrixD FitVector2 = TMatrixD(NumberTemplateBins,1);
807 
808  //Calculate Chi Sq
809  for(int i =0; i < NumberTemplateBins; i++){
810  data = DataValues[i];
811  signal = SignalLike_Values[i];
812  numu = NumuLike_Values[i];
813  nc = NC_Values[i];
814  other = Other_Values[i];
815 
816  if( (data < 30) || (signal + numu + nc + other < 30)){
817  FitVector[0][i] = 0;
818  FitVector2[i][0] = 0;
819  continue;
820  }
821 
822  FitVector[0][i] = ( data - (nue_par*signal + bkgd_par*(numu + nc)
823  + other));
824 
825  FitVector2[i][0] = ( data -(nue_par*signal + bkgd_par*(numu + nc)
826  + other));
827  }
828 
829  TMatrixD Chi = FitVector * VijInv;
830  TMatrixD Chi2 = Chi * FitVector2;
831  float chisquare = Chi2(0,0);
832  return chisquare;
833  }
834 
836  {
837  //Weights are good here
838  double signal_weis = nue_wei.first;
839  double nummu_weis = numu_wei.first;
840  double nc_weis = nc_wei.first;
841 
842  double signal_err = nue_wei.second;
843  double nummu_err = numu_wei.second;
844  double nc_err = nc_wei.second;
845 
846  hSignalWeights2D->SetBinContent(binx,biny,signal_weis);
847  hNumuCCWeights2D->SetBinContent(binx,biny,nummu_weis);
848  hNCWeights2D->SetBinContent(binx,biny,nc_weis);
849 
850  hSignalError2D->SetBinContent(binx,biny,signal_err);
851  hNumuCCError2D->SetBinContent(binx,biny,nummu_err);
852  hNCError2D->SetBinContent(binx,biny,nc_err);
853  hChiSquaredResults2D->SetBinContent(binx,biny,chisquared);
854  }
855 
857  {
858  for(int numChn = 0; numChn < (int)hAnalysis3D.size(); numChn++){
859  for(int zbin = 1;
860  zbin <= hAnalysis3D[0]->GetZaxis()->GetNbins();zbin++){
861  float nsig = hAnalysis3D[1]->GetBinContent(binx,biny,zbin);
862  float nnuebar = hAnalysis3D[2]->GetBinContent(binx,biny,zbin);
863  float nnumu = hAnalysis3D[3]->GetBinContent(binx,biny,zbin);
864  float nnumubar = hAnalysis3D[4]->GetBinContent(binx,biny,zbin);
865  float nnc = hAnalysis3D[5]->GetBinContent(binx,biny,zbin);
866  float nother = hAnalysis3D[6]->GetBinContent(binx,biny,zbin);
867  float nnf = hAnalysis3D[7]->GetBinContent(binx,biny,zbin);
868  if(numChn == 0){
869  //Add all Channels to total
870  float binCont = nsig * nue_wei.first;
871  binCont += nnuebar * nue_wei.first;
872  binCont += nnumu * numu_wei.first;
873  binCont += nnumubar * numu_wei.first;
874  binCont += nnc * nc_wei.first;
875  binCont += nother;
876  binCont += nnf * nue_wei.first;
877 
878  std::cout << hAnalysis3D[0]->GetBinContent(binx,biny,zbin)
879  << "\t" << binCont << std::endl;
880 
881  hAnalysis3D[0]->SetBinContent(binx,biny,zbin,binCont);
882 
883 
884 
885  //MINOS Error calculations already take the correlations between
886  //parameters into account
887  float binErr2 =
888  pow(fit_covariance_matrix[0][0],2)*pow((nnumu+nnumubar),2) +
889  pow(fit_covariance_matrix[1][1],2)* pow((nsig+nnuebar+nnf),2) +
890  pow(fit_covariance_matrix[2][2],2)* pow((nnc),2) +
891  nsig * pow(nue_wei.first,2) + nnuebar * pow(nue_wei.first,2) +
892  nnf * pow(nue_wei.first,2) + nnumu * pow(numu_wei.first,2) +
893  nnumubar * pow(numu_wei.first,2) + nnc * pow(nc_wei.first,2) +
894  nother;
895 
896  if(printVerbose)
897  std::cout << binx << "\t" << biny << "\t" << zbin << "\t"
898  << sqrt(binErr2) <<std::endl;
899 
900  if(fit_covariance_matrix[0][1] == 0 &&
901  fit_covariance_matrix[0][2] == 0 &&
902  fit_covariance_matrix[1][2] == 0 && nue_wei.second == 0 &&
903  numu_wei.second == 0 && nc_wei.second == 0){
904  binErr2 = nsig + nnuebar + nnf + nnumu + nnumubar +
905  nnc + nother;
906  }
907  hAnalysis3D[0]->SetBinError(binx,biny,zbin,sqrt(binErr2));
908  }
909  if(numChn == 1){
910  //Signal Estimate
911  float binCont = nsig * nue_wei.first;
912 
913  //Just uncertainty on Signal estimate
914  //Signal + Background (taken into account in TMinuit internally
915  //when calculating the internal error matrix)
916 
917  //Need to add wrong sign estimate uncertainty
918  //Lets just add this in quadrature to the
919  //signal estimate + background estimate uncertainties
920  //We add uncertainties on the signal template estimate and
921  //statistical uncertainties on the WS component
922  float binErr = sqrt( pow(nue_wei.second,2)*pow(nsig,2) +
923  nsig * pow(nue_wei.first,2) +
924  pow(nue_wei.second,2)*pow(nnuebar,2) +
925  nnuebar * pow(nue_wei.first,2));
926 
927  if(printVerbose)
928  std::cout << binx << "\t" << biny << "\t" << zbin << "\t"
929  << nsig << "\t" << "\t" << sqrt(nsig) << "\t"
930  << binErr << std::endl;
931 
932  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
933  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
934  }
935  if(numChn == 2){
936  float binCont = nnuebar * nue_wei.first;
937  float binErr = sqrt( pow(nue_wei.second,2)*pow(nnuebar,2) +
938  nnuebar * pow(nue_wei.first,2));
939  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
940  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
941  }
942  if(numChn == 3){
943  float binCont = nnumu * numu_wei.first;
944  float binErr = sqrt( pow(numu_wei.second,2)*pow(nnumu,2) +
945  nnumu * pow(numu_wei.first,2));
946  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
947  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
948  }
949  if(numChn == 4){
950  float binCont = nnumubar * numu_wei.first;
951  float binErr = sqrt( pow(numu_wei.second,2)*pow(nnumubar,2) +
952  nnumubar * pow(numu_wei.first,2));
953  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
954  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
955  }
956  if(numChn == 5){
957  float binCont = nnc * nc_wei.first;
958  float binErr = sqrt( pow(nc_wei.second,2)*pow(nnc,2) +
959  nnc * pow(nc_wei.first,2));
960  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
961  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
962  }
963  if(numChn == 7){
964  float binCont = nnf * nue_wei.first;
965  float binErr = sqrt( pow(nue_wei.second,2)*pow(nnf,2) +
966  nnf * pow(nue_wei.first,2));
967  hAnalysis3D[numChn]->SetBinContent(binx,biny,zbin,binCont);
968  hAnalysis3D[numChn]->SetBinError(binx,biny,zbin,binErr);
969  }
970  }//loop over all Z axis bins (not electron kinematic bin)
971  }//loop over all interaction channels
972  }
973 
974  std::vector<std::pair<double,double>>
976  {
977  std::vector<std::pair<double,double>> vResult;
978 
979  vResult.push_back(std::make_pair(hSignalWeights2D->
980  GetBinContent(binx,biny),
982  GetBinContent(binx,biny)));
983  vResult.push_back(std::make_pair(hNumuCCWeights2D->
984  GetBinContent(binx,biny),
986  GetBinContent(binx,biny)));
987  vResult.push_back(std::make_pair(hNCWeights2D->GetBinContent(binx,biny),
988  hNCError2D->GetBinContent(binx,biny)));
989 
990  return vResult;
991  }
992 
993  std::vector<TH1F*> NueCCIncTemplateFitter::GetTemplates(int binx, int biny,
995  {
996  std::vector<TH1F*> hResults;
997  for(uint i = 0; i < hTemplates3D.size(); i++){
998  hResults.push_back((TH1F*)hTemplates3D[i]->ProjectionZ
999  (ana::UniqueName().c_str(),binx,binx,biny,biny));
1000  }
1001  return hResults;
1002  }
1003 
1005  {
1006  //Make Sure this bin isn't excluded by the phase space cuts
1007  //If it is, don't worrky about calculating anything
1008  //Just checking one of the inner bins
1009  TH2F* hResultBad = new TH2F(ana::UniqueName().c_str(),";",1,0,1,1,0,1);
1010  if(hCV3D->Integral(binx,binx,biny,biny,1,hCV3D->GetZaxis()->GetNbins()) ==
1011  0) return hResultBad;
1012  //Systematic Samples
1013  //Turn Systamtic samples that are not generated from a multiverse
1014  //into a "one sigma" shift that we can randomly draw from in a moment
1015  std::vector<TH1F*> hSystsShifts =
1017 
1018  std::vector<std::pair<std::string,std::vector<double>>> hSysts_Pairs;
1019 
1020  std::vector<std::string> names = {"Cherenkov","Cal. Shape",
1021  "Light","Cal. Norm.",
1022  "Horn Current", "Beam Size",
1023  "Beam Shift X", "Beam Shift Y",
1024  "Horn 1 Shift X", "Horn 1 Shift Y",
1025  "Horn 2 Shift X", "Horn 2 Shift Y",
1026  "Target Z Shift", "Magnetic Field",
1027  "Horn Wather Thickness",
1028  "PPFX","GENIE"};
1029 
1030  TRandom *r0 = new TRandom();
1031 
1032  //Randomly draw from these one sigma shift samples
1033  //Draw 1000 from each sample
1034  for(int j = 0; j < (int)hSystsShifts.size(); j++){
1035  std::vector<double> vPseudoExp;
1036  for(int i = 0; i < 1000; i++){
1037  TH1F* hClone = (TH1F*)hCV3D->ProjectionZ(ana::UniqueName().c_str(),
1038  binx,binx,biny,biny);
1039  double wei = r0->Gaus(0,1);
1040  hClone->Add(hSystsShifts[j],wei);
1041  vPseudoExp.push_back(hClone->Integral(1,-1));
1042  delete hClone;
1043  }//loop over systematic samples
1044  hSysts_Pairs.push_back(std::make_pair(names[j],vPseudoExp));
1045  }//loop over pseudo experiments
1046 
1047  std::vector<double> vPPFX_Values;
1048  std::vector<double> vGenie_Values;
1049 
1050  for(int systNum = 0; systNum < (int)hSystsMultiverse.size();systNum++){
1051  double integral =
1052  hSystsMultiverse[systNum]->Integral(binx,binx,biny,biny,1,-1);
1053  if(systNum < 300)vGenie_Values.push_back(integral);
1054  else vPPFX_Values.push_back(integral);
1055  }//loop over multiverse universes
1056 
1057  hSysts_Pairs.push_back(std::make_pair("ppfx", vPPFX_Values));
1058  hSysts_Pairs.push_back(std::make_pair("genie", vGenie_Values));
1059 
1060  //Calculate Covariance
1061  //Multiverses are done first
1062  //Followed by the "multiverse" we created from the other shifts
1063  //above
1064  const int systNum = (int)hSysts_Pairs.size();
1065  double syst_covariance[systNum][systNum];
1066  for(int i = 0; i < systNum; i++){
1067  for(int j = 0; j < systNum; j++){
1068  double numerator_term = 0;
1069  double counter = hSysts_Pairs[i].second.size();
1070  for(int nUniv = 0; nUniv < counter; nUniv++){
1071  double term_i = hSysts_Pairs[i].second[nUniv] -
1072  hCV3D->Integral(binx,binx,biny,biny,1,-1);
1073  double term_j = hSysts_Pairs[j].second[nUniv] -
1074  hCV3D->Integral(binx,binx,biny,biny,1,-1);
1075  numerator_term += term_i*term_j * (1.0/(counter - 1));
1076  }
1077  double covariance = numerator_term;
1078  syst_covariance[i][j] = covariance;
1079  }
1080  }
1081 
1082 
1083  TH2F* hResult = new TH2F(ana::UniqueName().c_str(),";",
1084  systNum,0,systNum,systNum,0,systNum);
1085 
1086  //Move the covariance matrix to the final result
1087  for(int i = 0; i < systNum; i++){
1088  for(int j = 0; j < systNum; j++){
1089  hResult->SetBinContent(i+1,j+1,syst_covariance[i][j]);
1090  hResult->GetXaxis()->SetBinLabel(i+1,names[i].c_str());
1091  hResult->GetYaxis()->SetBinLabel(j+1,names[j].c_str());
1092  }
1093  }
1094 
1095  hResult->GetZaxis()->SetTitle("Covariance (Events^{2})");
1096 
1097  return hResult;
1098  }
1099 
1100  //////////// //////////// //////////// //////////// ////////////
1101  //Public Functions//////////////////////////////////////////////////
1102  //////////// //////////// //////////// //////////// ////////////
1104  {
1105  Fitter_obj = this;
1106  }
1107 
1109  {
1110  std::cout << "Values For Bins: " << binx << "\t" << biny << std::endl;
1112 
1113  for(unsigned int i = 0; i < numPIDBins; i++){
1114  std::cout << "Data: " << DataValues[i] << "\tSignal: " <<
1115  SignalLike_Values[i] << "\tNumu: " << NumuLike_Values[i] <<
1116  "\tNC: " << NC_Values[i] << "\tOther: " <<
1117  Other_Values[i] << std::endl;
1118  }
1119  }
1120 
1121 
1123  {
1124  std::vector<TH3F*> badentry;
1125  bool CheckSysts = NueCCIncTemplateFitter::CheckSystsSize();
1126  if(!CheckSysts){
1127  std::cerr << "Unequal up and down shifts. Exiting" << std::endl;
1128  return badentry;
1129  }
1130 
1131  std::vector<TH3F*> hResults =
1133 
1134  //Loop over all analysis bins (x and y axis bins)
1135  for(int xbin = 1; xbin <= hResults[0]->GetXaxis()->GetNbins(); xbin++){
1136  for(int ybin = 1; ybin <= hResults[0]->GetYaxis()->GetNbins(); ybin++){
1137 
1138  std::string caption =
1140 
1142 
1143  bool b_fit_this_bin =
1145 
1146  bool check_high_lo_energy = true;
1147  if(hResults[0]->GetYaxis()->GetBinLowEdge(ybin) == 6.0 ||
1148  hResults[0]->GetXaxis()->GetBinLowEdge(xbin) == 0.75)
1149  check_high_lo_energy = false;
1150 
1151  if(b_fit_this_bin && check_high_lo_energy){
1152  std::cout << "Doing Fit in this bin" << std::endl;
1156 
1157  //Now Perform the Fit
1158  //Fit 1: Get Good Starting Fit Values, Quickly
1159  //Three Parameter fit
1160 
1161  //Prepare Starting Seeds
1162  std::vector<float> nue_starting_points;
1163  std::vector<float> numu_starting_points;
1164  std::vector<float> nc_starting_points;
1165 
1166  for(float fillVal = 0.6; fillVal < 1.8; fillVal+=0.15){
1167  nue_starting_points.push_back(fillVal);
1168  numu_starting_points.push_back(fillVal);
1169  nc_starting_points.push_back(fillVal);
1170  }
1171 
1172 
1173  std::random_shuffle(nue_starting_points.begin(),
1174  nue_starting_points.end());
1175  std::random_shuffle(numu_starting_points.begin(),
1176  numu_starting_points.end());
1177  std::random_shuffle(nc_starting_points.begin(),
1178  nc_starting_points.end());
1179 
1180  //Try 3 Parameter Fit First
1181  std::vector<std::pair<float,std::vector<float>>>
1182  pairFitResultStartPts3Var;
1183 
1184  std::vector<std::vector<double>> vUncertainties;
1185  std::vector<std::vector<double>> vNormalizations;
1186 
1187  std::cout << nue_starting_points.size() <<
1188  "\t" << nc_starting_points.size() << "\t"
1189  << numu_starting_points.size() << std::endl;
1190 
1191  int counter = 0;
1192 
1193  for(auto nuestart : nue_starting_points){
1194  for(auto numustart : numu_starting_points){
1195  for(auto ncstart : nc_starting_points){
1196  counter++;
1197  TMinuit *gMinuit = new TMinuit(3);
1198  gMinuit->SetPrintLevel(-1);
1199  gMinuit->SetFCN(fcn3Var);
1200 
1201  Double_t arglist[4];
1202  Int_t ierflg = 0;
1203  arglist[0] = 1;
1204  gMinuit->mnexcm("SET ERR", arglist,1,ierflg);
1205 
1206  //Define Fit Parameters
1207  //0.15 is the approximate error on the parameter
1208  //0,0 refers to the lower and upper bounds on the parameters
1209  gMinuit->mnparm(0, "Numu-Like Scaling",
1210  numustart, 0.15, 0, 0, ierflg);
1211  gMinuit->mnparm(1, "Nue-Like Scaling",
1212  nuestart, 0.15, 0, 0, ierflg);
1213  gMinuit->mnparm(2, "NC Scaling",
1214  ncstart, 0.15, 0, 0, ierflg);
1215 
1216  arglist[0] = 0; //number of iterations 0 == no limit
1217  arglist[1] = 1; //1 == standard minimization strategy
1218  //SIMPLEX should provide a good starting point the next fits
1219  gMinuit->Command("SIMPLEX");
1220 
1221  //Fit 2: Improved Fit + Better Eror Estimation
1222  arglist[1] = 1; //2 == try to improve minimum (slower)
1223  gMinuit->mnexcm("SET STR", arglist,1,ierflg);
1224 
1225  gMinuit->Command("SIMPLEX");
1226  gMinuit->Command("HESSE");
1227  gMinuit->Command("MIGRAD");
1228  gMinuit->Command("MINOS");
1229 
1230 
1231  //Grab Final Results
1232  double vstart[3];
1233  double verror[3];
1234  double fParNewStart;
1235  double fParNewErr;
1236  gMinuit->GetParameter(0,fParNewStart,fParNewErr);
1237  vstart[0] = fParNewStart;
1238  verror[0] = fParNewErr;
1239  gMinuit->GetParameter(1,fParNewStart,fParNewErr);
1240  vstart[1] = fParNewStart;
1241  verror[1] = fParNewErr;
1242  gMinuit->GetParameter(2,fParNewStart,fParNewErr);
1243  vstart[2] = fParNewStart;
1244  verror[2] = fParNewErr;
1245 
1246  //Before Ending Fits,
1247  //Make Sure There Was No Issue With the Fit
1248  bool troubleFitting = false;
1249  Double_t fmin = 0,fedm = 0,ferrdef = 0;
1250  Int_t npari =0, nparx = 0, istat = 0;
1251  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1252 
1253  Double_t corr_mat_test[3][3];
1254  gMinuit->mnemat(*corr_mat_test,3);
1255 
1256  //Turn this test error matix into a correlation matrix
1257  double corr_01 =
1258  corr_mat_test[0][1]/(pow(corr_mat_test[0][0],0.5) *
1259  pow(corr_mat_test[1][1],0.5));
1260  double corr_02=
1261  corr_mat_test[0][2]/((pow(corr_mat_test[0][0],0.5) *
1262  pow(corr_mat_test[2][2],0.5)));
1263  double corr_12=
1264  corr_mat_test[1][2]/((pow(corr_mat_test[1][1],0.5) *
1265  pow(corr_mat_test[2][2],0.5)));
1266 
1267  if(istat != 3 ||
1268  vstart[0] < 0 || vstart[1] < 0 || vstart[2] < 0 ||
1269  vstart[0] > 2 || vstart[1] > 2 || vstart[2] > 2 ||
1270  fabs(corr_01) > 0.97 || fabs(corr_02) > 0.97 ||
1271  fabs(corr_12) > 0.97)
1272  troubleFitting = true;
1273 
1274  std::vector<double> vUncert = {verror[0]/vstart[0],
1275  verror[1]/vstart[1],
1276  verror[2]/vstart[2]};
1277  std::vector<double> vNormal = {vstart[0],vstart[1],vstart[2]};
1278 
1279  std::vector<float> vStarts = {nuestart,numustart,ncstart};
1280 
1281  if(!troubleFitting){
1282  pairFitResultStartPts3Var.push_back(std::make_pair
1283  (fmin,vStarts));
1284  vUncertainties.push_back(vUncert);
1285  vNormalizations.push_back(vNormal);
1286  }
1287 
1288  } //loop over nc starting values
1289  } //loop over numu starting values
1290  }//loop over nue starting values
1291 
1292  if(printVerbose){
1293  for(uint i = 0; i < pairFitResultStartPts3Var.size(); i++){
1294  std::cout << pairFitResultStartPts3Var[i].first << "\t"
1295  << vNormalizations[i][0] << "\t"
1296  << vUncertainties[i][0] << "\t"
1297  << vNormalizations[i][1] << "\t"
1298  << vUncertainties[i][1] << "\t"
1299  << vNormalizations[i][0] << "\t"
1300  << vUncertainties[i][2] << "\t" << std::endl;
1301  }
1302  }
1303 
1304  bool shouldSwitchTo2VarFit = false;
1305  //Find the best fits
1306  int saveIndex = 0;
1307  if(pairFitResultStartPts3Var.size() > 0){
1308  //now we can loop through all of the items in the vector
1309  for(uint i = 0; i < pairFitResultStartPts3Var.size(); i++)
1310  {
1311  if(floor(pairFitResultStartPts3Var[saveIndex].first*1e3) >
1312  floor(pairFitResultStartPts3Var[i].first*1e3))
1313  {
1314  saveIndex = i;
1315  chisquared = pairFitResultStartPts3Var[saveIndex].first;
1316  }
1317  if(floor(pairFitResultStartPts3Var[saveIndex].first*1e3) ==
1318  floor(pairFitResultStartPts3Var[i].first*1e3))
1319  {
1320  if(vUncertainties[i][0] <
1321  vUncertainties[saveIndex][0] &&
1322  vUncertainties[i][1] <
1323  vUncertainties[saveIndex][1] &&
1324  vUncertainties[i][2] <
1325  vUncertainties[saveIndex][2])
1326  {
1327  saveIndex = i;
1328  chisquared = pairFitResultStartPts3Var[saveIndex].first;
1329  }
1330  }
1331  }
1332  }
1333  else shouldSwitchTo2VarFit = true;
1334 
1335 
1336  //Put 3 Var Fit Here
1337  if(!shouldSwitchTo2VarFit){
1338  //Fit 1: Get Good Starting Fit Values, Quickly
1339  //Three Parameter fit
1340  TMinuit *gMinuit = new TMinuit(3);
1341  gMinuit->SetPrintLevel(-1);
1342  if(printVerbose)gMinuit->SetPrintLevel(1);
1343  gMinuit->SetFCN(fcn3Var);
1344  Double_t arglist[4];
1345  Int_t ierflg = 0;
1346  arglist[0] = 1;
1347  gMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1348 
1349 
1350  double vstart[3] = {0.8,0.9,1.15};
1351  double verror[3] = {0,0,0};
1352  double vstep[3] = {0.15,0.15,0.15};
1353 
1354  //Define Fit Parameters
1355  gMinuit->mnparm(0, "Numu-Like Scaling",
1356  pairFitResultStartPts3Var[saveIndex].second[1],
1357  vstep[0], 0, 0, ierflg);
1358  gMinuit->mnparm(1, "Nue-Like Scaling",
1359  pairFitResultStartPts3Var[saveIndex].second[0],
1360  vstep[1], 0, 0, ierflg);
1361  gMinuit->mnparm(2, "NC Scaling",
1362  pairFitResultStartPts3Var[saveIndex].second[2],
1363  vstep[2], 0, 0, ierflg);
1364 
1365  arglist[0] = 0; //number of iterations 0 == no limit
1366  arglist[1] = 1; //1 == standard minimization strategy
1367 
1368  //SIMPLEX should provide a good starting point for further fits
1369  gMinuit->Command("SIMPLEX");
1370 
1371  //Fit 2: Improved Fit + Better Eror Estimation
1372  arglist[1] = 2; //2 == try to improve minimum (slower)
1373  gMinuit->mnexcm("SET STR", arglist,1,ierflg);
1374  gMinuit->Command("SIMPLEX");
1375 
1376  Double_t fmin = 0,fedm = 0,ferrdef = 0;
1377  Int_t npari =0, nparx = 0, istat = 0;
1378  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1379 
1380  gMinuit->Command("SIMPLEX");
1381  gMinuit->Command("HESSE");
1382  gMinuit->Command("MIGRAD");
1383  gMinuit->Command("MINOS");
1384 
1385  //Grab Final Results
1386  double fParNewStart;
1387  double fParNewErr;
1388  gMinuit->GetParameter(0,fParNewStart,fParNewErr);
1389  vstart[0] = fParNewStart;
1390  verror[0] = fParNewErr;
1391  gMinuit->GetParameter(1,fParNewStart,fParNewErr);
1392  vstart[1] = fParNewStart;
1393  verror[1] = fParNewErr;
1394  gMinuit->GetParameter(2,fParNewStart,fParNewErr);
1395  vstart[2] = fParNewStart;
1396  verror[2] = fParNewErr;
1397 
1398  //Before Ending Fits, Make Sure There Was No Issue With the Fit
1399  gMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1400 
1401  Double_t corr_mat_test[3][3];
1402  gMinuit->mnemat(*corr_mat_test,3);
1403 
1404  //Turn this test error matix into a correlation matrix
1405  double corr_01 =
1406  corr_mat_test[0][1]/(pow(corr_mat_test[0][0],0.5) *
1407  pow(corr_mat_test[1][1],0.5));
1408  double corr_02=
1409  corr_mat_test[0][2]/((pow(corr_mat_test[0][0],0.5) *
1410  pow(corr_mat_test[2][2],0.5)));
1411  double corr_12=
1412  corr_mat_test[1][2]/((pow(corr_mat_test[1][1],0.5) *
1413  pow(corr_mat_test[2][2],0.5)));
1414 
1415 
1416  if(istat != 3 ||
1417  vstart[0] < 0 || vstart[1] < 0 || vstart[2] < 0 ||
1418  vstart[0] > 2. || vstart[1] > 2 || vstart[2] > 2.0 ||
1419  fabs(corr_01) > 0.98 || fabs(corr_02) > 0.98 ||
1420  fabs(corr_12) > 0.98)
1421  shouldSwitchTo2VarFit = true;
1422 
1423  //Get the covariance matrix
1424  Double_t emat[3][3];
1425  gMinuit->mnemat(*emat,3);
1426 
1427  nue_wei = std::make_pair (vstart[1],verror[1]);
1428  numu_wei = std::make_pair (vstart[0],verror[0]);
1429  nc_wei = std::make_pair (vstart[2],verror[2]);
1430  fit_covariance_matrix[0][0] = emat[0][0];
1431  fit_covariance_matrix[0][1] = emat[0][1];
1432  fit_covariance_matrix[0][2] = emat[0][2];
1433  fit_covariance_matrix[1][0] = emat[1][0];
1434  fit_covariance_matrix[1][1] = emat[1][1];
1435  fit_covariance_matrix[1][2] = emat[1][2];
1436  fit_covariance_matrix[2][0] = emat[2][0];
1437  fit_covariance_matrix[2][1] = emat[2][1];
1438  fit_covariance_matrix[2][2] = emat[2][2];
1439  std::cout << shouldSwitchTo2VarFit <<
1440  "\tMade it Here" << std::endl;
1441  }
1442 
1443  //Now If there were issues with the 3 Parameter Fit,
1444  //switch to a 2 ParameterFit
1445  std::vector<std::pair<float,std::vector<float>>>
1446  pairFitResultStartPts2Var;
1447 
1448  std::vector<std::vector<double>> vUncertainties2Var;
1449  std::vector<std::vector<double>> vNormalizations2Var;
1450 
1451  std::vector<float> signal_starting_points;
1452  std::vector<float> background_starting_points;
1453 
1454  for(float fillVal = 0.6; fillVal < 1.8; fillVal+=0.15){
1455  signal_starting_points.push_back(fillVal);
1456  background_starting_points.push_back(fillVal);
1457  }
1458 
1459  if(shouldSwitchTo2VarFit){
1460  //There Was trouble with the original fits,
1461  //most likely this was due to extremely
1462  //high correlation (approx 1)
1463  //between two of the fitting parameters
1464  //Switch to a fit that adjusts all backgrounds
1465  for(auto nuestart : signal_starting_points){
1466  for(auto bkgstart : background_starting_points){
1467  bool isGoodResult = true;
1468  Double_t fmin = 0,fedm = 0,ferrdef = 0;
1469  Int_t npari =0, nparx = 0, istat = 0;
1470  Double_t arglist[4];
1471  Int_t ierflg = 0;
1472  TMinuit *jMinuit = new TMinuit(2);
1473  jMinuit->SetPrintLevel(-1);
1474  jMinuit->SetFCN(fcn2Var);
1475  arglist[0] = 1;
1476  jMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1477 
1478  jMinuit->mnparm(0, "Background Scaling",bkgstart,
1479  0.15, 0,0, ierflg);
1480  jMinuit->mnparm(1, "Nue CC Scaling",
1481  nuestart, 0.15, 0,0, ierflg);
1482 
1483  arglist[0] = 0; //number of iterations 0 == no limit
1484  arglist[1] = 1;
1485  jMinuit->Command("SIMPLEX");
1486  jMinuit->Command("HESSE");
1487 
1488  //Fit 2: Improved Fit + Better Eror Estimation
1489  arglist[1] = 2; //2 == try to improve minimum (slower)
1490  jMinuit->mnexcm("SET STR", arglist,1,ierflg);
1491 
1492  jMinuit->Command("SIMPLEX");
1493  jMinuit->Command("HESSE");
1494  jMinuit->Command("MIGRAD");
1495  //jMinuit->Command("MINOS");
1496 
1497 
1498  //Before Ending Fits,Make Sure There Was No Issue With the Fit
1499  jMinuit->mnstat(fmin,fedm,ferrdef,npari,nparx,istat);
1500 
1501 
1502  //errors are
1503  // ISTAT:a status integerindicating how good is the covariance
1504  //*-* matrix:0= not calculated at all
1505  //*-* 1= approximation only, not accurate
1506  //*-* 2= full matrix, but forced positive-def
1507  //*-* 3= full accurate covariance matrix
1508  if(istat != 3)
1509  isGoodResult = false;
1510 
1511 
1512  //Check for high correlations between fit parameters,
1513  //can show an issue with the fit in this bin
1514  Double_t corr_mat_test[2][2];
1515  jMinuit->mnemat(*corr_mat_test,2);
1516 
1517 
1518  //Turn this test error matix into a correlation matrix
1519  double corr_01 =
1520  corr_mat_test[0][1]/(pow(corr_mat_test[0][0],0.5) *
1521  pow(corr_mat_test[1][1],0.5));
1522 
1523  if(fabs(corr_01) > 0.97)
1524  isGoodResult = false;
1525 
1526 
1527 
1528  std::vector<double> vUncert;
1529  std::vector<double> vNormal;
1530  for(int i = 0; i < 2; i++){
1531  double fParNormalization,fParError;
1532  jMinuit->GetParameter(i,fParNormalization,fParError);
1533  vUncert.push_back(fParError);
1534  vNormal.push_back(fParNormalization);
1535  if(fParNormalization < 0 || fParNormalization > 2 ||
1536  fParError/fParNormalization > 1)
1537  isGoodResult = false;
1538  if(!isGoodResult)break;
1539  }
1540 
1541  std::vector<float> vStart = {nuestart,bkgstart};
1542 
1543  if(isGoodResult)
1544  {pairFitResultStartPts2Var.push_back(std::make_pair
1545  (fmin,vStart));
1546  vUncertainties2Var.push_back(vUncert);
1547  vNormalizations2Var.push_back(vNormal);
1548 
1549  }
1550 
1551 
1552  }//background starting values loop
1553  } //nue starting values loop
1554  }//do 2Var Fit
1555 
1556  if(printVerbose){
1557  for(uint i = 0; i < pairFitResultStartPts2Var.size(); i++){
1558  std::cout << pairFitResultStartPts2Var[i].first << "\t"
1559  << vNormalizations2Var[i][0] << "\t"
1560  << vUncertainties2Var[i][0] << "\t"
1561  << vNormalizations2Var[i][1] << "\t"
1562  << vUncertainties2Var[i][1] << std::endl;
1563  }
1564  }
1565 
1566  //Grab best fit by ChiSquared Value
1567  //If those are equal sort by predicted uncertainties
1568  //All Uncertainties need to be better to be considered the best fit
1569 
1570  TH1F* hUncertaintyAverage = new TH1F("hUncertaintyAverage",
1571  ";Signal Fit Uncertainty;Count",
1572  100,0,2.0);
1573 
1574  saveIndex = 0;
1575 
1576 
1577  if(shouldSwitchTo2VarFit){
1578  for(uint i = 0; i < pairFitResultStartPts2Var.size(); i++)
1579  {
1580  hUncertaintyAverage->Fill(vUncertainties2Var[i][1]);
1581  if(floor(pairFitResultStartPts2Var[saveIndex].first*1e5) >
1582  floor(pairFitResultStartPts2Var[i].first*1e5))
1583  {
1584  saveIndex = i;
1585  chisquared = pairFitResultStartPts2Var[saveIndex].first;
1586  shouldSwitchTo2VarFit = true;
1587  }
1588  }
1589  }
1590 
1591  std::cout << "SaveIndex " << saveIndex
1592  << "************" << std::endl;
1593  std::cout << pairFitResultStartPts2Var.size() << std::endl;
1594 
1595 
1596  if(shouldSwitchTo2VarFit){
1597  //There Was trouble with the original fits,
1598  //most likely this was due to extremely high correlation (approx 1)
1599  //between two of the fitting parameters
1600  //Switch to a fit that adjusts all backgrounds
1601  double vstart[3] = {0.8,0.9,1.15};
1602  double verror[3] = {0,0,0};
1603  double vstep[3] = {0.15,0.15,0.15};//{0.01,0.01,0.01
1604 
1605  Double_t arglist[4];
1606  Int_t ierflg = 0;
1607  TMinuit *jMinuit = new TMinuit(2);
1608  jMinuit->SetPrintLevel(-1);
1609  if(printVerbose)jMinuit->SetPrintLevel(1);
1610  jMinuit->SetFCN(fcn2Var);
1611  arglist[0] = 1;
1612  jMinuit->mnexcm("SET ERR", arglist, 1, ierflg);
1613 
1614  double nue_start = 1.0;
1615  double bkgd_start = 1.0;
1616 
1617  if(pairFitResultStartPts2Var.size() > 0){
1618  nue_start = pairFitResultStartPts2Var[saveIndex].second[1];
1619  bkgd_start = pairFitResultStartPts2Var[saveIndex].second[0];
1620  std::cout << "Made it here" << std::endl;
1621  std::cout << nue_start << "\t" << bkgd_start << std::endl;
1622  }
1623 
1624 
1625  jMinuit->mnparm(0, "Background Scaling",
1626  bkgd_start,
1627  vstep[0], 0, 0, ierflg);
1628  jMinuit->mnparm(1, "Nue CC Scaling",
1629  nue_start,
1630  vstep[1], 0, 0, ierflg);
1631 
1632  arglist[0] = 0; //number of iterations 0 == no limit
1633  arglist[1] = 1;
1634 
1635  jMinuit->Command("SIMPLEX");
1636  jMinuit->Command("HESSE");
1637 
1638  //Fit 2: Improved Fit + Better Eror Estimation
1639  arglist[1] = 2; //2 == try to improve minimum (slower)
1640  jMinuit->mnexcm("SET STR", arglist,1,ierflg);
1641 
1642  jMinuit->Command("SIMPLEX");
1643  for(int i = 0; i < 2; i++){
1644  double fParNewStart = 0;
1645  double fParNewErr = 0;
1646  jMinuit->GetParameter(i,fParNewStart,fParNewErr);
1647  std::cout << "SIMPLEX" << "\t" << fParNewStart << "\t"
1648  << fParNewErr << std::endl;
1649  }
1650  jMinuit->Command("MIGRAD");
1651  for(int i = 0; i < 2; i++){
1652  double fParNewStart = 0;
1653  double fParNewErr = 0;
1654  jMinuit->GetParameter(i,fParNewStart,fParNewErr);
1655  std::cout << "MIGRAD" << "\t" << fParNewStart << "\t"
1656  << fParNewErr << std::endl;
1657  }
1658  //jMinuit->mnexcm("SEEk", arglist, 3, ierflg);
1659  jMinuit->Command("HESSE");
1660  for(int i = 0; i < 2; i++){
1661  double fParNewStart = 0;
1662  double fParNewErr = 0;
1663  jMinuit->GetParameter(i,fParNewStart,fParNewErr);
1664  std::cout << "SEEK" << "\t" << fParNewStart << "\t"
1665  << fParNewErr << std::endl;
1666  }
1667  //jMinuit->Command("MINOS");
1668  for(int i = 0; i < 2; i++){
1669  double fParNewStart = 0;
1670  double fParNewErr = 0;
1671  jMinuit->GetParameter(i,fParNewStart,fParNewErr);
1672  std::cout << "MINOS" << "\t" << fParNewStart << "\t"
1673  << fParNewErr << std::endl;
1674  }
1675 
1676  double fParNewStart;
1677  double fParNewErr;
1678  jMinuit->GetParameter(0,fParNewStart,fParNewErr);
1679  vstart[0] = fParNewStart;
1680  verror[0] = fParNewErr;
1681  jMinuit->GetParameter(0,fParNewStart,fParNewErr);
1682  vstart[2] = fParNewStart;
1683  verror[2] = fParNewErr;
1684  jMinuit->GetParameter(1,fParNewStart,fParNewErr);
1685  vstart[1] = fParNewStart;
1686  verror[1] = fParNewErr;
1687  //Get the covariance matrix
1688  Double_t emat[2][2];
1689  jMinuit->mnemat(*emat,2);
1690 
1691  nue_wei = std::make_pair (vstart[1],verror[1]);
1692  numu_wei = std::make_pair (vstart[0],verror[0]);
1693  nc_wei = std::make_pair (vstart[2],verror[2]);
1694  fit_covariance_matrix[0][0] = emat[0][0];
1695  fit_covariance_matrix[0][1] = emat[0][1];
1696  fit_covariance_matrix[0][2] = 0;
1697  fit_covariance_matrix[1][0] = emat[1][0];
1698  fit_covariance_matrix[1][1] = emat[1][1];
1699  fit_covariance_matrix[1][2] = 0;
1700  fit_covariance_matrix[2][0] = emat[0][0];
1701  fit_covariance_matrix[2][1] = emat[0][1];
1702  fit_covariance_matrix[2][2] = emat[0][0];
1703  std::cout << shouldSwitchTo2VarFit <<
1704  "\tMade it Here, 2Var Fit true" << std::endl;
1705  }//2var results
1706 
1707 
1708  }//fit this bin
1709  else{
1710  std::cout << "Skipping Fit in " << caption << std::endl;
1711 
1712  nue_wei = std::make_pair (1,0);
1713  numu_wei = std::make_pair (1,0);
1714  nc_wei = std::make_pair (1,0);
1715  chisquared = 0;
1716  fit_covariance_matrix[0][0] = 0;
1717  fit_covariance_matrix[0][1] = 0;
1718  fit_covariance_matrix[0][2] = 0;
1719  fit_covariance_matrix[1][0] = 0;
1720  fit_covariance_matrix[1][1] = 0;
1721  fit_covariance_matrix[1][2] = 0;
1722  fit_covariance_matrix[2][0] = 0;
1723  fit_covariance_matrix[2][1] = 0;
1724  fit_covariance_matrix[2][2] = 0;
1725  }
1726 
1727  std::cout << "*********Fit for bin " << caption << " : " << std::endl;
1728  std::cout << "NumuCC Scaling Factor= " << numu_wei.first << " +- " <<
1729  numu_wei.second << std::endl;
1730  std::cout << "NueCC Scaling Factor= " << nue_wei.first << " +- " <<
1731  nue_wei.second << std::endl;
1732  std::cout << "NC Scaling Factor= " << nc_wei.first << " +- " <<
1733  nc_wei.second << std::endl;
1734 
1737  }//loop over ybins
1738  }//loop over xbins
1739 
1740  //Now That we have gone over every Bin and propagated uncertainty
1741  //Return the final results
1742  std::vector<TH3F*> hResultsPostFit =
1744 
1745  for(uint i = 0; i < hResults.size(); i++){
1746  std::cout << i << "\t" << hResults[i]->Integral() << "\t"
1747  << hResultsPostFit[i]->Integral() << std::endl;
1748  }
1749 
1750 
1751  return hResultsPostFit;
1752  }
1753 
1754 
1755  std::vector<TH1F*> NueCCIncTemplateFitter::getPIDFitTemplates(int binx,
1756  int biny)
1757  {
1758  std::vector<TH1F*> hResults =
1759  NueCCIncTemplateFitter::GetTemplates(binx, biny,"fit");
1760  std::string caption =
1762  std::vector<std::pair<double,double>> vWeights =
1764  for(int bin = 1; bin <= hResults[0]->GetXaxis()->GetNbins(); bin++){
1765  double nue = hResults[1]->GetBinContent(bin) *vWeights[0].first;
1766  double nuebar = hResults[2]->GetBinContent(bin) *vWeights[0].first;
1767  double numu = hResults[3]->GetBinContent(bin) *vWeights[1].first;
1768  double numubar= hResults[4]->GetBinContent(bin) *vWeights[1].first;
1769  double nc = hResults[5]->GetBinContent(bin) *vWeights[2].first;
1770  double other = hResults[6]->GetBinContent(bin);
1771  double nuenf = hResults[7]->GetBinContent(bin) *vWeights[0].first;
1772 
1773  hResults[0]->SetBinContent(bin,nue+nuebar+numu+numubar+nc+other+nuenf);
1774  hResults[1]->SetBinContent(bin,nue);
1775  hResults[2]->SetBinContent(bin,nuebar);
1776  hResults[3]->SetBinContent(bin,numu);
1777  hResults[4]->SetBinContent(bin,numubar);
1778  hResults[5]->SetBinContent(bin,nc);
1779  hResults[6]->SetBinContent(bin,other);
1780  hResults[7]->SetBinContent(bin,nuenf);
1781  }//loop over template bins
1782  hResults[0]->SetTitle(caption.c_str());
1783  return hResults;
1784  }
1785 
1786  std::vector<std::pair<TH2F*,TH2F*>>
1788  {
1789  std::vector<std::pair<TH2F*,TH2F*>> hResults;
1790 
1791  std::vector<TH3F*> vAnalysisBinning =
1793 
1794  TH2F* hTemplate = (TH2F*)vAnalysisBinning.front()->Project3D("yx");
1795  hTemplate->SetName("hFitWeightsErrorTemplate");
1796 
1797  TH2F* hSignalNorm = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1798  TH2F* hSignalError = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1799 
1800  TH2F* hNumuNorm = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1801  TH2F* hNumuError = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1802 
1803  TH2F* hNCNorm = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1804  TH2F* hNCError = (TH2F*)hTemplate->Clone(ana::UniqueName().c_str());
1805 
1806  for(int xbin = 1; xbin <= hTemplate->GetXaxis()->GetNbins(); xbin++){
1807  for(int ybin = 1; ybin <= hTemplate->GetYaxis()->GetNbins(); ybin++){
1808  std::vector<std::pair<double,double>> fit_results =
1810  hSignalNorm->SetBinContent(xbin,ybin,fit_results[0].first);
1811  hSignalError->SetBinContent(xbin,ybin,fit_results[0].second);
1812  hNumuNorm->SetBinContent(xbin,ybin,fit_results[1].first);
1813  hNumuError->SetBinContent(xbin,ybin,fit_results[1].second);
1814  hNCNorm->SetBinContent(xbin,ybin,fit_results[2].first);
1815  hNCError->SetBinContent(xbin,ybin,fit_results[2].second);
1816  }//loop over ybins
1817  }//loop over xbins
1818 
1819  hResults.push_back(std::make_pair(hSignalNorm,hSignalError));
1820  hResults.push_back(std::make_pair(hNumuNorm,hNumuError));
1821  hResults.push_back(std::make_pair(hNCNorm,hNCError));
1822 
1823  return hResults;
1824  }
1825 
1826 
1827  std::vector<TH1F*> NueCCIncTemplateFitter::getPIDTemplates(int binx,
1828  int biny)
1829  {
1830  std::vector<TH1F*> hResults =
1831  NueCCIncTemplateFitter::GetTemplates(binx, biny,"nominal");
1832  std::string caption =
1834  hResults[0]->SetTitle(caption.c_str());
1835  return hResults;
1836  }
1837 
1839  int biny)
1840  {
1842 
1843  std::vector<TH1F*> hPIDAxis =
1844  NueCCIncTemplateFitter::GetTemplates(binx,biny,"nominal");
1845 
1846  TH2F* hCovariance =
1847  new TH2F(ana::UniqueName().c_str(),//"hCovariance_holder",
1848  ";Template Bins;Template Bins;Covariance (Events^{2})",
1849  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1850  hPIDAxis[0]->GetXaxis()->GetNbins(),
1851  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1852  hPIDAxis[0]->GetXaxis()->GetNbins());
1853 
1854  for(int i = 1; i <= hCovariance->GetXaxis()->GetNbins(); i++){
1855  for(int j = 1; j <= hCovariance->GetYaxis()->GetNbins(); j++){
1856  double val =
1858  hCovariance->SetBinContent(i,j,val);
1859  }//loop over y bins
1860  }//loop over x bins
1861 
1862  std::string caption =
1864  hCovariance->SetTitle(caption.c_str());
1865 
1866  return ((TH2F*)hCovariance->Clone("hCovariance"));
1867  //ana::UniqueName().c_str()));
1868  }
1869 
1871  int biny)
1872  {
1874  std::vector<TH1F*> hPIDAxis =
1875  NueCCIncTemplateFitter::GetTemplates(binx,biny,"nominal");
1876 
1877  TH2F* hCovariance =
1878  new TH2F(ana::UniqueName().c_str(),
1879  ";Template Bins;Template Bins;Covariance (Events^{2})",
1880  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1881  hPIDAxis[0]->GetXaxis()->GetNbins(),
1882  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1883  hPIDAxis[0]->GetXaxis()->GetNbins());
1884 
1885  TH2F* hCorrelation =
1886  new TH2F(ana::UniqueName().c_str(),
1887  ";Template Bins;Template Bins;Correlation",
1888  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1889  hPIDAxis[0]->GetXaxis()->GetNbins(),
1890  hPIDAxis[0]->GetXaxis()->GetNbins(),1,
1891  hPIDAxis[0]->GetXaxis()->GetNbins());
1892 
1893  for(int i = 1; i <= hCovariance->GetXaxis()->GetNbins(); i++){
1894  for(int j = 1; j <= hCovariance->GetYaxis()->GetNbins(); j++){
1895  double val =
1897  hCovariance->SetBinContent(i,j,val);
1898  }//loop over y bins
1899  }//loop over x bins
1900 
1901  for(int i = 1; i <= hCovariance->GetXaxis()->GetNbins(); i++){
1902  for(int j = 1; j <= hCovariance->GetYaxis()->GetNbins(); j++){
1903  double val = hCovariance->GetBinContent(i,j);
1904  double dia_i = sqrt(hCovariance->GetBinContent(i,i));
1905  double dia_j = sqrt(hCovariance->GetBinContent(j,j));
1906 
1907  hCorrelation->SetBinContent(i,j,val/(dia_i*dia_j));
1908  }//loop over y bins
1909  }//loop over x bins
1910 
1911  std::string caption =
1913  hCorrelation->SetTitle(caption.c_str());
1914 
1915  return ((TH2F*)hCorrelation->Clone("hCorrelation"));
1916  }
1917 
1918 
1920  {
1921  TH2F* hResult =
1923 
1924 
1925  return hResult;
1926 
1927  }
1928 
1930  {
1931  TH2F* hCovariance =
1933 
1934  TH2F* hResult = (TH2F*)hCovariance->Clone(ana::UniqueName().c_str());
1935 
1936  for(int i = 1; i <= hCovariance->GetXaxis()->GetNbins(); i++){
1937  for(int j = 1; j <= hCovariance->GetYaxis()->GetNbins(); j++){
1938  double val = hCovariance->GetBinContent(i,j);
1939  double dia_i = sqrt(hCovariance->GetBinContent(i,i));
1940  double dia_j = sqrt(hCovariance->GetBinContent(j,j));
1941 
1942  hResult->SetBinContent(i,j,val/(dia_i*dia_j));
1943  }//loop over y bins
1944  }//loop over x bins
1945 
1946  hResult->GetZaxis()->SetTitle("Correlation");
1947 
1948  return hResult;
1949  }
1950 
1951  }//end namespace nueccinc
1952 }//end namespace ana
const XML_Char * name
Definition: expat.h:151
NueCCIncTemplateFitter(TH3F *da, std::vector< TH3F * > templates, std::vector< TH3F * > analysis, TH3F *cv, std::vector< TH3F * > systs_up, std::vector< TH3F * > systs_down, std::vector< TH3F * > systs_oneside, std::vector< TH3F * > multiverse, int numbins)
ofstream output
fvar< T > fmin(const fvar< T > &x1, const fvar< T > &x2)
Definition: fmin.hpp:14
TH2F * getCorrelationMatrixTemplateBins(int binx, int biny)
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
correl_xv GetYaxis() -> SetDecimals()
Calculate bin to bin correlation matrix from the universe technique.
const bool useSlidingSignalRegion
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
static void fcn3Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
Int_t par
Definition: SimpleIterate.C:24
double covariance(Eigen::VectorXd x, Eigen::VectorXd y)
OStream cerr
Definition: OStream.cxx:7
std::string GetAnalysisBinCaption(int binx, int biny)
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
double GetStatisticCovarianceMatrixValue(int binx, int biny)
std::vector< TH1F * > CalculateOneSigmaShift(int binx, int biny)
double CovarianceMatrix[numPIDBins][numPIDBins]
const std::string covarianceDir
int isnan(const stan::math::var &a)
Definition: std_isnan.hpp:18
double myFunction3Vars(double numu_par, double nue_par, double nc_par)
const XML_Char const XML_Char * data
Definition: expat.h:268
Double_t ymax
Definition: plot.C:25
TH2F * getCovarianceMatrixTemplateBins(int binx, int biny)
std::vector< std::pair< TH2F *, TH2F * > > getFitNormalizationAndErrors()
std::vector< std::pair< double, double > > GetTemplateWeightsAndErrors(int binx, int biny)
TH2F * getCorrelationBetweenSysts(int binx, int biny)
TH2D * GetCorrelationMatrix()
Get correlation matrix.
correl_xv GetXaxis() -> SetDecimals()
const std::string cv[Ncv]
TH2F * getCovarianceBetweenSysts(int binx, int biny)
const double j
Definition: BetheBloch.cxx:29
const float sliding_signal_region[13]
TMatrixT< double > TMatrixD
Definition: Utilities.h:18
float bin[41]
Definition: plottest35.C:14
void FillCovarianceMatrix(int binx, int biny, bool makePlots)
OStream cout
Definition: OStream.cxx:6
std::vector< TH1F * > GetTemplates(int binx, int biny, std::string name)
TH2F * FillCovarianceBetweenSysts(int binx, int biny)
std::vector< TH1F * > getPIDFitTemplates(int binx, int biny)
void FillStatisticCovarianceMatrix(int binx, int biny)
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
static void fcn2Var(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
enum BeamMode nc
c1
Definition: demo5.py:24
double GetCovarianceMatrixValue(int binx, int biny)
TLegend * AutoPlaceLegend(double dx, double dy, double yPin)
Create a legend, maximizing distance from all histograms.
Definition: Plots.cxx:715
std::vector< TH1F * > getPIDTemplates(int binx, int biny)
NueCCIncTemplateFitter * Fitter_obj
enum BeamMode kBlue
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
double StatCovarianceMatrix[numPIDBins][numPIDBins]
double myFunction2Vars(double bkgd_par, double nue_par)
enum BeamMode string
unsigned int uint