Functions | Variables
make_surfprof_sensitivity.C File Reference
#include "includes.h"

Go to the source code of this file.

Functions

void make_surfprof_sensitivity (const std::string anaName="fhc", const std::string extrapName="extrap", bool usecosmics=false, bool nh=true, bool systematics=true, bool calc3a=true)
 

Variables

std::string dataName = "fakedata"
 
std::string horn_one = "rhc"
 
std::string horn_two = "fhc"
 
std::string period_one = "full"
 
std::string period_two = "full"
 
const double pot_fhc = kAna2018FHCPOT
 
const double pot_rhc = kAna2018RHCPOT
 
const double livetime_fhc = kAna2018FHCLivetime
 
const double livetime_rhc = kAna2018RHCLivetime
 
double pot = 0
 
double pot_one = pot_rhc
 
double pot_two = pot_fhc
 
double livetime = 0
 
double livetime_one = livetime_rhc
 
double livetime_two = livetime_fhc
 
int totquant = 8
 
bool twobeams = true
 

Function Documentation

void make_surfprof_sensitivity ( const std::string  anaName = "fhc",
const std::string  extrapName = "extrap",
bool  usecosmics = false,
bool  nh = true,
bool  systematics = true,
bool  calc3a = true 
)

Definition at line 33 of file make_surfprof_sensitivity.C.

References std::asin(), calc, make_associated_cosmic_defs::cosmics, om::cout, dataName, ana::DefaultOscCalc(), deltamax, deltamin, allTimeWatchdog::endl, extrapName, FakeData(), plot_validation_datamc::fdata, submit_syst::fout, ana::getAllAna2018Systs(), ana::getNumuMarginalisedOscParam(), horn_one, horn_two, ana::kFitDmSq32Scaled, ana::kFitSinSqTheta23, ana::kNuMu, ana::kNumuAna2018, ana::kSolarConstraintsPDG2017, livetime, livetime_fhc, livetime_one, livetime_two, period_one, pot, pot_fhc, pot_one, pot_two, ana::FrequentistSurface::SaveTo(), osc::_IOscCalcAdjustable< T >::SetDmsq32(), osc::_IOscCalcAdjustable< T >::SetTh23(), std::sqrt(), ana::SqrtProfile(), string, systs, totquant, twobeams, ana::WorldReactorConstraint2017(), and Write().

