47 double xq[1] = {percentile};
49 hist->GetQuantiles(1, yq, xq);
55 std::ifstream
infile(fileName);
61 bool corrSysts =
false,
63 TString
options=
"joint_realData_both_onlyIH",
64 bool dmsqSurf =
false,
65 bool th13Surf =
false,
67 bool fillContour =
false)
71 bool nueOnly =
options.Contains(
"nueOnly");
72 bool numuOnly =
options.Contains(
"numuOnly");
73 bool joint =
options.Contains(
"joint");
74 assert (nueOnly || numuOnly || joint);
76 bool FHCOnly =
options.Contains(
"FHCOnly");
77 bool RHCOnly =
options.Contains(
"RHCOnly");
78 bool both =
options.Contains(
"both");
79 assert (FHCOnly || RHCOnly || both);
81 bool fake2017 =
options.Contains(
"fake2017");
82 bool realData =
options.Contains(
"realData");
84 bool onlyNH =
options.Contains(
"onlyNH");
85 bool onlyIH =
options.Contains(
"onlyIH");
88 if(dmsqSurf) suffix+=
"_dmsq";
89 if(th13Surf) suffix+=
"_th13";
90 if(corrSysts) suffix+=
"_systs";
91 else suffix+=
"_stats10";
92 string syst = (corrSysts?
"syst":
"stat");
101 string outdir2018 =
"";
102 if(th13Surf) outdir2018 =
"th13_delta/";
103 else if (dmsqSurf) outdir2018 =
"th23_dmsq/";
104 else outdir2018 =
"delta_th23/";
105 TString ana2018 (
"/nova/ana/nu_e_ana/Ana2018/Results/RHC_and_FHC/contours/" + outdir2018 +
"stat/hist_contours_2018_" +
options +
"combo_"+(dmsqSurf?
"dmsq_":
"")+
"stats");
106 TFile * infile2018ana =
new TFile (ana2018 +
".root",
"read");
107 auto mins2018 =* (
TVectorD*)infile2018ana->Get(
"overall_minchi");
108 double minchi23_2018ana = mins2018[0];
110 TString
ana2019 (
"/nova/ana/nu_e_ana/Ana2019/Results/RHC_and_FHC/contours/" + outdir2018 +
"stat/hist_contours_2019_" +
options +
"_"+(dmsqSurf?
"dmsq_":
"")+
"stats");
111 TFile * infile2019ana =
new TFile (ana2019 +
".root",
"read");
112 auto mins2019 =* (
TVectorD*)infile2019ana->Get(
"overall_minchi");
113 double minchi23_2019ana = mins2019[0];
115 TString ana2020 (
"/nova/ana/3flavor/Ana2020/Fits/RealData/root_contours/hist_contours_2020_" +
options +
"_"+(dmsqSurf?
"dmsq_":
"")+
"stats");
116 TFile * infile2020ana =
new TFile (ana2020 +
".root",
"read");
117 auto mins2020 =* (
TVectorD*)infile2020ana->Get(
"overall_minchi");
118 double minchi23_2020ana = mins2020[0];
121 std::vector <TString> surfNames;
122 if(!th13Surf && !dmsqSurf)surfNames = {
"delta_th23"};
123 if (dmsqSurf) surfNames.push_back(
"th23_dm32");
124 if (th13Surf) surfNames.push_back(
"th13_delta");
125 for(TString surfName:surfNames){
126 for (
int hie:{-1,+1}){
127 if((onlyNH && hie<0) || (onlyIH && hie>0))
continue;
133 auto g1_2018ana = surf_2018ana.GetGraphs(surf_1Sigma_2018ana, minchi23_2018ana);
134 auto g2_2018ana = surf_2018ana.GetGraphs(surf_2Sigma_2018ana, minchi23_2018ana);
135 auto g3_2018ana = surf_2018ana.GetGraphs(surf_3Sigma_2018ana, minchi23_2018ana);
141 auto g1_2019ana = surf_2019ana.GetGraphs(surf_1Sigma_2019ana, minchi23_2019ana);
142 auto g2_2019ana = surf_2019ana.GetGraphs(surf_2Sigma_2019ana, minchi23_2019ana);
143 auto g3_2019ana = surf_2019ana.GetGraphs(surf_3Sigma_2019ana, minchi23_2019ana);
149 auto g1_2020ana = surf_2020ana.GetGraphs(surf_1Sigma_2020ana, minchi23_2020ana);
150 auto g2_2020ana = surf_2020ana.GetGraphs(surf_2Sigma_2020ana, minchi23_2020ana);
151 auto g3_2020ana = surf_2020ana.GetGraphs(surf_3Sigma_2020ana, minchi23_2020ana);
153 TColor *c_coral2 =
new TColor(102, 255/255.0, 169/255.0, 155/255.0);
154 TColor *c_coral =
new TColor(100, 255/255.0, 228/255.0, 218/255.0);
155 TColor *c_wave =
new TColor(101, 0/255.0, 167/255.0, 200/255.0);
156 Color_t kFillTestColor = ((hie>0)? 101:102);
157 Int_t kFillColorCont = TColor::GetColorTransparent(kFillTestColor, 0.3);
158 if (onlyIH) kFillColorCont = 100;
161 Int_t kFillColor1 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.75);
162 Int_t kFillColor2 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.38);
163 Int_t kFillColor3 = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.18);
164 Int_t kDarkColor = kGray+3;
167 std::vector <TH2 *> hhist;
168 std::vector <double> chis;
170 for(
int j=0;
j< universes;
j++){
173 TFile *
infile =
new TFile (outfilename,
"read");
174 for (
int k=0; k<10; k++){
177 double minchi23 = mins[0];
180 hhist.push_back(
surf.ToTH2(minchi23));
181 chis.push_back(minchi23);
190 if (surfName==
"delta_th23"){
191 median =
new TH2D(
"me",
"", 30, 0-1/30.0, 2+1/30.0, 30, 0.3-0.2/30, 0.7+0.2/30);
192 plsigma =
new TH2D(
"p",
"", 30, 0-1/30.0, 2+1/30.0, 30, 0.3-0.2/30, 0.7+0.2/30);
193 minsigma =
new TH2D(
"mi",
"", 30, 0-1/30.0, 2+1/30.0, 30, 0.3-0.2/30, 0.7+0.2/30);}
194 if (surfName==
"th23_dm32" && onlyNH){
195 median =
new TH2D(
"me",
"", 20, 0.35-0.35/20, 0.7+0.35/20, 20, 2.1-0.6/20.0, 2.7+0.6/20.0);
196 plsigma =
new TH2D(
"p",
"", 20, 0.35-0.35/20, 0.7+0.35/20, 20, 2.1-0.6/20.0, 2.7+0.6/20.0);
197 minsigma =
new TH2D(
"mi",
"", 20, 0.35-0.35/20, 0.7+0.35/20, 20, 2.1-0.6/20.0, 2.7+0.6/20.0);}
198 if (surfName==
"th23_dm32" && onlyIH){
199 median =
new TH2D(
"me",
"", 20, 0.35-0.35/20, 0.7+0.35/20, 20, -2.7-0.5/20.0, -2.2+0.5/20.0);
200 plsigma =
new TH2D(
"p",
"", 20, 0.35-0.35/20, 0.7+0.35/20, 20, -2.7-0.5/20.0, -2.2+0.5/20.0);
201 minsigma =
new TH2D(
"mi",
"", 20, 0.35-0.35/20, 0.7+0.35/20, 20, -2.7-0.5/20.0, -2.2+0.5/20.0);}
203 for (
int i = 0;
i<=31;
i++){
204 for (
int k = 0; k<=31; k++){
206 TH1F*
m =
new TH1F(
"tmp",
"tmp", 500, 0, 200);
207 for(
int j=0;
j< hhist.size();
j++){
208 auto tmp = hhist[
j]->GetBinContent(
i, k);
211 auto med =
sigma(m, 0.5);
212 median->SetBinContent(
i, k, med);
214 plsigma->SetBinContent(
i, k,
ps);
216 minsigma->SetBinContent(
i, k,
ms);
221 cout<<
"test "<<median->GetMinimum()<<
" "<<plsigma->GetMinimum()<<
" "<<minsigma->GetMinimum()<<
endl;
223 median->Draw(
"colz");
224 cout<<
"minimum "<< median->GetMinimum()<<
endl;
225 gPad->Print((
"median_"+suffix +
"_" +
folder +
".pdf"));
228 plsigma->Draw(
"colz");
229 cout<<
"minimum "<< plsigma->GetMinimum()<<
endl;
230 gPad->Print((
"plsigma_"+suffix +
"_" +
folder +
".pdf"));
233 minsigma->Draw(
"colz");
234 cout<<
"minimum "<< minsigma->GetMinimum()<<
endl;
235 gPad->Print((
"minsigma_"+suffix +
"_" +
folder +
".pdf"));
237 for(
auto &
sigma:{1,2,3}){
239 if(surfName.Contains(
"delta_th23")) {
240 if (zoomIn)
DrawContourCanvas(surfName, 0., 2., (RHCOnly?0.15:0.275), (RHCOnly?0.8:0.725));
243 if(surfName.Contains(
"th23_dm32")) {
244 if (hie>0)
DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?1.8:2.05), (RHCOnly?3.3:2.85));
245 else DrawContourCanvas(surfName, ((nueOnly||RHCOnly)?0.2:0.3), ((nueOnly||RHCOnly)?0.8:0.7), (RHCOnly?-3.5:-2.95), (RHCOnly?-2.0:-2.25));
247 if(surfName.Contains(
"th13_delta")) {
252 kFillColor = TColor::GetColorTransparent(k2SigmaConfidenceColor, 0.75);
255 if (
sigma == 1) value = 2.3;
256 if (
sigma == 2) value = 6.18;
257 if (
sigma == 3) value = 11.83;
259 TVirtualPad*
save = gPad;
261 median->SetContour(1, &value);
262 median->Draw(
"CONT Z LIST");
265 tmp.Print((
"median_cont_" +
folder +
".pdf").c_str());
266 std::vector<TGraph*> gr_m;
267 TCollection*
specs = gROOT->GetListOfSpecials();
268 TIter nextSpec(specs);
269 while(TObject* spec = nextSpec()){
270 if(!spec->InheritsFrom(TObjArray::Class()))
continue;
271 TObjArray* conts = (TObjArray*)spec;
272 if(conts->IsEmpty())
continue;
273 if(!conts->At(0)->InheritsFrom(TList::Class()))
continue;
274 TList* cont = (TList*)conts->At(0);
276 while(TObject* obj = nextObj()){
277 if(!obj->InheritsFrom(TGraph::Class()))
continue;
278 gr_m.push_back((TGraph*)obj->Clone(
UniqueName().c_str()));
283 plsigma->SetContour(1, &value);
284 plsigma->Draw(
"CONT Z LIST");
287 tmp2.Print((
"plsigma_cont_" +
folder +
".pdf").c_str());
288 std::vector<TGraph*> gr_p;
289 TCollection* specs_p = gROOT->GetListOfSpecials();
290 TIter nextSpec_p(specs_p);
291 while(TObject* spec_p = nextSpec_p()){
292 if(!spec_p->InheritsFrom(TObjArray::Class()))
continue;
293 TObjArray* conts_p = (TObjArray*)spec_p;
294 if(conts_p->IsEmpty())
continue;
295 if(!conts_p->At(0)->InheritsFrom(TList::Class()))
continue;
296 TList* cont_p = (TList*)conts_p->At(0);
297 TIter nextObj_p(cont_p);
298 while(TObject* obj_p = nextObj_p()){
299 if(!obj_p->InheritsFrom(TGraph::Class()))
continue;
300 gr_p.push_back((TGraph*)obj_p->Clone(
UniqueName().c_str()));
305 minsigma->SetContour(1, &value);
306 minsigma->Draw(
"CONT Z LIST");
309 tmp3.Print((
"minsigma_cont_" +
folder +
".pdf").c_str());
310 std::vector<TGraph*> gr_min;
311 TCollection* specs_min = gROOT->GetListOfSpecials();
312 TIter nextSpec_min(specs_min);
313 while(TObject* spec_min = nextSpec_min()){
314 if(!spec_min->InheritsFrom(TObjArray::Class()))
continue;
315 TObjArray* conts = (TObjArray*)spec_min;
316 if(conts->IsEmpty())
continue;
317 if(!conts->At(0)->InheritsFrom(TList::Class()))
continue;
318 TList* cont = (TList*)conts->At(0);
320 while(TObject* obj = nextObj()){
321 if(!obj->InheritsFrom(TGraph::Class()))
continue;
322 gr_min.push_back((TGraph*)obj->Clone(
UniqueName().c_str()));
326 cout<<
"size of gr "<<gr_m.size()<<
" p "<<gr_p.size()<<
" min "<<gr_min.size()<<
endl;
329 if(surfName.Contains(
"delta_th23") && onlyNH &&
sigma > 2) {
330 JoinGraphs(gr_min[0], gr_min[1], kFillColorCont)->Draw(
"f");
331 for(TGraph*
g: gr_min){
333 g->SetLineStyle(kSolid);
334 g->SetLineColor(kFillColorCont);
339 for(TGraph*
g: gr_min){
341 g->SetLineStyle(kSolid);
342 g->SetLineColor(kFillColorCont);
343 g->SetFillColor(kFillColorCont);
348 for(TGraph*
g: gr_p){
350 g->SetLineStyle(kSolid);
351 g->SetLineColor(kFillColorCont);
352 g->SetFillColor(kWhite);
359 if(surfName.Contains(
"delta_th23") &&
sigma >3 && onlyIH || surfName.Contains(
"delta_th23") && onlyNH) {
360 JoinGraphs(gr_min[0], gr_min[1], kFillColorCont)->Draw(
"f");
361 for(TGraph*
g: gr_min){
363 g->SetLineStyle(kSolid);
364 g->SetLineColor(kFillColorCont);
369 for(TGraph*
g: gr_min){
371 g->SetLineStyle(kSolid);
372 g->SetLineColor(kFillColorCont);
373 g->SetFillColor(kFillColorCont);
378 if (surfName.Contains(
"delta_th23") &&
sigma >2 && onlyNH){
380 for(TGraph*
g: gr_p){
382 g->SetLineStyle(kSolid);
383 g->SetLineColor(kFillColorCont);
388 for(TGraph*
g: gr_p){
390 g->SetLineStyle(kSolid);
391 g->SetLineColor(kFillColorCont);
392 g->SetFillColor(kWhite);
399 for(TGraph*
g: gr_m){
401 g->SetLineStyle(kSolid);
402 g->SetLineColor(kFillColor);
406 auto gr2018 = g1_2018ana;
407 auto gr2019 = g1_2019ana;
408 auto gr2020 = g1_2020ana;
409 if (
sigma==1)
cout<<
"2018 size "<<gr2018.size()<<
"2019 size "<<gr2019.size()<<
" 2020 size "<<gr2020.size()<<
endl;
410 if (
sigma==2) {gr2018 = g2_2018ana; gr2019 = g2_2019ana; gr2020 = g2_2020ana;
cout<<
"2018 size "<<gr2018.size()<<
" 2019 size "<<gr2019.size() <<
" 2020 size "<<gr2020.size() <<
endl;}
411 if (
sigma==3) {gr2018 = g3_2018ana; gr2019 = g3_2019ana; gr2020 = g3_2020ana;
cout<<
"2018 size "<<gr2018.size()<<
" 2019 size "<<gr2019.size() <<
" 2020 size "<<gr2020.size() <<
endl; }
413 for (
auto &
g:gr2018) {
g->SetLineColor(kGray);
g->Draw(
"l");}
414 for (
auto &
g:gr2019) {
g->SetLineColor(kGray+2);
g->Draw(
"l");}
415 for (
auto &
g:gr2020) {
g->SetLineColor(
kViolet-6);
g->Draw(
"l");}
417 TLegend *
leg =
new TLegend(0.16,0.17,0.8,0.28);
419 leg->SetFillStyle(0);
420 leg->SetMargin(1.3*leg->GetMargin());
422 leg->SetMargin(1.4*leg->GetMargin());
423 TGraph *
dummy =
new TGraph;
424 dummy->SetLineWidth(3);
425 dummy->SetLineColor(kFillColor);
426 leg->AddEntry(dummy->Clone(),
"median",
"l");
427 dummy->SetFillColor(kFillColorCont);
428 leg->AddEntry(dummy->Clone(),
"68\% exp",
"f");
431 TLegend * leg2 =
new TLegend(0.16,0.8,0.9,0.85);
433 leg2->SetNColumns(3);
434 dummy->SetLineColor(kGray);
435 leg2->AddEntry(dummy->Clone(),
"2018 stat.fit",
"l");
436 dummy->SetLineColor(kGray+2);
437 leg2->AddEntry(dummy->Clone(),
"2019 stat fit",
"l");
438 dummy->SetLineColor(
kViolet-6);
439 leg2->AddEntry(dummy->Clone(),
"2020 stat fit",
"l");
447 gPad->Print(outdir +
"contour_" + suffix +
"_" +
folder +
"_" +
450 (zoomIn?
"_zoom":
"") +
460 TLatex * ltxNH =
new TLatex(0.16,0.89,
"#bf{NH}");
463 ltxNH->SetTextAlign(22);
465 if(th13) ltxNH->Draw();
466 else if(!vert) ltxNH->DrawLatex(.85,.20,
"#bf{NH}");
467 else ltxNH->DrawLatex(.85,.85,
"#bf{NH}");
470 TLatex * ltxIH =
new TLatex (0.16,0.46,
"#bf{IH}");
473 ltxIH->SetTextAlign(22);
475 if(th13) ltxIH->Draw();
476 else if(!vert) ltxIH->DrawLatex(.85,.20,
"#bf{IH}");
477 else ltxIH->DrawLatex(.85,.85,
"#bf{IH}");
const Color_t kPrimColorIH
const Color_t kPrimColorNH
static constexpr Double_t ms
Cuts and Vars for the 2020 FD DiF Study.
TH2 * Gaussian68Percent2D(const FrequentistSurface &s)
Up-value surface for 68% confidence in 2D in gaussian approximation.
const Color_t k2SigmaConfidenceColorNH
TH2 * Gaussian2Sigma2D(const FrequentistSurface &s)
Up-value surface for 2 sigma confidence in 2D in gaussian approximation.
string outfilename
knobs that need extra care
const XML_Char int const XML_Char * value
TH2 * Gaussian3Sigma2D(const FrequentistSurface &s)
Up-value surface for 3 sigma confidence in 2D in gaussian approximation.
void save(bool IsData, bool Bad)
void DrawContourCanvas(TString surfName, double minx=0, double maxx=2, double miny=0, double maxy=1)
TGraph * JoinGraphs(TGraph *a, TGraph *b, int fillcol)
Join graphs and set a fill color. Useful for contours.
assert(nhit_max >=nhit_nbins)
const Float_t kBlessedLabelFontSize
std::string to_string(ModuleType mt)
static std::unique_ptr< FrequentistSurface > LoadFrom(TDirectory *dir, const std::string &name)
std::string UniqueName()
Return a different string each time, for creating histograms.
const Color_t k2SigmaConfidenceColorIH