mzfunctions.cxx File Reference

#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <sstream>
#include <cmath>
#include "mzparameters.h"
#include "mzfunctions.h"
#include "ranlxs.h"
#include "ranlxd.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)
double	mzdelta (int i, int k)
double	mzvscalar (int n, double *p, double *q)
double	mzvmod2 (int n, double *p)
void	mzboost (int flag, double *p, double *q, double *q_prime)
double	mzenergy (double m, double p_x, double p_y, double p_z)
double	mzpolar (double *p)
double	mzazimuthal (double *p)
double	mzcospolar (double *p)
double	mzangle (double *p, double *q)
double	mzrnd (double a, double b)
double	mzcomplexmod (double Re_z, double Im_z)
double	mz_E_to_s (double E)
double	mz_legendre_polynomial (int n, double x)
double	mz_linear_extrapolation (double x1, double y1, double x2, double y2, double x)

Function Documentation

double mz_E_to_s

(

double

E

)

function mz_E_to_s

arguments:

E: pbar energy in lab (GeV)

returns:

cm-energy square s of pbar p system

Definition at line 410 of file mzfunctions.cxx.

References exit(), i, mp, mzvmod2(), P, printf(), and s.

                          {
  double s;

  double M=mp;

  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 P;                //lab frame: pz pbar


  //check argument range
  if(E<M){
    printf("ERROR: function mz_E_to_s called with argument E out of range \n");
    printf("    => E should be larger than the proton mass M = %5.3f GeV\n",M);
    printf("\n\n");
    exit(1);
  }


  //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)
  P=pow((E*E-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];
  }

  s=mzvmod2(4,pbarp_p);

  return s;
}

double mz_legendre_polynomial	(	int	n,
		double	x
	)

function mz_legendre_polynomial

INPUT:

n -> polynomial degree, [0,10] x -> independent variable, [-1,1]

returns:

P_{n}(x)

Definition at line 470 of file mzfunctions.cxx.

References exit(), printf(), x, and y.

Referenced by mz_pp_to_pipi_vandewi_sigma_legendre().

                                              {
  double y;

  //checking parameter range...
  if( !( (n>=0)&&(n<=10) ) ){
    printf("ERROR: function  mz_legendre_polynomial called with n out of range \n");
    printf("    => n should be contained in the interval [0,10]\n");
    printf("\n\n");
    exit(1);
  }

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


  if(n==0){
    y=1.0;
  }

  if(n==1){
    y=x;
  }

  if(n==2){
    y=(1.0/2.0)*(3.0*pow(x,2.0) -1.0);
  }

  if(n==3){
    y=(1.0/2.0)*(5.0*pow(x,3.0) -3.0*x);
  }

  if(n==4){
    y=(1.0/8.0)*(35.0*pow(x,4.0) -30.0*pow(x,2.0) +3.0);
  }

  if(n==5){
    y=(1.0/8.0)*(63.0*pow(x,5.0) -70.0*pow(x,3.0) +15.0*x);
  }

  if(n==6){
    y=(1.0/16.0)*(231.0*pow(x,6.0) -315.0*pow(x,4.0) 
                 +105.0*pow(x,2.0) -5.0);
  }

  if(n==7){
    y=(1.0/16.0)*(429.0*pow(x,7.0) -693.0*pow(x,5.0) 
                  +315.0*pow(x,3.0) -35.0*x);
  }

  if(n==8){
    y=(1.0/128.0)*(6435.0*pow(x,8.0) -12012.0*pow(x,6.0)
                   +6930.0*pow(x,4.0) -1260.0*pow(x,2.0) 
                   +35.0);
  }

  if(n==9){
    y=(1.0/128.0)*(12155.0*pow(x,9.0) -25740.0*pow(x,7.0) 
                   +18018.0*pow(x,5.0) -4620.0*pow(x,3.0) 
                   +315.0*x);
  }

  if(n==10){
    y=(1.0/256.0)*(46189.0*pow(x,10.0) -109395.0*pow(x,8.0)
                   +90090.0*pow(x,6.0) -30030*pow(x,4.0)
                   +3465.0*pow(x,2.0) -63.0);
  }

  return y;
}

double mz_linear_extrapolation	(	double	x1,
		double	y1,
		double	x2,
		double	y2,
		double	x
	)

funtion mz_linear_extrapolation

INPUT:

(x1,y1) first point in (x,y)-plane (x2,y2) second point in (x,y)-plane

x independent variable

RETURNS:

y=f(x), with f the linear function atisfying y1=f(x1) and y2=f(x2)

Definition at line 558 of file mzfunctions.cxx.

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

Referenced by mz_pp_to_pipi_sigma().

{
  double y;
  double a,b;

  //checking parameter range...
  if(x2==x1){
    printf("ERROR: function  mz_linear_extrapolation called with x_2=x_1 \n");
    printf("    => argument x_2 should be different from argument x_1 \n");
    printf("\n\n");
    exit(1);
  }

  a = (y1-y2)/(x1-x2);
  b = y1-a*x1;

  y = a*x+b;

  return y;
}

