Belle II Software development
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, WeightInfoWeightMap
 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 >, DecayStringParticleDecayString
 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.

30 {
31 // Output
32 StoreArray<ContinuumSuppression> qqArray(maskName);
33 // Create ContinuumSuppression object
34 ContinuumSuppression* qqVars = qqArray.appendNew();
35
36 // Create relation: Particle <-> ContinuumSuppression
37 particle->addRelationTo(qqVars);
38
39 std::vector<ROOT::Math::XYZVector> p3_cms_sigB, p3_cms_roe, p3_cms_all;
40
41 std::vector<std::pair<ROOT::Math::XYZVector, int>> p3_cms_q_sigA;
42 std::vector<std::pair<ROOT::Math::XYZVector, int>> p3_cms_q_sigB;
43 std::vector<std::pair<ROOT::Math::XYZVector, int>> p3_cms_q_roe;
44
45 std::vector<float> ksfwFS0;
46 std::vector<float> ksfwFS1;
47
48 std::vector<float> cleoConesAll;
49 std::vector<float> cleoConesROE;
50
51 double et[2];
52
53 ROOT::Math::XYZVector thrustB;
54 ROOT::Math::XYZVector thrustO;
55
56 float thrustBm = -1;
57 float thrustOm = -1;
58 float cosTBTO = -1;
59 float cosTBz = -1;
60 float R2 = -1;
61
62
63 // -- B Cand --------------------------------------------------------------------------
64 PCmsLabTransform T;
65 double BeamEnergy = T.getCMSEnergy() / 2;
66
67 ROOT::Math::PxPyPzEVector p_cms_missA(0, 0, 0, 2 * BeamEnergy);
68 ROOT::Math::PxPyPzEVector p_cms_missB(0, 0, 0, 2 * BeamEnergy);
69 et[0] = et[1] = 0;
70
71 // -- SIG A --- Use B primary daughters - (Belle: use_finalstate_for_sig == 0) --------
72 std::vector<Belle2::Particle*> signalDaughters = particle->getDaughters();
73
74 for (const Belle2::Particle* sigFS0 : signalDaughters) {
75 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * sigFS0->get4Vector();
76
77 p3_cms_q_sigA.emplace_back(p_cms.Vect(), sigFS0->getCharge());
78
79 p_cms_missA -= p_cms;
80 et[0] += p_cms.Pt();
81 }
82
83 // -- SIG B --- Use B final-state daughters - (Belle: use_finalstate_for_sig == 1) ----
84 std::vector<const Belle2::Particle*> signalFSParticles = particle->getFinalStateDaughters();
85
86 for (const Belle2::Particle* sigFS1 : signalFSParticles) {
87 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * sigFS1->get4Vector();
88
89 p3_cms_all.push_back(p_cms.Vect());
90 p3_cms_sigB.push_back(p_cms.Vect());
91
92 p3_cms_q_sigB.emplace_back(p_cms.Vect(), sigFS1->getCharge());
93
94 p_cms_missB -= p_cms;
95 et[1] += p_cms.Pt();
96 }
97
98 // -- ROE -----------------------------------------------------------------------------
99 const RestOfEvent* roe = particle->getRelated<RestOfEvent>();
100
101 if (roe) {
102
103 // Charged tracks
104 //
105 std::vector<const Particle*> chargedROEParticles = roe->getChargedParticles(maskName);
106
107 for (const Particle* chargedROEParticle : chargedROEParticles) {
108
109 // TODO: Add helix and KVF with IpProfile once available. Port from L163-199 of:
110 // /belle/b20090127_0910/src/anal/ekpcontsuppress/src/ksfwmoments.cc
111
112 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * chargedROEParticle->get4Vector();
113
114 p3_cms_all.push_back(p_cms.Vect());
115 p3_cms_roe.push_back(p_cms.Vect());
116
117 p3_cms_q_roe.emplace_back(p_cms.Vect(), chargedROEParticle->getCharge());
118
119 p_cms_missA -= p_cms;
120 p_cms_missB -= p_cms;
121 et[0] += p_cms.Pt();
122 et[1] += p_cms.Pt();
123 }
124
125 // ECLCluster
126 //
127 std::vector<const Particle*> roePhotons = roe->getPhotons(maskName);
128
129 for (const Particle* photon : roePhotons) {
130
131 if (photon->getECLClusterEHypothesisBit() == ECLCluster::EHypothesisBit::c_nPhotons) {
132
133 ROOT::Math::PxPyPzEVector p_cms = T.rotateLabToCms() * photon->get4Vector();
134 p3_cms_all.push_back(p_cms.Vect());
135 p3_cms_roe.push_back(p_cms.Vect());
136
137 p3_cms_q_roe.emplace_back(p_cms.Vect(), photon->getCharge());
138
139 p_cms_missA -= p_cms;
140 p_cms_missB -= p_cms;
141 et[0] += p_cms.Pt();
142 et[1] += p_cms.Pt();
143 }
144 }
145
146 // Thrust variables
147 thrustB = Thrust::calculateThrust(p3_cms_sigB);
148 thrustO = Thrust::calculateThrust(p3_cms_roe);
149 thrustBm = thrustB.R();
150 thrustOm = thrustO.R();
151 cosTBTO = fabs(thrustB.Unit().Dot(thrustO.Unit()));
152 cosTBz = fabs(cos(thrustB.Theta()));
153
154 // Cleo Cones
155 CleoCones cc(p3_cms_all, p3_cms_roe, thrustB, true, true);
156 cleoConesAll = cc.cleo_cone_with_all();
157 cleoConesROE = cc.cleo_cone_with_roe();
158
159 // Fox-Wolfram Moments: Uses all final-state tracks (= sigB + ROE)
160 FoxWolfram FW(p3_cms_all);
161 FW.calculateBasicMoments();
162 R2 = FW.getR(2);
163
164 // KSFW moments
165 ROOT::Math::PxPyPzEVector p_cms_B = T.rotateLabToCms() * particle->get4Vector();
166 double Hso0_max(2 * (2 * BeamEnergy - p_cms_B.E()));
167 KsfwMoments KsfwM(Hso0_max,
168 p3_cms_q_sigA,
169 p3_cms_q_sigB,
170 p3_cms_q_roe,
171 p_cms_missA,
172 p_cms_missB,
173 et);
174 // use_finalstate_for_sig == 0
175 KsfwM.usefinal(0);
176 ksfwFS0.push_back(KsfwM.mm2());
177 ksfwFS0.push_back(KsfwM.et());
178 ksfwFS0.push_back(KsfwM.Hso(0, 0));
179 ksfwFS0.push_back(KsfwM.Hso(0, 1));
180 ksfwFS0.push_back(KsfwM.Hso(0, 2));
181 ksfwFS0.push_back(KsfwM.Hso(0, 3));
182 ksfwFS0.push_back(KsfwM.Hso(0, 4));
183 ksfwFS0.push_back(KsfwM.Hso(1, 0));
184 ksfwFS0.push_back(KsfwM.Hso(1, 2));
185 ksfwFS0.push_back(KsfwM.Hso(1, 4));
186 ksfwFS0.push_back(KsfwM.Hso(2, 0));
187 ksfwFS0.push_back(KsfwM.Hso(2, 2));
188 ksfwFS0.push_back(KsfwM.Hso(2, 4));
189 ksfwFS0.push_back(KsfwM.Hoo(0));
190 ksfwFS0.push_back(KsfwM.Hoo(1));
191 ksfwFS0.push_back(KsfwM.Hoo(2));
192 ksfwFS0.push_back(KsfwM.Hoo(3));
193 ksfwFS0.push_back(KsfwM.Hoo(4));
194 // use_finalstate_for_sig == 1
195 KsfwM.usefinal(1);
196 ksfwFS1.push_back(KsfwM.mm2());
197 ksfwFS1.push_back(KsfwM.et());
198 ksfwFS1.push_back(KsfwM.Hso(0, 0));
199 ksfwFS1.push_back(KsfwM.Hso(0, 1));
200 ksfwFS1.push_back(KsfwM.Hso(0, 2));
201 ksfwFS1.push_back(KsfwM.Hso(0, 3));
202 ksfwFS1.push_back(KsfwM.Hso(0, 4));
203 ksfwFS1.push_back(KsfwM.Hso(1, 0));
204 ksfwFS1.push_back(KsfwM.Hso(1, 2));
205 ksfwFS1.push_back(KsfwM.Hso(1, 4));
206 ksfwFS1.push_back(KsfwM.Hso(2, 0));
207 ksfwFS1.push_back(KsfwM.Hso(2, 2));
208 ksfwFS1.push_back(KsfwM.Hso(2, 4));
209 ksfwFS1.push_back(KsfwM.Hoo(0));
210 ksfwFS1.push_back(KsfwM.Hoo(1));
211 ksfwFS1.push_back(KsfwM.Hoo(2));
212 ksfwFS1.push_back(KsfwM.Hoo(3));
213 ksfwFS1.push_back(KsfwM.Hoo(4));
214
215 // TODO: The following is from the original belle ksfwmoments.cc module.
216 // Not sure if necessary here (i.e., will we be using rooksfw in belle II in the same way?).
217 // printf("rooksfw::rooksfw: mm2=%f et=%f hoo2=%f hso02=%f\n",
218 // m_mm2[0], et[0], m_Hoo[0][2], m_Hso[0][0][2]);
219 }
220
221 // Fill ContinuumSuppression with content
222 qqVars->addThrustB(thrustB);
223 qqVars->addThrustO(thrustO);
224 qqVars->addThrustBm(thrustBm);
225 qqVars->addThrustOm(thrustOm);
226 qqVars->addCosTBTO(cosTBTO);
227 qqVars->addCosTBz(cosTBz);
228 qqVars->addR2(R2);
229 qqVars->addKsfwFS0(ksfwFS0);
230 qqVars->addKsfwFS1(ksfwFS1);
231 qqVars->addCleoConesALL(cleoConesAll);
232 qqVars->addCleoConesROE(cleoConesROE);
233 }
Class to store reconstructed particles.
Definition: Particle.h:75

