TH1::ComputeIntegral()

From: Jamil Edgemir <unclejamil_at_gmail.com>
Date: Mon, 27 Jun 2005 21:40:10 -0700


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 Tue Jun 28 2005 - 06:40:18 MEST

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