]> git.uio.no Git - u/mrichter/AliRoot.git/blame - GEANT321/gheisha/nbanh.F
README file from R.Barbera
[u/mrichter/AliRoot.git] / GEANT321 / gheisha / nbanh.F
CommitLineData
fe4da5cc 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)
13C
14C *** ANTI NEUTRON ANNIHILATION AT REST ***
15C *** NVE 04-MAR-1988 CERN GENEVA ***
16C
17C ORIGIN : H.FESEFELDT (09-JULY-1987)
18C
19C NOPT=0 NO ANNIHILATION
20C NOPT=1 ANNIH.IN PI+ PI-
21C NOPT=2 ANNIH.IN PI0 PI0
22C NOPT=3 ANNIH.IN PI+ PI0
23C NOPT=4 ANNIH.IN GAMMA GAMMA
24C
25#include "geant321/s_defcom.inc"
26C
27 DIMENSION BRR(3)
28 DIMENSION RNDM(3)
29 DATA BRR/0.125,0.25,0.50/
30C
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
49C**
50C** EVAPORATION
51C**
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.
114C 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
179C**
180C** THEN ALSO DEUTERONS, TRITONS AND ALPHAS
181C**
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)
2371001 FORMAT(1H0,'NB ANNIHILATION AT REST POSITION',3(1X,F8.2),1X,
238 * 'PI MOMENTA,CHARGES',2(1X,F8.4,1X,F4.1))
239C
240 9999 CONTINUE
241 END