5 * Revision 1.1.1.1 1995/10/24 10:19:57 cernlib
9 #include "geant321/pilot.h"
10 *CMZ : 3.21/02 29/03/94 15.41.43 by S.Giani
15 *=== nucevv ===========================================================*
17 SUBROUTINE NUCEVV ( NNHAD, KPROJ, PPPROJ, EKPROJ, TXI, TYI, TZI )
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
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
29 *----------------------------------------------------------------------*
32 C**************************************************************
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
41 C IBARNU = BARYONIC NUMBER
42 C ANNU = NAME OF THE PARTICLE
43 C NRENU = TYPE NUMBER OF THE PARTICLE
44 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)
63 DATA PTHRSH / 16*5.D+00,2*2.5D+00,5.D+00,3*2.5D+00,8*5.D+00,
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 *----------------------------------------------------------------------*
94 EPROJ = EKPROJ + AMPROJ
99 * Set Atemp and Ztemp to their initial values
100 ATEMP = IBTAR - IGREYP - IGREYN
101 ZTEMP = ICHTAR - IGREYP
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)
110 * We have now sampled Nsea, the number of quark-antiquark pairs
111 * interacting besides the valence interaction (total interactions
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
120 EM = EKPROJ * ( XSEA(I) + XASEA(I) )
121 * | +----------------------------------------------------------------*
122 * | | Selection of the ddbar or qqbar composition
124 IF ( RNDM (1) .LT. 0.5D+00 ) THEN
127 * | +----------------------------------------------------------------*
133 * | +----------------------------------------------------------------*
136 * | +----------------------------------------------------------------*
137 * | | Energy and/or momentum is too low!!!
138 IF ( EKM .LT. ETHSEA .OR. PPROJ .LT. PTHRSH (KPTOIP(KPROJ)) )
140 * | | +-------------------------------------------------------------*
142 IF ( I .LT. NSEA ) THEN
144 XSEA(II) = XSEA(II) + XSEA(I)
145 XASEA(II)= XASEA(II) + XASEA(I)
147 * | | +-------------------------------------------------------------*
152 * | | +-------------------------------------------------------------*
154 * | | +-------------------------------------------------------------*
155 * | | | Kt is the index of the target nucleon (1=proton,8=neutron)
156 IF ( KT .EQ. 1 ) THEN
159 ZNCOLL = ZNCOLL - 1.D+00
161 * | | +-------------------------------------------------------------*
167 * | | +-------------------------------------------------------------*
169 ANCOLL = ANCOLL - 1.D+00
173 * | +----------------------------------------------------------------*
174 AMPIO2 = AMPIO * AMPIO
175 EPROLD = EKPROJ + AMPROJ
177 * | After selection of the meson energy, Ekproj is updated
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
186 * | | +-------------------------------------------------------------*
190 * | | | +----------------------------------------------------------*
191 * | | | | Kt is the index of the target nucleon (1=proton,8=neutron)
192 IF ( KT .EQ. 1 ) THEN
196 * | | | +----------------------------------------------------------*
202 * | | | +----------------------------------------------------------*
206 * | | +-------------------------------------------------------------*
210 * | +----------------------------------------------------------------*
211 * | Now the wounded nucleus is selected inside Corevt
213 ATEMP = ATEMP - 1.D+00
214 IF ( KT .EQ. 1 ) ZTEMP = ZTEMP - 1.D+00
215 * | +---------------------------------------------------------------*
217 IF ( NINT ( ATEMP ) .EQ. 1 ) THEN
220 KTINC = IJTARG ( NSEA + 1 )
221 PM = SQRT ( EM**2 - AMPIO2 )
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
231 EKPROJ = EKPROJ - 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 * | | | +----------------------------------------------------------*
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
249 ATEMP = ATEMP + 1.D+00
250 IF ( KT .EQ. 1 ) ZTEMP = ZTEMP + 1.D+00
251 * | | | | | +----------------------------------------------------*
255 * | | | | | | +-------------------------------------------------*
256 * | | | | | | | Kt is the index of the target nucleon
257 * | | | | | | | (1=proton,8=neutron)
258 IF ( KT .EQ. 1 ) THEN
262 * | | | | | | +-------------------------------------------------*
268 * | | | | | | +-------------------------------------------------*
272 * | | | | | +----------------------------------------------------*
273 PPROJ = SQRT ( EKPROJ * ( EKPROJ + 2.D+00
277 * | | | | +-------------------------------------------------------*
279 ELSE IF ( DELTAE .GE. EKPROJ ) THEN
280 WRITE ( LUNERR,* )' Nucevv: sea call impossible ',
281 & ' to get Umin2, go to resampling'
286 * | | | | +-------------------------------------------------------*
289 * | | +-<|--<--<--<--< We need more invariant mass!!
292 * | | | +----------------------------------------------------------*
294 * | | +-------------------------------------------------------------*
295 * | | Now we will divide Eleft and Pleft between the two
298 ELACMS = 0.5D+00 * ( UMO2 + AMLAST**2 - AMINC**2 )
300 EINCMS = UMO - ELACMS
301 PCMS = SQRT ( ELACMS**2 - AMLAST**2 )
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
326 EINCT = EINCP + EINCN
327 IF ( EINCT .GT. 0.D+00 ) THEN
328 EINCP = EINCP - EINCP / EINCT * TVGREY
329 EINCN = EINCN - EINCN / EINCT * TVGREY
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
338 CALL FEREVV ( NHAD, IM, KT, PM, EKM, TXX, TYY, TZZ, ATEMP,
340 LDIFFR (KPTOIP(IM)) = LDSAVE
343 * | +---------------------------------------------------------------*
349 PPROJ = SQRT ( EKPROJ * ( EKPROJ + 2.D+00 * AMPROJ ) )
350 IF ( PPROJ .GT. PPROLD ) THEN
351 PPROJ = PPROLD * EPROJ / EPROLD
353 CALL FERSEA ( NHAD, IM, KT, EKM, TXX, TYY, TZZ, ATEMP,
354 & ZTEMP, KPROJ, EPROJ, PPROJ, EPROLD, PPROLD )
358 * | +----------------------------------------------------------------*
360 * | +----------------------------------------------------------------*
362 IF ( NNHAD + NHAD .GT. MXPNUC ) THEN
363 WRITE (LUNOUT,1101) NNHAD
364 1101 FORMAT(' NHAD IN NUCEVT TOO BIG CHANGE DIMENSION',
369 * | +----------------------------------------------------------------*
371 * | +----------------------------------------------------------------*
372 * | | Looping over the produced particles!!!
385 * | | Updating energy and momentum accumulators
394 * | +----------------------------------------------------------------*
397 * +-------------------------------------------------------------------*
400 * Computing the actual pproj
402 EPROJ = EKPROJ + AMPROJ
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 * | +---------------------------------------------------------------*
420 IF ( PPROJ - 1.5D+00 * SIGN ( ONEONE, DKPR30 )
422 WRITE ( LUNERR,* )' Nucevt: Pproj < 4 GeV/c!!!', PPROJ,
426 * | +---------------------------------------------------------------*
429 * +------------------------------------------------------------------*
430 * Now the wounded nucleus is selected inside Corevt
431 KT = IJTARG ( NSEA + 1 )
432 * Update the Nsea value
434 *or IF (INIT.EQ.1) WRITE(LUNOUT,162)NEVT,NHAD,IM,KPROJ,KT,PM,
435 *or &EM,EKM,PM,EKPROJ,PPROJ,EPROJ
437 * Now ferevt performs the interaction: no longer a meson but the actual
440 ATEMP = ATEMP - 1.D+00
441 IF ( KT .EQ. 1 ) ZTEMP = ZTEMP - 1.D+00
442 * +-------------------------------------------------------------------*
444 IF ( NINT ( ATEMP ) .EQ. 1 ) THEN
447 IF ( ZNOW .GT. 0.D+00 ) THEN
454 UMIN2 = ( AMINC + AMLAST )**2
455 ITJ = MIN ( KTINC, 2 )
456 * | +----------------------------------------------------------------*
457 * | | Selection of the invariant mass of the two last nucleon
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 * | | +-------------------------------------------------------------*
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 * | | | +----------------------------------------------------------*
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'
480 * | | | +----------------------------------------------------------*
481 EKPROJ = EKPROJ - DELTAE
483 PPROJ = SQRT ( EKPROJ * ( EKPROJ + 2.D+00 * AMPROJ ) )
486 PXFRM = PXFRM + TXX * PLA
487 PYFRM = PYFRM + TYY * PLA
488 PZFRM = PZFRM + TZZ * PLA
491 * | +-<|--<--<--<--< We need more invariant mass!!
494 * | | +-------------------------------------------------------------*
496 * | +----------------------------------------------------------------*
497 * | Now we will divide Eleft and Pleft between the two
500 ELACMS = 0.5D+00 * ( UMO2 + AMLAST**2 - AMINC**2 ) / UMO
501 EINCMS = UMO - ELACMS
502 PCMS = SQRT ( ELACMS**2 - AMLAST**2 )
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
527 EINCT = EINCP + EINCN
528 IF ( EINCT .GT. 0.D+00 ) THEN
529 EINCP = EINCP - EINCP / EINCT * TVGREY
530 EINCN = EINCN - EINCN / EINCT * TVGREY
535 * +-------------------------------------------------------------------*
536 * The last argument for ferevv is set to 1 to signal that it
537 * is a true valence call
539 CALL FEREVV ( NHAD, KPROJ, KT, PPROJ, EKPROJ, TXX, TYY, TZZ,
540 & ATEMP, ZTEMP, IVFLAG )
543 * +-------------------------------------------------------------------*
545 IF ( NNHAD + NHAD .GT. MXPNUC ) THEN
546 WRITE (LUNOUT,1101) NNHAD
550 * +-------------------------------------------------------------------*
552 * +-------------------------------------------------------------------*
553 * | Looping over the produced particles
566 * | Updating energy and momentum accumulators
575 * +-------------------------------------------------------------------*
577 * IF (PPROJ .LT. 4.D0) INIT=1
578 **** print and test energy conservation
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
585 PMIN = 1.D-10 * PPPROJ
586 * +-------------------------------------------------------------------*
588 IF (ICTOT .NE. ICU .OR. IBTOT .NE. IBU) THEN
590 WRITE(LUNERR,*)' NUCEVT-FEREVT CHARGE CONSERVATION FAILURE:',
591 & ' ICU,ICTOT,IBU,IBTOT',ICU,ICTOT,IBU,IBTOT
594 * +-------------------------------------------------------------------*
595 PXCHCK = PPPROJ * TXI + PSEA * TXX
596 PYCHCK = PPPROJ * TYI + PSEA * TYY
597 PZCHCK = PPPROJ * TZI + PSEA * TZZ
598 * +-------------------------------------------------------------------*
600 IF ( ABS ( PUX - PXFRM - PXCHCK ) .GT. PMIN ) THEN
602 WRITE(LUNERR,*)' NUCEVT-FEREVT PX CONSERVATION FAILURE:',
603 & ' PUX,PXFRM+PXCHCK',PUX,PXFRM+PXCHCK
606 * +-------------------------------------------------------------------*
607 * +-------------------------------------------------------------------*
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
615 * +-------------------------------------------------------------------*
616 * +-------------------------------------------------------------------*
618 IF ( ABS ( PUZ - PZFRM - PZCHCK ) .GT. PMIN ) THEN
620 WRITE(LUNERR,*)' NUCEVT-FEREVT PZ CONSERVATION FAILURE:',
621 & ' PUZ,PZFRM+PZCHCK',PUZ,PZFRM+PZCHCK
624 * +-------------------------------------------------------------------*
625 * +-------------------------------------------------------------------*
627 IF ( ABS ( ETOT - EUZ ) .GT. EMIN ) THEN
629 WRITE(LUNERR,*)' NUCEVT-FEREVT E CONSERVATION FAILURE:',
630 & ' EUZ,ETOT',EUZ,ETOT
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,
638 12 FORMAT (1X,2I5,6F12.6,3I5)
639 * +-------------------------------------------------------------------*
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)
647 * +-------------------------------------------------------------------*