make_fc_surfaces_2020.C
Go to the documentation of this file.
1 #include <fstream>
2 
4 #include "CAFAna/Core/Spectrum.h"
6 
7 #include "CAFAna/Fit/Fit.h"
9 
12 
13 #include "CAFAna/Vars/FitVars.h"
14 
15 #include "CAFAna/FC/FCPoint.h"
16 #include "CAFAna/FC/FCCollection.h"
17 #include "CAFAna/FC/FCSurface.h"
18 
22 
24 
25 #include "CAFAna/Analysis/Calcs.h"
26 #include "OscLib/IOscCalc.h"
27 
28 #include "TRandom3.h"
29 #include "TH1.h"
30 
31 namespace fc {
32  enum kBeam_t {
35  };
36  enum kFlavor_t {
39  };
40 }
41 
42 using namespace ana;
43 
44 void SetSeeds(osc::IOscCalcAdjustable* calc, double delta, double th23, double dmsq)
45 {
47  kFitDeltaInPiUnits.SetValue(calc, delta);
48  kFitSinSqTheta23.SetValue(calc, th23);
49  kFitDmSq32.SetValue(calc, dmsq);
50 }
51 
52 void make_fc_surfaces_2020(int NPts, int startidx, int endidx, bool nh,
53  int N, int jid, int npts_flush = 0,
54  std::string plot = "deltassth23" ,
55  std::string options = "joint_fake2019" ,
56  std::string decomp = "",
57  TString beam = "both", bool corrSysts = false,
58  bool readFromFile = false, bool runOnGrid=false)
59 {
60  assert(plot == "deltath13" || plot=="ssth23dmsq32" ||
61  plot == "deltassth23" && "Plot must be ssth23dmsq32 or deltassth23");
62 
63  assert(beam == "FHCOnly" || beam == "RHCOnly" || beam == "both");
64 
65  if(readFromFile && !corrSysts){
66  std::cout << "No bin list for stats only contours" << endl;
67  // exit(1);
68  }
69 
70  //also treat the decomp option to get the right FCInput name
71 
72  bool fake2019 = options.find("fake2019") != std::string::npos;
73  bool realData = options.find("realData") != std::string::npos;
74 
75  bool RHCOnly = beam.Contains("RHCOnly");
76  bool FHCOnly = beam.Contains("FHCOnly");
77  bool both = beam.Contains("both");
78 
79  bool nueOnly = options.find("nueOnly") != std::string::npos;
80  bool numuOnly = options.find("numuOnly") != std::string::npos;
81  bool joint = options.find("joint") != std::string::npos;
82 
83  bool isFHC = true;
84  bool isRHC = true;
85 
86  // ssth23dmsq32 has 400 total bins while deltassth23 has 900
87  int TotBins = 0;
88  TotBins = (plot == "deltassth23") ? 900 : 400;
89  int bins[TotBins];
90  std::memset(bins, 0, TotBins);
91  std::string hier = nh? "_NH_" : "_IH_";
92  std::string stat = corrSysts? "_systs": "_stats";
93 
94 
95  // If we don't want to run full corrections, we'll read from a text file which holds
96  // a list of bins to be corrected based on deltachisq in each bin. This list is
97  // generated by generate_fc_binlists.C
98 
99  ifstream binlist;
100  // std::string path = getenv("FCHELPERANA2019_LIB_PATH");
101  std::string path = "/nova/app/users/scalvez/fchelper/";
102  path += "/binlists/";
103  // TString binlist_fname = path+plot+hier+beam+stat+"_binlist.txt";
104  TString binlist_fname = path+plot+hier+"binlist.txt";
105 
106  binlist.open(binlist_fname);
107  int j = 0;
108  // check if the binlist exists
109  if(!binlist.good()){
110  std::cout << "Could not find "<< binlist_fname << endl;
111  exit(1);
112  }
113  if(readFromFile){
114  do
115  {
116  binlist >> bins[j];
117  j++;
118  }while(!binlist.eof());
119  binlist.close();
120  }
121  else for (int ibin = 0; ibin < TotBins; ibin++) bins[ibin]=ibin;
122 
123  int ibin = bins[startidx];
124  int fbin = bins[endidx];
125 
126  std::string binname = "bins_"+std::to_string(ibin) + "_to_" + std::to_string(fbin);
127  std::string number = std::to_string(N);
128  std::string fcout = "FCCol_"+plot+"_"+hier+"_";
129  if(both) fcout+="both";
130  else fcout += RHCOnly? "RHCOnly":"FHCOnly";
131  fcout += corrSysts? "_systs_":"_stats_";
132  fcout += binname+"_v"+number+'_'+std::to_string(jid)+".root";
133 
134  std::string hierStr = nh ? "NH" : "IH";
135 
136  FCCollection fccol;
137 
138  // Relevant exposures, to be used later
139  double thisPOT;
140  double thisLivetime;
141 
143  if(fake2019) SetFakeCalc(calc, 0.564956, 2.48304e-3, 1.9997);
144  // if(fake2019) SetFakeCalc(calc, 0.565, 2.48e-3, 0);
145  else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}
146 
147  // struct predictions {
148  // const string name;
149  // const IPrediction * pred;
150  // std::pair <Spectrum*, double> cos;
151  // double pot;
152  // double livetime;
153  // const fc::kBeam_t beam;
154  // const fc::kFlavor_t flavor;
155  // };
156 
157  // std::vector <predictions> preds;
158  // std::vector <predictions> preds_numu_only;
159  // std::vector <Spectrum * > data;
160  // std::vector <Spectrum * > data_numu_only;
161 
162  std::vector <ana::predictions> preds = LoadPredictions (corrSysts, runOnGrid, decomp, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
163  std::vector <ana::predictions> preds_numu_only = LoadPredictions (corrSysts, runOnGrid, decomp, false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
164  // std::vector <Spectrum * > data = LoadRealData (runOnGrid, nueOnly || joint, numuOnly || joint, FHCOnly || both, RHCOnly || both);
165  // std::vector <Spectrum * > data_numu_only = LoadRealData (runOnGrid, false, numuOnly || joint, FHCOnly || both, RHCOnly || both);
166 
167  // We're setting up the requisite inputs to the nue experiment
168  /*
169  //////////////////////////////////////////////////////////////
170  // Nue Predictions
171  /////////////////////////////////////////////////////////////
172 
173  if(!numuOnly) {
174  if(FHCOnly || both ) {
175  preds.push_back({"nue_fhc",
176  GetNuePrediction2020(decomp+"combo", DefaultOscCalc(), corrSysts, "fhc",
177  true, runOnGrid),
178  GetNueCosmics2020("fhc", runOnGrid),
179  GetPOT(), GetLT(),
180  fc::kFHC, fc::kNue});
181  }
182  if(RHCOnly || both ) {
183  preds.push_back({"nue_rhc",
184  GetNuePrediction2020("prop", DefaultOscCalc(), corrSysts, "rhc",
185  false, runOnGrid),
186  GetNueCosmics2020("rhc", runOnGrid),
187  GetPOT(false), GetLT(false),
188  fc::kRHC, fc::kNue});
189  }
190  }
191  std::cout << "\n------------------Done Getting nue preds ---------------------\n";
192 
193  ///////////////////////////////////////////////////////////////
194  // Numu predictions
195  //////////////////////////////////////////////////////////////
196 
197  // Construct numu inputs
198  int nnumu = 4;
199  if(!nueOnly) {
200  for(fc::kBeam_t _beam : {fc::kFHC, fc::kRHC}) {
201  if(((_beam == fc::kRHC) && FHCOnly) ||
202  ((_beam == fc::kFHC) && RHCOnly))
203  continue;
204  double _pot = GetPOT(_beam == fc::kFHC);
205  double _livetime = GetLT (_beam == fc::kFHC);
206  std::string beam_name = _beam == fc::kFHC ? "fhc" : "rhc";
207  auto numu_preds = GetNumuPredictions2020(nnumu, decomp, DefaultOscCalc(), corrSysts, beam_name, false, false);
208  auto numu_cosmics = GetNumuCosmics2020 (nnumu, beam_name, false);
209  for (int i = 0; i< nnumu; i++ ){
210  preds.push_back({"numu_" + beam_name + "_" +std::to_string(i+1),
211  numu_preds[i],
212  numu_cosmics[i],
213  _pot,
214  _livetime,
215  _beam,
216  fc::kNumu});
217  preds_numu_only.push_back({"numu_" + beam_name + "_"+std::to_string(i+1),
218  numu_preds[i],
219  numu_cosmics[i],
220  _pot,
221  _livetime,
222  _beam,
223  fc::kNumu});
224  }
225  }
226  }
227  */
228  // Grab the systs
229  bool ptExtrap = false;
230  // std::vector<const ISyst*> systs = {};
231  auto systs = GetJointFitSystematicList(corrSysts, nueOnly, numuOnly, FHCOnly||both, RHCOnly||both, true, ptExtrap);
232  // if(corrSysts) systs = GetJointFitSystematicList(true, !nueOnly,
233  // !numuOnly, isFHC, isRHC);
234  // Equivalent for 2020?
235  // need to grab different systs for nue rhc and fhc
236  std::vector<const ISyst*> snue_fhc = {};
237  // if(corrSysts) snue_fhc = getAllAna2018Systs(kNueAna2018, true, kFHC);
238  std::vector<const ISyst*> snue_rhc = {};
239  // if(corrSysts) snue_rhc = getAllAna2018Systs(kNueAna2018, true, kRHC);
240 
241  // numu systs are independent of beam mode
242  std::vector<const ISyst*> snumu = {};
243  // snumu = GetJointFitSystematicList(true, !nueOnly,
244  // numuOnly,
245  // isFHC,
246  // isRHC);
247 
248  // std::cout << "Systematics for the numu only fit:\n";
249  auto systs_numu_only = GetJointFitSystematicList(corrSysts, false, true, true, true, true, ptExtrap);
250 
251  // Grab analysis contours to find bin center
252  std::string p = ""; //delta && !corrSysts
253  // if(plot == "ssth23dmsq32" && corrSysts) p = "ssth23dmsq32_";
254  // if(plot == "deltassth23" && corrSysts) p = "deltassth23_";
255  if(plot == "ssth23dmsq32") p = "dmsq_";
256  if(plot == "deltath13") p = "th13_";
257 
258  // std::string helperupsname = std::string(getenv("FCHELPERANA2019_LIB_PATH"));
259  std::string helperupsname = "/nova/app/users/scalvez/fchelper/";
260  TString anasurfname = helperupsname
261  + "/Surfaces/hist_contours_2019_joint_realData_"+beam+"_only";
262  anasurfname += nh? "NH_" : "IH_";
263  anasurfname += p;
264  anasurfname += corrSysts? "systs.root": "stats.root";
265 
266  TFile *anasurffile = TFile::Open(anasurfname);
267  anasurffile->cd();
268  std::string plotShortName="delta_th23";
269  if(plot == "ssth23dmsq32") plotShortName = "th23_dm32";
270  if(nh) plotShortName += "_NH";
271  else plotShortName += "_IH";
272  TH2* anasurf = (TH2F*)anasurffile->Get(plotShortName.c_str());
273 
274  // Grab the FCHelper2018 file
275  //std::string helperupsname = std::string(getenv("FCHELPERANA2017_LIB_PATH"));
276  TString helpername =
277  helperupsname+"FCInputs/contours_FCInput_2019_joint_realData_"+beam+"_only";
278  helpername += nh? "NH_" : "IH_";
279  helpername += p;
280  helpername += corrSysts? "systs.root":"stats.root";
281 
282  TFile *fchelp = TFile::Open(helpername);
283 
284  // Find the bin center for this job, probably should be arguments
285  // dCP / pi
286 
287  // The surfaces are made with ana::ExpandedHistrogram() which sets the bin edges so that
288  // the input min/max x and y fall in the center of the bins. Grab these boundaries.
289  int NX = anasurf->GetXaxis()->GetNbins();
290  int NY = anasurf->GetYaxis()->GetNbins();
291  double minX = anasurf->GetXaxis()->GetBinCenter(1); // root bins start at 1
292  double maxX = anasurf->GetXaxis()->GetBinCenter(NX);
293  double minY = anasurf->GetYaxis()->GetBinCenter(1);
294  double maxY = anasurf->GetYaxis()->GetBinCenter(NY);
295 
296  double widthX = anasurf->GetXaxis()->GetBinWidth(1);
297  double widthY = anasurf->GetYaxis()->GetBinWidth(1);
298  std::cout << "widthX = " << widthX << endl;
299  std::cout << "widthY = " << widthY << endl;
300  assert(fbin < NX*NY && "Invalid bin number");
301 
302  int bin = 0;
303  for (int ibin = startidx; ibin <= endidx; ibin++) {
304  bin = bins[ibin];
305  std::cout<<"- "<<ibin<<std::endl;
306  std::cout << "bin " << bin << endl;
307  // Once we run out of important bins, get out.
308  if(bins[ibin] == 0 && ibin != 0) continue;
309 
310  // bin = bins[i];
311 
312  int binY = bin/NX;
313  int binX = bin - NX*binY;
314  std::cout << "binX: " << binX << endl;
315  std::cout << "binY: " << binY << endl;
316  anasurffile->cd();
317  double centerX = anasurf->GetXaxis()->GetBinCenter(binX+1);
318  double centerY = anasurf->GetYaxis()->GetBinCenter(binY+1);
319  std::cout << "centerX = " << centerX << endl;
320  std::cout << "centerY = " << centerY << endl;
321  fchelp->cd();
322 
323  // Pull out hists and get values at current bin
324  double seedssth23 = 0;
325  if(plot == "deltath13"){
326  TH2* hssth23 = (TH2F*) fchelp->Get((hierStr+"_SinSqTheta23").c_str());
327  seedssth23 = hssth23->GetBinContent(binX+1,binY+1);
328  delete hssth23;
329  }
330  double seeddmsq32 = 0;
331  if (plot == "deltassth23"){
332  TH2* hdmsq = (TH2F*)fchelp->Get((hierStr+"_DmSq32").c_str());
333  seeddmsq32 = hdmsq->GetBinContent(binX+1,binY+1); // root bins off-by-one
334  delete hdmsq;
335  }
336  double seeddelta = 0;
337  if (plot == "ssth23dmsq32"){
338  TH2* hdelta = (TH2F*)fchelp->Get((hierStr+"_DeltaCPpi").c_str());
339  seeddelta = hdelta->GetBinContent(binX+1,binY+1);
340  delete hdelta;
341  }
342  TH2* hss2th13 = (TH2F*)fchelp->Get((hierStr+"_SinSq2Theta13").c_str());
343  double seedss2th13 = hss2th13->GetBinContent(binX+1,binY+1);
344  delete hss2th13;
345 
346  // Get current syst shifts
347  std::map<std::string,double> seedsyst;
348  for (const ISyst* syst : systs){
349  TH2* h = (TH2F*)fchelp->Get((hierStr+"_"+syst->ShortName()).c_str());
350  if (!h){
351  std::cout << "Don't see a prof hist for " << syst->ShortName() << ". Continuing, but check your ups version." << std::endl;
352  continue;
353  }
354  double seedval = h->GetBinContent(binX+1,binY+1);
355  seedsyst.insert(std::pair<std::string,double>(syst->ShortName(),seedval));
356  }
357 
358  std::cout << "Parameters from anasurf " << endl;
359  std::cout << "ssth23 = " << seedssth23 << endl;
360  std::cout << "dmsq = " << seeddmsq32 << endl;
361  std::cout << "delta = " << seeddelta << endl;
362  std::cout << "th13 = " << seedss2th13 << endl;
363 
364  // Need random numbers to sample the bin
365  TRandom3 rnd(0);
366 
367  for (int i = 0; i < NPts; i++){
368  ResetOscCalcToDefault(calc);
369 
370  const IFitVar* dmsq_var;
371  if(nh) dmsq_var = &kFitDmSq32ScaledNH;
372  else dmsq_var = &kFitDmSq32ScaledIH;
373 
374  // backwards compatability
375  double trueX = centerX;
376  double trueY = centerY;
377  std::cout << "trueX: " << trueX << endl;
378  std::cout << "trueY: " << trueY << endl;
379  // Plot = deltassth23: x=delta, y=ssth23
380  // Plot = ssth23dmsq: x=ssth23, y=dmsq32
381  // Plot = deltassth13: x=th13, y=delta
382  if(plot != "deltath13"){
383  kFitDeltaInPiUnits.SetValue(calc, plot=="deltassth23"?trueX:seeddelta);
384  kFitSinSqTheta23.SetValue(calc, plot=="deltassth23"?trueY:trueX);
385  dmsq_var->SetValue(calc, plot=="deltassth23"?seeddmsq32:trueY);
386  kFitSinSq2Theta13.SetValue(calc, seedss2th13);
387  }
388  else{
389  kFitDeltaInPiUnits.SetValue(calc, trueY);
390  kFitSinSqTheta23.SetValue(calc, seedssth23);
391  dmsq_var->SetValue(calc, seeddmsq32);
392  kFitSinSq2Theta13.SetValue(calc, trueX);
393  }
394 
395  // Need to seed psuedo-experiments at profiled syst shifts read from
396  // FCInput
397  SystShifts seednue_fhc;
398  for (const ISyst *syst : snue_fhc){
399  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
400  seednue_fhc.SetShift(syst,seedsyst[syst->ShortName()]);
401  }
402  SystShifts seednue_rhc;
403  for (const ISyst *syst : snue_rhc) {
404  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
405  seednue_rhc.SetShift(syst, seedsyst[syst->ShortName()]);
406 
407  }
408  SystShifts seednumu;
409  for (const ISyst *syst : snumu){
410  if (seedsyst.find(syst->ShortName()) == seedsyst.end()) continue;
411  seednumu.SetShift(syst,seedsyst[syst->ShortName()]);
412  }
413 
414 
415  // Make a mock data spectrum for nue
416  // Use mock data + cosmic as predictions
417 
418  std::cout << "----------Making mock data at-----------------\n";
419  std::cout << "ssth23 = " << kFitSinSqTheta23.GetValue(calc);
420  std::cout << ", delta = " << kFitDeltaInPiUnits.GetValue(calc);
421  std::cout << ", dmsq = " << dmsq_var->GetValue(calc);
422  std::cout << ", th13 = " << kFitSinSq2Theta13.GetValue(calc);
423  std::cout << endl;
424 
425  // std::vector<Spectrum> mockdata;
426  // std::vector<Spectrum> mockdata_numu;
427 
428  // Vestige?
429  // SystShifts this_seednue;
430 
431  // for (auto ipred = 0; ipred < preds.size(); ipred++) {
432 
433  // // SystShifts tmp_systs;
434  // // if(preds[ipred].flavor == fc::kNue) {
435  // // if(preds[ipred].beam == fc::kFHC) tmp_systs = seednue_fhc;
436  // // else tmp_systs = seednue_rhc;
437  // // }
438  // // else tmp_systs = seednumu;
439  // std::cout << "___debug ipred " << ipred << std::endl;
440  // Spectrum fakedata = preds[ipred].pred->Predict (calc ); // stats only
441  // // if(corrSysts) fakedata = preds[ipred].pred->PredictSyst(calc, tmp_systs); // with systs
442  // std::cout << "___debug cos first " << preds[ipred].cos.first << std::endl;
443 
444  // fakedata += *preds[ipred].cos.first;
445  // std::cout << "___debug mock pot " << preds[ipred].pot << std::endl;
446  // mockdata.push_back(fakedata.MockData(preds[ipred].pot));
447  // }
448 
449  std::vector<const IExperiment*> expts;
450  std::vector<const IExperiment*> numu_expts;
451 
452  int nnumu = 4;
453  // auto calc = DefaultOscCalc();
454  //if(fake2019) SetFakeCalc(calc, 0.564956, 2.48304e-3, 1.9997);
455  // if(fake2019) SetFakeCalc(calc, 0.565, 2.48e-3, 0);
456  // else if(!realData) {std::cerr << "need setting for data\n"; exit(1);}
457 
458  for(int ipred = 0; ipred < int(preds.size()); ++ipred) {
459  double POT = preds[ipred].pot;
460  // auto fakedata = GetFakeData (preds[ipred].pred, calc, POT, preds[ipred].cos.first, preds[ipred].livetime);
461  auto fakedata = preds[ipred].pred->Predict(calc);
462  auto hcos = preds[ipred].cos.first->ToTH1(preds[i].livetime, kBlack, kSolid, kLivetime);
463  fakedata += Spectrum(hcos,POT,0);
464  // cout<<i<<" "<<preds[i].name<<" POT "<<POT<<" tot MC "<<preds[i].pred->Predict(calc).Integral(POT)<<
465  // " cos "<<hcos->Integral()<<" cos er "<<preds[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
466  auto mockdata = fakedata.MockData(POT);
467  expts.push_back(new SingleSampleExperiment(preds[ipred].pred, mockdata, *preds[ipred].cos.first, preds[ipred].cos.second));
468  }
469 
470  for(int ipred = 0; ipred < (int) preds_numu_only.size(); ++ipred){
471  double POT = preds_numu_only[ipred].pot;
472  // auto fakedata = GetFakeData (preds_numu_only[ipred].pred, calc, POT, preds_numu_only[ipred].cos.first);
473  auto fakedata = preds_numu_only[ipred].pred->Predict(calc);
474  auto hcos = preds[ipred].cos.first->ToTH1(preds[i].livetime, kBlack, kSolid, kLivetime);
475  // cout<<i<<" "<<preds_numu_only[i].name<<" POT "<<POT<<" tot MC "<<preds_numu_only[i].pred->Predict(calc).Integral(POT)<<
476  // " cos "<<hcos->Integral()<<" cos er "<<preds_numu_only[i].cos.second<<" analyze data "<<thisdata->Integral(POT)<<endl;
477  fakedata += Spectrum(hcos,POT,0);
478  auto mockdata = fakedata.MockData(POT);
479  numu_expts.push_back(new SingleSampleExperiment(preds_numu_only[ipred].pred, mockdata, *preds_numu_only[ipred].cos.first, preds_numu_only[ipred].cos.second));
480  }
481 
482  // Set up experiment, nova nue + nova numu + 2017 reactor
483  // Get true chi-square before fit
484 
485  if (plot != "th13") {
486  expts .push_back( WorldReactorConstraint2019() );
487  numu_expts.push_back( WorldReactorConstraint2019() );
488  }
489  std::cout << "Creating multiexperiment from " << expts.size()
490  << " experiments" << std::endl;
491 
492  MultiExperiment expt(expts);
493  MultiExperiment numu_expt(numu_expts);
494 
495 
496  // Needed?
497  std::cout<<"Clean verctors"<<std::endl;
498  preds.clear();
499  preds_numu_only.clear();
500 
501  // int nnumu = 4;
502  std::cout << "\n\nSystematic correlations...\n";
503  if(!nueOnly && ! numuOnly && corrSysts){
504  if(FHCOnly){
505  expt.SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
506  auto notfornumu = GetCorrelations(false, true, ptExtrap);
507  for(int i =0; i < nnumu; ++i) expt.SetSystCorrelations(i+1, notfornumu);
508  }
509  if(RHCOnly){
510  expt.SetSystCorrelations(0, GetCorrelations(true, false, ptExtrap));
511  auto notfornumu = GetCorrelations(false, true, ptExtrap);
512  for(int i =0; i < nnumu; ++i) expt.SetSystCorrelations(i+1, notfornumu);
513  }
514  if(both){
515  expt.SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
516  expt.SetSystCorrelations(1, GetCorrelations(true, false, ptExtrap));
517  auto notfornumu = GetCorrelations(false, true, ptExtrap);
518  for(int i =0; i < 8; ++i) expt.SetSystCorrelations(i+2, notfornumu);
519  }
520  }
521  if (nueOnly && corrSysts){
522  if (FHCOnly) expt.SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
523  if (RHCOnly) expt.SetSystCorrelations(0, GetCorrelations(true, false, ptExtrap));
524  if (both) {
525  expt.SetSystCorrelations(0, GetCorrelations(true, true, ptExtrap));
526  expt.SetSystCorrelations(1, GetCorrelations(true, false, ptExtrap));
527  }
528  }
529 
530 
531  // To fix when systs are available
532  // MultiExperiment expt(expts);
533  // MultiExperiment numu_expt(numu_expts);
534  // // Correlate your systs
535  // // TODO GetCorrelation function calls are different
536  // for(int j = 0 ; j < (int)nue_preds.size(); j++){
537  // if(j==0 && !RHCOnly)
538  // expt.SetSystCorrelations(j, GetCorrelations(true, true)); // nue FHC
539  // else
540  // expt.SetSystCorrelations(j, GetCorrelations(true, false)); // nue RHC
541  // auto notfornumu = GetCorrelations(false, false); // numu not dependent on beam
542  // for(int j = 0; j < (int)numu_preds.size(); ++j)
543  // expt.SetSystCorrelations(j+nue_preds.size(), notfornumu);
544  // }
545 
546 
547 
548  /////////////////////////////////////////////////////////////////
549  // Numu only fit to figure out dmsq seeds
550  // (may be able to do without this for contours)
551  /////////////////////////////////////////////////////////////////
552  std::cout << "Fitting numu to find pre seeds" << endl;
553 
554  double maxmixing = .514;
556  // osc::IOscCalcAdjustable* numu_calc = DefaultOscCalc();
557  SystShifts auxShifts = SystShifts::Nominal();
558  ResetOscCalcToDefault(calc);
559 
560  MinuitFitter fitnumu_only(&numu_expt, {&kFitSinSqTheta23, &kFitDmSq32});
561 
562  double minchi_numu =
563  fitnumu_only.Fit(calc, auxShifts,
564  SeedList({{&kFitSinSqTheta23, {0.4}},
565  {&kFitDmSq32, {nh? 2.35e-3 : -2.35e-3}}}),
566  {},
567  IFitter::kQuiet);
568 
569  double pre_seed_dmsq = dmsq_var->GetValue(calc);
570  double pre_seed_th23 = kFitSinSqTheta23.GetValue(calc);
571 
572  double pre_seed_th23_LO = ( pre_seed_th23 > maxmixing ) ?
573  (2*maxmixing - pre_seed_th23) : pre_seed_th23;
574  double pre_seed_th23_UO = ( pre_seed_th23 > maxmixing ) ?
575  pre_seed_th23 : (2*maxmixing - pre_seed_th23);
576  if(pre_seed_th23_LO > .5) pre_seed_th23_LO = .45;
577 
578  std::cout << "Preseeds ---> ssth23_LO = " << pre_seed_th23_LO;
579  std::cout << "\tdmsq_NH = " << fabs(pre_seed_dmsq) << endl;
580  // done getting pre seeds
581  // --------------------------------------------------------------
582  ////////////////////////////////////////////////////////////////////
583  // True fit
584  ////////////////////////////////////////////////////////////////////
585 
586  std::vector<const IFitVar*> trueVars;
587  if (plot=="deltassth23"){
588  trueVars.push_back(dmsq_var);
589  trueVars.push_back(&kFitSinSq2Theta13);
590  }
591  if (plot=="ssth23dmsq32"){
592  trueVars.push_back(&kFitDeltaInPiUnits);
593  trueVars.push_back(&kFitSinSq2Theta13);
594  }
595 
596 
597  MinuitFitter fittrue(&expt, trueVars, systs);
598  double tb_th13 = 0;
599  double tb_dmsq = 0;
600  double tb_ssth23 = 0;
601  double tb_delta = 0;
602  SystShifts systshift = {};
603  // if not correcting deltassth23, then put one element in this
604  // vector so the following loop will run
605  std::vector <double> tf_delta_seeds = {0, .5, 1., 1.5};
606  std::vector <double> tf_ssth23_seeds = {pre_seed_th23_LO,
607  maxmixing,
608  pre_seed_th23_UO};
609 
610  // make a map of all of the seeds for the true fit
611  std::map<const IFitVar*, std::vector<double>> tf_seedPts;
612  if(plot == "deltassth23")
613  tf_seedPts = {{dmsq_var, {nh? fabs(pre_seed_dmsq) : -1*fabs(pre_seed_dmsq)}},
614  {&kFitSinSq2Theta13, {0.082}}};
615  else if(plot == "ssth23dmsq32")
616  tf_seedPts = {{&kFitDeltaInPiUnits, tf_delta_seeds},
617  {&kFitSinSq2Theta13, {0.082}}};
618  else if(plot == "deltath13")
619  tf_seedPts = {{dmsq_var, {nh? fabs(pre_seed_dmsq) : -1*fabs(pre_seed_dmsq)}},
620  {&kFitSinSqTheta23, tf_ssth23_seeds}};
621 
622  std::cout << "True fit parameters (pre fit): ssth23 = " << kFitSinSqTheta23.GetValue(calc);
623  std::cout << ", delta = " << kFitDeltaInPiUnits.GetValue(calc);
624  std::cout << ", ss2th13 = " << kFitSinSq2Theta13.GetValue(calc);
625  std::cout << ", dmsq = " << dmsq_var->GetValue(calc) << " x10^-3 eV" << endl;
626 
627 
628  std::cout << "Starting true fit --->" << endl;
629  auxShifts.ResetToNominal();
630  double trueChi = fittrue.Fit(calc,
631  auxShifts,
632  tf_seedPts,
633  {},
634  IFitter::kQuiet);
635 
636  // save best parameters for posterity
637  tb_delta = kFitDeltaInPiUnits.GetValue(calc);
638  tb_ssth23 = kFitSinSqTheta23.GetValue(calc);
639  tb_dmsq = dmsq_var->GetValue(calc);
640  tb_th13 = kFitSinSq2Theta13.GetValue(calc);
641 
642  std::cout << "True fit parameters (post fit): ssth23 = " << tb_ssth23;
643  std::cout << ", delta = " << tb_delta;
644  std::cout << ", ss2th13 = " << tb_th13;
645  std::cout << ", dmsq = " << tb_dmsq << " x10^-3 eV" << endl;
646 
647  // Set up best fitter
648  MinuitFitter fit(&expt,
651  systs);
652 
653  std::vector <double> bf_dmsq_seeds = {fabs(pre_seed_dmsq),
654  -1*fabs(pre_seed_dmsq)};
655  std::vector <double> bf_delta_seeds = {0., 0.5, 1., 1.5};
656  std::vector <double> bf_th23_seeds = {pre_seed_th23_LO, pre_seed_th23_UO,
657  maxmixing};
658  std::vector <double> bf_th13_seeds = {0.082};
659  std::map <const IFitVar*, std::vector<double>> bf_seedPts =
660  {{&kFitDeltaInPiUnits, bf_delta_seeds},
661  {&kFitSinSqTheta23, bf_th23_seeds},
662  {&kFitSinSq2Theta13, bf_th13_seeds},
663  {&kFitDmSq32Scaled, bf_dmsq_seeds}};
664 
665  double bestdelta = 0;
666  double bestssth23 = 0;
667  double bestth13 = 0;
668  double bestdmsq = 0;
669  std::cout << "Starting best fit ----> " << endl;
670  ResetOscCalcToDefault(calc);
671  auxShifts.ResetToNominal();
672  double bestChi = fit.Fit(calc,
673  auxShifts,
674  bf_seedPts,
675  {},
676  IFitter::kQuiet);
677 
678  bestdelta = kFitDeltaInPiUnits.GetValue(calc);
679  bestssth23 = kFitSinSqTheta23.GetValue(calc);
680  bestth13 = kFitSinSq2Theta13.GetValue(calc);
681  bestdmsq = kFitDmSq32Scaled.GetValue(calc);
682  std::cout << "Best fit parameters: ssth23 = " << bestssth23;
683  std::cout << ", delta = " << bestdelta;
684  std::cout << ", ss2th13 = " << bestth13;
685  std::cout << ", dmsq = " << bestdmsq << " x10^-3 eV"<< endl;
686 
687 
688  // Make pseudoexperiment dst
689  // with FCPoint in its current state (r31771) the names of the branches are
690  // not correct. We want to save best fit information
691  FCPoint pt(seeddelta, seedssth23, seedss2th13, seeddmsq32,
692  pre_seed_th23, pre_seed_dmsq,
693  tb_delta, tb_ssth23, tb_th13, tb_dmsq,
694  trueChi,
695  bestdelta, bestssth23, bestth13, bestdmsq,
696  bestChi,
697  0, 0);
698 
699  fccol.AddPoint(pt);
700 
701  unsigned int nfc_pts = fccol.NPoints();
702  if(npts_flush != 0) {
703  if(nfc_pts % npts_flush == 0) {
704  std::cout << "Saving " + std::to_string(nfc_pts) + "th point to " << fcout << endl;
705  fccol.SaveToFile(fcout);
706  }
707  }
708  }
709  }// end loop over bins
710  // Save to file
711  fccol.SaveToFile(fcout);
712 }
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
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool intersection=true)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var&#39;s GetValue()
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
const FitDmSq32ScaledNH kFitDmSq32ScaledNH
Definition: FitVars.cxx:21
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
const char * p
Definition: xmltok.h:285
void AddPoint(const FCPoint &p)
Definition: FCCollection.h:17
OStream cerr
Definition: OStream.cxx:7
double maxY
Definition: plot_hist.C:10
std::vector< predictions > LoadPredictions(bool corrSysts=false, bool runOnGrid=false, std::string decomp="", bool nue=true, bool numu=true, bool fhc=true, bool rhc=true, bool NERSC=false)
void make_fc_surfaces_2020(int NPts, int startidx, int endidx, bool nh, int N, int jid, int npts_flush=0, std::string plot="deltassth23", std::string options="joint_fake2019", std::string decomp="", TString beam="both", bool corrSysts=false, bool readFromFile=false, bool runOnGrid=false)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:48
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
unsigned int NPoints() const
Definition: FCCollection.h:26
Represents the results of a single Feldman-Cousins pseudo-experiment.
Definition: FCPoint.h:8
expt
Definition: demo5.py:34
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
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 double j
Definition: BetheBloch.cxx:29
const Binning bins
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
void SetSeeds(osc::IOscCalcAdjustable *calc, double delta, double th23, double dmsq)
float bin[41]
Definition: plottest35.C:14
std::vector< float > Spectrum
Definition: Constants.h:610
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
static bool isFHC
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
const std::string path
Definition: plot_BEN.C:43
Combine multiple component experiments.
const FitDmSq32ScaledIH kFitDmSq32ScaledIH
Definition: FitVars.cxx:22
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const =0
double livetime
Definition: saveFDMCHists.C:21
void ResetToNominal()
Definition: SystShifts.cxx:144
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
exit(0)
T cos(T number)
Definition: d0nt_math.hpp:78
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
Definition: IFitVar.h:16
const ReactorExperiment * WorldReactorConstraint2019()
Reactor constraint from PDG 2019.
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue)
const FitSinSq2Theta13 kFitSinSq2Theta13
Definition: FitVars.cxx:13
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
void SaveToFile(const std::string &fname) const
Float_t e
Definition: plot.C:35
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
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)
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
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const =0
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string