MaCh3 2.2.1
Reference Guide
Loading...
Searching...
No Matches
Functions
Sample Handler Getters

Functions

virtual M3::int_t SampleHandlerBase::GetNsamples ()
 
virtual std::string SampleHandlerBase::GetTitle () const
 
virtual std::string SampleHandlerBase::GetSampleName (int Sample) const =0
 
virtual double SampleHandlerBase::GetSampleLikelihood (const int isample)
 
MaCh3ModesSampleHandlerBase::GetMaCh3Modes () const
 Return pointer to MaCh3 modes.
 
virtual double SampleHandlerBase::GetLikelihood ()=0
 
unsigned int SampleHandlerBase::GetNEvents () const
 
virtual int SampleHandlerBase::GetNMCSamples ()
 
virtual int SampleHandlerBase::GetNOscChannels ()
 
double SampleHandlerBase::GetTestStatLLH (double data, double mc) const
 Calculate test statistic for a single bin using Poisson.
 
double SampleHandlerBase::GetTestStatLLH (const double data, const double mc, const double w2) const
 Calculate test statistic for a single bin. Calculation depends on setting of fTestStatistic. Data and mc -> 0 cut-offs are defined in M3::LOW_MC_BOUND.
 
void SampleHandlerBase::SetTestStatistic (TestStatistic testStat)
 Set the test statistic to be used when calculating the binned likelihoods.
 
int SampleHandlerFD::GetNDim () const
 DB Function to differentiate 1D or 2D binning.
 
std::string SampleHandlerFD::GetSampleName (int iSample=0) const override
 
std::string SampleHandlerFD::GetTitle () const override
 
std::string SampleHandlerFD::GetXBinVarName ()
 
std::string SampleHandlerFD::GetYBinVarName ()
 
TH1 * SampleHandlerFD::GetMCHist (const int Dimension)
 Get MC histogram.
 
TH1 * SampleHandlerFD::GetW2Hist (const int Dimension)
 Get W2 histogram.
 
TH1 * SampleHandlerFD::GetDataHist (const int Dimension)
 Get Data histogram.
 
int SampleHandlerFD::GetNOscChannels () override
 
std::string SampleHandlerFD::GetFlavourName (const int iChannel)
 
TH1 * SampleHandlerFD::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 >())
 
void SampleHandlerFD::Fill1DSubEventHist (TH1D *_h1DVar, const std::string &ProjectionVar, const std::vector< KinematicCut > &SubEventSelectionVec=std::vector< KinematicCut >(), int WeightStyle=0)
 
TH1 * SampleHandlerFD::Get1DVarHistByModeAndChannel (const std::string &ProjectionVar_Str, int kModeToFill=-1, int kChannelToFill=-1, int WeightStyle=0, TAxis *Axis=nullptr)
 
TH2 * SampleHandlerFD::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 * SampleHandlerFD::GetModeHist1D (int s, int m, int style=0)
 
TH2 * SampleHandlerFD::GetModeHist2D (int s, int m, int style=0)
 
std::vector< TH1 * > SampleHandlerFD::ReturnHistsBySelection1D (std::string KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *Axis=0)
 
std::vector< TH2 * > SampleHandlerFD::ReturnHistsBySelection2D (std::string KinematicProjectionX, std::string KinematicProjectionY, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *XAxis=0, TAxis *YAxis=0)
 
THStack * SampleHandlerFD::ReturnStackedHistBySelection1D (std::string KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *Axis=0)
 
TLegend * SampleHandlerFD::ReturnStackHistLegend ()
 
int SampleHandlerFD::ReturnKinematicParameterFromString (const std::string &KinematicStr) const
 ETA function to generically convert a string from xsec cov to a kinematic type.
 
std::string SampleHandlerFD::ReturnStringFromKinematicParameter (const int KinematicVariable) const
 ETA function to generically convert a kinematic type from xsec cov to a string.
 

Detailed Description

Group of functions to get various parameters, names, and values.

Function Documentation

◆ Fill1DSubEventHist()

void SampleHandlerFD::Fill1DSubEventHist ( TH1D *  _h1DVar,
const std::string &  ProjectionVar,
const std::vector< KinematicCut > &  SubEventSelectionVec = std::vector< KinematicCut >(),
int  WeightStyle = 0 
)

Definition at line 1553 of file SampleHandlerFD.cpp.

