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

#include "isajet/pilot.h"
      SUBROUTINE QFUNC
C
C          Find approximate QMW and QTW dependence for DRELLYAN.
C          Set up /WGEN/ to generate QMW and QTW. Fit is
C          Non-resonant:
C            SIGMA=ANORM/(Q2/QMAX**2)**QPOW/(PT**2+RNU2)**PTPOW
C          Resonant:
C            SIGMA=ANORM/((Q**2-M**2)**2+M**2*GAM**2)
C          with appropriate M and GAM.
C
C          Ver. 6.23: Remove extension of region 1 under region 2
C                     to avoid discontinuity in d(sigma)/d(M)
C          Ver. 6.40: Scale Q**2 fit by QMAX**2 to avoid underflow
C                     problems. Must also change DRLLYN
C
#if defined(CERNLIB_IMPNONE)
      IMPLICIT NONE
#endif
#include "isajet/itapes.inc"
#include "isajet/dypar.inc"
#include "isajet/dylim.inc"
#include "isajet/jetpar.inc"
#include "isajet/jetlim.inc"
#include "isajet/q1q2.inc"
#include "isajet/wcon.inc"
#include "isajet/wgen.inc"
#include "isajet/jetsig.inc"
#include "isajet/keys.inc"
#include "isajet/hcon.inc"
#include "isajet/tcpar.inc"
#include "isajet/xmssm.inc"
C
      REAL QT2CUT,DPT,QMN,QMX,EM,GAM,DELM,QSTOR,SUMS,DQ,ETAX,ETA,
     $Q2,XI,ALI,SIGSAV,T1,T2,T3,T4,T5,DET,DEVMAX,PTNU,ALPTNU,ALQ2,FIT,
     $DEV,DY3,DYW,SIG00,FACTOR,FAC1,C1,B1,SUM,AL1,QMAX2
      INTEGER NDIV1,NDIV2,K,I,NQS,J,N,NDIV3,NDIV4,IW,I3,II
      DIMENSION SUMS(9)
      DIMENSION QMN(3),QMX(3)
      DIMENSION SIGSAV(20,20)
C
C          QT cutoff function
      QT2CUT(QMW)=CUTOFF*QMW**CUTPOW
C
C          Entry
C
      IF(FIXQM) THEN
        NDIV1=1
      ELSE
        NDIV1=20
      ENDIF
      IF(FIXQT) THEN
        NDIV2=1
      ELSE
        NDIV2=20
      ENDIF
C
      DPT=(PTMAX(3)-PTMIN(3))/NDIV2
      YJ(3)=0
      YW=0.
      CTH(3)=0.
      STH(3)=1.
      IF(GODY(4)) JWTYP=4
      NKL=1
      NKH=1
      QMN(1)=QMIN
      QMX(1)=QMAX
      QMAX2=QMAX**2
C
C          Define resonance region
C
      IF(KEYS(3)) THEN
        IF(JWTYP.EQ.1) GO TO 99
        EM=WMASS(JWTYP)
        GAM=WGAM(JWTYP)
        DELM=20.
      ELSEIF(KEYS(7)) THEN
        EM=HMASS
        GAM=HGAM
        DELM=.201357*EM
        DELM=AMIN1(DELM,1.5*HGAM)
        DELM=AMAX1(DELM,.1*EM)
      ELSEIF(KEYS(9)) THEN
        EM=TCMRHO
        GAM=TCGRHO
        DELM=.201357*EM
        DELM=AMIN1(DELM,1.5*TCGRHO)
        DELM=AMAX1(DELM,.1*EM)
C          No resonance region for KKG
      ELSEIF(KEYS(11)) THEN
        EM=QMAX
        GAM=0.
        DELM=0.
      ENDIF
      EMGAM=EM*GAM
      EMSQ=EM**2
