NcRandom.cxx

Go to the documentation of this file.

 
#include "NcRandom.h"
#include "Riostream.h"
 
ClassImp(NcRandom); // Class implementation to enable ROOT I/O
 
NcRandom::NcRandom()
{
 
 Int_t seed=53310452; // Default seed
 Start(seed,0,0,0);   // Start the sequence for this seed from scratch
}

NcRandom::NcRandom(Int_t seed,NcTimestamp* ts)
{
 
 Start(seed,0,0,ts); // Start the sequence for this seed from scratch
}

NcRandom::NcRandom(Int_t seed,Int_t cnt1,Int_t cnt2,NcTimestamp* ts)
{
 
 Start(seed,cnt1,cnt2,ts); // Start the sequence from a user defined point
}

NcRandom::~NcRandom()
{
 
 if (fXa) delete [] fXa;
 fXa=0;
 if (fYa) delete [] fYa;
 fYa=0;
 if (fIbins) delete [] fIbins;
 fIbins=0;
}

void NcRandom::Start(Int_t seed,Int_t cnt1,Int_t cnt2,NcTimestamp* ts)
{
 
 if (seed>921350143) seed=53310452;
 
 // Use the provided c.q. current timestamp to create a seed value and start the sequence from scratch
 if (seed<0)
 {
  cnt1=0;
  cnt2=0;
  Int_t jd,sec,ns;
  if (ts)
  {
   ts->GetJD(jd,sec,ns);
  }
  else
  {
   NcTimestamp tx;
   tx.GetJD(jd,sec,ns);
  }
  Int_t jd100=jd%100;
  Int_t sec100=ns/10000000;
  seed=10000*sec+100*sec100+jd100;
 }
 
// Reset the area function
 fNa=0;
 fXa=0;
 fYa=0;
 fIbins=0;
 
// Clipping parameter to prevent overflow of the counting system
 fClip=1000000;
 
// Set the lags for the Fibonacci sequence of the first part
// The sequence is set to F(97,33,*) (see article)
 fI=97;
 fJ=33;
 
// Unpack the seed value and determine i, j, k and l
 fSeed=seed;
 Int_t i,j,k,l;
 Unpack(seed,i,j,k,l);
 
// Reset counters
 fCnt1=0;
 fCnt2=0;
 
// Fill the starting table for the first part of the combination
 Float_t s,t;
 Int_t m;
 for (Int_t ii=0; ii<97; ii++)
 {
  s=0.;
  t=0.5;
 
  for (Int_t jj=1; jj<=24; jj++)
  {
   m=(((i*j)%179)*k)%179;
   i=j;
   j=k;
   k=m;
   l=((53*l)+1)%169;
   if ((l*m)%64 >= 32) s+=t;
   t=0.5*t;
  }
  fU[ii]=s;
 }
 
// Initialise the second part of the combination
 fC=362436./16777216.;
 fCd=7654321./16777216.;
 fCm=16777213./16777216.;
 
// Generate random numbers upto the user selected starting point
// on basis of the counting system
 if (cnt1 > 0) Uniform(cnt1);
 if (cnt2 > 0)
 {
  for (Int_t n=1; n<=cnt2; n++)
  {
   Uniform(fClip);
  }
 }
}

void NcRandom::Unpack(Int_t seed,Int_t& i,Int_t& j,Int_t& k,Int_t& l)
{
 
 if ((seed >= 0) && (seed <= 921350143)) // Check seed value
 {
  Int_t idum=seed;
  Int_t imin2=idum/(176*176*169);
  idum=idum%(176*176*169);
  Int_t jmin2=idum/(176*169);
  idum=idum%(176*169);
  Int_t kmin2=idum/169;
 
  i=imin2+2;
  j=jmin2+2;
  k=kmin2+2;
  l=seed%169;
 }
 else
 {
  cout << " *NcRandom::unpack()* Unallowed seed value encountered."
       << " seed = " << seed << endl;
  cout << " Seed will be set to default value." << endl;
 
  seed=53310452;     // Default seed
  Start(seed,0,0,0); // Start the sequence for this seed from scratch
 }
}

