Fix examples and comments

README.rst

@@ -2,9 +2,9 @@
 BxDecay0 - C++ port of the legacy Decay0 FORTRAN library
 ============================================================================
 
-:authors: François Mauger, Vladimir Tretyak
-:date: 2021-12-10
-:copyright: Copyright (C) 2017-2021 the BxCppDev group
+:authors: François Mauger, Volodymyr Tretyak
+:date: 2022-04-04
+:copyright: Copyright (C) 2017-2022 the BxCppDev group
 
 .. image:: resources/images/logo/bxdecay0_logo.png
    :width: 20%
@@ -13,7 +13,7 @@ BxDecay0 - C++ port of the legacy Decay0 FORTRAN library
 The **BxDecay0** C++  library provides a set of  classes and functions
 for the random generation of  simulated nuclear decays. It consists in
 a C++ port of the original Decay0 Fortran program (also known as GENBB
-in  some other  context) written  and maintained  by Vladimir  Tretyak
+in  some other  context) written  and maintained  by Volodymyr  Tretyak
 (KINR_).  Decay0 was created to  address the Monte Carlo generation of
 nuclear decays  in the context  of double  beta decay and  dark matter
 experimental research.
@@ -40,7 +40,7 @@ History
 =======
 
 The first version of the Decay0  program was written in the early 90's
-by  Vladimir   Tretyak  and   collaborators,  using  the   Fortran  77
+by Volodymyr   Tretyak  and   collaborators,  using  the   Fortran  77
 programming language and  the CERNLIB library. It  has been maintained
 and improved up to now (2020).
 
@@ -61,7 +61,7 @@ Fortran  and  CERNLIB.   This  code  was integrated  in  2011  as  the
 
 This release of the BxDecay0 C++  library is extracted from the Bayeux
 ``genbb_help`` module with some changes from the 2017-03-01 release of
