17 #include "CAFAna/Analysis/CovMxSurface.h" 20 #include "CAFAna/Analysis/Fit.h" 43 std::vector<std::string> polarity = {
"fhc",
"rhc" };
45 std::vector<std::string>
detector = {
"neardet",
"fardet" };
50 std::vector<const ISyst*> isysts;
51 if (opt.Contains(
"covmxsysts") or opt.Contains(
"pullsyst")) {
59 size_t pspos = optstring.find(
"pullsyst=");
60 if (pspos != std::string::npos) {
61 if (optstring.substr(pspos+8, 1) !=
"=")
62 assert(
false and
"The dm41 option must take the form \"dm41=xxx\".");
63 size_t start = pspos+9;
64 size_t end = optstring.find(
"_", pspos);
65 pullSyst = optstring.substr(start, end-start);
68 if (pullSyst ==
"all" and covmxSysts)
assert (
false and
69 "Treating systematics as both pull terms and covariance matrix is not permitted.");
73 size_t universe = 999999;
74 size_t fdpos = optstring.find(
"fakedata=");
75 if (optstring.substr(fdpos+8, 1) !=
"=")
76 assert(
false and
"The data option must take the form \"fakedata=xxx\".");
77 if (fdpos != std::string::npos) {
78 size_t start = fdpos + 9;
79 if (optstring.substr(start, 6) ==
"asimov") {
82 else if (optstring.substr(start, 5) ==
"stats") {
84 size_t end = optstring.find(
"_", fdpos);
85 universe = std::stoi(optstring.substr(start+5, end-(start+5)));
87 else if (optstring.substr(start, 5) ==
"systs") {
89 size_t end = optstring.find(
"_", fdpos);
90 universe = std::stoi(optstring.substr(start+5, end-(start+5)));
93 assert(
false and
"Fake data type not recognised!");
99 size_t dmpos = optstring.find(
"dm41=");
100 if (dmpos != std::string::npos) {
101 if (optstring.substr(dmpos+4, 1) !=
"=")
102 throw std::runtime_error(
"The dm41 option must take the form \"dm41=xxx\".");
103 size_t start = dmpos+5;
104 size_t end = optstring.find(
"_", dmpos);
105 if (end != std::string::npos)
106 dm41 =
std::stod(optstring.substr(start, end-start));
108 dm41 =
std::stod(optstring.substr(start));
113 bool profDelta24 = not
CheckOption(opt,
"noprofdelta24");
114 bool profVars = not
CheckOption(opt,
"noprofvars");
117 bool th24vsth34 =
false;
118 bool th24vsdm41 =
false;
119 bool th34vsdm41 =
false;
120 if (
CheckOption(opt,
"th24vsth34")) th24vsth34 =
true;
121 else if (
CheckOption(opt,
"th24vsdm41")) th24vsdm41 =
true;
122 else if (
CheckOption(opt,
"th34vsdm41")) th34vsdm41 =
true;
123 else throw std::runtime_error(
"The options string must contain plot type: \"th24vsth34\", \"th24vsdm41\" or \"th34vsdm41\".");
126 if (th24vsth34 && dm41 == -1)
127 throw std::runtime_error(
"If you want to run a th24 vs th34 surface, you must specify a value for dm41.");
140 std::vector<covmx::Sample> samples;
141 std::vector<IPrediction*> preds;
142 for (
size_t d = 0;
d < 2; ++
d) {
143 for (
size_t p = 0;
p < 2; ++
p) {
144 if (!usePolarity[
p] || !useDetector[
d])
continue;
152 std::vector<Spectrum*>
cosmic;
153 for (
auto sample : samples)
157 std::vector<Spectrum*>
data;
158 if (fakedata ==
"asimov") {
160 for (
size_t i = 0;
i < preds.size(); ++
i) {
164 data.back()->OverrideLivetime(livetimeFD);
165 if (cosmic[
i]) *data[
i] += *cosmic[
i];
169 std::cout <<
"Loading universe " << universe <<
" of " << fakedata
175 std::vector<const ISyst*> pullSysts = {};
176 SystConcat* sc =
nullptr;
177 if (pullSyst ==
"all") {
179 if (samples.size() == 1)
184 pullSysts = sc->GetActiveSysts();
187 else if (pullSyst !=
"none") {
189 if (samples.size() == 1) {
190 for (
const ISyst* syst : isysts) {
192 pullSysts = { syst };
196 if (pullSysts.empty()) {
197 std::ostringstream err;
198 err <<
"Couldn't find requested pull systematic " << pullSyst;
199 assert (
false and err.str().c_str());
204 std::vector<SystConcat> all_scs =
GetSystConcats(samples, preds, isysts);
205 for (
auto it_sc : all_scs) {
206 if (
GetSystBaseName(it_sc.GetActiveSysts()[0]->ShortName()) == pullSyst) {
207 sc =
new SystConcat(it_sc);
208 pullSysts = sc->GetActiveSysts();
216 CovarianceMatrix* mx =
nullptr;
218 std::vector<Binning>
bins;
219 for (
auto d : data) bins.push_back(
d->GetBinnings()[0]);
220 mx =
new CovarianceMatrix(bins);
224 mx->AddMatrix(it_mx);
227 std::cout <<
"Omitting pull term " << pullSyst
242 if (samples.size() > 1 and sc) expt->ConfigureSysts(sc);
244 for (
int i = 1;
i <= 4; ++
i) {
245 noosccalc->SetDm(
i, 0);
247 expt->DisableOscillations({
true,
false}, noosccalc);
251 expts.push_back(expt);
255 if (th24vsth34) calc->
SetDm(4, dm41);
262 std::map<const IFitVar*, std::vector<double> > seedValues;
270 std::vector<const IFitVar*> xyVars;
278 else if (th24vsdm41) {
286 else if (th34vsdm41) {
311 preds, data, cosmic);
315 TFile*
f =
new TFile(outfile.c_str(),
"recreate");
316 surf.
SaveTo(f->mkdir(optstring.c_str()));
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
size_t JobNumber()
What's the process number for a grid job?
Represent the binning of a Spectrum's x-axis.
Cuts and Vars for the 2020 FD DiF Study.
std::string GetSystBaseName(std::string name)
Get the last part of a systematic's ShortName.
void SetNus20Params(osc::OscCalcSterile *calc, std::string type="3flav")
double stod(const std::string &val)
bool RunningOnGrid()
Is this a grid (condor) job?
void SetDelta(int i, int j, double delta)
Adapt the PMNS_Sterile calculator to standard interface.
const FitDmSq32Sterile kFitDmSq32Sterile
const double kAna2018SensitivityFHCNDPOT
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
static SystShifts Nominal()
Encapsulate code to systematically shift a caf::SRProxy.
void SaveTo(TDirectory *dir) const
Standard SaveTo implementation.
Representation of a spectrum in any variable, with associated POT.
const XML_Char const XML_Char * data
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Sum up livetimes from individual cosmic triggers.
std::vector< const IExperiment * > GetConstraints()
Gaussian constrains on atmospheric parameters.
std::vector< bool > GetOptionConfig(TString opt, std::vector< std::string > options)
Compare a single data spectrum to the MC + cosmics expectation.
const covmx::Sample kNusFHCFarDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kFarDet)
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
SystConcat GetSystConcat(std::vector< covmx::Sample > samples, std::vector< const ISyst * > systs, std::vector< IPrediction * > preds)
Get correctly concatenated systs for ND and FD.
bool CheckOption(TString opts, TString opt)
Combine multiple component experiments.
void SetDm(int i, double dm)
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
static Binning LogUniform(int n, double lo, double hi)
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
std::map< std::string, ana::Spectrum > LoadFakeData(const std::string &fakeDataFilename, const std::string &path)
const double kAna2018FHCPOT
Prevent histograms being added to the current directory.
void MakeSurfaceNoNDOsc(TString opt)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
const FitTheta23Sterile kFitTheta23Sterile