Classes | Functions
BlessedPlotsAnaByPeriod.C File Reference
#include "TCanvas.h"
#include "TDirectory.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TLatex.h"
#include "TLegend.h"
#include "TLine.h"
#include "TPad.h"
#include "CAFAna/Analysis/Calcs.h"
#include "CAFAna/Fit/Fit.h"
#include "CAFAna/Analysis/Style.h"
#include "CAFAna/Fit/FrequentistSurface.h"
#include "CAFAna/Core/IFitVar.h"
#include "3FlavorAna/Cuts/NumuCuts.h"
#include "CAFAna/Cuts/TimingCuts.h"
#include "CAFAna/Decomp/ProportionalDecomp.h"
#include "CAFAna/Experiment/GaussianConstraint.h"
#include "CAFAna/Experiment/MultiExperiment.h"
#include "CAFAna/Experiment/CountingExperiment.h"
#include "CAFAna/Experiment/SingleSampleExperiment.h"
#include "CAFAna/Extrap/ExtrapSterile.h"
#include "CAFAna/Prediction/PredictionExtrap.h"
#include "CAFAna/Prediction/PredictionNoExtrap.h"
#include "NuXAna/Prediction/PredictionSterile.h"
#include "CAFAna/Prediction/PredictionCombinePeriods.h"
#include "NuXAna/Vars/FitVarsSterile.h"
#include "OscLib/OscCalcSterile.h"
#include "NuXAna/macros/NuSPlotFunctions.h"

Go to the source code of this file.

Classes

struct  AngleValues
 

Functions

void 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)
 
void Plot1DSlices (IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, TDirectory *rootOut, AngleValues a, Color_t color)
 
void Plot2DSlices (IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, TDirectory *rootOut, AngleValues a)
 
void Plot2DSlices (FrequentistSurface *s3424, FrequentistSurface *s2324, FrequentistSurface *s3423, TDirectory *rootOut, std::string prof="")
 
void PlotStack (Spectrum spectra[], TDirectory *rootOut, FILE *textOFS, std::string name, std::string title, std::string det, Spectrum *sterile=nullptr, Spectrum *sterileNC=nullptr)
 
void PlotStack (IDecomp *decomp, TDirectory *rootOut, FILE *textOFS, std::string name, std::string title)
 
void PlotStack (IPrediction *pred, Spectrum &sCos, osc::IOscCalc *calc, osc::IOscCalc *calcSt, TDirectory *rootOut, FILE *textOFS, std::string name, std::string title)
 
void PlotPurEff (IPrediction *pFDNum, IPrediction *pFDPurDen, IPrediction *pFDEffDen, Spectrum &sNDNum, Spectrum &sNDPurDen, Spectrum &sNDEffDen, osc::IOscCalc *calc, Spectrum &sCos, TDirectory *rootOut, FILE *textOFS, std::string name, std::string title)
 
void ResetAngles (osc::OscCalcSterile *calc)
 
void BlessedPlotsAnaByPeriod ()
 

Function Documentation

void BlessedPlotsAnaByPeriod ( )

Definition at line 156 of file BlessedPlotsAnaByPeriod.C.

References ana::DefaultSterileCalc(), demo5::expt, fclose(), ana::IFitter::Fit(), PandAna.Demos.tute_pid_validation::folder, ana::kFitTheta23InDegreesSterile, ana::kFitTheta24InDegreesSterile, ana::kFitTheta34InDegreesSterile, kRed, ana::PredictionCombinePeriods::LoadFrom(), ana::ProportionalDecomp::LoadFrom(), ana::Spectrum::LoadFrom(), AngleValues::max23, AngleValues::max24, AngleValues::max34, AngleValues::min23, AngleValues::min24, AngleValues::min34, AngleValues::nbins23, AngleValues::nbins24, AngleValues::nbins34, Plot1DSlices(), Plot2DSlices(), ana::PlotStack(), plot_validation_datamc::pred, ResetAngles(), osc::OscCalcSterile::SetAngle(), osc::OscCalcSterile::SetDm(), osc::OscCalcSterile::SetNFlavors(), and string.

