diff --git a/EvtGenBase/EvtParticle.hh b/EvtGenBase/EvtParticle.hh
index a98473d..24e6505 100644
--- a/EvtGenBase/EvtParticle.hh
+++ b/EvtGenBase/EvtParticle.hh
@@ -1,503 +1,503 @@
 
 /***********************************************************************
 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
 *                                                                      *
 * This file is part of EvtGen.                                         *
 *                                                                      *
 * EvtGen is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU General Public License as published by *
 * the Free Software Foundation, either version 3 of the License, or    *
 * (at your option) any later version.                                  *
 *                                                                      *
 * EvtGen is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
 * GNU General Public License for more details.                         *
 *                                                                      *
 * You should have received a copy of the GNU General Public License    *
 * along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
 ***********************************************************************/
 
 #ifndef EVTPARTICLE_HH
 #define EVTPARTICLE_HH
 
 #include "EvtGenBase/EvtId.hh"
 #include "EvtGenBase/EvtSpinDensity.hh"
 #include "EvtGenBase/EvtSpinType.hh"
 #include "EvtGenBase/EvtVector4R.hh"
 
 #include <assert.h>
 #include <map>
 #include <string>
 #include <vector>
 
 class EvtDiracSpinor;
 class EvtVector4C;
 class EvtTensor4C;
 class EvtStdHep;
 class EvtSecondary;
 class EvtRaritaSchwinger;
 
 const int MAX_DAUG = 100;
 const int MAX_LEVEL = 10;
 const int MAX_TRIES = 10000;
 
 class EvtParticle {
   public:
     /**
   * Default constructor.
   */
     EvtParticle();
 
     /**
   * Destructor.
   */
     virtual ~EvtParticle();
 
     /**
   * Returns polarization vector in the parents restframe.
   */
     virtual EvtVector4C epsParent( int i ) const;
 
     /**
   * Returns polarization vector in the particles own restframe.
   */
     virtual EvtVector4C eps( int i ) const;
 
     /**
   * Returns polarization vector in the parents restframe for a photon.
   */
     virtual EvtVector4C epsParentPhoton( int i ) const;
 
     /**
   * Returns polarization vector in the particles own restframe for a photon.
   */
     virtual EvtVector4C epsPhoton( int i ) const;
 
     /**
   * Returns Dirac spinor in the parents restframe for a Dirac particle.
   */
     virtual EvtDiracSpinor spParent( int ) const;
 
     /**
   * Returns Dirac spinor in the particles own restframe for a Dirac particle.
   */
     virtual EvtDiracSpinor sp( int ) const;
 
     /**
   * Returns Dirac spinor in the parents restframe for a Neutrino particle.
   */
     virtual EvtDiracSpinor spParentNeutrino() const;
 
     /**
   * Returns Dirac spinor in the particles own restframe for a
   * Neutrino particle.
   */
     virtual EvtDiracSpinor spNeutrino() const;
 
     /**
   * Returns tensor in the parents restframe for a spin 2 particle.
   */
     virtual EvtTensor4C epsTensorParent( int i ) const;
 
     /**
   * Returns tensor in the particles own restframe for a spin 2 particle.
   */
     virtual EvtTensor4C epsTensor( int i ) const;
 
     /**
    * Returns Rarita-Schwinger spinor in the parents restframe for a 
    * Rarita-Schwinger particle.
    */
     virtual EvtRaritaSchwinger spRSParent( int ) const;
 
     /**
    * Returns Rarita-Schwinger spinor in the particles own restframe for a 
    * Rarita-Schwinger particle.
    */
     virtual EvtRaritaSchwinger spRS( int ) const;
 
     /**
   * Initialiaze particle with id and 4momentum.
   */
     virtual void init( EvtId part_n, const EvtVector4R& p4 ) = 0;
 
     /**
   * Add another daughter to the particle
   */
     void addDaug( EvtParticle* node );
 
     /**
   * Decay particle
   */
     void decay();
 
     /** 
   * Delete a decay chain
   */
     void deleteTree();
     void deleteDaughters( bool keepChannel = false );
 
     /**
   * Should only be used internally.
   */
     void setChannel( int i );
 
     /**
   * Creates the daughters in the list of ids and 
   * adds them to the parent. Note that momentum
   * is left uninitialized, this is _only_ creation.
   */
     void makeDaughters( unsigned int ndaug, EvtId* id );
 
     /**
   * Creates the daughters in the list of ids and 
   * adds them to the parent. Note that momentum
   * is left uninitialized, this is _only_ creation.
   */
     void makeDaughters( unsigned int ndaug, std::vector<EvtId> idVector );
 
     /**
   * Similar to the routine above except that here 
   * momentum is generated according to phase space 
   * daughters are filled with this momentum.
   */
     double initializePhaseSpace( unsigned int numdaughter, EvtId* daughters,
                                  bool forceResetMasses = false,
                                  double poleSize = -1., int whichTwo1 = 0,
                                  int whichTwo2 = 1 );
 
     /**
   * Get pointer the the i:th daugther.
   */
     EvtParticle* getDaug( const int i ) { return _daug[i]; }
 
     /**
   * Get const pointer the the i:th daugther.
   */
     const EvtParticle* getDaug( const int i ) const { return _daug[i]; }
 
     /**
   * Iterates over the particles in a decay chain.
   */
     EvtParticle* nextIter( EvtParticle* rootOfTree = nullptr );
 
     /**
   * Makes stdhep list
   */
     void makeStdHep( EvtStdHep& stdhep, EvtSecondary& secondary,
                      EvtId* stable_parent_ihep );
     void makeStdHep( EvtStdHep& stdhep );
 
     /**
   * Gets 4vector in the labframe, i.e., the frame in which the root
   * particles momentum is measured.
   */
     EvtVector4R getP4Lab() const;
 
     /**
   * Gets 4vector in the labframe for the 4-momentum before FSR was 
   * generated in the parents decay. The lab frame is where the root
   * particles momentum is measured.
   */
-    EvtVector4R getP4LabBeforeFSR();
+    EvtVector4R getP4LabBeforeFSR() const;
 
     /**
   * Gets 4vector in the particles restframe, i.e. this functiont will
   * return (m,0,0,0)
   */
     EvtVector4R getP4Restframe() const;
 
     /**
   * Returns the 4position of the particle in the lab frame.
   */
     EvtVector4R get4Pos() const;
 
     /**
   * Returns pointer to parent particle.
   */
     EvtParticle* getParent() const;
 
     /**
   * Makes partptr the idaug:th daugther.
   */
     void insertDaugPtr( int idaug, EvtParticle* partptr )
     {
         _daug[idaug] = partptr;
         partptr->_parent = this;
     }
     /**
   * Returns mass of particle.
   */
     double mass() const;
 
     /**
   * Used internally to decide if first time particle is decayed.
   */
     int firstornot() const;
     void setFirstOrNot();
     void resetFirstOrNot();
 
     /**
   * Returns Id of particle.
   */
     EvtId getId() const;
 
     /**
   * Returns the PDG id of the particle
   */
 
     int getPDGId() const;
 
     /** 
   * Returns particle type.
   */
 
     EvtSpinType::spintype getSpinType() const;
 
     /**
   * Returns number of spin states of the particle.
   */
     int getSpinStates() const;
 
     /**
   * Returns 4momentum in parents restframe.
   */
     const EvtVector4R& getP4() const;
 
     /**
   * Sets the 4momentum in the parents restframe.
   */
     void setP4( const EvtVector4R& p4 )
     {
         _p = p4;
         _pBeforeFSR = p4;
     }
 
     void setP4WithFSR( const EvtVector4R& p4 ) { _p = p4; }
 
     void setFSRP4toZero() { _pBeforeFSR.set( 0.0, 0.0, 0.0, 0.0 ); }
 
     /**
   * Retunrs the decay channel.
   */
     int getChannel() const;
 
     /**
   * Returns number of daugthers.
   */
     size_t getNDaug() const;
     void resetNDaug()
     {
         _ndaug = 0;
         return;
     }
 
     /**
   * Prints out the particle "tree" of a given particle.  The
   * tree consists of all daughters (and their daughters, etc)
   * and their properties.
   */
     void printTree() const;
 
     void printTreeRec( unsigned int level ) const;
 
     std::string treeStr() const;
     std::string treeStrRec( unsigned int level ) const;
 
     /**
   * Prints information for the particle.
   */
     void printParticle() const;
 
     /**
   * Set lifetime of the particle in parents restframe.
   */
     void setLifetime( double tau );
 
     /**
   * Generate lifetime according to pure exponential.
   */
     void setLifetime();
 
     /**
   * Returns the lifetime.
   */
     double getLifetime() const;
 
     /** 
   * Set diagonal spindensity matrix.
   */
     void setDiagonalSpinDensity();
 
     /** 
   * Set spindensity matrix for e+e- -> V
   */
     void setVectorSpinDensity();
 
     /**
   * Set forward spin density matrix.
   */
     void setSpinDensityForward( const EvtSpinDensity& rho )
     {
         _rhoForward = rho;
     }
 
     /**
   * Set forward spin density matrix according to the density matrix
   * rho in the helicity amplitude basis.
   */
     void setSpinDensityForwardHelicityBasis( const EvtSpinDensity& rho );
     void setSpinDensityForwardHelicityBasis( const EvtSpinDensity& rho,
                                              double alpha, double beta,
                                              double gamma );
 
     /**
   * Returns a rotation matrix need to rotate the basis state
   * to the helicity basis. The EvtSpinDensity matrix is just use
   * as a matrix here. This function is to be implemented in each
   * derived class.
   */
     virtual EvtSpinDensity rotateToHelicityBasis() const = 0;
     virtual EvtSpinDensity rotateToHelicityBasis( double alpha, double beta,
                                                   double gamma ) const = 0;
 
     /**
   * Get forward spin density matrix.
   */
     EvtSpinDensity getSpinDensityForward() { return _rhoForward; }
 
     /**
   * Set backward spin density matrix.
   */
     void setSpinDensityBackward( const EvtSpinDensity& rho )
     {
         _rhoBackward = rho;
     }
 
     /**
   * Get backward spin density matrix.
   */
     EvtSpinDensity getSpinDensityBackward() { return _rhoBackward; }
 
     //Hacks will be removed when better solutions are thought of!
     //This is used to suppress use of random numbers when doing initialization
     //of some models.
     void noLifeTime() { _genlifetime = 0; }
 
     //lange - April 29, 2002
     void setId( EvtId id ) { _id = id; }
     void initDecay( bool useMinMass = false );
     bool generateMassTree();
 
     double compMassProb();
 
     //setMass will blow away any existing 4vector
     void setMass( double m ) { _p = EvtVector4R( m, 0.0, 0.0, 0.0 ); }
 
     //void setMixed() {_mix=true;}
     //void setUnMixed() {_mix=false;}
     //bool getMixed() {return _mix;}
 
     //void takeCConj() {EvtGenReport(EVTGEN_INFO,"EvtGen") << "should take conj\n";}
 
     //this means that the particle has gone through initDecay
     // and thus has a mass
     bool isInitialized() { return _isInit; }
     bool hasValidP4() { return _validP4; }
     bool isDecayed() { return _isDecayed; }
 
     // decay prob - only relevent if already decayed
     // and is a scalar particle
     // returned is a double* that should be prob/probMax
     // FIXME - this should probably be changed to std::optional
     const double* decayProb() const { return _decayProb; }
     void setDecayProb( double p );
 
     // Return the name of the particle (from the EvtId number)
     std::string getName() const;
 
     // Specify whether the particle has a special named attribute with
     // a set value. By default, nothing is set, but derived classes
     // can set this to mean something specific, e.g. if a photon is FSR
     void setAttribute( std::string attName, int attValue )
     {
         _intAttributes[attName] = attValue;
     }
 
     // Retrieve the integer value for the given attribute name
     int getAttribute( std::string attName ) const;
 
     // Specify if the particle has a double attribute value, e.g. amplitude weight.
     // By default, nothing is set, but derived classes can set this to mean something specific
     void setAttributeDouble( std::string attName, double attValue )
     {
         _dblAttributes[attName] = attValue;
     }
 
     // Retrieve the double value for the given attribute name
     double getAttributeDouble( std::string attName ) const;
 
   protected:
     void setp( double e, double px, double py, double pz )
     {
         _p.set( e, px, py, pz );
         _pBeforeFSR = _p;
     }
 
     void setp( const EvtVector4R& p4 )
     {
         _p = p4;
         _pBeforeFSR = _p;
     }
 
     void setpart_num( EvtId particle_number )
     {
         assert( _channel == -10 || _id.getId() == particle_number.getId() ||
                 _id.getId() == -1 );
         _id = particle_number;
     }
     bool _validP4;
 
     // A typedef to define the attribute (name, integer) map
     typedef std::map<std::string, int> EvtAttIntMap;
     EvtAttIntMap _intAttributes;
 
     // A typedef to define the attribute (name, double) map
     typedef std::map<std::string, double> EvtAttDblMap;
     EvtAttDblMap _dblAttributes;
 
   private:
     EvtParticle* _daug[MAX_DAUG];
     size_t _ndaug;
     EvtParticle* _parent;
     int _channel;
     int _first;
     EvtId _id;
     EvtVector4R _p;
     EvtVector4R _pBeforeFSR;
     double _t;
     bool _isInit;
     bool _isDecayed;
 
     //bool _mix;
 
     EvtSpinDensity _rhoForward;
     EvtSpinDensity _rhoBackward;
 
     void makeStdHepRec( int firstparent, int lastparent, EvtStdHep& stdhep,
                         EvtSecondary& secondary, EvtId* stable_parent_ihep );
     void makeStdHepRec( int firstparent, int lastparent, EvtStdHep& stdhep );
 
     //This is a hack until things gets straightened out. (Ryd)
     int _genlifetime;
 
     //should never be used, therefor is private.
     //these does _not_ have an implementation
     EvtParticle& operator=( const EvtParticle& p );
     EvtParticle( const EvtParticle& p );
 
     double* _decayProb;
 };
 
 #endif
