/*
This is a proprietary part of DataMelt package.
License: the Creative Commons Attribution-Share Alike License; either version 3.0 of the License,
or any later version.
http://creativecommons.org/licenses/by-sa/3.0/
Permission to use, copy and modify this software and its documentation for
NON-COMMERCIAL purposes is granted,
without fee, provided that an acknowledgment to the author, Dr Sergei Chekanov
appears in all copies and associated documentation or publications.
Redistributions of the source code, or parts of the source codes,
must retain the copyright notice, this list of conditions and the following disclaimer
(all at the top of each source code) and requires the written permission of Sergei Chekanov
Redistribution in binary form of all or parts of a class must reproduce the
copyright notice, this list of conditions and the following disclaimer in the documentation and/or other
materials provided with the distribution and requires the written permission of Sergei Chekanov
Details about license agreement or commercial use are: http://jwork.org/jhepwork/license.php
*/
package jhpro.stat.limit;
import java.util.ArrayList;
import java.util.Random;
import org.apache.commons.math3.util.FastMath;
/**
*
* This is a class to give confidence intervals on various measurements
* calculated from the number of events observed, the acceptance factor, the
* background estimate and the associated errors. Results from different
* experiments (or different channels from the same experiment) can be combined.
* A complete description is given Computer Physics Communications 149 (2002)
* 97–102 by Roger Barlow.
*
* After initializing this class, add data (addData()). Then run "run()" method.
* You can get probabilities using getProbabilyty method.
*
*
* The user can iterate towards whatever limit type they choose, e.g., for a 90%
* upper limit they adjust the limit guess until the probability is shown as
* 0.10; for a 90% lower limit they would aim for 0.90.
*
*
* @author Sergei Chekanov
*
*/
public class StatConfidence {
static final double EPSILON = 0.0001;
double sensitivity; // common acceptance-sensitivity factor
double sigma; // absolute error on S
double limitguess; // limit guess
double NMCevents; // Number of MC events
ArrayList expdata;
private Random rrr;
private double R1, R2, R3, Q1, Q2, Q3;
private boolean allzero;
private double pole; // where Eq. 6 may blow up
private double outSWG;
private double outCH;
private double outJ;
/**
* Initialize the limit calculations
*
* @param sensitivit (or efficiency)
* @param sigma (uncertainty on efficiency)
* @param limitguess initial guess on limit
* @param NMCevents Number of Monte Carlo events
*/
public StatConfidence(double sensitivity, double sigma, double limitguess,
int NMCevents) {
this.sensitivity = sensitivity;
this.sigma = sigma;
this.limitguess = limitguess;
this.NMCevents = NMCevents;
rrr = new Random();
expdata = new ArrayList();
R3 = 0;
}
/**
* Initialize the limit calculations, assuming efficiency 1.
*
* @param limitguess initial guess on limit
* @param NMCevents Number of Monte Carlo events
*/
public StatConfidence(double limitguess, int NMCevents) {
this(1.0,0.0,limitguess,NMCevents);
}
/**
* Add experimental data
*
* @param data
* exp data
*/
public void addData(ExpData data) {
expdata.add(data);
}
/**
* Get number of experiments
*/
public void getDataSize() {
expdata.size();
}
/**
* Run limit calculations for upper limit. Similar to run(true).
* @return error
*/
public boolean runUpper(){
return run(true);
}
/**
* Run limit calculations for lower limit. Similar to run(false).
* @return error
*/
public boolean runLower(){
return run(false);
}
/**
* Run limit calculations. You can retrieve probabilities using getProbability* methods.
* The probabilities depend on guessed limit (which is input).
*
*
* @param isUpper
* set to true for the upper limit. Set to false for the lower limit.
* @return
*/
public boolean run(boolean isUpper) {
boolean tmp = true;
double epsilon = 0;
allzero = false;
if (isUpper) {
epsilon = EPSILON;
} else
epsilon = -EPSILON;
pole = 1; // used as flag to ensure it's updated first time through
for (int i = 0; i < expdata.size(); i++) { // read values from screen
if (expdata.get(i).getNobs() != 0) {
allzero = false;
double p = -expdata.get(i).getBack()
/ expdata.get(i).getFraction();
if (p > pole || pole > 0)
pole = p;
}
}
int NMAX = expdata.size();
int[] nobs = new int[NMAX];
for (int i = 0; i < NMAX; i++)
nobs[i] = expdata.get(i).getNobs();
double p = pole / sensitivity; // fold in global sensitivity
R1 = 0.0;
while ((Q1 = q(R1, nobs, sensitivity)) < 0) { // lower edge cannot go
// below pole or Eq 6
// explodes
R1 = (p + R1) / 2.;
}
R2 = 1.0;
while ((Q2 = q(R2, nobs, sensitivity)) > 0) { // upper edge can be as
// high as it likes
R2 *= 2;
}
// R1 and R2 now bracket the solution. Use a boring binary split
int safetycount = 0;
do {
R3 = 0.5 * (R1 + R2);
Q3 = q(R3, nobs, sensitivity);
if (Q1 * Q3 > 0) {
R1 = R3;
Q1 = Q3;
} else {
R2 = R3;
Q2 = Q3;
}
safetycount++;
} while (!((FastMath.abs(Q3) <= 10. * Float.MIN_VALUE)
|| (safetycount > 100) || (FastMath.abs(R1 - R2) <= 10. * Float.MIN_VALUE)));
// R3 is now the critical value - unless allzero is set
int scoreSWG = 0;
int scoreCH = 0;
int scoreJ = 0;
// generate events
int[] resultSWG = new int[NMAX];
int[] resultCH = new int[NMAX];
int[] resultJ = new int[NMAX];
int i, j;
for (j = 0; j < NMCevents; j++) { // repeat many times
boolean gotallzeroSWG = true;
boolean gotallzeroCH = true;
boolean gotallzeroJ = true;
double thing = sigma * rrr.nextGaussian() / sensitivity; // proportional
// error
// on
double aCH = sensitivity * (1 + thing); // S using Gaussian smearing
double aSWG = sensitivity / (1 + thing);
double aJ = sensitivity * FastMath.exp(thing);
for (i = 0; i < NMAX; i++) { // generate a set of results
thing = expdata.get(i).getErrFraction() * rrr.nextGaussian()
/ expdata.get(i).getFraction();
double cCH = expdata.get(i).getFraction() * (1 + thing); // smear
// individual
// factors
double cSWG = expdata.get(i).getFraction() / (1 + thing);
double cJ = expdata.get(i).getFraction() * FastMath.exp(thing);
// smear background
double b = expdata.get(i).getBack()
+ expdata.get(i).getErrBack() * rrr.nextGaussian();
// now the actual Poisson generation
resultSWG[i] = prand(limitguess * aSWG * cSWG + b);
if (resultSWG[i] != 0)
gotallzeroSWG = false;
resultCH[i] = prand(limitguess * aCH * cCH + b);
if (resultCH[i] != 0)
gotallzeroCH = false;
resultJ[i] = prand(limitguess * aJ * cJ + b);
if (resultJ[i] != 0)
gotallzeroJ = false;
}
if (allzero) {// complete null is a special case
if (gotallzeroSWG)
scoreSWG++;
if (gotallzeroCH)
scoreCH++;
if (gotallzeroJ)
scoreJ++;
} else {
if (q(R3, resultSWG, sensitivity) < epsilon)
scoreSWG++;
if (q(R3, resultCH, sensitivity) < epsilon)
scoreCH++;
if (q(R3, resultJ, sensitivity) < epsilon)
scoreJ++;
}
} // end of long computation loop. Present results
outSWG = (double) scoreSWG / NMCevents;
outCH = (double) scoreCH / NMCevents;
outJ = (double) scoreJ / NMCevents;
return tmp;
}
/**
* calculate equation 6 of the paper
*
* @param R
* @param c
* @return
*/
private double q(double R, int[] n, double c) {
double result = 0;
for (int i = 0; i < expdata.size(); i++) {
double cs = c * expdata.get(i).getFraction();
result += (n[i] * cs / (R * cs + expdata.get(i).getBack()) - cs);
}
return result;
}
/**
* Poisson generator
*
* @param f
* lambda
* @return
*/
private int prand(double f) { // Poisson generator
if (f < 0)
return 0; // just in case
double prob = rrr.nextDouble(); // throw the dice
int n = 0;
double probi = FastMath.exp(-f);
while (n < 1000) // just in case
{
prob -= probi;
if (prob <= 0)
return n;
probi *= f / (++n);
}
return 0; // should never happen
}
/**
* get probability (Cousins+Highland) for the given limit.
*
* @return
*/
public double getProbabilityCH() {
return outCH;
}
/**
* Get probability (BaBar SWG) for the given limit.
*
* @return
*/
public double getProbabilitySWG() {
return outSWG;
}
/**
* Probability (Jeffreys) for the given limit.
*
* @return
*/
public double getProbabilityJ() {
return outJ;
}
/**
* Set new limit
*
* @param limitguess
*/
public void setLimitGuess(double limitguess) {
this.limitguess = limitguess;
}
/**
* Set number of MC events
*
* @param NMC
* number of MC events
*/
public void setNMC(double NMC) {
NMCevents = NMC;
}
/**
* Set new sensitivity
*
* @param sensitivity
*/
public void setSensitivity(double sensitivity) {
this.sensitivity = sensitivity;
}
/**
* Set error on sensitivity
*
* @param sigma
*/
public void setSensitivityError(double sigma) {
this.sigma = sigma;
}
/**
* Get limit guess
*
* @return
*/
public double getLimitGuess() {
return limitguess;
}
/**
* Get number of MC events
*
* @return
*/
public double getNMC() {
return NMCevents;
}
/**
* Get experimental data
*
* @return
*/
public ArrayList getExpData() {
return expdata;
}
/**
* Get sensitivity
*
* @return
*/
public double getSensitivity() {
return sensitivity;
}
/**
* Get error on sensitivity
*
* @return
*/
public double getSensitivityError() {
return sigma;
}
}
Ads help maintain this website.