C C ------------------ SOURCE CODE for RIDI 2.0 ----------------------- C C --------------------------------------------------------------------- C C +KEEP kept for documentation C Commons included by means of INCLUDE statement C C --------------------------------------------------------------------- C+KEEP,RIDSEQ. C C logical duejet,trejet,fragm,wei,two,longit,transv C + ,user ! not in ZDIS C logical lightf,thlim,charm,strange,checkout,onlych C COMMON/LOGICA/duejet,trejet,fragm,wei,two,longit,transv, C + user, ! not in ZDIS C + lightf,thlim,charm,strange,checkout,onlych C COMMON/EVENTS/nevwei,ievent,jevent,twojets,iqqbar C COMMON/WEI2/nevent2,itoteve2,itwoev2,iwei2,nmass,fk_2,weimax2, C + sig2t,sig2l,weit2,weil2,twint2(7),twinte2(7), C + cs2(2),cst2(2),csl2(2),tcs2(2),tcse2(2), C + tcst2(2),tcste2(2),tcsl2(2),tcsle2(2), C + tcsb2,tcsbe2,tcsi2,tcsie2,tcsa2,tcsae2,tcsf2,tcsfe2 C COMMON/WEI3/nevent3,ITOTEVE3,itwoev3,IWEI3,nskip,fk_3,weimax3, C + weit3,weil3,TWINT3(7),TWINTE3(7), C + cs3(2),cst3(2),csl3(2),tcs3(2),tcse3(2), C + tcst3(2),tcste3(2),tcsl3(2),tcsle3(2), C + tcsb3,tcsbe3,tcsi3,tcsie3,tcsa3,tcsae3,tcsf3,tcsfe3 C C COMMON/MASSES/AMP,AME,AMG,AMPSQ,AMESQ,AMGSQ,CMCBAR,CNBARN,PI, C + HBARC,HBARC2,CM2DIM,ALPHA,ALPHA2,ALPSCUT, C + effmas(4),cthre,sthre C COMMON/BEAMS/EBEAM,PBEAM,EB_I(5),PB(5),STOT(5),LFLAG C COMMON/KCUTS/QSQC(2),TCUTS(2),AMX2C(2),XPOMC(2), C + AKT2C(2),VCUTS(12),xbjddc,betamx C COMMON/DIFFR/BSLOPE,CTH2,CTH C COMMON/LUND/IPS,idq C COMMON/GENERA/PP,RR,SS,AA,BB,BBB,CC,CCC,EEE,FF C COMMON/CHECK/laacheck,nptcms,PT2CMS, C + IAKT2,IQSQ,IEXBJ,IAMX2,IZ,IAMF2,ITHJCM,ialps, C + akt2lim,xpomlim(2),amx2lim(2),chmin,stmin C COMMON/CHOICE/akt2cut,ikfact,iwg,fks(2) C COMMON/GLUON/igluon,gval(2),gsfk02,gpar(4) C C logical radcor,radcor_i,radcor_f,radcor_a C common/radlog/radcor,radcor_i,radcor_f,radcor_a C common/radcon/ephmin,thphmn C c c NOT IN ZDIS c C COMMON/EVENTI/NEVENT,NEVPRI C COMMON/SEME/ISEED C C+KEEP,RIDKIN. C C COMMON/KINGEN/QSQ,QSQMIN,QSQMAX,FACQSQ,PFQSQ, C + XBJ,XBJMIN,XBJMAX,FACXBJ,PFXBJ, C + EXBJ,EXBJMN,EXBJMX,FACEXBJ,PFEXBJ, C + T,TMIN,TMAX,FACT,PFT,WTT, C + PT2,PT2MIN,PT2MAX,FACPT2,PFPT2,WTPT2, C + AMX2,AMX2MN,AMX2MX,FACMX2,PFMX2,WTMX2, C + XF,XFMIN,XFMAX,FACXF,PFXF,WTXF,xpommn, C + AKT2,AKT2MN,AKT2MX,FACKT2,PFKT2,WTKT2, C + Z,ZMIN,ZMAX,FACZ,PFZ,WTZ, C + THJCM,PHJCM,FACTHJ,PFTHJC,PFPHJC,WTTHJ, C + THJMN2,THJMX2,WTQINT C COMMON/VECTS/EB(5),ES(5),PS(5),GAM(5),POM(5),XVECT(5), C + AK(5),GLUON(5),YVECT(5),QQBAR(5),GLUECM(5), C + Q1CM(5),QINT(5),Q1(5),Q2(5),Q(5),QBAR(5) C COMMON/KINEM/S,WSQ,Y,EY,VU,PT,XL,XL_CMS,epoms, C + AMX,XPOM,XG,XK,AKT,AMKT2, C + THETAP,PHIP,PHIPCM,PFPHIE,PFPHIP, C + THETAE,PHIE, C + PHIK,THETAK,PHIAK,THEKCM,PHIKCM, C + EJET,PJET,THEQ,PHIQ,THEQB,PHIQB,IFLAV,q_kt2 C COMMON/WEIGHT/WEIGHT,WPHASE,GAMT,EPSILON C common/radkin/zph,phiph,theph,es_f(5),phot(5),elec_f(5),srad(5) C C+KEEP,RBOOST. C C IMPLICIT DOUBLE PRECISION (D) C COMMON/BOOST/PA(4,5),DBETA(1,3),DTHETA(1),DPHI(1),BETA(1,3), C + THETA(1),PHI(1),ibskip C C --------------------------------------------------------------------- C BLOCK DATA MCRIDI C --------------------------------------------------------------------- C Block data for RIDI C C Author : Ada Solano C --------------------------------------------------------------------- C include 'ridseq.inc' C C Constants C DATA AMP / 0.938272 / DATA AME / 0.000511 / DATA AMG / 0.75 / DATA CNBARN / 10.0E32 / DATA CMCBAR / 1.0E30 / DATA HBARC / 1.97E-14 / DATA CTH2 / 3.5 / DATA ALPSCUT / 0.7 / C C Default beam energies C DATA EBEAM / 27.5 / DATA PBEAM / 820. / C C Default kinematical cuts C DATA QSQC / 0.1, 10000. / DATA TCUTS / -1., -0.00001 / DATA AMX2C / 3., 10000. / DATA XPOMC / 0.00001, 0.1 / DATA AKT2C / 0.001, 100. / C C Choose the maximum x_Bjorken for diffractive events C ( xbjddc < 0.1 in HERA kinematics! ) C data xbjddc / 0.1 / C C Default setup C DATA LFLAG / 1 / ! 1=positron, -1=electron data duejet /.true./ ! q-qbar final states data trejet /.true./ ! q-qbar-gluon final states data transv /.true./ ! trasverse cross section data longit /.true./ ! longitudinal cross section data lightf /.false./ ! light flavours only data thlim /.false./ ! high k_t only data fragm /.true./ ! fragmentation data wei /.false./ ! generation of weighted events data nevwei / 1000 / ! events for calculating the xsect's data checkout /.false./ ! more output for check data radcor /.false./ ! radiative corrections: data radcor_i /.false./ ! initial state radiation data radcor_f /.false./ ! final state radiation data radcor_a /.false./ ! vertex correction data onlych /.false./ ! charm only C C down, up, strange, charm data effmas / 0.0099, 0.0056, 0.199, 1.35 / ! current mass (JETSET) data cthre / 2.05 / data sthre / 0.97 / C C Choose type of k-factors: 0 = constant k-factors C 1 = exp[alphas(MX^2)*pi*C_F(A)] C 2 = exp[alphas(MX^2/4)*pi*C_F] for qqbar C 2 = exp[alphas(MX^2/4)*pi*C_A]/3+ for qqbarg C exp[alphas(MX^2/4)*pi*C_F]*2/3 C 3 = (1.+alphas(MX^2/4)*pi*C_F(A)/2.)**2 data ikfact / 2 / data fk_2 / 2.1 / ! constant k-factors data fk_3 / 4. / C C Constants C data bslope / 6. / ! slope of the p_t^2 distribution data akt2cut / 0.06 / ! infrared cutoff for k_t^2' data ephmin /0.03/ ! minimum energy and theta angle (FSR) data thphmn /0.0001/ ! of a radiated photon C C Set-up C data igluon / 1 / ! 0 = simple xg(x); 1 = use PDFLIB C C If igluon=0: xg(x)=gpar(1)*(1-x)**gpar(2)*(x/gpar(3))**gpar(4) C data gpar / 3., 5., 0.01, -0.4 / C C Gluon structure function choice according to PDFLIB: C gval(1) indicates Ngroup (3=MRS, 5=GRV) C gval(2) allows the choice of the available Nset for each Ngroup C data gval / 5., 6. / data gsfk02 / 1. / ! cutoff for linear extrapolation C C w_g=(1-z')^3*(2z'+1)^2 with z'=x_Bj/((ex_Bj/x_gluon)*xpom) C takes into account a more precise form of the triple pomeron vertex C at not too small values of z' (Levin&Wusthoff 1994) C data iwg / 1 / ! 1: w_g in the cross section C C Type of parton shower for the q-qbar-gluon final state: C 0 -> parton shower only for the q-qbar system :qmax=sqrt(m_f^2+p_t^2) C 1 -> parton shower separately for the gluon jet | C (qmax=sqrt(m_g^2+k_t^2)) and the q-qbar system <-| C DATA IPS / 1 / C C Constants and powers for variable generation (unweighted events) C C q-qbar diagram data pp / 0.5 / ! 1/[x_bj^(PP+1)] data rr / 1. / ! 1/[Q^2^(RR+1)] data ss / 1. / ! 1/[(Mx^2+Q^2)^(SS+1)] C q-qbar-gluon diagram DATA AA / 0.05 / ! 1/[(1-x_F)^(AA+1)] DATA BB / 1.5 / ! 1/[(k'_t^2+m_g^2/BBB)^(BB+1)] DATA BBB / 0.025 / DATA CC / 1.5 / ! 1/[(Q^2+CCC)^(CC+1)] DATA CCC / 5. / DATA EEE / 0.03 / ! 1/[exbj*(exbj+EEE)] DATA FF / 0. / ! (1-z)^FF, where z=exbj/x_gluon C C Used for checks C data two /.false./ data chmin / 100. / data stmin / 100. / data akt2lim / 100. / data xpomlim / 1., 0. / data amx2lim / 10000., 0. / C C NOT IN ZDIS C DATA NEVENT / 10000 / data nevpri / 0 / ! events to be printed data user /.false./ ! user output C END C C ------------------------------------------------------------------- C SUBROUTINE RIDINI(JErr) C ------------------------------------------------------------------- C Initialization of RIDI C Called by : MAIN C C Author : Ada Solano C ------------------------------------------------------------------- C JErr = 0 c WRITE(6,*) WRITE(6,*)' *****************************************************' WRITE(6,*)' ** **' WRITE(6,*)' ** RIDI Version 2.0 **' WRITE(6,*)' ** **' WRITE(6,*)' ** A MONTE CARLO GENERATOR **' WRITE(6,*)' ** FOR DIFFRACTIVE DISSOCIATION **' WRITE(6,*)' ** OF VIRTUAL PHOTONS INTO Q-QBAR AND **' WRITE(6,*)' ** Q-QBAR-GLUON FINAL STATES **' WRITE(6,*)' ** BASED ON THE RYSKIN MODEL **' WRITE(6,*)' ** **' WRITE(6,*)' ** Written by Ada Solano **' WRITE(6,*)' ** Comments and questions to solano@vxdesy.desy.de **' WRITE(6,*)' ** **' WRITE(6,*)' *****************************************************' WRITE(6,*) c c 1) Initialize c CALL INIT(JErr) if (JErr.ne.0) return c c 2) Use the first NEVWEI events to calculate the maximum weight c CALL WEIMAX c c 3) Use NEVWEI events to determine the 2jet and 3jet cross sections c CALL JETNUM c RETURN END C C ------------------------------------------------------------------- C SUBROUTINE INIT(IErr) C ------------------------------------------------------------------- C Initialization of variables and kinematics C Called by : RIDINI C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' common/ludat2/kchg(500,3),pmas(500,4),parf(2000),vckm(4,4) EXTERNAL LUDATA C integer kc C IErr = 0 C C Initialize constants and variables C AMPSQ = AMP*AMP AMESQ = AME*AME AMGSQ = AMG*AMG PI = 4.0*ATAN(1.0) ALPHA = 1./134.5 ! alpha_em at about Q^2=10 GeV^2 ALPHA2 = ALPHA*ALPHA HBARC2 = HBARC*HBARC CM2DIM = HBARC2*CNBARN C CTH = SQRT(CTH2) C C Read random number seed from external file (unit 10) C read(10,*,ERR=10) iseed rewind(10) 10 if(iseed.eq.0) then write(6,*) ' ' write(6,*) 'WARNING-INIT: Missing file RIDI.RND' write(6,*) ' random number seed = 0' write(6,*) ' ' endif C C Read/write data cards (RIDI and JETSET) C CALL READCARDS CALL RIDICARDS C if(ikfact.eq.0) then fk_2 = fks(1) fk_3 = fks(2) endif C if(duejet) then weimax2 = 0 ITOTEVE2 = 0 ITWOEV2 = 0 IWEI2 = 0 NEVENT2 = 0 nmass = 0 laacheck = 0 nptcms = 0 tcsb2 = 0. tcsi2 = 0. tcsa2 = 0. tcsf2 = 0. tcsbe2 = 0. tcsie2 = 0. tcsae2 = 0. tcsfe2 = 0. do i=1,2 tcs2(i) = 0. tcse2(i) = 0. tcst2(i) = 0. tcste2(i) = 0. tcsl2(i) = 0. tcsle2(i) = 0. enddo do i=1,7 twint2(i) = 0. twinte2(i) = 0. enddo endif C if(trejet) then weimax3 = 0 ITOTEVE3 = 0 ITWOEV3 = 0 IWEI3 = 0 NEVENT3 = 0 IAKT2 = 0 IQSQ = 0 IEXBJ = 0 IAMX2 = 0 IZ = 0 IAMF2 = 0 ITHJCM = 0 IALPS = 0 tcsb3 = 0. tcsi3 = 0. tcsa3 = 0. tcsf3 = 0. tcsbe3 = 0. tcsie3 = 0. tcsae3 = 0. tcsfe3 = 0. do i=1,2 tcs3(i) = 0. tcse3(i) = 0. tcst3(i) = 0. tcste3(i) = 0. tcsl3(i) = 0. tcsle3(i) = 0. enddo do i=1,7 twint3(i) = 0. twinte3(i) = 0. enddo endif C if(akt2cut.gt.0.) then cthre = sqrt(4.16+akt2cut) sthre = sqrt(0.9+akt2cut) endif C if(fragm) then c c Change here JETSET current algebra masses c kc = lucomp(4) pmas(kc,1) = effmas(4) c endif C C Initialize kinematics C CALL CUTINI(IErr) if (IErr.ne.0) return C C Analysis initialization and histogram booking (not in ZDIS) C if(user) CALL USERINI C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE READCARDS C ------------------------------------------------------------------- C Read data cards C Called by : INIT C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' COMMON/CFREAD/SPACE(10000) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) EXTERNAL LUDATA C C ** Read RIDI generator data cards from external file (unit 7) C call ffinit(10000) call ffset('SIZE',6) call ffset('LINP',7) C C .true. = Radiative corrections on C call ffkey('RADCOR',radcor,1,'LOGI') C C .true. = Generation with weitghs C call ffkey('WEIGHT',wei,1,'LOGI') C C .true. = q-qbar and/or q-qbar-gluon final state C call ffkey('QQBAR',duejet,1,'LOGI') call ffkey('QQBARG',trejet,1,'LOGI') C C .true. = Transverse and/or longitudinal cross section C call ffkey('TRANSV',transv,1,'LOGI') call ffkey('LONGIT',longit,1,'LOGI') C C .true. = Light flavours only C call ffkey('LIGHTF',lightf,1,'LOGI') C C .true. = Charm flavour only C call ffkey('ONLYCH',onlych,1,'LOGI') C C .true. = High k_t only in q-qbar final state C call ffkey('THLIM',thlim,1,'LOGI') C C .true. = Fragmentation on C call ffkey('FRAGM',fragm,1,'LOGI') C C .true. = More output with check variables C call ffkey('CHECK',checkout,1,'LOGI') C C .true. = User output on (not in ZDIS) C call ffkey('USER',user,1,'LOGI') C C Number of final events (not in ZDIS) C call ffkey('NEVENT',nevent,1,'INTE') C C Number of events to be printed (not in ZDIS) C call ffkey('NEVPRI',nevpri,1,'INTE') C C Number of events to be used for the cross section determination C (in case of generation of unweighted events) C call ffkey('NEVWEI',nevwei,1,'INTE') C C K-factor choice: 0 = constant k-factors C 1 = exp[alphas(MX^2)*pi*C_F(A)] C 2 = exp[alphas(MX^2/4)*pi*C_F] for qqbar C 2 = exp[alphas(MX^2/4)*pi*C_A]/3+ for qqbarg C exp[alphas(MX^2/4)*pi*C_F]*2/3 C 3 = (1.+alphas(MX^2/4)*pi*C_F(A)/2.)**2 C call ffkey('IKFACT',ikfact,1,'INTE') C C Values of the constant K-factors if IKFACT = 0 C (default: fks(1) = 2.1; fks(2) = 4. ) C call ffkey('KFS',fks,2,'REAL') C C Type of parton shower for the q-qbar-gluon final state: C 0 -> parton shower only for the q-qbar system :qmax=sqrt(m_f^2+p_t^2) C 1 -> parton shower separately for the gluon jet | C (qmax=sqrt(m_g^2+k_t^2)) and the q-qbar system <-| C call ffkey('IPSTYP',ips,1,'INTE') C C w_g=(1-z')^3*(2z'+1)^2 with z'=x_Bj/((ex_Bj/x_gluon)*xpom) C takes into account a more precise form of the triple pomeron vertex C at not too small values of z' (Levin&Wusthoff 1994) C 0 -> w_g NOT in the cross section C 1 -> w_g in the cross section C call ffkey('IWG',iwg,1,'INTE') C C Lepton choice: 1=positron; -1=electron C call ffkey('LFLAG',lflag,1,'INTE') C C Lepton beam momentum (positive!) C call ffkey('EBEAM',ebeam,1,'REAL') C C Proton beam momentum (positive!) C call ffkey('PBEAM',pbeam,1,'REAL') C C Kinematical limits on variables: Q^2, t, Mx^2, x_pom and k_t^2 C call ffkey('Q2LIM',qsqc,2,'REAL') call ffkey('TLIM',tcuts,2,'REAL') call ffkey('MX2LIM',amx2c,2,'REAL') call ffkey('XPOLIM',xpomc,2,'REAL') call ffkey('KT2LIM',akt2c,2,'REAL') C C Maximum x_Bjorken for diffractive events: xbjddc < xpom_max < 0.1 C call ffkey('XBJDDC',xbjddc,1,'REAL') C C Slope of the p_t^2 distribution of the scattered proton C call ffkey('BSLOPE',bslope,1,'REAL') C C Charm mass (default: m_c = 1.35 GeV) C call ffkey('CHMASS',effmas(4),1,'REAL') C C Gluon effective mass (default: m_g = 0.75 GeV) C call ffkey('GMASS',amg,1,'REAL') C C Infrared cutoff for k'_t^2 (default: 0.06 GeV^2) C call ffkey('KT2IC',akt2cut,1,'REAL') C C Cut on alpha_s: if(alpha_s.gt.alpscut) alpha_s = alpscut C call ffkey('ALPSC',alpscut,1,'REAL') C C Gluon structure function choice: 0 = simple xg(x); 1 = PDFLIB C call ffkey('IGLUON',igluon,1,'INTE') C C In case of igluon=0 choose parameters of C xg(x)=gpar(1)*(1-x)**gpar(2)*(x/gpar(3))**gpar(4) C call ffkey('GPAR',gpar,4,'REAL') C C In case of igluon=1 choose gluon distribution from PDFLIB C (default GRV HO) C call ffkey('GVAL',gval,2,'REAL') C C Cutoff for the gluon structure function: C linear extrapolation for k'_t^2 < k_0^2 (default: 1 GeV^2) C call ffkey('GSFK02',gsfk02,1,'REAL') C C Minimum energy of the radiated photon in case of radiative events C call ffkey('EPHMIN',ephmin,1,'REAL') C C Minimum angle of the radiated photon w.r.t. to the incoming lepton C (FSR events only) C call ffkey('THPHMN',thphmn,1,'REAL') C C To optimaze variable generation in case of unweighted events C choose here the generation functions: C C q-qbar diagram: C C 1/[x_bj^(PP+1)] C call ffkey('PP',pp,1,'REAL') C C 1/[Q^2^(RR+1)] C call ffkey('RR',rr,1,'REAL') C C 1/[(Mx^2+Q^2)^(SS+1)] C call ffkey('SS',ss,1,'REAL') C C q-qbar-gluon diagram: C C 1/[(1-x_F)^(AA+1)], where x_F = p_z'/p_z C call ffkey('AA',aa,1,'REAL') C C 1/[(k'_t^2+m_g^2/BBB)^(BB+1)] C call ffkey('BB',bb,1,'REAL') call ffkey('BBB',bbb,1,'REAL') C C 1/[(Q^2+CCC)^(CC+1)] C call ffkey('CC',cc,1,'REAL') call ffkey('CCC',ccc,1,'REAL') C C 1/[exbj*(exbj+EEE)], where exbj=x_bj in the electron-pomeron cms C call ffkey('EEE',eee,1,'REAL') C C (1-z)^FF, where z=exbj/x_gluon C call ffkey('FF',ff,1,'REAL') C C C ** Initialize standard JETSET default values before reading the cards C MSTU(12) = 0 ! Avoid the header message C C ** Read JETSET data cards C call ffkey('MSTU',mstu,200,'INTE') call ffkey('MSTJ',mstj,200,'INTE') call ffkey('PARU',paru,200,'REAL') call ffkey('PARJ',parj,200,'REAL') C C ** End data card reading C call ffgo C return end C C ------------------------------------------------------------------- C SUBROUTINE RIDICARDS C ------------------------------------------------------------------- C Write data cards C Called by : INIT C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' C write(6,*)' ' write(6,*)'RIDI data cards:' write(6,*)' ' write(6,*)'RIDI-RADCOR = ',radcor write(6,*)'RIDI-WEIGHT = ',wei write(6,*)'RIDI-QQBAR = ',duejet write(6,*)'RIDI-QQBARG = ',trejet write(6,*)'RIDI-TRANSV = ',transv write(6,*)'RIDI-LONGIT = ',longit write(6,*)'RIDI-LIGHTF = ',lightf write(6,*)'RIDI-ONLYCH = ',onlych write(6,*)'RIDI-THLIM = ',thlim write(6,*)'RIDI-FRAGM = ',fragm write(6,*)'RIDI-CHECK = ',checkout C C Not in ZDIS C write(6,*)'RIDI-USER = ',user write(6,*)'RIDI-NEVENT = ',nevent write(6,*)'RIDI-NEVPRI = ',nevpri C write(6,*)'RIDI-NEVWEI = ',nevwei write(6,*)'RIDI-IKFACT = ',ikfact if(ikfact.eq.0) then write(6,*)'RIDI-KFS = ',fks(1),' ',fks(2) endif write(6,*)'RIDI-IWG = ',iwg write(6,*)'RIDI-IPSTYP = ',ips write(6,*)'RIDI-LFLAG = ',lflag write(6,*)'RIDI-EBEAM = ',ebeam write(6,*)'RIDI-PBEAM = ',pbeam write(6,*)'RIDI-XBJDDC = ',xbjddc write(6,*)'RIDI-Q2LIM = ',qsqc(1),' ',qsqc(2) write(6,*)'RIDI-TLIM = ',tcuts(1),' ',tcuts(2) write(6,*)'RIDI-MX2LIM = ',amx2c(1),' ',amx2c(2) write(6,*)'RIDI-XPOLIM = ',xpomc(1),' ',xpomc(2) write(6,*)'RIDI-KT2LIM = ',akt2c(1),' ',akt2c(2) write(6,*)'RIDI-BSLOPE = ',bslope write(6,*)'RIDI-CHMASS = ',effmas(4) write(6,*)'RIDI-GMASS = ',amg write(6,*)'RIDI-ALPSC = ',alpscut write(6,*)'RIDI-KT2IC = ',akt2cut write(6,*)'RIDI-IGLUON = ',igluon write(6,*)'RIDI-GPAR = ',gpar(1),gpar(2),gpar(3),gpar(4) write(6,*)'RIDI-GVAL = ',gval(1),' ',gval(2) write(6,*)'RIDI-GSFK02 = ',gsfk02 write(6,*)'RIDI-EPHMIN = ',ephmin write(6,*)'RIDI-THPHMN = ',thphmn write(6,*)'RIDI-PP = ',pp write(6,*)'RIDI-RR = ',rr write(6,*)'RIDI-SS = ',ss write(6,*)'RIDI-AA = ',aa write(6,*)'RIDI-BB = ',bb write(6,*)'RIDI-BBB = ',bbb write(6,*)'RIDI-CC = ',cc write(6,*)'RIDI-CCC = ',ccc write(6,*)'RIDI-EEE = ',eee write(6,*)'RIDI-FF = ',ff write(6,*)' ' C return end C C ------------------------------------------------------------------- C SUBROUTINE CUTINI(KErr) C ------------------------------------------------------------------- C Evaluation of the kinematical limits C Called by : INIT C C Author : Ada Solano C ------------------------------------------------------------------- C Implicit Double Precision (D) include 'ridseq.inc' include 'ridkin.inc' C real betamn C KErr = 0 C C Vectors of incident particles in the Zeus reference system C C Vector (5) = ( p_x, p_y, p_z, energy, mass**2 ) C EB_I(1) = 0. EB_I(2) = 0. EB_I(3) = -EBEAM EB_I(5) = AME*AME EB_I(4) = SQRT(EB_I(3)*EB_I(3)+EB_I(5)) PB(1) = 0. PB(2) = 0. PB(3) = PBEAM PB(5) = AMP*AMP PB(4) = SQRT(PB(3)*PB(3)+PB(5)) c do 1 i=1,5 1 eb(i) = eb_i(i) c c Total center of mass energy [ s = 4.*ebeam*pbeam ] c [ s = ampsq + 2.*amp*ebeam ] c call add5(eb_i,pb,stot) s = stot(5) c c Variable cuts (user vs. kinematics of diffractive events) c tmin = max(tcuts(1),-1.) tmax = min(tcuts(2),-0.00001) if(tmin.ge.tmax)goto 999 c c Avoid at least the resonance region ( Mx2>3 ) c A lower limit can be used since JETSET does not complain with c new PARJ(32), but still the fragmentation can be wrong (15.5.96) c amx2c(1) = max(amx2c(1),3.) c c Limits on Q**2 as function of s and m_lepton only ... c ... and the Mx**2 minimum which has to be generated c amlsq = amesq amxmsq = (amp+sqrt(amx2c(1)))*(amp+sqrt(amx2c(1))) daa = dble(s + amlsq - ampsq) daaasq = daa*daa-4.d0*s*amlsq daaa = dsqrt(daaasq) dbb = dble(s + amlsq - amxmsq) dbbbsq = dbb*dbb-4.d0*s*amlsq dbbb = dsqrt(dbbbsq) dqqmin = -2.d0*amlsq + (daa*dbb-daaa*dbbb)/s/2. dqqmax = -2.d0*amlsq + (daa*dbb+daaa*dbbb)/s/2. qqmin = sngl(dqqmin) qqmax = sngl(dqqmax) c c Combined Q**2 limits + Bjorken_x and Nu limits c qsqc(1) = max(qqmin,qsqc(1),0.1) qsqc(2) = min(qqmax,qsqc(2)) if(qsqc(1).ge.qsqc(2))goto 999 xbjmin = qsqc(1)/(s-ampsq-amesq) xbjmax = xbjddc if(xbjmin.ge.xbjmax)goto 999 c c xpom_min cannot be lower than xbj_min since beta<=1 c betamx = qsqc(1)/(qsqc(1)+amx2c(1)-tmax) betamn = qsqc(1)/(qsqc(1)+amx2c(2)-tmin) xpomc(1) = max(xpomc(1),xbjmin/betamx) xpomc(2) = min(xpomc(2),xbjddc/betamn) if(xpomc(1).ge.xpomc(2))goto 999 xpommn = xpomc(1) c if(trejet) then XFMIN = 1. - xpomc(2) XFMAX = 1. - xpomc(1) IF(XFMIN.GE.XFMAX)GOTO 999 PT2MIN = max(-tcuts(2),0.00001) PT2MAX = min(-tcuts(1),1.) IF(PT2MIN.GE.PT2MAX)GOTO 999 akt2c(1) = max(akt2c(1),0.001) AKT2MN = AKT2C(1) IF(AKT2MN.GE.AKT2C(2))GOTO 999 THJMX2 = PI*PI endif c c Fill vector VCUTS contening the final kinematical cuts c vcuts(1) = xbjmin vcuts(2) = xbjmax vcuts(3) = qsqc(1) vcuts(4) = qsqc(2) vcuts(5) = tmin vcuts(6) = tmax vcuts(7) = amx2c(1) vcuts(8) = amx2c(2) vcuts(9) = xpomc(1) vcuts(10)= xpomc(2) if(trejet) then vcuts(11)= AKT2C(1) vcuts(12)= AKT2C(2) endif C RETURN C 999 WRITE(6,*)'ERROR-INIT: WRONG VARIABLE CUT! STOP GENERATION' KErr = 1 C return end C C ------------------------------------------------------------------- C SUBROUTINE WEIMAX C ------------------------------------------------------------------- C Evaluation of the maximum weight C Called by : RIDINI C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' c if(wei) return ievent = 0 jevent = 1 DO 1 WHILE(JEVENT.LE.NEVWEI) c c Generate the partonic system and calculate the event weight c if(duejet) then CALL EVENT2(wei2) if(wei2.gt.weimax2) weimax2 = wei2 endif if(trejet) then CALL EVENT3(wei3) if(wei3.gt.weimax3) weimax3 = wei3 endif c jevent = jevent+1 1 continue c if(duejet) weimax2 = 1.5*weimax2 if(trejet) weimax3 = 1.5*weimax3 c return end C C ------------------------------------------------------------------- C SUBROUTINE JETNUM C ------------------------------------------------------------------- C q-qbar C Evaluation of the cross section ratio ---------------- C q-qbar-gluon C Called by : RIDINI C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' c nmass = 0 nskip = 0 ievent = 0 jevent = 1 c if(duejet.and.(.not.trejet)) then twojets = 1. goto 10 elseif(trejet.and.(.not.duejet)) then twojets = 0. goto 10 endif c DO 2 WHILE(JEVENT.LE.NEVWEI) c c Generate the partonic system and count the event cross section c if(duejet) then 3 CALL EVENT2(wei2) CALL CROSS2(jagain) if(jagain.ne.0) goto 3 endif c if(trejet) then 4 CALL EVENT3(wei3) CALL CROSS3(jagain) if(jagain.ne.0) goto 4 endif c jevent = jevent+1 2 continue c c 2jet and 3jet cross sections c if(.not.wei) then c if(duejet) NEVENT2 = ITOTEVE2+nmass if(trejet) NEVENT3 = ITOTEVE3+NSKIP c else c if(duejet) NEVENT2 = NEVWEI+nmass if(trejet) NEVENT3 = NEVWEI+NSKIP c endif c if(duejet) TCS2(1) = TCS2(1)/NEVENT2 if(trejet) TCS3(1) = TCS3(1)/NEVENT3 c c Fraction of 2jet events to be generated c twojets = TCS2(1)/( TCS2(1)+TCS3(1) ) c 10 continue if(duejet) then ITOTEVE2 = 0 ITWOEV2 = 0 IWEI2 = 0 nmass = 0 NEVENT2 = 0 tcsb2 = 0. tcsi2 = 0. tcsa2 = 0. tcsf2 = 0. tcsbe2 = 0. tcsie2 = 0. tcsae2 = 0. tcsfe2 = 0. do i=1,2 tcs2(i) = 0. tcse2(i) = 0. tcst2(i) = 0. tcste2(i) = 0. tcsl2(i) = 0. tcsle2(i) = 0. enddo endif if(trejet) then ITOTEVE3 = 0 ITWOEV3 = 0 IWEI3 = 0 NSKIP = 0 NEVENT3 = 0 tcsb3 = 0. tcsi3 = 0. tcsa3 = 0. tcsf3 = 0. tcsbe3 = 0. tcsie3 = 0. tcsae3 = 0. tcsfe3 = 0. do i=1,2 tcs3(i) = 0. tcse3(i) = 0. tcst3(i) = 0. tcste3(i) = 0. tcsl3(i) = 0. tcsle3(i) = 0. enddo endif c ievent = 1 two = .false. c RETURN END C C ------------------------------------------------------------------- C SUBROUTINE RIDIEVT C ------------------------------------------------------------------- C Main routine of RIDI C Called by : MAIN C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' c real choice Jerr = 0 c if(two) then c c weight.gt.weimax for the previous generated event c -> count the event two times in accordance with probability c two = .false. if(iqqbar.eq.1) then c if(fragm) CALL FRAGM2 c c User output of event variables (not in ZDIS) c if(user) CALL USEREVT c IEVENT = IEVENT+1 NEVENT2 = NEVENT2+1 c else c if(fragm) CALL FRAGM3 c c User output of event variables (not in ZDIS) c if(user) CALL USEREVT c IEVENT = IEVENT+1 NEVENT3 = NEVENT3+1 c endif goto 99 endif c c Start event generation c iqqbar = 0 choice = userran(iseed) if(choice.lt.twojets) then iqqbar = 1 c c Generate a q-qbar partonic system c 11 CALL EVENT2(wei2) c c Event cross section c CALL CROSS2(jagain) if(jagain.ne.0) goto 11 c c Insert initial and scattered electrons and protons c + fragmentation of the partonic system c if(fragm) CALL FRAGM2 c c User output of event variables (not in ZDIS) c if(user) CALL USEREVT c IEVENT = IEVENT+1 NEVENT2 = NEVENT2+1 c else c c Generate a q-qbar-gluon partonic system c 13 CALL EVENT3(wei3) c c Event cross section c CALL CROSS3(jagain) if(jagain.ne.0) goto 13 c c Insert initial and scattered electrons and protons c + fragmentation of the partonic system c if(fragm) CALL FRAGM3 c c User output of event variables (not in ZDIS) c if(user) CALL USEREVT c IEVENT = IEVENT+1 NEVENT3 = NEVENT3+1 c endif c 99 continue c Not in ZDIS if(mod(ievent,100).eq.0) + write(6,*)'Event number ',ievent,' Last random number ',iseed c return end C C ------------------------------------------------------------------- C SUBROUTINE RADCORR(jtype) C ------------------------------------------------------------------- C Routine for the generation of radiative events C Called by : RIDIEVT C C jtype = 0 : initial state radiation C jtype = 1 : final state radiation C jtype = 2 : absorption of initial state radiation C C Author : Ada Solano C ------------------------------------------------------------------- C include 'rboost.inc' include 'ridseq.inc' include 'ridkin.inc' c real elec_i(5) c radcor_a = .false. if(jtype.eq.2) then if(theph.ge.thetae) then radcor_a = .true. endif return endif c radcor_i = .false. radcor_f = .false. c if(jtype.eq.0) then do 10 i=1,5 10 elec_i(i) = eb_i(i) else if(thetae.lt.thphmn) return do 20 i=1,5 20 elec_i(i) = es(i) endif c c Generate the fraction of photon energy w.r.t. the electron: c zph = E_photon/E_e c zphmax = 1. faczph = alog(zphmax*elec_i(4)/ephmin) zph = (ephmin/elec_i(4))*exp(faczph*userran(iseed)) c c Generate the azimuthal angle of the radiated photon c phiph = 2*pi*userran(iseed) c c Generate the theta angle of the radiated photon c if(jtype.eq.0) then factph = alog((elec_i(4)*elec_i(4)*pi*pi+amesq)/amesq) theph = sqrt(exp(factph*userran(iseed))-1.)*ame/elec_i(4) else factph = alog( (elec_i(4)*elec_i(4)*thetae*thetae+amesq)/ & (elec_i(4)*elec_i(4)*thphmn*thphmn+amesq) ) theph = sqrt( (elec_i(4)*elec_i(4)*thphmn*thphmn+amesq)* & exp(factph*userran(iseed)) - amesq )/elec_i(4) endif c if(theph.lt.0.01) then vv = ((amesq/elec_i(4)/elec_i(4) + theph*theph)/2.)**2 else vv = (1.-cos(theph))**2 endif c if(theph.lt.1.e-10) then write(6,*) & 'WARNING-RADCORR: Theta_photon is about 0 ',ievent,theph theph = 1.e-10 endif ww = faczph*factph*(alpha/pi/2.)*(sin(theph)**3/vv)* & (amesq/elec_i(4)/elec_i(4) + theph*theph)/theph/2. if(ww.gt.1.) write(6,*)'WARNING-RADCORR: ww > 1 ',ievent,ww c radcheck = userran(iseed) if(radcheck.ge.ww) return c c Generate a radiative event c if(jtype.eq.0) then radcor_i = .true. else radcor_f = .true. endif c c Fill vector of radiated photon (Zeus reference system) c phot(5) = 0. phot(4) = zph*elec_i(4) phot(1) = phot(4)*sin(theph)*cos(phiph) phot(2) = phot(4)*sin(theph)*sin(phiph) phot(3) = -phot(4)*cos(theph) if(jtype.eq.1) then DO 1 I=1,5 PA(1,I) = elec_i(i) PA(2,I) = 0. PA(3,I) = 0. PA(4,I) = 0. 1 CONTINUE DPHI(1)=ULANGL(PA(1,1),PA(1,2)) CALL ROTATE2(0.D+0,-DPHI(1)) DTHETA(1)=ULANGL(PA(1,3),PA(1,1)) CALL ROTATE2(-DTHETA(1),0.D+0) DO 2 I=1,5 2 PA(3,I) = phot(i) pa(3,3) = -pa(3,3) CALL ROTATE2(DTHETA(1),0.D+0) CALL ROTATE2(0.D+0,DPHI(1)) do 3 i=1,5 3 phot(i) = pa(3,i) endif c c New vector for the incoming, off-mass shell, electron c if(jtype.eq.0) call dif5(elec_i,phot,elec_f) c return end C C ------------------------------------------------------------------- C SUBROUTINE EVTINI(Ierr) C ------------------------------------------------------------------- C Event initialization: kinematics and limits C Called by : RIDIEVT C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' include 'ridkin.inc' c real w2min Ierr = 0 c if(.not.radcor_i) then do 1 i=1,5 1 eb(i) = eb_i(i) s = stot(5) else do 2 i=1,5 2 eb(i) = elec_f(i) call add5(eb,pb,srad) s = srad(5) w2min = (amx2c(1)+qsqc(1))/xpomc(2) - qsqc(1) if(s.le.w2min) goto 99 endif c xbjmin = qsqc(1)/(s-ampsq-eb(5)) if(xbjmin.ge.xbjmax) goto 99 c xpommn = max(xpomc(1),xbjmin/betamx) if(xpommn.ge.xpomc(2)) goto 99 c xfmax = 1.-xpommn if(xfmin.ge.xfmax) goto 99 c return c 99 Ierr = 1 c return end C C ------------------------------------------------------------------- C SUBROUTINE EVENT2(sig2) C ------------------------------------------------------------------- C Main routine for the generation of a q-qbar final state C Called by : WEIMAX, JETNUM, RIDIEVT C C Author : Ada Solano C ------------------------------------------------------------------- C include 'rboost.inc' include 'ridseq.inc' include 'ridkin.inc' c c Loop until the generated event fits the kinematical requirements c real psmod,escat,thetae2,costhee,sinthee,phe real efmod,am2prime,esmod c 1000 continue c c Check for events with initial state radiation c if(.not.radcor) goto 1 1001 itype = 0 call radcorr(itype) c c Set kinematics and limits c call evtini(Jerr) if(Jerr.ne.0) then nmass = nmass+1 goto 1001 endif c 1 charm = .false. strange = .false. C C Generate Bjorken_x variable C if(wei) then ! like 1/x C FACXBJ = ALOG(XBJMAX/XBJMIN) XBJ = XBJMIN*EXP(FACXBJ*userran(iseed)) PFXBJ = 1.0/(FACXBJ*XBJ) C else ! like 1/x**(PP+1) C FACXBJ = (XBJMIN/XBJMAX)**PP - 1. XBJ = XBJMIN/((1.+FACXBJ*userran(iseed))**(1./PP)) PFXBJ = -PP*XBJMIN**PP/FACXBJ/(XBJ**(PP+1.)) C endif C C Generate QSQ C QSQMAX = MIN((S-AMPSQ-eb(5))*XBJ,QSQC(2)) C if(wei) then ! like 1/Q2**2 C FACQSQ = 1. - (QSQc(1)/QSQMAX) QSQ = QSQc(1)/(1.-FACQSQ*userran(iseed)) PFQSQ = QSQc(1)/(FACQSQ*QSQ*QSQ) C else ! like 1/Q2**(RR+1) C FACQSQ = 1. - (QSQc(1)/QSQMAX)**RR QSQ = QSQc(1)/((1.-FACQSQ*userran(iseed))**(1./RR)) PFQSQ = RR*QSQc(1)**RR/FACQSQ/(QSQ**(RR+1.)) C endif C C Bjorken_y variable [ y = 2.*amp*vu/(s-ampsq-m_e^2) ] C y = qsq/xbj/(s-ampsq-eb(5)) C C Generate the scattered electron, differently if initial state C radiation occurs or not C if(radcor_i) then C C Transform to CMS of electron and proton. Proton along +z' axis C DO 10 I=1,5 PA(1,I) = PB(I) PA(2,I) = EB(I) PA(3,I) = 0. PA(4,I) = 0. 10 CONTINUE DO 20 J=1,3 20 DBETA(1,J)=DBLE(PA(1,J)+PA(2,J))/DBLE(PA(1,4)+PA(2,4)) CALL DBOOST2(-DBETA(1,1),-DBETA(1,2),-DBETA(1,3)) DPHI(1)=ULANGL(PA(1,1),PA(1,2)) CALL ROTATE2(0.D+0,-DPHI(1)) DTHETA(1)=ULANGL(PA(1,3),PA(1,1)) CALL ROTATE2(-DTHETA(1),0.D+0) C C Energy and theta of the scattered electron (from Q2 and y) C taking into account the virtuality of the incoming electron C escat = (1.-y)*(pa(1,4)*pa(2,4)+pa(1,3)*pa(1,3)) + + (qsq+pa(2,5)+amesq)/2. escat = escat/(pa(1,4)+pa(2,4)) costhee = pa(2,4)*escat - (qsq+pa(2,5)+amesq)/2. efmod = escat*escat-amesq if(efmod.lt.0.) then write(6,*)'WARNING-EVENT2: efmod < 0 ',ievent,efmod nmass = nmass+1 goto 1001 endif costhee = costhee/pa(1,3)/sqrt(efmod) if(costhee.gt.1.) then if(costhee.lt.1.005) then costhee = 1. else nmass = nmass+1 goto 1001 endif endif sinthee = sqrt(1.-costhee*costhee) PHE = 2.*PI*userran(iseed) PFPHIE = 1.0 C C Fill vector of scattered electron C PA(3,1) = ESCAT*SINTHEE*COS(PHE) PA(3,2) = ESCAT*SINTHEE*SIN(PHE) PA(3,3) = -ESCAT*COSTHEE PA(3,4) = ESCAT PA(3,5) = amesq C C Back to the Zeus reference system C CALL ROTATE2(DTHETA(1),0.D+0) CALL ROTATE2(0.D+0,DPHI(1)) CALL DBOOST2(DBETA(1,1),DBETA(1,2),DBETA(1,3)) C C Fill vector of the scattered electron. Theta and phi angles C DO 30 I=1,5 30 ES(I) = PA(3,I) PHIE = ULANGL(ES(1),ES(2)) THETAE = ATAN2(SQRT(ES(1)*ES(1)+ES(2)*ES(2)),ES(3)) thetae = pi-thetae C C Check if the initial state radiation is absorbed by the C scattered electron (interference term) C itype = 2 call radcorr(itype) if(radcor_a) then C C Recalculate angles and momentum of the scattered electron C do 40 i=1,5 40 eb(i) = eb_i(i) s = stot(5) y = qsq/xbj/(s-ampsq-eb(5)) C ESCAT = (QSQ+2*AMESQ)/(4.*EB(4)) + (1.-Y)*EB(4) THETAE2 = (QSQ+2*AMESQ)/(4.*EB(4)*ESCAT) IF(THETAE2.GT.1.)THEN THETAE = PI/2. ELSE THETAE = 2.*ASIN(SQRT(THETAE2)) ENDIF PHIE = 2.*PI*userran(iseed) PFPHIE = 1.0 C C Fill vector of scattered electron C ES(1) = ESCAT*SIN(THETAE)*COS(PHIE) ES(2) = ESCAT*SIN(THETAE)*SIN(PHIE) ES(3) = -ESCAT*COS(THETAE) ES(4) = ESCAT ES(5) = amesq C endif C else C C Angles and momentum of scattered electron C ESCAT = (QSQ+2*AMESQ)/(4.*EB(4)) + (1.-Y)*EB(4) THETAE2 = (QSQ+2*AMESQ)/(4.*EB(4)*ESCAT) IF(THETAE2.GT.1.)THEN THETAE = PI/2. ELSE THETAE = 2.*ASIN(SQRT(THETAE2)) ENDIF PHIE = 2.*PI*userran(iseed) PFPHIE = 1.0 C C Fill vector of scattered electron C ES(1) = ESCAT*SIN(THETAE)*COS(PHIE) ES(2) = ESCAT*SIN(THETAE)*SIN(PHIE) ES(3) = -ESCAT*COS(THETAE) ES(4) = ESCAT ES(5) = amesq C if(.not.radcor) goto 100 C C Check if final state radiation occurs in this event C itype = 1 call radcorr(itype) if(radcor_f) then C C Recalculate the mass of the scattered, off-mass shell, electron C efmod = escat*escat*(1.-zph)*(1.-zph) - amesq if(efmod.lt.0.) then nmass = nmass+1 goto 1001 endif efmod = sqrt(efmod) if((escat*(1.-zph)).lt.0.2) then am2prime = amesq + + 2.*zph*escat*(escat*(1.-zph)-efmod*cos(theph)) else am2prime = amesq + 2.*zph*(1.-zph)*escat*sin(theph/2.)**2 + + zph*amesq*cos(theph)/(1.-zph) endif if(am2prime.lt.amesq) am2prime = amesq C C Redefine es vector C es(5) = am2prime esmod = escat*escat - es(5) if(esmod.lt.0.) then nmass = nmass+1 goto 1001 endif esmod = sqrt(esmod) es(1) = esmod*sin(thetae)*cos(phie) es(2) = esmod*sin(thetae)*sin(phie) es(3) = -esmod*cos(thetae) es(4) = escat C C Rescale Q2 and y in accordance with new mass am2prime C costhee = cos(thetae) qsq = 2*eb(4)*es(4)-2*abs(eb(3))*esmod*costhee-amesq-am2prime y = 1. - (pb(4)*es(4)+pb(3)*esmod*costhee)/ + (pb(4)*eb(4)+pb(3)*abs(eb(3))) C C Fill vector of the final, on-mass shell, electron C call dif5(es,phot,es_f) es_f(5) = amesq PHIE = ULANGL(ES_f(1),ES_f(2)) THETAE = ATAN2(SQRT(ES_f(1)*ES_f(1)+ES_f(2)*ES_f(2)),ES_f(3)) thetae = pi-thetae C C Recalculate xbj C xbj = qsq/y/(s-ampsq-eb(5)) C endif C 100 continue endif C C Nu variable [ vu = p.q/amp ] C VU = QSQ/(2.*AMP*XBJ) C C Mass hadronic final state [ wsq = ampsq + qsq*(1.-xbj)/xbj ] C wsq = (2.0*amp*vu)+ampsq-qsq C C Gamt = #photons / dxdQ2 (Hand_convention) C adu = qsq*ampsq/(s*s) one_eps = (y*y/2.+2.*adu)/(1.0-y+(y*y/2.0)+adu) aak = vu - qsq/(2.*amp) gamt = 4.*alpha*ampsq*aak*vu/(2.*pi*qsq*s*s*xbj*one_eps) C C Virtual photon vector. GAM(5) should be -QSQ generated above. C CALL DIF5(EB,ES,GAM) C C Compute quantities in photon-proton reference system. C CALL GAMPRO(Jmass) C IF(JMASS.NE.0) goto 1000 C C Fill vector of scattered proton C DO 50 I=1,5 50 PS(I) = PA(3,I) C C x_L = |p'|/|p| : standard LPS definition C PSMOD = SQRT(PS(1)*PS(1)+PS(2)*PS(2)+PS(3)*PS(3)) if(pbeam.ne.0.) XL = PSMOD/PBEAM C C Transverse momentun of scattered proton C PT = SQRT(PS(1)*PS(1)+PS(2)*PS(2)) C C Theta and phi angles between p and p' C thetap = atan2(pt,ps(3)) PHIP = ATAN2(PS(2),PS(1)) IF(PS(2).LT.0.) PHIP = PHIP + 2*PI IF((PS(2).EQ.0.).AND.PS(1).LT.0.) PHIP = PI C C Pomeron vector. POM(5) should be T generated above. C CALL DIF5(PB,PS,POM) C C Hadronic final state X. XVECT(5) should be AMX2 generated above. C CALL ADD5(GAM,POM,XVECT) C C Photon - Pomeron interaction ----> q - qbar C CALL GAMPOM2(Kmass) C IF(KMASS.NE.0) goto 1000 C C Fill vectors of q and qbar. Theta and phi angles C DO I=1,5 Q(I) = PA(3,I) QBAR(I) = PA(4,I) ENDDO C PHIQ = ULANGL(Q(1),Q(2)) THEQ = ATAN2(SQRT(Q(1)*Q(1)+Q(2)*Q(2)),Q(3)) PHIQB = ULANGL(QBAR(1),QBAR(2)) THEQB = ATAN2(SQRT(QBAR(1)*QBAR(1)+QBAR(2)*QBAR(2)),QBAR(3)) C C *** dsigma(gamma*+p->p+X) / dMx2 dt dk_t^2 *** C CALL SMXSQ2(WTMX2) CALL SIGT2(WTT) CALL SIGTHJ2(WTTHJ) C C Total phase space and total event weight *** in nanobarn *** C WPHASE = PFXBJ*PFQSQ*PFT*PFMX2*PFTHJC c c K-factor with running alpha_s c if(ikfact.eq.1) then alps = alphas(amx2) fk_2 = exp(alps*pi*4./3.) elseif(ikfact.eq.2) then alps = alphas(amx2/4) fk_2 = exp(alps*pi*4./3.) elseif(ikfact.eq.3) then alps = alphas(amx2/4) fk_2 = (1.+alps*pi*4./6.)**2 endif C cs2(1) = 1000.*fk_2*gamt*wtmx2*wtt*wtthj/wphase cs2(2) = 1000.*fk_2*wtmx2*wtt*wtthj/(pfmx2*pft*pfthjc) c cst2(1) = 1000.*fk_2*gamt*sig2t*wtt*wtthj/wphase cst2(2) = 1000.*fk_2*sig2t*wtt*wtthj/(pfmx2*pft*pfthjc) c csl2(1) = 1000.*fk_2*gamt*sig2l*wtt*wtthj/wphase csl2(2) = 1000.*fk_2*sig2l*wtt*wtthj/(pfmx2*pft*pfthjc) C sig2 = cs2(1) C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE CROSS2(IAGAIN) C ------------------------------------------------------------------- C Evaluation of the cross section for a q-qbar final state C Called by : JETNUM, RIDIEVT C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' include 'ridkin.inc' C real clara,bella c IAGAIN = 0 c c Transverse and longitudinal cross sections *** in nanobarn *** c do i = 1,2 tcst2(i) = tcst2(i) + cst2(i) ! sig_T(e+p->e+p+jj) tcste2(i) = tcste2(i) + cst2(i)*cst2(i) ! sig_T(gamma*+p->p+jj) enddo c do i = 1,2 tcsl2(i) = tcsl2(i) + csl2(i) ! sig_L(e+p->e+p+jj) tcsle2(i) = tcsle2(i) + csl2(i)*csl2(i) ! sig_L(gamma*+p->p+jj) enddo c weit2 = cst2(1) weil2 = csl2(1) C if(.not.wei) then c c Use weighted events (all) to calculate the gamma*-p c cross section *** in nanobarn *** c NOTE: the gamma*-p cross section is evaluated at the mean c values of Q^2 and W (xbj) in the kinematic c region of generation. Since these two variables are c not integrated over, the phase space does not c compensate for different functions of generation. c Different functions determine different mean values of c Q^2 and W and hence different values of cross sections. c tcs2(2) = tcs2(2) + cs2(2) tcse2(2) = tcse2(2) + cs2(2)**2 ! sig(gamma*+p->p+jj) c c Count all the generated events c itoteve2 = itoteve2+1 C C To produce unweighted events check that Rndm < CS(1)/WEIMAX2 C clara = userran(iseed) if(clara.gt.cs2(1)/weimax2) then iagain = 1 return endif C C Total e-p cross section *** in nanobarn *** C tcs2(1) = tcs2(1) + weimax2 ! sig(e+p->e+p+jj) tcse2(1) = tcse2(1) + weimax2**2 C C Cross sections for different classes of radiative events C if(radcor_i) then if(.not.radcor_a) then tcsi2 = tcsi2 + weimax2 tcsie2 = tcsie2 + weimax2**2 else tcsa2 = tcsa2 + weimax2 tcsae2 = tcsae2 + weimax2**2 endif elseif(radcor_f) then tcsf2 = tcsf2 + weimax2 tcsfe2 = tcsfe2 + weimax2**2 else tcsb2 = tcsb2 + weimax2 tcsbe2 = tcsbe2 + weimax2**2 endif C C If weight.gt.weimax count the event two times in C accordance with probability C if(cs2(1).gt.weimax2) then iwei2 = iwei2+1 bella = userran(iseed) if(bella.le.(cs2(1)/weimax2 - 1.)) then itoteve2 = itoteve2+1 itwoev2 = itwoev2+1 two = .true. tcs2(1) = tcs2(1) + weimax2 tcse2(1) = tcse2(1) + weimax2**2 if(radcor_i) then if(.not.radcor_a) then tcsi2 = tcsi2 + weimax2 tcsie2 = tcsie2 + weimax2**2 else tcsa2 = tcsa2 + weimax2 tcsae2 = tcsae2 + weimax2**2 endif elseif(radcor_f) then tcsf2 = tcsf2 + weimax2 tcsfe2 = tcsfe2 + weimax2**2 else tcsb2 = tcsb2 + weimax2 tcsbe2 = tcsbe2 + weimax2**2 endif endif endif C WEIGHT = WEIMAX2 C else C C Total e-p cross section *** in nanobarn *** C if(cs2(1).eq.0.) then nmass = nmass+1 iagain = 1 return endif do i = 1,2 tcs2(i) = tcs2(i) + cs2(i) ! sig(e+p->e+p+jj) tcse2(i) = tcse2(i) + cs2(i)*cs2(i) ! sig(gamma*+p->p+jj) enddo C C Cross sections for different classes of radiative events C if(radcor_i) then if(.not.radcor_a) then tcsi2 = tcsi2 + cs2(1) tcsie2 = tcsie2 + cs2(1)*cs2(1) else tcsa2 = tcsa2 + cs2(1) tcsae2 = tcsae2 + cs2(1)*cs2(1) endif elseif(radcor_f) then tcsf2 = tcsf2 + cs2(1) tcsfe2 = tcsfe2 + cs2(1)*cs2(1) else tcsb2 = tcsb2 + cs2(1) tcsbe2 = tcsbe2 + cs2(1)*cs2(1) endif C if(ievent.ne.0) WEIGHT = CS2(1)/twojets C C Check phase space of generated variables C if(ievent.ne.0) CALL GENERA2 C endif C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE GAMPRO(Imass) C ------------------------------------------------------------------- C To compute quantities in photon-proton reference system C Called by : EVENT2 C C Author : Ada Solano C ------------------------------------------------------------------- C include 'rboost.inc' include 'ridseq.inc' include 'ridkin.inc' C real p2_cms,p_cms,pz_cms,pt2_cms,pt_cms,pin_cms c imass = 0 C C Transform to CMS of photon and proton. Proton along +z' axis. C DO 1 I=1,5 PA(1,I) = PB(I) PA(2,I) = GAM(I) PA(3,I) = 0. PA(4,I) = 0. 1 CONTINUE DO 10 J=1,3 10 DBETA(1,J)=DBLE(PA(1,J)+PA(2,J))/DBLE(PA(1,4)+PA(2,4)) CALL DBOOST2(-DBETA(1,1),-DBETA(1,2),-DBETA(1,3)) DPHI(1)=ULANGL(PA(1,1),PA(1,2)) CALL ROTATE2(0.D+0,-DPHI(1)) DTHETA(1)=ULANGL(PA(1,3),PA(1,1)) CALL ROTATE2(-DTHETA(1),0.D+0) C C Generate t=(p-p')**2 like exp(b*t) : -----|---------|------|---> C tmin tmax 0 t C xpmin = (amx2c(1)+qsq)/(wsq+qsq-ampsq) if(xpmin.gt.xpomc(2)) then imass = 1 nmass = nmass + 1 return endif ttmax = -xpmin*xpmin*ampsq/(1.-xpmin) ttmax = min(ttmax,tmax) if(ttmax.le.tmin) then imass = 1 nmass = nmass + 1 return endif FACT = EXP(BSLOPE*ttmax) - EXP(BSLOPE*TMIN) T = ALOG(EXP(BSLOPE*ttmax)-FACT*userran(iseed))/BSLOPE PFT = -BSLOPE*EXP(BSLOPE*T)/FACT PFT = ABS(PFT) c c t = 0. C C M_X**2 = (1.-xl_cms)*(W**2+Q**2-Mp**2)+t-Q**2 C = (1.-xl_cms)*(Q**2/x)+t-Q**2 C xpom_max = (1.-xl_cms)_max = 0.1 : definition of diffractive event C Since beta<=1 --> xpom_min>=xbj C xpmax = (t+sqrt(t*t-4*ampsq*t))/ampsq/2. xpmax = min(xpmax,xpomc(2)) AMX2MN = max(xpommn*(WSQ+QSQ-AMPSQ)+T-QSQ,amx2c(1)) AMX2MX = min(xpmax*(WSQ+QSQ-AMPSQ)+T-QSQ,amx2c(2)) if(amx2mx.le.amx2mn) then imass = 1 nmass = nmass + 1 return endif C C Generate M_X**2 C if(wei) then ! like 1/(M_X**2+Q**2) c FACMX2 = ALOG((AMX2MX+QSQ)/(AMX2MN+QSQ)) AMX2 = (AMX2MN+QSQ)*EXP(FACMX2*userran(iseed)) - QSQ PFMX2 = 1./FACMX2/(AMX2+QSQ) C else ! like 1/(M_X**2+Q**2)^(SS+1) C FACMX2 = 1. - ((AMX2MN+QSQ)/(AMX2MX+QSQ))**SS AMX2 = (AMX2MN+QSQ)/((1.-FACMX2*userran(iseed))**(1./SS))-QSQ PFMX2 = SS*(AMX2MN+QSQ)**SS/FACMX2/((AMX2+QSQ)**(SS+1.)) c endif c c x_pom = fraction of proton momentum carried out by the pomeron c ( it should be equal to 1-x_L ) c xpom = (amx2+qsq-t)/(wsq-ampsq+qsq) if(xpom.lt.xpomlim(1)) xpomlim(1) = xpom if(xpom.gt.xpomlim(2)) xpomlim(2) = xpom C C Use energy and momentum conservation to calculate the C momentum of the scattered proton: in the c.m.s. C E_p' + E_X = W and |p'| = |p_X| C p2_cms = ( (ampsq+amx2-wsq)**2 - 4*ampsq*amx2 )/4./wsq p_cms = sqrt(p2_cms) C C Use t and |p'| to calculate the transverse momentum C of the scattered proton C pin_cms = sqrt(pa(1,4)**2-ampsq) pz_cms = ( t-2*ampsq+2*pa(1,4)*sqrt(p2_cms+ampsq) )/2./pin_cms pt2_cms = p2_cms-pz_cms**2 pt2cms = 0. c c Only to avoid problems from rounding errors (checked) c if(pt2_cms.lt.0) then nptcms = nptcms+1 pt2cms = pt2_cms pt2_cms = 0.000001 endif pt_cms = sqrt(pt2_cms) C xl_cms = p_cms/pin_cms C C Generate Phi angle flat C PHIPCM = 2.*PI*userran(iseed) PFPHIP = 1.0 C C Fill vectors of scattered proton ( photon-proton cms ) C PA(3,1) = PT_CMS*COS(PHIPCM) PA(3,2) = PT_CMS*SIN(PHIPCM) c pa(3,3) = pz_cms pa(3,4) = sqrt(p2_cms+ampsq) c PA(3,5) = PB(5) C C Back to the Zeus reference system C CALL ROTATE2(DTHETA(1),0.D+0) CALL ROTATE2(0.D+0,DPHI(1)) CALL DBOOST2(DBETA(1,1),DBETA(1,2),DBETA(1,3)) C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE GAMPOM2(Imass) C ------------------------------------------------------------------- C To compute quantities in gamma-pomeron reference system C Called by : EVENT2 C C Author : Ada Solano C ------------------------------------------------------------------- C include 'rboost.inc' include 'ridseq.inc' include 'ridkin.inc' C REAL PIPPO,PLUTO,MINNI REAL THJCM2,THEOSQ,THEO,THJLIM,ABSTHJ,COMTHJ,ALPTH,ATH C imass = 0 C C Transform to CMS of photon and pomeron. Pomeron along +z' axis. C DO 1 I=1,5 PA(1,I) = POM(I) PA(2,I) = GAM(I) PA(3,I) = 0. PA(4,I) = 0. 1 CONTINUE DO 10 J=1,3 10 DBETA(1,J)=DBLE(PA(1,J)+PA(2,J))/DBLE(PA(1,4)+PA(2,4)) CALL DBOOST2(-DBETA(1,1),-DBETA(1,2),-DBETA(1,3)) DPHI(1)=ULANGL(PA(1,1),PA(1,2)) CALL ROTATE2(0.D+0,-DPHI(1)) DTHETA(1)=ULANGL(PA(1,3),PA(1,1)) CALL ROTATE2(-DTHETA(1),0.D+0) C C Generate theta distribution of the 2 jets like C 1/((theta_jet**2+thetao**2/2)**2) C THEOSQ = 1./(XVECT(5)/4.) THEO = SQRT(THEOSQ) THJLIM = CTH*THEO ATH = EXP(-CTH2)*(SIN(CTH*THEO)**4)/(THEOSQ*THEOSQ) if(.not.thlim) then FACTHJ = PI*PI/(2.*THEOSQ*(THEOSQ+PI*PI/2.)) COMTHJ = 1./THEOSQ - FACTHJ*userran(iseed) else alpth = 1./(2.*thjlim*thjlim+theosq) facthj = alpth - 2./(pi*pi+2.*theosq) comthj = alpth - facthj*userran(iseed) endif THJCM2 = (1./COMTHJ - THEOSQ)/2. ABSTHJ = SQRT(THJCM2) PIPPO = userran(iseed) IF(PIPPO.LT.0.5)THEN THJCM = ABSTHJ ELSE THJCM = PI - ABSTHJ ENDIF PFTHJC = 2./FACTHJ/((2.*THJCM2 + THEOSQ)**2) ! theta**2 PFTHJC = 2.*absthj*pfthjc ! |theta| pfthjc = pfthjc/2. ! theta C PHJCM = 2.*PI*userran(iseed) PFPHJC = 1. C EJET = SQRT(XVECT(5))/2. C C Use all the four flavours for Ryskin's cross section: C Flavour of the initiators ( Lund ): 1 = down, 2 = up C 3 = strange, 4 = charm C C Check if charm or strange production is possible C if(onlych) then if(ejet.gt.cthre) then charm = .true. iflav = 4 goto 100 else imass = 1 nmass = nmass + 1 return endif endif c if(.not.lightf) then if(ejet.gt.cthre) then charm = .true. else if(ejet.gt.sthre) strange = .true. endif endif c IFLAV = 1 PLUTO = userran(iseed) IF(PLUTO.GT.0.2) IFLAV = IFLAV+1 if(.not.lightf) then MINNI = userran(iseed) if(charm) then IF(MINNI.GT.0.5) IFLAV = IFLAV+2 elseif(strange) then if(minni.lt.1./6.) iflav = 3 endif endif C C Fill vectors of the 2 jets ( photon-pomeron cms ) C The mass of the jets is NOT taken to be zero C 100 continue PA(3,5) = ULMASS(IFLAV)**2 PA(3,4) = EJET PJET = SQRT( EJET*EJET-PA(3,5) ) PA(3,1) = PJET*SIN(THJCM)*COS(PHJCM) PA(3,2) = PJET*SIN(THJCM)*SIN(PHJCM) PA(3,3) = PJET*COS(THJCM) PA(4,1) = -PA(3,1) PA(4,2) = -PA(3,2) PA(4,3) = -PA(3,3) PA(4,4) = PA(3,4) PA(4,5) = PA(3,5) C C Back to the Zeus reference system (for the four vectors C containing the photon,the pomeron and the two partons). C CALL ROTATE2(DTHETA(1),0.D+0) CALL ROTATE2(0.D+0,DPHI(1)) CALL DBOOST2(DBETA(1,1),DBETA(1,2),DBETA(1,3)) C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE FRAGM2 C ------------------------------------------------------------------- C Parton shower + hadronization of the q-qbar system (JETSET) C Called by : RIDIEVT C C Author : Ada Solano C ------------------------------------------------------------------- C include 'rboost.inc' include 'ridseq.inc' include 'ridkin.inc' C COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) EXTERNAL LUDATA C integer id,ide,idp,ijoin2(2) REAL QMAX,angphi,angthe C MSTU(1) = 0 MSTU(2) = 0 MSTU(3) = 0 C C Initial and scattered beam particles + radiated photons C id = 1 ide = id call lu1ent(id,-11*LFLAG,eb_i(4),paru(1),0.) k(id,1) = 21 k(id,2) = -11*LFLAG k(id,3) = 0 do 1 i=1,4 1 p(id,i) = eb_i(i) p(id,5) = ame c id = id+1 idp = id call lu1ent(id,2212,pb(4),0.,0.) k(id,1) = 21 k(id,2) = 2212 k(id,3) = 0 do 2 i=1,4 2 p(id,i) = pb(i) p(id,5) = amp c if(radcor_i.and.(.not.radcor_a)) then id = id+1 angphi = ULANGL(phot(1),phot(2)) angthe = ATAN2(SQRT(phot(1)*phot(1)+phot(2)*phot(2)),phot(3)) call lu1ent(id,22,phot(4),angthe,angphi) k(id,1) = 21 k(id,2) = 22 k(id,3) = ide do 3 i=1,4 3 p(id,i) = phot(i) p(id,5) = 0. c id = id+1 angphi = ULANGL(eb(1),eb(2)) angthe = ATAN2(SQRT(eb(1)*eb(1)+eb(2)*eb(2)),eb(3)) call lu1ent(id,-11*LFLAG,eb(4),angthe,angphi) k(id,1) = 21 k(id,2) = -11*LFLAG k(id,3) = ide do 4 i=1,4 4 p(id,i) = eb(i) p(id,5) = -sqrt(abs(eb(5))) ide = id endif c id = id+1 angphi = ULANGL(es(1),es(2)) angthe = ATAN2(SQRT(es(1)*es(1)+es(2)*es(2)),es(3)) call lu1ent(id,-11*LFLAG,es(4),angthe,angphi) k(id,1) = 21 k(id,2) = -11*LFLAG k(id,3) = ide do 5 i=1,4 5 p(id,i) = es(i) p(id,5) = sqrt(es(5)) ide = id c if(radcor_f) then id = id+1 angphi = ULANGL(phot(1),phot(2)) angthe = ATAN2(SQRT(phot(1)*phot(1)+phot(2)*phot(2)),phot(3)) call lu1ent(id,22,phot(4),angthe,angphi) k(id,1) = 21 k(id,2) = 22 k(id,3) = ide do 6 i=1,4 6 p(id,i) = phot(i) p(id,5) = 0. c id = id+1 angphi = ULANGL(es_f(1),es_f(2)) angthe = ATAN2(SQRT(es_f(1)*es_f(1)+es_f(2)*es_f(2)),es_f(3)) call lu1ent(id,-11*LFLAG,es_f(4),angthe,angphi) k(id,1) = 21 k(id,2) = -11*LFLAG k(id,3) = ide do 7 i=1,4 7 p(id,i) = es_f(i) p(id,5) = ame endif c id = id+1 call lu1ent(id,2212,ps(4),thetap,phip) k(id,1) = 21 k(id,2) = 2212 k(id,3) = idp do 8 i=1,4 8 p(id,i) = ps(i) p(id,5) = amp c do 9 i=1,id k(i,4) = 0 k(i,5) = 0 9 continue C C The fragmentation is done in the photon-pomeron cms C id = id+1 idq = id CALL LU1ENT(-id,IFLAV,EJET,THJCM,PHJCM) CALL LU1ENT(id+1,-IFLAV,EJET,PARU(1)-THJCM,PHJCM+PARU(1)) C C Connect the colour C ijoin2(1) = id ijoin2(2) = id+1 CALL LUJOIN(2,IJOIN2) C C Parton shower ( final state gluon radiation ) C if(mstj(41).gt.0)then QMAX = plu(id,12) CALL LUSHOW(id,id+1,QMAX) DO 10 J=1,N IF(K(J,1).EQ.3)CALL LUPREP(J) 10 CONTINUE endif C CALL LUEXEC if(mstu(28).ne.0) then write(6,*)'JETSET Warning: error code, ievent, iflav, ejet' write(6,*)mstu(28),ievent,iflav,ejet else if(iflav.eq.3) then if(ejet.lt.stmin) stmin = ejet elseif(iflav.eq.4) then if(ejet.lt.chmin) chmin = ejet endif endif C C Back to the Zeus reference system ( for all the entries of C the COMMON/LUJETS/ ) C theta(1) = sngl(dtheta(1)) phi(1) = sngl(dphi(1)) DO 100 I=1,3 100 BETA(1,I) = SNGL(DBETA(1,I)) c c Avoid rotation of initial and scattered beam particles c mstu(1) = id CALL LUROBO(theta(1),phi(1),BETA(1,1),BETA(1,2),BETA(1,3)) mstu(1) = 0 c c Scattered electron and proton + radiatied photon as final states c k(id-1,1) = 1 k(id-2,1) = 1 if(radcor_i.and.(.not.radcor_a)) k(id-4,1) = 1 if(radcor_f) k(id-3,1) = 1 c c Not in ZDIS if(ievent.le.nevpri) call lulist(1) c RETURN END C C ------------------------------------------------------------------- C SUBROUTINE SMXSQ2(SIG) C ------------------------------------------------------------------- C To compute the differential cross section in Mx**2 C Called by : EVENT2 C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' include 'ridkin.inc' C real z,amq2,q02,q002,amqkt2,q2l,aal C C dsigma_DD(gamma*+p->p+X) / dMx2 ( in microbarn! ) C z = xbj/xpom akt2 = (pjet*sin(thjcm))**2 if(akt2.lt.akt2c(1)) then c c Avoid to small k_t^2. The event is not taken. c laacheck = laacheck+1 sig2t = 0. sig2l = 0. goto 10 endif c if(akt2.lt.akt2lim) akt2lim = akt2 c amq2 = effmas(iflav)**2 amqkt2 = akt2+amq2 q02 = amqkt2/(1-z) c c Only to avoid log(number<1)<0 !! c q002 = q02 alps = alphas(q002) if(alps.gt.alpscut) alps = alpscut alps2 = alps*alps c c Note that the factor 1/2 wrt the paper comes from the generation c a` la Kolya of the theta distribution for the q-qbar pair (1 dec '92) c if(lightf) then sig = (1./bslope)*(xpom*GN(xpom,q02))**2* + (5./9.)*(2.*pi*pi*alpha/qsq/3.)*hbarc2*cmcbar*alps2* + (1.-z)**2*(z*z)* + qsq/((amx2+qsq-t)**2) elseif(onlych) then sig = (1./bslope)*(xpom*GN(xpom,q02))**2* + (4./9.)*(2.*pi*pi*alpha/qsq/3.)*hbarc2*cmcbar*alps2* + (1.-z)**2*(z*z)* + qsq/((amx2+qsq-t)**2) else if(charm) then sig = (1./bslope)*(xpom*GN(xpom,q02))**2* + (10./9.)*(2.*pi*pi*alpha/qsq/3.)*hbarc2*cmcbar*alps2* + (1.-z)**2*(z*z)* + qsq/((amx2+qsq-t)**2) elseif(strange) then sig = (1./bslope)*(xpom*GN(xpom,q02))**2* + (6./9.)*(2.*pi*pi*alpha/qsq/3.)*hbarc2*cmcbar*alps2* + (1.-z)**2*(z*z)* + qsq/((amx2+qsq-t)**2) else sig = (1./bslope)*(xpom*GN(xpom,q02))**2* + (5./9.)*(2.*pi*pi*alpha/qsq/3.)*hbarc2*cmcbar*alps2* + (1.-z)**2*(z*z)* + qsq/((amx2+qsq-t)**2) endif endif c c Longitudinal cross section [Ryskin*Wusthoff(7/2/96)*Levin(24/4/96)] c q2l = amqkt2/qsq aal = q2l*z/(1.-z) c c Note that since aal=(k_t^2+m_q^2)/M_x^2 we need to check c that (k_t^2+m_q^2) q - qbar C CALL GAMGLU(Jskip) IF(JSKIP.EQ.1) GOTO 1000 C C Fill vectors of q and qbar. Theta and phi angles C DO 30 I=1,5 Q1(I) = PA(3,I) Q2(I) = PA(4,I) 30 CONTINUE C C q1=q or q1=qbar ? C MINNI = userran(iseed) IF(MINNI.GT.0.5) THEN DO 40 I=1,5 QBAR(I) = Q1(I) Q(I) = Q2(I) 40 CONTINUE ELSE DO 50 I=1,5 Q(I) = Q1(I) QBAR(I) = Q2(I) 50 CONTINUE ENDIF C PHIQ = ULANGL(Q(1),Q(2)) THEQ = ATAN2(SQRT(Q(1)*Q(1)+Q(2)*Q(2)),Q(3)) PHIQB = ULANGL(QBAR(1),QBAR(2)) THEQB = ATAN2(SQRT(QBAR(1)*QBAR(1)+QBAR(2)*QBAR(2)),QBAR(3)) C C Total phase space and total event weight ( in nanobarn! ) C WPHASE = PFXF*PFPT2*PFKT2*PFQSQ*PFEXBJ*PFZ*PFTHJC c c k-factor with running alpha_s c if(ikfact.eq.1) then alps = alphas(amx2) fk_3 = exp(alps*pi*3.) elseif(ikfact.eq.2) then alps = alphas(amx2/4) c fk_3 = exp(alps*pi*3.) fk_3 = exp(alps*pi*3.)/3.+exp(alps*pi*4./3.)*2./3. elseif(ikfact.eq.3) then alps = alphas(amx2/4) fk_3 = (1.+alps*pi*3./2.)**2 endif C if(transv) then cst3(1) = fk_3*GAMT*WTXF*WTPT2*WTKT2*WTZ*WTQINT*CM2DIM/WPHASE cst3(2) = fk_3*WTXF*WTPT2*WTKT2*WTZ*WTQINT*CM2DIM / + (PFXF*PFPT2*PFKT2*PFZ*PFTHJC) else cst3(1) = 0. cst3(2) = 0. endif c if(longit) then c c R=sigma_L/sigma_T from Levin et al., hep-ph/9606443 (14.2.97) c erre = effmas(iflav)**2/qsq bet2 = 1.-4*erre*z/(1.-z) if(bet2.ge.0.) then bet = sqrt(bet2) else write(6,*)'WARNING-EVENT3: bet2<0 ',ievent,bet2 csl3(1) = 0. csl3(2) = 0. goto 999 endif aln = alog((1.+bet)*(1.+bet)) - alog(4*erre*z/(1.-z)) cl = 4*bet*z*(1.-z)-8*z*z*erre*aln ct = aln*(z*z+(1.-z)*(1.-z)+4*z*(1.-3*z)*erre+ + 8*z*z*(erre-erre*erre)) + bet*(4*z*(1.-z)*(1.-erre)-1.) csl3(1) = cst3(1)*cl/ct csl3(2) = cst3(2)*cl/ct else csl3(1) = 0. csl3(2) = 0. endif c 999 epsilon = 2*(1.-y)/(2.-2*y+y*y) cs3(1) = cst3(1)+epsilon*csl3(1) cs3(2) = cst3(2)+epsilon*csl3(2) C sig3 = cs3(1) C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE CROSS3(IAGAIN) C ------------------------------------------------------------------- C Evaluation of the cross section for a q-qbar-gluon final state C Called by : JETNUM, RIDIEVT C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' include 'ridkin.inc' C real clara,bella c IAGAIN = 0 c c Transverse and logitudinal cross sections *** in nanobarn *** c do i = 1,2 tcst3(i) = tcst3(i) + cst3(i) ! sig_T(e+p->e+p+3jets) tcste3(i) = tcste3(i) + cst3(i)*cst3(i) ! sig_T(gamma*+p->p+3jets) enddo c do i = 1,2 tcsl3(i) = tcsl3(i) + csl3(i) ! sig_L(e+p->e+p+3jets) tcsle3(i) = tcsle3(i) + csl3(i)*csl3(i) ! sig_L(gamma*+p->p+3jets) enddo c weit3 = cst3(1) weil3 = csl3(1) c if(.not.wei) then c c Use weighted events (all) to calculate the gamma*-p cross section c NOTE: the gamma*-p cross section is evaluated at the mean c values of Q^2 and W (xbj or vu) in the kinematic c region of generation. Since these two variables are c not integrated over, the phase space does not c compensate for different functions of generation. c Different functions determine different mean values of c Q^2 and W and hence different values of cross sections. c tcs3(2) = tcs3(2) + cs3(2) ! sig(gamma*+p -> p+3jets) tcse3(2) = tcse3(2) + cs3(2)*cs3(2) c c Count all the generated events c ITOTEVE3 = ITOTEVE3+1 C C To produce unweighted events check that Rndm < cs(1)/WEIMAX3 C clara = userran(iseed) if(clara.gt.cs3(1)/weimax3) then iagain = 1 return endif C C Total e-p cross section *** in nanobarn *** C tcs3(1) = tcs3(1) + WEIMAX3 ! sig(e+p -> e+p+3jets) tcse3(1) = tcse3(1) + WEIMAX3**2 C C Cross sections for different classes of radiative events C if(radcor_i) then if(.not.radcor_a) then tcsi3 = tcsi3 + weimax3 tcsie3 = tcsie3 + weimax3**2 else tcsa3 = tcsa3 + weimax3 tcsae3 = tcsae3 + weimax3**2 endif elseif(radcor_f) then tcsf3 = tcsf3 + weimax3 tcsfe3 = tcsfe3 + weimax3**2 else tcsb3 = tcsb3 + weimax3 tcsbe3 = tcsbe3 + weimax3**2 endif C C If weight.gt.weimax count the event two times in C accordance with probability C if(cs3(1).gt.weimax3) then iwei3 = iwei3+1 bella = userran(iseed) if(bella.le.(cs3(1)/weimax3 - 1.)) then itoteve3 = itoteve3+1 itwoev3 = itwoev3+1 two = .true. tcs3(1) = tcs3(1) + WEIMAX3 tcse3(1) = tcse3(1) + WEIMAX3**2 if(radcor_i) then if(.not.radcor_a) then tcsi3 = tcsi3 + weimax3 tcsie3 = tcsie3 + weimax3**2 else tcsa3 = tcsa3 + weimax3 tcsae3 = tcsae3 + weimax3**2 endif elseif(radcor_f) then tcsf3 = tcsf3 + weimax3 tcsfe3 = tcsfe3 + weimax3**2 else tcsb3 = tcsb3 + weimax3 tcsbe3 = tcsbe3 + weimax3**2 endif endif endif c WEIGHT = WEIMAX3 C else C C Total e-p cross section *** in nanobarn *** C if(cs3(1).eq.0.) then nskip = nskip+1 iagain = 1 return endif do i=1,2 tcs3(i) = tcs3(i) + cs3(i) ! sig(e+p->e+p+3jets) tcse3(i) = tcse3(i) + cs3(i)*cs3(i) ! sig(gamma*+p->p+3jets) enddo C C Cross sections for different classes of radiative events C if(radcor_i) then if(.not.radcor_a) then tcsi3 = tcsi3 + cs3(1) tcsie3 = tcsie3 + cs3(1)*cs3(1) else tcsa3 = tcsa3 + cs3(1) tcsae3 = tcsae3 + cs3(1)*cs3(1) endif elseif(radcor_f) then tcsf3 = tcsf3 + cs3(1) tcsfe3 = tcsfe3 + cs3(1)*cs3(1) else tcsb3 = tcsb3 + cs3(1) tcsbe3 = tcsbe3 + cs3(1)*cs3(1) endif C if(ievent.ne.0) WEIGHT = CS3(1)/(1.-twojets) C C Check phase space of generated variables C if(ievent.ne.0) CALL GENERA3 C endif C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE ELEPOM(Iskip) C ------------------------------------------------------------------- C To compute quantities in lepton-pomeron reference system C Called by : EVENT3 C C Author : Ada Solano C ------------------------------------------------------------------- C include 'rboost.inc' include 'ridseq.inc' include 'ridkin.inc' C real escat,costhee,sinthee,phe,thee,thee2 real qsq1,qsq2,qsq3,qsq4,exbj1,exbj2,exbj3,eeee ISKIP = 0 C C Transform to CMS of electron and pomeron. Pomeron along +z' axis. C DO 1 I=1,5 PA(1,I) = POM(I) PA(2,I) = EB(I) PA(3,I) = 0. PA(4,I) = 0. 1 CONTINUE DO 10 J=1,3 10 DBETA(1,J)=DBLE(PA(1,J)+PA(2,J))/DBLE(PA(1,4)+PA(2,4)) dbe = dsqrt(dbeta(1,1)**2+dbeta(1,2)**2+dbeta(1,3)**2) if(dbe.gt.1.) then nskip = nskip+1 iskip = 1 return endif CALL DBOOST3(-DBETA(1,1),-DBETA(1,2),-DBETA(1,3)) DPHI(1)=ULANGL(PA(1,1),PA(1,2)) CALL ROTATE3(0.D+0,-DPHI(1)) DTHETA(1)=ULANGL(PA(1,3),PA(1,1)) CALL ROTATE3(-DTHETA(1),0.D+0) C C Total center of mass energy C EPOMS = (PA(1,4)+PA(2,4))**2 - (PA(1,3)+PA(2,3))**2 C C Generate the transverse momentum of the hard gluon C inside the pomeron C AKT2MX = MIN(AKT2C(2),EPOMS/4-amgsq) IF(AKT2MN.GE.AKT2MX) THEN if(ievent.ge.1) IAKT2 = IAKT2+1 nskip = nskip+1 ISKIP = 1 RETURN ENDIF C if(wei) then ! like 1/(K_t**2+amg**2)**2 C FACKT2 = 1. - (AKT2MN+AMGSQ)/(AKT2MX+AMGSQ) AKT2 = (AKT2MN+AMGSQ)/(1.-FACKT2*userran(iseed)) - AMGSQ AMKT2 = AKT2+AMGSQ PFKT2 = (AKT2MN+AMGSQ)/(FACKT2*AMKT2*AMKT2) C else ! like 1/(K_t**2+amg**2/BBB)**(BB+1) C FACKT2 = ((AKT2MN+AMGSQ/BBB)/(AKT2MX+AMGSQ/BBB))**BB - 1. AKT2 = (AKT2MN+AMGSQ/BBB)/ + ((1.+FACKT2*userran(iseed))**(1./BB)) - AMGSQ/BBB PFKT2 = -BB*(AKT2MN+AMGSQ/BBB)**BB / + (FACKT2*(AKT2+AMGSQ/BBB)**(BB+1.)) AMKT2 = AKT2+AMGSQ C endif C alps = alphas(amkt2) if(alps.gt.alpscut) then if(ievent.ge.1) ialps = ialps+1 alps = alpscut endif C WTKT2 = (9.*ALPS*ALPS*ALPS/(64.*PI)/(AMKT2*AMKT2)) & *akt2*akt2/(akt2+akt2cut)/(akt2+akt2cut) C C Generate QSQ of the photon C QSQ2 = AMKT2/(1.-SQRT(4.*AMKT2/(EPOMS-T-eb(5)))) QSQMIN = MAX(QSQC(1),AKT2,QSQ2) QSQ3 = (EPOMS-T-eb(5))*(1.-SQRT(4.*AMKT2/(EPOMS-T-eb(5)))) qsq4 = EPOMS-4*AMKT2 QSQMAX = MIN(qsq4,QSQ3,qsqc(2)) IF(QSQMIN.GE.QSQMAX) THEN if(ievent.ge.1) IQSQ = IQSQ+1 nskip = nskip+1 ISKIP = 1 RETURN ENDIF C if(wei) then ! like 1/Q**4 C FACQSQ = 1. - (QSQMIN/QSQMAX) QSQ = QSQMIN/(1.-FACQSQ*userran(iseed)) PFQSQ = QSQMIN/(FACQSQ*QSQ*QSQ) C else ! like 1/(Q**2+CCC)**(CC+1) C FACQSQ = ((QSQMIN+CCC)/(QSQMAX+CCC))**CC - 1. QSQ = (QSQMIN+CCC)/((1.+FACQSQ*userran(iseed))**(1./CC)) - CCC PFQSQ = -CC*(QSQMIN+CCC)**CC/FACQSQ/((QSQ+CCC)**(CC+1.)) C endif C C Generate the effective Bjorken_x variable in the CMS of C electron and pomeron C exbj1 = qsq/(amx2c(2)+qsq-t) exbj1 = exbj1*0.95 EXBJMN = MAX(QSQ/(EPOMS-T-eb(5)),4*AMKT2/(EPOMS-QSQ),exbj1) EXBJ2 = (4.*AMKT2+2.*QSQ-T - + SQRT((4.*AMKT2+2.*QSQ-T)**2-4.*QSQ*(QSQ-T)))/(2.*(QSQ-T)) exbj3 = qsq/(amx2c(1)+qsq-t) exbj3 = exbj3*1.05 EXBJMX = MIN((1.-AMKT2/QSQ)**2,QSQ/(QSQ+4*AMKT2),EXBJ2,exbj3) IF (EXBJMN.GE.EXBJMX) THEN if(ievent.ge.1) IEXBJ = IEXBJ+1 nskip = nskip+1 ISKIP = 1 RETURN ENDIF C if(wei) then ! like 1/x C FACEXBJ = ALOG(EXBJMX/EXBJMN) EXBJ = EXBJMN*EXP(FACEXBJ*userran(iseed)) PFEXBJ = 1./(FACEXBJ*EXBJ) C else ! like 1/x(x+EEE) C EEEE = (EXBJMN+EEE)/EXBJMN if(EEEE*EXBJMX/(EXBJMX+EEE).eq.1.) then if(ievent.ge.1) IEXBJ = IEXBJ+1 nskip = nskip+1 ISKIP = 1 RETURN endif FACEXBJ = ALOG(EEEE*EXBJMX/(EXBJMX+EEE)) EXBJ = EEE/(EEEE*EXP(-FACEXBJ*userran(iseed))-1.) PFEXBJ = EEE/(FACEXBJ*EXBJ*(EXBJ+EEE)) C endif C C Effective Bjorken_y variable C EY = QSQ/EXBJ/(EPOMS-T-eb(5)) C C Generate the scattered electron, differently if initial state C radiation occurs or not C if(radcor_i) then C C Energy and theta of the scattered electron (from Q2 and y) C taking into account the virtuality of the incoming electron C escat = (1.-ey)*(pa(1,4)*pa(2,4)+pa(1,3)*pa(1,3)) + + (qsq+pa(2,5)+amesq)/2. escat = escat/(pa(1,4)+pa(2,4)) costhee = pa(2,4)*escat - (qsq+pa(2,5)+amesq)/2. efmod = escat*escat-amesq if(efmod.lt.0.) then write(6,*)'WARNING-ELEPOM: efmod < 0 ',ievent,efmod nskip = nskip+1 iskip = 1 return endif costhee = costhee/pa(1,3)/sqrt(efmod) if(costhee.gt.1.) then if(costhee.lt.1.005) then costhee = 1. else nskip = nskip+1 iskip = 1 return endif endif sinthee = sqrt(1.-costhee*costhee) PHE = 2.*PI*userran(iseed) PFPHIE = 1.0 C C Fill vector of scattered electron C PA(3,1) = ESCAT*SINTHEE*COS(PHE) PA(3,2) = ESCAT*SINTHEE*SIN(PHE) PA(3,3) = -ESCAT*COSTHEE PA(3,4) = ESCAT PA(3,5) = amesq C C Back to the Zeus reference system C CALL ROTATE3(DTHETA(1),0.D+0) CALL ROTATE3(0.D+0,DPHI(1)) CALL DBOOST3(DBETA(1,1),DBETA(1,2),DBETA(1,3)) C THETAE = ATAN2(SQRT(PA(3,1)*PA(3,1)+PA(3,2)*PA(3,2)),PA(3,3)) thetae = pi-thetae C C Check if the initial state radiation is absorbed by the C scattered electron (interference term) C itype = 2 call radcorr(itype) if(radcor_a) then C C Recalculate angles and momentum of the scattered electron C (in electron-pomeron c.m.s.) C do 20 i=1,5 20 eb(i) = eb_i(i) s = stot(5) C DO 30 I=1,5 PA(1,I) = POM(I) PA(2,I) = EB(I) PA(3,I) = 0. PA(4,I) = 0. 30 CONTINUE DO 40 J=1,3 40 DBETA(1,J)=DBLE(PA(1,J)+PA(2,J))/DBLE(PA(1,4)+PA(2,4)) CALL DBOOST3(-DBETA(1,1),-DBETA(1,2),-DBETA(1,3)) DPHI(1)=ULANGL(PA(1,1),PA(1,2)) CALL ROTATE3(0.D+0,-DPHI(1)) DTHETA(1)=ULANGL(PA(1,3),PA(1,1)) CALL ROTATE3(-DTHETA(1),0.D+0) C EPOMS = (PA(1,4)+PA(2,4))**2 - (PA(1,3)+PA(2,3))**2 EY = QSQ/EXBJ/(EPOMS-T-eb(5)) C ESCAT = (QSQ+2*AMESQ)/(4.*PA(2,4)) + (1.-EY)*PA(2,4) THEE2 = (QSQ+2*AMESQ)/(4.*PA(2,4)*ESCAT) IF(THEE2.GT.1.) THEN THEE = PI/2. ELSE THEE = 2.*ASIN(SQRT(THEE2)) ENDIF PHE = 2.*PI*userran(iseed) PFPHIE = 1. C C Fill vector of scattered electron C PA(3,1) = ESCAT*SIN(THEE)*COS(PHE) PA(3,2) = ESCAT*SIN(THEE)*SIN(PHE) PA(3,3) = -ESCAT*COS(THEE) PA(3,4) = ESCAT PA(3,5) = amesq C C Back to the Zeus reference system C CALL ROTATE3(DTHETA(1),0.D+0) CALL ROTATE3(0.D+0,DPHI(1)) CALL DBOOST3(DBETA(1,1),DBETA(1,2),DBETA(1,3)) C endif C else C C Angles and momentum of scattered electron (electron-pomeron CMS!!) C THEE --> pi - THEE (electron along -z' axis) C ESCAT = (QSQ+2*AMESQ)/(4.*PA(2,4)) + (1.-EY)*PA(2,4) THEE2 = (QSQ+2*AMESQ)/(4.*PA(2,4)*ESCAT) IF(THEE2.GT.1.) THEN THEE = PI/2. ELSE THEE = 2.*ASIN(SQRT(THEE2)) ENDIF PHE = 2.*PI*userran(iseed) PFPHIE = 1. C C Fill vector of scattered electron C PA(3,1) = ESCAT*SIN(THEE)*COS(PHE) PA(3,2) = ESCAT*SIN(THEE)*SIN(PHE) PA(3,3) = -ESCAT*COS(THEE) PA(3,4) = ESCAT PA(3,5) = amesq C C Back to the Zeus reference system C CALL ROTATE3(DTHETA(1),0.D+0) CALL ROTATE3(0.D+0,DPHI(1)) CALL DBOOST3(DBETA(1,1),DBETA(1,2),DBETA(1,3)) C endif C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE GAMPOM3(Iskip) C ------------------------------------------------------------------- C To compute quantities in gamma-pomeron reference system C Called by : EVENT3 C C Author : Ada Solano C ------------------------------------------------------------------- C include 'rboost.inc' include 'ridseq.inc' include 'ridkin.inc' C real aq,ag,amcsq,amssq,pippo,pluto,check real xgmin,xgmax,akplus,akminus ISKIP = 0 C C Transform to CMS of photon and pomeron. Pomeron along +z' axis. C DO 1 I=1,5 PA(1,I) = POM(I) PA(2,I) = GAM(I) PA(3,I) = 0. PA(4,I) = 0. 1 CONTINUE DO 10 J=1,3 10 DBETA(1,J)=DBLE(PA(1,J)+PA(2,J))/DBLE(PA(1,4)+PA(2,4)) CALL DBOOST3(-DBETA(1,1),-DBETA(1,2),-DBETA(1,3)) DPHI(1)=ULANGL(PA(1,1),PA(1,2)) CALL ROTATE3(0.D+0,-DPHI(1)) DTHETA(1)=ULANGL(PA(1,3),PA(1,1)) CALL ROTATE3(-DTHETA(1),0.D+0) c c Check if charm or strange production is possible c if(onlych) then amcsq = 1.94*1.94 amssq = 0.56*0.56 ag = (amgsq+akt2)/amx2 aq = (4*amcsq+akt2)/amx2 check=1.-exbj-(sqrt(aq)+sqrt(ag))**2 if(check.gt.0) then charm = .true. iflav = 4 goto 100 else iskip = 1 nskip = nskip + 1 return endif endif c if(.not.lightf) then amcsq = 1.94*1.94 amssq = 0.56*0.56 ag = (amgsq+akt2)/amx2 aq = (4*amcsq+akt2)/amx2 check=1.-exbj-(sqrt(aq)+sqrt(ag))**2 if(check.gt.0) then charm = .true. else aq = (4*amssq+akt2)/amx2 check=1.-exbj-(sqrt(aq)+sqrt(ag))**2 if(check.gt.0) strange = .true. endif endif C C Flavour of the initiators ( Lund ): 1 = down, 2 = up C 3 = strange, 4 = charm C IFLAV = 1 PIPPO = userran(iseed) IF(PIPPO.GT.0.2) IFLAV = IFLAV+1 if(.not.lightf) then PLUTO = userran(iseed) if(charm) then IF(PLUTO.GT.0.5) IFLAV = IFLAV+2 elseif(strange) then if(pluto.lt.1./6.) iflav = 3 endif endif C if(iflav.eq.3) aq = (4*amssq+akt2)/amx2 if(iflav.eq.4) aq = (4*amcsq+akt2)/amx2 C 100 continue if(iflav.le.2) then if((1.-EXBJ)*((1.-EXBJ)/4.-AMKT2/AMX2).lt.0)then write(6,*)'WARNING-GAMPOM3: No phase space for XG! ',ievent endif xxx = max( 0.,(1.-EXBJ)*((1.-EXBJ)/4.-AMKT2/AMX2) ) XGMAX = (1.+EXBJ)/2. + sqrt(xxx) XGMIN = (1.+EXBJ)/2. - sqrt(xxx) else XGMax = (1.+exbj+aq-ag)/2. + SQRT((1.-exbj-aq-ag)**2/4.-aq*ag) XGMin = (1.+exbj+aq-ag)/2. - SQRT((1.-exbj-aq-ag)**2/4.-aq*ag) endif C C Generate z = exbj/x_g as (1-z)**FF , where x_g is the C momentum fraction of the hard gluon inside the pomeron C ZMIN = EXBJ/XGMAX ZMAX = EXBJ/XGMIN IF(ZMIN.GE.ZMAX) THEN if(ievent.ge.1) IZ = IZ+1 nskip = nskip+1 ISKIP = 1 RETURN ENDIF FACZ = (1-ZMAX)**(FF+1) - (1-ZMIN)**(FF+1) Z = 1 - ( userran(iseed)*FACZ + (1-ZMIN)**(FF+1) )**(1./(FF+1)) PFZ = ( -(FF+1)*(1-Z)**FF )/FACZ C if(lightf) then WTZ = (Z*Z+(1.-Z)*(1.-Z))*10.*4*PI*PI*ALPHA/9./QSQ elseif(onlych) then WTZ = (Z*Z+(1.-Z)*(1.-Z))*8.*4*PI*PI*ALPHA/9./QSQ else if(charm) then WTZ = (Z*Z+(1.-Z)*(1.-Z))*20.*4*PI*PI*ALPHA/9./QSQ elseif(strange) then WTZ = (Z*Z+(1.-Z)*(1.-Z))*12.*4*PI*PI*ALPHA/9./QSQ else WTZ = (Z*Z+(1.-Z)*(1.-Z))*10.*4*PI*PI*ALPHA/9./QSQ endif endif C C Momentum fraction of the hard gluon inside the pomeron C XG = EXBJ/Z C C Momentum fraction of the hard gluon remnant from the pomeron: C x_k = 1-x_g = E_k+k_z / E_pom+pom_z = k+/pom+ C XK = 1. - XG C AKPLUS = XK*(PA(1,4)+PA(1,3)) AKMINUS = (AMGSQ+AKT2)/AKPLUS C C Hard gluon remnant C PHIK = 2*PI*userran(iseed) AKT = SQRT(AKT2) PA(3,1) = AKT*COS(PHIK) PA(3,2) = AKT*SIN(PHIK) PA(3,4) = 0.5*(AKPLUS+AKMINUS) PA(3,3) = AKPLUS - PA(3,4) PA(3,5) = AMGSQ C C Back to the Zeus reference system C CALL ROTATE3(DTHETA(1),0.D+0) CALL ROTATE3(0.D+0,DPHI(1)) CALL DBOOST3(DBETA(1,1),DBETA(1,2),DBETA(1,3)) C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE GAMGLU(Iskip) C ------------------------------------------------------------------- C To compute quantities in gamma-gluon reference system C Called by : EVENT3 C C Author : Ada Solano C ------------------------------------------------------------------- C include 'rboost.inc' include 'ridseq.inc' include 'ridkin.inc' C REAL EMF,THJCM2,THJ0SQ ISKIP = 0 C C Transform to CMS of photon and gluon. Gluon along +z' axis. C DO 1 I=1,5 PA(1,I) = GLUON(I) PA(2,I) = GAM(I) PA(3,I) = 0. PA(4,I) = 0. 1 CONTINUE DO 10 J=1,3 10 DBETA(1,J)=DBLE(PA(1,J)+PA(2,J))/DBLE(PA(1,4)+PA(2,4)) CALL DBOOST3(-DBETA(1,1),-DBETA(1,2),-DBETA(1,3)) if(ibskip.eq.1) then nskip = nskip+1 iskip = 1 return endif DPHI(1)=ULANGL(PA(1,1),PA(1,2)) CALL ROTATE3(0.D+0,-DPHI(1)) DTHETA(1)=ULANGL(PA(1,3),PA(1,1)) CALL ROTATE3(-DTHETA(1),0.D+0) C EJET = SQRT(QQBAR(5))/2. C C Skip the event if EJET**2 is less than the mass**2 of the quark C emf = effmas(iflav) if(iflav.eq.4) emf = 1.94 ! to avoid problems with JETSET IF(emf.GT.EJET) THEN if(ievent.ge.1) IAMF2 = IAMF2+1 nskip = nskip+1 ISKIP = 1 RETURN ENDIF C C Generate theta distribution of the 2 jets like C 1/(theta**2 + kt**2/ejet**2) ( 0-pi range ) C THJMN2 = AKT2/EJET/EJET IF(THJMN2.GE.THJMX2) THEN if(ievent.ge.1) ITHJCM = ITHJCM+1 nskip = nskip+1 ISKIP = 1 RETURN ENDIF THJ0SQ = AKT2/EJET/EJET FACTHJ = ALOG( (THJMX2+THJ0SQ)/(THJMN2+THJ0SQ) ) THJCM2 = (THJMN2+THJ0SQ)*EXP(FACTHJ*userran(iseed)) - THJ0SQ THJCM = SQRT(THJCM2) PFTHJC = 1./FACTHJ/(THJCM2+THJ0SQ) ! theta**2 PFTHJC = PFTHJC*(1.-Z)/EJET/EJET ! q_int**2 C C Fill vectors of the 2 jets ( photon-gluon CMS!! ) C The mass of the jets is NOT taken to be zero C PHJCM = 2.*PI*userran(iseed) PA(3,5) = ULMASS(IFLAV)**2 PA(3,4) = EJET PJET = SQRT( EJET*EJET-PA(3,5) ) PA(3,1) = PJET*SIN(THJCM)*COS(PHJCM) PA(3,2) = PJET*SIN(THJCM)*SIN(PHJCM) PA(3,3) = PJET*COS(THJCM) PA(4,1) = -PA(3,1) PA(4,2) = -PA(3,2) PA(4,3) = -PA(3,3) PA(4,4) = PA(3,4) PA(4,5) = PA(3,5) C q_kt2 = pa(3,1)**2+pa(3,2)**2 C C Invariant mass of the off-shell quark C DO 5 I=1,5 GLUECM(I) = PA(1,I) Q1CM(I) = PA(3,I) 5 CONTINUE CALL DIF5(GLUECM,Q1CM,QINT) C C Evaluate that part of cross section which is function of q_int**2 C WTQINT = 1./(ABS(QINT(5))+effmas(iflav)**2) C C Back to the Zeus reference system C CALL ROTATE3(DTHETA(1),0.D+0) CALL ROTATE3(0.D+0,DPHI(1)) CALL DBOOST3(DBETA(1,1),DBETA(1,2),DBETA(1,3)) C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE FRAGM3 C -------------------------------------------------------------------- C Parton shower + hadronization of the q-qbar-gluon system (JETSET) C Called by : RIDIEVT C C Author : Ada Solano C -------------------------------------------------------------------- C include 'rboost.inc' include 'ridseq.inc' include 'ridkin.inc' C COMMON/LUJETS/N,K(4000,5),P(4000,5),V(4000,5) COMMON/LUDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) common/ludat3/mdcy(500,3),mdme(2000,2),brat(2000),kfdp(2000,5) EXTERNAL LUDATA C integer id,ide,idp,ijoin3(3) REAL QMAX,SCFPX,SCFPY,SCFPZ,ENEW,VAREN,angphi,angthe C MSTU(1) = 0 MSTU(2) = 0 MSTU(3) = 0 C C Initial and scattered beam particles + radiated photons C id = 1 ide = id call lu1ent(id,-11*LFLAG,eb_i(4),paru(1),0.) k(id,1) = 21 k(id,2) = -11*LFLAG k(id,3) = 0 do 1 i=1,4 1 p(id,i) = eb_i(i) p(id,5) = ame c id = id+1 idp = id call lu1ent(id,2212,pb(4),0.,0.) k(id,1) = 21 k(id,2) = 2212 k(id,3) = 0 do 2 i=1,4 2 p(id,i) = pb(i) p(id,5) = amp c if(radcor_i.and.(.not.radcor_a)) then id = id+1 angphi = ULANGL(phot(1),phot(2)) angthe = ATAN2(SQRT(phot(1)*phot(1)+phot(2)*phot(2)),phot(3)) call lu1ent(id,22,phot(4),angthe,angphi) k(id,1) = 21 k(id,2) = 22 k(id,3) = ide do 3 i=1,4 3 p(id,i) = phot(i) p(id,5) = 0. c id = id+1 angphi = ULANGL(eb(1),eb(2)) angthe = ATAN2(SQRT(eb(1)*eb(1)+eb(2)*eb(2)),eb(3)) call lu1ent(id,-11*LFLAG,eb(4),angthe,angphi) k(id,1) = 21 k(id,2) = -11*LFLAG k(id,3) = ide do 4 i=1,4 4 p(id,i) = eb(i) p(id,5) = -sqrt(abs(eb(5))) ide = id endif c id = id+1 angphi = ULANGL(es(1),es(2)) angthe = ATAN2(SQRT(es(1)*es(1)+es(2)*es(2)),es(3)) call lu1ent(id,-11*LFLAG,es(4),angthe,angphi) k(id,1) = 21 k(id,2) = -11*LFLAG k(id,3) = ide do 5 i=1,4 5 p(id,i) = es(i) p(id,5) = sqrt(es(5)) ide = id c if(radcor_f) then id = id+1 angphi = ULANGL(phot(1),phot(2)) angthe = ATAN2(SQRT(phot(1)*phot(1)+phot(2)*phot(2)),phot(3)) call lu1ent(id,22,phot(4),angthe,angphi) k(id,1) = 21 k(id,2) = 22 k(id,3) = ide do 6 i=1,4 6 p(id,i) = phot(i) p(id,5) = 0. c id = id+1 angphi = ULANGL(es_f(1),es_f(2)) angthe = ATAN2(SQRT(es_f(1)*es_f(1)+es_f(2)*es_f(2)),es_f(3)) call lu1ent(id,-11*LFLAG,es_f(4),angthe,angphi) k(id,1) = 21 k(id,2) = -11*LFLAG k(id,3) = ide do 7 i=1,4 7 p(id,i) = es_f(i) p(id,5) = ame endif c id = id+1 call lu1ent(id,2212,ps(4),thetap,phip) k(id,1) = 21 k(id,2) = 2212 k(id,3) = idp do 8 i=1,4 8 p(id,i) = ps(i) p(id,5) = amp c do 9 i=1,id k(i,4) = 0 k(i,5) = 0 9 continue C id = id+1 idq = id ijoin3(1) = id ijoin3(2) = id+2 ijoin3(3) = id+1 C IF(IPS.EQ.0) THEN C C Parton shower and fragmentation are done in the Zeus ref. system C CALL LU1ENT(-id,IFLAV,Q(4),THEQ,PHIQ) CALL LU1ENT(-id-1,-IFLAV,QBAR(4),THEQB,PHIQB) C C Enter the gluon mass ( the default is mg=0 ) C MSTU(10) = 1 P(id+2,5) = AMG CALL LU1ENT(id+2,21,AK(4),THETAK,PHIAK) MSTU(10) = 2 C C Connect the colour C CALL LUJOIN(3,IJOIN3) C if(mstj(41).gt.0) then C C Parton shower: only for the q-qbar system C qmax = plu(id,12) CALL LUSHOW(id,id+1,QMAX) C C Rearrange parton shower endproducts C DO 10 J=1,N IF(K(J,1).EQ.3)CALL LUPREP(J) 10 CONTINUE C endif C CALL LUEXEC C ELSEIF(IPS.EQ.1) THEN C C Parton shower and fragmentation are done in the q-qbar CMS system, C quark along +z' axis. C DO 100 I=1,5 PA(1,I) = Q(I) PA(2,I) = QBAR(I) PA(3,I) = AK(I) PA(4,I) = 0. 100 CONTINUE DO 110 J=1,3 110 DBETA(1,J)=DBLE(PA(1,J)+PA(2,J))/DBLE(PA(1,4)+PA(2,4)) CALL DBOOST3(-DBETA(1,1),-DBETA(1,2),-DBETA(1,3)) DPHI(1)=ULANGL(PA(1,1),PA(1,2)) CALL ROTATE3(0.D+0,-DPHI(1)) DTHETA(1)=ULANGL(PA(1,3),PA(1,1)) CALL ROTATE3(-DTHETA(1),0.D+0) C CALL LU1ENT(-id,IFLAV,PA(1,4),0.,0.) CALL LU1ENT(-id-1,-IFLAV,PA(2,4),PARU(1),0.) C C Enter the gluon mass ( the default is mg=0 ) C MSTU(10) = 1 P(id+2,5) = AMG PHIKCM = ULANGL(PA(3,1),PA(3,2)) THEKCM = ATAN2(SQRT(PA(3,1)*PA(3,1)+PA(3,2)*PA(3,2)),PA(3,3)) CALL LU1ENT(id+2,21,PA(3,4),THEKCM,PHIKCM) MSTU(10) = 2 C C Connect the colour C CALL LUJOIN(3,IJOIN3) C if(mstj(41).gt.0) then C C Parton shower: first consider the gluon alone C QMAX = MIN( SQRT(AMKT2), + SQRT((2.*P(id,4)+P(id+2,4)-2.*P(id,5))**2 - P(id+2,4)**2)) 115 CALL LUSHOW(id+2,0,QMAX) IF(N.EQ.id+2) GOTO 150 C C Since LUSHOW does not conserve the momentum, scale all the partons C from the gluon to obtain momentum conservation again C SCFPX = P(id+2,1)/P(id+3,1) SCFPY = P(id+2,2)/P(id+3,2) SCFPZ = P(id+2,3)/P(id+3,3) ENEW = 0. DO 120 I=id+4,N IF(K(I,1).EQ.3) THEN P(I,1) = P(I,1)*SCFPX P(I,2) = P(I,2)*SCFPY P(I,3) = P(I,3)*SCFPZ C C Calculate the new energies C P(I,4) = SQRT( P(I,1)*P(I,1)+P(I,2)*P(I,2)+P(I,3)*P(I,3) ) ENEW = ENEW + P(I,4) ENDIF 120 CONTINUE VAREN = ENEW - P(id+2,4) C C Rearrange energy and momentum of the q-qbar system to get C total energy and momentum conservation C P(id,4) = P(id,4) - VAREN/2. P(id+1,4) = P(id+1,4) - VAREN/2. IF( (P(id,4).LT.P(id,5)).OR.(P(id+1,4).LT.P(id+1,5)) ) THEN QMAX = 0.8*QMAX DO 125 I=id+3,N DO 126 J=1,5 P(I,J) = 0. K(I,J) = 0. 126 CONTINUE 125 CONTINUE N = id+2 P(id,4) = PA(1,4) P(id+1,4) = PA(2,4) K(id+2,1) = 3 K(id+2,4) = 60000 K(id+2,5) = 50000 GOTO 115 ENDIF P(id,3) = SQRT( P(id,4)*P(id,4)-P(id,5)*P(id,5) ) P(id+1,3) = -SQRT( P(id+1,4)*P(id+1,4)-P(id+1,5)*P(id+1,5) ) C DO 130 I=1,5 PA(1,I) = P(id,I) PA(2,I) = P(id+1,I) PA(3,I) = P(id+3,I) PA(4,I) = 0. 130 CONTINUE C 150 CONTINUE C C Parton shower: now consider the q-qbar system C QMAX = plu(id,12) CALL LUSHOW(id,id+1,QMAX) C DO 160 J=1,N IF(K(J,1).EQ.3)CALL LUPREP(J) 160 CONTINUE C endif C CALL LUEXEC C C Back to the Zeus reference system ( for all the entries of C the COMMON/LUJETS/ ) C theta(1) = sngl(dtheta(1)) phi(1) = sngl(dphi(1)) MSTU(1) = id DO 170 I=1,3 170 BETA(1,I) = SNGL(DBETA(1,I)) CALL LUROBO(theta(1),phi(1),BETA(1,1),BETA(1,2),BETA(1,3)) MSTU(1) = 0 C C Back to the Zeus reference system for the new three vectors C of the remnant gluon and the q-qbar system C CALL ROTATE3(DTHETA(1),0.D+0) CALL ROTATE3(0.D+0,DPHI(1)) CALL DBOOST3(DBETA(1,1),DBETA(1,2),DBETA(1,3)) C ENDIF c c Scattered electron and proton + radiatied photon as final states c k(id-1,1) = 1 k(id-2,1) = 1 if(radcor_i.and.(.not.radcor_a)) k(id-4,1) = 1 if(radcor_f) k(id-3,1) = 1 C C not in ZDIS if(ievent.le.nevpri) call lulist(1) C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE GENERA3 C ------------------------------------------------------------------- C To check generated variable integrals in case of weighted events C Called by : CROSS3 C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' include 'ridkin.inc' C C Test on generated variables: the integrals should be 1! C real vol3(7),wint3(7) C C Volumes: C vol3(1) = xfmax - xfmin vol3(2) = pt2max - pt2min vol3(3) = akt2mx - akt2mn vol3(4) = qsqmax - qsqmin vol3(5) = exbjmx - exbjmn vol3(6) = zmax - zmin vol3(7) = thjmx2 - thjmn2 C C Generation integrals: C wint3(1) = 1./pfxf/vol3(1) wint3(2) = 1./pfpt2/vol3(2) wint3(3) = 1./pfkt2/vol3(3) wint3(4) = 1./pfqsq/vol3(4) wint3(5) = 1./pfexbj/vol3(5) wint3(6) = 1./pfz/vol3(6) wint3(7) = 1./(pfthjc/((1.-Z)/EJET/EJET))/vol3(7) C do i = 1,7 twint3(i) = twint3(i) + wint3(i) twinte3(i) = twinte3(i) + wint3(i)*wint3(i) enddo C RETURN END C C ------------------------------------------------------------------- C SUBROUTINE RIDIOUT C ------------------------------------------------------------------- C Run output C Called by : MAIN C C Author : Ada Solano C ------------------------------------------------------------------- C include 'ridseq.inc' C integer nev2,nev3 C C Write last random number back to RIDI.Rnd for future use C (not in ZDIS) write(10,*) iseed C WRITE(6,*) WRITE(6,*) WRITE(6,*)'*****************************************************', + '******************' WRITE(6,*)'*****************************************************', + '******************' WRITE(6,*)'* ' WRITE(6,*)'* ' WRITE(6,*)'* RIDI 2.0 ' WRITE(6,*)'* ' WRITE(6,*)'* ' WRITE(6,*)'* MONTE CARLO CONFIGURATION FOR THIS RUN ' WRITE(6,*)'* ' WRITE(6,*)'* ' if(radcor) then WRITE(6,*)'* Radiative corrections on.' else WRITE(6,*)'* Radiative corrections off.' endif if(wei) then WRITE(6,*)'* Weighted events have been generated ' else WRITE(6,*)'* Unweighted events have been generated ' endif if(duejet.and.trejet) then WRITE(6,*)'* in q-qbar and q-qbar-gluon final states. ' elseif(duejet) then WRITE(6,*)'* in q-qbar final states only. ' else WRITE(6,*)'* in q-qbar-gluon final states only. ' endif WRITE(6,*)'* ' if(transv.and.longit) then WRITE(6,*)'* Transverse and longitudinal cross sections ' WRITE(6,*)'* have been simulated. ' elseif(transv) then WRITE(6,*)'* Only the transverse cross section ' WRITE(6,*)'* has been simulated. ' else WRITE(6,*)'* Only the longitudinal cross section ' WRITE(6,*)'* has been simulated. ' endif WRITE(6,*)'* ' if(lightf) + WRITE(6,*)'* Light flavours only! ' if(onlych) + WRITE(6,*)'* Charm flavour only! ' C WRITE(6,*)'* ' WRITE(6,*)'* ' WRITE(6,*)'* KINEMATICAL CUTS ' WRITE(6,*)'* ' WRITE(6,*)'* Q^2 cuts ',vcuts(3),vcuts(4) WRITE(6,*)'* x_Bj cuts ',vcuts(1),vcuts(2) WRITE(6,*)'* t cuts ',vcuts(5),vcuts(6) WRITE(6,*)'* Mx^2 cuts ',vcuts(7),vcuts(8) WRITE(6,*)'* x_pom cuts ',vcuts(9),vcuts(10) if(trejet) + WRITE(6,*)'* k_t^2 cuts ',vcuts(11),vcuts(12) C WRITE(6,*)'* ' WRITE(6,*)'* During generation the cut on M_x^2 or ' WRITE(6,*)'* x_pom which is more restrictive is used: ' if(duejet) then WRITE(6,*)'* ' WRITE(6,*)'* Minimum generated value of x_pom',XPOMLIM(1) WRITE(6,*)'* Maximum generated value of x_pom',XPOMLIM(2) endif if(trejet) then WRITE(6,*)'* ' WRITE(6,*)'* Minimum generated value of Mx^2',AMX2LIM(1) WRITE(6,*)'* Maximum generated value of Mx^2',AMX2LIM(2) endif C if(duejet) then WRITE(6,*)'* ' WRITE(6,*)'* ' akm=sqrt(akt2lim) WRITE(6,*)'* Minimum kt in q-qbar final states = ',akm,' GeV' endif C WRITE(6,*)'* ' WRITE(6,*)'* ' write(6,*)'* FRACTION OF Q-QBAR EVENTS =',twojets WRITE(6,*)'* ' c if(duejet) then c if(nevent2.eq.0) goto 33 C WRITE(6,*)'* ' WRITE(6,*)'*****************************************************', + '******************' WRITE(6,*)'* ' WRITE(6,*)'* STATISTICS OF Q-QBAR EVENTS ' WRITE(6,*)'* ' WRITE(6,*)'* Number of final events',NEVENT2 if(.not.wei) then WRITE(6,*)'* Number of generated events',ITOTEVE2 write(6,*)'* Maximun weight',weimax2 WRITE(6,*)'* Number of events with WEIGHT.GT.WEIMAX',IWEI2 WRITE(6,*)'* Number of events counted twice',ITWOEV2 endif WRITE(6,*)'* Number of skipped events',nmass c if(checkout) then WRITE(6,*)'* ' c WRITE(6,*)'* Number of skipped events with imass.ne.0', c + nmass WRITE(6,*)'* Number of generated events with aal>0.25', + laacheck WRITE(6,*)'* Number of generated events with pt2_cms<0.', + nptcms if(fragm) then WRITE(6,*)'* ' write(6,*)'* Minimum strange mass',stmin write(6,*)'* Minimum charm mass',chmin endif endif c nev2 = nevent2 if(.not.wei) NEVENT2 = ITOTEVE2 C WRITE(6,*)'* ' WRITE(6,*)'* ' WRITE(6,*)'* CROSS SECTIONS FOR THE Q-QBAR DIAGRAM ' WRITE(6,*)'* ' C DO I = 1,2 TCS2(I) = TCS2(I)/(NEVENT2+nmass) TCSE2(I) = SQRT(TCSE2(I))/(NEVENT2+nmass) ENDDO c DO I = 1,2 TCST2(I) = TCST2(I)/(NEVENT2+nmass) TCSTE2(I) = SQRT(TCSTE2(I))/(NEVENT2+nmass) ENDDO c DO I = 1,2 TCSL2(I) = TCSL2(I)/(NEVENT2+nmass) TCSLE2(I) = SQRT(TCSLE2(I))/(NEVENT2+nmass) ENDDO c TCSB2 = TCSB2/(NEVENT2+nmass) TCSBE2 = SQRT(TCSBE2)/(NEVENT2+nmass) TCSI2 = TCSI2/(NEVENT2+nmass) TCSIE2 = SQRT(TCSIE2)/(NEVENT2+nmass) TCSA2 = TCSA2/(NEVENT2+nmass) TCSAE2 = SQRT(TCSAE2)/(NEVENT2+nmass) TCSF2 = TCSF2/(NEVENT2+nmass) TCSFE2 = SQRT(TCSFE2)/(NEVENT2+nmass) c c ************ Output in nanobarn ********************************* c WRITE(6,*)'* SIG(e+p -> e+p+2jets) ' + ,TCS2(1),'+/-',TCSE2(1),' nb' WRITE(6,*)'* SIG_T(e+p -> e+p+2jets) ' + ,TCST2(1),'+/-',TCSTE2(1),' nb' WRITE(6,*)'* SIG_L(e+p -> e+p+2jets) ' + ,TCSL2(1),'+/-',TCSLE2(1),' nb' WRITE(6,*)'* ' WRITE(6,*)'* SIG_BORN(e+p -> e+p+2jets)' + ,TCSB2,'+/-',TCSBE2,' nb' WRITE(6,*)'* SIG_ISR(e+p -> e+p+2jets) ' + ,TCSI2,'+/-',TCSIE2,' nb' WRITE(6,*)'* SIG_ABS(e+p -> e+p+2jets) ' + ,TCSA2,'+/-',TCSAE2,' nb' WRITE(6,*)'* SIG_FSR(e+p -> e+p+2jets) ' + ,TCSF2,'+/-',TCSFE2,' nb' C if(wei) then WRITE(6,*)'* ' WRITE(6,*)'* SIG(gamma*+p -> p+2jets) ' + ,TCS2(2),'+/-',TCSE2(2),' nb' WRITE(6,*)'* SIG_T(gamma*+p -> p+2jets)' + ,TCST2(2),'+/-',TCSTE2(2),' nb' WRITE(6,*)'* SIG_L(gamma*+p -> p+2jets)' + ,TCSL2(2),'+/-',TCSLE2(2),' nb' endif C if(.not.wei.and.checkout) then WRITE(6,*)'* ' WRITE(6,*)'* These values might be meaningless, depending on ' WRITE(6,*)'* generation functions of Q^2 and x_Bj!! ' WRITE(6,*)'* SIG(gamma*+p -> p+2jets) ' + ,TCS2(2),'+/-',TCSE2(2),' nb' WRITE(6,*)'* SIG_T(gamma*+p -> p+2jets)' + ,TCST2(2),'+/-',TCSTE2(2),' nb' WRITE(6,*)'* SIG_L(gamma*+p -> p+2jets)' + ,TCSL2(2),'+/-',TCSLE2(2),' nb' endif C if(wei.and.checkout) then WRITE(6,*)'* ' WRITE(6,*)'* ' WRITE(6,*)'* GENERATION INTEGRALS FOR Q-QBAR EVENTS ' WRITE(6,*)'* ' DO I = 1,7 TWINT2(I) = TWINT2(I)/NEVENT2 TWINTE2(I) = SQRT(TWINTE2(I))/NEVENT2 WRITE(6,*)'* ',TWINT2(I),'+/-',TWINTE2(I) ENDDO endif C endif C 33 if(trejet) then c if(nevent3.eq.0) goto 99 c WRITE(6,*)'* ' WRITE(6,*)'*****************************************************', + '******************' WRITE(6,*)'* ' WRITE(6,*)'* STATISTICS OF Q-QBAR-GLUON EVENTS ' WRITE(6,*)'* ' WRITE(6,*)'* Number of final events',NEVENT3 if(.not.wei) then WRITE(6,*)'* Number of generated events',ITOTEVE3 write(6,*)'* Maximun weight',weimax3 WRITE(6,*)'* Number of events with WEIGHT.GT.WEIMAX',IWEI3 write(6,*)'* Number of events counted twice',ITWOEV3 endif WRITE(6,*)'* Number of skipped events',NSKIP C if(checkout) then WRITE(6,*)'* ' WRITE(6,*)'* IAKT2 =',IAKT2 WRITE(6,*)'* IQSQ =',IQSQ WRITE(6,*)'* IEXBJ =',IEXBJ WRITE(6,*)'* IAMX2 =',IAMX2 WRITE(6,*)'* IZ =',IZ WRITE(6,*)'* IAMF2 =',IAMF2 WRITE(6,*)'* ITHJCM =',ITHJCM WRITE(6,*)'* IALPS =',IALPS endif C nev3 = nevent3 if(.not.wei) NEVENT3 = ITOTEVE3 C WRITE(6,*)'* ' WRITE(6,*)'* ' WRITE(6,*)'* CROSS SECTIONS FOR THE Q-QBAR-GLUON DIAGRAM' WRITE(6,*)'* ' C DO I = 1,2 TCS3(I) = TCS3(I)/(NEVENT3+NSKIP) TCSE3(I) = SQRT(TCSE3(I))/(NEVENT3+NSKIP) ENDDO c DO I = 1,2 TCST3(I) = TCST3(I)/(NEVENT3+NSKIP) TCSTE3(I) = SQRT(TCSTE3(I))/(NEVENT3+NSKIP) ENDDO c DO I = 1,2 TCSL3(I) = TCSL3(I)/(NEVENT3+NSKIP) TCSLE3(I) = SQRT(TCSLE3(I))/(NEVENT3+NSKIP) ENDDO c TCSB3 = TCSB3/(NEVENT3+NSKIP) TCSBE3 = SQRT(TCSBE3)/(NEVENT3+NSKIP) TCSI3 = TCSI3/(NEVENT3+NSKIP) TCSIE3 = SQRT(TCSIE3)/(NEVENT3+NSKIP) TCSA3 = TCSA3/(NEVENT3+NSKIP) TCSAE3 = SQRT(TCSAE3)/(NEVENT3+NSKIP) TCSF3 = TCSF3/(NEVENT3+NSKIP) TCSFE3 = SQRT(TCSFE3)/(NEVENT3+NSKIP) c c ************ Output in nanobarn ********************************* c WRITE(6,*)'* SIG(e+p -> e+p+3jets) ' + ,TCS3(1),'+/-',TCSE3(1),' nb' WRITE(6,*)'* SIG_T(e+p -> e+p+3jets) ' + ,TCST3(1),'+/-',TCSTE3(1),' nb' WRITE(6,*)'* SIG_L(e+p -> e+p+3jets) ' + ,TCSL3(1),'+/-',TCSLE3(1),' nb' WRITE(6,*)'* ' WRITE(6,*)'* SIG_BORN(e+p -> e+p+3jets)' + ,TCSB3,'+/-',TCSBE3,' nb' WRITE(6,*)'* SIG_ISR(e+p -> e+p+3jets) ' + ,TCSI3,'+/-',TCSIE3,' nb' WRITE(6,*)'* SIG_ABS(e+p -> e+p+3jets) ' + ,TCSA3,'+/-',TCSAE3,' nb' WRITE(6,*)'* SIG_FSR(e+p -> e+p+3jets) ' + ,TCSF3,'+/-',TCSFE3,' nb' C if(wei) then WRITE(6,*)'* ' WRITE(6,*)'* SIG(gamma*+p -> p+3jets) ' + ,TCS3(2),'+/-',TCSE3(2),' nb' WRITE(6,*)'* SIG_T(gamma*+p -> p+3jets)' + ,TCST3(2),'+/-',TCSTE3(2),' nb' WRITE(6,*)'* SIG_L(gamma*+p -> p+3jets)' + ,TCSL3(2),'+/-',TCSLE3(2),' nb' endif C if(.not.wei.and.checkout) then WRITE(6,*)'* ' WRITE(6,*)'* These values might be meaningless, depending on ' WRITE(6,*)'* generation functions of Q^2 and x_Bj!! ' WRITE(6,*)'* SIG(gamma*+p -> p+3jets) ' + ,TCS3(2),'+/-',TCSE3(2),' nb' WRITE(6,*)'* SIG_T(gamma*+p -> p+3jets)' + ,TCST3(2),'+/-',TCSTE3(2),' nb' WRITE(6,*)'* SIG_L(gamma*+p -> p+3jets)' + ,TCSL3(2),'+/-',TCSLE3(2),' nb' endif C if(wei.and.checkout) then WRITE(6,*)'* ' WRITE(6,*)'* ' WRITE(6,*)'* GENERATION INTEGRALS FOR Q-QBAR-GLUON EVENTS' WRITE(6,*)'* ' DO I = 1,7 TWINT3(I) = TWINT3(I)/NEVENT3 TWINTE3(I) = SQRT(TWINTE3(I))/NEVENT3 WRITE(6,*)'* ',TWINT3(I),'+/-',TWINTE3(I) ENDDO endif c endif c WRITE(6,*)'* ' WRITE(6,*)'*****************************************************', + '******************' WRITE(6,*)'* ' WRITE(6,*)'* NUMBER OF FINAL EVENTS: ',NEV2+NEV3 WRITE(6,*)'* ' WRITE(6,*)'* TOTAL CROSS SECTION WITHIN THE CUTS: ' WRITE(6,*)'* ' err = sqrt(TCSE2(1)*TCSE2(1)+TCSE3(1)*TCSE3(1)) errt = sqrt(TCSTE2(1)*TCSTE2(1)+TCSTE3(1)*TCSTE3(1)) errl = sqrt(TCSLE2(1)*TCSLE2(1)+TCSLE3(1)*TCSLE3(1)) errb = sqrt(TCSBE2*TCSBE2+TCSBE3*TCSBE3) erri = sqrt(TCSIE2*TCSIE2+TCSIE3*TCSIE3) erra = sqrt(TCSAE2*TCSAE2+TCSAE3*TCSAE3) errf = sqrt(TCSFE2*TCSFE2+TCSFE3*TCSFE3) WRITE(6,*)'* SIG(e+p -> e+p+X) ' + ,TCS2(1)+TCS3(1),'+/-',err,' nb' WRITE(6,*)'* SIG_T(e+p -> e+p+X) ' + ,TCST2(1)+TCST3(1),'+/-',errt,' nb' WRITE(6,*)'* SIG_L(e+p -> e+p+X) ' + ,TCSL2(1)+TCSL3(1),'+/-',errl,' nb' WRITE(6,*)'* ' WRITE(6,*)'* SIG_BORN(e+p -> e+p+X)' + ,TCSB2+TCSB3,'+/-',errb,' nb' WRITE(6,*)'* SIG_ISR(e+p -> e+p+X) ' + ,TCSI2+TCSI3,'+/-',erri,' nb' WRITE(6,*)'* SIG_ABS(e+p -> e+p+X) ' + ,TCSA2+TCSA3,'+/-',erra,' nb' WRITE(6,*)'* SIG_FSR(e+p -> e+p+X) ' + ,TCSF2+TCSF3,'+/-',errf,' nb' WRITE(6,*)'* ' WRITE(6,*)'* ' WRITE(6,*)'*****************************************************', + '******************' WRITE(6,*)'*****************************************************', + '******************' C C NOT IN ZDIS C if(user) CALL USEROUT C 99 RETURN END C C ------------------------------------------------------------------- C REAL FUNCTION ALPHAS(SCALE2) C C *** ALPHA STRONG C common/strong/alpsm(9),alpst(9) logical first save first c data first /.true./ data alpsm /0.5,1.,2.,4.,8.,16.,32.,64.,91.2/ c c amc=charm mass, amb=beauty mass, a2=alpha_s pi=3.141592654 p2=2.*pi amc=1.35 amb=4.6 c if(first) then first = .false. np = 10000 h=6.9/np do 10 i=1,9 qq = alpsm(i) a2 = 0.120 do 11 ip=1,np nf=3 pt=exp(-ip*h)*91.187 if(pt.lt.0.6) goto 11 if(pt.lt.qq) goto 11 if(pt.gt.amc) nf=4 if(pt.gt.amb) nf=5 b00=(11.-(2./3.)*nf)/p2 b1=(51.-19.*nf/3.)/(p2*p2) b2=(2857.-5033*nf/9.-(325./27.)*nf**2)/(8.*p2*p2*p2) a2=a2+h*a2*a2*(b00+a2*(b1+a2*b2)) 11 continue alpst(i) = a2 10 continue endif c qq = sqrt(scale2) i = int(alog(qq)/alog(2.)) + 3 if(qq.lt.1.) i = i-1 if(i.lt.1) i = 1 if(i.gt.9) i = 9 a2 = alpst(i) np=100 h=alog(2.)/np do 1 ip=1,np nf=3 pt=exp(-ip*h)*alpsm(i) if(pt.lt.0.6) goto 2 if(pt.lt.qq) goto 2 if(pt.gt.amc) nf=4 if(pt.gt.amb) nf=5 b00=(11.-(2./3.)*nf)/p2 b1=(51.-19.*nf/3.)/(p2*p2) b2=(2857.-5033*nf/9.-(325./27.)*nf**2)/(8.*p2*p2*p2) a2=a2+h*a2*a2*(b00+a2*(b1+a2*b2)) 1 continue 2 continue alphas=a2 C RETURN END C C ------------------------------------------------------------------- C REAL FUNCTION GN(XX,XKT2) C C *** GLUON DENSITY DISTRIBUTION C include 'ridseq.inc' logical first real scale DOUBLE PRECISION Up,Dw,Usea,Dsea,St,Ch,Bo,Top,xgx,val(20) CHARACTER*20 PARM(20) save first c data first/.true./ if(first)then first = .false. PARM(1) = 'Nptype' PARM(2) = 'Ngroup' PARM(3) = 'Nset' VAL(1) = 1.D+0 ! proton VAL(2) = DBLE(gval(1)) ! 3=MRS; 5=GRV VAL(3) = DBLE(gval(2)) ! 43=MRS(A) 4=LO; 3,6=HO CALL PDFSET(PARM,VAL) endif c if(igluon.eq.0) then GN = (gpar(1)*((1.-xx)**gpar(2))*((xx/gpar(3))**gpar(4)))/xx else if(gsfk02.le.0.) gsfk02 = 1. if(xkt2.lt.gsfk02) then scale = sqrt(gsfk02) CALL STRUCTM + (DBLE(xx),DBLE(scale),Up,Dw,Usea,Dsea,St,Ch,Bo,Top,XGX) xgx = xgx*xkt2/gsfk02 else scale = sqrt(xkt2) CALL STRUCTM + (DBLE(xx),DBLE(scale),Up,Dw,Usea,Dsea,St,Ch,Bo,Top,XGX) endif GN = sngl(xgx)/xx endif c RETURN END C C ------------------------------------------------------------------- C SUBROUTINE ADD5(PLA,PLB,PLR) C C *** PLA,PLB ARE (MOMENTUM) FIVE-VECTORS OF TRACKS C *** PLR IS THE RESULTANT FIVE VECTOR OBTAINED BY ADDING PLA AND PLB C *** THEREFORE ............. C *** PLR(5) IS INVARIANT MASS SQUARED OF COMBINATION C DIMENSION PLA(5),PLB(5),PLR(5) C DO 10 K=1,4 10 PLR(K)=PLA(K)+PLB(K) PLR(5)=PLR(4)**2 DO 20 K=1,3 20 PLR(5)=PLR(5)-PLR(K)**2 RETURN END C C ------------------------------------------------------------------- C SUBROUTINE DIF5(PLA,PLB,PLR) C DIMENSION PLA(5),PLB(5),PLR(5) C DO 10 K=1,4 10 PLR(K)=PLA(K)-PLB(K) PLR(5)=PLR(4)**2 DO 20 K=1,3 20 PLR(5)=PLR(5)-PLR(K)**2 RETURN END C C ------------------------------------------------------------------- C FUNCTION ULANGL(X,Y) C C *** ROUTINE TO DETERMINE THE ANGLE FROM X AND Y COORDINATES C REAL PIGREC PIGREC = 4.0*ATAN(1.0) ULANGL=0. R=SQRT(X**2+Y**2) IF(R.LT.1E-20) RETURN IF(ABS(X)/R.LT.0.8) THEN ULANGL=SIGN(ACOS(X/R),Y) ELSE ULANGL=ASIN(Y/R) IF(X.LT.0..AND.ULANGL.GE.0.) THEN ULANGL=PIGREC-ULANGL ELSEIF(X.LT.0.) THEN ULANGL=-PIGREC-ULANGL ENDIF ENDIF RETURN END C C ------------------------------------------------------------------- C SUBROUTINE DBOOST2(BEX,BEY,BEZ) C C *** ROUTINE TO MAKE BOOSTS IN DOUBLE PRECISION. C include 'rboost.inc' C DOUBLE PRECISION BEX,BEY,BEZ,GA,BEP,GABEP C GA=1./DSQRT(1.D+0-BEX**2-BEY**2-BEZ**2) DO 140 I=1,4 BEP=BEX*DBLE(PA(I,1))+BEY*DBLE(PA(I,2))+BEZ*DBLE(PA(I,3)) GABEP=GA*(GA*BEP/(1.D+0+GA)+DBLE(PA(I,4))) PA(I,1)=SNGL(DBLE(PA(I,1))+GABEP*BEX) PA(I,2)=SNGL(DBLE(PA(I,2))+GABEP*BEY) PA(I,3)=SNGL(DBLE(PA(I,3))+GABEP*BEZ) PA(I,4)=SNGL(GA*(DBLE(PA(I,4))+BEP)) 140 CONTINUE RETURN END C C ------------------------------------------------------------------- C SUBROUTINE ROTATE2(DTH,DPH) C C *** ROUTINE TO ROTATE IN DOUBLE PRECISION. C include 'rboost.inc' C DIMENSION DROT(3,3),DPV(3) DROT(1,1)=DCOS(DTH)*DCOS(DPH) DROT(1,2)=-DSIN(DPH) DROT(1,3)=DSIN(DTH)*DCOS(DPH) DROT(2,1)=DCOS(DTH)*DSIN(DPH) DROT(2,2)=DCOS(DPH) DROT(2,3)=DSIN(DTH)*DSIN(DPH) DROT(3,1)=-DSIN(DTH) DROT(3,2)=0.D+0 DROT(3,3)=DCOS(DTH) DO 120 I=1,4 DO 100 J=1,3 100 DPV(J)=DBLE(PA(I,J)) DO 110 J=1,3 PA(I,J)=SNGL(DROT(J,1)*DPV(1)+DROT(J,2)*DPV(2)+ + DROT(J,3)*DPV(3)) 110 CONTINUE 120 CONTINUE RETURN END C C ------------------------------------------------------------------- C SUBROUTINE DBOOST3(BEX,BEY,BEZ) C C *** ROUTINE TO MAKE BOOSTS IN DOUBLE PRECISION. C include 'rboost.inc' C DOUBLE PRECISION BEX,BEY,BEZ,GA,BEP,GABEP C ibskip = 0 if((1.D+0-BEX**2-BEY**2-BEZ**2).lt.0) then ibskip = 1 return endif GA=1./DSQRT(1.D+0-BEX**2-BEY**2-BEZ**2) DO 140 I=1,4 BEP=BEX*DBLE(PA(I,1))+BEY*DBLE(PA(I,2))+BEZ*DBLE(PA(I,3)) GABEP=GA*(GA*BEP/(1.D+0+GA)+DBLE(PA(I,4))) PA(I,1)=SNGL(DBLE(PA(I,1))+GABEP*BEX) PA(I,2)=SNGL(DBLE(PA(I,2))+GABEP*BEY) PA(I,3)=SNGL(DBLE(PA(I,3))+GABEP*BEZ) PA(I,4)=SNGL(GA*(DBLE(PA(I,4))+BEP)) 140 CONTINUE RETURN END C C ------------------------------------------------------------------- C SUBROUTINE ROTATE3(DTH,DPH) C C *** ROUTINE TO ROTATE IN DOUBLE PRECISION. C include 'rboost.inc' C DIMENSION DROT(3,3),DPV(3) DROT(1,1)=DCOS(DTH)*DCOS(DPH) DROT(1,2)=-DSIN(DPH) DROT(1,3)=DSIN(DTH)*DCOS(DPH) DROT(2,1)=DCOS(DTH)*DSIN(DPH) DROT(2,2)=DCOS(DPH) DROT(2,3)=DSIN(DTH)*DSIN(DPH) DROT(3,1)=-DSIN(DTH) DROT(3,2)=0.D+0 DROT(3,3)=DCOS(DTH) DO 120 I=1,4 DO 100 J=1,3 100 DPV(J)=DBLE(PA(I,J)) DO 110 J=1,3 PA(I,J)=SNGL(DROT(J,1)*DPV(1)+DROT(J,2)*DPV(2)+ + DROT(J,3)*DPV(3)) 110 CONTINUE 120 CONTINUE RETURN END C C ------------------------------------------------------------------- C REAL FUNCTION USERRAN(iseme) C C *** USER RANDOM NUMBER GENERATION FUNCTION C userran = ran(iseme) C RETURN END C