diff --git a/History.md b/History.md
index bb5da84..34fb214 100644
--- a/History.md
+++ b/History.md
@@ -1,856 +1,859 @@
 # History file for EvtGen
 
 From version 2.0.0,
 Tabc labels refer to [Maniphest tasks](https://phab.hepforge.org/maniphest/query/nkBRd9OhPCBN/), while
 Dxyz labels refer to [Differential code reviews](https://phab.hepforge.org/differential/query/YDY8EgjNGd.e/) on HepForge:
 
 https://phab.hepforge.org/Tabc
 
 https://phab.hepforge.org/Dxyz
 
 ===
 ## R02-0X-00
 
+9 Apr 2024 Fernando Abudinen
+* D111: Added tests for decays with FSR and implemented const correctness for getP4LabBeforeFSR() function.
+
 19 Feb 2024 Fernando Abudinen
 * D108: Bugfix protect polarisation vector for vector and tensor particles against wrong indices
 
 19 Feb 2024 Fernando Abudinen
 * D107: Implement const correctness for epsPhoton and epsParentPhoton
 
 3 Feb 2024 John Back
 * T230: Add EvtBcVPPHad model for Bc to Jpsi p pbar pi decays and generalised particle ordering
   for EvtBcVHad decay files, courtesy of Aleksei Luchinsky, Dmitrii Pereima & Vanya Belyaev (LHCb).
   Updated EvtPhiDalitz model to use helicity amplitudes and fixed the indices used in
   EvtVector3R::dot(), courtesy of Arnau Brossa Gonzalo & Antonio Romero Vidal (LHCb).
 
 31 Oct 2023 Fernando Abudinen
 * D102: Bugfix probmax issue for X38722-+\_PSI\_GAMMA model
   - Calculation of amplitude moved to dedicated calcAmp function and fixed initialisation.
   - Introduced weighting to compensate for phase-space changes due to resonance width.
   - Fixed bugs in rho couplings and loops over photon and vector states.
 
 16 Oct 2023 Fernando Abudinen
 * D99: Add tests for PHOTOS
   - EvtParticle::getAttribute function made const
   - Added variables for testing of FSR to testing framework
   - Added dedicated tests for FSR with Photos
   - Turned off Photos for all other tests
 
 16 Oct 2023 Thomas Latham
 * D98: Modernise EvtIdSet and other improvements
   - Modernise and greatly simplify EvtIdSet implementation
   - Fixes in EvtPropSLPole to avoid unnecessary dynamic allocations
   - Other minor fixes and tidy-ups
 * Credit to Heather Ratcliffe and Chris Brady for providing
   and/or inspiring several of these improvements
 
 22 Aug 2023 Fernando Abudinen
 * D97: Bugfix probmax issue and introduced pole compensation for VTOSLL model
   - Calculation of amplitude moved to dedicated calcAmp function and fixed initialisation.
 
 22 Aug 2023 Tom Latham
 * CMake updates
   - Update default C++ standard to 17
   - Suppress 'up-to-date' messages during install
   - Add 'EvtGen' prefix to file names of custom CMake modules
   - Remove unused CMake module
 * Update CI config
   - Update default LCG version to 103
   - Add builds for el9 OS using LCG 104
   - Allow switching on/off building against each external
   - Test stage: attempt to improve selection of commits to diff in different cases
   - Update CVMFS tags as per https://cern.service-now.com/service-portal?id=outage&n=OTG0079356
 * Fix script for generating src dependencies
   - Adapted to behaviour of newer CMake version in LCG 103
 * Apply clang-format
   - Few small changes after update of CI clang version to 12
   - Document version of clang-format to be used
 * Updates to install script
   - Update to latest Pythia8 and HepMC3 versions
   - Fix to ensure dependencies are picked up from the local install
 
 22 Aug 2023 Andrii Verbytskyi
 * Patch for finding Pythia8 xmldoc path
 
 21 Aug 2023 Tom Latham
 * D96: Work around change of interface for setRndmEnginePtr in Pythia8 310
 
 22 Jun 2023 Ludovico Massaccesi
 * T219: Fixed amplitudes for D_DALITZ model of D0 -> pi+ pi- pi0.
 
 1 March 2023 Fernando Abudinen
 * D92: Bugfix probmax issue for TENSOR daughter in EvtSSD_DirectCP model.
        Calculation of amplitude in EvtSSD_DirectCP model moved to dedicated calcAmp function.
        Got rid of a few static variables along the way.
 
 8 Feb 2023 Alexei Sibidanov and Tom Latham
 * D90: Apply clang format and enable checking of formatting in gitlab CI
 
 3 Feb 2023 John Back
 * D91: Check for non-zero momentum for EvtSLBaryonAmp parent spin density matrix.
        Print out integrals of JSON test histograms.
 
 16 Dec 2022 John Back
 * D89: Added probabilities for B_c -> V pi+ and V pi+ pi0 EvtBcVHad modes.
 
 16 Dec 2022 Alexei Sibidanov
 * D88: Applied clang-tidy modernize-use-nullptr
 
 13 Dec 2022 Fernando Abudinen and Tom Latham
 * Various improvements to testing framework
 
 16 Nov 2022 Tom Latham
 * T123: Provide documention of how to contribute bug reports, feature requests, new/modified code
 
 8 Sep 2022 Fernando Abudinen, John Back, Michal Kreps, Tom Latham, Alex Ward
 * T108: Implement JSON test framework for all decay models
 
 9 June 2022 Michal Kreps
 * D84: Improve efficiency of RareLbToLll decay model for final states with e+e- pair.
 
 ===
 ## R02-02-00
 
 12 May 2022 Michal Kreps
 * D83: Fix double counting of charmonia with K0 decays for B0.
 
 11 May 2022 Tom Latham
 * D80: Add CMake options for enabling clang-tidy static analysis checks during build
 
 14 Apr 2022 John Back
 * D82: EvtDecayProb: initialise data members and remove empty destructor
 
 14 Apr 2022 Michal Kreps
 * D81: Derive EvtFlatSqDalitz from EvtDecayIncoherent since we directly provide
        final kinematics
 
 2nd Mar 2022 John Back
 * D78: Add Bc -> J/psi K+ pi- pi+ pi- pi+, Bc -> J/psi K+ K- pi+ pi- pi+ &
        Bc -> J/psi 4pi+ 3pi- decay modes to the BC_VHAD model, courtesy of
        Aleksei Luchinsky, Anna Danilina, Dmitrii Pereima & Vanya Belyaev (LHCb)
 
 ===
 ## R02-01-01
 
 8th Sep 2021 Michal Kreps
 * Update location of web page for Pythia8 download in setup script.
 
 8th Sep 2021 Markus Prim, Lu Cao, Chaoyi Lyu and Michel De Cian (Michal Kreps)
 * D73: Add new model for semileptonic B decays with BCL and BGL form-factors
 
 8th June 2021 Michal Kreps
 * T110, D71: Fix B+ --> eta' l nu BF which was order of magnitude too high.
   Balance the decrease by increasing B+ --> D0 l nu, which is after change
   still bit smaller than PDG 2021.
 
 8th Jun 2021 Michal Kreps
 * D71: Fix B+ --> eta' l nu BF.
 
 20th Apr 2021 Tom Lathem
 * D68: Fix compilation with Pythia 8.304
 
 17th Mar 2021 Michal Kreps
 * D62: Improve PI0DALITZ model to dynamically work out maximum probability to
   make it usuable also for eta --> llgamma decays.
   Remove ETA2MUMUGAMMA since it is a pure one-to-one copy of PI0DALITZ.
 
 15th Jan 2021 Michal Kreps
 * D47: Model to generate 4-body phase-space decays in restricted part of the m12-m34 space
 
 12th Jan 2021 Michal Kreps
 * D48: Fix bug in calculation of the daughter momentum in decay model EvtBsMuMuKK
 
 7th Jan 2021 Michal Kreps
 * D43: Allow to pass particle properties table in form of stringstream to constructor
   of EvtGen for use case where these are created on fly.
 
 10th Dec 2020 Michal Kreps
 * D36: EvtFlatSqDalitz model to be more efficient and to avoid cut-off around the edges
 
 21st Aug 2020 John Back
 * T109: Add EvtEtaLLPiPi model for eta' -> l l pi pi decays (l = e or mu),
   - courtesy of Aleksei Luchinsky (LHCb).
 
 29th Jun 2020 Michal Kreps
 * T103: Add missing include in EvtGenBase/EvtMatrix.hh.
 
 15th May 2020 Michal Kreps
 * D28: Add EvtThreeBodyPhsp (rectangular Dalitz plot selection) and
     EvtPhspDecaytimeCut (minimum decay time requirement) models.
 * D27: Fix initialisation of constants in EvtBTo3hCP model.
 
 ===
 ## R02-00-00
 
 24th Apr 2020 Michal Kreps
 * Update particle properties file evt.pdl to 2019 version of RPP by PDG.
 
 23rd Apr 2020 Tom Latham
 * Apply copyright and licence notices to all relevant files.
 
 17th Apr 2020 Tom Latham
 * Add text of GNU GPLv3 in COPYING, add (preliminary) list of authors in
   AUTHORS, and clean-up all source files, prior to applying copyright and
   licence notices.
 
 9th Apr 2020 Tom Latham
 * Improve, and document use of, standalone installation script.
 * Apply clang-format formatting to all C++ source files.
 
 8th Apr 2020 Tom Latham
 * One small modernisation change to EvtPhotosEngine to match that already
   applied in EvtTauolaEngine.
 
 8th Apr 2020 John Back
 * Fixed NonReson amplitude and the 4-momentum boosts used for angles in
   EvtLambdacPHH,
   - courtesy of Elisabeth Niel (LHCb).
 
 7th Apr 2020 Gerhard Raven, Tom Latham, Michal Kreps and John Back
 * Incorporate C++ modernization changes from Gerhard Raven (LHCb).
   - Merged modernize branch into master.
 
 9th Mar 2020 John Back
 * Add EvtAbsExternalGen::getDecayProb() to allow external generators to
   return a probability that can be used in EvtDecayProb (default = 1).
 
 6th Mar 2020 Andrii Verbytskyi and Tom Latham
 * Implement HepMC3 support: EvtHepMCEvent, external engines & cmake files.
 
 15th Nov 2019 John Back
 * Added EvtPsi2JpsiPiPi model for psi2S -> J/psi pi+ pi- decays based on hep-ph/1507.07985,
   - courtesy of Aleksei Luchinsky (LHCb).
 
 21st August 2019 Michal Kreps
 * Added the EvtDToKpienu decay model for D -> K pi e nu decays from BESIII,
   - courtesy of Liaoyuan Dong.
 
 16th July 2019 John Back
 * Correct imaginary sign for EvtComplex /= (EvtComplex c) operator.
 
 3rd July 2019 John Back
 * Added the EvtLambdacPHH decay model for Lc -> p K pi decays with K*(890),
   Delta++(1232) and Lambda(1520) resonances, based on the Fermilab E791
   analysis hep-ex/9912003v1,
   - courtesy of Elisabeth Niel and Patrick Robbe (LHCb).
 * Modified EvtResonance2 to accept other barrier factor radii.
 
 3rd July 2019 Michal Kreps
 * Make sure minimum mass for resonances with non-zero widths is larger than
   1e-4 GeV in EvtRelBreitWignerBarrierFact.
 
 3rd May 2019 John Back
 * Corrected EvtSLDiBaryonAmp bugs/issues in the BToDiBaryonlnupQCD model:
   - parity, amplitude terms and B momentum reference frame variables.
 * Corrected treament of spinor indices in EvtRareLb2Lll,
   - courtesy of Tom Blake and Michal Kreps (LHCb).
 * Updated the EvtBcVHad model to also handle Bc -> psi Ks K decays,
   - courtesy of Aleksei Luchinsky (LHCb).
 * Add new decay model EvtBsMuMuKK (BS_MUMUKK) for Bs to J/psi (mu+mu-) K+K-,
   - courtesy of Veronika Chobanova, Jeremy Dalseno, Diego Martinez Santos and
     Marcos Romero Lamas (LHCb).
 * Fix infinite loop during initialisation of the EvtBTo3hCP model via
   EvtCBTo3piP00 and EvtCBTo3piMPP,
   - courtesy of Peter Richardson (Durham).
 
 15th March 2019 Tom Latham
 * Implement cmake build system, replacing the old config method.
 
 30th Jan 2019 John Back
 * Fix modernization compiler errors and warnings.
 
 29th Jan 2019 Michal Kreps
 * Allow reading decay files which are missing end-of-line before end-of-file.
 
 21st December 2018 John Back
 * Imported C++ modernization changes from Gerhard Raven (LHCb).
 
 7th December 2018 John Back
 * Added the EvtBLLNuL (BLLNUL) model that generates rare B -> ell ell nu ell
   decays, where ell = e or mu,
   - courtesy of Anna Danilina and Nikolai Nikitin (LHCb).
 * Removed the EvtB2MuMuMuNu (BUTOMMMN) model, since its now replaced
   by the more general BLLNuL one.
 
 5th November 2018 John Back
 * Added the BToDiBaryonlnupQCD model for generating B to p N* l nu decays,
   where N can be any (exited) charged baryon (spin 1/2 or 3/2),
   - courtesy of Mark Smith and Ryan Newcombe (LHCb), with added code optimisations.
 
 17th October 2018 John Back
 * Added various decay models from LHCb EvtGenExtras package:
   - EvtBcVHad ("BC_VHAD"),
   - Evtbs2llGammaMNT ("BSTOGLLMNT"),
   - Evtbs2llGammaISRFSR ("BSTOGLLISRFSR"),
   - EvtbTosllMS ("BTOSLLMS"),
   - EvtbTosllMSExt ("BTOSLLMSEXT"),
   - EvtLb2Baryonlnu ("Lb2Baryonlnu"),
   - EvtLb2plnuLCSR ("Lb2plnuLCSR"),
   - EvtLb2plnuLQCD ("Lb2plnuLQCD"),
   - EvtFlatSqDalitz ("FLATSQDALITZ"),
   - EvtPhspFlatLifetime ("PHSPFLATLIFETIME").
 
 5th October 2018 John Back
 * Updated setupEvtGen.sh to work with the new HepForge Phabricator site.
 
 13th March 2018 John Back
 * Updated EvtPythiaEngine to correctly handle updates of various particle
   properties so that Pythia uses the same information as EvtGen (evt.pdl)
   for the generic and alias PYTHIA decay model.
 
 12th March 2018 John Back
 * Updated EvtBcXMuNu models (X = Scalar, Vector, Tensor) to generate
   Bc to D0(star) mu nu decays, with associated form factors in EvtBCXFF,
   - courtesy of Aleksei Luchinsky (LHCb).
 * Also generalised the calculation
   of their maximum probabilities by reusing the CalcMaxProb method in
   EvtSemiLeptonicAmp, which now allows for different Q^2 binning
   (default remains at 25 bins).
 
 ===
 ## R01-07-00
 
 13th December 2017 John Back
 * New tag incorporating all changes below.
 * Recommended external packages are
   (as used in the setupEvtGen.sh script):
   - HepMC 2.06.09,
   - pythia 8.230,
   - Photos++ 3.61
   - Tauola++ 1.1.6c.
 
 12th December 2017 John Back
 * Changed Pythia validation example decay files to use Pythia8 codes.
 
 6th December 2017 John Back
 * Modified the examples to use DECAY.DEC (see 25th April 2016) instead of
   DECAY_2010.DEC. Changed EvtExternalGenList to assume Pythia8 codes are
   used in decay files by default, which is the case for DECAY.DEC. Also
   updated the setupEvtGen.sh script to work with Pythia 8.2x versions.
 
 29th November 2017 John Back
 * Modified EvtSVP, EvtVVP and EvtTVP models to handle both radiative and
   two-lepton decays,
   - courtesy of Aleksei Luchinsky (LHCb).
 
 14th July 2017 John Back
 * Only create external generator objects if they don't already exist in
   EvtExternalGenFactory.
 * Modified configure script to work with Pythia 8.2x
 
 5th July 2017 Michal Kreps
 * Register the VTOSLL model.
 
 14th June 2017 John Back
 * Add isNeutralKaon() boolean function and corrected comments in EvtDDalitz.
 
 8th May 2017 Michal Kreps
 * Fix bug in EvtbTosllVectorAmp to recognise Bs --> K*bar mu mu decay as
   b --> d ll transition.
 
 8th May 2017 Michal Kreps
 * Significantly simplify way how we decide on decay mode and daughters
   ordering in DDalitz model.
   - With new code by definition all orderings of
     daughters in the decay file will yield same output.
 
 4th May 2017 John Back
 * Further fixes to DDalitz particle ordering (including charge-conjugates):
   - Mode 5:  D0 -> K- K0bar K+ and K+ K- K0bar
   - Mode 12: D0 -> pi0 pi- pi+ and pi+ pi0 pi-
   - Removed unneeded index ordering checks for mode 10 (D+ -> pi- pi+ pi+)
     and mode 11 (Ds+ -> pi- pi+ pi+)
 
 27th April 2017 John Back
 * Fixed DDalitz particle ordering for mode 7: D+ -> pi+ K- K+ and K+ pi+ K-
   and their charge-conjugates
 
 7th April 2017 John Back
 * Modified EvtGenExternal/EvtPythiaEngine to ensure that the EvtGen-based
   instances of Pythia8 (for generic and alias decays) use the same
   particle properties as defined by EvtGen,
   - courtesy Patrick Robbe (LHCb).
 
 5th April 2017 Michal Kreps
 * Fixed indexing in copy constructor of Evt3Rank3C, which would otherwise
   produce an infinite loop;
   - bug report from David Grellscheid.
 
 3rd November 2016 John Back
 * Modified EvtFlatQ2 model to work for all B -> X lepton lepton modes, as
   well as adding an extra phase space factor to correct for the dip at low
   q^2, which is enabled by using "FLATQ2 1" instead of just "FLATQ2" in the decay file,
   - courtesy of Marcin Chrzaszcz & Thomas Blake (LHCb).
 
 
 13th October 2016 John Back
 * Added the TauolaCurrentOption decfile keyword to select the hadronic
   current in Tauola; default is the BaBar-tuned current option (int = 1).
 * EvtParticles can store double attributes using the functions
   setAttributeDouble(name, double) and getAttributeDouble(name), which can
   be useful for storing and retrieving amplitude weights, for example.
   - The analogous EvtParticle integer attribute interface remains unchanged:
     setAttribute(name, int) and getAttribute(name).
 
 13th September 2016 John Back
 * Modified EvtTauolaEngine to use internal Tauola spin matrices for
   tau pair events by temporarily setting the PDG id of the mother as a
   boson, keeping the same 4-momentum.
 * The BaBar hadronic currents are now used by default.
 * Also added the ability to change some Tauola parameters
   using the "Define" keyword in decay files.
 * Added an example decay file
   illustrating the new features: validation/TauolaFiles/Btautau.dec
 
 9th September 2016 Michal Kreps
 * Reimplement code in EvtBTo3pi.F, EvtBTo3piMPP.F, EvtBTo3piP00.F and
   EvtBToKpipi.F in C++ in order to remove dependence on Fortran compiler.
   - With this, there is no internal Fortran code in EvtGen.
 
 ===
 ## R01-06-00
 
 1st June 2016 John Back
 * New tag incorporating all changes below.
 * Recommended external packages are
   - HepMC 2.06.09
   - pythia 8.186
   - Photos++ 3.61
   - Tauola++ 1.1.5
 
 28th April 2016 Michal Kreps
 * For Ds+ --> 2pi+ pi- there was double counting of branching fraction
   resulting in total branching fraction being 1.5 times larger than measured
   one.
   - Fix by revisiting submodes, which now fill total Ds --> 3pi.
 
 25th April 2016 Michal Kreps
 * Added DECAY.DEC/XML, which contain updated semileptonic charm and beauty
   branching fractions using the 2014 PDG, tuned to minimize disagreements
   between measurements and EvtGen for both inclusive and exclusive decays.
 * Updated the evt.pdl particle properties file to the PDG 2014 edition.
 * Implemented new LQCD form factors for Lb --> L mu mu from arXiv paper
   1602.01399 (EvtRareLbToLllFFlQCD); old LQCD form factors are removed.
 
 18th March 2016 John Back
 * Fixed incorrect spinor algebra used in S1 -> 1/2 S2, 1/2 -> S3 S4 decays
   in EvtDiracParticle::rotateToHelicityBasis() functions,
   - courtesy of Luis Miguel Garcia Martin and the IFIC Valencia LHCb group.
 
 19th Feburary 2016 John Back
 * Fixed bug in the definition of the initial spinor term Sinit in
   EvtRareLbToLll::HadronicAmpRS(),
   - from Tom Blake (LHCb).
 
 12th February 2016 John Back
 * From LHCb, added extensions to the EvtHQET2(FF) model for semitauonic
   decays from Brian Hamilton, which needs a patch to EvtSemiLeptonicAmp
   from Jack Wimberley to ensure that the q^2 range is physical when
   finding the maximum amplitude probability.
 
 2nd December 2015 John Back
 * From LHCb, added EvtKStopizmumu model for KS -> pi0 mu mu decays based on
   JHEP08(1998)004,
   - courtesy of Veronika Chobanova, Diego Martinez Santos and Jeremy Dalseno.
 * Added EvtConst::Fermi for Fermi coupling constant.
 
 ===
 ## R01-05-00
 
 21st October 2015 John Back
 * New tag incorporating all changes below.
 * Recommended external packages are
   - HepMC 2.06.09
   - pythia 8.186
   - Photos++ 3.61
   - Tauola++ 1.1.5
 
 * Added the EvtB2MuMuMuNu model for simulating the very rare four-leptonic
   decays B- -> mu+ mu- anti-nu_mu mu-,
   - courtesy Nikolai Nikitin.
 
 16th October 2015 John Back
 * Updated the configure script to automatically select the library names
   for PHOTOS++; version 3.56 and below uses Fortran, version 3.61 and above
   uses C++ only (default). Avoid using v3.60, since it does not work.
   This needs the PHOTOS libraries built before EvtGen is configured.
   Modified setupEvtGen.sh to use Photos++ v3.61.
 
 7th October 2015 John Back
 * Updated EvtGenExternal/EvtPhotosEngine to check that additional particles
   from the outgoing vertex are indeed (FSR) photons, since later versions of
   PHOTOS introduce pair emission, where particles may not always be photons.
 * Added the genRootDecayChain.cc validation program to create ROOT files
   containing information about the complete decay tree. Two example test
   decay files BKstarGamma.dec and BuDst0rhop.dec can be used with this; the
   first tests PHOTOS, the second looks at sequential decay chain storage.
   The plotBKstarGamma.C ROOT macro can be used for B -> K* gamma plots.
 
 2nd October 2015 John Back
 * Modified EvtSVPHelAmp and added a new EvtSVPHelCPMix model, implementing
   the complete mixing phenomenology of Bs to vector gamma decays,
   - courtesy of Clara Remon (LHCb).
 * EvtD0mixDalitz code: cleanup, inverted q/p for decays of D0bar (simplifies
   user decay files) and fixed y parameter bug,
   - courtesy of Jordi Tico (LHCb).
 * Changed the initialisation order of the infrared cut-off in EvtPhotosEngine.
   This actually has no effect, since the exponentiation function sets it to the
   same 1e-7 value, but it is now in the correct order if we need to update it.
 * Removed all remaining obsolete pragma (Win32) warnings from some classes.
 
 23rd September 2015 Michal Kreps
 * Reimplement the real Spence function in C++ and removed its fortran
   implementation.
 
 15th September 2015 Michal Kreps
 * Fixed accessed uninitialised memory in EvtPDL.cpp, line 213.
 * Modified the configure and setupEvtGen.sh scripts to work on Mac; needed
   Mac compilation patch files added to the new "platform" subdirectory.
 
 10th September 2015 John Back
 * Updated setupEvtGen.sh to use the recommended external packages:
   - HepMC 2.06.09, pythia 8.186, Photos++ 3.56 and Tauola++ 1.1.5.
 * Fixed form-factor calculations for the BTOSLLBALL model 6 used to generate
   b -> sll decays,
   - courtesy of Christoph Langenbruch and David Loh (LHCb).
   - Affects B->K*ll, B->rholl and B->omegall, particularly the electron modes.
 * In the validation directory, added runPhotosTest.sh for testing FSR in
   Upsilon(4S) -> e+ e- decays, and changed the plot comparison scripts to
   use the 2nd directory "oldRootFiles" (which could be a soft-link) for
   including ROOT histograms made from a previous version of EvtGen.
 
 27th August 2015 John Back
 * Added Mersenne-Twister random number generator (RNG) EvtMTRandomEngine.
   - It requires c++11 compiler features (>= gcc 4.7), which should
     automatically be enabled by the configure script.
   - Introduced the preprocessor environment variable EVTGEN_CPP11 for c++11
     features.
   - EvtMTRandomEngine is the default RNG for the validation and test examples
     if c++11 features are enabled.
 * Added a phase-space test validation/genPHSP.sh and PhaseSpacePlots.C to
   visually check the flatness of Dalitz plots in order to ensure that the
   RNG is not producing biased results that depend on particle ordering.
 * Added the models EvtbsToLLLLAmp and EvtbsToLLLLHyperCP for
   B0_q -> l+ l- l+ l- decays (SM and one supersymmetric scenario),
   - courtesy of Nikolai Nikitin and Konstantin Toms.
   - Documentation provided in doc/evt_BQTOLLLL_model.pdf and
     doc/evt_BQTOLLLLHYPERCP_model.pdf.
 
 * Changed the installation and set-up script name to be just setupEvtGen.sh;
   it uses the VERSION variable to specify the required tag. List of tags
   are available using either "svn ls -v http://svn.cern.ch/guest/evtgen/tags"
   or by going to http://svn.cern.ch/guest/evtgen/tags in a web browser.
 
 12th June 2015 John Back
 * Changed the width of chi_b1 in evt.pdl from 9.8928 GeV (!) to zero.
 
 1st May 2015 John Back
 * Added Bc -> scalar ell nu (EvtBcSMuNu) and Bc -> tensor ell nu
   (EvtBcTMuNu) decays,
   - courtesy of Jack Wimberley, LHCb.
   - Also included the chi_c1 mode for EvtBcVMuNu.
 
 ===
 ## R01-04-00
 
 2nd April 2015 John Back
 * Removed the EvtStdlibRandomEngine class since this can produce biases
   to kinematic distributions when one or more of the daughters is a
   resonance, such as B0 -> K pi psi
   - (thanks to Antonio Augusto Alves Jr who discovered this issue).
   - EvtSimpleRandomEngine is now the default random number generator in the
     validation and test examples.
 * Incorporated several additions and modifications from LHCb:
   a) From Michal Kreps, Tom Blake & Christoph Langenbruch,
      added rare Lb --> Lambda^(*) ell ell models described in
      arXiv:1108.6129, with various form factors from Gutsche et al.
      (arXiv:1301.3737) and lattice QCD  (arXiv:1212.4827)
   b) From Andrew Crocombe, implemented Bs --> K* form factors
      from Ball-Zwicky and z-parametrization form factors from
      arXiv:1006.4945 for EvtbTosllBallFF
   c) Christoph Langenbruch fixed the Bs -> phi ll form factors in
      EvtbTosllBallFF; T3 showed a non-physical pole at very low
      q2 which significantly affected the electron mode
   d) From Michal Kreps, removed semicolons from wrong places to
      clear warnings when compiled with the "-pedantic" option.
 
 9th October 2014 John Back
 * Change svnweb.cern.ch to svn.cern.ch in the setup script.
 
 1st April 2014 John Back
 * In EvtReport, modified the logging output severity status flags
   to have the "EVTGEN_" prefix, e.g. INFO becomes EVTGEN_INFO.
 * The global report() function has been renamed to EvtGenReport().
 
 31st March 2014 John Back
 * Added the ability to store named attributes for EvtParticles in the
   form of a map<string, int>. The setAttribute(name, value) stores the
   required value, while getAttribute(name) retrieves the integer value.
   This is used in EvtPhotosEngine to specify the final-state radiation
   "FSR" attribute to 1 for any additional photons (EvtPhotonParticles)
   created by Photos++. It also stores the "ISR" attribute, but this
   is always set to zero, since only FSR photons are created.
   If the named attribute does not exist, then getAttribute() returns zero.
 
 29th January 2014 Daniel Craik
 * Removed mass assertion on GS shape in EvtDalitzReso to allow it to also
   be used for charged rho resonances.
 
 27th January 2014 John Back
 * Minor corrections to Vub models to remove further gcc 4.8 warnings.
 * Updated configure script to work for MacOS clang (from Genser team).
 
 ===
 ## R01-03-00
 
 9th January 2014 John Back
 * New tag version "1.3.0", incorporating all changes below.
 * Replaced auto-install script to work with this version as well as
   the latest versions of all external generator packages.
 * Updated README to mention the new CERN-based web pages for Photos++
   and Tauola++.
 
 8th January 2014 John Back
 * Fix gcc 4.6 and 4.8 compilation warnings,
   - courtesy of Patrick Robbe (LHCb);
   - main changes are removal of unused variables.
 * Changed the EvtPythiaEngine class and configure script to use new
   Pythia 8 header locations; Pythia 8.180 or above is now required.
 
 7th January 2014 John Back
 * Modified EvtBCVFF to correct the Kiselev form factors
   - from Jack Wimberley (LHCb).
 
 9th October 2013 Daniel Craik
 * Added Gounaris-Sakurai and Gaussian shapes to EvtGenericDalitz
   and set sensible defaults for LASS parameters.
 
 19th September 2013 John Back
 * Modified EvtGenExternal/EvtPythiaEngine to keep track of any new
   particles that are added to the default Pythia database to avoid
   duplicating particle/anti-particle entries that could override
   previously defined Pythia decay chains.
 
 18th September 2013 John Back
 * Added Mac OS flags for the configure script and src/Makefile.
 
 15th July 2013 Daniel Craik
 * Added flag to turn on scaling of LASS amplitude by M/q in EvtDalitzReso
 
 15th July 2013 Daniel Craik
 * EvtParserXML now accepts file names containing environment variables,
   exponential non-resonant shape in EvtDalitzReso now defined as exp(-alpha*m^2),
   LASS shape in EvtDalitzReso now takes a cutoff parameter
 
 4th July 2013 Daniel Craik
 * Added LASS, exponential non-resonant and linear non-resonant shapes to EvtGenericDalitz.
 
 3rd July 2013 Daniel Craik
 * Fixed auto-install script for R01-02-00.
 
 1st July 2013 Daniel Craik
 * Added auto-install script for R01-02-00.
 
 ===
 ## R01-02-00
 
 15th May 2013 John Back
 * New tag, version "1.2.0", incorporating all changes below.
 
 14th May 2013 Michal Kreps
 * Added Blatt-Weisskopf barrier factors up to L=5 in
   EvtGenBase/EvtBlattWeisskopf::compute().
 
 14th May 2013 John Back
 * Added additional entries (appended at the end) to the evt.pdl particle
   data file
   - courtesy of Romulus Godang and Belle II colleagues.
 
 14th March 2013 John Back
 * Added the method EvtParticle::getPDGId() to get the PDG integer for a
   particle directly (which just calls EvtPDL::getStdHep()).
 * Added a check in EvtPhotosEngine::doDecay to skip Photos if a given
   particle has too many daughters (>= 10) to avoid a problem with a
   hard coded upper limit in the Photos PHOENE subroutine.
 
 2nd February 2013 Daniel Craik
 * Updated EvtDalitzTable to estimate probMax when it is missing from a
   Dalitz model.
 
 1st  February 2013 John Back
 * Added the ability to read in Pythia 6 commands in ascii decay files in
   EvtDecayTable::readDecayFile (this was already possible in xml files).
 * Modified the Photos++ engine default settings to be more suited to B
   decays (from LHCb defaults).
 
 31st January 2013 John Back
 * Added the ability to read in Pythia 8 commands in ascii decay files
   in EvtDecayTable::readDecayFile. They can be set using the syntax:
   "PythiaTypeParam module:variable=value", where Type = Generic, Alias or
   Both for specifying whether the parameter is for the generic or alias
   Pythia decay engines (or both). The 2nd argument must not contain any
   spaces. Fixed the list of commands strings being used in the
   EvtPythiaEngine class (i.e. Pythia parameters that can be set via decay
   files).
 
 31st January 2013 Daniel Craik
 * Added named parameters to various decay models.
 
 30th January 2013 John Back
 * Fixed some of the parameter arguments used in the EvtSVSCPiso model.
 
 24th January 2013 John Back
 * Set the Photos++ and Tauola++ models to use the EvtGen random number
   engine when useEvtGenRandom is set to true in the EvtExternalGenList
   constructor.
 
 23rd January 2013 John Back
 * Added EvtGenExternal/EvtPythiaRandom to allow the use of the EvtGen
   random number engine to also be used for the random engine for Pythia 8.
 * Added a boolean (useEvtGenRandom, default = true) within the
   EvtExternalGenList constructor to use this feature.
 
 18th December 2012 John Back
 * Corrected some wrong daughter ordering assignments for decay modes 5 and
   12 in EvtDDalitz. Updated validation/DalitzDecays.xml to also contain
   D decay mode 12, as well as various modes that may use K_S0 and K_L0.
 * Added validation/genDDalitzModes.sh and updated validation/compareDalitz.C to
   do a complete comparison of the D Dalitz modes with the xml versions.
 
 11th December 2012 Daniel Craik
 * Updated the Xml parser to support named model parameters.
 * Updated the generic Dalitz model to use named model parameters as an example.
 
 15th October 2012 John Back
 * Make EvtSimpleRandomEngine inherit from EvtRandomEngine to avoid
   crash in EvtGen.cpp when no random engine is defined
   - (from Bjoern Spruck).
 
 ===
 ## R01-01-00
 
 4th October 2012 John Back
 * New tag, version "1.1.0", incorporating all changes below.
 
 * Provide proper default constructors for EvtVector4R and
   EvtPhotonParticle. Modified the validation and test code to also
   compile/link in the case of no external generators being included.
 
 3rd October 2012 John Back
 * Corrected the t3 vector form factor values for the Ball-Zwicky 2005
   model (modelId = 6) in EvtbTosllBallFF::getVectorFF(), which
   were set to t3tilde instead.
 
 18th September 2012 John Back
 * Moved the external generator engines to a new sub-directory
   EvtGenExternal. Building the code now creates 2 libraries:
   libEvtGen.so (Base+Models) and libEvtGenExternal.so.
   - This allows anyone to ignore using the new external generators
     if required (by not creating/loading the 2nd library).
 * Added prefix option to the configure script/Makefile to allow the user
   to specify an installation directory for the include files, libraries,
   DECAY.DEC and evt.pdl files (for Genser).
 
 14th September 2012 Michal Kreps
 * Fixed the calculation of the angle between decay planes in the function
   EvtKine::EvtDecayAngleChi. Fixed typo in EvtLb2Lll decay model. Only
   some NP scenarious could be affected, SM one is definitely unaffected.
 
 13th September 2012 John Back
 * Added the use of the environment variables EVTGEN_PHOTOS, EVTGEN_PYTHIA
   and EVTGEN_TAUOLA to specify if the Photos, Pythia and/or Tauola engine
   classes are used or not. These variables are set by the configure script,
   depending if the library paths are specified for these generators.
 
 ===
 ## R01-00-01
 
 12th September 2012 John Back
 * New tag incorporating all changes below, since R01-00-00.
 
 11th September 2012 John Back
 * Modified the Photos and Tauola engine classes to use the new
   Photospp and Tauolapp namespaces that are present in the latest
   versions of Photos++(3.5) and Tauola++(1.0.7).
 * Updated the configure file to get the correct location of the Tauola++
   include files.
 * Added the D0->pi+pi-pi0 decay mode in EvtDDalitz
   - from Marco Gersabeck and Frederic Dreyer (LHCb).
 * Added new decay models/classes from Alexey Luchinsky (LHCb):
   EvtBcVMuNu, EvtTVP, EvtWnPi, EvtSVP, EvtXPsiGamma, EvtBcVNpi
 
 29th June 2012 John Back
 * Corrected mass(squared) variables filled in the Dalitz TTree in
   validation/genExampleRootFiles.
 
 15th May 2012 Daniel Craik
 * Updated EvtD0gammaDalitz to deal with D mesons from neutral B->DK
 * Added save function to validation/compareDalitz.C.
 
 11th May 2012 Daniel Craik
 * Replaced BaBar specific configuration for BlattWeisskopf birth factors.
 * Updated XML conversion script to handle new configuration.
 * Fixed some bugs in the XML conversion script related to particle
   modifications.
 
 9th May 2012 Daniel Craik
 * Added latex documentation for xml decay files.
 
 2nd May 2012 Daniel Craik
 * Added resDaughters attribute to the Dalitz resonance xml tag to
   simplify defining symmetric resonances. Updated validation xml files to
   use the new functionality.
 
 27th April 2012 Daniel Craik
 * Upgraded EvtGenericDalitz to use EvtDalitzReso for resonances.
 * Added validation to compare EvtGenericDalitz to all 11 EvtDDalitz modes.
 * Added a root macro to quickly compare two Dalitz decays for validation.
 
 24th April 2012 John Back
 * Solved two bugs in the EvtD0gammaDalitz model (from Jordi Tico, LHCb):
   configuration of the conjugated model, and using only the B charge
   to determine the model used, not the D flavour.
 
 17th April 2012 Daniel Craik
 * Updated the GenericDalitz validation code to use the same probMax
   values as DDalitz.
 * Added XML decay file parsing to EvtGen::readUDecay.
   - Dec files are still the default.
 
 30th March 2012 John Back
 * Update maximum probability values in EvtDDalitz::initProbMax()
   for all DDalitz modes.
 
 23rd March 2012 John Back
 * Added the EvtEta2MuMuGamma decay model from LHCb.
 
 21st March 2012 John Back
 * Added EvtD0gammaDalitz decay model from LHCb.
 
 20th March 2012 Daniel Craik
 * Added backwards compatibility for Pythia 6 commands in the XML configuration.
 * Updated decay file conversion tool to convert JetSetPar lines to pythia6Param
   tags.
 
 19th March 2012 Daniel Craik
 * Added infrastructure to pass commands to external generators.
 * XML config now takes Pythia8 configuration commands.
 
 16th March 2012 Daniel Craik
 * Added the ability to define particles from the PDL for Dalitz decay
   resonances instead of defining mass, width and spin seperately.
 * Renamed the lifetime attribute of Dalitz decay resonaces to width to avoid
   confusion.
 * Added further validation code for the generic Dalitz model.
 
 15th March 2012 Daniel Craik
 * Added validation code for xml decay files and the generic Dalitz model.
 
 ===
 ## R01-00-00
 
 6th March 2012 John Back
 * First official version for Genser (evtgen 1.0.0) that includes
  support for external generators: Pythia8, Photos++ and Tauola++.
 * This also includes a preliminary version of creating Dalitz plot
  decay models using EvtGenericDalitz.
