FitParamEffectsAna.C
Go to the documentation of this file.
3 #include "CAFAna/Fit/Fit.h"
12 
13 #include "OscLib/OscCalcSterile.h"
14 
15 #include "TArrow.h"
16 #include "TCanvas.h"
17 #include "TDirectory.h"
18 #include "TFile.h"
19 #include "TH1.h"
20 #include "TKey.h"
21 #include "TLatex.h"
22 #include "TLegend.h"
23 #include "TLine.h"
24 #include "TPad.h"
25 
26 using namespace ana;
27 
28 struct AngleValues{
29  int nbins34;
30  double min34;
31  double max34;
32 
33  int nbins24;
34  double min24;
35  double max24;
36 };
37 
38 double CalcSystEffect(double statonly, double withsyst);
39 
40 double FindAngleFromChi2(double chi2, TH1* h);
41 
42 void LoadMaps(TDirectory* dir,
43  std::map<std::string, PredictionCombinePeriods*>* preds);
44 
46  const IFitVar* var, std::vector<const IFitVar*> profVars,
47  TDirectory* rootOut, int nbins, double min, double max,
48  std::string name, Color_t color,
49  std::vector<const ISyst*> systs);
50 
52  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24, const IFitVar* vard24,
53  TDirectory* rootOut, AngleValues a,
54  std::string years, Color_t color,
55  std::string samplelabel,
56  std::map<std::string, TH1*>* h24s,
57  std::map<std::string, TH1*>* h34s);
58 
60  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24, const IFitVar* vard24,
61  const IFitVar* extra,
62  TDirectory* rootOut, AngleValues a,
63  std::string years, Color_t color,
64  std::string samplelabel,
65  std::map<std::string, TH1*>* h24s,
66  std::map<std::string, TH1*>* h34s);
67 
68 void PrintAngles(FILE* text, std::string label,
69  double th2468, double th2490,
70  double th3468, double th3490);
71 
72 void PrintSystEffect(FILE* text,
73  double th24stat68, double th2468,
74  double th24stat90, double th2490,
75  double th34stat68, double th3468,
76  double th34stat90, double th3490);
77 
78 // After running a fit, reset the values in a calculator to preset defaults
80 
82 {
83  TH1::AddDirectory(0);
84 
85  std::map<std::string, PredictionCombinePeriods*> preds;
86 
87  // Set up path to file with filled CAFAna objects and for output of analysis
88  std::string folder = "";
89  std::string filenm = "FitParamEffects";
90 
91  std::string loadLocation = "pred_fitsysteffects_ana01_fix.root";
92  std::string saveLocation = folder + filenm + "Ana.root";
93  std::string textLocation = folder + filenm + "Ana.txt";
94 
95  // Load filled objects from file
96  TFile* rootL = new TFile(loadLocation.c_str(), "READ");
97  TDirectory* tmpL = gDirectory;
98  TDirectory* loadDir = gDirectory;
99  loadDir->cd((loadLocation + ":/pred").c_str());
101 
102  loadDir->cd((loadLocation + ":/sCosmic").c_str());
103  Spectrum sCosmic = *Spectrum::LoadFrom(gDirectory);
104  loadDir->cd((loadLocation + ":/sData").c_str());
105  Spectrum sData = *Spectrum::LoadFrom(gDirectory);
106 
107  rootL->Close(); // Don't forget to close the file!
108 
109  TH1* hCosmic = sCosmic.ToTH1(sCosmic.Livetime(), kLivetime);
110  Spectrum sCos(hCosmic, sData.POT(), sData.Livetime());
111 
112  int nData = sData.Integral(sData.POT());
113  std::cout << "Data spectrum: " << sData.POT()
114  << " and " << nData << " events" << std::endl;
115  double nCos = sCos.Integral(sCos.POT());
116  std::cout << "Cos spectrum: " << sCos.POT()
117  << " and " << nCos << " events" << std::endl;
118 
119  // Set up oscillation calculator that uses default 3 flavor parameters
121 
123  calc4f->SetNFlavors(4);
124  ResetAngles(calc4f);
125 
126  // Set up values for slice and surface plots
127  AngleValues avals;
128  avals.nbins34 = 20;
129  avals.min34 = 0.;
130  avals.max34 = 50.;
131  avals.nbins24 = 18;
132  avals.min24 = 0.;
133  avals.max24 = 45.;
134 
135  const Color_t kFitColor = kRed;
136 
137  std::map<std::string, std::pair<const IFitVar*, std::pair<double, double> > > parammap;
138  parammap["t13"] = {&kFitTheta13InDegreesSterile, {8.475, 0.0353}};
139  parammap["d13"] = {&kFitDelta13InPiUnitsSterile, {1.39, 0.325}};
140  parammap["m32"] = {&kFitDmSq32Sterile, {0.00244, 0.00006}};
141 
142  std::map<std::string, TH1*> h24s;
143  std::map<std::string, TH1*> h34s;
144 
145  // Open files for saving analysis output
146  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
147  FILE* textF;
148  textF = fopen(textLocation.c_str(), "w");
149 
150  if(true) {
151  std::cout << "NoParams" << std::endl;
152  ResetAngles(calc4f);
153 
154  GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile, 45.8, 3.2);
155  CountingExperiment expt(&pred, sData, sCos);
156  MultiExperiment multi({&th23Constraint, &expt});
157 
158  // Plot 1D fits for 1 data yr experiment
159  Plot1DSlices(
160  &multi, calc4f,
165  rootF, avals, "1", kFitColor,
166  "NoParams", &h24s, &h34s
167  );
168  }
169 
170  double th24stat68 = FindAngleFromChi2(1., h24s["NoParams"]);
171  double th24stat90 = FindAngleFromChi2(2.706, h24s["NoParams"]);
172  double th34stat68 = FindAngleFromChi2(1., h34s["NoParams"]);
173  double th34stat90 = FindAngleFromChi2(2.706, h34s["NoParams"]);
174  PrintAngles(textF, "NoParams", th24stat68, th24stat90, th34stat68, th34stat90);
175 
176  for(const auto& parampair : parammap) {
177  std::string label = parampair.first;
178  std::cout << "Label: " << label << std::endl;
179  ResetAngles(calc4f);
180 
181  GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile, 45.8, 3.2);
182  GaussianConstraint extraConstraint(parampair.second.first,
183  parampair.second.second.first,
184  parampair.second.second.second);
185  CountingExperiment expt(&pred, sData, sCos);
186  MultiExperiment multi({&extraConstraint, &th23Constraint, &expt});
187 
188  // Plot 1D fits for 1 data yr experiment
189  Plot1DSlices(
190  &multi, calc4f,
195  parampair.second.first,
196  rootF, avals, "1", kFitColor,
197  label, &h24s, &h34s
198  );
199 
200  double th2468 = FindAngleFromChi2(1., h24s[label]);
201  double th2490 = FindAngleFromChi2(2.706, h24s[label]);
202  double th3468 = FindAngleFromChi2(1., h34s[label]);
203  double th3490 = FindAngleFromChi2(2.706, h34s[label]);
204  PrintAngles(textF, label, th2468, th2490, th3468, th3490);
205  PrintSystEffect(textF,
206  th24stat68, th2468, th24stat90, th2490,
207  th34stat68, th3468, th34stat90, th3490);
208  }
209 
210  fclose(textF); // Text file is done now, close it
211  rootF->Close(); // Close the output root file
212 }
213 
214 //------------------------------------------------------------------------------
215 double CalcSystEffect(double statonly, double withsyst)
216 {
217  return 100.*(withsyst - statonly)/statonly;
218 }
219 
220 //------------------------------------------------------------------------------
221 double FindAngleFromChi2(double chi2, TH1* h)
222 {
223  for(int i = 1, n = h->GetNbinsX(); i < n; ++i) {
224  if(h->GetBinContent(i, 1) < chi2 &&
225  h->GetBinContent(i+1, 1) >= chi2) {
226  return h->GetXaxis()->GetBinCenter(i+1);
227  }
228  }
229 
230  return h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
231 }
232 
233 //------------------------------------------------------------------------------
234 /// Load CAFAna objects from a TDirectory
235 /// The maps to contain these objects must already exist
236 void LoadMaps(TDirectory* dir,
237  std::map<std::string, PredictionCombinePeriods*>* preds)
238 {
239  TDirectory* tmp = gDirectory;
240  TDirectory* keyDir = gDirectory;
241 
242  // Loop over all keys in this directory
243  dir->cd();
244  TIter nextkey(dir->GetListOfKeys());
245  TKey* key;
246  TKey* oldkey = 0;
247 
248  std::string sep = "__"; // This is the token separator
249  std::string cafanaType = "";
250  std::string sampleLabel = "";
251 
252  while((key = (TKey*)nextkey())) {
253  // Keep only the highest cycle number for each key
254  if(oldkey && !strcmp(oldkey->GetName(), key->GetName())) { continue; }
255  if(key->ReadObj()->IsFolder()) {
256  keyDir = (TDirectory*)key->ReadObj();
257  }
258  else {
259  continue;
260  }
261 
262  std::string name(key->GetName());
263  std::vector<std::string> tokens;
264  std::size_t pos = 0, prev = 0;
265  while(pos != std::string::npos) { // Loop over the TObject name
266  pos = name.find(sep, prev); // Find the next separator
267 
268  if(pos != std::string::npos) {
269  tokens.push_back(name.substr(prev, pos-prev));
270  prev = pos+2; // Move past the current separator
271  }
272  else {
273  tokens.push_back(name.substr(prev, name.length()-prev));
274  }
275  }
276 
277  if(tokens.size() > 1) {
278  cafanaType = tokens[0];
279  sampleLabel = tokens[1];
280  if(cafanaType.compare("pred") == 0) {
281  (*preds)[sampleLabel] = PredictionCombinePeriods::LoadFrom(keyDir).release();
282  }
283  }
284  else { continue; }
285  }
286 
287  tmp->cd();
288  return;
289 }
290 
291 //------------------------------------------------------------------------------
293  const IFitVar* var, std::vector<const IFitVar*> profVars,
294  TDirectory* rootOut, int nbins, double min, double max,
295  std::string name, Color_t color)
296 {
297  // Create slice histogram with/without profiling based on input profVars vector size
298 
299  TH1* h = Profile(expt, calc, var, nbins, min, max, -1, profVars);
300 
301  ResetAngles(calc); // Don't forget to reset calculator...
302 
303  // Formatting:
304  h->SetLineColor(color); // Set the line color
305  h->SetMaximum(5); // Want the plot to go from 0 to 5
306  h->SetMinimum(0);
307  h->SetName(name.c_str()); // Set the object name for saving
308  rootOut->WriteTObject(h); // Save the raw histogram
309 
310  return h;
311 }
312 
313 //------------------------------------------------------------------------------
315  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24, const IFitVar* vard24,
316  TDirectory* rootOut, AngleValues a,
317  std::string years, Color_t color,
318  std::string samplelabel,
319  std::map<std::string, TH1*>* h24s,
320  std::map<std::string, TH1*>* h34s)
321 {
322  // General form of the object name for saving
323  std::string name = "h" + years + "yrTh";
324  std::string fullname; // fullname adds 2 characters, the theta subscript indices
325 
326  // Plot profiled slices in same order as above
327  // When slice is of one angle, the other two are profiled
328  // I.e., when the slice is theta 23, theta 34 and 24 are profiled
329  fullname = name + "34Prof" + samplelabel;
330  (*h34s)[samplelabel] = Plot1DSlice(
331  expt, calc, var34, {var23, var24, vard24}, rootOut,
332  a.nbins34, a.min34, a.max34, fullname, color
333  );
334  fullname = name + "24Prof" + samplelabel;
335  (*h24s)[samplelabel] = Plot1DSlice(
336  expt, calc, var24, {var23, var34, vard24}, rootOut,
337  a.nbins24, a.min24, a.max24, fullname, color
338  );
339 
340  return;
341 }
342 
343 //------------------------------------------------------------------------------
345  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24, const IFitVar* vard24,
346  const IFitVar* extra,
347  TDirectory* rootOut, AngleValues a,
348  std::string years, Color_t color,
349  std::string samplelabel,
350  std::map<std::string, TH1*>* h24s,
351  std::map<std::string, TH1*>* h34s)
352 {
353  // General form of the object name for saving
354  std::string name = "h" + years + "yrTh";
355  std::string fullname; // fullname adds 2 characters, the theta subscript indices
356 
357  // Plot profiled slices in same order as above
358  // When slice is of one angle, the other two are profiled
359  // I.e., when the slice is theta 23, theta 34 and 24 are profiled
360  fullname = name + "34Prof" + samplelabel;
361  (*h34s)[samplelabel] = Plot1DSlice(
362  expt, calc, var34, {var23, var24, vard24, extra}, rootOut,
363  a.nbins34, a.min34, a.max34, fullname, color
364  );
365  fullname = name + "24Prof" + samplelabel;
366  (*h24s)[samplelabel] = Plot1DSlice(
367  expt, calc, var24, {var23, var34, vard24, extra}, rootOut,
368  a.nbins24, a.min24, a.max24, fullname, color
369  );
370 
371  return;
372 }
373 
374 //------------------------------------------------------------------------------
376  double th2468, double th2490,
377  double th3468, double th3490)
378 {
379  fprintf(text, "Angle limits for %s:\n", label.c_str());
380  fprintf(text, "Th24: 68%%: %4.2f, 90%%: %4.2f\n", th2468, th2490);
381  fprintf(text, "Th34: 68%%: %4.2f, 90%%: %4.2f\n", th3468, th3490);
382  fprintf(text, "\n");
383  return;
384 }
385 
386 //------------------------------------------------------------------------------
388  double th24stat68, double th2468,
389  double th24stat90, double th2490,
390  double th34stat68, double th3468,
391  double th34stat90, double th3490)
392 {
393  fprintf(text, "Diff24: 68%% %4.2f, 90%%: %4.2f\n",
394  CalcSystEffect(th24stat68, th2468),
395  CalcSystEffect(th24stat90, th2490));
396  fprintf(text, "Diff34: 68%% %4.2f, 90%%: %4.2f\n",
397  CalcSystEffect(th34stat68, th3468),
398  CalcSystEffect(th34stat90, th3490));
399  fprintf(text, "\n");
400  return;
401 }
402 
403 //------------------------------------------------------------------------------
405 {
407  calc->SetAngle(2, 3, M_PI/4); // 45 degrees (maximal mixing)
408  calc->SetAngle(3, 4, 0.);
409  calc->SetAngle(2, 4, 0.);
410 
411  return;
412 }
413 
const XML_Char * name
Definition: expat.h:151
enum BeamMode kRed
const FitTheta13InDegreesSterile kFitTheta13InDegreesSterile
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetNFlavors(int nflavors)
cons_index_list< index_multi, nil_index_list > multi
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
static std::unique_ptr< PredictionCombinePeriods > LoadFrom(TDirectory *dir, const std::string &name)
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.
double CalcSystEffect(double statonly, double withsyst)
Adapt the PMNS_Sterile calculator to standard interface.
const FitTheta23InDegreesSterile kFitTheta23InDegreesSterile
const FitDmSq32Sterile kFitDmSq32Sterile
TGraph * Profile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double input_minchi, const std::vector< const IFitVar * > &profVars, const std::vector< const ISyst * > &profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &profSystsMap, MinuitFitter::FitOpts opts)
scan in one variable, profiling over all others
Definition: Fit.cxx:48
Float_t tmp
Definition: plot.C:36
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:249
osc::OscCalcDumb calc
double chi2()
#define M_PI
Definition: SbMath.h:34
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void ResetAngles(osc::OscCalcSterile *calc)
void PrintSystEffect(FILE *text, double th24stat68, double th2468, double th24stat90, double th2490, double th34stat68, double th3468, double th34stat90, double th3490)
double FindAngleFromChi2(double chi2, TH1 *h)
const char * label
fclose(fg1)
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
const int nbins
Definition: cellShifts.C:15
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
expt
Definition: demo5.py:34
void prev()
Definition: show_event.C:91
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
const double a
void FitParamEffectsAna()
const FitDelta13InPiUnitsSterile kFitDelta13InPiUnitsSterile
void ResetSterileCalcToDefault(osc::OscCalcSterile *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:76
double POT() const
Definition: Spectrum.h:227
const char sep
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Combine multiple component experiments.
void SetAngle(int i, int j, double th)
TDirectory * dir
Definition: macro.C:5
TH1 * Plot1DSlice(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var, std::vector< const IFitVar * > profVars, TDirectory *rootOut, int nbins, double min, double max, std::string name, Color_t color, std::vector< const ISyst * > systs)
Base class defining interface for experiments.
Definition: IExperiment.h:14
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
Interface definition for fittable variables.
Definition: IFitVar.h:16
void LoadMaps(TDirectory *dir, std::map< std::string, IDecomp * > *nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > *shifted)
Definition: PPFXHelper.h:273
Sum MC predictions from different periods scaled according to data POT targets.
void Plot1DSlices(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, const IFitVar *vard24, TDirectory *rootOut, AngleValues a, std::string years, Color_t color, std::string samplelabel, std::map< std::string, TH1 * > *h24s, std::map< std::string, TH1 * > *h34s)
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
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.
void PrintAngles(FILE *text, std::string label, double th2468, double th2490, double th3468, double th3490)
enum BeamMode string