Home

QCDNUM16: A fast QCD evolution program

1. C TTTTTTLETTITTTT 7 o a eV a NE KA b 5 f 4 f 5 Moco Nocp oO 0 o OV C3 EX CY UC CX 9 rrmqeemmmPTIPTTETITTTITTITT TT O o NOUO B YD e prm Figure 1 a The evolution of a in NLO in the MS scheme starting from a M 0 110 0 120 0 130 with the quark mass thresholds taken to be me 1 5 5 GeV b The corresponding values of the QCD scale parameter A to evolve down to m with f 4 flavors and so on In this way a is continuous at the flavor thresholds but its derivative is not Fig la shows the evolution of a for three input values of a M2 0 110 0 120 and 0 130 The quark mass thresholds were set to me mp 1 5 5 GeV The values of the QCD scale parameter A as defined in eq 2 are shown in fig 1b Notice that A is discontinuous at the flavor thresholds Notice also from eq 2 that a cannot be evolved to values of u at or below the scale A for instance not below u 0 44 GeV if a M 0 130 and if the flavor thresholds are taken to be those given above 2 2 QCD evolution of parton densities The Altarelli Parisi evolution equations 6 can be written as oun ae S p Pig a ln u may Here u u is the mass factorization scale h x u are the parton density distributions P x are the QCD splitting functions and f is the number of active flavors When the index i gt 0 h denotes the quark distribution q x
2. The patchy directive tuse gheavy t exe makes the heavy quark code available This directive may be omitted or commented out producing a smaller executable but then the heavy quark structure functions cannot be calculated anymore 4 The Qcdnum code is written in standard Fortran77 and uses only a few CERN library routines namely DGAUSS D103 DDILOG C304 FLPSOR M103 and for character string manipulation LENOCC M507 and CHPACK M432 Qcdnum does not use a dynamic memory manager so that all common blocks are defined at compilation time The size of these common blocks is governed by the two parameters mxx and mq2 mxx 1 and mq2 1 are the maximum number of grid points in x and Q the user can define This can be set directly in the source code or through the patchy directive see table 1 e g 14 An attempt to do so will cause Qcdnum to abend with an error message 17 bin csh fx cern pro bin ypatchy qcdprog tty qcdprog list go lt lt OPTION MAPASM UREF USE DOUBLE USE P QUCOMM QUPROG T EXE user program USE QCDCOM QCDNUM QCUTIL T EXE Qcdnum code TUSE HEAVY TSEXE 54 RR Heavy quark code TREPL QCDCOM QCDCOM 6 7 lees Max grid points PARAMETER MXX 100 PARAMETER MQ2 40 EXE P CRA PAM 11 T ATTACH T CARDS qcduser car PAM 12 T ATTACH T CARDS qcdnumi6 car QUIT 77 f 77ldflags 0 qcdprog f o qcdprog L cern pro lib l
3. Doing so invalidates the current weight tables An attempt to use them afterwards in an evolution or structure function routine will cause a fatal error 4 4 4 4 4 var typl deflt value description 4 4 4 4 4 W1ANA L T T Analytical LO weight calculation WiNUM L F F Numerical LO weight calculation W2NUM L T T Numerical NLO weight calculation W2STF L1 T T Structure function NLO weights WTF2C L F F F2 charm weight calculation WTF2B L F F F2 bottom weight calculation WFLC LI F F FL charm weight calculation WTFLB L F F FL bottom weight calculation 4 4 4 4 4 Table 4 Logical variables which steer the weight calculations 2f The CPU time needed for the calculation of the light quark weights is modest about 60 seconds for 100 x grid points but will increase quadratically with the number of grid points The logical variables WTF2C WTF2B WTFLC and WTFLB enable the calculation of weight tables needed for the heavy flavor structure functions These tables depend on the definition of both the x and the Q grid which should not be changed after the call to QNFILW They also depend on the charm or bottom mass which should not be changed either The calculation of the hea
4. Some details on how to control accuracy and speed are given in section 5 whereas section 6 reference section contains a short description of all subroutine calls 2 The QCD analysis of structure functions 2 1 QCD evolution of a The renormalization scale u u dependence of the strong coupling constant a is governed by the renormalization group equation which reads in next to leading order NLO Ja5 11 In we ga u z fias u 1 where a as 4r and the QCD beta functions are given by fy 11 2 3 and 3 102 38f 3 with f the number of active quark flavors The first entry in the Qcdnum history record dates from January 1987 4 The solution of eq 1 can be written as I B n fas u T ra as u fs nous Pol 55 2 where we have introduced the QCD scale parameter A which is defined here as the scale where the left hand side of eq 2 vanishes There are however other definitions of A which are widely used we prefer to remove this ambiguous scale parameter by writing eq 2 in terms of the value of a at some input scale nos Lo tg HP Pu fas balu mu a S 5 i eri D S The beta functions depend on the number of flavors f 3 4 5 The flavor thresholds are usually set at u2 m2 Thus evolving downwards from a M7 for instance Qcdnum calculates a m from eq 3 with f 5 flavors This value of a is then taken as an input 25 E 9 Q MZ 0 110
5. this change being calculated at order n 1 Clearly if n is large enough so that AF P AF we have FO Q i3 e F 34 AFP Q 2 FP Qn so that the u2 dependence vanishes for large n From eq 23 we obtain for the scale dependence of the structure functions in LO see eqs 18 21 M FA in vasi 24 Thus in LO one takes for the structure functions F5 and xF the corresponding quark distri bution at 2 instead of at Q For the scale dependence in NLO we get FO Q 2 2 O s us co PO Q 2 25 23 Q uu salum ka T n im Q zq uy 25 FP Qu FP Notice that Fr when calculated in NLO still has a large LO type scale uncertainty Next we set u Q and vary p13 This leads to a change in the value of a being taken at u instead of Q and to a modification of the NLO splitting functions Oq Q LR as UR 0 as u Bo oln Q 2T n Y Ua pu Apam B do e The structure functions are calculated from the evolved parton distributions with 2 F Q uR 1 0kr ma Q we aar amp zq Q uR 27 The scale uncertainties are usually assigned to o and the parton distributions These uncer tainties are obtained by varying the scale and re determining a and the parton densities in a fit where the structure functions are kept fixed For determinations of the scale errors on a see 2 11 and 12 8To keep the notation simple we give in this section only the expres
6. PU f caq ci x PaG2 P amp 2 c f PPL2 P x f PMI2 PP a f Table 9 Selection options in QNSPLF 41 The input character variable opt selects the splitting or coefficient function as listed in table 9 As indicated in this table some functions depend on z only in which case the value of q2 and nf are irrelevant whilst others depend in addition on either f then q2 is irrelevant or Q then nf is irrelevant The function QNSPLF is not affected by grid definitions and may thus be called at any time after QNINIT Several splitting coefficient functions are not defined for x 1 only their integral over an interval containing x 1 as an upper limit these functions are set to zero at x 1 4 17 Qcdnum error handling Qcdnum is robust or at least attempts to be it should either give the correct answer numer ically or grind to a halt printing an error message An example of such a message is given in table 10 here we have attempted to calculate the proton F structure function outside the grid boundaries or cuts Fatal errors might have their cause upstream For example if weights are calculated and the x Q grid is changed afterwards the weight tables which are defined on the grid are invalidated The first call to a routine which uses them will cause an abend complaining that the weights are not available To avoid such errors we recommend to stick to the flow of action as described in section 4 2 In
7. W 6 W 7 W8 W 9 W 10 ID NAME nf SING UMIN SPLU CPLU BPLU XQVA FREE FREE FREE FREE 4 4 O GLUON 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 1 SINGL 1 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 2 UMIND 0 00 1 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 3 SPLUS 0 00 0 00 1 00 0 00 0 00 0 00 0 00 0 00 0 00 0 00 4 CPLUS 0 00 0 00 0 00 1 00 0 00 0 00 0 00 0 00 0 00 0 00 5 BPLUS 0 00 0 00 0 00 0 00 1 00 0 00 0 00 0 00 0 00 0 00 6 XQVAL 0 00 0 00 0 00 0 00 0 00 1 00 0 00 0 00 0 00 0 00 4 Re ee a Se a a r 11 PROTO 3 0 22 0 17 0 17 0 00 0 00 0 00 0 00 0 00 0 00 0 00 4 0 28 0 17 0 17 0 17 0 00 0 00 0 00 0 00 0 00 0 00 5 0 24 0 17 0 17 0 17 0 17 0 00 0 00 0 00 0 00 0 00 4 Table 5 Example of a list of parton distributions 3l 4 10 QCD evolution parameters and options Before evolving parton distributions and or calculating structure functions a few parame ters options can be set as follows 1 The order of all calculations in Qcdnum is governed by the value of iorder 1 LO or 2 NLO The default is set to 2 NLO which can be changed at any time by call QNISET ORDER 1 and vice versa 2 The QCD evolution of a is internally calculated by solving the renormalization group equation eq
8. a cut is applied the region outside it becomes inaccessible to the user GRCUTS can be called at any time to release one or more cuts or not apply it at all the corresponding parameter should be set to a value lt 0 The function ifail IFAILXQ x q2 can be used to investigate if a given z Q fails or passes the cuts 0 1 no yes fail QMINA cut see below 0 1 no yes fail ROOTS cut 0 1 no yes fail QMAX cut 0 1 no yes fail QMIN cut 0 1 no yes fail XMIN cut ifail ijklm B O wo e Thus ifail 1001 means that r Q is below the xmin cut and above the kinematic limit as defined by roots The corresponding function for the grid point ix iq is ifail IFAILIJ ix iq The QMINA cut mentioned above acts as a Q cut and is automatically set by Qcdnum min to the lowest Q grid point above the current value of A 1n QCD fits for instance you can set roots to the highest center of mass energy of your datasets e g 300 GeV if HERA data are included 33 4 11 QCD evolution of parton distributions The QCD evolution of parton distributions consists of two steps 1 Input of the parton distribution for all grid points in x at some fixed value of Q i e at a fixed Q grid point iq0 2 QCD evolution over the whole z Q grid or part of it when cuts are applied see sec tion 4 10 by calling one of the evolution routines described below The procedure is straight forward for the singlet gluon evoluti
9. in addition QADDSI cannot be called for the gluon and the singlet quark distribution 50 call QNPNUL name Set the patron distribution name to zero A fatal error occurs if name does not correspond to a memory resident distribution call EVOLSG iq0 iqmin iqmax Evolve the gluon and the singlet quark distribution over the range iqmin iqmax starting from the grid point iq0 A fatal error occurs if iqO lies outside the range iqmin iqmax This range might have been adjusted internally by Qcdnum if Q cuts are made For further details see section 4 11 call EVOLNP name iq0 iqmin iqmax As above but now evolve a non singlet xq or xA quark distribution A fatal error will occur if name does not correspond to a non singlet memory resident distribution For further details see section 4 11 call EVPLUS name iq0 iqmin iqmax As above but check that the number of flavors does not change in the range iqmin iqmax To be used for xq evolution which is discontinuous at the flavor thresholds For further details see section 4 11 val QPDFIJ name ix iq iflag Returns the value of the parton distribution name or linear combination at the grid point ix iq If ix or iq are outside the grid boundaries or cuts val is set to zero and iflag to 1 0 otherwise Qcdnum will then abend with an error message unless the f
10. of the user that the input argument name refers to a correct linear combination for instance to an object like the proton as defined in section 4 9 for f 3 flavors e Heavy quark structure functions can be calculated for all Q within the grid or cuts but might be numerically inaccurate for Q lt 1 5 GeV in this version of Qcdnum Attempts to calculate below Q 1 5 GeV result in a fatal error unless you switch this check off by call QNLSET CLOWQ false 4 13 The renormalization and factorization scale In Qcdnum the renormalization and mass factorization scales are defined as HR o aRUM be VON 2 HM amQ S bm 41 The default values for the parameters agr m and be m are given in table 6 notice that the mass factorization scale can be defined independently for the light and the heavy quark structure functions Qcdnum insists that all scales are within the Q grid boundaries and above the current value of A In addition the scale 4 2 is required to lie within the cuts if any In practice only one scale at the time is varied For instance one may identify the factorization scale 3 with Q and vary ww in the range Q 4 lt uh lt 4Q say This variation affects the value of a being taken at u instead of Q and also the evolution of the parton distributions through a modification of the NLO splitting functions In the second case 4 2 is identified with u3 which is varied in the range Q
11. points Each pass consisted of one singlet gluon evolution and four non singlet evolutions Furthermore F gt structure functions were calculated for 500 points randomly distributed below the kinematic limit Q xs with ys set to 300 GeV The amount of CPU time increases quadratically with the number of grid points in x and linearly with that in Q not shown but is still only about 8 minutes for 2500 evolutions and 25 x 10 structure function evaluations on a 150 x 60 grid Also shown in fig 3 is the effect of a ys cut which speeds up the evolutions by more than a factor of two fig 3b Only a quarter of the CPU time is spent when STFAST is used to calculate the structure functions instead of QSTFXQ fig 3c In practical cases however the gain will be less factor 2 typically since structure function data are usually ordered in bins of x and Q instead of being randomly distributed as in this example To summarize a good initial choice is about 100 grid points in z and about 40 in Q This gives an executable size of about 3 Mbyte an accuracy of well within 196 over a large kinematic range and execution times of the order of 10 minutes for several thousand evolutions and a few hundred thousand structure function evaluations 34If heavy quark structure functions are calculated the number of Q grid points should however not be taken too small 35In QCD fits you can still have the evolution on the whole z Q grid by releasing the cut o
12. slope Oh z u 801n u changes when crossing a flavor threshold This is also true for X q7 and A but not for qt this distribution is per definition discontinuous see eq 7 It is important to notice that Qcdnum can only evolve the singlet gluon and non singlet dis tributions Only such distributions should therefore be passed as input to the program Any linear combination of quark and anti quark distributions can be decomposed into a linear com bination of the singlet density X and one or more non singlet densities f f f it bi i bi 1 with c 2 5 and Ser As an example we take the proton distribution in charged lepton scattering which can be written as q x uel qx ai x 12 e D z gp x To be more precise one should pass the parton momentum densities xg xu xq etc instead of the parton number densities g used in this and the next sections 7 where e is the quark charge of flavor i in units of the electron charge and lt e gt 1 f Xe is the average of the square of the quark charges lt e gt 4 18 5 18 11 45 for f 3 4 5 flavors The distribution qg is a pure q type non singlet dp 2 J eia 2 13 f Because q and X are continuous functions of u it follows that the change in e must be compensated for by a discontinuity in qg when crossing a flavor threshold We find 1 d 4 df 3 4E 14 dU 9 du orzx In this way
13. spite of the error checking mechanism there is still room to go astray with a Qcdnum based analysis As a reminder we summarize below the most important do s and dont s described in the previous sections e QNINIT must be called before any other Qcdnum routine e When run in double precision mode floating point variables passed to Qcdnum should Input Opt F2 Name PROTO x 0 63000E 04 Q2 0 35000E 01 X Q2 or mu2 outside grid or cuts x 0 63000E 04 xmin 0 10000E 04 pass Q2 0 35000E 01 Qmin 0 75000E 00 pass Q2 0 35000E 01 Qmax 0 60000E 04 pass s 0 55556E 05 Smax 0 44100E 05 fail Q2 0 35000E 01 Qmin_alphas 0 16818E 00 pass 4 Table 10 Example of a Qcdnum error message 42 be declared double precision in the calling routine and actual values should be written in double precision format e g 0 188D0 Qcdnum floating point functions must also be declared double precision e Do not change Qcdnum parameters in between calculations such as evolving in NLO and calculating structure functions in LO e Input parton distributions must be the gluon the singlet or non singlet momentum den sities and not a mixture of these Qcdnum assumes that the gluon is stored in memory location ID 0 the singlet in ID 1 and the non singlet distributions in the remaining locations ID 2 10 e Qcdnum does not verify if starting values of parton distributions are suppli
14. structure functions however are given on the scale Q As an example let us assume that 3 Q 4 Then value QPDFXQ name x 8 DO iflag returns the parton distribution at u2 8 GeV whereas value QSTFXQ opt name x 8 DO iflag gives the structure function at Q 8 GeV This structure function is calculated from the parton distributions at u2 Q 4 2 GeV Indeed the reader may verify that in LO the Fy structure function given by QSTFXQ at 8 GeV is equal to the quark distribution given by QPDFXQ at 2 GeV 4 14 Fast structure function calculation The calculation of structure functions is expensive because they require in NLO the evaluation of convolution integrals This is particularly true for QSTFXQ x q2 where in general the result is obtained from linear interpolation of four structure functions calculated at the four grid points surrounding x and Q The use of QSTFXQ is inefficient if one or more of these grid points coincide with those of a previous call This overhead can be eliminated if Qcdnum knows beforehand for which values of x and Q it has to calculate the structure functions The structure function calculations might thus be speeded up considerably as follows 1 Pass a list of x Q points where structure functions are to be calculated This is done by calling the subroutine QFMARK e g do i 1 ndata x xdata i q2 qdata i call QFMARK x q2 38 This rou
15. the Q value is returned by q2 QFROMIQ iq The inverse function ix IXFROMX x returns the grid index ix for given x If x does not coincide with a grid point ix is negative and ix corresponds to the closest grid point which lies below x If x is outside the grid boundary ix is set to zero Thus ABS IXFROMX x can be thought of as a function which returns for given x the bin number i with the grid point x being the lower edge of the bin The corresponding function for the Q grid is iq IQFROMQ q2 The function ix IXNEARX x is similar to IXFROMX but here the grid index ix corre sponds to the grid point closest to x instead of that below x Again if x coincides with a grid point ix is positive otherwise negative ABS IXNEARX x can be thought of as a function which returns for given x the bin number i with x being the center of the bin The corresponding function for Q is iq IQNEARQ q2 The current grids including the logarithmic heavy quark x grid can be written to logical unit number lun by a call to QPRINT lun Xqgrid Finally when defining the z Q grid one should bear in mind that see section 3 The grid should be dense enough so that the linear quadratic interpolation in x In Q used by Qcdnum is a sufficiently good approximation The CPU time increases quadratically linearly with the number of grid points in x Q The starting scale of the evolutions Q2 must be included in the g
16. these distributions come from they may be the result of a Q evolution or for instance be read in directly from Pdflib and passed to memory by the appropriate calls to QNPSET see section 4 15 The current value of a parton density is returned by value QPDFIJ name ix iq iflag value QPDFXQ name x q2 iflag where name is a valid name of a memory resident distribution or of a linear combination e g proton see section 4 9 The function QPDFXQ gives by linear interpolation the parton distribution at any value of x and Q within the grid If ix iq or x q2 are outside the grid boundaries or cuts a zero value is returned and iflag 1 0 otherwise To protect against use of undefined parton distributions or structure functions Qcdnum will abend with an error message This behavior can be switched off by setting the logical variable LIMCK to false call QNLSET LIMCK false The structure function corresponding to a quark density or linear combination is calculated by value QSTFIJ opt name ix iq iflag value QSTFXQ opt name x q2 iflag where opt can be set to F2 FL or XF3 see below for the heavy quark structure func tions The calculation is in LO or in NLO depending on the value of order see section 4 10 It is the responsibility of the user to evaluate the structure function from the appropriate distri bution e g xF
17. to 20 different combinations of opt and name If there are more than 20 different combinations STFAST acts as a do nothing call STFCLR Clear the memory allocation made by STFAST Previous results calculated by STFAST are inval idated call QPRINT lun opt Print Qcdnum info on logical unit number lun A logical unit number other than 6 should be opened by the user The input parameter opt can be set to P Qcdnum parameters and options UY x Q2 grid B Booking list of parton distributions 8 Stats on structure function calls T Timelog 7A All of the above A printout of the timelog option T causes a halt of the logging If the log is to be continued please call QNTIME Start or QNTIME Continue after QNPRINT 52 7 Acknowledgments I am of course much indebted to the original authors of Qcdnum M Virchaux and A Ouraou I would like to thank M Virchaux for making the code available to me in 1991 and for his contributions during the first steps on the tortuous road to bring the program up to its present version J O Mara and M Vreeswijk contributed to the Qcdnum upgrades for the QCD analysis of the ZEUS 1993 F data The NLO heavy quark coefficient functions were taken from a program kindly provided by S Riemersma Qcdnum results were regularly checked against those obtained from a similar program developed independently by M Vreeswijk Detailed comparison with several ot
18. to be unequal 8 For clarity we will drop the argument Q in the following and write for the non singlet contri butions to Fk Qs FPE anle 2s f eo 2 o 19 2T z zZ where q is a shorthand notation for q7 q or any linear combination of these These structure functions do not depend on the gluon distribution The first term on the r h s of eq 18 and 19 is the leading order contribution to F and F7 respectively From eqs 18 and 19 the structure function for any linear combination q c E q can be written as Qs 1 dz x x F x Spo xq a P fog g z 4 c Ce s 20 The structure function zF is a pure non singlet Qs shla zale Sta f Zo 2 eto 21 2m zZ E The coefficient functions Go x a x and Ey x are calculated in the MS scheme in 7 and are given in appendix A As an example we consider the proton F gt gz structure functions in charged lepton scattering which can be calculated from eq 20 using the quark distributions given in eqs 12 and 13 in the previous section The neutron distribution is obtained by assuming isospin symmetry qu and qq are interchanged in eq 12 The deuteron quark density is usually defined as q q 2 so that q z d Auala 22 The distribution A a is thus constrained by the difference of proton and deuteron F structure functions To calculate Ff we need to evolve one more non singlet dist
19. 12 e ORDER 1 2 selects LO NLO calculations Its value can be set at any time Default NLO e SCAXO and SCAQO are related to the logarithmic linear scale of the x and Q grid see section 4 6 15See section 4 1 on how to change this 23 ic ecw ii typ deflt Ta el Nodd Huu 3g ATH S SJ gd dgdgdgdgdggdurrrr r rntebetttrt SU UJ SJ DU GU DJ DW T T Oo OO OO OO OO OO OO One AHA HAA A AHA TAHA 2 20000E 00 10000E 11 15000E 01 50000E 01 15000E 01 50000E 01 18800E 03 13000E 00 83152E 04 LOOOOE 01 00000E 00 10000E 01 00000E 00 10000E 01 00000E 00 10000E 11 10000E 11 10000E 01 10000E 01 LOOOOE 01 10000E 01 A7568E 00 Table 3 Default E descriptio Analytical Numerical Numerical Structure F2_charm F2_bottom FL_charm FL_bottom Check x Q Heavy F2 FL only for Q2 gt 1 5 GeV2 LO 1 or x grid sc Q2 grid sc C mass for mass for mass for mass for B C B T mass for n LO weight calculation LO weight calculation NLO weight calculation function NLO weights weight calculation weight c
20. 3 in section 2 1 Input to the calculation is the value of a at a given Q default 0 180 and 50 GeV respectively These values can be changed at any time by e g Mz 91 188 As 0 116 call QNRSET ALFQO Mz Mz call QNRSET ALFAS As As described in section 2 1 the flavor thresholds in the a evolution are taken to be equal to the quark masses which can be set by calls to QNRSET see section 4 5 Notice that you can also set the top mass default 188 GeV so that a can be evolved with up to f 6 flavors However there is no top threshold in the parton distribution evolutions or in the structure function calculations so that there Qcdnum uses only up to f 5 flavors see below The function alfas QALFAS Q2 qlambda nf ierr gives the value of a for any Q The function also returns the value of A qlambda and the number of active flavors nf If Q lt A the value of a is arbitrarily set to 100 A to zero and the flag ierr 1 0 otherwise 3 The flavor thresholds used in the evolution and structure function routines are not related to the quark masses Qcdnum uses instead the thresholds Q and Q such that f 3 for Q lt Q f A for Q lt Q lt Q and f 5 for Q lt Q At initialization Qcdnum sets Q Q at minus plus a very large value so that by default f 4 for all Q You can change this at any time by a call to call QTHRES Q2c Q2b Thus if you want to use thre
21. 4 lt pi lt 4Q for instance This variation does not affect the evolution but changes the coefficient functions and therefore relation between the structure functions given at the scale Q and the parton distributions given at the scale p 3 It is probably safe to calculate the structure functions down to Q 0 5 GeV or so 5 Changing the scales u or 4 does not require a re calculation of the weight tables 37 4 4 4 4 4 var typl deflt value description 4 4 4 4 4 AAAR2 R 1 0 0 10000E 01 R2 A M2 B ren scale BBBR2 R 0 0 0 00000E 00 R2 A M2 B ren scale AAM2L R 1 0 0 10000E 01 M2 A Q2 B light fact scale BBM2L R 0 0 0 00000E 00 M2 A Q2 B light fact scale AAM2H R 1 0 0 10000E 01 M2 A Q2 B heavy fact scale BBM2H R 0 0 0 00000E 00 M2 A Q2 B heavy fact scale 4 4 4 4 4 2 Table 6 Default settings for the renormalization and mass factorization scale Because in the previous sections no reference is made to the different scales they were all identified with Q some confusion may arise when 3 is varied with respect to Q The rule is that parton distribution input evolution and output are defined on the scale u34 The
22. BOOK 6 xqval Up to 20 linear combinations of quark distributions with identifiers id 11 to id 30 can be defined by call QNLINC id name nf factors where factors is a 10 dimensional array which contains the multiplication factor for the singlet distribution in factors 1 and that for the non singlet distributions in factors 2 through factors 10 Because such linear combinations often depend on the number of fla vors Qcdnum maintains for each identifier three factor arrays for f 3 4 or 5 For instance the proton quark distribution in charged lepton scattering for f 3 flavors can be written as a linear combination of the singlet and non singlet distributions 4 1 1 d ecce Ae qf 3 18 F 6 ud 6 ds Taking the non singlet distributions from the booking list given above the proton can be defined by the following code double precision factors 10 data factors 10 0 D0 Define proton distribution for 3 flavors factors 1 4 18 factors 2 1 6 factors 3 1 6 call QNLINC 11 proton 3 factors Likewise we have for four flavors 5 1 1 1 y MEE eNOS Ny a ENT dj jg t Gud gh gc factors 1 5 18 factors 2 1 6 factors 3 1 6 factors 4 1 6 call QNLINC 11 proton 4 factors Notice that only the factor arrays are stored in memory and not the linear combinations themselves Each time a linear combination is addressed Qcdnum evalua
23. RR ree SS RISE ESOS 3 8 Evolution in Q 4 Qcdnum user guide 4 1 The Qcdnum program 4 2 User program 4 3 Example program 4 4 Initialization 4 5 Qodnum variables and options cos curri we ates ees ee eee A RET RS 4 6 Definition of the grid 4 7 Weight calculation 4 8 Disk dump read of weight tables oa aeree Wane a Ke LA 4 9 Definition of parton distributions exei eu 4d Pe ew ee Se 4 10 QCD evolution parameters and options 2 2 2058 4 11 QCD evolution of parton distributions 204 4 12 Access to parton distributions and structure functions 4 13 The renormalization and factorization scale 0 0 0 0008 8 4 14 Fast structure function calculation n 4 15 Structure functions from external parton distribution sets 4 16 Access to more information 0 000 ss 4 17 Qcdnum error handling 5 Size accuracy and speed 11 13 13 14 15 16 16 19 21 23 23 25 27 28 29 32 34 36 37 38 40 40 42 43 6 Qcdnum subroutine calls reference section 6 1 Description of subroutines ox fee8 h Oh P atat eoe ep dr eoe deve 0B e ee 7 Acknowledgments A QCD splitting and coefficient functions 1 Introduction Qcdnum is a fast QCD evolution program which provides Calculation of the Q evolution of a up to NLO e Q evolution of the gluon singlet and non singlet distribu
24. The code uses the predefined names singlet gluon to access the singlet quark gluon distribution The quark and gluon distributions at Q are given by the user defined functions quarks x and gluons x The singlet F5 structure function is calculated from the evolved distributions patch QUCOMM user defined sequences if any patch QUPROG implicit double precision A H 0 Z x 1 Initialization call QNINIT call QNTIME start x 2 Setup grid in x and Q2 q0 7 call GRXDEF 80 4 D 4 call GRQDEF 59 4 D0 5000 D0 call GRQINP q0 1 call GRGIVE nx xmi xma ng qmi qma 21 x 3 Create weight tables call QNFILW 0 0 x 4 Set charm bottom threshold at Q2 4 25 GeV2 call QTHRES 4 D0 25 D0 5 Set input value for alphas Mz 91 2 call QNRSET alfas 0 118D0 call QNRSET alfQ0O Mz Mz 6 Input singlet and gluon at Q2 7 GeV2 quarks x and gluons x are user defined functions iq0 IQFROMQ q0 do ix i nx x XFROMIX ix call QNPSET singlet ix iqO quarks x call QNPSET gluon ix iq0 gluons x enddo 7 Singlet gluon evolution call EVOLSG iq0 1 nq x 8 Get singlet gluon F2 etc at any x and Q2 X 0 005 q2 20 qval QPDFXQ singlet x q2 iflag gval QPDFXQ gluon x q2 iflag fval QSTFXQ F2 singlet x q2 iflag x 9 Some printout call QPRINT 6 all end The user defined function gluons x for instance might look a
25. Zeus Note 97 066 QCDNUMI16 A fast QCD evolution program M A J Botje NIKHEF PO Box 41882 1009DB Amsterdam the Netherlands August 12 1998 QCDNUM version 16 10 12 Abstract Qcdnum is a program which numerically evolves parton distributions using the Altarelli Parisi QCD evolution equations in LO or in NLO in the MS scheme From the evolved distributions the structure functions P5 Fr and xF3 can be calculated The program can handle flavor thresholds at fixed values of Q light quark variable flavor number scheme or alternatively evolve the three light quarks only and calculate the charm and bottom contributions to the structure functions from the photon gluon fusion process heavy quark fixed flavor number scheme To study the scale uncertainties it is possible to independently vary the renormalization scale and the mass factorization scale In this report we describe in detail how to use Qcdnum in a QCD analysis of unpolarized structure functions Contents 1 Introduction 2 The QCD analysis of structure functions 2 1 QCD evolution of a 2 2 QCD evolution of parton densities x ace woe uS es ur Ee we BES 2 3 Structure function evaluation 0 00 bee eee 2 4 Renormalization and factorization scale dependence 2 5 Heavy flavor contributions oen ta cR SNS aa kd alate are Ie dA PES ene du 3 Numerical method 3 1 Overview 3 2 Evaluation of convolution integrals vou yore
26. alculation weight calculation weight calculation 2 limits and cuts NLO 2 calculations ale from log gt linear ale from log gt linear F2c FLc GeV F2b FLb GeV alpha s evolution GeV alpha s evolution GeV alpha s evolution GeV Value of alpha s Q2 R2 R2 M2 M2 M2 M2 A M2 A M2 A Q2 A Q2 A Q2 A Q2 Charm thre Bottom thr Xmin cut Qmin cut Qmax cut Roots cut Lowest Q2 where alpha_s is given GeV2 ren scale ren scale light fact light fact heavy fact heavy fact scale scale scale scale shold eshold 1e 0 1e 0 1e 0 no cut 1e 0 no cut with valid alpha s cut cut no no settings in Qcdnum 24 Al i M M fa pi al lc i A ey Si Na gi ig in Sea i i ec id a ase Sp a Sa ee sh ey a M H e MCSTF and MBSTF are the charm and bottom mass in GeV used in the heavy quark structure function calculations see section 4 12 e MCALF MBALF and MTALF are the charm bottom and top mass GeV which determine the flavor thresholds in the a evolution e ALFAS and ALFQO are the value of a and the Q where it is given These values may be changed at any time e AAAR2 and BBBR2 define the renormalization scale u aj b see section 4 13 e The mass factorization scale 4 aQ b can be defined for the light heavy
27. call QNFILW 0 O 48 Generate weight tables needed for the calculation of convolution integrals in the QCD evolution and structure function routines Should be called after the z Q grid is defined For a complete description of QNFILW see section 4 7 call QNDUMP lun Write weight tables to logical unit number lun which should have been opened before by the user e g OPEN unit 24 file weights file form unformatted status unknown Written to the file are i the Qcdnum version number ii the parameters mxx and mq2 see section 4 1 iii the z and Q grid iv the current value of the charm and bottom mass and v those weight tables which are available when QNDUMP is called See section 4 8 for details call QNREAD lun istop ierr Read weight tables from logical unit number lun which should have been opened before by the user see above This routine will not only read the weight tables but also overwrite the z Q grid if any and the charm bottom mass in memory See section 4 8 for a description of the input parameter istop and the output flag ierr call QNBOOK id name Associate memory location id with the name of a parton density e id 0 reserved for the gluon distribution with the predefined name gluon e id 1 reserved for the singlet quark distribution with the predefined name singl e id 2 10 reserved for non singlet dist
28. e back use OPEN unit 24 file weights file form unformatted status unknown CALL QNREAD 24 istop ierr However this will not only read the weight tables but also overwrite the current definition of the z grid if any and in case the heavy quark tables are stored that of the Q grid together with the value of the charm and or bottom mass To avoid unwanted side effects QNREAD compares the grid definition if any and the quark masses in memory with those on the file and sets the error flag ierr as follows 28 ierr 0 all OK ierr 1 xgrid exist in memory ne that on the file ierr 2 file contains heavy quark weight tables and qgrid exist in memory ne that on the file ierr 3 file contains charm weight tables and c mass in memory ne that on the file ierr 4 file contains bottom weight tables and b mass in memory ne that on the file The action taken by QNREAD depends on the input variable istop istop 0 read the file whatever the value of ierr istop 1 read only when ierr 0 do nothing otherwise istop 2 stop the program when ierr ne O Thus setting istop 1 provides the possibility to inspect the value of ierr take the appro priate action if any and then call QNREAD again 4 9 Definition of parton distributions Qcdnum can hold up to 11 parton momentum densities in memory which internally are ad dressed by an identifier running from id 0 to id 10 The user addresses these mem
29. e flavors throughout then set both Q and Q beyond the highest Q grid point to use 5 flavors set them both below the lowest Q grid point As an example we give below the code to set the thresholds to the square of the quark masses as in the a evolution 2Thus qlambda and nf are output variables and not input variables as the syntax may suggest 3 However always keep Q2 lt Q in the call to QTHRES otherwise a fatal error will result 32 call QNRGET MCALF Cmass call QNRGET MBALF Bmass call QTHRES Cmass Cmass Bmass Bmass The function nf NFLGET iq returns the number of flavors as used by Qcdnum for any grid point iq 4 As explained in section 4 8 weight tables can be calculated once and for all and dumped to disk for the next passes of your analysis It is convenient to do this on a large grid say 120 x grid points covering 1079 lt x lt 1 and 60 Q points covering 0 1 lt Q lt 10 GeV However one might not always be interested in such a large kinematic range Instead of creating a new weight file on a smaller grid you can set cuts so that parton distributions are evolved over part of the z Q grid only This may result in a considerable reduction in CPU time The cuts are passed to Qcdnum by call GRCUTS xmin qmin qmax roots where the meaning of xmin qmin and qmax should be obvious roots y s restricts the z Q domain to within the kinematic limit Q lt zs Once
30. ed it checks if they are booked though e Evolution routines do not clear the memory before the evolution Thus if you evolve over a restricted range iqmin iqmax in Q or apply cuts then outside this range the value of the parton distribution will remain to be whatever it was before the evolution e Non singlet xq distributions are discontinuous at the flavor thresholds if any The discontinuities are not automatically taken into account by the evolution routines e xF is a non singlet structure function and should therefore be calculated from a linear combination of non singlet densities only e Heavy quark F gt z structure functions should only be calculated from linear combinations representing the proton deuteron or neutron quark distributions in charged lepton scat tering The parton distributions must have been evolved with f 3 flavors for all Q 5 Size accuracy and speed The size of the Qcdnum executable can be controlled by setting the parameters mxx and mq2 to the maximum number of grid points you want to use in your analysis see section 4 1 For mxx 100 150 and mq2 40 the size is 3 3 5 2 Mbyte as measured with the Unix command size The accuracy of Qcdnum depends on the choice of the z Q grid as illustrated in fig 2 Here we take as a reference parton distributions evolved from Q 4 up to Q 10 GeV on a 400 x 60 x Q grid covering a large range 10 lt x lt 1 For such a grid Qcd
31. en value of x e If x is outside the grid boundaries ix is set to zero 47 e If x coincides with a grid point ix is positive and corresponds to that grid point e If x does not coincide with a grid point ix is negative and ix corresponds to the closest grid point which lies below x iq IQFROMQ q2 As above but now for the Q grid ix IXNEARX x As IXFROMX but ix now corresponds to the grid point closest to x instead of that below x iq IQNEARQ q2 As above but now for the Q grid call GRCUTS xmi qmi qma roots Pass cuts to Qcdnum such that evolutions are only performed in the kinematic domain defined by zi x lt 1 Q2 lt Q lt Qar and Q lt zs where roots 4 s is the center of mass energy To release a cut or not apply it at all set the corresponding parameter to a value lt 0 GRCUTS can be called at any time even before the z Q grid is defined istat IFAILIJ ix iq Returns a non zero value if ix or iq falls outside the grid boundaries or cuts 0 1 no yes fail QMINA cut 0 1 no yes fail ROOTS cut 0 1 no yes fail QMAX cut 0 1 no yes fail QMIN cut 0 1 no yes fail XMIN cut ifail ijklm Bor WU e The QMINA cut is set automatically by Qcdnum to the lowest Q grid point above the current value of the QCD scale parameter A istat IFAILXQ x q2 As above but now for z Q
32. entum densities H x Q zq x Q Extension of the method to singlet gluon evolution is straight forward For convenience we introduce the variable t In Q and write for a given grid point x t the evolution equation for H as cf eq 35 Baty 02 Ys Plena S Pana Hent 6 1 with H 0H dt This slope cannot be calculated directly because the first term in the sum on the r h s of eq 36 contains H x tj which is a priori unknown Separating this term from the sum and denoting the remaining sum which runs from k i 1 to n by S we write eq 36 in shorthand notation as H wH TS 37 Next we assume that the following condition is fulfilled Condition i In previous steps of the evolution H x t has been calculated for all values of x gt zj ie for k gt 1 With this assumption the sum S can be evaluated and as a consequence eq 37 becomes a linear equation with two unknowns H j and H A second linear equation is obtained as follows the density H z t is quadratically interpolated between the grid points _ and t H z tji lt t lt t at bt c With this interpolation H z t is related to H xj t j 1 H x tj and H t 1 through 1 Ht Hit 2 H zi t5 H zi tj 1 At 38 with At t t 1 If we assume that the following condition is satisfied 13The weight tables for the splitting and light quark coefficient function
33. genlib lpacklib lkernlib qcdprog gt qcdprog log lt lt EOF User datacards e g MINUIT if any EOF Table 1 Example of a unix script running Qcdnum 18 REPL QCDCOM QCDCOM 6 7 PARAMETER MXX 150 PARAMETER MQ2 60 See section 5 for further details on program size accuracy and speed A unix script which performs a patchy run compiles the resulting fortran code makes the executable and runs it might look like that shown in table 1 In this example the Qcdnum code is stored in the file qcdnum16 car and user code in qcduser car The user code contains two patches qucomm and quprog An example of such code is given in section 4 3 4 2 User program An analysis program based on Qcdnum should contain the following steps 1 2 Initialization Definition of the z Q grid Calculation of weights The weight calculations depend on the grid definition which therefore cannot be changed afterwards By default weight tables are computed for the LO and NLO splitting functions as well as for the NLO F gt Fr and zF coefficient functions Tables which depend on the number of flavors are generated for f 3 4 and 5 Optionally weight tables are calculated for the heavy quark coefficient functions Definition booking of the parton distributions Internally Qcdnum reserves space for the gluon and the quark singlet momentum density and for up to nine user defined non singlet densities Furthermore up
34. hecks that you do not evolve over the thresholds without adding the discontinuities i e call EVPLUS splus ig0 1 nq2 produces in this example a fatal error message The routine call QADDSI name iq factor adds for all x at the grid point iq factor times the singlet distribution to a non singlet distribution Likewise the charm distribution can be evolved from zq x Q2 1 AxX x Q2 by Notice that we first set aq to zero QNPNUL for all x and Q before supplying the starting value with QADDSI Otherwise in QCD fits for instance Qcdnum would continue to add 1 4X to the starting value of zq at each Minuit iteration 35 call QNPNUL cplus factor 1 4 call QADDSI cplus iqc factor call EVPLUS cplus iqc igc iqb factor 1 4 1 5 call QADDSI cplus igb factor call EVPLUS cplus iqb iqb nq2 Finally we remark that although Qcdnum knows that it deals with a non singlet distribution it has no information on what kind of distribution it actually is xq xA or xq It is therefore the responsibility of the user to call the proper evolution routine and in case of zq evolution to handle the threshold discontinuities correctly 4 12 Access to parton distributions and structure functions In this section we describe how to access parton distributions and how to calculate structure functions from these Qcdnum does of course not care where
35. her NLO evolution codes was made during the HERA 1995 1996 workshop here I much enjoyed the collaboration with J Bl mlein C Pascaud S Riemersma A Vogt and F Zomer I am grateful to A Vogt for many discussions and for further checks beyond those made during the HERA workshop Finally I would like to thank my colleagues from NMC and ZEUS for help and encouragement A QCD splitting and coefficient functions The leading order splitting functions are 14 P 2 lu Saq S 2 PO f a zy 42 Poe 5 0 2 pO 6 lacu c z 1 2 4 5 4 i12 where the so called prescription is defined by ols Fle 0 2 f fona 43 Notice that fols f FO Eed 10 alea x and in particular that ae dz f 2 E D scstiltstt ess The coefficient functions are given by 7 cQ 1 hs m 2 19 o 3 1 z c 4 l r S at ay 8 Ma Il e anh 53 1 4 8z 1 2 44 8 Coa gf CMa Afz 1 x 4 C CN t 45 Notice the prescription in pr a non singlet P5 structure function for instance is thus given by F 0 aq a 5 t 3 ya y a EN RE od du w 46 2T References 11 12 13 14 A Ouraou Ph D Thesis Universit de Paris XI 1988 M Virchaux Ph D Thesis Universit de Paris VII 1988 M Virchaux and A Milsztajn Phys Lett B274 1992 221 NMC M Arneodo et al Ph
36. hich is usually set to u Q or wi Q 4m The kinematic domain where the heavy quarks contribute is restricted by the requirement that the square of the y p center of mass energy must be sufficient to produce the hh pair W M Q 1 2 x gt M 4m i e the lower integration limit az lt 1 in eq 29 For the heavy quark coefficient functions C and D in eq 28 we refer to 13 Qcdnum however will continue to evolve q as a non singlet q distribution below the charm threshold if you ask the program to do so l0Notice that the proton distribution is continuous at the flavor thresholds although it is build up from distributions some of which are discontinuous In fact a powerful check on the consistency of a Qcdnum based analysis is provided by verifying that no discontinuities show up in distributions which should be continuous 11In the LO and the first two NLO terms the virtual photon couples to the heavy quark hence the factor e in eq 28 The last NLO term describes the process where the virtual photon couples to a light quark which subsequently branches into a hh pair via an intermediate gluon hence the appearance of the charge weighted sum q of light quark distributions I Some of these coefficient functions are given as interpolation tables since they are too complex to be cast into analytical form All heavy quark coefficient functions in Qcdnum are taken from code provided by S Riemersma 12 A QCD ana
37. i xma nq qmi qma GRXOUT xgrid GRQOUT qgrid xx XFROMIX ix q2 QFROMIQ iq ix IXFROMX xx iq IQFROMQ q2 ix IXNEARX xx iq IQNEARQ q2 GRCUTS xmi qmi qma roots istat IFAILIJ ix iq istat IFAILXQ x q2 QNFILWC 0 O QNDUMP lun QNREAD lun istop ierr QNBOOK id name QNLINC id name nf factors id IPDFID name QTHRES t34 t45 nf NFLGET iq as QALFAS q2 qlam nf ierr QNPSET name ix iq value QADDSI name iq factor QNPNUL name EVOLSG iq0 iqmin iqmax EVOLNP name iq0 iqmin iqmax EVPLUS name iqQ iqmin iqmax EVOLNM name iq0 iqmin iqmax val QPDFIJ name ix iq ifl val QPDFXQ name x q2 ifl NAI NA wes we val QFMARK x q2 QFMNUL STFAST opt name STFCLR QPRINT lun opt USTELI C opt name ix iq ifl val QSTFXQ opt name x q2 ifl Clear the z Q grid Give maximum number of grid points Get definition of the current grid Copy z Q grid to local array Get z Q for given grid index Get grid index for given x Q Get grid index for given x Q Set cuts ix iq fails passes cuts x Q fails passes cuts Calculate weight tables Write weight tables to disk Read weight tables from disk Book non singlet distribution Define linear combination Get id of quark distribution Set flavor
38. ication program which in addition to Qcdnum might use Pdflib for parton distribution input or Minuit for QCD fits 4 1 The Qcdnum program The Qcdnum code is available as a patchy car file The patches to be selected in the patchy run are 16 USE DOUBLE Select double precision version TUSE QCDCOM T EXE Qcdnum common blocks USE QCDNUM T EXE Qcdnum code light quarks USE QCUTIL T EXE Utility routines USE QHEAVY T EXE Qcdnum code heavy quarks The switch use double selects the double precision version of Qcdnum To obtain the single precision version you simply omit the line use double or comment it out c use double It is recommended to always compile the program in double precision except on machines with 60 bit arithmetic like the CRAY When the double precision version is used e floating point variables in Qcdnum subroutine calls and Qcdnum floating point functions should be declared double precision in the calling routine e g double precision val x q2 QPDFXQ val QPDFXQ name x q2 iflag The easiest is of course to simply declare all floating point variables as double precision in the calling routine implicit double precision A H 0 Z e Actual values passed to subroutines must be written in double precision format e g 1 3D0 instead of 1 3 in for instance call QNRSET MCSTF 1 3D0 In the remainder of this writeup we assume that Qcdnum is run in double precision mode
39. iqmin iqmax This feature is of course irrelevant for the singlet gluon evolution but useful in the context of xq evolution with discontinuities at the flavor thresholds see below A non singlet distribution of the type xq valence is evolved by 34 call EVOLNM name iq0 iqmin iqmax and non singlet density of the type xA by call EVOLNP name iq0 iqmin iqmax where name is a valid name of a non singlet distribution that is not gluon or singlet or whatever names are associated with id O and 1 A non singlet zq distribution can also be evolved by a call to EVOLNP but here one has to be careful with discontinuities at the flavor thresholds see section 2 2 and 2 5 As an example we assume that 1 lt iqc lt iq0 lt iqb lt nq2 where iqc and iqb are the grid indices corresponding to the charm and bottom thresholds respectively To evolve for instance the vst distribution you might code Evolve from iq0 down to iqc and up to iqb call EVPLUS splus iq0 iqc iqb Downward splus jumps at iqc because f 4 gt 3 factor 1 4 1 3 call QADDSI splus iqc factor call EVPLUS splus iqc 1 iqc Upward splus jumps at iqb because f 4 gt 5 factor 1 4 1 5 call QADDSI splus iqb factor call EVPLUS splus iqb iqb nq2 Here we have introduced the routine EVPLUS which is identical to EVOLNP except that Qcdnum c
40. lag LIMCK has been set to false val QPDFXQ name x q2 iflag As above but now for x and Q val QSTFIJ opt name ix iq iflag As QPDFIJ but calculate the structure function opt of parton distribution or linear combina tion name The input variable opt can be set to F2 FL xF3 F2C FLC F2B or FLB Qcdnum accepts every combination of opt and name except the name associated with the gluon distri bution The output variable iflag is set to e 1 ix and or iq outside grid or cuts In this case val is set to zero Qcdnum abends with a fatal error unless LIMCK is set to false 51 e 0 calculation successful e 1 structure function obtained from the results of STFAST see below For further details on the use of QSTFIJ see section 4 12 4 14 and 4 15 val QSTFXQ opt name x q2 iflag As above but now for z and Q call QFMARK x q2 Mark the four grid points surrounding z and Q for fast structure function calculation Qcdnum abends with a fatal error message if x and or Q are outside the grid boundaries call QFMNUL Clear all the marks previously set by QFMARK call STFAST opt name Calculate the structure function opt see above of parton distribution name for all grid points marked by QFMARK The results are stored internally for later interpolation by QSTFIJ or QSTFXQ for up
41. lysis of charged lepton structure functions now might proceed as follows the parton distributions given at some scale Q Q2 are evolved in Q with f 3 flavors The light quark contribution u d s to the proton quark density is assembled from the evolved distributions and the structure functions F gt z are calculated from eq 20 with f 3 flavors To these are added the heavy flavor contributions F and F as calculated from eq 28 above The resulting structure functions are then compared to data Notice that such heavy flavor corrections are process dependent the method described above applies to charged lepton deep inelastic scattering only Finally we remark that Qcdnum will calculate F from any quark distribution e g singlet non singlet but that the result only makes sense when obtained for the proton deuteron or neutron distributions as defined in section 2 3 for f 3 flavors 3 Numerical method 3 1 Overview The Qcdnum program calculates the Q evolution of parton momentum densities in LO or in NLO in the MS scheme on a user defined grid in z and Q The value of these parton densities at fixed Q Q has to be specified for each grid point in x This implies that Qcdnum does not use nor cares about parameterizations describing the x dependence of the input distributions The calculation of the logarithmic slopes in Q is based on the computation of convolution integrals see eqs 6 and 8 These are calc
42. may serve as a check on the accuracy of the numerical integration Inserting eq 32 into eqs 6 and 8 the evolution equations for the parton momentum densities can be written as weighted sums e g for the the singlet quark momentum density we have dx xo Qs us dlnQ E 24 2 wh zi 0 S Up rss zo zX x i 0 As wie riso gt wore o igle 35 Likewise the expressions for the structure functions F gt Fr and x F3 in NLO see eqs 18 21 can be written as weighted sums The weight tables for the heavy quark coefficient functions depend on Q which makes them 3 dimensional Their size is however much reduced by calculating the heavy quark structure 14 functions on an equidistant logarithmic grid in x this grid is automatically maintained by Qcdnum For such a grid the weights depend only on the difference x xo thus saving a factor n 2 in storage In this way fast calculation of the logarithmic slopes and of structure functions is achieved entirely based on lookup tables which can be pre calculated with high and controlled accuracy The z grid must of course be dense enough so that linear interpolation in x eq 31 is a good approximation Since parton densities are smoothly varying functions of x this can be achieved with relatively few grid points typically n 100 3 3 Evolution in Q In this section we describe the calculation of the Q evolution of non singlet quark mom
43. n the evolution limits First we have to assemble the proton quark distribution using eqs 15 16 and 17 for the regions Q lt Q Q Q Q and Q lt Q respectively F is then calculated from eq 20 Qcdnum provides a mechanism to define objects like the proton as linear combinations of quark distributions for f 3 4 5 flavors The outline given above should be considered as an example to illustrate the handling of the flavor thresholds in practice it would be easier and more efficient to use the proton non singlet distribution defined in eq 13 in section 2 2 instead of the Aya q gf and qj used above An alternative treatment of heavy flavor threshold effects in charged lepton scattering is based on the calculation of the heavy flavor contributions to F5 or Fr through in NLO 13 Os Qs Fe z Q fagot on ge Ci re ESC a DE j 28 where e is the charge of the heavy quark in units of the electron charge and q denotes the proton quark distribution for f 3 light flavors see for instance eq 15 in section 2 2 The first term in eq 28 is the LO contribution from the photon gluon fusion process 7 g gt hh The last three terms correspond to the NLO subprocess 4 g hhg and y q hhq In eq 28 we have introduced the shorthand notation S S a 1 Z 2 fec e 9030 fE Ee e Cale nm 29 where a 1 4m Q and u is the factorization equals renormalization scale w
44. nce the fit has converged iflag 3 in Minuit language 44 CPUsec 600 _ Standard L a Total time 400 L SCU 200 B ee pieces 202 Stfast 0 400 b Evolution _ Standard 300 200 E KSCU 100 0 c 150 c Structure functions Standard 100 50 eig F Stfast o La pay 4 tio tio iio or 40 60 80 100 120 140 160 N ly Figure 3 CPU time spent for 2500 evolutions and 250000 F calculations as function of the number of z grid points for a 60 point Q grid a Total time b time spent in QCD evolutions c time spent in structure function calculations The dashed curves correspond to evolutions on part of the grid and the dotted curves to invoking fast structure function calculation in addition 6 Qcdnum subroutine calls reference section In this section we give a short description of all Qcdnum subroutine calls see also table 2 in section 4 2 We use the convention described in section 4 2 to denote logical integer floating point and character variables Output variables are indicated with an asterisk 6 1 Description of subroutines call QNINIT To be called before anything else Initializes a large set of constants and sets reasonable defaults call QNVERS version Ldouble nxmax nqmax Output e version current version number Should be declared character 8 in the calling program 45 e Ldouble returned as tr
45. num results were shown to be in agreement with several other NLO evolution codes to within 0 0596 5 Choosing n 200 100 70 the accuracy on F and on the quark distributions is within 0 1 0 5 1 at the highest Q as shown in fig 2a This is slightly worse for the gluon distribution fig 2b in particular at high x but here the gluon density vanishes The deviations are of course zero at the input scale 4 GeV and increase linearly with In Q to the values at Q 104 GeV shown in fig 2 33The accuracy deteriorates rapidly if the z grid has less than 70 points 43 j Z Olu CS on Se EN M gt OQ UI Q E a Q 10 GeV E b Q 10 GeV 4 L 4 L Be ees N 70 E 3 N 100 3 E N 200 C 2r 2 1 F Tee ae C 7 l 4 n spe Rm EU d Cul LAJIIN LL uul Ettir il J titi L LLLLU im L uuu saul L uu LL ul PEC 4 1 4 3 1 cum xg 49 49 xX iO 40 qd 497 X8 X Figure 2 a The deviation AF gt F and b Axg zg at Q 10 GeV as a function of x for a 70 100 and 200 point x grid and a 60 point Q grid The deviations are taken with respect to parton distributions evolved on a 400 x 60 z Q grid It turns out that Qcdnum is rather insensitive to the number of Q grid points the change in xg is not more than 0 0696 if n is reduced from 60 to 20 Figure 3 shows the CPU time spent for 500 passes through a Qcdnum based calculation as function of the number of z grid
46. oints to be added to the grid 16A logical unit number other than 6 should be opened by the user before a call to QPRINT The default zo 0 20 can be changed by call QNRSET SCAXO xO There is a similar parameter for the Q grid SCAQO which at initialization is set to a very high value 25 call GRXINP xarray nx call GRQINP qarray nq where xarray qarray are arrays containing nx nq points in x Q which may be given in arbitrary order One may or may not specify the point x 1 Qcdnum will always generate it As mentioned in section 3 2 the heavy quark structure functions are calculated on an equidistant logarithmic grid in x This grid is automatically maintained by Qcdnum with the same range xmin 1 and the same number of grid points nx as the grid defined by the calls to GRXDEF and GRXINP A call to GRXNUL GRQNUL sets the x Q grid to zero There are several ways to access the grid GRMXMQ nxmax nqmax returns the maximum number of grid points the user can define GRGIVE nx xmi xma nq qmi qma returns the number of grid points and the limits of the current grid GRXOUT xgrid and GRQOUT qgrid return the list of current grid points in the arrays xgrid qgrid which should be dimensioned in the calling routine to at least nx nq The value of x corresponding to a given grid index is returned by the function x XFROMIX ix where x is set to zero if ix is outside the grid boundaries Likewise
47. on and for the non singlet evolu tions of xq valence and zA However the non singlet evolution of a xq distribution is more complicated because rq is per definition discontinuous at the flavor thresholds The value of a parton distribution can be set for a given grid point ix iq by call QNPSET name ix iq value where name is a valid name of a memory resident distribution To set the distribution to zero for all z and Q call QNPNUL name The following call evolves the gluon and the singlet quark density in Q call EVOLSG iq0 iqmin iqmax where iq0 is the starting point and iqmin iqmax the range of the evolution You must take care that iqmin lt iq0 lt iqmax failing to do so will result in a fatal error and also that the gluon and singlet quark distributions are set to their proper starting values at iq0 by previous calls to QNPSET The calculation is in LO or in NLO depending on the value of order see section 4 10 For an example of user code evolving the gluon and singlet distributions we refer to the program listing given in section 4 3 We mention at this point that Qcdnum evolution routines do not clear the memory before the evolution This allows to chain the evolutions i e if iqmin lt iq0 lt iq1 lt iq2 lt iqmax then the sequence call EVOLSG iq0 iqmin igi call EVOLSG iqi iqi iq2 call EVOLSG iq2 iq2 iqmax is equivalent to call EVOLSG iq0
48. order The routine will always generate the grid point x 1 even if it is not included in the list xarray call GRQINP qarray nq As GRXINP but now add nq points to the Q grid 46 call GRtNUL e GRXNUL clear the x grid e GRQNUL clear the Q grid call GRMXMQ nxmax nqmax Returns the maximum number of grid points the user can define See section 4 1 on how to set these in the Qcdnum source code call GRGIVE nx xmi xma nq qmi qma Output e nx xmi xma number of grid points and limits of the current z grid The grid point x 1 is not included in the count nx GRGIVE should always return xma 1 e nq qmi qma number of grid points and limits of the current Q grid call GRXOUT xgrid Copy the current z grid to a local array xgrid which must be dimensioned to at least nx in the calling routine The array xgrid does not contain the point x 1 call GRQOUT qgrid Copy the current Q grid to a local array qgrid which must be dimensioned to at least nq in the calling routine x XFROMIX ix Returns the value corresponding to the grid point ix If ix is outside the grid boundaries or if the x grid is not defined XFROMIX returns zero QFROMIQ iq q2 As above but now for the Q grid IXFROMX x ix Returns the grid index ix corresponding to a giv
49. ory locations by a name The identifier id 0 is reserved for the gluon and id 1 for the quark singlet distribution with the predefined names gluon and singlet respectively The re maining 9 memory locations which should contain non singlet quark distributions can be associated to a user defined name by call QNBOOK id name with 2 lt id lt 10 The following rules apply e name can be given in upper lower or mixed case it is translated to upper case inter nally e name can be of any length but the first five characters must be unique Qcdnum trun cates name to five characters or adds trailing blanks if it is shorter To avoid confusion please refrain from defining names with leading blanks e An identifier cannot be booked twice unless it is released previously by setting the name to free Thus if one does not like the predefined name gluon for id 0 it can be renamed by call QNBOOK 0 free call QNBOOK 0 xglue 20They can only be evolved as non singlets and can only enter as non singlet contributions in the calculation of structure functions 29 As an example we give below the calling sequence to book the non singlet distributions zA 4 rqi tq xq and xq ru xd these distributions will be used in examples later on call QNBOOK 2 umind call QNBOOK 3 splus call QNBOOK 4 cplus call QNBOOK 5 bplus call QN
50. quark structure functions by setting the variables AAM2L BBM2L AAM2H BBM2H see section 4 13 Logical integer and floating point variables can be set by calls to QNLSET QNISET and QNRSET respectively e g call QNLSET WiNUM true call QNISET ORDER iorder call QNRSET ALFAS 0 187D0 The actual value of a Qcdnum parameter option can be retrieved by similar calls to QNLGET QNIGET and QNRGET The list as given in table 3 can be written to logical unit number lun by a call to QPRINT lun Parameters Also given in table 3 are the thresholds Q TCHRM Q TBOTT and four cut parameters XMINC ROOTS These cannot be set by QNRSET but by special routines described in section 4 10 The cut QMINA is set automatically by Qcdnum to delimit the region Q 7 A You can retrieve the values of the cuts with calls to QNRGET 4 6 Definition of the grid There are two ways in which the z Q grid can be specified e Specify the kinematic range and the number of grid points to be generated call GRXDEF nx xmin call GRQDEF nq qmin qmax GRXDEF clears the z grid if any and generates nx grid points in the range min x lt 1 Qcdnum sets the nx 1 th grid point to one The scale is logarithmic for x lt 0 20 and linear for x gt 0 20 1 The routine GRQDEF clears the Q grid and generates a logarithmic grid of nq points in the range Q2 lt Q lt Q min max e Pass a list of n grid p
51. r q x valence distribution 7 In this section we set the renormalization scale u3 see section 2 1 equal to y y Here and in the following we will suppress the argument u and write a for as u q x for q x u etc 6 g a ale 72l The evolution of these distributions does not depend on the gluon iS x Th a2 8 This equation is linear in the quark density so that any linear combination of the qf q is again a non singlet distribution of the type q q which obeys eq 8 in the following we denote by A the non singlet combination Ay z q x qj 2 lax ai 2 la 2 G 9 Notice that the sum of the q distributions vanishes 57 q7 0 The splitting functions can be expanded in a perturbative series in o which presently is cal culated up to next to leading order NLO 8 9 Pis PRla ZPP a 10 Pap x PiR x 5 PAR 2 where AB stands for FF FG GF or GG The expressions for the LO splitting functions P are given in appendix A Those for the NLO splitting functions P are complicated results calculated the MS renormalization scheme are given in refs 8 non singlet case and 9 singlet case Like for the beta functions introduced in section 2 1 the number of active flavors f enters as a parameter in several LO and NLO splitting functions The parton distributions h x u are continuous functions of jj but the
52. r instance wrong to store the charged weighted sum of quark distributions and calculate structure functions from this e When NLO structure functions are calculated from the singlet quark distribution or from a linear combination containing the singlet the gluon distribution must be supplied as well For non singlet structure functions the gluon distribution is not needed e When heavy flavor structure functions are calculated the parton distribution set must have been evolved with f 3 flavors only e When the parton densities are read in for part of the z Q grid only one may delimit this region by appropriate cuts so that Qcdnum can verify that no undefined distributions are used e Once the parton distributions are read in one can study the factorization scale dependence of the structure functions by varying 3 A variation of the renormalization scale does not make sense because this requires an evolution of the parton distributions see section 2 4 4 16 Access to more information Qcdnum info is printed on logical unit number lun by a call to call QPRINT lun opt opt P Qcdnum parameters and options X x Q2 grid B Booking list of parton distributions S Stats on structure function calls T Timelog 7A All of the above A logical unit number other than 6 should be opened by the user An example of a timelog is given in table 8 It gives for the evolution routines the number of calls the number of evolution
53. ribution Aya in addition to q7 defined in eq 13 The structure functions in neutrino scattering are usually given as the sums and differences of vN and PN results on isoscalar targets Here the relevant quark distribution is zX for FIN and zq z u x d d for zFY To calculate the latter we need to evolve a third non singlet parton distribution xq 2 4 Renormalization and factorization scale dependence In the previous section we have assumed that the renormalization scale u and the mass factorization scale 44 are equal to Q Here we discuss two cases where these scales are chosen to be different i u u Q and ii u A u Q We neglect here electroweak contributions from y Z interference and Z exchange which become important above Q 5000 GeV For a full description of charged and neutral current cross sections and structure functions in deep inelastic scattering see 10 Let us first vary the mass factorization scale with respect to Q for instance Q 4 lt pi Um lt 4Q This variation affects the relation between the parton distributions and the structure functions through a modification of the NLO coefficient functions The 14 dependence of a structure function F calculated at a given Q is defined as FO aya qm ts eA S NOS o 23 where F My X u r is the structure function falce up to order n at the scale i2 and AF TEN is the change in F when moving from 42 to Q
54. ributions as defined by the user At initialization the names of these locations are set to free e Qcdnum translates all names to upper case and truncates them to the first five characters which must be unique e An identifier cannot be booked twice unless the name has been previously set to free call QNLINC id name nf factors Define a linear combination of quark distributions e id 11 30 input identifier e name name of the linear combination e nf 3 5 the number of flavors For each combination of id and name you can define three different linear combinations for f 3 4 or 5 flavors 49 e factors 10 dimensional floating point input array containing the multiplication factor for the singlet distribution in factors 1 and that for the non singlet distributions in factors 2 through 10 id IPDFID name Returns for a given name the identifier id The input name should correspond to that of a memory resident quark distribution i e not gluon or the name of a linear combination Qcdnum abends with a fatal error message if it cannot find an identifier associated with the input name call QTHRES q2c q2b Set the flavor thresholds Q and Q Q lt Q Can be called at any time even before the z Q grid is defined At initialization Q Q is set to minus plus a very large value so that by default f 4 in all Qcdnum calculations Actual
55. rid see section 4 11 26 e It is convenient but not strictly necessary to also put the flavor thresholds Q and Q into the grid A good initial choice is about 100 grid points in x and about 40 in Q It may be worthwhile to spent some effort in optimizing the grid for QCD fit applications where the program is likely to be run many times on the same dataset A modest decrease in the number of x grid points by say 30 will increase the program speed by a factor of two For further details see section 5 4 7 Weight calculation Weight tables are calculated by a call to call QNFILW 0 0 By default all weights corresponding to the LO and NLO splitting functions and those for the NLO light quark coefficient functions are evaluated for f 3 4 and 5 flavors The weight calculations depend on the x grid which thus must have been defined before the call to QNFILW and cannot be changed afterwards The logical variables listed in table 4 define which tables are calculated in QNFILW These variables can be set by the appropriate calls to QNLSET e g call QNLSET Wi1NUM TRUE forces Qcdnum to calculate the LO weight tables numerically instead of analytically default If one is not interested in NLO calculations the generation of NLO weight tables may be switched off by setting W2NUM and W2STF to false 18The two integer arguments 0 0 in QNFILW are irrelevant and kept only for reasons of backward compati bility
56. ring a large kinematic range down to x 107 and 4 lt Q lt 10000 GeV see section 5 It is important to notice that the Qcdnum program evolves parton momentum densities in contrast to the formalism presented in the previous sections which is given in terms of parton number densities The user has thus to supply starting values for zg z zX z instead of g x E x etc 13 3 2 Evaluation of convolution integrals When written in terms of the parton momentum density H x xh x the convolution integrals in eqs 6 8 and 18 21 take the following form 1 EP 2 nt 30 z2 zZ where P x is a QCD splitting or coefficient function The distribution H x is sampled on a grid Xo lt T1 lt lt En lt Tu 1 with H x444 H 1 0 This distribution is linearly interpolated between grid points t x zi epo Ti Inserting eq 31 into eq 30 the convolution integral can be written as a weighted sum n I xo 2 Uleta 32 with w o o S s So Ww Ti o S1 i41 Si S2 Si Si 1 33 where s x and Si u v XA z u P d 34 Pd So u v fF z v a UV uU Ju In Qcdnum the weights are calculated at initialization from Gauss integration of eq 34 for all LO and NLO QCD splitting and light quark coefficient functions and are stored in 2 dimensional tables In addition tables based on analytical expressions for the LO weights are available they
57. s 9 Notice that the authors of the set may have used a different a evolution algorithm than Qcdnum This can lead to small inconsistencies 40 CPU sec CPU evol Routine calls evols EVOLNM 1000 431 6 EVOLNP 1000 431 6 EVOLSG 500 215 8 AP total 2500 1079 0 STFAST 500 250000 0 QNFILW 0 Other Total Table 8 Example of a Qcdnum timelog normalized to the size of the z Q grid the CPU time spent and the number of CPU seconds per evolution Also given is the time spent in fast structure function calculations and that in creating the weight tables none in this example since the weights were taken from disk The entry other gives the time spent in any other Qcdnum routine e g QSTFXQ or user routine e g Minuit between the moment QNTIME was called and the moment of the printou Splitting and coefficient functions can be accessed by val QNSPLF opt x q2 nf 30Parton distributions are often evolved over part of the grid only 31The second column gives the number of structure function calculations f you want to continue logging after the printout please call QNTIME with argument Start start a new log or Continue continue with the present log Option Function Option Function PFFi P cig CiD s PFG1 Pad a f cie Cha f PGF1 PR coQ CL z pect P x cac O s PFF2 P f cLQ CY 2 PFG2 PER f cre CR a f PGF2
58. s can of course also be reduced in size but we prefer to give complete freedom to choose the grid density at both low and high x by allowing a non equidistant grid in In zx An equidistant logarithmic r grid with larger spacing at high x does not much affect the heavy quark structure functions since these are strongly suppressed by threshold effects at large x 15 Condition ii In previous steps of the evolution both H x t 1 and H a t 1 have been calculated for tj lt tj then also eq 38 is a linear equation with two unknowns H and H 4 in shorthand notation 1 Hj Hat HG Ay 39 Solving eqs 37 and 39 for H we get 2H 4 Hi S A 40 Because of conditions i and ii above all quantities on the r h s of eq 40 are known The value of H is found by substituting eq 40 back into eq 37 In this way both H z t and H t are calculated for a given grid point x t by solving two linear equations with two unknowns Next we show that it is possible to define starting from a suitable grid point a progression in x and t such that the conditions i and ii mentioned above are satisfied at each step of the evolution First condition i is satisfied for all t at the highest grid point x since we have defined n 1 1 and H z 41 t H 1 t 0 see section 3 2 Second at the starting value of the evolution tg H x tg is known input by the user and H x to can be calculated direc
59. s follows implicit double precision A H 0 Z dimension a 3 data a 1 01 0 35 5 22 GLUONS a 1 x a 2 1 x a 3 22 return end 4 4 Initialization Qcdnum is initialized by call QNINIT This routine which must be called before anything else initializes a large amount of constants and sets reasonable defaults call QNVERS version Ldouble nxmax nqmax returns the version number in the character 8 variable version and sets the logical variable Ldouble to true if the current program version is in double precision and to false if it is in single precision The variables nxmax and nqmax are set to the maximum number of grid points in z and Q the user can define A timelog of Qcdnum subroutine calls is started by call QNTIME Start At any moment you can interrupt the logging by a call to QVTIME Hold and resume it again by calling QNTIME Continue 4 5 Qcdnum variables and options A list of Qcdnum variables options is shown in table 3 The meaning of these variables is as follows e W1ANA WTFLB these eight logical options steer the weight calculations as described in section 4 7 e LIMCK if set to true one cannot access parton distributions or structure functions outside the grid boundaries or cuts See section 4 12 e CLOWQ if set to true one can access heavy quark structure functions only for Q gt 1 5 GeV If set to false no such check is made See section 4
60. should be calculated from a non singlet quark density and not from something else The variable iflag is set to 1 if Q or any other scale like 1 3 or u see section 4 13 6 An attempt to calculate any structure function from the gluon distribution will result in a fatal error Apart from this Qcdnum accepts any combination of opt and name 36 runs outside the grid boundaries or cuts to 0 if the structure function calculation is successful and to 1 if the result is obtained from fast calculations see section 4 14 The heavy quark contribution to Fy and FT can be calculated with QSTFIJ or QSTFXQ by setting the input option opt to F2C FLC F2B or FLB provided of course that the heavy quark weight tables are calculated beforehand There are some restrictions however we refer to section 2 5 and make a few remarks below e The quark distributions making up the target nucleon and the gluon distribution must have been evolved with f 3 flavors for all Q In fact a QCD analysis which uses heavy quark structure functions must have Q and Q set beyond the highest Q grid point to guarantee that Qcdnum uses f 3 in all calculations except as e The heavy quark structure functions can only be calculated for the neutral current pro cesses 7y A X where A stands for the target nucleon i e p proton n neutron or d n p 2 deuteron Qcdnum cannot check this so that it is the responsibility
61. sions for non singlet structure functions and quark densities Convolution integrals are indicated by the symbol 10 2 5 Heavy flavor contributions In Qcdnum there are two ways in which heavy flavors can be taken into account 1 Treat the heavy quarks h c b as light quarks but allow them to contribute only above a given threshold in Q In the following we will denote these thresholds by Q and Q for charmed and bottom quarks respectively This prescription is known as the light quark variable flavor number scheme 2 Do not introduce any thresholds but consider only three light flavors the parton distri butions are evolved with f 3 and the structure functions F gt or Fr are calculated from eq 20 The heavy quark contributions F or FP are calculated from the photon gluon fusion process including NLO corrections see below This prescription is known as the heavy quark fixed flavor number scheme In the following we will describe both schemes in more detail and point out a few Qcdnum do s and dont s To describe the first approach let us assume that we want to calculate the P5 structure function in charged lepton proton scattering from parton distributions given at an input scale Q2 We choose for the singlet non singlet decomposition of the proton quark distribution those given in section 2 2 eqs 15 16 17 for f 3 4 5 flavors The range of the evolution and the input scale Q are set such
62. tes it at any given z and Q If Q lt Q Qcdnum uses the factors defined for f 3 if Q lt Q Q those for f 4 and if Q lt Q those for f 5 see section 4 10 on how the flavor thresholds are set 30 And for five flavors 11 1 1 1 quc geo wear er Ld factors 1 11 45 factors 2 1 6 factors 3 1 6 factors 4 1 6 factors 5 1 6 call QNLINC 11 proton 5 factors The function id IPDFID name returns the identifier id for a given name so that you do not have to remember all identifiers e g Define proton distribution for 3 flavors factors IPDFID singl 4 18 factors IPDFID umind 1 6 factors IPDFID splus 1 6 call QNLINC 11 proton 3 factors It is guaranteed that 1 lt id lt 10 because Qcdnum abends with a fatal error message if the name passed to IPDFID does not correspond to a memory resident quark distribution A list of all parton distributions can be written to logical unit number lun by a call to QPRINT lun Booklist This gives for the example given above the list as shown in ta ble 5 Notice from table 5 that Qcdnum has assigned trivial multiplication factors to the memory resident distributions gluon xqval and also that all names are truncated to five characters and converted to upper case 4 4 Ww 1 W2W 3W 4 W5
63. that Qu lt Q lt Qi lt Qi lt Ce Input to the calculation are the parton distributions g x QS D x i Ault Qj and qi v QU Furthermore since qe 0 for Q lt Q and q 0 for Q lt Q we set the starting values of qt and qj at aH Q2 39 Q2 and qj v QD 9 QD The analysis now proceeds as follows e Evolve the singlet gluon distribution from Q2 down to Q2 and up to Q2 In the upward downward evolution the number of flavors changes from 4 to 5 3 at Q Q2 Qcdnum takes automatically care of this e Evolve A 4 here also Qcdnum takes care of the number of flavors e Evolve q this is more complicated because q7 is discontinuous at the thresholds We therefore evolve from Q up to Q add 1 4 1 5 X x Q2 to qf x Q2 and continue the evolution from Q up to Q Likewise for downward evolution where we add 1 4 1 3 X x Q2 at the threshold Q e Evolve q7 here we evolve from Q upwards to Q taking into account the discontinuity at Q as described above Downward evolution does not make sense because we have 11 assumed that charm does not contribute below the threshold Q Thus you should either take care that q below the charm threshold is not used in subsequent calculations or you have to set it explicitly to d sQ imemqQ Q Q e Evolve qj similar to q7 except that we evolve upwards starting from Q e Calculate the F structure function for any x and Q withi
64. the proton distribution for all j can be obtained from a singlet gluon evolution together with that of only one q type non singlet provided the threshold discontinuities in the latter are taken properly into account See section 2 5 for further details on heavy flavor thresholds Below we give another decomposition of the proton quark distribution 4 1 1 lp pa EE Deum 15 I Bo p 6 15 5 1 1 1 ME E y EE ges erates 16 11 1 1 1 1 M a ES a a q 17 qi ip ghut gh Ge gd 17 Such a decomposition is not very economic since it requires up to four non singlet evolutions It is given here only because it will be used in examples later on 2 3 Structure function evaluation In NLO in the MS scheme the relation between the singlet quark distribution and the structure functions F and F F 2aF is given by 7 F a Q k2 qm Q oe dz of 2 26 9 c 2 se 95 a8 z Z z x with k 2 L Here we have introduced the physical scale Q which in deep inelastic scattering is the negative square of the four momentum transferred from the initial lepton to the target In eq 18 the renormalization scale u see section 2 1 the mass factorization scale 2 see section 2 2 and Q are assumed to be equal uw piy Q We have a b in eq 11 so that the non singlet valence distributions q7 do not contribute to q See section 2 4 for when they are chosen
65. thresholds Get number of flavors Get a Q Set value of parton density Add factor times singlet Set parton distribution to zero Singlet gluon evolution Non singlet evolution rA Non singlet evolution xq Non singlet evolution xq Get value of parton density Get value of parton density Calculate Fo Fr or xF Calculate P5 F or xF Mark for fast F calculation Clear marks Fast P calculation Clear memory allocation Print Qcdnum info Table 2 Subroutine and function calls in Qcdnum 20 A variable with a name starting with the letter L is of type logical A variable with a name starting with the letter I N excluding L is of type integer e We denote character variables by giving the name in quotes e g name Significant characters of input variables are given in capitals as in e g call QPRINT 6 Timelog Qcdnum is case insensitive so that T t time and Timelog are all valid inputs e All other variables are either of type real or of type double precision depending on the Qcdnum version which is used 4 3 Example program In this section we give an example of a simple application program The program evolves the singlet quark and gluon distribution given at a starting value of Q2 7 GeV in NLO NLO is Qcdnum default on a 80 x 60 z Q grid The grid covers the kinematic range 4x 1074 lt z lt 1 and 4 lt Q lt 5000 GeV The flavor thresholds are set to Q Q 4 25 GeV
66. tine just marks the four grid points surrounding each data point x Q The rou tine QFMNUL clears all the marks The grid points have to be marked only once somewhere in the initialization phase of the user program but of course after the z Q grid has been defined 2 The call call STFAST opt name calculates the structure function for all marked grid points As before opt can be set to F2 FL and name to a valid name of a memory resident distribution or linear combination see section 4 12 STFAST does of course the same as QSTFIJ but the code is more efficient STFAST stores the results for later interpolation in 20 memory locations reserved for this purpose If more than 20 different structure functions are calculated i e more than 20 different combinations of opt and name STFAST runs out of storage space and acts as a do nothing Qcdnum then switches automatically to slow structure function calculations 3 Finally the value of the structure function can be obtained by calling QSTFIJ or QSTFXQ as described in section 4 12 These functions will if possible pick up the results calculated by STFAST and perform the interpolation if any if not possible the structure function is calculated from scratch 4 In the unlikely event that you want to calculate more than 20 structure functions with STFAST you can after step 3 clear the memory allocation by calling STFCLR and repeat steps 2 and 3 for a ne
67. tions up to NLO Calculation of F5 Fr and zF up to NLO Calculation of the heavy quark contributions to F and Fr up to NLO Possibility to vary the mass factorization or renormalization scale In NLO the calculations are carried out in the MS scheme Qcdnum has a long history The program was originally developed by members of the BCDMS collaboration 1 and was used for instance in the QCD analysis of SLAC BCDMS F structure functions 2 Qcdnum was later on upgraded for use at low x by the NMC 3 At this stage the program could only run on vectorized machines like the CRAY A complete revision for the QCD analysis of ZEUS F data 4 led to fast devectorized code which can run on any machine During the 1995 1996 Hera workshop results from Qcdnum were compared to those from several other NLO evolution codes and found to be in agreement to within 0 0596 5 The present version of Qcdnum can be regarded as a black box which performs numerical calculations with reasonable speed and accuracy A user interface allows to interact with the program A QCD analysis is therefore not performed by Qcdnum itself but by user code instructing Qcdnum what actions to take This report is organized as follows in section 2 we give a brief outline of the formalism un derlying a QCD analysis of structure functions Section 3 is devoted to the numerical method used by Qcdnum Section 4 user manual describes how to use Qcdnum in a QCD analysis
68. tly from the QCD evolution eq 36 At the next grid point t gt t condition ii is thus satisfied for all values of x It follows that both conditions i and ii are met at the grid point z t1 which serves as the starting point of the evolution It is easy to see that one can subsequently progress to larger t at fixed values of x or to smaller x at fixed values of t and continue the evolution over the whole z t grid Here we have tacitly assumed upward evolution in t gt to the method applies to downward evolution as well provided tj is replaced by t in the above It is straight forward to extend the evolution algorithm to the singlet gluon case by solving four linear equations with four unknowns at each evolution step Finally we remark that the method outlined in this and the previous section yields parton distributions which are represented by second order splines in t at fixed values of x and by first order splines in x at fixed values of t The evolution is fully numerical in the sense that no use is made of specific parametric forms of the parton densities 4 Qcdnum user guide The Qcdnum program consists of a set subroutines which perform the QCD evolution of parton distributions and o and which calculate structure functions from the evolved distributions Interface routines allow to interact with the program which from a user point of view can be regarded as a black box A QCD analysis is performed by an appl
69. to 20 linear combinations of the quark densities can be defined Setting of options parameters e g e input of the flavor thresholds e input value for o e LO or NLO calculations e etc Input of the parton densities for each grid point in x at some fixed value of Q Q followed by the QCD evolution of these distributions Routines are provided for the singlet gluon evolution and for the evolution of non singlet quark densities xq xA and Eq Evaluation of structure functions and in fitting applications the calculation of the x7 Qcdnum provides a routine for the calculation of F5 Fr and zF3 The contribution to F gt and Fr from charm and bottom can also be calculated In the following subsections we give a detailed description of the steps outlined above For quick reference a list of Qcdnum subroutines is given in table 2 In the remainder of this write up we use the following notation 19 Subroutine or function Description QNINIT Initialization QNVERS vers Ldoubl nxmax nqmax Get info about current version QNTIME option Start stop CPU timelog QNiSET var ival i L I R Set Qcdnum variable QNiGET var ival i L I R Get Qcdnum variable GRXDEF nx xmi Define the x grid GRQDEF nq qmi qma Define the Q grid GRXINP xarray nx Add x grid points GRQINP qarray nq Add Q grid points GRXNUL GRQNUL GRMXMQ nxmax nqmax GRGIVE nx xm
70. u of flavor i d u s for i lt 0 it denotes the corresponding anti quark distribution q r u and for i 0 it denotes the gluon distribution g x u In the quark parton model and also in leading order LO perturbative QCD the parton density distributions are defined such that h x u dx is at a given u the number of partons which carry a fraction of the nucleon momentum between x and z4 dx The distribution zh x u is then the fractional momentum density Beyond LO QCD no such intuitive interpretation of parton distributions is possible they become theoretical constructs and their definition depends on the renormalization and factorization scheme in which the calculations are carried out MS in Qcdnum If the x dependences of the parton densities are known at some fixed value of pz ug they can be evolved to any value of u using eq 4 which consists of 2f 1 coupled integro differential equations Using symmetries in the splitting functions see e g 7 a much more simple set of equations can be written for the evolution of the flavor singlet and flavor non singlet quark distributions The singlet quark distribution is defined as the sum of all quark and anti quark densities f gt ai x Gi w 5 i 1 and obeys an evolution equation which is coupled to that of the gluon mel z t mem Qe A EE a Far D z Paa s There are two types of non singlet distributions qd x qi
71. ue false if the current version is in double single precision Should be declared logical in the calling program e nxmax nqmax the maximum number of grid points in z and Q the user can define call QNTIME option Start Qcdnum timelog option Start stop the log Hold or continue with the current log Continue A printout see QPRINT stops the current timelog call QNtSET tvar tval e QNLSET Lvar Lval set the logical Qcdnum variable Lvar to Lval e QNISET Ivar Ival as above but for integer variables e QNRSET C Rvar Rval as above but for real double precision variables See table 3 in section 4 5 for a list of variables you can set call QNtGET tvar tval As QNtGET but now retrieve the value tval corresponding to the Qcdnum variable tvar call GRXDEF nx xmi Clear the current x grid and generate a grid covering the range min lt x lt 1 with nx points The grid is logarithmic linear below above the value scaxo default 0 2 can be reset by a call to QNRSET call GRQDEF nq qmi qma 2 Clear the current Q grid and generate a logarithmic grid covering the range Q2 lt Q with nq points call GRXINP xarray nx Add nx points as specified in the input array xarray to the x grid The list of grid points may be given in arbitrary
72. ulated with the assumption that the parton densities can be linearly interpolated at all Q from one z grid point to the next With this assumption the convolution integrals can be evaluated as weighted sums The weights which are essentially integrals of the QCD splitting functions are calculated to high precision at program initialization From the value of a given input parton distribution and the calculated slopes at QZ the dis tribution can be evaluated at the next grid point Q gt Q2 or Q lt Q3 This distribution then serves as an input to calculate the slopes at Q and so on In this way the evolution can be continued over the whole z Q grid In Qcdnum the evolution is based on a second order spline interpolation in In Q In earlier versions of the program 1 this quadratic interpolation was achieved iteratively starting from a linear evolution in In Q from one Q grid point to the next In the present program this is achieved without iteration by solving linear equations at each evolution step The number of grid points in x and Q has to be chosen with some care A larger number of grid points improves the numerical accuracy the interpolations underlying the evolution algorithm become better approximations but cannot be taken too large because the CPU time increases quadratically linearly with the number of points in x Q In practice good balance between accuracy and speed has been achieved on a 100 x 40 z Q grid cove
73. values can be retrieved by QNRGET tchrm q2c and QNRGET tbott q2b respectively nf NFLGET iq Get the number of flavors associated with the grid point iq If ig is outside the grid boundaries NFLGET returns with a fatal error alphas QALFAS q2 qlambda nf ierr Calculate o Q The calculation depends on the input value of o Q2 at some reference scale Q the quark mass thresholds and the order of the Qcdnum calculations All these parameters are set at initialisation and can be changed by calls to QNISET QNRSET see section 4 5 Apart from a the function also returns value of the QCD scale parameter A the number of active flavours f and an error flag ierr which is set to 0 1 if Q is above below A If Q A ierr 1 then arbirarily QALFAS 100 and qlambda 0 call QNPSET name ix iq value Set the value of the parton distribution name at the grid point ix iq A fatal error occurs if name does not correspond to that of a memory resident distribution if the z Q grid is not defined or if ix iq is outside the grid boundaries call QADDSI name iq factor Add factor times the singlet distribution to a non singlet distribution The addition is done at fixed iq for all grid points in x See section 4 11 on how QADDSI is used to add discontinuities of xq distributions at the flavor thresholds Error messages as for QNPSET
74. vy quark weight tables typically takes a few CPU minutes only 4 8 Disk dump read of weight tables To avoid the weight initialization for every Qcdnum run there is the possibility to write them to disk QNDUMP and read them back again QNREAD in a next pass of your QCD analysis The rules are as follows e Unformatted read write this implies that weight files cannot be exchanged across differ ent machines e Only those tables which are available at the point where QNDUMP is called are written to disk e For each table the full common block is dumped This implies that the size of the output file does not depend on how many z Q grid points are actually defined but on the size of the common block itself This depends on the parameters mxx and mq2 which are also written to the file Weight files cannot be exchanged across Qcdnum versions which have different values of mxx and or mq2 see also section 4 1 e All weight tables depend on the current x grid definition whereas the heavy quark tables depend in addition on the Q grid these grid definitions are also written to disk e The heavy quark tables depend on the charm and or bottom mass which are dumped too The weight tables are written to disk by a call to QNDUMP e g OPEN unit 24 file weights file form unformatted status unknown CALL QNDUMP 24 Files containing all weight tables including heavy quarks typically require about 5 Mbyte of storage To read the fil
75. xt batch of 20 structure functions A call to QPRINT 1un Stats prints a log as shown in table 7 From this log it is seen that almost all structure functions are obtained from STFAST and that no structure functions are calculated outside the grid or cuts To summarize you can speed up the structure function calculations by simply marking grid points using QFMARK and calling STFAST before the calls to QSTFIJ or QSTFXQ See section 5 on the performance Structure function calls F2 FL xF3 F2h FLh Slow calculation 577 0 10 624 0 Fast calculation 1423961 0 93241 1425183 0 Outside grid or cuts 0 0 0 0 0 Total 1424538 0 93251 1425807 0 Table 7 Qcdnum log of the number of structure function calls 39 4 15 Structure functions from external parton distribution sets As mentioned in section 4 12 Qcdnum can be used to calculate structure functions from any given parton distribution set These parton distributions can be obtained from e g Pdflib and stored in Qcdnum memory using the routine QNPSET Doing so you should keep the following in mind e The parton distribution set should have been defined in the MS scheme e The order LO or NLO the flavor thresholds and a in Qcdnum should be set identical to those of the parton distribution set e The quark distributions stored in Qcdnum memory should either be singlet xX or non singlet Composite objects like the proton can be defined using the routine QNLINC see section 4 9 It is fo
76. ys Lett B309 1993 222 ZEUS Collab M Derrick et al Phys Lett B345 1995 576 J Bl mlein et al Proc workshop Future Physics at HERA eds G Ingelman A de Roeck and R Klanner Vol 1 p 23 1996 V N Gribov and L N Lipatov Sov J Nucl Phys 15 1972 438 675 L N Lipatov Sov J Nucl Phys 20 1974 181 G Altarelli and G Parisi Nucl Phys B126 1977 297 W Furmanski and R Petronzio Z Phys C11 1982 293 G Gurci W Furmanski and R Petronzio Nucl Phys B175 1980 27 W Furmanski and R Petronzio Phys Lett 97B 1980 437 G Ingelman and R R ckl Z Phys C44 1989 291 J Bl mlein et al Z Phys C45 1990 501 R G Roberts The Structure of the Proton Cambridge University Press 1990 A D Martin W J Stirling and R G Roberts Phys Lett B266 1991 173 J Bl mlein et al Nucl Phys B Proc Suppl 51C 1996 97 S Riemersma J Smith and W L van Neerven Phys Lett B347 1995 143 and references therein A Buras Rev Mod Phys 52 1980 199 54

QCDNUM16: A fast QCD evolution program

Contents

Download Pdf Manuals

Related Search

Related Contents