make_fc_oct_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 
34 void make_fc_oct_nersc_2019(int NPts, int N, int jid, int npts_flush = 0)
35 {
36  std::string number = std::to_string(N);
37  std::string fcout = "FCCol_octant_number_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  nue_cosmics.push_back(GetNueCosmics2019("fhc", false, true));
81  nue_cosmics.push_back(GetNueCosmics2019("rhc", false, true));
82 
83  // We need to store cosmic hists as spectra to add into our fake data
84  // FHC and RHC have different POT and Livetimes
85  for(int i = 0; i < 2; i++) {
86  if(i==0) {
87  thisPOT = kAna2019FHCPOT;
88  thisLivetime = kAna2019FHCLivetime;
89  }
90  else {
91  thisPOT = kAna2019RHCPOT;
92  thisLivetime = kAna2019RHCLivetime;
93  }
94  nue_cosmics_spectra.push_back(Spectrum(nue_cosmics[i].first,
95  thisPOT, thisLivetime));
96  }
97 
98 
99  //////////////////////////////////////////////////////////////////////
100  // Numu Prediction
101  //////////////////////////////////////////////////////////////////////
102  std::vector<const IPrediction*> numu_preds;
103  std::vector< std::pair<TH1D*, double> > numu_cosmics;
104  std::vector<Spectrum> numu_cosmics_spectra;
105 
106  std::vector<const IPrediction*> this_numu_preds;
107  std::vector< std::pair<TH1D*, double> > this_numu_cosmics;
108 
109  // Construct numu inputs
110  for(std::string beam : {"fhc", "rhc"}){
111  this_numu_preds = GetNumuPredictions2019(4,true, beam, false, kNuMu, minimizeMemory, true);
112  this_numu_cosmics = GetNumuCosmics2019(4, beam, false, true);
113  numu_preds.insert(numu_preds.end(),
114  this_numu_preds.begin(),
115  this_numu_preds.end());
116  numu_cosmics.insert(numu_cosmics.end(),
117  this_numu_cosmics.begin(),
118  this_numu_cosmics.end());
119 
120  for(std::pair<TH1D*, double> h : this_numu_cosmics){
121  thisPOT = kAna2019FHCPOT;
122  thisLivetime = kAna2019FHCLivetime;
123  if(beam == "rhc") {
124  thisPOT = kAna2019RHCPOT;
125  thisLivetime = kAna2019RHCLivetime;
126  }
127  numu_cosmics_spectra.push_back(Spectrum(h.first, thisPOT, thisLivetime));
128  }
129 
130  }
131 
132  // Grab analysis slices to find the bin centers and the best
133  // fit point
134  std::string helperupsname = std::string(getenv("FCHELPERANA2019_LIB_PATH"));
135  std::string anaslicename = helperupsname +
136  "/Surfaces/hist_slices_2019_joint_realData_both_systs_dmsq.root";
137 
138  // need the best fit point in the lower octant to seed pseudoexperiments
139  TFile* anaslicefile = TFile::Open(anaslicename.c_str());
140  TH1* sliceNHLO = (TH1F*) anaslicefile->Get("slice_dmsq_NHLO");
141  TH1* sliceIHLO = (TH1F*) anaslicefile->Get("slice_dmsq_IHLO");
142  int bestNHLO = sliceNHLO->GetMinimumBin();
143  int bestIHLO = sliceIHLO->GetMinimumBin();
144 
145  int best_bin = bestNHLO;
146  std::string option = "NHLO";
147  double true_dmsq = sliceNHLO->GetBinCenter(best_bin);
148  if( sliceNHLO->GetBinContent(bestNHLO) > sliceIHLO->GetBinContent(bestIHLO)){
149  best_bin = bestIHLO;
150  option = "IHLO";
151  true_dmsq = sliceIHLO->GetBinCenter(best_bin);
152  }
153 
154  // Get profiled parameters at best fit
155  std::string helpername =
156  helperupsname
157  + "/FCInputs/"
158  + "slices_FCInput_2019_joint_realData_both_systs_dmsq.root";
159 
160  TFile* fchelp = TFile::Open(helpername.c_str());
161 
162  TGraph* hssth23 = (TGraph*)fchelp->Get((option+"_SinSqTheta23").c_str());
163  TGraph *hdelta = (TGraph*)fchelp->Get((option+"_DeltaCPpi").c_str());
164  TGraph *hss2th13 = (TGraph*)fchelp->Get((option+"_SinSq2Theta13").c_str());
165 
166  double seedssth23 = hssth23->Eval(true_dmsq);
167  double seeddelta = hdelta->Eval(true_dmsq);
168  double seedss2th13 = hss2th13->Eval(true_dmsq);
169  std::cout << "Got all of the parameters to make mock data" << endl;
170 
171  delete sliceNHLO;
172  delete sliceIHLO;
173  delete hssth23;
174  delete hdelta;
175  delete hss2th13;
176 
177  // Now, look for all the systematics
178  std::map<std::string,double> seedsyst;
179  for(const ISyst* syst :systs) {
180  TGraph* h = (TGraph*) fchelp->Get((option+"_"+syst->ShortName()).c_str());
181  if(!h){
182  std::cout << "Don't see a prof hist for " << syst->ShortName() << ". Continuing, but check your ups version." << std::endl;
183  continue;
184  }
185  double seedval = h->Eval(true_dmsq);
186 
187  seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
188  }
189 
190  // Need random numbers to sample the bin
191  TRandom3 rnd(0);
192 
193  // Now everything is set up, we can start throwing FC mock experiments
194 
195  for(int i = 0; i < NPts; i++){
196  ResetOscCalcToDefault(calc);
197 
198 
199  // The variable you're plotting isn't a nuisance parameter, should
200  // be set to the trueX that you threw for this experiment.
201  // Other oscillation parameters are nuisance parameters and
202  // should be seeded at their best fit
203  kFitDeltaInPiUnits.SetValue(calc, seeddelta);
204  kFitSinSq2Theta13.SetValue(calc, seedss2th13);
205  kFitDmSq32.SetValue(calc, true_dmsq);
206  kFitSinSqTheta23LowerOctant->SetValue(calc, seedssth23);
207 
208  // Set syst shift - they're nuisance parameters, so set them
209  // to the best fit. Nue and Numu predictions are taught about
210  // slightly different sets of systematics, so we need to throw
211  // the predictions with different seeds
212  SystShifts seednue_fhc;
213  for(const ISyst *syst : snue_fhc){
214  if(seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
215  seednue_fhc.SetShift(syst, seedsyst[syst->ShortName()]);
216  }
217  SystShifts seednue_rhc;
218  for(const ISyst *syst : snue_rhc){
219  if(seedsyst.find(syst->ShortName()) ==seedsyst.end()) continue;
220  seednue_rhc.SetShift(syst, seedsyst[syst->ShortName()]);
221  }
222  SystShifts seednumu;
223  for(const ISyst *syst : snumu){
224  if(seedsyst.find(syst->ShortName()) ==seedsyst.end()) continue;
225  seednumu.SetShift(syst, seedsyst[syst->ShortName()]);
226  }
227 
228  // Make a mock data spectrum for nue
229  // use mock data + cosmic as predictions
230  double mock_integral_fhc = 0;
231  double mock_integral_rhc = 0;
232  assert(nue_preds.size()==2 && "Something is wrong with nue preds. Size !=2");
233  std::vector<Spectrum> mocknue;
234  SystShifts this_seednue;
235  for (int j = 0; j < (int) nue_preds.size(); j++){
236 
237  double thisPOT = j == 0? kAna2019FHCPOT : kAna2019RHCPOT;
238 
239  if(j==0) this_seednue = seednue_fhc;
240  else this_seednue = seednue_rhc;
241  Spectrum fakenue = nue_preds[j]->PredictSyst(calc,this_seednue);
242 
243  // Add in cosmic prediction here, make sure it gets Poisson fluctuated
244  fakenue += nue_cosmics_spectra[j];
245  mocknue.push_back(fakenue.MockData(thisPOT));
246  // Save this, for posterity's sake
247  if(j==0) mock_integral_fhc += mocknue[j].Integral(thisPOT);
248  if(j==1) mock_integral_rhc += mocknue[j].Integral(thisPOT);
249  }
250 
251 
252  std::vector<Spectrum> mocknumu;
253  // Make a mock data spectrum for numu
254  assert(numu_cosmics_spectra.size()==8 && "Issue with N quantiles in cosmics");
255  assert(numu_preds.size()==8 && "Issue with N quantiles in predictions");
256  for (int j = 0; j < 8; j++){
257  double thisPOT = kAna2019FHCPOT;
258  if(j > 3) thisPOT = kAna2019RHCPOT;
259  Spectrum fakenm = numu_preds[j]->PredictSyst(calc,seednumu);
260  fakenm += numu_cosmics_spectra[j]; // Again, add cosmic spectra pre-fluctuation
261  mocknumu.push_back(fakenm.MockData(thisPOT)); // Fluctuate
262  }
263 
264 
265 
266  // Set up experiment, nova nue + nova numu (4 quantiles) + PDG reactor
267  std::vector<const IExperiment*> expts;
268  std::vector<const IExperiment*> numu_expts;
269  // Nue experiments
270  for(int j = 0; j < 2; j++) {
271  expts.push_back(new SingleSampleExperiment(nue_preds[j], mocknue[j],
272  nue_cosmics[j].first,
273  nue_cosmics[j].second)); // Nue
274  }
275  for (int j = 0; j < 8; j++) {
276  expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
277  numu_cosmics[j].first,
278  numu_cosmics[j].second));
279  numu_expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
280  numu_cosmics[j].first,
281  numu_cosmics[j].second));
282 
283  }
284 
285  // Add reactor
286  expts.push_back( WorldReactorConstraint2017() );
287  numu_expts.push_back(WorldReactorConstraint2017());
288  MultiExperiment expt(expts);
289  MultiExperiment numu_expt(numu_expts);
290 
291  // Correlate your systs
292  // A bit tricky since nue and numu predictions know about a
293  // different list of systs, but there are functions to do this
294  // in joint_fit_2018_tools.h
295  expt.SetSystCorrelations(0, GetCorrelations(true,true));
296  expt.SetSystCorrelations(1, GetCorrelations(true, false));
297 
298  auto notfornumu = GetCorrelations(false, false);
299 
300  for(int j = 0; j < 8; j++) {
301  expt.SetSystCorrelations(j+2, notfornumu);
302  }
303 
304  // Finally, we get to do some FC
305  //////////////////////////////////////////////////////////////////
306  // Use numu only fit to find th23 and dmsq seeds for the joint fit
307  //////////////////////////////////////////////////////////////////
308 
310  SystShifts auxShifts = SystShifts::Nominal();
311  double ssth23mirror = 0.514; // the point about which
312  // the disappearance oscillation is degenerate
313  Fitter fitnumu_only(&numu_expt, {kFitSinSqTheta23LowerOctant, &kFitDmSq32});
314 
315  std::cout << "\nFinding the seeds from numu only fit" << endl;
316  double minchi_numu = fitnumu_only.Fit(numu_calc, auxShifts,
317  {{&kFitDmSq32, {-2.35e-3}},
318  {kFitSinSqTheta23LowerOctant, {0.4}}},
319  {},
320  Fitter::kQuiet);
321  double pre_seed_th23 = kFitSinSqTheta23LowerOctant->GetValue(numu_calc);
322  double pre_seed_dmsq = kFitDmSq32.GetValue(numu_calc);
323 
324  double pre_seed_th23_LO = ( pre_seed_th23 > ssth23mirror ) ?
325  (2*ssth23mirror - pre_seed_th23) : pre_seed_th23;
326  double pre_seed_th23_UO = ( pre_seed_th23 > ssth23mirror ) ?
327  pre_seed_th23 : (2*ssth23mirror-pre_seed_th23);
328  if(pre_seed_th23_LO > 0.5) pre_seed_th23_LO = 0.45;
329 
330  // done getting pre seeds
331  //-----------------------------------------------------------------
332 
333  // Set up fitter for chi best
334  // make a vector of oscillation parameters that are our
335  // nuisance parameters
336  std::vector<const IFitVar*> LO_fitVars = {&kFitDeltaInPiUnits,
339  &kFitDmSq32};
340 
341  Fitter LO_fit(&expt, LO_fitVars, systs);
342  SystShifts seedShifts = {};
343 
344  double UO_delta = 0;
345  double UO_ssth23 = 0;
346  double UO_th13 = 0;
347  double UO_dmsq = 0;
348  double LO_dmsq = 0;
349  double LO_delta = 0;
350  double LO_ssth23 = 0;
351  double LO_th13 = 0;
352 
353  std::vector <double> LO_delta_seeds = {0., 0.5, 1.0, 1.5};
354  std::vector <double> LO_ssth23_seeds = {pre_seed_th23_LO};
355  // if below max mixing, add ssth23mirror to seeds in LO
356  if(ssth23mirror < 0.5) LO_ssth23_seeds.push_back(ssth23mirror);
357 
358  std::map <const IFitVar*, std::vector<double>> LO_seedPts =
359  {{kFitSinSqTheta23LowerOctant, LO_ssth23_seeds},
360  {&kFitDeltaInPiUnits, LO_delta_seeds},
361  {&kFitDmSq32, {-1*fabs(pre_seed_dmsq),fabs(pre_seed_dmsq)}}};
362  std::cout << "Starting LO fit ----->" << endl;
363  auxShifts.ResetToNominal();
364  double LO_chi = LO_fit.Fit(calc,
365  auxShifts,
366  LO_seedPts,
367  {},
368  Fitter::kQuiet);
369  LO_dmsq = kFitDmSq32.GetValue(calc);
370  LO_delta = kFitDeltaInPiUnits.GetValue(calc);
371  LO_ssth23 = kFitSinSqTheta23LowerOctant->GetValue(calc);
372  LO_th13 = kFitSinSq2Theta13.GetValue(calc);
373 
374  std::vector <double> UO_delta_seeds = {0., 0.5, 1.0, 1.5};
375  std::vector <double> UO_ssth23_seeds = {pre_seed_th23_UO};
376  // if above max mxing, add ssth23mirror to seeds in UO
377  if(ssth23mirror > 0.5) UO_ssth23_seeds.push_back(ssth23mirror);
378 
379  std::map <const IFitVar*, std::vector<double>> UO_seedPts =
380  {{kFitSinSqTheta23UpperOctant, UO_ssth23_seeds},
381  {&kFitDeltaInPiUnits, UO_delta_seeds},
382  {&kFitDmSq32, {fabs(pre_seed_dmsq), -1*fabs(pre_seed_dmsq)}},
383  {&kFitSinSq2Theta13, {0.082}}};
384 
385  std::vector <const IFitVar*> UO_vars =
387  Fitter UO_fit(&expt, UO_vars, systs);
388 
389 
390  std::cout << "Starting the UO fit ----->" << endl;
391  auxShifts.ResetToNominal();
392  double UO_chi = UO_fit.Fit(calc,
393  auxShifts,
394  UO_seedPts,
395  {},
396  Fitter::kQuiet);
397  UO_dmsq = kFitDmSq32.GetValue(calc);
398  UO_ssth23 = kFitSinSqTheta23.GetValue(calc);
399  UO_delta = kFitDeltaInPiUnits.GetValue(calc);
400  UO_th13 = kFitSinSq2Theta13.GetValue(calc);
401 
402 
403  FCPoint pt(seeddelta, seedssth23, seedss2th13, true_dmsq,
404  pre_seed_th23, pre_seed_dmsq,
405  LO_delta, LO_ssth23, LO_th13, LO_dmsq,
406  LO_chi,
407  UO_delta, UO_ssth23, UO_th13, UO_dmsq,
408  UO_chi);
409 
410  fccol.AddPoint(pt);
411  unsigned int nfc_pts = fccol.NPoints();
412 
413  if(npts_flush != 0){
414  if(nfc_pts % npts_flush == 0){
415  std::cout << "Saving " + std::to_string(nfc_pts) +
416  " points to " << fcout << endl;
417  fccol.SaveToFile(fcout);
418  }
419  }
420  } // end loop over N points
421  fccol.SaveToFile(fcout);
422 }
423 
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
Definition: FitVars.cxx:141
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
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)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:171
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
void AddPoint(const FCPoint &p)
Definition: FCCollection.h:17
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:333
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)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
unsigned int NPoints() const
Definition: FCCollection.h:26
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
std::string getenv(std::string const &name)
void make_fc_oct_nersc_2019(int NPts, int N, int jid, int npts_flush=0)
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const double j
Definition: BetheBloch.cxx:29
const double kAna2019FHCLivetime
Definition: Exposures.h:226
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:147
static bool isFHC
OStream cout
Definition: OStream.cxx:6
const double kAna2019RHCPOT
Definition: Exposures.h:224
Combine multiple component experiments.
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Spectrum MockData(double pot, int idx=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:384
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2019(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:141
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
void SaveToFile(const std::string &fname) const
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
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:78
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