# Heavy ion PDF in MG5_aMC_v2_6_7

Asked by Wenyi Song on 2019-11-27

Hello,

I would like to understand the implementation of heavy-ion photon flux in MG5_aMC_v2_6_7, as I am trying to benchmark the cross section of an elastic photon-fusion process, for example, a a > b b~, with a constrained invariant mass for the final state.

In Template/LO/Source/PDF/PhotonFlux.f, the heavy-ion photon flux is given by the PDF up to a multiplication of the total number of nucleons.

double precision function get_ion_pdf(pdf, pdg, nb_proton, nb_neutron)
C***********************************************************************
C computing (heavy) ion contribution from proton PDF
C***********************************************************************
double precision pdf(-7:7)
double precision tmppdf(-2:2)
integer pdg
integer nb_proton
integer nb_neutron
double precision tmp1, tmp2

if (nb_proton.eq.1.and.nb_neutron.eq.0)then
get_ion_pdf = pdf(pdg)
return
endif

if (pdg.eq.1.or.pdg.eq.2) then
tmp1 = pdf(1)
tmp2 = pdf(2)
tmppdf(1) = nb_proton * tmp1 + nb_neutron * tmp2
tmppdf(2) = nb_proton * tmp2 + nb_neutron * tmp1
get_ion_pdf = tmppdf(pdg)
else if (pdg.eq.-1.or.pdg.eq.-2) then
tmp1 = pdf(-1)
tmp2 = pdf(-2)
tmppdf(-1) = nb_proton * tmp1 + nb_neutron * tmp2
tmppdf(-2) = nb_proton * tmp2 + nb_neutron * tmp1
get_ion_pdf = tmppdf(pdg)
else
get_ion_pdf = pdf(pdg)*(nb_proton+nb_neutron)
endif

C set correct PDF normalisation
get_ion_pdf = get_ion_pdf * (nb_proton+nb_neutron)
return
end

Hence the same process will have a higher cross section in heavy ion collisions than in pp collisions. However, when I launched Pb collisions with the following updates to the run card, MadGraph returned zero cross section.

set mmbb 100.
set mmbbmax 140.
set lpp1 2
set lpp2 2
set ebeam1 7000
set ebeam2 7000
set nb_proton1 82
set nb_neutron1 126
set nb_proton2 82
set nb_neutron2 126

In order to understand the reason, I kept everything the same but the number of nucleons, for which I tried a series of numbers between 2 and 50. I noticed that the cross section decreases as a function of the total nucleon number and vanishes at 50. This is clearly not what the corresponding photon flux is supposed to do, but I have no clue why this is the case.

Any help would be appreciated.

Cheers,
Wenyi

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
2019-12-04
2019-12-04
 Olivier Mattelaer (olivier-mattelaer) said on 2019-11-27: #1

PhotonFlux for heavy ion is not supported.

This being said ebeam1 is the energy of the beam and not the energy of the proton/neutron.
So I do not think that your comparison is meaningful.

Cheers,

Olivier

