34 double dmsq32 = hier ? 2.44e-3 : -2.44e-3;
40 double& bestChi,
double& bestX,
double& bestY,
44 if (currentChiSq < bestChi){
45 bestChi = currentChiSq;
57 assert(bin >= 0 &&
"Invalid bin number");
81 Spectrum cosmicnue(hcosmicnue.first,pot,livetime);
84 std::vector<const IPrediction*> ps_numu =
86 std::vector< std::pair<TH1D*,double> > hcosmics_nm =
88 std::vector<Spectrum> cosmics_nm;
89 for (std::pair<TH1D*,double>
h : hcosmics_nm)
90 cosmics_nm.push_back(
Spectrum(
h.first,pot,livetime));
98 if (plot ==
"ssth23dmsq32"){
103 double widthX = (maxX-minX)/NX;
108 if (plot ==
"ssth23dmsq32"){
109 minY = nh ? 0.002 : -0.003;
110 maxY = nh ? 0.003 : -0.002;
113 double widthY = (maxY-
minY)/NY;
115 assert(bin < NX*NY &&
"Invalid bin number");
118 int binX = bin - NX*binY;
120 double centerX = minX + (binX+0.5)*widthX;
121 double centerY = minY + (binY+0.5)*widthY;
127 helperupsname+
"/contours_"+plot+
"_FCInput_2017.root";
128 TFile *fchelp = TFile::Open(helpername.c_str());
132 double seeddmsq32 = 0;
133 if (plot ==
"deltassth23"){
134 TH2* hdmsq = (TH2F*)fchelp->Get((hierStr+
"_DmSq32").c_str());
135 seeddmsq32 = hdmsq->GetBinContent(binX+1,binY+1);
138 double seeddelta = 0;
139 if (plot ==
"ssth23dmsq32"){
140 TH2* hdelta = (TH2F*)fchelp->Get((hierStr+
"_DeltaInPiUnits").c_str());
141 seeddelta = hdelta->GetBinContent(binX+1,binY+1);
144 TH2* hss2th13 = (TH2F*)fchelp->Get((hierStr+
"_SinSq2Theta13").c_str());
145 double seedss2th13 = hss2th13->GetBinContent(binX+1,binY+1);
147 std::map<std::string,double> seedsyst;
148 for (
const ISyst* syst : systs){
149 TH2*
h = (TH2F*)fchelp->Get((hierStr+
"_"+syst->ShortName()).c_str());
151 std::cout <<
"Don't see a prof hist for " << syst->ShortName() <<
". Continuing, but check your ups version." <<
std::endl;
154 double seedval = h->GetBinContent(binX+1,binY+1);
155 seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
161 for (
int i = 0;
i < NPts;
i++){
164 if (100*
i % NPts == 0){
170 while (trueX < minX || trueX > maxX)
171 trueX = rnd.Uniform(centerX-widthX/2,centerX+widthX/2);
173 while (trueY < minY || trueY > maxY)
174 trueY = rnd.Uniform(centerY-widthY/2,centerY+widthY/2);
178 kFitDmSq32. SetValue(calc, plot==
"deltassth23"?seeddmsq32:trueY);
182 for (
const ISyst *syst : snue){
183 if (seedsyst.find(syst->ShortName()) == seedsyst.end())
continue;
184 seednue.
SetShift(syst,seedsyst[syst->ShortName()]);
187 for (
const ISyst *syst : snumu){
188 if (seedsyst.find(syst->ShortName()) == seedsyst.end())
continue;
189 seednumu.
SetShift(syst,seedsyst[syst->ShortName()]);
197 fakenue += cosmicnue;
200 int mock_integral = mocknue.
Integral(pot);
203 std::vector<Spectrum> mocknumu;
204 assert(hcosmics_nm.size()==4 &&
"Issue with N quantiles in cosmics");
205 assert(ps_numu.size()==4 &&
"Issue with N quantiles in predictions");
206 for (
int i = 0;
i < 4;
i++){
207 Spectrum fakenm = ps_numu[
i]->PredictSyst(calc,seednumu);
208 fakenm += cosmics_nm[
i];
209 mocknumu.push_back(fakenm.
MockData(pot));
214 std::vector<const IExperiment*> expts;
217 for (
int i = 0;
i < 4;
i++)
219 hcosmics_nm[i].first,
230 std::vector<const IFitVar*> trueVars;
231 if (plot==
"deltassth23"){
235 if (plot==
"ssth23dmsq32"){
243 double trueChi = 1e6;
246 for (
double mecq0: {-1,+1}){
254 double dcp = calc->
GetdCP();
261 double bestChi = 1e6;
263 int bestHie = 0;
int bestDel = 0;
int bestOct = 0;
265 for (
bool hier : {
false,
true}){
266 for (
bool oct : {
false,
true}){
267 for (
int i = 0; i <=1; i++){
268 for (
double mecq0 : {-1,+1}){
274 bool isBest =
RunFitter(fit, calc, bestChi, bestX, bestY,
278 bestHie=hier; bestDel=
i; bestOct=oct;
285 if (trueChi < bestChi){
287 std::cout <<
"Pseudo-experiment bestChi, " << bestChi <<
" is less than the trueChi, " << trueChi <<
std::endl;
293 double fitloc = mock_integral*1e3+bestOct*1
e2+bestHie*1
e1+bestDel;
297 FCPoint pt(trueX, trueY, trueChi, bestX, bestY, bestChi, fitloc);
299 if (100*i % NPts == 0){
303 std::cout <<
"ChiSqaures: chitrue = " << trueChi <<
", chibest = " << bestChi <<
", deltachisquare = " << trueChi-bestChi <<
std::endl;
312 fccol.
SaveToFile(
"FCCol_"+plot+
"_"+binname+hier+number+
".root");
void joint_fit_2017_make_fc_surf(int NPts, int bin, bool nh, int N, std::string plot)
Cuts and Vars for the 2020 FD DiF Study.
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
bool RunFitter(MinuitFitter fit, osc::IOscCalcAdjustable *calc, double &bestChi, double &bestX, double &bestY, SystShifts seed=kNoShift)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var's GetValue()
const FitDmSq32 kFitDmSq32
Simple record of shifts applied to systematic parameters.
void AddPoint(const FCPoint &p)
std::vector< double > Spectrum
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var's SetValue()
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Encapsulate code to systematically shift a caf::SRProxy.
virtual void SetDmsq32(const T &dmsq32)=0
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
Representation of a spectrum in any variable, with associated POT.
Represents the results of a single Feldman-Cousins pseudo-experiment.
const double kAna2017Livetime
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::string getenv(std::string const &name)
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 Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const SystShifts kNoShift
Combine multiple component experiments.
const NOvARwgtSyst kMECq0ShapeSyst2017("MECq0Shape","MEC q_{0} shape", novarwgt::kMECq0ShapeSyst2017)
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
Standard interface to all prediction techniques.
bool SetHierarchy(osc::IOscCalcAdjustable *calc, bool hier)
const FitSinSq2Theta13 kFitSinSq2Theta13
T min(const caf::Proxy< T > &a, T b)
void SaveToFile(const std::string &fname) const
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)
void SetShift(const ISyst *syst, double shift, bool force=false)
std::string to_string(ModuleType const mt)
void SetSystCorrelations(int idx, const std::vector< std::pair< const ISyst *, const ISyst * >> &corrs)
virtual void SetdCP(const T &dCP)=0
Compare a single data spectrum to the MC + cosmics expectation.
Collection of FCPoint. Serializable to/from a file.
Perform MINUIT fits in one or two dimensions.