46 double deltaCP = 0.18,
54 bool corrSysts =
false,
55 TString
options=
"joint_mockData_both_onlyIH",
56 bool dmsqSurf =
false,
57 bool th13Surf =
false)
60 cout<<
"\nLOAD PREDICTIONS\n\n"<<
endl;
61 bool nueOnly =
options.Contains(
"nueOnly");
62 bool numuOnly =
options.Contains(
"numuOnly");
63 bool joint =
options.Contains(
"joint");
64 assert (nueOnly || numuOnly || joint);
66 bool FHCOnly =
options.Contains(
"FHCOnly");
67 bool RHCOnly =
options.Contains(
"RHCOnly");
68 bool both =
options.Contains(
"both");
69 assert (FHCOnly || RHCOnly || both);
71 bool onlyNH =
options.Contains(
"onlyNH");
72 bool onlyIH =
options.Contains(
"onlyIH");
75 if(dmsqSurf) suffix+=
"_dmsq";
76 if(th13Surf) suffix+=
"_th13";
77 if(corrSysts) suffix+=
"_systs";
78 else suffix+=
"_stats";
81 string syst = (corrSysts?
"syst":
"stat");
83 TString
outfilename (syst +
"_hist_contours_" + suffix +
"_" + uni_tag);
84 TFile *
outfile =
new TFile(outfilename+
".root",
"recreate");
93 std::vector <ana::predictions> preds =
LoadPredictions (corrSysts,
true, decomp, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
94 std::vector <ana::predictions> preds_numu_only =
LoadPredictions (corrSysts,
true, decomp,
false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
95 std::vector <Spectrum >
data;
96 std::vector <Spectrum > data_numu_only;
97 std::vector <const IExperiment * > expts;
98 std::vector <const IExperiment * > expts_numu_only;
109 osc.Write(
"osc_params");
111 for (
int j = 0;
j< universes;
j++){
113 cout<<
"clear expts and data vectors"<<
endl;
115 expts_numu_only.clear();
117 data_numu_only.clear();
121 if(RHCOnly || both ) data.push_back(
GenerateFutureData(calc_fake, preds[1].pred, preds[1].
cos.first, RHCexp, RHClive));
124 if(FHCOnly || both ){
125 std::vector <Spectrum> numu_data;
127 numu_data.push_back(
GenerateFutureData(calc_fake, preds[3].pred, preds[3].
cos.first, FHCexp, FHClive));
128 numu_data.push_back(
GenerateFutureData(calc_fake, preds[4].pred, preds[4].
cos.first, FHCexp, FHClive));
129 numu_data.push_back(
GenerateFutureData(calc_fake, preds[5].pred, preds[5].
cos.first, FHCexp, FHClive));
130 data.insert(data.end(),numu_data.begin(), numu_data.end());
131 data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());
133 if(RHCOnly || both ){
134 std::vector <Spectrum> numu_data;
136 numu_data.push_back(
GenerateFutureData(calc_fake, preds[7].pred, preds[7].
cos.first, RHCexp, RHClive));
137 numu_data.push_back(
GenerateFutureData(calc_fake, preds[8].pred, preds[8].
cos.first, RHCexp, RHClive));
138 numu_data.push_back(
GenerateFutureData(calc_fake, preds[9].pred, preds[9].
cos.first, RHCexp, RHClive));
139 data.insert(data.end(),numu_data.begin(), numu_data.end());
140 data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());
144 cout<<
"check all vector sizes: "<<preds.size()<<
" data "<<data.size()<<
endl;
149 for(
int i = 0;
i <
int(preds.size()); ++
i){
150 auto thisdata = data[
i];
151 double POT = thisdata.POT();
152 double LT = thisdata.Livetime();
154 auto hcos = preds[
i].cos.first->ToTH1(LT, kBlack, kSolid,
kLivetime);
155 cout<<
i<<
" "<<preds[
i].name<<
" POT "<<POT<<
" tot MC "<<preds[
i].pred->Predict(calc_fake).Integral(POT)<<
156 " cos "<<hcos->Integral()<<
" cos er "<<preds[
i].cos.second<<
" analyze data "<<thisdata.Integral(POT)<<
endl;
160 cout<<
"Make numu only experiment to get the seeds"<<
endl;
162 for(
int i = 0;
i < (
int) preds_numu_only.size(); ++
i){
163 auto thisdata = data_numu_only[
i];
164 double POT = thisdata.POT();
165 double LT = thisdata.Livetime();
166 auto hcos = preds_numu_only[
i].cos.first->ToTH1(LT, kBlack, kSolid,
kLivetime);
167 cout<<
i<<
" "<<preds_numu_only[
i].name<<
" POT "<<POT<<
" tot MC "<<preds_numu_only[
i].pred->Predict(calc_fake).Integral(POT)<<
168 " cos "<<hcos->Integral()<<
" cos er "<<preds_numu_only[
i].cos.second<<
" analyze data "<<thisdata.Integral(POT)<<
endl;
169 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,
true));
177 std::cout <<
"Adding WorldReactorConstraint2017\n";
181 std::cout <<
"Adding Dmsq32ConstraintPDG2017\n";
184 std::cout <<
"Creating Multiexperiment with a total of " 185 << expts.size() <<
" experiments\n\n" <<
std::endl;
188 std::cout <<
"Creating Multiexperiment of numu only SimpleExp with a total of " 189 << expts_numu_only.size() <<
" experiments\n\n" <<
std::endl;
196 std::cout <<
"Systematics for the fit:\n";
198 bool ptExtrap =
true;
202 std::cout <<
"\n\nSystematic correlations...\n";
203 if(!nueOnly && ! numuOnly && corrSysts){
205 exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
true, ptExtrap));
207 for(
int i =0;
i < nnumu; ++
i) exptThis->SetSystCorrelations(
i+1, notfornumu);
210 exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
false, ptExtrap));
212 for(
int i =0;
i < nnumu; ++
i) exptThis->SetSystCorrelations(
i+1, notfornumu);
215 exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
true, ptExtrap));
216 exptThis->SetSystCorrelations(1,
GetCorrelations(
true,
false, ptExtrap));
218 for(
int i =0;
i < 4; ++
i) exptThis->SetSystCorrelations(
i+2, notfornumufhc);
220 for(
int i =0;
i < 4; ++
i) exptThis->SetSystCorrelations(
i+6, notfornumurhc);
223 if (nueOnly && corrSysts){
224 if (FHCOnly) exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
true, ptExtrap));
225 if (RHCOnly) exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
false, ptExtrap));
227 exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
true, ptExtrap));
228 exptThis->SetSystCorrelations(1,
GetCorrelations(
true,
false, ptExtrap));
232 std::cout <<
"Systematics for the numu only fit:\n";
235 if (corrSysts && numuOnly && both){
237 for(
int i =0;
i < 4; ++
i) {
238 exptThis->SetSystCorrelations(
i, notfornumufhc);
239 exptThis_numu_only->SetSystCorrelations(
i, notfornumufhc);
242 for(
int i =0;
i < 4; ++
i) {
243 exptThis->SetSystCorrelations(
i+4, notfornumurhc);
244 exptThis_numu_only->SetSystCorrelations(
i+4, notfornumurhc);
257 std::vector <SystShifts> seedShifts = {};
264 cout<<
"------------------- Start prestage seeds --------------------------"<<
endl;
266 double minchi_numu = 1E20;
267 double pre_seed_th23;
268 double pre_seed_dmsq;
272 double maxmix = 0.514;
273 double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
276 std::cout <<
"\n\nFinding test best fit " << (hie>0?
"NH " :
"IH ") <<
"\n\n";
279 auto thisminchi = fitnumu_only.
Fit(calc, auxShifts ,
284 if(thisminchi < minchi_numu) minchi_numu = thisminchi;
289 numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
290 numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));
295 cout<<
"------------------- End prestage seeds --------------------------"<<
endl;
298 std::vector<double> seeds;
304 std::vector <th23helper> th23seeds;
314 std::vector<double> th23_all_seeds_NH = {numu_pre_seedLONH, 0.5, numu_pre_seedUONH};
316 std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
317 std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
323 double minchi23 = 1E20;
double thisdmsq = -1;
double thisth23 = -1;
double thisdcp = -1;
double thisth13 = -1;
324 double curdmsq, curth23, curdcp, curth13;
327 for (
auto & thseed: th23seeds){
330 << (hie>0?
" NH " :
" IH ")
332 for (
auto const &
s:thseed.seeds)
std::cout <<
s <<
", ";
334 cout<<
"pre seed from numu is "<<pre_seed_dmsq<<
endl;
341 auto thisminchi = fit23.
Fit(calc, auxShifts,
SeedList({
343 {thseed.var, thseed.seeds},
349 curth13 = kFitSinSq2Theta13.GetValue(calc);
350 curth23 = thseed.var->GetValue(calc);
358 string hietag = (hie>0?
"NH":
"IH");
360 if(thisminchi<minchi23){
362 thisth13 = kFitSinSq2Theta13.GetValue(calc);
363 thisth23 = thseed.var->GetValue(calc);
365 cout <<
"this values (dmsq, th13, th23, dcp): "<<thisdmsq<<
" "<<thisth13<<
" "<<thisth23<<
" "<<thisdcp<<
endl;
367 minchi23=
min(minchi23, thisminchi);
380 cout<<
"minchi wrote into the file "<<minchi23<<
endl;
385 std::cerr <<
"\n WARNING Using 20 bins for dmsq surface\n\n";
389 for(
int hie: {-1, +1}){
391 if((onlyNH && hie<0) || (onlyIH && hie>0))
continue;
393 std::cout <<
"Starting surface " << (hie>0?
"NH ":
"IH") <<
"\n\n";
394 cout<<
"the delta seeds are ";
395 for (
auto const &
s:delta_seeds)
std::cout <<
s <<
", ";
398 if(!th13Surf && ! dmsqSurf){
401 double low = (nueOnly)?0:0.3;
402 double high = (nueOnly)?1:0.7;
403 cout<<
"low "<<low<<
" high "<<high<<
endl;
411 auto surf1=surf23.
ToTH2(minchi23);
413 surf23.SaveTo(d, (TString)
"surf_delta_th23_univ"+
std::to_string(
j)+
"_"+hieStr);
415 if(!th13Surf && dmsqSurf){
418 double lowth23 = (RHCOnly)?0.2:0.35;
419 double highth23 = (RHCOnly)?1:0.7;
420 cout<<
"lowth23 "<<lowth23<<
" highth23 "<<highth23<<
endl;
421 double lowdmsq = (hie>0?(RHCOnly?2.0:2.1):(RHCOnly?-3.0:-2.7));
422 double highdmsq = (hie>0?(RHCOnly?3.0:2.7):(RHCOnly?-2.0:-2.2));
428 auto surf6=surf23m.
ToTH2(minchi23);
429 surf6->Write((TString)
"th23_dm32_univ"+
std::to_string(
j)+
"_"+(hie>0?
"NH":
"IH"));
442 auto surf4 = surf13.
ToTH2(minchi23);
444 surf4->Write((TString)
"th13_delta_univ"+
std::to_string(
j)+
"_"+(hie>0?
"NH":
"IH"));
Cuts and Vars for the 2020 FD DiF Study.
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
fvar< T > fabs(const fvar< T > &x)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var's GetValue()
Simple record of shifts applied to systematic parameters.
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
const Dmsq32Constraint kDmsq32ConstraintPDG2019(2.444e-3, 0.034e-3, 2.55e-3, 0.04e-3)
std::vector< predictions > LoadPredictions(bool corrSysts=false, bool runOnGrid=false, std::string decomp="", bool nue=true, bool numu=true, bool fhc=true, bool rhc=true, bool NERSC=false)
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
Log-likelihood scan across two parameters.
string outfilename
knobs that need extra care
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.
const double kAna2020FHCLivetime
virtual T GetDmsq32() const
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
const double kAna2020FHCPOT
void joint_fit_future_contour_univ(int universes=2, TString uni_tag="", double deltaCP=0.18, double dmsq32=2.51e-3, double th23=0.58, string sign="plus", double FHCexp=kAna2020FHCPOT, double RHCexp=kAna2020RHCPOT, double FHClive=kAna2020FHCLivetime, double RHClive=kAna2020RHCLivetime, bool corrSysts=false, TString options="joint_mockData_both_onlyIH", bool dmsqSurf=false, bool th13Surf=false)
const double kAna2020RHCPOT
Oscillation probability calculators.
std::vector< double > POT
static float min(const float a, const float b, const float c)
TH2 * ToTH2(double minchi=-1) const
Combine multiple component experiments.
void SaveTo(TDirectory *dir, const std::string &name) const
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
const double kAna2020RHCLivetime
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
const ReactorExperiment * WorldReactorConstraint2019()
Reactor constraint from PDG 2019.
const FitSinSq2Theta13 kFitSinSq2Theta13
std::string to_string(ModuleType mt)
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)
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.