44 bool runOnGrid =
false,
45 TString
options=
"joint_fake2019_both_onlyIH",
46 bool dmsqSurf =
false,
47 bool th13Surf =
false)
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 fake2019 =
options.Contains(
"fake2019");
61 bool realData =
options.Contains(
"realData");
63 bool onlyNH =
options.Contains(
"onlyNH");
64 bool onlyIH =
options.Contains(
"onlyIH");
67 if(FHCOnly) tag=
"FHCOnly/";
68 if(RHCOnly) tag=
"RHCOnly/";
69 if(both) tag =
"RHC_and_FHC/";
72 if(dmsqSurf) suffix+=
"_dmsq";
73 if(th13Surf) suffix+=
"_th13";
74 if(corrSysts) suffix+=
"_systs";
75 else suffix+=
"_stats";
77 TString
outdir =
"/nova/ana/nu_e_ana/Ana2020/prefits/";
78 if (runOnGrid) outdir =
"./";
79 TString outFCfilename (outdir +
"contours_FCInput_2020_" + suffix);
87 TString
outfilename (outdir +
"hist_contours_2020_" + suffix);
96 std::vector <ana::predictions> preds =
LoadPredictions (corrSysts, runOnGrid, decomp, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
97 std::vector <ana::predictions> preds_numu_only =
LoadPredictions (corrSysts, runOnGrid, decomp,
false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
98 std::vector <Spectrum * >
data =
LoadRealData (runOnGrid, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
99 std::vector <Spectrum * > data_numu_only =
LoadRealData (runOnGrid,
false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
100 std::vector <const IExperiment * > expts;
101 std::vector <const IExperiment * > expts_numu_only;
107 else if(!realData) {
std::cerr <<
"need setting for data\n";
exit(1);}
109 for(
int i = 0;
i <
int(preds.size()); ++
i){
110 double POT = preds[
i].pot;
111 auto thisdata =
GetFakeData (preds[
i].
pred, calc_fake, POT, preds[
i].
cos.first, preds[i].livetime);
112 if(realData) thisdata = data[
i];
114 cout<<i<<
" "<<preds[
i].name<<
" POT "<<POT<<
" tot MC "<<preds[
i].pred->Predict(calc_fake).Integral(POT)<<
115 " cos "<<hcos->Integral()<<
" cos er "<<preds[
i].cos.second<<
" analyze data "<<thisdata->Integral(thisdata->POT())<<
endl;
120 hcos->SetMarkerStyle(34);
122 hcos->Draw(
"hist p same");
123 gPad->Print(outdir + (corrSysts?
"syst/":
"stat/")+
"debug_predictions_" + suffix +
"_" +
std::to_string(i) +
".pdf");
127 cout<<
"Make numu only experiment to get the seeds"<<
endl;
130 for(
int i = 0;
i < (
int) preds_numu_only.size(); ++
i){
131 double POT = preds_numu_only[
i].pot;
132 auto thisdata =
GetFakeData (preds_numu_only[
i].
pred, calc_fake, POT, preds_numu_only[
i].
cos.first, preds_numu_only[i].livetime);
133 if(realData) thisdata = data_numu_only[
i];
134 auto hcos = preds_numu_only[
i].cos.first->ToTH1(preds_numu_only[i].
livetime, kBlack, kSolid,
kLivetime);
135 cout<<i<<
" "<<preds_numu_only[
i].name<<
" POT "<<POT<<
" tot MC "<<preds_numu_only[
i].pred->Predict(calc_fake).Integral(POT)<<
136 " cos "<<hcos->Integral()<<
" cos er "<<preds_numu_only[
i].cos.second<<
" analyze data "<<thisdata->Integral(thisdata->POT())<<endl;
137 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));
145 std::cout <<
"Adding WorldReactorConstraint2019\n";
149 std::cout <<
"Adding Dmsq32ConstraintPDG2019\n";
152 std::cout <<
"Creating Multiexperiment with a total of " << expts.size() <<
" experiments\n\n" <<
std::endl;
155 std::cout <<
"Creating Multiexperiment of numu only SimpleExp with a total of " 156 << expts_numu_only.size() <<
" experiments\n\n" <<
std::endl;
161 preds_numu_only.clear();
167 bool ptExtrap =
true;
168 if (decomp ==
"noPt") ptExtrap =
false;
171 std::cout <<
"\n\nSystematic correlations...\n";
172 if(!nueOnly && ! numuOnly && corrSysts){
174 exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
true, ptExtrap));
176 for(
int i =0;
i < nnumu; ++
i) exptThis->SetSystCorrelations(
i+1, notfornumu);
179 exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
false, ptExtrap));
181 for(
int i =0;
i < nnumu; ++
i) exptThis->SetSystCorrelations(
i+1, notfornumu);
184 exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
true, ptExtrap));
185 exptThis->SetSystCorrelations(1,
GetCorrelations(
true,
false, ptExtrap));
187 for(
int i =0;
i < 4; ++
i) exptThis->SetSystCorrelations(
i+2, notfornumufhc);
189 for(
int i =0;
i < 4; ++
i) exptThis->SetSystCorrelations(
i+6, notfornumurhc);
192 if (nueOnly && corrSysts){
193 if (FHCOnly) exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
true, ptExtrap));
194 if (RHCOnly) exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
false, ptExtrap));
196 exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
true, ptExtrap));
197 exptThis->SetSystCorrelations(1,
GetCorrelations(
true,
false, ptExtrap));
201 std::cout <<
"Systematics for the numu only fit:\n";
204 if (corrSysts && numuOnly && both){
206 for(
int i =0;
i < 4; ++
i) {
207 exptThis->SetSystCorrelations(
i, notfornumufhc);
208 exptThis_numu_only->SetSystCorrelations(
i, notfornumufhc);
211 for(
int i =0;
i < 4; ++
i) {
212 exptThis->SetSystCorrelations(
i+4, notfornumurhc);
213 exptThis_numu_only->SetSystCorrelations(
i+4, notfornumurhc);
226 std::vector <SystShifts> seedShifts = {};
228 for (
double systshift:{-2, +2}){
230 seedShifts.push_back(tempShifts);
238 cout<<
"------------------- Start prestage seeds --------------------------"<<
endl;
240 double minchi_numu = 1E20;
241 double pre_seed_th23;
242 double pre_seed_dmsq;
246 double maxmix = 0.514;
247 double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
250 std::cout <<
"\n\nFinding test best fit " << (hie>0?
"NH " :
"IH ") <<
"\n\n";
253 auto thisminchi = fitnumu_only.
Fit(calc, auxShifts ,
258 if(thisminchi < minchi_numu) minchi_numu = thisminchi;
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;
272 std::vector<double> seeds;
277 std::vector <th23helper> th23seeds;
287 std::vector<double> th23_all_seeds_NH = {numu_pre_seedLONH, 0.5, numu_pre_seedUONH};
289 std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
290 std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
296 double minchi23 = 1E20;
298 for (
auto & thseed: th23seeds){
299 std::cout <<
"\n\nFinding best fit " << thseed.label
300 << (hie>0?
" NH " :
" IH ") <<
"\n"<<endl;
301 for (
auto const &
s:thseed.seeds)
std::cout <<
s <<
", ";
303 cout<<
"pre seed from numu is "<<pre_seed_dmsq<<
endl;
312 auto thisminchi = fit23.
Fit(calc, auxShifts,
SeedList({
314 {thseed.var, thseed.seeds},
318 minchi23=
min(minchi23, thisminchi);
325 cout<<
"Check with oldstyle seeds when we didn't split into octant"<<
endl;
326 double minchi23test = 1E20;
328 std::cout <<
"\n\nFinding test best fit " << (hie>0?
"NH " :
"IH ") <<
"\n\n";
338 auto thisminchi = fit23.
Fit(calc,auxShifts ,
SeedList({
342 }),seedShifts)->EvalMetricVal();
343 minchi23test=
min(minchi23test, thisminchi);
349 TFile *
outfile =
new TFile(outfilename+
".root",
"recreate");
350 TFile * outFCfile =
new TFile(outFCfilename+
".root",
"recreate");
355 cout<<
"minchi wrote into the file "<<minchi23<<
endl;
356 v.Write(
"overall_minchi");
360 std::cerr <<
"\n WARNING Using 20 bins for dmsq surface\n\n";
364 for(
int hie: {-1, +1}){
366 if((onlyNH && hie<0) || (onlyIH && hie>0))
continue;
368 std::cout <<
"Starting surface " << (hie>0?
"NH ":
"IH") <<
"\n\n";
369 cout<<
"the delta seeds are ";
370 for (
auto const &
s:delta_seeds)
std::cout <<
s <<
", ";
373 if(!th13Surf && ! dmsqSurf){
376 double low = (nueOnly)?0:0.3;
377 double high = (nueOnly)?1:0.7;
378 cout<<
"low "<<low<<
" high "<<high<<
endl;
386 auto surf1=surf23.
ToTH2(minchi23);
388 surf1->Write((TString)
"delta_th23_"+hieStr);
389 surf23.SaveTo(outfile,
"surf_delta_th23_"+hieStr);
392 std::vector<TH2*> profhists = surf23.GetProfiledHists();
394 profhists[0]->Write((hieStr+
"_DmSq32").c_str());
395 profhists[1]->Write((hieStr+
"_SinSq2Theta13").c_str());
396 for (
int i = 0;
i <
int(systs.size());
i++){
397 profhists[
i+2]->Write((hieStr+
"_"+systs[
i]->ShortName()).c_str());
400 if(!th13Surf && dmsqSurf){
403 double lowth23 = (RHCOnly)?0.2:0.35;
404 double highth23 = (RHCOnly)?1:0.7;
405 cout<<
"lowth23 "<<lowth23<<
" highth23 "<<highth23<<
endl;
406 double lowdmsq = (hie>0?(RHCOnly?2.0:2.1):(RHCOnly?-3.0:-2.7));
407 double highdmsq = (hie>0?(RHCOnly?3.0:2.7):(RHCOnly?-2.0:-2.2));
408 cout<<
"lowdmsq "<<lowdmsq<<
" highdmsq "<<highdmsq<<
endl;
414 auto surf6=surf23m.
ToTH2(minchi23);
416 surf6->Write((TString)
"th23_dm32_"+(hie>0?
"NH":
"IH"));
422 profhists[0]->Write((hieStr+
"_SinSq2Theta13").c_str());
423 profhists[1]->Write((hieStr+
"_DeltaCPpi").c_str());
424 for (
int i = 0;
i <
int(systs.size());
i++){
425 profhists[
i+2]->Write((hieStr+
"_"+systs[
i]->ShortName()).c_str());
438 auto surf4 = surf13.
ToTH2(minchi23);
440 surf4->Write((TString)
"th13_delta_"+(hie>0?
"NH":
"IH"));
446 profhists[0]->Write((hieStr+
"_SinSqTheta23").c_str());
447 profhists[1]->Write((hieStr+
"_DmSq32").c_str());
448 for (
int i = 0;
i <
int(systs.size());
i++){
449 profhists[
i+2]->Write((hieStr+
"_"+systs[
i]->ShortName()).c_str());
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
const Color_t kCosmicBackgroundColor
void run_joint_fit_2020_contours(bool corrSysts=false, bool runOnGrid=false, TString options="joint_fake2019_both_onlyIH", bool dmsqSurf=false, bool th13Surf=false)
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
std::vector< TH2 * > GetProfiledHists()
Maps of the values taken on by the profiled parameters.
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.
virtual T GetDmsq32() const
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
std::vector< double > POT
static float min(const float a, const float b, const float c)
TH2 * ToTH2(double minchi=-1) const
std::vector< Spectrum * > LoadRealData(bool runOnGrid=false, bool nue=true, bool numu=true, bool fhc=true, bool rhc=true)
Combine multiple component experiments.
void SaveTo(TDirectory *dir, const std::string &name) const
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
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
TH1 * DataMCComparison(const Spectrum &data, const Spectrum &mc, EBinType bintype)
std::string to_string(ModuleType mt)
const DummyAnaSyst kAnaLightlevelNDSyst("Light_Level_ND","Light_Level_ND")
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.