An issue of MadGraph when calculating neutrino cross sections

Asked by Qiji Xin

Dear MadGraph experts,

It seems that I noticed a common issue of MadGraph when calculating neutrino cross sections.

Following the discussion on https://answers.launchpad.net/mg5amcnlo/+question/678937 . From the last couple posts, the neutrino DIS (p + νμ → μ + X) cross section calculated by MadGraph is only half of standard values (e.g., https://arxiv.org/abs/1106.3723).

The DIS is complicated. So then I've tried several processes between neutrinos and elementary particles. Surprisingly, I found MadGraph's results is usually ~half of standard value.

First,
I've tried the ve e- > ve e-. The xsec from MadGraph is very close to 1/1.5 of standard value, and the shape is slightly different at lowest energies.
For the runcards, I set e- to be at rest and removed all the cuts:
 set lpp1 0
 set lpp2 0
 set ebeam1 0.
 set no_parton_cut

Second,
I've tried vm e- > ve m-. The xsec from MadGraph is very close to 1/2 of standard value, the shape is slight different near threshold.
(The standard value should be sigma = GF^2*s/pi*(1- m_mu^2/s), where s = 2*Eν*m_e for electron at rest. )
For the runcards, same settings as above.

Third,
I've tried ve a > e- W+. The xsec from MadGraph is very close to 1/2 of standard value, in https://arxiv.org/abs/hep-ph/9709290 . Also, the threshold should be s=(m_W+m_e)^2≃ 6462 GeV^2, but MadGraph's threshold is ≃ 6468 GeV^2.
For the runcards,
 set lpp1 0
 set lpp2 0
 set ebeam1 0.25
 set no_parton_cut
(Here I set photon energy to be 0.25 GeV just for convenience because the value of s will be same as that of ebeam2.)

So any possible explanation for above? Since the difference is usually very close to 1/2 of standard value, maybe it's because MadGraph did a helicity average, 1/2, to neutrino (but neutrino should only have one helicity)? Or maybe there are some remaining cuts I didn't remove after "set no_parton_cut"?

Thanks so much!

Question information

Language:
English Edit question
Status:
Solved
For:
MadGraph5_aMC@NLO Edit question
Assignee:
No assignee Edit question
Solved by:
Qiji Xin
Solved:
Last query:
Last reply:
Revision history for this message
Olivier Mattelaer (olivier-mattelaer) said :
#1

By default beam are not polarised.

You can fully polarised one beam (in your case the neutrino one) by setting
beampol1 to +100 or -100 (depending of the polarization that you are looking for)

Cheers,

Olivier

> On 10 Jun 2019, at 19:27, Qiji Xin <email address hidden> wrote:
>
> New question #681325 on MadGraph5_aMC@NLO:
> https://answers.launchpad.net/mg5amcnlo/+question/681325
>
> Dear MadGraph experts,
>
> It seems that I noticed a common issue of MadGraph when calculating neutrino cross sections.
>
> Following the discussion on https://answers.launchpad.net/mg5amcnlo/+question/678937 . From the last couple posts, the neutrino DIS (p + νμ → μ + X) cross section calculated by MadGraph is only half of standard values (e.g., https://arxiv.org/abs/1106.3723).
>
> The DIS is complicated. So then I've tried several processes between neutrinos and elementary particles. Surprisingly, I found MadGraph's results is usually ~half of standard value.
>
> First,
> I've tried the ve e- > ve e-. The xsec from MadGraph is very close to 1/1.5 of standard value, and the shape is slightly different at lowest energies.
> For the runcards, I set e- to be at rest and removed all the cuts:
> set lpp1 0
> set lpp2 0
> set ebeam1 0.
> set no_parton_cut
>
> Second,
> I've tried vm e- > ve m-. The xsec from MadGraph is very close to 1/2 of standard value, the shape is slight different near threshold.
> (The standard value should be sigma = GF^2*s/pi*(1- m_mu^2/s), where s = 2*Eν*m_e for electron at rest. )
> For the runcards, same settings as above.
>
> Third,
> I've tried ve a > e- W+. The xsec from MadGraph is very close to 1/2 of standard value, in https://arxiv.org/abs/hep-ph/9709290 . Also, the threshold should be s=(m_W+m_e)^2≃ 6462 GeV^2, but MadGraph's threshold is ≃ 6468 GeV^2.
> For the runcards,
> set lpp1 0
> set lpp2 0
> set ebeam1 0.25
> set no_parton_cut
> (Here I set photon energy to be 0.25 GeV just for convenience because the value of s will be same as that of ebeam2.)
>
> So any possible explanation for above? Since the difference is usually very close to 1/2 of standard value, maybe it's because MadGraph did a helicity average, 1/2, to neutrino (but neutrino should only have one helicity)? Or maybe there are some remaining cuts I didn't remove after "set no_parton_cut"?
>
> Thanks so much!
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Revision history for this message
Qiji Xin (xinqijisuper) said :
#2

