Different particle counts between W1 and W2

[Juergen Reuter]

When looking at four jets (in that case coming from WW) then the LC generator group reported that the number of charged particles, photons and neutral particles in general is a lot higher in W2 then in W1. At the WHIZARD meeting we came up with a list of possible errors: Currently our best bet is actually 7) and would ask you to rerun Whizard2 with the suggested setting.

Possible issues:

1) Transferral of color flows from LHE to color ordering in PYJETS

2) Cuts

3) The exact partonic process: flavors, cross sections, parameters (especially

CKM matrix elements)

4) splitting into a showering / hadronization step

5) way of PYTHIA call (PY4FERM)

6) transferal of PYTHIA settings

7) Different scales in W1 vs W2 related to 5)

================================================================================ We observe:

1)

• Unlikely as the LHE Pythia6 interface is standard and used in countless runs (also by Madgraph).

2)

• W2: M > 10 GeV on all Jets. Looks simple enough to be the fine but we don't have the W1 input file to verify.

3)

• CKM Matrix elements in W2 have been updated and they are not changed from their defaults. Flavor definition of the process is weird but seems consistent between W1 and W2.

4)

• Unlikely as it works for 2jets.

5)

• We do not fully understand what is happening in py4ferm. It looks as if it is setting the Ws as mothers which could change the color reconnection that is done. See MSTP(115). In our setup with LHE there would be no color reconnection done for the 4 fermion process because there is no W as a mother. This is only an issue if the used tune actually uses MSTP(115). The setting of the probabilities in Tim's code according to the color flows looks fine if there is no bug hidden.

6)

• Unlikely as it works for 2jets.

7)

• Py4frm sets the initial showering scale to mW. This would also be a reasonable scale for W2

================================================================================ Things still to be done:

1)

• Verify that the color flow in the Whizard particle set is the same as in the Pythia common block with --debug2 shower.

2)

• Individual cross sections of each of the 5 components have to agree at the percent level

3)

• All SM parameters have to be checked between DBD samples and W2. Easily accessible in sindarin with show(model).

5)

• Verify in W1 that MSTP(115) is off and not activated by another setting (just form looking at the tune it seems to be off).

7)

• In W2 the scale should be set to mW with:

scale = mW

Per default the scale is sqrts. This could actually easily explain a lot more radiated particles in W2.

[Bijan Chokoufe Nejad]

Can you please add the sindarin file, including the commands to produce the plots, to reproduce this?

From Akiya's mail:

The whizard1.95 code used for the DBD production is available at ​https://svnweb.cern.ch/trac/lcgentools/browser/trunk/Whizard_1.95

They can be downloaded by svn co ​http://svn.cern.ch/guest/lcgentools/trunk/Whizard_1.95

As long as I see, PY4FRM is called at ​https://svnweb.cern.ch/trac/lcgentools/browser/trunk/Whizard_1.95/a6f/include/ilc_fragment_call.f90#L182

[Juergen Reuter]

I added the LC 1.95 version as branches/1.95_LC to our repo, starting to fix issues to get this running again.

[Juergen Reuter]

SINDARIN file to produce the distributions for charged particles and photons.

[Bijan Chokoufe Nejad]

Can you please add the sindarin file, including the commands to produce the plots, to reproduce this?

Another idea would be to do a factorized production such that pythia knows the parents of the jets similar to py4ferm.

[Bijan Chokoufe Nejad]

Apparently you still use the hacked model. For the record, to use the above attached sindarin one has to remove all lines involving pi0 (and thus also neutral). Still there is pythia6-parameters.sin missing. Do you have the file?

[Juergen Reuter]

Here is a file of the PYTHIA parameters truncated to the essentials.

[Juergen Reuter]

Replying to bchokoufe:

Apparently you still use the hacked model. For the record, to use the above attached sindarin one has to remove all lines involving pi0 (and thus also neutral). Still there is pythia6-parameters.sin missing. Do you have the file?

Partially yes, you should still use the 'hacked' model. 'neutral' should work with the actual trunk, however.

[Bijan Chokoufe Nejad]

AFAICT, We can cross off 1). I generated 5 events and the color flow in Whizard (c-sbar, ubar-d) was every time correctly recognized by Pythia. After shower there are usually only gluons in between but sometimes qqbar can break the color string (as expected).

[Bijan Chokoufe Nejad]

A run with 50k events can be found here:

​http://desy.de/~bcho/whizard_analysis.pdf

But I think only the photon distribution makes sense. The others likely didn't include the hadrons due to the model problem.

So an average of 50 photons per event is still a factor of 2 too many?

[Bijan Chokoufe Nejad]

