SampleStructs.h

Go to the documentation of this file.

#pragma once
 
// C++ includes
#include <set>
#include <list>
#include <unordered_map>
 
// MaCh3 includes
#include "Manager/MaCh3Exception.h"
#include "Manager/MaCh3Logger.h"
#include "Manager/Core.h"
 
_MaCh3_Safe_Include_Start_ //{
// ROOT include
#include "TSpline.h"
#include "TObjString.h"
#include "TFile.h"
#include "TF1.h"
#include "TH2Poly.h"
#include "TH1.h"
// NuOscillator includes
#include "Constants/OscillatorConstants.h"
_MaCh3_Safe_Include_End_ //}
 
 
 
// *******************
constexpr unsigned int str2int(const char* str, const int h = 0) {
// *******************
  return !str[h] ? 5381 : (str2int(str, h+1) * 33) ^ str[h];
}
 
 
// *****************
enum TargetMat {
// *****************
  kTarget_H  = 1,    
  kTarget_C  = 12,   
  kTarget_N  = 14,   
  kTarget_O  = 16,   
  kTarget_Al = 27,   
  kTarget_Ar = 40,   
  kTarget_Ti = 48,   
  kTarget_Fe = 56,   
  kTarget_Pb = 207   
};
 
// *****************
inline std::string TargetMat_ToString(const TargetMat i) {
// *****************
  std::string name;
 
  switch(i) {
    case kTarget_H:
      name = "Hydrogen";
      break;
    case kTarget_C:
      name = "Carbon";
      break;
    case kTarget_N:
      name = "Nitrogen";
      break;
    case kTarget_O:
      name = "Oxygen";
      break;
    case kTarget_Al:
      name = "Aluminium";
      break;
    case kTarget_Ar:
      name = "Argon";
      break;
    case kTarget_Ti:
      name = "Titanium";
      break;
    case kTarget_Fe:
      name = "Iron";
      break;
    case kTarget_Pb:
      name = "Lead";
      break;
    default:
      name = "TargetMat_Undefined";
      break;
  }
 
  return name;
}
 
// *****************
enum NuPDG {
// *****************
  kNue = 12,
  kNumu = 14,
  kNutau = 16,
  kNueBar = -12,
  kNumuBar = -14,
  kNutauBar = -16
};
 
enum TestStatistic {
  kPoisson,                 
  kBarlowBeeston,           
  kIceCube,                 
  kPearson,                 
  kDembinskiAbdelmotteleb,  
  kNTestStatistics          
};
 
