Merging the VirtualMC branch to the main development branch (HEAD)

[u/mrichter/AliRoot.git] / ISAJET / code / sigh3.F

#include "isajet/pilot.h"
      SUBROUTINE SIGH3
C
C          Calculate angular distributions for W decays from Higgs,
C          d(sigma)/d(qmw**2)d(yw)d(omega)d(omega1)d(omega2)
C
C          Ver 7.14: Only modification needed for MSSM is to check
C          GOMSSM flag instead of INITYP
C
#include "isajet/itapes.inc"
#include "isajet/qcdpar.inc"
#include "isajet/jetpar.inc"
#include "isajet/pjets.inc"
#include "isajet/primar.inc"
#include "isajet/q1q2.inc"
#include "isajet/jetsig.inc"
#include "isajet/wsig.inc"
#include "isajet/wwsig.inc"
#include "isajet/qsave.inc"
#include "isajet/wcon.inc"
#include "isajet/const.inc"
#include "isajet/wwpar.inc"
#include "isajet/hcon.inc"
#include "isajet/xmssm.inc"
C
      EQUIVALENCE (S,SHAT),(T,THAT),(U,UHAT)
      DIMENSION IDADDR(4),IW(2),LAM(3),LISTJ(29)
     $,T12(3,3),T34(3,3),FTERM(4),FR(3,3),FI(3,3)
     $,CPHI12(3),SPHI12(3),CPHI34(3),SPHI34(3)
      DIMENSION PFCM(5,4),PWCM(5,2)
#if defined(CERNLIB_DOUBLE)
      DOUBLE PRECISION TERM,FTERM,ZCM
#endif
      DATA LAM/0,1,-1/
      DATA LISTJ/
     $9,1,-1,2,-2,3,-3,4,-4,5,-5,6,-6,
     $11,-11,12,-12,13,-13,14,-14,15,-15,16,-16,
     $10,80,-80,90/
C
C          FUNCTIONS
      DOTP(I,J)=PPAIR(4,I)*PPAIR(4,J)-PPAIR(1,I)*PPAIR(1,J)
     $-PPAIR(2,I)*PPAIR(2,J)-PPAIR(3,I)*PPAIR(3,J)
C
C          ENTRY
      IF(NPAIR.NE.4) RETURN
C
C          RECONSTRUCT W-->FF DECAY ANGLES
C
C          INITIALIZE PFCM AND PWCM
      DO 10 I=1,4
      DO 10 K=1,5
10    PFCM(K,I)=PPAIR(K,I)
      DO 11 I=1,2
      DO 11 K=1,5
11    PWCM(K,I)=PJETS(K,I)
C
C          Z BOOST TO WW CENTER OF MASS
      CHWW=QWJET(4)/QWJET(5)
      SHWW=QWJET(3)/QWJET(5)
      DO 20 I=1,4
      TMP=CHWW*PFCM(4,I)-SHWW*PFCM(3,I)
      PFCM(3,I)=-SHWW*PFCM(4,I)+CHWW*PFCM(3,I)
20    PFCM(4,I)=TMP
      DO 21 I=1,2
      TMP=CHWW*PWCM(4,I)-SHWW*PWCM(3,I)
      PWCM(3,I)=-SHWW*PWCM(4,I)+CHWW*PWCM(3,I)
21    PWCM(4,I)=TMP
C
C          ROTATE W1 TO +Z AXIS
      PTW1=SQRT(PWCM(1,1)**2+PWCM(2,1)**2)
      CPHIW1=PWCM(1,1)/PTW1
      SPHIW1=PWCM(2,1)/PTW1
      PW1=SQRT(PTW1**2+PWCM(3,1)**2)
      CTHW1=PWCM(3,1)/PW1
      STHW1=PTW1/PW1
C          Z ROTATION
      DO 30 I=1,4
      TMP=CPHIW1*PFCM(1,I)+SPHIW1*PFCM(2,I)
      PFCM(2,I)=-SPHIW1*PFCM(1,I)+CPHIW1*PFCM(2,I)
30    PFCM(1,I)=TMP
C          Y ROTATION
      DO 31 I=1,4
      TMP=CTHW1*PFCM(1,I)-STHW1*PFCM(3,I)
      PFCM(3,I)=STHW1*PFCM(1,I)+CTHW1*PFCM(3,I)
31    PFCM(1,I)=TMP
C
C          BOOST TO W REST FRAMES
      CHW1=PWCM(4,1)/PWCM(5,1)
      SHW1=PW1/PWCM(5,1)
      DO 40 I=1,4
      IF(I.LE.2) THEN
        SHWI=SHW1
      ELSE
        SHWI=-SHW1
      ENDIF
      TMP=CHW1*PFCM(4,I)-SHWI*PFCM(3,I)
      PFCM(3,I)=-SHWI*PFCM(4,I)+CHW1*PFCM(3,I)
