mzfunctions_pp_to_leplep_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 "mzfunctions_pp_to_leplep_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_leplep_vandewi_init (int seed, int particle_flag, double P, double GE_to_GM, double cos_theta_min, double cos_theta_max)
double	mz_pp_to_leplep_vandewi_sigma (int particle_flag, double P, double GE, double GM, double cos_theta)
double	mz_pp_to_leplep_vandewi_sigma_nonorm (int particle_flag, double P, double GE_to_GM, double cos_theta)
void	mz_pp_to_leplep_vandewi_event (int particle_flag, double P, double GE_to_GM, double cos_theta_min, double cos_theta_max, double *lepplus_p, double *lepminus_p)

Function Documentation

void mz_pp_to_leplep_vandewi_event	(	int	particle_flag,
		double	P,
		double	GE_to_GM,
		double	cos_theta_min,
		double	cos_theta_max,
		double *	lepplus_p,
		double *	lepminus_p
	)

function mz_pp_to_leplep_vandewi_event

INPUT: particle_flag: 0=electron, 1=muon, 2=tau P: anti-proton momentum (lab frame) [GeV] GE_to_GM: ratio |G_E|/|G_M| cos_theta_min: cos(theta*) minimum; theta*=angle(lep-,pbar) CMS frame cos_theta_max: cos(theta*) maximum

OUTPUT:

lepplus_p: 4-momentum lep+ lab frame lepminus_p: 4-momentum lep- lab frame

Definition at line 331 of file mzfunctions_pp_to_leplep_vandewi.cxx.

References acos(), cos(), fabs(), i, max(), me, mmu, mp, mtau, mz_pp_to_leplep_vandewi_sigma_nonorm(), mzboost(), mzrnd(), mzvmod2(), P, phi, PI, R, sin(), theta, x, and y.

Referenced by PndLepLepGenerator::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

  double mlepton;          //lepton mass

  double tau;              // q^2/(4M^2)                M: proton mass
  double beta_l;           // sqrt{1-4m^2/s}            m: lepton mass
  double R;                // |G_E|/|G_M|
  double A;                // beta_l^2[1-R^2/tau] / ( 2-beta_l^2 + R^2/tau ) 

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

  double lepplus_pcm[4];    //CMS frame: 4-mom lepton+
  double lepminus_pcm[4];   //CMS frame: 4-mom lepton-


  //set the lepton mass
  if(particle_flag==0) mlepton=me;
  if(particle_flag==1) mlepton=mmu;
  if(particle_flag==2) mlepton=mtau;


 //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 non-normalized cross section
  tau=S/(4.0*M*M);
  beta_l=pow(1.0-4.0*mlepton*mlepton/S,0.5);
  R=GE_to_GM;
  A=( beta_l*beta_l*(1.0-R*R/tau) )/( 2.0-beta_l*beta_l + R*R/tau );
  max=1.0 + fabs(A);


  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_leplep_vandewi_sigma_nonorm(particle_flag,P,GE_to_GM,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* lepton- in [0,PI]
      phi=mzrnd(0.0,2.0*PI);                //phi lepton-    in [0,2.0*PI]
 
      //build lepton- 4-momentum in CMS frame
      lepminus_pcm[0]=sqrt_S/2.0;
      lepminus_pcm[1]=pow(S/4.0-mlepton*mlepton,0.5)*sin(theta)*cos(phi);
      lepminus_pcm[2]=pow(S/4.0-mlepton*mlepton,0.5)*sin(theta)*sin(phi);
      lepminus_pcm[3]=pow(S/4.0-mlepton*mlepton,0.5)*cos(theta);

      //build lepton+ 4-momentum in CMS frame
      lepplus_pcm[0]=lepminus_pcm[0];
      lepplus_pcm[1]=-1.0*lepminus_pcm[1];
      lepplus_pcm[2]=-1.0*lepminus_pcm[2];
      lepplus_pcm[3]=-1.0*lepminus_pcm[3];

      //boost lep+ and lep- 4-momenta to lab frame
      mzboost(1,pbarp_p,lepplus_pcm,lepplus_p);
      mzboost(1,pbarp_p,lepminus_pcm,lepminus_p);
    }

  }


}

void mz_pp_to_leplep_vandewi_init	(	int	seed,
		int	particle_flag,
		double	P,
		double	GE_to_GM,
		double	cos_theta_min,
		double	cos_theta_max
	)

function mz_pp_to_leplep_vandewi_init

INPUT:

