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"
22 
26 
27 #include <chrono>
28 
29 using namespace ana;
30 
31 using std::string;
32 using std::vector;
33 using std::ostringstream;
34 using std::cout;
35 using std::endl;
36 using std::exception;
37 
38 using std::chrono::high_resolution_clock;
39 using std::chrono::duration_cast;
40 using std::chrono::seconds;
41 
42 using covmx::Sample;
44 
45 void MakeSurface(TString opt, string systConfig="allsysts") {
46 
47  try {
48 
49  DontAddDirectory guard;
50 
51  /////////////////////////////////////////////////////////////////////
52  // //
53  // ARGUMENT PARSING //
54  // //
55  /////////////////////////////////////////////////////////////////////
56 
57  FitConfig cfg = GetConfig(opt);
58  if (cfg.samples.size() == 0) {
59  throw std::runtime_error("No samples found! Exiting.");
60  }
61 
62  for (Sample& s : cfg.samples) {
63  cout << "Setting prediction for sample " << s.GetName() << "..." << endl;
64  SetPrediction(s, false);
65  }
66 
67  /////////////////////////////////////////////////////////////////////
68  // //
69  // PREPARE INPUTS //
70  // //
71  /////////////////////////////////////////////////////////////////////
72 
73  // Oscillation calculator
75  SetNus20Params(calc, cfg.nullHyp);
76  if (CheckOption(opt, "th24vsth34")) calc->SetDm(4, 0.5);
77 
78  // Load cosmics
79  for (Sample& sample : cfg.samples) SetCosmic(sample);
80 
81  cout << cfg.xAxis->var->ShortName() << " " << cfg.yAxis->var->ShortName() << endl;
82 
83  // Load data
84  if (cfg.fakeDataType == kAsimov) {
85  cout << "Fitting with Asimov fake data with \"" << cfg.nullHyp << "\" null hypothesis:" << endl;
86  PrintOscParams(calc);
87  for (Sample& sample: cfg.samples) {
88  Spectrum* data = new Spectrum(sample.GetPrediction()->Predict(calc).FakeData(sample.GetPOT()));
89  if (sample.detector == covmx::kFarDet) data->OverrideLivetime(sample.GetLivetime());
90  if (sample.HasCosmic()) *data += *sample.GetCosmic();
91  sample.SetData(data);
92  }
93  SetNus20Params(calc, cfg.nullHyp);
94  if (CheckOption(opt, "th24vsth34")) calc->SetDm(4, 0.5);
95  } else if (cfg.fakeDataType == kPseudoexperiment) {
96  cout << "Loading fake data, type " << cfg.fakeDataSet << ", pseudoexperiment "
97  << cfg.fakeDataUniverse << endl;
98  ostringstream filepath;
99  filepath << InputPath() << "/fakedata/fakedata_" << cfg.fakeDataSet << ".root";
100  TFile* f = TFile::Open(filepath.str().c_str(), "read");
101  // Get universe and spectra
102  TDirectory* dir = f->GetDirectory(Form("universe_%u", cfg.fakeDataUniverse));
103  for (Sample& s : cfg.samples) {
104  s.SetData(Spectrum::LoadFrom(dir, s.GetTag()).release());
105  }
106  delete f; // Close fake data file
107  }
108 
109  // Load covariance matrix
110  CovarianceMatrix* mx = nullptr;
111  if (cfg.systType == kCovMxSysts) {
112  TFile* f = TFile::Open((InputPath()+"/covmx/mx.root").c_str(), "read");
113  ostringstream tag;
114  for (size_t i = 0; i < cfg.samples.size(); ++i) {
115  if (i > 0) tag << "_";
116  tag << cfg.samples[i].GetID();
117  }
118  mx = LoadFrom<CovarianceMatrix>(f, tag.str()+"/"+systConfig).release();
119  delete f;
120  }
121 
122  /////////////////////////////////////////////////////////////////////
123  // //
124  // MAKE SURFACE //
125  // //
126  /////////////////////////////////////////////////////////////////////
127 
128  // Create experiment
129  IExperiment* expt(nullptr);
130  if (cfg.fitType == kPoissonLL) {
132  expt = ll;
133  ll->SaveHists();
134  if (!opt.Contains("quiet")) ll->SetVerbose(true);
135  ll->SetResetShifts(false);
136  } else {
137  expt = new OscCovMxExperiment(cfg.samples, mx);
138  if (cfg.fitType == kCNP) dynamic_cast<OscCovMxExperiment*>(expt)->SetCNP();
139  }
140 
141  // Set Gaussian constraints, make MultiExperiment
142  vector<const IExperiment*> expts = { expt };
143  if (!opt.Contains("numusel")) {
144  cout << "No numu samples included! Adding Gaussian constraints on atmospheric parameters." << endl;
145  for (const IExperiment* e : GetConstraints()) expts.push_back(e);
146  } else {
147  cout << "Numu samples included! Omitting Gaussian constraints on atmospheric parameters." << endl;
148  }
149  MultiExperiment multiExpt(expts);
150  TH1D hTime("hTime", "", 1, 0, 1);
152  FrequentistSurface surf(&multiExpt, calc, *cfg.xAxis, *cfg.yAxis,
153  cfg.fitVars, {}, SeedList(cfg.fitSeeds));
155  hTime.SetBinContent(1, duration_cast<seconds>(end-start).count());
156 
157  // Open output file and save surface
158  string optstring = opt.Data();
159  TFile* outFile = new TFile((optstring+".root").c_str(), "recreate");
160  surf.SaveTo(outFile, "surface");
161  outFile->WriteTObject(&hTime, "time");
162  delete outFile;
163 
164  } catch(const exception& e) {
165  cout << endl << "Exception occurred! " << e.what() << endl;
166  }
167 
168 } // 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:35
void SetNus20Params(osc::OscCalcSterile *calc, std::string type="3flav")
Definition: Utilities.h:619
double lambda
Definition: FitUtilities.h:52
void SetPrediction(covmx::Sample &sample, bool systs=true)
Definition: Utilities.h:503
const IFitVar * var
Definition: FitAxis.h:17
Adapt the PMNS_Sterile calculator to standard interface.
unsigned int fakeDataUniverse
Definition: FitUtilities.h:39
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
void PrintOscParams(osc::OscCalcSterile *calc)
osc::OscCalcDumb calc
FakeDataType fakeDataType
Definition: FitUtilities.h:37
std::map< const IFitVar *, std::vector< double > > fitSeeds
Definition: FitUtilities.h:47
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const XML_Char const XML_Char * data
Definition: expat.h:268
FitAxis * xAxis
Definition: FitUtilities.h:43
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
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
void SetCosmic(covmx::Sample &sample)
Set cosmics in a given sample.
Definition: Utilities.h:488
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:236
std::vector< float > Spectrum
Definition: Constants.h:759
FitType fitType
Definition: FitUtilities.h:49
bool CheckOption(TString opts, TString opt)
Compare a single data spectrum to the MC + cosmics expectation.
OStream cout
Definition: OStream.cxx:6
string InputPath()
double epsilon
Definition: FitUtilities.h:51
SystType systType
Definition: FitUtilities.h:31
Combine multiple component experiments.
TString GetConfig(int dbcfg)
Definition: PlotDB_Web.C:247
void SetDm(int i, double dm)
const std::string & ShortName() const
Definition: IFitVar.h:36
TDirectory * dir
Definition: macro.C:5
void SaveTo(TDirectory *dir, const std::string &name) const
Base class defining interface for experiments.
Definition: IExperiment.h:14
FitAxis * yAxis
Definition: FitUtilities.h:44
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:41
std::string fakeDataSet
Definition: FitUtilities.h:38
std::vector< const IFitVar * > fitVars
Definition: FitUtilities.h:46
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...
surf
Definition: demo5.py:35
enum BeamMode string