53 int N,
int jid,
int npts_flush = 0,
57 TString
beam =
"both",
bool corrSysts =
false,
58 bool readFromFile =
false,
bool runOnGrid=
false)
61 plot ==
"deltassth23" &&
"Plot must be ssth23dmsq32 or deltassth23");
65 if(readFromFile && !corrSysts){
72 bool fake2019 =
options.find(
"fake2019") != std::string::npos;
73 bool realData =
options.find(
"realData") != std::string::npos;
75 bool RHCOnly =
beam.Contains(
"RHCOnly");
76 bool FHCOnly =
beam.Contains(
"FHCOnly");
77 bool both =
beam.Contains(
"both");
79 bool nueOnly =
options.find(
"nueOnly") != std::string::npos;
80 bool numuOnly =
options.find(
"numuOnly") != std::string::npos;
81 bool joint =
options.find(
"joint") != std::string::npos;
88 TotBins = (
plot ==
"deltassth23") ? 900 : 400;
90 std::memset(bins, 0, TotBins);
102 path +=
"/binlists/";
104 TString binlist_fname = path+
plot+hier+
"binlist.txt";
106 binlist.open(binlist_fname);
118 }
while(!binlist.eof());
123 int ibin = bins[startidx];
124 int fbin = bins[endidx];
129 if(both) fcout+=
"both";
130 else fcout += RHCOnly?
"RHCOnly":
"FHCOnly";
131 fcout += corrSysts?
"_systs_":
"_stats_";
143 if(fake2019)
SetFakeCalc(calc, 0.564956, 2.48304
e-3, 1.9997);
145 else if(!realData) {
std::cerr <<
"need setting for data\n";
exit(1);}
162 std::vector <ana::predictions> preds =
LoadPredictions (corrSysts, runOnGrid,
decomp, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
163 std::vector <ana::predictions> preds_numu_only =
LoadPredictions (corrSysts, runOnGrid,
decomp,
false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
174 if(FHCOnly || both ) {
175 preds.push_back({"nue_fhc",
176 GetNuePrediction2020(decomp+"combo", DefaultOscCalc(), corrSysts, "fhc",
178 GetNueCosmics2020("fhc", runOnGrid),
180 fc::kFHC, fc::kNue});
182 if(RHCOnly || both ) {
183 preds.push_back({"nue_rhc",
184 GetNuePrediction2020("prop", DefaultOscCalc(), corrSysts, "rhc",
186 GetNueCosmics2020("rhc", runOnGrid),
187 GetPOT(false), GetLT(false),
188 fc::kRHC, fc::kNue});
191 std::cout << "\n------------------Done Getting nue preds ---------------------\n";
197 // Construct numu inputs
200 for(fc::kBeam_t _beam : {fc::kFHC, fc::kRHC}) {
201 if(((_beam == fc::kRHC) && FHCOnly) ||
202 ((_beam == fc::kFHC) && RHCOnly))
204 double _pot = GetPOT(_beam == fc::kFHC);
205 double _livetime = GetLT (_beam == fc::kFHC);
206 std::string beam_name = _beam == fc::kFHC ? "fhc" : "rhc";
207 auto numu_preds = GetNumuPredictions2020(nnumu, decomp, DefaultOscCalc(), corrSysts, beam_name, false, false);
208 auto numu_cosmics = GetNumuCosmics2020 (nnumu, beam_name, false);
209 for (int i = 0; i< nnumu; i++ ){
210 preds.push_back({"numu_" + beam_name + "_" +std::to_string(i+1),
217 preds_numu_only.push_back({"numu_" + beam_name + "_"+std::to_string(i+1),
229 bool ptExtrap =
false;
236 std::vector<const ISyst*> snue_fhc = {};
238 std::vector<const ISyst*> snue_rhc = {};
242 std::vector<const ISyst*> snumu = {};
255 if(
plot ==
"ssth23dmsq32") p =
"dmsq_";
256 if(
plot ==
"deltath13") p =
"th13_";
259 std::string helperupsname =
"/nova/app/users/scalvez/fchelper/";
260 TString anasurfname = helperupsname
261 +
"/Surfaces/hist_contours_2019_joint_realData_"+
beam+
"_only";
262 anasurfname += nh?
"NH_" :
"IH_";
264 anasurfname += corrSysts?
"systs.root":
"stats.root";
266 TFile *anasurffile = TFile::Open(anasurfname);
269 if(
plot ==
"ssth23dmsq32") plotShortName =
"th23_dm32";
270 if(nh) plotShortName +=
"_NH";
271 else plotShortName +=
"_IH";
272 TH2* anasurf = (TH2F*)anasurffile->Get(plotShortName.c_str());
277 helperupsname+
"FCInputs/contours_FCInput_2019_joint_realData_"+beam+
"_only";
278 helpername += nh?
"NH_" :
"IH_";
280 helpername += corrSysts?
"systs.root":
"stats.root";
282 TFile *fchelp = TFile::Open(helpername);
289 int NX = anasurf->GetXaxis()->GetNbins();
290 int NY = anasurf->GetYaxis()->GetNbins();
291 double minX = anasurf->GetXaxis()->GetBinCenter(1);
292 double maxX = anasurf->GetXaxis()->GetBinCenter(NX);
293 double minY = anasurf->GetYaxis()->GetBinCenter(1);
294 double maxY = anasurf->GetYaxis()->GetBinCenter(NY);
296 double widthX = anasurf->GetXaxis()->GetBinWidth(1);
297 double widthY = anasurf->GetYaxis()->GetBinWidth(1);
300 assert(fbin < NX*NY &&
"Invalid bin number");
303 for (
int ibin = startidx; ibin <= endidx; ibin++) {
308 if(bins[ibin] == 0 && ibin != 0)
continue;
313 int binX = bin - NX*binY;
317 double centerX = anasurf->GetXaxis()->GetBinCenter(binX+1);
318 double centerY = anasurf->GetYaxis()->GetBinCenter(binY+1);
324 double seedssth23 = 0;
325 if(
plot ==
"deltath13"){
326 TH2* hssth23 = (TH2F*) fchelp->Get((hierStr+
"_SinSqTheta23").c_str());
327 seedssth23 = hssth23->GetBinContent(binX+1,binY+1);
330 double seeddmsq32 = 0;
331 if (
plot ==
"deltassth23"){
332 TH2* hdmsq = (TH2F*)fchelp->Get((hierStr+
"_DmSq32").c_str());
333 seeddmsq32 = hdmsq->GetBinContent(binX+1,binY+1);
336 double seeddelta = 0;
337 if (
plot ==
"ssth23dmsq32"){
338 TH2* hdelta = (TH2F*)fchelp->Get((hierStr+
"_DeltaCPpi").c_str());
339 seeddelta = hdelta->GetBinContent(binX+1,binY+1);
342 TH2* hss2th13 = (TH2F*)fchelp->Get((hierStr+
"_SinSq2Theta13").c_str());
343 double seedss2th13 = hss2th13->GetBinContent(binX+1,binY+1);
347 std::map<std::string,double> seedsyst;
349 TH2*
h = (TH2F*)fchelp->Get((hierStr+
"_"+syst->ShortName()).c_str());
351 std::cout <<
"Don't see a prof hist for " << syst->ShortName() <<
". Continuing, but check your ups version." <<
std::endl;
354 double seedval = h->GetBinContent(binX+1,binY+1);
355 seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
367 for (
int i = 0;
i < NPts;
i++){
375 double trueX = centerX;
376 double trueY = centerY;
382 if(
plot !=
"deltath13"){
385 dmsq_var->
SetValue(calc,
plot==
"deltassth23"?seeddmsq32:trueY);
391 dmsq_var->
SetValue(calc, seeddmsq32);
398 for (
const ISyst *syst : snue_fhc){
399 if (seedsyst.find(syst->ShortName()) == seedsyst.end())
continue;
400 seednue_fhc.
SetShift(syst,seedsyst[syst->ShortName()]);
403 for (
const ISyst *syst : snue_rhc) {
404 if (seedsyst.find(syst->ShortName()) == seedsyst.end())
continue;
405 seednue_rhc.
SetShift(syst, seedsyst[syst->ShortName()]);
409 for (
const ISyst *syst : snumu){
410 if (seedsyst.find(syst->ShortName()) == seedsyst.end())
continue;
411 seednumu.
SetShift(syst,seedsyst[syst->ShortName()]);
418 std::cout <<
"----------Making mock data at-----------------\n";
449 std::vector<const IExperiment*> expts;
450 std::vector<const IExperiment*> numu_expts;
458 for(
int ipred = 0; ipred <
int(preds.size()); ++ipred) {
459 double POT = preds[ipred].pot;
461 auto fakedata = preds[ipred].pred->Predict(calc);
462 auto hcos = preds[ipred].cos.first->ToTH1(preds[
i].
livetime, kBlack, kSolid,
kLivetime);
466 auto mockdata = fakedata.MockData(POT);
470 for(
int ipred = 0; ipred < (
int) preds_numu_only.size(); ++ipred){
471 double POT = preds_numu_only[ipred].pot;
473 auto fakedata = preds_numu_only[ipred].pred->Predict(calc);
474 auto hcos = preds[ipred].cos.first->ToTH1(preds[
i].
livetime, kBlack, kSolid,
kLivetime);
478 auto mockdata = fakedata.MockData(POT);
479 numu_expts.push_back(
new SingleSampleExperiment(preds_numu_only[ipred].
pred, mockdata, *preds_numu_only[ipred].
cos.first, preds_numu_only[ipred].cos.second));
485 if (
plot !=
"th13") {
489 std::cout <<
"Creating multiexperiment from " << expts.size()
499 preds_numu_only.clear();
502 std::cout <<
"\n\nSystematic correlations...\n";
503 if(!nueOnly && ! numuOnly && corrSysts){
521 if (nueOnly && corrSysts){
554 double maxmixing = .514;
563 fitnumu_only.
Fit(calc, auxShifts,
569 double pre_seed_dmsq = dmsq_var->
GetValue(calc);
572 double pre_seed_th23_LO = ( pre_seed_th23 > maxmixing ) ?
573 (2*maxmixing - pre_seed_th23) : pre_seed_th23;
574 double pre_seed_th23_UO = ( pre_seed_th23 > maxmixing ) ?
575 pre_seed_th23 : (2*maxmixing - pre_seed_th23);
576 if(pre_seed_th23_LO > .5) pre_seed_th23_LO = .45;
578 std::cout <<
"Preseeds ---> ssth23_LO = " << pre_seed_th23_LO;
586 std::vector<const IFitVar*> trueVars;
587 if (
plot==
"deltassth23"){
588 trueVars.push_back(dmsq_var);
591 if (
plot==
"ssth23dmsq32"){
600 double tb_ssth23 = 0;
605 std::vector <double> tf_delta_seeds = {0, .5, 1., 1.5};
606 std::vector <double> tf_ssth23_seeds = {pre_seed_th23_LO,
611 std::map<const IFitVar*, std::vector<double>> tf_seedPts;
612 if(
plot ==
"deltassth23")
613 tf_seedPts = {{dmsq_var, {nh?
fabs(pre_seed_dmsq) : -1*
fabs(pre_seed_dmsq)}},
615 else if(
plot ==
"ssth23dmsq32")
618 else if(
plot ==
"deltath13")
619 tf_seedPts = {{dmsq_var, {nh?
fabs(pre_seed_dmsq) : -1*
fabs(pre_seed_dmsq)}},
630 double trueChi = fittrue.
Fit(calc,
642 std::cout <<
"True fit parameters (post fit): ssth23 = " << tb_ssth23;
653 std::vector <double> bf_dmsq_seeds = {
fabs(pre_seed_dmsq),
654 -1*
fabs(pre_seed_dmsq)};
655 std::vector <double> bf_delta_seeds = {0., 0.5, 1., 1.5};
656 std::vector <double> bf_th23_seeds = {pre_seed_th23_LO, pre_seed_th23_UO,
658 std::vector <double> bf_th13_seeds = {0.082};
659 std::map <const IFitVar*, std::vector<double>> bf_seedPts =
665 double bestdelta = 0;
666 double bestssth23 = 0;
672 double bestChi = fit.Fit(calc,
680 bestth13 = kFitSinSq2Theta13.GetValue(calc);
682 std::cout <<
"Best fit parameters: ssth23 = " << bestssth23;
691 FCPoint pt(seeddelta, seedssth23, seedss2th13, seeddmsq32,
692 pre_seed_th23, pre_seed_dmsq,
693 tb_delta, tb_ssth23, tb_th13, tb_dmsq,
695 bestdelta, bestssth23, bestth13, bestdmsq,
701 unsigned int nfc_pts = fccol.
NPoints();
702 if(npts_flush != 0) {
703 if(nfc_pts % npts_flush == 0) {
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()
const FitDmSq32 kFitDmSq32
const FitDmSq32ScaledNH kFitDmSq32ScaledNH
Simple record of shifts applied to systematic parameters.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
void AddPoint(const FCPoint &p)
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)
void make_fc_surfaces_2020(int NPts, int startidx, int endidx, bool nh, int N, int jid, int npts_flush=0, std::string plot="deltassth23", std::string options="joint_fake2019", std::string decomp="", TString beam="both", bool corrSysts=false, bool readFromFile=false, bool runOnGrid=false)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var's SetValue()
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Encapsulate code to systematically shift a caf::SRProxy.
unsigned int NPoints() const
Represents the results of a single Feldman-Cousins pseudo-experiment.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
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 ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
void SetSeeds(osc::IOscCalcAdjustable *calc, double delta, double th23, double dmsq)
std::vector< float > Spectrum
std::vector< double > POT
Combine multiple component experiments.
const FitDmSq32ScaledIH kFitDmSq32ScaledIH
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const =0
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
const ReactorExperiment * WorldReactorConstraint2019()
Reactor constraint from PDG 2019.
const FitSinSq2Theta13 kFitSinSq2Theta13
std::string to_string(ModuleType mt)
void SaveToFile(const std::string &fname) const
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
void SetShift(const ISyst *syst, double shift, bool force=false)
void SetSystCorrelations(int idx, const std::vector< std::pair< const ISyst *, const ISyst * >> &corrs)
Compare a single data spectrum to the MC + cosmics expectation.
void plot(std::string label, std::map< std::string, std::map< std::string, Spectrum * >> &plots, bool log)
Collection of FCPoint. Serializable to/from a file.
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const =0
Perform MINUIT fits in one or two dimensions.