mzfunctions_pp_to_pipi_vandewi.cxx File Reference

#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include "ranlxs.h"
#include "ranlxd.h"
#include "mzparameters.h"
#include "mzfunctions.h"
#include "pipisigmas.h"
#include "mzfunctions_pp_to_pipi_vandewi.h"

Go to the source code of this file.

Functions
void	ranlxs (float r[], int n)
void	ranlxd (double r[], int n)
void	rlxs_init (int level, int seed)
void	rlxd_init (int level, int seed)
void	mz_pp_to_pipi_vandewi_init (int seed, double P, double cos_theta_min, double cos_theta_max)
void	mz_pp_to_pipi_get_p_lattice_site (double P, int &i_min, int &iflag)
void	mz_pp_to_pipi_get_costheta_lattice_site (double cos_theta, int &k_min, int &kflag)
double	mz_pp_to_pipi_vandewi_maximum_sigma (double P)
double	mz_pp_to_pipi_vandewi_sigma_legendre (int i, double cos_theta)
double	mz_pp_to_pipi_sigma (double P, double cos_theta)
void	mz_pp_to_pipi_vandewi_event (double P, double cos_theta_min, double cos_theta_max, double *piplus_p, double *piminus_p)

Function Documentation

void mz_pp_to_pipi_get_costheta_lattice_site	(	double	cos_theta,
		int &	k_min,
		int &	kflag
	)

function mz_pp_to_pipi_get_costheta_lattice_site

INPUT:

cos_theta -> cos(theta*)

OUTPUT:

k_min -> lower site in cos(theta*) lattice for input cos(theta*) kflag -> cos(theta*) extrapolation flag kflag=0: input cos(theta*) at lattice site, i.e. no extrapolation needed iflag=1: input cos(theta*) NOT at lattice site, i.e. extrapolation needed

Definition at line 216 of file mzfunctions_pp_to_pipi_vandewi.cxx.

References a, exit(), and printf().

Referenced by mz_pp_to_pipi_sigma().

{
  double a=(2.0)/(200.0);        //lattice spacing in cos(theta*) 
  double diff;

  //checking arguments range...

  if( !((cos_theta>=-1.0)&&(cos_theta<=1.0)) ){
    printf("ERROR: function mz_pp_to_pipi_get_costheta_lattice_site\n");
    printf("       called with argument cos_theta out of range \n");
    printf("    => cos_theta should be contained in the interval [-1.0,1.0]\n");
    printf("\n\n");
    exit(1);
  }


  //calculate lower cos(theta*) lattice site for input cos_theta
  k_min=int( (cos_theta+1.0)/(a) );

  //decide whether or not extrapolation is needed
  kflag=0;
  diff=cos_theta -1.0*(-1.0 +a*(double)(k_min));
  if(diff>0.0) kflag=1;
}

void mz_pp_to_pipi_get_p_lattice_site	(	double	P,
		int &	i_min,
		int &	iflag
	)

function mz_pp_to_pipi_get_p_lattice_site

INPUT:

P -> pbar momentum lab frame [GeV]

OUTPUT:

i_min -> lower site in momentum lattice for input momentum P iflag -> momentum extrapolation flag iflag=0: input momentum P at lattice site, i.e. no extrapolation needed iflag=1: input momentum P NOT at lattice site, i.e. extrapolation needed

Definition at line 183 of file mzfunctions_pp_to_pipi_vandewi.cxx.

References dim_pz, i, and pz.

Referenced by mz_pp_to_pipi_sigma(), mz_pp_to_pipi_vandewi_init(), and mz_pp_to_pipi_vandewi_maximum_sigma().

{
  for(int i=0; i < dim_pz; i++){
    if(P==pz[i]){ 
      iflag=0;
      i_min=i;
    }
    if(i<(dim_pz-1)){      
      if( (P>pz[i])&&(P<pz[i+1])  ){ 
        iflag=1;
        i_min=i;
      }
    }
  }

}