diff --git a/src/EvtGenBase/EvtParticle.cpp b/src/EvtGenBase/EvtParticle.cpp
index 09cffee..7eed7dc 100644
--- a/src/EvtGenBase/EvtParticle.cpp
+++ b/src/EvtGenBase/EvtParticle.cpp
@@ -1,1317 +1,1317 @@
 
 /***********************************************************************
 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
 *                                                                      *
 * This file is part of EvtGen.                                         *
 *                                                                      *
 * EvtGen is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU General Public License as published by *
 * the Free Software Foundation, either version 3 of the License, or    *
 * (at your option) any later version.                                  *
 *                                                                      *
 * EvtGen is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
 * GNU General Public License for more details.                         *
 *                                                                      *
 * You should have received a copy of the GNU General Public License    *
 * along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
 ***********************************************************************/
 
 #include "EvtGenBase/EvtParticle.hh"
 
 #include "EvtGenBase/EvtCPUtil.hh"
 #include "EvtGenBase/EvtDecayTable.hh"
 #include "EvtGenBase/EvtDiracParticle.hh"
 #include "EvtGenBase/EvtGenKine.hh"
 #include "EvtGenBase/EvtId.hh"
 #include "EvtGenBase/EvtIdSet.hh"
 #include "EvtGenBase/EvtNeutrinoParticle.hh"
 #include "EvtGenBase/EvtPDL.hh"
 #include "EvtGenBase/EvtParticleFactory.hh"
 #include "EvtGenBase/EvtPatches.hh"
 #include "EvtGenBase/EvtPhotonParticle.hh"
 #include "EvtGenBase/EvtRadCorr.hh"
 #include "EvtGenBase/EvtRandom.hh"
 #include "EvtGenBase/EvtRaritaSchwingerParticle.hh"
 #include "EvtGenBase/EvtReport.hh"
 #include "EvtGenBase/EvtScalarParticle.hh"
 #include "EvtGenBase/EvtSecondary.hh"
 #include "EvtGenBase/EvtStatus.hh"
 #include "EvtGenBase/EvtStdHep.hh"
 #include "EvtGenBase/EvtStringParticle.hh"
 #include "EvtGenBase/EvtTensorParticle.hh"
 #include "EvtGenBase/EvtVectorParticle.hh"
 
 #include <sys/stat.h>
 
 #include <iostream>
 #include <math.h>
 #include <stdio.h>
 #include <stdlib.h>
 
 using std::endl;
 
 EvtParticle::~EvtParticle()
 {
     delete _decayProb;
 }
 
 EvtParticle::EvtParticle()
 {
     _ndaug = 0;
     _parent = nullptr;
     _channel = -10;
     _t = 0.0;
     _genlifetime = 1;
     _first = 1;
     _isInit = false;
     _validP4 = false;
     _isDecayed = false;
     _decayProb = nullptr;
     _intAttributes.clear();
     _dblAttributes.clear();
     //   _mix=false;
 }
 
 void EvtParticle::setFirstOrNot()
 {
     _first = 0;
 }
 void EvtParticle::resetFirstOrNot()
 {
     _first = 1;
 }
 
 void EvtParticle::setChannel( int i )
 {
     _channel = i;
 }
 
 EvtParticle* EvtParticle::getParent() const
 {
     return _parent;
 }
 
 void EvtParticle::setLifetime( double tau )
 {
     _t = tau;
 }
 
 void EvtParticle::setLifetime()
 {
     if ( _genlifetime ) {
         _t = -log( EvtRandom::Flat() ) * EvtPDL::getctau( getId() );
     }
 }
 
 double EvtParticle::getLifetime() const
 {
     return _t;
 }
 
 void EvtParticle::addDaug( EvtParticle* node )
 {
     node->_daug[node->_ndaug++] = this;
     _ndaug = 0;
     _parent = node;
 }
 
 int EvtParticle::firstornot() const
 {
     return _first;
 }
 
 EvtId EvtParticle::getId() const
 {
     return _id;
 }
 
 int EvtParticle::getPDGId() const
 {
     return EvtPDL::getStdHep( _id );
 }
 
 EvtSpinType::spintype EvtParticle::getSpinType() const
 {
     return EvtPDL::getSpinType( _id );
 }
 
 int EvtParticle::getSpinStates() const
 {
     return EvtSpinType::getSpinStates( EvtPDL::getSpinType( _id ) );
 }
 
 const EvtVector4R& EvtParticle::getP4() const
 {
     return _p;
 }
 
 int EvtParticle::getChannel() const
 {
     return _channel;
 }
 
 size_t EvtParticle::getNDaug() const
 {
     return _ndaug;
 }
 
 double EvtParticle::mass() const
 {
     return _p.mass();
 }
 
 void EvtParticle::setDiagonalSpinDensity()
 {
     _rhoForward.setDiag( getSpinStates() );
 }
 
 void EvtParticle::setVectorSpinDensity()
 {
     if ( getSpinStates() != 3 ) {
         EvtGenReport( EVTGEN_ERROR, "EvtGen" )
             << "Error in EvtParticle::setVectorSpinDensity" << endl;
         EvtGenReport( EVTGEN_ERROR, "EvtGen" )
             << "spin_states:" << getSpinStates() << endl;
         EvtGenReport( EVTGEN_ERROR, "EvtGen" )
             << "particle:" << EvtPDL::name( _id ).c_str() << endl;
         ::abort();
     }
 
     EvtSpinDensity rho;
 
     //Set helicity +1 and -1 to 1.
     rho.setDiag( getSpinStates() );
     rho.set( 1, 1, EvtComplex( 0.0, 0.0 ) );
 
     setSpinDensityForwardHelicityBasis( rho );
 }
 
 void EvtParticle::setSpinDensityForwardHelicityBasis( const EvtSpinDensity& rho )
 {
     EvtSpinDensity R = rotateToHelicityBasis();
 
     assert( R.getDim() == rho.getDim() );
 
     int n = rho.getDim();
 
     _rhoForward.setDim( n );
 
     int i, j, k, l;
 
     for ( i = 0; i < n; i++ ) {
         for ( j = 0; j < n; j++ ) {
             EvtComplex tmp = 0.0;
             for ( k = 0; k < n; k++ ) {
                 for ( l = 0; l < n; l++ ) {
                     tmp += R.get( l, i ) * rho.get( l, k ) *
                            conj( R.get( k, j ) );
                 }
             }
             _rhoForward.set( i, j, tmp );
         }
     }
 }
 
 void EvtParticle::setSpinDensityForwardHelicityBasis( const EvtSpinDensity& rho,
                                                       double alpha, double beta,
                                                       double gamma )
 {
     EvtSpinDensity R = rotateToHelicityBasis( alpha, beta, gamma );
 
     assert( R.getDim() == rho.getDim() );
 
     int n = rho.getDim();
 
     _rhoForward.setDim( n );
 
     int i, j, k, l;
 
     for ( i = 0; i < n; i++ ) {
         for ( j = 0; j < n; j++ ) {
             EvtComplex tmp = 0.0;
             for ( k = 0; k < n; k++ ) {
                 for ( l = 0; l < n; l++ ) {
                     tmp += R.get( l, i ) * rho.get( l, k ) *
                            conj( R.get( k, j ) );
                 }
             }
             _rhoForward.set( i, j, tmp );
         }
     }
 }
 
 void EvtParticle::initDecay( bool useMinMass )
 {
     EvtParticle* p = this;
     // carefull - the parent mass might be fixed in stone..
     EvtParticle* par = p->getParent();
     double parMass = -1.;
     if ( par ) {
         if ( par->hasValidP4() )
             parMass = par->mass();
         for ( size_t i = 0; i < par->getNDaug(); i++ ) {
             EvtParticle* tDaug = par->getDaug( i );
             if ( p != tDaug )
                 parMass -= EvtPDL::getMinMass( tDaug->getId() );
         }
     }
 
     if ( _isInit ) {
         //we have already been here - just reroll the masses!
         if ( _ndaug > 0 ) {
             for ( size_t ii = 0; ii < _ndaug; ii++ ) {
                 if ( _ndaug == 1 ||
                      EvtPDL::getWidth( p->getDaug( ii )->getId() ) > 0.0000001 )
                     p->getDaug( ii )->initDecay( useMinMass );
                 else
                     p->getDaug( ii )->setMass(
                         EvtPDL::getMeanMass( p->getDaug( ii )->getId() ) );
             }
         }
 
         EvtId* dauId = nullptr;
         double* dauMasses = nullptr;
         if ( _ndaug > 0 ) {
             dauId = new EvtId[_ndaug];
             dauMasses = new double[_ndaug];
             for ( size_t j = 0; j < _ndaug; j++ ) {
                 dauId[j] = p->getDaug( j )->getId();
                 dauMasses[j] = p->getDaug( j )->mass();
             }
         }
         EvtId* parId = nullptr;
         EvtId* othDauId = nullptr;
         EvtParticle* tempPar = p->getParent();
         if ( tempPar ) {
             parId = new EvtId( tempPar->getId() );
             if ( tempPar->getNDaug() == 2 ) {
                 if ( tempPar->getDaug( 0 ) == this )
                     othDauId = new EvtId( tempPar->getDaug( 1 )->getId() );
                 else
                     othDauId = new EvtId( tempPar->getDaug( 0 )->getId() );
             }
         }
         if ( p->getParent() && _validP4 == false ) {
             if ( !useMinMass ) {
                 p->setMass( EvtPDL::getRandMass( p->getId(), parId, _ndaug, dauId,
                                                  othDauId, parMass, dauMasses ) );
             } else
                 p->setMass( EvtPDL::getMinMass( p->getId() ) );
         }
         if ( parId )
             delete parId;
         if ( othDauId )
             delete othDauId;
         if ( dauId )
             delete[] dauId;
         if ( dauMasses )
             delete[] dauMasses;
         return;
     }
 
     //Will include effects of mixing here
     //added by Lange Jan4,2000
     static EvtId BS0 = EvtPDL::getId( "B_s0" );
     static EvtId BSB = EvtPDL::getId( "anti-B_s0" );
     static EvtId BD0 = EvtPDL::getId( "B0" );
     static EvtId BDB = EvtPDL::getId( "anti-B0" );
     static EvtId D0 = EvtPDL::getId( "D0" );
     static EvtId D0B = EvtPDL::getId( "anti-D0" );
     static EvtId U4S = EvtPDL::getId( "Upsilon(4S)" );
     static EvtIdSet borUps{ BS0, BSB, BD0, BDB, U4S };
 
     //only makes sense if there is no parent particle which is a B or an Upsilon
     bool hasBorUps = false;
     if ( getParent() && borUps.contains( getParent()->getId() ) )
         hasBorUps = true;
     //    if ( (getNDaug()==0)&&(getParent()==0) && (getId()==BS0||getId()==BSB||getId()==BD0||getId()==BDB)){
     EvtId thisId = getId();
     // remove D0 mixing for now.
     //  if ( (getNDaug()==0 && !hasBorUps) && (thisId==BS0||thisId==BSB||thisId==BD0||thisId==BDB||thisId==D0||thisId==D0B)){
     if ( ( getNDaug() == 0 && !hasBorUps ) &&
          ( thisId == BS0 || thisId == BSB || thisId == BD0 || thisId == BDB ) ) {
         double t;
         int mix;
         EvtCPUtil::getInstance()->incoherentMix( getId(), t, mix );
         setLifetime( t );
 
         if ( mix ) {
             EvtScalarParticle* scalar_part;
 
             scalar_part = new EvtScalarParticle;
             if ( getId() == BS0 ) {
                 EvtVector4R p_init( EvtPDL::getMass( BSB ), 0.0, 0.0, 0.0 );
                 scalar_part->init( EvtPDL::chargeConj( getId() ), p_init );
             } else if ( getId() == BSB ) {
                 EvtVector4R p_init( EvtPDL::getMass( BS0 ), 0.0, 0.0, 0.0 );
                 scalar_part->init( EvtPDL::chargeConj( getId() ), p_init );
             } else if ( getId() == BD0 ) {
                 EvtVector4R p_init( EvtPDL::getMass( BDB ), 0.0, 0.0, 0.0 );
                 scalar_part->init( EvtPDL::chargeConj( getId() ), p_init );
             } else if ( getId() == BDB ) {
                 EvtVector4R p_init( EvtPDL::getMass( BD0 ), 0.0, 0.0, 0.0 );
                 scalar_part->init( EvtPDL::chargeConj( getId() ), p_init );
             } else if ( getId() == D0 ) {
                 EvtVector4R p_init( EvtPDL::getMass( D0B ), 0.0, 0.0, 0.0 );
                 scalar_part->init( EvtPDL::chargeConj( getId() ), p_init );
             } else if ( getId() == D0B ) {
                 EvtVector4R p_init( EvtPDL::getMass( D0 ), 0.0, 0.0, 0.0 );
                 scalar_part->init( EvtPDL::chargeConj( getId() ), p_init );
             }
 
             scalar_part->setLifetime( 0 );
             scalar_part->setDiagonalSpinDensity();
 
             insertDaugPtr( 0, scalar_part );
 
             _ndaug = 1;
             _isInit = true;
             p = scalar_part;
             p->initDecay( useMinMass );
             return;
         }
     }
 
     EvtDecayBase* decayer;
     decayer = EvtDecayTable::getInstance()->getDecayFunc( p );
 
     if ( decayer ) {
         p->makeDaughters( decayer->nRealDaughters(), decayer->getDaugs() );
         //then loop over the daughters and init their decay
         for ( size_t i = 0; i < p->getNDaug(); i++ ) {
             //      std::cout << EvtPDL::name(p->getDaug(i)->getId()) << " " << i << " " << p->getDaug(i)->getSpinType() << " " << EvtPDL::name(p->getId()) << std::endl;
             if ( EvtPDL::getWidth( p->getDaug( i )->getId() ) > 0.0000001 )
                 p->getDaug( i )->initDecay( useMinMass );
             else
                 p->getDaug( i )->setMass(
                     EvtPDL::getMeanMass( p->getDaug( i )->getId() ) );
         }
     }
 
     int j;
     EvtId* dauId = nullptr;
     double* dauMasses = nullptr;
     int nDaugT = p->getNDaug();
     if ( nDaugT > 0 ) {
         dauId = new EvtId[nDaugT];
         dauMasses = new double[nDaugT];
         for ( j = 0; j < nDaugT; j++ ) {
             dauId[j] = p->getDaug( j )->getId();
             dauMasses[j] = p->getDaug( j )->mass();
         }
     }
 
     EvtId* parId = nullptr;
     EvtId* othDauId = nullptr;
     EvtParticle* tempPar = p->getParent();
     if ( tempPar ) {
         parId = new EvtId( tempPar->getId() );
         if ( tempPar->getNDaug() == 2 ) {
             if ( tempPar->getDaug( 0 ) == this )
                 othDauId = new EvtId( tempPar->getDaug( 1 )->getId() );
             else
                 othDauId = new EvtId( tempPar->getDaug( 0 )->getId() );
         }
     }
     if ( p->getParent() && p->hasValidP4() == false ) {
         if ( !useMinMass ) {
             p->setMass( EvtPDL::getRandMass( p->getId(), parId, p->getNDaug(),
                                              dauId, othDauId, parMass,
                                              dauMasses ) );
         } else {
             p->setMass( EvtPDL::getMinMass( p->getId() ) );
         }
     }
     if ( parId )
         delete parId;
     if ( othDauId )
         delete othDauId;
     if ( dauId )
         delete[] dauId;
     if ( dauMasses )
         delete[] dauMasses;
     _isInit = true;
 }
 
 void EvtParticle::decay()
 {
     //P is particle to decay, typically 'this' but sometime
     //modified by mixing
     EvtParticle* p = this;
     //Did it mix?
     //if ( p->getMixed() ) {
     //should take C(p) - this should only
     //happen the first time we call decay for this
     //particle
     //p->takeCConj();
     // p->setUnMixed();
     //}
 
     EvtDecayBase* decayer;
     decayer = EvtDecayTable::getInstance()->getDecayFunc( p );
     //  if ( decayer ) {
     //    EvtGenReport(EVTGEN_INFO,"EvtGen") << "calling decay for " << EvtPDL::name(p->getId()) << " " << p->mass() << " " << p->getP4() << " " << p->getNDaug() << " " << p << endl;
     //    EvtGenReport(EVTGEN_INFO,"EvtGen") << "NDaug= " << decayer->getNDaug() << endl;
     //    int ti;
     //    for ( ti=0; ti<decayer->getNDaug(); ti++)
     //      EvtGenReport(EVTGEN_INFO,"EvtGen") << "Daug " << ti << " " << EvtPDL::name(decayer->getDaug(ti)) << endl;
     //  }
     //if (p->_ndaug>0) {
     //      EvtGenReport(EVTGEN_INFO,"EvtGen") <<"Is decaying particle with daughters!!!!!"<<endl;
     //     ::abort();
     //return;
     //call initdecay first - April 29,2002 - Lange
     //}
 
     //if there are already daughters, then this step is already done!
     // figure out the masses
     bool massTreeOK( true );
     if ( _ndaug == 0 ) {
         massTreeOK = generateMassTree();
     }
 
     if ( massTreeOK == false ) {
         EvtGenReport( EVTGEN_INFO, "EvtGen" )
             << "Could not decay " << EvtPDL::name( p->getId() ) << " with mass "
             << p->mass() << " to decay channel number " << _channel << endl;
         _isDecayed = false;
         return;
     }
 
     static EvtId BS0 = EvtPDL::getId( "B_s0" );
     static EvtId BSB = EvtPDL::getId( "anti-B_s0" );
     static EvtId BD0 = EvtPDL::getId( "B0" );
     static EvtId BDB = EvtPDL::getId( "anti-B0" );
     // static EvtId D0=EvtPDL::getId("D0");
     // static EvtId D0B=EvtPDL::getId("anti-D0");
 
     EvtId thisId = getId();
     // remove D0 mixing for now..
     //  if ( _ndaug==1 &&  (thisId==BS0||thisId==BSB||thisId==BD0||thisId==BDB||thisId==D0||thisId==D0B) ) {
     if ( _ndaug == 1 &&
          ( thisId == BS0 || thisId == BSB || thisId == BD0 || thisId == BDB ) ) {
         p = p->getDaug( 0 );
         decayer = EvtDecayTable::getInstance()->getDecayFunc( p );
     }
     //now we have accepted a set of masses - time
     if ( decayer != nullptr ) {
         decayer->makeDecay( p );
     } else {
         p->_rhoBackward.setDiag( p->getSpinStates() );
     }
 
     _isDecayed = true;
     return;
 }
 
 bool EvtParticle::generateMassTree()
 {
     bool isOK( true );
 
     double massProb = 1.;
     double ranNum = 2.;
     int counter = 0;
     EvtParticle* p = this;
     while ( massProb < ranNum ) {
         //check it out the first time.
         p->initDecay();
         massProb = p->compMassProb();
         ranNum = EvtRandom::Flat();
         counter++;
 
         if ( counter > 10000 ) {
             if ( counter == 10001 ) {
                 EvtGenReport( EVTGEN_INFO, "EvtGen" )
                     << "Too many iterations to determine the mass tree. Parent mass= "
                     << p->mass() << " " << massProb << endl;
                 p->printTree();
                 EvtGenReport( EVTGEN_INFO, "EvtGen" )
                     << "will take next combo with non-zero likelihood\n";
             }
             if ( massProb > 0. )
                 massProb = 2.0;
             if ( counter > 20000 ) {
                 // one last try - take the minimum masses
                 p->initDecay( true );
                 p->printTree();
                 massProb = p->compMassProb();
                 if ( massProb > 0. ) {
                     massProb = 2.0;
                     EvtGenReport( EVTGEN_INFO, "EvtGen" )
                         << "Taking the minimum mass of all particles in the chain\n";
                 } else {
                     EvtGenReport( EVTGEN_INFO, "EvtGen" )
                         << "Sorry, no luck finding a valid set of masses.  This may be a pathological combo\n";
                     isOK = false;
                     break;
                 }
             }
         }
     }
 
     return isOK;
 }
 
 double EvtParticle::compMassProb()
 {
     EvtParticle* p = this;
     double mass = p->mass();
     double parMass = 0.;
     if ( p->getParent() ) {
         parMass = p->getParent()->mass();
     }
 
     int nDaug = p->getNDaug();
     double* dMasses = nullptr;
 
     int i;
     if ( nDaug > 0 ) {
         dMasses = new double[nDaug];
         for ( i = 0; i < nDaug; i++ )
             dMasses[i] = p->getDaug( i )->mass();
     }
 
     double temp = 1.0;
     temp = EvtPDL::getMassProb( p->getId(), mass, parMass, nDaug, dMasses );
 
     //If the particle already has a mass, we dont need to include
     //it in the probability calculation
     if ( ( !p->getParent() || _validP4 ) && temp > 0.0 )
         temp = 1.;
 
     delete[] dMasses;
     for ( i = 0; i < nDaug; i++ ) {
         temp *= p->getDaug( i )->compMassProb();
     }
     return temp;
 }
 
 void EvtParticle::deleteDaughters( bool keepChannel )
 {
     for ( size_t i = 0; i < _ndaug; i++ ) {
         _daug[i]->deleteTree();
     }
 
     _ndaug = 0;
     if ( !keepChannel )
         _channel = -10;
     _first = 1;
     _isInit = false;
 }
 
 void EvtParticle::deleteTree()
 {
     this->deleteDaughters();
 
     delete this;
 }
 
 EvtVector4C EvtParticle::epsParent( int i ) const
 {
     EvtVector4C temp;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" )
         << "and you have asked for the:" << i << "th polarization vector."
         << " I.e. you thought it was a"
         << " vector particle!" << endl;
     ::abort();
     return temp;
 }
 
 EvtVector4C EvtParticle::eps( int i ) const
 {
     EvtVector4C temp;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" )
         << "and you have asked for the:" << i << "th polarization vector."
         << " I.e. you thought it was a"
         << " vector particle!" << endl;
     ::abort();
     return temp;
 }
 
 EvtVector4C EvtParticle::epsParentPhoton( int i ) const
 {
     EvtVector4C temp;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "and you have asked for the:" << i
                                            << "th polarization vector of photon."
                                            << " I.e. you thought it was a"
                                            << " photon particle!" << endl;
     ::abort();
     return temp;
 }
 
 EvtVector4C EvtParticle::epsPhoton( int i ) const
 {
     EvtVector4C temp;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" )
         << "and you have asked for the:" << i
         << "th polarization vector of a photon."
         << " I.e. you thought it was a"
         << " photon particle!" << endl;
     ::abort();
     return temp;
 }
 
 EvtDiracSpinor EvtParticle::spParent( int i ) const
 {
     EvtDiracSpinor tempD;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" )
         << "and you have asked for the:" << i << "th dirac spinor."
         << " I.e. you thought it was a"
         << " Dirac particle!" << endl;
     ::abort();
     return tempD;
 }
 
 EvtDiracSpinor EvtParticle::sp( int i ) const
 {
     EvtDiracSpinor tempD;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" )
         << "and you have asked for the:" << i << "th dirac spinor."
         << " I.e. you thought it was a"
         << " Dirac particle!" << endl;
     ::abort();
     return tempD;
 }
 
 EvtDiracSpinor EvtParticle::spParentNeutrino() const
 {
     EvtDiracSpinor tempD;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "and you have asked for the "
                                            << "dirac spinor."
                                            << " I.e. you thought it was a"
                                            << " neutrino particle!" << endl;
     ::abort();
     return tempD;
 }
 
 EvtDiracSpinor EvtParticle::spNeutrino() const
 {
     EvtDiracSpinor tempD;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "and you have asked for the "
                                            << "dirac spinor."
                                            << " I.e. you thought it was a"
                                            << " neutrino particle!" << endl;
     ::abort();
     return tempD;
 }
 
 EvtTensor4C EvtParticle::epsTensorParent( int i ) const
 {
     EvtTensor4C tempC;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" )
         << "and you have asked for the:" << i << "th tensor."
         << " I.e. you thought it was a"
         << " Tensor particle!" << endl;
     ::abort();
     return tempC;
 }
 
 EvtTensor4C EvtParticle::epsTensor( int i ) const
 {
     EvtTensor4C tempC;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" )
         << "and you have asked for the:" << i << "th tensor."
         << " I.e. you thought it was a"
         << " Tensor particle!" << endl;
     ::abort();
     return tempC;
 }
 
 EvtRaritaSchwinger EvtParticle::spRSParent( int i ) const
 {
     EvtRaritaSchwinger tempD;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" )
         << "and you have asked for the:" << i << "th Rarita-Schwinger spinor."
         << " I.e. you thought it was a"
         << " RaritaSchwinger particle!" << std::endl;
     ::abort();
     return tempD;
 }
 
 EvtRaritaSchwinger EvtParticle::spRS( int i ) const
 {
     EvtRaritaSchwinger tempD;
     printParticle();
     EvtGenReport( EVTGEN_ERROR, "EvtGen" )
         << "and you have asked for the:" << i << "th Rarita-Schwinger spinor."
         << " I.e. you thought it was a"
         << " RaritaSchwinger particle!" << std::endl;
     ::abort();
     return tempD;
 }
 
 EvtVector4R EvtParticle::getP4Lab() const
 {
     EvtVector4R temp, mom;
     const EvtParticle* ptemp;
 
     temp = this->getP4();
     ptemp = this;
 
     while ( ptemp->getParent() ) {
         ptemp = ptemp->getParent();
         mom = ptemp->getP4();
         temp = boostTo( temp, mom );
     }
     return temp;
 }
 
