calculateWrongSignNue.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////
2 // //
3 // WRONG SIGN NORMALIZATION STUDY SELECTION //
4 // NUE ONLY! //
5 // //
6 /////////////////////////////////////////////////////////////
7 
10 #include "CAFAna/Fit/Fit.h"
12 #include "CAFAna/Analysis/Style.h"
13 #include "CAFAna/Vars/FitVars.h"
14 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Core/Utilities.h"
17 #include "CAFAna/Analysis/Calcs.h"
24 #include "CAFAna/Core/Progress.h"
31 #include "OscLib/OscCalcPMNSOpt.h"
32 #include "OscLib/OscCalcDumb.h"
33 #include "Utilities/func/MathUtil.h"
35 
36 using namespace ana;
37 
38 #include "TFile.h"
39 #include "TH2.h"
40 #include "TMath.h"
41 #include "TCanvas.h"
42 #include "TGraph.h"
43 #include "TLatex.h"
44 #include "TLegend.h"
45 #include "TBox.h"
46 #include "TLine.h"
47 #include "TSystem.h"
48 #include "TColor.h"
49 #include "TPaletteAxis.h"
50 #include "THStack.h"
51 #include "TString.h"
52 #include "TEllipse.h"
53 
54 #include <iostream>
55 #include <iomanip>
56 #include <string>
57 #include <vector>
58 #include "Utilities/rootlogon.C"
59 
60 TCanvas *calcAlphaBetaEachBin(Spectrum *templateSpec, std::string anaString, TH1D *wsS[], TH1D *wsB[], TH1D *wsD,
61  TH1D *rsS[], TH1D *rsB[], TH1D *rsD, bool hasSplit=false, bool treatSplit=false){
62  TH1D *alphaHist = templateSpec->ToTH1( 8.0e20 );
63  TH1D *betaHist = templateSpec->ToTH1( 8.0e20 );
64 
65  int nBinsHist = alphaHist->GetXaxis()->GetNbins();
66 
67  float bNumuMCWS[2];
68  float bNumuMCOther[2];
69  float bNumuDataCt[2];
70 
71  // Loop through bins and calculate alpha/beta
72  for( int iBinHist=1; iBinHist<=nBinsHist; ++iBinHist ){
73  bNumuMCWS[0] = 0.;
74  bNumuMCWS[1] = 0.;
75  bNumuMCOther[0] = 0.;
76  bNumuMCOther[1] = 0.;
77  bNumuDataCt[0] = 0.;
78  bNumuDataCt[1] = 0.;
79 
80  if( !hasSplit && !treatSplit ){
81  bNumuMCWS[0] = rsS[0]->GetBinContent(iBinHist) + rsS[1]->GetBinContent(iBinHist);
82  bNumuMCOther[0] = rsB[0]->GetBinContent(iBinHist) + rsB[1]->GetBinContent(iBinHist) + rsB[2]->GetBinContent(iBinHist);
83  }
84  else if( hasSplit && !treatSplit ){
85  bNumuMCWS[0] = rsS[0]->GetBinContent(iBinHist) + rsS[1]->GetBinContent(iBinHist);
86  bNumuMCOther[0] = rsB[0]->GetBinContent(iBinHist) + rsB[1]->GetBinContent(iBinHist) +
87  rsB[2]->GetBinContent(iBinHist) + rsS[2]->GetBinContent(iBinHist);
88  }
89  else if( hasSplit && treatSplit ){
90  bNumuMCWS[0] = rsS[0]->GetBinContent(iBinHist) + rsS[1]->GetBinContent(iBinHist) + rsS[2]->GetBinContent(iBinHist);
91  bNumuMCOther[0] = rsB[0]->GetBinContent(iBinHist) + rsB[1]->GetBinContent(iBinHist) + rsB[2]->GetBinContent(iBinHist);
92  }
93  else{
94  std::cout << "Doing something wrong in bin-by-bin method...";
95  bNumuMCWS[0] = 0.;
96  bNumuMCOther[0] = 0.;
97  }
98  bNumuDataCt[0] = rsD->GetBinContent(iBinHist);
99 
100  if( !hasSplit && !treatSplit ){
101  bNumuMCWS[1] = wsS[0]->GetBinContent(iBinHist) + wsS[1]->GetBinContent(iBinHist);
102  bNumuMCOther[1] = wsB[0]->GetBinContent(iBinHist) + wsB[1]->GetBinContent(iBinHist) + wsB[2]->GetBinContent(iBinHist);
103  }
104  else if( hasSplit && !treatSplit ){
105  bNumuMCWS[1] = wsS[0]->GetBinContent(iBinHist) + wsS[1]->GetBinContent(iBinHist);
106  bNumuMCOther[1] = wsB[0]->GetBinContent(iBinHist) + wsB[1]->GetBinContent(iBinHist) +
107  wsB[2]->GetBinContent(iBinHist) + wsS[2]->GetBinContent(iBinHist);
108  }
109  else if( hasSplit && treatSplit ){
110  bNumuMCWS[1] = wsS[0]->GetBinContent(iBinHist) + wsS[1]->GetBinContent(iBinHist) + wsS[2]->GetBinContent(iBinHist);
111  bNumuMCOther[1] = wsB[0]->GetBinContent(iBinHist) + wsB[1]->GetBinContent(iBinHist) + wsB[2]->GetBinContent(iBinHist);
112  }
113  else{
114  std::cout << "Doing something wrong in bin-by-bin method...";
115  bNumuMCWS[1] = 0.;
116  bNumuMCOther[1] = 0.;
117  }
118  bNumuDataCt[1] = wsD->GetBinContent(iBinHist);
119 
120  float numerator = bNumuDataCt[0] - ((bNumuMCWS[0]/bNumuMCWS[1])*bNumuDataCt[1]);
121  float denominator = bNumuMCOther[0] - ((bNumuMCWS[0]/bNumuMCWS[1])*bNumuMCOther[1]);
122  float alpha = numerator / denominator;
123  float beta = (bNumuDataCt[1] - (alpha*bNumuMCOther[1]))/bNumuMCWS[1];
124 
125  alphaHist->SetBinContent( iBinHist, alpha );
126  alphaHist->SetBinError( iBinHist, 0. );
127  betaHist->SetBinContent( iBinHist, beta );
128  betaHist->SetBinError( iBinHist, 0. );
129  }
130 
131  alphaHist->GetYaxis()->SetTitle("#alpha parameter");
132  betaHist->GetYaxis()->SetTitle("#beta parameter");
133 
134  // Print a canvas and return
135  TCanvas *tC = new TCanvas( TString::Format("binsAlphaBeta_%s",anaString.c_str()) );
136  tC->Divide(1,2);
137  tC->cd(1);
138  CenterTitles(alphaHist);
139  alphaHist->Draw();
140  Preliminary();
141  tC->cd(2);
142  CenterTitles(betaHist);
143  betaHist->Draw();
144  return tC;
145 }
146 
148 {
149  // with Data File!!
150  bool hasData = true;
151  bool splitNCe[2] = {true,true}; // first is if you HAVE split NC, second is if you want to use it
152  bool splitNCmu[2] = {true,true}; // first is if you HAVE split NC, second is if you want to use it
153  bool eachBinMu = true;
154  bool quantileMC = true;
155 
156  if( (!splitNCe[0] && splitNCe[1]) || (!splitNCmu[0] && splitNCmu[1]) ){
157  std::cout << "NC split parameter mismatch. Quitting..." << std::endl;
158  return;
159  }
160 
161  std::string strOutFile = "./WrongSignNormalizationStudy_Nue";
162  std::string strSplitNCe = "_plotNCsplitE";
163  std::string strSplitNCeMd = "_treatNCsplitE";
164  std::string strSplitNCmu = "_plotNCsplitMu";
165  std::string strSplitNCmuMd = "_treatNCsplitMu";
166  std::string strEachBinMu = "_withEachBinMu";
167  std::string strQuantileMC = "_QuantilesInMC";
168  if( splitNCe[0] && !splitNCe[1] ) strOutFile += strSplitNCe;
169  if( splitNCe[0] && splitNCe[1] ) strOutFile += strSplitNCeMd;
170  if( splitNCmu[0] && !splitNCmu[1] ) strOutFile += strSplitNCmu;
171  if( splitNCmu[0] && splitNCmu[1] ) strOutFile += strSplitNCmuMd;
172  if( eachBinMu ) strOutFile += strEachBinMu;
173  if( quantileMC ) strOutFile += strQuantileMC;
174  TFile *outFile = new TFile( TString::Format("%s.root",strOutFile.c_str()), "recreate" );
175 
176  TFile *inFile = new TFile("/nova/ana/users/howard/Nue2018Ana/development/Predictions_And_Data-WSBinning-NCSplit-26March2018/predictions_and_data_WrongsignRHC_splitNC.root","read");
177 
179  std::vector< std::string > anaTitles = {"FSProton","CVNProng","BDT3"};
180  std::vector< std::string > nueModes = {"HiCVN","LoCVN"};
181 
182  // ----- Load files -----
183  PredictionNoOsc *pRHCNue[3][2];
184  Spectrum *sNCNue[3][2];
185  Spectrum *sNCAntiNue[3][2];
186 
187  Spectrum *sRHCNueData[3][2];
188 
189  for( int iAna=0; iAna<3; ++iAna ){
190  for( int iMd=0; iMd<2; ++iMd ){
191  pRHCNue[iAna][iMd] = PredictionNoOsc::LoadFrom(inFile, TString::Format("pNue%s%s",
192  anaTitles[iAna].c_str(),
193  nueModes[iMd].c_str()).Data() ).release();
194  if( splitNCe[0] ){
195  sNCNue[iAna][iMd] = Spectrum::LoadFrom(inFile, TString::Format("sNue%s%sNC",
196  anaTitles[iAna].c_str(),
197  nueModes[iMd].c_str()).Data() ).release();
198  sNCAntiNue[iAna][iMd] = Spectrum::LoadFrom(inFile, TString::Format("sNue%s%sNCAnti",
199  anaTitles[iAna].c_str(),
200  nueModes[iMd].c_str()).Data() ).release();
201  }
202  }
203  if( hasData ){
204  for( int iMd=0; iMd<2; ++iMd ){
205  sRHCNueData[iAna][iMd] = Spectrum::LoadFrom(inFile, TString::Format("sNueData_%s%s",
206  anaTitles[iAna].c_str(),
207  nueModes[iMd].c_str()).Data() ).release();
208  }
209  }
210  }
211 
212  delete inFile;
213 
214  // Get POT and if pulled from data check it against other spectra
215  float pot = 8e20; // replace this with data POT when we have that...
216  if( hasData ){
217  for( int iAna=0; iAna<3; ++iAna ){
218  if( int(sRHCNueData[iAna][0]->POT()) != int(sRHCNueData[iAna][1]->POT()) ){
219  std::cout << "Low and High CVN electron score bins don't match" << std::endl;
220  return;
221  }
222  if( iAna>0 && int(sRHCNueData[iAna][0]->POT()) != int(sRHCNueData[iAna-1][0]->POT()) ){
223  std::cout << "POT information between " << iAna << " and " << iAna-1 << " methods did not match." << std::endl;
224  return;
225  };
226  }
227  pot = sRHCNueData[0][0]->POT();
228  }
229 
230  std::cout << "----------------------------" << std::endl;
231  std::cout << " POT: " << pot << std::endl;
232  std::cout << "----------------------------" << std::endl;
233 
234  // ----- Nue -----
235  TCanvas *cResults[3][2];
236  THStack *stackAnaBins[3][2];
237  TH1D *nueS[3][3][2]; // make this also have 3 categories, so I can handle NC as part of signal (and NCAnti alone as part of BG)
238  TH1D *nueB[3][3][2];
239  TH1D *nueAll[3][2];
240 
241  TH1D *nueDataResults[3][2];
242  TLegend *nueLegend[3][2];
243 
244  for( int iAna=0; iAna<3; ++iAna ){
245  for( int iMd=0; iMd<2; ++iMd ){
246  // Get predictions
247  nueS[iAna][0][iMd] = pRHCNue[iAna][iMd]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUE
248  nueS[iAna][0][iMd]->SetName( TString::Format("Nue%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()) );
249  nueS[iAna][1][iMd] = pRHCNue[iAna][iMd]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUMU
250  nueS[iAna][1][iMd]->SetName( TString::Format("Numu%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()) );
251  nueB[iAna][0][iMd] = pRHCNue[iAna][iMd]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUE
252  nueB[iAna][0][iMd]->SetName( TString::Format("Nuebar%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()) );
253  nueB[iAna][1][iMd] = pRHCNue[iAna][iMd]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUMU
254  nueB[iAna][1][iMd]->SetName( TString::Format("Numubar%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()) );
255  if( !splitNCe[0] ){
256  nueB[iAna][2][iMd] = pRHCNue[iAna][iMd]->PredictComponent(&calc,Flavors::kAll,Current::kNC,Sign::kBoth).ToTH1(pot); // NC
257  nueB[iAna][2][iMd]->SetName( TString::Format("Nc%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()) );
258  }
259  // if splitNCe then handle differently
260  else{
261  nueB[iAna][2][iMd] = sNCAntiNue[iAna][iMd]->ToTH1(pot); // NCAnti
262  nueB[iAna][2][iMd]->SetName( TString::Format("NcAnti%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()) );
263  nueS[iAna][2][iMd] = sNCNue[iAna][iMd]->ToTH1(pot); // NC
264  nueS[iAna][2][iMd]->SetName( TString::Format("Nc%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()) );
265  }
266  // Set Colors
267  nueS[iAna][0][iMd]->SetLineColor(kGreen);
268  nueS[iAna][0][iMd]->SetFillColor(kGreen);
269  nueS[iAna][1][iMd]->SetLineColor(kOrange);
270  nueS[iAna][1][iMd]->SetFillColor(kOrange);
271  nueB[iAna][0][iMd]->SetLineColor(kGreen+2);
272  nueB[iAna][0][iMd]->SetFillColor(kGreen+2);
273  nueB[iAna][1][iMd]->SetLineColor(kRed+1);
274  nueB[iAna][1][iMd]->SetFillColor(kRed+1);
275  nueB[iAna][2][iMd]->SetLineColor(kMagenta);
276  nueB[iAna][2][iMd]->SetFillColor(kMagenta);
277  if( splitNCe[0] ){
278  nueS[iAna][2][iMd]->SetLineColor(kMagenta);
279  nueS[iAna][2][iMd]->SetFillColor(kMagenta);
280  nueB[iAna][2][iMd]->SetLineColor(kGray+2);
281  nueB[iAna][2][iMd]->SetFillColor(kGray+2);
282  }
283  // Fill up the sum histogram
284  nueAll[iAna][iMd] = new TH1D( TString::Format("All%s_%s",
285  anaTitles[iAna].c_str(),
286  nueModes[iMd].c_str()),";Reco Energy / Ana Bin;Events",18,0,18 );
287  nueAll[iAna][iMd]->Sumw2();
288  nueAll[iAna][iMd]->Add(nueS[iAna][0][iMd]);
289  nueAll[iAna][iMd]->Add(nueS[iAna][1][iMd]);
290  nueAll[iAna][iMd]->Add(nueB[iAna][0][iMd]);
291  nueAll[iAna][iMd]->Add(nueB[iAna][1][iMd]);
292  nueAll[iAna][iMd]->Add(nueB[iAna][2][iMd]);
293  // Get Data!
294  if( hasData ) nueDataResults[iAna][iMd] = sRHCNueData[iAna][iMd]->ToTH1( pot );
295  }
296  }
297 
298  double maxY[2] = {1.2*std::max( std::max(nueAll[0][0]->GetMaximum(),nueAll[1][0]->GetMaximum()), nueAll[2][0]->GetMaximum()),
299  1.2*std::max( std::max(nueAll[0][1]->GetMaximum(),nueAll[1][1]->GetMaximum()), nueAll[2][1]->GetMaximum())};
300 
301  if( hasData ){
302  maxY[0] = 1.2*std::max( std::max( std::max( std::max( std::max(
303  nueAll[0][0]->GetMaximum(),nueAll[1][0]->GetMaximum()),
304  nueAll[2][0]->GetMaximum()),
305  nueDataResults[0][0]->GetMaximum()),
306  nueDataResults[1][0]->GetMaximum()),
307  nueDataResults[2][0]->GetMaximum());
308  maxY[1] = 1.2*std::max( std::max( std::max( std::max( std::max(
309  nueAll[0][1]->GetMaximum(),nueAll[1][1]->GetMaximum()),
310  nueAll[2][1]->GetMaximum()),
311  nueDataResults[0][1]->GetMaximum()),
312  nueDataResults[1][1]->GetMaximum()),
313  nueDataResults[2][1]->GetMaximum());
314  }
315 
316  TH2F *blankH[2];
317  blankH[0] = new TH2F("blankH0",";NuE Energy / Analysis Bin;Events",18,0,18,int(maxY[0])+1,0,float(int(maxY[0])+1));
318  blankH[1] = new TH2F("blankH1",";NuE Energy / Analysis Bin;Events",18,0,18,int(maxY[1])+1,0,float(int(maxY[1])+1));
319  CenterTitles(blankH[0]);
320  CenterTitles(blankH[1]);
321 
322  TLine *bin1[2];
323  bin1[0] = new TLine(9,0.,9,float(int(maxY[0])+1));
324  bin1[1] = new TLine(9,0.,9,float(int(maxY[1])+1));
325  bin1[0]->SetLineStyle(7);
326  bin1[0]->SetLineColor(kBlack);
327  bin1[1]->SetLineStyle(7);
328  bin1[1]->SetLineColor(kBlack);
329 
330  // dataType[method 0-2][mode 0-1 (high/lo CVN)][WS suppressed (0), WS enhanced (1)]
331  float NueMCWS[3][2][2];
332  float NueMCOther[3][2][2];
333  float NueDataCt[3][2][2];
334 
335  outFile->cd();
336 
337  for( int iAna=0; iAna<3; ++iAna ){
338  for( int iMd=0; iMd<2; ++iMd ){
339  cResults[iAna][iMd] = new TCanvas();
340  cResults[iAna][iMd]->SetRightMargin(0.02);
341  stackAnaBins[iAna][iMd] = new THStack(TString::Format("stackAnaBins%i_%i",iAna,iMd),"");
342  stackAnaBins[iAna][iMd]->Add(nueB[iAna][2][iMd]);
343  if( splitNCe[0] ){
344  stackAnaBins[iAna][iMd]->Add(nueS[iAna][2][iMd]);
345  }
346  stackAnaBins[iAna][iMd]->Add(nueB[iAna][1][iMd]);
347  stackAnaBins[iAna][iMd]->Add(nueS[iAna][1][iMd]);
348  stackAnaBins[iAna][iMd]->Add(nueS[iAna][0][iMd]);
349  stackAnaBins[iAna][iMd]->Add(nueB[iAna][0][iMd]);
350  blankH[iMd]->Draw("");
351  stackAnaBins[iAna][iMd]->Draw("hist same");
352  if( hasData ) nueDataResults[iAna][iMd]->Draw("same");
353  nueLegend[iAna][iMd] = new TLegend(0.7,0.63,0.88,0.88);
354  nueLegend[iAna][iMd]->AddEntry( TString::Format("Nue%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()),"#nu_{e}","lf" );
355  nueLegend[iAna][iMd]->AddEntry( TString::Format("Numu%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()),"#nu_{#mu}","lf" );
356  nueLegend[iAna][iMd]->AddEntry( TString::Format("Nuebar%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()),"#bar{#nu}_{e}","lf" );
357  nueLegend[iAna][iMd]->AddEntry( TString::Format("Numubar%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()),"#bar{#nu}_{#mu}","lf" );
358  nueLegend[iAna][iMd]->AddEntry( TString::Format("Nc%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()),"#nu NC","lf" );
359  if( splitNCe[0] )
360  nueLegend[iAna][iMd]->AddEntry( TString::Format("NcAnti%s_%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()),"#bar{#nu} NC","lf" );
361  nueLegend[iAna][iMd]->Draw();
362  bin1[iMd]->Draw();
363  if( !hasData ) Simulation();
364  if( hasData ) Preliminary();
365  cResults[iAna][iMd]->Write( TString::Format("cRHCNue%s%s",anaTitles[iAna].c_str(),nueModes[iMd].c_str()) );
366 
367  // Values
368  if( !splitNCe[0] && !splitNCe[1] ){
369  NueMCWS[iAna][iMd][0] = nueS[iAna][0][iMd]->Integral(1,9) + nueS[iAna][1][iMd]->Integral(1,9);
370  NueMCOther[iAna][iMd][0] = nueB[iAna][0][iMd]->Integral(1,9) + nueB[iAna][1][iMd]->Integral(1,9) + nueB[iAna][2][iMd]->Integral(1,9);
371  }
372  else if( splitNCe[0] && !splitNCe[1] ){
373  NueMCWS[iAna][iMd][0] = nueS[iAna][0][iMd]->Integral(1,9) + nueS[iAna][1][iMd]->Integral(1,9);
374  NueMCOther[iAna][iMd][0] = nueB[iAna][0][iMd]->Integral(1,9) + nueB[iAna][1][iMd]->Integral(1,9) +
375  nueB[iAna][2][iMd]->Integral(1,9) + nueS[iAna][2][iMd]->Integral(1,9);
376  }
377  else if( splitNCe[0] && splitNCe[1] ){
378  NueMCWS[iAna][iMd][0] = nueS[iAna][0][iMd]->Integral(1,9) + nueS[iAna][1][iMd]->Integral(1,9) + nueS[iAna][2][iMd]->Integral(1,9);
379  NueMCOther[iAna][iMd][0] = nueB[iAna][0][iMd]->Integral(1,9) + nueB[iAna][1][iMd]->Integral(1,9) + nueB[iAna][2][iMd]->Integral(1,9);
380  }
381  else{
382  std::cout << "IMPROPER NUE NC MODE SELECTED. WILL FILL WITH 0s..." << std::endl;
383  NueMCWS[iAna][iMd][0] = 0.;
384  NueMCOther[iAna][iMd][0] = 0.;
385  }
386  if( hasData ) NueDataCt[iAna][iMd][0] = nueDataResults[iAna][iMd]->Integral(1,9);
387 
388  if( !splitNCe[0] && !splitNCe[1] ){
389  NueMCWS[iAna][iMd][1] = nueS[iAna][0][iMd]->Integral(10,18) + nueS[iAna][1][iMd]->Integral(10,18);
390  NueMCOther[iAna][iMd][1] = nueB[iAna][0][iMd]->Integral(10,18) + nueB[iAna][1][iMd]->Integral(10,18) +
391  nueB[iAna][2][iMd]->Integral(10,18);
392  }
393  else if( splitNCe[0] && !splitNCe[1] ){
394  NueMCWS[iAna][iMd][1] = nueS[iAna][0][iMd]->Integral(10,18) + nueS[iAna][1][iMd]->Integral(10,18);
395  NueMCOther[iAna][iMd][1] = nueB[iAna][0][iMd]->Integral(10,18) + nueB[iAna][1][iMd]->Integral(10,18) +
396  nueB[iAna][2][iMd]->Integral(10,18) + nueS[iAna][2][iMd]->Integral(10,18);
397  }
398  else if( splitNCe[0] && splitNCe[1] ){
399  NueMCWS[iAna][iMd][1] = nueS[iAna][0][iMd]->Integral(10,18) + nueS[iAna][1][iMd]->Integral(10,18) +
400  nueS[iAna][2][iMd]->Integral(10,18);
401  NueMCOther[iAna][iMd][1] = nueB[iAna][0][iMd]->Integral(10,18) + nueB[iAna][1][iMd]->Integral(10,18) +
402  nueB[iAna][2][iMd]->Integral(10,18);
403  }
404  else{
405  std::cout << "IMPROPER NUE NC MODE SELECTED. WILL FILL WITH 0s..." << std::endl;
406  NueMCWS[iAna][iMd][1] = 0.;
407  NueMCOther[iAna][iMd][1] = 0.;
408  }
409  if( hasData ) NueDataCt[iAna][iMd][1] = nueDataResults[iAna][iMd]->Integral(10,18);
410 
411  std::cout << "Nue " << anaTitles[iAna] << "(" << nueModes[iMd] << "):" << std::endl;
412  std::cout << " __Right-sign enhanced__" << std::endl;
413  std::cout << " signal (WS) = " << NueMCWS[iAna][iMd][0] << std::endl;
414  std::cout << " bgd (RS+NC) = " << NueMCOther[iAna][iMd][0] << std::endl;
415  if( hasData ) std::cout << " DATA = " << NueDataCt[iAna][iMd][0] << std::endl;
416  std::cout << std::endl;
417  std::cout << " __Wrong-sign enhanced__" << std::endl;
418  std::cout << " signal (WS) = " << NueMCWS[iAna][iMd][1] << std::endl;
419  std::cout << " bgd (RS+NC) = " << NueMCOther[iAna][iMd][1] << std::endl;
420  if( hasData ) std::cout << " DATA = " << NueDataCt[iAna][iMd][1] << std::endl;
421  std::cout << std::endl;
422  if( splitNCe[0] && splitNCe[1] ){
423  std::cout << "WSCC: " << nueS[iAna][0][iMd]->Integral() + nueS[iAna][1][iMd]->Integral() << std::endl;
424  std::cout << "WSNC: " << nueS[iAna][2][iMd]->Integral() << std::endl;
425  std::cout << "RSCC: " << nueB[iAna][0][iMd]->Integral() + nueB[iAna][1][iMd]->Integral() << std::endl;
426  std::cout << "RSNC: " << nueB[iAna][2][iMd]->Integral() << std::endl;
427  std::cout << std::endl;
428  std::cout << "Nue (RS+WS): " << nueS[iAna][0][iMd]->Integral() + nueB[iAna][0][iMd]->Integral() << std::endl;
429  std::cout << "Numu (RS+WS): " << nueS[iAna][1][iMd]->Integral() + nueB[iAna][1][iMd]->Integral() << std::endl;
430  std::cout << "NC (RS+WS): " << nueS[iAna][2][iMd]->Integral() + nueB[iAna][2][iMd]->Integral() << std::endl;
431  }
432  std::cout << std::endl;
433  std::cout << std::endl;
434  }
435  }
436 
437  // Solve System of equations for all three methods (if hasData is turned on)
438  TEllipse *loCVNEllipse[2];
439  TEllipse *hiCVNEllipse[2];
440  TEllipse *loCVNEllipsePt[2];
441  TEllipse *hiCVNEllipsePt[2];
442 
443  if( hasData ){
444  // Low CVN
445  // FSProton
446  float numeratorFSProtonLoCVN = NueDataCt[0][1][0] - ((NueMCWS[0][1][0]/NueMCWS[0][1][1])*NueDataCt[0][1][1]);
447  float denominatorFSProtonLoCVN = NueMCOther[0][1][0] - ((NueMCWS[0][1][0]/NueMCWS[0][1][1])*NueMCOther[0][1][1]);
448  float alphaFSProtonLoCVN = numeratorFSProtonLoCVN / denominatorFSProtonLoCVN;
449  float betaFSProtonLoCVN = (NueDataCt[0][1][1] - (alphaFSProtonLoCVN*NueMCOther[0][1][1]))/NueMCWS[0][1][1];
450  float denomCovarFSProtonLoCVN = (NueMCOther[0][1][0]*NueMCWS[0][1][1]) - (NueMCOther[0][1][1]*NueMCWS[0][1][0]);
451  float A11FSProtonLoCVN = NueMCWS[0][1][1] / denomCovarFSProtonLoCVN;
452  float A12FSProtonLoCVN = -NueMCWS[0][1][0] / denomCovarFSProtonLoCVN;
453  float A21FSProtonLoCVN = -NueMCOther[0][1][1] / denomCovarFSProtonLoCVN;
454  float A22FSProtonLoCVN = NueMCOther[0][1][0] / denomCovarFSProtonLoCVN;
455  float C11FSProtonLoCVN = (A11FSProtonLoCVN*A11FSProtonLoCVN*NueDataCt[0][1][0]) + (A12FSProtonLoCVN*A12FSProtonLoCVN*NueDataCt[0][1][1]);
456  float C12FSProtonLoCVN = (A11FSProtonLoCVN*A21FSProtonLoCVN*NueDataCt[0][1][0]) + (A12FSProtonLoCVN*A22FSProtonLoCVN*NueDataCt[0][1][1]);
457  float C21FSProtonLoCVN = C12FSProtonLoCVN;
458  float C22FSProtonLoCVN = (A21FSProtonLoCVN*A21FSProtonLoCVN*NueDataCt[0][1][0]) + (A22FSProtonLoCVN*A22FSProtonLoCVN*NueDataCt[0][1][1]);
459  float majorFSProtonLoCVN = sqrt(0.5*(C11FSProtonLoCVN + C22FSProtonLoCVN +
460  sqrt( (C11FSProtonLoCVN+C22FSProtonLoCVN)*(C11FSProtonLoCVN+C22FSProtonLoCVN)
461  -(4*( C11FSProtonLoCVN*C22FSProtonLoCVN - C12FSProtonLoCVN*C21FSProtonLoCVN )))));
462  float minorFSProtonLoCVN = sqrt(0.5*(C11FSProtonLoCVN + C22FSProtonLoCVN -
463  sqrt( (C11FSProtonLoCVN+C22FSProtonLoCVN)*(C11FSProtonLoCVN+C22FSProtonLoCVN)
464  -(4*( C11FSProtonLoCVN*C22FSProtonLoCVN - C12FSProtonLoCVN*C21FSProtonLoCVN )))));
465  float rotBetaXFSProtonLoCVN = 0.5*(180./TMath::Pi())*atan2(2.0*C12FSProtonLoCVN,C22FSProtonLoCVN - C11FSProtonLoCVN);
466  float wsfracFSProtonLoCVN =
467  (betaFSProtonLoCVN*(NueMCWS[0][1][0]+NueMCWS[0][1][1])) /
468  (betaFSProtonLoCVN*(NueMCWS[0][1][0]+NueMCWS[0][1][1]) + alphaFSProtonLoCVN*(NueMCOther[0][1][0]+NueMCOther[0][1][1]));
469  float wsfracFSProtonLoCVNErrDenom =
470  (betaFSProtonLoCVN*(NueMCWS[0][1][0]+NueMCWS[0][1][1]) + alphaFSProtonLoCVN*(NueMCOther[0][1][0]+NueMCOther[0][1][1])) *
471  (betaFSProtonLoCVN*(NueMCWS[0][1][0]+NueMCWS[0][1][1]) + alphaFSProtonLoCVN*(NueMCOther[0][1][0]+NueMCOther[0][1][1])) *
472  (betaFSProtonLoCVN*(NueMCWS[0][1][0]+NueMCWS[0][1][1]) + alphaFSProtonLoCVN*(NueMCOther[0][1][0]+NueMCOther[0][1][1])) *
473  (betaFSProtonLoCVN*(NueMCWS[0][1][0]+NueMCWS[0][1][1]) + alphaFSProtonLoCVN*(NueMCOther[0][1][0]+NueMCOther[0][1][1])); // term^4
474  float wsfracFSProtonLoCVNErrAB2 =
475  ( (NueMCWS[0][1][0]+NueMCWS[0][1][1])*(NueMCOther[0][1][0]+NueMCOther[0][1][1]) ) *
476  ( (NueMCWS[0][1][0]+NueMCWS[0][1][1])*(NueMCOther[0][1][0]+NueMCOther[0][1][1]) ); //AB^2
477  float wsfracFSProtonLoCVNErr =
478  sqrt( C11FSProtonLoCVN*wsfracFSProtonLoCVNErrAB2*betaFSProtonLoCVN*betaFSProtonLoCVN/wsfracFSProtonLoCVNErrDenom +
479  C22FSProtonLoCVN*wsfracFSProtonLoCVNErrAB2*alphaFSProtonLoCVN*alphaFSProtonLoCVN/wsfracFSProtonLoCVNErrDenom -
480  2*C12FSProtonLoCVN*wsfracFSProtonLoCVNErrAB2*alphaFSProtonLoCVN*betaFSProtonLoCVN/wsfracFSProtonLoCVNErrDenom );
481  float betaAlphaFSProtonLoCVN = betaFSProtonLoCVN / alphaFSProtonLoCVN;
482  float betaAlphaFSProtonLoCVNErr = sqrt(
483  C11FSProtonLoCVN*betaFSProtonLoCVN*betaFSProtonLoCVN/(alphaFSProtonLoCVN*alphaFSProtonLoCVN*alphaFSProtonLoCVN*alphaFSProtonLoCVN)+
484  C22FSProtonLoCVN*(1./alphaFSProtonLoCVN)*(1./alphaFSProtonLoCVN)-
485  2*C12FSProtonLoCVN*betaFSProtonLoCVN/(alphaFSProtonLoCVN*alphaFSProtonLoCVN*alphaFSProtonLoCVN) );
486  // CVNProng
487  float numeratorCVNProngLoCVN = NueDataCt[1][1][0] - ((NueMCWS[1][1][0]/NueMCWS[1][1][1])*NueDataCt[1][1][1]);
488  float denominatorCVNProngLoCVN = NueMCOther[1][1][0] - ((NueMCWS[1][1][0]/NueMCWS[1][1][1])*NueMCOther[1][1][1]);
489  float alphaCVNProngLoCVN = numeratorCVNProngLoCVN / denominatorCVNProngLoCVN;
490  float betaCVNProngLoCVN = (NueDataCt[1][1][1] - (alphaCVNProngLoCVN*NueMCOther[1][1][1]))/NueMCWS[1][1][1];
491  float denomCovarCVNProngLoCVN = (NueMCOther[1][1][0]*NueMCWS[1][1][1]) - (NueMCOther[1][1][1]*NueMCWS[1][1][0]);
492  float A11CVNProngLoCVN = NueMCWS[1][1][1] / denomCovarCVNProngLoCVN;
493  float A12CVNProngLoCVN = -NueMCWS[1][1][0] / denomCovarCVNProngLoCVN;
494  float A21CVNProngLoCVN = -NueMCOther[1][1][1] / denomCovarCVNProngLoCVN;
495  float A22CVNProngLoCVN = NueMCOther[1][1][0] / denomCovarCVNProngLoCVN;
496  float C11CVNProngLoCVN = (A11CVNProngLoCVN*A11CVNProngLoCVN*NueDataCt[1][1][0]) + (A12CVNProngLoCVN*A12CVNProngLoCVN*NueDataCt[1][1][1]);
497  float C12CVNProngLoCVN = (A11CVNProngLoCVN*A21CVNProngLoCVN*NueDataCt[1][1][0]) + (A12CVNProngLoCVN*A22CVNProngLoCVN*NueDataCt[1][1][1]);
498  float C21CVNProngLoCVN = C12CVNProngLoCVN;
499  float C22CVNProngLoCVN = (A21CVNProngLoCVN*A21CVNProngLoCVN*NueDataCt[1][1][0]) + (A22CVNProngLoCVN*A22CVNProngLoCVN*NueDataCt[1][1][1]);
500  float majorCVNProngLoCVN = sqrt(0.5*(C11CVNProngLoCVN + C22CVNProngLoCVN +
501  sqrt( (C11CVNProngLoCVN+C22CVNProngLoCVN)*(C11CVNProngLoCVN+C22CVNProngLoCVN)
502  -(4*( C11CVNProngLoCVN*C22CVNProngLoCVN - C12CVNProngLoCVN*C21CVNProngLoCVN )))));
503  float minorCVNProngLoCVN = sqrt(0.5*(C11CVNProngLoCVN + C22CVNProngLoCVN -
504  sqrt( (C11CVNProngLoCVN+C22CVNProngLoCVN)*(C11CVNProngLoCVN+C22CVNProngLoCVN)
505  -(4*( C11CVNProngLoCVN*C22CVNProngLoCVN - C12CVNProngLoCVN*C21CVNProngLoCVN )))));
506  float rotBetaXCVNProngLoCVN = 0.5*(180./TMath::Pi())*atan2(2.0*C12CVNProngLoCVN,C22CVNProngLoCVN - C11CVNProngLoCVN);
507  float wsfracCVNProngLoCVN =
508  (betaCVNProngLoCVN*(NueMCWS[1][1][0]+NueMCWS[1][1][1])) /
509  (betaCVNProngLoCVN*(NueMCWS[1][1][0]+NueMCWS[1][1][1]) + alphaCVNProngLoCVN*(NueMCOther[1][1][0]+NueMCOther[1][1][1]));
510  float wsfracCVNProngLoCVNErrDenom =
511  (betaCVNProngLoCVN*(NueMCWS[1][1][0]+NueMCWS[1][1][1]) + alphaCVNProngLoCVN*(NueMCOther[1][1][0]+NueMCOther[1][1][1])) *
512  (betaCVNProngLoCVN*(NueMCWS[1][1][0]+NueMCWS[1][1][1]) + alphaCVNProngLoCVN*(NueMCOther[1][1][0]+NueMCOther[1][1][1])) *
513  (betaCVNProngLoCVN*(NueMCWS[1][1][0]+NueMCWS[1][1][1]) + alphaCVNProngLoCVN*(NueMCOther[1][1][0]+NueMCOther[1][1][1])) *
514  (betaCVNProngLoCVN*(NueMCWS[1][1][0]+NueMCWS[1][1][1]) + alphaCVNProngLoCVN*(NueMCOther[1][1][0]+NueMCOther[1][1][1])); // term^4
515  float wsfracCVNProngLoCVNErrAB2 =
516  ( (NueMCWS[1][1][0]+NueMCWS[1][1][1])*(NueMCOther[1][1][0]+NueMCOther[1][1][1]) ) *
517  ( (NueMCWS[1][1][0]+NueMCWS[1][1][1])*(NueMCOther[1][1][0]+NueMCOther[1][1][1]) ); //AB^2
518  float wsfracCVNProngLoCVNErr =
519  sqrt( C11CVNProngLoCVN*wsfracCVNProngLoCVNErrAB2*betaCVNProngLoCVN*betaCVNProngLoCVN/wsfracCVNProngLoCVNErrDenom +
520  C22CVNProngLoCVN*wsfracCVNProngLoCVNErrAB2*alphaCVNProngLoCVN*alphaCVNProngLoCVN/wsfracCVNProngLoCVNErrDenom -
521  2*C12CVNProngLoCVN*wsfracCVNProngLoCVNErrAB2*alphaCVNProngLoCVN*betaCVNProngLoCVN/wsfracCVNProngLoCVNErrDenom );
522  float betaAlphaCVNProngLoCVN = betaCVNProngLoCVN / alphaCVNProngLoCVN;
523  float betaAlphaCVNProngLoCVNErr = sqrt(
524  C11CVNProngLoCVN*betaCVNProngLoCVN*betaCVNProngLoCVN/(alphaCVNProngLoCVN*alphaCVNProngLoCVN*alphaCVNProngLoCVN*alphaCVNProngLoCVN)+
525  C22CVNProngLoCVN*(1./alphaCVNProngLoCVN)*(1./alphaCVNProngLoCVN)-
526  2*C12CVNProngLoCVN*betaCVNProngLoCVN/(alphaCVNProngLoCVN*alphaCVNProngLoCVN*alphaCVNProngLoCVN) );
527  // BDT3
528  float numeratorBDT3LoCVN = NueDataCt[2][1][0] - ((NueMCWS[2][1][0]/NueMCWS[2][1][1])*NueDataCt[2][1][1]);
529  float denominatorBDT3LoCVN = NueMCOther[2][1][0] - ((NueMCWS[2][1][0]/NueMCWS[2][1][1])*NueMCOther[2][1][1]);
530  float alphaBDT3LoCVN = numeratorBDT3LoCVN / denominatorBDT3LoCVN;
531  float betaBDT3LoCVN = (NueDataCt[2][1][1] - (alphaBDT3LoCVN*NueMCOther[2][1][1]))/NueMCWS[2][1][1];
532  float denomCovarBDT3LoCVN = (NueMCOther[2][1][0]*NueMCWS[2][1][1]) - (NueMCOther[2][1][1]*NueMCWS[2][1][0]);
533  float A11BDT3LoCVN = NueMCWS[2][1][1] / denomCovarBDT3LoCVN;
534  float A12BDT3LoCVN = -NueMCWS[2][1][0] / denomCovarBDT3LoCVN;
535  float A21BDT3LoCVN = -NueMCOther[2][1][1] / denomCovarBDT3LoCVN;
536  float A22BDT3LoCVN = NueMCOther[2][1][0] / denomCovarBDT3LoCVN;
537  float C11BDT3LoCVN = (A11BDT3LoCVN*A11BDT3LoCVN*NueDataCt[2][1][0]) + (A12BDT3LoCVN*A12BDT3LoCVN*NueDataCt[2][1][1]);
538  float C12BDT3LoCVN = (A11BDT3LoCVN*A21BDT3LoCVN*NueDataCt[2][1][0]) + (A12BDT3LoCVN*A22BDT3LoCVN*NueDataCt[2][1][1]);
539  float C21BDT3LoCVN = C12BDT3LoCVN;
540  float C22BDT3LoCVN = (A21BDT3LoCVN*A21BDT3LoCVN*NueDataCt[2][1][0]) + (A22BDT3LoCVN*A22BDT3LoCVN*NueDataCt[2][1][1]);
541  float majorBDT3LoCVN = sqrt(0.5*(C11BDT3LoCVN + C22BDT3LoCVN +
542  sqrt( (C11BDT3LoCVN+C22BDT3LoCVN)*(C11BDT3LoCVN+C22BDT3LoCVN)
543  -(4*( C11BDT3LoCVN*C22BDT3LoCVN - C12BDT3LoCVN*C21BDT3LoCVN )))));
544  float minorBDT3LoCVN = sqrt(0.5*(C11BDT3LoCVN + C22BDT3LoCVN -
545  sqrt( (C11BDT3LoCVN+C22BDT3LoCVN)*(C11BDT3LoCVN+C22BDT3LoCVN)
546  -(4*( C11BDT3LoCVN*C22BDT3LoCVN - C12BDT3LoCVN*C21BDT3LoCVN )))));
547  float rotBetaXBDT3LoCVN = 0.5*(180./TMath::Pi())*atan2(2.0*C12BDT3LoCVN,C22BDT3LoCVN - C11BDT3LoCVN);
548  float wsfracBDT3LoCVN =
549  (betaBDT3LoCVN*(NueMCWS[2][1][0]+NueMCWS[2][1][1])) /
550  (betaBDT3LoCVN*(NueMCWS[2][1][0]+NueMCWS[2][1][1]) + alphaBDT3LoCVN*(NueMCOther[2][1][0]+NueMCOther[2][1][1]));
551  float wsfracBDT3LoCVNErrDenom =
552  (betaBDT3LoCVN*(NueMCWS[2][1][0]+NueMCWS[2][1][1]) + alphaBDT3LoCVN*(NueMCOther[2][1][0]+NueMCOther[2][1][1])) *
553  (betaBDT3LoCVN*(NueMCWS[2][1][0]+NueMCWS[2][1][1]) + alphaBDT3LoCVN*(NueMCOther[2][1][0]+NueMCOther[2][1][1])) *
554  (betaBDT3LoCVN*(NueMCWS[2][1][0]+NueMCWS[2][1][1]) + alphaBDT3LoCVN*(NueMCOther[2][1][0]+NueMCOther[2][1][1])) *
555  (betaBDT3LoCVN*(NueMCWS[2][1][0]+NueMCWS[2][1][1]) + alphaBDT3LoCVN*(NueMCOther[2][1][0]+NueMCOther[2][1][1])); // term^4
556  float wsfracBDT3LoCVNErrAB2 =
557  ( (NueMCWS[2][1][0]+NueMCWS[2][1][1])*(NueMCOther[2][1][0]+NueMCOther[2][1][1]) ) *
558  ( (NueMCWS[2][1][0]+NueMCWS[2][1][1])*(NueMCOther[2][1][0]+NueMCOther[2][1][1]) ); //AB^2
559  float wsfracBDT3LoCVNErr =
560  sqrt( C11BDT3LoCVN*wsfracBDT3LoCVNErrAB2*betaBDT3LoCVN*betaBDT3LoCVN/wsfracBDT3LoCVNErrDenom +
561  C22BDT3LoCVN*wsfracBDT3LoCVNErrAB2*alphaBDT3LoCVN*alphaBDT3LoCVN/wsfracBDT3LoCVNErrDenom -
562  2*C12BDT3LoCVN*wsfracBDT3LoCVNErrAB2*alphaBDT3LoCVN*betaBDT3LoCVN/wsfracBDT3LoCVNErrDenom );
563  float betaAlphaBDT3LoCVN = betaBDT3LoCVN / alphaBDT3LoCVN;
564  float betaAlphaBDT3LoCVNErr = sqrt(
565  C11BDT3LoCVN*betaBDT3LoCVN*betaBDT3LoCVN/(alphaBDT3LoCVN*alphaBDT3LoCVN*alphaBDT3LoCVN*alphaBDT3LoCVN)+
566  C22BDT3LoCVN*(1./alphaBDT3LoCVN)*(1./alphaBDT3LoCVN)-
567  2*C12BDT3LoCVN*betaBDT3LoCVN/(alphaBDT3LoCVN*alphaBDT3LoCVN*alphaBDT3LoCVN) );
568  std::cout << "==========================================" << std::endl;
569  std::cout << " SYSTEM OF EQUATIONS RESULTS: NUE " << std::endl;
570  std::cout << " LO CVN BIN (0.89 <= CVNSSe < 0.98) " << std::endl;
571  std::cout << "------------------------------------------" << std::endl;
572  std::cout << "alpha*BKG_mode + beta*SIG_mode = DATA_mode" << std::endl;
573  std::cout << "==========================================" << std::endl;
574  std::cout << "FSProton:" << std::endl;
575  std::cout << " alpha = " << alphaFSProtonLoCVN << std::endl;
576  std::cout << " beta = " << betaFSProtonLoCVN <<std::endl;
577  std::cout << " Cov Matrix: | " << C11FSProtonLoCVN << " " << C12FSProtonLoCVN << " |" << std::endl;
578  std::cout << " | " << C21FSProtonLoCVN << " " << C22FSProtonLoCVN << " |" << std::endl;
579  std::cout << " sqrt(Eigvl): " << majorFSProtonLoCVN << std::endl;
580  std::cout << " " << minorFSProtonLoCVN << std::endl;
581  std::cout << " wsfrc = " << wsfracFSProtonLoCVN << " +/- " << wsfracFSProtonLoCVNErr << std::endl;
582  std::cout << " b/a = " << betaAlphaFSProtonLoCVN << " +/- " << betaAlphaFSProtonLoCVNErr << std::endl;
583  std::cout << std::endl;
584  std::cout << "CVNProng:" << std::endl;
585  std::cout << " alpha = " << alphaCVNProngLoCVN << std::endl;
586  std::cout << " beta = " << betaCVNProngLoCVN << std::endl;
587  std::cout << " Cov Matrix: | " << C11CVNProngLoCVN << " " << C12CVNProngLoCVN << " |" << std::endl;
588  std::cout << " | " << C21CVNProngLoCVN << " " << C22CVNProngLoCVN << " |" << std::endl;
589  std::cout << " sqrt(Eigvl): " << majorCVNProngLoCVN << std::endl;
590  std::cout << " " << minorCVNProngLoCVN << std::endl;
591  std::cout << " wsfrc = " << wsfracCVNProngLoCVN << " +/- " << wsfracCVNProngLoCVNErr << std::endl;
592  std::cout << " b/a = " << betaAlphaCVNProngLoCVN << " +/- " << betaAlphaCVNProngLoCVNErr << std::endl;
593  std::cout << std::endl;
594  std::cout << "BDT3:" << std::endl;
595  std::cout << " alpha = " << alphaBDT3LoCVN << std::endl;
596  std::cout << " beta = " << betaBDT3LoCVN << std::endl;
597  std::cout << " Cov Matrix: | " << C11BDT3LoCVN << " " << C12BDT3LoCVN << " |" << std::endl;
598  std::cout << " | " << C21BDT3LoCVN << " " << C22BDT3LoCVN << " |" << std::endl;
599  std::cout << " sqrt(Eigvl): " << majorBDT3LoCVN << std::endl;
600  std::cout << " " << minorBDT3LoCVN << std::endl;
601  std::cout << " wsfrc = " << wsfracBDT3LoCVN << " +/- " << wsfracBDT3LoCVNErr << std::endl;
602  std::cout << " b/a = " << betaAlphaBDT3LoCVN << " +/- " << betaAlphaBDT3LoCVNErr << std::endl;
603  std::cout << "==========================================" << std::endl;
604  std::cout << std::endl;
605  std::cout << std::endl;
606  }
607  if( hasData ){
608  // High CVN
609  // FSProton
610  float numeratorFSProtonHiCVN = NueDataCt[0][0][0] - ((NueMCWS[0][0][0]/NueMCWS[0][0][1])*NueDataCt[0][0][1]);
611  float denominatorFSProtonHiCVN = NueMCOther[0][0][0] - ((NueMCWS[0][0][0]/NueMCWS[0][0][1])*NueMCOther[0][0][1]);
612  float alphaFSProtonHiCVN = numeratorFSProtonHiCVN / denominatorFSProtonHiCVN;
613  float betaFSProtonHiCVN = (NueDataCt[0][0][1] - (alphaFSProtonHiCVN*NueMCOther[0][0][1]))/NueMCWS[0][0][1];
614  float denomCovarFSProtonHiCVN = (NueMCOther[0][0][0]*NueMCWS[0][0][1]) - (NueMCOther[0][0][1]*NueMCWS[0][0][0]);
615  float A11FSProtonHiCVN = NueMCWS[0][0][1] / denomCovarFSProtonHiCVN;
616  float A12FSProtonHiCVN = -NueMCWS[0][0][0] / denomCovarFSProtonHiCVN;
617  float A21FSProtonHiCVN = -NueMCOther[0][0][1] / denomCovarFSProtonHiCVN;
618  float A22FSProtonHiCVN = NueMCOther[0][0][0] / denomCovarFSProtonHiCVN;
619  float C11FSProtonHiCVN = (A11FSProtonHiCVN*A11FSProtonHiCVN*NueDataCt[0][0][0]) + (A12FSProtonHiCVN*A12FSProtonHiCVN*NueDataCt[0][0][1]);
620  float C12FSProtonHiCVN = (A11FSProtonHiCVN*A21FSProtonHiCVN*NueDataCt[0][0][0]) + (A12FSProtonHiCVN*A22FSProtonHiCVN*NueDataCt[0][0][1]);
621  float C21FSProtonHiCVN = C12FSProtonHiCVN;
622  float C22FSProtonHiCVN = (A21FSProtonHiCVN*A21FSProtonHiCVN*NueDataCt[0][0][0]) + (A22FSProtonHiCVN*A22FSProtonHiCVN*NueDataCt[0][0][1]);
623  float majorFSProtonHiCVN = sqrt(0.5*(C11FSProtonHiCVN + C22FSProtonHiCVN +
624  sqrt( (C11FSProtonHiCVN+C22FSProtonHiCVN)*(C11FSProtonHiCVN+C22FSProtonHiCVN)
625  -(4*( C11FSProtonHiCVN*C22FSProtonHiCVN - C12FSProtonHiCVN*C21FSProtonHiCVN )))));
626  float minorFSProtonHiCVN = sqrt(0.5*(C11FSProtonHiCVN + C22FSProtonHiCVN -
627  sqrt( (C11FSProtonHiCVN+C22FSProtonHiCVN)*(C11FSProtonHiCVN+C22FSProtonHiCVN)
628  -(4*( C11FSProtonHiCVN*C22FSProtonHiCVN - C12FSProtonHiCVN*C21FSProtonHiCVN )))));
629  float rotBetaXFSProtonHiCVN = 0.5*(180/3.14159)*atan2(2.0*C12FSProtonHiCVN,C22FSProtonHiCVN - C11FSProtonHiCVN);
630  float wsfracFSProtonHiCVN =
631  (betaFSProtonHiCVN*(NueMCWS[0][0][0]+NueMCWS[0][0][1])) /
632  (betaFSProtonHiCVN*(NueMCWS[0][0][0]+NueMCWS[0][0][1]) + alphaFSProtonHiCVN*(NueMCOther[0][0][0]+NueMCOther[0][0][1]));
633  float wsfracFSProtonHiCVNErrDenom =
634  (betaFSProtonHiCVN*(NueMCWS[0][0][0]+NueMCWS[0][0][1]) + alphaFSProtonHiCVN*(NueMCOther[0][0][0]+NueMCOther[0][0][1])) *
635  (betaFSProtonHiCVN*(NueMCWS[0][0][0]+NueMCWS[0][0][1]) + alphaFSProtonHiCVN*(NueMCOther[0][0][0]+NueMCOther[0][0][1])) *
636  (betaFSProtonHiCVN*(NueMCWS[0][0][0]+NueMCWS[0][0][1]) + alphaFSProtonHiCVN*(NueMCOther[0][0][0]+NueMCOther[0][0][1])) *
637  (betaFSProtonHiCVN*(NueMCWS[0][0][0]+NueMCWS[0][0][1]) + alphaFSProtonHiCVN*(NueMCOther[0][0][0]+NueMCOther[0][0][1])); // term^4
638  float wsfracFSProtonHiCVNErrAB2 =
639  ( (NueMCWS[0][0][0]+NueMCWS[0][0][1])*(NueMCOther[0][0][0]+NueMCOther[0][0][1]) ) *
640  ( (NueMCWS[0][0][0]+NueMCWS[0][0][1])*(NueMCOther[0][0][0]+NueMCOther[0][0][1]) ); //AB^2
641  float wsfracFSProtonHiCVNErr =
642  sqrt( C11FSProtonHiCVN*wsfracFSProtonHiCVNErrAB2*betaFSProtonHiCVN*betaFSProtonHiCVN/wsfracFSProtonHiCVNErrDenom +
643  C22FSProtonHiCVN*wsfracFSProtonHiCVNErrAB2*alphaFSProtonHiCVN*alphaFSProtonHiCVN/wsfracFSProtonHiCVNErrDenom -
644  2*C12FSProtonHiCVN*wsfracFSProtonHiCVNErrAB2*alphaFSProtonHiCVN*betaFSProtonHiCVN/wsfracFSProtonHiCVNErrDenom );
645  float betaAlphaFSProtonHiCVN = betaFSProtonHiCVN / alphaFSProtonHiCVN;
646  float betaAlphaFSProtonHiCVNErr = sqrt(
647  C11FSProtonHiCVN*betaFSProtonHiCVN*betaFSProtonHiCVN/(alphaFSProtonHiCVN*alphaFSProtonHiCVN*alphaFSProtonHiCVN*alphaFSProtonHiCVN)+
648  C22FSProtonHiCVN*(1./alphaFSProtonHiCVN)*(1./alphaFSProtonHiCVN)-
649  2*C12FSProtonHiCVN*betaFSProtonHiCVN/(alphaFSProtonHiCVN*alphaFSProtonHiCVN*alphaFSProtonHiCVN) );
650  // CVNProng
651  float numeratorCVNProngHiCVN = NueDataCt[1][0][0] - ((NueMCWS[1][0][0]/NueMCWS[1][0][1])*NueDataCt[1][0][1]);
652  float denominatorCVNProngHiCVN = NueMCOther[1][0][0] - ((NueMCWS[1][0][0]/NueMCWS[1][0][1])*NueMCOther[1][0][1]);
653  float alphaCVNProngHiCVN = numeratorCVNProngHiCVN / denominatorCVNProngHiCVN;
654  float betaCVNProngHiCVN = (NueDataCt[1][0][1] - (alphaCVNProngHiCVN*NueMCOther[1][0][1]))/NueMCWS[1][0][1];
655  float denomCovarCVNProngHiCVN = (NueMCOther[1][0][0]*NueMCWS[1][0][1]) - (NueMCOther[1][0][1]*NueMCWS[1][0][0]);
656  float A11CVNProngHiCVN = NueMCWS[1][0][1] / denomCovarCVNProngHiCVN;
657  float A12CVNProngHiCVN = -NueMCWS[1][0][0] / denomCovarCVNProngHiCVN;
658  float A21CVNProngHiCVN = -NueMCOther[1][0][1] / denomCovarCVNProngHiCVN;
659  float A22CVNProngHiCVN = NueMCOther[1][0][0] / denomCovarCVNProngHiCVN;
660  float C11CVNProngHiCVN = (A11CVNProngHiCVN*A11CVNProngHiCVN*NueDataCt[1][0][0]) + (A12CVNProngHiCVN*A12CVNProngHiCVN*NueDataCt[1][0][1]);
661  float C12CVNProngHiCVN = (A11CVNProngHiCVN*A21CVNProngHiCVN*NueDataCt[1][0][0]) + (A12CVNProngHiCVN*A22CVNProngHiCVN*NueDataCt[1][0][1]);
662  float C21CVNProngHiCVN = C12CVNProngHiCVN;
663  float C22CVNProngHiCVN = (A21CVNProngHiCVN*A21CVNProngHiCVN*NueDataCt[1][0][0]) + (A22CVNProngHiCVN*A22CVNProngHiCVN*NueDataCt[1][0][1]);
664  float majorCVNProngHiCVN = sqrt(0.5*(C11CVNProngHiCVN + C22CVNProngHiCVN +
665  sqrt( (C11CVNProngHiCVN+C22CVNProngHiCVN)*(C11CVNProngHiCVN+C22CVNProngHiCVN)
666  -(4*( C11CVNProngHiCVN*C22CVNProngHiCVN - C12CVNProngHiCVN*C21CVNProngHiCVN )))));
667  float minorCVNProngHiCVN = sqrt(0.5*(C11CVNProngHiCVN + C22CVNProngHiCVN -
668  sqrt( (C11CVNProngHiCVN+C22CVNProngHiCVN)*(C11CVNProngHiCVN+C22CVNProngHiCVN)
669  -(4*( C11CVNProngHiCVN*C22CVNProngHiCVN - C12CVNProngHiCVN*C21CVNProngHiCVN )))));
670  float rotBetaXCVNProngHiCVN = 0.5*(180/3.14159)*atan2(2.0*C12CVNProngHiCVN,C22CVNProngHiCVN - C11CVNProngHiCVN);
671  float wsfracCVNProngHiCVN =
672  (betaCVNProngHiCVN*(NueMCWS[1][0][0]+NueMCWS[1][0][1])) /
673  (betaCVNProngHiCVN*(NueMCWS[1][0][0]+NueMCWS[1][0][1]) + alphaCVNProngHiCVN*(NueMCOther[1][0][0]+NueMCOther[1][0][1]));
674  float wsfracCVNProngHiCVNErrDenom =
675  (betaCVNProngHiCVN*(NueMCWS[1][0][0]+NueMCWS[1][0][1]) + alphaCVNProngHiCVN*(NueMCOther[1][0][0]+NueMCOther[1][0][1])) *
676  (betaCVNProngHiCVN*(NueMCWS[1][0][0]+NueMCWS[1][0][1]) + alphaCVNProngHiCVN*(NueMCOther[1][0][0]+NueMCOther[1][0][1])) *
677  (betaCVNProngHiCVN*(NueMCWS[1][0][0]+NueMCWS[1][0][1]) + alphaCVNProngHiCVN*(NueMCOther[1][0][0]+NueMCOther[1][0][1])) *
678  (betaCVNProngHiCVN*(NueMCWS[1][0][0]+NueMCWS[1][0][1]) + alphaCVNProngHiCVN*(NueMCOther[1][0][0]+NueMCOther[1][0][1])); // term^4
679  float wsfracCVNProngHiCVNErrAB2 =
680  ( (NueMCWS[1][0][0]+NueMCWS[1][0][1])*(NueMCOther[1][0][0]+NueMCOther[1][0][1]) ) *
681  ( (NueMCWS[1][0][0]+NueMCWS[1][0][1])*(NueMCOther[1][0][0]+NueMCOther[1][0][1]) ); //AB^2
682  float wsfracCVNProngHiCVNErr =
683  sqrt( C11CVNProngHiCVN*wsfracCVNProngHiCVNErrAB2*betaCVNProngHiCVN*betaCVNProngHiCVN/wsfracCVNProngHiCVNErrDenom +
684  C22CVNProngHiCVN*wsfracCVNProngHiCVNErrAB2*alphaCVNProngHiCVN*alphaCVNProngHiCVN/wsfracCVNProngHiCVNErrDenom -
685  2*C12CVNProngHiCVN*wsfracCVNProngHiCVNErrAB2*alphaCVNProngHiCVN*betaCVNProngHiCVN/wsfracCVNProngHiCVNErrDenom );
686  float betaAlphaCVNProngHiCVN = betaCVNProngHiCVN / alphaCVNProngHiCVN;
687  float betaAlphaCVNProngHiCVNErr = sqrt(
688  C11CVNProngHiCVN*betaCVNProngHiCVN*betaCVNProngHiCVN/(alphaCVNProngHiCVN*alphaCVNProngHiCVN*alphaCVNProngHiCVN*alphaCVNProngHiCVN)+
689  C22CVNProngHiCVN*(1./alphaCVNProngHiCVN)*(1./alphaCVNProngHiCVN)-
690  2*C12CVNProngHiCVN*betaCVNProngHiCVN/(alphaCVNProngHiCVN*alphaCVNProngHiCVN*alphaCVNProngHiCVN) );
691  // BDT3
692  float numeratorBDT3HiCVN = NueDataCt[2][0][0] - ((NueMCWS[2][0][0]/NueMCWS[2][0][1])*NueDataCt[2][0][1]);
693  float denominatorBDT3HiCVN = NueMCOther[2][0][0] - ((NueMCWS[2][0][0]/NueMCWS[2][0][1])*NueMCOther[2][0][1]);
694  float alphaBDT3HiCVN = numeratorBDT3HiCVN / denominatorBDT3HiCVN;
695  float betaBDT3HiCVN = (NueDataCt[2][0][1] - (alphaBDT3HiCVN*NueMCOther[2][0][1]))/NueMCWS[2][0][1];
696  float denomCovarBDT3HiCVN = (NueMCOther[2][0][0]*NueMCWS[2][0][1]) - (NueMCOther[2][0][1]*NueMCWS[2][0][0]);
697  float A11BDT3HiCVN = NueMCWS[2][0][1] / denomCovarBDT3HiCVN;
698  float A12BDT3HiCVN = -NueMCWS[2][0][0] / denomCovarBDT3HiCVN;
699  float A21BDT3HiCVN = -NueMCOther[2][0][1] / denomCovarBDT3HiCVN;
700  float A22BDT3HiCVN = NueMCOther[2][0][0] / denomCovarBDT3HiCVN;
701  float C11BDT3HiCVN = (A11BDT3HiCVN*A11BDT3HiCVN*NueDataCt[2][0][0]) + (A12BDT3HiCVN*A12BDT3HiCVN*NueDataCt[2][0][1]);
702  float C12BDT3HiCVN = (A11BDT3HiCVN*A21BDT3HiCVN*NueDataCt[2][0][0]) + (A12BDT3HiCVN*A22BDT3HiCVN*NueDataCt[2][0][1]);
703  float C21BDT3HiCVN = C12BDT3HiCVN;
704  float C22BDT3HiCVN = (A21BDT3HiCVN*A21BDT3HiCVN*NueDataCt[2][0][0]) + (A22BDT3HiCVN*A22BDT3HiCVN*NueDataCt[2][0][1]);
705  float majorBDT3HiCVN = sqrt(0.5*(C11BDT3HiCVN + C22BDT3HiCVN +
706  sqrt( (C11BDT3HiCVN+C22BDT3HiCVN)*(C11BDT3HiCVN+C22BDT3HiCVN)
707  -(4*( C11BDT3HiCVN*C22BDT3HiCVN - C12BDT3HiCVN*C21BDT3HiCVN )))));
708  float minorBDT3HiCVN = sqrt(0.5*(C11BDT3HiCVN + C22BDT3HiCVN -
709  sqrt( (C11BDT3HiCVN+C22BDT3HiCVN)*(C11BDT3HiCVN+C22BDT3HiCVN)
710  -(4*( C11BDT3HiCVN*C22BDT3HiCVN - C12BDT3HiCVN*C21BDT3HiCVN )))));
711  float rotBetaXBDT3HiCVN = 0.5*(180/3.14159)*atan2(2.0*C12BDT3HiCVN,C22BDT3HiCVN - C11BDT3HiCVN);
712  float wsfracBDT3HiCVN =
713  (betaBDT3HiCVN*(NueMCWS[2][0][0]+NueMCWS[2][0][1])) /
714  (betaBDT3HiCVN*(NueMCWS[2][0][0]+NueMCWS[2][0][1]) + alphaBDT3HiCVN*(NueMCOther[2][0][0]+NueMCOther[2][0][1]));
715  float wsfracBDT3HiCVNErrDenom =
716  (betaBDT3HiCVN*(NueMCWS[2][0][0]+NueMCWS[2][0][1]) + alphaBDT3HiCVN*(NueMCOther[2][0][0]+NueMCOther[2][0][1])) *
717  (betaBDT3HiCVN*(NueMCWS[2][0][0]+NueMCWS[2][0][1]) + alphaBDT3HiCVN*(NueMCOther[2][0][0]+NueMCOther[2][0][1])) *
718  (betaBDT3HiCVN*(NueMCWS[2][0][0]+NueMCWS[2][0][1]) + alphaBDT3HiCVN*(NueMCOther[2][0][0]+NueMCOther[2][0][1])) *
719  (betaBDT3HiCVN*(NueMCWS[2][0][0]+NueMCWS[2][0][1]) + alphaBDT3HiCVN*(NueMCOther[2][0][0]+NueMCOther[2][0][1])); // term^4
720  float wsfracBDT3HiCVNErrAB2 =
721  ( (NueMCWS[2][0][0]+NueMCWS[2][0][1])*(NueMCOther[2][0][0]+NueMCOther[2][0][1]) ) *
722  ( (NueMCWS[2][0][0]+NueMCWS[2][0][1])*(NueMCOther[2][0][0]+NueMCOther[2][0][1]) ); //AB^2
723  float wsfracBDT3HiCVNErr =
724  sqrt( C11BDT3HiCVN*wsfracBDT3HiCVNErrAB2*betaBDT3HiCVN*betaBDT3HiCVN/wsfracBDT3HiCVNErrDenom +
725  C22BDT3HiCVN*wsfracBDT3HiCVNErrAB2*alphaBDT3HiCVN*alphaBDT3HiCVN/wsfracBDT3HiCVNErrDenom -
726  2*C12BDT3HiCVN*wsfracBDT3HiCVNErrAB2*alphaBDT3HiCVN*betaBDT3HiCVN/wsfracBDT3HiCVNErrDenom );
727  float betaAlphaBDT3HiCVN = betaBDT3HiCVN / alphaBDT3HiCVN;
728  float betaAlphaBDT3HiCVNErr = sqrt(
729  C11BDT3HiCVN*betaBDT3HiCVN*betaBDT3HiCVN/(alphaBDT3HiCVN*alphaBDT3HiCVN*alphaBDT3HiCVN*alphaBDT3HiCVN)+
730  C22BDT3HiCVN*(1./alphaBDT3HiCVN)*(1./alphaBDT3HiCVN)-
731  2*C12BDT3HiCVN*betaBDT3HiCVN/(alphaBDT3HiCVN*alphaBDT3HiCVN*alphaBDT3HiCVN) );
732  std::cout << "==========================================" << std::endl;
733  std::cout << " SYSTEM OF EQUATIONS RESULTS: NUE " << std::endl;
734  std::cout << " HI CVN BIN (CVNSSe >=0.98) " << std::endl;
735  std::cout << "------------------------------------------" << std::endl;
736  std::cout << "alpha*BKG_mode + beta*SIG_mode = DATA_mode" << std::endl;
737  std::cout << "==========================================" << std::endl;
738  std::cout << "FSProton:" << std::endl;
739  std::cout << " alpha = " << alphaFSProtonHiCVN << std::endl;
740  std::cout << " beta = " << betaFSProtonHiCVN << std::endl;
741  std::cout << " Cov Matrix: | " << C11FSProtonHiCVN << " " << C12FSProtonHiCVN << " |" << std::endl;
742  std::cout << " | " << C21FSProtonHiCVN << " " << C22FSProtonHiCVN << " |" << std::endl;
743  std::cout << " sqrt(Eigvl): " << majorFSProtonHiCVN << std::endl;
744  std::cout << " " << minorFSProtonHiCVN << std::endl;
745  std::cout << " wsfrac = " << wsfracFSProtonHiCVN << " +/- " << wsfracFSProtonHiCVNErr << std::endl;
746  std::cout << " b/a = " << betaAlphaFSProtonHiCVN << " +/- " << betaAlphaFSProtonHiCVNErr << std::endl;
747  std::cout << std::endl;
748  std::cout << "CVNProng:" << std::endl;
749  std::cout << " alpha = " << alphaCVNProngHiCVN << std::endl;
750  std::cout << " beta = " << betaCVNProngHiCVN << std::endl;
751  std::cout << " Cov Matrix: | " << C11CVNProngHiCVN << " " << C12CVNProngHiCVN << " |" << std::endl;
752  std::cout << " | " << C21CVNProngHiCVN << " " << C22CVNProngHiCVN << " |" << std::endl;
753  std::cout << " sqrt(Eigvl): " << majorCVNProngHiCVN << std::endl;
754  std::cout << " " << minorCVNProngHiCVN << std::endl;
755  std::cout << " wsfrac = " << wsfracCVNProngHiCVN << " +/- " << wsfracCVNProngHiCVNErr << std::endl;
756  std::cout << " b/a = " << betaAlphaCVNProngHiCVN << " +/- " << betaAlphaCVNProngHiCVNErr << std::endl;
757  std::cout << std::endl;
758  std::cout << "BDT3:" << std::endl;
759  std::cout << " alpha = " << alphaBDT3HiCVN << std::endl;
760  std::cout << " beta = " << betaBDT3HiCVN << std::endl;
761  std::cout << " Cov Matrix: | " << C11BDT3HiCVN << " " << C12BDT3HiCVN << " |" << std::endl;
762  std::cout << " | " << C21BDT3HiCVN << " " << C22BDT3HiCVN << " |" << std::endl;
763  std::cout << " sqrt(Eigvl): " << majorBDT3HiCVN << std::endl;
764  std::cout << " " << minorBDT3HiCVN << std::endl;
765  std::cout << " wsfrac = " << wsfracBDT3HiCVN << " +/- " << wsfracBDT3HiCVNErr << std::endl;
766  std::cout << " b/a = " << betaAlphaBDT3HiCVN << " +/- " << betaAlphaBDT3HiCVNErr << std::endl;
767  std::cout << std::endl;
768  std::cout << "==========================================" << std::endl;
769  std::cout << std::endl;
770  std::cout << std::endl;
771  }
772 
773  std::cout << "Done." << std::endl;
774 
775  outFile->Close();
776 }
void Simulation()
Definition: tools.h:16
T max(const caf::Proxy< T > &a, T b)
enum BeamMode kOrange
Pass neutrinos through unchanged.
Definition: IOscCalc.h:45
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Antineutrinos-only.
Definition: IPrediction.h:50
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
T sqrt(T number)
Definition: d0nt_math.hpp:156
void calculateWrongSignNue()
TCanvas * calcAlphaBetaEachBin(Spectrum *templateSpec, std::string anaString, TH1D *wsS[], TH1D *wsB[], TH1D *wsD, TH1D *rsS[], TH1D *rsB[], TH1D *rsD, bool hasSplit=false, bool treatSplit=false)
static std::unique_ptr< PredictionNoOsc > LoadFrom(TDirectory *dir, const std::string &name)
double maxY
Definition: plot_hist.C:10
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
Double_t beta
osc::OscCalcDumb calc
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
ifstream inFile
Definition: AnaPlotMaker.h:34
Charged-current interactions.
Definition: IPrediction.h:39
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
TFile * outFile
Definition: PlotXSec.C:135
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
#define pot
double POT() const
Definition: Spectrum.h:227
void Preliminary()
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
All neutrinos, any flavor.
Definition: IPrediction.h:26
enum BeamMode kGreen
Prediction that wraps a simple Spectrum.
T atan2(T number)
Definition: d0nt_math.hpp:72
enum BeamMode string