double mz_pp_to_pipi_sigma	(	double	P,
		double	cos_theta
	)

function mz_pp_to_pipi_sigma

INPUT:

P -> anti-proton momentum (lab frame) [GeV] cos_theta -> cos(theta*)

RETURNS:

value of cross section for the pbar momentum P at the point cos_theta

Definition at line 346 of file mzfunctions_pp_to_pipi_vandewi.cxx.

References a, exit(), mz_linear_extrapolation(), mz_pp_to_pipi_get_costheta_lattice_site(), mz_pp_to_pipi_get_p_lattice_site(), mz_pp_to_pipi_vandewi_sigma_legendre(), printf(), pz, sigma, and sigmas_dt.

Referenced by mz_pp_to_pipi_vandewi_event().

{
  double sigma;

  int i_min, iflag;
  int k_min, kflag;
  double x1,y1,x2,y2;
  double factor=2.77207297; //conversion dsigma/d(cos\theta*) -> dsigma/dt at p=2.43 GeV 

  double sigma_1,sigma_2;
  double a=(2.0)/(200.0);        //lattice spacing in cos(theta*)

  //checking arguments range...

  if( !((cos_theta>=-1.0)&&(cos_theta<=1.0)) ){
    printf("ERROR: function mz_pp_to_pipi_vandewi_sigma\n");
    printf("       called with argument cos_theta out of range \n");
    printf("    => cos_theta should be contained in the interval [-1.0,1.0]\n");
    printf("\n\n");
    exit(1);
  }


  //get momentum-lattice site
  mz_pp_to_pipi_get_p_lattice_site(P, i_min, iflag);

  //get cos(theta*) lattice site
  mz_pp_to_pipi_get_costheta_lattice_site(cos_theta, k_min, kflag);



  //legendre zone
  if( (P>=pz[0])&&(P<=pz[19]) ){
    if(iflag==0){                      //NO p-extrapolation
      sigma=mz_pp_to_pipi_vandewi_sigma_legendre(i_min,cos_theta);
    }
    if(iflag==1){                      //YES p_extrapolation
      x1=pz[i_min];
      y1=mz_pp_to_pipi_vandewi_sigma_legendre(i_min,cos_theta);
      x2=pz[i_min+1];
      y2=mz_pp_to_pipi_vandewi_sigma_legendre(i_min+1,cos_theta);
      sigma=mz_linear_extrapolation(x1,y1,x2,y2,P);
    }
  }


  //legendre/vandewi extrapolation zone
  if( (P>pz[19])&&(P<pz[20]) ){
    if(kflag==0){                 //NO cos(theta*) extrapolation
      x1=pz[19];
      y1=mz_pp_to_pipi_vandewi_sigma_legendre(19,cos_theta)*factor;
      x2=pz[20];
      y2=sigmas_dt[0][k_min];
      sigma=mz_linear_extrapolation(x1,y1,x2,y2,P);
    }
    //
    if(kflag==1){                 //YES cos(theta*) extrapolation
      sigma_1=mz_pp_to_pipi_vandewi_sigma_legendre(19,cos_theta)*factor;
      //
      x1=-1.0+a*(double)(k_min);
      y1=sigmas_dt[0][k_min];
      x2=-1.0+a*(double)(k_min+1);
      y2=sigmas_dt[0][k_min+1];
      sigma_2=mz_linear_extrapolation(x1,y1,x2,y2,cos_theta);
      //
      x1=pz[19];
      y1=sigma_1;
      x2=pz[20];
      y2=sigma_2;
      sigma=mz_linear_extrapolation(x1,y1,x2,y2,P);
    }
  }


  //vandewi zone
  if( (P>=pz[20])&&(P<=pz[38]) ){
    if(  (iflag==0)&&(kflag==0)  ){       //NO p-extrapolation NO cos(theta*)-extrapolation
      sigma=sigmas_dt[i_min-20][k_min];
    }
    //
    if(  (iflag==0)&&(kflag==1)  ){       //NO p-extrapolation YES cos(theta*)-extrapolation
      x1=-1.0+a*(double)(k_min);
      y1=sigmas_dt[i_min-20][k_min];
      x2=-1.0+a*(double)(k_min+1);
      y2=sigmas_dt[i_min-20][k_min+1];
      sigma=mz_linear_extrapolation(x1,y1,x2,y2,cos_theta);
    }
    //
    if(  (iflag==1)&&(kflag==0)  ){       //YES p-extrapolation NO cos(theta*)-extrapolation
      x1=pz[i_min];
      y1=sigmas_dt[i_min-20][k_min];
      x2=pz[i_min+1];
      y2=sigmas_dt[i_min+1-20][k_min];
      sigma=mz_linear_extrapolation(x1,y1,x2,y2,P);
    }
    //
    if(  (iflag==1)&&(kflag==1)  ){       //YES p-extrapolation YES cos(theta*)-extrapolation
      x1=-1.0+a*(double)(k_min);
      y1=sigmas_dt[i_min-20][k_min];
      x2=-1.0+a*(double)(k_min+1);
      y2=sigmas_dt[i_min-20][k_min+1];
      sigma_1=mz_linear_extrapolation(x1,y1,x2,y2,cos_theta);
      //
      x1=-1.0+a*(double)(k_min);
      y1=sigmas_dt[i_min+1-20][k_min];
      x2=-1.0+a*(double)(k_min+1);
      y2=sigmas_dt[i_min+1-20][k_min+1];
      sigma_2=mz_linear_extrapolation(x1,y1,x2,y2,cos_theta);
      //
      x1=pz[i_min];
      y1=sigma_1;
      x2=pz[i_min+1];
      y2=sigma_2;
      sigma=mz_linear_extrapolation(x1,y1,x2,y2,P);
    }
  }

  return sigma;
}

