sensitivity2018.C
Go to the documentation of this file.
8 #include "CAFAna/Core/Spectrum.h"
9 #include "CAFAna/FC/FCSurface.h"
14 #include "CAFAna/Core/Spectrum.h"
15 #include "CAFAna/Core/Registry.h"
20 #include "CAFAna/Systs/Systs.h"
22 #include "CAFAna/Vars/FitVars.h"
24 #include "OscLib/IOscCalc.h"
25 #include "Utilities/rootlogon.C"
26 
27 #include "TFile.h"
28 #include "TStyle.h"
29 #include "TGraph.h"
30 #include "TCanvas.h"
31 #include "TH2.h"
32 
33 using namespace ana;
34 
35 
36 void sensitivity2018(bool MakeFile = false, std::string ana = "joint", bool isSyst = false, bool MakeSurf = false, std::string whatsurf = "th23dcp", bool isNH = true, bool MakeSens = false, std::string whatplot = "mh")
37 {
38 
39 std::string version = "v13_";
40 std::string outfilename = "run_sensitivity_"+version;
41 
42 // ana=="joint", "fhc"
43 // isSyst - use syst, not only nominals
44 // whatplot: {"mh", "dcp", "oct", "maxmix"}
45 // MakeFile - write to file or not, in case of not - read from file and plot
46 // whatsurf: th23dcp dmsqth23
47 // isNH true false
48 
49 ////////////////////////////////////////////////////////////////
50 /////////////////////// Load nue and numu stuff ///////////////
51 ///////////////////////////////////////////////////////////////
52 
53  cout<<"Greetings!"<<endl;
54  cout<<"You requested: "<<" with "<<(isSyst?" systematic ones ":" nominals ")<<(MakeFile?" make files ":" plot ")<<endl;
55  cout<<"It will be "<<ana<<" analysis for nue+numu, plot: "<<(MakeSurf?"contours ":"no contours ")<<(isNH?" NH ":" IH ")<<" ; "<<whatsurf<<(MakeSens?" slices ":" no slices")<<" in case of slices it will be "<<whatplot<<endl;
56 
57 
58  std::string SystTag = isSyst?"with_syst_":"just_stat_";
59 
60  std::string TagName = SystTag+ana;
61 
63  calc->SetdCP(1.21*M_PI);
64  calc->SetTh23(asin(sqrt(0.558)));
65  calc->SetDmsq32(2.44e-3);
66 
67  std::vector<double> pot;
68  std::vector<double> livetime;
69 
70  const double potFD_rhc = kAna2018RHCPOT;
71  const double potFD_fhc = kAna2018FHCPOT;
72 
73  const double cosFD_rhc = kAna2018RHCLivetime;
74  const double cosFD_fhc = kAna2018FHCLivetime;
75 
76 
77  std::vector<const IPrediction*> predsvector;
78  std::vector<std::pair <TH1D*, double>> cosmicsvector;
79 
80  std::vector<const IPrediction*> predictionsNue;
81  std::vector<const IPrediction*> predictionsNumu;
82  std::vector<std::pair <TH1D*, double>> cosmicsNue;
83  std::vector<std::pair <TH1D*, double>> cosmicsNumu;
84 
85 
86  string decomp_rhc = "prop";
87  string decomp_fhc = "combo"; //should be combo
88 
89  cosmicsNue.push_back(GetNueCosmics2018("fhc"));
90  if(ana=="joint") cosmicsNue.push_back(GetNueCosmics2018("rhc"));
91 
92  int nnumu = 4;
93 
94  auto nue_preds_fhc = GetNuePrediction2018(decomp_fhc, DefaultOscCalc(), isSyst, "fhc", false);
95  auto nue_preds_rhc = GetNuePrediction2018(decomp_rhc, DefaultOscCalc(), isSyst, "rhc", false);
96 
97  predictionsNue.push_back(nue_preds_fhc);
98  if(ana=="joint") predictionsNue.push_back(nue_preds_rhc);
99 
100  auto numu_preds_fhc = GetNumuPredictions2018(nnumu, isSyst, "fhc");
101  auto numu_preds_rhc = GetNumuPredictions2018(nnumu, isSyst, "rhc");
102 
103  predictionsNumu.insert(predictionsNumu.end(), numu_preds_fhc.begin(), numu_preds_fhc.end());
104  if(ana=="joint") predictionsNumu.insert(predictionsNumu.end(), numu_preds_rhc.begin(), numu_preds_rhc.end());
105 
106  predsvector.insert(predsvector.end(), predictionsNue.begin(), predictionsNue.end());
107  predsvector.insert(predsvector.end(), predictionsNumu.begin(), predictionsNumu.end());
108 
109  auto numu_cosmics_fhc = GetNumuCosmics2018(nnumu, "fhc");
110  auto numu_cosmics_rhc = GetNumuCosmics2018(nnumu, "rhc");
111 
112  cosmicsNumu.insert(cosmicsNumu.end(),numu_cosmics_fhc.begin(), numu_cosmics_fhc.end());
113  if(ana=="joint") cosmicsNumu.insert(cosmicsNumu.end(),numu_cosmics_rhc.begin(), numu_cosmics_rhc.end());
114 
115  cosmicsvector.insert(cosmicsvector.end(), cosmicsNue.begin(), cosmicsNue.end());
116  cosmicsvector.insert(cosmicsvector.end(), cosmicsNumu.begin(), cosmicsNumu.end());
117 
118 
119  cout<<"Good job"<<endl;
120 
121 
122  //////////////////////////////////////////////////////////////
123  ///////////////////// Check Predictions at first /////////////
124  /////////////////////////////////////////////////////////////
125 
126 
127  pot.push_back(potFD_fhc);
128  pot.push_back(potFD_rhc);
129 
130  livetime.push_back(cosFD_fhc);
131  livetime.push_back(cosFD_rhc);
132 
133  std::vector<TH1*> hNue; // nue signal
134  std::vector<TH1*> hMC; //total MC
135  std::vector<TH1*> hTotbkg; //total sum of all bkg
136  std::vector<TH1*> hBkg; // all bkg except cosmics
137  std::vector<TH1*> hNC; //1
138  std::vector<TH1*> hCC; //2
139  std::vector<TH1*> hBeam; //3
140  std::vector<TH1*> hTau; //4
141  std::vector<TH1*> hCos; //5
142  std::vector<TH1*> hRock; //6 rock events
143  std::vector<TH1*> hNue_bg; //7 wrong sign nue
144 
145  std::vector<TH1*> hNumu;
146  std::vector<TH1*> hCosNumu;
147  std::vector<TH1*> hMCNumu;
148 
149  for(int i = 0; i < (int)predictionsNue.size(); i++){
150 
151  if (i==0){
152  hNue.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu).ToTH1(pot[i])); //signal FHC
153  hNue_bg.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot[i])); //wrong sign bkg for FHC
154  }
155  if (i==1){
156  hNue.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kAntiNu).ToTH1(pot[i])); //signal RHC
157  hNue_bg.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuMuToNuE, Current::kCC, Sign::kNu).ToTH1(pot[i])); //wrong sign bkg for RHC
158  }
159 
160  hTau.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kAllNuTau, Current::kCC, Sign::kBoth).ToTH1(pot[i])); // tau CC bkg (2)
161  hBeam.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kNuEToNuE, Current::kCC, Sign::kBoth).ToTH1(pot[i])); // beam bkg (3)
162  hNC.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kAll, Current::kNC, Sign::kBoth).ToTH1(pot[i])); // NC bkg (4)
163  hCC.push_back(predictionsNue[i]->PredictComponent(calc, Flavors::kAllNuMu, Current::kCC, Sign::kBoth).ToTH1(pot[i])); // numu CC bkg (5)
164  hMC.push_back(predictionsNue[i]->Predict(calc).ToTH1(pot[i]));
165  TH1* temp = (TH1*)hNue_bg[i]->Clone(UniqueName().c_str());
166  temp->Add(hTau[i]);
167  temp->Add(hBeam[i]);
168  temp->Add(hNC[i]);
169  temp->Add(hCC[i]);
170  hBkg.push_back((TH1*)temp->Clone(UniqueName().c_str()));
171  hCos.push_back((TH1D*)cosmicsNue[i].first->Clone(UniqueName().c_str())); //cosmic bkg (6)
172  temp->Add(hCos[i]);
173  hMC[i]->Add(hCos[i]);
174  hTotbkg.push_back(temp);
175  }
176 
177  cout<<">>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>"<<endl;
178  cout<<"Here are the predictions for 2018 nue ana with 2017 Best fit values: "<<endl;
179 
180  for(int i = 0; i < (int)predictionsNue.size(); i++){
181 
182  if (i==0) {cout<<" FHC mode:"<<endl;}
183  if (i==1) {cout<<" RHC mode:"<<endl;}
184 
185  cout<<"With "<<pot[i]<<" POT "<<livetime[i]<<" s livetime"<<endl;
186  cout<<"total MC "<<hMC[i]->Integral()<<" signal "<<hNue[i]->Integral()<<" total bkg (beam+cosmic) "<<hTotbkg[i]->Integral()<<endl;
187  cout<<"total beam bkg "<<hBkg[i]->Integral()<<" wrong sign bkg "<<hNue_bg[i]->Integral()<<" beam nue "<<hBeam[i]->Integral()<<" NC "<<hNC[i]->Integral()<<endl;
188  cout<<"numu CC "<<hCC[i]->Integral()<<" tau CC "<<hTau[i]->Integral()<<" cosmic "<<hCos[i]->Integral()<<endl;
189 
190  }
191 
192 
193  for(int i = 0; i < (int)predictionsNumu.size(); i++){
194  double pot, livetime;
195  if(i<4) {pot=potFD_fhc; livetime = cosFD_fhc;}
196  else {pot=potFD_rhc; livetime = cosFD_rhc;}
197 
198  hNumu.push_back(predictionsNumu[i]->PredictComponent(calc, Flavors::kNuMuToNuMu, Current::kCC, Sign::kBoth).ToTH1(pot));
199  hMCNumu.push_back(predictionsNumu[i]->Predict(calc).ToTH1(pot));
200  hCosNumu.push_back((TH1D*)cosmicsNumu[i].first->Clone(UniqueName().c_str()));
201 
202  }
203 
204  cout<<"Here are the predictions for 2018 numu ana with 2017 Best fit values, for checks: "<<endl;
205  for(int i = 0; i < (int)predictionsNumu.size(); i++){
206 
207  if (i<4) {cout<<" FHC mode:"<<endl;}
208  else {cout<<" RHC mode:"<<endl;}
209 
210  cout<<"total MC "<<hMCNumu[i]->Integral()<<" signal "<<hNumu[i]->Integral()<<" cosmic "<<hCosNumu[i]->Integral()<<endl;
211 
212  }
213 
214  /////////////////////////////////////////////////////////////
215  /////////////////// Sensitivities ///////////////////////////
216  ////////////////////////////////////////////////////////////
217 
218 
219  // TODO Extrap, Extrap + Syst; also nue, numu, joint; RHC, FHC, joint
220  // species from 2017 TN: th23 vs delta; th23 vs dmsq23; slices: hie, dcp, octant, max mixing
221 
222 
223  //make fakedata at 2017 BF(make choosable)
224 
225  ResetOscCalcToDefault(calc);
226  calc->SetdCP(1.21*M_PI);
227  calc->SetTh23(asin(sqrt(0.558)));
228  calc->SetDmsq32(2.44e-3);
229 
230  std::vector <const IExperiment * > expts;
231  cout<<"Fake data with systematic pred. and multiexp"<<endl;
232  for(int i = 0; i < (int) predsvector.size(); ++i){
233  double pot;
234  if(ana=="joint" && (i==1 || i>5)) pot=potFD_rhc; //sigh
235  else pot=potFD_fhc;
236  cout<<pot<<endl;
237 
238  auto thisdata = GetFakeData (predsvector[i],calc, pot, cosmicsvector[i].first);
239  cout<<"fake data: "<<thisdata->Integral(pot)<<endl;
240  expts.push_back(new SingleSampleExperiment(predsvector[i], *thisdata, cosmicsvector[i].first, cosmicsvector[i].second));
241 
242  }
243 
244  expts.push_back(WorldReactorConstraint2017());
245  auto exptsAll = new MultiExperiment(expts);
246  cout<<"Experiments are ready, number of items in vector is "<<expts.size()<<endl;
247 
248  std::vector< const ISyst * > systs;
249  SystShifts auxShifts;
250 
251  if(isSyst){
252  systs = GetJointFitSystematicList(true, false, false, true, true, true);
253  auto notfornueFHC = GetCorrelations(true, true);
254  auto notfornueRHC = GetCorrelations(true, false);
255  auto notfornumu = GetCorrelations(false, true); //it will never enter the FHC, RHC ifs, so doesn't matter
256  if(ana=="fhc"){
257  exptsAll->SetSystCorrelations(0, notfornueFHC);
258  for(int i=1; i<5;i++) {
259  exptsAll->SetSystCorrelations(i, notfornumu);
260  }
261  }
262  if(ana=="joint"){
263  exptsAll->SetSystCorrelations(0, notfornueFHC);
264  exptsAll->SetSystCorrelations(1, notfornueRHC);
265  for(int i=2; i<10;i++) {exptsAll->SetSystCorrelations(i, notfornumu);}
266  }
267 
268  auxShifts = SystShifts::Nominal();
269  //Registry::Print();
270  }
271  else { cout<<"no systematic"<<endl; systs={}; auxShifts = junkShifts;}
272 
273  // find BF here:
274 
275  std::vector <const IFitVar *> fitvars = {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32Scaled, &kFitSinSq2Theta13};
276  MinuitFitter fit(exptsAll, fitvars, systs);
277 
278  double minchi23 = 1E20;
279  ResetOscCalcToDefault(calc);
280  for(double thseed:{0.45,0.55}){
281  for(int hie:{-1,1}){
282  std::cout << "Finding best fit for " << (thseed<0.5?"LO ":"UO ") << (hie>0? "NH " : "IH ") <<endl;
283  auto thisminchi = fit.Fit( calc, auxShifts ,
284  {{&kFitDmSq32, {hie*fabs(calc->GetDmsq32())}},
285  {&kFitSinSqTheta23, {thseed}},
286  {&kFitDeltaInPiUnits, {0,0.5,1,1.5}}},
288  )->EvalMetricVal();
289  minchi23= min(minchi23, thisminchi);
290  ResetOscCalcToDefault(calc);
291  if(isSyst){auxShifts = SystShifts::Nominal(); cout<<"reset syst as well"<<endl;}
292  }//end hie
293  }//end th23
294  std::cout << "Found overall min chi " << minchi23 << std::endl;
295 
296 
297  //contours from 2017 TN
298  if(MakeSurf){
299  if(MakeFile){
300  TFile * outfile = new TFile(("files/"+outfilename+TagName+"_Surfonly_"+whatsurf+"_"+(isNH?"NH":"IH")+".root").c_str(),"recreate");
301  outfile->cd();
302  TVectorD v(1);
303  v[0] = minchi23;
304  v.Write("overall_minchi");
305  int hie;
306  if(isNH) hie = 1;
307  else hie = -1;
308 
309  cout<<"clean verctors"<<endl;
310  predsvector.clear();
311 
312  ResetOscCalcToDefault(calc);
313  if(isSyst) {auxShifts = SystShifts::Nominal();}
314  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
315 
316  int Ysteps = 30; int Xsteps = 30;
317  if (whatsurf=="dmsqth23") {Ysteps = 26; Xsteps = 26; cout<<"minimize with grid "<< Ysteps<<" * "<<Xsteps<<" points"<<endl;}
318  if (whatsurf=="th23dcp") {FrequentistSurface surf_th23_dcp (exptsAll, calc,
319  &kFitDeltaInPiUnits, Xsteps, 0, 2,
320  &kFitSinSqTheta23, Ysteps, 0.2, 0.7,
321  {&kFitDmSq32Scaled, &kFitSinSq2Theta13},
322  systs
323  //{{&kFitSinSq2Theta13, {calc->GetTh13()}}}
324  //seedShifts
325  );
326 
327  outfile->cd();
328  surf_th23_dcp.SaveTo(outfile, (TString)"surf_delta_th23_"+(hie>0?"NH":"IH"));}
329 
330  ResetOscCalcToDefault(calc);
331  if(isSyst) {auxShifts = SystShifts::Nominal();}
332  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
333 
334  if (whatsurf=="dmsqth23") {FrequentistSurface surf_th23_dmsq23 (exptsAll, calc,
335  &kFitSinSqTheta23, Xsteps,0.3,0.7,
336  &kFitDmSq32Scaled, Ysteps, (hie>0?2.1:-2.8), (hie>0?2.8:-2.1),
338  systs,
339  {{&kFitDeltaInPiUnits,{0,0.5,1,1.5}}//, {&kFitSinSq2Theta13, {calc->GetTh13()}}
340  }
341  //seedShifts
342  );
343 
344 
345  outfile->cd();
346  auto surf_test=surf_th23_dmsq23.ToTH2(minchi23);
347  surf_test->Write((TString)"test_th23_dm32_"+(hie>0?"NH":"IH"));
348  surf_th23_dmsq23.SaveTo(outfile, (TString)"surf_th23_dmsq23_"+(hie>0?"NH":"IH"));}
349 
350  outfile->Close();
351  cout<<"Done with surfaces"<<endl;
352  }// end Make File
353 
354  else {
355 
356  TFile * infile = new TFile (("files/"+outfilename+TagName+"_Surfonly_"+whatsurf+".root").c_str(),"read");
357  auto mins =* (TVectorD*)infile->Get("overall_minchi");
358  double minchi23 = mins[0];
359  cout<<"Min chi: "<<minchi23<<endl;
360  int hie;
361  if(isNH) hie = +1;
362  else hie = -1;
363 
364  Color_t k2SigmaConfidenceColor = ((hie>0)? k2SigmaConfidenceColorNH:k2SigmaConfidenceColorIH);
365  Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.70);
366  Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.35);
367  Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.20);
368 
369  auto leg = ContourLegend(hie, false, false, kFillColor1,kFillColor2,kFillColor3,kBlack, false);
370 
371  if (whatsurf=="th23dcp"){
372  auto surf_th23_dcp = *ISurface::LoadFrom(infile, (TString)"surf_delta_th23_" + (hie > 0 ? "NH" : "IH"));
373 
374  TH2 * surf_th23_dcp_1Sigma = Gaussian68Percent2D(surf_th23_dcp);
375  TH2 * surf_th23_dcp_2Sigma = Gaussian2Sigma2D(surf_th23_dcp);
376  TH2 * surf_th23_dcp_3Sigma = Gaussian3Sigma2D(surf_th23_dcp);
377  DrawContourCanvas("delta_th23", 0., 2., 0.225, 0.725);
378  surf_th23_dcp.DrawContour(surf_th23_dcp_1Sigma, kSolid, kFillColor1, minchi23);
379  surf_th23_dcp.DrawContour(surf_th23_dcp_2Sigma, k2SigmaConfidenceStyle, kFillColor2, minchi23);
380  surf_th23_dcp.DrawContour(surf_th23_dcp_3Sigma, kSolid, kFillColor3, minchi23);
381 
382  leg->Draw();
383  Simulation();
384 
385  gPad->Print((TString)"plots/"+version+"surf_th23_dcp_"+TagName+(hie>0 ? "nh" : "ih")+".pdf");
386  }
387 
388  if (whatsurf=="dmsqth23") {
389  auto surf_th23_dmsq23 = *ISurface::LoadFrom(infile, (TString)"surf_th23_dmsq23_" + (hie > 0 ? "NH" : "IH"));
390 
391  TH2 * surf_th23_dmsq23_1Sigma = Gaussian68Percent2D(surf_th23_dmsq23);
392  TH2 * surf_th23_dmsq23_2Sigma = Gaussian2Sigma2D(surf_th23_dmsq23);
393  TH2 * surf_th23_dmsq23_3Sigma = Gaussian3Sigma2D(surf_th23_dmsq23);
394 
395  if (hie>0) DrawContourCanvas("th23_dm32", 0.3, 0.7, 2.05, 2.75);
396  else DrawContourCanvas("th23_dm32", 0.3, 0.7, -2.95, -2.25);
397 
398  surf_th23_dmsq23.DrawContour(surf_th23_dmsq23_1Sigma, kSolid, kFillColor1, minchi23);
399  surf_th23_dmsq23.DrawContour(surf_th23_dmsq23_2Sigma, k2SigmaConfidenceStyle, kFillColor2, minchi23);
400  surf_th23_dmsq23.DrawContour(surf_th23_dmsq23_3Sigma, kSolid, kFillColor3, minchi23);
401 
402  leg->Draw();
403  Simulation();
404  gPad->Print((TString)"plots/"+version+"surf_th23_dmsq23_"+TagName+(hie>0 ? "nh" : "ih")+".pdf");
405  }
406 
407  } // end ! MakeFile = plot
408  } //MakeSurf end
409 
410  if(MakeSens){
411 
412 
413  std::vector<std::pair<TGraph*, TString>> slice;
414  slice.clear();
415 
416  struct th23helper{
417  std::vector<double> seeds;
418  const IFitVar * var;
419  TString label;
420  };
421 
422 
423  std::vector <th23helper> th23seeds;
424 
425  th23seeds.push_back({ {0.4,0.6}, &kFitSinSqTheta23, "2017" });
426  th23seeds.push_back( { {0.45}, kFitSinSqTheta23LowerOctant, "LO SA" });
427  th23seeds.push_back( { {0.55}, kFitSinSqTheta23UpperOctant, "UO SA"});
428 
429  int id, revid;
430 
431  if(MakeFile){
432 
433  TFile * outfile = new TFile(("files/"+outfilename+TagName+whatplot+"_Slicesonly.root").c_str(),"recreate");
434 
435  outfile->cd();
436  TVectorD v(1);
437  v[0] = minchi23;
438  v.Write("overall_minchi");
439 
440  for (double th23 : {0.558}){ //2017 BF, SA LO, SA UO : 0.558, 0.404, 0.623
441 
442  th23helper th23seed;
443 
444  if(th23 == 0.558) {th23seed = th23seeds[0];}
445  if(th23 == 0.404) {th23seed = th23seeds[1];}
446  if(th23 == 0.623) {th23seed = th23seeds[2];}
447 
448  for(int hie : {-1,1}){
449 
450  int steps = 25; double chi2r[steps], chi2w[steps], chi[steps]; double delta[steps]; double delt= 0;
451  double step = 2.0/(steps-1);
452  delt=-step;
453  for (int i =0; i<steps; i++){
454 
455  delt = delt+step;
456 
457  ResetOscCalcToDefault(calc);
458  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
459  calc->SetTh23(asin(sqrt(th23)));
460  calc->SetdCP(delt*M_PI);
461 
462  std::vector <const IExperiment * > expts_temp;
463  for(int i = 0; i < (int) predsvector.size(); ++i){
464 
465  double pot;
466  if(ana=="joint" && (i==1 || i>5)) pot=potFD_rhc; //sigh
467  else pot=potFD_rhc;
468  auto thisdata = GetFakeData (predsvector[i],calc, pot, cosmicsvector[i].first);
469  expts_temp.push_back(new SingleSampleExperiment(predsvector[i], *thisdata,
470  cosmicsvector[i].first, cosmicsvector[i].second));
471  }
472 
473  expts_temp.push_back(WorldReactorConstraint2017());
474 
475  auto exptsAll_temp = new MultiExperiment(expts_temp);
476 
477  if(isSyst){
478  systs = GetJointFitSystematicList(true, false, false, true, true, true);
479  auto notfornueFHC = GetCorrelations(true, true);
480  auto notfornueRHC = GetCorrelations(true, false);
481  auto notfornumu = GetCorrelations(false, true); //it will never enter the FHC, RHC ifs, so doesn't matter
482  if(ana=="fhc"){
483  exptsAll_temp->SetSystCorrelations(0, notfornueFHC);
484  for(int i=1; i<5;i++) {
485  exptsAll_temp->SetSystCorrelations(i, notfornumu);
486  }
487  }
488  if(ana=="joint"){
489  exptsAll_temp->SetSystCorrelations(0, notfornueFHC);
490  exptsAll_temp->SetSystCorrelations(1, notfornueRHC);
491  for(int i=2; i<10;i++) {exptsAll_temp->SetSystCorrelations(i, notfornumu);}
492  }
493 
494  auxShifts = SystShifts::Nominal();
495  //Registry::Print();
496  }
497  if(whatplot=="mh"){
498 
499  MinuitFitter fit(exptsAll_temp, {&kFitDeltaInPiUnits, &kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13}, systs);
500 
501  ResetOscCalcToDefault(calc);
502  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
503  chi2r[i] = fit.Fit(calc, auxShifts, {{ &kFitDeltaInPiUnits, {0.,0.5, 1., 1.5} }, {&kFitSinSqTheta23, {0.4,0.6}}},
504  {}, IFitter::kQuiet)->EvalMetricVal();
505 
506  ResetOscCalcToDefault(calc);
507  calc->SetDmsq32(-hie*fabs(calc->GetDmsq32()));
508  chi2w[i] = fit.Fit(calc, auxShifts, {{ &kFitDeltaInPiUnits, {0.,0.5, 1., 1.5} }, {&kFitSinSqTheta23, {0.4,0.6}}},
509  {}, IFitter::kQuiet)->EvalMetricVal();
510 
511  }
512 
513  if(whatplot=="dcp"){
514 
515  MinuitFitter fit(exptsAll_temp, {&kFitSinSqTheta23, &kFitDmSq32, &kFitSinSq2Theta13}, systs);
516 
517  ResetOscCalcToDefault(calc);
518  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
519  calc->SetdCP(delt*M_PI);
520  chi2r[i] = fit.Fit(calc, auxShifts,
521  {{&kFitSinSqTheta23, {0.4, 0.6}}, {&kFitDmSq32,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())}}},
522  {}, IFitter::kQuiet);
523 
524  ResetOscCalcToDefault(calc);
525  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
526  calc->SetdCP(0);
527  double chi2w0 = fit.Fit(calc, auxShifts,
528  {{&kFitSinSqTheta23, {0.4, 0.6}}, {&kFitDmSq32,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())}}}, {}, IFitter::kQuiet)->EvalMetricVal();
529 
530  ResetOscCalcToDefault(calc);
531  calc->SetdCP(M_PI);
532  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
533  double chi2wpi = fit.Fit(calc, auxShifts,
534  {{&kFitSinSqTheta23, {0.4, 0.6}}, {&kFitDmSq32,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())}}},
535  {}, IFitter::kQuiet)->EvalMetricVal();
536 
537  chi2w[i] = std::min(chi2w0, chi2wpi);
538 
539  }
540 
541  if(whatplot=="maxmix"){
542 
543  MinuitFitter fit(exptsAll_temp, {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13}, systs); //for 0.5
544  MinuitFitter fitth23(exptsAll_temp, {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13, &kFitSinSqTheta23}, systs); // free th23 for fake data
545 
546  ResetOscCalcToDefault(calc);
547  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
548  calc->SetTh23(asin(sqrt(th23)));
549  chi2r[i] = fit.Fit(calc, auxShifts,
550  {{ &kFitDeltaInPiUnits, {0.,0.5, 1., 1.5} }, {&kFitDmSq32,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())}},
551  {&kFitSinSqTheta23, {0.4,0.6}}},
552  {}, IFitter::kQuiet)->EvalMetricVal();
553 
554  ResetOscCalcToDefault(calc);
555  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
556  calc->SetTh23(asin(sqrt(0.5)));
557  chi2w[i] = fit.Fit(calc, auxShifts,
558  {{ &kFitDeltaInPiUnits, {0.,0.5, 1., 1.5}}, {&kFitDmSq32,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())}}},
559  {}, IFitter::kQuiet)->EvalMetricVal();
560 
561  }
562 
563  if(whatplot=="oct"){
564 
565  if(th23==0.404) {id=1; revid=2;} // primitive flip between octants: minimize LO with kFitSinSqTheta23LowerOctant, then use wrong octant and minimize there with kFitSinSqTheta23UpperOctant
566  if(th23==0.623 || th23==0.558) {id=2; revid=1;} // vice versa
567 
568  MinuitFitter fit(exptsAll_temp, {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13, th23seeds[id].var}, systs); // right oct
569  MinuitFitter fitwo(exptsAll_temp, {&kFitDeltaInPiUnits, &kFitDmSq32, &kFitSinSq2Theta13, th23seeds[revid].var}, systs); //wrong oct
570 
571  ResetOscCalcToDefault(calc);
572  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
573  calc->SetTh23(asin(sqrt(th23)));
574  chi2r[i] = fit.Fit(calc, auxShifts,
575  {{ &kFitDeltaInPiUnits, {0.,0.5, 1., 1.5} }, {&kFitDmSq32,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())}},
576  {th23seeds[id].var, th23seeds[id].seeds}},
577  {}, IFitter::kQuiet)->EvalMetricVal();
578 
579  ResetOscCalcToDefault(calc);
580  calc->SetDmsq32(hie*fabs(calc->GetDmsq32()));
581  calc->SetTh23(asin(sqrt(1-th23)));
582  chi2w[i] = fitwo.Fit(calc, auxShifts,
583  {{ &kFitDeltaInPiUnits, {0.,0.5, 1., 1.5}}, {&kFitDmSq32,{fabs(calc->GetDmsq32()), -fabs(calc->GetDmsq32())}},
584  {th23seeds[revid].var, th23seeds[revid].seeds}},
585  {}, IFitter::kQuiet)->EvalMetricVal();
586 
587  }
588 
589  delta[i] = delt;
590  chi[i] = sqrt(fabs(chi2w[i]-chi2r[i]));
591 
592  //cout<<delt<<" chi2 "<<chi[i]<<" right "<<chi2r[i]<<" wrong assumption :"<<chi2w[i]<<endl;
593  }
594  TString name = th23seed.label+" "+(hie>0?"NH":"IH");
595  slice.push_back({new TGraph(steps, delta, chi), name}); //0 is ih, 1 is nh
596 
597  } // end hie
598  } //end th23
599  outfile->cd();
600 
601  TVectorD ss(1);
602  ss[0] = slice.size();
603  ss.Write("numsl");
604 
605  for(int i=0; i<(int) slice.size(); i++){
606 
607  slice[i].first->Write(slice[i].second);
608 
609  }
610 
611  outfile->Close();
612  cout<<"Done with slices"<<endl;
613  }//end Make plots
614 
615  else{
616 
617  TFile * infile = new TFile (("files/"+outfilename+TagName+whatplot+"_Slicesonly.root").c_str(),"read");
618  auto mins =* (TVectorD*)infile->Get("overall_minchi");
619  double minchi23 = mins[0];
620  cout<<"Min chi: "<<minchi23<<endl;
621 
622  auto ss =* (TVectorD*)infile->Get("numsl");
623  double numsl = ss[0];
624  cout<<"Num slices: "<<numsl<<endl;
625 
626  std::vector<std::pair<TGraph*, TString>> slice;
627  std::vector<std::pair<TGraph*, TString>> slice2017;
628 
629  for (double th23 : {0.558}){ //2017 BF, SA LO, SA UO : 0.558, 0.404, 0.623
630 
631  th23helper th23seed;
632  if(th23 == 0.558) {th23seed = th23seeds[0];}
633  if(th23 == 0.404) {th23seed = th23seeds[1];}
634  if(th23 == 0.623) {th23seed = th23seeds[2];}
635 
636  for(int hie : {-1,1}){
637 
638  TString name = th23seed.label+" "+(hie>0?"NH":"IH");
639  TGraph* get_slice;
640  infile->GetObject(name, get_slice);
641  if(name.Contains("2017")) slice2017.push_back({get_slice,name});
642  else slice.push_back({get_slice,name});
643 
644  }
645 
646  }
647 
648  double max;
649  if(whatplot=="dcp") max=3.1;
650  else max = 5.1;
651  /* DrawSliceCanvas("delta", max, 0, 2);
652  for(int i =0; i < slice.size(); i++){
653 
654  if(slice[i].second.Contains("IH")) slice[i].first->SetLineColor(kPrimColorIH);
655  if(slice[i].second.Contains("NH")) slice[i].first->SetLineColor(k3SigmaConfidenceColorNH);
656  if(slice[i].second.Contains("UO")) slice[i].first->SetLineStyle(7);
657  // if(slice[i].second.Contains("2017")) slice[i].first->SetLineStyle(10);
658  slice[i].first->SetLineWidth(3);
659  slice[i].first->Draw("same");
660  }
661 
662  auto leg = SliceLegend(slice, true, 0.53);
663  leg->Draw();
664  gPad->SetGrid();
665  Simulation();
666 
667  gPad->Print((TString)"plots/"+version+"sens"+TagName+"_"+whatplot+"_SA.pdf");*/
668 
669  new TCanvas;
670  DrawSliceCanvas("delta", max, 0, 2);
671  for(int i =0; i < slice2017.size(); i++){
672 
673  if(slice2017[i].second.Contains("IH")) slice2017[i].first->SetLineColor(kPrimColorIH);
674  if(slice2017[i].second.Contains("NH")) slice2017[i].first->SetLineColor(k3SigmaConfidenceColorNH);
675 
676  slice2017[i].first->SetLineWidth(3);
677  slice2017[i].first->Draw("same");
678  }
679 
680  auto leg2 = SliceLegend(slice2017, true, 0.53);
681  leg2->Draw();
682  gPad->SetGrid();
683  Simulation();
684 
685  gPad->Print((TString)"plots/"+version+"sens"+TagName+"_"+whatplot+"_2017.pdf");
686 
687  }//end !Make plots = figures
688  }//end Makeslice
689 
690 }
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
::xsd::cxx::tree::id< char, ncname > id
Definition: Database.h:165
const Color_t kPrimColorIH
Definition: Style.h:64
void Simulation()
Definition: tools.h:16
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2.)
const XML_Char * name
Definition: expat.h:151
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
Antineutrinos-only.
Definition: IPrediction.h:50
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
(&#39; appearance&#39;)
Definition: IPrediction.h:18
double th23
Definition: runWimpSim.h:98
double delta
Definition: runWimpSim.h:98
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
(&#39;beam &#39;)
Definition: IPrediction.h:15
T sqrt(T number)
Definition: d0nt_math.hpp:156
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
Float_t ss
Definition: plot.C:24
const Color_t k2SigmaConfidenceColorNH
Definition: Style.h:77
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
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)
static SystShifts Nominal()
Definition: SystShifts.h:34
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
virtual void SetDmsq32(const T &dmsq32)=0
#define M_PI
Definition: SbMath.h:34
const double kAna2018RHCPOT
Definition: Exposures.h:208
const char * label
Log-likelihood scan across two parameters.
string outfilename
knobs that need extra care
Charged-current interactions.
Definition: IPrediction.h:39
string infile
const Color_t k3SigmaConfidenceColorNH
Definition: Style.h:78
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
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
virtual T GetDmsq32() const
Definition: IOscCalc.h:91
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta)
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2018(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
Neutrinos-only.
Definition: IPrediction.h:49
void DrawContourCanvas(TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1)
OStream cout
Definition: OStream.cxx:6
(&#39; survival&#39;)
Definition: IPrediction.h:19
const XML_Char * version
Definition: expat.h:187
static float min(const float a, const float b, const float c)
Definition: absgeo.cxx:45
TH2 * ToTH2(double minchi=-1) const
Definition: ISurface.cxx:161
TLegend * ContourLegend(int hie, bool fillContour, bool fccorr, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor, bool bestFit)
Combine multiple component experiments.
const Style_t k2SigmaConfidenceStyle
Definition: Style.h:71
std::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool GetFromUPS=false, bool NERSC=false)
double livetime
Definition: saveFDMCHists.C:21
void SaveTo(TDirectory *dir, const std::string &name) const
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
Neutral-current interactions.
Definition: IPrediction.h:40
Interface definition for fittable variables.
Definition: IFitVar.h:16
Both neutrinos and antineutrinos.
Definition: IPrediction.h:52
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
T min(const caf::Proxy< T > &a, T b)
All neutrinos, any flavor.
Definition: IPrediction.h:26
Float_t e
Definition: plot.C:35
void sensitivity2018(bool MakeFile=false, std::string ana="joint", bool isSyst=false, bool MakeSurf=false, std::string whatsurf="th23dcp", bool isNH=true, bool MakeSens=false, std::string whatplot="mh")
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
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)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
const double kAna2018FHCLivetime
Definition: Exposures.h:213
const Color_t k2SigmaConfidenceColorIH
Definition: Style.h:82
FILE * outfile
Definition: dump_event.C:13
const double kAna2018RHCLivetime
Definition: Exposures.h:214
virtual void SetdCP(const T &dCP)=0
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string