23 #include "Utilities/func/MathUtil.h" 47 TString
subdir = par.Contains(
"ssth23dmsq32")?
"th23_dmsq" :
"delta_th23";
48 TString gaussDir =
"/nova/ana/nu_e_ana/Ana2018/Results/RHC_and_FHC/contours/"+subdir+
"/syst/";
49 TString gaussName =
"hist_contours_2018_joint_realData_both_only";
51 if(par.Contains(
"ssth23dmsq32"))
52 gaussName +=
"combo_dmsq_systs.root";
54 gaussName +=
"combo_systs.root";
58 TFile*
ret = TFile::Open(gaussDir + gaussName);
65 TH2*
temp = (TH2*)h->Clone();
67 int NX = h->GetXaxis()->GetNbins();
68 int NY = h->GetYaxis()->GetNbins();
69 double minX = h->GetXaxis()->GetBinLowEdge(1);
70 double maxX = h->GetXaxis()->GetBinUpEdge(NX);
71 double widthX = h->GetXaxis()->GetBinWidth(1);
72 double minY = h->GetYaxis()->GetBinLowEdge(1);
73 double maxY = h->GetYaxis()->GetBinUpEdge(NY);
74 double widthY = h->GetYaxis()->GetBinWidth(1);
81 int foldedBin = 1 + NX/2;
84 for(
int i = 1;
i <= NX/2;
i++) {
85 for(
int j = 1;
j <= NY/2;
j++) {
89 int x2 =
i%(NX/2+1)+2;
92 double z1 = h->GetBinContent(x1, y1);
93 double z2 = h->GetBinContent(x2, y2);
94 temp->SetBinContent(x1, y1, z2);
95 temp->SetBinContent(x2, y2, z1);
103 for(
int y = 1;
y <= NY;
y++) {
104 double sZ1 = temp->GetBinContent(foldedBin,
y);
105 double sZ2 = temp->GetBinContent(foldedBin,
y);
106 h->SetBinContent(1,
y, sZ1);
107 h->SetBinContent(NX,
y, sZ2);
114 if(beam ==
"rhc") isFHC =
false;
118 if(!NERSC)
return "";
123 +
"/Predictions/provisional_singles";
124 if(isFHC && !isFakeData)
125 ret = dir +
"/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_FHCAllSyst_nueconcat_realNDData_2018-05-22.root";
126 else if(!isFHC && !isFakeData)
127 ret = dir +
"/NueProd4SystematicsOnRealNDData_2018-05-22/pred_allxp_Nue2018Axis_full_RHCAllSyst_nueconcat_realNDData_2018-05-22.root";
128 else if(isFHC && isFakeData)
129 ret = dir +
"/NueProd4SystematicsOnFakeNDData_2018-05-03/pred_allxp_Nue2018Axis_full_FHCAllSysts_nueconcat_fakeNDData_05_03_2018_merged.root";
130 else if(!isFHC & isFakeData)
131 ret = dir +
"/NueProd4SystematicsOnFakeNDData_2018-05-03/pred_allxp_Nue2018Axis_full_RHCAllSysts_nueconcat_fakeNDData_05_03_2018_merged.root";
141 bool GetFromUPS=
false,
142 bool minimizeMemory =
false,
146 if(decomp ==
"noextrap") extrap =
kNoExtrap;
148 if(decomp ==
"combo") extrap =
kCombo;
149 if(decomp ==
"fakexp") extrap =
kFake;
156 std::cout <<
"Getting Nue Prediction from " + nue_pred_path <<
endl;
159 beam, isFakeData,
true,
160 nue_pred_path, minimizeMemory);
163 beam, isFakeData,
true,
164 nue_pred_path, minimizeMemory);
168 true, nue_pred_path, minimizeMemory);
177 beam, isFakeData,
true,
178 upsname+name, minimizeMemory);
181 beam, isFakeData,
true,
182 upsname+name, minimizeMemory);
187 true, upsname+
name, minimizeMemory);
201 double nonsScale = 1./10.45;
202 double swapScale = 1./12.91;
205 if(NERSC) dirName =
std::string(
getenv(
"FCHELPERANA2018_LIB_PATH"))+
"/Predictions/FDRock/";
209 if(beam==
"rhc") rName =
"pred_FDRock_merg_rhc";
210 else rName =
"pred_FDRock_merg_fhc";
211 if( decomp ==
"noextrap"){
212 if(beam==
"rhc") rName =
"pred_FDRock_rhc";
213 else rName =
"pred_FDRock_fhc";
218 fName2 =
"/fdrock_nue2017.root";
219 rName =
"pred_FDRock";
222 auto rock = LoadFromFile<PredictionNoExtrap>(dirName + fName2, rName).
release();
224 std::cout <<
"Nue prediction (rock only - not scaled) " 229 std::cout <<
"Nue prediction (rock added) " 230 << nue_pred->Predict(calc).Integral(pot) <<
std::endl;
250 std::string anaDir =
"/nova/ana/nu_e_ana/Ana2018/";
252 if(beam==
"rhc") name =
"cosmic_spect_rhc";
253 else name=
"cosmic_spect_fhc";
254 outtime = LoadFromFile<Spectrum>
255 (anaDir+
"/Predictions/cosmic/cosmic_prediction_real_data.root",
name ).
release();
266 <<
" cosmics from out-of-time." <<
std::endl;
285 intime = LoadFromFile<Spectrum>
286 (
"/nova/ana/nu_e_ana/Ana2018/Data2018/spectra.root",
"data_spect_fhc").
release();
288 else if(beam==
"rhc"){
289 intime = LoadFromFile<Spectrum>
290 (
"/nova/ana/nu_e_ana/Ana2018/Data2018/spectra.root",
"data_spect_rhc").
release();
298 <<
"POT " << intime->
POT()
299 <<
" (expected " << POT <<
")\n" 300 <<
"Livetime " << intime->
Livetime()
301 <<
" (expected " << livetime <<
")\n" 311 bool useSysts =
true,
313 bool GetFromUPS=
false,
315 bool minimizeMemory =
false,
320 std::vector <const IPrediction*> numu_preds;
323 std::string anaDir =
"/nova/ana/nu_mu_ana/Ana2018/";
325 if(
beam==
"rhc") filename=anaDir +
"/Predictions/pred_numuconcat_rhc__numu2018.root";
326 else filename= anaDir +
"/Predictions/pred_numuconcat_fhc__numu2018.root";
328 std::cout <<
"\nLoading Numu Predictions\n" ;
329 for (
int quant = 0; quant < nq; ++quant){
356 filename, minimizeMemory));
365 bool GetFromUPS=
false,
bool NERSC =
false)
373 fcosm = TFile::Open((upsname+
"/cos_bkg_hists_v2.root").c_str());
379 if (
beam==
"rhc") filename=
"cosmics_rhc__numu2018.root";
380 else filename=
"cosmics_fhc__numu2018.root";
381 fcosm = TFile::Open((dir+filename).c_str());
384 if(fcosm->IsZombie()) {
std::cerr<<
"bad cosmics\n";
exit(1);}
385 std::vector <std::pair <TH1D*, double > > numu_cosmics;
387 for (
int i = 1;
i <=nq; ++
i) {
392 numu_cosmics.push_back({
h,1/
sqrt(
h->Integral())});
393 std::cout <<
"Cosmics all periods quantile " <<
i <<
" " 397 << 1/
sqrt(
h->Integral())
407 if(
beam ==
"fhc") {file=
"/nova/ana/nu_mu_ana/Ana2018/Data/fd_dataspectrum_fhc_full__numu2018.root";}
408 if(
beam ==
"rhc") {file=
"/nova/ana/nu_mu_ana/Ana2018/Data/fd_dataspectrum_rhc_full__numu2018.root";}
410 std::vector <Spectrum *> numu_data;
412 for (
int i = 1;
i <=nq; ++
i) {
414 auto f=
new TFile(file.c_str());
443 double dmsq32 = -5,
double dCP_Pi = -5)
452 std::vector<const ISyst*>
GetJointFitSystematicList(
bool corrSysts,
bool nueExclusive =
false,
bool numuExclusive=
false,
bool isFHC =
true,
bool isRHC =
true ,
bool intersection =
true){
453 std::vector<const ISyst*>
ret ={} ;
484 std::cout <<
"Return list of systematics: " ;
485 for(
auto &
s:ret)
std::cout <<
s->ShortName() <<
", ";
490 std::vector<std::pair<const ISyst*,const ISyst*> >
ret ={};
491 std::cout<<(isNue?
"not nue ":
"not numu ")<<(isFHC&&!isNue?
" ":(isFHC?
" FHC ":
" RHC "));
495 for (
auto &
s:badsyst) ret.push_back ({
s,NULL});
540 c1->SetLeftMargin(0.14);
541 c1->SetBottomMargin(0.15);
543 if(slicename.Contains(
"delta")) {
544 title =
";delta;Significance (#sqrt{#Delta#chi^{2}})";
546 if(slicename.Contains(
"th23")) {
547 title =
";sin^{2}#theta_{23};Significance (#sqrt{#Delta#chi^{2}})";
549 if(slicename.Contains(
"dmsq")) {
550 title =
";|#Deltam^{2}_{32}| (10^{-3} eV^{2});Significance (#sqrt{#Delta#chi^{2}})";
552 if(slicename.Contains(
"th13")) {
553 title =
";sin^{2}2#theta_{13};Significance (#sqrt{#Delta#chi^{2}})";
555 auto axes =
new TH2F(
"",title, 100,
xmin,
xmax, 100, 0, ymax);
556 if(slicename.Contains(
"fccorr") ||
overlay)
axes->GetYaxis()->SetTitle(
"Significance (#sigma)");
558 axes->GetXaxis()->CenterTitle();
559 axes->GetYaxis()->CenterTitle();
565 axes->GetXaxis()->SetTitleOffset(0.8);
566 axes->GetYaxis()->SetTitleOffset(0.75);
571 TLegend *
SliceLegend(std::vector <std::pair <TGraph*, TString > > &
graphs,
bool isDelta,
double low=0.7,
bool isDmsq =
false,
bool overlay =
false)
574 if(isDelta) leg =
new TLegend(0.6, 0.6, 0.9,0.89);
576 leg =
new TLegend(0.18, 0.2, 0.4, 0.49);
579 leg->SetFillStyle(0);
582 if(gr.second.Contains(
"O")){
583 if(isDelta) leg->SetX1(0.57);
584 if(isDmsq) leg->SetX1(0.18);
585 if(!isDelta && !isDmsq) {leg->SetX1(0.6); leg->SetX2(0.89);}
586 leg->SetMargin(0.18);
587 str += gr.second.Contains(
"NH") ?
"NH":
"IH";
588 str += gr.second.Contains(
"LO") ?
" Lower":
" Upper";
593 str += gr.second.Contains(
"NH") ?
"Normal":
"Inverted";
594 str +=
"}{hierarchy}";
597 str += gr.second.Contains(
"NH") ?
"NH" :
"IH";
600 leg->AddEntry(gr.first, str,
"l");
607 double maxx = 2,
double miny = 0,
double maxy = 1){
610 c1->SetLeftMargin(0.14);
611 c1->SetBottomMargin(0.15);
612 TH2*
axes =
new TH2F();
614 if(surfName.Contains(
"delta_th23")) title=
";#delta_{CP};sin^{2}#theta_{23}";
615 if(surfName.Contains(
"th13_delta")) title=
";sin^{2}2#theta_{13};#delta_{CP}/#pi";
616 if(surfName.Contains(
"th23_dm32")) title=
";sin^{2}#theta_{23};#Deltam^{2}_{32} (10^{-3}eV^{2})";
617 axes =
new TH2F(
"",title,100,minx,maxx,100,miny,
maxy);
619 axes->GetXaxis()->CenterTitle();
620 axes->GetYaxis()->CenterTitle();
626 axes->GetXaxis()->SetTitleOffset(0.8);
627 axes->GetYaxis()->SetTitleOffset(0.8);
628 TGaxis::SetMaxDigits(3);
629 if(surfName.Contains(
"th23_dm32")) {
630 axes->GetYaxis()->SetDecimals() ;
631 axes->GetYaxis()->SetTitleOffset(0.85);
632 axes->GetYaxis()->SetLabelOffset(0.002);
633 axes->GetYaxis()->SetLabelSize(0.058);
634 axes->GetYaxis()->SetTitleSize(0.078);
636 axes->GetXaxis()->SetTitleOffset(0.86);
644 Int_t kFillColor1,Int_t kFillColor2, Int_t kFillColor3,Int_t kDarkColor,
646 TLegend *
leg =
new TLegend(0.16,fccorr?0.19:0.17,0.65,0.28);
648 leg->SetFillStyle(0);
649 leg->SetMargin(1.3*leg->GetMargin());
652 TGraph *
dummy =
new TGraph;
653 dummy->SetLineWidth(3);
656 dummy->SetLineColor(kDarkColor);
657 dummy->SetFillColor(kFillColor1);
658 leg->AddEntry(dummy->Clone(),
"1#sigma",
"f");
659 dummy->SetFillColor(kFillColor2);
660 leg->AddEntry(dummy->Clone(),
"2#sigma",
"f");
661 dummy->SetFillColor(kFillColor3);
662 leg->AddEntry(dummy->Clone(),
"3#sigma",
"f");
665 leg->SetMargin(1.4*leg->GetMargin());
666 dummy->SetLineColor(kFillColor1);
667 dummy->SetLineStyle(kSolid);
668 leg->AddEntry(dummy->Clone(),
"1#sigma",
"l");
670 dummy->SetLineColor(kFillColor2);
671 leg->AddEntry(dummy->Clone(),
"2#sigma",
"l");
672 dummy->SetLineStyle(kSolid);
673 dummy->SetLineColor(kFillColor3);
674 leg->AddEntry(dummy->Clone(),
"3#sigma",
"l");
681 dummy->SetMarkerStyle(43);
682 dummy->SetMarkerSize(2);
683 dummy->SetMarkerColor(kDarkColor);
684 dummy->SetLineColor(kDarkColor);
685 dummy->SetLineStyle(7);
686 if(hie>0) leg->AddEntry(dummy->Clone(),
"Best Fit",
"p");
687 if(hie<0) leg->AddEntry((TObject*)0,
"#color[0]{Best Fit}",
"");
688 dummy->SetLineStyle(1);
690 if(!fccorr) leg->SetHeader(
"No Feldman-Cousins");
691 else leg->SetHeader(
"Feldman-Cousins");
697 Int_t kFillColor1,Int_t kFillColor2, Int_t kFillColor3,Int_t kDarkColor)
699 TLegend *
leg =
new TLegend(0.16,hie>0?.78:0.77,0.65,0.86);
701 leg->SetFillStyle(0);
702 leg->SetMargin(1.3*leg->GetMargin());
705 TGraph *
dummy =
new TGraph;
706 dummy->SetLineWidth(3);
709 dummy->SetLineColor(kDarkColor);
710 dummy->SetFillColor(kFillColor1);
711 leg->AddEntry(dummy->Clone(),
"1#sigma",
"f");
712 dummy->SetFillColor(kFillColor2);
713 leg->AddEntry(dummy->Clone(),
"2#sigma",
"f");
714 dummy->SetFillColor(kFillColor3);
715 leg->AddEntry(dummy->Clone(),
"3#sigma",
"f");
718 leg->SetMargin(1.4*leg->GetMargin());
719 dummy->SetLineColor(kFillColor1);
720 dummy->SetLineStyle(kSolid);
721 leg->AddEntry(dummy->Clone(),
"1#sigma",
"l");
723 dummy->SetLineColor(kFillColor2);
724 leg->AddEntry(dummy->Clone(),
"2#sigma",
"l");
725 dummy->SetLineStyle(kSolid);
726 dummy->SetLineColor(kFillColor3);
727 leg->AddEntry(dummy->Clone(),
"3#sigma",
"l");
730 leg->SetHeader(
"No Feldman-Cousins");
TCut intime("tave>=217.0 && tave<=227.8")
const NOvARwgtSyst k2ndClassCurrs("2ndclasscurr","Second class currents", novarwgt::kSimpleSecondClassCurrentsSystKnob)
Second-class current syst. See documentation in NOvARwgt.
Cuts and Vars for the 2020 FD DiF Study.
std::map< std::string, double > xmax
Float_t y1[n_points_granero]
TH1D * ToTH1(double exposure, Color_t col=kBlack, Style_t style=kSolid, EExposureType expotype=kPOT, EBinType bintype=kBinContent) const
Histogram made from this Spectrum, scaled to some exposure.
Loads shifted spectra from files.
const FitDmSq32 kFitDmSq32
virtual Spectrum Predict(osc::IOscCalc *calc) const =0
Float_t x1[n_points_granero]
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
void XAxisDeltaCPLabels(TH1 *axes, bool t2kunits)
Label the x-axis with fractions of pi.
const FitSinSqTheta23UpperOctant kFitSinSqTheta23UpperOctant
General interface to oscillation calculators.
virtual double HighLimit() const
virtual T GetTh23() const
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var's SetValue()
void AddJointAna2018NotCorrelSysts(std::vector< const ISyst * > &systs, const EAnaType2018 ana, const BeamType2018 beam)
double Integral(double exposure, double *err=0, EExposureType expotype=kPOT) const
Return total number of events scaled to pot.
virtual double HighLimit() const
virtual double LowLimit() const
osc::IOscCalcAdjustable * DefaultOscCalc()
Create a new calculator with default assumptions for all parameters.
T sqr(T x)
More efficient square function than pow(x,2)
Spectrum MockData(double pot, int seed=0) const
Mock data is FakeData with Poisson fluctuations applied.
Representation of a spectrum in any variable, with associated POT.
const double kAna2018RHCPOT
const NOvARwgtSyst kRadCorrNuebar("radcorrnuebar","Radiative corrections for #bar{#nu}_{e}", novarwgt::kSimpleRadiativeCorrNuebarXsecSystKnob)
Radiative corrections syst (nuebars). See documentation in NOvARwgt.
std::vector< const ISyst * > getAllAna2018Systs(const EAnaType2018 ana, const bool smallgenie, const BeamType2018 beam, bool isFit)
static std::unique_ptr< Spectrum > LoadFrom(TDirectory *dir, const std::string &name)
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::string getenv(std::string const &name)
std::vector< float > Spectrum
Oscillation probability calculators.
std::vector< double > POT
virtual std::string LatexName() const
Spectrum Predict(osc::IOscCalc *calc) const override
const Style_t k2SigmaConfidenceStyle
const DummyAnaSyst kAna2018NormRHC("NormRHC2018","RHC. Norm.")
virtual double LowLimit() const
const FitSinSqTheta23LowerOctant kFitSinSqTheta23LowerOctant
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
void AddNonLoadable2018Systs(std::vector< const ISyst * > &systs, const EAnaType2018 ana)
TCut outtime("(tave>0.0&&tave<200.0) || (tave>250.0&&tave<450.0)")
std::vector< std::pair< const ISyst *, const ISyst * > > GetCorrelations(bool isNue, bool isFHC, bool ptExtrap)
Standard interface to all prediction techniques.
const Float_t kBlessedLabelFontSize
const double kAna2018FHCPOT
std::string to_string(ModuleType mt)
virtual std::string ShortName() const
const Float_t kBlessedTitleFontSize
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)
const NOvARwgtSyst kRadCorrNue("radcorrnue","Radiative corrections for #nu_{e}", novarwgt::kSimpleRadiativeCorrNueXsecSystKnob)
Radiative corrections syst (nues). See documentation in NOvARwgt.
double Livetime() const
Seconds. For informational purposes only. No calculations use this.
std::string UniqueName()
Return a different string each time, for creating histograms.
void rock(std::string suffix="full")
const DummyAnaSyst kAna2018NormFHC("NormFHC2018","FHC. Norm.")
const double kAna2018FHCLivetime
virtual std::string LatexName() const
virtual double GetValue(const osc::IOscCalcAdjustable *osc) const
const double kAna2018RHCLivetime
virtual void SetValue(osc::IOscCalcAdjustable *osc, double val) const
virtual std::string ShortName() const
FitSinSqTheta23LowerOctant()