analysis

Collaboration diagram for analysis:

Modules
	analysis data objects
	analysis modules

Namespaces
namespace	Belle2::ParticleListName
	Helper to deal with ParticleList names.
namespace	Belle2::DistanceTools
	This namespace contains a collection of function that are useful to compute distances between tracks and vertices.
namespace	Belle2::EvtPDLUtil
	Utilities for converting PDG codes into particle names.
namespace	Belle2::ParticleCopy
	Functions that create copies of Particles.

Classes
class	ClusterUtils
	Class to provide momentum-related information from ECLClusters. More...
class	CleoCones
	Class to calculate the Cleo clone variables. More...
class	FoxWolfram
	Class to calculate the Fox-Wolfram moments up to order 8. More...
class	HarmonicMoments
	Class to calculate the Harmonic moments up to order 8 with respect to a given axis. More...
class	KsfwMoments
	Moment-calculation of the k_sfw improved Super-Fox-Wolfram moments. More...
class	SphericityEigenvalues
	Class to calculate the Sphericity tensor eigenvalues and eigenvectors starting from an array of 3-momenta The tensor itself is not stored, only its eigenvalues and eigenvectors are. More...
class	Thrust
	Class to calculate the thrust axis. More...
class	ChargedPidMVAWeights
	Class to contain the payload of MVA weightfiles needed for charged particle identification. More...
class	ECLPhotonEnergyResolution
	Class to hold the information ECL energy resolution derived from PERC. More...
class	ParticleWeightingAxis
	Class for handling LookUp tables. More...
class	ParticleWeightingBinLimits
	Just pair of numbers - min and max values of bin border. More...
class	ParticleWeightingKeyMap
	Class for handling KeyMap. More...
class	ParticleWeightingLookUpTable
	Class for handling LookUp tables. More...
class	PIDCalibrationWeight
	Class for handling the PID calibration weight matrix. More...
class	PIDDetectorWeights
	Class for handling the PID weights per detector, used to calculate the track helix isolation score per particle. More...
class	PIDNeuralNetworkParameters
	Class for handling the parameters for the neural-network PID. More...
class	PIDPriors
	A database class to hold the prior probability for the particle identification. More...
class	PIDPriorsTable
	This class holds the prior distribution for a single particle species. More...
class	DecayDescriptor
	The DecayDescriptor stores information about a decay tree or parts of a decay tree. More...
class	DecayDescriptorParticle
	Represents a particle in the DecayDescriptor. More...
struct	DecayStringDecay
	Holds the information of a decay. More...
struct	DecayStringGrammar< Iterator >
	This class describes the grammar and the syntax elements of decay strings. More...
struct	DecayStringParticle
	Holds the information of a particle in the decay string. More...
class	ParticleListHelper
	Class to help managing creation and adding to ParticleLists. More...
class	ParticleIndexGenerator
	ParticleIndexGenerator is a generator for all the combinations of the particle indices stored in the particle lists. More...
class	ListIndexGenerator
	ListIndexGenerator is a generator for all the combinations of the sublists (FlavorSpecificParticle = 0, SelfConjugatedParticle = 1) of a set of particle lists. More...
class	ParticleGenerator
	ParticleGenerator is a generator for all the particles combined from the given ParticleLists. More...
class	ChargedParticleIdentificatorTest
	Test the MVA-based charged PID. More...
class	eventShapeCoreAlgorithmTest
class	PIDPriorsTest
class	TrackIsoScoreCalculatorTest
	Test the calculation of the track helix-based isolation score per particle. More...
class	AnalysisConfiguration
	Singleton class keeping configurables of analysis components. More...
class	DecayForest
	Contains several DecayTree objects, which belong all to the same candidate. More...
class	DecayNode
	DecayNode describes the decay of a particle identified by its pdg code, into list of daughters. More...
class	DecayTree< T >
	This is a helper class for the MCDecayFinderModule. More...
class	DetSurfCylBoundaries
	Simple class to encapsulate a detector surface's boundaries in cylindrical coordinates. More...
struct	DetectorSurface
	Detector surfaces information. More...
class	GenB0Tag
	Class to determine generated decay mode of B0 and B0bar. More...
class	GenBplusTag
	Class to determine generated decay modes of B+ and B-. More...
class	GenBsTag
	Class to determine generated decay mode of Bs0 and Bs0bar. More...
class	GenDTag
	Class to determine generated decay mode of D*, Ds, D0, D+, and their anti-particles. More...
class	GenTauTag
	Class to determine generated decay mode of tau+ and tau-. More...
struct	KlongCalculatorUtils
	Utility class to calculate the Klong kinematics. More...
struct	MCMatching
	Functions to perform Monte Carlo matching for reconstructed Particles. More...
class	ParticleSubset
	Specialised SelectSubset<Particle> that also fixes daughter indices and all ParticleLists. More...
class	PCmsLabTransform
	Class to hold Lorentz transformations from/to CMS and boost vector. More...
class	PIDCalibrationWeightUtil
	Class to call calibration weight matrix. More...
class	PIDNeuralNetwork
	Class to call PID neural network. More...
class	PostProcessingParticleWeighting
	Post-processing particle weighting. More...
class	ReferenceFrame
	Abstract base class of all reference frames. More...
class	RestFrame
	Rest frame of a particle. More...
class	LabFrame
	Lab frame. More...
class	CMSFrame
	CMS frame. More...
class	RotationFrame
	Rotation frame around vector. More...
class	CMSRotationFrame
	Stack frame for cms and Rotation frame. More...
class	UseReferenceFrame< T >
	A template class to apply the reference frame. More...
struct	VariableFormulaConstructor
	Struct to construct new variable function objects from a name or a double value or to apply operations on these variable function objects. More...

Macros
#define	MAKE_DEPRECATED(name, make_fatal, version, description)    static DeprecateProxy VARMANAGER_MAKE_UNIQUE(_deprecateproxy)(std::string(name), bool(make_fatal), std::string(version), std::string(description));
	Registers a variable as deprecated.

Typedefs
typedef std::vector< std::pair< double, double > >	Binning
	Bin holder as vector for bin limit pairs: [energy limits, theta limits, phi limits].
typedef std::map< int, ParticleWeightingBinLimits * >	BinMap
	Map of keys with bin limits.
typedef std::pair< std::vector< int >, int >	MultiDimBin
	Multidimensional bin: first element contains combination of bin IDs from 1D axis, second elements contain ID ("key") associated with this combination.
typedef std::map< std::string, ParticleWeightingBinLimits * >	NDBin
	N-dim bin: pairs of bin limits with name of the axis variable.
typedef std::map< std::string, double >	WeightInfo
	Weight information: a line from the weight lookup table.
typedef std::map< int, WeightInfo >	WeightMap
	Weight map: the whole lookup table
typedef std::vector< std::vector< double > >	WeightMatrix
	PID calibration weight matrix, 6 (particle type) x 6 (detectors).
typedef std::vector< std::tuple< size_t, float > >	PIDNNMissingInputs
	Stores information on how to handle missing inputs, i.e.
typedef std::vector< std::tuple< size_t, size_t, double, double, double > >	PIDNNInputsToCut
	Stores information on whether and how to overwrite certain inputs.
typedef boost::variant< boost::recursive_wrapper< DecayStringDecay >, DecayStringParticle >	DecayString
	The DecayStringElement can be either a DecayStringDecay or a vector of mother particles.

Functions
void	addContinuumSuppression (const Particle *particle, const std::string &maskName)
	Adds continuum suppression variables.
double	legendre (const double z, const int i)
	Legendre polynomials.
	TEST_F (ChargedParticleIdentificatorTest, TestDBRep)
	Test correct storage of weightfiles in the database representation inner structure.
	TEST_F (eventShapeCoreAlgorithmTest, Thrust)
	Test the calculation of a thrust axis.
	TEST_F (eventShapeCoreAlgorithmTest, CleoCones)
	Test the calculation of the CleoClones variables.
	TEST_F (eventShapeCoreAlgorithmTest, FoxWolfram)
	Test the calculation of the Fox-Wolfram moments.
	TEST_F (eventShapeCoreAlgorithmTest, HarmonicMoments)
	Test the calculation of the Harmonic moments.
	TEST_F (eventShapeCoreAlgorithmTest, Sphericity)
	Test the calculation of the Sphericity eigenvalues and eigenvectors.
	TEST_F (PIDPriorsTest, PIDPriorsTableTest)
	Test of the PIDPriorsTable class.
	TEST_F (PIDPriorsTest, PIDPriorTest)
	Test the PIDPriors dbobject.
	TEST_F (TrackIsoScoreCalculatorTest, TestDBRep)
	Test correct retrieval of information from the database representation inner structure.
bool	operator== (const DecayNode &node1, const DecayNode &node2)
	Compare two Decay Nodes: They are equal if All daughter decay nodes are equal or one of the daughter lists is empty (which means "inclusive")
bool	operator!= (const DecayNode &node1, const DecayNode &node2)
	Not equal: See operator==.
	CleoCones (const std::vector< ROOT::Math::XYZVector > &p3_cms_all, const std::vector< ROOT::Math::XYZVector > &p3_cms_roe, const ROOT::Math::XYZVector &thrustB, bool calc_CleoCones_with_all, bool calc_CleoCones_with_roe)
	Constructor.
	KsfwMoments (double Hso0_max, std::vector< std::pair< ROOT::Math::XYZVector, int > > p3_cms_q_sigA, std::vector< std::pair< ROOT::Math::XYZVector, int > > p3_cms_q_sigB, std::vector< std::pair< ROOT::Math::XYZVector, int > > p3_cms_q_roe, const ROOT::Math::PxPyPzEVector &p_cms_missA, const ROOT::Math::PxPyPzEVector &p_cms_missB, const double et[2])
	Constructor.
void	registerList (const std::string &listname, bool save=true)
	Register a list by name.
void	registerList (const DecayDescriptor &decay, bool save=true)
	Register a list using a decay descriptor.
void	create ()
	Create the list objects.
void	init (const std::vector< unsigned int > &_sizes)
	Initialises the generator to produce combinations with the given sizes of each particle list.
bool	loadNext ()
	Loads the next combination.
const std::vector< unsigned int > &	getCurrentIndices () const
	Returns theindices of the current loaded combination.
void	init (unsigned int _numberOfLists)
	Initialises the generator to produce the given type of sublist.
bool	loadNext ()
	Loads the next combination.
const std::vector< ParticleList::EParticleType > &	getCurrentIndices () const
	Returns the type of the sublist of the current loaded combination.
	ParticleGenerator (const std::string &decayString, const std::string &cutParameter="")
	Initialises the generator to produce the given type of sublist.
	ParticleGenerator (const DecayDescriptor &decaydescriptor, const std::string &cutParameter="")
	Initialises the generator to produce the given type of sublist.
void	init ()
	Initialises the generator to produce the given type of sublist.