◆ 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.

24 {
27
28 // ----------------------------------------------------------------------
29 // Calculate momentum flow in 9 cones for all particles in event
30 // ----------------------------------------------------------------------
31 if (calc_CleoCones_with_all == true) {
32 for (int i = 1; i <= 9; i++) {
33 float momentum_flow_all = 0;
34 for (auto& iter0 : p3_cms_all) {
35
36 /* Use the following intervals
37 0*10<= <1*10 0- 10 170-180 180-1*10< <=180-0*10
38 1*10<= <2*10 10- 20 160-170 180-2*10< <=180-1*10
39 2*10<= <3*10 20- 30 150-160 180-3*10< <=180-2*10
40 3*10<= <4*10 30- 40 140-150 180-4*10< <=180-3*10
41 4*10<= <5*10 40- 50 130-140 180-5*10< <=180-4*10
42 5*10<= <6*10 50- 60 120-130 180-6*10< <=180-5*10
43 6*10<= <7*10 60- 70 110-120 180-7*10< <=180-6*10
44 7*10<= <8*10 70- 80 100-110 180-8*10< <=180-7*10
45 8*10<= <9*10 80- 90 90-100 180-9*10< <=180-8*10
46 ==90 */
47
48 float angle = ((180 * acos(thrustB.Unit().Dot(iter0.Unit()))) / M_PI);
49 if (((((i - 1) * 10) <= angle) && (angle < (i * 10))) || (((180 - (i * 10)) < angle) && (angle <= (180 - ((i - 1) * 10))))) {
50 momentum_flow_all += iter0.R();
51 // 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 );
52 }
53 if ((i == 9) && (angle == 90)) {
54 momentum_flow_all += iter0.R();
55 }
56 }
57 m_cleo_cone_with_all.push_back(momentum_flow_all);
58 }
59 }
60
61 // ----------------------------------------------------------------------
62 // Calculate momentum flow in 9 cones for all particles in rest of event
63 // ----------------------------------------------------------------------
64 if (calc_CleoCones_with_roe == true) {
65 for (int i = 1; i <= 9; i++) {
66 float momentum_flow_roe = 0;
67 for (auto& iter1 : p3_cms_roe) {
68 float angle = ((180 * acos(thrustB.Unit().Dot(iter1.Unit()))) / M_PI);
69 if (((((i - 1) * 10) <= angle) && (angle < (i * 10))) || (((180 - (i * 10)) < angle) && (angle <= (180 - ((i - 1) * 10))))) {
70 momentum_flow_roe += iter1.R();
71 }
72 if ((i == 9) && (angle == 90)) {
73 momentum_flow_roe += iter1.R();
74 }
75 }
76 m_cleo_cone_with_roe.push_back(momentum_flow_roe);
77 }
78 }
79 }
std::vector< float > m_cleo_cone_with_roe
Cleo Cones calculated from only ROE tracks.
Definition: CleoCones.h:50
std::vector< float > m_cleo_cone_with_all
Cleo Cones calculated from all tracks.
Definition: CleoCones.h:49

◆ create()

void create ( )

Create the list objects.

Warning
This function should be called in Every event

Definition at line 40 of file ParticleListHelper.cc.

41 {
42 m_list.create();
43 m_list->initialize(m_pdg, m_list.getName());
44 if (m_antiList) {
45 m_antiList->create();
46 (*m_antiList)->initialize(-m_pdg, m_antiList->getName());
47 m_list->bindAntiParticleList(**m_antiList);
48 }
49 }
StoreObjPtr< ParticleList > m_list
Store object for the list.
int m_pdg
PDG code of the list.
std::optional< StoreObjPtr< ParticleList > > m_antiList
Optional store object for the conjugated list if that exists.

◆ createCurrentParticle()

Particle createCurrentParticle ( ) const
private

Create current particle object.

Definition at line 439 of file ParticleCombiner.cc.

