forW-ME.f

00001 
00002 c ###### 
00003 c  standard routine !!!
00004 c ######
00005 c$$$      FUNCTION PHOCOR(MPASQR,MCHREN,ME)
00006 c$$$C.----------------------------------------------------------------------
00007 c$$$C.
00008 c$$$C.    PHOTOS:   PHOton radiation in decays CORrection weight from
00009 c$$$C.              matrix elements
00010 c$$$C.
00011 c$$$C.    Purpose:  Calculate  photon  angle.  The reshaping functions  will
00012 c$$$C.              have  to  depend  on the spin S of the charged particle.
00013 c$$$C.              We define:  ME = 2 * S + 1 !
00014 c$$$C.
00015 c$$$C.    Input Parameters:  MPASQR:  Parent mass squared,
00016 c$$$C.                       MCHREN:  Renormalised mass of charged system,
00017 c$$$C.                       ME:      2 * spin + 1 determines matrix element
00018 c$$$C.
00019 c$$$C.    Output Parameter:  Function value.
00020 c$$$C.
00021 c$$$C.    Author(s):  Z. Was, B. van Eijk             Created at:  26/11/89
00022 c$$$C.                                                Last Update: 21/03/93
00023 c$$$C.
00024 c$$$C.----------------------------------------------------------------------
00025 c$$$      IMPLICIT NONE
00026 c$$$      DOUBLE PRECISION MPASQR,MCHREN,BETA,XX,YY,DATA
00027 c$$$      INTEGER ME
00028 c$$$      REAL*8 PHOCOR,PHOFAC,WT1,WT2,WT3
00029 c$$$      DOUBLE PRECISION MCHSQR,MNESQR
00030 c$$$      REAL*8 PNEUTR
00031 c$$$      COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
00032 c$$$      DOUBLE PRECISION COSTHG,SINTHG
00033 c$$$      REAL*8 XPHMAX,XPHOTO
00034 c$$$      COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
00035 c$$$      INTEGER IREP
00036 c$$$      REAL*8 PROBH,CORWT,XF
00037 c$$$      COMMON/PHOPRO/PROBH,CORWT,XF,IREP
00038 c$$$C--
00039 c$$$C--   Shaping (modified by ZW)...
00040 c$$$      XX=4.D0*MCHSQR/MPASQR*(1.D0-XPHOTO)/(1.D0-XPHOTO+(MCHSQR-MNESQR)/
00041 c$$$     &MPASQR)**2
00042 c$$$      IF (ME.EQ.1) THEN
00043 c$$$        YY=1.D0
00044 c$$$        WT3=(1.D0-XPHOTO/XPHMAX)/((1.D0+(1.D0-XPHOTO/XPHMAX)**2)/2.D0)
00045 c$$$      ELSEIF (ME.EQ.2) THEN
00046 c$$$        YY=0.5D0*(1.D0-XPHOTO/XPHMAX+1.D0/(1.D0-XPHOTO/XPHMAX))
00047 c$$$        WT3=1.D0
00048 c$$$      ELSEIF ((ME.EQ.3).OR.(ME.EQ.4).OR.(ME.EQ.5)) THEN
00049 c$$$        YY=1.D0
00050 c$$$        WT3=(1.D0+(1.D0-XPHOTO/XPHMAX)**2-(XPHOTO/XPHMAX)**3)/
00051 c$$$     &  (1.D0+(1.D0-XPHOTO/XPHMAX)** 2)
00052 c$$$      ELSE
00053 c$$$        DATA=(ME-1.D0)/2.D0
00054 c$$$        CALL PHOERR(6,'PHOCOR',DATA)
00055 c$$$        YY=1.D0
00056 c$$$        WT3=1.D0
00057 c$$$      ENDIF
00058 c$$$      BETA=SQRT(1.D0-XX)
00059 c$$$      WT1=(1.D0-COSTHG*SQRT(1.D0-MCHREN))/(1.D0-COSTHG*BETA)
00060 c$$$      WT2=(1.D0-XX/YY/(1.D0-BETA**2*COSTHG**2))*(1.D0+COSTHG*BETA)/2.D0
00061 c$$$      WT2=WT2*PHOFAC(1)
00062 c$$$      PHOCOR=WT1*WT2*WT3
00063 c$$$      CORWT=PHOCOR
00064 c$$$      IF (PHOCOR.GT.1.D0) THEN
00065 c$$$        DATA=PHOCOR
00066 c$$$        CALL PHOERR(3,'PHOCOR',DATA)
00067 c$$$      ENDIF
00068 c$$$      RETURN
00069 c$$$      END
00070 c$$$
00071 
00072 
00073       FUNCTION PHOCORN(MPASQR,MCHREN,ME)
00074 c ###### 
00075 c  replace with, 
00076 c ######
00077 
00078 C.----------------------------------------------------------------------
00079 C.
00080 C.    PHOTOS:   PHOton radiation in decays CORrection weight from
00081 C.              matrix elements This version for spin 1/2 is verified for
00082 C.              W decay only
00083 C.    Purpose:  Calculate  photon  angle.  The reshaping functions  will
00084 C.              have  to  depend  on the spin S of the charged particle.
00085 C.              We define:  ME = 2 * S + 1 !
00086 C.
00087 C.    Input Parameters:  MPASQR:  Parent mass squared,
00088 C.                       MCHREN:  Renormalised mass of charged system,
00089 C.                       ME:      2 * spin + 1 determines matrix element
00090 C.
00091 C.    Output Parameter:  Function value.
00092 C.
00093 C.    Author(s):  Z. Was, B. van Eijk, G. Nanava  Created at:  26/11/89
00094 C.                                                Last Update: 01/11/12
00095 C.
00096 C.----------------------------------------------------------------------
00097       IMPLICIT NONE
00098       DOUBLE PRECISION MPASQR,MCHREN,BETA,BETA0,BETA1,XX,YY,DATA
00099       INTEGER ME
00100       REAL*8 PHOCOR,PHOFAC,WT1,WT2,WT3,PHOTRI,S1,PHOCORN
00101       DOUBLE PRECISION MCHSQR,MNESQR
00102       REAL*8 PNEUTR
00103       COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
00104       DOUBLE PRECISION COSTHG,SINTHG,phocorWT3,phocorWT2,phocorWT1
00105       REAL*8 XPHMAX,XPHOTO
00106       COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
00107       common/phocorWT/phocorWT3,phocorWT2,phocorWT1
00108       INTEGER IREP
00109       REAL*8 PROBH,CORWT,XF
00110       COMMON/PHOPRO/PROBH,CORWT,XF,IREP
00111 
00112 C--
00113 C--   Shaping (modified by ZW)...
00114       XX=4.D0*MCHSQR/MPASQR*(1.D0-XPHOTO)/(1.D0-XPHOTO+(MCHSQR-MNESQR)/
00115      &MPASQR)**2
00116       IF (ME.EQ.1) THEN
00117         S1=MPASQR  * (1.D0-XPHOTO)
00118         BETA0=2*PHOTRI(1D0,dsqrt(MCHSQR/MPASQR),dsqrt(MNESQR/MPASQR))
00119         BETA1=2*PHOTRI(1D0,dsqrt(MCHSQR/S1),dsqrt(MNESQR/S1))
00120         WT1= (1.D0-COSTHG*SQRT(1.D0-MCHREN))
00121      $      /((1D0+(1D0-XPHOTO/XPHMAX)**2)/2.D0)*XPHOTO          ! de-presampler
00122      $     
00123         WT2= beta1/beta0*XPHOTO                                  !phase space jacobians
00124         WT3=  beta1**2* (1D0-COSTHG**2) * (1D0-XPHOTO)/XPHOTO**2 ! matrix element
00125      $    /(1D0 +MCHSQR/S1-MNESQR/S1-BETA1*COSTHG)**2/2D0 
00126       ELSEIF (ME.EQ.2) THEN
00127         S1=MPASQR  * (1.D0-XPHOTO)
00128         BETA0=2*PHOTRI(1D0,dsqrt(MCHSQR/MPASQR),dsqrt(MNESQR/MPASQR))
00129         BETA1=2*PHOTRI(1D0,dsqrt(MCHSQR/S1),dsqrt(MNESQR/S1))
00130         WT1= (1.D0-COSTHG*SQRT(1.D0-MCHREN))
00131      $      /((1D0+(1D0-XPHOTO/XPHMAX)**2)/2.D0)*XPHOTO          ! de-presampler
00132          
00133         WT2= beta1/beta0*XPHOTO                                  ! phase space jacobians
00134 
00135         WT3= beta1**2* (1D0-COSTHG**2) * (1D0-XPHOTO)/XPHOTO**2  ! matrix element
00136      $       /(1D0 +MCHSQR/S1-MNESQR/S1-BETA1*COSTHG)**2/2D0 
00137         WT3=WT3*(1-xphoto/xphmax+0.5*(xphoto/xphmax)**2)/(1-xphoto/xphmax)
00138 c       print*,"WT3=",wt3
00139         phocorWT3=WT3
00140         phocorWT2=WT2
00141         phocorWT1=WT1
00142 
00143 c       YY=0.5D0*(1.D0-XPHOTO/XPHMAX+1.D0/(1.D0-XPHOTO/XPHMAX))
00144 c       BETA=SQRT(1.D0-XX)
00145 c       WT1=(1.D0-COSTHG*SQRT(1.D0-MCHREN))/(1.D0-COSTHG*BETA)
00146 c       WT2=(1.D0-XX/YY/(1.D0-BETA**2*COSTHG**2))*(1.D0+COSTHG*BETA)/2.D0
00147 c       WT3=1.D0
00148       ELSEIF ((ME.EQ.3).OR.(ME.EQ.4).OR.(ME.EQ.5)) THEN
00149         YY=1.D0
00150         BETA=SQRT(1.D0-XX)
00151         WT1=(1.D0-COSTHG*SQRT(1.D0-MCHREN))/(1.D0-COSTHG*BETA)
00152         WT2=(1.D0-XX/YY/(1.D0-BETA**2*COSTHG**2))*(1.D0+COSTHG*BETA)/2.D0
00153         WT3=(1.D0+(1.D0-XPHOTO/XPHMAX)**2-(XPHOTO/XPHMAX)**3)/
00154      &  (1.D0+(1.D0-XPHOTO/XPHMAX)** 2)
00155       ELSE
00156         DATA=(ME-1.D0)/2.D0
00157         CALL PHOERR(6,'PHOCOR',DATA)
00158         YY=1.D0
00159         BETA=SQRT(1.D0-XX)
00160         WT1=(1.D0-COSTHG*SQRT(1.D0-MCHREN))/(1.D0-COSTHG*BETA)
00161         WT2=(1.D0-XX/YY/(1.D0-BETA**2*COSTHG**2))*(1.D0+COSTHG*BETA)/2.D0
00162         WT3=1.D0
00163       ENDIF
00164       WT2=WT2*PHOFAC(1)
00165       PHOCOR=WT1*WT2*WT3
00166       PHOCORN=PHOCOR
00167       CORWT=PHOCOR
00168       IF (PHOCOR.GT.1.D0) THEN
00169         DATA=PHOCOR
00170         CALL PHOERR(3,'PHOCOR',DATA)
00171       ENDIF
00172       RETURN
00173       END
00174       SUBROUTINE PHOBWnlo(WT)
00175 C.----------------------------------------------------------------------
00176 C.
00177 C.    PHOTOS:   PHOtos Boson W correction weight
00178 C.
00179 C.    Purpose:  calculates correction weight due to amplitudes of 
00180 C.              emission from W boson. It is ecact, but not verified
00181 C.              for exponentiation yet.
00182 C.              
00183 C.              
00184 C.              
00185 C.
00186 C.    Input Parameters:  Common /PHOEVT/, with photon added.
00187 C.                       wt  to be corrected
00188 C.                       
00189 C.                       
00190 C.                       
00191 C.    Output Parameters: wt
00192 C.
00193 C.    Author(s):  G. Nanava, Z. Was               Created at:  13/03/03
00194 C.                                                Last Update: 01/11/12
00195 C.
00196 C.----------------------------------------------------------------------
00197       IMPLICIT NONE
00198       DOUBLE PRECISION WT
00199       INTEGER NMXPHO
00200       PARAMETER (NMXPHO=10000)
00201       INTEGER IDPHO,ISTPHO,JDAPHO,JMOPHO,NEVPHO,NPHO
00202       REAL*8 PPHO,VPHO
00203       INTEGER PHLUN
00204       COMMON/PHOLUN/PHLUN
00205       REAL*8 ALPHA,XPHCUT
00206       COMMON/PHOCOP/ALPHA,XPHCUT
00207       COMMON/PHOEVT/NEVPHO,NPHO,ISTPHO(NMXPHO),IDPHO(NMXPHO),JMOPHO(2,NMXPHO),
00208      &              JDAPHO(2,NMXPHO),PPHO(5,NMXPHO),VPHO(4,NMXPHO)
00209       INTEGER I,JJ,II,I3,I4,IJ
00210       DOUBLE PRECISION EMU,MCHREN,BETA,COSTHG,MPASQR,XPH,
00211      &                 PW(0:3),PMU(0:3),PPHOT(0:3),PNE(0:3),
00212      &                 B_PW(0:3),B_PNE(0:3),B_PMU(0:3),AMPSQR,SANC_MC_INIT
00213       DOUBLE PRECISION WDecayAmplitudeSqrKS_1ph,WDecayEikonalSqrKS_1ph,
00214      &                 WDecayBornAmpSqrKS_1ph,EIKONALFACTOR
00215       EXTERNAL WDecayAmplitudeSqrKS_1ph,WDecayEikonalSqrKS_1ph,WDecayBornAmpSqrKS_1ph
00216 
00217       INTEGER  d_h_nmxhep         ! maximum number of particles
00218       PARAMETER (d_h_NMXHEP=10000)
00219       REAL*8  d_h_phep,  d_h_vhep ! to be real*4/ *8  depending on host
00220       INTEGER d_h_nevhep,d_h_nhep,d_h_isthep,d_h_idhep,d_h_jmohep,
00221      $        d_h_jdahep
00222       COMMON /ph_hepevt/
00223      $      d_h_nevhep,               ! serial number
00224      $      d_h_nhep,                 ! number of particles
00225      $      d_h_isthep(d_h_nmxhep),   ! status code
00226      $      d_h_idhep(d_h_nmxhep),    ! particle ident KF
00227      $      d_h_jmohep(2,d_h_nmxhep), ! parent particles
00228      $      d_h_jdahep(2,d_h_nmxhep), ! childreen particles
00229      $      d_h_phep(5,d_h_nmxhep),   ! four-momentum, mass [GeV]
00230      $      d_h_vhep(4,d_h_nmxhep)    ! vertex [mm]
00231 
00232       SAVE SANC_MC_INIT
00233       DATA SANC_MC_INIT /-123456789D0/
00234       INTEGER           mcLUN
00235       DOUBLE PRECISION  spV(0:3),bet(0:3)
00236       DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
00237       COMMON /Kleiss_Stirling/spV,bet
00238       COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN        
00239 !      write(*,*) 'IDPHOs=',IDPHO(1),IDPHO(2),IDPHO(3),IDPHO(4),IDPHO(5)
00240 !      write(*,*) 'IDPHOs=',JDAPHO(1,1),npho
00241 !      write(*,*) 'd_h_IDPHOs=',d_h_IDhep(1),d_h_IDhep(2),d_h_IDhep(3),d_h_IDhep(4),d_h_IDhep(5)
00242 
00243 C--
00244         IF (ABS(IDPHO(1)).EQ.24.AND.
00245      $     ABS(IDPHO(JDAPHO(1,1)  )).GE.11.AND.
00246      $     ABS(IDPHO(JDAPHO(1,1)  )).LE.16.AND.
00247      $     ABS(IDPHO(JDAPHO(1,1)+1)).GE.11.AND.
00248      $     ABS(IDPHO(JDAPHO(1,1)+1)).LE.16     ) THEN
00249 
00250            IF(
00251      $      ABS(IDPHO(JDAPHO(1,1)  )).EQ.11.OR.
00252      $      ABS(IDPHO(JDAPHO(1,1)  )).EQ.13.OR.
00253      $      ABS(IDPHO(JDAPHO(1,1)  )).EQ.15    ) THEN 
00254               I=JDAPHO(1,1)
00255            ELSE
00256               I=JDAPHO(1,1)+1
00257            ENDIF
00258 c..        muon energy   
00259            EMU=PPHO(4,I)
00260 c..        muon mass square
00261            MCHREN=ABS(PPHO(4,I)**2-PPHO(3,I)**2
00262      $               -PPHO(2,I)**2-PPHO(1,I)**2)
00263            BETA=SQRT(1- MCHREN/ PPHO(4,I)**2)
00264            COSTHG=(PPHO(3,I)*PPHO(3,NPHO)+PPHO(2,I)*PPHO(2,NPHO)
00265      $                                   +PPHO(1,I)*PPHO(1,NPHO))/
00266      $            SQRT(PPHO(3,I)**2+PPHO(2,I)**2+PPHO(1,I)**2)   /
00267      $            SQRT(PPHO(3,NPHO)**2+PPHO(2,NPHO)**2+PPHO(1,NPHO)**2)
00268            MPASQR=PPHO(4,1)**2
00269            XPH=PPHO(4,NPHO)
00270 
00271 c...       Initialization of the W->l\nu\gamma 
00272 c...       decay Matrix Element parameters 
00273            IF (SANC_MC_INIT.EQ.-123456789D0) THEN
00274               SANC_MC_INIT=1D0
00275 
00276               PI=4*datan(1d0)
00277               QF1=0d0                           ! neutrino charge
00278               MF1=1d-10                         ! newutrino mass
00279               VF=1d0                            ! V&A couplings
00280               AF=1d0
00281               alphaI=1d0/ALPHA
00282               CW=0.881731727d0                  ! Weak Weinberg angle
00283               SW=0.471751166d0
00284            
00285 
00286 c...          An auxilary K&S vectors
00287               bet(0)= 1d0
00288               bet(1)= 0.0722794881816159d0
00289               bet(2)=-0.994200045099866d0
00290               bet(3)= 0.0796363353729248d0 
00291 
00292               spV(0)= 0d0 
00293               spV(1)= 7.22794881816159d-002
00294               spV(2)=-0.994200045099866d0     
00295               spV(3)= 7.96363353729248d-002
00296 
00297               mcLUN = PHLUN
00298            ENDIF 
00299 
00300 
00301            MB=PPHO(4,1)                      ! W boson mass
00302            MF2=dsqrt(MCHREN)                 ! muon mass
00303 
00304            DO IJ=1,d_h_nhep
00305             IF(ABS(d_h_idhep(IJ)).EQ.24) I3=IJ ! position of W 
00306            ENDDO
00307            IF(
00308      $      ABS(d_h_idhep(d_h_jdahep(1,I3)  )).EQ.11.OR.
00309      $      ABS(d_h_idhep(d_h_jdahep(1,I3)  )).EQ.13.OR.
00310      $      ABS(d_h_idhep(d_h_jdahep(1,I3)  )).EQ.15    ) THEN 
00311               I4=d_h_jdahep(1,I3)              ! position of lepton
00312            ELSE
00313               I4=d_h_jdahep(1,I3)+1            ! position of lepton
00314            ENDIF
00315 
00316 
00317               IF (d_h_idhep(I3).EQ.-24) QB=-1D0  ! W boson charge
00318               IF (d_h_idhep(I3).EQ.+24) QB=+1D0    
00319               IF (d_h_idhep(I4).GT.0d0) QF2=-1D0 ! lepton charge
00320               IF (d_h_idhep(I4).LT.0d0) QF2=+1D0
00321 
00322 
00323 c...          Particle momenta before foton radiation; effective Born level
00324               DO JJ=1,4
00325                 B_PW(mod(JJ,4))=d_h_phep(JJ,I3)  ! W boson
00326                 B_PNE(mod(JJ,4))=d_h_phep(JJ,I3)-d_h_phep(JJ,I4) ! neutrino
00327                 B_PMU(mod(JJ,4))=d_h_phep(JJ,I4) ! muon
00328               ENDDO
00329 
00330 c..        Particle monenta after photon radiation
00331            DO JJ=1,4
00332              PW(mod(JJ,4))=PPHO(JJ,1)
00333              PMU(mod(JJ,4))=PPHO(JJ,I)
00334              PPHOT(mod(JJ,4))=PPHO(JJ,NPHO)
00335              PNE(mod(JJ,4))=PPHO(JJ,1)-PPHO(JJ,I)-PPHO(JJ,NPHO)
00336            ENDDO
00337 C two options of calculating neutrino (spectator) mass
00338            MF1=SQRT(ABS(B_PNE(0)**2-B_PNE(1)**2-B_PNE(2)**2-B_PNE(3)**2))
00339            MF1=SQRT(ABS(  PNE(0)**2-  PNE(1)**2-  PNE(2)**2-  PNE(3)**2))
00340 c..        Exact amplitude square      
00341            AMPSQR=WDecayAmplitudeSqrKS_1ph(PW,PNE,PMU,PPHOT)
00342 
00343            EIKONALFACTOR=WDecayBornAmpSqrKS_1ph(B_PW,B_PNE,B_PMU)
00344      &                  *WDecayEikonalSqrKS_1ph(PW,PNE,PMU,PPHOT)
00345       
00346 c..        New weight
00347 !           write(*,*) 'B_pne=',B_PNE
00348 !           write(*,*) 'B_PMU=',B_PMU
00349 !           write(*,*) 'bornie=',WDecayBornAmpSqrKS_1ph(B_PW,B_PNE,B_PMU)
00350 
00351 !           write(*,*) ' '
00352 !           write(*,*) '  pne=',pne
00353 !           write(*,*) '  pmu=',pmu
00354 !           write(*,*) 'pphot=',pphot
00355 !           write(*,*) ' '
00356 !           write(*,*) '  b_pw=',B_PW
00357 !           write(*,*) '  b_pne=',B_PNE
00358 !           write(*,*) 'b_pmu=',B_PMU
00359  
00360  !          write(*,*) 'cori=',AMPSQR/EIKONALFACTOR,AMPSQR,EIKONALFACTOR
00361            WT=WT*AMPSQR/EIKONALFACTOR
00362 c           
00363 c          WT=WT*(1-8*EMU*XPH*(1-COSTHG*BETA)*     
00364 c     $           (MCHREN+2*XPH*SQRT(MPASQR))/
00365 c     $            MPASQR**2/(1-MCHREN/MPASQR)/(4-MCHREN/MPASQR)) 
00366         ENDIF
00367 !      write(*,*)   'AMPSQR/EIKONALFACTOR= ',   AMPSQR/EIKONALFACTOR
00368       END
00369 
00370 c$$$
00371 c$$$
00372 c$$$
00373 c$$$
00374 c$$$
00375 c$$$
00376 c$$$
00377 c$$$
00378 c$$$
00379 c$$$
00380 c$$$
00381 c$$$
00382 c$$$
00383 c$$$
00384 c$$$
00385 c$$$   WffGammaME.f  follows now: 
00386 c$$$
00387 c$$$
00388 c$$$
00389 c$$$
00390 
00391       FUNCTION WDecayAmplitudeSqrKS_1ph(p3,p1,p2,k)
00392 c========================================================
00393 c        The squared amplitude for W decay              =
00394 c        into fermion pair and one photon               =
00395 c  INPUT :                                              =
00396 c                                                       =
00397 c  OUTPUT:                                              =
00398 c========================================================       
00399          IMPLICIT NONE    
00400          INTEGER          l1,l2,l3,s,flag
00401          DOUBLE PRECISION k(0:3),p1(0:3),p2(0:3),p3(0:3)
00402          DOUBLE PRECISION spinSumAvrg,WDecayAmplitudeSqrKS_1ph
00403          DOUBLE COMPLEX   WDecayAmplitudeKS_1ph,wDecAmp
00404          EXTERNAL         WDecayAmplitudeKS_1ph 
00405          INTEGER           mcLUN
00406          DOUBLE PRECISION  spV(0:3),bet(0:3)
00407          DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
00408          COMMON /Kleiss_Stirling/spV,bet
00409          COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN        
00410 
00411          spinSumAvrg = 0d0
00412          DO l3 = -1,1        
00413            DO l1 = -1,1,2  
00414              DO l2 = -1,1,2  
00415                DO s = -1,1,2  
00416                  wDecAmp     = WDecayAmplitudeKS_1ph(p3,l3,p1,l1,p2,l2,k,s)
00417                  spinSumAvrg = spinSumAvrg + wDecAmp*CONJG(wDecAmp) 
00418                ENDDO
00419              ENDDO
00420            ENDDO
00421          ENDDO        
00422          WDecayAmplitudeSqrKS_1ph = spinSumAvrg
00423       
00424       RETURN
00425       END    
00426 
00427       FUNCTION WDecayBornAmpSqrKS_1ph(p3,p1,p2)
00428 c========================================================
00429 c        The squared eikonal factor for W decay         =
00430 c        into fermion pair and one photon               =
00431 c  INPUT :                                              =
00432 c                                                       =
00433 c  OUTPUT:                                              =
00434 c========================================================       
00435          IMPLICIT NONE    
00436          INTEGER          l1,l2,l3
00437          DOUBLE PRECISION k(0:3),p1(0:3),p2(0:3),p3(0:3)
00438          DOUBLE PRECISION spinSumAvrg,WDecayBornAmpSqrKS_1ph
00439          DOUBLE COMPLEX   WDecayBornAmpKS_1ph,wDecAmp
00440          EXTERNAL         WDecayBornAmpKS_1ph
00441          INTEGER           mcLUN
00442          DOUBLE PRECISION  spV(0:3),bet(0:3)
00443          DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
00444          COMMON /Kleiss_Stirling/spV,bet
00445          COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN        
00446 
00447          spinSumAvrg = 0d0
00448          DO l3 = -1,1      
00449            DO l1 = -1,1,2  
00450              DO l2 = -1,1,2  
00451                 wDecAmp     = WDecayBornAmpKS_1ph(p3,l3,p1,l1,p2,l2)
00452                 spinSumAvrg = spinSumAvrg + wDecAmp*CONJG(wDecAmp) 
00453              ENDDO
00454            ENDDO
00455          ENDDO        
00456          WDecayBornAmpSqrKS_1ph = spinSumAvrg
00457 
00458        RETURN
00459        END   
00460 
00461       FUNCTION WDecayEikonalSqrKS_1ph(p3,p1,p2,k)
00462 c========================================================
00463 c        The squared eikonal factor for W decay         =
00464 c        into fermion pair and one photon               =
00465 c  INPUT :                                              =
00466 c                                                       =
00467 c  OUTPUT:                                              =
00468 c========================================================       
00469          IMPLICIT NONE    
00470          INTEGER          s
00471          DOUBLE PRECISION k(0:3),p1(0:3),p2(0:3),p3(0:3)
00472          DOUBLE PRECISION spinSumAvrg,WDecayEikonalSqrKS_1ph
00473          DOUBLE COMPLEX   WDecayEikonalKS_1ph,wDecAmp
00474          EXTERNAL         WDecayEikonalKS_1ph
00475          INTEGER           mcLUN
00476          DOUBLE PRECISION  spV(0:3),bet(0:3)
00477          DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
00478          COMMON /Kleiss_Stirling/spV,bet
00479          COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN         
00480 
00481          spinSumAvrg = 0d0
00482          DO s = -1,1,2  
00483            wDecAmp     = WDecayEikonalKS_1ph(p3,p1,p2,k,s)
00484            spinSumAvrg = spinSumAvrg + wDecAmp*CONJG(wDecAmp) 
00485          ENDDO
00486          WDecayEikonalSqrKS_1ph = spinSumAvrg
00487 
00488        RETURN
00489        END   
00490               
00491 
00492        FUNCTION WDecayEikonalKS_1ph(p3,p1,p2,k,s)
00493 c======================================================================                 
00494 c     
00495 c Eikonal factor of decay W->l_1+l_2+\gamma in terms of K&S objects !
00496 c 
00497 c   EikFactor = q1*eps.p1/k.p1 + q2*eps.p2/k.p2 - q3*eps.p3/k.p3
00498 c
00499 c   indices 1,2 are for charged decay products
00500 c   index 3 is for W
00501 c   
00502 c   q - charge
00503 c    
00504 c======================================================================
00505          IMPLICIT NONE    
00506          INTEGER          s
00507          DOUBLE PRECISION k(0:3),p1(0:3),p2(0:3),p3(0:3)
00508          DOUBLE PRECISION scalProd1,scalProd2,scalProd3
00509          DOUBLE COMPLEX   WDecayEikonalKS_1ph,BsFactor,BSoft1,BSoft2  
00510          EXTERNAL         BsFactor
00511          INTEGER           mcLUN
00512          DOUBLE PRECISION  spV(0:3),bet(0:3)
00513          DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
00514          COMMON /Kleiss_Stirling/spV,bet
00515          COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN    
00516 
00517          scalProd1 = p1(0)*k(0)-p1(1)*k(1)-p1(2)*k(2)-p1(3)*k(3)
00518          scalProd2 = p2(0)*k(0)-p2(1)*k(1)-p2(2)*k(2)-p2(3)*k(3)
00519          scalProd3 = p3(0)*k(0)-p3(1)*k(1)-p3(2)*k(2)-p3(3)*k(3)
00520 
00521          BSoft1  = BsFactor(s,k,p1,mf1)
00522          BSoft2  = BsFactor(s,k,p2,mf2)
00523  
00524         WDecayEikonalKS_1ph = 
00525      &       sqrt(pi/alphaI)*(-(qf1/scalProd1+qb/scalProd3)*BSoft1   
00526      &                        +(qf2/scalProd2-qb/scalProd3)*BSoft2)
00527 
00528        RETURN
00529        END
00530 
00531        FUNCTION WDecayBornAmpKS_1ph(p3,l3,p1,l1,p2,l2)
00532 c======================================================================                 
00533 c                                                                     =
00534 c                            p1,mf1,l1                                =
00535 c                           /                                         =
00536 c                         \/_                                         = 
00537 c                         /                                           =
00538 c        p3,mb,l3        /                                            =
00539 c              \/\/\/\/\/\      ------> g_(mu,1)*(1+g5_(1))           =
00540 c                         \                                           =
00541 c                         _\/                                         = 
00542 c                           \                                         =
00543 c                            p2,mf2,l2                                =
00544 c INPUT : p1,m1,l1; p2,m2,l2; p3,m3,l3  -- momenta,mass and helicity  =
00545 c                                                                     =
00546 c OUTPUT: value of functions            -- decay amplitude            =
00547 c                                                                     =
00548 c======================================================================
00549          IMPLICIT NONE    
00550          INTEGER          l1,l2,l3
00551          DOUBLE PRECISION p1(0:3),p2(0:3),p3(0:3),coeff
00552          DOUBLE COMPLEX   WDecayBornAmpKS_1ph,TrMatrix_mass
00553          EXTERNAL         TrMatrix_mass
00554          INTEGER           mcLUN
00555          DOUBLE PRECISION  spV(0:3),bet(0:3)
00556          DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
00557          COMMON /Kleiss_Stirling/spV,bet
00558          COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN 
00559 
00560          coeff = Dsqrt(pi/alphaI/2d0)/sw      ! vertex: g/2/sqrt(2)
00561 
00562          WDecayBornAmpKS_1ph = coeff*TrMatrix_mass(p2,mf2,l2,p3,mb,l3,p1,-mf1,-l1)
00563          
00564 
00565        RETURN
00566        END
00567 
00568        FUNCTION WDecayAmplitudeKS_1ph(p3,l3,p1,l1,p2,l2,k,s)
00569 c======================================================================                 
00570 c                k,0,l                                                =
00571 c                   \        p1,mf1,l1                                =
00572 c                   /       /                                         =
00573 c                   \     \/_                                         = 
00574 c                   /     /                                           =
00575 c        p3,mb,l3   \    /                                            =
00576 c              \/\/\/\/\/\      ------> g_(mu,1)*(1+g5_(1))           =
00577 c                         \                                           =
00578 c                         _\/                                         = 
00579 c                           \                                         =
00580 c                            p2,mf2,l2                                =
00581 c           { + }                                                     =
00582 c                            p1,mf1,l1                                =
00583 c                           /                                         = 
00584 c                         \/_~~~~~~~ k,0,s                            = 
00585 c                         /                                           =
00586 c        p3,mb,l3        /                                            =
00587 c              \/\/\/\/\/\      ------> g_(mu,1)*(1+g5_(1))           =
00588 c                         \                                           =
00589 c                         _\/                                         = 
00590 c                           \                                         =
00591 c                            p2,mf2,l2                                =
00592 c           { + }                                                     =
00593 c                            p1,mf1,l1                                =
00594 c                           /                                         =
00595 c                         \/_                                         =
00596 c                         /                                           =
00597 c        p3,mb,l3        /                                            =
00598 c              \/\/\/\/\/\      ------> g_(mu,1)*(1+g5_(1))           =
00599 c                         \                                           =
00600 c                         _\/ ~~~~~~~ k,0,s                           =
00601 c                           \                                         =
00602 c                             p2,mf2,l2                               =
00603 c                                                                     =
00604 c                   all momentas, exept k are incoming !!!            =
00605 c                                                                     =
00606 c This function culculates The W-ff\gamma decay amplitude into permion=
00607 c pair and one photon using Kleisse&Stirling method for helicity      =
00608 c amplitudes, which includes three above feynman diagramms..          = 
00609 c                                                                     =
00610 c INPUT : p1,m1,l1; p2,m2,l2; p3,m3,l3  -- momenta,mass and helicity  =
00611 c                                                                     =
00612 c OUTPUT: value of functions            -- decay amplitude            =
00613 c                                                                     =
00614 c======================================================================
00615          IMPLICIT NONE    
00616          INTEGER          s,l1,l2,l3,la1
00617          DOUBLE PRECISION k(0:3),p1(0:3),p2(0:3),p3(0:3)
00618          DOUBLE PRECISION scalProd1,scalProd2,scalProd3,coeff,theta3,ph3
00619          DOUBLE COMPLEX   bornAmp,TrMx1,TrMx2,WDecayAmplitudeKS_1ph,eikAmp
00620          DOUBLE COMPLEX   BsFactor,TrMatrix_zero,TrMatrix_mass,BSoft1,BSoft2  
00621          EXTERNAL         BsFactor,TrMatrix_zero,TrMatrix_mass
00622          INTEGER           mcLUN
00623          DOUBLE PRECISION  spV(0:3),bet(0:3)
00624          DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
00625          COMMON /Kleiss_Stirling/spV,bet
00626          COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN
00627 
00628          coeff = Dsqrt(2d0)*pi/sw/alphaI      ! vertex: g/2/sqrt(2) * e
00629 
00630          scalProd1 = p1(0)*k(0)-p1(1)*k(1)-p1(2)*k(2)-p1(3)*k(3)
00631          scalProd2 = p2(0)*k(0)-p2(1)*k(1)-p2(2)*k(2)-p2(3)*k(3)
00632          scalProd3 = p3(0)*k(0)-p3(1)*k(1)-p3(2)*k(2)-p3(3)*k(3)
00633 
00634          BSoft1  = BsFactor(s,k,p1,mf1)
00635          BSoft2  = BsFactor(s,k,p2,mf2)
00636          bornAmp = TrMatrix_mass(p2,mf2,l2,p3,mb,l3,p1,-mf1,-l1)
00637          TrMx1   = Dcmplx(0d0,0d0)  
00638          TrMx2   = Dcmplx(0d0,0d0)  
00639  
00640          DO la1=-1,1,2            
00641            TrMx1 = TrMx1 + TrMatrix_zero(k,0d0,-la1,k,s,p1,-mf1,-l1)*
00642      &                     TrMatrix_mass(p2,mf2,l2,p3,mb,l3,k,0d0,-la1)
00643            TrMx2 = TrMx2 + TrMatrix_zero(p2,mf2,l2,k,s,k,0d0,la1)*
00644      &                     TrMatrix_mass(k,0d0,la1,p3,mb,l3,p1,-mf1,-l1)
00645          ENDDO
00646 
00647          WDecayAmplitudeKS_1ph = coeff * (        
00648      &        + (-(qf1/scalProd1+qb/scalProd3)*BSoft1              ! IR-divergent part of amplitude      
00649      &           +(qf2/scalProd2-qb/scalProd3)*BSoft2)/2d0*bornAmp
00650 
00651      &        - (qf1/ScalProd1+qb/scalProd3)*TrMx1/2d0             ! IR-finite part of amplitude            
00652      &        + (qf2/ScalProd2-qb/scalProd3)*TrMx2/2d0
00653      &        ) 
00654 
00655        RETURN
00656        END
00657  
00658        FUNCTION TrMatrix_mass(p1,m1,l1,k,m,s,p2,m2,l2)
00659 c//////////////////////////////////////////////////////////////
00660 c          transition matrix for massive boson               //
00661 c                                                            // 
00662 c                                                            //
00663 c                         \ eps(k,m,s)                       //
00664 c                         /                                  // 
00665 c                        _\                                  //
00666 c                         /\ k                               // 
00667 c                         \                                  //
00668 c             <-- p1      /         <-- p2                   //                       
00669 c           ---<----------\----------<---                    //
00670 c       Ub(p1,m1,l1)                  U(p2,m2,l2)            //
00671 c                                                            // 
00672 c//////////////////////////////////////////////////////////////                         
00673            IMPLICIT NONE
00674            INTEGER          l1,l2,s,i
00675            DOUBLE PRECISION forSqrt1,forSqrt2
00676            DOUBLE PRECISION m,m1,m2
00677            DOUBLE PRECISION k_1(0:3),k_2(0:3),p1(0:3),p2(0:3),k(0:3)
00678            DOUBLE PRECISION forSqrt3,forSqrt4,sqrt3,sqrt1,sqrt2,sqrt4
00679            DOUBLE COMPLEX   InProd_zero,inPr1,inPr2,inPr3,inPr4
00680            DOUBLE COMPLEX   TrMatrix_mass
00681            EXTERNAL         InProd_zero
00682            INTEGER           mcLUN
00683            DOUBLE PRECISION  spV(0:3),bet(0:3)
00684            DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
00685            COMMON /Kleiss_Stirling/spV,bet
00686            COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN      
00687 
00688            DO i=0,3
00689              k_1(i) = 1d0/2d0*(k(i) - m*spV(i))
00690              k_2(i) = 1d0/2d0*(k(i) + m*spV(i))                                
00691            ENDDO
00692 
00693            IF ( (l1.EQ.+1).AND.(l2.EQ.+1).AND.(s.EQ.0) ) THEN 
00694                 
00695               inPr1 = InProd_zero(p1,+1,k_2,-1)
00696               inPr2 = InProd_zero(p2,-1,k_2,+1)
00697               inPr3 = InProd_zero(p1,+1,k_1,-1)
00698               inPr4 = InProd_zero(p2,-1,k_1,+1)
00699               sqrt1 = Dsqrt(p1(0)-p1(1))
00700               sqrt2 = Dsqrt(p2(0)-p2(1))
00701               sqrt3 = m1*m2/sqrt1/sqrt2
00702 
00703               TrMatrix_mass =                 
00704      &                      (inPr1*inPr2-inPr3*inPr4)*(vf+af)/m 
00705      &          + (k_1(0)-k_2(0)-k_1(1)+k_2(1))*sqrt3*(vf-af)/m         
00706                  
00707            ELSEIF ( (l1.EQ.+1).AND.(l2.EQ.-1).AND.(s.EQ.0) ) THEN 
00708 
00709               inPr1 = InProd_zero(p1,+1,k_1,-1)
00710               inPr2 = InProd_zero(p1,+1,k_2,-1)
00711               inPr3 = InProd_zero(p2,+1,k_2,-1)
00712               inPr4 = InProd_zero(p2,+1,k_1,-1)
00713 
00714               forSqrt1 = (k_1(0)-k_1(1))/(p2(0)-p2(1))
00715               forSqrt2 = (k_2(0)-k_2(1))/(p2(0)-p2(1))
00716               forSqrt3 = (k_2(0)-k_2(1))/(p1(0)-p1(1))
00717               forSqrt4 = (k_1(0)-k_1(1))/(p1(0)-p1(1))
00718               sqrt1 = Dsqrt(forSqrt1)
00719               sqrt2 = Dsqrt(forSqrt2)
00720               sqrt3 = Dsqrt(forSqrt3)
00721               sqrt4 = Dsqrt(forSqrt4)     
00722 
00723               TrMatrix_mass = 
00724      &                 (inPr1*sqrt1 - inPr2*sqrt2)*(vf+af)*m2/m
00725      &               + (inPr3*sqrt3 - inPr4*sqrt4)*(vf-af)*m1/m
00726                         
00727            ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.+1).AND.(s.EQ.0) ) THEN 
00728 
00729               inPr1 = InProd_zero(p1,-1,k_1,+1)
00730               inPr2 = InProd_zero(p1,-1,k_2,+1)
00731               inPr3 = InProd_zero(p2,-1,k_2,+1)
00732               inPr4 = InProd_zero(p2,-1,k_1,+1)
00733 
00734               forSqrt1 = (k_1(0)-k_1(1))/(p2(0)-p2(1))
00735               forSqrt2 = (k_2(0)-k_2(1))/(p2(0)-p2(1))
00736               forSqrt3 = (k_2(0)-k_2(1))/(p1(0)-p1(1))
00737               forSqrt4 = (k_1(0)-k_1(1))/(p1(0)-p1(1))
00738               sqrt1 = Dsqrt(forSqrt1)
00739               sqrt2 = Dsqrt(forSqrt2)
00740               sqrt3 = Dsqrt(forSqrt3)
00741               sqrt4 = Dsqrt(forSqrt4)     
00742         
00743               TrMatrix_mass = 
00744      &                 (inPr1*sqrt1 - inPr2*sqrt2)*(vf-af)*m2/m
00745      &               + (inPr3*sqrt3 - inPr4*sqrt4)*(vf+af)*m1/m
00746 
00747            ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.-1).AND.(s.EQ.0) ) THEN 
00748 
00749               inPr1 = InProd_zero(p2,+1,k_2,-1)
00750               inPr2 = InProd_zero(p1,-1,k_2,+1)
00751               inPr3 = InProd_zero(p2,+1,k_1,-1)
00752               inPr4 = InProd_zero(p1,-1,k_1,+1)
00753               sqrt1 = Dsqrt(p1(0)-p1(1))
00754               sqrt2 = Dsqrt(p2(0)-p2(1))
00755               sqrt3 = m1*m2/sqrt1/sqrt2
00756 
00757               TrMatrix_mass =                    
00758      &                      (inPr1*inPr2 - inPr3*inPr4)*(vf-af)/m  
00759      &            + (k_1(0)-k_2(0)-k_1(1)+k_2(1))*sqrt3*(vf+af)/m
00760         
00761            ELSEIF ( (l1.EQ.+1).AND.(l2.EQ.+1).AND.(s.EQ.+1) ) THEN 
00762 
00763               inPr1 = InProd_zero(p1,+1,k_1,-1)
00764               inPr2 = InProd_zero(k_2,-1,p2,+1)
00765               inPr3 = inPr1*inPr2
00766 
00767               forSqrt1 = (k_1(0)-k_1(1))/(p1(0)-p1(1))                       
00768               forSqrt2 = (k_2(0)-k_2(1))/(p2(0)-p2(1))  
00769               sqrt1 = Dsqrt(forSqrt1)                   
00770               sqrt2 = Dsqrt(forSqrt2)                   
00771               sqrt3 = m1*m2*sqrt1*sqrt2
00772 
00773               TrMatrix_mass =
00774      &                 Dsqrt(2d0)/m*(inPr3*(vf+af)+sqrt3*(vf-af))
00775 
00776            ELSEIF ( (l1.EQ.+1).AND.(l2.EQ.-1).AND.(s.EQ.+1) ) THEN 
00777 
00778               inPr1 = InProd_zero(p1,+1,k_1,-1)
00779               inPr2 = InProd_zero(p2,+1,k_1,-1) 
00780 
00781               forSqrt1 = (k_2(0)-k_2(1))/(p2(0)-p2(1))                       
00782               forSqrt2 = (k_2(0)-k_2(1))/(p1(0)-p1(1))                       
00783               sqrt1 = m2*Dsqrt(forSqrt1)                   
00784               sqrt2 = m1*Dsqrt(forSqrt2)                                     
00785                      
00786               TrMatrix_mass =
00787      &                 Dsqrt(2d0)/m*( + inPr1*sqrt1*(vf+af)
00788      &                                - inPr2*sqrt2*(vf-af)
00789      &                                )
00790 
00791            ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.+1).AND.(s.EQ.+1) ) THEN 
00792 
00793               inPr1 = InProd_zero(k_2,-1,p2,+1)
00794               inPr2 = InProd_zero(k_2,-1,p1,+1)
00795 
00796               forSqrt1 = (k_1(0)-k_1(1))/(p1(0)-p1(1))                       
00797               forSqrt2 = (k_1(0)-k_1(1))/(p2(0)-p2(1))                       
00798               sqrt1 = m1*Dsqrt(forSqrt1)                   
00799               sqrt2 = m2*Dsqrt(forSqrt2)                                     
00800                      
00801               TrMatrix_mass =
00802      &                 Dsqrt(2d0)/m*( + inPr1*sqrt1*(vf+af)
00803      &                                - inPr2*sqrt2*(vf-af)
00804      &                                )
00805 
00806            ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.-1).AND.(s.EQ.+1) ) THEN 
00807 
00808               inPr1 = InProd_zero(p2,+1,k_1,-1)
00809               inPr2 = InProd_zero(k_2,-1,p1,+1)
00810               inPr3 = inPr1*inPr2
00811 
00812               forSqrt1 = (k_1(0)-k_1(1))/(p1(0)-p1(1))                       
00813               forSqrt2 = (k_2(0)-k_2(1))/(p2(0)-p2(1))  
00814               sqrt1 = Dsqrt(forSqrt1)                   
00815               sqrt2 = Dsqrt(forSqrt2)                   
00816               sqrt3 = m1*m2*sqrt1*sqrt2
00817 
00818               TrMatrix_mass = 
00819      &                 Dsqrt(2d0)/m*(inPr3*(vf-af)+sqrt3*(vf+af))
00820 
00821            ELSEIF ( (l1.EQ.+1).AND.(l2.EQ.+1).AND.(s.EQ.-1) ) THEN 
00822 
00823               inPr1 = InProd_zero(p2,-1,k_1,+1)
00824               inPr2 = InProd_zero(k_2,+1,p1,-1)
00825               inPr3 = inPr1*inPr2
00826 
00827               forSqrt1 = (k_1(0)-k_1(1))/(p1(0)-p1(1))                       
00828               forSqrt2 = (k_2(0)-k_2(1))/(p2(0)-p2(1))  
00829               sqrt1 = Dsqrt(forSqrt1)                   
00830               sqrt2 = Dsqrt(forSqrt2)                   
00831               sqrt3 = m1*m2*sqrt1*sqrt2
00832 
00833               TrMatrix_mass =               
00834      &                 Dsqrt(2d0)/m*(inPr3*(vf+af)+sqrt3*(vf-af))
00835 
00836            ELSEIF ( (l1.EQ.+1).AND.(l2.EQ.-1).AND.(s.EQ.-1) ) THEN 
00837 
00838               inPr1 = InProd_zero(k_2,+1,p2,-1)
00839               inPr2 = InProd_zero(k_2,+1,p1,-1)
00840 
00841               forSqrt1 = (k_1(0)-k_1(1))/(p1(0)-p1(1))                       
00842               forSqrt2 = (k_1(0)-k_1(1))/(p2(0)-p2(1))                       
00843               sqrt1 = m1*Dsqrt(forSqrt1)                   
00844               sqrt2 = m2*Dsqrt(forSqrt2)                                     
00845                      
00846               TrMatrix_mass =
00847      &                 Dsqrt(2d0)/m*(+ inPr1*sqrt1*(vf-af)
00848      &                               - inPr2*sqrt2*(vf+af)
00849      &                               )
00850 
00851            ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.+1).AND.(s.EQ.-1) ) THEN 
00852 
00853               inPr1 = InProd_zero(p1,-1,k_1,+1)
00854               inPr2 = InProd_zero(p2,-1,k_1,+1)
00855 
00856               forSqrt1 = (k_2(0)-k_2(1))/(p2(0)-p2(1))                       
00857               forSqrt2 = (k_2(0)-k_2(1))/(p1(0)-p1(1))                       
00858               sqrt1 = m2*Dsqrt(forSqrt1)                   
00859               sqrt2 = m1*Dsqrt(forSqrt2)                                     
00860                      
00861               TrMatrix_mass =
00862      &                 Dsqrt(2d0)/m*(+ inPr1*sqrt1*(vf-af)
00863      &                               - inPr2*sqrt2*(vf+af) 
00864      &                               )
00865     
00866            ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.-1).AND.(s.EQ.-1) ) THEN 
00867 
00868               inPr1 = InProd_zero(p1,-1,k_1,+1)
00869               inPr2 = InProd_zero(k_2,+1,p2,-1)
00870               inPr3 = inPr1*inPr2
00871 
00872               forSqrt1 = (k_1(0)-k_1(1))/(p1(0)-p1(1))                       
00873               forSqrt2 = (k_2(0)-k_2(1))/(p2(0)-p2(1))  
00874               sqrt1 = Dsqrt(forSqrt1)                   
00875               sqrt2 = Dsqrt(forSqrt2)                   
00876               sqrt3 = m1*m2*sqrt1*sqrt2
00877 
00878               TrMatrix_mass =  
00879      &                 Dsqrt(2d0)/m*(inPr3*(vf-af)+sqrt3*(vf+af))
00880 
00881            ELSE
00882 
00883               WRITE(MCLUN,*) " "             
00884               WRITE(MCLUN,*) " TrMatrix_mass: Wrong values for l1,l2,s:"
00885               WRITE(MCLUN,*) "          l1,l2 = -1,+1; s = -1,0,1 "
00886               WRITE(MCLUN,*) " "             
00887               STOP
00888 
00889            ENDIF 
00890          
00891         RETURN
00892         END
00893 
00894         FUNCTION TrMatrix_zero(p1,m1,l1,k,s,p2,m2,l2)
00895 c############################################################################# 
00896 c                                                                            #
00897 c                         \ eps(k,0,s)                                       # 
00898 c                         /                                                  #   
00899 c                        _\                                                  # 
00900 c                         /\                                                 #
00901 c                         \                                                  #
00902 c                         /                                                  #
00903 c           ---<----------\-------------<---                                 #
00904 c       Ub(p1,m1,l1)                  U(p2,m2,l2)                            #
00905 c                                                                            #
00906 c                                                                            #
00907 c             definition of arbitrary light-like vector beta!!               #
00908 c                                                                            #
00909 c              bet(0) = 1.d0                                                 #
00910 c              bet(1) = 1.d0                                                 #
00911 c              bet(2) = 0.d0      <==> bet == k0  expression becomes easy!!  #
00912 c              bet(3) = 0.d0                                                 #
00913 c#############################################################################
00914               IMPLICIT NONE
00915               INTEGER          l1,l2,s,i
00916               DOUBLE PRECISION m1,m2,forSqrt1,forSqrt2,p1(0:3)
00917               DOUBLE PRECISION p1_1(0:3),p2_1(0:3),k(0:3),p2(0:3)
00918               DOUBLE PRECISION sqrt1,sqrt2,scalProd1,scalProd2
00919               DOUBLE COMPLEX   InProd_zero,inPr1,inPr2,inPr3,TrMatrix_zero
00920               LOGICAL          equal
00921               EXTERNAL         InProd_zero             
00922               INTEGER           mcLUN
00923               DOUBLE PRECISION  spV(0:3),bet(0:3)
00924               DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
00925               COMMON /Kleiss_Stirling/spV,bet
00926               COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN               
00927               equal = .TRUE.    
00928               DO i=0,3
00929                 IF (p1(i).NE.p2(i))  equal = equal.AND..FALSE.
00930               ENDDO
00931                     
00932 
00933               IF ( (m1.EQ.m2).AND.(equal) ) THEN
00934 c..          
00935 c..             when:  p1=p2=p <=> m1=m2 TrMatrix_zero is diagonal
00936 c..               
00937                  IF ( (l1.EQ.+1).AND.(l2.EQ.+1) ) THEN 
00938 
00939                     inPr1    = InProd_zero(k,+s,p1,-s)
00940                     forSqrt1 = (p1(0)-p1(1))/(k(0)-k(1)) 
00941                     sqrt1    = Dsqrt(2d0*forSqrt1)
00942 
00943                     TrMatrix_zero = sqrt1*inPr1
00944                     GOTO 111  
00945  
00946                  ELSEIF ( (l1.EQ.+1).AND.(l2.EQ.-1) ) THEN                
00947 
00948                     TrMatrix_zero = Dcmplx(0d0,0d0)
00949                     GOTO 111  
00950 
00951                  ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.+1) ) THEN                
00952 
00953                     TrMatrix_zero = Dcmplx(0d0,0d0)
00954                     GOTO 111  
00955 
00956                  ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.-1) ) THEN                
00957 
00958                     inPr1    = InProd_zero(k,+s,p1,-s)
00959                     forSqrt1 = (p1(0)-p1(1))/(k(0)-k(1)) 
00960                     sqrt1    = Dsqrt(2d0*forSqrt1)
00961 
00962                     TrMatrix_zero = sqrt1*inPr1
00963                     GOTO 111  
00964           
00965                  ELSE 
00966         
00967                     WRITE(MCLUN,*) ""             
00968                     WRITE(MCLUN,*) " ERROR IN  TrMatrix_zero: "
00969                     WRITE(MCLUN,*) "       WRONG VALUES FOR l1,l2,s" 
00970                     WRITE(MCLUN,*) ""             
00971                     STOP
00972 
00973                  ENDIF       
00974 
00975               ENDIF
00976 
00977               IF ( (l1.EQ.+1).AND.(l2.EQ.+1).AND.(s.EQ.+1) ) THEN 
00978 
00979                  inPr1    = InProd_zero(k,+1,p1,-1)
00980                  forSqrt1 = (p2(0)-p2(1))/(k(0)-k(1))
00981                  sqrt1    = Dsqrt(2d0*forSqrt1)                   
00982  
00983                  TrMatrix_zero = sqrt1*inPr1
00984 
00985               ELSEIF ( (l1.EQ.+1).AND.(l2.EQ.-1).AND.(s.EQ.+1) ) THEN 
00986  
00987                  TrMatrix_zero = Dcmplx(0d0,0d0)
00988 
00989               ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.+1).AND.(s.EQ.+1) ) THEN 
00990   
00991                  forSqrt1 = (p1(0)-p1(1))/(p2(0)-p2(1))             
00992                  forSqrt2 = 1d0/forSqrt1
00993                  sqrt1    = Dsqrt(2d0*forSqrt1)                   
00994                  sqrt2    = Dsqrt(2d0*forSqrt2)                   
00995                      
00996                  TrMatrix_zero = Dcmplx(m2*sqrt1-m1*sqrt2,0d0)
00997 
00998               ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.-1).AND.(s.EQ.+1) ) THEN 
00999 
01000                  inPr1    = InProd_zero(k,+1,p2,-1)
01001                  forSqrt1 = (p1(0)-p1(1))/(k(0)-k(1))
01002                  sqrt1    = Dsqrt(2d0*forSqrt1)                   
01003   
01004                  TrMatrix_zero = inPr1*sqrt1
01005 
01006               ELSEIF ( (l1.EQ.+1).AND.(l2.EQ.+1).AND.(s.EQ.-1) ) THEN 
01007  
01008                  inPr1    = -InProd_zero(k,-1,p2,+1)
01009                  forSqrt1 = (p1(0)-p1(1))/(k(0)-k(1))
01010                  sqrt1    = Dsqrt(2d0*forSqrt1)                   
01011  
01012                  TrMatrix_zero = -sqrt1*inPr1
01013 
01014 
01015               ELSEIF ( (l1.EQ.+1).AND.(l2.EQ.-1).AND.(s.EQ.-1) ) THEN 
01016            
01017                  forSqrt1 = (p1(0)-p1(1))/(p2(0)-p2(1))     
01018                  forSqrt2 = 1d0/forSqrt1
01019                  sqrt1    = Dsqrt(2d0*forSqrt1)                   
01020                  sqrt2    = Dsqrt(2d0*forSqrt2)                   
01021                      
01022                  TrMatrix_zero = Dcmplx(m2*sqrt1-m1*sqrt2,0d0)
01023 
01024               ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.+1).AND.(s.EQ.-1) ) THEN 
01025 
01026                  TrMatrix_zero = Dcmplx(0d0,0d0)
01027 
01028               ELSEIF ( (l1.EQ.-1).AND.(l2.EQ.-1).AND.(s.EQ.-1) ) THEN 
01029 
01030                  inPr1    = -InProd_zero(k,-1,p1,+1)
01031                  forSqrt1 = (p2(0)-p2(1))/(k(0)-k(1))
01032                  sqrt1    = Dsqrt(2d0*forSqrt1)                   
01033   
01034                  TrMatrix_zero = -inPr1*sqrt1
01035 
01036               ELSE
01037 
01038                  WRITE(MCLUN,*) ""
01039                  WRITE(MCLUN,*) " ERROR IN TrMatrix_zero: "
01040                  WRITE(MCLUN,*) "    WRONG VALUES FOR l1,l2,s"
01041                  WRITE(MCLUN,*) ""             
01042                  STOP
01043 
01044               ENDIF 
01045 
01046  111          CONTINUE           
01047 
01048        RETURN
01049        END
01050 
01051 
01052 
01053        DOUBLE PRECISION FUNCTION InSqrt(p,q)
01054 
01055            DOUBLE PRECISION q(0:3),p(0:3)
01056             
01057            InSqrt = Dsqrt( (p(0)-p(1)) / (q(0)-q(1)) )
01058 
01059        RETURN
01060        END    
01061 
01062 
01063        DOUBLE COMPLEX FUNCTION BsFactor(s,k,p,m)
01064 c///////////////////////////////////////////////////////////////////
01065 c                                                                 //
01066 c  this is small B_{s}(k,p) function when TrMartix is diaagonal!! //
01067 c                                                                 //
01068 c///////////////////////////////////////////////////////////////////
01069           IMPLICIT NONE
01070           INTEGER          s
01071           DOUBLE PRECISION p_1(0:3),p(0:3),k(0:3),m
01072           DOUBLE PRECISION forSqrt1,sqrt1
01073           DOUBLE COMPLEX   InProd_zero,inPr1
01074           EXTERNAL         InProd_zero             
01075           INTEGER           mcLUN
01076           DOUBLE PRECISION  spV(0:3),bet(0:3)
01077           DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
01078           COMMON /Kleiss_Stirling/spV,bet
01079           COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN            
01080 
01081           IF ( s.EQ.+1 ) THEN 
01082 
01083              inPr1    = InProd_zero(k,+1,p,-1)
01084              forSqrt1 = (p(0)-p(1))/(k(0)-k(1))
01085              sqrt1    = Dsqrt(2d0*forSqrt1)  
01086              BsFactor = inPr1*sqrt1
01087 
01088           ELSEIF ( s.EQ.-1 ) THEN 
01089 
01090              inPr1    = InProd_zero(k,-1,p,+1)
01091              forSqrt1 = (p(0)-p(1))/(k(0)-k(1))
01092              sqrt1    = Dsqrt(2d0*forSqrt1)  
01093              BsFactor = inPr1*sqrt1
01094 
01095          ELSE
01096 
01097             WRITE(MCLUN,*) " "             
01098             WRITE(MCLUN,*) " ERROR IN BsFactor: "
01099             WRITE(MCLUN,*) "       WRONG VALUES FOR s : s = -1,+1"
01100             WRITE(MCLUN,*) " "             
01101             STOP
01102 
01103          ENDIF 
01104 
01105        RETURN
01106        END   
01107 
01108        DOUBLE COMPLEX FUNCTION SoftFactor(s,k,p1,m1,p2,m2,Gmass2)
01109 
01110 c
01111 c       Gauge invariant soft factor for decay!!
01112 c       Gmass2 -- photon mass square       
01113 c 
01114           IMPLICIT NONE
01115           INTEGER          s
01116           DOUBLE PRECISION p1(0:3),p2(0:3),k(0:3),Gmass2
01117           DOUBLE PRECISION m1,m2,ScalProd1,ScalProd2
01118           DOUBLE COMPLEX   BsFactor2,BsFactor1,BsFactor
01119           EXTERNAL         BsFactor             
01120 
01121           ScalProd1 = k(0)*p1(0)-k(1)*p1(1)-k(2)*p1(2)-k(3)*p1(3)
01122           ScalProd2 = k(0)*p2(0)-k(1)*p2(1)-k(2)*p2(2)-k(3)*p2(3)
01123           
01124           BsFactor1 = BsFactor(s,k,p1,m1)
01125           BsFactor2 = BsFactor(s,k,p2,m2)
01126 
01127           SoftFactor= + BsFactor2/2.d0/(ScalProd2-Gmass2)
01128      &                - BsFactor1/2.d0/(ScalProd1-Gmass2)
01129 
01130        RETURN
01131        END
01132        DOUBLE COMPLEX FUNCTION InProd_mass(p1,m1,l1,p2,m2,l2)
01133 c////////////////////////////////////////////////////////////////
01134 c                                                              //
01135 c  Inner product for massive spinors: Ub(p1,m1,l1)*U(p2,m2,l2) //
01136 c                                                              //
01137 c////////////////////////////////////////////////////////////////
01138               
01139           IMPLICIT NONE
01140           INTEGER          l1,l2
01141           DOUBLE PRECISION p1(0:3),p2(0:3),m1,m2
01142           DOUBLE PRECISION sqrt1,sqrt2,forSqrt1
01143           DOUBLE COMPLEX   InProd_zero
01144           EXTERNAL         InProd_zero
01145           INTEGER           mcLUN
01146           DOUBLE PRECISION  spV(0:3),bet(0:3)
01147           DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
01148           COMMON /Kleiss_Stirling/spV,bet
01149           COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN        
01150 
01151           IF ((l1.EQ.+1).AND.(l2.EQ.+1)) THEN                
01152              forSqrt1    = (p1(0)-p1(1))/(p2(0)-p2(1))
01153              sqrt1       = dsqrt(forSqrt1)
01154              sqrt2       = 1d0/sqrt1
01155              InProd_mass = Dcmplx(m1*sqrt2+m2*sqrt1,0d0)
01156 
01157           ELSEIF ((l1.EQ.+1).AND.(l2.EQ.-1)) THEN                             
01158              InProd_mass = InProd_zero(p1,+1,p2,-1)
01159 
01160           ELSEIF ((l1.EQ.-1).AND.(l2.EQ.+1)) THEN                        
01161              InProd_mass = InProd_zero(p1,-1,p2,+1)               
01162 
01163           ELSEIF ((l1.EQ.-1).AND.(l2.EQ.-1)) THEN                             
01164              forSqrt1    = (p1(0)-p1(1))/(p2(0)-p2(1))
01165              sqrt1       = dsqrt(forSqrt1)
01166              sqrt2       = 1d0/sqrt1
01167              InProd_mass = Dcmplx(m1*sqrt2+m2*sqrt1,0d0)
01168 
01169           ELSE        
01170              WRITE(MCLUN,*) " "             
01171              WRITE(MCLUN,*) " ERROR IN InProd_mass.."
01172              WRITE(MCLUN,*) "       WRONG VALUES FOR l1,l2"
01173              WRITE(MCLUN,*) " "             
01174              STOP
01175           ENDIF       
01176 
01177        RETURN
01178        END
01179 
01180        DOUBLE COMPLEX FUNCTION InProd_zero(p1,l1,p2,l2)
01181 c////////////////////////////////////////////////////////////////
01182 c         small s_{+,-}(p1,p2) for massless case:              //
01183 c                 p1^2 = p2^2 = 0                              // 
01184 c                                                              //
01185 c     k0(0) = 1.d0                                             //
01186 c     k0(1) = 1.d0                                             //
01187 c     k0(2) = 0.d0  Kleisse_Stirling k0 points to X-axis       // 
01188 c     k0(3) = 0.d0                                             //
01189 c                                                              //
01190 c////////////////////////////////////////////////////////////////
01191            IMPLICIT NONE
01192            INTEGER           l1,l2,i    
01193            DOUBLE PRECISION  p1(0:3),p2(0:3)
01194            DOUBLE PRECISION  forSqrt1,forSqrt2,sqrt1,sqrt2
01195            DOUBLE COMPLEX    i_,Dcmplx
01196            LOGICAL           equal
01197            PARAMETER         ( i_=(0.d0,1.d0) )
01198            INTEGER           mcLUN
01199            DOUBLE PRECISION  spV(0:3),bet(0:3)
01200            DOUBLE PRECISION  pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af
01201            COMMON /Kleiss_Stirling/spV,bet
01202            COMMON /mc_parameters/pi,sw,cw,alphaI,qb,mb,mf1,mf2,qf1,qf2,vf,af,mcLUN                          
01203 
01204            equal = .TRUE.    
01205            DO i=0,3
01206              IF (p1(i).NE.p2(i))  equal = equal.AND..FALSE.                
01207            ENDDO
01208 
01209            IF ( (l1.EQ.l2) .OR. equal ) THEN
01210 
01211               InProd_zero = Dcmplx(0d0,0d0)
01212 
01213            ELSEIF ( (l1.EQ.+1) .AND. (l2.EQ.-1) ) THEN
01214 
01215               forSqrt1 = (p1(0)-p1(1))/(p2(0)-p2(1))
01216               forSqrt2 = 1.0d0/forSqrt1
01217               sqrt1    = Dsqrt(forSqrt2)
01218               sqrt2    = Dsqrt(forSqrt1)
01219 
01220               InProd_zero = (p1(2)+i_*p1(3))*sqrt1 -
01221      &                      (p2(2)+i_*p2(3))*sqrt2 
01222 
01223            ELSEIF ( (l1.EQ.-1) .AND. (l2.EQ.+1) ) THEN
01224 
01225               forSqrt1 = (p1(0)-p1(1))/(p2(0)-p2(1))
01226               forSqrt2 = 1.0d0/forSqrt1
01227               sqrt1    = Dsqrt(forSqrt2)
01228               sqrt2    = Dsqrt(forSqrt1)
01229 
01230               InProd_zero = (p2(2)-i_*p2(3))*sqrt2 -
01231      &                      (p1(2)-i_*p1(3))*sqrt1 
01232                                     
01233            ELSE
01234                  
01235               WRITE(MCLUN,*) " "             
01236               WRITE(MCLUN,*) " ERROR IN InProd_zero:"
01237               WRITE(MCLUN,*) "    WRONG VALUES FOR l1,l2: l1,l2 = -1,+1"
01238               WRITE(MCLUN,*) " "             
01239               STOP
01240 
01241            ENDIF   
01242 
01243         RETURN
01244         END