FCContour.C
Go to the documentation of this file.
6 #include "CAFAna/Core/Spectrum.h"
14 #include "NuXAna/Systs/NusSysts.h"
16 
17 #include "OscLib/OscCalcSterile.h"
18 
19 #include "TFile.h"
20 #include "TCanvas.h"
21 #include "TH1.h"
22 #include "TH2D.h"
23 #include "TMath.h"
24 
25 using namespace ana;
26 
27 double CalcSignificance(const std::vector<double>& num,
28  const std::vector<double>& den,
29  int nData,
30  double nPred);
31 
32 void ResetAngles(osc::OscCalcSterile* calc); //should just make this globally somewhere
33 
35 {
36  // Get the data spectrum
37  // Load in the NuMI data spectrum
38  Spectrum* data = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Ana01/fout_nus_ana01_box_opening_restricted.root", "spec_nus_cale_numi").release();
39  int nData = data->Integral(data->POT());
40  std::cout << "Data spectrum: " << data->POT()
41  << " and " << nData << " events" << std::endl;
42 
43  // Load in the out-of-NuMI-time weighted spectrum (cosmics)
44  Spectrum* cosmic_scale = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Ana01/AnaResults.root", "cosmic__CalE").release();
45  TH1D* hCosmic = cosmic_scale->ToTH1(cosmic_scale->Livetime(), kLivetime);
46  Spectrum sCos(hCosmic, data->POT(), data->Livetime());
47  double nCos = sCos.ToTH1(data->POT())->Integral(1, hCosmic->GetNbinsX());
48  std::cout << "Total cosmic count: " << nCos << std::endl;
49 
50  // Load in the MC Prediction
51  // TFile* fileMC = new TFile("$NOVA_ANA/steriles/Ana01/nus_ana01_prediction.root", "READ");
52  TFile* fileMC = new TFile("/nova/app/users/gsdavies/nus/ana.a/NuXAna/macros/nus_pred_osc_syst_incl_ana01.root");
53 
54  std::unique_ptr<PredictionInterp> pred123b =
55  PredictionInterp::LoadFrom((TDirectory*)fileMC->Get("nus_pred_p1p2e3b"));
56 
57  std::unique_ptr<PredictionInterp> pred3c3d =
58  PredictionInterp::LoadFrom((TDirectory*)fileMC->Get("nus_pred_e3c"));
59 
60  fileMC->Close();
61 
62  const double pot123b =
66 
67  const double pot3c3d = kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT;
68 
69  // Combine period1/period2/epoch3b + epoch3c, weighted appropriately
71  {&(*pred123b), pot123b},
72  {&(*pred3c3d), pot3c3d}
73  });
74 
75  // Initialize FC histogram
76  std::string title = ";#theta_{34} (Deg.); #theta_{24} (Deg.)";
77  int nbins34 = 200;
78  double min34 = 0.;
79  double max34 = 50.;
80  int nbins24 = 180;
81  double min24 = 0.;
82  double max24 = 45.;
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);
89 
91  calc4f->SetNFlavors(4);
92  calc4f->SetDm(4, 0.5); // #Delta m^{2}_{41} = 0.5 eV^2
93  calc4f->SetAngle(2, 4, 0.); // #theta_{24} = 0 degrees
94  calc4f->SetAngle(3, 4, 0.); // #theta_{34} = 0 degrees
95 
96  std::vector<const IFitVar*> fitvars_num;
97  fitvars_num.push_back(&kFitTheta23InDegreesSterile);
98  fitvars_num.push_back(&kFitDelta24InPiUnitsSterile);
99 
100  std::vector<const IFitVar*> fitvars_den;
101  fitvars_den.push_back(&kFitTheta23InDegreesSterile);
102  fitvars_den.push_back(&kFitTheta24InDegreesSterile);
103  fitvars_den.push_back(&kFitTheta34InDegreesSterile);
104  fitvars_den.push_back(&kFitDelta24InPiUnitsSterile);
105 
106  std::vector<double> rank_den;
107 
108 // std::cout << "Rank denominator: " << std::endl;
109  for(int i = 0; i <= nData; ++i) {
110 
111  ResetAngles(calc4f);
112 
113  // Create a fake data spectrum that only contains i events
114  TH1* h = data->ToTH1(kSecondAnaPOT);
115  h->Reset();
116  h->SetBinContent(1, i);
117  Spectrum fake(h, data->POT(), data->Livetime());
118 
120  45.8, 3.2);
121  CountingExperiment nus_expt(&pred, fake, sCos);
122  MultiExperiment expt({&th23Constraint, &nus_expt});
123 
124  std::vector<const ISyst*> systs = getAllNusSysts();
125  SystShifts shifts;
126  for(const auto& syst : systs) {
127  shifts.SetShift(syst, 0);
128  }
129 
130  MinuitFitter fit(&expt, fitvars_den, systs);
131  double chi2 = fit.Fit(calc4f, shifts, IFitter::kQuiet)->EvalMetricVal();
132 
133  double nPred = pred.PredictSyst(calc4f, shifts).Integral(kSecondAnaPOT);
134  nPred += sCos.Integral(kSecondAnaPOT);
135 
136  rank_den.push_back(TMath::Poisson(i, nPred));
137  std::cout << "Data value: " << i << std::endl;
138 // std::cout << "Data value: " << i << ", " << nPred << ", "
139 // << chi2 << ", " << rank_den[i] << std::endl;
140  }
141  std::cout << "Finished calculating R denominator values." << std::endl;
142 
143  for(int i34 = 1, n34 = hFCSig->GetNbinsX(); i34 <= n34; ++i34) {
144  std::cout << "On loop over theta_34 = " << i34 << std::endl;
145 
146  for(int i24 = 1, n24 = hFCSig->GetNbinsY(); i24 <= n24; ++i24) {
147  std::cout << "On loop over theta_24 = " << i24 << std::endl;
148 
149  ResetAngles(calc4f);
150  calc4f->SetAngle(3, 4, hFCSig->GetXaxis()->GetBinCenter(i34)*M_PI/180.);
151  calc4f->SetAngle(2, 4, hFCSig->GetYaxis()->GetBinCenter(i24)*M_PI/180.);
152  std::cout << "t34: " << hFCSig->GetXaxis()->GetBinCenter(i34)*M_PI/180. << std::endl;
153  std::cout << "t24: " << hFCSig->GetXaxis()->GetBinCenter(i24)*M_PI/180. << std::endl;
154 
155  std::vector<double> rank_num;
156 
158  45.8, 3.2);
159  CountingExperiment nus_expt(&pred, *data, sCos);
160  MultiExperiment expt({&th23Constraint, &nus_expt});
161 
162  std::vector<const ISyst*> systs = getAllNusSysts();
163  SystShifts shifts;
164  for(const auto& syst : systs) {
165  shifts.SetShift(syst, 0);
166  }
167 
168  MinuitFitter fit(&expt, fitvars_num, systs);
169  double chi2 = fit.Fit(calc4f, shifts, IFitter::kQuiet)->EvalMetricVal();
170  std::cout << " chi2: " << chi2 << std::endl;
171  double nPred = pred.PredictSyst(calc4f, shifts).Integral(kSecondAnaPOT);
172  std::cout << "npred: " << nPred << std::endl;
173  nPred += sCos.Integral(kSecondAnaPOT);
174  std::cout << "new nPred: " << nPred << std::endl;
175  for(int i = 0; i <= nData; ++i) {
176  chi2 += LogLikelihood(nPred, i);
177  rank_num.push_back(TMath::Poisson(i, nPred));
178 // std::cout << "Data value: " << i << ", " << nPred << ", "
179 // << chi2 << ", " << rank_num[i] << ", "
180 // << rank_num[i]/rank_den[i] << std::endl;
181  chi2 -= LogLikelihood(nPred, i);
182  }
183  std::cout << "new chi2: " << chi2 << std::endl;
184  double bin_significance = CalcSignificance(rank_num, rank_den, nData, nPred);
185  std::cout << "bin_sig: " << bin_significance << std::endl;
186  // Converting to “Gaussian” sigma
187  const double sigma = PValueToSigma(1.0-bin_significance);
188 
189 // std::cout << "Bin significance: " << bin_significance << ", "
190 // << "Sigma: " << sigma << ", "
191 // << "Chi2: " << sigma*sigma << std::endl;
192  hFCSig->SetBinContent(i34,i24, bin_significance);
193  hFCChi->SetBinContent(i34,i24, sigma*sigma);
194  } // end of loop over theta 24 bins
195  } // end of loop over theta 34 bins
196 
197  TFile* rootF = new TFile(outfile.c_str(), "RECREATE");
198 
199  rootF->WriteTObject(hFCSig);
200  rootF->WriteTObject(hFCChi);
201 
202  TCanvas* cSig = new TCanvas("cFCSig", "cFCSig", 800, 500);
203  //CenterTitles(hFCSig);
204  hFCSig->Draw("colz");
205  gPad->RedrawAxis();
206  //Preliminary();
207 
208  gPad->Print("FC2DSignificance_profd24.png");
209  rootF->WriteTObject(cSig);
210 
211  TCanvas* cChi = new TCanvas("cFCChi", "cFCChi", 800, 500);
212  //CenterTitles(hFCChi);
213  hFCChi->Draw("colz");
214  gPad->RedrawAxis();
215  //Preliminary();
216 
217  gPad->Print("FC2DChiSq_profd24.png");
218  rootF->WriteTObject(cChi);
219 
220  rootF->Close();
221 
222  std::cout << "Finished. Ctrl+Z should be safe if the macro didn't quit." << std::endl;
223 }
224 
225 //------------------------------------------------------------------------------
226 double CalcSignificance(const std::vector<double>& num,
227  const std::vector<double>& den,
228  int nData,
229  double nPred)
230 {
231  double prob_tot = 0.;
232 
233  std::map<double, std::vector<int> > rank;
234 
235  assert(num.size() == den.size());
236  for(int i = 0, n = num.size(); i < n; ++i) {
237  double r = num[i]/den[i];
238 
239  if(rank.find(r) == rank.end()) {
240  rank[r] = std::vector<int>();
241  }
242 
243  rank[r].push_back(i);
244  }
245 
246  double rData = num[nData]/den[nData];
247  //std::cout << "rData: " << rData << std::endl;
248 
249  std::map<double, std::vector<int> >::reverse_iterator top_rank = rank.rbegin();
250  double rCurr = top_rank->first;
251 
252  while(rCurr >= rData && top_rank != rank.rend()) {
253 // std::cout << "rCurr: " << rCurr << std::endl;
254  for(const auto& nObs : top_rank->second) {
255  prob_tot += TMath::Poisson(nObs, nPred);
256  }
257 
258  rank.erase(rCurr);
259  if(rank.size() == 0) { break; }
260  top_rank = rank.rbegin();
261  rCurr = top_rank->first;
262  }
263 
264 // std::cout << "Rank Size after: " << rank.size() << std::endl;
265 
266  return prob_tot;
267 }
268 
269 //------------------------------------------------------------------------------
271 {
273  calc->SetDm(4, 0.5); // #deltam41 = 0.5 eV^2
274  calc->SetAngle(2, 3, M_PI/4);
275  calc->SetAngle(3, 4, 0.);
276  calc->SetAngle(2, 4, 0.);
277 }
int rank(const C &v, int s)
Definition: rank.hpp:20
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
const double kSecondAnaPeriod2POT
Definition: Exposures.h:74
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.
Definition: Spectrum.cxx:209
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.
Definition: SystShifts.h:20
Float_t den
Definition: plot.C:36
const double kSecondAnaEpoch3dPOT
Definition: Exposures.h:77
const double kSecondAnaEpoch3bPOT
Definition: Exposures.h:75
Adapt the PMNS_Sterile calculator to standard interface.
const FitTheta23InDegreesSterile kFitTheta23InDegreesSterile
const double kSecondAnaPeriod1POT
Definition: Exposures.h:73
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:333
osc::OscCalcDumb calc
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
double chi2()
#define M_PI
Definition: SbMath.h:34
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const XML_Char const XML_Char * data
Definition: expat.h:268
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.
Definition: Calcs.cxx:97
expt
Definition: demo5.py:34
std::vector< const ISyst * > getAllNusSysts()
Get a vector of all the nus group systs.
Definition: NusSysts.cxx:202
double CalcSignificance(const std::vector< double > &num, const std::vector< double > &den, int nData, double nPred)
Definition: FCContour.C:226
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.
Definition: IFitter.cxx:64
void FCContour(std::string outfile)
Definition: FCContour.C:34
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:76
double POT() const
Definition: Spectrum.h:231
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
void SetAngle(int i, int j, double th)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir, const std::string &name)
void SetDm(int i, double dm)
int num
Definition: f2_nu.C:119
assert(nhit_max >=nhit_nbins)
TRandom3 r(0)
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
const double kSecondAnaEpoch3cPOT
Definition: Exposures.h:76
Sum MC predictions from different periods scaled according to data POT targets.
const double kSecondAnaPOT
Definition: Exposures.h:87
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:234
void ResetAngles(osc::OscCalcSterile *calc)
Definition: FCContour.C:270
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:78
Compare a data spectrum to MC expectation using only the event count.
FILE * outfile
Definition: dump_event.C:13
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17