forZ-ME.f

00001       function phwtnlo(xdumm) 
00002 C.----------------------------------------------------------------------
00003 C.
00004 C.    PHWTNLO:   PHotosWTatNLO
00005 C.
00006 C.    Purpose:  calculates instead of interference weight
00007 C.              complete NLO weight for vector boson decays
00008 C.              of pure vector (or pseudovector) couplings
00009 C.              Proper orientation of beams required.
00010 C.              This is not standard in PHOTOS.
00011 C.              At NLO more tuning than in standard is needed.
00012 C.              Works with KORALZ and KKMC. 
00013 C.              Note some commented out commons from MUSTAAL, KORALZ
00014 C.
00015 C.    Input Parameters:   Common /PHOEVT/
00016 C.
00017 C.    Output Parameters:  Function value
00018 C.
00019 C.    Author(s):  Z. Was                          Created at:  08/12/05
00020 C.                                                Last Update: 
00021 C.
00022 C.----------------------------------------------------------------------
00023 
00024       IMPLICIT NONE 
00025       INTEGER NMXHEP
00026       PARAMETER (NMXHEP=10000)
00027       COMMON/PHOEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
00028      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
00029       real*8 phep,vhep
00030       INTEGER NEVHEP,NHEP,ISTHEP,IDHEP,JMOHEP,JDAHEP
00031       SAVE  /PHOEVT/
00032 C      COMMON / UTIL3 / xFIG,xFI,IT,IBREX
00033       DOUBLE PRECISION COSTHG,SINTHG
00034       REAL*8 XPHMAX,XPHOTO
00035       COMMON/PHOPHS/XPHMAX,XPHOTO,COSTHG,SINTHG
00036 C fi3 orientation of photon, fi1,th1 orientation of neutral
00037       REAL*8 FI3,fi1,th1
00038       COMMON /PHOREST/ FI3,fi1,th1
00039       REAL*8 BETA,WT1,WT2,WT3
00040       COMMON /PHWT/ BETA,WT1,WT2,WT3
00041       REAL*8 PROBH,CORWT,XF
00042       COMMON/PHOPRO/PROBH,CORWT,XF,IREP
00043       REAL*8 PI
00044       DATA PI /3.141592653589793238462643D0/
00045       REAL*8  QP(4),QM(4),PH(4),QQ(4),PP(4),PM(4),QQS(4)
00046       REAL*8 s,c,svar,xkaM,xkaP,xk,phwtnlo,xdumm,PHINT
00047       REAL*8 ENE,a,t,u,t1,u1,wagan2,waga,PHASYZ,BT,BU,ENEB
00048       INTEGER IBREM,K,L,IREP,IDUM
00049       integer icont,ide,idf
00050       REAL*8 delta
00051       REAL*8 PNEUTR,MCHSQR,MNESQR
00052       COMMON/PHOMOM/MCHSQR,MNESQR,PNEUTR(5)
00053 
00054 */////////////////////
00055 !        call phlupa(299 500)
00056 
00057 
00058 */////////////////////
00059 !        call phlupa(299 500)
00060 
00061         XK=2*PHEP(4,nhep)/PHEP(4,1)
00062 
00063 !        XK=2*PHEP(4,nhep)/PHEP(4,1)/xphmax  ! it is not used becuse here
00064                                              ! order of emissions is meaningless
00065       IF(NHEP.LE.4) XK=0.0D0
00066 C the mother must be Z or gamma*  !!!!
00067       
00068       IF (XK.GT.1d-10.AND.(IDHEP(1).EQ.22.OR.IDHEP(1).EQ.23)) THEN
00069 
00070 !        write(*,*) 'nhep=',nhep
00071 !      DO K=1,3 ENDDO
00072 !      IF (K.EQ.1) IBREM= 1
00073 !      IF (K.EQ.2) IBREM=-1
00074 !      ICONT=ICONT+1
00075 !      IBREM=IBREX        ! that will be input parameter.
00076 !      IBREM=IBREY        ! that IS now   input parameter.
00077 
00078 C We initialize twice 4-vectors, here and again later after boost 
00079 C must be the same way. Important is how the reduction procedure will work.
00080 C It seems at present that the beams must be translated to be back to back.
00081 C this may be done after initialising, thus on 4-vectors.
00082 
00083         DO K=1,4
00084          PP(K)=PHEP(K,1)
00085          PM(K)=PHEP(K,2)
00086          QP(K)=PHEP(K,3)
00087          QM(K)=PHEP(K,4)
00088          PH(K)=PHEP(K,nhep)
00089          QQ(K)=0.0
00090          QQS(K)=QP(K)+QM(K)
00091         ENDDO
00092 
00093 
00094         PP(4)=(PHEP(4,1)+PHEP(4,2))/2
00095         PM(4)=(PHEP(4,1)+PHEP(4,2))/2
00096         PP(3)= PP(4)
00097         PM(3)=-PP(4)
00098         
00099          DO L=5,NHEP-1
00100          DO K=1,4
00101           QQ(K)=QQ(K)+ PHEP(K,L)
00102          QQS(K)=QQS(K)+ PHEP(K,L)
00103          ENDDO
00104         ENDDO        
00105 
00106 C go to the restframe of 3        
00107         CALL PHOB(1,QQS,qp)
00108         CALL PHOB(1,QQS,qm)
00109         CALL PHOB(1,QQS,QQ)
00110         ENE=(qp(4)+qm(4)+QQ(4))/2
00111 
00112 C preserve direction of emitting particle and wipeout QQ 
00113         IF (IREP.EQ.1) THEN
00114          a=sqrt(ENE**2-PHEP(5,3)**2)/sqrt(qm(4)**2-PHEP(5,3)**2)
00115          qm(1)= qm(1)*a
00116          qm(2)= qm(2)*a
00117          qm(3)= qm(3)*a
00118          qp(1)=-qm(1)
00119          qp(2)=-qm(2)
00120          qp(3)=-qm(3)
00121         ELSE
00122          a=sqrt(ENE**2-PHEP(5,3)**2)/sqrt(qp(4)**2-PHEP(5,3)**2)
00123          qP(1)= qP(1)*a
00124          qP(2)= qP(2)*a
00125          qP(3)= qP(3)*a
00126          qm(1)=-qP(1)
00127          qm(2)=-qP(2)
00128          qm(3)=-qP(3)
00129         ENDIF
00130          qp(4)=ENE
00131          qm(4)=ENE
00132 C go back to reaction frame (QQ eliminated) 
00133         CALL PHOB(-1,QQS,qp)
00134         CALL PHOB(-1,QQS,qm)
00135         CALL PHOB(-1,QQS,QQ)
00136 C IBREM is spurious but it numbers branches of MUSTRAAL
00137 
00138         IBREM=1
00139        IF (IREP.EQ.1)  IBREM=-1
00140 
00141         svar=PHEP(4,1)**2
00142 
00143 C we calculate C and S, note that TH1 exists in MUSTRAAL as well. 
00144 
00145         C=COS(TH1)
00146 C from off line application we had:
00147         IF(IBREM.EQ.-1) C=-C
00148 C ... we need to check it. 
00149         s=sqrt(1D0-C**2)
00150 
00151         IF (IBREM.EQ.1) THEN
00152          xkaM=(qp(4)*PH(4)-qp(3)*PH(3)-qp(2)*PH(2)-qp(1)*PH(1))/XK
00153          xkaP=(qM(4)*PH(4)-qM(3)*PH(3)-qM(2)*PH(2)-qM(1)*PH(1))/XK
00154         ELSE
00155          xkap=(qp(4)*PH(4)-qp(3)*PH(3)-qp(2)*PH(2)-qp(1)*PH(1))/XK
00156          xkaM=(qM(4)*PH(4)-qM(3)*PH(3)-qM(2)*PH(2)-qM(1)*PH(1))/XK
00157         ENDIF    
00158 
00159 !        XK=2*PHEP(4,nhep)/PHEP(4,1)/xphmax   ! it is not used becuse here
00160                                               ! order of emissions is meaningless
00161 !
00162 !        DELTA=2*PHEP(5,4)**2/svar/(1+(1-XK)**2)*(xKAP/xKAM+xKAM/xKAP)
00163 !        waga=SVAR/4./xkap
00164 !        waga=waga*(1.D0-COSTHG*BETA) ! sprawdzone 1= svar/xKAp/4   * (1.D0-COSTHG*BETA)
00165 !        waga=waga*(1-delta) /wt2 ! sprawdzone ze to jest =2/(1.D0+COSTHG*BETA)
00166 !                                 ! czyli ubija de-interferencje
00167  
00168 
00169 C this is true only for intermediate resonances with afb=0!
00170         t =2*(qp(4)*PP(4)-qp(3)*PP(3)-qp(2)*PP(2)-qp(1)*PP(1))
00171         u =2*(qM(4)*PP(4)-qM(3)*PP(3)-qM(2)*PP(2)-qM(1)*PP(1))
00172         u1=2*(qp(4)*Pm(4)-qp(3)*Pm(3)-qp(2)*Pm(2)-qp(1)*Pm(1))
00173         t1=2*(qM(4)*Pm(4)-qM(3)*Pm(3)-qM(2)*Pm(2)-qM(1)*Pm(1))
00174 
00175 C basically irrelevant lines  ...
00176         t =t - (qp(4)**2-qp(3)**2-qp(2)**2-qp(1)**2)
00177         u =u - (qm(4)**2-qm(3)**2-qm(2)**2-qm(1)**2)
00178         u1=u1- (qp(4)**2-qp(3)**2-qp(2)**2-qp(1)**2)
00179         t1=t1- (qm(4)**2-qm(3)**2-qm(2)**2-qm(1)**2)
00180 
00181 
00182 
00183 
00184       call GETIDEIDF(IDE,IDF)   ! we adjust to what is f-st,s-nd beam flavour 
00185        IF (IDE*IDHEP(3).GT.0) THEN
00186         BT=1+PHASYZ(SVAR)
00187         BU=1-PHASYZ(SVAR)
00188        ELSE
00189         BT=1-PHASYZ(SVAR)
00190         BU=1+PHASYZ(SVAR)
00191        ENDIF  
00192         wagan2=2*(BT*t**2+BU*u**2+BT*t1**2+BU*u1**2)
00193      $        /(1+(1-xk)**2)* 2.0/(BT*(1-c)**2+BU*(1+c)**2)/svar**2
00194 
00195 !        waga=waga*wagan2
00196 !        waga=waga*(1-delta) /wt2 ! sprawdzone ze to jest =2/(1.D0+COSTHG*BETA)
00197         waga=2/(1.D0+COSTHG*BETA)*wagan2  
00198 !     %       * SVAR/4./xkap*(1.D0-COSTHG*BETA)*sqrt(1.0-xk)
00199 
00200         phwtnlo=waga
00201         IF(wagan2.gt.3.8) THEN
00202 !         write(*,*) 'phwtnlo= ',phwtnlo
00203          write(*,*) 'idhepy= ',IDHEP(1),IDHEP(2),IDHEP(3),IDHEP(4),IDHEP(5)
00204          write(*,*) 'IDE=',IDE,'  IDF=',IDF
00205          write(*,*) 'bt,bu,bt+bu= ',bt,bu,bt+bu
00206          call PHODMP
00207          write(*,*) ' '
00208          write (*,*) IREP,IBREM, '<-- IREP,IBREM '
00209          write(*,*) 'pneutr= ',pneutr
00210          write(*,*) 'qp    = ',qp
00211          write(*,*) 'qm    = ',qm
00212          write(*,*) ' '
00213          write(*,*) 'ph    = ',ph
00214          write(*,*) 'p1= ',PHEP(1,1),PHEP(2,1),PHEP(3,1),PHEP(4,1)
00215          write(*,*) 'p2= ',PHEP(1,2),PHEP(2,2),PHEP(3,2),PHEP(4,2)
00216          write(*,*) 'p3= ',PHEP(1,3),PHEP(2,3),PHEP(3,3),PHEP(4,3)
00217          write(*,*) 'p4= ',PHEP(1,4),PHEP(2,4),PHEP(3,4),PHEP(4,4)
00218          write(*,*) 'p5= ',PHEP(1,5),PHEP(2,5),PHEP(3,5),PHEP(4,5)
00219 
00220          write (*,*) ' c= ',c,' theta= ',th1
00221 !         write(*,*)  'photos waga daje ... IBREM=',IBREM,' waga=',waga
00222 !         write(*,*) 'xk,COSTHG,c',xk,COSTHG,c
00223 !         write(*,*) SVAR/4./xkap*(1.D0-COSTHG*BETA), 
00224 !     $   (1-delta) /wt2 *(1.D0+COSTHG*BETA)/2, wagan2
00225 !         write(*,*) ' delta, wt2,beta',  delta, wt2,beta
00226          write(*,*) '  -  '
00227          write(*,*) 't,u       = ',t,u
00228          write(*,*) 't1,u1     = ',t1,u1
00229          write(*,*) 'sredniaki = ',svar*(1-c)/2
00230      $                            ,svar*(1+c)/2
00231 !         write(*,*) 'xk= ',xk,' c= ',c ,' COSTHG= ' ,COSTHG
00232          write(*,*) 'PHASYZ(SVAR)=',PHASYZ(SVAR),' svar= ',svar,' waga=',waga
00233          write(*,*) '  -  '
00234          write(*,*) 'BT-part= ',2*(BT*t**2+BT*t1**2)
00235      $        /(1+(1-xk)**2)* 2.0/(BT*(1-c)**2)/svar**2,
00236      $ ' BU-part= ',2*(BU*u**2+BU*u1**2)
00237      $        /(1+(1-xk)**2)* 2.0/(BU*(1+c)**2)/svar**2
00238          write(*,*) 'BT-part*BU-part= ',2*(BT*t**2+BT*t1**2)
00239      $        /(1+(1-xk)**2)* 2.0/(BT*(1-c)**2)/svar**2
00240      $         *2*(BU*u**2+BU*u1**2)
00241      $        /(1+(1-xk)**2)* 2.0/(BU*(1+c)**2)/svar**2, 'wagan2=',wagan2
00242 
00243          write(*,*) 'wagan2=',wagan2
00244          write(*,*) ' ###################  '
00245          wagan2=3.8  ! overwrite 
00246          waga=2/(1.D0+COSTHG*BETA)*wagan2  
00247 !     %       * SVAR/4./xkap*(1.D0-COSTHG*BETA)*sqrt(1.0-xk)
00248 
00249         phwtnlo=waga
00250 
00251         ENDIF
00252       ELSE
00253 C in other cases we just use default setups.
00254         phwtnlo= PHINT(IDUM)
00255       ENDIF
00256       end
00257       
00258       FUNCTION PHASYZ(SVAR)
00259 C.----------------------------------------------------------------------
00260 C.
00261 C.    PHWTNLO:   PHotosASYmmetryofZ
00262 C.
00263 C.    Purpose:  Calculates born level asymmetry for Z
00264 C.              between distributions (1-C)**2 and (1+C)**2
00265 C.              At present dummy, requrires effective Z and gamma 
00266 C.              Couplings and also spin polarization states
00267 C.              For initial and final states.
00268 C.              To be correct this function need to be tuned
00269 C.              to host generator. Axes orientation polarisation
00270 C.              conventions etc etc. 
00271 C.              Modularity of PHOTOS would break. 
00272 C.
00273 C.    Input Parameters:   SVAR
00274 C.
00275 C.    Output Parameters:  Function value
00276 C.
00277 C.    Author(s):  Z. Was                          Created at:  10/12/05
00278 C.                                                Last Update: 
00279 C.
00280 C.----------------------------------------------------------------------
00281 
00282       IMPLICIT NONE
00283       real*8 PHASYZ,SVAR,AFB,AFBCALC
00284       INTEGER IDE,IDF,IDEE,IDFF,GETIDEE
00285       call GETIDEIDF(IDE,IDF)
00286       IDEE=abs(GETIDEE(IDE))
00287       IDFF=abs(GETIDEE(IDF))
00288       AFB= -AFBCALC(SVAR,IDEE,IDFF)
00289 C      AFB=0
00290       PHASYZ=4.D0/3.D0*AFB
00291 C      write(*,*) 'IDE=',IDE,'  IDF=',IDF,'  SVAR=',SVAR,'AFB=',AFB
00292       END
00293 
00294       FUNCTION GETIDEE(IDE)
00295       IMPLICIT NONE
00296       INTEGER IDEE,IDE, GETIDEE
00297       IF(IDE.EQ.11.OR.IDE.EQ.13.OR.IDE.EQ.15) THEN
00298         IDEE=2
00299       ELSEIF(IDE.EQ.-11.OR.IDE.EQ.-13.OR.IDE.EQ.-15) THEN
00300         IDEE=-2
00301       ELSEIF(IDE.EQ.12.OR.IDE.EQ.14.OR.IDE.EQ.16) THEN
00302          IDEE=1
00303       ELSEIF(IDE.EQ.-12.OR.IDE.EQ.-14.OR.IDE.EQ.-16) THEN
00304          IDEE=-1
00305       ELSEIF(IDE.EQ.1.OR.IDE.EQ.3.OR.IDE.EQ.5) THEN
00306          IDEE=4
00307       ELSEIF(IDE.EQ.-1.OR.IDE.EQ.-3 .OR.IDE.EQ.-5 ) THEN
00308          IDEE=-4
00309       ELSEIF(IDE.EQ. 2.OR.IDE.EQ. 4.OR.IDE.EQ. 6) THEN
00310          IDEE=3
00311       ELSEIF(IDE.EQ.- 2.OR.IDE.EQ.- 4.OR.IDE.EQ.- 6) THEN
00312          IDEE=-3
00313       ENDIF
00314       getidee=idee
00315       END
00316 
00317       SUBROUTINE PHOB(MODE,PBOOS1,VEC)
00318 C.----------------------------------------------------------------------
00319 C.
00320 C.
00321 C.    PHOB:     PHotosBoost
00322 C.
00323 C.    Purpose:  Boosts VEC to (MODE=1)  rest frame of PBOOS1;  
00324 C.              or back (MODE=1)
00325 C.
00326 C.    Input Parameters:   MODE,PBOOS1,VEC
00327 C.
00328 C.    Output Parameters:  VEC
00329 C.
00330 C.    Author(s):                                  Created at:  08/12/05
00331 C.                Z. Was                          Last Update: 
00332 C.
00333 C.----------------------------------------------------------------------
00334       IMPLICIT NONE
00335       DOUBLE PRECISION BET1(3),GAM1,PB
00336       INTEGER I,J,MODE
00337       REAL*8 PBOOS1(4),vec(4)
00338 
00339       PB=sqrt(PBOOS1(4)**2-PBOOS1(3)**2-PBOOS1(2)**2-PBOOS1(1)**2)
00340       DO 10 J=1,3
00341         IF (MODE.EQ.1) THEN
00342           BET1(J)=-PBOOS1(J)/PB
00343         ELSE
00344           BET1(J)= PBOOS1(J)/PB
00345         ENDIF 
00346   10  CONTINUE
00347       GAM1=PBOOS1(4)/PB
00348 
00349 C--
00350 C--   Boost vector 
00351 
00352         PB=BET1(1)*vec(1)+BET1(2)*vec(2)+BET1(3)*vec(3)
00353         
00354          DO 30 J=1,3
00355    30    vec(J)=vec(J)+BET1(J)*(vec(4)+PB/(GAM1+1.D0))
00356          vec(4)=GAM1*vec(4)+PB
00357 
00358 C--
00359 
00360       RETURN
00361       END
00362 
00363 
00364 
00365       SUBROUTINE GETIDEIDF(IDE,IDF)
00366       IMPLICIT NONE
00367 c should provide flavour of first incoming beam, and first tau
00368       INTEGER NMXHEP
00369       PARAMETER (NMXHEP=10000)
00370       INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
00371       REAL*8 PHEP,VHEP
00372       COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
00373      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
00374       LOGICAL QEDRAD
00375       COMMON/PH_PHOQED/QEDRAD(NMXHEP)
00376       INTEGER   IDE,IDF   
00377 
00378   
00379       IDE=IDHEP(1)
00380       IDF=IDHEP(4)
00381       IF (ABS(IDHEP(4)).EQ.ABS(IDHEP(3))) IDF=IDHEP(3)
00382       END
00383 
00384       SUBROUTINE GETBM1BM2(IDE,IDF)
00385       IMPLICIT NONE
00386 c should provide flavour of first incoming beam, and first tau
00387       INTEGER NMXHEP
00388       PARAMETER (NMXHEP=10000)
00389       INTEGER IDHEP,ISTHEP,JDAHEP,JMOHEP,NEVHEP,NHEP
00390       REAL*8 PHEP,VHEP
00391       COMMON/PH_HEPEVT/NEVHEP,NHEP,ISTHEP(NMXHEP),IDHEP(NMXHEP),
00392      &JMOHEP(2,NMXHEP),JDAHEP(2,NMXHEP),PHEP(5,NMXHEP),VHEP(4,NMXHEP)
00393       LOGICAL QEDRAD
00394       COMMON/PH_PHOQED/QEDRAD(NMXHEP)
00395       INTEGER   IDE,IDF   
00396   
00397       IDE=IDHEP(1)
00398       IDF=IDHEP(4)
00399       IF (ABS(IDHEP(4)).EQ.ABS(IDHEP(3))) IDF=IDHEP(3)
00400       END
00401       FUNCTION AFBCALC(SVAR,IDEE,IDFF)
00402 C ----------------------------------------------------------------------
00403 C THIS ROUTINE CALCULATES  BORN ASYMMETRY.
00404 C IT EXPLOITS THE FACT THAT BORN X. SECTION = A + B*C + D*C**2
00405 C
00406 C     called by : EVENTM
00407 C ----------------------------------------------------------------------
00408       IMPLICIT REAL*8 (A-H,O-Z)
00409 C
00410       CALL GIVIZO(IDEE,-1,T3e,qe,KOLOR)
00411       CALL GIVIZO(IDFF,-1,T3f,qf,KOLOR1)
00412 
00413       A=PHBORNM(svar,0.5D0,T3e,qe,T3f,qf,KOLOR*KOLOR1)
00414       B=PHBORNM(svar,-0.5D0,T3e,qe,T3f,qf,KOLOR*KOLOR1)
00415       AFBCALC= (A-B)/(A+B)*5.0/2.0 *3.0/8.0
00416       END
00417       SUBROUTINE GIVIZO(IDFERM,IHELIC,SIZO3,CHARGE,KOLOR)
00418 C ----------------------------------------------------------------------
00419 C PROVIDES ELECTRIC CHARGE AND WEAK IZOSPIN OF A FAMILY FERMION
00420 C IDFERM=1,2,3,4 DENOTES NEUTRINO, LEPTON, UP AND DOWN QUARK
00421 C NEGATIVE IDFERM=-1,-2,-3,-4, DENOTES ANTIPARTICLE
00422 C IHELIC=+1,-1 DENOTES RIGHT AND LEFT HANDEDNES ( CHIRALITY)
00423 C SIZO3 IS THIRD PROJECTION OF WEAK IZOSPIN (PLUS MINUS HALF)
00424 C AND CHARGE IS ELECTRIC CHARGE IN UNITS OF ELECTRON CHARGE
00425 C KOLOR IS A QCD COLOUR, 1 FOR LEPTON, 3 FOR QUARKS
00426 C
00427 C     called by : EVENTE, EVENTM, FUNTIH, .....
00428 C ----------------------------------------------------------------------
00429       IMPLICIT REAL*8(A-H,O-Z)
00430 C
00431       IF(IDFERM.EQ.0.OR.IABS(IDFERM).GT.4) GOTO 901
00432       IF(IABS(IHELIC).NE.1)                GOTO 901
00433       IH  =IHELIC
00434       IDTYPE =IABS(IDFERM)
00435       IC  =IDFERM/IDTYPE
00436       LEPQUA=INT(IDTYPE*0.4999999D0)
00437       IUPDOW=IDTYPE-2*LEPQUA-1
00438       CHARGE  =(-IUPDOW+2D0/3D0*LEPQUA)*IC
00439       SIZO3   =0.25D0*(IC-IH)*(1-2*IUPDOW)
00440       KOLOR=1+2*LEPQUA
00441 C** NOTE THAT CONVENTIONALY Z0 COUPLING IS
00442 C** XOUPZ=(SIZO3-CHARGE*SWSQ)/SQRT(SWSQ*(1-SWSQ))
00443       RETURN
00444  901    PRINT *,' STOP IN GIVIZO: WRONG PARAMS.'
00445       STOP
00446       END
00447  
00448 
00449       DOUBLE PRECISION FUNCTION PHBORNM(svar,costhe,T3e,qe,T3f,qf,Ncf)
00450 *///////////////////////////////////////////////////////////////////////////
00451 *//                                                                       //
00452 *// This routine provides unsophisticated Born differential cross section //
00453 *// at the crude x-section level, with Z and gamma s-chanel exchange.     //
00454 *///////////////////////////////////////////////////////////////////////////
00455       IMPLICIT NONE
00456       DOUBLE PRECISION   svar,costhe
00457       DOUBLE PRECISION   s,t,Sw2,MZ,MZ2,GammZ,MW,MW2,AlfInv,GFermi
00458       DOUBLE PRECISION   sum,T3e,t3f,qf,Qe,deno,Ve,Ae,thresh
00459       DOUBLE PRECISION   xe,yf,xf,ye,ff0,ff1,amx2,amfin,vf,af
00460       DOUBLE PRECISION   ReChiZ,SqChiZ,RaZ,RaW,ReChiW,SqChiW
00461       DOUBLE PRECISION   Born, BornS
00462       INTEGER  KeyZet,HadMin,KFbeam
00463       INTEGER  i,ke,KFfin,ncf,kf,IsGenerated,iKF
00464       INTEGER  KeyWidFix
00465       REAL*8 PI,TWOPI
00466       COMMON/PHPICO/PI,TWOPI
00467       AlfInv= 137.0359895d0
00468       GFermi=1.16639d-5
00469 
00470 *--------------------------------------------------------------------
00471       s = svar
00472 *------------------------------
00473 *     EW paratemetrs taken from BornV
00474       MZ=91.187D0
00475       GammZ=2.50072032D0
00476       Sw2=.22276773D0
00477 *------------------------------
00478 * Z and gamma couplings to beams (electrons)
00479 * Z and gamma couplings to final fermions
00480 * Loop over all flavours defined in m_xpar(400+i)
00481 
00482 
00483 *------ incoming fermion
00484             Ve=  2*T3e -4*Qe*Sw2
00485             Ae=  2*T3e
00486 *------ final fermion couplings
00487             amfin = 0.000511D0 !  m_xpar(kf+6)
00488             Vf =  2*T3f -4*Qf*Sw2
00489             Af =  2*T3f
00490             IF(abs(costhe) .GT. 1d0) THEN
00491                WRITE(*,*) '+++++STOP in PHBORN: costhe>0 =',costhe
00492                STOP
00493             ENDIF
00494             MZ2  = MZ**2
00495             RaZ  = (GFermi *MZ2 *AlfInv  )/( DSQRT(2d0) *8d0 *pi) !
00496             RaZ  = 1/(16D0*Sw2*(1d0-Sw2))
00497             KeyWidFix = 1       ! fixed width
00498             KeyWidFix = 0       ! variable width
00499             IF( KeyWidFix .EQ. 0 ) THEN
00500                ReChiZ=(s-MZ2)*s/((s-MZ2)**2+(GammZ*s/MZ)**2) *RaZ    ! variable width
00501                SqChiZ=     s**2/((s-MZ2)**2+(GammZ*s/MZ)**2) *RaZ**2 ! variable width
00502             ELSE
00503                ReChiZ=(s-MZ2)*s/((s-MZ2)**2+(GammZ*MZ)**2) *RaZ    ! fixed width
00504                SqChiZ=     s**2/((s-MZ2)**2+(GammZ*MZ)**2) *RaZ**2 ! fixed width
00505             ENDIF
00506             xe= Ve**2 +Ae**2
00507             xf= Vf**2 +Af**2
00508             ye= 2*Ve*Ae
00509             yf= 2*Vf*Af
00510             ff0= qe**2*qf**2 +2*ReChiZ*qe*qf*Ve*Vf +SqChiZ*xe*xf
00511             ff1=             +2*ReChiZ*qe*qf*Ae*Af +SqChiZ*ye*yf
00512             Born    = (1d0+ costhe**2)*ff0 +2d0*costhe*ff1
00513 * Colour factor
00514             Born = NCf*Born
00515 * Crude method of correcting threshold, cos(theta) depencence incorrect!!!
00516             IF(    svar .LE.  4d0*amfin**2) THEN
00517                thresh=0d0
00518             ELSEIF(svar .LE. 16d0*amfin**2) THEN
00519                amx2=4d0*amfin**2/svar
00520                thresh=sqrt(1d0-amx2)*(1d0+amx2/2d0)
00521             ELSE
00522                thresh=1d0
00523             ENDIF
00524             Born= Born*thresh
00525       PHBORNM = Born
00526       END
00527