440 {
441 double px = 0;
442 double py = 0;
443 double pz = 0;
444 double E = 0;
445 for (const Particle* d : m_particles) {
446 px += d->getPx();
447 py += d->getPy();
448 pz += d->getPz();
449 E += d->getEnergy();
450 }
451 const ROOT::Math::PxPyPzEVector vec(px, py, pz, E);
452
453 switch (m_iParticleType) {
460 case 2: return Particle(vec, m_pdgCode, Particle::c_Unflavored, m_indices,
463 default: B2FATAL("You called getCurrentParticle although loadNext should have returned false!");
464 }
465
466 return Particle(); // This should never happen
467 }
R E
internal precision of FFTW codelets
unsigned int m_iParticleType
The type of particle which is currently generated.
std::vector< Particle * > m_particles
Pointers to the particle objects of the current combination.
const StoreArray< Particle > m_particleArray
Global list of particles.
std::vector< int > m_daughterProperties
Daughter's particle properties.
int m_pdgCode
PDG Code of the particle which is combined.
int m_properties
Particle property.
std::vector< int > m_indices
Indices stored in the ParticleLists of the current combination.
bool m_isSelfConjugated
True if the combined particle is self-conjugated.
@ c_Unflavored
Is its own antiparticle or we don't know whether it is a particle/antiparticle.
Definition: Particle.h:95
@ c_Flavored
Is either particle or antiparticle.
Definition: Particle.h:96
TClonesArray * getPtr() const
Raw access to the underlying TClonesArray.
Definition: StoreArray.h:311

◆ 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.

475 {
476 std::vector<Particle*> stack = m_particles;
477 static std::vector<int> sources; // stack for particle sources
478 sources.clear();
479
480 // recursively check all daughters and daughters of daughters
481 while (!stack.empty()) {
482 Particle* p = stack.back();
483 stack.pop_back();
484 const std::vector<int>& daughters = p->getDaughterIndices();
485
486 if (daughters.empty()) {
487 int source = p->getMdstSource();
488 for (int i : sources) {
489 if (source == i) return false;
490 }
491 sources.push_back(source);
492 } else {
493 for (int daughter : daughters) stack.push_back(m_particleArray[daughter]);
494 }
495 }
496 return true;
497 }

◆ 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.

523 {
524 unsigned nECLSource = 0;
525 std::vector<Particle*> stack = m_particles;
526 static std::vector<int> connectedregions;
527 static std::vector<ECLCluster::EHypothesisBit> hypotheses;
528 connectedregions.clear();
529 hypotheses.clear();
530
531 // recursively check all daughters and daughters of daughters
532 while (!stack.empty()) {
533 Particle* p = stack.back();
534 stack.pop_back();
535 const std::vector<int>& daughters = p->getDaughterIndices();
536
537 if (daughters.empty()) {
538 // Only test if the particle was created from an ECLCluster at source.
539 // This CAN CHANGE if we change the cluster <--> track matching,
540 // (currently we match nPhotons clusters to all track type particles,
541 // i.e. electrons, pions, kaons etc. We might gain by matching
542 // neutralHadron hypothesis clusters to kaons, for example).
543 //
544 // In the above case one would need to extend the check to ALL
545 // particles with an associated ECLCluster. Then replace the following
546 // active line with two lines...
547 //
548 // auto cluster = p->getECLCluster();
549 // if (cluster) { // then do stuff
550 if (p->getParticleSource() == Particle::EParticleSourceObject::c_ECLCluster) {
551 nECLSource++;
552 auto* cluster = p->getECLCluster();
553 connectedregions.push_back(cluster->getConnectedRegionId());
554 hypotheses.push_back(p->getECLClusterEHypothesisBit());
555 }
556 } else {
557 for (int daughter : daughters) stack.push_back(m_particleArray[daughter]);
558 }
559 }
560
561 // less than two particles from an ECL source is fine
562 // (unless cluster <--> track matching changes)
563 if (nECLSource < 2) return true;
564
565 // yes this is a nested for loop but it's fast, we promise
566 for (unsigned icr = 0; icr < connectedregions.size(); ++icr)
567 for (unsigned jcr = icr + 1; jcr < connectedregions.size(); ++jcr)
568 if (connectedregions[icr] == connectedregions[jcr])
569 if (hypotheses[icr] != hypotheses[jcr]) return false;
570
571 return true;
572 }

◆ 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.

501 {
502 std::set<int> indexSet;
503 if (not m_inputListsCollide)
504 indexSet.insert(m_indices.begin(), m_indices.end());
505 else
506 for (unsigned int i = 0; i < m_numberOfLists; i++)
507 indexSet.insert(m_indicesToUniqueIDs.at(m_indices[i]));
508
509 bool elementInserted = m_usedCombinations.insert(indexSet).second;
510 return elementInserted;
511 }
std::unordered_set< std::set< int > > m_usedCombinations
already used combinations (as sets of indices or unique IDs).
unsigned int m_numberOfLists
Number of lists which are combined.
std::unordered_map< int, int > m_indicesToUniqueIDs
map of store array indices of input Particles to their unique IDs.
bool m_inputListsCollide
True if the daughter lists can contain copies of Particles.

◆ 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.

288 {
289 const Particle* A, *B;
290 bool copies = false;
291 for (int i : listA) {
292 bool aIsAlreadyIn = m_indicesToUniqueIDs.count(i) ? true : false;
293
294 if (not aIsAlreadyIn)
295 m_indicesToUniqueIDs[ i ] = uniqueID++;
296
297 for (int j : listB) {
298 bool bIsAlreadyIn = m_indicesToUniqueIDs.count(j) ? true : false;
299
300 if (bIsAlreadyIn)
301 continue;
302
303 // are these two particles copies
304 A = m_particleArray[ i ];
305 B = m_particleArray[ j ];
306 copies = B->isCopyOf(A);
307
308 if (copies)
310 }
311 }
312 }

◆ 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.

316 {
317 for (int i : listA) {
318 bool aIsAlreadyIn = m_indicesToUniqueIDs.count(i) ? true : false;
319
320 if (not aIsAlreadyIn)
321 m_indicesToUniqueIDs[ i ] = uniqueID++;
322 }
323 }

◆ find_decay()

bool find_decay ( const DecayNode to_find) const

Check if the decay node contains the given decay tree.

Parameters
to_findDecayNode object describing the decay

Definition at line 19 of file DecayNode.cc.

20 {
21 if (to_find == (*this))
22 return true;
23 for (const auto& node : daughters) {
24 if (node.find_decay(to_find))
25 return true;
26 }
27 return false;
28 }
std::vector< DecayNode > daughters
daughter decay nodes
Definition: DecayNode.h:48

◆ 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.

76 {
77 return indices;
78 }
std::vector< unsigned int > indices
The indices of the current loaded combination.

◆ 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.

106 {
107 return m_types;
108 }
std::vector< ParticleList::EParticleType > m_types
The current types of sublist of the ParticleLists for this combination.

◆ getCurrentParticle()

Particle getCurrentParticle ( ) const

Returns the particle.

Definition at line 469 of file ParticleCombiner.cc.

470 {
471 return m_current_particle;
472 }
Particle m_current_particle
The current Particle object generated by this combiner.

◆ 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.

575 {
576 if (not m_inputListsCollide)
577 return 0;
578
579 if (not m_indicesToUniqueIDs.count(index))
580 return -1;
581
582 return m_indicesToUniqueIDs.at(index);
583 }

◆ init() [1/3]

void init ( )

Initialises the generator to produce the given type of sublist.

Definition at line 221 of file ParticleCombiner.cc.

222 {
223 m_iParticleType = 0;
224
225 m_particleIndexGenerator.init(std::vector<unsigned int> {}); // ParticleIndexGenerator will be initialised on first call
226 m_listIndexGenerator.init(m_numberOfLists); // ListIndexGenerator must be initialised here!
227 m_usedCombinations.clear();
230
233 }
ListIndexGenerator m_listIndexGenerator
listIndexGenerator makes the combinations of the types of sublists of the ParticleLists
ParticleIndexGenerator m_particleIndexGenerator
particleIndexGenerator makes the combinations of indices stored in the sublists of the ParticleLists
void init(const std::vector< unsigned int > &_sizes)
Initialises the generator to produce combinations with the given sizes of each particle list.
void initIndicesToUniqueIDMap()
In the case input daughter particle lists collide (two or more lists contain copies of Particles) the...
void init(unsigned int _numberOfLists)
Initialises the generator to produce the given type of sublist.

◆ 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
_sizesthe sizes of the particle lists to combine

Definition at line 28 of file ParticleCombiner.cc.

29 {
30
31 m_numberOfLists = _sizes.size();
32 sizes = _sizes;
33
35
37 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
38 indices[i] = 0;
39 }
40
42 for (unsigned int i = 0; i < m_numberOfLists; ++i)
44
45 if (m_numberOfLists == 0) {
47 }
48
49 }
unsigned int m_iCombination
The current position of the combination.
unsigned int m_numberOfLists
Number of lists which are combined.
std::vector< unsigned int > sizes
The sizes of the particle lists which are combined.
unsigned int m_nCombinations
The total amount of combinations.

◆ init() [3/3]

void init ( unsigned int  _numberOfLists)

Initialises the generator to produce the given type of sublist.

Parameters
_numberOfListsNumber of Particle Lists which shall be combined

Definition at line 80 of file ParticleCombiner.cc.

81 {
82 m_numberOfLists = _numberOfLists;
84
87
88 }
unsigned int m_iCombination
The current position of the combination.
unsigned int m_numberOfLists
Number of lists which are combined.
unsigned int m_nCombinations
The total amount of combinations.

◆ 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.

236 {
237 m_indicesToUniqueIDs.clear();
238
239 unsigned inputParticlesCount = 0;
240 for (unsigned int i = 0; i < m_numberOfLists; ++i)
241 inputParticlesCount += m_plists[i]->getListSize();
242
243 m_indicesToUniqueIDs.reserve(inputParticlesCount);
244
245 int uniqueID = 1;
246
247 for (auto& collidingList : m_collidingLists) {
248 StoreObjPtr<ParticleList> listA = m_plists[collidingList.first];
249 StoreObjPtr<ParticleList> listB = m_plists[collidingList.second];
250
251 bool sameSign = (listA->getPDGCode() == listB->getPDGCode());
252
253 // if sameSign == true then
254 // 1. compare FS to FS particles in lists A and B
255 // 2. compare anti-FS to anti-FS particles in lists A and B
256 // else
257 // 1. compare FS to anti-FS particles in lists A and B
258 // 2. compare anti-FS to FS particles in lists A and B
259 // and in either case compare
260 // 3. compare SC to SC particles in lists A and B
261 if (sameSign) {
262 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false),
263 listB->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
264 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true),
265 listB->getList(ParticleList::c_FlavorSpecificParticle, true), uniqueID);
266 } else {
267 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false),
268 listB->getList(ParticleList::c_FlavorSpecificParticle, true), uniqueID);
269 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true),
270 listB->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
271 }
272
273 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_SelfConjugatedParticle),
274 listB->getList(ParticleList::c_SelfConjugatedParticle), uniqueID);
275 }
276
277 // assign unique ids to others as well
278 for (unsigned i = 0; i < m_numberOfLists; i++) {
279 StoreObjPtr<ParticleList> listA = m_plists[i];
280
281 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, true), uniqueID);
282 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_FlavorSpecificParticle, false), uniqueID);
283 fillIndicesToUniqueIDMap(listA->getList(ParticleList::c_SelfConjugatedParticle), uniqueID);
284 }
285 }
std::vector< StoreObjPtr< ParticleList > > m_plists
particle lists
std::vector< std::pair< unsigned, unsigned > > m_collidingLists
pairs of lists that can contain copies.
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.
bool sameSign(double expected, double actual)
Predicate checking that two values have the same sign.
Definition: TestHelpers.cc:76

◆ 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.

514 {
515 for (const auto& collidingList : m_collidingLists)
516 if (pair == collidingList)
517 return true;
518
519 return false;
520 }

◆ 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.

