CountingExperimentStatistics.cxx
Go to the documentation of this file.
2 
4 
5 #include <cmath>
6 
7 #include "TMath.h"
8 
9 namespace ana
10 {
11  /// There doesn't seem to be a standard header for this. Implement our own.
12  double factorial(int x)
13  {
14  double ret = 1;
15  for(int i = 1; i <= x; ++i) ret *= i;
16  return ret;
17  }
18 
19  //----------------------------------------------------------------------
20  double PValueToSigma(double p)
21  {
22  return sqrt(2)*TMath::ErfcInverse(p);
23  }
24 
25  //----------------------------------------------------------------------
26  double CountingExperimentPValue(int N, double b0, double sigma)
27  {
28  double p = 0;
29 
30  // We compute the sum from 0 to N-1 rather than the infinite sum from N to
31  // infinity.
32  for(int x = 0; x < N; ++x){
33 
34  if(sigma == 0){
35  // If there's no systematic, can short-circuit the loop below
36  p += pow(b0, x) * exp(-b0) / factorial(x);
37  continue;
38  }
39 
40  // Small steps in the number of sigmas
41  const double ds = .001;
42  // Scan over +/-10 sigma range
43  for(double s = -10; s <= +10; s += ds){
44  // Background fluctuated this many sigma, cannot be below zero
45  const double b = std::max(0., b0 + s*sigma*b0);
46 
47  // Poission probability of seeing x given b
48  const double poi = pow(b, x) * exp(-b) / factorial(x);
49 
50  // Gaussian probability of seeing a fluctuation this large
51  const double gaus = exp(-s*s/2)/sqrt(2*M_PI);
52 
53  // Overall probability of x given b0 and sigma
54  p += poi * gaus * ds;
55  }
56  }
57 
58  // Probability of seeing N or more
59  return 1-p;
60  }
61 
62 
63  //----------------------------------------------------------------------
64  double CountingExperimentSigma(int N, double b0, double sigma)
65  {
66  return PValueToSigma(CountingExperimentPValue(N, b0, sigma));
67  }
68 
69 
70  //----------------------------------------------------------------------
71  double CountingExperimentPValueByLL(int N, double b0, double sigma)
72  {
73  double p = 0;
74 
75  const double LL = LogLikelihood(b0, N);
76 
77  // We compute the sum from 0 to N-1 rather than the infinite sum from N to
78  // infinity.
79  for(int x = 0; x < N+10*b0*sigma+1; ++x){
80 
81  if(sigma == 0){
82  // If there's no systematic, can short-circuit the loop below
83  if(LogLikelihood(b0, x) < LL-1e-6){
84  p += pow(b0, x) * exp(-b0) / factorial(x);
85  }
86  continue;
87  }
88 
89  // Small steps in the number of sigmas
90  const double ds = .001;
91  // Scan over +/-10 sigma range
92  for(double s = -10; s <= +10; s += ds){
93  // Background fluctuated this many sigma, cannot be below zero
94  const double b = std::max(0., b0 + s*sigma*b0);
95 
96  // Poission probability of seeing x given b
97  const double poi = pow(b, x) * exp(-b) / factorial(x);
98 
99  // Gaussian probability of seeing a fluctuation this large
100  const double gaus = exp(-s*s/2)/sqrt(2*M_PI);
101 
102  if(LogLikelihood(b0, x) < LL-1e-6){
103  // Overall probability of x given b0 and sigma
104  p += poi * gaus * ds;
105  }
106  } // end for s
107  }
108 
109  // Probability of seeing N or more
110  return 1-p;
111  }
112 
113  //----------------------------------------------------------------------
114  double CountingExperimentSigmaByLL(int N, double b0, double sigma)
115  {
116  return PValueToSigma(CountingExperimentPValueByLL(N, b0, sigma));
117  }
118 }
T max(const caf::Proxy< T > &a, T b)
Cuts and Vars for the 2020 FD DiF Study.
Definition: vars.h:6
const char * p
Definition: xmltok.h:285
T sqrt(T number)
Definition: d0nt_math.hpp:156
constexpr T pow(T x)
Definition: pow.h:75
double CountingExperimentSigmaByLL(int N, double b0, double sigma)
double LogLikelihood(const Eigen::ArrayXd &ea, const Eigen::ArrayXd &oa, bool useOverflow)
The log-likelihood formula from the PDG.
Definition: EigenUtils.cxx:36
#define M_PI
Definition: SbMath.h:34
double PValueToSigma(double p)
Compute the equivalent number of gaussian sigma for a p-value.
const XML_Char * s
Definition: expat.h:262
fvar< T > exp(const fvar< T > &x)
Definition: exp.hpp:10
double CountingExperimentPValueByLL(int N, double b0, double sigma)
A very simple service to remember what detector we&#39;re working in.
double sigma(TH1F *hist, double percentile)
double CountingExperimentSigma(int N, double b0, double sigma)
double CountingExperimentPValue(int N, double b0, double sigma)
Compute the probability of seeing N or more events.
double factorial(int x)
There doesn&#39;t seem to be a standard header for this. Implement our own.
const hit & b
Definition: hits.cxx:21
Float_t e
Definition: plot.C:35