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(TestStatistic i) {
// **************************************************
  std::string name = "";
 
  switch(i) {
    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>(i));
      throw MaCh3Exception(__FILE__ , __LINE__ );
  }
  return name;
}
 
// ***************************
struct KinematicCut {
// ***************************
  int ParamToCutOnIt = M3::_BAD_INT_;
  double LowerBound = M3::_BAD_DOUBLE_;
  double UpperBound = M3::_BAD_DOUBLE_;
};
 
// ***************************
struct SampleBinningInfo {
// ***************************
  std::vector<double> XBinEdges;
  std::vector<double> YBinEdges;
 
  size_t nXBins = M3::_BAD_INT_;
  size_t nYBins = M3::_BAD_INT_;
 
  std::vector<double> rw_lower_xbinedge;
  std::vector<double> rw_lower_lower_xbinedge;
  std::vector<double> rw_upper_xbinedge;
  std::vector<double> rw_upper_upper_xbinedge;
 
  int FindXBin(const double XVar, const int NomXBin) const {
    //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 < rw_upper_xbinedge[NomXBin] && XVar >= rw_lower_xbinedge[NomXBin]) {
      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 < rw_lower_xbinedge[NomXBin] && XVar >= rw_lower_lower_xbinedge[NomXBin]) {
      return NomXBin-1;
    }
    //Shifted up one bin from the event bin at nominal
    else if (XVar < rw_upper_upper_xbinedge[NomXBin] && XVar >= rw_upper_xbinedge[NomXBin]) {
      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);
    }
  }
};
 
 
// ***************************
// 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(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