Belle II Software development
TrackBuilder.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#include <tracking/trackFitting/trackBuilder/factories/TrackBuilder.h>
9
10
11#include <framework/datastore/StoreArray.h>
12#include <mdst/dataobjects/MCParticle.h>
13#include <mdst/dataobjects/HitPatternVXD.h>
14#include <mdst/dataobjects/HitPatternCDC.h>
15#include <mdst/dataobjects/Track.h>
16#include <mdst/dataobjects/TrackFitResult.h>
17#include <cdc/dataobjects/CDCRecoHit.h>
18#include <pxd/reconstruction/PXDRecoHit.h>
19#include <svd/reconstruction/SVDRecoHit.h>
20#include <svd/reconstruction/SVDRecoHit2D.h>
21#include <tracking/dataobjects/RecoTrack.h>
22
23#include <TMatrixDSym.h>
24
25#include <genfit/FitStatus.h>
26#include <genfit/KalmanFitterInfo.h>
27#include <genfit/Track.h>
28#include <genfit/TrackPoint.h>
29#include <genfit/MeasuredStateOnPlane.h>
30#include <genfit/FieldManager.h>
31
32using namespace Belle2;
33
35 const bool useClosestHitToIP, const bool useBFieldAtHit)
36{
39
40 const auto& trackReps = recoTrack.getRepresentations();
41 B2DEBUG(27, trackReps.size() << " track representations available.");
42 Track newTrack(recoTrack.getQualityIndicator());
43
44 bool repAlreadySet = false;
45 unsigned int repIDPlusOne = 0;
46 for (const auto& trackRep : trackReps) {
47 repIDPlusOne++;
48
49 // Check if the fitted particle type is in our charged stable set.
50 const Const::ParticleType particleType(std::abs(trackRep->getPDG()));
51 if (not Const::chargedStableSet.contains(particleType)) {
52 B2DEBUG(27, "Track fitted with hypothesis that is not a ChargedStable (PDG code = " << particleType.getPDGCode() << ")");
53 continue;
54 }
55
56 // Check if the fit worked.
57 if (not recoTrack.wasFitSuccessful(trackRep)) {
58 B2DEBUG(27, "The fit with the given track representation (" << std::abs(trackRep->getPDG()) <<
59 ") was not successful. Skipping ...");
60 continue;
61 }
62
63 if (not repAlreadySet) {
64 RecoTrackGenfitAccess::getGenfitTrack(recoTrack).setCardinalRep(repIDPlusOne - 1);
65 repAlreadySet = true;
66 }
67
68 // Extrapolate the tracks to the perigee.
69 genfit::MeasuredStateOnPlane msop;
70 try {
71 if (useClosestHitToIP) {
72 msop = recoTrack.getMeasuredStateOnPlaneClosestTo(ROOT::Math::XYZVector(0, 0, 0), trackRep);
73 } else {
74 msop = recoTrack.getMeasuredStateOnPlaneFromFirstHit(trackRep);
75 }
76 } catch (genfit::Exception& exception) {
77 B2WARNING(exception.what());
78 continue;
79 } catch (const std::runtime_error& er) {
80 B2WARNING("Runtime error encountered: " << er.what());
81 continue;
82 } catch (...) {
83 B2WARNING("Undefined exception encountered.");
84 continue;
85 }
86
87 genfit::MeasuredStateOnPlane extrapolatedMSoP = msop;
88 try {
89 extrapolatedMSoP.extrapolateToLine(m_beamSpot, m_beamAxis);
90 } catch (...) {
91 B2DEBUG(29, "Could not extrapolate the fit result for pdg " << particleType.getPDGCode() <<
92 " to the perigee point.");
93 continue;
94 }
95
96 // Build track fit result.
97
98 TVector3 poca(0., 0., 0.);
99 TVector3 dirInPoca(0., 0., 0.);
100 TMatrixDSym cov(6);
101 extrapolatedMSoP.getPosMomCov(poca, dirInPoca, cov);
102 B2DEBUG(29, "Point of closest approach: " << poca.x() << " " << poca.y() << " " << poca.z());
103 B2DEBUG(29, "Track direction in POCA: " << dirInPoca.x() << " " << dirInPoca.y() << " " << dirInPoca.z());
104
105 const int charge = recoTrack.getTrackFitStatus(trackRep)->getCharge();
106 const double pValue = recoTrack.getTrackFitStatus(trackRep)->getPVal();
107 const double nDF = recoTrack.getTrackFitStatus(trackRep)->getNdf();
108
109 double Bx, By, Bz; // In cgs units
110 if (useBFieldAtHit) {
111 const B2Vector3D& hitPosition = msop.getPos();
112 genfit::FieldManager::getInstance()->getFieldVal(hitPosition.X(), hitPosition.Y(), hitPosition.Z(), Bx, By, Bz);
113 } else {
114 genfit::FieldManager::getInstance()->getFieldVal(poca.X(), poca.Y(), poca.Z(), Bx, By, Bz);
115 }
116 Bz = Bz / 10.; // In SI-Units
117
118 const uint64_t hitPatternCDCInitializer = getHitPatternCDCInitializer(recoTrack, trackRep);
119 const uint32_t hitPatternVXDInitializer = getHitPatternVXDInitializer(recoTrack, trackRep);
120
121 const auto newTrackFitResult = trackFitResults.appendNew(
122 ROOT::Math::XYZVector(poca), ROOT::Math::XYZVector(dirInPoca), cov, charge, particleType, pValue, Bz,
123 hitPatternCDCInitializer, hitPatternVXDInitializer, nDF
124 );
125
126 const int newTrackFitResultArrayIndex = newTrackFitResult->getArrayIndex();
127 newTrack.setTrackFitResultIndex(particleType, newTrackFitResultArrayIndex);
128 }
129
130 B2DEBUG(27, "Number of fitted hypothesis = " << newTrack.getNumberOfFittedHypotheses());
131 if (newTrack.getNumberOfFittedHypotheses() > 0) {
132 Track* addedTrack = tracks.appendNew(newTrack);
133 addedTrack->addRelationTo(&recoTrack);
134 return true;
135 } else {
136 B2DEBUG(28, "No valid fit for any given hypothesis. No Track is added to the Tracks StoreArray.");
137 }
138 return true;
139}
140
141
142uint32_t TrackBuilder::getHitPatternVXDInitializer(const RecoTrack& recoTrack, const genfit::AbsTrackRep* representation)
143{
144 HitPatternVXD hitPatternVXD;
145
146 const auto& hitPointsWithMeasurements = recoTrack.getHitPointsWithMeasurement();
147 int nNotFittedVXDhits = 0;
148
149 for (const auto& trackPoint : hitPointsWithMeasurements) { // Loop on TrackPoint
150
151 genfit::KalmanFitterInfo* kalmanInfo = trackPoint->getKalmanFitterInfo(representation);
152
153 for (size_t measurementId = 0; measurementId < trackPoint->getNumRawMeasurements(); measurementId++) { //Loop on raw measurement
154
155 genfit::AbsMeasurement* absMeas = trackPoint->getRawMeasurement(measurementId);
156
157 PXDRecoHit* pxdHit = dynamic_cast<PXDRecoHit*>(absMeas);
158 SVDRecoHit* svdHit = dynamic_cast<SVDRecoHit*>(absMeas);
159 SVDRecoHit2D* svdHit2D = dynamic_cast<SVDRecoHit2D*>(absMeas);
160
161 if (!pxdHit && !svdHit2D && !svdHit)
162 continue; // consider only VXD hits
163
164 if (kalmanInfo) {
165
166 if (kalmanInfo->getNumMeasurements() > 1)
167 B2WARNING("VXD TrackPoint contains more than one KalmanFitterInfo: only the first will be considered");
168
169 const double weight = kalmanInfo->getWeights().at(0); // only 1st KalmanFitterInfo considered
170 if (weight < 1.e-9)
171 continue; // skip kfinfo with negligible weight
172
173 if (pxdHit) {
174 const int layerNumber = pxdHit->getSensorID().getLayerNumber();
175 const int currentHits = hitPatternVXD.getPXDLayer(layerNumber, HitPatternVXD::PXDMode::normal);
176 hitPatternVXD.setPXDLayer(layerNumber, currentHits + 1, HitPatternVXD::PXDMode::normal);
177 }
178
179 if (svdHit2D) {
180 const int layerNumber = svdHit2D->getSensorID().getLayerNumber();
181 const auto& currentHits = hitPatternVXD.getSVDLayer(layerNumber);
182 hitPatternVXD.setSVDLayer(layerNumber, currentHits.first + 1, currentHits.second + 1);
183 } else if (svdHit) {
184 const int layerNumber = svdHit->getSensorID().getLayerNumber();
185 const auto& currentHits = hitPatternVXD.getSVDLayer(layerNumber);
186
187 if (svdHit->isU())
188 hitPatternVXD.setSVDLayer(layerNumber, currentHits.first + 1, currentHits.second);
189 else
190 hitPatternVXD.setSVDLayer(layerNumber, currentHits.first, currentHits.second + 1);
191 }
192
193 } // end of if kalmanInfo
194 else {
195 // i.e. if !kalmanInfo)
196 ++nNotFittedVXDhits; // counts TrackPoints with VXD hits without KalmanFitterInfo
197 continue;
198 }
199
200 } // end of loop on raw measurements
201 } // end of loop on TrackPoint
202
203 if (nNotFittedVXDhits > 0) {
204 B2DEBUG(27, " No KalmanFitterInfo associated to some TrackPoints with VXD hits, not filling the HitPatternVXD");
205 B2DEBUG(27, nNotFittedVXDhits << " had no FitterInfo");
206 }
207 return hitPatternVXD.getInteger();
208}
209
210
211uint64_t TrackBuilder::getHitPatternCDCInitializer(const RecoTrack& recoTrack, const genfit::AbsTrackRep* representation)
212{
213 HitPatternCDC hitPatternCDC;
214
215 int nCDChits = 0;
216 int nNotFittedCDChits = 0;
217
218 const auto& hitPointsWithMeasurements = recoTrack.getHitPointsWithMeasurement();
219
220 for (const auto& trackPoint : hitPointsWithMeasurements) { // Loop on TrackPoint
221
222 genfit::KalmanFitterInfo* kalmanInfo = trackPoint->getKalmanFitterInfo(representation);
223
224 for (size_t measurementId = 0; measurementId < trackPoint->getNumRawMeasurements(); measurementId++) { //Loop on raw measurement
225
226 genfit::AbsMeasurement* absMeas = trackPoint->getRawMeasurement(measurementId);
227 CDCRecoHit* cdcHit = dynamic_cast<CDCRecoHit*>(absMeas);
228
229 if (!cdcHit)
230 continue; // consider only CDC hits
231
232 if (kalmanInfo) {
233 bool isValidWeight = false;
234 for (unsigned int kfinfoId = 0; kfinfoId < kalmanInfo->getNumMeasurements(); kfinfoId++) { // Loop on KalmanFitterInfo
235 const double weight = kalmanInfo->getWeights().at(kfinfoId);
236 if (weight < 1.e-9)
237 continue; // skip kfinfo with negligible weight
238 else {
239 isValidWeight = true;
240 B2DEBUG(27, "CDC: " << nCDChits << "\t" << cdcHit->getWireID().getEWire() << "\t" << kfinfoId << "\t" << weight);
241 break;
242 }
243 } // end of KalmanFitterInfo loop
244
245 if (isValidWeight) { // fill nCDChits only one time per each raw measurement
246 WireID wire = cdcHit->getWireID();
247 hitPatternCDC.setLayer(wire.getICLayer());
248 nCDChits++; // counts CDC hits where there is KalmanFitterInfo and not negligible weight
249 }
250
251 } // end of if kalmanInfo
252 else {
253 // i.e. if !kalmanInfo)
254 ++nNotFittedCDChits; // counts TrackPoints with CDC hits without KalmanFitterInfo
255 continue;
256 }
257
258 } // end of loop on raw measurements
259 } // end of loop on TrackPoint
260
261 if (nNotFittedCDChits > 0) {
262 B2DEBUG(27, " No KalmanFitterInfo associated to some TrackPoints with CDC hits, not filling the HitPatternCDC");
263 B2DEBUG(27, nNotFittedCDChits << " out of " << nCDChits << " had no FitterInfo");
264 }
265 hitPatternCDC.setNHits(nCDChits);
266
267 return hitPatternCDC.getInteger();
268
269}
270
271
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
This class is used to transfer CDC information to the track fit.
Definition: CDCRecoHit.h:32
WireID getWireID() const
Getter for WireID object.
Definition: CDCRecoHit.h:49
The ParticleType class for identifying different particle types.
Definition: Const.h:408
int getPDGCode() const
PDG code.
Definition: Const.h:473
static const ParticleSet chargedStableSet
set of charged stable particles
Definition: Const.h:618
Hit pattern of CDC hits within a track.
Definition: HitPatternCDC.h:35
void setNHits(unsigned short nHits)
Sets the 8 MSBs to the total number of hits in the CDC.
void setLayer(const unsigned short layer)
Set bit corresponding to layer to true.
ULong64_t getInteger() const
Getter for underlying integer type.
Hit pattern of the VXD within a track.
Definition: HitPatternVXD.h:37
unsigned int getInteger() const
Getter for the underlying integer.
Definition: HitPatternVXD.h:58
void setSVDLayer(const unsigned short layerId, unsigned short uHits, unsigned short vHits)
Set the number of hits in a specific layer of the SVD.
void setPXDLayer(const unsigned short layerId, unsigned short nHits, const PXDMode &mode=PXDMode::normal)
Set the number of hits in a specific layer of the PXD.
unsigned short getPXDLayer(const unsigned short layerId, const PXDMode &mode=PXDMode::normal) const
Get the number of hits in a specific layer of the PXD.
std::pair< const unsigned short, const unsigned short > getSVDLayer(const unsigned short layerId) const
Get the number of hits in a specific layer of the SVD.
PXDRecoHit - an extended form of PXDCluster containing geometry information.
Definition: PXDRecoHit.h:49
VxdID getSensorID() const
Get the compact ID.
Definition: PXDRecoHit.h:101
static genfit::Track & getGenfitTrack(RecoTrack &recoTrack)
Give access to the RecoTrack's genfit::Track.
Definition: RecoTrack.cc:404
This is the Reconstruction Event-Data Model Track.
Definition: RecoTrack.h:79
const std::vector< genfit::AbsTrackRep * > & getRepresentations() const
Return a list of track representations. You are not allowed to modify or delete them!
Definition: RecoTrack.h:638
bool wasFitSuccessful(const genfit::AbsTrackRep *representation=nullptr) const
Returns true if the last fit with the given representation was successful.
Definition: RecoTrack.cc:336
const genfit::MeasuredStateOnPlane & getMeasuredStateOnPlaneClosestTo(const ROOT::Math::XYZVector &closestPoint, const genfit::AbsTrackRep *representation=nullptr)
Return genfit's MasuredStateOnPlane, that is closest to the given point useful for extrapolation of m...
Definition: RecoTrack.cc:426
float getQualityIndicator() const
Get the quality index attached to this RecoTrack given by one of the reconstruction algorithms....
Definition: RecoTrack.h:841
const genfit::MeasuredStateOnPlane & getMeasuredStateOnPlaneFromFirstHit(const genfit::AbsTrackRep *representation=nullptr) const
Return genfit's MeasuredStateOnPlane for the first hit in a fit useful for extrapolation of measureme...
Definition: RecoTrack.cc:605
const genfit::FitStatus * getTrackFitStatus(const genfit::AbsTrackRep *representation=nullptr) const
Return the track fit status for the given representation or for the cardinal one. You are not allowed...
Definition: RecoTrack.h:621
const std::vector< genfit::TrackPoint * > & getHitPointsWithMeasurement() const
Return a list of measurements and track points, which can be used e.g. to extrapolate....
Definition: RecoTrack.h:708
void addRelationTo(const RelationsInterface< BASE > *object, float weight=1.0, const std::string &namedRelation="") const
Add a relation from this object to another object (with caching).
SVDRecoHit - an extended form of SVDHit containing geometry information.
Definition: SVDRecoHit2D.h:46
VxdID getSensorID() const
Get the compact ID.
Definition: SVDRecoHit2D.h:100
SVDRecoHit - an extended form of SVDHit containing geometry information.
Definition: SVDRecoHit.h:47
bool isU() const
Is the coordinate u or v?
Definition: SVDRecoHit.h:91
VxdID getSensorID() const
Get the compact ID.
Definition: SVDRecoHit.h:82
Accessor to arrays stored in the data store.
Definition: StoreArray.h:113
T * appendNew()
Construct a new T object at the end of the array.
Definition: StoreArray.h:246
static uint32_t getHitPatternVXDInitializer(const RecoTrack &recoTrack, const genfit::AbsTrackRep *representation=nullptr)
Get the HitPattern in the VXD.
bool storeTrackFromRecoTrack(RecoTrack &recoTrack, const bool useClosestHitToIP=false, const bool useBFieldAtHit=false)
Stores a Belle2 Track from a Reco Track.
Definition: TrackBuilder.cc:34
B2Vector3D m_beamSpot
Extrapolation target, origin.
Definition: TrackBuilder.h:91
std::string m_trackColName
TrackColName (output).
Definition: TrackBuilder.h:87
static uint64_t getHitPatternCDCInitializer(const RecoTrack &recoTrack, const genfit::AbsTrackRep *representation=nullptr)
Get the HitPattern in the CDC.
std::string m_trackFitResultColName
TrackFitResultColName (output).
Definition: TrackBuilder.h:89
B2Vector3D m_beamAxis
Extrapolation target, positive z direction.
Definition: TrackBuilder.h:93
Class that bundles various TrackFitResults.
Definition: Track.h:25
void setTrackFitResultIndex(const Const::ChargedStable &chargedStable, short index)
Set an index (for positive values) or unavailability-code (index = -1) for a specific mass hypothesis...
Definition: Track.h:174
unsigned int getNumberOfFittedHypotheses() const
Returns the number of fitted hypothesis which are stored in this track.
Definition: Track.cc:36
baseType getLayerNumber() const
Get the layer id.
Definition: VxdID.h:96
Class to identify a wire inside the CDC.
Definition: WireID.h:34
unsigned short getICLayer() const
Getter for continuous layer numbering.
Definition: WireID.cc:24
unsigned short getEWire() const
Getter for encoded wire number.
Definition: WireID.h:154
Abstract base class for different kinds of events.