nova_official_data_release.C
Go to the documentation of this file.
1 #include "TGraph.h"
2 #include "TFile.h"
3 #include "TH2.h"
4 #include "TVector.h"
5 #include "TCanvas.h"
6 #include "TLegendEntry.h"
7 
8 #include "CAFAna/FC/FCSurface.h"
12 #include "Utilities/rootlogon.C"
13 
14 
15 #include <fstream>
16 
17 using namespace ana;
18 
19 //copied from CAFAna/nue/Ana2018/nue/FitandFC/joint_fit_2018_slices.C
20 TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName,
21  bool fourierFit, std::string plot_name,
22  ofstream& summary, TString par, TString hier, TString oct,
23  bool useRoot = false, bool fccol=false);
24 
25 TGraph * SmoothWithFourierFit(TGraph* gr0,
26  bool isDelta = true,
27  bool useRoot = false,
28  int fourierFitN = 2);
29 
30 void SmoothWithFourierFit(TH1* hist,
31  int nbins = 80, double minX = 0, double maxX = 2,
32  int fourierFitN=3);
33 
34 
35 TGraph* HistoToCyclicGraph(TH1* h, bool useRoot);
36 
37 void AverageCyclicPoints(TGraph* g);
38 
39 void AddCyclicPoints(TGraph* g);
40 // end from joint_fit_2018_slices.C
41 
42 TH1D* ProfileAxes(TString par);
43 
44 TFile* GetGaussianProfile(TString par)
45 {
46  TString gaussDir = "/nova/ana/nu_e_ana/Ana2018/Results/RHC_and_FHC/slices/syst/";
47  TString baseName = "hist_slices_2017_joint_realData_bothcombo_systs";
48  TString shortName;
49  if(par.Contains("ssth23"))
50  shortName = "_th23_noOct";
51  else if(par.Contains("dmsq"))
52  shortName = "_dmsq";
53  else
54  shortName = "";
55  TFile* ret = TFile::Open(gaussDir + baseName + shortName + ".root");
56  return ret;
57 }
58 
59 TString GetFCPath(TString par, TString hier)
60 {
61  TString fcFolder = "/nova/ana/nu_e_ana/Ana2018/FC/";
62  TString fcFileName = par + "_";
63  fcFileName += hier.Contains("NH")? "nh":"ih";
64  fcFileName += ".root";
65  return fcFolder + fcFileName;
66 }
67 
68 TString GetFCPath(TString par, TString hier, TString oct)
69 {
70  TString fcFolder = "/nova/ana/nu_e_ana/Ana2018/FC/";
71  hier.ToLower();
72  if(par.Contains("ssth23"))
73  return fcFolder + TString::Format("%s_%s.root", par.Data(), hier.Data());
74  else
75  return fcFolder + TString::Format("%s_%s_%s.root", par.Data(), oct.Data(), hier.Data());
76 }
77 
78 TString MakeSummaryTable(FCSurface* fsurf, TTree* tfc,
79  TString par, TString hier, TString oct = "");
80 
81 void nova_official_data_release(bool contours = true, bool slices = true)
82 {
83  ofstream summary;
84  summary.open("fc_summary_2018/summary_table.tex");
85 
86  TString outdir = "NOvA_2018_official_data_release/";
87  system(TString::Format("mkdir %s", outdir.Data()).Data());
88  // Let's do the contours first
89  if(contours) {
90  double bf_delta = 0.16576;
91  double bf_th23 = 0.584509;
92  double bf_dmsq32 = 2.50515;
93  for(TString par : {"deltassth23", "ssth23dmsq32"}) {
94  TMarker* bf = new TMarker(bf_delta, bf_th23, 43);
95  if(par.Contains("ssth23dmsq32")) {
96  bf->SetX(bf_th23);
97  bf->SetY(bf_dmsq32);
98  }
99  TFile* official = new TFile(outdir + TString::Format("NOvA_2018_official_contour_%s.root", par.Data()), "recreate");
100  bf->Write("NH_best_fit");
101 
102  for(TString hier: {"NH", "IH"}) {
103 
104  TFile* fGauss = GetGaussianSurface(par, hier);
105  TFile* ffc = TFile::Open(GetFCPath(par, hier));
106  TTree* tfc = (TTree*)ffc->Get("fc");
107  TString rootDir = par.Contains("ssth23dmsq32")? "th23_dm32_" : "delta_th23_";
108  auto surf = ISurface::LoadFrom(fGauss, "surf_" + rootDir + hier);
109  TH2* hGauss = (TH2*) fGauss->Get(rootDir+hier);
110  auto mins =* (TVectorD*) fGauss->Get("overall_minchi");
111  double minchi23 = mins[0];
112 
113  int NX = hGauss->GetXaxis()->GetNbins();
114  int NY = hGauss->GetYaxis()->GetNbins();
115  double minX = hGauss->GetXaxis()->GetBinLowEdge(1);
116  double maxX = hGauss->GetXaxis()->GetBinUpEdge(NX);
117  double minY = hGauss->GetYaxis()->GetBinLowEdge(1);
118  double maxY = hGauss->GetYaxis()->GetBinUpEdge(NY);
119 
120  std::cout << "Making FCSurface with " << std::endl;
121  std::cout << "------------------------------------" << std::endl;
122  std::cout << "NX: " << NX << "\tNY: " << NY << std::endl;
123  std::cout << "minX: " << minX << "\tmaxX" << maxX << std::endl;
124  std::cout << "minY: " << minY << "\tmaxY" << maxY << std::endl;
125 
126  FCSurface* fc = new FCSurface(NX, minX, maxX, NY, minY, maxY);
127 
128  auto fccol = FCCollection::LoadFromFile(GetFCPath(par, hier).Data()).release();
129  fc->Add(*fccol, par.Data());
130 
131  // output TeX formatted summary table
132  summary << MakeSummaryTable(fc, tfc , par, hier) << std::endl;
133 
134 
135  TH2* fc_surf_1sig = fc->SurfaceForSignificance(0.6827); fc_surf_1sig->Smooth();
136  TH2* fc_surf_2sig = fc->SurfaceForSignificance(0.9545); fc_surf_2sig->Smooth();
137  TH2* fc_surf_3sig = fc->SurfaceForSignificance(0.9973);
138  TH2* fc_surf_90CL = fc->SurfaceForSignificance(0.90); fc_surf_90CL->Smooth();
139  if(par.Contains("deltassth23"))
140  CylindricalSmoothing(fc_surf_3sig);
141  else
142  fc_surf_3sig->Smooth();
143 
144  std::vector<TGraph*> g1 = surf->GetGraphs(fc_surf_1sig, minchi23);
145  std::vector<TGraph*> g2 = surf->GetGraphs(fc_surf_2sig, minchi23);
146  std::vector<TGraph*> g3 = surf->GetGraphs(fc_surf_3sig, minchi23);
147  std::vector<TGraph*> g90 = surf->GetGraphs(fc_surf_90CL, minchi23);
148 
149  // Start writing stuff out
150  official->cd();
151 
152  // canvas axes are different than what is returned from the surface object
153  // as a style choice. Values are pulled from
154  // CAFAna/nue/Ana2018/FitandFC/joint_fit_2018_contours.C
155  double axesMaxX; double axesMaxY;
156  double axesMinX; double axesMinY;
157  if(par.Contains("deltassth23")){
158  axesMinX = 0.;
159  axesMaxX = 2.;
160  axesMinY = 0.275;
161  axesMaxY = 0.725;
162  }
163  else {
164  axesMinX = 0.3;
165  axesMaxX = 0.7;
166  if(hier.Contains("NH")) {
167  axesMinY = 2.05;
168  axesMaxY = 2.85;
169  }
170  else {
171  axesMinY = -2.95;
172  axesMaxY = -2.25;
173  }
174  }
175 
176  TH2F* axes = new TH2F("axes", "", 100, axesMinX, axesMaxX, 100, axesMinY, axesMaxY);
177  TString xLabel = "#delta_{CP}/#pi";
178  TString yLabel = "sin^{2}#theta_{23}";
179  if(par.Contains("ssth23dmsq32")){
180  xLabel = "sin^{2}#theta_{23}";
181  yLabel = "#Deltam^{2}_{32} (10^{-3} eV^{2})";
182  }
183  axes->GetXaxis()->SetTitle(xLabel.Data());
184  axes->GetYaxis()->SetTitle(yLabel.Data());
185  CenterTitles(axes);
186 
187  // For consistency, only save one axis for deltassth23 as just 'axes'
188  TString axesName = hier + "_axes";
189  if(par.Contains("deltassth23"))
190  axesName = "axes";
191  if(!(par.Contains("deltassth23") && hier.Contains("IH"))){
192  axes->Write(axesName.Data());
193  }
194  std::vector<std::pair< TString, std::vector<TGraph*> > > graphs = {{"1sigma",g1},
195  {"2sigma", g2},
196  {"3sigma", g3},
197  {"90CL", g90}};
198  for(std::pair<TString, std::vector<TGraph*>> gn : graphs){
199  int i = 0;
200  for(TGraph* g : gn.second){
201  g->Write(TString::Format("%s_%s_%d",
202  hier.Data(),
203  gn.first.Data(),
204  i));
205  i++;
206  }
207  }
208  }
209  }
210  }
211 
212  // 1D profiles
213  if(slices){
214  std::map< TString, double > maxY = {{"delta", 5.0},
215  {"ssth23", 3.0},
216  {"dmsq", 3.0}};
217  for(TString par : {"delta", "ssth23", "dmsq"}) {
218  TString parShortName = par;
219  if(par.Contains("ssth23")) parShortName = "th23";
220  TFile* official = new TFile(outdir + TString::Format("NOvA_2018_official_profile_%s.root", par.Data()), "recreate");
221  TH1D* axes = ProfileAxes(par);
222  axes->SetMaximum(maxY[par]);
223  axes->Write("axes");
224 
225  TFile* fGauss = GetGaussianProfile(par);
226  // this will hold the profiles of each parameter
227 
228 
229  for(TString hier : {"IH", "NH"}) {
230  for(TString oct : {"LO", "UO"}) {
231  if(par.Contains("ssth23") && oct.Contains("LO"))
232  continue; // only do ssth23 once for each hierarchy
233 
234  if(par.Contains("ssth23"))
235  oct = "";
236  TString histName = "slice_" + parShortName + "_" + hier + oct;
237  std::cout << "Looking for " << histName << " in " << fGauss->GetName() << std::endl;
238  TH1D* h = (TH1D*) fGauss->Get(histName);
239  TGraph* graph = new TGraph(h);
240 
241  // only use fourier fit for period boundary conditions
242  // Fourier fit doesn't like delta_UO_nh. Would need
243  // lots of harmonics to fit it correctly anyway
244  bool useFourier = (par.Contains("delta") &&
245  !(hier.Contains("NH") && oct.Contains("UO")));
246 
247  bool smooth = true;
248  graph = FCCorrectSlice(graph, GetFCPath(par, hier, oct),
249  smooth, par.Data(),
250  summary, par, hier, oct,
251  !useFourier, true);
252 
253  if(par.Contains("dmsq")){ //HACK
254  for (int i =0; i < graph->GetN(); ++i){
255  graph->SetPoint(i,1000*fabs(graph->GetX()[i]),graph->GetY()[i]);
256  }
257  }
258  official->cd();
259  graph->Write(hier+oct);
260  }
261  }
262  }
263  }
264 
265  summary.close();
266 }
267 
268 TString MakeSummaryTable(FCSurface* fcsurf, TTree* tfc,
269  TString par, TString hier, TString oct)
270 {
271  int maxNPts = 0;
272  int minNPts = 100000;
273  int totNPts = tfc->GetEntries();
274  int NBins = 0;
275  int empty;
276  double avgNPts;
277  for(int i = 1; i < fcsurf->Binning()->GetXaxis()->GetNbins()+1; i++){
278  for(int j = 1; j < fcsurf->Binning()->GetYaxis()->GetNbins()+1; j++){
279  if(fcsurf->GetBin(i, j)->Empty()) {
280  empty++;
281  continue;
282  }
283  int NPts = fcsurf->GetBin(i, j)->NPts();
284  maxNPts = std::max(maxNPts, NPts);
285  minNPts = std::min(minNPts, NPts);
286  NBins++;
287  }
288  }
289  avgNPts = totNPts / double(NBins);
290 
291  TString label = "";
292  if(par == "ssth23dmsq32")
293  label = "$\\sin^2\\theta_{23} vs \\Deltam^2_{32}$";
294  if(par == "deltassth23")
295  label = "$\\delta_{CP} vs \\sin^2\\theta_{23}$";
296  if(par == "delta")
297  label = "$\\delta_{CP}$";
298  if(par == "ssth23")
299  label = "$\\sin^2\\theta_{23}$";
300  if(par == "dmsq")
301  label = "$\\Deltam^2_{32}$";
302  label += " ";
303  label += oct;
304  label += hier;
305 
306  TString row(TString::Format("%s &",label.Data()));
307  row += TString::Format("%d & ", NBins);
308  row += TString::Format("%d & ", totNPts);
309  row += TString::Format("%d & ", maxNPts);
310  row += TString::Format("%d & ", minNPts);
311  row += TString::Format("%.2f \\\\ ", avgNPts);
312 
313  std::cout << row << std::endl;
314 
315  return row;
316 
317 }
318 
319 TH1D* ProfileAxes(TString par)
320 {
321  double minX = 0;
322  double maxX = 2;
323  TString xLabel = "#delta_{CP}/#pi";
324  TString yLabel = "Significance (#sigma)";
325  // min and max values of X axis are copied from
326  // CAFAna/nue/Ana2018/FitandFC/joint_fit_2018_slices.C
327  if(par.Contains("dmsq")){
328  minX = 2.15;
329  maxX = 2.75;
330  xLabel = "|#Deltam^{2}_{32}| (10^{-3} eV^{2})";
331  }
332  else if(par.Contains("th23")){
333  minX = 0.32;
334  maxX = 0.72;
335  xLabel = "sin^{2}#theta_{23}";
336  }
337 
338  TH1D* ret = new TH1D("axes","", 100, minX, maxX);
339  ret->GetXaxis()->SetTitle(xLabel);
340  ret->GetYaxis()->SetTitle(yLabel);
341  CenterTitles(ret);
342 
343  return ret;
344 }
345 
346 TGraph* FCCorrectSlice(TGraph* sqrtslice, TString fcFileName,
347  bool fourierFit, std::string plot_name,
348  ofstream& summary,
349  TString par, TString hier, TString oct,
350  bool useRoot, bool fccol)
351 {
352  bool th23plot = fcFileName.Contains("th23");
353  bool dmsqplot = fcFileName.Contains("dmsq");
354  TString hiername = fcFileName.Contains("nh")? "NH":"IH";
355  double minX = 0;
356  double maxX = 2;
357  if (th23plot) {minX = 0.3; maxX = 0.7;}
358  if (dmsqplot) {
359  minX = hiername=="NH" ? 2e-3 : -3e-3;
360  maxX = hiername=="NH" ? 3e-3 : -2e-3;
361  }
362  //We want slice, not sqrt
363  auto n= sqrtslice->GetN();
364  auto sliceXH = new TGraph(n);
365  for (int i = 0; i < n;i++ ){
366  sliceXH->SetPoint(i,sqrtslice->GetX()[i],util::sqr(sqrtslice->GetY()[i]));
367  }
368  FCSurface *fcXH = 0;
369  if(!fccol) fcXH = FCSurface::LoadFromFile(fcFileName.Data()).release();
370  else {
371  auto fccol = FCCollection::LoadFromFile(fcFileName.Data()).release();
372  //Only one Y bin required for slices.
373  // All slices this year had 60 bins. Need to find a better way
374  // next year
375  fcXH = new FCSurface(60,minX,maxX,1,0,1);
376  fcXH->Add(*fccol, plot_name);
377  }
378  TFile* ffc = TFile::Open(fcFileName.Data());
379  TTree* tfc = (TTree*) ffc->Get("fc");
380 
381  // output TeX formatted summary table
382  summary << MakeSummaryTable(fcXH, tfc , par, hier, oct) << std::endl;
383 
384  TGraph* sigXH = new TGraph;
385  TGraph* pcXH = new TGraph;
386  for(int i = 0; i < n; ++i){
387  const double delta = sliceXH->GetX()[i];
388  const double upXH = sliceXH->GetY()[i];
389  // Fall back to gaussian up if dchisq is too large or FCBin is empty.
390  // Harder to believe these dchisq are covered with our
391  // pseudoexperiments anyway
392  if(upXH > 14 || fcXH->GetBin(i+1,1)->Empty()){
393  sigXH->SetPoint(sigXH->GetN(), delta, sqrt(upXH));
394  }
395  else {
396  double pXH;
397  if(upXH > 0) pXH = fcXH->GetBin(i+1, 1)->SignificanceForUpValue(upXH);
398  else pXH = 0; //HACK HACK
399  sigXH->SetPoint(sigXH->GetN(), delta, PValueToSigma(1-pXH));
400  }
401  }
402  if(fourierFit) {
403  sigXH = SmoothWithFourierFit(sigXH, fcFileName.Contains("delta"), useRoot);
404  }
405  return sigXH;
406 }
407 
408 
409 TGraph * SmoothWithFourierFit(TGraph* gr0, bool isDelta, bool useRoot, int fourierFitN)
410 {
411  int nbins = gr0->GetN();
412  double minX = gr0->GetX()[0];
413  double maxX = gr0->GetX()[nbins-1];
414  //TGraph to center bin
415  minX = minX - (gr0->GetX()[1] - gr0->GetX()[0])/2;
416  maxX = maxX + (gr0->GetX()[1] - gr0->GetX()[0])/2;
417  //Fourier fit relies on periodicity, ignore extra points
418  if(isDelta && !useRoot){
419  nbins = 60;
420  minX = 0;
421  maxX = 2;
422  }
423  std::cout << "Smooth " << nbins << " nbins " << minX << ", " << maxX << "\n";
424  TH1D * hist = new TH1D(UniqueName().c_str(),"",nbins,minX,maxX);
425  for (int i = 0; i < nbins; ++i){
426  double val = gr0->Eval(hist->GetBinCenter(i+1));
427  hist->SetBinContent(i + 1, val * val); //use dchisq not chi2 to smooth
428  }
429  if(useRoot) hist->Smooth(2);
430  else SmoothWithFourierFit(hist, nbins, minX, maxX, fourierFitN);
431  //now return sqrt
432  for (int i = 1; i <= nbins; ++i){
433  hist->SetBinContent(i, sqrt(hist->GetBinContent(i)));
434  }
435  if(isDelta) return HistoToCyclicGraph(hist, useRoot);
436  else return new TGraph(hist);
437 }
438 
439 
440 TGraph* HistoToCyclicGraph(TH1* h, bool useRoot = false)
441 {
442  TGraph* g = new TGraph;
443  g->SetPoint(0, 0, 0); // temporary
444  for(int i = 1; i <= h->GetNbinsX(); ++i)
445  g->SetPoint(i, h->GetXaxis()->GetBinCenter(i), h->GetBinContent(i));
446  g->SetLineColor(h->GetLineColor());
447  g->SetLineWidth(2);
448 
449  if(useRoot)
451  else
452  AddCyclicPoints(g);
453 
454  return g;
455 }
456 void SmoothWithFourierFit(TH1* hist, int nbins, double minX, double maxX, int fourierFitN)
457 {
458  FitToFourier fitF (hist, minX, maxX, fourierFitN);
459  auto f1= fitF.Fit();
460  hist->Reset();
461  hist->Add(f1);
462 }
463 
464 void AddCyclicPoints(TGraph* g)
465 {
466  double x, y;
467  g->GetPoint(g->GetN()-1, x, y);
468  g->SetPoint(0, x-2, y);
469  g->GetPoint(1, x, y);
470  g->SetPoint(g->GetN(), x+2, y);
471  std::cout << "\n Added cyclic points " << 0 << " " << g->GetX()[0]<< " " << g->GetY()[0]
472  << ", " << g->GetN() << x+2 << y << "\n\n";
473 }
474 void AverageCyclicPoints(TGraph* g)
475 {
476  double x0, y0;
477  double x2pi, y2pi;
478  g->GetPoint(g->GetN()-1, x2pi, y2pi);
479  g->GetPoint(1, x0, y0);
480  std::cout << "(x0, y0) = (" << x0 << ", " << y0 << ")\n";
481  std::cout << "(x2pi, y2pi) = (" << x2pi << ", " << y2pi << ")\n";
482 
483  double avgY = ( y0 + y2pi ) / 2;
484  std::cout << "Average = " << avgY << std::endl;
485  g->SetPoint(0, 0, avgY);
486  g->SetPoint(g->GetN(), 2, avgY);
487 
488  std::cout << "\nFound average cyclic point: (0, " << avgY << ")\n";
489  std::cout << "\nAdded average cyclic points " << 0 << ": " << g->GetX()[0] << " " << g->GetY()[0];
490  std::cout << ", " << g->GetN() << ": " << g->GetX()[g->GetN()-1] << " " << g->GetY()[g->GetN()-1] << "\n\n";
491 }
T max(const caf::Proxy< T > &a, T b)
TGraph * FCCorrectSlice(TGraph *sqrtslice, TString fcFileName, bool fourierFit, std::string plot_name, ofstream &summary, TString par, TString hier, TString oct, bool useRoot=false, bool fccol=false)
void AddCyclicPoints(TGraph *g)
double minY
Definition: plot_hist.C:9
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
system("rm -rf microbeam.root")
fvar< T > fabs(const fvar< T > &x)
Definition: fabs.hpp:15
TH1D * ProfileAxes(TString par)
void nova_official_data_release(bool contours=true, bool slices=true)
double delta
Definition: runWimpSim.h:98
TGraph * SmoothWithFourierFit(TGraph *gr0, bool isDelta=true, bool useRoot=false, int fourierFitN=2)
T sqrt(T number)
Definition: d0nt_math.hpp:156
Int_t par
Definition: SimpleIterate.C:24
TString MakeSummaryTable(FCSurface *fsurf, TTree *tfc, TString par, TString hier, TString oct="")
double maxY
Definition: plot_hist.C:10
void CenterTitles(TH1 *histo)
Definition: Plots.cxx:1481
void CylindricalSmoothing(TH2 *h)
const FCBin * GetBin(int x, int y) const
Definition: FCSurface.cxx:152
T sqr(T x)
More efficient square function than pow(x,2)
Definition: MathUtil.h:23
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
Definition: LoadFromFile.h:17
double SignificanceForUpValue(double up) const
Definition: FCBin.cxx:59
const char * label
static std::unique_ptr< FCSurface > LoadFromFile(const std::string &fname)
Definition: FCSurface.cxx:170
double PValueToSigma(double p)
Compute the equivalent number of gaussian sigma for a p-value.
const int nbins
Definition: cellShifts.C:15
Float_t f1
const double j
Definition: BetheBloch.cxx:29
void Add(const FCPoint &pt, std::string plot)
Definition: FCSurface.cxx:39
TGraph * HistoToCyclicGraph(TH1 *h, bool useRoot)
TF1 * Fit() const
Definition: Utilities.cxx:422
OStream cout
Definition: OStream.cxx:6
TString GetFCPath(TString par, TString hier)
ofstream summary
::xsd::cxx::tree::string< char, simple_type > string
Definition: Database.h:154
TFile * GetGaussianProfile(TString par)
const TH2F * Binning() const
Definition: FCSurface.h:30
string shortName
THUMBNAIL BLOCK: We need to make a thumbnail for each.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
Definition: Style.cxx:154
T min(const caf::Proxy< T > &a, T b)
TH2 * SurfaceForSignificance(double sig)
Definition: FCSurface.cxx:103
Float_t e
Definition: plot.C:35
bool Empty() const
Definition: FCBin.h:23
void AverageCyclicPoints(TGraph *g)
static std::unique_ptr< FCCollection > LoadFromFile(const std::string &wildcard)
std::string UniqueName()
Return a different string each time, for creating histograms.
Definition: Utilities.cxx:29
const std::string outdir
TFile * GetGaussianSurface(TString par, TString hier)
Pseudo-experiments, binned to match a log-likelihood surface.
Definition: FCSurface.h:14
surf
Definition: demo5.py:35