Reenable shower and Pythia

[kilian]

said in the (new) title

[sschmidt]

Status of the Parton Shower Implementation: My implementation so far is in branches/schmidtboschmann/firstFSR. In whizard.nw I implemented an additional Fortran file shower_interactions.f90 with a subroutine event_apply_shower_particle_set that extracts the partons that can undergo showering from event%particle_set, transfers these partons to my shower program (located in src/miisr) and transfers the results back to the particle_set of the event. It should be working, although I can't guarantee it. The next steps will be:

• Up to now all parameters are hardcoded. The shower should use given parameters, probably given in the sin file. Moreover the switches whether the showers are to be executed at all are missing.
• The interface to PYTHIA is to be implemented.

[sschmidt]

About the implementation of the PYTHIA interface: There are as I see it, three ways of doing this, and I'm not sure which on to prefer:

• Implementation via LHEF

It could be implemented by letting W2 write events out in Les Houche Event Files, then letting PYTHIA read these back in, shower (and hadronize) them and transfer them back to W2. The advantage is obviously that this is according to accepted standards, the disadvantage being that the user needs to have PYTHIA >= 6.404 (Some experimentalists might still use PYTHIA 6.2xx?) and that every event needs to be written to a file which might make the whole process quite slow.

• Implementation like pythia_interface.f90 in WHIZARD 1

In W1 the transfer was done by transfering the event to a /HEPEVT/ COMMON BLOCK, then importing this COMMON BLOCK into PYTHIA by a PYTHIA routine (PYHEPC). Transfering the results back ist just the same steps in reverse order. The advantage is that this is according to standard(s) too, while the disadvantage is the appearance of the otherwise unused COMMON BLOCK /HEPEVT/.

• Implementation like my shower_to_pythia.f90

For my standalone PS program I implemented an interface to PYTHIA that writes directly into the PYTHIA COMMON BLOCKS. I used it with the PYTHIA versions 6.228 and 6.409 without problems, although as it is no standard incompatibilities might arise with different versions of PYTHIA.

Any preferences?

[dwiesler]

My 2 cents on this last point:

The first implementation is the one I realized using F77 code, as in my branch dwiesler/mlm.

The second and third point are quite similar in principle. In some Freiburg discussion I also suggested to go w/ point 2, for which I could then rewrite the recent version of my interface using the pythia_interface.f90 of W1 as a starting point. The thread safety regarding the usage of common blocks can then be restored using a dummy interface, which is used if PYTHIA is not linked.

Regarding the version of PYTHIA: I think it's fair to require at most V6.41x (dated end of 2007), since other tools are required to be far more uptodate (gfortran, nagfor, ocaml etc...)

[Christian Speckner]

To me, version 2 seems to make most sense. Using an external file if the data can also be passed internally seems a bit awkward, and if PYTHIA and its common blocks are hidden behind an interface which communicates with WHIZARD via derived types, 2 can be turned into 1 easily if preferred at some later point. Also, thread safety could be achieved without resorting to an external program, e.g. by forking and using pipes and signals to communicate between the two different processes.

[kilian]

Hmm.

I'm actually inclined towards version 1. It's not so much the COMMON blocks, it is the complication of linking PYTHIA directly, which opens the same can of worms in the configure/make step as WHIZARD 1 suffered from. And we are literally linked to PYTHIA, while in version 1 PYTHIA can trivially be replaced by any program that reads LHEF, including potentially WHIZARD itself.

The other complication is the fact that the interface requires WHIZARD (as an external process) to be called as a subprogram from PYTHIA as main program, while WHIZARD again should have the overall control. In W1, this was awkward, and I would not want to implement it again. In version 1, PYTHIA, HERWIG, or whatever would be called standalone (fork) or as a subprogram, and we just have to respect initialization/calling conventions.

Writing external files has no performance penalty on current OS, since everything is held in memory anyway and written to cache/disk only occasionally. The Fortran compilers know about this, and we may even use status=scratch files that leave no trace in the file system. I use this mechanism for the phase-space file; the 'cascades' module constructs phase space and transfers it to the main program not as an argument, but as a scratch file. I also use scratch files for MD5 sum calculation.

[Juergen Reuter]

Well, there is one reason which might favour version 2. Namely event files can become really, really large, especially for showered LHC events. So reading in and out files might become problematic.

[sschmidt]

Replying to kilian:

The other complication is the fact that the interface requires WHIZARD (as an external process) to be called as a subprogram from PYTHIA as main program, while WHIZARD again should have the overall control. In W1, this was awkward, and I would not want to implement it again. In version 1, PYTHIA, HERWIG, or whatever would be called standalone (fork) or as a subprogram, and we just have to respect initialization/calling conventions.

If you are referring to implementing WHIZARD events as external processes for PYTHIA, I agree.

But it should be possible to do the reverse. IMHO the shower subroutine using PYTHIA could look something like this

... definitions of common blocks etc...

call pyinit(...)                               ! initialize PYTHIA
call pyhepc(1)                                 ! transfer data from /HEPEVT/ fo PYTHIA
call pyevt()  or call pysspa(); call pyshow()  ! perform event (here only showering)
call pyhepc(2)                                  ! transfer data back

so that PYTHIA could be called by WHIZARD. That should be the way the interface was realized in W1 in pythia_interface.f90.

[kilian]

Forgot about that interface. IIRC, the external-process interface was supposed to be the preferred mode of interaction with PYTHIA. If the simpler version still works, we could implement it.

Still, configuring is complicated by the need for linking PYTHIA. It used to require parts of the CERNLIB and even PDFLIB, both long obsolete. What are the needs of the more recent PYTHIAs?

[sschmidt]

According to Daniel: PYTHIA, at least the version Daniel uses, doesn't require anything to be linked.

[sschmidt]

In my latest commit (r2554) i committed my implementation of solution #1 (lhe files). Originally I favoured solution #2 (HEPEVT) but it turned out that starting with 6.40? PYTHIA reads in lhef when invoking PYINIT, so for not using lhef one would have to replace PYTHIA's UPINIT and UPEVNT routines, which would have been messy.

The next steps/problems are:

1. Rewrite to use scratch files (maybe keep external files as a means of debugging)
   
2. Decide how to handle beam-remnants. As lhe is intended for parton-level events, beam-remnants are not considered in the first place. For now I just keep them as ordinary particles, which breaks W2 when reading them back in as they are no SM particles.
   
3. Implement a SINDARIN-interface, to read in parameter settings from the .sin file.
   
4. Implement a Hadronization interface? As the PYTHIA interface is there, Hadronization should be just some switches...

Can't think of anything more right now. If anybody has an opinion about these points...

[sschmidt]

In Reply to the next steps mentioned above:

1. done
   
2. Still don't know how to handle beam-remnants. When beam remnants are read back in, they are labelled to be gluons as the kT-codes used for them are not present in the models.
   
3. mostly done, the problem is in an extra ticket -> #340
   
4. Didn't implement a hadronization interface as this leads to the same problem as point 2. To read the results of the hadronization back in, you would need a model file with all the hadrons in it because otherwise W2 complains that it doesn't know the particle.

[sschmidt]

A comment on the status of the shower:

1) the interface to PYTHIA works with the exception of beam remnants mentioned in the last post.
2) my shower "works", it produces output. Comparing the results of the FSR with data (e+e- => qqbar at LEP) shows good agreement. The agreement is not perfect, but the comparison was made using a hadronisation that was fitted to PYTHIA, so there should be no bigger problem here. For the ISR, the results (pT-distribution of Z in Drell-Yan) isn't as good. Recently I discovered that I use the Sudakov factor in a slightly different interpretation with respect to PYTHIA, I don't know if the differences stem from this difference.

[sschmidt]

An update:

1) the hadronization interface to PYTHIA is implemented in my branch.
2) a handling of beam-remnants is implemented, although so far it's only a minimal version just to be able to interface the hadronization and it only works for one hard interaction so far.
3) A seperate model file SM_hadrons.mdl is implemented to be able to read the hadronization results back in.
4) I found a mistake in one of my subroutines and comitted a corrected version in r2917. The problems with the pT-spectrum of Z-Bosons at Tevatron seem to have been resolved, I made a quick comparison with some Tevatron data, ​http://james.physik.uni-freiburg.de/~sschmidt/Zptcomparison.png. However the problem remains, that there too few events at a high transverse momentum (pT > 30 GeV) are generated.

Currently the algorithm is much too slow, I will try some rearrangement to make it faster.

[Juergen Reuter]

Ok, I leave a reminder here for the reweighting of the MLM matched samples in the HepMC format. So S.S. can comment this as 'done' here when it is done.

[Juergen Reuter]

Isn't it time to close that ticket and maybe open up more specific tickets for open problems connected to that topic!?

[Juergen Reuter]

That is somehow the definition of 2.1.0.

[Juergen Reuter]

After the experience from tonight we should add a shower test to the WHIZARD self-tests before officially releasing the shower which is a light-weight version of one of the examples w/ and w\ matching.