void mz_pp_to_pipi_vandewi_event	(	double	P,
		double	cos_theta_min,
		double	cos_theta_max,
		double *	piplus_p,
		double *	piminus_p
	)

function mz_pp_to_pipi_vandewi_event

INPUT:

P: anti-proton momentum (lab frame) [GeV] cos_theta_min: cos(theta*) minimum; theta*=angle(pi-,pbar) CMS frame cos_theta_max: cos(theta*) maximum

OUTPUT:

piplus_p: 4-momentum pi- lab frame piminus_p: 4-momentum pi+ lab frame

Definition at line 480 of file mzfunctions_pp_to_pipi_vandewi.cxx.

References acos(), cos(), i, max(), mp, mpi, mz_pp_to_pipi_sigma(), mz_pp_to_pipi_vandewi_maximum_sigma(), mzboost(), mzrnd(), mzvmod2(), P, phi, PI, sin(), theta, x, and y.

Referenced by PndPiPiGenerator::ReadEvent().

{
  double M=mp;             //proton mass
  double PI=acos(-1.0);

  double p_p[4];           //lab frame: 4-mom p    (proton target)
  double pbar_p[4];        //lab frame: 4-mom pbar (anti-proton beam)
  double pbarp_p[4];       //lab frame: 4-mom pbarpsystem
  double E;                //lab frame: energy pbar
  double S;                //invariant mass square pbarpsystem
  double sqrt_S;           //invariant mass pbarpsystem

  int flag;
  double x,y;                              
  double max;
  double prob;
  double costheta, theta, phi;

  double piplus_pcm[4];    //CMS frame: 4-mom pi+
  double piminus_pcm[4];   //CMS frame: 4-mom pi-



 //initial state kinematic

  //4-mom of p (lab frame)
  p_p[0]=M;
  p_p[1]=0.;
  p_p[2]=0.;
  p_p[3]=0.;

  //4-mom of pbar (lab frame)
  E=pow((P*P+M*M),0.5);

  pbar_p[0]=E;
  pbar_p[1]=0.;
  pbar_p[2]=0.;
  pbar_p[3]=P;

  //4-mom pbarpsystem (lab frame)
  for(int i=0; i<4; i++){
    pbarp_p[i]=p_p[i]+pbar_p[i];
  }

  //invariant mass pbarpsystem (cms energy)
  S=mzvmod2(4,pbarp_p);
  sqrt_S=pow(S,0.5);


  //get upper bound to cross section for pbar momentum P
  max=mz_pp_to_pipi_vandewi_maximum_sigma(P);


  flag=0;

  while(flag==0){
    x=mzrnd(cos_theta_min,cos_theta_max);   //random cos(theta*) in range [cos_theta_min,cos_theta_max]  
    y=mzrnd(0.0,max);                       //random in range [0.0,upper bound to the cross section]
    prob=mz_pp_to_pipi_sigma(P,x);          //cross section at generated cos(theta*)

    if(y<prob){                             //accept event
      flag=1;

      costheta=x;
      if(costheta > 1.0) costheta=1.0;      //just make sure about range...
      if(costheta < -1.0) costheta=-1.0;
      theta=acos(costheta);                 //theta* pi-
      phi=mzrnd(0.0,2.0*PI);                //phi pi-
 
      //build pi- 4-momentum in CMS frame
      piminus_pcm[0]=sqrt_S/2.0;
      piminus_pcm[1]=pow(S/4.0-mpi*mpi,0.5)*sin(theta)*cos(phi);
      piminus_pcm[2]=pow(S/4.0-mpi*mpi,0.5)*sin(theta)*sin(phi);
      piminus_pcm[3]=pow(S/4.0-mpi*mpi,0.5)*cos(theta);

      //build pi+ 4-momentum in CMS frame
      piplus_pcm[0]=piminus_pcm[0];
      piplus_pcm[1]=-1.0*piminus_pcm[1];
      piplus_pcm[2]=-1.0*piminus_pcm[2];
      piplus_pcm[3]=-1.0*piminus_pcm[3];

      //boost pi- and pi+ 4-momenta to lab frame
      mzboost(1,pbarp_p,piplus_pcm,piplus_p);
      mzboost(1,pbarp_p,piminus_pcm,piminus_p);
    }

  }

}

