Belle II Software  release-05-02-19
RaveVertexFitter.cc
1 /**************************************************************************
2  * BASF2 (Belle Analysis Framework 2) *
3  * Copyright(C) 2012 - Belle II Collaboration *
4  * *
5  * Author: The Belle II Collaboration *
6  * Contributors: Moritz Nadler, Luigi Li Gioi *
7  * *
8  * This software is provided "as is" without any warranty. *
9  **************************************************************************/
10 
11 
12 #include <analysis/VertexFitting/RaveInterface/RaveVertexFitter.h>
13 #include <analysis/VertexFitting/RaveInterface/RaveSetup.h>
14 #include <mdst/dataobjects/Track.h>
15 
16 #include <rave/VertexFactory.h>
17 
18 //root
19 #include <Math/ProbFunc.h>
20 //stl
21 using std::string;
22 #include <vector>
23 using std::vector;
24 
25 using namespace Belle2;
26 using namespace analysis;
27 
28 
29 
30 RaveVertexFitter::RaveVertexFitter(): m_useBeamSpot(false)
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 
55 void RaveVertexFitter::addTrack(const TrackFitResult* const aTrackPtr)
56 {
57  m_raveTracks.push_back(TrackFitResultToRaveTrack(aTrackPtr));
58 }
59 
60 rave::Track RaveVertexFitter::TrackFitResultToRaveTrack(const TrackFitResult* const aTrackPtr) const
61 {
62  const int id = m_raveTracks.size();
63 
64  TVector3 pos = aTrackPtr->getPosition();
65  TVector3 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 
87 void RaveVertexFitter::addTrack(const Particle* aParticlePtr)
88 {
89  const int id = m_raveTracks.size();
90  const TVector3& pos = aParticlePtr->getVertex();
91  const TVector3& 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 
110 void 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 
121 int 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) {
141  const TVector3& bsPos = RaveSetup::getRawInstance()->m_beamSpot;
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 
155 void RaveVertexFitter::isVertexIdValid(const VecSize vertexId) const
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 
169 TVector3 RaveVertexFitter::getPos(VecSize vertexId) const
170 {
171  isVertexIdValid(vertexId);
172 
173  return TVector3(m_raveVertices[vertexId].position().x(), m_raveVertices[vertexId].position().y(),
174  m_raveVertices[vertexId].position().z());
175 
176 }
177 
178 double 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( "returing weight for track with x coord: " <<weightedTracks[i].second.state().x() );
187  return weightedTracks[i].first;
188  }
189  }
190  return -1;
191 
192 }
193 
194 std::vector<int> RaveVertexFitter::getTrackIdsForOneVertex(VecSize vertexId) const
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 
206 double RaveVertexFitter::getPValue(VecSize vertexId) const
207 {
208  isVertexIdValid(vertexId);
209 
210  return ROOT::Math::chisquared_cdf_c(m_raveVertices[vertexId].chiSquared(), m_raveVertices[vertexId].ndf());
211 
212 }
213 
214 double RaveVertexFitter::getNdf(VecSize vertexId) const
215 {
216  isVertexIdValid(vertexId);
217 
218  return m_raveVertices[vertexId].ndf();
219 
220 }
221 
222 double RaveVertexFitter::getChi2(VecSize vertexId) const
223 {
224  isVertexIdValid(vertexId);
225 
226  return m_raveVertices[vertexId].chiSquared();
227 
228 }
229 
230 TMatrixDSym 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 sigle 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  TVector3 x3(fittedV.x(), fittedV.y(), fittedV.z());
271  TLorentzVector p4;
272  double fittedE = sqrt(fittedP.mag2() + (m_belleDaughters[i]->getMass() * (m_belleDaughters[i]->getMass())));
273  p4.SetXYZT(fittedP.x(), fittedP.y(), fittedP.z(), fittedE);
274 
275  TMatrixDSym fitted7CovPart = m_belleDaughters[i]->getMomentumVertexErrorMatrix() ;
276 
277  //px,py,pz,E,x,y,z
278  TMatrixDSym fitted7CovM(7);
279  fitted7CovM(0, 0) = fittedCov.dpxpx(); fitted7CovM(0, 0) = fitted7CovM(0, 0);
280  fitted7CovM(0, 1) = fittedCov.dpxpy(); fitted7CovM(1, 0) = fitted7CovM(0, 1);
281  fitted7CovM(0, 2) = fittedCov.dpxpz(); fitted7CovM(2, 0) = fitted7CovM(0, 2);
282  fitted7CovM(0, 3) = fitted7CovPart(0, 3); fitted7CovM(3, 0) = fitted7CovM(0, 3);
283  fitted7CovM(0, 4) = fittedCov.dxpx(); fitted7CovM(4, 0) = fitted7CovM(0, 4);
284  fitted7CovM(0, 5) = fittedCov.dypx(); fitted7CovM(5, 0) = fitted7CovM(0, 5);
285  fitted7CovM(0, 6) = fittedCov.dzpx(); fitted7CovM(6, 0) = fitted7CovM(0, 6);
286 
287  fitted7CovM(1, 1) = fittedCov.dpypy(); fitted7CovM(1, 1) = fitted7CovM(1, 1);
288  fitted7CovM(1, 2) = fittedCov.dpypz(); fitted7CovM(2, 1) = fitted7CovM(1, 2);
289  fitted7CovM(1, 3) = fitted7CovPart(1, 3); fitted7CovM(3, 1) = fitted7CovM(1, 3);
290  fitted7CovM(1, 4) = fittedCov.dxpy(); fitted7CovM(4, 1) = fitted7CovM(1, 4);
291  fitted7CovM(1, 5) = fittedCov.dypy(); fitted7CovM(5, 1) = fitted7CovM(1, 5);
292  fitted7CovM(1, 6) = fittedCov.dzpy(); fitted7CovM(6, 1) = fitted7CovM(1, 6);
293 
294  fitted7CovM(2, 2) = fittedCov.dpzpz(); fitted7CovM(2, 2) = fitted7CovM(2, 2);
295  fitted7CovM(2, 3) = fitted7CovPart(2, 3); fitted7CovM(3, 2) = fitted7CovM(2, 3);
296  fitted7CovM(2, 4) = fittedCov.dxpz(); fitted7CovM(4, 2) = fitted7CovM(2, 4);
297  fitted7CovM(2, 5) = fittedCov.dypz(); fitted7CovM(5, 2) = fitted7CovM(2, 5);
298  fitted7CovM(2, 6) = fittedCov.dzpz(); fitted7CovM(6, 2) = fitted7CovM(2, 6);
299 
300  fitted7CovM(3, 3) = fitted7CovPart(3, 3); fitted7CovM(3, 3) = fitted7CovM(3, 3);
301  fitted7CovM(3, 4) = fitted7CovPart(3, 4); fitted7CovM(4, 3) = fitted7CovM(3, 4);
302  fitted7CovM(3, 5) = fitted7CovPart(3, 5); fitted7CovM(5, 3) = fitted7CovM(3, 5);
303  fitted7CovM(3, 6) = fitted7CovPart(3, 6); fitted7CovM(6, 3) = fitted7CovM(3, 6);
304 
305  fitted7CovM(4, 4) = fittedCov.dxx(); fitted7CovM(4, 4) = fitted7CovM(4, 4);
306  fitted7CovM(4, 5) = fittedCov.dxy(); fitted7CovM(5, 4) = fitted7CovM(4, 5);
307  fitted7CovM(4, 6) = fittedCov.dxz(); fitted7CovM(6, 4) = fitted7CovM(4, 6);
308 
309  fitted7CovM(5, 5) = fittedCov.dyy(); fitted7CovM(5, 5) = fitted7CovM(5, 5);
310  fitted7CovM(5, 6) = fittedCov.dyz(); fitted7CovM(6, 5) = fitted7CovM(5, 6);
311 
312  fitted7CovM(6, 6) = fittedCov.dzz(); fitted7CovM(6, 6) = fitted7CovM(6, 6);
313 
314 
315  float pValDau = m_belleDaughters[i]->getPValue();
316 
317  m_belleDaughters[i]->updateMomentum(p4, x3, fitted7CovM, pValDau);
318 
319  }
320 
321 }
Belle2::analysis::RaveVertexFitter::isVertexIdValid
void isVertexIdValid(const VecSize vertexId) const
checks if the vertex the user requested via one of the getters it actually there
Definition: RaveVertexFitter.cc:155
Belle2::Vector3D
HepGeom::Vector3D< double > Vector3D
3D Vector
Definition: Cell.h:35
Belle2::analysis::RaveVertexFitter::getPValue
double getPValue(VecSize vertexId=0) const
get the p value of the fitted vertex.
Definition: RaveVertexFitter.cc:206
Belle2::Particle::getVertex
TVector3 getVertex() const
Returns vertex position (POCA for charged, IP for neutral FS particles)
Definition: Particle.h:529
Belle2::analysis::RaveVertexFitter::getPos
TVector3 getPos(VecSize vertexId=0) const
get the position of the fitted vertex.
Definition: RaveVertexFitter.cc:169
Belle2::TrackFitResult::getMomentum
TVector3 getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:116
Belle2::analysis::RaveVertexFitter::m_belleDaughters
std::vector< Particle * > m_belleDaughters
Belle Particle pointers input.
Definition: RaveVertexFitter.h:153
Belle2::analysis::RaveVertexFitter::getChi2
double getChi2(VecSize vertexId=0) const
get the χ² of the fitted vertex.
Definition: RaveVertexFitter.cc:222
Belle2::analysis::RaveSetup::m_useBeamSpot
bool m_useBeamSpot
flag determines if beam spot information should be used for vertex fit.
Definition: RaveSetup.h:74
Belle2::analysis::RaveVertexFitter::TrackFitResultToRaveTrack
rave::Track TrackFitResultToRaveTrack(const TrackFitResult *const aTrackPtr) const
converts a track from Belle2::TrackFitResult format to rave::Track format
Definition: RaveVertexFitter.cc:60
Belle2::TrackFitResult
Values of the result of a track fit with a given particle hypothesis.
Definition: TrackFitResult.h:59
Belle2::analysis::RaveVertexFitter::m_raveVertices
std::vector< rave::Vertex > m_raveVertices
holds the fitted vertices after fit() was called in the format used by Rave
Definition: RaveVertexFitter.h:150
Belle2::TrackFitResult::getPosition
TVector3 getPosition() const
Getter for vector of position at closest approach of track in r/phi projection.
Definition: TrackFitResult.h:109
Belle2::analysis::RaveVertexFitter::m_raveTracks
std::vector< rave::Track > m_raveTracks
holds the tracks that were added to a RaveVertexFitter object in the format used by Rave
Definition: RaveVertexFitter.h:147
Belle2::analysis::RaveVertexFitter::m_useBeamSpot
bool m_useBeamSpot
flag determines if the beam spot will be used or not.
Definition: RaveVertexFitter.h:141
Belle2::Particle::getCharge
float getCharge(void) const
Returns particle charge.
Definition: Particle.cc:586
Belle2::analysis::RaveSetup::m_raveVertexFactory
rave::VertexFactory * m_raveVertexFactory
The RAVE vertex factory is the principal interface offered by the RAVE vertex fitting library.
Definition: RaveSetup.h:79
Belle2::Particle::getDaughters
std::vector< Belle2::Particle * > getDaughters() const
Returns a vector of pointers to daughter particles.
Definition: Particle.cc:601
Belle2
Abstract base class for different kinds of events.
Definition: MillepedeAlgorithm.h:19
Belle2::analysis::RaveVertexFitter::~RaveVertexFitter
~RaveVertexFitter()
Destructor.
Belle2::Particle::getMomentumVertexErrorMatrix
TMatrixFSym getMomentumVertexErrorMatrix() const
Returns 7x7 error matrix.
Definition: Particle.cc:386
Belle2::analysis::RaveVertexFitter::getNdf
double getNdf(VecSize vertexId=0) const
get the number of degrees of freedom (NDF) of the fitted vertex.
Definition: RaveVertexFitter.cc:214
Belle2::analysis::RaveVertexFitter::addTrack
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
Definition: RaveVertexFitter.cc:87
Belle2::analysis::RaveVertexFitter::RaveVertexFitter
RaveVertexFitter()
The constructor.
Definition: RaveVertexFitter.cc:30
Belle2::analysis::RaveVertexFitter::getCov
TMatrixDSym getCov(VecSize vertexId=0) const
get the covariance matrix (3x3) of the of the fitted vertex position.
Definition: RaveVertexFitter.cc:230
Belle2::analysis::RaveVertexFitter::initBeamSpotMember
void initBeamSpotMember()
Initialize m_useBeamSpot.
Definition: RaveVertexFitter.cc:41
Belle2::analysis::RaveVertexFitter::getTrackIdsForOneVertex
std::vector< int > getTrackIdsForOneVertex(VecSize vertexId=0) const
get the ids of the tracks Rave associated with a specific vertex.
Definition: RaveVertexFitter.cc:194
Belle2::Particle::getMomentum
TVector3 getMomentum() const
Returns momentum vector.
Definition: Particle.h:475
Belle2::analysis::RaveVertexFitter::updateDaughters
void updateDaughters()
update the Daughters particles
Definition: RaveVertexFitter.cc:250
Belle2::Particle
Class to store reconstructed particles.
Definition: Particle.h:77
Belle2::analysis::RaveVertexFitter::addMother
void addMother(const Particle *const aMotherParticlePtr)
All daughters of the argument of this function will be used as input for the vertex fit.
Definition: RaveVertexFitter.cc:110
Belle2::analysis::RaveSetup::getRawInstance
static RaveSetup * getRawInstance()
Same as getInstance(), but no check if the instance is initialised.
Definition: RaveSetup.cc:29
Belle2::analysis::RaveVertexFitter::VecSize
std::vector< int >::size_type VecSize
I am using std::vector<int>::size_type because it is the official return value of ....
Definition: RaveVertexFitter.h:48
Belle2::analysis::RaveVertexFitter::fit
int fit(std::string options="default")
do the vertex fit with all tracks previously added with the addTrack or addMother function.
Definition: RaveVertexFitter.cc:121
Belle2::analysis::RaveSetup::m_beamSpot
TVector3 m_beamSpot
beam spot position.
Definition: RaveSetup.h:75
Belle2::TrackFitResult::getCovariance6
TMatrixDSym getCovariance6() const
Position and Momentum Covariance Matrix.
Definition: TrackFitResult.h:147
Belle2::Point3D
HepGeom::Point3D< double > Point3D
3D point
Definition: Cell.h:33
Belle2::analysis::RaveVertexFitter::getWeight
double getWeight(int trackId, VecSize vertexId=0) const
get the weight Rave assigned to a specific input track.
Definition: RaveVertexFitter.cc:178
Belle2::TrackFitResult::getChargeSign
short getChargeSign() const
Return track charge (1 or -1).
Definition: TrackFitResult.h:160
Belle2::analysis::RaveSetup::m_beamSpotCov
TMatrixDSym m_beamSpotCov
beam spot position covariance matrix.
Definition: RaveSetup.h:76