Upper limit on a 0 result

[faca87]

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 [0] TEfficiency::FeldmanCousins(300000000, 0,0.90,true)

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

root [1] TFeldmanCousins fc;
root [2] fc.SetCL(0.90);
root [3] double upperLimitFL = fc.CalculateUpperLimit(0, 0.5);
root [4] upperLimitFL

Maybe the second one could be ok for me?

Maybe @moneta would know the best method?

---

Please read tips for efficient and successful posting and posting code

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

---

[moneta]

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

[faca87]

Hello @moneta thank you for your help
I added the function

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

[moneta]

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

[faca87]

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

root [3]  int n =0;
root [4] double alpha=1.-0.9;
root [5] ROOT::Math::gamma_quantile_c(alpha,n+1,1);
root [6]

But i didn’t get any value

[moneta]

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;

[faca87]

Thank you @moneta!