20 #include "../joint_fit_2020_loader_tools.h" 21 #include "../joint_fit_2020_datarelease_tools.h" 38 #include "TGraphSmooth.h" 54 void POTtag(
bool both,
bool FHCOnly,
bool RHCOnly);
59 int nbins,
double minX,
double maxX,
63 bool useRoot,
int fourierFitN);
66 bool useRoot,
bool fccol);
70 std::ifstream
infile(fileName);
76 bool runOnGrid =
false,
77 TString
options=
"joint_realData_both",
78 bool th23slice =
false,
79 bool octantSlice =
true,
80 bool dmsqSlice =
false,
81 bool th13Slice =
false,
83 bool fourierFit =
true)
86 bool nueOnly =
options.Contains(
"nueOnly");
87 bool numuOnly =
options.Contains(
"numuOnly");
88 bool joint =
options.Contains(
"joint");
89 assert (nueOnly || numuOnly || joint);
91 bool FHCOnly =
options.Contains(
"FHCOnly");
92 bool RHCOnly =
options.Contains(
"RHCOnly");
93 bool both =
options.Contains(
"both");
94 assert (FHCOnly || RHCOnly || both);
96 bool fake2019 =
options.Contains(
"fake2019");
97 bool realData =
options.Contains(
"realData");
100 if(corrSysts) suffix+=
"_systs";
101 else suffix+=
"_stats";
103 if(th23slice) suffix+=
"_th23";
104 if(dmsqSlice) suffix+=
"_dmsq";
105 if(th13Slice) suffix+=
"_th13";
106 if(!octantSlice)suffix+=
"_noOct";
108 TString
outdir =
"/nova/ana/3flavor/Ana2020/Fits/RealData/WithSeeds/";
109 if (runOnGrid) outdir =
"./";
111 TString outFCfilename (outdir +
"slices_FCInput_2020_" + suffix);
112 if(joint && corrSysts && !th23slice && ! dmsqSlice && !th13Slice) outFCfilename = outdir +
"slices_FCInput_2020_" +suffix;
114 TString
outfilename (outdir +
"hist_slices_2020_" + suffix);
116 TString infilename =
outdir;
118 infilename +=
"hist_slices_2020_";
120 infilename += corrSysts?
"_systs":
"_stats";
121 if(th23slice) infilename+=
"_th23";
122 if(dmsqSlice) infilename+=
"_dmsq";
123 if(th13Slice) infilename+=
"_th13";
124 if(!octantSlice)infilename+=
"_noOct";
125 infilename +=
".root";
128 TFile *
infile =
new TFile (infilename,
"read");
130 TString officialpreffix =
"./";
131 bool makedatarelease =
true;
132 if(makedatarelease &&
133 !(corrSysts &&
options.Contains(
"joint") && fccorr )) {
135 <<
"makedatarelease is set to true" 136 <<
"for a weird combination of options\n";
138 if(makedatarelease) {
140 th23slice, dmsqSlice);
142 auto mins =* (
TVectorD*)infile->Get(
"overall_minchi");
143 double minchi23 = mins[0];
145 std::vector <TString> sliceNames = {};
146 if (!th23slice && ! dmsqSlice && !th13Slice) sliceNames = {
"slice_delta_"};
147 if (th23slice) sliceNames = {
"slice_th23_"};
148 if (dmsqSlice) sliceNames = {
"slice_dmsq_"};
149 if (th13Slice) sliceNames = {
"slice_th13_"};
154 std::vector<double> seeds;
158 std::vector <th23helper> th23seeds;
167 for(TString sliceName:sliceNames){
168 std::vector < std::pair <TGraph*, TString> >
graphs;
170 if(sliceName.Contains(
"delta")) {
171 for (TString hie:{
"IH",
"NH"}){
172 for(
auto const & thseed:th23seeds){
174 auto gr = (TGraph*) infile ->
Get( sliceName + hie + thseed.label);
176 std::cerr <<
"Problem " << sliceName + hie + thseed.label <<
"\n";
180 graphs.push_back({gr, hie +
" " + thseed.label});
188 for (TString hie:{
"NH",
"IH"}) {
189 for(
auto const & thseed:th23seeds){
190 auto gr = (TGraph*) infile ->
Get( sliceName + hie + thseed.label);
191 if(!gr) {
std::cerr <<
"Problem " << sliceName + hie + thseed.label <<
"\n"; }
192 graphs.push_back({gr, hie +
" " + thseed.label});
196 if(sliceName.Contains(
"th23")) {
198 (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
201 if(sliceName.Contains(
"dmsq")) {
204 if(sliceName.Contains(
"th13")) {
209 TString fcFolder=
"/nova/ana/nu_e_ana/Ana2020/FC/slices/";
212 plot_name =
"ssth23";
225 fcFolder +=
"delta/";
229 if(!joint || !corrSysts)
233 for (
auto &gr:graphs) {
235 fcSliceN += gr.second.Contains(
"NH")?
"_nh" :
"_ih";
237 fcSliceN += gr.second.Contains(
"UO")?
"uo" :
"lo";
245 bool useFourier = (fcSliceN.Contains(
"delta") &&
246 !fcSliceN.Contains(
"delta_nhuo"));
249 fourierFit, plot_name, !useFourier,
true);
255 for (
auto &gr:graphs) {
258 if(gr.second.Contains(
"IH")) gr.first->SetLineColor(
cih_dark);
259 if(gr.second.Contains(
"NH")) gr.first->SetLineColor(
cnh);
260 if(gr.second.Contains(
"LO")) gr.first->SetLineStyle(7);
261 gr.first->SetLineWidth(3);
263 for (
int i =0;
i < gr.first->GetN(); ++
i){
264 gr.first->SetPoint(
i,1000*
fabs(gr.first->GetX()[
i]),gr.first->GetY()[
i]);
267 TGraphSmooth *
test =
new TGraphSmooth(
"normal");
268 TGraph * test_sm = test->SmoothKern(gr.first,
"normal", 0.01);
269 test_sm->SetLineColor(gr.first->GetLineColor());
270 test_sm->SetLineStyle(gr.first->GetLineStyle());
271 test_sm->SetLineWidth(3);
272 test_sm->Draw(
"l same");
276 POTtag(both, FHCOnly, RHCOnly);
280 auto leg =
SliceLegend(graphs, !dmsqSlice && !th23slice && !th13Slice, 0.7, dmsqSlice);
285 for(TString
ext: {
".pdf",
".root"}) {
286 gPad->Print(
"./slice_" + suffix +
288 (fccorr?
"_fccorr":
"") +
289 (fourierFit?
"_smooth":
"") +
296 for (
auto &gr:graphs) {
300 gPad->Print(
"./points_slice_" + suffix +
302 (fccorr?
"_fccorr":
"") +
303 (fourierFit?
"_smooth":
"") +
314 TFile * inFCfile =
new TFile(outFCfilename+
".root",
"read");
317 std::vector <TString> strprof = {
"DmSq32",
"SinSq2Theta13",
"SinSqTheta23"};
318 if(th23slice) strprof = {
"DmSq32",
"SinSq2Theta13",
"DeltaCPpi"};
319 if(th13Slice) strprof = {
"DmSq32",
"DeltaCPpi",
"SinSqTheta23"};
320 if(dmsqSlice) strprof = {
"SinSqTheta23",
"SinSq2Theta13",
"DeltaCPpi"};
325 strprof.push_back(
s->ShortName());
328 for(
auto const &
str:strprof){
330 double miny=-1.5,
maxy=1.5;
331 if(
str==
"DeltaCPpi") {miny=0;
maxy=2;}
332 if(
str==
"DmSq32") {miny=2.2;
maxy=2.9;}
333 if(
str==
"SinSq2Theta13") {miny=0.080;
maxy=0.086;}
334 if(
str==
"SinSqTheta23") {miny=0.3;
maxy=0.8;}
336 auto axes =
new TH2F(
"",
";#delta_{CP}/#pi;" +
str, 100, 0, 2, 100, miny,
maxy);
337 if(dmsqSlice)
axes =
new TH2F(
"",
";|#Deltam^{2}_{32}| (10^{-3} eV^{2});" +
str, 100, 2.2, 2.8, 100, miny,
maxy);
338 if(th23slice)
axes =
new TH2F(
"",
";sin^{2}#theta_{23};" +
str, 100, 0.3, 0.7, 100, miny,
maxy);
339 if(th13Slice)
axes =
new TH2F(
"",
";sin^{2}#theta_{13};" +
str, 100, 0.01, 0.2, 100, miny,
maxy);
342 for (TString hie:{
"NH",
"IH"}){
343 for(
auto const & thseed:th23seeds){
344 auto gr = (TGraph*) inFCfile->Get(hie + thseed.label +
"_" +
str);
346 std::cerr <<
"Problem " << hie + thseed.label +
"_" +
str <<
"\n";
350 for (
int i = 0;
i < (
int) gr->GetN(); ++
i){
351 gr->SetPoint(
i, gr->GetX()[
i], 1000*
fabs(gr->GetY()[
i]));
355 for (
int i = 0;
i < (
int) gr->GetN(); ++
i){
356 gr->SetPoint(
i, 1000*
fabs(gr->GetX()[
i]), gr->GetY()[
i]);
361 if(thseed.label.Contains(
"LO")) gr->SetLineStyle(7);
363 gr->SetMarkerColor(gr->GetLineColor());
368 gPad->Print(
"./debug_slice_prof_" + suffix +
"_"+
str +
".pdf");
380 if(t[0] > 0 &&
abs(y - t[0]) > 0.05)
return 999.;
381 return abs(y - t[0]);
385 TMarker *
m =
new TMarker (x, y,29);
386 m->SetMarkerColor(color);
389 TLatex * ltx =
new TLatex(x,y,pstr.Data());
390 ltx->SetTextAngle(45);
391 ltx->SetTextColor(color);
397 double minX = gr->GetX()[0];
398 double maxX = gr->GetX()[gr->GetN()];
400 TF1 *
f1 =
new TF1(
"sign",
signFunc,minX,maxX,1);
401 f1->SetParameter(0,sigma);
402 double xmin = f1->GetMinimumX(lowx,highx);
416 std::vector <SPoint > sig_range;
418 if(grname.Contains(
"dmsq")){
419 sig_range = {{0.,2.3,2.6},
420 {1,2.3,2.4},{1,2.4,2.6},
421 {2,2.2,2.5},{2,2.5,2.7}};
423 if(grname.Contains(
"th23")){
424 sig_range = {{0.,0,0.5},{0,0.5,0.7},{0,0.5,0.5001},
425 {1,0.3,0.45},{1,0.45,0.5},{1,0.5,0.515},{1,0.515,0.55},{1,0.54,0.7},
426 {2,0.3,0.5},{2,0.5,0.7}};
428 if(grname.Contains(
"th13")){
429 sig_range = {{0.,0,0.1},{0,0.1,0.15},{0,0.15,0.2},
430 {1,0.0,0.1},{1,0.1,0.15},{1,0.15,0.2},
431 {2,0.0,0.1},{2,0.1,0.2}};
433 if(grname.Contains(
"delta")){
434 sig_range = {{0.,1., 2.}, {0.,0., 1.},
435 {1., 0, 0.5},{1.,0.5, 1.},{1.,1, 1.7},{1.,1.7, 2.},
436 {2., 0, 0.5},{2.,0.5, 1.},{2.,1, 1.3},{2.,1.5, 2.},
437 {3., 0, 0.5}, {3,0.5,1.5}, {3,1.7,2.}};
440 int color = kMagenta;
441 for (
auto &
sr:sig_range){
442 if(
sr.sigma == 1 ) color =
kGreen + 1;
443 if(
sr.sigma == 2 ) color = kAzure ;
449 void POTtag(
bool both,
bool FHCOnly,
bool RHCOnly)
451 TLatex* ltx =
new TLatex();
457 ltx1->SetTextAlign(11);
461 ltx->SetTextAlign(11);
470 g->GetPoint(g->GetN()-1, x2pi, y2pi);
471 g->GetPoint(1, x0, y0);
472 std::cout <<
"(x0, y0) = (" << x0 <<
", " << y0 <<
")\n";
473 std::cout <<
"(x2pi, y2pi) = (" << x2pi <<
", " << y2pi <<
")\n";
475 double avgY = ( y0 + y2pi ) / 2;
477 g->SetPoint(0, 0, avgY);
478 g->SetPoint(g->GetN(), 2, avgY);
480 std::cout <<
"\nFound average cyclic point: (0, " << avgY <<
")\n";
481 std::cout <<
"\nAdded average cyclic points " << 0 <<
": " << g->GetX()[0] <<
" " << g->GetY()[0];
482 std::cout <<
", " << g->GetN() <<
": " << g->GetX()[g->GetN()-1] <<
" " << g->GetY()[g->GetN()-1] <<
"\n\n";
490 g->GetPoint(g->GetN()-1,
x,
y);
491 g->SetPoint(0, x-2, y);
492 g->GetPoint(1, x, y);
493 g->SetPoint(g->GetN(), x+2,
y);
495 std::cout <<
"\n Added cyclic points " << 0 <<
" " << g->GetX()[0]<<
" " << g->GetY()[0]
496 <<
", " << g->GetN()-1 << x+2 << y <<
"\n\n";
520 int nbins = 80,
double minX = 0,
double maxX = 2,
531 bool useRoot =
false,
534 int nbins = gr0->GetN();
535 double minX = gr0->GetX()[0];
536 double maxX = gr0->GetX()[nbins-1];
538 minX = minX - (gr0->GetX()[1] - gr0->GetX()[0])/2;
539 maxX = maxX + (gr0->GetX()[1] - gr0->GetX()[0])/2;
542 if(isDelta && !useRoot){
547 std::cout <<
"Smooth " << nbins <<
" nbins " << minX <<
", " << maxX <<
"\n";
552 double val = gr0->Eval(hist->GetBinCenter(
i+1));
553 hist->SetBinContent(
i + 1, val * val);
555 if(useRoot) hist->Smooth(2);
559 hist->SetBinContent(
i,
sqrt(hist->GetBinContent(
i)));
563 return new TGraph(hist);
568 bool useRoot =
false,
bool fccol=
false)
570 bool th23plot = fcFileName.Contains(
"th23");
571 bool dmsqplot = fcFileName.Contains(
"dmsq");
572 bool dmsqnhuo = fcFileName.Contains(
"dmsq_nhuo");
573 TString hiername = fcFileName.Contains(
"nh")?
"NH":
"IH";
577 if (th23plot) {minX = 0.3; maxX = 0.7;}
579 minX = hiername==
"NH" ? 2
e-3 : -3
e-3;
580 maxX = hiername==
"NH" ? 3
e-3 : -2
e-3;
584 auto n= sqrtslice->GetN();
586 auto sliceXH =
new TGraph(
n);
587 for (
int i = 0;
i <
n;
i++ ){
588 sliceXH->SetPoint(
i,sqrtslice->GetX()[
i],
util::sqr(sqrtslice->GetY()[
i]));
597 fcXH =
new FCSurface(60,minX,maxX,1,0,1);
598 fcXH->
Add(*fccol, plot_name);
602 TGraph* sigXH =
new TGraph;
603 TGraph* pcXH =
new TGraph;
605 TFile* bestfitfile =
new TFile(
"/nova/ana/nu_e_ana/Ana2019/Results/BestFits/bestfits_joint_realData_both_systs.root",
"read");
610 for(
int i = 0;
i <
n; ++
i){
611 const double delta = sliceXH->GetX()[
i];
612 const double upXH = sliceXH->GetY()[
i];
616 std::cout <<
" Falling back on gaussian up value: " << upXH
617 <<
" for bin " <<
i <<
" , delta = " << delta <<
std::endl;
618 sigXH->SetPoint(sigXH->GetN(),
delta,
sqrt(upXH));
626 if(dmsqnhuo && (
i+1 < n) && (delta < dmsq32bf) && (dmsq32bf < sliceXH->
GetX()[
i+1])){
627 sigXH->SetPoint(sigXH->GetN(), dmsq32bf, 0.);
629 if ( special_smoothing && (
i+1 < n) && (delta < deltabf) && (deltabf < sliceXH->
GetX()[
i+1])){
631 sigXH->SetPoint(sigXH->GetN(), deltabf, 0.);
633 if ( fcFileName.Contains(
"th23_nh") && (
i+1 <
n) && (delta < th23bf) && (th23bf < sliceXH->GetX()[
i+1])){
635 sigXH->SetPoint(sigXH->GetN(), th23bf, 0.);
638 if(fourierFit && !th23plot) {
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2.)
Cuts and Vars for the 2020 FD DiF Study.
void HighlightSignPoints(TGraph *gr, TString grname)
fvar< T > fabs(const fvar< T > &x)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
double GetValue(const osc::IOscCalcAdjustable *osc) const override
Forward to wrapped Var's GetValue()
const FitDmSq32 kFitDmSq32
double GetValue(const osc::IOscCalcAdjustable *osc) const override
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
const FCBin * GetBin(int x, int y) const
T sqr(T x)
More efficient square function than pow(x,2)
Encapsulate code to systematically shift a caf::SRProxy.
const double kAna2020FHCFullDetEquivPOT
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.
string outfilename
knobs that need extra care
bool is_file_exist(const char *fileName)
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta)
void Add(const FCPoint &pt, std::string plot)
const double kAna2020FHCPOT
void AverageCyclicPoints(TGraph *g)
void DrawSignMarker(double x, double y, int color)
const double kAna2020RHCPOT
Double_t signFunc(Double_t *x, Double_t *t)
And supporting functions:
void AddCyclicPoints(TGraph *g)
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
TGraph * FCCorrectSlice(TGraph *sqrtslice, TString fcFileName, bool fourierFit, std::string plot_name, bool useRoot, bool fccol)
void SmoothWithFourierFit(TH1 *hist, int nbins, double minX, double maxX, int fourierFitN)
assert(nhit_max >=nhit_nbins)
void FindSignPoint(TGraph *gr, double sigma, double lowx, double highx, int color)
Interface definition for fittable variables.
void plot_joint_fit_2020_slices(bool corrSysts=false, bool runOnGrid=false, TString options="joint_realData_both", bool th23slice=false, bool octantSlice=true, bool dmsqSlice=false, bool th13Slice=false, bool fccorr=false, bool fourierFit=true)
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
const Float_t kBlessedLabelFontSize
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::vector< const ISyst * > GetJointFitSystematicList(bool corrSysts, bool nueExclusive=false, bool numuExclusive=false, bool isFHC=true, bool isRHC=true, bool intersection=true, bool ptExtrap=true)
static std::unique_ptr< FCCollection > LoadFromFile(const std::string &wildcard)
std::string UniqueName()
Return a different string each time, for creating histograms.
void POTtag(bool both, bool FHCOnly, bool RHCOnly)
TGraph * HistoToCyclicGraph(TGraph *g, bool useRoot)
double GetX(int ndb=14, int db=1, int feb=0, int pix=0)
Pseudo-experiments, binned to match a log-likelihood surface.