40    PFCM(4,I)=TMP
C
C          COMPUTE ANGLES
      TH12=ACOS(PFCM(3,1)/SQRT(PFCM(1,1)**2+PFCM(2,1)**2+PFCM(3,1)**2))
      PHI12=ATAN2(PFCM(2,1),PFCM(1,1))
      TH34=ACOS(PFCM(3,3)/SQRT(PFCM(1,3)**2+PFCM(2,3)**2+PFCM(3,3)**2))
      PHI34=ATAN2(PFCM(2,3),PFCM(1,3))
C
C          COMPUTE DECAY ANGULAR DISTRIBUTIONS.
C
      DO 100 I=1,4
      IDADDR(I)=IABS(IDPAIR(I))
100   IF(IDADDR(I).GE.11) IDADDR(I)=IDADDR(I)-4
      IF(GOMSSM) THEN
        IW(1)=JETTYP(1)-76
        IW(2)=JETTYP(2)-76
      ELSE
        IW(1)=JETTYP(1)-25
        IW(2)=JETTYP(2)-25
      ENDIF
C
      AMV=PJETS(5,1)
      GAMV=WGAM(IW(1))
      QMH=QMW
C          COUPLINGS
      A12=AQ(IDADDR(1),IW(1))
      B12=BQ(IDADDR(1),IW(1))
      A34=AQ(IDADDR(3),IW(2))
      B34=BQ(IDADDR(3),IW(2))
C          DECAY DISTRIBUTIONS
      TVV12=8.*PI*ALFA*(A12**2+B12**2)
      TVA12=16.*PI*ALFA*A12*B12
      COS12=COS(TH12)
      SIN12=SIN(TH12)
      T12(1,1)=TVV12*SIN12**2
      T12(1,2)=TVV12*SIN12*COS12/SQRT2+TVA12*SIN12/SQRT2
      T12(1,3)=-TVV12*SIN12*COS12/SQRT2+TVA12*SIN12/SQRT2
      T12(2,1)=T12(1,2)
      T12(2,2)=TVV12*(.5+.5*COS12**2)+TVA12*COS12
      T12(2,3)=TVV12*.5*SIN12**2
      T12(3,1)=T12(1,3)
      T12(3,2)=T12(2,3)
      T12(3,3)=TVV12*(.5+.5*COS12**2)-TVA12*COS12
C
      TVV34=8.*PI*ALFA*(A34**2+B34**2)
      TVA34=16.*PI*ALFA*A34*B34
      COS34=COS(TH34)
      SIN34=SIN(TH34)
      T34(1,1)=TVV34*SIN34**2
      T34(1,2)=TVV34*SIN34*COS34/SQRT2+TVA34*SIN34/SQRT2
      T34(1,3)=-TVV34*SIN34*COS34/SQRT2+TVA34*SIN34/SQRT2
      T34(2,1)=T34(1,2)
      T34(2,2)=TVV34*(.5+.5*COS34**2)+TVA34*COS34
      T34(2,3)=TVV34*.5*SIN34**2
      T34(3,1)=T34(1,3)
      T34(3,2)=T34(2,3)
      T34(3,3)=TVV34*(.5+.5*COS34**2)-TVA34*COS34
C
      CPHI12(1)=1.
      CPHI12(2)=COS(PHI12)
      CPHI12(3)=COS(2.*PHI12)
      SPHI12(1)=0.
      SPHI12(2)=SIN(PHI12)
      SPHI12(3)=SIN(2.*PHI12)
      CPHI34(1)=1.
      CPHI34(2)=COS(PHI34)
      CPHI34(3)=COS(2.*PHI34)
      SPHI34(1)=0.
      SPHI34(2)=SIN(PHI34)
      SPHI34(3)=SIN(2.*PHI34)
C
      TCPHI=CPHI12(2)*CPHI34(2)-SPHI12(2)*SPHI34(2)
      TSPHI=SPHI12(2)*CPHI34(2)+CPHI12(2)*SPHI34(2)
      TC2PHI=CPHI12(3)*CPHI34(3)-SPHI12(3)*SPHI34(3)
      TS2PHI=SPHI12(3)*CPHI34(3)+CPHI12(3)*SPHI34(3)
