# Cross section calculation abnormal when the dark photon mass is extreme small

Asked by Yao Yao on 2020-09-24

Dear expert,

I am working on a dark photon analysis with signal: p p > m+ m- zp using the CMS detector. I am currently generating the signals with HAHM (Hidden Abelian Higgs Model). So far I've been analyzing signal samples with dark photon (zp) mass equals to 0.5 MeV, and things worked fine at that point. Now we are trying lower down the mass of dark photon, to see how far our analysis can go.

The problem of cross section calculation appears when I am trying to lower down the mass of dark photon to 0.5 KeV and even further. We realize that the total cross section for the processes becomes dramatically huge and some diagrams dominate the process in a scale that I feel quite unreasonable. (Of course please feel free to debate if you think the calculation result is reasonable.)

The cuts that I specified for all my signal samples:
1) pT of dark photon > 10 GeV,
2) dilepton invariant mass (l+ l- same flavor) > 10 GeV,
3) deltaR (dark photon, muon) > 0.4

The total cross section of the signal with:
1) dark photon mass = 0.5 MeV, xsec = 0.16 pb, dominated by P1_qq_mpmmzp (0.1459 pb), diagram #2 and #12 dominates.
2) dark photon mass = 0.5 KeV, xsec = 1.3e+31 pb, dominated by P1_bbx_mpmmzp (1.049e+31 pb), and P1_ccx_mpmmzp (2.194e+30 pb), while P1_qq_mpmmzp is small compared to them but still quite huge (7.233e+07 pb), and in the P1_qq_mpmmzp, diagram #9 and #12 dominates.
3) dark photon mass = 5e-6 eV, xec = 4.7e+79 pb, dominated by P1_bbx_mpmmzp (3.905e+79 pb), and P1_ccx_mpmmzp (8.22e+78), and P1_qq_mpmmzp (6.262e+52 pb), and in the P1_qq_mpmmzp, diagram #9 and #12 dominates.

Diagrams #9 and #12 are t (or u) channels that one branch is zp, the other branch is zp > m+ m-.
Diagram #2 is s channel that q q > z > m+ m-, and zp recoil from m+.

I cannot figure out a reason why the cross section changes in such a big scale by only changing the mass of the dark photon, and why certain diagrams suddenly dominate and blow up. We are considering that dark photon should be very similar to a SM photon, but with a bit mass. Please let us know if you have any thinking on this or anything we can check.

Thank you,
Yao

## Question information

Language:
English Edit question
Status:
For:
Assignee:
No assignee Edit question
Last query:
2020-10-01
2020-10-02
 Olivier Mattelaer (olivier-mattelaer) said on 2020-09-24: #1

Hi,

I have a few questions:
1) How did you implement cut #3?
2) What is the mass of the muon?
3) Do you have a width to your dark photon?
4) How many events are generated for those runs? What is the statistical error of such runs?

Cheers,

Olivier

 Yao Yao (yyao1) said on 2020-09-24: #2

Hi Olivier,

1) for cut #3, I asked this question to you before in [Question #686216]: and solved it by writing some code in dummy_fct.f. (Please see the quote)
2) The mass of the muon should be standard, about 105 MeV. It should be set in HAHM model, and I didn't touch it.
3) I don’t set a width for dark photon, or maybe I should?
4) I generated 50000 events for each of them. The statistical error of the cross section for each is :

0.1601 +/- 0.04335 pb
1.268e+31 +/- 1.792e+28 pb
4.727e+79 +/- 5.962e+86 pb

Yao

> Hi Olivier,
>
> I got Madgraph running fast now when generating the events! I write the following code in dummy_fct.f. It implies that a dark photon needs to be deltaR>0.4 from any muon. Please help me check if the code seems good?
>
>
> include 'maxamps.inc'
> integer iproc,i,j
> integer idup(nexternal,maxproc,maxsproc)
> integer mothup(2,nexternal)
> integer icolup(2,nexternal,maxflow,maxsproc)
> include 'leshouche.inc'
>
> double precision p1(0:3),p2(0:3)
> double precision R2
> double precision deltaR
> double precision deltaRmin
> deltaR=1d3
> deltaRmin=1d3
>
> do iproc=1,maxsproc
>
> do i=nincoming+1,nexternal-1
> if (idup(i,1,iproc).eq.13.or.idup(i,1,iproc).eq.-13)then
> do j=i+1,nexternal
> if (idup(j,1,iproc).eq.1023)then
> p1(0:3)=P(0:3,i)
> p2(0:3)=P(0:3,j)
> deltaR=dsqrt(R2(p1,p2))
> if (deltaR.lt.deltaRmin)then
> deltaRmin=deltaR
> endif
> endif
> enddo
> endif
> enddo
> do i=nincoming+1,nexternal-1
> if(idup(i,1,iproc).eq.1023)then
> do j=i+1,nexternal
> if(idup(j,1,iproc).eq.13.or.idup(j,1,iproc).eq.-13)then
> p1(0:3)=P(0:3,i)
> p2(0:3)=P(0:3,j)
> deltaR=dsqrt(R2(p1,p2))
> if (deltaR.lt.deltaRmin)then
> deltaRmin=deltaR
> endif
> endif
> enddo
> endif
> enddo
> enddo
> dummy_cuts=.true.
>
> if(deltaRmin.lt.4.0d-1) then
> dummy_cuts=.false.
> endif
> return
> end

