MakeSurfaceLLTest.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"
20 
22 
23 #include <ctime>
24 
25 using namespace ana;
26 
27 void MakeSurfaceLLTest(TString opt) {
28 
29  DontAddDirectory guard;
30 
31  // Define exposures
32  double potND = kAna2018SensitivityFHCNDPOT;
33  double potFD = kAna2018FHCPOT;
34  double livetimeFD = kAna2018FHCLivetime;
35 
36  /////////////////////////////////////////////////////////////////////
37  // //
38  // ARGUMENT PARSING //
39  // //
40  /////////////////////////////////////////////////////////////////////
41 
42  std::string optstring(opt.Data());
43 
44  // Set parallelising options
45  size_t njobs = 999;
46  size_t pos = optstring.find("njobs=");
47  if (pos != std::string::npos) {
48  size_t start = pos + 6;
49  size_t end = optstring.find("_", pos);
50  njobs = std::stoi(optstring.substr(start, end-start));
51  std::cout << "Splitting the surface into " << njobs << " jobs." << std::endl;
52  }
53 
54  // Get pull and covariance matrix syst options
55  bool covmxSysts = CheckOption(opt, "covmxsysts");
56  std::string pullSyst = "none";
57  pos = optstring.find("pullsyst=");
58  if (pos != std::string::npos) {
59  size_t start = pos + 9;
60  size_t end = optstring.find("_", pos);
61  pullSyst = optstring.substr(start, end-start);
62  std::cout << "Treating syst " << pullSyst << " as a pull term." << std::endl;
63  }
64  if (pullSyst == "all" and covmxSysts) assert (false and
65  "Treating systematics as both pull terms and covariance matrix is not permitted.");
66 
67  // Get samples, predictions and systematics
68  std::vector<covmx::Sample> samples = GetSamplesFromOptString(opt);
69  std::vector<IPrediction*> preds;
70  std::vector<std::vector<const ISyst*>> isysts;
71  for (covmx::Sample sample : samples) {
72  isysts.push_back(LoadSysts(sample, kUseCVMFS));
73  preds.push_back(LoadPrediction(sample, pullSyst == "none" ? "nosysts" : "systs", kUseCVMFS));
74  }
75 
76  // Load saved fake data, if requested
77  std::string fakedata;
78  size_t universe = 999999;
79  pos = optstring.find("fakedata=");
80  if (pos != std::string::npos) {
81  size_t start = pos + 9;
82  if (optstring.substr(start, 6) == "asimov") {
83  fakedata = "asimov";
84  }
85  else if (optstring.substr(start, 5) == "stats") {
86  fakedata = "stats";
87  size_t end = optstring.find("_", pos);
88  universe = std::stoi(optstring.substr(start+5, end-(start+5)));
89  }
90  else if (optstring.substr(start, 5) == "systs") {
91  fakedata = "systs";
92  size_t end = optstring.find("_", pos);
93  universe = std::stoi(optstring.substr(start+5, end-(start+5)));
94  }
95  else {
96  assert(false and "Fake data type not recognised!");
97  }
98  }
99 
100  // Set delta m, if requested
101  double dm41 = -1;
102  pos = optstring.find("dm41=");
103  if (pos != std::string::npos) {
104  size_t start = pos+5;
105  size_t end = optstring.find("_", start);
106  dm41 = std::stod(optstring.substr(start, end-start));
107  std::cout << "Using dm41 value " << dm41 << std::endl;
108  }
109 
110  // Don't profile delta24
111  bool profDelta24 = not CheckOption(opt, "noprofdelta24");
112  bool profVars = not CheckOption(opt, "noprofvars");
113 
114  // Zoom y axis
115  std::tuple<bool, double, double> yrange = std::make_tuple(false, -1, -1);
116  pos = optstring.find("zoomyaxis");
117  if (pos != std::string::npos) {
118  // Using manually defined y axis range
119  std::get<0>(yrange) = true;
120  // Get lower bound
121  size_t start = pos+10;
122  size_t end = optstring.find("-", start);
123  // std::cout << "Lower bound: " << optstring.substr(start, end-start) << std::endl;
124  std::get<1>(yrange) = std::stod(optstring.substr(start, end-start));
125  // Get upper bound
126  start = end + 1;
127  end = optstring.find("_", start);
128  // std::cout << "Upper bound: " << optstring.substr(start, end-start) << std::endl;
129  std::get<2>(yrange) = std::stod(optstring.substr(start, end-start));
130  }
131 
132 
133  // Check what type of plot we're making
134  bool th24vsth34 = false;
135  bool th24vsdm41 = false;
136  bool th34vsdm41 = false;
137  if (CheckOption(opt, "th24vsth34")) th24vsth34 = true;
138  else if (CheckOption(opt, "th24vsdm41")) th24vsdm41 = true;
139  else if (CheckOption(opt, "th34vsdm41")) th34vsdm41 = true;
140  else throw std::runtime_error("The options string must contain plot type: \"th24vsth34\", \"th24vsdm41\" or \"th34vsdm41\".");
141 
142  // Make sure dm41 is set for th24 vs th34 case
143  if (th24vsth34 && dm41 == -1)
144  throw std::runtime_error("If you want to run a th24 vs th34 surface, you must specify a value for dm41.");
145 
146  /////////////////////////////////////////////////////////////////////
147  // //
148  // PREPARE INPUTS //
149  // //
150  /////////////////////////////////////////////////////////////////////
151 
152  // Oscillation calculator
154  SetNus20Params(calc);
155  PrintOscParams(calc);
156 
157  // Load cosmics
158  std::vector<Spectrum*> cosmic;
159  for (auto sample : samples) {
160  cosmic.push_back(LoadCosmic(sample, kUseCVMFS).release());
161  }
162 
163  // Load data
164  std::vector<Spectrum*> data;
165  if (fakedata == "asimov") {
166  std::cout << "Using Asimov fake data." << std::endl;
167  for (size_t i = 0; i < preds.size(); ++i) {
168  double pot = (samples[i].detector == covmx::Detector::kNearDet) ? potND : potFD;
169  data.push_back(new Spectrum(preds[i]->Predict(calc).FakeData(pot)));
170  if(samples[i].detector == covmx::Detector::kFarDet)
171  data.back()->OverrideLivetime(livetimeFD);
172  if (cosmic[i]) *data[i] += *cosmic[i];
173  }
174  }
175  else {
176  std::cout << "Loading universe " << universe << " of " << fakedata
177  << " fake data." << std::endl;
178  data = LoadFakeData(samples, universe, (fakedata=="systs"), kUseCVMFS);
179  }
180 
181  // Pull term systematics
182  std::vector<const ISyst*> pullSysts = {};
183  SystConcat* sc = nullptr;
184  if (pullSyst == "all") {
185  // All pull systs, single detector
186  if (samples.size() == 1)
187  pullSysts = isysts[0];
188  // All pull systs, both detectors
189  else {
190  sc = new SystConcat(GetSystConcat(samples, isysts, preds));
191  pullSysts = sc->GetActiveSysts();
192  }
193  }
194  else if (pullSyst != "none") {
195  // Single pull syst, single detector
196  if (samples.size() == 1) {
197  for (std::vector<const ISyst*> detSysts : isysts) {
198  for (const ISyst* syst : detSysts) {
199  if (GetSystBaseName(syst->ShortName()) == pullSyst) {
200  pullSysts = { syst };
201  break;
202  }
203  }
204  }
205  if (pullSysts.empty()) {
206  std::ostringstream err;
207  err << "Couldn't find requested pull systematic " << pullSyst;
208  assert (false and err.str().c_str());
209  }
210  }
211  // Single pull syst, both detectors
212  else {
213  std::vector<SystConcat> allSCs = GetSystConcats(samples, preds, isysts);
214  for (auto it_sc : allSCs) {
215  if (GetSystBaseName(it_sc.GetActiveSysts()[0]->ShortName()) == pullSyst) {
216  sc = new SystConcat(it_sc);
217  pullSysts = sc->GetActiveSysts();
218  break;
219  }
220  }
221  }
222  }
223 
224  // Load covariance matrix
225  CovarianceMatrix* mx = nullptr;
226  if (covmxSysts) {
227  std::vector<Binning> bins;
228  for (auto d : data) bins.push_back(d->GetBinnings()[0]);
229  mx = new CovarianceMatrix(bins);
230  std::cout << std::endl;
231  for (CovarianceMatrix* it_mx : GetMatrices(samples, kUseCVMFS)) {
232  if (GetSystBaseName(it_mx->GetName()) != pullSyst) {
233  mx->AddMatrix(it_mx);
234  }
235  else
236  std::cout << "Omitting pull term " << pullSyst
237  << " from matrix." << std::endl;
238  }
239  }
240 
241  /////////////////////////////////////////////////////////////////////
242  // //
243  // MAKE SURFACE //
244  // //
245  /////////////////////////////////////////////////////////////////////
246 
247  // Create experiment
248  OscCovMxExperiment* expt = nullptr;
249  if (covmxSysts) expt = new OscCovMxExperiment(preds, mx, data, cosmic);
250  else expt = new OscCovMxExperiment(preds, data, cosmic);
251  if (samples.size() > 1 and sc) expt->ConfigureSysts(sc);
252  expt->SetLogLikelihood(); // Enable log likelihood correction
253  // Set Gaussian constraints, make MultiExperiment
254 
255  std::vector<const IChiSqExperiment*> expts = GetConstraints();
256  expts.push_back(expt);
257  MultiExperiment multiExpt(expts);
258 
259  // Set fixed dm41 for th24 vs th34 case
260  if (th24vsth34) calc->SetDm(4, dm41);
261 
262  // Binning, fit vars and seeds
263  std::vector<const IFitVar*> fitVars = { &kFitTheta23Sterile,
265  if (profDelta24) fitVars.push_back(&kFitDelta24InPiUnitsSterile);
266  else calc->SetDelta(2, 4, 0.5);
267  std::map<const IFitVar*, std::vector<double> > seedValues;
268  seedValues[&kFitTheta23Sterile] = { kNus20Th23 };
269  seedValues[&kFitDmSq32Sterile] = { kNus20Dm32 };
270  if (profDelta24) seedValues[&kFitDelta24InPiUnitsSterile] = { 0.5 };
271 
272  Binning* xbins = nullptr;
273  Binning* ybins = nullptr;
274 
275  std::vector<const IFitVar*> xyVars;
276 
277  if (th24vsth34) {
278  xbins = new Binning(Binning::Simple(50, 0, 50));
279  ybins = new Binning(Binning::Simple(50, 0, 45));
280  xyVars.push_back(&kFitTheta34InDegreesSterile);
281  xyVars.push_back(&kFitTheta24InDegreesSterile);
282  }
283  else if (th24vsdm41) {
284  xbins = new Binning(Binning::LogUniform(50, 1e-4, 1));
285  // Zoom in y axis if requested
286  if (std::get<0>(yrange)) {
287  ybins = new Binning(Binning::Simple(50, std::get<1>(yrange), std::get<2>(yrange)));
288  }
289  else {
290  ybins = new Binning(Binning::LogUniform(50, 1e-2, 100));
291  }
292  xyVars.push_back(&kFitSinSqTheta24Sterile);
293  xyVars.push_back(&kFitDmSq41Sterile);
294  fitVars.push_back(&kFitSinSqTheta34Sterile);
295  seedValues[&kFitSinSqTheta34Sterile] = { 0., 1. };
296  }
297  else if (th34vsdm41) {
298  xbins = new Binning(Binning::Simple(50, 0, 1));
299  ybins = new Binning(Binning::LogUniform(50, 1e-2, 100));
300  xyVars.push_back(&kFitSinSqTheta34Sterile);
301  xyVars.push_back(&kFitDmSq41Sterile);
302  fitVars.push_back(&kFitSinSqTheta24Sterile);
303  seedValues[&kFitSinSqTheta24Sterile] = { 0., 1. };
304  }
305 
306  // Set up parallelising options
307  size_t part = 999;
308  if (njobs != 999) {
309  if (RunningOnGrid()) {
310  part = JobNumber();
311  }
312  else {
313  // Assume interactive test, set patch to the centre of the surface
314  part = njobs / 2;
315  }
316  }
317 
318  // Override fit variables if requested
319  if (not profVars) {
320  fitVars.clear();
321  seedValues.clear();
322  }
323 
324  CovMxSurface surf(*xbins, *ybins, &multiExpt, calc, xyVars[0], xyVars[1],
325  fitVars, seedValues, pullSysts, SystShifts::Nominal(), true, nullptr,
326  preds, sc, data, cosmic, part, njobs);
327 
328  // Open output file and save surface
329  std::string fitName = optstring.substr(0, optstring.find("_njobs"));
330  std::string outFileName = fitName + ".root";
331  TFile* outFile = new TFile(outFileName.c_str(), "recreate");
332  surf.SaveTo(outFile);
333  delete outFile;
334 
335 }
std::vector< SystConcat > GetSystConcats(std::vector< covmx::Sample > samples, std::vector< IPrediction * > preds, std::vector< const ISyst * > systs)
Get correctly concatenated systs for ND and FD.
Definition: Utilities.h:294
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
Represent the binning of a Spectrum&#39;s x-axis.
Definition: Binning.h:16
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::string GetSystBaseName(std::string name)
Get the last part of a systematic&#39;s ShortName.
Definition: Utilities.h:232
double stod(const std::string &val)
bool RunningOnGrid()
Is this a grid (condor) job?
Definition: Utilities.cxx:363
void SetDelta(int i, int j, double delta)
Adapt the PMNS_Sterile calculator to standard interface.
const FitDmSq32Sterile kFitDmSq32Sterile
const bool kUseCVMFS
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
void PrintOscParams(osc::OscCalcSterile *calc)
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
void SaveTo(TDirectory *dir) const
Standard SaveTo implementation.
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
void MakeSurfaceLLTest(TString opt)
expt
Definition: demo5.py:34
Sum up livetimes from individual cosmic triggers.
TFile * outFile
Definition: PlotXSec.C:135
void SetNus20Params(osc::OscCalcSterile *calc)
Definition: Utilities.h:682
std::vector< const IExperiment * > GetConstraints()
Gaussian constrains on atmospheric parameters.
Definition: Utilities.h:416
const int xbins
Definition: MakePlots.C:82
#define pot
Compare a single data spectrum to the MC + cosmics expectation.
Float_t d
Definition: plot.C:236
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
Definition: Utilities.h:350
std::vector< float > Spectrum
Definition: Constants.h:570
SystConcat GetSystConcat(std::vector< covmx::Sample > samples, std::vector< const ISyst * > systs, std::vector< IPrediction * > preds)
Get correctly concatenated systs for ND and FD.
Definition: Utilities.h:249
bool CheckOption(TString opts, TString opt)
OStream cout
Definition: OStream.cxx:6
const int ybins
const Binning bins
Definition: NumuCC_CPiBin.h:8
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
void SetDm(int i, double dm)
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::vector< covmx::Sample > GetSamplesFromOptString(TString optString)
Function to take an option TString and return a vector of associated covmx::Samples.
Definition: Utilities.h:384
std::unique_ptr< Spectrum > LoadCosmic(covmx::Sample sample, bool cvmfs=true)
Get cosmics for a given sample.
Definition: Utilities.h:145
assert(nhit_max >=nhit_nbins)
std::vector< const ISyst * > LoadSysts(covmx::Sample sample)
Get systematics for a given sample.
Definition: Utilities.h:524
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
std::map< std::string, ana::Spectrum > LoadFakeData(const std::string &fakeDataFilename, const std::string &path)
const double kAna2018FHCPOT
Definition: Exposures.h:207
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const FitTheta23Sterile kFitTheta23Sterile
surf
Definition: demo5.py:35