Better pruning of events from Pythia8 and enable in HepMC (#12275)

Fixes the pruning of events in GeneratorPythia8. The old implementation had a number of problems

- Beam particles (HepMC status = 4) was not kept
- The event structure was not preserved.

That meant that the event would be split into many disjoint trees. Thus, subsequent passes over the event would not be able to tie things together. This was especially a problem when converting the event to HepMC (via AOD) since the events would not be complete.

Added the option to prune events read in via GeneratorHepMC. The configuration key HepMC.prune=true will enable this mode (off by default).

By pruning, I mean that all particles that are not

- beam particles (HepMC status = 4),
- decayed (HepMC status = 2), nor
- final state (HepMC status = 1)

are removed from the event. When we do that, it is extremely important to connect up particles so that we have a consistent event structure. For example, we may remove a gluon from the tree. In that case, we must take care to move all child particles to the proper node in the tree, or the children will be disjoint from the rest of the event.

This patch also adds Doxygen documentation to GeneratorPythia8, as well as a new run/SimExamples/Pythia8/README.md file.

Author: cholmcc

Generators/include/Generators/GeneratorHepMC.h

@@ -31,6 +31,7 @@ namespace HepMC3
 class Reader;
 class GenEvent;
 class FourVector;
+class GenParticle;
 } // namespace HepMC3
 #endif
 
@@ -101,12 +102,22 @@ class GeneratorHepMC : public Generator, public GeneratorFileOrCmd
   /** Make our reader */
   bool makeReader();
 
+  /** Type of function to select particles to keep when pruning
+   * events */
+  typedef bool (*Select)(std::shared_ptr<const HepMC3::GenParticle>);
+  /** Prune event of particles that are not selected by passed
+   * function.  The event structure is preserved. */
+  void pruneEvent(Select select);
+
   /** HepMC interface **/
   uint64_t mEventsToSkip = 0;
+  /** HepMC event record version to expected.  Deprecated. */
   int mVersion = 0;
   std::shared_ptr<HepMC3::Reader> mReader;
   /** Event structure */
   HepMC3::GenEvent* mEvent = nullptr;
+  /** Option whether to prune event */
+  bool mPrune; //!
 
   ClassDefOverride(GeneratorHepMC, 1);
 

Generators/include/Generators/GeneratorHepMCParam.h

@@ -30,9 +30,27 @@ namespace eventgen
  **/
 
 struct GeneratorHepMCParam : public o2::conf::ConfigurableParamHelper<GeneratorHepMCParam> {
+  /** Version number of event structure to decode.  Note, when reading
+   *  from a file, this key is ignored.  The interface will figure out
+   *  the version automatically.  When reading from a pipe, and the
+   *  command line writes HepMC version 2 format, then this parameter
+   *  should be set to 2 */
   int version = 0;
+  /** Number of events to skip at the start of the input */
   uint64_t eventsToSkip = 0;
+  /** Deprecated.  Set the name of the file to read from.  Use
+   * GeneratorFileOrCmd.fileNames instead. */
   std::string fileName = "";
+  /** Wether to prune the event tree.  If true, then only particles that are
+   *
+   * - beam particles (status == 4)
+   * - decayed (status == 2), or
+   * - final state (status == 1)
+   *
+   * are kept.  This reduces the event size.  How much depends on the
+   * event generator producing the event.  Use with caution, as it may
+   * corrupt the event record. */
+  bool prune = false;
   O2ParamDef(GeneratorHepMCParam, "HepMC");
 };
 

Generators/include/Generators/GeneratorPythia8.h

@@ -24,7 +24,63 @@ namespace eventgen
 
 /*****************************************************************/
 /*****************************************************************/