157 {
158  TH1::AddDirectory(0);
159 
160  // Set up path to file with filled CAFAna objects and for output of analysis
161  //std::string folder = "/nova/ana/steriles/Ana01/";
162  std::string folder = "";
163  std::string filenm = "test_veto";
164 
165  std::string loadLocation = folder + filenm + ".root";
166  std::string saveLocation = folder + filenm + "Ana.root";
167  std::string textLocation = folder + filenm + "Ana.txt";
168 
169  // Load filled objects from file
170  TFile* rootL = new TFile(loadLocation.c_str(), "UPDATE");
171  TDirectory* tmpL = gDirectory;
172  TDirectory* loadDir = gDirectory;
173 
174  loadDir->cd((loadLocation + ":/decompNC").c_str());
175  ProportionalDecomp decompNC = *ProportionalDecomp::LoadFrom(gDirectory);
176 
177  loadDir->cd((loadLocation + ":/prediction").c_str());
179 
180  loadDir->cd((loadLocation + ":/sData").c_str());
181  Spectrum sFDData = *Spectrum::LoadFrom(gDirectory);
182 
183  loadDir->cd((loadLocation + ":/sCosSA").c_str());
184  Spectrum sCosmic = *Spectrum::LoadFrom(gDirectory);
185 
186  tmpL->cd();
187  rootL->Close(); // Don't forget to close the file!
188 
189  // Set up oscillation calculator that uses default 3 flavor parameters
191 
192  // Set up oscillation calculator that uses default 3 flavor parameters
194  calc4f->SetNFlavors(4);
195  calc4f->SetDm(4, 0.5); //Dm41 = 0.5 eV^2
196  calc4f->SetAngle(2, 4, 0.);
197  calc4f->SetAngle(3, 4, 0.);
198 
199  std::string labelEReco = "Calorimetric Energy (GeV)";
200 
201  // Open files for saving analysis output
202  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
203  FILE* textF;
204  textF = fopen(textLocation.c_str(), "w");
205 
206  // Choose some angle values for making 4 flavor predictions
207  calc4f->SetAngle(2, 4, 0.17453293); // 10 degrees
208  calc4f->SetAngle(3, 4, 0.61086524); // 35 degrees
209 
210  // Plot spectra for Ana01 prediction
211  std::string titleHelper = " Spectra for 6.05e20 POT-equiv.";
212  PlotStack(&pred, sCosmic, calc3f, calc4f, rootF, textF, "FDSpectra", titleHelper);
213  PlotStack(&decompNC, rootF, textF, "NDSpectra", titleHelper);
214 
215  fclose(textF); // Text file is done now, close it
216  ResetAngles(calc4f);
217 
218  // Create Experiment objects
219  GaussianConstraint th23Constraint(&kFitTheta23InDegreesSterile, 45., 8.73);
220  CountingExperiment expt(&pred, sFDData, sCosmic);
221  MultiExperiment multi({&th23Constraint, &expt});
222 
223  // Vector of parameters to fit: theta 23, 34, 24
224  std::vector<const IFitVar*> fitvars;
225  fitvars.push_back(&kFitTheta34InDegreesSterile);
226  fitvars.push_back(&kFitTheta24InDegreesSterile);
227  fitvars.push_back(&kFitTheta23InDegreesSterile);
228 
229  // Create fit objects
230  MinuitFitter fit(&multi, fitvars);
231 
232  // Run fits!
233  fit.Fit(calc4f);
234  ResetAngles(calc4f);
235 
236  // Set up values for slice and surface plots
237  AngleValues avals;
238  avals.nbins23 = 200;
239  avals.min23 = 20.;
240  avals.max23 = 70.;
241  avals.nbins34 = 200;
242  avals.min34 = 0.;
243  avals.max34 = 50.;
244  avals.nbins24 = 180;
245  avals.min24 = 0.;
246  avals.max24 = 45.;
247 
248  const Color_t kFitColor = kRed;
249 
250  // Plot 1D fits
251  Plot1DSlices(
252  &multi, calc4f,
256  rootF, avals, kFitColor
257  );
258 
259  // Lower the number of bins for surfaces
260  avals.nbins23 = 50;
261  avals.nbins34 = 50;
262  avals.nbins24 = 45;
263 
264  // Plot the surfaces
265  Plot2DSlices(
266  &multi, calc4f,
270  rootF, avals
271  );
272 
273  rootF->Close(); // Close the output root file
274 }
enum BeamMode kRed
void SetNFlavors(int nflavors)
cons_index_list< index_multi, nil_index_list > multi
A simple Gaussian constraint on an arbitrary IFitVar.
Adapt the PMNS_Sterile calculator to standard interface.
const FitTheta23InDegreesSterile kFitTheta23InDegreesSterile
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
fclose(fg1)
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 ResetAngles(osc::OscCalcSterile *calc)
void PlotStack(Spectrum spectra[], TDirectory *rootOut, FILE *textOFS, std::string name, std::string title, std::string det, Spectrum *sterile=nullptr, Spectrum *sterileNC=nullptr)
void Plot1DSlices(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, TDirectory *rootOut, AngleValues a, Color_t color)
Combine multiple component experiments.
void SetAngle(int i, int j, double th)
Splits Data proportionally according to MC.
void SetDm(int i, double dm)
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
Sum MC predictions from different periods scaled according to data POT targets.
void Plot2DSlices(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, TDirectory *rootOut, AngleValues a)
Compare a data spectrum to MC expectation using only the event count.
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string
void 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 
)

Definition at line 277 of file BlessedPlotsAnaByPeriod.C.

References plot_validation_datamc::c, om::cout, allTimeWatchdog::endl, make_syst_table_plots::h, ana::Profile(), ResetAngles(), Simulation(), ana::Slice(), string, and tex.

Referenced by Plot1DSlices().