void	initIndicesToUniqueIDMap ()
	In the case input daughter particle lists collide (two or more lists contain copies of Particles) the Particle's Store Array index can not be longer used as its unique identifier, which is needed to check for uniqueness of accepted combinations.
void	fillIndicesToUniqueIDMap (const std::vector< int > &listA, const std::vector< int > &listB, int &uniqueID)
	Assigns unique IDs to all particles in list A, which do not have the unique ID already assigned.
void	fillIndicesToUniqueIDMap (const std::vector< int > &listA, int &uniqueID)
	Assigns unique IDs to all particles in list A, which do not have the unique ID already assigned.
bool	loadNext (bool loadAntiParticle=true)
	Loads the next combination.
bool	loadNextParticle (bool useAntiParticle)
	Loads the next combination.
bool	loadNextSelfConjugatedParticle ()
	Loads the next combination.
Particle	createCurrentParticle () const
	Create current particle object.
Particle	getCurrentParticle () const
	Returns the particle.
bool	currentCombinationHasDifferentSources ()
	Check that all FS particles of a combination differ.
bool	currentCombinationIsUnique ()
	Check that the combination is unique.
bool	inputListsCollide (const std::pair< unsigned, unsigned > &pair) const
	True if the pair of input lists collide.
bool	currentCombinationIsECLCRUnique ()
	Check that: if the current combination has at least two particles from an ECL source, then they are from different connected regions or from the same connected region but have the same hypothesis.
int	getUniqueID (int index) const
	Returns the unique ID assigned to Particle with given index from the IndicesToUniqueID map.
bool	find_decay (const DecayNode &to_find) const
	Check if the decay node contains the given decay tree.
std::string	print_node (unsigned int indent=0) const
	Output a single node.

Detailed Description

Macro Definition Documentation

MAKE_DEPRECATED

#define MAKE_DEPRECATED	(		name,
			make_fatal,
			version,
			description
	)		static DeprecateProxy VARMANAGER_MAKE_UNIQUE(_deprecateproxy)(std::string(name), bool(make_fatal), std::string(version), std::string(description));

Registers a variable as deprecated.

Definition at line 443 of file Manager.h.

Typedef Documentation

BinMap

typedef std::map<int, ParticleWeightingBinLimits*> BinMap

Map of keys with bin limits.

Definition at line 21 of file ParticleWeightingAxis.h.

Binning

typedef std::vector<std::pair<double, double> > Binning

Bin holder as vector for bin limit pairs: [energy limits, theta limits, phi limits].

Definition at line 29 of file ECLPhotonEnergyResolution.h.

DecayString

typedef boost::variant< boost::recursive_wrapper<DecayStringDecay>, DecayStringParticle > DecayString

The DecayStringElement can be either a DecayStringDecay or a vector of mother particles.

User documentation is located at analysis/doc/DecayDescriptor.rst Please modify in accordingly to introduced changes.

Definition at line 23 of file DecayString.h.

MultiDimBin

typedef std::pair<std::vector<int>, int> MultiDimBin

Multidimensional bin: first element contains combination of bin IDs from 1D axis, second elements contain ID ("key") associated with this combination.

Definition at line 27 of file ParticleWeightingKeyMap.h.

NDBin

typedef std::map<std::string, ParticleWeightingBinLimits*> NDBin

N-dim bin: pairs of bin limits with name of the axis variable.

Definition at line 32 of file ParticleWeightingKeyMap.h.

PIDNNInputsToCut

typedef std::vector< std::tuple<size_t, size_t, double, double, double> > PIDNNInputsToCut

Stores information on whether and how to overwrite certain inputs.

It is a vector of tuples of 5 elements:

• Element 0: Index i of input to be overwritten
• Element 1: Index j of input that defines if i is overwritten
• [Element1, Element2] give the range in j in which i is overwritten
• Element 4: Value to which i will be set For example, if input 0 should be overwritten to -1.0 if input 1 is within [1.0, 2.0]: (0, 1, 1.0, 2.0, -1.0)

Definition at line 40 of file PIDNeuralNetworkParameters.h.

PIDNNMissingInputs

typedef std::vector<std::tuple<size_t, float> > PIDNNMissingInputs

Stores information on how to handle missing inputs, i.e.

inputs that are NaN It is a vector of tuples of two elements, where the first element is the index i of the input variable and the second element is the value to which input i is set if it is NaN.

Definition at line 28 of file PIDNeuralNetworkParameters.h.

WeightInfo

typedef std::map<std::string, double> WeightInfo

Weight information: a line from the weight lookup table.

Definition at line 21 of file ParticleWeightingLookUpTable.h.

WeightMap

typedef std::map<int, WeightInfo> WeightMap

Weight map: the whole lookup table

Definition at line 23 of file ParticleWeightingLookUpTable.h.

WeightMatrix

typedef std::vector<std::vector<double> > WeightMatrix

PID calibration weight matrix, 6 (particle type) x 6 (detectors).

The particle types are sorted by their invariant mass (e, mu, pi, K, p, d), which is inherited from Const::chargedStableSet. Each vector<double> has 6 weights for different sub-detectors for a particle type. The order is as follows, [SVD, CDC, TOP, ARICH, ECL. KLM] that is inherited from Const::PIDDetectors.

Definition at line 30 of file PIDCalibrationWeight.h.

Function Documentation

addContinuumSuppression()

void addContinuumSuppression	(	const Particle *	particle,
		const std::string &	maskName
	)

Adds continuum suppression variables.

Definition at line 29 of file ContinuumSuppression.cc.

  {
    // Output
    StoreArray<ContinuumSuppression> qqArray(maskName);
    // Create ContinuumSuppression object
    ContinuumSuppression* qqVars = qqArray.appendNew();
 
    // Create relation: Particle <-> ContinuumSuppression
    particle->addRelationTo(qqVars);
 
    std::vector<ROOT::Math::XYZVector> p3_cms_sigB, p3_cms_roe, p3_cms_all;
 
    std::vector<std::pair<ROOT::Math::XYZVector, int>> p3_cms_q_sigA;
    std::vector<std::pair<ROOT::Math::XYZVector, int>> p3_cms_q_sigB;
    std::vector<std::pair<ROOT::Math::XYZVector, int>> p3_cms_q_roe;
 
    std::vector<float> ksfwFS0;
    std::vector<float> ksfwFS1;
 
    std::vector<float> cleoConesAll;
    std::vector<float> cleoConesROE;
 
    double et[2];
 
    ROOT::Math::XYZVector thrustB;
    ROOT::Math::XYZVector thrustO;
 
    float thrustBm = -1;
    float thrustOm = -1;
    float cosTBTO  = -1;
    float cosTBz   = -1;
    float R2       = -1;
 
 
    // -- B Cand --------------------------------------------------------------------------
    PCmsLabTransform T;
    double BeamEnergy = T.getCMSEnergy() / 2;
 
    ROOT::Math::PxPyPzEVector p_cms_missA(0, 0, 0, 2 * BeamEnergy);
    ROOT::Math::PxPyPzEVector p_cms_missB(0, 0, 0, 2 * BeamEnergy);
    et[0] = et[1] = 0;
 
    // -- SIG A --- Use B primary daughters - (Belle: use_finalstate_for_sig == 0) --------
    std::vector<Belle2::Particle*> signalDaughters = particle->getDaughters();
 
    for (const Belle2::Particle* sigFS0 : signalDaughters) {
      ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * sigFS0->get4Vector();
 
      p3_cms_q_sigA.emplace_back(p_cms.Vect(), sigFS0->getCharge());
 
      p_cms_missA -= p_cms;
      et[0] += p_cms.Pt();
    }
 
    // -- SIG B --- Use B final-state daughters - (Belle: use_finalstate_for_sig == 1) ----
    std::vector<const Belle2::Particle*> signalFSParticles = particle->getFinalStateDaughters();
 
    for (const Belle2::Particle* sigFS1 : signalFSParticles) {
      ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * sigFS1->get4Vector();
 
      p3_cms_all.push_back(p_cms.Vect());
      p3_cms_sigB.push_back(p_cms.Vect());
 
      p3_cms_q_sigB.emplace_back(p_cms.Vect(), sigFS1->getCharge());
 
      p_cms_missB -= p_cms;
      et[1] += p_cms.Pt();
    }
 
    // -- ROE -----------------------------------------------------------------------------
    const RestOfEvent* roe = particle->getRelated<RestOfEvent>();
 
    if (roe) {
 
      // Charged tracks
      //
      std::vector<const Particle*> chargedROEParticles = roe->getChargedParticles(maskName);
 
      for (const Particle* chargedROEParticle : chargedROEParticles) {
 
        // TODO: Add helix and KVF with IpProfile once available. Port from L163-199 of:
        // /belle/b20090127_0910/src/anal/ekpcontsuppress/src/ksfwmoments.cc
 
        ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * chargedROEParticle->get4Vector();
 
        p3_cms_all.push_back(p_cms.Vect());
        p3_cms_roe.push_back(p_cms.Vect());
 
        p3_cms_q_roe.emplace_back(p_cms.Vect(), chargedROEParticle->getCharge());
 
        p_cms_missA -= p_cms;
        p_cms_missB -= p_cms;
        et[0] += p_cms.Pt();
        et[1] += p_cms.Pt();
      }
 
      // ECLCluster
      //
      std::vector<const Particle*> roePhotons = roe->getPhotons(maskName);
 
      for (const Particle* photon : roePhotons) {
 
        if (photon->getECLClusterEHypothesisBit() == ECLCluster::EHypothesisBit::c_nPhotons) {
 
          ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * photon->get4Vector();
          p3_cms_all.push_back(p_cms.Vect());
          p3_cms_roe.push_back(p_cms.Vect());
 
          p3_cms_q_roe.emplace_back(p_cms.Vect(), photon->getCharge());
 
          p_cms_missA -= p_cms;
          p_cms_missB -= p_cms;
          et[0] += p_cms.Pt();
          et[1] += p_cms.Pt();
        }
      }
 
      // Thrust variables
      thrustB = Thrust::calculateThrust(p3_cms_sigB);
      thrustO = Thrust::calculateThrust(p3_cms_roe);
      thrustBm = thrustB.R();
      thrustOm = thrustO.R();
      cosTBTO  = fabs(thrustB.Unit().Dot(thrustO.Unit()));
      cosTBz   = fabs(cos(thrustB.Theta()));
 
      // Cleo Cones
      CleoCones cc(p3_cms_all, p3_cms_roe, thrustB, true, true);
      cleoConesAll = cc.cleo_cone_with_all();
      cleoConesROE = cc.cleo_cone_with_roe();
 
      // Fox-Wolfram Moments: Uses all final-state tracks (= sigB + ROE)
      FoxWolfram FW(p3_cms_all);
      FW.calculateBasicMoments();
      R2 = FW.getR(2);
 
      // KSFW moments
      ROOT::Math::PxPyPzEVector p_cms_B = T.rotateLabToCms() * particle->get4Vector();
      double Hso0_max(2 * (2 * BeamEnergy - p_cms_B.E()));
      KsfwMoments KsfwM(Hso0_max,
                        p3_cms_q_sigA,
                        p3_cms_q_sigB,
                        p3_cms_q_roe,
                        p_cms_missA,
                        p_cms_missB,
                        et);
      // use_finalstate_for_sig == 0
      KsfwM.usefinal(0);
      ksfwFS0.push_back(KsfwM.mm2());
      ksfwFS0.push_back(KsfwM.et());
      ksfwFS0.push_back(KsfwM.Hso(0, 0));
      ksfwFS0.push_back(KsfwM.Hso(0, 1));
      ksfwFS0.push_back(KsfwM.Hso(0, 2));
      ksfwFS0.push_back(KsfwM.Hso(0, 3));
      ksfwFS0.push_back(KsfwM.Hso(0, 4));
      ksfwFS0.push_back(KsfwM.Hso(1, 0));
      ksfwFS0.push_back(KsfwM.Hso(1, 2));
      ksfwFS0.push_back(KsfwM.Hso(1, 4));
      ksfwFS0.push_back(KsfwM.Hso(2, 0));
      ksfwFS0.push_back(KsfwM.Hso(2, 2));
      ksfwFS0.push_back(KsfwM.Hso(2, 4));
      ksfwFS0.push_back(KsfwM.Hoo(0));
      ksfwFS0.push_back(KsfwM.Hoo(1));
      ksfwFS0.push_back(KsfwM.Hoo(2));
      ksfwFS0.push_back(KsfwM.Hoo(3));
      ksfwFS0.push_back(KsfwM.Hoo(4));
      // use_finalstate_for_sig == 1
      KsfwM.usefinal(1);
      ksfwFS1.push_back(KsfwM.mm2());
      ksfwFS1.push_back(KsfwM.et());
      ksfwFS1.push_back(KsfwM.Hso(0, 0));
      ksfwFS1.push_back(KsfwM.Hso(0, 1));
      ksfwFS1.push_back(KsfwM.Hso(0, 2));
      ksfwFS1.push_back(KsfwM.Hso(0, 3));
      ksfwFS1.push_back(KsfwM.Hso(0, 4));
      ksfwFS1.push_back(KsfwM.Hso(1, 0));
      ksfwFS1.push_back(KsfwM.Hso(1, 2));
      ksfwFS1.push_back(KsfwM.Hso(1, 4));
      ksfwFS1.push_back(KsfwM.Hso(2, 0));
      ksfwFS1.push_back(KsfwM.Hso(2, 2));
      ksfwFS1.push_back(KsfwM.Hso(2, 4));
      ksfwFS1.push_back(KsfwM.Hoo(0));
      ksfwFS1.push_back(KsfwM.Hoo(1));
      ksfwFS1.push_back(KsfwM.Hoo(2));
      ksfwFS1.push_back(KsfwM.Hoo(3));
      ksfwFS1.push_back(KsfwM.Hoo(4));
 
      // TODO: The following is from the original belle ksfwmoments.cc module.
      //       Not sure if necessary here (i.e., will we be using rooksfw in belle II in the same way?).
      //       printf("rooksfw::rooksfw: mm2=%f et=%f hoo2=%f hso02=%f\n",
      //       m_mm2[0], et[0], m_Hoo[0][2], m_Hso[0][0][2]);
    }
 
    // Fill ContinuumSuppression with content
    qqVars->addThrustB(thrustB);
    qqVars->addThrustO(thrustO);
    qqVars->addThrustBm(thrustBm);
    qqVars->addThrustOm(thrustOm);
    qqVars->addCosTBTO(cosTBTO);
    qqVars->addCosTBz(cosTBz);
    qqVars->addR2(R2);
    qqVars->addKsfwFS0(ksfwFS0);
    qqVars->addKsfwFS1(ksfwFS1);
    qqVars->addCleoConesALL(cleoConesAll);
    qqVars->addCleoConesROE(cleoConesROE);
  }