56 {
57 // Private member needs to be initialized. Here it is initialized to -1 (illegal value).
58 m_uf = -1;
59
60 m_et[0] = et[0];
61 m_et[1] = et[1];
62
63 m_mm2[0] = p_cms_missA.E() > 0
64 ? p_cms_missA.mag2()
65 : -p_cms_missA.E() * p_cms_missA.E() - p_cms_missA.P2();
66 m_mm2[1] = p_cms_missB.E() > 0
67 ? p_cms_missB.mag2()
68 : -p_cms_missB.E() * p_cms_missB.E() - p_cms_missB.P2();
69
70 //========================
71 // Calculate discriminants
72 //========================
73 std::vector<std::pair<ROOT::Math::XYZVector, int>>::iterator pqi, pqj;
74
75 // Calculate Hso components
76 for (int i = 0; i < 3; i++) {
77 for (int k = 0; k < 5; k++) {
78 m_Hso[0][i][k] = m_Hso[1][i][k] = 0;
79 }
80 }
81
82 // Signal A (use_finalstate_for_sig == 0)
83 for (pqi = p3_cms_q_sigA.begin(); pqi != p3_cms_q_sigA.end(); ++pqi) {
84 const double pi_mag((pqi->first).R());
85 for (pqj = p3_cms_q_roe.begin(); pqj != p3_cms_q_roe.end(); ++pqj) {
86 const double pj_mag((pqj->first).R());
87 const double ij_cos((pqi->first).Dot(pqj->first) / pi_mag / pj_mag);
88 const int c_or_n(0 == (pqj->second) ? 1 : 0); // 0: charged 1: neutral
89 for (int k = 0; k < 5; k++) {
90 m_Hso[0][c_or_n][k] += (k % 2)
91 ? (pqi->second) * (pqj->second) * pj_mag * legendre(ij_cos, k)
92 : pj_mag * legendre(ij_cos, k);
93 }
94 }
95 const double p_miss_mag(p_cms_missA.P());
96 const double i_miss_cos((pqi->first).Dot(p_cms_missA.Vect()) / pi_mag / p_miss_mag);
97 for (int k = 0; k < 5; k++) {
98 m_Hso[0][2][k] += (k % 2) ? 0 : p_miss_mag * legendre(i_miss_cos, k);
99 }
100 }
101
102 // Signal B (use_finalstate_for_sig == 1)
103 for (pqi = p3_cms_q_sigB.begin(); pqi != p3_cms_q_sigB.end(); ++pqi) {
104 const double pi_mag((pqi->first).R());
105 for (pqj = p3_cms_q_roe.begin(); pqj != p3_cms_q_roe.end(); ++pqj) {
106 const double pj_mag((pqj->first).R());
107 const double ij_cos((pqi->first).Dot(pqj->first) / pi_mag / pj_mag);
108 const int c_or_n(0 == (pqj->second) ? 1 : 0); // 0: charged 1: neutral
109 for (int k = 0; k < 5; k++) {
110 m_Hso[1][c_or_n][k] += (k % 2)
111 ? (pqi->second) * (pqj->second) * pj_mag * legendre(ij_cos, k)
112 : pj_mag * legendre(ij_cos, k);
113 }
114 }
115 const double p_miss_mag(p_cms_missB.P());
116 const double i_miss_cos((pqi->first).Dot(p_cms_missB.Vect()) / pi_mag / p_miss_mag);
117 for (int k = 0; k < 5; k++) {
118 m_Hso[1][2][k] += (k % 2) ? 0 : p_miss_mag * legendre(i_miss_cos, k);
119 }
120 }
121
122 // Add missing to the lists
123 std::vector<std::pair<ROOT::Math::XYZVector, int>> p3_cms_q_roeA(p3_cms_q_roe), p3_cms_q_roeB(p3_cms_q_roe);
124 p3_cms_q_roeA.emplace_back(p_cms_missA.Vect(), 0);
125 p3_cms_q_roeB.emplace_back(p_cms_missB.Vect(), 0);
126
127 // Calculate Hoo components
128 for (int k = 0; k < 5; k++) {
129 m_Hoo[0][k] = m_Hoo[1][k] = 0;
130 }
131 for (pqi = p3_cms_q_roeA.begin(); pqi != p3_cms_q_roeA.end(); ++pqi) {
132 const double pi_mag((pqi->first).R());
133 for (pqj = p3_cms_q_roeA.begin(); pqj != pqi; ++pqj) {
134 const double pj_mag((pqj->first).R());
135 const double ij_cos((pqi->first).Dot(pqj->first) / pi_mag / pj_mag);
136 for (int k = 0; k < 5; k++) {
137 m_Hoo[0][k] += (k % 2)
138 ? (pqi->second) * (pqj->second) * pi_mag * pj_mag * legendre(ij_cos, k)
139 : pi_mag * pj_mag * legendre(ij_cos, k);
140 }
141 }
142 }
143 for (pqi = p3_cms_q_roeB.begin(); pqi != p3_cms_q_roeB.end(); ++pqi) {
144 const double pi_mag((pqi->first).R());
145 for (pqj = p3_cms_q_roeB.begin(); pqj != pqi; ++pqj) {
146 const double pj_mag((pqj->first).R());
147 const double ij_cos((pqi->first).Dot(pqj->first) / pi_mag / pj_mag);
148 for (int k = 0; k < 5; k++) {
149 m_Hoo[1][k] += (k % 2)
150 ? (pqi->second) * (pqj->second) * pi_mag * pj_mag * legendre(ij_cos, k)
151 : pi_mag * pj_mag * legendre(ij_cos, k);
152 }
153 }
154 }
155
156 // Normalize so that it does not dependent on delta_e
157 for (int k = 0; k < 5; k++) {
158 for (int j = 0; j < ((k % 2) ? 1 : 3); j++) {
159 m_Hso[0][j][k] /= Hso0_max;
160 m_Hso[1][j][k] /= Hso0_max;
161 }
162 m_Hoo[0][k] /= (Hso0_max * Hso0_max);
163 m_Hoo[1][k] /= (Hso0_max * Hso0_max);
164 }
165
166 }
double et(int uf=-1) const
Returns calculated transverse energy.
Definition: KsfwMoments.h:93
double m_et[2]
Transverse energy.
Definition: KsfwMoments.h:127
int m_uf
Flag that specifiies we are using the final state for signal.
Definition: KsfwMoments.h:124
double m_Hso[2][3][5]
KSFW moments.
Definition: KsfwMoments.h:125
double m_mm2[2]
Missing mass squared.
Definition: KsfwMoments.h:128
double m_Hoo[2][5]
KSFW moments.
Definition: KsfwMoments.h:126
double legendre(const double z, const int i)
Legendre polynomials.
Definition: KsfwMoments.cc:33

◆ legendre()

double legendre ( const double  z,
const int  i 
)
inline

Legendre polynomials.

Definition at line 33 of file KsfwMoments.cc.

34 {
35 switch (i) {
36 case 0: return 1.;
37 case 1: return z;
38 case 2: return 1.5 * z * z - 0.5;
39 case 3: return z * (2.5 * z * z - 1.5);
40 case 4: return (4.375 * z * z * z * z - 3.75 * z * z + 0.375);
41 default: return 0;
42 }
43 }

◆ 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.

52 {
53
55 return false;
56 }
57
58 if (m_iCombination > 0) {
59
60 // TF SC this does not yet account for double counting so will produce:
61 // { 000, 100, 200, ... 010, 110, .... } even if the first and second
62 // place are the same particle list
63 for (unsigned int i = 0; i < m_numberOfLists; i++) {
64 indices[i]++;
65 if (indices[i] < sizes[i]) break;
66 indices[i] = 0;
67 }
68
69 }
70
72 return true;
73 }

◆ 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.

91 {
92
94 return false;
95
96 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
97 bool useSelfConjugatedDaughterList = (m_iCombination & (1 << i));
98 m_types[i] = useSelfConjugatedDaughterList ? ParticleList::c_SelfConjugatedParticle : ParticleList::c_FlavorSpecificParticle;
99 }
100
102 return true;
103 }

◆ 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.

327 {
328
329 bool loadedNext = false;
336 while (true) {
337 if (m_iParticleType == 0) {
338 loadedNext = loadNextParticle(false);
339 } else if (m_iParticleType == 1 and loadAntiParticle) {
340 loadedNext = loadNextParticle(true);
341 } else if (m_iParticleType == 2) {
342 loadedNext = loadNextSelfConjugatedParticle();
343 } else {
344 return false;
345 }
346
347 if (loadedNext) return true;
348 else ++m_iParticleType;
349
350 if (m_iParticleType == 2) {
351 std::vector<unsigned int> sizes(m_numberOfLists);
352 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
353 sizes[i] = m_plists[i]->getList(ParticleList::c_SelfConjugatedParticle, false).size();
354 }
356 } else {
358 }
359
360 }
361 }
bool loadNextSelfConjugatedParticle()
Loads the next combination.
bool loadNextParticle(bool useAntiParticle)
Loads the next combination.