281 {
282  // Create slice histogram with/without profiling based on input profVars vector size
283 
284  TH1* h = (profVars.size() > 0 ?
285  Profile(expt, calc, var, nbins, min, max, -1, profVars) :
286  Slice(expt, calc, var, nbins, min, max));
287 
288  ResetAngles(calc); // Don't forget to reset calculator...
289 
290  // Formatting:
291  h->SetLineColor(color); // Set the line color
292  h->SetMaximum(5); // Want the plot to go from 0 to 5
293  h->SetMinimum(0);
294  h->SetName(name.c_str()); // Set the object name for saving
295  rootOut->WriteTObject(h); // Save the raw histogram
296 
297  // Create an object for writing LaTeX to canvas
298  TLatex *tex = new TLatex();
299  tex->SetNDC();
300  tex->SetTextFont(42);
301  tex->SetTextSize(0.037);
302  tex->SetTextAlign(11);
303 
304  // Coordinates for LaTeX object
305  double xTex = 0.92, y68 = 0.25, y90 = 0.52;
306 
307  // Dotted lines representing 1 and 2 sigma
308  TLine* line68 = new TLine(min, 1, max, 1);
309  TLine* line90 = new TLine(min, 2.71, max, 2.71);
310 
311  // Sometimes histogram has horizontal artifact at max canvas value
312  // This line hides it
313  TLine* line5 = new TLine(min, 5, max, 5);
314 
315  // Dashed line for 1 and 2 sigma lines
316  line68->SetLineStyle(2);
317  line90->SetLineStyle(2);
318 
319  // Formatting
320  line68->SetLineWidth(2);
321  line90->SetLineWidth(2);
322  line5 ->SetLineWidth(2);
323 
324  // Input "name" should have "h" for histogram at the beginning
325  // Make canvas name with input "name" but remove that "h"
326  std::cout << name << std::endl;
327  std::string cname = "c1DSlice" + name.substr(1, name.size()-1);
328 
329  // Theta subscripts should be at position 4, 5, get them
330  std::string indices = name.substr(4, 2);
331 
332  // Create POT string
333  std::string cpotstr = "6.05e20 POT-equiv.";
334 
335  // Create a title for canvas
336  std::string ctitle = "Delta chi^2 for Theta" + indices;
337 
338  // Add info on profiling if applicable
339  ctitle += " (Profiling Over Other Angles) for " + cpotstr;
340 
341  // Create the canvas and draw!
342  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 600, 500);
343  h->Draw("]["); // Don't show vertical drop to 0 at edge of histogram
344  line68->DrawClone("same");
345  line90->DrawClone("same");
346  line5 ->DrawClone("same");
347  tex->DrawLatex(xTex, y68, "68%"); // Add label for 1 sigma line
348  tex->DrawLatex(xTex, y90, "90%"); // Add label for 2 sigma line
349  Simulation();
350  rootOut->WriteTObject(c); // Save full canvas
351 
352  return;
353 }
void Simulation()
Definition: tools.h:16
const XML_Char * name
Definition: expat.h:151
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
const int nbins
Definition: cellShifts.C:15
void ResetAngles(osc::OscCalcSterile *calc)
TGraph * Slice(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, MinuitFitter::FitOpts opts)
scan in one variable, holding all others constant
Definition: Fit.cxx:169
OStream cout
Definition: OStream.cxx:6
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
TLatex * tex
Definition: f2_nu.C:499
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
enum BeamMode string
void Plot1DSlices ( IExperiment expt,
osc::OscCalcSterile calc,
const IFitVar var23,
const IFitVar var34,
const IFitVar var24,
TDirectory *  rootOut,
AngleValues  a,
Color_t  color 
)

Definition at line 356 of file BlessedPlotsAnaByPeriod.C.

References plot_validation_datamc_2018::color, AngleValues::max23, AngleValues::max24, AngleValues::max34, AngleValues::min23, AngleValues::min24, AngleValues::min34, AngleValues::nbins23, AngleValues::nbins24, AngleValues::nbins34, Plot1DSlice(), and string.

Referenced by BlessedPlotsAnaByPeriod().

360 {
361  // General form of the object name for saving
362  std::string name = "hTh";
363  std::string fullname; // fullname adds 2 characters, the theta subscript indices
364 
365  // Plot unprofiled slices, for theta 23, then 34, then 24
366  fullname = name + "23";
367  Plot1DSlice(expt, calc, var23, {},
368  rootOut, a.nbins23, a.min23, a.max23, fullname, color);
369  fullname = name + "34";
370  Plot1DSlice(expt, calc, var34, {},
371  rootOut, a.nbins34, a.min34, a.max34, fullname, color);
372  fullname = name + "24";
373  Plot1DSlice(expt, calc, var24, {},
374  rootOut, a.nbins24, a.min24, a.max24, fullname, color);
375 
376 
377  // Plot profiled slices in same order as above
378  // When slice is of one angle, the other two are profiled
379  // I.e., when the slice is theta 23, theta 34 and 24 are profiled
380  fullname = name + "23Prof";
381  Plot1DSlice(expt, calc, var23, {var34, var24},
382  rootOut, a.nbins23, a.min23, a.max23, fullname, color);
383  fullname = name + "34Prof";
384  Plot1DSlice(expt, calc, var34, {var23, var24},
385  rootOut, a.nbins34, a.min34, a.max34, fullname, color);
386  fullname = name + "24Prof";
387  Plot1DSlice(expt, calc, var24, {var23, var34},
388  rootOut, a.nbins24, a.min24, a.max24, fullname, color);
389 
390  return;
391 }
const XML_Char * name
Definition: expat.h:151
void 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)
enum BeamMode string
void Plot2DSlices ( IExperiment expt,
osc::OscCalcSterile calc,
const IFitVar var23,
const IFitVar var34,
const IFitVar var24,
TDirectory *  rootOut,
AngleValues  a 
)

Definition at line 394 of file BlessedPlotsAnaByPeriod.C.

References AngleValues::max23, AngleValues::max24, AngleValues::max34, AngleValues::min23, AngleValues::min24, AngleValues::min34, AngleValues::nbins23, AngleValues::nbins24, AngleValues::nbins34, and ResetAngles().

