/*
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/scavis
*/
package jhpro.stat;
import java.awt.Color;
import java.util.*;
import jhplot.*;
import jhplot.math.SortUtils;
import jhplot.gui.*;
import org.apache.commons.math3.util.FastMath;
/**
* Confidence level calculations.
*
* For discoveries, 1-CLb indicates the probability that the background
* fluctuates to produce a distribution of candidates at least
* as signal-like as those observed in the data. For discovery,
* 1-CLb is required to be no more than 2.87x10-7, or twice that,
* depending on how one interprets what is meant by “five sigma,” including just one side of a
* Gaussian tail or both. A “three sigma” excess is defined to be 1-CLb = 1.3x10-3 or twice that.
* But forming discovery p-values, we must compute 1-CLb values of the order
* of 10-7. This computation involves generating of the order of 10^8
* pseudoexperiments, just to be on the safe side.
*
* Read Reference: HEP-EX/9902006. see: Tom Junk,NIM A434, p. 435-443,
*
* @author S.Chekanov
*
*/
public class ConfidenceLevel {
double fgMCLM2S = 0.025;
double fgMCLM1S = 0.16;
double fgMCLMED = 0.5;
double fgMCLP1S = 0.84;
double fgMCLP2S = 0.975;
// LHWG "one-sided" definition
double fgMCL3S1S = 2.6998E-3;
double fgMCL5S1S = 5.7330E-7;
// the other definition (not chosen by the LHWG)
double fgMCL3S2S = 1.349898E-3;
double fgMCL5S2S = 2.866516E-7;
double fStot;
double fBtot;
double fDtot;
double fTSD;
int fNMC;
int fNNMC;
double fMCL3S;
double fMCL5S;
double[] fTSB;
double[] fTSS;
double[] fLRS;
double[] fLRB;
int[] fISS;
int[] fISB;
/**
* Default constructor.
*/
public ConfidenceLevel() {
fStot = 0;
fBtot = 0;
fDtot = 0;
fTSD = 0;
fTSB = null;
fTSS = null;
fLRS = null;
fLRB = null;
fNMC = 0;
fNNMC = 0;
fISS = new int[1];
fISB = new int[1];
fMCL3S = fgMCL3S1S;
fMCL5S = fgMCL5S1S;
}
/**
* Construct ConfLevel
* @param mc number of MonteCarlo experiments
*/
public ConfidenceLevel(int mc) {
this(mc, true);
}
/**
* Build confidence level.
*
* @param mc
* is the number of Monte Carlo experiments
* @param onesided
* specifies if the intervals are one-sided or not.
*/
public ConfidenceLevel(int mc, boolean onesided) {
fStot = 0;
fBtot = 0;
fDtot = 0;
fTSD = 0;
fTSB = null;
fTSS = null;
fLRS = null;
fLRB = null;
fNMC = mc;
fNNMC = mc;
fISS = new int[mc];
fISB = new int[mc];
fMCL3S = onesided ? fgMCL3S1S : fgMCL3S2S;
fMCL5S = onesided ? fgMCL5S1S : fgMCL5S2S;
}
/**
* Get the expected statistic value in the background only hypothesis
* @param sigma between -2 and 2
* @return
*/
public double getExpectedStatistic_b(int sigma) {
switch (sigma) {
case -2:
return (-2 * ((fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLP2S)))]]) - fStot));
case -1:
return (-2 * ((fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLP1S)))]]) - fStot));
case 0:
return (-2 * ((fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLMED)))]]) - fStot));
case 1:
return (-2 * ((fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLM1S)))]]) - fStot));
case 2:
return (-2 * ((fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLM2S)))]]) - fStot));
default:
return 0;
}
}
/**
* Get the expected statistic value in the signal plus background hypothesis
*
* @param sigma
* @return
*/
public double getExpectedStatistic_sb(int sigma) {
switch (sigma) {
case -2:
return (-2 * ((fTSS[fISS[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLP2S)))]]) - fStot));
case -1:
return (-2 * ((fTSS[fISS[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLP1S)))]]) - fStot));
case 0:
return (-2 * ((fTSS[fISS[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLMED)))]]) - fStot));
case 1:
return (-2 * ((fTSS[fISS[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLM1S)))]]) - fStot));
case 2:
return (-2 * ((fTSS[fISS[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLM2S)))]]) - fStot));
default:
return 0;
}
}
/**
* Get the Confidence Level for the background only.
* This confidence level quantifies the confidence of a potential
* discovery, as it expresses the probability that background processes would
* give fewer than or equal to the number of candidates observed.
* 1-CLb is the probability that the null hypothesis will give an outcome
* that looks at least as signal-like as the one observed. For discovery,
* 1-CLb is required to be no more than 2.87x10-7, or twice that,
* depending on how one interprets what is meant by “five sigma,” including just one side of a
* Gaussian tail or both. A “three sigma” excess is defined to be 1-CLb=1.3x10-3 or twice that.
* But forming discovery p-values, we must compute 1-CLb values of the order
* of 10-7. This computation involves generating of the order of 10^8 pseudoexperiments,
* just to be on the safe side.
*
* @return Confidence Level for the background only.
*/
public double getCLb() {
return getCLb(false);
}
/**
* Get the Confidence Level for the background only.
* This confidence level quantifies the confidence of a potential
* discovery, as it expresses the probability that background processes would
* give fewer than or equal to the number of candidates observed.
* 1-CLb is the probability that the null hypothesis will give an outcome
* that looks at least as signal-like as the one observed. For discovery,
* 1-CLb is required to be no more than 2.87*10-7, or twice that,
* depending on how one interprets what is meant by “five sigma,” including just one side of a
* Gaussian tail or both. A “three sigma” excess is defined to be 1-CLb = 1.3*10^-3 or twice that.
*
* But forming discovery p-values, we must compute 1-CLb values of the order
* of 10-7. This computation involves generating of the order of 10^8 pseudoexperiments,
* just to be on the safe side.
*
* @param use_sMC
* @return Confidence Level for the background only.
*/
public double getCLb(boolean use_sMC) {
int use = 0;
if (use_sMC==false)
use = 1;
if (use_sMC==true)
use = 2;
double result = 0;
switch (use) {
case 1: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] < fTSD)
result = ((double) (i + 1)) / fNMC;
return result;
}
case 2: {
for (int i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] < fTSD)
result += (1 / (fLRS[fISS[i]] * fNMC));
return result;
}
}
return result;
}
/**
* Get the Confidence Level for the signal plus background hypothesis
* The confidence level for excluding the possibility of simultaneous presence
* of new particle production and background (the s + b hypothesis)
*
* @return
*/
public double getCLsb() {
return getCLsb(false);
}
/**
* Get the Confidence Level for the signal plus background hypothesis
*
* @param use_sMC
* @return
*/
public double getCLsb(boolean use_sMC) {
int use = 0;
if (use_sMC == false)
use = 1;
if (use_sMC == true)
use = 2;
double result = 0;
switch (use) {
case 1: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSD)
result += (fLRB[fISB[i]]) / fNMC;
return result;
}
case 2: {
for (int i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSD)
result = i / fNMC;
return result;
}
}
return result;
}
/**
* Get the Confidence Level defined by CLs = CLsb/CLb.
* This quantity is
* stable w.r.t. background fluctuations.
*
* This hypothesis is excluded at the 95% CL if CLs = 0.05,
* and at more than the 95% CL if CLs < 0.05, assuming that signal is present.
*
*
* @return
*/
public double getCLs() {
return getCLs(false);
}
/**
* Get the Confidence Level defined by CLs = CLsb/CLb. This quantity is
* stable w.r.t. background fluctuations.
* This hypothesis is excluded at the 95% CL if CLs = 0.05,
* and at more than the 95% CL if CLs < 0.05, assuming that signal is present.
*
*
* @param use_sMC use or not MC.
* @return
*/
public double getCLs(boolean use_sMC) {
double clb = getCLb(false);
double clsb = getCLsb(use_sMC);
if (clb == 0) {
System.out.println("Warning: clb = 0 !");
return 0;
} else
return clsb / clb;
}
/**
* Get the expected Confidence Level for the signal plus background
* hypothesis if there is only background.
*
* @param sigma
* @return
*/
public double getExpectedCLsb_b(int sigma) {
double result = 0;
switch (sigma) {
case -2: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLP2S)))]])
result += fLRB[fISB[i]] / fNMC;
return result;
}
case -1: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLP1S)))]])
result += fLRB[fISB[i]] / fNMC;
return result;
}
case 0: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLMED)))]]) {
result += fLRB[fISB[i]] / fNMC;
// System.out.println(fISB[i]);
// System.out.println(fLRB[fISB[i]]); // sergei
}
return result;
}
case 1: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLM1S)))]])
result += fLRB[fISB[i]] / fNMC;
return result;
}
case 2: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLM2S)))]])
result += fLRB[fISB[i]] / fNMC;
return result;
}
default:
return 0;
}
}
/**
* Get the expected Confidence Level for the background only if there is
* signal and background.
*
* @param sigma
* @return
*/
public double getExpectedCLb_sb(int sigma) {
double result = 0;
switch (sigma) {
case 2: {
for (int i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSS[fISS[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLP2S)))]])
result += fLRS[fISS[i]] / fNMC;
return result;
}
case 1: {
for (int i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSS[fISS[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLP1S)))]])
result += fLRS[fISS[i]] / fNMC;
return result;
}
case 0: {
for (int i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSS[fISS[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLMED)))]])
result += fLRS[fISS[i]] / fNMC;
return result;
}
case -1: {
for (int i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSS[fISS[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLM1S)))]])
result += fLRS[fISS[i]] / fNMC;
return result;
}
case -2: {
for (int i = 0; i < fNMC; i++)
if (fTSS[fISS[i]] <= fTSS[fISS[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLM2S)))]])
result += fLRS[fISS[i]] / fNMC;
return result;
}
default:
return 0;
}
}
/**
* Get the expected Confidence Level for the background only if there is
* only background.
*
* @param sigma
* @return
*/
public double getExpectedCLb_b(int sigma) {
double result = 0;
switch (sigma) {
case 2: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLM2S)))]])
result = (i + 1) / (double) fNMC;
return result;
}
case 1: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLM1S)))]])
result = (i + 1) / (double) fNMC;
return result;
}
case 0: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLMED)))]])
result = (i + 1) / (double) fNMC;
return result;
}
case -1: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLP1S)))]])
result = (i + 1) / (double) fNMC;
return result;
}
case -2: {
for (int i = 0; i < fNMC; i++)
if (fTSB[fISB[i]] <= fTSB[fISB[FastMath.min((int) fNMC,
(int) FastMath.max((int) 1, (int) (fNMC * fgMCLP2S)))]])
result = (i + 1) / (double) fNMC;
return result;
}
}
return result;
}
/**
* Get average CLsb.
*
* @return
*/
public double getAverageCLsb() {
double result = 0;
double psumsb = 0;
for (int i = 0; i < fNMC; i++) {
psumsb += fLRB[fISB[i]] / fNMC;
result += psumsb / fNMC;
}
return result;
}
/**
* Get average CLs.
*/
public double getAverageCLs() {
double result = 0;
double psumsb = 0;
for (int i = 0; i < fNMC; i++) {
psumsb += fLRB[fISB[i]] / fNMC;
result += ((psumsb / fNMC) / ((i + 1) / fNMC));
}
return result;
}
/**
* Get 3s probability.
*/
public double get3sProbability() {
double result = 0;
double psumbs = 0;
for (int i = 0; i < fNMC; i++) {
psumbs += 1 / (double) (fLRS[(fISS[(int) (fNMC - i)])] * fNMC);
if (psumbs <= fMCL3S)
result = i / fNMC;
}
return result;
}
/**
* Get 5s probability.
*/
public double get5sProbability() {
double result = 0;
double psumbs = 0;
for (int i = 0; i < fNMC; i++) {
psumbs += 1 / (double) (fLRS[(fISS[(int) (fNMC - i)])] * fNMC);
if (psumbs <= fMCL5S)
result = i / fNMC;
}
return result;
}
/**
* Display sort of a "canonical" -2lnQ plot. This results in a plot with 2
* elements: // - The histogram of -2lnQ for background hypothesis (full) -
* The histogram of -2lnQ for signal and background hypothesis (dashed) The
* 2 histograms are respectively named b_hist and sb_hist.
*
* @param Option
* @return
*/
public ArrayList getResults(String Option) {
H1D h = new H1D("TConfidenceLevel_Draw", 50, 0, 0);
int i;
for (i = 0; i < fNMC; i++) {
h.fill(-2 * (fTSB[i] - fStot));
h.fill(-2 * (fTSS[i] - fStot));
}
H1D b_hist = new H1D("-2lnQ", 50, h.getMin(), h.getMax());
H1D sb_hist = new H1D("-2lnQ", 50, h.getMin(), h.getMax());
for (i = 0; i < fNMC; i++) {
b_hist.fill(-2 * (fTSB[i] - fStot));
sb_hist.fill(-2 * (fTSS[i] - fStot));
}
sb_hist.setColor(Color.red);
ArrayList arr = new ArrayList();
arr.add(b_hist);
arr.add(sb_hist);
return arr;
}
public void setTSD(double in) {
fTSD = in;
}
public void setLRS(double[] in) {
fLRS = in;
}
public void setLRB(double[] in) {
fLRB = in;
}
public void setBtot(double in) {
fBtot = in;
}
public void setStot(double in) {
fStot = in;
}
public void setDtot(int in) {
fDtot = in;
}
public double getStatistic() {
return -2 * (fTSD - fStot);
}
public void setTSB(double[] in) {
// Set the TSB.
fTSB = in;
// make a clone and sort in increasing order
Double[] clone = new Double[in.length];
for (int i=0; i
Ads help maintain this website.