Heavy ions in MG5 for MM production whithin photon fusion process

Hi,

The first issue that I have to comment is the meaning of lpp=2. Which makes the field "pdlabel" and "lhaid" to not be used.
So you need to decide to either run:
set lpp 2
or do
set lpp 1
set pdlabel lhapdf
set lhaid 901300

For lpp 2 you will get elastic photon, while for lpp1 you will get in principle inelastic photon but you need to check the paper of the according PDF set since different collaboration are handling differently photon PDF some collaboration include inelastic as well.

By default lpp 2 is set for proton beam and not for heavy ion beam. Starting with version 2.7.0, we started the support for heavy beam in that mode.
For that you need to specify in the run_card the mass of the ion and the number of proton.

#*********************************************************************
# Heavy ion PDF / rescaling of PDF *
#*********************************************************************
  1 = nb_proton1 # number of proton for the first beam
  0 = nb_neutron1 # number of neutron for the first beam
  -1.0 = mass_ion1 # mass of the heavy ion (first beam)
# Note that seting differently the two beams only work if you use
# group_subprocess=False when generating your matrix-element
  1 = nb_proton2 # number of proton for the second beam
  0 = nb_neutron2 # number of neutron for the second beam
  -1.0 = mass_ion2 # mass of the heavy ion (second beam)

So in short, I would not trust any of your result, even those who are not crashing.

For the crash per se, since the code is very old, the best is to test if you face the same issue in the latest version of the code, which is anyway needed in your case if you want to use the feature that you are looking for.

Cheers,

Olivier

Dear Olivier

First of all i want tot thank you for your help.
This time I used these commands as input for MadGraph5 v 2.8.2:

import model mono_spinzero/
generate a a > mm+ mm-
output simulari/monopoles_13tev_spinzeroPb
launch
analysis=madanalysis5+
set PbPb
set lpp1 1
set lpp2 1
set dynamical_scale_choice -1
set fixed_couplings False
set param_card mass 25 125
set decay 4110000 0.000000e+0
set gch 1 1.0
set pdlabel lhapdf
set lhaid 901300
set MMM 1000
set ebeam1 = 2750*208
set ebeam2 = 2750*208
set nevents = 10000

The output is this:
INFO: Update the dependent parameter of the param_card.dat
Generating 10000 events with run name run_01
survey run_01
INFO: compile directory
INFO: Using LHAPDF v6.3.0 interface for PDFs
compile Source Directory
Using random number seed offset = 21
INFO: Running Survey
Creating Jobs
Working on SubProcesses
INFO: P1_aa_mmpmmm
INFO: Idle: 1, Running: 0, Completed: 0 [ current time: 14h29 ]
INFO: Idle: 0, Running: 0, Completed: 1 [ 13.3s ]
INFO: Idle: 0, Running: 0, Completed: 1 [ 13.3s ]
INFO: End survey
refine 10000
Creating Jobs
INFO: Refine results to 10000
INFO: Generating 10000.0 unweighted events.
sum of cpu time of last step: 13 seconds
INFO: Effective Luminosity 1.2e+103 pb^-1
INFO: need to improve 0 channels
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.
   Please check/correct your param_card and/or your run_card.
Zero result detected: See https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/FAQ-General-14
INFO:
quit

I don't have much experience in MadGraph and I'm trying to learn it more and more. Thanks again for your help

Best regards,
Horea

Hi,

If you look at the info-file of your PDF set:
https://lhapdfsets.web.cern.ch/current/EPPS16nlo_CT14nlo_Pb208/EPPS16nlo_CT14nlo_Pb208.info
or to the associated paper:
https://arxiv.org/pdf/1612.05741.pdf
you will see that they are not photon PDF in that pdf set, so according to that PDf set, the probability to have a photon as part of the proton is zero. Consequently you are in the second case here, your pDF forbids such process to occur and therefore the cross-section is zero.

