make_surfprof.C
Go to the documentation of this file.
1 #pragma once
2 
5 
6 using namespace ana;
7 
8 
9 std::string dataName = "realdata";
10 
11 
13 {
15  calc->SetTh23(asin(sqrt(0.558)));
16  calc->SetDmsq32(0.00244);
17 }
18 
19 
20 void make_surfprof(const std::string anaName = "rhcfhc",
21  const std::string extrapName = "extrap",
22  const int quantN = 0,
23  bool usecosmics = true,
24  bool systematics = false,
25  bool calc3a = true,
26  bool nh = true)
27 {
28 
29  if(anaName=="rhc"){
30  twobeams = false;
31  totquant = 4;
32  }
33  if(anaName=="fhc"){
34  horn_one = "fhc";
35  period_one = "full";
36  pot_one = pot_fhc;
38  twobeams = false;
39  totquant = 4;
40  }
41 
42 
43  // systematics
44  std::string systName = "";
45  std::vector< const ISyst* > systs;
46  systs = getAllAna2018Systs(kNumuAna2018,true);
47  if(systematics){
48  systName = "systs";
49  }
50  if(!systematics){
51  systName = "stats";
52  systs = {};
53  }
54  for (auto & sys:systs) std::cout << sys->ShortName() << std::endl;
55 
56 
57  // calculator
58  double sinmax = 0.7, sinmin = 0.3;
59  double deltamax = 3.0, deltamin = 2.0;
60  std::string calcName = "2018calc";
61  std::string hierarchy = "";
63  if(nh){
64  hierarchy = "nh";
65  calc->SetTh23(asin(sqrt(0.558)));
66  calc->SetDmsq32(0.00244);
67  }
68  if(!nh){
69  hierarchy = "ih";
70  calc->SetTh23(asin(sqrt(0.558)));
71  calc->SetDmsq32(-0.00244);
72  deltamax = -2.0; deltamin = -3.0;
73  }
74 
75  std::vector<const IFitVar*> fit_sin23delta32 = getNumuMarginalisedOscParam();
76  std::vector<const IFitVar*> fit_delta32 = getNumuMarginalisedOscParam();
77  std::vector<const IFitVar*> fit_sin23 = getNumuMarginalisedOscParam();
78  fit_sin23.push_back(&kFitSinSqTheta23);
79  fit_delta32.push_back(&kFitDmSq32Scaled);
80 
81 
82  // print on screen /////////////////////////////////////////////////////////
83  std::cout << "\n\nmake numu surface and profiles\n" << std::endl;
84  std::cout << "ana name: " << anaName << ", use cosmics: " << usecosmics << ", extrapolation: " << extrapName << std::endl;
85  std::cout << "systematics: " << systName << ", calculator: " << calcName << ", hierarchy: " << hierarchy << std::endl;
86 
87 
88  //inputs
89  //TString dout_name = "../results/";
90  TString dout_name = "";
91  std::string ddata_name = "/nova/ana/nu_mu_ana/Ana2018/Data/";
92  std::string dcosm_name = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
93  //std::string dpred_name = "/nova/ana/nu_mu_ana/Ana2018/Predictions/";
94  std::string dpred_name = "/pnfs/nova/scratch/users/dmendez/Numu/RPAFix_noDIS/";
95  std::string fpred_name = "";
96  std::string fpred_suffix = "";
97 
98 
99  // this is not needed for sensitivity
100  // but using ana2107 selected fd data for testing
101  TFile* fdata_one = new TFile((ddata_name + "fd_dataspectrum_" + horn_one + "_full__numu2018.root").c_str());
102  TFile* fdata_two = new TFile((ddata_name + "fd_dataspectrum_" + horn_two + "_full__numu2018.root").c_str());
103  std::vector<Spectrum*> s_realdata;
104  std::cout << "\nloading data files and spectra" << std::endl;
105  for(int quant=1; quant<totquant+1; quant++){
106  if(quant<5) s_realdata.push_back(LoadFrom<Spectrum>(fdata_one, Form("fd_data_q%d", quant) ).release());
107  else s_realdata.push_back(LoadFrom<Spectrum>(fdata_two, Form("fd_data_q%d", quant-4) ).release());
108  }
109 
110 
111  //predictions systs framework
112  const std::string dpred_name_fhc= dpred_name+"pred_numuconcat_fhc__numu2018.root";
113  const std::string dpred_name_rhc= dpred_name+"pred_numuconcat_rhc__numu2018.root";
114  std::string dpred_name_one = dpred_name_rhc;
115  std::string dpred_name_two = dpred_name_fhc;
116  if(anaName=="fhc") dpred_name_one = dpred_name_fhc;
117 
118 
119  /*
120  std::cout << "\nloading predictions" << std::endl;
121  std::vector<PredictionSystJoint2018*> predictions;
122  ENu2018ExtrapType extrap = kNuMu;
123  for(int quant = 1; quant < totquant+1; quant++){
124  std::cout << "quantile " << quant << std::endl;
125  if(quant<5) predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_one, quant, systs, dpred_name_one));
126  else predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_two, quant-4, systs, dpred_name_two));
127  }
128  */
129 
130  //predictions old style
131  std::vector<PredictionInterp*> predictions;
132  std::cout << "\nloading predictions" << std::endl;
133  for(int quant = 1; quant < totquant+1; quant++){
134  std::cout << "quantile " << quant << std::endl;
135  if(quant<5){
136  fpred_suffix = extrapName + "_" + systName + "__" + horn_one + "_" + calcName + "_" + hierarchy + "__" + period_one + "__numu2018.root";
137  fpred_name = Form("prediction_quant%d__", quant);
138  std::string input = dpred_name + fpred_name + fpred_suffix;
139  TFile* file = new TFile(input.c_str());
140  predictions.push_back(LoadFrom<PredictionInterp>(file, "prediction")).release();
141  file->Close();
142  }
143  else{
144  fpred_suffix = extrapName + "_" + systName + "__" + horn_two + "_" + calcName + "_" + hierarchy + "__" + period_two + "__numu2018.root";
145  fpred_name = Form("prediction_quant%d__", quant-4);
146  std::string input = dpred_name + fpred_name + fpred_suffix;
147  TFile* file = new TFile(input.c_str());
148  predictions.push_back(LoadFrom<PredictionInterp>(file, "prediction")).release();
149  file->Close();
150  }
151  }//loop quant
152 
153 
154 
155  /*
156  // predictions and data (real or fake)
157  std::cout << "\nfilling prediction" << std::endl;
158  std::vector<Spectrum> s_predictions;
159  for(int quant = 0; quant < totquant; quant++){
160  std::cout << "quantile " << quant+1 << std::endl;
161  s_predictions.push_back(predictions[quant]->Predict(calc));
162  }
163  */
164 
165  // cosmics
166  std::string cosmics = "nocosmics";
167  std::vector<Spectrum> s_cosmics;
168  if(usecosmics){
170  cosmics = "cosmics";
171  std::cout << "\nloading cosmics" << std::endl;
172  std::cout << "" << horn_one << std::endl;
173  std::string fcosmics_name_one = dcosm_name + "cosmics_" + horn_one + "__numu2018.root";
174  TFile fcosmics_one(fcosmics_name_one.c_str());
175  auto *cosmics_quant1_one = (TH1*)fcosmics_one.Get("cosmics_q1");
176  auto *cosmics_quant2_one = (TH1*)fcosmics_one.Get("cosmics_q2");
177  auto *cosmics_quant3_one = (TH1*)fcosmics_one.Get("cosmics_q3");
178  auto *cosmics_quant4_one = (TH1*)fcosmics_one.Get("cosmics_q4");
179  s_cosmics.push_back(Spectrum (cosmics_quant1_one, pot, livetime));
180  s_cosmics.push_back(Spectrum (cosmics_quant2_one, pot, livetime));
181  s_cosmics.push_back(Spectrum (cosmics_quant3_one, pot, livetime));
182  s_cosmics.push_back(Spectrum (cosmics_quant4_one, pot, livetime));
183 
184  if(twobeams){
186  std::cout << "" << horn_two << std::endl;
187  std::string fcosmics_name_two = dcosm_name + "cosmics_" + horn_two + "__numu2018.root";
188  TFile fcosmics_two(fcosmics_name_two.c_str());
189  auto *cosmics_quant1_two = (TH1*)fcosmics_two.Get("cosmics_q1");
190  auto *cosmics_quant2_two = (TH1*)fcosmics_two.Get("cosmics_q2");
191  auto *cosmics_quant3_two = (TH1*)fcosmics_two.Get("cosmics_q3");
192  auto *cosmics_quant4_two = (TH1*)fcosmics_two.Get("cosmics_q4");
193  s_cosmics.push_back(Spectrum (cosmics_quant1_two, pot, livetime));
194  s_cosmics.push_back(Spectrum (cosmics_quant2_two, pot, livetime));
195  s_cosmics.push_back(Spectrum (cosmics_quant3_two, pot, livetime));
196  s_cosmics.push_back(Spectrum (cosmics_quant4_two, pot, livetime));
197  }
198  }
199 
200 
201  std::vector <const IExperiment*> experiments;
202  if(quantN==0){// all quantiles
203  for(int quant = 0; quant < totquant; quant++){
204  if(usecosmics){
205  experiments.push_back(new SingleSampleExperiment(predictions[quant], *s_realdata[quant], s_cosmics[quant], 0.1));
206  }
207  else experiments.push_back(new SingleSampleExperiment(predictions[quant], *s_realdata[quant]));
208  }
209  }
210  else{//specific quantile
211  if(usecosmics){
212  experiments.push_back(new SingleSampleExperiment(predictions[quantN-1], *s_realdata[quantN-1], s_cosmics[quantN-1], 0.1));
213  }
214  else experiments.push_back(new SingleSampleExperiment(predictions[quantN-1], *s_realdata[quantN-1]));
215  }
216 
217  predictions.clear();
218  experiments.push_back(&kSolarConstraintsPDG2017);
219  experiments.push_back(WorldReactorConstraint2017());
220  MultiExperiment multiexpt(experiments);
221 
222  TFile* fout;
223  if(quantN==0) fout = new TFile(dout_name + "surfprof_" + dataName + "_" + cosmics + "_" + anaName + "__" + extrapName + "_" + systName + "_" + calcName + "_" + hierarchy + "__numu2018.root","RECREATE");
224  else fout = new TFile(dout_name + "surfprof_" + dataName + "_" + cosmics + "_" + anaName + "__" + extrapName + "_" + systName + "_" + calcName + "_" + hierarchy + Form("__numu2018__q%d.root", quantN),"RECREATE");
225 
226  RestartCalculator(calc);
227  FrequentistSurface* surface = new FrequentistSurface(&multiexpt, calc, &kFitSinSqTheta23, 40, sinmin, sinmax, &kFitDmSq32Scaled, 50, deltamin, deltamax, fit_sin23delta32, systs);
228  RestartCalculator(calc);
229  TH1* hprofile_delta32 = SqrtProfile(&multiexpt, calc, &kFitDmSq32Scaled, 50, deltamin, deltamax, -1, fit_sin23, systs);
230  RestartCalculator(calc);
231  TH1* hprofile_sin23 = SqrtProfile(&multiexpt, calc, &kFitSinSqTheta23, 40, sinmin, sinmax, -1, fit_delta32, systs);
232 
233  surface->SaveTo(fout, "surface");
234  hprofile_sin23->Write("profile_sin23");
235  hprofile_delta32->Write("profile_delta32");
236 
237  experiments.clear();
238 
239 
240 
241 } // end make_surfprof
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
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
double pot_two
Definition: saveFDMCHists.C:20
int totquant
Definition: saveFDMCHists.C:40
T sqrt(T number)
Definition: d0nt_math.hpp:156
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
bool twobeams
Definition: saveFDMCHists.C:39
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::string extrapName
Log-likelihood scan across two parameters.
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)
void make_surfprof(const std::string anaName="rhcfhc", const std::string extrapName="extrap", const int quantN=0, bool usecosmics=true, bool systematics=false, bool calc3a=true, bool nh=true)
Definition: make_surfprof.C:20
void RestartCalculator(osc::IOscCalcAdjustable *calc)
Definition: make_surfprof.C:12
#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)
std::string dataName
Definition: make_surfprof.C:9
const double pot_fhc
Definition: saveFDMCHists.C:10
double livetime
Definition: saveFDMCHists.C:21
void SaveTo(TDirectory *dir, const std::string &name) const
double livetime_one
Definition: saveFDMCHists.C:22
double deltamax
Definition: plot_shifts.C:10
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
General interface to any calculator that lets you set the parameters.
TFile * file
Definition: cellShifts.C:17
bool systematics
Definition: fnexvscaf.C:31
std::string period_two
Definition: saveFDMCHists.C:36
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.
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