CleoCones()

CleoCones	(	const std::vector< ROOT::Math::XYZVector > &	p3_cms_all,
		const std::vector< ROOT::Math::XYZVector > &	p3_cms_roe,
		const ROOT::Math::XYZVector &	thrustB,
		bool	calc_CleoCones_with_all,
		bool	calc_CleoCones_with_roe
	)

Constructor.

Definition at line 18 of file CleoCones.cc.

  {
    m_cleo_cone_with_all.clear();
    m_cleo_cone_with_roe.clear();
 
    // ----------------------------------------------------------------------
    // Calculate momentum flow in 9 cones for all particles in event
    // ----------------------------------------------------------------------
    if (calc_CleoCones_with_all == true) {
      for (int i = 1; i <= 9; i++) {
        float momentum_flow_all = 0;
        for (auto& iter0 : p3_cms_all) {
 
          /* Use the following intervals
             0*10<= <1*10  0- 10   170-180  180-1*10< <=180-0*10
             1*10<= <2*10 10- 20   160-170  180-2*10< <=180-1*10
             2*10<= <3*10 20- 30   150-160  180-3*10< <=180-2*10
             3*10<= <4*10 30- 40   140-150  180-4*10< <=180-3*10
             4*10<= <5*10 40- 50   130-140  180-5*10< <=180-4*10
             5*10<= <6*10 50- 60   120-130  180-6*10< <=180-5*10
             6*10<= <7*10 60- 70   110-120  180-7*10< <=180-6*10
             7*10<= <8*10 70- 80   100-110  180-8*10< <=180-7*10
             8*10<= <9*10 80- 90    90-100  180-9*10< <=180-8*10
             ==90 */
 
          float angle = ((180 * acos(thrustB.Unit().Dot(iter0.Unit()))) / M_PI);
          if (((((i - 1) * 10) <= angle) && (angle < (i * 10))) || (((180 - (i * 10)) < angle) && (angle <= (180 - ((i - 1) * 10))))) {
            momentum_flow_all += iter0.R();
            // B2DEBUG(19, "interval " << ((i-1)*10) << " to " << (i*10) << " and " << (180-(i*10)) << " to " << (180-((i-1)*10)) << " has value " << (180*(thrustB.angle(*iter0)))/M_PI << ", momentum flow is " << momentum_flow );
          }
          if ((i == 9) && (angle == 90)) {
            momentum_flow_all += iter0.R();
          }
        }
        m_cleo_cone_with_all.push_back(momentum_flow_all);
      }
    }
 
    // ----------------------------------------------------------------------
    // Calculate momentum flow in 9 cones for all particles in rest of event
    // ----------------------------------------------------------------------
    if (calc_CleoCones_with_roe == true) {
      for (int i = 1; i <= 9; i++) {
        float momentum_flow_roe = 0;
        for (auto& iter1 : p3_cms_roe) {
          float angle = ((180 * acos(thrustB.Unit().Dot(iter1.Unit()))) / M_PI);
          if (((((i - 1) * 10) <= angle) && (angle < (i * 10))) || (((180 - (i * 10)) < angle) && (angle <= (180 - ((i - 1) * 10))))) {
            momentum_flow_roe += iter1.R();
          }
          if ((i == 9) && (angle == 90)) {
            momentum_flow_roe += iter1.R();
          }
        }
        m_cleo_cone_with_roe.push_back(momentum_flow_roe);
      }
    }
  }

create()

void create

(

)

Create the list objects.

Warning

This function should be called in Every event

Definition at line 40 of file ParticleListHelper.cc.

  {
    m_list.create();
    m_list->initialize(m_pdg, m_list.getName());
    if (m_antiList) {
      m_antiList->create();
      (*m_antiList)->initialize(-m_pdg, m_antiList->getName());
      m_list->bindAntiParticleList(**m_antiList);
    }
  }

createCurrentParticle()

Particle createCurrentParticle ( ) const

private

Create current particle object.

Definition at line 439 of file ParticleCombiner.cc.

  {
    double px = 0;
    double py = 0;
    double pz = 0;
    double E = 0;
    for (const Particle* d : m_particles) {
      px += d->getPx();
      py += d->getPy();
      pz += d->getPz();
      E += d->getEnergy();
    }
    const ROOT::Math::PxPyPzEVector vec(px, py, pz, E);
 
    switch (m_iParticleType) {
      case 0: return Particle(vec, m_pdgCode, m_isSelfConjugated ? Particle::c_Unflavored : Particle::c_Flavored, m_indices,
                                m_properties, m_daughterProperties,
                                m_particleArray.getPtr());
      case 1: return Particle(vec, -m_pdgCode, m_isSelfConjugated ? Particle::c_Unflavored : Particle::c_Flavored, m_indices,
                                m_properties, m_daughterProperties,
                                m_particleArray.getPtr());
      case 2: return Particle(vec, m_pdgCode, Particle::c_Unflavored, m_indices,
                                m_properties, m_daughterProperties,
                                m_particleArray.getPtr());
      default: B2FATAL("You called getCurrentParticle although loadNext should have returned false!");
    }
 
    return Particle(); // This should never happen
  }

currentCombinationHasDifferentSources()

bool currentCombinationHasDifferentSources ( )

private

Check that all FS particles of a combination differ.

The comparison is made at the MDST objects level. If for example a kaon and a pion Particle's are created from the same MSDT Track object, then these two Particles have the same source and the function will return false.

Returns

true if all FS particles of a combination differ

Definition at line 474 of file ParticleCombiner.cc.

  {
    std::vector<Particle*> stack = m_particles;
    static std::vector<int> sources; // stack for particle sources
    sources.clear();
 
    // recursively check all daughters and daughters of daughters
    while (!stack.empty()) {
      Particle* p = stack.back();
      stack.pop_back();
      const std::vector<int>& daughters = p->getDaughterIndices();
 
      if (daughters.empty()) {
        int source = p->getMdstSource();
        for (int i : sources) {
          if (source == i) return false;
        }
        sources.push_back(source);
      } else {
        for (int daughter : daughters) stack.push_back(m_particleArray[daughter]);
      }
    }
    return true;
  }

currentCombinationIsECLCRUnique()

bool currentCombinationIsECLCRUnique ( )

private

Check that: if the current combination has at least two particles from an ECL source, then they are from different connected regions or from the same connected region but have the same hypothesis.

Returns

true if indices not found in the stack; if true indices pushed to stack