Referenced by BlessedPlotsAnaByPeriod().

397 {
398  // Create unprofiled surfaces
399  // 34 vs 24
400  FrequentistSurface s3424(expt, calc,
401  var34, a.nbins34, a.min34, a.max34,
402  var24, a.nbins24, a.min24, a.max24);
403  ResetAngles(calc);
404 
405  // 23 vs 24
406  FrequentistSurface s2324(expt, calc,
407  var23, a.nbins23, a.min23, a.max23,
408  var24, a.nbins24, a.min24, a.max24);
409  ResetAngles(calc);
410 
411  // 34 vs 23
412  FrequentistSurface s3423(expt, calc,
413  var34, a.nbins34, a.min34, a.max34,
414  var23, a.nbins23, a.min23, a.max23);
415  ResetAngles(calc);
416 
417  // Create profiled surfaces in same order as above
418  // The third angle is the one that is profiled
419  // I.e., the 34 vs 24 surface profiles 23
420  FrequentistSurface s3424Prof(expt, calc,
421  var34, a.nbins34, a.min34, a.max34,
422  var24, a.nbins24, a.min24, a.max24,
423  {var23});
424  ResetAngles(calc);
425 
426  FrequentistSurface s2324Prof(expt, calc,
427  var23, a.nbins23, a.min23, a.max23,
428  var24, a.nbins24, a.min24, a.max24,
429  {var34});
430  ResetAngles(calc);
431 
432  FrequentistSurface s3423Prof(expt, calc,
433  var34, a.nbins34, a.min34, a.max34,
434  var23, a.nbins23, a.min23, a.max23,
435  {var24});
436  ResetAngles(calc);
437 
438  // Plot all of the surfaces
439  Plot2DSlices(&s3424, &s2324, &s3423, rootOut);
440  Plot2DSlices(&s3424Prof, &s2324Prof, &s3423Prof, rootOut, "Prof");
441 
442  return;
443 }
Log-likelihood scan across two parameters.
void ResetAngles(osc::OscCalcSterile *calc)
void Plot2DSlices(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, TDirectory *rootOut, AngleValues a)
void Plot2DSlices ( FrequentistSurface s3424,
FrequentistSurface s2324,
FrequentistSurface s3423,
TDirectory *  rootOut,
std::string  prof = "" 
)

Definition at line 446 of file BlessedPlotsAnaByPeriod.C.

References plot_validation_datamc::c, ana::ISurface::Draw(), ana::ISurface::DrawBestFit(), ana::ISurface::DrawContour(), ana::Gaussian68Percent2D(), ana::Gaussian90Percent2D(), Simulation(), string, and ana::ISurface::ToTH2().

448 {
449  // General format for the name for saving
450  std::string name = "hTh";
451  std::string fullname;
452 
453  // Add indices for first theta, add second theta, and profile string if applicaable
454  fullname = name + "34Th24" + prof;
455  TH1* h3424 = s3424->ToTH2(); // Create a histogram from the surface
456  h3424->SetName(fullname.c_str()); // Set the name for saving
457  rootOut->WriteTObject(h3424); // Save the raw histogram
458 
459  fullname = name + "23Th24" + prof;
460  TH1* h2324 = s2324->ToTH2();
461  h2324->SetName(fullname.c_str());
462  rootOut->WriteTObject(h2324);
463 
464  fullname = name + "34Th23" + prof;
465  TH1* h3423 = s3423->ToTH2();
466  h3423->SetName(fullname.c_str());
467  rootOut->WriteTObject(h3423);
468 
469  // Create name for canvas
470  std::string cname = "c2DSlice" + prof;
471  // Create string for POT based on whether name has 1 or 3 in it
472  std::string potstr = "6.05e20 POT-equiv.";
473  // Create relevant title for the canvas
474  std::string ctitle = (
475  prof.size() > 0 ?
476  "Surfaces (Profiling Over Other Angles) for " + potstr :
477  "Surfaces for " + potstr
478  );
479 
480  // Create triptych canvas and draw!
481  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 800, 800);
482  c->Divide(2, 2); // Split into 2 x 2 grid
483  c->cd(1); // Go to top left
484  s3424->Draw(); // Draw histogram
485  s3424->DrawContour(Gaussian68Percent2D(*s3424), 7, kBlack); // Draw 1 sigma contour
486  s3424->DrawContour(Gaussian90Percent2D(*s3424), 1, kBlack); // Draw 2 sigma contour
487  s3424->DrawBestFit(kBlack); // Draw best fit point
488  c->cd(2); // Go to top right
489  s2324->Draw();
490  s2324->DrawContour(Gaussian68Percent2D(*s2324), 7, kBlack);
491  s2324->DrawContour(Gaussian90Percent2D(*s2324), 1, kBlack);
492  s2324->DrawBestFit(kBlack);
493  Simulation();
494  c->cd(3); // Go to bottom left
495  s3423->Draw();
496  s3423->DrawContour(Gaussian68Percent2D(*s3423), 7, kBlack);
497  s3423->DrawContour(Gaussian90Percent2D(*s3423), 1, kBlack);
498  s3423->DrawBestFit(kBlack);
499  rootOut->WriteTObject(c); // Save the root canvas
500 
501  return;
502 }
void Simulation()
Definition: tools.h:16
const XML_Char * name
Definition: expat.h:151
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:37
void Draw() const
Draw the surface itself.
Definition: ISurface.cxx:21
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:50
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:170
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
enum BeamMode string
void PlotPurEff ( IPrediction pFDNum,
IPrediction pFDPurDen,
IPrediction pFDEffDen,
Spectrum sNDNum,
Spectrum sNDPurDen,
Spectrum sNDEffDen,
osc::IOscCalc calc,
Spectrum sCos,
TDirectory *  rootOut,
FILE *  textOFS,
std::string  name,
std::string  title 
)