double mzangle	(	double *	p,
		double *	q
	)

function mzangle

input: 4-mom p, 4-mom q output: angle between 3-dim components of p, q in range [0,180]

Definition at line 310 of file mzfunctions.cxx.

References acos(), alpha, mzvmod2(), mzvscalar(), and PI.

{
 double PI=acos(-1.0);

 double alpha;

 double scalar, pmod, qmod, cosalpha;

 pmod=pow(mzvmod2(3,p),0.5);
 qmod=pow(mzvmod2(3,q),0.5);

 scalar=mzvscalar(3,p,q);

 if(pmod==0.0){
    cerr << "WARNING: mzangle called with p[1]=p[2]=p[3]=0.0 " << endl;
    return 9999.;
 }

 if(qmod==0.0){
    cerr << "WARNING: mzangle called with q[1]=q[2]=q[3]=0.0 " << endl;
    return 9999.;
 }


 cosalpha=scalar/(pmod*qmod);

 if(cosalpha > 1.0) cosalpha=1.0;
 if(cosalpha < -1.0) cosalpha=-1.0;

 alpha=acos(cosalpha) * 180.0 / PI;


 return alpha;
}

double mzazimuthal

(

double *

p

)

function mzazimuthal

input: 4-mom p output: azimuthal angle phi in range [0,180]

Definition at line 258 of file mzfunctions.cxx.

References acos(), atan2(), phi, and PI.

{
 double PI=acos(-1.0);

 double phi;

 if ( (p[1]==0.) && (p[2]==0.0) ) {
   cerr << "WARNING: mzazimuthal called with p[1]=p[2]=0.0 " << endl;
   return 9999.;
 }


 phi=atan2(p[2],p[1])*180.0/PI;


 if (phi<0.0) phi = phi + 360.0;

 return phi;
}

void mzboost	(	int	flag,
		double *	p,
		double *	q,
		double *	q_prime
	)

function mzboost

input: flag=0 lab->CM, flag=1 CM->lab p 4-mom of CM particle (M!=0) in lab frame q 4-mom of test particle in lab frame (flag=0) / CM frame (flag=1)

output: q_prime 4-mom of test particle in the CM frame (flag=0) / lab frame (flag=1)

Definition at line 123 of file mzfunctions.cxx.

References i, L, m, mzdelta(), and mzvmod2().

Referenced by mz_pp_to_leplep_vandewi_event(), and mz_pp_to_pipi_vandewi_event().

{
  int i,k;

  double m;
  double gamma;
  double p_uni[4];

  double L[4][4];
  double LI[4][4];

  double mztemp;



  //computing matrix L(p)

  m=pow(mzvmod2(4,p),0.5);
  gamma=p[0]/m;


  for (i=1; i<=3; i++) {
    p_uni[i]=p[i]/pow(mzvmod2(3,p),0.5);
  }


  L[0][0]=gamma;


  for (i=1; i<=3; i++) {
    L[0][i] = p_uni[i] * pow((gamma*gamma - 1.0),0.5);
    L[i][0] = L[0][i];
  }


  for (i=1; i<=3; i++) {
    for (k=1; k<=3; k++) {
      L[i][k] = mzdelta(i,k) + (gamma -1.0)*p_uni[i]*p_uni[k];
    }
  }




  //inverting matrix L(p)

  LI[0][0]=L[0][0];

  for (i=1; i<=3; i++) {
    LI[0][i] = -1.0*L[i][0];
    LI[i][0] = -1.0*L[0][i];
  }

  for (i=1; i<=3; i++) {
    for (k=1; k<=3; k++) {
      LI[i][k] = L[k][i];
    }
  }



  //boosting 4-mom q

  for (i=0; i<=3; i++) {
    mztemp = 0.0;
    for (k=0; k<=3; k++) {
      if (flag==0) mztemp = mztemp + LI[i][k]*q[k];
      if (flag==1) mztemp = mztemp + L[i][k]*q[k];
    }
    q_prime[i] =mztemp;
  }

}

double mzcomplexmod	(	double	Re_z,
		double	Im_z
	)

function mzcomplexmod

arguments:

Re_z = real part of complex number z Im_z = imaginary part of complex number z

returns:

modulus of complex number z

Definition at line 389 of file mzfunctions.cxx.

References y.

                                             {
  double y;

  y=pow(Re_z*Re_z + Im_z*Im_z,0.5);

  return y;
}

double mzcospolar

(

double *

p

)

function mzcospolar

input: 4-mom p output: cos of polar angle theta

Definition at line 285 of file mzfunctions.cxx.

References mzvmod2().

{
  double costheta;
  
  double pmod;

  pmod=pow(mzvmod2(3,p),0.5);
  
  if(pmod==0.0){
    cerr << "WARNING: mzcospolar called with p[1]=p[2]=p[3]=0.0 " << endl;
    return 9999.;
  }

  costheta  = p[3]/pmod;

  return costheta;
}