Definition at line 522 of file ParticleCombiner.cc.

  {
    unsigned nECLSource = 0;
    std::vector<Particle*> stack = m_particles;
    static std::vector<int> connectedregions;
    static std::vector<ECLCluster::EHypothesisBit> hypotheses;
    connectedregions.clear();
    hypotheses.clear();
 
    // recursively check all daughters and daughters of daughters
    while (!stack.empty()) {
      Particle* p = stack.back();
      stack.pop_back();
      const std::vector<int>& daughters = p->getDaughterIndices();
 
      if (daughters.empty()) {
        // Only test if the particle was created from an ECLCluster at source.
        // This CAN CHANGE if we change the cluster <--> track matching,
        // (currently we match nPhotons clusters to all track type particles,
        // i.e. electrons, pions, kaons etc. We might gain by matching
        // neutralHadron hypothesis clusters to kaons, for example).
        //
        // In the above case one would need to extend the check to ALL
        // particles with an associated ECLCluster. Then replace the following
        // active line with two lines...
        //
        // auto cluster = p->getECLCluster();
        // if (cluster) { // then do stuff
        if (p->getParticleSource() == Particle::EParticleSourceObject::c_ECLCluster) {
          nECLSource++;
          auto* cluster = p->getECLCluster();
          connectedregions.push_back(cluster->getConnectedRegionId());
          hypotheses.push_back(p->getECLClusterEHypothesisBit());
        }
      } else {
        for (int daughter : daughters) stack.push_back(m_particleArray[daughter]);
      }
    }
 
    // less than two particles from an ECL source is fine
    // (unless cluster <--> track matching changes)
    if (nECLSource < 2) return true;
 
    // yes this is a nested for loop but it's fast, we promise
    for (unsigned icr = 0; icr < connectedregions.size(); ++icr)
      for (unsigned jcr = icr + 1; jcr < connectedregions.size(); ++jcr)
        if (connectedregions[icr] == connectedregions[jcr])
          if (hypotheses[icr] != hypotheses[jcr]) return false;
 
    return true;
  }

currentCombinationIsUnique()

bool currentCombinationIsUnique ( )

private

Check that the combination is unique.

Especially in the case of reconstructing self conjugated decays we get combinations like M -> A B and M -> B A. These two particles are the same and hence only one combination of the two should be kept. This function takes care of this. It keeps track of all combinations that were already accepted (ordered set of unique IDs of daughter particles) and if the current combination is already found in the set it is discarded. The unique ID of daughter particles are either their StoreArray indices (if input lists do not collide) or specially set unique IDs during initialization phase (if input lists collide).

Returns

true if indices not found in the stack; if true indices pushed to stack

Definition at line 500 of file ParticleCombiner.cc.

  {
    std::set<int> indexSet;
    if (not m_inputListsCollide)
      indexSet.insert(m_indices.begin(), m_indices.end());
    else
      for (unsigned int i = 0; i < m_numberOfLists; i++)
        indexSet.insert(m_indicesToUniqueIDs.at(m_indices[i]));
 
    bool elementInserted = m_usedCombinations.insert(indexSet).second;
    return elementInserted;
  }

fillIndicesToUniqueIDMap() [1/2]

void fillIndicesToUniqueIDMap ( const std::vector< int > &  listA, const std::vector< int > &  listB, int &  uniqueID  )

private

Assigns unique IDs to all particles in list A, which do not have the unique ID already assigned.

The same unique ID is assigned to copies of particles from list A found in the list B. This function has to be executed first.

Definition at line 287 of file ParticleCombiner.cc.

  {
    const Particle* A, *B;
    bool copies = false;
    for (int i : listA) {
      bool aIsAlreadyIn = m_indicesToUniqueIDs.count(i) ? true : false;
 
      if (not aIsAlreadyIn)
        m_indicesToUniqueIDs[ i ] = uniqueID++;
 
      for (int j : listB) {
        bool bIsAlreadyIn = m_indicesToUniqueIDs.count(j) ? true : false;
 
        if (bIsAlreadyIn)
          continue;
 
        // are these two particles copies
        A = m_particleArray[ i ];
        B = m_particleArray[ j ];
        copies = B->isCopyOf(A);
 
        if (copies)
          m_indicesToUniqueIDs[ j ] = m_indicesToUniqueIDs[ i ];
      }
    }
  }

fillIndicesToUniqueIDMap() [2/2]

void fillIndicesToUniqueIDMap ( const std::vector< int > &  listA, int &  uniqueID  )

private

Assigns unique IDs to all particles in list A, which do not have the unique ID already assigned.

Definition at line 315 of file ParticleCombiner.cc.

  {
    for (int i : listA) {
      bool aIsAlreadyIn = m_indicesToUniqueIDs.count(i) ? true : false;
 
      if (not aIsAlreadyIn)
        m_indicesToUniqueIDs[ i ] = uniqueID++;
    }
  }

find_decay()

bool find_decay

(

const DecayNode &

to_find

)

const

Check if the decay node contains the given decay tree.

Parameters

to_find

DecayNode object describing the decay

Definition at line 19 of file DecayNode.cc.

  {
    if (to_find == (*this))
      return true;
    for (const auto& node : daughters) {
      if (node.find_decay(to_find))
        return true;
    }
    return false;
  }

getCurrentIndices() [1/2]

const std::vector< unsigned int > & getCurrentIndices

(

)

const

Returns theindices of the current loaded combination.

Definition at line 75 of file ParticleCombiner.cc.

  {
    return indices;
  }

getCurrentIndices() [2/2]

const std::vector< ParticleList::EParticleType > & getCurrentIndices

(

)

const

Returns the type of the sublist of the current loaded combination.

Definition at line 105 of file ParticleCombiner.cc.

  {
    return m_types;
  }

getCurrentParticle()

Particle getCurrentParticle

(

)

const

Returns the particle.

Definition at line 469 of file ParticleCombiner.cc.

  {
    return m_current_particle;
  }

getUniqueID()

int getUniqueID

(

int

index

)

const

Returns the unique ID assigned to Particle with given index from the IndicesToUniqueID map.

If the Particle's index is not found in the map -1 is returned instead. If the map was never filled (inputListsCollide is false) 0 is returned.

Needed only for tests.

Definition at line 574 of file ParticleCombiner.cc.

  {
    if (not m_inputListsCollide)
      return 0;
 
    if (not m_indicesToUniqueIDs.count(index))
      return -1;
 
    return m_indicesToUniqueIDs.at(index);
  }

init() [1/3]

void init

(

)

Initialises the generator to produce the given type of sublist.

Definition at line 221 of file ParticleCombiner.cc.

  {
    m_iParticleType = 0;
 
    m_particleIndexGenerator.init(std::vector<unsigned int> {}); // ParticleIndexGenerator will be initialised on first call
    m_listIndexGenerator.init(m_numberOfLists); // ListIndexGenerator must be initialised here!
    m_usedCombinations.clear();
    m_indices.resize(m_numberOfLists);
    m_particles.resize(m_numberOfLists);
 
    if (m_inputListsCollide)
      initIndicesToUniqueIDMap();
  }

init() [2/3]

void init

(

const std::vector< unsigned int > &

_sizes

)

Initialises the generator to produce combinations with the given sizes of each particle list.

Parameters

_sizes

the sizes of the particle lists to combine

Definition at line 28 of file ParticleCombiner.cc.

  {
 
    m_numberOfLists = _sizes.size();
    sizes = _sizes;
 
    m_iCombination = 0;
 
    indices.resize(m_numberOfLists);
    for (unsigned int i = 0; i < m_numberOfLists; ++i) {
      indices[i] = 0;
    }
 
    m_nCombinations = 1;
    for (unsigned int i = 0; i < m_numberOfLists; ++i)
      m_nCombinations *= sizes[i];
 
    if (m_numberOfLists == 0) {
      m_nCombinations = 0;
    }
 
  }

init() [3/3]

void init

(

unsigned int

_numberOfLists

)

Initialises the generator to produce the given type of sublist.

Parameters

_numberOfLists

Number of Particle Lists which shall be combined

Definition at line 80 of file ParticleCombiner.cc.

  {
    m_numberOfLists = _numberOfLists;
    m_types.resize(m_numberOfLists);
 
    m_iCombination = 0;
    m_nCombinations = (1 << m_numberOfLists) - 1;
 
  }

initIndicesToUniqueIDMap()

void initIndicesToUniqueIDMap ( )

private

In the case input daughter particle lists collide (two or more lists contain copies of Particles) the Particle's Store Array index can not be longer used as its unique identifier, which is needed to check for uniqueness of accepted combinations.

Instead unique identifier is created for all particles in the input particle lists and use those when checking for uniqueness of current combination.

Definition at line 235 of file ParticleCombiner.cc.

  {
    m_indicesToUniqueIDs.clear();
 
    unsigned inputParticlesCount = 0;
    for (unsigned int i = 0; i < m_numberOfLists; ++i)
      inputParticlesCount +=  m_plists[i]->getListSize();
 
    m_indicesToUniqueIDs.reserve(inputParticlesCount);
 
    int uniqueID = 1;
 
    for (auto& collidingList : m_collidingLists) {
      StoreObjPtr<ParticleList> listA =  m_plists[collidingList.first];
      StoreObjPtr<ParticleList> listB =  m_plists[collidingList.second];
 
      bool sameSign = (listA->getPDGCode() == listB->getPDGCode());
 
      // if sameSign == true then
      // 1. compare FS to FS particles in lists A and B
      // 2. compare anti-FS to anti-FS particles in lists A and B
      // else
      // 1. compare FS to anti-FS particles in lists A and B
      // 2. compare anti-FS to FS particles in lists A and B
      // and in either case compare
      // 3. compare SC to SC particles in lists A and B
      if (sameSign) {
        fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false),
                                 listB->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
        fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true),
                                 listB->getList(ParticleList::c_FlavorSpecificParticle, true),  uniqueID);
      } else {
        fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false),
                                 listB->getList(ParticleList::c_FlavorSpecificParticle, true),  uniqueID);
        fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true),
                                 listB->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
      }
 
      fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_SelfConjugatedParticle),
                               listB->getList(ParticleList::c_SelfConjugatedParticle),  uniqueID);
    }
 
    // assign unique ids to others as well
    for (unsigned i = 0; i < m_numberOfLists; i++) {
      StoreObjPtr<ParticleList> listA =  m_plists[i];
 
      fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true),  uniqueID);
      fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
      fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_SelfConjugatedParticle), uniqueID);
    }
  }

inputListsCollide()

bool inputListsCollide

(

const std::pair< unsigned, unsigned > &

pair

)

const

True if the pair of input lists collide.

Needed only for tests.

Definition at line 513 of file ParticleCombiner.cc.

  {
    for (const auto& collidingList : m_collidingLists)
      if (pair == collidingList)
        return true;
 
    return false;
  }

KsfwMoments()

KsfwMoments	(	double	Hso0_max,
		std::vector< std::pair< ROOT::Math::XYZVector, int > >	p3_cms_q_sigA,
		std::vector< std::pair< ROOT::Math::XYZVector, int > >	p3_cms_q_sigB,
		std::vector< std::pair< ROOT::Math::XYZVector, int > >	p3_cms_q_roe,
		const ROOT::Math::PxPyPzEVector &	p_cms_missA,
		const ROOT::Math::PxPyPzEVector &	p_cms_missB,
		const double	et[2]
	)