Int_t NcRandom::GetSeed() const
{
 
 return fSeed;
}

Int_t NcRandom::GetCnt1() const
{
 
 return fCnt1;
}

Int_t NcRandom::GetCnt2() const
{
 
 return fCnt2;
}

void NcRandom::Data() const
{
 
 cout << " *NcRandom* seed = " << fSeed
      << " cnt1 = " << fCnt1 << " cnt2 = " << fCnt2 << endl;
}

Float_t NcRandom::Uniform()
{
 
// First part of the combination : F(97,33,*) (see article for explanation)
 Float_t unirnu=fU[fI-1]-fU[fJ-1];
 if (unirnu < 0) unirnu+=1.;
 fU[fI-1]=unirnu;
 fI-=1;
 if (fI == 0) fI=97;
 fJ-=1;
 if (fJ == 0) fJ=97;
 
// Second part of the combination (see article for explanation)
 fC-=fCd;
 if (fC < 0.) fC+=fCm;
 unirnu-=fC;
 if (unirnu < 0.) unirnu+=1.;
 
// Update the counting system to enable sequence continuation
// at an arbitrary starting position.
// Two counters have been introduced to avoid overflow
// fCnt1 : Counter which goes up to fClip
//         and is reset when fClip is reached
// fCnt2 : Counts the number of times fClip has been reached
 fCnt1+=1;
 if (fCnt1 >= fClip)
 {
  fCnt1=0;
  fCnt2+=1;
 }
 
 if (unirnu <= 0.) unirnu=Uniform(); // Exclude 0. from the range
 
 return unirnu;
}

Float_t NcRandom::Uniform(Float_t a,Float_t b)
{
 
 Float_t rmin=a;
 if (a > b) rmin=b;
 
 Float_t rndm=Uniform();
 rndm=rmin+fabs(rndm*(a-b));
 
 return rndm;
}

void NcRandom::Uniform(Float_t* vec,Int_t n,Float_t a,Float_t b)
{
 
// Determine the minimum of a and b
 Float_t rmin=a;
 if (a > b) rmin=b;
 
// First generate random numbers within <0,1>
 if (n > 0) // Check n value
 {
  for (Int_t jvec=0; jvec<n; jvec++)
  {
   // First part of the combination : F(97,33,*)
   Float_t unirnu=fU[fI-1]-fU[fJ-1];
   if (unirnu < 0) unirnu+=1.;
   fU[fI-1]=unirnu;
   fI-=1;
   if (fI == 0) fI=97;
   fJ-=1;
   if (fJ == 0) fJ=97;
 
   // Second part of the combination
   fC-=fCd;
   if (fC < 0.) fC+=fCm;
   unirnu-=fC;
   if (unirnu < 0.) unirnu+=1.;
 
   // Update the counting system to enable sequence continuation
   // at an arbitrary starting position.
   // Two counters have been introduced to avoid overflow
   // fCnt1 : Counter which goes up to fClip
   //         and is reset when fClip is reached
   // fCnt2 : Counts the number of times fClip has been reached
   fCnt1+=1;
   if (fCnt1 >= fClip)
   {
    fCnt1=0;
    fCnt2+=1;
   }
 
   if (unirnu <= 0.) unirnu=Uniform(); // Exclude 0. from the range
 
   // Fill the vector within the selected range
   vec[jvec]=rmin+fabs(unirnu*(a-b));
  }
 }
 else
 {
  cout << " *NcRandom::Uniform* Invalid value n = " << n << endl;
 }
}

void NcRandom::Uniform(Float_t* vec,Int_t n)
{
 
 Uniform(vec,n,0.,1.);
}

