Calculating amplitude in python

Asked by Li on 2021-04-10

Hi,
    I followed https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/FAQ-General-4 to generate the amplitude for pp->ttbar process, as the "Full library of process in python".

1) I get the same result as the example, however after using some different momentum, like

p = np.array([[ 53.35411266, 0. , 0. , 53.35411266],
       [ 1146.1299388 , -0. , -0. , -1146.1299388 ],
       [ 261.02808531, -70.24765479, 58.6271451 , -196.95877203],
       [ 938.45596613, 70.24765479, -58.6271451 , -895.81705408]])

    The result "ans" is 0.0. This four momentum is from the lhe file of a pp->tt process.
    The zero result occurs for around half of the events for a 1 million event lhe file.
     Why the result is 0?

2)nhel = -1 # means sum over all helicity
    In the code the nhel is defined as above. Does it sum over the inital plus final state particle, or only final state particles?

3)
    In the matrix.f file,

      IF (NHEL.NE.0) THEN
        CALL M0_SMATRIXHEL(P, NHEL, ANS)
      ELSE
        CALL M0_SMATRIX(P, ANS)

    And in the subroutine M0_SMATRIXHEL,

      USERHEL=HEL
      CALL M0_SMATRIX(P,ANS)
      USERHEL=-1

    And in the subroutine M0_SMATRIX

      DATA USERHEL/-1/

    What does "DATA USERHEL/-1/" means? Is it a selection of the amplitude from the following:

      DATA (NHEL(I, 1),I=1,4) /-1,-1,-1, 1/
      DATA (NHEL(I, 2),I=1,4) /-1,-1,-1,-1/
      DATA (NHEL(I, 3),I=1,4) /-1,-1, 1, 1/
      DATA (NHEL(I, 4),I=1,4) /-1,-1, 1,-1/
      DATA (NHEL(I, 5),I=1,4) /-1, 1,-1, 1/
      DATA (NHEL(I, 6),I=1,4) /-1, 1,-1,-1/
      DATA (NHEL(I, 7),I=1,4) /-1, 1, 1, 1/
      DATA (NHEL(I, 8),I=1,4) /-1, 1, 1,-1/
      DATA (NHEL(I, 9),I=1,4) / 1,-1,-1, 1/
      DATA (NHEL(I, 10),I=1,4) / 1,-1,-1,-1/
      DATA (NHEL(I, 11),I=1,4) / 1,-1, 1, 1/
      DATA (NHEL(I, 12),I=1,4) / 1,-1, 1,-1/
      DATA (NHEL(I, 13),I=1,4) / 1, 1,-1, 1/
      DATA (NHEL(I, 14),I=1,4) / 1, 1,-1,-1/
      DATA (NHEL(I, 15),I=1,4) / 1, 1, 1, 1/
      DATA (NHEL(I, 16),I=1,4) / 1, 1, 1,-1/
      DATA IDEN/256/

4) In DATA IDEN/256/, why is it 256 but not 16?

Thank you!

Best,
Li

Question information

Language:
English Edit question
Status:
Solved
For:
MadGraph5_aMC@NLO Edit question
Assignee:
No assignee Edit question
Solved by:
Li
Solved:
Last query:
Last reply:

Hi,

I have tested the script below with version 2.9.3 and did not get any issue for that phase-space point.
Can you check the same script (adapted since it seems that you use an old version of MG5amC. (so you likely need to remove proc_id parameter). Note that I run with python3.9

2)nhel = -1 # means sum over all helicity
    In the code the nhel is defined as above. Does it sum over the inital plus final state particle, or only final state particles?

Sum over all possible helicities according to your process. so for g g > t t~ this is 2**4=16 combination
if you use polarised syntax like
generate g g{+} > t t~
it will sum over 8

3) What does "DATA USERHEL/-1/" means?
This is a fortran syntax to initialise that parameter USERHEL to -1, so this is the default value as long as that parameter is not modified explicitly.

4) DATA IDEN/256/, why is it 256 but not 16?
256= 8 (number of color of the gluon initial state) * 8 (number of color of the second gluon) * 2 (number of helicity of the first gluon) *2 (number of helicity of the second gluon)

This factor also include symmetry factor if you have identical final states, so that number is
512 for g g > g g and 1536 for g g > 3g

Note that one optimization include is that vanishing helicity are automatically removed from the computation when summing over all helicity. This could lead to weird behavior when changing benchmark or if the filtering is not done in the correct way.
You can force goodhel array to stay always on True to avoid that filtering and see if that resolve your zero result.

Here is the exact script that i have run:

import allmatrix2py
allmatrix2py.initialise('../Cards/param_card.dat')

def invert_momenta(p):
        """ fortran/C-python do not order table in the same order"""
        new_p = []
        for i in range(len(p[0])): new_p.append([0]*len(p))
        for i, onep in enumerate(p):
            for j, x in enumerate(onep):
                new_p[j][i] = x
        return new_p
p =[[ 0.5000000E+03, 0.0000000E+00, 0.0000000E+00, 0.5000000E+03],
    [0.5000000E+03, 0.0000000E+00, 0.0000000E+00, -0.5000000E+03],
    [ 0.5000000E+03, 0.1040730E+03, 0.4173556E+03, -0.1872274E+03],
    [ 0.5000000E+03, -0.1040730E+03, -0.4173556E+03, 0.1872274E+03]
    ]
for i in range(100):
        p2 =invert_momenta(p)
        nhel = -1 # means sum over all helicity
        procid = -1 #
        pdgs = [21,21,6,-6] #which pdg is requested
        scale2 = 0. #only used for loop matrix element. should be set to 0 for tree-level
        alphas=0.13

        ans = allmatrix2py.smatrixhel(pdgs,procid,p2,alphas,scale2,nhel)
        print('ans', ans)