1553 {
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}
bool IsEventSelected(const int iEvent)
DB Function which determines if an event is selected, where Selection double looks like {{ND280Kinema...
virtual std::vector< double > ReturnKinematicVector(std::string KinematicParameter, int iEvent)
M3::float_t GetEventWeight(const int iEntry) const
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.
int ReturnKinematicVectorFromString(const std::string &KinematicStr) const
unsigned int GetNEvents() const

◆ Get1DVarHist()

TH1 * SampleHandlerFD::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 >() 
)

Definition at line 1500 of file SampleHandlerFD.cpp.

1501 {
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}
std::vector< double > ReturnKinematicParameterBinning(const std::string &KinematicParameter)
Return the binning used to draw a kinematic parameter.
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 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.
void Fill1DSubEventHist(TH1D *_h1DVar, const std::string &ProjectionVar, const std::vector< KinematicCut > &SubEventSelectionVec=std::vector< KinematicCut >(), int WeightStyle=0)
int ReturnKinematicParameterFromString(const std::string &KinematicStr) const
ETA function to generically convert a string from xsec cov to a kinematic type.

◆ Get1DVarHistByModeAndChannel()

TH1 * SampleHandlerFD::Get1DVarHistByModeAndChannel ( const std::string &  ProjectionVar_Str,
int  kModeToFill = -1,
int  kChannelToFill = -1,
int  WeightStyle = 0,
TAxis *  Axis = nullptr 
)

Definition at line 1790 of file SampleHandlerFD.cpp.

1791 {
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}
#define MACH3LOG_ERROR
Definition: MaCh3Logger.h:25
Custom exception class for MaCh3 errors.
std::unique_ptr< MaCh3Modes > Modes
Holds information about used Generator and MaCh3 modes.
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 >())
int GetNOscChannels() override
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.

◆ Get2DVarHistByModeAndChannel()

TH2 * SampleHandlerFD::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 
)

Definition at line 1840 of file SampleHandlerFD.cpp.

1841 {
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}
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 >())

◆ GetDataHist()

TH1 * SampleHandlerFD::GetDataHist ( const int  Dimension)

Get Data histogram.

Definition at line 1136 of file SampleHandlerFD.cpp.

1136 {
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}

◆ GetFlavourName()

std::string SampleHandlerFD::GetFlavourName ( const int  iChannel)
inline

Definition at line 83 of file SampleHandlerFD.h.

83 {
84 if (iChannel < 0 || iChannel > static_cast<int>(OscChannels.size())) {
85 MACH3LOG_ERROR("Invalid Channel Requested: {}", iChannel);
86 throw MaCh3Exception(__FILE__ , __LINE__);
87 }
88 return OscChannels[iChannel].flavourName;
89 }
std::vector< OscChannelInfo > OscChannels

◆ GetLikelihood()

virtual double SampleHandlerBase::GetLikelihood ( )
pure virtual

Implemented in PySampleHandlerBase, and SampleHandlerFD.

◆ GetMaCh3Modes()

MaCh3Modes * SampleHandlerBase::GetMaCh3Modes ( ) const
inline

Return pointer to MaCh3 modes.

Definition at line 53 of file SampleHandlerBase.h.

53{ return Modes.get(); }

◆ GetMCHist()

TH1 * SampleHandlerFD::GetMCHist ( const int  Dimension)

Get MC histogram.

Definition at line 1121 of file SampleHandlerFD.cpp.

1121 {
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}
void Fill1DHist()
Fill a 1D histogram with the event-level information used in the fit.
void Fill2DHist()
Fill a 2D histogram with the event-level information used in the fit.

◆ GetModeHist1D()

TH1 * SampleHandlerFD::GetModeHist1D ( int  s,
int  m,
int  style = 0 
)
inline

Definition at line 111 of file SampleHandlerFD.h.

111 {
112 return Get1DVarHistByModeAndChannel(XVarStr,m,s,style);
113 }
std::string XVarStr
the strings associated with the variables used for the binning e.g. "RecoNeutrinoEnergy"
TH1 * Get1DVarHistByModeAndChannel(const std::string &ProjectionVar_Str, int kModeToFill=-1, int kChannelToFill=-1, int WeightStyle=0, TAxis *Axis=nullptr)

◆ GetModeHist2D()

TH2 * SampleHandlerFD::GetModeHist2D ( int  s,
int  m,
int  style = 0 
)
inline

Definition at line 115 of file SampleHandlerFD.h.

115 {
117 }
std::string YVarStr
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)

◆ GetNDim()

int SampleHandlerFD::GetNDim ( ) const
inline

DB Function to differentiate 1D or 2D binning.

