systematics_summary_from_pred_interp.C
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Macro for making the nue systematics bar plots
4 //
5 ////////////////////////////////////////////////////////////////////////////////
6 
7 //#include "CAFAna/Analysis/Calcs.h"
8 //#include "CAFAna/Analysis/Exposures.h"
9 //#include "CAFAna/Analysis/Plots.h"
10 
11 //#include "OscLib/IOscCalc.h"
12 //#include "OscLib/OscCalcDumb.h"
13 
14 //#include "3FlavorAna/Ana2020/joint_fit_style_tools.h"
15 //#include "3FlavorAna/Ana2020/joint_fit_2020_loader_tools.h"
16 //#include "CAFAna/Analysis/3FlavorAna2020Systs.h"
17 
18 //#include "3FlavorAna/Systs/DummySystStorage.h"
20 
21 #include "TCanvas.h"
22 #include "TGraphAsymmErrors.h"
23 
24 #include <algorithm>
25 #include <string>
26 #include <iostream>
27 #include <iomanip>
28 
29 using namespace ana;
30 
31 void systematics_summary_from_pred_interp(std::string beam = "rhc", bool IsNuMu = false, bool sum = true, bool extrapComp = false, bool quantsComp = false)
32 {
33  std::vector<ColumnDef> cols;
34  double pot = 0;
35  double livetime = 0;
36  const bool beamLab = beam == "fhc" ? true : false;
37  std::string sampLab = "";
38  std::vector<std::vector<const ISyst *>> systs;
39  std::vector<std::string> extrapStr;
40  std::vector<std::string> extrapLab;
41 
42  // Do you want to make extrap/noextrap comparisons?
43  if(extrapComp)
44  {
45  extrapStr.push_back("noextrap_noPt");
46  extrapLab.push_back("Not Extrapolated");
47  }
48  std::string sNuMu = IsNuMu? "NuMu":"NuE";
49 
50  // Definitions of table columns follow, use the convention or you'll get wrong numbers:
51  // 0 = Total Pred
52  // 1 = Signal
53  // 2 = Total Bkg
54  // 3 = Beam Bkg
55  // last = Cosmics
56  if( beam == "fhc" ) {
57  pot = kAna2020FHCPOT;
58  livetime = kAna2020FHCLivetime;
59  if (IsNuMu) {
60  cols = {
61  {"Total Pred", Flavors::kAll, Current::kBoth, Sign::kBoth},
62  {"$\\nu_\\mu$ Signal", Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth},
63  {"Total Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth},
64  {"Beam Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth},
65  {"$\\nu_\\mu$ App.", Flavors::kNuEToNuMu,Current::kCC, Sign::kBoth},
66  {"NC", Flavors::kAll, Current::kNC, Sign::kBoth},
67  {"$\\nu_e$ CC", Flavors::kAllNuE, Current::kCC, Sign::kBoth},
68  {"$\\nu_\\tau$ CC", Flavors::kAllNuTau, Current::kCC, Sign::kBoth},
69  {"Cosmics", Flavors::kAll, Current::kBoth, Sign::kBoth}
70 
71  };
72  if(extrapComp) systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNumuAna, true, BeamType2020::kFHC, false, false ) );
73  systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNumuAna, true, BeamType2020::kFHC, false, true ) );
74  extrapStr.push_back("extrap");
75  extrapLab.push_back("Extrapolated");
76  } else {
77  cols = {
78  {"Total Pred", Flavors::kAll, Current::kBoth, Sign::kBoth},
79  {"$\\nu_e$ Signal", Flavors::kNuMuToNuE, Current::kCC, Sign::kNu},
80  {"Total Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth},
81  {"Beam Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth},
82  {"$\\bar\\nu_e$ WS", Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu},
83  {"beam $\\nu_e$ CC", Flavors::kNuEToNuE, Current::kCC, Sign::kBoth},
84  {"NC", Flavors::kAll, Current::kNC, Sign::kBoth},
85  {"$\\nu_\\mu$ CC", Flavors::kAllNuMu, Current::kCC, Sign::kBoth},
86  {"$\\nu_\\tau$ CC", Flavors::kAllNuTau, Current::kCC, Sign::kBoth},
87  {"Cosmics", Flavors::kAll, Current::kBoth, Sign::kBoth}
88  };
89  if(extrapComp)
90  {
91  systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNueAna, true, BeamType2020::kFHC, false, false ) );
92  systs.back().push_back(&kRockScaleSyst);
93  }
94  systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNueAna, true, BeamType2020::kFHC, false, true ) );
95  systs.back().push_back(&kRockScaleSyst);
96  extrapStr.push_back("combo");
97  extrapLab.push_back("Extrapolated");
98  }
99  } else { // if RHC
100  pot = kAna2020RHCPOT;
101  livetime = kAna2020RHCLivetime;
102  if (IsNuMu) {
103  cols = {
104  {"Total Pred", Flavors::kAll, Current::kBoth, Sign::kBoth},
105  {"$\\bar\\nu_\\mu$ Signal", Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth},
106  {"Total Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth},
107  {"Beam Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth},
108  {"$\\nu_\\mu$ App.", Flavors::kNuEToNuMu,Current::kCC, Sign::kBoth},
109  {"NC", Flavors::kAll, Current::kNC, Sign::kBoth},
110  {"$\\nu_e$ CC", Flavors::kAllNuE, Current::kCC, Sign::kBoth},
111  {"$\\nu_\\tau$ CC", Flavors::kAllNuTau, Current::kCC, Sign::kBoth},
112  {"Cosmics", Flavors::kAll, Current::kBoth, Sign::kBoth}
113  };
114  if(extrapComp) systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNumuAna, true, BeamType2020::kRHC, false, false ) );
115  systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNumuAna, true, BeamType2020::kRHC, false, true ) );
116  extrapStr.push_back("extrap");
117  extrapLab.push_back("Extrapolated");
118  } else {
119  cols = {
120  {"Total Pred", Flavors::kAll, Current::kBoth, Sign::kBoth},
121  {"$\\bar\\nu_e$ Signal",Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu},
122  {"Total Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth},
123  {"Beam Bkg", Flavors::kAll, Current::kBoth, Sign::kBoth},
124  {"$\\nu_e$ WS", Flavors::kNuMuToNuE, Current::kCC, Sign::kNu},
125  {"beam $\\nu_e$ CC", Flavors::kNuEToNuE, Current::kCC, Sign::kBoth},
126  {"NC", Flavors::kAll, Current::kNC, Sign::kBoth},
127  {"$\\nu_\\mu$ CC", Flavors::kAllNuMu, Current::kCC, Sign::kBoth},
128  {"$\\nu_\\tau$ CC", Flavors::kAllNuTau, Current::kCC, Sign::kBoth},
129  {"Cosmics", Flavors::kAll, Current::kBoth, Sign::kBoth}
130  };
131  if(extrapComp)
132  {
133  systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNueAna, true, BeamType2020::kRHC, false, false ) );
134  systs.back().push_back(&kRockScaleSyst);
135  }
136  systs.push_back( get3FlavorAna2020AllSysts( EAnaType2020::kNueAna, true, BeamType2020::kRHC, false, true ) );
137  systs.back().push_back(&kRockScaleSyst);
138  extrapStr.push_back("prop");
139  extrapLab.push_back("Extrapolated");
140  }
141  }
142  gSystem->MakeDirectory( "plots" );
143 
144  const unsigned int cosmicsIdx = cols.size()-1;
145 
146  std::cout << std::fixed;
147  std::cout << std::setprecision(2);
148 
149  // 2020 Best Fit
150  const double dCP = .82;
151  const double ssth23 = .568;
152  const double ss2th13 = .085;
153  const double dmsq32 = 2.406e-3;
154 
156  calcosc->SetdCP ( dCP*(TMath::Pi()) );
157  calcosc->SetTh23 ( asin(sqrt(ssth23)) );
158  calcosc->SetDmsq32( dmsq32 );
159  calcosc->SetTh13 ( .5*asin(sqrt(ss2th13)) );
160 
161  std::vector<std::map<std::string,std::pair<double,double>>> tots, sigs, bkgs;
162  std::vector<std::pair<double,double>> totQuad, sigQuad, bkgQuad;
163 
164  for(unsigned int ext = 0; ext < extrapStr.size(); ++ext){
165  std::vector< std::vector<Spectrum> > noms;
166  std::vector<const IPrediction*> pred;
167  std::vector< std::pair <Spectrum*, double> > cosmics;
168  if (IsNuMu) {
169  pred = GetNumuPredictions2020( 4, extrapStr[ext], calcosc, true, beam );
170  cosmics = GetNumuCosmics2020( 4, beam );
171  } else {
172  pred.push_back( GetNuePrediction2020( extrapStr[ext], calcosc, true, beam, false ) );
173  cosmics.push_back( GetNueCosmics2020( beam ) );
174  }
175  systs[ext].push_back(&kCosmicBkgScaleSyst);
176 
177  for (size_t pr=0; pr<pred.size(); ++pr) {
178  std::vector<Spectrum> Tempnom;
179  for(ColumnDef col: cols) {
180  if(col.name == "Cosmics")
181  {
182  Tempnom.push_back(*cosmics[pr].first);
183  std::cerr << "Cosmics integral is " << Tempnom.back().Integral(livetime, 0, kLivetime) << std::endl;
184  }
185  else
186  {
187  Tempnom.push_back(pred[pr]->PredictComponent(calcosc, col.flav, col.curr, col.sign));
188  std::cerr << "Tempnom integral is " << Tempnom.back().Integral(pot) << std::endl;
189  }
190  }
191  noms.push_back(Tempnom);
192  }
193 
194  // Sum of numu quartiles
195  if (IsNuMu)
196  {
197  std::cerr << "Summing numu nominal quantiles!" << std::endl;
198  std::vector<Spectrum> Tempnom;
199  for(ColumnDef col: cols) {
200  if(col.name == "Cosmics")
201  {
202  Spectrum sumCosmics = *cosmics[0].first + *cosmics[1].first + *cosmics[2].first + *cosmics[3].first;
203  Tempnom.push_back(sumCosmics);
204  std::cerr << "Cosmics integral is " << Tempnom.back().Integral(livetime, 0, kLivetime) << std::endl;
205  }
206  else
207  {
208  Spectrum sumBeam = pred[0]->PredictComponent(calcosc, col.flav, col.curr, col.sign) +
209  pred[1]->PredictComponent(calcosc, col.flav, col.curr, col.sign) +
210  pred[2]->PredictComponent(calcosc, col.flav, col.curr, col.sign) +
211  pred[3]->PredictComponent(calcosc, col.flav, col.curr, col.sign);
212  Tempnom.push_back(sumBeam);
213  std::cerr << "Tempnom integral is " << Tempnom.back().Integral(pot) << std::endl;
214  }
215  }
216  noms.push_back(Tempnom);
217  }
218 
220  std::cout<<"\\section{Best Fit Point}"<<std::endl;
221  std::cout<<"$\\delta / \\pi = " << dCP << "$\\\\" << std::endl;
222  std::cout<<"$\\sin^{2} \\theta_{23} = " << ssth23 << "$\\\\" << std::endl;
223  std::cout<<"$\\Delta m_{32}^{2} = " << dmsq32 << "$\\\\" << std::endl;
224  std::cout<<"$\\sin^{2} 2\\theta_{13} = " << ss2th13 << "$\\\\"<<std::endl;
225  std::cout<<"\\pagebreak"<<std::endl;
226 
227  unsigned int predSize = pred.size();
228  std::cerr << "predSize = " << predSize << std::endl;
229  if(IsNuMu) predSize += 1; // for summed quartiles prediction, add one more summing step
230  if(!IsNuMu && quantsComp) predSize += 3; // for low/high PID and peripheral bin, mimic the three additional bins
231  // --- And now I want to loop through my vectors...
232  for (unsigned int prIdx=0; prIdx<predSize; ++prIdx) {
233  unsigned int pr = prIdx;
234  if(!IsNuMu && quantsComp) pr = 0;
235 
236  TH1* hNomTot = noms[pr][0].ToTH1(pot);
237  // Add cosmics
238  hNomTot->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
239 
240  TH1* hNomSig = noms[pr][1].ToTH1(pot);
241 
242  if(IsNuMu || prIdx==0)
243  {
244  // Total bkg - signal
245  noms[pr][2] -= noms[pr][1];
246  // Beam bkg - signal
247  noms[pr][3] -= noms[pr][1];
248  }
249 
250  TH1* hNomBkg = noms[pr][2].ToTH1(pot);
251  hNomBkg->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
252 
253  std::map<std::string, TH1*> hUpsTot, hUpsSig, hUpsBkg;
254  std::map<std::string, TH1*> hDnsTot, hDnsSig, hDnsBkg;
255 
256  // need to add cosmic to total prediction and bkg prediction
257  double totErr = 100/sqrt(Integral(noms[pr][0], pot)+Integral(noms[pr][cosmicsIdx], livetime, 0, kLivetime));
258  double sigErr = 100/sqrt(Integral(noms[pr][1], pot));
259  double bkgErr = 100/sqrt(Integral(noms[pr][2], pot)+Integral(noms[pr][cosmicsIdx], livetime, 0, kLivetime));
260 
261  std::string sSamp = beam + " " + sNuMu;
262 
263  std::cout<<"\n\n"<<std::endl;
264  // -- Table Header.
265  if (IsNuMu) {
266  if(pr==4) {
267  sSamp += " All Quartiles";
268  std::cout<< "\\section{All Quartiles}"<<std::endl;
269  }
270  else {
271  sSamp += " Quartile " + std::to_string(pr+1);
272  std::cout<< "\\section{Ehad Quartile "<<std::to_string(pr+1)<<"}"<<std::endl;
273  }
274  std::cout<<"\\hspace*{-1.5cm}\\resizebox{1.25\\textwidth}{!}{"<<std::endl;
275  std::cout<<"\\small"<<std::endl;
276  std::cout<<"\\begin{tabular}{l c c c c c c c c c c c c c c c c c c}"<<std::endl;
277  std::cout<<"\\multicolumn{19}{c}{\\bf{Systematic Summary for " << sSamp << " (\\%)}} \\\\"<<std::endl;
278  } else {
279  std::string sBin = "All Prediction";
280  if(prIdx==1) sBin = "Low PID";
281  else if(prIdx==2) sBin = "High PID";
282  else if(prIdx==3) sBin = "Peripheral";
283  sSamp += " "+sBin;
284  std::cout<<"\\section{"<< sBin <<"}"<<std::endl;
285  std::cout<<"\\hspace*{-1.5cm}\\resizebox{1.25\\textwidth}{!}{"<<std::endl;
286  std::cout<<"\\small"<<std::endl;
287  std::cout<<"\\begin{tabular}{l c c c c c c c c c c c c c c c c c c c c}"<<std::endl;
288  std::cout<<"\\multicolumn{21}{c}{\\bf{Systematic Summary for " << sSamp << " (\\%)}} \\\\"<<std::endl;
289  }
290  // -- Top line of table
291  std::cout << "{Systematic}";
292  for(ColumnDef col: cols) std::cout << " & \\multicolumn{2}{c}{" << col.name << "}";
293  std::cout << " \\\\" << std::endl;
294  std::cout << "\\hline" << std::endl;
295 
296  std::pair<double,double> quads[cols.size()];
297  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx) quads[colIdx] = {0,0};
298 
299  std::map<std::string, std::pair<double,double>> barsTot, barsSig, barsBkg;
300 
301  for(const ISyst* syst: systs[ext]){
302 
303  std::cout << FixSystName(syst->ShortName());
304 
305  std::vector<Spectrum> shiftsUp, shiftsDn;
306  for(ColumnDef col: cols){
307  // not summing the quantiles (nue or not the last pr)
308  if(!IsNuMu || pr != pred.size())
309  {
310  if(col.name != "Cosmics")
311  {
312  if(syst == &kCosmicBkgScaleSyst)
313  {
314  shiftsUp.push_back(pred[pr]->PredictComponent(calcosc, col.flav, col.curr, col.sign));
315  shiftsDn.push_back(pred[pr]->PredictComponent(calcosc, col.flav, col.curr, col.sign));
316  }
317  else
318  {
319  shiftsUp.push_back(pred[pr]->PredictComponentSyst(calcosc, SystShifts(syst, +1), col.flav, col.curr, col.sign));
320  shiftsDn.push_back(pred[pr]->PredictComponentSyst(calcosc, SystShifts(syst, -1), col.flav, col.curr, col.sign));
321  }
322  }
323  else
324  {
325  shiftsUp.push_back(*cosmics[pr].first);
326  shiftsDn.push_back(*cosmics[pr].first);
327 
328  if(syst == &kCosmicBkgScaleSyst)
329  {
330  shiftsUp.back() = *ShiftedCosmics(cosmics[pr].first, +1);
331  shiftsDn.back() = *ShiftedCosmics(cosmics[pr].first, -1);
332  }
333  }
334  }
335  // summing the numu quantiles
336  else
337  {
338  if(col.name != "Cosmics")
339  {
340  if(syst == &kCosmicBkgScaleSyst)
341  {
342  shiftsUp.push_back(pred[0]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
343  pred[1]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
344  pred[2]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
345  pred[3]->PredictComponent(calcosc, col.flav, col.curr, col.sign));
346  shiftsDn.push_back(pred[0]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
347  pred[1]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
348  pred[2]->PredictComponent(calcosc, col.flav, col.curr, col.sign)+
349  pred[3]->PredictComponent(calcosc, col.flav, col.curr, col.sign));
350  }
351  else
352  {
353  shiftsUp.push_back(pred[0]->PredictComponentSyst(calcosc, SystShifts(syst, +1), col.flav, col.curr, col.sign)+
354  pred[1]->PredictComponentSyst(calcosc, SystShifts(syst, +1), col.flav, col.curr, col.sign)+
355  pred[2]->PredictComponentSyst(calcosc, SystShifts(syst, +1), col.flav, col.curr, col.sign)+
356  pred[3]->PredictComponentSyst(calcosc, SystShifts(syst, +1), col.flav, col.curr, col.sign));
357  shiftsDn.push_back(pred[0]->PredictComponentSyst(calcosc, SystShifts(syst, -1), col.flav, col.curr, col.sign)+
358  pred[1]->PredictComponentSyst(calcosc, SystShifts(syst, -1), col.flav, col.curr, col.sign)+
359  pred[2]->PredictComponentSyst(calcosc, SystShifts(syst, -1), col.flav, col.curr, col.sign)+
360  pred[3]->PredictComponentSyst(calcosc, SystShifts(syst, -1), col.flav, col.curr, col.sign));
361  }
362  }
363  else
364  {
365  if(syst == &kCosmicBkgScaleSyst)
366  {
367  Spectrum cosmicsUp = *ShiftedCosmics(cosmics[0].first, +1) +
368  *ShiftedCosmics(cosmics[1].first, +1) +
369  *ShiftedCosmics(cosmics[2].first, +1) +
370  *ShiftedCosmics(cosmics[3].first, +1);
371 
372  Spectrum cosmicsDn = *ShiftedCosmics(cosmics[0].first, -1) +
373  *ShiftedCosmics(cosmics[1].first, -1) +
374  *ShiftedCosmics(cosmics[2].first, -1) +
375  *ShiftedCosmics(cosmics[3].first, -1);
376 
377  shiftsUp.push_back(cosmicsUp);
378  shiftsDn.push_back(cosmicsDn);
379  }
380  else
381  {
382  shiftsUp.push_back(*cosmics[0].first+*cosmics[1].first+*cosmics[2].first+*cosmics[3].first);
383  shiftsDn.push_back(*cosmics[0].first+*cosmics[1].first+*cosmics[2].first+*cosmics[3].first);
384  }
385  }
386  }
387  }
388 
389  hUpsTot[ (std::string)FixSystName(syst->ShortName()) ] = shiftsUp[0].ToTH1(pot);
390  hDnsTot[ (std::string)FixSystName(syst->ShortName()) ] = shiftsDn[0].ToTH1(pot);
391  hUpsSig[ (std::string)FixSystName(syst->ShortName()) ] = shiftsUp[1].ToTH1(pot);
392  hDnsSig[ (std::string)FixSystName(syst->ShortName()) ] = shiftsDn[1].ToTH1(pot);
393  hUpsBkg[ (std::string)FixSystName(syst->ShortName()) ] = (shiftsUp[0] - shiftsUp[1]).ToTH1(pot);
394  hDnsBkg[ (std::string)FixSystName(syst->ShortName()) ] = (shiftsDn[0] - shiftsDn[1]).ToTH1(pot);
395  if(syst == &kCosmicBkgScaleSyst)
396  {
397  TH1D* cosUp;
398  TH1D* cosDn;
399  if(IsNuMu && pr==4)
400  {
401  cosUp = ShiftedCosmics(cosmics[0].first,+1)->ToTH1(livetime,kLivetime);
402  cosUp -> Add(ShiftedCosmics(cosmics[1].first,+1)->ToTH1(livetime,kLivetime));
403  cosUp -> Add(ShiftedCosmics(cosmics[2].first,+1)->ToTH1(livetime,kLivetime));
404  cosUp -> Add(ShiftedCosmics(cosmics[3].first,+1)->ToTH1(livetime,kLivetime));
405 
406  cosDn = ShiftedCosmics(cosmics[0].first,-1)->ToTH1(livetime,kLivetime);
407  cosDn -> Add(ShiftedCosmics(cosmics[1].first,-1)->ToTH1(livetime,kLivetime));
408  cosDn -> Add(ShiftedCosmics(cosmics[2].first,-1)->ToTH1(livetime,kLivetime));
409  cosDn -> Add(ShiftedCosmics(cosmics[3].first,-1)->ToTH1(livetime,kLivetime));
410  } else {
411  cosUp = ShiftedCosmics(cosmics[pr].first,+1)->ToTH1(livetime,kLivetime);
412  cosDn = ShiftedCosmics(cosmics[pr].first,-1)->ToTH1(livetime,kLivetime);
413  }
414  hUpsTot[ (std::string)FixSystName(syst->ShortName()) ]->Add(cosUp);
415  hDnsTot[ (std::string)FixSystName(syst->ShortName()) ]->Add(cosDn);
416  hUpsBkg[ (std::string)FixSystName(syst->ShortName()) ]->Add(cosUp);
417  hDnsBkg[ (std::string)FixSystName(syst->ShortName()) ]->Add(cosDn);
418  }
419  else
420  {
421  hUpsTot[ (std::string)FixSystName(syst->ShortName()) ]->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
422  hDnsTot[ (std::string)FixSystName(syst->ShortName()) ]->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
423  hUpsBkg[ (std::string)FixSystName(syst->ShortName()) ]->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
424  hDnsBkg[ (std::string)FixSystName(syst->ShortName()) ]->Add(noms[pr][cosmicsIdx].ToTH1(livetime, kLivetime));
425  }
426 
427  // Shifted Total bkg = Total pred with cosmics - signal
428  shiftsUp[2] -= shiftsUp[1];
429  shiftsDn[2] -= shiftsDn[1];
430 
431  // Shifted Beam bkg = Total pred - signal
432  shiftsUp[3] -= shiftsUp[1];
433  shiftsDn[3] -= shiftsDn[1];
434 
435  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx){
436  const Spectrum& nom = noms[pr][colIdx];
437  const Spectrum& cosNom = noms[pr][cosmicsIdx];
438 
439  double ratioUp, errUp, ratioDn, errDn;
440  if(colIdx == 0 || colIdx == 2)
441  {
442  ratioUp = 1.0*( Integral(shiftsUp[colIdx], pot) + Integral(shiftsUp[cosmicsIdx], livetime, 0, kLivetime) ) / (1.0*(Integral(nom, pot)+Integral(cosNom,livetime,0,kLivetime)));
443  errUp = 100.0*(ratioUp-1.0);
444 
445  ratioDn = 1.0*( Integral(shiftsDn[colIdx], pot) + Integral(shiftsDn[cosmicsIdx], livetime, 0, kLivetime) ) / (1.0*(Integral(nom, pot)+Integral(cosNom,livetime,0,kLivetime)));
446  errDn = 100.0*(ratioDn-1.0);
447  }
448  else if(colIdx == cosmicsIdx)
449  {
450  ratioUp = 1.0*( Integral(shiftsUp[colIdx],livetime,0,kLivetime)) / (1.0*(Integral(nom,livetime,0,kLivetime)));
451  errUp = 100.0*(ratioUp-1.0);
452 
453  ratioDn = 1.0*( Integral(shiftsDn[colIdx],livetime,0,kLivetime)) / (1.0*(Integral(nom,livetime,0,kLivetime)));
454  errDn = 100.0*(ratioDn-1.0);
455  }
456  else
457  {
458  ratioUp = 1.0*Integral(shiftsUp[colIdx], pot) / (1.0*Integral(nom, pot));
459  errUp = 100.0*(ratioUp-1.0);
460 
461  ratioDn = 1.0*Integral(shiftsDn[colIdx], pot) / (1.0*Integral(nom, pot));
462  errDn = 100.0*(ratioDn-1.0);
463  }
464 
465  // if investigating the individual nue bins, we need to integrate over the particular bins of the predicted spectrum
466  if(!IsNuMu && quantsComp && prIdx != 0)
467  {
468  std::pair<int,int> Int = {1,9}; // low PID
469  if(prIdx==2) Int = {10,18}; // high PID
470  if(prIdx==3) Int = {19,hNomTot->GetNbinsX()}; // peripheral
471 
472  if(colIdx == 0 || colIdx == 2)
473  {
474  ratioUp = 1.0*( shiftsUp[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) + shiftsUp[cosmicsIdx].ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(pot)->Integral(Int.first,Int.second) + cosNom.ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) );
475  errUp = 100.0*(ratioUp-1.0);
476 
477  ratioDn = 1.0*( shiftsDn[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) + shiftsDn[cosmicsIdx].ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(pot)->Integral(Int.first,Int.second) + cosNom.ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) );
478  errDn = 100.0*(ratioDn-1.0);
479  }
480  else if(colIdx == cosmicsIdx)
481  {
482  ratioUp = 1.0*( shiftsUp[colIdx].ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) );
483  errUp = 100.0*(ratioUp-1.0);
484 
485  ratioDn = 1.0*( shiftsDn[colIdx].ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(livetime,kLivetime)->Integral(Int.first,Int.second) ) );
486  errDn = 100.0*(ratioDn-1.0);
487  }
488  else
489  {
490  ratioUp = 1.0*( shiftsUp[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(pot)->Integral(Int.first,Int.second) ) );
491  errUp = 100.0*(ratioUp-1.0);
492 
493  ratioDn = 1.0*( shiftsDn[colIdx].ToTH1(pot)->Integral(Int.first,Int.second) ) / ( 1.0*( nom.ToTH1(pot)->Integral(Int.first,Int.second) ) );
494  errDn = 100.0*(ratioDn-1.0);
495  }
496  }
497 
498  double min = std::min(errUp,errDn);
499  double max = std::max(errUp,errDn);
500  if(max < 0.01) max = 0;
501  if(min > -0.01) min = 0;
502 
503  quads[colIdx].first += min*min;
504  quads[colIdx].second += max*max;
505 
506 
507  std::string minstr = "white";
508  std::string maxstr = "white";
509  /*
510  if( fabs(min) >= 1 && fabs(min) < 2) minstr="Apricot";
511  if( fabs(min) >= 2 && fabs(min) < 3) minstr="YellowOrange";
512  if( fabs(min) >= 3 ) minstr="RedOrange";
513 
514  if( fabs(max) >= 1 && fabs(max) < 2) maxstr="Apricot";
515  if( fabs(max) >= 2 && fabs(max) < 3) maxstr="YellowOrange";
516  if( fabs(max) >= 3 ) maxstr="RedOrange";
517  */
518 
519  std::cout << " & \\cellcolor{" << minstr << "}" << (min == 0 ? "-" : "") << min << " & \\cellcolor{"<< maxstr << "}+" << max;
520 
521  min = fabs(min);
522  max = fabs(max);
523  if(cols[colIdx].name == "$\\nu_e$ Signal" ||
524  cols[colIdx].name == "$\\bar\\nu_e$ Signal" ||
525  cols[colIdx].name == "$\\nu_\\mu$ Signal" ||
526  cols[colIdx].name == "$\\bar\\nu_\\mu$ Signal")
527  barsSig[ (std::string)FixSystName(syst->ShortName()) ] = {min,max};
528 
529  if(cols[colIdx].name == "Total Pred")
530  barsTot[ (std::string)FixSystName(syst->ShortName()) ] = {min,max};
531 
532  if(cols[colIdx].name == "Total Bkg")
533  barsBkg[ (std::string)FixSystName(syst->ShortName()) ] = {min,max};
534  } // end for colIdx
535 
536  std::cout << " \\\\" << std::endl;
537  } // end for syst
538 
539  std::cout << "\\hline" << std::endl;
540 
541  std::cout << "Quadrature";
542  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx) std::cout << " & -" << sqrt(quads[colIdx].first) << " & +" << sqrt(quads[colIdx].second);
543  std::cout << std::endl;
544  std::cout<<"\\end{tabular}"<<std::endl;
545  std::cout<<"}"<<std::endl;
546 
547  if(sum){
548  barsTot = SumSysts(barsTot);
549  barsSig = SumSysts(barsSig);
550  barsBkg = SumSysts(barsBkg);
551 
553  std::cout<<"\\vfill" << std::endl;
554  std::cout<<"Summing..."<<std::endl;
556 
557  std::cout<<"\\begin{tabular}{l c c c c c c}"<<std::endl;
558  std::cout<<"\\multicolumn{7}{c}{\\bf{Summed Systematic for " << sSamp << " (\\%)}} \\\\"<<std::endl;
559  std::cout<<"{Systematic}";
560  for(int colIdx = 0; colIdx < 3; ++colIdx) std::cout << " & \\multicolumn{2}{c}{" << cols[colIdx].name << "}";
561  std::cout<<" \\\\"<<std::endl;
562  std::cout << "\\hline" << std::endl;
563 
564  for(int i = 0;i < kNumSystTypes; ++i)
565  {
566  std::cout << systTypes[i] << " & "
567  << "-" << barsTot[systTypes[i]].first << " & +" << barsTot[systTypes[i]].second << " & "
568  << "-" << barsSig[systTypes[i]].first << " & +" << barsSig[systTypes[i]].second << " & "
569  << "-" << barsBkg[systTypes[i]].first << " & +" << barsBkg[systTypes[i]].second <<" \\\\"
570  << std::endl;
571  }
572 
573  std::cout << "\\hline" << std::endl;
574  std::cout << "Quadrature";
575  for(int colIdx = 0; colIdx < 3; ++colIdx) std::cout << " & -" << sqrt(quads[colIdx].first) << " & +" << sqrt(quads[colIdx].second);
577  std::cout<<"\\end{tabular}"<<std::endl;
579  }
580 
581  // tots, sigs and bkg have following struct:
582  // nue: extrap1_all, extrap1_lowPID, extrap1_highPID, extrap1_peripheral, extrap2_all...
583  // numu: extrap1_Q1, extrap1_Q2, extrap1_Q3, extrap1_Q4, extrap1_allQuants, extrap2_Q1...
584  tots.push_back(barsTot);
585  sigs.push_back(barsSig);
586  bkgs.push_back(barsBkg);
587 
588  for(unsigned int colIdx = 0; colIdx < cols.size(); ++colIdx){
589  quads[colIdx].first = sqrt(quads[colIdx].first);
590  quads[colIdx].second = sqrt(quads[colIdx].second);
591  }
592 
593  totQuad.push_back(quads[0]);
594  sigQuad.push_back(quads[1]);
595  bkgQuad.push_back(quads[2]);
596 
597  std::string sumStr = sum == true ? "short" : "full";
598  std::string sQuant = "";
599  if (IsNuMu) {
600  sQuant = "_Q"+std::to_string(pr+1);
601  if(pr==pred.size()) sQuant = "_AllQuartiles";
602  }
603  else
604  {
605  sQuant = "_AllPrediction";
606  if(prIdx==1) sQuant = "_lowPID";
607  else if(prIdx==2) sQuant = "_highPID";
608  else if(prIdx==3) sQuant = "_peripheral";
609  }
610 
611  std::cout << "\n\n\n And now to save things..." << std::endl;
612 
613  sampLab = "#nu_{e} Selection";
614  if(beam == "rhc") sampLab = "#bar{#nu}_{e} Selection";
615  if(IsNuMu)
616  {
617  std::string quart = "All Quartiles";
618  if(pr<4) quart ="Quartile "+std::to_string(pr+1);
619  sampLab = "#nu_{#mu} Selection "+quart;
620  if(beam == "rhc") sampLab = "#bar{#nu}_{#mu} Selection "+quart;
621  }
622  else
623  {
624  if(prIdx==1) sampLab += " Low PID";
625  if(prIdx==2) sampLab += " High PID";
626  if(prIdx==3) sampLab += " Peripheral";
627  }
628 
629  // Energy prediction ratios
630  TH1* hD = (TH1*)hNomTot->Clone();
631  TH1* hNomTotRatio = (TH1*)hNomTot->Clone();
632  hNomTotRatio->Divide(hNomTot);
633  std::vector<TH1*> tempUp, tempDn;
634  for(auto it: hUpsTot) tempUp.push_back(it.second);
635  for(auto it: hDnsTot) tempDn.push_back(it.second);
636 
637  std::vector<TH1*> hUpsTotRatios = MakeRatios(hNomTot, tempUp);
638  std::vector<TH1*> hDnsTotRatios = MakeRatios(hNomTot, tempDn);
639 
640  std::map<std::string, std::vector<TH1*>> hUpsTotType = SortSystHists(hUpsTot);
641  std::map<std::string, std::vector<TH1*>> hDnsTotType = SortSystHists(hDnsTot);
642 
643  std::vector<int> types = {0,1,2,4,5};
644  if(IsNuMu) types = {1,2,4,5,6};
645  TLegend *lTypes = new TLegend(.1,.1,.9,.9);
646  lTypes->SetFillStyle(0);
647 
648  TCanvas *c = new TCanvas(TString("c")+TString(sQuant),"c");
649  MyPlotWithSystErrorBand(hNomTotRatio, hUpsTotRatios, hDnsTotRatios, kBlack, kTotalMCErrorBandColor, 2, true);
650  hD->SetLineColor(kTotalMCErrorBandColor);
651  lTypes->AddEntry(hD,"Total syst. uncertainty","l");
652 
653  for(auto j: types)
654  {
655  std::vector<TH1*> hUpsRatios = MakeRatios(hNomTot, hUpsTotType[systTypes[j]]);
656  std::vector<TH1*> hDnsRatios = MakeRatios(hNomTot, hDnsTotType[systTypes[j]]);
657 
658  MyPlotWithSystErrorBand(hNomTotRatio, hUpsRatios, hDnsRatios, kBlack, systCols[j], 1, false);
659  TH1* hDummy = (TH1*)hNomTot->Clone();
660  hDummy->SetLineColor(systCols[j]);
661  lTypes->AddEntry(hDummy,systTypes[j].c_str(),"l");
662  }
663  MyPlotWithSystErrorBand(hNomTotRatio, hUpsTotRatios, hDnsTotRatios, kBlack, kTotalMCErrorBandColor, 2, false);
664  drawBeamLabel(beamLab,0.1);
665  drawSampleLabel(sampLab);
666  preliminary();
667  gPad->Print(("plots/ratioband_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_tot.pdf").c_str());
668 
669  TCanvas *cLeg = new TCanvas(TString("cLeg")+TString(sQuant),"cLeg");
670 
671  lTypes->Draw();
672  gPad->Print(("plots/ratioband_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_legend.pdf").c_str());
673 
674  // Classic bar-plots
675  TCanvas *cTot = new TCanvas(TString("cTot")+TString(sQuant),"cTot");
676  ErrorBarChart(barsTot,{totErr,totErr},"Total Prediction Uncertainty (%)");
677  gPad->SetLeftMargin(.35);
678  drawBeamLabel(beamLab);
679  drawSampleLabel(sampLab);
680  preliminary();
681  gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_tot.pdf").c_str());
682 
683  TCanvas *cTot2 = new TCanvas(TString("cTot2")+TString(sQuant),"cTot2");
684  ErrorBarChart(barsTot,{0,0},"Total Prediction Uncertainty (%)");
685  gPad->SetLeftMargin(.35);
686  drawBeamLabel(beamLab);
687  drawSampleLabel(sampLab);
688  preliminary();
689  gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_tot_nostat.pdf").c_str());
690 
691  TCanvas *cSig = new TCanvas(TString("cSig")+TString(sQuant),"cSig");
692  ErrorBarChart(barsSig,{sigErr,sigErr},"Signal Uncertainty (%)");
693  gPad->SetLeftMargin(.35);
694  drawBeamLabel(beamLab);
695  drawSampleLabel(sampLab);
696  preliminary();
697  gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_sig.pdf").c_str());
698 
699  TCanvas *cSig2 = new TCanvas(TString("cSig2")+TString(sQuant),"cSig2");
700  ErrorBarChart(barsSig,{0,0},"Signal Uncertainty (%)");
701  gPad->SetLeftMargin(.35);
702  drawBeamLabel(beamLab);
703  drawSampleLabel(sampLab);
704  preliminary();
705  gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_sig_nostat.pdf").c_str());
706 
707  TCanvas *cBkg = new TCanvas(TString("cBkg")+TString(sQuant),"cBkg");
708  ErrorBarChart(barsBkg,{bkgErr,bkgErr},"Background Uncertainty (%)");
709  gPad->SetLeftMargin(.35);
710  drawBeamLabel(beamLab);
711  drawSampleLabel(sampLab);
712  preliminary();
713  gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_bkg.pdf").c_str());
714 
715  TCanvas *cBkg2 = new TCanvas(TString("cBkg2")+TString(sQuant),"cBkg2");
716  ErrorBarChart(barsBkg,{0,0},"Background Uncertainty (%)");
717  gPad->SetLeftMargin(.35);
718  drawBeamLabel(beamLab);
719  drawSampleLabel(sampLab);
720  preliminary();
721  gPad->Print(("plots/syst_summary_"+sNuMu+"_"+extrapStr[ext]+"_"+sumStr+"_"+beam+sQuant+"_bkg_nostat.pdf").c_str());
722  } // Pred loop
723  std::cout<<"\n"<<std::endl;
724 
725  if (!IsNuMu) std::cout<<"\\end{landscape}"<<std::endl;
726  std::cout<<"\\end{document}"<<std::endl;
727  } // extrapStr loop
728 
729  // numu quartiles comparison
730  if(IsNuMu && quantsComp)
731  {
732  unsigned int nQuants = 4; // 4
733 
734  for(unsigned int i = 0; i < extrapStr.size(); ++i)
735  {
736  std::vector<std::string> quantLab = {"Quartile 1", "Quartile 2", "Quartile 3", "Quartile 4"};
737  std::string sumStr = "full";
738  std::vector<std::string> labels;
739  std::vector<std::vector<std::pair<double,double>>> fTot;
740  std::vector<std::vector<std::pair<double,double>>> fSig;
741  std::vector<std::vector<std::pair<double,double>>> fBkg;
742  std::vector<std::pair<double,double>> fTotQuad;
743  std::vector<std::pair<double,double>> fSigQuad;
744  std::vector<std::pair<double,double>> fBkgQuad;
745 
746  for(unsigned int q = 0; q < nQuants; ++q)
747  {
748  std::vector<std::pair<double,double>> tempTot;
749  std::vector<std::pair<double,double>> tempSig;
750  std::vector<std::pair<double,double>> tempBkg;
751  int idx = i*(nQuants+1) + q;
752  std::cerr << "Quantile idx: " << idx << " " << extrapStr[i] << std::endl;
753 
754  if(sum)
755  {
756  sumStr = "short";
757  for(int j = 0; j < kNumSystTypes; ++j){
758  if(q == 0) labels.push_back(systTypes[j]);
759  tempTot.push_back(tots[idx][systTypes[j]]);
760  tempSig.push_back(sigs[idx][systTypes[j]]);
761  tempBkg.push_back(bkgs[idx][systTypes[j]]);
762  }
763  }
764  else
765  {
766  for(const ISyst* syst: systs[0]){
767  if(q == 0) labels.push_back( (std::string)FixSystName(syst->ShortName()) );
768  tempTot.push_back(tots[idx][ (std::string)FixSystName(syst->ShortName()) ]);
769  tempSig.push_back(sigs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
770  tempBkg.push_back(bkgs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
771  }
772  }
773 
774 
775  fTot.push_back(tempTot);
776  fSig.push_back(tempSig);
777  fBkg.push_back(tempBkg);
778  fTotQuad.push_back(totQuad[idx]);
779  fSigQuad.push_back(sigQuad[idx]);
780  fBkgQuad.push_back(bkgQuad[idx]);
781  }
782 
783  sampLab = "#nu_{#mu} Selection";
784  if(beam == "rhc") sampLab = "#bar{#nu}_{#mu} Selection";
785 
786  std::vector<int> colors = {kRed-7,kPink+6,kMagenta-6,kViolet-3};
787  new TCanvas;
788  myBarChart("Total Prediction Uncertainty (%)", labels,
789  fTot, fTotQuad, quantLab, colors);
790  gPad->SetLeftMargin(.35);
791  drawBeamLabel(beamLab);
792  drawSampleLabel(sampLab);
793  preliminary();
794 
795  gPad->Print(("plots/syst_quant_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_tot_"+extrapStr[i]+".pdf").c_str());
796 
797  new TCanvas;
798  myBarChart("Signal Uncertainty (%)", labels,
799  fSig, fSigQuad, quantLab, colors);
800  gPad->SetLeftMargin(.35);
801  drawBeamLabel(beamLab);
802  drawSampleLabel(sampLab);
803  preliminary();
804 
805  gPad->Print(("plots/syst_quant_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_sig_"+extrapStr[i]+".pdf").c_str());
806 
807  new TCanvas;
808  myBarChart("Background Uncertainty (%)", labels,
809  fBkg, fBkgQuad, quantLab, colors);
810  gPad->SetLeftMargin(.35);
811  drawBeamLabel(beamLab);
812  drawSampleLabel(sampLab);
813  preliminary();
814 
815  gPad->Print(("plots/syst_quant_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_bkg_"+extrapStr[i]+".pdf").c_str());
816  }
817  }
818 
819  // nue bins comparison
820  if(!IsNuMu && quantsComp)
821  {
822  unsigned int nBins = 3; // 4
823 
824  for(unsigned int i = 0; i < extrapStr.size(); ++i)
825  {
826  std::vector<std::string> quantLab = {"Low PID", "High PID", "Peripheral"};
827  std::string sumStr = "full";
828  std::vector<std::string> labels;
829  std::vector<std::vector<std::pair<double,double>>> fTot;
830  std::vector<std::vector<std::pair<double,double>>> fSig;
831  std::vector<std::vector<std::pair<double,double>>> fBkg;
832  std::vector<std::pair<double,double>> fTotQuad;
833  std::vector<std::pair<double,double>> fSigQuad;
834  std::vector<std::pair<double,double>> fBkgQuad;
835 
836  for(unsigned int b = 1; b <= nBins; ++b)
837  {
838  std::vector<std::pair<double,double>> tempTot;
839  std::vector<std::pair<double,double>> tempSig;
840  std::vector<std::pair<double,double>> tempBkg;
841  int idx = i*(nBins+1) + b;
842  std::cerr << "Bin idx: " << idx << " " << extrapStr[i] << std::endl;
843 
844  if(sum)
845  {
846  sumStr = "short";
847  for(int j = 0; j < kNumSystTypes; ++j){
848  if(b == 1) labels.push_back(systTypes[j]);
849  tempTot.push_back(tots[idx][systTypes[j]]);
850  tempSig.push_back(sigs[idx][systTypes[j]]);
851  tempBkg.push_back(bkgs[idx][systTypes[j]]);
852  }
853  }
854  else
855  {
856  for(const ISyst* syst: systs[0]){
857  if(b == 1) labels.push_back( (std::string)FixSystName(syst->ShortName()) );
858  tempTot.push_back(tots[idx][ (std::string)FixSystName(syst->ShortName()) ]);
859  tempSig.push_back(sigs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
860  tempBkg.push_back(bkgs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
861  }
862  }
863 
864 
865  fTot.push_back(tempTot);
866  fSig.push_back(tempSig);
867  fBkg.push_back(tempBkg);
868  fTotQuad.push_back(totQuad[idx]);
869  fSigQuad.push_back(sigQuad[idx]);
870  fBkgQuad.push_back(bkgQuad[idx]);
871  }
872 
873  sampLab = "#nu_{e} Selection";
874  if(beam == "rhc") sampLab = "#bar{#nu}_{e} Selection";
875 
876  std::vector<int> colors = {kViolet-4,kRed-7,kViolet+8};
877  new TCanvas;
878  myBarChart("Total Prediction Uncertainty (%)", labels,
879  fTot, fTotQuad, quantLab, colors);
880  gPad->SetLeftMargin(.35);
881  drawBeamLabel(beamLab);
882  drawSampleLabel(sampLab);
883  preliminary();
884 
885  gPad->Print(("plots/syst_bins_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_tot_"+extrapStr[i]+".pdf").c_str());
886 
887  new TCanvas;
888  myBarChart("Signal Uncertainty (%)", labels,
889  fSig, fSigQuad, quantLab, colors);
890  gPad->SetLeftMargin(.35);
891  drawBeamLabel(beamLab);
892  drawSampleLabel(sampLab);
893  preliminary();
894 
895  gPad->Print(("plots/syst_bins_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_sig_"+extrapStr[i]+".pdf").c_str());
896 
897  new TCanvas;
898  myBarChart("Background Uncertainty (%)", labels,
899  fBkg, fBkgQuad, quantLab, colors);
900  gPad->SetLeftMargin(.35);
901  drawBeamLabel(beamLab);
902  drawSampleLabel(sampLab);
903  preliminary();
904 
905  gPad->Print(("plots/syst_bins_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_bkg_"+extrapStr[i]+".pdf").c_str());
906  }
907  }
908 
909  if(extrapComp)
910  {
911  unsigned int nQuants = 1; //nue
912  if(!IsNuMu && quantsComp) nQuants = 4; // low, high PID, peripheral
913  if(IsNuMu) nQuants = 5; // 4 + summed
914 
915  for(unsigned int q = 0; q < nQuants; ++q)
916  {
917  std::string sumStr = "full";
918  std::vector<std::string> labels;
919  std::vector<std::vector<std::pair<double,double>>> fTot;
920  std::vector<std::vector<std::pair<double,double>>> fSig;
921  std::vector<std::vector<std::pair<double,double>>> fBkg;
922  std::vector<std::pair<double,double>> fTotQuad;
923  std::vector<std::pair<double,double>> fSigQuad;
924  std::vector<std::pair<double,double>> fBkgQuad;
925 
926  for(unsigned int i = 0; i < extrapStr.size(); ++i)
927  {
928  std::vector<std::pair<double,double>> tempTot;
929  std::vector<std::pair<double,double>> tempSig;
930  std::vector<std::pair<double,double>> tempBkg;
931  int idx = i*nQuants + q;
932  std::cerr << "Quartile idx: " << idx << " " << extrapStr[i] << std::endl;
933 
934  if(sum)
935  {
936  sumStr = "short";
937  for(int j = 0; j < kNumSystTypes; ++j){
938  if(i == 0) labels.push_back(systTypes[j]);
939  tempTot.push_back(tots[idx][systTypes[j]]);
940  tempSig.push_back(sigs[idx][systTypes[j]]);
941  tempBkg.push_back(bkgs[idx][systTypes[j]]);
942  }
943  }
944  else
945  {
946  // just for now hardcoded because the noextrap and normal extrap does not have Lepton Angle Syst and different accept syst
947  for(const ISyst* syst: systs[1]){
948  if(i == 0) labels.push_back( (std::string)FixSystName(syst->ShortName()) );
949  if ( i==0 && FixSystName(syst->ShortName()).Contains("LeptonAngle") )
950  {
951  tempTot.push_back({.0,.0});
952  tempSig.push_back({.0,.0});
953  tempBkg.push_back({.0,.0});
954  }
955  else if ( i==0 && FixSystName(syst->ShortName()).Contains("accept") )
956  {
957  if (beam == "fhc")
958  {
959  tempTot.push_back(tots[idx][ (std::string)"accept signalkin FHC 2020" ]);
960  tempSig.push_back(sigs[idx][ (std::string)"accept signalkin FHC 2020" ]);
961  tempBkg.push_back(bkgs[idx][ (std::string)"accept signalkin FHC 2020" ]);
962  }
963  else
964  {
965  tempTot.push_back(tots[idx][ (std::string)"accept signalkin RHC 2020" ]);
966  tempSig.push_back(sigs[idx][ (std::string)"accept signalkin RHC 2020" ]);
967  tempBkg.push_back(bkgs[idx][ (std::string)"accept signalkin RHC 2020" ]);
968  }
969  }
970  else
971  {
972  tempTot.push_back(tots[idx][ (std::string)FixSystName(syst->ShortName()) ]);
973  tempSig.push_back(sigs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
974  tempBkg.push_back(bkgs[idx][ (std::string)FixSystName(syst->ShortName()) ]);
975  }
976  }
977  }
978 
979  fTot.push_back(tempTot);
980  fSig.push_back(tempSig);
981  fBkg.push_back(tempBkg);
982  fTotQuad.push_back(totQuad[idx]);
983  fSigQuad.push_back(sigQuad[idx]);
984  fBkgQuad.push_back(bkgQuad[idx]);
985  }
986 
987  std::string sQuant = "";
988  if (IsNuMu) {
989  sQuant = "_Q"+std::to_string(q+1);
990  if(q==4) sQuant = "_AllQuartiles";
991  }
992  else {
993  sQuant = "_AllPrediction";
994  if(q==1) sQuant = "_lowPID";
995  if(q==2) sQuant = "_highPID";
996  if(q==3) sQuant = "_peripheral";
997  }
998 
999  sampLab = "#nu_{e} Selection";
1000  if(beam == "rhc") sampLab = "#bar{#nu}_{e} Selection";
1001  if(IsNuMu)
1002  {
1003  std::string quart = "All Quartiles";
1004  if(q<4) quart ="Quartile "+std::to_string(q+1);
1005  sampLab = "#nu_{#mu} Selection "+quart;
1006  if(beam == "rhc") sampLab = "#bar{#nu}_{#mu} Selection "+quart;
1007  }
1008  else
1009  {
1010  if(q==1) sampLab += " Low PID";
1011  if(q==2) sampLab += " High PID";
1012  if(q==3) sampLab += " Peripheral";
1013  }
1014 
1015  new TCanvas;
1016  myBarChart("Total Prediction Uncertainty (%)", labels,
1017  fTot, fTotQuad, extrapLab, {-1});
1018  gPad->SetLeftMargin(.35);
1019  drawBeamLabel(beamLab);
1020  drawSampleLabel(sampLab);
1021  preliminary();
1022 
1023  gPad->Print(("plots/syst_extrap_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_tot"+sQuant+".pdf").c_str());
1024 
1025  new TCanvas;
1026  myBarChart("Signal Uncertainty (%)", labels,
1027  fSig, fSigQuad, extrapLab, {-1});
1028  gPad->SetLeftMargin(.35);
1029  drawBeamLabel(beamLab);
1030  drawSampleLabel(sampLab);
1031  preliminary();
1032 
1033  gPad->Print(("plots/syst_extrap_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_sig"+sQuant+".pdf").c_str());
1034 
1035  new TCanvas;
1036  myBarChart("Background Uncertainty (%)", labels,
1037  fBkg, fBkgQuad, extrapLab, {-1});
1038  gPad->SetLeftMargin(.35);
1039  drawBeamLabel(beamLab);
1040  drawSampleLabel(sampLab);
1041  preliminary();
1042 
1043  gPad->Print(("plots/syst_extrap_comp_"+sNuMu+"_"+sumStr+"_"+beam+"_bkg"+sQuant+".pdf").c_str());
1044  }
1045  }
1046 
1047 }
1048 
1049 
1050 
1051 
1052 
1053 
1054 
1055 
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
void systematics_summary_from_pred_interp(std::string beam="rhc", std::string extrapStr="prop", bool sum=true)
int nBins
Definition: plotROC.py:16
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void drawBeamLabel(bool isFHC, double x=0.35)
double ssth23
set< int >::iterator it
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
std::map< std::string, std::vector< TH1 * > > SortSystHists(std::map< std::string, TH1 * > shifts)
(&#39; appearance&#39;)
Definition: IPrediction.h:18
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
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
(&#39;beam &#39;)
Definition: IPrediction.h:15
std::string systTypes[kNumSystTypes]
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual void SetTh13(const T &th13)=0
OStream cerr
Definition: OStream.cxx:7
std::map< std::string, std::pair< double, double > > SumSysts(std::map< std::string, std::pair< double, double >> bars)
Spectrum * ShiftedCosmics(Spectrum *cosmic, double sigma)
std::vector< const IPrediction * > GetNumuPredictions2020(const int nq=4, std::string decomp="noPt", osc::IOscCalc *calc=DefaultOscCalc(), bool useSysts=true, std::string beam="fhc", bool isFakeData=false, bool GetFromUPS=false, bool minimizeMemory=true, bool NERSC=false)
const Color_t kTotalMCErrorBandColor
Definition: Style.h:17
const CosmicBkgScaleSyst kCosmicBkgScaleSyst
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
virtual void SetDmsq32(const T &dmsq32)=0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::pair< Spectrum *, double > GetNueCosmics2020(std::string beam, bool GetFromUPS=false, bool NERSC=false)
void documentHeader(bool IsNuMu)
bool IsNuMu(int pdgc)
Definition: PDGUtils.cxx:155
Charged-current interactions.
Definition: IPrediction.h:39
int colors[6]
Definition: tools.h:1
Interactions of both types.
Definition: IPrediction.h:42
std::vector< TH1 * > MakeRatios(TH1 *nom, std::vector< TH1 * > shifts)
const int cols[3]
double dCP
Int_t col[ntarg]
Definition: Style.C:29
#define pot
const IPrediction * GetNuePrediction2020(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=true, bool NERSC=false)
const double kAna2020FHCLivetime
Definition: Exposures.h:237
const double j
Definition: BetheBloch.cxx:29
const double kAna2020FHCPOT
Definition: Exposures.h:233
const double kAna2020RHCPOT
Definition: Exposures.h:235
Neutrinos-only.
Definition: IPrediction.h:49
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
void myBarChart(std::string xlabel, std::vector< std::string > ylabels, std::vector< std::vector< std::pair< double, double >>> errs, std::vector< std::pair< double, double >> sum, std::vector< std::string > legend)
std::vector< std::pair< Spectrum *, double > > GetNumuCosmics2020(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
void drawSampleLabel(std::string s, double x=.915, double y=.9)
void MyPlotWithSystErrorBand(TH1 *&nom, std::vector< TH1 * > &ups, std::vector< TH1 * > &dns, int col, int errCol, float linewidth, bool newaxis)
double dmsq32
double ss2th13
double livetime
Definition: saveFDMCHists.C:21
double Integral(const Spectrum &s, double pot)
const double kAna2020RHCLivetime
Definition: Exposures.h:238
enum BeamMode kViolet
virtual void SetTh23(const T &th23)=0
TString FixSystName(TString mystring)
const hit & b
Definition: hits.cxx:21
Neutral-current interactions.
Definition: IPrediction.h:40
std::vector< const ISyst * > get3FlavorAna2020AllSysts(const EAnaType2020 ana, const bool smallgenie, const BeamType2020 beam, const bool isFit, const bool ptExtrap)
void ErrorBarChart(const std::map< std::string, std::pair< double, double >> &systErrors, const std::pair< double, double > &statErr, const std::string &label)
Make a simple plot of relative size of different errors.
Definition: Plots.cxx:798
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
int systCols[kNumSystTypes]
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
(&#39; appearance&#39;)
Definition: IPrediction.h:16
Double_t sum
Definition: plot.C:31
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
float fSig[xbins]
Definition: MakePlots.C:86
const DummyRockScaleSyst kRockScaleSyst
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string