FDCHQHP:A Fortran Package for Heavy Quarkonium HadroProduction

FDCHQHP is a Fortran package to calculate the transverse momentum (p_t) distribution of yield and polarization for heavy quarkonium hadroproduction at next-to-leading-order (NLO) based on non-relativistic QCD(NRQCD) framework. It contains the complete color-singlet and color-octet intermediate states in present theoretical level, and is available to calculate different polarization parameters in different frames. As the LHC running now and in the future, it supplies a very useful tool to obtain theoretical prediction on the heavy quarkonium hadroproduction.


I. INTRODUCTION
From the discovery of heavy quarkonium in the 1970s, both the theoretical and experimental study on production and decay of J/ψ and Υ, plays an important role in the research on the nonperturbative and perturbative aspects of QCD. The factorization framework of non-relativistic QCD (NRQCD) had been proposed to study the production and decay of heavy quarkonium [1]. By extracting the NRQCD long-distance matrix elements (LDMEs) from the matching between theoretical prediction and experimental data, the NRQCD calculation had given a good description on the transverse momentum distribution (p t ) of heavy quarkonium production at hadron colliders. But the polarization of heavy quarkonium hadroproduction had been an open question for many years.
In the past, many works has been done by applying the Fortran codes generated by our Feynman Diagram Calculation(FDC) package [2] which is written in REDUCE and RLisp and can generate Fortran codes automatically. In the recent two works [3,4], we have presented the complete next-to-leading-order(NLO) predictions to the p t distribution of yield and polarization for J/ψ and Υ prompt hadroproduction. It took us a lot of efforts to manipulate running of more than ten thousands the numerical calculation jobs and control the final numerical precision when more than hundred numerical results from sub-processes must be sum up for each polarization result at each p t point. Finally, based on the collection of all the Fortran codes generated by using FDC, the efforts and experience on job submission and numerical precision control are implemented in an independent package FDCHQHP, it can be easily run by inexperienced user to obtain any desirable results on heavy quarkonium hadroproduction at any given experimental condition. It is well known that only prompt experimental measurements on heavy quarkonium hadroproduction are available up to now. Therefore, to fit the experimental measurement, the theoretical prediction must be on the prompt production. FDCHQHP supplys the complete calculation on the p t distribution of yield and polarization for prompt hadroproduction, which includes direct production part and feed-down contribution form excited states production.
From another aspect, as the LHC running at the current energy, the p t distribution of yield are measured up to large p t region (e.g. the measurement of Υ by ATLAS are up to 70GeV [5]), and there are more polarization parameters measurement e.g. λ θ , λ θφ , λ φ in different reference frame for heavy quarkonium polarization, and also the frame-independent polarization parameter λ [6]. What's more, in the next years, the LHC running energy will upgrade to 14 TeV, and there will be more high-precision measurements to be published in the future. Although we have already finished a few theoretical prediction in comparison with the experimental measurements in past two work. There are still a lot of experimental measurements in future need the theoretical predictions which needed large scale calculation in supercomputer environment. There are already many request for theoretical prediction in different experimental condition, we are ready to make the package public so that it can be used to obtained theoretical prediction by other people.
This package includes 6 channels and 76 sub-processes as summarized in above  [8,9] are included in the package for convenience. Please note that, if the file makefile was changed, you should execute make clean first, and then execute make, otherwise, just make is OK.

III. THE DIRECTORY STRUCTURES AND PROGRAMS
A.