double mzdelta	(	int	i,
		int	k
	)

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

function mzdelta

input: i,k integers output: kroneker delta(i,k)

Definition at line 30 of file mzfunctions.cxx.

References d.

Referenced by mzboost().

{
  double d;

  if( k==i )
  {
    d=1.0;
  }
  else{
    d=0.0;
  }

  return d;
}

double mzenergy	(	double	m,
		double	p_x,
		double	p_y,
		double	p_z
	)

input: mass m, 3-mom components px, py, pz of a particle

output: energy of the particle

Definition at line 207 of file mzfunctions.cxx.

{
  double E;

  E=pow(m*m +p_x*p_x +p_y*p_y +p_z*p_z,0.5);

  return E;
}

double mzpolar

(

double *

p

)

function mzpolar

input: 4-mom p output: polar angle theta in range [0,180]

Definition at line 223 of file mzfunctions.cxx.

References acos(), mzvmod2(), PI, and theta.

{
  double PI=acos(-1.0);

  double theta;

  double pmod;
  double costheta;

  pmod=pow(mzvmod2(3,p),0.5);
  
  if(pmod==0.0){
    cerr << "WARNING: mzpolar called with p[1]=p[2]=p[3]=0.0 " << endl;
    return 9999.;
  }

  costheta  = p[3]/pmod;

  if(costheta > 1.0) costheta=1.0;
  if(costheta < -1.0) costheta=-1.0;

  theta = acos(costheta) * 180.0 / PI ;


  return theta;
}

double mzrnd	(	double	a,
		double	b
	)

function mzrnd

arguments:

a = inf[a,b] b = sup[a,b]

returns:

random number in interval [a,b]

Definition at line 353 of file mzfunctions.cxx.

References a, exit(), printf(), ranlxd(), x, and y.

Referenced by mz_pp_to_leplep_vandewi_event(), and mz_pp_to_pipi_vandewi_event().

                                {
  double y;

  int dim=1;
  double rvec[dim];
  double x;

  //check arguments
  if(b<=a){
    printf("ERROR: function mzrnd called with arguments b<=a \n");
    printf("    => argument b should be larger than argument a \n");
    printf("\n\n");
    exit(1);
  }

  ranlxd(rvec,dim);  //generate random number in interval [0,1]
  x=rvec[0];
  y=(b-a)*x + a; 

  return y;
}

double mzvmod2	(	int	n,
		double *	p
	)

function mzvmod2

input: int n, 4-vector p

output:

n=2 -> returns p[1]*p[1] + p[2]*p[2] n=3 -> returns p[1]*p[1] + p[2]*p[2] + p[3]*p[3] n=4 -> returns p[0]*p[0] -p[1]*p[1] -p[2]*p[2] -p[3]*p[3]

Definition at line 93 of file mzfunctions.cxx.

References exit(), m2(), mzvscalar(), and printf().

Referenced by mz_E_to_s(), mz_pp_to_leplep_vandewi_event(), mz_pp_to_pipi_vandewi_event(), mz_pp_to_pipi_vandewi_init(), mzangle(), mzboost(), mzcospolar(), and mzpolar().

{
  double m2;


  if( !(n==2||n==3||n==4) ){
    printf("ERROR: function mzvmod2 called with flag n out of range \n");
    printf("    => n should take the value 2,3 or 4 \n");
    printf("\n\n");
    exit(1);
  }


  m2 = mzvscalar(n,p,p);

  return m2;
}

double mzvscalar	(	int	n,
		double *	p,
		double *	q
	)

function mzvscalar

input: int n, 4-vector p, 4-vector q

output:

n=2 -> returns p[1]*q[1] + p[2]*q[2] n=3 -> returns p[1]*q[1] + p[2]*q[2] + p[3]*q[3] n=4 -> returns p[0]*q[0] -p[1]*q[1] -p[2]*q[2] -p[3]*q[3]

Definition at line 53 of file mzfunctions.cxx.

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

Referenced by mzangle(), and mzvmod2().

{
  double s;


  if( !(n==2||n==3||n==4) ){
    printf("ERROR: function mzvscalar called with flag n out of range \n");
    printf("    => n should take the value 2,3 or 4 \n");
    printf("\n\n");
    exit(1);
  }


  if (n==2) {
    s = p[1]*q[1] + p[2]*q[2];
  }

  if (n==3) {
    s = p[1]*q[1] + p[2]*q[2] + p[3]*q[3];
  }

  if (n==4) {
    s = p[0]*q[0] -1.0*p[1]*q[1] -1.0*p[2]*q[2] -1.0*p[3]*q[3];
  }

  return s;
}

void ranlxd	(	double	r[],
		int	n
	)

Definition at line 571 of file ranlxd.cxx.

Referenced by mzrnd().

{
   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.

Referenced by mz_pp_to_leplep_vandewi_init(), mz_pp_to_pipi_vandewi_init(), and ranlxd().

{
   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.

Referenced by ranlxs().

{
   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;
}