-the Decay0 program  by Vladimir Tretyak.  It has then  been updated to
+the Decay0 program  by Volodymyr Tretyak.  It has then  been updated to
 the 2018-12-05 and 2020-04-20 releases.
 
 BxDecay0 is  a standalone library  with very few  dependencies (mostly
@@ -182,6 +182,7 @@ GSL_, ROOT_ or whatever).
 - Release 1.0.13 (in preparation):
 
   - Fix G4String issue in Geant4 extension (issue #26)
+  - Fix bugs concerning the manipulation of particles' times
 
 .. _SuperNEMO: http://supernemo.org/
 
@@ -900,7 +901,7 @@ functionalities.
 Authors and contributors
 ========================
 
-* Vladimir  Tretyak  (KINR_,  Kiev  Institute  for  Nuclear  Research,
+* Volodymyr  Tretyak  (KINR_,  Kiev  Institute  for  Nuclear  Research,
   Ukraine) is the original author  and maintener of the legacy Fortran
   Decay0 package.
 * François Mauger (`LPC Caen`_,  Laboratoire de Physique Corpusculaire
@@ -916,7 +917,7 @@ Authors and contributors
 * Malak Hoballah  and Laurent Simard (IJCLab,  Orsay, France) (support
   for DBD generator with gA process).
 * Rastislav Dvornicky (Comenius  University, Bratislava, Slovakia) has
-  calculated specific phase space factors for the DBD gA process.
+  calculated specific phase space factors for some DBD gA process.
 * Pierre Lasorak  (University of  Sussex, Oxford, United  Kingdom) has
   added  simplified versions  of Po210,  Po218, Th230,  U234 and  U238
   alpha emitters.
@@ -948,7 +949,7 @@ of their respective authors and/or scientific communities.
 References
 ===========
 
-* Vladimir Tretyak, *DECAY0 event  generator for initial kinematics of
+* Volodymyr Tretyak, *DECAY0 event  generator for initial kinematics of
   particles in  alpha, beta  and double beta  decays*, talk_  given at
   Laboratori Nazionali del Gran Sasso, 17 March 2015 :
 * O.A.Ponkratenko, V.I.Tretyak, Yu.G.Zdesenko,  Event Generator DECAY4

bxdecay0/event.cc

@@ -217,9 +217,17 @@ namespace bxdecay0 {
     double pz = p * ctet;
     tdlev_    = tclev_;
     if (thlev_ > 0.) {
+      // Randomize the emission time of the particle:
       tdlev_ = tclev_ - thlev_ / std::log(2.) * std::log(prng_());
     }
     part.set_code(np_);
+
+    // Here we  cumulate the decay  time of the new  particle (tdlev_)
+    // and the generation  time of the previous particle  in the event
+    // (last_time). This  is different of  what is done in  the legacy
+    // Decay0 program  which stores the time  delay between successive
+    // particles in the event:
+
     part.set_time(last_time + tdlev_);
     part.set_momentum(px, py, pz);
     event_.add_particle(part);

bxdecay0/event.h

@@ -43,16 +43,18 @@ namespace bxdecay0 {
   /// - generator : an optional generator string which may identifies the event generator responsible
   ///               of the creation of the current event,
   /// - time: an optional decay time reference which can be set externally by some activiy model
-  ///         (default value: unset==NaN)
+  ///         (default value: unset==NaN, unit: second)
   /// - particles : a list of generated particles, ordered in time.
   ///
   /// Each particle has the following attributes:
   ///
   /// - code : the particle identifier code (using GEANT3 particle code for alpha, electron, positron,
   ///          gamma, neutron or proton),
-  /// - time : the time elapsed from the event time reference (expressed in second),
+  /// - time : the particle creation time relative to the event time reference (unit: second);
+  ///          this is a major change compared to the original DECAY0 program which stores
+  ///          the times elapsed between successive particles within the event.
   /// - momentum: an array of 3 double precision reals representing the coordinates (px, py, pz)
-  ///             of the particle's momentum (expressed in MeV/c)
+  ///             of the particle's momentum (unit: MeV/c)
   ///
   class event
   {
@@ -94,7 +96,7 @@ namespace bxdecay0 {
     /// Check the validity of the event
     bool is_valid() const;
 
-    /// Shift the time of the particles by a delay (unit: second), starting from a given rank
+    /// Shift the times of the particles by a delay (unit: second), starting from a given rank
     void shift_particles_time(double delta_time_, const int from_ = 0);
 
     /// Reset the event
@@ -136,16 +138,19 @@ namespace bxdecay0 {
   // time of the last particle in the event (or 0 if there is no particle yet),
   // shifted from tclev + tdlev, where tdlev is randomized from thlev (halflife)
   //
-  //    event
+  // \code
+  //                                           last particle
+  //    event                                     /
   //  reference     particle#0          particle#1         thlev
   //    time            :                   :         :---------------------->: 
   // -----+-------------*-------------------*---------:----------------*-------------------> time
   //   t0 :------------>:                   :-------->:                :\
-  //   t1 :-------------------------------->:  tclev     tdlev(random) : new particle
+  //   t1 :-------------------------------->:  tclev     tdlev(random) : new particle #2
   //                                        :------------------------->:                             
   //   t2 :----------------------------------------------------------->:
   //                        time of the new particle
   //
+  // \endcode
   void randomize_particle(i_random & prng_,
                           event & event_,
                           const particle_code np_,

examples/ex02/ex02.cxx

@@ -1,7 +1,7 @@
 /** ex02.cxx
  *
- * Copyright 2017 François Mauger <mauger@lpccaen.in2p3.fr>
- * Copyright 2017 Normandie Université
+ * Copyright 2017-2022 François Mauger <mauger@lpccaen.in2p3.fr>
+ * Copyright 2017-2022 Normandie Université
  *
  * This file is part of BxDecay0.
  *
@@ -171,9 +171,10 @@ int main()
       gendecay.set_time(evtime);
 
       // Debug dump:
-      if (debug)
+      if (debug) {
         gendecay.print(std::cerr, "DBD event:", "[debug] ");
-
+      }
+
       // Build a Hep MC event:
       static const double C_LIGHT_MM_PER_SEC = 3e11;
       std::shared_ptr<HepMC3::GenEvent> genEvtPtr
@@ -211,8 +212,8 @@ int main()
           break;
         }
         if (particle.has_time()) {
-          // Cumulative particle time:
-          part_time += particle.get_time();
+          // Extract the particle time relatice to some time reference:
+          part_time = particle.get_time();
         }
         genPartPtr->set_pid(pid);
         genPartPtr->set_momentum(HepMC3::FourVector(particle.get_px(),

examples/ex03/ex03.cxx

@@ -133,9 +133,10 @@ int main()
       gendecay.set_time(evtime);
 
       // Debug dump:
-      if (debug)
+      if (debug) {
         gendecay.print(std::cerr, "DBD event:", "[debug] ");
-
+      }
+
       // Store the event:
       gendecay.store(fout, store_flags);
 

examples/ex04/README.rst

@@ -18,7 +18,7 @@ events for various isotopes and decay modes.
 The  ``bxdecay0::genbbsub``  function  is  the low  level  engine  for
 generating  decay  events  from  BxDecay0. Its  usage  is  mainly  for
 experts.  Normal users  should use  the ``bxdecay0::decay0_generator``
-class   which    wraps   in    a   OOP    approach   the    magic   of
+class   which    wraps   in    a   OOP    approach   the   *magic*  of
 ``bxdecay0::genbbsub``  inherited from  the Decay0  Fortran code  (see
 example ``ex01``, ``ex02`` and ``ex03``).
 

examples/ex04/ex04.cxx

@@ -65,7 +65,7 @@ int main()
     // Isotope:
     std::string chnuclide = "Mo100";
 
-    // Daughter's energy level (ground state):
+    // Daughter's energy level (0 = ground state):
     int ilevel = 0;
 
     // DBD mode (neutrinoless):
@@ -129,8 +129,9 @@ int main()
       gendecay.set_time(evtime);
 
       // Debug dump:
-      if (debug)
+      if (debug) {
         gendecay.print(std::cerr, "DBD event:", "[debug] ");
+      }
 
       // Clear the event:
       gendecay.reset();

extensions/bxdecay0_g4/bxdecay0_g4/primary_generator_action.cc

@@ -702,7 +702,11 @@ namespace bxdecay0_g4 {
     if (IsTrace()) std::cerr << "[trace] bxdecay0_g4::PrimaryGeneratorAction::GeneratePrimaries: Entering..." << '\n';
     bxdecay0::event gendecay;
     _pimpl_->get_decay0().shoot(_pimpl_->get_prng(), gendecay);
-    gendecay.set_time(0.0);
+    double event_time = 0.0 * CLHEP::second;
+    // Force event reference time to be zero (even if the generator tells something else):
+    // if (gendecay.has_time()) {
+    //    event_time = gendecay.get_time() * CLHEP::second;
+    // }
     double particle_time = 0.0;
     const auto & particles = gendecay.get_particles();
     if (IsDebug()) std::cerr << "[debug] bxdecay0_g4::PrimaryGeneratorAction::GeneratePrimaries: Nb particles=" << particles.size() << '\n';
@@ -738,8 +742,12 @@ namespace bxdecay0_g4 {
       }
       // Time:
       if (particle.has_time()) {
-        particle_time += gendecay.get_time() * CLHEP::second;
+        particle_time = particle.get_time() * CLHEP::second;
+      } else {
+        particle_time = 0.0 * CLHEP::second;
       }
+      // Shift particle time by the event reference time:
+      particle_time += event_time;
       _particle_gun_->SetParticleTime(particle_time);
       // Momentum:
       G4ThreeVector momentum(particle.get_px() * CLHEP::MeV,