C          Region limits
      QMN(2)=EM-DELM
      QMN(3)=EM+DELM
      QMX(1)=QMN(2)
      QMX(2)=QMN(3)
      NKL=1
      NKH=3
      IF(QMAX.LE.QMN(3)) NKH=2
      IF(QMAX.LE.QMN(2)) NKH=1
      IF(QMIN.GE.QMN(2)) NKL=2
      IF(QMIN.GE.QMN(3)) NKL=3
      QMX(NKH)=QMAX
      QMN(NKL)=QMIN
   99 CONTINUE
C
C            Fit over regions NKL to NKH
C            Region 1 is below resonance
C            Region 2 is inside resonance
C            Region 3 is above resonance
C            FIT=ANORM/(Q2/QMAX**2)**QPOW/(PT**2+RNU2)**PTPOW
C
      DO 100 K=1,3
        ANORM(K)=0.
        PTPOW(K)=0.
        QPOW(K)=0.
        RNU2(K)=QT2CUT(QMIN)
100   CONTINUE
C
C          Loop over regions
C
      DO 200 K=NKL,NKH
        DO 210 I=1,9
210     SUMS(I)=0
        DQ=(QMX(K)-QMN(K))/NDIV1
        NQS=NDIV1
        DO 220 I=1,NDIV2
          PT(3)=PTMIN(3)+(I-1)*DPT
          QTW=PT(3)
          P(3)=PT(3)
          RNU2(K)=QT2CUT(QMN(K))
          ETAX=PT(3)**2+RNU2(K)
          ETA=ALOG(ETAX)
          DO 230 J=1,NQS
            QMW=QMN(K)+(J-1)*DQ
            Q2=QMW*QMW
            XI=ALOG(Q2/QMAX2)
            SUMS(1)=SUMS(1)+XI
            SUMS(2)=SUMS(2)+ETA
            SUMS(5)=SUMS(5)+ETA*ETA
            SUMS(4)=SUMS(4)+XI**2
            SUMS(7)=SUMS(7)+XI*ETA
C            Cross section
            IF(KEYS(3)) THEN
              CALL SIGDY
            ELSEIF(KEYS(7).AND..NOT.GOMSSM) THEN
              CALL SIGH
            ELSEIF(KEYS(7).AND.GOMSSM) THEN
              CALL SIGHSS
            ELSEIF(KEYS(9)) THEN
              CALL SIGTC
            ELSEIF(KEYS(11)) THEN
              CALL SIGKKG
            ENDIF
            IF(SIGMA.EQ.0.) GO TO 999
            AL1=ALOG(SIGMA)
            SIGSAV(I,J)=AL1
            IF(K.EQ.2) AL1=AL1+ALOG((Q2-EM**2)**2+EMGAM**2)
            SUMS(3)=SUMS(3)+AL1
            SUMS(8)=SUMS(8)+AL1*XI
            SUMS(9)=SUMS(9)+AL1*ETA
230         CONTINUE
220      CONTINUE
C
C            Find coefficients minimizing chisq
C
        N=NQS*NDIV2
        T1=N*SUMS(7)-SUMS(1)*SUMS(2)
        T2=N*SUMS(5)-SUMS(2)**2
        T3=N*SUMS(4)-SUMS(1)**2
        T4=N*SUMS(8)-SUMS(1)*SUMS(3)
        T5=N*SUMS(9)-SUMS(2)*SUMS(3)
        IF((FIXQM.OR.K.EQ.2).AND.FIXQT) THEN
          PTPOW(K)=0.
          QPOW(K)=0.
        ELSEIF(FIXQT) THEN
          PTPOW(K)=0.
          QPOW(K)=-T4/T3
        ELSEIF(FIXQM.OR.K.EQ.2) THEN
          PTPOW(K)=-T5/T2
          QPOW(K)=0.
        ELSE
          DET=T1**2-T2*T3
          PTPOW(K)=(T5*T3-T4*T1)/DET
          QPOW(K)=(T4*T2-T1*T5)/DET
        ENDIF
        ANORM(K)=(QPOW(K)*SUMS(1)+PTPOW(K)*SUMS(2)+SUMS(3))/N
