check_predinterp_numu.C
Go to the documentation of this file.
3 #include "CAFAna/Fit/Fit.h"
11 #include "CAFAna/FC/FCSurface.h"
14 #include "CAFAna/Vars/FitVars.h"
15 
16 #include "Utilities/rootlogon.C"
17 #include "OscLib/IOscCalc.h"
18 
19 #include "TCanvas.h"
20 #include "TBox.h"
21 #include "TColor.h"
22 #include "TGraph.h"
23 #include "TVectorD.h"
24 #include "TF1.h"
25 #include "TLegend.h"
26 #include "TLine.h"
27 #include "TMarker.h"
28 #include "TStyle.h"
29 #include "TSystem.h"
30 #include "TGaxis.h"
31 
32 #include <algorithm>
33 #include <string>
34 
35 using namespace ana;
36 
37 void check_predinterp_numu(std::string decomp = "noextrap", int WhichQuant=1, bool NewBestFit = false)
38 {
39 
40  std::cout << "\n\n NewBestFit is " << NewBestFit << "\n\n" << std::endl;
41 
42  TString outBase = "/nova/ana/users/karlwarb/Systematic_Interp/";
43 // if (NewBestFit) outBase = "/nova/ana/users/karlwarb/Systematic_Interp_NewBestFit/";
44  if (NewBestFit) outBase = "/nova/ana/users/ecatanom/Ana2017/numu_predinterp_debug/"+decomp + "/";
45  TString outdir = outBase + decomp + "_Quant"+std::to_string(WhichQuant)+"/";
46  gSystem -> MakeDirectory(outBase);
47  gSystem -> MakeDirectory(outdir );
48 
50  if(decomp == "numu") extrap = kNuMu;
51 
52  auto systs = kNumu2017Systs;
53  for (auto s:systs){
54  std::cout << s->ShortName() << ",";
55  }
57 
58  auto calc = DefaultOscCalc();
59  if (NewBestFit) {
60  calc->SetDmsq21(7.53011e-05);
61  calc->SetTh12 (asin( sqrt( 0.851012) ) );
62  calc->SetTh23 (asin( sqrt( 0.471974) ) );
63  calc->SetDmsq32(2.43071e-3 );
64  calc->SetdCP (M_PI*0.915869);
65  }
66  PredictionSystNumu2017 pred(extrap, calc, WhichQuant);
67  pred.DebugPlots(calc, (outdir+"plot_debug_predinterp_numu2017_"+decomp+"_Quant"+std::to_string(WhichQuant)+"_%s.pdf").Data());
68  //pred.DebugPlots(calc, "debug/plot_debug_predinterp_nueSA_sig_Quant"+std::to_string(WhichQuant)+"_%s.pdf" ,Flavors::kNuMuToNuE,Current::kCC);
69  //pred.DebugPlots(calc, "debug/plot_debug_predinterp_nueSA_beam_Quant"+std::to_string(WhichQuant)+"_%s.pdf",Flavors::kNuEToNuE ,Current::kCC);
70  //pred.DebugPlots(calc, "debug/plot_debug_predinterp_nueSA_nc_Quant"+std::to_string(WhichQuant)+"_%s.pdf" ,Flavors::kAll ,Current::kNC);
71  //pred.DebugPlots(calc, "debug/plot_debug_predinterp_nueSA_numu_Quant"+std::to_string(WhichQuant)+"_%s.pdf",Flavors::kAllNuMu ,Current::kCC);
72 
73  struct prettyshift{
74  int sigma;
75  TString name;
76  int color;
77  int style;
78  };
79 
80  std::vector <prettyshift> prettyshifts= {{+3,"+3 #sigma",kYellow+1,3},
81  {-3,"-3 #sigma",kGreen+1,3},
82  {+2,"+2 #sigma",kOrange,9},
83  {-2,"-2 #sigma",kAzure+10,7},
84  {+1,"+1 #sigma",kRed,1},
85  {-1,"-1 #sigma",kAzure,1},
86  {0,"Nominal " ,kBlack,1}};
87 
88  double fdpot = kSecondAnaPOT;
89  for (int i=0; i<1/*5*/; ++i){
90  std::string pname ="";
91  TString plotname = "Numu 2017 - total prediction";
92  switch(i){
93  case 1: pname = "_sig" ; plotname = "Numu 2017 - signal " ;break;
94  case 2: pname = "_beam"; plotname = "Numu 2017 - Beam Nue CC"; break;
95  case 3: pname = "_nc" ; plotname = "Numu 2017 - NC" ; break;
96  case 4: pname = "_numu"; plotname = "Numu 2017 - Numu App" ; break;
97  }
98  for (auto const & syst:systs){
99  new TCanvas;
100  auto leg = new TLegend(0.1,0.65,0.5,0.89);
101  leg->SetNColumns(2);
102  for(auto const & sh:prettyshifts){
103  TH1 * h1;
104  SystShifts shift(syst,sh.sigma);
105  if (i==0) h1 = pred.PredictSyst(calc,shift).ToTH1(fdpot,sh.color, sh.style);
106  if (i==1) h1 = pred.PredictComponentSyst(calc,shift,Flavors::kNuMuToNuE,Current::kCC,Sign::kBoth).ToTH1(fdpot,sh.color, sh.style);
107  if (i==2) h1 = pred.PredictComponentSyst(calc,shift,Flavors::kNuEToNuE ,Current::kCC,Sign::kBoth).ToTH1(fdpot,sh.color, sh.style);
108  if (i==3) h1 = pred.PredictComponentSyst(calc,shift,Flavors::kAll ,Current::kNC,Sign::kBoth).ToTH1(fdpot,sh.color, sh.style);
109  if (i==4) h1 = pred.PredictComponentSyst(calc,shift,Flavors::kAllNuMu ,Current::kCC,Sign::kBoth).ToTH1(fdpot,sh.color, sh.style);
110  h1->SetTitle(plotname);
111  h1->Scale(0.1,"width");
112  h1->GetYaxis()->SetTitle("Events / 0.1 GeV");
113  h1->GetYaxis()->SetRangeUser(0,2);
114  h1->Draw("hist same");
115  leg->AddEntry(h1,sh.name,"l");
116  }
117  leg->SetHeader((syst->LatexName()).c_str());
118  leg->SetFillStyle(0);
119  leg->SetTextSize(0.06);
120  leg->Draw();
121  gPad->Print((outdir+"plot_debug_predinterp_numu2017_"+decomp+"_Quant"+std::to_string(WhichQuant)+"_spec"+pname+"_"+syst->ShortName()+"_calc.pdf"));
122  }
123  }
124 }
const XML_Char * name
Definition: expat.h:151
enum BeamMode kOrange
void check_predinterp_numu(std::string decomp="noextrap", int WhichQuant=1, bool NewBestFit=false)
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
(&#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
T sqrt(T number)
Definition: d0nt_math.hpp:156
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
#define M_PI
Definition: SbMath.h:34
const std::vector< const DummyNumu2017Syst * > kNumu2017Systs
const XML_Char * s
Definition: expat.h:262
Charged-current interactions.
Definition: IPrediction.h:39
Loads shifted spectra from files.
Spectrum PredictComponentSyst(osc::IOscCalc *calc, const SystShifts &shift, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const override
string pname
Definition: eplot.py:33
double sigma(TH1F *hist, double percentile)
void DebugPlots(osc::IOscCalc *calc, const std::string &savePattern="", Flavors::Flavors_t flav=Flavors::kAll, Current::Current_t curr=Current::kBoth, Sign::Sign_t sign=Sign::kBoth) const
OStream cout
Definition: OStream.cxx:6
TH1F * h1
Neutral-current interactions.
Definition: IPrediction.h:40
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t e
Definition: plot.C:35
const double kSecondAnaPOT
Definition: Exposures.h:87
enum BeamMode kGreen
const std::string outdir
T asin(T number)
Definition: d0nt_math.hpp:60
enum BeamMode string