Belle II Software light-2406-ragdoll
ParticleBase.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * External Contributor: Wouter Hulsbergen *
5 * *
6 * See git log for contributors and copyright holders. *
7 * This file is licensed under LGPL-3.0, see LICENSE.md. *
8 **************************************************************************/
9
10//Creates a particle. Base class for all other particle classes.
11
12#include <TDatabasePDG.h>
13
14#include <analysis/VertexFitting/TreeFitter/ParticleBase.h>
15#include <analysis/VertexFitting/TreeFitter/InternalParticle.h>
16#include <analysis/VertexFitting/TreeFitter/Composite.h>
17#include <analysis/VertexFitting/TreeFitter/RecoResonance.h>
18#include <analysis/VertexFitting/TreeFitter/RecoTrack.h>
19
20#include <analysis/VertexFitting/TreeFitter/RecoPhoton.h>
21#include <analysis/VertexFitting/TreeFitter/RecoKlong.h>
22#include <analysis/VertexFitting/TreeFitter/Resonance.h>
23#include <analysis/VertexFitting/TreeFitter/Origin.h>
24#include <analysis/VertexFitting/TreeFitter/FitParams.h>
25
26#include <framework/geometry/BFieldManager.h>
27
28namespace TreeFitter {
29
31 m_particle(particle),
32 m_mother(mother),
33 m_isStronglyDecayingResonance(false),
34 m_config(config),
35 m_index(0),
36 m_name("Unknown")
37 {
38 if (particle) {
40 const int pdgcode = particle->getPDGCode();
41 if (pdgcode) { // PDG code != 0
43 }
44 }
45 }
46
47 ParticleBase::ParticleBase(const std::string& name) :
48 m_particle(nullptr),
49 m_mother(nullptr),
50 m_isStronglyDecayingResonance(false),
51 m_config(nullptr),
52 m_index(0),
53 m_name(name)
54 {}
55
56
58 {
59 for (auto daughter : m_daughters) {
60 delete daughter;
61 }
62 m_daughters.clear();
63 }
64
66 {
67 auto newDaughter = ParticleBase::createParticle(cand, this, config, forceFitAll);
68 m_daughters.push_back(newDaughter);
69 return m_daughters.back();
70 }
71
72
74 {
75 auto iter = std::find(m_daughters.begin(), m_daughters.end(), pb);
76 if (iter != m_daughters.end()) {
77 delete *iter;
78 m_daughters.erase(iter);
79 } else {
80 B2ERROR("Cannot remove particle, because not found ...");
81 }
82 }
83
84 void ParticleBase::updateIndex(int& offset)
85 {
86 for (auto* daughter : m_daughters) {
87 daughter->updateIndex(offset);
88 }
89 m_index = offset;
90 offset += dim();
91 }
92
94 Belle2::Particle* daughter,
95 const ConstraintConfiguration& config,
96 bool forceFitAll
97 )
98 {
99 return new Origin(daughter, config, forceFitAll);
100 }
101
103 const ConstraintConfiguration& config, bool forceFitAll)
104 {
105 ParticleBase* rc = nullptr;
106
107 if (!mother) { // If there is no mother, this is the 'head of tree' particle (is never a resonance)
108 rc = new InternalParticle(particle, nullptr, config, forceFitAll);
109 } else if (particle->hasExtraInfo("bremsCorrected") // Has Bremsstrahlungs-recovery
110 && particle->getExtraInfo("bremsCorrected") != 0) { // and gammas are attached
111 rc = new Composite(particle, mother, config, true);
112 // if no gamma is attached, it is treated as a RecoTrack
113 } else if (particle->hasExtraInfo("treeFitterTreatMeAsInvisible")
114 && particle->getExtraInfo("treeFitterTreatMeAsInvisible") == 1) { // dummy particles with invisible flag
115 rc = new RecoResonance(particle, mother, config);
116 } else if (particle->getTrack()) { // external reconstructed track
117 rc = new RecoTrack(particle, mother);
118 } else if (particle->getECLCluster()) { // external reconstructed photon
119 rc = new RecoPhoton(particle, mother);
120 } else if (particle->getKLMCluster()) { // external reconstructed klong
121 rc = new RecoKlong(particle, mother);
122 } else if (particle->getMdstArrayIndex()) { // external composite e.g. V0
123 rc = new InternalParticle(particle, mother, config, forceFitAll);
124 } else { // 'internal' particles
125 if (isAResonance(particle)) {
126 rc = new Resonance(particle, mother, config, forceFitAll);
127 } else {
128 rc = new InternalParticle(particle, mother, config, forceFitAll);
129 }
130 }
131 return rc;
132 }
133
135 {
136 bool rc = false ;
137 const int pdgcode = std::abs(particle->getPDGCode());
138
139 if (pdgcode && !(particle->getMdstArrayIndex())) {
140 switch (pdgcode) {
141 case 22: //photon conversion
142 rc = false;
143 break ;
144
145 case -11: //bremsstrahlung
146 case 11:
147 rc = true ;
148 break ;
149 default: //everything with boosted flight length less than 1 micrometer
150 rc = (pdgcode && particle->getPDGLifetime() < 1e-14);
151 }
152 }
153 return rc ;
154 }
155
156 void ParticleBase::collectVertexDaughters(std::vector<ParticleBase*>& particles, int posindex)
157 {
158 if (mother() && mother()->posIndex() == posindex) {
159 particles.push_back(this);
160 }
161
162 for (auto* daughter : m_daughters) {
163 daughter->collectVertexDaughters(particles, posindex);
164 }
165 }
166
168 {
169 // this is very sensitive and can heavily affect the efficiency of the fit
170 ErrCode status;
171
172 const int posindex = posIndex();
173 if (posindex >= 0) {
174 for (int i = 0; i < 3; ++i) {
175 fitparams.getCovariance()(posindex + i, posindex + i) = 1;
176 }
177 }
178
179 const int momindex = momIndex();
180 if (momindex >= 0) {
181 const int maxrow = hasEnergy() ? 4 : 3;
182 for (int i = 0; i < maxrow; ++i) {
183 fitparams.getCovariance()(momindex + i, momindex + i) = 10.;
184 }
185 }
186
187 const int tauindex = tauIndex();
188 if (tauindex >= 0) {
189 fitparams.getCovariance()(tauindex, tauindex) = 1.;
190 }
191 return status;
192 }
193
194 std::string ParticleBase::parname(int thisindex) const
195 {
196 std::string rc = m_name;
197 switch (thisindex) {
198 case 0: rc += "_x "; break;
199 case 1: rc += "_y "; break;
200 case 2: rc += "_z "; break;
201 case 3: rc += "_tau"; break;
202 case 4: rc += "_px "; break;
203 case 5: rc += "_py "; break;
204 case 6: rc += "_pz "; break;
205 case 7: rc += "_E "; break;
206 default:;
207 }
208 return rc;
209 }
210
212 {
213 const ParticleBase* rc = (m_particle == particle) ? this : nullptr;
214 if (!rc) {
215 for (auto* daughter : m_daughters) {
216 rc = daughter->locate(particle);
217 if (rc) {break;}
218 }
219 }
220 return rc;
221 }
222
224 {
225 map.push_back(std::pair<const ParticleBase*, int>(this, index()));
226 for (auto* daughter : m_daughters) {
227 daughter->retrieveIndexMap(map);
228 }
229 }
230
231 double ParticleBase::chiSquare(const FitParams& fitparams) const
232 {
233 double rc = 0;
234 for (auto* daughter : m_daughters) {
235 rc += daughter->chiSquare(fitparams);
236 }
237 return rc;
238 }
239
241 {
242 int rc = 0;
243 for (auto* daughter : m_daughters) {
244 rc += daughter->nFinalChargedCandidates();
245 }
246 return rc;
247 }
248
250 {
251 assert(m_config);
252 // only allow 2d for head of tree particles that are beam constrained
253 const int dim = m_config->m_originDimension == 2 && std::abs(m_particle->getPDGCode()) == m_config->m_headOfTreePDG ? 2 : 3;
254 const int posindexmother = mother()->posIndex();
255 const int posindex = posIndex();
256 const int tauindex = tauIndex();
257 const int momindex = momIndex();
258
259 const double tau = fitparams.getStateVector()(tauindex);
260 Eigen::Matrix < double, 1, -1, 1, 1, 3 > x_vec = fitparams.getStateVector().segment(posindex, dim);
261 Eigen::Matrix < double, 1, -1, 1, 1, 3 > x_m = fitparams.getStateVector().segment(posindexmother, dim);
262 Eigen::Matrix < double, 1, -1, 1, 1, 3 > p_vec = fitparams.getStateVector().segment(momindex, dim);
263 const double mom = p_vec.norm();
264 const double mom3 = mom * mom * mom;
265
266 if (3 == dim) {
267 // we can already set these
268 //diagonal momentum
269 p.getH()(0, momindex) = tau * (p_vec(1) * p_vec(1) + p_vec(2) * p_vec(2)) / mom3 ;
270 p.getH()(1, momindex + 1) = tau * (p_vec(0) * p_vec(0) + p_vec(2) * p_vec(2)) / mom3 ;
271 p.getH()(2, momindex + 2) = tau * (p_vec(0) * p_vec(0) + p_vec(1) * p_vec(1)) / mom3 ;
272
273 //offdiagonal momentum
274 p.getH()(0, momindex + 1) = - tau * p_vec(0) * p_vec(1) / mom3 ;
275 p.getH()(0, momindex + 2) = - tau * p_vec(0) * p_vec(2) / mom3 ;
276
277 p.getH()(1, momindex + 0) = - tau * p_vec(1) * p_vec(0) / mom3 ;
278 p.getH()(1, momindex + 2) = - tau * p_vec(1) * p_vec(2) / mom3 ;
279
280 p.getH()(2, momindex + 0) = - tau * p_vec(2) * p_vec(0) / mom3 ;
281 p.getH()(2, momindex + 1) = - tau * p_vec(2) * p_vec(1) / mom3 ;
282
283 } else if (2 == dim) {
284
285 // NOTE THAT THESE ARE DIFFERENT IN 2d
286 p.getH()(0, momindex) = tau * (p_vec(1) * p_vec(1)) / mom3 ;
287 p.getH()(1, momindex + 1) = tau * (p_vec(0) * p_vec(0)) / mom3 ;
288
289 //offdiagonal momentum
290 p.getH()(0, momindex + 1) = - tau * p_vec(0) * p_vec(1) / mom3 ;
291 p.getH()(1, momindex + 0) = - tau * p_vec(1) * p_vec(0) / mom3 ;
292 } else {
293 B2FATAL("Dimension of Geometric constraint is not 2 or 3. This will crash many things. You should feel bad.");
294 }
295
296 for (int row = 0; row < dim; ++row) {
297
298 double posxmother = x_m(row);
299 double posx = x_vec(row);
300 double momx = p_vec(row);
301
305 p.getResiduals()(row) = posxmother + tau * momx / mom - posx ;
306 p.getH()(row, posindexmother + row) = 1;
307 p.getH()(row, posindex + row) = -1;
308 p.getH()(row, tauindex) = momx / mom;
309 }
310
311 return ErrCode(ErrCode::Status::success);
312 }
313
315 Projection& p) const
316 {
317 const double mass = particle()->getPDGMass();
318 const double mass2 = mass * mass;
319 double px = 0;
320 double py = 0;
321 double pz = 0;
322 double E = 0;
323
324 // the parameters of the daughters must be used otherwise the mass constraint does not have an effect on the extracted daughter momenta
325 for (const auto* daughter : m_daughters) {
326 const int momindex = daughter->momIndex();
327 // in most cases the daughters will be final states so we cache the value to use it in the energy column
328 const double px_daughter = fitparams.getStateVector()(momindex);
329 const double py_daughter = fitparams.getStateVector()(momindex + 1);
330 const double pz_daughter = fitparams.getStateVector()(momindex + 2);
331
332 px += px_daughter;
333 py += py_daughter;
334 pz += pz_daughter;
335 if (daughter->hasEnergy()) {
336 E += fitparams.getStateVector()(momindex + 3);
337 } else {
338 // final states dont have an energy index
339 const double m = daughter->particle()->getPDGMass();
340 E += std::sqrt(m * m + px_daughter * px_daughter + py_daughter * py_daughter + pz_daughter * pz_daughter);
341 }
342 }
343
347 p.getResiduals()(0) = mass2 - E * E + px * px + py * py + pz * pz;
348
349 for (const auto* daughter : m_daughters) {
350 //dr/dx = d/dx m2-{E1+E2+...}^2+{p1+p2+...}^2 = 2*x (x= E or p)
351 const int momindex = daughter->momIndex();
352 p.getH()(0, momindex) = 2.0 * px;
353 p.getH()(0, momindex + 1) = 2.0 * py;
354 p.getH()(0, momindex + 2) = 2.0 * pz;
355
356 if (daughter->hasEnergy()) {
357 p.getH()(0, momindex + 3) = -2.0 * E;
358 } else {
359 const double px_daughter = fitparams.getStateVector()(momindex);
360 const double py_daughter = fitparams.getStateVector()(momindex + 1);
361 const double pz_daughter = fitparams.getStateVector()(momindex + 2);
362 const double m = daughter->particle()->getPDGMass();
363
364 const double E_daughter = std::sqrt(m * m + px_daughter * px_daughter + py_daughter * py_daughter + pz_daughter * pz_daughter);
365 const double E_by_E_daughter = E / E_daughter;
366 p.getH()(0, momindex) -= 2.0 * E_by_E_daughter * px_daughter;
367 p.getH()(0, momindex + 1) -= 2.0 * E_by_E_daughter * py_daughter;
368 p.getH()(0, momindex + 2) -= 2.0 * E_by_E_daughter * pz_daughter;
369 }
370
371 }
372 return ErrCode(ErrCode::Status::success);
373 }
374
376 Projection& p) const
377 {
378 const double mass = particle()->getPDGMass();
379 const double mass2 = mass * mass;
380 const int momindex = momIndex();
381 const double px = fitparams.getStateVector()(momindex);
382 const double py = fitparams.getStateVector()(momindex + 1);
383 const double pz = fitparams.getStateVector()(momindex + 2);
384 const double E = fitparams.getStateVector()(momindex + 3);
385
389 p.getResiduals()(0) = mass2 - E * E + px * px + py * py + pz * pz;
390
391 p.getH()(0, momindex) = 2.0 * px;
392 p.getH()(0, momindex + 1) = 2.0 * py;
393 p.getH()(0, momindex + 2) = 2.0 * pz;
394 p.getH()(0, momindex + 3) = -2.0 * E;
395
396 // TODO 0 in most cases -> needs special treatment if width=0 to not crash chi2 calculation
397 // const double width = TDatabasePDG::Instance()->GetParticle(particle()->getPDGCode())->Width();
398 // transport measurement uncertainty into residual system
399 // f' = sigma_x^2 * (df/dx)^2
400 // p.getV()(0) = width * width * 4 * mass2;
401
402 return ErrCode(ErrCode::Status::success);
403 }
404
406 Projection& p) const
407 {
408 assert(m_config);
409 if (m_config->m_massConstraintType == 0) {
410 return projectMassConstraintParticle(fitparams, p);
411 } else {
412 return projectMassConstraintDaughters(fitparams, p);
413 }
414 }
415
417 {
418 if (type == Constraint::mass) {
419 return projectMassConstraint(fitparams, p);
420 } else {
421 B2FATAL("Trying to project constraint of ParticleBase type. This is undefined.");
422 }
423 return ErrCode(ErrCode::Status::badsetup);
424 }
425
427 {
428 const int tauindex = tauIndex();
429 if (tauindex >= 0 && hasPosition()) {
430
431 const int posindex = posIndex();
432 const int mother_ps_index = mother()->posIndex();
433 const int dim = m_config->m_originDimension; // TODO can we configure this to be particle specific?
434
435 // tau has different meaning depending on the dimension of the constraint
436 // 2-> use x-y projection
437 const Eigen::Matrix < double, 1, -1, 1, 1, 3 > vertex_dist =
438 fitparams.getStateVector().segment(posindex, dim) - fitparams.getStateVector().segment(mother_ps_index, dim);
439 const Eigen::Matrix < double, 1, -1, 1, 1, 3 >
440 mom = fitparams.getStateVector().segment(posindex, dim) - fitparams.getStateVector().segment(mother_ps_index, dim);
441
442 // if an intermediate vertex is not well defined by a track or so it will be initialised with 0
443 // same for the momentum of for example B0, it might be initialised with 0
444 // in those cases use pdg value
445 const double mom_norm = mom.norm();
446 const double dot = std::abs(vertex_dist.dot(mom));
447 const double tau = dot / mom_norm;
448 if (0 == mom_norm || 0 == dot) {
449 const double mass = m_particle->getPDGMass();
450 if (mass > 0)
451 fitparams.getStateVector()(tauindex) = m_particle->getPDGLifetime() * 1e9 * Belle2::Const::speedOfLight / mass;
452 else
453 fitparams.getStateVector()(tauindex) = 0;
454 } else {
455 fitparams.getStateVector()(tauindex) = tau;
456 }
457 }
458
459 return ErrCode(ErrCode::Status::success);
460 }
461} //end namespace TreeFitter
462
static const double speedOfLight
[cm/ns]
Definition: Const.h:695
Class to store reconstructed particles.
Definition: Particle.h:75
const KLMCluster * getKLMCluster() const
Returns the pointer to the KLMCluster object that was used to create this Particle (ParticleType == c...
Definition: Particle.cc:926
const Track * getTrack() const
Returns the pointer to the Track object that was used to create this Particle (ParticleType == c_Trac...
Definition: Particle.cc:845
const ECLCluster * getECLCluster() const
Returns the pointer to the ECLCluster object that was used to create this Particle (if ParticleType =...
Definition: Particle.cc:891
std::string getName() const override
Return name of this particle.
Definition: Particle.cc:1165
bool hasExtraInfo(const std::string &name) const
Return whether the extra info with the given name is set.
Definition: Particle.cc:1266
unsigned getMdstArrayIndex(void) const
Returns 0-based index of MDST store array object (0 for composite particles)
Definition: Particle.h:487
int getPDGCode(void) const
Returns PDG code.
Definition: Particle.h:454
double getPDGMass(void) const
Returns uncertainty on the invariant mass (requires valid momentum error matrix)
Definition: Particle.cc:604
double getPDGLifetime() const
Returns particle nominal lifetime.
Definition: Particle.cc:613
double getExtraInfo(const std::string &name) const
Return given value if set.
Definition: Particle.cc:1289
A class for composite particles, where the daughters must be ignored by the fitter.
Definition: Composite.h:17
constraint configuration class
const int m_originDimension
dimension of the origin constraint and ALL geometric gcosntraints
int m_headOfTreePDG
PDG code of the head particle.
const bool m_massConstraintType
const flag for the type of the mass constraint
Type
type of constraints the order of these constraints is important: it is the order in which they are ap...
Definition: Constraint.h:27
abstract errorocode be aware that the default is success
Definition: ErrCode.h:14
Class to store and manage fitparams (statevector)
Definition: FitParams.h:20
Eigen::Matrix< double, -1, 1, 0, MAX_MATRIX_SIZE, 1 > & getStateVector()
getter for the fit parameters/statevector
Definition: FitParams.h:65
Eigen::Matrix< double, -1, -1, 0, MAX_MATRIX_SIZE, MAX_MATRIX_SIZE > & getCovariance()
getter for the states covariance
Definition: FitParams.h:53
another unnecessary layer of abstraction
representation of the beamspot as a particle
Definition: Origin.h:19
base class for all particles
Definition: ParticleBase.h:25
virtual void updateIndex(int &offset)
this sets the index for momentum, position, etc.
Definition: ParticleBase.cc:84
virtual ErrCode projectMassConstraintParticle(const FitParams &, Projection &) const
project mass constraint using the particles parameters
Belle2::Particle * particle() const
get basf2 particle
Definition: ParticleBase.h:92
ErrCode initTau(FitParams &par) const
initialises tau as a length
virtual ParticleBase * addDaughter(Belle2::Particle *, const ConstraintConfiguration &config, bool forceFitAll=false)
add daughter
Definition: ParticleBase.cc:65
static ParticleBase * createOrigin(Belle2::Particle *daughter, const ConstraintConfiguration &config, bool forceFitAll)
create a custom origin particle or a beamspot
Definition: ParticleBase.cc:93
static ParticleBase * createParticle(Belle2::Particle *particle, const ParticleBase *mother, const ConstraintConfiguration &config, bool forceFitAll=false)
create the according treeFitter particle obj for a basf2 particle type
virtual int dim() const =0
get dimension of constraint
virtual int nFinalChargedCandidates() const
number of charged candidates
virtual void retrieveIndexMap(indexmap &anindexmap) const
get index map
virtual ErrCode projectMassConstraint(const FitParams &, Projection &) const
project mass constraint abstract
const ConstraintConfiguration * m_config
has all the constraint config
Definition: ParticleBase.h:204
virtual std::string parname(int index) const
get name of parameter i
static bool isAResonance(Belle2::Particle *particle)
controls if a particle is treated as a resonance(lifetime=0) or a particle that has a finite lifetime...
virtual ErrCode projectConstraint(Constraint::Type, const FitParams &, Projection &) const
project constraint.
virtual ErrCode projectMassConstraintDaughters(const FitParams &, Projection &) const
project mass constraint using the parameters of the daughters
const ParticleBase * locate(Belle2::Particle *particle) const
get particle base from basf2 particle
virtual void removeDaughter(const ParticleBase *pb)
remove daughter
Definition: ParticleBase.cc:73
bool m_isStronglyDecayingResonance
decay length less than 1 micron
Definition: ParticleBase.h:201
virtual ErrCode initCovariance(FitParams &) const
init covariance matrix
virtual int type() const =0
get particle type
Belle2::Particle * m_particle
pointer to framework type
Definition: ParticleBase.h:192
ParticleBase(Belle2::Particle *particle, const ParticleBase *mother, const ConstraintConfiguration *config=nullptr)
default constructor
Definition: ParticleBase.cc:30
virtual ErrCode projectGeoConstraint(const FitParams &, Projection &) const
project geometrical constraint
virtual double chiSquare(const FitParams &) const
get chi2
virtual int posIndex() const
get vertex index (in statevector!)
Definition: ParticleBase.h:122
std::vector< std::pair< const ParticleBase *, int > > indexmap
alias
Definition: ParticleBase.h:55
virtual int momIndex() const
get momentum index
Definition: ParticleBase.h:128
int index() const
get index
Definition: ParticleBase.h:95
virtual bool hasEnergy() const
get momentum dimension
Definition: ParticleBase.h:132
std::vector< ParticleBase * > m_daughters
daughter container
Definition: ParticleBase.h:198
std::string m_name
name
Definition: ParticleBase.h:211
virtual ~ParticleBase()
destructor, actually does something
Definition: ParticleBase.cc:57
void collectVertexDaughters(std::vector< ParticleBase * > &particles, int posindex)
get vertex daughters
virtual int tauIndex() const
get tau index
Definition: ParticleBase.h:125
const ParticleBase * mother() const
getMother() / hasMother()
Definition: ParticleBase.h:98
virtual bool hasPosition() const
get false
Definition: ParticleBase.h:135
class to store the projected residuals and the corresponding jacobian as well as the covariance matri...
Definition: Projection.h:18
representation of the Klong constraint
Definition: RecoKlong.h:16
representation of the photon constraint
Definition: RecoPhoton.h:16
A class for resonances.
Definition: RecoResonance.h:17
representation of all charged final states FIXME rename since this name is taken in tracking
Definition: RecoTrack.h:18
class for resonances as internal particles
Definition: Resonance.h:17