Constructor.

Definition at line 48 of file KsfwMoments.cc.

  {
    // Private member needs to be initialized. Here it is initialized to -1 (illegal value).
    m_uf = -1;
 
    m_et[0] = et[0];
    m_et[1] = et[1];
 
    m_mm2[0] = p_cms_missA.E() > 0
               ? p_cms_missA.mag2()
               : -p_cms_missA.E() * p_cms_missA.E() - p_cms_missA.P2();
    m_mm2[1] = p_cms_missB.E() > 0
               ? p_cms_missB.mag2()
               : -p_cms_missB.E() * p_cms_missB.E() - p_cms_missB.P2();
 
    //========================
    // Calculate discriminants
    //========================
    std::vector<std::pair<ROOT::Math::XYZVector, int>>::iterator pqi, pqj;
 
    // Calculate Hso components
    for (int i = 0; i < 3; i++) {
      for (int k = 0; k < 5; k++) {
        m_Hso[0][i][k] = m_Hso[1][i][k] = 0;
      }
    }
 
    // Signal A (use_finalstate_for_sig == 0)
    for (pqi = p3_cms_q_sigA.begin(); pqi != p3_cms_q_sigA.end(); ++pqi) {
      const double pi_mag((pqi->first).R());
      for (pqj = p3_cms_q_roe.begin(); pqj != p3_cms_q_roe.end(); ++pqj) {
        const double pj_mag((pqj->first).R());
        const double ij_cos((pqi->first).Dot(pqj->first) / pi_mag / pj_mag);
        const int c_or_n(0 == (pqj->second) ? 1 : 0);  // 0: charged 1: neutral
        for (int k = 0; k < 5; k++) {
          m_Hso[0][c_or_n][k] += (k % 2)
                                 ? (pqi->second) * (pqj->second) * pj_mag * legendre(ij_cos, k)
                                 : pj_mag * legendre(ij_cos, k);
        }
      }
      const double p_miss_mag(p_cms_missA.P());
      const double i_miss_cos((pqi->first).Dot(p_cms_missA.Vect()) / pi_mag / p_miss_mag);
      for (int k = 0; k < 5; k++) {
        m_Hso[0][2][k] += (k % 2) ? 0 : p_miss_mag * legendre(i_miss_cos, k);
      }
    }
 
    // Signal B (use_finalstate_for_sig == 1)
    for (pqi = p3_cms_q_sigB.begin(); pqi != p3_cms_q_sigB.end(); ++pqi) {
      const double pi_mag((pqi->first).R());
      for (pqj = p3_cms_q_roe.begin(); pqj != p3_cms_q_roe.end(); ++pqj) {
        const double pj_mag((pqj->first).R());
        const double ij_cos((pqi->first).Dot(pqj->first) / pi_mag / pj_mag);
        const int c_or_n(0 == (pqj->second) ? 1 : 0);  // 0: charged 1: neutral
        for (int k = 0; k < 5; k++) {
          m_Hso[1][c_or_n][k] += (k % 2)
                                 ? (pqi->second) * (pqj->second) * pj_mag * legendre(ij_cos, k)
                                 : pj_mag * legendre(ij_cos, k);
        }
      }
      const double p_miss_mag(p_cms_missB.P());
      const double i_miss_cos((pqi->first).Dot(p_cms_missB.Vect()) / pi_mag / p_miss_mag);
      for (int k = 0; k < 5; k++) {
        m_Hso[1][2][k] += (k % 2) ? 0 : p_miss_mag * legendre(i_miss_cos, k);
      }
    }
 
    // Add missing to the lists
    std::vector<std::pair<ROOT::Math::XYZVector, int>> p3_cms_q_roeA(p3_cms_q_roe), p3_cms_q_roeB(p3_cms_q_roe);
    p3_cms_q_roeA.emplace_back(p_cms_missA.Vect(), 0);
    p3_cms_q_roeB.emplace_back(p_cms_missB.Vect(), 0);
 
    // Calculate Hoo components
    for (int k = 0; k < 5; k++) {
      m_Hoo[0][k] = m_Hoo[1][k] = 0;
    }
    for (pqi = p3_cms_q_roeA.begin(); pqi != p3_cms_q_roeA.end(); ++pqi) {
      const double pi_mag((pqi->first).R());
      for (pqj = p3_cms_q_roeA.begin(); pqj != pqi; ++pqj) {
        const double pj_mag((pqj->first).R());
        const double ij_cos((pqi->first).Dot(pqj->first) / pi_mag / pj_mag);
        for (int k = 0; k < 5; k++) {
          m_Hoo[0][k] += (k % 2)
                         ? (pqi->second) * (pqj->second) * pi_mag * pj_mag * legendre(ij_cos, k)
                         : pi_mag * pj_mag * legendre(ij_cos, k);
        }
      }
    }
    for (pqi = p3_cms_q_roeB.begin(); pqi != p3_cms_q_roeB.end(); ++pqi) {
      const double pi_mag((pqi->first).R());
      for (pqj = p3_cms_q_roeB.begin(); pqj != pqi; ++pqj) {
        const double pj_mag((pqj->first).R());
        const double ij_cos((pqi->first).Dot(pqj->first) / pi_mag / pj_mag);
        for (int k = 0; k < 5; k++) {
          m_Hoo[1][k] += (k % 2)
                         ? (pqi->second) * (pqj->second) * pi_mag * pj_mag * legendre(ij_cos, k)
                         : pi_mag * pj_mag * legendre(ij_cos, k);
        }
      }
    }
 
    // Normalize so that it does not dependent on delta_e
    for (int k = 0; k < 5; k++) {
      for (int j = 0; j < ((k % 2) ? 1 : 3); j++) {
        m_Hso[0][j][k] /= Hso0_max;
        m_Hso[1][j][k] /= Hso0_max;
      }
      m_Hoo[0][k] /= (Hso0_max * Hso0_max);
      m_Hoo[1][k] /= (Hso0_max * Hso0_max);
    }
 
  }

legendre()

double legendre ( const double  z, const int  i  )

inline

Legendre polynomials.

Definition at line 33 of file KsfwMoments.cc.

  {
    switch (i) {
      case 0:  return 1.;
      case 1:  return z;
      case 2:  return 1.5 * z * z - 0.5;
      case 3:  return z * (2.5 * z * z - 1.5);
      case 4:  return (4.375 * z * z * z * z - 3.75 * z * z + 0.375);
      default: return 0;
    }
  }

loadNext() [1/3]

bool loadNext

(

)

Loads the next combination.

Returns false if there is no next combination

Definition at line 51 of file ParticleCombiner.cc.

  {
 
    if (m_iCombination == m_nCombinations) {
      return false;
    }
 
    if (m_iCombination > 0) {
 
      // TF SC this does not yet account for double counting so will produce:
      // { 000, 100, 200, ... 010, 110, .... } even if the first and second
      // place are the same particle list
      for (unsigned int i = 0; i < m_numberOfLists; i++) {
        indices[i]++;
        if (indices[i] < sizes[i]) break;
        indices[i] = 0;
      }
 
    }
 
    m_iCombination++;
    return true;
  }

loadNext() [2/3]

bool loadNext

(

)

Loads the next combination.

Returns false if there is no next combination

Definition at line 90 of file ParticleCombiner.cc.

  {
 
    if (m_iCombination == m_nCombinations)
      return false;
 
    for (unsigned int i = 0; i < m_numberOfLists; ++i) {
      bool useSelfConjugatedDaughterList = (m_iCombination & (1 << i));
      m_types[i] = useSelfConjugatedDaughterList ? ParticleList::c_SelfConjugatedParticle : ParticleList::c_FlavorSpecificParticle;
    }
 
    ++m_iCombination;
    return true;
  }

loadNext() [3/3]

bool loadNext

(

bool

loadAntiParticle = true

)

Loads the next combination.

Returns false if there is no next combination

Three cases are distinguished: First, particles matching the flavor specified in the decay string are used to form combinations. Secondly, the anti-particles of flavored particles are used, but only if requested. Lastly, self-conjugated particles are handled specifically.

Definition at line 326 of file ParticleCombiner.cc.

  {
 
    bool loadedNext = false;
    while (true) {
      if (m_iParticleType == 0) {
        loadedNext = loadNextParticle(false);
      } else if (m_iParticleType == 1 and loadAntiParticle) {
        loadedNext = loadNextParticle(true);
      } else if (m_iParticleType == 2) {
        loadedNext = loadNextSelfConjugatedParticle();
      } else {
        return false;
      }
 
      if (loadedNext) return true;
      else ++m_iParticleType;
 
      if (m_iParticleType == 2) {
        std::vector<unsigned int> sizes(m_numberOfLists);
        for (unsigned int i = 0; i < m_numberOfLists; ++i) {
          sizes[i] = m_plists[i]->getList(ParticleList::c_SelfConjugatedParticle, false).size();
        }
        m_particleIndexGenerator.init(sizes);
      } else {
        m_listIndexGenerator.init(m_numberOfLists);
      }
 
    }
  }

loadNextParticle()

bool loadNextParticle ( bool  useAntiParticle )

private

Loads the next combination.

Returns false if there is no next combination

Definition at line 363 of file ParticleCombiner.cc.

  {
    while (true) {
 
      // Load next index combination if available
      if (m_particleIndexGenerator.loadNext()) {
 
        const auto& m_types = m_listIndexGenerator.getCurrentIndices();
        const auto& indices = m_particleIndexGenerator.getCurrentIndices();
 
        for (unsigned int i = 0; i < m_numberOfLists; i++) {
          const auto& list = m_plists[i]->getList(m_types[i], m_types[i] == ParticleList::c_FlavorSpecificParticle ? useAntiParticle : false);
          m_indices[i] =  list[ indices[i] ];
          m_particles[i] = m_particleArray[ m_indices[i] ];
        }
 
        if (not currentCombinationHasDifferentSources()) continue;
 
        m_current_particle = createCurrentParticle();
        if (!m_cut->check(&m_current_particle))
          continue;
 
        if (not currentCombinationIsUnique()) continue;
 
        if (not currentCombinationIsECLCRUnique()) continue;
 
        return true;
      }
 
      // Load next list combination if available and reset indexCombiner
      if (m_listIndexGenerator.loadNext()) {
        const auto& m_types = m_listIndexGenerator.getCurrentIndices();
        std::vector<unsigned int> sizes(m_numberOfLists);
        for (unsigned int i = 0; i < m_numberOfLists; ++i) {
          sizes[i] = m_plists[i]->getList(m_types[i], m_types[i] == ParticleList::c_FlavorSpecificParticle ? useAntiParticle : false).size();
        }
        m_particleIndexGenerator.init(sizes);
        continue;
      }
      return false;
    }
 
  }

loadNextSelfConjugatedParticle()

bool loadNextSelfConjugatedParticle ( )

