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

Last change on this file 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: 9.3 KB
Line 
1 SUBROUTINE ISOBAR( E,KIND,AMASS,ASMASS,NOPI )
2
3C-----------------------------------------------------------------------
4C ISOBAR
5C
6C THREE AND FOUR PION DECAY OF HEAVY ISOBAR, DECAY IS PICKED AT
7C RANDOM FROM A UNIFORM DISTRIBUTION WITH EQUAL PROBABILITY
8C CHARGE IS DISTRIBUTED AT RANDOM WITH EQUAL PROBABILITY
9C DECAYS ARE COMPUTED VIA MOMENTA, HAVING UNIFORM DISTRIBUTION
10C UPPER LIMIT OF MOMENTUM DISTRIBUTIONS ARE SPECIFIED BY INPUT DATA
11C ENERGY IS STRICTLY CONSERVED, MOMENTA ONLY ON AVERAGE
12C THIS SUBROUTINE IS CALLED FROM MANY BOX ROUTINES
13C ARGUMENTS:
14C E = AVAILABLE ENERGY IN CM
15C KIND = 1 BACKWARD ISOBAR
16C = 0 FORWARD ISOBAR
17C AMASS = MASS OF HEAVY MESON
18C ASMASS = MASS TO BE LEFT OVER FOR OTHER PARTICLES
19C NOPI = NUMBER OF PIONS TO BE GENERATED
20C-----------------------------------------------------------------------
21
22 IMPLICIT DOUBLE PRECISION (A-H,O-Z)
23*KEEP,BAL.
24 COMMON /BAL/ EBAL
25 DOUBLE PRECISION EBAL(10)
26*KEEP,CONST.
27 COMMON /CONST/ PI,PI2,OB3,TB3,ENEPER
28 DOUBLE PRECISION PI,PI2,OB3,TB3,ENEPER
29*KEEP,ELASTY.
30 COMMON /ELASTY/ ELAST,IELIS,IELHM,IELNU,IELPI
31 DOUBLE PRECISION ELAST
32 INTEGER IELIS(20),IELHM(20),IELNU(20),IELPI(20)
33*KEEP,MULT.
34 COMMON /MULT/ EKINL,MSMM,MULTMA,MULTOT
35 DOUBLE PRECISION EKINL
36 INTEGER MSMM,MULTMA(37,13),MULTOT(37,13)
37*KEEP,PAM.
38 COMMON /PAM/ PAMA,SIGNUM
39 DOUBLE PRECISION PAMA(6000),SIGNUM(6000)
40*KEEP,PARPAR.
41 COMMON /PARPAR/ CURPAR,SECPAR,PRMPAR,OUTPAR,C,
42 * E00,E00PN,PTOT0,PTOT0N,THICKH,ITYPE,LEVL
43 DOUBLE PRECISION CURPAR(14),SECPAR(14),PRMPAR(14),OUTPAR(14),
44 * C(50),E00,E00PN,PTOT0,PTOT0N,THICKH
45 INTEGER ITYPE,LEVL
46*KEEP,PARPAE.
47 DOUBLE PRECISION GAMMA,COSTHE,PHI,H,T,X,Y,CHI,BETA,GCM,ECM
48 EQUIVALENCE (CURPAR(2),GAMMA), (CURPAR(3),COSTHE),
49 * (CURPAR(4), PHI ), (CURPAR(5), H ),
50 * (CURPAR(6), T ), (CURPAR(7), X ),
51 * (CURPAR(8), Y ), (CURPAR(9), CHI ),
52 * (CURPAR(10),BETA), (CURPAR(11),GCM ),
53 * (CURPAR(12),ECM )
54*KEEP,RANDPA.
55 COMMON /RANDPA/ FAC,U1,U2,RD,NSEQ,ISEED,KNOR
56 DOUBLE PRECISION FAC,U1,U2
57 REAL RD(3000)
58 INTEGER ISEED(103,10),NSEQ
59 LOGICAL KNOR
60*KEEP,RUNPAR.
61 COMMON /RUNPAR/ FIXHEI,THICK0,HILOECM,HILOELB,
62 * STEPFC,NRRUN,NSHOW,PATAPE,MONIIN,
63 * MONIOU,MDEBUG,NUCNUC,
64 * CETAPE,
65 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
66 * N1STTR,MDBASE,
67 * DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
68 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
69 * ,GHEISH,GHESIG
70 COMMON /RUNPAC/ DSN,HOST,USER
71 DOUBLE PRECISION FIXHEI,THICK0,HILOECM,HILOELB
72 REAL STEPFC
73 INTEGER NRRUN,NSHOW,PATAPE,MONIIN,MONIOU,MDEBUG,NUCNUC,
74 * SHOWNO,ISHW,NOPART,NRECS,NBLKS,MAXPRT,NDEBDL,
75 * N1STTR,MDBASE
76 INTEGER CETAPE
77 CHARACTER*79 DSN
78 CHARACTER*20 HOST,USER
79
80 LOGICAL DEBDEL,DEBUG,FDECAY,FEGS,FIRSTI,FIXINC,FIXTAR,
81 * FIX1I,FMUADD,FNKG,FPRINT,FDBASE
82 * ,GHEISH,GHESIG
83*KEEP,VKIN.
84 COMMON /VKIN/ BETACM
85 DOUBLE PRECISION BETACM
86*KEND.
87
88C-----------------------------------------------------------------------
89
90 IF ( DEBUG ) WRITE(MDEBUG,201)E,KIND,AMASS,ASMASS,NOPI
91 201 FORMAT(' ISOBAR: E,KIND,AMASS,ASMASS,NOPI=',1P,E10.4,I3,2E10.4,I3)
92
93C RETURN 1 KILLS PARTICLE
94
95 ISET = 1
96 EDI = 0.D0
97 PACC = 0.D0
98 RESTE = AMASS
99
100 IF ( KIND .NE. 0 ) GOTO 100
101
102C-----------------------------------------------------------------------
103C FORWARD ISOBAR
104C PIONS FROM FORWARD ISOBAR
105 IF ( NOPI .EQ. 1 ) THEN
106 INOPI = NOPI
107 A = C(36)
108 W = 0.38D0
109 ELSE
110 W = 1.5D0
111 CALL RMMAR( RD,1,1 )
112 IF ( RD(1) .LE. 0.5 ) THEN
113 INOPI = 3
114 A = C(36)
115 ELSE
116 INOPI = 4
117 A = C(35)
118 ENDIF
119 ENDIF
120
121 GIFCM = (E**2+AMASS**2-ASMASS**2) * 0.5D0 /(E*AMASS)
122 BEIFCM = SQRT(GIFCM**2 - 1.D0) / GIFCM
123 GIFLAB = GCM * GIFCM * (1.D0+BETACM*BEIFCM)
124 BEIFL = SQRT(GIFLAB**2 - 1.D0) / GIFLAB
125
126C PION LOOP FOR FORWARD ISOBAR DECAY
127 13 CONTINUE
128 P = PCL(A,W)
129 PT = PTRANS(DUMMY)
130 GPIIF = SQRT( P**2 / PAMA(8)**2 + 1.D0 )
131 BEPIIF = SQRT(GPIIF**2 - 1.D0) / GPIIF
132 EDI = EDI + SQRT( PAMA(8)**2 + P**2 + PT**2 )
133 RESTE = RESTE - EDI
134
135 CALL RMMAR( RD,3,1 )
136 IF ( RD(1) .GT. 0.5 ) THEN
137 GPILAB = GPIIF*GIFLAB*(1.D0-BEIFL*BEPIIF)
138 PACC = PACC - P
139 ELSE
140 GPILAB = GPIIF*GIFLAB*(1.D0+BEIFL*BEPIIF)
141 PACC = PACC + P
142 ENDIF
143
144C CORRECTIVE ACTION IF PLLAB2 LE 1.E-6
145 PLLAB2 = PAMA(8)**2 *(GPILAB**2 - 1.D0)
146 PLLAB2 = MAX( 1.D-6, PLLAB2 )
147 CTHETA = SQRT( PLLAB2 / (PT**2+PLLAB2) )
148 IF ( CTHETA .GE. C(27) ) THEN
149 CALL ADDANG( COSTHE,PHI, CTHETA,RD(2)*PI2, SECPAR(3),SECPAR(4) )
150 IF ( SECPAR(3) .GE. C(29) ) THEN
151 SECPAR(2) = GPILAB
152C CHARGE ASSIGNMENT
153 IF ( RD(3) .LE. OB3 ) THEN
154 SECPAR(1) = 7.D0
155 ELSEIF ( RD(3) .LE. TB3 ) THEN
156 SECPAR(1) = 8.D0
157 ELSE
158 SECPAR(1) = 9.D0
159 ENDIF
160 DO 5 J = 5,8
161 SECPAR(J) = CURPAR(J)
162 5 CONTINUE
163
164 CALL TSTACK
165 ENDIF
166 ENDIF
167 IF ( RESTE .LE. PAMA(14) .OR. ISET .EQ. INOPI ) GOTO 14
168 ISET = ISET + 1
169 GOTO 13
170
171C NUCLEON FROM FORWARD ISOBAR
172 14 CONTINUE
173 PT = PTRANS(DUMMY)
174 MSMM = MSMM + ISET
175 EPT = SQRT( PAMA(14)**2 + PT**2 )
176 RESTE = RESTE - EPT
177
178 IF ( RESTE .LE. 0.D0 ) THEN
179 GNFLAB = GIFLAB
180 EBAL(1) = EBAL(1) + RESTE
181 ELSE
182 GNIF = (RESTE+PAMA(14)) / PAMA(14)
183 BENIF = SQRT(GNIF**2 - 1.D0) / GNIF
184 IF ( PACC .LE. 0.D0 ) THEN
185 GNFLAB = GIFLAB * GNIF * (1.D0 + BENIF*BEIFL)
186 ELSE
187 GNFLAB = GIFLAB * GNIF * (1.D0 - BENIF*BEIFL)
188 ENDIF
189 ENDIF
190
191 PLLAB2 = PAMA(14)**2 * (GNFLAB**2 - 1.D0)
192 PLLAB2 = MAX( 1.D-6, PLLAB2 )
193 CTHETA = SQRT( PLLAB2 / (PT**2+PLLAB2) )
194 IF ( CTHETA .LT. C(27) ) RETURN
195 CALL RMMAR( RD,2,1 )
196 CALL ADDANG( COSTHE,PHI, CTHETA,RD(1)*PI2, SECPAR(3),SECPAR(4) )
197 IF ( SECPAR(3) .LT. C(29) ) RETURN
198 SECPAR(2) = GNFLAB
199
200 IF ( RD(2) .LT. 0.5 ) THEN
201 IADD = 1
202 ELSE
203 IADD = 0
204 ENDIF
205
206 IF ( ITYPE .EQ. 13 .OR. ITYPE .EQ. 14 ) THEN
207 SECPAR(1) = 14 - IADD
208 ELSE
209 SECPAR(1) = 15 + 10*IADD
210 ENDIF
211
212C CHARGE ASSIGNMENT
213 DO 9 J = 5,8
214 SECPAR(J) = CURPAR(J)
215 9 CONTINUE
216 CALL TSTACK
217
218C FILL HISTOGRAM
219 IN = 1.D0 + SECPAR(2) / GAMMA * 20.D0
220 IN = MIN( IN, 20 )
221 IELIS(IN) = IELIS(IN) + 1
222
223 RETURN
224
225C-----------------------------------------------------------------------
226C BACKWARD ISOBAR
227C PIONS FROM BACKWARD ISOBAR
228 100 CONTINUE
229
230 IF ( NOPI .EQ. 1 ) THEN
231 INOPI = NOPI
232 ELSE
233 CALL RMMAR( RD,1,1 )
234 IF ( RD(1) .LE. 0.5 ) THEN
235 INOPI = 3
236 ELSE
237 INOPI = 4
238 ENDIF
239 ENDIF
240 WORK = MIN( C(11), GAMMA*0.5D0 )
241 MSMM = MSMM + INOPI
242 DO 101 J = 1,INOPI
243 CALL RMMAR( RD,3,1 )
244 GPI = RD(1)*(WORK-1.D0) + 1.D0
245 PT = PTRANS(DUMMY)
246 EDI = EDI+SQRT( PAMA(8)**2+PAMA(8)**2*(GPI**2-1.D0)+PT**2 )
247 RESTE = ASMASS - EDI
248 PLLAB2 = PAMA(8)**2 * (GPI**2 - 1.D0)
249 PLLAB2 = MAX( 1.D-6, PLLAB2 )
250 CTHETA = SQRT( PLLAB2 / (PT**2+PLLAB2) )
251 IF ( CTHETA .GE. C(27) ) THEN
252 CALL ADDANG( COSTHE,PHI, CTHETA,RD(2)*PI2,
253 * SECPAR(3),SECPAR(4) )
254 IF ( SECPAR(3) .GE. C(29) ) THEN
255 SECPAR(2) = GPI
256C CHARGE ASSIGNMENT
257 RR = RD(3)
258 IF ( RR .LE. OB3 ) THEN
259 SECPAR(1) = 7.D0
260 ELSEIF ( RR .LE. TB3 ) THEN
261 SECPAR(1) = 8.D0
262 ELSE
263 SECPAR(1) = 9.D0
264 ENDIF
265 DO 104 I = 5,8
266 SECPAR(I) = CURPAR(I)
267 104 CONTINUE
268
269 CALL TSTACK
270 ENDIF
271 ENDIF
272 IF ( RESTE .LE. PAMA(14) ) GOTO 110
273 101 CONTINUE
274
275C NUCLEON FROM BACKWARD ISOBAR
276
277 110 CONTINUE
278 WORK = MIN( C(10), GAMMA*0.5D0 )
279 CALL RMMAR( RD,3,1 )
280 GNRLAB = RD(1) * (WORK-1.D0) + 1.D0
281 PT = PTRANS(DUMMY)
282 EDI = EDI+SQRT( PAMA(14)**2 + PAMA(14)**2*(GNRLAB**2-1.D0)+PT**2)
283 RESTE = ASMASS - EDI
284 EBAL(2) = EBAL(2) + RESTE
285 PLLAB2 = PAMA(14)**2 * (GNRLAB**2 - 1.D0)
286 PLLAB2 = MAX( 1.D-6, PLLAB2 )
287 CTHETA = SQRT( PLLAB2 / (PT**2+PLLAB2) )
288 IF ( CTHETA .LT. C(27) ) RETURN
289 CALL ADDANG( COSTHE,PHI, CTHETA,RD(2)*PI2, SECPAR(3),SECPAR(4) )
290 IF ( SECPAR(3) .LT. C(29) ) RETURN
291 SECPAR(2) = GNRLAB
292C CHARGE ASSIGNMENT
293 IF ( RD(3) .LT. 0.5 ) THEN
294 SECPAR(1) = 13.D0
295 ELSE
296 SECPAR(1) = 14.D0
297 ENDIF
298 DO 113 J = 5,8
299 SECPAR(J) = CURPAR(J)
300 113 CONTINUE
301 CALL TSTACK
302
303 RETURN
304 END
Note: See TracBrowser for help on using the repository browser.