]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 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 |