void NcRandom::Uniform(Int_t n)
{
 
 if (n > 0) // Check n value
 {
  for (Int_t jvec=0; jvec<n; jvec++)
  {
   // First part of the combination : F(97,33,*)
   Float_t unirnu=fU[fI-1]-fU[fJ-1];
   if (unirnu < 0) unirnu+=1.;
   fU[fI-1]=unirnu;
   fI-=1;
   if (fI == 0) fI=97;
   fJ-=1;
   if (fJ == 0) fJ=97;
 
   // Second part of the combination
   fC-=fCd;
   if (fC < 0.) fC+=fCm;
   unirnu-=fC;
   if (unirnu < 0.) unirnu+=1.;
 
   // Update the counting system to enable sequence continuation
   // at an arbitrary starting position.
   // Two counters have been introduced to avoid overflow
   // fCnt1 : Counter which goes up to fClip
   //         and is reset when fClip is reached
   // fCnt2 : Counts the number of times fClip has been reached
   fCnt1+=1;
   if (fCnt1 >= fClip)
   {
    fCnt1=0;
    fCnt2+=1;
   }
  }
 }
 else
 {
  cout << " *NcRandom::Uniform* Invalid value n = " << n << endl;
 }
}

Float_t NcRandom::Gauss(Float_t mean,Float_t sigma)
{
 
// Generate the two uniform random numbers q1 and q2 in <0,1>
 Float_t q1,q2;
 q1=Uniform();
 q2=Uniform();
 
// Use q1 and q2 to get the gaussian distributed random number
 Float_t pi=acos(-1.);
 Float_t unirng=mean+cos(2.*pi*q2)*sigma*sqrt(-2.*log(q1));
 
 return unirng;
}

Float_t NcRandom::Gauss()
{
 
 return Gauss(0.,1.);
}

void NcRandom::Gauss(Float_t* vec,Int_t n,Float_t mean,Float_t sigma)
{
 
 if (n > 0) // Check n value
 {
  // The vector to hold the q1 and q2 values.
  // Two subsequent q[] values are used for q1 and q2
  // in order to obtain identical random numbers in the vector
  // as when generating n single ones.
  Int_t m=2*n;
  Float_t* q=new Float_t[m];
  Uniform(q,m);
 
  // Fill the vector with gaussian random numbers
  Float_t pi=acos(-1.);
  Float_t q1,q2;
  for (Int_t jvec=0; jvec<n; jvec++)
  {
   q1=q[jvec*2];     // use two subsequent q[] values
   q2=q[(jvec*2)+1];
   vec[jvec]=mean+cos(2.*pi*q2)*sigma*sqrt(-2.*log(q1));
  }
  delete [] q;
 }
 else
 {
  cout << " *NcRandom::Gauss* Invalid value n = " << n << endl;
 }
}

void NcRandom::Gauss(Float_t* vec,Int_t n)
{
 
 Gauss(vec,n,0.,1.);
}

Float_t NcRandom::Poisson(Float_t mean)
{
 
 Float_t unirnp=0.; // Initialise the random number value
 
 if (mean <= 0.) // Check the mean value
 {
  cout << " *NcRandom::Poisson* Invalid value mean = " << mean
       << " Should be positive." << endl;
 }
 else
 {
  if (mean > 80.) // Use gaussian dist. for high mean values to save time
  {
   Float_t grndm=Gauss();
   Float_t rpois=mean+grndm*sqrt(mean);
   Int_t npois=int(rpois);
   if ((rpois-float(npois)) > 0.5) npois++;
   unirnp=float(npois);
  }
  else // Construct a Poisson random number from a uniform one
  {
   Int_t npois=-1;
   Float_t expxm=exp(-mean);
   Float_t poitst=1.;
   while (poitst > expxm)
   {
    Float_t rndm=Uniform();
    npois++;
    poitst=poitst*rndm;
   }
   unirnp=float(npois);
  } // end of check for using Gauss method
 }  // end of mean validity checkn
 return unirnp;
}

