6 #include "TLegendEntry.h" 22 ofstream&
summary, TString
par, TString hier, TString oct,
23 bool useRoot =
false,
bool fccol=
false);
31 int nbins = 80,
double minX = 0,
double maxX = 2,
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";
49 if(par.Contains(
"ssth23"))
50 shortName =
"_th23_noOct";
51 else if(par.Contains(
"dmsq"))
55 TFile*
ret = TFile::Open(gaussDir + baseName + shortName +
".root");
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;
70 TString fcFolder =
"/nova/ana/nu_e_ana/Ana2018/FC/";
72 if(par.Contains(
"ssth23"))
73 return fcFolder +
TString::Format(
"%s_%s.root", par.Data(), hier.Data());
75 return fcFolder +
TString::Format(
"%s_%s_%s.root", par.Data(), oct.Data(), hier.Data());
79 TString
par, TString hier, TString oct =
"");
84 summary.open(
"fc_summary_2018/summary_table.tex");
86 TString
outdir =
"NOvA_2018_official_data_release/";
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")) {
99 TFile* official =
new TFile(outdir +
TString::Format(
"NOvA_2018_official_contour_%s.root",
par.Data()),
"recreate");
100 bf->Write(
"NH_best_fit");
102 for(TString hier: {
"NH",
"IH"}) {
106 TTree* tfc = (TTree*)ffc->Get(
"fc");
107 TString rootDir =
par.Contains(
"ssth23dmsq32")?
"th23_dm32_" :
"delta_th23_";
109 TH2* hGauss = (TH2*) fGauss->Get(rootDir+hier);
110 auto mins =* (
TVectorD*) fGauss->Get(
"overall_minchi");
111 double minchi23 = mins[0];
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);
129 fc->
Add(*fccol,
par.Data());
139 if(
par.Contains(
"deltassth23"))
142 fc_surf_3sig->Smooth();
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);
155 double axesMaxX;
double axesMaxY;
156 double axesMinX;
double axesMinY;
157 if(
par.Contains(
"deltassth23")){
166 if(hier.Contains(
"NH")) {
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})";
183 axes->GetXaxis()->SetTitle(xLabel.Data());
184 axes->GetYaxis()->SetTitle(yLabel.Data());
188 TString axesName = hier +
"_axes";
189 if(
par.Contains(
"deltassth23"))
191 if(!(
par.Contains(
"deltassth23") && hier.Contains(
"IH"))){
192 axes->Write(axesName.Data());
194 std::vector<std::pair< TString, std::vector<TGraph*> > >
graphs = {{
"1sigma",g1},
198 for(std::pair<TString, std::vector<TGraph*>> gn : graphs){
200 for(TGraph*
g : gn.second){
214 std::map< TString, double >
maxY = {{
"delta", 5.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");
222 axes->SetMaximum(maxY[
par]);
229 for(TString hier : {
"IH",
"NH"}) {
230 for(TString oct : {
"LO",
"UO"}) {
231 if(par.Contains(
"ssth23") && oct.Contains(
"LO"))
234 if(par.Contains(
"ssth23"))
236 TString
histName =
"slice_" + parShortName +
"_" + hier + oct;
238 TH1D*
h = (TH1D*) fGauss->Get(histName);
239 TGraph* graph =
new TGraph(h);
244 bool useFourier = (par.Contains(
"delta") &&
245 !(hier.Contains(
"NH") && oct.Contains(
"UO")));
253 if(par.Contains(
"dmsq")){
254 for (
int i =0;
i < graph->GetN(); ++
i){
255 graph->SetPoint(
i,1000*fabs(graph->GetX()[
i]),graph->GetY()[
i]);
259 graph->Write(hier+oct);
269 TString
par, TString hier, TString oct)
272 int minNPts = 100000;
273 int totNPts = tfc->GetEntries();
277 for(
int i = 1;
i < fcsurf->
Binning()->GetXaxis()->GetNbins()+1;
i++){
278 for(
int j = 1;
j < fcsurf->
Binning()->GetYaxis()->GetNbins()+1;
j++){
283 int NPts = fcsurf->
GetBin(
i,
j)->NPts();
289 avgNPts = totNPts / double(NBins);
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}$";
297 label =
"$\\delta_{CP}$";
299 label =
"$\\sin^2\\theta_{23}$";
301 label =
"$\\Deltam^2_{32}$";
323 TString xLabel =
"#delta_{CP}/#pi";
324 TString yLabel =
"Significance (#sigma)";
327 if(par.Contains(
"dmsq")){
330 xLabel =
"|#Deltam^{2}_{32}| (10^{-3} eV^{2})";
332 else if(par.Contains(
"th23")){
335 xLabel =
"sin^{2}#theta_{23}";
338 TH1D*
ret =
new TH1D(
"axes",
"", 100, minX, maxX);
339 ret->GetXaxis()->SetTitle(xLabel);
340 ret->GetYaxis()->SetTitle(yLabel);
349 TString
par, TString hier, TString oct,
350 bool useRoot,
bool fccol)
352 bool th23plot = fcFileName.Contains(
"th23");
353 bool dmsqplot = fcFileName.Contains(
"dmsq");
354 TString hiername = fcFileName.Contains(
"nh")?
"NH":
"IH";
357 if (th23plot) {minX = 0.3; maxX = 0.7;}
359 minX = hiername==
"NH" ? 2
e-3 : -3
e-3;
360 maxX = hiername==
"NH" ? 3
e-3 : -2
e-3;
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]));
375 fcXH =
new FCSurface(60,minX,maxX,1,0,1);
376 fcXH->
Add(*fccol, plot_name);
378 TFile* ffc = TFile::Open(fcFileName.Data());
379 TTree* tfc = (TTree*) ffc->Get(
"fc");
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];
393 sigXH->SetPoint(sigXH->GetN(),
delta,
sqrt(upXH));
411 int nbins = gr0->GetN();
412 double minX = gr0->GetX()[0];
413 double maxX = gr0->GetX()[nbins-1];
415 minX = minX - (gr0->GetX()[1] - gr0->GetX()[0])/2;
416 maxX = maxX + (gr0->GetX()[1] - gr0->GetX()[0])/2;
418 if(isDelta && !useRoot){
423 std::cout <<
"Smooth " << nbins <<
" nbins " << minX <<
", " << maxX <<
"\n";
426 double val = gr0->Eval(hist->GetBinCenter(
i+1));
427 hist->SetBinContent(
i + 1, val * val);
429 if(useRoot) hist->Smooth(2);
433 hist->SetBinContent(
i,
sqrt(hist->GetBinContent(
i)));
436 else return new TGraph(hist);
442 TGraph*
g =
new TGraph;
443 g->SetPoint(0, 0, 0);
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());
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";
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";
483 double avgY = ( y0 + y2pi ) / 2;
485 g->SetPoint(0, 0, avgY);
486 g->SetPoint(g->GetN(), 2, avgY);
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";
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)
Cuts and Vars for the 2020 FD DiF Study.
system("rm -rf microbeam.root")
TH1D * ProfileAxes(TString par)
void nova_official_data_release(bool contours=true, bool slices=true)
TGraph * SmoothWithFourierFit(TGraph *gr0, bool isDelta=true, bool useRoot=false, int fourierFitN=2)
TString MakeSummaryTable(FCSurface *fsurf, TTree *tfc, TString par, TString hier, TString oct="")
void CenterTitles(TH1 *histo)
const FCBin * GetBin(int x, int y) const
T sqr(T x)
More efficient square function than pow(x,2)
std::unique_ptr< T > LoadFrom(TDirectory *dir, const std::string &label)
double SignificanceForUpValue(double up) const
static std::unique_ptr< FCSurface > LoadFromFile(const std::string &fname)
double PValueToSigma(double p)
Compute the equivalent number of gaussian sigma for a p-value.
void Add(const FCPoint &pt, std::string plot)
TGraph * HistoToCyclicGraph(TH1 *h, bool useRoot)
TString GetFCPath(TString par, TString hier)
TFile * GetGaussianProfile(TString par)
const TH2F * Binning() const
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)
T min(const caf::Proxy< T > &a, T b)
TH2 * SurfaceForSignificance(double sig)
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.
Pseudo-experiments, binned to match a log-likelihood surface.