starPlot.C
Go to the documentation of this file.
6 #include "CAFAna/Core/Loaders.h"
8 #include "CAFAna/Core/Progress.h"
10 #include "CAFAna/Core/Spectrum.h"
12 #include "CAFAna/Cuts/Cuts.h"
14 #include "CAFAna/Cuts/SpillCuts.h"
15 #include "CAFAna/Cuts/TimingCuts.h"
30 #include "CAFAna/Vars/FitVars.h"
34 #include "CAFAna/Vars/Vars.h"
37 #include "OscLib/OscCalcPMNSOpt.h"
38 
39 using namespace ana;
40 
41 #include "TCanvas.h"
42 #include "TFile.h"
43 #include "TGraph.h"
44 #include "TGraphAsymmErrors.h"
45 #include "TH1.h"
46 #include "TH2.h"
47 #include "TLatex.h"
48 #include "TLegend.h"
49 #include "TStyle.h"
50 
51 void Preliminary() // For some reason it doesn't include it automatically
52 {
53  TLatex* prelim = new TLatex(.9, .95, "NOvA Preliminary");
54  prelim->SetTextColor(kBlue);
55  prelim->SetNDC();
56  prelim->SetTextSize(2/30.);
57  prelim->SetTextAlign(32);
58  prelim->Draw();
59 }
60 
62 {
64  calc->SetDmsq32(2.44006e-3);
65  calc->SetTh23(asin(sqrt(0.558965)));
66  calc->SetL(810);
67  calc->SetRho(0); // No matter effects
68  calc->SetDmsq21(7.59e-5);
69  calc->SetTh12(.601);
70  calc->SetTh13(.1567);
71  calc->SetdCP(0);
72 
73 }
74 
75 
76 // These two clamp th23 to the uo / lo
78 {
79 public:
80  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const
81  {
82  return util::sqr(sin(osc->GetTh23()));
83  };
84  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const
85  {
86  osc->SetTh23(asin(sqrt(Clamp(val))));
87  };
88  virtual std::string ShortName() const {return "sinsqth23_UO";}
89  virtual std::string LatexName() const {return "sin^{2}#theta_{23} Constrained to Upper Octant";}
90 
91  // And now set max/min values for upper octant boundaries
92  virtual double LowLimit() const {return 0.514;}
93  virtual double HighLimit() const {return 1.;}
94 };
96 
98 {
99 public:
100  virtual double GetValue(const osc::IOscCalcAdjustable* osc) const
101  {
102  return util::sqr(sin(osc->GetTh23()));
103  };
104  virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const
105  {
106  osc->SetTh23(asin(sqrt(Clamp(val))));
107  };
108  virtual std::string ShortName() const {return "sinsqth23_LO";}
109  virtual std::string LatexName() const {return "sin^{2}#theta_{23} Constraint to Lower Octant";}
110 
111  // And now set max/min values for lower octant boundaries
112  virtual double LowLimit() const {return 0.;}
113  //virtual double HighLimit() const {return 0.514;}
114  virtual double HighLimit() const {return 0.49;}
115 };
117 
118 
119 
120 void starPlot()
121 {
122 
123  gStyle->SetPadLeftMargin(0.12);
124  gStyle->SetTitleOffset(1.15,"x");
125  gStyle->SetTitleOffset(0.95,"y");
126 
127 
130 
131  // make syst list of the largest systs
132  std::vector<const ISyst*> systs_XSec;
133  std::vector<const ISyst*> systs_Beam;
134  std::vector<const ISyst*> systs_FileSysts;
135  std::vector<const ISyst*> systs_Others;
136  AddJointAna2017XSecSysts(systs_XSec, kNumuAna2017, false);
138  AddJointAna2017FileSysts(systs_FileSysts, kNumuAna2017);
140 
141  std::map< TString, std::vector<const ISyst*> > mapSystLists;
142  mapSystLists["XSec"] = systs_XSec;
143  mapSystLists["Beam"] = systs_Beam;
144  mapSystLists["FileSysts"] = systs_FileSysts;
145  mapSystLists["Others"] = systs_Others;
146 
147  ENumu2017ExtrapType extrap = kNuMu;
148  PredictionSystNumu2017 prediction_Q1(extrap, &calc, 1);
149  PredictionSystNumu2017 prediction_Q2(extrap, &calc, 2);
150  PredictionSystNumu2017 prediction_Q3(extrap, &calc, 3);
151  PredictionSystNumu2017 prediction_Q4(extrap, &calc, 4);
152 
153  // get the contour
154  TFile* fSurfSens = new TFile("resultCont__allSysts.root"); // result
155  FrequentistSurface* surfaceSens = LoadFrom< FrequentistSurface >(fSurfSens, "surface").release();
156  fSurfSens -> Close();
157  double bestFitDm = surfaceSens -> fBestFitY;
158  double bestFitSsq = surfaceSens -> fBestFitX;
159  std::cout<< "Best fit from result contour: ssqth23: " << bestFitSsq
160  << ", dm^2_32: " << bestFitDm <<std::endl;
161  delete fSurfSens;
162 
163 
164  for(auto const& thisSystList: mapSystLists){
165 
166  std::string NameZoom = "starPlotZoom" + (std::string)thisSystList.first;
167  TCanvas* cZoom = new TCanvas("cZoom",(NameZoom).c_str());
168  cZoom->SetBottomMargin(0.15);
169  TH2* h2d = (TH2*) surfaceSens -> ToTH2();
170  for(int binx = 1; binx<h2d->GetNbinsX(); binx++)
171  for(int biny = 1; biny<h2d->GetNbinsY(); biny++)
172  h2d -> SetBinContent(binx, biny, 0);
173  h2d -> GetYaxis()->SetRangeUser(2.4,2.48);
174  h2d -> GetXaxis()->SetRangeUser(0.455,0.48);
175  h2d -> GetYaxis()->CenterTitle();
176  h2d -> GetXaxis()->CenterTitle();
177  h2d -> Draw();
178 
179  std::string Name = "starPlot" + (std::string)thisSystList.first;
180  TCanvas* c = new TCanvas("c",(Name).c_str());
181  c->SetBottomMargin(0.15);
182  surfaceSens->DrawContour(Gaussian90Percent2D(*surfaceSens), kSolid, kBlack);
183  surfaceSens->DrawBestFit(kBlack);
184 
185 
186  TLegend *leg = new TLegend(0.15, 0.65, 0.60, 0.85);
187  int nOfGraphs = 0;
188 
189 
190  std::vector<const ISyst*> systList = thisSystList.second;
191 
192  for(const ISyst* s: systList){
193  Progress prog("Making star plot for "+s->ShortName());
194  TGraph* g = new TGraph;
195  TGraph* gLO = new TGraph;
196  for(double sigma = -1; sigma < +1.05; sigma += 0.2){
197 
198  SystShifts systListShift;
199  systListShift.SetShift(s, sigma);
200  RestartCalculator(&calc); // Restart &calculator after fitting
201 
202  // Systematically shifted fake data to fit to
203  Spectrum fakeData_Q1 = prediction_Q1.PredictSyst(&calc, systListShift).FakeData(kAna2017POT);
204  Spectrum fakeData_Q2 = prediction_Q2.PredictSyst(&calc, systListShift).FakeData(kAna2017POT);
205  Spectrum fakeData_Q3 = prediction_Q3.PredictSyst(&calc, systListShift).FakeData(kAna2017POT);
206  Spectrum fakeData_Q4 = prediction_Q4.PredictSyst(&calc, systListShift).FakeData(kAna2017POT);
207 
208  SingleSampleExperiment expt1(&prediction_Q1, fakeData_Q1);
209  SingleSampleExperiment expt2(&prediction_Q2, fakeData_Q2);
210  SingleSampleExperiment expt3(&prediction_Q3, fakeData_Q3);
211  SingleSampleExperiment expt4(&prediction_Q4, fakeData_Q4);
212 
213  std::vector<const IExperiment*> vConstraints;
214  vConstraints.push_back(&expt1);
215  vConstraints.push_back(&expt2);
216  vConstraints.push_back(&expt3);
217  vConstraints.push_back(&expt4);
218  vConstraints.push_back(&kSolarConstraintsPDG2017);
219  vConstraints.push_back(WorldReactorConstraint2017());
220  MultiExperiment exptSyst(vConstraints);
221 
222  MinuitFitter fitter(&exptSyst, {&kFitSinSqTheta23UO, &kFitDmSq32Scaled});
223  fitter.Fit(&calc, IFitter::kQuiet);
224 
225  g->SetPoint(g->GetN(),
226  kFitSinSqTheta23UO.GetValue(&calc),
227  kFitDmSq32Scaled.GetValue(&calc));
228 
229  std::cout<< "\nsigma: " << sigma <<std::endl;
231 
232 
233  MinuitFitter fitterLO(&exptSyst, {&kFitSinSqTheta23LO, &kFitDmSq32Scaled});
234  fitterLO.Fit(&calc, IFitter::kQuiet);
235 
236  gLO->SetPoint(gLO->GetN(),
237  kFitSinSqTheta23LO.GetValue(&calc),
238  kFitDmSq32Scaled.GetValue(&calc));
239 
240  std::cout << kFitSinSqTheta23LO.GetValue(&calc) << " " << kFitDmSq32Scaled.GetValue(&calc) << std::endl;
241 
242 
243 
244  prog.SetProgress((sigma+1)/2.);
245 
246  }
247  if(nOfGraphs<=12)
248  {
249  // This would still need some playing with the output to make it look pretty,
250  // but at least it's not a bunch of graphs all looking the same
251  g->SetLineColor(7-nOfGraphs/2); // Hopefully these are different enough
252  g->SetLineStyle(nOfGraphs%2 + 1); // Alternating dashed and normal line
253 
254  gLO->SetLineColor(7-nOfGraphs/2); // Hopefully these are different enough
255  gLO->SetLineStyle(nOfGraphs%2 + 1); // Alternating dashed and normal line
256  }
257  else
258  {
259  g->SetLineColor(nOfGraphs%10);
260  g->SetLineStyle((nOfGraphs-10)%10);
261 
262  gLO->SetLineColor(nOfGraphs%10);
263  gLO->SetLineStyle((nOfGraphs-10)%10);
264  }
265 
266  nOfGraphs++;
267 
268 
269  g->SetTitle((s->LatexName()).c_str());
270  g->SetFillColor(10);
271  g->SetMarkerStyle(0);
272  g->SetLineWidth(3);
273  gLO->SetTitle((s->LatexName()).c_str());
274  gLO->SetFillColor(10);
275  gLO->SetMarkerStyle(0);
276  gLO->SetLineWidth(3);
277 
278  leg->AddEntry(g, (s->LatexName()).c_str(), "l");
279  leg -> SetFillColor(0);
280  leg -> SetFillStyle(0);
281 
282  c -> cd();
283  g->Draw("l");
284  gLO->Draw("l");
285  cZoom -> cd();
286  g->Draw("l");
287  gLO->Draw("l");
288 
289  }
290 
291  c -> cd();
292  Preliminary();
293  leg->Draw();
294 
295  c->SaveAs(("plots/"+Name+".root").c_str());
296  c->SaveAs(("plots/"+Name+".pdf").c_str());
297  c->SaveAs(("plots/"+Name+".png").c_str());
298  c->Close();
299 
300 
301  cZoom -> cd();
302  Preliminary();
303  leg->Draw();
304 
305  cZoom->SaveAs(("plots/"+NameZoom+".root").c_str());
306  cZoom->SaveAs(("plots/"+NameZoom+".pdf").c_str());
307  cZoom->SaveAs(("plots/"+NameZoom+".png").c_str());
308  cZoom->Close();
309 
310  }
311 
312 } // End of function
tree Draw("slc.nhit")
void Preliminary()
Definition: starPlot.C:51
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
Definition: starPlot.C:84
virtual void SetL(double L)=0
bin1_2sigma SetFillColor(3)
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
void starPlot()
Definition: starPlot.C:120
void AddJointAna2017OtherSysts(std::vector< const ISyst * > &systs, const EAnaType2017 ana)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
correl_xv GetYaxis() -> SetDecimals()
virtual void SetDmsq21(const T &dmsq21)=0
Simple record of shifts applied to systematic parameters.
Definition: SystShifts.h:20
T sqrt(T number)
Definition: d0nt_math.hpp:156
Spectrum PredictSyst(osc::IOscCalc *calc, const SystShifts &shift) const override
virtual void SetTh13(const T &th13)=0
virtual std::string LatexName() const
Definition: starPlot.C:109
virtual std::string ShortName() const
Definition: starPlot.C:108
virtual T GetTh23() const
Definition: IOscCalc.h:89
const double kAna2017POT
Definition: Exposures.h:174
osc::OscCalcDumb calc
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
virtual std::string LatexName() const
Definition: starPlot.C:89
int Clamp(int x)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
Log-likelihood scan across two parameters.
const XML_Char * s
Definition: expat.h:262
void DrawBestFit(Color_t color, Int_t marker=kFullCircle) const
Draw the best fit point.
Definition: ISurface.cxx:32
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:402
virtual double LowLimit() const
Definition: starPlot.C:92
correl_xv GetXaxis() -> SetDecimals()
void DrawContour(TH2 *fc, Style_t style, Color_t color, double minchi=-1)
Definition: ISurface.cxx:44
Loads shifted spectra from files.
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:64
Optimized version of OscCalcPMNS.
Definition: StanTypedefs.h:31
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
void RestartCalculator(osc::IOscCalcAdjustable *calc)
Definition: starPlot.C:61
void AddJointAna2017XSecSysts(std::vector< const ISyst * > &systs, const EAnaType2017 ana, bool smallgenie)
const ConstrainedFitVarWithPrior fitDmSq32Scaled_UniformPrior & kFitDmSq32Scaled
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
Definition: starPlot.C:80
TLatex * prelim
Definition: Xsec_final.C:133
c1 Close()
Oscillation probability calculators.
Definition: Calcs.h:5
double sigma(TH1F *hist, double percentile)
OStream cout
Definition: OStream.cxx:6
cout<< "--"<< endl;for(Int_t iP=1;iP<=hyz->GetNbinsX();iP++){for(Int_t iC=1;iC<=hyz->GetNbinsY();iC++){if(hyv->GetBinContent(iP, iC)>-999){goal_hyv-> SetBinContent(iP, iC,-(dy[iP-1][iC-1]))
virtual double HighLimit() const
Definition: starPlot.C:93
virtual double LowLimit() const
Definition: starPlot.C:112
Combine multiple component experiments.
const FitSinSqTheta23UO kFitSinSqTheta23UO
Definition: starPlot.C:95
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
virtual void SetRho(double rho)=0
const SolarConstraints kSolarConstraintsPDG2017(7.53e-5, 0.18e-5, 0.851, 0.020)
T sin(T number)
Definition: d0nt_math.hpp:132
void AddJointAna2017FileSysts(std::vector< const ISyst * > &systs, const EAnaType2017 ana)
void AddJointAna2017BeamSysts(std::vector< const ISyst * > &systs, const EAnaType2017 ana)
virtual std::string ShortName() const
Definition: starPlot.C:88
TH2 * Gaussian90Percent2D(const FrequentistSurface &s)
Up-value surface for 90% confidence in 2D in gaussian approximation.
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
A simple ascii-art progress bar.
Definition: Progress.h:9
const FitSinSqTheta23LO kFitSinSqTheta23LO
Definition: starPlot.C:116
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
Definition: starPlot.C:100
Float_t e
Definition: plot.C:35
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
Definition: starPlot.C:104
void SetShift(const ISyst *syst, double shift, bool force=false)
Definition: SystShifts.cxx:78
virtual double HighLimit() const
Definition: starPlot.C:114
c cd(1)
TH2 * ToTH2(const Spectrum &s, double exposure, ana::EExposureType expotype, const Binning &binsx, const Binning &binsy, ana::EBinType bintype)
For use with Var2D.
Definition: Utilities.cxx:256
virtual void SetTh12(const T &th12)=0
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