TEvtGen/PHOTOS/phodo.F

   1       SUBROUTINE PHODO(IP,NCHARB,NEUDAU)
   2 C.----------------------------------------------------------------------
   3 C.
   4 C.    PHOTOS:   PHOton radiation in  decays DOing of KINematics
   5 C.
   6 C.    Purpose:  Starting  from   the  charged  particle energy/momentum,
   7 C.              PNEUTR, photon  energy  fraction and photon  angle  with
   8 C.              respect  to  the axis formed by charged particle energy/
   9 C.              momentum  vector  and PNEUTR, scale the energy/momentum,
  10 C.              keeping the original direction of the neutral system  in
  11 C.              the lab. frame untouched.
  12 C.
  13 C.    Input Parameters:   IP:      Pointer  to   decaying  particle   in
  14 C.                                 /PHOEVT/  and   the   common   itself
  15 C.                        NCHARB:  pointer to the charged radiating
  16 C.                                 daughter in /PHOEVT/.
  17 C.                        NEUDAU:  pointer to the first neutral daughter
  18 C.    Output Parameters:  Common /PHOEVT/, with photon added.
  19 C.
  20 C.    Author(s):  Z. Was, B. van Eijk             Created at:  26/11/89
  21 C.                                                Last Update: 27/05/93
  22 C.
  23 C.----------------------------------------------------------------------
  24       IMPLICIT NONE
  25       DOUBLE PRECISION PHOAN1,PHOAN2,ANGLE,FI1,FI3,FI4,FI5,TH1,TH3,TH4
  26       DOUBLE PRECISION PARNE,QNEW,QOLD,DATA
  27       INTEGER IP,FI3DUM,I,J,NEUDAU,FIRST,LAST
  28       INTEGER NCHARB
  29       REAL*8 EPHOTO,PMAVIR,PHOTRI
  30       REAL*8 GNEUT,PHORAN,CCOSTH,SSINTH,PVEC(4)
  31       INTEGER NMXPHO
  32       PARAMETER (NMXPHO=10000)
  33       INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
  34       REAL*8 PPHO,VPHO
  35       COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),
  36      &JMOPHO(2,NMXPHO),JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
  37       DOUBLE PRECISION MCHSQR,MNESQR
  38       REAL*8 PNEUTR
  39       COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
  40       DOUBLE PRECISION COSTHG,SINTHG
  41       REAL*8 XPHMAX,XPHOTO
  42       COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
  43       REAL*8 PI,TWOPI
  44       COMMON/PHPICO/PI,TWOPI
  45 C--
  46       EPHOTO=XPHOTO*PPHO(5,IP)/2.D0
  47       PMAVIR=SQRT(PPHO(5,IP)*(PPHO(5,IP)-2.D0*EPHOTO))
  48 C--
  49 C--   Reconstruct  kinematics  of  charged particle  and  neutral system
  50       FI1=PHOAN1(PNEUTR(1),PNEUTR(2))
  51 C--
  52 C--   Choose axis along  z of  PNEUTR, calculate  angle  between x and y
  53 C--   components  and z  and x-y plane and  perform Lorentz transform...
  54       TH1=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)**2+PNEUTR(2)**2))
  55       CALL PHORO3(-FI1,PNEUTR(1))
  56       CALL PHORO2(-TH1,PNEUTR(1))
  57 C--
  58 C--   Take  away  photon energy from charged particle and PNEUTR !  Thus
  59 C--   the onshell charged particle  decays into virtual charged particle
  60 C--   and photon.  The virtual charged  particle mass becomes:
  61 C--   SQRT(PPHO(5,IP)*(PPHO(5,IP)-2*EPHOTO)).  Construct  new PNEUTR mo-
  62 C--   mentum in the rest frame of the parent:
  63 C--   1) Scaling parameters...
  64       QNEW=PHOTRI(PMAVIR,PNEUTR(5),PPHO(5,NCHARB))
  65       QOLD=PNEUTR(3)
  66       GNEUT=(QNEW**2+QOLD**2+MNESQR)/(QNEW*QOLD+SQRT((QNEW**2+MNESQR)*
  67      &(QOLD**2+MNESQR)))
  68       IF (GNEUT.LT.1.D0) THEN
  69         DATA=0.D0
  70         CALL PHOERR(4,'PHOKIN',DATA)
  71       ENDIF
  72       PARNE=GNEUT-SQRT(MAX(GNEUT**2-1.0D0,0.D0))
  73 C--
  74 C--   2) ...reductive boost...
  75       CALL PHOBO3(PARNE,PNEUTR)
  76 C--
  77 C--   ...calculate photon energy in the reduced system...
  78       NPHO=NPHO+1
  79       ISTPHO(NPHO)=1
  80       IDPHO(NPHO) =22
  81 C--   Photon mother and daughter pointers !
  82       JMOPHO(1,NPHO)=IP
  83       JMOPHO(2,NPHO)=0
  84       JDAPHO(1,NPHO)=0
  85       JDAPHO(2,NPHO)=0
  86       PPHO(4,NPHO)=EPHOTO*PPHO(5,IP)/PMAVIR
  87 C--
  88 C--   ...and photon momenta
  89       CCOSTH=-COSTHG
  90       SSINTH=SINTHG
  91       TH3=PHOAN2(CCOSTH,SSINTH)
  92       FI3=TWOPI*PHORAN(FI3DUM)
  93       PPHO(1,NPHO)=PPHO(4,NPHO)*SINTHG*COS(FI3)
  94       PPHO(2,NPHO)=PPHO(4,NPHO)*SINTHG*SIN(FI3)
  95 C--
  96 C--   Minus sign because axis opposite direction of charged particle !
  97       PPHO(3,NPHO)=-PPHO(4,NPHO)*COSTHG
  98       PPHO(5,NPHO)=0.D0
  99 C--
 100 C--   Rotate in order to get photon along z-axis
 101       CALL PHORO3(-FI3,PNEUTR(1))
 102       CALL PHORO3(-FI3,PPHO(1,NPHO))
 103       CALL PHORO2(-TH3,PNEUTR(1))
 104       CALL PHORO2(-TH3,PPHO(1,NPHO))
 105       ANGLE=EPHOTO/PPHO(4,NPHO)
 106 C--
 107 C--   Boost to the rest frame of decaying particle
 108       CALL PHOBO3(ANGLE,PNEUTR(1))
 109       CALL PHOBO3(ANGLE,PPHO(1,NPHO))
 110 C--
 111 C--   Back in the parent rest frame but PNEUTR not yet oriented !
 112       FI4=PHOAN1(PNEUTR(1),PNEUTR(2))
 113       TH4=PHOAN2(PNEUTR(3),SQRT(PNEUTR(1)**2+PNEUTR(2)**2))
 114       CALL PHORO3(FI4,PNEUTR(1))
 115       CALL PHORO3(FI4,PPHO(1,NPHO))
 116 C--
 117         DO 60 I=2,4
 118    60   PVEC(I)=0.D0
 119         PVEC(1)=1.D0
 120         CALL PHORO3(-FI3,PVEC)
 121         CALL PHORO2(-TH3,PVEC)
 122         CALL PHOBO3(ANGLE,PVEC)
 123         CALL PHORO3(FI4,PVEC)
 124         CALL PHORO2(-TH4,PNEUTR)
 125         CALL PHORO2(-TH4,PPHO(1,NPHO))
 126         CALL PHORO2(-TH4,PVEC)
 127         FI5=PHOAN1(PVEC(1),PVEC(2))
 128 C--
 129 C--   Charged particle restores original direction
 130         CALL PHORO3(-FI5,PNEUTR)
 131         CALL PHORO3(-FI5,PPHO(1,NPHO))
 132         CALL PHORO2(TH1,PNEUTR(1))
 133         CALL PHORO2(TH1,PPHO(1,NPHO))
 134         CALL PHORO3(FI1,PNEUTR)
 135         CALL PHORO3(FI1,PPHO(1,NPHO))
 136 C--   See whether neutral system has multiplicity larger than 1...
 137       IF ((JDAPHO(2,IP)-JDAPHO(1,IP)).GT.1) THEN
 138 C--   Find pointers to components of 'neutral' system
 139 C--
 140         FIRST=NEUDAU
 141         LAST=JDAPHO(2,IP)
 142         DO 70 I=FIRST,LAST
 143           IF (I.NE.NCHARB.AND.(JMOPHO(1,I).EQ.IP)) THEN
 144 C--
 145 C--   Reconstruct kinematics...
 146             CALL PHORO3(-FI1,PPHO(1,I))
 147             CALL PHORO2(-TH1,PPHO(1,I))
 148 C--
 149 C--   ...reductive boost
 150             CALL PHOBO3(PARNE,PPHO(1,I))
 151 C--
 152 C--   Rotate in order to get photon along z-axis
 153             CALL PHORO3(-FI3,PPHO(1,I))
 154             CALL PHORO2(-TH3,PPHO(1,I))
 155 C--
 156 C--   Boost to the rest frame of decaying particle
 157             CALL PHOBO3(ANGLE,PPHO(1,I))
 158 C--
 159 C--   Back in the parent rest-frame but PNEUTR not yet oriented.
 160             CALL PHORO3(FI4,PPHO(1,I))
 161             CALL PHORO2(-TH4,PPHO(1,I))
 162 C--
 163 C--   Charged particle restores original direction
 164             CALL PHORO3(-FI5,PPHO(1,I))
 165             CALL PHORO2(TH1,PPHO(1,I))
 166             CALL PHORO3(FI1,PPHO(1,I))
 167           ENDIF
 168    70   CONTINUE
 169       ELSE
 170 C--
 171 C--   ...only one 'neutral' particle in addition to photon!
 172         DO 80 J=1,4
 173    80   PPHO(J,NEUDAU)=PNEUTR(J)
 174       ENDIF
 175 C--
 176 C--   All 'neutrals' treated, fill /PHOEVT/ for charged particle...
 177       DO 90 J=1,3
 178    90 PPHO(J,NCHARB)=-(PPHO(J,NPHO)+PNEUTR(J))
 179       PPHO(4,NCHARB)=PPHO(5,IP)-(PPHO(4,NPHO)+PNEUTR(4))
 180 C--
 181       END