void NcRandom::Poisson(Float_t* vec,Int_t n,Float_t mean)
{
 
 if (n <= 0) // Check n value
 {
  cout << " *NcRandom::Poisson* Invalid value n = " << n << endl;
 }
 else
 {
  if (mean <= 0.) // Check the mean value
  {
   cout << " *NcRandom::Poisson* Invalid value mean = " << mean
        << " Should be positive." << endl;
  }
  else
  {
   if (mean > 80.) // Use gaussian dist. for high mean values to save time
   {
    Float_t* grndm=new Float_t[n];
    Gauss(grndm,n);
    Int_t npois;
    Float_t rpois;
    for (Int_t jvec=0; jvec<n; jvec++)
    {
     rpois=mean+grndm[jvec]*sqrt(mean);
     npois=int(rpois);
     if ((rpois-float(npois)) > 0.5) npois++;
     vec[jvec]=float(npois);
    }
    delete [] grndm;
   }
   else // Construct Poisson random numbers from a uniform ones
   {
    Float_t expxm=exp(-mean);
    Int_t npois;
    Float_t poitst;
    for (Int_t jvec=0; jvec<n; jvec++)
    {
     npois=-1;
     poitst=1.;
     while (poitst > expxm)
     {
      Float_t rndm=Uniform();
      npois++;
      poitst=poitst*rndm;
     }
     vec[jvec]=float(npois);
    }
   } // end of check for using Gauss method
  }  // end of mean validity check
 }   // end of n validity check
}

void NcRandom::SetUser(Float_t a,Float_t b,Int_t n,Float_t (*f)(Float_t))
{
 
 fNa=n+1;             // The number of bins for the area function
 fXa=new Float_t[fNa];  // The binned x values of the area function
 fYa=new Float_t[fNa];  // The corresponding summed f(x) values
 fIbins=new Int_t[fNa]; // The bin numbers of the random x candidates
 
 Float_t xmin=a;
 if (a > b) xmin=b;
 Float_t step=fabs(a-b)/float(n);
 
 Float_t x;
 Float_t extreme=0;
 for (Int_t i=0; i<fNa; i++) // Fill bins of the area function
 {
  x=xmin+float(i)*step;
  fXa[i]=x;
  fYa[i]=f(x);
  if (i > 0) fYa[i]+=fYa[i-1];
  if (fabs(fYa[i]) > extreme) extreme=fabs(fYa[i]);
 }
 fYamin=fYa[0]/extreme;
 fYamax=fYa[0]/extreme;
 for (Int_t j=0; j<fNa; j++) // Normalise the area function
 {
  fYa[j]=fYa[j]/extreme;
  if (fYa[j] < fYamin) fYamin=fYa[j];
  if (fYa[j] > fYamax) fYamax=fYa[j];
 }
}

void NcRandom::SetUser(Float_t* x,Float_t* y,Int_t n)
{
 
 fNa=n;               // The number of bins for the area function
 fXa=new Float_t[fNa];  // The binned x values of the area function
 fYa=new Float_t[fNa];  // The corresponding summed y values
 fIbins=new Int_t[fNa]; // The bin numbers of the random x candidates
 
// Order input data with increasing x
 fXa[0]=x[0];
 fYa[0]=y[0];
 for (Int_t i=1; i<fNa; i++) // Loop over x[]
 {
  for (Int_t j=0; j<i; j++) // Loop over xa[]
  {
   if (x[i] < fXa[j])
   {
    for (Int_t k=i; k>=j; k--) // Create empty position
    {
     fXa[k+1]=fXa[k];
     fYa[k+1]=fYa[k];
    }
    fXa[j]=x[i]; // Put x[] value in empty position
    fYa[j]=y[i]; // Put y[] value in empty position
    break; // Go for next x[] value
   }
   if (j == i-1) // This x[] value is the largest so far
   {
    fXa[i]=x[i]; // Put x[] value at the end of x[]
    fYa[i]=y[i]; // Put y[] value at the end of y[]
   }
  } // End loop over fXa[]
 }  // End loop over x[]
 
 Float_t extreme=0;
 for (Int_t l=0; l<fNa; l++) // Fill area function
 {
  if (l > 0) fYa[l]+=fYa[l-1];
  if (fabs(fYa[l]) > extreme) extreme=fabs(fYa[l]);
 }
 fYamin=fYa[0]/extreme;
 fYamax=fYa[0]/extreme;
 for (Int_t m=0; m<fNa; m++) // Normalise the area function
 {
  fYa[m]=fYa[m]/extreme;
  if (fYa[m] < fYamin) fYamin=fYa[m];
  if (fYa[m] > fYamax) fYamax=fYa[m];
 }
}