private

Loads the next combination.

Returns false if there is no next combination

Definition at line 407 of file ParticleCombiner.cc.

  {
    while (true) {
 
      // Load next index combination if available
      if (m_particleIndexGenerator.loadNext()) {
 
        const auto& indices = m_particleIndexGenerator.getCurrentIndices();
 
        for (unsigned int i = 0; i < m_numberOfLists; i++) {
          m_indices[i] = m_plists[i]->getList(ParticleList::c_SelfConjugatedParticle, false) [ indices[i] ];
          m_particles[i] = m_particleArray[ m_indices[i] ];
        }
 
        if (not currentCombinationHasDifferentSources()) continue;
 
        m_current_particle = createCurrentParticle();
        if (!m_cut->check(&m_current_particle))
          continue;
 
        if (not currentCombinationIsUnique()) continue;
 
        if (not currentCombinationIsECLCRUnique()) continue;
 
        return true;
      }
 
      return false;
    }
 
  }

operator!=()

bool operator!=	(	const DecayNode &	node1,
		const DecayNode &	node2
	)

Not equal: See operator==.

Definition at line 65 of file DecayNode.cc.

  {
    return not(node1 == node2);
  }

operator==()

bool operator==	(	const DecayNode &	node1,
		const DecayNode &	node2
	)

Compare two Decay Nodes: They are equal if All daughter decay nodes are equal or one of the daughter lists is empty (which means "inclusive")

In particular the order of the daughter particles matter! 423 (--> 421 13) is not the same as 423 (--> 13 421)

Parameters

node1	first node
node2	second node

Definition at line 48 of file DecayNode.cc.

  {
    if (node1.pdg == node2.pdg) {
      if (node1.daughters.size() > 0 and node2.daughters.size() > 0) {
        if (node1.daughters.size() != node2.daughters.size())
          return false;
        for (unsigned int i = 0; i < node1.daughters.size(); ++i) {
          if (node1.daughters[i] != node2.daughters[i])
            if (node1.daughters[i] != 0 and node2.daughters[i] != 0)
              return false;
        }
      }
      return true;
    }
    return false;
  }

ParticleGenerator() [1/2]

ParticleGenerator ( const DecayDescriptor &  decaydescriptor, const std::string &  cutParameter = ""  )

explicit

Initialises the generator to produce the given type of sublist.

Parameters

decaydescriptor
cutParameter

Definition at line 167 of file ParticleCombiner.cc.

                                                                                                            : m_iParticleType(0),
    m_listIndexGenerator(),
    m_particleIndexGenerator()
  {
    bool valid = decaydescriptor.isInitOK();
    if (!valid)
      B2ERROR("Given decaydescriptor failed to initialized");
 
    m_properties = decaydescriptor.getProperty();
 
    // Mother particle
    const DecayDescriptorParticle* mother = decaydescriptor.getMother();
    m_pdgCode = mother->getPDGCode();
    m_properties |= mother->getProperty();
 
    // Daughters
    m_numberOfLists = decaydescriptor.getNDaughters();
    for (unsigned int i = 0; i < m_numberOfLists; ++i) {
      // Get the mother of the subdecaystring of the ith daughter
      // eg. "B -> [D -> K pi] [tau -> pi pi pi]". The 0th daughter is the
      // *decaystring* D -> K pi whose mother is the D.
      const DecayDescriptorParticle* daughter = decaydescriptor.getDaughter(i)->getMother();
      StoreObjPtr<ParticleList> list(daughter->getFullName());
      m_plists.push_back(list);
 
      int daughterProperty = daughter->getProperty();
      m_daughterProperties.push_back(daughterProperty);
    }
 
    m_cut = Variable::Cut::compile(cutParameter);
 
    m_isSelfConjugated = decaydescriptor.isSelfConjugated();
 
    // check if input lists can contain copies
    // remember (list_i, list_j) pairs, where list_i and list_j can contain copies
    m_inputListsCollide = false;
    m_collidingLists.clear();
    for (unsigned int i = 0; i < m_numberOfLists; ++i) {
      const DecayDescriptorParticle* daughter_i = decaydescriptor.getDaughter(i)->getMother();
      for (unsigned int j = i + 1; j < m_numberOfLists; ++j) {
        const DecayDescriptorParticle* daughter_j = decaydescriptor.getDaughter(j)->getMother();
 
        if (abs(daughter_i->getPDGCode()) != abs(daughter_j->getPDGCode()))
          continue;
 
        if (daughter_i->getLabel() != daughter_j->getLabel()) {
          m_inputListsCollide = true;
          m_collidingLists.emplace_back(i, j);
        }
      }
    }
 
  }

ParticleGenerator() [2/2]

ParticleGenerator ( const std::string &  decayString, const std::string &  cutParameter = ""  )

explicit

Initialises the generator to produce the given type of sublist.

Parameters

decayString
cutParameter

Definition at line 110 of file ParticleCombiner.cc.

                                                                                                  : m_iParticleType(0),
    m_listIndexGenerator(),
    m_particleIndexGenerator()
  {
 
    DecayDescriptor decaydescriptor;
    bool valid = decaydescriptor.init(decayString);
    if (!valid)
      B2ERROR("Invalid input DecayString: " << decayString);
 
    m_properties = decaydescriptor.getProperty();
 
    // Mother particle
    const DecayDescriptorParticle* mother = decaydescriptor.getMother();
    m_pdgCode = mother->getPDGCode();
    m_properties |= mother->getProperty();
 
    // Daughters
    m_numberOfLists = decaydescriptor.getNDaughters();
    for (unsigned int i = 0; i < m_numberOfLists; ++i) {
      // Get the mother of the subdecaystring of the ith daughter
      // eg. "B -> [D -> K pi] [tau -> pi pi pi]". The 0th daughter is the
      // *decaystring* D -> K pi whose mother is the D.
      const DecayDescriptorParticle* daughter = decaydescriptor.getDaughter(i)->getMother();
      StoreObjPtr<ParticleList> list(daughter->getFullName());
      m_plists.push_back(list);
 
      int daughterProperty = daughter->getProperty();
      m_daughterProperties.push_back(daughterProperty);
    }
 
    m_cut = Variable::Cut::compile(cutParameter);
 
    m_isSelfConjugated = decaydescriptor.isSelfConjugated();
 
    // check if input lists can contain copies
    // remember (list_i, list_j) pairs, where list_i and list_j can contain copies
    m_inputListsCollide = false;
    m_collidingLists.clear();
    for (unsigned int i = 0; i < m_numberOfLists; ++i) {
      const DecayDescriptorParticle* daughter_i = decaydescriptor.getDaughter(i)->getMother();
      for (unsigned int j = i + 1; j < m_numberOfLists; ++j) {
        const DecayDescriptorParticle* daughter_j = decaydescriptor.getDaughter(j)->getMother();
 
        if (abs(daughter_i->getPDGCode()) != abs(daughter_j->getPDGCode()))
          continue;
 
        if (daughter_i->getLabel() != daughter_j->getLabel()) {
          m_inputListsCollide = true;
          m_collidingLists.emplace_back(i, j);
        }
      }
    }
 
  }

print_node()

std::string print_node

(

unsigned int

indent = 0

)

const

Output a single node.

Used by to_string to convert tree into a string

Parameters

indent

-ion level

Definition at line 30 of file DecayNode.cc.

  {
    std::stringstream output;
 
    for (unsigned int i = 0; i < indent; ++i) {
      output << "    ";
    }
 
    output << std::to_string(pdg);
    output << std::endl;
 
    for (const auto& d : daughters) {
      output << d.print_node(indent + 1);
    }
 
    return output.str();
  }

registerList() [1/2]

void registerList	(	const DecayDescriptor &	decay,
		bool	save = true
	)

Register a list using a decay descriptor.

Warning

should be called in initialize

Parameters

decay	decay descriptor to use for the list name
save	if true the list will be saved in the output file

Definition at line 27 of file ParticleListHelper.cc.

  {
    auto listname = decay.getMother()->getFullName();
    m_pdg = decay.getMother()->getPDGCode();
    auto storeFlags = save ? DataStore::c_WriteOut : DataStore::c_DontWriteOut;
    m_particles.registerInDataStore(storeFlags);
    m_list.registerInDataStore(listname, storeFlags);
    auto antiListName = ParticleListName::antiParticleListName(listname);
    if (antiListName != listname) {
      m_antiList.emplace(antiListName);
      m_antiList->registerInDataStore(storeFlags);
    }
  }

registerList() [2/2]

void registerList	(	const std::string &	listname,
		bool	save = true
	)

Register a list by name.

Warning

should be called in initialize

Parameters

listname	full list name in the form '{particlename}:{label}'
save	if true the list will be saved in the output file

Definition at line 18 of file ParticleListHelper.cc.

  {
    DecayDescriptor decayDescriptor;
    bool valid = decayDescriptor.init(listname);
    if (!valid)
      throw std::runtime_error("Invalid ParticleList name: '" + listname + "' Should be EVTPDLNAME[:LABEL], e.g. B+ or B+:mylist.");
    ParticleListHelper::registerList(decayDescriptor, save);
  }

TEST_F() [1/9]

TEST_F	(	ChargedParticleIdentificatorTest	,
		TestDBRep
	)

Test correct storage of weightfiles in the database representation inner structure.

Definition at line 129 of file chargedParticleIdentificator.cc.

  {
 
    // Pick a random index in the clusterTheta, p, charge bins arrays.
    // Exclude underflow and overflow.
    std::random_device rd; // non-deterministic uniform int rand generator.
    std::uniform_int_distribution<int> binx_idx_distr(1, m_thetabins.size() - 1);
    std::uniform_int_distribution<int> biny_idx_distr(1, m_pbins.size() - 1);
    std::uniform_int_distribution<int> binz_idx_distr(1, m_chbins.size() - 1);
    int binx = binx_idx_distr(rd);
    int biny = biny_idx_distr(rd);
    int binz = binz_idx_distr(rd);
 
    // Pick each axis' bin centre as a test value for (clusterTheta, p, charge).
    auto theta = m_dbrep.getWeightCategories()->GetXaxis()->GetBinCenter(binx);
    auto p = m_dbrep.getWeightCategories()->GetYaxis()->GetBinCenter(biny);
    auto charge = m_dbrep.getWeightCategories()->GetZaxis()->GetBinCenter(binz);
 
    int jth, ip, kch;
    auto jik = m_dbrep.getMVAWeightIdx(theta, p, charge, jth, ip, kch);
 
    EXPECT_EQ(jth, binx);
    EXPECT_EQ(ip, biny);
    EXPECT_EQ(kch, binz);
 
    auto thisfname = m_basename
                     + "__clusterTheta__" + std::to_string(m_thetabins.at(jth - 1)) + "_" + std::to_string(m_thetabins.at(jth))
                     + "__p__" + std::to_string(m_pbins.at(ip - 1)) + "_" + std::to_string(m_pbins.at(ip))
                     + "__charge__" + std::to_string(charge);
    std::replace(thisfname.begin(), thisfname.end(), '.', '_');
    thisfname += ".xml";
 
    EXPECT_EQ(thisfname, m_dummyfiles.at(jik));
 
    auto matchitr = std::find(m_dummyfiles.begin(), m_dummyfiles.end(), thisfname);
    auto thisidx = std::distance(m_dummyfiles.begin(), matchitr);
 
    EXPECT_EQ(thisidx, jik);
 
  }

