joint_fit_style_tools.h
Go to the documentation of this file.
1 #pragma once
2 
8 //#include "CAFAna/Prediction/PredictionNoExtrap.h"
14 #include "CAFAna/Systs/XSecSysts.h"
15 #include "CAFAna/Vars/FitVars.h"
16 #include "CAFAna/Analysis/Style.h"
18 
19 #include "Utilities/func/MathUtil.h"
20 
21 #include "TCanvas.h"
22 #include "TBox.h"
23 #include "TLatex.h"
24 #include "TColor.h"
25 #include "TGraph.h"
26 #include "TVectorD.h"
27 #include "TF1.h"
28 #include "TH2.h"
29 #include "TLegend.h"
30 #include "TLine.h"
31 #include "TMarker.h"
32 #include "TStyle.h"
33 #include "TSystem.h"
34 #include "TGaxis.h"
35 
36 using namespace ana;
37 
38 #include "OscLib/IOscCalc.h"
39 
40 #include "TFile.h"
41 
43  double dmsq32 = -5, double dCP_Pi = -5)
44  {
45  if(dCP_Pi > -1) kFitDeltaInPiUnits.SetValue(calc, dCP_Pi);
46  if(ssth23 > -1) kFitSinSqTheta23.SetValue(calc, ssth23);
47  if(dmsq32 > -1) kFitDmSq32.SetValue(calc, dmsq32);
48  }
49 
51 {
52  TH2* temp = (TH2*)h->Clone();
53 
54  int NX = h->GetXaxis()->GetNbins();
55  int NY = h->GetYaxis()->GetNbins();
56 
57  // assert(NX % 2 == 0 &&
58  // "Histogram must have an even number of X bins for this to work");
59 
60  // Save the location we expect the end column of the original
61  // TH2 to be after we do the folding
62  int foldedBin = 1 + NX/2;
63 
64  // Fold the initial histogram to join the ends of the X axis
65  for(int i = 1; i <= NX/2; i++) {
66  for(int j = 1; j <= NY/2; j++) {
67  // x1, y1 denotes a position in the left-hand side half plane
68  // x2, y2 denotes a position in the right-hand side half plane
69  int x1 = i+NX/2;
70  int x2 = i%(NX/2+1)+2;
71  int y1 = j;
72  int y2 = j;
73  double z1 = h->GetBinContent(x1, y1);
74  double z2 = h->GetBinContent(x2, y2);
75  temp->SetBinContent(x1, y1, z2);
76  temp->SetBinContent(x2, y2, z1);
77  }
78  }
79  // smooth both histograms
80  temp->Smooth();
81  h->Smooth();
82 
83  // Replace end columns
84  for(int y = 1; y <= NY; y++) {
85  double sZ1 = temp->GetBinContent(foldedBin, y);
86  double sZ2 = temp->GetBinContent(foldedBin, y);
87  h->SetBinContent(1, y, sZ1);
88  h->SetBinContent(NX, y, sZ2);
89  }
90 }
91 
92 //////////////////////////////////////////////////////////////////////
93 
95  const double pot, Spectrum* cosmics = 0, const double livetime = 0)
96 {
97  // TODO - why does this need to be a pointer?
98  Spectrum* ret = new Spectrum(pred->Predict(calc).FakeData(pot, livetime));
99  if(cosmics) *ret += *cosmics;
100  return ret;
101 }
102 
104  const double pot, TH1D* cosmics = 0)
105 {
106  // TODO - why does this need to be a pointer?
107  Spectrum ret = pred->Predict(calc);
108  if(cosmics){
109  std::cout << "joint_fit_style_tools.h: GetMockData() - no way to set cosmics livetime. Can we convert this to Spectrum too please? And in any case, the cosmics should be fluctuated too, so the old code was definitely wrong" << std::endl;
110  abort();
111  // ret += *cosmics;
112  }
113  return new Spectrum(ret.MockData(pot)); // TODO why a pointer?
114 }
115 
117  const IPrediction* pred,
118  Spectrum* cosmics,
119  double futurePOT,
120  double futureLT,
121  bool randomSystShifts = false,
122  bool asimov=false)
123 {
124  Spectrum total = pred->Predict(calc).FakeData(futurePOT, futureLT) + *cosmics;
125  if (randomSystShifts){
126  std::cout << "GenerateFutureData(): randomSystShifts not implemented" << std::endl;
127  abort();
128  }
129 
130  if(asimov) return total.FakeData(futurePOT); else return total.MockData(futurePOT);
131 }
132 // Moved these to CAFAna/Vars/FitVars.*
133 ////////////////////////////////////////////////////////////////////////
134 //
135 //// For slices plots, plot by true octant, and make separate plots of each
136 //class FitSinSqTheta23LowerOctant: public IConstrainedFitVar
137 //{
138 //public:
139 // FitSinSqTheta23LowerOctant() : IConstrainedFitVar("ssth23LO", "sin^{2}#theta_{23} LO") {}
140 // virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
141 // return util::sqr(sin(osc->GetTh23()));};
142 // virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
143 // osc->SetTh23(asin(sqrt(Clamp(val))));};
144 // //virtual std::string ShortName() const {return "ssth23LO";}
145 // //virtual std::string LatexName() const {return "sin^{2}#theta_{23} LO";}
146 // virtual double LowLimit() const {return 0;}
147 // virtual double HighLimit() const {return 0.5;}
148 //};
149 //
150 //const IConstrainedFitVar * kFitSinSqTheta23LowerOctant = new FitSinSqTheta23LowerOctant();
151 //
152 //class FitSinSqTheta23UpperOctant: public IConstrainedFitVar
153 //{
154 //public:
155 // FitSinSqTheta23UpperOctant() : IConstrainedFitVar("ssth23UO", "sin^{2}#theta_{23} UO") {}
156 // virtual double GetValue(const osc::IOscCalcAdjustable* osc) const {
157 // return util::sqr(sin(osc->GetTh23()));};
158 // virtual void SetValue(osc::IOscCalcAdjustable* osc, double val) const {
159 // osc->SetTh23(asin(sqrt(Clamp(val))));};
160 // //virtual std::string ShortName() const {return "ssth23UO";}
161 // //virtual std::string LatexName() const {return "sin^{2}#theta_{23} UO";}
162 // virtual double LowLimit() const {return 0.5;}
163 // virtual double HighLimit() const {return 1;}
164 //};
165 //const IConstrainedFitVar * kFitSinSqTheta23UpperOctant = new FitSinSqTheta23UpperOctant();
166 //
167 
168 ////////////////////////////////////////////////////////////////////////////
169 //////////////////////// Slice, Contours Style /////////////////////////////
170 ///////////////////////////////////////////////////////////////////////////
171 //move functions from plotting script here
172 
173 void DrawSliceCanvas(TString slicename, const double ymax,
174  const double xmin = 0, const double xmax = 2., bool overlay = false, bool t2kunits = false, bool isFuture = false){
175  auto c1 = new TCanvas(ana::UniqueName().c_str());
176  c1->SetFillStyle(0);
177  c1->SetLeftMargin(0.14);
178  c1->SetBottomMargin(0.15);
179  TString title;
180  if(slicename.Contains("delta")) {
181  title = ";delta;Significance (#sqrt{#Delta#chi^{2}})";
182  }
183  if(slicename.Contains("th23")) {
184  title = ";sin^{2}#theta_{23};Significance (#sqrt{#Delta#chi^{2}})";
185  }
186  if(slicename.Contains("dmsq")) {
187  title = ";|#Deltam^{2}_{32}| (10^{-3} eV^{2});Significance (#sqrt{#Delta#chi^{2}})";
188  }
189  if(slicename.Contains("th13")) {
190  title = ";sin^{2}2#theta_{13};Significance (#sqrt{#Delta#chi^{2}})";
191  }
192 
193  if (isFuture) {
194  if ( slicename.Contains("dcp" ) ) { title = ";#delta_{CP}/#pi;Significance of CP Violation (#sigma)"; }
195  else if ( slicename.Contains("mh" ) ) { title = ";#delta_{CP}/#pi;Significance of hierarchy resolution (#sigma)"; }
196  else if ( slicename.Contains("maxmix") ) { title = ";sin^{2}#theta_{23};Significance of maximal mixing rejection (#sigma)"; }
197  else if ( slicename.Contains("oct" ) ) { title = ";sin^{2}#theta_{23};Significance of octant determination (#sigma)"; }
198  }
199 
200  auto axes = new TH2F("",title, 100, xmin, xmax, 100, 0, ymax);
201  if(slicename.Contains("fccorr") || overlay) axes->GetYaxis()->SetTitle("Significance (#sigma)");
202  // CenterTitles(axes);
203  axes->GetXaxis()->CenterTitle();
204  axes->GetYaxis()->CenterTitle();
205  axes->Draw();
206  if (!isFuture) {
207  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
208  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
209  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
210  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
211  } else {
212  if ( slicename.Contains("maxmix") || slicename.Contains("oct") )
213  axes->GetYaxis()->SetTitleSize( 0.8*axes->GetYaxis()->GetTitleSize() );
214  }
215  axes->GetXaxis()->SetTitleOffset(0.8);
216  axes->GetYaxis()->SetTitleOffset(0.75);
217  if(slicename.Contains("delta")) XAxisDeltaCPLabels(axes, t2kunits);
218  c1->RedrawAxis();
219  }
220 
221 TLegend * SliceLegend(std::vector <std::pair <TGraph*, TString > > & graphs, bool isDelta, double low=0.7, bool isDmsq = false, bool overlay = false,
222  bool isSyst = false, bool isFuture = false) {
223  TLegend * leg;
224  double xmin = 0.18, ymin = 0.2, xmax = 0.4 , ymax = 0.49;
225  if(isDelta) { xmin = 0.6 , ymin = 0.6, xmax = 0.9 , ymax = 0.89; }
226  if(isFuture){ xmin = 0.45, ymin = 0.6, xmax = 0.65, ymax = 0.85; }
227  leg = new TLegend( xmin, ymin, xmax, ymax );
228 
229  leg->SetTextSize(0.8*kBlessedLabelFontSize);
230  leg->SetFillStyle(0);
231  if (isFuture) {
232  TString str="#splitline{Full-dataset Sensitivity}";
233  if (isSyst) str+="{2020 Syst.}";
234  else str+="{Stat. Only}";
235  leg->SetHeader(str,"C");
236  }
237  for (auto &gr:graphs){
238  TString str="";
239  if(gr.second.Contains("O")){
240  if(isDelta) leg->SetX1(0.57);
241  if(isDmsq) leg->SetX1(0.18);
242  if(!isDelta && !isDmsq) {leg->SetX1(0.6); leg->SetX2(0.89);}
243  leg->SetMargin(0.18);
244  str += gr.second.Contains("NH") ? "NH":"IH";
245  str += gr.second.Contains("LO") ? " Lower": " Upper";
246  str += " octant";
247  }
248  else{
249  if (!isFuture) str="#splitline{";
250  str += gr.second.Contains("NH") ? "Normal":"Inverted";
251  if (!isFuture) str += "}{hierarchy}";
252  if(overlay) {
253  str = "Uncorr. ";
254  str += gr.second.Contains("NH") ? "NH" : "IH";
255  }
256  }
257  leg->AddEntry(gr.first, str, "l");
258  }
259  return leg;
260 }
261 
262 
263 void DrawContourCanvas (TString surfName, double minx = 0,
264  double maxx = 2, double miny = 0, double maxy = 1, bool t2kunits = false){
265  auto c1 = new TCanvas(ana::UniqueName().c_str());
266  c1->SetFillStyle(0);
267  c1->SetLeftMargin(0.14);
268  c1->SetBottomMargin(0.15);
269  TH2* axes = new TH2F();
270  TString title;
271  if(surfName.Contains("delta_th23")) title=";#delta_{CP};sin^{2}#theta_{23}";
272  if(surfName.Contains("th13_delta")) title=";sin^{2}2#theta_{13};#delta_{CP}/#pi";
273  if(surfName.Contains("th23_dm32")) title=";sin^{2}#theta_{23};#Deltam^{2}_{32} (10^{-3}eV^{2})";
274  axes = new TH2F("",title,100,minx,maxx,100,miny,maxy);
275 // CenterTitles(axes);
276  axes->GetXaxis()->CenterTitle();
277  axes->GetYaxis()->CenterTitle();
278  axes->Draw();
279  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize);
280  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize);
281  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize);
282  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize);
283  axes->GetXaxis()->SetTitleOffset(0.8);
284  axes->GetYaxis()->SetTitleOffset(0.8);
285  TGaxis::SetMaxDigits(3);
286  if(surfName.Contains("th23_dm32")) {
287  axes->GetYaxis()->SetDecimals() ;
288  axes->GetYaxis()->SetTitleOffset(0.85);
289  axes->GetYaxis()->SetLabelOffset(0.002);
290  axes->GetYaxis()->SetLabelSize(0.058);
291  axes->GetYaxis()->SetTitleSize(0.078);
292  axes->GetXaxis()->SetTitleSize(.95*kBlessedTitleFontSize);
293  axes->GetXaxis()->SetTitleOffset(0.86);
294  }
295  if(surfName.Contains("delta_th23")) XAxisDeltaCPLabels(axes, t2kunits);
296  c1->RedrawAxis();
297 }
298 
299 
300 TLegend * ContourLegend(int hie, bool fillContour, bool fccorr,
301  Int_t kFillColor1,Int_t kFillColor2, Int_t kFillColor3,Int_t kDarkColor,
302  bool bestFit){
303  TLegend * leg = new TLegend(0.16,fccorr?0.19:0.17,0.65,0.28);
304  leg->SetTextSize(kBlessedLabelFontSize);
305  leg->SetFillStyle(0);
306  leg->SetMargin(1.3*leg->GetMargin());
307  leg->SetNColumns(3);
308 
309  TGraph * dummy = new TGraph;
310  dummy->SetLineWidth(3);
311 
312  if(fillContour){
313  dummy->SetLineColor(kDarkColor);
314  dummy->SetFillColor(kFillColor1);
315  leg->AddEntry(dummy->Clone(),"1#sigma","f");
316  dummy->SetFillColor(kFillColor2);
317  leg->AddEntry(dummy->Clone(),"2#sigma","f");
318  dummy->SetFillColor(kFillColor3);
319  leg->AddEntry(dummy->Clone(),"3#sigma","f");
320  }
321  else{
322  leg->SetMargin(1.4*leg->GetMargin());
323  dummy->SetLineColor(kFillColor1);
324  dummy->SetLineStyle(kSolid);
325  leg->AddEntry(dummy->Clone(),"1#sigma","l");
326  dummy->SetLineStyle(k2SigmaConfidenceStyle);
327  dummy->SetLineColor(kFillColor2);
328  leg->AddEntry(dummy->Clone(),"2#sigma","l");
329  dummy->SetLineStyle(kSolid);
330  dummy->SetLineColor(kFillColor3);
331  leg->AddEntry(dummy->Clone(),"3#sigma","l");
332  }
333  if(bestFit){
334  leg->SetNColumns(4);
335  //leg->SetTextSize(0.8*kBlessedLabelFontSize);
336 // leg->SetColumnSeparation(0.08);
337  leg->SetX2(0.75);
338  dummy->SetMarkerStyle(43);
339  dummy->SetMarkerSize(2);
340  dummy->SetMarkerColor(kDarkColor);
341  dummy->SetLineColor(kDarkColor);
342  dummy->SetLineStyle(7);
343  if(hie>0) leg->AddEntry(dummy->Clone(),"Best Fit","p");
344  if(hie<0) leg->AddEntry((TObject*)0, "#color[0]{Best Fit}", "");
345  dummy->SetLineStyle(1);
346  }
347  if(!fccorr) leg->SetHeader("No Feldman-Cousins");
348  else leg->SetHeader("Feldman-Cousins");
349  return leg;
350 }
351 
352 
353 TLegend * FCOverlayContourLegend(int hie, bool fillContour,
354  Int_t kFillColor1,Int_t kFillColor2, Int_t kFillColor3,Int_t kDarkColor)
355 {
356  TLegend * leg = new TLegend(0.16,hie>0?.78:0.77,0.65,0.86);
357  leg->SetTextSize(kBlessedLabelFontSize);
358  leg->SetFillStyle(0);
359  leg->SetMargin(1.3*leg->GetMargin());
360  leg->SetNColumns(3);
361 
362  TGraph * dummy = new TGraph;
363  dummy->SetLineWidth(3);
364 
365  if(fillContour){
366  dummy->SetLineColor(kDarkColor);
367  dummy->SetFillColor(kFillColor1);
368  leg->AddEntry(dummy->Clone(),"1#sigma","f");
369  dummy->SetFillColor(kFillColor2);
370  leg->AddEntry(dummy->Clone(),"2#sigma","f");
371  dummy->SetFillColor(kFillColor3);
372  leg->AddEntry(dummy->Clone(),"3#sigma","f");
373  }
374  else{
375  leg->SetMargin(1.4*leg->GetMargin());
376  dummy->SetLineColor(kFillColor1);
377  dummy->SetLineStyle(kSolid);
378  leg->AddEntry(dummy->Clone(),"1#sigma","l");
379  dummy->SetLineStyle(k2SigmaConfidenceStyle);
380  dummy->SetLineColor(kFillColor2);
381  leg->AddEntry(dummy->Clone(),"2#sigma","l");
382  dummy->SetLineStyle(kSolid);
383  dummy->SetLineColor(kFillColor3);
384  leg->AddEntry(dummy->Clone(),"3#sigma","l");
385  }
386 
387  leg->SetHeader("No Feldman-Cousins");
388  return leg;
389 }
390 
391 void latex(double x, double y, std::string txt, double ang = 0, int align = 22, double size = 2/30.) {
392  TLatex* ltx = new TLatex(x, y, txt.c_str());
393  ltx->SetNDC();
394  ltx->SetTextSize(size);
395  ltx->SetTextAlign(align);
396  ltx->SetTextAngle(ang);
397  cout<<"i can draw!"<<endl;
398  ltx->Draw();
399 }
400 
401  int cnh = TColor::GetFreeColorIndex();
402  int cih = TColor::GetFreeColorIndex() + 1;
403  TColor *color_ih = new TColor(cih, 255/255.0, 134/255.0, 122/255.0);
404  TColor *color_nh = new TColor(cnh, 0/255.0, 167/255.0, 200/255.0);
405 
406  int cnh_dark = TColor::GetFreeColorIndex();
407  int cih_dark = TColor::GetFreeColorIndex() + 1;
408  TColor *color_ih_dark = new TColor(cih_dark, 243/255.0, 91/255.0, 83/255.0);
409  TColor *color_nh_dark = new TColor(cnh_dark, 0/255.0, 167/255.0, 200/255.0);
410 
Simple analyzer module to print to text file information required for c++ based block alignment study...
int cih_dark
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
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
TColor * color_nh_dark
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Definition: FitVars.cxx:177
void XAxisDeltaCPLabels(TH1 *axes, bool t2kunits)
Label the x-axis with fractions of pi.
Spectrum * GetFakeData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, Spectrum *cosmics=0, const double livetime=0)
int cnh_dark
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var&#39;s SetValue()
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2., bool overlay=false, bool t2kunits=false, bool isFuture=false)
osc::OscCalcDumb calc
Spectrum * GetMockData(const IPrediction *pred, osc::IOscCalc *calc, const double pot, TH1D *cosmics=0)
Spectrum GenerateFutureData(osc::IOscCalc *calc, const IPrediction *pred, Spectrum *cosmics, double futurePOT, double futureLT, bool randomSystShifts=false, bool asimov=false)
Representation of a spectrum in any variable, with associated POT.
Definition: Spectrum.h:33
Double_t ymax
Definition: plot.C:25
void DrawContourCanvas(TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1, bool t2kunits=false)
Spectrum FakeData(double pot) const
Fake data is a MC spectrum scaled to the POT expected in the data.
Definition: Spectrum.cxx:402
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
#define pot
TColor * color_ih_dark
const double j
Definition: BetheBloch.cxx:29
TColor * color_nh
TLegend * ContourLegend(int hie, bool fillContour, bool fccorr, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor, bool bestFit)
OStream cout
Definition: OStream.cxx:6
void SetFakeCalc(osc::IOscCalcAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
Spectrum MockData(double pot, int idx=0) const
Mock data is FakeData with Poisson fluctuations applied.
Definition: Spectrum.cxx:384
TColor * color_ih
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta, double low=0.7, bool isDmsq=false, bool overlay=false, bool isSyst=false, bool isFuture=false)
const Style_t k2SigmaConfidenceStyle
Definition: Style.h:71
void CylindricalSmoothing(TH2 *h)
double dmsq32
double livetime
Definition: saveFDMCHists.C:21
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
TLegend * FCOverlayContourLegend(int hie, bool fillContour, Int_t kFillColor1, Int_t kFillColor2, Int_t kFillColor3, Int_t kDarkColor)
c1
Definition: demo5.py:24
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
Double_t ymin
Definition: plot.C:24
const Float_t kBlessedLabelFontSize
Definition: Style.h:90
const Float_t kBlessedTitleFontSize
Definition: Style.h:89
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
void latex(double x, double y, std::string txt, double ang=0, int align=22, double size=2/30.)