HOW TO RUN: FOURJPHACT A. Ballestrero Vers 0.2 (March 2000) In this note flags, parameters and options of the program FOURJPHACT are examined. Other explanations can be found in the file readme.ex in the subdirectory /example The program is intended for studies of e+e- -> four jets (qcd) and in particular for interfacing massive 4 parton matrix elements with PYTHIA and its routine PY4JET. For this use the program must be linked to PYTHIA version 6.137 or higher. Input Parameters The parameters ppmax, irad, icom for the call to Pythia subroutine py4jet are set at 0.d0, 0, 0 in a data statement in fxn function. They can of course be changed if needed Other parameters like quark masses, coupling constants etc. are given in the data statement of the main program and one may eventually change them. The quark masses have been fixed to the same value as in PYTHIA Flags and inputs The possibilities given by the program as well as its use can be understood by an analysis of the read statements and of the possible answers to them. We report and comment them below Every parameter whose initial is i or n is integer*4. All others are real*8. When a variable has a yes/no option the value 1 corresponds to YES, 0 to NO. All energies and momenta must be expressed in GeV. Read statements: read*, e_cm * Energy in GeV read*, ioneshot ! yes/no all final states at once * ioneshot=0 for choosing only one particular final state * ioneshot=1 for computing all final state at once if (ioneshot.eq.1) then read*, icompprob !yes/no compute first every channel's probab. * icompprob=1 : all 20 final states are computed one by one * for determining the cross sections for the * chosen energy and cuts. The corresponding * probabilities of the various final states * are then stored in the file abprob.dat * In it, the nth entry represents the sum of * the probabilities of the first n final * states. These probabilities are then used * to compute all final states at once. * icompprob=0 : the program assumes that in the same * directory is present the file abprob.dat, * so that it just computes all final states in * one shot. else read*, if1 ! flavour index 3 4 (1=d,2=u,3=s,4=c,5=b,6=t) * if oneshot is 0 one has to choose the final state, * and this is done here by choosing first the type of * quarks corresponding to particles 3 and 4 read*, if2 !flavour index 5 6 (1=d,2=u,3=s,4=c,5=b,6=t,21=gg) * here the identity of particles 5 and 6 is chosen. They * can be quarks of various flavour or gluons endif read*, ialfasrun ! * ialfasrun=0 the value of alpha strong is taken from data * statement ( where it can be modified if * needed) * ialfasrun=1 alpha strong running is used. For each event * the alpha_s(M) is computed and used. * The scale M is chosen to be the invariant * mass of the two gluons for final states * qq~gg, the invariant mass of the gluon * propagator for each diagram of 4q processes. READ*,isr ! yes/no ISR * isr=0 no initial state radiation * isr=1 collinear initial state radiation implemented via * structure funcion method. READ*,iycut ! iycut=0 no ycut; =1 durham; =2 jade; =3 cambridge * ycut at parton level can be implemented according to * Durham (iycut=1), Jade (iycut=2) or Cambridge (iycut=3) * algorithm. If iycut=0 is chosen no ycut at parton level * is requested. IF (iycut.gt.0) then READ*, ycut ENDIF * the actual value of ycut is specified here if iycut=1 has * been chosen READ*,icut ! yes/no cuts * Cuts specified in the following may ( icut=1) or not * (icut=0}) be implemented. If icut=1 all default cuts * of the list must be specified. IF(icut.EQ.1)THEN READ*,e_min ! 4 energy lower cuts (GeV) READ*,e_max ! 4 energy upper cuts (GeV) * e_min and e_max correspond to the 4 lower and upper * energies for particle 3, 4, 5, 6 respectively. READ*,rm_min ! 6 invariant mass lower limits (GeV) ! (34, 35, 36, 45, 46, 56) READ*,rm_max ! 6 invariant mass upper limits (GeV) * rm_min and rm_max are the 6 invariant mass lower and * upper limits respectively. They must be given in the * following order: * m(34), m(35), m(36), m(45), m(46), m(56). READ*,pt_min ! 4 transverse momenta lower cuts (GeV) READ*,pt_max ! 4 transverse momenta upper cuts (GeV) * pt_min and pt_max are the 4 lower and upper values of * the transverse momenta. The order is as before 3, 4, 5, 6 READ*,icos ! angular cuts in deg (0) or cos (1) * icos = 1 implies that the following angular cuts must be * expressed in terms of the cosines of the angles. With * icos = 0 one must instead specify the angles in degrees. READ*,thbeam_min! 4 particle-beam angle lower (in degrees) cuts READ*,thbeam_max! 4 particle-beam angle upper (in degrees) cuts * thbeam_min and thbeam_max are the 4 lower and the * 4 upper limits for the angle that particles 3,4,5,6 * produce with the beam (e^+). READ*,thsep_min ! 6 particle-particle angle lower (in degs) cuts READ*,thsep_max ! 6 particle-particle angle upper (in degs) cuts * thsep_min and thsep_max are the 6 lower and the 6 * upper limits for particle-particle angles: the order is * again (3 4), (3 5), (3 6), (4 5), (4 6), (5 6). ENDIF ************************************************************* c this part has to be repeated in input file if ioneshot=1 c and icompprob =1 has been chosen ************************************************************* * As indicated in the above comment from the code, the * following read statements have to be answered twice if * ioneshot=1 and icompprob =1 has been chosen. The first * time the input is used in the calculation of the various * final states one by one, the second in that of all final * states simultaneously. One may of course vary the input * for the two steps. For instance one may require a * higher number of calls or unweighted events generation * only for the second step where all processes are * computed simultaneously. * * When ioneshot=0 only one channel will be computed. * When ioneshot=1 and icompprob=0 all channels in one shot * are computed. In both cases only one cross section is * computed and one specifies only the inputs for this one. 100 READ*,idistr ! yes/no distributions * Distributions at the parton level can be implemented * very easily. If one wants to use this possibility the * flag idistr must be =1, =0 otherwise. * In the program there is the line: * ' include * abdis.dis' * which must be uncommented before compiling if one chooses * idistr=1. In this case the file abdis.dis must be written * for implementing distributions, as in the following * example: * * string(1)=' Distribution: bb~ invariant mass' * distr_var(1)=sqrt((p3(0)+p4(0))**2-(p3(1)+p4(1))**2 * & -(p3(2)+p4(2))**2-(p3(3)+p4(3))**2) * string(2)='Distribution: Charged lepton energy ' * distr_var(2)=p5(0) * * In it, the title of the ith distribution is given in the * character string(i) and the ith quantity to be * distributed in bins, computed from the particle momenta, * is assigned to distr_var(i). * The resulting cross sections, corresponding to every * single bin will be stored in the file abdis.dat, where * each line will contain 3 numbers: the value of the * central point of the bin, the distribution for the bin * (cross section divided by the width of the bin, in order * to reproduce d sigma/dx for a distribution of the * variable x), and the estimated statistical error. IF(idistr.EQ.1)THEN READ*,ndistr ! number of distributions * ndistr is the number of distributions defined in abdis.dis DO i=1,ndistr READ*,nsubint(i) ! number of sub-intervals with different ! binning READ*,(distr_estrinf(i,j),j=1,nsubint(i)+1) !lower limits of ! each subint (which coincide with the upper limit !of the previous one) + upper limit of the last subint READ*,(nbin_number(i,j),j=1,nsubint(i)) !number of bins in !each subint ENDDO !i ENDIF * For each distribution i, one must specify: * nsubint(i), the number of sub-intervals with different * binning in the ith distribution (=1 when all bins are * of the same length). The subintervals must be * contiguous. * (distr_estrinf(i,j),j=1,nsubint(i)+1), the lower limits * of each subinterval (which coincide with the upper * limit of the previous one as they must be contiguous) * and, as last entry, the upper limit of the last * subinterval. In case all bins are of the same length, * this corresponds only to the lower and upper limit of * the interval for the distribution. * (nbin_number(i,j),j=1,nsubint(i)), the number of bins in * each subinterval. In case all bins are of the same * length, this corresponds only to the number of bins READ*,iflat ! yes/no flat event generation * One may choose ( iflat=1) or not ( iflat =0) to generate * unweighted events. In the first case, the number of * iterations ( itmx) which must be specified in the * following must be 2. The integration routine will * perform the requested number of iterations for * thermalization (see iterm, ncall_term,itmx_term below) and * then the two iterations in which the integral is * evaluated. In the first iteration the maximum for the * hit-or-miss procedure will be determined and used in the * second iteration where the unweighted generation will * take place. After the run the output file will report as * usual the result of the integration and its error. * It will also report the maximum used for the hit-or-miss * procedure, the maximum found in the second iteration, and * the number of events which were greater than the maximum * used. There is also the possibility to repeat the * generation just starting directly from the second * iteration. This might be useful if too many events * exceeded the maximum chosen. Another useful possibility * is to rerun the program asking for a predetermined number * of unweighted events to generate. So one may for instance * run the first time just to find an appropriate maximum * and a correct cross section (and eventually to produce * partonic distributions). One controls then if it is * better to multiply the maximum for an appropriate factor * and with a second run one generates the desired number * of events. To save time, the hadronization procedure with * the link to Pythia may be activated only in this second * run. IF(iflat.EQ.1)THEN READ*,scalemax ! scale factor for the maximum * scalemax is the coefficient by which the maximum of the * first iteration can be multiplied, in order to vary the * efficiency of the hit-or-miss procedure or in order to * avoid values exceeding the maximum. READ*,istorvegas ! yes/no VEGAS data stored * VEGAS data are (if istorvegas=1) or not (if istorvegas=0) * stored in ABVEGAS.DAT after the first iteration. Stored * Vegas data are necessary if one wants to rerun the program * to generate again unweighted events. When the program is * rerun using Vegas data stored, the maximum of the second * iteration will be automatically used as the new maximum. READ*,irepeat ! 0=normal;1=repeat 2nd only;2=rep. with ! nflevts fixed ! 1=repeat only second iteration using vegas ! data stored in previous run, ! 2=repeat, generating an exact number ! nflevts of flat events using vegas data ! stored in previous run * irepeat has to be set to 0 for the first run * It has to be set to 1 if one wants to rerun exactly with * the same input starting from the second iteration. In this * case the same weighted points (if scalemax is not varied) * will be reproduced. * irepeat=2 has to be chosen if one wants to rerun with * the same imput and grid as before, but letting the * program run until a requested number of events nfltevts * is reached. For both cases irepeat=1 and 2 one might of * course vary scalemax, ipyth or istormom. IF(irepeat.eq.2)THEN READ*, nflevts ! number of events to be generated as ! described above END IF READ*,istormom ! yes/no momenta of flat events written ! in .dat files * The momenta of the unweighted events are written in * ABMOM.DAT file if istormom=1, they are not written if * istormom=0. READ*,ipyth ! yes/no call to Pythia * Every unweighted event is passed to the standard COMMON * HEPTV if ipyth=1 and Pythia subroutine PY4JET is then * called. Pythia is not used if ipyth=0. ENDIF * Parameters that specify how the integration will be performed by * VEGAS: READ*,acc ! integration accuracy * acc is the integration accuracy. When this accuracy is * reached after a certain integration iteration, the * remaining iterations are not performed. READ*,iterm ! yes/no thermalization * If iterm=1 a certain number of integration iterations * ( =itmx_term) are used only for adapting the integration * grid. * Their result is not used for the final integral. Each * thermalizing iteration makes use of a maximum of * ncall_term evaluations of the amplitude. If iterm=0 * these iterations are skipped. One or two itmx_term with * few ncall_term are often useful. READ*,ncall_term ! thermalization calls per iteration READ*,itmx_term ! thermalization iterations READ*,ncall ! integration calls per iteration READ*,itmx ! integration iterations ! (2 for flat event generation) * itmx is the maximum number of iterations used to evaluate * the integral. Each iteration makes use of a maximum of * ncall evaluations of the amplitude. Normally it is better * not to use more than about five iterations. If higher * precision is requested it is convenient to increase ncall * and not itmx. * As a final remark about the choiche of these parameters, * one must be aware of the fact that final results with a * chi^2 much greater than the number of iterations are * not to be trusted. When this happens, one has to increase * ncall.