C
C          Shift fit to obtain envelope for SIGDY
C
        DEVMAX=0.
        DO 240 I=1,NDIV2
          PT(3)=PTMIN(3)+(I-1)*DPT
          PTNU=PT(3)**2+RNU2(K)
          DO 250 J=1,NDIV1
            QMW=QMN(K)+(J-1)*DQ
            Q2=QMW**2
            ALPTNU=ALOG(PTNU)
            ALQ2=ALOG(Q2/QMAX2)
            IF(K.EQ.2) THEN
              FIT=EXP(ANORM(K)-PTPOW(K)*ALPTNU
     $        -ALOG((Q2-EM**2)**2+EMGAM**2))
            ELSE
              FIT=EXP(ANORM(K)-PTPOW(K)*ALPTNU-QPOW(K)*ALQ2)
            ENDIF
            DEV=SIGSAV(I,J)-ALOG(FIT)
            IF(DEV.GT.DEVMAX) DEVMAX=DEV
250       CONTINUE
240     CONTINUE
        ANORM(K)=ANORM(K)+DEVMAX
200   CONTINUE
C
C          Shift fit to obtain envelope in YW
      NDIV3=20
      IF(STDDY) THEN
        NDIV4=1
        DY3=0.
      ELSE
        NDIV4=20
        DY3=(YJMAX(3)-YJMIN(3))/(NDIV4-1)
      ENDIF
      DYW=(YWMAX-YWMIN)/(NDIV3-1)
C
      DO 300 K=NKL,NKH
        QMW=QMN(K)
        Q2=QMW**2
        QTW=QTMIN
        PT(3)=QTW
        P(3)=PT(3)
        YW=0.
        YJ(3)=0.
        CTH(3)=0.
        STH(3)=1.
        IF(KEYS(3)) THEN
          CALL SIGDY
        ELSEIF(KEYS(7).AND..NOT.GOMSSM) THEN
          CALL SIGH
        ELSEIF(KEYS(7).AND.GOMSSM) THEN
          CALL SIGHSS
        ELSEIF(KEYS(9)) THEN
          CALL SIGTC
        ELSEIF(KEYS(11)) THEN
          CALL SIGKKG
        ENDIF
        SIG00=SIGMA
        FACTOR=1.
        DO 310 IW=1,NDIV3
          YW=YWMIN+(IW-1)*DYW
          DO 320 I3=1,NDIV4
            IF(.NOT.STDDY) THEN
              YJ(3)=YJMIN(3)+(I3-1)*DY3
              CTH(3)=TANH(YJ(3))
              STH(3)=SQRT(1.-CTH(3)**2)
              IF(STH(3).EQ.0.) GO TO 320
              TH(3)=ACOS(CTH(3))
              P(3)=PT(3)/STH(3)
            ENDIF
            IF(KEYS(3)) THEN
              CALL SIGDY
            ELSEIF(KEYS(7).AND..NOT.GOMSSM) THEN
              CALL SIGH
            ELSEIF(KEYS(7).AND.GOMSSM) THEN
              CALL SIGHSS
            ELSEIF(KEYS(9)) THEN
              CALL SIGTC
            ELSEIF(KEYS(11)) THEN
              CALL SIGKKG
            ENDIF
            FAC1=SIGMA/SIG00
            FACTOR=AMAX1(FACTOR,FAC1)
320       CONTINUE
310     CONTINUE
        ANORM(K)=ALOG(FACTOR)+ANORM(K)
300   CONTINUE
C
C          Set up generating constants for PT**2 and QMW**2
C
      DO 400 K=NKL,NKH
        C1=1.-PTPOW(K)
        PTGN(1,K)=(PTMIN(3)**2+RNU2(K))**C1
        PTGN(2,K)=(PTMAX(3)**2+RNU2(K))**C1-PTGN(1,K)
        PTGN(3,K)=1./C1
        IF(K.EQ.2) THEN
          QGEN(1,2)=ATAN((QMN(2)**2-EMSQ)/EMGAM)
          QGEN(2,2)=ATAN((QMX(2)**2-EMSQ)/EMGAM)-QGEN(1,2)
          QGEN(3,2)=EMGAM
        ELSE
          B1=1.-QPOW(K)
          QGEN(1,K)=(QMN(K)/QMAX)**(2.*B1)
          QGEN(2,K)=(QMX(K)/QMAX)**(2.*B1)-QGEN(1,K)
          QGEN(3,K)=1./B1
        ENDIF
