Jpp 21.0.0-rc.1
the software that should make you happy
Loading...
Searching...
No Matches
Classes | Public Member Functions | Static Public Member Functions | Static Public Attributes | Protected Member Functions | Protected Attributes | Static Private Member Functions | List of all members
JASTRONOMY::JRealExperiment Struct Reference

Real experiment using PDF of signal and background. More...

#include <JRealExperiment.hh>

Inheritance diagram for JASTRONOMY::JRealExperiment:
JASTRONOMY::JExperiment JASTRONOMY::JAspera std::vector< double >

Classes

struct  remnant_type
 

Public Member Functions

 JRealExperiment ()
 Default constructor.
 
template<class H_t >
 JRealExperiment (const H_t &hd, const H_t &hs, const H_t &hb)
 Constructor.
 
void add (const TObject *pd, const TObject *ps, const TObject *pb)
 Add objects with data and PDFs of signal and background.
 
void add (const TH1 &hd, const TH1 &hs, const TH1 &hb)
 Add histograms with data and PDFs of signal and background.
 
void add (const TH2 &hd, const TH2 &hs, const TH2 &hb)
 Add histograms with data and PDFs of signal and background.
 
void add (const TH3 &hd, const TH3 &hs, const TH3 &hb)
 Add histograms with data and PDFs of signal and background.
 
void add (const size_t n, const double s, const double b)
 Add data, signal and background.
 
void add ()
 Add remnant data, signal and background.
 
void put (const double s, const double b)
 Put signal and background to list of pre-computed N/S values.
 
void put (const size_t n, const double s, const double b)
 Put signal and background to list of pre-computed N/S values.
 
double getLikelihood (const double p) const
 Get likelihood for given signal strength.
 
double getDerivative (const double p) const
 Get derivative of likelihood for given signal strength.
 
fit_type operator() (const bool ns=false) const
 Fit signal strength.
 
double getTestStatisticForUpperLimit (const double ps) const
 Get test statistic for given signal strength.
 
double getSignal () const
 Get total signal strength.
 
void setSignal (const double wS)
 Set signal strength.
 
void addSignal (const double wS)
 Add signal strength.
 

Static Public Member Functions

static double getSNR ()
 Get minimal signa-to-noise ratio.
 
static void setSNR (const double value)
 Set minimal signa-to-noise ratio.
 
static bool check (const double s, const double b)
 Check validity of signal and background.
 
static bool check (const TAxis *ha, const TAxis *hb)
 Check histogram bins.
 
static bool check (const TH1 &ha, const TH1 &hb)
 Check histogram bins.
 
static bool check (const TH2 &ha, const TH2 &hb)
 Check histogram bins.
 
static bool check (const TH3 &ha, const TH3 &hb)
 Check histogram bins.
 

Static Public Attributes

static constexpr double EPSILON = 1.0e-3
 precision determination of signal strength
 

Protected Member Functions

template<class H_t >
bool add (const TObject *pd, const TObject *ps, const TObject *pb)
 Add objects with data and PDFs of signal and background.
 

Protected Attributes

struct JASTRONOMY::JRealExperiment::remnant_type remnant
 
double ws = 0.0
 total signal strength
 

Static Private Member Functions

static double & get_snr ()
 Get minimal signa-to-noise ratio.
 

Detailed Description

Real experiment using PDF of signal and background.

Definition at line 29 of file JRealExperiment.hh.

Constructor & Destructor Documentation

◆ JRealExperiment() [1/2]

JASTRONOMY::JRealExperiment::JRealExperiment ( )
inline

Default constructor.

Definition at line 36 of file JRealExperiment.hh.

37 {}

◆ JRealExperiment() [2/2]

template<class H_t >
JASTRONOMY::JRealExperiment::JRealExperiment ( const H_t & hd,
const H_t & hs,
const H_t & hb )
inline

Constructor.

Parameters
hdhistogram with data
hshistogram with PDF of signal
hbhistogram with PDF of background

Definition at line 48 of file JRealExperiment.hh.

51 {
52 add(hd, hs, hb);
53 }
JASTRONOMY::JRealExperiment::add
void add()
Add remnant data, signal and background.
Definition JRealExperiment.hh:158

Member Function Documentation

◆ add() [1/7]

void JASTRONOMY::JRealExperiment::add ( const TObject * pd,
const TObject * ps,
const TObject * pb )
inline

Add objects with data and PDFs of signal and background.

