MakeSurfaceNoNDOsc.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"
17 #include "CAFAna/Analysis/CovMxSurface.h"
20 #include "CAFAna/Analysis/Fit.h"
21 
23 
24 #include <ctime>
25 
26 using namespace ana;
27 
28 void MakeSurfaceNoNDOsc(TString opt) {
29 
30  DontAddDirectory guard;
31 
32  // Define exposures
33  double potND = kAna2018SensitivityFHCNDPOT;
34  double potFD = kAna2018FHCPOT;
35  double livetimeFD = kAna2018FHCLivetime;
36 
37  /////////////////////////////////////////////////////////////////////
38  // //
39  // ARGUMENT PARSING //
40  // //
41  /////////////////////////////////////////////////////////////////////
42 
43  std::vector<std::string> polarity = { "fhc", "rhc" };
44  std::vector<bool> usePolarity = GetOptionConfig(opt, polarity);
45  std::vector<std::string> detector = { "neardet", "fardet" };
46  std::vector<bool> useDetector = GetOptionConfig(opt, detector);
47 
48  std::string optstring(opt.Data());
49 
50  std::vector<const ISyst*> isysts;
51  if (opt.Contains("covmxsysts") or opt.Contains("pullsyst")) {
52  if (useDetector[1]) isysts = LoadSysts(kNusFHCFarDet, kUseCVMFS);
53  if (useDetector[0]) isysts = LoadSysts(kNusFHCNearDet, kUseCVMFS);
54  }
55 
56  // Get pull and covariance matrix syst options
57  bool covmxSysts = CheckOption(opt, "covmxsysts");
58  std::string pullSyst = "none";
59  size_t pspos = optstring.find("pullsyst=");
60  if (pspos != std::string::npos) {
61  if (optstring.substr(pspos+8, 1) != "=")
62  assert(false and "The dm41 option must take the form \"dm41=xxx\".");
63  size_t start = pspos+9;
64  size_t end = optstring.find("_", pspos);
65  pullSyst = optstring.substr(start, end-start);
66  std::cout << "Treating syst " << pullSyst << " as a pull term." << std::endl;
67  }
68  if (pullSyst == "all" and covmxSysts) assert (false and
69  "Treating systematics as both pull terms and covariance matrix is not permitted.");
70 
71  // Load saved fake data, if requested
72  std::string fakedata;
73  size_t universe = 999999;
74  size_t fdpos = optstring.find("fakedata=");
75  if (optstring.substr(fdpos+8, 1) != "=")
76  assert(false and "The data option must take the form \"fakedata=xxx\".");
77  if (fdpos != std::string::npos) {
78  size_t start = fdpos + 9;
79  if (optstring.substr(start, 6) == "asimov") {
80  fakedata = "asimov";
81  }
82  else if (optstring.substr(start, 5) == "stats") {
83  fakedata = "stats";
84  size_t end = optstring.find("_", fdpos);
85  universe = std::stoi(optstring.substr(start+5, end-(start+5)));
86  }
87  else if (optstring.substr(start, 5) == "systs") {
88  fakedata = "systs";
89  size_t end = optstring.find("_", fdpos);
90  universe = std::stoi(optstring.substr(start+5, end-(start+5)));
91  }
92  else {
93  assert(false and "Fake data type not recognised!");
94  }
95  }
96 
97  // Set delta m, if requested
98  double dm41 = -1;
99  size_t dmpos = optstring.find("dm41=");
100  if (dmpos != std::string::npos) {
101  if (optstring.substr(dmpos+4, 1) != "=")
102  throw std::runtime_error("The dm41 option must take the form \"dm41=xxx\".");
103  size_t start = dmpos+5;
104  size_t end = optstring.find("_", dmpos);
105  if (end != std::string::npos)
106  dm41 = std::stod(optstring.substr(start, end-start));
107  else
108  dm41 = std::stod(optstring.substr(start));
109  std::cout << "Using dm41 value " << dm41 << std::endl;
110  }
111 
112  // Don't profile delta24
113  bool profDelta24 = not CheckOption(opt, "noprofdelta24");
114  bool profVars = not CheckOption(opt, "noprofvars");
115 
116  // Check what type of plot we're making
117  bool th24vsth34 = false;
118  bool th24vsdm41 = false;
119  bool th34vsdm41 = false;
120  if (CheckOption(opt, "th24vsth34")) th24vsth34 = true;
121  else if (CheckOption(opt, "th24vsdm41")) th24vsdm41 = true;
122  else if (CheckOption(opt, "th34vsdm41")) th34vsdm41 = true;
123  else throw std::runtime_error("The options string must contain plot type: \"th24vsth34\", \"th24vsdm41\" or \"th34vsdm41\".");
124 
125  // Make sure dm41 is set for th24 vs th34 case
126  if (th24vsth34 && dm41 == -1)
127  throw std::runtime_error("If you want to run a th24 vs th34 surface, you must specify a value for dm41.");
128 
129  /////////////////////////////////////////////////////////////////////
130  // //
131  // PREPARE INPUTS //
132  // //
133  /////////////////////////////////////////////////////////////////////
134 
135  // Oscillation calculator
137  SetNus20Params(calc);
138 
139  // Load predictions
140  std::vector<covmx::Sample> samples;
141  std::vector<IPrediction*> preds;
142  for (size_t d = 0; d < 2; ++d) { // for detector
143  for (size_t p = 0; p < 2; ++p) { // for beam polarity
144  if (!usePolarity[p] || !useDetector[d]) continue;
145  samples.push_back(covmx::Sample(covmx::Analysis::kNC,
147  preds.push_back(LoadPrediction(samples.back(), pullSyst == "none" ? "nosysts" : "systs", kUseCVMFS));
148  } // for beam polarity
149  } // for detector
150 
151  // Load cosmics
152  std::vector<Spectrum*> cosmic;
153  for (auto sample : samples)
154  cosmic.push_back(LoadCosmic(sample, kUseCVMFS).release());
155 
156  // Load data
157  std::vector<Spectrum*> data;
158  if (fakedata == "asimov") {
159  std::cout << "Using Asimov fake data." << std::endl;
160  for (size_t i = 0; i < preds.size(); ++i) {
161  double pot = (samples[i].detector == covmx::Detector::kNearDet) ? potND : potFDl
162  data.push_back(new Spectrum(preds[i]->Predict(calc).FakeData(pot)));
163  if(samples[i].detector == covmx::Detector::kFarDet)
164  data.back()->OverrideLivetime(livetimeFD);
165  if (cosmic[i]) *data[i] += *cosmic[i];
166  }
167  }
168  else {
169  std::cout << "Loading universe " << universe << " of " << fakedata
170  << " fake data." << std::endl;
171  data = LoadFakeData(samples, universe, (fakedata=="systs"), kUseCVMFS);
172  }
173 
174  // Pull term systematics
175  std::vector<const ISyst*> pullSysts = {};
176  SystConcat* sc = nullptr;
177  if (pullSyst == "all") {
178  // All pull systs, single detector
179  if (samples.size() == 1)
180  pullSysts = isysts;
181  // All pull systs, both detectors
182  else {
183  sc = new SystConcat(GetSystConcat(samples, isysts, preds));
184  pullSysts = sc->GetActiveSysts();
185  }
186  }
187  else if (pullSyst != "none") {
188  // Single pull syst, single detector
189  if (samples.size() == 1) {
190  for (const ISyst* syst : isysts) {
191  if (GetSystBaseName(syst->ShortName()) == pullSyst) {
192  pullSysts = { syst };
193  break;
194  }
195  }
196  if (pullSysts.empty()) {
197  std::ostringstream err;
198  err << "Couldn't find requested pull systematic " << pullSyst;
199  assert (false and err.str().c_str());
200  }
201  }
202  // Single pull syst, both detectors
203  else {
204  std::vector<SystConcat> all_scs = GetSystConcats(samples, preds, isysts);
205  for (auto it_sc : all_scs) {
206  if (GetSystBaseName(it_sc.GetActiveSysts()[0]->ShortName()) == pullSyst) {
207  sc = new SystConcat(it_sc);
208  pullSysts = sc->GetActiveSysts();
209  break;
210  }
211  }
212  }
213  }
214 
215  // Load covariance matrix
216  CovarianceMatrix* mx = nullptr;
217  if (covmxSysts) {
218  std::vector<Binning> bins;
219  for (auto d : data) bins.push_back(d->GetBinnings()[0]);
220  mx = new CovarianceMatrix(bins);
221  std::cout << std::endl;
222  for (CovarianceMatrix* it_mx : GetMatrices(samples, kUseCVMFS)) {
223  if (GetSystBaseName(it_mx->GetName()) != pullSyst) {
224  mx->AddMatrix(it_mx);
225  }
226  else
227  std::cout << "Omitting pull term " << pullSyst
228  << " from matrix." << std::endl;
229  }
230  }
231 
232  /////////////////////////////////////////////////////////////////////
233  // //
234  // MAKE SURFACE //
235  // //
236  /////////////////////////////////////////////////////////////////////
237 
238  // Create experiment
239  OscCovMxExperiment* expt = nullptr;
240  if (covmxSysts) expt = new OscCovMxExperiment(preds, mx, data, cosmic);
241  else expt = new OscCovMxExperiment(preds, data, cosmic);
242  if (samples.size() > 1 and sc) expt->ConfigureSysts(sc);
243  auto noosccalc = DefaultSterileCalc(4);
244  for (int i = 1; i <= 4; ++i) {
245  noosccalc->SetDm(i, 0);
246  }
247  expt->DisableOscillations({true, false}, noosccalc);
248 
249  // Set Gaussian constraints, make MultiExperiment
250  std::vector<const IChiSqExperiment*> expts = GetConstraints();
251  expts.push_back(expt);
252  MultiExperiment multiExpt(expts);
253 
254  // Set fixed dm41 for th24 vs th34 case
255  if (th24vsth34) calc->SetDm(4, dm41);
256 
257  // Binning, fit vars and seeds
258  std::vector<const IFitVar*> fitVars = { &kFitTheta23Sterile,
260  if (profDelta24) fitVars.push_back(&kFitDelta24InPiUnitsSterile);
261  else calc->SetDelta(2, 4, 0.5);
262  std::map<const IFitVar*, std::vector<double> > seedValues;
263  seedValues[&kFitTheta23Sterile] = { kNus19Th23 };
264  seedValues[&kFitDmSq32Sterile] = { kNus19Dm32 };
265  if (profDelta24) seedValues[&kFitDelta24InPiUnitsSterile] = { 0.5 };
266 
267  Binning* xbins = nullptr;
268  Binning* ybins = nullptr;
269 
270  std::vector<const IFitVar*> xyVars;
271 
272  if (th24vsth34) {
273  xbins = new Binning(Binning::Simple(50, 0, 50));
274  ybins = new Binning(Binning::Simple(50, 0, 45));
275  xyVars.push_back(&kFitTheta34InDegreesSterile);
276  xyVars.push_back(&kFitTheta24InDegreesSterile);
277  }
278  else if (th24vsdm41) {
279  xbins = new Binning(Binning::LogUniform(50, 1e-4, 1));
280  ybins = new Binning(Binning::LogUniform(50, 1e-2, 100));
281  xyVars.push_back(&kFitSinSqTheta24Sterile);
282  xyVars.push_back(&kFitDmSq41Sterile);
283  fitVars.push_back(&kFitSinSqTheta34Sterile);
284  seedValues[&kFitSinSqTheta34Sterile] = { 1. };
285  }
286  else if (th34vsdm41) {
287  xbins = new Binning(Binning::Simple(50, 0, 1));
288  ybins = new Binning(Binning::LogUniform(50, 1e-2, 100));
289  xyVars.push_back(&kFitSinSqTheta34Sterile);
290  xyVars.push_back(&kFitDmSq41Sterile);
291  fitVars.push_back(&kFitSinSqTheta24Sterile);
292  seedValues[&kFitSinSqTheta24Sterile] = { 1. };
293  }
294 
295  // Create surface
296  size_t part = 127;
297  size_t njobs = 250;
298  if (RunningOnGrid()){
299  part=JobNumber();
300  njobs=10;
301  }
302 
303  // Override fit variables if requested
304  if (not profVars) {
305  fitVars.clear();
306  seedValues.clear();
307  }
308 
309  CovMxSurface surf(*xbins, *ybins, &multiExpt, calc, xyVars[0], xyVars[1],
310  fitVars, seedValues, pullSysts, SystShifts::Nominal(), true, nullptr,
311  preds, data, cosmic);//, part, njobs);
312 
313  // Open output file and save surface
314  std::string outfile = optstring + ".root";
315  TFile* f = new TFile(outfile.c_str(), "recreate");
316  surf.SaveTo(f->mkdir(optstring.c_str()));
317  delete f;
318 
319 }
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
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
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
const char * p
Definition: xmltok.h:285
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
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.
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
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.
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
std::vector< bool > GetOptionConfig(TString opt, std::vector< std::string > options)
#define pot
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< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
Definition: Utilities.h:350
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.
const double kNus19Th23
Definition: Utilities.h:33
::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::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
void MakeSurfaceNoNDOsc(TString opt)
static Binning Simple(int n, double lo, double hi, const std::vector< std::string > &labels={})
Definition: Binning.cxx:114
const covmx::Sample kNusFHCNearDet(covmx::Analysis::kNC, covmx::Polarity::kFHC, covmx::Detector::kNearDet)
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
Definition: Exposures.h:213
FILE * outfile
Definition: dump_event.C:13
const FitTheta23Sterile kFitTheta23Sterile
surf
Definition: demo5.py:35