◆ 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.

364 {
365 while (true) {
366
367 // Load next index combination if available
369
370 const auto& m_types = m_listIndexGenerator.getCurrentIndices();
371 const auto& indices = m_particleIndexGenerator.getCurrentIndices();
372
373 for (unsigned int i = 0; i < m_numberOfLists; i++) {
374 const auto& list = m_plists[i]->getList(m_types[i], m_types[i] == ParticleList::c_FlavorSpecificParticle ? useAntiParticle : false);
375 m_indices[i] = list[ indices[i] ];
377 }
378
379 if (not currentCombinationHasDifferentSources()) continue;
380
382 if (!m_cut->check(&m_current_particle))
383 continue;
384
385 if (not currentCombinationIsUnique()) continue;
386
387 if (not currentCombinationIsECLCRUnique()) continue;
388
389 return true;
390 }
391
392 // Load next list combination if available and reset indexCombiner
394 const auto& m_types = m_listIndexGenerator.getCurrentIndices();
395 std::vector<unsigned int> sizes(m_numberOfLists);
396 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
397 sizes[i] = m_plists[i]->getList(m_types[i], m_types[i] == ParticleList::c_FlavorSpecificParticle ? useAntiParticle : false).size();
398 }
400 continue;
401 }
402 return false;
403 }
404
405 }
std::unique_ptr< Variable::Cut > m_cut
cut object which performs the cuts
bool currentCombinationHasDifferentSources()
Check that all FS particles of a combination differ.
bool currentCombinationIsECLCRUnique()
Check that: if the current combination has at least two particles from an ECL source,...
const std::vector< unsigned int > & getCurrentIndices() const
Returns theindices of the current loaded combination.
bool currentCombinationIsUnique()
Check that the combination is unique.
const std::vector< ParticleList::EParticleType > & getCurrentIndices() const
Returns the type of the sublist of the current loaded combination.
Particle createCurrentParticle() const
Create current particle object.
bool loadNext()
Loads the next combination.

◆ loadNextSelfConjugatedParticle()

bool loadNextSelfConjugatedParticle ( )
private

Loads the next combination.

Returns false if there is no next combination

Definition at line 407 of file ParticleCombiner.cc.

408 {
409 while (true) {
410
411 // Load next index combination if available
413
414 const auto& indices = m_particleIndexGenerator.getCurrentIndices();
415
416 for (unsigned int i = 0; i < m_numberOfLists; i++) {
417 m_indices[i] = m_plists[i]->getList(ParticleList::c_SelfConjugatedParticle, false) [ indices[i] ];
419 }
420
421 if (not currentCombinationHasDifferentSources()) continue;
422
424 if (!m_cut->check(&m_current_particle))
425 continue;
426
427 if (not currentCombinationIsUnique()) continue;
428
429 if (not currentCombinationIsECLCRUnique()) continue;
430
431 return true;
432 }
433
434 return false;
435 }
436
437 }

◆ operator!=()

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

Not equal: See operator==.

Definition at line 65 of file DecayNode.cc.

66 {
67 return not(node1 == node2);
68 }

◆ 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
node1first node
node2second node

Definition at line 48 of file DecayNode.cc.

49 {
50 if (node1.pdg == node2.pdg) {
51 if (node1.daughters.size() > 0 and node2.daughters.size() > 0) {
52 if (node1.daughters.size() != node2.daughters.size())
53 return false;
54 for (unsigned int i = 0; i < node1.daughters.size(); ++i) {
55 if (node1.daughters[i] != node2.daughters[i])
56 if (node1.daughters[i] != 0 and node2.daughters[i] != 0)
57 return false;
58 }
59 }
60 return true;
61 }
62 return false;
63 }

◆ 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.

167 : m_iParticleType(0),
170 {
171 bool valid = decaydescriptor.isInitOK();
172 if (!valid)
173 B2ERROR("Given decaydescriptor failed to initialized");
174
175 m_properties = decaydescriptor.getProperty();
176
177 // Mother particle
178 const DecayDescriptorParticle* mother = decaydescriptor.getMother();
179 m_pdgCode = mother->getPDGCode();
180 m_properties |= mother->getProperty();
181
182 // Daughters
183 m_numberOfLists = decaydescriptor.getNDaughters();
184 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
185 // Get the mother of the subdecaystring of the ith daughter
186 // eg. "B -> [D -> K pi] [tau -> pi pi pi]". The 0th daughter is the
187 // *decaystring* D -> K pi whose mother is the D.
188 const DecayDescriptorParticle* daughter = decaydescriptor.getDaughter(i)->getMother();
189 StoreObjPtr<ParticleList> list(daughter->getFullName());
190 m_plists.push_back(list);
191
192 int daughterProperty = daughter->getProperty();
193 m_daughterProperties.push_back(daughterProperty);
194 }
195
196 m_cut = Variable::Cut::compile(cutParameter);
197
198 m_isSelfConjugated = decaydescriptor.isSelfConjugated();
199
200 // check if input lists can contain copies
201 // remember (list_i, list_j) pairs, where list_i and list_j can contain copies
202 m_inputListsCollide = false;
203 m_collidingLists.clear();
204 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
205 const DecayDescriptorParticle* daughter_i = decaydescriptor.getDaughter(i)->getMother();
206 for (unsigned int j = i + 1; j < m_numberOfLists; ++j) {
207 const DecayDescriptorParticle* daughter_j = decaydescriptor.getDaughter(j)->getMother();
208
209 if (abs(daughter_i->getPDGCode()) != abs(daughter_j->getPDGCode()))
210 continue;
211
212 if (daughter_i->getLabel() != daughter_j->getLabel()) {
213 m_inputListsCollide = true;
214 m_collidingLists.emplace_back(i, j);
215 }
216 }
217 }
218
219 }
static std::unique_ptr< GeneralCut > compile(const std::string &cut)
Creates an instance of a cut and returns a unique_ptr to it, if you need a copy-able object instead y...
Definition: GeneralCut.h:84

◆ 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.

110 : m_iParticleType(0),
113 {
114
115 DecayDescriptor decaydescriptor;
116 bool valid = decaydescriptor.init(decayString);
117 if (!valid)
118 B2ERROR("Invalid input DecayString: " << decayString);
119
120 m_properties = decaydescriptor.getProperty();
121
122 // Mother particle
123 const DecayDescriptorParticle* mother = decaydescriptor.getMother();
124 m_pdgCode = mother->getPDGCode();
125 m_properties |= mother->getProperty();
126
127 // Daughters
128 m_numberOfLists = decaydescriptor.getNDaughters();
129 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
130 // Get the mother of the subdecaystring of the ith daughter
131 // eg. "B -> [D -> K pi] [tau -> pi pi pi]". The 0th daughter is the
132 // *decaystring* D -> K pi whose mother is the D.
133 const DecayDescriptorParticle* daughter = decaydescriptor.getDaughter(i)->getMother();
134 StoreObjPtr<ParticleList> list(daughter->getFullName());
135 m_plists.push_back(list);
136
137 int daughterProperty = daughter->getProperty();
138 m_daughterProperties.push_back(daughterProperty);
139 }
140
141 m_cut = Variable::Cut::compile(cutParameter);
142
143 m_isSelfConjugated = decaydescriptor.isSelfConjugated();
144
145 // check if input lists can contain copies
146 // remember (list_i, list_j) pairs, where list_i and list_j can contain copies
147 m_inputListsCollide = false;
148 m_collidingLists.clear();
149 for (unsigned int i = 0; i < m_numberOfLists; ++i) {
150 const DecayDescriptorParticle* daughter_i = decaydescriptor.getDaughter(i)->getMother();
151 for (unsigned int j = i + 1; j < m_numberOfLists; ++j) {
152 const DecayDescriptorParticle* daughter_j = decaydescriptor.getDaughter(j)->getMother();
153
154 if (abs(daughter_i->getPDGCode()) != abs(daughter_j->getPDGCode()))
155 continue;
156
157 if (daughter_i->getLabel() != daughter_j->getLabel()) {
158 m_inputListsCollide = true;
159 m_collidingLists.emplace_back(i, j);
160 }
161 }
162 }
163
164 }

◆ 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.

31 {
32 std::stringstream output;
33
34 for (unsigned int i = 0; i < indent; ++i) {
35 output << " ";
36 }
37
38 output << std::to_string(pdg);
39 output << std::endl;
40
41 for (const auto& d : daughters) {
42 output << d.print_node(indent + 1);
43 }
44
45 return output.str();
46 }
int pdg
pdg code of the particle
Definition: DecayNode.h:47

