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 covmx analysis
9  */
10 
12 #include "CAFAna/Core/Progress.h"
13 #include "CAFAna/Core/Sample.h"
21 
22 #include "NuXAna/Core/Utilities.h"
26 
27 #include <chrono>
28 
29 using namespace ana;
30 
31 using std::cout;
32 using std::endl;
33 using std::ostringstream;
34 using std::string;
35 using std::vector;
36 
37 using std::chrono::high_resolution_clock;
38 using std::chrono::duration_cast;
39 using std::chrono::seconds;
40 
41 using covmx::Sample;
43 
44 void MakeSurface(TString opt, string systConfig="allsysts") {
45 
46  DontAddDirectory guard;
47 
48  /////////////////////////////////////////////////////////////////////
49  // //
50  // ARGUMENT PARSING //
51  // //
52  /////////////////////////////////////////////////////////////////////
53 
54  FitConfig cfg = GetConfig(opt);
55  if (cfg.samples.size() == 0) {
56  throw std::runtime_error("No samples found! Exiting.");
57  }
58 
59  for (Sample& s : cfg.samples) {
60  cout << "Setting prediction for sample " << s.GetName() << "..." << endl;
61  bool usePredInterp = cfg.fakeDataType == kPseudoexperiment && cfg.systType == kCovMxSysts;
62  SetPrediction(s, usePredInterp);
63  SetCosmic(s);
64  if (usePredInterp) SetSystAlias(s);
65  }
66 
67  /////////////////////////////////////////////////////////////////////
68  // //
69  // PREPARE INPUTS //
70  // //
71  /////////////////////////////////////////////////////////////////////
72 
73  cout << cfg.xAxis->var->ShortName() << " " << cfg.yAxis->var->ShortName() << endl;
74 
75  // Oscillation calculator
77  SetNus20Params(calc, cfg.nullHyp);
78  PrintOscParams(calc);
79 
80  // Load data
81  if (cfg.fakeDataType == kAsimov) {
82  cout << "Fitting with Asimov fake data with \"" << cfg.nullHyp << "\" null hypothesis:" << endl;
83  for (Sample& sample: cfg.samples) {
84  Spectrum* data = new Spectrum(sample.GetPrediction()->Predict(calc).FakeData(sample.GetPOT()));
85  if (sample.detector == covmx::kFarDet) data->OverrideLivetime(sample.GetLivetime());
86  if (sample.HasCosmic()) *data += *sample.GetCosmic();
87  sample.SetData(data);
88  }
89  } else if (cfg.fakeDataType == kPseudoexperiment) {
90  cout << "Fitting with fluctuated fake data \"" << cfg.nullHyp
91  << "\" null hypothesis, pseudoexperiment " << cfg.fakeDataUniverse << endl;
93  }
94  if (CheckOption(opt, "th24vsth34")) calc->SetDm(4, 0.5);
95 
96  // Load covariance matrix
97  CovarianceMatrix* mx = nullptr;
98  if (cfg.systType == kCovMxSysts) {
99  TFile* f = TFile::Open((InputPath()+"/covmx/mx.root").c_str(), "read");
100  ostringstream tag;
101  for (size_t i = 0; i < cfg.samples.size(); ++i) {
102  if (i > 0) tag << "_";
103  tag << cfg.samples[i].GetID();
104  }
105  mx = LoadFrom<CovarianceMatrix>(f, tag.str()+"/"+systConfig).release();
106  delete f;
107  }
108 
109  /////////////////////////////////////////////////////////////////////
110  // //
111  // MAKE SURFACE //
112  // //
113  /////////////////////////////////////////////////////////////////////
114 
115  // Create experiment
116  IExperiment* expt(nullptr);
117  if (cfg.fitType == kPoissonLL) {
119  expt = ll;
120  if (!opt.Contains("quiet")) ll->SetVerbose(true);
121  } else {
122  expt = new OscCovMxExperiment(cfg.samples, mx);
123  if (cfg.fitType == kCNP) dynamic_cast<OscCovMxExperiment*>(expt)->SetCNP();
124  }
125 
126  // Set Gaussian constraints, make MultiExperiment
127  vector<const IExperiment*> expts = { expt };
128  if (!opt.Contains("numusel")) {
129  cout << "No numu samples included! Adding Gaussian constraints on atmospheric parameters." << endl;
130  for (const IExperiment* e : GetConstraints()) expts.push_back(e);
131  } else {
132  cout << "Numu samples included! Adding loose constraint on dmsq32." << endl;
133  expts.push_back(GetConstraintLoose());
134  }
135  MultiExperiment multiExpt(expts);
136  TH1D hTime("hTime", "", 1, 0, 1);
138  FrequentistSurface surf(&multiExpt, calc, *cfg.xAxis, *cfg.yAxis,
139  cfg.fitVars, {}, SeedList(cfg.fitSeeds));
141  hTime.SetBinContent(1, duration_cast<seconds>(end-start).count());
142 
143  // Open output file and save surface
144  string optstring = opt.Data();
145  TFile* outFile = new TFile((optstring+".root").c_str(), "recreate");
146  surf.SaveTo(outFile, "surface");
147  outFile->WriteTObject(&hTime, "time");
148  delete outFile;
149 
150 } // main function
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void MakeSurface()
Definition: MakeSurface.C:31
std::vector< covmx::Sample > samples
Definition: FitUtilities.h:31
virtual void SetDm(int i, double dm)=0
double lambda
Definition: FitUtilities.h:47
void SetSystAlias(Sample &sample)
void SetPrediction(covmx::Sample &sample, bool systs=true)
Definition: Utilities.h:411
const IFitVar * var
Definition: FitAxis.h:17
std::vector< double > Spectrum
Definition: Constants.h:746
void SetNus20Params(osc::IOscCalcSterile *calc, std::string type)
Definition: Calcs.cxx:29
unsigned int fakeDataUniverse
Definition: FitUtilities.h:34
osc::OscCalcDumb calc
FakeDataType fakeDataType
Definition: FitUtilities.h:33
std::map< const IFitVar *, std::vector< double > > fitSeeds
Definition: FitUtilities.h:42
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:39
const XML_Char const XML_Char * data
Definition: expat.h:268
const IExperiment * GetConstraintLoose()
Loose Gaussian constraint on atmospheric mass splitting.
Definition: Utilities.h:469
FitAxis * xAxis
Definition: FitUtilities.h:38
Log-likelihood scan across two parameters.
const XML_Char * s
Definition: expat.h:262
osc::OscCalcSterile * DefaultSterileCalc(int nflavors)
Create a sterile calculator with default assumptions for all parameters.
Definition: Calcs.cxx:97
expt
Definition: demo5.py:34
TFile * outFile
Definition: PlotXSec.C:135
std::vector< const IExperiment * > GetConstraints()
Gaussian constrains on atmospheric parameters.
Definition: Utilities.h:416
void PrintOscParams(osc::IOscCalcSterile *calc)
Definition: Calcs.cxx:12
void SetCosmic(covmx::Sample &sample)
Set cosmics in a given sample.
Definition: Utilities.h:396
Compare a single data spectrum to the MC + cosmics expectation.
void OverrideLivetime(double newlive)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:235
base class for sterile oscillation calculators In the context of a sterile oscillation calculator...
FitType fitType
Definition: FitUtilities.h:44
bool CheckOption(TString opts, TString opt)
Compare a single data spectrum to the MC + cosmics expectation.
OStream cout
Definition: OStream.cxx:6
string InputPath()
Definition: Utilities.cxx:12
double epsilon
Definition: FitUtilities.h:46
SystType systType
Definition: FitUtilities.h:27
Combine multiple component experiments.
TString GetConfig(int dbcfg)
Definition: PlotDB_Web.C:247
const std::string & ShortName() const
Definition: IFitVar.h:37
void SaveTo(TDirectory *dir, const std::string &name) const
Base class defining interface for experiments.
Definition: IExperiment.h:14
FitAxis * yAxis
Definition: FitUtilities.h:39
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
std::string nullHyp
Definition: FitUtilities.h:36
std::vector< const IFitVar * > fitVars
Definition: FitUtilities.h:41
void GeneratePseudoexpt(std::vector< covmx::Sample > &samples, osc::IOscCalcSterile *calc, int seed, bool systs)
Create fluctuated pseudoexperiment.
Definition: Utilities.h:499
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...
surf
Definition: demo5.py:35
enum BeamMode string