Phase space efficiency, precision and speed ticket

[Juergen Reuter]

The setup for VBF does not integrate with the same result as in WHIZARD 1, Helac or Sherpa. The process is e1, E1 => e1, E1, u, U, d, D, the result from WHIZARD 1 is 0.6521(55) fb. The WHIZARD 2 result however is 0.866(35) fb where the integration is stuck and does not improve at all:

| Initializating integration for process xxxxx:
| Phase space file 'xxxxx.phs' not found.
| Generating phase space configuration ...
| ... done.
| ... found 795 phase space channels, collected in 55 groves.
| Phase space: found 1051 equivalences between channels.
| Wrote phase-space configuration file 'xxxxx.phs'.
| iterations = 5:100000, 25:200000
| Applying user-defined cuts.
| Creating VAMP integration grids:
| Using phase-space channel equivalences.
| 100000 calls, 795 channels, 14 dimensions, 12 bins, stratified = T
| Integrating process 'arsch':
|=============================================================================|
| It      Calls  Integral[fb]  Error[fb]   Err[%]    Acc  Eff[%]   Chi2 N[It] |
|=============================================================================|
   1     100000  5.7038854E-01  9.58E-02   16.80   53.13*   1.21
   2     100000  5.4088309E-01  5.46E-02   10.09   31.90*   0.78
   3     100000  6.9953138E-01  6.56E-02    9.38   29.66*   0.24
   4     100000  9.0424999E-01  1.75E-01   19.33   61.12    0.17