◆ 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
decaydecay descriptor to use for the list name
saveif true the list will be saved in the output file

Definition at line 27 of file ParticleListHelper.cc.

28 {
29 auto listname = decay.getMother()->getFullName();
30 m_pdg = decay.getMother()->getPDGCode();
31 auto storeFlags = save ? DataStore::c_WriteOut : DataStore::c_DontWriteOut;
33 m_list.registerInDataStore(listname, storeFlags);
34 auto antiListName = ParticleListName::antiParticleListName(listname);
35 if (antiListName != listname) {
36 m_antiList.emplace(antiListName);
37 m_antiList->registerInDataStore(storeFlags);
38 }
39 }
@ c_WriteOut
Object/array should be saved by output modules.
Definition: DataStore.h:70
@ c_DontWriteOut
Object/array should be NOT saved by output modules.
Definition: DataStore.h:71
StoreArray< Particle > m_particles
Store array for the particles.
bool registerInDataStore(DataStore::EStoreFlags storeFlags=DataStore::c_WriteOut)
Register the object/array in the DataStore.
std::string antiParticleListName(const std::string &listName)
Returns name of anti-particle-list corresponding to listName.

◆ registerList() [2/2]

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

Register a list by name.

Warning
should be called in initialize
Parameters
listnamefull list name in the form '{particlename}:{label}'
saveif true the list will be saved in the output file

Definition at line 18 of file ParticleListHelper.cc.

19 {
20 DecayDescriptor decayDescriptor;
21 bool valid = decayDescriptor.init(listname);
22 if (!valid)
23 throw std::runtime_error("Invalid ParticleList name: '" + listname + "' Should be EVTPDLNAME[:LABEL], e.g. B+ or B+:mylist.");
24 ParticleListHelper::registerList(decayDescriptor, save);
25 }
void registerList(const std::string &listname, bool save=true)
Register a list by name.

◆ 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.

130 {
131
132 // Pick a random index in the clusterTheta, p, charge bins arrays.
133 // Exclude underflow and overflow.
134 std::random_device rd; // non-deterministic uniform int rand generator.
135 std::uniform_int_distribution<int> binx_idx_distr(1, m_thetabins.size() - 1);
136 std::uniform_int_distribution<int> biny_idx_distr(1, m_pbins.size() - 1);
137 std::uniform_int_distribution<int> binz_idx_distr(1, m_chbins.size() - 1);
138 int binx = binx_idx_distr(rd);
139 int biny = biny_idx_distr(rd);
140 int binz = binz_idx_distr(rd);
141
142 // Pick each axis' bin centre as a test value for (clusterTheta, p, charge).
143 auto theta = m_dbrep.getWeightCategories()->GetXaxis()->GetBinCenter(binx);
144 auto p = m_dbrep.getWeightCategories()->GetYaxis()->GetBinCenter(biny);
145 auto charge = m_dbrep.getWeightCategories()->GetZaxis()->GetBinCenter(binz);
146
147 int jth, ip, kch;
148 auto jik = m_dbrep.getMVAWeightIdx(theta, p, charge, jth, ip, kch);
149
150 EXPECT_EQ(jth, binx);
151 EXPECT_EQ(ip, biny);
152 EXPECT_EQ(kch, binz);
153
154 auto thisfname = m_basename
155 + "__clusterTheta__" + std::to_string(m_thetabins.at(jth - 1)) + "_" + std::to_string(m_thetabins.at(jth))
156 + "__p__" + std::to_string(m_pbins.at(ip - 1)) + "_" + std::to_string(m_pbins.at(ip))
157 + "__charge__" + std::to_string(charge);
158 std::replace(thisfname.begin(), thisfname.end(), '.', '_');
159 thisfname += ".xml";
160
161 EXPECT_EQ(thisfname, m_dummyfiles.at(jik));
162
163 auto matchitr = std::find(m_dummyfiles.begin(), m_dummyfiles.end(), thisfname);
164 auto thisidx = std::distance(m_dummyfiles.begin(), matchitr);
165
166 EXPECT_EQ(thisidx, jik);
167
168 }
double charge(int pdgCode)
Returns electric charge of a particle with given pdg code.
Definition: EvtPDLUtil.cc:44

◆ TEST_F() [2/9]

Test the calculation of the CleoClones variables.

Definition at line 57 of file eventShapeCoreAlgorithm.cc.

58 {
59 const bool use_all = true;
60 const bool use_roe = true;
61 std::vector<ROOT::Math::XYZVector> momenta;
62 std::vector<ROOT::Math::XYZVector> sig_side_momenta;
63 std::vector<ROOT::Math::XYZVector> roe_side_momenta;
64
65 // "Signal Side"
66 sig_side_momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
67 sig_side_momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
68 sig_side_momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
69 sig_side_momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
70 sig_side_momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
71
72 // "ROE Side"
73 roe_side_momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
74 roe_side_momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
75 roe_side_momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
76 roe_side_momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
77 roe_side_momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
78
79 // "All momenta"
80 momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
81 momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
82 momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
83 momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
84 momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
85 momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
86 momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
87 momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
88 momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
89 momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
90
91
92 // Calculate thrust from "Signal Side"
93 const ROOT::Math::XYZVector thrustB = Thrust::calculateThrust(sig_side_momenta);
94
95 CleoCones myCleoCones(momenta, roe_side_momenta, thrustB, use_all, use_roe);
96
97 const auto& cleo_cone_with_all = myCleoCones.cleo_cone_with_all();
98 const auto& cleo_cone_with_roe = myCleoCones.cleo_cone_with_roe();
99
100 EXPECT_FLOAT_EQ(0.823567, cleo_cone_with_all[0]);
101 EXPECT_FLOAT_EQ(4.7405558, cleo_cone_with_all[1]);
102 EXPECT_FLOAT_EQ(1.7517139, cleo_cone_with_all[2]);
103 EXPECT_FLOAT_EQ(0.37677661, cleo_cone_with_all[3]);
104 EXPECT_FLOAT_EQ(0.622467, cleo_cone_with_all[4]);
105 EXPECT_FLOAT_EQ(0, cleo_cone_with_all[5]);
106 EXPECT_FLOAT_EQ(0, cleo_cone_with_all[6]);
107 EXPECT_FLOAT_EQ(0, cleo_cone_with_all[7]);
108 EXPECT_FLOAT_EQ(0, cleo_cone_with_all[8]);
109
110 EXPECT_FLOAT_EQ(0.823567, cleo_cone_with_roe[0]);
111 EXPECT_FLOAT_EQ(1.9106253, cleo_cone_with_roe[1]);
112 EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[2]);
113 EXPECT_FLOAT_EQ(0.37677661, cleo_cone_with_roe[3]);
114 EXPECT_FLOAT_EQ(0.622467, cleo_cone_with_roe[4]);
115 EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[5]);
116 EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[6]);
117 EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[7]);
118 EXPECT_FLOAT_EQ(0, cleo_cone_with_roe[8]);
119 }

◆ TEST_F() [3/9]

Test the calculation of the Fox-Wolfram moments.

Definition at line 122 of file eventShapeCoreAlgorithm.cc.

123 {
124 std::vector<ROOT::Math::XYZVector> momenta;
125
126 momenta.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
127 momenta.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
128 momenta.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
129 momenta.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
130 momenta.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
131 momenta.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
132 momenta.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
133 momenta.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
134 momenta.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
135 momenta.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
136
137 FoxWolfram FW(momenta);
138 FW.calculateBasicMoments();
139 EXPECT_FLOAT_EQ(0.63011014, FW.getR(2));
140
141 // Tests using the analytic values provided by Fox & Wolfram
142 // for simple distributions
143 // first: Back-to-back versors
144 momenta.clear();
145 momenta.emplace_back(1., 0., 0.);
146 momenta.emplace_back(-1., 0., 0.);
147
148 FW.setMomenta(momenta);
149 FW.calculateBasicMoments();
150
151 EXPECT_TRUE(TMath::Abs(FW.getR(0) - 1) < 0.001);
152 EXPECT_TRUE(TMath::Abs(FW.getR(1)) < 0.001);
153 EXPECT_TRUE(TMath::Abs(FW.getR(2) - 1) < 0.001);
154 EXPECT_TRUE(TMath::Abs(FW.getR(3)) < 0.001);
155 EXPECT_TRUE(TMath::Abs(FW.getR(4) - 1) < 0.001);
156
157 // second: equatorial line
158 momenta.clear();
159 TRandom3 rnd;
160 for (int i = 0; i < 10000; i++) {
161 double phi = rnd.Uniform(0., 2 * TMath::Pi());
162 momenta.emplace_back(TMath::Cos(phi), TMath::Sin(phi), 0.);
163 }
164 FW.setMomenta(momenta);
165 FW.calculateBasicMoments();
166
167 EXPECT_TRUE(TMath::Abs(FW.getR(0) - 1.) < 0.001);
168 EXPECT_TRUE(TMath::Abs(FW.getR(1)) < 0.001);
169 EXPECT_TRUE(TMath::Abs(FW.getR(2) - 0.25) < 0.001);
170 EXPECT_TRUE(TMath::Abs(FW.getR(3)) < 0.001);
171 EXPECT_TRUE(TMath::Abs(FW.getR(4) - 0.14) < 0.001);
172
173 }

