make_fc_mh_nersc_2019.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
4 
5 #include "CAFAna/Fit/Fit.h"
7 
10 
11 #include "CAFAna/Vars/FitVars.h"
12 
13 #include "CAFAna/FC/FCPoint.h"
14 #include "CAFAna/FC/FCCollection.h"
15 #include "CAFAna/FC/FCSurface.h"
16 
20 
21 //#include "3FlavorAna/Ana2018/nue/joint_fit_2018_tools.h"
23 
24 #include "CAFAna/Analysis/Calcs.h"
25 #include "OscLib/IOscCalc.h"
26 
27 #include "TRandom3.h"
28 #include "TH1.h"
29 #include "TGraph.h"
30 
31 using namespace ana;
32 
33 void make_fc_mh_nersc_2019(bool nh, int NPts, int N, int jid, int npts_flush = 0)
34 {
35  std::string hier = nh ? "nh" : "ih";
36  std::string number = std::to_string(N);
37  std::string fcout = "FCCol_mass_hierarchy_number_" + hier + "_v"
38  + number
39  + "_"
40  + std::to_string(jid)
41  +".root";
42 
43  FCCollection fccol;
45 
46  double thisPOT;
47  double thisLivetime;
48 
49  // Grab the systs
50  // bools to make the function calls a little more readable
51  bool nueExclusive = false;
52  bool numuExclusive = false;
53  bool isFHC = true;
54  bool isRHC = true;
55 
56  // joint
57  std::vector<const ISyst*> systs = GetJointFitSystematicList(true, nueExclusive, numuExclusive, true, true, true);
58 
59  // nue fhc
60  std::vector<const ISyst*> snue_fhc = getAllAna2018Systs(kNueAna2018, true, kFHC);
61 
62  // nue rhc
63  std::vector<const ISyst*> snue_rhc = getAllAna2018Systs(kNueAna2018, true, kRHC);
64 
65  // numu
66  std::vector<const ISyst*> snumu = GetJointFitSystematicList(true, nueExclusive, !numuExclusive, isFHC, isRHC, true);
67 
68  // We're setting up the requisite inputs to the nue experiment
69  // nue prediction
70 
71  std::vector <const IPrediction*> nue_preds;
72  std::vector <std::pair <TH1D*, double > > nue_cosmics;
73  std::vector <Spectrum> nue_cosmics_spectra;
74  bool minimizeMemory = true;
75 
76  // Grab the nue predictions for fhc and rhc
77  nue_preds.push_back(GetNuePrediction2019("combo", calc, true, "fhc", false, false, minimizeMemory, true));
78  nue_preds.push_back(GetNuePrediction2019("prop", calc, true, "rhc", false, false, minimizeMemory, true));
79 
80  std::cout << "\nDone getting Nue Predictions\n\n" ;
81 
82  nue_cosmics.push_back(GetNueCosmics2019("fhc", false, true));
83  nue_cosmics.push_back(GetNueCosmics2019("rhc", false, true));
84 
85  // We need to store cosmic hists as spectra to add into our fake data
86  // FHC and RHC have different POT and Livetimes
87  for(int i = 0; i < 2; i++) {
88  if(i==0) {
89  thisPOT = kAna2019FHCPOT;
90  thisLivetime = kAna2019FHCLivetime;
91  }
92  else {
93  thisPOT = kAna2019RHCPOT;
94  thisLivetime = kAna2019RHCLivetime;
95  }
96  nue_cosmics_spectra.push_back(Spectrum(nue_cosmics[i].first,
97  thisPOT, thisLivetime));
98  }
99 
100 
101  //////////////////////////////////////////////////////////////////////
102  // Numu Prediction
103  //////////////////////////////////////////////////////////////////////
104  std::vector<const IPrediction*> numu_preds;
105  std::vector< std::pair<TH1D*, double> > numu_cosmics;
106  std::vector<Spectrum> numu_cosmics_spectra;
107 
108  std::vector<const IPrediction*> this_numu_preds;
109  std::vector< std::pair<TH1D*, double> > this_numu_cosmics;
110 
111  // Construct numu inputs
112  for(std::string beam : {"fhc", "rhc"}){
113  this_numu_preds = GetNumuPredictions2019(4,true, beam, false, kNuMu, minimizeMemory, true);
114  this_numu_cosmics = GetNumuCosmics2019(4, beam, false, true);
115  numu_preds.insert(numu_preds.end(),
116  this_numu_preds.begin(),
117  this_numu_preds.end());
118  numu_cosmics.insert(numu_cosmics.end(),
119  this_numu_cosmics.begin(),
120  this_numu_cosmics.end());
121 
122  for(std::pair<TH1D*, double> h : this_numu_cosmics){
123  thisPOT = kAna2019FHCPOT;
124  thisLivetime = kAna2019FHCLivetime;
125  if(beam == "rhc") {
126  thisPOT = kAna2019RHCPOT;
127  thisLivetime = kAna2019RHCLivetime;
128  }
129  numu_cosmics_spectra.push_back(Spectrum(h.first, thisPOT, thisLivetime));
130  }
131 
132  }
133 
134  // Grab analysis slices to find the bin centers and the best
135  // fit point
136  std::string helperupsname = std::string(getenv("FCHELPERANA2019_LIB_PATH"));
137  // incorrectly named 2017, but don't worry
138  std::string anaslicename = helperupsname +
139  "/Surfaces/hist_slices_2019_joint_realData_both_systs_dmsq_noOct.root";
140 
141  TFile* anaslicefile = TFile::Open(anaslicename.c_str());
142 
143  std::string HIER = nh ? "NH" : "IH";
144 
145  std::string slice_dmsq_hier = "slice_dmsq_" + HIER;
146  HIER+="_";
147 
148  TH1* anaslice = (TH1F*) anaslicefile->Get(slice_dmsq_hier.c_str());
149 
150  // best fit bin
151  // get upper and lower edge to get a range of dmsq
152  // to throw mock experiments in
153  int best_bin = anaslice->GetMinimumBin();
154  double true_dmsq = anaslice->GetBinCenter(best_bin);
155  double best_bin_low = anaslice->GetBinLowEdge(best_bin);
156  double best_bin_high = anaslice->GetBinLowEdge(best_bin+1);
157  double widthX = anaslice->GetBinWidth(best_bin);
158  std::cout << "Got best fit point bin info" << endl;
159 
160  std::string helpername =
161  helperupsname
162  + "/FCInputs/"
163  + "/slices_FCInput_2019_joint_realData_both_systs_dmsq_noOct.root";
164 
165  TFile* fchelp = TFile::Open(helpername.c_str());
166 
167  // Get profiled parameters at best fit
168  TGraph *hdelta = (TGraph*)fchelp->Get((HIER + "DeltaCPpi").c_str());
169  TGraph *hssth23 = (TGraph*)fchelp->Get((HIER + "SinSqTheta23").c_str());
170  TGraph *hss2th13 = (TGraph*)fchelp->Get((HIER + "SinSq2Theta13").c_str());
171 
172  double seeddelta = hdelta->Eval(true_dmsq);
173  double seedssth23 = hssth23->Eval(true_dmsq);
174  double seedss2th13 = hss2th13->Eval(true_dmsq);
175 
176  std::cout << "Got all of the parameters to make mock data" << endl;
177 
178  delete hdelta;
179  delete hssth23;
180  delete hss2th13;
181 
182  // Now, look for all the systematics
183  std::map<std::string,double> seedsyst;
184  for(const ISyst* syst :systs) {
185  TGraph* h = (TGraph*) fchelp->Get((HIER + syst->ShortName()).c_str());
186  if(!h){
187  std::cout << "Don't see a prof hist for " << syst->ShortName() << ". Continuing, but check your ups version." << std::endl;
188  continue;
189  }
190  double seedval = h->Eval(true_dmsq);
191 
192  seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
193  }
194 
195 
196  // Need random numbers to sample the bin
197  TRandom3 rnd(0);
198 
199  // Now everything is set up, we can start throwing FC mock experiments
200 
201  for(int i = 0; i < NPts; i++){
202  ResetOscCalcToDefault(calc);
203 
204 
205  // The variable you're plotting isn't a nuisance parameter, should
206  // be set to the trueX that you threw for this experiment.
207  // Other oscillation parameters are nuisance parameters and
208  // should be seeded at their best fit
209  kFitDeltaInPiUnits.SetValue(calc, seeddelta);
210  kFitSinSqTheta23.SetValue(calc, seedssth23);
211  kFitSinSq2Theta13.SetValue(calc, seedss2th13);
212 
213  nh ?
214  kFitDmSq32NormalHierarchy.SetValue(calc, true_dmsq) :
215  kFitDmSq32InvertedHierarchy.SetValue(calc, true_dmsq);
216 
217  // Set syst shift - they're nuisance parameters, so set them
218  // to the best fit. Nue and Numu predictions are taught about
219  // slightly different sets of systematics, so we need to throw
220  // the predictions with different seeds
221  SystShifts seednue_fhc;
222  for(const ISyst *syst : snue_fhc){
223  if(seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
224  seednue_fhc.SetShift(syst, seedsyst[syst->ShortName()]);
225  }
226  SystShifts seednue_rhc;
227  for(const ISyst *syst : snue_rhc){
228  if(seedsyst.find(syst->ShortName()) ==seedsyst.end()) continue;
229  seednue_rhc.SetShift(syst, seedsyst[syst->ShortName()]);
230  }
231  SystShifts seednumu;
232  for(const ISyst *syst : snumu){
233  if(seedsyst.find(syst->ShortName()) ==seedsyst.end()) continue;
234  seednumu.SetShift(syst, seedsyst[syst->ShortName()]);
235  }
236 
237  // Make a mock data spectrum for nue
238  // use mock data + cosmic as predictions
239  double mock_integral_fhc = 0;
240  double mock_integral_rhc = 0;
241  assert(nue_preds.size()==2 && "Something is wrong with nue preds. Size !=2");
242  std::vector<Spectrum> mocknue;
243  SystShifts this_seednue;
244  for (int j = 0; j < (int) nue_preds.size(); j++){
245  std::cout << "\n\n Lets make some nue fake data \n\n";
246  double thisPOT = j == 0? kAna2019FHCPOT : kAna2019RHCPOT;
247  std::cout << "Getting fake data for j=" << std::to_string(j);
248  if(j==0) this_seednue = seednue_fhc;
249  else this_seednue = seednue_rhc;
250  Spectrum fakenue = nue_preds[j]->PredictSyst(calc,this_seednue);
251  std::cout << "\n\n Thats some nice nue fake data \n\n";
252  // Add in cosmic prediction here, make sure it gets Poisson fluctuated
253  fakenue += nue_cosmics_spectra[j];
254  mocknue.push_back(fakenue.MockData(thisPOT));
255  // Save this, for posterity's sake
256  if(j==0) mock_integral_fhc += mocknue[j].Integral(thisPOT);
257  if(j==1) mock_integral_rhc += mocknue[j].Integral(thisPOT);
258  }
259 
260  std::cout << "\n\n Got Nue Preds \n\n";
261  std::vector<Spectrum> mocknumu;
262  // Make a mock data spectrum for numu
263  assert(numu_cosmics_spectra.size()==8 && "Issue with N quantiles in cosmics");
264  assert(numu_preds.size()==8 && "Issue with N quantiles in predictions");
265  for (int j = 0; j < 8; j++){
266  double thisPOT = kAna2019FHCPOT;
267  if(j > 3) thisPOT = kAna2019RHCPOT;
268  Spectrum fakenm = numu_preds[j]->PredictSyst(calc,seednumu);
269  fakenm += numu_cosmics_spectra[j]; // Again, add cosmic spectra pre-fluctuation
270  mocknumu.push_back(fakenm.MockData(thisPOT)); // Fluctuate
271  }
272  std::cout << "Got Numu Preds" << endl;
273 
274 
275  // Set up experiment, nova nue + nova numu (4 quantiles) + PDG reactor
276  std::vector<const IExperiment*> expts;
277  std::vector<const IExperiment*> numu_expts;
278  // Nue experiments
279  for(int j = 0; j < 2; j++) {
280  expts.push_back(new SingleSampleExperiment(nue_preds[j], mocknue[j],
281  nue_cosmics[j].first,
282  nue_cosmics[j].second)); // Nue
283  }
284  for (int j = 0; j < 8; j++) {
285  expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
286  numu_cosmics[j].first,
287  numu_cosmics[j].second));
288  numu_expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
289  numu_cosmics[j].first,
290  numu_cosmics[j].second));
291 
292  }
293 
294  // Add reactor
295  expts.push_back( WorldReactorConstraint2017() );
296  numu_expts.push_back(WorldReactorConstraint2017());
297  MultiExperiment expt(expts);
298  MultiExperiment numu_expt(numu_expts);
299 
300  // Correlate your systs
301  // A bit tricky since nue and numu predictions know about a
302  // different list of systs, but there are functions to do this
303  // in joint_fit_2018_tools.h
304  expt.SetSystCorrelations(0, GetCorrelations(true,true));
305  expt.SetSystCorrelations(1, GetCorrelations(true, false));
306 
307  auto notfornumu = GetCorrelations(false, false);
308 
309  for(int j = 0; j < 8; j++) {
310  expt.SetSystCorrelations(j+2, notfornumu);
311  }
312 
313  // Finally, we get to do some FC
314  //////////////////////////////////////////////////////////////////
315  // Use numu only fit to find th23 and dmsq seeds for the joint fit
316  //////////////////////////////////////////////////////////////////
317 
319  SystShifts auxShifts = SystShifts::Nominal();
320  double maxmixing = 0.514;
321 
322  Fitter fitnumu_only_ih(&numu_expt, {&kFitSinSqTheta23, &kFitDmSq32InvertedHierarchy});
323  Fitter fitnumu_only_nh(&numu_expt, {&kFitSinSqTheta23, &kFitDmSq32NormalHierarchy});
324 
325  std::cout << "\nFinding the seeds from numu only fit" << endl;
326  double minchi_numu = nh ?
327  fitnumu_only_nh.Fit(numu_calc, auxShifts,
328  {{&kFitDmSq32NormalHierarchy, {-2.35e-3}},
329  {&kFitSinSqTheta23, {0.4}}},
330  {},
331  Fitter::kQuiet):
332  fitnumu_only_ih.Fit(numu_calc, auxShifts,
333  {{&kFitDmSq32InvertedHierarchy, {-2.35e-3}},
334  {&kFitSinSqTheta23, {0.4}}},
335  {},
336  Fitter::kQuiet);
337 
338  double pre_seed_th23 = kFitSinSqTheta23.GetValue(numu_calc);
339  double pre_seed_dmsq = kFitDmSq32.GetValue(numu_calc);
340 
341  double pre_seed_th23_LO = ( pre_seed_th23 > maxmixing ) ?
342  (2*maxmixing - pre_seed_th23) : pre_seed_th23;
343  double pre_seed_th23_UO = ( pre_seed_th23 > maxmixing ) ?
344  pre_seed_th23 : (2*maxmixing-pre_seed_th23);
345  if(pre_seed_th23_LO > 0.5) pre_seed_th23_LO = 0.45;
346 
347  // done getting pre seeds
348  //-----------------------------------------------------------------
349 
350  // Set up fitter for chi best
351  // make a vector of oscillation parameters that are our
352  // nuisance parameters
353  std::vector<const IFitVar*> ih_fitVars = {&kFitDeltaInPiUnits,
357 
358  Fitter ih_fit(&expt, ih_fitVars, systs);
359  SystShifts seedShifts = {};
360 
361  double nh_delta = 0;
362  double nh_ssth23 = 0;
363  double nh_th13 = 0;
364  double nh_dmsq = 0;
365  double ih_delta = 0;
366  double ih_ssth23 = 0;
367  double ih_th13 = 0;
368  double ih_dmsq = 0;
369 
370  std::vector <double> ih_delta_seeds = {0., 0.5, 1.0, 1.5};
371  std::vector <double> ih_ssth23_seeds = {pre_seed_th23_LO, pre_seed_th23_UO, maxmixing};
372  std::map <const IFitVar*, std::vector<double>> ih_seedPts =
373  {{&kFitSinSqTheta23, ih_ssth23_seeds},
374  {&kFitDeltaInPiUnits, ih_delta_seeds},
375  {&kFitDmSq32InvertedHierarchy, {-1*fabs(pre_seed_dmsq)}}};
376  std::cout << "Starting IH fit ----->" << endl;
377  auxShifts.ResetToNominal();
378  double ih_chi = ih_fit.Fit(calc,
379  auxShifts,
380  ih_seedPts,
381  {},
382  Fitter::kQuiet);
384  ih_delta = kFitDeltaInPiUnits.GetValue(calc);
385  ih_ssth23 = kFitSinSqTheta23.GetValue(calc);
386  ih_th13 = kFitSinSq2Theta13.GetValue(calc);
387 
388  std::vector <double> nh_delta_seeds = {0., 0.5, 1.0, 1.5};
389  std::vector <double> nh_ssth23_seeds = {pre_seed_th23_LO, pre_seed_th23_UO, maxmixing};
390  std::map <const IFitVar*, std::vector<double>> nh_seedPts =
391  {{&kFitSinSqTheta23, nh_ssth23_seeds},
392  {&kFitDeltaInPiUnits, nh_delta_seeds},
393  {&kFitDmSq32NormalHierarchy, {fabs(pre_seed_dmsq)}},
394  {&kFitSinSq2Theta13, {0.082}}};
395 
396  std::vector <const IFitVar*> nh_vars =
398  Fitter nh_fit(&expt, nh_vars, systs);
399 
400  std::cout << "Starting the NH fit ----->" << endl;
401  auxShifts.ResetToNominal();
402  double nh_chi = nh_fit.Fit(calc,
403  auxShifts,
404  nh_seedPts,
405  {},
406  Fitter::kQuiet);
407  nh_dmsq = kFitDmSq32.GetValue(calc);
408  nh_ssth23 = kFitSinSqTheta23.GetValue(calc);
409  nh_delta = kFitDeltaInPiUnits.GetValue(calc);
410  nh_th13 = kFitSinSq2Theta13.GetValue(calc);
411 
412  // Keep same order regardless of hierarchy:
413  // true = IH
414  // best = NH
415  FCPoint pt(seeddelta, seedssth23, seedss2th13, true_dmsq,
416  pre_seed_th23, pre_seed_dmsq,
417  ih_delta, ih_ssth23, ih_th13, ih_dmsq,
418  ih_chi,
419  nh_delta, nh_ssth23, nh_th13, nh_dmsq,
420  nh_chi);
421 
422  fccol.AddPoint(pt);
423  unsigned int nfc_pts = fccol.NPoints();
424 
425  if(npts_flush != 0){
426  if(nfc_pts % npts_flush == 0){
427  std::cout << "Saving " + std::to_string(nfc_pts) +
428  " points to " << fcout << endl;
429  fccol.SaveToFile(fcout);
430  }
431  }
432  } // end loop over N points
433  fccol.SaveToFile(fcout);
434 }
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
Definition: Calcs.cxx:23
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double GetValue(const osc::IOscCalcAdjustable *osc) const override
const double kAna2019RHCLivetime
Definition: Exposures.h:227
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var&#39;s GetValue()
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:250
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
std::vector< const IPrediction * > GetNumuPredictions2019(const int nq=4, bool useSysts=true, std::string beam="fhc", bool GetFromUPS=false, ENu2018ExtrapType numuExtrap=kNuMu, bool minimizeMemory=false, bool NERSC=false)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:171
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:48
static SystShifts Nominal()
Definition: SystShifts.h:34
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:42
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
const IPrediction * GetNuePrediction2019(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:300
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
Represents the results of a single Feldman-Cousins pseudo-experiment.
Definition: FCPoint.h:8
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
expt
Definition: demo5.py:34
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::string getenv(std::string const &name)
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const double j
Definition: BetheBloch.cxx:29
const double kAna2019FHCLivetime
Definition: Exposures.h:226
std::vector< float > Spectrum
Definition: Constants.h:610
static bool isFHC
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:214
OStream cout
Definition: OStream.cxx:6
const double kAna2019RHCPOT
Definition: Exposures.h:224
Combine multiple component experiments.
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2019(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
const double kAna2019FHCPOT
Definition: Exposures.h:223
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
const FitDmSq32InvertedHierarchy kFitDmSq32InvertedHierarchy
Definition: FitVars.cxx:23
void make_fc_mh_nersc_2019(bool nh, int NPts, int N, int jid, int npts_flush=0)
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:244
const FitDmSq32NormalHierarchy kFitDmSq32NormalHierarchy
Definition: FitVars.cxx:20
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool isFHC=true, bool isRHC=true, bool intersection=true, bool ptExtrap=true)
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:80
void SetSystCorrelations(int idx, const std::vector< std::pair< const ISyst *, const ISyst * >> &corrs)
std::pair< TH1D *, double > GetNueCosmics2019(std::string beam, bool GetFromUPS=false, bool NERSC=false)
Compare a single data spectrum to the MC + cosmics expectation.
Collection of FCPoint. Serializable to/from a file.
Definition: FCCollection.h:12
enum BeamMode string