Dear Olivier,

Thanks for your reply!
Yes, adding "set polbeam2 = -100" for the neutrino fixed the cross sections at most of the energies.

However, there is still difference from the standard values near the threshold. And MadGraph's result is lower, or have higher threshold.
For example.
+) for e- vm > ve m-, threshold from MadGraph is correct but near threshold it gives lower cross sections than standard values.
+) for ve a > e- W+, the threshold should be s = (m_W+m_e)^2 ≃ 6462 GeV^2, but MadGraph's result is ≃ 6468 GeV^2.
So maybe there are some remaining cuts I didn't remove after "set no_parton_cut"?

I can send you the figures of comparing MadGraph's results with standard values if you need.

Thank you so much!

Revision history for this message
Qiji Xin (xinqijisuper) said :
#3

Dear Olivier and all,

Any possible solutions for that follow-up question?
The model that I use is sm-full, which includes the mass of every particle involved in my processes, so this couldn't be due to massless final state.

Thanks so much!

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) said :
#4

Hi,

> for ve a > e- W+, the threshold should be s = (m_W+m_e)^2 ≃ 6462 GeV^2, but MadGraph's result is ≃ 6468 GeV^2.

So if I look at your formula and use the W and electron mass: (80.419000**2 + 5.110000e-04**2) = 6467.2155612611205
So it is in agreement with the value that you quote from madgraph.

So I guess that your issue is either with the value of the W mass or the electron mass.
The handling of the Wmass is tricky in sm-full since the W mass is not a free parameter in this model but fixed by gauge invariance.
(and therefore fixed by aeW, Gf and MZ). If you want to have to be able to change directly MW,

You can either do
set complex_mass_scheme
which in addition to setting all the masses to complex number within the evaluation of the matrix-element (m+iGamma)
also changes the EW scheme to have aw to depend of MW and MZ (and therefore aW will be a complex number now)

You can also decide to fully define the EW yourself (but then you have the responsability to keep this coherent)
by setting
set ewscheme external

Obviously here if your break the gauge invariance relation, no one knows what will be the exact consequence in term of physics.
(and in that case, the result can change from one version of the code to the next)

Cheers,

Olivier

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) said :
#6

> Could you please refer me to a literature/document about how the W mass is obtained from aeW, Gf and MZ?

For the reference any text book about SM and/or fermi theory should have such formula.
Since most of them are not free, I will simply recomend wikipedia which certainly have such information.

>
But how does MadGraph use 80.419 GeV?

You can tweak the value of each parameter in the param_card.
As said the W mass is not a free parameter by default, so you need to change MZ or aEW or Gf to change his value.