> On 27 Nov 2019, at 18:37, Wenyi Song <email address hidden> wrote:
>
> New question #686313 on MadGraph5_aMC@NLO:
>
> Hello,
>
> I would like to understand the implementation of heavy-ion photon flux in MG5_aMC_v2_6_7, as I am trying to benchmark the cross section of an elastic photon-fusion process, for example, a a > b b~, with a constrained invariant mass for the final state.
>
> In Template/LO/Source/PDF/PhotonFlux.f, the heavy-ion photon flux is given by the PDF up to a multiplication of the total number of nucleons.
>
> double precision function get_ion_pdf(pdf, pdg, nb_proton, nb_neutron)
> C***********************************************************************
> C computing (heavy) ion contribution from proton PDF
> C***********************************************************************
> double precision pdf(-7:7)
> double precision tmppdf(-2:2)
> integer pdg
> integer nb_proton
> integer nb_neutron
> double precision tmp1, tmp2
>
> if (nb_proton.eq.1.and.nb_neutron.eq.0)then
> get_ion_pdf = pdf(pdg)
> return
> endif
>
> if (pdg.eq.1.or.pdg.eq.2) then
> tmp1 = pdf(1)
> tmp2 = pdf(2)
> tmppdf(1) = nb_proton * tmp1 + nb_neutron * tmp2
> tmppdf(2) = nb_proton * tmp2 + nb_neutron * tmp1
> get_ion_pdf = tmppdf(pdg)
> else if (pdg.eq.-1.or.pdg.eq.-2) then
> tmp1 = pdf(-1)
> tmp2 = pdf(-2)
> tmppdf(-1) = nb_proton * tmp1 + nb_neutron * tmp2
> tmppdf(-2) = nb_proton * tmp2 + nb_neutron * tmp1
> get_ion_pdf = tmppdf(pdg)
> else
> get_ion_pdf = pdf(pdg)*(nb_proton+nb_neutron)
> endif
>
> C set correct PDF normalisation
> get_ion_pdf = get_ion_pdf * (nb_proton+nb_neutron)
> return
> end
>
> Hence the same process will have a higher cross section in heavy ion collisions than in pp collisions. However, when I launched Pb collisions with the following updates to the run card, MadGraph returned zero cross section.
>
> set mmbb 100.
> set mmbbmax 140.
> set lpp1 2
> set lpp2 2
> set ebeam1 7000
> set ebeam2 7000
> set nb_proton1 82
> set nb_neutron1 126
> set nb_proton2 82
> set nb_neutron2 126
>
> In order to understand the reason, I kept everything the same but the number of nucleons, for which I tried a series of numbers between 2 and 50. I noticed that the cross section decreases as a function of the total nucleon number and vanishes at 50. This is clearly not what the corresponding photon flux is supposed to do, but I have no clue why this is the case.
>
> Any help would be appreciated.
>
> Cheers,
> Wenyi
>
> --

 Wenyi Song (wenyisong) said on 2019-11-27: #2

Hi Olivier,

If this version only supports photon fluxes for electrons and protons, one would have to scale proton's photon flux to make it work for heavy ion. Am I correct?

Regarding your point about the meaningfulness of the test, I just wanted to figure out what has the zero cross section. I thought if I observed a change in the cross section, then the cause would have to be PDF-related, thereby rule out the other possibilities, i.e., phase space integration, narrow decay width...

Assuming heavy ion collision is supported, I need to be careful with the ebeam setting. If ebeam is the energy of the beam, i.e., 7 TeV as in this case, then I should just leave it as 7000 as I did above?

Many thanks,
Wenyi

 Olivier Mattelaer (olivier-mattelaer) said on 2019-11-27: #3

Hi Wenyi,

> If this version only supports photon fluxes for electrons and protons,
> one would have to scale proton's photon flux to make it work for heavy
> ion. Am I correct?

Well, the first question is what do expect for this photon flux.
From the syntax that you tried, I would expect that you want the coherent emission
of the photon from the heavy ion and not the sum of the coherent emission of the proton.

Assuming no interference, the second can be indeed obtained by rescaling from the one of a single proton. For the first a patch technically exists but I do not think that it works (issue of phase-space integration)

> Assuming heavy ion collision is supported, I need to be careful with the
> ebeam setting. If ebeam is the energy of the beam, i.e., 7 TeV as in
> this case, then I should just leave it as 7000 as I did above?

If you keep 7000 per beam but increase the number of proton/neutron
then the energy of each proton/neutron will decrease (obviously).
For scaling test, we typically keep the energy per proton constant.

Cheers,

Olivier

> On 27 Nov 2019, at 21:27, Wenyi Song <email address hidden> wrote:
>
> Question #686313 on MadGraph5_aMC@NLO changed:
>
>
> Wenyi Song is still having a problem:
> Hi Olivier,
>
>
> If this version only supports photon fluxes for electrons and protons,
> one would have to scale proton's photon flux to make it work for heavy
> ion. Am I correct?
>
> Regarding your point about the meaningfulness of the test, I just wanted
> to figure out what has the zero cross section. I thought if I observed a
> change in the cross section, then the cause would have to be PDF-
> related, thereby rule out the other possibilities, i.e., phase space
> integration, narrow decay width...
>
> Assuming heavy ion collision is supported, I need to be careful with the
> ebeam setting. If ebeam is the energy of the beam, i.e., 7 TeV as in
> this case, then I should just leave it as 7000 as I did above?
>
> Many thanks,
> Wenyi
>
> --

 Wenyi Song (wenyisong) said on 2019-12-03: #4

HI Olivier,

