]>
Commit | Line | Data |
---|---|---|
fe4da5cc | 1 | * |
2 | * $Id$ | |
3 | * | |
4 | * $Log$ | |
5 | * Revision 1.1.1.1 1995/10/24 10:19:57 cernlib | |
6 | * Geant | |
7 | * | |
8 | * | |
9 | #include "geant321/pilot.h" | |
10 | *CMZ : 3.21/02 29/03/94 15.41.43 by S.Giani | |
11 | *-- Author : | |
12 | *$ CREATE NUCEVV.FOR | |
13 | *COPY NUCEVV | |
14 | * | |
15 | *=== nucevv ===========================================================* | |
16 | * | |
17 | SUBROUTINE NUCEVV ( NNHAD, KPROJ, PPPROJ, EKPROJ, TXI, TYI, TZI ) | |
18 | ||
19 | #include "geant321/dblprc.inc" | |
20 | #include "geant321/dimpar.inc" | |
21 | #include "geant321/iounit.inc" | |
22 | *----------------------------------------------------------------------* | |
23 | * Nucevt90: new version by A. Ferrari INFN-Milan and CERN-SPS | |
24 | * | |
25 | * Besides code cleaning, some changes have been made to extract | |
26 | * the quantities needed for a correct energy, momentum, electric | |
27 | * charge and baryonic charge conservation. This quantities are put | |
28 | * in /balanc/ common | |
29 | *----------------------------------------------------------------------* | |
30 | * | |
31 | C | |
32 | C************************************************************** | |
33 | C /NUCPAR/ | |
34 | C PARTICLES CREATED IN NUCEVT | |
35 | C PXNU,PYNU,PZNU = X-, Y- AND Z-COMPONENTS OF THE | |
36 | C LAB MOMENTUM OF THE SECONDARY. COORDINATE SYSTEM | |
37 | C DEFINED BY THE PRIMARY PARTICLE. | |
38 | C HEPNU = TOTAL ENERGY IN THE LAB | |
39 | C AMNU = REST MASS | |
40 | C ICHNU = CHARGE | |
41 | C IBARNU = BARYONIC NUMBER | |
42 | C ANNU = NAME OF THE PARTICLE | |
43 | C NRENU = TYPE NUMBER OF THE PARTICLE | |
44 | C************************************************************** | |
45 | C | |
46 | #include "geant321/balanc.inc" | |
47 | #include "geant321/corinc.inc" | |
48 | #include "geant321/depnuc.inc" | |
49 | #include "geant321/hadpar.inc" | |
50 | #include "geant321/inpdat2.inc" | |
51 | #include "geant321/nucdat.inc" | |
52 | #include "geant321/nucpar.inc" | |
53 | #include "geant321/parevt.inc" | |
54 | #include "geant321/part.inc" | |
55 | #include "geant321/resnuc.inc" | |
56 | COMMON /FKPRIN/ IPRI, INIT | |
57 | * The following dimension statement is now obsolete | |
58 | * DIMENSION XQT(20), XDQT(20), IFQT(3,20) | |
59 | REAL RNDM(1) | |
60 | LOGICAL LDSAVE | |
61 | DIMENSION PTHRSH (39) | |
62 | SAVE PTHRSH | |
63 | DATA PTHRSH / 16*5.D+00,2*2.5D+00,5.D+00,3*2.5D+00,8*5.D+00, | |
64 | & 9*2.5D+00 / | |
65 | * | |
66 | * | |
67 | * | |
68 | *----------------------------------------------------------------------* | |
69 | * Eproj = total energy of the projectile | |
70 | * Amproj = mass energy of the projectile | |
71 | * Ekproj = kinetic energy of the projectile | |
72 | * Pproj = momentum of the projectile | |
73 | * Eeproj = original total energy of the projectile | |
74 | * Ppproj = original momentum of the projectile | |
75 | * Ibproj = barionic charge of the projectile | |
76 | * Icproj = electric charge of the projectile | |
77 | * Nnhad = number of particles produced in Nucevt | |
78 | * Nevt = number of high energy collisions | |
79 | * Atemp = local variable with the mass number of the residual | |
80 | * nucleus : it is updated at any interaction and alrea- | |
81 | * dy includes contributions from cascade particles | |
82 | * (even though actually they are produced after the | |
83 | * high energy interactions: since there is no real | |
84 | * correlation except for the one with Nsea, this does | |
85 | * not represent a problem) | |
86 | * Ztemp = same as Atemp for the atomic number | |
87 | *----------------------------------------------------------------------* | |
88 | * | |
89 | TXX = TXI | |
90 | TYY = TYI | |
91 | TZZ = TZI | |
92 | PPROJ = PPPROJ | |
93 | AMPROJ= AM (KPROJ) | |
94 | EPROJ = EKPROJ + AMPROJ | |
95 | EEPROJ= EPROJ | |
96 | C | |
97 | IBPROJ= IBAR (KPROJ) | |
98 | ICPROJ= ICH (KPROJ) | |
99 | * Set Atemp and Ztemp to their initial values | |
100 | ATEMP = IBTAR - IGREYP - IGREYN | |
101 | ZTEMP = ICHTAR - IGREYP | |
102 | NNHAD = 0 | |
103 | NEVT = 0 | |
104 | * Now Nsea is sampled inside Corrin and passed through common /corrinc | |
105 | IF ( NSEA .EQ. 0 ) GO TO 1000 | |
106 | IF ( INIT .EQ. 1 ) WRITE(LUNOUT,1)NNHAD,KPROJ,PPROJ, | |
107 | * AMPROJ,EPROJ,ANUAV,NSEA | |
108 | 1 FORMAT (1X,2I5,4F12.5,I10) | |
109 | * | |
110 | * We have now sampled Nsea, the number of quark-antiquark pairs | |
111 | * interacting besides the valence interaction (total interactions | |
112 | * = Nsea+1 | |
113 | * | |
114 | *or IF (INIT.EQ.1) WRITE(LUNOUT,4)XO,(XSEA(I),XASEA(I),I=1,NSEA) | |
115 | *or 4 FORMAT (10F12.6) | |
116 | * +-------------------------------------------------------------------* | |
117 | * | Sample effective meson nucleus collisions | |
118 | * | Em, im, ekm, pm refer to the meson quantities | |
119 | DO 6 I = 1, NSEA | |
120 | EM = EKPROJ * ( XSEA(I) + XASEA(I) ) | |
121 | * | +----------------------------------------------------------------* | |
122 | * | | Selection of the ddbar or qqbar composition | |
123 | CALL GRNDM(RNDM,1) | |
124 | IF ( RNDM (1) .LT. 0.5D+00 ) THEN | |
125 | IM = 23 | |
126 | * | | | |
127 | * | +----------------------------------------------------------------* | |
128 | * | | | |
129 | ELSE | |
130 | IM = 26 | |
131 | END IF | |
132 | * | | | |
133 | * | +----------------------------------------------------------------* | |
134 | AMPIO = AM (IM) | |
135 | EKM = EM - AMPIO | |
136 | * | +----------------------------------------------------------------* | |
137 | * | | Energy and/or momentum is too low!!! | |
138 | IF ( EKM .LT. ETHSEA .OR. PPROJ .LT. PTHRSH (KPTOIP(KPROJ)) ) | |
139 | & THEN | |
140 | * | | +-------------------------------------------------------------* | |
141 | * | | | | |
142 | IF ( I .LT. NSEA ) THEN | |
143 | II = I + 1 | |
144 | XSEA(II) = XSEA(II) + XSEA(I) | |
145 | XASEA(II)= XASEA(II) + XASEA(I) | |
146 | * | | | | |
147 | * | | +-------------------------------------------------------------* | |
148 | * | | | | |
149 | ELSE | |
150 | END IF | |
151 | * | | | | |
152 | * | | +-------------------------------------------------------------* | |
153 | KT = IJTARG (I) | |
154 | * | | +-------------------------------------------------------------* | |
155 | * | | | Kt is the index of the target nucleon (1=proton,8=neutron) | |
156 | IF ( KT .EQ. 1 ) THEN | |
157 | ZNOW = ZNOW + 1.D+00 | |
158 | KTARP = KTARP - 1 | |
159 | ZNCOLL = ZNCOLL - 1.D+00 | |
160 | * | | | | |
161 | * | | +-------------------------------------------------------------* | |
162 | * | | | | |
163 | ELSE | |
164 | KTARN = KTARN - 1 | |
165 | END IF | |
166 | * | | | | |
167 | * | | +-------------------------------------------------------------* | |
168 | ANOW = ANOW + 1.D+00 | |
169 | ANCOLL = ANCOLL - 1.D+00 | |
170 | GO TO 6 | |
171 | END IF | |
172 | * | | | |
173 | * | +----------------------------------------------------------------* | |
174 | AMPIO2 = AMPIO * AMPIO | |
175 | EPROLD = EKPROJ + AMPROJ | |
176 | PPROLD = PPROJ | |
177 | * | After selection of the meson energy, Ekproj is updated | |
178 | EKPROJ = EKPROJ - EM | |
179 | EPROJ = EKPROJ + AMPROJ | |
180 | * | +----------------------------------------------------------------* | |
181 | * | | All ekproj energy is transferred to em if ekproj < 4 GeV | |
182 | * | | Now modified since it causes troubles to energy conservation by | |
183 | * | | A. Ferrari: force them to have only a valence interaction | |
184 | IF ( EKPROJ .LT. 3.99D0 ) THEN | |
185 | EKPROJ = EKPROJ + EM | |
186 | * | | +-------------------------------------------------------------* | |
187 | * | | | | |
188 | DO 666 IN = I, NSEA | |
189 | KT = IJTARG (IN) | |
190 | * | | | +----------------------------------------------------------* | |
191 | * | | | | Kt is the index of the target nucleon (1=proton,8=neutron) | |
192 | IF ( KT .EQ. 1 ) THEN | |
193 | ZNOW = ZNOW + 1.D0 | |
194 | KTARP = KTARP - 1 | |
195 | * | | | | | |
196 | * | | | +----------------------------------------------------------* | |
197 | * | | | | | |
198 | ELSE | |
199 | KTARN = KTARN - 1 | |
200 | END IF | |
201 | * | | | | | |
202 | * | | | +----------------------------------------------------------* | |
203 | ANOW = ANOW + 1.D0 | |
204 | 666 CONTINUE | |
205 | * | | | | |
206 | * | | +-------------------------------------------------------------* | |
207 | GO TO 66 | |
208 | END IF | |
209 | * | | | |
210 | * | +----------------------------------------------------------------* | |
211 | * | Now the wounded nucleus is selected inside Corevt | |
212 | KT = IJTARG (I) | |
213 | ATEMP = ATEMP - 1.D+00 | |
214 | IF ( KT .EQ. 1 ) ZTEMP = ZTEMP - 1.D+00 | |
215 | * | +---------------------------------------------------------------* | |
216 | * | | | |
217 | IF ( NINT ( ATEMP ) .EQ. 1 ) THEN | |
218 | LLASTN = .TRUE. | |
219 | KTLAST = KT | |
220 | KTINC = IJTARG ( NSEA + 1 ) | |
221 | PM = SQRT ( EM**2 - AMPIO2 ) | |
222 | AMLAST = AM (KTLAST) | |
223 | AMINC = AM (KTINC) | |
224 | UMIN2 = ( AMINC + AMLAST )**2 | |
225 | ITJ = MIN ( KTINC, 2 ) | |
226 | DELTAE = V0WELL (ITJ) - EFRMAV (ITJ) | |
227 | * | | +------------------------------------------------------------* | |
228 | * | | | Selection of the invariant mass of the two last nucleon | |
229 | * | | | system | |
230 | 2020 CONTINUE | |
231 | EKPROJ = EKPROJ - DELTAE | |
232 | EFRM = EFRM - DELTAE | |
233 | ELEFT = ETTOT - EINTR - EUZ - EKPROJ - AMPROJ - EM | |
234 | PPROJ = SQRT ( EKPROJ * ( EKPROJ + 2.D+00 * AMPROJ ) ) | |
235 | PXLEFT = PXTTOT - PXINTR - PUX - ( PPROJ + PM ) * TXX | |
236 | PYLEFT = PYTTOT - PYINTR - PUY - ( PPROJ + PM ) * TYY | |
237 | PZLEFT = PZTTOT - PZINTR - PUZ - ( PPROJ + PM ) * TZZ | |
238 | UMO2 = ELEFT**2 - PXLEFT**2 - PYLEFT**2 - PZLEFT**2 | |
239 | * | | | +----------------------------------------------------------* | |
240 | * | | | | | |
241 | IF ( UMO2 .LE. UMIN2 ) THEN | |
242 | PPDTPL = PXLEFT * TXX + PYLEFT * TYY + PZLEFT * TZZ | |
243 | DELTAE = 0.51D+00 * ( UMIN2 - UMO2 ) / ( ELEFT - | |
244 | & PPDTPL * ( EKPROJ + AMPROJ ) / PPROJ ) | |
245 | * | | | | +-------------------------------------------------------* | |
246 | * | | | | | Skip the sea interaction if deltae < 0 | |
247 | IF ( DELTAE .LT. 0.D+00 ) THEN | |
248 | EKPROJ = EKPROJ + EM | |
249 | ATEMP = ATEMP + 1.D+00 | |
250 | IF ( KT .EQ. 1 ) ZTEMP = ZTEMP + 1.D+00 | |
251 | * | | | | | +----------------------------------------------------* | |
252 | * | | | | | | | |
253 | DO 777 IN = I, NSEA | |
254 | KT = IJTARG (IN) | |
255 | * | | | | | | +-------------------------------------------------* | |
256 | * | | | | | | | Kt is the index of the target nucleon | |
257 | * | | | | | | | (1=proton,8=neutron) | |
258 | IF ( KT .EQ. 1 ) THEN | |
259 | ZNOW = ZNOW + 1.D0 | |
260 | KTARP = KTARP - 1 | |
261 | * | | | | | | | | |
262 | * | | | | | | +-------------------------------------------------* | |
263 | * | | | | | | | | |
264 | ELSE | |
265 | KTARN = KTARN - 1 | |
266 | END IF | |
267 | * | | | | | | | | |
268 | * | | | | | | +-------------------------------------------------* | |
269 | ANOW = ANOW + 1.D0 | |
270 | 777 CONTINUE | |
271 | * | | | | | | | |
272 | * | | | | | +----------------------------------------------------* | |
273 | PPROJ = SQRT ( EKPROJ * ( EKPROJ + 2.D+00 | |
274 | & * AMPROJ ) ) | |
275 | GO TO 66 | |
276 | * | | | | | | |
277 | * | | | | +-------------------------------------------------------* | |
278 | * | | | | | | |
279 | ELSE IF ( DELTAE .GE. EKPROJ ) THEN | |
280 | WRITE ( LUNERR,* )' Nucevv: sea call impossible ', | |
281 | & ' to get Umin2, go to resampling' | |
282 | LRESMP = .TRUE. | |
283 | RETURN | |
284 | END IF | |
285 | * | | | | | | |
286 | * | | | | +-------------------------------------------------------* | |
287 | GO TO 2020 | |
288 | * | | | | |
289 | * | | +-<|--<--<--<--< We need more invariant mass!! | |
290 | END IF | |
291 | * | | | | | |
292 | * | | | +----------------------------------------------------------* | |
293 | * | | | | |
294 | * | | +-------------------------------------------------------------* | |
295 | * | | Now we will divide Eleft and Pleft between the two | |
296 | * | | left nucleons! | |
297 | UMO = SQRT ( UMO2 ) | |
298 | ELACMS = 0.5D+00 * ( UMO2 + AMLAST**2 - AMINC**2 ) | |
299 | & / UMO | |
300 | EINCMS = UMO - ELACMS | |
301 | PCMS = SQRT ( ELACMS**2 - AMLAST**2 ) | |
302 | GAMCM = ELEFT / UMO | |
303 | ETAX = PXLEFT / UMO | |
304 | ETAY = PYLEFT / UMO | |
305 | ETAZ = PZLEFT / UMO | |
306 | CALL RACO ( CXXINC, CYYINC, CZZINC ) | |
307 | PCMSX = PCMS * CXXINC | |
308 | PCMSY = PCMS * CYYINC | |
309 | PCMSZ = PCMS * CZZINC | |
310 | * | | Now go back from the CMS frame to the lab frame!!! | |
311 | * | | First the "inc" nucleon: | |
312 | ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ | |
313 | EKINC = GAMCM * EINCMS + ETAPCM - AMINC | |
314 | PHELP = ETAPCM / (GAMCM + 1.D+00) + EINCMS | |
315 | PXXINC = PCMSX + ETAX * PHELP | |
316 | PYYINC = PCMSY + ETAY * PHELP | |
317 | PZZINC = PCMSZ + ETAZ * PHELP | |
318 | * | | Now the "last" nucleon | |
319 | EKLAST = GAMCM * ELACMS - ETAPCM - AMLAST | |
320 | PHELP = - ETAPCM / (GAMCM + 1.D+00) + ELACMS | |
321 | PXLAST = - PCMSX + ETAX * PHELP | |
322 | PYLAST = - PCMSY + ETAY * PHELP | |
323 | PZLAST = - PCMSZ + ETAZ * PHELP | |
324 | TVEUZ = 0.D+00 | |
325 | TVGRE0 = 0.D+00 | |
326 | EINCT = EINCP + EINCN | |
327 | IF ( EINCT .GT. 0.D+00 ) THEN | |
328 | EINCP = EINCP - EINCP / EINCT * TVGREY | |
329 | EINCN = EINCN - EINCN / EINCT * TVGREY | |
330 | END IF | |
331 | TVGREY = 0.D+00 | |
332 | PSEA = PM + PPROJ - PPROLD | |
333 | LDSAVE = LDIFFR (KPTOIP(IM)) | |
334 | LDIFFR (KPTOIP(IM)) = .FALSE. | |
335 | * | | The last argument for ferevv is set to 0 to signal that it | |
336 | * | | is not a true valence call | |
337 | IVFLAG = 0 | |
338 | CALL FEREVV ( NHAD, IM, KT, PM, EKM, TXX, TYY, TZZ, ATEMP, | |
339 | & ZTEMP, IVFLAG ) | |
340 | LDIFFR (KPTOIP(IM)) = LDSAVE | |
341 | IF ( LRESMP ) RETURN | |
342 | * | | | |
343 | * | +---------------------------------------------------------------* | |
344 | * | | | |
345 | ELSE | |
346 | TXXOLD = TXX | |
347 | TYYOLD = TYY | |
348 | TZZOLD = TZZ | |
349 | PPROJ = SQRT ( EKPROJ * ( EKPROJ + 2.D+00 * AMPROJ ) ) | |
350 | IF ( PPROJ .GT. PPROLD ) THEN | |
351 | PPROJ = PPROLD * EPROJ / EPROLD | |
352 | END IF | |
353 | CALL FERSEA ( NHAD, IM, KT, EKM, TXX, TYY, TZZ, ATEMP, | |
354 | & ZTEMP, KPROJ, EPROJ, PPROJ, EPROLD, PPROLD ) | |
355 | IF ( LRESMP ) RETURN | |
356 | END IF | |
357 | * | | | |
358 | * | +----------------------------------------------------------------* | |
359 | NEVT = NEVT + 1 | |
360 | * | +----------------------------------------------------------------* | |
361 | * | | | |
362 | IF ( NNHAD + NHAD .GT. MXPNUC ) THEN | |
363 | WRITE (LUNOUT,1101) NNHAD | |
364 | 1101 FORMAT(' NHAD IN NUCEVT TOO BIG CHANGE DIMENSION', | |
365 | * ' NNHAD=',I10) | |
366 | STOP | |
367 | END IF | |
368 | * | | | |
369 | * | +----------------------------------------------------------------* | |
370 | ||
371 | * | +----------------------------------------------------------------* | |
372 | * | | Looping over the produced particles!!! | |
373 | DO 7 JJ = 1, NHAD | |
374 | NNHAD = NNHAD + 1 | |
375 | J = NNHAD | |
376 | PXNU(J) = PXH(JJ) | |
377 | PYNU(J) = PYH(JJ) | |
378 | PZNU(J) = PZH(JJ) | |
379 | HEPNU(J) = HEPH(JJ) | |
380 | AMNU(J) = AMH(JJ) | |
381 | IBARNU(J)= IBARH(JJ) | |
382 | ICHNU(J) = ICHH(JJ) | |
383 | NRENU(J) = NREH(JJ) | |
384 | ANNU(J) = ANH(JJ) | |
385 | * | | Updating energy and momentum accumulators | |
386 | PUX = PUX+PXNU(J) | |
387 | PUY = PUY+PYNU(J) | |
388 | PUZ = PUZ+PZNU(J) | |
389 | EUZ = EUZ+HEPNU(J) | |
390 | ICU = ICU+ICHNU(J) | |
391 | IBU = IBU+IBARNU(J) | |
392 | 7 CONTINUE | |
393 | * | | | |
394 | * | +----------------------------------------------------------------* | |
395 | 6 CONTINUE | |
396 | * | | |
397 | * +-------------------------------------------------------------------* | |
398 | 66 CONTINUE | |
399 | * | |
400 | * Computing the actual pproj | |
401 | * | |
402 | EPROJ = EKPROJ + AMPROJ | |
403 | 1000 CONTINUE | |
404 | * +-------------------------------------------------------------------+ | |
405 | * | Now modified because of troubles to energy conservation by | |
406 | * | A. Ferrari: if the incident projectile is a proton or a neutron | |
407 | * | maybe it was stopped in the nucleus, else it is added to the stack, | |
408 | * | other particles are also added to the secondary stack. | |
409 | * | May be is not physical, but better than nothing!!! A different | |
410 | * | solution maybe to let meson and other baryons decay at rest or | |
411 | * | also force them to have only a valence interaction or also | |
412 | * | convert a charged meson plus a residual neutron/proton to a | |
413 | * | proton/neutron and add the extra energy to the sea meson and | |
414 | * | so on ???????????????????? | |
415 | * | Now the condition Pproj .ge. .4 is always fulfilled!!!!! | |
416 | IF (PPROJ .LT. 4.D0) THEN | |
417 | * | +---------------------------------------------------------------* | |
418 | * | | | |
419 | DKPR30 = KPROJ-30 | |
420 | IF ( PPROJ - 1.5D+00 * SIGN ( ONEONE, DKPR30 ) | |
421 | & .LT. 4.D0 ) THEN | |
422 | WRITE ( LUNERR,* )' Nucevt: Pproj < 4 GeV/c!!!', PPROJ, | |
423 | & KPROJ | |
424 | END IF | |
425 | * | | | |
426 | * | +---------------------------------------------------------------* | |
427 | END IF | |
428 | * | | |
429 | * +------------------------------------------------------------------* | |
430 | * Now the wounded nucleus is selected inside Corevt | |
431 | KT = IJTARG ( NSEA + 1 ) | |
432 | * Update the Nsea value | |
433 | NSEA = NEVT | |
434 | *or IF (INIT.EQ.1) WRITE(LUNOUT,162)NEVT,NHAD,IM,KPROJ,KT,PM, | |
435 | *or &EM,EKM,PM,EKPROJ,PPROJ,EPROJ | |
436 | * | |
437 | * Now ferevt performs the interaction: no longer a meson but the actual | |
438 | * projectile | |
439 | * | |
440 | ATEMP = ATEMP - 1.D+00 | |
441 | IF ( KT .EQ. 1 ) ZTEMP = ZTEMP - 1.D+00 | |
442 | * +-------------------------------------------------------------------* | |
443 | * | | |
444 | IF ( NINT ( ATEMP ) .EQ. 1 ) THEN | |
445 | LLASTN = .TRUE. | |
446 | KTLAST = KT | |
447 | IF ( ZNOW .GT. 0.D+00 ) THEN | |
448 | KTINC = 1 | |
449 | ELSE | |
450 | KTINC = 8 | |
451 | END IF | |
452 | AMLAST = AM (KTLAST) | |
453 | AMINC = AM (KTINC) | |
454 | UMIN2 = ( AMINC + AMLAST )**2 | |
455 | ITJ = MIN ( KTINC, 2 ) | |
456 | * | +----------------------------------------------------------------* | |
457 | * | | Selection of the invariant mass of the two last nucleon | |
458 | * | | system | |
459 | 2040 CONTINUE | |
460 | ELEFT = ETTOT - EINTR - EUZ - EKPROJ - AMPROJ | |
461 | PXLEFT = PXTTOT - PXINTR - PUX - PPROJ * TXX | |
462 | PYLEFT = PYTTOT - PYINTR - PUY - PPROJ * TYY | |
463 | PZLEFT = PZTTOT - PZINTR - PUZ - PPROJ * TZZ | |
464 | UMO2 = ELEFT**2 - PXLEFT**2 - PYLEFT**2 - PZLEFT**2 | |
465 | * | | +-------------------------------------------------------------* | |
466 | * | | | | |
467 | IF ( UMO2 .LE. UMIN2 ) THEN | |
468 | PPDTPL = PXLEFT * TXX + PYLEFT * TYY + PZLEFT * TZZ | |
469 | DELTAE = 0.51D+00 * ( UMIN2 - UMO2 ) / ( ELEFT - | |
470 | & PPDTPL * ( EKPROJ + AMPROJ ) / PPROJ ) | |
471 | * | | | +----------------------------------------------------------* | |
472 | * | | | | | |
473 | IF ( DELTAE .GE. EKPROJ .OR. DELTAE .LT. 0.D+00 ) THEN | |
474 | WRITE ( LUNERR,* )' Nucevv: valence call impossible ', | |
475 | & ' to get Umin2, go to resampling' | |
476 | LRESMP = .TRUE. | |
477 | RETURN | |
478 | END IF | |
479 | * | | | | | |
480 | * | | | +----------------------------------------------------------* | |
481 | EKPROJ = EKPROJ - DELTAE | |
482 | PLA = PPROJ | |
483 | PPROJ = SQRT ( EKPROJ * ( EKPROJ + 2.D+00 * AMPROJ ) ) | |
484 | PLA = PPROJ - PLA | |
485 | EFRM = EFRM - DELTAE | |
486 | PXFRM = PXFRM + TXX * PLA | |
487 | PYFRM = PYFRM + TYY * PLA | |
488 | PZFRM = PZFRM + TZZ * PLA | |
489 | GO TO 2040 | |
490 | * | | | |
491 | * | +-<|--<--<--<--< We need more invariant mass!! | |
492 | END IF | |
493 | * | | | | |
494 | * | | +-------------------------------------------------------------* | |
495 | * | | | |
496 | * | +----------------------------------------------------------------* | |
497 | * | Now we will divide Eleft and Pleft between the two | |
498 | * | left nucleons! | |
499 | UMO = SQRT ( UMO2 ) | |
500 | ELACMS = 0.5D+00 * ( UMO2 + AMLAST**2 - AMINC**2 ) / UMO | |
501 | EINCMS = UMO - ELACMS | |
502 | PCMS = SQRT ( ELACMS**2 - AMLAST**2 ) | |
503 | GAMCM = ELEFT / UMO | |
504 | ETAX = PXLEFT / UMO | |
505 | ETAY = PYLEFT / UMO | |
506 | ETAZ = PZLEFT / UMO | |
507 | CALL RACO ( CXXINC, CYYINC, CZZINC ) | |
508 | PCMSX = PCMS * CXXINC | |
509 | PCMSY = PCMS * CYYINC | |
510 | PCMSZ = PCMS * CZZINC | |
511 | * | Now go back from the CMS frame to the lab frame!!! | |
512 | * | First the "inc" nucleon: | |
513 | ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ | |
514 | EKINC = GAMCM * EINCMS + ETAPCM - AMINC | |
515 | PHELP = ETAPCM / (GAMCM + 1.D+00) + EINCMS | |
516 | PXXINC = PCMSX + ETAX * PHELP | |
517 | PYYINC = PCMSY + ETAY * PHELP | |
518 | PZZINC = PCMSZ + ETAZ * PHELP | |
519 | * | Now the "last" nucleon | |
520 | EKLAST = GAMCM * ELACMS - ETAPCM - AMLAST | |
521 | PHELP = - ETAPCM / (GAMCM + 1.D+00) + ELACMS | |
522 | PXLAST = - PCMSX + ETAX * PHELP | |
523 | PYLAST = - PCMSY + ETAY * PHELP | |
524 | PZLAST = - PCMSZ + ETAZ * PHELP | |
525 | TVEUZ = 0.D+00 | |
526 | TVGRE0 = 0.D+00 | |
527 | EINCT = EINCP + EINCN | |
528 | IF ( EINCT .GT. 0.D+00 ) THEN | |
529 | EINCP = EINCP - EINCP / EINCT * TVGREY | |
530 | EINCN = EINCN - EINCN / EINCT * TVGREY | |
531 | END IF | |
532 | TVGREY = 0.D+00 | |
533 | END IF | |
534 | * | | |
535 | * +-------------------------------------------------------------------* | |
536 | * The last argument for ferevv is set to 1 to signal that it | |
537 | * is a true valence call | |
538 | IVFLAG = 1 | |
539 | CALL FEREVV ( NHAD, KPROJ, KT, PPROJ, EKPROJ, TXX, TYY, TZZ, | |
540 | & ATEMP, ZTEMP, IVFLAG ) | |
541 | IF ( LRESMP ) RETURN | |
542 | NEVT = NEVT + 1 | |
543 | * +-------------------------------------------------------------------* | |
544 | * | | |
545 | IF ( NNHAD + NHAD .GT. MXPNUC ) THEN | |
546 | WRITE (LUNOUT,1101) NNHAD | |
547 | STOP | |
548 | END IF | |
549 | * | | |
550 | * +-------------------------------------------------------------------* | |
551 | ||
552 | * +-------------------------------------------------------------------* | |
553 | * | Looping over the produced particles | |
554 | DO 8 JJ=1,NHAD | |
555 | NNHAD = NNHAD+1 | |
556 | J = NNHAD | |
557 | PXNU(J) = PXH(JJ) | |
558 | PYNU(J) = PYH(JJ) | |
559 | PZNU(J) = PZH(JJ) | |
560 | HEPNU(J) = HEPH(JJ) | |
561 | AMNU(J) = AMH(JJ) | |
562 | IBARNU(J)= IBARH(JJ) | |
563 | ICHNU(J) = ICHH(JJ) | |
564 | NRENU(J) = NREH(JJ) | |
565 | ANNU(J) = ANH(JJ) | |
566 | * | Updating energy and momentum accumulators | |
567 | PUX = PUX+PXNU(J) | |
568 | PUY = PUY+PYNU(J) | |
569 | PUZ = PUZ+PZNU(J) | |
570 | EUZ = EUZ+HEPNU(J) | |
571 | ICU = ICU+ICHNU(J) | |
572 | IBU = IBU+IBARNU(J) | |
573 | 8 CONTINUE | |
574 | * | | |
575 | * +-------------------------------------------------------------------* | |
576 | 888 CONTINUE | |
577 | * IF (PPROJ .LT. 4.D0) INIT=1 | |
578 | **** print and test energy conservation | |
579 | * | |
580 | ETOT = EEPROJ + KTARP * ( AMNUCL (1) - EBNDNG (1) ) | |
581 | & + KTARN * ( AMNUCL (2) - EBNDNG (2) ) + EFRM | |
582 | ICTOT = ICH (KPROJ) + KTARP | |
583 | IBTOT = IBAR (KPROJ) + KTARP + KTARN | |
584 | EMIN = 1.D-10 * ETOT | |
585 | PMIN = 1.D-10 * PPPROJ | |
586 | * +-------------------------------------------------------------------* | |
587 | * | | |
588 | IF (ICTOT .NE. ICU .OR. IBTOT .NE. IBU) THEN | |
589 | LRESMP = .TRUE. | |
590 | WRITE(LUNERR,*)' NUCEVT-FEREVT CHARGE CONSERVATION FAILURE:', | |
591 | & ' ICU,ICTOT,IBU,IBTOT',ICU,ICTOT,IBU,IBTOT | |
592 | END IF | |
593 | * | | |
594 | * +-------------------------------------------------------------------* | |
595 | PXCHCK = PPPROJ * TXI + PSEA * TXX | |
596 | PYCHCK = PPPROJ * TYI + PSEA * TYY | |
597 | PZCHCK = PPPROJ * TZI + PSEA * TZZ | |
598 | * +-------------------------------------------------------------------* | |
599 | * | | |
600 | IF ( ABS ( PUX - PXFRM - PXCHCK ) .GT. PMIN ) THEN | |
601 | LRESMP = .TRUE. | |
602 | WRITE(LUNERR,*)' NUCEVT-FEREVT PX CONSERVATION FAILURE:', | |
603 | & ' PUX,PXFRM+PXCHCK',PUX,PXFRM+PXCHCK | |
604 | END IF | |
605 | * | | |
606 | * +-------------------------------------------------------------------* | |
607 | * +-------------------------------------------------------------------* | |
608 | * | | |
609 | IF ( ABS ( PUY - PYFRM - PYCHCK ) .GT. PMIN ) THEN | |
610 | WRITE(LUNERR,*)' NUCEVT-FEREVT PY CONSERVATION FAILURE:', | |
611 | & ' PUY,PYFRM+PYCHCK*TYY',PUY,PYFRM+PYCHCK | |
612 | LRESMP = .TRUE. | |
613 | END IF | |
614 | * | | |
615 | * +-------------------------------------------------------------------* | |
616 | * +-------------------------------------------------------------------* | |
617 | * | | |
618 | IF ( ABS ( PUZ - PZFRM - PZCHCK ) .GT. PMIN ) THEN | |
619 | LRESMP = .TRUE. | |
620 | WRITE(LUNERR,*)' NUCEVT-FEREVT PZ CONSERVATION FAILURE:', | |
621 | & ' PUZ,PZFRM+PZCHCK',PUZ,PZFRM+PZCHCK | |
622 | END IF | |
623 | * | | |
624 | * +-------------------------------------------------------------------* | |
625 | * +-------------------------------------------------------------------* | |
626 | * | | |
627 | IF ( ABS ( ETOT - EUZ ) .GT. EMIN ) THEN | |
628 | LRESMP = .TRUE. | |
629 | WRITE(LUNERR,*)' NUCEVT-FEREVT E CONSERVATION FAILURE:', | |
630 | & ' EUZ,ETOT',EUZ,ETOT | |
631 | END IF | |
632 | * | | |
633 | * +-------------------------------------------------------------------* | |
634 | EUZ0 = EUZ-NEVT*AM(1) | |
635 | IF (INIT.NE.1) GO TO 11 | |
636 | WRITE(LUNOUT,12)NNHAD,KPROJ,PPPROJ,EPROJ,PUX,PUY,PUZ,EUZ0, | |
637 | *ICU,IBU,NNHAD | |
638 | 12 FORMAT (1X,2I5,6F12.6,3I5) | |
639 | * +-------------------------------------------------------------------* | |
640 | * | | |
641 | DO 13 I=1,NNHAD | |
642 | WRITE(LUNOUT,14)I,NRENU(I),ICHNU(I),IBARNU(I),ANNU(I),PXNU(I), | |
643 | * PYNU(I),PZNU(I),HEPNU(I),AMNU(I) | |
644 | 14 FORMAT (1X,4I5,A8,6F12.6) | |
645 | 13 CONTINUE | |
646 | * | | |
647 | * +-------------------------------------------------------------------* | |
648 | INIT=0 | |
649 | 11 CONTINUE | |
650 | RETURN | |
651 | END |