Float_t NcRandom::User()
{
 
 Float_t unirnf=0;
 
 Float_t ra=Uniform(fYamin,fYamax);     // Random value of the area function
 Float_t dmin=100.*fabs(fYamax-fYamin); // Init. the min. distance encountered
 Int_t ncand=0;
 Float_t dist;
 for (Int_t i=0; i<fNa; i++) // Search for fYa[] value(s) closest to ra
 {
  dist=fabs(ra-fYa[i]);
  if (dist <= dmin) // fYa[i] within smallest distance --> extra candidate
  {
   ncand++;
   if (dist < dmin) ncand=1; // Smallest distance so far --> THE candidate
   dmin=dist;
   fIbins[ncand-1]=i+1;
  }
 }
 
 Int_t jbin=0; // The bin number to hold the required x value
 if (ncand == 1) jbin=fIbins[0];
 
 if (ncand > 1) // Multiple x value candidates --> pick one randomly
 {
  Float_t cand=Uniform(1.,float(ncand));
  Int_t jcand=int(cand);
  if ((cand-float(jcand)) > 0.5) jcand++;
  jbin=fIbins[jcand-1];
 }
 
 if (jbin > 0) // Pick randomly a value in this x-bin
 {
  Float_t xlow=fXa[jbin-1];
  if (jbin > 1) xlow=fXa[jbin-2];
  Float_t xup=fXa[jbin-1];
  if (jbin < fNa-1) xup=fXa[jbin];
  unirnf=Uniform(xlow,xup);
 }
 
 if (ncand == 0) cout << " *NcRandom::User* No candidate found." << endl;
 
 return unirnf;
}

void NcRandom::User(Float_t* vec,Int_t n)
{
 
 Float_t unirnf,ra,dmin,dist;
 Int_t ncand,jbin;
 for (Int_t jvec=0; jvec<n; jvec++)
 {
  unirnf=0;
  ra=Uniform(fYamin,fYamax);     // Random value of the area function
  dmin=100.*fabs(fYamax-fYamin); // Init. the min. distance encountered
  ncand=0;
  for (Int_t i=0; i<fNa; i++) // Search for fYa[] value(s) closest to ra
  {
   dist=fabs(ra-fYa[i]);
   if (dist <= dmin) // fYa[i] within smallest distance --> extra candidate
   {
    ncand++;
    if (dist < dmin) ncand=1; // Smallest distance so far --> THE candidate
    dmin=dist;
    fIbins[ncand-1]=i+1;
   }
  }
 
  jbin=0; // The bin number to hold the required x value
  if (ncand == 1) jbin=fIbins[0];
 
  if (ncand > 1) // Multiple x value candidates --> pick one randomly
  {
   Float_t cand=Uniform(1.,float(ncand));
   Int_t jcand=int(cand);
   if ((cand-float(jcand)) > 0.5) jcand++;
   jbin=fIbins[jcand-1];
  }
 
  if (jbin > 0) // Pick randomly a value in this x-bin
  {
   Float_t xlow=fXa[jbin-1];
   if (jbin > 1) xlow=fXa[jbin-2];
   Float_t xup=fXa[jbin-1];
   if (jbin < fNa-1) xup=fXa[jbin];
   unirnf=Uniform(xlow,xup);
  }
 
  if (ncand == 0) cout << " *NcRandom::User* No candidate found." << endl;
 
  vec[jvec]=unirnf;
 
 }
}

