MakeSurfaceBinningStudy.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 
24 
25 #include <ctime>
26 
27 using namespace ana;
28 
29 using std::back_inserter;
30 using std::begin;
31 using std::endl;
32 using std::end;
33 using std::endl;
34 using std::exception;
35 using std::ostringstream;
36 using std::runtime_error;
37 using std::string;
38 using std::transform;
39 using std::vector;
40 
42 
44 using covmx::Sample;
45 
46 void MakeSurfaceBinningStudy(TString opt, size_t ndPred, size_t fdPred) {
47 
48  try {
49 
50  DontAddDirectory guard;
51 
52  /////////////////////////////////////////////////////////////////////
53  // //
54  // ARGUMENT PARSING //
55  // //
56  /////////////////////////////////////////////////////////////////////
57 
58  FitConfig cfg = GetConfig(opt);
59 
60  if (ndPred < 1 || ndPred > 6) {
61  ostringstream err;
62  err << "Near detector binning option must be an integer between 1 and 6!"
63  << " Value provided was " << ndPred;
64  throw runtime_error(err.str());
65  }
66  if (fdPred < 1 || fdPred > 6) {
67  ostringstream err;
68  err << "Far detector binning option must be an integer between 1 and 6!"
69  << " Value provided was " << fdPred;
70  throw runtime_error(err.str());
71  }
72 
73 
74  //string path = "/cvmfs/nova.osgstorage.org/analysis/nux/binning_studies/";
75  string path = pnfs2xrootd("/pnfs/nova/persistent/analysis/nux/binning_studies/");
76  vector<string> suffix = { "", "_res10", "_res20", "_res30", "_res40", "_res50" };
77  for (Sample& s : cfg.samples) {
79  if (s.detector == covmx::kNearDet) {
80  string predFileName = path+"pred/pred_"+s.GetTag()+suffix[ndPred-1]+".root";
81  cout << "Loading prediction from " << predFileName << endl;
82  TFile* f = TFile::Open(predFileName.c_str(), "read");
83  pred = LoadFromFile<NDPredictionSterile>(f, "pred_"+s.GetTag()).release();
84  delete f;
85  SetAxis(s, 10*(ndPred-1));
86  } else {
87  string predFileName = path+"pred/pred_"+s.GetTag()+suffix[fdPred-1]+".root";
88  cout << "Loading prediction from " << predFileName << endl;
89  TFile* f = TFile::Open(predFileName.c_str(), "read");
90  pred = LoadFromFile<FDPredictionSterile>(f, "pred_"+s.GetTag()).release();
91  delete f;
92  SetAxis(s, 10*(fdPred-1));
93  }
94  s.SetPrediction(pred);
95  }
96 
97  /////////////////////////////////////////////////////////////////////
98  // //
99  // PREPARE INPUTS //
100  // //
101  /////////////////////////////////////////////////////////////////////
102 
103  // Oscillation calculator
105  SetNus20Params(calc);
106 
107  // Load cosmics
108  for (Sample& s : cfg.samples) {
109  if (s.detector == covmx::kFarDet) {
110  string cosmicFileName = path+"cosmics/"+s.GetTag()+suffix[fdPred-1]+".root";
111  cout << "Loading cosmic spectrum from " << cosmicFileName << endl;
112  TFile* f = TFile::Open(cosmicFileName.c_str(), "read");
113  s.SetCosmic(LoadFromFile<Spectrum>(f, s.GetTag()).release());
114  delete f;
115  }
116  }
117 
118  cout << cfg.xVar->ShortName() << " " << cfg.yVar->ShortName() << endl;
119 
120  // Load data
121  if (cfg.fakeDataType == kAsimov || cfg.fakeDataType == kFakeSignal) {
122  if (cfg.fakeDataType == kFakeSignal) {
123  cout << "Using fake signal oscillation parameter set " << cfg.fakeDataSet << endl;
124  SetFakeSignalParams(calc, cfg.fakeDataSet);
125  } else {
126  cout << "Using Asimov fake data." << endl;
127  }
128  cout << "Oscillation parameters used for fake data:" << endl;
129  PrintOscParams(calc);
130  for (Sample& sample: cfg.samples) {
131  Spectrum* data = new Spectrum(sample.GetPrediction()->Predict(calc).FakeData(sample.GetPOT()));
132  if (sample.detector == covmx::kFarDet) data->OverrideLivetime(sample.GetLivetime());
133  TH1D* tmp = data->ToTH1(9e20);
134  if (sample.HasCosmic()) *data += *sample.GetCosmic();
135  tmp = data->ToTH1(9e20);
136  sample.SetData(data);
137  }
138  SetNus20Params(calc);
139  }
140 
141  else if (cfg.fakeDataType == kFakeSignalFluctuated) {
142 
143  cout << "Loading fake signal data, type " << cfg.fakeDataSet << ", pseudoexperiment "
144  << cfg.fakeDataUniverse << endl;
145 
146  ostringstream filepath;
147  filepath << kNus20Path << "/fakedata/fhc_fake_signal_" << cfg.fakeDataSet << ".root";
148  TFile* f = TFile::Open(filepath.str().c_str(), "read");
149 
150  // Get universe and spectra
151  TDirectory* dir = f->GetDirectory(Form("universe_%u", cfg.fakeDataUniverse));
152  for (Sample& s : cfg.samples) {
153  s.SetData(Spectrum::LoadFrom(dir->GetDirectory(s.GetTag().c_str())).release());
154  }
155  delete f; // Close fake data file
156  }
157 
158  // Load covariance matrix
159  CovarianceMatrix* mx = nullptr;
160  if (cfg.systType == kCovMxSysts) {
161  mx = new CovarianceMatrix(cfg.samples);
162  for (CovarianceMatrix* itMx : GetMatrices(cfg.samples, kUseCVMFS)) {
163  bool use = false;
164  if (cfg.systConfig == kAllSysts) use = true;
165  else if (cfg.systConfig == kDetectorSysts) {
166  for (string name : kDetectorSystNames) {
167  if (itMx->GetName() == name) {
168  use = true;
169  break;
170  }
171  } // for detector syst name
172  } else if (cfg.systConfig == kBeamSysts) {
173  for (string name : kBeamSystNames) {
174  if (itMx->GetName() == name) {
175  use = true;
176  break;
177  }
178  } // for beam syst name
179  } else if (cfg.systConfig == kXsecSysts) {
180  for (const ISyst* syst : getAllXsecSysts_2018_RPAFix()) {
181  if (itMx->GetName() == syst->ShortName()) {
182  use = true;
183  break;
184  }
185  } // for xsec isyst
186  }
187  if (use) mx->AddMatrix(itMx);
188  delete itMx;
189  }
190  } else if (cfg.systType == kCovMxSingleSyst) {
191  for (CovarianceMatrix* itMx : GetMatrices(cfg.samples, kUseCVMFS)) {
192  if (itMx->GetName() == cfg.systName) {
193  mx = itMx;
194  } else {
195  delete itMx;
196  }
197  }
198  if (!mx) { // Throw an error if we didn't find the right matrix
199  ostringstream err;
200  err << "Couldn't find matrix with name \"" << cfg.systName << "\"!";
201  throw std::runtime_error(err.str());
202  }
203  } else if (cfg.systType == kCovMxWithSinglePullSyst) {
204  mx = new covmx::CovarianceMatrix(cfg.samples);
205  for (CovarianceMatrix* itMx : GetMatrices(cfg.samples, kUseCVMFS)) {
206  if (itMx->GetName() != cfg.systName) {
207  mx->AddMatrix(itMx);
208  }
209  delete itMx;
210  }
211  }
212  /////////////////////////////////////////////////////////////////////
213  // //
214  // MAKE SURFACE //
215  // //
216  /////////////////////////////////////////////////////////////////////
217 
218  // Create experiment
219  IChiSqExperiment* expt(nullptr);
220  if (cfg.fitType == kPoissonLL) {
221  expt = new LikelihoodCovMxExperiment(cfg.samples, mx, cfg.epsilon);
222  } else {
223  expt = new OscCovMxExperiment(cfg.samples, mx);
224  if (cfg.fitType == kCNP) dynamic_cast<OscCovMxExperiment*>(expt)->SetCNP();
225  }
226 
227  // Set Gaussian constraints, make MultiExperiment
228  vector<const IChiSqExperiment*> expts = GetConstraints();
229  expts.push_back(expt);
230  MultiExperiment multiExpt(expts);
231 
232  CovMxSurface surf(*cfg.xBins, *cfg.yBins, &multiExpt, calc, cfg.xVar, cfg.yVar,
233  cfg.fitVars, cfg.fitSeeds, {} /*FIX THIS!!*/, SystShifts::Nominal(),
234  cfg.jobNumber, cfg.nJobs);
235 
236  // Open output file and save surface
237  string optstring = opt.Data();
238  string fitName = optstring.substr(0, optstring.find("_njobs"));
239  string outFileName = fitName + ".root";
240  TFile* outFile = TFile::Open(outFileName.c_str(), "recreate");
241  surf.SaveTo(outFile);
242 
243  } catch(const exception& e) {
244  cout << endl << "Exception occurred! " << e.what() << endl;
245  }
246 
247 } // main function
248 
const IFitVar * yVar
Definition: FitUtilities.h:52
const XML_Char * name
Definition: expat.h:151
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const std::vector< std::string > kBeamSystNames
Definition: SystUtilities.h:15
std::vector< covmx::Sample > samples
Definition: FitUtilities.h:44
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
Adapt the PMNS_Sterile calculator to standard interface.
unsigned int fakeDataUniverse
Definition: FitUtilities.h:48
::xsd::cxx::tree::exception< char > exception
Definition: Database.h:225
const bool kUseCVMFS
Float_t tmp
Definition: plot.C:36
std::string pnfs2xrootd(std::string loc, bool unauth)
Definition: UtilsExt.cxx:237
void PrintOscParams(osc::OscCalcSterile *calc)
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
FakeDataType fakeDataType
Definition: FitUtilities.h:46
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
void SaveTo(TDirectory *dir) const
Standard SaveTo implementation.
std::map< const IFitVar *, std::vector< double > > fitSeeds
Definition: FitUtilities.h:56
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
void SetAxis(TH1 *h, bool visible, int color, int style=1, bool scale=true)
const XML_Char const XML_Char * data
Definition: expat.h:268
std::string systName
Definition: FitUtilities.h:40
const XML_Char * s
Definition: expat.h:262
void SetFakeSignalParams(osc::OscCalcSterile *calc, size_t type)
Definition: Utilities.h:702
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
void SetNus20Params(osc::OscCalcSterile *calc)
Definition: Utilities.h:682
std::vector< const IExperiment * > GetConstraints()
Gaussian constrains on atmospheric parameters.
Definition: Utilities.h:416
SystConfig systConfig
Definition: FitUtilities.h:39
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:546
const std::string kNus20Path
void MakeSurfaceBinningStudy(TString opt, size_t ndPred, size_t fdPred)
Compare a single data spectrum to the MC + cosmics expectation.
unsigned int nJobs
Definition: FitUtilities.h:37
void OverrideLivetime(double newlive)
DO NOT USE UNLESS YOU ARE 110% CERTAIN THERE ISN&#39;T A BETTER WAY!
Definition: Spectrum.h:228
std::vector< CovarianceMatrix * > GetMatrices(std::vector< covmx::Sample > samples, bool cvmfs=true)
Definition: Utilities.h:350
unsigned int fakeDataSet
Definition: FitUtilities.h:47
std::vector< float > Spectrum
Definition: Constants.h:570
FitType fitType
Definition: FitUtilities.h:58
Compare a single data spectrum to the MC + cosmics expectation.
unsigned int jobNumber
Definition: FitUtilities.h:36
OStream cout
Definition: OStream.cxx:6
double epsilon
Definition: FitUtilities.h:42
const std::string path
Definition: plot_BEN.C:43
SystType systType
Definition: FitUtilities.h:38
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TString GetConfig(int dbcfg)
Definition: PlotDB_Web.C:247
std::vector< const ISyst * > getAllXsecSysts_2018_RPAFix()
const std::string & ShortName() const
Definition: IFitVar.h:36
TDirectory * dir
Definition: macro.C:5
const IFitVar * xVar
Definition: FitUtilities.h:50
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Binning * xBins
Definition: FitUtilities.h:51
Binning * yBins
Definition: FitUtilities.h:53
Prevent histograms being added to the current directory.
Definition: UtilsExt.h:46
Float_t e
Definition: plot.C:35
std::vector< const IFitVar * > fitVars
Definition: FitUtilities.h:55
A class for generating a covariance matrices as a function of oscillation parameters and systematics ...
const std::vector< std::string > kDetectorSystNames
Definition: SystUtilities.h:19
surf
Definition: demo5.py:35