> On Sep 24, 2020, at 4:25 PM, Olivier Mattelaer <email address hidden> wrote:
>
> Your question #693074 on MadGraph5_aMC@NLO changed:
>
> Status: Open => Answered
>
> Olivier Mattelaer proposed the following answer:
> Hi,
>
> I have a few questions:
> 1) How did you implement cut #3?
> 2) What is the mass of the muon?
> 3) Do you have a width to your dark photon?
> 4) How many events are generated for those runs? What is the statistical error of such runs?
>
> Cheers,
>
> Olivier
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
>
> You received this question notification because you asked the question.

 Yao Yao (yyao1) said on 2020-09-24: #3

Oops, sorry a type on the last cross section, it should be e+76 for the error.
> On Sep 24, 2020, at 4:40 PM, Yao Yao <email address hidden> wrote:
>
> Hi Olivier,
>
> 1) for cut #3, I asked this question to you before in [Question #686216]: and solved it by writing some code in dummy_fct.f. (Please see the quote)
> 2) The mass of the muon should be standard, about 105 MeV. It should be set in HAHM model, and I didn't touch it.
> 3) I don’t set a width for dark photon, or maybe I should?
> 4) I generated 50000 events for each of them. The statistical error of the cross section for each is :
>
> 0.1601 +/- 0.04335 pb
> 1.268e+31 +/- 1.792e+28 pb
> 4.727e+79 +/- 5.962e+86 pb
>
> Yao
>
>> Hi Olivier,
>>
>> I got Madgraph running fast now when generating the events! I write the following code in dummy_fct.f. It implies that a dark photon needs to be deltaR>0.4 from any muon. Please help me check if the code seems good?
>>
>>
>> include 'maxamps.inc'
>> integer iproc,i,j
>> integer idup(nexternal,maxproc,maxsproc)
>> integer mothup(2,nexternal)
>> integer icolup(2,nexternal,maxflow,maxsproc)
>> include 'leshouche.inc'
>>
>> double precision p1(0:3),p2(0:3)
>> double precision R2
>> double precision deltaR
>> double precision deltaRmin
>> deltaR=1d3
>> deltaRmin=1d3
>>
>> do iproc=1,maxsproc
>>
>> do i=nincoming+1,nexternal-1
>> if (idup(i,1,iproc).eq.13.or.idup(i,1,iproc).eq.-13)then
>> do j=i+1,nexternal
>> if (idup(j,1,iproc).eq.1023)then
>> p1(0:3)=P(0:3,i)
>> p2(0:3)=P(0:3,j)
>> deltaR=dsqrt(R2(p1,p2))
>> if (deltaR.lt.deltaRmin)then
>> deltaRmin=deltaR
>> endif
>> endif
>> enddo
>> endif
>> enddo
>> do i=nincoming+1,nexternal-1
>> if(idup(i,1,iproc).eq.1023)then
>> do j=i+1,nexternal
>> if(idup(j,1,iproc).eq.13.or.idup(j,1,iproc).eq.-13)then
>> p1(0:3)=P(0:3,i)
>> p2(0:3)=P(0:3,j)
>> deltaR=dsqrt(R2(p1,p2))
>> if (deltaR.lt.deltaRmin)then
>> deltaRmin=deltaR
>> endif
>> endif
>> enddo
>> endif
>> enddo
>> enddo
>> dummy_cuts=.true.
>>
>> if(deltaRmin.lt.4.0d-1) then
>> dummy_cuts=.false.
>> endif
>> return
>> end
>
>
>> On Sep 24, 2020, at 4:25 PM, Olivier Mattelaer <<email address hidden> <mailto:<email address hidden>>> wrote:
>>
>> Your question #693074 on MadGraph5_aMC@NLO changed:
>>
>> Status: Open => Answered
>>
>> Olivier Mattelaer proposed the following answer:
>> Hi,
>>
>> I have a few questions:
>> 1) How did you implement cut #3?
>> 2) What is the mass of the muon?
>> 3) Do you have a width to your dark photon?
>> 4) How many events are generated for those runs? What is the statistical error of such runs?
>>
>> Cheers,
>>
>> Olivier
>>
>> --
>> If this answers your question, please go to the following page to let us
>> know that it is solved:
>>
>> If you still need help, you can reply to this email or go to the
>> following page to enter your feedback:
>>
>> You received this question notification because you asked the question.
>

 Yao Yao (yyao1) said on 2020-09-24: #4

I check the param_card.dat, I set DECAY 1023 0.000000e+00 # wzp, which is the width of the dark photon.