[Juergen Reuter]

Replying to jr_reuter:

After the experience from tonight we should add a shower test to the WHIZARD self-tests before officially releasing the shower which is a light-weight version of one of the examples w/ and w\ matching.

Ok, doing so at the moment.

[Juergen Reuter]

I will not enable the test for the ISR at the moment, because it is hard-coded to use LHAPDF. We should have an option to use puilt-in PDFs here. Postpone this until 2.0.8. or 2.1.0.

[Juergen Reuter]

The release will be done in the trunk in commit r3846. Documentation and more validation is still pending.

[Juergen Reuter]

We change the scope of this ticket to reenabling the parton shower again after restructuring the code from 2.1 -> 2.2. But it clearly comes after #468. The documentation part of the ticket silently goes to the general manual and documentation ticket.

[kilian]

likely tbd after decays/polarization

[Juergen Reuter]

First steps in r4960: Shower_interface is compiling again. Minimal-invasive:

1. Removed strfun and event_formats dep., added hep_common dep.
2. commented out old STRF_LHAPDF, STRF_PDF_BUILTIN flags.

[Juergen Reuter]

Next steps: 1) reinstate rt_data options for shower/hadronization/matching 2) making shower_settings routines type-bound 3) steer shower settings through dispatch function (incl. unit test)

[Juergen Reuter]

The steps 1)-3) have been mostly done. I will look into this for now, and then eventually pass back to WK for event_transform sugar.

[kilian]

As suggested by JR, I'll take this over for a while, trying to adapt the event transform workflow where necessary.

[kilian]

I don't change the severity, but priority is after current bug fixes

[kilian]

Trying to do this before the Würzburg meeting.

[kilian]

The shower module is now part of the event-transform facility and works, at least superficially. Done in r5326.

I don't close the ticket yet because there are a few issues remaining:

• The PYTHIA shower appears to have a problem with e+ e- events without explicit beams. Maybe this was already in 2.1. I don't include a test for this variant, yet.
• I have to find out the exact meaning of the PDF integer IDs that the shower routine wants. Probably the ISR shower depends on it, also not yet tested.
• Should also include tests for matching etc.
• Where does the shower get its random numbers?

[Juergen Reuter]

In the src/shower, there is a routine randomseed which is however never used, using the system clock to initialize TAO random seed. Otherwise there are many different tao_random_number calls in shower_modul and shower_parton_module.

[Juergen Reuter]

So shouldn't the shower use the rng_tao_factory as well?

[Juergen Reuter]

shower_basics_module.f90 contains partially of auxiliary functions (like Altarelli-Parisi splitting functions), and partially initializers for global variables for the shower setup, shouldn't this be cleaned up?

[kilian]

I think the shower code needs a thorough revision regarding the points that we just raised, but I don't see this happen for the current release. Presently, we have to make sure that all allowed modes do something sensible, include a basic test for each of those, and make sure that system dependency doesn't break the test suite. We *may* be able to achieve that in WÜ.

[Juergen Reuter]

As discussed in Würzburg, WK please first continue (a bit) with this one here ...

[kilian]

That's all I wanted to do for now (r5363). I added a test for the ISR shower (parton_shower_2) which uses the PDF assignment.

Unfortunately, this test takes forever for a single event (actually, a couple of seconds). Maybe this is initialization, I didn't check. It is not possible to stop the program there, we must include signal_is_pending tests and return breakpoints in the shower routine, wherever it hangs. Furthermore, numerics is unstable (at low level, but still), so the output is not checked.

JR, can you take over? Ideally, we should have

1. discover the bottleneck
2. signal catching
3. more checks, in particular PYTHIA shower
4. validation of the results.

[Juergen Reuter]

Will take over. #612 is part of it, in principle.

[Juergen Reuter]

As of r5496, there is also a complete overhaul of the shower code: I didn't change any contents, and made sure step by step that the shower events are the same as before, but the whole shower code including the dummy is now in a single noweb file. This is included in the documented source code in share/doc and might serve as an example how things work when splitting up the monster .nw (#499). I even transferred shower routines into type-bound procedures. Though this doesn't change anything, it might help BC for further development. Pythia showering and hadronization do work (though I might have to check WHIZARD internal analysis using SM_hadrons). There are two more tests for the Pythia shower now. All further verification and profling are left either to the profiling ticket or ticket to be newly created. Some small technical issues have special tickets (#555 and #612). I will also go and check the MLM matching next (#545). Closing (rather burying forever) this one ...