void mz_pp_to_pipi_vandewi_init	(	int	seed,
		double	P,
		double	cos_theta_min,
		double	cos_theta_max
	)

Author

Manuel Zambrana zambr.nosp@m.ana@.nosp@m.kph.u.nosp@m.ni-m.nosp@m.ainz..nosp@m.de and Dmitry Khaneftd khane.nosp@m.ftd@.nosp@m.kph.u.nosp@m.ni-m.nosp@m.ainz..nosp@m.de

Mainz, May 2011

Definition at line 34 of file mzfunctions_pp_to_pipi_vandewi.cxx.

References exit(), i, mp, mz_pp_to_pipi_get_p_lattice_site(), mzvmod2(), P, printf(), pz, and rlxd_init().

Referenced by PndPiPiGenerator::Init().

{

  double p_a=pz[0];        //momentum range
  double p_b=pz[19];
  double p_c=pz[20];
  double p_d=pz[38];

  double P_PANDA_THRESHOLD=1.5;    //pbar momentum threshold in PANDA

  double M=mp;             //proton mass

  double p_p[4];           //lab frame: 4-mom p    (proton target)
  double pbar_p[4];        //lab frame: 4-mom pbar (anti-proton beam)
  double pbarp_p[4];       //lab frame: 4-mom pbarpsystem
  double E;                //lab frame: energy pbar
  double S;                //invariant mass square pbarpsystem
  double sqrt_S;           //invariant mass pbarpsystem

  int p_flag;              //momentum extrapolation flag
  int i_min;               //lower site in p-lattice for input momentum P

  //checking parameter range...
  printf("\n\n");
  printf("###########################################################\n");

  if( !( (P>=p_a)&&(P<=p_d) ) ){
    printf("ERROR: P (anti-proton momentum) out of range \n");
    printf("    => P should be contained in the interval [%4.2f,%5.2f] GeV\n",p_a,p_d);
    printf("\n\n");
    exit(1);
  }


  if(P<P_PANDA_THRESHOLD){
    printf("WARNING: P (anti-proton momentum) smaller than PANDA threshold %4.2f GeV \n", P_PANDA_THRESHOLD);
    printf("\n\n");
  }


  if( !( (cos_theta_min>=-1.0)&&(cos_theta_max>cos_theta_min)&&(cos_theta_max<=1.0) ) ){
    printf("ERROR: cos(theta*) range NOT allowed \n");
    printf("    => cos(theta*) range [cos(theta*) min, cos(theta*) max]\n");
    printf("       should be contained within [-1.0,1.0]\n");
    printf("\n\n");
    exit(1);
  }


  if( ((P>=p_a)&&(P<p_c))&&( (cos_theta_min<-0.94)||(cos_theta_max>0.94)) ){
    printf("WARNING: in the momentum range [%4.2f,%5.2f) GeV,\n",p_a,p_c);
    printf("         it is recommended to keep cos(theta*) range within the interval [-0.94,0.94]\n");
    printf("         to avoid potentially large extrapolations in the cross section\n");
    printf("\n\n");
  }



  //initial state kinematic

  //4-mom of p (lab frame)
  p_p[0]=M;
  p_p[1]=0.;
  p_p[2]=0.;
  p_p[3]=0.;

  //4-mom of pbar (lab frame)
  E=pow((P*P+M*M),0.5);

  pbar_p[0]=E;
  pbar_p[1]=0.;
  pbar_p[2]=0.;
  pbar_p[3]=P;

  //4-mom pbarpsystem (lab frame)
  for(int i=0; i<4; i++){
    pbarp_p[i]=p_p[i]+pbar_p[i];
  }

  //invariant mass pbarpsystem (cms energy)
  S=mzvmod2(4,pbarp_p);
  sqrt_S=pow(S,0.5);

  //print out kinematics
  printf("Initial state kinematics\n");
  printf("------------------------\n");
  printf("\n");                                
  printf("=> CM energy square:        S = %5.2f GeV^2 \n", S);
  printf("   CM energy:         sqrt(S) = %5.2f GeV \n", sqrt_S);
  printf("\n\n\n\n");



  //determining kinematic range and cross section parametrization...

  printf("Kinematic regime and cross section parametrization\n");
  printf("--------------------------------------------------\n");
  printf("\n");

  if( (P>=p_a)&&(P<=p_b)  ){
    printf("%5.2f < P < %5.2f GeV => low energy regime\n", p_a, p_b);
    printf("cross section parametrization: Legendre polynomial series\n");
    printf("\n\n");
  }

  if( (P>p_b)&&(P<p_c)  ){
    printf("%5.2f < P < %5.2f GeV => GAP energy regime\n", p_b, p_c);
    printf("cross section parametrization: extrapolation Legendre polynomial series / Regge description\n");
    printf("\n\n");
  }

  if( (P>=p_c)&&(P<=pz[24])  ){
    printf("%5.2f < P < %5.2f GeV => transition energy regime\n", p_c, pz[24]);
    printf("cross section parametrization: Regge description\n");
    printf("\n\n");
  }


  if( (P>pz[24])&&(P<=p_d)  ){
    printf("%5.2f < P < %5.2f GeV => high energy regime\n", pz[24],p_d);
    printf("cross section parametrization: Regge description\n");
    printf("\n\n");
  }


  //determining if some cross section extrapolation is needed...

  mz_pp_to_pipi_get_p_lattice_site(P, i_min, p_flag);


  if(p_flag==0){
    printf("P at lattice site: no cross section extrapolation needed \n");
    printf("\n\n\n\n");
  }

  if(p_flag==1){
    printf("P NOT at lattice site: extrapolating cross section between momenta %5.2f and %5.2f GeV\n", pz[i_min], pz[i_min+1]);
    printf("\n\n\n\n");
  }

  printf("###########################################################\n");
  printf("\n\n");

  //initialize RANLUX
  rlxd_init(1,seed);
}