> On Sep 24, 2020, at 4:45 PM, Yao Yao <email address hidden> wrote:
>
> Your question #693074 on MadGraph5_aMC@NLO changed:
>
> You gave more information on the question:
> Oops, sorry a type on the last cross section, it should be e+76 for the error.
>> On Sep 24, 2020, at 4:40 PM, Yao Yao <<email address hidden> <mailto:<email address hidden>>> wrote:
>>
>> Hi Olivier,
>>
>> 1) for cut #3, I asked this question to you before in [Question #686216]: and solved it by writing some code in dummy_fct.f. (Please see the quote)
>> 2) The mass of the muon should be standard, about 105 MeV. It should be set in HAHM model, and I didn't touch it.
>> 3) I don’t set a width for dark photon, or maybe I should?
>> 4) I generated 50000 events for each of them. The statistical error of the cross section for each is :
>>
>> 0.1601 +/- 0.04335 pb
>> 1.268e+31 +/- 1.792e+28 pb
>> 4.727e+79 +/- 5.962e+86 pb
>>
>> Yao
>>
>>> Hi Olivier,
>>>
>>> I got Madgraph running fast now when generating the events! I write the following code in dummy_fct.f. It implies that a dark photon needs to be deltaR>0.4 from any muon. Please help me check if the code seems good?
>>>
>>>
>>> include 'maxamps.inc'
>>> integer iproc,i,j
>>> integer idup(nexternal,maxproc,maxsproc)
>>> integer mothup(2,nexternal)
>>> integer icolup(2,nexternal,maxflow,maxsproc)
>>> include 'leshouche.inc'
>>>
>>> double precision p1(0:3),p2(0:3)
>>> double precision R2
>>> double precision deltaR
>>> double precision deltaRmin
>>> deltaR=1d3
>>> deltaRmin=1d3
>>>
>>> do iproc=1,maxsproc
>>>
>>> do i=nincoming+1,nexternal-1
>>> if (idup(i,1,iproc).eq.13.or.idup(i,1,iproc).eq.-13)then
>>> do j=i+1,nexternal
>>> if (idup(j,1,iproc).eq.1023)then
>>> p1(0:3)=P(0:3,i)
>>> p2(0:3)=P(0:3,j)
>>> deltaR=dsqrt(R2(p1,p2))
>>> if (deltaR.lt.deltaRmin)then
>>> deltaRmin=deltaR
>>> endif
>>> endif
>>> enddo
>>> endif
>>> enddo
>>> do i=nincoming+1,nexternal-1
>>> if(idup(i,1,iproc).eq.1023)then
>>> do j=i+1,nexternal
>>> if(idup(j,1,iproc).eq.13.or.idup(j,1,iproc).eq.-13)then
>>> p1(0:3)=P(0:3,i)
>>> p2(0:3)=P(0:3,j)
>>> deltaR=dsqrt(R2(p1,p2))
>>> if (deltaR.lt.deltaRmin)then
>>> deltaRmin=deltaR
>>> endif
>>> endif
>>> enddo
>>> endif
>>> enddo
>>> enddo
>>> dummy_cuts=.true.
>>>
>>> if(deltaRmin.lt.4.0d-1) then
>>> dummy_cuts=.false.
>>> endif
>>> return
>>> end
>>
>>
>>> On Sep 24, 2020, at 4:25 PM, Olivier Mattelaer <<email address hidden> <mailto:<email address hidden>><mailto:<email address hidden> <mailto:<email address hidden>>>> wrote:
>>>
>>> Your question #693074 on MadGraph5_aMC@NLO changed:
>>>
>>> Status: Open => Answered
>>>
>>> Olivier Mattelaer proposed the following answer:
>>> Hi,
>>>
>>> I have a few questions:
>>> 1) How did you implement cut #3?
>>> 2) What is the mass of the muon?
>>> 3) Do you have a width to your dark photon?
>>> 4) How many events are generated for those runs? What is the statistical error of such runs?
>>>
>>> Cheers,
>>>
>>> Olivier
>>>
>>> --
>>> If this answers your question, please go to the following page to let us
>>> know that it is solved:
>>>
>>> If you still need help, you can reply to this email or go to the
>>> following page to enter your feedback:
>>>
>>> You received this question notification because you asked the question.
>>
>
> --
> You received this question notification because you asked the question.

 Olivier Mattelaer (olivier-mattelaer) said on 2020-09-28: #5

Looks like we are not using the same model since in my HAHM model, the second benchmark point is crashing at the model level wit h the following error:

Error detected in "generate_events run_01"
write debug file /Users/omattelaer/Documents/workspace/2.7.2_alternate/PROC_HAHM_darkphoton_iDM_new_UFO_0/run_01_tag_1_debug.log
str : Unable to evaluate mdl_th = (-mdl_MHinput__exp__2 + mdl_MHSinput__exp__2 + (2*cmath.atan(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*(mdl_MHinput - mdl_MHSinput))*cmath.sqrt((mdl_MHinput__exp__2 - mdl_MHSinput__exp__2)**2 - 4*mdl_kap__exp__2*mdl_v__exp__2*mdl_xi__exp__2))/cmath.pi)/(2.*mdl_kap*mdl_v*mdl_xi): raise error: complex division by zero