Definition at line 31 of file SampleHandlerFD.h.

31{ return nDimensions; }
int nDimensions
Keep track of the dimensions of the sample binning.

◆ GetNEvents()

unsigned int SampleHandlerBase::GetNEvents ( ) const
inline

Definition at line 60 of file SampleHandlerBase.h.

60{return nEvents;}
unsigned int nEvents
Number of MC events are there.

◆ GetNMCSamples()

virtual int SampleHandlerBase::GetNMCSamples ( )
inlinevirtual

Definition at line 62 of file SampleHandlerBase.h.

62{ return nSamples; }
M3::int_t nSamples
Contains how many samples we've got.

◆ GetNOscChannels() [1/2]

virtual int SampleHandlerBase::GetNOscChannels ( )
inlinevirtual

Reimplemented in SampleHandlerFD.

Definition at line 64 of file SampleHandlerBase.h.

64{ return 1; }

◆ GetNOscChannels() [2/2]

int SampleHandlerFD::GetNOscChannels ( )
inlineoverridevirtual

Reimplemented from SampleHandlerBase.

Definition at line 80 of file SampleHandlerFD.h.

80{return static_cast<int>(OscChannels.size());}

◆ GetNsamples()

virtual M3::int_t SampleHandlerBase::GetNsamples ( )
inlinevirtual

Definition at line 40 of file SampleHandlerBase.h.

40{ return nSamples; };

◆ GetSampleLikelihood()

virtual double SampleHandlerBase::GetSampleLikelihood ( const int  isample)
inlinevirtual

Definition at line 46 of file SampleHandlerBase.h.

46{(void) isample; return GetLikelihood();};
virtual double GetLikelihood()=0

◆ GetSampleName() [1/2]

std::string SampleHandlerFD::GetSampleName ( int  iSample = 0) const
overridevirtual

Implements SampleHandlerBase.

Definition at line 1341 of file SampleHandlerFD.cpp.

1341 {
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}
std::unique_ptr< manager > SampleManager
The manager object used to read the sample yaml file.
std::string SampleName
A unique ID for each sample based on powers of two for quick binary operator comparisons.

◆ GetSampleName() [2/2]

virtual std::string SampleHandlerBase::GetSampleName ( int  Sample) const
pure virtual

Implemented in PySampleHandlerBase, and SampleHandlerFD.

◆ GetTestStatLLH() [1/2]

double SampleHandlerBase::GetTestStatLLH ( const double  data,
const double  mc,
const double  w2 
) const

Calculate test statistic for a single bin. Calculation depends on setting of fTestStatistic. Data and mc -> 0 cut-offs are defined in M3::LOW_MC_BOUND.

