make_fc_oct_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_oct_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_octant_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 
137  std::string helpername =
138  helperupsname
139  + "/FCInput/"
140  + "octant_rejection_FCInput_2018.root";
141 
142  // helperupsname
143  // + "/FCInput/"
144  // + "/slices_FCInput_2018_joint_realData_bothcombo_systs_th23_noOct.root";
145 
146  TFile* fchelp = TFile::Open(helpername.c_str());
147 
148  TGraph* hssth23 = (TGraph*)fchelp->Get("ssth23");
149  double true_th23 = hssth23->Eval(0);
150 
151  // Get profiled parameters at best fit
152  TGraph *hdelta = (TGraph*)fchelp->Get("delta(pi)");
153  TGraph *hdmsq32 = (TGraph*)fchelp->Get("dmsq32");
154  TGraph *hss2th13 = (TGraph*)fchelp->Get("ss2th13");
155 
156  double seeddmsq32 = hdmsq32->Eval(0.);
157  double seeddelta = hdelta->Eval(0.);
158  double seedss2th13 = hss2th13->Eval(0.);
159  std::cout << "Got all of the parameters to make mock data" << endl;
160 
161  delete hssth23;
162  delete hdmsq32;
163  delete hdelta;
164  delete hss2th13;
165 
166  // Now, look for all the systematics
167  std::map<std::string,double> seedsyst;
168  for(const ISyst* syst :systs) {
169  TGraph* h = (TGraph*) fchelp->Get((syst->ShortName()).c_str());
170  if(!h){
171  std::cout << "Don't see a prof hist for " << syst->ShortName() << ". Continuing, but check your ups version." << std::endl;
172  continue;
173  }
174  double seedval = h->Eval(0);
175 
176  seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
177  }
178 
179  // Need random numbers to sample the bin
180  TRandom3 rnd(0);
181 
182  // Now everything is set up, we can start throwing FC mock experiments
183 
184  for(int i = 0; i < NPts; i++){
185  ResetOscCalcToDefault(calc);
186 
187 
188  // The variable you're plotting isn't a nuisance parameter, should
189  // be set to the trueX that you threw for this experiment.
190  // Other oscillation parameters are nuisance parameters and
191  // should be seeded at their best fit
192  kFitDeltaInPiUnits.SetValue(calc, seeddelta);
193  kFitSinSq2Theta13.SetValue(calc, seedss2th13);
194  kFitDmSq32.SetValue(calc, seeddmsq32);
195  kFitSinSqTheta23LowerOctant->SetValue(calc, true_th23);
196 
197  // Set syst shift - they're nuisance parameters, so set them
198  // to the best fit. Nue and Numu predictions are taught about
199  // slightly different sets of systematics, so we need to throw
200  // the predictions with different seeds
201  SystShifts seednue_fhc;
202  for(const ISyst *syst : snue_fhc){
203  if(seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
204  seednue_fhc.SetShift(syst, seedsyst[syst->ShortName()]);
205  }
206  SystShifts seednue_rhc;
207  for(const ISyst *syst : snue_rhc){
208  if(seedsyst.find(syst->ShortName()) ==seedsyst.end()) continue;
209  seednue_rhc.SetShift(syst, seedsyst[syst->ShortName()]);
210  }
211  SystShifts seednumu;
212  for(const ISyst *syst : snumu){
213  if(seedsyst.find(syst->ShortName()) ==seedsyst.end()) continue;
214  seednumu.SetShift(syst, seedsyst[syst->ShortName()]);
215  }
216 
217  // Make a mock data spectrum for nue
218  // use mock data + cosmic as predictions
219  double mock_integral_fhc = 0;
220  double mock_integral_rhc = 0;
221  assert(nue_preds.size()==2 && "Something is wrong with nue preds. Size !=2");
222  std::vector<Spectrum> mocknue;
223  SystShifts this_seednue;
224  for (int j = 0; j < (int) nue_preds.size(); j++){
225  std::cout << "\n\n Lets make some nue fake data \n\n";
226  double thisPOT = j == 0? kAna2018FHCPOT : kAna2018RHCPOT;
227  std::cout << "Getting fake data for j=" << std::to_string(j);
228  if(j==0) this_seednue = seednue_fhc;
229  else this_seednue = seednue_rhc;
230  Spectrum fakenue = nue_preds[j]->PredictSyst(calc,this_seednue);
231  std::cout << "\n\n Thats some nice nue fake data \n\n";
232  // Add in cosmic prediction here, make sure it gets Poisson fluctuated
233  fakenue += nue_cosmics_spectra[j];
234  mocknue.push_back(fakenue.MockData(thisPOT));
235  // Save this, for posterity's sake
236  if(j==0) mock_integral_fhc += mocknue[j].Integral(thisPOT);
237  if(j==1) mock_integral_rhc += mocknue[j].Integral(thisPOT);
238  }
239 
240  std::cout << "\n\n Got Nue Preds \n\n";
241  std::vector<Spectrum> mocknumu;
242  // Make a mock data spectrum for numu
243  assert(numu_cosmics_spectra.size()==8 && "Issue with N quantiles in cosmics");
244  assert(numu_preds.size()==8 && "Issue with N quantiles in predictions");
245  for (int j = 0; j < 8; j++){
246  double thisPOT = kAna2018FHCPOT;
247  if(j > 3) thisPOT = kAna2018RHCPOT;
248  Spectrum fakenm = numu_preds[j]->PredictSyst(calc,seednumu);
249  fakenm += numu_cosmics_spectra[j]; // Again, add cosmic spectra pre-fluctuation
250  mocknumu.push_back(fakenm.MockData(thisPOT)); // Fluctuate
251  }
252  std::cout << "Got Numu Preds" << endl;
253 
254 
255  // Set up experiment, nova nue + nova numu (4 quantiles) + PDG reactor
256  std::vector<const IExperiment*> expts;
257  std::vector<const IExperiment*> numu_expts;
258  // Nue experiments
259  for(int j = 0; j < 2; j++) {
260  expts.push_back(new SingleSampleExperiment(nue_preds[j], mocknue[j],
261  nue_cosmics[j].first,
262  nue_cosmics[j].second)); // Nue
263  }
264  for (int j = 0; j < 8; j++) {
265  expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
266  numu_cosmics[j].first,
267  numu_cosmics[j].second));
268  numu_expts.push_back(new SingleSampleExperiment(numu_preds[j], mocknumu[j],
269  numu_cosmics[j].first,
270  numu_cosmics[j].second));
271 
272  }
273 
274  // Add reactor
275  expts.push_back( WorldReactorConstraint2017() );
276  numu_expts.push_back(WorldReactorConstraint2017());
277  MultiExperiment expt(expts);
278  MultiExperiment numu_expt(numu_expts);
279 
280  // Correlate your systs
281  // A bit tricky since nue and numu predictions know about a
282  // different list of systs, but there are functions to do this
283  // in joint_fit_2018_tools.h
284  expt.SetSystCorrelations(0, GetCorrelations(true,true));
285  expt.SetSystCorrelations(1, GetCorrelations(true, false));
286 
287  auto notfornumu = GetCorrelations(false, false);
288 
289  for(int j = 0; j < 8; j++) {
290  expt.SetSystCorrelations(j+2, notfornumu);
291  }
292 
293  // Finally, we get to do some FC
294  //////////////////////////////////////////////////////////////////
295  // Use numu only fit to find th23 and dmsq seeds for the joint fit
296  //////////////////////////////////////////////////////////////////
297 
299  SystShifts auxShifts = SystShifts::Nominal();
300  double maxmixing = 0.514;
301  MinuitFitter fitnumu_only(&numu_expt, {kFitSinSqTheta23LowerOctant, &kFitDmSq32});
302 
303  std::cout << "\nFinding the seeds from numu only fit" << endl;
304  double minchi_numu = fitnumu_only.Fit(numu_calc, auxShifts,
305  {{&kFitDmSq32, {-2.35e-3}},
306  {kFitSinSqTheta23LowerOctant, {0.4}}},
307  {},
308  IFitter::kQuiet)->EvalMetricVal();
309  double pre_seed_th23 = kFitSinSqTheta23LowerOctant->GetValue(numu_calc);
310  double pre_seed_dmsq = kFitDmSq32.GetValue(numu_calc);
311 
312  double pre_seed_th23_LO = ( pre_seed_th23 > maxmixing ) ?
313  (2*maxmixing - pre_seed_th23) : pre_seed_th23;
314  double pre_seed_th23_UO = ( pre_seed_th23 > maxmixing ) ?
315  pre_seed_th23 : (2*maxmixing-pre_seed_th23);
316  if(pre_seed_th23_LO > 0.5) pre_seed_th23_LO = 0.45;
317 
318  // done getting pre seeds
319  //-----------------------------------------------------------------
320 
321  // Set up fitter for chi best
322  // make a vector of oscillation parameters that are our
323  // nuisance parameters
324  std::vector<const IFitVar*> LO_fitVars = {&kFitDeltaInPiUnits,
327  &kFitDmSq32};
328 
329  MinuitFitter LO_fit(&expt, LO_fitVars, systs);
330  SystShifts seedShifts = {};
331 
332  double bestdelta = 0;
333  double bestssth23 = 0;
334  double bestth13 = 0;
335  double bestdmsq = 0;
336  double bestchi = 1e10;
337  double LO_dmsq = 0;
338  double LO_delta = 0;
339  double LO_ssth23 = 0;
340  double LO_th13 = 0;
341 
342  std::vector <double> LO_delta_seeds = {0., 0.5, 1.0, 1.5};
343  std::vector <double> LO_ssth23_seeds = {pre_seed_th23_LO, pre_seed_th23_UO, maxmixing};
344  std::map <const IFitVar*, std::vector<double>> LO_seedPts =
345  {{kFitSinSqTheta23LowerOctant, {pre_seed_th23_LO}},
346  {&kFitDeltaInPiUnits, LO_delta_seeds},
347  {&kFitDmSq32, {-1*fabs(pre_seed_dmsq),fabs(pre_seed_dmsq)}}};
348  std::cout << "Starting LO fit ----->" << endl;
349  auxShifts.ResetToNominal();
350  double chi_LO = LO_fit.Fit(calc,
351  auxShifts,
352  LO_seedPts,
353  {},
354  IFitter::kQuiet)->EvalMetricVal();
355  LO_dmsq = kFitDmSq32.GetValue(calc);
356  LO_delta = kFitDeltaInPiUnits.GetValue(calc);
357  LO_ssth23 = kFitSinSqTheta23LowerOctant->GetValue(calc);
358  LO_th13 = kFitSinSq2Theta13.GetValue(calc);
359 
360  std::vector <double> bf_delta_seeds = {0., 0.5, 1.0, 1.5};
361  std::vector <double> bf_ssth23_seeds = {pre_seed_th23_LO, pre_seed_th23_UO, maxmixing};
362  std::map <const IFitVar*, std::vector<double>> bf_seedPts =
363  {{&kFitSinSqTheta23, bf_ssth23_seeds},
364  {&kFitDeltaInPiUnits, bf_delta_seeds},
365  {&kFitDmSq32, {fabs(pre_seed_dmsq), -1*fabs(pre_seed_dmsq)}},
366  {&kFitSinSq2Theta13, {0.082}}};
367 
368  std::vector <const IFitVar*> bestVars =
370  MinuitFitter bestfit(&expt, bestVars, systs);
371 
372  std::cout << "Starting the best fit ----->" << endl;
373  auxShifts.ResetToNominal();
374  double chi_UO = bestfit.Fit(calc,
375  auxShifts,
376  bf_seedPts,
377  {},
378  IFitter::kQuiet)->EvalMetricVal();
379  bestdmsq = kFitDmSq32.GetValue(calc);
380  bestssth23 = kFitSinSqTheta23.GetValue(calc);
381  bestdelta = kFitDeltaInPiUnits.GetValue(calc);
382  bestth13 = kFitSinSq2Theta13.GetValue(calc);
383  bestchi = std::min(chi_LO, chi_UO);
384 
385  FCPoint pt(seeddelta, true_th23, seedss2th13, seeddmsq32,
386  pre_seed_th23, pre_seed_dmsq,
387  LO_delta, LO_ssth23, LO_th13, LO_dmsq,
388  chi_LO,
389  bestdelta, bestssth23, bestth13, bestdmsq,
390  bestchi);
391 
392  fccol.AddPoint(pt);
393  unsigned int nfc_pts = fccol.NPoints();
394 
395  if(npts_flush != 0){
396  if(nfc_pts % npts_flush == 0){
397  std::cout << "Saving " + std::to_string(nfc_pts) +
398  " points to " << fcout << endl;
399  fccol.SaveToFile(fcout);
400  }
401  }
402  } // end loop over N points
403  fccol.SaveToFile(fcout);
404 }
405 
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
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
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
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 make_fc_oct_nersc_2018(int NPts, int N, int jid, int npts_flush=0)
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 SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:147
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)
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
void ResetToNominal()
Definition: SystShifts.cxx:144
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 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:81
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