-EvtVector4R EvtParticle::getP4LabBeforeFSR()
+EvtVector4R EvtParticle::getP4LabBeforeFSR() const
 {
     EvtVector4R temp, mom;
-    EvtParticle* ptemp;
+    const EvtParticle* ptemp;
 
     temp = this->_pBeforeFSR;
     ptemp = this;
 
     while ( ptemp->getParent() ) {
         ptemp = ptemp->getParent();
         mom = ptemp->getP4();
         temp = boostTo( temp, mom );
     }
     return temp;
 }
 
 EvtVector4R EvtParticle::getP4Restframe() const
 {
     return EvtVector4R( mass(), 0.0, 0.0, 0.0 );
 }
 
 EvtVector4R EvtParticle::get4Pos() const
 {
     EvtVector4R temp, mom;
     EvtParticle* ptemp;
 
     temp.set( 0.0, 0.0, 0.0, 0.0 );
     ptemp = getParent();
 
     if ( !ptemp ) {
         return temp;
     }
 
     temp = ( ptemp->_t / ptemp->mass() ) * ( ptemp->getP4() );
 
     while ( ptemp->getParent() ) {
         ptemp = ptemp->getParent();
         mom = ptemp->getP4();
         temp = boostTo( temp, mom );
         temp = temp + ( ptemp->_t / ptemp->mass() ) * ( ptemp->getP4() );
     }
 
     return temp;
 }
 
 EvtParticle* EvtParticle::nextIter( EvtParticle* rootOfTree )
 {
     EvtParticle* bpart;
     EvtParticle* current;
 
     current = this;
     size_t i;
 
     if ( _ndaug != 0 )
         return _daug[0];
 
     do {
         bpart = current->_parent;
         if ( !bpart ) {
             return nullptr;
         }
         i = 0;
         while ( bpart->_daug[i] != current ) {
             i++;
         }
 
         if ( bpart == rootOfTree ) {
             if ( i + 1 == bpart->_ndaug ) {
                 return nullptr;
             }
         }
 
         i++;
         current = bpart;
 
     } while ( i >= bpart->_ndaug );
 
     return bpart->_daug[i];
 }
 
 void EvtParticle::makeStdHep( EvtStdHep& stdhep, EvtSecondary& secondary,
                               EvtId* list_of_stable )
 {
     //first add particle to the stdhep list;
     stdhep.createParticle( getP4Lab(), get4Pos(), -1, -1,
                            EvtPDL::getStdHep( getId() ) );
 
     int ii = 0;
 
     //lets see if this is a longlived particle and terminate the
     //list building!
 
     while ( list_of_stable[ii] != EvtId( -1, -1 ) ) {
         if ( getId() == list_of_stable[ii] ) {
             secondary.createSecondary( 0, this );
             return;
         }
         ii++;
     }
 
     for ( size_t i = 0; i < _ndaug; i++ ) {
         stdhep.createParticle( _daug[i]->getP4Lab(), _daug[i]->get4Pos(), 0, 0,
                                EvtPDL::getStdHep( _daug[i]->getId() ) );
     }
 
     for ( size_t i = 0; i < _ndaug; i++ ) {
         _daug[i]->makeStdHepRec( 1 + i, 1 + i, stdhep, secondary, list_of_stable );
     }
     return;
 }
 
 void EvtParticle::makeStdHep( EvtStdHep& stdhep )
 {
     //first add particle to the stdhep list;
     stdhep.createParticle( getP4Lab(), get4Pos(), -1, -1,
                            EvtPDL::getStdHep( getId() ) );
 
     for ( size_t i = 0; i < _ndaug; i++ ) {
         stdhep.createParticle( _daug[i]->getP4Lab(), _daug[i]->get4Pos(), 0, 0,
                                EvtPDL::getStdHep( _daug[i]->getId() ) );
     }
 
     for ( size_t i = 0; i < _ndaug; i++ ) {
         _daug[i]->makeStdHepRec( 1 + i, 1 + i, stdhep );
     }
     return;
 }
 
 void EvtParticle::makeStdHepRec( int firstparent, int lastparent,
                                  EvtStdHep& stdhep, EvtSecondary& secondary,
                                  EvtId* list_of_stable )
 {
     int ii = 0;
 
     //lets see if this is a longlived particle and terminate the
     //list building!
 
     while ( list_of_stable[ii] != EvtId( -1, -1 ) ) {
         if ( getId() == list_of_stable[ii] ) {
             secondary.createSecondary( firstparent, this );
             return;
         }
         ii++;
     }
 
     int parent_num = stdhep.getNPart();
     for ( size_t i = 0; i < _ndaug; i++ ) {
         stdhep.createParticle( _daug[i]->getP4Lab(), _daug[i]->get4Pos(),
                                firstparent, lastparent,
                                EvtPDL::getStdHep( _daug[i]->getId() ) );
     }
 
     for ( size_t i = 0; i < _ndaug; i++ ) {
         _daug[i]->makeStdHepRec( parent_num + i, parent_num + i, stdhep,
                                  secondary, list_of_stable );
     }
     return;
 }
 
 void EvtParticle::makeStdHepRec( int firstparent, int lastparent,
                                  EvtStdHep& stdhep )
 {
     int parent_num = stdhep.getNPart();
     for ( size_t i = 0; i < _ndaug; i++ ) {
         stdhep.createParticle( _daug[i]->getP4Lab(), _daug[i]->get4Pos(),
                                firstparent, lastparent,
                                EvtPDL::getStdHep( _daug[i]->getId() ) );
     }
 
     for ( size_t i = 0; i < _ndaug; i++ ) {
         _daug[i]->makeStdHepRec( parent_num + i, parent_num + i, stdhep );
     }
     return;
 }
 
 void EvtParticle::printTreeRec( unsigned int level ) const
 {
     size_t newlevel, i;
     newlevel = level + 1;
 
     if ( _ndaug != 0 ) {
         if ( level > 0 ) {
             for ( i = 0; i < ( 5 * level ); i++ ) {
                 EvtGenReport( EVTGEN_INFO, "" ) << " ";
             }
         }
         EvtGenReport( EVTGEN_INFO, "" ) << EvtPDL::name( _id ).c_str();
         EvtGenReport( EVTGEN_INFO, "" ) << " -> ";
         for ( i = 0; i < _ndaug; i++ ) {
             EvtGenReport( EVTGEN_INFO, "" )
                 << EvtPDL::name( _daug[i]->getId() ).c_str() << " ";
         }
         for ( i = 0; i < _ndaug; i++ ) {
             EvtGenReport( EVTGEN_INFO, "" )
                 << _daug[i]->mass() << " " << _daug[i]->getP4() << " "
                 << _daug[i]->getSpinStates() << "; ";
         }
         EvtGenReport( EVTGEN_INFO, "" ) << endl;
         for ( i = 0; i < _ndaug; i++ ) {
             _daug[i]->printTreeRec( newlevel );
         }
     }
 }
 
 void EvtParticle::printTree() const
 {
     EvtGenReport( EVTGEN_INFO, "EvtGen" )
         << "This is the current decay chain" << endl;
     EvtGenReport( EVTGEN_INFO, "" )
         << "This top particle is " << EvtPDL::name( _id ).c_str() << " "
         << this->mass() << " " << this->getP4() << endl;
 
     this->printTreeRec( 0 );
     EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "End of decay chain." << endl;
 }
 
 std::string EvtParticle::treeStrRec( unsigned int level ) const
 {
     size_t newlevel, i;
     newlevel = level + 1;
 
     std::string retval = "";
 
     for ( i = 0; i < _ndaug; i++ ) {
         retval += EvtPDL::name( _daug[i]->getId() );
         if ( _daug[i]->getNDaug() > 0 ) {
             retval += " (";
             retval += _daug[i]->treeStrRec( newlevel );
             retval += ") ";
         } else {
             if ( i + 1 != _ndaug )
                 retval += " ";
         }
     }
 
     return retval;
 }
 
 std::string EvtParticle::treeStr() const
 {
     std::string retval = EvtPDL::name( _id );
     retval += " -> ";
 
     retval += treeStrRec( 0 );
 
     return retval;
 }
 
 void EvtParticle::printParticle() const
 {
     switch ( EvtPDL::getSpinType( _id ) ) {
         case EvtSpinType::SCALAR:
             EvtGenReport( EVTGEN_INFO, "EvtGen" )
                 << "This is a scalar particle:" << EvtPDL::name( _id ).c_str()
                 << "\n";
             break;
         case EvtSpinType::VECTOR:
             EvtGenReport( EVTGEN_INFO, "EvtGen" )
                 << "This is a vector particle:" << EvtPDL::name( _id ).c_str()
                 << "\n";
             break;
         case EvtSpinType::TENSOR:
             EvtGenReport( EVTGEN_INFO, "EvtGen" )
                 << "This is a tensor particle:" << EvtPDL::name( _id ).c_str()
                 << "\n";
             break;
         case EvtSpinType::DIRAC:
             EvtGenReport( EVTGEN_INFO, "EvtGen" )
                 << "This is a dirac particle:" << EvtPDL::name( _id ).c_str()
                 << "\n";
             break;
         case EvtSpinType::PHOTON:
             EvtGenReport( EVTGEN_INFO, "EvtGen" )
                 << "This is a photon:" << EvtPDL::name( _id ).c_str() << "\n";
             break;
         case EvtSpinType::NEUTRINO:
             EvtGenReport( EVTGEN_INFO, "EvtGen" )
                 << "This is a neutrino:" << EvtPDL::name( _id ).c_str() << "\n";
             break;
         case EvtSpinType::STRING:
             EvtGenReport( EVTGEN_INFO, "EvtGen" )
                 << "This is a string:" << EvtPDL::name( _id ).c_str() << "\n";
             break;
         default:
             EvtGenReport( EVTGEN_INFO, "EvtGen" )
                 << "Unknown particle type in EvtParticle::printParticle()"
                 << endl;
             break;
     }
     EvtGenReport( EVTGEN_INFO, "EvtGen" )
         << "Number of daughters:" << _ndaug << "\n";
 }
 
 void init_vector( EvtParticle** part )
 {
     *part = new EvtVectorParticle;
 }
 
 void init_scalar( EvtParticle** part )
 {
     *part = new EvtScalarParticle;
 }
 
 void init_tensor( EvtParticle** part )
 {
     *part = new EvtTensorParticle;
 }
 
 void init_dirac( EvtParticle** part )
 {
     *part = new EvtDiracParticle;
 }
 
 void init_photon( EvtParticle** part )
 {
     *part = new EvtPhotonParticle;
 }
 
 void init_neutrino( EvtParticle** part )
 {
     *part = new EvtNeutrinoParticle;
 }
 
 void init_string( EvtParticle** part )
 {
     *part = new EvtStringParticle;
 }
 
 double EvtParticle::initializePhaseSpace( unsigned int numdaughter,
                                           EvtId* daughters,
                                           bool forceDaugMassReset,
                                           double poleSize, int whichTwo1,
                                           int whichTwo2 )
 {
     double m_b;
     unsigned int i;
     //lange
     //  this->makeDaughters(numdaughter,daughters);
 
     static EvtVector4R p4[100];
     static double mass[100];
 
     m_b = this->mass();
 
     //lange - Jan2,2002 - Need to check to see if the daughters of the parent
     // have changed. If so, delete them and start over.
     //EvtGenReport(EVTGEN_INFO,"EvtGen") << "the parent is\n";
     //if ( this->getParent() ) {
     //  if ( this->getParent()->getParent() ) this->getParent()->getParent()->printTree();
     //    this->getParent()->printTree();
     //}
     //EvtGenReport(EVTGEN_INFO,"EvtGen") << "and this is\n";
     //if ( this) this->printTree();
     bool resetDaughters = false;
 
     if ( numdaughter != this->getNDaug() && this->getNDaug() > 0 )
         resetDaughters = true;
     if ( numdaughter == this->getNDaug() )
         for ( i = 0; i < numdaughter; i++ ) {
             if ( this->getDaug( i )->getId() != daughters[i] )
                 resetDaughters = true;
             //EvtGenReport(EVTGEN_INFO,"EvtGen") << EvtPDL::name(this->getDaug(i)->getId())
             //		    << " " << EvtPDL::name(daughters[i]) << endl;
         }
 
     if ( resetDaughters || forceDaugMassReset ) {
         bool t1 = true;
         //but keep the decay channel of the parent.
         this->deleteDaughters( t1 );
         this->makeDaughters( numdaughter, daughters );
         bool massTreeOK = this->generateMassTree();
         if ( massTreeOK == false ) {
             return 0.0;
         }
     }
 
     double weight = 0.;
     for ( i = 0; i < numdaughter; i++ ) {
         mass[i] = this->getDaug( i )->mass();
     }
 
     if ( poleSize < -0.1 ) {
         //special case to enforce 4-momentum conservation in 1->1 decays
         if ( numdaughter == 1 ) {
             this->getDaug( 0 )->init( daughters[0],
                                       EvtVector4R( m_b, 0.0, 0.0, 0.0 ) );
         } else {
             EvtGenKine::PhaseSpace( numdaughter, mass, p4, m_b );
             for ( i = 0; i < numdaughter; i++ ) {
                 this->getDaug( i )->init( daughters[i], p4[i] );
             }
         }
     } else {
         if ( numdaughter != 3 ) {
             EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                 << "Only can generate pole phase space "
                 << "distributions for 3 body final states" << endl
                 << "Will terminate." << endl;
             ::abort();
         }
         bool ok = false;
         if ( ( whichTwo1 == 1 && whichTwo2 == 0 ) ||
              ( whichTwo1 == 0 && whichTwo2 == 1 ) ) {
             weight = EvtGenKine::PhaseSpacePole( m_b, mass[0], mass[1], mass[2],
                                                  poleSize, p4 );
             this->getDaug( 0 )->init( daughters[0], p4[0] );
             this->getDaug( 1 )->init( daughters[1], p4[1] );
             this->getDaug( 2 )->init( daughters[2], p4[2] );
             ok = true;
         }
         if ( ( whichTwo1 == 1 && whichTwo2 == 2 ) ||
              ( whichTwo1 == 2 && whichTwo2 == 1 ) ) {
             weight = EvtGenKine::PhaseSpacePole( m_b, mass[2], mass[1], mass[0],
                                                  poleSize, p4 );
             this->getDaug( 0 )->init( daughters[0], p4[2] );
             this->getDaug( 1 )->init( daughters[1], p4[1] );
             this->getDaug( 2 )->init( daughters[2], p4[0] );
             ok = true;
         }
         if ( ( whichTwo1 == 0 && whichTwo2 == 2 ) ||
              ( whichTwo1 == 2 && whichTwo2 == 0 ) ) {
             weight = EvtGenKine::PhaseSpacePole( m_b, mass[1], mass[0], mass[2],
                                                  poleSize, p4 );
             this->getDaug( 0 )->init( daughters[0], p4[1] );
             this->getDaug( 1 )->init( daughters[1], p4[0] );
             this->getDaug( 2 )->init( daughters[2], p4[2] );
             ok = true;
         }
         if ( !ok ) {
             EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                 << "Invalid pair of particle to generate a pole dist "
                 << whichTwo1 << " " << whichTwo2 << endl
                 << "Will terminate." << endl;
             ::abort();
         }
     }
 
     return weight;
 }
 
 void EvtParticle::makeDaughters( unsigned int ndaugstore,
                                  std::vector<EvtId> idVector )
 {
     // Convert the STL vector method to use the array method for now, since the
     // array method pervades most of the EvtGen code...
 
     unsigned int nVector = idVector.size();
     if ( nVector < ndaugstore ) {
         EvtGenReport( EVTGEN_ERROR, "EvtGen" )
             << "Asking to make " << ndaugstore << " daughters when there "
             << "are only " << nVector << " EvtId values available" << endl;
         return;
     }
 
     EvtId* idArray = new EvtId[ndaugstore];
     unsigned int i;
     for ( i = 0; i < ndaugstore; i++ ) {
         idArray[i] = idVector[i];
     }
 
     this->makeDaughters( ndaugstore, idArray );
 
     delete[] idArray;
 }
 
 void EvtParticle::makeDaughters( unsigned int ndaugstore, EvtId* id )
 {
     unsigned int i;
     if ( _channel < 0 ) {
         setChannel( 0 );
     }
     EvtParticle* pdaug;
     if ( _ndaug != 0 ) {
         if ( _ndaug != ndaugstore ) {
             EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                 << "Asking to make a different number of "
                 << "daughters than what was previously created." << endl;
             EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                 << "Original parent:" << EvtPDL::name( _id ) << endl;
             for ( size_t i = 0; i < _ndaug; i++ ) {
                 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                     << "Original daugther:"
                     << EvtPDL::name( getDaug( i )->getId() ) << endl;
             }
             for ( size_t i = 0; i < ndaugstore; i++ ) {
                 EvtGenReport( EVTGEN_ERROR, "EvtGen" )
                     << "New Daug:" << EvtPDL::name( id[i] ) << endl;
             }
             EvtGenReport( EVTGEN_ERROR, "EvtGen" ) << "Will terminate." << endl;
             ::abort();
         }
     } else {
         for ( i = 0; i < ndaugstore; i++ ) {
             pdaug = EvtParticleFactory::particleFactory(
                 EvtPDL::getSpinType( id[i] ) );
             pdaug->setId( id[i] );
             pdaug->addDaug( this );
         }
 
     }    //else
 }    //makeDaughters
 
 void EvtParticle::setDecayProb( double prob )
 {
     if ( _decayProb == nullptr )
         _decayProb = new double;
     *_decayProb = prob;
 }
 
 std::string EvtParticle::getName() const
 {
     std::string theName = _id.getName();
     return theName;
 }
 
 int EvtParticle::getAttribute( std::string attName ) const
 {
     // Retrieve the attribute integer if the name exists.
     // Otherwise, simply return 0
 
     int attValue = 0;
 
     EvtAttIntMap::const_iterator mapIter;
 
     if ( ( mapIter = _intAttributes.find( attName ) ) != _intAttributes.end() ) {
         attValue = mapIter->second;
     }
 
     return attValue;
 }
 
 double EvtParticle::getAttributeDouble( std::string attName ) const
 {
     // Retrieve the attribute double if the name exists.
     // Otherwise, simply return 0.0
 
     double attValue = 0.0;
 
     EvtAttDblMap::const_iterator mapIter;
 
     if ( ( mapIter = _dblAttributes.find( attName ) ) != _dblAttributes.end() ) {
         attValue = mapIter->second;
     }
 
     return attValue;
 }