After some literature search, I implemented an analytical expression for the heavy ion pdf into PhotonFlux.f, but I am still getting zero cross sections. I don't think this is caused by the implemented pdf being evaluated to zero, because I also passed non-zero constants to it and ended up with the same cross section. Do you know why this is happening? Here is how I defined the pdf.

real*8 function get_ion_pdf(x,xi,beamid)
integer i
integer beamid
real*8 x
real*8 xin,bmin,xi
real*8 alpha
real*8 f,y1,y2
real*8 PI
integer nb_proton(2), nb_neutron(2)
common/to_heavyion_pdg/ nb_proton, nb_neutron
double precision mass_ion(2)
common/to_heavyion_mass/mass_ion

data PI/3.14159265358979323846/

data xin/0.9315/ ! nucleus mass in GeV

alpha = .0072992701

! Ion option turned on by !1 proton number or non-zero neutron number
if (nb_proton(beamid).ne.1.or.nb_neutron(beamid).ne.0) then
! impact parameter in GeV^-1 (1fm ~ 5 GeV^-1)
bmin = 5.*1.2*mass_ion(beamid)**(1/3)
if (x.lt.1) then
xi= x*xin*bmin
y1=2*(-log(xi/2)+0.5772)
y2=xi*xi*(1/xi**2-(-log(xi/2)+0.5772)**2)
f = alpha/PI*nb_proton(beamid)**2/x*(y1-y2)
else
f= 0.
endif
if (f .lt. 0) f = 0
get_ion_pdf= f
endif
end

Many thanks,
Wenyi

 Olivier Mattelaer (olivier-mattelaer) said on 2019-12-03: #5

Hi,

My guess is that this is related to the phase-space integration issue.
Photon-PDF are quite soft and are the phase-space integrator does not know this.
The phase-space integrator is actually already quite border line for the proton case.
If you still decrease the value of allowed x, this can indeed make the code to badly behaves.

Cheers,

Olivier

> On 3 Dec 2019, at 19:58, Wenyi Song <email address hidden> wrote:
>
> Question #686313 on MadGraph5_aMC@NLO changed:
>
>
> Wenyi Song is still having a problem:
> HI Olivier,
>
> After some literature search, I implemented an analytical expression for
> the heavy ion pdf into PhotonFlux.f, but I am still getting zero cross
> sections. I don't think this is caused by the implemented pdf being
> evaluated to zero, because I also passed non-zero constants to it and
> ended up with the same cross section. Do you know why this is happening?
> Here is how I defined the pdf.
>
> real*8 function get_ion_pdf(x,xi,beamid)
> integer i
> integer beamid
> real*8 x
> real*8 xin,bmin,xi
> real*8 alpha
> real*8 f,y1,y2
> real*8 PI
> integer nb_proton(2), nb_neutron(2)
> common/to_heavyion_pdg/ nb_proton, nb_neutron
> double precision mass_ion(2)
> common/to_heavyion_mass/mass_ion
>
> data PI/3.14159265358979323846/
>
> data xin/0.9315/ ! nucleus mass in GeV
>
> alpha = .0072992701
>
> ! Ion option turned on by !1 proton number or non-zero neutron number
> if (nb_proton(beamid).ne.1.or.nb_neutron(beamid).ne.0) then
> ! impact parameter in GeV^-1 (1fm ~ 5 GeV^-1)
> bmin = 5.*1.2*mass_ion(beamid)**(1/3)
> if (x.lt.1) then
> xi= x*xin*bmin
> y1=2*(-log(xi/2)+0.5772)
> y2=xi*xi*(1/xi**2-(-log(xi/2)+0.5772)**2)
> f = alpha/PI*nb_proton(beamid)**2/x*(y1-y2)
> else
> f= 0.
> endif
> if (f .lt. 0) f = 0
> get_ion_pdf= f
> endif
> end
>
> Many thanks,
> Wenyi
>
> --

 Wenyi Song (wenyisong) said on 2019-12-04: #6

Hi Olivier,

I don't quite understand the issue here. I think I understand that the soft photon PDF makes the code badly behave, but when I replaced the PDF with a constant, i.e., 1, I still ended up with the same zero cross sections. In addition, what would you suggest to resolve this issue?

Cheers,
Wenyi

 Olivier Mattelaer (olivier-mattelaer) said on 2019-12-04: #7