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 NuS contours for 2018 analysis
9  */
10 
11 #ifdef __CINT__
12 
13 #include <iostream>
14 void MakeSurface(std::string polarity_in, std::string type_in,
15  std::string systs_in, std::string detectors_in, double dm41 = -1,
16  std::string option)
17 {
18  std::cout << "This macro must be run in compiled mode." << std::endl;
19 }
20 
21 #else
22 
28 #include "CAFAna/Fit/Fit.h"
30 
31 using namespace ana;
32 
33 std::vector<double> LogBins(int nbins,double xlo, double xhi)
34 {
35  double log_spacing = exp( (log(xhi) - log(xlo))/(nbins) );
36  vector<double> binning;
37  for (int i = 0; i <= nbins; ++i) {
38  if (i == 0) binning.push_back(xlo);
39  else binning.push_back(log_spacing*binning[i-1]);
40  }
41  return binning;
42 }
43 
44 void MakeSurface(std::string polarity_in, std::string type_in,
45  std::string systs_in, std::string detectors_in, int universe = 1,
46  double dm41 = -1, std::string option="")
47 {
48  DontAddDirectory guard;
49 
50  bool rhc;
51  if (polarity_in == "fhc") rhc = false;
52  else if (polarity_in == "rhc") rhc = true;
53  else throw std::runtime_error("The polarity_in string must be set to \"fhc\" or \"rhc\".");
54 
55  // Check whether we're using systematics
56  bool use_systs;
57  if (systs_in == "stats") use_systs = false;
58  else if (systs_in == "systs") use_systs = true;
59  else throw std::runtime_error("The systs_in string must be \"stats\" or \"systs\".");
60 
61  // Check what type of plot we're making
62  bool th24vsth34 = false;
63  bool th24vsdm41 = false;
64  bool th34vsdm41 = false;
65  if (type_in == "th24vsth34") th24vsth34 = true;
66  else if (type_in == "th24vsdm41") th24vsdm41 = true;
67  else if (type_in == "th34vsdm41") th34vsdm41 = true;
68  else throw std::runtime_error("You must set the type_in string as \"th24vsth34\", \"th24vsdm41\" or \"th34vsdm41\".");
69 
70  // Check on detectors
71  bool use_nd;
72  bool use_fd;
73  if (detectors_in == "nd") {
74  use_nd = true;
75  use_fd = false;
76  }
77  else if (detectors_in == "fd") {
78  use_nd = false;
79  use_fd = true;
80  }
81  else if (detectors_in == "both") {
82  use_nd = true;
83  use_fd = true;
84  }
85  else throw std::runtime_error("The detectors_in string must be \"nd\", \"fd\" or \"both\".");
86 
87  // Make sure dm41 is set for th24 vs th34 case
88  if (th24vsth34 && dm41 == -1)
89  throw std::runtime_error("If you want to run a th24 vs th34 surface, you must specify a value for dm41.");
90 
91  // Check any extra options
92  bool fake_signal = false;
93  bool old_ppfx = false;
94  if (option.size() > 0) {
95  if (option == "fake_signal") fake_signal = true;
96  else if (option == "old_ppfx") old_ppfx = true;
97  else throw std::runtime_error("Surface option \""+option+"\" was not understood.");
98  }
99 
100  // Set up grid environment
101  char* job_c = std::getenv("PROCESS");
102  int Ix, Iy, N;
103  if (job_c) {
104  int job = std::stoi(job_c);
105  Ix = (10 * std::floor((double)job/5.)) + 1;
106  Iy = (10 * (job%5)) + 1;
107  N = 10;
108  }
109  else {
110  std::cout << "PROCESS variable not set, running local test job." << std::endl;
111  Ix = 1;
112  Iy = 1;
113  N = 1;
114  }
115 
116  // Polarity flag & exposures
117  std::string polarity = rhc? "RHC" : "FHC";
119  double fd_pot = rhc ? kAna2018RHCPOT : kAna2018FHCPOT;
120  double fd_livetime = rhc ? kAna2018RHCLivetime : kAna2018FHCLivetime;
121 
122  // Now put it all together
123  PredictionConcat* sim = GetDefaultSimulation(rhc, detectors_in);
124  double pot = use_fd? fd_pot : nd_pot;
125  double livetime = use_fd? fd_livetime : 0;
126 
127  // Merged cosmic spectrum
128  std::vector<TH1D*> cosmics = GetCosmics(rhc, detectors_in);
129  Spectrum cosmic(sim->MergeHists(cosmics), pot, livetime);
130 
131  // Create sterile oscillation calculator
133  calc->SetNFlavors(4);
134  SetAngles(calc);
135 
136  // Load fake data spectrum and add in cosmics
137  std::string fakedata_filename = rhc? "rhc" : "fhc";
138  fakedata_filename += "_fake_data";
139  Spectrum fake_data(LoadFakeData(fakedata_filename, universe, pot, livetime));
140 
141  // No point adding cosmics if there's no far detector
142  if (use_fd) fake_data += cosmic;
143 
144  if (sim->GetNBins() != (unsigned int)fake_data.ToTH1(1)->GetNbinsX())
145  throw std::runtime_error("Number of bins in prediction does not match number of bins in fake data!");
146 
147  ana::SystShifts shifts;
148 
149  // Load the covariance matrix generator, make experiment and Gaussian constraints
150  CovarianceMatrix* mx = GetFullCovMx(sim, calc, pot);
151  OscCovMxExperiment expt(sim, mx, fake_data, cosmic);
152 
153  // Set Gaussian constraints, make MultiExperiment
154  std::vector<const IExperiment*> expts = GetNus18Constraints();
155  expts.push_back(&expt);
156  MultiExperiment multiExpt(expts);
157 
158  // Set seed values and marg vars
159  std::vector<double> fSeedValues = GetNus18SeedValues(calc);
160  std::vector<const IFitVar*> margVars = GetNus18MargVars();
161  if (th24vsdm41) {
162  fSeedValues.push_back(kFitSinSqTheta34Sterile.GetValue(calc));
163  margVars.push_back(&kFitSinSqTheta34Sterile);
164  }
165  else if (th34vsdm41) {
166  fSeedValues.push_back(kFitSinSqTheta24Sterile.GetValue(calc));
167  margVars.push_back(&kFitSinSqTheta24Sterile);
168  }
169 
170  // Set fixed dm41 for th24 vs th34 case
171  if (th24vsth34) calc->SetDm(4, dm41);
172 
173  TH2D* fMargHists[margVars.size()];
174 
175  // Instantiate surfaces depending on type (as they all have different names & binning)
176  TH2D* surf = nullptr;
177  if (th24vsth34) {
178  surf = new TH2D("th34 vs th24",";#theta_{34};#theta_{24}",50,0.,50.,50,0.,45.);
179  for (unsigned int i = 0; i < margVars.size(); ++i)
180  fMargHists[i] = new TH2D("",margVars[i]->ShortName().c_str(),50,0.,50.,50,0.,45.);
181  }
182  else if (th24vsdm41) {
183  std::vector<double> th24Bins = LogBins(50, 1.0e-6, 1);
184  std::vector<double> dm2Bins = LogBins(50, 1.0e-2, 100);
185  surf = new TH2D("th24 vs dm41", ";sin^{2}#theta_{24};#Deltam^{2}_{41}",
186  th24Bins.size()-1, &th24Bins[0], dm2Bins.size()-1, &dm2Bins[0]);
187  for (unsigned int i = 0; i < margVars.size(); ++i)
188  fMargHists[i] = new TH2D("", margVars[i]->ShortName().c_str(),
189  th24Bins.size()-1, &th24Bins[0], dm2Bins.size()-1, &dm2Bins[0]);
190  }
191  else if (th34vsdm41) {
192  std::vector<double> th34Bins = LogBins(50, 1.0e-6, 1);
193  std::vector<double> dm2Bins = LogBins(50, 1.0e-2, 100);
194  surf = new TH2D("th34 vs dm41", ";sin^{2}#theta_{34};#Deltam^{2}_{41}",
195  th34Bins.size()-1, &th34Bins[0], dm2Bins.size()-1, &dm2Bins[0]);
196  for(unsigned int i = 0; i < margVars.size(); ++i)
197  fMargHists[i] = new TH2D("",margVars[i]->ShortName().c_str(),
198  th34Bins.size()-1, &th34Bins[0], dm2Bins.size()-1, &dm2Bins[0]);
199  }
200 
201  for(int i=Ix; i<Ix+N; ++i) {
202  for(int j=Iy; j<Iy+N; ++j) {
203 
204  std::cout << " *********************************************" << std::endl;
205  std::cout << " Bin X Index: " << i
206  << " Bin Y Index: " << j
207  << " Bin Center X: " << surf->GetXaxis()->GetBinCenter(i) << " "
208  << " Bin Center Y: " << surf->GetYaxis()->GetBinCenter(j) << std::endl;
209 
210  double xvar=surf->GetXaxis()->GetBinCenter(i);
211  double yvar=surf->GetYaxis()->GetBinCenter(j);
212 
213  // Set binning variables according to surface type
214  if (th24vsth34) {
217  }
218  else if (th24vsdm41) {
219  kFitSinSqTheta24Sterile.SetValue(calc, xvar);
220  kFitDmSq41Sterile.SetValue(calc, yvar);
221  }
222  else if (th34vsdm41) {
223  kFitSinSqTheta34Sterile.SetValue(calc, xvar);
224  kFitDmSq41Sterile.SetValue(calc, yvar);
225  }
226 
227  for(unsigned int k = 0; k<margVars.size(); ++k) {
228  margVars[k]->SetValue(calc, fSeedValues[k]);
229  std::cout << " Marg values : " << fSeedValues[k] << std::endl;
230  }
231 
232  double chi;
233 
234  if(margVars.size()==0){
235  chi = multiExpt.ChiSq(calc);
236  }
237  else {
238  Fitter fitter(&multiExpt, margVars, {});
239  SystShifts systSeed = SystShifts::Nominal();
240  chi = fitter.Fit(calc, systSeed, Fitter::kQuiet)->EvalMetricVal();
241 
242  for(unsigned int k = 0; k < margVars.size(); ++k){
243  fMargHists[k]->SetBinContent(i, j, margVars[k]->GetValue(calc));
244  std::cout << " margVar Value: " << margVars[k]->GetValue(calc) << " for " << k << std::endl;
245  }
246  }
247 
248  std::cout << " chi from first fit: " << chi << std::endl;
249  surf->SetBinContent(i,j, chi);
250  }
251  }
252 
253  // Write surface to file
254  std::string output_filename = type_in + "_" + polarity_in + "_" + systs_in + "_" + detectors_in + ".root";
255  TFile* f = new TFile(output_filename.c_str(), "recreate");
256  surf->Write("surface");
257 
258  // Write marginalisation pull surfaces to file
259  for(unsigned int i=0;i<margVars.size();++i) {
260  fMargHists[i]->Write(margVars[i]->ShortName().c_str());
261  }
262 
263  // Bye o/
264  f->Close();
265 
266 }
267 
268 #endif
std::vector< double > GetNus18SeedValues(osc::IOscCalcAdjustable *calc)
Seed values for Nus18 marginalisation parameters.
const FitSinSqTheta24Sterile kFitSinSqTheta24Sterile
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void MakeSurface()
Definition: MakeSurface.C:31
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
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
std::vector< double > LogBins(int nbins, double xlo, double xhi)
Definition: MakeSurface.C:33
Adapt the PMNS_Sterile calculator to standard interface.
const double kAna2018SensitivityFHCNDPOT
Definition: Exposures.h:210
const FitTheta24InDegreesSterile kFitTheta24InDegreesSterile
std::vector< const IExperiment * > GetNus18Constraints()
Gaussian constrains on atmospheric parameters.
const FitSinSqTheta34Sterile kFitSinSqTheta34Sterile
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
void SetAngles(osc::OscCalcSterile *calc)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const double kAna2018RHCPOT
Definition: Exposures.h:208
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
const int nbins
Definition: cellShifts.C:15
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.
std::string getenv(std::string const &name)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
#define pot
Compare a single data spectrum to the MC + cosmics expectation.
virtual double ChiSq(osc::IOscCalcAdjustable *osc, const SystShifts &syst=SystShifts::Nominal()) const override
const double j
Definition: BetheBloch.cxx:29
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Example routine for calculating the "ultimate e-m mother" of a particle in a simulated event...
Definition: FillTruth.h:16
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
double livetime
Definition: saveFDMCHists.C:21
fvar< T > floor(const fvar< T > &x)
Definition: floor.hpp:11
const double kAna2018SensitivityRHCNDPOT
Definition: Exposures.h:211
double GetValue(const osc::IOscCalcAdjustable *osc) const override
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::vector< TH1D * > GetCosmics(bool rhc, std::string detector)
Function to get default cosmics.
const FitTheta34InDegreesSterile kFitTheta34InDegreesSterile
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
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
const FitDmSq41Sterile kFitDmSq41Sterile
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const double kAna2018RHCLivetime
Definition: Exposures.h:214
surf
Definition: demo5.py:35
enum BeamMode string