How to speed up process and event generation in a complicated model?

Asked by Sara Vanini

Hi,

I'm studying a BSM signal, and I have listed the processes I'm interested in. Since one single generation with all the 6 processes took 1 month to generate the diagrams, I decided it was not possible to handle. therefore I've produced one LHE file for each of the 6 processes, for each of the mass point.

Nevertheless, to submit the request to the CMS production system I need to merge all the processes for each mass point, so I need to create ONE single LHE file with all the generated events from the 6 processes. Is this possible? could I edit the header, the cross sections, and paste all the events in one single file? is there a standard script available to do so?

Please let me know. Thanks a lot, and best regards

Sara

Question information

Language:
English Edit question
Status:
Solved
For:
MadGraph5_aMC@NLO Edit question
Assignee:
No assignee Edit question
Solved by:
Johan Alwall
Solved:
Last query:
Last reply:
Revision history for this message
Johan Alwall (johan-alwall) said :
#1

Hello Sara,

The main problem is that if you merge the LHE files, you probably want the events in correct proportion according to their relative cross section? So you'll need to make sure to generate the correct number of events for each process. If you have done that, and use different process labels for the different processes, you can just combine the event files (and make sure to combine also the <init> blocks).

It sounds strange that it would take quite that long to generate the diagrams. If you want some advice how to generate your processes in a more efficient way, feel free to send me an email on jalwall_at_ntu.edu.tw.

All the best,
Johan

Revision history for this message
Andrea Gozzelino (andrea-gozzelino) said :
#2

Hi Johan,

I'm working with Sara.
I provide below the process cards
import model_v4 typeIIIseesaw140
# Define multiparticle labels
define l+ = e+ mu+ ta+
define l- = e- mu- ta-
define vll = v1 v2 v3
# Specify process(es) to run
generate p p > tr+ tr0 > l- w+ w+ vll , w+ > l+ vll
#generate p p > tr+ tr0 > l+ w- w+ vll , w+ > l+ vll, w- > l- vll
#generate p p > tr+ tr0 > l+ l- z w+ , z > vll vll , w+ > l+ vll
#generate p p > tr+ tr0 > l+ l+ z w- , z > vll vll , w- > l- vll
#generate p p > tr+ tr0 > l+ l+ z w- , z > j j , w- > l- vll
#generate p p > tr+ tr0 > l+ l- z w+ , z > j j , w+ > l+ vll
# Output processes to MadEvent directory
output -f

We processed a single card with all 6 process, but it did NOT produce LHE file. It crashed after week of work.

Below I write the run card, too:
#*******************
# Running parameters
#*******************
#
#*********************************************************************
# Tag name for the run (one word) *
#*********************************************************************
  'seesaw' = run_tag ! name of the run
#*********************************************************************
# Run to generate the grid pack *
#*********************************************************************
  .true. = gridpack !True = setting up the grid pack
#*********************************************************************
# Number of events and rnd seed *
# Warning: Do not generate more than 100K event in a single run *
#*********************************************************************
  30000 = nevents ! Number of unweighted events requested
      0 = iseed ! rnd seed (0=assigned automatically=default))
#*********************************************************************
# Collider type and energy *
#*********************************************************************
        1 = lpp1 ! beam 1 type (0=NO PDF)
        1 = lpp2 ! beam 2 type (0=NO PDF)
     4000 = ebeam1 ! beam 1 energy in GeV
     4000 = ebeam2 ! beam 2 energy in GeV
#*********************************************************************
# Beam polarization from -100 (left-handed) to 100 (right-handed) *
#*********************************************************************
        0 = polbeam1 ! beam polarization for beam 1
        0 = polbeam2 ! beam polarization for beam 2
#*********************************************************************
# PDF CHOICE: this automatically fixes also alpha_s and its evol. *
#*********************************************************************
 'cteq6l1' = pdlabel ! PDF set
#*********************************************************************
# 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
#*********************************************************************
 1 = 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
#*********************************************************************
# Automatic ptj and mjj cuts if xqcut > 0
# (turn off for VBF and single top processes)
#**********************************************************
   T = auto_ptj_mjj
#**********************************************************
#
#**********************************
# BW cutoff (M+/-bwcutoff*Gamma)
#**********************************
# 10000 = bwcutoff
   15 = bwcutoff
#**********************************************************
# Apply pt/E/eta/dr/mij cuts on decay products or not
# (note that etmiss/ptll/ptheavy/sorted cuts always apply)
#**********************************************************
   T = cut_decays
