joint_fit_2017_make_fc_surf.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 
29 using namespace ana;
30 
31 // Needed for running on the OSG in a tag prior to r20638
33 {
34  double dmsq32 = hier ? 2.44e-3 : -2.44e-3;
35  calc->SetDmsq32(dmsq32);
36 }
37 
38 
40  double& bestChi, double& bestX, double& bestY,
42 {
43  double currentChiSq = fit.Fit(calc, seed, IFitter::kQuiet)->EvalMetricVal();
44  if (currentChiSq < bestChi){
45  bestChi = currentChiSq;
46  bestX = kFitDeltaInPiUnits.GetValue(calc);
47  bestY = kFitSinSqTheta23. GetValue(calc);
48  return true;
49  }
50  return false;
51 }
52 
53 
54 void joint_fit_2017_make_fc_surf(int NPts, int bin, bool nh, int N,
56 {
57  assert(bin >= 0 && "Invalid bin number");
58 
59  std::string hierStr = nh ? "NH" : "IH";
60 
61  FCCollection fccol;
63 
64  // Relevant exposures, to be used later
65  double pot = kAna2017POT;
66  double livetime = kAna2017Livetime;
67 
68  // Grab the systs
69  std::vector<const ISyst*> systs = GetJointFitSystematicList(true, true,
70  true, true);
71  std::vector<const ISyst*> snue = GetJointFitSystematicList(true, true,
72  false,true);
73  std::vector<const ISyst*> snumu = GetJointFitSystematicList(true, false,
74  true, true);
75 
76  // We're setting up the requisite inputs to the nue experiment
77  // Nue Prediction
78  const IPrediction *p = GetNuePrediction2017("combo",calc,true,true);
79  std::pair<TH1D*,double> hcosmicnue = GetNueCosmics2017(true);
80  // We need to store cosmic hists as spectra to add into the our fake data
81  Spectrum cosmicnue(hcosmicnue.first,pot,livetime);
82 
83  // Construct numu inputs
84  std::vector<const IPrediction*> ps_numu =
85  GetNumuPredictions2017(4,true,true);
86  std::vector< std::pair<TH1D*,double> > hcosmics_nm =
87  GetNumuCosmics2017(4,true);
88  std::vector<Spectrum> cosmics_nm;
89  for (std::pair<TH1D*,double> h : hcosmics_nm)
90  cosmics_nm.push_back(Spectrum(h.first,pot,livetime));
91 
92 
93 
94  // Find the bin center for this job, probably should be arguments
95  // dCP / pi
96  double minX = 0;
97  double maxX = 2;
98  if (plot == "ssth23dmsq32"){
99  minX = 0.3;
100  maxX = 0.7;
101  }
102  int NX = 20;
103  double widthX = (maxX-minX)/NX;
104 
105  // ss theta23
106  double minY = 0.3;
107  double maxY = 0.7;
108  if (plot == "ssth23dmsq32"){
109  minY = nh ? 0.002 : -0.003;
110  maxY = nh ? 0.003 : -0.002;
111  }
112  double NY = 20;
113  double widthY = (maxY-minY)/NY;
114 
115  assert(bin < NX*NY && "Invalid bin number");
116 
117  int binY = bin/NX;
118  int binX = bin - NX*binY;
119 
120  double centerX = minX + (binX+0.5)*widthX;
121  double centerY = minY + (binY+0.5)*widthY;
122 
123 
124  // Grab the FCHelper2017 file
125  std::string helperupsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));
126  std::string helpername =
127  helperupsname+"/contours_"+plot+"_FCInput_2017.root";
128  TFile *fchelp = TFile::Open(helpername.c_str());
129  fchelp->cd();
130 
131  // Pull out hists and get values at current bin
132  double seeddmsq32 = 0;
133  if (plot == "deltassth23"){
134  TH2* hdmsq = (TH2F*)fchelp->Get((hierStr+"_DmSq32").c_str());
135  seeddmsq32 = hdmsq->GetBinContent(binX+1,binY+1); // root bins off-by-one
136  delete hdmsq;
137  }
138  double seeddelta = 0;
139  if (plot == "ssth23dmsq32"){
140  TH2* hdelta = (TH2F*)fchelp->Get((hierStr+"_DeltaInPiUnits").c_str());
141  seeddelta = hdelta->GetBinContent(binX+1,binY+1);
142  delete hdelta;
143  }
144  TH2* hss2th13 = (TH2F*)fchelp->Get((hierStr+"_SinSq2Theta13").c_str());
145  double seedss2th13 = hss2th13->GetBinContent(binX+1,binY+1);
146  delete hss2th13;
147  std::map<std::string,double> seedsyst;
148  for (const ISyst* syst : systs){
149  TH2* h = (TH2F*)fchelp->Get((hierStr+"_"+syst->ShortName()).c_str());
150  if (!h){
151  std::cout << "Don't see a prof hist for " << syst->ShortName() << ". Continuing, but check your ups version." << std::endl;
152  continue;
153  }
154  double seedval = h->GetBinContent(binX+1,binY+1);
155  seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
156  }
157 
158  // Need random numbers to sample the bin
159  TRandom3 rnd(0);
160 
161  for (int i = 0; i < NPts; i++){
162  ResetOscCalcToDefault(calc);
163 
164  if (100*i % NPts == 0){
165  std::cout << i << " / " << NPts << " " << std::endl;
166  }
167  // Throw a true X,Y value, randomly accross current bin
168  double trueX = -1;
169  // Edge bins centered on a physical edge :(
170  while (trueX < minX || trueX > maxX)
171  trueX = rnd.Uniform(centerX-widthX/2,centerX+widthX/2);
172  double trueY = -1;
173  while (trueY < minY || trueY > maxY)
174  trueY = rnd.Uniform(centerY-widthY/2,centerY+widthY/2);
175 
176  kFitDeltaInPiUnits.SetValue(calc, plot=="deltassth23"?trueX:seeddelta);
177  kFitSinSqTheta23. SetValue(calc, plot=="deltassth23"?trueY:trueX);
178  kFitDmSq32. SetValue(calc, plot=="deltassth23"?seeddmsq32:trueY);
179  kFitSinSq2Theta13.SetValue(calc,seedss2th13);
180 
181  SystShifts seednue;
182  for (const ISyst *syst : snue){
183  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
184  seednue.SetShift(syst,seedsyst[syst->ShortName()]);
185  }
186  SystShifts seednumu;
187  for (const ISyst *syst : snumu){
188  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
189  seednumu.SetShift(syst,seedsyst[syst->ShortName()]);
190  }
191 
192 
193  // Make a mock data spectrum for nue
194  // Use mock data + cosmic as predictions
195  Spectrum fakenue = p->PredictSyst(calc,seednue);
196  // Add in cosmic prediction here, make sure it gets Poisson fluctuated
197  fakenue += cosmicnue;
198  Spectrum mocknue = fakenue.MockData(pot);
199  // Save this, for posterity's sake
200  int mock_integral = mocknue.Integral(pot);
201 
202  // Make a mock data spectrum for numu
203  std::vector<Spectrum> mocknumu;
204  assert(hcosmics_nm.size()==4 && "Issue with N quantiles in cosmics");
205  assert(ps_numu.size()==4 && "Issue with N quantiles in predictions");
206  for (int i = 0; i < 4; i++){
207  Spectrum fakenm = ps_numu[i]->PredictSyst(calc,seednumu);
208  fakenm += cosmics_nm[i];
209  mocknumu.push_back(fakenm.MockData(pot));
210  }
211 
212  // Set up experiment, nova nue + nova numu + 2017 reactor
213  // Get true chi-square before fit
214  std::vector<const IExperiment*> expts;
215  expts.push_back(new SingleSampleExperiment(p, mocknue, hcosmicnue.first,
216  hcosmicnue.second)); // Nue
217  for (int i = 0; i < 4; i++)
218  expts.push_back(new SingleSampleExperiment(ps_numu[i], mocknumu[i],
219  hcosmics_nm[i].first,
220  hcosmics_nm[i].second));
221  expts.push_back( WorldReactorConstraint2017() );
222  MultiExperiment expt(expts);
223 
224  // Correlate your systs
225  expt.SetSystCorrelations(0, GetCorrelations(true));
226  auto notfornumu = GetCorrelations(false);
227  for(int i = 0; i < 4; ++i) expt.SetSystCorrelations(i+1, notfornumu);
228 
229 
230  std::vector<const IFitVar*> trueVars;
231  if (plot=="deltassth23"){
232  trueVars.push_back(&kFitDmSq32);
233  trueVars.push_back(&kFitSinSq2Theta13);
234  }
235  if (plot=="ssth23dmsq32"){
236  trueVars.push_back(&kFitDeltaInPiUnits);
237  trueVars.push_back(&kFitSinSq2Theta13);
238  }
239 
240  // Finally, we're get to do some FC
241  // Set up fitter for #chi^{2}_{true}
242  MinuitFitter fittrue(&expt, trueVars, systs);
243  double trueChi = 1e6;
244  SystShifts mecq0seed;
245  for (int i: {0,1}){
246  for (double mecq0: {-1,+1}){
247  mecq0seed.SetShift(&kMECq0ShapeSyst2017,mecq0);
248  if (plot != "deltassth23") kFitDeltaInPiUnits.SetValue(calc,i);
249  trueChi = std::min(trueChi,fittrue.Fit(calc,mecq0seed,IFitter::kQuiet)->EvalMetricVal());
250  }
251  }
252 
253  // Set up for best fit
254  double dcp = calc->GetdCP();
255  double bestX = kFitDeltaInPiUnits.GetValue(calc);
256  double bestY = kFitSinSqTheta23. GetValue(calc);
257 
258  // Set up fitter
260 
261  double bestChi = 1e6;
262  // Save best fit xy for seeding with #chi^{2}_{best}
263  int bestHie = 0; int bestDel = 0; int bestOct = 0;
264  // Seed fit at dif delta / oct values to deal with local minima
265  for (bool hier : {false, true}){
266  for (bool oct : {false, true}){
267  for (int i = 0; i <=1; i++){
268  for (double mecq0 : {-1,+1}){
269  mecq0seed.SetShift(&kMECq0ShapeSyst2017,mecq0);
270  SetHierarchy(calc, hier);
271  // Test at best fit from fitter, and the value reflected about 0.5
272  kFitSinSqTheta23.SetValue(calc,oct?bestY:2*0.5-bestY);
273  calc->SetdCP(dcp + i*M_PI);
274  bool isBest = RunFitter(fit, calc, bestChi, bestX, bestY,
275  mecq0seed);
276  // Set fit location if this is global best chisquare
277  if (isBest){
278  bestHie=hier; bestDel=i; bestOct=oct;
279  }
280  }
281  }
282  }
283  }
284  // If trueChi is lower than bestChi, we know bestChi is lower than we found
285  if (trueChi < bestChi){
286  std::cout << "Warning:" << std::endl;
287  std::cout << "Pseudo-experiment bestChi, " << bestChi << " is less than the trueChi, " << trueChi << std::endl;
288  bestChi = trueChi;
289  }
290 
291 
292  // These are all useful numbers, encode them and save as FCPoint seed
293  double fitloc = mock_integral*1e3+bestOct*1e2+bestHie*1e1+bestDel;
294 
295 
296  // Make pseudoexperiment dst
297  FCPoint pt(trueX, trueY, trueChi, bestX, bestY, bestChi, fitloc);
298  // Every once in a while, print out values we're getting
299  if (100*i % NPts == 0){
300  std::cout << i << " / " << NPts << std::endl;
301  std::cout << "Thrown parameters: " << trueX << " " << trueY << std::endl;
302  std::cout << "Best fit parameters: " << bestX << " " << bestY << std::endl;
303  std::cout << "ChiSqaures: chitrue = " << trueChi << ", chibest = " << bestChi << ", deltachisquare = " << trueChi-bestChi << std::endl;
304  }
305  fccol.AddPoint(pt);
306  }
307 
308  // Save to file
309  std::string binname = std::to_string(bin);
310  std::string hier = nh ? "_nh_" : "_ih_";
311  std::string number = std::to_string(N);
312  fccol.SaveToFile("FCCol_"+plot+"_"+binname+hier+number+".root");
313 }
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2017(const int nq=4, bool GetFromUPS=false)
void joint_fit_2017_make_fc_surf(int NPts, int bin, bool nh, int N, std::string plot)
double minY
Definition: plot_hist.C:9
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
bool RunFitter(MinuitFitter fit, osc::IOscCalcAdjustable *calc, double &bestChi, double &bestX, double &bestY, SystShifts seed=kNoShift)
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
double maxY
Definition: plot_hist.C:10
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
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
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
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:69
#define pot
virtual Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &syst) const
Definition: IPrediction.cxx:49
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
float bin[41]
Definition: plottest35.C:14
std::vector< float > Spectrum
Definition: Constants.h:610
virtual T GetdCP() const
Definition: IOscCalc.h:95
const SystShifts kNoShift
Definition: SystShifts.cxx:22
OStream cout
Definition: OStream.cxx:6
Combine multiple component experiments.
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
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
bool SetHierarchy(osc::IOscCalcAdjustable *calc, bool hier)
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
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
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
enum BeamMode string