28 const std::vector<double>&
den,
38 Spectrum*
data = LoadFromFile<Spectrum>(
"$NOVA_ANA/steriles/Ana01/fout_nus_ana01_box_opening_restricted.root",
"spec_nus_cale_numi").
release();
41 <<
" and " << nData <<
" events" <<
std::endl;
44 Spectrum* cosmic_scale = LoadFromFile<Spectrum>(
"$NOVA_ANA/steriles/Ana01/AnaResults.root",
"cosmic__CalE").
release();
52 TFile* fileMC =
new TFile(
"/nova/app/users/gsdavies/nus/ana.a/NuXAna/macros/nus_pred_osc_syst_incl_ana01.root");
54 std::unique_ptr<PredictionInterp> pred123b =
57 std::unique_ptr<PredictionInterp> pred3c3d =
62 const double pot123b =
71 {&(*pred123b), pot123b},
72 {&(*pred3c3d), pot3c3d}
83 TH2D* hFCSig =
new TH2D(
"hFCSig", title.c_str(),
84 nbins34, min34, max34,
85 nbins24, min24, max24);
86 TH2D* hFCChi =
new TH2D(
"hFCChi", title.c_str(),
87 nbins34, min34, max34,
88 nbins24, min24, max24);
91 calc4f->SetNFlavors(4);
92 calc4f->SetDm(4, 0.5);
93 calc4f->SetAngle(2, 4, 0.);
94 calc4f->SetAngle(3, 4, 0.);
96 std::vector<const IFitVar*> fitvars_num;
100 std::vector<const IFitVar*> fitvars_den;
106 std::vector<double> rank_den;
109 for(
int i = 0;
i <= nData; ++
i) {
116 h->SetBinContent(1,
i);
126 for(
const auto& syst : systs) {
136 rank_den.push_back(TMath::Poisson(
i, nPred));
143 for(
int i34 = 1, n34 = hFCSig->GetNbinsX(); i34 <= n34; ++i34) {
146 for(
int i24 = 1, n24 = hFCSig->GetNbinsY(); i24 <= n24; ++i24) {
150 calc4f->SetAngle(3, 4, hFCSig->GetXaxis()->GetBinCenter(i34)*
M_PI/180.);
151 calc4f->SetAngle(2, 4, hFCSig->GetYaxis()->GetBinCenter(i24)*
M_PI/180.);
155 std::vector<double> rank_num;
164 for(
const auto& syst : systs) {
175 for(
int i = 0;
i <= nData; ++
i) {
177 rank_num.push_back(TMath::Poisson(
i, nPred));
184 double bin_significance =
CalcSignificance(rank_num, rank_den, nData, nPred);
192 hFCSig->SetBinContent(i34,i24, bin_significance);
193 hFCChi->SetBinContent(i34,i24, sigma*sigma);
197 TFile* rootF =
new TFile(outfile.c_str(),
"RECREATE");
199 rootF->WriteTObject(hFCSig);
200 rootF->WriteTObject(hFCChi);
202 TCanvas* cSig =
new TCanvas(
"cFCSig",
"cFCSig", 800, 500);
204 hFCSig->Draw(
"colz");
208 gPad->Print(
"FC2DSignificance_profd24.png");
209 rootF->WriteTObject(cSig);
211 TCanvas* cChi =
new TCanvas(
"cFCChi",
"cFCChi", 800, 500);
213 hFCChi->Draw(
"colz");
217 gPad->Print(
"FC2DChiSq_profd24.png");
218 rootF->WriteTObject(cChi);
222 std::cout <<
"Finished. Ctrl+Z should be safe if the macro didn't quit." <<
std::endl;
227 const std::vector<double>&
den,
231 double prob_tot = 0.;
233 std::map<double, std::vector<int> >
rank;
235 assert(num.size() == den.size());
236 for(
int i = 0,
n = num.size();
i <
n; ++
i) {
237 double r = num[
i]/den[
i];
239 if(rank.find(r) == rank.end()) {
240 rank[
r] = std::vector<int>();
243 rank[
r].push_back(
i);
246 double rData = num[nData]/den[nData];
249 std::map<double, std::vector<int> >::reverse_iterator top_rank = rank.rbegin();
250 double rCurr = top_rank->first;
252 while(rCurr >= rData && top_rank != rank.rend()) {
254 for(
const auto& nObs : top_rank->second) {
255 prob_tot += TMath::Poisson(nObs, nPred);
259 if(rank.size() == 0) {
break; }
260 top_rank = rank.rbegin();
261 rCurr = top_rank->first;
int rank(const C &v, int s)
Cuts and Vars for the 2020 FD DiF Study.
const double kSecondAnaPeriod2POT
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
A simple Gaussian constraint on an arbitrary IFitVar.
double Integral(const Spectrum &s, const double pot, int cvnregion)
Simple record of shifts applied to systematic parameters.
const double kSecondAnaEpoch3dPOT
const double kSecondAnaEpoch3bPOT
Adapt the PMNS_Sterile calculator to standard interface.
const FitTheta23InDegreesSterile kFitTheta23InDegreesSterile
const double kSecondAnaPeriod1POT
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Representation of a spectrum in any variable, with associated POT.
const XML_Char const XML_Char * data
double PValueToSigma(double p)
Compute the equivalent number of gaussian sigma for a p-value.
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
std::vector< const ISyst * > getAllNusSysts()
Get a vector of all the nus group systs.
double CalcSignificance(const std::vector< double > &num, const std::vector< double > &den, int nData, double nPred)
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 FCContour(std::string outfile)
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
Combine multiple component experiments.
void SetAngle(int i, int j, double th)
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir, const std::string &name)
void SetDm(int i, double dm)
assert(nhit_max >=nhit_nbins)
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
const double kSecondAnaEpoch3cPOT
Sum MC predictions from different periods scaled according to data POT targets.
const double kSecondAnaPOT
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
void ResetAngles(osc::OscCalcSterile *calc)
void SetShift(const ISyst *syst, double shift, bool force=false)
Compare a data spectrum to MC expectation using only the event count.
Perform MINUIT fits in one or two dimensions.