For comparison, I can confirm that with scale = 500 GeV one gets only a bit more, cf. here

​http://desy.de/~bcho/whizard_analysis_500GeVscale.pdf

(53.8 instead of 50.4). So also 7) can be crossed off.

[Juergen Reuter]

There is the following line in shower.nw which looks suspicious: call pygive ("MSTP(111)=1") !!! Allow hadronization and decays as it seems to switch on hadronization even during the shower stage. However, the following test SINDARIN does not produce hadrons in any of the 4 events, so there hadronization seems to be switched off.

model = SM
process eeqq = e1, E1 => u, U
sqrts = 500 GeV
integrate (eeqq)
?ps_fsr_active = true
sample_format = lhef
n_events = 4
simulate (eeqq)

[Juergen Reuter]

I did the following debug procedure using this SINDARIN file:

model = SM
process eeqq = e1, E1 => u, U
sqrts = 500 GeV
integrate (eeqq)
?ps_fsr_active = true
$shower_method = "PYTHIA6"
?hadronization_active = true
sample_format = lhef
n_events = 4
simulate (eeqq)

with --rebuild-events --debug2 shower --debug2 transforms. What looks funny is, that for the second event onwards in the line 'before pyevnt, before boosting' there seems to be the old PYTHIA event still in the common blocks. Then, after all but the first event's statement 'After pyevnt, after boosting' there are hadronic particles in the event. So shouldn't we somehow clear the PYTHIA6 event record before calling pyevnt again?

[Juergen Reuter]

eeqq debug output, please look for the 'before/after pyevnt' statements

[Juergen Reuter]

Ok, the listings in debug.001 above look ok, all events are different, particle numbers also become smaller, and there seem to be no left-overs from previous events in the new events.

[Bijan Chokoufe Nejad]

Well yes it seems that 'before pyevnt' still has the old particle in the PYTHIA common blocks as pyevnt has not been called yet. But as you observe it should be no problem because otherwise particle number would never get smaller. But what I still find suspicous are the 'after pyevnt' listings. Only in the first event, the event is as expected (unhadronized). Afterwards it is not (and also not the old one). Maybe the settings are being overwritten?

[Bijan Chokoufe Nejad]

Replying to jr_reuter:

There is the following line in shower.nw which looks suspicious: call pygive ("MSTP(111)=1") !!! Allow hadronization and decays as it seems to switch on hadronization even during the shower stage.

The comment is not super clear. MSTP(111) is the master switch to allow hadronization and decays in principle. But in the two lines directly thereafter you find

       call pygive ("MSTJ(1)=0")       !!! No jet fragmentation
       call pygive ("MSTJ(21)=1")      !!! Allow decays but no jet fragmentation

[Bijan Chokoufe Nejad]

Replying to bchokoufe:

.. But what I still find suspicous are the 'after pyevnt' listings. Only in the first event, the event is as expected (unhadronized). Afterwards it is not (and also not the old one). Maybe the settings are being overwritten?

Well at least the settings stay fixed. A pystat(5) at the 'after pyevnt' shows that MSTP and PARP blocks are identical

[Juergen Reuter]

It really looks like that from the 2nd event on the hadronization would happen already at the stage of the shower, but then the complete event seems to stay unchanged through the hadronization event transform stage.

[Bijan Chokoufe Nejad]

Further possible issues:

8) Polarization. Everything is set to polarized in the given sindarin. Disable it and see if the number of photons changes. Pythia throws some warnings concerning polarization:

   I IST     ID Mothers   Colours    p_x      p_y      p_z       E        m
   1 -1      11   0   0    0    0    0.000    0.000  252.800  252.800    0.000
   2 -1     -11   0   0    0    0    1.341  -60.388 -325.710  325.710    0.000
   3  1       4   1   2  501    0  -24.574  -30.417 -199.219  203.021    0.000
   4  1      -3   1   2    0  501   21.452  -53.392 -109.592  123.779    0.000
   5  1      -2   1   2    0  502    1.917   25.519   12.852   28.636    0.000
   6  1       1   1   2  502    0    2.547   -2.097  223.050  223.074    0.000
Warning: Stored helicity information is wrong: 9for ip = 162.

          Event listing of user process at input (simplified)

   I IST     ID Mothers   Colours    p_x      p_y      p_z       E        m
   1 -1      11   0   0    0    0    0.000    0.000  254.898  254.898    0.000
   2 -1     -11   0   0    0    0  -22.142  -74.258 -306.423  306.423    0.000
   3  1      -4   1   2    0  502  -11.079  -39.319   74.752   85.186    0.000
   4  1      -1   1   2    0  501  -92.085  -11.238 -208.572  228.272    0.000
   5  1       2   1   2  501    0  -13.839   40.416  -54.645   69.362    0.000
   6  1       3   1   2  502    0   94.862  -64.118  136.940  178.501    0.000
