flyopunta/src/src_main/random.cpp

/***************************************************************************
                          grandom.cpp  -  description
                             -------------------
    begin                : Sun Nov 26 2000
    copyright            : (C) 2000 by Giuseppe m Pierazzini
    email                : pierazzini@unipi.it
 ***************************************************************************
 *                                                                         *
 *   NA48  simulation program.                                             *
 *                                                                         *
 ***************************************************************************/

// grandom  random generators    GMP 1.10.2000
#include "parm.h"
#include "evento.h"
#include "particella.h"
#include "flyorea.h"

int Rdm1=314,Rdm2=625,idum1=0,idum2=0,iy_rand=0,iv_rand[32];
//------------------------------------------------------------
//-------------------------R A N D O M -----------------------
#define IM1 2147483563
#define IM2 2147483399
#define AM (1.0/IM1)
#define IMM1 (IM1-1)
#define IA1 40014
#define IA2 40692
#define IQ1 53668
#define IQ2 52774
#define IR1 12211
#define IR2 3791
#define NTAB 32
#define NDIV (1+IMM1/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)


//-------------------------------------------------------------
//             G a u s s r
//-------------------------------------------------------------
//
double  Pgauss ( double s, double c )

{
        int i = 0;
        double a = 0.0;
        while ( i++ < 12 )
                a += Pran();
        return ( s * ( a - 6. ) + c );
}
//
double Pran()
//
// Long period random number generator. Returns a uniform random
// deviate between 0.0 and 1.0 (both exclusive).
// It should be called with  idum1 negative, or with  (idum1,idum2)
// from a previous call.
// The obtained random number is smaller than   RNMX
// pr->prandom contains the last random value.
{
        double prandom;
        int j;
        int k;
        //      int idum1=314,idum2=625;
        //      int iy_rand=0;
        //      int iv_rand[NTAB];
    
        //============================================================
        if ( Rdm1<1 )
        {
                idum1=Rdm1; idum2=Rdm2; Rdm1=100;  //attenzione Rdm1 ridefinito..
                if ( idum1 <= 0 )
                {
                        if ( - ( idum1 ) < 1 ) idum1=1; //avoid idum1=0
                        else idum1 = - ( idum1 );
                        //    idum2=(idum1);  original!!!
                        idum2= ( idum2 );
                        for ( j=NTAB+7;j>=0;j-- )
                        {
                                k= ( idum1 ) /IQ1;
                                idum1=IA1* ( idum1-k*IQ1 )-k*IR1;
                                if ( idum1 < 0 ) idum1 +=IM1;
                                if ( j < NTAB ) iv_rand[j]=idum1;
                        }
                        iy_rand=iv_rand[0];
                }
        }
        //===========================================================

        // We start from here, when not initializing ...
        k= ( idum1 ) /IQ1;
        idum1=IA1* ( idum1-k*IQ1 )-k*IR1;
        if ( idum1 < 0 ) idum1 +=IM1;
        k= ( idum2 ) /IQ2;
        idum2=IA2* ( idum2-k*IQ2 )-k*IR2;
        if ( idum2 < 0 ) idum2 +=IM2;
        j=iy_rand/NDIV;
        iy_rand=iv_rand[j]- ( idum2 );
        iv_rand[j]=idum1;
        if ( iy_rand < 1 ) iy_rand +=IMM1;
        if ( ( prandom=AM*iy_rand ) >RNMX ) prandom=RNMX;
        return prandom;

}
//================================================================
//              S e t _ R a n d o m   X run
//=================================================================
void    Set_Pran()
{
        // correlazione con il run
        Rdm2=evento_.Gen.Run;
        Rdm1=-31415;
        Pran();
}

void Set_PranEv()
{
        // per  correlazione.anche con l' evento
        Rdm2=evento_.Gen.Run;
        Rdm1=-evento_.Gen.Event*512;
        Pran();
}