MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
OscProcessor.cpp
Go to the documentation of this file.
1#include "OscProcessor.h"
2
4// ROOT includes
5#include "TArrow.h"
7
8#pragma GCC diagnostic ignored "-Wfloat-conversion"
9
10// ****************************
11OscProcessor::OscProcessor(const std::string &InputFile) : MCMCProcessor(InputFile) {
12// ****************************
13 //KS: WARNING this only work when you project from Chain, will nor work when you try SetBranchAddress etc. Turn it on only if you know how to use it
14 PlotJarlskog = false;
15
22}
23
24// ****************************
25// The destructor
27// ****************************
28
29}
30
31// ***************
32// Read the Osc cov file and get the input central values and errors
34// ***************
35 // KS: Check if OscParams were enabled, in future we will also get
36 for(size_t i = 0; i < ParameterGroup.size(); i++) {
37 if(ParameterGroup[i] == "Osc"){
38 OscEnabled = true;
39 break;
40 }
41 }
42
43 if(OscEnabled)
44 {
45 for (int i = 0; i < nDraw; ++i)
46 {
47 //Those keep which parameter type we run currently and relative number
48 const int ParamEnum = ParamType[i];
49 const int ParamNo = i - ParamTypeStartPos[ParameterEnum(ParamEnum)];
50 const std::string CurrentName = ParamNames[ParamEnum][ParamNo].Data();
51
53 if (CurrentName == "sin2th_13") {
55 Sin2Theta13Name = CurrentName;
56 } else if (CurrentName == "sin2th_12") {
58 Sin2Theta12Name = CurrentName;
59 } else if (CurrentName == "sin2th_23") {
61 Sin2Theta23Name = CurrentName;
62 } else if (CurrentName == "delta_cp") {
63 DeltaCPIndex = i;
64 DeltaCPName = CurrentName;
65 } else if (CurrentName == "delm2_23") {
67 DeltaM2_23Name = CurrentName;
68 }
69 }
70 } else{
71 MACH3LOG_WARN("Didn't find oscillation parameters");
72 }
73
75 {
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)");
77 BranchNames.push_back("J_cp");
78 ParamType.push_back(kXSecPar);
80 nDraw++;
81
83 ParamNom[kXSecPar].push_back( 0. );
84 ParamCentral[kXSecPar].push_back( 0. );
85 ParamErrors[kXSecPar].push_back( 1. );
86 // Push back the name
87 ParamNames[kXSecPar].push_back("J_cp");
88 ParamFlat[kXSecPar].push_back( false );
89 } else if(PlotJarlskog && !OscEnabled) {
90 MACH3LOG_ERROR("Trying to enable Jarlskog without oscillations");
91 throw MaCh3Exception(__FILE__,__LINE__);
92 }
93}
94
95// ***************
96// Calculate Jarlskog Invariant using oscillation parameters
97double OscProcessor::CalcJarlskog(const double s2th13, const double s2th23, const double s2th12, const double dcp) const {
98// ***************
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);
106
107 const double j = s13*c13*c13*s12*c12*s23*c23*sdcp;
108
109 return j;
110}
111
112// ***************
113double OscProcessor::SamplePriorForParam(const int paramIndex, const std::unique_ptr<TRandom3>& randGen, const std::vector<double>& FlatBounds) const {
114// ***************
115 TString Title = "";
116 double Prior = 1.0, PriorError = 1.0;
117 bool FlatPrior = false;
118
119 // Get info for this parameter
120 GetNthParameter(paramIndex, Prior, PriorError, Title);
121
122 ParameterEnum ParType = ParamType[paramIndex];
123 int ParamTemp = paramIndex - ParamTypeStartPos[ParType];
124 FlatPrior = ParamFlat[ParType][ParamTemp];
125
126 if (FlatPrior) {
127 return randGen->Uniform(FlatBounds[0], FlatBounds[1]);
128 } else {
129 // Gaussian prior centered at Prior with width PriorError
130 return randGen->Gaus(Prior, PriorError);
131 }
132}
133
134// ***************
135// Perform Several Jarlskog Plotting
137// ***************
138 if(!OscEnabled ||
144 {
145 MACH3LOG_WARN("Will not {}, as oscillation parameters are missing", __func__);
146 return;
147 }
148 MACH3LOG_INFO("Starting {}", __func__);
149
150 TDirectory *JarlskogDir = OutputFile->mkdir("Jarlskog");
151 JarlskogDir->cd();
152
153 bool DoReweight = false;
154
155 double s2th13, s2th23, s2th12, dcp, dm2 = M3::_BAD_DOUBLE_;
156 double weight = 1.;
157 int step = M3::_BAD_INT_;
158 Chain->SetBranchStatus("*", false);
159
160 Chain->SetBranchStatus(Sin2Theta13Name.c_str(), true);
161 Chain->SetBranchAddress(Sin2Theta13Name.c_str(), &s2th13);
162
163 Chain->SetBranchStatus(Sin2Theta23Name.c_str(), true);
164 Chain->SetBranchAddress(Sin2Theta23Name.c_str(), &s2th23);
165
166 Chain->SetBranchStatus(Sin2Theta12Name.c_str(), true);
167 Chain->SetBranchAddress(Sin2Theta12Name.c_str(), &s2th12);
168
169 Chain->SetBranchStatus(DeltaCPName.c_str(), true);
170 Chain->SetBranchAddress(DeltaCPName.c_str(), &dcp);
171
172 Chain->SetBranchStatus(DeltaM2_23Name.c_str(), true);
173 Chain->SetBranchAddress(DeltaM2_23Name.c_str(), &dm2);
174
175 Chain->SetBranchStatus("step", true);
176 Chain->SetBranchAddress("step", &step);
177
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);
182 DoReweight = true;
183 } else {
184 MACH3LOG_WARN("Not applying reweighting weight");
185 weight = 1.0;
186 }
187 // Original histograms
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);
194
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");
197
198 // Clones
199 auto jarl_IH = M3::Clone(jarl.get(), "jarl_IH");
200 auto jarl_NH = M3::Clone(jarl.get(), "jarl_NH");
201
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");
204
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");
207
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");
211
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");
215
216 auto jarl_prior = M3::Clone(jarl.get(), "jarl_prior");
217
218 std::unique_ptr<TH1D> jarl_wRC_prior, jarl_wRC_prior_flatsindcp, jarl_wRC_prior_t2kth23;
219 // Only use this if chain has reweigh weight [mostly coming from Reactor Constrains]
220 if(DoReweight){
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");
224 }
225
226 // to apply a prior that is flat in sin(dcp) intead of dcp
227 auto prior3 = std::make_unique<TF1>("prior3", "TMath::Abs(TMath::Cos(x))");
228
229 // T2K prior is flat (and uncorrelated) in dcp, sin^2(th13), sin^2(th23)
230 auto randGen = std::make_unique<TRandom3>(0);
231 const Long64_t countwidth = nEntries/5;
232
233 for(int i = 0; i < nEntries; ++i) {
234 if (i % countwidth == 0) {
237 } else {
238 Chain->GetEntry(i);
239 }
240
241 if(step < BurnInCut) continue; // burn-in cut
242
243 const double j = CalcJarlskog(s2th13, s2th23, s2th12, dcp);
244 const double prior_weight = prior3->Eval(dcp);
245
246 jarl->Fill(j, weight);
247 jarl_th23->Fill(j, s2th23,weight);
248 jarl_dcp->Fill(j, dcp,weight);
249
250 jarl_flatsindcp->Fill(j, prior_weight*weight);
251 jarl_th23_flatsindcp->Fill(j, s2th23,prior_weight*weight);
252
253 const double prior_s2th13 = SamplePriorForParam(Sin2Theta13Index, randGen, {0.,1.});
254 const double prior_s2th23 = SamplePriorForParam(Sin2Theta23Index, randGen, {0.,1.});
255 const double prior_s2th12 = SamplePriorForParam(Sin2Theta12Index, randGen, {0.,1.});
256 const double prior_dcp = SamplePriorForParam(DeltaCPIndex, randGen, {-1.*TMath::Pi(),TMath::Pi()});
257 // KS: This is hardcoded but we always assume flat in delta CP so probably fine
258 const double prior_sindcp = randGen->Uniform(-1.,1.);
259
260 const double prior_s13 = std::sqrt(prior_s2th13);
261
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;
269
270 jarl_prior->Fill(prior_j);
271
272 if(DoReweight){
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;
279
280 const double s23 = std::sqrt(s2th23);
281 const double c23 = std::sqrt(1.-s2th23);
282
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);
286 }
287
288
289 if(dm2 > 0.) {
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);
295 }
296 else if(dm2 < 0.) {
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);
302 }
303 }
304
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");
311
312 jarl_dcp->Write("jarlskog_dcp_both");
313 jarl_dcp_NH->Write("jarlskog_dcp_NH");
314 jarl_dcp_IH->Write("jarlskog_dcp_IH");
315
316
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");
323
324 jarl_prior->Write("jarl_prior");
325
326 if(DoReweight) {
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");
330 }
331
332 MakeJarlskogPlot(jarl, jarl_flatsindcp,
333 jarl_NH, jarl_NH_flatsindcp,
334 jarl_IH, jarl_IH_flatsindcp);
335
336 JarlskogDir->Close();
337 delete JarlskogDir;
338
339 Chain->SetBranchStatus("*", true);
340 OutputFile->cd();
341}
342
343
344// ***************
345void OscProcessor::MakeJarlskogPlot(const std::unique_ptr<TH1D>& jarl,
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) {
351// ***************
352 MACH3LOG_INFO("Starting {}", __func__);
353 int originalErrorLevel = gErrorIgnoreLevel;
354 gErrorIgnoreLevel = kFatal;
355
356 // 1-->NH, 0-->both, -1-->IH
357 for(int hierarchy = -1; hierarchy <= 1; hierarchy++)
358 {
359 std::unique_ptr<TH1D> j_hist;
360 std::unique_ptr<TH1D> j_hist_sdcp;
361 if(hierarchy == 1) {
362 j_hist = M3::Clone(jarl_NH.get(), "");
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) {
366 j_hist = M3::Clone(jarl.get(), "");
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) {
370 j_hist = M3::Clone(jarl_IH.get(), "");
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");
373 } else {
374 MACH3LOG_ERROR("Invalid hierarchy option. 1 for NH, 0 for both, -1 for IH");
375 throw MaCh3Exception(__FILE__ , __LINE__ );
376 }
377
378 j_hist->Rebin(7);
379 j_hist_sdcp->Rebin(7);
380
381 j_hist->SetLineColor(kAzure-2);
382 j_hist_sdcp->SetLineColor(kOrange+1);
383 j_hist->SetLineWidth(2);
384 j_hist_sdcp->SetLineWidth(2);
385
386 auto StyleAxis = [](TH1* h) {
387 auto xAxis = h->GetXaxis();
388 auto yAxis = h->GetYaxis();
389
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);
397
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);
405 };
406
407 StyleAxis(j_hist.get());
408
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());
412
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");
417
418 //upper and lower edges
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.;
426
427
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");
432
433 //upper and lower edges
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.;
440
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.;
446
447 integral = j_hist_copy->Integral();
448
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.;
461 }
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.;
468 }
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.;
474 }
475 tsum+=tmax;
476 }
477
478 integral = j_hist_sdcp_copy->Integral();
479 tsum = 0.;
480
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.;
493 }
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.;
500 }
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.;
506 }
507 tsum+=tmax;
508 }
509
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);
516
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);
529
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);
534
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);
538
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);
542
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);
546
547 double arrowLength = 0.003;
548 double arrowHeight = vertUp;
549
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);
554 return arrow;
555 };
556
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());
560
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());
564
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);
567
568 auto CopyLineStyle = [](const TH1D* src, TLine* dst) {
569 dst->SetLineColor(src->GetLineColor());
570 dst->SetLineStyle(src->GetLineStyle());
571 dst->SetLineWidth(src->GetLineWidth());
572 };
573
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());
580
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());
587
588 auto leg = std::make_unique<TLegend>(0.45, 0.60, 0.75, 0.90);
589 leg->SetTextSize(0.05);
590 leg->SetFillStyle(0);
591 leg->SetNColumns(1);
592 leg->SetTextFont(132);
593 leg->SetBorderSize(0);
594
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");
600
601 j_hist->GetYaxis()->SetRangeUser(0., j_hist->GetMaximum()*1.15);
602 j_hist->Draw("h");
603 j_hist_sdcp->Draw("same h");
604
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");
611
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();
618 leg->Draw("same");
619
620 auto ttext = std::make_unique<TText>();
621 ttext->SetNDC(); // Use normalized device coordinates
622 ttext->SetTextSize(0.03); // Adjust size as needed
623 ttext->SetTextAlign(13); // Align left-top
624
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");
628
629 gPad->RedrawAxis();
630 Posterior->Update();
631 gPad->Update();
632
633 Posterior->Print(CanvasName);
634
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");
638 }
639
640 gErrorIgnoreLevel = originalErrorLevel;
641}
#define _MaCh3_Safe_Include_Start_
KS: Avoiding warning checking for headers.
Definition: Core.h:106
#define _MaCh3_Safe_Include_End_
KS: Restore warning checking after including external headers.
Definition: Core.h:117
ParameterEnum
KS: Enum for different covariance classes.
Definition: MCMCProcessor.h:52
@ kXSecPar
Definition: MCMCProcessor.h:53
#define MACH3LOG_CRITICAL
Definition: MaCh3Logger.h:26
#define MACH3LOG_DEBUG
Definition: MaCh3Logger.h:22
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:24
Class responsible for processing MCMC chains, performing diagnostics, generating plots,...
Definition: MCMCProcessor.h:65
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.
Definition: OscProcessor.h:75
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 .
Definition: OscProcessor.h:58
std::string Sin2Theta12Name
Name of the parameter representing .
Definition: OscProcessor.h:60
int Sin2Theta12Index
Index of in the parameter list.
Definition: OscProcessor.h:71
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.
Definition: OscProcessor.h:77
std::string DeltaCPName
Name of the parameter representing (the CP-violating phase).
Definition: OscProcessor.h:64
bool PlotJarlskog
Will plot Jarlskog Invariant using information in the chain.
Definition: OscProcessor.h:52
int Sin2Theta13Index
Index of in the parameter list.
Definition: OscProcessor.h:69
std::string Sin2Theta23Name
Name of the parameter representing .
Definition: OscProcessor.h:62
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.
Definition: OscProcessor.h:73
virtual ~OscProcessor()
Destroys the OscProcessor object.
std::string DeltaM2_23Name
Name of the parameter representing (mass-squared difference).
Definition: OscProcessor.h:66
bool OscEnabled
Will plot Jarlskog Invariant using information in the chain.
Definition: OscProcessor.h:55
void PerformJarlskogAnalysis()
Perform Several Jarlskog Plotting.
static constexpr const double _BAD_DOUBLE_
Default value used for double initialisation.
Definition: Core.h:43
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.
Definition: Core.h:45
void PrintProgressBar(const Long64_t Done, const Long64_t All)
KS: Simply print progress bar.
Definition: Monitor.cpp:212
void EstimateDataTransferRate(TChain *chain, const Long64_t entry)
KS: Check what CPU you are using.
Definition: Monitor.cpp:195