39 {
40 
41  if(anaName=="rhc"){
42  twobeams = false;
43  totquant = 4;
44  }
45  if(anaName=="fhc"){
46  horn_one = "fhc";
47  period_one = "full";
48  pot_one = pot_fhc;
50  twobeams = false;
51  totquant = 4;
52  }
53 
54 
55  // systematics
56  std::string systName = "";
58  if(systematics){
59  systName = "systs";
60  }
61  if(!systematics){
62  systName = "stats";
63  systs = {};
64  }
65  for (auto & sys:systs) std::cout << sys->ShortName() << std::endl;
66 
67 
68  // calculator
69  double sinmax = 0.7, sinmin = 0.3;
70  double deltamax = 3.0, deltamin = 2.0;
71  std::string calcName = "2018calc";
72  std::string hierarchy = "";
74  if(nh){
75  hierarchy = "nh";
76  calc->SetTh23(asin(sqrt(0.558)));
77  calc->SetDmsq32(0.00244);
78  }
79  if(!nh){
80  hierarchy = "ih";
81  calc->SetTh23(asin(sqrt(0.558)));
82  calc->SetDmsq32(-0.00244);
83  deltamax = -2.0; deltamin = -3.0;
84  }
85 
86 
87  std::vector<const IFitVar*> fit_sin23delta32 = getNumuMarginalisedOscParam();
88  std::vector<const IFitVar*> fit_delta32 = getNumuMarginalisedOscParam();
89  std::vector<const IFitVar*> fit_sin23 = getNumuMarginalisedOscParam();
90  fit_sin23.push_back(&kFitSinSqTheta23);
91  fit_delta32.push_back(&kFitDmSq32Scaled);
92 
93 
94  // print on screen /////////////////////////////////////////////////////////
95  std::cout << "\n\nmake numu surface and profiles\n" << std::endl;
96  std::cout << "ana name: " << anaName << ", use cosmics: " << usecosmics << ", extrapolation: " << extrapName << std::endl;
97  std::cout << "systematics: " << systName << ", calculator: " << calcName << ", hierarchy: " << hierarchy << std::endl;
98 
99 
100  //inputs
101  TString dout_name = "results/";
102  std::string ddata_name = "data/";
103  std::string dcosm_name = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
104  std::string dpred_name = "/nova/ana/nu_mu_ana/Ana2018/Predictions/";
105  std::string fpred_name = "";
106  std::string fpred_suffix = "";
107 
108 
109  // this is not needed for sensitivity
110  // but using ana2107 selected fd data for testing
111  std::vector<TFile*> fdata;
112  std::cout << "\nloading data files" << std::endl;
113  std::cout << "" << horn_one << std::endl;
114  fdata.push_back(new TFile( (ddata_name + "fd_data_quant1_" + horn_one + ".root").c_str() ));
115  fdata.push_back(new TFile( (ddata_name + "fd_data_quant2_" + horn_one + ".root").c_str() ));
116  fdata.push_back(new TFile( (ddata_name + "fd_data_quant3_" + horn_one + ".root").c_str() ));
117  fdata.push_back(new TFile( (ddata_name + "fd_data_quant3_" + horn_one + ".root").c_str() ));
118  if(twobeams){
119  std::cout << "" << horn_two << std::endl;
120  fdata.push_back(new TFile( (ddata_name + "fd_data_quant1_" + horn_two + ".root").c_str() ));
121  fdata.push_back(new TFile( (ddata_name + "fd_data_quant2_" + horn_two + ".root").c_str() ));
122  fdata.push_back(new TFile( (ddata_name + "fd_data_quant3_" + horn_two + ".root").c_str() ));
123  fdata.push_back(new TFile( (ddata_name + "fd_data_quant3_" + horn_two + ".root").c_str() ));
124  }
125 
126  // const int totquant = fdata.size();
127 
128 
129  std::cout << "\nloading predictions" << std::endl;
130  std::vector<PredictionInterp*> predictions;
131  ENumu2018ExtrapType extrap = kNuMu;
132 
133  //predictions systs framework
134  const std::string dpred_name_fhc= dpred_name+"pred_numuconcat_fhc__numu2018.root";
135  const std::string dpred_name_rhc= dpred_name+"pred_numuconcat_rhc__numu2018.root";
136  std::string dpred_name_one = dpred_name_rhc;
137  std::string dpred_name_two = dpred_name_fhc;
138  if(anaName=="fhc") dpred_name_one = dpred_name_fhc;
139 
140  for(int quant = 1; quant < totquant+1; quant++){
141  std::cout << "quantile " << quant << std::endl;
142  if(quant<5) predictions.push_back(new PredictionSystNumu2018(extrap, calc, quant, kNumu2018Systs, dpred_name_one));
143  else predictions.push_back(new PredictionSystNumu2018(extrap, calc, quant-4, kNumu2018Systs, dpred_name_two));
144  }
145 
146 
147  /*
148  //predictions old style
149  std::cout << "\nloading predictions" << std::endl;
150  for(int quant = 1; quant < totquant+1; quant++){
151  std::cout << "quantile " << quant << std::endl;
152  if(quant<5){
153  fpred_suffix = extrapName + "_" + systName + "__" + horn_one + "_" + calcName + "_" + hierarchy + "__" + period_one + "__numu2018.root";
154  fpred_name = Form("prediction_quant%d__", quant);
155  std::string input = dpred_name + fpred_name + fpred_suffix;
156  TFile* file = new TFile(input.c_str());
157  predictions.push_back(LoadFrom<PredictionInterp>(file, "prediction")).release();
158  file->Close();
159  }
160  else{
161  fpred_suffix = extrapName + "_" + systName + "__" + horn_two + "_" + calcName + "_" + hierarchy + "__" + period_two + "__numu2018.root";
162  fpred_name = Form("prediction_quant%d__", quant-4);
163  std::string input = dpred_name + fpred_name + fpred_suffix;
164  TFile* file = new TFile(input.c_str());
165  predictions.push_back(LoadFrom<PredictionInterp>(file, "prediction")).release();
166  file->Close();
167  }
168  }//loop quant
169  */
170 
171 
172 
173  // predictions and data (real or fake)
174  std::cout << "\nfilling prediction" << std::endl;
175  std::vector<Spectrum> s_predictions;
176  std::vector<Spectrum> s_fakedata;
177  for(int quant = 0; quant < totquant; quant++){
178  std::cout << "quantile " << quant+1 << std::endl;
179  s_predictions.push_back(predictions[quant]->Predict(calc));
180  }
181  std::cout << "\nfilling fake data" << std::endl;
182  for(int quant = 0; quant < totquant; quant++){
183  std::cout << "quantile " << quant+1 << std::endl;
184  if(quant<4){ pot = pot_one; livetime = livetime_one;}
185  else {pot = pot_two; livetime = livetime_two;}
186  s_fakedata.push_back(s_predictions[quant].FakeData(pot));
187  }
188 
189 
190  // cosmics
191  std::string cosmics = "nocosmics";
192  std::vector<Spectrum> s_cosmics;
193  if(usecosmics){
195  cosmics = "cosmics";
196  std::cout << "\nloading cosmics" << std::endl;
197  std::cout << "" << horn_one << std::endl;
198  std::string fcosmics_name_one = dcosm_name + "cosmics_" + horn_one + "__numu2018.root";
199  TFile fcosmics_one(fcosmics_name_one.c_str());
200  auto *cosmics_quant1_one = (TH1*)fcosmics_one.Get("cosmics_q1");
201  auto *cosmics_quant2_one = (TH1*)fcosmics_one.Get("cosmics_q2");
202  auto *cosmics_quant3_one = (TH1*)fcosmics_one.Get("cosmics_q3");
203  auto *cosmics_quant4_one = (TH1*)fcosmics_one.Get("cosmics_q4");
204  s_cosmics.push_back(Spectrum (cosmics_quant1_one, pot, livetime));
205  s_cosmics.push_back(Spectrum (cosmics_quant2_one, pot, livetime));
206  s_cosmics.push_back(Spectrum (cosmics_quant3_one, pot, livetime));
207  s_cosmics.push_back(Spectrum (cosmics_quant4_one, pot, livetime));
208 
209  if(twobeams){
211  std::cout << "" << horn_two << std::endl;
212  std::string fcosmics_name_two = dcosm_name + "cosmics_" + horn_two + "__numu2018.root";
213  TFile fcosmics_two(fcosmics_name_two.c_str());
214  auto *cosmics_quant1_two = (TH1*)fcosmics_two.Get("cosmics_q1");
215  auto *cosmics_quant2_two = (TH1*)fcosmics_two.Get("cosmics_q2");
216  auto *cosmics_quant3_two = (TH1*)fcosmics_two.Get("cosmics_q3");
217  auto *cosmics_quant4_two = (TH1*)fcosmics_two.Get("cosmics_q4");
218  s_cosmics.push_back(Spectrum (cosmics_quant1_two, pot, livetime));
219  s_cosmics.push_back(Spectrum (cosmics_quant2_two, pot, livetime));
220  s_cosmics.push_back(Spectrum (cosmics_quant3_two, pot, livetime));
221  s_cosmics.push_back(Spectrum (cosmics_quant4_two, pot, livetime));
222  }
223  }
224 
225 
226 
227  std::vector <const IExperiment*> experiments;
228  for(int quant = 0; quant < totquant; quant++){
229  if(usecosmics){
230  s_fakedata[quant] += s_cosmics[quant];
231  experiments.push_back(new SingleSampleExperiment(predictions[quant], s_fakedata[quant], s_cosmics[quant], 0.1));
232  }
233  else experiments.push_back(new SingleSampleExperiment(predictions[quant], s_fakedata[quant]));
234  }
235 
236 
237  predictions.clear();
238  experiments.push_back(&kSolarConstraintsPDG2017);
239  experiments.push_back(WorldReactorConstraint2017());
240  MultiExperiment multiexpt(experiments);
241 
242 
243  TFile* fout = new TFile(dout_name + "surfprof_" + dataName + "_" + cosmics + "_" + anaName + "__" + extrapName + "_" + systName + "_" + calcName + "_" + hierarchy + "__numu2018.root","RECREATE");
244  std::cout<< "save fake data histograms:" <<std::endl;
245  for(int quant = 0; quant < totquant; quant++){
246  std::cout << "quantile " << quant+1 << std::endl;
247  double pot_temp = s_fakedata[quant].POT();
248  TH1D* h_fakedata = s_fakedata[quant].ToTH1(pot_temp);
249  TH1D* h_prediction = predictions[quant]->Predict(calc).ToTH1(pot_temp);
250  h_fakedata -> Write(Form("fakedata_quant%d",quant+1));
251  h_prediction -> Write(Form("prediction_quant%d",quant+1));
252  }
253 
254 
255  FrequentistSurface* surface = new FrequentistSurface(&multiexpt, calc, &kFitSinSqTheta23, 40, sinmin, sinmax, &kFitDmSq32Scaled, 50, deltamin, deltamax, fit_sin23delta32, systs);
256  TH1* hprofile_delta32 = SqrtProfile(&multiexpt, calc, &kFitDmSq32Scaled, 50, deltamin, deltamax, -1, fit_sin23, systs);
257  TH1* hprofile_sin23 = SqrtProfile(&multiexpt, calc, &kFitSinSqTheta23, 40, sinmin, sinmax, -1, fit_delta32, systs);
258 
259  surface->SaveTo(fout, "surface");
260  hprofile_sin23->Write("profile_sin23");
261  hprofile_delta32->Write("profile_delta32");
262 
263  experiments.clear();
264 
265 
266 
267 } // end make_surfprof
double pot_one
std::string period_one
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
void FakeData()
Definition: rootlogon.C:156
std::string dataName
double livetime
T sqrt(T number)
Definition: d0nt_math.hpp:156
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 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)
void SaveTo(TDirectory *dir, const std::string &name) const
double pot
double deltamax
Definition: plot_shifts.C:10
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
double pot_two
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

