Awful runtime scaling in ZZjj event generation

Hi Philipp,

2.3.3 should improve the speed of various processes (and in particular the scaling with the requested precision/number of event).
But in any case, we do not care of scaling time for such few number of event.

Cheers,

Olivier

On 10 Oct 2015, at 13:46, Philipp Pigard <email address hidden> wrote:

> New question #272287 on MadGraph5_aMC@NLO:
> https://answers.launchpad.net/mg5amcnlo/+question/272287
>
> Dear experts,
>
> I am using MadGraph5_aMC@NLO 2.3.2.2 to generate events of the electroweak production of ZZjj. The time needed to generate one event takes about 2 seconds for 100 events but increases to a minute when I request 500 events. These numbers are based on running on a single core.
>
> I am now wondering whether this kind of bad scaling is expected or if I made a mistake in setting up the production. The process and run card are copied below.
>
> Best regards
> Philipp
>
> process card:
>
> set group_subprocesses Auto
> set ignore_six_quark_processes False
> set loop_optimized_output True
> set complex_mass_scheme False
> import model sm-ckm_no_b_mass
> define p = g u c d s b u~ c~ d~ s~ b~
> define j = p
> define l+ = e+ mu+ ta+
> define l- = e- mu- ta-
> define vl = ve vm vt
> define vl~ = ve~ vm~ vt~
>
> generate p p > z z j j QED=4 QCD=0, z > e+ e-, z > mu+ mu- @ 1
> output ZZjj_ewk
>
> run card extract:
>
> 'lhapdf' = pdlabel ! PDF set
> 263000 = lhaid ! if pdlabel=lhapdf, this is the lhapdf number
> #*********************************************************************
> # Renormalization and factorization scales *
> #*********************************************************************
> F = fixed_ren_scale ! if .true. use fixed ren scale
> F = fixed_fac_scale ! if .true. use fixed fac scale
> 91.1880 = scale ! fixed ren scale
> 91.1880 = dsqrt_q2fact1 ! fixed fact scale for pdf1
> 91.1880 = dsqrt_q2fact2 ! fixed fact scale for pdf2
> 1 = scalefact ! scale factor for event-by-event scales
> #*********************************************************************
> # Matching - Warning! ickkw > 1 is still beta
> #*********************************************************************
> 0 = ickkw ! 0 no matching, 1 MLM, 2 CKKW matching
> 1 = highestmult ! for ickkw=2, highest mult group
> 1 = ktscheme ! for ickkw=1, 1 Durham kT, 2 Pythia pTE
> 1 = alpsfact ! scale factor for QCD emission vx
> F = chcluster ! cluster only according to channel diag
> T = pdfwgt ! for ickkw=1, perform pdf reweighting
> 5 = asrwgtflavor ! highest quark flavor for a_s reweight
> T = clusinfo ! include clustering tag in output
> 3.0 = lhe_version ! Change the way clustering information pass to shower.
>
> F = auto_ptj_mjj ! Automatic setting of ptj and mjj
> #**********************************************************
> #
> #**********************************
> # BW cutoff (M+/-bwcutoff*Gamma)
> #**********************************
> 15000 = bwcutoff ! (M+/-bwcutoff*Gamma)
> #**********************************************************
> # Apply pt/E/eta/dr/mij cuts on decay products or not
> # (note that etmiss/ptll/ptheavy/ht/sorted cuts always apply)
> #**********************************************************
> T = cut_decays ! Cut decay products
> #*************************************************************
> # Number of helicities to sum per event (0 = all helicities)
> # 0 gives more stable result, but longer run time (needed for
> # long decay chains e.g.).
> # Use >=2 if most helicities contribute, e.g. pure QCD.
> #*************************************************************
> 0 = nhel ! Number of helicities used per event
> #*******************
> # Standard Cuts
> #*******************
> #
> #*********************************************************************
> # Minimum and maximum pt's (for max, -1 means no cut) *
> #*********************************************************************
> 10 = ptj ! minimum pt for the jets
> 0 = ptb ! minimum pt for the b
> 0 = pta ! minimum pt for the photons
> 0 = ptl ! minimum pt for the charged leptons
> 0 = misset ! minimum missing Et (sum of neutrino's momenta)
> 0 = ptheavy ! minimum pt for one heavy final state
> 1.0 = ptonium ! minimum pt for the quarkonium states
> -1 = ptjmax ! maximum pt for the jets
> -1 = ptbmax ! maximum pt for the b
> -1 = ptamax ! maximum pt for the photons
> -1 = ptlmax ! maximum pt for the charged leptons
> -1 = missetmax ! maximum missing Et (sum of neutrino's momenta)
> #*********************************************************************
> # Minimum and maximum E's (in the center of mass frame) *
> #*********************************************************************
> 0 = ej ! minimum E for the jets
> 0 = eb ! minimum E for the b
> 0 = ea ! minimum E for the photons
> 0 = el ! minimum E for the charged leptons
> -1 = ejmax ! maximum E for the jets
> -1 = ebmax ! maximum E for the b
> -1 = eamax ! maximum E for the photons
> -1 = elmax ! maximum E for the charged leptons
> #*********************************************************************
> # Maximum and minimum absolute rapidity (for max, -1 means no cut) *
> #*********************************************************************
> -1 = etaj ! max rap for the jets
> -1 = etab ! max rap for the b
> 2.5 = etaa ! max rap for the photons
> -1 = etal ! max rap for the charged leptons
> 0.6 = etaonium ! max rap for the quarkonium states
> 0 = etajmin ! min rap for the jets
> 0 = etabmin ! min rap for the b
> 0 = etaamin ! min rap for the photons
> 0 = etalmin ! main rap for the charged leptons
> #*********************************************************************
> # 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.0 = draa ! min distance between gammas
> 0.0 = drbj ! min distance between b and jet
> 0.0 = 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.0 = dral ! min distance between gamma and lepton
> -1 = drjjmax ! max distance between jets
> -1 = drbbmax ! max distance between b's
> -1 = drllmax ! max distance between leptons
> -1 = draamax ! max distance between gammas
> -1 = drbjmax ! max distance between b and jet
> -1 = drajmax ! max distance between gamma and jet
> -1 = drjlmax ! max distance between jet and lepton
> -1 = drabmax ! max distance between gamma and b
> -1 = drblmax ! max distance between b and lepton
> -1 = 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! *
> #*********************************************************************
> 100 = mmjj ! min invariant mass of a jet pair
> 0 = mmbb ! min invariant mass of a b pair
> 0 = mmaa ! min invariant mass of gamma gamma pair
> 0 = mmll ! min invariant mass of l+l- (same flavour) lepton pair
> -1 = mmjjmax ! max invariant mass of a jet pair
> -1 = mmbbmax ! max invariant mass of a b pair
> -1 = mmaamax ! max invariant mass of gamma gamma pair
> -1 = mmllmax ! max invariant mass of l+l- (same flavour) lepton pair
> #*********************************************************************
> # Minimum and maximum invariant mass for all letpons *
> #*********************************************************************
> 0 = mmnl ! min invariant mass for all letpons (l+- and vl)
> -1 = mmnlmax ! max invariant mass for all letpons (l+- and vl)
> #*********************************************************************
> # Minimum and maximum pt for 4-momenta sum of leptons *
> #*********************************************************************
> 0 = ptllmin ! Minimum pt for 4-momenta sum of leptons(l and vl)
> -1 = ptllmax ! Maximum pt for 4-momenta sum of leptons(l and vl)
> #*********************************************************************
> # Inclusive cuts *
> #*********************************************************************
> 0 = xptj ! minimum pt for at least one jet
> 0 = xptb ! minimum pt for at least one b
> 0 = xpta ! minimum pt for at least one photon
> 0 = xptl ! minimum pt for at least one charged lepton
> #*********************************************************************
> # Control the pt's of the jets sorted by pt *
> #*********************************************************************
> 0 = ptj1min ! minimum pt for the leading jet in pt
> 0 = ptj2min ! minimum pt for the second jet in pt
> 0 = ptj3min ! minimum pt for the third jet in pt
> 0 = ptj4min ! minimum pt for the fourth jet in pt
> -1 = ptj1max ! maximum pt for the leading jet in pt
> -1 = ptj2max ! maximum pt for the second jet in pt
> -1 = ptj3max ! maximum pt for the third jet in pt
> -1 = 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 = ptl1min ! minimum pt for the leading lepton in pt
> 0 = ptl2min ! minimum pt for the second lepton in pt
> 0 = ptl3min ! minimum pt for the third lepton in pt
> 0 = ptl4min ! minimum pt for the fourth lepton in pt
> -1 = ptl1max ! maximum pt for the leading lepton in pt
> -1 = ptl2max ! maximum pt for the second lepton in pt
> -1 = ptl3max ! maximum pt for the third lepton in pt
> -1 = ptl4max ! maximum pt for the fourth lepton in pt
> #*********************************************************************
> # Control the Ht(k)=Sum of k leading jets *
> #*********************************************************************
> 0 = htjmin ! minimum jet HT=Sum(jet pt)
> -1 = htjmax ! maximum jet HT=Sum(jet pt)
> 0 = ihtmin !inclusive Ht for all partons (including b)
> -1 = ihtmax !inclusive Ht for all partons (including b)
> 0 = ht2min ! minimum Ht for the two leading jets
> 0 = ht3min ! minimum Ht for the three leading jets
> 0 = ht4min ! minimum Ht for the four leading jets
> -1 = ht2max ! maximum Ht for the two leading jets
> -1 = ht3max ! maximum Ht for the three leading jets
> -1 = 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 = 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 = xetamin ! minimum rapidity for two jets in the WBF case
> 0.0 = deltaeta ! minimum rapidity for two jets in the WBF case
> #*********************************************************************
> # KT DURHAM CUT *
> #*********************************************************************
> -1 = ktdurham
> 0.4 = dparameter
> #*********************************************************************
> # maximal pdg code for quark to be considered as a light jet *
> # (otherwise b cuts are applied) *
> #*********************************************************************
> 5 = maxjetflavor ! Maximum jet pdg code
> #*********************************************************************
> # Jet measure cuts *
> #*********************************************************************
> 0 = xqcut ! minimum kt jet measure between partons
> #*********************************************************************
> # Store info for systematics studies *
> # WARNING: If use_syst is T, matched Pythia output is *
> # meaningful ONLY if plotted taking matchscale *
> # reweighting into account! *
> #*********************************************************************
> T = use_syst ! Enable systematics studies
> #
> #**************************************
> # Parameter of the systematics study
> # will be used by SysCalc (if installed)
> #**************************************
> #
> 0.5 1 2 = sys_scalefact # factorization/renormalization scale factor
> 0.5 1 2 = sys_alpsfact # \alpha_s emission scale factors
> 30 50 = sys_matchscale # variation of merging scale
> # PDF sets and number of members (0 or none for all members).
> NNPDF30_lo_as_0130.LHgrid = sys_pdf # matching scales
> # MSTW2008nlo68cl.LHgrid 1 = sys_pdf
>
>
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Hello Olivier,

thank you for you quick response! Please let me ask a follow up question :

I understand that the number of events I am generating is really small, but is this scaling expected? Because it is and it is not the result of an improper or idiotic configuration, it will put constraints on generating samples with proper statistics, as we would be limited to running production with less than 1k events per job.

I will explore the possibility of using a more recent version, but my current workflow is based on the version we have in CMSSW.

Again thank you for your support
Philipp

Dear Philip,

Do I understand that your plan is to submit a huge number of runs (running on one core).

Then we have a dedicated mode for that called the gridpack (which is the standard method used by both CMS and ATLAS MC group).

See details here:
https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/GridDevelopment

For the details of the difference with a standard jobs, see:
https://cp3.irmp.ucl.ac.be/projects/madgraph/wiki/IntroGrid

Cheers,

Olivier

On 10 Oct 2015, at 15:22, Philipp Pigard <email address hidden> wrote:

> Question #272287 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/272287
>
> Status: Answered => Open
>
> Philipp Pigard is still having a problem:
> Hello Olivier,
>
> thank you for you quick response! Please let me ask a follow up question
> :
>
> I understand that the number of events I am generating is really small,
> but is this scaling expected? Because it is and it is not the result of
> an improper or idiotic configuration, it will put constraints on
> generating samples with proper statistics, as we would be limited to
> running production with less than 1k events per job.
>
> I will explore the possibility of using a more recent version, but my
> current workflow is based on the version we have in CMSSW.
>
> Again thank you for your support
> Philipp
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Dear Olivier,

why isn't the scaling mentioned by Philipp worrisome?
Naively, with more events generated the average time decreases instead of increasing, and the change (from 2s/event to 60s/event) is quite large indeed.

I think that it's important to know whether such a change is expected or not, at least as an indirect validation of Philipp's local setup.

Best regards,

pietro

Dear Pietro,

> Naively, with more events generated the average time decreases instead of increasing,

The best possible program can technically have a linear scaling for a large number of event —not too large to avoid RAM/disk problem—.
This is possible only if you know in advance the best way to integrate (i.e. if you can keep your grid of integration fixed). This is one of the idea of the gridpack: prepare all
grid to have a better scaling.

Now in practise (without gridpack mode) you need to iterate on those grid to find on the flight the correct one.
This step force you to run your code in a series of iteration which breaks the linearity.
In our code, we double the number of probed events in each iteration up to reach the target of precision/number of events.
Asking 1 more event/changing the seed, can therefore (if you are unlucky) need one more iteration which double the total times for that channel by a factor of two.

Asking more event, can also require to make more computation since small channel can then contribute. Having more channel slows down the computation even more.
Since it is a new channel of integration, the convergence of such channel is independent of the other ones/… and the dependence in the total time is a priori arbitrarily large.
Also for security, we always request three iterations for a channel which again break the linearity at low number of events.

So I do not have any a priori expectation on the scaling of the running time for a such “small” number of events.
For large number of events, I would more expect to have a scaling as N^2 just because of the re-ordering of the event in a random order.
I have implemented a new method scaling like N in the latest version but this is not yet used for all type of generation (but for most of them) so you can still face this N^2 dependencies.

Now to have an idea of how much you are dominated by those above point, you can look at the luminosity obtained channel by channel.
The perfect behaviour should be a constant number which should be the sign of a perfect integration (which a linear scaling in the generation phase).
But in practise, you never generate enough event in order for us to reach such ideal.

Cheers,

Olivier

On 14 Oct 2015, at 13:22, pietro govoni <email address hidden> wrote:

> Question #272287 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/272287
>
> pietro govoni posted a new comment:
> Dear Olivier,
>
> why isn't the scaling mentioned by Philipp worrisome?
> Naively, with more events generated the average time decreases instead of increasing, and the change (from 2s/event to 60s/event) is quite large indeed.
>
> I think that it's important to know whether such a change is expected or
> not, at least as an indirect validation of Philipp's local setup.
>
> Best regards,
>
> pietro
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Dear Olivier,

thanks a lot for the detailed explanation.
So it seems to be that having a gridpack is not only in general faster, but even safer.
In fact, from this sentence of yours:

"Asking more event, can also require to make more computation since small channel can then contribute."

is seems to be that generating a small set of events might lead to neglecting some channels, and therefore that a produciton of 1M events might give a different outcome than lumping together 1000 productions of 1K events.
Is my interpretation correct?

Is it correct also to state that this problem does not exist in case of gridpacks?

Thanks a lot!!

pietro

Dear Pietro,

> is seems to be that generating a small set of events might lead to neglecting some channels, and therefore that a produciton of 1M events might give a different outcome than lumping together 1000 productions of 1K events.

First 1000 production of 1K events do not correspond to 1M events but typically less. (each production has a different cross-section and therefore all the event do not have the same weight, if you combine the events correctly you have to “loose” some event such that you can have all the events with the same weight).
This explains why making 1M production of 1 events do not work.
Obviously this is something which is only important if you generate a lot of small sample of events.

> Is my interpretation correct?

Not really, we still keep the (few) events generated for that channel when we decide what to keep/discard and they participate to the final unweighting. They typically have a small weight and are therefore typically discarded by the unweighting but it can happen that they are kept such that the sum of multiple sample should have that contribution present. Now the integration error is quite large on those channels and the weight associated is therefore not precise. But this effect should be averaged between the multiple production (some time the weight will be too high and sometimes too low) such that the contribution in the sum of the sample should occur at the correct rate.

> Is it correct also to state that this problem does not exist in case of
> gridpacks?

In gridpack, all the contribution are computed with high precision during the creation of the gridpack. Then for a given number of requested event, we know how much event is required for each channel. if this number is smaller than 1 (you can change the threshold), then a random number is thrown to decide if we have to run that channel or not. So this is not the same strategy but you basically have the same effect.

Cheers,

Olivier

On 15 Oct 2015, at 13:03, pietro govoni <email address hidden> wrote:

> Question #272287 on MadGraph5_aMC@NLO changed:
> https://answers.launchpad.net/mg5amcnlo/+question/272287
>
> pietro govoni posted a new comment:
> Dear Olivier,
>
> thanks a lot for the detailed explanation.
> So it seems to be that having a gridpack is not only in general faster, but even safer.
> In fact, from this sentence of yours:
>
> "Asking more event, can also require to make more computation since
> small channel can then contribute."
>
> is seems to be that generating a small set of events might lead to neglecting some channels, and therefore that a produciton of 1M events might give a different outcome than lumping together 1000 productions of 1K events.
> Is my interpretation correct?
>
> Is it correct also to state that this problem does not exist in case of
> gridpacks?
>
> Thanks a lot!!
>
> pietro
>
> --
> You received this question notification because you are an answer
> contact for MadGraph5_aMC@NLO.

Thanks Olivier, everything is clear.

Cheers,

pietro