Also note that your command
set PbPb is doing the following:
INFO: modify parameter lpp1 of the run_card.dat to 1
INFO: modify parameter lpp2 of the run_card.dat to 1
INFO: modify parameter nb_proton1 of the run_card.dat to 82
INFO: modify parameter nb_neutron1 of the run_card.dat to 126
INFO: modify parameter mass_ion1 of the run_card.dat to 195.0820996698
INFO: modify parameter nb_proton2 of the run_card.dat to 82
INFO: modify parameter nb_neutron2 of the run_card.dat to 126
INFO: modify parameter mass_ion2 of the run_card.dat to 195.0820996698

When combine with a PDf set, we will do PDf rescaling (lpp=1) assuming a PDF from the proton.
So this is not designed to be used with the above PDf set where all the effect of many proton/neutron are handle internally to the PDF set.

Sorry that such part of the computation is not that user friendly. Lead collision is not our first focus since in general Argantur is more designed for such type of physics.

Cheers,

Olivier

PS: so I guess that the syntax that should fit your need is:
import model mono_spinzero
generate a a > mm+ mm-
output simulari/monopoles_13tev_spinzeroPb
launch
analysis=madanalysis5
set PbPb
set lpp1 2
set lpp2 2
set fixed_scale 2000
set fixed_couplings False
set param_card mass 25 125
set decay 4110000 0.000000e+0
set gch 1 1.0
set MMM 1000
set ebeam1 = 2750*208
set ebeam2 = 2750*208
set nevents = 10000

Hi Olivier.

Thanks for the help. I managed to get the cross sections with the pdf set LUXqed (set lhaid 82000) for Pb Pb. At this point I am trying to find out how I could change the impact parameter for heavy ion collisions into peripheral and ultra peripheral collisions. How can I do this in MadGraph5?

Best regrads,
Horea

You do not have any control on that in MG5aMC.

Cheers,

Olivier

> On 24 Mar 2021, at 17:10, Branzas Horea Florin <email address hidden> wrote:
>
> Question #696091 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/696091
>
> Branzas Horea Florin posted a new comment:
> Hi Olivier.
>
> Thanks for the help. I managed to get the cross sections with the pdf
> set LUXqed (set lhaid 82000) for Pb Pb. At this point I am trying to
> find out how I could change the impact parameter for heavy ion
> collisions into peripheral and ultra peripheral collisions. How can I do
> this in MadGraph5?
>
> Best regrads,
> Horea
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Dear Olivier

How could I get besides monopoles, the total number of particles produced in this type of interactions (PbPb)? I try to obtain total multiplicities that include monopoles. Many thanks for your help!!!

Best regards,
Horea

Hi,

I do not think that our tool can correctly handle such type of prediction.
The only tool that I know that can handle such type of prediction (as far as I know only in the SM) is argantur.
From the pythia8 team.

Cheers,

Olivier

> On 19 Jul 2021, at 11:20, Branzas Horea Florin <email address hidden> wrote:
>
> Question #696091 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/696091
>
> Branzas Horea Florin posted a new comment:
> Dear Olivier
>
> How could I get besides monopoles, the total number of particles
> produced in this type of interactions (PbPb)? I try to obtain total
> multiplicities that include monopoles. Many thanks for your help!!!
>
> Best regards,
> Horea
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Dear Olivier, first of all I want to wish you a Happy New Year!

After some searching, I want to modify the Photonflux.f file to get monopoles in ultrapheripheral (heavy ion PbPB)collisions but I can't do the right thing. I tried to modify the Phtonflux.f file by implementing the Weizsacker-Williams formula. The PhotonFlux.f file looks like that:

c/* ********************************************************* */
c/* Equivalent photon approximation structure function. * */
c/* V.M.Budnev et al., Phys.Rep. 15C (1975) 181 * */
c/* Improved Weizsaecker-Williams formula * */
c/* http://inspirehep.net/record/359425 * */
c/* ********************************************************* */
c provided by Tomasz Pierzchala - UCL

      real*8 function epa_lepton(x,q2max, mode)
      implicit none
      integer i, mode, imode
