make_fc_mh_nersc_2018.C
Go to the documentation of this file.
2 #include "CAFAna/Core/Spectrum.h"
4 
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 
22 
23 #include "CAFAna/Analysis/Calcs.h"
24 #include "OscLib/IOscCalc.h"
25 
26 #include "TRandom3.h"
27 #include "TH1.h"
28 #include "TGraph.h"
29 
30 using namespace ana;
31 
32 
33 void make_fc_mh_nersc_2018(int NPts, int N, int jid, int npts_flush = 0)
34 {
35  std::string number = std::to_string(N);
36  std::string fcout = "FCCol_mass_hierarchy_number_v"
37  + number
38  + "_"
39  + std::to_string(jid)
40  +".root";
41 
42  FCCollection fccol;
44 
45  double thisPOT;
46  double thisLivetime;
47 
48  // Grab the systs
49  // bools to make the function calls a little more readable
50  bool nueOnly = true;
51  bool numuOnly= true;
52  bool isFHC = true;
53  bool isRHC = true;
54 
55  // joint
56  std::vector<const ISyst*> systs = GetJointFitSystematicList(true, !nueOnly, !numuOnly, true, true, true);
57 
58  // nue fhc
59  std::vector<const ISyst*> snue_fhc = getAllAna2018Systs(kNueAna2018, true, kFHC);
60 
61  // nue rhc
62  std::vector<const ISyst*> snue_rhc = getAllAna2018Systs(kNueAna2018, true, kRHC);
63 
64  // numu
65  std::vector<const ISyst*> snumu = GetJointFitSystematicList(true, !nueOnly, numuOnly, isFHC, isRHC, true);
66 
67  // We're setting up the requisite inputs to the nue experiment
68  // nue prediction
69 
70  std::vector <const IPrediction*> nue_preds;
71  std::vector <std::pair <TH1D*, double > > nue_cosmics;
72  std::vector <Spectrum> nue_cosmics_spectra;
73  bool minimizeMemory = true;
74 
75  // Grab the nue predictions for fhc and rhc
76  nue_preds.push_back(GetNuePrediction2018("combo", calc, true, "fhc", false, false, minimizeMemory, true));
77  nue_preds.push_back(GetNuePrediction2018("prop", calc, true, "rhc", false, false, minimizeMemory, true));
78 
79  std::cout << "\nDone getting Nue Predictions\n\n" ;
80 
81  nue_cosmics.push_back(GetNueCosmics2018("fhc", false, true));
82  nue_cosmics.push_back(GetNueCosmics2018("rhc", false, true));
83 
84  // We need to store cosmic hists as spectra to add into our fake data
85  // FHC and RHC have different POT and Livetimes
86  for(int i = 0; i < 2; i++) {
87  if(i==0) {
88  thisPOT = kAna2018FHCPOT;
89  thisLivetime = kAna2018FHCLivetime;
90  }
91  else {
92  thisPOT = kAna2018RHCPOT;
93  thisLivetime = kAna2018RHCLivetime;
94  }
95  nue_cosmics_spectra.push_back(Spectrum(nue_cosmics[i].first,
96  thisPOT, thisLivetime));
97  }
98 
99 
100  //////////////////////////////////////////////////////////////////////
101  // Numu Prediction
102  //////////////////////////////////////////////////////////////////////
103  std::vector<const IPrediction*> numu_preds;
104  std::vector< std::pair<TH1D*, double> > numu_cosmics;
105  std::vector<Spectrum> numu_cosmics_spectra;
106 
107  std::vector<const IPrediction*> this_numu_preds;
108  std::vector< std::pair<TH1D*, double> > this_numu_cosmics;
109 
110  // Construct numu inputs
111  for(std::string beam : {"fhc", "rhc"}){
112  this_numu_preds = GetNumuPredictions2018(4,true, beam, false, kNuMu, minimizeMemory, true);
113  this_numu_cosmics = GetNumuCosmics2018(4, beam, false, true);
114  numu_preds.insert(numu_preds.end(),
115  this_numu_preds.begin(),
116  this_numu_preds.end());
117  numu_cosmics.insert(numu_cosmics.end(),
118  this_numu_cosmics.begin(),
119  this_numu_cosmics.end());
120 
121  for(std::pair<TH1D*, double> h : this_numu_cosmics){
122  thisPOT = kAna2018FHCPOT;
123  thisLivetime = kAna2018FHCLivetime;
124  if(beam == "rhc") {
125  thisPOT = kAna2018RHCPOT;
126  thisLivetime = kAna2018RHCLivetime;
127  }
128  numu_cosmics_spectra.push_back(Spectrum(h.first, thisPOT, thisLivetime));
129  }
130 
131  }
132 
133  // Grab analysis slices to find the bin centers and the best
134  // fit point
135  std::string helperupsname = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
136  // incorrectly named 2017, but don't worry
137  std::string anaslicename = helperupsname +
138  "/slices/hist_slices_2017_joint_realData_bothcombo_systs_dmsq_noOct.root";
139 
140  TFile* anaslicefile = TFile::Open(anaslicename.c_str());
141  TH1* anaslice = (TH1F*) anaslicefile->Get("slice_dmsq_IH");
142 
143  // best fit bin
144  // get upper and lower edge to get a range of dmsq
145  // to throw mock experiments in
146  int best_bin = anaslice->GetMinimumBin();
147  double true_dmsq = anaslice->GetBinCenter(best_bin);
148  double best_bin_low = anaslice->GetBinLowEdge(best_bin);
149  double best_bin_high = anaslice->GetBinLowEdge(best_bin+1);
150  double widthX = anaslice->GetBinWidth(best_bin);
151  std::cout << "Got best fit point bin info" << endl;
152 
153  std::string helpername =
154  helperupsname
155  + "/FCInput/"
156  + "/slices_FCInput_2018_joint_realData_bothcombo_systs_dmsq_noOct.root";
157 
158  TFile* fchelp = TFile::Open(helpername.c_str());
159 
160  // Get profiled parameters at best fit
161  TGraph *hdelta = (TGraph*)fchelp->Get("IH_DeltaCPpi");
162  TGraph *hssth23 = (TGraph*)fchelp->Get("IH_SinSqTheta23");
163  TGraph *hss2th13 = (TGraph*)fchelp->Get("IH_SinSq2Theta13");
164 
165  double seeddelta = hdelta->Eval(true_dmsq);
166  double seedssth23 = hssth23->Eval(true_dmsq);
167  double seedss2th13 = hss2th13->Eval(true_dmsq);
168  std::cout << "Got all of the parameters to make mock data" << endl;
169 
170  delete hdelta;
171  delete hssth23;
172  delete hss2th13;
173 
174  // Now, look for all the systematics
175  std::map<std::string,double> seedsyst;
176  for(const ISyst* syst :systs) {
177  TGraph* h = (TGraph*) fchelp->Get(("IH_"+syst->ShortName()).c_str());
178  if(!h){
179  std::cout << "Don't see a prof hist for " << syst->ShortName() << ". Continuing, but check your ups version." << std::endl;
180  continue;
181  }
182  double seedval = h->Eval(true_dmsq);
183 
184  seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
185  }
186 
187 
188  // Need random numbers to sample the bin
189  TRandom3 rnd(0);
190 
191  // Now everything is set up, we can start throwing FC mock experiments
192 
193  for(int i = 0; i < NPts; i++){
194  ResetOscCalcToDefault(calc);
195 
196 
197  // The variable you're plotting isn't a nuisance parameter, should
198  // be set to the trueX that you threw for this experiment.
199  // Other oscillation parameters are nuisance parameters and
200  // should be seeded at their best fit
201  kFitDeltaInPiUnits.SetValue(calc, seeddelta);
202  kFitSinSqTheta23.SetValue(calc, seedssth23);
203  kFitSinSq2Theta13.SetValue(calc, seedss2th13);
204  kFitDmSq32InvertedHierarchy.SetValue(calc, true_dmsq);
205 
206  // Set syst shift - they're nuisance parameters, so set them
207  // to the best fit. Nue and Numu predictions are taught about
208  // slightly different sets of systematics, so we need to throw
209  // the predictions with different seeds
210  SystShifts seednue_fhc;
211  for(const ISyst *syst : snue_fhc){
212  if(seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
213  seednue_fhc.SetShift(syst, seedsyst[syst->ShortName()]);
214  }
215  SystShifts seednue_rhc;
216  for(const ISyst *syst : snue_rhc){
217  if(seedsyst.find(syst->ShortName()) ==seedsyst.end()) continue;
218  seednue_rhc.SetShift(syst, seedsyst[syst->ShortName()]);
219  }
220  SystShifts seednumu;
221  for(const ISyst *syst : snumu){
222  if(seedsyst.find(syst->ShortName()) ==seedsyst.end()) continue;
223  seednumu.SetShift(syst, seedsyst[syst->ShortName()]);
224  }
225 
226  // Make a mock data spectrum for nue
227  // use mock data + cosmic as predictions
228  double mock_integral_fhc = 0;
229  double mock_integral_rhc = 0;
230  assert(nue_preds.size()==2 && "Something is wrong with nue preds. Size !=2");
231  std::vector<Spectrum> mocknue;
232  SystShifts this_seednue;
233  for (int j = 0; j < (int) nue_preds.size(); j++){
234  std::cout << "\n\n Lets make some nue fake data \n\n";
235  double thisPOT = j == 0? kAna2018FHCPOT : kAna2018RHCPOT;
236  std::cout << "Getting fake data for j=" << std::to_string(j);
237  if(j==0) this_seednue = seednue_fhc;
238  else this_seednue = seednue_rhc;
239  Spectrum fakenue = nue_preds[j]->PredictSyst(calc,this_seednue);
240  std::cout << "\n\n Thats some nice nue fake data \n\n";
241  // Add in cosmic prediction here, make sure it gets Poisson fluctuated
242  fakenue += nue_cosmics_spectra[j];
243  mocknue.push_back(fakenue.MockData(thisPOT));
244  // Save this, for posterity's sake
245  if(j==0) mock_integral_fhc += mocknue[j].Integral(thisPOT);
246  if(j==1) mock_integral_rhc += mocknue[j].Integral(thisPOT);
247  }
248 
249  std::cout << "\n\n Got Nue Preds \n\n";
250  std::vector<Spectrum> mocknumu;
251  // Make a mock data spectrum for numu
252  assert(numu_cosmics_spectra.size()==8 && "Issue with N quantiles in cosmics");
253  assert(numu_preds.size()==8 && "Issue with N quantiles in predictions");
254  for (int j = 0; j < 8; j++){
255  double thisPOT = kAna2018FHCPOT;
256  if(j > 3) thisPOT = kAna2018RHCPOT;
257  Spectrum fakenm = numu_preds[j]->PredictSyst(calc,seednumu);
258  fakenm += numu_cosmics_spectra[j]; // Again, add cosmic spectra pre-fluctuation
259  mocknumu.push_back(fakenm.MockData(thisPOT)); // Fluctuate
260  }
261  std::cout << "Got Numu Preds" << endl;
262 
263 
264  // Set up experiment, nova nue + nova numu (4 quantiles) + PDG reactor
265  std::vector<const IExperiment*> expts;
266  std::vector<const IExperiment*> numu_expts;
267  // Nue experiments
268  for(int j = 0; j < 2; j++) {
269  expts.push_back(new SingleSampleExperiment(nue_preds[j], mocknue[j],
270  nue_cosmics[j].first,
271  nue_cosmics[j].second)); // Nue
272  }
273  for (int j = 0; j < 8; j++) {
274  expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
275  numu_cosmics[j].first,
276  numu_cosmics[j].second));
277  numu_expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
278  numu_cosmics[j].first,
279  numu_cosmics[j].second));
280 
281  }
282 
283  // Add reactor
284  expts.push_back( WorldReactorConstraint2017() );
285  numu_expts.push_back(WorldReactorConstraint2017());
286  MultiExperiment expt(expts);
287  MultiExperiment numu_expt(numu_expts);
288 
289  // Correlate your systs
290  // A bit tricky since nue and numu predictions know about a
291  // different list of systs, but there are functions to do this
292  // in joint_fit_2018_tools.h
293  expt.SetSystCorrelations(0, GetCorrelations(true,true));
294  expt.SetSystCorrelations(1, GetCorrelations(true, false));
295 
296  auto notfornumu = GetCorrelations(false, false);
297 
298  for(int j = 0; j < 8; j++) {
299  expt.SetSystCorrelations(j+2, notfornumu);
300  }
301 
302  // Finally, we get to do some FC
303  //////////////////////////////////////////////////////////////////
304  // Use numu only fit to find th23 and dmsq seeds for the joint fit
305  //////////////////////////////////////////////////////////////////
306 
308  SystShifts auxShifts = SystShifts::Nominal();
309  double maxmixing = 0.514;
310  MinuitFitter fitnumu_only(&numu_expt, {&kFitSinSqTheta23, &kFitDmSq32InvertedHierarchy});
311 
312  std::cout << "\nFinding the seeds from numu only fit" << endl;
313  double minchi_numu = fitnumu_only.Fit(numu_calc, auxShifts,
314  {{&kFitDmSq32InvertedHierarchy, {-2.35e-3}},
315  {&kFitSinSqTheta23, {0.4}}},
316  {},
317  IFitter::kQuiet)->EvalMetricVal();
318  double pre_seed_th23 = kFitSinSqTheta23.GetValue(numu_calc);
319  double pre_seed_dmsq = kFitDmSq32.GetValue(numu_calc);
320 
321  double pre_seed_th23_LO = ( pre_seed_th23 > maxmixing ) ?
322  (2*maxmixing - pre_seed_th23) : pre_seed_th23;
323  double pre_seed_th23_UO = ( pre_seed_th23 > maxmixing ) ?
324  pre_seed_th23 : (2*maxmixing-pre_seed_th23);
325  if(pre_seed_th23_LO > 0.5) pre_seed_th23_LO = 0.45;
326 
327  // done getting pre seeds
328  //-----------------------------------------------------------------
329 
330  // Set up fitter for chi best
331  // make a vector of oscillation parameters that are our
332  // nuisance parameters
333  std::vector<const IFitVar*> ih_fitVars = {&kFitDeltaInPiUnits,
337 
338  MinuitFitter ih_fit(&expt, ih_fitVars, systs);
339  SystShifts seedShifts = {};
340 
341  double bestdelta = 0;
342  double bestssth23 = 0;
343  double bestth13 = 0;
344  double bestdmsq = 0;
345  double bestchi = 1e10;
346  double ih_dmsq = 0;
347  double ih_delta = 0;
348  double ih_ssth23 = 0;
349  double ih_th13 = 0;
350 
351  std::vector <double> ih_delta_seeds = {0., 0.5, 1.0, 1.5};
352  std::vector <double> ih_ssth23_seeds = {pre_seed_th23_LO, pre_seed_th23_UO, maxmixing};
353  std::map <const IFitVar*, std::vector<double>> ih_seedPts =
354  {{&kFitSinSqTheta23, ih_ssth23_seeds},
355  {&kFitDeltaInPiUnits, ih_delta_seeds},
356  {&kFitDmSq32InvertedHierarchy, {-1*fabs(pre_seed_dmsq)}}};
357  std::cout << "Starting IH fit ----->" << endl;
358  auxShifts.ResetToNominal();
359  double chi_IH = ih_fit.Fit(calc,
360  auxShifts,
361  ih_seedPts,
362  {},
363  IFitter::kQuiet)->EvalMetricVal();
364  ih_dmsq = kFitDmSq32InvertedHierarchy.GetValue(calc);
365  ih_delta = kFitDeltaInPiUnits.GetValue(calc);
366  ih_ssth23 = kFitSinSqTheta23.GetValue(calc);
367  ih_th13 = kFitSinSq2Theta13.GetValue(calc);
368 
369  std::vector <double> bf_delta_seeds = {0., 0.5, 1.0, 1.5};
370  std::vector <double> bf_ssth23_seeds = {pre_seed_th23_LO, pre_seed_th23_UO, maxmixing};
371  std::map <const IFitVar*, std::vector<double>> bf_seedPts =
372  {{&kFitSinSqTheta23, bf_ssth23_seeds},
373  {&kFitDeltaInPiUnits, bf_delta_seeds},
374  {&kFitDmSq32, {fabs(pre_seed_dmsq), -1*fabs(pre_seed_dmsq)}},
375  {&kFitSinSq2Theta13, {0.082}}};
376 
377  std::vector <const IFitVar*> bestVars =
379  MinuitFitter bestfit(&expt, bestVars, systs);
380 
381  std::cout << "Starting the best fit ----->" << endl;
382  auxShifts.ResetToNominal();
383  double chi_NH = bestfit.Fit(calc,
384  auxShifts,
385  bf_seedPts,
386  {},
387  IFitter::kQuiet)->EvalMetricVal();
388  bestdmsq = kFitDmSq32.GetValue(calc);
389  bestssth23 = kFitSinSqTheta23.GetValue(calc);
390  bestdelta = kFitDeltaInPiUnits.GetValue(calc);
391  bestth13 = kFitSinSq2Theta13.GetValue(calc);
392  bestchi = std::min(chi_NH, chi_IH);
393 
394  FCPoint pt(seeddelta, seedssth23, seedss2th13, true_dmsq,
395  pre_seed_th23, pre_seed_dmsq,
396  ih_delta, ih_ssth23, ih_th13, ih_dmsq,
397  chi_IH,
398  bestdelta, bestssth23, bestth13, bestdmsq,
399  bestchi);
400 
401  fccol.AddPoint(pt);
402  unsigned int nfc_pts = fccol.NPoints();
403 
404  if(npts_flush != 0){
405  if(nfc_pts % npts_flush == 0){
406  std::cout << "Saving " + std::to_string(nfc_pts) +
407  " points to " << fcout << endl;
408  fccol.SaveToFile(fcout);
409  }
410  }
411  } // end loop over N points
412  fccol.SaveToFile(fcout);
413 }
414 
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
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
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Definition: FitVars.cxx:171
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()
std::vector< const IPrediction * > GetNumuPredictions2018(const int nq=4, bool useSysts=true, std::string beam="fhc", bool GetFromUPS=false, ENu2018ExtrapType numuExtrap=kNuMu, bool minimizeMemory=false, bool NERSC=false)
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
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
const double kAna2018RHCPOT
Definition: Exposures.h:208
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
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::string getenv(std::string const &name)
virtual std::unique_ptr< IFitSummary > Fit(osc::IOscCalcAdjustable *seed, SystShifts &bestSysts=junkShifts, const SeedList &seedPts=SeedList(), const std::vector< SystShifts > &systSeedPts={}, Verbosity verb=kVerbose) const
Master fitting method. Depends on FitHelper and FitHelperSeeded.
Definition: IFitter.cxx:69
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const double j
Definition: BetheBloch.cxx:29
void make_fc_mh_nersc_2018(int NPts, int N, int jid, int npts_flush=0)
std::vector< float > Spectrum
Definition: Constants.h:610
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2018(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
static bool isFHC
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
std::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool GetFromUPS=false, bool NERSC=false)
void ResetToNominal()
Definition: SystShifts.cxx:143
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
assert(nhit_max >=nhit_nbins)
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
const IPrediction * GetNuePrediction2018(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
const double kAna2018FHCPOT
Definition: Exposures.h:207
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
T min(const caf::Proxy< T > &a, T b)
void SaveToFile(const std::string &fname) const
const FitDmSq32InvertedHierarchy kFitDmSq32InvertedHierarchy
Definition: FitVars.cxx:23
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:80
const double kAna2018FHCLivetime
Definition: Exposures.h:213
void SetSystCorrelations(int idx, const std::vector< std::pair< const ISyst *, const ISyst * >> &corrs)
const double kAna2018RHCLivetime
Definition: Exposures.h:214
Compare a single data spectrum to the MC + cosmics expectation.
Collection of FCPoint. Serializable to/from a file.
Definition: FCCollection.h:12
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string