8#pragma GCC diagnostic ignored "-Wfloat-conversion"
45 for (
int i = 0; i <
nDraw; ++i)
50 const std::string CurrentName =
ParamNames[ParamEnum][ParamNo].Data();
53 if (CurrentName ==
"sin2th_13") {
56 }
else if (CurrentName ==
"sin2th_12") {
59 }
else if (CurrentName ==
"sin2th_23") {
62 }
else if (CurrentName ==
"delta_cp") {
65 }
else if (CurrentName ==
"delm2_23") {
76 Chain->SetAlias(
"J_cp",
"TMath::Sqrt(sin2th_13)*TMath::Sqrt(1.-sin2th_13)*TMath::Sqrt(1.-sin2th_13)*TMath::Sqrt(sin2th_12)*TMath::Sqrt(1.-sin2th_12)*TMath::Sqrt(sin2th_23)*TMath::Sqrt(1.-sin2th_23)*TMath::Sin(delta_cp)");
99 const double s13 = std::sqrt(s2th13);
100 const double s23 = std::sqrt(s2th23);
101 const double s12 = std::sqrt(s2th12);
102 const double sdcp = std::sin(dcp);
103 const double c13 = std::sqrt(1.-s2th13);
104 const double c12 = std::sqrt(1.-s2th12);
105 const double c23 = std::sqrt(1.-s2th23);
107 const double j = s13*c13*c13*s12*c12*s23*c23*sdcp;
116 double Prior = 1.0, PriorError = 1.0;
117 bool FlatPrior =
false;
124 FlatPrior =
ParamFlat[ParType][ParamTemp];
127 return randGen->Uniform(FlatBounds[0], FlatBounds[1]);
130 return randGen->Gaus(Prior, PriorError);
145 MACH3LOG_WARN(
"Will not {}, as oscillation parameters are missing", __func__);
150 TDirectory *JarlskogDir =
OutputFile->mkdir(
"Jarlskog");
153 bool DoReweight =
false;
158 Chain->SetBranchStatus(
"*",
false);
175 Chain->SetBranchStatus(
"step",
true);
176 Chain->SetBranchAddress(
"step", &step);
178 if (
Chain->GetBranch(
"Weight")) {
179 MACH3LOG_CRITICAL(
"Found Weight in chain, using terrible hardcoding for RC prior...");
180 Chain->SetBranchStatus(
"Weight",
true);
181 Chain->SetBranchAddress(
"Weight", &weight);
188 auto jarl = std::make_unique<TH1D>(
"jarl",
"jarl", 1000, -0.05, 0.05);
189 jarl->SetDirectory(
nullptr);
190 auto jarl_th23 = std::make_unique<TH2D>(
"jarl_th23",
"jarl_th23", 500, -0.05, 0.05, 500, 0.3, 0.7);
191 jarl_th23->SetDirectory(
nullptr);
192 auto jarl_dcp = std::make_unique<TH2D>(
"jarl_dcp",
"jarl_dcp", 500, -0.05, 0.05, 500, -1. * TMath::Pi(), TMath::Pi());
193 jarl_dcp->SetDirectory(
nullptr);
195 jarl->SetTitle(
"Jarlskog Invariant;J #equiv s_{13}c_{13}^{2}s_{12}c_{12}s_{23}c_{23}sin#delta;Posterior probability");
196 jarl_th23->SetTitle(
"Jarlskog Invariant;J #equiv s_{13}c_{13}^{2}s_{12}c_{12}s_{23}c_{23}sin#delta;Posterior probability");
199 auto jarl_IH =
M3::Clone(jarl.get(),
"jarl_IH");
200 auto jarl_NH =
M3::Clone(jarl.get(),
"jarl_NH");
202 auto jarl_th23_IH =
M3::Clone(jarl_th23.get(),
"jarl_th23_IH");
203 auto jarl_th23_NH =
M3::Clone(jarl_th23.get(),
"jarl_th23_NH");
205 auto jarl_dcp_IH =
M3::Clone(jarl_dcp.get(),
"jarl_dcp_IH");
206 auto jarl_dcp_NH =
M3::Clone(jarl_dcp.get(),
"jarl_dcp_NH");
208 auto jarl_flatsindcp =
M3::Clone(jarl.get(),
"jarl_flatsindcp");
209 auto jarl_IH_flatsindcp =
M3::Clone(jarl.get(),
"jarl_IH_flatsindcp");
210 auto jarl_NH_flatsindcp =
M3::Clone(jarl.get(),
"jarl_NH_flatsindcp");
212 auto jarl_th23_flatsindcp =
M3::Clone(jarl_th23.get(),
"jarl_th23_flatsindcp");
213 auto jarl_th23_IH_flatsindcp =
M3::Clone(jarl_th23.get(),
"jarl_th23_IH_flatsindcp");
214 auto jarl_th23_NH_flatsindcp =
M3::Clone(jarl_th23.get(),
"jarl_th23_NH_flatsindcp");
216 auto jarl_prior =
M3::Clone(jarl.get(),
"jarl_prior");
218 std::unique_ptr<TH1D> jarl_wRC_prior, jarl_wRC_prior_flatsindcp, jarl_wRC_prior_t2kth23;
221 jarl_wRC_prior =
M3::Clone(jarl.get(),
"jarl_wRC_prior");
222 jarl_wRC_prior_flatsindcp =
M3::Clone(jarl.get(),
"jarl_wRC_prior_flatsindcp");
223 jarl_wRC_prior_t2kth23 =
M3::Clone(jarl.get(),
"jarl_wRC_prior_flatsindcp");
227 auto prior3 = std::make_unique<TF1>(
"prior3",
"TMath::Abs(TMath::Cos(x))");
230 auto randGen = std::make_unique<TRandom3>(0);
231 const Long64_t countwidth =
nEntries/5;
234 if (i % countwidth == 0) {
243 const double j =
CalcJarlskog(s2th13, s2th23, s2th12, dcp);
244 const double prior_weight = prior3->Eval(dcp);
246 jarl->Fill(j, weight);
247 jarl_th23->Fill(j, s2th23,weight);
248 jarl_dcp->Fill(j, dcp,weight);
250 jarl_flatsindcp->Fill(j, prior_weight*weight);
251 jarl_th23_flatsindcp->Fill(j, s2th23,prior_weight*weight);
258 const double prior_sindcp = randGen->Uniform(-1.,1.);
260 const double prior_s13 = std::sqrt(prior_s2th13);
262 const double prior_s23 = std::sqrt(prior_s2th23);
263 const double prior_s12 = std::sqrt(prior_s2th12);
264 const double prior_sdcp = std::sin(prior_dcp);
265 const double prior_c13 = std::sqrt(1.-prior_s2th13);
266 const double prior_c12 = std::sqrt(1.-prior_s2th12);
267 const double prior_c23 = std::sqrt(1.-prior_s2th23);
268 const double prior_j = prior_s13*prior_c13*prior_c13*prior_s12*prior_c12*prior_s23*prior_c23*prior_sdcp;
270 jarl_prior->Fill(prior_j);
274 const double prior_wRC_s2th13 = randGen->Gaus(0.0220,0.0007);
275 const double prior_wRC_s13 = std::sqrt(prior_wRC_s2th13);
276 const double prior_wRC_c13 = std::sqrt(1.-prior_wRC_s2th13);
277 const double prior_wRC_j = prior_wRC_s13*prior_wRC_c13*prior_wRC_c13*prior_s12*prior_c12*prior_s23*prior_c23*prior_sdcp;
278 const double prior_wRC_flatsindcp_j = prior_wRC_s13*prior_wRC_c13*prior_wRC_c13*prior_s12*prior_c12*prior_s23*prior_c23*prior_sindcp;
280 const double s23 = std::sqrt(s2th23);
281 const double c23 = std::sqrt(1.-s2th23);
283 jarl_wRC_prior->Fill(prior_wRC_j);
284 jarl_wRC_prior_flatsindcp->Fill(prior_wRC_flatsindcp_j);
285 jarl_wRC_prior_t2kth23->Fill(prior_wRC_s13*prior_wRC_c13*prior_wRC_c13*prior_s12*prior_c12*s23*c23*prior_sdcp);
290 jarl_NH->Fill(j,weight);
291 jarl_th23_NH->Fill(j,s2th23,weight);
292 jarl_dcp_NH->Fill(j,dcp,weight);
293 jarl_NH_flatsindcp->Fill(j,prior_weight*weight);
294 jarl_th23_NH_flatsindcp->Fill(j,s2th23,prior_weight*weight);
297 jarl_IH->Fill(j,weight);
298 jarl_th23_IH->Fill(j,s2th23,weight);
299 jarl_dcp_IH->Fill(j,dcp,weight);
300 jarl_IH_flatsindcp->Fill(j, prior_weight*weight);
301 jarl_th23_IH_flatsindcp->Fill(j, s2th23, prior_weight*weight);
305 jarl->Write(
"jarlskog_both");
306 jarl_NH->Write(
"jarlskog_NH");
307 jarl_IH->Write(
"jarlskog_IH");
308 jarl_th23->Write(
"jarlskog_th23_both");
309 jarl_th23_NH->Write(
"jarlskog_th23_NH");
310 jarl_th23_IH->Write(
"jarlskog_th23_IH");
312 jarl_dcp->Write(
"jarlskog_dcp_both");
313 jarl_dcp_NH->Write(
"jarlskog_dcp_NH");
314 jarl_dcp_IH->Write(
"jarlskog_dcp_IH");
317 jarl_flatsindcp->Write(
"jarlskog_both_flatsindcp");
318 jarl_NH_flatsindcp->Write(
"jarlskog_NH_flatsindcp");
319 jarl_IH_flatsindcp->Write(
"jarlskog_IH_flatsindcp");
320 jarl_th23_flatsindcp->Write(
"jarlskog_th23_both_flatsindcp");
321 jarl_th23_NH_flatsindcp->Write(
"jarlskog_th23_NH_flatsindcp");
322 jarl_th23_IH_flatsindcp->Write(
"jarlskog_th23_IH_flatsindcp");
324 jarl_prior->Write(
"jarl_prior");
327 jarl_wRC_prior->Write(
"jarl_wRC_prior");
328 jarl_wRC_prior_flatsindcp->Write(
"jarl_wRC_prior_flatsindcp");
329 jarl_wRC_prior_t2kth23->Write(
"jarl_wRC_prior_t2kth23");
333 jarl_NH, jarl_NH_flatsindcp,
334 jarl_IH, jarl_IH_flatsindcp);
336 JarlskogDir->Close();
339 Chain->SetBranchStatus(
"*",
true);
346 const std::unique_ptr<TH1D>& jarl_flatsindcp,
347 const std::unique_ptr<TH1D>& jarl_NH,
348 const std::unique_ptr<TH1D>& jarl_NH_flatsindcp,
349 const std::unique_ptr<TH1D>& jarl_IH,
350 const std::unique_ptr<TH1D>& jarl_IH_flatsindcp) {
353 int originalErrorLevel = gErrorIgnoreLevel;
354 gErrorIgnoreLevel = kFatal;
357 for(
int hierarchy = -1; hierarchy <= 1; hierarchy++)
359 std::unique_ptr<TH1D> j_hist;
360 std::unique_ptr<TH1D> j_hist_sdcp;
363 j_hist_sdcp =
M3::Clone(jarl_NH_flatsindcp.get(),
"");
364 j_hist->SetTitle(
";J_{CP} #equiv s_{13}c^{2}_{13}s_{12}c_{12}s_{23}c_{23}sin#delta_{CP};Posterior probability");
365 }
else if(hierarchy == 0) {
367 j_hist_sdcp =
M3::Clone(jarl_flatsindcp.get(),
"");
368 j_hist->SetTitle(
";J_{CP} #equiv s_{13}c^{2}_{13}s_{12}c_{12}s_{23}c_{23}sin#delta_{CP};Posterior probability");
369 }
else if(hierarchy == -1) {
371 j_hist_sdcp =
M3::Clone(jarl_IH_flatsindcp.get(),
"");
372 j_hist->SetTitle(
";J_{CP} #equiv s_{13}c^{2}_{13}s_{12}c_{12}s_{23}c_{23}sin#delta_{CP};Posterior probability");
374 MACH3LOG_ERROR(
"Invalid hierarchy option. 1 for NH, 0 for both, -1 for IH");
379 j_hist_sdcp->Rebin(7);
381 j_hist->SetLineColor(kAzure-2);
382 j_hist_sdcp->SetLineColor(kOrange+1);
383 j_hist->SetLineWidth(2);
384 j_hist_sdcp->SetLineWidth(2);
386 auto StyleAxis = [](TH1* h) {
387 auto xAxis = h->GetXaxis();
388 auto yAxis = h->GetYaxis();
390 xAxis->SetLabelSize(0.04);
391 xAxis->SetLabelFont(132);
392 xAxis->SetTitleSize(0.04);
393 xAxis->SetTitleOffset(0.80);
394 xAxis->SetTitleFont(132);
395 xAxis->SetNdivisions(505);
396 xAxis->SetTickSize(0.04);
398 yAxis->SetLabelSize(0.04);
399 yAxis->SetLabelFont(132);
400 yAxis->SetTitleSize(0.04);
401 yAxis->SetTitleOffset(1.2);
402 yAxis->SetTitleFont(132);
403 yAxis->SetNdivisions(505);
404 yAxis->SetTickSize(0.04);
407 StyleAxis(j_hist.get());
409 j_hist->GetXaxis()->SetRangeUser(-0.04,0.04);
410 j_hist->Scale(1./j_hist->Integral());
411 j_hist_sdcp->Scale(1./j_hist_sdcp->Integral());
413 std::unique_ptr<TH1D> j_hist_copy =
M3::Clone(j_hist.get(),
"j_hist_copy");
414 std::unique_ptr<TH1D> j_hist_1sig =
M3::Clone(j_hist.get(),
"j_hist_1sig");
415 std::unique_ptr<TH1D> j_hist_2sig =
M3::Clone(j_hist.get(),
"j_hist_2sig");
416 std::unique_ptr<TH1D> j_hist_3sig =
M3::Clone(j_hist.get(),
"j_hist_3sig");
419 double j_bf = j_hist_copy->GetXaxis()->GetBinCenter(j_hist_copy->GetMaximumBin());
420 double j_1sig_low = 9999999.;
421 double j_1sig_up = -9999999.;
422 double j_2sig_low = 9999999.;;
423 double j_2sig_up = -9999999.;
424 double j_3sig_low = 9999999.;;
425 double j_3sig_up = -9999999.;
428 std::unique_ptr<TH1D> j_hist_sdcp_copy =
M3::Clone(j_hist_sdcp.get(),
"j_hist_sdcp_copy");
429 std::unique_ptr<TH1D> j_hist_sdcp_1sig =
M3::Clone(j_hist_sdcp.get(),
"j_hist_sdcp_1sig");
430 std::unique_ptr<TH1D> j_hist_sdcp_2sig =
M3::Clone(j_hist_sdcp.get(),
"j_hist_sdcp_2sig");
431 std::unique_ptr<TH1D> j_hist_sdcp_3sig =
M3::Clone(j_hist_sdcp.get(),
"j_hist_sdcp_3sig");
434 double j_sdcp_1sig_low = 9999999.;
435 double j_sdcp_1sig_up = -9999999.;
436 double j_sdcp_2sig_low = 9999999.;;
437 double j_sdcp_2sig_up = -9999999.;
438 double j_sdcp_3sig_low = 9999999.;;
439 double j_sdcp_3sig_up = -9999999.;
441 double contlevel1 = 0.68;
442 double contlevel2 = 0.90;
443 double contlevel4 = 0.99;
444 double contlevel5 = 0.9973;
445 double integral, tsum = 0.;
447 integral = j_hist_copy->Integral();
449 while((tsum/integral)<contlevel5) {
450 double tmax = j_hist_copy->GetMaximum();
451 int bin = j_hist_copy->GetMaximumBin();
452 double xval = j_hist_copy->GetXaxis()->GetBinCenter(bin);
453 double xwidth = j_hist_copy->GetXaxis()->GetBinWidth(bin);
454 if((tsum/integral)<contlevel1) {
455 j_hist_copy->SetBinContent(bin,-1.0);
456 j_hist_1sig->SetBinContent(bin,0.);
457 j_hist_2sig->SetBinContent(bin,0.);
458 j_hist_3sig->SetBinContent(bin,0.);
459 if(xval<j_1sig_low && xval<j_bf) j_1sig_low = xval - xwidth/2.;
460 if(xval>j_1sig_up && xval>j_bf) j_1sig_up = xval + xwidth/2.;
462 if((tsum/integral)<contlevel2 && (tsum / integral > contlevel1) ) {
463 j_hist_copy->SetBinContent(bin,-5.0);
464 j_hist_2sig->SetBinContent(bin,0.);
465 j_hist_3sig->SetBinContent(bin,0.);
466 if(xval<j_2sig_low && xval<j_bf) j_2sig_low = xval - xwidth/2.;
467 if(xval>j_2sig_up && xval>j_bf) j_2sig_up = xval + xwidth/2.;
469 if((tsum/integral)<contlevel4 && (tsum / integral > contlevel1) ) {
470 j_hist_copy->SetBinContent(bin,-9.0);
471 j_hist_3sig->SetBinContent(bin,0.);
472 if(xval < j_3sig_low && xval <j_bf) j_3sig_low = xval - xwidth/2.;
473 if(xval > j_3sig_up && xval > j_bf) j_3sig_up = xval + xwidth/2.;
478 integral = j_hist_sdcp_copy->Integral();
481 while((tsum/integral)<contlevel5) {
482 double tmax = j_hist_sdcp_copy->GetMaximum();
483 int bin = j_hist_sdcp_copy->GetMaximumBin();
484 double xval = j_hist_sdcp_copy->GetXaxis()->GetBinCenter(bin);
485 double xwidth = j_hist_sdcp_copy->GetXaxis()->GetBinWidth(bin);
486 if((tsum/integral)<contlevel1) {
487 j_hist_sdcp_copy->SetBinContent(bin,-1.0);
488 j_hist_sdcp_1sig->SetBinContent(bin,0.);
489 j_hist_sdcp_2sig->SetBinContent(bin,0.);
490 j_hist_sdcp_3sig->SetBinContent(bin,0.);
491 if(xval<j_sdcp_1sig_low && xval<j_bf) j_sdcp_1sig_low = xval - xwidth/2.;
492 if(xval>j_sdcp_1sig_up && xval>j_bf) j_sdcp_1sig_up = xval + xwidth/2.;
494 if((tsum/integral)<contlevel2 && (tsum / integral > contlevel1) ) {
495 j_hist_sdcp_copy->SetBinContent(bin,-5.0);
496 j_hist_sdcp_2sig->SetBinContent(bin,0.);
497 j_hist_sdcp_3sig->SetBinContent(bin,0.);
498 if(xval<j_sdcp_2sig_low && xval<j_bf) j_sdcp_2sig_low = xval - xwidth/2.;
499 if(xval>j_sdcp_2sig_up && xval>j_bf) j_sdcp_2sig_up = xval + xwidth/2.;
501 if((tsum/integral)<contlevel4 && (tsum / integral > contlevel1) ) {
502 j_hist_sdcp_copy->SetBinContent(bin,-9.0);
503 j_hist_sdcp_3sig->SetBinContent(bin,0.);
504 if(xval<j_sdcp_3sig_low && xval<j_bf) j_sdcp_3sig_low = xval - xwidth/2.;
505 if(xval>j_sdcp_3sig_up && xval>j_bf) j_sdcp_3sig_up = xval + xwidth/2.;
510 j_hist_1sig->SetLineStyle(9);
511 j_hist_sdcp_1sig->SetLineStyle(9);
512 j_hist_2sig->SetLineStyle(7);
513 j_hist_sdcp_2sig->SetLineStyle(7);
514 j_hist_3sig->SetLineStyle(2);
515 j_hist_sdcp_3sig->SetLineStyle(2);
517 auto ldash = std::make_unique<TH1D>(
"ldash",
"solid", 10, -0.04, 0.04);
518 auto sdash = std::make_unique<TH1D>(
"sdash",
"dashed", 10, -0.04, 0.04);
519 auto fdash = std::make_unique<TH1D>(
"fdash",
"fdashed",10, -0.04, 0.04);
520 ldash->SetLineColor(kBlack);
521 sdash->SetLineColor(kBlack);
522 fdash->SetLineColor(kBlack);
523 ldash->SetLineWidth(2);
524 sdash->SetLineWidth(2);
525 fdash->SetLineWidth(2);
526 ldash->SetLineStyle(9);
527 sdash->SetLineStyle(7);
528 fdash->SetLineStyle(2);
530 double vertUp = 0.5 * j_hist->GetMaximum();
531 auto jline_1sig_low = std::make_unique<TLine>(j_1sig_low, 0., j_1sig_low, vertUp);
532 auto jline_2sig_low = std::make_unique<TLine>(j_2sig_low, 0., j_2sig_low, vertUp);
533 auto jline_3sig_low = std::make_unique<TLine>(j_3sig_low, 0., j_3sig_low, vertUp);
535 auto jline_1sig_up = std::make_unique<TLine>(j_1sig_up, 0., j_1sig_up,vertUp);
536 auto jline_2sig_up = std::make_unique<TLine>(j_2sig_up, 0., j_2sig_up,vertUp);
537 auto jline_3sig_up = std::make_unique<TLine>(j_3sig_up, 0., j_3sig_up,vertUp);
539 auto jline_sdcp_1sig_low = std::make_unique<TLine>(j_sdcp_1sig_low, 0., j_sdcp_1sig_low, vertUp);
540 auto jline_sdcp_2sig_low = std::make_unique<TLine>(j_sdcp_2sig_low, 0., j_sdcp_2sig_low, vertUp);
541 auto jline_sdcp_3sig_low = std::make_unique<TLine>(j_sdcp_3sig_low, 0., j_sdcp_3sig_low, vertUp);
543 auto jline_sdcp_1sig_up = std::make_unique<TLine>(j_sdcp_1sig_up, 0., j_sdcp_1sig_up, vertUp);
544 auto jline_sdcp_2sig_up = std::make_unique<TLine>(j_sdcp_2sig_up, 0., j_sdcp_2sig_up, vertUp);
545 auto jline_sdcp_3sig_up = std::make_unique<TLine>(j_sdcp_3sig_up, 0., j_sdcp_3sig_up, vertUp);
547 double arrowLength = 0.003;
548 double arrowHeight = vertUp;
550 auto MakeArrow = [&](
double x, Color_t color, Width_t width) -> std::unique_ptr<TArrow> {
551 auto arrow = std::make_unique<TArrow>(x, arrowHeight, x - arrowLength, arrowHeight, 0.02,
">");
552 arrow->SetLineColor(color);
553 arrow->SetLineWidth(width);
557 auto j_arrow_1sig_up = MakeArrow(j_1sig_up, j_hist_1sig->GetLineColor(), j_hist_1sig->GetLineWidth());
558 auto j_arrow_2sig_up = MakeArrow(j_2sig_up, j_hist_2sig->GetLineColor(), j_hist_2sig->GetLineWidth());
559 auto j_arrow_3sig_up = MakeArrow(j_3sig_up, j_hist_3sig->GetLineColor(), j_hist_3sig->GetLineWidth());
561 auto j_sdcp_arrow_1sig_up = MakeArrow(j_sdcp_1sig_up, j_hist_sdcp_1sig->GetLineColor(), j_hist_sdcp_1sig->GetLineWidth());
562 auto j_sdcp_arrow_2sig_up = MakeArrow(j_sdcp_2sig_up, j_hist_sdcp_2sig->GetLineColor(), j_hist_sdcp_2sig->GetLineWidth());
563 auto j_sdcp_arrow_3sig_up = MakeArrow(j_sdcp_3sig_up, j_hist_sdcp_3sig->GetLineColor(), j_hist_sdcp_3sig->GetLineWidth());
565 MACH3LOG_DEBUG(
"j_1sig_low = {}, j_2sig_low = {}, j_3sig_low = {}", j_1sig_low, j_2sig_low, j_3sig_low);
566 MACH3LOG_DEBUG(
"j_1sig_up = {}, j_2sig_up = {}, j_3sig_up = {}", j_1sig_up, j_2sig_up, j_3sig_up);
568 auto CopyLineStyle = [](
const TH1D* src, TLine* dst) {
569 dst->SetLineColor(src->GetLineColor());
570 dst->SetLineStyle(src->GetLineStyle());
571 dst->SetLineWidth(src->GetLineWidth());
574 CopyLineStyle(j_hist_1sig.get(), jline_1sig_low.get());
575 CopyLineStyle(j_hist_1sig.get(), jline_1sig_up.get());
576 CopyLineStyle(j_hist_2sig.get(), jline_2sig_low.get());
577 CopyLineStyle(j_hist_2sig.get(), jline_2sig_up.get());
578 CopyLineStyle(j_hist_3sig.get(), jline_3sig_low.get());
579 CopyLineStyle(j_hist_3sig.get(), jline_3sig_up.get());
581 CopyLineStyle(j_hist_sdcp_1sig.get(), jline_sdcp_1sig_low.get());
582 CopyLineStyle(j_hist_sdcp_1sig.get(), jline_sdcp_1sig_up.get());
583 CopyLineStyle(j_hist_sdcp_2sig.get(), jline_sdcp_2sig_low.get());
584 CopyLineStyle(j_hist_sdcp_2sig.get(), jline_sdcp_2sig_up.get());
585 CopyLineStyle(j_hist_sdcp_3sig.get(), jline_sdcp_3sig_low.get());
586 CopyLineStyle(j_hist_sdcp_3sig.get(), jline_sdcp_3sig_up.get());
588 auto leg = std::make_unique<TLegend>(0.45, 0.60, 0.75, 0.90);
589 leg->SetTextSize(0.05);
590 leg->SetFillStyle(0);
592 leg->SetTextFont(132);
593 leg->SetBorderSize(0);
595 leg->AddEntry(j_hist.get(),
"Prior flat in #delta_{CP}",
"l");
596 leg->AddEntry(j_hist_sdcp.get(),
"Prior flat in sin#delta_{CP}",
"l");
597 leg->AddEntry(ldash.get(),
"68% CI",
"l");
598 leg->AddEntry(sdash.get(),
"90% CI",
"l");
599 leg->AddEntry(fdash.get(),
"99% CI",
"l");
601 j_hist->GetYaxis()->SetRangeUser(0., j_hist->GetMaximum()*1.15);
603 j_hist_sdcp->Draw(
"same h");
605 jline_sdcp_1sig_up->Draw(
"same");
606 jline_sdcp_2sig_up->Draw(
"same");
607 jline_sdcp_3sig_up->Draw(
"same");
608 jline_1sig_up->Draw(
"same");
609 jline_2sig_up->Draw(
"same");
610 jline_3sig_up->Draw(
"same");
612 j_arrow_1sig_up->Draw();
613 j_arrow_2sig_up->Draw();
614 j_arrow_3sig_up->Draw();
615 j_sdcp_arrow_1sig_up->Draw();
616 j_sdcp_arrow_2sig_up->Draw();
617 j_sdcp_arrow_3sig_up->Draw();
620 auto ttext = std::make_unique<TText>();
622 ttext->SetTextSize(0.03);
623 ttext->SetTextAlign(13);
625 if (hierarchy == 1) ttext->DrawText(0.15, 0.85,
"Normal Ordering");
626 else if (hierarchy == 0) ttext->DrawText(0.15, 0.85,
"Both Orderings");
627 else if (hierarchy == -1) ttext->DrawText(0.15, 0.85,
"Inverted Ordering");
635 if(hierarchy == 1)
Posterior->Write(
"jarl1D_NH_comp");
636 else if(hierarchy == 0)
Posterior->Write(
"jarl1D_both_comp");
637 else if(hierarchy == -1)
Posterior->Write(
"jarl1D_IH_comp");
640 gErrorIgnoreLevel = originalErrorLevel;
#define _MaCh3_Safe_Include_Start_
KS: Avoiding warning checking for headers.
#define _MaCh3_Safe_Include_End_
KS: Restore warning checking after including external headers.
ParameterEnum
KS: Enum for different covariance classes.
#define MACH3LOG_CRITICAL
Class responsible for processing MCMC chains, performing diagnostics, generating plots,...
void GetNthParameter(const int param, double &Prior, double &PriorError, TString &Title) const
Get properties of parameter by passing it number.
int BurnInCut
Value of burn in cut.
std::vector< std::vector< double > > ParamCentral
Parameters central values which we are going to analyse.
std::vector< std::vector< double > > ParamErrors
Uncertainty on a single parameter.
std::vector< int > nParam
Number of parameters per type.
std::vector< std::vector< bool > > ParamFlat
Whether Param has flat prior or not.
std::vector< std::vector< double > > ParamNom
std::unique_ptr< TCanvas > Posterior
Fancy canvas used for our beautiful plots.
TFile * OutputFile
The output file.
TChain * Chain
Main chain storing all steps etc.
int nDraw
Number of all parameters used in the analysis.
std::vector< ParameterEnum > ParamType
Make an enum for which class this parameter belongs to so we don't have to keep string comparing.
TString CanvasName
Name of canvas which help to save to the sample pdf.
std::vector< std::string > ParameterGroup
std::vector< TString > BranchNames
std::vector< std::vector< TString > > ParamNames
Name of parameters which we are going to analyse.
std::vector< int > ParamTypeStartPos
int nEntries
KS: For merged chains number of entries will be different from nSteps.
Custom exception class for MaCh3 errors.
double SamplePriorForParam(const int paramIndex, const std::unique_ptr< TRandom3 > &randGen, const std::vector< double > &FlatBounds) const
Draw Prior value.
int DeltaCPIndex
Index of in the parameter list.
double CalcJarlskog(const double s2th13, const double s2th23, const double s2th12, const double dcp) const
Calculate Jarlskog Invariant using oscillation parameters.
std::string Sin2Theta13Name
Name of the parameter representing .
std::string Sin2Theta12Name
Name of the parameter representing .
int Sin2Theta12Index
Index of in the parameter list.
OscProcessor(const std::string &InputFile)
Constructs an OscProcessor object with the specified input file and options.
int DeltaM2_23Index
Index of in the parameter list.
std::string DeltaCPName
Name of the parameter representing (the CP-violating phase).
bool PlotJarlskog
Will plot Jarlskog Invariant using information in the chain.
int Sin2Theta13Index
Index of in the parameter list.
std::string Sin2Theta23Name
Name of the parameter representing .
void LoadAdditionalInfo() override
Read the Osc cov file and get the input central values and errors Here we allow Jarlskog Shenanigans.
void MakeJarlskogPlot(const std::unique_ptr< TH1D > &jarl, const std::unique_ptr< TH1D > &jarl_flatsindcp, const std::unique_ptr< TH1D > &jarl_NH, const std::unique_ptr< TH1D > &jarl_NH_flatsindcp, const std::unique_ptr< TH1D > &jarl_IH, const std::unique_ptr< TH1D > &jarl_IH_flatsindcp)
Perform Jarlskog Plotting.
int Sin2Theta23Index
Index of in the parameter list.
virtual ~OscProcessor()
Destroys the OscProcessor object.
std::string DeltaM2_23Name
Name of the parameter representing (mass-squared difference).
bool OscEnabled
Will plot Jarlskog Invariant using information in the chain.
void PerformJarlskogAnalysis()
Perform Several Jarlskog Plotting.
static constexpr const double _BAD_DOUBLE_
Default value used for double initialisation.
std::unique_ptr< ObjectType > Clone(const ObjectType *obj, const std::string &name="")
KS: Creates a copy of a ROOT-like object and wraps it in a smart pointer.
static constexpr const int _BAD_INT_
Default value used for int initialisation.
void PrintProgressBar(const Long64_t Done, const Long64_t All)
KS: Simply print progress bar.
void EstimateDataTransferRate(TChain *chain, const Long64_t entry)
KS: Check what CPU you are using.