joint_fit_2018_tools.h
Go to the documentation of this file.
1 #pragma once
2 
3 // TODO: find a better home for this header file; functions are duplicated
4 
9 //#include "CAFAna/Prediction/PredictionSystJoint2018.h"
17 #include "CAFAna/Systs/XSecSysts.h"
18 #include "CAFAna/Vars/FitVars.h"
19 
20 #include "CAFAna/Analysis/Style.h"
22 
23 #include "Utilities/func/MathUtil.h"
24 
25 #include "TCanvas.h"
26 #include "TBox.h"
27 #include "TLatex.h"
28 #include "TColor.h"
29 #include "TGraph.h"
30 #include "TVectorD.h"
31 #include "TF1.h"
32 #include "TLegend.h"
33 #include "TLine.h"
34 #include "TMarker.h"
35 #include "TStyle.h"
36 #include "TSystem.h"
37 #include "TGaxis.h"
38 
39 using namespace ana;
40 
41 #include "OscLib/IOscCalc.h"
42 
43 #include "TFile.h"
44 
45 TFile* GetGaussianSurface(TString par, TString hier)
46 {
47  TString subdir = par.Contains("ssth23dmsq32")? "th23_dmsq" : "delta_th23";
48  TString gaussDir = "/nova/ana/nu_e_ana/Ana2018/Results/RHC_and_FHC/contours/"+subdir+"/syst/";
49  TString gaussName = "hist_contours_2018_joint_realData_both_only";
50  gaussName += hier;
51  if(par.Contains("ssth23dmsq32"))
52  gaussName += "combo_dmsq_systs.root";
53  else
54  gaussName += "combo_systs.root";
55 
56 
57  std::cout << "Opening " << gaussDir + gaussName << std::endl;
58  TFile* ret = TFile::Open(gaussDir + gaussName);
59  return ret;
60 }
61 
62 
64 {
65  TH2* temp = (TH2*)h->Clone();
66 
67  int NX = h->GetXaxis()->GetNbins();
68  int NY = h->GetYaxis()->GetNbins();
69  double minX = h->GetXaxis()->GetBinLowEdge(1);
70  double maxX = h->GetXaxis()->GetBinUpEdge(NX);
71  double widthX = h->GetXaxis()->GetBinWidth(1);
72  double minY = h->GetYaxis()->GetBinLowEdge(1);
73  double maxY = h->GetYaxis()->GetBinUpEdge(NY);
74  double widthY = h->GetYaxis()->GetBinWidth(1);
75 
76  // assert(NX % 2 == 0 &&
77  // "Histogram must have an even number of X bins for this to work");
78 
79  // Save the location we expect the end column of the original
80  // TH2 to be after we do the folding
81  int foldedBin = 1 + NX/2;
82 
83  // Fold the initial histogram to join the ends of the X axis
84  for(int i = 1; i <= NX/2; i++) {
85  for(int j = 1; j <= NY/2; j++) {
86  // x1, y1 denotes a position in the left-hand side half plane
87  // x2, y2 denotes a position in the right-hand side half plane
88  int x1 = i+NX/2;
89  int x2 = i%(NX/2+1)+2;
90  int y1 = j;
91  int y2 = j;
92  double z1 = h->GetBinContent(x1, y1);
93  double z2 = h->GetBinContent(x2, y2);
94  temp->SetBinContent(x1, y1, z2);
95  temp->SetBinContent(x2, y2, z1);
96  }
97  }
98  // smooth both histograms
99  temp->Smooth();
100  h->Smooth();
101 
102  // Replace end columns
103  for(int y = 1; y <= NY; y++) {
104  double sZ1 = temp->GetBinContent(foldedBin, y);
105  double sZ2 = temp->GetBinContent(foldedBin, y);
106  h->SetBinContent(1, y, sZ1);
107  h->SetBinContent(NX, y, sZ2);
108  }
109 }
110 
111 std::string GetNuePredPath(std::string beam, bool isFakeData, bool NERSC) {
113  bool isFHC = true;
114  if(beam == "rhc") isFHC = false;
115 
116  // default path name for Nue Preds constructor
117  // to point to official files
118  if(!NERSC) return "";
119 
120  // NERSC paths
121  else {
122  std::string dir = std::string(getenv("FCHELPERANA2018_LIB_PATH"))
123  + "/Predictions/provisional_singles";
124  if(isFHC && !isFakeData)
125  ret = dir + "/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_FHCAllSyst_nueconcat_realNDData_2018-05-22.root";
126  else if(!isFHC && !isFakeData)
127  ret = dir + "/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_RHCAllSyst_nueconcat_realNDData_2018-05-22.root";
128  else if(isFHC && isFakeData)
129  ret = dir + "/NueProd4SystematicsOnFakeNDData_2018-05-03/pred_allxp_Nue2018Axis_full_FHCAllSysts_nueconcat_fakeNDData_05_03_2018_merged.root";
130  else if(!isFHC & isFakeData)
131  ret = dir + "/NueProd4SystematicsOnFakeNDData_2018-05-03/pred_allxp_Nue2018Axis_full_RHCAllSysts_nueconcat_fakeNDData_05_03_2018_merged.root";
132  }
133  return ret;
134 }
135 
138  bool corrSysts,
139  std::string beam,
140  bool isFakeData,
141  bool GetFromUPS=false,
142  bool minimizeMemory = false,
143  bool NERSC = false){
144 
146  if(decomp == "noextrap") extrap = kNoExtrap;
147  if(decomp == "prop") extrap = kProportional;
148  if(decomp == "combo") extrap = kCombo;
149  if(decomp == "fakexp") extrap = kFake;//Hack
150 
151  PredictionSystJoint2018 * pred_fid;
152  std::string nue_pred_path = GetNuePredPath(beam, isFakeData, NERSC);
153  // Interactive running
154  if (!GetFromUPS){
155  if(corrSysts){
156  std::cout << "Getting Nue Prediction from " + nue_pred_path << endl;
157  if(beam=="rhc") pred_fid = new PredictionSystJoint2018 (extrap, calc,
159  beam, isFakeData, true,
160  nue_pred_path, minimizeMemory);
161  else if(beam=="fhc") pred_fid = new PredictionSystJoint2018 (extrap, calc,
163  beam, isFakeData, true,
164  nue_pred_path, minimizeMemory);
165  else{ std::cout << "Beam mode needs to be either fhc or rhc" << endl; exit(1);}
166  }
167  else pred_fid = new PredictionSystJoint2018 (extrap, calc,{}, beam, isFakeData,
168  true, nue_pred_path, minimizeMemory);
169  }
170  // Grid running for FC
171  else{
172  std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH")); // need an update
173  std::string name = "/pred_nue_reduced_v9.root"; //need an update
174  if(corrSysts){
175  if(beam=="rhc") pred_fid = new PredictionSystJoint2018(extrap, calc,
177  beam, isFakeData, true,
178  upsname+name, minimizeMemory);
179  else if(beam=="fhc") pred_fid = new PredictionSystJoint2018(extrap, calc,
181  beam, isFakeData, true,
182  upsname+name, minimizeMemory);
183  else{ std::cout << "Beam mode should be either fhc or rhc" << endl; exit(1);}
184  }
185  else
186  pred_fid = new PredictionSystJoint2018(extrap, calc, {},beam, isFakeData,
187  true, upsname+name, minimizeMemory);
188  }
189 
190 
191 
192  double pot;
193  if(beam=="rhc") pot = kAna2018RHCPOT;
194  else pot = kAna2018FHCPOT;
195  std::cout << "Scale to "<<pot<<std::endl;
196 
197  std::cout << "Nue prediction (no rock)"
198  << pred_fid->Predict(calc).Integral(pot) <<std::endl;
199  // Load rock
200 
201  double nonsScale = 1./10.45;
202  double swapScale = 1./12.91;
203  //fix to the normal path later
204  std::string dirName = "/nova/ana/nu_e_ana/Ana2018/Predictions/FDRock/";
205  if(NERSC) dirName = std::string(getenv("FCHELPERANA2018_LIB_PATH"))+"/Predictions/FDRock/";
206  std::string fName2 = "fdrock_nue2018.root";
207  std::string rName;
208 
209  if(beam=="rhc") rName = "pred_FDRock_merg_rhc";
210  else rName = "pred_FDRock_merg_fhc";
211  if( decomp == "noextrap"){
212  if(beam=="rhc") rName = "pred_FDRock_rhc";
213  else rName = "pred_FDRock_fhc";
214  }
215  // Point to the rock file in the ups product for FC'ing
216  if (GetFromUPS){
217  dirName = std::string(getenv("FCHELPERANA2017_LIB_PATH")); //need an update
218  fName2 = "/fdrock_nue2017.root"; //need an update
219  rName = "pred_FDRock";
220  }
221 
222  auto rock = LoadFromFile<PredictionNoExtrap>(dirName + fName2, rName).release();
223 
224  std::cout << "Nue prediction (rock only - not scaled) "
225  << rock->Predict(calc).Integral(pot) <<std::endl;
226  std::cout << "Adding Rock" << endl;
227  auto nue_pred = new PredictionAddRock (pred_fid, rock, nonsScale, swapScale);
228 
229  std::cout << "Nue prediction (rock added) "
230  << nue_pred->Predict(calc).Integral(pot) << std::endl;
231  std::cout << "Pure rock events: "
232  << nue_pred->Predict(calc).Integral(pot) - pred_fid->Predict(calc).Integral(pot) << std::endl;
233 
234 
235  return nue_pred;
236 }
237 
238 //////////////////////////////////////////////////////////////////////
239 
240 std::pair <TH1D*, double > GetNueCosmics2018 (std::string beam, bool GetFromUPS=false, bool NERSC=false){
241  Spectrum* outtime;
242  if (GetFromUPS){
243  std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH")); //update
244  std::string fname = "/nue2017_cospred.root"; //update
245  std::string name = "cos";
246  outtime = LoadFromFile<Spectrum>(upsname+fname,name).release();
247  }
248  else{
250  std::string anaDir = "/nova/ana/nu_e_ana/Ana2018/";
251  if(NERSC) anaDir = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
252  if(beam=="rhc") name = "cosmic_spect_rhc";
253  else name="cosmic_spect_fhc";
254  outtime = LoadFromFile<Spectrum>
255  (anaDir+"/Predictions/cosmic/cosmic_prediction_real_data.root", name ).release();
256  }
257 
258  double livetime;
259  if(beam=="rhc") livetime = kAna2018RHCLivetime ;
260  else livetime = kAna2018FHCLivetime ;
261 
262  //double outN = outtime->Integral(outtime->Livetime(), 0, kLivetime);
263  double outN = outtime->Integral(livetime, 0, kLivetime);
264 
265  std::cout << "\nAdding " << outtime->Integral(livetime, 0, kLivetime)
266  << " cosmics from out-of-time." << std::endl;
267  std::cout << "Which scale to " << outtime->Integral(livetime, 0,
268  kLivetime)
269  << " in the spill window." << std::endl;
270 
271  std::cout << "uncertainty " << 1/sqrt(outN) << std::endl;
272 
273  TH1 *h = outtime->ToTH1(livetime, kBlack, kSolid, kLivetime);
274 // assert( h->Integral() > 4 && "Nue cosmic prediction < 4. Something is wrong"); //RHC has ~0.6 cosmic events
275 
276  return {outtime->ToTH1(livetime, kBlack, kSolid, kLivetime), 1/sqrt(outN)};
277 }
278 
279 //////////////////////////////////////////////////////////////////////
280 
282  Spectrum* intime;
283 
284  if(beam=="fhc"){
285  intime = LoadFromFile<Spectrum>
286  ("/nova/ana/nu_e_ana/Ana2018/Data2018/spectra.root", "data_spect_fhc").release();
287  }
288  else if(beam=="rhc"){
289  intime = LoadFromFile<Spectrum>
290  ("/nova/ana/nu_e_ana/Ana2018/Data2018/spectra.root", "data_spect_rhc").release();
291  }
292  else{ std::cout << "Beam mode != fhc or rhc" << std::endl; exit(1);}
293  double POT, livetime;
294  if(beam == "rhc") {livetime = kAna2018RHCLivetime; POT = kAna2018RHCPOT;}
295  if(beam == "fhc") {livetime = kAna2018FHCLivetime; POT = kAna2018FHCPOT;}
296 
297  std::cout << "Loaded data spectrum:\n"
298  << "POT " << intime->POT()
299  << " (expected " << POT << ")\n"
300  << "Livetime " << intime->Livetime()
301  << " (expected " << livetime << ")\n"
302  << "Integral " << intime->Integral(intime->POT())
303  << std::endl;
304 
305  return intime;
306 }
307 
308 //////////////////////////////////////////////////////////////////////
309 
310 std::vector <const IPrediction*> GetNumuPredictions2018(const int nq = 4,
311  bool useSysts = true,
312  std::string beam="fhc",
313  bool GetFromUPS=false,
314  ENu2018ExtrapType numuExtrap = kNuMu,
315  bool minimizeMemory = false,
316  bool NERSC = false)
317 {
318  auto calc = DefaultOscCalc();
319 
320  std::vector <const IPrediction*> numu_preds;
321 
323  std::string anaDir = "/nova/ana/nu_mu_ana/Ana2018/";
324  if(NERSC) anaDir = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
325  if(beam=="rhc") filename=anaDir + "/Predictions/pred_numuconcat_rhc__numu2018.root";
326  else filename= anaDir + "/Predictions/pred_numuconcat_fhc__numu2018.root";
327 
328  std::cout << "\nLoading Numu Predictions\n" ;
329  for (int quant = 0; quant < nq; ++quant){
330  if (GetFromUPS){
331  std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH")); //update
332  std::string fname = "/pred_numu_reduced_v5.root"; //update
333  if(!useSysts) {
334  numu_preds.push_back(new PredictionSystJoint2018(numuExtrap,
335  calc, beam,
336  quant+1,{},
337  upsname+fname,
338  minimizeMemory)); //update
339  }
340  else{
341  numu_preds.push_back(new PredictionSystJoint2018(numuExtrap, calc, beam,
342  quant+1,
344  upsname+fname,
345  minimizeMemory)); //update
346  }
347  }
348  else{
349  if(!useSysts) {
350  numu_preds.push_back(new PredictionSystJoint2018(numuExtrap, calc, beam, quant+1,
351  {}, filename, minimizeMemory));
352  }
353  else {
354  numu_preds.push_back(new PredictionSystJoint2018(numuExtrap, calc, beam, quant+1,
356  filename, minimizeMemory));
357  }
358  }
359  }
360  return numu_preds;
361 }
362 //////////////////////////////////////////////////////////////////////
363 
364 std::vector <std::pair <TH1D*, double> > GetNumuCosmics2018(const int nq = 4, std::string beam="fhc",
365  bool GetFromUPS=false, bool NERSC = false)
366 {
367  cout<<"Beam is "<<beam<<endl;
368  //cosmics
369  //the versionn in cafana has float not double and it's a pain
370  TFile *fcosm;
371  if (GetFromUPS){
372  std::string upsname = std::string(getenv("FCHELPERANA2017_LIB_PATH")); //update
373  fcosm = TFile::Open((upsname+"/cos_bkg_hists_v2.root").c_str()); //update
374  }
375  else{
376  std::string dir = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
377  if(NERSC) dir = std::string(getenv("FCHELPERANA2018_LIB_PATH"));
379  if (beam=="rhc") filename="cosmics_rhc__numu2018.root";
380  else filename="cosmics_fhc__numu2018.root";
381  fcosm = TFile::Open((dir+filename).c_str());
382  }
383 
384  if(fcosm->IsZombie()) {std::cerr<< "bad cosmics\n"; exit(1);}
385  std::vector <std::pair <TH1D*, double > > numu_cosmics;
386 
387  for (int i = 1; i <=nq; ++i) {
388  auto h = (TH1D*) fcosm->Get((TString)"cosmics_q"+std::to_string(i))
389  ->Clone(UniqueName().c_str());
390  //auto hs = (TH1D*) fcosm->Get((TString)"q"+std::to_string(i)+"allsyst"); //update. no such stuff in numu files
391  //numu_cosmics.push_back({h,hs->GetBinContent(1)});
392  numu_cosmics.push_back({h,1/sqrt(h->Integral())});
393  std::cout << "Cosmics all periods quantile " << i << " "
394  << h->Integral()
395  << " scale error "
396  //<< hs->GetBinContent(1)
397  << 1/sqrt(h->Integral())
398  << "\n";
399  }//end quantiles
400  return numu_cosmics;
401 }
402 //////////////////////////////////////////////////////////////////////
403 
404 std::vector <Spectrum * > GetNumuData2018(const int nq = 4, std::string beam = "fhc"){
405 
407  if(beam == "fhc") {file="/nova/ana/nu_mu_ana/Ana2018/Data/fd_dataspectrum_fhc_full__numu2018.root";}
408  if(beam == "rhc") {file="/nova/ana/nu_mu_ana/Ana2018/Data/fd_dataspectrum_rhc_full__numu2018.root";}
409 
410  std::vector <Spectrum *> numu_data;
411 
412  for (int i = 1; i <=nq; ++i) {
413 
414  auto f= new TFile(file.c_str());
415  if(f->IsZombie()) {std::cerr<< "bad data\n"; exit(1);}
416  numu_data.push_back(Spectrum::LoadFrom(f, ("fd_data_q"+std::to_string(i)).c_str())).release();
417 
418  }//end quantiles
419  return numu_data;
420 }
421 
422 //////////////////////////////////////////////////////////////////////
423 
425  const double pot, TH1D* cosmics = 0)
426 {
427  auto temp = pred->Predict(calc).FakeData(pot);
428  if(cosmics) temp += Spectrum(cosmics, pot,0);
429  return new Spectrum(temp);
430 }
431 
433  const double pot, TH1D* cosmics = 0)
434 {
435  auto temp = pred->Predict(calc).MockData(pot);
436  if(cosmics) temp += Spectrum(cosmics, pot,0);
437  return new Spectrum(temp);
438 }
439 
440 //////////////////////////////////////////////////////////////////////
441 
443  double dmsq32 = -5, double dCP_Pi = -5)
444 {
445  if(dCP_Pi > -1) kFitDeltaInPiUnits.SetValue(calc, dCP_Pi);
446  if(ssth23 > -1) kFitSinSqTheta23.SetValue(calc, ssth23);
447  if(dmsq32 > -1) kFitDmSq32.SetValue(calc, dmsq32);
448 }
449 
450 //////////////////////////////////////////////////////////////////////
451 
452 std::vector<const ISyst*> GetJointFitSystematicList(bool corrSysts, bool nueExclusive = false, bool numuExclusive=false, bool isFHC = true, bool isRHC = true , bool intersection = true){
453  std::vector<const ISyst*> ret ={} ;
454  if(corrSysts){
455  if(! intersection){
456  if(nueExclusive){
457 
458  ret.push_back(&kRadCorrNue);
459  ret.push_back(&kRadCorrNuebar);
460  ret.push_back(&k2ndClassCurrs);
463 
464  }
465  if(numuExclusive){
466 
470  if(isFHC) ret.push_back(&kAna2018NormRHC);
471  if(isRHC) ret.push_back(&kAna2018NormFHC);
473 
474  }
475  }
476  else{
477  if(nueExclusive) ret = getAllAna2018Systs(kNueAna2018);
478  if(numuExclusive) ret = getAllAna2018Systs(kNumuAna2018);
479  if(!nueExclusive && !numuExclusive) ret = getAllAna2018Systs(kJointAna2018, true, kBoth, true);
480  //std::cout <<"I'm in else" <<std::endl;
481  //if(!nueExclusive && !numuExclusive) ret = getAllAna2018Systs(kJointAna2018, true, true);
482  }
483  }
484  std::cout << "Return list of systematics: " ;
485  for(auto & s:ret) std::cout << s->ShortName() << ", ";
486  std::cout << std::endl;
487  return ret;
488 }
489 std::vector<std::pair<const ISyst*,const ISyst*> > GetCorrelations (bool isNue, bool isFHC){
490  std::vector<std::pair<const ISyst*,const ISyst*> > ret ={};
491  std::cout<<(isNue?"not nue ":"not numu ")<<(isFHC&&!isNue?" ":(isFHC?" FHC ":" RHC "));
492  //to first approximation if they are not the same they are not correlated.
493  //might have to be more complicated later
494  auto badsyst = GetJointFitSystematicList(true, !isNue, isNue, isFHC, !isFHC, false);
495  for (auto & s:badsyst) ret.push_back ({s,NULL});
496 
497  return ret;
498 }
499 
500 // For slices plots, plot by true octant, and make separate plots of each
502 {
503 public:
504  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
505  return util::sqr(sin(osc->GetTh23()));};
506  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
507  osc->SetTh23(asin(sqrt(Clamp(val))));};
508  virtual std::string ShortName() const {return "ssth23LO";}
509  virtual std::string LatexName() const {return "sin^{2}#theta_{23} LO";}
510  virtual double LowLimit() const {return 0;}
511  virtual double HighLimit() const {return 0.5;}
512 };
513 
515 
517 {
518 public:
519  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
520  return util::sqr(sin(osc->GetTh23()));};
521  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
522  osc->SetTh23(asin(sqrt(Clamp(val))));};
523  virtual std::string ShortName() const {return "ssth23UO";}
524  virtual std::string LatexName() const {return "sin^{2}#theta_{23} UO";}
525  virtual double LowLimit() const {return 0.5;}
526  virtual double HighLimit() const {return 1;}
527 };
529 
530 
531 ////////////////////////////////////////////////////////////////////////////
532 //////////////////////// Slice, Contours Style /////////////////////////////
533 ///////////////////////////////////////////////////////////////////////////
534 //move functions from plotting script here
535 
536 void DrawSliceCanvas(TString slicename, const double ymax,
537  const double xmin = 0, const double xmax = 2., bool overlay = false){
538  auto c1 = new TCanvas(ana::UniqueName().c_str());
539  c1->SetFillStyle(0);
540  c1->SetLeftMargin(0.14);
541  c1->SetBottomMargin(0.15);
542  TString title;
543  if(slicename.Contains("delta")) {
544  title = ";delta;Significance (#sqrt{#Delta#chi^{2}})";
545  }
546  if(slicename.Contains("th23")) {
547  title = ";sin^{2}#theta_{23};Significance (#sqrt{#Delta#chi^{2}})";
548  }
549  if(slicename.Contains("dmsq")) {
550  title = ";|#Deltam^{2}_{32}| (10^{-3} eV^{2});Significance (#sqrt{#Delta#chi^{2}})";
551  }
552  if(slicename.Contains("th13")) {
553  title = ";sin^{2}2#theta_{13};Significance (#sqrt{#Delta#chi^{2}})";
554  }
555  auto axes = new TH2F("",title, 100, xmin, xmax, 100, 0, ymax);
556  if(slicename.Contains("fccorr") || overlay) axes->GetYaxis()->SetTitle("Significance (#sigma)");
557  // CenterTitles(axes);
558  axes->GetXaxis()->CenterTitle();
559  axes->GetYaxis()->CenterTitle();
560  axes->Draw();
561  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
562  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
563  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
564  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
565  axes->GetXaxis()->SetTitleOffset(0.8);
566  axes->GetYaxis()->SetTitleOffset(0.75);
567  if(slicename.Contains("delta")) XAxisDeltaCPLabels(axes);
568  c1->RedrawAxis();
569  }
570 
571 TLegend * SliceLegend(std::vector <std::pair <TGraph*, TString > > & graphs, bool isDelta, double low=0.7, bool isDmsq = false, bool overlay = false)
572 {
573  TLegend * leg;
574  if(isDelta) leg = new TLegend(0.6, 0.6, 0.9,0.89);
575  else {
576  leg = new TLegend(0.18, 0.2, 0.4, 0.49);
577  }
578  leg->SetTextSize(0.8*kBlessedLabelFontSize);
579  leg->SetFillStyle(0);
580  for (auto &gr:graphs){
581  TString str="";
582  if(gr.second.Contains("O")){
583  if(isDelta) leg->SetX1(0.57);
584  if(isDmsq) leg->SetX1(0.18);
585  if(!isDelta && !isDmsq) {leg->SetX1(0.6); leg->SetX2(0.89);}
586  leg->SetMargin(0.18);
587  str += gr.second.Contains("NH") ? "NH":"IH";
588  str += gr.second.Contains("LO") ? " Lower": " Upper";
589  str += " octant";
590  }
591  else{
592  str="#splitline{";
593  str += gr.second.Contains("NH") ? "Normal":"Inverted";
594  str += "}{hierarchy}";
595  if(overlay) {
596  str = "Uncorr. ";
597  str += gr.second.Contains("NH") ? "NH" : "IH";
598  }
599  }
600  leg->AddEntry(gr.first, str, "l");
601  }
602  return leg;
603 }
604 
605 
606 void DrawContourCanvas (TString surfName, double minx = 0,
607  double maxx = 2, double miny = 0, double maxy = 1){
608  auto c1 = new TCanvas(ana::UniqueName().c_str());
609  c1->SetFillStyle(0);
610  c1->SetLeftMargin(0.14);
611  c1->SetBottomMargin(0.15);
612  TH2* axes = new TH2F();
613  TString title;
614  if(surfName.Contains("delta_th23")) title=";#delta_{CP};sin^{2}#theta_{23}";
615  if(surfName.Contains("th13_delta")) title=";sin^{2}2#theta_{13};#delta_{CP}/#pi";
616  if(surfName.Contains("th23_dm32")) title=";sin^{2}#theta_{23};#Deltam^{2}_{32} (10^{-3}eV^{2})";
617  axes = new TH2F("",title,100,minx,maxx,100,miny,maxy);
618 // CenterTitles(axes);
619  axes->GetXaxis()->CenterTitle();
620  axes->GetYaxis()->CenterTitle();
621  axes->Draw();
622  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
623  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
624  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
625  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
626  axes->GetXaxis()->SetTitleOffset(0.8);
627  axes->GetYaxis()->SetTitleOffset(0.8);
628  TGaxis::SetMaxDigits(3);
629  if(surfName.Contains("th23_dm32")) {
630  axes->GetYaxis()->SetDecimals() ;
631  axes->GetYaxis()->SetTitleOffset(0.85);
632  axes->GetYaxis()->SetLabelOffset(0.002);
633  axes->GetYaxis()->SetLabelSize(0.058);
634  axes->GetYaxis()->SetTitleSize(0.078);
635  axes->GetXaxis()->SetTitleSize(.95*kBlessedTitleFontSize);
636  axes->GetXaxis()->SetTitleOffset(0.86);
637  }
638  if(surfName.Contains("delta_th23")) XAxisDeltaCPLabels(axes);
639  c1->RedrawAxis();
640 }
641 
642 
643 TLegend * ContourLegend(int hie, bool fillContour, bool fccorr,
644  Int_t kFillColor1,Int_t kFillColor2, Int_t kFillColor3,Int_t kDarkColor,
645  bool bestFit){
646  TLegend * leg = new TLegend(0.16,fccorr?0.19:0.17,0.65,0.28);
647  leg->SetTextSize(kBlessedLabelFontSize);
648  leg->SetFillStyle(0);
649  leg->SetMargin(1.3*leg->GetMargin());
650  leg->SetNColumns(3);
651 
652  TGraph * dummy = new TGraph;
653  dummy->SetLineWidth(3);
654 
655  if(fillContour){
656  dummy->SetLineColor(kDarkColor);
657  dummy->SetFillColor(kFillColor1);
658  leg->AddEntry(dummy->Clone(),"1#sigma","f");
659  dummy->SetFillColor(kFillColor2);
660  leg->AddEntry(dummy->Clone(),"2#sigma","f");
661  dummy->SetFillColor(kFillColor3);
662  leg->AddEntry(dummy->Clone(),"3#sigma","f");
663  }
664  else{
665  leg->SetMargin(1.4*leg->GetMargin());
666  dummy->SetLineColor(kFillColor1);
667  dummy->SetLineStyle(kSolid);
668  leg->AddEntry(dummy->Clone(),"1#sigma","l");
669  dummy->SetLineStyle(k2SigmaConfidenceStyle);
670  dummy->SetLineColor(kFillColor2);
671  leg->AddEntry(dummy->Clone(),"2#sigma","l");
672  dummy->SetLineStyle(kSolid);
673  dummy->SetLineColor(kFillColor3);
674  leg->AddEntry(dummy->Clone(),"3#sigma","l");
675  }
676  if(bestFit){
677  leg->SetNColumns(4);
678  //leg->SetTextSize(0.8*kBlessedLabelFontSize);
679 // leg->SetColumnSeparation(0.08);
680  leg->SetX2(0.75);
681  dummy->SetMarkerStyle(43);
682  dummy->SetMarkerSize(2);
683  dummy->SetMarkerColor(kDarkColor);
684  dummy->SetLineColor(kDarkColor);
685  dummy->SetLineStyle(7);
686  if(hie>0) leg->AddEntry(dummy->Clone(),"Best Fit","p");
687  if(hie<0) leg->AddEntry((TObject*)0, "#color[0]{Best Fit}", "");
688  dummy->SetLineStyle(1);
689  }
690  if(!fccorr) leg->SetHeader("No Feldman-Cousins");
691  else leg->SetHeader("Feldman-Cousins");
692  return leg;
693 }
694 
695 
696 TLegend * FCOverlayContourLegend(int hie, bool fillContour,
697  Int_t kFillColor1,Int_t kFillColor2, Int_t kFillColor3,Int_t kDarkColor)
698 {
699  TLegend * leg = new TLegend(0.16,hie>0?.78:0.77,0.65,0.86);
700  leg->SetTextSize(kBlessedLabelFontSize);
701  leg->SetFillStyle(0);
702  leg->SetMargin(1.3*leg->GetMargin());
703  leg->SetNColumns(3);
704 
705  TGraph * dummy = new TGraph;
706  dummy->SetLineWidth(3);
707 
708  if(fillContour){
709  dummy->SetLineColor(kDarkColor);
710  dummy->SetFillColor(kFillColor1);
711  leg->AddEntry(dummy->Clone(),"1#sigma","f");
712  dummy->SetFillColor(kFillColor2);
713  leg->AddEntry(dummy->Clone(),"2#sigma","f");
714  dummy->SetFillColor(kFillColor3);
715  leg->AddEntry(dummy->Clone(),"3#sigma","f");
716  }
717  else{
718  leg->SetMargin(1.4*leg->GetMargin());
719  dummy->SetLineColor(kFillColor1);
720  dummy->SetLineStyle(kSolid);
721  leg->AddEntry(dummy->Clone(),"1#sigma","l");
722  dummy->SetLineStyle(k2SigmaConfidenceStyle);
723  dummy->SetLineColor(kFillColor2);
724  leg->AddEntry(dummy->Clone(),"2#sigma","l");
725  dummy->SetLineStyle(kSolid);
726  dummy->SetLineColor(kFillColor3);
727  leg->AddEntry(dummy->Clone(),"3#sigma","l");
728  }
729 
730  leg->SetHeader("No Feldman-Cousins");
731  return leg;
732 }
733 
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
const XML_Char * name
Definition: expat.h:151
TCut intime("tave>=217.0 && tave<=227.8")
double minY
Definition: plot_hist.C:9
const NOvARwgtSyst k2ndClassCurrs("2ndclasscurr","Second class currents", novarwgt::kSimpleSecondClassCurrentsSystKnob)
Second-class current syst. See documentation in NOvARwgt.
Definition: XSecSysts.h:71
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
std::map< std::string, double > xmax
double ssth23
Float_t y1[n_points_granero]
Definition: compare.C:5
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Definition: Spectrum.cxx:148
subdir
Definition: cvnie.py:7
Loads shifted spectra from files.
const FitDmSq32 kFitDmSq32
Definition: FitVars.cxx:18
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
double maxy
Float_t x1[n_points_granero]
Definition: compare.C:5
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
std::string GetNuePredPath(std::string beam, bool isFakeData, bool NERSC)
void XAxisDeltaCPLabels(TH1 *axes, bool t2kunits)
Label the x-axis with fractions of pi.
T sqrt(T number)
Definition: d0nt_math.hpp:156
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
Definition: FitVars.cxx:16
General interface to oscillation calculators.
Definition: StanTypedefs.h:22
Int_t par
Definition: SimpleIterate.C:24
virtual double HighLimit() const
std::vector< Spectrum * > GetNumuData2018(const int nq=4, std::string beam="fhc")
OStream cerr
Definition: OStream.cxx:7
virtual T GetTh23() const
Definition: IOscCalc.h:94
double maxY
Definition: plot_hist.C:10
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
string filename
Definition: shutoffs.py:106
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
void AddJointAna2018NotCorrelSysts(std::vector< const ISyst * > &systs, const EAnaType2018 ana, const BeamType2018 beam)
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)
void CylindricalSmoothing(TH2 *h)
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:249
virtual double HighLimit() const
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta, double low=0.7, bool isDmsq=false, bool overlay=false)
virtual double LowLimit() const
osc::OscCalcDumb calc
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:300
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
const double kAna2018RHCPOT
Definition: Exposures.h:208
const NOvARwgtSyst kRadCorrNuebar("radcorrnuebar","Radiative corrections for #bar{#nu}_{e}", novarwgt::kSimpleRadiativeCorrNuebarXsecSystKnob)
Radiative corrections syst (nuebars). See documentation in NOvARwgt.
Definition: XSecSysts.h:67
Spectrum * GetMockData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
Double_t ymax
Definition: plot.C:25
const XML_Char * s
Definition: expat.h:262
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Definition: Spectrum.cxx:535
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
Definition: Spectrum.cxx:349
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2., bool overlay=false)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::string getenv(std::string const &name)
#define pot
TLegend * ContourLegend(int hie, bool fillContour, bool fccorr, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor, bool bestFit)
const double j
Definition: BetheBloch.cxx:29
T Clamp(T val) const
Definition: IFitVar.h:60
std::vector< float > Spectrum
Definition: Constants.h:610
TLegend * FCOverlayContourLegend(int hie, bool fillContour, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor)
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
double POT() const
Definition: Spectrum.h:227
std::vector< std::pair< TH1D *, double > > GetNumuCosmics2018(const int nq=4, std::string beam="fhc", bool GetFromUPS=false, bool NERSC=false)
static bool isFHC
Oscillation probability calculators.
Definition: Calcs.h:5
OStream cout
Definition: OStream.cxx:6
std::vector< double > POT
virtual std::string LatexName() const
Spectrum Predict(osc::IOscCalc *calc) const override
const Style_t k2SigmaConfidenceStyle
Definition: Style.h:71
T sin(T number)
Definition: d0nt_math.hpp:132
void DrawContourCanvas(TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1)
const DummyAnaSyst kAna2018NormRHC("NormRHC2018","RHC. Norm.")
std::pair< TH1D *, double > GetNueCosmics2018(std::string beam, bool GetFromUPS=false, bool NERSC=false)
std::string dirName
Definition: PlotSpectra.h:47
virtual double LowLimit() const
double dmsq32
TDirectory * dir
Definition: macro.C:5
double livetime
Definition: saveFDMCHists.C:21
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
Definition: FitVars.cxx:17
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
void AddNonLoadable2018Systs(std::vector< const ISyst * > &systs, const EAnaType2018 ana)
exit(0)
TCut outtime("(tave>0.0&&tave<200.0) || (tave>250.0&&tave<450.0)")
TFile * file
Definition: cellShifts.C:17
c1
Definition: demo5.py:24
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 Float_t kBlessedLabelFontSize
Definition: Style.h:90
const IPrediction * GetNuePrediction2018(std::string decomp, osc::IOscCalc *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
const double kAna2018FHCPOT
Definition: Exposures.h:207
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
virtual std::string ShortName() const
const Float_t kBlessedTitleFontSize
Definition: Style.h:89
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)
const NOvARwgtSyst kRadCorrNue("radcorrnue","Radiative corrections for #nu_{e}", novarwgt::kSimpleRadiativeCorrNueXsecSystKnob)
Radiative corrections syst (nues). See documentation in NOvARwgt.
Definition: XSecSysts.h:64
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
Definition: Spectrum.h:230
Spectrum * GetNueData2018(std::string beam)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
void rock(std::string suffix="full")
Definition: rock.C:28
const DummyAnaSyst kAna2018NormFHC("NormFHC2018","FHC. Norm.")
const double kAna2018FHCLivetime
Definition: Exposures.h:213
virtual std::string LatexName() const
TFile * GetGaussianSurface(TString par, TString hier)
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
const double kAna2018RHCLivetime
Definition: Exposures.h:214
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
T asin(T number)
Definition: d0nt_math.hpp:60
virtual std::string ShortName() const
enum BeamMode string