-
+/** Interface to the Pythia8 event generator.  This generator is
+ * configured by configuration files (e.g.,
+ * `Generators/share/egconfig/pythia8_inel.cfg` for the type of
+ * events to generate.
+ *
+ * The above file, for example, contains
+ *
+ * @verbatim
+ * ### Beams, proton-proton collisions at sqrt(s)=14TeV
+ * Beams:idA 2212                  # proton
+ * Beams:idB 2212                  # proton
+ * Beams:eCM 14000.                # GeV
+ *
+ * ### Processes, min-bias inelastic events
+ * SoftQCD:inelastic on            # all inelastic processes
+ *
+ * ### Decays. Only decay until 1cm/c - corresponding to (physical) primaries
+ * ParticleDecays:limitTau0 on
+ * ParticleDecays:tau0Max 10.
+ * @endverbatim
+ *
+ * The configuration file to read is initially set in
+ * `GeneratorFactory`, but an additional configuration file can be
+ * specified with the configuration key `GeneratorPythia8::config`.
+ *
+ * If the configuration key `GeneratorPythia8::includePartonEvent` is
+ * set to false (default), then the event is pruned.  That is, all
+ * particles that are not
+ *
+ * - beam particles (HepMC status = 4),
+ * - decayed particles (HepMC status = 2), nor
+ * - final state partcles (HepMC status = 1)
+ *
+ * are removed from the event before passing on to the simulation.
+ * The event structure is kept, so that we have a well-formed event
+ * structure.  This reduces the event size by roughly 30%.
+ *
+ * If the configuration key `GeneratorPythia8::includePartonEvent` is
+ * true, then the full event is kept, including intermediate partons
+ * such as gluons, pomerons, diffractive hadrons, and so on.
+ *
+ * In the future, the way to prune events may become more flexible.
+ * For example, we could have the configuration keys
+ *
+ * - GeneratorPythia8::onlyStatus a list of HepMC status codes to accept
+ * - GeneratorPythia8::onlyPDGs a list of PDG particle codes to accept
+ *
+ * The configuration key `GeneratorPythia8::hooksFileName` allows the
+ * definition of a Pythia8 user hook.  See for example
+ * `Generators/share/egconfig/pythia8_userhooks_charm.C`.  The file
+ * specified is interpreted via ROOT (i.e., a ROOT script), and the
+ * function name set via the configuration key
+ * `GeneratorPythia8::hooksFuncName` (default `pythia8_user_hooks`) is
+ * executed.  That function must return a pointer to a
+ * `Pythia8::UserHooks` object (see the Pythia8 manual for more on
+ * this).
+ */
 class GeneratorPythia8 : public Generator
 {
 
@@ -36,47 +92,69 @@ class GeneratorPythia8 : public Generator
   /** destructor **/
   ~GeneratorPythia8() override = default;
 
+  /** @{
+      @name methods to override **/
   /** Initialize the generator if needed **/
   Bool_t Init() override;
-
-  /** methods to override **/
+  /** Generate a single event */
   Bool_t generateEvent() override;
+  /** Import particles from Pythia onto the simulation event stack */
   Bool_t importParticles() override { return importParticles(mPythia.event); };
+  /** @} */
 
-  /** setters **/
+  /** @{
+   * @name setters **/
+  /** Set Pythia8 configuration file to read */
   void setConfig(std::string val) { mConfig = val; };
+  /**
+   * Set the ROOT script file name that defines a Pythia8::UserHooks
+   * object */
   void setHooksFileName(std::string val) { mHooksFileName = val; };
+  /** Function in ROOT script that returns a pointer to a
+   * Pythia8::UserHooks object */
   void setHooksFuncName(std::string val) { mHooksFuncName = val; };
-  void setUserHooks(Pythia8::UserHooks* hooks)
-  {
-#if PYTHIA_VERSION_INTEGER < 8300
-    mPythia.setUserHooksPtr(hooks);
-#else
-    mPythia.setUserHooksPtr(std::shared_ptr<Pythia8::UserHooks>(hooks));
-#endif
-  }
-
-  /** methods **/
+  /** Set the user hooks (defined in a Pythia8::UserHooks object) for
+   * the event generator. */
+  void setUserHooks(Pythia8::UserHooks* hooks);
+  /** @} */
+
+  /** @{
+   * @name Configuration methods **/
+  /** Read a Pythia8 configuration string */
   bool readString(std::string val) { return mPythia.readString(val, true); };
+  /** Read a Pythia8 configuration file */
   bool readFile(std::string val) { return mPythia.readFile(val, true); };
+  /** @} */
 
-  /** utilities **/
+  /** @{
+   * @name Utilities **/
+  /** Get number of binary collisions.  Note that this method deviates
+   * from how the Pythia authors count number of binary collisions */
   void getNcoll(int& nColl)
   {
     getNcoll(mPythia.info, nColl);
   };
+  /** Get number of participants.  Note that this method deviates
+   * from how the Pythia authors count number of participants */
   void getNpart(int& nPart)
   {
     getNpart(mPythia.info, nPart);
   };
+  /** Get number of participants, split by nucleon type and origin.
+   * Note that this method deviates from how the Pythia authors count
+   * number of participants */
   void getNpart(int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg)
   {
     getNpart(mPythia.info, nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg);
   };
+  /** Get number of nuclei remnants, split by nucleon type and origin. */
   void getNremn(int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg)
   {
     getNremn(mPythia.event, nProtonProj, nNeutronProj, nProtonTarg, nNeutronTarg);
   };
+  /** Get number of spectators, split by nucleon type and origin.
+   * Note that this method deviates from how the Pythia authors count
+   * number of spectators */
   void getNfreeSpec(int& nFreenProj, int& nFreepProj, int& nFreenTarg, int& nFreepTarg)
   {
     getNfreeSpec(mPythia.info, nFreenProj, nFreepProj, nFreenTarg, nFreepTarg);
@@ -88,27 +166,101 @@ class GeneratorPythia8 : public Generator
   /** operator= **/
   GeneratorPythia8& operator=(const GeneratorPythia8&);
 
-  /** methods that can be overridded **/
-  void updateHeader(o2::dataformats::MCEventHeader* eventHeader) override;
+  /** @{
+   * @name Some function definitions
+   */
+  /** Select particles when pruning event */
+  typedef bool (*Select)(const Pythia8::Particle& particle);
+  /** Get relatives (mothers or daughters) of a particle */
+  typedef std::vector<int> (*GetRelatives)(const Pythia8::Particle&);
+  /** Set relatives (mothers or daughters) of a particle */
+  typedef void (*SetRelatives)(Pythia8::Particle&, int, int);
+  /** Get range of relatives (mothers or daughters) of a particle */
+  typedef std::pair<int, int> (*FirstLastRelative)(const Pythia8::Particle&);
+  /** @} */
 
-  /** internal methods **/
-  Bool_t importParticles(Pythia8::Event const& event);
+  /** @{
+   * @name Methods that can be overridded **/
+  /** Update the event header.  This propagates all sorts of
+   * information from Pythia8 to the simulation event header,
+   * including parton distribution function parameters, event weight,
+   * cross-section information, heavy-ion collision information, and
+   * so on. */
+  void updateHeader(o2::dataformats::MCEventHeader* eventHeader) override;
+  /** @} */
 
-  /** utilities **/
+  /** @{
+   * @name Internal methods **/
+  /** Import particles from Pythia onto the simulation stack */
+  Bool_t importParticles(Pythia8::Event& event);
+  /** Prune an event.  Only particles for which the function select
+   * returns true are kept in the event record.  The structure of the
+   * event is preserved, meaning that particles will point back to
+   * their ultimate (selected) mothers and, if select preserves the
+   * beam particles, the ultimate collision interaction. */
+  void pruneEvent(Pythia8::Event& event, Select select);
+  /** Investigate relatives (mothers or daughters) for particles to
+   * keep when pruning an event.  This checks the current particle,
+   * identified by index, if any of its relatives (either mothers or
+   * daughters) are to be kept.  If a relative is to be kept, then
+   * that is added to an internal list.  If a relative is _not_ to be
+   * kept, then that relatives relatives are queried (recursive call).
+   * The result of the recursive call is a list of relatives to the
+   * current particle which are to be kept.  These are then also added
+   * to the internal list.  The relatives that are found to be kept
+   * are then set to be relatives of the current particle.  Note that
+   * this member function modifies the relatives of the passed
+   * particle, and thus modifies the passed event structure.
+   * Calculations are cached. */
+  void investigateRelatives(Pythia8::Event& event,           // Event
+                            const std::vector<int>& old2New, // Map from old to new idx
+                            size_t index,                    // Current particle
+                            std::vector<bool>& done,         // cache flag
+                            GetRelatives getter,             // get relatives
+                            SetRelatives setter,             // set relatives
+                            FirstLastRelative firstLast,     // get first and last relative
+                            const std::string& what,         // what are we looking for
+                            const std::string& ind = "");    // logging indent
+  /** @{
+   * @name utilities **/
+  /** Select from ancestor. Fills the output event with all particles
+   * related to an ancestor of the input event
+   *
+   * Starting from ancestor, select all daughters (and their daughters
+   * recursively), and store them in the output event.
+   **/
   void selectFromAncestor(int ancestor, Pythia8::Event& inputEvent, Pythia8::Event& outputEvent);
+  /** Get number of binary collisions.  Note that this method deviates
+   * from how the Pythia authors count number of binary collisions */
   void getNcoll(const Pythia8::Info& info, int& nColl);
+  /** Get number of participants.  Note that this method deviates
+   * from how the Pythia authors count number of participants */
   void getNpart(const Pythia8::Info& info, int& nPart);
+  /** Get number of participants, split by nucleon type and origin.
+   * Note that this method deviates from how the Pythia authors count
+   * number of participants */
   void getNpart(const Pythia8::Info& info, int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg);
+  /** Get number of nuclei remnants, split by nucleon type and origin. */
   void getNremn(const Pythia8::Event& event, int& nProtonProj, int& nNeutronProj, int& nProtonTarg, int& nNeutronTarg);
+  /** Get number of spectators, split by nucleon type and origin.
+   * Note that this method deviates from how the Pythia authors count
+   * number of spectators */
   void getNfreeSpec(const Pythia8::Info& info, int& nFreenProj, int& nFreepProj, int& nFreenTarg, int& nFreepTarg);
+  /** @} */
 
   /** Pythia8 **/
   Pythia8::Pythia mPythia; //!
 
-  /** configuration **/
+  /** @{
+   * @name Configurations */
+  /** configuration file to read **/
   std::string mConfig;
+  /** ROOT script defining a Pythia8::UserHooks object */
   std::string mHooksFileName;
+  /** Function in `mHooksFileName` to execute to return pointer to
+   * Pythia8::UserHooks object */
   std::string mHooksFuncName;
+  /** @} */
 
   ClassDefOverride(GeneratorPythia8, 1);
 
@@ -121,3 +273,6 @@ class GeneratorPythia8 : public Generator
 } // namespace o2
 
 #endif /* ALICEO2_EVENTGEN_GENERATORPYTHIA8_H_ */
+//
+// EOF
+//