47 TString
options=
"joint_realData_both")
50 bool nueOnly =
options.Contains(
"nueOnly");
51 bool numuOnly =
options.Contains(
"numuOnly");
52 bool joint =
options.Contains(
"joint");
53 assert (nueOnly || numuOnly || joint);
55 bool FHCOnly =
options.Contains(
"FHCOnly");
56 bool RHCOnly =
options.Contains(
"RHCOnly");
57 bool both =
options.Contains(
"both");
58 assert (FHCOnly || RHCOnly || both);
60 bool fake2018 =
options.Contains(
"fake2018");
61 bool realData =
options.Contains(
"realData");
64 if(corrSysts) suffix+=
"_systs";
65 else suffix+=
"_stats";
69 TString
outfilename (outdir +
"bestfits_" + suffix);
79 std::pair <TH1D*, double>
cos;
83 std::vector <predictions> preds;
84 std::vector <predictions> preds_numu_only;
85 std::vector <Spectrum * >
data;
86 std::vector <Spectrum * > data_numu_only;
87 std::vector <const IExperiment * > expts;
88 std::vector <const IExperiment * > expts_numu_only;
91 if(fake2018)
SetFakeCalc(calc_fake, 0.58, 2.51
e-3, 0.17);
92 else if(!realData) {
std::cerr <<
"need setting for data\n";
exit(1);}
95 if(FHCOnly || both ) {
98 if(RHCOnly || both ) {
105 if(FHCOnly || both ) {
108 for (
int i = 0;
i< nnumu;
i++ ){
113 if(RHCOnly || both ) {
116 for (
int i = 0;
i< nnumu;
i++ ){
130 if(FHCOnly || both ){
132 data.insert(data.end(),numu_data.begin(), numu_data.end());
133 data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());}
134 if(RHCOnly || both ){
136 data.insert(data.end(),numu_data.begin(), numu_data.end());
137 data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());}
141 for(
int i = 0;
i <
int(preds.size()); ++
i){
142 double POT = preds[
i].pot;
144 if(realData) thisdata = data[
i];
145 cout<<i<<
" "<<preds[
i].name<<
" POT "<<POT<<
" tot MC "<<preds[
i].pred->Predict(calc_fake).Integral(POT)<<
" cos "<<preds[
i].cos.first->Integral()<<
" cos er "<<preds[
i].cos.second<<
" analyze data "<<thisdata->Integral(POT)<<
endl;
149 preds[
i].cos.first->SetMarkerStyle(34);
151 preds[
i].cos.first->Draw(
"hist p same");
152 gPad->Print(outdir +
"debug_predictions_" + suffix +
"_" +
std::to_string(i) +
".pdf");
155 cout<<
"Make numu only experiment to get the seeds"<<
endl;
158 for(
int i = 0;
i < (
int) preds_numu_only.size(); ++
i){
159 double POT = preds_numu_only[
i].pot;
160 auto thisdata =
GetFakeData (preds_numu_only[
i].
pred, calc_fake, POT, preds_numu_only[
i].
cos.first);
161 if(realData) thisdata = data_numu_only[
i];
162 cout<<i<<
" "<<preds_numu_only[
i].name<<
" POT "<<POT<<
" tot MC "<<preds_numu_only[
i].pred->Predict(calc_fake).Integral(POT)<<
" cos "<<preds_numu_only[
i].cos.first->Integral()<<
" cos er "<<preds_numu_only[
i].cos.second<<
" analyze data "<<thisdata->Integral(POT)<<
endl;
163 expts_numu_only.push_back(
new SingleSampleExperiment(preds_numu_only[i].pred, *thisdata, preds_numu_only[i].
cos.first, preds_numu_only[i].cos.second));
173 std::cout <<
"Adding Dmsq32ConstraintPDG2017\n";
176 std::cout <<
"Creating Multiexperiment with a total of " 177 << expts.size() <<
" experiments\n\n" <<
std::endl;
180 std::cout <<
"Creating Multiexperiment of numu only SimpleExp with a total of " 181 << expts_numu_only.size() <<
" experiments\n\n" <<
std::endl;
187 std::cout <<
"Systematics for the fit:\n";
190 std::cout <<
"\n\nSystematic correlations...\n";
191 if(!nueOnly && ! numuOnly && corrSysts){
195 for(
int i =0;
i < nnumu; ++
i) exptThis->SetSystCorrelations(
i+1, notfornumu);
200 for(
int i =0;
i < nnumu; ++
i) exptThis->SetSystCorrelations(
i+1, notfornumu);
206 for(
int i =0;
i < 8; ++
i) exptThis->SetSystCorrelations(
i+2, notfornumu);
209 if (nueOnly && corrSysts){
210 if (FHCOnly) exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
true));
211 if (RHCOnly) exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
false));
218 std::cout <<
"Systematics for the numu only fit:\n";
231 std::vector <SystShifts> seedShifts = {};
237 cout<<
"------------------- Start prestage seeds --------------------------"<<
endl;
239 double minchi_numu = 1E20;
240 double pre_seed_th23;
241 double pre_seed_dmsq;
245 double maxmix = 0.514;
246 double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
249 std::cout <<
"\n\nFinding test best fit " << (hie>0?
"NH " :
"IH ") <<
"\n\n";
252 auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
257 if(thisminchi < minchi_numu) minchi_numu = thisminchi;
260 pre_seed_dmsq = kFitDmSq32.GetValue(calc);
263 numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
264 numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));
269 cout<<
"------------------- End prestage seeds --------------------------"<<
endl;
273 std::vector<double> seeds;
279 std::vector <th23helper> th23seeds;
290 std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
291 std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
298 TFile *
outfile =
new TFile(outfilename+
".root",
"recreate");
301 double minchi23 = 1E20;
302 for(
int hie:{-1, 1}){
303 for (
auto & thseed:th23seeds){
305 std::cout <<
"\n\nFinding best fit " << (hie > 0 ?
"NH " :
"IH ")
306 << thseed.label <<
", " 316 Fitter fit23(exptThis, fitvars,
systs);
317 cout<<
" pre dmsq seed is "<<pre_seed_dmsq<<
endl;
318 auto thisminchi = fit23.Fit(calc, auxShifts,
321 { thseed.var, thseed.seeds },
326 minchi23=
min(minchi23, thisminchi);
331 TString
str =
"Best fit " ;
332 for (
auto &
v:fitvars){
336 shifts->SetTitle(str);
338 TLine *
l=
new TLine(gPad->GetUxmin(),0,gPad->GetUxmax(),0);
340 gPad->Print(outdir +
"debug_slice_shifts_bestfit_" + suffix +
341 (hie>0?
"NH":
"IH") + thseed.label +
".pdf");
343 TString dirname = (hie > 0 ?
"NH_" :
"IH_") + thseed.label;
345 auxShifts.
SaveTo(outfile, dirname+
"_systs");
350 std::cout <<
"\nFound overall minchi " << minchi23 <<
"\n\n";
352 cout<<
"Check with oldstyle seeds "<<
endl;
353 double minchi23test = 1E20;
355 std::cout <<
"\n\nFinding test best fit " 356 << (hie>0?
"NH " :
"IH ")
363 Fitter fit23(exptThis, fitvars,
systs);
364 auto thisminchi = fit23.Fit(calc,auxShifts ,{
370 minchi23test=
min(minchi23test, thisminchi);
373 TString
str =
"Best fit " ;
374 for (
auto &
v:fitvars){
378 shifts->SetTitle(str);
379 gPad->Print(outdir +
"debug_slice_shifts_bestfit_" + suffix +
380 (hie>0?
"NH":
"IH") +
"test_chisq_with_no_octant_differentiation.pdf");
389 cout<<
"<<<<<<<<<<<<<<<<<<<<<<<<<< Test with only systematics free"<<
endl;
413 Fitter fit23(exptThis, {},
systs);
414 auto thisminchi = fit23.Fit(calc,auxShifts ,{}, seedShifts,
416 ofstream txtfile(outdir +
"debug_slice_shifts_bestfit_" + suffix +
417 +
"test_chisq_with_osc_par_fixed.txt");
419 txtfile<<
v->ShortName().c_str()<<
"="<<auxShifts.
GetShift(
v)<<
endl;
422 minchi23test=
min(minchi23test, thisminchi);
424 TString
str =
"Best fit " ;
426 shifts->SetTitle(str);
427 gPad->Print(outdir +
"debug_slice_shifts_bestfit_" + suffix +
428 +
"test_chisq_with_osc_par_fixed.pdf");
431 cout<<
"<<<<<<<<<<<< This is not true BF : "<< minchi23test<<
endl;
Cuts and Vars for the 2020 FD DiF Study.
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
TH1 * PlotSystShifts(const SystShifts &shifts, bool sortName)
fvar< T > fabs(const fvar< T > &x)
void SaveTo(const osc::IOscCalc &x, TDirectory *dir, const std::string &name)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
const Color_t kCosmicBackgroundColor
const FitDmSq32 kFitDmSq32
Simple record of shifts applied to systematic parameters.
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
General interface to oscillation calculators.
virtual void SetTh13(const T &th13)=0
static SystShifts Nominal()
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
virtual void SetDmsq32(const T &dmsq32)=0
const XML_Char const XML_Char * data
string outfilename
knobs that need extra care
void joint_fit_2019_bestfit(bool corrSysts=true, TString options="joint_realData_both")
T GetShift(const ISyst *syst) const
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
std::vector< double > POT
static float min(const float a, const float b, const float c)
Combine multiple component experiments.
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
double GetPOT(bool isFHC=true)
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
Standard interface to all prediction techniques.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
const FitSinSq2Theta13 kFitSinSq2Theta13
void SaveTo(TDirectory *dir, const std::string &name) const
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
std::string to_string(ModuleType mt)
std::pair< Spectrum *, double > cos
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool isFHC=true, bool isRHC=true, bool intersection=true, bool ptExtrap=true)
virtual void SetdCP(const T &dCP)=0
Compare a single data spectrum to the MC + cosmics expectation.