c mode is +3/-3 for electron and +4/-4 for muon
      real*8 x,phi_f
      real*8 xin(3:4)
      real*8 alpha
      real*8 f, q2min,q2max
      real*8 PI
      data PI/3.14159265358979323846/

      data xin/0.511d-3, 0.105658d0/ !electron mass in GeV

      alpha = .0072992701
      imode = abs(mode)

C // x = omega/E = (E-E')/E
      if (x.lt.1) then
         q2min= xin(imode)*xin(imode)*x*x/(1-x)
         if(q2min.lt.q2max) then
             f = alpha/2d0/PI*
     & (2d0*xin(imode)*xin(imode)*x*(-1/q2min+1/q2max)+
     & (2-2d0*x+x*x)/x*dlog(q2max/q2min))

         else
           f = 0.
         endif
      else
         f= 0.
      endif
c write (*,*) x,dsqrt(q2min),dsqrt(q2max),f
      if (f .lt. 0) f = 0
      epa_lepton = f

      end

      real*8 function epa_proton(x,q2max,beamid)
      integer i
      integer beamid
      real*8 x,phi_f
      real*8 xin
      real*8 alpha,qz
      real*8 f, qmi,qma, q2max
      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.938/ ! proton mass in GeV

      alpha = .0072992701
      qz = 0.71

      if (nb_proton(beamid).ne.1.or.nb_neutron(beamid).ne.0)then
         xin = mass_ion(beamid)
         alpha = alpha * nb_proton(beamid)
      endif

C // x = omega/E = (E-E')/E
      if (x.lt.1) then
         qmi= xin*xin*x*x/(1-x)
         if(qmi.lt.q2max) then
            f = alpha/PI*(phi_f(x,q2max/qz)-phi_f(x,qmi/qz))*(1-x)/x
         else
            f=0.
         endif
      else
         f= 0.
      endif
      if (f .lt. 0) f = 0
      epa_proton= f
      end

      real*8 function phi_f(x,qq)
      real*8 x, qq
      real*8 y,qq1,f
      real*8 a,b,c

       a = 7.16
       b = -3.96
       c = .028

      qq1=1+qq
      y= x*x/(1-x)
      f=(1+a*y)*(-log(qq1/qq)+1/qq1+1/(2*qq1*qq1)+1/(3*qq1*qq1*qq1))
      f=f + (1-b)*y/(4*qq*qq1*qq1*qq1);
      f=f+ c*(1+y/4)*(log((qq1-b)/qq1)+b/qq1+b*b/(2*qq1*qq1)+
     $b*b*b/(3*qq1*qq1*qq1))
      phi_f= f
      end
C ****Horea****
      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.938/ ! nucleus mass in GeV

      alpha = .0072992701

      if (nb_proton(beamid).ne.1.or.nb_neutron(beamid).ne.0) then
         bmin = 6.64
         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

where b is the impact parameter . When I run an example with the input data:

import model mono_spinzero/
generate a a > mm+ mm-
output simulari/monopoles_14tev_spinzeroGG
launch
analysis = madanalysis5
set lpp1 2
set lpp2 2
set nb_proton1 82
set nb_neutron1 126
set mass_ion1 = 195
set nb_proton2 82
set nb_neutron2 126
set mass_ion2 = 195
set fixed_scale 500
set fixed_fac_scale1 False
set fixed_fac_scale2 False
set dynamical_scale_choice -1
set fixed_couplings False
set param_card mass 25 125
set decay 4110000 0.000000e+0
set gch 1 1.0
set pdlabel lhapdf
set lhaid 82000
set MMM 10
set ebeam1 = 2560*208
set ebeam2 = 2560*208
set nevents = 100
I don't get any errors but if I change the impact parameter I get the same cross sections every time. Regarding the epa_proton function, how should it be modified in this case? I also mention that I use MadGraph5 version 3.2.0. Many thanks!!!

Best regards,
Horea

Hi,