Definition at line 785 of file BlessedPlotsAnaByPeriod.C.

References plot_validation_datamc::c, ana::CenterTitles(), kBlue, kRed, MECModelEnuComparisons::leg, ana::NCSCALE, ana::IPrediction::Predict(), Simulation(), and ana::Spectrum::ToTH1().

790 {
791  Spectrum sFDNum = pFDNum ->Predict(calc);
792  Spectrum sFDPurDen = pFDPurDen->Predict(calc) + sCos;
793  Spectrum sFDEffDen = pFDEffDen->Predict(calc);
794 
795  TH1* hFDNum = sFDNum .ToTH1(NCSCALE);
796  TH1* hFDPurDen = sFDPurDen.ToTH1(NCSCALE);
797  TH1* hFDEffDen = sFDEffDen.ToTH1(NCSCALE);
798 
799  TH1* hNDNum = sNDNum .ToTH1(NCSCALE);
800  TH1* hNDPurDen = sNDPurDen.ToTH1(NCSCALE);
801  TH1* hNDEffDen = sNDEffDen.ToTH1(NCSCALE);
802 
803  TH1* hFDEff = (TH1*)hFDNum->Clone();
804  hFDEff->Divide(hFDEffDen);
805 
806  TH1* hNDEff = (TH1*)hNDNum->Clone();
807  hNDEff->Divide(hNDEffDen);
808 
809  TH1* hFDPur = (TH1*)hFDNum->Clone();
810  hFDPur->Divide(hFDPurDen);
811 
812  TH1* hNDPur = (TH1*)hNDNum->Clone();
813  hNDPur->Divide(hNDPurDen);
814 
815  hFDEff->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
816  hFDPur->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
817  hNDEff->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
818  hNDPur->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
819 
820  hFDEff->SetLineColor(kBlue);
821  hNDEff->SetLineColor(kBlue);
822 
823  hFDPur->SetLineColor(kRed);
824  hNDPur->SetLineColor(kRed);
825 
826  hNDEff->SetLineStyle(2);
827  hNDPur->SetLineStyle(2);
828 
829  hFDEff->SetMaximum(1.);
830  hFDEff->SetMinimum(0.);
831  hFDEff->GetXaxis()->SetNdivisions(-406);
832 
833  hNDEff->SetMaximum(1.);
834  hNDEff->SetMinimum(0.);
835 
836  hFDPur->SetMaximum(1.);
837  hFDPur->SetMinimum(0.);
838 
839  hNDPur->SetMaximum(1.);
840  hNDPur->SetMinimum(0.);
841 
842  CenterTitles(hFDEff);
843  CenterTitles(hNDEff);
844  CenterTitles(hFDPur);
845  CenterTitles(hNDPur);
846 
847  TLegend* leg = new TLegend(0.6, 0.35, 0.85, 0.55);
848  leg->AddEntry(hFDEff, "FD MC Efficiency", "L");
849  leg->AddEntry(hNDEff, "ND MC Efficiency", "L");
850  leg->AddEntry(hFDPur, "FD MC Purity", "L");
851  leg->AddEntry(hNDPur, "ND MC Purity", "L");
852 
853  TCanvas* c = new TCanvas(name.c_str(), "Purity and Efficiency", 800, 500);
854  hFDEff->Draw("hist");
855  hNDEff->Draw("hist same");
856  hFDPur->Draw("hist same");
857  hNDPur->Draw("hist same");
858  leg->Draw();
859 
860  gPad->RedrawAxis();
861  Simulation();
862 
863  rootOut->WriteTObject(c);
864 }
void Simulation()
Definition: tools.h:16
const XML_Char * name
Definition: expat.h:151
enum BeamMode kRed
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
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void CenterTitles(TH1 *h)
static double NCSCALE
enum BeamMode kBlue
void PlotStack ( Spectrum  spectra[],
TDirectory *  rootOut,
FILE *  textOFS,
std::string  name,
std::string  title,
std::string  det,
Spectrum sterile = nullptr,
Spectrum sterileNC = nullptr 
)

Definition at line 505 of file BlessedPlotsAnaByPeriod.C.

References plot_validation_datamc::c, ana::CenterTitles(), MECModelEnuComparisons::i, ana::kCosmicBackgroundColor, ana::kNCBackgroundColor, ana::kNueSignalColor, ana::kNumuBackgroundColor, kOrange, ana::kTotalMCColor, MECModelEnuComparisons::leg, std::max(), getGoodRuns4SAM::n, POT, ana::Spectrum::POT(), Simulation(), std::sqrt(), string, tex, and ana::Spectrum::ToTH1().

