GEANT321/fluka/betrst.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:19:54  cernlib
   6 * Geant
   7 *
   8 *
   9 #include "geant321/pilot.h"
  10 *CMZ :  3.21/02 29/03/94  15.41.41  by  S.Giani
  11 *-- Author :
  12 *$ CREATE BETRST.FOR
  13 *COPY BETRST
  14 *
  15 *=== betrst ===========================================================*
  16 *
  17       FUNCTION BETRST ( GAM, ETA, X0, X1 )
  18
  19 #include "geant321/dblprc.inc"
  20 #include "geant321/dimpar.inc"
  21 #include "geant321/iounit.inc"
  22 *
  23 *----------------------------------------------------------------------*
  24 *                                                                      *
  25 *     New version:                                                     *
  26 *     Created on 20 february 1991  by    Alfredo Ferrari & Paola Sala  *
  27 *                                                   Infn - Milan       *
  28 *                                                                      *
  29 *     Last change on 20-feb-91     by    Alfredo Ferrari               *
  30 *                                                                      *
  31 *     Sampling from beta distribution in [X0,X1) :                     *
  32 *                                                                      *
  33 *     P(X) = X**(GAM-1.D0)*(1.D0-X)**(ETA-1)*GAMM(ETA+GAM)             *
  34 *          / (GAMM(GAM*GAMM(ETA))                                      *
  35 *                                                                      *
  36 *----------------------------------------------------------------------*
  37 *
  38       REAL RNDM(2)
  39 *
  40 *  +-------------------------------------------------------------------*
  41 *  |
  42       IF ( X1 .LT. X0 ) THEN
  43          WRITE (LUNOUT,*)' Betrst: x1<x0, gam, eta', X1, X0, GAM, ETA
  44          WRITE (LUNERR,*)' Betrst: x1<x0, gam, eta', X1, X0, GAM, ETA
  45          X0 = MIN ( X0, ONEONE-0.05 )
  46          X1 = MIN ( ONEONE, X0/TWOTHI )
  47       END IF
  48 *  |
  49 *  +-------------------------------------------------------------------*
  50       GAMI  = 1.D+00 / GAM
  51       X0GAM = X0**GAM
  52       X1GAM = X1**GAM
  53       X0TO1 = 1.D+00 - X0
  54       CNORM = X1GAM - X0GAM
  55       ETAM1 = ETA - 1.D+00
  56       NTAM1 = NINT (ETA - 1.D+00)
  57 *  +-------------------------------------------------------------------*
  58 *  |
  59       IF ( ETAM1 - NTAM1 .NE. 0.D+00 ) THEN
  60 *  |  +----------------------------------------------------------------*
  61 *  |  |  First sample from X**(gam-1) and then reject according to
  62 *  |  |  (1-X)**(eta-1)
  63   100    CONTINUE
  64             CALL GRNDM(RNDM,2)
  65             BETRST = ( RNDM (1) * CNORM + X0GAM )**GAMI
  66             REJE   = ( ( 1.D+00 - BETRST ) / X0TO1 )**ETAM1
  67          IF ( RNDM (2) .GE. REJE ) GO TO 100
  68 *  |  |
  69 *  |  +----------------------------------------------------------------*
  70 *  |
  71 *  +-------------------------------------------------------------------*
  72 *  |
  73       ELSE
  74 *  |  +----------------------------------------------------------------*
  75 *  |  |  First sample from X**(gam-1) and then reject according to
  76 *  |  |  (1-X)**(eta-1)
  77   200    CONTINUE
  78             CALL GRNDM(RNDM,2)
  79             BETRST = ( RNDM (1) * CNORM + X0GAM )**GAMI
  80             REJE   = ( ( 1.D+00 - BETRST ) / X0TO1 )**NTAM1
  81          IF ( RNDM (2) .GE. REJE ) GO TO 200
  82 *  |  |
  83 *  |  +----------------------------------------------------------------*
  84       END IF
  85 *  |
  86 *  +-------------------------------------------------------------------*
  87       RETURN
  88 *=== End of function betrst ===========================================*
  89       END