Could you create a bzr branch such that I can use all the tools of launchpad to understand what you did change.

From what I see, you modified the function "get_ion_pdf("
which is not used for lpp=2 type of collision,
you need to modify directly the epa_proton function ( I guess)

Cheers,

Olivier

I'm sorry but I don't know exactly how to create and upload a branch with bzr branch ....

The original Photonflux.f file was:

      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

And then I changed to:

      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.938/ ! nucleus mass in GeV

      alpha = .0072992701

      if (nb_proton(beamid).ne.1.or.nb_neutron(beamid).ne.0) then
         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

The only thing I changed was the get_ion_pdf function and I don't know how to change the epa_proton function. Many thanks for your answer!

The only thing I changed was the get_ion_pdf function.

But this function is not used in this case.
In order for that function to be used, you need to set lpp=1 in your run.

 > I don't know how to change the epa_proton function.

i do not understand the physics that you try to achieve, so I can not comment on that,
but basically the change that you did in get_ion_pdf should be implemented in epa_proton if you want to use lpp=2

Cheers,

Olivier

Dear Olivier,

I'm trying to get magnetic monopoles through the process of photon fusion in the ultrapheripheral collisions of heavy ions(PbPb). More precisely, I am trying to obtain what is similar as in this article, chapter V. point D and chapter VI :

https://arxiv.org/pdf/2107.10789.pdf

I modified the epa_proton function and it works for lpp = 2, there are also some minor differences when I change the impact parameter. I suspect that lpp2 is a photon that comes from a proton and not from an given ion. For this reason I think the get_ion_pdf function should be modified and not epa_proton. If i use lpp1 (modified get_ion_pdf) i obtained:

The following floating-point exceptions are signalling: IEEE_INVALID_FLAG
....

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.
   Please check/correct your param_card and/or your run_card.

Many many thanks for your help!!!

Best regards,
Horea

Good that you succesfully hacked madgraph to do the physics that interest you.

Dear Olivier

After conducting some photon flux studies in the ultrapheripheral collisions of heavy ions, I learned that photons cannot have an energy greater than 80 GeV. How can I parameterize the fraction of energy that photons can have (x) from heavy ion collisions and and where can I make this change? setscale.f must be modified or this can be done in photonflux.f?. The idea is that at the moment I can get monopoles in ultra-peripheral collisions with mass much higher than 80GeV. When in fact it should have a maximum mass of 80 GeV (because photons cannot have an energy of more than 80 GeV in ultra-peripheral collisions of PbPb). I forgot to mention that energy of photon in heavy ion collision (X )is defined as :Lorentz factor/radius of heavy ion. Thank you very much !!!!

Hi,

I do not understand your question,
is madgraph has a cut that you do not want, or you want to impose one additional cut?
If this is the second this should be the same as changing the cutoff of the IWW formula no?
(so passing in fixed_scaled computation and changing that scale?)

Cheers,

Olivier

Dear Olivier

 I want to add a new cut to the fraction of energy carried by the photon (x). I am not sure why is so high at this moment. According to theory x can take a maximum value of 82 GeV in PbPb collision at c.m 5.5 TeV. It is as if the value of x is not in relation to an entire ion. In pp collision the energy of x can take much higher values than in PbPb collisions an here i think i must to add a cut somehow.

>If this is the second this should be the same as changing the cutoff of the IWW formula no?

Yeah something like that

Thanks for your time!!!

So yes changing to fix scale and set the scale to 80 would likely be the solution

> On 13 Apr 2022, at 13:50, Branzas Horea Florin <email address hidden> wrote:
>
> Question #696091 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/696091
>
> Status: Answered => Open
>
> Branzas Horea Florin is still having a problem:
> Dear Olivier
>
> I want to add a new cut to the fraction of energy carried by the photon
> (x). I am not sure why is so high at this moment. According to theory x
> can take a maximum value of 82 GeV in PbPb collision at c.m 5.5 TeV. It
> is as if the value of x is not in relation to an entire ion. In pp
> collision the energy of x can take much higher values than in PbPb
> collisions an here i think i must to add a cut somehow.
>
>> If this is the second this should be the same as changing the cutoff of
> the IWW formula no?
>
> Yeah something like that
>
> Thanks for your time!!!
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Is this command the right one?

