Belle II Software light-2406-ragdoll
RaveVertexFitter.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9
10#include <analysis/VertexFitting/RaveInterface/RaveVertexFitter.h>
11#include <analysis/VertexFitting/RaveInterface/RaveSetup.h>
12#include <mdst/dataobjects/Track.h>
13
14#include <rave/VertexFactory.h>
15
16//root
17#include <Math/ProbFunc.h>
18#include <Math/Vector3D.h>
19#include <Math/Vector4D.h>
20//stl
21using std::string;
22#include <vector>
23using std::vector;
24
25using namespace Belle2;
26using namespace analysis;
27
28
29
31{
32 //B2WARNING( "RaveVertexFitter::RaveVertexFitter()" );
33 if (RaveSetup::getRawInstance() == nullptr) {
34 B2FATAL("RaveSetup::initialize was not called. It has to be called before RaveSetup or RaveVertexFitter are used");
35 }
36 //B2WARNING "m_useBeamSpot " << m_useBeamSpot );
38 //B2WARNING("m_useBeamSpot " << m_useBeamSpot );
39}
40
42{
43 m_useBeamSpot = false;
44 if (RaveSetup::getRawInstance() == nullptr) {
45 B2FATAL("RaveSetup::initialize was not called. It has to be called before RaveSetup or RaveVertexFitter are used");
46 }
48}
49
50
51
53
54
55void RaveVertexFitter::addTrack(const TrackFitResult* const aTrackPtr)
56{
57 m_raveTracks.push_back(TrackFitResultToRaveTrack(aTrackPtr));
58}
59
60rave::Track RaveVertexFitter::TrackFitResultToRaveTrack(const TrackFitResult* const aTrackPtr) const
61{
62 const int id = m_raveTracks.size();
63
64 B2Vector3D pos = aTrackPtr->getPosition();
65 B2Vector3D mom = aTrackPtr->getMomentum();
66 TMatrixF cov(aTrackPtr->getCovariance6());
67
68
69 // state
70 rave::Vector6D ravestate(pos.X(), pos.Y(), pos.Z(),
71 mom.X(), mom.Y(), mom.Z());
72
73 rave::Covariance6D ravecov(cov(0, 0), cov(1, 0), cov(2, 0),
74 cov(1, 1), cov(2, 1), cov(2, 2),
75 cov(3, 0), cov(4, 0), cov(5, 0),
76 cov(3, 1), cov(4, 1), cov(5, 1),
77 cov(3, 2), cov(4, 2), cov(5, 2),
78 cov(3, 3), cov(4, 3), cov(5, 3),
79 cov(4, 4), cov(5, 4), cov(5, 5));
80
81 return rave::Track(id, ravestate, ravecov, rave::Charge(aTrackPtr->getChargeSign()), 1,
82 1); //the two 1s are just dummy values. They are not used by Rave anyway
83
84}
85
86
87void RaveVertexFitter::addTrack(const Particle* aParticlePtr)
88{
89 const int id = m_raveTracks.size();
90 const ROOT::Math::XYZVector& pos = aParticlePtr->getVertex();
91 const ROOT::Math::XYZVector& mom = aParticlePtr->getMomentum();
92 rave::Vector6D ravestate(pos.X(), pos.Y(), pos.Z(), mom.X(), mom.Y(), mom.Z());
93
94 const TMatrixFSym& cov = aParticlePtr->getMomentumVertexErrorMatrix();
95
96 rave::Covariance6D ravecov(cov(4, 4), cov(4, 5), cov(4, 6),
97 cov(5, 5), cov(5, 6), cov(6, 6),
98 cov(0, 4), cov(1, 4), cov(2, 4),
99 cov(0, 5), cov(1, 5), cov(2, 5),
100 cov(0, 6), cov(1, 6), cov(2, 6),
101 cov(0, 0), cov(0, 1), cov(0, 2),
102 cov(1, 1), cov(1, 2), cov(2, 2));
103
104 // 1 and 1 are dummy values for chi2 and ndf. the are not used for the vertex fit
105 m_raveTracks.emplace_back(id, ravestate, ravecov, rave::Charge(aParticlePtr->getCharge() + 0.1), 1, 1);
106
107 m_belleDaughters.push_back(const_cast<Particle*>(aParticlePtr));
108}
109
110void RaveVertexFitter::addMother(const Particle* aMotherParticlePtr)
111{
112 vector<Particle*> daughters = aMotherParticlePtr->getDaughters();
113
114 int nDaughters = daughters.size();
115 for (int i = 0; i not_eq nDaughters; ++i) {
116 addTrack(daughters[i]);
117 }
118
119}
120
121int RaveVertexFitter::fit(string options)
122{
123 //B2WARNING("RaveVertexFitter::fit(string options)" );
124 //B2WARNING("m_useBeamSpot " << m_useBeamSpot );
125 if (options == "default") {
126 options = "kalman";
127 }
128 int ndf = 0;
129
130 ndf = 2 * m_raveTracks.size();
131
132 if (m_useBeamSpot == true) {
133 ndf += 3;
134 }
135 if (ndf < 4) {
136 return -1;
137 }
138 int nOfVertices = -100;
139
140 if (m_useBeamSpot == true) {
142 const TMatrixDSym& bsCov = RaveSetup::getRawInstance()->m_beamSpotCov;
143 const rave::Covariance3D bsCovRave(bsCov(0, 0), bsCov(0, 1), bsCov(0, 2), bsCov(1, 1), bsCov(1, 2), bsCov(2, 2));
144 RaveSetup::getRawInstance()->m_raveVertexFactory->setBeamSpot(rave::Ellipsoid3D(rave::Point3D(bsPos.X(), bsPos.Y(), bsPos.Z()),
145 bsCovRave));
146 }
147 //B2WARNING( "now fitting with m_raveVertexFactory" );
148 RaveSetup::getRawInstance()->m_raveVertexFactory->setDefaultMethod(options);
150 nOfVertices = m_raveVertices.size();
151
152 return nOfVertices;
153}
154
156{
157
158 if (m_raveVertices.empty()) {
159 B2ERROR("There is no fitted Vertex. Maybe you did not call fit() or maybe the fit was not successful");
160 throw;
161 }
162 if (vertexId >= m_raveVertices.size()) {
163 B2ERROR("The Vertex id " << vertexId << " does not correspond to any fitted vertex");
164 throw;
165 }
166
167}
168
170{
171 isVertexIdValid(vertexId);
172
173 return B2Vector3D(m_raveVertices[vertexId].position().x(), m_raveVertices[vertexId].position().y(),
174 m_raveVertices[vertexId].position().z());
175
176}
177
178double RaveVertexFitter::getWeight(int trackId, VecSize vertexId)const
179{
180 isVertexIdValid(vertexId);
181
182
183 const std::vector < std::pair < float, rave::Track > >& weightedTracks = m_raveVertices[vertexId].weightedTracks();
184 for (unsigned int i = 0; i not_eq weightedTracks.size(); ++i) {
185 if (weightedTracks[i].second.id() == trackId) {
186// B2WARNING( "returning weight for track with x coord: " <<weightedTracks[i].second.state().x() );
187 return weightedTracks[i].first;
188 }
189 }
190 return -1;
191
192}
193
195{
196
197 const std::vector < std::pair < float, rave::Track > >& weightedTracks = m_raveVertices[vertexId].weightedTracks();
198 const int n = weightedTracks.size();
199 vector<int> trackIds(n);
200 for (int i = 0; i not_eq n; ++i) {
201 trackIds[i] = weightedTracks[i].second.id();
202 }
203 return trackIds;
204}
205
207{
208 isVertexIdValid(vertexId);
209
210 return ROOT::Math::chisquared_cdf_c(m_raveVertices[vertexId].chiSquared(), m_raveVertices[vertexId].ndf());
211
212}
213
214double RaveVertexFitter::getNdf(VecSize vertexId) const
215{
216 isVertexIdValid(vertexId);
217
218 return m_raveVertices[vertexId].ndf();
219
220}
221
222double RaveVertexFitter::getChi2(VecSize vertexId) const
223{
224 isVertexIdValid(vertexId);
225
226 return m_raveVertices[vertexId].chiSquared();
227
228}
229
230TMatrixDSym RaveVertexFitter::getCov(VecSize vertexId) const
231{
232 isVertexIdValid(vertexId);
233
234 TMatrixDSym Cov(3);
235 Cov(0, 0) = m_raveVertices[vertexId].error().dxx();
236 Cov(0, 1) = m_raveVertices[vertexId].error().dxy();
237 Cov(0, 2) = m_raveVertices[vertexId].error().dxz();
238 Cov(1, 1) = m_raveVertices[vertexId].error().dyy();
239 Cov(1, 2) = m_raveVertices[vertexId].error().dyz();
240 Cov(2, 2) = m_raveVertices[vertexId].error().dzz();
241 Cov(1, 0) = Cov(0, 1);
242 Cov(2, 1) = Cov(1, 2);
243 Cov(2, 0) = Cov(0, 2);
244
245 return Cov;
246
247}
248
249
251{
252 if (m_raveVertices.size() != 1) {
253 B2ERROR("RaveVertexFitter: Daughters update works only with a single vertex");
254 return;
255 }
256
257 if (!m_raveVertices[0].hasRefittedTracks()) {
258 B2WARNING("RaveVertexFitter: Fitted vertex has no refitted tracks");
259 return;
260 }
261
262 std::vector < rave::Track > rTracks = m_raveVertices[0].tracks(); //< the original tracks
263
264 for (unsigned int i = 0; i < rTracks.size(); i++) {
265 rave::Track rtrk = m_raveVertices[0].refittedTrack(rTracks[i]);
266 const rave::Point3D fittedV = rtrk.position();
267 const rave::Vector3D fittedP = rtrk.momentum();
268 const rave::Covariance6D& fittedCov = rtrk.error();
269
270 ROOT::Math::XYZVector x3(fittedV.x(), fittedV.y(), fittedV.z());
271 double fittedE = sqrt(fittedP.mag2() + (m_belleDaughters[i]->getMass() * (m_belleDaughters[i]->getMass())));
272 ROOT::Math::PxPyPzEVector p4(fittedP.x(), fittedP.y(), fittedP.z(), fittedE);
273
274 TMatrixDSym fitted7CovPart = m_belleDaughters[i]->getMomentumVertexErrorMatrix() ;
275
276 //px,py,pz,E,x,y,z
277 TMatrixDSym fitted7CovM(7);
278 fitted7CovM(0, 0) = fittedCov.dpxpx(); fitted7CovM(0, 0) = fitted7CovM(0, 0);
279 fitted7CovM(0, 1) = fittedCov.dpxpy(); fitted7CovM(1, 0) = fitted7CovM(0, 1);
280 fitted7CovM(0, 2) = fittedCov.dpxpz(); fitted7CovM(2, 0) = fitted7CovM(0, 2);
281 fitted7CovM(0, 3) = fitted7CovPart(0, 3); fitted7CovM(3, 0) = fitted7CovM(0, 3);
282 fitted7CovM(0, 4) = fittedCov.dxpx(); fitted7CovM(4, 0) = fitted7CovM(0, 4);
283 fitted7CovM(0, 5) = fittedCov.dypx(); fitted7CovM(5, 0) = fitted7CovM(0, 5);
284 fitted7CovM(0, 6) = fittedCov.dzpx(); fitted7CovM(6, 0) = fitted7CovM(0, 6);
285
286 fitted7CovM(1, 1) = fittedCov.dpypy(); fitted7CovM(1, 1) = fitted7CovM(1, 1);
287 fitted7CovM(1, 2) = fittedCov.dpypz(); fitted7CovM(2, 1) = fitted7CovM(1, 2);
288 fitted7CovM(1, 3) = fitted7CovPart(1, 3); fitted7CovM(3, 1) = fitted7CovM(1, 3);
289 fitted7CovM(1, 4) = fittedCov.dxpy(); fitted7CovM(4, 1) = fitted7CovM(1, 4);
290 fitted7CovM(1, 5) = fittedCov.dypy(); fitted7CovM(5, 1) = fitted7CovM(1, 5);
291 fitted7CovM(1, 6) = fittedCov.dzpy(); fitted7CovM(6, 1) = fitted7CovM(1, 6);
292
293 fitted7CovM(2, 2) = fittedCov.dpzpz(); fitted7CovM(2, 2) = fitted7CovM(2, 2);
294 fitted7CovM(2, 3) = fitted7CovPart(2, 3); fitted7CovM(3, 2) = fitted7CovM(2, 3);
295 fitted7CovM(2, 4) = fittedCov.dxpz(); fitted7CovM(4, 2) = fitted7CovM(2, 4);
296 fitted7CovM(2, 5) = fittedCov.dypz(); fitted7CovM(5, 2) = fitted7CovM(2, 5);
297 fitted7CovM(2, 6) = fittedCov.dzpz(); fitted7CovM(6, 2) = fitted7CovM(2, 6);
298
299 fitted7CovM(3, 3) = fitted7CovPart(3, 3); fitted7CovM(3, 3) = fitted7CovM(3, 3);
300 fitted7CovM(3, 4) = fitted7CovPart(3, 4); fitted7CovM(4, 3) = fitted7CovM(3, 4);
301 fitted7CovM(3, 5) = fitted7CovPart(3, 5); fitted7CovM(5, 3) = fitted7CovM(3, 5);
302 fitted7CovM(3, 6) = fitted7CovPart(3, 6); fitted7CovM(6, 3) = fitted7CovM(3, 6);
303
304 fitted7CovM(4, 4) = fittedCov.dxx(); fitted7CovM(4, 4) = fitted7CovM(4, 4);
305 fitted7CovM(4, 5) = fittedCov.dxy(); fitted7CovM(5, 4) = fitted7CovM(4, 5);
306 fitted7CovM(4, 6) = fittedCov.dxz(); fitted7CovM(6, 4) = fitted7CovM(4, 6);
307
308 fitted7CovM(5, 5) = fittedCov.dyy(); fitted7CovM(5, 5) = fitted7CovM(5, 5);
309 fitted7CovM(5, 6) = fittedCov.dyz(); fitted7CovM(6, 5) = fitted7CovM(5, 6);
310
311 fitted7CovM(6, 6) = fittedCov.dzz(); fitted7CovM(6, 6) = fitted7CovM(6, 6);
312
313
314 float pValDau = m_belleDaughters[i]->getPValue();
315
316 m_belleDaughters[i]->updateMomentum(p4, x3, fitted7CovM, pValDau);
317
318 }
319
320}
DataType Z() const
access variable Z (= .at(2) without boundary check)
Definition: B2Vector3.h:435
DataType X() const
access variable X (= .at(0) without boundary check)
Definition: B2Vector3.h:431
DataType Y() const
access variable Y (= .at(1) without boundary check)
Definition: B2Vector3.h:433
Class to store reconstructed particles.
Definition: Particle.h:75
ROOT::Math::XYZVector getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:631
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:637
double getCharge(void) const
Returns particle charge.
Definition: Particle.cc:622
ROOT::Math::XYZVector getMomentum() const
Returns momentum vector.
Definition: Particle.h:560
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:420
Values of the result of a track fit with a given particle hypothesis.
short getChargeSign() const
Return track charge (1 or -1).
TMatrixDSym getCovariance6() const
Position and Momentum Covariance Matrix.
ROOT::Math::XYZVector getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
ROOT::Math::XYZVector getPosition() const
Getter for vector of position at closest approach of track in r/phi projection.
TMatrixDSym m_beamSpotCov
beam spot position covariance matrix.
Definition: RaveSetup.h:74
bool m_useBeamSpot
flag determines if beam spot information should be used for vertex fit.
Definition: RaveSetup.h:72
static RaveSetup * getRawInstance()
Same as getInstance(), but no check if the instance is initialised.
Definition: RaveSetup.cc:27
B2Vector3D m_beamSpot
beam spot position.
Definition: RaveSetup.h:73
rave::VertexFactory * m_raveVertexFactory
The RAVE vertex factory is the principal interface offered by the RAVE vertex fitting library.
Definition: RaveSetup.h:77
std::vector< Particle * > m_belleDaughters
Belle Particle pointers input.
TMatrixDSym getCov(VecSize vertexId=0) const
get the covariance matrix (3x3) of the of the fitted vertex position.
bool m_useBeamSpot
flag determines if the beam spot will be used or not.
int fit(std::string options="default")
do the vertex fit with all tracks previously added with the addTrack or addMother function.
std::vector< int >::size_type VecSize
I am using std::vector<int>::size_type because it is the official return value of ....
void addTrack(const Particle *const aParticlePtr)
add a track (in the format of a Belle2::Particle) to set of tracks that should be fitted to a vertex
std::vector< int > getTrackIdsForOneVertex(VecSize vertexId=0) const
get the ids of the tracks Rave associated with a specific vertex.
double getNdf(VecSize vertexId=0) const
get the number of degrees of freedom (NDF) of the fitted vertex.
rave::Track TrackFitResultToRaveTrack(const TrackFitResult *const aTrackPtr) const
converts a track from Belle2::TrackFitResult format to rave::Track format
double getChi2(VecSize vertexId=0) const
get the χ² of the fitted vertex.
B2Vector3D getPos(VecSize vertexId=0) const
get the position of the fitted vertex.
std::vector< rave::Vertex > m_raveVertices
holds the fitted vertices after fit() was called in the format used by Rave
void isVertexIdValid(const VecSize vertexId) const
checks if the vertex the user requested via one of the getters it actually there
void initBeamSpotMember()
Initialize m_useBeamSpot.
void updateDaughters()
update the Daughters particles
double getWeight(int trackId, VecSize vertexId=0) const
get the weight Rave assigned to a specific input track.
double getPValue(VecSize vertexId=0) const
get the p value of the fitted vertex.
std::vector< rave::Track > m_raveTracks
holds the tracks that were added to a RaveVertexFitter object in the format used by Rave
void addMother(const Particle *const aMotherParticlePtr)
All daughters of the argument of this function will be used as input for the vertex fit.
B2Vector3< double > B2Vector3D
typedef for common usage with double
Definition: B2Vector3.h:516
Abstract base class for different kinds of events.
Definition: ClusterUtils.h:24