Heavy ion init incorrectly reported in LHE file

[Zachary Marshall]

Dear authors,

We've seen what is apparently an issue in the LHE output of MadGraph when trying to run with a heavy ion initial state (as described in https://arxiv.org/pdf/2207.03012.pdf). The process is photon-induced:

generate a a > y2 > a a

and the card is set accordingly to use ions:

2 = lpp1 ! beam 1 type
2 = lpp2 ! beam 2 type
557440.0 = ebeam1 ! beam 1 total energy in GeV
557440.0 = ebeam2 ! beam 2 total energy in GeV
edff = pdlabel ! PDF set
230000 = lhaid ! if pdlabel=lhapdf, this is the lhapdf number
82 = nb_proton1
126 = nb_neutron1
82 = nb_proton2
126 = nb_neutron2

The LHE init line looks like this:

<init>
2212 2212 6.500000e+03 6.500000e+03 0 0 247000 247000 -4 1
   +1.5507000e-07 +4.7623000e-10 +1.5507000e-07 1
<generator name='MadGraph5_aMC@NLO' version='3.5.1'>please cite 1405.0301 </gen\
erator>
</init>

which seems wrong — we're back to 13 TeV center of mass proton-proton collisions somehow. Can you help us understand what might be going wrong?

Thanks,
Zach

[Olivier Mattelaer (olivier-mattelaer)]

Can you check 3.5.3?
Are you sure that this is not an Athena post-processing issue (if you use that framework?)

I do have this in the file:

   399 <init>
    400 2212 2212 5.574400e+05 5.574400e+05 0 0 0 0 -4 1
    401 1.461900e+05 3.975100e+02 1.461900e+05 1
    402 <generator name='MadGraph5_aMC@NLO' version='3.5.3'>please cite 1405.0301 </generator>
    403 </init>

for the following script:

generate a a > j j
output
launch
set lpp 2
set ebeam 557440.0
set pdlabel edff
set nb_proton1 82
set nb_neutron1 126
set nb_proton2 82
set nb_neutron2 126

Cheers,

Olivier

[Zachary Marshall (zach-marshall)]

Thanks Olivier! I do see when running stand-alone (both in 3.5.1 and 3.5.3) that the center of mass energy is respected, which suggests there's something fishy in our launch configuration. With that said, for sure even in your example the particles are wrong (the beams aren't protons), right?

Thanks,
Zach

[Olivier Mattelaer (olivier-mattelaer)]

Indeed I never tried to set up those number following the 10 digit code that can be computed for heavy ion state. I can fix it but this will increase the number of digit in that section and some code might not parse it (or refuse to shower it /...).

I will look at a patch for that.

Olivier

