Heavy ion PDF in MG5_aMC_v2_6_7

PhotonFlux for heavy ion is not supported.
This will lead to a crash in the next version (I did not think about this special case)

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:
> https://answers.launchpad.net/mg5amcnlo/+question/686313
>
> 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
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi Olivier,

Thank you for your reply.

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

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:
> https://answers.launchpad.net/mg5amcnlo/+question/686313
>
> Status: Answered => Open
>
> Wenyi Song is still having a problem:
> Hi Olivier,
>
> Thank you for your reply.
>
> 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
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

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

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:
> https://answers.launchpad.net/mg5amcnlo/+question/686313
>
> Status: Answered => Open
>
> 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
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hi Olivier,

Thank you for your reply.

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

Hi,

please look at https://answers.launchpad.net/mg5amcnlo/+question/685493

Cheers,

Olivier

Hi Olivier,

Thank you for pointing me to this. I was able to follow the instruction and get it to work.

Now I would like to come up with a better parametrization for the photon PDF. I modified PhotonFlux.f but ended up with an error:

"DiscreteSampler:: Error, no point could be picked with random variable 0.7649418711662292 using upper bound found of NaN.
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO
STOP 1".

Part of the error message asks me to file a report to report this with a log file attached. I think it's better for me to first understand what this error means.

Thanks,
Wenyi

The real error is
"Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG IEEE_DIVIDE_BY_ZERO"
meaning that you have issue with your function (likely returning NAN or something as bad as that)

Cheers,

Olivier

Hi Olivier,

Thank you for your comment. I see what was causing this problem now.

I am running into a new error now:

"Survey return zero cross section.
   Typical reasons are the following:
   1) A massive s-channel particle has a width set to zero.
   2) The pdf are zero for at least one of the initial state particles
      or you are using maxjetflavor=4 for initial state b:s.
   3) The cuts are too strong."

This points me to the modified epa_proton function, which now provides an effective QED coupling (alpha times ion charge) for heavy ion collisions. The fraction of energy carried by the photon, x, is still with respect to one proton, instead of the entire ion. However, the photon flux parametrization I am using requires x to be with respect to the entire ion. I think this is why I am getting the zero cross section. What would you suggest for a redefinition of x?

Thanks for your continuous help!

-Wenyi

Here I do not know...

Cheers,

Olivier

Hi Olivier,

I see. Maybe I should look for a way around it. I truly appreciate all your time spent on this.

-Wenyi