◆ TEST_F() [4/9]

Test the calculation of the Harmonic moments.

Definition at line 178 of file eventShapeCoreAlgorithm.cc.

179 {
180
181 float dummySqrtS = 10.;
182 std::vector<XYZVector> partMom;
183
184 partMom.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
185 partMom.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
186 partMom.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
187 partMom.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
188 partMom.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
189 partMom.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
190 partMom.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
191 partMom.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
192 partMom.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
193 partMom.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
194
195 XYZVector axis(0., 0., 1.);
196 // repeats the calculation
197 double moment[9] = {0.};
198 for (auto& p : partMom) {
199 double pMag = p.R();
200 double cTheta = p.Dot(axis) / pMag;
201 for (short i = 0; i < 9; i++)
202 moment[i] += pMag * boost::math::legendre_p(i, cTheta) / dummySqrtS;
203 }
204
205 HarmonicMoments MM(partMom, axis);
206 MM.calculateAllMoments();
207
208 for (short i = 0; i < 9; i++)
209 EXPECT_FLOAT_EQ(moment[i], MM.getMoment(i, dummySqrtS));
210 }

◆ 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.

217 {
218 std::vector<XYZVector> partMom;
219
220 partMom.emplace_back(0.5429965262452898, 0.37010582077332344, 0.0714978744529432);
221 partMom.emplace_back(0.34160659934755344, 0.6444967896760643, 0.18455766323674105);
222 partMom.emplace_back(0.9558442475237068, 0.3628892505037786, 0.545225050633818);
223 partMom.emplace_back(0.8853521332124835, 0.340704481181513, 0.34728211023189237);
224 partMom.emplace_back(0.3155615844988947, 0.8307541128801257, 0.45701302024212986);
225 partMom.emplace_back(0.6100164897524695, 0.5077455724845565, 0.06639458334119974);
226 partMom.emplace_back(0.5078972239903029, 0.9196504908351234, 0.3710366834603026);
227 partMom.emplace_back(0.06252858849289977, 0.4680168989606487, 0.4056055050148607);
228 partMom.emplace_back(0.61672460498333, 0.4472311336875816, 0.31288581834261064);
229 partMom.emplace_back(0.18544654870476218, 0.0758107751704592, 0.31909701462121065);
230
231
232 SphericityEigenvalues Sph(partMom);
233 EXPECT_FLOAT_EQ(0., Sph.getEigenvalue(0));
234 EXPECT_FLOAT_EQ(0., Sph.getEigenvalue(1));
235 EXPECT_FLOAT_EQ(0., Sph.getEigenvalue(2));
236
237 Sph.calculateEigenvalues();
238
239 EXPECT_FLOAT_EQ(0.87272972, Sph.getEigenvalue(0));
240 EXPECT_FLOAT_EQ(0.095489956, Sph.getEigenvalue(1));
241 EXPECT_FLOAT_EQ(0.031780295, Sph.getEigenvalue(2));
242 EXPECT_FLOAT_EQ(1., Sph.getEigenvalue(0) + Sph.getEigenvalue(1) + Sph.getEigenvalue(2));
243 }

◆ TEST_F() [6/9]

TEST_F ( eventShapeCoreAlgorithmTest  ,
Thrust   
)

Test the calculation of a thrust axis.

Definition at line 38 of file eventShapeCoreAlgorithm.cc.

39 {
40 std::vector<ROOT::Math::XYZVector> momenta;
41 // random generated numbers
42 momenta.emplace_back(0.5935352844151847, 0.28902324918117417, 0.9939000705771412);
43 momenta.emplace_back(0.7097025137911714, 0.5118418422879152, 0.44501044145648994);
44 momenta.emplace_back(0.6005771199332856, 0.12366608492454145, 0.7541373665256832);
45 momenta.emplace_back(0.8548902083897467, 0.6887268865943484, 0.34301136659215437);
46 momenta.emplace_back(0.26863521039211535, 0.011148638667487942, 0.96186334199901);
47
48 const ROOT::Math::XYZVector thrustB = Thrust::calculateThrust(momenta);
49
50 EXPECT_FLOAT_EQ(0.925838, thrustB.R());
51 EXPECT_FLOAT_EQ(0.571661, thrustB.X());
52 EXPECT_FLOAT_EQ(0.306741, thrustB.Y());
53 EXPECT_FLOAT_EQ(0.660522, thrustB.Z());
54 }

◆ TEST_F() [7/9]

TEST_F ( PIDPriorsTest  ,
PIDPriorsTableTest   
)

Test of the PIDPriorsTable class.

Definition at line 35 of file PIDPriors.cc.

36 {
37
38 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.};
39 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.};
40
41 PIDPriorsTable PIDTable = PIDPriorsTable();
42 PIDTable.setBinEdges(edgesX, edgesY);
43
44 PIDTable.printPrior();
45
46 PIDTable.setPriorValue(1.1, 0.54, 0.1);
47 PIDTable.setPriorValue(9.5, -0.67, 0.3);
48
49 std::cout << "after setting (1.1, 0.54) amd (9.5, -0.67)" << std::endl;
50 PIDTable.printPrior();
51
52 EXPECT_B2WARNING(PIDTable.setPriorValue(11.0, -10.01, 3.));
53 EXPECT_B2WARNING(PIDTable.setPriorValue(4.23, -0.01, 4.));
54
55 std::cout << " ----- Setting the priors from outside ------ " << std::endl;
56
57 std::vector<float> prior;
58
59 for (int j = 0; j < (int)edgesY.size() - 1; j++) {
60 for (int i = 0; i < (int)edgesX.size() - 1; i++) {
61 prior.push_back(0.001 * (i + j * (edgesX.size() - 1)));
62 }
63 }
64
65
66 PIDTable.setPriorsTable(prior);
67
68 PIDTable.printPrior();
69
70 EXPECT_FLOAT_EQ(PIDTable.getPriorValue(0.1, -0.98), 0.001 * (0 + 0 * (edgesX.size() - 1)));
71 EXPECT_FLOAT_EQ(PIDTable.getPriorValue(0.71, -0.65), 0.001 * (5 + 4 * (edgesX.size() - 1)));
72 EXPECT_FLOAT_EQ(PIDTable.getPriorValue(0.71, -1.65), 0.);
73 EXPECT_FLOAT_EQ(PIDTable.getPriorValue(11.71, 0.876), 0.);
74 EXPECT_B2WARNING(PIDTable.getPriorValue(0.71, -1.65));
75 EXPECT_B2WARNING(PIDTable.getPriorValue(11.71, 0.876));
76 EXPECT_FLOAT_EQ(0., PIDTable.getPriorValue(11.71, 0.876));
77
78
79
80 // now let's try to construct a table with no or on edge only
81 std::vector<float> badEdges1 = {0.};
82 std::vector<float> badEdges2 = {};
83
84
85 PIDPriorsTable PIDTable2 = PIDPriorsTable();
86
87 EXPECT_B2WARNING(PIDTable2.setBinEdges(edgesX, badEdges1));
88 EXPECT_B2WARNING(PIDTable2.setBinEdges(edgesX, badEdges2));
89
90
91 PIDTable2.setBinEdges(edgesX, badEdges2);
92
93 PIDTable2.printPrior();
94
95 PIDTable2.setPriorValue(1.1, 0.54, 0.3);
96 PIDTable2.setPriorValue(9.5, -0.67, 0.3);
97
98 std::cout << "after setting (1.1, 0.54) and (9.5, -0.67)" << std::endl;
99 PIDTable2.printPrior();
100
101 EXPECT_B2WARNING(PIDTable2.setPriorValue(11.0, -10.01, 0.3));
102 EXPECT_FLOAT_EQ(PIDTable2.getPriorValue(9.516, 0.), 0.3);
103
104 }

◆ TEST_F() [8/9]

TEST_F ( PIDPriorsTest  ,
PIDPriorTest   
)

Test the PIDPriors dbobject.

Definition at line 108 of file PIDPriors.cc.

