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