joint_fit_2017_make_fc_slice.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 
34  double& bestChi, double& bestX, double& bestY,
36 {
37  double currentChiSq = fit.Fit(calc, seed, IFitter::kQuiet)->EvalMetricVal();
38  if (currentChiSq < bestChi){
39  bestChi = currentChiSq;
40  bestX = kFitDeltaInPiUnits.GetValue(calc);
41  bestY = kFitSinSqTheta23. GetValue(calc);
42  return true;
43  }
44  return false;
45 }
46 
47 // Needed for running on the OSG in a tag prior to r20638
49 {
50  double dmsq32 = hier ? 2.44e-3 : -2.44e-3;
51  calc->SetDmsq32(dmsq32);
52 }
53 
54 
55 // Conventionally, we run with NPts = 100 for std production of FC points
56 // with just a few for testing.
57 // We've usually set 80 bins for slice plots, so bin runs from 0-79
58 // N is just an identifier marker to add to the job output, increment if you
59 // ever need to top-up stats anywhere so you don't overwrite same-named files
60 // plot is either delta, dmsq32, or ssth23 - the three variables we make
61 // a slice plto for
62 void joint_fit_2017_make_fc_slice(int NPts, int bin, bool nh, int N,
63  std::string plot="delta")
64 {
65  assert(bin >= 0 && "Invalid bin number");
66 
67  FCCollection fccol;
69 
70 
71  // We're setting up the requisite inputs to the nue experiment
72  // Nue Prediction
73  // Relevant exposures, to be used later
74  double pot = kAna2017POT;
75  double livetime = kAna2017Livetime;
76 
77  // Grab the systs
78  std::vector<const ISyst*> systs = GetJointFitSystematicList(true, true,
79  true, true);
80  std::vector<const ISyst*> snue = GetJointFitSystematicList(true, true,
81  false,true);
82  std::vector<const ISyst*> snumu = GetJointFitSystematicList(true, false,
83  true, true);
84 
85  // We're setting up the requisite inputs to the nue experiment
86  // Nue Prediction
87  const IPrediction *p = GetNuePrediction2017("combo",calc,true,true);
88  std::pair<TH1D*,double> hcosmicnue = GetNueCosmics2017(true);
89  // We need to store cosmic hists as spectra to add into the our fake data
90  Spectrum cosmicnue(hcosmicnue.first,pot,livetime);
91 
92  // Construct numu inputs
93  std::vector<const IPrediction*> ps_numu =
94  GetNumuPredictions2017(4,true,true);
95  std::vector< std::pair<TH1D*,double> > hcosmics_nm =
96  GetNumuCosmics2017(4,true);
97  std::vector<Spectrum> cosmics_nm;
98  for (std::pair<TH1D*,double> h : hcosmics_nm)
99  cosmics_nm.push_back(Spectrum(h.first,pot,livetime));
100 
101 
102 
103 
104  // Find the extent of the bin for this job, probably should be arguments
105  // This job will randomly throw points in this bin
106  double minX = 0; // assuming plot = "delta"
107  if (plot == "ssth23") minX = 0.3;
108  if (plot == "dmsq32") minX = nh ? 2e-3 : -3e-3;
109  double maxX = 2;
110  if (plot == "ssth23") maxX = 0.7;
111  if (plot == "dmsq32") maxX = nh ? 3e-3 : -2e-3;
112  int NX = 80;
113  double widthX = (maxX-minX)/(NX-1);
114 
115  assert(bin < NX && "Invalid bin number");
116 
117  double centerX = minX + bin*widthX;
118 
119  // A note on throwing experiments - we've always set the nuisance parameters
120  // for plots, the parameters in the fit that don't show up in the plot, to
121  // the best fit at each given value of the intersting physics parameter
122  // you're plotting (the frequentist way of doing this). To do that, we
123  // need to store the best fit values of all nuisance parameters as a function
124  // of the physics parameter - that's done with TGraph's are a function of the
125  // physics parameter and there's one stored for each hierarchy / nuisance
126  // parameter combo
127 
128  // Grab the FCHelper2017 file
129  // This file saved the best fit parameters at every value of the variable
130  // we're plotting
131  std::string helperupsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));
132  std::string helpername =
133  helperupsname+"/slices_"+plot+"_FCInput_2017.root";
134  TFile *fchelp = TFile::Open(helpername.c_str());
135  fchelp->cd();
136 
137 
138  // Need to look up the right hierarchy when searching for best fits
139  std::string hiername = nh ? "NH" : "IH";
140 
141  // Pull out hists and get values at current bin
142  // Do this for the oscillation parameters first
143  // Be sure to not find a seed for the variable you're FC'ing!!!
144  double seeddelta = 0;
145  if (plot != "delta"){
146  TGraph *hdelta=(TGraph*)fchelp->Get((hiername+"_DeltaInPiUnits").c_str());
147  seeddelta = hdelta->Eval(centerX);
148  delete hdelta;
149  }
150  double seedss2th13 = 0.085;
151  TGraph *hss2th13=(TGraph*)fchelp->Get((hiername+"_SinSq2Theta13").c_str());
152  seedss2th13 = hss2th13->Eval(centerX);
153  delete hss2th13;
154  double seedssth23 = 0.6;
155  if (plot != "ssth23"){
156  TGraph *hssth23=(TGraph*)fchelp->Get((hiername+"_SinSqTheta23").c_str());
157  seedssth23 = hssth23->Eval(centerX);
158  delete hssth23;
159  }
160  double seeddmsq32 = hiername == "NH" ? 0.244 : -0.244;;
161  if (plot != "dmsq32"){
162  TGraph *hdmsq32=(TGraph*)fchelp->Get((hiername+"_DmSq32").c_str());
163  seeddmsq32 = hdmsq32->Eval(centerX);
164  delete hdmsq32;
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((hiername+"_"+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(centerX);
175  std::cout << "seeding " << syst->ShortName() << " at " << seedval << std::endl;
176  seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
177  }
178 
179 
180  // Need random numbers to sample the bin
181  TRandom3 rnd(0);
182 
183  // Now everything is set up, we can start throwing FC mock experiments
184  for (int i = 0; i < NPts; i++){
185  ResetOscCalcToDefault(calc);
186 
187  if (100*i % NPts == 0){
188  std::cout << i << " / " << NPts << " " << std::endl;
189  }
190 
191 
192 
193  // Throw a true X,Y value, randomly accross current bin
194  double trueX = -1;
195  // Edge bins centered on a physical edge :(
196  // Keep throwing until we find a value in the physical region
197  while (trueX < minX || trueX > maxX)
198  trueX = rnd.Uniform(centerX-widthX/2,centerX+widthX/2);
199 
200  // The variable you're plotting isn't a nuisance parameter, should
201  // be set to the trueX that you threw for this exeriment.
202  // Other oscillation parameters are nuisance parameters and should be
203  // seeded at their best fit
204  kFitDeltaInPiUnits.SetValue(calc, plot=="delta" ?trueX:seeddelta);
205  kFitSinSqTheta23. SetValue(calc, plot=="ssth23"?trueX:seedssth23);
206  kFitSinSq2Theta13. SetValue(calc, seedss2th13);
207  kFitDmSq32. SetValue(calc, plot=="dmsq32"?trueX:seeddmsq32);
208 
209 
210 
211  // Set syst shift - they're nuisance parameters, so set to the best fit
212  // Nue and numu predictions are taught about slightly different sets of
213  // systematics, so we need to throw the predictions with different seeds
214  SystShifts seednue;
215  for (const ISyst *syst : snue){
216  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
217  seednue.SetShift(syst,seedsyst[syst->ShortName()]);
218  }
219  SystShifts seednumu;
220  for (const ISyst *syst : snumu){
221  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
222  seednumu.SetShift(syst,seedsyst[syst->ShortName()]);
223  }
224 
225 
226 
227  // Make a mock data spectrum for nue
228  // Use mock data + cosmic as predictions
229  Spectrum fakenue = p->PredictSyst(calc,seednue);
230  // Add in cosmic prediction here, make sure it gets Poisson fluctuated
231  fakenue += cosmicnue;
232  Spectrum mocknue = fakenue.MockData(pot);
233  // Save this, for posterity's sake
234  int mock_integral = mocknue.Integral(pot);
235 
236  // Make a mock data spectrum for numu
237  std::vector<Spectrum> mocknumu;
238  assert(hcosmics_nm.size()==4 && "Issue with N quantiles in cosmics");
239  assert(ps_numu.size()==4 && "Issue with N quantiles in predictions");
240  for (int i = 0; i < 4; i++){
241  Spectrum fakenm = ps_numu[i]->PredictSyst(calc,seednumu);
242  fakenm += cosmics_nm[i]; // Again, add cosmic spectra pre-fluctuation
243  mocknumu.push_back(fakenm.MockData(pot)); // Fluctuate
244  }
245 
246  // Set up experiment, nova nue + nova numu (4 quantiles) + PDG reactor
247  std::vector<const IExperiment*> expts;
248  expts.push_back(new SingleSampleExperiment(p, mocknue, hcosmicnue.first,
249  hcosmicnue.second)); // Nue
250  for (int i = 0; i < 4; i++)
251  expts.push_back(new SingleSampleExperiment(ps_numu[i], mocknumu[i],
252  hcosmics_nm[i].first,
253  hcosmics_nm[i].second));
254  expts.push_back( WorldReactorConstraint2017() );
255  MultiExperiment expt(expts);
256 
257  // Correlate your systs
258  // A bit tricky since nue and numu predictions know about a different list
259  // of systs, but there are functions to do this in joint_fit_2017_tools.h
260  expt.SetSystCorrelations(0, GetCorrelations(true));
261  auto notfornumu = GetCorrelations(false);
262  for(int i = 0; i < 4; ++i) expt.SetSystCorrelations(i+1, notfornumu);
263 
264 
265  // Finally, we're get to do some FC
266  // Set up fitter for #chi^{2}_{best}
267 
268  // make a vector of oscillation parameters that are our nuisance parameters
269  std::vector<const IFitVar*> trueVars;
270  if (plot == "delta"){
271  trueVars.push_back(&kFitDmSq32);
272  trueVars.push_back(&kFitSinSqTheta23);
273  trueVars.push_back(&kFitSinSq2Theta13); // no delta included in this case
274  }
275  if (plot == "ssth23"){
276  trueVars.push_back(&kFitDmSq32);
277  trueVars.push_back(&kFitDeltaInPiUnits);
278  trueVars.push_back(&kFitSinSq2Theta13);
279  }
280  if (plot == "dmsq32"){
281  trueVars.push_back(&kFitDeltaInPiUnits);
282  trueVars.push_back(&kFitSinSqTheta23);
283  trueVars.push_back(&kFitSinSq2Theta13);
284  }
285  // Set up a fitter to find the best fits constrained to true value
286  MinuitFitter fittrue(&expt, trueVars, systs);
287  double trueChi = 1e6; // start the trueChi at something unphysically large
288  SystShifts mecq0seed; // related to MEC seeding problem, need to throw at
289  // +- 1 sigma, see DocDB 24987
290  // Test best fit at several seeds
291  for (double mecq0 : {-1,+1}){
292  for (double ssth23 : {0.4,0.6}){
293  for (double delta : {0,1}){
294  // Only apply multiple ssth23 / delta if you're not FC'ing that var
295  if (plot != "ssth23") kFitSinSqTheta23.SetValue(calc,ssth23);
296  if (plot != "delta") kFitDeltaInPiUnits.SetValue(calc,delta);
297  mecq0seed.SetShift(kMECq0ShapeSyst2017, mecq0);
298  double trueChiCur = fittrue.Fit(calc, mecq0seed, IFitter::kQuiet)->EvalMetricVal();
299  trueChi = std::min(trueChi,trueChiCur);
300  }
301  }
302  }
303 
304  // Save our seeds for best fit
305  double dcp = calc->GetdCP();
306  double bestX = kFitDeltaInPiUnits.GetValue(calc);
307  double bestY = kFitSinSqTheta23 .GetValue(calc);
308 
309  // Set up fitter for calculating the global best fit
312 
313  double bestChi = 1e6;
314  int bestHie = 0; int bestDel = 0; int bestOct = 0;
315  // Seed fit at dif delta / oct values to deal with local minima
316  for (bool hier : {false, true}){
317  for (bool oct : {false, true}){
318  for (int i = 0; i <=1; i++){
319  for (double mecq0 : {-1,+1}){
320  mecq0seed.SetShift(kMECq0ShapeSyst2017,mecq0);
321  SetHierarchy(calc, hier);
322  // Test at best fit from fitture, and the value reflected about 0.5
323  kFitSinSqTheta23.SetValue(calc,oct?bestY:2*0.5-bestY);
324  calc->SetdCP(dcp + i*M_PI);
325  bool isBest = RunFitter(fit, calc, bestChi, bestX, bestY,
326  mecq0seed);
327  // Set fit location if this is global best chisquare
328  if (isBest){
329  bestHie=hier; bestDel=i; bestOct=oct;
330  }
331  }
332  }
333  }
334  }
335  // If trueChi is lower than bestChi, we know bestChi is lower than we found
336  if (trueChi < bestChi){
337  std::cout << "Warning:" << std::endl;
338  std::cout << "Pseudo-experiment bestChi, " << bestChi << " is less than the trueChi, " << trueChi << std::endl;
339  bestChi = trueChi;
340  }
341 
342 
343  // These are all useful numbers, encode them and save as FCPoint seed
344  double fitloc = mock_integral*1e3+bestOct*1e2+bestHie*1e1+bestDel;
345 
346 
347  // Make pseudoexperiment dst
348  FCPoint pt(trueX, seedssth23, trueChi, bestX, bestY, bestChi, fitloc);
349  // Every once in a while, print out values we're getting
350  if (100*i % NPts == 0){
351  std::cout << i << " / " << NPts << std::endl;
352  std::cout << "Thrown parameters: " << trueX << " " << seedssth23 << std::endl;
353  std::cout << "Best fit parameters: " << bestX << " " << bestY << std::endl;
354  std::cout << "ChiSqaures: chitrue = " << trueChi << ", chibest = " << bestChi << ", deltachisquare = " << trueChi-bestChi << std::endl;
355  }
356  fccol.AddPoint(pt);
357  }
358 
359  // Save to file
360  // Do our very best to insure we never have a same-named file and over-write
361  // good FC points.
362  std::string binname = std::to_string(bin)+"_";
363  std::string hier = nh ? "_nh_" : "_ih_";
364  std::string number = std::to_string(N);
365  fccol.SaveToFile("FCCol_"+plot+hier+binname+number+".root");
366 }
void joint_fit_2017_make_fc_slice(int NPts, int bin, bool nh, int N, std::string plot="delta")
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2017(const int nq=4, bool GetFromUPS=false)
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
bool SetHierarchy(osc::IOscCalcAdjustable *calc, bool hier)
double ssth23
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double delta
Definition: runWimpSim.h:98
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
const char * p
Definition: xmltok.h:285
void AddPoint(const FCPoint &p)
Definition: FCCollection.h:17
const IPrediction * GetNuePrediction2017(std::string decomp, osc::IOscCalc *calc, bool corrSysts, bool GetFromUPS=false)
const double kAna2017POT
Definition: Exposures.h:174
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
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
virtual void SetDmsq32(const T &dmsq32)=0
#define M_PI
Definition: SbMath.h:34
bool RunFitter(MinuitFitter fit, osc::IOscCalcAdjustable *calc, double &bestChi, double &bestX, double &bestY, SystShifts seed=kNoShift)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
Represents the results of a single Feldman-Cousins pseudo-experiment.
Definition: FCPoint.h:8
expt
Definition: demo5.py:34
std::pair< TH1D *, double > GetNueCosmics2017(bool GetFromUPS=false)
const double kAna2017Livetime
Definition: Exposures.h:200
std::vector< const IPrediction * > GetNumuPredictions2017(const int nq=4, bool useSysts=true, bool GetFromUPS=false)
unsigned int seed
Definition: runWimpSim.h:102
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:64
#define pot
virtual Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:72
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
float bin[41]
Definition: plottest35.C:14
std::vector< float > Spectrum
Definition: Constants.h:527
virtual T GetdCP() const
Definition: IOscCalc.h:90
const SystShifts kNoShift
Definition: SystShifts.h:115
OStream cout
Definition: OStream.cxx:6
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
double dmsq32
double livetime
Definition: saveFDMCHists.C:21
const NOvARwgtSyst kMECq0ShapeSyst2017("MECq0Shape","MEC q_{0} shape", novarwgt::kMECq0ShapeSyst2017)
Definition: MECSysts.h:41
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
General interface to any calculator that lets you set the parameters.
assert(nhit_max >=nhit_nbins)
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
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
Float_t e
Definition: plot.C:35
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)
virtual void SetdCP(const T &dCP)=0
Compare a single data spectrum to the MC + cosmics expectation.
void plot(std::string label, std::map< std::string, std::map< std::string, Spectrum * >> &plots, bool log)
Collection of FCPoint. Serializable to/from a file.
Definition: FCCollection.h:12
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17