Page MenuHomeHEPForge
Files F9501186

No OneTemporary
Actions

View Options

diff --git a/EvtGenBase/EvtMTRandomEngine.hh b/EvtGenBase/EvtMTRandomEngine.hh
index 4d713e4..32ce16b 100644
--- a/EvtGenBase/EvtMTRandomEngine.hh
+++ b/EvtGenBase/EvtMTRandomEngine.hh
@@ -1,43 +1,44 @@
/***********************************************************************
* 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 EVTMTRANDOMENGINE_HH
#define EVTMTRANDOMENGINE_HH
#include "EvtGenBase/EvtRandomEngine.hh"
#include <random>
class EvtMTRandomEngine : public EvtRandomEngine {
public:
- EvtMTRandomEngine( unsigned int seed = 1430957218 );
+ EvtMTRandomEngine( unsigned long int seed = 1430957218 );
- virtual double random() override;
+ double random() override;
- virtual void setSeed( unsigned int seed ) override;
+ void setSeed( unsigned long int seed ) override;
+
+ unsigned long int lastSeed() const override { return m_lastSeed; }
private:
std::mt19937 m_engine;
-
- typedef std::uniform_real_distribution<double> URDist;
- URDist m_distribution;
+ std::uniform_real_distribution<double> m_distribution;
+ unsigned long int m_lastSeed;
};
#endif
diff --git a/EvtGenBase/EvtRandom.hh b/EvtGenBase/EvtRandom.hh
index 17582f9..12228b3 100644
--- a/EvtGenBase/EvtRandom.hh
+++ b/EvtGenBase/EvtRandom.hh
@@ -1,48 +1,50 @@
/***********************************************************************
* 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 EVTRANDOM_HH
#define EVTRANDOM_HH
class EvtRandomEngine;
class EvtRandom {
public:
static double Flat();
static double Flat( double max );
static double Flat( double min, double max );
//generate unit Gaussian
static double Gaussian();
static double random();
- static void setSeed( unsigned int seed );
+ static void setSeed( unsigned long int seed );
+
+ static unsigned long int lastSeed();
//This class does not take ownership of the random engine;
//the caller needs to make sure that the engine is not
//destroyed.
static void setRandomEngine( EvtRandomEngine* randomEngine );
private:
static thread_local EvtRandomEngine* m_randomEngine;
};
#endif
diff --git a/EvtGenBase/EvtRandomEngine.hh b/EvtGenBase/EvtRandomEngine.hh
index b4edb89..cbe7c41 100644
--- a/EvtGenBase/EvtRandomEngine.hh
+++ b/EvtGenBase/EvtRandomEngine.hh
@@ -1,39 +1,41 @@
/***********************************************************************
* 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 EVTRANDOMENGINE_HH
#define EVTRANDOMENGINE_HH
// Description:Class to generate random numbers. Single member
// function random is expected to return a random
// number in the range ]0..1[.
class EvtRandomEngine {
public:
virtual ~EvtRandomEngine(){};
virtual double random() = 0;
- virtual void setSeed( unsigned int seed ) = 0;
+ virtual void setSeed( unsigned long int seed ) = 0;
+
+ virtual unsigned long int lastSeed() const = 0;
private:
};
#endif
diff --git a/EvtGenBase/EvtSimpleRandomEngine.hh b/EvtGenBase/EvtSimpleRandomEngine.hh
index e9fb9bc..0730b8a 100644
--- a/EvtGenBase/EvtSimpleRandomEngine.hh
+++ b/EvtGenBase/EvtSimpleRandomEngine.hh
@@ -1,40 +1,42 @@
/***********************************************************************
* 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 EVTSIMPLERANDOMENGINE_HH
#define EVTSIMPLERANDOMENGINE_HH
#include "EvtGenBase/EvtRandomEngine.hh"
class EvtSimpleRandomEngine : public EvtRandomEngine {
public:
EvtSimpleRandomEngine() { m_next = 1; }
void reset() { m_next = 1; }
double random() override;
- void setSeed( unsigned int seed ) override;
+ void setSeed( unsigned long int seed ) override;
+
+ unsigned long int lastSeed() const override { return m_next; }
private:
unsigned long int m_next;
};
#endif
diff --git a/EvtGenExternal/EvtExternalGenFactory.hh b/EvtGenExternal/EvtExternalGenFactory.hh
index 4bc19f4..eb6ce21 100644
--- a/EvtGenExternal/EvtExternalGenFactory.hh
+++ b/EvtGenExternal/EvtExternalGenFactory.hh
@@ -1,69 +1,70 @@
/***********************************************************************
* 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 EVTEXTERNALGENFACTORY_HH
#define EVTEXTERNALGENFACTORY_HH
#include "EvtGenModels/EvtAbsExternalGen.hh"
#include <map>
#include <memory>
// Description: A factory type method to create engines for external physics
// generators like Pythia.
class EvtExternalGenFactory final {
public:
enum class GenId
{
PythiaGenId = 0,
TauolaGenId
};
static EvtExternalGenFactory& getInstance();
EvtAbsExternalGen* getGenerator( const GenId genId );
void initialiseAllGenerators();
void definePythiaGenerator( std::string xmlDir, bool convertPhysCodes,
bool useEvtGenRandom = true );
- void defineTauolaGenerator( bool useEvtGenRandom = true );
+ void defineTauolaGenerator( bool useEvtGenRandom = true,
+ bool seedTauolaFortran = true );
//methods to add configuration commands to the pythia generators
//void addPythiaCommand( std::string generator, std::string module, std::string param, std::string value);
//void addPythia6Command(std::string generator, std::string module, std::string param, std::string value);
private:
EvtExternalGenFactory() = default;
~EvtExternalGenFactory() = default;
EvtExternalGenFactory( const EvtExternalGenFactory& ) = delete;
EvtExternalGenFactory( EvtExternalGenFactory&& ) = delete;
EvtExternalGenFactory& operator=( const EvtExternalGenFactory& ) = delete;
EvtExternalGenFactory& operator=( EvtExternalGenFactory&& ) = delete;
typedef std::map<GenId, std::unique_ptr<EvtAbsExternalGen>> ExtGenMap;
//typedef std::map<GenId, std::map<std::string, std::vector<std::string>>> ExtGenCommandMap;
ExtGenMap m_extGenMap;
//ExtGenCommandMap m_extGenCommandMap;
};
#endif
diff --git a/EvtGenExternal/EvtExternalGenList.hh b/EvtGenExternal/EvtExternalGenList.hh
index 6fbd63e..b85182b 100644
--- a/EvtGenExternal/EvtExternalGenList.hh
+++ b/EvtGenExternal/EvtExternalGenList.hh
@@ -1,57 +1,58 @@
/***********************************************************************
* 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 EVTEXTERNALGENLIST_HH
#define EVTEXTERNALGENLIST_HH
#include "EvtGenBase/EvtAbsRadCorr.hh"
#include "EvtGenBase/EvtDecayBase.hh"
#include <list>
// Description: A factory type method to create engines for external physics
// generators like Pythia.
class EvtExternalGenList {
public:
EvtExternalGenList( bool convertPythiaCodes = false,
std::string pythiaXmlDir = "",
std::string photonType = "gamma",
- bool useEvtGenRandom = true );
+ bool useEvtGenRandom = true,
+ bool seedTauolaFortran = true );
virtual ~EvtExternalGenList();
std::list<EvtDecayBase*> getListOfModels();
EvtAbsRadCorr* getPhotosModel( const double infraredCutOff = 1.0e-7,
const double maxWtInterference = 64.0 );
EvtAbsRadCorr* getSherpaPhotonsModel( const double infraredCutOff = 1.0e-7,
const int mode = 2,
const int useME = 0 );
protected:
private:
std::string m_photonType;
bool m_useEvtGenRandom;
};
#endif
diff --git a/EvtGenExternal/EvtTauolaEngine.hh b/EvtGenExternal/EvtTauolaEngine.hh
index d86d1f2..8eedbdc 100644
--- a/EvtGenExternal/EvtTauolaEngine.hh
+++ b/EvtGenExternal/EvtTauolaEngine.hh
@@ -1,86 +1,87 @@
/***********************************************************************
* 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/>. *
***********************************************************************/
#ifdef EVTGEN_TAUOLA
#ifndef EVTTAUOLAENGINE_HH
#define EVTTAUOLAENGINE_HH
#include "EvtGenBase/EvtDecayBase.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "EvtGenModels/EvtAbsExternalGen.hh"
#ifdef EVTGEN_HEPMC3
#include "HepMC3/Relatives.h"
#include "HepMC3/Units.h"
#include "Tauola/TauolaHepMC3Event.h"
#include "Tauola/TauolaHepMC3Particle.h"
#else
#include "Tauola/TauolaHepMCEvent.h"
#include "Tauola/TauolaHepMCParticle.h"
#include "Tauola/TauolaParticle.h"
#endif
#include "EvtGenBase/EvtHepMCEvent.hh"
#include <mutex>
// Description: Interface to the TAUOLA external generator
class EvtTauolaEngine : public EvtAbsExternalGen {
public:
- EvtTauolaEngine( bool useEvtGenRandom = true );
+ EvtTauolaEngine( bool useEvtGenRandom = true, bool seedTauolaFortran = true );
bool doDecay( EvtParticle* theMother ) override;
void initialise() override;
protected:
private:
GenParticlePtr createGenParticle( const EvtParticle* theParticle ) const;
void setUpPossibleTauModes();
void setOtherParameters();
int getModeInt( EvtDecayBase* decayModel ) const;
void decayTauEvent( EvtParticle* tauParticle );
bool m_useEvtGenRandom{ true };
+ bool m_seedTauolaFortran{ true };
// PDG standard code integer ID for tau particle
static constexpr int m_tauPDG{ 15 };
// Number of possible decay modes in Tauola
static constexpr int m_nTauolaModes{ 22 };
// Neutral and charged spin propagator choices
static int m_neutPropType;
static int m_posPropType;
static int m_negPropType;
static bool m_initialised;
static std::mutex m_tauola_mutex;
};
#endif
#endif
diff --git a/src/EvtGenBase/EvtMTRandomEngine.cpp b/src/EvtGenBase/EvtMTRandomEngine.cpp
index 9433c51..6c176e2 100644
--- a/src/EvtGenBase/EvtMTRandomEngine.cpp
+++ b/src/EvtGenBase/EvtMTRandomEngine.cpp
@@ -1,43 +1,44 @@
/***********************************************************************
* 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/EvtMTRandomEngine.hh"
#include "EvtGenBase/EvtReport.hh"
#include <iostream>
-EvtMTRandomEngine::EvtMTRandomEngine( unsigned int seed ) :
- m_engine( seed ), m_distribution( URDist( 0.0, 1.0 ) )
+EvtMTRandomEngine::EvtMTRandomEngine( unsigned long int seed ) :
+ m_engine{ seed }, m_distribution{ 0.0, 1.0 }, m_lastSeed{ seed }
{
EvtGenReport( EVTGEN_INFO, "EvtMTRandomEngine" )
<< "Mersenne-Twister random number generator with seed = " << seed
<< std::endl;
}
double EvtMTRandomEngine::random()
{
return m_distribution( m_engine );
}
-void EvtMTRandomEngine::setSeed( unsigned int seed )
+void EvtMTRandomEngine::setSeed( unsigned long int seed )
{
m_engine.seed( seed );
+ m_lastSeed = seed;
}
diff --git a/src/EvtGenBase/EvtRandom.cpp b/src/EvtGenBase/EvtRandom.cpp
index 5d9e8da..c6725b6 100644
--- a/src/EvtGenBase/EvtRandom.cpp
+++ b/src/EvtGenBase/EvtRandom.cpp
@@ -1,94 +1,106 @@
/***********************************************************************
* 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/EvtRandom.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtRandomEngine.hh"
#include "EvtGenBase/EvtReport.hh"
#include <iostream>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
using std::endl;
thread_local EvtRandomEngine* EvtRandom::m_randomEngine = nullptr;
void EvtRandom::setRandomEngine( EvtRandomEngine* randomEngine )
{
m_randomEngine = randomEngine;
}
double EvtRandom::random()
{
if ( m_randomEngine == nullptr ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "No random engine available in "
<< "EvtRandom::random()." << endl;
::abort();
}
return m_randomEngine->random();
}
-void EvtRandom::setSeed( unsigned int seed )
+void EvtRandom::setSeed( unsigned long int seed )
{
if ( m_randomEngine == nullptr ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "No random engine available in "
<< "EvtRandom::random()." << endl;
::abort();
}
m_randomEngine->setSeed( seed );
}
+unsigned long int EvtRandom::lastSeed()
+{
+ if ( m_randomEngine == nullptr ) {
+ EvtGenReport( EVTGEN_ERROR, "EvtGen" )
+ << "No random engine available in "
+ << "EvtRandom::random()." << endl;
+ ::abort();
+ }
+
+ return m_randomEngine->lastSeed();
+}
+
// Random number routine to generate numbers between
// min and max. By djl on July 27, 1995.
double EvtRandom::Flat( double min, double max )
{
if ( min > max ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "min>max in EvtRandom::Flat(" << min << "," << max << ")" << endl;
::abort();
}
return EvtRandom::random() * ( max - min ) + min;
}
double EvtRandom::Flat( double max )
{
return max * EvtRandom::random();
}
double EvtRandom::Flat()
{
return EvtRandom::random();
}
double EvtRandom::Gaussian()
{
double x = EvtRandom::random();
double y = EvtRandom::random();
return cos( x * EvtConst::twoPi ) * sqrt( -2.0 * log( 1 - y ) );
}
diff --git a/src/EvtGenBase/EvtSimpleRandomEngine.cpp b/src/EvtGenBase/EvtSimpleRandomEngine.cpp
index 3fcecec..51e665a 100644
--- a/src/EvtGenBase/EvtSimpleRandomEngine.cpp
+++ b/src/EvtGenBase/EvtSimpleRandomEngine.cpp
@@ -1,38 +1,38 @@
/***********************************************************************
* 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/EvtSimpleRandomEngine.hh"
#include <iostream>
#include <math.h>
#include <stdio.h>
double EvtSimpleRandomEngine::random()
{
m_next = m_next * 1103515245 + 123345;
unsigned temp = (unsigned)( m_next / 65536 ) % 32768;
return ( temp + 1.0 ) / 32769.0;
}
-void EvtSimpleRandomEngine::setSeed( unsigned int seed )
+void EvtSimpleRandomEngine::setSeed( unsigned long int seed )
{
m_next = seed;
}
diff --git a/src/EvtGenExternal/EvtExternalGenFactory.cpp b/src/EvtGenExternal/EvtExternalGenFactory.cpp
index 9d11ad4..6b94d76 100644
--- a/src/EvtGenExternal/EvtExternalGenFactory.cpp
+++ b/src/EvtGenExternal/EvtExternalGenFactory.cpp
@@ -1,112 +1,113 @@
/***********************************************************************
* 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 "EvtGenExternal/EvtExternalGenFactory.hh"
#include "EvtGenBase/EvtReport.hh"
#ifdef EVTGEN_PYTHIA
#include "EvtGenExternal/EvtPythiaEngine.hh"
#endif
#ifdef EVTGEN_TAUOLA
#include "EvtGenExternal/EvtTauolaEngine.hh"
#endif
#include <iostream>
using std::endl;
EvtExternalGenFactory& EvtExternalGenFactory::getInstance()
{
static thread_local EvtExternalGenFactory theFactory;
return theFactory;
}
// Only define the generator if we have the external ifdef variable set
#ifdef EVTGEN_PYTHIA
void EvtExternalGenFactory::definePythiaGenerator( std::string xmlDir,
bool convertPhysCodes,
bool useEvtGenRandom )
{
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Defining EvtPythiaEngine: data tables defined in " << xmlDir << endl;
if ( convertPhysCodes == true ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Pythia 6 codes in decay files will be converted to Pythia 8 codes"
<< endl;
} else {
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Pythia 8 codes need to be used in decay files" << endl;
}
if ( useEvtGenRandom == true ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Using EvtGen random engine for Pythia 8 as well" << endl;
}
m_extGenMap[GenId::PythiaGenId] = std::make_unique<EvtPythiaEngine>(
xmlDir, convertPhysCodes, useEvtGenRandom );
}
#else
void EvtExternalGenFactory::definePythiaGenerator( std::string, bool, bool )
{
}
#endif
#ifdef EVTGEN_TAUOLA
-void EvtExternalGenFactory::defineTauolaGenerator( bool useEvtGenRandom )
+void EvtExternalGenFactory::defineTauolaGenerator( bool useEvtGenRandom,
+ bool seedTauolaFortran )
{
EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Defining EvtTauolaEngine." << endl;
- m_extGenMap[GenId::TauolaGenId] = std::make_unique<EvtTauolaEngine>(
- useEvtGenRandom );
+ m_extGenMap[GenId::TauolaGenId] =
+ std::make_unique<EvtTauolaEngine>( useEvtGenRandom, seedTauolaFortran );
}
#else
-void EvtExternalGenFactory::defineTauolaGenerator( bool )
+void EvtExternalGenFactory::defineTauolaGenerator( bool, bool )
{
}
#endif
EvtAbsExternalGen* EvtExternalGenFactory::getGenerator( const GenId genId )
{
ExtGenMap::iterator iter = m_extGenMap.find( genId );
if ( iter == m_extGenMap.end() ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "EvtAbsExternalGen::getGenerator: could not find generator for genId = "
//FIXME C++23 use std::to_underlying
<< static_cast<std::underlying_type_t<GenId>>( genId ) << endl;
return nullptr;
}
// Retrieve the external generator engine
auto& theGenerator = iter->second;
return theGenerator.get();
}
void EvtExternalGenFactory::initialiseAllGenerators()
{
for ( auto& [id, theGenerator] : m_extGenMap ) {
if ( theGenerator ) {
theGenerator->initialise();
}
}
}
diff --git a/src/EvtGenExternal/EvtExternalGenList.cpp b/src/EvtGenExternal/EvtExternalGenList.cpp
index cc7c093..27c4e31 100644
--- a/src/EvtGenExternal/EvtExternalGenList.cpp
+++ b/src/EvtGenExternal/EvtExternalGenList.cpp
@@ -1,120 +1,121 @@
/***********************************************************************
* 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 "EvtGenExternal/EvtExternalGenList.hh"
#include "EvtGenModels/EvtNoRadCorr.hh"
#include "EvtGenExternal/EvtExternalGenFactory.hh"
#include "EvtGenExternal/EvtPHOTOS.hh"
#include "EvtGenExternal/EvtPythia.hh"
#include "EvtGenExternal/EvtSherpaPhotons.hh"
#include "EvtGenExternal/EvtTauola.hh"
EvtExternalGenList::EvtExternalGenList( bool convertPythiaCodes,
std::string pythiaXmlDir,
std::string photonType,
- bool useEvtGenRandom ) :
+ bool useEvtGenRandom,
+ bool seedTauolaFortran ) :
m_photonType{ photonType }, m_useEvtGenRandom{ useEvtGenRandom }
{
// Instantiate the external generator factory
EvtExternalGenFactory& extFactory = EvtExternalGenFactory::getInstance();
if ( pythiaXmlDir.size() < 1 ) {
// If we have no string defined, check the value of the
// PYTHIA8DATA environment variable which should be set to the
// xmldoc Pythia directory
char* pythiaDataDir = getenv( "PYTHIA8DATA" );
if ( pythiaDataDir != nullptr ) {
pythiaXmlDir = pythiaDataDir;
}
}
extFactory.definePythiaGenerator( pythiaXmlDir, convertPythiaCodes,
useEvtGenRandom );
- extFactory.defineTauolaGenerator( useEvtGenRandom );
+ extFactory.defineTauolaGenerator( useEvtGenRandom, seedTauolaFortran );
}
EvtExternalGenList::~EvtExternalGenList()
{
}
#ifdef EVTGEN_PHOTOS
EvtAbsRadCorr* EvtExternalGenList::getPhotosModel( const double infraredCutOff,
const double maxWtInterference )
{
// Define the Photos model, which uses the EvtPHOTOS class.
EvtPHOTOS* photosModel = new EvtPHOTOS( m_photonType, m_useEvtGenRandom,
infraredCutOff, maxWtInterference );
return photosModel;
}
#else
EvtAbsRadCorr* EvtExternalGenList::getPhotosModel(
const double /*infraredCutOff*/, const double /*maxWtInterference*/ )
{
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< " PHOTOS generator has been called for FSR simulation, but it was not switched on during compilation."
<< std::endl;
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< " The simulation will be generated without FSR." << std::endl;
return new EvtNoRadCorr{};
}
#endif
#ifdef EVTGEN_SHERPA
EvtAbsRadCorr* EvtExternalGenList::getSherpaPhotonsModel(
const double infraredCutOff, const int mode, const int useME )
{
// Define the Photos model, which uses the EvtSherpaPhotonsEngine class.
EvtSherpaPhotons* sherpaPhotonsModel =
new EvtSherpaPhotons( m_useEvtGenRandom, infraredCutOff, mode, useME );
return sherpaPhotonsModel;
}
#else
EvtAbsRadCorr* EvtExternalGenList::getSherpaPhotonsModel(
const double /*infraredCutOff*/, const int /*mode*/, const int /*useME*/ )
{
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< " Sherpa's PHOTONS++ generator has been called for FSR simulation, but Sherpa was not switched on during compilation."
<< std::endl;
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< " The simulation will be generated without FSR." << std::endl;
return new EvtNoRadCorr{};
}
#endif
std::list<EvtDecayBase*> EvtExternalGenList::getListOfModels()
{
// Create the Pythia and Tauola models, which use their own engine classes.
EvtPythia* pythiaModel = new EvtPythia();
EvtTauola* tauolaModel = new EvtTauola();
std::list<EvtDecayBase*> extraModels;
extraModels.push_back( pythiaModel );
extraModels.push_back( tauolaModel );
// Return the list of models
return extraModels;
}
diff --git a/src/EvtGenExternal/EvtTauolaEngine.cpp b/src/EvtGenExternal/EvtTauolaEngine.cpp
index b444289..ecc7e87 100644
--- a/src/EvtGenExternal/EvtTauolaEngine.cpp
+++ b/src/EvtGenExternal/EvtTauolaEngine.cpp
@@ -1,591 +1,601 @@
/***********************************************************************
* 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/>. *
***********************************************************************/
#ifdef EVTGEN_TAUOLA
#include "EvtGenExternal/EvtTauolaEngine.hh"
#include "EvtGenBase/EvtDecayTable.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtReport.hh"
#include "EvtGenBase/EvtSymTable.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include "Tauola/Log.h"
#include "Tauola/Tauola.h"
#include <array>
#include <cmath>
#include <iostream>
+#include <limits>
#include <memory>
#include <sstream>
#include <string>
using std::endl;
// Mutex Tauola as it is not thread safe.
int EvtTauolaEngine::m_neutPropType = 0;
int EvtTauolaEngine::m_posPropType = 0;
int EvtTauolaEngine::m_negPropType = 0;
bool EvtTauolaEngine::m_initialised = false;
std::mutex EvtTauolaEngine::m_tauola_mutex;
-EvtTauolaEngine::EvtTauolaEngine( bool useEvtGenRandom ) :
- m_useEvtGenRandom{ useEvtGenRandom }
+EvtTauolaEngine::EvtTauolaEngine( bool useEvtGenRandom, bool seedTauolaFortran ) :
+ m_useEvtGenRandom{ useEvtGenRandom }, m_seedTauolaFortran{ seedTauolaFortran }
{
}
void EvtTauolaEngine::initialise()
{
const std::lock_guard<std::mutex> lock( m_tauola_mutex );
// Set up all possible tau decay modes.
// This should be done just before the first doDecay() call,
// since we want to make sure that any decay.dec files are processed
// first to get lists of particle modes and their alias definitions
// (for creating EvtParticles with the right history information).
if ( m_initialised ) {
return;
}
EvtGenReport( EVTGEN_INFO, "EvtGen" ) << "Setting up TAUOLA." << endl;
// These three lines are not really necessary since they are the default.
// But they are here so that we know what the initial conditions are.
// tau PDG code
Tauolapp::Tauola::setDecayingParticle( m_tauPDG );
// all modes allowed
Tauolapp::Tauola::setSameParticleDecayMode( Tauolapp::Tauola::All );
// all modes allowed
Tauolapp::Tauola::setOppositeParticleDecayMode( Tauolapp::Tauola::All );
// Limit the number of warnings printed out.
// Can't choose zero here, unfortunately.
Tauolapp::Log::SetWarningLimit( 1 );
// Initial the Tauola external generator
if ( m_useEvtGenRandom ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Using EvtGen random number engine also for Tauola++" << endl;
Tauolapp::Tauola::setRandomGenerator( EvtRandom::Flat );
}
// Use the BaBar-tuned chiral current calculations by default. Can be changed using the
// TauolaCurrentOption keyword in decay files
Tauolapp::Tauola::setNewCurrents( 1 );
Tauolapp::Tauola::initialize();
// Set-up possible decay modes _after_ we have read the (user) decay file
this->setUpPossibleTauModes();
this->setOtherParameters();
m_initialised = true;
}
void EvtTauolaEngine::setUpPossibleTauModes()
{
// Get the decay table list defined by the decay.dec files.
// Only look for the first tau particle decay mode definitions with the Tauola name,
// since that generator only allows the same BFs for both tau+ and tau- decays.
// We can not choose a specific tau decay event-by-event, since this is
// only possible before we call Tauola::initialize().
// Otherwise, we could have selected a random mode ourselves for tau- and tau+
// separately (via selecting a random number and comparing it to be less than
// the cumulative BF) for each event.
const int nPDL = EvtPDL::entries();
bool gotAnyTauolaModes( false );
for ( int iPDL = 0; iPDL < nPDL; iPDL++ ) {
const EvtId particleId = EvtPDL::getEntry( iPDL );
const int PDGId = EvtPDL::getStdHep( particleId );
if ( abs( PDGId ) == m_tauPDG && gotAnyTauolaModes == false ) {
const int aliasInt = particleId.getAlias();
// Get the list of decay modes for this tau particle (alias)
const int nModes = EvtDecayTable::getInstance().getNModes( aliasInt );
// Vector to store tau mode branching fractions.
// The size of this vector equals the total number of possible
// Tauola decay modes. Initialise all BFs to zero.
std::vector<double> tauolaModeBFs;
tauolaModeBFs.assign( m_nTauolaModes, 0.0 );
double totalTauModeBF( 0.0 );
int nNonTauolaModes( 0 );
// Loop through each decay mode
for ( int iMode = 0; iMode < nModes; iMode++ ) {
EvtDecayBase* decayModel =
EvtDecayTable::getInstance().findDecayModel( aliasInt, iMode );
if ( decayModel ) {
// Check that the decay model name matches TAUOLA
std::string modelName = decayModel->getName();
if ( modelName == "TAUOLA" ) {
gotAnyTauolaModes = true;
// Extract the decay mode integer type and branching fraction
double BF = decayModel->getBranchingFraction();
int modeArrayInt = this->getModeInt( decayModel ) - 1;
if ( modeArrayInt >= 0 && modeArrayInt < m_nTauolaModes ) {
tauolaModeBFs[modeArrayInt] = BF;
totalTauModeBF += BF;
}
} else {
nNonTauolaModes++;
}
} // Decay mode exists
} // Loop over decay models
if ( gotAnyTauolaModes == true && nNonTauolaModes > 0 ) {
EvtGenReport( EVTGEN_ERROR, "EvtGen" )
<< "Please remove all non-TAUOLA decay modes for particle "
<< EvtPDL::name( particleId ) << endl;
::abort();
}
// Normalise all (non-zero) tau mode BFs to sum up to 1.0, and
// let Tauola know about these normalised branching fractions
if ( totalTauModeBF > 0.0 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Setting TAUOLA BF modes using the definitions for the particle "
<< EvtPDL::name( particleId ) << endl;
for ( int iTauMode = 0; iTauMode < m_nTauolaModes; iTauMode++ ) {
tauolaModeBFs[iTauMode] /= totalTauModeBF;
double modeBF = tauolaModeBFs[iTauMode];
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Setting TAUOLA BF for mode " << iTauMode + 1
<< " = " << modeBF << endl;
Tauolapp::Tauola::setTauBr( iTauMode + 1, modeBF );
}
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "Any other TAUOLA BF modes for other tau particle decay mode definitions will be ignored!"
<< endl;
}
} // Got tau particle and have yet to get a TAUOLA mode
} // Loop through PDL entries
}
int EvtTauolaEngine::getModeInt( EvtDecayBase* decayModel ) const
{
int modeInt( 0 );
if ( decayModel ) {
int nVars = decayModel->getNArg();
if ( nVars > 0 ) {
modeInt = static_cast<int>( decayModel->getArg( 0 ) );
}
}
return modeInt;
}
void EvtTauolaEngine::setOtherParameters()
{
// Set other Tauola parameters using the "Defined" keyword in the decay file. If any of
// these are not found in the decay file, then default values are assumed/kept
// 1) TauolaNeutralProp: Specify the neutral propagator type used for spin matrix calculations
// "Z" (default), "Gamma", "Higgs" (H0), "PseudoHiggs" (A0), "MixedHiggs" (A0/H0)
int iErr( 0 );
std::string neutPropName = EvtSymTable::get( "TauolaNeutralProp", iErr );
if ( neutPropName == "Z0" || neutPropName == "Z" ) {
m_neutPropType = Tauolapp::TauolaParticle::Z0;
} else if ( neutPropName == "Gamma" ) {
m_neutPropType = Tauolapp::TauolaParticle::GAMMA;
} else if ( neutPropName == "Higgs" ) {
m_neutPropType = Tauolapp::TauolaParticle::HIGGS;
} else if ( neutPropName == "PseudoHiggs" ) {
m_neutPropType = Tauolapp::TauolaParticle::HIGGS_A;
} else if ( neutPropName == "MixedHiggs" ) {
m_neutPropType = Tauolapp::Tauola::getHiggsScalarPseudoscalarPDG();
}
if ( m_neutPropType != 0 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "TAUOLA neutral spin propagator PDG id set to " << m_neutPropType
<< endl;
}
// 2) TauolaChargedProp: Specify the charged propagator type used for spin matrix calculations
// "W" (default), "Higgs" (H+/H-)
std::string chargedPropName = EvtSymTable::get( "TauolaChargedProp", iErr );
if ( chargedPropName == "W" ) {
m_negPropType = Tauolapp::TauolaParticle::W_MINUS;
m_posPropType = Tauolapp::TauolaParticle::W_PLUS;
} else if ( chargedPropName == "Higgs" ) {
m_negPropType = Tauolapp::TauolaParticle::HIGGS_MINUS;
m_posPropType = Tauolapp::TauolaParticle::HIGGS_PLUS;
}
if ( m_negPropType != 0 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "TAUOLA negative charge spin propagator PDG id set to "
<< m_negPropType << endl;
}
if ( m_posPropType != 0 ) {
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "TAUOLA positive charge spin propagator PDG id set to "
<< m_posPropType << endl;
}
// 3) TauolaHiggsMixingAngle: Specify the mixing angle between the neutral scalar & pseudoscalar Higgs
// A0/H0; the default mixing angle is pi/4 radians
std::string mixString = EvtSymTable::get( "TauolaHiggsMixingAngle", iErr );
// If the definition name is not found, get() just returns the first argument string
if ( mixString != "TauolaHiggsMixingAngle" ) {
double mixAngle = std::atof( mixString.c_str() );
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "TAUOLA Higgs mixing angle set to " << mixAngle << " radians"
<< endl;
Tauolapp::Tauola::setHiggsScalarPseudoscalarMixingAngle( mixAngle );
}
// 4) TauolaBRi, where i = 1,2,3,4: Redefine sub-channel branching fractions using the setTaukle
// function, after initialized() has been called. Default values = 0.5, 0.5, 0.5 and 0.6667
std::array<double, 4> BRVect{ 0.5, 0.5, 0.5, 0.6667 };
for ( int j = 0; j < 4; j++ ) {
std::ostringstream o;
o << j + 1;
std::string BRName = "TauolaBR" + o.str();
std::string stringBR = EvtSymTable::get( BRName, iErr );
// If the definition name is not found, get() just returns the first argument string
if ( stringBR != BRName ) {
BRVect[j] = std::atof( stringBR.c_str() );
}
}
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "TAUOLA::setTaukle values are " << BRVect[0] << ", " << BRVect[1]
<< ", " << BRVect[2] << ", " << BRVect[3] << endl;
Tauolapp::Tauola::setTaukle( BRVect[0], BRVect[1], BRVect[2], BRVect[3] );
// 5) Specify the hadronic current option, e.g. orig CLEO = 0, BaBar-tuned = 1 (default), ...
// No check is made by EvtGen on valid integer options - its just passed to Tauola
std::string currentOption = EvtSymTable::get( "TauolaCurrentOption", iErr );
// If the definition name is not found, get() just returns the first argument string
if ( currentOption != "TauolaCurrentOption" ) {
int currentOpt = std::atoi( currentOption.c_str() );
EvtGenReport( EVTGEN_INFO, "EvtGen" )
<< "TAUOLA current option = " << currentOpt << endl;
Tauolapp::Tauola::setNewCurrents( currentOpt );
}
}
bool EvtTauolaEngine::doDecay( EvtParticle* tauParticle )
{
if ( !m_initialised ) {
this->initialise();
}
if ( !tauParticle ) {
return false;
}
// Check that we have a tau particle.
EvtId partId = tauParticle->getId();
if ( abs( EvtPDL::getStdHep( partId ) ) != m_tauPDG ) {
return false;
}
int nTauDaug = tauParticle->getNDaug();
// If the number of tau daughters is not zero, then we have already decayed
// it using Tauola/another decay algorithm.
if ( nTauDaug > 0 ) {
return true;
}
this->decayTauEvent( tauParticle );
return true;
}
void EvtTauolaEngine::decayTauEvent( EvtParticle* tauParticle )
{
// Either we have a tau particle within a decay chain, or a single particle.
// Create a dummy HepMC event & vertex for the parent particle, containing the tau as
// one of the outgoing particles. If we have a decay chain, the parent will be the
// incoming particle, while the daughters, including the tau, are outgoing particles.
// For the single particle case, the incoming particle is null, while the single tau
// is the only outgoing particle.
// We can then pass this event to Tauola which should then decay the tau particle.
// We also consider all other tau particles from the parent decay in the logic below.
// Create the dummy event.
auto theEvent = std::make_unique<GenEvent>( Units::GEV, Units::MM );
// Create the decay "vertex".
GenVertexPtr theVertex = newGenVertexPtr();
theEvent->add_vertex( theVertex );
// Get the parent of this tau particle
EvtParticle* theParent = tauParticle->getParent();
GenParticlePtr hepMCParent( nullptr );
// Assign the parent particle as the incoming particle to the vertex.
if ( theParent ) {
hepMCParent = this->createGenParticle( theParent );
theVertex->add_particle_in( hepMCParent );
} else {
// The tau particle has no parent. Set "itself" as the incoming particle for the first vertex.
// This is needed, otherwise Tauola warns of momentum non-conservation for this (1st) vertex.
GenParticlePtr tauGenInit = this->createGenParticle( tauParticle );
theVertex->add_particle_in( tauGenInit );
}
// Find all daughter particles and assign them as outgoing particles to the vertex.
// This will include the tau particle we are currently processing.
// If the parent decay has more than one tau particle, we need to include them as well.
// This is important since Tauola needs the correct physics correlations: we do not
// want Tauola to decay each particle separately if they are from tau pair combinations.
// Tauola will process the event, and we will create EvtParticles from all tau decay
// products, i.e. the tau particle we currently have and any other tau particles.
// EvtGen will try to decay the other tau particle(s) by calling EvtTauola and therefore
// this function. However, we check to see if the tau candidate has any daughters already.
// If it does, then we have already set the tau decay products from Tauola.
// Map to store (HepMC,EvtParticle) pairs for each tau candidate from the parent
// decay. This is needed to find out what EvtParticle corresponds to a given tau HepMC
// candidate: we do not want to recreate existing EvtParticle pointers.
std::map<GenParticlePtr, EvtParticle*> tauMap;
// Keep track of the original EvtId of the parent particle, since we may need to set
// the equivalent HepMCParticle has a gauge boson to let Tauola calculate spin effects
EvtId origParentId( -1, -1 );
if ( theParent ) {
// Original parent id
origParentId = EvtPDL::getId( theParent->getName() );
// Find all tau particles in the decay tree and store them in the map.
// Keep track of how many tau daughters this parent particle has
int nTaus( 0 );
int nDaug( theParent->getNDaug() );
int iDaug( 0 );
for ( iDaug = 0; iDaug < nDaug; iDaug++ ) {
EvtParticle* theDaughter = theParent->getDaug( iDaug );
if ( theDaughter ) {
GenParticlePtr hepMCDaughter = this->createGenParticle(
theDaughter );
theVertex->add_particle_out( hepMCDaughter );
EvtId theId = theDaughter->getId();
int PDGInt = EvtPDL::getStdHep( theId );
if ( abs( PDGInt ) == m_tauPDG ) {
// Delete any siblings for the tau particle
if ( theDaughter->getNDaug() > 0 ) {
theDaughter->deleteDaughters( false );
}
tauMap[hepMCDaughter] = theDaughter;
nTaus++;
} else {
// Treat all other particles as "stable"
hepMCDaughter->set_status( Tauolapp::TauolaParticle::STABLE );
}
} // theDaughter != 0
} // Loop over daughters
// For the parent particle, artifically set the PDG to a boson with the same 4-momentum
// so that spin correlations are calculated inside Tauola.
// This leaves the original parent m_EvtParticle_ unchanged
if ( nTaus > 0 && hepMCParent ) {
int parCharge = EvtPDL::chg3( origParentId ) /
3; // (3*particle charge)/3 = particle charge
if ( parCharge == 0 && m_neutPropType != 0 ) {
hepMCParent->set_pdg_id( m_neutPropType );
} else if ( parCharge == -1 && m_negPropType != 0 ) {
hepMCParent->set_pdg_id( m_negPropType );
} else if ( parCharge == 1 && m_posPropType != 0 ) {
hepMCParent->set_pdg_id( m_posPropType );
}
}
} else {
// We only have the one tau particle. Store only this in the map.
GenParticlePtr singleTau = this->createGenParticle( tauParticle );
theVertex->add_particle_out( singleTau );
tauMap[singleTau] = tauParticle;
}
{
const std::lock_guard<std::mutex> lock( m_tauola_mutex );
+ if ( m_useEvtGenRandom && m_seedTauolaFortran ) {
+ static thread_local auto lastSeed{
+ std::numeric_limits<unsigned long int>::max() };
+ if ( lastSeed != EvtRandom::lastSeed() ) {
+ lastSeed = EvtRandom::lastSeed();
+ Tauolapp::Tauola::setSeed( lastSeed, 0, 0 );
+ }
+ }
+
// Now pass the event to Tauola for processing
// Create a Tauola event object
#ifdef EVTGEN_HEPMC3
Tauolapp::TauolaHepMC3Event tauolaEvent( theEvent.get() );
#else
Tauolapp::TauolaHepMCEvent tauolaEvent( theEvent.get() );
#endif
// Run the Tauola algorithm
tauolaEvent.decayTaus();
}
// Loop over all tau particles in the HepMC event and create their EvtParticle daughters.
// Store all final "stable" descendent particles as the tau daughters, i.e.
// let Tauola decay any resonances such as a_1 or rho.
// If there is more than one tau particle in the event, then also create the
// corresponding EvtParticles for their daughters as well. They will not be
// re-decayed since we check at the start of this function if the tau particle has
// any daughters before running Tauola decayTaus().
#ifdef EVTGEN_HEPMC3
for ( auto aParticle : theEvent->particles() ) {
#else
HepMC::GenEvent::particle_iterator eventIter;
for ( eventIter = theEvent->particles_begin();
eventIter != theEvent->particles_end(); ++eventIter ) {
// Check to see if we have a tau particle
HepMC::GenParticle* aParticle = ( *eventIter );
#endif
if ( aParticle && abs( aParticle->pdg_id() ) == m_tauPDG ) {
// Find out what EvtParticle corresponds to the HepMC particle.
// We need this to create and attach EvtParticle daughters.
EvtParticle* tauEvtParticle = tauMap[aParticle];
if ( tauEvtParticle ) {
// Get the tau 4-momentum in the lab (first mother) frame. We need to boost
// all the tau daughters to this frame, such that daug.getP4() is in the tau restframe.
EvtVector4R tauP4CM = tauEvtParticle->getP4Lab();
tauP4CM.set( tauP4CM.get( 0 ), -tauP4CM.get( 1 ),
-tauP4CM.get( 2 ), -tauP4CM.get( 3 ) );
// Get the decay vertex for the tau particle
GenVertexPtr endVertex = aParticle->end_vertex();
std::vector<EvtId> daugIdVect;
std::vector<EvtVector4R> daugP4Vect;
// Loop through all descendants
#ifdef EVTGEN_HEPMC3
for ( auto tauDaug :
HepMC3::Relatives::DESCENDANTS( endVertex ) ) {
#else
HepMC::GenVertex::particle_iterator tauIter;
// Loop through all descendants
for ( tauIter = endVertex->particles_begin( HepMC::descendants );
tauIter != endVertex->particles_end( HepMC::descendants );
++tauIter ) {
HepMC::GenParticle* tauDaug = ( *tauIter );
#endif
// Check to see if this descendant has its own decay vertex, e.g. rho resonance.
// If so, skip this daughter and continue looping through the descendant list
// until we reach the final "stable" products (e.g. pi pi from rho -> pi pi).
GenVertexPtr daugDecayVtx = tauDaug->end_vertex();
if ( daugDecayVtx ) {
continue;
}
// Store the particle id and 4-momentum
int tauDaugPDG = tauDaug->pdg_id();
EvtId daugId = EvtPDL::evtIdFromStdHep( tauDaugPDG );
daugIdVect.push_back( daugId );
FourVector tauDaugP4 = tauDaug->momentum();
double tauDaug_px = tauDaugP4.px();
double tauDaug_py = tauDaugP4.py();
double tauDaug_pz = tauDaugP4.pz();
double tauDaug_E = tauDaugP4.e();
EvtVector4R daugP4( tauDaug_E, tauDaug_px, tauDaug_py,
tauDaug_pz );
daugP4Vect.push_back( daugP4 );
} // Loop over HepMC tau daughters
// Create the tau EvtParticle daughters and assign their ids and 4-mtm
int nDaug = daugIdVect.size();
tauEvtParticle->makeDaughters( nDaug, daugIdVect );
int iDaug( 0 );
for ( iDaug = 0; iDaug < nDaug; iDaug++ ) {
EvtParticle* theDaugPart = tauEvtParticle->getDaug( iDaug );
if ( theDaugPart ) {
EvtId theDaugId = daugIdVect[iDaug];
EvtVector4R theDaugP4 = daugP4Vect[iDaug];
theDaugP4.applyBoostTo(
tauP4CM ); // Boost the 4-mtm to the tau rest frame
theDaugPart->init( theDaugId, theDaugP4 );
}
} // Loop over tau daughters
}
} // We have a tau HepMC particle in the event
}
theEvent->clear();
}
GenParticlePtr EvtTauolaEngine::createGenParticle( const EvtParticle* theParticle ) const
{
// Method to create an HepMC::GenParticle version of the given EvtParticle.
if ( theParticle == nullptr ) {
return nullptr;
}
// Get the 4-momentum (E, px, py, pz) for the EvtParticle
EvtVector4R p4 = theParticle->getP4Lab();
// Convert this to the HepMC 4-momentum
double E = p4.get( 0 );
double px = p4.get( 1 );
double py = p4.get( 2 );
double pz = p4.get( 3 );
FourVector hepMC_p4( px, py, pz, E );
int PDGInt = EvtPDL::getStdHep( theParticle->getId() );
// Set the status flag for the particle.
int status = Tauolapp::TauolaParticle::HISTORY;
GenParticlePtr genParticle = newGenParticlePtr( hepMC_p4, PDGInt, status );
return genParticle;
}
#endif
diff --git a/test/testDecayModel.cc b/test/testDecayModel.cc
index 41d50be..933ae7f 100644
--- a/test/testDecayModel.cc
+++ b/test/testDecayModel.cc
@@ -1,1828 +1,1829 @@
/***********************************************************************
* 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 "TROOT.h"
#include "tbb/tbb.h"
#include <chrono>
#include <fstream>
#include <future>
#include <iostream>
#include <list>
#include <memory>
using nlohmann::json;
std::once_flag TestDecayModel::m_createDecFile_threadlock;
TestHistos::TestHistos( const std::string& parentName, const json& config )
{
// Histogram information
const std::size_t nHistos{ config.size() };
m_1DhistVect.reserve( nHistos );
m_2DhistVect.reserve( nHistos );
for ( const auto& hInfo : config ) {
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" ) ) {
auto hist = std::make_unique<TH1D>( histName.c_str(),
varTitle.c_str(), nBins, xmin,
xmax );
m_1DhistVect.emplace_back(
std::make_pair( HistInfo{ varName, d1, d2 }, std::move( hist ) ) );
} 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 );
}
auto hist = std::make_unique<TH2D>( histName.c_str(),
varTitle.c_str(), nBins, xmin,
xmax, nBinsY, ymin, ymax );
m_2DhistVect.emplace_back(
std::make_pair( HistInfo{ varName, d1, d2, varNameY, d1Y, d2Y },
std::move( hist ) ) );
}
}
// Add a mixed/unmixed histogram
// Useful for the case where the parent is a neutral K, D or B
const std::array<std::string, 14> parentsThatMix{
"B_s0", "anti-B_s0", "B_s0L", "B_s0H", "B0", "anti-B0", "B0L",
"B0H", "D0", "anti-D0", "K0", "anti-K0", "K_S0", "K_L0" };
if ( std::find( parentsThatMix.begin(), parentsThatMix.end(), parentName ) !=
parentsThatMix.end() ) {
const std::string varTitle{ parentName + " mixed" };
m_mixedHist = std::make_unique<TH1D>( "mixed", varTitle.c_str(), 2, 0.0,
2.0 );
// TODO maybe set bin labels?
}
}
TestHistos::TestHistos( const TestHistos& rhs )
{
m_1DhistVect.reserve( rhs.m_1DhistVect.size() );
for ( auto& [info, hist] : rhs.m_1DhistVect ) {
auto newHist = std::unique_ptr<TH1>{ static_cast<TH1*>( hist->Clone() ) };
m_1DhistVect.push_back( std::make_pair( info, std::move( newHist ) ) );
}
m_2DhistVect.reserve( rhs.m_2DhistVect.size() );
for ( auto& [info, hist] : rhs.m_2DhistVect ) {
auto newHist = std::unique_ptr<TH2>{ static_cast<TH2*>( hist->Clone() ) };
m_2DhistVect.push_back( std::make_pair( info, std::move( newHist ) ) );
}
if ( rhs.m_mixedHist ) {
m_mixedHist.reset( static_cast<TH1*>( rhs.m_mixedHist->Clone() ) );
}
}
TestHistos::TestHistos( TestHistos&& rhs ) noexcept
{
this->swap( rhs );
}
TestHistos& TestHistos::operator=( const TestHistos& rhs )
{
TestHistos tmp{ rhs };
this->swap( tmp );
return *this;
}
TestHistos& TestHistos::operator=( TestHistos&& rhs ) noexcept
{
this->swap( rhs );
return *this;
}
void TestHistos::swap( TestHistos& rhs ) noexcept
{
m_1DhistVect.swap( rhs.m_1DhistVect );
m_2DhistVect.swap( rhs.m_2DhistVect );
std::swap( m_mixedHist, rhs.m_mixedHist );
}
void TestHistos::add( const TestHistos& rhs )
{
// handle the special case where we have been default constructed and the rhs has not
if ( m_1DhistVect.empty() && m_2DhistVect.empty() && !m_mixedHist ) {
( *this ) = rhs;
return;
}
// TODO - should really check that the sets of histograms are the same between left and right
const std::size_t n1DHists{ rhs.m_1DhistVect.size() };
for ( std::size_t i{ 0 }; i < n1DHists; ++i ) {
m_1DhistVect[i].second->Add( rhs.m_1DhistVect[i].second.get() );
}
const std::size_t n2DHists{ rhs.m_2DhistVect.size() };
for ( std::size_t i{ 0 }; i < n2DHists; ++i ) {
m_2DhistVect[i].second->Add( rhs.m_2DhistVect[i].second.get() );
}
if ( m_mixedHist && rhs.m_mixedHist ) {
m_mixedHist->Add( rhs.m_mixedHist.get() );
}
}
void TestHistos::normalise()
{
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 );
}
}
if ( m_mixedHist ) {
const double area{ m_mixedHist->Integral() };
if ( area > 0.0 ) {
m_mixedHist->Scale( 1.0 / area );
}
}
}
void TestHistos::save( TFile* outputFile )
{
outputFile->cd();
for ( auto& [info, hist] : m_1DhistVect ) {
hist->SetDirectory( outputFile );
hist->Write();
hist.release();
}
for ( auto& [info, hist] : m_2DhistVect ) {
hist->SetDirectory( outputFile );
hist->Write();
hist.release();
}
if ( m_mixedHist ) {
m_mixedHist->SetDirectory( outputFile );
m_mixedHist->Write();
m_mixedHist.release();
}
}
TestDecayModel::TestDecayModel( const json& config ) :
m_config{ checkMandatoryFields( config )
? readConfig( config )
: throw std::runtime_error{
"ERROR : json does not contain all required fields" } }
{
}
bool TestDecayModel::checkMandatoryFields( const json& config )
{
const std::array<std::string, 7> mandatoryFields{ "parent", "daughters",
"models", "parameters",
"outfile", "events",
"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 ( !config.contains( field ) ) {
std::cerr << "ERROR : json does not contain required field: " << field
<< std::endl;
allMandatoryFields = false;
continue;
}
if ( field == "histograms" ) {
const json& jHistos{ 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;
}
TestConfig TestDecayModel::readConfig( const json& config )
{
TestConfig cfg;
// Get all the mandatory fields first
cfg.parentName = config.at( "parent" ).get<std::string>();
cfg.daughterNames = config.at( "daughters" ).get<std::vector<std::string>>();
cfg.modelNames = config.at( "models" ).get<std::vector<std::string>>();
cfg.modelParameters =
config.at( "parameters" ).get<std::vector<std::vector<std::string>>>();
cfg.nEvents = config.at( "events" ).get<std::size_t>();
// Histogram information
cfg.testHistograms = TestHistos{ cfg.parentName, config.at( "histograms" ) };
// Then check for optional fields, setting default values if not present
if ( config.contains( "grand_daughters" ) &&
config.at( "grand_daughters" ).is_array() ) {
cfg.grandDaughterNames = config.at( "grand_daughters" )
.get<std::vector<std::vector<std::string>>>();
}
if ( config.contains( "extras" ) && config.at( "extras" ).is_array() ) {
cfg.extraCommands = config.at( "extras" ).get<std::vector<std::string>>();
}
// Set the number of threads to use, 1 by default
cfg.nThreads = 1;
if ( config.contains( "threads" ) ) {
cfg.nThreads = config.at( "threads" ).get<std::size_t>();
}
// Set the FSR generator. Use PHOTOS by default.
cfg.fsrGenerator = FSRGenerator::PHOTOS;
if ( config.contains( "fsr_generator" ) ) {
cfg.fsrGenerator = config.at( "fsr_generator" ).get<FSRGenerator>();
}
// Set the type of threading to use, stdlib by default
cfg.threadModel = ThreadModel::StdLib;
if ( config.contains( "thread_model" ) ) {
cfg.threadModel = config.at( "thread_model" ).get<ThreadModel>();
}
// Set the RNG seed base (defaults to zero), to which the event number is added
cfg.rngSeed = 0;
if ( config.contains( "rng_seed" ) ) {
cfg.rngSeed = config.at( "rng_seed" ).get<std::size_t>();
}
// Set reference and output file names
// Insert fsrGenerator name if FSR simulation is not deactivated
const bool noFSR = std::find( cfg.extraCommands.begin(),
cfg.extraCommands.end(),
"noFSR" ) != cfg.extraCommands.end();
const std::string fileNameEnd =
noFSR ? ".root" : "_" + to_string( cfg.fsrGenerator ) + ".root";
const auto outFileStrSize = config.at( "outfile" ).get<std::string>().size() -
5;
cfg.outFileName =
config.at( "outfile" ).get<std::string>().substr( 0, outFileStrSize ) +
fileNameEnd;
cfg.refFileName = "Ref/" + cfg.outFileName;
cfg.decFileName = cfg.outFileName.substr( 0, cfg.outFileName.size() - 5 ) +
".dec";
cfg.debugFlag = ( config.contains( "debug_flag" ) &&
config.at( "debug_flag" ).is_boolean() )
? config.at( "debug_flag" ).get<bool>()
: false;
if ( config.contains( "do_conjugate_decay" ) &&
config.at( "do_conjugate_decay" ).is_array() ) {
cfg.doConjDecay =
config.at( "do_conjugate_decay" ).get<std::vector<bool>>();
}
if ( cfg.doConjDecay.size() != cfg.modelNames.size() ) {
cfg.doConjDecay.resize( cfg.modelNames.size(), false );
}
return cfg;
}
void TestDecayModel::run()
{
TestHistos theHistos;
const auto start{ std::chrono::steady_clock::now() };
if ( m_config.nThreads > 1 ) {
// Run multi-threaded using the specified thread model
switch ( m_config.threadModel ) {
case ThreadModel::StdLib:
theHistos = runStdThreads();
break;
case ThreadModel::TBB:
theHistos = runTBBThreads();
break;
}
} else {
// Run in the main thread
theHistos = runDecayBody( 0, m_config.nEvents );
}
const auto end{ std::chrono::steady_clock::now() };
const std::chrono::duration<double, std::milli> elapsed_ms{ ( end - start ) };
const std::chrono::duration<double, std::milli> elapsed_ms_per_event{
elapsed_ms / m_config.nEvents };
std::cout << "Took " << elapsed_ms.count() << " ms to generate "
<< m_config.nEvents << " events using " << m_config.nThreads
<< " " << to_string( m_config.threadModel ) << " threads ("
<< elapsed_ms_per_event.count() << " ms per event)" << std::endl;
// Normalise histograms.
theHistos.normalise();
// Compare with reference histograms
compareHistos( theHistos, m_config.refFileName );
// Create the root output file and write the histograms to it
// Only save the mixed/unmixed histogram for neutral K, D, B
std::unique_ptr<TFile> outFile{
TFile::Open( m_config.outFileName.c_str(), "recreate" ) };
theHistos.save( outFile.get() );
outFile->Close();
std::cout << "Created output file: " << m_config.outFileName.c_str()
<< std::endl;
}
TestHistos TestDecayModel::runStdThreads() const
{
// Determine the number of threads and the number of events per thread
const std::size_t nThreads{ std::min( m_config.nThreads, m_config.nEvents ) };
const std::size_t nEventsPerThread{ ( m_config.nEvents % nThreads )
? m_config.nEvents / nThreads + 1
: m_config.nEvents / nThreads };
// Create the store for the results from each thread
std::vector<std::future<TestHistos>> allHistos;
allHistos.reserve( nThreads );
// Launch the threads
std::size_t firstEvent{ 0 };
std::size_t nEvents{ nEventsPerThread };
for ( std::size_t iThread{ 0 }; iThread < nThreads; ++iThread ) {
// The last thread may need to generate slightly fewer events to get the required total
if ( ( firstEvent + nEvents ) >= m_config.nEvents ) {
nEvents = m_config.nEvents - firstEvent;
}
std::cout << "Thread " << iThread << " will generate " << nEvents
<< " events" << std::endl;
allHistos.emplace_back(
std::async( std::launch::async, [this, firstEvent, nEvents]() {
return runDecayBody( firstEvent, nEvents );
} ) );
firstEvent += nEvents;
}
// Now wait for each thread to finish
bool complete{ false };
do {
// Set the flag to completed and set it back if we find incomplete threads
complete = true;
for ( auto& future : allHistos ) {
auto status = future.wait_for( std::chrono::seconds( 10 ) );
if ( status != std::future_status::ready ) {
complete = false;
}
}
} while ( !complete );
// Accumulate the histograms from all threads
TestHistos theHistos{ allHistos[0].get() };
for ( std::size_t iThread{ 1 }; iThread < nThreads; ++iThread ) {
theHistos.add( allHistos[iThread].get() );
}
return theHistos;
}
TestHistos TestDecayModel::runTBBThreads() const
{
tbb::global_control gc{ tbb::global_control::parameter::max_allowed_parallelism,
m_config.nThreads };
TestHistos init;
return tbb::parallel_reduce(
tbb::blocked_range<std::size_t>( 0, m_config.nEvents ), init,
[this]( const tbb::blocked_range<std::size_t>& range,
const TestHistos& init ) -> TestHistos {
std::cout << "Thread "
<< tbb::this_task_arena::current_thread_index()
<< " will generate " << range.size()
<< " events: " << range.begin() << " - " << range.end()
<< std::endl;
TestHistos tmp{ init };
tmp.add( runDecayBody( range.begin(), range.size() ) );
return tmp;
},
[]( const TestHistos& lhs, const TestHistos& rhs ) -> TestHistos {
TestHistos tmp{ lhs };
tmp.add( rhs );
return tmp;
} );
}
TestHistos TestDecayModel::runDecayBody( const std::size_t firstEvent,
const std::size_t nEvents ) const
{
// 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
static thread_local auto randomEngine{ std::make_unique<EvtMTRandomEngine>() };
// TODO - need to streamline the extra models stuff
static thread_local EvtAbsRadCorr* radCorrEngine{ nullptr };
std::list<EvtDecayBase*> extraModels;
static thread_local bool initialised{ false };
#ifdef EVTGEN_EXTERNAL
if ( !initialised ) {
bool convertPythiaCodes( false );
bool useEvtGenRandom( true );
+ bool seedTauolaFortran( true );
EvtExternalGenList genList( convertPythiaCodes, "", "gamma",
- useEvtGenRandom );
+ useEvtGenRandom, seedTauolaFortran );
switch ( m_config.fsrGenerator ) {
case FSRGenerator::PHOTOS:
radCorrEngine = genList.getPhotosModel();
break;
case FSRGenerator::SherpaPhotons1:
radCorrEngine = genList.getSherpaPhotonsModel( 1e-7, 1, 0 );
break;
case FSRGenerator::SherpaPhotons20:
radCorrEngine = genList.getSherpaPhotonsModel( 1e-7, 2, 0 );
break;
case FSRGenerator::SherpaPhotons21:
radCorrEngine = genList.getSherpaPhotonsModel( 1e-7, 2, 1 );
break;
}
extraModels = genList.getListOfModels();
}
#endif
static thread_local EvtGen theGen{ "../DECAY.DEC", "../evt.pdl",
randomEngine.get(), radCorrEngine,
&extraModels };
// Creates a decay file based on json file input
// We don't want this to be called by every thread!
std::call_once( m_createDecFile_threadlock, [this]() { createDecFile(); } );
// Read the decay file
if ( !initialised ) {
theGen.readUDecay( m_config.decFileName.c_str() );
initialised = true;
}
// Define the histograms to be saved
TestHistos theHistos{ m_config.testHistograms };
// Generate events and fill histograms
generateEvents( theGen, theHistos, firstEvent, nEvents );
return theHistos;
}
void TestDecayModel::createDecFile() const
{
// Create (or overwrite) the decay file
std::ofstream decFile{ m_config.decFileName };
// Create daughter aliases if needed
std::vector<std::string> aliasPrefix;
for ( std::size_t daughter_index{ 0 };
daughter_index < m_config.daughterNames.size(); daughter_index++ ) {
if ( !m_config.grandDaughterNames.empty() &&
!m_config.grandDaughterNames[daughter_index].empty() ) {
decFile << "Alias My" << m_config.daughterNames[daughter_index] << " "
<< m_config.daughterNames[daughter_index] << std::endl;
if ( m_config.doConjDecay[daughter_index + 1] ) {
const EvtId daugID{
EvtPDL::getId( m_config.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( m_config.daughterNames ),
std::end( m_config.daughterNames ),
daugConjID.getName() ) !=
std::end( m_config.daughterNames ) ) {
conjName_Alias = conjName_Alias + "_" + daughter_index;
}
decFile << "Alias My" << conjName_Alias << " " << conjName
<< std::endl;
decFile << "ChargeConj My"
<< m_config.daughterNames[daughter_index] << " My"
<< conjName_Alias << std::endl;
} else if ( m_config.doConjDecay[0] ) {
decFile << "ChargeConj My"
<< m_config.daughterNames[daughter_index] << " My"
<< m_config.daughterNames[daughter_index] << std::endl;
}
aliasPrefix.push_back( "My" );
} else {
aliasPrefix.push_back( "" );
}
}
for ( const auto& iExtra : m_config.extraCommands ) {
decFile << iExtra << std::endl;
}
// Parent decay
decFile << "Decay " << m_config.parentName << std::endl;
decFile << "1.0";
for ( std::size_t daughter_index{ 0 };
daughter_index < m_config.daughterNames.size(); daughter_index++ ) {
decFile << " " << aliasPrefix[daughter_index]
<< m_config.daughterNames[daughter_index];
}
decFile << " " << m_config.modelNames[0];
for ( const auto& par : m_config.modelParameters[0] ) {
decFile << " " << par;
}
decFile << ";" << std::endl;
decFile << "Enddecay" << std::endl;
if ( m_config.doConjDecay[0] ) {
EvtId parID{ EvtPDL::getId( m_config.parentName ) };
EvtId parConjID{ EvtPDL::chargeConj( parID ) };
decFile << "CDecay " << parConjID.getName() << std::endl;
}
// Daughter decays into granddaughters
for ( std::size_t daughter_index{ 0 };
daughter_index < m_config.grandDaughterNames.size(); daughter_index++ ) {
if ( m_config.grandDaughterNames[daughter_index].empty() )
continue;
decFile << "Decay " << aliasPrefix[daughter_index]
<< m_config.daughterNames[daughter_index] << std::endl;
decFile << "1.0";
for ( std::size_t grandDaughter_index{ 0 };
grandDaughter_index <
m_config.grandDaughterNames[daughter_index].size();
grandDaughter_index++ ) {
decFile << " "
<< m_config.grandDaughterNames[daughter_index][grandDaughter_index];
}
decFile << " " << m_config.modelNames[daughter_index + 1];
for ( const auto& par : m_config.modelParameters[daughter_index + 1] ) {
decFile << " " << par;
}
decFile << ";" << std::endl;
decFile << "Enddecay" << std::endl;
if ( m_config.doConjDecay[daughter_index + 1] ) {
EvtId daugID{
EvtPDL::getId( m_config.daughterNames[daughter_index] ) };
EvtId daugConjID{ EvtPDL::chargeConj( daugID ) };
std::string conjName_Alias{ daugConjID.getName() };
if ( std::find( std::begin( m_config.daughterNames ),
std::end( m_config.daughterNames ),
daugConjID.getName() ) !=
std::end( m_config.daughterNames ) ) {
conjName_Alias = conjName_Alias + "_" + daughter_index;
}
decFile << "CDecay " << aliasPrefix[daughter_index]
<< conjName_Alias << std::endl;
}
}
decFile << "End" << std::endl;
decFile.close();
}
void TestDecayModel::generateEvents( EvtGen& theGen, TestHistos& theHistos,
const std::size_t firstEvent,
const std::size_t nEvents ) const
{
// Generate the decays
const EvtId parId{ EvtPDL::getId( m_config.parentName.c_str() ) };
const EvtId conjId{ m_config.doConjDecay[0] ? EvtPDL::chargeConj( parId )
: parId };
for ( std::size_t i{ firstEvent }; i < ( firstEvent + nEvents ); i++ ) {
if ( i % 1000 == 0 ) {
std::cout << "Event " << firstEvent + nEvents - i << std::endl;
}
// seed the RNG based on the event number
EvtRandom::setSeed( m_config.rngSeed + i );
// 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 };
TH1* mixedHist{ theHistos.getMixedHist() };
if ( mixedHist ) {
if ( parent->getNDaug() == 1 ) {
prodParent = parent;
parent = prodParent->getDaug( 0 );
mixedHist->Fill( 1 );
} else {
mixedHist->Fill( 0 );
}
}
// To debug
if ( m_config.debugFlag ) {
std::cout << "Parent PDG code: " << parent->getPDGId()
<< " has daughters " << parent->getNDaug() << std::endl;
for ( std::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 ( std::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";
const std::string perDaughter = "_perDaughter";
// Store information
for ( auto& [info, hist] : theHistos.get1DHistos() ) {
if ( !hist ) {
continue;
}
const std::string varName = info.getName();
std::string reducedVarName = varName;
const std::string::size_type findFSRstr = varName.find( fsrStr );
const std::string::size_type findPerDaugStr = varName.find(
perDaughter );
// If the variable name does not have the substrings '_perDaughter' and '_FSRPhotons', then add an entry per event and continue
if ( findFSRstr == std::string::npos &&
findPerDaugStr == std::string::npos ) {
const double value{
getValue( parent, varName, info.getd1(), info.getd2() ) };
hist->Fill( value );
continue;
}
// If the variable name has the substring '_perDaughter' and does not have the substring '_FSRPhotons', then add an entry per daughter and continue
else if ( findFSRstr == std::string::npos ) {
reducedVarName.erase( findPerDaugStr, perDaughter.length() );
// Figure out whether the variable indicates to be saved for a requested particle type
const std::string::size_type findIsStr = reducedVarName.find(
"_is" );
const std::string requestedType = findIsStr == std::string::npos
? ""
: reducedVarName.substr(
findIsStr + 3 );
if ( !requestedType.empty() ) {
reducedVarName.erase( findIsStr,
reducedVarName.length() - findIsStr );
}
for ( int iDaug{ 0 }; iDaug < (int)parent->getNDaug(); iDaug++ ) {
const double value{
getValue( parent, reducedVarName, iDaug + 1, 0 ) };
const int partGroup = getPartGroup(
parent->getDaug( iDaug )->getPDGId() );
if ( requestedType.empty() ||
isPartType( partGroup, requestedType ) ) {
hist->Fill( value );
}
}
continue;
}
// If otherwise the variable contains the substring '_FSRPhotons', then add an entry per photon
if ( leadingChargedDaughter == -1 ) {
leadingChargedDaughter = findChargedDaughterWithMaxE( parent );
}
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] : theHistos.get2DHistos() ) {
if ( !hist ) {
continue;
}
const double valueX{ getValue( parent,
info.getName( HistInfo::Axis::X ),
info.getd1( HistInfo::Axis::X ),
info.getd2( HistInfo::Axis::X ) ) };
const double valueY{ getValue( parent,
info.getName( HistInfo::Axis::Y ),
info.getd1( HistInfo::Axis::Y ),
info.getd2( HistInfo::Axis::Y ) ) };
hist->Fill( valueX, valueY );
}
if ( m_config.debugFlag ) {
if ( prodParent ) {
prodParent->printTree();
} else {
parent->printTree();
}
}
// Cleanup
if ( prodParent ) {
prodParent->deleteTree();
} else {
parent->deleteTree();
}
}
}
int TestDecayModel::findChargedDaughterWithMaxE( 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( "particleType" ) ) {
if ( par1 ) {
value = getPartGroup( par1->getPDGId() );
} else {
value = getPartGroup( selectedParent->getPDGId() );
}
} else if ( !selectedVarName.compare( "nDaug" ) ) {
// Number of daughters
value = sel_NDaugMax;
} 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" ) ||
!selectedVarName.compare( "cosHelParent" ) ) {
// Cosine of helicity angle
EvtVector4R p12;
EvtVector4R p1Res;
if ( !selectedVarName.compare( "cosHelParent" ) ) {
p12 = selectedParent->getP4Lab();
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 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 p4_Res{ res->getP4Lab() };
const EvtVector4R p4_p4{ bac != nullptr ? bac->getP4Lab()
: EvtVector4R() };
const EvtVector4R p4_p1{ daug1 != nullptr ? daug1->getP4Lab()
: EvtVector4R() };
const EvtVector4R p4_p3{ daug3 != nullptr ? daug3->getP4Lab()
: EvtVector4R() };
// Boost 4-vector for resonance frame
const EvtVector4R boost{ p4_Res.get( 0 ), -p4_Res.get( 1 ),
-p4_Res.get( 2 ), -p4_Res.get( 3 ) };
// Momentum of p1 and p3 in resonance frame
const EvtVector4R p1Res{ boostTo( p4_p1, boost ) };
const EvtVector4R p3Res{ boostTo( p4_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_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 p4_p1{ daug1 != nullptr ? daug1->getP4Lab()
: EvtVector4R() };
const EvtVector4R p4_p2{ daug2 != nullptr ? daug2->getP4Lab()
: EvtVector4R() };
const EvtVector4R p4_p3{ daug3 != nullptr ? daug3->getP4Lab()
: EvtVector4R() };
// p1 + p2
const EvtVector4R p12{ p4_p1 + p4_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( p4_p1, boost ) };
const EvtVector4R p3_12{ boostTo( p4_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 ( std::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 ( std::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 TestHistos& theHistos,
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?
// TODO - add comparison of mixedHist
for ( auto& [_, hist] : theHistos.get1DHistos() ) {
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.get(), 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] : theHistos.get2DHistos() ) {
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.get(), 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();
}
int TestDecayModel::getPartGroup( const int PDGId ) const
{
int group( -1 );
const int absPDGId = std::abs( PDGId );
if ( absPDGId >= 11 && absPDGId <= 16 ) {
group = 0; // leptons
} else if ( absPDGId == 22 ) {
group = 1; // photon
} else if ( absPDGId == 211 ) {
group = 2; // pi+-
} else if ( absPDGId == 111 ) {
group = 3; // pi0
} else if ( absPDGId == 321 ) {
group = 4; // K+-
} else if ( absPDGId == 311 || absPDGId == 130 || absPDGId == 310 ) {
group = 5; // K0
} else if ( absPDGId == 411 ) {
group = 6; // D+-
} else if ( absPDGId == 421 ) {
group = 7; // D0
} else if ( absPDGId == 2212 || absPDGId == 2112 || absPDGId == 2224 ||
absPDGId == 2214 || absPDGId == 2114 || absPDGId == 1114 ) {
group = 8; // light baryons
} else if ( absPDGId >= 3112 && absPDGId <= 3334 ) {
group = 9; // strange baryons
} else if ( absPDGId != 0 ) {
group = 10; // other particles
}
return group;
}
bool TestDecayModel::isPartType( const int group,
const std::string& particleType ) const
{
if ( ( !particleType.compare( "LeptonOrPhoton" ) &&
( group == 0 || group == 1 ) ) ||
( !particleType.compare( "Pion" ) && ( group == 2 || group == 3 ) ) ||
( !particleType.compare( "Kaon" ) && ( group == 4 || group == 5 ) ) ||
( !particleType.compare( "Dmeson" ) && ( group == 6 || group == 7 ) ) ||
( !particleType.compare( "Baryon" ) && ( group == 8 || group == 9 ) ) ||
( !particleType.compare( "Other" ) && group == 10 ) ) {
return true;
} else {
return false;
}
}
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 || argc > 3 ) {
std::cerr << "Expecting at least one argument: test configuration json file"
<< "\nAdditional argument supported for: general configuration json file"
<< std::endl;
return 1;
}
// Tweak ROOT behaviour
ROOT::EnableThreadSafety();
TH1::AddDirectory( kFALSE );
// Load input file in json format
const std::string testConfigFileName{ argv[1] };
const std::string generalConfigFileName{
( argc > 2 ) ? argv[2] : "jsonFiles/config/default.json" };
json generalConfig;
json testConfig;
{
std::ifstream inputStr{ generalConfigFileName };
inputStr >> generalConfig;
}
{
std::ifstream inputStr{ testConfigFileName };
inputStr >> testConfig;
}
if ( testConfig.is_array() ) {
for ( auto& cc : testConfig ) {
cc.merge_patch( generalConfig );
TestDecayModel test{ cc };
test.run();
}
} else {
testConfig.merge_patch( generalConfig );
TestDecayModel test{ testConfig };
test.run();
}
return 0;
}

File Metadata

Mime Type
text/x-diff
Expires
Sun, Feb 23, 2:01 PM (2 h, 17 m)
Storage Engine
blob
Storage Format
Raw Data
Storage Handle
4486435
Default Alt Text
(134 KB)

Event Timeline