nova_t2k_comparison.C
Go to the documentation of this file.
3 #include "CAFAna/Fit/Fit.h"
11 #include "CAFAna/FC/FCSurface.h"
16 #include "CAFAna/Vars/FitVars.h"
17 #include "Utilities/rootlogon.C"
18 #include "OscLib/IOscCalc.h"
19 
20 #include "../joint_fit_2020_loader_tools.h"
21 #include "../joint_fit_2020_datarelease_tools.h"
22 
23 #include "TCanvas.h"
24 #include "TMarker.h"
25 #include "TBox.h"
26 #include "TLatex.h"
27 #include "TColor.h"
28 #include "TGraph.h"
29 #include "TVectorD.h"
30 #include "TF1.h"
31 #include "TLegend.h"
32 #include "TLine.h"
33 #include "TMarker.h"
34 #include "TStyle.h"
35 #include "TSystem.h"
36 #include "TGaxis.h"
37 
38 #include <algorithm>
39 #include <vector>
40 #include <string>
41 
42 using namespace ana;
43 
44 TGraph *JointGraphs(TGraph *a, TGraph *b)
45 {
46  TGraph* ret = (TGraph*)a->Clone(UniqueName().c_str());
47  for(int i = 0; i < b->GetN(); ++i){
48  double x, y;
49  b->GetPoint(i, x, y);
50  ret->SetPoint(ret->GetN(), x, y);
51  }
52  return ret;
53 }
54 
55 TGraph* VectorsToGraph(std::vector<double> vx, std::vector<double> vy){
56  double x[vx.size()]; double y[vy.size()];
57  std::copy(vx.begin(), vx.end(), x);
58  std::copy(vy.begin(), vy.end(), y);
59  TGraph *gr = new TGraph(vx.size(), x, y);
60  return gr;
61 }
62 
63 std::vector<TGraph*> MoveDeltaToT2KUnits (TGraph* gr){
64  std::vector<double> x_less, y_less, x_more, y_more;
65  std::vector<double> x_less1, y_less1, x_less2, y_less2;
66  std::vector<TGraph*> new_gr;
67  for (int i = 0 ; i< gr->GetN(); i++){
68  double x; double y;
69  gr->GetPoint(i, x , y);
70  if (x>=0.85) {x_more.push_back(x-2); y_more.push_back(y);}
71  if (x<=1.15) {x_less.push_back(x); y_less.push_back(y);}
72  }
73  if (x_more.size()>0) {
74  new_gr.push_back(VectorsToGraph(x_more, y_more));
75  if (x_less.size()>0){
76  double max = *max_element(x_less.begin(), x_less.end());
77  int n = std::distance(x_less.begin(), std::find(x_less.begin(), x_less.end(), max));
78  for (int i = 0; i<x_less.size(); i++){
79  if (i <= n ) {x_less1.push_back(x_less[i]); y_less1.push_back(y_less[i]);}
80  else {x_less2.push_back(x_less[i]); y_less2.push_back(y_less[i]);}
81  }
82  if(x_less1.size()!=0 ) new_gr.push_back(VectorsToGraph(x_less1, y_less1));
83  if(x_less2.size()!=0 ) new_gr.push_back(VectorsToGraph(x_less2, y_less2));
84  }
85  }
86  else if (x_less.size()>0){
87  new_gr.push_back(VectorsToGraph(x_less, y_less));
88  }
89  return new_gr;
90 }
91 
92 TGraph* RotateT2K(TGraph* gr){
93  double x[gr->GetN()], y[gr->GetN()];
94  for (int i = 0 ; i< gr->GetN(); i++){
95  double x1; double y1;
96  gr->GetPoint(i, x1 , y1);
97  x[i] = x1; y[i] = y1/M_PI;
98  }
99  TGraph *gr_new = new TGraph(gr->GetN(), y, x);
100  return gr_new;
101 }
102 
103 TGraph* MoveT2KtoNOvAUnits(TGraph* gr){
104  double x[gr->GetN()], y[gr->GetN()];
105  for (int i = 0 ; i< gr->GetN(); i++){
106  double x1; double y1;
107  gr->GetPoint(i, x1 , y1);
108  x[i] = x1+2; y[i] = y1;
109  }
110  TGraph *gr_new = new TGraph(gr->GetN(), x, y);
111  return gr_new;
112 }
113 
114 void nova_t2k_comparison(bool contour = true)
115 {
116 
117  // Loading NOvA contours
118  TFile * infile = new TFile ("NOvA_2020_official_contour_deltassth23.root","read");
119  DrawContourCanvas("delta_th23", 0., 2., 0.275, 0.725);
120  std::vector<TGraph*> saveallgr;
121  for (auto &str:{"1sigma", "90CL"}){
122  for ( int k = 0; k < 2; k++ ){
123  TGraph * gr = (TGraph*)infile->Get((TString)"NH_"+str+"_"+std::to_string(k));
124  gr->Draw("l same");
125  saveallgr.push_back(gr);
126  }
127  }
128  TMarker *novabf = (TMarker*)infile->Get("NH_best_fit");
129  gPad->Print( "test_nova_contours.pdf" );
130  // Load T2K contours
131  TFile * t2k = new TFile ("t2k_th23_deltacp_wrc_fits.root","read");
132  TGraph * t2k_tmp68 = (TGraph*)t2k->Get("g68_nh_0");
133  TGraph * t2k_tmp90 = (TGraph*)t2k->Get("g90_nh_0");
134  TMarker *t2kbf_tmp = (TMarker*)t2k->Get("best_fit_nh");
135  TMarker *t2kbf = new TMarker (t2kbf_tmp->GetY()/M_PI, t2kbf_tmp->GetX(), 21);
136  TMarker *t2kbf_novau = new TMarker (t2kbf_tmp->GetY()/M_PI + 2, t2kbf_tmp->GetX(), 21);
137  new TCanvas;
138  t2k_tmp90->Draw();
139  t2k_tmp68->Draw("same");
140  gPad->Print( "test_t2k_contours.pdf" );
141 
142  // Work with input curves to make them in the same format
143  auto gr_t2k68 = RotateT2K(t2k_tmp68); // in T2K files delta is Y axis and th23 is X axis + make delta in pi units
144  auto gr_t2k90 = RotateT2K(t2k_tmp90);
145  // NOvA has two 1sigma graphs
146  auto vec68_1 = MoveDeltaToT2KUnits(saveallgr[0]); // move [0, 2pi] to [-pi, pi]: slice graphs that cross pi etc
147  auto vec68_2 = MoveDeltaToT2KUnits(saveallgr[1]);
148  auto vec90_1 = MoveDeltaToT2KUnits(saveallgr[2]); // move [0, 2pi] to [-pi, pi]: slice graphs that cross pi etc
149  auto vec90_2 = MoveDeltaToT2KUnits(saveallgr[3]);
150 
151  auto gr_t2k68_novau = MoveT2KtoNOvAUnits(gr_t2k68);
152  auto gr_t2k90_novau = MoveT2KtoNOvAUnits(gr_t2k90);
153 
154  t2kbf->SetMarkerStyle(21);
155  //t2kbf->SetMarkerColor(kBlue+3);
156  novabf->SetMarkerStyle(43);
157  novabf->SetMarkerColor(kBlue+3);
158  novabf->SetMarkerSize(2);
159  t2kbf_novau->SetMarkerStyle(21);
160 
161  // Draw 68% CL
162  DrawContourCanvas("delta_th23", -1., 1., 0.275, 0.725, true);
163  for (auto &gr: {vec68_1[0], vec68_1[1], vec68_1[2], vec68_2[0]}){
164  gr->SetLineColor(cnh);
165  gr->SetLineWidth(3);
166  }
167  vec68_1[0]->Draw("l same");
168  vec68_1[1]->Draw("l same");
169  vec68_1[2]->Draw("l same");
170  vec68_2[0]->Draw("l same");
171  //gr_t2k68->SetLineColor(kBlue+3);
172  gr_t2k68->SetLineStyle(2);
173  gr_t2k68->SetLineWidth(3);
174  gr_t2k68->Draw("l same");
175  t2kbf->Draw("same");
176  novabf->Draw("same");
177  TText * txtNH = new TText(0.4,0.85,"Normal Hierarchy, 68% CL");
178  txtNH->SetNDC();
179  txtNH->SetTextAlign(22);
180  txtNH->SetTextSize(kBlessedLabelFontSize);
181  txtNH->Draw();
182  Preliminary();
183  TLegend *leg = new TLegend(0.18,0.18,0.5,0.35);
184  leg->SetTextSize(kBlessedLabelFontSize);
185  leg->AddEntry(gr_t2k68, "T2K 2020", "l");
186  leg->AddEntry(vec68_1[0], "NOvA", "l");
187  leg->Draw();
188  gPad->Print( "nova-t2k-1sigma.pdf" );
189 
190  // Draw 90% CL
191  DrawContourCanvas("delta_th23", -1., 1., 0.275, 0.725, true);
192  for (auto &gr: {vec90_1[0], vec90_1[1], vec90_2[1], vec90_2[0], vec90_2[2]}){
193  gr->SetLineColor(cnh);
194  gr->SetLineWidth(3);
195  gr->Draw("l same");
196  }
197  //gr_t2k90->SetLineColor(kBlue+3);
198  gr_t2k90->SetLineWidth(3);
199  gr_t2k90->Draw("l same");
200  t2kbf->Draw("same");
201  novabf->Draw("same");
202  TText * txtNH2 = new TText(0.4,0.85,"Normal Hierarchy, 90% CL");
203  txtNH2->SetNDC();
204  txtNH2->SetTextAlign(22);
205  txtNH2->SetTextSize(kBlessedLabelFontSize);
206  txtNH2->Draw();
207  Preliminary();
208  TLegend *leg2 = new TLegend(0.18,0.18,0.5,0.35);
209  leg2->SetTextSize(kBlessedLabelFontSize);
210  leg2->AddEntry(gr_t2k90, "T2K 2020", "l");
211  leg2->AddEntry(vec90_1[0], "NOvA", "l");
212  leg2->Draw();
213  gPad->Print( "nova-t2k-90cl.pdf" );
214 
215  //All together
216  TGraph* nova_90_1 = JointGraphs(vec90_1[0], vec90_1[1]);
217  TGraph* nova_90_2 = JointGraphs(vec90_2[2], vec90_2[0]);
218  TGraph* nova_90 = JointGraphs(nova_90_1, nova_90_2);
219  TGraph* nova_68_1_tmp1 = JointGraphs(vec68_2[0], vec68_1[1]);
220  TGraph* nova_68_1 = JointGraphs(vec68_1[2], nova_68_1_tmp1);
221  Int_t kFillColor90 = TColor::GetColorTransparent(cnh, 0.25);
222  Int_t kFillColor68 = TColor::GetColorTransparent(cnh, 0.9);
223  Int_t kLineColor90 = TColor::GetColorTransparent(cnh, 0.7);
224  DrawContourCanvas("delta_th23", -1., 1., 0.275, 0.725, true);
225  nova_90->SetFillColor(kFillColor90);
226  nova_90->SetLineColor(kLineColor90);
227  nova_90->Draw("f same");
228  for (auto &gr: {nova_68_1, vec68_1[0]}) {
229  gr->SetFillColor(kWhite);
230  gr->Draw("f same");
231  gr->SetFillColor(kFillColor68);
232  gr->SetLineColor(cnh_dark);
233  gr->Draw("f same");
234  gr->Draw("l same");
235  }
236  nova_90->Draw("l same");
237  gr_t2k90->Draw("l same");
238  gr_t2k68->Draw("l same");
239  t2kbf->Draw("same");
240  novabf->Draw("same");
241  TText * txtNH3 = new TText(0.35,0.85,"Normal Hierarchy");
242  txtNH3->SetNDC();
243  txtNH3->SetTextAlign(22);
244  txtNH3->SetTextSize(kBlessedLabelFontSize*0.9);
245  txtNH3->Draw();
246  Preliminary();
247  TLegend *leg3 = new TLegend(0.13,0.25,0.85,0.35);
248  leg3->SetFillStyle(0);
249  leg3->SetNColumns(4);
250  leg3->SetTextSize(kBlessedLabelFontSize*0.7);
251  TMarker* tmp = new TMarker(10, 1, 1);
252  tmp->SetMarkerColor(kWhite);
253  leg3->AddEntry(tmp, "T2K, Nature 580:", "p");
254  leg3->AddEntry(t2kbf, "BF ", "p");
255  leg3->AddEntry(gr_t2k90, " #leq 90% CL ", "l");
256  leg3->AddEntry(gr_t2k68, " #leq 68% CL", "l");
257  leg3->Draw();
258  TLegend *leg4 = new TLegend(0.265,0.18,0.87,0.25);
259  leg4->SetFillStyle(0);
260  leg4->SetNColumns(4);
261  leg4->SetTextSize(kBlessedLabelFontSize*0.7);
262  leg4->AddEntry(tmp, "NOvA:", "p");
263  leg4->AddEntry(novabf, " BF ", "p");
264  leg4->AddEntry(nova_90, " #leq 90% CL ", "f");
265  leg4->AddEntry(nova_68_1, " #leq 68% CL", "f");
266  leg4->Draw();
267  gPad->Update();
268  gPad->RedrawAxis();
269  gPad->Print( "nova-t2k-68+90.pdf" );
270 
271  // 1d DeltaCP slice
272  TFile * infile_slice = new TFile ("NOvA_2020_official_slice_delta.root","read");
273  TGraph * nova = (TGraph *)infile_slice->Get("NH UO");
274  auto vecnova = MoveDeltaToT2KUnits(nova);
275  DrawSliceCanvas("delta", 5, -1, 1, false, true);
276  vecnova[0]->Draw("l same");
277  vecnova[1]->Draw("l same");
278  gPad->Print( "test_slice.pdf" );
279  TFile * t2k_slice = new TFile ("t2k_deltacp_wrc_fits.root","read");
280  TGraph * gr_t2k = (TGraph *)t2k_slice->Get("dchisq_map_nh");
281  new TCanvas;
282  gr_t2k->Draw();
283  gPad->Print( "input_t2k_slice.pdf" );
284 
285  DrawContourCanvas("delta_th23", 0., 2., 0.275, 0.725);
286  TGraph* nova_90_usual = JointGraphs(saveallgr[2], saveallgr[3]);
287  nova_90_usual->SetFillColor(kFillColor90);
288  nova_90_usual->SetLineColor(kLineColor90);
289  nova_90_usual->SetLineWidth(3);
290  nova_90_usual->Draw("f same");
291  for (auto &gr: {saveallgr[0], saveallgr[1]}) {
292  gr->SetFillColor(kWhite);
293  gr->Draw("f same");
294  gr->SetLineWidth(3);
295  gr->SetFillColor(kFillColor68);
296  gr->SetLineColor(cnh_dark);
297  gr->Draw("f same");
298  gr->Draw("l same");
299  }
300  nova_90_usual->Draw("l same");
301  novabf->Draw("same");
302  gr_t2k68_novau->SetLineStyle(2);
303  gr_t2k68_novau->SetLineWidth(3);
304  gr_t2k90_novau->SetLineWidth(3);
305  gr_t2k68_novau->Draw("l same");
306  gr_t2k90_novau->Draw("l same");
307  t2kbf_novau->Draw("same");
308  txtNH3->Draw();
309  Preliminary();
310  leg3->Draw();
311  leg4->Draw();
312  gPad->Update();
313  gPad->RedrawAxis();
314  gPad->Print( "nova-t2k-68+90_nova_units.pdf" );
315 
316 
317 }
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2.)
TGraph * MoveT2KtoNOvAUnits(TGraph *gr)
void nova_t2k_comparison(bool contour=true)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
Float_t y1[n_points_granero]
Definition: compare.C:5
TGraph * RotateT2K(TGraph *gr)
Float_t x1[n_points_granero]
Definition: compare.C:5
std::vector< TGraph * > MoveDeltaToT2KUnits(TGraph *gr)
Float_t tmp
Definition: plot.C:36
unsigned distance(const T &t1, const T &t2)
TGraphAsymmErrors * t2k
Definition: Xsec_final.C:30
#define M_PI
Definition: SbMath.h:34
TGraph * VectorsToGraph(std::vector< double > vx, std::vector< double > vy)
string infile
const double a
Hold drift constants from current run.
Definition: DriftCache.h:17
void Preliminary()
void DrawContourCanvas(TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1)
const hit & b
Definition: hits.cxx:21
TGraph * JointGraphs(TGraph *a, TGraph *b)
const Float_t kBlessedLabelFontSize
Definition: Style.h:90
std::string to_string(ModuleType mt)
Definition: ModuleType.h:32
T max(sqlite3 *const db, std::string const &table_name, std::string const &column_name)
Definition: statistics.h:68
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29