Which seems to indicate that the model (or some version of it) can be subject to some numerical issue of some mass are close to zero. Can you check if this is the case of your model? and if not could you provide me a link to your model --maybe you did in another thread but I have so many thread that I will not be able to find it back--

Cheers,

Olivier

 Yao Yao (yyao1) said on 2020-09-28: #6

Hi Olivier,

I am using the model downloaded from this webpage: http://insti.physics.sunysb.edu/~curtin/hahm_mg.html <http://insti.physics.sunysb.edu/~curtin/hahm_mg.html>
The UFO file that I use is within the HVHM_MG5model_v3, and is called HAHM_variableMV_v3_UFO.

I also attach the param_card, proc_card and run_card.

Thank you so much for helping me checking this.
Yao

————————————————————— proc_card—————————————————————————
import model --modelname HAHM_variableMW_v3_UFO
define p = g u c d s b u~ c~ d~ s~ b~
define higgs = h hs
generate p p > m+ m- zp /higgs
output darkPhoton_m5e-19_df0p01_dR0p4

————————————————————— param_card—————————————————————————
######################################################################
## PARAM_CARD AUTOMATICALY GENERATED BY MG5 ####
######################################################################
###################################
## INFORMATION FOR CKMBLOCK
###################################
BLOCK CKMBLOCK #
1 4.880000e-01 # cabi
###################################
## INFORMATION FOR GAUGEMASS
###################################
BLOCK GAUGEMASS #
1 9.118800e+01 # mzinput
###################################
## INFORMATION FOR HIDDEN
###################################
BLOCK HIDDEN #
1 5.000000e-15 # mzdinput
2 1.000000e+04 # mhsinput
3 1.000000e-01 # epsilon
4 1.000000e-10 # kap
5 1.279000e+02 # axm1
###################################
## INFORMATION FOR HIGGS
###################################
BLOCK HIGGS #
1 1.250000e+02 # mhinput
###################################
## INFORMATION FOR MASS
###################################
BLOCK MASS #
4 1.420000e+00 # mc
5 4.700000e+00 # mb
6 1.743000e+02 # mt
11 5.110000e-04 # me
13 1.057000e-01 # mm
15 1.777000e+00 # mta
1 0.000000e+00 # d : 0.0
2 0.000000e+00 # u : 0.0
3 0.000000e+00 # s : 0.0
12 0.000000e+00 # ve : 0.0
14 0.000000e+00 # vm : 0.0
16 0.000000e+00 # vt : 0.0
21 0.000000e+00 # g : 0.0
22 0.000000e+00 # a : 0.0
23 9.118800e+01 # z : mzinput
24 8.015873e+01 # w+ : cw*mz0
25 1.250000e+02 # h : mhinput
35 1.000000e+04 # hs : mhsinput
1023 5.000000e-15 # zp : mzdinput
###################################
## INFORMATION FOR SMINPUTS
###################################
BLOCK SMINPUTS #
1 2.250000e-01 # swsq
2 1.279000e+02 # aewm1
3 1.166390e-05 # gf
4 1.180000e-01 # as
###################################
## INFORMATION FOR YUKAWA
###################################
BLOCK YUKAWA #
4 1.420000e+00 # ymc
5 4.700000e+00 # ymb
6 1.743000e+02 # ymt
11 5.110000e-04 # ymel
13 1.057000e-01 # ymmu
15 1.777000e+00 # ymtau
###################################
## INFORMATION FOR DECAY
###################################
DECAY 6 1.508336e+00 # wt
DECAY 23 2.441404e+00 # wz
DECAY 24 2.047600e+00 # ww
DECAY 25 2.822990e-03 # wh
DECAY 35 5.237950e+00 # whs
DECAY 1023 auto # wzp
DECAY 1 0.000000e+00 # d : 0.0
DECAY 2 0.000000e+00 # u : 0.0
DECAY 3 0.000000e+00 # s : 0.0
DECAY 4 0.000000e+00 # c : 0.0
DECAY 5 0.000000e+00 # b : 0.0
DECAY 11 0.000000e+00 # e- : 0.0
DECAY 12 0.000000e+00 # ve : 0.0
DECAY 13 0.000000e+00 # m- : 0.0
DECAY 14 0.000000e+00 # vm : 0.0
DECAY 15 0.000000e+00 # tt- : 0.0
DECAY 16 0.000000e+00 # vt : 0.0
DECAY 21 0.000000e+00 # g : 0.0
DECAY 22 0.000000e+00 # a : 0.0
###################################
## INFORMATION FOR QNUMBERS 1023
###################################
BLOCK QNUMBERS 1023 # zp
1 0 # 3 times electric charge
2 3 # number of spin states (2s+1)
3 1 # colour rep (1: singlet, 3: triplet, 8: octet)
4 0 # particle/antiparticle distinction (0=own anti)
###################################
## INFORMATION FOR QNUMBERS 35
###################################
BLOCK QNUMBERS 35 # hs
1 0 # 3 times electric charge
2 1 # number of spin states (2s+1)
3 1 # colour rep (1: singlet, 3: triplet, 8: octet)
4 0 # particle/antiparticle distinction (0=own anti)

