From: Jamil Edgemir <unclejamil_at_gmail.com>

Date: Thu, 30 Jun 2005 05:14:01 -0400

Date: Thu, 30 Jun 2005 05:14:01 -0400

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 -------------------------------------------------------------Received on Thu Jun 30 2005 - 11:14:10 MEST

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