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