Double_t NcRandom::RanBm(Double_t nr,Double_t n,Int_t m,Double_t* p,Double_t* na,Double_t* psia,Double_t psi0,Int_t f,TH1* psih,Int_t ncut,Double_t* nrx)
{
 
 if (nrx) *nrx=-1;
 
 if (nr<0 || fabs(n)<1 || m<1 || nr>1e19 || fabs(n)>1e19) return -1;
 if (n<0 && !na) return -1;
 
 Double_t retval=-1;
 
 Double_t psi=-1;
 Double_t psimin=-1;
 Double_t npsi=0;
 Double_t pk=1./double(m);
 Double_t psum=0;
 Double_t* nk=new Double_t[m];
 Double_t* nsig=0;
 
 NcMath math;
 
 ULong64_t nrep=ULong64_t(nr);
 if (!nrep)
 {
  if (ncut)
  {
   nrep=ULong64_t(1.e19);
  }
  else
  {
   return -1;
  }
 }
 
 ULong64_t ntrial=ULong64_t(fabs(n));
 ULong64_t jrep,jtrial;
 Int_t jm;
 
 // Store the signal configuration, if specified.
 if (na && n<0)
 {
  nsig=new Double_t[m];
  ULong64_t isig;
  for (jm=0; jm<m; jm++)
  {
   isig=0;
   if (na[jm]>0) isig=ULong64_t(na[jm]);
   nsig[jm]=isig;
  }
 }
 
 Float_t rndm;
 for (jrep=0; jrep<nrep; jrep++) // Repetition loop
 {
  if (psia) psia[jrep]=-1;
  for (jm=0; jm<m; jm++)
  {
   nk[jm]=0;
   if (na && jrep==0) na[jm]=0;
  }
 
  for (jtrial=0; jtrial<ntrial; jtrial++) // Random trial loop
  {
   // Selecting randomly (according to p) a certain outcome
   rndm=Uniform();
   psum=0;
   for (jm=0; jm<m; jm++)
   {
    if (p) pk=p[jm];
    psum+=pk;
    if (rndm<psum) // Record this outcome jm
    {
     nk[jm]+=1;
     if (na) na[jm]+=1;
     break;
    }
   } // End of "m" loop
  } // End of "n" loop
 
  // Now we have generated n independent random trials of B_m
 
  // Superimpose a possibly specified signal configuration
  if (nsig) 
  {
   for (jm=0; jm<m; jm++)
   {
    nk[jm]+=nsig[jm];
    na[jm]+=nsig[jm];
   }
  }
 
  // Calculate corresponding statistics
  psi=math.PsiValue(m,nk,p,f);
  if (psi<psimin || psimin<0) psimin=psi;
  if (psi0>=0 && psi>=psi0) npsi+=1;
  if (psia) psia[jrep]=psi;
  if (psih) psih->Fill(psi);
 
  // Check for cutting short the repetiton loop to save CPU time
  if (!ncut) continue;
  if (npsi>=ncut) break;
 
 } // End of repetition loop
 
 if (nrep==1) retval=psi;
 if (nrep>1 && psi0<0) retval=psimin;
 if (nrep>1 && psi0>=0) retval=npsi;
 
 if (nrx)
 {
  *nrx=nrep;
  if (jrep<nrep) *nrx=jrep+1;
 }
 
 if (nk) delete [] nk;
 if (nsig) delete [] nsig;
 
 return retval;
}