TEST_F() [2/9]

TEST_F	(	eventShapeCoreAlgorithmTest	,
		CleoCones
	)

Test the calculation of the CleoClones variables.

Definition at line 57 of file eventShapeCoreAlgorithm.cc.

  {
    const bool use_all = true;
    const bool use_roe = true;
    std::vector<ROOT::Math::XYZVector> momenta;
    std::vector<ROOT::Math::XYZVector> sig_side_momenta;
    std::vector<ROOT::Math::XYZVector> roe_side_momenta;
 
    // "Signal Side"
    sig_side_momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
    sig_side_momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
    sig_side_momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
    sig_side_momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
    sig_side_momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
 
    // "ROE Side"
    roe_side_momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
    roe_side_momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
    roe_side_momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
    roe_side_momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
    roe_side_momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
 
    // "All momenta"
    momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
    momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
    momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
    momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
    momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
    momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
    momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
    momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
    momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
    momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
 
 
    // Calculate thrust from "Signal Side"
    const ROOT::Math::XYZVector thrustB = Thrust::calculateThrust(sig_side_momenta);
 
    CleoCones myCleoCones(momenta, roe_side_momenta, thrustB, use_all, use_roe);
 
    const auto& cleo_cone_with_all = myCleoCones.cleo_cone_with_all();
    const auto& cleo_cone_with_roe = myCleoCones.cleo_cone_with_roe();
 
    EXPECT_FLOAT_EQ(0.823567, cleo_cone_with_all[0]);
    EXPECT_FLOAT_EQ(4.7405558, cleo_cone_with_all[1]);
    EXPECT_FLOAT_EQ(1.7517139, cleo_cone_with_all[2]);
    EXPECT_FLOAT_EQ(0.37677661, cleo_cone_with_all[3]);
    EXPECT_FLOAT_EQ(0.622467, cleo_cone_with_all[4]);
    EXPECT_FLOAT_EQ(0, cleo_cone_with_all[5]);
    EXPECT_FLOAT_EQ(0, cleo_cone_with_all[6]);
    EXPECT_FLOAT_EQ(0, cleo_cone_with_all[7]);
    EXPECT_FLOAT_EQ(0, cleo_cone_with_all[8]);
 
    EXPECT_FLOAT_EQ(0.823567, cleo_cone_with_roe[0]);
    EXPECT_FLOAT_EQ(1.9106253, cleo_cone_with_roe[1]);
    EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[2]);
    EXPECT_FLOAT_EQ(0.37677661, cleo_cone_with_roe[3]);
    EXPECT_FLOAT_EQ(0.622467, cleo_cone_with_roe[4]);
    EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[5]);
    EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[6]);
    EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[7]);
    EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[8]);
  }

TEST_F() [3/9]

TEST_F	(	eventShapeCoreAlgorithmTest	,
		FoxWolfram
	)

Test the calculation of the Fox-Wolfram moments.

Definition at line 122 of file eventShapeCoreAlgorithm.cc.

  {
    std::vector<ROOT::Math::XYZVector> momenta;
 
    momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
    momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
    momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
    momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
    momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
    momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
    momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
    momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
    momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
    momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
 
    FoxWolfram FW(momenta);
    FW.calculateBasicMoments();
    EXPECT_FLOAT_EQ(0.63011014, FW.getR(2));
 
    // Tests using the analytic values provided by Fox & Wolfram
    // for simple distributions
    // first: Back-to-back versors
    momenta.clear();
    momenta.emplace_back(1., 0., 0.);
    momenta.emplace_back(-1., 0., 0.);
 
    FW.setMomenta(momenta);
    FW.calculateBasicMoments();
 
    EXPECT_TRUE(TMath::Abs(FW.getR(0) - 1) < 0.001);
    EXPECT_TRUE(TMath::Abs(FW.getR(1)) < 0.001);
    EXPECT_TRUE(TMath::Abs(FW.getR(2) - 1) < 0.001);
    EXPECT_TRUE(TMath::Abs(FW.getR(3)) < 0.001);
    EXPECT_TRUE(TMath::Abs(FW.getR(4) - 1) < 0.001);
 
    // second: equatorial line
    momenta.clear();
    TRandom3 rnd;
    for (int i = 0; i < 10000; i++) {
      double phi = rnd.Uniform(0., 2 * TMath::Pi());
      momenta.emplace_back(TMath::Cos(phi), TMath::Sin(phi), 0.);
    }
    FW.setMomenta(momenta);
    FW.calculateBasicMoments();
 
    EXPECT_TRUE(TMath::Abs(FW.getR(0) - 1.) < 0.001);
    EXPECT_TRUE(TMath::Abs(FW.getR(1)) < 0.001);
    EXPECT_TRUE(TMath::Abs(FW.getR(2) - 0.25) < 0.001);
    EXPECT_TRUE(TMath::Abs(FW.getR(3)) < 0.001);
    EXPECT_TRUE(TMath::Abs(FW.getR(4) - 0.14) < 0.001);
 
  }

TEST_F() [4/9]

TEST_F	(	eventShapeCoreAlgorithmTest	,
		HarmonicMoments
	)

Test the calculation of the Harmonic moments.

Definition at line 178 of file eventShapeCoreAlgorithm.cc.

  {
 
    float dummySqrtS = 10.;
    std::vector<XYZVector> partMom;
 
    partMom.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
    partMom.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
    partMom.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
    partMom.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
    partMom.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
    partMom.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
    partMom.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
    partMom.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
    partMom.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
    partMom.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
 
    XYZVector axis(0., 0., 1.);
    // repeats the calculation
    double moment[9] = {0.};
    for (auto& p : partMom) {
      double pMag = p.R();
      double cTheta = p.Dot(axis) / pMag;
      for (short i = 0; i < 9; i++)
        moment[i] += pMag * boost::math::legendre_p(i, cTheta) / dummySqrtS;
    }
 
    HarmonicMoments MM(partMom, axis);
    MM.calculateAllMoments();
 
    for (short i = 0; i < 9; i++)
      EXPECT_FLOAT_EQ(moment[i], MM.getMoment(i, dummySqrtS));
  }

TEST_F() [5/9]

TEST_F	(	eventShapeCoreAlgorithmTest	,
		Sphericity
	)

Test the calculation of the Sphericity eigenvalues and eigenvectors.

Definition at line 216 of file eventShapeCoreAlgorithm.cc.

  {
    std::vector<XYZVector> partMom;
 
    partMom.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
    partMom.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
    partMom.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
    partMom.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
    partMom.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
    partMom.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
    partMom.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
    partMom.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
    partMom.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
    partMom.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
 
 
    SphericityEigenvalues Sph(partMom);
    EXPECT_FLOAT_EQ(0., Sph.getEigenvalue(0));
    EXPECT_FLOAT_EQ(0., Sph.getEigenvalue(1));
    EXPECT_FLOAT_EQ(0., Sph.getEigenvalue(2));
 
    Sph.calculateEigenvalues();
 
    EXPECT_FLOAT_EQ(0.87272972, Sph.getEigenvalue(0));
    EXPECT_FLOAT_EQ(0.095489956, Sph.getEigenvalue(1));
    EXPECT_FLOAT_EQ(0.031780295, Sph.getEigenvalue(2));
    EXPECT_FLOAT_EQ(1., Sph.getEigenvalue(0) + Sph.getEigenvalue(1) + Sph.getEigenvalue(2));
  }

TEST_F() [6/9]

TEST_F	(	eventShapeCoreAlgorithmTest	,
		Thrust
	)

Test the calculation of a thrust axis.

Definition at line 38 of file eventShapeCoreAlgorithm.cc.

  {
    std::vector<ROOT::Math::XYZVector> momenta;
    // random generated numbers
    momenta.emplace_back(0.5935352844151847, 0.28902324918117417, 0.9939000705771412);
    momenta.emplace_back(0.7097025137911714, 0.5118418422879152, 0.44501044145648994);
    momenta.emplace_back(0.6005771199332856, 0.12366608492454145, 0.7541373665256832);
    momenta.emplace_back(0.8548902083897467, 0.6887268865943484, 0.34301136659215437);
    momenta.emplace_back(0.26863521039211535, 0.011148638667487942, 0.96186334199901);
 
    const ROOT::Math::XYZVector thrustB = Thrust::calculateThrust(momenta);
 
    EXPECT_FLOAT_EQ(0.925838, thrustB.R());
    EXPECT_FLOAT_EQ(0.571661, thrustB.X());
    EXPECT_FLOAT_EQ(0.306741, thrustB.Y());
    EXPECT_FLOAT_EQ(0.660522, thrustB.Z());
  }

TEST_F() [7/9]

TEST_F	(	PIDPriorsTest	,
		PIDPriorsTableTest
	)

Test of the PIDPriorsTable class.