#*************************************************************
# 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
#*******************
# Standard Cuts
#*******************
#
#*********************************************************************
# Minimum and maximum pt's *
#*********************************************************************
 20 = ptj ! minimum pt for the jets
 20 = 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
 1d5 = ptjmax ! maximum pt for the jets
 1d5 = ptbmax ! maximum pt for the b
 1d5 = ptamax ! maximum pt for the photons
 1d5 = ptlmax ! maximum pt for the charged leptons
 1d5 = missetmax ! maximum missing Et (sum of neutrino's momenta)
#*********************************************************************
# Minimum and maximum E's (in the lab 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
 1d5 = ejmax ! maximum E for the jets
 1d5 = ebmax ! maximum E for the b
 1d5 = eamax ! maximum E for the photons
 1d5 = elmax ! maximum E for the charged leptons
#*********************************************************************
# Maximum and minimum rapidity *
#*********************************************************************
   5 = etaj ! max rap for the jets
   5 = etab ! max rap for the b
 2d5 = etaa ! max rap for the photons
 2d5 = etal ! max rap for the charged leptons
 0.6 = etaonium ! max rap for the quarkonium states
 0d0 = etajmin ! min rap for the jets
 0d0 = etabmin ! min rap for the b
 0d0 = etaamin ! min rap for the photons
 0d0 = etalmin ! main rap for the charged leptons
#*********************************************************************
# Minimum and maximum DeltaR distance *
#*********************************************************************
 0.001 = drjj ! min distance between jets
 0.001 = drbb ! min distance between b's
 0 = drll ! min distance between leptons
 0 = draa ! min distance between gammas
 0.001 = drbj ! min distance between b and jet
 0 = draj ! min distance between gamma and jet
 0 = drjl ! min distance between jet and lepton
 0 = drab ! min distance between gamma and b
 0 = drbl ! min distance between b and lepton
 0 = dral ! min distance between gamma and lepton
 1d2 = drjjmax ! max distance between jets
 1d2 = drbbmax ! max distance between b's
 1d2 = drllmax ! max distance between leptons
 1d2 = draamax ! max distance between gammas
 1d2 = drbjmax ! max distance between b and jet
 1d2 = drajmax ! max distance between gamma and jet
 1d2 = drjlmax ! max distance between jet and lepton
 1d2 = drabmax ! max distance between gamma and b
 1d2 = drblmax ! max distance between b and lepton
 1d2 = dralmax ! maxdistance between gamma and lepton
#*********************************************************************
# Minimum and maximum invariant mass for pairs *
#*********************************************************************
 0 = 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
 1d5 = mmjjmax ! max invariant mass of a jet pair
 1d5 = mmbbmax ! max invariant mass of a b pair
 1d5 = mmaamax ! max invariant mass of gamma gamma pair
 1d5 = 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)
 1d5 = 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)
 1d5 = 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
 1d5 = ptj1max ! maximum pt for the leading jet in pt
 1d5 = ptj2max ! maximum pt for the second jet in pt
 1d5 = ptj3max ! maximum pt for the third jet in pt
 1d5 = ptj4max ! maximum pt for the fourth jet in pt
 0 = cutuse ! reject event if fails any (0) / all (1) jet pt cuts
#*********************************************************************
# Control the Ht(k)=Sum of k leading jets *
#*********************************************************************
 0 = htjmin ! minimum jet HT=Sum(jet pt)
 1d5 = htjmax ! maximum jet HT=Sum(jet pt)
 0 = ihtmin !inclusive Ht for all partons (including b)
 1d5 = 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
 1d5 = ht2max ! maximum Ht for the two leading jets
 1d5 = ht3max ! maximum Ht for the three leading jets
 1d5 = ht4max ! maximum Ht for the four leading jets
#*********************************************************************
# WBF cuts *
#*********************************************************************
 0 = xetamin ! minimum rapidity for two jets in the WBF case
 0 = deltaeta ! minimum rapidity for two jets in the WBF case
#*********************************************************************
# maximal pdg code for quark to be considered as a light jet *
# (otherwise b cuts are applied) *
#*********************************************************************
 5 = maxjetflavor
#*********************************************************************
# Jet measure cuts *
#*********************************************************************
 20 = xqcut ! minimum kt jet measure between partons (for P1 P11 P2 P3)
#10 = xqcut ! minimum kt jet measure between partons (for P4 P5)
#*********************************************************************

Have you some suggestions about the strategies?

Thanks,
Andrea

Revision history for this message
Olivier Mattelaer (olivier-mattelaer) said :
#3

Hi,

Could you send to me (and Johan) the model such that we can play a bit with it.
I suspect that this is so slow due to some inefficiency in the model.
My email is omatt AT illinois.edu

Cheers,

Olivier

Revision history for this message
Andrea Gozzelino (andrea-gozzelino) said :
#4

Hi Olivier and Johan,

I put all 5 models and cards (5 param, 1 run, 1 proc) in my personal web page.
The file is typeIIIseesaw.tar.gz (the last on the page)
link = http://www.pd.infn.it/~agozzeli/

The ideal situation will compact 6 process in 1 process cards and finally we will have 5 mass points * 2 signs = 10 process cards/10 LHE files.
As it is now, we will produce 60 LHE files.