double mz_pp_to_pipi_vandewi_maximum_sigma

(

double

P

)

funtion mz_pp_to_pipi_vandewi_maximum_sigma

INPUT:

P -> pbar input momentum lab frame [GeV]

RETURNS:

cross section maximum for momentum P in full cos(theta*) range

Definition at line 257 of file mzfunctions_pp_to_pipi_vandewi.cxx.

References max(), mz_pp_to_pipi_get_p_lattice_site(), pz, and sigmas_dt.

Referenced by mz_pp_to_pipi_vandewi_event().

{
  double max;
  int i_min, p_flag;

  if( (P>=pz[0])&&(P<pz[1]) )    max=300.0;               // legendre zone
  if( (P>=pz[1])&&(P<pz[2]) )    max=300.0;
  if( (P>=pz[2])&&(P<pz[3]) )    max=150.0;
  if( (P>=pz[3])&&(P<pz[4]) )    max=150.0;
  if( (P>=pz[4])&&(P<pz[5]) )    max=120.0;
  if( (P>=pz[5])&&(P<pz[6]) )    max=160.0;
  if( (P>=pz[6])&&(P<pz[7]) )    max=160.0;
  if( (P>=pz[7])&&(P<pz[8]) )    max=100.0;
  if( (P>=pz[8])&&(P<pz[9]) )    max=130.0;
  if( (P>=pz[9])&&(P<pz[10]) )   max=130.0;
  if( (P>=pz[10])&&(P<pz[11]) )  max=90.0;
  if( (P>=pz[11])&&(P<pz[12]) )  max=90.0;
  if( (P>=pz[12])&&(P<pz[13]) )  max=65.0;
  if( (P>=pz[13])&&(P<pz[14]) )  max=65.0;
  if( (P>=pz[14])&&(P<pz[15]) )  max=65.0;
  if( (P>=pz[15])&&(P<pz[16]) )  max=70.0;
  if( (P>=pz[16])&&(P<pz[17]) )  max=70.0;
  if( (P>=pz[17])&&(P<pz[18]) )  max=65.0;
  if( (P>=pz[18])&&(P<=pz[19]) ) max=23.0;
  //
  if( (P>pz[19])&&(P<pz[20]) ) max=50.0;                   // legendre/vandewi extrapolation zone
  //
  if( (P>=pz[20])&&(P<=pz[38]) ){                           // vandewi zone
       mz_pp_to_pipi_get_p_lattice_site(P, i_min, p_flag);
       max=sigmas_dt[i_min-20][200]*1.10;
  }

  return max;
}

