MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
SampleSummary.cpp
Go to the documentation of this file.
2
3//this file is choc full of usage of a root interface that only takes floats, turn this warning off for this CU for now
4#pragma GCC diagnostic ignored "-Wfloat-conversion"
5#pragma GCC diagnostic ignored "-Wuseless-cast"
6
7// *******************
8// The constructor
9SampleSummary::SampleSummary(const int n_Samples, const std::string &Filename, SampleHandlerBase* const sample, const int nSteps) {
10// *******************
11 MACH3LOG_DEBUG("Making sample summary class...");
12 #ifdef MULTITHREAD
13 MACH3LOG_DEBUG("With OpenMP and {} threads", omp_get_max_threads());
14 #endif
15
17
18 if(StandardFluctuation) MACH3LOG_INFO("Using standard method of statistical fluctuation");
19 else MACH3LOG_INFO("Using alternative method of statistical fluctuation, which is much slower");
20
21 //KS: If true it will print posterior predictive for every beta parameter it is quite useful but make root big number of plots
22 DoBetaParam = true;
23 if(DoBetaParam) MACH3LOG_INFO("I will calculate #beta parameters from Barlow-Beeston formalism");
24
25 //If true code will normalise each histogram, this way you can calculate shape only error. etc. pvalue will be completely wrong unfortunately
26 doShapeOnly = false;
27
28 nChainSteps = nSteps;
29 //KS: nChainSteps == 0 means we run PriorPredcitive
30 if(nChainSteps == 0) isPriorPredictive = true;
31 else isPriorPredictive = false;
32
33 OutputName = Filename;
34 nSamples = n_Samples;
35 SampleHandler = sample;
36
37 //Get mach3 modes from manager
39
40 nThrows = 0;
41 first_pass = true;
42 Outputfile = nullptr;
43 OutputTree = nullptr;
44 rnd = std::make_unique<TRandom3>();
45
46 DataHist.resize(nSamples);
49 NominalHist.resize(nSamples);
50 PosteriorHist.resize(nSamples);
51 W2NomHist.resize(nSamples);
52 w2Hist.resize(nSamples);
53
56
57 if(DoBetaParam) BetaHist.resize(nSamples);
58
59 maxBins.resize(nSamples);
60
61 lnLHist_Mean.resize(nSamples);
62 lnLHist_Mode.resize(nSamples);
64 MeanHist.resize(nSamples);;
66 ModeHist.resize(nSamples);
67 W2MeanHist.resize(nSamples);
68 W2ModeHist.resize(nSamples);
74
75 //KS: When a histogram is created with an axis lower limit greater or equal to its upper limit ROOT will automatically adjust histogram range
76 // https://root.cern.ch/doc/master/classTH1.html#auto-bin
77 lnLHist = std::make_unique<TH1D>("lnLHist_predpredfluc", "lnLHist_predpredfluc", 100, 1, -1);
78 lnLHist->SetDirectory(nullptr);
79 lnLHist->GetXaxis()->SetTitle("-2LLH (Pred Fluc, Pred)");
80 lnLHist->GetYaxis()->SetTitle("Counts");
81
82 lnLHist_drawdata = std::make_unique<TH1D>("lnLHist_drawdata", "lnLHist_drawdata", 100, 1, -1);
83 lnLHist_drawdata->SetDirectory(nullptr);
84 lnLHist_drawdata->GetXaxis()->SetTitle("-2LLH (Data, Draw)");
85 lnLHist_drawdata->GetYaxis()->SetTitle("Counts");
86
87 lnLHist_drawfluc = std::make_unique<TH1D>("lnLHist_drawpredfluc", "lnLHist_drawpredfluc", 100, 1, -1);
88 lnLHist_drawfluc->SetDirectory(nullptr);
89 lnLHist_drawfluc->GetXaxis()->SetTitle("-2LLH (Pred Fluc, Draw)");
90 lnLHist_drawfluc->GetYaxis()->SetTitle("Counts");
91
92 lnLHist_drawflucdraw = std::make_unique<TH1D>("lnLHist_drawflucdraw", "lnLHist_drawflucdraw", 100, 1, -1);
93 lnLHist_drawflucdraw->SetDirectory(nullptr);
94 lnLHist_drawflucdraw->GetXaxis()->SetTitle("-2LLH (Draw Fluc, Draw)");
95 lnLHist_drawflucdraw->GetYaxis()->SetTitle("Counts");
96
97 lnLDrawHist = std::make_unique<TH2D>("lnLDrawHist", "lnLDrawHist", 50, 1, -1, 50, 1, -1);
98 lnLDrawHist->SetDirectory(nullptr);
99 lnLDrawHist->GetXaxis()->SetTitle("-2LLH_{Pred Fluc, Draw}");
100 lnLDrawHist->GetYaxis()->SetTitle("-2LLH_{Data, Draw}");
101
102 lnLFlucHist = std::make_unique<TH2D>("lnLFlucHist", "lnLFlucHist", 50, 1, -1, 50, 1, -1);
103 lnLFlucHist->SetDirectory(nullptr);
104 lnLFlucHist->GetXaxis()->SetTitle("-2LLH_{Draw Fluc, Draw}");
105 lnLFlucHist->GetYaxis()->SetTitle("-2LLH_{Data, Draw}");
106
107 lnLDrawHistRate = std::make_unique<TH2D>("lnLDrawHistRate", "lnLDrawHistRate", 50, 1, -1, 50, 1, -1);
108 lnLDrawHistRate->SetDirectory(nullptr);
109 lnLDrawHistRate->GetXaxis()->SetTitle("-2LLH_{Pred Fluc, Draw}");
110 lnLDrawHistRate->GetYaxis()->SetTitle("-2LLH_{Data, Draw}");
111
112 //KS: This is silly as it assumes all samples uses same kinematics
113 lnLFlucHist_ProjectX = std::make_unique<TH2D>("lnLFlucHist_ProjectX", "lnLFlucHist_ProjectX", 50, 1, -1, 50, 1, -1);
114 lnLFlucHist_ProjectX->SetDirectory(nullptr);
115 lnLFlucHist_ProjectX->GetXaxis()->SetTitle(("-2LLH_{Draw Fluc, Draw} for " + SampleHandler->GetKinVarLabel(0, 0)).c_str());
116 lnLFlucHist_ProjectX->GetYaxis()->SetTitle(("-2LLH_{Data, Draw} for " + SampleHandler->GetKinVarLabel(0, 0)).c_str());
117
118 // Holds the hist of random number draws, only works for posterior predictive
120 {
121 RandomHist = std::make_unique<TH1D>("RandomHist", "RandomHist", 100, 0, nChainSteps);
122 RandomHist->SetDirectory(nullptr);
123 RandomHist->GetXaxis()->SetTitle("Step");
124 const double binwidth = nChainSteps/RandomHist->GetNbinsX();
125 std::stringstream ss;
126 ss << "Draws/" << binwidth;
127 RandomHist->GetYaxis()->SetTitle(ss.str().c_str());
128 RandomHist->SetLineWidth(2);
129 }
130 else RandomHist = nullptr;
131
132 for (int i = 0; i < nSamples; ++i)
133 {
134 DataHist[i] = nullptr;
135 DataHist_ProjectX[i] = nullptr;
136 DataHist_ProjectY[i] = nullptr;
137 NominalHist[i] = nullptr;
138
139 MeanHist[i] = nullptr;
140 if(DoBetaParam) MeanHistCorrected[i] = nullptr;
141 W2MeanHist[i] = nullptr;
142 W2ModeHist[i] = nullptr;
143 lnLHist_Mean[i] = nullptr;
144 lnLHist_Mode[i] = nullptr;
145 lnLHist_Mean_ProjectX[i] = nullptr;
146 lnLHist_Mean1D[i] = nullptr;
147 lnLHist_Mode1D[i] = nullptr;
148 lnLHist_Sample_DrawData[i] = nullptr;
149 lnLHist_Sample_DrawflucDraw[i] = nullptr;
150 lnLHist_Sample_PredflucDraw[i] = nullptr;
151 }//end loop over samples
152
153 DoByModePlots = false;
154 PosteriorHist_ByMode = nullptr;
155
156 nModelParams = 0;
157
158 Debug = 0;
159}
160
161// *******************
162// Destructor
164// *******************
165 Outputfile->cd();
166
167 //ROOT is weird and once you write TFile claim ownership of histograms. Best is to first delete histograms and then close file
168 Outputfile->Close();
169 delete Outputfile;
170
171 if(DoByModePlots)
172 {
173 for (int i = 0; i < nSamples; ++i)
174 {
175 if(DataHist[i] == nullptr) continue;
176 for (int j = 0; j < Modes->GetNModes()+1; j++)
177 {
178 for (int k = 1; k <= maxBins[i]; ++k)
179 {
180 if(PosteriorHist_ByMode[i][j][k] != nullptr) delete PosteriorHist_ByMode[i][j][k];
181 }
182 delete[] PosteriorHist_ByMode[i][j];
183 if(MeanHist_ByMode[i][j] != nullptr) delete MeanHist_ByMode[i][j];
184 }
185 delete[] PosteriorHist_ByMode[i];
186 }
187 delete[] PosteriorHist_ByMode;
188 }
189
190 for (int i = 0; i < nSamples; ++i)
191 {
192 if(DataHist[i] == nullptr) continue;
193 if(DataHist[i] != nullptr) delete DataHist[i];
194 if(NominalHist[i] != nullptr) delete NominalHist[i];
195 if(MeanHist[i] != nullptr) delete MeanHist[i];
196 if(ModeHist[i] != nullptr) delete ModeHist[i];
197 if(DoBetaParam && MeanHistCorrected[i] != nullptr) delete MeanHistCorrected[i];
198 if(W2MeanHist[i] != nullptr) delete W2MeanHist[i];
199 if(W2ModeHist[i] != nullptr) delete W2ModeHist[i];
200
201 if(ViolinHists_ProjectX[i] != nullptr) delete ViolinHists_ProjectX[i];
202 if(ViolinHists_ProjectY[i] != nullptr) delete ViolinHists_ProjectY[i];
203
204 if(lnLHist_Mean[i] != nullptr) delete lnLHist_Mean[i];
205 if(lnLHist_Mode[i] != nullptr) delete lnLHist_Mode[i];
206 if(lnLHist_Mean_ProjectX[i] != nullptr) delete lnLHist_Mean_ProjectX[i];
207 if(lnLHist_Mean1D[i] != nullptr) delete lnLHist_Mean1D[i];
208 if(lnLHist_Mode1D[i] != nullptr) delete lnLHist_Mode1D[i];
209 if(lnLHist_Sample_DrawData[i] != nullptr) delete lnLHist_Sample_DrawData[i];
212 }
213}
214
215// *******************
216// Check size of sample against size of vectors
217bool SampleSummary::CheckSamples(const int Length) {
218// *******************
219 bool ok = (nSamples == Length);
220 if (!ok) {
221 MACH3LOG_ERROR("Size of SampleVector input != number of defined samples");
222 MACH3LOG_ERROR("Size of SampleVector: {}", Length);
223 MACH3LOG_ERROR("Size of defined samples: {}", nSamples);
224 MACH3LOG_ERROR("Something has gone wrong with making the Samples");
225 throw MaCh3Exception(__FILE__ , __LINE__ );
226 }
227 return ok;
228}
229
230// *******************
231// Add a data histogram to the list (will have N_samples of these)
232// Since the data doesn't change with varying the MC
233void SampleSummary::AddData(std::vector<TH2Poly*> &Data) {
234// *******************
235 const int Length = int(Data.size());
236 // Check length of samples are OK
237 if (!CheckSamples(Length)) throw MaCh3Exception(__FILE__ , __LINE__ );
238 for (int i = 0; i < Length; ++i) {
239 if (Data[i] == nullptr) {
240 DataHist[i] = nullptr;
241 DataHist_ProjectX[i] = nullptr;
242 DataHist_ProjectY[i] = nullptr;
243 maxBins[i] = 0;
244 } else {
245 std::string classname = std::string(DataHist[i]->Class_Name());
246 if(classname == "TH2Poly")
247 {
248 DataHist[i] = static_cast<TH2Poly*>(Data[i]->Clone());
250 DataHist_ProjectX[i] = ProjectPoly(DataHist[i], true, i);
251 DataHist_ProjectY[i] = ProjectPoly(DataHist[i], false, i);
252 maxBins[i] = DataHist[i]->GetNumberOfBins();
253 } else {
254 MACH3LOG_ERROR("Somehow sample {} doesn't use TH2Poly", SampleHandler->GetSampleName(i));
255 MACH3LOG_ERROR("Right now I only support TH2Poly but I am ambitious piece of code and surely will have more support in the future");
256 throw MaCh3Exception(__FILE__ , __LINE__ );
257 }
258 }
259 }
260}
261
262// *******************
263// Add the nominal histograms to the list (will have N_samples of these)
264void SampleSummary::AddNominal(std::vector<TH2Poly*> &Nominal, std::vector<TH2Poly*> &NomW2) {
265// *******************
266 const int Length = int(Nominal.size());
267 if (!CheckSamples(Length)) throw MaCh3Exception(__FILE__ , __LINE__ );
268
269 //KS: ROOT is super annoying and you cannot use clone with openMP, hence we have another loop below
270 for (int i = 0; i < Length; ++i)
271 {
272 if (Nominal[i] == nullptr) {
273 NominalHist[i] = nullptr;
274 W2NomHist[i] = nullptr;
275 lnLHist_Mean[i] = nullptr;
276 lnLHist_Mode[i] = nullptr;
277 lnLHist_Mean_ProjectX[i] = nullptr;
278 MeanHist[i] = nullptr;
279 if(DoBetaParam) MeanHistCorrected[i] = nullptr;
280 ModeHist[i] = nullptr;
281 W2MeanHist[i] = nullptr;
282 W2ModeHist[i] = nullptr;
283 lnLHist_Sample_DrawData[i] = nullptr;
284 lnLHist_Sample_DrawflucDraw[i] = nullptr;
285 lnLHist_Sample_PredflucDraw[i] = nullptr;
286 // If not nullptr it indicates the selection was turned on, so initialise the privates
287 } else {
288 NominalHist[i] = static_cast<TH2Poly*>(Nominal[i]->Clone());
290 W2NomHist[i] = static_cast<TH2Poly*>(NomW2[i]->Clone());
291
292 lnLHist_Mean[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
293 lnLHist_Mean[i]->SetDirectory(nullptr);
294 lnLHist_Mode[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
295 lnLHist_Mode[i]->SetDirectory(nullptr);
296 lnLHist_Mean_ProjectX[i] = static_cast<TH1D*>(DataHist_ProjectX[i]->Clone());
297 MeanHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
298 if(DoBetaParam) MeanHistCorrected[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
299 ModeHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
300 W2MeanHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
301 W2ModeHist[i] = static_cast<TH2Poly*>(NominalHist[i]->Clone());
302 }
303 }
304
305 // Loop over the length of nominal and set the initial distributions up
306 //KS: Don't multithread, mostly due to fact that we initialise histograms
307 for (int i = 0; i < Length; ++i) {
308 // If NULL it indicates the selection was turned off, so initialise all the hists to NULL
309 if (Nominal[i] != nullptr)
310 {
311 std::string name = std::string(NominalHist[i]->GetName());
312 name = name.substr(0, name.find("_nom"));
313
314 PosteriorHist[i].resize(maxBins[i]+1);
315 w2Hist[i].resize(maxBins[i]+1);
316
317 if(DoBetaParam) BetaHist[i].resize(maxBins[i]+1);
318
319 for (int j = 0; j <= maxBins[i]; ++j)
320 {
321 PosteriorHist[i][j] = nullptr;
322 }
323 lnLHist_Mean[i]->SetNameTitle((name+"_MeanlnL").c_str(), (name+"_MeanlnL").c_str());
324 lnLHist_Mean[i]->Reset("");
325 lnLHist_Mean[i]->GetZaxis()->SetTitle("-2lnL_{sample} #times sign(MC-data)");
326
327 lnLHist_Mode[i]->SetNameTitle((name+"_ModelnL").c_str(), (name+"_ModelnL").c_str());
328 lnLHist_Mode[i]->Reset("");
329 lnLHist_Mode[i]->GetZaxis()->SetTitle("-2lnL_{sample} #times sign(MC-data)");
330
331 lnLHist_Mean_ProjectX[i]->SetNameTitle((name+"_MeanlnL_ProjectX").c_str(), (name+"_MeanlnL_ProjectX").c_str());
332 lnLHist_Mean_ProjectX[i]->Reset("");
333 lnLHist_Mean_ProjectX[i]->GetYaxis()->SetTitle("-2lnL_{sample} #times sign(MC-data)");
334
335 MeanHist[i]->SetNameTitle((name+"_mean").c_str(), (name+"_mean").c_str());
336 MeanHist[i]->Reset("");
337 MeanHist[i]->GetZaxis()->SetTitle("Mean");
338
339 if(DoBetaParam)
340 {
341 MeanHistCorrected[i]->SetNameTitle((name+"_mean_corrected").c_str(), (name+"_mean_corrected").c_str());
342 MeanHistCorrected[i]->Reset("");
343 MeanHistCorrected[i]->GetZaxis()->SetTitle("Mean");
344 }
345 std::vector<double> xbins;
346 std::vector<double> ybins;
347
348 SampleHandler->SetupBinning(M3::int_t(i), xbins, ybins);
349
350 //KS: Y axis is number of events to get estimate of maximal number we use integral
351 const int MaxBinning = doShapeOnly ? 1 : int(NoOverflowIntegral(NominalHist[i])/4);
352 ViolinHists_ProjectX[i] = new TH2D((name+"_Violin_ProjectX").c_str(), (name+"_Violin_ProjectX").c_str(), int(xbins.size()-1), &xbins[0] , 400, 0, MaxBinning);
353 ViolinHists_ProjectX[i]->GetYaxis()->SetTitle("Events");
354 ViolinHists_ProjectX[i]->GetXaxis()->SetTitle(std::string(NominalHist[i]->GetXaxis()->GetTitle()).c_str() );
355 ViolinHists_ProjectX[i]->SetDirectory(nullptr);
356
357 ViolinHists_ProjectY[i] = new TH2D((name+"_Violin_ProjectY").c_str(), (name+"_Violin_ProjectY").c_str(), int(ybins.size()-1), &ybins[0] , 400, 0, MaxBinning);
358 ViolinHists_ProjectY[i]->GetYaxis()->SetTitle("Events");
359 ViolinHists_ProjectY[i]->GetXaxis()->SetTitle(std::string(NominalHist[i]->GetYaxis()->GetTitle()).c_str());
360 ViolinHists_ProjectY[i]->SetDirectory(nullptr);
361
362 ModeHist[i]->SetNameTitle((name+"_mode").c_str(), (name+"_mode").c_str());
363 ModeHist[i]->Reset("");
364 ModeHist[i]->GetZaxis()->SetTitle("Mode");
365
366 W2MeanHist[i]->SetNameTitle((name+"_w2mean").c_str(), (name+"_w2mean").c_str());
367 W2MeanHist[i]->Reset("");
368 W2MeanHist[i]->GetZaxis()->SetTitle("W2Mean");
369
370 W2ModeHist[i]->SetNameTitle((name+"_w2mode").c_str(), (name+"_w2mode").c_str());
371 W2ModeHist[i]->Reset("");
372 W2ModeHist[i]->GetZaxis()->SetTitle("W2Mode");
373
374 // Declare the lnL histograms
375 lnLHist_Mean1D[i] = new TH1D((name+"_MeanlnL1D").c_str(),(name+"_MeanlnL1D").c_str(), 50, 1, -1);
376 lnLHist_Mean1D[i]->GetXaxis()->SetTitle("-2LLH (Data, Pred)");
377 lnLHist_Mean1D[i]->GetYaxis()->SetTitle("Counts");
378
379 lnLHist_Mode1D[i] = new TH1D((name+"_ModelnL1D").c_str(),(name+"_ModelnL1D").c_str(), 50, 1, -1);
380 lnLHist_Mode1D[i]->GetXaxis()->SetTitle("-2LLH (Data, Pred)");
381 lnLHist_Mode1D[i]->GetYaxis()->SetTitle("Counts");
382
383 lnLHist_Sample_DrawData[i] = new TH1D((name+"_lnLdrawdata").c_str(),(name+"_lnL").c_str(), 100, 1, -1);
384 lnLHist_Sample_DrawData[i]->GetXaxis()->SetTitle("-2LLH (Data, Draw)");
385 lnLHist_Sample_DrawData[i]->GetYaxis()->SetTitle("Counts");
386
387 lnLHist_Sample_DrawflucDraw[i] = new TH1D((name+"_lnLdrawfluc").c_str(),(name+"_lnL").c_str(), 100, 1, -1);
388 lnLHist_Sample_DrawflucDraw[i]->GetXaxis()->SetTitle("-2LLH (Draw Fluc, Draw)");
389 lnLHist_Sample_DrawflucDraw[i]->GetYaxis()->SetTitle("Counts");
390
391 lnLHist_Sample_PredflucDraw[i] = new TH1D((name+"_lnLpredfluc").c_str(),(name+"_lnL").c_str(), 100, 1, -1);
392 lnLHist_Sample_PredflucDraw[i]->GetXaxis()->SetTitle("-2LLH (Pred Fluc, Draw)");
393 lnLHist_Sample_PredflucDraw[i]->GetYaxis()->SetTitle("Counts");
394 }
395 }
396 //KS: Separate loop for thread safe reasons
397 for (int i = 0; i < Length; ++i)
398 {
399 //KS: We copy histograms so delete original
400 delete Nominal[i];
401 delete NomW2[i];
402 }
403}
404
405// *******************
406// Add a throw from the MCMC to the posterior predictive
407// The input here is nSamples long
408void SampleSummary::AddThrow(std::vector<TH2Poly*> &SampleVector, std::vector<TH2Poly*> &W2Vec, const double LLHPenalty, const double Weight, const int DrawNumber) {
409// *******************
410 nThrows++;
411 //KS: Only make sense for PosteriorPredictive
412 if( !isPriorPredictive )RandomHist->Fill(DrawNumber);
413
414 const int size = int(SampleVector.size());
415 if (!CheckSamples(size)) throw MaCh3Exception(__FILE__ , __LINE__ );
416
417 // Push back the throw
418 MCVector.push_back(SampleVector);
419 LLHPenaltyVector.push_back(LLHPenalty);
420 WeightVector.push_back(Weight);
421 W2MCVector.push_back(W2Vec);
422
423 // Initialise the posterior hist
424 if (first_pass)
425 {
426 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
427 {
428 const int nXBins = 500;
429 //Initialise TH1D which corresponds to each bin in the sample's th2poly
430 std::string name = std::string(SampleVector[SampleNum]->GetName());
431 for (int i = 1; i <= maxBins[SampleNum]; ++i)
432 {
433 //Get PolyBin
434 TH2PolyBin* bin = static_cast<TH2PolyBin*>(SampleVector[SampleNum]->GetBins()->At(i-1));
435
436 // Just make a little fancy name
437 std::stringstream ss2;
438 ss2 << name << "_";
439 ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
440 ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
441
442 PosteriorHist[SampleNum][i] = std::make_unique<TH1D>(ss2.str().c_str(), ss2.str().c_str(),nXBins, 1, -1);
443 PosteriorHist[SampleNum][i]->SetDirectory(nullptr);
444 w2Hist[SampleNum][i] = std::make_unique<TH1D>(("w2_"+ss2.str()).c_str(), ("w2_"+ss2.str()).c_str(),nXBins, 1, -1);
445 w2Hist[SampleNum][i]->SetDirectory(nullptr);
446 if(DoBetaParam)
447 {
448 std::string betaName = "#beta_param_";
449 BetaHist[SampleNum][i] = std::make_unique<TH1D>((betaName + ss2.str()).c_str(), (betaName + ss2.str()).c_str(), 70, 1, -1);
450 BetaHist[SampleNum][i]->SetDirectory(nullptr);
451 BetaHist[SampleNum][i]->GetXaxis()->SetTitle("#beta parameter value");
452 BetaHist[SampleNum][i]->GetYaxis()->SetTitle("Counts");
453 }
454 } //end loop over bins
455 }//end loop over samples
456 }
457 first_pass = false;
458
459 // Loop over the samples
460 #ifdef MULTITHREAD
461 #pragma omp parallel for
462 #endif
463 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
464 {
465 if (SampleVector[SampleNum] == nullptr) continue;
466 if(doShapeOnly) NormaliseTH2Poly(SampleVector[SampleNum]);
467 // Loop over the distribution and fill the prior/posterior predictive
468 for (int i = 1; i <= maxBins[SampleNum]; ++i) {
469 const double Content = SampleVector[SampleNum]->GetBinContent(i);
470 PosteriorHist[SampleNum][i]->Fill(Content, Weight);
471 const double w2 = W2Vec[SampleNum]->GetBinContent(i);
472 w2Hist[SampleNum][i]->Fill(w2, Weight);
473 if(DoBetaParam)
474 {
475 const double data = DataHist[SampleNum]->GetBinContent(i);
476 const double BetaParam = GetBetaParameter(data, Content, w2, likelihood);
477 BetaHist[SampleNum][i]->Fill(BetaParam, Weight);
478 }
479 } // end bin loop
480 } // end samples loop
481} // end AddThrow
482
483// *******************
484// Add a throw from the MCMC to the posterior predictive
485// The input here is has dimension [nsample][nMaCh3Modes]
486void SampleSummary::AddThrowByMode(std::vector<std::vector<TH2Poly*>> &SampleVector_ByMode) {
487// *******************
488 MCVectorByMode.push_back(SampleVector_ByMode);
489
490 //KS: This means this is first time
491 if(!DoByModePlots)
492 {
493 MACH3LOG_INFO("Turning reaction breadkwon mode, brum brum");
494 PosteriorHist_ByMode = new TH1D***[nSamples];
496 for (int SampleNum = 0; SampleNum < nSamples; SampleNum++)
497 {
498 if (DataHist[SampleNum] == nullptr) continue;
499
500 PosteriorHist_ByMode[SampleNum] = new TH1D**[Modes->GetNModes()+1];
501 MeanHist_ByMode[SampleNum].resize(Modes->GetNModes()+1);
502 for (int j = 0; j < Modes->GetNModes()+1; j++)
503 {
504 PosteriorHist_ByMode[SampleNum][j] = new TH1D*[maxBins[SampleNum]+1];
505 constexpr int nXBins = 500;
506
507 std::string name = std::string(NominalHist[SampleNum]->GetName());
508 name = name.substr(0, name.find("_nom"));
509 name = name + "_"+Modes->GetMaCh3ModeName(j);
510
511 for (int i = 1; i <= maxBins[SampleNum]; i++)
512 {
513 //Get PolyBin
514 TH2PolyBin* bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(i-1));
515
516 // Just make a little fancy name
517 std::stringstream ss2;
518 ss2 << name << "_";
519 ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
520 ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
521
522 //Initialise TH1D which corresponds to each bin in the sample's th2poly
523 PosteriorHist_ByMode[SampleNum][j][i] = new TH1D((name+ss2.str()).c_str(),(name+ss2.str()).c_str(),nXBins, 1, -1);
524 }
525 MeanHist_ByMode[SampleNum][j] = static_cast<TH2Poly*>(NominalHist[SampleNum]->Clone());
526 MeanHist_ByMode[SampleNum][j]->SetNameTitle((name+"_mean").c_str(), (name+"_mean").c_str());
527 MeanHist_ByMode[SampleNum][j]->Reset("");
528 MeanHist_ByMode[SampleNum][j]->GetZaxis()->SetTitle("Mean");
529 }
530 }
531 }
532 DoByModePlots = true;
533 // Loop over the sameples
534 #ifdef MULTITHREAD
535 #pragma omp parallel for
536 #endif
537 for (int SampleNum = 0; SampleNum < nSamples; SampleNum++)
538 {
539 if (DataHist[SampleNum] == nullptr) continue;
540
541 for (int j = 0; j < Modes->GetNModes()+1; j++)
542 {
543 if(doShapeOnly) NormaliseTH2Poly(SampleVector_ByMode[SampleNum][j]);
544 // Loop over the distribution and fill the prior/posterior predictive
545 for (int i = 1; i <= maxBins[SampleNum]; ++i)
546 {
547 const double Content = SampleVector_ByMode[SampleNum][j]->GetBinContent(i);
548 const int Entries = int(PosteriorHist_ByMode[SampleNum][j][i]->GetEntries());
549 PosteriorHist_ByMode[SampleNum][j][i]->Fill(Content, WeightVector[Entries]);
550 }
551 }
552 }
553} // end AddThrowByMode
554
555// **********************
557// **********************
558
559 // Make the output file (MakePosterioPredictive call writes to this)
560 std::string TempString = OutputName;
561 TempString.replace(TempString.find(".root"), 5, std::string("_procsW2.root"));
562 Outputfile = new TFile(TempString.c_str(), "RECREATE");
563
564 // The array of doubles we write to the TTree
565 // Data vs Draw
566 llh_data_draw.resize(nSamples);
567 // Fluctuated Draw vs Draw
569 // Fluctuated Predicitve vs Draw
571
572 // Data vs Draw using Rate
574 // Data vs Fluctuated Predictive using Rate
576
577 // Data vs Fluctuated Draw
579 // Data vs Fluctuated Predictive
581 // Draw vs Predictive
582 llh_draw_pred.resize(nSamples);
583 // Fluctuated Draw vs Predictive
585 // Fluctuated Draw vs Fluctuated Predictive
587
588 // Fluctuated Predictive vs Predictive
590 // Fluctuated Data vs Draw
592
593 // Data vs Draw for 1D projection
596
597 // The output tree we're going to write to
598 OutputTree = new TTree("LLH_draws", "LLH_draws");
599 SampleNames.resize(nSamples);
600 // Loop over the samples and set the addresses of the variables to write to file
601 for (int i = 0; i < nSamples; ++i)
602 {
603 // Get the name
604 std::string SampleName = SampleHandler->GetSampleName(i);
605 // Strip out spaces
606 while (SampleName.find(" ") != std::string::npos) {
607 SampleName.replace(SampleName.find(" "), 1, std::string("_"));
608 }
609 SampleNames[i] = SampleName;
610 //CW: Also strip out - signs because it messes up TBranches
611 while (SampleName.find("-") != std::string::npos) {
612 SampleName.replace(SampleName.find("-"), 1, std::string("_"));
613 }
614// All LLH below are used for actual p-value calculations
615 OutputTree->Branch((SampleName+"_data_draw").c_str(), &llh_data_draw[i]);
616 OutputTree->Branch((SampleName+"_drawfluc_draw").c_str(), &llh_drawfluc_draw[i]);
617 OutputTree->Branch((SampleName+"_predfluc_draw").c_str(), &llh_predfluc_draw[i]);
618
619// All LLH below are used for actual p-value calculations however using rate only
620 OutputTree->Branch((SampleName+"_rate_data_draw").c_str(), &llh_rate_data_draw[i]);
621 OutputTree->Branch((SampleName+"_rate_predfluc_draw").c_str(), &llh_rate_predfluc_draw[i]);
622
623// All LLH below are for validation reason but not used for final P-Value
624 OutputTree->Branch((SampleName+"_data_drawfluc").c_str(), &llh_data_drawfluc[i]);
625 OutputTree->Branch((SampleName+"_data_predfluc").c_str(), &llh_data_predfluc[i]);
626 OutputTree->Branch((SampleName+"_draw_pred").c_str(), &llh_draw_pred[i]);
627 OutputTree->Branch((SampleName+"_drawfluc_pred").c_str(), &llh_drawfluc_pred[i]);
628 OutputTree->Branch((SampleName+"_drawfluc_predfluc").c_str(), &llh_drawfluc_predfluc[i]);
629 OutputTree->Branch((SampleName+"_predfluc_pred").c_str(), &llh_predfluc_pred[i]);
630 OutputTree->Branch((SampleName+"_datafluc_draw").c_str(), &llh_datafluc_draw[i]);
631
632// All LLH below are used for calcauting P-Value but using 1D projections
633 OutputTree->Branch((SampleName+"_data_draw_ProjectX").c_str(), &llh_data_draw_ProjectX[i]);
634 OutputTree->Branch((SampleName+"_drawfluc_draw_ProjectX").c_str(), &llh_drawfluc_draw_ProjectX[i]);
635 }
636//All LLH below are used for actual p-value calculations
637 OutputTree->Branch("LLH_Penalty", &llh_penalty);
638 OutputTree->Branch("Total_LLH_Data_Draw", &total_llh_data_draw);
639 OutputTree->Branch("Total_LLH_DrawFluc_Draw", &total_llh_drawfluc_draw);
640 OutputTree->Branch("Total_LLH_PredFluc_Draw", &total_llh_predfluc_draw);
641
642// All LLH below are used for actual p-value calculations however using rate only
643 OutputTree->Branch("Total_LLH_Rate_PredFluc_Draw", &total_llh_rate_predfluc_draw);
644
645//All LLH below are for validation reason but not used for final P-Value
646 OutputTree->Branch("Total_LLH_Data_DrawFluc", &total_llh_data_drawfluc);
647 OutputTree->Branch("Total_LLH_Data_PredFluc", &total_llh_data_predfluc);
648 OutputTree->Branch("Total_LLH_Draw_Pred", &total_llh_draw_pred);
649 OutputTree->Branch("Total_LLH_DrawFluc_Pred", &total_llh_drawfluc_pred);
650 OutputTree->Branch("Total_LLH_DrawFluc_PredFluc", &total_llh_drawfluc_predfluc);
651 OutputTree->Branch("Total_LLH_PredFluc_Pred", &total_llh_predfluc_pred);
652 OutputTree->Branch("Total_LLH_DataFluc_Draw", &total_llh_datafluc_draw);
653
654//All LLH below are used for calcauting P-Value but 1D projections
655 OutputTree->Branch("total_llh_data_draw_ProjectX", &total_llh_data_draw_ProjectX);
656 OutputTree->Branch("total_llh_drawfluc_draw_ProjectX", &total_llh_drawfluc_draw_ProjectX);
657
658 Outputfile->cd();
659 Dir.resize(nSamples);
660 for (int i = 0; i < nSamples; ++i)
661 {
662 // Make a new directory
663 Dir[i] = Outputfile->mkdir((SampleNames[i]).c_str());
664 }
665}
666
667// *******************
668// Write the contents to the file
670// *******************
671 // Prepare the output tree
673
674 MACH3LOG_INFO("Summarising {} throws...", nThrows);
675 // After all the throws are added finalise the sample
676 TStopwatch timer;
677 timer.Start();
679 timer.Stop();
680 MACH3LOG_INFO("Made Prior/Posterior Predictive, it took {:.2f}s, now writing...", timer.RealTime());
681
682 // Studying information criterion
684
685 OutputTree->Write();
686
687 // Make the various distributions
688 lnLHist->Write();
689 lnLHist_drawfluc->Write();
690 lnLHist_drawflucdraw->Write();
691 lnLHist_drawdata->Write();
692 lnLDrawHist->Write();
693 lnLFlucHist->Write();
694 lnLDrawHistRate->Write();
695 //KS: Only available for Posterior Predictive
696 if(!isPriorPredictive) RandomHist->Write();
697
698 lnLFlucHist_ProjectX->Write();
699
700 // Loop over each sample and write to file
701 //KS: Multithreading is tempting here but we also write to ROOT file, separating all LLH and poly projections from write could work well
702 for (int i = 0; i < nSamples; ++i)
703 {
704 // Skip the null histograms
705 if (DataHist[i] == nullptr || NoOverflowIntegral(DataHist[i]) == 0) continue;
706 Dir[i]->cd();
707
708 // Make the data/MC ratio histogram
709 TH2Poly *RatioHistMean = RatioPolys(DataHist[i], MeanHist[i]);
710 RatioHistMean->GetZaxis()->SetTitle("Data/Mean");
711 TH2Poly *RatioHistMode = RatioPolys(DataHist[i], ModeHist[i]);
712 RatioHistMode->GetZaxis()->SetTitle("Data/Mode");
713 TH2Poly *RatioHistNom = RatioPolys(DataHist[i], NominalHist[i]);
714 RatioHistNom->GetZaxis()->SetTitle("Data/Nom");
715
716 // And the normalised data histogram
717 TH2Poly *DataNormHist = NormalisePoly(DataHist[i]);
718 // Last true refers to if project along x or y
719 TH2Poly *MeanNormHist = NormalisePoly(MeanHist[i]);
720 TH2Poly *ModeNormHist = NormalisePoly(ModeHist[i]);
721 TH1D *MeanProjectX = ProjectPoly(MeanHist[i], true, i, true);
722 TH1D *MeanProjectY = ProjectPoly(MeanHist[i], false, i, true);
723 TH1D *ModeProjectX = ProjectPoly(ModeHist[i], true, i, true);
724 TH1D *ModeProjectY = ProjectPoly(ModeHist[i], false, i, true);
725
726 TH1D *MeanHistCorrectedProjectX = nullptr;
727 if(DoBetaParam) MeanHistCorrectedProjectX = ProjectPoly(MeanHistCorrected[i], true, i, true);
728 TH1D *MeanHistCorrectedProjectY = nullptr;
729 if(DoBetaParam) MeanHistCorrectedProjectY = ProjectPoly(MeanHistCorrected[i], false, i, true);
730
731 TH1D *W2MeanProjectX = ProjectPoly(W2MeanHist[i], true, i);
732 TH1D *W2MeanProjectY = ProjectPoly(W2MeanHist[i], false, i);
733 TH1D *W2ModeProjectX = ProjectPoly(W2ModeHist[i], true, i);
734 TH1D *W2ModeProjectY = ProjectPoly(W2ModeHist[i], false, i);
735
736 TH2Poly *NomNormHist = NormalisePoly(NominalHist[i]);
737 TH1D *NomProjectX = ProjectPoly(NominalHist[i], true, i);
738 TH1D *NomProjectY = ProjectPoly(NominalHist[i], false, i);
739
740 TH1D *W2NomProjectX = ProjectPoly(W2NomHist[i], true, i);
741 TH1D *W2NomProjectY = ProjectPoly(W2NomHist[i], false, i);
742
743 // Same for the TH2Ds
747
748 // Calculate the log likelihood for the 1D dists
749 // Sets the title of the second TH1D to the -2LLH
750 CalcLLH(DataHist_ProjectX[i], NomProjectX, W2NomProjectX);
751 CalcLLH(DataHist_ProjectX[i], MeanProjectX, W2MeanProjectX);
752 CalcLLH(DataHist_ProjectX[i], ModeProjectX, W2ModeProjectX);
753 CalcLLH(DataHist_ProjectY[i], NomProjectY, W2NomProjectY);
754 CalcLLH(DataHist_ProjectY[i], MeanProjectY, W2MeanProjectY);
755 CalcLLH(DataHist_ProjectY[i], ModeProjectY, W2ModeProjectY);
756
757 std::string SampleName = SampleNames[i];
758 // Also strip out - signs because it messes up TBranches
759 while (SampleName.find("-") != std::string::npos) {
760 SampleName.replace(SampleName.find("-"), 1, std::string("_"));
761 }
762 OutputTree->Draw((SampleName+"_data_draw:"+SampleName+"_drawfluc_draw>>htemp").c_str());
763 TH2D *TempHistogram = static_cast<TH2D*>(gDirectory->Get("htemp")->Clone());
764 TempHistogram->GetXaxis()->SetTitle("-2LLH(Draw Fluc, Draw)");
765 TempHistogram->GetYaxis()->SetTitle("-2LLH(Data, Draw)");
766 TempHistogram->SetNameTitle((SampleNames[i]+"_drawfluc_draw").c_str(), (SampleNames[i]+"_drawfluc_draw").c_str());
767 Get2DBayesianpValue(TempHistogram);
768 TempHistogram->Write();
769 delete TempHistogram;
770
771 // Also write the 2D histograms for the p-value
772 OutputTree->Draw((SampleName+"_data_draw:"+SampleName+"_predfluc_draw>>htemp2").c_str());
773 TH2D *TempHistogram2 = static_cast<TH2D*>(gDirectory->Get("htemp2")->Clone());
774 TempHistogram2->GetXaxis()->SetTitle("-2LLH(Pred Fluc, Draw)");
775 TempHistogram2->GetYaxis()->SetTitle("-2LLH(Data, Draw)");
776 TempHistogram2->SetNameTitle((SampleNames[i]+"_predfluc_draw").c_str(), (SampleNames[i]+"_predfluc_draw").c_str());
777 Get2DBayesianpValue(TempHistogram2);
778 TempHistogram2->Write();
779 delete TempHistogram2;
780
781 // finally p-value for 1D projection
782 OutputTree->Draw((SampleName+"_rate_data_draw:"+SampleName+"_rate_predfluc_draw>>htemp3").c_str());
783 TH2D *TempHistogram3 = static_cast<TH2D*>(gDirectory->Get("htemp3")->Clone());
784 TempHistogram3->GetXaxis()->SetTitle("-2LLH(Pred Fluc, Draw)");
785 TempHistogram3->GetYaxis()->SetTitle("-2LLH(Data, Draw)");
786 TempHistogram3->SetNameTitle((SampleNames[i]+"_rate_predfluc_draw").c_str(), (SampleNames[i]+"_rate_predfluc_draw").c_str());
787 Get2DBayesianpValue(TempHistogram3);
788 TempHistogram3->Write();
789 delete TempHistogram3;
790
791 // finally p-value for 1D projection
792 OutputTree->Draw((SampleName+"_data_draw_ProjectX:"+SampleName+"_drawfluc_draw_ProjectX>>htemp4").c_str());
793 TH2D *TempHistogram4 = static_cast<TH2D*>(gDirectory->Get("htemp4")->Clone());
794 TempHistogram4->GetXaxis()->SetTitle(("-2LLH_{Draw Fluc, Draw} for " + SampleHandler->GetKinVarLabel(i, 0)).c_str());
795 TempHistogram4->GetYaxis()->SetTitle(("-2LLH_{Data, Draw} for " + SampleHandler->GetKinVarLabel(i, 0)).c_str());
796 TempHistogram4->SetNameTitle((SampleNames[i]+"_drawfluc_draw_ProjectX").c_str(), (SampleNames[i]+"_drawfluc_draw_ProjectX").c_str());
797 Get2DBayesianpValue(TempHistogram4);
798 TempHistogram4->Write();
799 delete TempHistogram4;
800
801 // Write the Histograms to each folder
802 DataHist[i]->Write();
803 NominalHist[i]->Write();
804 MeanHist[i]->Write();
805 ModeHist[i]->Write();
806 RatioHistMean->Write();
807 RatioHistMode->Write();
808 RatioHistNom->Write();
809 if(DoBetaParam) MeanHistCorrected[i]->Write();
810
811 W2NomHist[i]->Write();
812 W2MeanHist[i]->Write();
813 W2ModeHist[i]->Write();
814
815 DataNormHist->Write();
816 NomNormHist->Write();
817 MeanNormHist->Write();
818 ModeNormHist->Write();
819
820 DataHist_ProjectX[i]->Write();
821 NomProjectX->Write();
822 MeanProjectX->Write();
823 ModeProjectX->Write();
824 if(DoBetaParam) MeanHistCorrectedProjectX->Write();
825 ViolinHists_ProjectX[i]->Write();
826
827 DataHist_ProjectY[i]->Write();
828 NomProjectY->Write();
829 MeanProjectY->Write();
830 ModeProjectY->Write();
831 if(DoBetaParam) MeanHistCorrectedProjectY->Write();
832 ViolinHists_ProjectY[i]->Write();
833
834 W2NomProjectX->Write();
835 W2MeanProjectX->Write();
836 W2ModeProjectX->Write();
837
838 W2NomProjectY->Write();
839 W2MeanProjectY->Write();
840 W2ModeProjectY->Write();
841
842 //KS: This will dump lots of hists, use it only for debugging
843 if(Debug > 0)
844 {
845 TDirectory* DebugDir = Dir[i]->mkdir("Debug");
846 DebugDir->cd();
847 for (int b = 1; b <= maxBins[i]; ++b)
848 {
849 PosteriorHist[i][b]->Write();
850 std::string Title = PosteriorHist[i][b]->GetName();
851
852 auto TempLine = std::make_unique<TLine>(NominalHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMinimum(),
853 NominalHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMaximum());
854 TempLine->SetLineColor(kRed);
855 TempLine->SetLineWidth(2);
856
857 auto TempLineData = std::make_unique<TLine>(DataHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMinimum(),
858 DataHist[i]->GetBinContent(b), PosteriorHist[i][b]->GetMaximum());
859 TempLineData->SetLineColor(kGreen);
860 TempLineData->SetLineWidth(2);
861
862 // Also fit a Gaussian because why not?
863 TF1 *Fitter = new TF1("Fit", "gaus", PosteriorHist[i][b]->GetBinLowEdge(1), PosteriorHist[i][b]->GetBinLowEdge(PosteriorHist[i][b]->GetNbinsX()+1));
864 PosteriorHist[i][b]->Fit(Fitter, "RQ");
865 Fitter->SetLineColor(kRed-5);
866
867 auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
868 Legend->SetFillColor(0);
869 Legend->SetFillStyle(0);
870 Legend->SetLineWidth(0);
871 Legend->SetLineColor(0);
872 Legend->AddEntry(TempLineData.get(), Form("Data #mu=%.2f", DataHist[i]->GetBinContent(b)), "l");
873 Legend->AddEntry(TempLine.get(), Form("Prior #mu=%.2f", NominalHist[i]->GetBinContent(b)), "l");
874 Legend->AddEntry(PosteriorHist[i][b].get(), Form("Post, #mu=%.2f#pm%.2f", PosteriorHist[i][b]->GetMean(), PosteriorHist[i][b]->GetRMS()), "l");
875 Legend->AddEntry(Fitter, Form("Gauss, #mu=%.2f#pm%.2f", Fitter->GetParameter(1), Fitter->GetParameter(2)), "l");
876 std::string TempTitle = std::string(PosteriorHist[i][b]->GetName());
877
878 TempTitle += "_canv";
879 TCanvas *TempCanvas = new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
880 TempCanvas->SetGridx();
881 TempCanvas->SetGridy();
882 TempCanvas->SetRightMargin(0.03);
883 TempCanvas->SetBottomMargin(0.08);
884 TempCanvas->SetLeftMargin(0.10);
885 TempCanvas->SetTopMargin(0.06);
886 TempCanvas->cd();
887 PosteriorHist[i][b]->Draw();
888 TempLine->Draw("same");
889 TempLineData->Draw("same");
890 Fitter->Draw("same");
891 Legend->Draw("same");
892 TempCanvas->Write();
893
894 delete TempCanvas;
895 delete Fitter;
896 //This isn't useful check only in desperation
897 if(Debug > 1) w2Hist[i][b]->Write();
898 }
899 DebugDir->Close();
900 delete DebugDir;
901 Dir[i]->cd();
902 }
903 lnLHist_Mean[i]->Write();
904 lnLHist_Mode[i]->Write();
905
906 lnLHist_Mean_ProjectX[i]->Write();
907
908 lnLHist_Mean1D[i]->Write();
909 lnLHist_Mode1D[i]->Write();
910
912 lnLHist_Sample_DrawData[i]->Write();
914 lnLHist_Sample_DrawflucDraw[i]->Write();
916 lnLHist_Sample_PredflucDraw[i]->Write();
917
918 if(DoByModePlots)
919 {
920 for (int j = 0; j < Modes->GetNModes()+1; ++j)
921 {
922 MeanHist_ByMode[i][j]->Write();
923 TH1D *MeanProjectX_ByMode = ProjectPoly(MeanHist_ByMode[i][j], true, i, true);
924 TH1D *MeanProjectY_ByMode = ProjectPoly(MeanHist_ByMode[i][j], false, i, true);
925 MeanProjectX_ByMode->Write();
926 MeanProjectY_ByMode->Write();
927 //KS: This will dump lots of hists, use it only for debugging
928 if(Debug > 0)
929 {
930 for (int b = 1; b <= maxBins[i]; ++b)
931 {
932 PosteriorHist_ByMode[i][j][b]->Write();
933 }
934 }
935 delete MeanProjectX_ByMode;
936 delete MeanProjectY_ByMode;
937 } // End loop over bins
938 }
939 // Delete temporary objects
940 delete RatioHistMean;
941 delete RatioHistMode;
942 delete RatioHistNom;
943
944 delete DataNormHist;
945 delete MeanNormHist;
946 delete ModeNormHist;
947 delete NomNormHist;
948
949 delete DataHist_ProjectX[i];
950 delete MeanProjectX;
951 delete ModeProjectX;
952 if(DoBetaParam) delete MeanHistCorrectedProjectX;
953 delete NomProjectX;
954
955 delete DataHist_ProjectY[i];
956 delete MeanProjectY;
957 delete ModeProjectY;
958 if(DoBetaParam) delete MeanHistCorrectedProjectY;
959 delete NomProjectY;
960
961 delete W2NomProjectX;
962 delete W2MeanProjectX;
963 delete W2ModeProjectX;
964
965 delete W2NomProjectY;
966 delete W2MeanProjectY;
967 delete W2ModeProjectY;
968 MACH3LOG_INFO("");
969 } //end loop over samples
971
973 MACH3LOG_INFO("Wrote to {}", Outputfile->GetName());
974}
975
976// *******************
977// Make the posterior predictive distributions: fit Poisson etc
979// *******************
980 // First make the projection on the z axis of the TH3D* for every pmu cosmu bin
981 double llh_total_temp = 0.0;
982
983 // Loop over the samples
984 #ifdef MULTITHREAD
985 #pragma omp parallel for reduction(+:llh_total_temp)
986 #endif
987 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
988 {
989 // Skip disabled samples
990 if (DataHist[SampleNum] == nullptr || NoOverflowIntegral(DataHist[SampleNum]) == 0) continue;
991
992 // Count the -2LLH for each histogram
993 double negLogL_Mean = 0.0;
994 double negLogL_Mode = 0.0;
995
996 // Loop over each pmu cosmu bin
997 for (int j = 1; j < maxBins[SampleNum]+1; ++j)
998 {
999 TH1D *Projection = PosteriorHist[SampleNum][j].get();
1000 TH1D *W2Projection = w2Hist[SampleNum][j].get();
1001
1002 // Data content for the j,kth bin
1003 const double nData = DataHist[SampleNum]->GetBinContent(j);
1004
1005 // Get the mean for this projection for all the samples
1006 // This is the mean prediction for this given j,k bin
1007 const double nMean = Projection->GetMean();
1008 const double nMeanError = Projection->GetRMS();
1009 const double nMode = Projection->GetBinCenter(Projection->GetMaximumBin());
1010 const double nModeError = GetModeError(Projection);
1011
1012 const double nW2Mean = W2Projection->GetMean();
1013 const double nW2Mode = W2Projection->GetBinCenter(W2Projection->GetMaximumBin());
1014
1015 double TempLLH_Mean = 0.0;
1016 double TempLLH_Mode = 0.0;
1017
1018 //KS:Get LLH contribution getTestStatLLH can calculate Barlow Beeston/IceCube or Poisson
1019 TempLLH_Mean = SampleHandler->GetTestStatLLH(nData, nMean, nW2Mean);
1020 TempLLH_Mode = SampleHandler->GetTestStatLLH(nData, nMode, nW2Mode);
1021
1022 // Increment -2LLH
1023 //KS: do times 2 because banff reports chi2
1024 negLogL_Mean += 2*TempLLH_Mean;
1025 negLogL_Mode += 2*TempLLH_Mode;
1026
1027 // Set the content and error to the mean in the bin
1028 MeanHist[SampleNum]->SetBinContent(j, MeanHist[SampleNum]->GetBinContent(j)+nMean);
1029 // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
1030 MeanHist[SampleNum]->SetBinError(j, nMeanError);
1031
1032 if(DoBetaParam)
1033 {
1034 TH1D *BetaTemp = BetaHist[SampleNum][j].get();
1035 const double nBetaMean = BetaTemp->GetMean();
1036 const double nBetaMeanError = BetaTemp->GetRMS();
1037 //KS: Here we modify predictions by beta parameter from Barlow-Beeston
1038 MeanHistCorrected[SampleNum]->SetBinContent(j, MeanHistCorrected[SampleNum]->GetBinContent(j)+nMean*nBetaMean);
1039 //KS: Use error propagation to calcuate error
1040 const double ErrorTemp = std::sqrt( (nBetaMean*nMeanError) * (nBetaMean*nMeanError) + (nMean*nBetaMeanError) * (nMean*nBetaMeanError));
1041 // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
1042 MeanHistCorrected[SampleNum]->SetBinError(j, ErrorTemp);
1043 }
1044 // Set the content to the mode in the bin
1045 ModeHist[SampleNum]->SetBinContent(j, ModeHist[SampleNum]->GetBinContent(j)+nMode);
1046 // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
1047 ModeHist[SampleNum]->SetBinError(j, nModeError);
1048 // Set the content to the mean in the bin
1049 W2MeanHist[SampleNum]->SetBinContent(j, W2MeanHist[SampleNum]->GetBinContent(j)+nW2Mean);
1050 // Set the content to the mode in the bin
1051 W2ModeHist[SampleNum]->SetBinContent(j, W2ModeHist[SampleNum]->GetBinContent(j)+nW2Mode);
1052
1053 // Set the mean and average LLH for this given bin
1054 // Can use these hists to see where the largest -2LLH hists come from
1055 lnLHist_Mean[SampleNum]->SetBinContent(j, 2.0*TempLLH_Mean);
1056 lnLHist_Mode[SampleNum]->SetBinContent(j, 2.0*TempLLH_Mode);
1057
1058 lnLHist_Mean1D[SampleNum]->Fill(2.0*TempLLH_Mean);
1059 lnLHist_Mode1D[SampleNum]->Fill(2.0*TempLLH_Mode);
1060 } // End loop over bins
1061 if(DoByModePlots)
1062 {
1063 for (int j = 0; j < Modes->GetNModes()+1; j++)
1064 {
1065 // Loop over each pmu cosmu bin
1066 for (int i = 1; i < maxBins[SampleNum]+1; ++i)
1067 {
1068 // Make the posterior/prior predictive projection on z
1069 // The z axis of Predictive is the bin content
1070 // Essentially zooming in on one bin and looking at the mean and mode of that bin
1071 TH1D *Projection = PosteriorHist_ByMode[SampleNum][j][i];
1072
1073 // Get the mean for this projection for all the samples
1074 const double nMean = Projection->GetMean();
1075 const double nMeanError = Projection->GetRMS();
1076
1077 // Set the content and error to the mean in the bin
1078 MeanHist_ByMode[SampleNum][j]->SetBinContent(i, MeanHist_ByMode[SampleNum][j]->GetBinContent(i)+nMean);
1079 // KS: This -1 is only needed for root older than 6.18 for more see https://t2k-experiment.slack.com/archives/CU9CBG6NS/p1714551365661589
1080 MeanHist_ByMode[SampleNum][j]->SetBinError(i, nMeanError);
1081 } // End loop over bins
1082 }
1083 }
1084 llh_total_temp += negLogL_Mean;
1085 } // End loop over samples
1086
1087 // This is not multithreaded as due to ProjectPoly it is not safe
1088 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1089 {
1090 // Skip disabled samples
1091 if (DataHist[SampleNum] == nullptr || NoOverflowIntegral(DataHist[SampleNum]) == 0) continue;
1092
1093 //KS:: Might consider caching it as we use it once agian much later
1094 TH1D *MeanProjectX = ProjectPoly(MeanHist[SampleNum], true, SampleNum, true);
1095 TH1D *W2MeanProjectX = ProjectPoly(W2MeanHist[SampleNum], true, SampleNum);
1096 // Loop over each pmu bin for 1D projection
1097 for (int j = 1; j <= lnLHist_Mean_ProjectX[SampleNum]->GetXaxis()->GetNbins(); ++j)
1098 {
1099 // Data content for the j,kth bin
1100 const double nData = DataHist_ProjectX[SampleNum]->GetBinContent(j);
1101 const double nMean = MeanProjectX->GetBinContent(j);
1102 const double nW2Mean = W2MeanProjectX->GetBinContent(j);
1103
1104 double TempLLH_Mean = 0.0;
1105 TempLLH_Mean = SampleHandler->GetTestStatLLH(nData, nMean, nW2Mean);
1106
1107 //KS: do times 2 because banff reports chi2
1108 lnLHist_Mean_ProjectX[SampleNum]->SetBinContent(j, 2.0*TempLLH_Mean);
1109 }// End loop over bins
1110
1111 delete MeanProjectX;
1112 delete W2MeanProjectX;
1113 } // End loop over samples
1114
1115 llh_total = llh_total_temp;
1116 // Now we have our posterior predictive histogram and it's LLH
1117 MACH3LOG_INFO("Prior/Posterior predictive LLH mean (sample only) = {:.2f}", llh_total);
1118 std::stringstream ss;
1119 ss << llh_total;
1120 lnLHist->SetTitle((std::string(lnLHist->GetTitle())+"_"+ss.str()).c_str());
1121
1122 // Now make the fluctuated hists of the MeanHist and ModeHist
1123 MakeChi2Hists();
1124
1125 // Get the 1D LLH dists
1126 MakeCutLLH();
1127} // End MakePredictive() function
1128
1129// *******************
1130// Make the fluctuated histograms (2D and 1D) for the chi2s
1131// Essentially taking the MCMC draws and calculating their LLH to the Posterior predictive distribution
1132// And additionally taking the data histogram and calculating the LLH to the predictive distribution
1133// Additionally we calculate the chi2 of the draws (fluctuated) of the MC with the prior/posterior predictive and plot it vs the chi2 from the draws of MCMC and the data
1135// *******************
1136 MACH3LOG_INFO("Making the chi2 histograms...");
1137 // Have this to signify if we're doing the first pass
1138 first_pass = true;
1139
1140 double AveragePenalty = 0;
1141
1142 // Vectors to hold exact LLH
1143 std::vector<double> LLH_PredFluc_V(nThrows);
1144 std::vector<double> LLH_DataDraw_V(nThrows);
1145 std::vector<double> LLH_DrawFlucDraw_V(nThrows);
1146
1147 // Loop over the draws
1148 // Should look into multi-threading this. Would require temporary THxx structures from arrays
1149 //KS: Update above would be ideal but currently we loop over samples (see loop below) which isn't as efficient as loop over throws but much much easier to implement
1150 for (unsigned int i = 0; i < nThrows; ++i)
1151 {
1152 if (i % (nThrows/10) == 0) {
1154 }
1155
1156 // Set the total LLH to zero to initialise
1157 double total_llh_data_draw_temp = 0.0;
1158 double total_llh_drawfluc_draw_temp = 0.0;
1159 double total_llh_predfluc_draw_temp = 0.0;
1160
1161 double total_llh_rate_data_draw_temp = 0.0;
1162 double total_llh_rate_predfluc_draw_temp = 0.0;
1163
1164 double total_llh_data_drawfluc_temp = 0.0;
1165 double total_llh_data_predfluc_temp = 0.0;
1166 double total_llh_draw_pred_temp = 0.0;
1167 double total_llh_drawfluc_pred_temp = 0.0;
1168 double total_llh_drawfluc_predfluc_temp = 0.0;
1169 double total_llh_predfluc_pred_temp = 0.0;
1170 double total_llh_datafluc_draw_temp = 0.0;
1171
1172 double total_llh_data_draw_ProjectX_temp = 0.0;
1173 double total_llh_drawfluc_draw_ProjectX_temp = 0.0;
1174
1175 // Save the double that gets written to file
1177 AveragePenalty += llh_penalty;
1178
1179 // Make the Poisson fluctuated hist
1180 std::vector<TH2Poly*> FluctHist(nSamples);
1181 // Also Poisson fluctuate the drawn MCMC hist
1182 std::vector<TH2Poly*> FluctDrawHist(nSamples);
1183 // Finally Poisson fluctuate the data histogram
1184 std::vector<TH2Poly*> DataFlucHist(nSamples);
1185
1186 // Finally Poisson fluctuate the data histogram
1187 std::vector<TH1D*> FluctDrawHistProjectX(nSamples);
1188 std::vector<TH1D*> DrawHistProjectX(nSamples);
1189 std::vector<TH1D*> DrawHistProjectY(nSamples);
1190 std::vector<TH1D*> DrawW2HistProjectX(nSamples);
1191
1192 //KS: We have to clone histograms here to avoid cloning in MP loop, we have to make sure binning matches, content doesn't have to
1193 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1194 {
1195 FluctHist[SampleNum] = static_cast<TH2Poly*>(MeanHist[SampleNum]->Clone());
1196 FluctDrawHist[SampleNum] = static_cast<TH2Poly*>(MeanHist[SampleNum]->Clone());
1197 DataFlucHist[SampleNum] = static_cast<TH2Poly*>(MeanHist[SampleNum]->Clone());
1198
1199 FluctDrawHistProjectX[SampleNum] = static_cast<TH1D*>(DataHist_ProjectX[SampleNum]->Clone());
1200
1201 // Get the ith draw for the jth sample
1202 TH2Poly *DrawHist = MCVector[i][SampleNum];
1203 TH2Poly *DrawW2Hist = W2MCVector[i][SampleNum];
1204
1205 //ProjectPoly calls new TH1D under the hood, never define new ROOT object under MP...
1206 DrawHistProjectX[SampleNum] = ProjectPoly(DrawHist, true, SampleNum);
1207 DrawW2HistProjectX[SampleNum] = ProjectPoly(DrawW2Hist, true, SampleNum);
1208 DrawHistProjectY[SampleNum] = ProjectPoly(DrawHist, false, SampleNum);
1209 }
1210 #ifdef MULTITHREAD
1211 //KS: might be most obscure OMP reduction I have ever made...
1212 #pragma omp parallel for reduction(+:total_llh_data_draw_temp, total_llh_drawfluc_draw_temp, total_llh_predfluc_draw_temp, total_llh_rate_data_draw_temp, total_llh_rate_predfluc_draw_temp, total_llh_data_drawfluc_temp, total_llh_data_predfluc_temp, total_llh_draw_pred_temp, total_llh_drawfluc_pred_temp, total_llh_drawfluc_predfluc_temp, total_llh_predfluc_pred_temp, total_llh_datafluc_draw_temp, total_llh_data_draw_ProjectX_temp, total_llh_drawfluc_draw_ProjectX_temp)
1213 #endif
1214 // Loop over the samples
1215 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1216 {
1217 // Get the ith draw for the jth sample
1218 TH2Poly *DrawHist = MCVector[i][SampleNum];
1219 TH2Poly *DrawW2Hist = W2MCVector[i][SampleNum];
1220 // Skip empty samples
1221 if (DrawHist == nullptr) continue;
1222
1223 // Add LLH penalties from the systematics to the LLH that use the drawn histogram
1224 // Data vs Draw
1225 llh_data_draw[SampleNum] = llh_penalty;
1226 // Fluctuated Draw vs Draw
1227 llh_drawfluc_draw[SampleNum] = llh_penalty;
1228 // Fluctuated Predicitve vs Draw
1229 llh_predfluc_draw[SampleNum] = llh_penalty;
1230
1231 // Data vs Draw using rate
1232 llh_rate_data_draw[SampleNum] = llh_penalty;
1233 // Fluctuated Predicitve vs Draw using rate
1235
1236 // Data vs Fluctuated Draw
1237 llh_data_drawfluc[SampleNum] = llh_penalty;
1238 // Draw vs Predictive
1239 llh_draw_pred[SampleNum] = llh_penalty;
1240 // Fluctuated Draw vs Predictive
1241 llh_drawfluc_pred[SampleNum] = llh_penalty;
1242 // Fluctuated Draw vs Fluctuated Predictive
1244 // Fluctuated Data vs Draw
1245 llh_datafluc_draw[SampleNum] = llh_penalty;
1246
1247 //Some LLH for 1D projections
1250
1251 //Other get 0 penalty term
1252 // Fluctuated Predictive vs Predictive
1253 llh_predfluc_pred[SampleNum] = 0.0;
1254 // Data vs Fluctuated Predictive
1255 llh_data_predfluc[SampleNum] = 0.0;
1256
1257 // Make the Poisson fluctuated hist
1258 MakeFluctuatedHistogram(FluctHist[SampleNum], MeanHist[SampleNum]);
1259 // Also Poisson fluctuate the drawn MCMC hist
1260 MakeFluctuatedHistogram(FluctDrawHist[SampleNum], DrawHist);
1261 // Finally Poisson fluctuate the data histogram
1262 MakeFluctuatedHistogram(DataFlucHist[SampleNum], DataHist[SampleNum]);
1263
1264 // Likelihood between the drawn histogram and the data
1265 const double DataDrawLLH = GetLLH(DataHist[SampleNum], DrawHist, DrawW2Hist);
1266 llh_data_draw[SampleNum] += DataDrawLLH;
1267 total_llh_data_draw_temp += DataDrawLLH;
1268
1269 // Likelihood between drawn histogram and fluctuated drawn histogram
1270 const double DrawFlucDrawLLH = GetLLH(FluctDrawHist[SampleNum], DrawHist, DrawW2Hist);
1271 llh_drawfluc_draw[SampleNum] += DrawFlucDrawLLH;
1272 total_llh_drawfluc_draw_temp += DrawFlucDrawLLH;
1273
1274 // Likelihood between drawn histogram and fluctuated posterior predictive distribution
1275 const double PredFlucDrawLLH = GetLLH(FluctHist[SampleNum], DrawHist, DrawW2Hist);
1276 llh_predfluc_draw[SampleNum] += PredFlucDrawLLH;
1277 total_llh_predfluc_draw_temp += PredFlucDrawLLH;
1278
1279//Rate Based p-value
1280 // Likelihood between the drawn histogram and the data
1281 const double RateDataDrawLLH = SampleHandler->GetTestStatLLH(NoOverflowIntegral(DataHist[SampleNum]), NoOverflowIntegral(DrawHist), NoOverflowIntegral(DrawW2Hist));
1282 llh_rate_data_draw[SampleNum] += RateDataDrawLLH;
1283 total_llh_rate_data_draw_temp += RateDataDrawLLH;
1284
1285 // Likelihood between drawn histogram and fluctuated posterior predictive distribution using rate
1286 const double RatePredFlucDrawLLH = SampleHandler->GetTestStatLLH(NoOverflowIntegral(FluctHist[SampleNum]), NoOverflowIntegral(DrawHist), NoOverflowIntegral(DrawW2Hist));
1287 llh_rate_predfluc_draw[SampleNum] += RatePredFlucDrawLLH;
1288 total_llh_rate_predfluc_draw_temp += RatePredFlucDrawLLH;
1289
1290// All LLH below are for validation reason but not used for final P-Value
1291 // Likelihood between the fluctuated drawn histogram and the data
1292 const double DataDrawFlucLLH = GetLLH(DataHist[SampleNum], FluctDrawHist[SampleNum], DrawW2Hist);
1293 llh_data_drawfluc[SampleNum] += DataDrawFlucLLH;
1294 total_llh_data_drawfluc_temp += DataDrawFlucLLH;
1295
1296 // Likelihood between the drawn histogram and the data
1297 const double DataPredFlucLLH = GetLLH(DataHist[SampleNum], FluctHist[SampleNum], W2MeanHist[SampleNum]);
1298 llh_data_predfluc[SampleNum] += DataPredFlucLLH;
1299 total_llh_data_predfluc_temp += DataPredFlucLLH;
1300
1301 // Likelihood between the drawn hist and the Posterior Predictive
1302 const double DrawPredLLH = GetLLH(DrawHist, MeanHist[SampleNum], W2MeanHist[SampleNum]);
1303 llh_draw_pred[SampleNum] += DrawPredLLH;
1304 total_llh_draw_pred_temp += DrawPredLLH;
1305
1306 // Likelihood between fluctuated drawn and predictive
1307 const double DrawFlucPredLLH = GetLLH(FluctDrawHist[SampleNum], MeanHist[SampleNum], W2MeanHist[SampleNum]);
1308 llh_drawfluc_pred[SampleNum] += DrawFlucPredLLH;
1309 total_llh_drawfluc_pred_temp += DrawFlucPredLLH;
1310
1311 // Likelihood between drawn histogram and fluctuated drawn histogram
1312 const double DrawFlucPredFlucLLH = GetLLH(FluctDrawHist[SampleNum], FluctHist[SampleNum], W2MeanHist[SampleNum]);
1313 llh_drawfluc_predfluc[SampleNum] += DrawFlucPredFlucLLH;
1314 total_llh_drawfluc_predfluc_temp += DrawFlucPredFlucLLH;
1315
1316 // Likelihood between the fluctuated drawn histogram and the posterior predictive
1317 const double PredFlucPredLLH = GetLLH(FluctHist[SampleNum], MeanHist[SampleNum], W2MeanHist[SampleNum]);
1318 llh_predfluc_pred[SampleNum] += PredFlucPredLLH;
1319 total_llh_predfluc_pred_temp += PredFlucPredLLH;
1320
1321 // Likelihood between fluctuated data histogram and drawn histogram
1322 const double DataFlucDrawLLH = GetLLH(DataFlucHist[SampleNum], DrawHist, DrawW2Hist);
1323 llh_datafluc_draw[SampleNum] += DataFlucDrawLLH;
1324 total_llh_datafluc_draw_temp += DataFlucDrawLLH;
1325
1326 lnLHist_Sample_DrawData[SampleNum]->Fill(DataDrawLLH);
1327 lnLHist_Sample_DrawflucDraw[SampleNum]->Fill(DrawFlucDrawLLH);
1328 lnLHist_Sample_PredflucDraw[SampleNum]->Fill(PredFlucDrawLLH);
1329
1330// At the end we leave LLH for 1D projections
1331 MakeFluctuatedHistogram(FluctDrawHistProjectX[SampleNum], DrawHistProjectX[SampleNum]);
1332
1333 // Likelihood between the drawn histogram and the data for muon momentum
1334 const double DataDrawLLH_ProjectX = GetLLH(DataHist_ProjectX[SampleNum], DrawHistProjectX[SampleNum], DrawW2HistProjectX[SampleNum]);
1335 llh_data_draw_ProjectX[SampleNum] += DataDrawLLH_ProjectX;
1336 total_llh_data_draw_ProjectX_temp += DataDrawLLH_ProjectX;
1337
1338 const double DrawFlucDrawLLH_ProjectX = GetLLH(FluctDrawHistProjectX[SampleNum], DrawHistProjectX[SampleNum], DrawW2HistProjectX[SampleNum]);
1339 llh_drawfluc_draw_ProjectX[SampleNum] += DrawFlucDrawLLH_ProjectX;
1340 total_llh_drawfluc_draw_ProjectX_temp += DrawFlucDrawLLH_ProjectX;
1341
1342 //KS: This might seem complicated but we make X and Y projection for each sample. Then we add this to the Violin plot making nice Gaussian in each kineatmical bin of x and y axis
1343 FastViolinFill(ViolinHists_ProjectX[SampleNum], DrawHistProjectX[SampleNum]);
1344 FastViolinFill(ViolinHists_ProjectY[SampleNum], DrawHistProjectY[SampleNum]);
1345 } // End loop over samples (still looping throws)
1346
1347 // Delete the temporary histograms
1348 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1349 {
1350 delete FluctHist[SampleNum];
1351 delete FluctDrawHist[SampleNum];
1352 delete DataFlucHist[SampleNum];
1353 delete FluctDrawHistProjectX[SampleNum];
1354 delete DrawHistProjectX[SampleNum];
1355 delete DrawHistProjectY[SampleNum];
1356 delete DrawW2HistProjectX[SampleNum];
1357 }
1358
1359 total_llh_data_draw = total_llh_data_draw_temp;
1360 total_llh_drawfluc_draw = total_llh_drawfluc_draw_temp;
1361 total_llh_predfluc_draw = total_llh_predfluc_draw_temp;
1362
1363 total_llh_rate_data_draw = total_llh_rate_data_draw_temp;
1364 total_llh_rate_predfluc_draw = total_llh_rate_predfluc_draw_temp;
1365
1366 total_llh_data_drawfluc = total_llh_data_drawfluc_temp;
1367 total_llh_data_predfluc = total_llh_data_predfluc_temp;
1368 total_llh_draw_pred = total_llh_draw_pred_temp;
1369 total_llh_drawfluc_pred = total_llh_drawfluc_pred_temp;
1370 total_llh_drawfluc_predfluc = total_llh_drawfluc_predfluc_temp;
1371 total_llh_predfluc_pred = total_llh_predfluc_pred_temp;
1372 total_llh_datafluc_draw = total_llh_datafluc_draw_temp;
1373
1374 total_llh_data_draw_ProjectX = total_llh_data_draw_ProjectX_temp;
1375 total_llh_drawfluc_draw_ProjectX = total_llh_drawfluc_draw_ProjectX_temp;
1376
1377 // Add LLH penalties from the systematics to the LLH that use the drawn histogram
1381 //Rate based
1384
1389
1392
1397
1400
1402
1404
1405 // Also save to arrays to make sure we have the utmost super accuracy
1406 LLH_PredFluc_V[i] = total_llh_predfluc_draw;
1407 LLH_DataDraw_V[i] = total_llh_data_draw;
1408 LLH_DrawFlucDraw_V[i] = total_llh_drawfluc_draw;
1409
1410 // Write to the output tree
1411 OutputTree->Fill();
1412 } // End loop over throws
1413
1414 AveragePenalty = AveragePenalty/double(nThrows);
1415 MACH3LOG_INFO("Average LLH penalty over toys is {:.2f}", AveragePenalty);
1416 // Calculate exact p-value instead of binned
1417 unsigned int Accept_PredFluc = 0;
1418 unsigned int Accept_DrawFluc = 0;
1419 for (unsigned int i = 0; i < nThrows; ++i)
1420 {
1421 if (LLH_DataDraw_V[i] > LLH_DrawFlucDraw_V[i]) Accept_DrawFluc++;
1422 if (LLH_DataDraw_V[i] > LLH_PredFluc_V[i]) Accept_PredFluc++;
1423 }
1424 const double pvalue_DrawFluc = double(Accept_DrawFluc)/double(nThrows);
1425 const double pvalue_PredFluc = double(Accept_PredFluc)/double(nThrows);
1426
1427 MACH3LOG_INFO("Calculated exact p-value using Fluctuation of Draw: {:.2f}", pvalue_DrawFluc);
1428 MACH3LOG_INFO("Calculated exact p-value using Fluctuation of Prediction: {:.2f}", pvalue_PredFluc);
1429}
1430
1431// *******************
1432// Make the cut LLH histogram
1434// *******************
1435 Outputfile->cd();
1436 MakeCutLLH1D(lnLHist.get());
1440
1445}
1446
1447// ****************
1448// Make the 1D cut distribution and give the 1D p-value
1449void SampleSummary::MakeCutLLH1D(TH1D *Histogram, double llh_ref) {
1450// ****************
1451 const double TotalIntegral = Histogram->Integral();
1452 double Above = 0.0;
1453 // Get the LLH reference from total llh or some reference histogram
1454 double llh_reference = 0.0;
1455 if (llh_ref >= 0) {
1456 llh_reference = llh_ref;
1457 } else {
1458 llh_reference = llh_total;
1459 }
1460 for (int i = 0; i < Histogram->GetXaxis()->GetNbins(); ++i) {
1461 const double xvalue = Histogram->GetBinCenter(i+1);
1462 if (xvalue >= llh_reference) {
1463 Above += Histogram->GetBinContent(i+1);
1464 }
1465 }
1466 const double pvalue = Above/TotalIntegral;
1467 std::stringstream ss;
1468 ss << int(Above) << "/" << int(TotalIntegral) << "=" << pvalue;
1469 Histogram->SetTitle((std::string(Histogram->GetTitle())+"_"+ss.str()).c_str());
1470
1471 // Write a TCanvas and make a line and a filled histogram
1472 auto TempLine = std::make_unique<TLine>(llh_reference , Histogram->GetMinimum(), llh_reference, Histogram->GetMaximum());
1473 TempLine->SetLineColor(kBlack);
1474 TempLine->SetLineWidth(2);
1475
1476 // Make the fill histogram
1477 TH1D *TempHistogram = static_cast<TH1D*>(Histogram->Clone());
1478 TempHistogram->SetFillStyle(1001);
1479 TempHistogram->SetFillColor(kRed);
1480 for (int i = 0; i < TempHistogram->GetNbinsX(); ++i) {
1481 if (TempHistogram->GetBinCenter(i+1) < llh_reference) {
1482 TempHistogram->SetBinContent(i+1, 0.0);
1483 }
1484 }
1485
1486 auto Legend = std::make_unique<TLegend>(0.6, 0.6, 0.9, 0.9);
1487 Legend->SetFillColor(0);
1488 Legend->SetFillStyle(0);
1489 Legend->SetLineWidth(0);
1490 Legend->SetLineColor(0);
1491 Legend->AddEntry(TempLine.get(), Form("Reference LLH, %.0f, p-value=%.2f", llh_reference, pvalue), "l");
1492 Legend->AddEntry(Histogram, Form("LLH, #mu=%.1f#pm%.1f", Histogram->GetMean(), Histogram->GetRMS()), "l");
1493 std::string Title = Histogram->GetName();
1494 Title += "_canv";
1495 TCanvas *TempCanvas = new TCanvas(Title.c_str(), Title.c_str(), 1024, 1024);
1496 TempCanvas->SetGridx();
1497 TempCanvas->SetGridy();
1498 Histogram->Draw();
1499 TempHistogram->Draw("same");
1500 TempLine->Draw("same");
1501 Legend->Draw("same");
1502
1503 TempCanvas->Write();
1504
1505 delete TempHistogram;
1506 delete TempCanvas;
1507}
1508
1509// ****************
1510// Make the 1D Event Rate Hist
1511void SampleSummary::MakeCutEventRate(TH1D *Histogram, const double DataRate) {
1512// ****************
1513 // For the event rate histogram add a TLine to the data rate
1514 auto TempLine = std::make_unique<TLine>(DataRate, Histogram->GetMinimum(), DataRate, Histogram->GetMaximum());
1515 TempLine->SetLineColor(kRed);
1516 TempLine->SetLineWidth(2);
1517 // Also fit a Gaussian because why not?
1518 TF1 *Fitter = new TF1("Fit", "gaus", Histogram->GetBinLowEdge(1), Histogram->GetBinLowEdge(Histogram->GetNbinsX()+1));
1519 Histogram->Fit(Fitter, "RQ");
1520 Fitter->SetLineColor(kRed-5);
1521 // Calculate a p-value
1522 double Above = 0.0;
1523 for (int z = 0; z < Histogram->GetNbinsX(); ++z) {
1524 const double xvalue = Histogram->GetBinCenter(z+1);
1525 if (xvalue >= DataRate) {
1526 Above += Histogram->GetBinContent(z+1);
1527 }
1528 }
1529 const double pvalue = Above/Histogram->Integral();
1530 auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
1531 Legend->SetFillColor(0);
1532 Legend->SetFillStyle(0);
1533 Legend->SetLineWidth(0);
1534 Legend->SetLineColor(0);
1535 Legend->AddEntry(TempLine.get(), Form("Data, %.0f, p-value=%.2f", DataRate, pvalue), "l");
1536 Legend->AddEntry(Histogram, Form("MC, #mu=%.1f#pm%.1f", Histogram->GetMean(), Histogram->GetRMS()), "l");
1537 Legend->AddEntry(Fitter, Form("Gauss, #mu=%.1f#pm%.1f", Fitter->GetParameter(1), Fitter->GetParameter(2)), "l");
1538 std::string TempTitle = std::string(Histogram->GetName());
1539 TempTitle += "_canv";
1540 TCanvas *TempCanvas = new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
1541 TempCanvas->SetGridx();
1542 TempCanvas->SetGridy();
1543 TempCanvas->SetRightMargin(0.03);
1544 TempCanvas->SetBottomMargin(0.08);
1545 TempCanvas->SetLeftMargin(0.10);
1546 TempCanvas->SetTopMargin(0.06);
1547 TempCanvas->cd();
1548 Histogram->Draw();
1549 TempLine->Draw("same");
1550 Fitter->Draw("same");
1551 Legend->Draw("same");
1552 TempCanvas->Write();
1553 Histogram->Write();
1554
1555 delete TempCanvas;
1556 delete Fitter;
1557}
1558
1559// ****************
1560// Calculate the LLH for TH1D, set the LLH to title of MCHist
1561void SampleSummary::CalcLLH(TH1D * const & DatHist, TH1D * const & MCHist, TH1D * const & W2Hist) {
1562// ****************
1563 const double llh = GetLLH(DatHist, MCHist, W2Hist);
1564 std::stringstream ss;
1565 ss << "_2LLH=" << llh;
1566 MCHist->SetTitle((std::string(MCHist->GetTitle())+ss.str()).c_str());
1567 MACH3LOG_INFO("{:<55} {:<10.2f} {:<10.2f} {:<10.2f}", MCHist->GetName(), DatHist->Integral(), MCHist->Integral(), llh);
1568}
1569
1570// ****************
1571// Calculate the LLH for TH1D, set the LLH to title of MCHist
1572void SampleSummary::CalcLLH(TH2Poly * const & DatHist, TH2Poly * const & MCHist, TH2Poly * const & W2Hist) {
1573// ****************
1574 const double llh = GetLLH(DatHist, MCHist, W2Hist);
1575 std::stringstream ss;
1576 ss << "_2LLH=" << llh;
1577 MCHist->SetTitle((std::string(MCHist->GetTitle())+ss.str()).c_str());
1578 MACH3LOG_INFO("{:<55} {:<10.2f} {:<10.2f} {:<10.2f}", MCHist->GetName(), NoOverflowIntegral(DatHist), NoOverflowIntegral(MCHist), llh);
1579}
1580
1581// ****************
1582double SampleSummary::GetLLH(TH2Poly * const & DatHist, TH2Poly * const & MCHist, TH2Poly * const & W2Hist) {
1583// ****************
1584 double llh = 0.0;
1585 for (int i = 1; i < DatHist->GetNumberOfBins()+1; ++i)
1586 {
1587 const double data = DatHist->GetBinContent(i);
1588 const double mc = MCHist->GetBinContent(i);
1589 const double w2 = W2Hist->GetBinContent(i);
1590 llh += SampleHandler->GetTestStatLLH(data, mc, w2);
1591 }
1592 //KS: do times 2 because banff reports chi2
1593 return 2*llh;
1594}
1595
1596// ****************
1597double SampleSummary::GetLLH(TH1D * const & DatHist, TH1D * const & MCHist, TH1D * const & W2Hist) {
1598// ****************
1599 double llh = 0.0;
1600 for (int i = 1; i <= DatHist->GetXaxis()->GetNbins(); ++i)
1601 {
1602 const double data = DatHist->GetBinContent(i);
1603 const double mc = MCHist->GetBinContent(i);
1604 const double w2 = W2Hist->GetBinContent(i);
1605 llh += SampleHandler->GetTestStatLLH(data, mc, w2);
1606 }
1607 //KS: do times 2 because banff reports chi2
1608 return 2*llh;
1609}
1610
1611// ****************
1613// ****************
1614 // Make a new directory
1615 TDirectory *BetaDir = Outputfile->mkdir("BetaParameters");
1616 BetaDir->cd();
1617
1618 int originalErrorLevel = gErrorIgnoreLevel;
1619
1620 //To avoid Warning in <Fit>: Fit data is empty
1621 gErrorIgnoreLevel = kFatal;
1622
1623 MACH3LOG_INFO("Writing Beta parameters");
1624 std::vector<TDirectory *> DirBeta(nSamples);
1625 for (int i = 0; i < nSamples; ++i)
1626 {
1627 // Make a new directory
1628 DirBeta[i] = BetaDir->mkdir((SampleNames[i]).c_str());
1629 DirBeta[i]->cd();
1630
1631 // Loop over each pmu cosmu bin
1632 for (int j = 1; j < maxBins[i]+1; ++j)
1633 {
1634 const double data = DataHist[i]->GetBinContent(j);
1635 const double mc = NominalHist[i]->GetBinContent(j);
1636 const double w2 = W2NomHist[i]->GetBinContent(j);
1637
1638 const double BetaPrior = GetBetaParameter(data, mc, w2, likelihood);
1639
1640 auto TempLine = std::unique_ptr<TLine>(new TLine(BetaPrior, BetaHist[i][j]->GetMinimum(), BetaPrior, BetaHist[i][j]->GetMaximum()));
1641 TempLine->SetLineColor(kRed);
1642 TempLine->SetLineWidth(2);
1643
1644 // Also fit a Gaussian because why not?
1645 TF1 *Fitter = new TF1("Fit", "gaus", BetaHist[i][j]->GetBinLowEdge(1), BetaHist[i][j]->GetBinLowEdge(BetaHist[i][j]->GetNbinsX()+1));
1646 BetaHist[i][j]->Fit(Fitter, "RQ");
1647 Fitter->SetLineColor(kRed-5);
1648
1649 auto Legend = std::make_unique<TLegend>(0.4, 0.75, 0.98, 0.90);
1650 Legend->SetFillColor(0);
1651 Legend->SetFillStyle(0);
1652 Legend->SetLineWidth(0);
1653 Legend->SetLineColor(0);
1654 Legend->AddEntry(TempLine.get(), Form("Prior #mu=%.4f, N_{data}=%.0f", BetaPrior, data), "l");
1655 Legend->AddEntry(BetaHist[i][j].get(), Form("Post, #mu=%.4f#pm%.4f", BetaHist[i][j]->GetMean(), BetaHist[i][j]->GetRMS()), "l");
1656 Legend->AddEntry(Fitter, Form("Gauss, #mu=%.4f#pm%.4f", Fitter->GetParameter(1), Fitter->GetParameter(2)), "l");
1657 std::string TempTitle = std::string(BetaHist[i][j]->GetName());
1658
1659 TempTitle += "_canv";
1660 TCanvas *TempCanvas = new TCanvas(TempTitle.c_str(), TempTitle.c_str(), 1024, 1024);
1661 TempCanvas->SetGridx();
1662 TempCanvas->SetGridy();
1663 TempCanvas->SetRightMargin(0.03);
1664 TempCanvas->SetBottomMargin(0.08);
1665 TempCanvas->SetLeftMargin(0.10);
1666 TempCanvas->SetTopMargin(0.06);
1667 TempCanvas->cd();
1668 BetaHist[i][j]->Draw();
1669 TempLine->Draw("same");
1670 Fitter->Draw("same");
1671 Legend->Draw("same");
1672 TempCanvas->Write();
1673 BetaHist[i][j]->Write();
1674
1675 delete TempCanvas;
1676 delete Fitter;
1677 }
1678 DirBeta[i]->Write();
1679 delete DirBeta[i];
1680 }
1681 BetaDir->Write();
1682 delete BetaDir;
1683
1684 gErrorIgnoreLevel = originalErrorLevel;
1685 Outputfile->cd();
1686}
1687
1688// ****************
1689// Make a projection
1691// ****************
1692 MACH3LOG_INFO("Calculating Correlations");
1693 TStopwatch timer;
1694 timer.Start();
1695
1696 // Data vs Draw for 1D projection
1697 std::vector<double> NEvents_Sample(nSamples);
1698 double event_rate = 0.;
1699
1700 // The output tree we're going to write to
1701 TTree* Event_Rate_Tree = new TTree("Event_Rate_draws", "Event_Rate_draws");
1702 Event_Rate_Tree->Branch("Event_Rate", &event_rate);
1703 // Loop over the samples and set the addresses of the variables to write to file
1704 for (int i = 0; i < nSamples; ++i)
1705 {
1706 // Get the name
1707 std::string SampleName = SampleNames[i];
1708 //CW: Also strip out - signs because it messes up TBranches
1709 while (SampleName.find("-") != std::string::npos) {
1710 SampleName.replace(SampleName.find("-"), 1, std::string("_"));
1711 }
1712 Event_Rate_Tree->Branch((SampleName+"_Event_Rate").c_str(), &NEvents_Sample[i]);
1713 }
1714
1715 // Holds the total event rate
1716 auto EventHist = std::make_unique<TH1D>("EventHist", "Total Event Rate", 100, 1, -1);
1717 EventHist->SetDirectory(nullptr);
1718 EventHist->GetXaxis()->SetTitle("Total event rate");
1719 EventHist->GetYaxis()->SetTitle("Counts");
1720 EventHist->SetLineWidth(2);
1721
1722 // Holds the event rate for the distribution
1723 std::vector<std::unique_ptr<TH1D>> SumHist(nSamples);
1724 for (int i = 0; i < nSamples; ++i)
1725 {
1726 std::string name = std::string(NominalHist[i]->GetName());
1727 name = name.substr(0, name.find("_nom"));
1728
1729 SumHist[i] = std::make_unique<TH1D>((name+"_sum").c_str(),(name+"_sum").c_str(), 100, 1, -1);
1730 SumHist[i]->GetXaxis()->SetTitle("N_{events}");
1731 SumHist[i]->GetYaxis()->SetTitle("Counts");
1732 double Integral = NoOverflowIntegral(DataHist[i]);
1733 std::stringstream ss;
1734 ss << Integral;
1735 SumHist[i]->SetTitle((std::string(SumHist[i]->GetTitle())+"_"+ss.str()).c_str());
1736 }
1737
1738 for (unsigned int it = 0; it < nThrows; ++it)
1739 {
1740 double event_rate_temp = 0.;
1741 // Loop over the samples
1742 #ifdef MULTITHREAD
1743 #pragma omp parallel for reduction(+:event_rate_temp)
1744 #endif
1745 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1746 {
1747 NEvents_Sample[SampleNum] = NoOverflowIntegral(MCVector[it][SampleNum]);
1748 // Fill the sum histogram with the integral of the sampled distribution
1749 SumHist[SampleNum]->Fill(NEvents_Sample[SampleNum], WeightVector[it]);
1750
1751 event_rate_temp += NEvents_Sample[SampleNum];
1752 } // end samples loop
1753 event_rate = event_rate_temp;
1754 EventHist->Fill(event_rate);
1755 Event_Rate_Tree->Fill();
1756 } //end loops over throws
1757 Event_Rate_Tree->Write();
1758 delete Event_Rate_Tree;
1759
1760 double DataRate = 0.0;
1761 #ifdef MULTITHREAD
1762 #pragma omp parallel for reduction(+:DataRate)
1763 #endif
1764 for (int i = 0; i < nSamples; ++i)
1765 {
1766 DataRate += NoOverflowIntegral(DataHist[i]);
1767 }
1768 MakeCutEventRate(EventHist.get(), DataRate);
1769
1770 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1771 {
1772 Dir[SampleNum]->cd();
1773 //Make fancy event rate histogram
1774 MakeCutEventRate(SumHist[SampleNum].get(), NoOverflowIntegral(DataHist[SampleNum]));
1775 }
1776
1777 // Make a new directory
1778 TDirectory *CorrDir = Outputfile->mkdir("Correlations");
1779 CorrDir->cd();
1780
1781 TMatrixDSym* SampleCorrelation = new TMatrixDSym(nSamples);
1782 std::vector<std::vector<std::unique_ptr<TH2D>>> SamCorr(nSamples);
1783 for (int i = 0; i < nSamples; ++i)
1784 {
1785 SamCorr[i].resize(nSamples);
1786
1787 (*SampleCorrelation)(i,i) = 1.0;
1788 const double Min_i = SumHist[i]->GetXaxis()->GetBinLowEdge(1);
1789 const double Max_i = SumHist[i]->GetXaxis()->GetBinUpEdge(SumHist[i]->GetNbinsX()+1);
1790 for (int j = 0; j < nSamples; ++j)
1791 {
1792 const double Min_j = SumHist[j]->GetXaxis()->GetBinLowEdge(1);
1793 const double Max_j = SumHist[j]->GetXaxis()->GetBinUpEdge(SumHist[j]->GetNbinsX()+1);
1794
1795 // TH2D to hold the Correlation
1796 SamCorr[i][j] = std::make_unique<TH2D>(Form("SamCorr_%i_%i", i, j), Form("SamCorr_%i_%i", i, j), 70, Min_i, Max_i, 70, Min_j, Max_j);
1797 SamCorr[i][j]->SetDirectory(nullptr);
1798 SamCorr[i][j]->SetMinimum(0);
1799 SamCorr[i][j]->GetXaxis()->SetTitle(SampleNames[i].c_str());
1800 SamCorr[i][j]->GetYaxis()->SetTitle(SampleNames[j].c_str());
1801 SamCorr[i][j]->GetZaxis()->SetTitle("Events");
1802 }
1803 }
1804
1805 // Now we are sure we have the diagonal elements, let's make the off-diagonals
1806 #ifdef MULTITHREAD
1807 #pragma omp parallel for
1808 #endif
1809 for (int i = 0; i < nSamples; ++i)
1810 {
1811 for (int j = 0; j <= i; ++j)
1812 {
1813 // Skip the diagonal elements which we've already done above
1814 if (j == i) continue;
1815
1816 for (unsigned int it = 0; it < nThrows; ++it)
1817 {
1818 SamCorr[i][j]->Fill(NoOverflowIntegral(MCVector[it][i]), NoOverflowIntegral(MCVector[it][j]));
1819 }
1820 SamCorr[i][j]->Smooth();
1821
1822 // Get the Covariance for these two parameters
1823 (*SampleCorrelation)(i,j) = SamCorr[i][j]->GetCorrelationFactor();
1824 (*SampleCorrelation)(j,i) = (*SampleCorrelation)(i,j);
1825 }// End j loop
1826 }// End i loop
1827
1828 auto hSamCorr = std::make_unique<TH2D>("Sample Correlation", "Sample Correlation", nSamples, 0, nSamples, nSamples, 0, nSamples);
1829 hSamCorr->SetDirectory(nullptr);
1830 hSamCorr->GetZaxis()->SetTitle("Correlation");
1831 hSamCorr->SetMinimum(-1);
1832 hSamCorr->SetMaximum(1);
1833 hSamCorr->GetXaxis()->SetLabelSize(0.015);
1834 hSamCorr->GetYaxis()->SetLabelSize(0.015);
1835
1836 // Loop over the Covariance matrix entries
1837 for (int i = 0; i < nSamples; ++i)
1838 {
1839 hSamCorr->GetXaxis()->SetBinLabel(i+1, SampleNames[i].c_str());
1840
1841 for (int j = 0; j < nSamples; ++j)
1842 {
1843 hSamCorr->GetYaxis()->SetBinLabel(j+1, SampleNames[j].c_str());
1844 // The value of the Covariance
1845 const double corr = (*SampleCorrelation)(i,j);
1846 hSamCorr->SetBinContent(i+1, j+1, corr);
1847 }
1848 }
1849 hSamCorr->Draw("colz");
1850 hSamCorr->Write("Sample_Corr");
1851
1852 SampleCorrelation->Write("Sample_Correlation");
1853 delete SampleCorrelation;
1854
1855 for (int i = 0; i < nSamples; ++i)
1856 {
1857 for (int j = 0; j <= i; ++j)
1858 {
1859 // Skip the diagonal elements which we've already done above
1860 if (j == i) continue;
1861 SamCorr[i][j]->Write();
1862 }// End j loop
1863 }// End i loop
1864
1865 //KS: This can take ages so better turn it off by default
1866 bool DoPerKinemBin = false;
1867 if(DoPerKinemBin)
1868 {
1869 //KS: Now the same but for kinematic bin of each sample
1870 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
1871 {
1872 TMatrixDSym* KinCorrelation = new TMatrixDSym(maxBins[SampleNum]);
1873 std::vector<std::vector<std::unique_ptr<TH2D>>> KinCorr(maxBins[SampleNum]);
1874 for (int i = 0; i < maxBins[SampleNum]; ++i)
1875 {
1876 KinCorr[i].resize(maxBins[SampleNum]);
1877 (*KinCorrelation)(i,i) = 1.0;
1878
1879 const double Min_i = PosteriorHist[SampleNum][i+1]->GetXaxis()->GetBinLowEdge(1);
1880 const double Max_i = PosteriorHist[SampleNum][i+1]->GetXaxis()->GetBinUpEdge(PosteriorHist[SampleNum][i+1]->GetNbinsX()+1);
1881
1882 //Get PolyBin
1883 TH2PolyBin* bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(i));
1884 // Just make a little fancy name
1885 std::stringstream ss2;
1886 ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
1887 ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
1888
1889 for (int j = 0; j < maxBins[SampleNum]; ++j)
1890 {
1891 const double Min_j = PosteriorHist[SampleNum][j+1]->GetXaxis()->GetBinLowEdge(1);
1892 const double Max_j = PosteriorHist[SampleNum][j+1]->GetXaxis()->GetBinUpEdge(PosteriorHist[SampleNum][j+1]->GetNbinsX()+1);
1893
1894 // TH2D to hold the Correlation
1895 KinCorr[i][j] = std::make_unique<TH2D>( Form("Kin_%i_%i_%i", SampleNum, i, j),
1896 Form("Kin_%i_%i_%i", SampleNum, i, j), 70, Min_i, Max_i, 70, Min_j, Max_j);
1897 KinCorr[i][j]->SetDirectory(nullptr);
1898 KinCorr[i][j]->SetMinimum(0);
1899
1900 KinCorr[i][j]->GetXaxis()->SetTitle(ss2.str().c_str());
1901
1902 bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(j));
1903 // Just make a little fancy name
1904 std::stringstream ss3;
1905 ss3 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
1906 ss3 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
1907 KinCorr[i][j]->GetYaxis()->SetTitle(ss3.str().c_str());
1908 KinCorr[i][j]->GetZaxis()->SetTitle("Events");
1909 }
1910 }
1911 // Now we are sure we have the diagonal elements, let's make the off-diagonals
1912 #ifdef MULTITHREAD
1913 #pragma omp parallel for
1914 #endif
1915 for (int i = 0; i < maxBins[SampleNum]; ++i)
1916 {
1917 for (int j = 0; j <= i; ++j)
1918 {
1919 // Skip the diagonal elements which we've already done above
1920 if (j == i) continue;
1921
1922 for (unsigned int it = 0; it < nThrows; ++it)
1923 {
1924 KinCorr[i][j]->Fill(MCVector[it][SampleNum]->GetBinContent(i+1), MCVector[it][SampleNum]->GetBinContent(j+1));
1925 }
1926 KinCorr[i][j]->Smooth();
1927
1928 // Get the Covariance for these two parameters
1929 (*KinCorrelation)(i,j) = KinCorr[i][j]->GetCorrelationFactor();
1930 (*KinCorrelation)(j,i) = (*KinCorrelation)(i,j);
1931 }// End j loop
1932 }// End i loop
1933
1934 auto hKinCorr = std::make_unique<TH2D>(SampleNames[SampleNum].c_str(), SampleNames[SampleNum].c_str(),
1935 maxBins[SampleNum], 0, maxBins[SampleNum], maxBins[SampleNum], 0, maxBins[SampleNum]);
1936 hKinCorr->SetDirectory(nullptr);
1937 hKinCorr->GetZaxis()->SetTitle("Correlation");
1938 hKinCorr->SetMinimum(-1);
1939 hKinCorr->SetMaximum(1);
1940 hKinCorr->GetXaxis()->SetLabelSize(0.015);
1941 hKinCorr->GetYaxis()->SetLabelSize(0.015);
1942
1943 // Loop over the Covariance matrix entries
1944 for (int i = 0; i < maxBins[SampleNum]; ++i)
1945 {
1946 //Get PolyBin
1947 TH2PolyBin* bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(i));
1948 // Just make a little fancy name
1949 std::stringstream ss2;
1950 ss2 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
1951 ss2 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
1952 hKinCorr->GetXaxis()->SetBinLabel(i+1, ss2.str().c_str());
1953
1954 for (int j = 0; j < maxBins[SampleNum]; ++j)
1955 {
1956 bin = static_cast<TH2PolyBin*>(NominalHist[SampleNum]->GetBins()->At(j));
1957 // Just make a little fancy name
1958 std::stringstream ss3;
1959 ss3 << "p_{#mu} (" << bin->GetXMin() << "-" << bin->GetXMax() << ")";
1960 ss3 << " cos#theta_{#mu} (" << bin->GetYMin() << "-" << bin->GetYMax() << ")";
1961 KinCorr[i][j]->GetYaxis()->SetTitle(ss3.str().c_str());
1962
1963 hKinCorr->GetYaxis()->SetBinLabel(j+1, ss3.str().c_str());
1964 // The value of the Covariance
1965 const double corr = (*KinCorrelation)(i,j);
1966 hKinCorr->SetBinContent(i+1, j+1, corr);
1967 }
1968 }
1969 hKinCorr->Draw("colz");
1970 hKinCorr->Write((SampleNames[SampleNum] + "_Corr").c_str());
1971
1972 KinCorrelation->Write((SampleNames[SampleNum] + "_Correlation").c_str());
1973 delete KinCorrelation;
1974
1975 /*
1976 for (int i = 0; i < maxBins[SampleNum]; ++i)
1977 {
1978 for (int j = 0; j <= i; ++j)
1979 {
1980 // Skip the diagonal elements which we've already done above
1981 if (j == i) continue;
1982 KinCorr[i][j]->Write();
1983 }// End j loop
1984 }// End i loop
1985 */
1986 }//end loop over samples
1987 }//end if DoPerKinemBin
1988 else
1989 {
1990 MACH3LOG_INFO("Not calculating correlations per each kinematic bin");
1991 }
1992
1993 if(DoByModePlots)
1994 {
1995 // Holds the total event rate by mode
1996 std::vector<TH1D*> EventHist_ByMode(Modes->GetNModes()+1);
1997 for (int j = 0; j < Modes->GetNModes()+1; j++)
1998 {
1999 std::string ModeName = Modes->GetMaCh3ModeName(j);
2000 EventHist_ByMode[j] = new TH1D(Form("EventHist_%s", ModeName.c_str()), Form("Total Event Rate %s", ModeName.c_str()), 100, 1, -1);
2001 EventHist_ByMode[j]->GetXaxis()->SetTitle("Total event rate");
2002 EventHist_ByMode[j]->GetYaxis()->SetTitle("Counts");
2003 EventHist_ByMode[j]->SetLineWidth(2);
2004 }
2005
2006 //KS: Here we calculate total event rates for each mode, maybe not most efficient but can be improved in the future
2007 for (unsigned int it = 0; it < nThrows; ++it)
2008 {
2009 for (int j = 0; j < Modes->GetNModes()+1; j++)
2010 {
2011 double event_rate_temp = 0.;
2012 #ifdef MULTITHREAD
2013 #pragma omp parallel for reduction(+:event_rate_temp)
2014 #endif
2015 for (int SampleNum = 0; SampleNum < nSamples; SampleNum++)
2016 {
2017 event_rate_temp += NoOverflowIntegral(MCVectorByMode[it][SampleNum][j]);
2018 }
2019 EventHist_ByMode[j]->Fill(event_rate_temp);
2020 }
2021 }
2022
2023 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2024 {
2025 MakeCutEventRate(EventHist_ByMode[i], DataRate);
2026 }
2027
2028 TMatrixDSym* ModeCorrelation = new TMatrixDSym(Modes->GetNModes()+1);
2029
2030 TH2D*** ModeCorr = new TH2D**[Modes->GetNModes()+1]();
2031 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2032 {
2033 ModeCorr[i] = new TH2D*[Modes->GetNModes()+1]();
2034
2035 (*ModeCorrelation)(i,i) = 1.0;
2036
2037 const double Min_i = EventHist_ByMode[i]->GetXaxis()->GetBinLowEdge(1);
2038 const double Max_i = EventHist_ByMode[i]->GetXaxis()->GetBinUpEdge(EventHist_ByMode[i]->GetNbinsX()+1);
2039 for (int j = 0; j < Modes->GetNModes()+1; ++j)
2040 {
2041 const double Min_j = EventHist_ByMode[j]->GetXaxis()->GetBinLowEdge(1);
2042 const double Max_j = EventHist_ByMode[j]->GetXaxis()->GetBinUpEdge(EventHist_ByMode[j]->GetNbinsX()+1);
2043
2044 // TH2D to hold the Correlation
2045 ModeCorr[i][j] = new TH2D(Form("ModeCorr_%i_%i",i,j), Form("ModeCorr_%i_%i",i,j), 70, Min_i, Max_i, 70, Min_j, Max_j);
2046 ModeCorr[i][j]->SetDirectory(nullptr);
2047 ModeCorr[i][j]->SetMinimum(0);
2048 ModeCorr[i][j]->GetXaxis()->SetTitle(Modes->GetMaCh3ModeName(i).c_str());
2049 ModeCorr[i][j]->GetYaxis()->SetTitle(Modes->GetMaCh3ModeName(j).c_str());
2050 ModeCorr[i][j]->GetZaxis()->SetTitle("Events");
2051 }
2052 }
2053
2054 // Now we are sure we have the diagonal elements, let's make the off-diagonals
2055 #ifdef MULTITHREAD
2056 #pragma omp parallel for
2057 #endif
2058 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2059 {
2060 for (int j = 0; j <= i; ++j)
2061 {
2062 // Skip the diagonal elements which we've already done above
2063 if (j == i) continue;
2064
2065 for (unsigned int it = 0; it < nThrows; ++it)
2066 {
2067 double Integral_X = 0.;
2068 double Integral_Y = 0.;
2069 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
2070 {
2071 Integral_X += NoOverflowIntegral(MCVectorByMode[it][SampleNum][i]);
2072 Integral_Y += NoOverflowIntegral(MCVectorByMode[it][SampleNum][j]);
2073 }
2074 ModeCorr[i][j]->Fill(Integral_X, Integral_Y);
2075 }
2076 ModeCorr[i][j]->Smooth();
2077
2078 // Get the Covariance for these two parameters
2079 (*ModeCorrelation)(i,j) = ModeCorr[i][j]->GetCorrelationFactor();
2080 (*ModeCorrelation)(j,i) = (*ModeCorrelation)(i,j);
2081 }// End j loop
2082 }// End i loop
2083
2084 TH2D* hModeCorr = new TH2D("Mode Correlation", "Mode Correlation", Modes->GetNModes()+1, 0, Modes->GetNModes()+1, Modes->GetNModes()+1, 0, Modes->GetNModes()+1);
2085 hModeCorr->SetDirectory(nullptr);
2086 hModeCorr->GetZaxis()->SetTitle("Correlation");
2087 hModeCorr->SetMinimum(-1);
2088 hModeCorr->SetMaximum(1);
2089 hModeCorr->GetXaxis()->SetLabelSize(0.015);
2090 hModeCorr->GetYaxis()->SetLabelSize(0.015);
2091
2092 // Loop over the Covariance matrix entries
2093 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2094 {
2095 hModeCorr->GetXaxis()->SetBinLabel(i+1, Modes->GetMaCh3ModeName(i).c_str());
2096
2097 for (int j = 0; j < Modes->GetNModes()+1; ++j)
2098 {
2099 hModeCorr->GetYaxis()->SetBinLabel(j+1, Modes->GetMaCh3ModeName(j).c_str());
2100 // The value of the Covariance
2101 const double corr = (*ModeCorrelation)(i,j);
2102 hModeCorr->SetBinContent(i+1, j+1, corr);
2103 }
2104 }
2105 hModeCorr->Draw("colz");
2106 hModeCorr->Write("Mode_Corr");
2107 delete hModeCorr;
2108
2109 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2110 {
2111 for (int j = 0; j <= i; ++j)
2112 {
2113 // Skip the diagonal elements which we've already done above
2114 if (j == i) continue;
2115 ModeCorr[i][j]->Write();
2116 }// End j loop
2117 }// End i loop
2118
2119 for (int i = 0; i < Modes->GetNModes()+1; ++i)
2120 {
2121 for (int j = 0; j < Modes->GetNModes()+1; ++j)
2122 {
2123 delete ModeCorr[i][j];
2124 }
2125 delete[] ModeCorr[i];
2126 }
2127 delete[] ModeCorr;
2128 ModeCorrelation->Write("Mode_Correlation");
2129 delete ModeCorrelation;
2130
2131 for (int j = 0; j < Modes->GetNModes()+1; j++)
2132 {
2133 delete EventHist_ByMode[j];
2134 }
2135 }
2136
2137 CorrDir->Close();
2138 delete CorrDir;
2139
2140 timer.Stop();
2141 MACH3LOG_INFO("Calculating correlations took {:.2f}s", timer.RealTime());
2142 Outputfile->cd();
2143}
2144
2145// ****************
2146// Make a projection
2147TH1D* SampleSummary::ProjectHist(TH2D* Histogram, const bool ProjectX) {
2148// ****************
2149 TH1D* Projection = nullptr;
2150 std::string name;
2151 if (ProjectX) {
2152 name = std::string(Histogram->GetName()) + "_x";
2153 Projection = Histogram->ProjectionX(name.c_str(), 1, Histogram->GetYaxis()->GetNbins(), "e");
2154 } else {
2155 name = std::string(Histogram->GetName()) + "_y";
2156 Projection = Histogram->ProjectionY(name.c_str(), 1, Histogram->GetXaxis()->GetNbins(), "e");
2157 }
2158 return Projection;
2159}
2160
2161// ****************
2162// Make a projection
2163TH1D* SampleSummary::ProjectPoly(TH2Poly* Histogram, const bool ProjectX, const int selection, const bool MakeErrorHist) {
2164// ****************
2165 std::vector<double> xbins;
2166 std::vector<double> ybins;
2167
2168 SampleHandler->SetupBinning(M3::int_t(selection), xbins, ybins);
2169 TH1D* Projection = nullptr;
2170 std::string name;
2171 if (ProjectX) {
2172 name = std::string(Histogram->GetName()) + "_x";
2173 Projection = PolyProjectionX(Histogram, name.c_str(), xbins, MakeErrorHist);
2174 } else {
2175 name = std::string(Histogram->GetName()) + "_y";
2176 Projection = PolyProjectionY(Histogram, name.c_str(), ybins, MakeErrorHist);
2177 }
2178 return Projection;
2179}
2180
2181// ****************
2182//KS: We have two methods how to apply statistical fluctuation standard is faster hence is default
2183void SampleSummary::MakeFluctuatedHistogram(TH1D *FluctHist, TH1D* PolyHist){
2184// ****************
2185 if(StandardFluctuation) MakeFluctuatedHistogramStandard(FluctHist, PolyHist, rnd.get());
2186 else MakeFluctuatedHistogramAlternative(FluctHist, PolyHist, rnd.get());
2187}
2188
2189// ****************
2190//KS: We have two methods how to apply statistical fluctuation standard is faster hence is default
2191void SampleSummary::MakeFluctuatedHistogram(TH2Poly *FluctHist, TH2Poly* PolyHist){
2192// ****************
2193 if(StandardFluctuation) MakeFluctuatedHistogramStandard(FluctHist, PolyHist, rnd.get());
2194 else MakeFluctuatedHistogramAlternative(FluctHist, PolyHist, rnd.get());
2195}
2196
2197
2198// ****************
2200// ****************
2201 MACH3LOG_INFO("******************************");
2202 switch(Criterion) {
2203 case M3::kInfCrit::kBIC:
2204 // Study Bayesian Information Criterion
2205 StudyBIC();
2206 break;
2207 case M3::kInfCrit::kDIC:
2208 // Study Deviance Information Criterion
2209 StudyDIC();
2210 break;
2212 // Study Watanabe-Akaike information criterion (WAIC)
2213 StudyWAIC();
2214 break;
2216 MACH3LOG_ERROR("kInfCrits is not a valid kInfCrit!");
2217 throw MaCh3Exception(__FILE__, __LINE__);
2218 default:
2219 MACH3LOG_ERROR("UNKNOWN Information Criterion SPECIFIED!");
2220 MACH3LOG_ERROR("You gave {}", static_cast<int>(Criterion));
2221 throw MaCh3Exception(__FILE__ , __LINE__ );
2222 }
2223 MACH3LOG_INFO("******************************");
2224}
2225
2226// ****************
2228// ****************
2229 //make fancy event rate histogram
2230 double DataRate = 0.0;
2231 double BinsRate = 0.0;
2232 #ifdef MULTITHREAD
2233 #pragma omp parallel for reduction(+:DataRate, BinsRate)
2234 #endif
2235 for (int i = 0; i < nSamples; ++i)
2236 {
2237 if (DataHist[i] == nullptr) continue;
2238 DataRate += NoOverflowIntegral(DataHist[i]);
2239 BinsRate += maxBins[i];
2240 }
2241
2242 const double EventRateBIC = GetBIC(llh_total, DataRate, nModelParams);
2243 const double BinBasedBIC = GetBIC(llh_total, BinsRate, nModelParams);
2244 MACH3LOG_INFO("Calculated Bayesian Information Criterion using global number of events: {:.2f}", EventRateBIC);
2245 MACH3LOG_INFO("Calculated Bayesian Information Criterion using global number of bins: {:.2f}", BinBasedBIC);
2246 MACH3LOG_INFO("Additional info: nModelParams {} DataRate: {:.2f} BinsRate: {:.2f}", nModelParams, DataRate, BinsRate);
2247}
2248
2249// ****************
2250// Get the Deviance Information Criterion (DIC)
2252// ****************
2253 //The posterior mean of the deviance
2254 double Dbar = 0.;
2255
2256 #ifdef MULTITHREAD
2257 #pragma omp parallel for reduction(+:Dbar)
2258 #endif
2259 for (unsigned int i = 0; i < nThrows; ++i)
2260 {
2261 double LLH_temp = 0.;
2262 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum)
2263 {
2264 // Get -2*log-likelihood
2265 LLH_temp += GetLLH(DataHist[SampleNum], MCVector[i][SampleNum], W2MCVector[i][SampleNum]);
2266 }
2267 Dbar += LLH_temp;
2268 }
2269 Dbar = Dbar / nThrows;
2270
2271 // A point estimate of the deviance
2272 const double Dhat = llh_total;
2273
2274 //Effective number of parameters
2275 const double p_D = std::fabs(Dbar - Dhat);
2276
2277 //Actual test stat
2278 const double DIC_stat = Dhat + 2 * p_D;
2279 MACH3LOG_INFO("Effective number of parameters following DIC formalism is equal to: {:.2f}", p_D);
2280 MACH3LOG_INFO("DIC test statistic = {:.2f}", DIC_stat);
2281}
2282
2283// ****************
2284// Get the Watanabe-Akaike information criterion (WAIC)
2286// ****************
2287 // log pointwise predictive density
2288 double lppd = 0.;
2289 // effective number of parameters
2290 double p_WAIC = 0.;
2291
2292 #ifdef MULTITHREAD
2293 #pragma omp parallel for reduction(+:lppd, p_WAIC)
2294 #endif
2295 for (int SampleNum = 0; SampleNum < nSamples; ++SampleNum) {
2296 int nBins = maxBins[SampleNum];
2297 for (int i = 1; i <= nBins; ++i) {
2298 double mean_llh = 0.;
2299 double sum_exp_llh = 0;
2300 double mean_llh_squared = 0.;
2301
2302 for (unsigned int s = 0; s < nThrows; ++s) {
2303 const double data = DataHist[SampleNum]->GetBinContent(i);
2304 const double mc = MCVector[s][SampleNum]->GetBinContent(i);
2305 const double w2 = W2MCVector[s][SampleNum]->GetBinContent(i);
2306
2307 // Get the -log-likelihood for this sample and bin
2308 double neg_LLH_temp = SampleHandler->GetTestStatLLH(data, mc, w2);
2309
2310 // Negate the negative log-likelihood to get the actual log-likelihood
2311 double LLH_temp = -neg_LLH_temp;
2312
2313 mean_llh += LLH_temp;
2314 mean_llh_squared += LLH_temp * LLH_temp;
2315 sum_exp_llh += std::exp(LLH_temp);
2316 }
2317
2318 // Compute the mean log-likelihood and the squared mean
2319 mean_llh /= nThrows;
2320 mean_llh_squared /= nThrows;
2321 sum_exp_llh /= nThrows;
2322 sum_exp_llh = std::log(sum_exp_llh);
2323
2324 // Log pointwise predictive density based on Eq. 4 in Gelman2014
2325 lppd += sum_exp_llh;
2326
2327 // Compute the effective number of parameters for WAIC
2328 p_WAIC += mean_llh_squared - (mean_llh * mean_llh);
2329 }
2330 }
2331
2332 // Compute WAIC, see Eq. 13 in Gelman2014
2333 double WAIC = -2 * (lppd - p_WAIC);
2334 MACH3LOG_INFO("Effective number of parameters following WAIC formalism is equal to: {:.2f}", p_WAIC);
2335 MACH3LOG_INFO("WAIC = {:.2f}", WAIC);
2336}
int size
TH1D * PolyProjectionX(TObject *poly, std::string TempName, const std::vector< double > &xbins, const bool computeErrors)
WP: Poly Projectors.
void MakeFluctuatedHistogramAlternative(TH1D *FluctHist, TH1D *PolyHist, TRandom3 *rand)
Make Poisson fluctuation of TH1D hist using slow method which is only for cross-check.
void NormaliseTH2Poly(TH2Poly *Histogram)
Helper to Normalise histograms.
void MakeFluctuatedHistogramStandard(TH1D *FluctHist, TH1D *PolyHist, TRandom3 *rand)
Make Poisson fluctuation of TH1D hist using default fast method.
void FastViolinFill(TH2D *violin, TH1D *hist_1d)
KS: Fill Violin histogram with entry from a toy.
TH1D * PolyProjectionY(TObject *poly, std::string TempName, const std::vector< double > &ybins, const bool computeErrors)
WP: Poly Projectors.
double NoOverflowIntegral(TH2Poly *poly)
WP: Helper function for calculating binned Integral of TH2Poly i.e not including overflow.
TH2Poly * RatioPolys(TH2Poly *NumHist, TH2Poly *DenomHist)
Helper to make ratio of TH2Polys.
TH2Poly * NormalisePoly(TH2Poly *Histogram)
WP: Helper to Normalise histograms.
#define MACH3LOG_DEBUG
Definition: MaCh3Logger.h:22
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
double GetBIC(const double llh, const int data, const int nPars)
Get the Bayesian Information Criterion (BIC) or Schwarz information criterion (also SIC,...
void Get2DBayesianpValue(TH2D *Histogram)
Calculates the 2D Bayesian p-value and generates a visualization.
double GetModeError(TH1D *hpost)
Get the mode error from a TH1D.
double GetBetaParameter(const double data, const double mc, const double w2, TestStatistic TestStat)
KS: Calculate Beta parameter which will be different based on specified test statistic.
Custom exception class for MaCh3 errors.
int GetNModes() const
KS: Get number of modes, keep in mind actual number is +1 greater due to unknown category.
Definition: MaCh3Modes.h:52
std::string GetMaCh3ModeName(const int Index)
KS: Get normal name of mode, if mode not known you will get UNKNOWN_BAD.
Definition: MaCh3Modes.cpp:138
Class responsible for handling implementation of samples used in analysis, reweighting and returning ...
virtual void SetupBinning(const M3::int_t Selection, std::vector< double > &BinningX, std::vector< double > &BinningY)
virtual std::string GetKinVarLabel(const int sample, const int Dimension)
void AddThrowByMode(std::vector< std::vector< TH2Poly * > > &SampleVector_ByMode)
KS: Add histograms for each mode.
void StudyInformationCriterion(M3::kInfCrit Criterion)
Information Criterion.
std::vector< TH2Poly * > lnLHist_Mode
The LLH distribution in pmu cosmu for using the mode in each bin.
std::vector< double > LLHPenaltyVector
Vector to hold the penalty term.
std::vector< double > llh_data_draw_ProjectX
Projection X (most likely muon momentum) of LLH.
void MakeCutEventRate(TH1D *Histogram, const double DataRate)
Make the 1D Event Rate Hist.
SampleHandlerBase * SampleHandler
Pointer to SampleHandler object, mostly used to get sample names, binning etc.
~SampleSummary()
Destructor.
std::vector< double > llh_predfluc_draw
Fluctuated Predictive vs Draw.
void AddThrow(std::vector< TH2Poly * > &MCHist, std::vector< TH2Poly * > &W2Hist, const double LLHPenalty=0.0, const double Weight=1.0, const int DrawNumber=0)
KS: Add histograms with throws.
std::vector< std::vector< TH2Poly * > > MCVector
Vector of vectors which holds the loaded MC histograms.
std::vector< TH2Poly * > lnLHist_Mean
The LLH distribution in pmu cosmu for using the mean in each bin.
void CalcLLH(TH2Poly *const &Data, TH2Poly *const &MC, TH2Poly *const &W2)
Helper functions to calculate likelihoods using TH2Poly, will modify MC hist tittle to include LLH.
double total_llh_predfluc_pred
Fluctuated Predictive vs Predictive.
double total_llh_drawfluc_draw_ProjectX
Fluctuated Draw vs Draw for projection X (most likely muon momentum)
std::vector< std::vector< TH2Poly * > > W2MCVector
Vector of vectors which holds the loaded W2 histograms.
std::vector< TH2Poly * > W2MeanHist
Pointer to the w2 histograms (for mean values).
double total_llh_data_draw
Data vs Draw.
std::unique_ptr< TH1D > RandomHist
Holds the history of which entries have been drawn in the MCMC file.
double total_llh_rate_predfluc_draw
Fluctuated Predictive vs Draw using Rate.
void MakeChi2Hists()
Make the fluctuated histograms (2D and 1D) for the chi2s Essentially taking the MCMC draws and calcul...
std::vector< double > llh_data_predfluc
Data vs Fluctuated Predictive.
void StudyBIC()
Study Bayesian Information Criterion (BIC) .
std::unique_ptr< TRandom3 > rnd
Random number generator.
std::vector< TH1D * > lnLHist_Mode1D
Holds the bin-by-bin LLH for the mode posterior predictive vs the data.
std::string OutputName
Output filename.
void MakePredictive()
Finalise the distributions from the thrown samples.
std::vector< double > llh_predfluc_pred
Fluctuated Predictive vs Predictive.
std::vector< std::vector< std::unique_ptr< TH1D > > > BetaHist
Distribution of beta parameters in Barlow Beeston formalisms.
int nModelParams
Number of parameters.
std::vector< TH2D * > ViolinHists_ProjectX
Posterior predictive but for X projection but as a violin plot.
bool first_pass
KS: Hacky flag to let us know if this is first toy.
std::vector< TH2Poly * > NominalHist
The nominal histogram for the selection.
std::vector< std::vector< std::unique_ptr< TH1D > > > PosteriorHist
The posterior predictive for the whole selection: this gets built after adding in the toys....
double total_llh_draw_pred
Draw vs Predictive.
std::vector< double > llh_drawfluc_draw
Fluctuated Draw vs Draw.
double total_llh_data_draw_ProjectX
Data vs Draw for projection X (most likely muon momentum)
std::vector< TH2Poly * > MeanHistCorrected
The posterior predictive distribution in pmu cosmu using the mean after applying Barlow-Beeston Corre...
std::vector< double > llh_drawfluc_predfluc
Fluctuated Draw vs Fluctuated Predictive.
MaCh3Modes * Modes
MaCh3 Modes.
std::vector< double > llh_data_draw
Data vs Draw.
double total_llh_data_drawfluc
Data vs Fluctuated Draw.
void AddData(std::vector< TH2Poly * > &DataHist)
KS: Add data histograms.
TH1D * ProjectPoly(TH2Poly *Histogram, const bool ProjectX, const int selection, const bool MakeErrorHist=false)
Helper to project TH2Poly onto axis.
std::unique_ptr< TH1D > lnLHist
The histogram containing the lnL for each throw.
std::vector< TH1D * > lnLHist_Mean1D
Holds the bin-by-bin LLH for the mean posterior predictive vs the data.
TTree * OutputTree
TTree which we save useful information to.
bool DoByModePlots
By mode variables.
std::vector< TH1D * > lnLHist_Mean_ProjectX
The LLH distribution in pmu using the mean in each bin.
void PrepareOutput()
KS: Prepare output tree and necessary variables.
SampleSummary(const int n_Samples, const std::string &Filename, SampleHandlerBase *const sample, const int nSteps)
Constructor.
double total_llh_drawfluc_draw
Fluctuated Draw vs Draw.
std::vector< TH2Poly * > DataHist
The data histogram for the selection.
std::vector< double > llh_draw_pred
Draw vs Predictive.
double total_llh_predfluc_draw
Fluctuated Predictive vs Draw.
int Debug
Tells Debug level to save additional histograms.
std::vector< double > llh_drawfluc_draw_ProjectX
bool CheckSamples(const int Length)
Check the length of samples agrees.
std::unique_ptr< TH2D > lnLDrawHist
The 2D lnLhist, showing (draw vs data) and (draw vs fluct), anything above y=x axis is the p-value.
void AddNominal(std::vector< TH2Poly * > &NominalHist, std::vector< TH2Poly * > &W2Nom)
KS: Add prior histograms.
void StudyDIC()
KS: Get the Deviance Information Criterion (DIC) .
std::vector< TH2D * > ViolinHists_ProjectY
Posterior predictive but for Y projection but as a violin plot.
std::vector< double > llh_data_drawfluc
Data vs Fluctuated Draw.
std::vector< TH2Poly * > ModeHist
The posterior predictive distribution in pmu cosmu using the mode.
std::vector< std::string > SampleNames
name for each sample
std::vector< TH2Poly * > W2NomHist
Pointer to the w2 histograms (for nominal values).
double total_llh_drawfluc_pred
Fluctuated Draw vs Predictive.
void StudyWAIC()
KS: Get the Watanabe-Akaike information criterion (WAIC) .
std::vector< TDirectory * > Dir
Directory for each sample.
std::vector< double > llh_rate_predfluc_draw
Fluctuated Predictive vs Draw using rate only.
unsigned int nThrows
Number of throws.
std::vector< TH1D * > lnLHist_Sample_DrawData
The histogram containing the lnL (draw vs data) for each throw for each sample.
std::vector< std::vector< std::unique_ptr< TH1D > > > w2Hist
The posterior predictive for the whole selection: this gets built after adding in the toys....
TH1D **** PosteriorHist_ByMode
Histogram which corresponds to each bin in the sample's th2poly.
std::unique_ptr< TH1D > lnLHist_drawflucdraw
The lnLhist for the draw vs draw fluctuated.
TFile * Outputfile
Output filename.
double llh_total
Total LLH for the posterior predictive distribution.
void Write()
KS: Write results into root file.
std::unique_ptr< TH2D > lnLDrawHistRate
The 2D lnLhist, showing (draw vs data) and (draw vs fluct), using rate, anything above y=x axis is th...
double total_llh_rate_data_draw
Rate Data vs Draw.
std::unique_ptr< TH1D > lnLHist_drawdata
The lnLhist for the draw vs data.
TestStatistic likelihood
Type of likelihood for example Poisson, Barlow-Beeston or Ice Cube.
int nSamples
Number of samples.
void MakeCutLLH()
Make the cut LLH histogram.
std::vector< std::vector< std::vector< TH2Poly * > > > MCVectorByMode
Vector of vectors which holds the loaded MC histograms for each mode.
std::vector< double > llh_drawfluc_pred
Fluctuated Draw vs Predictive.
std::vector< std::vector< TH2Poly * > > MeanHist_ByMode
The posterior predictive distribution in pmu cosmu using the mean.
std::vector< double > llh_datafluc_draw
Fluctuated Data vs Draw.
std::vector< TH1D * > DataHist_ProjectX
The data histogram for the selection X projection.
void MakeCutLLH1D(TH1D *Histogram, double llh_ref=-999)
std::vector< TH1D * > DataHist_ProjectY
The data histogram for the selection Y projection.
void PlotBetaParameters()
KS: In Barlow Beeston we have Beta Parameters which scale generated MC.
double GetLLH(TH2Poly *const &Data, TH2Poly *const &MC, TH2Poly *const &W2)
Helper functions to calculate likelihoods using TH2Poly.
std::vector< double > WeightVector
Vector holding weight.
void MakeFluctuatedHistogram(TH1D *FluctHist, TH1D *PolyHist)
Make Poisson fluctuation of TH1D hist.
std::vector< double > llh_rate_data_draw
Data vs Draw using rate only.
bool isPriorPredictive
bool whether we have Prior or Posterior Predictive
bool DoBetaParam
Are we making Beta Histograms.
std::vector< TH2Poly * > W2ModeHist
Pointer to the w2 histograms (for mode values).
std::unique_ptr< TH2D > lnLFlucHist
The 2D lnLHist, showing (draw vs data) and (draw vs draw fluct), anything above y=x axis is the p-val...
bool StandardFluctuation
KS: We have two methods for Poissonian fluctuation.
std::unique_ptr< TH2D > lnLFlucHist_ProjectX
The 2D lnLHist but for ProjectionX histogram (pmu), showing (draw vs data) and (draw vs draw fluct),...
std::unique_ptr< TH1D > lnLHist_drawfluc
The lnLhist for the draw vs MC fluctuated.
bool doShapeOnly
bool whether to normalise each toy to have shape based p-value and pos pred distribution
std::vector< TH1D * > lnLHist_Sample_PredflucDraw
The histogram containing the lnL (draw vs pred fluct) for each throw for each sample.
double total_llh_data_predfluc
Data vs Fluctuated Predictive.
double llh_penalty
LLH penalty for each throw.
std::vector< TH2Poly * > MeanHist
The posterior predictive distribution in pmu cosmu using the mean.
unsigned int nChainSteps
Number of throws by user.
std::vector< int > maxBins
Max Number of Bins per each sample.
TH1D * ProjectHist(TH2D *Histogram, const bool ProjectX)
Helper to project TH2D onto axis.
double total_llh_drawfluc_predfluc
Fluctuated Draw vs Fluctuated Predictive.
double total_llh_datafluc_draw
Fluctuated Data vs Draw.
void StudyKinematicCorrelations()
KS: Study how correlated are sample or kinematic bins.
std::vector< TH1D * > lnLHist_Sample_DrawflucDraw
The histogram containing the lnL (draw vs draw fluct) for each throw for each sample.
virtual std::string GetSampleName(int Sample) const =0
MaCh3Modes * GetMaCh3Modes() const
Return pointer to MaCh3 modes.
double GetTestStatLLH(double data, double mc) const
Calculate test statistic for a single bin using Poisson.
kInfCrit
KS: Different Information Criterion tests mostly based Gelman paper.
Definition: SampleSummary.h:10
@ kWAIC
Watanabe-Akaike information criterion.
Definition: SampleSummary.h:13
@ kInfCrits
This only enumerates.
Definition: SampleSummary.h:14
@ kBIC
Bayesian Information Criterion.
Definition: SampleSummary.h:11
@ kDIC
Deviance Information Criterion.
Definition: SampleSummary.h:12
int int_t
Definition: Core.h:29
void PrintProgressBar(const Long64_t Done, const Long64_t All)
KS: Simply print progress bar.
Definition: Monitor.cpp:212