MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
SampleHandlerFD.cpp
Go to the documentation of this file.
1#include "SampleHandlerFD.h"
4#include <cstddef>
5
6#include <algorithm>
7#include <memory>
8
9// ************************************************
10SampleHandlerFD::SampleHandlerFD(std::string ConfigFileName, ParameterHandlerGeneric* xsec_cov,
11 const std::shared_ptr<OscillationHandler>& OscillatorObj_) : SampleHandlerBase() {
12// ************************************************
13 MACH3LOG_INFO("-------------------------------------------------------------------");
14 MACH3LOG_INFO("Creating SampleHandlerFD object");
15
16 //ETA - safety feature so you can't pass a NULL xsec_cov
17 if(!xsec_cov){
18 MACH3LOG_ERROR("You've passed me a nullptr to a SystematicHandlerGeneric... I need this to setup splines!");
19 throw MaCh3Exception(__FILE__, __LINE__);
20 }
21 ParHandler = xsec_cov;
22
23 nSamples = 1;
25 if (OscillatorObj_ != nullptr) {
26 MACH3LOG_WARN("You have passed an Oscillator object through the constructor of a SampleHandlerFD object - this will be used for all oscillation channels");
27 Oscillator = OscillatorObj_;
28 }
29 dathist = nullptr;
30 dathist2d = nullptr;
31 _hPDF1D = nullptr;
32 _hPDF2D = nullptr;
33
34 KinematicParameters = nullptr;
36 KinematicVectors = nullptr;
38
39 SampleHandlerFD_array = nullptr;
40 SampleHandlerFD_data = nullptr;
42 SampleName = "";
43 SampleManager = std::unique_ptr<manager>(new manager(ConfigFileName.c_str()));
44
45 // Variables related to MC stat
46 FirstTimeW2 = true;
47 UpdateW2 = false;
48}
49
51 MACH3LOG_DEBUG("I'm deleting SampleHandlerFD");
52
53 for (unsigned int yBin=0;yBin<(Binning.YBinEdges.size()-1);yBin++) {
54 if(SampleHandlerFD_array != nullptr){delete[] SampleHandlerFD_array[yBin];}
55 delete[] SampleHandlerFD_array_w2[yBin];
56 //ETA - there is a chance that you haven't added any data...
57 if(SampleHandlerFD_data != nullptr){delete[] SampleHandlerFD_data[yBin];}
58 }
59
60 if(SampleHandlerFD_array != nullptr){delete[] SampleHandlerFD_array;}
62 //ETA - there is a chance that you haven't added any data...
63 if(SampleHandlerFD_data != nullptr){delete[] SampleHandlerFD_data;}
64
65 if(THStackLeg != nullptr) delete THStackLeg;
66
67 if(dathist != nullptr) delete dathist;
68 if(dathist2d != nullptr) delete dathist2d;
69 if(_hPDF1D != nullptr) delete _hPDF1D;
70 if(_hPDF2D != nullptr) delete _hPDF2D;
71}
72
74{
75 auto ModeName = Get<std::string>(SampleManager->raw()["MaCh3ModeConfig"], __FILE__ , __LINE__);
76 Modes = std::make_unique<MaCh3Modes>(ModeName);
77 //SampleTitle has to be provided in the sample yaml otherwise this will throw an exception
78 SampleTitle = Get<std::string>(SampleManager->raw()["SampleTitle"], __FILE__ , __LINE__);
79 //SampleName has to be provided in the sample yaml otherwise this will throw an exception
80 SampleName = Get<std::string>(SampleManager->raw()["SampleName"], __FILE__ , __LINE__);
81
82 fTestStatistic = static_cast<TestStatistic>(SampleManager->GetMCStatLLH());
83 if (CheckNodeExists(SampleManager->raw(), "LikelihoodOptions")) {
84 UpdateW2 = GetFromManager<bool>(SampleManager->raw()["LikelihoodOptions"]["UpdateW2"], false);
85 }
86 //Binning
87 nDimensions = 0;
88 XVarStr = GetFromManager(SampleManager->raw()["Binning"]["XVarStr"], std::string(""));
89 Binning.XBinEdges = GetFromManager(SampleManager->raw()["Binning"]["XVarBins"], std::vector<double>());
90 const auto& edgesx = Binning.XBinEdges;
91 if (!std::is_sorted(edgesx.begin(), edgesx.end())) {
92 MACH3LOG_ERROR("XVarBins must be in increasing order in sample config {}\n XVarBins: [{}]",
93 GetTitle(), fmt::join(edgesx, ", "));
94 throw MaCh3Exception(__FILE__, __LINE__);
95 }
96 if(XVarStr.length() > 0){
98 } else{
99 MACH3LOG_ERROR("Please specify an X-variable string in sample config {}", SampleManager->GetFileName());
100 throw MaCh3Exception(__FILE__, __LINE__);
101 }
102
103 YVarStr = GetFromManager(SampleManager->raw()["Binning"]["YVarStr"], std::string(""));
104 Binning.YBinEdges = GetFromManager(SampleManager->raw()["Binning"]["YVarBins"], std::vector<double>());
105 const auto& edgesy = Binning.YBinEdges;
106 if (!std::is_sorted(edgesy.begin(), edgesy.end())) {
107 MACH3LOG_ERROR("YBinEdges must be in increasing order in sample config {}\n YBinEdges: [{}]",
108 GetTitle(), fmt::join(edgesy, ", "));
109 throw MaCh3Exception(__FILE__, __LINE__);
110 }
111 if(YVarStr.length() > 0){
112 if(XVarStr.length() == 0){
113 MACH3LOG_ERROR("Please specify an X-variable string in sample config {}. I won't work only with a Y-variable", SampleManager->GetFileName());
114 throw MaCh3Exception(__FILE__, __LINE__);
115 }
116 nDimensions++;
117 }
118
119 if(nDimensions == 0){
120 MACH3LOG_ERROR("Error setting up the sample binning");
121 MACH3LOG_ERROR("Number of dimensions is {}", nDimensions);
122 MACH3LOG_ERROR("Check that an XVarStr has been given in the sample config");
123 throw MaCh3Exception(__FILE__, __LINE__);
124 } else{
125 MACH3LOG_INFO("Found {} dimensions for sample binning", nDimensions);
126 }
127
128 //Sanity check that some binning has been specified
129 if(Binning.XBinEdges.size() == 0 && Binning.YBinEdges.size() == 0){
130 MACH3LOG_ERROR("No binning specified for either X or Y of sample binning, please add some binning to the sample config {}", SampleManager->GetFileName());
131 throw MaCh3Exception(__FILE__, __LINE__);
132 }
133
134 if (!CheckNodeExists(SampleManager->raw(), "BinningFile")){
135 MACH3LOG_ERROR("BinningFile not given in for sample {}, ReturnKinematicParameterBinning will not work", GetTitle());
136 throw MaCh3Exception(__FILE__, __LINE__);
137 }
138
139 auto mtupleprefix = Get<std::string>(SampleManager->raw()["InputFiles"]["mtupleprefix"], __FILE__, __LINE__);
140 auto mtuplesuffix = Get<std::string>(SampleManager->raw()["InputFiles"]["mtuplesuffix"], __FILE__, __LINE__);
141 auto splineprefix = Get<std::string>(SampleManager->raw()["InputFiles"]["splineprefix"], __FILE__, __LINE__);
142 auto splinesuffix = Get<std::string>(SampleManager->raw()["InputFiles"]["splinesuffix"], __FILE__, __LINE__);
143
144 int NChannels = static_cast<M3::int_t>(SampleManager->raw()["SubSamples"].size());
145 OscChannels.reserve(NChannels);
146
147 for (auto const &osc_channel : SampleManager->raw()["SubSamples"]) {
148 std::string MTupleFileName = mtupleprefix+osc_channel["mtuplefile"].as<std::string>()+mtuplesuffix;
149
150 OscChannelInfo OscInfo;
151 OscInfo.flavourName = osc_channel["Name"].as<std::string>();
152 OscInfo.flavourName_Latex = osc_channel["LatexName"].as<std::string>();
153 OscInfo.InitPDG = static_cast<NuPDG>(osc_channel["nutype"].as<int>());
154 OscInfo.FinalPDG = static_cast<NuPDG>(osc_channel["oscnutype"].as<int>());
155 OscInfo.ChannelIndex = static_cast<int>(OscChannels.size());
156
157
158 OscChannels.push_back(std::move(OscInfo));
159
160 FileToInitPDGMap[MTupleFileName] = static_cast<NuPDG>(osc_channel["nutype"].as<int>());
161 FileToFinalPDGMap[MTupleFileName] = static_cast<NuPDG>(osc_channel["oscnutype"].as<int>());
162
163 mc_files.push_back(MTupleFileName);
164 spline_files.push_back(splineprefix+osc_channel["splinefile"].as<std::string>()+splinesuffix);
165 }
166
167 //Now grab the selection cuts from the manager
168 for ( auto const &SelectionCuts : SampleManager->raw()["SelectionCuts"]) {
169 auto TempBoundsVec = GetBounds(SelectionCuts["Bounds"]);
170 KinematicCut CutObj;
171 CutObj.LowerBound = TempBoundsVec[0];
172 CutObj.UpperBound = TempBoundsVec[1];
173 CutObj.ParamToCutOnIt = ReturnKinematicParameterFromString(SelectionCuts["KinematicStr"].as<std::string>());
174 MACH3LOG_INFO("Adding cut on {} with bounds {} to {}", SelectionCuts["KinematicStr"].as<std::string>(), TempBoundsVec[0], TempBoundsVec[1]);
175 StoredSelection.push_back(CutObj);
176 }
177
178 // EM: initialise the mode weight map
179 for( int iMode=0; iMode < Modes->GetNModes(); iMode++ ) {
180 _modeNomWeightMap[Modes->GetMaCh3ModeName(iMode)] = 1.0;
181 }
182
183 // EM: multiply by the nominal weight specified in the sample config file
184 if ( SampleManager->raw()["NominalWeights"] ) {
185 for( int iMode=0; iMode<Modes->GetNModes(); iMode++ ) {
186 std::string modeStr = Modes->GetMaCh3ModeName(iMode);
187 if( SampleManager->raw()["NominalWeights"][modeStr] ) {
188 double modeWeight = SampleManager->raw()["NominalWeights"][modeStr].as<double>();
189 _modeNomWeightMap[Modes->GetMaCh3ModeName(iMode)] *= modeWeight;
190 }
191 }
192 }
193
194 // EM: print em out
195 MACH3LOG_INFO(" Nominal mode weights to apply: ");
196 for(int iMode=0; iMode<Modes->GetNModes(); iMode++ ) {
197 std::string modeStr = Modes->GetMaCh3ModeName(iMode);
198 MACH3LOG_INFO(" - {}: {}", modeStr, _modeNomWeightMap.at(modeStr));
199 }
200}
201
203 //First grab all the information from your sample config via your manager
205
206 //Now initialise all the variables you will need
207 Init();
208
210
212 SetupFDMC();
213
214 MACH3LOG_INFO("=============================================");
215 MACH3LOG_INFO("Total number of events is: {}", GetNEvents());
216
218 if (OscParams.size() > 0) {
219 MACH3LOG_INFO("Setting up NuOscillator..");
220 if (Oscillator != nullptr) {
221 MACH3LOG_INFO("You have passed an OscillatorBase object through the constructor of a SampleHandlerFD object - this will be used for all oscillation channels");
222 if(Oscillator->isEqualBinningPerOscChannel() != true) {
223 MACH3LOG_ERROR("Trying to run shared NuOscillator without EqualBinningPerOscChannel, this will not work");
224 throw MaCh3Exception(__FILE__, __LINE__);
225 }
226
227 if(OscParams.size() != Oscillator->GetOscParamsSize()){
228 MACH3LOG_ERROR("Sample {} with {} has {} osc params, while shared NuOsc has {} osc params", GetTitle(), GetSampleName(),
229 OscParams.size(), Oscillator->GetOscParamsSize());
230 MACH3LOG_ERROR("This indicate misconfiguration in your Osc yaml");
231 throw MaCh3Exception(__FILE__, __LINE__);
232 }
233 } else {
235 }
237 } else{
238 MACH3LOG_WARN("Didn't find any oscillation params, thus will not enable oscillations");
239 }
240
241 MACH3LOG_INFO("Setting up Sample Binning..");
243 MACH3LOG_INFO("Setting up Splines..");
244 SetupSplines();
245 MACH3LOG_INFO("Setting up Normalisation Pointers..");
247 MACH3LOG_INFO("Setting up Functional Pointers..");
249 MACH3LOG_INFO("Setting up Weight Pointers..");
251 MACH3LOG_INFO("Setting up Kinematic Map..");
253 MACH3LOG_INFO("=======================================================");
254}
255
256// ************************************************
258// ************************************************
259 if(KinematicParameters == nullptr || ReversedKinematicParameters == nullptr) {
260 MACH3LOG_INFO("Map KinematicParameters or ReversedKinematicParameters hasn't been initialised");
261 throw MaCh3Exception(__FILE__, __LINE__);
262 }
263 // KS: Ensure maps exist correctly
264 for (const auto& pair : *KinematicParameters) {
265 const auto& key = pair.first;
266 const auto& value = pair.second;
267
268 auto it = ReversedKinematicParameters->find(value);
269 if (it == ReversedKinematicParameters->end() || it->second != key) {
270 MACH3LOG_ERROR("Mismatch found: {} -> {} but {} -> {}",
271 key, value, value, (it != ReversedKinematicParameters->end() ? it->second : "NOT FOUND"));
272 throw MaCh3Exception(__FILE__, __LINE__);
273 }
274 }
275}
276
278{
279 // DB Commented out by default - Code heading towards GetLikelihood using arrays instead of root objects
280 // Wouldn't actually need this for GetLikelihood as TH objects wouldn't be filled
281 _hPDF1D->Reset();
282 for (unsigned int yBin=0;yBin<(Binning.YBinEdges.size()-1);yBin++) {
283 for (unsigned int xBin=0;xBin<(Binning.XBinEdges.size()-1);xBin++) {
284 _hPDF1D->AddBinContent(xBin+1,SampleHandlerFD_array[yBin][xBin]);
285 }
286 }
287}
288
290{
291 // DB Commented out by default - Code heading towards GetLikelihood using arrays instead of root objects
292 // Wouldn't actually need this for GetLikelihood as TH objects wouldn't be filled
293 _hPDF2D->Reset();
294 for (unsigned int yBin=0;yBin<(Binning.YBinEdges.size()-1);yBin++) {
295 for (unsigned int xBin=0;xBin<(Binning.XBinEdges.size()-1);xBin++) {
296 _hPDF2D->SetBinContent(xBin+1,yBin+1,SampleHandlerFD_array[yBin][xBin]);
297 }
298 }
299}
300
301// ************************************************
308// ************************************************
309 MACH3LOG_INFO("Setting up Sample Binning");
310 TString histname1d = (XVarStr).c_str();
311 TString histname2d = (XVarStr+"_"+YVarStr).c_str();
312 TString histtitle = GetTitle();
313
314 //The binning here is arbitrary, now we get info from cfg so the
315 //set1DBinning and set2Dbinning calls below will make the binning
316 //to be what we actually want
317 _hPDF1D = new TH1D("h" + histname1d + GetTitle(),histtitle, 1, 0, 1);
318 dathist = new TH1D("d" + histname1d + GetTitle(),histtitle, 1, 0, 1);
319 _hPDF2D = new TH2D("h" + histname2d + GetTitle(),histtitle, 1, 0, 1, 1, 0, 1);
320 dathist2d = new TH2D("d" + histname2d + GetTitle(),histtitle, 1, 0, 1, 1, 0, 1);
321
322 _hPDF1D->GetXaxis()->SetTitle(XVarStr.c_str());
323 dathist->GetXaxis()->SetTitle(XVarStr.c_str());
324 _hPDF2D->GetXaxis()->SetTitle(XVarStr.c_str());
325 dathist2d->GetXaxis()->SetTitle(XVarStr.c_str());
326
327 //A string to store the binning for a nice print out
328 std::string XBinEdgesStr = "";
329 std::string YBinEdgesStr = "";
330
331 for(auto XBinEdge : Binning.XBinEdges){
332 XBinEdgesStr += std::to_string(XBinEdge);
333 XBinEdgesStr += ", ";
334 }
335 MACH3LOG_INFO("XBinning:");
336 MACH3LOG_INFO("{}", XBinEdgesStr);
337
338 //And now the YBin Edges
339 for(auto YBinEdge : Binning.YBinEdges){
340 YBinEdgesStr += std::to_string(YBinEdge);
341 YBinEdgesStr += ", ";
342 }
343 MACH3LOG_INFO("YBinning:");
344 MACH3LOG_INFO("{}", YBinEdgesStr);
345
346 //Check whether you are setting up 1D or 2D binning
347 if(nDimensions == 1){
348 MACH3LOG_INFO("Setting up 1D binning with {}", XVarStr);
350 }
351 else if(nDimensions == 2){
352 MACH3LOG_INFO("Setting up 2D binning with {} and {}", XVarStr, YVarStr);
354 }
355 else{
356 MACH3LOG_ERROR("Number of dimensions is not 1 or 2, this is unsupported at the moment");
357 throw MaCh3Exception(__FILE__, __LINE__);
358 }
359}
360
361// ************************************************
362bool SampleHandlerFD::IsEventSelected(const int iEvent) {
363// ************************************************
364 const int SelectionSize = static_cast<int>(Selection.size());
365 for (int iSelection = 0; iSelection < SelectionSize; ++iSelection) {
366 const auto& Cut = Selection[iSelection];
367 const double Val = ReturnKinematicParameter(Cut.ParamToCutOnIt, iEvent);
368 if ((Val < Cut.LowerBound) || (Val >= Cut.UpperBound)) {
369 return false;
370 }
371 }
372 //DB To avoid unnecessary checks, now return false rather than setting bool to true and continuing to check
373 return true;
374}
375
376// === JM Define function to check if sub-event is selected ===
377bool SampleHandlerFD::IsSubEventSelected(const std::vector<KinematicCut> &SubEventCuts, const int iEvent, const unsigned int iSubEvent, size_t nsubevents) {
378 for (unsigned int iSelection=0;iSelection < SubEventCuts.size() ;iSelection++) {
379 std::vector<double> Vec = ReturnKinematicVector(SubEventCuts[iSelection].ParamToCutOnIt, iEvent);
380 if (nsubevents != Vec.size()) {
381 MACH3LOG_ERROR("Cannot apply kinematic cut on {} as it is of different size to plotting variable");
382 throw MaCh3Exception(__FILE__, __LINE__);
383 }
384 const double Val = Vec[iSubEvent];
385 if ((Val < SubEventCuts[iSelection].LowerBound) || (Val >= SubEventCuts[iSelection].UpperBound)) {
386 return false;
387 }
388 }
389 //DB To avoid unnecessary checks, now return false rather than setting bool to true and continuing to check
390 return true;
391}
392// ===========================================================
393
394//************************************************
395// Reweight function
397//************************************************
398 //KS: Reset the histograms before reweight
400
401 //You only need to do these things if Oscillator has been initialised
402 //if not then you're not considering oscillations
403 if (Oscillator) Oscillator->Evaluate();
404
405 // Calculate weight coming from all splines if we initialised handler
406 if(SplineHandler) SplineHandler->Evaluate();
407
408 #ifdef MULTITHREAD
409 // Call entirely different routine if we're running with openMP
410 FillArray_MP();
411 #else
412 FillArray();
413 #endif
414
415 //KS: If you want to not update W2 wights then uncomment this line
416 if(!UpdateW2) FirstTimeW2 = false;
417}
418
419//************************************************
428//************************************************
429 //DB Reset which cuts to apply
431
433
434 for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
435 ApplyShifts(iEvent);
436
437 if (!IsEventSelected(iEvent)) {
438 continue;
439 }
440
441 FarDetectorCoreInfo* MCEvent = &MCSamples[iEvent];
442 M3::float_t splineweight = CalcWeightSpline(MCEvent);
443 //DB Catch negative spline weights and skip any event with a negative event. Previously we would set weight to zero and continue but that is inefficient. Do this on a spline-by-spline basis
444 if (splineweight <= 0.){
445 MCEvent->xsec_w = 0.;
446 continue;
447 }
448
449 //Loop over stored normalisation and function pointers
450 M3::float_t normweight = CalcWeightNorm(MCEvent);
451
452 // Virtual by default does nothing
453 CalcWeightFunc(iEvent);
454
455 MCEvent->xsec_w = splineweight*normweight;
456
457 //DB Total weight
458 M3::float_t totalweight = GetEventWeight(iEvent);
459 //DB Catch negative weights and skip any event with a negative event
460 if (totalweight <= 0.){
461 MCEvent->xsec_w = 0.;
462 continue;
463 }
464 //DB Switch on BinningOpt to allow different binning options to be implemented
465 //The alternative would be to have inheritance based on BinningOpt
466 const double XVar = *(MCEvent->x_var);
467
468 //DB Find the relevant bin in the PDF for each event
469 const int XBinToFill = Binning.FindXBin(XVar, MCEvent->NomXBin);
470 const int YBinToFill = MCEvent->NomYBin;
471
472 //DB Fill relevant part of thread array
473 if (XBinToFill != -1 && YBinToFill != -1) {
474 SampleHandlerFD_array[YBinToFill][XBinToFill] += totalweight;
475 if (FirstTimeW2) SampleHandlerFD_array_w2[YBinToFill][XBinToFill] += totalweight*totalweight;
476 }
477 }
478}
479
480#ifdef MULTITHREAD
481// ************************************************
483void SampleHandlerFD::FillArray_MP() {
484// ************************************************
485 //DB Reset which cuts to apply
487
489
490 //This is stored as [y][x] due to shifts only occurring in the x variable (Erec/Lep mom) - I believe this will help reduce cache misses
491 double** SampleHandlerFD_array_private = nullptr;
492 double** SampleHandlerFD_array_private_w2 = nullptr;
493 // Declare the omp parallel region
494 // The parallel region needs to stretch beyond the for loop!
495 #pragma omp parallel private(SampleHandlerFD_array_private, SampleHandlerFD_array_private_w2)
496 {
497 // private to each thread
498 // ETA - maybe we can use parallel firstprivate to initialise these?
499 SampleHandlerFD_array_private = new double*[Binning.nYBins];
500 SampleHandlerFD_array_private_w2 = new double*[Binning.nYBins];
501 for (size_t yBin=0;yBin<Binning.nYBins; ++yBin) {
502 SampleHandlerFD_array_private[yBin] = new double[Binning.nXBins];
503 SampleHandlerFD_array_private_w2[yBin] = new double[Binning.nXBins];
504
505 std::fill_n(SampleHandlerFD_array_private[yBin], Binning.nXBins, 0.0);
506 std::fill_n(SampleHandlerFD_array_private_w2[yBin], Binning.nXBins, 0.0);
507 }
508
509 //DB - Brain dump of speedup ideas
510 //
511 //Those relevant to reweighting
512 // 1. Don't bother storing and calculating NC signal events - Implemented and saves marginal s/step
513 // 2. Loop over spline event weight calculation in the following event loop - Currently done in splineSKBase->calcWeight() where multi-threading won't be optimised - Implemented and saves 0.3s/step
514 // 3. Inline getDiscVar or somehow include that calculation inside the multi-threading - Implemented and saves about 0.01s/step
515 // 4. Include isCC inside SKMCStruct so don't have to have several 'if' statements determine if oscillation weight needs to be set to 1.0 for NC events - Implemented and saves marginal s/step
516 // 5. Do explicit check on adjacent bins when finding event XBin instead of looping over all BinEdge indices - Implemented but doesn't significantly affect s/step
517 //
518 //Other aspects
519 // 1. Order minituples in Y-axis variable as this will *hopefully* reduce cache misses inside SampleHandlerFD_array_class[yBin][xBin]
520 //
521 // We will hit <0.1 s/step eventually! :D
522
523 const unsigned int NumberOfEvents = GetNEvents();
524 #pragma omp for
525 for (unsigned int iEvent = 0; iEvent < NumberOfEvents; ++iEvent) {
526 //ETA - generic functions to apply shifts to kinematic variables
527 // Apply this before IsEventSelected is called.
528 ApplyShifts(iEvent);
529
530 //ETA - generic functions to apply shifts to kinematic variable
531 //this is going to be slow right now due to string comps under the hood.
532 //Need to implement a more efficient version of event-by-event cut checks
533 if(!IsEventSelected(iEvent)){
534 continue;
535 }
536
537 //DB SKDet Syst
538 //As weights were skdet::fParProp, and we use the non-shifted erec, we might as well cache the corresponding fParProp index for each event and the pointer to it
539 FarDetectorCoreInfo* MCEvent = &MCSamples[iEvent];
540 const M3::float_t splineweight = CalcWeightSpline(MCEvent);
541 //DB Catch negative spline weights and skip any event with a negative event. Previously we would set weight to zero and continue but that is inefficient
542 if (splineweight <= 0.){
543 MCEvent->xsec_w = 0.;
544 continue;
545 }
546
547 const M3::float_t normweight = CalcWeightNorm(MCEvent);
548
549 // Virtual by default does nothing
550 CalcWeightFunc(iEvent);
551
552 MCEvent->xsec_w = splineweight*normweight;
553
554 const M3::float_t totalweight = GetEventWeight(iEvent);
555
556 //DB Catch negative weights and skip any event with a negative event
557 if (totalweight <= 0.){
558 MCEvent->xsec_w = 0.;
559 continue;
560 }
561
562 //DB Switch on BinningOpt to allow different binning options to be implemented
563 //The alternative would be to have inheritance based on BinningOpt
564 const double XVar = (*(MCEvent->x_var));
565
566 //DB Find the relevant bin in the PDF for each event
567 const int XBinToFill = Binning.FindXBin(XVar, MCEvent->NomXBin);
568 const int YBinToFill = MCEvent->NomYBin;
569
570 //ETA - we can probably remove this final if check on the -1?
571 //Maybe we can add an overflow bin to the array and assign any events to this bin?
572 //Might save us an extra if call?
573 //DB Fill relevant part of thread array
574 if (XBinToFill != -1 && YBinToFill != -1) {
575 SampleHandlerFD_array_private[YBinToFill][XBinToFill] += totalweight;
576 SampleHandlerFD_array_private_w2[YBinToFill][XBinToFill] += totalweight*totalweight;
577 }
578 }
579 //End of Calc Weights and fill Array
580 //==================================================
581 // DB Copy contents of 'SampleHandlerFD_array_private' into 'SampleHandlerFD_array' which can then be used in GetLikelihood
582 for (size_t yBin = 0; yBin < Binning.nYBins; ++yBin) {
583 for (size_t xBin = 0; xBin < Binning.nXBins; ++xBin) {
584 #pragma omp atomic
585 SampleHandlerFD_array[yBin][xBin] += SampleHandlerFD_array_private[yBin][xBin];
586 if(FirstTimeW2) {
587 #pragma omp atomic
588 SampleHandlerFD_array_w2[yBin][xBin] += SampleHandlerFD_array_private_w2[yBin][xBin];
589 }
590 }
591 }
592
593 for (size_t yBin = 0; yBin < Binning.nYBins; ++yBin) {
594 delete[] SampleHandlerFD_array_private[yBin];
595 delete[] SampleHandlerFD_array_private_w2[yBin];
596 }
597 delete[] SampleHandlerFD_array_private;
598 delete[] SampleHandlerFD_array_private_w2;
599 } //end of parallel region
600}
601#endif
602
603// **************************************************
604// Helper function to reset the data and MC histograms
606// **************************************************
607 //DB Reset values stored in PDF array to 0.
608 // Don't openMP this; no significant gain
609 for (size_t yBin = 0; yBin < Binning.nYBins; ++yBin) {
610 #ifdef MULTITHREAD
611 #pragma omp simd
612 #endif
613 for (size_t xBin = 0; xBin < Binning.nXBins; ++xBin) {
614 SampleHandlerFD_array[yBin][xBin] = 0.;
615 if(FirstTimeW2) SampleHandlerFD_array_w2[yBin][xBin] = 0.;
616 }
617 }
618} // end function
619
620void SampleHandlerFD::RegisterIndividualFunctionalParameter(const std::string& fpName, int fpEnum, FuncParFuncType fpFunc){
621 // Add protections to not add the same functional parameter twice
622 if (funcParsNamesMap.find(fpName) != funcParsNamesMap.end()) {
623 MACH3LOG_ERROR("Functional parameter {} already registered in funcParsNamesMap with enum {}", fpName, funcParsNamesMap[fpName]);
624 throw MaCh3Exception(__FILE__, __LINE__);
625 }
626 if (std::find(funcParsNamesVec.begin(), funcParsNamesVec.end(), fpName) != funcParsNamesVec.end()) {
627 MACH3LOG_ERROR("Functional parameter {} already in funcParsNamesVec", fpName);
628 throw MaCh3Exception(__FILE__, __LINE__);
629 }
630 if (funcParsFuncMap.find(fpEnum) != funcParsFuncMap.end()) {
631 MACH3LOG_ERROR("Functional parameter enum {} already registered in funcParsFuncMap", fpEnum);
632 throw MaCh3Exception(__FILE__, __LINE__);
633 }
634 funcParsNamesMap[fpName] = fpEnum;
635 funcParsNamesVec.push_back(fpName);
636 funcParsFuncMap[fpEnum] = fpFunc;
637}
638
641 // RegisterFunctionalParameters is implemented in experiment-specific code,
642 // which calls RegisterIndividualFuncPar to populate funcParsNamesMap, funcParsNamesVec, and funcParsFuncMap
644 funcParsMap.resize(funcParsNamesMap.size());
645 funcParsGrid.resize(GetNEvents());
646
647 // For every functional parameter in XsecCov that matches the name in funcParsNames, add it to the map
648 for (FunctionalParameter & fp : funcParsVec) {
649 for (std::string name : funcParsNamesVec) {
650 if (fp.name == name) {
651 MACH3LOG_INFO("Adding functional parameter: {} to funcParsMap with key: {}", fp.name, funcParsNamesMap[fp.name]);
652 fp.funcPtr = &funcParsFuncMap[funcParsNamesMap[fp.name]];
653 funcParsMap[static_cast<std::size_t>(funcParsNamesMap[fp.name])] = &fp;
654 continue;
655 }
656 }
657 // If we don't find a match, we need to throw an error
658 if (funcParsMap[static_cast<std::size_t>(funcParsNamesMap[fp.name])] == nullptr) {
659 MACH3LOG_ERROR("Functional parameter {} not found, did you define it in RegisterFunctionalParameters()?", fp.name);
660 throw MaCh3Exception(__FILE__, __LINE__);
661 }
662 }
663
664 // Mostly the same as CalcXsecNormsBins
665 // For each event, make a vector of pointers to the functional parameters
666 for (std::size_t iEvent = 0; iEvent < static_cast<std::size_t>(GetNEvents()); ++iEvent) {
667 // Now loop over the functional parameters and get a vector of enums corresponding to the functional parameters
668 for (std::vector<FunctionalParameter>::iterator it = funcParsVec.begin(); it != funcParsVec.end(); ++it) {
669 // Check whether the interaction modes match
670 bool ModeMatch = MatchCondition((*it).modes, static_cast<int>(std::round(*(MCSamples[iEvent].mode))));
671 if (!ModeMatch) {
672 MACH3LOG_TRACE("Event {}, missed Mode check ({}) for dial {}", iEvent, *(MCSamples[iEvent].mode), (*it).name);
673 continue;
674 }
675 // Now check whether within kinematic bounds
676 bool IsSelected = true;
677 if ((*it).hasKinBounds) {
678 const auto& kinVars = (*it).KinematicVarStr;
679 const auto& selection = (*it).Selection;
680
681 for (std::size_t iKinPar = 0; iKinPar < kinVars.size(); ++iKinPar) {
682 const double kinVal = ReturnKinematicParameter(kinVars[iKinPar], static_cast<int>(iEvent));
683
684 bool passedAnyBound = false;
685 const auto& boundsList = selection[iKinPar];
686
687 for (const auto& bounds : boundsList) {
688 if (kinVal > bounds[0] && kinVal <= bounds[1]) {
689 passedAnyBound = true;
690 break;
691 }
692 }
693
694 if (!passedAnyBound) {
695 MACH3LOG_TRACE("Event {}, missed kinematic check ({}) for dial {}",
696 iEvent, kinVars[iKinPar], (*it).name);
697 IsSelected = false;
698 break;
699 }
700 }
701 }
702 // Need to then break the event loop
703 if(!IsSelected){
704 MACH3LOG_TRACE("Event {}, missed Kinematic var check for dial {}", iEvent, (*it).name);
705 continue;
706 }
707 auto funcparenum = funcParsNamesMap[(*it).name];
708 funcParsGrid.at(iEvent).push_back(funcparenum);
709 }
710 }
711 MACH3LOG_INFO("Finished setting up functional parameters");
712}
713
715 // Given a sample and event, apply the shifts to the event based on the vector of functional parameter enums
716 // First reset shifted array back to nominal values
717 resetShifts(iEvent);
718
719 const std::vector<int>& fpEnums = funcParsGrid[iEvent];
720 const std::size_t nShifts = fpEnums.size();
721
722 for (std::size_t i = 0; i < nShifts; ++i) {
723 const int fpEnum = fpEnums[i];
724 FunctionalParameter *fp = funcParsMap[static_cast<std::size_t>(fpEnum)];
725 // if (fp->funcPtr) {
726 // (*fp->funcPtr)(fp->valuePtr, iEvent);
727 // } else {
728 // MACH3LOG_ERROR("Functional parameter function pointer for {} is null for event {} in sample {}", fp->name, iEvent);
729 // throw MaCh3Exception(__FILE__, __LINE__);
730 // }
731 (*fp->funcPtr)(fp->valuePtr, iEvent);
732 }
733}
734
735// ***************************************************************************
736// Calculate the spline weight for one event
738// ***************************************************************************
739 M3::float_t spline_weight = 1.0;
740 const int nSplines = static_cast<int>(MCEvent->xsec_spline_pointers.size());
741 //DB Xsec syst
742 //Loop over stored spline pointers
743 #ifdef MULTITHREAD
744 #pragma omp simd
745 #endif
746 for (int iSpline = 0; iSpline < nSplines; ++iSpline) {
747 spline_weight *= *(MCEvent->xsec_spline_pointers[iSpline]);
748 }
749 return spline_weight;
750}
751
752// ***************************************************************************
753// Calculate the normalisation weight for an event
755// ***************************************************************************
756 M3::float_t xsecw = 1.0;
757 const int nNorms = static_cast<int>(MCEvent->xsec_norm_pointers.size());
758 //Loop over stored normalisation and function pointers
759 #ifdef MULTITHREAD
760 #pragma omp simd
761 #endif
762 for (int iParam = 0; iParam < nNorms; ++iParam)
763 {
764#pragma GCC diagnostic push
765#pragma GCC diagnostic ignored "-Wuseless-cast"
766 xsecw *= static_cast<M3::float_t>(*(MCEvent->xsec_norm_pointers[iParam]));
767#pragma GCC diagnostic pop
768 #ifdef DEBUG
769 if (std::isnan(xsecw)) MACH3LOG_WARN("iParam= {} xsecweight=nan from norms", iParam);
770 #endif
771 }
772 return xsecw;
773}
774
775// ***************************************************************************
776// Setup the norm parameters
778// ***************************************************************************
779 std::vector< std::vector< int > > xsec_norms_bins(GetNEvents());
780
781 std::vector<NormParameter> norm_parameters = ParHandler->GetNormParsFromSampleName(GetSampleName());
782
783 if(!ParHandler){
784 MACH3LOG_ERROR("ParHandler is not setup!");
785 throw MaCh3Exception(__FILE__ , __LINE__ );
786 }
787
788 // Assign xsec norm bins in MCSamples tree
789 CalcNormsBins(norm_parameters, xsec_norms_bins);
790
791 //DB Attempt at reducing impact of SystematicHandlerGeneric::calcReweight()
792 int counter;
793 for (unsigned int iEvent = 0; iEvent < GetNEvents(); ++iEvent) {
794 counter = 0;
795
796 MCSamples[iEvent].xsec_norm_pointers.resize(xsec_norms_bins[iEvent].size());
797 for(auto const & norm_bin: xsec_norms_bins[iEvent]) {
798 MCSamples[iEvent].xsec_norm_pointers[counter] = ParHandler->RetPointer(norm_bin);
799 counter += 1;
800 }
801 }
802}
803
804// ************************************************
805//A way to check whether a normalisation parameter applies to an event or not
806void SampleHandlerFD::CalcNormsBins(std::vector<NormParameter>& norm_parameters, std::vector< std::vector< int > >& xsec_norms_bins) {
807// ************************************************
808 #ifdef DEBUG
809 std::vector<int> VerboseCounter(norm_parameters.size(), 0);
810 #endif
811 for(unsigned int iEvent = 0; iEvent < GetNEvents(); ++iEvent){
812 std::vector< int > NormBins = {};
813 if (ParHandler) {
814 // Skip oscillated NC events
815 // Not strictly needed, but these events don't get included in oscillated predictions, so
816 // no need to waste our time calculating and storing information about xsec parameters
817 // that will never be used.
818 if (MCSamples[iEvent].isNC && (*MCSamples[iEvent].nupdg != *MCSamples[iEvent].nupdgUnosc) ) {
819 MACH3LOG_TRACE("Event {}, missed NC/signal check", iEvent);
820 continue;
821 } //DB Abstract check on MaCh3Modes to determine which apply to neutral current
822 for (std::vector<NormParameter>::iterator it = norm_parameters.begin(); it != norm_parameters.end(); ++it) {
823 //Now check that the target of an interaction matches with the normalisation parameters
824 bool TargetMatch = MatchCondition((*it).targets, *(MCSamples[iEvent].Target));
825 if (!TargetMatch) {
826 MACH3LOG_TRACE("Event {}, missed target check ({}) for dial {}", iEvent, *(MCSamples[iEvent].Target), (*it).name);
827 continue;
828 }
829
830 //Now check that the neutrino flavour in an interaction matches with the normalisation parameters
831 bool FlavourMatch = MatchCondition((*it).pdgs, *(MCSamples[iEvent].nupdg));
832 if (!FlavourMatch) {
833 MACH3LOG_TRACE("Event {}, missed PDG check ({}) for dial {}", iEvent,(*MCSamples[iEvent].nupdg), (*it).name);
834 continue;
835 }
836
837 //Now check that the unoscillated neutrino flavour in an interaction matches with the normalisation parameters
838 bool FlavourUnoscMatch = MatchCondition((*it).preoscpdgs, *(MCSamples[iEvent].nupdgUnosc));
839 if (!FlavourUnoscMatch){
840 MACH3LOG_TRACE("Event {}, missed FlavourUnosc check ({}) for dial {}", iEvent,(*MCSamples[iEvent].nupdgUnosc), (*it).name);
841 continue;
842 }
843
844 //Now check that the mode of an interaction matches with the normalisation parameters
845 bool ModeMatch = MatchCondition((*it).modes, static_cast<int>(std::round(*(MCSamples[iEvent].mode))));
846 if (!ModeMatch) {
847 MACH3LOG_TRACE("Event {}, missed Mode check ({}) for dial {}", iEvent, *(MCSamples[iEvent].mode), (*it).name);
848 continue;
849 }
850
851 //Now check whether the norm has kinematic bounds
852 //i.e. does it only apply to events in a particular kinematic region?
853 // Now check whether within kinematic bounds
854 bool IsSelected = true;
855 if ((*it).hasKinBounds) {
856 const auto& kinVars = (*it).KinematicVarStr;
857 const auto& selection = (*it).Selection;
858
859 for (std::size_t iKinPar = 0; iKinPar < kinVars.size(); ++iKinPar) {
860 const double kinVal = ReturnKinematicParameter(kinVars[iKinPar], static_cast<int>(iEvent));
861
862 bool passedAnyBound = false;
863 const auto& boundsList = selection[iKinPar];
864
865 for (const auto& bounds : boundsList) {
866 if (kinVal > bounds[0] && kinVal <= bounds[1]) {
867 passedAnyBound = true;
868 break;
869 }
870 }
871
872 if (!passedAnyBound) {
873 MACH3LOG_TRACE("Event {}, missed kinematic check ({}) for dial {}",
874 iEvent, kinVars[iKinPar], (*it).name);
875 IsSelected = false;
876 break;
877 }
878 }
879 }
880 // Need to then break the event loop
881 if(!IsSelected){
882 MACH3LOG_TRACE("Event {}, missed Kinematic var check for dial {}", iEvent, (*it).name);
883 continue;
884 }
885 // Now set 'index bin' for each normalisation parameter
886 // All normalisations are just 1 bin for 2015, so bin = index (where index is just the bin for that normalisation)
887 int bin = (*it).index;
888
889 NormBins.push_back(bin);
890 MACH3LOG_TRACE("Event {}, will be affected by dial {}", iEvent, (*it).name);
891 #ifdef DEBUG
892 VerboseCounter[std::distance(norm_parameters.begin(), it)]++;
893 #endif
894 //}
895 } // end iteration over norm_parameters
896 } // end if (ParHandler)
897 xsec_norms_bins[iEvent] = NormBins;
898 }//end loop over events
899 #ifdef DEBUG
900 MACH3LOG_DEBUG("┌──────────────────────────────────────────────────────────┐");
901 for (std::size_t i = 0; i < norm_parameters.size(); ++i) {
902 const auto& norm = norm_parameters[i];
903 double eventRatio = static_cast<double>(VerboseCounter[i]) / static_cast<double>(GetNEvents());
904
905 MACH3LOG_DEBUG("│ Param {:<15}, affects {:<8} events ({:>6.2f}%) │",
906 ParHandler->GetParFancyName(norm.index), VerboseCounter[i], eventRatio);
907 }
908 MACH3LOG_DEBUG("└──────────────────────────────────────────────────────────┘");
909 #endif
910}
911
912// ************************************************
913void SampleHandlerFD::SetupReweightArrays(const size_t numberXBins, const size_t numberYBins) {
914// ************************************************
915 //Set the number of X and Y bins now
916 Binning.nXBins = numberXBins;
917 Binning.nYBins = numberYBins;
918
919 SampleHandlerFD_array = new double*[Binning.nYBins];
921 SampleHandlerFD_data = new double*[Binning.nYBins];
922 for (size_t yBin=0;yBin<Binning.nYBins;yBin++) {
923 SampleHandlerFD_array[yBin] = new double[Binning.nXBins];
924 SampleHandlerFD_array_w2[yBin] = new double[Binning.nXBins];
925 SampleHandlerFD_data[yBin] = new double[Binning.nXBins];
926 for (size_t xBin=0;xBin<Binning.nXBins;xBin++) {
927 SampleHandlerFD_array[yBin][xBin] = 0.;
928 SampleHandlerFD_array_w2[yBin][xBin] = 0.;
929 SampleHandlerFD_data[yBin][xBin] = 0.;
930 }
931 }
932}
933
934//ETA
935//New versions of set binning functions is SampleHandlerBase
936//so that we can set the values of the bin and lower/upper
937//edges in the skmc_base. Hopefully we can use this to make
938//fill1Dhist and fill2Dhist quicker
939void SampleHandlerFD::Set1DBinning(size_t nbins, double* boundaries)
940{
941 _hPDF1D->Reset();
942 _hPDF1D->SetBins(static_cast<int>(nbins),boundaries);
943 dathist->SetBins(static_cast<int>(nbins),boundaries);
944
945 Binning.YBinEdges = std::vector<double>(2);
946 Binning.YBinEdges[0] = -1e8;
947 Binning.YBinEdges[1] = 1e8;
948
949 double YBinEdges_Arr[2];
950 YBinEdges_Arr[0] = Binning.YBinEdges[0];
951 YBinEdges_Arr[1] = Binning.YBinEdges[1];
952
953 _hPDF2D->Reset();
954 _hPDF2D->SetBins(static_cast<int>(nbins),boundaries,1,YBinEdges_Arr);
955 dathist2d->SetBins(static_cast<int>(nbins),boundaries,1,YBinEdges_Arr);
956
957 //Set the number of X and Y bins now
958 SetupReweightArrays(Binning.XBinEdges.size() - 1, Binning.YBinEdges.size() - 1);
959
961}
962
964 for(unsigned int event_i = 0; event_i < GetNEvents(); event_i++){
965 //Set x_var and y_var values based on XVarStr and YVarStr
966 MCSamples[event_i].x_var = GetPointerToKinematicParameter(XVarStr, event_i);
967 if (std::isnan(*MCSamples[event_i].x_var) || std::isinf(*MCSamples[event_i].x_var)) {
968 MACH3LOG_ERROR("X var for event {} is ill-defined and equal to {}", event_i, *MCSamples[event_i].x_var);
969 throw MaCh3Exception(__FILE__, __LINE__);
970 }
971
972 //Give y_var M3::_BAD_DOUBLE_ value for the 1D case since this won't be used
973 MCSamples[event_i].y_var = &(M3::_BAD_DOUBLE_);
974 int bin = _hPDF1D->FindBin(*(MCSamples[event_i].x_var));
975
976 if ((bin-1) >= 0 && (bin-1) < int(Binning.XBinEdges.size()-1)) {
977 MCSamples[event_i].NomXBin = bin-1;
978 } else {
979 MCSamples[event_i].NomXBin = -1;
980 }
981 MCSamples[event_i].NomYBin = 0;
982 }
983
988 //Set rw_pdf_bin and rw_upper_xbinedge and rw_lower_xbinedge for each skmc_base
989 for(size_t bin_x = 0; bin_x < Binning.nXBins; bin_x++){
990 double low_lower_edge = M3::_DEFAULT_RETURN_VAL_;
991 double low_edge = Binning.XBinEdges[bin_x];
992 double upper_edge = Binning.XBinEdges[bin_x+1];
993 double upper_upper_edge = M3::_DEFAULT_RETURN_VAL_;
994
995 if (bin_x == 0) {
996 low_lower_edge = Binning.XBinEdges[0];
997 } else {
998 low_lower_edge = Binning.XBinEdges[bin_x-1];
999 }
1000
1001 if (bin_x + 2 < Binning.nXBins) {
1002 upper_upper_edge = Binning.XBinEdges[bin_x + 2];
1003 } else if (bin_x + 1 < Binning.nXBins) {
1004 upper_upper_edge = Binning.XBinEdges[bin_x + 1];
1005 }
1006
1007 Binning.rw_lower_xbinedge[bin_x] = low_edge;
1008 Binning.rw_upper_xbinedge[bin_x] = upper_edge;
1009 Binning.rw_lower_lower_xbinedge[bin_x] = low_lower_edge;
1010 Binning.rw_upper_upper_xbinedge[bin_x] = upper_upper_edge;
1011 }
1012}
1013
1014void SampleHandlerFD::Set2DBinning(size_t nbins1, double* boundaries1, size_t nbins2, double* boundaries2)
1015{
1016 _hPDF1D->Reset();
1017 _hPDF1D->SetBins(static_cast<int>(nbins1),boundaries1);
1018 dathist->SetBins(static_cast<int>(nbins1),boundaries1);
1019
1020 _hPDF2D->Reset();
1021 _hPDF2D->SetBins(static_cast<int>(nbins1),boundaries1,static_cast<int>(nbins2),boundaries2);
1022 dathist2d->SetBins(static_cast<int>(nbins1),boundaries1,static_cast<int>(nbins2),boundaries2);
1023
1024 //Set the number of X and Y bins now
1025 SetupReweightArrays(Binning.XBinEdges.size() - 1, Binning.YBinEdges.size() - 1);
1026
1028}
1029
1030// ************************************************
1032// ************************************************
1033 for(unsigned int event_i = 0 ; event_i < GetNEvents(); event_i++) {
1034 //Set x_var and y_var values based on XVarStr and YVarStr
1035 MCSamples[event_i].x_var = GetPointerToKinematicParameter(XVarStr, event_i);
1036 MCSamples[event_i].y_var = GetPointerToKinematicParameter(YVarStr, event_i);
1037
1038 if (std::isnan(*MCSamples[event_i].x_var) || std::isinf(*MCSamples[event_i].x_var)) {
1039 MACH3LOG_ERROR("X var for event {} is ill-defined and equal to {}", event_i, *MCSamples[event_i].x_var);
1040 throw MaCh3Exception(__FILE__, __LINE__);
1041 }
1042 if (std::isnan(*MCSamples[event_i].y_var) || std::isinf(*MCSamples[event_i].y_var)) {
1043 MACH3LOG_ERROR("Y var for event {} is ill-defined and equal to {}", event_i, *MCSamples[event_i].y_var);
1044 throw MaCh3Exception(__FILE__, __LINE__);
1045 }
1046 //Global bin number
1047 int bin = _hPDF2D->FindBin(*(MCSamples[event_i].x_var), *(MCSamples[event_i].y_var));
1048
1049 int bin_x = M3::_BAD_INT_;
1050 int bin_y = M3::_BAD_INT_;
1051 int bin_z = M3::_BAD_INT_;
1052 _hPDF2D->GetBinXYZ(bin, bin_x, bin_y, bin_z);
1053
1054 if ((bin_x-1) >= 0 && (bin_x-1) < int(Binning.XBinEdges.size()-1)) {
1055 MCSamples[event_i].NomXBin = bin_x-1;
1056 }
1057 MCSamples[event_i].NomYBin = bin_y-1;
1058 if(MCSamples[event_i].NomYBin < 0){
1059 MACH3LOG_WARN("Nominal YBin PROBLEM, y-bin is {}", MCSamples[event_i].NomYBin);
1060 }
1061 }
1062
1067 //Set rw_pdf_bin and rw_upper_xbinedge and rw_lower_xbinedge for each skmc_base
1068 for(size_t bin_x = 0; bin_x < Binning.nXBins; bin_x++){
1069 double low_lower_edge = M3::_DEFAULT_RETURN_VAL_;
1070 double low_edge = Binning.XBinEdges[bin_x];
1071 double upper_edge = Binning.XBinEdges[bin_x+1];
1072 double upper_upper_edge = M3::_DEFAULT_RETURN_VAL_;
1073
1074 if (bin_x == 0) {
1075 low_lower_edge = Binning.XBinEdges[0];
1076 } else {
1077 low_lower_edge = Binning.XBinEdges[bin_x-1];
1078 }
1079
1080 if (bin_x + 2 < Binning.nXBins) {
1081 upper_upper_edge = Binning.XBinEdges[bin_x + 2];
1082 } else if (bin_x + 1 < Binning.nXBins) {
1083 upper_upper_edge = Binning.XBinEdges[bin_x + 1];
1084 }
1085
1086 Binning.rw_lower_xbinedge[bin_x] = low_edge;
1087 Binning.rw_upper_xbinedge[bin_x] = upper_edge;
1088 Binning.rw_lower_lower_xbinedge[bin_x] = low_lower_edge;
1089 Binning.rw_upper_upper_xbinedge[bin_x] = upper_upper_edge;
1090 }
1091}
1092
1093// ************************************************
1094TH1* SampleHandlerFD::GetW2Hist(const int Dimension) {
1095// ************************************************
1096 if(Dimension == 1) {
1097 TH1D* W2Hist = dynamic_cast<TH1D*>(_hPDF1D->Clone((_hPDF1D->GetName() + std::string("_W2")).c_str()));
1098 W2Hist->Reset();
1099 for (unsigned int yBin = 0; yBin < (Binning.YBinEdges.size()-1); yBin++) {
1100 for (unsigned int xBin = 0; xBin < (Binning.XBinEdges.size()-1); xBin++) {
1101 W2Hist->AddBinContent(xBin+1, SampleHandlerFD_array_w2[yBin][xBin]);
1102 }
1103 }
1104 return W2Hist;
1105 } else if(Dimension == 2) {
1106 TH2D* W2Hist = dynamic_cast<TH2D*>(_hPDF2D->Clone((_hPDF2D->GetName() + std::string("_W2")).c_str()));
1107 W2Hist->Reset();
1108 for (unsigned int yBin = 0; yBin < (Binning.YBinEdges.size()-1); yBin++) {
1109 for (unsigned int xBin = 0; xBin < (Binning.XBinEdges.size()-1); xBin++) {
1110 W2Hist->SetBinContent(xBin+1, yBin+1, SampleHandlerFD_array_w2[yBin][xBin]);
1111 }
1112 }
1113 return W2Hist;
1114 } else{
1115 MACH3LOG_ERROR("Asking for {} with N Dimension = {}. This is not implemented", __func__, Dimension);
1116 throw MaCh3Exception(__FILE__, __LINE__);
1117 }
1118}
1119
1120// ************************************************
1121TH1* SampleHandlerFD::GetMCHist(const int Dimension) {
1122// ************************************************
1123 if(Dimension == 1) {
1124 Fill1DHist();
1125 return _hPDF1D;
1126 } else if(Dimension == 2) {
1127 Fill2DHist();
1128 return _hPDF2D;
1129 } else{
1130 MACH3LOG_ERROR("Asdking for {} with N Dimension = {}. This is not implemented", __func__, Dimension);
1131 throw MaCh3Exception(__FILE__, __LINE__);
1132 }
1133}
1134
1135// ************************************************
1136TH1* SampleHandlerFD::GetDataHist(const int Dimension) {
1137// ************************************************
1138 if(Dimension == 1) {
1139 return dathist;
1140 } else if(Dimension == 2) {
1141 return dathist2d;
1142 } else{
1143 MACH3LOG_ERROR("Asdking for {} with N Dimension = {}. This is not implemented", __func__, Dimension);
1144 throw MaCh3Exception(__FILE__, __LINE__);
1145 }
1146}
1147
1148void SampleHandlerFD::AddData(std::vector<double> &data) {
1149 dathist2d = nullptr;
1150 dathist->Reset();
1151
1152 if (GetNDim()!=1) {
1153 MACH3LOG_ERROR("Trying to set a 1D 'data' histogram when the number of dimensions for this sample is {}", GetNDim());
1154 MACH3LOG_ERROR("This won't work, please specify the correct dimensions in your sample config with the X and Y variables");
1155 throw MaCh3Exception(__FILE__, __LINE__);
1156 }
1157
1158 for (auto const& data_point : data){
1159 dathist->Fill(data_point);
1160 }
1161
1162 if(SampleHandlerFD_data == nullptr) {
1163 MACH3LOG_ERROR("SampleHandlerFD_data haven't been initialised yet");
1164 throw MaCh3Exception(__FILE__, __LINE__);
1165 }
1166 for (size_t yBin=0;yBin<Binning.nYBins;yBin++) {
1167 for (size_t xBin=0;xBin<Binning.nXBins;xBin++) {
1168 SampleHandlerFD_data[yBin][xBin] = dathist->GetBinContent(static_cast<int>(xBin+1));
1169 }
1170 }
1171}
1172
1173void SampleHandlerFD::AddData(std::vector< std::vector <double> > &data) {
1174 dathist = nullptr;
1175 dathist2d->Reset();
1176
1177 if (GetNDim()!=2) {
1178 MACH3LOG_ERROR("Trying to set a 2D 'data' histogram when the number of dimensions for this sample is {}", GetNDim());
1179 MACH3LOG_ERROR("This won't work, please specify the correct dimensions in your sample config with the X and Y variables");
1180 throw MaCh3Exception(__FILE__, __LINE__);
1181 }
1182
1183 //TODO: this assumes that std::vector is of length 2 and then both data.at(0) and data
1184 //ETA: I think this might just be wrong? We should probably just make this AddData(std::vector<double> data_x, std::vector<double> data_y)
1185 // or maybe something like AddData(std::vector<std::pair<double, double>> data)?
1186 for (int i = 0; i < int(data.size()); i++) {
1187 dathist2d->Fill(data.at(0)[i],data.at(1)[i]);
1188 }
1189
1190 if(SampleHandlerFD_data == nullptr) {
1191 MACH3LOG_ERROR("SampleHandlerFD_data haven't been initialised yet");
1192 throw MaCh3Exception(__FILE__, __LINE__);
1193 }
1194 for (size_t yBin=0;yBin<Binning.nYBins;yBin++) {
1195 for (size_t xBin=0;xBin<Binning.nXBins;xBin++) {
1196 //Need to cast to an int (Int_t) for ROOT
1197 //Need to do +1 for the bin, this is to be consistent with ROOTs binning scheme
1198 SampleHandlerFD_data[yBin][xBin] = dathist2d->GetBinContent(static_cast<int>(xBin+1),static_cast<int>(yBin+1));
1199 }
1200 }
1201}
1202
1204 MACH3LOG_INFO("Adding 1D data histogram: {} with {:.2f} events", Data->GetTitle(), Data->Integral());
1205 dathist2d = nullptr;
1206 dathist = static_cast<TH1D*>(Data->Clone());
1207
1208 if (GetNDim() != 1) {
1209 MACH3LOG_ERROR("Trying to set a 1D 'data' histogram in a 2D sample - Quitting");
1210 MACH3LOG_ERROR("The number of dimensions for this sample is {}", GetNDim());
1211 throw MaCh3Exception(__FILE__ , __LINE__ );
1212 }
1213
1214 if(SampleHandlerFD_data == nullptr) {
1215 MACH3LOG_ERROR("SampleHandlerFD_data haven't been initialised yet");
1216 throw MaCh3Exception(__FILE__, __LINE__);
1217 }
1218 for (size_t yBin=0;yBin<Binning.nYBins;yBin++) {
1219 for (size_t xBin=0;xBin<Binning.nXBins;xBin++) {
1220 //Need to cast to an int (Int_t) for ROOT
1221 //Need to do +1 for the bin, this is to be consistent with ROOTs binning scheme
1222 SampleHandlerFD_data[yBin][xBin] = Data->GetBinContent(static_cast<int>(xBin+1));
1223 }
1224 }
1225}
1226
1228 MACH3LOG_INFO("Adding 2D data histogram: {} with {:.2f} events", Data->GetTitle(), Data->Integral());
1229 dathist2d = static_cast<TH2D*>(Data->Clone());
1230 dathist = nullptr;
1231
1232 if (GetNDim() != 2) {
1233 MACH3LOG_ERROR("Trying to set a 2D 'data' histogram in a 1D sample - Quitting");
1234 throw MaCh3Exception(__FILE__ , __LINE__ );}
1235
1236 if(SampleHandlerFD_data == nullptr) {
1237 MACH3LOG_ERROR("SampleHandlerFD_data haven't been initialised yet");
1238 throw MaCh3Exception(__FILE__, __LINE__);
1239 }
1240 for (size_t yBin=0;yBin<Binning.nYBins;yBin++) {
1241 for (size_t xBin=0;xBin<Binning.nXBins;xBin++) {
1242 //Need to cast to an int (Int_t) for ROOT
1243 //Need to do +1 for the bin, this is to be consistent with ROOTs binning scheme
1244 SampleHandlerFD_data[yBin][xBin] = dathist2d->GetBinContent(static_cast<int>(xBin+1),static_cast<int>(yBin+1));
1245 }
1246 }
1247}
1248
1249// ************************************************
1251// ************************************************
1252 auto NuOscillatorConfigFile = Get<std::string>(SampleManager->raw()["NuOsc"]["NuOscConfigFile"], __FILE__ , __LINE__);
1253 auto EqualBinningPerOscChannel = Get<bool>(SampleManager->raw()["NuOsc"]["EqualBinningPerOscChannel"], __FILE__ , __LINE__);
1254
1255 // TN override the sample setting if not using binned oscillation
1256 if (EqualBinningPerOscChannel) {
1257 if (YAML::LoadFile(NuOscillatorConfigFile)["General"]["CalculationType"].as<std::string>() == "Unbinned") {
1258 MACH3LOG_WARN("Tried using EqualBinningPerOscChannel while using Unbinned oscillation calculation, changing EqualBinningPerOscChannel to false");
1259 EqualBinningPerOscChannel = false;
1260 }
1261 }
1262 std::vector<const double*> OscParams = ParHandler->GetOscParsFromSampleName(SampleName);
1263 if (OscParams.empty()) {
1264 MACH3LOG_ERROR("OscParams is empty for sample '{}'.", GetTitle());
1265 MACH3LOG_ERROR("This likely indicates an error in your oscillation YAML configuration.");
1266 throw MaCh3Exception(__FILE__, __LINE__);
1267 }
1268 Oscillator = std::make_shared<OscillationHandler>(NuOscillatorConfigFile, EqualBinningPerOscChannel, OscParams, static_cast<int>(OscChannels.size()));
1269
1270 if (!EqualBinningPerOscChannel) {
1271 for(int iChannel = 0; iChannel < static_cast<int>(OscChannels.size()); iChannel++) {
1272 std::vector<M3::float_t> EnergyArray;
1273 std::vector<M3::float_t> CosineZArray;
1274
1275 for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
1276 // KS: This is bit weird but we basically loop over all events and push to vector only these which are part of a given OscChannel
1277 const int Channel = GetOscChannel(OscChannels, (*MCSamples[iEvent].nupdgUnosc), (*MCSamples[iEvent].nupdg));
1278 //DB Remove NC events from the arrays which are handed to the NuOscillator objects
1279 if (!MCSamples[iEvent].isNC && Channel == iChannel) {
1280 EnergyArray.push_back(M3::float_t(*(MCSamples[iEvent].rw_etru)));
1281 }
1282 }
1283 std::sort(EnergyArray.begin(),EnergyArray.end());
1284
1285 //============================================================================
1286 //DB Atmospheric only part, can only happen if truecz has been initialised within the experiment specific code
1287 if (*(MCSamples[0].rw_truecz) != M3::_BAD_DOUBLE_) {
1288 for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
1289 // KS: This is bit weird but we basically loop over all events and push to vector only these which are part of a given OscChannel
1290 const int Channel = GetOscChannel(OscChannels, (*MCSamples[iEvent].nupdgUnosc), (*MCSamples[iEvent].nupdg));
1291 //DB Remove NC events from the arrays which are handed to the NuOscillator objects
1292 if (!MCSamples[iEvent].isNC && Channel == iChannel) {
1293 CosineZArray.push_back(M3::float_t(*(MCSamples[iEvent].rw_truecz)));
1294 }
1295 }
1296 std::sort(CosineZArray.begin(),CosineZArray.end());
1297 }
1298 Oscillator->SetOscillatorBinning(iChannel, EnergyArray, CosineZArray);
1299 } // end loop over channels
1300 }
1301}
1302
1304 for (unsigned int iEvent=0;iEvent<GetNEvents();iEvent++) {
1305 MCSamples[iEvent].osc_w_pointer = &M3::Unity;
1306 if (MCSamples[iEvent].isNC) {
1307 if (*MCSamples[iEvent].nupdg != *MCSamples[iEvent].nupdgUnosc) {
1308 MCSamples[iEvent].osc_w_pointer = &M3::Zero;
1309 } else {
1310 MCSamples[iEvent].osc_w_pointer = &M3::Unity;
1311 }
1312 } else {
1313 int InitFlav = M3::_BAD_INT_;
1314 int FinalFlav = M3::_BAD_INT_;
1315
1316 InitFlav = MaCh3Utils::PDGToNuOscillatorFlavour((*MCSamples[iEvent].nupdgUnosc));
1317 FinalFlav = MaCh3Utils::PDGToNuOscillatorFlavour((*MCSamples[iEvent].nupdg));
1318
1319 if (InitFlav == M3::_BAD_INT_ || FinalFlav == M3::_BAD_INT_) {
1320 MACH3LOG_ERROR("Something has gone wrong in the mapping between MCSamples.nutype and the enum used within NuOscillator");
1321 MACH3LOG_ERROR("MCSamples.nupdgUnosc: {}", (*MCSamples[iEvent].nupdgUnosc));
1322 MACH3LOG_ERROR("InitFlav: {}", InitFlav);
1323 MACH3LOG_ERROR("MCSamples.nupdg: {}", (*MCSamples[iEvent].nupdg));
1324 MACH3LOG_ERROR("FinalFlav: {}", FinalFlav);
1325 throw MaCh3Exception(__FILE__, __LINE__);
1326 }
1327
1328 const int OscIndex = GetOscChannel(OscChannels, (*MCSamples[iEvent].nupdgUnosc), (*MCSamples[iEvent].nupdg));
1329 //Can only happen if truecz has been initialised within the experiment specific code
1330 if (*(MCSamples[iEvent].rw_truecz) != M3::_BAD_DOUBLE_) {
1331 //Atmospherics
1332 MCSamples[iEvent].osc_w_pointer = Oscillator->GetNuOscillatorPointers(OscIndex, InitFlav, FinalFlav, FLOAT_T(*(MCSamples[iEvent].rw_etru)), FLOAT_T(*(MCSamples[iEvent].rw_truecz)));
1333 } else {
1334 //Beam
1335 MCSamples[iEvent].osc_w_pointer = Oscillator->GetNuOscillatorPointers(OscIndex, InitFlav, FinalFlav, FLOAT_T(*(MCSamples[iEvent].rw_etru)));
1336 }
1337 } // end if NC
1338 } // end loop over events
1339}
1340
1341std::string SampleHandlerFD::GetSampleName(int iSample) const {
1342 //ETA - this is just to suppress a warning for an unused variable
1343 (void)iSample;
1344
1345 //ETA - extra safety to make sure SampleName is actually set
1346 // probably unnecessary due to the requirement for it to be in the yaml config
1347 if(SampleName.length() == 0){
1348 MACH3LOG_ERROR("No sample name provided");
1349 MACH3LOG_ERROR("Please provide a SampleName in your configuration file: {}", SampleManager->GetFileName());
1350 throw MaCh3Exception(__FILE__, __LINE__);
1351 }
1352 return SampleName;
1353}
1354
1356 M3::float_t totalweight = 1.0;
1357 const int nParams = static_cast<int>(MCSamples[iEntry].total_weight_pointers.size());
1358 #ifdef MULTITHREAD
1359 #pragma omp simd
1360 #endif
1361 for (int iParam = 0; iParam < nParams; ++iParam) {
1362 totalweight *= *(MCSamples[iEntry].total_weight_pointers[iParam]);
1363 }
1364
1365 return totalweight;
1366}
1367
1372 //Now loop over events and get the spline bin for each event
1373 for (unsigned int j = 0; j < GetNEvents(); ++j) {
1374 const int OscIndex = GetOscChannel(OscChannels, (*MCSamples[j].nupdgUnosc), (*MCSamples[j].nupdg));
1375
1376 std::vector< std::vector<int> > EventSplines;
1377 switch(nDimensions){
1378 case 1:
1379 EventSplines = SplineHandler->GetEventSplines(GetSampleName(), OscIndex, int(*(MCSamples[j].mode)), *(MCSamples[j].rw_etru), *(MCSamples[j].x_var), 0.);
1380 break;
1381 case 2:
1382 EventSplines = SplineHandler->GetEventSplines(GetSampleName(), OscIndex, int(*(MCSamples[j].mode)), *(MCSamples[j].rw_etru), *(MCSamples[j].x_var), *(MCSamples[j].y_var));
1383 break;
1384 default:
1385 MACH3LOG_ERROR("Error in assigning spline bins because nDimensions = {}", nDimensions);
1386 MACH3LOG_ERROR("MaCh3 only supports splines binned in Etrue + the sample binning");
1387 MACH3LOG_ERROR("Please check the sample binning you specified in your sample config ");
1388 throw MaCh3Exception(__FILE__, __LINE__);
1389 break;
1390 }
1391 int NSplines = int(EventSplines.size());
1392 if(NSplines < 0){
1393 throw MaCh3Exception(__FILE__, __LINE__);
1394 }
1395 MCSamples[j].xsec_spline_pointers.resize(NSplines);
1396 for(size_t spline = 0; spline < MCSamples[j].xsec_spline_pointers.size(); spline++) {
1397 //Event Splines indexed as: sample name, oscillation channel, syst, mode, etrue, var1, var2 (var2 is a dummy 0 for 1D splines)
1398 MCSamples[j].xsec_spline_pointers[spline] = SplineHandler->retPointer(EventSplines[spline][0], EventSplines[spline][1],
1399 EventSplines[spline][2], EventSplines[spline][3],
1400 EventSplines[spline][4], EventSplines[spline][5],
1401 EventSplines[spline][6]);
1402 }
1403 }
1404}
1405
1406// ************************************************
1408// ************************************************
1409 if (SampleHandlerFD_data == nullptr) {
1410 MACH3LOG_ERROR("Data sample is empty! Can't calculate a likelihood!");
1411 throw MaCh3Exception(__FILE__, __LINE__);
1412 }
1413
1414 double negLogL = 0.;
1415 #ifdef MULTITHREAD
1416 #pragma omp parallel for collapse(2) reduction(+:negLogL)
1417 #endif
1418 for (size_t yBin = 0; yBin < Binning.nYBins; ++yBin)
1419 {
1420 for (size_t xBin = 0; xBin < Binning.nXBins; ++xBin)
1421 {
1422 const double DataVal = SampleHandlerFD_data[yBin][xBin];
1423 const double MCPred = SampleHandlerFD_array[yBin][xBin];
1424 const double w2 = SampleHandlerFD_array_w2[yBin][xBin];
1425
1426 //KS: Calculate likelihood using Barlow-Beeston Poisson or even IceCube
1427 negLogL += GetTestStatLLH(DataVal, MCPred, w2);
1428 }
1429 }
1430 return negLogL;
1431}
1432
1434 MCSamples.resize(nEvents);
1435}
1436
1437
1438// ************************************************
1440// ************************************************
1441 Dir->cd();
1442
1443 YAML::Node Config = SampleManager->raw();
1444 TMacro ConfigSave = YAMLtoTMacro(Config, (std::string("Config_") + GetTitle()));
1445 ConfigSave.Write();
1446
1447 std::unique_ptr<TH1> data_hist;
1448
1449 if (GetNDim() == 1) {
1450 data_hist = M3::Clone<TH1D>(dynamic_cast<TH1D*>(GetDataHist(1)), "data_" + GetTitle());
1451 data_hist->GetXaxis()->SetTitle(XVarStr.c_str());
1452 data_hist->GetYaxis()->SetTitle("Number of Events");
1453 } else if (GetNDim() == 2) {
1454 data_hist = M3::Clone<TH2D>(dynamic_cast<TH2D*>(GetDataHist(2)), "data_" + GetTitle());
1455 data_hist->GetXaxis()->SetTitle(XVarStr.c_str());
1456 data_hist->GetYaxis()->SetTitle(YVarStr.c_str());
1457 data_hist->GetZaxis()->SetTitle("Number of Events");
1458 } else {
1459 MACH3LOG_ERROR("Not implemented");
1460 throw MaCh3Exception(__FILE__, __LINE__);
1461 }
1462
1463 if (!data_hist) {
1464 MACH3LOG_ERROR("nullptr data hist :(");
1465 throw MaCh3Exception(__FILE__, __LINE__);
1466 }
1467
1468 data_hist->SetTitle(("data_" + GetTitle()).c_str());
1469 data_hist->Write();
1470}
1471
1473 std::vector<std::string> spline_filepaths;
1474 for(unsigned iChannel = 0 ; iChannel < OscChannels.size() ; iChannel++){
1475 spline_filepaths.push_back(spline_files[iChannel]);
1476 }
1477
1478 //Keep a track of the spline variables
1479 std::vector<std::string> SplineVarNames = {"TrueNeutrinoEnergy"};
1480 if(XVarStr.length() > 0){
1481 SplineVarNames.push_back(XVarStr);
1482 }
1483 if(YVarStr.length() > 0){
1484 SplineVarNames.push_back(YVarStr);
1485 }
1486
1487 SplineHandler->AddSample(SampleName, spline_filepaths, SplineVarNames);
1488 SplineHandler->CountNumberOfLoadedSplines(false, 1);
1489 SplineHandler->TransferToMonolith();
1490
1491 MACH3LOG_INFO("--------------------------------");
1492 MACH3LOG_INFO("Setup Far Detector splines");
1493
1495
1496 SplineHandler->cleanUpMemory();
1497}
1498
1499// === JM adjust GetNDVarHist functions to allow for subevent-level plotting ===
1500TH1* SampleHandlerFD::Get1DVarHist(const std::string& ProjectionVar_Str, const std::vector< KinematicCut >& EventSelectionVec,
1501 int WeightStyle, TAxis* Axis, const std::vector< KinematicCut >& SubEventSelectionVec) {
1502 //DB Need to overwrite the Selection member variable so that IsEventSelected function operates correctly.
1503 // Consequently, store the selection cuts already saved in the sample, overwrite the Selection variable, then reset
1504 std::vector< KinematicCut > tmp_Selection = Selection;
1505 std::vector< KinematicCut > SelectionVecToApply;
1506
1507 //DB Add all the predefined selections to the selection vector which will be applied
1508 for (size_t iSelec=0;iSelec<Selection.size();iSelec++) {
1509 SelectionVecToApply.emplace_back(Selection[iSelec]);
1510 }
1511
1512 //DB Add all requested cuts from the argument to the selection vector which will be applied
1513 for (size_t iSelec=0;iSelec<EventSelectionVec.size();iSelec++) {
1514 SelectionVecToApply.emplace_back(EventSelectionVec[iSelec]);
1515 }
1516
1517 //DB Set the member variable to be the cuts to apply
1518 Selection = SelectionVecToApply;
1519
1520 //DB Define the histogram which will be returned
1521 TH1D* _h1DVar;
1522 if (Axis) {
1523 _h1DVar = new TH1D("","",Axis->GetNbins(),Axis->GetXbins()->GetArray());
1524 } else {
1525 std::vector<double> xBinEdges = ReturnKinematicParameterBinning(ProjectionVar_Str);
1526 _h1DVar = new TH1D("", "", int(xBinEdges.size())-1, xBinEdges.data());
1527 }
1528
1529 if (IsSubEventVarString(ProjectionVar_Str)) {
1530 Fill1DSubEventHist(_h1DVar, ProjectionVar_Str, SubEventSelectionVec, WeightStyle);
1531 } else {
1532 //DB Grab the associated enum with the argument string
1533 int ProjectionVar_Int = ReturnKinematicParameterFromString(ProjectionVar_Str);
1534
1535 //DB Loop over all events
1536 for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
1537 if (IsEventSelected(iEvent)) {
1538 double Weight = GetEventWeight(iEvent);
1539 if (WeightStyle == 1) {
1540 Weight = 1.;
1541 }
1542 double Var = ReturnKinematicParameter(ProjectionVar_Int,iEvent);
1543 _h1DVar->Fill(Var,Weight);
1544 }
1545 }
1546 }
1547 //DB Reset the saved selection
1548 Selection = tmp_Selection;
1549
1550 return _h1DVar;
1551}
1552
1553void SampleHandlerFD::Fill1DSubEventHist(TH1D* _h1DVar, const std::string& ProjectionVar_Str, const std::vector< KinematicCut >& SubEventSelectionVec, int WeightStyle) {
1554 int ProjectionVar_Int = ReturnKinematicVectorFromString(ProjectionVar_Str);
1555
1556 //JM Loop over all events
1557 for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
1558 if (IsEventSelected(iEvent)) {
1559 double Weight = GetEventWeight(iEvent);
1560 if (WeightStyle == 1) {
1561 Weight = 1.;
1562 }
1563 std::vector<double> Vec = ReturnKinematicVector(ProjectionVar_Int,iEvent);
1564 size_t nsubevents = Vec.size();
1565 //JM Loop over all subevents in event
1566 for (unsigned int iSubEvent = 0; iSubEvent < nsubevents; iSubEvent++) {
1567 if (IsSubEventSelected(SubEventSelectionVec, iEvent, iSubEvent, nsubevents)) {
1568 double Var = Vec[iSubEvent];
1569 _h1DVar->Fill(Var,Weight);
1570 }
1571 }
1572 }
1573 }
1574}
1575
1576// ************************************************
1577TH2* SampleHandlerFD::Get2DVarHist(const std::string& ProjectionVar_StrX, const std::string& ProjectionVar_StrY,
1578 const std::vector< KinematicCut >& EventSelectionVec, int WeightStyle, TAxis* AxisX, TAxis* AxisY, const std::vector< KinematicCut >& SubEventSelectionVec) {
1579// ************************************************
1580 //DB Need to overwrite the Selection member variable so that IsEventSelected function operates correctly.
1581 // Consequently, store the selection cuts already saved in the sample, overwrite the Selection variable, then reset
1582 std::vector< KinematicCut > tmp_Selection = Selection;
1583 std::vector< KinematicCut > SelectionVecToApply;
1584
1585 //DB Add all the predefined selections to the selection vector which will be applied
1586 for (size_t iSelec=0;iSelec<Selection.size();iSelec++) {
1587 SelectionVecToApply.emplace_back(Selection[iSelec]);
1588 }
1589
1590 //DB Add all requested cuts from the argument to the selection vector which will be applied
1591 for (size_t iSelec=0;iSelec<EventSelectionVec.size();iSelec++) {
1592 SelectionVecToApply.emplace_back(EventSelectionVec[iSelec]);
1593 }
1594
1595 //DB Set the member variable to be the cuts to apply
1596 Selection = SelectionVecToApply;
1597
1598 //DB Define the histogram which will be returned
1599 TH2D* _h2DVar;
1600 if (AxisX && AxisY) {
1601 _h2DVar = new TH2D("","",AxisX->GetNbins(),AxisX->GetXbins()->GetArray(),AxisY->GetNbins(),AxisY->GetXbins()->GetArray());
1602 } else {
1603 std::vector<double> xBinEdges = ReturnKinematicParameterBinning(ProjectionVar_StrX);
1604 std::vector<double> yBinEdges = ReturnKinematicParameterBinning(ProjectionVar_StrY);
1605 _h2DVar = new TH2D("", "", int(xBinEdges.size())-1, xBinEdges.data(), int(yBinEdges.size())-1, yBinEdges.data());
1606 }
1607
1608 bool IsSubEventHist = IsSubEventVarString(ProjectionVar_StrX) || IsSubEventVarString(ProjectionVar_StrY);
1609 if (IsSubEventHist) Fill2DSubEventHist(_h2DVar, ProjectionVar_StrX, ProjectionVar_StrY, SubEventSelectionVec, WeightStyle);
1610 else {
1611 //DB Grab the associated enum with the argument string
1612 int ProjectionVar_IntX = ReturnKinematicParameterFromString(ProjectionVar_StrX);
1613 int ProjectionVar_IntY = ReturnKinematicParameterFromString(ProjectionVar_StrY);
1614
1615 //DB Loop over all events
1616 for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
1617 if (IsEventSelected(iEvent)) {
1618 double Weight = GetEventWeight(iEvent);
1619 if (WeightStyle == 1) {
1620 Weight = 1.;
1621 }
1622 double VarX = ReturnKinematicParameter(ProjectionVar_IntX, iEvent);
1623 double VarY = ReturnKinematicParameter(ProjectionVar_IntY, iEvent);
1624 _h2DVar->Fill(VarX,VarY,Weight);
1625 }
1626 }
1627 }
1628 //DB Reset the saved selection
1629 Selection = tmp_Selection;
1630
1631 return _h2DVar;
1632}
1633
1634void SampleHandlerFD::Fill2DSubEventHist(TH2D* _h2DVar, const std::string& ProjectionVar_StrX, const std::string& ProjectionVar_StrY,
1635 const std::vector< KinematicCut >& SubEventSelectionVec, int WeightStyle) {
1636 bool IsSubEventVarX = IsSubEventVarString(ProjectionVar_StrX);
1637 bool IsSubEventVarY = IsSubEventVarString(ProjectionVar_StrY);
1638
1639 int ProjectionVar_IntX, ProjectionVar_IntY;
1640 if (IsSubEventVarX) ProjectionVar_IntX = ReturnKinematicVectorFromString(ProjectionVar_StrX);
1641 else ProjectionVar_IntX = ReturnKinematicParameterFromString(ProjectionVar_StrX);
1642 if (IsSubEventVarY) ProjectionVar_IntY = ReturnKinematicVectorFromString(ProjectionVar_StrY);
1643 else ProjectionVar_IntY = ReturnKinematicParameterFromString(ProjectionVar_StrY);
1644
1645 //JM Loop over all events
1646 for (unsigned int iEvent = 0; iEvent < GetNEvents(); iEvent++) {
1647 if (IsEventSelected(iEvent)) {
1648 double Weight = GetEventWeight(iEvent);
1649 if (WeightStyle == 1) {
1650 Weight = 1.;
1651 }
1652 std::vector<double> VecX = {}, VecY = {};
1653 double VarX = M3::_BAD_DOUBLE_, VarY = M3::_BAD_DOUBLE_;
1654 size_t nsubevents = 0;
1655 // JM Three cases: subeventX vs eventY || eventX vs subeventY || subeventX vs subeventY
1656 if (IsSubEventVarX && !IsSubEventVarY) {
1657 VecX = ReturnKinematicVector(ProjectionVar_IntX, iEvent);
1658 VarY = ReturnKinematicParameter(ProjectionVar_IntY, iEvent);
1659 nsubevents = VecX.size();
1660 }
1661 else if (!IsSubEventVarX && IsSubEventVarY) {
1662 VecY = ReturnKinematicVector(ProjectionVar_IntY, iEvent);
1663 VarX = ReturnKinematicParameter(ProjectionVar_IntX, iEvent);
1664 nsubevents = VecY.size();
1665 }
1666 else {
1667 VecX = ReturnKinematicVector(ProjectionVar_IntX, iEvent);
1668 VecY = ReturnKinematicVector(ProjectionVar_IntY, iEvent);
1669 if (VecX.size() != VecY.size()) {
1670 MACH3LOG_ERROR("Cannot plot {} of size {} against {} of size {}", ProjectionVar_StrX, VecX.size(), ProjectionVar_StrY, VecY.size());
1671 throw MaCh3Exception(__FILE__, __LINE__);
1672 }
1673 nsubevents = VecX.size();
1674 }
1675 //JM Loop over all subevents in event
1676 for (unsigned int iSubEvent = 0; iSubEvent < nsubevents; iSubEvent++) {
1677 if (IsSubEventSelected(SubEventSelectionVec, iEvent, iSubEvent, nsubevents)) {
1678 if (IsSubEventVarX) VarX = VecX[iSubEvent];
1679 if (IsSubEventVarY) VarY = VecY[iSubEvent];
1680 _h2DVar->Fill(VarX,VarY,Weight);
1681 }
1682 }
1683 }
1684 }
1685}
1686// ================================================
1687
1688// ************************************************
1689int SampleHandlerFD::ReturnKinematicParameterFromString(const std::string& KinematicParameterStr) const {
1690// ************************************************
1691 auto it = KinematicParameters->find(KinematicParameterStr);
1692 if (it != KinematicParameters->end()) return it->second;
1693
1694 MACH3LOG_ERROR("Did not recognise Kinematic Parameter type: {}", KinematicParameterStr);
1695 throw MaCh3Exception(__FILE__, __LINE__);
1696
1697 return M3::_BAD_INT_;
1698}
1699
1700// ************************************************
1701std::string SampleHandlerFD::ReturnStringFromKinematicParameter(const int KinematicParameter) const {
1702// ************************************************
1703 auto it = ReversedKinematicParameters->find(KinematicParameter);
1704 if (it != ReversedKinematicParameters->end()) {
1705 return it->second;
1706 }
1707
1708 MACH3LOG_ERROR("Did not recognise Kinematic Parameter type: {}", KinematicParameter);
1709 throw MaCh3Exception(__FILE__, __LINE__);
1710
1711 return "";
1712}
1713
1714// === JM define KinematicVector-to-string mapping functions ===
1715// ************************************************
1716int SampleHandlerFD::ReturnKinematicVectorFromString(const std::string& KinematicVectorStr) const {
1717// ************************************************
1718 auto it = KinematicVectors->find(KinematicVectorStr);
1719 if (it != KinematicVectors->end()) return it->second;
1720
1721 MACH3LOG_ERROR("Did not recognise Kinematic Vector: {}", KinematicVectorStr);
1722 throw MaCh3Exception(__FILE__, __LINE__);
1723
1724 return M3::_BAD_INT_;
1725}
1726
1727// ************************************************
1728std::string SampleHandlerFD::ReturnStringFromKinematicVector(const int KinematicVector) const {
1729// ************************************************
1730 auto it = ReversedKinematicVectors->find(KinematicVector);
1731 if (it != ReversedKinematicVectors->end()) {
1732 return it->second;
1733 }
1734
1735 MACH3LOG_ERROR("Did not recognise Kinematic Vector: {}", KinematicVector);
1736 throw MaCh3Exception(__FILE__, __LINE__);
1737
1738 return "";
1739}
1740
1741// ************************************************
1742std::vector<double> SampleHandlerFD::ReturnKinematicParameterBinning(const std::string& KinematicParameter) {
1743// ************************************************
1744 // If x or y variable return used binning
1745 if(KinematicParameter == XVarStr) {
1746 return Binning.XBinEdges;
1747 } else if (KinematicParameter == YVarStr) {
1748 return Binning.YBinEdges;
1749 }
1750
1751 auto MakeBins = [](int nBins) {
1752 std::vector<double> bins(nBins + 1);
1753 for (int i = 0; i <= nBins; ++i)
1754 bins[i] = static_cast<double>(i) - 0.5;
1755 return bins;
1756 };
1757
1758 if (KinematicParameter == "OscillationChannel") {
1759 return MakeBins(GetNOscChannels());
1760 } else if (KinematicParameter == "Mode") {
1761 return MakeBins(Modes->GetNModes());
1762 }
1763
1764 // We first check if binning for a sample has been specified
1765 auto BinningConfig = M3OpenConfig(SampleManager->raw()["BinningFile"].as<std::string>());
1766 if(BinningConfig[GetTitle()] && BinningConfig[GetTitle()][KinematicParameter]){
1767 auto BinningVect = Get<std::vector<double>>(BinningConfig[GetTitle()][KinematicParameter], __FILE__, __LINE__);
1768 return BinningVect;
1769 } else {
1770 auto BinningVect = Get<std::vector<double>>(BinningConfig[KinematicParameter], __FILE__, __LINE__);
1771 return BinningVect;
1772 }
1773}
1774
1775
1776bool SampleHandlerFD::IsSubEventVarString(const std::string& VarStr) {
1777 if (KinematicVectors == nullptr) return false;
1778
1779 if (KinematicVectors->count(VarStr)) {
1780 if (!KinematicParameters->count(VarStr)) return true;
1781 else {
1782 MACH3LOG_ERROR("Attempted to plot kinematic variable {}, but it appears in both KinematicVectors and KinematicParameters", VarStr);
1783 throw MaCh3Exception(__FILE__,__LINE__);
1784 }
1785 }
1786 return false;
1787}
1788// ===============================================================
1789
1790TH1* SampleHandlerFD::Get1DVarHistByModeAndChannel(const std::string& ProjectionVar_Str,
1791 int kModeToFill, int kChannelToFill, int WeightStyle, TAxis* Axis) {
1792 bool fChannel;
1793 bool fMode;
1794
1795 if (kChannelToFill != -1) {
1796 if (kChannelToFill > GetNOscChannels()) {
1797 MACH3LOG_ERROR("Required channel is not available. kChannelToFill should be between 0 and {}", GetNOscChannels());
1798 MACH3LOG_ERROR("kChannelToFill given:{}", kChannelToFill);
1799 MACH3LOG_ERROR("Exiting.");
1800 throw MaCh3Exception(__FILE__, __LINE__);
1801 }
1802 fChannel = true;
1803 } else {
1804 fChannel = false;
1805 }
1806
1807 if (kModeToFill!=-1) {
1808 if (!(kModeToFill >= 0) && (kModeToFill < Modes->GetNModes())) {
1809 MACH3LOG_ERROR("Required mode is not available. kModeToFill should be between 0 and {}", Modes->GetNModes());
1810 MACH3LOG_ERROR("kModeToFill given:{}", kModeToFill);
1811 MACH3LOG_ERROR("Exiting..");
1812 throw MaCh3Exception(__FILE__, __LINE__);
1813 }
1814 fMode = true;
1815 } else {
1816 fMode = false;
1817 }
1818
1819 std::vector< KinematicCut > SelectionVec;
1820
1821 if (fMode) {
1822 KinematicCut SelecMode;
1824 SelecMode.LowerBound = kModeToFill;
1825 SelecMode.UpperBound = kModeToFill+1;
1826 SelectionVec.push_back(SelecMode);
1827 }
1828
1829 if (fChannel) {
1830 KinematicCut SelecChannel;
1831 SelecChannel.ParamToCutOnIt = ReturnKinematicParameterFromString("OscillationChannel");
1832 SelecChannel.LowerBound = kChannelToFill;
1833 SelecChannel.UpperBound = kChannelToFill+1;
1834 SelectionVec.push_back(SelecChannel);
1835 }
1836
1837 return Get1DVarHist(ProjectionVar_Str,SelectionVec,WeightStyle,Axis);
1838}
1839
1840TH2* SampleHandlerFD::Get2DVarHistByModeAndChannel(const std::string& ProjectionVar_StrX, const std::string& ProjectionVar_StrY,
1841 int kModeToFill, int kChannelToFill, int WeightStyle, TAxis* AxisX, TAxis* AxisY) {
1842 bool fChannel;
1843 bool fMode;
1844
1845 if (kChannelToFill!=-1) {
1846 if (kChannelToFill > GetNOscChannels()) {
1847 MACH3LOG_ERROR("Required channel is not available. kChannelToFill should be between 0 and {}", GetNOscChannels());
1848 MACH3LOG_ERROR("kChannelToFill given:{}", kChannelToFill);
1849 MACH3LOG_ERROR("Exiting.");
1850 throw MaCh3Exception(__FILE__, __LINE__);
1851 }
1852 fChannel = true;
1853 } else {
1854 fChannel = false;
1855 }
1856
1857 if (kModeToFill!=-1) {
1858 if (!(kModeToFill >= 0) && (kModeToFill < Modes->GetNModes())) {
1859 MACH3LOG_ERROR("Required mode is not available. kModeToFill should be between 0 and {}", Modes->GetNModes());
1860 MACH3LOG_ERROR("kModeToFill given:{}", kModeToFill);
1861 MACH3LOG_ERROR("Exiting..");
1862 throw MaCh3Exception(__FILE__, __LINE__);
1863 }
1864 fMode = true;
1865 } else {
1866 fMode = false;
1867 }
1868
1869 std::vector< KinematicCut > SelectionVec;
1870
1871 if (fMode) {
1872 KinematicCut SelecMode;
1874 SelecMode.LowerBound = kModeToFill;
1875 SelecMode.UpperBound = kModeToFill+1;
1876 SelectionVec.push_back(SelecMode);
1877 }
1878
1879 if (fChannel) {
1880 KinematicCut SelecChannel;
1881 SelecChannel.ParamToCutOnIt = ReturnKinematicParameterFromString("OscillationChannel");
1882 SelecChannel.LowerBound = kChannelToFill;
1883 SelecChannel.UpperBound = kChannelToFill+1;
1884 SelectionVec.push_back(SelecChannel);
1885 }
1886
1887 return Get2DVarHist(ProjectionVar_StrX,ProjectionVar_StrY,SelectionVec,WeightStyle,AxisX,AxisY);
1888}
1889
1890void SampleHandlerFD::PrintIntegral(TString OutputFileName, int WeightStyle, TString OutputCSVFileName) {
1891 int space = 14;
1892
1893 bool printToFile=false;
1894 if (OutputFileName.CompareTo("/dev/null")) {printToFile = true;}
1895
1896 bool printToCSV=false;
1897 if(OutputCSVFileName.CompareTo("/dev/null")) printToCSV=true;
1898
1899 std::ofstream outfile;
1900 if (printToFile) {
1901 outfile.open(OutputFileName.Data(), std::ios_base::app);
1902 outfile.precision(7);
1903 }
1904
1905 std::ofstream outcsv;
1906 if(printToCSV){
1907 outcsv.open(OutputCSVFileName, std::ios_base::app); // Appened to CSV
1908 outcsv.precision(7);
1909 }
1910
1911 double PDFIntegral = 0.;
1912
1913 std::vector< std::vector< TH1* > > IntegralList;
1914 IntegralList.resize(Modes->GetNModes());
1915
1916 std::vector<double> ChannelIntegral;
1917 ChannelIntegral.resize(GetNOscChannels());
1918 for (unsigned int i=0;i<ChannelIntegral.size();i++) {ChannelIntegral[i] = 0.;}
1919
1920 for (int i=0;i<Modes->GetNModes();i++) {
1921 if (GetNDim()==1) {
1922 IntegralList[i] = ReturnHistsBySelection1D(XVarStr,1,i,WeightStyle);
1923 } else {
1924 IntegralList[i] = CastVector<TH2, TH1>(ReturnHistsBySelection2D(XVarStr,YVarStr,1,i,WeightStyle));
1925 }
1926 }
1927
1928 MACH3LOG_INFO("-------------------------------------------------");
1929
1930 if (printToFile) {
1931 outfile << "\\begin{table}[ht]" << std::endl;
1932 outfile << "\\begin{center}" << std::endl;
1933 outfile << "\\caption{Integral breakdown for sample: " << GetTitle() << "}" << std::endl;
1934 outfile << "\\label{" << GetTitle() << "-EventRate}" << std::endl;
1935
1936 TString nColumns;
1937 for (int i=0;i<GetNOscChannels();i++) {nColumns+="|c";}
1938 nColumns += "|c|";
1939 outfile << "\\begin{tabular}{|l" << nColumns.Data() << "}" << std::endl;
1940 outfile << "\\hline" << std::endl;
1941 }
1942
1943 if(printToCSV){
1944 // HW Probably a better way but oh well, here I go making MaCh3 messy again
1945 outcsv<<"Integral Breakdown for sample :"<<GetTitle()<<"\n";
1946 }
1947
1948 MACH3LOG_INFO("Integral breakdown for sample: {}", GetTitle());
1949 MACH3LOG_INFO("");
1950
1951 if (printToFile) {outfile << std::setw(space) << "Mode:";}
1952 if(printToCSV) {outcsv<<"Mode,";}
1953
1954 std::string table_headings = fmt::format("| {:<8} |", "Mode");
1955 std::string table_footline = "------------"; //Scalable table horizontal line
1956 for (size_t i = 0;i < OscChannels.size(); i++) {
1957 table_headings += fmt::format(" {:<17} |", OscChannels[i].flavourName);
1958 table_footline += "--------------------";
1959 if (printToFile) {outfile << "&" << std::setw(space) << OscChannels[i].flavourName_Latex << " ";}
1960 if (printToCSV) {outcsv << OscChannels[i].flavourName << ",";}
1961 }
1962 if (printToFile) {outfile << "&" << std::setw(space) << "Total:" << "\\\\ \\hline" << std::endl;}
1963 if (printToCSV) {outcsv <<"Total\n";}
1964 table_headings += fmt::format(" {:<10} |", "Total");
1965 table_footline += "-------------";
1966
1967 MACH3LOG_INFO("{}", table_headings);
1968 MACH3LOG_INFO("{}", table_footline);
1969
1970 for (unsigned int i=0;i<IntegralList.size();i++) {
1971 double ModeIntegral = 0;
1972 if (printToFile) {outfile << std::setw(space) << Modes->GetMaCh3ModeName(i);}
1973 if(printToCSV) {outcsv << Modes->GetMaCh3ModeName(i) << ",";}
1974
1975 table_headings = fmt::format("| {:<8} |", Modes->GetMaCh3ModeName(i)); //Start string with mode name
1976
1977 for (unsigned int j=0;j<IntegralList[i].size();j++) {
1978 double Integral = IntegralList[i][j]->Integral();
1979
1980 if (Integral<1e-100) {Integral=0;}
1981
1982 ModeIntegral += Integral;
1983 ChannelIntegral[j] += Integral;
1984 PDFIntegral += Integral;
1985
1986 if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",Integral) << " ";}
1987 if (printToCSV) {outcsv << Form("%4.5f", Integral) << ",";}
1988
1989 table_headings += fmt::format(" {:<17.4f} |", Integral);
1990 }
1991 if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",ModeIntegral) << " \\\\ \\hline" << std::endl;}
1992 if (printToCSV) {outcsv << Form("%4.5f", ModeIntegral) << "\n";}
1993
1994 table_headings += fmt::format(" {:<10.4f} |", ModeIntegral);
1995
1996 MACH3LOG_INFO("{}", table_headings);
1997 }
1998
1999 if (printToFile) {outfile << std::setw(space) << "Total:";}
2000 if (printToCSV) {outcsv << "Total,";}
2001
2002 //Clear the table_headings to print last row of totals
2003 table_headings = fmt::format("| {:<8} |", "Total");
2004 for (unsigned int i=0;i<ChannelIntegral.size();i++) {
2005 if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",ChannelIntegral[i]) << " ";}
2006 if (printToCSV) {outcsv << Form("%4.5f", ChannelIntegral[i]) << ",";}
2007 table_headings += fmt::format(" {:<17.4f} |", ChannelIntegral[i]);
2008 }
2009 if (printToFile) {outfile << "&" << std::setw(space) << Form("%4.5f",PDFIntegral) << " \\\\ \\hline" << std::endl;}
2010 if (printToCSV) {outcsv << Form("%4.5f", PDFIntegral) << "\n\n\n\n";} // Let's have a few new lines!
2011
2012 table_headings += fmt::format(" {:<10.4f} |", PDFIntegral);
2013 MACH3LOG_INFO("{}", table_headings);
2014 MACH3LOG_INFO("{}", table_footline);
2015
2016 if (printToFile) {
2017 outfile << "\\end{tabular}" << std::endl;
2018 outfile << "\\end{center}" << std::endl;
2019 outfile << "\\end{table}" << std::endl;
2020 }
2021
2022 MACH3LOG_INFO("");
2023
2024 if (printToFile) {
2025 outfile << std::endl;
2026 outfile.close();
2027 }
2028}
2029
2030std::vector<TH1*> SampleHandlerFD::ReturnHistsBySelection1D(std::string KinematicProjection, int Selection1, int Selection2, int WeightStyle, TAxis* XAxis) {
2031 std::vector<TH1*> hHistList;
2032 std::string legendEntry;
2033
2034 if (THStackLeg != nullptr) {delete THStackLeg;}
2035 THStackLeg = new TLegend(0.1,0.1,0.9,0.9);
2036
2037 int iMax = -1;
2038 if (Selection1 == 0) {
2039 iMax = Modes->GetNModes();
2040 }
2041 if (Selection1 == 1) {
2042 iMax = GetNOscChannels();
2043 }
2044 if (iMax == -1) {
2045 MACH3LOG_ERROR("You've passed me a Selection1 which was not implemented in ReturnHistsBySelection1D. Selection1 and Selection2 are counters for different indexable quantities");
2046 throw MaCh3Exception(__FILE__, __LINE__);
2047 }
2048
2049 for (int i=0;i<iMax;i++) {
2050 if (Selection1==0) {
2051 hHistList.push_back(Get1DVarHistByModeAndChannel(KinematicProjection,i,Selection2,WeightStyle,XAxis));
2052 THStackLeg->AddEntry(hHistList[i],(Modes->GetMaCh3ModeName(i)+Form(" : (%4.2f)",hHistList[i]->Integral())).c_str(),"f");
2053
2054 hHistList[i]->SetFillColor(static_cast<Color_t>(Modes->GetMaCh3ModePlotColor(i)));
2055 hHistList[i]->SetLineColor(static_cast<Color_t>(Modes->GetMaCh3ModePlotColor(i)));
2056 }
2057 if (Selection1==1) {
2058 hHistList.push_back(Get1DVarHistByModeAndChannel(KinematicProjection,Selection2,i,WeightStyle,XAxis));
2059 THStackLeg->AddEntry(hHistList[i],(OscChannels[i].flavourName+Form(" | %4.2f",hHistList[i]->Integral())).c_str(),"f");
2060 }
2061 }
2062
2063 return hHistList;
2064}
2065
2066std::vector<TH2*> SampleHandlerFD::ReturnHistsBySelection2D(std::string KinematicProjectionX, std::string KinematicProjectionY,
2067 int Selection1, int Selection2, int WeightStyle,
2068 TAxis* XAxis, TAxis* YAxis) {
2069 std::vector<TH2*> hHistList;
2070
2071 int iMax = -1;
2072 if (Selection1 == 0) {
2073 iMax = Modes->GetNModes();
2074 }
2075 if (Selection1 == 1) {
2076 iMax = GetNOscChannels();
2077 }
2078 if (iMax == -1) {
2079 MACH3LOG_ERROR("You've passed me a Selection1 which was not implemented in ReturnHistsBySelection1D. Selection1 and Selection2 are counters for different indexable quantities");
2080 throw MaCh3Exception(__FILE__, __LINE__);
2081 }
2082
2083 for (int i=0;i<iMax;i++) {
2084 if (Selection1==0) {
2085 hHistList.push_back(Get2DVarHistByModeAndChannel(KinematicProjectionX,KinematicProjectionY,i,Selection2,WeightStyle,XAxis,YAxis));
2086 }
2087 if (Selection1==1) {
2088 hHistList.push_back(Get2DVarHistByModeAndChannel(KinematicProjectionX,KinematicProjectionY,Selection2,i,WeightStyle,XAxis,YAxis));
2089 }
2090 }
2091
2092 return hHistList;
2093}
2094
2095THStack* SampleHandlerFD::ReturnStackedHistBySelection1D(std::string KinematicProjection, int Selection1, int Selection2, int WeightStyle, TAxis* XAxis) {
2096 std::vector<TH1*> HistList = ReturnHistsBySelection1D(KinematicProjection, Selection1, Selection2, WeightStyle, XAxis);
2097 THStack* StackHist = new THStack((GetTitle()+"_"+KinematicProjection+"_Stack").c_str(),"");
2098 for (unsigned int i=0;i<HistList.size();i++) {
2099 StackHist->Add(HistList[i]);
2100 }
2101 return StackHist;
2102}
2103
2104
2105// ************************************************
2106const double* SampleHandlerFD::GetPointerToOscChannel(int iEvent) const {
2107// ************************************************
2108 const int Channel = GetOscChannel(OscChannels, (*MCSamples[iEvent].nupdgUnosc), (*MCSamples[iEvent].nupdg));
2109
2110 return &(OscChannels[Channel].ChannelIndex);
2111}
int GetOscChannel(const std::vector< OscChannelInfo > &OscChannel, const int InitFlav, const int FinalFlav)
KS: Get Osc Channel Index based on initial and final PDG codes.
int size
KS: Based on this https://github.com/gabime/spdlog/blob/a2b4262090fd3f005c2315dcb5be2f0f1774a005/incl...
#define MACH3LOG_DEBUG
Definition: MaCh3Logger.h:22
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
#define MACH3LOG_INFO
Definition: MaCh3Logger.h:23
#define MACH3LOG_WARN
Definition: MaCh3Logger.h:24
#define MACH3LOG_TRACE
Definition: MaCh3Logger.h:21
std::function< void(const double *, std::size_t)> FuncParFuncType
NuPDG
Enum to track the incoming neutrino species.
TestStatistic
Make an enum of the test statistic that we're using.
Type GetFromManager(const YAML::Node &node, Type defval, const std::string File="", const int Line=1)
Get content of config file if node is not found take default value specified.
Definition: YamlHelper.h:298
TMacro YAMLtoTMacro(const YAML::Node &yaml_node, const std::string &name)
Convert a YAML node to a ROOT TMacro object.
Definition: YamlHelper.h:161
bool CheckNodeExists(const YAML::Node &node, Args... args)
KS: Wrapper function to call the recursive helper.
Definition: YamlHelper.h:54
#define M3OpenConfig(filename)
Macro to simplify calling LoadYaml with file and line info.
Definition: YamlHelper.h:560
#define GetBounds(filename)
Definition: YamlHelper.h:561
Custom exception class for MaCh3 errors.
const double * RetPointer(const int iParam)
DB Pointer return to param position.
Class responsible for handling of systematic error parameters with different types defined in the con...
Class responsible for handling implementation of samples used in analysis, reweighting and returning ...
std::unique_ptr< MaCh3Modes > Modes
Holds information about used Generator and MaCh3 modes.
TestStatistic fTestStatistic
Test statistic tells what kind of likelihood sample is using.
M3::int_t nSamples
Contains how many samples we've got.
unsigned int nEvents
Number of MC events are there.
bool MatchCondition(const std::vector< T > &allowedValues, const T &value)
check if event is affected by following conditions, for example pdg, or modes etc
std::vector< std::vector< int > > funcParsGrid
HH - a grid of vectors of enums for each sample and event.
void CalcNormsBins(std::vector< NormParameter > &norm_parameters, std::vector< std::vector< int > > &xsec_norms_bins)
Check whether a normalisation systematic affects an event or not.
std::vector< FunctionalParameter > funcParsVec
HH - a vector that stores all the FuncPars struct.
virtual void resetShifts(int iEvent)
HH - reset the shifted values to the original values.
std::vector< FunctionalParameter * > funcParsMap
HH - a map that relates the funcpar enum to pointer of FuncPars struct HH - Changed to a vector of po...
virtual void Init()=0
Initialise any variables that your experiment specific SampleHandler needs.
std::string ReturnStringFromKinematicVector(const int KinematicVariable) const
std::vector< double > ReturnKinematicParameterBinning(const std::string &KinematicParameter)
Return the binning used to draw a kinematic parameter.
void AddData(TH1D *Data)
std::unique_ptr< BinnedSplineHandler > SplineHandler
Contains all your binned splines and handles the setup and the returning of weights from spline evalu...
virtual void CalcWeightFunc(int iEvent)
Calculate weights for function parameters.
SampleBinningInfo Binning
KS: This stores binning information, in future could be come vector to store binning for every used s...
virtual void SetupSplines()=0
initialise your splineXX object and then use InitialiseSplineObject to conviently setup everything up
bool IsEventSelected(const int iEvent)
DB Function which determines if an event is selected, where Selection double looks like {{ND280Kinema...
ParameterHandlerGeneric * ParHandler
std::string SampleTitle
the name of this sample e.g."muon-like"
bool FirstTimeW2
KS:Super hacky to update W2 or not.
std::unordered_map< std::string, NuPDG > FileToInitPDGMap
std::vector< KinematicCut > Selection
a way to store selection cuts which you may push back in the get1DVar functions most of the time this...
virtual int SetupExperimentMC()=0
Experiment specific setup, returns the number of events which were loaded.
const std::unordered_map< int, std::string > * ReversedKinematicVectors
double ** SampleHandlerFD_array_w2
KS Array used for MC stat.
std::vector< KinematicCut > StoredSelection
What gets pulled from config options, these are constant after loading in this is of length 3: 0th in...
virtual std::vector< double > ReturnKinematicVector(std::string KinematicParameter, int iEvent)
std::vector< std::string > spline_files
names of spline files associated associated with this object
double GetLikelihood() override
DB Multi-threaded GetLikelihood.
std::unordered_map< std::string, NuPDG > FileToFinalPDGMap
std::unordered_map< int, FuncParFuncType > funcParsFuncMap
HH - a map that relates the funcpar enum to pointer of the actual function.
const double * GetPointerToOscChannel(const int iEvent) const
Get pointer to oscillation channel associated with given event. Osc channel is const.
std::vector< OscChannelInfo > OscChannels
std::string XVarStr
the strings associated with the variables used for the binning e.g. "RecoNeutrinoEnergy"
const std::unordered_map< std::string, int > * KinematicVectors
void SetupNormParameters()
Setup the norm parameters by assigning each event with bin.
void SetupKinematicMap()
Ensure Kinematic Map is setup and make sure it is initialised correctly.
void FillSplineBins()
Finds the binned spline that an event should apply to and stored them in a a vector for easy evaluati...
virtual void PrepFunctionalParameters()
Update the functional parameter values to the latest proposed values. Needs to be called before every...
SampleHandlerFD(std::string ConfigFileName, ParameterHandlerGeneric *xsec_cov, const std::shared_ptr< OscillationHandler > &OscillatorObj_=nullptr)
Constructor.
void Fill2DSubEventHist(TH2D *_h2DVar, const std::string &ProjectionVarX, const std::string &ProjectionVarY, const std::vector< KinematicCut > &SubEventSelectionVec=std::vector< KinematicCut >(), int WeightStyle=0)
void Initialise()
Function which does a lot of the lifting regarding the workflow in creating different MC objects.
virtual ~SampleHandlerFD()
destructor
void FindNominalBinAndEdges1D()
std::vector< std::string > funcParsNamesVec
HH - a vector of string names for each functional parameter.
virtual void SetupFDMC()=0
Function which translates experiment struct into core struct.
std::string YVarStr
void SaveAdditionalInfo(TDirectory *Dir) override
Store additional info in a chan.
std::unique_ptr< manager > SampleManager
The manager object used to read the sample yaml file.
M3::float_t GetEventWeight(const int iEntry) const
void Fill1DHist()
Fill a 1D histogram with the event-level information used in the fit.
void InitialiseNuOscillatorObjects()
including Dan's magic NuOscillator
void Fill2DHist()
Fill a 2D histogram with the event-level information used in the fit.
void RegisterIndividualFunctionalParameter(const std::string &fpName, int fpEnum, FuncParFuncType fpFunc)
HH - a helper function for RegisterFunctionalParameter.
M3::float_t CalcWeightNorm(const FarDetectorCoreInfo *MCEvent) const
Calculate the norm weight for a given event.
bool UpdateW2
KS:Super hacky to update W2 or not.
std::shared_ptr< OscillationHandler > Oscillator
Contains all your binned splines and handles the setup and the returning of weights from spline evalu...
double ** SampleHandlerFD_array
DB Array to be filled after reweighting.
bool IsSubEventSelected(const std::vector< KinematicCut > &SubEventCuts, const int iEvent, unsigned const int iSubEvent, size_t nsubevents)
JM Function which determines if a subevent is selected.
const std::unordered_map< std::string, int > * KinematicParameters
Mapping between string and kinematic enum.
void Set2DBinning(size_t nbins1, double *boundaries1, size_t nbins2, double *boundaries2)
set the binning for 2D sample used for the likelihood calculation
void InitialiseSingleFDMCObject()
function to create the member of the FarDetectorInfo struct so they are the appropriate size.
TH2 * Get2DVarHist(const std::string &ProjectionVarX, const std::string &ProjectionVarY, const std::vector< KinematicCut > &EventSelectionVec=std::vector< KinematicCut >(), int WeightStyle=0, TAxis *AxisX=nullptr, TAxis *AxisY=nullptr, const std::vector< KinematicCut > &SubEventSelectionVec=std::vector< KinematicCut >())
virtual void RegisterFunctionalParameters()=0
HH - a experiment-specific function where the maps to actual functions are set up.
M3::float_t CalcWeightSpline(const FarDetectorCoreInfo *MCEvent) const
Calculate the spline weight for a given event.
std::vector< std::string > mc_files
names of mc files associated associated with this object
void ResetHistograms()
Helper function to reset histograms.
virtual double ReturnKinematicParameter(std::string KinematicParamter, int iEvent)=0
Return the value of an assocaited kinematic parameter for an event.
bool IsSubEventVarString(const std::string &VarStr)
JM Check if a kinematic parameter string corresponds to a subevent-level variable.
virtual void ApplyShifts(int iEvent)
ETA - generic function applying shifts.
std::unordered_map< std::string, double > _modeNomWeightMap
const std::unordered_map< int, std::string > * ReversedKinematicParameters
Mapping between kinematic enum and string.
int ReturnKinematicVectorFromString(const std::string &KinematicStr) const
double ** SampleHandlerFD_data
DB Array to be filled in AddData.
virtual void SetupWeightPointers()=0
DB Function to determine which weights apply to which types of samples pure virtual!...
void SetupSampleBinning()
wrapper to call set binning functions based on sample config info
void Reweight() override
void SetupNuOscillatorPointers()
virtual void SetupFunctionalParameters()
ETA - a function to setup and pass values to functional parameters where you need to pass a value to ...
void FillArray()
DB Nice new multi-threaded function which calculates the event weights and fills the relevant bins of...
int nDimensions
Keep track of the dimensions of the sample binning.
TLegend * THStackLeg
DB Miscellaneous Variables.
void PrintIntegral(TString OutputName="/dev/null", int WeightStyle=0, TString OutputCSVName="/dev/null")
std::vector< FarDetectorCoreInfo > MCSamples
void SetupReweightArrays(const size_t numberXBins, const size_t numberYBins)
Initialise data, MC and W2 histograms.
std::unordered_map< std::string, int > funcParsNamesMap
HH - a map that relates the name of the functional parameter to funcpar enum.
void Set1DBinning(size_t nbins, double *boundaries)
sets the binning used for the likelihood calculation, used for both data and MC
std::string SampleName
A unique ID for each sample based on powers of two for quick binary operator comparisons.
virtual const double * GetPointerToKinematicParameter(std::string KinematicParamter, int iEvent)=0
The manager class is responsible for managing configurations and settings.
Definition: Manager.h:16
const std::vector< NormParameter > GetNormParsFromSampleName(const std::string &SampleName) const
DB Get norm/func parameters depending on given SampleName.
std::string GetParFancyName(const int i) const
Get fancy name of the Parameter.
const std::vector< FunctionalParameter > GetFunctionalParametersFromSampleName(const std::string &SampleName) const
HH Get functional parameters for the relevant SampleName.
std::vector< const double * > GetOscParsFromSampleName(const std::string &SampleName)
Get pointers to Osc params from Sample name.
TH1 * Get1DVarHist(const std::string &ProjectionVar, const std::vector< KinematicCut > &EventSelectionVec=std::vector< KinematicCut >(), int WeightStyle=0, TAxis *Axis=nullptr, const std::vector< KinematicCut > &SubEventSelectionVec=std::vector< KinematicCut >())
std::string ReturnStringFromKinematicParameter(const int KinematicVariable) const
ETA function to generically convert a kinematic type from xsec cov to a string.
std::string GetSampleName(int iSample=0) const override
void Fill1DSubEventHist(TH1D *_h1DVar, const std::string &ProjectionVar, const std::vector< KinematicCut > &SubEventSelectionVec=std::vector< KinematicCut >(), int WeightStyle=0)
int GetNDim() const
DB Function to differentiate 1D or 2D binning.
unsigned int GetNEvents() const
std::vector< TH2 * > ReturnHistsBySelection2D(std::string KinematicProjectionX, std::string KinematicProjectionY, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *XAxis=0, TAxis *YAxis=0)
std::string GetTitle() const override
int GetNOscChannels() override
TH2 * Get2DVarHistByModeAndChannel(const std::string &ProjectionVar_StrX, const std::string &ProjectionVar_StrY, int kModeToFill=-1, int kChannelToFill=-1, int WeightStyle=0, TAxis *AxisX=nullptr, TAxis *AxisY=nullptr)
TH1 * GetW2Hist(const int Dimension)
Get W2 histogram.
TH1 * Get1DVarHistByModeAndChannel(const std::string &ProjectionVar_Str, int kModeToFill=-1, int kChannelToFill=-1, int WeightStyle=0, TAxis *Axis=nullptr)
TH1 * GetMCHist(const int Dimension)
Get MC histogram.
double GetTestStatLLH(double data, double mc) const
Calculate test statistic for a single bin using Poisson.
std::vector< TH1 * > ReturnHistsBySelection1D(std::string KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *Axis=0)
int ReturnKinematicParameterFromString(const std::string &KinematicStr) const
ETA function to generically convert a string from xsec cov to a kinematic type.
TH1 * GetDataHist(const int Dimension)
Get Data histogram.
THStack * ReturnStackedHistBySelection1D(std::string KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *Axis=0)
double float_t
Definition: Core.h:28
static constexpr const double _BAD_DOUBLE_
Default value used for double initialisation.
Definition: Core.h:43
static constexpr const double _DEFAULT_RETURN_VAL_
Definition: Core.h:46
static constexpr const float_t Unity
Definition: Core.h:54
static constexpr const float_t Zero
Definition: Core.h:61
int int_t
Definition: Core.h:29
static constexpr const int _BAD_INT_
Default value used for int initialisation.
Definition: Core.h:45
int PDGToNuOscillatorFlavour(int NuPdg)
Convert from PDG flavour to NuOscillator type beware that in the case of anti-neutrinos the NuOscilla...
constructors are same for all three so put in here
std::vector< const double * > xsec_norm_pointers
Pointers to normalisation weights which are being taken from Parameter Handler.
M3::float_t xsec_w
Total weight of norm and spline parameters.
std::vector< const M3::float_t * > xsec_spline_pointers
Pointers to spline weights which are being calculated by Splines Handler.
HH - Functional parameters Carrier for whether you want to apply a systematic to an event or not.
const double * valuePtr
Parameter value pointer.
FuncParFuncType * funcPtr
Function pointer.
KS: Small struct used for applying kinematic cuts.
double UpperBound
Upper bound on which we apply cut.
double LowerBound
Lower bound on which we apply cut.
int ParamToCutOnIt
Index or enum value identifying the kinematic variable to cut on.
KS: Store info about used osc channels.
std::string flavourName
Name of osc channel.
std::vector< double > YBinEdges
Vector to hold y-axis bin-edges.
std::vector< double > rw_upper_xbinedge
upper to check if Eb has moved the erec bin
size_t nXBins
Number of X axis bins in the histogram used for likelihood calculation.
std::vector< double > rw_lower_lower_xbinedge
lower to check if Eb has moved the erec bin
size_t nYBins
Number of Y axis bins in the histogram used for likelihood calculation.
std::vector< double > rw_lower_xbinedge
lower to check if Eb has moved the erec bin
std::vector< double > rw_upper_upper_xbinedge
upper to check if Eb has moved the erec bin
std::vector< double > XBinEdges
Vector to hold x-axis bin-edges.
int FindXBin(const double XVar, const int NomXBin) const
DB Find the relevant bin in the PDF for each event.