double mz_pp_to_pipi_vandewi_sigma_legendre	(	int	i,
		double	cos_theta
	)

function mz_pp_to_pipi_vandewi_sigma_legendre

INPUT:

i -> momentum lattice site cos_theta -> cos(theta*)

RETURNS:

cross section (legendre polinomial series parametrization) at cos(theta*) for the pbar momentum corresponding to the i-th lattice site

Definition at line 305 of file mzfunctions_pp_to_pipi_vandewi.cxx.

References exit(), fitlegendre, i, mz_legendre_polynomial(), printf(), and sigma.

Referenced by mz_pp_to_pipi_sigma().

{
  double sigma;
  double sum=0.0;

  //checking arguments range...

  if( !((i>=0)&&(i<=19)) ){
    printf("ERROR: function mz_pp_to_pipi_vandewi_sigma_legendre\n");
    printf("       called with argument i out of range \n");
    printf("    => i should be contained in the interval [0,19]\n");
    printf("\n\n");
    exit(1);
  }


  //calculate the cross section using the fit coeficients
  for(int k=0; k<11; k++){
    sum=sum+fitlegendre[i][k]*mz_legendre_polynomial(k,cos_theta);
  }

  sigma=sum;

  return sigma;
}

void ranlxd	(	double	r[],
		int	n
	)

