Re: TH1::ComputeIntegral()

Hi Jamil,

I am not sure to understand your problem. myHistogram->FillRandom("gaus", 1000) does not call ComputeIntegral or/and GetRandom. A short concrete example would help.

Concerning the error handling, you can set your own error handler. The DefaultErrorHandler prints on stdout/err the error messages. You can implement your own handler. See TError.h,cxx

Rene Brun

On Thu,
30 Jun 2005, Jamil
Edgemir wrote:

> 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];
>>> }
>>>
>>>
>>>
>>>
>>
>
>
>
Received on Wed Jul 06 2005 - 22:50:47 MEST