Thanks,
Andrea

Revision history for this message
Best Johan Alwall (johan-alwall) said :
#5

Hello Andrea,

I would suggest the following changes in your models and process definitions, in order to get (very) reasonable generation times:

1. Make sure that
 a. all light quark and lepton masses are set to zero, including tau mass. This allows MadGraph to combine as many processes as possible into fewer process directories, and speeds up the generation by orders of magnitude.
 b. You should also remove the distinction between neutrino flavor and gauge eigenstates (just use the regular ve, vm, vt), to avoid a huge number of unnecessary process distinctions. The flavor mixing of neutrinos is certainly not relevant for collider physics (unless what you're doing is simulating neutrino flavor experiments).
 c. Same thing with the CKM matrix - just use unit matrix for CKM, and remove the irrelevant interactions (different quark flavors can anyway not be distinguished at a collider experiment).
 d. I would finally suggest to also remove the couplings between Higgs and light quarks/leptons, to reduce the number of (completely irrelevant) diagrams.

2. You should change your process formulation to

generate p p > tr+ tr0, (tr0 > l- w+, w+ > l+ vll), (tr+ > w+ vll , w+ > l+ vll) @0
add process p p > tr+ tr0, (tr0 > l+ w-, w- > l- vll), (tr+ > w+ vll , w+ > l+ vll) @1
etc.

This will give very fast process generation (of the order minutes or even seconds), without any loss of generality (since your tr+ and tr0 widths are small enough that the narrow width approximation should work extremely well). Also your event generation should be reduced to minutes.

All the best,
Johan

Revision history for this message
Andrea Gozzelino (andrea-gozzelino) said :
#6

Hi,

I write down the process card:
import model_v4 typeIIIseesaw140bis
# Define multiparticle labels
define l+ = e+ mu+ ta+
define l- = e- mu- ta-
define vl = v1 v2 v3
# Specify process(es) to run
generate p p > tr+ tr0, (tr0 > l- w+, w+ > l+ vl), (tr+ > w+ vl , w+ > l+ vl) @0
add process p p > tr+ tr0, (tr0 > l+ w-, w- > l- vl), (tr+ > w+ vl , w+ > l+ vl) @1
add process p p > tr+ tr0, (tr0 > l- w+, w+ > l+ vl), (tr+ > z l+ , z > vl vl) @2
add process p p > tr+ tr0, (tr0 > l+ w-, w- > l- vl), (tr+ > z l+ , z > vl vl) @3
add process p p > tr+ tr0, (tr0 > l+ w-, w- > l- vl), (tr+ > z l+ , z > j j) @4
add process p p > tr+ tr0, (tr0 > l- w+, w+ > l+ vl), (tr+ > z l+ , z > j j) @5
# Output processes to MadEvent directory
output -f

I'm not familiar with models and parameter card. I try to build a new models, following your instructions, [typeIIIseesaw140bis]: I'm NOT able to manage that.
I will ask to my theoretician colleagues to produce new param cards.

Another problem is related with gridpack. I produce gridpack, but I can NOT run it on CRAB. Are there some examples?

Thanks,
Andrea

Revision history for this message
Johan Alwall (johan-alwall) said :
#7

Hello again Andrea,

Please note that what needs to be changed is NOT the param_cards (that would be trivial) but the actual model files (particles.dat and interactions.dat). Also note that in general, you only want one model and then have different param_cards corresponding to different parameter points, to avoid the risk of mismatch between different model points.

If this is a model created in FeynRules, there are simple ways to introduce model restrictions of the type I'm suggesting without changing the actual model files (see the FR manual). In that case, if you output the model as a UFO model (MG5 format) instead of an MG4 model, you can also use MG5 restriction cards to further reduce parameters for different scenarios (see the MG5 release paper).

Regarding the gridpack problems, I recommend that you open a new question for that.

All the best,
Johan

Revision history for this message
Andrea Gozzelino (andrea-gozzelino) said :
#8

Hi Johan,

I use UFO model in MG5format, but I'm really not expert on the other cards.
I contact directly people who made the UFO model.

Regarding the gripack problems, I open question 216203.

Thanks for very useful info,
Andrea

Revision history for this message
Andrea Gozzelino (andrea-gozzelino) said :
#9

Hi Johan,

I followed your suggestions.
I studied and modified model (interactions and particles data; parameters, processes and run cards).
Events generation is working properly and faster than before.

The issue is solved.
Thank you very much for your help,
Andrea

Revision history for this message
Sara Vanini (vanini-sara) said :
#10

Thanks a lot, following the suggestions from Johan we speed up the generation process indeed, without drawbacks. You've been very helpful. Regards Sara

Revision history for this message
Sara Vanini (vanini-sara) said :
#11

Thanks Johan Alwall, that solved my question.