508 {
509 
510 
511  int POT = 6.05;
512  std::string potStr ="6.05e20 POT-equiv.";
513 
514  // Get each histogram from the Spectrum array
515  TH1* hAll = spectra[0].ToTH1(spectra[0].POT());
516  TH1* hNC = spectra[1].ToTH1(spectra[1].POT());
517  TH1* hNumu = spectra[2].ToTH1(spectra[2].POT());
518  TH1* hNue = spectra[3].ToTH1(spectra[3].POT());
519  TH1* hCos = spectra[4].ToTH1(spectra[4].POT());
520  // The default for these inputs are nullptr, only get a TH1 if they exist
521  TH1* hStrl = (sterile ? sterile ->ToTH1(sterile->POT()) : nullptr);
522  TH1* hStrlNC = (sterileNC ? sterileNC->ToTH1(sterileNC->POT()) : nullptr);
523 
524  // Create the title and axis labels, which is an empty title,
525  // the x label from an input Spectrum, and Events/POT
526  std::string xLabel(hAll->GetXaxis()->GetTitle());
527  std::string yLabel = (det.compare("ND") == 0 ? "10^{3} Events" : "Events");
528  yLabel += "/" + potStr + " / 0.25 GeV";
529  std::string histTitle = ";" + xLabel + ";" + yLabel;
530 
531  // Print the column title and headers
532  fprintf(textOFS, "Event numbers at the %2.2s for %2de20 POT-equiv.:\n", det.c_str(), POT);
533  if(hStrlNC) { fprintf(textOFS, "%12.12s, ", "NC (Sterile)"); }
534  else { fprintf(textOFS, "%14.14s", ""); }
535  fprintf(textOFS, "%12.12s, %12.12s, %12.12s, %12.12s, %12.12s, %12.12s\n",
536  "NC (3 Flav)", "All", "Numu", "Nue", "Cosmic", "FOM");
537 
538  // Get the integral number of events, then print them
539  double nNC3F = hNC ->Integral();
540  double nAll = hAll ->Integral();
541  double nNumu = hNumu->Integral();
542  double nNue = hNue ->Integral();
543  double nCos = hCos ->Integral();
544  double fom = nNC3F/sqrt(nAll);
545  if(hStrlNC) {
546  double nNCSt = hStrlNC->Integral();
547  fprintf(textOFS, "%10E, ", nNCSt);
548  }
549  else {
550  fprintf(textOFS, "%14.14s", "");
551  }
552  fprintf(textOFS, "%10E, %10E, %10E, %10E, %10E, %10E\n\n",
553  nNC3F, nAll, nNumu, nNue, nCos, fom);
554 
555  // Stylize the unstacked plots
556  CenterTitles(hNue);
557  CenterTitles(hNumu);
558  CenterTitles(hNC);
559  CenterTitles(hCos);
560  CenterTitles(hAll);
561 
562  hAll ->SetMinimum(0);
563  hNC ->SetMinimum(0);
564  hNumu->SetMinimum(0);
565  hNue ->SetMinimum(0);
566  hCos ->SetMinimum(0);
567  if(hStrlNC) { hStrl->SetMinimum(0); }
568 
569  const Color_t kCosmicBackgroundColor = kOrange+7;
570  hAll ->SetLineColor(kTotalMCColor);
571  hNC ->SetLineColor(kNCBackgroundColor);
572  hNumu->SetLineColor(kNumuBackgroundColor);
573  hNue ->SetLineColor(kNueSignalColor);
574  hCos ->SetLineColor(kCosmicBackgroundColor);
575  if(hStrlNC) {
576  hStrlNC->SetLineColor(kNCBackgroundColor);
577  hStrlNC->SetLineStyle(2);
578  }
579 
580  hAll ->SetTitle(histTitle.c_str());
581  hNC ->SetTitle(histTitle.c_str());
582  hNumu->SetTitle(histTitle.c_str());
583  hNue ->SetTitle(histTitle.c_str());
584  hCos ->SetTitle(histTitle.c_str());
585  if(hStrlNC) { hStrlNC->SetTitle(histTitle.c_str()); }
586 
587  // Create "stacked" copies of flavored histograms
588  // Cos Stack = Cos
589  TH1* hCosStack = (TH1*)hCos->Clone();
590 
591  // Numu Stack = Cos Stack + Numu = Cos + Numu
592  TH1* hNumuStack = (TH1*)hCosStack->Clone();
593  hNumuStack->Add(hNumu);
594 
595  // Nue Stack = Numu Stack + Nue = Cos + Numu + Nue
596  TH1* hNueStack = (TH1*)hNumuStack->Clone();
597  hNueStack->Add(hNue);
598 
599  // NC Stack = Nue Stack + NC = Cos + Numu + Nue + NC
600  TH1* hNCStack = (TH1*)hNC->Clone();
601  //hNCStack->Add(hNC);
602 
603  // Stylize the stacked plots
604  CenterTitles(hNueStack);
605  CenterTitles(hNumuStack);
606  CenterTitles(hNCStack);
607  CenterTitles(hCosStack);
608  if(hStrlNC) { CenterTitles(hStrlNC); }
609 
610  hNCStack ->SetLineColor(kNCBackgroundColor);
611  hNumuStack->SetLineColor(kNumuBackgroundColor);
612  hNumuStack->SetFillColor(kNumuBackgroundColor);
613  hNueStack ->SetLineColor(kNueSignalColor);
614  hNueStack ->SetFillColor(kNueSignalColor);
615  hCosStack ->SetLineColor(kCosmicBackgroundColor);
616  hCosStack ->SetFillColor(kCosmicBackgroundColor);
617 
618  hNCStack ->SetMinimum(0);
619  hNumuStack->SetMinimum(0);
620  hNueStack ->SetMinimum(0);
621  hCosStack ->SetMinimum(0);
622 
623  // Manually set errors, and calculate the max value (including errors)
624  double maxval = 0.;
625  double maxvalstack = 0.;
626  for(int i = 1, n = hNC->GetNbinsX(); i <= n; ++i) {
627  double error = sqrt(hNC->GetBinContent(i));
628  hNC ->SetBinError(i, error);
629  hNCStack->SetBinError(i, error);
630  maxval = std::max(maxval, std::max(hAll->GetBinContent(i), hNC->GetBinContent(i) + error));
631  maxvalstack = std::max(maxvalstack, hNCStack->GetBinContent(i) + error);
632  }
633 
634  // Scale down the ND by 10e3 (don't forget the max values!)
635  if(det.compare("ND") == 0) {
636  double ndScale = 1./1000.;
637  hAll ->Scale(ndScale);
638  hNC ->Scale(ndScale);
639  hNue ->Scale(ndScale);
640  hNumu->Scale(ndScale);
641  hNCStack ->Scale(ndScale);
642  hNueStack ->Scale(ndScale);
643  hNumuStack->Scale(ndScale);
644  maxval *= ndScale;
645  maxvalstack *= ndScale;
646  }
647 
648  TLatex *tex = new TLatex();
649  tex->SetNDC();
650  tex->SetTextFont(42);
651  tex->SetTextSize(0.037);
652  tex->SetTextAlign(11);
653 
654  // Create a canvas for and draw the unstacked histograms
655  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 500);
656  hNue->SetMaximum(maxval*1.1);
657  hNue->GetXaxis()->SetNdivisions(-406);
658 
659  // Create a legend for the unstacked histograms
660  double xL = 0.5, xR = 0.7, yB = 0.59, yT = 0.84;
661  if(det.compare("ND") == 0) { yB += 0.05; }
662  TLegend* leg = new TLegend(xL, yB, xR, yT);
663  leg->SetBorderSize(0);
664  leg->SetFillColor(kWhite);
665  leg->SetFillStyle(0);
666  leg->SetFillStyle(0);
667  leg->SetTextFont(42);
668  leg->SetTextSize(0.037);
669  leg->AddEntry(hAll, "Total Prediction", "l");
670  leg->AddEntry(hNC, "NC 3 Flavor Prediction", "le");
671  leg->AddEntry(hNue, "#nu_{e} CC Background", "l");
672  leg->AddEntry(hNumu, "#nu_{#mu} CC Background", "l");
673  if(det.compare("FD") == 0) {
674  leg->AddEntry(hCos, "Cosmic Background", "l");
675  }
676 
677  hNue ->Draw("hist");
678  hNumu->Draw("hist same");
679  if(det.compare("FD") == 0) { hCos->Draw("hist same"); }
680  hAll ->Draw("hist same");
681  hNC ->Draw("same");
682  leg->Draw();
683  if(det.compare("FD") == 0) {
684  tex->DrawLatex(0.505, yB - 0.045, "#Deltam^{2}_{41} = 0.5 eV^{2}, #theta_{24} = 10#circ, #theta_{34} = 35#circ");
685  tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.37x10^{-3} eV^{2}, #theta_{13} = 8.5#circ, #theta_{23} = 45#circ");
686  tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
687  }
688 
689  gPad->RedrawAxis();
690  Simulation();
691 
692  // Save the canvas
693  rootOut->WriteTObject(c);
694 
695  // Create a canvas for and draw the stacked histograms
696  TCanvas* cStack = new TCanvas((name + "Stack").c_str(), title.c_str(), 800, 500);
697  hNueStack->SetMaximum(maxvalstack*1.1);
698  hNueStack->GetXaxis()->SetNdivisions(-406);
699 
700  if(det.compare("ND") == 0) { yB += 0.05; }
701  TLegend* legStack = new TLegend(xL, yB, xR, yT);
702  legStack->SetBorderSize(0);
703  legStack->SetFillColor(kWhite);
704  legStack->SetFillStyle(0);
705  legStack->SetFillStyle(0);
706  legStack->SetTextFont(42);
707  legStack->SetTextSize(0.037);
708  legStack->AddEntry(hNCStack, "NC 3 Flavor Prediction", "le");
709  if(hStrlNC) { legStack->AddEntry(hStrlNC, "NC 4 Flavor Prediction", "l"); }
710  legStack->AddEntry(hNueStack, "#nu_{e} CC Background", "f");
711  legStack->AddEntry(hNumuStack, "#nu_{#mu} CC Background", "f");
712  if(det.compare("FD") == 0) {
713  legStack->AddEntry(hCosStack, "Cosmic Background", "f");
714  }
715 
716  hNueStack ->Draw("hist");
717  hNumuStack->Draw("hist same");
718  if(det.compare("FD") == 0) { hCosStack->Draw("hist same"); }
719  hNCStack ->Draw("same");
720  if(hStrlNC) { hStrlNC->Draw("hist same"); }
721  legStack->Draw();
722  if(det.compare("FD") == 0) {
723  tex->DrawLatex(0.505, yB - 0.045, "#Deltam^{2}_{41} = 0.5 eV^{2}, #theta_{24} = 10#circ, #theta_{34} = 35#circ");
724  tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.37x10^{-3} eV^{2}, #theta_{13} = 8.5#circ, #theta_{23} = 45#circ");
725  tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
726  }
727 
728  gPad->RedrawAxis();
729  Simulation();
730 
731  rootOut->WriteTObject(cStack);
732 
733  return;
734 }
void Simulation()
Definition: tools.h:16
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
enum BeamMode kOrange
const Color_t kCosmicBackgroundColor
Definition: Style.h:41
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
T sqrt(T number)
Definition: d0nt_math.hpp:156
const Color_t kNumuBackgroundColor
Definition: Style.h:30
const Color_t kNueSignalColor
Definition: Style.h:19
void CenterTitles(TH1 *h)
double POT() const
Definition: Spectrum.h:227
std::vector< double > POT
TLatex * tex
Definition: f2_nu.C:499
const Color_t kTotalMCColor
Definition: Style.h:16
const Color_t kNCBackgroundColor
Definition: Style.h:22
enum BeamMode string
void PlotStack ( IDecomp decomp,
TDirectory *  rootOut,
FILE *  textOFS,
std::string  name,
std::string  title 
)