^[[1;5C   5     100000  8.3137084E-01  1.26E-01   15.17   47.97*   0.05
|-----------------------------------------------------------------------------|
   5     500000  6.3176528E-01  3.60E-02    5.69   40.26    0.05    2.30   5
|-----------------------------------------------------------------------------|
   6     200000  8.3751280E-01  1.54E-01   18.34   82.00    0.02
   7     200000  8.7231558E-01  1.45E-01   16.60   74.25*   0.02
   8     200000  8.2172036E-01  4.00E-02    4.87   21.78*   0.01
   9     200000  8.4897070E-01  3.68E-02    4.33   19.36*   0.01
  10     200000  8.6302096E-01  3.54E-02    4.10   18.32*   0.01
  11     200000  8.6603706E-01  3.53E-02    4.08   18.23*   0.02
  12     200000  8.6603852E-01  3.53E-02    4.08   18.23*   0.03
  13     200000  8.6604215E-01  3.53E-02    4.08   18.23*   0.02
  14     200000  8.6608253E-01  3.53E-02    4.08   18.23*   0.00
  15     200000  8.6612555E-01  3.53E-02    4.08   18.23*   0.01
  16     200000  8.6612643E-01  3.53E-02    4.08   18.23*   0.04
  17     200000  8.6612643E-01  3.53E-02    4.08   18.23*   0.02
  18     200000  8.6612645E-01  3.53E-02    4.08   18.23*    NaN
  19     200000  8.6612647E-01  3.53E-02    4.08   18.23*    NaN

I attach the corresponding SINDARIN input file as well as the WHIZARD 1 input and output files.

[Juergen Reuter]

W1 output

[Juergen Reuter]

Process definition

[Juergen Reuter]

WHIZARD 1 cut file

[Juergen Reuter]

WHIZARD 1 input file

[Juergen Reuter]

WHIZARD 2 SINDARIN file

[Juergen Reuter]

It looks like that this is indeed a wrong cut set-up, I did an angle cut between jet-jet which should not be there, and forgot the angle cut between lepton-lepton. Now it looks much better.

[Juergen Reuter]

There are loads of numbers completely off. Either the phase space setup of WHIZARD 2 is BY FAR too sensitive to small fluctuations from different channels or there are some problems with the cuts not properly transported from one process to the other or there are still some problems lingering around with collinear singularities. It seems that problems are inherent always then if there are electrons in the final state for electron-positron initial state. I'm stymied.

[Juergen Reuter]

The very problem seems to be the population of channels. It seems that according whether WHIZARD hits or misses a channel the outcome of the result is more or less arbitrary. So in principle for a 2->6 process we should always at least demand 15 iterations or even more (???).

[Juergen Reuter]

Started test runs for the same process ten times. It really looks like the result is more or less arbitrary. Waiting for the final results.

[Juergen Reuter]

I ran the process e1, E1 => e1, E1, u, U, d, D on 10 different cores. Here, I do post the results, the SINDARIN input file and the corresponding output file with the vamp history. (I do post one complete file in order for you to get all the info.)

[Juergen Reuter]

complete output file

[Juergen Reuter]

I also attach the complete log file for run 10 as this run is deviating a lot.

[Juergen Reuter]

complete output file run 10

[Juergen Reuter]

Run 10

[Juergen Reuter]

Run 09

[Juergen Reuter]

Run 08

[Juergen Reuter]

Run 07

[Juergen Reuter]

Run 06

[Juergen Reuter]

Run 05

[Juergen Reuter]

Run 04

[Juergen Reuter]

Run 03

[Juergen Reuter]

Run 02

[Juergen Reuter]

Run 01

[Juergen Reuter]

When you look into the run files, then the dominating channels/groves are the following: 38, 51, sometimes 54 (multi-peripheral with two t-channel props.), 21 with one t-channel electron or neutrino prop., sometimes 16 (Higgsstrahlung with H -> WW -> jjjj). Very rarely the groves 46,47 (VBF -> Higgs -> WW -> jjjj and VBF -> Higgs -> ZZ -> jjjj) are populated. Would a larger weight on this channel yield the Sherpa/Helac? results? Will test when switching off the Higgs exchange.

[Juergen Reuter]

As a prime example to test the speed issue for the generation of the phase space as well as the generation of the adaption and phase space integration, take the following process:

model = NMSSM_Hgg  
 
alias tau_b = e3:E3:b:B 
alias proton = u:d:s:c:U:D:S:C:g  
alias fjet = proton:tau_b

# 4 taus + XX
# process uu_2x_4t_nmssm  = u:U,u:U => fjet, fjet, e3, E3, e3, E3
# process dd_2x_4t_nmssm  = d:D,d:D => fjet, fjet, e3, E3, e3, E3
# process cc_2x_4t_nmssm  = c:C,c:C => fjet, fjet, e3, E3, e3, E3
# process ss_2x_4t_nmssm  = s:S,s:S => fjet, fjet, e3, E3, e3, E3
# process gg_2x_4t_nmssm  = g,g => fjet, fjet, e3, E3, e3, E3
process pp_4t = proton, proton => fjet, fjet, e3, E3, e3, E3
process pp_4b = proton, proton => fjet, fjet, b, B, b, B
process pp_2t2b = proton, proton => fjet, fjet, e3, E3, b, B

compile

?slha_read_decays = true
read_slha("nmssm.slha")

cuts = all Pt >= 10 GeV [fjet]
       and all -5 <= Eta <= 5 [fjet]
       and all Dist > 0.4 [fjet,fjet]

sqrts = 14 TeV 
beams = p,p => lhapdf 


# integrate (uu_2x_4t_nmssm,dd_2x_4t_nmssm,cc_2x_4t_nmssm,ss_2x_4t_nmssm,gg_2x_4t_nmssm) { iterations = 5:352080,5:600000 }
integrate (pp_4t,pp_4b,pp_2t2b) { iterations = 5:50000,15:100000 }

# simulate (uu_2x_4t_nmssm,dd_2x_4t_nmssm,cc_2x_4t_nmssm,ss_2x_4t_nmssm,gg_2x_4t_nmssm) { luminosity = 100 }
simulate (pp_4t,pp_4b,pp_2t2b) { luminosity = 100 }

Code generation takes 1-2 hours, compilation is negligible, while the generation of the phase space files takes oder 12 hours, and integration is really horribly slow.

[kilian]

Investigating this, finally.

Something's definitely fishy. Look at the first channel in the generate phs file. It looks like

 ! Channel #            1
   tree     3   12   15   48   63
   map   3 s_channel      23 ! Z
   map  12 s_channel      23 ! Z
   map  15 s_channel      25 ! H
   map  48 s_channel      23 ! Z
   map  63 s_channel      23 ! Z

This is kinematically impossible, they can't all be resonant simultaneously. I'll check this first.

[kilian]

The problem with the first process ee -> ee udUD is solved.

With a tiny modification in the phase-space code that eliminates the impossible resonances, and with the following parameters

   phs_off_shell = 2
   phs_t_channel = 3

I get a stable and reliable result (15x100k + 5x100k):

  20     500000  6.8627427E-01  4.42E-03    0.64    4.55    1.07    0.02   5

The efficiency is above 1%, which is quite good. Stability and efficiency are better than with W1. BTW, the central value 0.6863(44) is also slightly higher and thus in line with PHEGAS and SHERPA.

[Juergen Reuter]

An NMSSM process I tried ended with a core dump, possibly because of too much memory. Here is the corresponding SINDARIN file:

model = NMSSM_Hgg
alias tau_b = e3:E3:b:B
alias proton = u:d:s:c:U:D:S:C:g
alias fjet = proton:tau_b
process pp_4t = proton, proton => fjet, fjet, e3, E3, e3, E3
process pp_4b = proton, proton => fjet, fjet, b, B, b, B
process pp_2t2b = proton, proton => fjet, fjet, e3, E3, b, B
compile
?slha_read_decays = true
read_slha("nmssm.slha")
cuts = all Pt >= 10 GeV [fjet]
       and all -5 <= Eta <= 5 [fjet]
       and all Dist > 0.4 [fjet,fjet]
sqrts = 14 TeV
beams = p,p => lhapdf
{ iterations = 5:352080,5:600000 }
integrate (pp_4t,pp_4b,pp_2t2b) { iterations = 5:50000,15:100000 }

This should be monitored.

[kilian]

Some more comments about the process ee -> ee udUD (all with the modified phs code as mentioned above):

First, important subprocesses are Higgs-strahlung, ZZ->H, Ae -> WWe, and AA -> WW.

After several test runs, there is a systematic pattern:

• If I leave phs_off_shell as default (=1), I get a typical result 0.653(2).
• If I set phs_off_shell=2, I get typically 0.677(3) or more (see above: 0.686(4)).
• Modifying phase space manually, i.e., adding important offshell=2 channels, I can get intermediate results.

I consistently use t-channel=3 now, because of the AA -> WW subprocess.

The lower result appears to be wrong (compare HELAC, Sherpa) but is stable. The error estimate apparently doesn't recognize the missing channels and corresponding contributions.

The obvious solution is to have off=2 by default, at least for processes like this. Disadvantage is: More channels require more statistics to begin with, so low-statistic runs (which make sense for more favorable processes) may be less stable.

I have still to find some heuristics to determine when off=2 is useful and when not. Maybe for 2->6 and more as default?

[kilian]

Finally decided to set the defaults to phs_off_shell=2 (see above) and threshold_calls=10 (no channels are dropped during adaptation) in r3220. This appears to cover the most complicated cases, at the expense of some efficiency loss for simple cases.