FitSystEffectsAna.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  std::vector<const ISyst*> systs = {});
59 
60 void PrintAngles(FILE* text, std::string label,
61  double th2468, double th2490,
62  double th3468, double th3490);
63 
64 void PrintSystEffect(FILE* text,
65  double th24stat68, double th2468,
66  double th24stat90, double th2490,
67  double th34stat68, double th3468,
68  double th34stat90, double th3490);
69 
70 // After running a fit, reset the values in a calculator to preset defaults
72 
74 {
75  TH1::AddDirectory(0);
76 
77  std::map<std::string, PredictionCombinePeriods*> preds;
78 
79  // Set up path to file with filled CAFAna objects and for output of analysis
80  std::string folder = "";
81  std::string filenm = "pred_fitsysteffects_ana01_fix";
82 
83  std::string loadLocation = folder + filenm + ".root";
84  std::string saveLocation = folder + filenm + "Ana.root";
85  std::string textLocation = folder + filenm + "Ana.txt";
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  loadDir->cd((loadLocation + ":/pred").c_str());
93 
94  loadDir->cd((loadLocation + ":/sCosmic").c_str());
95  Spectrum sCosmic = *Spectrum::LoadFrom(gDirectory);
96  loadDir->cd((loadLocation + ":/sData").c_str());
97  Spectrum sData = *Spectrum::LoadFrom(gDirectory);
98 
99  LoadMaps(rootL, &preds);
100 
101  rootL->Close(); // Don't forget to close the file!
102 
103  TH1* hCosmic = sCosmic.ToTH1(sCosmic.Livetime(), kLivetime);
104  Spectrum sCos(hCosmic, sData.POT(), sData.Livetime());
105 
106  int nData = sData.Integral(sData.POT());
107  std::cout << "Data spectrum: " << sData.POT()
108  << " and " << nData << " events" << std::endl;
109  double nCos = sCos.Integral(sCos.POT());
110  std::cout << "Cos spectrum: " << sCos.POT()
111  << " and " << nCos << " events" << std::endl;
112 
113  // Set up oscillation calculator that uses default 3 flavor parameters
115 
117  calc4f->SetNFlavors(4);
118  calc4f->SetDm(4, 0.5); //Dm41 = 0.5 eV^2
119  calc4f->SetAngle(2, 4, 0.);
120  calc4f->SetAngle(3, 4, 0.);
121 
122  // Set up values for slice and surface plots
123  AngleValues avals;
124  avals.nbins34 = 500;
125  avals.min34 = 0.;
126  avals.max34 = 50.;
127  avals.nbins24 = 450;
128  avals.min24 = 0.;
129  avals.max24 = 45.;
130 
131  const Color_t kFitColor = kRed;
132 
133  std::map<std::string, std::vector<const ISyst*> > systmap;
134  systmap["OscParam"] = std::vector<const ISyst*>();
135  systmap["Beam"] = std::vector<const ISyst*>();
136  systmap["Birks"] = std::vector<const ISyst*>();
137  systmap["Calib"] = std::vector<const ISyst*>();
138  systmap["DataMC"] = std::vector<const ISyst*>();
139  systmap["GENIE"] = std::vector<const ISyst*>();
140  systmap["NDCont"] = std::vector<const ISyst*>();
141  systmap["NDRock"] = std::vector<const ISyst*>();
142  systmap["Norm"] = std::vector<const ISyst*>();
143  systmap["Stats"] = std::vector<const ISyst*>();
144 
145  systmap["OscParam"] .push_back(&kNusOscParamSysts);
146  systmap["Beam"] .push_back(&kNusBeamSysts);
147  systmap["Birks"] .push_back(&kNusBirksSyst);
148  systmap["Calib"] .push_back(&kNusCalibFlatSyst);
149  systmap["Calib"] .push_back(&kNusCalibRelSyst);
150  systmap["Calib"] .push_back(&kNusCalibSlopeXSyst);
151  systmap["Calib"] .push_back(&kNusCalibSlopeYSyst);
152  systmap["DataMC"].push_back(&kNusNueCCSyst);
153  systmap["DataMC"].push_back(&kNusNumuCCSyst);
154  systmap["GENIE"] .push_back(&kNusGENIESmallSysts);
155  systmap["NDCont"].push_back(&kNusNDContSyst);
156  systmap["NDRock"].push_back(&kNusNDRockSyst);
157  systmap["Norm"] .push_back(&kNusNormSyst);
158  systmap["Stats"] .push_back(&kNusMCStatsSyst);
159 
160  systmap["All"] = getAllNusSysts();
161 
162  std::map<std::string, TH1*> h24s;
163  std::map<std::string, TH1*> h34s;
164 
165  // Open files for saving analysis output
166  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
167  FILE* textF;
168  textF = fopen(textLocation.c_str(), "w");
169 
170  if(true) {
171  std::cout << "NoSysts" << std::endl;
172  ResetAngles(calc4f);
173 
174  GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile, 45.8, 3.2);
175  CountingExperiment expt(&pred, sData, sCos);
176  MultiExperiment multi({&th23Constraint, &expt});
177 
178  // Plot 1D fits for 1 data yr experiment
179  Plot1DSlices(
180  &multi, calc4f,
185  rootF, avals, "1", kFitColor,
186  "NoSysts", &h24s, &h34s
187  );
188  }
189 
190  double th24stat68 = FindAngleFromChi2(1., h24s["NoSysts"]);
191  double th24stat90 = FindAngleFromChi2(2.706, h24s["NoSysts"]);
192  double th34stat68 = FindAngleFromChi2(1., h34s["NoSysts"]);
193  double th34stat90 = FindAngleFromChi2(2.706, h34s["NoSysts"]);
194  PrintAngles(textF, "NoSysts", th24stat68, th24stat90, th34stat68, th34stat90);
195 
196  for(const auto& predpair : preds) {
197  std::string label = predpair.first;
198  std::cout << "Label: " << label << std::endl;
199  ResetAngles(calc4f);
200 
201  GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile, 45.8, 3.2);
202  CountingExperiment expt(predpair.second, sData, sCos);
203  MultiExperiment multi({&th23Constraint, &expt});
204 
205  // Plot 1D fits for 1 data yr experiment
206  Plot1DSlices(
207  &multi, calc4f,
212  rootF, avals, "1", kFitColor,
213  label, &h24s, &h34s, systmap[label]
214  );
215 
216  double th2468 = FindAngleFromChi2(1., h24s[label]);
217  double th2490 = FindAngleFromChi2(2.706, h24s[label]);
218  double th3468 = FindAngleFromChi2(1., h34s[label]);
219  double th3490 = FindAngleFromChi2(2.706, h34s[label]);
220  PrintAngles(textF, label, th2468, th2490, th3468, th3490);
221  PrintSystEffect(textF,
222  th24stat68, th2468, th24stat90, th2490,
223  th34stat68, th3468, th34stat90, th3490);
224  }
225 
226  fclose(textF); // Text file is done now, close it
227  rootF->Close(); // Close the output root file
228 }
229 
230 //------------------------------------------------------------------------------
231 double CalcSystEffect(double statonly, double withsyst)
232 {
233  return 100.*(withsyst - statonly)/statonly;
234 }
235 
236 //------------------------------------------------------------------------------
237 double FindAngleFromChi2(double chi2, TH1* h)
238 {
239  for(int i = 1, n = h->GetNbinsX(); i < n; ++i) {
240  if(h->GetBinContent(i, 1) < chi2 &&
241  h->GetBinContent(i+1, 1) >= chi2) {
242  return h->GetXaxis()->GetBinCenter(i+1);
243  }
244  }
245 
246  return h->GetXaxis()->GetBinUpEdge(h->GetNbinsX());
247 }
248 
249 //------------------------------------------------------------------------------
250 /// Load CAFAna objects from a TDirectory
251 /// The maps to contain these objects must already exist
252 void LoadMaps(TDirectory* dir,
253  std::map<std::string, PredictionCombinePeriods*>* preds)
254 {
255  TDirectory* tmp = gDirectory;
256  TDirectory* keyDir = gDirectory;
257 
258  // Loop over all keys in this directory
259  dir->cd();
260  TIter nextkey(dir->GetListOfKeys());
261  TKey* key;
262  TKey* oldkey = 0;
263 
264  std::string sep = "__"; // This is the token separator
265  std::string cafanaType = "";
266  std::string sampleLabel = "";
267 
268  while((key = (TKey*)nextkey())) {
269  // Keep only the highest cycle number for each key
270  if(oldkey && !strcmp(oldkey->GetName(), key->GetName())) { continue; }
271  if(key->ReadObj()->IsFolder()) {
272  keyDir = (TDirectory*)key->ReadObj();
273  }
274  else {
275  continue;
276  }
277 
278  std::string name(key->GetName());
279  std::vector<std::string> tokens;
280  std::size_t pos = 0, prev = 0;
281  while(pos != std::string::npos) { // Loop over the TObject name
282  pos = name.find(sep, prev); // Find the next separator
283 
284  if(pos != std::string::npos) {
285  tokens.push_back(name.substr(prev, pos-prev));
286  prev = pos+2; // Move past the current separator
287  }
288  else {
289  tokens.push_back(name.substr(prev, name.length()-prev));
290  }
291  }
292 
293  if(tokens.size() > 1) {
294  cafanaType = tokens[0];
295  sampleLabel = tokens[1];
296  if(cafanaType.compare("pred") == 0) {
297  (*preds)[sampleLabel] = PredictionCombinePeriods::LoadFrom(keyDir).release();
298  }
299  }
300  else { continue; }
301  }
302 
303  tmp->cd();
304  return;
305 }
306 
307 //------------------------------------------------------------------------------
309  const IFitVar* var, std::vector<const IFitVar*> profVars,
310  TDirectory* rootOut, int nbins, double min, double max,
311  std::string name, Color_t color,
312  std::vector<const ISyst*> systs)
313 {
314  // Create slice histogram with/without profiling based on input profVars vector size
315 
316  TH1* h = Profile(expt, calc, var, nbins, min, max, -1, profVars, systs);
317 
318  ResetAngles(calc); // Don't forget to reset calculator...
319 
320  // Formatting:
321  h->SetLineColor(color); // Set the line color
322  h->SetMaximum(5); // Want the plot to go from 0 to 5
323  h->SetMinimum(0);
324  h->SetName(name.c_str()); // Set the object name for saving
325  rootOut->WriteTObject(h); // Save the raw histogram
326 
327  return h;
328 }
329 
330 //------------------------------------------------------------------------------
332  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24, const IFitVar* vard24,
333  TDirectory* rootOut, AngleValues a,
334  std::string years, Color_t color,
335  std::string samplelabel,
336  std::map<std::string, TH1*>* h24s,
337  std::map<std::string, TH1*>* h34s,
338  std::vector<const ISyst*> systs)
339 {
340  // General form of the object name for saving
341  std::string name = "h" + years + "yrTh";
342  std::string fullname; // fullname adds 2 characters, the theta subscript indices
343 
344  // Plot profiled slices in same order as above
345  // When slice is of one angle, the other two are profiled
346  // I.e., when the slice is theta 23, theta 34 and 24 are profiled
347  fullname = name + "34Prof" + samplelabel;
348  (*h34s)[samplelabel] = Plot1DSlice(
349  expt, calc, var34, {var23, var24}, rootOut,
350  a.nbins34, a.min34, a.max34, fullname, color, systs
351  );
352  fullname = name + "24Prof" + samplelabel;
353  (*h24s)[samplelabel] = Plot1DSlice(
354  expt, calc, var24, {var23, var34}, rootOut,
355  a.nbins24, a.min24, a.max24, fullname, color, systs
356  );
357 
358  return;
359 }
360 
361 //------------------------------------------------------------------------------
363  double th2468, double th2490,
364  double th3468, double th3490)
365 {
366  fprintf(text, "Angle limits for %s:\n", label.c_str());
367  fprintf(text, "Th24: 68% (%): %4.2f, 90%%: %4.2f\n", th2468, th2490);
368  fprintf(text, "Th34: 68% (%): %4.2f, 90%%: %4.2f\n", th3468, th3490);
369  fprintf(text, "\n");
370  return;
371 }
372 
373 //------------------------------------------------------------------------------
375  double th24stat68, double th2468,
376  double th24stat90, double th2490,
377  double th34stat68, double th3468,
378  double th34stat90, double th3490)
379 {
380  fprintf(text, "Diff24: 68% (%) %4.2f, 90%%: %4.2f\n",
381  CalcSystEffect(th24stat68, th2468),
382  CalcSystEffect(th24stat90, th2490));
383  fprintf(text, "Diff34: 68% (%) %4.2f, 90%%: %4.2f\n",
384  CalcSystEffect(th34stat68, th3468),
385  CalcSystEffect(th34stat90, th3490));
386  fprintf(text, "\n");
387  return;
388 }
389 
390 //------------------------------------------------------------------------------
392 {
394  calc->SetAngle(2, 3, M_PI/4); // 45 degrees (maximal mixing)
395  calc->SetAngle(3, 4, 0.);
396  calc->SetAngle(2, 4, 0.);
397 
398  return;
399 }
const XML_Char * name
Definition: expat.h:151
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, std::vector< const ISyst * > systs={})
enum BeamMode kRed
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetNFlavors(int nflavors)
double CalcSystEffect(double statonly, double withsyst)
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.
const NusFlatSyst kNusNDContSyst("ndcont","ND Containment", 1.0, 0.6)
Definition: NusSysts.h:92
Adapt the PMNS_Sterile calculator to standard interface.
const FitTheta23InDegreesSterile kFitTheta23InDegreesSterile
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
void PrintSystEffect(FILE *text, double th24stat68, double th2468, double th24stat90, double th2490, double th34stat68, double th3468, double th34stat90, double th3490)
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()
const NusSystFromHist kNusCalibFlatSyst(kNusAna01SystFile,"EX","CalFlat","Flat Miscalibration")
Definition: NusSysts.h:71
#define M_PI
Definition: SbMath.h:34
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const char * label
fclose(fg1)
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
const NusSystFromHist kNusBeamSysts(kNusAna01SystFile,"EX","Beam","All Beam")
Definition: NusSysts.h:69
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
std::vector< const ISyst * > getAllNusSysts()
Get a vector of all the nus group systs.
Definition: NusSysts.cxx:202
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
const NusSystFromHist kNusBirksSyst(kNusAna01SystFile,"EX","Birks","Birks C")
Definition: NusSysts.h:70
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)
const NusFlatSyst kNusMCStatsSyst("mcstat","MC Stats", 2.0, 4.8)
Definition: NusSysts.h:91
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 NusSystFromHist kNusGENIESmallSysts(kNusAna01SystFile,"EX","GENIESm","Summed small GENIE Systs")
Definition: NusSysts.h:75
const NusSystFromHist kNusNueCCSyst(kNusAna01SystFile,"EX","NueCC","#nu_{e} CC Background")
Definition: NusSysts.h:77
const char sep
OStream cout
Definition: OStream.cxx:6
void ResetAngles(osc::OscCalcSterile *calc)
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
Combine multiple component experiments.
double FindAngleFromChi2(double chi2, TH1 *h)
void SetAngle(int i, int j, double th)
void SetDm(int i, double dm)
TDirectory * dir
Definition: macro.C:5
const NusSystFromHist kNusCalibRelSyst(kNusAna01SystFile,"EX","CalRel","Relative Detector Miscalibration")
Definition: NusSysts.h:72
const NusFlatSyst kNusNormSyst("normNus","Normalization", 4.9, 4.9)
Definition: NusSysts.h:93
Base class defining interface for experiments.
Definition: IExperiment.h:14
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
Interface definition for fittable variables.
Definition: IFitVar.h:16
const NusSystFromHist kNusNDRockSyst(kNusAna01SystFile,"EX","NDNDRock","ND Rock")
Definition: NusSysts.h:76
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 FitSystEffectsAna()
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
const NusSystFromHist kNusCalibSlopeXSyst(kNusAna01SystFile,"EX","CalSlopeX","Sloped Miscalibration, X")
Definition: NusSysts.h:73
const NusSystFromHist kNusCalibSlopeYSyst(kNusAna01SystFile,"EX","CalSlopeY","Sloped Miscalibration, Y")
Definition: NusSysts.h:74
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.
const NusSystFromHist kNusNumuCCSyst(kNusAna01SystFile,"EX","NumuCC","#nu_{#mu} CC Background")
Definition: NusSysts.h:78
void PrintAngles(FILE *text, std::string label, double th2468, double th2490, double th3468, double th3490)
enum BeamMode string