make_nus_fc_surfs.C
Go to the documentation of this file.
1 // Ran at Caltech
2 
4 #include "CAFAna/Cuts/Cuts.h"
7 #include "CAFAna/FC/FCPoint.h"
8 #include "CAFAna/FC/FCSurface.h"
10 #include "CAFAna/Vars/FitVars.h"
12 #include "CAFAna/Core/Loaders.h"
18 #include "CAFAna/Analysis/Plots.h"
19 #include "CAFAna/Core/Progress.h"
24 #include "NuXAna/Systs/NusSysts.h"
25 #include "CAFAna/Core/Utilities.h"
26 #include "CAFAna/Vars/Vars.h"
27 #include "CAFAna/Analysis/Calcs.h"
28 
31 
32 using namespace ana;
33 
34 #include "TCanvas.h"
35 #include "TFile.h"
36 #include "TH2.h"
37 #include "TLegend.h"
38 #include "TMarker.h"
39 #include "TRandom3.h"
40 
41 // Reset calculator values after each fit to preset defaults
43 
44 void make_nus_fc_surfs(int nTrials = 100, int bin = 0, int Idx = 0)
45 {
46  //std::cout << "Bin: " << bin << ", nTrials: " << nTrials << std::endl;
47  // Load in the NuMI data spectrum
48  Spectrum* intime = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Ana01/fout_nus_ana01_box_opening_restricted.root","spec_nus_cale_numi").release();
49  //Spectrum* intime = LoadFromFile<Spectrum>("fout_nus_ana01_box_opening_restricted.root","spec_nus_cale_numi").release();
50 
51  TH1 *hSpectrum = intime->ToTH1(intime->POT());
52  Spectrum data_spec(hSpectrum, intime->POT(), intime->Livetime());
53  std::cout << "Data spectrum: " << data_spec.POT()
54  <<" POT and "
55  << data_spec.ToTH1(intime->POT())->GetSumOfWeights()
56  << " events" << std::endl;
57 
58 
59  // Load in the out-of-NuMI-time weighted spectrum (cosmics)
60  Spectrum* cosmic_scale = LoadFromFile<Spectrum>("$NOVA_ANA/steriles/Ana01/AnaResults.root","cosmic__CalE").release();
61  //Spectrum* cosmic_scale = LoadFromFile<Spectrum>("AnaResults.root","cosmic__CalE").release();
62  TH1D* hcosmic = (TH1D*)cosmic_scale->ToTH1(cosmic_scale->Livetime(), kLivetime);
63  Spectrum cosmics(hcosmic, intime->POT(), intime->Livetime());
64  double totcosmics = cosmics.ToTH1(intime->POT())->GetSumOfWeights();
65  std::cout << "Total cosmic count: " << totcosmics << std::endl;
66 
67  // Load in MC prediction
68  TFile input("$NOVA_ANA/steriles/Ana01/nus_ana01_prediction.root", "READ");
69  //TFile input("nus_ana01_prediction.root", "READ");
70 
71  std::unique_ptr<PredictionInterp> pred_p1p2e3b =
72  PredictionInterp::LoadFrom((TDirectory*)input.Get("nus_pred_p1p2e3b"));
73 
74  std::unique_ptr<PredictionInterp> pred_e3c =
75  PredictionInterp::LoadFrom((TDirectory*)input.Get("nus_pred_e3c"));
76 
77  const double pot_p1p2e3b =
81 
82  const double pot_e3ce3d = kSecondAnaEpoch3cPOT + kSecondAnaEpoch3dPOT;
83 
84  // Combine period1/period2/epoch3b + epoch3c, weighted appropriately
86  {&(*pred_p1p2e3b),pot_p1p2e3b},
87  {&(*pred_e3c),pot_e3ce3d}
88  });
89 
90 
92  calc4f->SetNFlavors(4);
93  calc4f->SetDm(4,0.5); // #Delta m^{2}_{41} = 0.5 eV^2
94  calc4f->SetAngle(2, 4, 0.); // #theta_{24} = 0 degrees
95  calc4f->SetAngle(3, 4, 0.); // #theta_{34} = 0 degrees
96 
97  // theta_34
98  double minX = 0.;
99  double maxX = 50.;
100  int NX = 50;
101  double widthX = (maxX-minX)/(NX-1);
102 
103  // theta_24
104  double minY = 0.;
105  double maxY = 45.;
106  double NY = 45.;
107  double widthY = (maxY-minY)/(NY-1);
108 
109  assert(bin < NX*NY && "Invalid bin number");
110  int binY = bin/NX;
111  int binX = bin - NX*binY;
112  double centerX = minX + binX*widthX;
113  double centerY = minY + binY*widthY;
114 
115  // Vector of parameters to fit: #theta_{23}, #theta_{34}, #theta_{24}
116  std::vector<const IFitVar*> fitvars;
117  fitvars.push_back(&kFitTheta34InDegreesSterile);
118  fitvars.push_back(&kFitTheta24InDegreesSterile);
119  fitvars.push_back(&kFitTheta23InDegreesSterile);
120 
121  FCCollection fccol;
122  //Zero sets the seed to be a unique value in space and time
123  TRandom3 rnd(0); // zero seeds randomly
124 
125  //Progress prog("Performing mock experiments");
126 
127  for(int iTrial = 0; iTrial < nTrials; ++iTrial){
129  ResetAngles(calc4f);
130 
131  std::cout << "On trial: " << iTrial << std::endl;
132 
133 
134  //For each trial, getting random location in true X, Y space.
135  // Limits are the edges of 2D plot you want contours drawn on.
136  // Throw a true X,Y value, randomly accross current bin
137  double trueX = -1;
138  // Edge bins centered on a physical edge
139  while (trueX < minX || trueX > maxX)
140  trueX = rnd.Uniform(centerX-widthX/2,centerX+widthX/2);
141  double trueY = -1;
142  while (trueY < minY || trueY > maxY)
143  trueY = rnd.Uniform(centerY-widthY/2,centerY+widthY/2);
144 
145  //const double trueX = rnd.Uniform(0.0, 50.);
146  //const double trueY = rnd.Uniform(0.0, 45.);
147 
148  kFitTheta34InDegreesSterile.SetValue(calc4f, trueX);
149  kFitTheta24InDegreesSterile.SetValue(calc4f, trueY);
150 
151 
152  // SystShifts shifts;
153  std::vector<const ISyst*> systlist = getAllNusSysts();
154 
155  // for(const auto& syst:systlist) {
156  // shifts.SetShift(syst, rnd.Gaus(0, 1));
157  // }
158 
159  const double pot = kSecondAnaPOT;
160 
161  Spectrum obs = pred.Predict(calc4f);
162  obs += cosmics;
163  Spectrum mock = obs.MockData(pot);
164 
165  // Compares (mock)data sample with MC sample
167  45.8,
168  3.2);
169  CountingExperiment nus_expt(&pred, mock, cosmics);
170  MultiExperiment expt({&th23Constraint, &nus_expt});
171 
172  double trueChi;
173  kFitTheta23InDegreesSterile.SetValue(calc4f, rnd.Gaus(45.8,3.2));
174 
175  MinuitFitter fittrue(&expt, {&kFitTheta23InDegreesSterile}, systlist);
176  trueChi = fittrue.Fit(calc4f, IFitter::kQuiet)->EvalMetricVal();
177 
178  MinuitFitter fit(&expt, fitvars, systlist);
179  double bestChi = fit.Fit(calc4f, IFitter::kQuiet)->EvalMetricVal();
180 
181  // The fit set the best fit values in the calculator
182  // now we are getting them back
183  const double bestX = kFitTheta34InDegreesSterile.GetValue(calc4f);
184  const double bestY = kFitTheta24InDegreesSterile.GetValue(calc4f);
185 
186  // Now we are setting true oscillation values, chi squared value
187  // when fit with true points, best fit oscillation values, and
188  // chi squared value for the best fit point
189  if (trueChi < bestChi) {
190  std::cout << "True Chi less than Best Chi" << std::endl;
191  bestChi = trueChi;
192  }
193 
194 
195  FCPoint pt(trueX, trueY, trueChi,bestX, bestY, bestChi,0); // Don't know seed or rng state
196  if (100*iTrial % nTrials == 0){
197  std::cout << "Trial: " << iTrial << ", TrueX: " << trueX << ", Y: " << trueY << ", trueChi: " << trueChi << ", bestX: " << bestX << ", Y: " << bestY << ", bestChi: " << bestChi << ", true-best: " << trueChi-bestChi << std::endl;
198  }
199 
200  fccol.AddPoint(pt);
201  }
202 
203  std::string binname = std::to_string(bin);
204  std::string number = std::to_string(Idx);
205 
206  fccol.SaveToFile("FCCol_"+binname+"_"+number+".root");
207 }
208 
209 
211 {
212  calc->SetAngle(2, 3, M_PI/4); // 45 degress, maximal mixing
213  calc->SetAngle(3, 4, 0.); // 0 degrees
214  calc->SetAngle(2, 4, 0.); // 0 degrees
215 
216  return;
217 }
TCut intime("tave>=217.0 && tave<=227.8")
double minY
Definition: plot_hist.C:9
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
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:148
A simple Gaussian constraint on an arbitrary IFitVar.
const double kSecondAnaEpoch3dPOT
Definition: Exposures.h:77
void AddPoint(const FCPoint &p)
Definition: FCCollection.h:17
const double kSecondAnaEpoch3bPOT
Definition: Exposures.h:75
Adapt the PMNS_Sterile calculator to standard interface.
const FitTheta23InDegreesSterile kFitTheta23InDegreesSterile
double maxY
Definition: plot_hist.C:10
const double kSecondAnaPeriod1POT
Definition: Exposures.h:73
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
osc::OscCalcDumb calc
#define M_PI
Definition: SbMath.h:34
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:300
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Represents the results of a single Feldman-Cousins pseudo-experiment.
Definition: FCPoint.h:8
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
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 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.
Definition: IFitter.cxx:69
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
#define pot
float bin[41]
Definition: plottest35.C:14
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:76
double POT() const
Definition: Spectrum.h:227
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
void SetAngle(int i, int j, double th)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir, const std::string &name)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
void make_nus_fc_surfs(int nTrials=100, int bin=0, int Idx=0)
assert(nhit_max >=nhit_nbins)
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
void ResetAngles(osc::OscCalcSterile *calc)
const double kSecondAnaEpoch3cPOT
Definition: Exposures.h:76
Sum MC predictions from different periods scaled according to data POT targets.
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void SaveToFile(const std::string &fname) const
mock
Definition: demo4.py:28
const double kSecondAnaPOT
Definition: Exposures.h:87
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
Compare a data spectrum to MC expectation using only the event count.
Collection of FCPoint. Serializable to/from a file.
Definition: FCCollection.h:12
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string