109 {
110
111 // Creates 6 prior objects in different formats (array, PIDPriorsTable, TH2F) to test different setters
112
113 std::vector<float> edgesX = {0., 1, 2, 3, 4, 5, 6, 7, 8, 9, 10.};
114 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.};
115
116 std::vector<float> prior1;
117 std::vector<float> prior2;
118 std::vector<float> prior3;
119 std::vector<float> prior4;
120 std::vector<float> prior5;
121 std::vector<float> prior6;
122
123 TH2F* counts1 = new TH2F("counts1", "counts1", 10, 0, 10., 20, -1, 1);
124 TH2F* counts2 = new TH2F("counts2", "counts2", 10, 0, 10., 20, -1, 1);
125 TH2F* counts3 = new TH2F("counts3", "counts3", 10, 0, 10., 20, -1, 1);
126 TH2F* counts4 = new TH2F("counts4", "counts4", 10, 0, 10., 20, -1, 1);
127 TH2F* counts5 = new TH2F("counts5", "counts5", 10, 0, 10., 20, -1, 1);
128 TH2F* counts6 = new TH2F("counts6", "counts6", 10, 0, 10., 20, -1, 1);
129
130 TH2F* norms = new TH2F("norms", "norms", 10, 0, 10., 20, -1, 1);
131
132 TRandom3* rn = new TRandom3();
133
134
135
136 for (int i = 0; i < (int)edgesX.size() - 1; i++) {
137 for (int j = 0; j < (int)edgesY.size() - 1; j++) {
138
139 float n1 = (float)rn->Poisson(100);
140 float n2 = (float)rn->Poisson(100);
141 float n3 = (float)rn->Poisson(700);
142 float n4 = (float)rn->Poisson(200);
143 float n5 = (float)rn->Poisson(100);
144 float n6 = (float)rn->Poisson(10);
145
146 float n = n1 + n2 + n3 + n4 + n5 + n6;
147
148 prior1.push_back(n1 / n);
149 prior2.push_back(n2 / n);
150 prior3.push_back(n3 / n);
151 prior4.push_back(n4 / n);
152 prior5.push_back(n5 / n);
153 prior6.push_back(n6 / n);
154
155 counts1->SetBinContent(i + 1, j + 1, n1);
156 counts2->SetBinContent(i + 1, j + 1, n2);
157 counts3->SetBinContent(i + 1, j + 1, n3);
158 counts4->SetBinContent(i + 1, j + 1, n4);
159 counts5->SetBinContent(i + 1, j + 1, n5);
160 counts6->SetBinContent(i + 1, j + 1, n6);
161 norms->SetBinContent(i + 1, j + 1, n);
162
163 }
164 }
165
166
167 PIDPriorsTable PIDTable1 = PIDPriorsTable();
168 PIDTable1.setBinEdges(edgesX, edgesY);
169 PIDTable1.setPriorsTable(prior1);
170
171 PIDPriorsTable PIDTable2 = PIDPriorsTable();
172 PIDTable2.setBinEdges(edgesX, edgesY);
173 PIDTable2.setPriorsTable(prior2);
174
175 PIDPriorsTable PIDTable3 = PIDPriorsTable();
176 PIDTable3.setBinEdges(edgesX, edgesY);
177 PIDTable3.setPriorsTable(prior3);
178
179 PIDPriorsTable PIDTable4 = PIDPriorsTable();
180 PIDTable4.setBinEdges(edgesX, edgesY);
181 PIDTable4.setPriorsTable(prior4);
182
183 PIDPriorsTable PIDTable5 = PIDPriorsTable();
184 PIDTable5.setBinEdges(edgesX, edgesY);
185 PIDTable5.setPriorsTable(prior5);
186
187 PIDPriorsTable PIDTable6 = PIDPriorsTable();
188 PIDTable6.setBinEdges(edgesX, edgesY);
189 PIDTable6.setPriorsTable(prior6);
190
191 Const::ChargedStable part_e = Const::ChargedStable(11);
192 Const::ChargedStable part_mu = Const::ChargedStable(13);
193 Const::ChargedStable part_pi = Const::ChargedStable(211);
194 Const::ChargedStable part_K = Const::ChargedStable(321);
195 Const::ChargedStable part_p = Const::ChargedStable(2212);
196 Const::ChargedStable part_d = Const::ChargedStable(1000010020);
197
198
199 std::cout << "setting the tables " << std::endl;
200 // setter with the priortable
201 PIDPriors priors = PIDPriors();
202 priors.setPriors(part_e, PIDTable1);
203 priors.setPriors(part_mu, PIDTable2);
204 priors.setPriors(part_pi, PIDTable3);
205 priors.setPriors(part_K, PIDTable4);
206 priors.setPriors(part_p, PIDTable5);
207 priors.setPriors(part_d, PIDTable6);
208
209 // setter with the count histograms
210 PIDPriors priors2 = PIDPriors();
211 priors2.setPriors(part_e, counts1, norms);
212 priors2.setPriors(part_mu, counts2, norms);
213 priors2.setPriors(part_pi, counts3, norms);
214 priors2.setPriors(part_K, counts4, norms);
215 priors2.setPriors(part_p, counts5, norms);
216 priors2.setPriors(part_d, counts6, norms);
217
218 counts1->Divide(norms);
219 counts2->Divide(norms);
220 counts3->Divide(norms);
221 counts4->Divide(norms);
222 counts5->Divide(norms);
223 counts6->Divide(norms);
224
225 std::cout << "division done" << std::endl;
226 // setter with priors stored in a TH2F
227 std::cout << "setting the TH2 " << std::endl;
228 PIDPriors priors3 = PIDPriors();
229 priors3.setPriors(part_e, counts1);
230 priors3.setPriors(part_mu, counts2);
231 priors3.setPriors(part_pi, counts3);
232 priors3.setPriors(part_K, counts4);
233 priors3.setPriors(part_p, counts5);
234 priors3.setPriors(part_d, counts6);
235
236
237 std::vector<Const::ChargedStable> codes = {part_e, part_mu, part_pi, part_K, part_p, part_d};
238 for (auto code : codes) {
239 for (int i = 0; i < 3; i++) {
240 double p = rn->Uniform(0., 10);
241 double c = rn->Uniform(-1., 1);
242 EXPECT_FLOAT_EQ(priors.getPriorValue(code, p, c), priors2.getPriorValue(code, p, c));
243 EXPECT_FLOAT_EQ(priors.getPriorValue(code, p, c), priors3.getPriorValue(code, p, c));
244 }
245 }
246
247 EXPECT_FLOAT_EQ(0., priors.getPriorValue(part_pi, 11.1, 0.));
248 EXPECT_FLOAT_EQ(0., priors.getPriorValue(part_pi, -1, 0.5));
249 EXPECT_FLOAT_EQ(0., priors.getPriorValue(part_pi, 0.98, -1.001));
250
251 EXPECT_FLOAT_EQ(0., priors2.getPriorValue(part_pi, 11.1, 0.));
252 EXPECT_FLOAT_EQ(0., priors2.getPriorValue(part_pi, -1, 0.5));
253 EXPECT_FLOAT_EQ(0., priors2.getPriorValue(part_pi, 0.98, -1.001));
254
255 EXPECT_FLOAT_EQ(0., priors3.getPriorValue(part_pi, 11.1, 0.));
256 EXPECT_FLOAT_EQ(0., priors3.getPriorValue(part_pi, -1, 0.5));
257 EXPECT_FLOAT_EQ(0., priors3.getPriorValue(part_pi, 0.98, -1.001));
258
259 }
TString rn()
Get random string.
Definition: tools.h:38

◆ 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.

164 {
165
166 PIDDetectorWeights dbrep("tree", m_dummyFile);
167
168 // Test for correct filling of the RDataFrame.
169 auto pdgIds = dbrep.getWeightsRDF().Take<double>("pdgId").GetValue();
170 for (const auto& pdgId : pdgIds) {
171 EXPECT_EQ(pdgId, m_testHypo.getPDGCode());
172 }
173
174 // Retrieve weight for a (p, theta) pair in the available weights range.
175 auto p = 1.23; // GeV/c
176 auto theta = 1.34; // rad
177 auto weight = dbrep.getWeight(m_testHypo, m_detector, p, theta);
178 EXPECT_EQ(weight, -0.118);
179
180 // Test for weight in case of out-of-range p and/or theta values.
181 p = 1.46;
182 theta = 0.12;
183 weight = dbrep.getWeight(m_testHypo, m_detector, p, theta);
184 EXPECT_TRUE(std::isnan(weight));
185
186 // Trigger a FATAL if reading an ill-defined source table.
187 EXPECT_B2FATAL(PIDDetectorWeights("tree_broken", m_dummyFile));
188
189 }