plot_cont_expfriends.C
Go to the documentation of this file.
1 /////////////////////////////////////////////////////////////02
2 // //
3 // plot_cont_expfriends.C //
4 // picks up output root files from make_surfprof.C //
5 // and results from our friendds MINOS, IceCube, T2K and SK //
6 // to plot all contours overlayed //
7 // //
8 ///////////////////////////////////////////////////////////////
9 
12 #include "CAFAna/FC/FCSurface.h"
14 
15 std::string NoFC = "No Feldman-Cousins";
16 
18 
19 
20  bool preliminary = false;
21  bool lazyprint = true;
22  bool numuonly = true;
23  bool onlycf = true;
24  // print only confidence levels so it
25  // does or doesnt change contour colors
26  // already define colors
27  Color_t color68joint = kBlue-7;
28  Color_t color90joint = kViolet+2;
29  Color_t color99joint = kMagenta+1;
30  Color_t colorjoint = kBlack;
31  Color_t colornumu = kBlack;
32  Color_t colorsa = kGreen+1;
33  Color_t coloric = kBlue;
34  Color_t colorsk = kMagenta+1;
35  Color_t colort2k = kGreen+2;
36  Color_t colorminos = kRed;
37  Color_t colorsens = kGray+1;
38  gStyle->SetLineStyleString(11,"400 200");
39 
40  Style_t line68joint = kSolid;
41  Style_t line90joint = kSolid;
42  Style_t line99joint = kSolid;
43  Style_t linejoint = kSolid;
44  Style_t linenumu = kSolid;
45  Style_t linesa = 10;
46  Style_t lineic = kDotted; // dotted
47  Style_t linesk = 6; // dashed
48  Style_t linet2k = 7; // dashed
49  Style_t lineminos = 9; // dotted
50  Style_t linesens = kDashed;
51 
52  Style_t marker68joint = kFullCircle;
53  Style_t marker90joint = kFullCircle;
54  Style_t marker99joint = kFullCircle;
55  Style_t markerjoint = kFullCircle;
56  Style_t markernumu = kFullCircle;
57  Style_t markersa = kFullCircle;
58  Style_t markersens = kFullCircle;
59 
60 
61  // nova 2018 joint analysis
62  // Get Gaussian Stuff
63  TString hier = "NH";
64  TString gaussDir = "/nova/ana/nu_e_ana/Ana2018/Results/RHC_and_FHC/contours/th23_dmsq/syst/";
65  TString gaussName = "hist_contours_2018_joint_realData_both_only";
66  gaussName += hier;
67  gaussName += "combo_dmsq_systs.root";
68 
69  TFile* fGauss = new TFile(gaussDir + gaussName);
70  TH1* hGauss = (TH1*) fGauss->Get("th23_dm32_"+hier);
71  auto surfJoint = *ISurface::LoadFrom(fGauss, "surf_th23_dm32_" + hier);
72  auto mins =* (TVectorD*)fGauss->Get("overall_minchi");
73  double minchi23 = mins[0];
74 
75  int NX = hGauss->GetXaxis()->GetNbins();
76  int NY = hGauss->GetYaxis()->GetNbins();
77  double minX = hGauss->GetXaxis()->GetBinCenter(1);
78  double maxX = hGauss->GetXaxis()->GetBinCenter(NX);
79  double minY = hGauss->GetYaxis()->GetBinCenter(1);
80  double maxY = hGauss->GetYaxis()->GetBinCenter(NY);
81 
82  std::cout << "Making FCSurface with " << std::endl;
83  std::cout << "------------------------------------" << std::endl;
84  std::cout << "NX: " << NX << "\tNY: " << NY << std::endl;
85  std::cout << "minX: " << minX << "\tmaxX" << maxX << std::endl;
86  std::cout << "minY: " << minY << "\tmaxY" << maxY << std::endl;
87 
88  // Get FC Stuff
89  TString fcFolder = "/nova/ana/nu_e_ana/Ana2018/FC/";
90  TString fcFileName = "ssth23dmsq32_";
91  fcFileName += hier.Contains("NH")? "nh":"ih";
92  fcFileName += ".root";
93 
94  FCSurface *fcXH = 0;
95  auto fccol = FCCollection::LoadFromFile((fcFolder+fcFileName).Data()).release();
96  fcXH = new FCSurface(NX, minX, maxX, NY, minY, maxY);
97 
98  fcXH->Add(*fccol, "ssth23dmsq32");
99  TH2* fc_surf_90CL = fcXH->SurfaceForSignificance(0.90);
100  fc_surf_90CL->Smooth();
101  auto graphJoint = surfJoint.GetGraphs(fc_surf_90CL,minchi23);
102 
103 
104 
105  std::cout << "\n\n============================================" << std::endl;
106  std::cout << "============================================\n" << std::endl;
107  std::cout << "Now get your friends and hug them together!" << std::endl;
108  std::cout << "\n============================================" << std::endl;
109  std::cout << "============================================\n\n" << std::endl;
110  // get icecube results
111  TString fIC = FindCAFAnaDir().c_str() + (TString)"/data/expt/IceCubeDeepCore_2017/IC2017_90CL_FC__copy.dat";
112  std::ifstream inIC; inIC.open(fIC);
113  const int pointsIC = 81;
114  double valIC[pointsIC][2];
115  for(int row = 0; row < pointsIC; row ++){ // pair
116  for(int col = 0; col < 2; col ++){ // (sin,delta)
117  inIC >> valIC[row][col];
118  }
119  }
120  inIC.close();
121  TGraph *icplot = new TGraph;
122  for (int point = 0; point < pointsIC; point++) {
123  double xval, yval;
124  icplot->SetPoint(point,valIC[point][0],valIC[point][1]);
125  }
126  icplot->SetLineStyle(lineic);
127  icplot->SetLineColor(coloric);
128  icplot->SetLineWidth(2);
129 
130 
131  // get and rescale t2k results
132  TFile* fT2K = new TFile(FindCAFAnaDir().c_str() + (TString)"/data/expt/T2KJointNuNuBarOscillation_Run18_2017/sinsqth23Vsdmsq_Run18_Data_react_nh.root");
133  TGraph *t2kplot = (TGraph*)fT2K->Get("g90_0");
134  fT2K -> Close();
135  for (int point = 0; point < t2kplot->GetN(); point++) {
136  double xval, yval;
137  t2kplot->GetPoint(point,xval,yval);
138  yval = yval * 1000.0;
139  t2kplot->SetPoint(point,xval,yval);
140  }
141  t2kplot->SetLineStyle(linet2k);
142  t2kplot->SetLineColor(colort2k);
143  t2kplot->SetLineWidth(2);
144 
145 
146  // get and rescale superk results
147  TFile* fSK = new TFile(FindCAFAnaDir().c_str() + (TString)"/data/expt/SuperK_2017/sk1234.atmospheric.201803.nh.root");
148  TGraph *skplot = (TGraph*)fSK->Get("sk1234_NH_S23M23_90CL");
149  fSK -> Close();
150  for (int point = 0; point < skplot->GetN(); point++) {
151  double xval, yval;
152  skplot->GetPoint(point,xval,yval);
153  yval = yval * 1000.0;
154  skplot->SetPoint(point,xval,yval);
155  }
156  skplot->SetLineStyle(linesk);
157  skplot->SetLineColor(colorsk);
158  skplot->SetLineWidth(2);
159 
160 
161 
162  // get and rescale minos result
163  TFile* minosfile = new TFile(FindCAFAnaDir().c_str() + (TString)"/data/expt/MINOSx_Nu2018_dmsq32_sinsq23_contours.root");
164  TGraph *minosplot = (TGraph*)minosfile->Get("NH_90percent");
165  minosfile->Close();
166  for (int point = 0; point < minosplot->GetN(); point++) {
167  double xval, yval;
168  minosplot->GetPoint(point,xval,yval);
169  yval = yval * 1000.0;
170  minosplot->SetPoint(point,xval,yval);
171  }
172  minosplot->SetLineStyle(lineminos);
173  minosplot->SetLineColor(colorminos);
174  minosplot->SetLineWidth(2);
175 
176 
177 
178  if(lazyprint){
179  numuonly = false;
180  // canvas
181  // 90% systs + T2K + MINOS
182  std::ofstream file;
183  new TCanvas("canvas","canvas",1120,800);
184  TString out = "../plots_blessing/contours_nova_rhcfhc__withfriends__numu2018";
185  gPad->SetFillStyle(0);
186  gPad->SetTopMargin(0.08);
187  gPad->SetRightMargin(0.08);
188  gPad->SetBottomMargin(0.12);
189  gPad->SetLeftMargin(0.12);
190  gStyle->SetTitleOffset(0.95,"x"); //override style script so nothing is clipped or overlapping with axis labels
191  gStyle->SetTitleOffset(0.85,"y");
192  TH2F *axes = new TH2F("setaxes","" ,38 ,0.32001, 0.6999, 55, 1.90, 3.3);
193  axes->GetXaxis()->SetTitle("sin^{2}#theta_{23}");
194  axes->GetYaxis()->SetTitle("#Deltam^{2}_{32} (10^{-3} eV^{2})");
195  axes->GetXaxis()->SetTitleOffset(0.95);
196  axes->GetYaxis()->SetTitleOffset(0.85);
197  axes->GetXaxis()->SetTitleSize(0.055);
198  axes->GetYaxis()->SetTitleSize(0.055);
199  axes->GetXaxis()->SetLabelSize(0.04);
200  axes->GetYaxis()->SetLabelSize(0.04);
201  axes->GetXaxis()->CenterTitle();
202  axes->GetYaxis()->CenterTitle();
203  axes->SetTitle("");
204  axes->GetYaxis()->SetDecimals();
205  axes->Draw();
206  axes->Draw();
207  icplot->Draw("same");
208  skplot -> Draw("same");
209  t2kplot -> Draw("same");
210  minosplot-> Draw("same");
211  DrawLegendBFNull(surfJoint.GetMinX(), surfJoint.GetMinY(), kBlack, kFullCircle);
212  surfJoint.DrawBestFit(kBlack, kFullCircle);
213  graphJoint[0]->SetLineWidth(3);
214  graphJoint[0]->Draw("c");
216  // for nova only
217  TPaveText *pText1 = new TPaveText(0.23, 0.84, 0.43, 0.89, "brNDC");
218  TText *text1 = pText1->AddText("Normal Hierarchy 90% CL");
219  text1->SetTextSize(0.05);
220  pText1->SetBorderSize(0);
221  pText1->SetFillStyle(0);
222  pText1->Draw();
223  auto leg = new TLegend(0.15,0.67,0.70,0.82);
224  leg->SetNColumns(2);
225  leg->SetFillStyle(0);
226  leg->SetTextSize(0.05);
227  leg->SetMargin(1.3*leg->GetMargin());
228  TGraph * dummy = new TGraph();
229  dummy->SetLineColor(colorjoint);
230  dummy->SetLineStyle(linejoint);
231  dummy->SetLineWidth(3);
232  leg->AddEntry(dummy,"NOvA","l");
233  leg->AddEntry(minosplot,"MINOS+ 2018","l");
234  leg->AddEntry(t2kplot,"T2K 2018","l");
235  leg->AddEntry(icplot,"IceCube 2017","l");
236  leg->AddEntry(skplot,"SK 2017","l");
237  leg->Draw();
238  gPad->Modified();
239  gPad->Print(out + ".png");
240  gPad->Print(out + ".pdf");
241  gPad->Print(out + ".eps");
242  file.open (out + ".txt");
243 
244  file<< "90% confidence level contour from NOvA's 2018 joint fit with neutrino + antineutrino data with exposure to 8.85e20+6.89e20 POT-equivalent, with systematics and Feldman-Cousins corrections applied. The latest results from MINOS+(Neutrino 2018), T2K(Phys.Rev.Lett. 121, 171802 (2018), SuperKamiokande(Phys. Rev. D 97, 072001 (2018)) and Ice-Cube(Phys. Rev. Lett. 120, 071801 (2018)) are plotted for reference.";
245  file.close();
246 
247  }
248 
249 
250 
251 } // End of function
252 
tree Draw("slc.nhit")
enum BeamMode kRed
double minY
Definition: plot_hist.C:9
double maxY
Definition: plot_hist.C:10
std::string FindCAFAnaDir()
Definition: Utilities.cxx:204
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
Int_t col[ntarg]
Definition: Style.C:29
void plot_cont_expfriends()
Int_t preliminary
Definition: SimpleIterate.C:63
std::string NoFC
void DrawLegendBFNull(double bfSin, double bfDm, Color_t color, Style_t style)
Definition: numu_tools.h:80
void PreliminaryBoxOpening()
Definition: numu_tools.h:187
void Add(const FCPoint &pt, std::string plot)
Definition: FCSurface.cxx:39
c1 Close()
OStream cout
Definition: OStream.cxx:6
enum BeamMode kViolet
TFile * file
Definition: cellShifts.C:17
TH2 * SurfaceForSignificance(double sig)
Definition: FCSurface.cxx:103
enum BeamMode kGreen
enum BeamMode kBlue
FileContents LoadFromFile(const std::string &filename)
Pseudo-experiments, binned to match a log-likelihood surface.
Definition: FCSurface.h:14
enum BeamMode string