Definition at line 35 of file PIDPriors.cc.

  {
 
    std::vector<float> edgesX = {0., 0.2, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1., 1.1, 1.2, 1.3, 1.4, 1.5, 1.7, 1.9, 2.2, 2.5, 3., 3.5, 4., 5., 6., 7., 10.};
    std::vector<float> edgesY = { -1., -0.95, -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 0.95, 1.};
 
    PIDPriorsTable PIDTable = PIDPriorsTable();
    PIDTable.setBinEdges(edgesX, edgesY);
 
    PIDTable.printPrior();
 
    PIDTable.setPriorValue(1.1, 0.54, 0.1);
    PIDTable.setPriorValue(9.5, -0.67, 0.3);
 
    std::cout << "after setting (1.1, 0.54) amd (9.5, -0.67)" << std::endl;
    PIDTable.printPrior();
 
    EXPECT_B2WARNING(PIDTable.setPriorValue(11.0, -10.01, 3.));
    EXPECT_B2WARNING(PIDTable.setPriorValue(4.23, -0.01, 4.));
 
    std::cout << " ----- Setting the priors from outside ------ " << std::endl;
 
    std::vector<float>  prior;
 
    for (int j = 0; j < (int)edgesY.size() - 1; j++) {
      for (int i = 0; i < (int)edgesX.size() - 1; i++) {
        prior.push_back(0.001 * (i + j * (edgesX.size() - 1)));
      }
    }
 
 
    PIDTable.setPriorsTable(prior);
 
    PIDTable.printPrior();
 
    EXPECT_FLOAT_EQ(PIDTable.getPriorValue(0.1, -0.98), 0.001 * (0 + 0 * (edgesX.size() - 1)));
    EXPECT_FLOAT_EQ(PIDTable.getPriorValue(0.71, -0.65), 0.001 * (5 + 4 * (edgesX.size() - 1)));
    EXPECT_FLOAT_EQ(PIDTable.getPriorValue(0.71, -1.65), 0.);
    EXPECT_FLOAT_EQ(PIDTable.getPriorValue(11.71, 0.876), 0.);
    EXPECT_B2WARNING(PIDTable.getPriorValue(0.71, -1.65));
    EXPECT_B2WARNING(PIDTable.getPriorValue(11.71, 0.876));
    EXPECT_FLOAT_EQ(0., PIDTable.getPriorValue(11.71, 0.876));
 
 
 
    // now let's try to construct a table with no or on edge only
    std::vector<float> badEdges1 = {0.};
    std::vector<float> badEdges2 = {};
 
 
    PIDPriorsTable PIDTable2 = PIDPriorsTable();
 
    EXPECT_B2WARNING(PIDTable2.setBinEdges(edgesX, badEdges1));
    EXPECT_B2WARNING(PIDTable2.setBinEdges(edgesX, badEdges2));
 
 
    PIDTable2.setBinEdges(edgesX, badEdges2);
 
    PIDTable2.printPrior();
 
    PIDTable2.setPriorValue(1.1, 0.54, 0.3);
    PIDTable2.setPriorValue(9.5, -0.67, 0.3);
 
    std::cout << "after setting (1.1, 0.54) and (9.5, -0.67)" << std::endl;
    PIDTable2.printPrior();
 
    EXPECT_B2WARNING(PIDTable2.setPriorValue(11.0, -10.01, 0.3));
    EXPECT_FLOAT_EQ(PIDTable2.getPriorValue(9.516, 0.), 0.3);
 
  }

TEST_F() [8/9]

TEST_F	(	PIDPriorsTest	,
		PIDPriorTest
	)

Test the PIDPriors dbobject.

Definition at line 108 of file PIDPriors.cc.

  {
 
    // Creates 6 prior objects in different formats (array, PIDPriorsTable, TH2F) to test different setters
 
    std::vector<float> edgesX = {0., 1, 2, 3, 4, 5, 6, 7, 8, 9, 10.};
    std::vector<float> edgesY = { -1., -0.9, -0.8, -0.7, -0.6, -0.5, -0.4, -0.3, -0.2, -0.1, 0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.};
 
    std::vector<float>  prior1;
    std::vector<float>  prior2;
    std::vector<float>  prior3;
    std::vector<float>  prior4;
    std::vector<float>  prior5;
    std::vector<float>  prior6;
 
    TH2F* counts1 = new TH2F("counts1", "counts1",  10, 0, 10., 20, -1, 1);
    TH2F* counts2 = new TH2F("counts2", "counts2",  10, 0, 10., 20, -1, 1);
    TH2F* counts3 = new TH2F("counts3", "counts3",  10, 0, 10., 20, -1, 1);
    TH2F* counts4 = new TH2F("counts4", "counts4",  10, 0, 10., 20, -1, 1);
    TH2F* counts5 = new TH2F("counts5", "counts5",  10, 0, 10., 20, -1, 1);
    TH2F* counts6 = new TH2F("counts6", "counts6",  10, 0, 10., 20, -1, 1);
 
    TH2F* norms = new TH2F("norms", "norms",  10, 0, 10., 20, -1, 1);
 
    TRandom3* rn = new TRandom3();
 
 
 
    for (int i = 0; i < (int)edgesX.size() - 1; i++) {
      for (int j = 0; j < (int)edgesY.size() - 1; j++) {
 
        float n1 = (float)rn->Poisson(100);
        float n2 = (float)rn->Poisson(100);
        float n3 = (float)rn->Poisson(700);
        float n4 = (float)rn->Poisson(200);
        float n5 = (float)rn->Poisson(100);
        float n6 = (float)rn->Poisson(10);
 
        float n = n1 + n2 + n3 + n4 + n5 + n6;
 
        prior1.push_back(n1 / n);
        prior2.push_back(n2 / n);
        prior3.push_back(n3 / n);
        prior4.push_back(n4 / n);
        prior5.push_back(n5 / n);
        prior6.push_back(n6 / n);
 
        counts1->SetBinContent(i + 1, j + 1, n1);
        counts2->SetBinContent(i + 1, j + 1, n2);
        counts3->SetBinContent(i + 1, j + 1, n3);
        counts4->SetBinContent(i + 1, j + 1, n4);
        counts5->SetBinContent(i + 1, j + 1, n5);
        counts6->SetBinContent(i + 1, j + 1, n6);
        norms->SetBinContent(i + 1, j + 1, n);
 
      }
    }
 
 
    PIDPriorsTable PIDTable1 = PIDPriorsTable();
    PIDTable1.setBinEdges(edgesX, edgesY);
    PIDTable1.setPriorsTable(prior1);
 
    PIDPriorsTable PIDTable2 = PIDPriorsTable();
    PIDTable2.setBinEdges(edgesX, edgesY);
    PIDTable2.setPriorsTable(prior2);
 
    PIDPriorsTable PIDTable3 = PIDPriorsTable();
    PIDTable3.setBinEdges(edgesX, edgesY);
    PIDTable3.setPriorsTable(prior3);
 
    PIDPriorsTable PIDTable4 = PIDPriorsTable();
    PIDTable4.setBinEdges(edgesX, edgesY);
    PIDTable4.setPriorsTable(prior4);
 
    PIDPriorsTable PIDTable5 = PIDPriorsTable();
    PIDTable5.setBinEdges(edgesX, edgesY);
    PIDTable5.setPriorsTable(prior5);
 
    PIDPriorsTable PIDTable6 = PIDPriorsTable();
    PIDTable6.setBinEdges(edgesX, edgesY);
    PIDTable6.setPriorsTable(prior6);
 
    Const::ChargedStable part_e = Const::ChargedStable(11);
    Const::ChargedStable part_mu = Const::ChargedStable(13);
    Const::ChargedStable part_pi = Const::ChargedStable(211);
    Const::ChargedStable part_K = Const::ChargedStable(321);
    Const::ChargedStable part_p = Const::ChargedStable(2212);
    Const::ChargedStable part_d = Const::ChargedStable(1000010020);
 
 
    std::cout << "setting the tables " << std::endl;
    // setter with the priortable
    PIDPriors priors = PIDPriors();
    priors.setPriors(part_e, PIDTable1);
    priors.setPriors(part_mu, PIDTable2);
    priors.setPriors(part_pi, PIDTable3);
    priors.setPriors(part_K, PIDTable4);
    priors.setPriors(part_p, PIDTable5);
    priors.setPriors(part_d, PIDTable6);
 
    // setter with the count histograms
    PIDPriors priors2 = PIDPriors();
    priors2.setPriors(part_e, counts1, norms);
    priors2.setPriors(part_mu, counts2, norms);
    priors2.setPriors(part_pi, counts3, norms);
    priors2.setPriors(part_K, counts4, norms);
    priors2.setPriors(part_p, counts5, norms);
    priors2.setPriors(part_d, counts6, norms);
 
    counts1->Divide(norms);
    counts2->Divide(norms);
    counts3->Divide(norms);
    counts4->Divide(norms);
    counts5->Divide(norms);
    counts6->Divide(norms);
 
    std::cout << "division done" << std::endl;
    // setter with priors stored in a TH2F
    std::cout << "setting the TH2 " << std::endl;
    PIDPriors priors3 = PIDPriors();
    priors3.setPriors(part_e, counts1);
    priors3.setPriors(part_mu, counts2);
    priors3.setPriors(part_pi, counts3);
    priors3.setPriors(part_K, counts4);
    priors3.setPriors(part_p, counts5);
    priors3.setPriors(part_d, counts6);
 
 
    std::vector<Const::ChargedStable> codes = {part_e, part_mu, part_pi, part_K, part_p, part_d};
    for (auto code : codes) {
      for (int i = 0; i < 3; i++) {
        double p = rn->Uniform(0., 10);
        double c = rn->Uniform(-1., 1);
        EXPECT_FLOAT_EQ(priors.getPriorValue(code, p, c),  priors2.getPriorValue(code, p, c));
        EXPECT_FLOAT_EQ(priors.getPriorValue(code, p, c),  priors3.getPriorValue(code, p, c));
      }
    }
 
    EXPECT_FLOAT_EQ(0.,  priors.getPriorValue(part_pi, 11.1, 0.));
    EXPECT_FLOAT_EQ(0.,  priors.getPriorValue(part_pi, -1, 0.5));
    EXPECT_FLOAT_EQ(0.,  priors.getPriorValue(part_pi, 0.98, -1.001));
 
    EXPECT_FLOAT_EQ(0.,  priors2.getPriorValue(part_pi, 11.1, 0.));
    EXPECT_FLOAT_EQ(0.,  priors2.getPriorValue(part_pi, -1, 0.5));
    EXPECT_FLOAT_EQ(0.,  priors2.getPriorValue(part_pi, 0.98, -1.001));
 
    EXPECT_FLOAT_EQ(0.,  priors3.getPriorValue(part_pi, 11.1, 0.));
    EXPECT_FLOAT_EQ(0.,  priors3.getPriorValue(part_pi, -1, 0.5));
    EXPECT_FLOAT_EQ(0.,  priors3.getPriorValue(part_pi, 0.98, -1.001));
 
  }

TEST_F() [9/9]

TEST_F	(	TrackIsoScoreCalculatorTest	,
		TestDBRep
	)

Test correct retrieval of information from the database representation inner structure.

Definition at line 163 of file trackIsoScoreCalculator.cc.

  {
 
    PIDDetectorWeights dbrep("tree", m_dummyFile);
 
    // Test for correct filling of the RDataFrame.
    auto pdgIds = dbrep.getWeightsRDF().Take<double>("pdgId").GetValue();
    for (const auto& pdgId : pdgIds) {
      EXPECT_EQ(pdgId, m_testHypo.getPDGCode());
    }
 
    // Retrieve weight for a (p, theta) pair in the available weights range.
    auto p = 1.23; // GeV/c
    auto theta = 1.34; // rad
    auto weight = dbrep.getWeight(m_testHypo, m_detector, p, theta);
    EXPECT_EQ(weight, -0.118);
 
    // Test for weight in case of out-of-range p and/or theta values.
    p = 1.46;
    theta = 0.12;
    weight = dbrep.getWeight(m_testHypo, m_detector, p, theta);
    EXPECT_TRUE(std::isnan(weight));
 
    // Trigger a FATAL if reading an ill-defined source table.
    EXPECT_B2FATAL(PIDDetectorWeights("tree_broken", m_dummyFile));
 
  }