|     90.0           45000               5000              0m:52s
Warning: Stored helicity information is wrong: 9for ip = 99.

[Juergen Reuter]

There is also the following:

9) TAUOLA and PHOTOS: both are switched on. On the other hand, the discrepancy I guess was there already before TAUOLA and PHOTOS have been added to WHIZARD2.

[Bijan Chokoufe Nejad]

When allowing hadronization in the shower step (by comment out "MSTJ(1)=0") and setting hadronization_active = true, I am still getting 52 photons on average. So 4) can likely be crossed off.

[Bijan Chokoufe Nejad]

Replying to jr_reuter:

It really looks like that from the 2nd event on the hadronization would happen already at the stage of the shower, but then the complete event seems to stay unchanged through the hadronization event transform stage.

Indeed and I think I found the culprit. Due to this

    if (shower%initialized_for_NPRUP >= NPRUP) then
      call msg_debug (D_SHOWER, "calling upinit")
      call upinit ()
    else
      ... pygive statements ...
       shower%initialized_for_NPRUP = NPRUP
    end if

mechanism the settings are only set once for each process. What also confuses me is why upinit should not be called for the first event if it provides necessary information.

[Bijan Chokoufe Nejad]

Enforcing to set

       call pygive ("MSTP(111)=1")     !!! Allow hadronization and decays
       call pygive ("MSTJ(1)=0")       !!! No jet fragmentation
       call pygive ("MSTJ(21)=1")      !!! Allow decays but no jet fragmentation

for each event still gives 50.4 photons... debug output looks more like expected though (after pyevnt only partons and only later on hadrons)

[Bijan Chokoufe Nejad]

I also checked that incoming and outgoing momenta also coiincide with the script tests/check-debug-output.py. Not a single mismatch in 5000 events. Not even if I reduce the numerical thresholds.

[Bijan Chokoufe Nejad]

By deactivating the polarized statements, I get far less warnings from Pythia but the average number of photons is still 50.5. So cross off 8) as well. Same when I deactivate the manual tauola activation (<#photons> = 50.0).

---

Mini summary: We checked 1), 4), 7) 8), 9)

Leaving open 2) Cuts, 3) exact partonic process, 5) PY4FERM, 6) pythia settings transferral

[Bijan Chokoufe Nejad]

I checked all that ALL of the settings that are set in pythia6-parameters.sin are correctly landing in Pythia's common blocks and stay the same over 5 events. So cross off 6) as well.

[Bijan Chokoufe Nejad]

Mhh. You get a very different photon distribution when you let Pythia decay the W's, see here ​http://desy.de/~bcho/whizard_analysis_WW.pdf Mean is 28.4 (so could be the factor of 2 we are looking for). I suspect that py4ferm ensures that the 4j case is treated like a WW. Next I try factorized decay in Whizard.

[kilian]

Replying to bchokoufe:

Mhh. You get a very different photon distribution when you let Pythia decay the W's, see here

That would be a disaster for any generic shower approach, right? The shower has to start from the specific EW process, then, i.e., matched to the matrix element...

[Bijan Chokoufe Nejad]

Replying to kilian:

Replying to bchokoufe:

Mhh. You get a very different photon distribution when you let Pythia decay the W's, see here

That would be a disaster for any generic shower approach, right? The shower has to start from the specific EW process, then, i.e., matched to the matrix element...

Well I don't think we are looking at a physical observable, when I only histogram photons with at least 1 GeV, the histogram changes again strongly

​http://desy.de/~bcho/whizard_analysis_WW_whizard_cut.pdf - mean is 21

(so this is factorized WW decay in Whizard and count if E > 1 Gev [A]). This is still sensitive to collinear splittings. For a proper analysis I think I have to build jets out of all particles in the event (currently not possible in sindarin due to the model issue). Or at least demand that photons are well separated from the rest (also not clear to me how I can formulate this).

Nevertheless also with the full process I now get with a 1 GeV cut

​http://desy.de/~bcho/whizard_analysis_cut.pdf - mean of 22 !!!

[Juergen Reuter]

Here are the comparisons of the cross sections for the five processes denoted

P1 = ww_h0cuxx: e1, E1 -> c:cbar, d:s:b:D:S:B, u:U, d:s:b:D:S:B
P2 = ww_h0uubd: e1, E1 -> u:U, b:B, u:U, d:D
P3 = ww_h0uusd: e1, E1 -> u:U, s:S, u:U, d:D
P4 = ww_h0ccbs: e1, E1 -> c:C, b:B, c:C, s:S
P5 = ww_h0ccds: e1, E1 -> c:C, d:D, c:C, s:S

