GEANT321/gheisha/pimabs.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 PIMABS(NOPT)
  13 C
  14 C *** CHARGED PION ABSORPTION BY A NUCLEUS ***
  15 C *** NVE 04-MAR-1988 CERN GENEVA ***
  16 C
  17 C ORIGIN : H.FESEFELDT (09-JULY-1987)
  18 C
  19 C PANOFSKY RATIO (PI- P --> N PI0/PI- P --> N GAMMA) = 3/2
  20 C FOR CAPTURE ON PROTON (HYDROGEN),
  21 C STAR PRODUCTION FOR HEAVIER ELEMENTS
  22 C
  23 #include "geant321/s_defcom.inc"
  24       DIMENSION RNDM(4)
  25       SAVE NT
  26 C
  27       CALL COMPO
  28       PV(1,1)=0.
  29       PV(2,1)=0.
  30       PV(3,1)=0.
  31       PV(4,1)=RMASS(9)
  32       PV(5,1)=RMASS(9)
  33       PV(6,1)=-1.
  34       PV(7,1)=TOF
  35       PV(8,1)=IPART
  36       PV(9,1)=0.
  37       PV(10,1)=USERW
  38       IER(87)=IER(87)+1
  39       IF(ATNO2.GT.1.5) GOTO 30
  40       CALL GRNDM(RNDM,2)
  41       RAN=RNDM(1)
  42       ISW=1
  43       IF(RAN.LT.0.33) ISW=2
  44       NOPT=ISW
  45       RAN=RNDM(2)
  46       TOF1=-25.*LOG(RAN)
  47       TOF1=20.*TOF1
  48       IF(ISW.EQ.1) GOTO 20
  49       PCM=0.02
  50       CALL GRNDM(RNDM,2)
  51       PHI=2.*PI*RNDM(1)
  52       COST=-1.+2.*RNDM(2)
  53       SINT=SQRT(ABS(1.-COST*COST))
  54       PV(1,2)=PCM*SINT*COS(PHI)
  55       PV(2,2)=PCM*SINT*SIN(PHI)
  56       PV(3,2)=PCM*COST
  57       PV(4,2)=PCM
  58       PV(5,2)=0.
  59       PV(7,2)=TOF+TOF1
  60       PV(9,2)=0.
  61       PV(10,2)=0.
  62       PV(6,2)=0.
  63       PV(8,2)=1.
  64       GOTO 21
  65    20 PV(1,2)=PV(1,1)
  66       PV(2,2)=PV(2,1)
  67       PV(3,2)=PV(3,1)
  68       PV(4,2)=PV(4,1)
  69       PV(5,2)=PV(5,1)
  70       PV(6,2)=0.
  71       PV(7,2)=TOF+TOF1
  72       PV(8,2)=8.
  73       PV(9,2)=0.
  74       PV(10,2)=0.
  75    21 INTCT=INTCT+1.
  76       CALL SETCUR(2)
  77       NTK=NTK+1
  78       IF(NPRT(3))
  79      *WRITE(NEWBCD,1002) XEND,YEND,ZEND,P,NCH
  80 1002  FORMAT(1H0,'PION ABSORBTION   POSITION',3(2X,F8.2),2X,
  81      * 'PI0/GAMMA MOMENTUM,CHARGE',2X,F8.4,2X,F4.1)
  82       GO TO 9999
  83 C**
  84 C** STAR PRODUCTION FOR PION ABSORPTION IN HEAVY ELEMENTS
  85 C**
  86    30 ENP(1)=0.0135
  87       ENP(3)=0.0058
  88       NT=1
  89       TEX=ENP(1)
  90       BLACK=0.5*LOG(ATNO2)
  91       CALL POISSO(BLACK,NBL)
  92       IF(NBL.LE.0) NBL=1
  93       IF(NPRT(3))
  94      *WRITE(NEWBCD,3003) NBL,TEX
  95       IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
  96       EKIN=TEX/NBL
  97       EKIN2=0.
  98       DO 31 I=1,NBL
  99       IF(NT.EQ.MXGKPV-2) GOTO 31
 100       CALL GRNDM(RNDM,4)
 101       RAN2=RNDM(1)
 102       EKIN1=-EKIN*LOG(RAN2)
 103       EKIN2=EKIN2+EKIN1
 104       IPA1=16
 105       PNRAT=1.-ZNO2/ATNO2
 106       IF(RNDM(2).GT.PNRAT) IPA1=14
 107       NT=NT+1
 108       COST=-1.+RNDM(3)*2.
 109       SINT=SQRT(ABS(1.-COST*COST))
 110       PHI=TWPI*RNDM(4)
 111       IPA(NT)=-IPA1
 112       PV(5,NT)=ABS(RMASS(IPA1))
 113       PV(6,NT)=RCHARG(IPA1)
 114       PV(7,NT)=2.
 115       PV(4,NT)=EKIN1+PV(5,NT)
 116       PP=SQRT(ABS(PV(4,NT)**2-PV(5,NT)**2))
 117       PV(1,NT)=PP*SINT*SIN(PHI)
 118       PV(2,NT)=PP*SINT*COS(PHI)
 119       PV(3,NT)=PP*COST
 120       IF(EKIN2.GT.TEX) GOTO 33
 121    31 CONTINUE
 122    33 TEX=ENP(3)
 123       BLACK=0.50*LOG(ATNO2)
 124       CALL POISSO(BLACK,NBL)
 125       IF(NT+NBL.GT.MXGKPV-2) NBL=MXGKPV-2-NT
 126       IF(NBL.LE.0) NBL=1
 127       EKIN=TEX/NBL
 128       EKIN2=0.
 129       IF(NPRT(3))
 130      *WRITE(NEWBCD,3004) NBL,TEX
 131       DO 32 I=1,NBL
 132       IF(NT.EQ.MXGKPV-2) GOTO 32
 133       CALL GRNDM(RNDM,4)
 134       RAN2=RNDM(1)
 135       EKIN1=-EKIN*LOG(RAN2)
 136       EKIN2=EKIN2+EKIN1
 137       NT=NT+1
 138       COST=-1.+RNDM(2)*2.
 139       SINT=SQRT(ABS(1.-COST*COST))
 140       PHI=TWPI*RNDM(3)
 141       RAN=RNDM(4)
 142       IPA(NT)=-30
 143       IF(RAN.GT.0.60) IPA(NT)=-31
 144       IF(RAN.GT.0.90) IPA(NT)=-32
 145       INVE=ABS(IPA(NT))
 146       PV(5,NT)=RMASS(INVE)
 147       PV(6,NT)=RCHARG(INVE)
 148       PV(7,NT)=2.
 149       PV(4,NT)=PV(5,NT)+EKIN1
 150       PP=SQRT(ABS(PV(4,NT)**2-PV(5,NT)**2))
 151       PV(1,NT)=PP*SINT*SIN(PHI)
 152       PV(2,NT)=PP*SINT*COS(PHI)
 153       PV(3,NT)=PP*COST
 154       IF(EKIN2.GT.TEX) GOTO 40
 155    32 CONTINUE
 156 C**
 157 C** STORE ON EVENT COMMON
 158 C**
 159    40 CALL GRNDM(RNDM,1)
 160       RAN=RNDM(1)
 161       TOF1=-25.*LOG(RAN)
 162       TOF1=20.*TOF1
 163       DO 41 I=2,NT
 164       IF(PV(7,I).LT.0.) PV(5,I)=-PV(5,I)
 165       PV(7,I)=TOF+TOF1
 166       PV(8,I)=ABS(IPA(I))
 167       PV(9,I)=0.
 168    41 PV(10,I)=0.
 169       INTCT=INTCT+1.
 170       CALL SETCUR(2)
 171       NTK=NTK+1
 172       IF(NT.EQ.2) GO TO 9999
 173       DO 50 I=3,NT
 174       IF(NTOT.LT.NSIZE/12) GOTO 43
 175       GO TO 9999
 176    43 CALL SETTRK(I)
 177    50 CONTINUE
 178 C
 179  3003 FORMAT(1H ,I3,' BLACK TRACK PARTICLES PRODUCED WITH TOTAL KINETIC
 180      * ENERGY OF ',F8.3,' GEV')
 181  3004 FORMAT(1H ,I5,' HEAVY FRAGMENTS WITH TOTAL ENERGY OF',F8.4,' GEV')
 182 C
 183  9999 CONTINUE
 184       END