Structure of directories
All the Fortran codes are collected in the folder FDCHQHP. The folder's structure is like : • source/<channels>/<sub-processes>/fort/ • fdc2.0/f77/src/src64d ifort(src64d ifort p100, src64q ifort)/ , respectively. The <sub-processes> under each <channel> are the ones listed in Tab. I, they are named with 'g'or 'q', which denotes the partons in initial and final states except the heavy quarkonium. For example, gggg denotes g + g → J/ψ(Υ) + g + g, qggq denotes g + q(q) → J/ψ(Υ) + g + q(q), and so on. Without specification of quadruple precision, all the calculation are double precision calculation. Especially, in the channel named 3pj 8, the sub-process ggg loop32, denotes the one-loop part for g + g → J/ψ(Υ) + g(one-loop) with quadruple precision calculation Fortran source involved. A corresponding sub-process, ggg other denotes the remaining part aside from the loop part(i.e. ggg loop32 ) with double precision calculation.
We know that the physical processes are divided into several channels, and a channel are divided into a few sub-processes. Therefore, the final result in physics should be the summation of all of them. For any sub-processes, the Fortran source codes and the corresponding compiled executable program (i.e. int) are placed in the fort directory.

B. Programs
In the fort directory, there are several kinds of files, and the most important one for users is input.dat which present all the input to run the Fortran program. We describe them in the following.

Introduction to the file input.dat
For the direct production part, including 1 S 8 0 , 3 P 8 J , 3 S 1 1 , 3 S 8 1 , as well as for the feed-down contributions from χ c in the directories chic-3pj 8 and chic-3s1 8, the file input.dat are completely same, just as Tab.II.
As an example, when we calculated the p t distribution of polarization parameter(α) for J/ψ prompt production under a certain experimental measurement (here denotes the rapidity cut) at the LHC [3], we prepared two input.dat files which uses the format showing in Tab.II with me2 type set as 4 for longitude polarization and 1 for transverse polarization; Other parameters can be easily changed according to the detailed description in Tab.II.

Introduction to other files
Users can skip over this block directly, if they do not concern the details of other files.
• The files prefixed with amp, are the amplitude calculation part the ampN.f is the amplitude for the Nth tree-diagram, amplN.f is that for the Nth loop-diagram, amplNn2.f, amplNn3.f and amplNn4.f are the 1 ε IR , 1 In amps2.f file, it calls all the diagrams contribution and squares the matrix elements, there are different blocks for the tree-diagram part, the loop-diagram part, the remaining term from divergence canceling part. It also contains the different choices for polarization direction and frame.
• The files prefixed with incl contain the common blocks and the definition of data, used in the Fortran programs.
• makefile, list the rules, libraries with location to link and compiler for building the executable file int.
• parameter.f and parameter1.f are two files defining the physical parameters in calculation.
• func.f, generate a certain phase space point with a set of random number and calculate the cross section at the point.
• int.f, the main program, read all input and calls the Monte Carlo program for phase space integration.
• genppp.f, and coupling.f calculate all the Lorentz invariable and constants in the programs.
• bases.data and the files prefixed with plot, are generated automatically in calculation, which are usually not used.
• input.dat, lists the input parameters and provides an interface for users.
We make two tables in this paper to explain it as Tab. II.
• The files fresult.dat and convergence.dat, are the results we expect after programs running. The former is used in the data processing, and the latter is the precision analysis for the integration. They will be replaced by new ones if the executable int runs again. We will give some detailed introduction in the section OUTPUT FILES.
• Other files not listed.

A. Calculation
As we know, a complete prompt result should be the sum up of all the channels with fixed LDMEs. Users can also do the calculation for one sub-process independently. The following is the calculation introduction.

calculation for a sub-process
Once the programs compiled successfully, the executable file int will be found in the fort folder. Two calculation methods are available: It may take a few seconds to a few hours depending on the process if just one CPU is used. Usually, the loop diagrams will take more time. No matter which method is used, the input.dat and cteq6m.tbl(or other PDF used) files should be placed in the calculating folders correctly.

calculation for a channel or channels
Go into the source directory, execute ./run-all to do the calculation for all channel. It may take users several hours to a few days depending on the calculation precision and computer.
By modifying the variable $op in line 4 with any channel name in the file run-all, such as op = ′ 3s1 1/ * /f ort ′ , users can obtain the calculation result for the specified channel.
In the script run-all, we have processed the data collecting step. The results for different channels are listed in the corresponding folders, named as data *.