Parameters
pdpointer to object with data
pspointer to object with PDF of signal
pbpointer to object with PDF of background

Definition at line 63 of file JRealExperiment.hh.

66 {
67 if (add<TH3>(pd, ps, pb)) { return; }
68 if (add<TH2>(pd, ps, pb)) { return; }
69 if (add<TH1>(pd, ps, pb)) { return; }
70 }

◆ add() [2/7]

void JASTRONOMY::JRealExperiment::add ( const TH1 & hd,
const TH1 & hs,
const TH1 & hb )
inline

Add histograms with data and PDFs of signal and background.

Parameters
hdhistogram with data
hshistogram with PDF of signal
hbhistogram with PDF of background

Definition at line 80 of file JRealExperiment.hh.

83 {
84 for (Int_t ix = 1; ix <= hs.GetXaxis()->GetNbins(); ++ix) {
85 add(hd.GetBinContent(ix),
86 hs.GetBinContent(ix),
87 hb.GetBinContent(ix));
88 }
89 }

◆ add() [3/7]

void JASTRONOMY::JRealExperiment::add ( const TH2 & hd,
const TH2 & hs,
const TH2 & hb )
inline

Add histograms with data and PDFs of signal and background.

Parameters
hdhistogram with data
hshistogram with PDF of signal
hbhistogram with PDF of background

Definition at line 99 of file JRealExperiment.hh.

102 {
103 for (Int_t ix = 1; ix <= hs.GetXaxis()->GetNbins(); ++ix) {
104 for (Int_t iy = 1; iy <= hs.GetYaxis()->GetNbins(); ++iy) {
105 add(hd.GetBinContent(ix, iy),
106 hs.GetBinContent(ix, iy),
107 hb.GetBinContent(ix, iy));
108 }
109 }
110 }

◆ add() [4/7]

void JASTRONOMY::JRealExperiment::add ( const TH3 & hd,
const TH3 & hs,
const TH3 & hb )
inline

Add histograms with data and PDFs of signal and background.

Parameters
hdhistogram with data
hshistogram with PDF of signal
hbhistogram with PDF of background

Definition at line 120 of file JRealExperiment.hh.

123 {
124 for (Int_t ix = 1; ix <= hs.GetXaxis()->GetNbins(); ++ix) {
125 for (Int_t iy = 1; iy <= hs.GetYaxis()->GetNbins(); ++iy) {
126 for (Int_t iz = 1; iz <= hs.GetZaxis()->GetNbins(); ++iz) {
127 add(hd.GetBinContent(ix, iy, iz),
128 hs.GetBinContent(ix, iy, iz),
129 hb.GetBinContent(ix, iy, iz));
130 }
131 }
132 }
133 }

◆ add() [5/7]

void JASTRONOMY::JRealExperiment::add ( const size_t n,
const double s,
const double b )
inline

Add data, signal and background.

Parameters
ndata
ssignal
bbackground

Definition at line 143 of file JRealExperiment.hh.

146 {
147 if (check(s, b)) {
148 put(n, s, b);
149 } else {
150 remnant.add(n, s, b);
151 }
152 }
JASTRONOMY::JAspera::put
void put(const double s, const double b)
Put signal and background to list of pre-computed N/S values.
Definition JAspera.hh:44
JASTRONOMY::JExperiment::check
static bool check(const double s, const double b)
Check validity of signal and background.
Definition JExperiment.hh:58
JASTRONOMY::JRealExperiment::remnant_type::add
void add(const size_t n, const double s, const double b)
Definition JRealExperiment.hh:204
JASTRONOMY::JRealExperiment::remnant
struct JASTRONOMY::JRealExperiment::remnant_type remnant

◆ add() [6/7]

void JASTRONOMY::JRealExperiment::add ( )
inline

Add remnant data, signal and background.

Definition at line 158 of file JRealExperiment.hh.

159 {
160 if (remnant.n != 0) {
161
162 add(remnant.n, remnant.s, remnant.b);
163
164 remnant.reset();
165 }
166 }
JASTRONOMY::JRealExperiment::remnant_type::b
double b
Definition JRealExperiment.hh:220
JASTRONOMY::JRealExperiment::remnant_type::n
size_t n
Definition JRealExperiment.hh:218
JASTRONOMY::JRealExperiment::remnant_type::reset
void reset()
Definition JRealExperiment.hh:211
JASTRONOMY::JRealExperiment::remnant_type::s
double s
Definition JRealExperiment.hh:219