// **************************************************
inline std::string TestStatistic_ToString(const TestStatistic TestStat) {
// **************************************************
  std::string name = "";
 
  switch(TestStat) {
    case TestStatistic::kPoisson:
    name = "Poisson";
    break;
    case TestStatistic::kBarlowBeeston:
    name = "Barlow-Beeston";
    break;
    case TestStatistic::kIceCube:
    name = "IceCube";
    break;
    case TestStatistic::kPearson:
    name = "Pearson";
    break;
    case TestStatistic::kDembinskiAbdelmotteleb:
    name = "Dembinski-Abdelmotteleb";
    break;
    case TestStatistic::kNTestStatistics:
      MACH3LOG_ERROR("kNTestStatistics is not a valid TestStatistic!");
      throw MaCh3Exception(__FILE__, __LINE__);
    default:
      MACH3LOG_ERROR("UNKNOWN LIKELIHOOD SPECIFIED!");
      MACH3LOG_ERROR("You gave test-statistic {}", static_cast<int>(TestStat));
      throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  return name;
}
 
// ***************************
struct KinematicCut {
// ***************************
  int ParamToCutOnIt = M3::_BAD_INT_;
  double LowerBound = M3::_BAD_DOUBLE_;
  double UpperBound = M3::_BAD_DOUBLE_;
};
 
 
// ***************************
struct BinShiftLookup {
// ***************************
  double lower_binedge;
  double lower_lower_binedge;
  double upper_binedge;
  double upper_upper_binedge;
};
 
// ***************************
struct SampleBinningInfo {
// ***************************
  std::vector<double> XBinEdges;
  std::vector<double> YBinEdges;
 
  size_t nXBins = M3::_BAD_INT_;
  size_t nYBins = M3::_BAD_INT_;
  size_t nBins = M3::_BAD_INT_;
  size_t GlobalOffset = M3::_BAD_INT_;
  std::vector<BinShiftLookup> xBinLookup;
 
  int GetBin(const int xBin, const int yBin) const {
    return static_cast<int>(yBin * nXBins + xBin);
  }
 
  int GetBinSafe(const int xBin, const int yBin) const {
    if (xBin < 0 || yBin < 0 || static_cast<size_t>(xBin) >= nXBins || static_cast<size_t>(yBin) >= nYBins) {
      MACH3LOG_ERROR("GetBinSafe: Bin indices out of range: xBin={}, yBin={}, max xBin={}, max yBin={}",
                     xBin, yBin, nXBins - 1, nYBins - 1);
      throw MaCh3Exception(__FILE__, __LINE__);
    }
    return GetBin(xBin, yBin);
  }
 
  int GetBinGlobal(const int xBin, const int yBin) const {
    return static_cast<int>(GlobalOffset + GetBin(xBin, yBin));
  }
 
  int FindXBin(const double XVar, const int NomXBin) const {
    // KS: Get reference to avoid repeated indexing and help with performance
    const auto& xBin = xBinLookup[NomXBin];
 
    //DB Check to see if momentum shift has moved bins
    //DB - First , check to see if the event is outside of the binning range and skip event if it is
     if (XVar < XBinEdges[0] || XVar >= XBinEdges[nXBins]) {
      return -1;
    }
    //DB - Second, check to see if the event is still in the nominal bin
    else if (XVar < xBin.upper_binedge && XVar >= xBin.lower_binedge) {
      return NomXBin;
    }
    //DB - Thirdly, check the adjacent bins first as Eb+CC+EScale shifts aren't likely to move an Erec more than 1bin width
    //Shifted down one bin from the event bin at nominal
    else if (XVar < xBin.lower_binedge && XVar >= xBin.lower_lower_binedge) {
      return NomXBin-1;
    }
    //Shifted up one bin from the event bin at nominal
    else if (XVar < xBin.upper_upper_binedge && XVar >= xBin.upper_binedge) {
      return NomXBin+1;
    }
    //DB - If we end up in this loop, the event has been shifted outside of its nominal bin, but is still within the allowed binning range
    else {
      // KS: Perform binary search to find correct bin. We already checked if isn't outside of bounds
      return static_cast<int>(std::distance(XBinEdges.begin(), std::upper_bound(XBinEdges.begin(), XBinEdges.end(), XVar)) - 1);
    }
  }
 
  void InitialiseLookUpSingleDimension(std::vector<BinShiftLookup>& BinLookup, const std::vector<double>& BinEdges, const size_t TotBins) {
    BinLookup.resize(TotBins);
    //Set rw_pdf_bin and upper_binedge and lower_binedge for each skmc_base
    for(size_t bin_i = 0; bin_i < TotBins; bin_i++){
      double low_lower_edge = M3::_DEFAULT_RETURN_VAL_;
      double low_edge = BinEdges[bin_i];
      double upper_edge = BinEdges[bin_i+1];
      double upper_upper_edge = M3::_DEFAULT_RETURN_VAL_;
 
      if (bin_i == 0) {
        low_lower_edge = BinEdges[0];
      } else {
        low_lower_edge = BinEdges[bin_i-1];
      }
 
      if (bin_i + 2 < TotBins) {
        upper_upper_edge = BinEdges[bin_i + 2];
      } else if (bin_i + 1 < TotBins) {
        upper_upper_edge = BinEdges[bin_i + 1];
      }
 
      BinLookup[bin_i].lower_binedge = low_edge;
      BinLookup[bin_i].upper_binedge = upper_edge;
      BinLookup[bin_i].lower_lower_binedge = low_lower_edge;
      BinLookup[bin_i].upper_upper_binedge = upper_upper_edge;
    }
  }
 
  void InitialiseBinMigrationLookUp() {
    InitialiseLookUpSingleDimension(xBinLookup, XBinEdges, nXBins);
  }
};
 
inline int GetSampleFromGlobalBin(const std::vector<SampleBinningInfo>& BinningInfo, const size_t GlobalBin) {
  for (size_t iSample = 0; iSample < BinningInfo.size(); ++iSample) {
    const SampleBinningInfo& info = BinningInfo[iSample];
    if (GlobalBin >= info.GlobalOffset && GlobalBin < info.GlobalOffset + info.nBins) {
      return static_cast<int>(iSample);
    }
  }
 
  MACH3LOG_ERROR("Couldn't find sample corresponding to bin {}", GlobalBin);
  throw MaCh3Exception(__FILE__, __LINE__);
 
  // GlobalBin is out of range for all samples
  return M3::_BAD_INT_;
}
 
inline void SetGlobalBinNumbers(std::vector<SampleBinningInfo>& BinningInfo) {
  if (BinningInfo.empty()) {
    MACH3LOG_ERROR("No binning samples provided.");
    throw MaCh3Exception(__FILE__, __LINE__);
  }
 
  size_t GlobalOffsetCounter = 0;
  for(size_t iSample = 0; iSample < BinningInfo.size(); iSample++){
    BinningInfo[iSample].GlobalOffset = GlobalOffsetCounter;
    GlobalOffsetCounter += BinningInfo[iSample].nBins;
  }
}
 
// ***************************
// A handy namespace for variables extraction
namespace MaCh3Utils {
// ***************************
  // *****************************
  inline double GetMassFromPDG(const int PDG) {
  // *****************************
    switch (abs(PDG)) {
      // Leptons
      case 11: return 0.00051099895; // e
      case 13: return 0.1056583755;  // mu
      case 15: return 1.77693;       // tau
      // Neutrinos
      case 12:
      case 14:
      case 16:
        return 0.; 
      // Photon
      case 22: return 0.; 
      // Mesons
      case 211: return 0.13957039; // pi_+/-
      case 111: return 0.1349768;  // pi_0
      case 221: return 0.547862;   // eta
      case 311:                    // K_0
      case 130:                    // K_0_L
      case 310:                    // K_0_S
        return 0.497611;
      case 321: return 0.493677;   // K_+/-
      // Baryons
      case 2112: return 0.939565; // n
      case 2212: return 0.938272; // p
      case 3122: return 1.115683; // lambda
      case 3222: return 1.118937; // sig_+
      case 3112: return 1.197449; // sig_-
      case 3212: return 1.192642; // sig_0
      // Nuclei
      case 1000050110: return 10.255103; // Boron-11
      case 1000060120: return 11.177929; // Carbon-12
      case 1000070140: return 13.043781; // Nitrogen-14
      case 1000080160: return 14.899169; // Oxygen-16
      case 1000090190: return 17.696901; // Fluorine-19
      case 1000110230: return 21.414835; // Sodium-23
      case 1000130270: return 25.133144; // Aluminum-27
      case 1000140280: return 26.060342; // Silicon-28
      case 1000190390: return 36.294463; // Potassium-39
      case 1000180400: return 37.224724; // Argon-40
      case 1000220480: return 44.663224; // Titanium-48
      case 1000300640: return 59.549619; // Zinc-64
      default:
        MACH3LOG_ERROR("Haven't got a saved mass for PDG: {}", PDG);
        MACH3LOG_ERROR("Please implement me!");
        throw MaCh3Exception(__FILE__, __LINE__);
    } // End switch
    MACH3LOG_ERROR("Warning, didn't catch a saved mass");
    return 0;
  }
  // ***************************
 
  // ***************************
  inline int PDGToNuOscillatorFlavour(const int NuPdg){
    int NuOscillatorFlavour = M3::_BAD_INT_;
    switch(std::abs(NuPdg)){
      case NuPDG::kNue:
        NuOscillatorFlavour = NuOscillator::kElectron;
        break;
      case NuPDG::kNumu:
        NuOscillatorFlavour = NuOscillator::kMuon;
        break;
      case NuPDG::kNutau:
        NuOscillatorFlavour = NuOscillator::kTau;
        break;
      default:
        MACH3LOG_ERROR("Unknown Neutrino PDG {}, cannot convert to NuOscillator type", NuPdg);
        break;
    }
 
    //This is very cheeky but if the PDG is negative then multiply the PDG by -1
    // This is consistent with the treatment that NuOscillator expects as enums only
    // exist for the generic matter flavour and not the anti-matter version
    if(NuPdg < 0){NuOscillatorFlavour *= -1;}
 
    return NuOscillatorFlavour;
  }
  // ***************************
 
  // ***************************
  inline std::string FormatDouble(const double value, const int precision) {
  // ***************************
    std::ostringstream oss;
    oss << std::fixed << std::setprecision(precision) << value;
    return oss.str();
  }
 
} // end MaCh3Utils namespace