]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/fluka/eventv.F
git://git.uio.no / u / mrichter / AliRoot.git / blame
? search:
summary | shortlog | log | commit | commitdiff | tree
blob | blame (incremental) | history | HEAD
Open flukaaf.dat in the ALICE root
[u/mrichter/AliRoot.git] / GEANT321 / fluka / eventv.F
CommitLineData
fe4da5cc 1*
2* $Id$
3*
4* $Log$
5* Revision 1.1.1.1 1995/10/24 10:19:55 cernlib
6* Geant
7*
8*
9#include "geant321/pilot.h"
10*CMZ : 3.21/02 29/03/94 15.41.42 by S.Giani
11*-- Author :
12*$ CREATE EVENTV.FOR
13*COPY EVENTV
14*
15*=== eventv ===========================================================*
16*
17 SUBROUTINE EVENTV ( IJIJ, POO, EKE, TXI, TYI, TZI, WE, IMAT )
18
19#include "geant321/dblprc.inc"
20#include "geant321/dimpar.inc"
21#include "geant321/iounit.inc"
22*
23*----------------------------------------------------------------------*
24* *
25* Created on 7 september 1989 by Alfredo Ferrari *
26* INFN - Milan *
27* *
28* Last change on 15 april 1993 by Alfredo Ferrari *
29* *
30* This routine is a strongly modified and improved version of the *
31* original Eventq routine of the Fluka package: it requires modified *
32* versions of many other routines *
33* *
34* Eventv90: new version by A. Ferrari, INFN-Milan and CERN-SPS, 7-9-89*
35* main modifications: *
36* - new treatment for hydrogen (following Moehring's *
37* suggestion) *
38* - new commons and coding for energy, momentum, *
39* electric and baryonic charge conservation *
40* - new correlations both in the "low" (p < 5 GeV/c) *
41* and in the "high" energy models *
42* - complete revision of kinematics for both models *
43* - tentative treatment for Lambdas, Sigmas with the *
44* "high" energy model down to 2.5 GeV/c to overcome *
45* the impossibility of Hadrin to treat these parti- *
46* cles *
47* - tentative introduction of Xsi0, Xsi-, Omega and *
48* their antiparticles and of Sigmabars's *
49* - complete new treatment of cascade nucleons *
50* - accurate treatment of binding energy effects, *
51* using atomic mass tables and the proper isotopic *
52* composition of elements *
53* - possibility to use an evaporation module (derived *
54* from the Hetc-Kfa one) after the interaction *
55* - possibility to use a nuclear gamma deexcitation *
56* module (written by P.Sala, INFN-Milan) to produce *
57* deexcitation gammas after the evaporation step *
58* *
59*----------------------------------------------------------------------*
60*
61 PARAMETER ( COS45 = 0.7071067811865475D+00 )
62 PARAMETER ( SIN30 = 0.5D+00 )
63 PARAMETER ( COS30 = 0.8660254037844387D+00 )
64 PARAMETER ( CSNPRN = 50.D+00 * CSNNRM )
65#include "geant321/balanc.inc"
66#include "geant321/corinc.inc"
67#include "geant321/depnuc.inc"
68#include "geant321/fheavy.inc"
69#include "geant321/finlsp3.inc"
70#include "geant321/finuc.inc"
71#include "geant321/hadflg.inc"
72#include "geant321/hadpar.inc"
73#include "geant321/isotop.inc"
74#include "geant321/mapa.inc"
75#include "geant321/nucdat.inc"
76#include "geant321/nucpar.inc"
77#include "geant321/part.inc"
78#include "geant321/parevt.inc"
79#include "geant321/resnuc.inc"
80 COMMON / FKCASF / PKFRMI, COSTH, PKIN
81 COMMON /FKEVNT/ LNUCRI, LHADRI
82 LOGICAL LNUCRI, LHADRI, LINCCK, LACCEP, LCHTYP
83 REAL RNDM(1)
84* Note Pthrsh is also in Ferevv!!
85 DIMENSION SOPPP(2),EXSOP(2),PTHRSH(39),IJNUCR(39)
86* Use this "save" if common dbgdbg is commented
87 SAVE PTHRSH, IPRI, IEVT, IJNUCR
88 DATA PTHRSH / 16*5.D+00,2*2.5D+00,5.D+00,3*2.5D+00,8*5.D+00,
89 & 9*2.5D+00 /
90 DATA IJNUCR / 16*1,2*0,1,3*0,8*1,9*0 /
91 DATA IPRI / 0 /
92 DATA IEVT / 0 /
93*
94* +-------------------------------------------------------------------*
95* | Heavy ions are treated in eventa: now switched off!!
96 IF (IJIJ .EQ. 30) THEN
97 STOP 'HEAVY ION'
98 END IF
99* |
100* +-------------------------------------------------------------------*
101* First of all: get the "part" index of the incoming "paprop" ijij
102 IJ = IPTOKP (IJIJ)
103 IEVT = IEVT + 1
104 AMCHCK = 0.5D+00 * ( POO * POO - EKE * EKE ) / EKE
105 AHELP = ABS ( AMCHCK - AM ( IJ ) )
106* +-------------------------------------------------------------------*
107* |
108 IF ( AHELP .GT. ANGLGB * AM ( IJ ) ) THEN
109 POO = SQRT ( EKE * ( EKE + 2.D+00 * AM (IJ) ) )
110 END IF
111* |
112* +-------------------------------------------------------------------*
113*
114* Variable initialization for conservation checks
115*
116 PTTOT = POO
117 PXTTOT = PTTOT*TXI
118 PYTTOT = PTTOT*TYI
119 PZTTOT = PTTOT*TZI
120 ICHTAR = NINT(ZTAR(IMAT))
121 ICHTOT = ICHTAR + ICH(IJ)
122* +-------------------------------------------------------------------*
123* | Choice of the mass number of the target nucleus: use the input
124* | one if any
125 IF ( MSSNUM (IMAT) .GT. 0 ) THEN
126 IBTAR = MSSNUM (IMAT)
127* |
128* +-------------------------------------------------------------------*
129* | Choice of the mass number of the target nucleus according
130* | to the natural isotopic composition of element ichtar
131 ELSE
132 CALL GRNDM(RNDM,1)
133 RNDISO = RNDM (1)
134* | +----------------------------------------------------------------*
135* | | Loop on the stable isotopes
136 DO 25 IS = ISONDX (1,ICHTAR), ISONDX (2,ICHTAR) - 1
137 RNDISO = RNDISO - ABUISO (IS)
138 IF ( RNDISO .LT. 0.D+00 ) GO TO 30
139 25 CONTINUE
140* | |
141* | +----------------------------------------------------------------*
142 IS = ISONDX (2,ICHTAR)
143 30 CONTINUE
144 IBTAR = ISOMNM (IS)
145 END IF
146* |
147* +-------------------------------------------------------------------*
148 IBTOT = IBTAR + IBAR(IJ)
149 BBTAR = IBTAR
150 ZZTAR = ICHTAR
151 ZTO103 = ZZTAR**0.3333333333333333D+00
152* +-------------------------------------------------------------------*
153* |
154 IF ( IBTAR .EQ. 1 ) THEN
155 ITTA = 8 - 7 * ICHTAR
156 ETTOT = AM (ITTA) + EKE + AM(IJ)
157 AMNTAR = AM (ITTA)
158* |
159* +-------------------------------------------------------------------*
160* | The following should be done with the proper mass of the nuclide
161 ELSE
162* AMNTAR = BBTAR * AMUC12
163 AMMTAR = BBTAR * AMUAMU + 0.001D+00 * FKENER ( BBTAR, ZZTAR )
164 AMNTAR = AMMTAR - ZZTAR * AMELEC + ELBNDE ( ICHTAR )
165 ETTOT = AMNTAR + EKE + AM (IJ)
166 END IF
167* |
168* +-------------------------------------------------------------------*
169* +-------------------------------------------------------------------*
170* | Kaon-short and kaon-long are transformed into a kaon0 or a
171* | kaon0-bar
172 IF (IJ .EQ. 12 .OR. IJ. EQ. 19) THEN
173 IJ = 24
174 CALL GRNDM(RNDM,1)
175 IF (RNDM(1) .GT. 0.5D0) IJ = 25
176* |
177* +-------------------------------------------------------------------*
178* | Pi0 quark configurations are selected according to 50% uubar and
179* | 50% ddbar
180 ELSE IF ( IJ .EQ. 23 .OR. IJ .EQ. 26 ) THEN
181* | +----------------------------------------------------------------*
182* | |
183 CALL GRNDM(RNDM,1)
184 IF ( RNDM (1) .LT. 0.5D+00 ) THEN
185 IJ = 23
186* | |
187* | +----------------------------------------------------------------*
188* | |
189 ELSE
190 IJ = 26
191 END IF
192* | |
193* | +----------------------------------------------------------------*
194 END IF
195* |
196* +-------------------------------------------------------------------*
197 ETEPS = MAX ( 1.D-10 * EKE, 1.D-10 )
198 50 CONTINUE
199*-->-->-->-->--> Come here for resampling
200* Reset all the accumulators
201 NP = NP0
202 NPHEAV = 0
203 LRNFSS = .FALSE.
204 LEVDIF = .FALSE.
205 LHADRI = .FALSE.
206 ICHTOT = ICHTAR + ICH(IJ)
207 TV = 0.D+00
208 TVRECL = 0.D+00
209 TVHEAV = 0.D+00
210 TVCMS = 0.D+00
211 EOTEST = ETTOT
212 LRESMP = .FALSE.
213 LLASTN = .FALSE.
214 LLAST1 = .FALSE.
215 EUZ = 0.D+00
216 PUX = 0.D+00
217 PUY = 0.D+00
218 PUZ = 0.D+00
219 TVEUZ = 0.D+00
220 ICU = 0
221 IBU = 0
222 EINTR = 0.D+00
223 PXINTR = 0.D+00
224 PYINTR = 0.D+00
225 PZINTR = 0.D+00
226 ICINTR = 0
227 IBINTR = 0
228 ENUCR = 0.D+00
229 PXNUCR = 0.D+00
230 PYNUCR = 0.D+00
231 PZNUCR = 0.D+00
232 ICNUCR = 0
233 IBNUCR = 0
234 EFRM = 0.D+00
235 PXFRM = 0.D+00
236 PYFRM = 0.D+00
237 PZFRM = 0.D+00
238 PSEA = 0.D+00
239 TVGRE0 = 0.D+00
240 TVGREY = 0.D+00
241 IGREYP = 0
242 IGREYN = 0
243 KTARP = 0
244 KTARN = 0
245 IEVAPL = 0
246 IEVAPH = 0
247 IEVPRO = 0
248 IEVNEU = 0
249 IEVDEU = 0
250 IEVTRI = 0
251 IEV3HE = 0
252 IEV4HE = 0
253 IDEEXG = 0
254 ANOW = BBTAR
255 ZNOW = ZZTAR
256 DSOPP = 0.D+00
257 NNHAD = 0
258 RN1GSC = -1.D+00
259 RN2GSC = -1.D+00
260 IF ( POO .GE. PTHRSH (IJIJ) ) GO TO 200
261* Below 5 GeV/c use Nucrin:
262 100 CONTINUE
263 ANCOLL = 1.D+00
264* +-------------------------------------------------------------------*
265* | Check for Hydrogen
266 IF ( IBTAR .NE. 1 ) THEN
267* | +----------------------------------------------------------------*
268* | | Steering for Peanut: very temporary
269* | | Protons and Neutrons below 260 MeV:
270 IF ( EKE .LT. 0.260D+00 .AND. ( IJ .EQ. 1 .OR. IJ .EQ. 8 )
271 & .AND. IBTAR .GT. 4 ) THEN
272 CALL PEANUT ( IJ, EKE, POO, TXI, TYI, TZI, WE )
273 IF ( LRESMP ) GO TO 50
274 RETURN
275* | |
276* | +----------------------------------------------------------------*
277* | | Stopping pi-s:
278 ELSE IF ( IJ .EQ. 14 .AND. EKE .LT. 2.D+00 * GAMMIN .AND. IBTAR
279 & .GT. 4 ) THEN
280 CALL PEANUT ( IJ, EKE, POO, TXI, TYI, TZI, WE )
281 IF ( LRESMP ) GO TO 50
282 RETURN
283 END IF
284* | |
285* | +----------------------------------------------------------------*
286* | Now at first we sample the number of interactions
287 CALL CORRIN ( ZZTAR, BBTAR, IJ, POO, EKE )
288 END IF
289 LNUCRI = .TRUE.
290 NP = NP0
291 ELAB = EKE + AM(IJ)
292* Redefinition of particle types: the recovery of proper charge and
293* strangeness conservation after the interaction is not yet implemented
294 KPROJ = IJ
295 LCHTYP = .FALSE.
296* +-------------------------------------------------------------------*
297* |
298 IF ( KPROJ .LE. 30 ) THEN
299 GO TO ( 2017, 2018, 2999, 2020, 2021, 2022, 2999, 2999, 2999,
300 & 2026 ), KPROJ - 16
301 GO TO 2999
302* |
303* +-------------------------------------------------------------------*
304* |
305 ELSE
306 GO TO ( 2031, 2032, 2033, 2034, 2035, 2036, 2037, 2038, 2039 ),
307 & KPTOIP (KPROJ) - 30
308 GO TO 2999
309 END IF
310* |
311* +-------------------------------------------------------------------*
312* +-------------------------------------------------------------------*
313* | Lambdas are considered as neutrons: after the interaction
314* | a possible way to recover strangeness and charge conservation
315* | is to flip a pi- --> K-, or an eventual neutron into
316* | a Sigma0/Lambda,
317* | (we must flip a d quark into an s quark), or a pi0 --> K0bar,
318* | or a p --> Sigma+
319 2017 CONTINUE
320 LCHTYP = .TRUE.
321 KPROJ = 8
322* | Adjusting mass and momentum for lambdas (now n)
323 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
324 PTTOT = POO
325 PXTTOT = PTTOT*TXI
326 PYTTOT = PTTOT*TYI
327 PZTTOT = PTTOT*TZI
328 GO TO 2999
329* |
330* +-------------------------------------------------------------------*
331* | Alambdas are considered as antineutrons : after the interaction
332* | a possible way to recover strangeness and charge conservation
333* | is to flip a pi+ --> K+, or an eventual neutronbar into a
334* | Lambdabar, (we must flip a dbar quark into an sbar quark),
335* | or a pi0 --> K0, or a pbar --> Sigma+bar (ASigma-) ( there
336* | are also possible ways through eta, rho and omega mesons )
337 2018 CONTINUE
338 LCHTYP = .TRUE.
339 KPROJ = 9
340* | Adjusting mass and momentum for alambdas (now an)
341 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
342 PTTOT = POO
343 PXTTOT = PTTOT*TXI
344 PYTTOT = PTTOT*TYI
345 PZTTOT = PTTOT*TZI
346 GO TO 2999
347* |
348* +-------------------------------------------------------------------*
349* | Sigma-s are considered as neutrons : after the interaction
350* | a possible way to recover strangeness and charge conservation
351* | is to flip a pi+ --> K0bar, or an eventual neutron into a Sigma-,
352* | (we must flip a u quark into an s quark), or a pi0 --> K-,
353* | or a p --> Sigma0/Lambda
354 2020 CONTINUE
355 LCHTYP = .TRUE.
356 KPROJ = 8
357* | The following card must be commented following the strangeness
358* | and charge conservation recovery implementation
359 ICHTOT = ICHTOT + 1
360* | Adjusting mass and momentum for Sigma-s (now n)
361 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
362 PTTOT = POO
363 PXTTOT = PTTOT*TXI
364 PYTTOT = PTTOT*TYI
365 PZTTOT = PTTOT*TZI
366 GO TO 2999
367* |
368* +-------------------------------------------------------------------*
369* | Sigma+s are considered as protons : after the interaction
370* | a possible way to recover strangeness and charge conservation
371* | is to flip a pi- --> K-, or an eventual proton into a Sigma+,
372* | (we must flip a d quark into an s quark), or a pi0 --> K0bar,
373* | or a n --> Sigma0/Lambda
374 2021 CONTINUE
375 LCHTYP = .TRUE.
376 KPROJ = 1
377* | Adjusting mass and momentum for Sigma+s (now p)
378 PTTOT = POO
379 PXTTOT = PTTOT*TXI
380 PYTTOT = PTTOT*TYI
381 PZTTOT = PTTOT*TZI
382 GO TO 2999
383* |
384* +-------------------------------------------------------------------*
385* | Sigma0s are considered as neutrons: after the interaction
386* | a possible way to recover strangeness and charge conservation
387* | is to flip a pi- --> K-, or an eventual neutron into
388* | a Sigma0/Lambda,
389* | (we must flip a d quark into an s quark), or a pi0 --> K0bar,
390* | or a p --> Sigma+
391 2022 CONTINUE
392 LCHTYP = .TRUE.
393 KPROJ = 8
394* | Adjusting mass and momentum for Sigma0s (now n)
395 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
396 PTTOT = POO
397 PXTTOT = PTTOT*TXI
398 PYTTOT = PTTOT*TYI
399 PZTTOT = PTTOT*TZI
400 GO TO 2999
401* |
402* +-------------------------------------------------------------------*
403* | Pi0 must be 23 for Nucrin
404 2026 CONTINUE
405 KPROJ = 23
406 GO TO 2999
407* |
408* +-------------------------------------------------------------------*
409* | ASigma-s are considered as pbars: after the interaction
410* | a possible way to recover strangeness and charge conservation
411* | is to flip a pi+ --> K+, or an eventual pbar into a ASigma-,
412* | (we must flip a dbar quark into an sbar quark), or a pi0 --> K0,
413* | or a nbar --> ASigma0/ALambda
414 2031 CONTINUE
415 LCHTYP = .TRUE.
416 KPROJ = 2
417* | Adjusting mass and momentum for ASigma-s (now pbar)
418 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
419 PTTOT = POO
420 PXTTOT = PTTOT*TXI
421 PYTTOT = PTTOT*TYI
422 PZTTOT = PTTOT*TZI
423 GO TO 2999
424* |
425* +-------------------------------------------------------------------*
426* | ASigma0s are considered as nbar: after the interaction
427* | a possible way to recover strangeness and charge conservation
428* | is to flip a pi+ --> K+, or an eventual nbar into
429* | a ASigma0/ALambda,
430* | (we must flip a dbar quark into an sbar quark), or a pi0 --> K0,
431* | or a pbar --> ASigma-
432 2032 CONTINUE
433 LCHTYP = .TRUE.
434 KPROJ = 9
435* | Adjusting mass and momentum for ASigma0s (now nbar)
436 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
437 PTTOT = POO
438 PXTTOT = PTTOT*TXI
439 PYTTOT = PTTOT*TYI
440 PZTTOT = PTTOT*TZI
441 GO TO 2999
442* |
443* +-------------------------------------------------------------------*
444* | ASigma+s are considered as nbar : after the interaction
445* | a possible way to recover strangeness and charge conservation
446* | is to flip a pi- --> K0, or an eventual nbar into a ASigma+,
447* | (we must flip a ubar quark into an sbar quark), or a pi0 --> K+,
448* | or a pbar --> ASigma0/ALambda
449 2033 CONTINUE
450 LCHTYP = .TRUE.
451 KPROJ = 9
452* | The following card must be commented following the strangeness
453* | and charge conservation recovery implementation
454 ICHTOT = ICHTOT - 1
455* | Adjusting mass and momentum for ASigma+s (now nbar)
456 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
457 PTTOT = POO
458 PXTTOT = PTTOT*TXI
459 PYTTOT = PTTOT*TYI
460 PZTTOT = PTTOT*TZI
461 GO TO 2999
462* |
463* +-------------------------------------------------------------------*
464* | Xsi0s are considered as neutrons: after the interaction
465* | a possible way to recover strangeness and charge conservation
466* | is to flip an eventual neutron into a Xsi0,
467* | or a proton --> Sigma+ and a pi- --> K- (or a pi0 --> K0bar) etc.
468* | (we must flip two d quarks into s quarks)
469 2034 CONTINUE
470 LCHTYP = .TRUE.
471 KPROJ = 8
472* | Adjusting mass and momentum for Xsi0s (now n)
473 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
474 PTTOT = POO
475 PXTTOT = PTTOT*TXI
476 PYTTOT = PTTOT*TYI
477 PZTTOT = PTTOT*TZI
478 GO TO 2999
479* |
480* +-------------------------------------------------------------------*
481* | AXsi0s are considered as nbars: after the interaction
482* | a possible way to recover strangeness and charge conservation
483* | is to flip an eventual nbar into a AXsi0,
484* | or a pbar --> ASigma- and a pi+ --> K+ (or a pi0 --> K0) etc.
485* | (we must flip two dbar quarks into sbar quarks)
486 2035 CONTINUE
487 LCHTYP = .TRUE.
488 KPROJ = 9
489* | Adjusting mass and momentum for AXsi0s (now nbar)
490 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
491 PTTOT = POO
492 PXTTOT = PTTOT*TXI
493 PYTTOT = PTTOT*TYI
494 PZTTOT = PTTOT*TZI
495 GO TO 2999
496* |
497* +-------------------------------------------------------------------*
498* | Xsi-s are considered as protons: after the interaction
499* | a possible way to recover strangeness and charge conservation
500* | is to flip an eventual proton into a Xsi-,
501* | or a neutron --> Sigma- and a pi+ --> K0bar (or a pi0 --> K-) etc.
502* | (we must flip two u quarks into s quarks)
503 2036 CONTINUE
504 LCHTYP = .TRUE.
505 KPROJ = 1
506* | The following card must be commented following the strangeness
507* | and charge conservation recovery implementation
508 ICHTOT = ICHTOT + 2
509* | Adjusting mass and momentum for Xsi-s (now p)
510 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
511 PTTOT = POO
512 PXTTOT = PTTOT*TXI
513 PYTTOT = PTTOT*TYI
514 PZTTOT = PTTOT*TZI
515 GO TO 2999
516* |
517* +-------------------------------------------------------------------*
518* | AXsi+s are considered as pbars: after the interaction
519* | a possible way to recover strangeness and charge conservation
520* | is to flip an eventual pbar into a AXsi+,
521* | or a nbar --> ASigma+ and a pi- --> K0 (or a pi0 --> K+) etc.
522* | (we must flip two ubar quarks into sbar quarks)
523 2037 CONTINUE
524 LCHTYP = .TRUE.
525 KPROJ = 2
526* | The following card must be commented following the strangeness
527* | and charge conservation recovery implementation
528 ICHTOT = ICHTOT - 2
529* | Adjusting mass and momentum for AXsi+s (now pbar)
530 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
531 PTTOT = POO
532 PXTTOT = PTTOT*TXI
533 PYTTOT = PTTOT*TYI
534 PZTTOT = PTTOT*TZI
535 GO TO 2999
536* |
537* +-------------------------------------------------------------------*
538* | Omega-s are considered as protons: after the interaction there
539* | are many possible ways to recover strangeness and charge
540* | conservation (we must flip a d quark into an s quark and two u
541* | quarks into s quarks).
542 2038 CONTINUE
543 LCHTYP = .TRUE.
544 KPROJ = 1
545* | The following card must be commented following the strangeness
546* | and charge conservation recovery implementation
547 ICHTOT = ICHTOT + 2
548* | Adjusting mass and momentum for Omega-s (now p)
549 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
550 PTTOT = POO
551 PXTTOT = PTTOT*TXI
552 PYTTOT = PTTOT*TYI
553 PZTTOT = PTTOT*TZI
554 GO TO 2999
555* |
556* +-------------------------------------------------------------------*
557* | Omegabar+s are considered as pbars: after the interaction there
558* | are many possible ways to recover strangeness and charge
559* | conservation (we must flip a dbar quark into an sbar quark and
560* | two ubar quarks into sbar quarks).
561 2039 CONTINUE
562 LCHTYP = .TRUE.
563 KPROJ = 2
564* | The following card must be commented following the strangeness
565* | and charge conservation recovery implementation
566 ICHTOT = ICHTOT - 2
567* | Adjusting mass and momentum for AOmega+s (now pbar)
568 POO = SQRT ( ( ELAB + AM(KPROJ) ) * ( ELAB - AM(KPROJ) ) )
569 PTTOT = POO
570 PXTTOT = PTTOT*TXI
571 PYTTOT = PTTOT*TYI
572 PZTTOT = PTTOT*TZI
573 GO TO 2999
574* |
575* +-------------------------------------------------------------------*
576 2999 CONTINUE
577* +-------------------------------------------------------------------*
578* | Check for hydrogen!!!!
579 IF ( IBTAR .EQ. 1 ) THEN
580 TV = 0.D+00
581 NP = NP0
582 ITTA = 8 - 7 * ICHTAR
583 IF ( KPROJ .EQ. 14 .AND. EKE .LT. 2.D+00 * GAMMIN ) THEN
584 CALL PMPRAB ( KPROJ, EKE, POO, TXI, TYI, TZI, WE )
585 RETURN
586 END IF
587* | Set a flag to avoid elastic collision in Hadrin
588 IELFLG = +1
589* | Set a flag to fully exploit charge exchange
590 ICXFLG = 0
591 CALL HADRIV ( KPROJ, POO, ELAB, TXI, TYI, TZI, ITTA )
592 IELFLG = -1
593 ICXFLG = -1
594* | +----------------------------------------------------------------*
595* | |
596 IF ( IH .LE. 1 ) THEN
597* | | +-------------------------------------------------------------*
598* | | |
599 IF ( ITH (1) .EQ. KPROJ ) THEN
600 IH = IH + 1
601 ITH (2) = 1
602 CXH (2) = TXI
603 CYH (2) = TYI
604 CZH (2) = TZI
605 PLH (2) = 0.D+00
606 ELH (2) = AM (1)
607* | | |
608* | | +-------------------------------------------------------------*
609* | | |
610 ELSE IF ( ITH (1) .EQ. 1 ) THEN
611 IH = IH + 1
612 ITH (2) = KPROJ
613 CXH (2) = TXI
614 CYH (2) = TYI
615 CZH (2) = TZI
616 PLH (2) = 0.D+00
617 ELH (2) = AM (KPROJ)
618* | | |
619* | | +-------------------------------------------------------------*
620* | | |
621 ELSE
622 END IF
623* | | |
624* | | +-------------------------------------------------------------*
625 END IF
626* | |
627* | +----------------------------------------------------------------*
628
629* | +----------------------------------------------------------------*
630* | | Looping over the particles produced in hadrin
631 DO 110 J=1,IH
632 NP = NP + 1
633 TKI(NP) = ELH(J) - AM(ITH(J))
634* | | +-------------------------------------------------------------*
635* | | |
636 IF ( TKI(NP) .GE. 0.D+00 ) THEN
637 KPART(NP) = ITH(J)
638 CXR(NP) = CXH(J)
639 CYR(NP) = CYH(J)
640 CZR(NP) = CZH(J)
641 PLR(NP) = PLH(J)
642 WEI(NP) = WE
643* | | | Updating conservation counters
644 ENUCR = ENUCR + ELH(J)
645 PXNUCR = PXNUCR + PLR(NP)*CXR(NP)
646 PYNUCR = PYNUCR + PLR(NP)*CYR(NP)
647 PZNUCR = PZNUCR + PLR(NP)*CZR(NP)
648 ICNUCR = ICNUCR + ICH(KPART(NP))
649 IBNUCR = IBNUCR + IBAR(KPART(NP))
650* | | |
651* | | +-------------------------------------------------------------*
652* | | |
653 ELSE IF ( ABS ( TKI (NP) ) .LE. ANGLGB ) THEN
654* | | | Updating conservation counters
655 ENUCR = ENUCR + ELH(J)
656 ICNUCR = ICNUCR + ICH(ITH(J))
657 IBNUCR = IBNUCR + IBAR(ITH(J))
658 NP = NP - 1
659* | | |
660* | | +-------------------------------------------------------------*
661* | | |
662 ELSE
663 NP = NP - 1
664 END IF
665* | | |
666* | | +-------------------------------------------------------------*
667 110 CONTINUE
668* | |
669* | +----------------------------------------------------------------*
670 KTARP = 1
671 ANOW = 0.D+00
672 ZNOW = 0.D+00
673 ICRES = 0
674 IBRES = 0
675 ERES = 0.D+00
676 EOTEST = EOTEST - ENUCR
677* | +----------------------------------------------------------------*
678* | |
679 IF ( ABS (EOTEST) .GT. ETEPS ) THEN
680 WRITE (LUNERR,*)' Eventq: eotest failure with Hadrin',
681 & EOTEST
682 LRESMP = .TRUE.
683 GO TO 50
684* |
685* +-<|--<--<--<--<--< go to resampling
686 END IF
687* | |
688* | +----------------------------------------------------------------*
689 RETURN
690 END IF
691* | end hydrogen check
692* +-------------------------------------------------------------------*
693*
694* Now calling Nucrin to produce secondary particles, which are put
695* into the /Finuc/ common
696*
697 CALL NUCRIV ( KPROJ, ELAB, TXI, TYI, TZI, BBTAR,
698 & ZZTAR, RHO(IMAT) )
699 IF ( LRESMP ) GO TO 50
700*--<--<--<--<--< go to resampling if something was wrong
701 TVGRE0 = TV - TVGREY
702 TVEUZ = 0.D+00
703 TVTENT = TV
704 TV = 0.D+00
705 IF (NP .GT. MXP) GO TO 400
706* +-------------------------------------------------------------------*
707* | Looping over the particles produced in nucrin
708* | No need for Nucrin produced particles to
709* | change the numbering
710 DO 101 I=NP0+1,NP
711 I1 = KPART(I)
712 WEI(I) = WE
713 TKI(I) = TKI(I) - AM(I1)
714 TXYZ = CXR (I)**2 + CYR (I)**2 + CZR (I)**2
715 IF ( ABS (TXYZ-1.D0) .GT. CSNPRN ) THEN
716 WRITE ( LUNERR,* )
717 & ' **** Bad normalized cosines from Nucriv!!!',I,TXYZ,CXR(I),
718 & CYR(I),CZR(I)
719 TXYZ = TXYZ - 1.D+00
720 TXYZ = 1.D+00 - 0.5D+00 * TXYZ + 0.375D+00 * TXYZ * TXYZ
721 CXR (I) = CXR (I) * TXYZ
722 CYR (I) = CYR (I) * TXYZ
723 CZR (I) = CZR (I) * TXYZ
724 TXYZ = CXR (I)**2 + CYR (I)**2 + CZR (I)**2
725 IF ( ABS (TXYZ-1.D0) .GT. CSNPRN ) THEN
726 LRESMP = .TRUE.
727 GO TO 50
728 END IF
729 END IF
730 IF (TKI(I) .GE. 0.D+00) GO TO 101
731 TKI(I) = 0.D+00
732 PLR(I) = 0.D+00
733 101 CONTINUE
734* |
735* +-------------------------------------------------------------------*
736 GO TO 500
737C
738C********************************************************************
739C GENERATE FIRST THE LOW ENERGY NUCLEONS FROM THE INTRANUCLEAR
740C CASCADE - CHECK IF ACCEPTABLE.
741C NOTE THAT AINEL AND RAKEKA ARE USING THE OLD PARTICLE NUMBERING
742C********************************************************************
743C
744 200 CONTINUE
745 NP = NP0
746 LNUCRI = .FALSE.
747* +-------------------------------------------------------------------*
748* | Check for hydrogen!!!!
749 IF ( IBTAR .EQ. 1 ) THEN
750 ANCOLL = 1.D+00
751 NHAD = 0
752 KPROJ = IJ
753 KTARG = 1
754 EPROJ = EKE + AM (KPROJ)
755 UMO = SQRT ( AM (1)**2 + AM (KPROJ)**2 + 2.D+00 * AM (1) *
756 & EPROJ )
757* | +----------------------------------------------------------------*
758* | | Now diffractive events are taken into account for projectile
759* | | -proton collisions!!! A. Ferrari and J. Ranft, 25-3-90
760 IF ( LDIFFR (KPTOIP(KPROJ)) ) THEN
761 CALL GRNDM(RNDM,1)
762 IF ( RNDM (1) .LT. FRDIFF ) THEN
763 IF ( POO .GT. PTHDFF ) THEN
764 LEVDIF = .TRUE.
765 CALL DIFEVV ( NHAD, KPROJ, KTARG, POO, EPROJ, UMO )
766 ELSE
767 CALL HADEVV ( NHAD, KPROJ, KTARG, POO, EPROJ, UMO )
768 END IF
769 ELSE
770 CALL HADEVV ( NHAD, KPROJ, KTARG, POO, EPROJ, UMO )
771 END IF
772* | |
773* | +----------------------------------------------------------------*
774* | |
775 ELSE
776 CALL HADEVV ( NHAD, KPROJ, KTARG, POO, EPROJ, UMO )
777 END IF
778* | |
779* | +----------------------------------------------------------------*
780* | +----------------------------------------------------------------*
781* | | Loop over particles produced in hadevt/difevt
782 DO 207 J=1,NHAD
783 NP = NP + 1
784 CFE = 1.D+00
785 SFE = 0.D+00
786* | | Get the "paprop" index
787 KPART(NP) = KPTOIP (NREH(J))
788 IF ( KPART (NP) .EQ. 0 ) THEN
789 WRITE (LUNERR,*)
790 & ' **** Charmed particle produced in Hadevv/Difevv',
791 & NREH(J),HEPH(J),AMH(J),' ****'
792 KPART (NP) = NREH (J)
793 END IF
794 PLR(NP) = SQRT ( PXH(J)**2 + PYH(J)**2 + PZH(J)**2 )
795 PTH = SQRT( PXH(J)**2 + PYH(J)**2 )
796 CDE = PZH(J) / PLR(NP)
797 SDE = PTH / PLR(NP)
798 IF (SDE .GE. ANGLGB) THEN
799 CFE = PXH(J) / PTH
800 SFE = PYH(J) / PTH
801 END IF
802 CALL TTRANS ( TXI, TYI, TZI, CDE, SDE, SFE, CFE,
803 & CXR(NP), CYR(NP), CZR(NP) )
804 TKI(NP) = HEPH(J) - AMH(J)
805 WEI(NP) = WE
806* | | +-------------------------------------------------------------*
807* | | | Updating conservation counters
808 IF ( TKI(NP) .GT. 0.D+00 ) THEN
809 EUZ = EUZ + HEPH(J)
810 PUX = PUX + PLR(NP)*CXR(NP)
811 PUY = PUY + PLR(NP)*CYR(NP)
812 PUZ = PUZ + PLR(NP)*CZR(NP)
813 ICU = ICU + ICH(NREH(J))
814 IBU = IBU + IBAR(NREH(J))
815* | | |
816* | | +-------------------------------------------------------------*
817* | | |
818 ELSE
819 NP = NP-1
820 END IF
821* | | |
822* | | +-------------------------------------------------------------*
823 207 CONTINUE
824* | | end loop
825* | +----------------------------------------------------------------*
826 ANOW = 0.D+00
827 ZNOW = 0.D+00
828 KTARP = 1
829 ICRES = 0
830 IBRES = 0
831 ERES = 0.D+00
832 EOTEST = EOTEST - EUZ
833* | +----------------------------------------------------------------*
834* | |
835 IF ( ABS (EOTEST) .GT. ETEPS ) THEN
836 WRITE (LUNERR,*)' Eventq: eotest failure with',
837 & ' Hadevt/Diffevt',EOTEST
838 END IF
839* | |
840* | +----------------------------------------------------------------*
841 RETURN
842 END IF
843* | end hydrogen check
844* +-------------------------------------------------------------------*
845* Now at first we sample the number of interactions
846 CALL COREVT ( ZZTAR, BBTAR, IJ, POO, EKE )
847*
848* First decide how much energy will go into the intranuclear cascade
849*
850 EKIN = EKE
851* Here starts the cascade particle production module!!!
852 SOPPP (1) = EINCP
853 SOPPP (2) = EINCN
854 ARECL = MAX ( ANOW - 1.D+00, 1.D+00 )
855 KREJE = 0
856 NGREYT = NGREYP + NGREYN
857 TVTENT = ANCOLL * ( AV0WEL - AEFRMA ) + TVGRE0
858 TMP015 = 0.15D+00
859 EEXTRA = 2.D+00* MIN ( TWOTWO * ( EKUPNU (1) + EKUPNU (2) )
860 & , TMP015 *
861 & ( EKE - EINCP - EINCN - TVTENT ), EINCP + EINCN )
862 EEXTRA = MAX ( EEXTRA, EKUPNU (1), EKUPNU (2) )
863 LINCCK = .FALSE.
864 PXCASC = 0.D+00
865 PYCASC = 0.D+00
866 PZCASC = 0.D+00
867 PTXINT = 0.D+00
868 PTYINT = 0.D+00
869 PPINTR = 0.D+00
870 EKRECL = 0.D+00
871 IEXTRN = 0
872 IEXTRP = 0
873 LACCEP = .FALSE.
874 IGREYT = 0
875 IF ( NGREYT .EQ. 0 ) GO TO 206
876* +-------------------------------------------------------------------*
877* | Now looping until we reach the correct number of cascade
878* | nucleons
879 205 CONTINUE
880 IGREYT = IGREYP + IGREYN
881 DSOPP = SOPPP (1) + SOPPP (2)
882* | Now it is not possible to reduce too heavily the available energy
883 IF ( EKIN - TVGRE0 .LT. 0.5D+00 * EKE .OR. PZCASC .GE.
884 & 0.5D+00 * POO ) GO TO 206
885* | +----------------------------------------------------------------*
886* | |
887 IF ( IGREYT .GE. NGREYT ) THEN
888 AEXTRA = 2.D+00 * DSOPP / ( EKINAV (1) + EKINAV (2)
889 & + ESWELL (1) + ESWELL (2) )
890 & / ( IGREYT - NGREYT + 1 )
891 CALL GRNDM(RNDM,1)
892 RNDUMO = RNDM (1)
893* | | +-------------------------------------------------------------*
894* | | |
895 IF ( RNDUMO .LT. AEXTRA .AND. LACCEP .AND. NINT (ANOW)
896 & .GT. 0 ) THEN
897 REJE = ZNOW / ANOW
898* | | | +----------------------------------------------------------*
899* | | | |
900 CALL GRNDM(RNDM,1)
901 IF ( RNDM (1) .LT. REJE ) THEN
902 IEXTRP = MIN ( 1, NINT (ZNOW), NGREYP / 2 )
903 IF ( NGREYP + IEXTRP - IGREYP .LE. 0 ) GO TO 206
904* | | | |
905* | | | +----------------------------------------------------------*
906* | | | |
907 ELSE
908* IEXTRN = MIN ( IEXTRN + 1, NINT (ANOW - ZNOW), 2 )
909 IEXTRN = MIN ( 1, NINT (ANOW - ZNOW), NGREYN / 2 )
910 IF ( NGREYN + IEXTRN - IGREYN .LE. 0 ) GO TO 206
911 END IF
912* | | | |
913* | | | +----------------------------------------------------------*
914* | | |
915* | | +-------------------------------------------------------------*
916* | | |
917 ELSE
918 GO TO 206
919* | | |
920* | | |-->-->-->-->--> We have finished particles!!
921 END IF
922* | | |
923* | | +-------------------------------------------------------------*
924* | |
925* | +----------------------------------------------------------------*
926* | |
927 ELSE IF ( DSOPP .LE. - EEXTRA ) THEN
928 GO TO 206
929* | |
930* | |-->-->-->-->--> We have finished energy
931 END IF
932* | |
933* | +----------------------------------------------------------------*
934 LACCEP = .FALSE.
935* | +----------------------------------------------------------------*
936* | | Check how many nucleons are still available for cascade
937 IF ( ANOW .LT. 1.5D+00 ) THEN
938* | | +-------------------------------------------------------------*
939* | | | Check if at least two high energy collisions are foreseen
940* | | | if yes the proper energy / momentum balance for the last
941* | | | two nucleons will be done inside Ferevv
942 IF ( ANCOLL .GT. 1.5D+00 ) THEN
943* | | |
944* | | +-------------------------------------------------------------*
945* | | | Only the valence collision is foreseen: we are forced
946* | | | to perform here the balance for the 2 last nucleons!
947* | | | (One might be used here for cascading, the other
948* | | | will be used by Hadevv / Difevv)
949 ELSE
950 LLASTN = .TRUE.
951 KTLAST = IJTARG (1)
952 AMLAST = AM ( KTLAST )
953* | | | +----------------------------------------------------------*
954* | | | |
955 IF ( NINT ( ZNOW ) .GT. 0 ) THEN
956 KTINC = 1
957 IGREYP = IGREYP + 1
958* | | | |
959* | | | +----------------------------------------------------------*
960* | | | |
961 ELSE
962 KTINC = 8
963 IGREYN = IGREYN + 1
964 END IF
965* | | | |
966* | | | +----------------------------------------------------------*
967 AMINC = AM (KTINC)
968 UMIN2 = ( AMLAST + AMINC )**2
969 DSOPP = MAX ( DSOPP, AV0WEL - AEFRMA )
970 216 CONTINUE
971 EKIN = EKIN - DSOPP
972* | | | +----------------------------------------------------------*
973* | | | | Check if we can get enough invariant mass
974 IF ( EKIN .LE. 0.2D+00 * EKE ) THEN
975 WRITE ( LUNOUT, * )
976 & ' *** Eventv: impossible to get enough invariant',
977 & ' mass for the last two nucleons! ****'
978 WRITE ( LUNERR, * )
979 & ' *** Eventv: impossible to get enough invariant',
980 & ' mass for the last two nucleons! ****'
981 LRESMP = .TRUE.
982 GO TO 50
983* | | | |---<---<---< Go to resampling
984 END IF
985* | | | |
986* | | | +----------------------------------------------------------*
987 ELEFT = ETTOT - EINTR - EKIN - AM (IJ)
988 PLA = SQRT ( EKIN * ( EKIN + 2.D+00 * AM (IJ) ) )
989 PXLEFT = PXTTOT - PXINTR - PLA * TXI
990 PYLEFT = PYTTOT - PYINTR - PLA * TYI
991 PZLEFT = PZTTOT - PZINTR - PLA * TZI
992* | | | Now we will divide Eleft and Pleft between the two
993* | | | left nucleons!
994 UMO2 = ELEFT**2 - PXLEFT**2 - PYLEFT**2 - PZLEFT**2
995* | | | +----------------------------------------------------------*
996* | | | |
997 IF ( UMO2 .LE. UMIN2 ) THEN
998* | | | | +-------------------------------------------------------*
999* | | | | |
1000 IF ( KTINC .EQ. 1 ) THEN
1001 EINCP = EINCP + DSOPP
1002* | | | | |
1003* | | | | +-------------------------------------------------------*
1004* | | | | |
1005 ELSE
1006 EINCN = EINCN + DSOPP
1007 END IF
1008* | | | | |
1009* | | | | +-------------------------------------------------------*
1010 GO TO 216
1011 END IF
1012* | | | |
1013* | | | +----------------------------------------------------------*
1014 UMO = SQRT ( UMO2 )
1015 ELACMS = 0.5D+00 * ( UMO2 + AMLAST**2 - AMINC**2 )
1016 & / UMO
1017 EINCMS = UMO - ELACMS
1018 PCMS = SQRT ( ELACMS**2 - AMLAST**2 )
1019 GAMCM = ELEFT / UMO
1020 ETAX = PXLEFT / UMO
1021 ETAY = PYLEFT / UMO
1022 ETAZ = PZLEFT / UMO
1023 CALL RACO ( CXXINC, CYYINC, CZZINC )
1024 PCMSX = PCMS * CXXINC
1025 PCMSY = PCMS * CYYINC
1026 PCMSZ = PCMS * CZZINC
1027* | | | Now go back from the CMS frame to the lab frame!!!
1028* | | | First the inc nucleon:
1029 ETAPCM = PCMSX * ETAX + PCMSY * ETAY + PCMSZ * ETAZ
1030 NP = NP + 1
1031 IF (NP .GT. MXP) GO TO 400
1032 KPART(NP)= KTINC
1033 WEI (NP) = WE
1034 TKI (NP) = GAMCM * EINCMS + ETAPCM - AMINC
1035 PHELP = ETAPCM / (GAMCM + 1.D0) + EINCMS
1036 CXXINC = PCMSX + ETAX * PHELP
1037 CYYINC = PCMSY + ETAY * PHELP
1038 CZZINC = PCMSZ + ETAZ * PHELP
1039* | | | Updating conservation counters
1040 EINTR = EINTR + TKI(NP) + AMINC
1041 PXINTR = PXINTR + CXXINC
1042 PYINTR = PYINTR + CYYINC
1043 PZINTR = PZINTR + CZZINC
1044 ICINTR = ICINTR + ICH (KPART(NP))
1045 IBINTR = IBINTR + IBAR (KPART(NP))
1046 PLR (NP) = SQRT (CXXINC * CXXINC + CYYINC * CYYINC +
1047 & CZZINC * CZZINC)
1048 CXR (NP) = CXXINC / PLR (NP)
1049 CYR (NP) = CYYINC / PLR (NP)
1050 CZR (NP) = CZZINC / PLR (NP)
1051* | | | Now the high energy collision nucleon
1052 EKLAST = GAMCM * ELACMS - ETAPCM - AMLAST
1053 PHELP = - ETAPCM / (GAMCM + 1.D0) + ELACMS
1054 PXLAST = - PCMSX + ETAX * PHELP
1055 PYLAST = - PCMSY + ETAY * PHELP
1056 PZLAST = - PCMSZ + ETAZ * PHELP
1057* | | | Now we perform a Lorentz transformation to the rest fra-
1058* | | | me of the "last" nucleon, with the projectile momentum
1059* | | | along the +z axis
1060 AMPROJ = AM (IJ)
1061 EPROJ = EKIN + AMPROJ
1062 ECHCK = EPROJ + AMLAST + EKLAST
1063 PXCHCK = PLA * TXI + PXLAST
1064 PYCHCK = PLA * TYI + PYLAST
1065 PZCHCK = PLA * TZI + PZLAST
1066 UMO = SQRT ( ECHCK**2 - PXCHCK**2 - PYCHCK**2 -
1067 & PZCHCK**2 )
1068 EPROJX = 0.5D+00 * ( UMO**2 - AMPROJ**2 - AMLAST**2 )
1069 & / AMLAST
1070 PPROJX = SQRT ( EPROJX**2 - AMPROJ**2 )
1071 ETOTX = EPROJX + AMLAST
1072* | | | Now set the parameters for the Lorentz transformation
1073 AAFACT = ECHCK + ETOTX
1074 BBFACT = PPROJX - PZCHCK
1075 DDENOM = ETOTX * AAFACT - PPROJX * BBFACT
1076 GAMTRA = ( ECHCK * AAFACT + PPROJX * BBFACT ) / DDENOM
1077 ETAZ = - BBFACT * AAFACT / DDENOM
1078 ETAX = PXCHCK * ( GAMTRA + 1.D+00 ) / AAFACT
1079 ETAY = PYCHCK * ( GAMTRA + 1.D+00 ) / AAFACT
1080 PLABS = PPROJX
1081 ELABS = EPROJX
1082* | | | +----------------------------------------------------------*
1083* | | | |
1084 IF ( .NOT. LDIFFR (KPTOIP(IJ)) .OR. PLABS .LE. PTHDFF )
1085 & THEN
1086 RUUN = 1.D0
1087* | | | |
1088* | | | +----------------------------------------------------------*
1089* | | | |
1090 ELSE
1091 CALL GRNDM(RNDM,1)
1092 RUUN = RNDM (1)
1093 END IF
1094* | | | |
1095* | | | +----------------------------------------------------------*
1096 NHAD = 0
1097* | | | +----------------------------------------------------------*
1098* | | | | Diffractive event
1099 IF ( RUUN .LE. FRDIFF ) THEN
1100 LEVDIF = .TRUE.
1101 CALL DIFEVV ( NHAD, IJ, KTLAST, PLABS, ELABS, UMO )
1102* | | | |
1103* | | | +----------------------------------------------------------*
1104* | | | | Usual event
1105 ELSE
1106 CALL HADEVV ( NHAD, IJ, KTLAST, PLABS, ELABS, UMO )
1107 END IF
1108* | | | |
1109* | | | +----------------------------------------------------------*
1110* | | | Now we have the particles in the nucleon rest frame,
1111* | | | transform back in the lab frame
1112* | | | +----------------------------------------------------------*
1113* | | | | Looping over the produced particles
1114 DO 236 I = 1, NHAD
1115 NP = NP + 1
1116 IF (NP .GT. MXP) GO TO 400
1117 CALL ALTRA ( GAMTRA, ETAX, ETAY, ETAZ, PXH (I),
1118 & PYH (I), PZH (I), HEPH (I), PLR (NP),
1119 & CXR (NP), CYR (NP), CZR (NP), ELR )
1120* | | | | Updating conservation counters
1121 EUZ = EUZ + ELR
1122 PUX = PUX + CXR (NP)
1123 PUY = PUY + CYR (NP)
1124 PUZ = PUZ + CZR (NP)
1125 KPART (NP) = KPTOIP ( NREH (I) )
1126 IF ( KPART (NP) .EQ. 0 ) THEN
1127 WRITE (LUNERR,*)
1128 & ' **** Charmed particle produced in Hadevv',
1129 & NREH(I),ELR,AMH(I),' ****'
1130 KPART (NP) = NREH (I)
1131 END IF
1132 ICU = ICU + ICH (NREH(I))
1133 IBU = IBU + IBAR (NREH(I))
1134 CXR (NP) = CXR (NP) / PLR (NP)
1135 CYR (NP) = CYR (NP) / PLR (NP)
1136 CZR (NP) = CZR (NP) / PLR (NP)
1137 TKI (NP) = ELR - AMH (I)
1138 WEI (NP) = WE
1139 236 CONTINUE
1140* | | | |
1141* | | | +----------------------------------------------------------*
1142* | | | Now perform the standard tests for energy and momentum
1143* | | | conservation
1144 ERES = ETTOT - EUZ - ENUCR - EINTR
1145 PXRES = PXTTOT - PUX - PXNUCR - PXINTR
1146 PYRES = PYTTOT - PUY - PYNUCR - PYINTR
1147 PZRES = PZTTOT - PUZ - PZNUCR - PZINTR
1148 ICRES = ICHTOT - ICU - ICNUCR - ICINTR
1149 IBRES = IBTOT - IBU - IBNUCR - IBINTR
1150 PTRES = MAX ( ABS (PXRES), ABS (PYRES), ABS (PZRES) )
1151* | | | +----------------------------------------------------------*
1152* | | | |
1153 IF ( IBRES .NE. 0 .OR. ICRES .NE. 0 .OR. ABS (ERES)
1154 & .GT. 1.D-10*EPROJ .OR. PTRES .GT. 1.D-10*PTTOT )
1155 & THEN
1156 WRITE ( LUNERR, * )
1157 & ' Eventq: last nucleon failure!!', ICRES,
1158 & IBRES, REAL (ERES), REAL (PXRES),
1159 & REAL (PYRES), REAL (PZRES)
1160 LRESMP = .TRUE.
1161 GO TO 50
1162* | | |--|--<--<--<--<--< go to resampling
1163 END IF
1164* | | | |
1165* | | | +----------------------------------------------------------*
1166 PTRES = 0.D+00
1167 RETURN
1168 END IF
1169* | | |
1170* | | +-------------------------------------------------------------*
1171 END IF
1172* | |
1173* | +----------------------------------------------------------------*
1174
1175* | +----------------------------------------------------------------*
1176* | |
1177 IF ( NINT ( ANOW - ZNOW ) .LE. 0 ) THEN
1178 REJE = 1.D+00
1179* | |
1180* | +----------------------------------------------------------------*
1181* | |
1182 ELSE
1183 REJE = MAX ( ZNOW * ( NGREYP + IEXTRP - IGREYP ),
1184 & ZERZER )
1185 HELP = REJE + ( ANOW - ZNOW ) * ( NGREYN +
1186 & IEXTRN - IGREYN )
1187 IF ( HELP .GT. 0.D+00 ) THEN
1188 REJE = REJE / HELP
1189 ELSE
1190 REJE = 1.D+00
1191 WRITE(LUNERR,*)' EVENTV: HELP=0, ANOW,ZNOW',ANOW,ZNOW
1192 & ,' NGREYN,IGREYN,IEXTRN',
1193 & NGREYN,IGREYN,IEXTRN,
1194 & ' NGREYP,IGREYP,IEXTRP',
1195 & NGREYP,IGREYP,IEXTRP
1196 END IF
1197 END IF
1198* | |
1199* | +----------------------------------------------------------------*
1200
1201* | +----------------------------------------------------------------*
1202* | |
1203 CALL GRNDM(RNDM,1)
1204 IF ( RNDM (1) .LT. REJE ) THEN
1205 ILO = 1
1206 ILLO = 2
1207 KP = 1
1208* | |
1209* | +----------------------------------------------------------------*
1210* | |
1211 ELSE
1212 ILO = 2
1213 ILLO = 1
1214 KP = 8
1215 END IF
1216* | |
1217* | +----------------------------------------------------------------*
1218 CALL RAKEKV ( ILO, EXSOP, BBTAR, TKIN, TSTRCK, PLA, ARECL,
1219 & TKRECL, EFERMI, CDE, SDE )
1220 IF ( EKIN - TKIN .LT. 0.5D+00 * EKE ) GO TO 206
1221* |--|--<--<--<--<--< avoid to deplete too much the energy
1222* | +----------------------------------------------------------------*
1223* | | Now check if the nucleon energy is acceptable:
1224 IF ( SOPPP (ILO) .LT. TKIN ) THEN
1225* | | +-------------------------------------------------------------*
1226* | | |
1227 IF ( TKIN - SOPPP (ILO) .LT. SOPPP (ILLO) + 1.5D+00 *
1228 & EEXTRA ) THEN
1229 SOPPP (ILLO) = SOPPP (ILLO) + SOPPP (ILO) - TKIN
1230 & - EXSOP (ILO)
1231 SOPPP (ILO) = 0.D+00
1232* | | |
1233* | | +-------------------------------------------------------------*
1234* | | |
1235 ELSE
1236 KREJE = KREJE + 1
1237 IF ( KREJE .GT. 10 ) GO TO 206
1238 SOPPP (ILLO) = SOPPP (ILLO) + SOPPP (ILO)
1239 SOPPP (ILO) = 0.D+00
1240 GO TO 205
1241 END IF
1242* | | |
1243* | | +-------------------------------------------------------------*
1244* | |
1245* | +----------------------------------------------------------------*
1246* | |
1247 ELSE
1248 SOPPP (ILO) = SOPPP (ILO) - TKIN - EXSOP (ILO)
1249 END IF
1250* | |
1251* | +----------------------------------------------------------------*
1252 IF ( NINT ( ANOW ) .LE. 0 ) GO TO 206
1253* | +----------------------------------------------------------------*
1254* | | Update the current atomic and mass number
1255 IF ( ILO .EQ. 1 ) THEN
1256 IF ( NINT ( ZNOW ) .LE. 0 ) GO TO 206
1257 IGREYP = IGREYP + 1
1258 ZNOW = ZNOW - 1.D+00
1259* | |
1260* | +----------------------------------------------------------------*
1261* | |
1262 ELSE
1263 IGREYN = IGREYN + 1
1264 END IF
1265* | |
1266* | +----------------------------------------------------------------*
1267 ANOW = ANOW - 1.D+00
1268* | Now make the kinematical tests!!
1269 FRSURV = ARECL / BBTAR
1270 NP = NP + 1
1271 IF ( NP .GT. MXP ) GO TO 400
1272 EKIN = EKIN - TKIN
1273* Note the change!!
1274 LINCCK = ARECL .LT. 30.D+00 .AND. FRSURV .LE. 0.33D+00
1275* | +----------------------------------------------------------------*
1276* | |
1277 IF ( .NOT. LINCCK ) THEN
1278 CALL SFECFE ( SFE, CFE )
1279* | |
1280* | +----------------------------------------------------------------*
1281* | |
1282 ELSE
1283 IF ( CDE .GT. 0.D+00 ) THEN
1284 CDE = CDE * FRSURV / 0.33D+00
1285 SDE = SQRT ( 1.D+00 - CDE**2 )
1286 END IF
1287 CALL SFECFE ( SFE, CFE )
1288 END IF
1289* | |
1290* | +----------------------------------------------------------------*
1291 PTXINT = PTXINT + PLA * CFE * SDE
1292 PTYINT = PTYINT + PLA * SFE * SDE
1293 PPINTR = PPINTR + PLA * CDE
1294 CALL TTRANS ( TXI, TYI, TZI, CDE, SDE, SFE, CFE,
1295 & CXR (NP), CYR (NP), CZR (NP) )
1296* | Generated nucleon is acceptable - save it
1297 LACCEP = .TRUE.
1298 TVGRE0 = TVGRE0 + EXSOP (ILO)
1299 TVGREY = TVGREY + TKIN + TKRECL - TSTRCK - EBNDNG (ILO)
1300 EKRECL = EKRECL + TKRECL
1301 ARECL = MAX ( ARECL - 1.D+00, 1.D+00 )
1302 KPART(NP)= KP
1303 WEI(NP) = WE
1304 TKI(NP) = TSTRCK
1305 PLR(NP) = PLA
1306* | Updating the counters for spent momenta
1307* | of the cascade nucleons (the momentum
1308* | spent by the high energy particles)
1309 SINTH = SQRT ( 1.D+00 - COSTH * COSTH )
1310 PLA = SQRT ( ( EFERMI + TKIN )**2 - AMNUSQ (ILO) )
1311 PZLEFT = PLA * CDE - COSTH * PKFRMI
1312 IF ( PZLEFT + PZCASC .LT. EKE - EKIN + TVGRE0 ) THEN
1313 PZCASC = EKE - EKIN + TVGRE0
1314 PXCASC = PXCASC + CFE * ( PLA * SDE - SINTH * PKFRMI )
1315 PYCASC = PYCASC + SFE * ( PLA * SDE - SINTH * PKFRMI )
1316 ELSE
1317 PZCASC = PZCASC + PZLEFT
1318 PXCASC = PXCASC + CFE * ( PLA * SDE - SINTH * PKFRMI )
1319 PYCASC = PYCASC + SFE * ( PLA * SDE - SINTH * PKFRMI )
1320 END IF
1321* | Updating conservation counters
1322 EINTR = EINTR + TKI(NP) + AM (KPART(NP))
1323 PXINTR = PXINTR + PLR(NP) * CXR(NP)
1324 PYINTR = PYINTR + PLR(NP) * CYR(NP)
1325 PZINTR = PZINTR + PLR(NP) * CZR(NP)
1326 ICINTR = ICINTR + ICH (KPART(NP))
1327 IBINTR = IBINTR + IBAR (KPART(NP))
1328 GO TO 205
1329* |
1330* +--<--<--<--<--< go to sample another nucleon
1331 206 CONTINUE
1332* +-------------------------------------------------------------------*
1333* | First check charge and baryonic number conservation:
1334 IF ( IGREYP .NE. ICINTR .OR. ( IGREYP + IGREYN ) .NE.
1335 & IBINTR ) THEN
1336 WRITE ( LUNERR,* )' Eventq: charge or baryon number',
1337 & ' conservation failure in the Inc',
1338 & ' section!!', IGREYP, ICINTR,
1339 & IGREYP + IGREYN, IBINTR
1340 LRESMP = .TRUE.
1341 GO TO 50
1342*--|--<--<--<--<--< go to resampling
1343 END IF
1344* |
1345* +-------------------------------------------------------------------*
1346 TVTENT = TVGRE0
1347* Tvgrey is already inside Eincp and Eincn since Rakekv updates the
1348* estimated excitation energy (and Tvgre0, Eincp, Eincn and the
1349* Soppp array)
1350 EINCP = EINCP - SOPPP (1)
1351 EINCN = EINCN - SOPPP (2)
1352 EKIN = EKIN - TVTENT
1353* +-------------------------------------------------------------------*
1354* | Check if we have not spent too much energy!!!
1355 IF ( EKIN .LT. 0.5D+00 * EKE ) THEN
1356 LRESMP = .TRUE.
1357 WRITE ( LUNERR,* ) ' Eventq: Ekin after inc too low!! ',
1358 & EKIN, EKE, IJIJ
1359 END IF
1360* |
1361* +-------------------------------------------------------------------*
1362 IF ( LRESMP ) GO TO 50
1363*
1364*--<--<--<--<--<--< go to resampling if something was wrong
1365* Here the modification to take into account properly the
1366* cascade nucleon momenta
1367* New version: rotate the spent momentum in the lab frame
1368 IF ( IGREYT .GT. 0 ) THEN
1369 PZCASC = MAX ( PZCASC, ZERZER )
1370 PTH = PXCASC * PXCASC + PYCASC * PYCASC
1371 PLA = SQRT ( PTH + PZCASC * PZCASC )
1372 PTH = SQRT ( PTH )
1373 CDE = PZCASC / PLA
1374 SDE = PTH / PLA
1375 IF (SDE .GE. ANGLGB) THEN
1376 CFE = PXCASC / PTH
1377 SFE = PYCASC / PTH
1378 ELSE
1379 CFE = 1.D+00
1380 SFE = 0.D+00
1381 END IF
1382 CALL TTRANS ( TXI, TYI, TZI, CDE, SDE, SFE, CFE,
1383 & CXXINC, CYYINC, CZZINC )
1384 PXCASC = PLA * CXXINC
1385 PYCASC = PLA * CYYINC
1386 PZCASC = PLA * CZZINC
1387 END IF
1388 PXLEFT = PXTTOT - PXCASC - TXI * TVGRE0
1389 PYLEFT = PYTTOT - PYCASC - TYI * TVGRE0
1390 PZLEFT = PZTTOT - PZCASC - TZI * TVGRE0
1391 PLA = SQRT ( PXLEFT * PXLEFT + PYLEFT * PYLEFT + PZLEFT
1392 & * PZLEFT )
1393 DEKVSP = PLA - EKIN - AM (IJ)
1394* +-------------------------------------------------------------------*
1395* | Check the momentum versus the total energy
1396 IF ( DEKVSP .GE. 0.D+00 ) THEN
1397* | +----------------------------------------------------------------*
1398* | | There are good reasons to believe that the following attempt
1399* | | is dangerous rather than useful, leading to a Dp > DEk
1400 IF ( TVGRE0 .GT. 0.D+00 ) THEN
1401 TVTENT = TVTENT - TVGRE0
1402 EKIN = EKIN + TVGRE0
1403 TVGRE0 = 0.D+00
1404 PXLEFT = PXTTOT - PXCASC
1405 PYLEFT = PYTTOT - PYCASC
1406 PZLEFT = PZTTOT - PZCASC
1407 PLA = SQRT ( PXLEFT * PXLEFT + PYLEFT * PYLEFT + PZLEFT
1408 & * PZLEFT )
1409 IF ( EKIN + AM (IJ) .GT. PLA ) GO TO 300
1410* | | This part is new (1-10-91)
1411 DEEXTR = EKE - EKIN - EINTR + IGREYP * AM (1)
1412 & + IGREYN * AM (8) - 1.5D+00
1413 & * (IGREYP+IGREYN) * EBNDAV
1414 IF ( DEEXTR .GT. 0.D+00 ) THEN
1415 EKIN = EKIN + 0.5D+00 * DEEXTR
1416 EINCP = EINCP - 0.5D+00 * DEEXTR
1417 EINCN = EINCN - 0.5D+00 * DEEXTR
1418 PXLEFT = PXTTOT - PXCASC + 0.5D+00 * TXI * DEEXTR
1419 PYLEFT = PYTTOT - PYCASC + 0.5D+00 * TYI * DEEXTR
1420 PZLEFT = PZTTOT - PZCASC + 0.5D+00 * TZI * DEEXTR
1421 PLA = SQRT ( PXLEFT * PXLEFT + PYLEFT * PYLEFT + PZLEFT
1422 & * PZLEFT )
1423 IF ( EKIN + AM (IJ) .GT. PLA ) GO TO 300
1424 END IF
1425 END IF
1426* | |
1427* | +----------------------------------------------------------------*
1428 IF ( PLA - EKIN - AM (IJ) .LE. 1.D-04 * PLA ) GO TO 300
1429* | +----------------------------------------------------------------*
1430* | | Printing condition relaxed, A.F., 23-12-92
1431 IF ( PLA - EKIN - AM (IJ) .GT. 1.D-02 * PLA ) THEN
1432 WRITE ( LUNERR,* )' Eventv: ekin+am < pla,ij,igreyt',
1433 & EKIN+AM(IJ),PLA,IJ,IGREYT
1434 END IF
1435* | |
1436* | +----------------------------------------------------------------*
1437 PTH = PLA
1438 PLA = EKIN + AM(IJ)
1439 PXLEFT = PXLEFT * PLA / PTH
1440 PYLEFT = PYLEFT * PLA / PTH
1441 PZLEFT = PZLEFT * PLA / PTH
1442 END IF
1443* |
1444* +-------------------------------------------------------------------*
1445 300 CONTINUE
1446* end new version
1447* +-------------------------------------------------------------------*
1448* | Check if we have enough energy to enter the high energy module
1449* | if not reset the accumulators and go to nucrin
1450 IF ( PLA .LT. PTHRSH (IJIJ) ) THEN
1451 PREF = 0.45D+00 * POO
1452* | +----------------------------------------------------------------*
1453* | | Check if the momentum is much smaller than the original one
1454 IF ( PLA .LE. PREF ) THEN
1455* | | +-------------------------------------------------------------*
1456* | | |
1457 IF ( PLA .LE. 0.4D+00 * POO ) THEN
1458 WRITE ( LUNERR,* )' Eventv: Pla < 0.4 Poo',PLA,POO,IJIJ,
1459 & TKIN,IMAT
1460 LRESMP = .TRUE.
1461 GO TO 50
1462 END IF
1463* | | |
1464* | | +-------------------------------------------------------------*
1465 END IF
1466* | |
1467* | +----------------------------------------------------------------*
1468 PREF = MAX ( PREF, TWOTWO )
1469 PREF = MIN ( PREF, FOUFOU )
1470* | +----------------------------------------------------------------*
1471* | | Originally it was < 10, but with Fermi momentum sometimes
1472* | | we were calling Hadrin with p > 10 GeV/c
1473* IF ( POO .LE. 9.75D+00 ) THEN
1474 IF ( ( POO .LE. 9.75D+00 .AND. IJNUCR (IJIJ) .LE. 0 )
1475 & .OR. PLA .LT. PREF ) THEN
1476 IGREYP = 0
1477 IGREYN = 0
1478 KTARP = 0
1479 KTARN = 0
1480 TVGRE0 = 0.D+00
1481 TVGREY = 0.D+00
1482 EINTR = 0.D+00
1483 PXINTR = 0.D+00
1484 PYINTR = 0.D+00
1485 PZINTR = 0.D+00
1486 ICINTR = 0
1487 IBINTR = 0
1488 ANOW = BBTAR
1489 ZNOW = ZZTAR
1490 IF ( IJ .EQ. 26 ) IJ = 23
1491 GO TO 100
1492 END IF
1493* | |
1494* | +----------------------------------------------------------------*
1495 END IF
1496* |
1497* +-------------------------------------------------------------------*
1498 TXX = PXLEFT / PLA
1499 TYY = PYLEFT / PLA
1500 TZZ = PZLEFT / PLA
1501 ZTEMP = ZZTAR - IGREYP
1502 ATEMP = BBTAR - IGREYN
1503 ERES = ETTOT - EKIN - AM (IJ) - EINTR
1504 AMNRES = AMUAMU * ATEMP + 1.D-03 * FKENER ( ATEMP, ZTEMP ) -
1505 & ZTEMP * AMELEC + ELBNDE ( NINT (ZTEMP) )
1506C
1507C********************************************************************
1508C FOR MOMENTA ABOVE 5.0 GEV/C USE NUCEVT
1509C********************************************************************
1510C
1511*
1512* From here the high energy model....
1513*
1514 NNHAD = 0
1515 CALL NUCEVV ( NNHAD, IJ, PLA, EKIN, TXX, TYY, TZZ )
1516 IF ( LRESMP ) GO TO 50
1517*--<--<--<--<--< go to resampling if something was wrong
1518* +-------------------------------------------------------------------*
1519* |
1520 IF ( LLASTN .AND. NINT ( ANOW ) .EQ. 1 ) THEN
1521* | +----------------------------------------------------------------*
1522* | |
1523 IF ( KTINC .EQ. 1 ) THEN
1524 IGREYP = IGREYP + 1
1525 EINCP = EINCP + EKINC
1526 ZNOW = ZNOW - 1.D+00
1527* | |
1528* | +----------------------------------------------------------------*
1529* | |
1530 ELSE
1531 IGREYN = IGREYN + 1
1532 EINCN = EINCN + EKINC
1533 END IF
1534* | |
1535* | +----------------------------------------------------------------*
1536 ANOW = ANOW - 1.D+00
1537 NP = NP + 1
1538 IF (NP .GT. MXP) GO TO 400
1539 KPART(NP)= KPTOIP (KTINC)
1540 WEI (NP) = WE
1541 TKI (NP) = EKINC
1542* | Updating conservation counters
1543 EINTR = EINTR + TKI (NP) + AMINC
1544 PXINTR = PXINTR + PXXINC
1545 PYINTR = PYINTR + PYYINC
1546 PZINTR = PZINTR + PZZINC
1547 ICINTR = ICINTR + ICH (KTINC)
1548 IBINTR = IBINTR + IBAR (KTINC)
1549 PLR (NP) = SQRT (PXXINC * PXXINC + PYYINC * PYYINC +
1550 & PZZINC * PZZINC)
1551 CXR (NP) = PXXINC / PLR (NP)
1552 CYR (NP) = PYYINC / PLR (NP)
1553 CZR (NP) = PZZINC / PLR (NP)
1554 END IF
1555* |
1556* +-------------------------------------------------------------------*
1557
1558* +-------------------------------------------------------------------*
1559* | Loop over particles produced in nucevt
1560 DO 301 I=1,NNHAD
1561 NP = NP+1
1562 IF (NP .GT. MXP) GO TO 400
1563* | Get the "paprop" index for Nucevt produced particles
1564 KPART(NP) = KPTOIP (NRENU(I))
1565 IF ( KPART (NP) .EQ. 0 ) THEN
1566 WRITE (LUNERR,*)' **** Charmed particle produced in Nucevv',
1567 & NRENU(I),HEPNU(I),AMNU(I),' ****'
1568 KPART (NP) = NRENU (I)
1569 END IF
1570 PLR (NP) = SQRT ( PXNU (I)**2 + PYNU (I)**2 + PZNU (I)**2 )
1571 CXR (NP) = PXNU (I) / PLR (NP)
1572 CYR (NP) = PYNU (I) / PLR (NP)
1573 CZR (NP) = PZNU (I) / PLR (NP)
1574 TKI(NP) = HEPNU(I) - AMNU(I)
1575 WEI(NP) = WE
1576* | +----------------------------------------------------------------*
1577* | |
1578 IF (TKI(NP) .LE. 0.D+00) THEN
1579* | | Is this the only check on the generated particle energy??
1580 EUZ = EUZ - HEPNU(I)
1581 PUX = PUX - PXNU(I)
1582 PUY = PUY - PYNU(I)
1583 PUZ = PUZ - PZNU(I)
1584 ICU = ICU - ICHNU(I)
1585 IBU = IBU - IBARNU(I)
1586 WRITE (LUNERR,*) ' Eventq: Kin en. < 0 from nucevt',TKI(NP)
1587 NP = NP - 1
1588 END IF
1589* | |
1590* | +----------------------------------------------------------------*
1591 301 CONTINUE
1592* |
1593* +-------------------------------------------------------------------*
1594 500 CONTINUE
1595*
1596* Now we have to compute the excitation energy: we have used Eincp and
1597* Eincn for cascade protons and neutrons, Tvgre0 and Tvgrey have been
1598* generated by cascade nucleons under threshold, Tveuz has been gene-
1599* rated by the high energy collisions: however Tveuz and Tvgrey are
1600* only approximate since they have been computed starting from an
1601* average binding energy, furthermore Tvgrey is already accounted for
1602* inside Eincp and Eincn. The actual energy spent inside high energy
1603* collisions was Eke - Eincp - Eincn - Tvgre0 + Efrm, the one in ca-
1604* scade nucleons (approximately) Eincp + Eincn - Tvgrey: so first
1605* check the energy balance without excitation energy
1606* and then compute Tv!!!
1607*
1608* Now the balance!!!!!!!!!!!!!!!
1609*
1610 ERES = ETTOT - EUZ - ENUCR - EINTR
1611 PXRES = PXTTOT - PUX - PXNUCR - PXINTR
1612 PYRES = PYTTOT - PUY - PYNUCR - PYINTR
1613 PZRES = PZTTOT - PUZ - PZNUCR - PZINTR
1614 ICRES = ICHTOT - ICU - ICNUCR - ICINTR
1615 IBRES = IBTOT - IBU - IBNUCR - IBINTR
1616* +-------------------------------------------------------------------*
1617* |
1618 IF ( NINT (ZNOW) .NE. ICRES .OR. NINT (ANOW) .NE. IBRES ) THEN
1619* | +----------------------------------------------------------------*
1620* | |
1621 IF ( LNUCRI ) THEN
1622 WRITE (LUNERR,*)' Eventq: charge/baryon conservation',
1623 & ' failure with Nucrin',
1624 & NINT (ZNOW), ICRES, NINT (ANOW), IBRES
1625* | |
1626* | +----------------------------------------------------------------*
1627* | |
1628 ELSE
1629 WRITE (LUNERR,*)' Eventq: charge/baryon conservation',
1630 & ' failure with Nucevt',
1631 & NINT (ZNOW), ICRES, NINT (ANOW), IBRES
1632 END IF
1633* | |
1634* | +----------------------------------------------------------------*
1635 LRESMP = .TRUE.
1636 GO TO 50
1637* |--<--<--<--<--< go to resampling
1638 END IF
1639* |
1640* +-------------------------------------------------------------------*
1641* +-------------------------------------------------------------------*
1642* |
1643 IF ( IBRES .GT. 0 ) THEN
1644 AMMRES = ANOW * AMUAMU + 0.001D+00 * FKENER ( ANOW, ZNOW )
1645 AMNRES = AMMRES - ZNOW * AMELEC + ELBNDE ( ICRES )
1646* AMNRES = ANOW * AMUC12
1647 EKR0 = ERES - AMNRES
1648* | Now switch to atomic masses:
1649 ERES = ERES + AMMTAR - AMNTAR - ( ZZTAR - ZNOW ) * AMELEC
1650 EKRES = ERES - AMMRES
1651 TVTENT = TVGRE0 + TVGREY + TVEUZ
1652 TVCOMP = ERES - ANOW * AMUC12
1653 IF ( LNUCRI ) TVCOMP = TVCOMP + ( AMNTAR - BBTAR * AMUC12 )
1654* |
1655* +-------------------------------------------------------------------*
1656* |
1657 ELSE
1658 AMMRES = 0.D+00
1659 AMNRES = 0.D+00
1660 ERES = 0.D+00
1661 EKR0 = 0.D+00
1662 EKRES = 0.D+00
1663 TVTENT = 0.D+00
1664 GO TO 600
1665 END IF
1666* |
1667* +-------------------------------------------------------------------*
1668* Check that the kinetic energy of the residual nuclei is not much
1669* different from our prevision and kinematically consistent with
1670* its momentum
1671* +-------------------------------------------------------------------*
1672* |
1673 IF ( EKRES .LE. 0.D+00 ) THEN
1674 WRITE ( LUNERR,* )' Eventq: negative kinetic energy for',
1675 & ' the residual nucleus!!',ICRES,IBRES,
1676 & REAL (EKRES), LNUCRI
1677* | +----------------------------------------------------------------*
1678* | |
1679 IF ( EKRES .LT. -3.D-3 ) THEN
1680 LRESMP = .TRUE.
1681 GO TO 50
1682*--|--|--<--<--<--<--< go to resampling
1683 END IF
1684* | |
1685* | +----------------------------------------------------------------*
1686 EKRES = 0.D+00
1687 TVRECL = 0.D+00
1688 AMSTAR = AMMRES
1689 TVCMS = 0.D+00
1690 PTRES2 = 0.D+00
1691 PXRES = 0.D+00
1692 PYRES = 0.D+00
1693 PZRES = 0.D+00
1694* |
1695* +-------------------------------------------------------------------*
1696* |
1697 ELSE
1698 PTRES2 = PXRES**2 + PYRES**2 + PZRES**2
1699 AMSTAR = ERES**2 - PTRES2
1700* | +----------------------------------------------------------------*
1701* | |
1702 IF ( AMSTAR .GE. AMMRES**2 ) THEN
1703 AMSTAR = SQRT ( AMSTAR )
1704 TVCMS = AMSTAR - AMMRES
1705* | |
1706* | +----------------------------------------------------------------*
1707* | | If the following condition is satisfied it is only
1708* | | a rounding problem, set the excitation energy to 0
1709* | | and continue
1710 ELSE IF ( AMMRES**2 - AMSTAR .LT. 2.D+00 * AMSTAR * TVEPSI
1711 & ) THEN
1712 AMSTAR = AMMRES
1713 TVCMS = 0.D+00
1714* | |
1715* | +----------------------------------------------------------------*
1716* | |
1717 ELSE IF ( AMSTAR .LE. 0.D+00 ) THEN
1718 WRITE ( LUNERR,* )' Eventq: immaginary mass for',
1719 & ' the residual nucleus!!',ICRES,IBRES,
1720 & REAL (AMSTAR)
1721 LRESMP = .TRUE.
1722 GO TO 50
1723*--|--|--<--<--<--<--< go to resampling
1724* AMSTAR = AMMRES
1725* TVCMS = 0.D+00
1726* | |
1727* | +----------------------------------------------------------------*
1728* | |
1729 ELSE
1730 AMSTAR = SQRT ( AMSTAR )
1731* | | +-------------------------------------------------------------*
1732* | | | If the following condition is satisfied it is only
1733* | | | a rounding problem, set the excitation energy to 0
1734* | | | and continue
1735 IF ( AMMRES - AMSTAR .LT. TVEPSI ) THEN
1736 AMSTAR = AMMRES
1737 TVCMS = 0.D+00
1738 TVRECL = ERES - AMSTAR
1739 GO TO 550
1740 END IF
1741* | | |
1742* | | +-------------------------------------------------------------*
1743 WRITE ( LUNERR,* )' Eventq: negative excitation energy for',
1744 & ' the residual nucleus!!',ICRES,IBRES,
1745 & REAL ( AMSTAR - AMMRES ), LNUCRI
1746* | | +-------------------------------------------------------------*
1747* | | |
1748 IF ( AMSTAR - AMMRES .LT. -3.D-3 ) THEN
1749 LRESMP = .TRUE.
1750 GO TO 50
1751*--|--|--|--<--<--<--<--< go to resampling
1752 END IF
1753* | | |
1754* | | +-------------------------------------------------------------*
1755 AMSTAR = AMMRES
1756 TVCMS = 0.D+00
1757 AMSTAR = AMMRES
1758 TVCMS = 0.D+00
1759 HELP = SQRT ( ( ERES - AMMRES ) * ( ERES + AMMRES )
1760 & / PTRES2 )
1761 PXRES = PXRES * HELP
1762 PYRES = PYRES * HELP
1763 PZRES = PZRES * HELP
1764 PTRES2 = PTRES2 * HELP * HELP
1765 END IF
1766* | |
1767* | +----------------------------------------------------------------*
1768 TVRECL = ERES - AMSTAR
1769 END IF
1770* |
1771* +-------------------------------------------------------------------*
1772 550 CONTINUE
1773* The following two cards are equivalent providing the kinematical
1774* limits are ok and we use for both amnres or ammres! Now Tv is left
1775* = 0
1776 TV = 0.D+00
1777* TV = EKRES
1778* TV = TVRECL + TVCMS
1779* +-------------------------------------------------------------------*
1780* |
1781 IF ( .NOT. LEVPRT .AND. ABS ( ( TVTENT - TVCOMP ) / TVCOMP )
1782 & .GT. 30000000.D+00 ) THEN
1783* | +----------------------------------------------------------------*
1784* | |
1785 IF ( LINCTV ) THEN
1786 WRITE ( LUNERR,* )
1787 & ' Eventq: excitation energy very different',
1788 & ' from the approximate one!!', ICRES, IBRES,
1789 & REAL (TVTENT), REAL (TVCOMP), LNUCRI
1790 END IF
1791* | |
1792* | +----------------------------------------------------------------*
1793 END IF
1794* |
1795* +-------------------------------------------------------------------*
1796 EKRES = TVRECL
1797 600 CONTINUE
1798 EOTEST = EOTEST - EUZ - ENUCR - EINTR - EKR0 - AMNRES
1799* +-------------------------------------------------------------------*
1800* |
1801 IF ( ABS (EOTEST) .GT. ETEPS ) THEN
1802* | +----------------------------------------------------------------*
1803* | |
1804 IF ( LNUCRI ) THEN
1805 WRITE (LUNERR,*)' Eventq: eotest failure with Nucrin',
1806 & EOTEST
1807* | |
1808* | +----------------------------------------------------------------*
1809* | |
1810 ELSE
1811 WRITE (LUNERR,*)' Eventq: eotest failure with Nucevt',
1812 & EOTEST
1813 END IF
1814* | |
1815* | +----------------------------------------------------------------*
1816 LRESMP = .TRUE.
1817 GO TO 50
1818* |--|--<--<--<--<--< go to resampling
1819 END IF
1820* |
1821* +-------------------------------------------------------------------*
1822 IF ( IBRES .EQ. 0 ) RETURN
1823*-->-->-->-->--> Here for the evaporation step!!!
1824* +-------------------------------------------------------------------*
1825* | Check if the evaporation is requested
1826 IF ( LEVPRT ) THEN
1827 PTRES = SQRT ( PTRES2 )
1828 CALL EVEVAP ( WE )
1829 IF ( LRESMP ) GO TO 50
1830*
1831*--|--<--<--<--<--<--< go to resampling
1832* |
1833* +-------------------------------------------------------------------*
1834* | No evaporation
1835 ELSE
1836 TVHEAV = 0.D+00
1837 END IF
1838* |
1839* +-------------------------------------------------------------------*
1840 IF (IPRI.NE.1) RETURN
1841C
1842C********************************************************************
1843C TEST PRINTOUT
1844C********************************************************************
1845C
1846 WRITE(LUNOUT,590)NP-NP0,NNHAD,IJ,IMAT,POO,EKE,TXI,TYI,TZI,WE
1847 590 FORMAT (4I7,6F12.6)
1848 DO 501 I=NP0+1,NP
1849 WRITE(LUNOUT,591)I,KPART(I),CXR(I),CYR(I),CZR (I),TKI(I),PLR (I),
1850 * WEI(I)
1851 591 FORMAT (2I5,6F12.6)
1852 501 CONTINUE
1853 RETURN
1854C
1855C********************************************************************
1856C FINUC FLOWS OVER - THIS IS FATAL - INCREASE THE SIZE OF IT
1857C********************************************************************
1858C
1859 400 CONTINUE
1860 WRITE(LUNOUT,490)
1861 490 FORMAT(1X,'OVERFLOW IN EVENTQ - INCREASE THE SIZE OF THE',
1862 1' COMMON BLOCK FINUC.')
1863 STOP
1864 END