Unified treatment of beam particles

[jncots]

Fixes #101.

The patch prepend to event 2 beam particles if the generator doesn't contain them, or refill 2 first entry if generator doesn't produce correct one (case of nuclei for DpmJet). These 2 beam particles are calculated from EventKinematics and take into account reference frames of generator and EventKinematics to provide correct px, py, pz, en values. 2 beam particles have status == 4.

The event entries for beam particles are calculated only once for EventKinematics object and then applied for each event generated for this EventKinematics. Therefore, the beam particles event entries are kept in EventKinematics as one of the attribute.

For the particles which have mothers = [-1, -1] and some other similar combinations ([0, -1], [1, -1], [1, 1]), the mothers are reset to [0,1], so all such particles are daughters of beam particles. Phojet193 has mothers of the kind [6, 4]. I don't know what it means, but when such event are stored, pyhepmc produces error. I swapped number to [4, 6] to fix this. Svg graphs after these fixes have a single origin vertex where beam particles meet.

[jncots]

BeamData class contains beam particles in the format of EventData. Each event kinematics contains an instance of BeamData ready to apply to the event produced by generator. It is done because of efficiency considerations, otherwise one should create a new BeamData for each new event. However, it is still needed to create a new BeamData in case of CompositeTarget.

Initially implementation used EventData, but there is a circular dependence + there is a need for ref_frame attribute. However, if there is an idea how to use EventData instead of BeamData I will change it back.

Some tests were changed to reflect changes in MCEvent class.

[jncots]

It seems that all tests pass. It is a part of #165 issue. Together with #167 it is a part of release #164.

[jncots]

@HDembinski, thanks for comments, the implementation became much simpler. Now it is consist of one additional method for EventKinematics: get_beam_data and 2 methods in EventData. I needed an inverse boost, so I added 1 parameter for apply_boost and used SimpleNamespace to have en, pz attributes.

[jncots]

There have been introduced couple of changes.

• Because models are more different that it was assumed, it is more convenient to introduce _history_zero_indexing and _repair_initial_beam for MCEvent and override them for each model separately.
• Sophia history is added
• Phojet and Dpmjet have a complex history that hepmc cannot understand, therefore to_hepmc3 was overriden to output only final particles attached to beam particles. The original history (mothers and daughters) is left unchanged except shift to -1 to follow zero-indexing.
• Epos-LHC history changes when target is nuclei. It is only possible to add target ahead of the stack to be valid history for hepmc. Not sure if this is correct.
• Test for writing and reading hepmc format is added to tests correctness of event history.
• Test for beam particles is added. Epos-LHC is checked only for pp because pO switches beam particles.

[jncots]

Following the discussion in #164

• PhoJet and DpmJet history for hepmc has been changed. Now it is done similar of what @HDembinski did for Pythia8
• new version of pyhepmc was introduced. To make sure that only pyhepmc>=2.13.2 is used, Runtime error raised if the version of pyhepmc < 2.13.2. Let me know if it should be done more gently.
• main was merged into this branch

[jncots]

I noticed 2 following problems with Phojet and Dpmjet.

• In Phojet all decayed particles disappear from history, i.e. only gammas appear and not pi0, and so on. It is a problem of Phojet, because the printout of original stack has the same problem
• In DpmJet only final particles are boosted to lab frame. I found the place where it could be changed:

C in SUBROUTINE DT_LT2LAB

C This

      DO i = NPOint(4) , NHKk
         IF ( (ABS(ISThkk(i)).EQ.1) .OR. (ISThkk(i).EQ.1000) .OR. 
     &        (ISThkk(i).EQ.1001) ) THEN
            CALL DT_LTNUC(PHKk(3,i),PHKk(4,i),pz,pe,-3)
            PHKk(3,i) = pz
            PHKk(4,i) = pe
         END IF
      END DO

C should be changed to

      DO i = 1 , NHKk
            CALL DT_LTNUC(PHKk(3,i),PHKk(4,i),pz,pe,-3)
            PHKk(3,i) = pz
            PHKk(4,i) = pe
      END DO

@afedynitch, could you change it?

[afedynitch]

Cool! Did you try the proposed source code? It may need to be behind a flag, maybe. Need to have a deeper look

[jncots]

Cool! Did you try the proposed source code? It may need to be behind a flag, maybe. Need to have a deeper look

Yes, I tried it. It works good.

I also tried to fix the energy using only chromo:

self._lib.dtflg1.iframe = 2
self._frame = EventFrame.CENTER_OF_MASS

but I've seen events where the energy of decayed particle (rho+ if I remember it correctly) is larger than the sum of its products.

[jncots]

The issues noticed by @afedynitch have been resolved.

Problems with Phojet (removed history) and Dpmjet (Lorentz boost only for final particles) should be resolved.

One strage behaviour of EposLHC is noticed in connection with test_final_state: it produce different particles from run to run for the same seed. Because of that it is not possible to define known_issues to pass the test.

[afedynitch]

EposLHC is noticed in connection with test_final_state: it produce different particles from run to run for the same seed.

Does it mean that it doesn't conserve the random number sequence? It would be strange because it runs in CORSIKA 7.

[HDembinski]

Does it mean that it doesn't conserve the random number sequence? It would be strange because it runs in CORSIKA 7.

I worked on the rng instrumentation in EPOS-LHC. It is confusing because they use several seeds. Maybe some of them are not properly initialized. I believe it is rather a bug in our RNG handling than in EPOS-LHC itself.