Re: TH1::ComputeIntegral()

Hi Rene,

I have written code to perform error propagation using Monte Carlo and while doing so I ran into the following problem: When calling the TH1::GetRandom() as shown below;

    TRandom* generator = new TRandom();
    generator->SetSeed(0);

    for (int it = 0; it < counts; it++) {       myOtherHistogram->Fill(myHistogram->GetRandom() );     }

there is a call made to TH1::ComputeIntegral() . That in itself is not a problem unless there is a problem with the TH1. In my case, there was a logical error beforehand and I ended up with myHistogram being constructed with the upper and lower limits the same:

    TH1D* myHistogram = new TH1D("ex2", "ex2", 100, 10, 10);

On such a histogram you can still,

    myHistogram->FillRandom("gaus", 1000);

but the subsequent call to GetRandom() then will then produce an error:

    <TH1D::FillRandom>: Integral = zero

and then continue.

Perhaps this is more of a C++ design question rather than a ROOT issue but why not throw something like a TException rather than just send the Error to the screen?? If I were to call TH1::ComputeIntegral() before TH1::GetRandom() then I could test for a possible zero return value but then the way TH1::GetRandom() is designed means I end up calling it twice when it turns out everything is okay and the GetRandom() should procede.

Is there a reason the call to Error(...) was chosen instead of throwing an exception or somehow separating the two methods? Maybe a TException?

-j

On 6/28/05, Rene Brun <brun_at_pcroot.cern.ch> wrote:
> 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];
> > }
> >
> >
> >
> >
>

-- 
-------------------------------------------------------------
Jamil Egdemir               jamil_at_grad.physics.sunysb.edu
SUNY Stony Brook    http://grad.physics.sunysb.edu/~jamil
Physics C-120
(631) 632-4482
-------------------------------------------------------------