◆ add() [7/7]

template<class H_t >
bool JASTRONOMY::JRealExperiment::add ( const TObject * pd,
const TObject * ps,
const TObject * pb )
inlineprotected

Add objects with data and PDFs of signal and background.

Parameters
pdpointer to object with data
pspointer to object with PDF of signal
pbpointer to object with PDF of background
Returns
true if added; else false

Definition at line 178 of file JRealExperiment.hh.

181 {
182 if (dynamic_cast<const H_t*>(pd) != NULL &&
183 dynamic_cast<const H_t*>(ps) != NULL &&
184 dynamic_cast<const H_t*>(pb) != NULL) {
185
186 const H_t& hd = dynamic_cast<const H_t&>(*pd);
187 const H_t& hs = dynamic_cast<const H_t&>(*ps);
188 const H_t& hb = dynamic_cast<const H_t&>(*pb);
189
190 if (check(hd, hs) && check(hs, hb)) {
191
192 add(hd, hs, hb);
193
194 return true;
195 }
196 }
197
198 return false;
199 }

◆ getSNR()

static double JASTRONOMY::JExperiment::getSNR ( )
inlinestaticinherited

Get minimal signa-to-noise ratio.

Returns
signa-to-noise ratio.

Definition at line 33 of file JExperiment.hh.

34 {
35 return get_snr();
36 }
JASTRONOMY::JExperiment::get_snr
static double & get_snr()
Get minimal signa-to-noise ratio.
Definition JExperiment.hh:125

◆ setSNR()

static void JASTRONOMY::JExperiment::setSNR ( const double value)
inlinestaticinherited

Set minimal signa-to-noise ratio.

Parameters
valuesigna-to-noise ratio.

Definition at line 44 of file JExperiment.hh.

45 {
46 get_snr() = value;
47 }

◆ check() [1/5]

static bool JASTRONOMY::JExperiment::check ( const double s,
const double b )
inlinestaticinherited

Check validity of signal and background.

Parameters
ssignal
bbackground compare significance expectation for 1 bin with cut value
Returns
true if signal and backgroud are valid; else false

Definition at line 58 of file JExperiment.hh.

59 {
60 return (s > 0.0 && b > 0.0 && (s+s*s)/(b+b*b) >= getSNR()); //E(S^2) / E(B^2) for Poisson distribution
61 }
JASTRONOMY::JExperiment::getSNR
static double getSNR()
Get minimal signa-to-noise ratio.
Definition JExperiment.hh:33

◆ check() [2/5]

static bool JASTRONOMY::JExperiment::check ( const TAxis * ha,
const TAxis * hb )
inlinestaticinherited

Check histogram bins.

Parameters
hahistogram axis
hbhistogram axis
Returns
true if same binning; else false

Definition at line 71 of file JExperiment.hh.

72 {
73 return (ha->GetNbins() == hb->GetNbins() &&
74 ha->GetXmin() == hb->GetXmin() &&
75 ha->GetXmax() == hb->GetXmax());
76 }

◆ check() [3/5]

static bool JASTRONOMY::JExperiment::check ( const TH1 & ha,
const TH1 & hb )
inlinestaticinherited

Check histogram bins.

Parameters
hahistogram
hbhistogram
Returns
true if same binning; else false

Definition at line 85 of file JExperiment.hh.

86 {
87 return check(ha.GetXaxis(), hb.GetXaxis());
88 }

◆ check() [4/5]

static bool JASTRONOMY::JExperiment::check ( const TH2 & ha,
const TH2 & hb )
inlinestaticinherited

Check histogram bins.

Parameters
hahistogram
hbhistogram
Returns
true if same binning; else false

Definition at line 98 of file JExperiment.hh.

99 {
100 return (check(ha.GetXaxis(), hb.GetXaxis()) &&
101 check(ha.GetYaxis(), hb.GetYaxis()));
102 }

◆ check() [5/5]

static bool JASTRONOMY::JExperiment::check ( const TH3 & ha,
const TH3 & hb )
inlinestaticinherited

Check histogram bins.

Parameters
hahistogram
hbhistogram
Returns
true if same binning; else false

Definition at line 112 of file JExperiment.hh.

113 {
114 return (check(ha.GetXaxis(), hb.GetXaxis()) &&
115 check(ha.GetYaxis(), hb.GetYaxis()) &&
116 check(ha.GetZaxis(), hb.GetZaxis()));
117 }

◆ get_snr()

