 Upper limit on a 0 result

Hello I simulated by geant4 a beam hitting a target positioned in a vacuum chamber. out of the vacuum chamber I positioned a scintillator to estimate which particles hit the scintillator and the radiation absorbed by it. This is the setup:

I’m looking for electrons, positrons and photons hitting the scinttillator, but there are some runs in which I obtained only electrons and photons but I didn’t obtained positrons (i.e. no positrons hitted the scintillator).
Instead of to write in the results that I got 0 positrons I would set an UL
Is there some UL for this case that I can easily estimate by ROOT? I simulated 3*10^8 primary electrons.

I know this method (FC) to calculate the selected events Ns on a total of events NTot

root  TEfficiency::FeldmanCousins(300000000, 0,0.90,true)

and this one to estimate the UL on an integer number of events

root  TFeldmanCousins fc;
root  fc.SetCL(0.90);
root  double upperLimitFL = fc.CalculateUpperLimit(0, 0.5);
root  upperLimitFL

Maybe the second one could be ok for me?

Maybe @moneta would know the best method?

ROOT Version: Not Provided
Platform: Not Provided
Compiler: Not Provided

Hi,
I think the method you are using are not correct for your problem. From my understanding the observation of positrons is a Poisson process without any background. TFeldmanCousins should be used for cases where there is a background, and TEfficiency for binomial processes.
You can use the classic frequentist Poisson confidence interval, see for example table 40.3 of https://pdg.lbl.gov/2021/reviews/rpp2020-rev-statistics.pdf.
ROOT implements the 1-sigma (66%) or 2-sigma (95%) intervals for the Poisson in TH1::GetBInErrorUp/Low after setting TH1::SetBinErrorOption(TH1::kPoisson) or ``TH1::SetBinErrorOption(TH1::kPoisson2)`.
See ROOT: hist/hist/src/TH1.cxx Source File.

Best regards
Lorenzo

Hello @moneta thank you for your help

Double_t TH1::GetBinErrorUp(Int_t bin) const
{
if (fBinStatErrOpt == kNormal) return GetBinError(bin);
// in case of weighted histogram check if it is really weighted
if (fSumw2.fN && fTsumw != fTsumw2) return GetBinError(bin);
if (bin < 0) bin = 0;
if (bin >= fNcells) bin = fNcells-1;
if (fBuffer) ((TH1*)this)->BufferEmpty();

Double_t alpha = 1.- 0.682689492;
if (fBinStatErrOpt == kPoisson2) alpha = 0.05;

Double_t c = RetrieveBinContent(bin);
Int_t n = int(c);
if (n < 0) {
Warning("GetBinErrorUp","Histogram has negative bin content-force usage to normal errors");
((TH1*)this)->fBinStatErrOpt = kNormal;
return GetBinError(bin);
}

// for N==0 return an upper limit at 0.68 or (1-alpha)/2 ?
// decide to return always (1-alpha)/2 upper interval
//if (n == 0) return ROOT::Math::gamma_quantile_c(alpha,n+1,1);
return ROOT::Math::gamma_quantile_c( alpha/2, n+1, 1) - c;
}

in my macro…but I didn’t understand how to get the value of the UL

scintillator.cpp (4.2 KB)

I runned a simulation with few events if it can help WeTransfer - Send Large Files & Share Photos Online - Up to 2GB Free

Hi,
For the upper limit when you observe 0 events, you just need to call:

int n = 0;  // number of observed events
double alpha = 1. - 0.9; // 0.9 is the CL
ROOT::Math::gamma_quantile_c( alpha, n+1, 1);

For varying n you can reproduce the numbers in the table 40.3 of the documents linked above.

Lorenzo

Hello @moneta thank you, I wrote the code in a root session,

root   int n =0;
root  double alpha=1.-0.9;
root  ROOT::Math::gamma_quantile_c(alpha,n+1,1);
root 

But i didn’t get any value

If you remove ; at the end of the line, ROOT will print the result.
Otherwise do:

double value =   ROOT::Math::gamma_quantile_c(alpha,n+1,1);
std::cout << value << std::endl;

Thank you @moneta!

This topic was automatically closed 14 days after the last reply. New replies are no longer allowed.