p = [[ 53.35411266, 0. , 0. , 53.35411266],
      [ 1146.1299388 , -0. , -0. , -1146.1299388 ],
      [ 261.02808531, -70.24765479, 58.6271451 , -196.95877203],
      [ 938.45596613, 70.24765479, -58.6271451 , -895.81705408]]

p =invert_momenta(p)
nhel = -1 # means sum over all helicity
procid = -1 #
pdgs = [21,21,6,-6] #which pdg is requested
scale2 = 0. #only used for loop matrix element. should be set to 0 for tree-level
alphas=0.13

ans = allmatrix2py.smatrixhel(pdgs,procid,p,alphas,scale2,nhel)
print(ans)

p = [[ 53.35411266, 0. , 0. , 53.35411266],
      [ 1146.1299388 , -0. , -0. , -1146.1299388 ],
      [ 261.02808531, -70.24765479, 58.6271451 , -196.95877203],
      [ 938.45596613, 70.24765479, -58.6271451 , -895.81705408]]

p =invert_momenta(p)
nhel = -1 # means sum over all helicity
procid = -1 #
pdgs = [21,21,6,-6] #which pdg is requested
scale2 = 0. #only used for loop matrix element. should be set to 0 for tree-level
alphas=0.13

ans = allmatrix2py.smatrixhel(pdgs,procid,p,alphas,scale2,nhel)
print(ans)

# output
100 times ans 0.7192889343295779
followed by
1.271062614780256
1.271062614780256

Li (huangli-itp) said : #2

Hi Olivier,
    Thank you very much!

1) I have tested the script below with version 2.9.3 and did not get any issue for that phase-space point.

    You're right, this point gives result for pdg=[21, 21, 6, -6]. And I found that my mistake is that I forget to convert the momentum in this case and this momentum array happens to be a 4x4 matrix. And when I use the momentum from lhe file, I converted the momentum but there are around half events with different pdg order. When the pdg is [-3, 3, 6, -6] or [21, 21, -6, 6], the result is zero. I found that the Subprocess name in the mg5 output file is P1_gg_ttx and P1_uux_ttx, does it mean I should follow the subprocess name order?

2)Sum over all possible helicities according to your process. so for g g > t t~ this is 2**4=16 combination

     I see. So nhel is a flag? I thought it's value with physical meaning.

3)This is a fortran syntax to initialise that parameter USERHEL to -1, so this is the default value as long as that parameter is not modified explicitly.

     I see. I am not farmiliar with Fortran. So does
           USERHEL=HEL
          CALL M0_SMATRIX(P,ANS)
     change the default value -1 to the value HEL in the subroutine M0_SMATRIX?

Thank you very much again! This python version amplitude calculation is really useful!

Best,
Li

Hi,

> I found that the Subprocess
> name in the mg5 output file is P1_gg_ttx and P1_uux_ttx, does it mean I
> should follow the subprocess name order?

Yes the handling of the ordering needs to be done within your python script.
(doing that inside fortran is really annoying).
You should have a Fortran API providing you information about the supporting order for the particle
(as well as associated process id).
You actually also have an API to get the relation between the value of nhel and the helicity/polarization of each particles.

> I see. So nhel is a flag? I thought it's value with physical
> meaning.

if set negative then you sum over all polarization/helicity
if set positive you only compute one entry in the list of polarization.
You do have an api function to know which nhel value means in term of helicity for each initial/final state.

> I see. I am not farmiliar with Fortran. So does
> USERHEL=HEL
> CALL M0_SMATRIX(P,ANS)
> change the default value -1 to the value HEL in the subroutine M0_SMATRIX?

USERHEL is indeed a global variable, so yes.

> Thank you very much again! This python version amplitude calculation is really useful!

Glad that it is not only usefull for us via the re-weighting module.

Cheers,

Olivier

> On 11 Apr 2021, at 02:50, Li <email address hidden> wrote:
>
> Question #696494 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/696494
>
> Li posted a new comment:
> Hi Olivier,
> Thank you very much!
>
> 1) I have tested the script below with version 2.9.3 and did not get any
> issue for that phase-space point.
>
> You're right, this point gives result for pdg=[21, 21, 6, -6]. And I
> found that my mistake is that I forget to convert the momentum in this
> case and this momentum array happens to be a 4x4 matrix. And when I use
> the momentum from lhe file, I converted the momentum but there are
> around half events with different pdg order. When the pdg is [-3, 3, 6,
> -6] or [21, 21, -6, 6], the result is zero. I found that the Subprocess
> name in the mg5 output file is P1_gg_ttx and P1_uux_ttx, does it mean I
> should follow the subprocess name order?
>
> 2)Sum over all possible helicities according to your process. so for g g
>> t t~ this is 2**4=16 combination
>
> I see. So nhel is a flag? I thought it's value with physical
> meaning.
>
> 3)This is a fortran syntax to initialise that parameter USERHEL to -1,
> so this is the default value as long as that parameter is not modified
> explicitly.
>
> I see. I am not farmiliar with Fortran. So does
> USERHEL=HEL
> CALL M0_SMATRIX(P,ANS)
> change the default value -1 to the value HEL in the subroutine M0_SMATRIX?
>
> Thank you very much again! This python version amplitude calculation is really useful!
>
> Best,
> Li
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Li (huangli-itp) said : #4

Hi Olivier,
    This solved my problem.
     Thank you very much! This gives me another interactive way to access MG5 and gives me a chance to do the phase space integration by myself.

Best,
Li