seed -> RANLUX seed particle_flag -> 0=electron, 1=muon, 2=tau P -> anti-proton momentum (lab frame) [GeV] GE_to_GM -> ratio |G_E|/|G_M| cos_theta_min -> cos(theta*) minimum; theta*=angle(lep-,pbar) CMS frame cos_theta_max -> cos_theta_max: cos(theta*) maximum

OUTPUT:

without output

DESCRIPTION:

function mz_pp_to_leplep_vandewi_init checks the values of all arguments, killing the job if at least one value lies outside its allowed range. In case of tau production, it determines compatibility with the input antiproton momentum, killing the job if the proton-antiproton center-of-mass energy does not reach the threshold for tau production. Kinematic information is printed out to the log file. RANLUX is initialised.

Definition at line 30 of file mzfunctions_pp_to_leplep_vandewi.cxx.

References exit(), mp, mtau, printf(), and rlxd_init().

Referenced by PndLepLepGenerator::Init().

{
  double P_PANDA_THRESHOLD=1.5;        //pbar momentum threshold in PANDA
  double P_PANDA_UP=15.0;              //pbar momentum upper limit in PANDA

  double M=mp;             //proton mass

  double E;                //lab frame: pbar energy
  double S;                //center-of-mass energy squared of pbarpsystem
  double sqrt_S;           //center-of-mass energy of pbarpsystem
  double S_min;            //threshold for tau+tau- production: center-of-mass energy squared of pbarpsystem 
  double E_min;            //threshold for tau+tau- production: (lab frame) pbar energy 
  double P_min;            //threshold for tau+tau- production: (lab frame) pbar momentum

  //checking parameter range...

  if( (seed<1)||(seed>2147483647) ){   
   printf("ERROR: RANLUX seed out of range: \n");
   printf("    => RANLUX seed should be contained within [1,2^31-1]\n");
   printf("    => job killed!!! \n");
   printf("\n\n");
   exit(1);
  }


  if( !( ( particle_flag==0)||( particle_flag==1)||( particle_flag==2) ) ){
    printf("ERROR: particle flag out of range: \n");
    printf("    => particle flag should be 0 (=electron), 1 (=muon) or 2 (=tau)\n");
    printf("    => job killed!!! \n");
    printf("\n\n");
    exit(1);
  }


  if(P<0.0){
    printf("ERROR: P (anti-proton momentum) out of range: \n");
    printf("    => P should be a positive real number \n");
    printf("    => job killed!!! \n");
    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(P>P_PANDA_UP){
    printf("WARNING: P (anti-proton momentum) larger than PANDA upper limit %4.2f GeV \n", 
           P_PANDA_UP);
    printf("\n\n");
  }


  if(GE_to_GM<0.0){
    printf("ERROR: |G_E|/|G_M|  out of range: \n");
    printf("    => |G_E|/|G_M|  should be a positive real number \n");
    printf("    => job killed!!! \n");
    printf("\n\n");
    exit(1);
  }


  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("    => job killed!!! \n");
    printf("\n\n");
    exit(1);
  }



  //initial state kinematics


  //antiproton energy
  E=pow((P*P+M*M),0.5);


  //center-of-mass energy squared and center-of-mass energy of pbarpsystem
  S=2.0*M*(M+E);
  sqrt_S=pow(S,0.5);


  //kinematic thershold for tau production
  S_min=4.0*mtau*mtau;
  E_min=S_min/(2.0*M)-M;
  P_min=pow((E_min*E_min-M*M),0.5);


  //check if tau production is possible
  if( (particle_flag==2)&&(S<=S_min) ){
   printf("ERROR: center of mass energy below threshold for tau+tau- production \n");
   printf("    => tau+tau- production requires: \n");
   printf("           center-of-mass energy squared   s > 4.0*m_{tau}^2 = %4.2f GeV^2, i.e. \n", S_min);
   printf("           antiproton energy               E > %4.2f GeV, i.e. \n", E_min);
   printf("           antiproton momentum             P > %4.2f GeV \n", P_min);
   printf("    => job killed!!! \n");
   exit(1);
   printf("\n\n");
  }

  //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");


  //initialize RANLUX
  rlxd_init(1,seed);

}

double mz_pp_to_leplep_vandewi_sigma	(	int	particle_flag,
		double	P,
		double	GE,
		double	GM,
		double	cos_theta
	)

function mz_pp_to_leplep_vandewi_sigma

INPUT: particle_flag: 0=electron, 1=muon, 2=tau P: anti-proton momentum (lab frame) [GeV] GE: |G_E| GM: |G_M| cos_theta: cos(theta*), theta*=angle(lep-,pbar) CMS frame

RETURNS:

value of the cross section d / d cos(theta*) [fb] for antiproton momentum P, electric and magnetic form factors GE_to_GM at the point cos_theta