V. LARGE SCALE CALCULATION IN SUPERCOMPUTERS
To simplify the tedious manual work, we have developed a bash-shell script package WBIN to do the numerical calculation automatically on supercomputers. This package provides a very simply interface and some commands for users. It has been realized to automatically parallize the calculation, generate and submit jobs, analyze the accuracy of the results, fix the suitable number of CPU to each job, separately calculate the positive and negative part of a sub-process when needed, realize the test calculation(for 3 p t points) and the real calculation(for all p t points), check the state of each job, and finally collect all the data. It is customized for our package, but also could be adapted in any similar application which use BASES MC integral program.

A. Installation
The package can be downloaded together with the Fortran package, named as WBIN.taz.
There are two steps to install it: • execute tar -xzvf WBIN.taz ; • execute cd WBIN/ • add the absolute directory of WBIN to the default search path, e.g. add the line 'set path=($path $HOME/WBIN)' in ∼/.cshrc .

B. Usage
In the calculation of transverse momentum distribution of yield and polarization parameters of heavy quarkonium hadroproduction, all the results used the p t as the their x-coordinate and there are usually many p t points. The result at these p t points can be calculated independently. Instead of the sequential calculation mode for all the p t points at our personal computer environment, the calculation could be parallelled for all he p t points at the same time in supercomputer environment. Therefore, the parallelization scale will become very large dramatically. Finally, there are many jobs have to be submitted to finished in the parallel computing environment. And each job is for calculation on one polarization component, one sub-process at one p t point with one experimental condition.
We need to do a three p t points test running with low accuracy requirement for the result of each channel which is obtained by summing up all the sub-directories in this channel. In the test calculation, the number of sample point in the Monte Carlo integration is adjusted to obtained the required accuracy, as well as the number of CPU cores is adjusted to keep the calculation finish at two hours. With all the information obtained in the successful test calculation, we can arrange how to run program in each sub-process so that the accuracy requirement for the final result of each channel can be achieved. For each jobs, the number of CUP cores and number of sample point in the Monte Carlo integration are fixed by applying the information obtained in test calculation so that every jobs will be finished within two hours by themselves. Then the results for all p t points can be obtained in the real calculation stage.

preparation
The most important preparation is the method of submitting jobs, which depends on the job management system in supercomputers. In WBIN package, two files refer to this method: the file Submit and sub job. It needs users themselves to understand the method of submitting jobs in their supercomputer environment to modify it properly. Here we provide two examples which was applied by us at the Shenteng7000 and TH-1A in China.
• For Shenteng7000, we give the example as files Submit shenteng and sub job shenteng.
It adapts the LSF(Load Sharing Facility) HPC7.0 job management system at the Shenteng7000. There is no upper limit for the total number of CPU cores applied by users, but every job submitted is added in task queues together with other users' jobs.
• For TH-1A, we give the example as files Submit th-1a and sub job th-1a.
In TH-1A, it adapts the method of batch jobs. We can apply a certain number of CPU nodes within a upper limit, in these CPU nodes, we can arrange our small scale jobs.
Once the CPU nodes are allocated to a user, the user's jobs will monopolize the CPU nodes until all the jobs are over. It is called a non-preemptive scheduling system.

