This commit was generated by cvs2svn to compensate for changes in r2,
[u/mrichter/AliRoot.git] / GEANT321 / gheisha / nbanh.F
1 *
2 * $Id$
3 *
4 * $Log$
5 * Revision 1.1.1.1  1995/10/24 10:20:58  cernlib
6 * Geant
7 *
8 *
9 #include "geant321/pilot.h"
10 *CMZ :  3.21/02 29/03/94  15.41.38  by  S.Giani
11 *-- Author :
12       SUBROUTINE NBANH(NOPT)
13 C
14 C *** ANTI NEUTRON ANNIHILATION AT REST ***
15 C *** NVE 04-MAR-1988 CERN GENEVA ***
16 C
17 C ORIGIN : H.FESEFELDT (09-JULY-1987)
18 C
19 C NOPT=0    NO ANNIHILATION
20 C NOPT=1    ANNIH.IN PI+ PI-
21 C NOPT=2    ANNIH.IN PI0 PI0
22 C NOPT=3    ANNIH.IN PI+ PI0
23 C NOPT=4    ANNIH.IN GAMMA GAMMA
24 C
25 #include "geant321/s_defcom.inc"
26 C
27       DIMENSION BRR(3)
28       DIMENSION RNDM(3)
29       DATA BRR/0.125,0.25,0.50/
30 C
31       PV(1,1)=0.
32       PV(2,1)=0.
33       PV(3,1)=0.
34       PV(4,1)=ABS(RMASS(17))
35       PV(5,1)=RMASS(17)
36       PV(6,1)=0.
37       PV(7,1)=TOF
38       PV(8,1)=IPART
39       PV(9,1)=0.
40       PV(10,1)=USERW
41       IER(85)=IER(85)+1
42       ISW=1
43       CALL GRNDM(RNDM,1)
44       RAN=RNDM(1)
45       IF(RAN.GT.BRR(1)) ISW=2
46       IF(RAN.GT.BRR(2)) ISW=3
47       IF(RAN.GT.BRR(3)) ISW=4
48       NOPT=ISW
49 C**
50 C**  EVAPORATION
51 C**
52       CALL COMPO
53       RMNVE1=RMASS(7)
54       RMNVE2=RMASS(9)
55       IF (ISW .EQ. 2) RMNVE1=RMASS(8)
56       IF (ISW .EQ. 2) RMNVE2=RMASS(8)
57       IF (ISW .EQ. 3) RMNVE2=RMASS(8)
58       IF (ISW .EQ. 4) RMNVE1=RMASS(1)
59       IF (ISW .EQ. 4) RMNVE2=RMASS(1)
60       EK=RMASS(16)+ABS(RMASS(17))-RMNVE1-RMNVE2
61       TKIN=EXNU(EK)
62       EK=EK-TKIN
63       IF(EK.LT.0.0001) EK=0.0001
64       EK=EK/2.
65       EN=EK+(RMNVE1+RMNVE2)/2.0
66       S=AMAS*AMAS+RMASS(16)**2+2.0*RMASS(16)*EN
67       RS=SQRT(S)
68       PCM=SQRT(ABS(EN*EN-RMNVE1*RMNVE2))
69       CALL GRNDM(RNDM,2)
70       PHI=2.*PI*RNDM(1)
71       COST=-1.+2.*RNDM(2)
72       SINT=SQRT(ABS(1.-COST*COST))
73       PV(1,2)=PCM*SINT*COS(PHI)
74       PV(2,2)=PCM*SINT*SIN(PHI)
75       PV(3,2)=PCM*COST
76       DO 1 I=1,3
77     1 PV(I,3)=-PV(I,2)
78       PV(5,2)=RMNVE1
79       PV(5,3)=RMNVE2
80       IF(ISW.LE.3) GOTO 2
81       PV(5,2)=0.
82       PV(5,3)=0.
83     2 PV(4,2)=SQRT(PV(5,2)**2+PCM**2)
84       PV(4,3)=SQRT(PV(5,3)**2+PCM**2)
85       PV(7,2)=TOF
86       PV(7,3)=TOF
87       PV(9,2)=0.
88       PV(9,3)=0.
89       PV(10,2)=0.
90       PV(10,3)=0.
91       GOTO (21,22,23,24), ISW
92    21 PV(6,2)=1.
93       PV(6,3)=-1.
94       PV(8,2)=7.
95       PV(8,3)=9.
96       GOTO 25
97    22 PV(6,2)=0.
98       PV(6,3)=0.
99       PV(8,2)=8.
100       PV(8,3)=8.
101       GOTO 25
102    23 PV(6,2)=1.
103       PV(6,3)=0.
104       PV(8,2)=7.
105       PV(8,3)=8.
106       GOTO 25
107    24 PV(6,2)=0.
108       PV(6,3)=0.
109       PV(8,2)=1.
110       PV(8,3)=1.
111    25 NT=3
112       IF(ATNO2.LT.1.5) GOTO 40
113       AFC=0.312+0.200*LOG(LOG(S))+S**1.5/6000.
114 C     TARG=AFC*(ATNO2**0.33 -1.0)
115       CFA=0.025*((ATNO2-1.)/120.)*EXP(-(ATNO2-1.)/120.)
116       TARG=1.
117       TEX=ENP(1)
118       IF(TEX.LT.0.001) GOTO 445
119       BLACK=(1.5+1.25*TARG)*ENP(1)/(ENP(1)+ENP(3))
120       CALL POISSO(BLACK,NBL)
121       IF(IFIX(TARG)+NBL.GT.ATNO2) NBL=ATNO2-TARG
122       IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
123       IF(NBL.LE.0) GOTO 445
124       EKIN=TEX/NBL
125       EKIN2=0.
126       CALL STEEP(XX)
127       DO 441 I=1,NBL
128       IF(NT.EQ.MXGKPV-2) GOTO 441
129       IF(EKIN2.GT.TEX) GOTO 443
130       CALL GRNDM(RNDM,1)
131       RAN1=RNDM(1)
132       CALL NORMAL(RAN2)
133       EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
134       IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1)
135       EKIN1=EKIN1*XX
136       EKIN2=EKIN2+EKIN1
137       IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
138       IF(EKIN1.LT.0.) EKIN1=0.001
139       IPA1=16
140       PNRAT=1.-ZNO2/ATNO2
141       CALL GRNDM(RNDM,3)
142       IF(RNDM(1).GT.PNRAT) IPA1=14
143       NT=NT+1
144       COST=-1.+RNDM(2)*2.
145       SINT=SQRT(ABS(1.-COST*COST))
146       PHI=TWPI*RNDM(3)
147       IPA(NT)=-IPA1
148       PV(5,NT)=ABS(RMASS(IPA1))
149       PV(6,NT)=RCHARG(IPA1)
150       PV(7,NT)=TOF
151       PV(8,NT)=IPA1
152       PV(9,NT)=0.
153       PV(10,NT)=0.
154       PV(4,NT)=EKIN1+PV(5,NT)
155       PP=SQRT(ABS(PV(4,NT)**2-PV(5,NT)**2))
156       PV(1,NT)=PP*SINT*SIN(PHI)
157       PV(2,NT)=PP*SINT*COS(PHI)
158       PV(3,NT)=PP*COST
159   441 CONTINUE
160   443 IF(ATNO2.LT.230.) GOTO 445
161       IF(EK.GT.2.0) GOTO 445
162       II=NT+1
163       KK=0
164       EKA=EK
165       IF(EKA.GT.1.) EKA=EKA*EKA
166       IF(EKA.LT.0.1) EKA=0.1
167       IKA=IFIX(3.6/EKA)
168       DO 444 I=1,NT
169       II=II-1
170       IF(IPA(II).NE.-14) GOTO 444
171       IPA(II)=-16
172       IPA1  = 16
173       PV(5,II)=ABS(RMASS(IPA1))
174       PV(6,II)=RCHARG(IPA1)
175       PV(8,II)=IPA1
176       KK=KK+1
177       IF(KK.GT.IKA) GOTO 445
178   444 CONTINUE
179 C**
180 C** THEN ALSO DEUTERONS, TRITONS AND ALPHAS
181 C**
182   445 TEX=ENP(3)
183       IF(TEX.LT.0.001) GOTO 40
184       BLACK=(1.5+1.25*TARG)*ENP(3)/(ENP(1)+ENP(3))
185       CALL POISSO(BLACK,NBL)
186       IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
187       IF(NBL.LE.0) GOTO 40
188       EKIN=TEX/NBL
189       EKIN2=0.
190       CALL STEEP(XX)
191       DO 442 I=1,NBL
192       IF(NT.EQ.MXGKPV-2) GOTO 442
193       IF(EKIN2.GT.TEX) GOTO 40
194       CALL GRNDM(RNDM,1)
195       RAN1=RNDM(1)
196       CALL NORMAL(RAN2)
197       EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
198       IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1)
199       EKIN1=EKIN1*XX
200       EKIN2=EKIN2+EKIN1
201       IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
202       IF(EKIN1.LT.0.) EKIN1=0.001
203       CALL GRNDM(RNDM,3)
204       COST=-1.+RNDM(1)*2.
205       SINT=SQRT(ABS(1.-COST*COST))
206       PHI=TWPI*RNDM(2)
207       RAN=RNDM(3)
208       IPA(NT+1)=-30
209       IF(RAN.GT.0.60) IPA(NT+1)=-31
210       IF(RAN.GT.0.90) IPA(NT+1)=-32
211       INVE=ABS(IPA(NT+1))
212       PV(5,NT+1)=RMASS(INVE)
213       NT=NT+1
214       PV(6,NT)=RCHARG(INVE)
215       PV(7,NT)=TOF
216       PV(8,NT)=ABS(IPA(NT))
217       PV(9,NT)=0.
218       PV(10,NT)=0.
219       PV(4,NT)=PV(5,NT)+EKIN1
220       PP=SQRT(ABS(PV(4,NT)**2-PV(5,NT)**2))
221       PV(1,NT)=PP*SINT*SIN(PHI)
222       PV(2,NT)=PP*SINT*COS(PHI)
223       PV(3,NT)=PP*COST
224   442 CONTINUE
225    40 INTCT=INTCT+1.
226       CALL SETCUR(2)
227       NTK=NTK+1
228       IF(NT.EQ.2) GO TO 9999
229       DO 50 I=3,NT
230       IF(NTOT.LT.NSIZE/12) GOTO 43
231       GO TO 9999
232    43 CALL SETTRK(I)
233    50 CONTINUE
234       CALL LENGTX(3,PP)
235       IF(NPRT(3))
236      *WRITE(NEWBCD,1001) XEND,YEND,ZEND,P,NCH,PP,PV(6,3)
237 1001  FORMAT(1H0,'NB ANNIHILATION AT REST  POSITION',3(1X,F8.2),1X,
238      * 'PI MOMENTA,CHARGES',2(1X,F8.4,1X,F4.1))
239 C
240  9999 CONTINUE
241       END