TEvtGen/PHOTOS/phopre.F

   1
   2
   3       SUBROUTINE PHOPRE(IPARR,WT,NEUDAU,NCHARB)
   4 C.----------------------------------------------------------------------
   5 C.
   6 C.    PHOTOS:   Photon radiation in decays
   7 C.
   8 C.    Purpose:  Order (alpha) radiative corrections  are  generated  in
   9 C.              the decay of the IPPAR-th particle in the HEP-like
  10 C.              common /PHOEVT/.  Photon radiation takes place from one
  11 C.              of the charged daughters of the decaying particle IPPAR
  12 C.              WT is calculated, eventual rejection will be performed
  13 C.              later after inclusion of interference weight.
  14 C.
  15 C.    Input Parameter:    IPPAR:  Pointer   to   decaying  particle  in
  16 C.                                /PHOEVT/ and the common itself,
  17 C.
  18 C.    Output Parameters:  Common  /PHOEVT/, either  with  or  without a
  19 C.                                photon(s) added.
  20 C.                        WT      weight of the configuration
  21 C.
  22 C.    Author(s):  Z. Was, B. van Eijk             Created at:  26/11/89
  23 C.                                                Last Update: 26/05/93
  24 C.
  25 C.----------------------------------------------------------------------
  26       IMPLICIT NONE
  27       DOUBLE PRECISION MINMAS,MPASQR,MCHREN
  28       DOUBLE PRECISION BETA,EPS,DEL1,DEL2,DATA
  29       REAL*8 PHOCHA,PHOSPI,PHORAN,PHOCOR,MASSUM
  30       INTEGER IP,IPARR,IPPAR,I,J,ME,NCHARG,NEUPOI,NLAST,THEDUM
  31       INTEGER IDABS,IDUM
  32       INTEGER NCHARB,NEUDAU
  33       REAL*8 WT
  34       INTEGER NMXPHO
  35       PARAMETER (NMXPHO=10000)
  36       INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
  37       REAL*8 PPHO,VPHO
  38       COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
  39      &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
  40       LOGICAL CHKIF
  41       COMMON/PHOIF/CHKIF(NMXPHO)
  42       INTEGER CHAPOI(NMXPHO)
  43       DOUBLE PRECISION MCHSQR,MNESQR
  44       REAL*8 PNEUTR
  45       COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
  46       DOUBLE PRECISION COSTHG,SINTHG
  47       REAL*8 XPHMAX,XPHOTO
  48       COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
  49       REAL*8 ALPHA,XPHCUT
  50       COMMON/PHOCOP/ALPHA,XPHCUT
  51       INTEGER IREP
  52       REAL*8 PROBH,CORWT,XF
  53       COMMON/PHOPRO/PROBH,CORWT,XF,IREP
  54 C--
  55       IPPAR=IPARR
  56 C--   Store pointers for cascade treatement...
  57       IP=IPPAR
  58       NLAST=NPHO
  59       IDUM=1
  60 C--
  61 C--   Check decay multiplicity..
  62       IF (JDAPHO(1,IP).EQ.0) RETURN
  63 C--
  64 C--   Loop over daughters, determine charge multiplicity
  65    10 NCHARG=0
  66       IREP=0
  67       MINMAS=0.D0
  68       MASSUM=0.D0
  69       DO 20 I=JDAPHO(1,IP),JDAPHO(2,IP)
  70 C--
  71 C--
  72 C--   Exclude marked particles, quarks and gluons etc...
  73         IDABS=ABS(IDPHO(I))
  74         IF (CHKIF(I-JDAPHO(1,IP)+3)) THEN
  75           IF (PHOCHA(IDPHO(I)).NE.0) THEN
  76             NCHARG=NCHARG+1
  77             IF (NCHARG.GT.NMXPHO) THEN
  78               DATA=NCHARG
  79               CALL PHOERR(1,'PHOTOS',DATA)
  80             ENDIF
  81             CHAPOI(NCHARG)=I
  82           ENDIF
  83           MINMAS=MINMAS+PPHO(5,I)**2
  84         ENDIF
  85         MASSUM=MASSUM+PPHO(5,I)
  86    20 CONTINUE
  87       IF (NCHARG.NE.0) THEN
  88 C--
  89 C--   Check that sum of daughter masses does not exceed parent mass
  90         IF ((PPHO(5,IP)-MASSUM)/PPHO(5,IP).GT.2.D0*XPHCUT) THEN
  91 C--
  92 C--   Order  charged  particles  according  to decreasing mass, this  to
  93 C--   increase efficiency (smallest mass is treated first).
  94           IF (NCHARG.GT.1) CALL PHOOMA(1,NCHARG,CHAPOI)
  95 C--
  96    30       CONTINUE
  97             DO 70 J=1,3
  98    70       PNEUTR(J)=-PPHO(J,CHAPOI(NCHARG))
  99             PNEUTR(4)=PPHO(5,IP)-PPHO(4,CHAPOI(NCHARG))
 100 C--
 101 C--   Calculate  invariant  mass of 'neutral' etc. systems
 102           MPASQR=PPHO(5,IP)**2
 103           MCHSQR=PPHO(5,CHAPOI(NCHARG))**2
 104           IF ((JDAPHO(2,IP)-JDAPHO(1,IP)).EQ.1) THEN
 105             NEUPOI=JDAPHO(1,IP)
 106             IF (NEUPOI.EQ.CHAPOI(NCHARG)) NEUPOI=JDAPHO(2,IP)
 107             MNESQR=PPHO(5,NEUPOI)**2
 108             PNEUTR(5)=PPHO(5,NEUPOI)
 109           ELSE
 110             MNESQR=PNEUTR(4)**2-PNEUTR(1)**2-PNEUTR(2)**2-PNEUTR(3)**2
 111             MNESQR=MAX(MNESQR,MINMAS-MCHSQR)
 112             PNEUTR(5)=SQRT(MNESQR)
 113           ENDIF
 114 C--
 115 C--   Determine kinematical limit...
 116           XPHMAX=(MPASQR-(PNEUTR(5)+PPHO(5,CHAPOI(NCHARG)))**2)/MPASQR
 117 C--
 118 C--   Photon energy fraction...
 119           CALL PHOENE(MPASQR,MCHREN,BETA,IDPHO(CHAPOI(NCHARG)))
 120 C--
 121 C--   Energy fraction not too large (very seldom) ? Define angle.
 122           IF ((XPHOTO.LT.XPHCUT).OR.(XPHOTO.GT.XPHMAX)) THEN
 123 C--
 124 C--   No radiation was accepted, check  for more daughters  that may ra-
 125 C--   diate and correct radiation probability...
 126             NCHARG=NCHARG-1
 127             IF (NCHARG.GT.0) THEN
 128               IREP=IREP+1
 129               GOTO 30
 130             ENDIF
 131           ELSE
 132 C--
 133 C--   Angle is generated  in  the  frame defined  by  charged vector and
 134 C--   PNEUTR, distribution is taken in the infrared limit...
 135             EPS=MCHREN/(1.D0+BETA)
 136 C--
 137 C--   Calculate sin(theta) and cos(theta) from interval variables
 138             DEL1=(2.D0-EPS)*(EPS/(2.D0-EPS))**PHORAN(THEDUM)
 139             DEL2=2.D0-DEL1
 140             COSTHG=(1.D0-DEL1)/BETA
 141             SINTHG=SQRT(DEL1*DEL2-MCHREN)/BETA
 142 C--
 143 C--   Determine spin of  particle and construct code  for matrix element
 144             ME=2.D0*PHOSPI(IDPHO(CHAPOI(NCHARG)))+1.D0
 145 C--
 146 C--   Weighting procedure with 'exact' matrix element, reconstruct kine-
 147 C--   matics for photon, neutral and charged system and update /PHOEVT/.
 148 C--   Find pointer to the first component of 'neutral' system
 149       DO  I=JDAPHO(1,IP),JDAPHO(2,IP)
 150         IF (I.NE.CHAPOI(NCHARG)) THEN
 151           NEUDAU=I
 152           GOTO 51
 153         ENDIF
 154       ENDDO
 155 C--
 156 C--   Pointer not found...
 157       DATA=NCHARG
 158       CALL PHOERR(5,'PHOKIN',DATA)
 159  51   CONTINUE
 160       NCHARB=CHAPOI(NCHARG)
 161       NCHARB=NCHARB-JDAPHO(1,IP)+3
 162       NEUDAU=NEUDAU-JDAPHO(1,IP)+3
 163         WT=PHOCOR(MPASQR,MCHREN,ME)
 164
 165           ENDIF
 166         ELSE
 167           DATA=PPHO(5,IP)-MASSUM
 168           CALL PHOERR(10,'PHOTOS',DATA)
 169         ENDIF
 170       ENDIF
 171 C--
 172       RETURN
 173       END