MakeExtrapSurface.C
Go to the documentation of this file.
1 //////////////////////////////////////////////////////////////
2 // MakeExtrapSurface.C
3 // M Wallbank (wallbank@fnal.gov)
4 //
5 // Make 1D profiles and 2D surface in th24-th34 parameter
6 // space, using extrapolated predictions.
7 //////////////////////////////////////////////////////////////
8 
9 // framework
10 #include "CAFAna/Analysis/Calcs.h"
12 #include "CAFAna/Fit/Fit.h"
13 #include "CAFAna/Core/IFitVar.h"
23 
24 // root
25 #include "TFile.h"
26 #include "TH1.h"
27 
28 // stl
29 #include <sstream>
30 
31 using namespace ana;
32 
33 void MakeExtrapSurface(TString opt) {
34 
35  TH1::AddDirectory(0);
36 
37  // Get running options
38  bool fhc = opt.Contains("FHC", TString::kIgnoreCase);
39  bool rhc = opt.Contains("RHC", TString::kIgnoreCase);
40  assert(fhc || rhc);
41 
42  bool oned = opt.Contains("1D", TString::kIgnoreCase);
43  bool twod = opt.Contains("2D", TString::kIgnoreCase);
44  assert(oned || twod);
45 
46  bool systematics = opt.Contains("syst", TString::kIgnoreCase);
47  bool cos = opt.Contains("cos", TString::kIgnoreCase);
48  bool data = opt.Contains("data", TString::kIgnoreCase);
49 
50  bool scale = opt.Contains("scale", TString::kIgnoreCase);
51 
52  // Get the prediction
55 
56  // Get POT & livetimes
57  double pot, livetime;
58  if (scale) {
59  double scale_pot = opt.Contains("125")? 12.5 : opt.Contains("12")? 12 : opt.Contains("35")? 35 : -1;
60  assert(scale_pot > 0);
61  pot = scale_pot * 1e20;
63  }
64  else {
65  pot = rhc? kAna2018RHCPOT : kAna2018FHCPOT;
66  livetime = rhc? kAna2018RHCLivetime : kAna2018FHCLivetime;
67  }
68 
69  std::cout << "POT is " << pot << ", livetime is " << livetime << std::endl;
70 
71  // Get the systs
72  std::vector<const ISyst*> systs = {};
73  if (systematics)
74  systs = LoadISysts();
75  //systs = rhc? GetNus18RHCExtrapSysts() : GetNus18FHCFDSysts("all");
76 
77  // Set up oscillation calculator that uses default 3 flavor parameters
79  calc->SetNFlavors(4);
80  SetAngles(calc);
81 
82  // Get the far detector data
83  // -- make sure to scale appropriately
84  Spectrum prediction = pred->Predict(calc);
85  cosmic.Scale(livetime/cosmic.Livetime());
86  cosmic.OverridePOT(pot);
87  Spectrum cos_spec = cosmic.FakeData(pot);
88  Spectrum fddata_spec = data ? LoadData(0, rhc) : (prediction.FakeData(pot) + cos_spec);
89 
90  // Create experiment objects
91  // Gaussian constraint on th23, from NOvA Jan 2018 results (joint numu+nue)
92  GaussianConstraint th23Constraint(&kFitSinSqTheta23Sterile, 0.558, 0.042);
93 
94  // Make experiment
95  SingleSampleExperiment experiment_ss(pred, fddata_spec, cos_spec);
96  MultiExperiment experiment({&th23Constraint, &experiment_ss});
97 
98  // Variables to profile over
99  std::vector<const IFitVar*> prof_vars = {&kFitSinSqTheta23Sterile, &kFitDelta24InPiUnitsSterile};
100  std::vector<const IFitVar*> prof_vars_th24 = prof_vars, prof_vars_th34 = prof_vars;
101  prof_vars_th24.push_back(&kFitTheta34InDegreesSterile);
102  prof_vars_th34.push_back(&kFitTheta24InDegreesSterile);
103 
104  // Binnings
105  int nbins34_1d = 240, nbins34_2d = 60;
106  float min34 = 0.;
107  float max34 = 60.;
108  float nbins24_1d = 20, nbins24_2d = 40; //160
109  float min24 = 0.;
110  float max24 = 40.;
111 
112  // 1D profiles
113  TH1* ProfileTh24, *ProfileTh34;
114  if (oned) {
115  SetAngles(calc);
116  ProfileTh24 = Profile(&experiment, calc,
117  &kFitTheta24InDegreesSterile, nbins24_1d, min24, max24,
118  -1, prof_vars_th24, systs);
119  SetAngles(calc);
120  ProfileTh34 = Profile(&experiment, calc,
121  &kFitTheta34InDegreesSterile, nbins34_1d, min34, max34,
122  -1, prof_vars_th34, systs);
123  std::cout << "Profiles made! Mean is " << ProfileTh24->GetMean() << std::endl;
124  std::cout << "Profiles made! NEntries is " << ProfileTh24->GetEntries() << std::endl;
125  }
126 
127  // 2D surface
128  FrequentistSurface* SurfaceTh34Th24;
129  if (twod) {
130  SetAngles(calc);
131  std::map<const IFitVar*, std::vector<double> > seedPoints;
132  //seedPoints[&kFitDelta24InPiUnitsSterile] = {TMath::Pi()/2};
133  SurfaceTh34Th24 = new FrequentistSurface(&experiment, calc,
134  &kFitTheta34InDegreesSterile, nbins34_2d, min34, max34,
135  &kFitTheta24InDegreesSterile, nbins24_2d, min24, max24,
136  prof_vars, systs, seedPoints);
137  }
138 
139  // Open files for saving analysis output
140  TFile* outFile = new TFile("ExtrapSurface.root", "RECREATE");
141  if (oned) {
142  ProfileTh24->Write("ProfileTh24");
143  ProfileTh34->Write("ProfileTh34");
144  }
145  if (twod)
146  SurfaceTh34Th24->SaveTo(outFile->mkdir("SurfaceTh34Th24"));
147  outFile->Close();
148  delete outFile;
149 
150 }
TGraph * Profile(const IChiSqExperiment *expt, osc::IOscCalculatorAdjustable *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
IPrediction * LoadExtrapPrediction(TString opt)
Spectrum LoadCosmicSpectrum(bool rhc=false)
Function to return cosmic spectrum.
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
std::vector< SystGroupDef > systs
Definition: syst_header.h:384
Adapt the PMNS_Sterile calculator to standard interface.
void MakeExtrapSurface(TString opt)
A simple Gaussian constraint on an arbitrary IFitVar.
void OverridePOT(double newpot)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:269
void SaveTo(TDirectory *dir) const
std::vector< const ISyst * > LoadISysts()
Function to load ISysts.
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
void Scale(double c)
Multiply this spectrum by a constant c.
Definition: Spectrum.cxx:734
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:37
const double kAna2018RHCPOT
Definition: Exposures.h:208
const XML_Char const XML_Char * data
Definition: expat.h:268
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
Log-likelihood scan across two parameters.
Double_t scale
Definition: plot.C:25
Sum up livetimes from individual cosmic triggers.
TFile * outFile
Definition: PlotXSec.C:135
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:819
Spectrum LoadData(bool nd, bool rhc=false)
Function to load NuS data file.
void SetAngles(osc::OscCalculatorSterile *calc)
#define pot
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
const FitSinSqTheta23Sterile kFitSinSqTheta23Sterile
T cos(T number)
Definition: d0nt_math.hpp:78
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
bool systematics
Definition: fnexvscaf.C:31
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const double kAna2018FHCPOT
Definition: Exposures.h:207
osc::OscCalculatorSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:266
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const double kAna2018RHCLivetime
Definition: Exposures.h:214
Compare a single data spectrum to the MC + cosmics expectation.