400   CONTINUE
C
      DO 410 K=1,3
410   QSELWT(K)=0.
      SUM=0.
C
      DO 420 K=NKL,NKH
        QSELWT(K)=1.
        IF(.NOT.FIXQT) QSELWT(K)=QSELWT(K)*PTGN(2,K)*PTGN(3,K)
        IF(.NOT.FIXQM) THEN
          IF(K.EQ.2) THEN
            QSELWT(K)=QSELWT(K)*QGEN(2,K)/EMGAM
          ELSE
            QSELWT(K)=QMAX**2*QSELWT(K)*QGEN(2,K)*QGEN(3,K)
          ENDIF
        ENDIF
        QSELWT(K)=EXP(ALOG(QSELWT(K))+ANORM(K))
        SUM=SUM+QSELWT(K)
420   CONTINUE
C
      DO 430 K=1,3
        QSELWT(K)=QSELWT(K)/SUM
430   CONTINUE
C
C          Write fit to output
C
      WRITE(ITLIS,4301)
4301  FORMAT(//10X,' QT AND Q FIRST GENERATED BY--'/)
      DO 440 K=NKL,NKH
        WRITE(ITLIS,4402) K,QMN(K),QMX(K)
4402    FORMAT(//5X,' REGION',I2,5X,E11.4,' < Q < ',E11.5)
        WRITE(ITLIS,4403) (PTGN(II,K),II=1,3),RNU2(K)
4403    FORMAT(/' QT**2 = (',E11.4,' + ',E11.4,' * RANF) ** ',E11.4,
     $  ' - ',E11.4)
        IF(K.NE.2) THEN
          WRITE(ITLIS,4404) QMAX2,(QGEN(II,K),II=1,3)
4404      FORMAT(/' Q**2  = ',E11.4,' * (',E11.4,' + ',E11.4,
     $    ' * RANF) ** ',E11.4)
        ELSE
          WRITE(ITLIS,4505) QGEN(3,K),QGEN(1,K),QGEN(2,K),EMSQ
4505       FORMAT(/' Q**2  = ',E11.4,' * TAN(',E11.4,' + ',E11.4,
     $    ' * RANF) + ',E11.4)
        ENDIF
        WRITE(ITLIS,4506) QSELWT(K)
4506    FORMAT(/' WEIGHT = ',E11.4)
440   CONTINUE
C
C          Set fixed limits if any
C
      IF(FIXQT) THEN
        PTMAX(3)=PTMIN(3)
        PT(3)=PTMIN(3)
        QTW=PT(3)
      ENDIF
      IF(FIXQM) THEN
        QMAX=QMIN
        QMW=QMIN
      ENDIF
      RETURN
C
C          Fit fails if SIGMA=0 in allowed range
C
999   WRITE(ITLIS,9990) QMW,QTW
9990  FORMAT(//' ERROR IN QFUNC...SIGMA=0 FOR QMW = ',E12.4,' , QTW = ',
     1E12.4/' CHECK YOUR LIMITS')
      STOP 99
      END
Commit	Line	Data
0795afa3	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/QMAX2)QPOW/(PT2+RNU2)PTPOW
	8	C Resonant:
	9	C SIGMA=ANORM/((Q2-M2)2+M2GAM*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 Q2 fit by QMAX2 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)=CUTOFFQMW*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/QMAX2)QPOW/(PT2+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-EM2)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=NSUMS(7)-SUMS(1)SUMS(2)
174	T2=NSUMS(5)-SUMS(2)*2
175	T3=NSUMS(4)-SUMS(1)*2
176	T4=NSUMS(8)-SUMS(1)SUMS(3)
177	T5=NSUMS(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-T2T3
189	PTPOW(K)=(T5T3-T4T1)/DET
190	QPOW(K)=(T4T2-T1T5)/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-EM2)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 PT2 and QMW2
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*2QSELWT(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