ISAJET/code/qfunc.F

   1 #include "isajet/pilot.h"
   2       SUBROUTINE QFUNC
   3 C
   4 C          Find approximate QMW and QTW dependence for DRELLYAN.
   5 C          Set up /WGEN/ to generate QMW and QTW. Fit is
   6 C          Non-resonant:
   7 C            SIGMA=ANORM/(Q2/QMAX**2)**QPOW/(PT**2+RNU2)**PTPOW
   8 C          Resonant:
   9 C            SIGMA=ANORM/((Q**2-M**2)**2+M**2*GAM**2)
  10 C          with appropriate M and GAM.
  11 C
  12 C          Ver. 6.23: Remove extension of region 1 under region 2
  13 C                     to avoid discontinuity in d(sigma)/d(M)
  14 C          Ver. 6.40: Scale Q**2 fit by QMAX**2 to avoid underflow
  15 C                     problems. Must also change DRLLYN
  16 C
  17 #if defined(CERNLIB_IMPNONE)
  18       IMPLICIT NONE
  19 #endif
  20 #include "isajet/itapes.inc"
  21 #include "isajet/dypar.inc"
  22 #include "isajet/dylim.inc"
  23 #include "isajet/jetpar.inc"
  24 #include "isajet/jetlim.inc"
  25 #include "isajet/q1q2.inc"
  26 #include "isajet/wcon.inc"
  27 #include "isajet/wgen.inc"
  28 #include "isajet/jetsig.inc"
  29 #include "isajet/keys.inc"
  30 #include "isajet/hcon.inc"
  31 #include "isajet/tcpar.inc"
  32 #include "isajet/xmssm.inc"
  33 C
  34       REAL QT2CUT,DPT,QMN,QMX,EM,GAM,DELM,QSTOR,SUMS,DQ,ETAX,ETA,
  35      $Q2,XI,ALI,SIGSAV,T1,T2,T3,T4,T5,DET,DEVMAX,PTNU,ALPTNU,ALQ2,FIT,
  36      $DEV,DY3,DYW,SIG00,FACTOR,FAC1,C1,B1,SUM,AL1,QMAX2
  37       INTEGER NDIV1,NDIV2,K,I,NQS,J,N,NDIV3,NDIV4,IW,I3,II
  38       DIMENSION SUMS(9)
  39       DIMENSION QMN(3),QMX(3)
  40       DIMENSION SIGSAV(20,20)
  41 C
  42 C          QT cutoff function
  43       QT2CUT(QMW)=CUTOFF*QMW**CUTPOW
  44 C
  45 C          Entry
  46 C
  47       IF(FIXQM) THEN
  48         NDIV1=1
  49       ELSE
  50         NDIV1=20
  51       ENDIF
  52       IF(FIXQT) THEN
  53         NDIV2=1
  54       ELSE
  55         NDIV2=20
  56       ENDIF
  57 C
  58       DPT=(PTMAX(3)-PTMIN(3))/NDIV2
  59       YJ(3)=0
  60       YW=0.
  61       CTH(3)=0.
  62       STH(3)=1.
  63       IF(GODY(4)) JWTYP=4
  64       NKL=1
  65       NKH=1
  66       QMN(1)=QMIN
  67       QMX(1)=QMAX
  68       QMAX2=QMAX**2
  69 C
  70 C          Define resonance region
  71 C
  72       IF(KEYS(3)) THEN
  73         IF(JWTYP.EQ.1) GO TO 99
  74         EM=WMASS(JWTYP)
  75         GAM=WGAM(JWTYP)
  76         DELM=20.
  77       ELSEIF(KEYS(7)) THEN
  78         EM=HMASS
  79         GAM=HGAM
  80         DELM=.201357*EM
  81         DELM=AMIN1(DELM,1.5*HGAM)
  82         DELM=AMAX1(DELM,.1*EM)
  83       ELSEIF(KEYS(9)) THEN
  84         EM=TCMRHO
  85         GAM=TCGRHO
  86         DELM=.201357*EM
  87         DELM=AMIN1(DELM,1.5*TCGRHO)
  88         DELM=AMAX1(DELM,.1*EM)
  89 C          No resonance region for KKG
  90       ELSEIF(KEYS(11)) THEN
  91         EM=QMAX
  92         GAM=0.
  93         DELM=0.
  94       ENDIF
  95       EMGAM=EM*GAM
  96       EMSQ=EM**2
  97 C          Region limits
  98       QMN(2)=EM-DELM
  99       QMN(3)=EM+DELM
 100       QMX(1)=QMN(2)
 101       QMX(2)=QMN(3)
 102       NKL=1
 103       NKH=3
 104       IF(QMAX.LE.QMN(3)) NKH=2
 105       IF(QMAX.LE.QMN(2)) NKH=1
 106       IF(QMIN.GE.QMN(2)) NKL=2
 107       IF(QMIN.GE.QMN(3)) NKL=3
 108       QMX(NKH)=QMAX
 109       QMN(NKL)=QMIN
 110    99 CONTINUE
 111 C
 112 C            Fit over regions NKL to NKH
 113 C            Region 1 is below resonance
 114 C            Region 2 is inside resonance
 115 C            Region 3 is above resonance
 116 C            FIT=ANORM/(Q2/QMAX**2)**QPOW/(PT**2+RNU2)**PTPOW
 117 C
 118       DO 100 K=1,3
 119         ANORM(K)=0.
 120         PTPOW(K)=0.
 121         QPOW(K)=0.
 122         RNU2(K)=QT2CUT(QMIN)
 123 100   CONTINUE
 124 C
 125 C          Loop over regions
 126 C
 127       DO 200 K=NKL,NKH
 128         DO 210 I=1,9
 129 210     SUMS(I)=0
 130         DQ=(QMX(K)-QMN(K))/NDIV1
 131         NQS=NDIV1
 132         DO 220 I=1,NDIV2
 133           PT(3)=PTMIN(3)+(I-1)*DPT
 134           QTW=PT(3)
 135           P(3)=PT(3)
 136           RNU2(K)=QT2CUT(QMN(K))
 137           ETAX=PT(3)**2+RNU2(K)
 138           ETA=ALOG(ETAX)
 139           DO 230 J=1,NQS
 140             QMW=QMN(K)+(J-1)*DQ
 141             Q2=QMW*QMW
 142             XI=ALOG(Q2/QMAX2)
 143             SUMS(1)=SUMS(1)+XI
 144             SUMS(2)=SUMS(2)+ETA
 145             SUMS(5)=SUMS(5)+ETA*ETA
 146             SUMS(4)=SUMS(4)+XI**2
 147             SUMS(7)=SUMS(7)+XI*ETA
 148 C            Cross section
 149             IF(KEYS(3)) THEN
 150               CALL SIGDY
 151             ELSEIF(KEYS(7).AND..NOT.GOMSSM) THEN
 152               CALL SIGH
 153             ELSEIF(KEYS(7).AND.GOMSSM) THEN
 154               CALL SIGHSS
 155             ELSEIF(KEYS(9)) THEN
 156               CALL SIGTC
 157             ELSEIF(KEYS(11)) THEN
 158               CALL SIGKKG
 159             ENDIF
 160             IF(SIGMA.EQ.0.) GO TO 999
 161             AL1=ALOG(SIGMA)
 162             SIGSAV(I,J)=AL1
 163             IF(K.EQ.2) AL1=AL1+ALOG((Q2-EM**2)**2+EMGAM**2)
 164             SUMS(3)=SUMS(3)+AL1
 165             SUMS(8)=SUMS(8)+AL1*XI
 166             SUMS(9)=SUMS(9)+AL1*ETA
 167 230         CONTINUE
 168 220      CONTINUE
 169 C
 170 C            Find coefficients minimizing chisq
 171 C
 172         N=NQS*NDIV2
 173         T1=N*SUMS(7)-SUMS(1)*SUMS(2)
 174         T2=N*SUMS(5)-SUMS(2)**2
 175         T3=N*SUMS(4)-SUMS(1)**2
 176         T4=N*SUMS(8)-SUMS(1)*SUMS(3)
 177         T5=N*SUMS(9)-SUMS(2)*SUMS(3)
 178         IF((FIXQM.OR.K.EQ.2).AND.FIXQT) THEN
 179           PTPOW(K)=0.
 180           QPOW(K)=0.
 181         ELSEIF(FIXQT) THEN
 182           PTPOW(K)=0.
 183           QPOW(K)=-T4/T3
 184         ELSEIF(FIXQM.OR.K.EQ.2) THEN
 185           PTPOW(K)=-T5/T2
 186           QPOW(K)=0.
 187         ELSE
 188           DET=T1**2-T2*T3
 189           PTPOW(K)=(T5*T3-T4*T1)/DET
 190           QPOW(K)=(T4*T2-T1*T5)/DET
 191         ENDIF
 192         ANORM(K)=(QPOW(K)*SUMS(1)+PTPOW(K)*SUMS(2)+SUMS(3))/N
 193 C
 194 C          Shift fit to obtain envelope for SIGDY
 195 C
 196         DEVMAX=0.
 197         DO 240 I=1,NDIV2
 198           PT(3)=PTMIN(3)+(I-1)*DPT
 199           PTNU=PT(3)**2+RNU2(K)
 200           DO 250 J=1,NDIV1
 201             QMW=QMN(K)+(J-1)*DQ
 202             Q2=QMW**2
 203             ALPTNU=ALOG(PTNU)
 204             ALQ2=ALOG(Q2/QMAX2)
 205             IF(K.EQ.2) THEN
 206               FIT=EXP(ANORM(K)-PTPOW(K)*ALPTNU
 207      $        -ALOG((Q2-EM**2)**2+EMGAM**2))
 208             ELSE
 209               FIT=EXP(ANORM(K)-PTPOW(K)*ALPTNU-QPOW(K)*ALQ2)
 210             ENDIF
 211             DEV=SIGSAV(I,J)-ALOG(FIT)
 212             IF(DEV.GT.DEVMAX) DEVMAX=DEV
 213 250       CONTINUE
 214 240     CONTINUE
 215         ANORM(K)=ANORM(K)+DEVMAX
 216 200   CONTINUE
 217 C
 218 C          Shift fit to obtain envelope in YW
 219       NDIV3=20
 220       IF(STDDY) THEN
 221         NDIV4=1
 222         DY3=0.
 223       ELSE
 224         NDIV4=20
 225         DY3=(YJMAX(3)-YJMIN(3))/(NDIV4-1)
 226       ENDIF
 227       DYW=(YWMAX-YWMIN)/(NDIV3-1)
 228 C
 229       DO 300 K=NKL,NKH
 230         QMW=QMN(K)
 231         Q2=QMW**2
 232         QTW=QTMIN
 233         PT(3)=QTW
 234         P(3)=PT(3)
 235         YW=0.
 236         YJ(3)=0.
 237         CTH(3)=0.
 238         STH(3)=1.
 239         IF(KEYS(3)) THEN
 240           CALL SIGDY
 241         ELSEIF(KEYS(7).AND..NOT.GOMSSM) THEN
 242           CALL SIGH
 243         ELSEIF(KEYS(7).AND.GOMSSM) THEN
 244           CALL SIGHSS
 245         ELSEIF(KEYS(9)) THEN
 246           CALL SIGTC
 247         ELSEIF(KEYS(11)) THEN
 248           CALL SIGKKG
 249         ENDIF
 250         SIG00=SIGMA
 251         FACTOR=1.
 252         DO 310 IW=1,NDIV3
 253           YW=YWMIN+(IW-1)*DYW
 254           DO 320 I3=1,NDIV4
 255             IF(.NOT.STDDY) THEN
 256               YJ(3)=YJMIN(3)+(I3-1)*DY3
 257               CTH(3)=TANH(YJ(3))
 258               STH(3)=SQRT(1.-CTH(3)**2)
 259               IF(STH(3).EQ.0.) GO TO 320
 260               TH(3)=ACOS(CTH(3))
 261               P(3)=PT(3)/STH(3)
 262             ENDIF
 263             IF(KEYS(3)) THEN
 264               CALL SIGDY
 265             ELSEIF(KEYS(7).AND..NOT.GOMSSM) THEN
 266               CALL SIGH
 267             ELSEIF(KEYS(7).AND.GOMSSM) THEN
 268               CALL SIGHSS
 269             ELSEIF(KEYS(9)) THEN
 270               CALL SIGTC
 271             ELSEIF(KEYS(11)) THEN
 272               CALL SIGKKG
 273             ENDIF
 274             FAC1=SIGMA/SIG00
 275             FACTOR=AMAX1(FACTOR,FAC1)
 276 320       CONTINUE
 277 310     CONTINUE
 278         ANORM(K)=ALOG(FACTOR)+ANORM(K)
 279 300   CONTINUE
 280 C
 281 C          Set up generating constants for PT**2 and QMW**2
 282 C
 283       DO 400 K=NKL,NKH
 284         C1=1.-PTPOW(K)
 285         PTGN(1,K)=(PTMIN(3)**2+RNU2(K))**C1
 286         PTGN(2,K)=(PTMAX(3)**2+RNU2(K))**C1-PTGN(1,K)
 287         PTGN(3,K)=1./C1
 288         IF(K.EQ.2) THEN
 289           QGEN(1,2)=ATAN((QMN(2)**2-EMSQ)/EMGAM)
 290           QGEN(2,2)=ATAN((QMX(2)**2-EMSQ)/EMGAM)-QGEN(1,2)
 291           QGEN(3,2)=EMGAM
 292         ELSE
 293           B1=1.-QPOW(K)
 294           QGEN(1,K)=(QMN(K)/QMAX)**(2.*B1)
 295           QGEN(2,K)=(QMX(K)/QMAX)**(2.*B1)-QGEN(1,K)
 296           QGEN(3,K)=1./B1
 297         ENDIF
 298 400   CONTINUE
 299 C
 300       DO 410 K=1,3
 301 410   QSELWT(K)=0.
 302       SUM=0.
 303 C
 304       DO 420 K=NKL,NKH
 305         QSELWT(K)=1.
 306         IF(.NOT.FIXQT) QSELWT(K)=QSELWT(K)*PTGN(2,K)*PTGN(3,K)
 307         IF(.NOT.FIXQM) THEN
 308           IF(K.EQ.2) THEN
 309             QSELWT(K)=QSELWT(K)*QGEN(2,K)/EMGAM
 310           ELSE
 311             QSELWT(K)=QMAX**2*QSELWT(K)*QGEN(2,K)*QGEN(3,K)
 312           ENDIF
 313         ENDIF
 314         QSELWT(K)=EXP(ALOG(QSELWT(K))+ANORM(K))
 315         SUM=SUM+QSELWT(K)
 316 420   CONTINUE
 317 C
 318       DO 430 K=1,3
 319         QSELWT(K)=QSELWT(K)/SUM
 320 430   CONTINUE
 321 C
 322 C          Write fit to output
 323 C
 324       WRITE(ITLIS,4301)
 325 4301  FORMAT(//10X,' QT AND Q FIRST GENERATED BY--'/)
 326       DO 440 K=NKL,NKH
 327         WRITE(ITLIS,4402) K,QMN(K),QMX(K)
 328 4402    FORMAT(//5X,' REGION',I2,5X,E11.4,' < Q < ',E11.5)
 329         WRITE(ITLIS,4403) (PTGN(II,K),II=1,3),RNU2(K)
 330 4403    FORMAT(/' QT**2 = (',E11.4,' + ',E11.4,' * RANF) ** ',E11.4,
 331      $  ' - ',E11.4)
 332         IF(K.NE.2) THEN
 333           WRITE(ITLIS,4404) QMAX2,(QGEN(II,K),II=1,3)
 334 4404      FORMAT(/' Q**2  = ',E11.4,' * (',E11.4,' + ',E11.4,
 335      $    ' * RANF) ** ',E11.4)
 336         ELSE
 337           WRITE(ITLIS,4505) QGEN(3,K),QGEN(1,K),QGEN(2,K),EMSQ
 338 4505       FORMAT(/' Q**2  = ',E11.4,' * TAN(',E11.4,' + ',E11.4,
 339      $    ' * RANF) + ',E11.4)
 340         ENDIF
 341         WRITE(ITLIS,4506) QSELWT(K)
 342 4506    FORMAT(/' WEIGHT = ',E11.4)
 343 440   CONTINUE
 344 C
 345 C          Set fixed limits if any
 346 C
 347       IF(FIXQT) THEN
 348         PTMAX(3)=PTMIN(3)
 349         PT(3)=PTMIN(3)
 350         QTW=PT(3)
 351       ENDIF
 352       IF(FIXQM) THEN
 353         QMAX=QMIN
 354         QMW=QMIN
 355       ENDIF
 356       RETURN
 357 C
 358 C          Fit fails if SIGMA=0 in allowed range
 359 C
 360 999   WRITE(ITLIS,9990) QMW,QTW
 361 9990  FORMAT(//' ERROR IN QFUNC...SIGMA=0 FOR QMW = ',E12.4,' , QTW = ',
 362      1E12.4/' CHECK YOUR LIMITS')
 363       STOP 99
 364       END