source: trunk/MagicSoft/Simulation/Corsika/Mmcs/ptram.f

      SUBROUTINE PTRAM( ZN,FACT,PTX,PTY )
 
C-----------------------------------------------------------------------
C  TRA(NSVERSE MOMENTUM FROM) M(IMI EXPERIMENT)
C
C  GENERATION OF TRANSVERSE MOMENTUM FOR PARTICLES IN HDPM  GENERATOR
C  SEE RESULTS FROM UA1/MIMI/96
C  SOME CONSTANTS CHANGED FROM MATHEMATICAL SOLUTION BY DICHOTOMY TO
C  TO TAKE INTO ACCOUNT EFFECT OF REJECTIONS. (TESTIFIED AT VS=630 GEV
C  ONLY)  SEE J.N. CAPDEVIELLE, 24TH ICRC, ROMA 1995
C         AND J.N. CAPDEVIELLE, 9TH ISVHECRI, KARLSRUHE 1996
C  THIS SUBROUTINE IS CALLED FROM PPARAM
C  ARGUMENTS:
C   ZN    = POWER OF TRANSV. MOMENTUM FUNCTION, DEP. ON CENT.RAP.DENSITY
C   FACT  = FACTOR TAKING INTO ACCOUNT PARTICLE SPECIFIC TRANSV.MOMENTUM
C   PTX   = TRANSVERSE MOMENTUM IN X DIRECTION
C   PTY   = TRANSVERSE MOMENTUM IN Y DIRECTION
C
C  DESIGN  : J.N. CAPDEVIELLE CDF PARIS
C  CHANGES : D. HECK    IK3  FZK KARLSRUHE
C-----------------------------------------------------------------------
 
      IMPLICIT DOUBLE PRECISION (A-H,O-Z)
*KEEP,CONST.
      COMMON /CONST/   PI,PI2,OB3,TB3,ENEPER
      DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
*KEEP,RANDPA.
      COMMON /RANDPA/  FAC,U1,U2,RD,NSEQ,ISEED,KNOR
      DOUBLE PRECISION FAC,U1,U2
      REAL             RD(3000)
      INTEGER          ISEED(103,10),NSEQ
      LOGICAL          KNOR
*KEEP,RUNPAR.
      COMMON /RUNPAR/  FIXHEI,THICK0,HILOECM,HILOELB,
     *                 STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
     *                 MONIOU,MDEBUG,NUCNUC,
     *                 CETAPE,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE,
     *                 DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
      COMMON /RUNPAC/  DSN,HOST,USER
      DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
      REAL             STEPFC
      INTEGER          NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
     *                 SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
     *                 N1STTR,MDBASE
      INTEGER          CETAPE
      CHARACTER*79     DSN
      CHARACTER*20     HOST,USER
 
      LOGICAL          DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
     *                 FIX1I,FMUADD,FNKG,FPRINT,FDBASE
     *                ,GHEISH,GHESIG
*KEND.
 
C-----------------------------------------------------------------------
 
CC    IF ( DEBUG ) WRITE(MDEBUG,*) 'PTRAM : ZN=',SNGL(ZN)
 
C  TWO RANDOM NUMBERS
      CALL RMMAR( RD,2,1 )
C  GENERATE <P_T>
      ALFA = -0.05D0
      B = ZN-1.D0
      A = RD(1)/B
      U = 0.D0
      DO 5 J = 1,1000
        F1 =  A * (U+1.D0)**B - 1.D0/B
        IF ( F1 .GE. U ) GOTO 15
        U  = U + 0.05D0
 5    CONTINUE
 15   BETA = U
      ALFA = U - 0.05D0
      IF ( F1-U .EQ. 0.D0 ) GOTO 30
      I = 0
 14   U = 0.5D0 * (ALFA+BETA)
      I = I + 1
      F = A * (U+1.D0)**B - 1.D0/B - U
      IF ( F .EQ. 0.D0 ) GO TO 30
      IF ( ABS(U-BETA) .LE. 1.D-4 ) GOTO 30
      FA = A * (ALFA+1.D0)**B - ALFA - 1.D0/B
      FB = B * (BETA+1.D0)**B - BETA - 1.D0/B
      IF ( F*FA .GE. 0.D0 ) THEN
        ALFA = U
      ELSE
        BETA = U
      ENDIF
      GO TO 14
 
 30   XPT = 0.9154D0 * U
C  2*PI*RANDOM NUMBER FOR ANGLE PHI
      Z   = PI2 * RD(2)
      PTX = XPT * FACT * COS(Z)
      PTY = XPT * FACT * SIN(Z)
 
CC    IF ( DEBUG ) WRITE(MDEBUG,*) 'PTRAM : RD(1,2),XPT=',
CC   *                                      RD(1),RD(2),SNGL(XPT)
 
      RETURN
      END