> I've also tried another process in which there is no W final state, and MadGraph gives different threshold.
> a vm > vm e e
> ( I've set the photon energy to be 0.25 GeV, so s = Eν * 1 GeV. )
> The threshold is Eν = (0.511e-3+0.511e-3)^2 / 1 GeV = 1.04448e-6 GeV
> But what MadGraph gives is the threshold to be ~2.3e-6 GeV.
> Below is my script for this process. So anything wrong with what I did below?

 I will look at your mistake here on next wednesday. (but if you find it yourself in the meantime)

Cheers,

Olivier

> On 14 Jun 2019, at 22:42, Qiji Xin <email address hidden> wrote:
>
> Question #681325 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/681325
>
> Status: Answered => Open
>
> Qiji Xin is still having a problem:
> Dear Olivier,
>
> Thanks so much for your answer!
>
> Could you please refer me to a literature/document about how the W mass is obtained from aeW, Gf and MZ?
> And, for the (m_W+m_e)^2 ≃ 6462 GeV^2, I am using PDF value, which is experimentally measured, that m_W = 80.379 (± 0.012) GeV.
But how does MadGraph use 80.419 GeV?
>
> I've also tried another process in which there is no W final state, and MadGraph gives different threshold.
> a vm > vm e e
> ( I've set the photon energy to be 0.25 GeV, so s = Eν * 1 GeV. )
> The threshold is Eν = (0.511e-3+0.511e-3)^2 / 1 GeV = 1.04448e-6 GeV
> But what MadGraph gives is the threshold to be ~2.3e-6 GeV.
> Below is my script for this process. So anything wrong with what I did below?
>
> Thanks a lot!
> ======================================================================
> import model sm-full #!
>
> generate a vm > vm e- e+ #!channel
> display diagrams
>
> output
> launch
>
> set lpp1 0
> set lpp2 0
> set ebeam1 0.25
> set polbeam2 = -100
> set no_parton_cut
> set nevents 1000.
> set ebeam2 1.0445e-6
> launch
> set ebeam2 1.3e-6
> launch
> set ebeam2 1.6e-6
> launch
> set ebeam2 2e-6
> launch
> set ebeam2 2.3e-6
> launch
> set ebeam2 2.7e-6
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) said :
#7

Hi,

Sorry this is actually an issue of our code and not a mistake on your side.
The phase-space is so constrained that the code fails to find any phase-space points.

One of the reason is that the code choose a wrong ordering for the phase-space generation and therefore miss to correctly generate the phase-space.
This explains why you do not have that issue for
>> vm a > vm e e

Additionally, the code fails to identify such threshold in advance, which decrease the efficiency of the code.
I will fix that and hope that it will improve the efficiency of generation in the first case.

Cheers and thanks,

Olivier

> On 15 Jun 2019, at 11:06, Olivier Mattelaer <email address hidden> wrote:
>
>> Could you please refer me to a literature/document about how the W mass is obtained from aeW, Gf and MZ?
>
> For the reference any text book about SM and/or fermi theory should have such formula.
> Since most of them are not free, I will simply recomend wikipedia which certainly have such information.
>
>>
But how does MadGraph use 80.419 GeV?
>
>
> You can tweak the value of each parameter in the param_card.
> As said the W mass is not a free parameter by default, so you need to change MZ or aEW or Gf to change his value.
>
>> I've also tried another process in which there is no W final state, and MadGraph gives different threshold.
>> a vm > vm e e
>> ( I've set the photon energy to be 0.25 GeV, so s = Eν * 1 GeV. )
>> The threshold is Eν = (0.511e-3+0.511e-3)^2 / 1 GeV = 1.04448e-6 GeV
>> But what MadGraph gives is the threshold to be ~2.3e-6 GeV.
>> Below is my script for this process. So anything wrong with what I did below?
>
>
>
> I will look at your mistake here on next wednesday. (but if you find it yourself in the meantime)
>
> Cheers,
>
> Olivier
>
>> On 14 Jun 2019, at 22:42, Qiji Xin <email address hidden> wrote:
>>
>> Question #681325 on MadGraph5_aMC@NLO changed:
>> https://answers.launchpad.net/mg5amcnlo/+question/681325
>>
>> Status: Answered => Open
>>
>> Qiji Xin is still having a problem:
>> Dear Olivier,
>>
>> Thanks so much for your answer!
>>
>> Could you please refer me to a literature/document about how the W mass is obtained from aeW, Gf and MZ?
>> And, for the (m_W+m_e)^2 ≃ 6462 GeV^2, I am using PDF value, which is experimentally measured, that m_W = 80.379 (± 0.012) GeV.
But how does MadGraph use 80.419 GeV?
>>
>> I've also tried another process in which there is no W final state, and MadGraph gives different threshold.
>> a vm > vm e e
>> ( I've set the photon energy to be 0.25 GeV, so s = Eν * 1 GeV. )
>> The threshold is Eν = (0.511e-3+0.511e-3)^2 / 1 GeV = 1.04448e-6 GeV
>> But what MadGraph gives is the threshold to be ~2.3e-6 GeV.
>> Below is my script for this process. So anything wrong with what I did below?
>>
>> Thanks a lot!
>> ======================================================================
>> import model sm-full #!
>>
>> generate a vm > vm e- e+ #!channel
>> display diagrams
>>
>> output
>> launch
>>
>> set lpp1 0
>> set lpp2 0
>> set ebeam1 0.25
>> set polbeam2 = -100
>> set no_parton_cut
>> set nevents 1000.
>> set ebeam2 1.0445e-6
>> launch
>> set ebeam2 1.3e-6
>> launch
>> set ebeam2 1.6e-6
>> launch
>> set ebeam2 2e-6
>> launch
>> set ebeam2 2.3e-6
>> launch
>> set ebeam2 2.7e-6
>>
>> --
>> You received this question notification because you are an answer
>> contact for MadGraph5_aMC@NLO.
>

Revision history for this message
Qiji Xin (xinqijisuper) said :
#8

Hi Olivier,

Thanks a lot for your hard work for this!

A quick follow-up question.
As you said, "The phase-space is so constrained that the code fails to find any phase-space points".
So does it help if I "set nevents" to a much larger number? I guess no (after some tries)?

Best,
Qiji

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) said :
#9

Hi,

Indeed that parameter does not have any impact on the situation.
If you want to increase the number of phase-point generated before the code gives up.
You need to edit the file
Source/dsample.f
and replace the line
      integer max_events
      parameter (max_events=50000000) !Maximum # events before get non_zero
by a bigger number.

I have tried by increasing it by a factor 100 and this was not solving the issue for the energy that I was trying.
I did not test it , but you should be able to see the effective threshold to decrease when doing that.
But for production, I would advise to run with
generate vm a > vm e e
which does not have such issue.

Cheers,

Olivier

> On 17 Jun 2019, at 03:47, Qiji Xin <email address hidden> wrote:
>
> Question #681325 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/681325
>
> Status: Answered => Open
>
> Qiji Xin is still having a problem:
> Hi Olivier,
>
> Thanks a lot for your hard work for this!
>
> A quick follow-up question.
> As you said, "The phase-space is so constrained that the code fails to find any phase-space points".
> So does it help if I "set nevents" to a much larger number? I guess no (after some tries)?
>
> Best,
> Qiji
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) said :
#11

Hi,

In this case the ordering should have no impact.
The ordering has only an impact when you have at least two propagators in t-channel (and therefore for at least 2>3 process).

The reason is linked to the fact that we probe the phase-space by generating the invariant mass of the sum of two particles. In your previous example, we have to pick the two electron for that in order to probe efficiently the phase-space close to the threshold and the pick is not "smart" but just based on the ordering.

Cheers,

Olivier

> On 26 Jun 2019, at 00:33, Qiji Xin <email address hidden> wrote:
>
> Question #681325 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/681325
>
> Status: Answered => Open
>
> Qiji Xin is still having a problem:
> Dear Olivier,
>
> Another follow-up question.
> When I am doing the fix-target DIS cross section calculation, like "νμ p → μ- j", do I also have to also put the neutrino as beam1 while p as beam2 in order to avoid the phase space issue?
> Does this issue happen only when there is a neutrino as the initial state?
>
> Thanks!
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Revision history for this message
Qiji Xin (xinqijisuper) said :
#12

Hi Olivier,

These answers solved me problem.
Thanks you so much for your help!