51 void MakeMaps ( std::vector <const IFitVar*> profvars,
52 std::map<const IFitVar*, TGraph*> &profMap,
53 std::vector <const ISyst* >
systs,
54 std::map<const ISyst*, TGraph*> &systMap);
57 std::vector <const IFitVar*> profvars,
58 std::vector <TString > profnames,
59 std::map<const IFitVar*, TGraph*> &profMap,
60 std::vector <const ISyst* >
systs,
61 std::map<const ISyst*, TGraph*> &systMap);
62 void POTtag(
bool both,
bool FHCOnly,
bool RHCOnly);
67 int nbins,
double minX,
double maxX,
71 bool useRoot,
int fourierFitN);
74 bool useRoot,
bool fccol);
79 bool corrSysts =
true,
80 bool runOnGrid =
false,
81 TString
options=
"joint_realData_both",
82 bool th23slice =
false,
83 bool octantSlice =
true,
84 bool dmsqSlice =
false,
85 bool th13Slice =
false,
87 bool fourierFit =
true)
90 bool nueOnly =
options.Contains(
"nueOnly");
91 bool numuOnly =
options.Contains(
"numuOnly");
92 bool joint =
options.Contains(
"joint");
93 assert (nueOnly || numuOnly || joint);
95 bool FHCOnly =
options.Contains(
"FHCOnly");
96 bool RHCOnly =
options.Contains(
"RHCOnly");
97 bool both =
options.Contains(
"both");
98 assert (FHCOnly || RHCOnly || both);
100 bool fake2018 =
options.Contains(
"fake2018");
101 bool realData =
options.Contains(
"realData");
104 if(corrSysts) suffix+=
"_systs";
105 else suffix+=
"_stats";
108 if(FHCOnly) tag=
"FHCOnly/";
109 if(RHCOnly) tag=
"RHCOnly/";
110 if(both) tag =
"RHC_and_FHC/";
112 if(th23slice) suffix+=
"_th23";
113 if(dmsqSlice) suffix+=
"_dmsq";
114 if(th13Slice) suffix+=
"_th13";
115 if(!octantSlice)suffix+=
"_noOct";
118 TString
outdir =
"/nova/ana/nu_e_ana/Ana2019/Results/"+tag+
"slices/FCInputs/";
119 if (runOnGrid) outdir =
"./";
121 TString outFCfilename (outdir +
"slices_FCInput_2019_" + suffix);
122 if(joint && corrSysts && !th23slice && ! dmsqSlice && !th13Slice) outFCfilename = outdir +
"slices_FCInput_2019_" +suffix;
125 outdir=
"./results/blessed_plots/slices/";
126 outdir += corrSysts?
"syst/":
"stat/";
127 if (runOnGrid) outdir =
"./";
129 TString
outfilename (outdir +
"hist_slices_2019_" + suffix);
141 std::pair <TH1D*, double>
cos;
145 std::vector <predictions> preds;
146 std::vector <predictions> preds_numu_only;
147 std::vector <Spectrum * >
data;
148 std::vector <Spectrum * > data_numu_only;
149 std::vector <const IExperiment * > expts;
150 std::vector <const IExperiment * > expts_numu_only;
153 if(fake2018)
SetFakeCalc(calc_fake, 0.58, 2.51
e-3, 0.17);
154 else if(!realData) {
std::cerr <<
"need setting for data\n";
exit(1);}
157 if(FHCOnly || both ) {
160 if(RHCOnly || both ) {
167 if(FHCOnly || both ) {
170 for (
int i = 0;
i< nnumu;
i++ ){
175 if(RHCOnly || both ) {
178 for (
int i = 0;
i< nnumu;
i++ ){
188 if(FHCOnly || both ) data.push_back(
GetNueData2019(
"fhc", runOnGrid));
189 if(RHCOnly || both ) data.push_back(
GetNueData2019(
"rhc", runOnGrid));
192 if(FHCOnly || both ){
194 data.insert(data.end(),numu_data.begin(), numu_data.end());
195 data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());}
196 if(RHCOnly || both ){
198 data.insert(data.end(),numu_data.begin(), numu_data.end());
199 data_numu_only.insert(data_numu_only.end(),numu_data.begin(), numu_data.end());}
203 for(
int i = 0;
i <
int(preds.size()); ++
i){
204 double POT = preds[
i].pot;
206 if(realData) thisdata = data[
i];
207 cout<<i<<
" "<<preds[
i].name<<
" POT "<<POT<<
" tot MC "<<preds[
i].pred->Predict(calc_fake).Integral(POT)<<
" cos "<<preds[
i].cos.first->Integral()<<
" cos er "<<preds[
i].cos.second<<
" analyze data "<<thisdata->Integral(POT)<<
endl;
211 cout<<
"Make numu only experiment to get the seeds"<<
endl;
214 for(
int i = 0;
i < (
int) preds_numu_only.size(); ++
i){
215 double POT = preds_numu_only[
i].pot;
216 auto thisdata =
GetFakeData (preds_numu_only[
i].
pred, calc_fake, POT, preds_numu_only[
i].
cos.first);
217 if(realData) thisdata = data_numu_only[
i];
218 cout<<i<<
" "<<preds_numu_only[
i].name<<
" POT "<<POT<<
" tot MC "<<preds_numu_only[
i].pred->Predict(calc_fake).Integral(POT)<<
" cos "<<preds_numu_only[
i].cos.first->Integral()<<
" cos er "<<preds_numu_only[
i].cos.second<<
" analyze data "<<thisdata->Integral(POT)<<
endl;
219 expts_numu_only.push_back(
new SingleSampleExperiment(preds_numu_only[i].pred, *thisdata, preds_numu_only[i].
cos.first, preds_numu_only[i].cos.second));
227 else {
std::cout <<
"No reactor constraint, that's th13 slice\n";}
229 std::cout <<
"Adding Dmsq32ConstraintPDG2017\n";
232 std::cout <<
"Creating Multiexperiment with a total of " 233 << expts.size() <<
" experiments\n\n" <<
std::endl;
236 std::cout <<
"Creating Multiexperiment of numu only SimpleExp with a total of " 237 << expts_numu_only.size() <<
" experiments\n\n" <<
std::endl;
244 std::cout <<
"Systematics for the fit:\n";
247 std::cout <<
"\n\nSystematic correlations...\n";
248 if(!nueOnly && ! numuOnly && corrSysts){
252 for(
int i =0;
i < nnumu; ++
i) exptThis->SetSystCorrelations(
i+1, notfornumu);
257 for(
int i =0;
i < nnumu; ++
i) exptThis->SetSystCorrelations(
i+1, notfornumu);
263 for(
int i =0;
i < 8; ++
i) exptThis->SetSystCorrelations(
i+2, notfornumu);
266 if (nueOnly && corrSysts){
267 if (FHCOnly) exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
true));
268 if (RHCOnly) exptThis->SetSystCorrelations(0,
GetCorrelations(
true,
false));
275 std::cout <<
"Systematics for the numu only fit:\n";
288 std::vector <SystShifts> seedShifts = {};
289 if(corrSysts && !th23slice && !dmsqSlice && !th13Slice && octantSlice && (both || FHCOnly)){
290 cout<<
"Add Cherenkov seeds \n";
291 for (
double systshift:{-0.5, +0.5}){
293 seedShifts.push_back(tempShifts);
296 cout<<
"Add neutron systematic seeds \n";
297 for (
double systshift:{-0.5, +0.5}){
299 seedShifts.push_back(tempShifts);
309 cout<<
"------------------- Start prestage seeds --------------------------"<<
endl;
311 double minchi_numu = 1E20;
312 double pre_seed_th23;
313 double pre_seed_dmsq;
317 double maxmix = 0.514;
318 double numu_pre_seedLONH, numu_pre_seedUONH, numu_pre_seedLOIH, numu_pre_seedUOIH, dmsq_numu_pre_seedNH, dmsq_numu_pre_seedIH;
321 std::cout <<
"\n\nFinding test best fit " << (hie>0?
"NH " :
"IH ") <<
"\n\n";
324 auto thisminchi = fitnumu_only.Fit(calc, auxShifts ,
329 if(thisminchi < minchi_numu) minchi_numu = thisminchi;
332 pre_seed_dmsq = kFitDmSq32.GetValue(calc);
335 numu_pre_seedLONH = ((pre_seed_th23>maxmix)?(2*maxmix-pre_seed_th23):pre_seed_th23);
336 numu_pre_seedUONH = ((pre_seed_th23>maxmix)?pre_seed_th23:(2*maxmix-pre_seed_th23));
341 cout<<
"------------------- End prestage seeds --------------------------"<<
endl;
345 std::vector<double> seeds;
351 std::vector <th23helper> th23seeds;
363 else th23seeds.push_back({ {numu_pre_seedLONH, 0.5, numu_pre_seedUONH}, &
kFitSinSqTheta23,
""});
365 std::vector<double> delta_seeds = {0, 0.5, 1., 1.5};
366 std::vector<double> th23_old_seeds = {0.45, 0.5, 0.55};
373 double minchi23 = 1E20;
374 for(
int hie:{-1, 1}){
375 for (
auto & thseed:th23seeds){
377 std::cout <<
"\n\nFinding best fit " << (hie > 0 ?
"NH " :
"IH ")
378 << thseed.label <<
", " 388 Fitter fit23(exptThis, fitvars,
systs);
389 cout<<
" pre dmsq seed is "<<pre_seed_dmsq<<
endl;
390 auto thisminchi = fit23.Fit(calc, auxShifts,
393 { thseed.var, thseed.seeds },
397 minchi23=
min(minchi23, thisminchi);
402 TString
str =
"Best fit " ;
403 for (
auto &
v:fitvars){
407 shifts->SetTitle(str);
409 TLine *
l=
new TLine(gPad->GetUxmin(),0,gPad->GetUxmax(),0);
411 gPad->Print(outdir +
"debug_slice_shifts_bestfit_" + suffix +
412 (hie>0?
"NH":
"IH") + thseed.label +
".pdf");
419 std::cout <<
"\nFound overall minchi " << minchi23 <<
"\n\n";
421 cout<<
"Check with oldstyle seeds "<<
endl;
422 double minchi23test = 1E20;
424 std::cout <<
"\n\nFinding test best fit " 425 << (hie>0?
"NH " :
"IH ")
432 Fitter fit23(exptThis, fitvars,
systs);
433 auto thisminchi = fit23.Fit(calc,auxShifts ,{
439 minchi23test=
min(minchi23test, thisminchi);
442 TString
str =
"Best fit " ;
443 for (
auto &
v:fitvars){
447 shifts->SetTitle(str);
448 gPad->Print(outdir +
"debug_slice_shifts_bestfit_" + suffix +
449 (hie>0?
"NH":
"IH") +
"test_chisq_with_no_octant_differentiation.pdf");
458 cout<<
"<<<<<<<<<<<<<<<<<<<<<<<<<< Test with only systematics free"<<
endl;
468 Fitter fit23(exptThis, {},
systs);
469 auto thisminchi = fit23.Fit(calc,auxShifts ,{}, seedShifts,
471 ofstream txtfile(outdir +
"debug_slice_shifts_bestfit_" + suffix +
472 +
"test_chisq_with_no_octant_test_with_osc_par_fixed.txt");
474 txtfile<<
v->ShortName().c_str()<<
"="<<auxShifts.
GetShift(
v)<<
endl;
477 minchi23test=
min(minchi23test, thisminchi);
479 TString
str =
"Best fit " ;
481 shifts->SetTitle(str);
482 gPad->Print(outdir +
"debug_slice_shifts_bestfit_" + suffix +
483 +
"test_chisq_with_no_octant_test_with_osc_par_fixed.pdf");
486 cout<<
"<<<<<<<<<<<< This is not true BF : "<< minchi23test<<
endl;
493 TFile *
outfile =
new TFile(outfilename+
".root",
"recreate");
494 TFile * outFCfile =
new TFile(outFCfilename+
".root",
"recreate");
499 v.Write(
"overall_minchi");
504 for(
int hie: {-1, +1}){
509 std::cout <<
"Starting profile " << (hie>0?
"NH ":
"IH") <<
"\n\n";
515 std::vector <const IFitVar * > profvars;
516 std::vector <TString> profvarnames;
517 std::map <const IFitVar *, std::vector <double> >profseeds;
521 for(
auto const & thseed:th23seeds){
523 const IFitVar* kCurTheta23 = thseed.var;
525 if(!th23slice && !dmsqSlice && !th13Slice){
530 profdef.profvarnames = {
"DmSq32",
"SinSq2Theta13",
"SinSqTheta23"};
532 profdef.shortname=
"delta";
538 profdef.minX = (nueOnly ? 0:0.3);
539 profdef.maxX = (nueOnly ? 1:0.7);
541 profdef.profvarnames = {
"DmSq32",
"SinSq2Theta13",
"DeltaCPpi"};
543 profdef.shortname=
"th23";
551 profdef.profvarnames = {
"DmSq32",
"SinSqTheta23",
"DeltaCPpi"};
553 profdef.shortname=
"th13";
558 profdef.minX = hie > 0 ? 2
e-3 : -3
e-3;
559 profdef.maxX = hie > 0 ? 3
e-3 : -2
e-3;
561 profdef.profvarnames = {
"SinSqTheta23",
"SinSq2Theta13",
"DeltaCPpi"};
562 profdef.profseeds = { {kCurTheta23, thseed.seeds}, {&
kFitDeltaInPiUnits,delta_seeds}};
563 profdef.shortname=
"dmsq";
566 std::map<const IFitVar*, TGraph*> profVarsMap;
567 std::map<const ISyst*, TGraph*> systMap;
568 MakeMaps (profdef.profvars, profVarsMap, systs, systMap);
570 std::cout <<
" Profile "<< thseed.label <<
"\n";
573 profdef.fitvar, steps, profdef.minX, profdef.maxX,
579 profVarsMap, systMap);
581 TString hieroctstr = (hie > 0 ?
"NH" :
"IH") + thseed.label;
584 slice->Write((TString )
"slice_" + profdef.shortname +
"_" + hieroctstr);
586 profdef.profvars, profdef.profvarnames, profVarsMap, systs, systMap);
600 TString infilename =
"/nova/ana/nu_e_ana/Ana2019/Results/RHC_and_FHC/slices/";
601 infilename += corrSysts?
"syst/":
"stat/";
602 infilename +=
"hist_slices_2019_";
604 infilename += corrSysts?
"_systs":
"_stats";
605 if(th23slice) infilename+=
"_th23";
606 if(dmsqSlice) infilename+=
"_dmsq";
607 if(th13Slice) infilename+=
"_th13";
608 if(!octantSlice)infilename+=
"_noOct";
609 infilename +=
".root";
612 TFile *
infile =
new TFile (infilename,
"read");
614 auto mins =* (
TVectorD*)infile->Get(
"overall_minchi");
615 double minchi23 = mins[0];
617 std::vector <TString> sliceNames = {};
618 if (!th23slice && ! dmsqSlice && !th13Slice) sliceNames = {
"slice_delta_"};
619 if (th23slice) sliceNames = {
"slice_th23_"};
620 if (dmsqSlice) sliceNames = {
"slice_dmsq_"};
621 if (th13Slice) sliceNames = {
"slice_th13_"};
626 std::vector<double> seeds;
630 std::vector <th23helper> th23seeds;
639 for(TString sliceName:sliceNames){
640 std::vector < std::pair <TGraph*, TString> >
graphs;
642 if(sliceName.Contains(
"delta")) {
643 for (TString hie:{
"NH",
"IH"}){
644 for(
auto const & thseed:th23seeds){
645 auto h = (TH1D*) infile ->
Get( sliceName + hie + thseed.label);
647 std::cerr <<
"Problem " << sliceName + hie + thseed.label <<
"\n";
651 graphs.push_back({
new TGraph(
h), hie +
" " + thseed.label});
659 for (TString hie:{
"NH",
"IH"}) {
660 for(
auto const & thseed:th23seeds){
661 auto h = (TH1D*) infile ->
Get( sliceName + hie + thseed.label);
662 if(!
h) {
std::cerr <<
"Problem " << sliceName + hie + thseed.label <<
"\n"; }
663 graphs.push_back({
new TGraph(
h),hie +
" " + thseed.label});
667 if(sliceName.Contains(
"th23")) {
669 (nueOnly ? 0 : 0.32), (nueOnly ? 1 : 0.72));
672 if(sliceName.Contains(
"dmsq")) {
675 if(sliceName.Contains(
"th13")) {
680 TString fcFolder=
"/nova/ana/nu_e_ana/Ana2019/FC/slices/";
683 plot_name =
"ssth23";
696 fcFolder +=
"delta/";
700 if(!joint || !corrSysts)
704 for (
auto &gr:graphs) {
706 fcSliceN += gr.second.Contains(
"NH")?
"_nh" :
"_ih";
708 fcSliceN += gr.second.Contains(
"UO")?
"uo" :
"lo";
714 bool useFourier = (fcSliceN.Contains(
"delta") &&
715 !fcSliceN.Contains(
"delta_nhuo"));
718 fourierFit, plot_name, !useFourier,
true);
723 for (
auto &gr:graphs) {
726 if(gr.second.Contains(
"IH")) gr.first->SetLineColor(
cih_dark);
727 if(gr.second.Contains(
"NH")) gr.first->SetLineColor(
cnh);
728 if(gr.second.Contains(
"LO")) gr.first->SetLineStyle(7);
729 gr.first->SetLineWidth(3);
731 for (
int i =0;
i < gr.first->GetN(); ++
i){
732 gr.first->SetPoint(
i,1000*
fabs(gr.first->GetX()[
i]),gr.first->GetY()[
i]);
735 gr.first->Draw(
"l same");
738 POTtag(both, FHCOnly, RHCOnly);
742 auto leg =
SliceLegend(graphs, !dmsqSlice && !th23slice && !th13Slice, 0.7, dmsqSlice);
747 for(TString
ext: {
".pdf",
".root"}) {
748 gPad->Print(outdir +
"slice_" + suffix +
750 (fccorr?
"_fccorr":
"") +
751 (fourierFit?
"_smooth":
"") +
758 for (
auto &gr:graphs) {
763 gPad->Print(outdir +
"points_slice_" + suffix +
765 (fccorr?
"_fccorr":
"") +
766 (fourierFit?
"_smooth":
"") +
776 TFile * inFCfile =
new TFile(outFCfilename+
".root",
"read");
779 std::vector <TString> strprof = {
"DmSq32",
"SinSq2Theta13",
"SinSqTheta23"};
780 if(th23slice) strprof = {
"DmSq32",
"SinSq2Theta13",
"DeltaCPpi"};
781 if(th13Slice) strprof = {
"DmSq32",
"DeltaCPpi",
"SinSqTheta23"};
782 if(dmsqSlice) strprof = {
"SinSqTheta23",
"SinSq2Theta13",
"DeltaCPpi"};
787 strprof.push_back(
s->ShortName());
790 for(
auto const &
str:strprof){
792 double miny=-1.5,
maxy=1.5;
793 if(
str==
"DeltaCPpi") {miny=0;
maxy=2;}
794 if(
str==
"DmSq32") {miny=2.2;
maxy=2.9;}
795 if(
str==
"SinSq2Theta13") {miny=0.080;
maxy=0.086;}
796 if(
str==
"SinSqTheta23") {miny=0.3;
maxy=0.8;}
798 auto axes =
new TH2F(
"",
";#delta_{CP}/#pi;" +
str, 100, 0, 2, 100, miny,
maxy);
799 if(dmsqSlice)
axes =
new TH2F(
"",
";|#Deltam^{2}_{32}| (10^{-3} eV^{2});" +
str, 100, 2.2, 2.8, 100, miny,
maxy);
800 if(th23slice)
axes =
new TH2F(
"",
";sin^{2}#theta_{23};" +
str, 100, 0.3, 0.7, 100, miny,
maxy);
801 if(th13Slice)
axes =
new TH2F(
"",
";sin^{2}#theta_{13};" +
str, 100, 0.01, 0.2, 100, miny,
maxy);
804 for (TString hie:{
"NH",
"IH"}){
805 for(
auto const & thseed:th23seeds){
806 auto gr = (TGraph*) inFCfile->Get(hie + thseed.label +
"_" +
str);
808 std::cerr <<
"Problem " << hie + thseed.label +
"_" +
str <<
"\n";
812 for (
int i = 0;
i < (
int) gr->GetN(); ++
i){
813 gr->SetPoint(
i, gr->GetX()[
i], 1000*
fabs(gr->GetY()[
i]));
817 for (
int i = 0;
i < (
int) gr->GetN(); ++
i){
818 gr->SetPoint(
i, 1000*
fabs(gr->GetX()[
i]), gr->GetY()[
i]);
823 if(thseed.label.Contains(
"LO")) gr->SetLineStyle(7);
825 gr->SetMarkerColor(gr->GetLineColor());
830 gPad->Print(outdir +
"debug_slice_prof_" + suffix +
"_"+
str +
".pdf");
843 if(t[0] > 0 &&
abs(y - t[0]) > 0.05)
return 999.;
844 return abs(y - t[0]);
848 TMarker *
m =
new TMarker (x, y,29);
849 m->SetMarkerColor(color);
852 TLatex * ltx =
new TLatex(x,y,pstr.Data());
853 ltx->SetTextAngle(45);
854 ltx->SetTextColor(color);
860 double minX = gr->GetX()[0];
861 double maxX = gr->GetX()[gr->GetN()];
863 TF1 *
f1 =
new TF1(
"sign",
signFunc,minX,maxX,1);
864 f1->SetParameter(0,sigma);
865 double xmin = f1->GetMinimumX(lowx,highx);
879 std::vector <SPoint > sig_range;
881 if(grname.Contains(
"dmsq")){
882 sig_range = {{0.,2.3,2.6},
883 {1,2.3,2.4},{1,2.4,2.6},
884 {2,2.2,2.5},{2,2.5,2.7}};
886 if(grname.Contains(
"th23")){
887 sig_range = {{0.,0,0.5},{0,0.5,0.7},{0,0.5,0.5001},
888 {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},
889 {2,0.3,0.5},{2,0.5,0.7}};
891 if(grname.Contains(
"th13")){
892 sig_range = {{0.,0,0.1},{0,0.1,0.15},{0,0.15,0.2},
893 {1,0.0,0.1},{1,0.1,0.15},{1,0.15,0.2},
894 {2,0.0,0.1},{2,0.1,0.2}};
896 if(grname.Contains(
"delta")){
897 sig_range = {{0.,1., 2.}, {0.,0., 1.},
898 {1., 0, 0.5},{1.,0.5, 1.},{1.,1, 1.7},{1.,1.7, 2.},
899 {2., 0, 0.5},{2.,0.5, 1.},{2.,1, 1.3},{2.,1.5, 2.},
900 {3., 0, 0.5}, {3,0.5,1.5}, {3,1.7,2.}};
903 int color = kMagenta;
904 for (
auto &
sr:sig_range){
905 if(
sr.sigma == 1 ) color =
kGreen + 1;
906 if(
sr.sigma == 2 ) color = kAzure ;
913 void MakeMaps ( std::vector <const IFitVar*> profvars,
914 std::map<const IFitVar*, TGraph*> &profMap,
915 std::vector <const ISyst* >
systs,
916 std::map<const ISyst*, TGraph*> &systMap)
920 profMap.insert(std::pair<const IFitVar*, TGraph*> (
var,
new TGraph()));
921 for (
const ISyst* syst : systs)
922 systMap.insert(std::pair<const ISyst*, TGraph*> (syst,
new TGraph()));
927 std::vector <const IFitVar*> profvars,
928 std::vector <TString > profnames,
929 std::map<const IFitVar*, TGraph*> &profMap,
930 std::vector <const ISyst* >
systs,
931 std::map<const ISyst*, TGraph*> &systMap)
933 TDirectory *
tmp = gDirectory;
935 for (
int i = 0;
i < (
int) profvars.size(); ++
i){
936 profMap[profvars[
i]]->Write((prefix +
"_" + profnames[
i]));
938 for (
const ISyst*
s : systs)
939 systMap[
s]->Write((prefix +
"_" +
s->ShortName()));
943 void POTtag(
bool both,
bool FHCOnly,
bool RHCOnly)
945 TLatex* ltx =
new TLatex();
951 ltx1->SetTextAlign(11);
955 ltx->SetTextAlign(11);
964 g->GetPoint(g->GetN()-1, x2pi, y2pi);
965 g->GetPoint(1, x0, y0);
966 std::cout <<
"(x0, y0) = (" << x0 <<
", " << y0 <<
")\n";
967 std::cout <<
"(x2pi, y2pi) = (" << x2pi <<
", " << y2pi <<
")\n";
969 double avgY = ( y0 + y2pi ) / 2;
971 g->SetPoint(0, 0, avgY);
972 g->SetPoint(g->GetN(), 2, avgY);
974 std::cout <<
"\nFound average cyclic point: (0, " << avgY <<
")\n";
975 std::cout <<
"\nAdded average cyclic points " << 0 <<
": " << g->GetX()[0] <<
" " << g->GetY()[0];
976 std::cout <<
", " << g->GetN() <<
": " << g->GetX()[g->GetN()-1] <<
" " << g->GetY()[g->GetN()-1] <<
"\n\n";
984 g->GetPoint(g->GetN()-1,
x,
y);
985 g->SetPoint(0, x-2, y);
986 g->GetPoint(1, x, y);
987 g->SetPoint(g->GetN(), x+2,
y);
989 std::cout <<
"\n Added cyclic points " << 0 <<
" " << g->GetX()[0]<<
" " << g->GetY()[0]
990 <<
", " << g->GetN()-1 << x+2 << y <<
"\n\n";
996 TGraph*
g =
new TGraph;
997 g->SetPoint(0, 0, 0);
998 for(
int i = 1;
i <= h->GetNbinsX(); ++
i)
999 g->SetPoint(
i, h->GetXaxis()->GetBinCenter(
i), h->GetBinContent(
i));
1001 g->SetLineColor(h->GetLineColor());
1014 int nbins = 80,
double minX = 0,
double maxX = 2,
1018 auto f1= fitF.
Fit();
1024 bool isDelta =
true,
1025 bool useRoot =
false,
1026 int fourierFitN = 2)
1028 int nbins = gr0->GetN();
1029 double minX = gr0->GetX()[0];
1030 double maxX = gr0->GetX()[nbins-1];
1032 minX = minX - (gr0->GetX()[1] - gr0->GetX()[0])/2;
1033 maxX = maxX + (gr0->GetX()[1] - gr0->GetX()[0])/2;
1036 if(isDelta && !useRoot){
1041 std::cout <<
"Smooth " << nbins <<
" nbins " << minX <<
", " << maxX <<
"\n";
1046 double val = gr0->Eval(hist->GetBinCenter(
i+1));
1047 hist->SetBinContent(
i + 1, val * val);
1049 if(useRoot) hist->Smooth(2);
1053 hist->SetBinContent(
i,
sqrt(hist->GetBinContent(
i)));
1057 else return new TGraph(hist);
1062 bool useRoot =
false,
bool fccol=
false)
1064 bool th23plot = fcFileName.Contains(
"th23");
1065 bool dmsqplot = fcFileName.Contains(
"dmsq");
1066 bool dmsqnhuo = fcFileName.Contains(
"dmsq_nhuo");
1067 TString hiername = fcFileName.Contains(
"nh")?
"NH":
"IH";
1071 if (th23plot) {minX = 0.3; maxX = 0.7;}
1073 minX = hiername==
"NH" ? 2
e-3 : -3
e-3;
1074 maxX = hiername==
"NH" ? 3
e-3 : -2
e-3;
1078 auto n= sqrtslice->GetN();
1080 auto sliceXH =
new TGraph(
n);
1081 for (
int i = 0;
i <
n;
i++ ){
1082 sliceXH->SetPoint(
i,sqrtslice->GetX()[
i],
util::sqr(sqrtslice->GetY()[
i]));
1091 fcXH =
new FCSurface(60,minX,maxX,1,0,1);
1092 fcXH->
Add(*fccol, plot_name);
1096 TGraph* sigXH =
new TGraph;
1097 TGraph* pcXH =
new TGraph;
1099 TFile* bestfitfile =
new TFile(
"/nova/ana/nu_e_ana/Ana2019/Results/BestFits/bestfits_joint_realData_both_systs.root",
"read");
1102 for(
int i = 0;
i <
n; ++
i){
1103 const double delta = sliceXH->GetX()[
i];
1104 const double upXH = sliceXH->GetY()[
i];
1108 std::cout <<
" Falling back on gaussian up value: " << upXH
1109 <<
" for bin " <<
i <<
" , delta = " << delta <<
std::endl;
1110 sigXH->SetPoint(sigXH->GetN(),
delta,
sqrt(upXH));
1118 if(dmsqnhuo && (
i+1 < n) && (delta < dmsq32bf) && (dmsq32bf < sliceXH->
GetX()[
i+1])){
1119 sigXH->SetPoint(sigXH->GetN(), dmsq32bf, 0.);
void DrawSliceCanvas(TString slicename, const double ymax, const double xmin=0, const double xmax=2.)
void joint_fit_2019_slices(bool createFile=false, bool corrSysts=true, 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)
Cuts and Vars for the 2020 FD DiF Study.
void ResetOscCalcToDefault(osc::IOscCalcAdjustable *calc)
Reset calculator to default assumptions for all parameters.
TH1 * PlotSystShifts(const SystShifts &shifts, bool sortName)
fvar< T > fabs(const fvar< T > &x)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
TGraph * HistoToCyclicGraph(TH1 *h, bool useRoot)
const FitDmSq32 kFitDmSq32
Simple record of shifts applied to systematic parameters.
const Dmsq32Constraint kDmsq32ConstraintPDG2017(2.45e-3, 0.05e-3, 2.52e-3, 0.05e-3)
double GetValue(const osc::IOscCalcAdjustable *osc) const override
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
void SaveMaps(TDirectory *out, std::map< std::string, IDecomp * > decomps_nominal, std::map< std::string, std::map< std::string, std::map< int, IDecomp * > > > decomps_shifted, std::map< std::string, PredictionNoExtrap * > predsNE_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionNoExtrap * > > > predsNE_shifted, std::map< std::string, PredictionSterile * > predsSt_nominal, std::map< std::string, std::map< std::string, std::map< int, PredictionSterile * > > > predsSt_shifted)
Save all of the objects in the input maps to the out directory/file.
void MakeMaps(std::vector< const IFitVar * > profvars, std::map< const IFitVar *, TGraph * > &profMap, std::vector< const ISyst * > systs, std::map< const ISyst *, TGraph * > &systMap)
virtual void SetTh13(const T &th13)=0
void SmoothWithFourierFit(TH1 *hist, int nbins, double minX, double maxX, int fourierFitN)
cout<< "Opened file "<< fin<< " ixs= "<< ixs<< endl;if(ixs==0) hhh=(TH1F *) fff-> Get("h1")
const NeutronVisEScaleSyst2018 kNeutronVisEScalePrimariesSyst2018(false)
void FindSignPoint(TGraph *gr, double sigma, double lowx, double highx, int color)
static SystShifts Nominal()
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
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.
virtual void SetDmsq32(const T &dmsq32)=0
double SignificanceForUpValue(double up) const
const XML_Char const XML_Char * data
static std::unique_ptr< FCSurface > LoadFromFile(const std::string &fname)
TGraph * FCCorrectSlice(TGraph *sqrtslice, TString fcFileName, bool fourierFit, std::string plot_name, bool useRoot, bool fccol)
double PValueToSigma(double p)
Compute the equivalent number of gaussian sigma for a p-value.
void DrawSignMarker(double x, double y, int color)
string outfilename
knobs that need extra care
const DummyAnaSyst kAnaCherenkovSyst("Cherenkov","Cherenkov")
T GetShift(const ISyst *syst) const
const ReactorExperiment * WorldReactorConstraint2017()
Reactor constraint from PDG 2017.
virtual T GetDmsq32() const
TLegend * SliceLegend(std::vector< std::pair< TGraph *, TString > > &graphs, bool isDelta)
void Add(const FCPoint &pt, std::string plot)
void HighlightSignPoints(TGraph *gr, TString grname)
std::vector< double > POT
void AverageCyclicPoints(TGraph *g)
static float min(const float a, const float b, const float c)
const double kAna2019RHCPOT
Combine multiple component experiments.
const double kAna2017FullDetEquivPOT
void AddCyclicPoints(TGraph *g)
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
assert(nhit_max >=nhit_nbins)
Interface definition for fittable variables.
double GetPOT(bool isFHC=true)
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
Standard interface to all prediction techniques.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
const Float_t kBlessedLabelFontSize
Double_t signFunc(Double_t *x, Double_t *t)
And supporting functions:
const FitSinSq2Theta13 kFitSinSq2Theta13
std::string to_string(ModuleType mt)
std::pair< Spectrum *, double > cos
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)
void POTtag(bool both, bool FHCOnly, bool RHCOnly)
static std::unique_ptr< FCCollection > LoadFromFile(const std::string &wildcard)
std::string UniqueName()
Return a different string each time, for creating histograms.
double GetX(int ndb=14, int db=1, int feb=0, int pix=0)
Pseudo-experiments, binned to match a log-likelihood surface.
virtual void SetdCP(const T &dCP)=0
Compare a single data spectrum to the MC + cosmics expectation.
TGraph * SqrtProfile(const IExperiment *expt, osc::IOscCalcAdjustable *calc, const IFitVar *v, int nbinsx, double minx, double maxx, double minchi, std::vector< const IFitVar * > profVars, std::vector< const ISyst * > profSysts, const SeedList &seedPts, const std::vector< SystShifts > &systSeedPts, std::map< const IFitVar *, TGraph * > &profVarsMap, std::map< const ISyst *, TGraph * > &systsMap, MinuitFitter::FitOpts opts)
Forward to Profile but sqrt the result for a crude significance.