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