MakeSurface.C
Go to the documentation of this file.
1 /*
2  * MakeSurface.C
3  *
4  * Created on: May 7, 2018
5  *
6  * Author: J. Hewes <jhewes15@fnal.gov>
7  *
8  * Create contours for 2019 covmx analysis
9  */
10 
12 #include "CAFAna/Core/Progress.h"
13 #include "CAFAna/Core/Sample.h"
14 
16 
19 
20 #include "CAFAna/Analysis/CovMxSurface.h"
23 #include "CAFAna/Analysis/Fit.h"
24 
25 #include "NuXAna/nus/Nus19FHC/Utilities.h"
26 
27 #include <ctime>
28 
29 using namespace ana;
30 
31 void MakeSurface() {
32 
33  DontAddDirectory guard;
34 
35  /////////////////////////////////////////////////////////////////////
36  // //
37  // PREPARE INPUTS //
38  // //
39  /////////////////////////////////////////////////////////////////////
40 
41  // Oscillation calculator
42  osc::OscCalculatorSterile* calc = DefaultSterileCalc(4);
43  SetAngles(calc);
44 
45  // Load predictions
46  std::vector<covmx::Sample> samples = { kNusFHCNearDet, kNusFHCFarDet };
47  std::vector<IPrediction*> preds;
48  for (auto sample : samples) preds.push_back(LoadPrediction(sample, "nosysts"));
49 
50  // Load cosmics
51  std::vector<Spectrum*> cosmic;
52  for (auto sample : samples)
53  cosmic.push_back(LoadCosmic(sample).release());
54 
55  // Load data
56  std::vector<Spectrum*> data;
57  for (size_t i = 0; i < preds.size(); ++i) {
58  data.push_back(new Spectrum(preds[i]->Predict(calc).FakeData(12e20)));
59  if(samples[i].detector == covmx::Detector::kFarDet)
60  data.back()->OverrideLivetime(kNus19SensLivetime);
61  if (cosmic[i]) *data[i] += *cosmic[i];
62  }
63 
64  // Create an empty matrix and then add in the ones you produced
65  CovarianceMatrix* mx = nullptr;
66  std::vector<Binning> bins;
67  for (auto d : data) bins.push_back(d->GetBinnings()[0]);
68  mx = new CovarianceMatrix(bins);
70  TFile* mxfile = TFile::Open("covmx.root", "read");
71  CovarianceMatrix* mx_cherenkov = CovarianceMatrix::LoadFrom(mxfile->GetDirectory("cherenkov")).release();
72  CovarianceMatrix* mx_abs_calib = CovarianceMatrix::LoadFrom(mxfile->GetDirectory("abs_calib")).release();
73  mx->AddMatrix(mx_cherenkov);
74  mx->AddMatrix(mx_abs_calib);
75 
76  /////////////////////////////////////////////////////////////////////
77  // //
78  // MAKE SURFACE //
79  // //
80  /////////////////////////////////////////////////////////////////////
81 
82  // Create experiment
83  OscCovMxExperiment* expt = new OscCovMxExperiment(preds, mx, data, cosmic);
84 
85  // Set Gaussian constraints, make MultiExperiment
86  std::vector<const IExperiment*> expts = GetConstraints();
87  expts.push_back(expt);
88  MultiExperiment multi_expt(expts);
89 
90  // Binning, fit vars and seeds
91  std::vector<const IFitVar*> fit_vars = { &kFitTheta23Sterile, &kFitDmSq32Sterile,
93  std::map<const IFitVar*, std::vector<double> > seed_values;
94  seed_values[&kFitTheta23Sterile] = { kNus19Th23 };
95  seed_values[&kFitDmSq32Sterile] = { kNus19Dm32 };
96  seed_values[&kFitDelta24InPiUnitsSterile] = { 1.5 };
97  seed_values[&kFitSinSqTheta24Sterile] = { 0 };
98 
99  // Create surface
100  size_t part = 1275;
101  size_t njobs = 2500;
102  if (RunningOnGrid()){
103  part=JobNumber();
104  njobs=10;
105  }
106 
107  // Just run a single point on the surface interactively
108  CovMxSurface surf(Binning::LogUniform(50, 1e-4, 1), Binning::LogUniform(50, 1e-2, 100),
109  &multi_expt, calc, &kFitSinSqTheta24Sterile, &kFitDmSq41Sterile, fit_vars, seed_values,
110  {}, SystShifts::Nominal(), true, nullptr, preds, data, cosmic, part, njobs);
111 
112  // Open output file and save surface
113  TFile* f = TFile::Open("surface.root", "recreate");
114  surf.SaveTo(f->mkdir("surface"));
115  delete f;
116 
117 }
const double kNus19Dm32
Definition: Utilities.h:28
void FakeData()
Definition: rootlogon.C:156
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
size_t JobNumber()
What&#39;s the process number for a grid job?
Definition: Utilities.cxx:370
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void MakeSurface()
Definition: MakeSurface.C:31
bool RunningOnGrid()
Is this a grid (condor) job?
Definition: Utilities.cxx:363
const FitDmSq32Sterile kFitDmSq32Sterile
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
void SetAngles(osc::OscCalcSterile *calc)
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
const XML_Char const XML_Char * data
Definition: expat.h:268
const FitDelta24InPiUnitsSterile kFitDelta24InPiUnitsSterile
TString part[npart]
Definition: Style.C:32
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
expt
Definition: demo5.py:34
Sum up livetimes from individual cosmic triggers.
std::vector< const IExperiment * > GetConstraints()
Gaussian constrains on atmospheric parameters.
Definition: Utilities.h:416
Compare a single data spectrum to the MC + cosmics expectation.
Float_t d
Definition: plot.C:236
const covmx::Sample kNusFHCFarDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kFarDet)
std::vector< float > Spectrum
Definition: Constants.h:570
OStream cout
Definition: OStream.cxx:6
const Binning bins
Definition: NumuCC_CPiBin.h:8
Combine multiple component experiments.
const double kNus19Th23
Definition: Utilities.h:33
IPrediction * LoadPrediction(std::string detector, bool rhc=false, std::string syst_type="all")
Function to load prediction object.
static Binning LogUniform(int n, double lo, double hi)
Definition: Binning.cxx:134
std::unique_ptr< Spectrum > LoadCosmic(covmx::Sample sample, bool cvmfs=true)
Get cosmics for a given sample.
Definition: Utilities.h:145
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
const FitDmSq41Sterile kFitDmSq41Sterile
const double kNus19SensLivetime
Definition: Utilities.h:39
const FitTheta23Sterile kFitTheta23Sterile
surf
Definition: demo5.py:35