Re: TH1::ComputeIntegral()

From: Rene Brun <brun_at_pcroot.cern.ch>
Date: Tue, 28 Jun 2005 08:08:01 +0200 (MEST)


Could you clearly explain what is your problem? I do not see anything wrong with this function.

Rene Brun

On Mon, 27 Jun 2005, Jamil
Edgemir wrote:

> Hi,
>
> I have a question about TH1::ComputeIntegral() . Why is it that we
> don't just crash if fIntegral is empty rather than just printing an
> error? Specifically, I'm talking about this part:
>
> if (fIntegral[nbins] == 0 ) {
> Error("ComputeIntegral", "Integral = zero"); return 0;
> }
>
> I guess I don't understand the logic. If you call GetRandom() many
> times as people will often do when sampling from a histogram then this
> error seems to leave things in a state that's sort of hard to debug.
> Or maybe I'll ask the question another way: What is the best way to
> deal with this and locate the problem when I find this in my log files
> where I piped my output?
>
> thanks,
> -jamil
>
> p.s. The function is shown below:
>
> Double_t TH1::ComputeIntegral()
> {
> // Compute integral (cumulative sum of bins)
> // The result stored in fIntegral is used by the GetRandom functions.
> // This function is automatically called by GetRandom when the fIntegral
> // array does not exist or when the number of entries in the histogram
> // has changed since the previous call to GetRandom.
> // The resulting integral is normalized to 1
>
> Int_t bin, binx, biny, binz, ibin;
>
> // delete previously computed integral (if any)
> if (fIntegral) delete [] fIntegral;
>
> // - Allocate space to store the integral and compute integral
> Int_t nbinsx = GetNbinsX();
> Int_t nbinsy = GetNbinsY();
> Int_t nbinsz = GetNbinsZ();
> Int_t nxy = nbinsx*nbinsy;
> Int_t nbins = nxy*nbinsz;
>
> fIntegral = new Double_t[nbins+2];
> ibin = 0;
> fIntegral[ibin] = 0;
> for (binz=1;binz<=nbinsz;binz++) {
> for (biny=1;biny<=nbinsy;biny++) {
> for (binx=1;binx<=nbinsx;binx++) {
> ibin++;
> bin = GetBin(binx, biny, binz);
> fIntegral[ibin] = fIntegral[ibin-1] + GetBinContent(bin);
> }
> }
> }
>
> // - Normalize integral to 1
> if (fIntegral[nbins] == 0 ) {
> Error("ComputeIntegral", "Integral = zero"); return 0;
> }
> for (bin=1;bin<=nbins;bin++) fIntegral[bin] /= fIntegral[nbins];
> fIntegral[nbins+1] = fEntries;
> return fIntegral[nbins];
> }
>
>
>
>
Received on Tue Jun 28 2005 - 08:08:07 MEST

This archive was generated by hypermail 2.2.0 : Tue Jan 02 2007 - 14:45:09 MET