make_starplots.C
Go to the documentation of this file.
1 #pragma once
2 
4 #include "../numu_tools.h"
5 #include "../settings.h"
6 #include "includes_starplots.h"
7 
8 #include "CAFAna/Core/Progress.h"
9 
10 using namespace ana;
11 
12 #include "TCanvas.h"
13 #include "TFile.h"
14 #include "TGraph.h"
15 #include "TGraphAsymmErrors.h"
16 #include "TH1.h"
17 #include "TH2.h"
18 #include "TLatex.h"
19 #include "TLegend.h"
20 #include "TStyle.h"
21 
22 std::string dataName = "fakedata";
23 
24 std::vector<Style_t> styleline = {kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid};
25 std::vector<Color_t> color = {kMagenta, kGreen+1, kOrange+1, kBlue-7, kCyan+1, kRed-4, kBlue+1, kGray+1, kOrange+3};
26 
27 
29  std::string anaName,
30  bool nh)
31 {
32  double sinHere = sin_rhcfhc_stats;
33  double deltaHere = delta_rhcfhc_stats;
34  // naively changing hierarchy
35  // this isnt totally right
36  double deltaSign = 1.;
37  if(nh==false) deltaSign = -1.;
38 
39  if(anaName=="rhc"){
40  sinHere = sin_rhc_stats;
41  deltaHere = delta_rhc_stats*deltaSign;
42  }
43  if(anaName=="fhc"){
44  sinHere = sin_fhc_stats;
45  deltaHere = delta_fhc_stats*deltaSign;
46  }
48  calc->SetTh23(asin(sqrt(sinHere)));
49  calc->SetDmsq32(deltaHere);
50 }
51 
52 
53 void make_starplots(const std::string anaName = "rhcfhc",
54  const int quantN = 0,
55  bool usecosmics = true,
56  bool nh = true)
57 {
58 
59  float minx = 0.35, maxx = 0.65;
60  float minx_zoom = 0.58, maxx_zoom = 0.60;
61  float miny = 2.30, maxy = 2.70;
62  float miny_zoom = 2.43, maxy_zoom = 2.50;
63 
64  if(anaName=="rhc"){
65  twobeams = false;
66  totquant = 4;
67  minx = 0.55; maxx = 0.70;
68  minx_zoom = 0.62; maxx_zoom = 0.68;
69  miny = 2.20; maxy = 3.00;
70  miny_zoom = 2.555;maxy_zoom = 2.65;
71  }
72  if(anaName=="fhc"){
73  horn_one = "fhc";
74  period_one = "full";
75  pot_one = pot_fhc;
77  twobeams = false;
78  totquant = 4;
79  minx = 0.35; maxx = 0.65;
80  minx_zoom = 0.48; maxx_zoom = 0.55;
81  miny = 2.20; maxy = 2.70;
82  miny_zoom = 2.40; maxy_zoom = 2.55;
83  }
84 
85 
86  // calculator
88  RestartCalculator(calc, anaName, nh);
89 
90  // full systematics
91  std::vector< const ISyst* > systs;
92  systs = getAllAna2018Systs(kNumuAna2018,true);
93 
94  std::map< TString, std::vector<const ISyst*> > mapSystLists;
95  mapSystLists["GenieRwgt"] = MyGenieRwgtSysts;
96  mapSystLists["GenieRwgtFr"] = MyGenieRwgtFrSysts;
97  mapSystLists["MEC"] = MyMECSysts;
98  mapSystLists["XSec"] = MyXSecSysts;
99  mapSystLists["GeniePCA"] = MyGeniePCASysts;
100  mapSystLists["FluxSysts"] = MyFluxSysts ;
101  mapSystLists["CalibSysts"] = MyCalibSysts;
102  mapSystLists["OtherSysts"] = MyOtherSysts ;
103 
104  // get contour
105  std::string inFile = "../results/surfprof_realdata_cosmics_"+anaName+"__extrap_stats_2018calc_nh__numu2018.root";
106  TFile* fSurfSens = new TFile(inFile.c_str());
107  FrequentistSurface* surfaceSens = LoadFrom< FrequentistSurface >(fSurfSens,"surface").release();
108  fSurfSens -> Close();
109 
110  double bestFitDm = surfaceSens->GetBestFitY();
111  double bestFitSsq = surfaceSens->GetBestFitX();
112  std::cout<< "Best fit from result contour: ssqth23: " << bestFitSsq
113  << ", dm^2_32: " << bestFitDm <<std::endl;
114  delete fSurfSens;
115 
116 
117  //inputs
118  TString dout_name = "../plots_blessing/";
119  std::string ddata_name = "/nova/ana/nu_mu_ana/Ana2018/Data/";
120  std::string dcosm_name = "/nova/ana/nu_mu_ana/Ana2018/Cosmics/";
121  std::string dpred_name = "/nova/ana/nu_mu_ana/Ana2018/Predictions/";
122  std::string fpred_name = "";
123  std::string fpred_suffix = "";
124 
125 
126  //predictions systs framework
127  const std::string dpred_name_fhc= dpred_name+"pred_numuconcat_fhc__numu2018.root";
128  const std::string dpred_name_rhc= dpred_name+"pred_numuconcat_rhc__numu2018.root";
129  std::string dpred_name_one = dpred_name_rhc;
130  std::string dpred_name_two = dpred_name_fhc;
131  if(anaName=="fhc") dpred_name_one = dpred_name_fhc;
132 
133  std::cout << "\nloading predictions" << std::endl;
134  std::vector<PredictionSystJoint2018*> predictions;
135  ENu2018ExtrapType extrap = kNuMu;
136  for(int quant = 1; quant < totquant+1; quant++){
137  std::cout << "quantile " << quant << std::endl;
138  if(quant<5) predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_one, quant, systs, dpred_name_one));
139  else predictions.push_back(new PredictionSystJoint2018(extrap, calc, horn_two, quant-4, systs, dpred_name_two));
140  }
141 
142 
143  // cosmics
144  std::string cosmics = "nocosmics";
145  std::vector<Spectrum> s_cosmics;
146  if(usecosmics){
148  cosmics = "cosmics";
149  std::cout << "\nloading cosmics" << std::endl;
150  std::cout << "" << horn_one << std::endl;
151  std::string fcosmics_name_one = dcosm_name + "cosmics_" + horn_one + "__numu2018.root";
152  TFile fcosmics_one(fcosmics_name_one.c_str());
153  auto *cosmics_quant1_one = (TH1*)fcosmics_one.Get("cosmics_q1");
154  auto *cosmics_quant2_one = (TH1*)fcosmics_one.Get("cosmics_q2");
155  auto *cosmics_quant3_one = (TH1*)fcosmics_one.Get("cosmics_q3");
156  auto *cosmics_quant4_one = (TH1*)fcosmics_one.Get("cosmics_q4");
157  s_cosmics.push_back(Spectrum (cosmics_quant1_one, pot, livetime));
158  s_cosmics.push_back(Spectrum (cosmics_quant2_one, pot, livetime));
159  s_cosmics.push_back(Spectrum (cosmics_quant3_one, pot, livetime));
160  s_cosmics.push_back(Spectrum (cosmics_quant4_one, pot, livetime));
161  if(twobeams){
163  std::cout << "" << horn_two << std::endl;
164  std::string fcosmics_name_two = dcosm_name + "cosmics_" + horn_two + "__numu2018.root";
165  TFile fcosmics_two(fcosmics_name_two.c_str());
166  auto *cosmics_quant1_two = (TH1*)fcosmics_two.Get("cosmics_q1");
167  auto *cosmics_quant2_two = (TH1*)fcosmics_two.Get("cosmics_q2");
168  auto *cosmics_quant3_two = (TH1*)fcosmics_two.Get("cosmics_q3");
169  auto *cosmics_quant4_two = (TH1*)fcosmics_two.Get("cosmics_q4");
170  s_cosmics.push_back(Spectrum (cosmics_quant1_two, pot, livetime));
171  s_cosmics.push_back(Spectrum (cosmics_quant2_two, pot, livetime));
172  s_cosmics.push_back(Spectrum (cosmics_quant3_two, pot, livetime));
173  s_cosmics.push_back(Spectrum (cosmics_quant4_two, pot, livetime));
174  }
175  }
176 
177 
178 
179  // actually do something
180  for(auto const& thisSystList: mapSystLists){
181 
182  std::string NameZoom = "starplot_"+anaName+"_"+(std::string)thisSystList.first+"__zoom";
183  TCanvas* canvas_zoom = new TCanvas("canvas_zoom", (NameZoom).c_str());
184  TString out_zoom = dout_name+NameZoom.c_str();
185 
186  std::string Name = "starplot_"+anaName+"_"+(std::string)thisSystList.first;
187  TCanvas* canvas = new TCanvas("canvas", (Name).c_str());
188  TString out = dout_name+Name.c_str();
189 
190  TH2F *axes_zoom = new TH2F("setaxes_zoom","", 30, minx_zoom, maxx_zoom, 50, miny_zoom, maxy_zoom);
191  axes_zoom->GetXaxis()->SetTitle("sin^{2}#theta_{23}");
192  axes_zoom->GetYaxis()->SetTitle("#Deltam^{2}_{32} (10^{-3} eV^{2})");
193  axes_zoom->GetXaxis()->SetTitleOffset(0.95);
194  axes_zoom->GetYaxis()->SetTitleOffset(0.85);
195  axes_zoom->GetXaxis()->SetTitleSize(0.055);
196  axes_zoom->GetYaxis()->SetTitleSize(0.055);
197  axes_zoom->GetXaxis()->SetLabelSize(0.04);
198  axes_zoom->GetYaxis()->SetLabelSize(0.04);
199  axes_zoom->GetXaxis()->CenterTitle();
200  axes_zoom->GetYaxis()->CenterTitle();
201  axes_zoom->SetTitle("");
202  axes_zoom->GetYaxis()->SetDecimals();
203 
204  TH2F *axes = new TH2F("setaxes","" ,38 ,minx, maxx, 52, miny, maxy);
205  axes->GetXaxis()->SetTitle("sin^{2}#theta_{23}");
206  axes->GetYaxis()->SetTitle("#Deltam^{2}_{32} (10^{-3} eV^{2})");
207  axes->GetXaxis()->SetTitleOffset(0.95);
208  axes->GetYaxis()->SetTitleOffset(0.85);
209  axes->GetXaxis()->SetTitleSize(0.055);
210  axes->GetYaxis()->SetTitleSize(0.055);
211  axes->GetXaxis()->SetLabelSize(0.04);
212  axes->GetYaxis()->SetLabelSize(0.04);
213  axes->GetXaxis()->CenterTitle();
214  axes->GetYaxis()->CenterTitle();
215  axes->SetTitle("");
216  axes->GetYaxis()->SetDecimals();
217  axes->Draw();
218  axes->Draw();
219 
220 
221  canvas_zoom->cd();
222  gPad->SetFillStyle(0);
223  gPad->SetTopMargin(0.08);
224  gPad->SetRightMargin(0.08);
225  gPad->SetBottomMargin(0.12);
226  gPad->SetLeftMargin(0.12);
227  gStyle->SetTitleOffset(0.95,"x");
228  gStyle->SetTitleOffset(0.85,"y");
229  axes_zoom->Draw();
230  surfaceSens->DrawContour(Gaussian90Percent2D(*surfaceSens), kSolid, kBlack);
231  surfaceSens->DrawBestFit(kBlack, kFullStar);
232  /*
233  TH2* h2d = (TH2*) surfaceSens -> ToTH2();
234  for(int binx = 1; binx<h2d->GetNbinsX(); binx++)
235  for(int biny = 1; biny<h2d->GetNbinsY(); biny++)
236  h2d-> SetBinContent(binx, biny, 0);
237  h2d-> GetYaxis()->SetRangeUser(miny_zoom,maxy_zoom);
238  h2d-> GetXaxis()->SetRangeUser(minx_zoom,maxx_zoom);
239  h2d-> GetYaxis()->CenterTitle();
240  h2d-> GetXaxis()->CenterTitle();
241  h2d-> Draw();
242  */
243 
244 
245  canvas->cd();
246  gPad->SetFillStyle(0);
247  gPad->SetTopMargin(0.08);
248  gPad->SetRightMargin(0.08);
249  gPad->SetBottomMargin(0.12);
250  gPad->SetLeftMargin(0.12);
251  gStyle->SetTitleOffset(0.95,"x");
252  gStyle->SetTitleOffset(0.85,"y");
253  axes->Draw();
254  surfaceSens->DrawContour(Gaussian90Percent2D(*surfaceSens), kSolid, kBlack);
255  surfaceSens->DrawBestFit(kBlack, kFullStar);
256  DrawLegendBFNull(surfaceSens->GetBestFitX(), surfaceSens->GetBestFitY(), kBlack, kFullStar);
257 
258  TLegend *leg = new TLegend(0.40, 0.65, 0.70, 0.85);
259  TLegendEntry *entry = leg->AddEntry("NULL","NOvA NH #nu_{#mu}+#bar{#nu}_{#mu} Systematic","h");
260  /*
261  entry = legend->AddEntry("zero","90% CL Sensitivity","L");
262  entry->SetLineColor(kBlack);
263  entry->SetLineStyle(kSolid);
264  entry->SetLineWidth(3);
265  */
266  int graphId = 0;
267 
268 
269  std::vector<const ISyst*> systList = thisSystList.second;
270 
271  for(const ISyst* s: systList){
272  Progress prog("Making star plot for "+s->ShortName());
273  TGraph* g = new TGraph;
274  TGraph* gLO = new TGraph;
275  for(double sigma = -1; sigma < +1.05; sigma += 0.2){
276 
277  SystShifts systListShift;
278  systListShift.SetShift(s, sigma);
279  RestartCalculator(calc, anaName, nh);
280 
281  // fake data
282  std::cout << "\nfilling fake data" << std::endl;
283  std::vector<Spectrum> s_fakedata;
284  for(int quant = 0; quant < totquant; quant++){
285  std::cout << "quantile " << quant+1 << std::endl;
286  if(quant<4){ pot = pot_one; livetime = livetime_one;}
287  else {pot = pot_two; livetime = livetime_two;}
288  s_fakedata.push_back( predictions[quant]->PredictSyst(calc, systListShift).FakeData(pot) );
289  }
290 
291  // single sample experiments
292  std::vector <const IExperiment*> experiments;
293  if(quantN==0){
294  for(int quant = 0; quant < totquant; quant++){
295  if(usecosmics){
296  s_fakedata[quant] += s_cosmics[quant];
297  experiments.push_back(new SingleSampleExperiment(predictions[quant], s_fakedata[quant], s_cosmics[quant], 0.1));
298  }
299  else experiments.push_back(new SingleSampleExperiment(predictions[quant], s_fakedata[quant]));
300  }
301  }
302  else{
303  if(usecosmics){
304  s_fakedata[quantN-1] += s_cosmics[quantN-1];
305  experiments.push_back(new SingleSampleExperiment(predictions[quantN-1], s_fakedata[quantN-1], s_cosmics[quantN-1], 0.1));
306  }
307  else experiments.push_back(new SingleSampleExperiment(predictions[quantN-1], s_fakedata[quantN-1]));
308  }
309 
310  experiments.push_back(&kSolarConstraintsPDG2017);
311  experiments.push_back(WorldReactorConstraint2017());
312  MultiExperiment exptSyst(experiments);
313 
314  std::vector <const IFitVar*> fitvarsUO = {&kFitSinSqTheta23UO, &kFitDmSq32Scaled};
315  std::vector <const IFitVar*> fitvarsLO = {&kFitSinSqTheta23LO, &kFitDmSq32Scaled};
316 
317  RestartCalculator(calc, anaName, nh);
318  MinuitFitter fitter(&exptSyst, fitvarsUO);
319  fitter.Fit(calc, IFitter::kQuiet);
320 
321  g->SetPoint(g->GetN(),
323  kFitDmSq32Scaled.GetValue(calc));
324 
325  std::cout<< "\nsigma: " << sigma <<std::endl;
326  std::cout << kFitSinSqTheta23.GetValue(calc) << " " << kFitDmSq32Scaled.GetValue(calc) << std::endl;
327 
328  RestartCalculator(calc, anaName, nh);
329  MinuitFitter fitterLO(&exptSyst, fitvarsLO);
330  fitterLO.Fit(calc, IFitter::kQuiet);
331 
332  gLO->SetPoint(gLO->GetN(),
334  kFitDmSq32Scaled.GetValue(calc));
335 
336  std::cout << kFitSinSqTheta23LO.GetValue(calc) << " " << kFitDmSq32Scaled.GetValue(calc) << std::endl;
337  prog.SetProgress((sigma+1)/2.);
338 
339  }
340 
341  g->SetLineColor(color[graphId]);
342  g->SetLineStyle(styleline[graphId]);
343  gLO->SetLineColor(color[graphId]);
344  gLO->SetLineStyle(styleline[graphId]);
345  graphId++;
346 
347  g->SetTitle((s->LatexName()).c_str());
348  g->SetFillColor(10);
349  g->SetMarkerStyle(0);
350  g->SetLineWidth(3);
351  gLO->SetTitle((s->LatexName()).c_str());
352  gLO->SetFillColor(10);
353  gLO->SetMarkerStyle(0);
354  gLO->SetLineWidth(3);
355 
356  leg-> AddEntry(g, (s->LatexName()).c_str(), "l");
357  leg-> SetFillColor(0);
358  leg-> SetFillStyle(0);
359 
360  canvas-> cd();
361  g ->Draw("l");
362  gLO->Draw("l");
363  surfaceSens->DrawBestFit(kBlack,kFullStar);
364 
365  canvas_zoom-> cd();
366  g ->Draw("l");
367  gLO->Draw("l");
368  surfaceSens->DrawBestFit(kBlack,kFullStar);
369  }
370 
371  canvas-> cd();
372  leg->Draw();
373  canvas->SaveAs((out+".root"));
374  canvas->SaveAs((out+".pdf"));
375  canvas->SaveAs((out+".png"));
376  std::ofstream file;
377  file.open (out + ".txt");
378  file << "90% confidence level sensitivity at NOvA's combined neutrino + antineutrino 2018 stats only best fit. Each point in the coloured lines represent the new best fit when a single systematic is applied to the fit: from start to end, a shift to the systematic is applied to cover a +-1$\sigma$ range.";
379  file.close();
380  canvas->Close();
381 
382  canvas_zoom-> cd();
383  leg->Draw();
384  canvas_zoom->SaveAs((out_zoom+".root"));
385  canvas_zoom->SaveAs((out_zoom+".pdf"));
386  canvas_zoom->SaveAs((out_zoom+".png"));
387  std::ofstream file_zoom;
388  file_zoom.open (out_zoom + ".txt");
389  file_zoom << "Zoom in the 90% confidence level sensitivity at NOvA's combined neutrino + antineutrino 2018 stats only best fit. Each point in the coloured lines represent the new best fit when a single systematic is applied to the fit: from start to end, a shift to the systematic is applied to cover a +-1$\sigma$ range.";
390  file_zoom.close();
391  canvas_zoom->Close();
392 
393  }
394 
395 
396 } // End of function
std::vector< const ISyst * > MyMECSysts
enum BeamMode kOrange
std::vector< const ISyst * > MyFluxSysts
std::vector< const ISyst * > MyGeniePCASysts
bin1_2sigma SetFillColor(3)
enum BeamMode kRed
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
double GetValue(const osc::IOscCalcAdjustable *osc) const override
void FakeData()
Definition: rootlogon.C:156
Loads shifted spectra from files.
double pot_two
Definition: saveFDMCHists.C:20
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
double maxy
int totquant
Definition: saveFDMCHists.C:40
T sqrt(T number)
Definition: d0nt_math.hpp:156
const double livetime_fhc
Definition: saveFDMCHists.C:12
std::vector< const ISyst * > MyCalibSysts
const double sin_rhc_stats
Definition: settings.h:9
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
ntuple SetFillStyle(1001)
double pot_one
Definition: saveFDMCHists.C:19
virtual void SetDmsq32(const T &dmsq32)=0
bool twobeams
Definition: saveFDMCHists.C:39
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:40
std::vector< const ISyst * > MyOtherSysts
Log-likelihood scan across two parameters.
ifstream inFile
Definition: AnaPlotMaker.h:34
const XML_Char * s
Definition: expat.h:262
double livetime_two
Definition: saveFDMCHists.C:23
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
std::vector< const ISyst * > MyGenieRwgtSysts
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
leg AddEntry(GRdata,"data","p")
#define pot
const double delta_rhc_stats
Definition: settings.h:12
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
void DrawLegendBFNull(double bfSin, double bfDm, Color_t color, Style_t style)
Definition: numu_tools.h:80
std::vector< const ISyst * > MyGenieRwgtFrSysts
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
Definition: starPlot.C:80
c1 Close()
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
void RestartCalculator(osc::IOscCalcAdjustable *calc, std::string anaName, bool nh)
Combine multiple component experiments.
const FitSinSqTheta23UO kFitSinSqTheta23UO
Definition: starPlot.C:95
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
std::vector< Style_t > styleline
const double pot_fhc
Definition: saveFDMCHists.C:10
double livetime
Definition: saveFDMCHists.C:21
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
double livetime_one
Definition: saveFDMCHists.C:22
const double delta_fhc_stats
Definition: settings.h:11
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
const double sin_rhcfhc_stats
Definition: settings.h:10
TFile * file
Definition: cellShifts.C:17
A simple ascii-art progress bar.
Definition: Progress.h:9
std::string dataName
std::vector< const ISyst * > MyXSecSysts
const FitSinSqTheta23LO kFitSinSqTheta23LO
Definition: starPlot.C:116
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
Definition: starPlot.C:100
enum BeamMode kGreen
enum BeamMode kBlue
std::string period_one
Definition: saveFDMCHists.C:35
def entry(str)
Definition: HTMLTools.py:26
const double delta_rhcfhc_stats
Definition: settings.h:13
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
c cd(1)
std::string horn_one
Definition: saveFDMCHists.C:33
const double sin_fhc_stats
Definition: settings.h:8
std::string horn_two
Definition: saveFDMCHists.C:34
std::vector< Color_t > color
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
void make_starplots(const std::string anaName="rhcfhc", const int quantN=0, bool usecosmics=true, bool nh=true)
Perform MINUIT fits in one or two dimensions.
Definition: MinuitFitter.h:17
enum BeamMode string