static double & JASTRONOMY::JExperiment::get_snr ( )
inlinestaticprivateinherited

Get minimal signa-to-noise ratio.

Returns
signa-to-noise ratio.

Definition at line 125 of file JExperiment.hh.

126 {
127 static double value = 0.0;
128
129 return value;
130 }

◆ put() [1/2]

void JASTRONOMY::JAspera::put ( const double s,
const double b )
inlineinherited

Put signal and background to list of pre-computed N/S values.

Parameters
ssignal
bbackground

Definition at line 44 of file JAspera.hh.

46 {
47 push_back(b/s);
48
49 ws += s;
50 }
JASTRONOMY::JAspera::ws
double ws
total signal strength
Definition JAspera.hh:378

◆ put() [2/2]

void JASTRONOMY::JAspera::put ( const size_t n,
const double s,
const double b )
inlineinherited

Put signal and background to list of pre-computed N/S values.

Parameters
ndata
ssignal
bbackground

Definition at line 60 of file JAspera.hh.

63 {
64 for (size_t i = 0; i != n; ++i) {
65 push_back(b/s);
66 }
67
68 ws += s;
69 }
JTOOLS::n
const int n
Definition JPolint.hh:791

◆ getLikelihood()

double JASTRONOMY::JAspera::getLikelihood ( const double p) const
inlineinherited

Get likelihood for given signal strength.

Parameters
psignal strength
Returns
likelihood

Definition at line 78 of file JAspera.hh.

79 {
80 double y = -p * ws;
81
82 for (const double i : static_cast<const std::vector<double>&>(*this)) {
83 y += log1p(p/i);
84 }
85
86 return y;
87 }

◆ getDerivative()

double JASTRONOMY::JAspera::getDerivative ( const double p) const
inlineinherited

Get derivative of likelihood for given signal strength.

Parameters
psignal strength
Returns
derivative of likelihood

Definition at line 96 of file JAspera.hh.

97 {
98 double y = -ws;
99
100 for (const double i : static_cast<const std::vector<double>&>(*this)) {
101 y += 1.0 / (p + i);
102 }
103
104 return y;
105 }

◆ operator()()

fit_type JASTRONOMY::JAspera::operator() ( const bool ns = false) const
inlineinherited

Fit signal strength.

Parameters
nsallow for negative signal
Returns
result

Definition at line 149 of file JAspera.hh.