> On 16 Jan 2024, at 19:54, Zachary Marshall <email address hidden> wrote:
>
> Thanks Olivier! I do see when running stand-alone (both in 3.5.1 and
> 3.5.3) that the center of mass energy is respected, which suggests
> there's something fishy in our launch configuration. With that said, for
> sure even in your example the particles are wrong (the beams aren't
> protons), right?
>
> Thanks,
> Zach
>
> --
> You received this bug notification because you are subscribed to
> MadGraph5_aMC@NLO.
> https://bugs.launchpad.net/bugs/2049278
>
> Title:
> Heavy ion init incorrectly reported in LHE file
>
> Status in MadGraph5_aMC@NLO:
> New
>
> Bug description:
> Dear authors,
>
> We've seen what is apparently an issue in the LHE output of MadGraph
> when trying to run with a heavy ion initial state (as described in
> https://arxiv.org/pdf/2207.03012.pdf). The process is photon-induced:
>
> generate a a > y2 > a a
>
> and the card is set accordingly to use ions:
>
> 2 = lpp1 ! beam 1 type
> 2 = lpp2 ! beam 2 type
> 557440.0 = ebeam1 ! beam 1 total energy in GeV
> 557440.0 = ebeam2 ! beam 2 total energy in GeV
> edff = pdlabel ! PDF set
> 230000 = lhaid ! if pdlabel=lhapdf, this is the lhapdf number
> 82 = nb_proton1
> 126 = nb_neutron1
> 82 = nb_proton2
> 126 = nb_neutron2
>
> The LHE init line looks like this:
>
> <init>
> 2212 2212 6.500000e+03 6.500000e+03 0 0 247000 247000 -4 1
> +1.5507000e-07 +4.7623000e-10 +1.5507000e-07 1
> <generator name='MadGraph5_aMC@NLO' version='3.5.1'>please cite 1405.0301 </gen\
> erator>
> </init>
>
> which seems wrong — we're back to 13 TeV center of mass proton-proton
> collisions somehow. Can you help us understand what might be going
> wrong?
>
> Thanks,
> Zach
>
> To manage notifications about this bug go to:
> https://bugs.launchpad.net/mg5amcnlo/+bug/2049278/+subscriptions
>

[Zachary Marshall (zach-marshall)]

Just a brief follow-up: yes, it looks like the more complex running here caused the center of mass change:

set group_subprocesses Auto
set ignore_six_quark_processes False
set low_mem_multicore_nlo_generation False
set complex_mass_scheme False
set include_lepton_initiated_processes False
set gauge unitary
set loop_optimized_output True
set loop_color_flows False
set max_npoint_for_channel 0
set default_unset_couplings 99
set max_t_for_channel 99
set zerowidth_tchannel True
set nlo_mixed_expansion True
import model DMspin2-full
define p = g u c d s u~ c~ d~ s~
define j = g u c d s u~ c~ d~ s~
define l+ = e+ mu+
define l- = e- mu-
define vl = ve vm vt
define vl~ = ve~ vm~ vt~
define p = 21 2 4 1 3 -2 -4 -1 -3 5 -5
define j = p
generate a a > y2 > a a
output aay2aa -nojpeg
launch aay2aa
update ion_pdf
launch aay2aa
update ion_pdf

and at the end the center of mass has been reset to 6500 GeV.

Cheers,
Zach

[Zachary Marshall (zach-marshall)]

Thanks Olivier! I think it's probably worth an email to the Pythia8 folks to see what they're going to try to do. I don't know how they're handling these things.

Best,
Zach

[Olivier Mattelaer (olivier-mattelaer)]

ok keep me updated. If needed I have a patch available:

diff --git a/Template/LO/Source/setrun.f b/Template/LO/Source/setrun.f
index ba99f895a..876d325ed 100644
--- a/Template/LO/Source/setrun.f
+++ b/Template/LO/Source/setrun.f
@@ -145,9 +145,21 @@ c banner.py is expected to correct such wrong run_card.
 C Fill common block for Les Houches init info
       do i=1,2
         if(lpp(i).eq.1.or.lpp(i).eq.2) then
- idbmup(i)=2212
+ if (nb_proton(i).eq.1.and.nb_neutron(i).eq.0) then
+ idbmup(i)=2212
+ elseif (nb_proton(i).eq.0.and.nb_neutron(i).eq.1) then
+ idbmup(i)=2112
+ else
+ idbmup(i) = 1000000000 + (nb_proton(i)+nb_neutron(i))*10
+ $ + nb_proton(i)*10000
+ endif
         elseif(lpp(i).eq.-1.or.lpp(i).eq.-2) then
- idbmup(i)=-2212
+ if (nb_proton(i).eq.1.and.nb_neutron(i).eq.0) then
+ idbmup(i)=-2212
+ else
+ idbmup(i) = -1*(1000000000 + (nb_proton(i)+nb_neutron(i))*10
+ $ + nb_proton(i)*10000)
+ endif
         elseif(lpp(i).eq.3) then
           idbmup(i)=11
         elseif(lpp(i).eq.-3) then
diff --git a/madgraph/various/banner.py b/madgraph/various/banner.py
index ece8d0a33..bee9a0061 100755
--- a/madgraph/various/banner.py
+++ b/madgraph/various/banner.py
@@ -2766,11 +2766,22 @@ class RunCard(ConfigFile):
                     fsock.writelines(line)
             fsock.close()

- @staticmethod
- def get_idbmup(lpp):
+ def get_idbmup(self, lpp, beam=1):
         """return the particle colliding pdg code"""
         if lpp in (1,2, -1,-2):
- return math.copysign(2212, lpp)
+ target = 2212
+ if 'nb_proton1' in self:
+ nbp = self['nb_proton%s' % beam]
+ nbn = self['nb_neutron%s' % beam]
+ if nbp == 1 and nbn ==0:
+ target = 2212
+ elif nbp==0 and nbn ==1:
+ target = 2112
+ else:
+ target = 1000000000
+ target += 10 * (nbp+nbn)
+ target += 10000 * nbp
+ return math.copysign(target, lpp)
         elif lpp in (3,-3):
             return math.copysign(11, lpp)
         elif lpp in (4,-4):
@@ -2786,8 +2797,8 @@ class RunCard(ConfigFile):
         the first line of the <init> block of the lhe file."""

         output = {}
- output["idbmup1"] = self.get_idbmup(self['lpp1'])
- output["idbmup2"] = self.get_idbmup(self['lpp2'])
+ output["idbmup1"] = self.get_idbmup(self['lpp1'], beam=1)
+ output["idbmup2"] = self.get_idbmup(self['lpp2'], beam=2)
         output["ebmup1"] = self["ebeam1"]
         output["ebmup2"] = self["ebeam2"]
         output["pdfgup1"] = 0

[Olivier Mattelaer (olivier-mattelaer)]

Hi Zach,

I was planing to release 3.5.4 today, and therefore checking any open issue and face this one?
Did you contact them? Should I push the patch in 3.5.4 to include the 10 digit information from the beam?

Now I also realized that we have not fix the issue with the energy.
In your script you use the command
update ion_pdf
Which only add a block in the run_card.dat (with default value for p p).
Like this such command is pointless since you do not specify any of the content (at least from the script provided).
So indeed from the script that you provide, I would expect p p @13TeV

[Zachary Marshall (zach-marshall)]

I don't know that we've heard back from the Pythia8 authors; I've just pinged the contacts on our side to check. I think it would be better to go ahead without the patch for now, just in case this breaks some assumption they're making.

I believe removing the extra junk in the launch command helped, but I'm not 100% sure if we're the whole way there.

Thanks for following up!

Best,
Zach