————————————————————— run_card—————————————————————————
#*********************************************************************
# *
# run_card.dat MadEvent *
# *
# This file is used to set the parameters of the run. *
# *
# Some notation/conventions: *
# *
# Lines starting with a '# ' are info or comments *
# *
# mind the format: value = variable ! comment *
# *
# To display more options, you can type the command: *
# update full_run_card *
#*********************************************************************
#
#*******************
# Running parameters
#*******************
#
#*********************************************************************
# Tag name for the run (one word) *
#*********************************************************************
tag_1 = run_tag ! name of the run
#*********************************************************************
# Number of events and rnd seed *
# Warning: Do not generate more than 1M events in a single run *
# If you want to run Pythia, avoid more than 50k events in a run. *
#*********************************************************************
50000 = nevents ! Number of unweighted events requested
0 = iseed ! rnd seed (0=assigned automatically=default))
#*********************************************************************
# Collider type and energy *
# lpp: 0=No PDF, 1=proton, -1=antiproton, 2=photon from proton, *
# 3=photon from electron *
#*********************************************************************
1 = lpp1 ! beam 1 type
1 = lpp2 ! beam 2 type
6500.0 = ebeam1 ! beam 1 total energy in GeV
6500.0 = ebeam2 ! beam 2 total energy in GeV
# To see polarised beam options: type "update beam_pol"
#*********************************************************************
# PDF CHOICE: this automatically fixes also alpha_s and its evol. *
#*********************************************************************
nn23lo1 = pdlabel ! PDF set
230000 = lhaid ! if pdlabel=lhapdf, this is the lhapdf number
# To see heavy ion options: type "update ion_pdf"
#*********************************************************************
# Renormalization and factorization scales *
#*********************************************************************
False = fixed_ren_scale ! if .true. use fixed ren scale
False = fixed_fac_scale ! if .true. use fixed fac scale
91.188 = scale ! fixed ren scale
91.188 = dsqrt_q2fact1 ! fixed fact scale for pdf1
91.188 = dsqrt_q2fact2 ! fixed fact scale for pdf2
-1 = dynamical_scale_choice ! Choose one of the preselected dynamical choices
1.0 = scalefact ! scale factor for event-by-event scales
#*********************************************************************
# Type and output format
#*********************************************************************
False = gridpack !True = setting up the grid pack
-1.0 = time_of_flight ! threshold (in mm) below which the invariant livetime is not written (-1 means not written)
3.0 = lhe_version ! Change the way clustering information pass to shower.
True = clusinfo ! include clustering tag in output
average = event_norm ! average/sum. Normalization of the weight in the LHEF

