24 #include "Utilities/func/MathUtil.h" 54 std::vector <std::pair <int, int> > potCombos,
55 Int_t
color, Int_t style);
67 bool stats = (systType==
"_stats");
68 double ssth23 = ((fakeOctant==
"UO")?0.625:0.403);
69 if(fakeOctant ==
"MAX") ssth23 = 0.5;
73 TString
filename =
"root_reach/hist_novareach_"+ suffix;
74 TString filenameCompare = (
"root_reach/hist_novareach_"+
MakeSuffix(
"CVN",fakeHieDelta,ssth23,potCompID,systType,
version));
75 TString plotfolder =
version+
"/";
76 gSystem->mkdir(
"root_reach");
77 gSystem->mkdir(plotfolder);
80 std::string hieDelStr =
"Normal #delta_{CP}=3#pi/2";
81 if(fakeHieDelta ==
"IH3pi2") hieDelStr =
"Inverted #delta_{CP}=3#pi/2";
83 TString plottitle =
TString::Format(
"%s, sin^{2}#theta_{23}=%.3f, #Deltam^{2}_{32}=2.5#times10^{-3}eV^{2}, sin^{2}#theta_{13}=0.022",hieDelStr.c_str(),
ssth23);
84 TString plottitle_split =
TString::Format(
"#splitline{%s, sin^{2}#theta_{23}=%.3f}{#Deltam^{2}_{32}=2.5#times10^{-3}eV^{2}, sin^{2}#theta_{13}=0.022}",hieDelStr.c_str(),
ssth23);
86 std::vector <int> potYears = {2016,2017, 2018, 2019, 2020, 2021, 2022, 2023, 2024};
87 std::vector <std::pair <int, int> > potComboT = {{6,0},{9,3}};
88 std::vector <std::pair <int, int> > potComboA = {{6,0},{9,3},{9,9},{12,12},{15,15},{18,18},{21,21},{24,24},{27,27}};
89 std::vector <std::pair <int, int> > potComboB = {{6,0},{12,0},{18,0},{24,0},{30,0},{36,0},{42,0},{48,0},{54,0}};
90 std::vector <std::pair <int, int> > potComboC = {{6,0},{9,3},{9,9},{15,9},{18,12},{21,15},{24,18},{27,21},{30,24}};
91 std::vector <std::pair <int, int> > potComboD = {{6,0},{9,3},{9,9},{13,12},{17,15},{21,18},{25,21},{29,24},{33,27}};
92 std::vector <std::pair <int, int> > potComboE = {{6,0},{9,3},{9,9},{13,12},{16,16},{20,19},{23,23},{27,26},{30,30}};
94 std::vector <std::pair <int, int> > potComboBase;
95 std::vector <std::pair <int, int> > potComboCompare;
97 if(potBaseID ==
"comboT") potComboBase = potComboT;
98 if(potBaseID ==
"comboA") potComboBase = potComboA;
99 if(potBaseID ==
"comboB") potComboBase = potComboB;
100 if(potBaseID ==
"comboC") potComboBase = potComboC;
101 if(potBaseID ==
"comboD") potComboBase = potComboD;
102 if(potBaseID ==
"comboE") potComboBase = potComboE;
104 if(potCompID ==
"comboA") potComboCompare = potComboA;
105 if(potCompID ==
"comboB") potComboCompare = potComboB;
106 if(potCompID ==
"comboC") potComboCompare = potComboC;
107 if(potCompID ==
"comboD") potComboCompare = potComboD;
108 if(potCompID ==
"comboE") potComboCompare = potComboE;
110 bool compare = (potCompID!=
"none");
113 TFile *
file =
new TFile (filename+
".root",
"recreate");
116 TH2D * sign_Hie =
new TH2D (
"hie",
"hie; Total POT (#times10^{20}); RHC POT (#times10^{20})",
117 67,5.5,72.5,67,-0.5,66.5);
118 TH2D * sign_Max =
new TH2D (
"max",
"max; Total POT (#times10^{20}); RHC POT (#times10^{20})",
119 67,5.5,72.5,67,-0.5,66.5);
120 TH2D * sign_Oct =
new TH2D (
"oct",
"oct; Total POT (#times10^{20}); RHC POT (#times10^{20})",
121 67,5.5,72.5,67,-0.5,66.5);
122 TH2D * sign_Cpv =
new TH2D (
"cpv",
"cpv; Total POT (#times10^{20}); RHC POT (#times10^{20})",
123 67,5.5,72.5,67,-0.5,66.5);
125 TString
dir =
"root_predictions/";
126 TFile * file1 =
new TFile( dir+
"predInterp_nue_FHC1_"+
version+
".root",
"read");
127 TFile * file2 =
new TFile( dir+
"predInterp_nue_FHC2_"+
version+
".root",
"read");
128 TFile * file3 =
new TFile( dir+
"predInterp_nue_RHC_"+
version+
".root",
"read");
129 TFile * file4 =
new TFile( dir+
"predInterp_numu_FHC1_full_"+
version+
".root",
"read");
130 TFile * file5 =
new TFile( dir+
"predInterp_numu_FHC2_full_"+
version+
".root",
"read");
131 TFile * file6 =
new TFile( dir+
"predInterp_numu_RHC_full_"+
version+
".root",
"read");
147 delete file1, file2, file3, file4, file5, file6;
153 int totpot = fhc+rhc;
155 std::cerr <<
"\nFinding best fit points for " << fhc <<
"p" << rhc
156 <<
"ssth23 " << ssth23 <<
"\n\n";
159 double p2POT = fhc*1E20 - 3.41E20;
160 double potFHC = p1POT + p2POT;
161 double potRHC = rhc*1E20;
169 if(fakeHieDelta ==
"IH3pi2") calc->SetDmsq32(-1*
fabs(calc->GetDmsq32()));
171 auto obsFHC = predFHC->Predict(calc).FakeData(potFHC);
172 auto obsRHC = predRHC->Predict(calc).FakeData(potRHC);
174 auto ubsRHC = prudRHC->Predict(calc).
FakeData(potRHC);
183 if(rhc > 0 ) expt =
new MultiExperiment({exptFHC, uxptFHC, exptRHC, uxptRHC,
193 std::vector< std::pair< const ISyst *, const ISyst * >> corr_nue_FHC;
194 std::vector< std::pair< const ISyst *, const ISyst * >> corr_numu_FHC;
195 std::vector< std::pair< const ISyst *, const ISyst * >> corr_nue_RHC;
196 std::vector< std::pair< const ISyst *, const ISyst * >> corr_numu_RHC;
200 for (
auto nosys:
systsNumuRHC) corr_nue_FHC.push_back({nosys,NULL});
203 for (
auto nosys:
systsNueRHC) corr_numu_FHC.push_back({nosys,NULL});
207 for (
auto nosys:
systsNumuFHC) corr_nue_RHC.push_back({nosys,NULL});
210 for (
auto nosys:
systsNueFHC) corr_numu_RHC.push_back({nosys,NULL});
214 if(stats) { mySysts.clear();
std::cout <<
"WARNING STATS ONLY ACTIVATED \n"; }
223 double minchi_NH_UO=1E20;
224 double minchi_NH_LO=1E20;
225 double minchi_IH_UO=1E20;
226 double minchi_IH_LO=1E20;
228 for(
int hie:{-1,+1}){
229 for(
double seed:{0.3,0.7}){
230 double minchi_temp = 1E20;
233 calc->SetDmsq32(hie*
fabs(calc->GetDmsq32()));
234 minchi_temp =
std::min(minchi_temp, fit23.Fit(calc)->EvalMetricVal());
235 const double dcp = calc->GetdCP();
236 for(
int n = 1;
n <= 3; ++
n){
237 calc->SetdCP(dcp+
n*
M_PI/2);
238 minchi_temp =
std::min(minchi_temp, fit23.Fit(calc)->EvalMetricVal());
240 if(hie>0 &&
seed<0.5) minchi_NH_LO = minchi_temp;
241 if(hie>0 &&
seed>0.5) minchi_NH_UO = minchi_temp;
242 if(hie<0 &&
seed<0.5) minchi_IH_LO = minchi_temp;
243 if(hie<0 && seed>0.5) minchi_IH_UO = minchi_temp;
245 std::cout <<
"\n \n Found minchi " << (hie>0?
"NH ":
"IH ")
246 << (
seed<0.5?
"LO ":
"UO ") << minchi_temp <<
"\n\n";
250 double minchi_NH =
std::min(minchi_NH_UO,minchi_NH_LO);
251 double minchi_IH =
std::min(minchi_IH_UO,minchi_IH_LO);
252 double minchi_UO =
std::min(minchi_NH_UO,minchi_IH_UO);
253 double minchi_LO =
std::min(minchi_NH_LO,minchi_IH_LO);
254 double minchi_best =
std::min(minchi_NH,minchi_IH);
258 double minchi_Max=1E20;
260 for(
int hie:{-1,+1}){
263 calc->SetDmsq32(hie*
fabs(calc->GetDmsq32()));
264 minchi_Max =
std::min(minchi_Max,fit23Max.Fit(calc)->EvalMetricVal());
265 const double dcp2 = calc->GetdCP();
266 for(
int n = 1;
n <= 3; ++
n){
267 calc->SetdCP(dcp2+
n*
M_PI/2);
268 minchi_Max =
std::min(minchi_Max, fit23Max.Fit(calc)->EvalMetricVal());
273 double minchi_CPV=1E20;
275 for(
int hie:{-1,+1}){
276 for(
double seed:{0.3,0.7}){
277 for(
double dCP:{0,1}){
281 calc->SetDmsq32(hie*
fabs(calc->GetDmsq32()));
282 minchi_CPV =
min(minchi_CPV, fit23CP.Fit(calc)->EvalMetricVal());
287 std::cout <<
"\n \n Found minchi NH " << minchi_NH
288 <<
" IH " << minchi_IH
289 <<
" UO " << minchi_UO
290 <<
" LO " << minchi_LO
291 <<
" Max " << minchi_Max
292 <<
" CPV " << minchi_CPV
295 sign_Hie->SetBinContent(sign_Hie->FindBin(totpot,rhc),
fabs(minchi_IH - minchi_NH));
296 sign_Oct->SetBinContent(sign_Oct->FindBin(totpot,rhc),
fabs(minchi_LO - minchi_UO));
297 sign_Max->SetBinContent(sign_Max->FindBin(totpot,rhc),minchi_Max - minchi_best);
298 sign_Cpv->SetBinContent(sign_Cpv->FindBin(totpot,rhc),minchi_CPV - minchi_best);
301 sign_Hie->SetTitle(
"#Delta #chi^{2} Wrong hie. rejection - "+ plottitle);
302 sign_Oct->SetTitle(
"#Delta #chi^{2} Octant resolution - "+plottitle);
303 sign_Max->SetTitle(
"#Delta #chi^{2} Maximal mix. rejection - "+plottitle);
304 sign_Cpv->SetTitle(
"#Delta #chi^{2} CP violation - "+plottitle);
306 sign_Hie->SetOption(
"colz");
307 sign_Oct->SetOption(
"colz");
308 sign_Max->SetOption(
"colz");
309 sign_Cpv->SetOption(
"colz");
323 std::vector <TString> sign_label={
"max",
"hie",
"oct",
"cpv"};
324 std::vector <TString> sign_name={
"Max. mixing",
"Hierarchy",
"Octant",
"CPV"};
327 TColor::GetColor(
"#6c71c4"), TColor::GetColor(
"#268bd2"),
328 TColor::GetColor(
"#2aa198"), TColor::GetColor(
"#859900"),
329 TColor::GetColor(
"#b58900"), TColor::GetColor(
"#cb4b16"),
330 TColor::GetColor(
"#dc322f"), TColor::GetColor(
"#d33682")
338 TLegend *
leg =
new TLegend (0.12, 0.52, 0.47, 0.85);
339 leg->SetTextSize(0.055);
340 leg->SetFillStyle(0);
341 leg->SetMargin(1.3*leg->GetMargin());
342 leg->SetHeader(
"NOvA joint #nu_{e}+#nu_{#mu}");
345 TFile *
file =
new TFile (filename+
".root",
"read");
347 for(
auto sign:sign_label){
349 auto gr =
GetSignGraph(
hist, potYears,potComboBase,potColors[sigIdx*2],lineStyle[sigIdx]);
351 leg->AddEntry(gr,sign_name[sigIdx],
"l");
356 TFile * fileCompare =
new TFile (filenameCompare+
".root",
"read");
357 if(fileCompare->IsZombie()) {
std::cout <<
"can t find " << fileCompare->GetName()
360 for(TString
sign:{
"hie",
"oct",
"max",
"cpv"}){
361 auto hist= (TH1D*) fileCompare->Get(
sign);
362 Int_t trans_color = TColor::GetColorTransparent(potColors[sigIdx*2], 0.3);
370 if(!compare) gPad->Print(plotfolder+
"sign_reach_"+suffix+
".pdf");
371 else gPad->Print(plotfolder+
"sign_reach_"+suffix+
"vs"+potCompID+
".pdf");
406 std::vector <std::pair <int, int> > potCombos,
407 Int_t
color, Int_t style){
408 TGraph * gr =
new TGraph();
413 int totpot = fhc+rhc;
414 double this_chisq = hist->GetBinContent(hist->FindBin(totpot,rhc));
416 gr->SetPoint(gr->GetN(),potYears[yearIdx], this_sign);
420 gr->SetLineColor(color);
421 gr->SetLineStyle(style);
427 gPad->SetTopMargin(0.14);
428 TH2D *
axes =
new TH2D (
"ax",
";Year;Significance (#sigma)",18,2015.5,2024.5,100,0,5.9);
429 axes->GetXaxis()->CenterTitle();
430 axes->GetYaxis()->CenterTitle();
431 axes->GetYaxis()->SetTitleOffset(0.6);
433 TLatex * ltx =
new TLatex();
434 ltx->SetTextAlign(13);
435 ltx->SetTextSize(axes->GetXaxis()->GetLabelSize());
436 ltx->DrawLatexNDC(0.12,0.99,plottitle_split.Data() );
438 TLatex* ltx2 =
new TLatex();
439 ltx2->SetTextSize(1.2*axes->GetXaxis()->GetLabelSize());
440 ltx2->DrawLatexNDC(0.35,0.16,
"#splitline{2016 analysis techniques with projected}{systematic uncertainty improvements}");
451 suffix+=
"_"+fakeHieDelta;
452 if(ssth23<1) suffix+=(
std::to_string(ssth23)).replace(0,2,
"_0p").erase(6,4);
453 else suffix+=
"_0pXX";
454 bool stats = (systType==
"_stats");
455 suffix+=(
"_"+ potID+(stats?
"":
"_syst")+systType);
T max(const caf::Proxy< T > &a, T b)
virtual void SetL(double L)=0
Cuts and Vars for the 2020 FD DiF Study.
_OscCalcPMNSOpt< double > OscCalcPMNSOpt
void futureSig_reach_singlePOTcombo_syst(const bool makeFile=false, const std::string version="v01", const std::string fakeHieDelta="NH3pi2", const std::string fakeOctant="UO", const std::string systType="_12", const std::string potBaseID="comboA", const std::string potCompID="none")
fvar< T > fabs(const fvar< T > &x)
const double kSecondAnaPeriod2POT
virtual void SetDmsq21(const T &dmsq21)=0
const FitDmSq32 kFitDmSq32
const ReactorExperiment * WorldReactorConstraint2015()
Weighted average of all experiments as of first nue paper writing.
std::vector< const ISyst * > systsNumuFHC
const double kSecondAnaEpoch3bPOT
osc::IOscCalcAdjustable * MyDefaultOscCalc()
virtual void SetTh13(const T &th13)=0
void PaintReachCanvas(TString plottitle_split)
void ResetOscCalcToMyDefault(osc::IOscCalcAdjustable *calc)
int stats(TString inputFilePath, Int_t firstRun, Int_t lastRun, Float_t thresh, TString myDet)
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
Forward to wrapped Var's SetValue()
const double kSecondAnaPeriod1POT
virtual void SetDmsq32(const T &dmsq32)=0
std::vector< const ISyst * > GetSystsList(TString systType)
make combinations of the systematics
TGraph * GetSignGraph(TH1D *hist, std::vector< int > potYears, std::vector< std::pair< int, int > > potCombos, Int_t color, Int_t style)
Spectrum FakeData(double pot) const
Synonymous with AsimovData(). Retained for compatibility.
std::vector< const ISyst * > systsNueFHC
void SetValue(osc::IOscCalcAdjustable *osc, double val) const override
std::vector< const ISyst * > systsNumuRHC
std::vector< const ISyst * > systsNueRHC
static float min(const float a, const float b, const float c)
Combine multiple component experiments.
virtual void SetRho(double rho)=0
static std::unique_ptr< PredictionInterp > LoadFrom(TDirectory *dir, const std::string &name)
const ConstrainedFitVarWithPrior fitSsqTh23_UniformPriorSsqTh23 & kFitSinSqTheta23
virtual void SetTh23(const T &th23)=0
bool compare(const GFluxGenerator &g1, const GFluxGenerator &g2)
General interface to any calculator that lets you set the parameters.
Sum MC predictions from different periods scaled according to data POT targets.
void Format(TGraph *gr, int lcol, int lsty, int lwid, int mcol, int msty, double msiz)
const FitSinSq2Theta13 kFitSinSq2Theta13
std::string to_string(ModuleType mt)
T min(const caf::Proxy< T > &a, T b)
const FitVarWithPrior fitDeltaInPiUnits_UniformPriordCP & kFitDeltaInPiUnits
std::string MakeSuffix(const std::string PID, const std::string fakeHieDelta, const double ssth23, const std::string potID, const std::string systType, const std::string version)
void SetSystCorrelations(int idx, const std::vector< std::pair< const ISyst *, const ISyst * >> &corrs)
virtual void SetTh12(const T &th12)=0
virtual void SetdCP(const T &dCP)=0
Spectrum Predict(osc::IOscCalc *calc) const override
Compare a single data spectrum to the MC + cosmics expectation.
Perform MINUIT fits in one or two dimensions.