C...Sample main program and subroutines illustrating how to set up C...an external process in Pythia. C...Example: q + g -> q + gamma at the Tevatron. This process C...already exists in standard Pythia, som comparisons are possible C...between the two implementations. C...Reminder: since a dummy copy of the PYUPEV routine is included C...in the Pythia standard distribution (to avoid problems with C...unresolved external references) it is necessary to remove C...(the linking of) this dummy to avoid clashes with the routine C...below. C----------------------------------------------------------------- C...Preamble: declarations. C...All real arithmetic in double precision. IMPLICIT DOUBLE PRECISION(A-H, O-Z) C...Three Pythia functions return integers, so need declaring. INTEGER PYK,PYCHGE,PYCOMP C...EXTERNAL statement links PYDATA on most machines. EXTERNAL PYDATA C...Commonblocks. C...The event record. COMMON/PYJETS/N,NPAD,K(4000,5),P(4000,5),V(4000,5) C...Parameters. COMMON/PYDAT1/MSTU(200),PARU(200),MSTJ(200),PARJ(200) C...Particle properties + some flavour parameters. COMMON/PYDAT2/KCHG(500,4),PMAS(500,4),PARF(2000),VCKM(4,4) C...Decay information. COMMON/PYDAT3/MDCY(500,3),MDME(4000,2),BRAT(4000),KFDP(4000,5) C...Selection of hard scattering subprocesses. COMMON/PYSUBS/MSEL,MSELPD,MSUB(500),KFIN(2,-40:40),CKIN(200) C...Parameters. COMMON/PYPARS/MSTP(200),PARP(200),MSTI(200),PARI(200) C...Information on user-selected processes. COMMON/PYUPPR/NUP,KUP(20,7),NFUP,IFUP(10,2),PUP(20,5),Q2UP(0:10) C...The user's own transfer of information. COMMON/MYCOMM/ECM,PTMIN,INIT C...Counter for number of generated events of each type. DIMENSION NCOUNT(2) DATA NCOUNT/2*0/ C----------------------------------------------------------------- C...First section: initialization. C...Number of events and cm energy. NEV=1000 ECM=1800D0 C...Define a pTmin cutoff for the process. PTMIN=10D0 CKIN(3)=PTMIN C...Estimate reasonable maximum of the cross section. C...(This is done by a cheating call to PYUPEV with special switch.) INIT=1 CALL PYUPEV(ISUB,SIGMAX) INIT=0 C...Use ISUB code 200 for external process q + g -> q + gamma. CALL PYUPIN(200,'q + g -> q + gamma',SIGMAX) C...Switch on both internal and external process (for comparison). MSEL=0 MSUB(29)=1 MSUB(200)=1 C...Switch off unnecessary aspects: initial and final state C...showers, multiple interactions, hadronization. C...(Optional for faster simulation of the parton-level C...processes only.) C MSTP(61)=0 C MSTP(71)=0 C MSTP(81)=0 C MSTP(111)=0 C...Initialize. CALL PYINIT('CMS','p','pbar',ECM) C...Book histograms. CALL PYBOOK(1,'pT spectrum proc 29',100,0D0,100D0) CALL PYBOOK(2,'pT spectrum proc 200',100,0D0,100D0) CALL PYBOOK(3,'n_ch distribution proc 29',100,-1D0,199D0) CALL PYBOOK(4,'n_ch distribution proc 200',100,-1D0,199D0) C----------------------------------------------------------------- C...Second section: event loop. C...Generate events and look at first few. DO 200 IEV=1,NEV CALL PYEVNT ISUB=MSTI(1) ICASE=1 IF(ISUB.EQ.200) ICASE=2 NCOUNT(ICASE)=NCOUNT(ICASE)+1 IF(NCOUNT(ICASE).LE.1) THEN WRITE(6,*) ' Following event is subprocess',ISUB CALL PYLIST(1) ENDIF C...Read out generated pT distribution and charged multiplicity. PTGEN=PARI(17) CALL PYFILL(ICASE,PTGEN,1D0) CALL PYEDIT(3) CALL PYFILL(ICASE+2,DBLE(N),1D0) C...End of loop over events. 200 CONTINUE C----------------------------------------------------------------- C...Third section: produce output and end. C...Cross section table. CALL PYSTAT(1) C...Histograms. CALL PYHIST END C********************************************************************* C...This routine samples the phase space, evaluates the process cross C...and sets the colour topology and parton shower scales. SUBROUTINE PYUPEV(ISUB,SIGEV) C...All real arithmetic in double precision. IMPLICIT DOUBLE PRECISION(A-H, O-Z) C...Cross-section information. COMMON/PYINT5/NGENPD,NGEN(0:500,3),XSEC(0:500,3) C...Information on user-selected processes. COMMON/PYUPPR/NUP,KUP(20,7),NFUP,IFUP(10,2),PUP(20,5),Q2UP(0:10) C...The user's own transfer of information. COMMON/MYCOMM/ECM,PTMIN,INIT C...Local arrays and parameters. DIMENSION XPP(-25:25),XPPBAR(-25:25),TERM(20) DATA PI/3.141592653589793D0/ DATA CONV/0.3894D0/ C----------------------------------------------------------------- C...Default return values. ISUB=200 SIGEV=0D0 C...Default flag for event rejection. IREJ=0 C----------------------------------------------------------------- C...The first step is to decide how the phase space should be sampled. C...The total cross section for a 2 -> 2 process can be written C... sigma = Integral d(x_1) Integral d(x_2) Integral d(t-hat) C... f_1(x_1,Q^2) f_2(x_2,Q^2) d(sigma-hat)/d(t-hat) C...(see manual section 7.2), so a simple recipe would be to pick C...x_1, x_2 and t_hat uniformly. However, the matrix element and C...parton distributions are peaked, so this would be inefficient. C...In particlar, the matrix element is proportional to 1/u-hat, C...so note that d(t-hat)= d(u-hat), and therefore rewrite as C... sigma = Integral d(x_1)/x_1 Integral d(x_2)/x_2 C... Integral d(u-hat)/u-hat x_1 f_1(x_1,Q^2) C... x_2 f_2(x_2,Q^2) u-hat d(sigma-hat)/d(t-hat). C...(More fancy methods are used in Pythia to improve efficiency, C...but here the principles are more important.) C...The pTmin cut translates into cuts on x_1, x_2 and u-hat. X1MIN=4D0*PTMIN**2/ECM**2 X2MIN=X1MIN UHTMIN=PTMIN**2 UHTMAX=ECM**2 IF(INIT.EQ.0) THEN C...Normally select a phase space point at random from C...d(x_1)/x_1 d(x_2)/x_2 d(u-hat)/u-hat. X1=X1MIN**PYR(0) X2=X2MIN**PYR(0) UHAT=-UHTMIN*(UHTMAX/UHTMIN)**PYR(0) ELSE C...For special initialization call pick lowest point as C...likely to represent maximum cross section. X1=2D0*PTMIN/ECM X2=X1 UHAT=-2D0*PTMIN**2 ENDIF C...Derive other kinematical quantities. SHAT=X1*X2*ECM**2 THAT=-SHAT-UHAT PT2=THAT*UHAT/SHAT C...By decoupling the x_1, x_2 and u-hat variables we overestimate C...the phase space region. We therefore have to reject events C...outside allowed phase space by putting the cross section = 0. C...However, to keep bookkeeping correct, one still needs to set C...variables as if there is an event, e.g. picking the same C...point as for the initialization call. IF(INIT.EQ.0.AND.PT2.LT.PTMIN**2) THEN IREJ=1 X1=2D0*PTMIN/ECM X2=X1 UHAT=-2D0*PTMIN**2 SHAT=X1*X2*ECM**2 THAT=-SHAT-UHAT PT2=THAT*UHAT/SHAT ENDIF C...The selection strategy for the phase space point relies on C...a phase space volume V inside which points are sampled. The C...expression sigma = Integral dV d(sigma)/dV = C... = < d(sigma)/dV > * Integral dV C...(with < ... > denoting (Monte Carlo sampled) mean value over V) C...implies we need Integral dV = C...Integral d(x_1)/x_1 Integral d(x_2)/x_2 Integral d(u-hat)/u-hat. PHSPV=(-LOG(X1MIN))*(-LOG(X2MIN))*LOG(UHTMAX/UHTMIN) C----------------------------------------------------------------- C...The second step is to evaluate the cross section in the C...phase space point selected above. C...Pick Q2 scale (which involves some arbitrariness) C...and evaluate alpha_em and alpha_s. Q2=PT2 AEM=PYALEM(Q2) ALPS=PYALPS(Q2) C...Evaluate matrix element d(sigma-hat)/d(t-hat), except for C...a factor e_q^2 that varies over the quark flavours. DSIGDT=(PI*AEM*ALPS/3D0)*(SHAT**2+UHAT**2)/(-SHAT**3*UHAT) C...Correction factor u-hat from our choice of phase space sampling. DSIGDT=(-UHAT)*DSIGDT C...Now need sum of charge-square-weighted quark*gluon distributions. C...The relative size of these terms is also needed for the C...Monte Carlo selection of one specific flavour in initial state, C...and therefore the terms are stored in an array. C...Parton distributions (multiplied by x) of p and pbar. CALL PYPDFU(2212,X1,Q2,XPP) CALL PYPDFU(-2212,X2,Q2,XPPBAR) C...Loop over quark flavours (up to b). SUM=0D0 DO 100 IQ=1,5 C...Squared quark charge. EQ2=1D0/9D0 IF(IQ.EQ.2.OR.IQ.EQ.4) EQ2=4D0/9D0 C...q from p, g from pbar. TERM(4*IQ-3)=EQ2*XPP(IQ)*XPPBAR(21) C...qbar from p, g from pbar. TERM(4*IQ-2)=EQ2*XPP(-IQ)*XPPBAR(21) C...g from p, q from pbar. TERM(4*IQ-1)=EQ2*XPP(21)*XPPBAR(IQ) C...g from p, qbar from pbar. TERM(4*IQ)=EQ2*XPP(21)*XPPBAR(-IQ) SUM=SUM+TERM(4*IQ-3)+TERM(4*IQ-2)+TERM(4*IQ-1)+TERM(4*IQ) 100 CONTINUE C...Now possible to define the SIGEV return value. C...The CONV factorconverts from GeV**(-2) to mb. IF(IREJ.EQ.0) SIGEV=CONV*PHSPV*DSIGDT*SUM C...Done when initializing for maximum cross section. IF(INIT.EQ.1) RETURN C----------------------------------------------------------------- C...The third step is to set up the partonic process that is selected. C...It is always to be defined in the cm frame of the two incoming C...particles (which here coincides with the "lab" frame). C...Select one of the possible terms at random. RANSUM=PYR(0)*SUM ISEL=0 110 ISEL=ISEL+1 RANSUM=RANSUM-TERM(ISEL) IF(ISEL.LT.20.AND.RANSUM.GT.0D0) GOTO 110 C...Read off quark flavour and quark side from selection. KFLQ=(ISEL+3)/4 IF(ISEL-4*KFLQ.EQ.-2.OR.ISEL-4*KFLQ.EQ.0) KFLQ=-KFLQ ISIDE=1 IF(ISEL-4*KFLQ.GE.-1) ISIDE=2 C...Define number of partons involved and reset to zero C...(except status code). NUP=4 DO 130 I=1,NUP KUP(I,1)=1 DO 120 J=2,7 120 KUP(I,J)=0 DO 130 J=1,5 130 PUP(I,J)=0D0 C...Flavour codes for entries. Note that definition of t-hat C...means quark either is 1 and 3 or 2 and 4. KUP(ISIDE,2)=KFLQ KUP(3-ISIDE,2)=21 KUP(ISIDE+2,2)=KFLQ KUP(5-ISIDE,2)=22 C...Masses of final state entries (initial assumed massless). PUP(3,5)=PYMASS(KUP(3,2)) PUP(4,5)=PYMASS(KUP(4,2)) C...Four-momenta of the incoming partons simple. PUP(1,4)=0.5D0*X1*ECM PUP(1,3)=PUP(1,4) PUP(2,4)=0.5D0*X2*ECM PUP(2,3)=-PUP(2,4) C...Energies and absolute momentum of the outgoing partons in C...the subsystem frame. RTSHAT=SQRT(X1*X2)*ECM PUP34=0.5D0*(RTSHAT**2+PUP(3,5)**2-PUP(4,5)**2)/RTSHAT PUP44=RTSHAT-PUP34 PABS=0.5D0*SQRT(MAX(0D0,(RTSHAT**2-PUP(3,5)**2- &PUP(4,5)**2)**2-4D0*PUP(3,5)**2*PUP(4,5)**2))/RTSHAT C...Subsystem scattering angle defined neglecting quark mass. COSTHE=(THAT-UHAT)/SHAT SINTHE=SQRT(MAX(0D0,1D0-COSTHE**2)) C...Azimuthal angle at random. PHI=2D0*PI*PYR(0) C...Momenta of outgoing partons in the subsystem frame. PUP(3,1)=PABS*SINTHE*COS(PHI) PUP(3,2)=PABS*SINTHE*SIN(PHI) PUP33=PABS*COSTHE PUP(4,1)=-PUP(3,1) PUP(4,2)=-PUP(3,2) PUP43=-PUP33 C...Longitudinal boost of outgoing partons to cm frame. BETA=(X1-X2)/(X1+X2) GAMMA=0.5D0*(X1+X2)/SQRT(X1*X2) PUP(3,3)=GAMMA*(PUP33+BETA*PUP34) PUP(3,4)=GAMMA*(PUP34+BETA*PUP33) PUP(4,3)=GAMMA*(PUP43+BETA*PUP44) PUP(4,4)=GAMMA*(PUP44+BETA*PUP43) C...Store colour flow. IF(KFLQ.GT.0) THEN C...Colour stretched from initial quark to gluon to final quark. KUP(ISIDE,6)=3-ISIDE KUP(3-ISIDE,7)=ISIDE KUP(3-ISIDE,6)=ISIDE+2 KUP(ISIDE+2,4)=3-ISIDE ELSE C...Ditto but with colour <--> anticolour. KUP(ISIDE,7)=3-ISIDE KUP(3-ISIDE,6)=ISIDE KUP(3-ISIDE,7)=ISIDE+2 KUP(ISiDE+2,5)=3-ISIDE ENDIF C...Pick Q2 scales for initial state shower evolution. Q2UP(0)=Q2 C...Allow final state shower of q and gamma. NFUP=1 IFUP(1,1)=3 IFUP(1,2)=4 Q2UP(1)=Q2 C...Done. RETURN END