C
C          PURE HIGGS --> W W. CALCULATE ANGULAR DISTRIBUTION FOR
C          HIGGS DECAY AND MULTIPLY BY CROSS SECTION.
C
      IF(INITYP(1).LE.25.OR.GOMSSM) THEN
        F0=.5*QMH**2/AMV**2-1.
        F1=1.
        TOTAL=(8.*PI/3.)**2*TVV12*TVV34*(F0**2+2.*F1**2)
        DIFF=F0**2*T12(1,1)*T34(1,1)
     $  +F0*F1*(2.*T12(1,2)*T34(1,2)+2.*T12(1,3)*T34(1,3))*TCPHI
     $  +F1**2*(T12(2,2)*T34(1,2)+T12(3,3)*T34(3,3)
     $    +2.*T12(2,3)*T34(2,3)*TC2PHI)
        WWSIG=SIGLLQ*DIFF/TOTAL
        RETURN
      ENDIF
C
C          W W FUSION. CALCULATE ANGULAR DISTRIUBTION FOR DECAY
C          INCLUDING ALL GRAPHS.
C
C          KINEMATICS
      IFL1=LISTJ(JETTYP(1))
      IFL2=LISTJ(JETTYP(2))
      IFIN1=LISTJ(INITYP(1))
      IFIN2=LISTJ(INITYP(2))
      WMF=AMASS(IFL1)
      WMI=AMASS(IFIN1)
      PINPF=SQRT((S-4.*WMI**2)*(S-4.*WMF**2))
      ZCM=(.5*S+T-WMI**2-WMF**2)/(.5*PINPF)
C          PRODUCTION AMPLITUDES. REMEMBER MISSING SIN(THETA)/SQRT(2)
      DO 110 L=1,4
      FTERM(L)=0.
      DO 120 J=1,4
      TERM=0.
      DO 130 I=1,4
130   TERM=TERM+ANWWWW(I,J,L)*ZCM**(I-1)
      TERM=TERM/(ADWWWW(1,J)+ADWWWW(2,J)*ZCM)
120   FTERM(L)=FTERM(L)+TERM
110   CONTINUE
      FTERM(4)=FTERM(4)*SQRT(ABS(1.-ZCM**2))/SQRT2
C          HELICITY AMPLITUDES. NOTATION IS 0,+,-
      FR(1,1)=FTERM(1)
      FI(1,1)=AIWWWW(1)
      FR(1,2)=FTERM(4)
      FI(1,2)=AIWWWW(4)
      FR(2,2)=FTERM(3)
      FI(2,2)=AIWWWW(3)
      FR(2,3)=FTERM(2)
      FI(2,3)=AIWWWW(2)
C
      FR(1,3)=FR(1,2)
      FI(1,3)=FI(1,2)
      FR(3,1)=FR(1,3)
      FI(3,1)=FI(1,3)
      FR(2,1)=FR(1,2)
      FI(2,1)=FI(1,2)
C
      FR(3,3)=FR(2,2)
      FI(3,3)=FI(2,2)
      FR(3,2)=FR(2,3)
      FI(3,2)=FI(2,3)
C
C          DIFFERENTIAL DISTRIBUTION FROM DENSITY MATRIX
      DIFF=0.
      DO 140 I1=1,3
      L1=LAM(I1)
      DO 140 I2=1,3
      L2=LAM(I2)
      DO 140 I3=1,3
      L3=LAM(I3)
      DO 140 I4=1,3
      L4=LAM(I4)
      L12=L1-L2
      I12=IABS(L12)+1
      IF(I12.EQ.0) I12=3
      L34=L3-L4
      I34=IABS(L34)+1
      IF(I34.EQ.0) I34=3
      C1234=CPHI12(I12)*CPHI34(I34)
     $-SPHI12(I12)*ISIGN(1,L12)*SPHI34(I34)*ISIGN(1,L34)
      S1234=SPHI12(I12)*ISIGN(1,L12)*CPHI34(I34)
     $+CPHI12(I12)*SPHI34(I34)*ISIGN(1,L34)
      DIFF=DIFF+(FR(I1,I2)*FR(I3,I4)+FI(I1,I2)*FI(I3,I4))
     $*T12(I3,I1)*T34(I4,I2)*C1234
     $+(FR(I1,I2)*FI(I3,I4)-FI(I1,I2)*FR(I3,I4))
     $*T12(I3,I1)*T34(I4,I2)*S1234
140   CONTINUE
C          INTEGRATED DISTRIBUTION
      TOTAL=0.
      DO 150 I1=1,3
      DO 150 I2=1,3
      TOTAL=TOTAL+FR(I1,I2)**2+FI(I1,I2)**2
150   CONTINUE
      FAC=(16.*PI/3.*4.*PI*ALFA)**2
      FAC=FAC*(A12**2+B12**2)*(A34**2+B34**2)
      TOTAL=TOTAL*FAC
      WWSIG=DIFF/TOTAL*SIGLLQ
      RETURN
      END