plot_shifts.C
Go to the documentation of this file.
1 #pragma once
2 
6 
7 using namespace ana;
8 
9 std::string dataName = "realdata";
10 double deltamax = 3.0, deltamin = 2.0;
11 
12 
13 
15  std::string anaName,
16  bool nh)
17 {
18  double sinHere = sin_rhcfhc;
19  double deltaHere = delta_rhcfhc;
20  // naively changing hierarchy
21  // this isnt totally right
22  double deltaSign = 1.;
23  if(nh==false) deltaSign = -1.;
24 
25  if(anaName=="rhc"){
26  sinHere = sin_rhc;
27  deltaHere = delta_rhc*deltaSign;
28  }
29  if(anaName=="fhc"){
30  sinHere = sin_fhc;
31  deltaHere = delta_fhc*deltaSign;
32  }
34  calc->SetTh23(asin(sqrt(sinHere)));
35  calc->SetDmsq32(deltaHere);
36 }
37 
38 
39 void plot_shifts(const std::string anaName = "fhc",
40  const int quantN = 0,
41  bool usecosmics = true,
42  bool calc3a = true,
43  bool nh = true)
44 {
45 
46  std::string formalLabel = "";
47  if(anaName=="rhcfhc") formalLabel = "#nu_{#mu} + #bar{#nu}_{#mu}";
48  if(anaName=="fhc")formalLabel = "#nu_{#mu}";
49  if(anaName=="rhc")formalLabel = "#bar{#nu}_{#mu}";
50 
51  if(anaName=="rhc"){
52  twobeams = false;
53  totquant = 4;
54  }
55  if(anaName=="fhc"){
56  horn_one = "fhc";
57  period_one = "full";
58  pot_one = pot_fhc;
60  twobeams = false;
61  totquant = 4;
62  }
63 
64 
65  // systematics
66  std::string systName = "systs";
67  std::vector< const ISyst* > systs;
68  systs = getAllAna2018Systs(kNumuAna2018,true);
69  for (auto & sys:systs) std::cout << sys->ShortName() << std::endl;
70 
71 
72  // calculator
73  std::string calcName = "2018calc";
74  std::string hierarchy = "";
76  RestartCalculator(calc, anaName, nh);
77 
78 
79  // print on screen /////////////////////////////////////////////////////////
80  std::cout << "\n\nmake numu surface and profiles\n" << std::endl;
81  std::cout << "ana name: " << anaName << ", use cosmics: " << usecosmics << std::endl;
82  std::cout << "systematics: " << systName << ", calculator: " << calcName << ", hierarchy: " << hierarchy << std::endl;
83 
84 
85  //inputs
86  TString dout_name = "";
87  std::string ddata_name = "/nova/ana/nu_mu_ana/Ana2018/Data/";
88  std::string dcosm_name = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
89  std::string dpred_name = "/nova/ana/nu_mu_ana/Ana2018/Predictions/";
90  std::string fpred_name = "";
91  std::string fpred_suffix = "";
92 
93 
94  // this is not needed for sensitivity
95  // but using ana2107 selected fd data for testing
96  TFile* fdata_one = new TFile((ddata_name + "fd_dataspectrum_" + horn_one + "_full__numu2018.root").c_str());
97  TFile* fdata_two = new TFile((ddata_name + "fd_dataspectrum_" + horn_two + "_full__numu2018.root").c_str());
98  std::vector<Spectrum*> s_realdata;
99  std::cout << "\nloading data files and spectra" << std::endl;
100  for(int quant=1; quant<totquant+1; quant++){
101  if(quant<5) s_realdata.push_back(LoadFrom<Spectrum>(fdata_one, Form("fd_data_q%d", quant) ).release());
102  else s_realdata.push_back(LoadFrom<Spectrum>(fdata_two, Form("fd_data_q%d", quant-4) ).release());
103  }
104 
105 
106 
107  //predictions systs framework
108  const std::string dpred_name_fhc= dpred_name+"pred_numuconcat_fhc__numu2018.root";
109  const std::string dpred_name_rhc= dpred_name+"pred_numuconcat_rhc__numu2018.root";
110  std::string dpred_name_one = dpred_name_rhc;
111  std::string dpred_name_two = dpred_name_fhc;
112  if(anaName=="fhc") dpred_name_one = dpred_name_fhc;
113 
114  std::cout << "\nloading predictions" << std::endl;
115  std::vector<PredictionSystJoint2018*> predictions;
116  ENu2018ExtrapType extrap = kNuMu;
117  for(int quant = 1; quant < totquant+1; quant++){
118  std::cout << "quantile " << quant << std::endl;
119  if(quant<5) predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_one, quant, systs, dpred_name_one));
120  else predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_two, quant-4, systs, dpred_name_two));
121  }
122 
123 
124  // cosmics
125  std::string cosmics = "nocosmics";
126  std::vector<Spectrum> s_cosmics;
127  if(usecosmics){
129  cosmics = "cosmics";
130  std::cout << "\nloading cosmics" << std::endl;
131  std::cout << "" << horn_one << std::endl;
132  std::string fcosmics_name_one = dcosm_name + "cosmics_" + horn_one + "__numu2018.root";
133  TFile fcosmics_one(fcosmics_name_one.c_str());
134  auto *cosmics_quant1_one = (TH1*)fcosmics_one.Get("cosmics_q1");
135  auto *cosmics_quant2_one = (TH1*)fcosmics_one.Get("cosmics_q2");
136  auto *cosmics_quant3_one = (TH1*)fcosmics_one.Get("cosmics_q3");
137  auto *cosmics_quant4_one = (TH1*)fcosmics_one.Get("cosmics_q4");
138  s_cosmics.push_back(Spectrum (cosmics_quant1_one, pot, livetime));
139  s_cosmics.push_back(Spectrum (cosmics_quant2_one, pot, livetime));
140  s_cosmics.push_back(Spectrum (cosmics_quant3_one, pot, livetime));
141  s_cosmics.push_back(Spectrum (cosmics_quant4_one, pot, livetime));
142 
143  if(twobeams){
145  std::cout << "" << horn_two << std::endl;
146  std::string fcosmics_name_two = dcosm_name + "cosmics_" + horn_two + "__numu2018.root";
147  TFile fcosmics_two(fcosmics_name_two.c_str());
148  auto *cosmics_quant1_two = (TH1*)fcosmics_two.Get("cosmics_q1");
149  auto *cosmics_quant2_two = (TH1*)fcosmics_two.Get("cosmics_q2");
150  auto *cosmics_quant3_two = (TH1*)fcosmics_two.Get("cosmics_q3");
151  auto *cosmics_quant4_two = (TH1*)fcosmics_two.Get("cosmics_q4");
152  s_cosmics.push_back(Spectrum (cosmics_quant1_two, pot, livetime));
153  s_cosmics.push_back(Spectrum (cosmics_quant2_two, pot, livetime));
154  s_cosmics.push_back(Spectrum (cosmics_quant3_two, pot, livetime));
155  s_cosmics.push_back(Spectrum (cosmics_quant4_two, pot, livetime));
156  }
157  }
158 
159 
160  std::vector <const IExperiment*> experiments;
161  if(quantN==0){// all quantiles
162  for(int quant = 0; quant < totquant; quant++){
163  if(usecosmics){
164  experiments.push_back(new SingleSampleExperiment(predictions[quant], *s_realdata[quant], s_cosmics[quant], 0.1));
165  }
166  else experiments.push_back(new SingleSampleExperiment(predictions[quant], *s_realdata[quant]));
167  }
168  }
169  else{//specific quantile
170  if(usecosmics){
171  experiments.push_back(new SingleSampleExperiment(predictions[quantN-1], *s_realdata[quantN-1], s_cosmics[quantN-1], 0.1));
172  }
173  else experiments.push_back(new SingleSampleExperiment(predictions[quantN-1], *s_realdata[quantN-1]));
174  }
175 
176  predictions.clear();
177  experiments.push_back(&kSolarConstraintsPDG2017);
178  experiments.push_back(WorldReactorConstraint2017());
179  MultiExperiment multiexpt(experiments);
180 
181  TH1D* hMinus = new TH1D("hMinus","hMinus",100,-50,50);
182  TH1D* hPlus = new TH1D("hPlus","hPlus",100,-50,50);
183  TH1D* hZero = new TH1D("hZero","hZero",100,-50,50);
184  for(int k=0; k<100;k++){
185  hMinus->SetBinContent(k,-0.5);
186  hPlus->SetBinContent(k,0.5);
187  hZero->SetBinContent(k,0.0);
188  }
189  hMinus->SetLineStyle(2);
190  hPlus->SetLineStyle(2);
191  hMinus->SetLineWidth(1);
192  hPlus->SetLineWidth(1);
193  hZero->SetLineWidth(1);
194 
195  std::vector <const IFitVar*> fitvars = {&kFitSinSqTheta23, &kFitDmSq32Scaled,
198  &kFitDmSq21,
200 
201  MinuitFitter fit(&multiexpt, fitvars, systs);
202  double chi = fit.Fit(calc)->EvalMetricVal();
203  std::cout<< "\n\nLL: " << chi << "\n\n" << std::endl;
204 
205  // PLOT
206  // plot the systematic pulls and label
207  const SystShifts fitShifts = fit.GetSystShifts();
208  auto shifts = PlotSystShifts(fitShifts);
209  TString str = "Best fit " ;
210  for (auto &v:fitvars){
211  str += TString::Format(" %s=%.3f ",v->LatexName().c_str(),v->GetValue(calc));
212  }
213  str+= TString::Format(" LL=%.3f", chi);
214  shifts->SetTitle(str);
215  hPlus->Draw("hist same");
216  hMinus->Draw("hist same");
217  hZero->Draw("hist same");
218  CornerLabel(formalLabel);
219  if(quantN==0) gPad->Print(dout_name+"plotshifts_bestfit_" + anaName + ".png");
220  else gPad->Print(dout_name+"plotshifts_bestfit_"+anaName+Form("__q%d.png", quantN));
221  experiments.clear();
222 
223 
224 
225 
226 } // end make_surfprof
const FitSinSq2Theta12 kFitSinSq2Theta12
Definition: FitVars.cxx:25
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
TH1 * PlotSystShifts(const SystShifts &shifts, bool sortName)
Definition: Plots.cxx:1495
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
std::string dataName
Definition: plot_shifts.C:9
Loads shifted spectra from files.
double pot_two
Definition: saveFDMCHists.C:20
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
int totquant
Definition: saveFDMCHists.C:40
T sqrt(T number)
Definition: d0nt_math.hpp:156
void CornerLabel(std::string &s)
Definition: numu_tools.h:145
const double delta_rhcfhc
const double livetime_fhc
Definition: saveFDMCHists.C:12
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
double pot_one
Definition: saveFDMCHists.C:19
virtual void SetDmsq32(const T &dmsq32)=0
const FitDmSq21 kFitDmSq21
Definition: FitVars.cxx:26
bool twobeams
Definition: saveFDMCHists.C:39
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
double deltamin
Definition: plot_shifts.C:10
double livetime_two
Definition: saveFDMCHists.C:23
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
const double delta_fhc
const double sin_rhc
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:69
#define pot
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
const double sin_fhc
const double pot_fhc
Definition: saveFDMCHists.C:10
double livetime
Definition: saveFDMCHists.C:21
double livetime_one
Definition: saveFDMCHists.C:22
const double delta_rhc
double deltamax
Definition: plot_shifts.C:10
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
std::unique_ptr< SystShifts > GetSystShifts() const
Definition: IFitter.h:82
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
void plot_shifts(const std::string anaName="fhc", const int quantN=0, bool usecosmics=true, bool calc3a=true, bool nh=true)
Definition: plot_shifts.C:39
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
const double sin_rhcfhc
void RestartCalculator(osc::IOscCalcAdjustable *calc, std::string anaName, bool nh)
Definition: plot_shifts.C:14
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::string period_one
Definition: saveFDMCHists.C:35
std::string horn_one
Definition: saveFDMCHists.C:33
std::string horn_two
Definition: saveFDMCHists.C:34
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string