4 #include "TDirectory.h" 25 #include "CAFAna/Analysis/CovMxSurface.h" 40 double log_spacing = exp( (
log(xhi) -
log(xlo))/(nbins) );
41 std::vector<double> binning;
43 if (
i == 0) binning.push_back(xlo);
44 else binning.push_back(log_spacing*binning[
i-1]);
51 double lin_spacing = (xhi -
xlo)/(nbins);
52 std::vector<double> binning;
54 if (
i == 0) binning.push_back(xlo);
55 else binning.push_back(lin_spacing+binning[
i-1]);
61 int N = 1,
int jid = 1,
int npts_flush = 1,
63 TString
beam =
"neardetfardetFHCcovmx",
bool corrSysts =
true,
64 bool readFromFile =
true,
double max_time_in_sec = 3600)
67 std::clock_t tstart = std::clock();
73 std::vector<std::string> polarity = {
"fhc",
"rhc" };
75 std::vector<std::string>
detector = {
"neardet",
"fardet" };
78 bool covmx =
beam.Contains(
"covmx", TString::kIgnoreCase);
79 bool extrap =
beam.Contains(
"extrap", TString::kIgnoreCase);
82 if (covmx && use_polarity[1])
83 throw std::runtime_error(
"RHC covariance matrix not supported yet.");
85 if (extrap && use_polarity[0])
86 throw std::runtime_error(
"Extrapolated fit not supported for FHC.");
90 if (covmx && corrSysts) {
92 std::cout <<
"Using covariance matrix, so disabling traditional systematics." <<
std::endl;
95 bool upperOctant =
false;
96 if(
plot ==
"th24vsth34_UO" ||
99 plot ==
"dm41vsth24_UO" ||
100 plot ==
"dm41vsth34_UO")
107 bool th24vsth34 =
false;
108 bool dm41vsth24 =
false;
109 bool dm41vsth34 =
false;
110 if(
plot ==
"th24_UO" ||
plot ==
"th24_LO"){
113 else if(
plot ==
"th34_UO" ||
plot ==
"th34_LO"){
116 else if(
plot ==
"th24vsth34_UO" ||
plot ==
"th24vsth34_LO") th24vsth34 =
true;
117 else if(
plot ==
"dm41vsth24_UO" ||
plot ==
"dm41vsth24_LO") dm41vsth24 =
true;
118 else if(
plot ==
"dm41vsth34_UO" ||
plot ==
"dm41vsth34_LO") dm41vsth34 =
true;
121 if(corrSysts)
beam +=
"syst";
131 std::vector<covmx::Sample> samples;
132 std::vector<IPrediction*> preds;
133 std::vector<Spectrum*>
cosmic;
134 std::vector<Spectrum*>
data;
135 std::vector<Spectrum*>::iterator iter;
137 CovarianceMatrix* mx =
nullptr;
140 std::vector<const ISyst*>
systs = {};
141 SystConcat* sc =
nullptr;
149 for (
size_t d = 0;
d < 2; ++
d) {
150 for (
size_t p = 0;
p < 2; ++
p) {
151 if (!use_polarity[
p] || !use_detector[
d])
continue;
154 preds.push_back(
LoadPrediction(samples.back(), pull_syst ==
"none" ?
"nosysts" :
"systs"));
158 for (
auto sample : samples)
163 if (pull_syst ==
"all") {
165 if (samples.size() == 1)
170 systs = sc->GetActiveSysts();
173 else if (pull_syst !=
"none") {
175 if (samples.size() == 1) {
176 for (
const ISyst* syst : isysts) {
183 std::ostringstream err;
184 err <<
"Couldn't find requested pull systematic " << pull_syst;
185 assert (
false and err.str().c_str());
190 std::vector<SystConcat> all_scs =
GetSystConcats(samples, preds, isysts);
191 for (
auto it_sc : all_scs) {
192 if (
GetSystBaseName(it_sc.GetActiveSysts()[0]->ShortName()) == pull_syst) {
193 sc =
new SystConcat(it_sc);
194 systs = sc->GetActiveSysts();
201 for (
size_t i = 0;
i < preds.size(); ++
i) {
205 if (cosmic[
i]) *data[
i] += *cosmic[
i];
208 std::vector<Binning>
bins;
209 for (
auto d : data) bins.push_back(
d->GetBinnings()[0]);
210 mx =
new CovarianceMatrix(bins);
212 for (CovarianceMatrix* it_mx :
GetMatrices(samples)) {
214 std::cout <<
"Adding " << it_mx->GetName()
215 <<
" to covariance matrix." <<
std::endl;
216 mx->AddMatrix(it_mx);
219 std::cout <<
"Omitting pull term " << pull_syst
244 std::vector<double> th24Bins;
245 std::vector<double> th34Bins;
246 std::vector<double> dm2Bins;
249 th24Bins =
LinBins(180, 0., 45.);
252 th34Bins =
LinBins(360, 0., 90.);
254 else if (th24vsth34) {
255 th24Bins =
LinBins(45, 0., 45.);
256 th34Bins =
LinBins(90, 0., 90.);
258 else if (dm41vsth24) {
262 else if (dm41vsth34) {
267 double binxlowedge = 0.;
268 double binxuppedge = 0.;
269 double binylowedge = 0.;
270 double binyuppedge = 0.;
271 double binxcenter = 0.;
272 double binycenter = 0.;
277 std::ifstream binlist;
280 binlist.open(binlist_fname);
292 }
while(!binlist.eof());
295 else for (
int i = 0;
i < TotBins;
i++) bins[
i]=
i;
297 int ibin = bins[startidx];
298 int fbin = bins[endidx];
304 fcout += use_polarity[0] ?
"FHC":
"RHC";
305 fcout += corrSysts ?
"_systs":
"_stats";
309 std::clock_t tload = std::clock();
310 double loadtime = ( tload - tstart ) / (
double) CLOCKS_PER_SEC;
314 for (
int i = startidx;
i <= endidx;
i++) {
316 int iter_bin = bins[
i];
319 binxlowedge = th34Bins.front();
320 binxuppedge = th34Bins.back();
321 binylowedge = th24Bins[iter_bin];
322 binyuppedge = th24Bins[iter_bin+1];
323 binxcenter = 0.5*(binxlowedge+binxuppedge);
324 binycenter = 0.5*(binylowedge+binyuppedge);
327 binxlowedge = th34Bins[iter_bin];
328 binxuppedge = th34Bins[iter_bin+1];
329 binylowedge = th24Bins.front();
330 binyuppedge = th24Bins.back();
331 binxcenter = 0.5*(binxlowedge+binxuppedge);
332 binycenter = 0.5*(binylowedge+binyuppedge);
334 else if (th24vsth34) {
335 int x_bins = th34Bins.size() - 1;
336 int x_idx = iter_bin % x_bins;
337 int y_idx = iter_bin / x_bins;
338 binxlowedge = th34Bins[x_idx];
339 binxuppedge = th34Bins[x_idx+1];
340 binylowedge = th24Bins[y_idx];
341 binyuppedge = th24Bins[y_idx+1];
342 binxcenter = 0.5*(binxlowedge+binxuppedge);
343 binycenter = 0.5*(binylowedge+binyuppedge);
345 else if (dm41vsth24) {
346 int x_bins = th24Bins.size() - 1;
347 int x_idx = iter_bin % x_bins;
348 int y_idx = iter_bin / x_bins;
349 binxlowedge = th24Bins[x_idx];
350 binxuppedge = th24Bins[x_idx+1];
351 binylowedge = dm2Bins[y_idx];
352 binyuppedge = dm2Bins[y_idx+1];
353 binxcenter = 0.5*(binxlowedge+binxuppedge);
354 binycenter = 0.5*(binylowedge+binyuppedge);
356 else if (dm41vsth34) {
357 int x_bins = th34Bins.size() - 1;
358 int x_idx = iter_bin % x_bins;
359 int y_idx = iter_bin / x_bins;
360 binxlowedge = th34Bins[x_idx];
361 binxuppedge = th34Bins[x_idx+1];
362 binylowedge = dm2Bins[y_idx];
363 binyuppedge = dm2Bins[y_idx+1];
364 binxcenter = 0.5*(binxlowedge+binxuppedge);
365 binycenter = 0.5*(binylowedge+binyuppedge);
368 std::cout <<
"Bin X (low-center-up): " << binxlowedge <<
"-" << binxcenter <<
"-" << binxuppedge <<
std::endl;
369 std::cout <<
"Bin Y (low-center-up): " << binylowedge <<
"-" << binycenter <<
"-" << binyuppedge <<
std::endl;
371 double profDelta13_PiUnits = 1.46;
373 double profssTh23 = 0.470;
375 double profDelta24_PiUnits = 1.;
380 double profTh24 = 0.;
381 double profTh34 = 0.;
382 double profDm41 = 0.;
383 if(th24 || th34 || th24vsth34) {
384 profTh24 = binycenter;
385 profTh34 = binxcenter;
388 else if(dm41vsth24) {
389 profTh24 = binxcenter;
391 profDm41 = binycenter;
393 else if(dm41vsth34) {
395 profTh34 = binxcenter;
396 profDm41 = binycenter;
400 std::vector<const IFitVar*> prof_vars;
411 else if(dm41vsth24) {
414 else if(dm41vsth34) {
418 std::vector<const IFitVar*> fitvars = prof_vars;
426 else if(th24vsth34) {
430 else if(dm41vsth24) {
434 else if(dm41vsth34) {
439 for(
int iTrial = 0; iTrial < NPts; ++iTrial){
441 std::clock_t tistart = std::clock();
442 if (max_time_in_sec > 0.0){
443 double elapstime = ( tistart - tstart ) / (
double) CLOCKS_PER_SEC;
444 if (elapstime >= max_time_in_sec)
break;
451 auto trueX = rnd.Uniform(binxlowedge, binxuppedge);
452 auto trueY = rnd.Uniform(binylowedge, binyuppedge);
453 auto trueDelta24 = rnd.Uniform(0.,2.);
455 double trueTh24 = 0.;
456 double trueTh34 = 0.;
457 double trueDm41 = 0.;
462 if(th24 || th34 || th24vsth34) {
469 else if(dm41vsth24) {
472 trueTh34 = rnd.Uniform(0.,1.);
476 else if(dm41vsth34) {
479 trueTh24 = rnd.Uniform(0.,1.);
489 for (iter=data.begin(); iter != data.end(); ++iter){
497 for (
size_t i = 0;
i < preds.size(); ++
i) {
498 data.push_back(
new Spectrum(preds[
i]->Predict(calc).MockData(12e20)));
501 if (cosmic[
i]) *data[
i] += *cosmic[
i];
526 expts.push_back(expt);
546 else if(dm41vsth24) {
549 else if(dm41vsth34) {
556 double bestChi = 1E20;
557 double bestDelta13 = 0;
558 double bestssTh13 = 0;
559 double bestssTh23 = 0;
561 double bestDelta24 = 0;
566 if(th24 || th34 || th24vsth34) {
567 for (
double th34_seed: {15, 30, 45, 60, 75}) {
568 for (
double th24_seed: {9, 18, 27, 36}) {
569 std::vector<double> th34_seed_vec; th34_seed_vec.push_back(th34_seed);
570 std::vector<double> th24_seed_vec; th24_seed_vec.push_back(th24_seed);
571 std::map<const IFitVar*, std::vector<double> > seedPts;
576 std::vector<SystShifts> seedShifts = {};
579 auto bestChi_mtplSeed = fit.
Fit(calc, auxShifts, seedPts, seedShifts,
IFitter::kQuiet)->EvalMetricVal();
581 if (bestChi_mtplSeed < bestChi) {
582 bestChi = bestChi_mtplSeed;
599 else if(dm41vsth24 || dm41vsth34) {
600 for (
double th_seed: {1.0e-5, 1.0e-4, 1.0e-3, 1.0e-2, 1.0e-1}) {
601 for (
double dm_seed: {5.0e-2, 5.0e-1, 5.0, 5.0e+1}) {
602 std::vector<double> th_seed_vec; th_seed_vec.push_back(th_seed);
603 std::vector<double> dm_seed_vec; dm_seed_vec.push_back(dm_seed);
604 std::map<const IFitVar*, std::vector<double> > seedPts;
610 else if(dm41vsth34) {
617 std::vector<SystShifts> seedShifts = {};
620 auto bestChi_mtplSeed = fit.
Fit(calc, auxShifts, seedPts, seedShifts,
IFitter::kQuiet)->EvalMetricVal();
622 if (bestChi_mtplSeed < bestChi) {
623 bestChi = bestChi_mtplSeed;
640 std::clock_t tiend = std::clock();
641 double calctime = ( tiend - tistart ) / (
double) CLOCKS_PER_SEC;
643 FCPoint pt(profDelta13_PiUnits, profssTh23, profssTh13, profDmSq,
644 profDelta24_PiUnits, profTh24, profTh34, profDm41,
646 trueDelta13, truessTh23, truessTh13, trueDmSq,
647 trueDelta24, trueTh24, trueTh34, trueDm41,
649 bestDelta13, bestssTh23, bestssTh13, bestDmSq,
650 bestDelta24, bestTh24, bestTh34, bestDm41,
652 trueChi-bestChi, loadtime, calctime);
657 if(th24 || th34 || th24vsth34) {
661 else if(dm41vsth24) {
665 else if(dm41vsth34) {
670 std::cout <<
"Trial: " << iTrial <<
", TrueX: " << trueX <<
", Y: " << trueY <<
", trueChi: " << trueChi <<
", bestX: " << bestX <<
", Y: " << bestY <<
", bestChi: " << bestChi <<
", true-best: " << trueChi-bestChi <<
", loadTime: " << loadtime <<
", calcTime: " << calctime <<
std::endl;
std::vector< SystConcat > GetSystConcats(std::vector< covmx::Sample > samples, std::vector< IPrediction * > preds, std::vector< const ISyst * > systs)
Get correctly concatenated systs for ND and FD.
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
IPrediction * LoadExtrapPrediction(TString opt)
Spectrum LoadCosmicSpectrum(bool rhc=false)
Function to return cosmic spectrum.
Cuts and Vars for the 2020 FD DiF Study.
std::string GetSystBaseName(std::string name)
Get the last part of a systematic's ShortName.
double GetValue(const osc::IOscCalcAdjustable *osc) const override
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Simple record of shifts applied to systematic parameters.
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN'T A BETTER WAY!
std::vector< const ISyst * > LoadISysts()
Function to load ISysts.
void AddPoint(const FCPoint &p)
std::vector< double > Spectrum
Adapt the PMNS_Sterile calculator to standard interface.
const FitDmSq32Sterile kFitDmSq32Sterile
std::vector< double > LogBins(int nbins, double xlo, double xhi)
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
static SystShifts Nominal()
void SetAngles(osc::OscCalcSterile *calc)
Encapsulate code to systematically shift a caf::SRProxy.
void Scale(double c)
Multiply this spectrum by a constant c.
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Representation of a spectrum in any variable, with associated POT.
const double kAna2018RHCPOT
const XML_Char const XML_Char * data
Represents the results of a single Feldman-Cousins pseudo-experiment.
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
const FitSinSqTheta13Sterile kFitSinSqTheta13Sterile
Sum up livetimes from individual cosmic triggers.
std::vector< const IExperiment * > GetConstraints()
Gaussian constrains on atmospheric parameters.
double GetValue(const osc::IOscCalcAdjustable *osc) 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.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::vector< bool > GetOptionConfig(TString opt, std::vector< std::string > options)
Compare a single data spectrum to the MC + cosmics expectation.
double GetValue(const osc::IOscCalcAdjustable *osc) const override
void OverrideLivetime(double newlive)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN'T A BETTER WAY!
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
const FitDelta13InPiUnitsSterile kFitDelta13InPiUnitsSterile
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
SystConcat GetSystConcat(std::vector< covmx::Sample > samples, std::vector< const ISyst * > systs, std::vector< IPrediction * > preds)
Get correctly concatenated systs for ND and FD.
std::vector< double > LinBins(int nbins, double xlo, double xhi)
Combine multiple component experiments.
void make_fc_nus_surfs_nersc_2019(int NPts=1, int startidx=0, int endidx=0, bool nh=true, int N=1, int jid=1, int npts_flush=1, std::string plot="th24vsth34_UO", TString beam="neardetfardetFHCcovmx", bool corrSysts=true, bool readFromFile=true, double max_time_in_sec=3600)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
const FitSinSqTheta23Sterile kFitSinSqTheta23Sterile
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Base class defining interface for experiments.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::unique_ptr< Spectrum > LoadCosmic(covmx::Sample sample, bool cvmfs=true)
Get cosmics for a given sample.
assert(nhit_max >=nhit_nbins)
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
double GetValue(const osc::IOscCalcAdjustable *osc) const override
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Standard interface to all prediction techniques.
const double kAna2018FHCPOT
void SaveToFile(const std::string &fname) const
Prevent histograms being added to the current directory.
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
std::string to_string(ModuleType const mt)
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
const double kNus19SensLivetime
const double kAna2018RHCLivetime
std::vector< const ISyst * > GetSystsFromFile(covmx::Sample sample)
Get systematics for a given sample.
Compare a single data spectrum to the MC + cosmics expectation.
const double kNus19Delta13
Collection of FCPoint. Serializable to/from a file.
Perform MINUIT fits in one or two dimensions.