Definition at line 571 of file ranlxd.cxx.

{
   int k;

   if (init==0)
      rlxd_init(1,1);

   for (k=0;k<n;k++) 
   {
      is=next[is];
      if (is==is_old)
         update();
      r[k]=one_bit*((double)(x.num[is+4])+one_bit*(double)(x.num[is]));      
   }
}

void ranlxs	(	float	r[],
		int	n
	)

Author

Manuel Zambrana zambr.nosp@m.ana@.nosp@m.kph.u.nosp@m.ni-m.nosp@m.ainz..nosp@m.de and Dmitry Khaneftd khane.nosp@m.ftd@.nosp@m.kph.u.nosp@m.ni-m.nosp@m.ainz..nosp@m.de

Mainz, May 2011

Definition at line 566 of file ranlxs.cxx.

{
   int k;

   if (init==0)
      rlxs_init(0,1);

   for (k=0;k<n;k++) 
   {
      is=next[is];
      if (is==is_old)
         update();
      r[k]=one_bit*(float)(x.num[is]);      
   }
}

void rlxd_init	(	int	level,
		int	seed
	)

Definition at line 501 of file ranlxd.cxx.

{
   int i,k,l;
   int ibit,jbit,xbit[31];
   int ix,iy;

   if ((INT_MAX<2147483647)||(FLT_RADIX!=2)||(FLT_MANT_DIG<24)||
       (DBL_MANT_DIG<48))
      error(0);

   define_constants();
   
   if (level==1)
      pr=202;
   else if (level==2)
      pr=397;
   else
      error(1);
   
   i=seed;

   for (k=0;k<31;k++) 
   {
      xbit[k]=i%2;
      i/=2;
   }

   if ((seed<=0)||(i!=0))
      error(2);

   ibit=0;
   jbit=18;

   for (i=0;i<4;i++)
   {
      for (k=0;k<24;k++)
      {
         ix=0;

         for (l=0;l<24;l++) 
         {
            iy=xbit[ibit];
            ix=2*ix+iy;
         
            xbit[ibit]=(xbit[ibit]+xbit[jbit])%2;
            ibit=(ibit+1)%31;
            jbit=(jbit+1)%31;
         }

         if ((k%4)!=i)
            ix=16777215-ix;

         x.num[4*k+i]=ix;
      }
   }

   carry.c1=0;
   carry.c2=0;
   carry.c3=0;
   carry.c4=0;

   ir=0;
   jr=7;
   is=91;
   is_old=0;
   prm=pr%12;
   init=1;
}

void rlxs_init	(	int	level,
		int	seed
	)

Definition at line 495 of file ranlxs.cxx.

{
   int i,k,l;
   int ibit,jbit,xbit[31];
   int ix,iy;

   if ((INT_MAX<2147483647)||(FLT_RADIX!=2)||(FLT_MANT_DIG<24))
      error(0);

   define_constants();
   
   if (level==0)
      pr=109;
   else if (level==1)
      pr=202;
   else if (level==2)
      pr=397;
   else
      error(1);
   
   i=seed;

   for (k=0;k<31;k++) 
   {
      xbit[k]=i%2;
      i/=2;
   }

   if ((seed<=0)||(i!=0))
      error(2);

   ibit=0;
   jbit=18;

   for (i=0;i<4;i++)
   {
      for (k=0;k<24;k++)
      {
         ix=0;

         for (l=0;l<24;l++) 
         {
            iy=xbit[ibit];
            ix=2*ix+iy;
         
            xbit[ibit]=(xbit[ibit]+xbit[jbit])%2;
            ibit=(ibit+1)%31;
            jbit=(jbit+1)%31;
         }

         if ((k%4)==i)
            ix=16777215-ix;

         x.num[4*k+i]=ix;
      }
   }

   carry.c1=0;
   carry.c2=0;
   carry.c3=0;
   carry.c4=0;

   ir=0;
   jr=7;
   is=95;
   is_old=0;
   prm=pr%12;
   init=1;
}