Implemented fTestStatistic are kPoisson (with Stirling's approx.), kBarlowBeeston (arXiv:1103.0354), kDembinskiAbdelmotteleb (arXiv:2206.12346), kIceCube (arxiv:1901.04645), and kPearson. Test statistics require mc > 0, therefore low mc and data values are treated with cut-offs based on M3::LOW_MC_BOUND = .00001 by default. For kPoisson, kBarlowBeeston, kDembinskiAbdelmotteleb, kPearson: data > LOW_MC_BOUND & mc <= LOW_MC_BOUND: returns GetTestStatLLH(data, LOW_MC_BOUND, w2), with Poisson(data,LOW_MC_BOUND) limit for mc->0, w2->0. mc < data <= LOW_MC_BOUND: returns 0 (as if any data <= LOW_MC_BOUND were effectively consistent with 0 data count), with a limit of 0 for mc->0. data = 0: returns mc (or mc/2. for kPearson), with a limit of 0 for mc->0. For kIceCube: mc < data returns the lower of IceCube(data,mc,w2) and Poisson(data,mc) penalties, with a Poisson(data,LOW_MC_BOUND) limit for mc->0, w2->0.

Parameters
datais data
mcis mc
w2is is Sum(w_{i}^2) (sum of weights squared), which is sigma^2_{MC stats}

Definition at line 36 of file SampleHandlerBase.cpp.

36 {
37// *************************
38 switch (fTestStatistic)
39 {
40 //CW: Not full Barlow-Beeston or what is referred to as "light": we're not introducing any more parameters
41 // Assume the MC has a Gaussian distribution around generated
42 // As in https://arxiv.org/abs/1103.0354 eq 10, 11
43 //CW: Calculate the Barlow-Beeston likelihood contribution from MC statistics
44 // Assumes the beta scaling parameters are Gaussian distributed
45 // Follows arXiv:1103.0354 section 5 and equation 8, 9, 10, 11 on page 4/5
46 // Essentially solves equation 11
47 case (kBarlowBeeston):
48 {
49 // The MC used in the likelihood calculation is allowed to be changed by Barlow Beeston beta parameters
50 double newmc = mc;
51
52 // If MC falls below the low MC bound, use low MC bound for newmc
53 if ( mc < M3::_LOW_MC_BOUND_ ) {
54 if ( data > M3::_LOW_MC_BOUND_ ) newmc = M3::_LOW_MC_BOUND_;
55 else if ( data >= mc ) return 0.;
56 }
57
58 // Barlow-Beeston uses fractional uncertainty on MC, so sqrt(sum[w^2])/mc
59 const double fractional = std::sqrt( w2 ) / newmc;
60 // fractional^2 to avoid doing same operation several times
61 const double fractional2 = fractional * fractional;
62 // b in quadratic equation
63 const double temp = newmc * fractional2 - 1;
64 // b^2 - 4ac in quadratic equation
65 const double temp2 = temp * temp + 4 * data * fractional2;
66 if ( temp2 < 0 ) {
67 MACH3LOG_ERROR("Negative square root in Barlow Beeston coefficient calculation!");
68 throw MaCh3Exception(__FILE__ , __LINE__ );
69 }
70 // Solve for the positive beta
71 const double beta = ( -1 * temp + sqrt( temp2 ) ) / 2.;
72
73 // If there is no data, test-stat shall return only MC*beta
74 double stat = mc * beta;
75 // With data, test-stat shall return LogL for newMC*beta which includes the low MC bound
76 if ( data > 0 ) {
77 newmc *= beta;
78 stat = newmc - data + data * std::log( data / newmc );
79 }
80
81 // Now, MC stat penalty
82 // The penalty from MC statistics using Conways approach (https://cds.cern.ch/record/1333496?)
83 double penalty = 0;
84 if ( fractional > 0 ) penalty = ( beta - 1 ) * ( beta - 1 ) / ( 2 * fractional2 );
85
86 // Returns test-stat plus the MC stat penalty
87 return stat+penalty;
88 }
89 break;
90 //KS: Alternative calculation of Barlow-Beeston following Hans Dembinski and Ahmed Abdelmottele arXiv:2206.12346v2
92 {
93 //KS: code follows authors implementation from:
94 //https://github.com/scikit-hep/iminuit/blob/059d06b00cae097ebf340b218b4eb57357111df8/src/iminuit/cost.py#L274-L300
95
96 // If there is no MC stat error for any reason, return Poisson LogL
97 if ( w2 == 0 ) return GetTestStatLLH(data,mc);
98
99 // The MC can be changed
100 double newmc = mc;
101
102 // If MC falls below the low MC bound, use low MC bound for newmc
103 if ( mc < M3::_LOW_MC_BOUND_ ) {
104 if ( data > M3::_LOW_MC_BOUND_ ) newmc = M3::_LOW_MC_BOUND_;
105 else if ( data >= mc ) return 0.;
106 }
107
108 //the so-called effective count
109 const double k = newmc * newmc / w2;
110 //Calculate beta which is scaling factor between true and generated MC
111 const double beta = ( data + k ) / ( newmc + k );
112
113 newmc *= beta;
114
115 // And penalise the movement in beta relative the mc uncertainty
116 const double penalty = k * beta - k + k * std::log( k / ( k * beta ) );
117
118 // If there are no data, this shall return newmc
119 double stat = newmc;
120 // All likelihood calculations may use the bare Poisson likelihood, so calculate here
121 // Only if there are some data
122 if ( data > 0 ) stat = newmc - data + data * std::log( data / newmc );
123
124 // Return the statistical contribution and penalty
125 return stat+penalty;
126 }
127 break;
128 //CW: Also try the IceCube likelihood
129 // It does not modify the MC content
130 // https://arxiv.org/abs/1901.04645
131 // Argüelles, C.A., Schneider, A. & Yuan, T. J. High Energ. Phys. (2019) 2019: 30. https://doi.org/10.1007/JHEP06(2019)030
132 // We essentially construct eq 3.16 and take the logarithm
133 // in eq 3.16, mu is MC, sigma2 is w2, k is data
134 case (kIceCube):
135 {
136 // IceCube low MC bound is implemented to return Poisson(data, _LOW_MC_BOUND_)
137 // up until the IceCube(data, mc) test-statistic is less than Poisson(data, _LOW_MC_BOUND_)
138 // The 0 MC limit is set to Poisson(data, 0.) as there is no way to get a non-diverging and reasonable guess on w2
139
140 // If there is 0 MC uncertainty (technically also when MC is 0) => Return Poisson(data, mc)
141 if ( w2 == 0 ) return GetTestStatLLH(data,mc);
142
143 // Auxiliary variables
144 const long double b = mc / w2;
145 const long double a = mc * b + 1;
146
147 // Use C99's implementation of log of gamma function to not be C++11 dependent
148 const double stat = double( -1 * ( a * logl( b ) + lgammal( data + a ) - lgammal( data + 1 ) - ( ( data + a ) * log1pl( b ) ) - lgammal( a ) ) );
149
150 // Check whether the stat is more than Poisson-like bound for low mc (mc < data)
151 // TN: I believe this might get some extra optimization
152 if ( mc <= data ) {
153 if ( data <= M3::_LOW_MC_BOUND_ ) return 0.;
154 const double poisson = GetTestStatLLH(data, M3::_LOW_MC_BOUND_);
155 if ( stat > poisson ) return poisson;
156 }
157
158 // Otherwise, return IceCube test-stat
159 return stat;
160 }
161 break;
162 //KS: Pearson works on assumption that event distribution in each bin is described by a Gaussian which in our case is not fulfilled for all bins, hence use it at your own risk
163 case (kPearson):
164 {
165 //KS: 2 is because this function returns -LLH not -2LLH
166 // With no data return the MC/2.
167 if ( data == 0 ) return mc/2.;
168
169 // If MC is lower than the low MC bound, return the test-stat at the bound
170 if ( mc < M3::_LOW_MC_BOUND_ ) {
171 if ( data > M3::_LOW_MC_BOUND_ ) return ( data - M3::_LOW_MC_BOUND_ ) * ( data - M3::_LOW_MC_BOUND_ ) / ( 2. * M3::_LOW_MC_BOUND_ );
172 else if ( data >= mc ) return 0.;
173 }
174
175 // Return the Pearson metric
176 return ( data - mc ) * ( data - mc ) / ( 2 * mc );
177 }
178 break;
179 case (kPoisson):
180 {
181 //Just call GetTestStatLLH which doesn't take in weights
182 //and is a Poisson likelihood comparison.
183 return GetTestStatLLH(data, mc);
184 }
185 break;
186 case TestStatistic::kNTestStatistics:
187 MACH3LOG_ERROR("kNTestStatistics is not a valid TestStatistic!");
188 throw MaCh3Exception(__FILE__, __LINE__);
189 default:
190 MACH3LOG_ERROR("Couldn't find TestStatistic {} exiting!", static_cast<int>(fTestStatistic));
191 throw MaCh3Exception(__FILE__ , __LINE__ );
192 } // end switch
193}
@ kPearson
Standard Pearson likelihood .
@ kBarlowBeeston
Barlow-Beeston () following Conway approximation ()
@ kIceCube
Based on .
@ kDembinskiAbdelmotteleb
Based on .
@ kPoisson
Standard Poisson likelihood .
TestStatistic fTestStatistic
Test statistic tells what kind of likelihood sample is using.
double GetTestStatLLH(double data, double mc) const
Calculate test statistic for a single bin using Poisson.
static constexpr const double _LOW_MC_BOUND_
MC prediction lower bound in bin to identify problematic binning definitions and handle LogL calculat...
Definition: Core.h:73

◆ GetTestStatLLH() [2/2]

double SampleHandlerBase::GetTestStatLLH ( double  data,
double  mc 
) const

Calculate test statistic for a single bin using Poisson.

Parameters
datais data
mcis mc

Definition at line 18 of file SampleHandlerBase.cpp.

18 {
19// ***************************************************************************
20 // Return MC if there are no data, returns 0 for data == 0 && mc == 0
21 if ( data == 0 ) return mc;
22
23 // If there are some data, but the prediction falls below the MC bound => return Poisson LogL for the low MC bound
24 if ( mc < M3::_LOW_MC_BOUND_ ) {
25 if ( data > M3::_LOW_MC_BOUND_ ) return ( M3::_LOW_MC_BOUND_ - data + data * std::log( data/M3::_LOW_MC_BOUND_ ) );
26 else if ( data >= mc ) return 0.;
27 }
28
29 // Otherwise, just return usual Poisson LogL using Stirling's approximation
30 // http://hyperphysics.phy-astr.gsu.edu/hbase/math/stirling.html
31 return ( mc - data + data * std::log( data / mc ) );
32}

◆ GetTitle() [1/2]

virtual std::string SampleHandlerBase::GetTitle ( ) const
inlinevirtual

Reimplemented in SampleHandlerFD.

Definition at line 42 of file SampleHandlerBase.h.

42{return "SampleHandler";};

◆ GetTitle() [2/2]

std::string SampleHandlerFD::GetTitle ( ) const
inlineoverridevirtual

Reimplemented from SampleHandlerBase.

Definition at line 35 of file SampleHandlerFD.h.

35{return SampleTitle;}
std::string SampleTitle
the name of this sample e.g."muon-like"

◆ GetW2Hist()

TH1 * SampleHandlerFD::GetW2Hist ( const int  Dimension)

Get W2 histogram.

Definition at line 1094 of file SampleHandlerFD.cpp.

1094 {
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}
SampleBinningInfo Binning
KS: This stores binning information, in future could be come vector to store binning for every used s...
double ** SampleHandlerFD_array_w2
KS Array used for MC stat.
std::vector< double > YBinEdges
Vector to hold y-axis bin-edges.
std::vector< double > XBinEdges
Vector to hold x-axis bin-edges.

◆ GetXBinVarName()

std::string SampleHandlerFD::GetXBinVarName ( )
inline

Definition at line 38 of file SampleHandlerFD.h.

38{return XVarStr;}

◆ GetYBinVarName()

std::string SampleHandlerFD::GetYBinVarName ( )
inline

Definition at line 40 of file SampleHandlerFD.h.

40{return YVarStr;}

◆ ReturnHistsBySelection1D()

std::vector< TH1 * > SampleHandlerFD::ReturnHistsBySelection1D ( std::string  KinematicProjection,
int  Selection1,
int  Selection2 = -1,
int  WeightStyle = 0,
TAxis *  Axis = 0 
)

Definition at line 2030 of file SampleHandlerFD.cpp.

2030 {
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}
TLegend * THStackLeg
DB Miscellaneous Variables.

◆ ReturnHistsBySelection2D()

std::vector< TH2 * > SampleHandlerFD::ReturnHistsBySelection2D ( std::string  KinematicProjectionX,
std::string  KinematicProjectionY,
int  Selection1,
int  Selection2 = -1,
int  WeightStyle = 0,
TAxis *  XAxis = 0,
TAxis *  YAxis = 0 
)

Definition at line 2066 of file SampleHandlerFD.cpp.

2068 {
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}

◆ ReturnKinematicParameterFromString()

int SampleHandlerFD::ReturnKinematicParameterFromString ( const std::string &  KinematicStr) const

ETA function to generically convert a string from xsec cov to a kinematic type.

Definition at line 1689 of file SampleHandlerFD.cpp.

1689 {
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}
const std::unordered_map< std::string, int > * KinematicParameters
Mapping between string and kinematic enum.
static constexpr const int _BAD_INT_
Default value used for int initialisation.
Definition: Core.h:45

◆ ReturnStackedHistBySelection1D()

THStack * SampleHandlerFD::ReturnStackedHistBySelection1D ( std::string  KinematicProjection,
int  Selection1,
int  Selection2 = -1,
int  WeightStyle = 0,
TAxis *  Axis = 0 
)

Definition at line 2095 of file SampleHandlerFD.cpp.

2095 {
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}
std::string GetTitle() const override
std::vector< TH1 * > ReturnHistsBySelection1D(std::string KinematicProjection, int Selection1, int Selection2=-1, int WeightStyle=0, TAxis *Axis=0)

◆ ReturnStackHistLegend()

TLegend * SampleHandlerFD::ReturnStackHistLegend ( )
inline

Definition at line 126 of file SampleHandlerFD.h.

126{return THStackLeg;}

◆ ReturnStringFromKinematicParameter()

std::string SampleHandlerFD::ReturnStringFromKinematicParameter ( const int  KinematicVariable) const

ETA function to generically convert a kinematic type from xsec cov to a string.

Definition at line 1701 of file SampleHandlerFD.cpp.

1701 {
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}
const std::unordered_map< int, std::string > * ReversedKinematicParameters
Mapping between kinematic enum and string.

◆ SetTestStatistic()

void SampleHandlerBase::SetTestStatistic ( TestStatistic  testStat)
inline

Set the test statistic to be used when calculating the binned likelihoods.

Parameters
testStatThe test statistic to use.

Definition at line 99 of file SampleHandlerBase.h.

99{ fTestStatistic = testStat; }