The three numbers below are unpol., pol. LR (1./1.), pol. LR (1./1.) with ISR:

P1
WHIZARD1    1602.2(2.2) fb     6385.3(9.0) fb    6787.41(7.30) fb
WHIZARD2    1597.0(3.3) fb     6373.1(12.9)fb    6747.32(4.06) fb
----------------------------------------------------------------------
P2
WHIZARD1      -----------       0.0947(140) fb    0.101194(110) fb
WHIZARD2      -----------       0.0947(2)   fb    0.100463(60)  fb
-----------------------------------------------------------------------
P3
WHIZARD1      ------------      311.70(44) fb     332.46(36) fb
WHIZARD2      ------------      312.09(62) fb     330.83(19) fb
-----------------------------------------------------------------------
P4
WHIZARD1      ------------      10.739(147) fb     11.432(13) fb
WHIZARD2      ------------      10.731(22) fb      11.402(7) fb
------------------------------------------------------------------------
P5
WHIZARD1       ------------     310.86(79) fb      332.36(37) fb
WHIZARD2       ------------     312.33(60) fb      329.67(20) fb

The only cuts applied are the invariant mass cuts on the pairs of all jets of 10 GeV each. Over all, the numbers look fine, while for the case with ISR there is a slight undershooting in the case of WHIZARD2, which is borderline, and most likely comes from the different mapping and integration strategies between WHIZARD1 and WHIZARD2.

[Juergen Reuter]

Update from Akiya Miyamoto: the discrepancy between WHIZARD2 and WHIZARD 1.95LCGenG still persists when switching beamstrahlung and ISR off. The partonic cross sections here conincide at the fifth of a permil level, so no chance that tails of events produce more charged or neutral particles or that different mappings of ISR structure functions play a role. It is most likely the details in the treatment of PY4FRM based on HEPEVT events vs. the exchange via LHE events in WHIZARD2.

[Bijan Chokoufe Nejad]

Parts of Akiyas Mail with references:

According to the LEP-1 and LEP-2 measurements, the averaged charged multiplicity, for Z is about 20 and those for WW is about 40. See for example, ​http://pdg.lbl.gov/2015/reviews/rpp2015-rev-frag-functions.pdf Figure 20.6 and DELPHI paper, ​https://www.researchgate.net/deref/http%3A%2F%2Farxiv.org%2Fabs%2Fhep-ex%2F0103031v2 Eq. 36 and 37.

Copy from my Mail to Akiya:

I think I could track down the issue to whether Pythia knows that the 4j come from a WW or not. As I didn't manage to build ROOT in my environment, I rewrote the basics of your analysis in Rivet. Results are here (plots are normalized to number of events) [ATTACHED] (numberOfChargedTracks.pdf) scale is mW, no ISR, no beamstrahlung, sqrts = 500 GeV for both simulations.

For the $4j$ I just ran the sindarin with these options. For $WW\to 4j$, I let Whizard do a factorized WW decay with

  process jjjj_ww = e1, E1 => "W+", "W-"
  process Wp_dec = "W+" => uq, down_type_q
  process Wm_dec = "W-" => cq, down_type_q
  unstable "W+" (Wp_dec)
  unstable "W-" (Wm_dec)

I think the difference is similar to "W1 noISR, noBS" to "W2 noISR, noBS" in your plots. However, you seem to get an average number of charged tracks in W2 of ~55 while I get 49.6 but maybe this is not significant. In the factorized case the mean is 39.5 so just on point with your references.

The question remains what the correct physics is. My guess is that Pythia will preserve the mass of the resonant Ws when it has the mother daughter relations that the jets come from the Ws. One could try to fudge a py4ferm by inserting the mother-daughter-relations. However, this is also not really correct as we can see from the diagramgs [ATTACHED] (jjjj_i1_diags.pdf) there are even for this flavor composition still quite a number of non-double-resonant diagrams. Maybe it is best to just simulate factorized WW when you want to reproduce WW results.

PS: To the Whizards: I think the best thing one could do is to write out the most probable resonance history a la Jezo / Nason. Most probable can be computed with Breit-Wigner-Mappings. I am not quite sure how we could use our infrastructure for the Born case.

Upshot: We hope that the generator group is happy that they can reproduce the W1 results with the factorized approach and will close this issue for now. A more rigorous approach to set resonance histories in Whizard in general for parton showers is very desirable as it changes particle multiplicities, jet counts, etc. more strongly than naively expected.

[Bijan Chokoufe Nejad]

Couldn't upload all diagrams but the others have the same topologies