150 {
151 using namespace std;
152
153 if (this->empty()) {
154
155 // nothing to be done
156
157 return { 0.0, 0.0 };
158
159 } else if (this->size() == 1 ) {
160
161 // analytical solution
162
163 const double x = 1.0/ws - (*this)[0];
164
165 if (x > 0.0 || ns)
166 return { getLikelihood(x), x };
167 else
168 return { 0.0, 0.0 };
169 }
170
171 double x1 = 0.0;
172 double x2 = 0.0;
173
174 double f1 = getDerivative(0.0); // discriminator between positive and negative signal
175 double f2 = 0.0;
176
177 if (f1 == 0.0) {
178
179 return { 0.0, 0.0 };
180
181 } else if (f1 > 0.0) { // positive signal
182
183 x1 = 0.0; // lower limit corresponds to no signal
184 x2 = (double) this->size() / ws; // upper limit corresponds to no background (i.e. all N/S = 0)
185
186 f2 = getDerivative(x2);
187
188 } else if (ns) { // negative signal
189
190 x2 = 0.0; // upper limit corresponds to no signal
191 x1 = -(*this)[0]; // lower limit corresponds to largest negated N/S ratio
192
193 for (const double x : static_cast<const std::vector<double>&>(*this)) {
194 if (-x > x1) {
195 x1 = -x;
196 }
197 }
198
199 x1 += 1.0 / ws; // offset
200
201 f2 = f1;
202 f1 = getDerivative(x1);
203
204 } else {
205
206 return { 0.0, 0.0 };
207 }
208
209#if defined(BISECTION)
210
211 // binary search
212
213 while (x2 - x1 > EPSILON) {
214
215 const double x = 0.5 * (x1 + x2);
216
217 if (getDerivative(x) > 0.0)
218 x1 = x;
219 else
220 x2 = x;
221 }
222
223 const double x = 0.5 * (x1 + x2);
224
225 return { getLikelihood(x), x };
226
227#elif defined(NEWTON_RAPHSON)
228
229 // Newton-Raphson method
230
231 double dx = x2 - x1;
232 double x = 0.5 * (x1 + x2);
233
234 derivatives_type y = getDerivatives(x);
235
236 struct {
237 double dx;
238 } old;
239
240 old.dx = dx;
241
242 for ( ; ; ) {
243
244 if ((y.fp + (x2 - x)*y.fpp)*(y.fp + (x1 - x)*y.fpp) > 0.0 || fabs(2.0*y.fp) > fabs(old.dx*y.fpp)) {
245
246 old.dx = dx;
247
248 dx = 0.5 * (x2 - x1);
249 x = x1 + dx;
250
251 } else {
252
253 old.dx = dx;
254
255 dx = y.fp/y.fpp;
256 x -= dx;
257 }
258
259 if (fabs(dx) <= EPSILON) {
260 break;
261 }
262
263 y = getDerivatives(x);
264
265 if (y.fp < 0.0)
266 x2 = x;
267 else
268 x1 = x;
269 }
270
271 return { getLikelihood(x), x };
272
273#else
274
275 // Ridder's method
276
277 while (x2 - x1 > EPSILON) {
278
279 const double xm = 0.5 * (x1 + x2);
280 const double fm = getDerivative(xm);
281
282 const double s = sqrt(fm*fm - f1*f2);
283
284 if (s == 0.0) {
285 break;
286 }
287
288 const double xn = xm + (xm - x1) * fm/s;
289 const double fn = getDerivative(xn);
290
291 if (fn == 0.0) {
292 return { getLikelihood(xn), xn };
293 }
294
295 if (signbit(fn) != signbit(fm)) {
296
297 x1 = xm;
298 f1 = fm;
299 x2 = xn;
300 f2 = fn;
301
302 } else {
303
304 if (signbit(fn)) {
305
306 x2 = xn;
307 f2 = fn;
308
309 } else {
310
311 x1 = xn;
312 f1 = fn;
313 }
314 }
315 }
316
317 const double x = 0.5 * (x1 + x2);
318
319 return { getLikelihood(x), x };
320#endif
321 }
JASTRONOMY::JAspera::getLikelihood
double getLikelihood(const double p) const
Get likelihood for given signal strength.
Definition JAspera.hh:78
JASTRONOMY::JAspera::EPSILON
static constexpr double EPSILON
precision determination of signal strength
Definition JAspera.hh:26
JASTRONOMY::JAspera::getDerivative
double getDerivative(const double p) const
Get derivative of likelihood for given signal strength.
Definition JAspera.hh:96

◆ getTestStatisticForUpperLimit()

double JASTRONOMY::JAspera::getTestStatisticForUpperLimit ( const double ps) const
inlineinherited

Get test statistic for given signal strength.

See formula 16 in this reference.

Parameters
pssignal strength
Returns
test statistic

Definition at line 331 of file JAspera.hh.

332 {
333 const fit_type result = (*this)(true);
334
335 if (result.signal <= 0.0)
336 return 0.0 - this->getLikelihood(ps);
337 else if (result.signal <= ps)
338 return result.likelihood - this->getLikelihood(ps);
339 else
340 return 0.0;
341 }

◆ getSignal()

double JASTRONOMY::JAspera::getSignal ( ) const
inlineinherited

Get total signal strength.

Returns
signal strength

Definition at line 349 of file JAspera.hh.

350 {
351 return ws;
352 }

◆ setSignal()

void JASTRONOMY::JAspera::setSignal ( const double wS)
inlineinherited

Set signal strength.

Parameters
wSsignal strength

Definition at line 360 of file JAspera.hh.

361 {
362 ws = wS;
363 }

◆ addSignal()

void JASTRONOMY::JAspera::addSignal ( const double wS)
inlineinherited

Add signal strength.

Parameters
wSsignal strength

Definition at line 371 of file JAspera.hh.

372 {
373 ws += wS;
374 }

Member Data Documentation

◆ remnant

struct JASTRONOMY::JRealExperiment::remnant_type JASTRONOMY::JRealExperiment::remnant
protected

◆ EPSILON

double JASTRONOMY::JAspera::EPSILON = 1.0e-3
staticconstexprinherited

precision determination of signal strength

Definition at line 26 of file JAspera.hh.

◆ ws

double JASTRONOMY::JAspera::ws = 0.0
protectedinherited

total signal strength

Definition at line 378 of file JAspera.hh.

The documentation for this struct was generated from the following file:
Generated by doxygen 1.12.0