]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/fluka/nucevv.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Default compile option changed to -g (Alpha)
[u/mrichter/AliRoot.git] / GEANT321 / fluka / nucevv.F
CommitLineData
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*
31C
32C**************************************************************
33C /NUCPAR/
34C PARTICLES CREATED IN NUCEVT
35C PXNU,PYNU,PZNU = X-, Y- AND Z-COMPONENTS OF THE
36C LAB MOMENTUM OF THE SECONDARY. COORDINATE SYSTEM
37C DEFINED BY THE PRIMARY PARTICLE.
38C HEPNU = TOTAL ENERGY IN THE LAB
39C AMNU = REST MASS
40C ICHNU = CHARGE
41C IBARNU = BARYONIC NUMBER
42C ANNU = NAME OF THE PARTICLE
43C NRENU = TYPE NUMBER OF THE PARTICLE
44C**************************************************************
45C
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
96C
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
2302020 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
3641101 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
4031000 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
4592040 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