diff --git a/test/jsonFiles/D_DALITZ__D+_Kpipi+_yesPhotos.json b/test/jsonFiles/D_DALITZ__D+_Kpipi+_yesPhotos.json
new file mode 100644
index 0000000..2f7ccf6
--- /dev/null
+++ b/test/jsonFiles/D_DALITZ__D+_Kpipi+_yesPhotos.json
@@ -0,0 +1,31 @@
+{
+    "parent" : "D+",
+    "daughters" : ["K-", "pi+", "pi+"],
+    "models" : ["D_DALITZ"],
+    "parameters" : [[]],
+    "do_conjugate_decay" : [true],
+    "extras" : ["yesPhotos"],
+    "events" : 10000,
+    "histograms" : [
+	{"variable" : "prob", "title" : "Probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 0.8 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.8},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta/#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+	{"variable" : "mass", "title" : "m(K^{-} #pi^{+})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.5, "xmax" : 1.8},
+        {"variable" : "mass", "title" : "m(K^{-} #pi^{+})", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 0.5, "xmax" : 1.8},
+        {"variable" : "mass", "title" : "m(#pi^{+} #pi^{+})", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.2, "xmax" : 1.5},
+        {"variable" : "massSq", "title" : "Dalitz m^{2}(K^{-} #pi^{+}) vs m^{2}(K^{-} #pi^{+})", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.3, "xmax" : 3.2, "variableY" : "massSq", "d1Y" : 1, "d2Y" : 3, "nbinsY" : 50, "ymin" : 0.3, "ymax" : 3.2},
+        {"variable" : "massSq", "title" : "Dalitz m^{2}(#pi^{+} #pi^{+}) vs m^{2}(K^{-} #pi^{+})", "d1" : 2, "d2" : 3, "nbins" : 50, "xmin" : 0.0, "xmax" : 2.0, "variableY" : "massSq", "d1Y" : 1, "d2Y" : 2, "nbinsY" : 50, "ymin" : 0.3, "ymax" : 3.2},
+        {"variable" : "mPrime", "title" : "sqDalitz", "d1" : 1, "d2" : 2, "nbins" : 25, "xmin" : 0.0, "xmax" : 1.00, "variableY" : "thetaPrime", "d1Y" : 1, "d2Y" : 2, "nbinsY" : 25, "ymin" : 0.0, "ymax" : 1.0},
+	{"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosHel", "title" : "cosHel23", "d1" : 2, "d2" : 3, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosHel", "title" : "cosHel13", "d1" : 1, "d2" : 3, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0}
+    ],
+    "outfile" : "D_DALITZ__D+_Kpipi+_yesPhotos.root",
+    "reference" : "RefD_DALITZ__D+_Kpipi+_yesPhotos.root"
+}
diff --git a/test/jsonFiles/PHSP__Bd_KKpipi_yesPhotos.json b/test/jsonFiles/PHSP__Bd_KKpipi_yesPhotos.json
new file mode 100644
index 0000000..17148ce
--- /dev/null
+++ b/test/jsonFiles/PHSP__Bd_KKpipi_yesPhotos.json
@@ -0,0 +1,33 @@
+{
+    "parent" : "B0",
+    "daughters" : ["K+", "K-", "pi+", "pi-"],
+    "models" : ["PHSP"],
+    "parameters" : [[]],
+    "do_conjugate_decay" : [true],    
+    "extras" : ["yesPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.5 },
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.0 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 2.0},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta/#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+	{"variable" : "parMass", "title" : "m(B^{0})", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+	{"variable" : "mass", "title" : "m(K^{+}K^{-})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+        {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},	
+        {"variable" : "E", "title" : "E(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "E", "title" : "E(K^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(K^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(K^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "cosTheta", "title" : "cosTheta(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(K^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+    ],
+    "outfile" : "PHSP__Bd_KKpipi_yesPhotos.root",
+    "reference" : "RefPHSP__Bd_KKpipi_yesPhotos.root"
+}
diff --git a/test/jsonFiles/PHSP__Bu_Kpi0_yesPhotos.json b/test/jsonFiles/PHSP__Bu_Kpi0_yesPhotos.json
new file mode 100644
index 0000000..0e24d92
--- /dev/null
+++ b/test/jsonFiles/PHSP__Bu_Kpi0_yesPhotos.json
@@ -0,0 +1,33 @@
+{
+    "parent" : "B+",
+    "daughters" : ["K+", "pi0"],
+    "models" : ["PHSP"],
+    "parameters" : [[]],
+    "do_conjugate_decay" : [true],    
+    "extras" : ["yesPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.5 },
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.0 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 2.0},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta/#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+	{"variable" : "parMass", "title" : "m(B^{+})", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+	{"variable" : "mass", "title" : "m(K^{+}#pi^{-})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+        {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},	
+        {"variable" : "E", "title" : "E(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "E", "title" : "E(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "cosTheta", "title" : "cosTheta(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+    ],
+    "outfile" : "PHSP__Bu_Kpi0_yesPhotos.root",
+    "reference" : "RefPHSP__Bu_Kpi0_yesPhotos.root"
+}
diff --git a/test/jsonFiles/PHSP__D0_KK_yesPhotos.json b/test/jsonFiles/PHSP__D0_KK_yesPhotos.json
new file mode 100644
index 0000000..9fb1ef0
--- /dev/null
+++ b/test/jsonFiles/PHSP__D0_KK_yesPhotos.json
@@ -0,0 +1,33 @@
+{
+    "parent" : "D0",
+    "daughters" : ["K+", "K-"],
+    "models" : ["PHSP"],
+    "parameters" : [[]],
+    "do_conjugate_decay" : [true],    
+    "extras" : ["yesPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.5 },
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 0.8 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.8},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta/#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+	{"variable" : "parMass", "title" : "m(D^{0})", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+	{"variable" : "mass", "title" : "m(K^{+}K^{-})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+        {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},	
+        {"variable" : "E", "title" : "E(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "E", "title" : "E(K^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(K^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(K^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "cosTheta", "title" : "cosTheta(K^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(K^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+    ],
+    "outfile" : "PHSP__D0_KK_yesPhotos.root",
+    "reference" : "RefPHSP__D0_KK_yesPhotos.root"
+}
diff --git a/test/jsonFiles/PHSP__D0_Kpi_yesPhotos.json b/test/jsonFiles/PHSP__D0_Kpi_yesPhotos.json
new file mode 100644
index 0000000..ab5d1dc
--- /dev/null
+++ b/test/jsonFiles/PHSP__D0_Kpi_yesPhotos.json
@@ -0,0 +1,33 @@
+{
+    "parent" : "D0",
+    "daughters" : ["K-", "pi+"],
+    "models" : ["PHSP"],
+    "parameters" : [[]],
+    "do_conjugate_decay" : [true],    
+    "extras" : ["yesPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.5 },
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 0.8 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.8},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta/#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+	{"variable" : "parMass", "title" : "m(D^{0})", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+	{"variable" : "mass", "title" : "m(K^{-}#pi^{+})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+        {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},	
+        {"variable" : "E", "title" : "E(K^{-})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "E", "title" : "E(#pi^{+})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(K^{-})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#pi^{+})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(K^{-})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#pi^{+})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "cosTheta", "title" : "cosTheta(K^{-})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(#pi^{+})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+    ],
+    "outfile" : "PHSP__D0_Kpi_yesPhotos.root",
+    "reference" : "RefPHSP__D0_Kpi_yesPhotos.root"
+}
diff --git a/test/jsonFiles/PHSP__D0_pipi_yesPhotos.json b/test/jsonFiles/PHSP__D0_pipi_yesPhotos.json
new file mode 100644
index 0000000..a4ae339
--- /dev/null
+++ b/test/jsonFiles/PHSP__D0_pipi_yesPhotos.json
@@ -0,0 +1,33 @@
+{
+    "parent" : "D0",
+    "daughters" : ["pi+", "pi-"],
+    "models" : ["PHSP"],
+    "parameters" : [[]],
+    "do_conjugate_decay" : [true],    
+    "extras" : ["yesPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.5 },
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 0.8 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.8},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta/#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+	{"variable" : "parMass", "title" : "m(D^{0})", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+	{"variable" : "mass", "title" : "m(#pi^{+}#pi^{-})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+        {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},	
+        {"variable" : "E", "title" : "E(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "E", "title" : "E(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "cosTheta", "title" : "cosTheta(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+    ],
+    "outfile" : "PHSP__D0_pipi_yesPhotos.root",
+    "reference" : "RefPHSP__D0_pipi_yesPhotos.root"
+}
diff --git a/test/jsonFiles/PHSP__Lambda0_ppi_yesPhotos.json b/test/jsonFiles/PHSP__Lambda0_ppi_yesPhotos.json
new file mode 100644
index 0000000..99ecd64
--- /dev/null
+++ b/test/jsonFiles/PHSP__Lambda0_ppi_yesPhotos.json
@@ -0,0 +1,33 @@
+{
+    "parent" : "Lambda0",
+    "daughters" : ["p+", "pi-"],
+    "models" : ["PHSP"],
+    "parameters" : [[]],
+    "do_conjugate_decay" : [true],    
+    "extras" : ["yesPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.5 },
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 6, "xmin" : 0, "xmax" : 6.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 50, "xmin" : 0, "xmax" : 0.05 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.05},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta/#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+	{"variable" : "parMass", "title" : "m(#Lambda)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+	{"variable" : "mass", "title" : "m(p^{+}#pi^{-})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+        {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},	
+        {"variable" : "E", "title" : "E(p^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "E", "title" : "E(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(p^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(p^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "cosTheta", "title" : "cosTheta(p^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+    ],
+    "outfile" : "PHSP__Lambda0_ppi_yesPhotos.root",
+    "reference" : "RefPHSP__Lambda0_ppi_yesPhotos.root"
+}
diff --git a/test/jsonFiles/PHSP__Lambdab0_ppi_yesPhotos.json b/test/jsonFiles/PHSP__Lambdab0_ppi_yesPhotos.json
new file mode 100644
index 0000000..6fc4c1e
--- /dev/null
+++ b/test/jsonFiles/PHSP__Lambdab0_ppi_yesPhotos.json
@@ -0,0 +1,33 @@
+{
+    "parent" : "Lambda_b0",
+    "daughters" : ["p+", "pi-"],
+    "models" : ["PHSP"],
+    "parameters" : [[]],
+    "do_conjugate_decay" : [true],    
+    "extras" : ["yesPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.5 },
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.0 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 2.0},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta/#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+	{"variable" : "parMass", "title" : "m(#Lambda_{b})", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+	{"variable" : "mass", "title" : "m(p^{+}#pi^{-})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+        {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},	
+        {"variable" : "E", "title" : "E(p^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "E", "title" : "E(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(p^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(p^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "cosTheta", "title" : "cosTheta(p^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+    ],
+    "outfile" : "PHSP__Lambdab0_ppi_yesPhotos.root",
+    "reference" : "RefPHSP__Lambdab0_ppi_yesPhotos.root"
+}
diff --git a/test/jsonFiles/SVP_HELAMP=VSS__Bu_Kstgamma_Kspi_yesPhotos.json b/test/jsonFiles/SVP_HELAMP=VSS__Bu_Kstgamma_Kspi_yesPhotos.json
new file mode 100644
index 0000000..f1e22dc
--- /dev/null
+++ b/test/jsonFiles/SVP_HELAMP=VSS__Bu_Kstgamma_Kspi_yesPhotos.json
@@ -0,0 +1,42 @@
+{
+    "parent" : "B+",
+    "daughters" : ["K*+", "gamma"],
+    "grand_daughters" : [["K_S0", "pi+"], []],
+    "models" : ["SVP_HELAMP", "VSS", ""],
+    "parameters" : [["1.0", "0.0", "1.0", "0.0"], [], []],
+    "extras" : ["yesPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "Prob",  "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.0 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 2.0},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "mass", "title" : "M(#it{B}^{+})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 5.0, "xmax" : 5.28},
+        {"variable" : "p", "title" : "p(#it{K}*^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.0, "xmax" : 3.0},
+        {"variable" : "pSq", "title" : "p^{2}(#it{K}*^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 2.0, "xmax" : 7.0},
+        {"variable" : "pz", "title" : "pz(#it{K}*^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -3.0, "xmax" : 3.0},
+        {"variable" : "p", "title" : "p(#it{#gamma})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.0, "xmax" : 3.0},
+        {"variable" : "pSq", "title" : "p^{2}(#it{#gamma})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 2.0, "xmax" : 7.0},
+        {"variable" : "pz", "title" : "pz(#it{#gamma})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -3.0, "xmax" : 3.0},
+        {"variable" : "cosHel", "title" : "cosHel(#it{K}*^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosHelDiTau", "title" : "cosHel(#it{K}*^{+})", "d1" : 1, "d2" : 1, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(#it{K}*^{+})", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(#it{#gamma})", "d1" : 2, "d2" : 0, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "phi", "title" : "phi(#it{K}*^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
+        {"variable" : "phi", "title" : "phi(#it{#gamma})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
+        {"variable" : "mass_daug1", "title" : "M(#it{K}*^{+})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.6, "xmax" : 1.5},
+        {"variable" : "pLab_daug1", "title" : "pLab(#it{K}^{0}_{S})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 3.0},
+        {"variable" : "cosTheta_daug1", "title" : "cosTheta(#it{K}^{0}_{S})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "phi_daug1", "title" : "phi(#it{K}^{0}_{S})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0},
+        {"variable" : "pLab_daug1", "title" : "pLab(#it{#pi}^{0})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 3.0},
+        {"variable" : "cosTheta_daug1", "title" : "cosTheta(#it{#pi}^{0})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "phi_daug1", "title" : "phi(#it{#pi}^{0})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -180.0, "xmax" : 180.0}
+    ],
+    "outfile" : "SVP_HELAMP=VSS__Bu_Kstgamma_Kspi_yesPhotos.root",
+    "reference" : "RefSVP_HELAMP=VSS__Bu_Kstgamma_Kspi_yesPhotos.root"
+}
diff --git a/test/jsonFiles/TAULNUNU__tau_enunu_yesPhotos.json b/test/jsonFiles/TAULNUNU__tau_enunu_yesPhotos.json
new file mode 100644
index 0000000..0334ee7
--- /dev/null
+++ b/test/jsonFiles/TAULNUNU__tau_enunu_yesPhotos.json
@@ -0,0 +1,32 @@
+{
+    "parent" : "tau+",
+    "daughters" : ["e+", "nu_e", "anti-nu_tau"],
+    "models" : ["TAULNUNU"],
+    "parameters" : [[]],
+    "extras" : ["yesPhotos"],
+    "events" : 100000,
+    "histograms" : [
+	{"variable" : "prob", "title" : "Prob",  "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.8, "xmax" : 1.2},
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 8, "xmin" : 0, "xmax" : 8.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 1.0 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta/#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+	{"variable" : "pSumSq", "title" : "qSq(e#nu_{e})", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : 0.0, "xmax" : 3.0},
+	{"variable" : "cosHel", "title" : "cosHel(e)", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "E_over_Eparent", "title" : "E(e)/E(tau)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.5},
+        {"variable" : "E_over_Eparent", "title" : "E(nue)/E(tau)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.5},
+        {"variable" : "E_over_Eparent", "title" : "E(nutau_{2})/E(tau)", "d1" : 3, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.5},
+        {"variable" : "cosTheta", "title" : "cosTheta(e)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(nue)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(nutau_{2})", "d1" : 3, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "mass", "title" : "m(e nue)", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.25, "xmax" : 2.0},
+        {"variable" : "mass", "title" : "m(e nutau_{2})", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 0.25, "xmax" : 2.0},
+        {"variable" : "mass", "title" : "m(nue nutau_{2})", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.25, "xmax" : 2.0}
+    ],
+    "outfile" : "TAULNUNU__tau_enunu_yesPhotos.root",
+    "reference" : "RefTAULNUNU__tau_enunu_yesPhotos.root"
+}
diff --git a/test/jsonFiles/TAUOLA__tau_pipipinu_yesPhotos.json b/test/jsonFiles/TAUOLA__tau_pipipinu_yesPhotos.json
new file mode 100644
index 0000000..0a33ebd
--- /dev/null
+++ b/test/jsonFiles/TAUOLA__tau_pipipinu_yesPhotos.json
@@ -0,0 +1,40 @@
+{
+    "parent" : "tau+",
+    "daughters" : [" "],
+    "models" : ["TAUOLA"],
+    "parameters" : [["5"]],
+    "extras" : ["Define TauolaCurrentOption 1", 
+                "Define TauolaBR1 1.0",
+                "yesPhotos"],
+    "events" : 100000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "Prob",  "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.8, "xmax" : 1.2},
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 6, "xmin" : 0, "xmax" : 6.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 0.8 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.8},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta/#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "E_over_Eparent", "title" : "E(nutau_{2})/E(tau)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.5},
+        {"variable" : "E_over_Eparent", "title" : "E(pi+_{1})/E(tau)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.5},
+        {"variable" : "E_over_Eparent", "title" : "E(pi+_{2})/E(tau)", "d1" : 3, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.5},
+        {"variable" : "E_over_Eparent", "title" : "E(pi-)/E(tau)", "d1" : 4, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 0.5},
+        {"variable" : "cosTheta", "title" : "cosTheta(nutau_{2})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(pi+_{1})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(pi+_{2})", "d1" : 3, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(pi-)", "d1" : 4, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosHel", "title" : "cosHel(nutau_{2})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "mass", "title" : "m(nutau_{2} pi+_{1})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 0.25, "xmax" : 2.0},
+        {"variable" : "mass", "title" : "m(nutau_{2} pi+_{2})", "d1" : 1, "d2" : 3, "nbins" : 100, "xmin" : 0.25, "xmax" : 2.0},
+        {"variable" : "mass", "title" : "m(nutau_{2} pi-)", "d1" : 1, "d2" : 4, "nbins" : 100, "xmin" : 0.25, "xmax" : 2.0},
+        {"variable" : "mass", "title" : "m(pi+_{1} pi+_{2})", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.25, "xmax" : 2.0},
+        {"variable" : "mass", "title" : "m(pi+_{1} pi-)", "d1" : 2, "d2" : 4, "nbins" : 100, "xmin" : 0.25, "xmax" : 2.0},
+        {"variable" : "mass", "title" : "m(pi+_{2} pi-)", "d1" : 3, "d2" : 4, "nbins" : 100, "xmin" : 0.25, "xmax" : 2.0},
+        {"variable" : "mass3", "title" : "m(pi+_{1} pi+_{2} pi-)", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : 0.2, "xmax" : 1.8},
+        {"variable" : "cosTheta3", "title" : "cosTheta3pi", "d1" : 2, "d2" : 3, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+    ],
+    "outfile" : "TAUOLA__tau_pipipinu_yesPhotos.root",
+    "reference" : "RefTAUOLA__tau_pipipinu_yesPhotos.root"
+}
diff --git a/test/jsonFiles/VLL__Jpsi_ee_yesPhotos.json b/test/jsonFiles/VLL__Jpsi_ee_yesPhotos.json
index 6fe069a..dafc540 100644
--- a/test/jsonFiles/VLL__Jpsi_ee_yesPhotos.json
+++ b/test/jsonFiles/VLL__Jpsi_ee_yesPhotos.json
@@ -1,26 +1,32 @@
 {
     "parent" : "J/psi",
     "daughters" : ["e+", "e-"],
     "models" : ["VLL"],
     "parameters" : [[]],
     "extras" : ["yesPhotos"],
     "events" : 100000,
     "histograms" : [
         {"variable" : "prob", "title" : "probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.5 },
-        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons(J/#psi)", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
-        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy(J/#psi)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.0 },  
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.0 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 2.0},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
 	{"variable" : "parMass", "title" : "m(J/#psi)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
 	{"variable" : "mass", "title" : "m(e^{+}e^{-})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
         {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},	
-        {"variable" : "E", "title" : "E(e+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
-        {"variable" : "E", "title" : "E(e-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
-        {"variable" : "pz", "title" : "p_{z}(e+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
-        {"variable" : "pz", "title" : "p_{z}(e-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
-        {"variable" : "pLab", "title" : "pLab(e+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
-        {"variable" : "pLab", "title" : "pLab(e-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
-        {"variable" : "cosTheta", "title" : "cosTheta(e+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
-        {"variable" : "cosTheta", "title" : "cosTheta(e-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+        {"variable" : "E", "title" : "E(e^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "E", "title" : "E(e^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(e^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(e^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(e^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(e^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "cosTheta", "title" : "cosTheta(e^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(e^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
     ],
     "outfile" : "VLL__Jpsi_ee_yesPhotos.root",
     "reference" : "RefVLL__Jpsi_ee_yesPhotos.root"
 }
diff --git a/test/jsonFiles/VLL__Jpsi_mumu_yesPhotos.json b/test/jsonFiles/VLL__Jpsi_mumu_yesPhotos.json
index 5c315c8..a4c2354 100644
--- a/test/jsonFiles/VLL__Jpsi_mumu_yesPhotos.json
+++ b/test/jsonFiles/VLL__Jpsi_mumu_yesPhotos.json
@@ -1,26 +1,32 @@
 {
     "parent" : "J/psi",
     "daughters" : ["mu+", "mu-"],
     "models" : ["VLL"],
     "parameters" : [[]],
     "extras" : ["yesPhotos"],
     "events" : 100000,
     "histograms" : [
         {"variable" : "prob", "title" : "probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.5 },
-        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons(J/#psi)", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
-        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy(J/#psi)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.0 },  
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.0 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 2.0},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
 	{"variable" : "parMass", "title" : "m(J/#psi)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
 	{"variable" : "mass", "title" : "m(#mu^{+}#mu^{-})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
         {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},	
-        {"variable" : "E", "title" : "E(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
-        {"variable" : "E", "title" : "E(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
-        {"variable" : "pz", "title" : "p_{z}(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
-        {"variable" : "pz", "title" : "p_{z}(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
-        {"variable" : "pLab", "title" : "pLab(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
-        {"variable" : "pLab", "title" : "pLab(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
-        {"variable" : "cosTheta", "title" : "cosTheta(mu+)", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
-        {"variable" : "cosTheta", "title" : "cosTheta(mu-)", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+        {"variable" : "E", "title" : "E(#mu^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "E", "title" : "E(#mu^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#mu^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#mu^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#mu^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#mu^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "cosTheta", "title" : "cosTheta(#mu^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(#mu^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
     ],
     "outfile" : "VLL__Jpsi_mumu_yesPhotos.root",
     "reference" : "RefVLL__Jpsi_mumu_yesPhotos.root"
 }
diff --git a/test/jsonFiles/VSS__rho0_pipi_yesPhotos.json b/test/jsonFiles/VSS__rho0_pipi_yesPhotos.json
new file mode 100644
index 0000000..c2375f7
--- /dev/null
+++ b/test/jsonFiles/VSS__rho0_pipi_yesPhotos.json
@@ -0,0 +1,32 @@
+{
+    "parent" : "rho0",
+    "daughters" : ["pi+", "pi-"],
+    "models" : ["VSS"],
+    "parameters" : [[]],
+    "extras" : ["yesPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.5 },
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 1.0 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+	{"variable" : "parMass", "title" : "m(#rho^{+})", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+	{"variable" : "mass", "title" : "m(#pi^{+}#pi^{-})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+        {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},	
+        {"variable" : "E", "title" : "E(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "E", "title" : "E(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "cosTheta", "title" : "cosTheta(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(#pi^{-})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+    ],
+    "outfile" : "VSS__rho0_pipi_yesPhotos.root",
+    "reference" : "RefVSS__rho0_pipi_yesPhotos.root"
+}
diff --git a/test/jsonFiles/VSS__rho_pipi0_yesPhotos.json b/test/jsonFiles/VSS__rho_pipi0_yesPhotos.json
new file mode 100644
index 0000000..ff86dd3
--- /dev/null
+++ b/test/jsonFiles/VSS__rho_pipi0_yesPhotos.json
@@ -0,0 +1,32 @@
+{
+    "parent" : "rho+",
+    "daughters" : ["pi+", "pi0"],
+    "models" : ["VSS"],
+    "parameters" : [[]],
+    "extras" : ["yesPhotos"],
+    "events" : 10000,
+    "histograms" : [
+        {"variable" : "prob", "title" : "probability", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 2.5 },
+        {"variable" : "nFSRPhotons", "title" : "nFSRPhotons", "d1" : 0, "d2" : 0, "nbins" : 10, "xmin" : 0, "xmax" : 10.0 },  
+        {"variable" : "totalFSREnergy", "title" : "totalFSREnergy", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0, "xmax" : 1.0 },  
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 0.0, "xmax" : 1.0},
+        {"variable" : "E_FSRPhotons", "title" : "FSR E(#gamma) cut off", "d1" : 1, "d2" : 0, "nbins" : 50, "xmin" : 0.0, "xmax" : 0.5e-6},
+        {"variable" : "cosTheta_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+        {"variable" : "cosTheta_EnergyWeight_FSRPhotons", "title" : "FSR cosTheta(#gamma)", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+	{"variable" : "openingAngle_EnergyWeight_FSRPhotons", "title" : "FSR opening angle", "d1" : 0, "d2" : 0, "nbins" : 90, "xmin" : 0, "xmax" : 180},
+	{"variable" : "parMass", "title" : "m(#rho^{+})", "d1" : 0, "d2" : 0, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+	{"variable" : "mass", "title" : "m(#pi^{+}#pi^{0})", "d1" : 1, "d2" : 2, "nbins" : 100, "xmin" : 3.0967, "xmax" : 3.09695},
+        {"variable" : "cosHel", "title" : "cosHel12", "d1" : 1, "d2" : 2, "nbins" : 50, "xmin" : -1.0, "xmax" : 1.0},	
+        {"variable" : "E", "title" : "E(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "E", "title" : "E(#pi^{0})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pz", "title" : "p_{z}(#pi^{0})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.6, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "pLab", "title" : "pLab(#pi^{0})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : 1.2, "xmax" : 1.6},
+        {"variable" : "cosTheta", "title" : "cosTheta(#pi^{+})", "d1" : 1, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0},
+        {"variable" : "cosTheta", "title" : "cosTheta(#pi^{0})", "d1" : 2, "d2" : 0, "nbins" : 100, "xmin" : -1.0, "xmax" : 1.0}
+    ],
+    "outfile" : "VSS__rho_pipi0_yesPhotos.root",
+    "reference" : "RefVSS__rho_pipi0_yesPhotos.root"
+}
diff --git a/test/testDecayModel.cc b/test/testDecayModel.cc
index 7f467ea..4d54cee 100644
--- a/test/testDecayModel.cc
+++ b/test/testDecayModel.cc
@@ -1,1331 +1,1413 @@
 
 /***********************************************************************
 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
 *                                                                      *
 * This file is part of EvtGen.                                         *
 *                                                                      *
 * EvtGen is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU General Public License as published by *
 * the Free Software Foundation, either version 3 of the License, or    *
 * (at your option) any later version.                                  *
 *                                                                      *
 * EvtGen is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
 * GNU General Public License for more details.                         *
 *                                                                      *
 * You should have received a copy of the GNU General Public License    *
 * along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
 ***********************************************************************/
 
 #include "testDecayModel.hh"
 
 #include "EvtGen/EvtGen.hh"
 
 #include "EvtGenBase/EvtAbsRadCorr.hh"
 #include "EvtGenBase/EvtConst.hh"
 #include "EvtGenBase/EvtDecayBase.hh"
 #include "EvtGenBase/EvtId.hh"
 #include "EvtGenBase/EvtKine.hh"
 #include "EvtGenBase/EvtMTRandomEngine.hh"
 #include "EvtGenBase/EvtPDL.hh"
 #include "EvtGenBase/EvtParticle.hh"
 #include "EvtGenBase/EvtParticleFactory.hh"
 #include "EvtGenBase/EvtRandom.hh"
 #include "EvtGenBase/EvtVector4R.hh"
 
 #ifdef EVTGEN_EXTERNAL
 #include "EvtGenExternal/EvtExternalGenList.hh"
 #endif
 
 #include <fstream>
 #include <iostream>
 #include <list>
 #include <memory>
 
 using nlohmann::json;
 
 TestDecayModel::TestDecayModel( const json& config ) : m_config{ config }
 {
 }
 
 bool TestDecayModel::checkMandatoryFields()
 {
     const std::array<std::string, 8> mandatoryFields{
         "parent",  "daughters", "models",    "parameters",
         "outfile", "events",    "reference", "histograms" };
     const std::array<std::string, 7> mandatoryHistoFields{
         "title", "variable", "d1", "d2", "nbins", "xmin", "xmax" };
     const std::array<std::string, 6> extra2DHistoFields{ "variableY", "d1Y",
                                                          "d2Y",       "nbinsY",
                                                          "ymin",      "ymax" };
 
     bool allMandatoryFields{ true };
 
     for ( const auto& field : mandatoryFields ) {
         if ( !m_config.contains( field ) ) {
             std::cerr << "ERROR : json does not contain required field: " << field
                       << std::endl;
             allMandatoryFields = false;
             continue;
         }
         if ( field == "histograms" ) {
             const json& jHistos{ m_config.at( "histograms" ) };
             for ( const auto& hInfo : jHistos ) {
                 for ( const auto& hField : mandatoryHistoFields ) {
                     if ( !hInfo.contains( hField ) ) {
                         std::cerr
                             << "ERROR : json does not contain required field for histogram definition: "
                             << hField << std::endl;
                         allMandatoryFields = false;
                     }
                 }
                 if ( hInfo.contains( extra2DHistoFields[0] ) ) {
                     for ( const auto& hField : extra2DHistoFields ) {
                         if ( !hInfo.contains( hField ) ) {
                             std::cerr
                                 << "ERROR : json does not contain required field for 2D histogram definition: "
                                 << hField << std::endl;
                             allMandatoryFields = false;
                         }
                     }
                 }
             }
         }
     }
 
     return allMandatoryFields;
 }
 
 bool TestDecayModel::run()
 {
     // Check that we have, and then get all the mandatory fields first
     if ( !checkMandatoryFields() ) {
         std::cerr << "ERROR : json does not contain all mandatory fields - dumping config for debugging:\n";
         std::cerr << m_config << std::endl;
         return false;
     }
 
     const auto parentName{ m_config.at( "parent" ).get<std::string>() };
     const auto daughterNames{
         m_config.at( "daughters" ).get<std::vector<std::string>>() };
     const auto modelNames{
         m_config.at( "models" ).get<std::vector<std::string>>() };
     const auto modelParameters{
         m_config.at( "parameters" ).get<std::vector<std::vector<std::string>>>() };
     const auto outFileName{ m_config.at( "outfile" ).get<std::string>() };
     const auto nEvents{ m_config.at( "events" ).get<int>() };
     const auto refFileName{ m_config.at( "reference" ).get<std::string>() };
 
     // Then check for optional fields, setting default values if not present
 
     const auto grandDaughterNames{
         ( m_config.contains( "grand_daughters" ) &&
           m_config.at( "grand_daughters" ).is_array() )
             ? m_config.at( "grand_daughters" )
                   .get<std::vector<std::vector<std::string>>>()
             : std::vector<std::vector<std::string>>{} };
 
     const auto extraCommands{
         ( m_config.contains( "extras" ) && m_config.at( "extras" ).is_array() )
             ? m_config.at( "extras" ).get<std::vector<std::string>>()
             : std::vector<std::string>{} };
 
     const auto debugFlag{ ( m_config.contains( "debug_flag" ) &&
                             m_config.at( "debug_flag" ).is_boolean() )
                               ? m_config.at( "debug_flag" ).get<bool>()
                               : false };
 
     std::vector<bool> doConjDecay;
     if ( m_config.contains( "do_conjugate_decay" ) &&
          m_config.at( "do_conjugate_decay" ).is_array() ) {
         doConjDecay = m_config.at( "do_conjugate_decay" ).get<std::vector<bool>>();
     }
     if ( doConjDecay.size() != modelNames.size() ) {
         doConjDecay.resize( modelNames.size(), false );
     }
 
     // Initialise the EvtGen object and hence the EvtPDL tables
     // The EvtGen object is used by generateEvents, while the
     // latter are also used within createDecFile
 
     // Define the random number generator
     auto randomEngine{ std::make_unique<EvtMTRandomEngine>() };
 
     EvtAbsRadCorr* radCorrEngine{ nullptr };
     std::list<EvtDecayBase*> extraModels;
 
 #ifdef EVTGEN_EXTERNAL
     bool convertPythiaCodes( false );
     bool useEvtGenRandom( true );
     EvtExternalGenList genList( convertPythiaCodes, "", "gamma", useEvtGenRandom );
     radCorrEngine = genList.getPhotosModel();
     extraModels = genList.getListOfModels();
 #endif
 
     EvtGen theGen( "../DECAY.DEC", "../evt.pdl", randomEngine.get(),
                    radCorrEngine, &extraModels );
 
     /*! Creates a decay file based on json file input. */
     const std::string decFileName{
         outFileName.substr( 0, outFileName.size() - 5 ) };
     const std::string decFile{ createDecFile( parentName, daughterNames,
                                               grandDaughterNames, modelNames,
                                               modelParameters, doConjDecay,
                                               extraCommands, decFileName ) };
 
     /*! Define the root output file and histograms to be saved. */
     std::unique_ptr<TFile> outFile{
         TFile::Open( outFileName.c_str(), "recreate" ) };
     defineHistos( outFile.get() );
 
     /*!  Generate events and fill histograms. */
     generateEvents( theGen, decFile, parentName, doConjDecay[0], nEvents,
                     debugFlag );
 
     // Normalize histograms.
     for ( auto& [_, hist] : m_1DhistVect ) {
         const double area{ hist->Integral() };
         if ( area > 0.0 ) {
             hist->Scale( 1.0 / area );
         }
     }
     for ( auto& [_, hist] : m_2DhistVect ) {
         const double area{ hist->Integral() };
         if ( area > 0.0 ) {
             hist->Scale( 1.0 / area );
         }
     }
 
     /*! Compare with reference histograms. */
     compareHistos( refFileName );
 
     // Write output.
     outFile->cd();
     for ( auto& [_, hist] : m_1DhistVect ) {
         hist->Write();
     }
     for ( auto& [_, hist] : m_2DhistVect ) {
         hist->Write();
     }
     if ( m_mixedHist ) {
         m_mixedHist->Write();
     }
     outFile->Close();
     std::cout << "Created output file: " << outFileName.c_str() << std::endl;
 
     return true;
 }
 
 std::string TestDecayModel::createDecFile(
     const std::string& parent, const std::vector<std::string>& daughterNames,
     const std::vector<std::vector<std::string>>& grandDaughterNames,
     const std::vector<std::string>& modelNames,
     const std::vector<std::vector<std::string>>& parameters,
     const std::vector<bool>& doConjDecay,
     const std::vector<std::string>& extras, const std::string decFileName ) const
 {
     // Create (or overwrite) the decay file
     std::string decName( decFileName + ".dec" );
     std::ofstream decFile;
     decFile.open( decName.c_str() );
 
     // Create daughter aliases if needed
     std::vector<std::string> aliasPrefix;
     for ( long unsigned int daughter_index{ 0 };
           daughter_index < daughterNames.size(); daughter_index++ ) {
         if ( !grandDaughterNames.empty() &&
              !grandDaughterNames[daughter_index].empty() ) {
             decFile << "Alias My" << daughterNames[daughter_index] << " "
                     << daughterNames[daughter_index] << std::endl;
             if ( doConjDecay[daughter_index + 1] ) {
                 const EvtId daugID{
                     EvtPDL::getId( daughterNames[daughter_index] ) };
                 const EvtId daugConjID{ EvtPDL::chargeConj( daugID ) };
                 const std::string conjName{ daugConjID.getName() };
                 std::string conjName_Alias{ daugConjID.getName() };
                 if ( std::find( std::begin( daughterNames ),
                                 std::end( daughterNames ), daugConjID.getName() ) !=
                      std::end( daughterNames ) ) {
                     conjName_Alias = conjName_Alias + "_" + daughter_index;
                 }
                 decFile << "Alias My" << conjName_Alias << " " << conjName
                         << std::endl;
                 decFile << "ChargeConj My" << daughterNames[daughter_index]
                         << "  My" << conjName_Alias << std::endl;
             } else if ( doConjDecay[0] ) {
                 decFile << "ChargeConj My" << daughterNames[daughter_index]
                         << "  My" << daughterNames[daughter_index] << std::endl;
             }
             aliasPrefix.push_back( "My" );
         } else {
             aliasPrefix.push_back( "" );
         }
     }
 
     for ( const auto& iExtra : extras ) {
         decFile << iExtra << std::endl;
     }
 
     // Parent decay
     decFile << "Decay " << parent << std::endl;
     decFile << "1.0";
 
     for ( long unsigned int daughter_index{ 0 };
           daughter_index < daughterNames.size(); daughter_index++ ) {
         decFile << " " << aliasPrefix[daughter_index]
                 << daughterNames[daughter_index];
     }
 
     decFile << " " << modelNames[0];
 
     for ( const auto& par : parameters[0] ) {
         decFile << " " << par;
     }
 
     decFile << ";" << std::endl;
     decFile << "Enddecay" << std::endl;
     if ( doConjDecay[0] ) {
         EvtId parID{ EvtPDL::getId( parent ) };
         EvtId parConjID{ EvtPDL::chargeConj( parID ) };
         decFile << "CDecay " << parConjID.getName() << std::endl;
     }
 
     // Daughter decays into granddaughters
     for ( long unsigned int daughter_index{ 0 };
           daughter_index < grandDaughterNames.size(); daughter_index++ ) {
         if ( grandDaughterNames[daughter_index].empty() )
             continue;
         decFile << "Decay " << aliasPrefix[daughter_index]
                 << daughterNames[daughter_index] << std::endl;
         decFile << "1.0";
         for ( long unsigned int grandDaughter_index{ 0 };
               grandDaughter_index < grandDaughterNames[daughter_index].size();
               grandDaughter_index++ ) {
             decFile << " "
                     << grandDaughterNames[daughter_index][grandDaughter_index];
         }
         decFile << " " << modelNames[daughter_index + 1];
         for ( const auto& par : parameters[daughter_index + 1] ) {
             decFile << " " << par;
         }
         decFile << ";" << std::endl;
         decFile << "Enddecay" << std::endl;
         if ( doConjDecay[daughter_index + 1] ) {
             EvtId daugID{ EvtPDL::getId( daughterNames[daughter_index] ) };
             EvtId daugConjID{ EvtPDL::chargeConj( daugID ) };
             std::string conjName_Alias{ daugConjID.getName() };
             if ( std::find( std::begin( daughterNames ),
                             std::end( daughterNames ), daugConjID.getName() ) !=
                  std::end( daughterNames ) ) {
                 conjName_Alias = conjName_Alias + "_" + daughter_index;
             }
             decFile << "CDecay " << aliasPrefix[daughter_index]
                     << conjName_Alias << std::endl;
         }
     }
 
     decFile << "End" << std::endl;
 
     decFile.close();
 
     return decName;
 }
 
 void TestDecayModel::defineHistos( TFile* outFile )
 {
     // Histogram information
     const json& jHistos{ m_config.at( "histograms" ) };
     const size_t nHistos{ jHistos.size() };
 
     m_1DhistVect.reserve( nHistos );
     m_2DhistVect.reserve( nHistos );
 
     for ( const auto& hInfo : jHistos ) {
         const auto varTitle{ hInfo.at( "title" ).get<std::string>() };
 
         const auto varName{ hInfo.at( "variable" ).get<std::string>() };
         // Integer values that define what particles need to be used
         // for invariant mass combinations or helicity angles etc
         const auto d1{ hInfo.at( "d1" ).get<int>() };
         const auto d2{ hInfo.at( "d2" ).get<int>() };
 
         const auto nBins{ hInfo.at( "nbins" ).get<int>() };
         const auto xmin{ hInfo.at( "xmin" ).get<double>() };
         const auto xmax{ hInfo.at( "xmax" ).get<double>() };
 
         std::string histName( varName.c_str() );
         if ( d1 != 0 ) {
             histName += "_";
             histName += std::to_string( d1 );
         }
         if ( d2 != 0 ) {
             histName += "_";
             histName += std::to_string( d2 );
         }
 
         if ( !hInfo.contains( "variableY" ) ) {
             TH1* hist{ new TH1D{ histName.c_str(), varTitle.c_str(), nBins,
                                  xmin, xmax } };
             hist->SetDirectory( outFile );
             m_1DhistVect.emplace_back(
                 std::make_pair( TestInfo( varName, d1, d2 ), hist ) );
             continue;
         } else {
             const auto varNameY{ hInfo.at( "variableY" ).get<std::string>() };
             const auto d1Y{ hInfo.at( "d1Y" ).get<int>() };
             const auto d2Y{ hInfo.at( "d2Y" ).get<int>() };
 
             const auto nBinsY{ hInfo.at( "nbinsY" ).get<int>() };
             const auto ymin{ hInfo.at( "ymin" ).get<double>() };
             const auto ymax{ hInfo.at( "ymax" ).get<double>() };
 
             histName += "_";
             histName += varNameY;
             if ( d1Y != 0 ) {
                 histName += "_";
                 histName += std::to_string( d1Y );
             }
             if ( d2Y != 0 ) {
                 histName += "_";
                 histName += std::to_string( d2Y );
             }
             TH2* hist{ new TH2D{ histName.c_str(), varTitle.c_str(), nBins,
                                  xmin, xmax, nBinsY, ymin, ymax } };
             hist->SetDirectory( outFile );
             m_2DhistVect.emplace_back( std::make_pair(
                 TestInfo( varName, d1, d2, varNameY, d1Y, d2Y ), hist ) );
         }
     }
 
     // For the case where the parent is either a neutral B or D add a mixed/unmixed histogram
     const auto parentName{ m_config.at( "parent" ).get<std::string>() };
     const int parentID{ abs( EvtPDL::getStdHep( EvtPDL::getId( parentName ) ) ) };
     if ( parentID == 511 || parentID == 531 || parentID == 421 ) {
         const std::string varTitle{ parentName + " mixed" };
         m_mixedHist = new TH1D{ "mixed", varTitle.c_str(), 2, 0.0, 2.0 };
         // TODO maybe set bin labels?
         m_mixedHist->SetDirectory( outFile );
     }
 }
 
 void TestDecayModel::generateEvents( EvtGen& theGen, const std::string& decFile,
                                      const std::string& parentName,
                                      const bool doConjDecay, const int nEvents,
                                      const bool debug_flag )
 {
     // Read the decay file
     theGen.readUDecay( decFile.c_str() );
 
     // Generate the decays
     EvtId parId{ EvtPDL::getId( parentName.c_str() ) };
     EvtId conjId{ doConjDecay ? EvtPDL::chargeConj( parId ) : parId };
     for ( int i{ 0 }; i < nEvents; i++ ) {
         if ( i % 1000 == 0 ) {
             std::cout << "Event " << nEvents - i << std::endl;
         }
 
         // Initial 4-momentum and particle
         EvtVector4R pInit( EvtPDL::getMass( parId ), 0.0, 0.0, 0.0 );
         EvtParticle* parent{ nullptr };
         if ( EvtRandom::Flat() < 0.5 ) {
             parent = EvtParticleFactory::particleFactory( parId, pInit );
         } else {
             parent = EvtParticleFactory::particleFactory( conjId, pInit );
         }
 
         theGen.generateDecay( parent );
 
         // Check for mixing (and fill histogram)
         EvtParticle* prodParent{ nullptr };
         if ( m_mixedHist ) {
             if ( parent->getNDaug() == 1 ) {
                 prodParent = parent;
                 parent = prodParent->getDaug( 0 );
 
                 m_mixedHist->Fill( 1 );
             } else {
                 m_mixedHist->Fill( 0 );
             }
         }
 
         // To debug
 
         if ( debug_flag ) {
             std::cout << "Parent PDG code: " << parent->getPDGId()
                       << " has daughters " << parent->getNDaug() << std::endl;
             for ( size_t iDaughter{ 0 }; iDaughter < parent->getNDaug();
                   iDaughter++ ) {
                 std::cout << "Parent PDG code of daughter " << iDaughter
                           << " : " << parent->getDaug( iDaughter )->getPDGId()
                           << " has daughters "
                           << parent->getDaug( iDaughter )->getNDaug()
                           << std::endl;
 
                 for ( size_t iGrandDaughter{ 0 };
                       iGrandDaughter < parent->getDaug( iDaughter )->getNDaug();
                       iGrandDaughter++ ) {
                     std::cout << "Parent PDG code of grand daughter "
                               << iGrandDaughter << " : "
                               << parent->getDaug( iDaughter )
                                      ->getDaug( iGrandDaughter )
                                      ->getPDGId()
                               << " has daughters "
                               << parent->getDaug( iDaughter )
                                      ->getDaug( iGrandDaughter )
                                      ->getNDaug()
                               << std::endl;
                 }
             }
         }
 
+        int leadingChargedDaughter = -1;
+        const std::string fsrStr = "_FSRPhotons";
+        const std::string ewStr = "_EnergyWeight";
+
         // Store information
         for ( auto& [info, hist] : m_1DhistVect ) {
-            const double value{ getValue( parent, info.getName(), info.getd1(),
-                                          info.getd2() ) };
+            if ( !hist ) {
+                continue;
+            }
+
+            const std::string varName = info.getName();
+            const std::string::size_type findFSRstr = varName.find( fsrStr );
+
+            // If the variable name does not have the substring '_FSRPhotons', then add an entry per event and continue
+            if ( findFSRstr == std::string::npos ) {
+                const double value{
+                    getValue( parent, varName, info.getd1(), info.getd2() ) };
 
-            if ( hist ) {
                 hist->Fill( value );
+                continue;
+            }
+
+            // If otherwise the variable contains the substring '_FSRPhotons', then add an entry per photon
+
+            if ( leadingChargedDaughter == -1 ) {
+                leadingChargedDaughter = findChargedDaugtherWithMaxE( parent );
+            }
+
+            std::string reducedVarName = varName;
+            reducedVarName.erase( findFSRstr, fsrStr.length() );
+
+            const std::string::size_type findEwStr = reducedVarName.find( ewStr );
+
+            // If variable name has substring '_EnergyWeight' then add an entry weighted by the photon energy
+            bool energyWeight = false;
+            if ( findEwStr != std::string::npos ) {
+                reducedVarName.erase( findEwStr, ewStr.length() );
+                energyWeight = true;
+            }
+
+            for ( int iDaug{ 0 }; iDaug < (int)parent->getNDaug(); iDaug++ ) {
+                const EvtParticle* iDaughter = parent->getDaug( iDaug );
+
+                if ( iDaughter->getAttribute( "FSR" ) == 1 ) {
+                    const double value{ getValue( parent, reducedVarName,
+                                                  iDaug + 1,
+                                                  leadingChargedDaughter + 1 ) };
+                    if ( energyWeight ) {
+                        const double photonEnergy = iDaughter->getP4().get( 0 );
+                        hist->Fill( value, photonEnergy );
+                    } else {
+                        hist->Fill( value );
+                    }
+                }
             }
         }
+
         for ( auto& [info, hist] : m_2DhistVect ) {
+            if ( !hist ) {
+                continue;
+            }
+
             const double valueX{ getValue( parent, info.getName(), info.getd1(),
                                            info.getd2() ) };
             const double valueY{ getValue( parent, info.getName( 2 ),
                                            info.getd1( 2 ), info.getd2( 2 ) ) };
-            if ( hist ) {
-                hist->Fill( valueX, valueY );
-            }
+            hist->Fill( valueX, valueY );
         }
 
         if ( debug_flag ) {
             if ( prodParent ) {
                 prodParent->printTree();
             } else {
                 parent->printTree();
             }
         }
 
         // Cleanup
         if ( prodParent ) {
             prodParent->deleteTree();
         } else {
             parent->deleteTree();
         }
     }
 }
 
+int TestDecayModel::findChargedDaugtherWithMaxE( const EvtParticle* parent ) const
+{
+    /* This function returns the index of the charged daughter with the highest energy 
+     * following the sign convention below. */
+
+    double max_E = 0;
+    double max_index = 0;
+
+    const int parentCh3 = EvtPDL::chg3( parent->getId() );
+
+    for ( int iDaug{ 0 }; iDaug < (int)parent->getNDaug(); iDaug++ ) {
+        const EvtParticle* iDaughter = parent->getDaug( iDaug );
+        const int dauCh3 = EvtPDL::chg3( iDaughter->getId() );
+        const double daugE = iDaughter->getP4LabBeforeFSR().get( 0 );
+        // Sign convention: take negative daughter if mother is neutral,
+        // otherwise the one with the same sign as the mother.
+        if ( ( ( parentCh3 == 0 && dauCh3 < 0 ) || ( parentCh3 * dauCh3 > 0 ) ) &&
+             daugE > max_E ) {
+            max_E = daugE;
+            max_index = iDaug;
+        }
+    }
+    return max_index;
+}
+
 double TestDecayModel::getValue( const EvtParticle* parent,
                                  const std::string& varName, const int d1,
                                  const int d2 ) const
 {
     double value{ std::numeric_limits<double>::quiet_NaN() };
     if ( !parent ) {
         return value;
     }
 
     const int NDaugMax( parent->getNDaug() );
     // If variable name contains "_daugX", we are interested in daughters of daughter X
     // Else we are interested in daughters
     const EvtParticle* selectedParent{ parent };
     std::string selectedVarName{ varName };
     if ( varName.find( "_daug" ) != std::string::npos ) {
         // Get daughter index from last character in string
         const int iDaughter{ varName.back() - '0' };
         selectedVarName = varName.substr( 0, varName.size() - 6 );
         if ( iDaughter > 0 && iDaughter <= NDaugMax ) {
             selectedParent = parent->getDaug( iDaughter - 1 );
         } else {
             return value;
         }
     }
 
     const int sel_NDaugMax( selectedParent->getNDaug() );
     const EvtParticle* par1{ d1 > 0 && d1 <= sel_NDaugMax
                                  ? selectedParent->getDaug( d1 - 1 )
                                  : nullptr };
     const EvtParticle* par2{ d2 > 0 && d2 <= sel_NDaugMax
                                  ? selectedParent->getDaug( d2 - 1 )
                                  : nullptr };
     const EvtParticle* par3{ sel_NDaugMax > 0
                                  ? selectedParent->getDaug( sel_NDaugMax - 1 )
                                  : nullptr };
 
     // 4-momenta in parent rest frame
     const EvtVector4R p1{ par1 != nullptr ? par1->getP4() : EvtVector4R() };
     const EvtVector4R p2{ par2 != nullptr ? par2->getP4() : EvtVector4R() };
     const EvtVector4R p3{ par3 != nullptr ? par3->getP4() : EvtVector4R() };
 
     // 4-momenta in lab frame (1st parent in decay tree)
     const EvtVector4R p1_lab{ par1 != nullptr ? par1->getP4Lab() : EvtVector4R() };
     const EvtVector4R p2_lab{ par2 != nullptr ? par2->getP4Lab() : EvtVector4R() };
     const EvtVector4R p3_lab{ par3 != nullptr ? par3->getP4Lab() : EvtVector4R() };
 
     if ( !selectedVarName.compare( "id" ) ) {
         // StdHep ID of one of the daughters (controlled by d1) or the parent
         if ( par1 ) {
             value = par1->getPDGId();
         } else {
             value = selectedParent->getPDGId();
         }
 
     } else if ( !selectedVarName.compare( "parMass" ) ) {
         // Parent invariant mass
         value = selectedParent->mass();
 
     } else if ( !selectedVarName.compare( "mass" ) ) {
         // Invariant mass
         if ( d2 != 0 ) {
             // Invariant 4-mass combination of particles d1 and d2
             value = ( p1 + p2 ).mass();
         } else {
             // Invariant mass of particle d1 only
             value = p1.mass();
         }
 
     } else if ( !selectedVarName.compare( "massSq" ) ) {
         // Invariant mass
         if ( d2 != 0 ) {
             // Invariant 4-mass combination of particles d1 and d2
             value = ( p1 + p2 ).mass2();
         } else {
             // Invariant mass of particle d1 only
             value = p1.mass2();
         }
     } else if ( !selectedVarName.compare( "mPrime" ) ) {
         // pick up first unused particle rather than last one
         if ( sel_NDaugMax != 3 ) {
             return -1;
         }
         int unused{ 0 };
         for ( int ii{ 1 }; ii <= sel_NDaugMax; ++ii ) {
             if ( ii != d1 && ii != d2 ) {
                 unused = ii;
                 break;
             }
         }
         if ( unused == 0 ) {
             unused = sel_NDaugMax;
         }
         const auto parL{ selectedParent->getDaug( unused - 1 ) };
         //            const auto& pL = parL->getP4();
 
         const double mB{ selectedParent->mass() };
         const double m1{ par1->mass() };
         const double m2{ par2->mass() };
         const double m3{ parL->mass() };
         const double m12{ ( p1 + p2 ).mass() };
         const double m12norm{
             2 * ( ( m12 - ( m1 + m2 ) ) / ( mB - ( m1 + m2 + m3 ) ) ) - 1 };
         value = acos( m12norm ) / EvtConst::pi;
 
     } else if ( !selectedVarName.compare( "thetaPrime" ) ) {
         // pick up first unused particle rather than last one
         if ( sel_NDaugMax != 3 ) {
             return -1;
         }
         int unused{ 0 };
         for ( int ii{ 1 }; ii <= sel_NDaugMax; ++ii ) {
             if ( ii != d1 && ii != d2 ) {
                 unused = ii;
                 break;
             }
         }
         if ( unused == 0 ) {
             unused = sel_NDaugMax;
         }
         const auto parL{ selectedParent->getDaug( unused - 1 ) };
         const auto& pL{ parL->getP4() };
 
         const double mB{ selectedParent->mass() };
         const double m1{ p1.mass() };
         const double m2{ p2.mass() };
         const double m3{ pL.mass() };
         double mBSq{ mB * mB };
         double m1Sq{ m1 * m1 };
         double m2Sq{ m2 * m2 };
         double m3Sq{ m3 * m3 };
         const double m12{ ( p1 + p2 ).mass() };
         const double m13{ ( p1 + p3 ).mass() };
         const double m12Sq{ m12 * m12 };
         const double m13Sq{ m13 * m13 };
         double en1{ ( m12Sq - m2Sq + m1Sq ) / ( 2.0 * m12 ) };
         double en3{ ( mBSq - m12Sq - m3Sq ) / ( 2.0 * m12 ) };
         double p1_12{ std::sqrt( en1 * en1 - m1Sq ) };
         double p3_12{ std::sqrt( en3 * en3 - m3Sq ) };
         double cosTheta{ ( -m13Sq + m1Sq + m3Sq + 2. * en1 * en3 ) /
                          ( 2. * p1_12 * p3_12 ) };
         value = acos( cosTheta ) / EvtConst::pi;
 
     } else if ( !selectedVarName.compare( "pSumSq" ) ) {
         // Invariant momentum sum squared
         value = ( p1 + p2 ).mass2();
 
     } else if ( !selectedVarName.compare( "pDiffSq" ) ) {
         // Invariant momentum difference squared
         value = ( p1 - p2 ).mass2();
 
     } else if ( !selectedVarName.compare( "mass3" ) ) {
         // Invariant mass of 3 daughters
         value = ( p1 + p2 + p3 ).mass();
 
     } else if ( !selectedVarName.compare( "mass3_specified" ) ) {
         // Invariant mass of 3 daughters, d1 is first daughter
         // second daughter is d1 + 1, d2 is last daughter
         const EvtParticle* daug2{ selectedParent->getDaug( d1 ) };
         const EvtParticle* daug3{ selectedParent->getDaug( d2 - 1 ) };
 
         const EvtVector4R p2_specified{ daug2 != nullptr ? daug2->getP4Lab()
                                                          : EvtVector4R() };
         const EvtVector4R p3_specified{ daug3 != nullptr ? daug3->getP4Lab()
                                                          : EvtVector4R() };
 
         value = ( p1 + p2_specified + p3_specified ).mass();
 
     } else if ( !selectedVarName.compare( "cosTheta3" ) ) {
         // Cosine of the polar angle of the momentum of d1 + d2 + d3
         const EvtVector4R p123{ p1 + p2 + p3 };
         if ( p123.d3mag() > 0.0 ) {
             value = p123.get( 3 ) / p123.d3mag();
         }
 
     } else if ( !selectedVarName.compare( "pLab" ) ) {
         // Momentum of particle d1 in lab frame
         value = p1_lab.d3mag();
 
     } else if ( !selectedVarName.compare( "p" ) ) {
         // Momentum of particle d1 in parent rest frame
         value = p1.d3mag();
 
     } else if ( !selectedVarName.compare( "pLabSq" ) ) {
         // Momentum squared of particle d1 (in lab frame)
         const double p1_lab_x{ p1_lab.get( 1 ) };
         const double p1_lab_y{ p1_lab.get( 2 ) };
         const double p1_lab_z{ p1_lab.get( 3 ) };
         value = p1_lab_x * p1_lab_x + p1_lab_y * p1_lab_y + p1_lab_z * p1_lab_z;
 
     } else if ( !selectedVarName.compare( "pSq" ) ) {
         // Momentum squared of particle d1 (in lab frame)
         const double p1_x{ p1.get( 1 ) };
         const double p1_y{ p1.get( 2 ) };
         const double p1_z{ p1.get( 3 ) };
         value = p1_x * p1_x + p1_y * p1_y + p1_z * p1_z;
 
     } else if ( !selectedVarName.compare( "pxLab" ) ) {
         // x momentum of particle d1 in lab frame
         value = p1_lab.get( 1 );
 
     } else if ( !selectedVarName.compare( "px" ) ) {
         // x momentum of particle d1 in parent frame
         value = p1.get( 1 );
 
     } else if ( !selectedVarName.compare( "pyLab" ) ) {
         // y momentum of particle d1 in lab frame
         value = p1_lab.get( 2 );
 
     } else if ( !selectedVarName.compare( "py" ) ) {
         // y momentum of particle d1 in parent frame
         value = p1.get( 2 );
 
     } else if ( !selectedVarName.compare( "pzLab" ) ) {
         // z momentum of particle d1 in lab frame
         value = p1_lab.get( 3 );
 
     } else if ( !selectedVarName.compare( "pz" ) ) {
         // z momentum of particle d1 in parent frame
         value = p1.get( 3 );
 
     } else if ( !selectedVarName.compare( "cosHel" ) ||
                 !selectedVarName.compare( "absCosHel" ) ) {
         // Cosine of helicity angle
         EvtVector4R p12;
         EvtVector4R p1Res;
 
         if ( d2 != 0 ) {
             // Resonance center-of-mass system (d1 and d2)
             p12 = p1_lab + p2_lab;
             // Boost vector
             const EvtVector4R boost{ p12.get( 0 ), -p12.get( 1 ), -p12.get( 2 ),
                                      -p12.get( 3 ) };
             // Momentum of particle d1 in resonance frame, p1Res
             p1Res = boostTo( p1_lab, boost );
         } else {
             // The resonance is d1
             p12 = p1;
 
             // Find its first daughter
             const EvtParticle* gpar{ par1 != nullptr ? par1->getDaug( 0 )
                                                      : nullptr };
 
             p1Res = gpar != nullptr ? gpar->getP4() : EvtVector4R();
         }
 
         // Cosine of angle between p1Res and momentum of resonance in parent frame
         const double p1ResMag{ p1Res.d3mag() };
         const double p12Mag{ p12.d3mag() };
         if ( p1ResMag > 0.0 && p12Mag > 0.0 ) {
             value = -p1Res.dot( p12 ) / ( p1ResMag * p12Mag );
         }
 
         if ( !selectedVarName.compare( "absCosHel" ) ) {
             value = std::abs( value );
         }
 
     } else if ( !selectedVarName.compare( "cosHelTau" ) ) {
         // Works only for B -> X nu_1 tau -> pi nu_2.
         // Cosine of helicity angle between pi momentum and opposite W momentum -(nu_1 + tau)
         // in the tau rest frame. Index d1 must match with tau.
         // p3 (momentum of last daughter) is taken as the neutrino momentum
 
         // W momentum
         const EvtVector4R p_W{ p1_lab + p3_lab };
 
         // Index d2 must match the index of the pion (daughter of tau)
         const EvtParticle* pion{
             selectedParent->getDaug( d1 - 1 )->getDaug( d2 - 1 ) };
         const EvtVector4R p_pion{ pion != nullptr ? pion->getP4Lab()
                                                   : EvtVector4R() };
 
         // Boost vector to tau frame
         const EvtVector4R boost{ p1_lab.get( 0 ), -p1_lab.get( 1 ),
                                  -p1_lab.get( 2 ), -p1_lab.get( 3 ) };
 
         // Boost both momenta to tau frame
         const EvtVector4R p_W_boosted{ boostTo( p_W, boost ) };
         const EvtVector4R p_pion_boosted{ boostTo( p_pion, boost ) };
         // Cosine of angle between opposite W momentum and pion momentum in tau frame
         const double p_W_boostedMag{ p_W_boosted.d3mag() };
         const double p_pion_boostedMag{ p_pion_boosted.d3mag() };
 
         if ( p_W_boostedMag > 0.0 && p_pion_boostedMag > 0.0 ) {
             value = -p_W_boosted.dot( p_pion_boosted ) /
                     ( p_W_boostedMag * p_pion_boostedMag );
         }
 
     } else if ( !selectedVarName.compare( "cosHelDiTau" ) ||
                 !selectedVarName.compare( "cosHelDiTau_over05" ) ||
                 !selectedVarName.compare( "cosHelDiTau_under05" ) ) {
         // Works for B -> tau (pi nu) tau -> (pi nu), B -> phi (-> KK) l l decays,
         // B -> 4 pions, or similar decays..
         // Cosine of helicity angle between pi momentum and opposite B momentum in tau rest frame.
         // Index d1 must match with tau
 
         if ( sel_NDaugMax < 2 || sel_NDaugMax > 4 ) {
             return value;
         }
 
         // B momentum
         const EvtVector4R p_B{ selectedParent->getP4Lab() };
 
         // Index d2 must match the index of the pion (daughter of tau)
         const EvtParticle* pion_1{
             sel_NDaugMax <= 3
                 ? selectedParent->getDaug( d1 - 1 )->getDaug( d2 - 1 )
                 : selectedParent->getDaug( d1 - 1 ) };
         const EvtVector4R p_pion_1{ pion_1 != nullptr ? pion_1->getP4Lab()
                                                       : EvtVector4R() };
 
         const EvtVector4R p_first_daughter{
             sel_NDaugMax <= 3 ? selectedParent->getDaug( d1 - 1 )->getP4Lab()
                               : selectedParent->getDaug( d1 - 1 )->getP4Lab() +
                                     selectedParent->getDaug( d1 )->getP4Lab() };
 
         // Boost vector to tau frame
         const EvtVector4R boost_1{ p_first_daughter.get( 0 ),
                                    -p_first_daughter.get( 1 ),
                                    -p_first_daughter.get( 2 ),
                                    -p_first_daughter.get( 3 ) };
 
         // Boost both momenta to tau frame
         const EvtVector4R p_B_boosted_1{ boostTo( p_B, boost_1 ) };
         const EvtVector4R p_pion_boosted_1{ boostTo( p_pion_1, boost_1 ) };
         // Cosine of angle between opposite W momentum and pion momentum in tau frame
         const double p_B_boosted_1_Mag{ p_B_boosted_1.d3mag() };
         const double p_pion_boosted_1_Mag{ p_pion_boosted_1.d3mag() };
 
         double hel1{ std::numeric_limits<double>::quiet_NaN() };
         double hel{ std::numeric_limits<double>::quiet_NaN() };
 
         if ( p_B_boosted_1_Mag > 0.0 && p_pion_boosted_1_Mag > 0.0 ) {
             hel1 = -p_B_boosted_1.dot( p_pion_boosted_1 ) /
                    ( p_B_boosted_1_Mag * p_pion_boosted_1_Mag );
         }
 
         if ( ( !selectedVarName.compare( "cosHelDiTau_over05" ) ) ||
              ( !selectedVarName.compare( "cosHelDiTau_under05" ) ) ) {
             // Works for B -> tau (pi nu) tau -> (pi nu) or similar decays.
             // Cosine of helicity angle between pi momentum and opposite B momentum in tau rest frame
             // Index d1 must match with tau; cosHelicity above +0.5 or below -0.5
 
             // Index d2 must match the index of the pion (daughter of tau)
             const EvtParticle* pion{
                 sel_NDaugMax <= 3 ? selectedParent->getDaug( 0 )->getDaug( d2 - 1 )
                                   : selectedParent->getDaug( 0 ) };
             const EvtVector4R p_pion{ pion != nullptr ? pion->getP4Lab()
                                                       : EvtVector4R() };
             // Boost vector to tau frame
             const EvtVector4R p_second_daughter{
                 sel_NDaugMax == 2
                     ? selectedParent->getDaug( 0 )->getP4Lab()
                     : selectedParent->getDaug( 0 )->getP4Lab() +
                           selectedParent->getDaug( 1 )->getP4Lab() };
 
             const EvtVector4R boost{ p_second_daughter.get( 0 ),
                                      -p_second_daughter.get( 1 ),
                                      -p_second_daughter.get( 2 ),
                                      -p_second_daughter.get( 3 ) };
 
             // Boost both momenta to tau frame
             const EvtVector4R p_B_boosted{ boostTo( p_B, boost ) };
             const EvtVector4R p_pion_boosted{ boostTo( p_pion, boost ) };
             // Cosine of angle between opposite W momentum and pion momentum in tau frame
             const double p_B_boostedMag{ p_B_boosted.d3mag() };
             const double p_pion_boostedMag{ p_pion_boosted.d3mag() };
 
             if ( p_B_boostedMag > 0.0 && p_pion_boostedMag > 0.0 ) {
                 hel = -p_B_boosted.dot( p_pion_boosted ) /
                       ( p_B_boostedMag * p_pion_boostedMag );
             }
         }
 
         if ( !selectedVarName.compare( "cosHelDiTau" ) ||
              ( hel > 0.5 && !selectedVarName.compare( "cosHelDiTau_over05" ) ) ||
              ( hel < -0.5 && !selectedVarName.compare( "cosHelDiTau_under05" ) ) ) {
             value = hel1;
         }
 
     } else if ( !selectedVarName.compare( "cosAcoplanarityAngle" ) ||
                 !selectedVarName.compare( "acoplanarityAngle" ) ) {
         // Acoplanarity angle or cosine for B -> tau (pi nu) tau -> (pi nu),
         // B -> phi (-> KK) l l decays, B -> 4 pions, or similar decays
 
         value = getCosAcoplanarityAngle( selectedParent, sel_NDaugMax, d1, d2 );
         if ( !selectedVarName.compare( "acoplanarityAngle" ) ) {
             value = std::acos( value ) * 180.0 / EvtConst::pi;    // degrees
         }
 
     } else if ( !selectedVarName.compare( "cosTheta" ) ) {
         // Cosine of polar angle of first daughter in lab frame
         const double p1_lab_mag{ p1_lab.d3mag() };
         if ( p1_lab_mag > 0.0 ) {
             value = p1_lab.get( 3 ) / p1_lab_mag;
         }
 
     } else if ( !selectedVarName.compare( "phi" ) ) {
         // Azimuthal angle of first daughter in lab frame (degrees)
         const double p1_lab_mag{ p1_lab.d3mag() };
         if ( p1_lab_mag > 0.0 ) {
             value = atan2( p1_lab.get( 1 ), p1_lab.get( 2 ) ) * 180.0 /
                     EvtConst::pi;
         }
 
+    } else if ( !selectedVarName.compare( "openingAngle" ) ) {
+        // Polar angle between first and second daughters in parent's frame
+
+        const double cost{ p1.dot( p2 ) / ( p1.d3mag() * p2.d3mag() ) };
+
+        value = acos( cost ) * 180.0 / EvtConst::pi;
+
     } else if ( !selectedVarName.compare( "decayangle" ) ) {
         // Polar angle between first and second daughters in lab frame
 
         const EvtVector4R p{ selectedParent->getP4() };
         const EvtVector4R q{ p1 + p2 };
         const EvtVector4R d{ p1 };
 
         const double cost{ EvtDecayAngle( p, q, d ) };
 
         value = acos( cost ) * 180.0 / EvtConst::pi;
 
     } else if ( !selectedVarName.compare( "decayangle3" ) ) {
         // Polar angle between combined first and second daughters
         // with the third daughter in lab frame
 
         const EvtVector4R p{ selectedParent->getP4() };
         const EvtVector4R q{ p1 + p2 + p3 };
         const EvtVector4R d{ p1 + p2 };
 
         const double cost{ EvtDecayAngle( p, q, d ) };
 
         value = acos( cost ) * 180.0 / EvtConst::pi;
 
     } else if ( !selectedVarName.compare( "decayangle_BTO4PI" ) ) {
         // Polar angle between first and second daughters in lab frame.
         // Used in PIPIPI (BTO4PI_CP) model
 
         const EvtVector4R p{ p1 + p2 + p3 };
         const EvtVector4R q{ p1 + p2 };
         const EvtVector4R d{ p1 };
 
         const double cost{ EvtDecayAngle( p, q, d ) };
 
         value = acos( cost ) * 180.0 / EvtConst::pi;
 
     } else if ( !selectedVarName.compare( "chi" ) ) {
         // Chi angle (degrees) using all 4 daughters and parent 4 mom
 
         const EvtParticle* daug1{ selectedParent->getDaug( d1 - 1 ) };
         const EvtParticle* daug2{ selectedParent->getDaug( d1 ) };
         const EvtParticle* daug3{ selectedParent->getDaug( d1 + 1 ) };
         const EvtParticle* daug4{ selectedParent->getDaug( d1 + 2 ) };
 
         const EvtVector4R p_parent{ selectedParent->getP4() };
         const EvtVector4R p_daug1{ daug1 != nullptr ? daug1->getP4()
                                                     : EvtVector4R() };
         const EvtVector4R p_daug2{ daug2 != nullptr ? daug2->getP4()
                                                     : EvtVector4R() };
         const EvtVector4R p_daug3{ daug3 != nullptr ? daug3->getP4()
                                                     : EvtVector4R() };
         const EvtVector4R p_daug4{ daug4 != nullptr ? daug4->getP4()
                                                     : EvtVector4R() };
 
         const double chi{
             EvtDecayAngleChi( p_parent, p_daug1, p_daug2, p_daug3, p_daug4 ) };
         value = chi * 180.0 / EvtConst::pi;
 
     } else if ( !selectedVarName.compare( "cosThetaResNorm" ) ) {
         // P -> R p4 -> (p1 p2 p3) p4, where the resonance R decays to p1 p2 p3.
         // Theta is the angle between the normal of the plane containing p1, p2 & p3
         // and the bachelor particle p4 in the rest frame of resonance R.
         // The normal vector is given by the cross product p3 x p1
 
         if ( sel_NDaugMax > 1 ) {
             const EvtParticle* res{ selectedParent->getDaug( 0 ) };
             const EvtParticle* bac{ selectedParent->getDaug( 1 ) };
 
             // Check resonance has 3 daughters
             if ( res != nullptr && res->getNDaug() == 3 ) {
                 const EvtParticle* daug1 = res->getDaug( 0 );
                 const EvtParticle* daug3 = res->getDaug( 2 );
 
                 // 4-momenta in base parent P lab frame
                 const EvtVector4R pRes{ res->getP4Lab() };
                 const EvtVector4R p4{ bac != nullptr ? bac->getP4Lab()
                                                      : EvtVector4R() };
                 const EvtVector4R p1{ daug1 != nullptr ? daug1->getP4Lab()
                                                        : EvtVector4R() };
                 const EvtVector4R p3{ daug3 != nullptr ? daug3->getP4Lab()
                                                        : EvtVector4R() };
 
                 // Boost 4-vector for resonance frame
                 const EvtVector4R boost{ pRes.get( 0 ), -pRes.get( 1 ),
                                          -pRes.get( 2 ), -pRes.get( 3 ) };
 
                 // Momentum of p1 and p3 in resonance frame
                 const EvtVector4R p1Res{ boostTo( p1, boost ) };
                 const EvtVector4R p3Res{ boostTo( p3, boost ) };
 
                 // Plane normal vector (just uses 3-momentum components)
                 const EvtVector4R norm{ p3Res.cross( p1Res ) };
 
                 // Momentum of p4 in resonance frame
                 const EvtVector4R p4Res{ boostTo( p4, boost ) };
 
                 // Cosine of the angle between the normal and p4 in the resonance frame
                 const double normMag{ norm.d3mag() };
                 const double p4ResMag{ p4Res.d3mag() };
                 if ( normMag > 0.0 && p4ResMag > 0.0 ) {
                     value = norm.dot( p4Res ) / ( normMag * p4ResMag );
                 }
             }
         }
 
     } else if ( !selectedVarName.compare( "cosBetaRes" ) ) {
         // For resonance P (parent) -> p1 p2 p3, beta is the
         // angle between p1 & p3 in the (p1 + p2) rest frame
         if ( sel_NDaugMax > 2 ) {
             const EvtParticle* daug1 = selectedParent->getDaug( 0 );
             const EvtParticle* daug2 = selectedParent->getDaug( 1 );
             const EvtParticle* daug3 = selectedParent->getDaug( 2 );
 
             // 4-momenta in base parent frame
             const EvtVector4R p1{ daug1 != nullptr ? daug1->getP4Lab()
                                                    : EvtVector4R() };
             const EvtVector4R p2{ daug2 != nullptr ? daug2->getP4Lab()
                                                    : EvtVector4R() };
             const EvtVector4R p3{ daug3 != nullptr ? daug3->getP4Lab()
                                                    : EvtVector4R() };
 
             // p1 + p2
             const EvtVector4R p12{ p1 + p2 };
 
             // Boost 4-vector for p12 frame
             const EvtVector4R boost{ p12.get( 0 ), -p12.get( 1 ), -p12.get( 2 ),
                                      -p12.get( 3 ) };
 
             // Momentum of p1 & p3 in p12 frame
             const EvtVector4R p1_12{ boostTo( p1, boost ) };
             const EvtVector4R p3_12{ boostTo( p3, boost ) };
 
             // Cosine of angle between p1 & p3 in p12 frame
             const double p1_12Mag{ p1_12.d3mag() };
             const double p3_12Mag{ p3_12.d3mag() };
             if ( p1_12Mag > 0.0 && p3_12Mag > 0.0 ) {
                 value = p1_12.dot( p3_12 ) / ( p1_12Mag * p3_12Mag );
             }
         }
 
     } else if ( !selectedVarName.compare( "E" ) ) {
         // Energy of first daughter in lab frame
         value = p1_lab.get( 0 );
 
     } else if ( !selectedVarName.compare( "E_over_Eparent" ) ) {
         // Energy of first daughter w.r.t parent energy (lab frame)
         value = p1_lab.get( 0 ) / selectedParent->getP4Lab().get( 0 );
 
     } else if ( !selectedVarName.compare( "E_over_Eparent_over05" ) ) {
         // First daughter E_over_Eparent (lab frame) if d2 granddaughter E ratio > 0.5
         const double E_over_Eparent_2{
             parent->getDaug( 0 )->getDaug( d2 - 1 )->getP4Lab().get( 0 ) /
             parent->getDaug( 0 )->getP4Lab().get( 0 ) };
 
         if ( E_over_Eparent_2 > 0.5 ) {
             value = p1_lab.get( 0 ) / selectedParent->getP4Lab().get( 0 );
         }
 
     } else if ( !selectedVarName.compare( "totE_over_Mparent" ) ) {
         // Energy of given daughters w.r.t parent mass (lab frame)
         value = ( p1_lab.get( 0 ) + p2_lab.get( 0 ) ) / selectedParent->mass();
 
     } else if ( !selectedVarName.compare( "prob" ) ) {
         // Decay probability
         const double* dProb{ selectedParent->decayProb() };
         if ( dProb ) {
             value = *dProb;
         }
 
     } else if ( !selectedVarName.compare( "lifetime" ) ) {
         // Lifetime of particle d1 (ps)
         if ( d1 ) {
             value = par1 != nullptr ? par1->getLifetime() * 1e12 / EvtConst::c
                                     : 0.0;
         } else {
             value = parent != nullptr
                         ? parent->getLifetime() * 1e12 / EvtConst::c
                         : 0.0;
         }
 
     } else if ( !selectedVarName.compare( "deltaT" ) ) {
         // Lifetime difference between particles d1 and d2
         const double t1{
             par1 != nullptr ? par1->getLifetime() * 1e12 / EvtConst::c : 0.0 };
         const double t2{
             par2 != nullptr ? par2->getLifetime() * 1e12 / EvtConst::c : 0.0 };
         value = t1 - t2;
 
     } else if ( !selectedVarName.compare( "decTime" ) ) {
         // Decay flight time of particle d1 in picoseconds.
         // Decay vertex = position of (1st) decay particle of d1
         const EvtParticle* gpar{ par1 != nullptr ? par1->getDaug( 0 ) : nullptr };
         const EvtVector4R vtxPos{ gpar != nullptr ? gpar->get4Pos()
                                                   : EvtVector4R() };
         const double p{ p1_lab.d3mag() };
         value = p > 0.0
                     ? 1e12 * vtxPos.d3mag() * p1_lab.mass() / ( p * EvtConst::c )
                     : 0.0;
     } else if ( !selectedVarName.compare( "nFSRPhotons" ) ) {
         // Loop over all daughters and get number of FSR photons
 
         double nFSRPhotons{ 0 };
 
         for ( size_t iDaughter{ 0 }; iDaughter < selectedParent->getNDaug();
               iDaughter++ ) {
             const EvtParticle* iDaug = selectedParent->getDaug( iDaughter );
 
             if ( iDaug->getAttribute( "FSR" ) == 1 )
                 nFSRPhotons += 1.0;
         }
 
         value = nFSRPhotons;
 
     } else if ( !selectedVarName.compare( "totalFSREnergy" ) ) {
         // Loop over all daughters and get number of FSR photons
 
         double totalFSREnergy{ 0 };
 
         for ( size_t iDaughter{ 0 }; iDaughter < selectedParent->getNDaug();
               iDaughter++ ) {
             const EvtParticle* iDaug = selectedParent->getDaug( iDaughter );
 
             if ( iDaug->getAttribute( "FSR" ) == 1 )
                 totalFSREnergy += iDaug->getP4Lab().get( 0 );
         }
 
         value = totalFSREnergy;
     } else {
         std::cerr << "Warning: Did not recognise variable name "
                   << selectedVarName << std::endl;
     }
 
     return value;
 }
 
 void TestDecayModel::compareHistos( const std::string& refFileName ) const
 {
     // Compare histograms with the same name, calculating the chi-squared
     std::unique_ptr<TFile> refFile{ TFile::Open( refFileName.c_str(), "read" ) };
 
     if ( !refFile ) {
         std::cerr << "Could not open reference file " << refFileName << std::endl;
         return;
     }
 
     // TODO - should we plot the (signed) chisq histogram? and save it as pdf/png?
 
     for ( auto& [_, hist] : m_1DhistVect ) {
         const std::string histName{ hist->GetName() };
         // Get equivalent reference histogram
         const TH1* refHist{
             dynamic_cast<TH1*>( refFile->Get( histName.c_str() ) ) };
 
         if ( refHist ) {
             double chiSq{ 0.0 };
             int nDof{ 0 };
             int iGood{ 0 };
             const double pValue{
                 refHist->Chi2TestX( hist, chiSq, nDof, iGood, "WW" ) };
             const double integral{ refHist->Integral() };
             std::cout << "Histogram " << histName << " chiSq/nDof = " << chiSq
                       << "/" << nDof << ", pValue = " << pValue
                       << ", integral = " << integral << std::endl;
 
         } else {
             std::cerr << "Could not find reference histogram " << histName
                       << std::endl;
         }
     }
 
     for ( auto& [_, hist] : m_2DhistVect ) {
         const std::string histName{ hist->GetName() };
         // Get equivalent reference histogram
         const TH2* refHist{
             dynamic_cast<TH2*>( refFile->Get( histName.c_str() ) ) };
 
         if ( refHist ) {
             double chiSq{ 0.0 };
             int nDof{ 0 };
             int iGood{ 0 };
             const double pValue{
                 refHist->Chi2TestX( hist, chiSq, nDof, iGood, "WW" ) };
             const double integral{ refHist->Integral() };
             std::cout << "Histogram " << histName << " chiSq/nDof = " << chiSq
                       << "/" << nDof << ", pValue = " << pValue
                       << ", integral = " << integral << std::endl;
 
         } else {
             std::cerr << "Could not find reference histogram " << histName
                       << std::endl;
         }
     }
 
     refFile->Close();
 }
 
 double TestDecayModel::getCosAcoplanarityAngle( const EvtParticle* selectedParent,
                                                 const int sel_NDaugMax,
                                                 const int d1, const int d2 ) const
 {
     // Given a two-body decay, the acoplanarity angle is defined as the angle between the
     // two decay planes (normal vectors) in the reference frame of the mother.
     // Each normal vector is the cross product of the momentum of one daughter (in the frame of the mother)
     // and the momentum of one of the granddaughters (in the reference frame of the daughter).
     // In case of 3 daughters, we build an intermediate daughter (2+3) from the last 2 daughters
     // (which are then treated as grand daughters), and in case of 4 daughters, we build
     // 2 intermediate daughters (1+2) and (3+4)
 
     double cosAco{ 0.0 };
 
     if ( sel_NDaugMax < 2 || sel_NDaugMax > 4 ) {
         return cosAco;
     }
 
     const EvtParticle* daughter1{ selectedParent->getDaug( 0 ) };
     const EvtParticle* daughter2{ selectedParent->getDaug( 1 ) };
     const EvtParticle* daughter3{ selectedParent->getDaug( 2 ) };
     const EvtParticle* daughter4{ selectedParent->getDaug( 3 ) };
 
     const EvtParticle* grandDaughter1{ daughter1->getDaug( d1 - 1 ) };
     const EvtParticle* grandDaughter2{ daughter2->getDaug( d2 - 1 ) };
 
     const EvtVector4R parent4Vector{ selectedParent->getP4Lab() };
 
     const EvtVector4R daughter4Vector1{
         sel_NDaugMax <= 3 ? daughter1->getP4Lab()
                           : daughter1->getP4Lab() + daughter2->getP4Lab() };
 
     const EvtVector4R daughter4Vector2{
         sel_NDaugMax == 2   ? daughter2->getP4Lab()
         : sel_NDaugMax == 3 ? daughter2->getP4Lab() + daughter3->getP4Lab()
                             : daughter3->getP4Lab() + daughter4->getP4Lab() };
 
     const EvtVector4R grandDaughter4Vector1{
         sel_NDaugMax <= 3 ? grandDaughter1->getP4Lab() : daughter1->getP4Lab() };
 
     const EvtVector4R grandDaughter4Vector2{
         sel_NDaugMax == 2   ? grandDaughter2->getP4Lab()
         : sel_NDaugMax == 3 ? daughter2->getP4Lab()
                             : daughter3->getP4Lab() };
 
     const EvtVector4R parentBoost{ parent4Vector.get( 0 ),
                                    -parent4Vector.get( 1 ),
                                    -parent4Vector.get( 2 ),
                                    -parent4Vector.get( 3 ) };
 
     const EvtVector4R daughter1Boost{ daughter4Vector1.get( 0 ),
                                       -daughter4Vector1.get( 1 ),
                                       -daughter4Vector1.get( 2 ),
                                       -daughter4Vector1.get( 3 ) };
 
     const EvtVector4R daughter2Boost{ daughter4Vector2.get( 0 ),
                                       -daughter4Vector2.get( 1 ),
                                       -daughter4Vector2.get( 2 ),
                                       -daughter4Vector2.get( 3 ) };
 
     // Boosting daughters to reference frame of the mother
     const EvtVector4R daughter4Vector1_boosted{
         boostTo( daughter4Vector1, parentBoost ) };
     const EvtVector4R daughter4Vector2_boosted{
         boostTo( daughter4Vector2, parentBoost ) };
 
     // Boosting each granddaughter to reference frame of its mother
     const EvtVector4R grandDaughter4Vector1_boosted{
         boostTo( grandDaughter4Vector1, daughter1Boost ) };
     const EvtVector4R grandDaughter4Vector2_boosted{
         boostTo( grandDaughter4Vector2, daughter2Boost ) };
 
     // We calculate the normal vectors of the decay two planes
     const EvtVector4R normalVector1{
         daughter4Vector1_boosted.cross( grandDaughter4Vector1_boosted ) };
     const EvtVector4R normalVector2{
         daughter4Vector2_boosted.cross( grandDaughter4Vector2_boosted ) };
 
     const double normalVector1Mag{ normalVector1.d3mag() };
     const double normalVector2Mag{ normalVector2.d3mag() };
 
     if ( normalVector1Mag > 0.0 && normalVector2Mag > 0.0 ) {
         cosAco = normalVector1.dot( normalVector2 ) /
                  ( normalVector1Mag * normalVector2Mag );
     }
 
     return cosAco;
 }
 
 int main( int argc, char* argv[] )
 {
     if ( argc != 2 ) {
         std::cerr << "Expecting one argument: json input file" << std::endl;
         return 1;
     }
 
     /*! Load input file in json format. */
     json config;
     std::ifstream inputStr{ argv[1] };
     inputStr >> config;
     inputStr.close();
 
     bool allOK{ true };
     if ( config.is_array() ) {
         for ( const auto& cc : config ) {
             TestDecayModel test{ cc };
             allOK &= test.run();
         }
     } else {
         TestDecayModel test{ config };
         allOK &= test.run();
     }
 
     return allOK ? 0 : 1;
 }
diff --git a/test/testDecayModel.hh b/test/testDecayModel.hh
index 7082024..f4804fc 100644
--- a/test/testDecayModel.hh
+++ b/test/testDecayModel.hh
@@ -1,102 +1,104 @@
 
 /***********************************************************************
 * Copyright 1998-2020 CERN for the benefit of the EvtGen authors       *
 *                                                                      *
 * This file is part of EvtGen.                                         *
 *                                                                      *
 * EvtGen is free software: you can redistribute it and/or modify       *
 * it under the terms of the GNU General Public License as published by *
 * the Free Software Foundation, either version 3 of the License, or    *
 * (at your option) any later version.                                  *
 *                                                                      *
 * EvtGen is distributed in the hope that it will be useful,            *
 * but WITHOUT ANY WARRANTY; without even the implied warranty of       *
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the        *
 * GNU General Public License for more details.                         *
 *                                                                      *
 * You should have received a copy of the GNU General Public License    *
 * along with EvtGen.  If not, see <https://www.gnu.org/licenses/>.     *
 ***********************************************************************/
 
 #ifndef TEST_DECAY_MODEL_HH
 #define TEST_DECAY_MODEL_HH
 
 #include "TFile.h"
 #include "TH1D.h"
 #include "TH2D.h"
 
 #include "nlohmann/json.hpp"
 
 #include <cmath>
 #include <string>
 #include <utility>
 #include <vector>
 
 class EvtGen;
 class EvtParticle;
 
 class TestInfo {
   public:
     TestInfo( const std::string& name, const int d1, const int d2,
               const std::string& nameY = "", const int d1Y = 0,
               const int d2Y = 0 ) :
         m_name{ name }, m_nameY{ nameY }, m_d1{ d1 }, m_d2{ d2 }, m_d1Y{ d1Y }, m_d2Y{ d2Y }
     {
     }
 
     const std::string getName( const int var = 1 ) const
     {
         return var == 2 ? m_nameY : m_name;
     }
     int getd1( const int var = 1 ) const { return var == 2 ? m_d1Y : m_d1; }
     int getd2( const int var = 1 ) const { return var == 2 ? m_d2Y : m_d2; }
 
   private:
     std::string m_name;
     std::string m_nameY;
     int m_d1;
     int m_d2;
     int m_d1Y;
     int m_d2Y;
 };
 
 class TestDecayModel {
   public:
     TestDecayModel( const nlohmann::json& config );
 
     bool run();
 
   private:
     bool checkMandatoryFields();
 
     std::string createDecFile(
         const std::string& parent, const std::vector<std::string>& daughterNames,
         const std::vector<std::vector<std::string>>& grandDaughterNames,
         const std::vector<std::string>& models,
         const std::vector<std::vector<std::string>>& parameters,
         const std::vector<bool>& doConjDecay,
         const std::vector<std::string>& extras,
         const std::string decFileName ) const;
 
     void defineHistos( TFile* theFile );
 
     void generateEvents( EvtGen& theGen, const std::string& decFile,
                          const std::string& parentName, const bool doConjDecay,
                          const int nEvents, const bool debugFlag );
 
+    int findChargedDaugtherWithMaxE( const EvtParticle* parent ) const;
+
     double getValue( const EvtParticle* rootPart, const std::string& varName,
                      const int d1, const int d2 ) const;
 
     void compareHistos( const std::string& refFileName ) const;
 
     double getCosAcoplanarityAngle( const EvtParticle* selectedParent,
                                     const int sel_NDaugMax, const int d1,
                                     const int d2 ) const;
 
     const nlohmann::json& m_config;
     std::vector<std::pair<TestInfo, TH1*>> m_1DhistVect;
     std::vector<std::pair<TestInfo, TH2*>> m_2DhistVect;
     TH1* m_mixedHist{ nullptr };
 };
 
 #endif