bicount_2018.C
Go to the documentation of this file.
1 
4 #include "CAFAna/Fit/Fit.h"
7 #include "CAFAna/Core/Spectrum.h"
16 #include "CAFAna/Core/Progress.h"
17 #include "CAFAna/FC/FCCollection.h"
18 #include "CAFAna/FC/FCPoint.h"
19 #include "CAFAna/FC/FCSurface.h"
20 
23 #include "Utilities/func/MathUtil.h"
24 #include "Utilities/rootlogon.C"
25 
26 #include "joint_fit_2018_tools.h"
27 
28 using namespace ana;
29 
30 #include "TFile.h"
31 #include "TH2.h"
32 #include "TRandom3.h"
33 #include "TMath.h"
34 #include "TMarker.h"
35 #include "TCanvas.h"
36 #include "TGraph.h"
37 #include "TLatex.h"
38 #include "TLegend.h"
39 #include "TColor.h"
40 #include "TLegend.h"
41 #include "TStyle.h"
42 #include "TGraphAsymmErrors.h"
43 #include "TFeldmanCousins.h"
44 #include "TPaveText.h"
45 
46 #include <iostream>
47 #include <iomanip>
48 
49 const double scale1 = 500./750. * 0.9;
50 TCanvas * BiCountCanvas(double minx=0, double maxx=40, double miny=0,double maxy=40,
51  TString title="", TString prelim_text=""){
52  auto c1= new TCanvas(UniqueName().c_str(), UniqueName().c_str(), 700, 750);
53  c1->SetLeftMargin(0.14);
54  c1->SetTopMargin(0.08);
55 // c1->SetBottomMargin(0.15);
56  TH2* axes = new TH2F(UniqueName().c_str(), ";P_{#mue};P_{#bar{#mu}#bar{e}}",
57  100, minx, maxx, 100, miny,maxy);
58  axes->SetTitle(title);
59 
60  axes->GetXaxis()->SetTitle("Total events - neutrino beam");
61  axes->GetYaxis()->SetTitle("Total events - antineutrino beam");
62  axes->GetXaxis()->CenterTitle();
63  axes->GetYaxis()->CenterTitle();
64 
65  axes->GetXaxis()->SetTitleSize(kBlessedTitleFontSize*scale1);
66  axes->GetYaxis()->SetTitleSize(kBlessedTitleFontSize*scale1);
67  axes->GetXaxis()->SetLabelSize(kBlessedLabelFontSize*scale1);
68  axes->GetYaxis()->SetLabelSize(kBlessedLabelFontSize*scale1);
69  // axes->GetXaxis()->SetTitleOffset(0.8);
70  axes->GetYaxis()->SetTitleOffset(0.95);
71 
72  axes->Draw();
73 
74  TLatex* prelim = new TLatex(.9, .95, ("NOvA " + prelim_text).Data());
75  prelim->SetTextColor(prelim_text=="Preliminary"? kBlue:kGray+1);
76  prelim->SetNDC();
77  prelim->SetTextSize(2/30.*scale1);
78  prelim->SetTextAlign(32);
79  if(prelim_text!="") prelim->Draw();
80 
81  return c1;
82 }
83 
84 TMarker * GetMarker_dCP (int delta, Int_t color)
85 {
86  TMarker * mark = new TMarker;
87  mark->SetMarkerColor(color);
88  Int_t style;
89  if(delta==0) style = kOpenCircle;
90  if(delta==1) style = kFullCircle;
91  if(delta==2) style = kOpenSquare;
92  if(delta==3) style = kFullSquare;
93  mark->SetMarkerStyle(style);
94  return mark;
95 }
96 
98  const IPrediction * predFHC, const IPrediction * predRHC,
99  double potFHC, double potRHC,
100  double cosmFHC=0, double cosmRHC=0,
101  TString label="", Int_t color=1, Int_t style = kSolid,
102  bool markers=true)
103 {
104  TGraph* g = new TGraph;
105  double avgx = 0, avgy = 0;
106  for(int dIdx = 0; dIdx <= 100; ++dIdx){
107  double delta = dIdx*0.02;
108  calc->SetdCP(delta*M_PI);
109  const double P = predFHC->Predict(calc).Integral(potFHC) + cosmFHC;
110  const double Pbar = predRHC->Predict(calc).Integral(potRHC) + cosmRHC;
111  g->SetPoint(g->GetN(), P, Pbar);
112 
113  if(markers){
114  if(dIdx%25==0)
115  GetMarker_dCP(dIdx/25,color)->DrawMarker(P,Pbar);
116  }
117  avgx += P;
118  avgy += Pbar;
119  }
120  avgx=avgx/100.;
121  avgy=avgy/100.;
122  TLatex* ltx = new TLatex(avgx, avgy,label);
123  ltx->SetTextAlign(22);
124  ltx->SetTextSize(0.04);
125  ltx->SetTextColor(color);
126  ltx->Draw();
127 
128  g->SetLineWidth(2);
129  g->SetLineColor(color);
130  g->SetLineStyle(style);
131  g->Draw("l");
132  return g;
133 }
134 
135 TGraph* GetDataPoint(double x, double y)
136 {
137  auto gr = new TGraphAsymmErrors();
138  gr->SetPoint(0,x,y);
139  gr->SetMarkerStyle(kFullCircle);
140  gr->SetLineWidth(2);
141  gr->SetMarkerColor(kBlack);
142  TFeldmanCousins fc(0.6827);//1 sigma
143  if ( y < 50 ) gr->SetPointEYlow(0, y-fc.CalculateLowerLimit(y,0));
144  else gr->SetPointEYlow(0,sqrt(y));
145  if ( y < 30 ) gr->SetPointEYhigh(0,fc.CalculateUpperLimit(y,0)-y);
146  else gr->SetPointEYhigh(0,sqrt(y));
147 
148  if ( x < 50 ) gr->SetPointEXlow(0, x-fc.CalculateLowerLimit(x,0));
149  else gr->SetPointEXlow(0,sqrt(x));
150  if ( x < 30 ) gr->SetPointEXhigh(0,fc.CalculateUpperLimit(x,0)-x);
151  else gr->SetPointEXhigh(0,sqrt(x));
152 
153  return gr;
154 }
155 
156 
157 TLegend * Legend_dCP(double x1, double y1,double x2, double y2)
158 {
159  auto leg = new TLegend(x1,y1,x2,y2);
160 // leg->SetFillStyle(1001);
161  // leg->SetFillColor(kWhite);
162  leg->SetFillStyle(0);
163  leg->SetTextSize(scale1*kBlessedLabelFontSize);
164  leg->AddEntry(GetMarker_dCP(0,1),"#delta_{CP}= 0","p");
165  leg->AddEntry(GetMarker_dCP(1,1),"#delta_{CP}= #pi/2","p");
166  leg->AddEntry(GetMarker_dCP(2,1),"#delta_{CP}= #pi","p");
167  leg->AddEntry(GetMarker_dCP(3,1),"#delta_{CP}= 3#pi/2","p");
168  leg->SetNColumns(2);
169  return leg;
170 }
171 
172 TLegend * Legend_BFV(double x1, double y1, double x2, double y2, TMarker* mark, TString text){
173  auto leg = new TLegend(x1,y1,x2,y2);
174  leg->SetFillStyle(0);
175  leg->SetTextSize(scale1*kBlessedLabelFontSize);
176  leg->SetTextColor(mark->GetMarkerColor());
177  leg->AddEntry(mark, text,"p");
178  return leg;
179 }
180 
181 TLatex * HieLabel(int hie, Int_t color, double x, double y)
182 {
183  auto ltx = new TLatex();
184  ltx->SetTextSize(scale1*kBlessedLabelFontSize);
185  ltx->SetTextColor(color);
186  ltx->SetTextAlign(22);
187  if(hie>0) ltx->DrawLatex(x,y,"#splitline{ NH}{#Deltam_{32}^{2}=#plus2.50#times10^{-3}eV^{2}}");
188  if(hie<0) ltx->DrawLatex(x,y,"#splitline{ IH}{#Deltam_{32}^{2}=#minus2.55#times10^{-3}eV^{2}}");
189 // if(hie>0) ltx->DrawLatex(x,y,"#splitline{ NH}{#Deltam_{32}^{2}=#plus2.45#times10^{-3}eV^{2}}");
190 // if(hie<0) ltx->DrawLatex(x,y,"#splitline{ IH}{#Deltam_{32}^{2}=#minus2.51#times10^{-3}eV^{2}}");
191 // if(hie<0) ltx->DrawLatex(x,y,"#Deltam_{32}^{2}<0");
192 
193  return ltx;
194 }
195 TLatex * OctLabel(int oct, Int_t color, double x, double y)
196 {
197  auto ltx = new TLatex();
198  ltx->SetTextSize(scale1*kBlessedLabelFontSize);
199  ltx->SetTextColor(color);
200  ltx->SetTextAlign(22);
201  if(oct>0) ltx->DrawLatex(x,y,"#splitline{ UO}{sin^{2}#theta_{23}=0.59}");
202  if(oct<0) ltx->DrawLatex(x,y,"#splitline{ LO}{sin^{2}#theta_{23}=0.46}");
203  // if(oct<0) ltx->DrawLatex(x,y,"#splitline{ IH}{#Deltam_{32}^{2}=#minus2.55#times10^{-3}eV^{2}}");
204 // if(hie>0) ltx->DrawLatex(x,y,"#splitline{ NH}{#Deltam_{32}^{2}=#plus2.45#times10^{-3}eV^{2}}");
205 // if(hie<0) ltx->DrawLatex(x,y,"#splitline{ IH}{#Deltam_{32}^{2}=#minus2.51#times10^{-3}eV^{2}}");
206 // if(hie<0) ltx->DrawLatex(x,y,"#Deltam_{32}^{2}<0");
207 
208  return ltx;
209 }
210 TPaveText * PaveParams (double x1, double y1, double x2, double y2)
211 {
212  auto pave = new TPaveText(x1,y1,x2,y2,"NDC NB");
213  pave->SetTextAlign(33); //right top
214  pave->SetFillStyle(0);
215  pave->SetLineColorAlpha(0,0);
216  pave->SetTextSize(scale1*kBlessedLabelFontSize);
217  //pave->AddText("|#Deltam_{32}^{2}|=2.67#times10^{-3}eV^{2}");
218  pave->AddText("sin^{2}2#theta_{13}=0.082");
219  // pave->AddText("sin^{2}#theta_{23}=0.46,0.59");
220 // pave->AddText("sin^{2}#theta_{23}=0.47,0.56");
221  return pave;
222 }
223 
224 TPaveText* PavePOT(double x1, double y1, double x2, double y2)
225 {
226  auto pave = new TPaveText(x1,y1,x2,y2,"NDC NB");
227  pave->SetTextAlign(13); //;eft top
228  pave->SetFillStyle(0);
229  pave->SetLineColorAlpha(0,0);
230  pave->SetTextSize(scale1*kBlessedLabelFontSize);
231  pave->AddText("NOvA FD");
232  pave->AddText("9.48#times10^{20} POT (#nu)");
233  pave->AddText("6.91#times10^{20} POT (#bar{#nu})");
234  return pave;
235 }
236 
237 
238 TMarker * BestFitMarker (double x, double y, int color = kViolet-5){
239  auto bfmarker = new TMarker (x,y,kStar);
240  bfmarker->SetMarkerStyle(kFullStar);
241  bfmarker->SetMarkerColor(color);
242  bfmarker->SetMarkerSize(3);
243  return bfmarker;
244 }
245 
246 
247 
249 {
250 
251  TString folder = "plots_prob/";
252 
253  double potFHC = kAna2018FHCPOT;
254  double potRHC = kAna2018RHCPOT;
255 
256  std::string decomp_FHC = "combo";
257  std::string decomp_RHC = "prop";
258 
259  auto predFHC = GetNuePrediction2018(decomp_FHC, DefaultOscCalc(), {}, "fhc", false);
260  auto predRHC = GetNuePrediction2018(decomp_RHC, DefaultOscCalc(), {}, "rhc", false);
261 
262  double cosmFHC = 3.33;
263  double cosmRHC = 0.71;
264 
265  for(TString opt : {"basic","bfv","data","bfv_data"}){
266  auto ca = BiCountCanvas(20,80,2,28,"", opt.Contains("data")?"Preliminary":"Simulation");
267 
268  auto calc= DefaultOscCalc();
269 
270  for(auto hie:{-1,+1}){
271  if(hie>0) calc->SetDmsq32( 2.502e-3);
272  if(hie<0) calc->SetDmsq32(-2.552e-3);
273 
274  Int_t theoryColor = (hie>0?k3SigmaConfidenceColorNH:kPrimColorIH);
275 
276  for(auto ssth:{0.461,0.586}){
277  calc->SetTh23(asin(sqrt(ssth)));
278  auto g1 = BicountEllipse_dCP(calc,
279  predFHC, predRHC, potFHC, potRHC,
280  3.33, 0.71,
281  /*TString::Format("%.1lf", ssth)*/"",theoryColor);
282  }
283  }
284 
285  Legend_dCP(0.14,0.12,0.52,0.23)->Draw();
286  HieLabel(-1,kPrimColorIH,35,12);
288  OctLabel(-1,kBlack,30,15);
289  OctLabel(+1,kBlack,72,17);
290  PaveParams(0.52,0.83,0.89,0.92)->Draw();
291  PavePOT(0.15,0.78,0.5,0.92)->Draw();
292 
293  if(opt.Contains("bfv")) {
294  SetFakeCalc(calc, 0.586, 2.502e-3, 0.186);
295  auto bfmarker = BestFitMarker(predFHC->Predict(calc).Integral(potFHC) + cosmFHC,
296  predRHC->Predict(calc).Integral(potRHC) + cosmRHC);
297 
298  bfmarker->Draw();
299  Legend_BFV(0.55, 0.12, 0.89,0.17, bfmarker, "2018 best fit")->Draw();
300  }
301 
302  if(opt.Contains("data")){
303  auto datapoint1 = GetDataPoint(58,18);
304  datapoint1->Draw("pe");
305  }
306  gPad->Print(folder + "bicount_"+opt + ".pdf");
307  }
308 
309 // SetFakeCalc(calc, 0.558, 2.445e-3, 1.21);
310 // auto bf2017 = BestFitMarker(predFHC->Predict(calc).Integral(potFHC) + cosmFHC,
311 // predRHC->Predict(calc).Integral(potRHC) + cosmRHC,
312 // kGray+1);
313 // // bf2017->Draw();
314 // // Legend_BFV(0.55, 0.17, 0.89,0.22, bf2017, "2017 best fit")->Draw();
315 //
316 // SetFakeCalc(calc, 0.538, 2.494e-3, 1.3);
317 // auto nufit = BestFitMarker(predFHC->Predict(calc).Integral(potFHC) + cosmFHC,
318 // predRHC->Predict(calc).Integral(potRHC) + cosmRHC,
319 // kGray+1);
320 // nufit->SetMarkerStyle(22);
321 // nufit->SetMarkerSize(2);
322 // // nufit->Draw();
323 // // Legend_BFV(0.55, 0.22, 0.89,0.27, nufit, "NuFit 2018")->Draw();
324 //
325 // gPad->Print(folder + "bicount_bfv_data.pdf");
326 // gPad->Print(folder + "bicount_bfv_data_2.pdf");
327  return;
328 }
329 
const Color_t kPrimColorIH
Definition: Style.h:55
osc::OscCalculatorDumb calc
Oscillation analysis framework, runs over CAF files outside of ART.
osc::IOscCalculatorAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
Definition: Calcs.cxx:49
TPaveText * PavePOT(double x1, double y1, double x2, double y2)
Definition: bicount_2018.C:224
Float_t y1[n_points_granero]
Definition: compare.C:5
TCanvas * BiCountCanvas(double minx=0, double maxx=40, double miny=0, double maxy=40, TString title="", TString prelim_text="")
Definition: bicount_2018.C:50
double delta
Definition: runWimpSim.h:98
double maxy
Float_t x1[n_points_granero]
Definition: compare.C:5
T sqrt(T number)
Definition: d0nt_math.hpp:156
TPaveText * PaveParams(double x1, double y1, double x2, double y2)
Definition: bicount_2018.C:210
virtual void SetdCP(const T &dCP)=0
TLegend * Legend_BFV(double x1, double y1, double x2, double y2, TMarker *mark, TString text)
Definition: bicount_2018.C:172
void bicount_2018()
Definition: bicount_2018.C:248
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
Definition: Spectrum.cxx:742
#define M_PI
Definition: SbMath.h:34
list markers
Definition: styles.py:7
virtual Spectrum Predict(osc::IOscCalculator *calc) const =0
const double kAna2018RHCPOT
Definition: Exposures.h:208
const char * label
#define P(a, b, c, d, e, x)
General interface to any calculator that lets you set the parameters.
const Color_t k3SigmaConfidenceColorNH
Definition: Style.h:69
const IPrediction * GetNuePrediction2018(std::string decomp, osc::IOscCalculator *calc, bool corrSysts, std::string beam, bool isFakeData, bool GetFromUPS=false, bool minimizeMemory=false, bool NERSC=false)
const double scale1
Definition: bicount_2018.C:49
TLatex * HieLabel(int hie, Int_t color, double x, double y)
Definition: bicount_2018.C:181
TLatex * prelim
Definition: Xsec_final.C:133
TMarker * BestFitMarker(double x, double y, int color=kViolet-5)
Definition: bicount_2018.C:238
TGraph * BicountEllipse_dCP(osc::IOscCalculatorAdjustable *calc, const IPrediction *predFHC, const IPrediction *predRHC, double potFHC, double potRHC, double cosmFHC=0, double cosmRHC=0, TString label="", Int_t color=1, Int_t style=kSolid, bool markers=true)
Definition: bicount_2018.C:97
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TGraph * GetDataPoint(double x, double y)
Definition: bicount_2018.C:135
c1
Definition: demo5.py:24
Standard interface to all prediction techniques.
Definition: IPrediction.h:57
void SetFakeCalc(osc::IOscCalculatorAdjustable *calc, double ssth23=-5, double dmsq32=-5, double dCP_Pi=-5)
TLatex * OctLabel(int oct, Int_t color, double x, double y)
Definition: bicount_2018.C:195
const Float_t kBlessedLabelFontSize
Definition: Style.h:81
const double kAna2018FHCPOT
Definition: Exposures.h:207
Float_t e
Definition: plot.C:35
const Float_t kBlessedTitleFontSize
Definition: Style.h:80
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:30
TLegend * Legend_dCP(double x1, double y1, double x2, double y2)
Definition: bicount_2018.C:157
T asin(T number)
Definition: d0nt_math.hpp:60
TMarker * GetMarker_dCP(int delta, Int_t color)
Definition: bicount_2018.C:84