source: trunk/MagicSoft/Simulation/Corsika/Mmcs/rnegbi.f@ 343

Last change on this file since 343 was 286, checked in by harald, 25 years ago
This is the start point for further developments of the Magic Monte Carlo Simulation written by Jose Carlos Gonzales. Now it is under control of one CVS repository for the whole collaboration. Everyone should use this CVS repository for further developments.
File size: 2.9 KB
Line 
1 SUBROUTINE RNEGBI( N,XN,ECM )
2
3C-----------------------------------------------------------------------
4C R(ANDOM NUMBER WITH) NEG(ATIVE) BI(NOMIAL DISTRIBUTION)
5C
6C RANDOM NUMBER GENERATOR FOR INTEGER NUMBERS DISTRIBUTED ACCORDING TO
7C A NEGATIVE BINOMIAL DISTRIBUTION WITH PARAMETERS <N> AND K
8C DELIVERS ONLY EVEN NUMBERS AS CHARGE MUST BE CONSERVED
9C THIS SUBROUTINE IS CALLED FROM HDPM
10C ARGUMENTS:
11C XN = <N> AVERAGE VALUE OF N
12C ECM = CENTER OF MASS ENERGY
13C N = RANDOM NUMBER DISTRIBUTED WITH NEG. BIN. DISTR.
14C-----------------------------------------------------------------------
15
16 IMPLICIT NONE
17*KEEP,RANDPA.
18 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
19 DOUBLE PRECISION FAC,U1,U2
20 REAL RD(3000)
21 INTEGER ISEED(103,10),NSEQ
22 LOGICAL KNOR
23*KEEP,RUNPAR.
24 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
25 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
26 * MONIOU,MDEBUG,NUCNUC,
27 * CETAPE,
28 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
29 * N1STTR,MDBASE,
30 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
31 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
32 * ,GHEISH,GHESIG
33 COMMON /RUNPAC/ DSN,HOST,USER
34 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
35 REAL STEPFC
36 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
37 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
38 * N1STTR,MDBASE
39 INTEGER CETAPE
40 CHARACTER*79 DSN
41 CHARACTER*20 HOST,USER
42
43 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
44 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
45 * ,GHEISH,GHESIG
46*KEND.
47
48 DOUBLE PRECISION ECM,P,PN,Q,R,SUM,XI,XK,XN
49 INTEGER N
50C-----------------------------------------------------------------------
51
52CC IF ( DEBUG ) WRITE(MDEBUG,*) 'RNEGBI: XN,ECM=',SNGL(XN),SNGL(ECM)
53
54C PARAMETRIZATION OF PARAMETER K OF NEG.BIN. DISTRIBUTION ACCORDING
55C TO UA5 COLLABORATION, PHYS. LETT. 167B (1986) 476
56 XK = 1.D0 / ( -0.104D0 + 0.058D0 * LOG(ECM) )
57C OTHER PARAMETERS
58 R = XN / XK
59 Q = 1.D0 / (1.D0 + R)
60 P = R * Q
61
62C VALUES FOR N EQUAL 0
63 1 CONTINUE
64 N = 0
65 PN = Q**XK
66 SUM = PN
67C GET UNIFORM RANDOM NUMBER
68 CALL RMMAR( RD,1,1 )
69 IF ( RD(1) .LE. SUM ) GOTO 100
70C COMPARE WITH SUM OVER P(N)
71 DO 2 XI = 1.D0, 1350.D0
72 PN = PN * P * (XK - 1.D0 + XI) / XI
73 SUM = SUM + PN
74 IF ( RD(1) .LE. SUM ) THEN
75 N = XI
76 GOTO 100
77 ENDIF
78 2 CONTINUE
79 N = 1350
80
81 100 CONTINUE
82 IF ( MOD(N,2) .NE. 0 .AND. N .NE. 1 ) GOTO 1
83CC IF (DEBUG) WRITE(MDEBUG,*)'RNEGBI: RD(1),N,<N>=',RD(1),N,SNGL(XN)
84
85 RETURN
86 END
Note: See TracBrowser for help on using the repository browser.