make_surfprof_sensitivity.C
Go to the documentation of this file.
1 #pragma once
2 
5 
6 using namespace ana;
7 
8 std::string dataName = "fakedata";
9 
10 
12  std::string anaName,
13  bool nh)
14 {
15 
16  double sinHere = sin_rhcfhc;
17  double deltaHere = delta_rhcfhc;
18 
19  // naively changing hierarchy
20  // this isnt totally right
21  double deltaSign = 1.;
22  if(nh==false) deltaSign = -1.;
23 
24  if(anaName=="rhc"){
25  sinHere = sin_rhc;
26  deltaHere = delta_rhc*deltaSign;
27  }
28  if(anaName=="fhc"){
29  sinHere = sin_fhc;
30  deltaHere = delta_fhc*deltaSign;
31  }
33  calc->SetTh23(asin(sqrt(sinHere)));
34  calc->SetDmsq32(deltaHere);
35 }
36 
37 
38 void make_surfprof_sensitivity(const std::string anaName = "rhc",
39  const std::string extrapName = "extrap",
40  const int quantN = 0,
41  bool usecosmics = true,
42  bool systematics = true,
43  bool calc3a = true,
44  bool nh = true)
45 {
46 
47  if(anaName=="rhc"){
48  twobeams = false;
49  totquant = 4;
50  }
51  if(anaName=="fhc"){
52  horn_one = "fhc";
53  period_one = "full";
54  pot_one = pot_fhc;
56  twobeams = false;
57  totquant = 4;
58  }
59 
60 
61  // systematics
62  std::string systName = "";
63  std::vector< const ISyst* > systs;
64  systs = getAllAna2018Systs(kNumuAna2018,true);
65  if(systematics){
66  systName = "systs";
67  }
68  if(!systematics){
69  systName = "stats";
70  systs = {};
71  }
72  for (auto & sys:systs) std::cout << sys->ShortName() << std::endl;
73 
74 
75  // calculator
76  double sinmax = 0.7, sinmin = 0.3;
77  double deltamax = 3.0, deltamin = 2.0;
78  std::string calcName = "2018calc";
79  std::string hierarchy = "";
80  if(nh) hierarchy = "nh";
81  if(!nh) hierarchy = "ih";
82 
83 
84 
86  RestartCalculator(calc, anaName, nh);
87 
88  std::vector<const IFitVar*> fit_sin23delta32 = getNumuMarginalisedOscParam();
89  std::vector<const IFitVar*> fit_delta32 = getNumuMarginalisedOscParam();
90  std::vector<const IFitVar*> fit_sin23 = getNumuMarginalisedOscParam();
91  fit_sin23.push_back(&kFitSinSqTheta23);
92  fit_delta32.push_back(&kFitDmSq32Scaled);
93 
94 
95  // print on screen /////////////////////////////////////////////////////////
96  std::cout << "\n\nmake numu surface and profiles\n" << std::endl;
97  std::cout << "ana name: " << anaName << ", use cosmics: " << usecosmics << ", extrapolation: " << extrapName << std::endl;
98  std::cout << "systematics: " << systName << ", calculator: " << calcName << ", hierarchy: " << hierarchy << std::endl;
99 
100 
101  //inputs
102  TString dout_name = "../results/";
103  std::string ddata_name = "/nova/ana/nu_mu_ana/Ana2018/Data/";
104  std::string dcosm_name = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
105  std::string dpred_name = "/nova/ana/nu_mu_ana/Ana2018/Predictions/";
106  std::string fpred_name = "";
107  std::string fpred_suffix = "";
108 
109 
110  //predictions systs framework
111  const std::string dpred_name_fhc= dpred_name+"pred_numuconcat_fhc__numu2018.root";
112  const std::string dpred_name_rhc= dpred_name+"pred_numuconcat_rhc__numu2018.root";
113  std::string dpred_name_one = dpred_name_rhc;
114  std::string dpred_name_two = dpred_name_fhc;
115  if(anaName=="fhc") dpred_name_one = dpred_name_fhc;
116 
117 
118  std::cout << "\nloading predictions" << std::endl;
119  std::vector<PredictionSystJoint2018*> predictions;
120  ENu2018ExtrapType extrap = kNuMu;
121  for(int quant = 1; quant < totquant+1; quant++){
122  std::cout << "quantile " << quant << std::endl;
123  if(quant<5) predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_one, quant, systs, dpred_name_one));
124  else predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_two, quant-4, systs, dpred_name_two));
125  }
126 
127 
128  // predictions and data (real or fake)
129  std::cout << "\nfilling prediction" << std::endl;
130  std::vector<Spectrum> s_predictions;
131  std::vector<Spectrum> s_fakedata;
132  for(int quant = 0; quant < totquant; quant++){
133  std::cout << "quantile " << quant+1 << std::endl;
134  s_predictions.push_back(predictions[quant]->Predict(calc));
135  }
136  std::cout << "\nfilling fake data" << std::endl;
137  for(int quant = 0; quant < totquant; quant++){
138  std::cout << "quantile " << quant+1 << std::endl;
139  if(quant<4){ pot = pot_one; livetime = livetime_one;}
140  else {pot = pot_two; livetime = livetime_two;}
141  s_fakedata.push_back(s_predictions[quant].FakeData(pot));
142  }
143 
144 
145  // cosmics
146  std::string cosmics = "nocosmics";
147  std::vector<Spectrum> s_cosmics;
148  if(usecosmics){
150  cosmics = "cosmics";
151  std::cout << "\nloading cosmics" << std::endl;
152  std::cout << "" << horn_one << std::endl;
153  std::string fcosmics_name_one = dcosm_name + "cosmics_" + horn_one + "__numu2018.root";
154  TFile fcosmics_one(fcosmics_name_one.c_str());
155  auto *cosmics_quant1_one = (TH1*)fcosmics_one.Get("cosmics_q1");
156  auto *cosmics_quant2_one = (TH1*)fcosmics_one.Get("cosmics_q2");
157  auto *cosmics_quant3_one = (TH1*)fcosmics_one.Get("cosmics_q3");
158  auto *cosmics_quant4_one = (TH1*)fcosmics_one.Get("cosmics_q4");
159  s_cosmics.push_back(Spectrum (cosmics_quant1_one, pot, livetime));
160  s_cosmics.push_back(Spectrum (cosmics_quant2_one, pot, livetime));
161  s_cosmics.push_back(Spectrum (cosmics_quant3_one, pot, livetime));
162  s_cosmics.push_back(Spectrum (cosmics_quant4_one, pot, livetime));
163 
164  if(twobeams){
166  std::cout << "" << horn_two << std::endl;
167  std::string fcosmics_name_two = dcosm_name + "cosmics_" + horn_two + "__numu2018.root";
168  TFile fcosmics_two(fcosmics_name_two.c_str());
169  auto *cosmics_quant1_two = (TH1*)fcosmics_two.Get("cosmics_q1");
170  auto *cosmics_quant2_two = (TH1*)fcosmics_two.Get("cosmics_q2");
171  auto *cosmics_quant3_two = (TH1*)fcosmics_two.Get("cosmics_q3");
172  auto *cosmics_quant4_two = (TH1*)fcosmics_two.Get("cosmics_q4");
173  s_cosmics.push_back(Spectrum (cosmics_quant1_two, pot, livetime));
174  s_cosmics.push_back(Spectrum (cosmics_quant2_two, pot, livetime));
175  s_cosmics.push_back(Spectrum (cosmics_quant3_two, pot, livetime));
176  s_cosmics.push_back(Spectrum (cosmics_quant4_two, pot, livetime));
177  }
178  }
179 
180 
181 
182  std::vector <const IExperiment*> experiments;
183  if(quantN==0){
184  for(int quant = 0; quant < totquant; quant++){
185  if(usecosmics){
186  s_fakedata[quant] += s_cosmics[quant];
187  experiments.push_back(new SingleSampleExperiment(predictions[quant], s_fakedata[quant], s_cosmics[quant], 0.1));
188  }
189  else experiments.push_back(new SingleSampleExperiment(predictions[quant], s_fakedata[quant]));
190  }
191  }
192  else{
193  if(usecosmics){
194  s_fakedata[quantN-1] += s_cosmics[quantN-1];
195  experiments.push_back(new SingleSampleExperiment(predictions[quantN-1], s_fakedata[quantN-1], s_cosmics[quantN-1], 0.1));
196  }
197  else experiments.push_back(new SingleSampleExperiment(predictions[quantN-1], s_fakedata[quantN-1]));
198  }
199 
200  predictions.clear();
201  experiments.push_back(&kSolarConstraintsPDG2017);
202  experiments.push_back(WorldReactorConstraint2017());
203  MultiExperiment multiexpt(experiments);
204 
205 
206  TFile* fout;
207  if(quantN==0) fout = new TFile(dout_name + "surfprof_" + dataName + "_" + cosmics + "_" + anaName + "__" + extrapName + "_" + systName + "_" + calcName + "_" + hierarchy + "__numu2018.root","RECREATE");
208  else fout = new TFile(dout_name + "surfprof_" + dataName + "_" + cosmics + "_" + anaName + "__" + extrapName + "_" + systName + "_" + calcName + "_" + hierarchy + Form("__numu2018__q%d.root", quantN),"RECREATE");
209 
210  std::cout<< "save fake data histograms:" <<std::endl;
211  for(int quant = 0; quant < totquant; quant++){
212  std::cout << "quantile " << quant+1 << std::endl;
213  double pot_temp = s_fakedata[quant].POT();
214  TH1D* h_fakedata = s_fakedata[quant].ToTH1(pot_temp);
215  TH1D* h_prediction = predictions[quant]->Predict(calc).ToTH1(pot_temp);
216  h_fakedata -> Write(Form("fakedata_quant%d",quant+1));
217  h_prediction -> Write(Form("prediction_quant%d",quant+1));
218  }
219 
220 
221  RestartCalculator(calc, anaName, nh);
222  FrequentistSurface* surface = new FrequentistSurface(&multiexpt, calc, &kFitSinSqTheta23, 40, sinmin, sinmax, &kFitDmSq32Scaled, 50, deltamin, deltamax, fit_sin23delta32, systs);
223  RestartCalculator(calc, anaName, nh);
224  TH1* hprofile_delta32 = SqrtProfile(&multiexpt, calc, &kFitDmSq32Scaled, 50, deltamin, deltamax, -1, fit_sin23, systs);
225  RestartCalculator(calc, anaName, nh);
226  TH1* hprofile_sin23 = SqrtProfile(&multiexpt, calc, &kFitSinSqTheta23, 40, sinmin, sinmax, -1, fit_delta32, systs);
227 
228  surface->SaveTo(fout, "surface");
229  hprofile_sin23->Write("profile_sin23");
230  hprofile_delta32->Write("profile_delta32");
231 
232  experiments.clear();
233 
234 
235 
236 } // end make_surfprof
void make_surfprof_sensitivity(const std::string anaName="rhc", const std::string extrapName="extrap", const int quantN=0, bool usecosmics=true, bool systematics=true, bool calc3a=true, bool nh=true)
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
double pot_one
std::string period_one
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void FakeData()
Definition: rootlogon.C:156
void RestartCalculator(osc::IOscCalcAdjustable *calc, std::string anaName, bool nh)
Loads shifted spectra from files.
double livetime
std::string dataName
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double delta_rhcfhc
const double pot_fhc
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
virtual void SetDmsq32(const T &dmsq32)=0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::string extrapName
double livetime_two
Log-likelihood scan across two parameters.
double deltamin
Definition: plot_shifts.C:10
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
const double livetime_fhc
const double delta_fhc
const double sin_rhc
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
std::string horn_two
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
double livetime_one
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
void SaveTo(TDirectory *dir, const std::string &name) const
double pot
const double delta_rhc
double deltamax
Definition: plot_shifts.C:10
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
double pot_two
bool systematics
Definition: fnexvscaf.C:31
const double sin_rhcfhc
std::string horn_one
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
gm Write()
TGraph * SqrtProfile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, std::vector< const IFitVar * > profVars, std::vector< const ISyst * > profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &systsMap, MinuitFitter::FitOpts opts)
Forward to Profile but sqrt the result for a crude significance.
Definition: Fit.cxx:143
std::vector< const IFitVar * > getNumuMarginalisedOscParam()
Definition: FitVarsNumu.cxx:9
enum BeamMode string