set fixed_scale 80

or

set fixed_fac_scale True
set scale 80

The first one should be correct.
You should see that such command will change multiple line in the run_card including
fixed_fac_scale (which will be set to True)

> set scale 80

I do not think that we have a parameter scale (but I did not double check)

In any case, please check the log, the impact of such command will be clearly printed

Cheers,

Olivier

> On 13 Apr 2022, at 17:25, Branzas Horea Florin <email address hidden> wrote:
>
> Question #696091 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/696091
>
> Status: Answered => Open
>
> Branzas Horea Florin is still having a problem:
> Is this command the right one?
>
> set fixed_scale 80
>
> or
>
> set fixed_fac_scale True
> set scale 80
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

set fixed_scale 80
INFO: modify parameter fixed_fac_scale of the run_card.dat to True
INFO: modify parameter fixed_ren_scale of the run_card.dat to True
INFO: modify parameter scale of the run_card.dat to 80
INFO: modify parameter dsqrt_q2fact1 of the run_card.dat to 80
INFO: modify parameter dsqrt_q2fact2 of the run_card.dat to 80

But it does not seem to affect the cross sections obtained in any way even if I set a value lower than 80 or higher.
The idea is that I shouldn't get magnetic monopoles in PbPb (trought photon fusion channel a a > mm+ mm-) collisions with a mass greater than 80 GeV because photons from PbPb collisions cannot have an energy greater than 80 GeV. But in my case I can get monopoles with a mass greater than 1 TeV. That's why I think that the fraction of energy carried by the photon, x, needs a certain cut, or a new parameterization.

https://inspirehep.net/files/14e301bc3365f02dd078451a66eb64cc
At Chapter 1 and Chapter 2 there is a better explication regarding this problem but in my case are Monopoles and not Higgs particles.
Thank you very much for your guidance

Hi,

From your paper
1) it is not a clear maximum
2) it is not that clear in which frame that statement hold (this is obviously a non lorentz invariant statement ...)

Cheers,

Olivier

Dear Olivier

I think I solved it with the implementation of the flux of photons from the ultrapheripheral collisions of heavy ions. It seems that the ion radius and the impact parameter play a very important role in the energy of photons coming from the magnetic fields of the ions. But I'm facing a new problem now. In the Geant4 simulator when I want to see how monopoles interact (i use photon fusion channel aa > mm + mm-) with the geometry of the detector everything is ok when I use pp collisions (default photonflux.f, Luxqed as PDF, lpp1 and lpp2 set to 1) . But when I use the data obtained for heavy ions (lpp1,2 set to 2 for elastic photons from ions and with modified photonflux.f) the simulation in Geant4 stops immediately because I receive certain warnings for unknown particles ( with code 1201 and 1203 ) that travel a finite distance. I checked the structure of the lhe file I got for heavy ions, everything seems fine. Also the cross sections obtained are as expected as values. Have you encountered this type of error before? Thank you very much!!!

Best regards,
Horea

Hi,

I have never run geant4 so can not really help here.
I will assign this thread to Hua-Sheng since he created a branch UPS (that I have to review actually)
and that might go in the same direction of what you did. So he might be of more help than me in such kind of physics.
(and know why pythia is adding such particles)

Cheers,

Olivier

Hi, yes, as Olivier said, we will have a UPC extension of photon-photon collisions in proton-proton, proton-nuclear and nuclear-nuclear collisions in MG5, which implements two different form factors and survival probability. However, regarding your last question about Geant4 (or even Pythia8), I have no clue how they interpret the UPC lhe events. Cheers, Hua-Sheng

Hello

Thank you very much for the answers !!!

Best regards,
Horea