GEANT321/gheisha/twob.F

   1 *
   2 * $Id$
   3 *
   4 * $Log$
   5 * Revision 1.1.1.1  1995/10/24 10:20:59  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 TWOB(IPPP,NFL,AVERN)
  13 C
  14 C *** GENERATION OF MOMENTA FOR ELAST. AND QUASI ELAST. 2 BODY REACT. ***
  15 C *** NVE 04-MAY-1988 CERN GENEVA ***
  16 C
  17 C ORIGIN : H.FESEFELDT 15-SEP-1987
  18 C
  19 C THE SIMPLE FORMULA DS/D|T| = S0* EXP(-B*|T|) IS USED.
  20 C THE B VALUES ARE PARAMETRIZATIONS FROM EXPERIMENTAL DATA .
  21 C NOT AVAILABLE VALUES ARE TAKEN FROM THOSE OF SIMILAR REACTIONS
  22 C
  23 #include "geant321/s_defcom.inc"
  24       DIMENSION RNDM(3)
  25 C
  26 C     DATA CB/3./
  27       DATA CB/0.01/
  28 C
  29 C --- STATEMENT FUNCTIONS ---
  30       BPP(X)=4.225+1.795*LOG(X)
  31 C
  32 C**
  33 C**  FOR DIFFRACTION SCATTERING ON HEAVY NUCLEI USE BETTER ROUTINE
  34 C**  "COSCAT"
  35 C
  36       TARMAS=RMASS(14)
  37       IF (NFL .EQ. 2) TARMAS=RMASS(16)
  38       ENP(8)=RMASS(IPPP)**2+TARMAS**2+2.0*TARMAS*ENP(6)
  39       ENP(9)=SQRT(ENP(8))
  40       EK=ENP(5)
  41       EN=ENP(6)
  42       P=ENP(7)
  43       S=ENP(8)
  44       RS=ENP(9)
  45       CFA=0.025*((ATNO2-1.)/120.)*EXP(-(ATNO2-1.)/120.)
  46       IF(ATNO2.LT.1.5) GOTO 500
  47       IPA1=ABS(IPA(1))
  48       IPA2=ABS(IPA(2))
  49       RMC=RMASS(IPA1)
  50       RMD=RMASS(IPA2)
  51       RCHC=RCHARG(IPA1)
  52       RCHD=RCHARG(IPA2)
  53       IF(ABS(RMC-AMAS).GT.0.001) GOTO 500
  54       RMNVE=RMASS(14)
  55       IF (NFL .EQ. 2) RMNVE=RMASS(16)
  56       IF(ABS(RMD-RMNVE).GT.0.001) GOTO 500
  57       IF(ABS(RCHC-NCH).GT.0.5) GOTO 500
  58       IF(NFL.EQ.1.AND.RCHD.LT.0.5) GOTO 500
  59       IF(NFL.EQ.2.AND.ABS(RCHD).GT.0.5) GOTO 500
  60       IF(ENP(1).GT.0.0001.OR.ENP(3).GT.0.0001) GOTO 500
  61       CALL COSCAT
  62       GO TO 9999
  63 C**
  64 C**  SET EFFECTIVE 4-MOMENTUM OF INITIAL PARTICLE
  65 C**
  66   500 PV( 1,MXGKPV-1)=P*PX
  67       PV( 2,MXGKPV-1)=P*PY
  68       PV( 3,MXGKPV-1)=P*PZ
  69       PV( 4,MXGKPV-1)=EN
  70       PV( 5,MXGKPV-1)=AMAS
  71       PV( 6,MXGKPV-1)=NCH
  72       PV( 7,MXGKPV-1)=TOF
  73       PV( 8,MXGKPV-1)=IPART
  74       PV( 9,MXGKPV-1)=0.
  75       PV(10,MXGKPV-1)=USERW
  76       IER(47)=IER(47)+1
  77       IF(NPRT(4))
  78      $  WRITE(NEWBCD,4001) (PV(J,MXGKPV-1),J=1,10),IPA(1),IPA(2)
  79       DO 2 J=1,6
  80     2 PV(J,1)=PV(J,MXGKPV-1)
  81       PV(7,1)=1.
  82       IF(PV(5,1).LT.0.) PV(7,1)=-1.
  83       PV(5,1)=ABS(PV(5,1))
  84       NT=1
  85 C**
  86 C** TWO-BODY SCATTERING POSSIBLE?? IF NOT, CONTINUE WITH ORIGINAL
  87 C** PARTICLE, BUT SPEND THE NUCLEAR EVAPORATION ENERGY
  88 C**
  89       IF(P.LT.0.1) GOTO 200
  90       IF(RS.LT.0.01) GOTO 200
  91 C**
  92 C** CALCULATE SLOPE B FOR ELASTIC SCATTERING ON PROTON/NEUTRON
  93 C**
  94       B=BPP(P)
  95       IF(B.LT.CB) B=CB
  96       IF(ABS(IPA(2)).GT.13) GOTO 9
  97       IPA(2)=14
  98       CALL GRNDM(RNDM,1)
  99       IF(RNDM(1).LT.0.5) IPA(2)=16
 100 C**
 101 C** SET MASSES AND MOMENTA FOR FINAL STATE PARTICLES
 102 C**
 103     9 RMC=RMASS(ABS(IPA(1)))
 104       RMD=RMASS(ABS(IPA(2)))
 105       PV(6,1)=RCHARG(ABS(IPA(1)))
 106       PV(6,2)=RCHARG(ABS(IPA(2)))
 107       PF=(S+RMD*RMD-RMC*RMC)**2 - 4*S*RMD*RMD
 108       IF(NPRT(4)) WRITE(NEWBCD,4002) RMC,RMD,PV(6,1),PV(6,2),RS,S,PF
 109       IF(PF.LT.0.001) GO TO 9999
 110       PF=SQRT(PF)/(2.*RS)
 111 C**
 112 C** SET BEAM AND TARGET IN CMS
 113 C**
 114       PV(1,3)=0.
 115       PV(2,3)=0.
 116       PV(3,3)=P
 117       PV(5,3)=ABS(AMAS)
 118       PV(4,3)=SQRT(P*P+AMAS*AMAS)
 119       PV(1,4)=0.
 120       PV(2,4)=0.
 121       PV(3,4)=0.
 122       RMNVE=RMASS(14)
 123       IF (NFL .EQ. 2) RMNVE=RMASS(16)
 124       PV(4,4)=RMNVE
 125       PV(5,4)=RMNVE
 126 C**
 127 C** TRANSFORM INTO CMS.
 128 C**
 129       CALL ADD(3,4,10)
 130       CALL LOR(3,10,3)
 131       CALL LOR(4,10,4)
 132 C**
 133 C** SET FINAL STATE MASSES AND ENERGIES IN CMS
 134 C**
 135       PV(5,1)=ABS(RMC)
 136       PV(5,2)=ABS(RMD)
 137       PV(7,1)=1.
 138       PV(7,2)=1.
 139       IF(RMC.LT.0.) PV(7,1)=-1.
 140       IF(RMD.LT.0.) PV(7,2)=-1.
 141       PV(4,1)=SQRT(PF*PF+PV(5,1)*PV(5,1))
 142       PV(4,2)=SQRT(PF*PF+PV(5,2)*PV(5,2))
 143 C**
 144 C** SET |T| AND |TMIN|
 145 C**
 146       CALL GRNDM(RNDM,2)
 147       CALL LENGTX(3,PIN)
 148       BTRANG=B*4.*PIN*PF
 149 C**
 150 C** SIMPLY A PROTECTION AGAINST EXPONENT OVERFLOW 1.E20 IS BIG ENOUGH
 151 C**
 152       EXINDT=-1.
 153       IF(BTRANG.LT.46) EXINDT=EXINDT+EXP(-BTRANG)
 154       TDN=LOG(1.+RNDM(1)*EXINDT)/BTRANG
 155 C**
 156 C** CACULATE (SIN(TETA/2.)**2 AND COS(TETA), SET AZIMUTH ANGLE PHI
 157 C**
 158       CTET=1.+2.*TDN
 159       IF(ABS(CTET).GT.1.) CTET=SIGN(1.,CTET)
 160       STET=SQRT((1.-CTET)*(1.+CTET))
 161       PHI=RNDM(2)*TWPI
 162 C**
 163 C** CALCULATE FINAL STATE MOMENTA IN CMS
 164 C**
 165       PV(1,1)=PF*STET*SIN(PHI)
 166       PV(2,1)=PF*STET*COS(PHI)
 167       PV(3,1)=PF*CTET
 168       PV(1,2)=-PV(1,1)
 169       PV(2,2)=-PV(2,1)
 170       PV(3,2)=-PV(3,1)
 171 C**
 172 C** TRANSFORM INTO LAB
 173 C**
 174       DO 11 I=1,2
 175       CALL LOR(I,4,I)
 176       CALL DEFS1(I,MXGKPV-1,I)
 177       IF(ATNO2.LT.1.5) GOTO 11
 178       CALL LENGTX(I,PP)
 179       IF(PP.LT.0.001) GOTO 11
 180       EKIN=PV(4,I)-ABS(PV(5,I))
 181       CALL NORMAL(RAN)
 182       EKIN=EKIN-CFA*(1.+0.5*RAN)
 183       IF(EKIN.LT.0.0001) EKIN=0.0001
 184       PP1=SQRT(EKIN*(EKIN+2.*ABS(PV(5,I))))
 185       PV(4,I)=EKIN+ABS(PV(5,I))
 186       PV(1,I)=PV(1,I)*PP1/PP
 187       PV(2,I)=PV(2,I)*PP1/PP
 188       PV(3,I)=PV(3,I)*PP1/PP
 189    11 CONTINUE
 190       NT=2
 191 C**
 192 C** ADD BLACK TRACK PARTICLES .
 193 C** HERE THE PROCEDURE IS SOMEWHAT DIFFERENT AS IN 'TWOCLU' AND 'GENXPT'
 194 C** THE REASON IS, THAT WE HAVE TO SIMULATE ALSO THE NUCLEAR REACTIONS
 195 C** AT LOW ENERGIES LIKE A(H,P)B, A(H,P P)B, A(H,N)B E.T.C.
 196 C**
 197   200 IF(ENP(1).LE.0.0001.AND.ENP(3).LE.0.0001) GOTO 40
 198       SPALL=0.
 199       TEX=ENP(1)
 200       IF(TEX.LT.0.0001) GOTO 445
 201       BLACK=TEX/0.02
 202       CALL POISSO(BLACK,NBL)
 203       IF(NBL.GT.ATNO2) NBL=ATNO2
 204       IF(ENP(1).GT.0.0001.AND.NBL.LE.0) NBL=1
 205       IF (NPRT(4)) WRITE(NEWBCD,3003) NBL,TEX
 206       IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
 207       IF(NBL.LE.0) GOTO 445
 208       EKIN=TEX/NBL
 209       EKIN2=0.
 210       CALL STEEP(XX)
 211       DO 441 I=1,NBL
 212       IF(NT.EQ.MXGKPV-2) GOTO 441
 213       IF(EKIN2.GT.TEX) GOTO 443
 214       CALL GRNDM(RNDM,1)
 215       RAN1=RNDM(1)
 216       CALL NORMAL(RAN2)
 217       EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
 218       IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1)
 219       EKIN1=EKIN1*XX
 220       EKIN2=EKIN2+EKIN1
 221       IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
 222       IF(EKIN1.LT.0.) EKIN1=0.0001
 223       IPA1=16
 224       PNRAT=1.-ZNO2/ATNO2
 225       CALL GRNDM(RNDM,3)
 226       IF(RNDM(1).GT.PNRAT) IPA1=14
 227       NT=NT+1
 228       SPALL=SPALL+1.
 229       COST=-1.+RNDM(2)*2.
 230       SINT=SQRT(ABS(1.-COST*COST))
 231       PHI=TWPI*RNDM(3)
 232       IPA(NT)=-IPA1
 233       PV(5,NT)=ABS(RMASS(IPA1))
 234       PV(6,NT)=RCHARG(IPA1)
 235       PV(7,NT)=2.
 236       PV(4,NT)=EKIN1+PV(5,NT)
 237       PP=SQRT(ABS(PV(4,NT)**2-PV(5,NT)**2))
 238       PV(1,NT)=PP*SINT*SIN(PHI)
 239       PV(2,NT)=PP*SINT*COS(PHI)
 240       PV(3,NT)=PP*COST
 241   441 CONTINUE
 242   443 IF(ATNO2.LT.10.) GOTO 445
 243       IF(EK.GT.2.0) GOTO 445
 244       II=NT+1
 245       KK=0
 246       EKA=EK
 247       IF(EKA.GT.1.) EKA=EKA*EKA
 248       IF(EKA.LT.0.1) EKA=0.1
 249       IKA=3.6*EXP((ZNO2**2/ATNO2-35.56)/6.45)/EKA
 250       IF(IKA.LE.0) GO TO 445
 251       DO 444 I=1,NT
 252       II=II-1
 253       IF(IPA(II).NE.-14) GOTO 444
 254       IPA(II)=-16
 255       IPA1  = 16
 256       PV(5,II)=ABS(RMASS(IPA1))
 257       PV(6,II)=RCHARG(IPA1)
 258       KK=KK+1
 259       IF(KK.GT.IKA) GOTO 445
 260   444 CONTINUE
 261   445 TEX=ENP(3)
 262       IF(TEX.LT.0.0001) GOTO 40
 263       NBL=IFIX(2.*LOG(ATNO2))
 264       IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
 265       IF(NBL.LE.0) GOTO 40
 266       EKIN=TEX/NBL
 267       EKIN2=0.
 268       CALL STEEP(XX)
 269       IF (NPRT(4)) WRITE(NEWBCD,3004) NBL,TEX
 270       DO 442 I=1,NBL
 271       IF(NT.EQ.MXGKPV-2) GOTO 442
 272       IF(EKIN2.GT.TEX) GOTO 40
 273       CALL GRNDM(RNDM,1)
 274       RAN1=RNDM(1)
 275       CALL NORMAL(RAN2)
 276       EKIN1=-EKIN*LOG(RAN1)-CFA*(1.+0.5*RAN2)
 277       IF(EKIN1.LT.0.0) EKIN1=-0.010*LOG(RAN1)
 278       EKIN1=EKIN1*XX
 279       EKIN2=EKIN2+EKIN1
 280       IF(EKIN2.GT.TEX) EKIN1=TEX-(EKIN2-EKIN1)
 281       IF(EKIN1.LT.0.) EKIN1=0.0001
 282       CALL GRNDM(RNDM,3)
 283       COST=-1.+RNDM(1)*2.
 284       SINT=SQRT(ABS(1.-COST*COST))
 285       PHI=TWPI*RNDM(2)
 286       RAN=RNDM(3)
 287       IPA(NT+1)=-30
 288       IF(RAN.GT.0.60) IPA(NT+1)=-31
 289       IF(RAN.GT.0.90) IPA(NT+1)=-32
 290       INVE=ABS(IPA(NT+1))
 291       PV(5,NT+1)=RMASS(INVE)
 292       SPALL=SPALL+PV(5,NT+1)*1.066
 293       IF(SPALL.GT.ATNO2) GOTO 40
 294       NT=NT+1
 295       PV(6,NT)=RCHARG(INVE)
 296       PV(7,NT)=2.
 297       PV(4,NT)=PV(5,NT)+EKIN1
 298       PP=SQRT(ABS(PV(4,NT)**2-PV(5,NT)**2))
 299       PV(1,NT)=PP*SINT*SIN(PHI)
 300       PV(2,NT)=PP*SINT*COS(PHI)
 301       PV(3,NT)=PP*COST
 302   442 CONTINUE
 303 C**
 304 C** STORE ON EVENT COMMON
 305 C**
 306    40 EKIN=PV(4,MXGKPV)-ABS(PV(5,MXGKPV))
 307       EKIN1=PV(4,MXGKPV-1)-ABS(PV(5,MXGKPV-1))
 308       EKIN2=0.
 309       CALL TDELAY(TOF1)
 310       CALL GRNDM(RNDM,1)
 311       RAN=RNDM(1)
 312       TOF=TOF-TOF1*LOG(RAN)
 313       DO 1 I=1,NT
 314       EKIN2=EKIN2+PV(4,I)-ABS(PV(5,I))
 315       IF(PV(7,I).LT.0.) PV(5,I)=-PV(5,I)
 316       PV(7,I)=TOF
 317       PV(8,I)=ABS(IPA(I))
 318       PV(9,I)=0.
 319     1 PV(10,I)=0.
 320       IF (NPRT(4)) WRITE(NEWBCD,1003) NT,EKIN,EKIN1,EKIN2
 321       INTCT=INTCT+1.
 322       CALL SETCUR(NT)
 323       NTK=NTK+1
 324       IF(NT.EQ.1) GO TO 9999
 325       DO 50 II=2,NT
 326       I=II-1
 327       IF(NTOT.LT.NSIZE/12) GOTO 43
 328       RETURN
 329    43 CALL SETTRK(I)
 330    50 CONTINUE
 331 C
 332  1002 FORMAT(' *TWOB* ',5F10.4,10X,5F10.4/1H ,7X,5F10.4,10X,5F10.4/
 333      $ ' LAB SYSTEM FINAL STATE FOUR VECTORS')
 334  1003 FORMAT(' *TWOB* COMPARISON',2X,I5,2X,3F10.4)
 335  4001 FORMAT(' *TWOB* ',10F10.4,2X,2I3)
 336  4002 FORMAT(' *TWOB* ',7F10.4)
 337  3003 FORMAT(' *TWOB* ',I3,' BLACK TRACK PARTICLES PRODUCED',
 338      $ ' WITH TOTAL KINETIC ENERGY OF ',F8.3,' GEV')
 339  3004 FORMAT(' *TWOB* ',I5,' HEAVY FRAGMENTS PRODUCED',
 340      $ ' WITH TOTAL ENERGY OF',F8.4,' GEV')
 341 C
 342  9999 CONTINUE
 343 C
 344       END