Definition at line 737 of file BlessedPlotsAnaByPeriod.C.

References ana::IDecomp::AntiNueComponent(), ana::IDecomp::AntiNumuComponent(), ana::Spectrum::Clear(), fillBadChanDBTables::det, MAXSPEC, ana::IDecomp::NCTotalComponent(), ana::IDecomp::NueComponent(), ana::IDecomp::NumuComponent(), ana::PlotStack(), string, and plotROC::title.

739 {
740  Spectrum empty = decomp->NCTotalComponent();
741  empty.Clear();
742 
743  Spectrum spectraND[MAXSPEC] = {
744  decomp->NCTotalComponent() +
745  decomp->NueComponent() + decomp->AntiNueComponent() +
746  decomp->NumuComponent() + decomp->AntiNumuComponent(),
747  decomp->NCTotalComponent(),
748  decomp->NumuComponent() + decomp->AntiNumuComponent(),
749  decomp->NueComponent() + decomp->AntiNueComponent(),
750  empty
751  };
752 
753  std::string det = "ND";
754  std::string fulltitle = det + title;
755  PlotStack(spectraND, rootOut, textOFS, name, fulltitle, det);
756 
757  return;
758 }
const XML_Char * name
Definition: expat.h:151
virtual Spectrum AntiNueComponent() const =0
virtual Spectrum NumuComponent() const =0
#define MAXSPEC
Definition: NusLoadProd3.h:18
void Clear()
Definition: Spectrum.cxx:361
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void PlotStack(Spectrum spectra[], TDirectory *rootOut, FILE *textOFS, std::string name, std::string title, std::string det, Spectrum *sterile=nullptr, Spectrum *sterileNC=nullptr)
virtual Spectrum AntiNumuComponent() const =0
virtual Spectrum NCTotalComponent() const
Definition: IDecomp.h:18
virtual Spectrum NueComponent() const =0
enum BeamMode string
void PlotStack ( IPrediction pred,
Spectrum sCos,
osc::IOscCalc calc,
osc::IOscCalc calcSt,
TDirectory *  rootOut,
FILE *  textOFS,
std::string  name,
std::string  title 
)