test calculation
For the usage of WBIN package, these steps will be done by manual operation or by shell script as the following shows.
• A: Make sure all of the sub-processes in folder FDCHQHP/source have been compiled successfully.
• B: Go into the directory source, execute echo $PWD > dir s.
• C: Make a directory any where which is specially used to do the calculation, e.g.
• D: Copy the files all-pwd, dir s, input.dat, cteq6m.tbl in source folder into the folder work.
• E: Go into the folder work, and use the template file input.dat to generate more input files.
In details, copy input.dat to other files named as inputX.dat with X can be any letters except '.', and then modify inputX.dat with the options you want according to Tab. II. The number of inputX.dat is not limited, such as inputlp.dat, inputtp.dat, inputdirectlp.dat ....
• F: If users want to calculate part of the total 76 sub-processes, the user can just keep the corresponding lines and delete others in file all-pwd.
• G: Execute init test job. The flow for this step: 0) generate the folders redo job, done job in work ; 1) generate the directories according to inputX.dat, all the X value will be the directory name in the folder work ; 2) generate the multiple sub-directories in directory X according to the file all-pwd, such as X/1s0 8/gggg/, and copy the inputX.dat file into X/1s0 8/gggg/input.dat; 3) generate the sub-directories in X/1s0 8/gggg/ according to p t points in the current input.dat and name them as tedir1, tedir2, tedir3..., for test calculation, and then generate input.dat in tedir1, tedir2, tedir3..., with just one p t point as the p t list at the end of input.dat; 4) link the file cteq6m.tbl in the work folder to tedir1, tedir2, tedir3...,;
• H: Submit jobs with executing Submit.
It will submit all the jobs listed in the file work/redo job/job list. Each line in the file denotes a job and looks like:  e) The step J: is also working for real calculation.

data collection
The test and real calculations make up the whole process for the calculation in supercomputers. Finally, execute final data collect to collect the data for all the p t points in all the sub-processes, and then give the final results for the 6 channels.
C. some detailed processing

control the calculation precision
Using the Monte Carlo integral program, the error estimate and the square root of integral points are in inverse ratio. Therefore, if we know the error for a certain number of sample points, it's easy to obtain the error for another number of sample points.
where N is the total number of sub-processes in a channel, c i is the central value for each sub-process, e i is the absolute error for each sub-process, and r i is the relative error for each sub-process. Therefore R is the relative error for a channel. Consuming the maximal relative error for a channel we desired is R 0 (the default value is 1%), when R > R 0 , we can easily obtain that in each sub-process, the relative error should be Using the new error r ′ i , the error r i and the number of sample points Ncall in the last calculation, we can get the new number of sample points Ncall ′ : Ncall ′ is the new number of sample points in the next calculation process. It may do several iterations until R < R 0 is achieved. This flow chart of these steps is shown in Fig. 1.
The calculation includes two steps. The first step is test calculation, in which we set a less precise error value r i = 1% for each sub-process, then do the calculation only at three p t points which are the minimal, intermediate, and the maximal number of p t points listed in input.dat. In this step, the relative error R maybe not meet the precision requirement for the result of a channel in the case of large number cancellation among different sub-processes in this channel. Then another iteration will work until R is small enough just as mentioned above. The second step is real calculation, in which all the p t points listed in input.dat will be calculated. Using the method of linear continuation for the number of sample points at the 3 p t points calculated in test step, we can obtain the numbers of sample points for all other p t points listed in input.dat. In this case, the calculation precision in the whole p t region could be guaranteed R 0 .

other details
• In some sub-processes, for the p t points in the region where the central value changes from positive to negative or in inverse, the precision will get worse in such region. In this case, we adapt to do the calculation for positive and negative part, respectively, and then sum up these two parts as the central value for the sub-process.
• In each sub-process, the number of sample points decides the run-time, so it is easy to estimate the number of CPU cores and the time needed in large number of sample points through the information obtained in small sample number case.
• When using the same number of CPU cores, run-time for different sub-processes are quite different. So we need to allocate the suitable number of CPU cores for each job to keep it can be finished in two hours.
• In WBIN package, we decide the priority to submit job by the number of CPU cores needed and then estimated run-time for the jobs.

