make_nus17_fc_surfs.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TFile.h"
3 #include "TH2.h"
4 #include "TLegend.h"
5 #include "TMarker.h"
6 #include "TRandom3.h"
7 #include "TDirectory.h"
8 #include "TH1.h"
9 #include "TLatex.h"
10 #include "TLine.h"
11 #include "TPad.h"
12 
13 #include "CAFAna/Cuts/Cuts.h"
14 #include "CAFAna/Cuts/SpillCuts.h"
15 #include "CAFAna/FC/FCCollection.h"
16 #include "CAFAna/FC/FCPoint.h"
17 #include "CAFAna/FC/FCSurface.h"
18 #include "CAFAna/Vars/FitVars.h"
20 #include "CAFAna/Vars/Vars.h"
28 #include "CAFAna/Core/Utilities.h"
29 #include "CAFAna/Core/Progress.h"
30 #include "CAFAna/Core/Loaders.h"
32 #include "CAFAna/Core/IFitVar.h"
33 #include "CAFAna/Analysis/Calcs.h"
36 #include "CAFAna/Analysis/Style.h"
38 #include "CAFAna/Analysis/Plots.h"
40 #include "OscLib/OscCalcSterile.h"
41 
43 
44 using namespace ana;
45 
46 // Global variables
47 bool kLowerOctant = true;
48 
49 // Reset calculator values after each fit to preset defaults
51 
52 void make_nus17_fc_surfs(int nTrials = 100, int bin = 0, int Idx = 0, bool lowerOctant = true) {
53  kLowerOctant = lowerOctant;
54  TH1::AddDirectory(0);
55 
56  // Set up path to file with filled CAFAna objects and for output of analysis
57  std::string folder = "/nova/ana/steriles/Nus17/";
58  std::string filenm = "predInterpSysts_nus17_v1";
59 
60  std::string loadLocation = "/nova/ana/steriles/Nus17/predInterpSysts_nus17_v1.root";
61  std::string saveLocation = "plot_" + filenm + "_shape_2D_lower.root";
62  std::string textLocation = "count_" + filenm + "_shape_2D_lower.txt";
63 
64  // Load in the NuMI data spectrum
65  Spectrum* intime = LoadFromFile<Spectrum>("/nova/ana/steriles/Nus17/fout_nus17_box_opening_restricted.root","spec_nus_vise_numi").release();
66 
67  TH1 *hSpectrum = intime->ToTH1(intime->POT());
68  Spectrum data_spec(hSpectrum, intime->POT(), intime->Livetime());
69  std::cout << "Data spectrum: " << data_spec.POT()
70  <<" POT and "
71  << data_spec.ToTH1(intime->POT())->GetSumOfWeights()
72  << " events" << std::endl;
73  delete hSpectrum;
74 
75  TFile *cos = new TFile("/nova/ana/steriles/Nus17/nus17_cosmic_background.root","read");
76  TH1 *hcosmic = (TH1*)cos->Get("cosall");
78  std::cout << "Cosmic Data spectrum: " << cosmic.POT()
79  <<" POT and "
80  << cosmic.ToTH1(cosmic.POT())->GetSumOfWeights()
81  << " events"
82  << ", Livetime: "
83  << cosmic.Livetime() << std::endl;
84  delete hcosmic;
85  cos->Close();
86 
87  // Load filled objects from file
88  TFile* rootL = new TFile(loadLocation.c_str(), "read");
89  TDirectory* tmpL = gDirectory;
90  TDirectory* loadDir = gDirectory;
91 
92  loadDir->cd((loadLocation + ":/nus17_predI").c_str());
93  std::unique_ptr<PredictionInterp> pred = PredictionInterp::LoadFrom(gDirectory);
94 
95  tmpL->cd();
96  rootL->Close(); // Don't forget to close the file!
97 
98  // OscCalc for Sterile 3+1 Model
100  calc4f->SetNFlavors(4);
101  ResetAngles(calc4f);
102 
103  // theta_34
104  double minX = 0.;
105  double maxX = 50.;
106  int NX = 50;
107  double widthX = (maxX-minX)/(NX-1);
108 
109  // theta_24
110  double minY = 0.;
111  double maxY = 45.;
112  double NY = 45.;
113  double widthY = (maxY-minY)/(NY-1);
114 
115  assert(bin < NX*NY && "Invalid bin number");
116  int binY = bin/NX;
117  int binX = bin - NX*binY;
118  double centerX = minX + binX*widthX;
119  double centerY = minY + binY*widthY;
120 
121  // Vector of parameters to fit: #theta_{23}, #theta_{34}, #theta_{24}, #deltam_{32}
122  std::vector<const IFitVar*> fitvars;
123  fitvars.push_back(&kFitTheta34InDegreesSterile);
124  fitvars.push_back(&kFitTheta24InDegreesSterile);
125  fitvars.push_back(&kFitSinSqTheta23Sterile); // to be constrained to PDG
126  fitvars.push_back(&kFitDmSq32Sterile); // to be constrained to PDG
127  fitvars.push_back(&kFitDelta24InPiUnitsSterile); // profile delta_24
128 
129  FCCollection fccol;
130  TRandom3 rnd(0); // zero seeds randomly
131 
132  for(int iTrial = 0; iTrial < nTrials; ++iTrial){
134  ResetAngles(calc4f);
135 
136  std::cout << "On trial: " << iTrial << std::endl;
137 
138  // Throw a true X,Y value, randomly accross current bin
139  double trueX = -1;
140  // Edge bins centered on a physical edge
141  while (trueX < minX || trueX > maxX)
142  trueX = rnd.Uniform(centerX-widthX/2,centerX+widthX/2);
143  double trueY = -1;
144  while (trueY < minY || trueY > maxY)
145  trueY = rnd.Uniform(centerY-widthY/2,centerY+widthY/2);
146 
147  kFitTheta34InDegreesSterile.SetValue(calc4f, trueX);
148  kFitTheta24InDegreesSterile.SetValue(calc4f, trueY);
149  kFitDelta24InPiUnitsSterile.SetValue(calc4f, rnd.Uniform(0, 2));
150 
151  // SystShifts shifts;
152  std::vector<const ISyst*> systlist = getAllNusFDSysts();
153 
154  // Update POT to 2017 Analysis POT
155  const double pot = kAna2017POT;
156 
157  Spectrum obs = pred.get()->Predict(calc4f);
158  obs += cosmic;
159  Spectrum mock = obs.MockData(pot);
160  TH1 *hMockSpectrum = mock.ToTH1(pot);
161  Spectrum mock_spec(hMockSpectrum, pot, intime->Livetime());
162  delete hMockSpectrum;
163 
164  // Compares (mock)data sample with MC sample
165  Double_t meanSinSqTheta23Sterile = 0;
166  Double_t sigmaSinSqTheta23Sterile = 0;
167 
168  if (kLowerOctant) {
169  meanSinSqTheta23Sterile = 0.404;
170  sigmaSinSqTheta23Sterile = 0.030;
171  } else {
172  meanSinSqTheta23Sterile = 0.603;
173  sigmaSinSqTheta23Sterile = 0.030;
174  }
175 
176  GaussianConstraint s2th23Constraint(&kFitSinSqTheta23Sterile, meanSinSqTheta23Sterile, sigmaSinSqTheta23Sterile);
177  GaussianConstraint dmsq32Constraint(&kFitDmSq32Sterile, 2.67e-3, 0.11e-3);
178 
179  SingleSampleExperiment nus_expt(pred.get(), mock_spec, cosmic);
180  MultiExperiment expt({&s2th23Constraint, &dmsq32Constraint, &nus_expt});
181 
182  if (kLowerOctant) {
183  kFitSinSqTheta23Sterile.SetValue(calc4f, rnd.Gaus(0.404, 0.030));
184  } else {
185  kFitSinSqTheta23Sterile.SetValue(calc4f, rnd.Gaus(0.603, 0.030));
186  }
187  kFitDmSq32Sterile.SetValue(calc4f, rnd.Gaus(2.67e-3, 0.11e-3));
188 
190  double trueChi = fittrue.Fit(calc4f, IFitter::kQuiet)->EvalMetricVal();
191 
192  if (kLowerOctant) {
193  kFitSinSqTheta23Sterile.SetValue(calc4f, rnd.Gaus(0.404, 0.030));
194  } else {
195  kFitSinSqTheta23Sterile.SetValue(calc4f, rnd.Gaus(0.603, 0.030));
196  }
197  kFitDmSq32Sterile.SetValue(calc4f, rnd.Gaus(2.67e-3, 0.11e-3));
198 
199  MinuitFitter fit(&expt, fitvars, systlist);
200  double bestChi = fit.Fit(calc4f, IFitter::kQuiet)->EvalMetricVal();
201 
202  // The fit set the best fit values in the calculator
203  // now we are getting them back
204  const double bestX = kFitTheta34InDegreesSterile.GetValue(calc4f);
205  const double bestY = kFitTheta24InDegreesSterile.GetValue(calc4f);
206 
207  // Now we are setting true oscillation values, chi squared value
208  // when fit with true points, best fit oscillation values, and
209  // chi squared value for the best fit point
210  if (trueChi < bestChi) {
211  std::cout << "True Chi less than Best Chi" << std::endl;
212  bestChi = trueChi;
213  }
214 
215  FCPoint pt(trueX, trueY, trueChi,bestX, bestY, bestChi,0); // Don't know seed or rng state
216  if (100*iTrial % nTrials == 0){
217  std::cout << "Trial: " << iTrial << ", TrueX: " << trueX << ", Y: " << trueY << ", trueChi: " << trueChi << ", bestX: " << bestX << ", Y: " << bestY << ", bestChi: " << bestChi << ", true-best: " << trueChi-bestChi << std::endl;
218  }
219 
220  fccol.AddPoint(pt);
221  }
222 
223  std::string binname = std::to_string(bin);
224  std::string number = std::to_string(Idx);
225  std::string octant = "_lowerOctant";
226  if (!kLowerOctant) {
227  octant = "_upperOctant";
228  }
229 
230  fccol.SaveToFile("FC_nus17_Col_"+binname+"_"+number+octant+".root");
231 }
232 
233 
235 {
237  calc->SetDm(4, 0.5);
238  calc->SetAngle(3, 4, 0); // 0, icecube
239  calc->SetAngle(2, 4, 0); // 0, sk
240  calc->SetDm(3, 7.53e-5 + 2.67e-3);
241  if (kLowerOctant) {
242  calc->SetAngle(2, 3, (asin(sqrt(0.404))));
243  calc->SetDelta(1, 3, 1.48*M_PI);
244  } else {
245  calc->SetAngle(2, 3, (asin(sqrt(0.603))));
246  calc->SetDelta(1, 3, 0.74*M_PI);
247  }
248  return;
249 }
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
void SetNFlavors(int nflavors)
void ResetAngles(osc::OscCalcSterile *calc)
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.
void make_nus17_fc_surfs(int nTrials=100, int bin=0, int Idx=0, bool lowerOctant=true)
T sqrt(T number)
Definition: d0nt_math.hpp:156
void SetDelta(int i, int j, double delta)
void AddPoint(const FCPoint &p)
Definition: FCCollection.h:17
Adapt the PMNS_Sterile calculator to standard interface.
const FitTheta23InDegreesSterile kFitTheta23InDegreesSterile
const FitDmSq32Sterile kFitDmSq32Sterile
const double kAna2017POT
Definition: Exposures.h:174
double maxY
Definition: plot_hist.C:10
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
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
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
const double kAna2017Livetime
Definition: Exposures.h:200
Sum up livetimes from individual cosmic triggers.
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
bool kLowerOctant
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 SetDm(int i, double dm)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
const FitSinSqTheta23Sterile kFitSinSqTheta23Sterile
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void SaveToFile(const std::string &fname) const
mock
Definition: demo4.py:28
Float_t e
Definition: plot.C:35
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
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