diff --git a/EvtGenBase/EvtParticle.hh b/EvtGenBase/EvtParticle.hh
index 953613a..084cb9f 100644
--- a/EvtGenBase/EvtParticle.hh
+++ b/EvtGenBase/EvtParticle.hh
@@ -1,502 +1,503 @@
/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
* *
* This file is part of EvtGen. *
* *
* EvtGen is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* EvtGen is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/
#ifndef EVTPARTICLE_HH
#define EVTPARTICLE_HH
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtSpinDensity.hh"
#include "EvtGenBase/EvtSpinType.hh"
#include "EvtGenBase/EvtVector4R.hh"
#include <assert.h>
#include <map>
#include <string>
#include <vector>
class EvtDiracSpinor;
class EvtVector4C;
class EvtTensor4C;
class EvtStdHep;
class EvtSecondary;
class EvtRaritaSchwinger;
const int MAX_DAUG = 100;
const int MAX_LEVEL = 10;
const int MAX_TRIES = 10000;
class EvtParticle {
public:
/**
* Default constructor.
*/
EvtParticle();
/**
* Destructor.
*/
virtual ~EvtParticle();
/**
* Returns polarization vector in the parents restframe.
*/
virtual EvtVector4C epsParent( int i ) const;
/**
* Returns polarization vector in the particles own restframe.
*/
virtual EvtVector4C eps( int i ) const;
/**
* Returns polarization vector in the parents restframe for a photon.
*/
virtual EvtVector4C epsParentPhoton( int i );
/**
* Returns polarization vector in the particles own restframe for a photon.
*/
virtual EvtVector4C epsPhoton( int i );
/**
* Returns Dirac spinor in the parents restframe for a Dirac particle.
*/
virtual EvtDiracSpinor spParent( int ) const;
/**
* Returns Dirac spinor in the particles own restframe for a Dirac particle.
*/
virtual EvtDiracSpinor sp( int ) const;
/**
* Returns Dirac spinor in the parents restframe for a Neutrino particle.
*/
virtual EvtDiracSpinor spParentNeutrino() const;
/**
* Returns Dirac spinor in the particles own restframe for a
* Neutrino particle.
*/
virtual EvtDiracSpinor spNeutrino() const;
/**
* Returns tensor in the parents restframe for a spin 2 particle.
*/
virtual EvtTensor4C epsTensorParent( int i ) const;
/**
* Returns tensor in the particles own restframe for a spin 2 particle.
*/
virtual EvtTensor4C epsTensor( int i ) const;
/**
* Returns Rarita-Schwinger spinor in the parents restframe for a
* Rarita-Schwinger particle.
*/
virtual EvtRaritaSchwinger spRSParent( int ) const;
/**
* Returns Rarita-Schwinger spinor in the particles own restframe for a
* Rarita-Schwinger particle.
*/
virtual EvtRaritaSchwinger spRS( int ) const;
/**
* Initialiaze particle with id and 4momentum.
*/
virtual void init( EvtId part_n, const EvtVector4R& p4 ) = 0;
/**
* Add another daughter to the particle
*/
void addDaug( EvtParticle* node );
/**
* Decay particle
*/
void decay();
/**
* Delete a decay chain
*/
void deleteTree();
void deleteDaughters( bool keepChannel = false );
/**
* Should only be used internally.
*/
void setChannel( int i );
/**
* Creates the daughters in the list of ids and
* adds them to the parent. Note that momentum
* is left uninitialized, this is _only_ creation.
*/
void makeDaughters( unsigned int ndaug, EvtId* id );
/**
* Creates the daughters in the list of ids and
* adds them to the parent. Note that momentum
* is left uninitialized, this is _only_ creation.
*/
void makeDaughters( unsigned int ndaug, std::vector<EvtId> idVector );
/**
* Similar to the routine above except that here
* momentum is generated according to phase space
* daughters are filled with this momentum.
*/
double initializePhaseSpace( unsigned int numdaughter, EvtId* daughters,
bool forceResetMasses = false,
double poleSize = -1., int whichTwo1 = 0,
int whichTwo2 = 1 );
/**
* Get pointer the the i:th daugther.
*/
EvtParticle* getDaug( const int i ) { return _daug[i]; }
/**
* Get const pointer the the i:th daugther.
*/
const EvtParticle* getDaug( const int i ) const { return _daug[i]; }
/**
* Iterates over the particles in a decay chain.
*/
EvtParticle* nextIter( EvtParticle* rootOfTree = 0 );
/**
* Makes stdhep list
*/
void makeStdHep( EvtStdHep& stdhep, EvtSecondary& secondary,
EvtId* stable_parent_ihep );
void makeStdHep( EvtStdHep& stdhep );
/**
* Gets 4vector in the labframe, i.e., the frame in which the root
* particles momentum is measured.
*/
EvtVector4R getP4Lab() const;
/**
* Gets 4vector in the labframe for the 4-momentum before FSR was
* generated in the parents decay. The lab frame is where the root
* particles momentum is measured.
*/
EvtVector4R getP4LabBeforeFSR();
/**
* Gets 4vector in the particles restframe, i.e. this functiont will
* return (m,0,0,0)
*/
EvtVector4R getP4Restframe() const;
/**
* Returns the 4position of the particle in the lab frame.
*/
EvtVector4R get4Pos() const;
/**
* Returns pointer to parent particle.
*/
EvtParticle* getParent() const;
/**
* Makes partptr the idaug:th daugther.
*/
void insertDaugPtr( int idaug, EvtParticle* partptr )
{
_daug[idaug] = partptr;
partptr->_parent = this;
}
/**
* Returns mass of particle.
*/
double mass() const;
/**
* Used internally to decide if first time particle is decayed.
*/
int firstornot() const;
void setFirstOrNot();
void resetFirstOrNot();
/**
* Returns Id of particle.
*/
EvtId getId() const;
/**
* Returns the PDG id of the particle
*/
int getPDGId() const;
/**
* Returns particle type.
*/
EvtSpinType::spintype getSpinType() const;
/**
* Returns number of spin states of the particle.
*/
int getSpinStates() const;
/**
* Returns 4momentum in parents restframe.
*/
const EvtVector4R& getP4() const;
/**
* Sets the 4momentum in the parents restframe.
*/
void setP4( const EvtVector4R& p4 )
{
_p = p4;
_pBeforeFSR = p4;
}
void setP4WithFSR( const EvtVector4R& p4 ) { _p = p4; }
void setFSRP4toZero() { _pBeforeFSR.set( 0.0, 0.0, 0.0, 0.0 ); }
/**
* Retunrs the decay channel.
*/
int getChannel() const;
/**
* Returns number of daugthers.
*/
size_t getNDaug() const;
void resetNDaug()
{
_ndaug = 0;
return;
}
/**
* Prints out the particle "tree" of a given particle. The
* tree consists of all daughters (and their daughters, etc)
* and their properties.
*/
void printTree() const;
void printTreeRec( unsigned int level ) const;
std::string treeStr() const;
std::string treeStrRec( unsigned int level ) const;
/**
* Prints information for the particle.
*/
void printParticle() const;
/**
* Set lifetime of the particle in parents restframe.
*/
void setLifetime( double tau );
/**
* Generate lifetime according to pure exponential.
*/
void setLifetime();
/**
* Returns the lifetime.
*/
double getLifetime() const;
/**
* Set diagonal spindensity matrix.
*/
void setDiagonalSpinDensity();
/**
* Set spindensity matrix for e+e- -> V
*/
void setVectorSpinDensity();
/**
* Set forward spin density matrix.
*/
void setSpinDensityForward( const EvtSpinDensity& rho )
{
_rhoForward = rho;
}
/**
* Set forward spin density matrix according to the density matrix
* rho in the helicity amplitude basis.
*/
void setSpinDensityForwardHelicityBasis( const EvtSpinDensity& rho );
void setSpinDensityForwardHelicityBasis( const EvtSpinDensity& rho,
double alpha, double beta,
double gamma );
/**
* Returns a rotation matrix need to rotate the basis state
* to the helicity basis. The EvtSpinDensity matrix is just use
* as a matrix here. This function is to be implemented in each
* derived class.
*/
virtual EvtSpinDensity rotateToHelicityBasis() const = 0;
virtual EvtSpinDensity rotateToHelicityBasis( double alpha, double beta,
double gamma ) const = 0;
/**
* Get forward spin density matrix.
*/
EvtSpinDensity getSpinDensityForward() { return _rhoForward; }
/**
* Set backward spin density matrix.
*/
void setSpinDensityBackward( const EvtSpinDensity& rho )
{
_rhoBackward = rho;
}
/**
* Get backward spin density matrix.
*/
EvtSpinDensity getSpinDensityBackward() { return _rhoBackward; }
//Hacks will be removed when better solutions are thought of!
//This is used to suppress use of random numbers when doing initialization
//of some models.
void noLifeTime() { _genlifetime = 0; }
//lange - April 29, 2002
void setId( EvtId id ) { _id = id; }
void initDecay( bool useMinMass = false );
bool generateMassTree();
double compMassProb();
//setMass will blow away any existing 4vector
void setMass( double m ) { _p = EvtVector4R( m, 0.0, 0.0, 0.0 ); }
//void setMixed() {_mix=true;}
//void setUnMixed() {_mix=false;}
//bool getMixed() {return _mix;}
//void takeCConj() {EvtGenReport(EVTGEN_INFO,"EvtGen") << "should take conj\n";}
//this means that the particle has gone through initDecay
// and thus has a mass
bool isInitialized() { return _isInit; }
bool hasValidP4() { return _validP4; }
bool isDecayed() { return _isDecayed; }
// decay prob - only relevent if already decayed
// and is a scalar particle
// returned is a double* that should be prob/probMax
- double* decayProb() { return _decayProb; }
+ // FIXME - this should probably be changed to std::optional
+ const double* decayProb() const { return _decayProb; }
void setDecayProb( double p );
// Return the name of the particle (from the EvtId number)
std::string getName() const;
// Specify whether the particle has a special named attribute with
// a set value. By default, nothing is set, but derived classes
// can set this to mean something specific, e.g. if a photon is FSR
void setAttribute( std::string attName, int attValue )
{
_intAttributes[attName] = attValue;
}
// Retrieve the integer value for the given attribute name
int getAttribute( std::string attName );
// Specify if the particle has a double attribute value, e.g. amplitude weight.
// By default, nothing is set, but derived classes can set this to mean something specific
void setAttributeDouble( std::string attName, double attValue )
{
_dblAttributes[attName] = attValue;
}
// Retrieve the double value for the given attribute name
double getAttributeDouble( std::string attName );
protected:
void setp( double e, double px, double py, double pz )
{
_p.set( e, px, py, pz );
_pBeforeFSR = _p;
}
void setp( const EvtVector4R& p4 )
{
_p = p4;
_pBeforeFSR = _p;
}
void setpart_num( EvtId particle_number )
{
assert( _channel == -10 || _id.getId() == particle_number.getId() ||
_id.getId() == -1 );
_id = particle_number;
}
bool _validP4;
// A typedef to define the attribute (name, integer) map
typedef std::map<std::string, int> EvtAttIntMap;
EvtAttIntMap _intAttributes;
// A typedef to define the attribute (name, double) map
typedef std::map<std::string, double> EvtAttDblMap;
EvtAttDblMap _dblAttributes;
private:
EvtParticle* _daug[MAX_DAUG];
size_t _ndaug;
EvtParticle* _parent;
int _channel;
int _first;
EvtId _id;
EvtVector4R _p;
EvtVector4R _pBeforeFSR;
double _t;
bool _isInit;
bool _isDecayed;
//bool _mix;
EvtSpinDensity _rhoForward;
EvtSpinDensity _rhoBackward;
void makeStdHepRec( int firstparent, int lastparent, EvtStdHep& stdhep,
EvtSecondary& secondary, EvtId* stable_parent_ihep );
void makeStdHepRec( int firstparent, int lastparent, EvtStdHep& stdhep );
//This is a hack until things gets straightened out. (Ryd)
int _genlifetime;
//should never be used, therefor is private.
//these does _not_ have an implementation
EvtParticle& operator=( const EvtParticle& p );
EvtParticle( const EvtParticle& p );
double* _decayProb;
};
#endif
diff --git a/test/testDecayModel.cc b/test/testDecayModel.cc
index bdea396..92e6cb1 100644
--- a/test/testDecayModel.cc
+++ b/test/testDecayModel.cc
@@ -1,1215 +1,1215 @@
/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
* *
* This file is part of EvtGen. *
* *
* EvtGen is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* EvtGen is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/
#include "testDecayModel.hh"
#include "EvtGen/EvtGen.hh"
#include "EvtGenBase/EvtAbsRadCorr.hh"
#include "EvtGenBase/EvtConst.hh"
#include "EvtGenBase/EvtDecayBase.hh"
#include "EvtGenBase/EvtId.hh"
#include "EvtGenBase/EvtKine.hh"
#include "EvtGenBase/EvtMTRandomEngine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtParticleFactory.hh"
#include "EvtGenBase/EvtRandom.hh"
#include "EvtGenBase/EvtVector4R.hh"
#ifdef EVTGEN_EXTERNAL
#include "EvtGenExternal/EvtExternalGenList.hh"
#endif
#include <fstream>
#include <iostream>
#include <list>
#include <memory>
using nlohmann::json;
TestDecayModel::TestDecayModel( const json& config ) : m_config{ config }
{
}
bool TestDecayModel::checkMandatoryFields()
{
const std::array<std::string, 8> mandatoryFields{
"parent", "daughters", "models", "parameters",
"outfile", "events", "reference", "histograms" };
const std::array<std::string, 7> mandatoryHistoFields{
"title", "variable", "d1", "d2", "nbins", "xmin", "xmax" };
const std::array<std::string, 6> extra2DHistoFields{ "variableY", "d1Y",
"d2Y", "nbinsY",
"ymin", "ymax" };
bool allMandatoryFields{ true };
for ( const auto& field : mandatoryFields ) {
if ( !m_config.contains( field ) ) {
std::cerr << "ERROR : json does not contain required field: " << field
<< std::endl;
allMandatoryFields = false;
continue;
}
if ( field == "histograms" ) {
- json jHistos = m_config["histograms"];
- for ( auto hInfo : jHistos ) {
+ const json& jHistos{ m_config.at( "histograms" ) };
+ for ( const auto& hInfo : jHistos ) {
for ( const auto& hField : mandatoryHistoFields ) {
if ( !hInfo.contains( hField ) ) {
std::cerr
<< "ERROR : json does not contain required field for histogram definition: "
<< hField << std::endl;
allMandatoryFields = false;
}
}
if ( hInfo.contains( extra2DHistoFields[0] ) ) {
for ( const auto& hField : extra2DHistoFields ) {
if ( !hInfo.contains( hField ) ) {
std::cerr
<< "ERROR : json does not contain required field for 2D histogram definition: "
<< hField << std::endl;
allMandatoryFields = false;
}
}
}
}
}
}
return allMandatoryFields;
}
bool TestDecayModel::run()
{
// Check that we have, and then get all the mandatory fields first
if ( !checkMandatoryFields() ) {
std::cerr << "ERROR : json does not contain all mandatory fields - dumping config for debugging:\n";
std::cerr << m_config << std::endl;
return false;
}
- const auto parentName = m_config["parent"].get<std::string>();
- const auto daughterNames =
- m_config["daughters"].get<std::vector<std::string>>();
- const auto modelNames = m_config["models"].get<std::vector<std::string>>();
- const auto modelParameters =
- m_config["parameters"].get<std::vector<std::vector<std::string>>>();
- const auto outFileName = m_config["outfile"].get<std::string>();
- const auto nEvents = m_config["events"].get<int>();
- const auto refFileName = m_config["reference"].get<std::string>();
+ const auto parentName{ m_config.at( "parent" ).get<std::string>() };
+ const auto daughterNames{
+ m_config.at( "daughters" ).get<std::vector<std::string>>() };
+ const auto modelNames{
+ m_config.at( "models" ).get<std::vector<std::string>>() };
+ const auto modelParameters{
+ m_config.at( "parameters" ).get<std::vector<std::vector<std::string>>>() };
+ const auto outFileName{ m_config.at( "outfile" ).get<std::string>() };
+ const auto nEvents{ m_config.at( "events" ).get<int>() };
+ const auto refFileName{ m_config.at( "reference" ).get<std::string>() };
// Then check for optional fields, setting default values if not present
- json grandDaughters = m_config["grand_daughters"];
- const auto grandDaughterNames =
- grandDaughters.is_array()
- ? grandDaughters.get<std::vector<std::vector<std::string>>>()
- : std::vector<std::vector<std::string>>{};
+ const auto grandDaughterNames{
+ ( m_config.contains( "grand_daughters" ) &&
+ m_config.at( "grand_daughters" ).is_array() )
+ ? m_config.at( "grand_daughters" )
+ .get<std::vector<std::vector<std::string>>>()
+ : std::vector<std::vector<std::string>>{} };
- json extras = m_config["extras"];
- const auto extraCommands = extras.is_array()
- ? extras.get<std::vector<std::string>>()
- : std::vector<std::string>{};
+ const auto extraCommands{
+ ( m_config.contains( "extras" ) && m_config.at( "extras" ).is_array() )
+ ? m_config.at( "extras" ).get<std::vector<std::string>>()
+ : std::vector<std::string>{} };
- json debug = m_config["debug_flag"];
- const auto debugFlag = debug.is_boolean() ? debug.get<bool>() : false;
+ const auto debugFlag{ ( m_config.contains( "debug_flag" ) &&
+ m_config.at( "debug_flag" ).is_boolean() )
+ ? m_config.at( "debug_flag" ).get<bool>()
+ : false };
- json conjugates = m_config["do_conjugate_decay"];
std::vector<bool> doConjDecay;
- if ( conjugates.is_array() ) {
- doConjDecay = conjugates.get<std::vector<bool>>();
+ if ( m_config.contains( "do_conjugate_decay" ) &&
+ m_config.at( "do_conjugate_decay" ).is_array() ) {
+ doConjDecay = m_config.at( "do_conjugate_decay" ).get<std::vector<bool>>();
}
if ( doConjDecay.size() != modelNames.size() ) {
doConjDecay.resize( modelNames.size(), false );
}
// Initialise the EvtGen object and hence the EvtPDL tables
// The EvtGen object is used by generateEvents, while the
// latter are also used within createDecFile
// Define the random number generator
- auto randomEngine = std::make_unique<EvtMTRandomEngine>();
+ auto randomEngine{ std::make_unique<EvtMTRandomEngine>() };
- EvtAbsRadCorr* radCorrEngine = nullptr;
+ EvtAbsRadCorr* radCorrEngine{ nullptr };
std::list<EvtDecayBase*> extraModels;
#ifdef EVTGEN_EXTERNAL
bool convertPythiaCodes( false );
bool useEvtGenRandom( true );
EvtExternalGenList genList( convertPythiaCodes, "", "gamma", useEvtGenRandom );
radCorrEngine = genList.getPhotosModel();
extraModels = genList.getListOfModels();
#endif
EvtGen theGen( "../DECAY.DEC", "../evt.pdl", randomEngine.get(),
radCorrEngine, &extraModels );
/*! Creates a decay file based on json file input. */
- const std::string decFileName = outFileName.substr( 0,
- outFileName.size() - 5 );
- const std::string decFile = createDecFile( parentName, daughterNames,
- grandDaughterNames, modelNames,
- modelParameters, doConjDecay,
- extraCommands, decFileName );
+ const std::string decFileName{
+ outFileName.substr( 0, outFileName.size() - 5 ) };
+ const std::string decFile{ createDecFile( parentName, daughterNames,
+ grandDaughterNames, modelNames,
+ modelParameters, doConjDecay,
+ extraCommands, decFileName ) };
/*! Define the root output file and histograms to be saved. */
- TFile* outFile = TFile::Open( outFileName.c_str(), "recreate" );
- defineHistos( outFile );
+ std::unique_ptr<TFile> outFile{
+ TFile::Open( outFileName.c_str(), "recreate" ) };
+ defineHistos( outFile.get() );
/*! Generate events and fill histograms. */
generateEvents( theGen, decFile, parentName, doConjDecay[0], nEvents,
debugFlag );
// Normalize histograms.
- for ( auto hist : m_1DhistVect ) {
- double area = hist.second->Integral();
+ for ( auto& [_, hist] : m_1DhistVect ) {
+ const double area{ hist->Integral() };
if ( area > 0.0 ) {
- hist.second->Scale( 1.0 / area );
+ hist->Scale( 1.0 / area );
}
}
- for ( auto hist : m_2DhistVect ) {
- double area = hist.second->Integral();
+ for ( auto& [_, hist] : m_2DhistVect ) {
+ const double area{ hist->Integral() };
if ( area > 0.0 ) {
- hist.second->Scale( 1.0 / area );
+ hist->Scale( 1.0 / area );
}
}
/*! Compare with reference histograms. */
compareHistos( refFileName );
// Write output.
outFile->cd();
- for ( auto& hist : m_1DhistVect ) {
- hist.second->Write();
+ for ( auto& [_, hist] : m_1DhistVect ) {
+ hist->Write();
}
- for ( auto& hist : m_2DhistVect ) {
- hist.second->Write();
+ for ( auto& [_, hist] : m_2DhistVect ) {
+ hist->Write();
}
if ( m_mixedHist ) {
m_mixedHist->Write();
}
outFile->Close();
std::cout << "Created output file: " << outFileName.c_str() << std::endl;
return true;
}
std::string TestDecayModel::createDecFile(
const std::string& parent, const std::vector<std::string>& daughterNames,
const std::vector<std::vector<std::string>>& grandDaughterNames,
const std::vector<std::string>& modelNames,
const std::vector<std::vector<std::string>>& parameters,
const std::vector<bool>& doConjDecay,
const std::vector<std::string>& extras, const std::string decFileName ) const
{
// Create (or overwrite) the decay file
std::string decName( decFileName + ".dec" );
std::ofstream decFile;
decFile.open( decName.c_str() );
// Create daughter aliases if needed
std::vector<std::string> aliasPrefix;
- for ( long unsigned int daughter_index = 0;
+ for ( long unsigned int daughter_index{ 0 };
daughter_index < daughterNames.size(); daughter_index++ ) {
if ( !grandDaughterNames.empty() &&
!grandDaughterNames[daughter_index].empty() ) {
decFile << "Alias My" << daughterNames[daughter_index] << " "
<< daughterNames[daughter_index] << std::endl;
if ( doConjDecay[daughter_index + 1] ) {
const EvtId daugID{
EvtPDL::getId( daughterNames[daughter_index] ) };
const EvtId daugConjID{ EvtPDL::chargeConj( daugID ) };
const std::string conjName{ daugConjID.getName() };
std::string conjName_Alias{ daugConjID.getName() };
if ( std::find( std::begin( daughterNames ),
std::end( daughterNames ), daugConjID.getName() ) !=
std::end( daughterNames ) ) {
conjName_Alias = conjName_Alias + "_" + daughter_index;
}
decFile << "Alias My" << conjName_Alias << " " << conjName
<< std::endl;
decFile << "ChargeConj My" << daughterNames[daughter_index]
<< " My" << conjName_Alias << std::endl;
} else if ( doConjDecay[0] ) {
decFile << "ChargeConj My" << daughterNames[daughter_index]
<< " My" << daughterNames[daughter_index] << std::endl;
}
aliasPrefix.push_back( "My" );
} else {
aliasPrefix.push_back( "" );
}
}
- for ( auto iExtra : extras ) {
+ for ( const auto& iExtra : extras ) {
decFile << iExtra << std::endl;
}
// Parent decay
decFile << "Decay " << parent << std::endl;
decFile << "1.0";
- for ( long unsigned int daughter_index = 0;
+ for ( long unsigned int daughter_index{ 0 };
daughter_index < daughterNames.size(); daughter_index++ ) {
decFile << " " << aliasPrefix[daughter_index]
<< daughterNames[daughter_index];
}
decFile << " " << modelNames[0];
- for ( auto par : parameters[0] ) {
+ for ( const auto& par : parameters[0] ) {
decFile << " " << par;
}
decFile << ";" << std::endl;
decFile << "Enddecay" << std::endl;
if ( doConjDecay[0] ) {
EvtId parID{ EvtPDL::getId( parent ) };
EvtId parConjID{ EvtPDL::chargeConj( parID ) };
decFile << "CDecay " << parConjID.getName() << std::endl;
}
// Daughter decays into granddaughters
- for ( long unsigned int daughter_index = 0;
+ for ( long unsigned int daughter_index{ 0 };
daughter_index < grandDaughterNames.size(); daughter_index++ ) {
if ( grandDaughterNames[daughter_index].empty() )
continue;
decFile << "Decay " << aliasPrefix[daughter_index]
<< daughterNames[daughter_index] << std::endl;
decFile << "1.0";
- for ( long unsigned int grandDaughter_index = 0;
+ for ( long unsigned int grandDaughter_index{ 0 };
grandDaughter_index < grandDaughterNames[daughter_index].size();
grandDaughter_index++ ) {
decFile << " "
<< grandDaughterNames[daughter_index][grandDaughter_index];
}
decFile << " " << modelNames[daughter_index + 1];
- for ( auto par : parameters[daughter_index + 1] ) {
+ for ( const auto& par : parameters[daughter_index + 1] ) {
decFile << " " << par;
}
decFile << ";" << std::endl;
decFile << "Enddecay" << std::endl;
if ( doConjDecay[daughter_index + 1] ) {
EvtId daugID{ EvtPDL::getId( daughterNames[daughter_index] ) };
EvtId daugConjID{ EvtPDL::chargeConj( daugID ) };
std::string conjName_Alias{ daugConjID.getName() };
if ( std::find( std::begin( daughterNames ),
std::end( daughterNames ), daugConjID.getName() ) !=
std::end( daughterNames ) ) {
conjName_Alias = conjName_Alias + "_" + daughter_index;
}
decFile << "CDecay " << aliasPrefix[daughter_index]
<< conjName_Alias << std::endl;
}
}
decFile << "End" << std::endl;
decFile.close();
return decName;
}
void TestDecayModel::defineHistos( TFile* outFile )
{
// Histogram information
- json jHistos = m_config["histograms"];
- size_t nHistos = jHistos.size();
+ const json& jHistos{ m_config.at( "histograms" ) };
+ const size_t nHistos{ jHistos.size() };
m_1DhistVect.reserve( nHistos );
m_2DhistVect.reserve( nHistos );
- for ( auto hInfo : jHistos ) {
- const auto varTitle = hInfo["title"].get<std::string>();
+ for ( const auto& hInfo : jHistos ) {
+ const auto varTitle{ hInfo.at( "title" ).get<std::string>() };
- const auto varName = hInfo["variable"].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["d1"].get<int>();
- const auto d2 = hInfo["d2"].get<int>();
+ const auto d1{ hInfo.at( "d1" ).get<int>() };
+ const auto d2{ hInfo.at( "d2" ).get<int>() };
- const auto nBins = hInfo["nbins"].get<int>();
- const auto xmin = hInfo["xmin"].get<double>();
- const auto xmax = hInfo["xmax"].get<double>();
+ 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" ) ) {
- TH1D* hist = new TH1D( histName.c_str(), varTitle.c_str(), nBins,
- xmin, xmax );
+ TH1* hist{ new TH1D{ histName.c_str(), varTitle.c_str(), nBins,
+ xmin, xmax } };
hist->SetDirectory( outFile );
m_1DhistVect.emplace_back(
std::make_pair( TestInfo( varName, d1, d2 ), hist ) );
continue;
} else {
- const auto varNameY = hInfo["variableY"].get<std::string>();
- const auto d1Y = hInfo["d1Y"].get<int>();
- const auto d2Y = hInfo["d2Y"].get<int>();
+ 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["nbinsY"].get<int>();
- const auto ymin = hInfo["ymin"].get<double>();
- const auto ymax = hInfo["ymax"].get<double>();
+ 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 );
}
- TH2D* hist = new TH2D( histName.c_str(), varTitle.c_str(), nBins,
- xmin, xmax, nBinsY, ymin, ymax );
+ TH2* hist{ new TH2D{ histName.c_str(), varTitle.c_str(), nBins,
+ xmin, xmax, nBinsY, ymin, ymax } };
hist->SetDirectory( outFile );
m_2DhistVect.emplace_back( std::make_pair(
TestInfo( varName, d1, d2, varNameY, d1Y, d2Y ), hist ) );
}
}
// For the case where the parent is either a neutral B or D add a mixed/unmixed histogram
- const auto parentName = m_config["parent"].get<std::string>();
- const int parentID = abs( EvtPDL::getStdHep( EvtPDL::getId( parentName ) ) );
+ const auto parentName{ m_config.at( "parent" ).get<std::string>() };
+ const int parentID{ abs( EvtPDL::getStdHep( EvtPDL::getId( parentName ) ) ) };
if ( parentID == 511 || parentID == 531 || parentID == 421 ) {
const std::string varTitle{ parentName + " mixed" };
- m_mixedHist = new TH1D( "mixed", varTitle.c_str(), 2, 0.0, 2.0 );
+ m_mixedHist = new TH1D{ "mixed", varTitle.c_str(), 2, 0.0, 2.0 };
// TODO maybe set bin labels?
m_mixedHist->SetDirectory( outFile );
}
}
void TestDecayModel::generateEvents( EvtGen& theGen, const std::string& decFile,
const std::string& parentName,
- bool doConjDecay, int nEvents,
- bool debug_flag )
+ const bool doConjDecay, const int nEvents,
+ const bool debug_flag )
{
// Read the decay file
theGen.readUDecay( decFile.c_str() );
// Generate the decays
- EvtId parId = EvtPDL::getId( parentName.c_str() );
- EvtId conjId = doConjDecay ? EvtPDL::chargeConj( parId ) : parId;
- for ( int i = 0; i < nEvents; i++ ) {
+ EvtId parId{ EvtPDL::getId( parentName.c_str() ) };
+ EvtId conjId{ doConjDecay ? EvtPDL::chargeConj( parId ) : parId };
+ for ( int i{ 0 }; i < nEvents; i++ ) {
if ( i % 1000 == 0 ) {
std::cout << "Event " << nEvents - i << std::endl;
}
// Initial 4-momentum and particle
EvtVector4R pInit( EvtPDL::getMass( parId ), 0.0, 0.0, 0.0 );
EvtParticle* parent{ nullptr };
if ( EvtRandom::Flat() < 0.5 ) {
parent = EvtParticleFactory::particleFactory( parId, pInit );
} else {
parent = EvtParticleFactory::particleFactory( conjId, pInit );
}
theGen.generateDecay( parent );
// Check for mixing (and fill histogram)
EvtParticle* prodParent{ nullptr };
if ( m_mixedHist ) {
if ( parent->getNDaug() == 1 ) {
prodParent = parent;
parent = prodParent->getDaug( 0 );
m_mixedHist->Fill( 1 );
} else {
m_mixedHist->Fill( 0 );
}
}
// To debug
if ( debug_flag ) {
std::cout << "Parent PDG code: " << parent->getPDGId()
<< " has daughters " << parent->getNDaug() << std::endl;
- for ( size_t iDaughter = 0; iDaughter < parent->getNDaug();
+ for ( size_t iDaughter{ 0 }; iDaughter < parent->getNDaug();
iDaughter++ ) {
std::cout << "Parent PDG code of daughter " << iDaughter
<< " : " << parent->getDaug( iDaughter )->getPDGId()
<< " has daughters "
<< parent->getDaug( iDaughter )->getNDaug()
<< std::endl;
- for ( size_t iGrandDaughter = 0;
+ for ( size_t iGrandDaughter{ 0 };
iGrandDaughter < parent->getDaug( iDaughter )->getNDaug();
iGrandDaughter++ ) {
std::cout << "Parent PDG code of grand daughter "
<< iGrandDaughter << " : "
<< parent->getDaug( iDaughter )
->getDaug( iGrandDaughter )
->getPDGId()
<< " has daughters "
<< parent->getDaug( iDaughter )
->getDaug( iGrandDaughter )
->getNDaug()
<< std::endl;
}
}
}
// Store information
- for ( auto& hist : m_1DhistVect ) {
- const auto& info = hist.first;
- const double value = getValue( parent, info.getName(), info.getd1(),
- info.getd2() );
+ for ( auto& [info, hist] : m_1DhistVect ) {
+ const double value{ getValue( parent, info.getName(), info.getd1(),
+ info.getd2() ) };
- if ( hist.second ) {
- hist.second->Fill( value );
+ if ( hist ) {
+ hist->Fill( value );
}
}
- for ( auto& hist : m_2DhistVect ) {
- const auto& info = hist.first;
- const double valueX = getValue( parent, info.getName(),
- info.getd1(), info.getd2() );
- const double valueY = getValue( parent, info.getName( 2 ),
- info.getd1( 2 ), info.getd2( 2 ) );
- if ( hist.second ) {
- hist.second->Fill( valueX, valueY );
+ for ( auto& [info, hist] : m_2DhistVect ) {
+ const double valueX{ getValue( parent, info.getName(), info.getd1(),
+ info.getd2() ) };
+ const double valueY{ getValue( parent, info.getName( 2 ),
+ info.getd1( 2 ), info.getd2( 2 ) ) };
+ if ( hist ) {
+ hist->Fill( valueX, valueY );
}
}
if ( debug_flag ) {
if ( prodParent ) {
prodParent->printTree();
} else {
parent->printTree();
}
}
// Cleanup
if ( prodParent ) {
prodParent->deleteTree();
} else {
parent->deleteTree();
}
}
}
-double TestDecayModel::getValue( EvtParticle* parent, const std::string& varName,
- const int d1, const int d2 ) const
+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();
+ 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
- EvtParticle* selectedParent = parent;
- std::string selectedVarName = varName;
+ 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';
+ 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;
+ 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();
+ 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();
+ const EvtVector4R p1_lab{ par1 != nullptr ? par1->getP4Lab() : EvtVector4R() };
+ const EvtVector4R p2_lab{ par2 != nullptr ? par2->getP4Lab() : EvtVector4R() };
+ const EvtVector4R p3_lab{ par3 != nullptr ? par3->getP4Lab() : EvtVector4R() };
if ( !selectedVarName.compare( "id" ) ) {
// StdHep ID of one of the daughters (controlled by d1) or the parent
if ( par1 ) {
value = par1->getPDGId();
} else {
value = selectedParent->getPDGId();
}
} else if ( !selectedVarName.compare( "parMass" ) ) {
// Parent invariant mass
value = selectedParent->mass();
} else if ( !selectedVarName.compare( "mass" ) ) {
// Invariant mass
if ( d2 != 0 ) {
// Invariant 4-mass combination of particles d1 and d2
value = ( p1 + p2 ).mass();
} else {
// Invariant mass of particle d1 only
value = p1.mass();
}
} else if ( !selectedVarName.compare( "massSq" ) ) {
// Invariant mass
if ( d2 != 0 ) {
// Invariant 4-mass combination of particles d1 and d2
value = ( p1 + p2 ).mass2();
} else {
// Invariant mass of particle d1 only
value = p1.mass2();
}
} else if ( !selectedVarName.compare( "mPrime" ) ) {
// pick up first unused particle rather than last one
if ( sel_NDaugMax != 3 ) {
return -1;
}
- int unused = 0;
- for ( int ii = 1; ii <= sel_NDaugMax; ++ii ) {
+ 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 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;
+ 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 ) {
+ 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 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();
+ 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 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 );
+ 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 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();
+ 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;
+ 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 );
+ 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 );
+ const double p1_x{ p1.get( 1 ) };
+ const double p1_y{ p1.get( 2 ) };
+ const double p1_z{ p1.get( 3 ) };
value = p1_x * p1_x + p1_y * p1_y + p1_z * p1_z;
} else if ( !selectedVarName.compare( "pxLab" ) ) {
// x momentum of particle d1 in lab frame
value = p1_lab.get( 1 );
} else if ( !selectedVarName.compare( "px" ) ) {
// x momentum of particle d1 in parent frame
value = p1.get( 1 );
} else if ( !selectedVarName.compare( "pyLab" ) ) {
// y momentum of particle d1 in lab frame
value = p1_lab.get( 2 );
} else if ( !selectedVarName.compare( "py" ) ) {
// y momentum of particle d1 in parent frame
value = p1.get( 2 );
} else if ( !selectedVarName.compare( "pzLab" ) ) {
// z momentum of particle d1 in lab frame
value = p1_lab.get( 3 );
} else if ( !selectedVarName.compare( "pz" ) ) {
// z momentum of particle d1 in parent frame
value = p1.get( 3 );
} else if ( !selectedVarName.compare( "cosHel" ) ||
!selectedVarName.compare( "absCosHel" ) ) {
// Cosine of helicity angle
EvtVector4R p12;
EvtVector4R p1Res;
if ( d2 != 0 ) {
// Resonance center-of-mass system (d1 and d2)
p12 = p1_lab + p2_lab;
// Boost vector
- const EvtVector4R boost( p12.get( 0 ), -p12.get( 1 ), -p12.get( 2 ),
- -p12.get( 3 ) );
+ 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;
+ 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();
+ 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;
+ 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();
+ 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 ) );
+ 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 );
+ 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();
+ 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();
+ const EvtVector4R p_B{ selectedParent->getP4Lab() };
// Index d2 must match the index of the pion (daughter of tau)
- const EvtParticle* pion_1 =
+ 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();
+ : selectedParent->getDaug( d1 - 1 ) };
+ const EvtVector4R p_pion_1{ pion_1 != nullptr ? pion_1->getP4Lab()
+ : EvtVector4R() };
- const EvtVector4R p_first_daughter =
+ const EvtVector4R p_first_daughter{
sel_NDaugMax <= 3 ? selectedParent->getDaug( d1 - 1 )->getP4Lab()
: selectedParent->getDaug( d1 - 1 )->getP4Lab() +
- selectedParent->getDaug( d1 )->getP4Lab();
+ selectedParent->getDaug( d1 )->getP4Lab() };
// Boost vector to tau frame
- const EvtVector4R boost_1( p_first_daughter.get( 0 ),
+ const EvtVector4R boost_1{ p_first_daughter.get( 0 ),
-p_first_daughter.get( 1 ),
-p_first_daughter.get( 2 ),
- -p_first_daughter.get( 3 ) );
+ -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 );
+ 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();
+ 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();
+ 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();
+ 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 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 ),
+ const EvtVector4R boost{ p_second_daughter.get( 0 ),
-p_second_daughter.get( 1 ),
-p_second_daughter.get( 2 ),
- -p_second_daughter.get( 3 ) );
+ -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 );
+ 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();
+ 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();
+ 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();
+ 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( "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 EvtVector4R p{ selectedParent->getP4() };
+ const EvtVector4R q{ p1 + p2 };
+ const EvtVector4R d{ p1 };
- const double cost = EvtDecayAngle( p, q, d );
+ 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 EvtVector4R p{ selectedParent->getP4() };
+ const EvtVector4R q{ p1 + p2 + p3 };
+ const EvtVector4R d{ p1 + p2 };
- const double cost = EvtDecayAngle( p, q, d );
+ 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 EvtVector4R p{ p1 + p2 + p3 };
+ const EvtVector4R q{ p1 + p2 };
+ const EvtVector4R d{ p1 };
- const double cost = EvtDecayAngle( p, q, d );
+ 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 );
+ 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( "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 =
+ const double E_over_Eparent_2{
parent->getDaug( 0 )->getDaug( d2 - 1 )->getP4Lab().get( 0 ) /
- parent->getDaug( 0 )->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
- double* dProb = selectedParent->decayProb();
+ 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;
+ 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();
+ 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 {
std::cerr << "Warning: Did not recognise variable name "
<< selectedVarName << std::endl;
}
return value;
}
void TestDecayModel::compareHistos( const std::string& refFileName ) const
{
// Compare histograms with the same name, calculating the chi-squared
- TFile* refFile = TFile::Open( refFileName.c_str(), "read" );
+ std::unique_ptr<TFile> refFile{ TFile::Open( refFileName.c_str(), "read" ) };
if ( !refFile ) {
std::cerr << "Could not open reference file " << refFileName << std::endl;
return;
}
- for ( auto& hist : m_1DhistVect ) {
- const std::string histName = hist.second->GetName();
+ // TODO - should we plot the (signed) chisq histogram? and save it as pdf/png?
+
+ for ( auto& [_, hist] : m_1DhistVect ) {
+ const std::string histName{ hist->GetName() };
// Get equivalent reference histogram
- const TH1* refHist = dynamic_cast<TH1*>(
- refFile->Get( histName.c_str() ) );
+ const TH1* refHist{
+ dynamic_cast<TH1*>( refFile->Get( histName.c_str() ) ) };
if ( refHist ) {
- double chiSq( 0.0 );
- int nDof( 0 );
- int iGood( 0 );
- double pValue = refHist->Chi2TestX( hist.second, chiSq, nDof, iGood,
- "WW" );
+ double chiSq{ 0.0 };
+ int nDof{ 0 };
+ int iGood{ 0 };
+ double pValue{ refHist->Chi2TestX( hist, chiSq, nDof, iGood, "WW" ) };
std::cout << "Histogram " << histName << " chiSq/nDof = " << chiSq
<< "/" << nDof << ", pValue = " << pValue << std::endl;
} else {
std::cerr << "Could not find reference histogram " << histName
<< std::endl;
}
}
- for ( auto& hist : m_2DhistVect ) {
- const std::string histName = hist.second->GetName();
+
+ for ( auto& [_, hist] : m_2DhistVect ) {
+ const std::string histName{ hist->GetName() };
// Get equivalent reference histogram
- const TH2* refHist = dynamic_cast<TH2*>(
- refFile->Get( histName.c_str() ) );
+ const TH2* refHist{
+ dynamic_cast<TH2*>( refFile->Get( histName.c_str() ) ) };
if ( refHist ) {
- double chiSq( 0.0 );
- int nDof( 0 );
- int iGood( 0 );
- double pValue = refHist->Chi2TestX( hist.second, chiSq, nDof, iGood,
- "WW" );
+ double chiSq{ 0.0 };
+ int nDof{ 0 };
+ int iGood{ 0 };
+ double pValue{ refHist->Chi2TestX( hist, chiSq, nDof, iGood, "WW" ) };
std::cout << "Histogram " << histName << " chiSq/nDof = " << chiSq
<< "/" << nDof << ", pValue = " << pValue << std::endl;
} else {
std::cerr << "Could not find reference histogram " << histName
<< std::endl;
}
}
refFile->Close();
}
-double TestDecayModel::getCosAcoplanarityAngle( EvtParticle* selectedParent,
- int sel_NDaugMax, int d1,
- int d2 ) const
+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 );
+ 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* 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 EvtParticle* grandDaughter1{ daughter1->getDaug( d1 - 1 ) };
+ const EvtParticle* grandDaughter2{ daughter2->getDaug( d2 - 1 ) };
- const EvtVector4R parent4Vector = selectedParent->getP4Lab();
+ const EvtVector4R parent4Vector{ selectedParent->getP4Lab() };
- const EvtVector4R daughter4Vector1 = sel_NDaugMax <= 3
- ? daughter1->getP4Lab()
- : daughter1->getP4Lab() +
- daughter2->getP4Lab();
+ const EvtVector4R daughter4Vector1{
+ sel_NDaugMax <= 3 ? daughter1->getP4Lab()
+ : daughter1->getP4Lab() + daughter2->getP4Lab() };
- const EvtVector4R daughter4Vector2 =
+ const EvtVector4R daughter4Vector2{
sel_NDaugMax == 2
? daughter2->getP4Lab()
: sel_NDaugMax == 3 ? daughter2->getP4Lab() + daughter3->getP4Lab()
- : daughter3->getP4Lab() + daughter4->getP4Lab();
+ : daughter3->getP4Lab() + daughter4->getP4Lab() };
- const EvtVector4R grandDaughter4Vector1 = sel_NDaugMax <= 3
- ? grandDaughter1->getP4Lab()
- : daughter1->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 grandDaughter4Vector2{
+ sel_NDaugMax == 2 ? grandDaughter2->getP4Lab()
+ : sel_NDaugMax == 3 ? daughter2->getP4Lab()
+ : daughter3->getP4Lab() };
- const EvtVector4R parentBoost( parent4Vector.get( 0 ),
+ const EvtVector4R parentBoost{ parent4Vector.get( 0 ),
-parent4Vector.get( 1 ),
-parent4Vector.get( 2 ),
- -parent4Vector.get( 3 ) );
+ -parent4Vector.get( 3 ) };
- const EvtVector4R daughter1Boost( daughter4Vector1.get( 0 ),
+ const EvtVector4R daughter1Boost{ daughter4Vector1.get( 0 ),
-daughter4Vector1.get( 1 ),
-daughter4Vector1.get( 2 ),
- -daughter4Vector1.get( 3 ) );
+ -daughter4Vector1.get( 3 ) };
- const EvtVector4R daughter2Boost( daughter4Vector2.get( 0 ),
+ const EvtVector4R daughter2Boost{ daughter4Vector2.get( 0 ),
-daughter4Vector2.get( 1 ),
-daughter4Vector2.get( 2 ),
- -daughter4Vector2.get( 3 ) );
+ -daughter4Vector2.get( 3 ) };
// Boosting daughters to reference frame of the mother
- EvtVector4R daughter4Vector1_boosted = boostTo( daughter4Vector1,
- parentBoost );
- EvtVector4R daughter4Vector2_boosted = boostTo( daughter4Vector2,
- parentBoost );
+ EvtVector4R daughter4Vector1_boosted{
+ boostTo( daughter4Vector1, parentBoost ) };
+ EvtVector4R daughter4Vector2_boosted{
+ boostTo( daughter4Vector2, parentBoost ) };
// Boosting each granddaughter to reference frame of its mother
- EvtVector4R grandDaughter4Vector1_boosted = boostTo( grandDaughter4Vector1,
- daughter1Boost );
- EvtVector4R grandDaughter4Vector2_boosted = boostTo( grandDaughter4Vector2,
- daughter2Boost );
+ EvtVector4R grandDaughter4Vector1_boosted{
+ boostTo( grandDaughter4Vector1, daughter1Boost ) };
+ 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 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();
+ const double normalVector1Mag{ normalVector1.d3mag() };
+ const double normalVector2Mag{ normalVector2.d3mag() };
if ( normalVector1Mag > 0.0 && normalVector2Mag > 0.0 ) {
cosAco = normalVector1.dot( normalVector2 ) /
( normalVector1Mag * normalVector2Mag );
}
return cosAco;
}
int main( int argc, char* argv[] )
{
if ( argc != 2 ) {
std::cerr << "Expecting one argument: json input file" << std::endl;
return 1;
}
/*! Load input file in json format. */
json config;
- std::ifstream inputStr( argv[1] );
+ std::ifstream inputStr{ argv[1] };
inputStr >> config;
inputStr.close();
bool allOK{ true };
if ( config.is_array() ) {
for ( const auto& cc : config ) {
- TestDecayModel test( cc );
+ TestDecayModel test{ cc };
allOK &= test.run();
}
} else {
- TestDecayModel test( config );
+ TestDecayModel test{ config };
allOK &= test.run();
}
return allOK ? 0 : 1;
}
diff --git a/test/testDecayModel.hh b/test/testDecayModel.hh
index 96b4b49..7082024 100644
--- a/test/testDecayModel.hh
+++ b/test/testDecayModel.hh
@@ -1,101 +1,102 @@
/***********************************************************************
* Copyright 1998-2020 CERN for the benefit of the EvtGen authors *
* *
* This file is part of EvtGen. *
* *
* EvtGen is free software: you can redistribute it and/or modify *
* it under the terms of the GNU General Public License as published by *
* the Free Software Foundation, either version 3 of the License, or *
* (at your option) any later version. *
* *
* EvtGen is distributed in the hope that it will be useful, *
* but WITHOUT ANY WARRANTY; without even the implied warranty of *
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the *
* GNU General Public License for more details. *
* *
* You should have received a copy of the GNU General Public License *
* along with EvtGen. If not, see <https://www.gnu.org/licenses/>. *
***********************************************************************/
#ifndef TEST_DECAY_MODEL_HH
#define TEST_DECAY_MODEL_HH
#include "TFile.h"
#include "TH1D.h"
#include "TH2D.h"
#include "nlohmann/json.hpp"
#include <cmath>
#include <string>
#include <utility>
#include <vector>
class EvtGen;
class EvtParticle;
class TestInfo {
public:
TestInfo( const std::string& name, const int d1, const int d2,
const std::string& nameY = "", const int d1Y = 0,
const int d2Y = 0 ) :
m_name{ name }, m_nameY{ nameY }, m_d1{ d1 }, m_d2{ d2 }, m_d1Y{ d1Y }, m_d2Y{ d2Y }
{
}
const std::string getName( const int var = 1 ) const
{
return var == 2 ? m_nameY : m_name;
}
int getd1( const int var = 1 ) const { return var == 2 ? m_d1Y : m_d1; }
int getd2( const int var = 1 ) const { return var == 2 ? m_d2Y : m_d2; }
private:
std::string m_name;
std::string m_nameY;
int m_d1;
int m_d2;
int m_d1Y;
int m_d2Y;
};
class TestDecayModel {
public:
TestDecayModel( const nlohmann::json& config );
bool run();
private:
bool checkMandatoryFields();
std::string createDecFile(
const std::string& parent, const std::vector<std::string>& daughterNames,
const std::vector<std::vector<std::string>>& grandDaughterNames,
const std::vector<std::string>& models,
const std::vector<std::vector<std::string>>& parameters,
const std::vector<bool>& doConjDecay,
const std::vector<std::string>& extras,
const std::string decFileName ) const;
void defineHistos( TFile* theFile );
void generateEvents( EvtGen& theGen, const std::string& decFile,
- const std::string& parentName, bool doConjDecay,
- int nEvents, bool debugFlag );
+ const std::string& parentName, const bool doConjDecay,
+ const int nEvents, const bool debugFlag );
- double getValue( EvtParticle* rootPart, const std::string& varName,
+ double getValue( const EvtParticle* rootPart, const std::string& varName,
const int d1, const int d2 ) const;
void compareHistos( const std::string& refFileName ) const;
- double getCosAcoplanarityAngle( EvtParticle* selectedParent,
- int sel_NDaugMax, int d1, int d2 ) const;
+ double getCosAcoplanarityAngle( const EvtParticle* selectedParent,
+ const int sel_NDaugMax, const int d1,
+ const int d2 ) const;
const nlohmann::json& m_config;
- std::vector<std::pair<TestInfo, TH1D*>> m_1DhistVect;
- std::vector<std::pair<TestInfo, TH2D*>> m_2DhistVect;
- TH1D* m_mixedHist{ nullptr };
+ std::vector<std::pair<TestInfo, TH1*>> m_1DhistVect;
+ std::vector<std::pair<TestInfo, TH2*>> m_2DhistVect;
+ TH1* m_mixedHist{ nullptr };
};
#endif
diff --git a/validation/genRootDecayChain.cc b/validation/genRootDecayChain.cc
index 88a8ba6..88d5066 100644
--- a/validation/genRootDecayChain.cc
+++ b/validation/genRootDecayChain.cc
@@ -1,406 +1,406 @@
/***********************************************************************
* 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/>. *
***********************************************************************/
// Program to create ROOT files for EvtGen validation plots,
// generalising storage of daughter information for each decay chain step.
// There are three integers to keep track of the decay tree history:
// dVtx = present decay vertex integer number
// pVtx = parent decay vertex integer number
// daug = daughter order number for the given particle
//
// Consider the following decay chain history, where Vj is vertex number j:
// V1
// A -> B + C <---- First chain level
// | |
// +--> B1 + B2 +--> C1 C2 C3 <---- Second chain level
// V2 | V5 |
// +--> B1a B1b +--> C1a C1b C1c <---- Third chain level
// V3 | V6
// +--> d1 d2 <---- Fourth chain level
// V4
//
// Since the decay tree is stored recursively, the information order follows
// the immediate, sequential, sub-decay vertex order V1, V2, V3, V4, V5, V6:
// A, B, B1, B1a, B1b, d1, d2, B2, C, C1, C1a, C1b, C1c, C2, C3
//
// The unique set of integer "attributes" for each particle is as follows:
//
// Particle dVtx pVtx daug
// A 0 0 0 The mother particle
// B 1 0 1
// B1 2 1 1
// B1a 3 2 1
// B1b 3 2 2
// d1 4 3 1
// d2 4 3 2
// B2 2 1 2
// C 1 0 2
// C1 5 1 1
// C1a 6 5 1
// C1b 6 5 2
// C1c 6 5 3
// C2 5 1 2
// C3 5 1 3
//
// The decay tree will be defined by the decay file. Additional photons from PHOTOS will
// only appear as appended, extra daughters on the RHS of a given vertex, e.g. if a FSR
// photon is added to the RHS of V4, there will be a "d3" particle with attributes (4,3,3)
#include "genRootDecayChain.hh"
#include "EvtGen/EvtGen.hh"
#include "EvtGenBase/EvtAbsRadCorr.hh"
#include "EvtGenBase/EvtDecayBase.hh"
#include "EvtGenBase/EvtMTRandomEngine.hh"
#include "EvtGenBase/EvtPDL.hh"
#include "EvtGenBase/EvtParticle.hh"
#include "EvtGenBase/EvtParticleFactory.hh"
#include "EvtGenBase/EvtPatches.hh"
#include "EvtGenBase/EvtRandomEngine.hh"
#include "EvtGenBase/EvtSimpleRandomEngine.hh"
#ifdef EVTGEN_EXTERNAL
#include "EvtGenExternal/EvtExternalGenList.hh"
#endif
#include "TROOT.h"
#include "TStyle.h"
#include <iostream>
#include <list>
#include <string>
using std::cout;
using std::endl;
using std::string;
genRootDecayChain::genRootDecayChain( const string& decayFileName,
const string& rootFileName,
const string& parentName, int nEvents,
bool storeMtmXYZ ) :
_decayFileName( decayFileName ),
_rootFileName( rootFileName ),
_parentName( parentName ),
_nEvents( nEvents ),
_storeMtmXYZ( storeMtmXYZ ),
_theFile( 0 ),
_theTree( 0 ),
_probHist( 0 ),
_theCanvas( 0 )
{
_theFile = new TFile( rootFileName.c_str(), "recreate" );
_theTree = new TTree( "Data", "Data" );
_theTree->SetDirectory( _theFile );
_probHist = new TH1D( "probHist", "probHist", 100, 0.0, 0.0 );
_probHist->SetXTitle( "Prob/MaxProb" );
_probHist->SetDirectory( _theFile );
gROOT->SetStyle( "Plain" );
gStyle->SetOptStat( 0 );
_theCanvas = new TCanvas( "theCanvas", "", 900, 700 );
_theCanvas->Clear();
_theCanvas->UseCurrentStyle();
}
genRootDecayChain::~genRootDecayChain()
{
_theFile->Close();
delete _theCanvas;
}
void genRootDecayChain::run()
{
this->initTree();
this->generateEvents();
this->writeTree();
}
void genRootDecayChain::initTree()
{
// Initialise internal variables
_eventId = 0;
_PDGId = 0;
_dVtx = 0, _pVtx = 0, _daug = 0;
_px = 0.0, _py = 0.0, _pz = 0.0;
_p = 0.0, _E = 0.0, _m = 0.0, _t = 0.0;
// Set-up the TTree
_theTree->Branch( "eventId", &_eventId, "eventId/I" );
_theTree->Branch( "PDGId", &_PDGId, "PDGId/I" );
_theTree->Branch( "dVtx", &_dVtx, "dVtx/I" );
_theTree->Branch( "pVtx", &_pVtx, "pVtx/I" );
_theTree->Branch( "daug", &_daug, "daug/I" );
_theTree->Branch( "p", &_p, "p/D" );
_theTree->Branch( "E", &_E, "E/D" );
_theTree->Branch( "pL", &_pL, "pL/D" );
_theTree->Branch( "EL", &_EL, "EL/D" );
_theTree->Branch( "m", &_m, "m/D" );
if ( _storeMtmXYZ ) {
// Store momentum components as well
_theTree->Branch( "px", &_px, "px/D" );
_theTree->Branch( "py", &_py, "py/D" );
_theTree->Branch( "pz", &_pz, "pz/D" );
_theTree->Branch( "pxL", &_pxL, "pxL/D" );
_theTree->Branch( "pyL", &_pyL, "pyL/D" );
_theTree->Branch( "pzL", &_pzL, "pzL/D" );
}
// Lifetime
_theTree->Branch( "t", &_t, "t/D" );
}
void genRootDecayChain::writeTree()
{
_theFile->cd();
_theTree->Write();
_probHist->Write();
}
void genRootDecayChain::generateEvents()
{
EvtRandomEngine* randomEngine = 0;
EvtAbsRadCorr* radCorrEngine = 0;
std::list<EvtDecayBase*> extraModels;
// Define the random number generator
#ifdef EVTGEN_CPP11
// Use the Mersenne-Twister generator (C++11 only)
randomEngine = new EvtMTRandomEngine();
#else
randomEngine = new EvtSimpleRandomEngine();
#endif
// Initialize the generator - read in the decay table and particle properties.
// For our validation purposes, we just want to read in one decay file
#ifdef EVTGEN_EXTERNAL
bool convertPythiaCodes( false );
bool useEvtGenRandom( true );
EvtExternalGenList genList( convertPythiaCodes, "", "gamma", useEvtGenRandom );
radCorrEngine = genList.getPhotosModel();
extraModels = genList.getListOfModels();
#endif
int mixingType( 1 );
bool useXml( false );
EvtGen evtGen( _decayFileName.c_str(), "../evt.pdl", randomEngine,
radCorrEngine, &extraModels, mixingType, useXml );
EvtParticle* theParent( 0 );
EvtId theId = EvtPDL::getId( _parentName );
if ( theId.getId() == -1 && theId.getAlias() == -1 ) {
cout << "Error. Could not find valid EvtId for " << _parentName << endl;
return;
}
// Loop to create nEvents
int i;
for ( i = 0; i < _nEvents; i++ ) {
// Parent particle 4-momentum
EvtVector4R pInit( EvtPDL::getMass( theId ), 0.0, 0.0, 0.0 );
_eventId = i;
if ( i % 1000 == 0 ) {
cout << "Event number = " << i + 1 << " out of " << _nEvents
<< std::endl;
}
theParent = EvtParticleFactory::particleFactory( theId, pInit );
if ( theParent->getSpinStates() == 3 ) {
theParent->setVectorSpinDensity();
}
// Generate the event
evtGen.generateDecay( theParent );
// Decay vertex index
theParent->setAttribute( "dVtx", 0 );
// Parent vertex index
theParent->setAttribute( "pVtx", 0 );
// Daughter index
theParent->setAttribute( "daug", 0 );
int nDaug = theParent->getNDaug();
int iDaug( 0 );
storeTreeInfo( theParent );
//cout<<"Event number = "<<i<<endl;
//theParent->printTree();
// Initialise the vertex number
_vertexNo = 1;
// Loop over the daughter tracks
for ( iDaug = 0; iDaug < nDaug; iDaug++ ) {
EvtParticle* daug = theParent->getDaug( iDaug );
// Decay vertex index
daug->setAttribute( "dVtx", 1 );
// Parent decay vertex
daug->setAttribute( "pVtx", 0 );
// Daughter index: 1,2,..,nDaug
daug->setAttribute( "daug", iDaug + 1 );
// Recursively store the daughter information
this->storeDaughterInfo( daug );
} // daughter loop
// Store probability/max probability
- double* dProb = theParent->decayProb();
+ const double* dProb = theParent->decayProb();
if ( dProb ) {
_probHist->Fill( *dProb );
}
theParent->deleteTree();
} // event loop
delete randomEngine;
}
void genRootDecayChain::storeDaughterInfo( EvtParticle* theParticle )
{
if ( !theParticle ) {
return;
}
// Store the particle information in the TTree
storeTreeInfo( theParticle );
// Loop over the particle's own daughter tracks and call this function
// recursively, keeping track of the decay chain order number
int nDaug = theParticle->getNDaug();
int iDaug( 0 );
// Increment the decay vertex integer if this particle has daughters
if ( nDaug > 0 ) {
_vertexNo++;
}
// The parent vertex index for the daughters is equal to the
// current particle decay vertex index
int parentVtx = theParticle->getAttribute( "dVtx" );
// First, we need to loop over the given particle's daughters and set the
// attributes so we keep track of the decay tree history at this chain level
for ( iDaug = 0; iDaug < nDaug; iDaug++ ) {
EvtParticle* daug = theParticle->getDaug( iDaug );
// Set the attributes for the daughters
daug->setAttribute( "dVtx", _vertexNo );
daug->setAttribute( "pVtx", parentVtx );
daug->setAttribute( "daug", iDaug + 1 );
}
// Now loop over the daughters again and recursively call this function to
// follow the rest of the decay tree history
for ( iDaug = 0; iDaug < nDaug; iDaug++ ) {
EvtParticle* daug = theParticle->getDaug( iDaug );
storeDaughterInfo( daug );
}
}
void genRootDecayChain::storeTreeInfo( EvtParticle* theParticle )
{
if ( !theParticle ) {
return;
}
// Store the particle information in the TTree by first setting the internal
// variables and then calling Fill()
_dVtx = theParticle->getAttribute( "dVtx" );
_pVtx = theParticle->getAttribute( "pVtx" );
_daug = theParticle->getAttribute( "daug" );
//cout<<"Particle "<<theParticle->getName()<<", dVtx = "<<_dVtx
//<<", pVtx = "<<_pVtx<<", daug = "<<_daug<<endl;
EvtVector4R p4Lab = theParticle->getP4Lab(); // lab frame = mother frame
EvtVector4R p4 = theParticle->getP4(); // frame = parents frame
// PDG id
_PDGId = EvtPDL::getStdHep( theParticle->getId() );
// 4-momenta in parent frame
_E = p4.get( 0 );
_px = p4.get( 1 );
_py = p4.get( 2 );
_pz = p4.get( 3 );
_p = sqrt( _px * _px + _py * _py + _pz * _pz );
// 4-momenta in lab frame
_EL = p4Lab.get( 0 );
_pxL = p4Lab.get( 1 );
_pyL = p4Lab.get( 2 );
_pzL = p4Lab.get( 3 );
_pL = sqrt( _pxL * _pxL + _pyL * _pyL + _pzL * _pzL );
// Rest mass and lifetime
_m = theParticle->mass();
_t = theParticle->getLifetime();
// Store the information in the TTree
_theTree->Fill();
}
int main( int argc, char** argv )
{
// Use the B0 -> K'*0 gamma, K'*0 -> K+ pi- decays as an example
string decayFileName( "BKstarGamma.dec" );
//string decayFileName("BuDst0rhop.dec");
if ( argc > 1 ) {
decayFileName = argv[1];
}
cout << "Decay file name is " << decayFileName << endl;
string rootFileName( "BKstarGamma.root" );
//string rootFileName("BuDst0rhop.root");
if ( argc > 2 ) {
rootFileName = argv[2];
}
cout << "Root file name is " << rootFileName << endl;
string parentName( "B0" );
//string parentName("B-");
if ( argc > 3 ) {
parentName = argv[3];
}
cout << "Parent name is " << parentName << endl;
int nEvents( 100000 );
if ( argc > 4 ) {
nEvents = atoi( argv[4] );
}
bool storeMtmXYZ = true;
genRootDecayChain myGen( decayFileName, rootFileName, parentName, nEvents,
storeMtmXYZ );
myGen.run();
return 0;
}