Variable Documentation

std::string dataName = "fakedata"

Definition at line 9 of file make_surfprof_sensitivity.C.

std::string horn_one = "rhc"

Definition at line 11 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

std::string horn_two = "fhc"

Definition at line 12 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

double livetime = 0

Definition at line 25 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

const double livetime_fhc = kAna2018FHCLivetime

Definition at line 19 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

double livetime_one = livetime_rhc

Definition at line 26 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

const double livetime_rhc = kAna2018RHCLivetime

Definition at line 20 of file make_surfprof_sensitivity.C.

double livetime_two = livetime_fhc

Definition at line 27 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

std::string period_one = "full"

Definition at line 13 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

std::string period_two = "full"

Definition at line 14 of file make_surfprof_sensitivity.C.

double pot = 0

Definition at line 22 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

const double pot_fhc = kAna2018FHCPOT

Definition at line 17 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

double pot_one = pot_rhc

Definition at line 23 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

const double pot_rhc = kAna2018RHCPOT

Definition at line 18 of file make_surfprof_sensitivity.C.

double pot_two = pot_fhc

Definition at line 24 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

int totquant = 8

Definition at line 29 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().

bool twobeams = true

Definition at line 30 of file make_surfprof_sensitivity.C.

Referenced by make_surfprof_sensitivity().