BlessedPlotsAna.C
Go to the documentation of this file.
1 #include "TCanvas.h"
2 #include "TDirectory.h"
3 #include "TFile.h"
4 #include "TH1.h"
5 #include "TH2.h"
6 #include "TLatex.h"
7 #include "TLegend.h"
8 #include "TLine.h"
9 #include "TPad.h"
10 
11 #include "CAFAna/Analysis/Calcs.h"
12 #include "CAFAna/Fit/Fit.h"
13 #include "CAFAna/Analysis/Style.h"
15 #include "CAFAna/Core/IFitVar.h"
17 #include "CAFAna/Cuts/TimingCuts.h"
28 
29 #include "OscLib/OscCalcSterile.h"
30 
32 
33 using namespace ana;
34 
35 // Definitions for organization
36 
37 // One struct to contain all values for angle fits
38 struct AngleValues{
39  int nbins23;
40  double min23;
41  double max23;
42 
43  int nbins34;
44  double min34;
45  double max34;
46 
47  int nbins24;
48  double min24;
49  double max24;
50 };
51 
52 // Plot a single 1D slice canvas -- this function does the actual plotting
53 // The macro does not call this function directly in the main function
54 // \param expt The CAFAna experiment object, contains FD prediction and "data"
55 // \param calc Sterile calculator for oscillating the spectra
56 // \param var Parameter to fit for!
57 // \param profVars Vector of parameters to profile over
58 // \param rootOut Already open file for root output
59 // \param name Base of the string for creating the name for saving the root objects
61  const IFitVar* var, std::vector<const IFitVar*> profVars,
62  TDirectory* rootOut, int nbins, double min, double max,
63  std::string name, Color_t color);
64 
65 // Plot 1D slice canvases
66 // This function sets up multiple calls to the single canvas version of Plot1DSlice
67 // \param expt The CAFAna experiment object, contains FD prediction and "data"
68 // \param calc Sterile calculator for oscillating the spectra
69 // \param rootOut Already open file for root output
70 // \param a Struct containing bin values for slice histograms
71 // \param years A string for the number of data years
73  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
74  TDirectory* rootOut, AngleValues a,
75  std::string years, Color_t color);
76 
77 // Plot 2D slice surfaces
78 // This function sets up multiple calls to the workhorse version of Plot2DSlices
79 // \param expt The CAFAna experiment object, contains FD prediction and "data"
80 // \param calc Sterile calculator for oscillating the spectra
81 // \param rootOut Already open file for root output
82 // \param a Struct containing bin values for slice histograms
83 // \param years A string for the number of data years
85  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
86  TDirectory* rootOut, AngleValues a, std::string years);
87 
88 // Plot 2D slice surfaces
89 // The macro does not call this function directly
90 // \param rootOut Already open file for root output
91 // \param prof String indicating whether surfaces are profiled or not (optional)
93  TDirectory* rootOut, std::string years, std::string prof = "");
94 
95 // Plot spectra and stack canvases -- the workhorse version
96 // The macro does not call this function directly
97 // \param spectra Array of spectra, should match the output of the SPECARR macro
98 // \param rootOut Already open file for root output
99 // \param textOFS Already open filestream for text output
100 // \param name A name for the output, only part of the saved output object's name
101 // \param title Title for the histograms. Typically this is the x axis label
102 // \param det The detector the spectra correspond to, should be 2 chars long for nice formatting
103 // \param POT INTEGER number of POT, in units of e20, i.e., a value of 6 means 6e20
104 // \param sterile Pointer to a spectrum with 4 flavor osc prediction, includes all event types (optional)
105 // \param sterileNC Pointer to a spectrum with 4 flavor osc prediction, only includes NC events (optional)
106 void PlotStack(Spectrum spectra[], TDirectory* rootOut, FILE* textOFS,
108  Spectrum* sterile = nullptr, Spectrum* sterileNC = nullptr);
109 
110 // Plot spectra and stack canvases for ND
111 // Sets up call to the workhorse version of PlotStack
112 // \param decomp The ND decomposition
113 // \param rootOut Already open file for root output
114 // \param textOFS Already open filestream for text output
115 // \param name A name for the output, only part of the saved output object's name
116 // \param title Title for the histograms. Typically this is the x axis label
117 // \param POT INTEGER number of POT, in units of e20, i.e., a value of 6 means 6e20
118 void PlotStack(IDecomp* decomp, TDirectory* rootOut, FILE* textOFS,
120 
121 // Plot spectra and stack canvases for FD
122 // Sets up call to the workhorse version of PlotStack
123 // \param pred The FD prediction object
124 // \param sCos The FD cosmic prediction spectrum
125 // \param calc Sterile calculator assuming 3 flavor oscillations
126 // \param calcSt Sterile calculator assuming 4 flavor oscillations
127 // \param rootOut Already open file for root output
128 // \param textOFS Already open filestream for text output
129 // \param name A name for the output, only part of the saved output object's name
130 // \param title Title for the histograms. Typically this is the x axis label
131 // \param POT INTEGER number of POT, in units of e20, i.e., a value of 6 means 6e20
133  Spectrum& sCos,
135  osc::IOscCalc* calcSt,
136  TDirectory* rootOut, FILE* textOFS,
138 
139 // Function to plot purity and efficiency at both FD and ND
140 // \param pFDNum Prediction object for FD, used as nominator for FD purity and efficiency
141 // \param pFDPurDen Prediction object for FD, used as denominator for FD purity
142 // \param pFDPurDen Prediction object for FD, used as denominator for FD efficiency
143 // \param sNDNum Spectrum object for ND, used as nominator for ND purity and efficiency
144 // \param sNDPurDen Spectrum object for ND, used as denominator for ND purity
145 // \param sNDPurDen Spectrum object for ND, used as denominator for ND efficiency
146 // \param calc Sterile calculator assuming 3 flavor oscillations
147 // \param rootOut Already open file for root output
148 // \param textOFS Already open filestream for text output
149 // \param name A name for the output, only part of the saved output object's name
150 // \param title Title for the histograms. Typically this is the x axis label
151 void PlotPurEff(IPrediction* pFDNum, IPrediction* pFDPurDen, IPrediction* pFDEffDen,
152  Spectrum& sNDNum, Spectrum& sNDPurDen, Spectrum& sNDEffDen,
153  osc::IOscCalc* calc, Spectrum& sCos,
154  TDirectory* rootOut, FILE* textOFS,
156 
157 // After running a fit, reset the values in a calculator to preset defaults
159 
161 {
162  TH1::AddDirectory(0);
163 
164  // Set up path to file with filled CAFAna objects and for output of analysis
165  //std::string folder = "/nova/ana/steriles/Ana01/";
166  std::string folder = "";
167  std::string filenm = "BlessedPlots14db";
168 
169  std::string loadLocation = folder + filenm + ".root";
170  std::string saveLocation = folder + filenm + "Ana.root";
171  std::string textLocation = folder + filenm + "Ana.txt";
172 
173  // Load filled objects from file
174  TFile* rootL = new TFile(loadLocation.c_str(), "UPDATE");
175  TDirectory* tmpL = gDirectory;
176  TDirectory* loadDir = gDirectory;
177 
178  loadDir->cd((loadLocation + ":/decompNC").c_str());
179  ProportionalDecomp decompNC = *ProportionalDecomp::LoadFrom(gDirectory);
180  loadDir->cd((loadLocation + ":/prediction").c_str());
182  loadDir->cd((loadLocation + ":/sData1Yr").c_str());
183  Spectrum sFDData1y = *Spectrum ::LoadFrom(gDirectory);
184  loadDir->cd((loadLocation + ":/sData3Yr").c_str());
185  Spectrum sFDData3y = *Spectrum ::LoadFrom(gDirectory);
186  loadDir->cd((loadLocation + ":/sCos1Yr").c_str());
187  Spectrum sCosmic1y = *Spectrum ::LoadFrom(gDirectory);
188  loadDir->cd((loadLocation + ":/sCos3Yr").c_str());
189  Spectrum sCosmic3y = *Spectrum ::LoadFrom(gDirectory);
190  loadDir->cd((loadLocation + ":/pFDPENum").c_str());
191  PredictionExtrap pFDPurEffNum = *PredictionExtrap ::LoadFrom(gDirectory);
192  loadDir->cd((loadLocation + ":/pFDPDen").c_str());
193  PredictionExtrap pFDPurityDen = *PredictionExtrap ::LoadFrom(gDirectory);
194  loadDir->cd((loadLocation + ":/pFDEDen").c_str());
195  PredictionExtrap pFDEfcncyDen = *PredictionExtrap ::LoadFrom(gDirectory);
196  loadDir->cd((loadLocation + ":/sNDPENum").c_str());
197  Spectrum sNDPurEffNum = *Spectrum ::LoadFrom(gDirectory);
198  loadDir->cd((loadLocation + ":/sNDPDen").c_str());
199  Spectrum sNDPurityDen = *Spectrum ::LoadFrom(gDirectory);
200  loadDir->cd((loadLocation + ":/sNDEDen").c_str());
201  Spectrum sNDEfcncyDen = *Spectrum ::LoadFrom(gDirectory);
202  tmpL->cd();
203  rootL->Close(); // Don't forget to close the file!
204 
205  // Set up oscillation calculator that uses default 3 flavor parameters
207 
208  // Set up oscillation calculator that uses default 3 flavor parameters
210  calc4f->SetNFlavors(4);
211  calc4f->SetDm(4, 0.5); //Dm41 = 0.5 eV^2
212  calc4f->SetAngle(2, 4, 0.);
213  calc4f->SetAngle(3, 4, 0.);
214 
215  std::string labelEReco = "Calorimetric Energy (GeV)";
216 
217  // Open files for saving analysis output
218  TFile* rootF = new TFile(saveLocation.c_str(), "RECREATE");
219  FILE* textF;
220  textF = fopen(textLocation.c_str(), "w");
221 
222  // Choose some angle values for making 4 flavor predictions
223  calc4f->SetAngle(2, 4, 0.17453293); // 10 degrees
224  calc4f->SetAngle(3, 4, 0.61086524); // 35 degrees
225 
226  // Plot spectra for 1 data yr
227  std::string titleHelper = " Spectra for 6e20 POT";
228  PlotStack(&pred, sCosmic1y, calc3f, calc4f, rootF, textF, "FDSpectra1yr", titleHelper, 6);
229  PlotStack(&decompNC, rootF, textF, "NDSpectra1yr", titleHelper, 6);
230 
231  // Plot spectra for 3 data yr
232  titleHelper = " Spectra for 18e20 POT";
233  PlotStack(&pred, sCosmic3y, calc3f, calc4f, rootF, textF, "FDSpectra3yr", titleHelper, 18);
234  PlotStack(&decompNC, rootF, textF, "NDSpectra3yr", titleHelper, 18);
235 
236  // Plot purity and efficiency
237  PlotPurEff(&pFDPurEffNum, &pFDPurityDen, &pFDEfcncyDen,
238  sNDPurEffNum, sNDPurityDen, sNDEfcncyDen,
239  calc3f, sCosmic1y, rootF, textF, "PurEff", labelEReco);
240 
241  fclose(textF); // Text file is done now, close it
242  ResetAngles(calc4f);
243 
244  // Create Experiment objects, one each for the 1 data yr and 3 data yr predictions
245  GaussianConstraint th23Constraint1y(&kFitTheta23InDegreesSterile, 45., 8.73);
246  CountingExperiment expt1y(&pred, sFDData1y, sCosmic1y);
247  MultiExperiment multi1y({&th23Constraint1y, &expt1y});
248 
249  GaussianConstraint th23Constraint3y(&kFitTheta23InDegreesSterile, 45., 8.73);
250  CountingExperiment expt3y(&pred, sFDData3y, sCosmic3y);
251  MultiExperiment multi3y({&th23Constraint3y, &expt3y});
252 
253  // Vector of parameters to fit: theta 23, 34, 24
254  std::vector<const IFitVar*> fitvars;
255  fitvars.push_back(&kFitTheta34InDegreesSterile);
256  fitvars.push_back(&kFitTheta24InDegreesSterile);
257  fitvars.push_back(&kFitTheta23InDegreesSterile);
258 
259  // Create fit objects
260  MinuitFitter fit1y(&multi1y, fitvars);
261  MinuitFitter fit3y(&multi3y, fitvars);
262 
263  // Run fits!
264  fit1y.Fit(calc4f);
265  ResetAngles(calc4f);
266  fit3y.Fit(calc4f);
267  ResetAngles(calc4f);
268 
269  // Set up values for slice and surface plots
270  AngleValues avals;
271  avals.nbins23 = 200;
272  avals.min23 = 20.;
273  avals.max23 = 70.;
274  avals.nbins34 = 200;
275  avals.min34 = 0.;
276  avals.max34 = 50.;
277  avals.nbins24 = 180;
278  avals.min24 = 0.;
279  avals.max24 = 45.;
280 
281  const Color_t kFitColor = kRed;
282 
283  // Plot 1D fits for 1 data yr experiment
284  Plot1DSlices(
285  &multi1y, calc4f,
289  rootF, avals, "1", kFitColor
290  );
291 
292  // Plot 1D fits for 3 data yr experiment
293  Plot1DSlices(
294  &multi3y, calc4f,
298  rootF, avals, "3", kFitColor
299  );
300 
301  // Lower the number of bins for surfaces
302  avals.nbins23 = 50;
303  avals.nbins34 = 50;
304  avals.nbins24 = 45;
305 
306  // Plot the surfaces for 1 data yr experiment
307  Plot2DSlices(
308  &multi1y, calc4f,
312  rootF, avals, "1"
313  );
314 
315  // Plot the surfaces for 3 data yr experiment
316  Plot2DSlices(
317  &multi3y, calc4f,
321  rootF, avals, "3"
322  );
323 
324  rootF->Close(); // Close the output root file
325 }
326 
327 //........................................................................
329  const IFitVar* var, std::vector<const IFitVar*> profVars,
330  TDirectory* rootOut, int nbins, double min, double max,
331  std::string name, Color_t color)
332 {
333  // Create slice histogram with/without profiling based on input profVars vector size
334 
335  TH1* h = (profVars.size() > 0 ?
336  Profile(expt, calc, var, nbins, min, max, -1, profVars) :
337  Slice(expt, calc, var, nbins, min, max));
338 
339  ResetAngles(calc); // Don't forget to reset calculator...
340 
341  // Formatting:
342  h->SetLineColor(color); // Set the line color
343  h->SetMaximum(5); // Want the plot to go from 0 to 5
344  h->SetMinimum(0);
345  h->SetName(name.c_str()); // Set the object name for saving
346  rootOut->WriteTObject(h); // Save the raw histogram
347 
348  // Create an object for writing LaTeX to canvas
349  TLatex *tex = new TLatex();
350  tex->SetNDC();
351  tex->SetTextFont(42);
352  tex->SetTextSize(0.037);
353  tex->SetTextAlign(11);
354 
355  // Coordinates for LaTeX object
356  double xTex = 0.92, y68 = 0.25, y90 = 0.52;
357 
358  // Dotted lines representing 1 and 2 sigma
359  TLine* line68 = new TLine(min, 1, max, 1);
360  TLine* line90 = new TLine(min, 2.71, max, 2.71);
361 
362  // Sometimes histogram has horizontal artifact at max canvas value
363  // This line hides it
364  TLine* line5 = new TLine(min, 5, max, 5);
365 
366  // Dashed line for 1 and 2 sigma lines
367  line68->SetLineStyle(2);
368  line90->SetLineStyle(2);
369 
370  // Formatting
371  line68->SetLineWidth(2);
372  line90->SetLineWidth(2);
373  line5 ->SetLineWidth(2);
374 
375  // Input "name" should have "h" for histogram at the beginning
376  // Make canvas name with input "name" but remove that "h"
377  std::string cname = "c1DSlice" + name.substr(1, name.size()-1);
378 
379  // Theta subscripts should be at position 6, 7, get them
380  std::string indices = name.substr(6, 2);
381 
382  // Create POT string based on whether "name" has 1 or 3 in it
383  std::string cpotstr = ( name.find("1") != std::string::npos ?
384  "6e20 POT" : "18e20 POT" );
385 
386  // Create a title for canvas
387  std::string ctitle = "Delta chi^2 for Theta" + indices;
388  // Add info on profiling if applicable
389  ctitle += ( name.find("Prof") != std::string::npos ?
390  " (Profiling Over Other Angles) for " + cpotstr :
391  " for " + cpotstr );
392 
393  // Create the canvas and draw!
394  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 600, 500);
395  h->Draw("]["); // Don't show vertical drop to 0 at edge of histogram
396  line68->DrawClone("same");
397  line90->DrawClone("same");
398  line5 ->DrawClone("same");
399  tex->DrawLatex(xTex, y68, "68%"); // Add label for 1 sigma line
400  tex->DrawLatex(xTex, y90, "90%"); // Add label for 2 sigma line
401  Simulation();
402  rootOut->WriteTObject(c); // Save full canvas
403 
404  return;
405 }
406 
407 //........................................................................
409  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
410  TDirectory* rootOut, AngleValues a,
411  std::string years, Color_t color)
412 {
413  // General form of the object name for saving
414  std::string name = "h" + years + "yrTh";
415  std::string fullname; // fullname adds 2 characters, the theta subscript indices
416 
417  // Plot unprofiled slices, for theta 23, then 34, then 24
418  fullname = name + "23";
419  Plot1DSlice(expt, calc, var23, {},
420  rootOut, a.nbins23, a.min23, a.max23, fullname, color);
421  fullname = name + "34";
422  Plot1DSlice(expt, calc, var34, {},
423  rootOut, a.nbins34, a.min34, a.max34, fullname, color);
424  fullname = name + "24";
425  Plot1DSlice(expt, calc, var24, {},
426  rootOut, a.nbins24, a.min24, a.max24, fullname, color);
427 
428 
429  // Plot profiled slices in same order as above
430  // When slice is of one angle, the other two are profiled
431  // I.e., when the slice is theta 23, theta 34 and 24 are profiled
432  fullname = name + "23Prof";
433  Plot1DSlice(expt, calc, var23, {var34, var24},
434  rootOut, a.nbins23, a.min23, a.max23, fullname, color);
435  fullname = name + "34Prof";
436  Plot1DSlice(expt, calc, var34, {var23, var24},
437  rootOut, a.nbins34, a.min34, a.max34, fullname, color);
438  fullname = name + "24Prof";
439  Plot1DSlice(expt, calc, var24, {var23, var34},
440  rootOut, a.nbins24, a.min24, a.max24, fullname, color);
441 
442  return;
443 }
444 
445 //........................................................................
447  const IFitVar* var23, const IFitVar* var34, const IFitVar* var24,
448  TDirectory* rootOut, AngleValues a, std::string years)
449 {
450  // Create unprofiled surfaces
451  // 34 vs 24
452  FrequentistSurface s3424(expt, calc,
453  var34, a.nbins34, a.min34, a.max34,
454  var24, a.nbins24, a.min24, a.max24);
455  ResetAngles(calc);
456 
457  // 23 vs 24
458  FrequentistSurface s2324(expt, calc,
459  var23, a.nbins23, a.min23, a.max23,
460  var24, a.nbins24, a.min24, a.max24);
461  ResetAngles(calc);
462 
463  // 34 vs 23
464  FrequentistSurface s3423(expt, calc,
465  var34, a.nbins34, a.min34, a.max34,
466  var23, a.nbins23, a.min23, a.max23);
467  ResetAngles(calc);
468 
469  // Create profiled surfaces in same order as above
470  // The third angle is the one that is profiled
471  // I.e., the 34 vs 24 surface profiles 23
472  FrequentistSurface s3424Prof(expt, calc,
473  var34, a.nbins34, a.min34, a.max34,
474  var24, a.nbins24, a.min24, a.max24,
475  {var23});
476  ResetAngles(calc);
477 
478  FrequentistSurface s2324Prof(expt, calc,
479  var23, a.nbins23, a.min23, a.max23,
480  var24, a.nbins24, a.min24, a.max24,
481  {var34});
482  ResetAngles(calc);
483 
484  FrequentistSurface s3423Prof(expt, calc,
485  var34, a.nbins34, a.min34, a.max34,
486  var23, a.nbins23, a.min23, a.max23,
487  {var24});
488  ResetAngles(calc);
489 
490  // Plot all of the surfaces
491  Plot2DSlices(&s3424, &s2324, &s3423, rootOut, years);
492  Plot2DSlices(&s3424Prof, &s2324Prof, &s3423Prof, rootOut, years, "Prof");
493 
494  return;
495 }
496 
497 //........................................................................
499  TDirectory* rootOut, std::string years, std::string prof)
500 {
501  // General format for the name for saving
502  std::string name = "h" + years + "yrTh";
503  std::string fullname;
504 
505  // Add indices for first theta, add second theta, and profile string if applicaable
506  fullname = name + "34Th24" + prof;
507  TH1* h3424 = s3424->ToTH2(); // Create a histogram from the surface
508  h3424->SetName(fullname.c_str()); // Set the name for saving
509  rootOut->WriteTObject(h3424); // Save the raw histogram
510 
511  fullname = name + "23Th24" + prof;
512  TH1* h2324 = s2324->ToTH2();
513  h2324->SetName(fullname.c_str());
514  rootOut->WriteTObject(h2324);
515 
516  fullname = name + "34Th23" + prof;
517  TH1* h3423 = s3423->ToTH2();
518  h3423->SetName(fullname.c_str());
519  rootOut->WriteTObject(h3423);
520 
521  // Create name for canvas
522  std::string cname = "c2DSlice" + years + "yr" + prof;
523  // Create string for POT based on whether name has 1 or 3 in it
524  std::string potstr = ( name.find("1") != std::string::npos ?
525  "6e20 POT" : "18e20 POT" );
526  // Create relevant title for the canvas
527  std::string ctitle = (
528  prof.size() > 0 ?
529  "Surfaces (Profiling Over Other Angles) for " + potstr :
530  "Surfaces for " + potstr
531  );
532 
533  // Create triptych canvas and draw!
534  TCanvas* c = new TCanvas(cname.c_str(), ctitle.c_str(), 800, 800);
535  c->Divide(2, 2); // Split into 2 x 2 grid
536  c->cd(1); // Go to top left
537  s3424->Draw(); // Draw histogram
538  s3424->DrawContour(Gaussian68Percent2D(*s3424), 7, kBlack); // Draw 1 sigma contour
539  s3424->DrawContour(Gaussian90Percent2D(*s3424), 1, kBlack); // Draw 2 sigma contour
540  s3424->DrawBestFit(kBlack); // Draw best fit point
541  c->cd(2); // Go to top right
542  s2324->Draw();
543  s2324->DrawContour(Gaussian68Percent2D(*s2324), 7, kBlack);
544  s2324->DrawContour(Gaussian90Percent2D(*s2324), 1, kBlack);
545  s2324->DrawBestFit(kBlack);
546  Simulation();
547  c->cd(3); // Go to bottom left
548  s3423->Draw();
549  s3423->DrawContour(Gaussian68Percent2D(*s3423), 7, kBlack);
550  s3423->DrawContour(Gaussian90Percent2D(*s3423), 1, kBlack);
551  s3423->DrawBestFit(kBlack);
552  rootOut->WriteTObject(c); // Save the root canvas
553 
554  return;
555 }
556 
557 //........................................................................
558 void PlotStack(Spectrum spectra[], TDirectory* rootOut, FILE* textOFS,
560  Spectrum* sterile, Spectrum* sterileNC)
561 {
562  // The input is in units of 1e20, turn this into the real scale
563  double realPOT = (double)POT*1e20;
564 
565  char potBuff[2];
566  sprintf(potBuff, "%d", POT);
567  std::string potStr = std::string(potBuff) + " x 10^{20} POT";
568 
569  // Get each histogram from the Spectrum array
570  TH1* hAll = spectra[0].ToTH1(realPOT);
571  TH1* hNC = spectra[1].ToTH1(realPOT);
572  TH1* hNumu = spectra[2].ToTH1(realPOT);
573  TH1* hNue = spectra[3].ToTH1(realPOT);
574  TH1* hCos = spectra[4].ToTH1(realPOT);
575  // The default for these inputs are nullptr, only get a TH1 if they exist
576  TH1* hStrl = (sterile ? sterile ->ToTH1(realPOT) : nullptr);
577  TH1* hStrlNC = (sterileNC ? sterileNC->ToTH1(realPOT) : nullptr);
578 
579  // Create the title and axis labels, which is an empty title,
580  // the x label from an input Spectrum, and Events/POT
581  std::string xLabel(hAll->GetXaxis()->GetTitle());
582  std::string yLabel = (det.compare("ND") == 0 ? "10^{3} Events" : "Events");
583  yLabel += "/" + potStr + " / 0.25 GeV";
584  std::string histTitle = ";" + xLabel + ";" + yLabel;
585 
586  // Print the column title and headers
587  fprintf(textOFS, "Event numbers at the %2.2s for %2de20 POT:\n", det.c_str(), POT);
588  if(hStrlNC) { fprintf(textOFS, "%12.12s, ", "NC (Sterile)"); }
589  else { fprintf(textOFS, "%14.14s", ""); }
590  fprintf(textOFS, "%12.12s, %12.12s, %12.12s, %12.12s, %12.12s, %12.12s\n",
591  "NC (3 Flav)", "All", "Numu", "Nue", "Cosmic", "FOM");
592 
593  // Get the integral number of events, then print them
594  double nNC3F = hNC ->Integral();
595  double nAll = hAll ->Integral();
596  double nNumu = hNumu->Integral();
597  double nNue = hNue ->Integral();
598  double nCos = hCos ->Integral();
599  double fom = nNC3F/sqrt(nAll);
600  if(hStrlNC) {
601  double nNCSt = hStrlNC->Integral();
602  fprintf(textOFS, "%10E, ", nNCSt);
603  }
604  else {
605  fprintf(textOFS, "%14.14s", "");
606  }
607  fprintf(textOFS, "%10E, %10E, %10E, %10E, %10E, %10E\n\n",
608  nNC3F, nAll, nNumu, nNue, nCos, fom);
609 
610  // Stylize the unstacked plots
611  CenterTitles(hNue);
612  CenterTitles(hNumu);
613  CenterTitles(hNC);
614  CenterTitles(hCos);
615  CenterTitles(hAll);
616 
617  hAll ->SetMinimum(0);
618  hNC ->SetMinimum(0);
619  hNumu->SetMinimum(0);
620  hNue ->SetMinimum(0);
621  hCos ->SetMinimum(0);
622  if(hStrlNC) { hStrl->SetMinimum(0); }
623 
624  const Color_t kCosmicBackgroundColor = kOrange+1;
625  hAll ->SetLineColor(kTotalMCColor);
626  hNC ->SetLineColor(kNCBackgroundColor);
627  hNumu->SetLineColor(kNumuBackgroundColor);
628  hNue ->SetLineColor(kNueSignalColor);
629  hCos ->SetLineColor(kCosmicBackgroundColor);
630  if(hStrlNC) {
631  hStrlNC->SetLineColor(kNCBackgroundColor);
632  hStrlNC->SetLineStyle(2);
633  }
634 
635  hAll ->SetTitle(histTitle.c_str());
636  hNC ->SetTitle(histTitle.c_str());
637  hNumu->SetTitle(histTitle.c_str());
638  hNue ->SetTitle(histTitle.c_str());
639  hCos ->SetTitle(histTitle.c_str());
640  if(hStrlNC) { hStrlNC->SetTitle(histTitle.c_str()); }
641 
642  // Create "stacked" copies of flavored histograms
643  // Cos Stack = Cos
644  TH1* hCosStack = (TH1*)hCos->Clone();
645 
646  // Numu Stack = Cos Stack + Numu = Cos + Numu
647  TH1* hNumuStack = (TH1*)hCosStack->Clone();
648  hNumuStack->Add(hNumu);
649 
650  // Nue Stack = Numu Stack + Nue = Cos + Numu + Nue
651  TH1* hNueStack = (TH1*)hNumuStack->Clone();
652  hNueStack->Add(hNue);
653 
654  // NC Stack = Nue Stack + NC = Cos + Numu + Nue + NC
655  TH1* hNCStack = (TH1*)hNC->Clone();
656  //hNCStack->Add(hNC);
657 
658  // Stylize the stacked plots
659  CenterTitles(hNueStack);
660  CenterTitles(hNumuStack);
661  CenterTitles(hNCStack);
662  CenterTitles(hCosStack);
663  if(hStrlNC) { CenterTitles(hStrlNC); }
664 
665  hNCStack ->SetLineColor(kNCBackgroundColor);
666  hNumuStack->SetLineColor(kNumuBackgroundColor);
667  hNumuStack->SetFillColor(kNumuBackgroundColor);
668  hNueStack ->SetLineColor(kNueSignalColor);
669  hNueStack ->SetFillColor(kNueSignalColor);
670  hCosStack ->SetLineColor(kCosmicBackgroundColor);
671  hCosStack ->SetFillColor(kCosmicBackgroundColor);
672 
673  hNCStack ->SetMinimum(0);
674  hNumuStack->SetMinimum(0);
675  hNueStack ->SetMinimum(0);
676  hCosStack ->SetMinimum(0);
677 
678  // Manually set errors, and calculate the max value (including errors)
679  double maxval = 0.;
680  double maxvalstack = 0.;
681  for(int i = 1, n = hNC->GetNbinsX(); i <= n; ++i) {
682  double error = sqrt(hNC->GetBinContent(i));
683  hNC ->SetBinError(i, error);
684  hNCStack->SetBinError(i, error);
685  maxval = std::max(maxval, std::max(hAll->GetBinContent(i), hNC->GetBinContent(i) + error));
686  maxvalstack = std::max(maxvalstack, hNCStack->GetBinContent(i) + error);
687  }
688 
689  // Scale down the ND by 10e3 (don't forget the max values!)
690  if(det.compare("ND") == 0) {
691  double ndScale = 1./1000.;
692  hAll ->Scale(ndScale);
693  hNC ->Scale(ndScale);
694  hNue ->Scale(ndScale);
695  hNumu->Scale(ndScale);
696  hNCStack ->Scale(ndScale);
697  hNueStack ->Scale(ndScale);
698  hNumuStack->Scale(ndScale);
699  maxval *= ndScale;
700  maxvalstack *= ndScale;
701  }
702 
703  TLatex *tex = new TLatex();
704  tex->SetNDC();
705  tex->SetTextFont(42);
706  tex->SetTextSize(0.037);
707  tex->SetTextAlign(11);
708 
709  // Create a canvas for and draw the unstacked histograms
710  TCanvas* c = new TCanvas(name.c_str(), title.c_str(), 800, 500);
711  hNue->SetMaximum(maxval*1.1);
712  hNue->GetXaxis()->SetNdivisions(-406);
713 
714  // Create a legend for the unstacked histograms
715  double xL = 0.5, xR = 0.7, yB = 0.59, yT = 0.84;
716  if(det.compare("ND") == 0) { yB += 0.05; }
717  TLegend* leg = new TLegend(xL, yB, xR, yT);
718  leg->SetBorderSize(0);
719  leg->SetFillColor(kWhite);
720  leg->SetFillStyle(0);
721  leg->SetFillStyle(0);
722  leg->SetTextFont(42);
723  leg->SetTextSize(0.037);
724  leg->AddEntry(hAll, "Total Prediction", "l");
725  leg->AddEntry(hNC, "NC 3 Flavor Prediction", "le");
726  leg->AddEntry(hNue, "#nu_{e} CC Background", "l");
727  leg->AddEntry(hNumu, "#nu_{#mu} CC Background", "l");
728  if(det.compare("FD") == 0) {
729  leg->AddEntry(hCos, "Cosmic Background", "l");
730  }
731 
732  hNue ->Draw("hist");
733  hNumu->Draw("hist same");
734  if(det.compare("FD") == 0) { hCos->Draw("hist same"); }
735  hAll ->Draw("hist same");
736  hNC ->Draw("same");
737  leg->Draw();
738  if(det.compare("FD") == 0) {
739  tex->DrawLatex(0.505, yB - 0.045, "#Deltam^{2}_{41} = 0.5 eV^{2}, #theta_{24} = 10#circ, #theta_{34} = 35#circ");
740  tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.37x10^{-3} eV^{2}, #theta_{13} = 8.5#circ, #theta_{23} = 45#circ");
741  tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
742  }
743 
744  gPad->RedrawAxis();
745  Simulation();
746 
747  // Save the canvas
748  rootOut->WriteTObject(c);
749 
750  // Create a canvas for and draw the stacked histograms
751  TCanvas* cStack = new TCanvas((name + "Stack").c_str(), title.c_str(), 800, 500);
752  hNueStack->SetMaximum(maxvalstack*1.1);
753  hNueStack->GetXaxis()->SetNdivisions(-406);
754 
755  if(det.compare("ND") == 0) { yB += 0.05; }
756  TLegend* legStack = new TLegend(xL, yB, xR, yT);
757  legStack->SetBorderSize(0);
758  legStack->SetFillColor(kWhite);
759  legStack->SetFillStyle(0);
760  legStack->SetFillStyle(0);
761  legStack->SetTextFont(42);
762  legStack->SetTextSize(0.037);
763  legStack->AddEntry(hNCStack, "NC 3 Flavor Prediction", "le");
764  if(hStrlNC) { legStack->AddEntry(hStrlNC, "NC 4 Flavor Prediction", "l"); }
765  legStack->AddEntry(hNueStack, "#nu_{e} CC Background", "f");
766  legStack->AddEntry(hNumuStack, "#nu_{#mu} CC Background", "f");
767  if(det.compare("FD") == 0) {
768  legStack->AddEntry(hCosStack, "Cosmic Background", "f");
769  }
770 
771  hNueStack ->Draw("hist");
772  hNumuStack->Draw("hist same");
773  if(det.compare("FD") == 0) { hCosStack->Draw("hist same"); }
774  hNCStack ->Draw("same");
775  if(hStrlNC) { hStrlNC->Draw("hist same"); }
776  legStack->Draw();
777  if(det.compare("FD") == 0) {
778  tex->DrawLatex(0.505, yB - 0.045, "#Deltam^{2}_{41} = 0.5 eV^{2}, #theta_{24} = 10#circ, #theta_{34} = 35#circ");
779  tex->DrawLatex(0.505, yB - 0.100, "#Deltam^{2}_{32} = 2.37x10^{-3} eV^{2}, #theta_{13} = 8.5#circ, #theta_{23} = 45#circ");
780  tex->DrawLatex(0.505, yB - 0.155, potStr.c_str());
781  }
782 
783  gPad->RedrawAxis();
784  Simulation();
785 
786  rootOut->WriteTObject(cStack);
787 
788  return;
789 }
790 
791 //........................................................................
792 void PlotStack(IDecomp* decomp, TDirectory* rootOut, FILE* textOFS,
794 {
795  Spectrum empty = decomp->NCTotalComponent();
796  empty.Clear();
797 
798  Spectrum spectraND[MAXSPEC] = {
799  decomp->NCTotalComponent() +
800  decomp->NueComponent() + decomp->AntiNueComponent() +
801  decomp->NumuComponent() + decomp->AntiNumuComponent(),
802  decomp->NCTotalComponent(),
803  decomp->NumuComponent() + decomp->AntiNumuComponent(),
804  decomp->NueComponent() + decomp->AntiNueComponent(),
805  empty
806  };
807 
808  std::string det = "ND";
809  std::string fulltitle = det + title;
810  PlotStack(spectraND, rootOut, textOFS, name, fulltitle, det, POT);
811 
812  return;
813 }
814 
815 //........................................................................
817  osc::IOscCalc* calc, osc::IOscCalc* calcSt,
818  TDirectory* rootOut, FILE* textOFS,
820 {
821  Spectrum spectraFD[MAXSPEC] = {
822  pred->Predict(calc) + sCos,
826  sCos
827  };
828 
829  Spectrum sterile = pred->Predict(calcSt) + sCos;
830  Spectrum sterileNC = pred->PredictComponent(calcSt, Flavors::kAll, Current::kNC, Sign::kBoth);
831 
832  std::string det = "FD";
833  std::string fulltitle = det + title;
834  PlotStack(spectraFD, rootOut, textOFS, name, fulltitle, det, POT, &sterile, &sterileNC);
835 
836  return;
837 }
838 
839 //........................................................................
840 void PlotPurEff(IPrediction* pFDNum, IPrediction* pFDPurDen, IPrediction* pFDEffDen,
841  Spectrum& sNDNum, Spectrum& sNDPurDen, Spectrum& sNDEffDen,
842  osc::IOscCalc* calc, Spectrum& sCos,
843  TDirectory* rootOut, FILE* textOFS,
845 {
846  Spectrum sFDNum = pFDNum ->Predict(calc);
847  Spectrum sFDPurDen = pFDPurDen->Predict(calc) + sCos;
848  Spectrum sFDEffDen = pFDEffDen->Predict(calc);
849 
850  TH1* hFDNum = sFDNum .ToTH1(NCSCALE);
851  TH1* hFDPurDen = sFDPurDen.ToTH1(NCSCALE);
852  TH1* hFDEffDen = sFDEffDen.ToTH1(NCSCALE);
853 
854  TH1* hNDNum = sNDNum .ToTH1(NCSCALE);
855  TH1* hNDPurDen = sNDPurDen.ToTH1(NCSCALE);
856  TH1* hNDEffDen = sNDEffDen.ToTH1(NCSCALE);
857 
858  TH1* hFDEff = (TH1*)hFDNum->Clone();
859  hFDEff->Divide(hFDEffDen);
860 
861  TH1* hNDEff = (TH1*)hNDNum->Clone();
862  hNDEff->Divide(hNDEffDen);
863 
864  TH1* hFDPur = (TH1*)hFDNum->Clone();
865  hFDPur->Divide(hFDPurDen);
866 
867  TH1* hNDPur = (TH1*)hNDNum->Clone();
868  hNDPur->Divide(hNDPurDen);
869 
870  hFDEff->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
871  hFDPur->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
872  hNDEff->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
873  hNDPur->SetTitle((";" + title + ";Efficiency, Purity / 0.25 GeV").c_str());
874 
875  hFDEff->SetLineColor(kBlue);
876  hNDEff->SetLineColor(kBlue);
877 
878  hFDPur->SetLineColor(kRed);
879  hNDPur->SetLineColor(kRed);
880 
881  hNDEff->SetLineStyle(2);
882  hNDPur->SetLineStyle(2);
883 
884  hFDEff->SetMaximum(1.);
885  hFDEff->SetMinimum(0.);
886  hFDEff->GetXaxis()->SetNdivisions(-406);
887 
888  hNDEff->SetMaximum(1.);
889  hNDEff->SetMinimum(0.);
890 
891  hFDPur->SetMaximum(1.);
892  hFDPur->SetMinimum(0.);
893 
894  hNDPur->SetMaximum(1.);
895  hNDPur->SetMinimum(0.);
896 
897  CenterTitles(hFDEff);
898  CenterTitles(hNDEff);
899  CenterTitles(hFDPur);
900  CenterTitles(hNDPur);
901 
902  TLegend* leg = new TLegend(0.6, 0.35, 0.85, 0.55);
903  leg->AddEntry(hFDEff, "FD MC Efficiency", "L");
904  leg->AddEntry(hNDEff, "ND MC Efficiency", "L");
905  leg->AddEntry(hFDPur, "FD MC Purity", "L");
906  leg->AddEntry(hNDPur, "ND MC Purity", "L");
907 
908  TCanvas* c = new TCanvas(name.c_str(), "Purity and Efficiency", 800, 500);
909  hFDEff->Draw("hist");
910  hNDEff->Draw("hist same");
911  hFDPur->Draw("hist same");
912  hNDPur->Draw("hist same");
913  leg->Draw();
914 
915  gPad->RedrawAxis();
916  Simulation();
917 
918  rootOut->WriteTObject(c);
919 }
920 
921 
922 //........................................................................
924 {
925  calc->SetAngle(2, 3, M_PI/4); // 45 degrees (maximal mixing)
926  calc->SetAngle(3, 4, 0.);
927  calc->SetAngle(2, 4, 0.);
928 
929  return;
930 }
void PlotStack(Spectrum spectra[], TDirectory *rootOut, FILE *textOFS, std::string name, std::string det, double POT, double potEquiv, PlotOptions opt, Spectrum *dataspec)
void Simulation()
Definition: tools.h:16
T max(const caf::Proxy< T > &a, T b)
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void SetNFlavors(int nflavors)
virtual Spectrum PredictComponent(osc::IOscCalc *calc, Flavors::Flavors_t flav, Current::Current_t curr, Sign::Sign_t sign) const =0
virtual Spectrum AntiNueComponent() const =0
void PlotPurEff(Spectrum pureff[], Spectrum effcyD[], TDirectory *out, std::string name, std::string title)
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:209
virtual Spectrum NumuComponent() const =0
#define MAXSPEC
Definition: NusLoadProd3.h:18
A simple Gaussian constraint on an arbitrary IFitVar.
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
T sqrt(T number)
Definition: d0nt_math.hpp:156
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
void Plot1DSlices(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, TDirectory *rootOut, AngleValues a, std::string years, Color_t color)
static std::unique_ptr< PredictionSterile > LoadFrom(TDirectory *dir, const std::string &name)
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
void Clear()
Definition: Spectrum.cxx:433
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 CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
osc::OscCalcDumb calc
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
const Color_t kNumuBackgroundColor
Definition: Style.h:30
#define M_PI
Definition: SbMath.h:34
static std::unique_ptr< ProportionalDecomp > LoadFrom(TDirectory *dir, const std::string &name)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
const Color_t kNueSignalColor
Definition: Style.h:19
void BlessedPlotsAna()
fclose(fg1)
Log-likelihood scan across two parameters.
const int nbins
Definition: cellShifts.C:15
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:32
Charged-current interactions.
Definition: IPrediction.h:39
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 Draw() const
Draw the surface itself.
Definition: ISurface.cxx:19
const double a
void Plot2DSlices(IExperiment *expt, osc::OscCalcSterile *calc, const IFitVar *var23, const IFitVar *var34, const IFitVar *var24, TDirectory *rootOut, AngleValues a, std::string years)
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
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:64
static std::unique_ptr< PredictionExtrap > LoadFrom(TDirectory *dir, const std::string &name)
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
virtual Spectrum AntiNumuComponent() const =0
std::vector< double > POT
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:161
Combine multiple component experiments.
void SetAngle(int i, int j, double th)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Splits Data proportionally according to MC.
void SetDm(int i, double dm)
void ResetAngles(osc::OscCalcSterile *calc)
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
virtual Spectrum NCTotalComponent() const
Definition: IDecomp.h:18
TLatex * tex
Definition: f2_nu.C:499
Standard interface to all decomposition techniques.
Definition: IDecomp.h:13
Base class defining interface for experiments.
Definition: IExperiment.h:14
Neutral-current interactions.
Definition: IPrediction.h:40
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
Interface definition for fittable variables.
Definition: IFitVar.h:16
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
A prediction object compatible with sterile oscillations.
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
static double NCSCALE
const Color_t kTotalMCColor
Definition: Style.h:16
Take the output of an extrapolation and oscillate it as required.
All neutrinos, any flavor.
Definition: IPrediction.h:26
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
const Color_t kNCBackgroundColor
Definition: Style.h:22
Compare a data spectrum to MC expectation using only the event count.
virtual Spectrum NueComponent() const =0
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17