demoStarPlots.C
Go to the documentation of this file.
1 ////////////////////////////////////////////////////////////////////////////////
2 //
3 // Another method of evaluating the impact of a systematic
4 //
5 // Taken largely from macros by Luke, and Diana
6 // I've trimmed it down a bit for simplicity
7 // if you want further examples like how to include cosmics
8 // they can be found here:
9 // CAFAna/numu/Analysis2018/BoxOpening/make_starplots.C
10 //
11 // The idea here is to apply a systematic shift then run a fit that is ignorant
12 // of the systematics. If we didn't account for it in our fit but it was there
13 // in data where would it pull the fit.
14 //
15 ////////////////////////////////////////////////////////////////////////////////
16 
17 #include "CAFAna/Analysis/Calcs.h"
19 #include "CAFAna/Fit/Fit.h"
20 #include "CAFAna/Analysis/Plots.h"
27 //#include "CAFAna/Prediction/PredictionSystJointDemo.h"
30 #include "CAFAna/Systs/Systs.h"
31 #include "CAFAna/Vars/FitVars.h"
32 #include "OscLib/IOscCalc.h"
33 
34 #include "TCanvas.h"
35 #include "TGraph.h"
36 #include "TLegend.h"
37 #include "TStyle.h"
38 
39 
40 
41 #include <fstream>
42 #include <iostream>
43 
44 using namespace ana;
45 
46 // style choices
47 std::vector<Style_t> styleline = {kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid, kSolid};
48 std::vector<Color_t> color = {kMagenta, kGreen+1, kOrange+1, kBlue-7, kCyan+1, kRed-4, kBlue+1, kGray+1, kOrange+3};
49 
50 // Helper function to organize the syst groups
51 std::map< std::string, std::vector<const ISyst*> > GetSystMap(){
52  std::map< std::string, std::vector<const ISyst*> > mapSystVecs;
53 
54  std::vector<const ISyst*> muESysts;
55  muESysts.push_back(&kMuEScaleSyst2017);
56  muESysts.push_back(&kRelMuEScaleSyst2017);
57  std::vector<const ISyst*> neutronSysts;
58  neutronSysts.push_back(&kNeutronVisEScalePrimariesSyst2018);
59 
60  mapSystVecs["MuESysts"] = muESysts;
61  mapSystVecs["NeutronSysts"] = neutronSysts;
62 
63  return mapSystVecs;
64 }
65 
66 // And another helper to reset the calculator
67 void ResetCalc (osc::IOscCalcAdjustable * calc, double ssth23 = -5, double dmsq32 = -5)
68 {
70  calc->SetTh23(asin(sqrt(ssth23)));
71  calc->SetDmsq32(dmsq32);
72 }
73 
74 // This should probably go in a header
75 // These two clamp th23 to the uo / lo
77 {
78  public:
79  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const
80  {
81  return util::sqr(sin(osc->GetTh23()));
82  };
83  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const
84  {
85  osc->SetTh23(asin(sqrt(Clamp(val))));
86  };
87  virtual std::string ShortName() const {return "sinsqth23_UO";}
88  virtual std::string LatexName() const {return "sin^{2}#theta_{23} Constrained to Upper Octant";}
89 
90  // And now set max/min values for upper octant boundaries
91  virtual double LowLimit() const {return 0.514;}
92  virtual double HighLimit() const {return 1.;}
93 };
95 
97 {
98  public:
99  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const
100  {
101  return util::sqr(sin(osc->GetTh23()));
102  };
103  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const
104  {
105  osc->SetTh23(asin(sqrt(Clamp(val))));
106  };
107  virtual std::string ShortName() const {return "sinsqth23_LO";}
108  virtual std::string LatexName() const {return "sin^{2}#theta_{23} Constraint to Lower Octant";}
109 
110  // And now set max/min values for lower octant boundaries
111  virtual double LowLimit() const {return 0.;}
112  //virtual double HighLimit() const {return 0.514;}
113  virtual double HighLimit() const {return 0.49;}
114 };
116 
117 
118 
119 
121 {
122  /////////////////////////////////////////////////////////////
123  //
124  // all the setup stuff
125  //
126  /////////////////////////////////////////////////////////////
127  double pot = kAna2018FHCPOT;
128  double livetime = kAna2018FHCLivetime;
129 
130  // some fit points
131  const double ssth23 = 0.585;
132  const double dmsq32 = 2.501e-3;
133 
134  // for later use when plotting
135  float minx = 0.35, maxx = 0.65;
136  float minx_zoom = 0.575, maxx_zoom = 0.595;
137  float miny = 2.30, maxy = 2.70;
138  float miny_zoom = 2.44, maxy_zoom = 2.52;
139 
140  auto calc = DefaultOscCalc();
141  ResetCalc(calc, ssth23, dmsq32);
142 
143  std::vector< const ISyst* > systs = {};
144  auto systmap = GetSystMap();
145  for(auto thissyst : systmap){
146  for(unsigned int sysVecIt = 0; sysVecIt < thissyst.second.size(); ++sysVecIt){
147  systs.push_back(thissyst.second[sysVecIt]);
148  }
149  }
150 
151  /////////////////////////////////////////////////////////////
152  //
153  // Loading in predictions and surface if you want it
154  //
155  /////////////////////////////////////////////////////////////
156  std::vector< PredictionSystJoint2018* > predictions;
157  for(int quant = 1; quant < 5; ++quant){
158  predictions.push_back(new PredictionSystJoint2018(kNuMu, calc, "fhc", quant, systs, "/nova/ana/nu_mu_ana/Ana2018/Predictions/pred_numuconcat_fhc__numu2018.root"));
159  }
160 
161 /*
162  // get contour and best fit points
163  TFile* fSurfSens = new TFile("/nova/ana/nu_mu_ana/Ana2018/Results/surfprof_fakedata_cosmics_fhc__extrap_stats_2018calc_nh__numu2018.root");
164  Surface* surfaceSens = LoadFrom< Surface >(fSurfSens -> GetDirectory("surface")).release();
165  fSurfSens -> Close();
166 
167  double bestFitDm = surfaceSens-> GetMinY();
168  double bestFitSsq = surfaceSens-> GetMinX();
169  delete fSurfSens;
170 */
171 
172  /////////////////////////////////////////////////////////////
173  //
174  // actually do something
175  //
176  /////////////////////////////////////////////////////////////
177  for(auto const& thisSystList: systmap){
178  std::string NameZoom = "sysfit_starplot_" + (std::string)thisSystList.first + "__zoom";
179  TCanvas* canvas_zoom = new TCanvas("canvas_zoom", (NameZoom).c_str());
180  TString out_zoom = NameZoom.c_str();
181 
182  std::string Name = "sysfit_starplot_"+(std::string)thisSystList.first;
183  TCanvas* canvas = new TCanvas("canvas", (Name).c_str());
184  TString out = Name.c_str();
185 
186  // lots of beautification stuff that you can skip over
187  TH2F *axes_zoom = new TH2F("setaxes_zoom","", 30, minx_zoom, maxx_zoom, 50, miny_zoom, maxy_zoom);
188  axes_zoom->GetXaxis()->SetTitle("sin^{2}#theta_{23}");
189  axes_zoom->GetYaxis()->SetTitle("#Deltam^{2}_{32} (10^{-3} eV^{2})");
190  axes_zoom->GetXaxis()->SetTitleOffset(0.95);
191  axes_zoom->GetYaxis()->SetTitleOffset(0.85);
192  axes_zoom->GetXaxis()->SetTitleSize(0.055);
193  axes_zoom->GetYaxis()->SetTitleSize(0.055);
194  axes_zoom->GetXaxis()->SetLabelSize(0.04);
195  axes_zoom->GetYaxis()->SetLabelSize(0.04);
196  axes_zoom->GetXaxis()->CenterTitle();
197  axes_zoom->GetYaxis()->CenterTitle();
198  axes_zoom->SetTitle("");
199  axes_zoom->GetYaxis()->SetDecimals();
200 
201  TH2F *axes = new TH2F("setaxes","" ,38 ,minx, maxx, 52, miny, maxy);
202  axes->GetXaxis()->SetTitle("sin^{2}#theta_{23}");
203  axes->GetYaxis()->SetTitle("#Deltam^{2}_{32} (10^{-3} eV^{2})");
204  axes->GetXaxis()->SetTitleOffset(0.95);
205  axes->GetYaxis()->SetTitleOffset(0.85);
206  axes->GetXaxis()->SetTitleSize(0.055);
207  axes->GetYaxis()->SetTitleSize(0.055);
208  axes->GetXaxis()->SetLabelSize(0.04);
209  axes->GetYaxis()->SetLabelSize(0.04);
210  axes->GetXaxis()->CenterTitle();
211  axes->GetYaxis()->CenterTitle();
212  axes->SetTitle("");
213  axes->GetYaxis()->SetDecimals();
214  axes->Draw();
215  axes->Draw();
216 
217 
218  canvas_zoom->cd();
219  gPad->SetFillStyle(0);
220  gPad->SetTopMargin(0.08);
221  gPad->SetRightMargin(0.08);
222  gPad->SetBottomMargin(0.12);
223  gPad->SetLeftMargin(0.12);
224  gStyle->SetTitleOffset(0.95,"x");
225  gStyle->SetTitleOffset(0.85,"y");
226  axes_zoom->Draw();
227 // surfaceSens->DrawContour(Gaussian90Percent2D(*surfaceSens), kSolid, kBlack);
228 // surfaceSens->DrawBestFit(kBlack, kFullStar);
229 
230  canvas->cd();
231  gPad->SetFillStyle(0);
232  gPad->SetTopMargin(0.08);
233  gPad->SetRightMargin(0.08);
234  gPad->SetBottomMargin(0.12);
235  gPad->SetLeftMargin(0.12);
236  gStyle->SetTitleOffset(0.95,"x");
237  gStyle->SetTitleOffset(0.85,"y");
238  axes->Draw();
239 // surfaceSens->DrawContour(Gaussian90Percent2D(*surfaceSens), kSolid, kBlack);
240 // surfaceSens->DrawBestFit(kBlack, kFullStar);
241 // DrawLegendBFNull(surfaceSens -> GetMinX(), surfaceSens -> GetMinY(), kBlack, kFullStar);
242 
243  TLegend *leg = new TLegend(0.125, 0.45, 0.60, 0.85);
244  TLegendEntry *entry = leg->AddEntry("NULL","NOvA NH #nu_{#mu}+#bar{#nu}_{#mu} Systematic","h");
245  int graphId = 0;
246 
247  std::vector<const ISyst*> systList = thisSystList.second;
248 
249 
250  /////////////////////////////////////////////////////////////
251  //
252  // now really do something
253  //
254  /////////////////////////////////////////////////////////////
255  for(const ISyst* s: systList){
256  // Two graphs one for upper octant and one for lower
257  TGraph* gUO = new TGraph;
258  TGraph* gLO = new TGraph;
259  for(double sigma = -1; sigma < +1.05; sigma += 0.2){
260  SystShifts systListShift;
261  systListShift.SetShift(s, sigma);
262  ResetCalc(calc, ssth23, dmsq32);
263 
264  // fake data and experiments
265  std::vector<Spectrum> s_fakedata;
266  std::vector <const IExperiment*> experiments;
267  for(int quant = 0; quant < 4; quant++){
268  s_fakedata.push_back( predictions[quant]->PredictSyst(calc, systListShift).FakeData(pot) );
269  experiments.push_back(new SingleSampleExperiment(predictions[quant], s_fakedata.back()));
270  }
271  experiments.push_back(&kSolarConstraintsPDG2017);
272  experiments.push_back(WorldReactorConstraint2017());
273  MultiExperiment exptSyst(experiments);
274 
275  std::vector <const IFitVar*> fitvarsUO = {&kFitSinSqTheta23UO, &kFitDmSq32Scaled};
276  std::vector <const IFitVar*> fitvarsLO = {&kFitSinSqTheta23LO, &kFitDmSq32Scaled};
277 
278  // Fit the UO and put in graph
279  ResetCalc(calc, ssth23, dmsq32);
280  Fitter fitter(&exptSyst, fitvarsUO/*, {}, Fitter::kCareful*/);
281  fitter.Fit(calc, Fitter::kQuiet);
282  gUO->SetPoint(gUO->GetN(), kFitSinSqTheta23UO.GetValue(calc), kFitDmSq32Scaled.GetValue(calc));
283 
284  // Fit the LO and put in graph
285  ResetCalc(calc, ssth23, dmsq32);
286  Fitter fitterLO(&exptSyst, fitvarsLO/*, {}, Fitter::kCareful*/);
287  fitterLO.Fit(calc, Fitter::kQuiet);
288  gLO->SetPoint(gLO->GetN(), kFitSinSqTheta23LO.GetValue(calc), kFitDmSq32Scaled.GetValue(calc));
289 
290 
291  }// end of loop through sigmas
292 
293  // More beautification
294  gUO->SetLineColor(color[graphId]);
295  gUO->SetLineStyle(styleline[graphId]);
296  gLO->SetLineColor(color[graphId]);
297  gLO->SetLineStyle(styleline[graphId]);
298  graphId++;
299 
300  gUO->SetTitle((s->LatexName()).c_str());
301  gUO->SetFillColor(10);
302  gUO->SetMarkerStyle(0);
303  gUO->SetLineWidth(3);
304  gLO->SetTitle((s->LatexName()).c_str());
305  gLO->SetFillColor(10);
306  gLO->SetMarkerStyle(0);
307  gLO->SetLineWidth(3);
308 
309  leg-> AddEntry(gUO, (s->LatexName()).c_str(), "l");
310  leg-> SetFillColor(0);
311  leg-> SetFillStyle(0);
312 
313  canvas-> cd();
314  gUO ->Draw("l");
315  gLO->Draw("l");
316 // surfaceSens->DrawBestFit(kBlack,kFullStar);
317 
318  canvas_zoom-> cd();
319  gUO ->Draw("l");
320  gLO->Draw("l");
321 // surfaceSens->DrawBestFit(kBlack,kFullStar);
322 
323  }// end of loop over individual systs
324 
325  // Draw and save
326  canvas-> cd();
327  leg->Draw();
328  canvas->SaveAs((Name+".png").c_str());
329  canvas->Close();
330 
331  canvas_zoom-> cd();
332  leg->Draw();
333  canvas_zoom->SaveAs((NameZoom+".png").c_str());
334  canvas_zoom->Close();
335 
336  }// end of loop over systmap
337 
338 }
enum BeamMode kOrange
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
Definition: demoStarPlots.C:83
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
double ssth23
std::vector< SystGroupDef > systs
Definition: syst_header.h:385
Loads shifted spectra from files.
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
double maxy
T sqrt(T number)
Definition: d0nt_math.hpp:156
virtual std::string LatexName() const
virtual std::string ShortName() const
virtual T GetTh23() const
Definition: IOscCalc.h:94
std::vector< Color_t > color
Definition: demoStarPlots.C:48
const FitSinSqTheta23LO kFitSinSqTheta23LO
TCanvas * canvas(const char *nm, const char *ti, const char *a)
Definition: pass1_plots.C:36
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
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
Encapsulate code to systematically shift a caf::SRProxy.
Definition: ISyst.h:14
ntuple SetFillStyle(1001)
virtual void SetDmsq32(const T &dmsq32)=0
void demoStarPlots()
virtual std::string LatexName() const
Definition: demoStarPlots.C:88
int Clamp(int x)
const MuEScaleSyst2017 kMuEScaleSyst2017(0.0074, 0.0012)
const XML_Char * s
Definition: expat.h:262
std::vector< Style_t > styleline
Definition: demoStarPlots.C:47
const RelMuEScaleSyst2017 kRelMuEScaleSyst2017(0.0045, 10.5)
std::map< std::string, std::vector< const ISyst * > > GetSystMap()
Definition: demoStarPlots.C:51
void FakeData()
Definition: rootlogon.C:156
virtual double LowLimit() const
Definition: demoStarPlots.C:91
const FitSinSqTheta23UO kFitSinSqTheta23UO
Definition: demoStarPlots.C:94
leg AddEntry(GRdata,"data","p")
#define pot
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
Definition: demoStarPlots.C:79
Oscillation probability calculators.
Definition: Calcs.h:5
double sigma(TH1F *hist, double percentile)
void ResetCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5)
Definition: demoStarPlots.C:67
virtual double HighLimit() const
Definition: demoStarPlots.C:92
virtual double LowLimit() const
Combine multiple component experiments.
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
T sin(T number)
Definition: d0nt_math.hpp:132
double dmsq32
virtual std::string ShortName() const
Definition: demoStarPlots.C:87
double livetime
Definition: saveFDMCHists.C:21
virtual void SetTh23(const T &th23)=0
const double kAna2018FHCPOT
Definition: Exposures.h:207
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
Definition: demoStarPlots.C:99
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
enum BeamMode kGreen
enum BeamMode kBlue
def entry(str)
Definition: HTMLTools.py:26
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:81
virtual double HighLimit() const
c cd(1)
const double kAna2018FHCLivetime
Definition: Exposures.h:213
T asin(T number)
Definition: d0nt_math.hpp:60
Compare a single data spectrum to the MC + cosmics expectation.
enum BeamMode string