Definition at line 183 of file mzfunctions_pp_to_leplep_vandewi.cxx.

References acos(), alpha_QED, hc2, me, mmu, mp, mtau, PI, and sigma.

                                                      {
  double sigma;

  double hc2_units=hc2*1000000000000.0;    //(hbar*c)^2 [fb*GeV^2]

  double PI=acos(-1.0);

  double M=mp;             //proton mass
  double mlepton;          //lepton mass
 
  double E;                //lab frame: pbar energy
  double S;                //center-of-mass-energy squared of pbarpsystem

  double tau;              // q^2/(4M^2)                M: proton mass
  double beta_l;           // sqrt{1-4m^2/s}            m: lepton mass
  double beta_p;           // sqrt{1-4M^2/s}

  double A_1, A_2, A_3, A_4;  //kinematic factors 

  double S_min;   //threshold for tau+tau- production: center-of-mass energy squared of pbarpsystem

  //set the lepton mass
  if(particle_flag==0) mlepton=me;
  if(particle_flag==1) mlepton=mmu;
  if(particle_flag==2) mlepton=mtau;

  //antiproton energy
  E=pow((P*P+M*M),0.5);

  //center-of-mass-energy squared of pbarpsystem
  S=2.0*M*(M+E);

  //kinematic factors
  tau=S/(4.0*M*M);
  beta_l=pow(1.0-4.0*mlepton*mlepton/S,0.5);
  beta_p=pow(1.0-4.0*M*M/S,0.5);

  A_1 = hc2_units*PI*alpha_QED*alpha_QED/2.0;

  A_2 = beta_l/(beta_p*S);

  A_3 = (2.0 -beta_l*beta_l +beta_l*beta_l*cos_theta*cos_theta)*GM*GM;

  A_4 = (1.0/tau)*(1.0 -beta_l*beta_l*cos_theta*cos_theta)*GE*GE;

  //calculate cross section
  sigma = A_1*A_2*(A_3+A_4);

  //return null value for tau production below threshold
  S_min=4.0*mtau*mtau;
  if( (particle_flag==2)&&(S<=S_min) ) sigma=0.0;

  return sigma;
}

double mz_pp_to_leplep_vandewi_sigma_nonorm	(	int	particle_flag,
		double	P,
		double	GE_to_GM,
		double	cos_theta
	)

function mz_pp_to_leplep_vandewi_sigma_nonorm

INPUT: particle_flag: 0=electron, 1=muon, 2=tau P: anti-proton momentum (lab frame) [GeV] GE_to_GM: ratio |G_E|/|G_M| cos_theta: cos(theta*), theta*=angle(lep-,pbar) CMS frame

RETURNS:

value of the non-normalized cross section d / d cos(theta*) for antiproton momentum P and and ratio GE_to_GM at the point cos_theta

the non-normalized cross section is

d sigma / d cos(theta*) ~ 1 + A*cos(theta*), with

A = beta_l^2[1-R^2/tau] / ( 2-beta_l^2 + R^2/tau ), where

beta_l = sqrt{1-4m^2/s} m: lepton mass tau = q^2/(4M^2) M: proton mass

Definition at line 263 of file mzfunctions_pp_to_leplep_vandewi.cxx.

References me, mmu, mp, mtau, R, and sigma.

Referenced by mz_pp_to_leplep_vandewi_event().

                                                             {
  double sigma;

  double M=mp;             //proton mass
  double mlepton;          //lepton mass
 
  double E;                //lab frame: pbar energy
  double S;                //center-of-mass-energy squared of pbarpsystem

  double tau;              // q^2/(4M^2)                M: proton mass
  double beta_l;           // sqrt{1-4m^2/s}            m: lepton mass
  double R;                // |G_E|/|G_M|
  double A;                // beta_l^2[1-R^2/tau] / ( 2-beta_l^2 + R^2/tau )


  //set the lepton mass
  if(particle_flag==0) mlepton=me;
  if(particle_flag==1) mlepton=mmu;
  if(particle_flag==2) mlepton=mtau;

  //antiproton energy
  E=pow((P*P+M*M),0.5);

  //center-of-mass-energy squared of pbarpsystem
  S=2.0*M*(M+E);

  //kinematic factors
  tau=S/(4.0*M*M);
  beta_l=pow(1.0-4.0*mlepton*mlepton/S,0.5);
  R=GE_to_GM;
  A=( beta_l*beta_l*(1.0-R*R/tau) )/( 2.0-beta_l*beta_l + R*R/tau );  

  //cross section
  sigma = 1.0 +A*cos_theta*cos_theta;

  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.

Referenced by mz_pp_to_leplep_vandewi_init().

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