calculateWrongSignNumuQ3.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////
2 // //
3 // WRONG SIGN NORMALIZATION STUDY SELECTION //
4 // NUMU Q3! //
5 // //
6 /////////////////////////////////////////////////////////////
7 
8 #include "CAFAna/Analysis/CountingExperimentStatistics.h"
9 #include "CAFAna/Analysis/Exposures.h"
10 #include "CAFAna/Fit/Fit.h"
11 #include "3FlavorAna/Plotting/NuePlotStyle.h"
12 #include "CAFAna/Analysis/Style.h"
13 #include "CAFAna/Vars/FitVars.h"
14 #include "CAFAna/Core/Spectrum.h"
15 #include "CAFAna/Core/LoadFromFile.h"
16 #include "CAFAna/Core/Utilities.h"
17 #include "CAFAna/Analysis/Calcs.h"
18 #include "CAFAna/Experiment/CountingExperiment.h"
19 #include "CAFAna/Experiment/SingleSampleExperiment.h"
20 #include "CAFAna/Experiment/MultiExperiment.h"
21 #include "CAFAna/Experiment/Dmsq32Constraint.h"
22 #include "CAFAna/Experiment/AtmConstraint.h"
23 #include "CAFAna/Experiment/ReactorExperiment.h"
24 #include "CAFAna/Core/Progress.h"
25 #include "CAFAna/Prediction/IPrediction.h"
26 #include "CAFAna/Prediction/PredictionNoExtrap.h"
27 #include "CAFAna/Prediction/PredictionInterp.h"
28 #include "CAFAna/Prediction/PredictionExtrap.h"
29 #include "CAFAna/Prediction/PredictionCombinePeriods.h"
30 #include "CAFAna/Systs/SystComponentScale.h"
31 #include "OscLib/OscCalcPMNSOpt.h"
32 #include "OscLib/OscCalcDumb.h"
33 #include "Utilities/func/MathUtil.h"
34 #include "CAFAna/Prediction/PredictionNoOsc.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 #include "TPaveText.h"
54 
55 #include <iostream>
56 #include <iomanip>
57 #include <string>
58 #include <vector>
59 #include "Utilities/rootlogon.C"
60 
61 TCanvas *calcAlphaBetaEachBin(Spectrum *templateSpec, std::string anaString, TH1D *wsS[], TH1D *wsB[], TH1D *wsD,
62  TH1D *rsS[], TH1D *rsB[], TH1D *rsD, bool hasSplit=false, bool treatSplit=false){
63  TH1D *alphaHist = templateSpec->ToTH1( 8.0e20 );
64  TH1D *betaHist = templateSpec->ToTH1( 8.0e20 );
65 
66  int nBinsHist = alphaHist->GetXaxis()->GetNbins();
67 
68  float bNumuMCWS[2];
69  float bNumuMCOther[2];
70  float bNumuDataCt[2];
71 
72  // Loop through bins and calculate alpha/beta
73  for( int iBinHist=1; iBinHist<=nBinsHist; ++iBinHist ){
74  bNumuMCWS[0] = 0.;
75  bNumuMCWS[1] = 0.;
76  bNumuMCOther[0] = 0.;
77  bNumuMCOther[1] = 0.;
78  bNumuDataCt[0] = 0.;
79  bNumuDataCt[1] = 0.;
80 
81  if( !hasSplit && !treatSplit ){
82  bNumuMCWS[0] = rsS[0]->GetBinContent(iBinHist) + rsS[1]->GetBinContent(iBinHist);
83  bNumuMCOther[0] = rsB[0]->GetBinContent(iBinHist) + rsB[1]->GetBinContent(iBinHist) + rsB[2]->GetBinContent(iBinHist);
84  }
85  else if( hasSplit && !treatSplit ){
86  bNumuMCWS[0] = rsS[0]->GetBinContent(iBinHist) + rsS[1]->GetBinContent(iBinHist);
87  bNumuMCOther[0] = rsB[0]->GetBinContent(iBinHist) + rsB[1]->GetBinContent(iBinHist) +
88  rsB[2]->GetBinContent(iBinHist) + rsS[2]->GetBinContent(iBinHist);
89  }
90  else if( hasSplit && treatSplit ){
91  bNumuMCWS[0] = rsS[0]->GetBinContent(iBinHist) + rsS[1]->GetBinContent(iBinHist) + rsS[2]->GetBinContent(iBinHist);
92  bNumuMCOther[0] = rsB[0]->GetBinContent(iBinHist) + rsB[1]->GetBinContent(iBinHist) + rsB[2]->GetBinContent(iBinHist);
93  }
94  else{
95  std::cout << "Doing something wrong in bin-by-bin method...";
96  bNumuMCWS[0] = 0.;
97  bNumuMCOther[0] = 0.;
98  }
99  bNumuDataCt[0] = rsD->GetBinContent(iBinHist);
100 
101  if( !hasSplit && !treatSplit ){
102  bNumuMCWS[1] = wsS[0]->GetBinContent(iBinHist) + wsS[1]->GetBinContent(iBinHist);
103  bNumuMCOther[1] = wsB[0]->GetBinContent(iBinHist) + wsB[1]->GetBinContent(iBinHist) + wsB[2]->GetBinContent(iBinHist);
104  }
105  else if( hasSplit && !treatSplit ){
106  bNumuMCWS[1] = wsS[0]->GetBinContent(iBinHist) + wsS[1]->GetBinContent(iBinHist);
107  bNumuMCOther[1] = wsB[0]->GetBinContent(iBinHist) + wsB[1]->GetBinContent(iBinHist) +
108  wsB[2]->GetBinContent(iBinHist) + wsS[2]->GetBinContent(iBinHist);
109  }
110  else if( hasSplit && treatSplit ){
111  bNumuMCWS[1] = wsS[0]->GetBinContent(iBinHist) + wsS[1]->GetBinContent(iBinHist) + wsS[2]->GetBinContent(iBinHist);
112  bNumuMCOther[1] = wsB[0]->GetBinContent(iBinHist) + wsB[1]->GetBinContent(iBinHist) + wsB[2]->GetBinContent(iBinHist);
113  }
114  else{
115  std::cout << "Doing something wrong in bin-by-bin method...";
116  bNumuMCWS[1] = 0.;
117  bNumuMCOther[1] = 0.;
118  }
119  bNumuDataCt[1] = wsD->GetBinContent(iBinHist);
120 
121  float numerator = bNumuDataCt[0] - ((bNumuMCWS[0]/bNumuMCWS[1])*bNumuDataCt[1]);
122  float denominator = bNumuMCOther[0] - ((bNumuMCWS[0]/bNumuMCWS[1])*bNumuMCOther[1]);
123  float alpha = numerator / denominator;
124  float beta = (bNumuDataCt[1] - (alpha*bNumuMCOther[1]))/bNumuMCWS[1];
125 
126  alphaHist->SetBinContent( iBinHist, alpha );
127  alphaHist->SetBinError( iBinHist, 0. );
128  betaHist->SetBinContent( iBinHist, beta );
129  betaHist->SetBinError( iBinHist, 0. );
130  }
131 
132  alphaHist->GetYaxis()->SetTitle("#alpha parameter");
133  betaHist->GetYaxis()->SetTitle("#beta parameter");
134 
135  // Print a canvas and return
136  TCanvas *tC = new TCanvas( TString::Format("binsAlphaBeta_%s",anaString.c_str()) );
137  tC->Divide(1,2);
138  tC->cd(1);
139  CenterTitles(alphaHist);
140  alphaHist->Draw();
141  Preliminary();
142  tC->cd(2);
143  CenterTitles(betaHist);
144  betaHist->Draw();
145  return tC;
146 }
147 
148 void calculateWrongSignNumuQ3()
149 {
150  // with Data File!!
151  bool hasData = true;
152  bool splitNCe[2] = {false,false}; // first is if you HAVE split NC, second is if you want to use it
153  bool splitNCmu[2] = {true,true}; // first is if you HAVE split NC, second is if you want to use it
154  bool eachBinMu = false;
155  bool quantileMC = false;
156 
157  if( (!splitNCe[0] && splitNCe[1]) || (!splitNCmu[0] && splitNCmu[1]) ){
158  std::cout << "NC split parameter mismatch. Quitting..." << std::endl;
159  return;
160  }
161 
162  std::string strOutFile = "./WrongSignNormalizationStudy_NumuQ3";
163  std::string strSplitNCe = "_plotNCsplitE";
164  std::string strSplitNCeMd = "_treatNCsplitE";
165  std::string strSplitNCmu = "_plotNCsplitMu";
166  std::string strSplitNCmuMd = "_treatNCsplitMu";
167  std::string strEachBinMu = "_withEachBinMu";
168  std::string strQuantileMC = "_QuantilesInMC";
169  if( splitNCe[0] && !splitNCe[1] ) strOutFile += strSplitNCe;
170  if( splitNCe[0] && splitNCe[1] ) strOutFile += strSplitNCeMd;
171  if( splitNCmu[0] && !splitNCmu[1] ) strOutFile += strSplitNCmu;
172  if( splitNCmu[0] && splitNCmu[1] ) strOutFile += strSplitNCmuMd;
173  if( eachBinMu ) strOutFile += strEachBinMu;
174  if( quantileMC ) strOutFile += strQuantileMC;
175  TFile *outFile = new TFile( TString::Format("%s.root",strOutFile.c_str()), "recreate" );
176 
177  TFile *inFile = new TFile("/nova/ana/users/howard/Nue2018Ana/development/Predictions_And_Data-WSBinning-NCSplit-NumuQ3OnlyDataStudy-31March2018/predictions_and_data_WrongsignRHC_splitNC_onlyNumuQ3.root","read");
178 
179  osc::NoOscillations calc;
180  std::vector< std::string > anaTitles = {"FSProton","CVNProng","BDT3"};
181  std::vector< std::string > nueModes = {"HiCVN","LoCVN"};
182 
183  // ----- Load files -----
184  PredictionNoOsc *pRHCNumuRS[3];
185  PredictionNoOsc *pRHCNumuWS[3];
186  Spectrum *sNCNumuRS[3];
187  Spectrum *sNCAntiNumuRS[3];
188  Spectrum *sNCNumuWS[3];
189  Spectrum *sNCAntiNumuWS[3];
190 
191  PredictionNoOsc *pQRHCNumuRS[4][3];
192  PredictionNoOsc *pQRHCNumuWS[4][3];
193 
194  Spectrum *sRHCNumuRSData[3];
195  Spectrum *sRHCNumuWSData[3];
196 
197  for( int iAna=0; iAna<3; ++iAna ){
198  pRHCNumuRS[iAna] = PredictionNoOsc::LoadFrom(inFile, TString::Format("pQ3Numu%sRS",anaTitles[iAna].c_str()) ).release();
199  pRHCNumuWS[iAna] = PredictionNoOsc::LoadFrom(inFile, TString::Format("pQ3Numu%sWS",anaTitles[iAna].c_str()) ).release();
200  if( splitNCmu[0] ){
201  sNCNumuRS[iAna] =
202  Spectrum::LoadFrom(inFile, TString::Format("sQ3Numu%sRSNC",anaTitles[iAna].c_str()) ).release();
203  sNCAntiNumuRS[iAna] =
204  Spectrum::LoadFrom(inFile, TString::Format("sQ3Numu%sRSNCAnti",anaTitles[iAna].c_str()) ).release();
205  sNCNumuWS[iAna] =
206  Spectrum::LoadFrom(inFile, TString::Format("sQ3Numu%sWSNC",anaTitles[iAna].c_str()) ).release();
207  sNCAntiNumuWS[iAna] =
208  Spectrum::LoadFrom(inFile, TString::Format("sQ3Numu%sWSNCAnti",anaTitles[iAna].c_str()) ).release();
209  }
210  if( hasData ){
211  sRHCNumuRSData[iAna] = Spectrum::LoadFrom(inFile, TString::Format("sNumuData_%sRS",anaTitles[iAna].c_str()) ).release();
212  sRHCNumuWSData[iAna] = Spectrum::LoadFrom(inFile, TString::Format("sNumuData_%sWS",anaTitles[iAna].c_str()) ).release();
213  }
214  }
215 
216  delete inFile;
217 
218  // Set POT, load POT from spectra if data
219  float pot = 8e20; // replace this with data POT when we have that...
220  if( hasData ){
221  pot = sRHCNumuRSData[0]->POT();
222  if( sRHCNumuRSData[0]->POT() != sRHCNumuWSData[0]->POT() ) std::cout << "You might want to check POT!!!!!!!!!!!" << std::endl;
223  }
224 
225  std::cout << "----------------------------" << std::endl;
226  std::cout << " POT: " << pot << std::endl;
227  std::cout << "----------------------------" << std::endl;
228 
229  // ----- Numu -----
230  TCanvas *cResultsNumu[3];
231  THStack *stackAnaBinsNumuRS[3];
232  THStack *stackAnaBinsNumuWS[3];
233  TLegend *numuLegend[3];
234  TH1D *numuRS_S[3][3]; // adding in 3rd category here for ability to split out NC
235  TH1D *numuRS_B[3][3];
236  TH1D *numuWS_S[3][3]; // adding in 3rd category here for ability to split out NC
237  TH1D *numuWS_B[3][3];
238  TH1D *numuRS_SCalc[3][3]; // adding in 3rd category here for ability to split out NC
239  TH1D *numuRS_BCalc[3][3];
240  TH1D *numuWS_SCalc[3][3]; // adding in 3rd category here for ability to split out NC
241  TH1D *numuWS_BCalc[3][3];
242  TH1D *numuRSDataResults[3];
243  TH1D *numuWSDataResults[3];
244 
245  // To be used later for error estimation
246  TH1D *numuErrorEstMC[3][2]; //RS, WS
247  TH1D *numuErrorEstData[3][2]; //RS, WS
248 
249  for( int iAna=0; iAna<3; ++iAna ){
250  // Get predictions
251  numuRS_S[iAna][0] = pRHCNumuRS[iAna]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUE
252  numuRS_S[iAna][0]->SetName( TString::Format("RSNue%s",anaTitles[iAna].c_str()) );
253  numuRS_S[iAna][1] = pRHCNumuRS[iAna]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUMU
254  numuRS_S[iAna][1]->SetName( TString::Format("RSNumu%s",anaTitles[iAna].c_str()) );
255  numuRS_B[iAna][0] = pRHCNumuRS[iAna]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUE
256  numuRS_B[iAna][0]->SetName( TString::Format("RSNuebar%s",anaTitles[iAna].c_str()) );
257  numuRS_B[iAna][1] = pRHCNumuRS[iAna]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUMU
258  numuRS_B[iAna][1]->SetName( TString::Format("RSNumubar%s",anaTitles[iAna].c_str()) );
259  if( !splitNCmu[0] ){
260  numuRS_B[iAna][2] = pRHCNumuRS[iAna]->PredictComponent(&calc,Flavors::kAll,Current::kNC,Sign::kBoth).ToTH1(pot); // NC
261  numuRS_B[iAna][2]->SetName( TString::Format("RSNc%s",anaTitles[iAna].c_str()) );
262  }
263  else{
264  numuRS_B[iAna][2] = sNCAntiNumuRS[iAna]->ToTH1(pot); // NCAnti
265  numuRS_B[iAna][2]->SetName( TString::Format("RSNcAnti%s",anaTitles[iAna].c_str()) );
266  numuRS_S[iAna][2] = sNCNumuRS[iAna]->ToTH1(pot); // NC
267  numuRS_S[iAna][2]->SetName( TString::Format("RSNc%s",anaTitles[iAna].c_str()) );
268  }
269  // Versions to use for calculations (since we add together the plot versions...)
270  numuRS_SCalc[iAna][0] = pRHCNumuRS[iAna]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUE
271  numuRS_SCalc[iAna][0]->SetName( TString::Format("RSCalcNue%s",anaTitles[iAna].c_str()) );
272  numuRS_SCalc[iAna][1] = pRHCNumuRS[iAna]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUMU
273  numuRS_SCalc[iAna][1]->SetName( TString::Format("RSCalcNumu%s",anaTitles[iAna].c_str()) );
274  numuRS_BCalc[iAna][0] = pRHCNumuRS[iAna]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUE
275  numuRS_BCalc[iAna][0]->SetName( TString::Format("RSCalcNuebar%s",anaTitles[iAna].c_str()) );
276  numuRS_BCalc[iAna][1] = pRHCNumuRS[iAna]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUMU
277  numuRS_BCalc[iAna][1]->SetName( TString::Format("RSCalcNumubar%s",anaTitles[iAna].c_str()) );
278  if( !splitNCmu[0] ){
279  numuRS_BCalc[iAna][2] = pRHCNumuRS[iAna]->PredictComponent(&calc,Flavors::kAll,Current::kNC,Sign::kBoth).ToTH1(pot); // NC
280  numuRS_BCalc[iAna][2]->SetName( TString::Format("RSCalcNc%s",anaTitles[iAna].c_str()) );
281  }
282  else{
283  numuRS_BCalc[iAna][2] = sNCAntiNumuRS[iAna]->ToTH1(pot); // NCAnti
284  numuRS_BCalc[iAna][2]->SetName( TString::Format("RSCalcNcAnti%s",anaTitles[iAna].c_str()) );
285  numuRS_SCalc[iAna][2] = sNCNumuRS[iAna]->ToTH1(pot); // NC
286  numuRS_SCalc[iAna][2]->SetName( TString::Format("RSCalcNc%s",anaTitles[iAna].c_str()) );
287  }
288  // FOR THE ERROR ESTIMATION
289  numuErrorEstMC[iAna][0] = pRHCNumuRS[iAna]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kNu).ToTH1(pot);
290  // Set Colors
291  numuRS_S[iAna][0]->SetLineColor(kGreen);
292  numuRS_S[iAna][0]->SetFillColor(kGreen);
293  numuRS_S[iAna][1]->SetLineColor(kOrange);
294  numuRS_S[iAna][1]->SetFillColor(kOrange);
295  numuRS_B[iAna][0]->SetLineColor(kGreen+2);
296  numuRS_B[iAna][0]->SetFillColor(kGreen+2);
297  numuRS_B[iAna][1]->SetLineColor(kRed+1);
298  numuRS_B[iAna][1]->SetFillColor(kRed+1);
299  numuRS_B[iAna][2]->SetLineColor(kMagenta);
300  numuRS_B[iAna][2]->SetFillColor(kMagenta);
301  if( splitNCmu[0] ){
302  numuRS_S[iAna][2]->SetLineColor(kMagenta);
303  numuRS_S[iAna][2]->SetFillColor(kMagenta);
304  numuRS_B[iAna][2]->SetLineColor(kGray+2);
305  numuRS_B[iAna][2]->SetFillColor(kGray+2);
306  }
307  // Fill up the sum histograms for numu "stack"
308  // To do this, simply add the things below to the one you want to show
309  // i.e. reverse drawing order, which is (bot->top): nc, nue RS, nue WS, numu WS, numu RS
310  if( splitNCmu[0] ) numuRS_B[iAna][2]->Add(numuRS_S[iAna][2]);
311  numuRS_B[iAna][0]->Add(numuRS_B[iAna][2]);
312  numuRS_S[iAna][0]->Add(numuRS_B[iAna][0]);
313  numuRS_S[iAna][1]->Add(numuRS_S[iAna][0]);
314  numuRS_B[iAna][1]->Add(numuRS_S[iAna][1]);
315  // Get Data!
316  if( hasData ) numuRSDataResults[iAna] = sRHCNumuRSData[iAna]->ToTH1( pot );
317  if( hasData ) numuErrorEstData[iAna][0] = sRHCNumuRSData[iAna]->ToTH1( pot );
318  }
319 
320  for( int iAna=0; iAna<3; ++iAna ){
321  // Get predictions
322  numuWS_S[iAna][0] = pRHCNumuWS[iAna]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUE
323  numuWS_S[iAna][0]->SetName( TString::Format("WSNue%s",anaTitles[iAna].c_str()) );
324  numuWS_S[iAna][1] = pRHCNumuWS[iAna]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUMU
325  numuWS_S[iAna][1]->SetName( TString::Format("WSNumu%s",anaTitles[iAna].c_str()) );
326  numuWS_B[iAna][0] = pRHCNumuWS[iAna]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUE
327  numuWS_B[iAna][0]->SetName( TString::Format("WSNuebar%s",anaTitles[iAna].c_str()) );
328  numuWS_B[iAna][1] = pRHCNumuWS[iAna]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUMU
329  numuWS_B[iAna][1]->SetName( TString::Format("WSNumubar%s",anaTitles[iAna].c_str()) );
330  if( !splitNCmu[0] ){
331  numuWS_B[iAna][2] = pRHCNumuWS[iAna]->PredictComponent(&calc,Flavors::kAll,Current::kNC,Sign::kBoth).ToTH1(pot); // NC
332  numuWS_B[iAna][2]->SetName( TString::Format("WSNc%s",anaTitles[iAna].c_str()) );
333  }
334  else{
335  numuWS_B[iAna][2] = sNCAntiNumuWS[iAna]->ToTH1(pot); // NCAnti
336  numuWS_B[iAna][2]->SetName( TString::Format("WSNcAnti%s",anaTitles[iAna].c_str()) );
337  numuWS_S[iAna][2] = sNCNumuWS[iAna]->ToTH1(pot); // NC
338  numuWS_S[iAna][2]->SetName( TString::Format("WSNc%s",anaTitles[iAna].c_str()) );
339  }
340  // Get predictions -- calc version as above
341  numuWS_SCalc[iAna][0] = pRHCNumuWS[iAna]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUE
342  numuWS_SCalc[iAna][0]->SetName( TString::Format("WSCalcNue%s",anaTitles[iAna].c_str()) );
343  numuWS_SCalc[iAna][1] = pRHCNumuWS[iAna]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kNu).ToTH1(pot); // WS NUMU
344  numuWS_SCalc[iAna][1]->SetName( TString::Format("WSCalcNumu%s",anaTitles[iAna].c_str()) );
345  numuWS_BCalc[iAna][0] = pRHCNumuWS[iAna]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUE
346  numuWS_BCalc[iAna][0]->SetName( TString::Format("WSCalcNuebar%s",anaTitles[iAna].c_str()) );
347  numuWS_BCalc[iAna][1] = pRHCNumuWS[iAna]->PredictComponent(&calc,Flavors::kAllNuMu,Current::kCC,Sign::kAntiNu).ToTH1(pot); // RS NUMU
348  numuWS_BCalc[iAna][1]->SetName( TString::Format("WSCalcNumubar%s",anaTitles[iAna].c_str()) );
349  if( !splitNCmu[0] ){
350  numuWS_BCalc[iAna][2] = pRHCNumuWS[iAna]->PredictComponent(&calc,Flavors::kAll,Current::kNC,Sign::kBoth).ToTH1(pot); // NC
351  numuWS_BCalc[iAna][2]->SetName( TString::Format("WSCalcNc%s",anaTitles[iAna].c_str()) );
352  }
353  else{
354  numuWS_BCalc[iAna][2] = sNCAntiNumuWS[iAna]->ToTH1(pot); // NCAnti
355  numuWS_BCalc[iAna][2]->SetName( TString::Format("WSCalcNcAnti%s",anaTitles[iAna].c_str()) );
356  numuWS_SCalc[iAna][2] = sNCNumuWS[iAna]->ToTH1(pot); // NC
357  numuWS_SCalc[iAna][2]->SetName( TString::Format("WSCalcNc%s",anaTitles[iAna].c_str()) );
358  }
359  // FOR ERROR ESTIMATION
360  numuErrorEstMC[iAna][1] = pRHCNumuWS[iAna]->PredictComponent(&calc,Flavors::kAllNuE,Current::kCC,Sign::kNu).ToTH1(pot);
361  // Set Colors
362  numuWS_S[iAna][0]->SetLineColor(kGreen);
363  numuWS_S[iAna][0]->SetFillColor(kGreen);
364  numuWS_S[iAna][1]->SetLineColor(kOrange);
365  numuWS_S[iAna][1]->SetFillColor(kOrange);
366  numuWS_B[iAna][0]->SetLineColor(kGreen+2);
367  numuWS_B[iAna][0]->SetFillColor(kGreen+2);
368  numuWS_B[iAna][1]->SetLineColor(kRed+1);
369  numuWS_B[iAna][1]->SetFillColor(kRed+1);
370  numuWS_B[iAna][2]->SetLineColor(kMagenta);
371  numuWS_B[iAna][2]->SetFillColor(kMagenta);
372  if( splitNCmu[0] ){
373  numuWS_S[iAna][2]->SetLineColor(kMagenta);
374  numuWS_S[iAna][2]->SetFillColor(kMagenta);
375  numuWS_B[iAna][2]->SetLineColor(kGray+2);
376  numuWS_B[iAna][2]->SetFillColor(kGray+2);
377  }
378  // Fill up the sum histograms for numu "stack"
379  // To do this, simply add the things below to the one you want to show
380  // i.e. reverse drawing order, which is (bot->top): nc, nue RS, nue WS, numu WS, numu RS
381  if( splitNCmu[0] ) numuWS_B[iAna][2]->Add(numuWS_S[iAna][2]);
382  numuWS_B[iAna][0]->Add(numuWS_B[iAna][2]);
383  numuWS_S[iAna][0]->Add(numuWS_B[iAna][0]);
384  numuWS_S[iAna][1]->Add(numuWS_S[iAna][0]);
385  numuWS_B[iAna][1]->Add(numuWS_S[iAna][1]);
386  // Get Data!
387  if( hasData ) numuWSDataResults[iAna] = sRHCNumuWSData[iAna]->ToTH1( pot );
388  if( hasData ) numuErrorEstData[iAna][1] = sRHCNumuWSData[iAna]->ToTH1( pot );
389  }
390 
391  double maxYmuRS = 1.2*std::max( std::max(numuRS_B[0][1]->GetMaximum(),numuRS_B[1][1]->GetMaximum()), numuRS_B[2][1]->GetMaximum());
392  double maxYmuWS = 1.2*std::max( std::max(numuWS_B[0][1]->GetMaximum(),numuWS_B[1][1]->GetMaximum()), numuWS_B[2][1]->GetMaximum());
393 
394  if( hasData ){
395  maxYmuRS = 1.2*std::max( std::max( std::max( std::max( std::max(
396  numuRS_B[0][1]->GetMaximum(),numuRS_B[1][1]->GetMaximum()),
397  numuRS_B[2][1]->GetMaximum()),
398  numuRSDataResults[0]->GetMaximum()),
399  numuRSDataResults[1]->GetMaximum()),
400  numuRSDataResults[2]->GetMaximum());
401  maxYmuWS = 1.2*std::max( std::max( std::max( std::max( std::max(
402  numuWS_B[0][1]->GetMaximum(),numuWS_B[1][1]->GetMaximum()),
403  numuWS_B[2][1]->GetMaximum()),
404  numuWSDataResults[0]->GetMaximum()),
405  numuWSDataResults[1]->GetMaximum()),
406  numuWSDataResults[2]->GetMaximum());
407  }
408 
409  double maxYmu = std::max(maxYmuRS,maxYmuWS);
410 
411  // dataType[method 0-2][WS suppressed (0), WS enhanced (1)]
412  float NumuMCWS[3][2];
413  float NumuMCOther[3][2];
414  float NumuDataCt[3][2];
415 
416  outFile->cd();
417 
418  for( int iAna=0; iAna<3; ++iAna ){
419  cResultsNumu[iAna] = new TCanvas();
420  cResultsNumu[iAna]->Divide(2,1);
421  cResultsNumu[iAna]->cd(1);
422  //numuRS_B[iAna][1]->GetYaxis()->SetRangeUser(0,maxYmu);
423  if( hasData )
424  numuRS_B[iAna][1]->GetYaxis()->SetRangeUser(0,1.2*std::max( numuRS_B[iAna][1]->GetMaximum(),numuRSDataResults[iAna]->GetMaximum() ));
425  numuRS_B[iAna][1]->GetYaxis()->SetLabelSize(0.03);
426  CenterTitles(numuRS_B[iAna][1]);
427  numuRS_B[iAna][1]->Draw("hist");
428  numuRS_S[iAna][1]->Draw("same hist");
429  numuRS_S[iAna][0]->Draw("same hist");
430  numuRS_B[iAna][0]->Draw("same hist");
431  numuRS_B[iAna][2]->Draw("same hist");
432  if( splitNCmu[0] ) numuRS_S[iAna][2]->Draw("same hist");
433  if( hasData ) numuRSDataResults[iAna]->Draw("same");
434  numuLegend[iAna] = new TLegend(0.7,0.63,0.88,0.88);
435  numuLegend[iAna]->AddEntry( TString::Format("RSNumu%s",anaTitles[iAna].c_str()), "#nu_{#mu}", "lf" );
436  numuLegend[iAna]->AddEntry( TString::Format("RSNue%s",anaTitles[iAna].c_str()), "#nu_{e}", "lf" );
437  numuLegend[iAna]->AddEntry( TString::Format("RSNumubar%s",anaTitles[iAna].c_str()), "#bar{#nu}_{#mu}", "lf" );
438  numuLegend[iAna]->AddEntry( TString::Format("RSNuebar%s",anaTitles[iAna].c_str()), "#bar{#nu}_{e}", "lf" );
439  numuLegend[iAna]->AddEntry( TString::Format("RSNc%s",anaTitles[iAna].c_str()), "NC", "lf" );
440  if( splitNCmu[0] ) numuLegend[iAna]->AddEntry( TString::Format("RSNcAnti%s",anaTitles[iAna].c_str()), "NCAnti", "lf" );
441  numuLegend[iAna]->Draw();
442  cResultsNumu[iAna]->cd(2);
443  //numuWS_B[iAna][1]->GetYaxis()->SetRangeUser(0,maxYmu);
444  if( hasData )
445  numuWS_B[iAna][1]->GetYaxis()->SetRangeUser(0,1.2*std::max( numuWS_B[iAna][1]->GetMaximum(),numuWSDataResults[iAna]->GetMaximum() ));
446  CenterTitles(numuWS_B[iAna][1]);
447  numuWS_B[iAna][1]->Draw("hist");
448  numuWS_S[iAna][1]->Draw("same hist");
449  numuWS_S[iAna][0]->Draw("same hist");
450  numuWS_B[iAna][0]->Draw("same hist");
451  numuWS_B[iAna][2]->Draw("same hist");
452  if( splitNCmu[0] ) numuWS_S[iAna][2]->Draw("same hist");
453  if( hasData ) numuWSDataResults[iAna]->Draw("same");
454  if( !hasData ) Simulation();
455  if( hasData ) Preliminary();
456  cResultsNumu[iAna]->Write( TString::Format("cRHCNumu%s",anaTitles[iAna].c_str()) );
457 
458  // Values
459  if( !splitNCmu[0] && !splitNCmu[1] ){
460  NumuMCWS[iAna][0] = numuRS_SCalc[iAna][0]->Integral() + numuRS_SCalc[iAna][1]->Integral();
461  NumuMCOther[iAna][0] = numuRS_BCalc[iAna][0]->Integral() + numuRS_BCalc[iAna][1]->Integral() + numuRS_BCalc[iAna][2]->Integral();
462  }
463  else if( splitNCmu[0] && !splitNCmu[1] ){
464  NumuMCWS[iAna][0] = numuRS_SCalc[iAna][0]->Integral() + numuRS_SCalc[iAna][1]->Integral();
465  NumuMCOther[iAna][0] = numuRS_BCalc[iAna][0]->Integral() + numuRS_BCalc[iAna][1]->Integral() +
466  numuRS_BCalc[iAna][2]->Integral() + numuRS_SCalc[iAna][2]->Integral();
467  }
468  else if( splitNCmu[0] && splitNCmu[1] ){
469  NumuMCWS[iAna][0] = numuRS_SCalc[iAna][0]->Integral() + numuRS_SCalc[iAna][1]->Integral() + numuRS_SCalc[iAna][2]->Integral();
470  NumuMCOther[iAna][0] = numuRS_BCalc[iAna][0]->Integral() + numuRS_BCalc[iAna][1]->Integral() + numuRS_BCalc[iAna][2]->Integral();
471  }
472  else{
473  std::cout << "IMPROPR NUMU SPLITTING SELECTED. FILLING WITH 0s" << std::endl;
474  NumuMCWS[iAna][0] = 0.;
475  NumuMCOther[iAna][0] = 0.;
476  }
477  if( hasData ) NumuDataCt[iAna][0] = numuRSDataResults[iAna]->Integral();
478 
479  if( !splitNCmu[0] && !splitNCmu[1] ){
480  NumuMCWS[iAna][1] = numuWS_SCalc[iAna][0]->Integral() + numuWS_SCalc[iAna][1]->Integral();
481  NumuMCOther[iAna][1] = numuWS_BCalc[iAna][0]->Integral() + numuWS_BCalc[iAna][1]->Integral() + numuWS_BCalc[iAna][2]->Integral();
482  }
483  else if( splitNCmu[0] && !splitNCmu[1] ){
484  NumuMCWS[iAna][1] = numuWS_SCalc[iAna][0]->Integral() + numuWS_SCalc[iAna][1]->Integral();
485  NumuMCOther[iAna][1] = numuWS_BCalc[iAna][0]->Integral() + numuWS_BCalc[iAna][1]->Integral() +
486  numuWS_BCalc[iAna][2]->Integral() + numuWS_SCalc[iAna][2]->Integral();
487  }
488  else if( splitNCmu[0] && splitNCmu[1] ){
489  NumuMCWS[iAna][1] = numuWS_SCalc[iAna][0]->Integral() + numuWS_SCalc[iAna][1]->Integral() + numuWS_SCalc[iAna][2]->Integral();
490  NumuMCOther[iAna][1] = numuWS_BCalc[iAna][0]->Integral() + numuWS_BCalc[iAna][1]->Integral() + numuWS_BCalc[iAna][2]->Integral();
491  }
492  else{
493  std::cout << "IMPROPR NUMU SPLITTING SELECTED. FILLING WITH 0s" << std::endl;
494  NumuMCWS[iAna][1] = 0.;
495  NumuMCOther[iAna][1] = 0.;
496  }
497  if( hasData ) NumuDataCt[iAna][1] = numuWSDataResults[iAna]->Integral();
498 
499  std::cout << "Numu " << anaTitles[iAna] << ":" << std::endl;
500  std::cout << " __Right-sign enhanced__" << std::endl;
501  std::cout << " signal (WS) = " << NumuMCWS[iAna][0] << std::endl;
502  std::cout << " bgd (RS+NC) = " << NumuMCOther[iAna][0] << std::endl;
503  if( hasData ) std::cout << " DATA = " << NumuDataCt[iAna][0] << std::endl;
504  std::cout << std::endl;
505  std::cout << " __Wrong-sign enhanced__" << std::endl;
506  std::cout << " signal (WS) = " << NumuMCWS[iAna][1] << std::endl;
507  std::cout << " bgd (RS+NC) = " << NumuMCOther[iAna][1] << std::endl;
508  if( hasData ) std::cout << " DATA = " << NumuDataCt[iAna][1] << std::endl;
509  if( splitNCmu[0] && splitNCmu[1] ){
510  // Print info on the right sign and wrong sign and composition of neutrino types
511  std::cout << std::endl;
512  std::cout << "WS: " << numuRS_SCalc[iAna][0]->Integral() +
513  numuRS_SCalc[iAna][1]->Integral() + numuWS_SCalc[iAna][0]->Integral() + numuWS_SCalc[iAna][1]->Integral() +
514  numuWS_SCalc[iAna][2]->Integral() + numuRS_SCalc[iAna][2]->Integral() << std::endl;
515  std::cout << "RS: " << numuRS_BCalc[iAna][0]->Integral() +
516  numuRS_BCalc[iAna][1]->Integral() + numuWS_BCalc[iAna][0]->Integral() + numuWS_BCalc[iAna][1]->Integral() +
517  numuWS_BCalc[iAna][2]->Integral() + numuRS_BCalc[iAna][2]->Integral() << std::endl;
518  std::cout << "Data: " << numuWSDataResults[iAna]->Integral() + numuRSDataResults[iAna]->Integral();
519  std::cout << std::endl;
520  std::cout << std::endl;
521  std::cout << "Nue (WS+RS): " << (numuRS_SCalc[iAna][0]->Integral() + numuRS_BCalc[iAna][0]->Integral() +
522  numuWS_SCalc[iAna][0]->Integral() + numuWS_BCalc[iAna][0]->Integral()) << std::endl;
523  std::cout << "Numu (WS+RS): " << (numuRS_SCalc[iAna][1]->Integral() + numuRS_BCalc[iAna][1]->Integral() +
524  numuWS_SCalc[iAna][1]->Integral() + numuWS_BCalc[iAna][1]->Integral()) << std::endl;
525  std::cout << "WS NC: " << numuWS_SCalc[iAna][2]->Integral() + numuRS_SCalc[iAna][2]->Integral() << std::endl;
526  std::cout << "RS NC: " << numuWS_BCalc[iAna][2]->Integral() + numuRS_BCalc[iAna][2]->Integral() << std::endl;
527  }
528  std::cout << std::endl;
529  std::cout << std::endl;
530  }
531 
532  float alphaNumu[] = {0., 0.};
533  float betaNumu[] = {0., 0.};
534 
535  // Solve numu system of equations, if hasData turned on
536  if( hasData ){
537  // High CVN
538  // FSProton
539  float numeratorFSProtonNumu = NumuDataCt[0][0] - ((NumuMCWS[0][0]/NumuMCWS[0][1])*NumuDataCt[0][1]);
540  float denominatorFSProtonNumu = NumuMCOther[0][0] - ((NumuMCWS[0][0]/NumuMCWS[0][1])*NumuMCOther[0][1]);
541  float alphaFSProtonNumu = numeratorFSProtonNumu / denominatorFSProtonNumu;
542  float betaFSProtonNumu = (NumuDataCt[0][1] - (alphaFSProtonNumu*NumuMCOther[0][1]))/NumuMCWS[0][1];
543  // FSProton error ellipses and WSFraction
544  float denomCovarFSProtonNumu = (NumuMCOther[0][0]*NumuMCWS[0][1]) - (NumuMCOther[0][1]*NumuMCWS[0][0]);
545  float A11FSProtonNumu = NumuMCWS[0][1] / denomCovarFSProtonNumu;
546  float A12FSProtonNumu = -NumuMCWS[0][0] / denomCovarFSProtonNumu;
547  float A21FSProtonNumu = -NumuMCOther[0][1] / denomCovarFSProtonNumu;
548  float A22FSProtonNumu = NumuMCOther[0][0] / denomCovarFSProtonNumu;
549  float C11FSProtonNumu = (A11FSProtonNumu*A11FSProtonNumu*NumuDataCt[0][0]) + (A12FSProtonNumu*A12FSProtonNumu*NumuDataCt[0][1]);
550  float C12FSProtonNumu = (A11FSProtonNumu*A21FSProtonNumu*NumuDataCt[0][0]) + (A12FSProtonNumu*A22FSProtonNumu*NumuDataCt[0][1]);
551  float C21FSProtonNumu = C12FSProtonNumu;
552  float C22FSProtonNumu = (A21FSProtonNumu*A21FSProtonNumu*NumuDataCt[0][0]) + (A22FSProtonNumu*A22FSProtonNumu*NumuDataCt[0][1]);
553  float majorFSProtonNumu = sqrt(0.5*(C11FSProtonNumu + C22FSProtonNumu +
554  sqrt( (C11FSProtonNumu+C22FSProtonNumu)*(C11FSProtonNumu+C22FSProtonNumu)
555  -(4*( C11FSProtonNumu*C22FSProtonNumu - C12FSProtonNumu*C21FSProtonNumu )))));
556  float minorFSProtonNumu = sqrt(0.5*(C11FSProtonNumu + C22FSProtonNumu -
557  sqrt( (C11FSProtonNumu+C22FSProtonNumu)*(C11FSProtonNumu+C22FSProtonNumu)
558  -(4*( C11FSProtonNumu*C22FSProtonNumu - C12FSProtonNumu*C21FSProtonNumu )))));
559  float rotBetaXFSProtonNumu = 0.5*(180./TMath::Pi())*atan2(2.0*C12FSProtonNumu,C22FSProtonNumu - C11FSProtonNumu);
560  float wsfracFSProtonNumu =
561  (betaFSProtonNumu*(NumuMCWS[0][0]+NumuMCWS[0][1])) /
562  (betaFSProtonNumu*(NumuMCWS[0][0]+NumuMCWS[0][1]) + alphaFSProtonNumu*(NumuMCOther[0][0]+NumuMCOther[0][1]));
563  float wsfracFSProtonNumuErrDenom =
564  (betaFSProtonNumu*(NumuMCWS[0][0]+NumuMCWS[0][1]) + alphaFSProtonNumu*(NumuMCOther[0][0]+NumuMCOther[0][1])) *
565  (betaFSProtonNumu*(NumuMCWS[0][0]+NumuMCWS[0][1]) + alphaFSProtonNumu*(NumuMCOther[0][0]+NumuMCOther[0][1])) *
566  (betaFSProtonNumu*(NumuMCWS[0][0]+NumuMCWS[0][1]) + alphaFSProtonNumu*(NumuMCOther[0][0]+NumuMCOther[0][1])) *
567  (betaFSProtonNumu*(NumuMCWS[0][0]+NumuMCWS[0][1]) + alphaFSProtonNumu*(NumuMCOther[0][0]+NumuMCOther[0][1])); // term^4
568  float wsfracFSProtonNumuErrAB2 =
569  ( (NumuMCWS[0][0]+NumuMCWS[0][1])*(NumuMCOther[0][0]+NumuMCOther[0][1]) ) *
570  ( (NumuMCWS[0][0]+NumuMCWS[0][1])*(NumuMCOther[0][0]+NumuMCOther[0][1]) ); //AB^2
571  float wsfracFSProtonNumuErr =
572  sqrt( C11FSProtonNumu*wsfracFSProtonNumuErrAB2*betaFSProtonNumu*betaFSProtonNumu/wsfracFSProtonNumuErrDenom +
573  C22FSProtonNumu*wsfracFSProtonNumuErrAB2*alphaFSProtonNumu*alphaFSProtonNumu/wsfracFSProtonNumuErrDenom -
574  2*C12FSProtonNumu*wsfracFSProtonNumuErrAB2*alphaFSProtonNumu*betaFSProtonNumu/wsfracFSProtonNumuErrDenom );
575  // CVNProng
576  float numeratorCVNProngNumu = NumuDataCt[1][0] - ((NumuMCWS[1][0]/NumuMCWS[1][1])*NumuDataCt[1][1]);
577  float denominatorCVNProngNumu = NumuMCOther[1][0] - ((NumuMCWS[1][0]/NumuMCWS[1][1])*NumuMCOther[1][1]);
578  float alphaCVNProngNumu = numeratorCVNProngNumu / denominatorCVNProngNumu;
579  float betaCVNProngNumu = (NumuDataCt[1][1] - (alphaCVNProngNumu*NumuMCOther[1][1]))/NumuMCWS[1][1];
580  // CVNProng error ellipses and WSFraction
581  float denomCovarCVNProngNumu = (NumuMCOther[1][0]*NumuMCWS[1][1]) - (NumuMCOther[1][1]*NumuMCWS[1][0]);
582  float A11CVNProngNumu = NumuMCWS[1][1] / denomCovarCVNProngNumu;
583  float A12CVNProngNumu = -NumuMCWS[1][0] / denomCovarCVNProngNumu;
584  float A21CVNProngNumu = -NumuMCOther[1][1] / denomCovarCVNProngNumu;
585  float A22CVNProngNumu = NumuMCOther[1][0] / denomCovarCVNProngNumu;
586  float C11CVNProngNumu = (A11CVNProngNumu*A11CVNProngNumu*NumuDataCt[1][0]) + (A12CVNProngNumu*A12CVNProngNumu*NumuDataCt[1][1]);
587  float C12CVNProngNumu = (A11CVNProngNumu*A21CVNProngNumu*NumuDataCt[1][0]) + (A12CVNProngNumu*A22CVNProngNumu*NumuDataCt[1][1]);
588  float C21CVNProngNumu = C12CVNProngNumu;
589  float C22CVNProngNumu = (A21CVNProngNumu*A21CVNProngNumu*NumuDataCt[1][0]) + (A22CVNProngNumu*A22CVNProngNumu*NumuDataCt[1][1]);
590  float majorCVNProngNumu = sqrt(0.5*(C11CVNProngNumu + C22CVNProngNumu +
591  sqrt( (C11CVNProngNumu+C22CVNProngNumu)*(C11CVNProngNumu+C22CVNProngNumu)
592  -(4*( C11CVNProngNumu*C22CVNProngNumu - C12CVNProngNumu*C21CVNProngNumu )))));
593  float minorCVNProngNumu = sqrt(0.5*(C11CVNProngNumu + C22CVNProngNumu -
594  sqrt( (C11CVNProngNumu+C22CVNProngNumu)*(C11CVNProngNumu+C22CVNProngNumu)
595  -(4*( C11CVNProngNumu*C22CVNProngNumu - C12CVNProngNumu*C21CVNProngNumu )))));
596  float rotBetaXCVNProngNumu = 0.5*(180./TMath::Pi())*atan2(2.0*C12CVNProngNumu,C22CVNProngNumu - C11CVNProngNumu);
597  float wsfracCVNProngNumu =
598  (betaCVNProngNumu*(NumuMCWS[1][0]+NumuMCWS[1][1])) /
599  (betaCVNProngNumu*(NumuMCWS[1][0]+NumuMCWS[1][1]) + alphaCVNProngNumu*(NumuMCOther[1][0]+NumuMCOther[1][1]));
600  float wsfracCVNProngNumuErrDenom =
601  (betaCVNProngNumu*(NumuMCWS[1][0]+NumuMCWS[1][1]) + alphaCVNProngNumu*(NumuMCOther[1][0]+NumuMCOther[1][1])) *
602  (betaCVNProngNumu*(NumuMCWS[1][0]+NumuMCWS[1][1]) + alphaCVNProngNumu*(NumuMCOther[1][0]+NumuMCOther[1][1])) *
603  (betaCVNProngNumu*(NumuMCWS[1][0]+NumuMCWS[1][1]) + alphaCVNProngNumu*(NumuMCOther[1][0]+NumuMCOther[1][1])) *
604  (betaCVNProngNumu*(NumuMCWS[1][0]+NumuMCWS[1][1]) + alphaCVNProngNumu*(NumuMCOther[1][0]+NumuMCOther[1][1])); // term^4
605  float wsfracCVNProngNumuErrAB2 =
606  ( (NumuMCWS[1][0]+NumuMCWS[1][1])*(NumuMCOther[1][0]+NumuMCOther[1][1]) ) *
607  ( (NumuMCWS[1][0]+NumuMCWS[1][1])*(NumuMCOther[1][0]+NumuMCOther[1][1]) ); //AB^2
608  float wsfracCVNProngNumuErr =
609  sqrt( C11CVNProngNumu*wsfracCVNProngNumuErrAB2*betaCVNProngNumu*betaCVNProngNumu/wsfracCVNProngNumuErrDenom +
610  C22CVNProngNumu*wsfracCVNProngNumuErrAB2*alphaCVNProngNumu*alphaCVNProngNumu/wsfracCVNProngNumuErrDenom -
611  2*C12CVNProngNumu*wsfracCVNProngNumuErrAB2*alphaCVNProngNumu*betaCVNProngNumu/wsfracCVNProngNumuErrDenom );
612  // BDT3
613  float numeratorBDT3Numu = NumuDataCt[2][0] - ((NumuMCWS[2][0]/NumuMCWS[2][1])*NumuDataCt[2][1]);
614  float denominatorBDT3Numu = NumuMCOther[2][0] - ((NumuMCWS[2][0]/NumuMCWS[2][1])*NumuMCOther[2][1]);
615  float alphaBDT3Numu = numeratorBDT3Numu / denominatorBDT3Numu;
616  float betaBDT3Numu = (NumuDataCt[2][1] - (alphaBDT3Numu*NumuMCOther[2][1]))/NumuMCWS[2][1];
617  // BDT3 error ellipses and WSFraction
618  float denomCovarBDT3Numu = (NumuMCOther[2][0]*NumuMCWS[2][1]) - (NumuMCOther[2][1]*NumuMCWS[2][0]);
619  float A11BDT3Numu = NumuMCWS[2][1] / denomCovarBDT3Numu;
620  float A12BDT3Numu = -NumuMCWS[2][0] / denomCovarBDT3Numu;
621  float A21BDT3Numu = -NumuMCOther[2][1] / denomCovarBDT3Numu;
622  float A22BDT3Numu = NumuMCOther[2][0] / denomCovarBDT3Numu;
623  float C11BDT3Numu = (A11BDT3Numu*A11BDT3Numu*NumuDataCt[2][0]) + (A12BDT3Numu*A12BDT3Numu*NumuDataCt[2][1]);
624  float C12BDT3Numu = (A11BDT3Numu*A21BDT3Numu*NumuDataCt[2][0]) + (A12BDT3Numu*A22BDT3Numu*NumuDataCt[2][1]);
625  float C21BDT3Numu = C12BDT3Numu;
626  float C22BDT3Numu = (A21BDT3Numu*A21BDT3Numu*NumuDataCt[2][0]) + (A22BDT3Numu*A22BDT3Numu*NumuDataCt[2][1]);
627  float majorBDT3Numu = sqrt(0.5*(C11BDT3Numu + C22BDT3Numu +
628  sqrt( (C11BDT3Numu+C22BDT3Numu)*(C11BDT3Numu+C22BDT3Numu)
629  -(4*( C11BDT3Numu*C22BDT3Numu - C12BDT3Numu*C21BDT3Numu )))));
630  float minorBDT3Numu = sqrt(0.5*(C11BDT3Numu + C22BDT3Numu -
631  sqrt( (C11BDT3Numu+C22BDT3Numu)*(C11BDT3Numu+C22BDT3Numu)
632  -(4*( C11BDT3Numu*C22BDT3Numu - C12BDT3Numu*C21BDT3Numu )))));
633  float rotBetaXBDT3Numu = 0.5*(180./TMath::Pi())*atan2(2.0*C12BDT3Numu,C22BDT3Numu - C11BDT3Numu);
634  float wsfracBDT3Numu =
635  (betaBDT3Numu*(NumuMCWS[2][0]+NumuMCWS[2][1])) /
636  (betaBDT3Numu*(NumuMCWS[2][0]+NumuMCWS[2][1]) + alphaBDT3Numu*(NumuMCOther[2][0]+NumuMCOther[2][1]));
637  float wsfracBDT3NumuErrDenom =
638  (betaBDT3Numu*(NumuMCWS[2][0]+NumuMCWS[2][1]) + alphaBDT3Numu*(NumuMCOther[2][0]+NumuMCOther[2][1])) *
639  (betaBDT3Numu*(NumuMCWS[2][0]+NumuMCWS[2][1]) + alphaBDT3Numu*(NumuMCOther[2][0]+NumuMCOther[2][1])) *
640  (betaBDT3Numu*(NumuMCWS[2][0]+NumuMCWS[2][1]) + alphaBDT3Numu*(NumuMCOther[2][0]+NumuMCOther[2][1])) *
641  (betaBDT3Numu*(NumuMCWS[2][0]+NumuMCWS[2][1]) + alphaBDT3Numu*(NumuMCOther[2][0]+NumuMCOther[2][1])); // term^4
642  float wsfracBDT3NumuErrAB2 =
643  ( (NumuMCWS[2][0]+NumuMCWS[2][1])*(NumuMCOther[2][0]+NumuMCOther[2][1]) ) *
644  ( (NumuMCWS[2][0]+NumuMCWS[2][1])*(NumuMCOther[2][0]+NumuMCOther[2][1]) ); //AB^2
645  float wsfracBDT3NumuErr =
646  sqrt( C11BDT3Numu*wsfracBDT3NumuErrAB2*betaBDT3Numu*betaBDT3Numu/wsfracBDT3NumuErrDenom +
647  C22BDT3Numu*wsfracBDT3NumuErrAB2*alphaBDT3Numu*alphaBDT3Numu/wsfracBDT3NumuErrDenom -
648  2*C12BDT3Numu*wsfracBDT3NumuErrAB2*alphaBDT3Numu*betaBDT3Numu/wsfracBDT3NumuErrDenom );
649  std::cout << "===========================================" << std::endl;
650  std::cout << " SYSTEM OF EQUATIONS RESULTS: NUMU " << std::endl;
651  std::cout << " All selected events " << std::endl;
652  std::cout << "-------------------------------------------" << std::endl;
653  std::cout << "alpha*BKG_mode + beta*SIG_mode = DATA_mode" << std::endl;
654  std::cout << "===========================================" << std::endl;
655  std::cout << "FSProton:" << std::endl;
656  std::cout << " alpha = " << alphaFSProtonNumu << std::endl;
657  std::cout << " beta = " << betaFSProtonNumu << std::endl;
658  std::cout << " Cov Matrix: | " << C11FSProtonNumu << " " << C12FSProtonNumu << " |" << std::endl;
659  std::cout << " | " << C21FSProtonNumu << " " << C22FSProtonNumu << " |" << std::endl;
660  std::cout << " sqrt(Eigvl): " << majorFSProtonNumu << std::endl;
661  std::cout << " " << minorFSProtonNumu << std::endl;
662  std::cout << " wsfrac = " << wsfracFSProtonNumu << " +/- " << wsfracFSProtonNumuErr << std::endl;
663  std::cout << std::endl;
664  std::cout << "CVNProng:" << std::endl;
665  std::cout << " alpha = " << alphaCVNProngNumu << std::endl;
666  std::cout << " beta = " << betaCVNProngNumu << std::endl;
667  std::cout << " Cov Matrix: | " << C11CVNProngNumu << " " << C12CVNProngNumu << " |" << std::endl;
668  std::cout << " | " << C21CVNProngNumu << " " << C22CVNProngNumu << " |" << std::endl;
669  std::cout << " sqrt(Eigvl): " << majorCVNProngNumu << std::endl;
670  std::cout << " " << minorCVNProngNumu << std::endl;
671  std::cout << " wsfrac = " << wsfracCVNProngNumu << " +/- " << wsfracCVNProngNumuErr << std::endl;
672  std::cout << std::endl;
673  std::cout << "BDT3:" << std::endl;
674  std::cout << " alpha = " << alphaBDT3Numu << std::endl;
675  std::cout << " beta = " << betaBDT3Numu << std::endl;
676  std::cout << " Cov Matrix: | " << C11BDT3Numu << " " << C12BDT3Numu << " |" << std::endl;
677  std::cout << " | " << C21BDT3Numu << " " << C22BDT3Numu << " |" << std::endl;
678  std::cout << " sqrt(Eigvl): " << majorBDT3Numu << std::endl;
679  std::cout << " " << minorBDT3Numu << std::endl;
680  std::cout << " wsfrac = " << wsfracBDT3Numu << " +/- " << wsfracBDT3NumuErr << std::endl;
681  std::cout << std::endl;
682  std::cout << std::endl;
683 
684  alphaNumu[0] = alphaFSProtonNumu;
685  alphaNumu[1] = alphaBDT3Numu;
686  betaNumu[0] = betaFSProtonNumu;
687  betaNumu[1] = betaBDT3Numu;
688  }
689 
690  // calculate chi2/dof in numu sample
691  TCanvas *cNumuErrorEst[2];
692  float errorScaleTerm[2];
693 
694  if( splitNCmu[0] && splitNCmu[1] ){
695  for( int iMode=0; iMode<2; ++iMode ){
696  int iAna=-1;
697  if( iMode==0 ) iAna=0;
698  if( iMode==1 ) iAna=2;
699  if( iAna==-1 ){
700  std::cout << "Skipping iMode==" << iMode << std::endl;
701  continue;
702  }
703  // Modified Predictions to compare to data
704  numuRS_SCalc[iAna][0]->Scale(betaNumu[iMode]);
705  numuErrorEstMC[iAna][0]->Scale(betaNumu[iMode]); // base of numuErrorEst is numuRS_SCalc[iAna][0], the wrong sign numu component
706  numuErrorEstMC[iAna][0]->SetName("rescaledMC_RSConsistent");
707  numuRS_SCalc[iAna][1]->Scale(betaNumu[iMode]);
708  numuRS_SCalc[iAna][2]->Scale(betaNumu[iMode]);
709  numuRS_BCalc[iAna][0]->Scale(alphaNumu[iMode]);
710  numuRS_BCalc[iAna][1]->Scale(alphaNumu[iMode]);
711  numuRS_BCalc[iAna][2]->Scale(alphaNumu[iMode]);
712  numuWS_SCalc[iAna][0]->Scale(betaNumu[iMode]);
713  numuErrorEstMC[iAna][1]->Scale(betaNumu[iMode]); // base of numuErrorEst is numuWS_SCalc[iAna][0], the wrong sign numu component
714  numuErrorEstMC[iAna][1]->SetName("rescaledMC_WSConsistent");
715  numuWS_SCalc[iAna][1]->Scale(betaNumu[iMode]);
716  numuWS_SCalc[iAna][2]->Scale(betaNumu[iMode]);
717  numuWS_BCalc[iAna][0]->Scale(alphaNumu[iMode]);
718  numuWS_BCalc[iAna][1]->Scale(alphaNumu[iMode]);
719  numuWS_BCalc[iAna][2]->Scale(alphaNumu[iMode]);
720  // Make plots! Summed up rescaled numu RS [0] and WS [1] components
721  numuErrorEstMC[iAna][0]->Add(numuRS_SCalc[iAna][1]);
722  numuErrorEstMC[iAna][0]->Add(numuRS_SCalc[iAna][2]);
723  numuErrorEstMC[iAna][0]->Add(numuRS_BCalc[iAna][0]);
724  numuErrorEstMC[iAna][0]->Add(numuRS_BCalc[iAna][1]);
725  numuErrorEstMC[iAna][0]->Add(numuRS_BCalc[iAna][2]);
726  numuErrorEstMC[iAna][1]->Add(numuWS_SCalc[iAna][1]);
727  numuErrorEstMC[iAna][1]->Add(numuWS_SCalc[iAna][2]);
728  numuErrorEstMC[iAna][1]->Add(numuWS_BCalc[iAna][0]);
729  numuErrorEstMC[iAna][1]->Add(numuWS_BCalc[iAna][1]);
730  numuErrorEstMC[iAna][1]->Add(numuWS_BCalc[iAna][2]);
731  // calculate chi2 / dof
732  int NbinsRS = numuRS_SCalc[iAna][0]->GetXaxis()->GetNbins();
733  int NbinsWS = numuWS_SCalc[iAna][0]->GetXaxis()->GetNbins();
734  float dof = float( NbinsRS + NbinsWS );
735  float chi2 = 0.;
736  // Loop over the RS first
737  for( int iBinChi2=1; iBinChi2<=NbinsRS; ++iBinChi2 ){
738  float modPred = numuErrorEstMC[iAna][0]->GetBinContent(iBinChi2);
739  float rsDataRes = numuRSDataResults[iAna]->GetBinContent(iBinChi2);
740  chi2 += (modPred-rsDataRes)*(modPred-rsDataRes)/rsDataRes;
741  }
742  // Loop over the WS now
743  for( int iBinChi2=1; iBinChi2<=NbinsWS; ++iBinChi2 ){
744  float modPred = numuErrorEstMC[iAna][1]->GetBinContent(iBinChi2);
745  float wsDataRes = numuWSDataResults[iAna]->GetBinContent(iBinChi2);
746  chi2 += (modPred-wsDataRes)*(modPred-wsDataRes)/wsDataRes;
747  }
748  std::cout << "-----------------------------------------" << std::endl;
749  std::cout << "ANA METHOD " << iAna << std::endl;
750  std::cout << "Chi2/dof = " << chi2 << "/" << dof << " = " << chi2/dof << std::endl;
751  std::cout << " so scale the error bars by " << sqrt(chi2/dof) << std::endl;
752  std::cout << "-----------------------------------------" << std::endl;
753  errorScaleTerm[iMode] = chi2/dof;
754  // Loop over the RS first, scale error bars
755  for( int iBinChi2=1; iBinChi2<=NbinsRS; ++iBinChi2 ){
756  numuErrorEstData[iAna][0]->SetBinError(iBinChi2,sqrt(chi2/dof)*numuErrorEstData[iAna][0]->GetBinError(iBinChi2));
757  }
758  // Loop over the WS now, scale error bars
759  for( int iBinChi2=1; iBinChi2<=NbinsWS; ++iBinChi2 ){
760  numuErrorEstData[iAna][1]->SetBinError(iBinChi2,sqrt(chi2/dof)*numuErrorEstData[iAna][1]->GetBinError(iBinChi2));
761  }
762  cNumuErrorEst[iMode] = new TCanvas();
763  cNumuErrorEst[iMode]->Divide(2,1);
764  cNumuErrorEst[iMode]->cd(1);
765  numuErrorEstMC[iAna][0]->SetLineColor(kRed);
766  // Format Titles
767  numuErrorEstMC[iAna][0]->Scale(0.1,"width"); // get it in Events/0.1 GeV (I think?)
768  numuErrorEstMC[iAna][0]->GetYaxis()->SetTitle("");
769  CenterTitles(numuErrorEstMC[iAna][0]);
770  numuErrorEstMC[iAna][1]->Scale(0.1,"width");
771  numuErrorEstMC[iAna][1]->GetYaxis()->SetTitle("");
772  numuErrorEstMC[iAna][1]->GetYaxis()->SetTitleOffset(1.1);
773  CenterTitles(numuErrorEstMC[iAna][1]);
774  // -------------------
775  // Y axis title
776  numuErrorEstMC[iAna][0]->Draw("hist");
777  double axisScaleFctr = 1.05*numuErrorEstMC[iAna][0]->GetMaximum();
778  TPaveText *yAxisBox = new TPaveText(3.75,axisScaleFctr/3.,4.75,2.*axisScaleFctr/3.,"NB");
779  TText *yTitle = yAxisBox->AddText("Events/0.1 GeV");
780  yTitle->SetTextAngle(90.);
781  yTitle->SetTextSize( numuErrorEstMC[iAna][0]->GetXaxis()->GetTitleSize() );
782  yAxisBox->SetFillColor(0);
783  yAxisBox->Draw();
784  numuErrorEstMC[iAna][0]->Draw("axissame");
785  // -------------------------
786  cNumuErrorEst[iMode]->cd(2);
787  numuErrorEstMC[iAna][1]->SetLineColor(kRed);
788  numuErrorEstMC[iAna][1]->Draw("hist");
789  cNumuErrorEst[iMode]->cd(1);
790  numuErrorEstData[iAna][0]->Scale(0.1,"width");
791  numuErrorEstData[iAna][0]->Draw("same");
792  cNumuErrorEst[iMode]->cd(2);
793  numuErrorEstData[iAna][1]->Scale(0.1,"width");
794  numuErrorEstData[iAna][1]->Draw("same");
795  Preliminary();
796  if( iAna==0 ) cNumuErrorEst[iMode]->Print("FSProton_RescaledResult_q3.pdf");
797  if( iAna==2 ) cNumuErrorEst[iMode]->Print("BDT_RescaledResult_q3.pdf");
798  cNumuErrorEst[iMode]->Write(TString::Format("RescaledResults_Ana%d",iAna));
799  }
800  // Re-calculate the terms
801  // FSProton
802  float numeratorFSProtonNumu = NumuDataCt[0][0] - ((NumuMCWS[0][0]/NumuMCWS[0][1])*NumuDataCt[0][1]);
803  float denominatorFSProtonNumu = NumuMCOther[0][0] - ((NumuMCWS[0][0]/NumuMCWS[0][1])*NumuMCOther[0][1]);
804  float alphaFSProtonNumu = numeratorFSProtonNumu / denominatorFSProtonNumu;
805  float betaFSProtonNumu = (NumuDataCt[0][1] - (alphaFSProtonNumu*NumuMCOther[0][1]))/NumuMCWS[0][1];
806  float denomCovarFSProtonNumu = (NumuMCOther[0][0]*NumuMCWS[0][1]) - (NumuMCOther[0][1]*NumuMCWS[0][0]);
807  float A11FSProtonNumu = NumuMCWS[0][1] / denomCovarFSProtonNumu;
808  float A12FSProtonNumu = -NumuMCWS[0][0] / denomCovarFSProtonNumu;
809  float A21FSProtonNumu = -NumuMCOther[0][1] / denomCovarFSProtonNumu;
810  float A22FSProtonNumu = NumuMCOther[0][0] / denomCovarFSProtonNumu;
811  float C11FSProtonNumu =
812  (A11FSProtonNumu*A11FSProtonNumu*(NumuDataCt[0][0]*errorScaleTerm[0])) +
813  (A12FSProtonNumu*A12FSProtonNumu*(NumuDataCt[0][1]*errorScaleTerm[0]));
814  float C12FSProtonNumu =
815  (A11FSProtonNumu*A21FSProtonNumu*(NumuDataCt[0][0]*errorScaleTerm[0])) +
816  (A12FSProtonNumu*A22FSProtonNumu*(NumuDataCt[0][1]*errorScaleTerm[0]));
817  float C21FSProtonNumu = C12FSProtonNumu;
818  float C22FSProtonNumu =
819  (A21FSProtonNumu*A21FSProtonNumu*(NumuDataCt[0][0]*errorScaleTerm[0])) +
820  (A22FSProtonNumu*A22FSProtonNumu*(NumuDataCt[0][1]*errorScaleTerm[0]));
821  // BDT3
822  float numeratorBDT3Numu = NumuDataCt[2][0] - ((NumuMCWS[2][0]/NumuMCWS[2][1])*NumuDataCt[2][1]);
823  float denominatorBDT3Numu = NumuMCOther[2][0] - ((NumuMCWS[2][0]/NumuMCWS[2][1])*NumuMCOther[2][1]);
824  float alphaBDT3Numu = numeratorBDT3Numu / denominatorBDT3Numu;
825  float betaBDT3Numu = (NumuDataCt[2][1] - (alphaBDT3Numu*NumuMCOther[2][1]))/NumuMCWS[2][1];
826  float denomCovarBDT3Numu = (NumuMCOther[2][0]*NumuMCWS[2][1]) - (NumuMCOther[2][1]*NumuMCWS[2][0]);
827  float A11BDT3Numu = NumuMCWS[2][1] / denomCovarBDT3Numu;
828  float A12BDT3Numu = -NumuMCWS[2][0] / denomCovarBDT3Numu;
829  float A21BDT3Numu = -NumuMCOther[2][1] / denomCovarBDT3Numu;
830  float A22BDT3Numu = NumuMCOther[2][0] / denomCovarBDT3Numu;
831  float C11BDT3Numu =
832  (A11BDT3Numu*A11BDT3Numu*(NumuDataCt[2][0]*errorScaleTerm[1])) + (A12BDT3Numu*A12BDT3Numu*(NumuDataCt[2][1]*errorScaleTerm[1]));
833  float C12BDT3Numu =
834  (A11BDT3Numu*A21BDT3Numu*(NumuDataCt[2][0]*errorScaleTerm[1])) + (A12BDT3Numu*A22BDT3Numu*(NumuDataCt[2][1]*errorScaleTerm[1]));
835  float C21BDT3Numu = C12BDT3Numu;
836  float C22BDT3Numu =
837  (A21BDT3Numu*A21BDT3Numu*(NumuDataCt[2][0]*errorScaleTerm[1])) + (A22BDT3Numu*A22BDT3Numu*(NumuDataCt[2][1]*errorScaleTerm[1]));
838  std::cout << "Mode\talpha\tbeta\tC11[mode]\tC22[mode]\tC12[mode]" << std::endl;
839  std::cout << "FSPr" <<"\t"<< alphaFSProtonNumu <<"\t"<< betaFSProtonNumu <<"\t"<< C11FSProtonNumu
840  <<"\t"<< C22FSProtonNumu <<"\t"<< C12FSProtonNumu << std::endl;
841  std::cout << "BDT3" <<"\t"<< alphaBDT3Numu <<"\t"<< betaBDT3Numu <<"\t"<< C11BDT3Numu
842  <<"\t"<< C22BDT3Numu <<"\t"<< C12BDT3Numu << std::endl;
843  }
844 
845  std::cout << "Done." << std::endl;
846 
847  outFile->Close();
848 }
Simulation
void Simulation()
Definition: tools.h:16
std::max
T max(const caf::Proxy< T > &a, T b)
Definition: BasicTypesProxy.h:397
osc::_NoOscillations
Pass neutrinos through unchanged.
Definition: IOscCalc.h:40
ana
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
ana::Sign::kAntiNu
Antineutrinos-only.
Definition: IPrediction.h:50
ana::Spectrum::ToTH1
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
calculateWrongSignNumuQ3
void calculateWrongSignNumuQ3()
Definition: calculateWrongSignNumuQ3.C:148
std::sqrt
T sqrt(T number)
Definition: d0nt_math.hpp:156
ana::PredictionNoOsc::LoadFrom
static std::unique_ptr< PredictionNoOsc > LoadFrom(TDirectory *dir, const std::string &name)
Definition: PredictionNoOsc.cxx:105
ana::CenterTitles
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
beta
Double_t beta
Definition: TG43_relative_dose.C:15
calc
osc::OscCalcDumb calc
Definition: make_nueFDprediction_kinematics.C:70
chi2
double chi2()
Definition: cana_time_peak.C:532
ana::Spectrum
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
inFile
ifstream inFile
Definition: AnaPlotMaker.h:34
ana::Current::kCC
Charged-current interactions.
Definition: IPrediction.h:39
ana::PredictionNoOsc::PredictComponent
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
Definition: PredictionNoOsc.cxx:47
check_time_usage.float
float
Definition: check_time_usage.py:6
outFile
TFile * outFile
Definition: PlotXSec.C:135
ana::Spectrum::LoadFrom
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
GetXaxis
correl_xv GetXaxis() -> SetDecimals()
Definition: plot_lowE_sig_spread_genie.C:84
pot
#define pot
Definition: template_basic.C:49
calcAlphaBetaEachBin
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)
Definition: calculateWrongSignNumuQ3.C:61
ana::Spectrum::POT
double POT() const
Definition: Spectrum.h:219
Preliminary
void Preliminary()
Definition: draw_decomp_plots.C:55
ana::Sign::kNu
Neutrinos-only.
Definition: IPrediction.h:49
om::cout
OStream cout
Definition: OStream.cxx:6
POT
std::vector< double > POT
Definition: neutKEsyst_plot.C:45
xml_schema::string
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
ana::Current::kNC
Neutral-current interactions.
Definition: IPrediction.h:40
make_mec_shifts_plots.GetMaximum
def GetMaximum(hist)
Definition: make_mec_shifts_plots.py:20
ana::Sign::kBoth
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
genie::utils::style::Format
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
ana::Flavors::kAll
All neutrinos, any flavor.
Definition: IPrediction.h:26
ana::PredictionNoOsc
Prediction that wraps a simple Spectrum.
Definition: PredictionNoOsc.h:11
std::atan2
T atan2(T number)
Definition: d0nt_math.hpp:72
Generated on Mon Oct 26 2020 05:13:18 for NOvA by   doxygen 1.8.11