Definition at line 761 of file BlessedPlotsAnaByPeriod.C.

References fillBadChanDBTables::det, ana::Flavors::kAll, ana::Flavors::kAllNuE, ana::Flavors::kAllNuMu, ana::Sign::kBoth, ana::Current::kCC, ana::Current::kNC, MAXSPEC, ana::PlotStack(), ana::IPrediction::Predict(), ana::IPrediction::PredictComponent(), string, and plotROC::title.

765 {
766  Spectrum spectraFD[MAXSPEC] = {
767  pred->Predict(calc) + sCos,
771  sCos
772  };
773 
774  Spectrum sterile = pred->Predict(calcSt) + sCos;
775  Spectrum sterileNC = pred->PredictComponent(calcSt, Flavors::kAll, Current::kNC, Sign::kBoth);
776 
777  std::string det = "FD";
778  std::string fulltitle = det + title;
779  PlotStack(spectraFD, rootOut, textOFS, name, fulltitle, det, &sterile, &sterileNC);
780 
781  return;
782 }
const XML_Char * name
Definition: expat.h:151
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
#define MAXSPEC
Definition: NusLoadProd3.h:18
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Charged-current interactions.
Definition: IPrediction.h:39
void PlotStack(Spectrum spectra[], TDirectory *rootOut, FILE *textOFS, std::string name, std::string title, std::string det, Spectrum *sterile=nullptr, Spectrum *sterileNC=nullptr)
All neutrinos, any flavor.
Definition: IPrediction.h:26
enum BeamMode string
void ResetAngles ( osc::OscCalcSterile calc)

Definition at line 868 of file BlessedPlotsAnaByPeriod.C.

References M_PI, and osc::OscCalcSterile::SetAngle().

Referenced by BlessedPlotsAnaByPeriod(), Plot1DSlice(), and Plot2DSlices().

869 {
870  calc->SetAngle(2, 3, M_PI/4); // 45 degrees (maximal mixing)
871  calc->SetAngle(3, 4, 0.);
872  calc->SetAngle(2, 4, 0.);
873 
874  return;
875 }
#define M_PI
Definition: SbMath.h:34
void SetAngle(int i, int j, double th)