a brief guide for the shell scripts
You can skip this subsection if you just want to run the package. It is presented here just for user who want to involved in modification of program in the package.
• Initialize the directories and the job list in test calculationinit test job, get test pt • Job Controllingjob control • Submit the jobssub job • Analyze the running state of jobanly job • Analyze the precisionanly precision In FDCHQHP package , the result for a channel "n" are summing up of all the subprocesses in this channel as: where i, j, k denotes the partons, σ a n is the result in the sub-process (corresponding to a) directories without the suffix '-m' in folder names, and σ b n is the result in the sub-process directories with '-m' suffix in folder names. The second part in the Eq.(5) only appears in the channels 3pj 8 and chic-3pj 1. It depends on the choice of µ Λ which usually take the value of m Q . Users can take different µ Λ values to study the dependence on NRQCD scheme, as we did in Ref. [4].
Considering the kinematics effect for the feed-down contribution from high excited states, could be performed after obtaining the result of a whole channel.

VII. OUTPUT FILES AND THE PRECISION ANALYSIS
As mentioned above, the files fresult.dat and convergence.dat are the output files needed.
The format offresult.dat is shown in Fig. 3, and the format of convergence.dat are shown in Fig.4 and Fig.5.
The following are some notes for convergence.dat: • In Fig.4, the first two lines are the beam energies for colliders; the 3rd line is the p t point, the 22nd shows the number of sample points (given and real used).
• The 38th and 41st lines show the maximal times for the grid optimization step and integration step, and the 39th and 42nd lines give the highest accuracy we expected in the two steps, usually they are set to be 0.01%.
• In Fig.5, the lines 45th to 49th lines give a statistical run-time for the calculation on one p t point. The lines from 60th to 69th are the result of each iteration in the grid optimization step. From left to right, they are the iteration number, the efficiency, the portion of the negative value, integration value, the relative error, the cumulative result for the data on the front(the absolute error), the total relative error, the time when each iteration finished.
• Lines 78th to 87th, are the reliable result in the integration step. In this step, after the grid optimization step, the results for each iteration are close to each other, and the final cumulative central value is the average of these data, the corresponding relative error(e.g. the 8 column of line 78, which is 0.077%) is the 1/ √ 10 times as the first line(e.g. the 8th column of line 87, 0.244% ).
There are two standards to ensure that the convergence behavior of integration is good enough: 1), the relative error in the last 10 lines of convergence.dat is average and small enough, no more than 1% (if the accuracy is as high as almost 0.01%, the less average result also could be accepted); 2), on the base of 1) achieved, the final cumulative accuracy in the last line, which is also the accuracy presented in fresult.dat in the form of absolute error, should meet the required precision when all the sub-processes are summed, since we know that large number cancellation usually appears in the summation of the sub-processes.
In fact, we also require that the relative error for each sub-process itself is no more than 1%( when it is larger than 1% , the relative error in each iteration must be more than 3%, that means the the integration is not reliable). Besides, for a certain sub-process, when some parameters are changed, such as the rapidity range, the p t points, the polarization direction, and so on, the convergence behavior of integration could change at the same time.
In this case, the number of sample points should be reset for the sub-processes to ensure a smooth result in the whole p t region, especially for the loop-diagram sub-processes such as 'g + g →J/ψ + g ', and the sub-process 'g + g →J/ψ + g + g ', since the value of these two processes are usually large and opposite, which could lead to large number canceling. Fig.4 and Fig.5 is a complete result for one p t point. If more p t points were listed in input.dat file, the results will be written into one convergence.dat in the order of p t points, as well as the file fresult.dat.

VIII. SUMMARY
We have presented the Fortran package FDCHQHP, which is used to do the calculation on the p t distribution of yield and polarization for heavy quarkonium prompt production in hadron colliders. The package contains all the sub-processes of CS and CO states up to NLO level, and it focus on the calculation of the polarization parameters.
This package has been fully developed in our past works with rather reliable performance and accuracy. We also provide a small package WBIN, which can help users to run FD-CHQHP in supercomputers easily. Some other improvement are in progress so that we can make it more automatically.