#*********************************************************************
# Matching parameter (MLM only)
#*********************************************************************
0 = ickkw ! 0 no matching, 1 MLM
1.0 = alpsfact ! scale factor for QCD emission vx
False = chcluster ! cluster only according to channel diag
5 = asrwgtflavor ! highest quark flavor for a_s reweight
False = auto_ptj_mjj ! Automatic setting of ptj and mjj if xqcut >0
! (turn off for VBF and single top processes)
0.0 = xqcut ! minimum kt jet measure between partons
#*********************************************************************
#
#*********************************************************************
# handling of the helicities:
# 0: sum over all helicities
# 1: importance sampling over helicities
#*********************************************************************
0 = nhel ! using helicities importance sampling or not.
#*********************************************************************
# Generation bias, check the wiki page below for more information: *
#*********************************************************************
None = bias_module ! Bias type of bias, [None, ptj_bias, -custom_folder-]
{} = bias_parameters ! Specifies the parameters of the module.
#
#*******************************
# Parton level cuts definition *
#*******************************
#
#
#*********************************************************************
# BW cutoff (M+/-bwcutoff*Gamma) ! Define on/off-shell for "\$" and decay
#*********************************************************************
15.0 = bwcutoff ! (M+/-bwcutoff*Gamma)
#*********************************************************************
# Apply pt/E/eta/dr/mij/kt_durham cuts on decay products or not
# (note that etmiss/ptll/ptheavy/ht/sorted cuts always apply)
#*********************************************************************
False = cut_decays ! Cut decay products
#*********************************************************************
# Standard Cuts *
#*********************************************************************
# Minimum and maximum pt's (for max, -1 means no cut) *
#*********************************************************************
20.0 = ptj ! minimum pt for the jets
0.0 = ptb ! minimum pt for the b
10.0 = pta ! minimum pt for the photons
10.0 = ptl ! minimum pt for the charged leptons
0.0 = misset ! minimum missing Et (sum of neutrino's momenta)
-1.0 = ptjmax ! maximum pt for the jets
-1.0 = ptbmax ! maximum pt for the b
-1.0 = ptamax ! maximum pt for the photons
-1.0 = ptlmax ! maximum pt for the charged leptons
-1.0 = missetmax ! maximum missing Et (sum of neutrino's momenta)
{1023: 10} = pt_min_pdg ! pt cut for other particles (use pdg code). Applied on particle and anti-particle
{} = pt_max_pdg ! pt cut for other particles (syntax e.g. {6: 100, 25: 50})
#*********************************************************************
# Minimum and maximum E's (in the center of mass frame) *
#*********************************************************************
0.0 = ej ! minimum E for the jets
0.0 = eb ! minimum E for the b
0.0 = ea ! minimum E for the photons
0.0 = el ! minimum E for the charged leptons
-1.0 = ejmax ! maximum E for the jets
-1.0 = ebmax ! maximum E for the b
-1.0 = eamax ! maximum E for the photons
-1.0 = elmax ! maximum E for the charged leptons
{} = e_min_pdg ! E cut for other particles (use pdg code). Applied on particle and anti-particle
{} = e_max_pdg ! E cut for other particles (syntax e.g. {6: 100, 25: 50})
#*********************************************************************
# Maximum and minimum absolute rapidity (for max, -1 means no cut) *
#*********************************************************************
5.0 = etaj ! max rap for the jets
-1.0 = etab ! max rap for the b
2.5 = etaa ! max rap for the photons
2.5 = etal ! max rap for the charged leptons
0.0 = etajmin ! min rap for the jets
0.0 = etabmin ! min rap for the b
0.0 = etaamin ! min rap for the photons
0.0 = etalmin ! main rap for the charged leptons
{} = eta_min_pdg ! rap cut for other particles (use pdg code). Applied on particle and anti-particle
{} = eta_max_pdg ! rap cut for other particles (syntax e.g. {6: 2.5, 23: 5})
#*********************************************************************
# Minimum and maximum DeltaR distance *
#*********************************************************************
0.4 = drjj ! min distance between jets
0.0 = drbb ! min distance between b's
0.4 = drll ! min distance between leptons
0.4 = draa ! min distance between gammas
0.0 = drbj ! min distance between b and jet
0.4 = draj ! min distance between gamma and jet
0.4 = drjl ! min distance between jet and lepton
0.0 = drab ! min distance between gamma and b
0.0 = drbl ! min distance between b and lepton
0.4 = dral ! min distance between gamma and lepton
-1.0 = drjjmax ! max distance between jets
-1.0 = drbbmax ! max distance between b's
-1.0 = drllmax ! max distance between leptons
-1.0 = draamax ! max distance between gammas
-1.0 = drbjmax ! max distance between b and jet
-1.0 = drajmax ! max distance between gamma and jet
-1.0 = drjlmax ! max distance between jet and lepton
-1.0 = drabmax ! max distance between gamma and b
-1.0 = drblmax ! max distance between b and lepton
-1.0 = dralmax ! maxdistance between gamma and lepton
#*********************************************************************
# Minimum and maximum invariant mass for pairs *
# WARNING: for four lepton final state mmll cut require to have *
# different lepton masses for each flavor! *
#*********************************************************************
0.0 = mmjj ! min invariant mass of a jet pair
0.0 = mmbb ! min invariant mass of a b pair
0.0 = mmaa ! min invariant mass of gamma gamma pair
10.0 = mmll ! min invariant mass of l+l- (same flavour) lepton pair
-1.0 = mmjjmax ! max invariant mass of a jet pair
-1.0 = mmbbmax ! max invariant mass of a b pair
-1.0 = mmaamax ! max invariant mass of gamma gamma pair
-1.0 = mmllmax ! max invariant mass of l+l- (same flavour) lepton pair
{} = mxx_min_pdg ! min invariant mass of a pair of particles X/X~ (e.g. {6:250})
{'default': False} = mxx_only_part_antipart ! if True the invariant mass is applied only
! to pairs of particle/antiparticle and not to pairs of the same pdg codes.
#*********************************************************************
# Minimum and maximum invariant mass for all letpons *
#*********************************************************************
0.0 = mmnl ! min invariant mass for all letpons (l+- and vl)
-1.0 = mmnlmax ! max invariant mass for all letpons (l+- and vl)
#*********************************************************************
# Minimum and maximum pt for 4-momenta sum of leptons *
#*********************************************************************
0.0 = ptllmin ! Minimum pt for 4-momenta sum of leptons(l and vl)
-1.0 = ptllmax ! Maximum pt for 4-momenta sum of leptons(l and vl)
#*********************************************************************
# Inclusive cuts *
#*********************************************************************
0.0 = ptheavy ! minimum pt for at least one heavy final state
0.0 = xptj ! minimum pt for at least one jet
0.0 = xptb ! minimum pt for at least one b
0.0 = xpta ! minimum pt for at least one photon
0.0 = xptl ! minimum pt for at least one charged lepton
#*********************************************************************
# Control the pt's of the jets sorted by pt *
#*********************************************************************
0.0 = ptj1min ! minimum pt for the leading jet in pt
0.0 = ptj2min ! minimum pt for the second jet in pt
0.0 = ptj3min ! minimum pt for the third jet in pt
0.0 = ptj4min ! minimum pt for the fourth jet in pt
-1.0 = ptj1max ! maximum pt for the leading jet in pt
-1.0 = ptj2max ! maximum pt for the second jet in pt
-1.0 = ptj3max ! maximum pt for the third jet in pt
-1.0 = ptj4max ! maximum pt for the fourth jet in pt
0 = cutuse ! reject event if fails any (0) / all (1) jet pt cuts
#*********************************************************************
# Control the pt's of leptons sorted by pt *
#*********************************************************************
0.0 = ptl1min ! minimum pt for the leading lepton in pt
0.0 = ptl2min ! minimum pt for the second lepton in pt
0.0 = ptl3min ! minimum pt for the third lepton in pt
0.0 = ptl4min ! minimum pt for the fourth lepton in pt
-1.0 = ptl1max ! maximum pt for the leading lepton in pt
-1.0 = ptl2max ! maximum pt for the second lepton in pt
-1.0 = ptl3max ! maximum pt for the third lepton in pt
-1.0 = ptl4max ! maximum pt for the fourth lepton in pt
#*********************************************************************
# Control the Ht(k)=Sum of k leading jets *
#*********************************************************************
0.0 = htjmin ! minimum jet HT=Sum(jet pt)
-1.0 = htjmax ! maximum jet HT=Sum(jet pt)
0.0 = ihtmin !inclusive Ht for all partons (including b)
-1.0 = ihtmax !inclusive Ht for all partons (including b)
0.0 = ht2min ! minimum Ht for the two leading jets
0.0 = ht3min ! minimum Ht for the three leading jets
0.0 = ht4min ! minimum Ht for the four leading jets
-1.0 = ht2max ! maximum Ht for the two leading jets
-1.0 = ht3max ! maximum Ht for the three leading jets
-1.0 = ht4max ! maximum Ht for the four leading jets
#***********************************************************************
# Photon-isolation cuts, according to hep-ph/9801442 *
# When ptgmin=0, all the other parameters are ignored *
# When ptgmin>0, pta and draj are not going to be used *
#***********************************************************************
0.0 = ptgmin ! Min photon transverse momentum
0.4 = R0gamma ! Radius of isolation code
1.0 = xn ! n parameter of eq.(3.4) in hep-ph/9801442
1.0 = epsgamma ! epsilon_gamma parameter of eq.(3.4) in hep-ph/9801442
True = isoEM ! isolate photons from EM energy (photons and leptons)
#*********************************************************************
# WBF cuts *
#*********************************************************************
0.0 = xetamin ! minimum rapidity for two jets in the WBF case
0.0 = deltaeta ! minimum rapidity for two jets in the WBF case
#***********************************************************************
# Turn on either the ktdurham or ptlund cut to activate *
# CKKW(L) merging with Pythia8 [arXiv:1410.3012, arXiv:1109.4829] *
#***********************************************************************
-1.0 = ktdurham
0.4 = dparameter
-1.0 = ptlund
1, 2, 3, 4, 5, 6, 21 = pdgs_for_merging_cut ! PDGs for two cuts above
#*********************************************************************
# maximal pdg code for quark to be considered as a light jet *
# (otherwise b cuts are applied) *
#*********************************************************************
5 = maxjetflavor ! Maximum jet pdg code
#*********************************************************************
#
#*********************************************************************
# Store info for systematics studies *
# WARNING: Do not use for interference type of computation *
#*********************************************************************
True = use_syst ! Enable systematics studies
#
systematics = systematics_program ! none, systematics [python], SysCalc [depreceted, C++]
['--mur=0.5,1,2', '--muf=0.5,1,2', '--pdf=errorset'] = systematics_arguments ! see: https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/Systematics#Systematicspythonmodule
# Syscalc is deprecated but to see the associate options type'update syscalc'

> On Sep 28, 2020, at 10:15 PM, Olivier Mattelaer <email address hidden> wrote:
>
> Your question #693074 on MadGraph5_aMC@NLO changed:
>
> Olivier Mattelaer posted a new comment:
> Looks like we are not using the same model since in my HAHM model, the
> second benchmark point is crashing at the model level wit h the
> following error:
>
> Error detected in "generate_events run_01"
> write debug file /Users/omattelaer/Documents/workspace/2.7.2_alternate/PROC_HAHM_darkphoton_iDM_new_UFO_0/run_01_tag_1_debug.log
> str : Unable to evaluate mdl_th = (-mdl_MHinput__exp__2 + mdl_MHSinput__exp__2 + (2*cmath.atan(10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000*(mdl_MHinput - mdl_MHSinput))*cmath.sqrt((mdl_MHinput__exp__2 - mdl_MHSinput__exp__2)**2 - 4*mdl_kap__exp__2*mdl_v__exp__2*mdl_xi__exp__2))/cmath.pi)/(2.*mdl_kap*mdl_v*mdl_xi): raise error: complex division by zero
>
> Which seems to indicate that the model (or some version of it) can be
> subject to some numerical issue of some mass are close to zero. Can you
> check if this is the case of your model? and if not could you provide me
> a link to your model --maybe you did in another thread but I have so
> many thread that I will not be able to find it back--
>
> Cheers,
>
> Olivier
>
> --
> You received this question notification because you asked the question.

 Olivier Mattelaer (olivier-mattelaer) said on 2020-09-29: #7

Hi,

For what I see the issue is that your issue is related to a very large contribution of
M_mumu > 1TeV

I can safely confirm that this is not an issue of the phase-space integrator.
I have compare the phase-space integrator of 2.8.0 with the quite different one that I'm building for 2.9.0 and they behave the same way.

I can see that they are some issue with the matrix-element evaluation in that part of the phase-space if i check the validity of the code for the first events that is generated at such high energy:

Lorentz invariance results:
Process Min element Max element Relative diff. Result
u~ u > m+ m- zp 6.1451618890e-11 6.1451619029e-11 2.2580989883e-09 Failed
JAMP 0 7.3741942668e-10 7.3741942834e-10 2.2580988480e-09 Failed
Summary: 0/1 passed, 1/1 failed
Failed processes: u~ u > m+ m- zp
Process permutation results:
Process Min element Max element Relative diff. Result
u~ u > m+ m- zp 6.1397114979e-11 6.1397114979e-11 4.2101968718e-16 Passed
Summary: 1/1 passed, 0/1 failed

You can see that the Lorentz invariance of the matrix-element is quite bad, which is typically related to a gauge/unitary issue of the model/process under consideration. So you are right to say that this is not physical but I do not think that the issue is on our side.

Cheers,

Olivier

 Yao Yao (yyao1) said on 2020-10-01: #8

Hi Olivier,

I just realized that if I remove the higgs exclusion from the process p p > m+ m- zp / higgs, (which means I am generating p p > m+ m- zp) the cross section becomes zero and there are errors generated as follows:

Generating 10000 events with run name run_05
survey run_05
INFO: compile directory
compile Source Directory
Using random number seed offset = 57
INFO: Running Survey
Creating Jobs
Working on SubProcesses
INFO: P1_ccx_mpmmzp
INFO: P1_bbx_mpmmzp
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG
...
INFO: End survey
refine 10000
Creating Jobs
INFO: Refine results to 10000
INFO: Generating 10000.0 unweigthed events.
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.
Zero result detected: See https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/FAQ-General-14

So do you think that excluding higgs and dark higgs destroys the gauge invariant in this situation, and anything that we probably can modify in the cards that can help solving this issue?

Best,
Yao

> On Sep 29, 2020, at 11:30 AM, Olivier Mattelaer <email address hidden> wrote:
>
> Your question #693074 on MadGraph5_aMC@NLO changed:
>
> Status: Open => Answered
>
> Olivier Mattelaer proposed the following answer:
> Hi,
>
> For what I see the issue is that your issue is related to a very large contribution of
> M_mumu > 1TeV
>
> I can safely confirm that this is not an issue of the phase-space integrator.
> I have compare the phase-space integrator of 2.8.0 with the quite different one that I'm building for 2.9.0 and they behave the same way.
>
> I can see that they are some issue with the matrix-element evaluation in
> that part of the phase-space if i check the validity of the code for the
> first events that is generated at such high energy:
>
> Lorentz invariance results:
> Process Min element Max element Relative diff. Result
> u~ u > m+ m- zp 6.1451618890e-11 6.1451619029e-11 2.2580989883e-09 Failed
> JAMP 0 7.3741942668e-10 7.3741942834e-10 2.2580988480e-09 Failed
> Summary: 0/1 passed, 1/1 failed
> Failed processes: u~ u > m+ m- zp
> Process permutation results:
> Process Min element Max element Relative diff. Result
> u~ u > m+ m- zp 6.1397114979e-11 6.1397114979e-11 4.2101968718e-16 Passed
> Summary: 1/1 passed, 0/1 failed
>
> You can see that the Lorentz invariance of the matrix-element is quite
> bad, which is typically related to a gauge/unitary issue of the
> model/process under consideration. So you are right to say that this is
> not physical but I do not think that the issue is on our side.
>
> Cheers,
>
> Olivier
>
> --
> If this answers your question, please go to the following page to let us
> know that it is solved:
>
> If you still need help, you can reply to this email or go to the
> following page to enter your feedback:
>
> You received this question notification because you asked the question.

 Olivier Mattelaer (olivier-mattelaer) said on 2020-10-02: #9

HI,

I do not know your model, so I can not comment on the gauge invariant issue of the model implementation/process.

Concerning:
Note: The following floating-point exceptions are signalling: IEEE_INVALID_FLAG
Those are typically raised in presence of square root of negative numbe, division by zero,... in the evaluation of some coupling.
Therefore, This typically points to an issue in the model implementation which leads to this zero result

Cheers,

Olivier