Belle II Software development
VXDDedxPIDModule.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#include <reconstruction/modules/VXDDedxPID/VXDDedxPIDModule.h>
10
11#include <framework/gearbox/Const.h>
12
13#include <vxd/geometry/GeoCache.h>
14
15#include <genfit/MaterialEffects.h>
16
17#include <Math/VectorUtil.h>
18#include <TH2F.h>
19
20#include <memory>
21#include <cassert>
22#include <cmath>
23#include <algorithm>
24
25using namespace Belle2;
26using namespace Dedx;
27
28REG_MODULE(VXDDedxPID);
29
31{
33
34 //Set module properties
35 setDescription("Extract dE/dx and corresponding log-likelihood from fitted tracks and hits in the SVD and PXD.");
36
37 //Parameter definitions
38 addParam("useIndividualHits", m_useIndividualHits,
39 "Use individual hits (true) or truncated mean (false) to determine likelihoods", false);
40 addParam("onlyPrimaryParticles", m_onlyPrimaryParticles,
41 "For MC only: if true, only save data for primary particles (as determined by MC truth)", false);
42 addParam("usePXD", m_usePXD, "Use PXDClusters for dE/dx calculation", false);
43 addParam("useSVD", m_useSVD, "Use SVDClusters for dE/dx calculation", true);
44
45 m_eventID = -1;
46 m_trackID = 0;
47}
48
50
52{
53 if (m_usePXD) {
54 if (not m_PXDDedxPDFs) B2FATAL("No PXD dE/dx PDF's available");
55 bool ok = m_PXDDedxPDFs->checkPDFs(not m_useIndividualHits);
56 if (not ok) B2FATAL("Binning or ranges of PXD dE/dx PDF's differ");
57 }
58 if (m_useSVD) {
59 if (not m_SVDDedxPDFs) B2FATAL("No SVD Dedx PDF's available");
60 bool ok = m_SVDDedxPDFs->checkPDFs(not m_useIndividualHits);
61 if (not ok) B2FATAL("Binning or ranges of SVD dE/dx PDF's differ");
62 }
63}
64
66{
67
68 // required inputs
69 m_tracks.isRequired();
70 m_recoTracks.isRequired();
71
72 //optional inputs
73 m_mcparticles.isOptional();
74 m_tracks.optionalRelationTo(m_mcparticles);
75
76 if (m_useSVD)
77 m_svdClusters.isRequired();
78 else
79 m_svdClusters.isOptional();
80 if (m_usePXD)
81 m_pxdClusters.isRequired();
82 else
83 m_pxdClusters.isOptional();
84
85 // register dE/dx data points
86 m_dedxTracks.registerInDataStore();
87 m_tracks.registerRelationTo(m_dedxTracks);
88
89 // register likelihoods
90 m_dedxLikelihoods.registerInDataStore();
91 m_tracks.registerRelationTo(m_dedxLikelihoods);
92
93 m_SVDDedxPDFs.addCallback([this]() {checkPDFs();});
94 m_PXDDedxPDFs.addCallback([this]() {checkPDFs();});
95 checkPDFs();
96
97
98 // create instances here to not confuse profiling
100
101 if (!genfit::MaterialEffects::getInstance()->isInitialized()) {
102 B2FATAL("Need to have SetupGenfitExtrapolationModule in path before this one.");
103 }
104}
105
107{
108 m_dedxTracks.clear();
109 m_dedxLikelihoods.clear();
110
111 // go through Tracks
112 // get fitresult and gftrack and do extrapolations, save corresponding dE/dx and likelihood values
113 // get genfit::TrackCand through genfit::Track::getCand()
114 // get hit indices through genfit::TrackCand::getHit(...)
115 // create one VXDDedxTrack per fitresult/gftrack
116 //create one DedkLikelihood per Track (plus rel)
117 m_eventID++;
118
119 // inputs
120 const int numMCParticles = m_mcparticles.getEntries();
121
122 // **************************************************
123 //
124 // LOOP OVER TRACKS
125 //
126 // **************************************************
127
128 for (const auto& track : m_tracks) {
129 m_trackID++;
130
131 std::shared_ptr<VXDDedxTrack> dedxTrack = std::make_shared<VXDDedxTrack>();
132 dedxTrack->m_eventID = m_eventID;
133 dedxTrack->m_trackID = m_trackID;
134
135 // load the pion fit hypothesis or the hypothesis which is the closest in mass to a pion
136 // the tracking will not always successfully fit with a pion hypothesis
137 const TrackFitResult* fitResult = track.getTrackFitResultWithClosestMass(Const::pion);
138 if (!fitResult) {
139 B2WARNING("No related fit for track ...");
140 continue;
141 }
142
143 if (numMCParticles != 0) {
144 // find MCParticle corresponding to this track
145 const MCParticle* mcpart = track.getRelatedTo<MCParticle>();
146
147 if (mcpart) {
149 continue; //not a primary particle, ignore
150 }
151
152 //add some MC truths to VXDDedxTrack object
153 dedxTrack->m_pdg = mcpart->getPDG();
154 const MCParticle* mother = mcpart->getMother();
155 dedxTrack->m_motherPDG = mother ? mother->getPDG() : 0;
156
157 const ROOT::Math::XYZVector trueMomentum = mcpart->getMomentum();
158 dedxTrack->m_pTrue = trueMomentum.R();
159 }
160 }
161
162 // get momentum (at origin) from fit result
163 const ROOT::Math::XYZVector& trackMom = fitResult->getMomentum();
164 dedxTrack->m_p = trackMom.R();
165 dedxTrack->m_cosTheta = cos(trackMom.Theta());
166 dedxTrack->m_charge = fitResult->getChargeSign();
167
168 // dE/dx values will be calculated using associated RecoTrack
169 const RecoTrack* recoTrack = track.getRelatedTo<RecoTrack>();
170 if (!recoTrack) {
171 B2WARNING("No related track for this fit...");
172 continue;
173 }
174
175 // Check to see if the track is pruned. We use the cardinal representation for this, as we do not expect that one
176 // representation is pruned whereas the other is not.
177 if (recoTrack->hasTrackFitStatus()) {
178 if (recoTrack->getTrackFitStatus()->isTrackPruned()) {
179 B2ERROR("GFTrack is pruned, please run VXDDedxPID only on unpruned tracks! Skipping this track.");
180 continue;
181 }
182 } else // The RecoTrack might not have a fit status for reasons: let's skip it
183 continue;
184
185 if (m_usePXD) {
186 const std::vector<PXDCluster*>& pxdClusters = recoTrack->getPXDHitList();
187 saveSiHits(dedxTrack.get(), pxdClusters, recoTrack);
188 }
189
190 if (m_useSVD) {
191 const std::vector<SVDCluster*>& svdClusters = recoTrack->getSVDHitList();
192 saveSiHits(dedxTrack.get(), svdClusters, recoTrack);
193 }
194
195 if (dedxTrack->dedx.empty()) {
196 B2DEBUG(50, "Found track with no hits, ignoring.");
197 continue;
198 }
199
200 // add a few last things to the VXDDedxTrack
201 const int numDedx = dedxTrack->dedx.size();
202 dedxTrack->m_nHits = numDedx;
203 // no need to define lowedgetruncated and highedgetruncated as we always remove the highest 2 dE/dx values from 8 dE/dx value
204 dedxTrack->m_nHitsUsed = numDedx - 2;
205
206 // calculate log likelihoods
207 dedxTrack->clearLogLikelihoods();
208 bool truncated = not m_useIndividualHits;
209 if (m_usePXD) dedxTrack->addLogLikelihoods(m_PXDDedxPDFs->getPDFs(truncated), Dedx::c_PXD, truncated);
210 if (m_useSVD) dedxTrack->addLogLikelihoods(m_SVDDedxPDFs->getPDFs(truncated), Dedx::c_SVD, truncated);
211
212 // save log likelihoods
213 if (dedxTrack->areLogLikelihoodsAvailable()) {
214 VXDDedxLikelihood* likelihoodObj = m_dedxLikelihoods.appendNew(dedxTrack->m_vxdLogl);
215 track.addRelationTo(likelihoodObj);
216 }
217
218 // save the VXDDedxTrack
219 VXDDedxTrack* newVXDDedxTrack = m_dedxTracks.appendNew(*dedxTrack);
220 track.addRelationTo(newVXDDedxTrack);
221
222 } // end of loop over tracks
223}
224
226{
227
228 B2DEBUG(50, "VXDDedxPIDModule exiting after processing " << m_trackID <<
229 " tracks in " << m_eventID + 1 << " events.");
230}
231
232
233// calculateMeans need some change as we always remove highest 2 dE/dx values
234void VXDDedxPIDModule::calculateMeans(double& mean, double& truncatedMean, double& truncatedMeanErr,
235 const std::vector<double>& dedx) const
236{
237 // Calculate the truncated average by skipping only highest two value
238 std::vector<double> sortedDedx = dedx;
239 std::sort(sortedDedx.begin(), sortedDedx.end());
240
241 double truncatedMeanTmp = 0.0;
242 double meanTmp = 0.0;
243 double sumOfSquares = 0.0;
244 const int numDedx = sortedDedx.size();
245
246
247 for (int i = 0; i < numDedx; i++) {
248 meanTmp += sortedDedx[i];
249 }
250 if (numDedx != 0) {
251 meanTmp /= numDedx;
252 }
253
254 for (int i = 0; i < numDedx - 2; i++) {
255 truncatedMeanTmp += sortedDedx[i];
256 sumOfSquares += sortedDedx[i] * sortedDedx[i];
257 }
258 if (numDedx - 2 != 0) {
259 truncatedMeanTmp /= numDedx - 2;
260 }
261
262 mean = meanTmp;
263 truncatedMean = truncatedMeanTmp;
264
265 if (numDedx - 2 > 1) {
266 truncatedMeanErr = sqrt(sumOfSquares / double(numDedx - 2) - truncatedMeanTmp * truncatedMeanTmp) / double((numDedx - 2) - 1);
267 } else {
268 truncatedMeanErr = 0;
269 }
270}
271
272template <class HitClass> double VXDDedxPIDModule::getTraversedLength(const HitClass* hit, const RecoTrack* recoTrack, double& p)
273{
275 const VXD::SensorInfoBase& sensor = geo.getSensorInfo(hit->getSensorID());
276 const ROOT::Math::XYZVector& sensorNormal = sensor.vectorToGlobal(ROOT::Math::XYZVector(0.0, 0.0, 1.0));
277
278 const auto* hitInfo = recoTrack->getRecoHitInformation(hit);
279 if (not hitInfo) return 0;
280
281 const auto* trackPoint = recoTrack->getCreatedTrackPoint(hitInfo);
282 if (not trackPoint) return 0;
283
284 const auto* fitterInfo = trackPoint->getFitterInfo();
285 if (not fitterInfo) return 0;
286
287 double cosTheta = 0;
288 try {
289 const genfit::MeasuredStateOnPlane& mop = fitterInfo->getFittedState();
290 auto pocaMom = ROOT::Math::XYZVector(mop.getMom());
291 p = pocaMom.R();
292 cosTheta = sensorNormal.Dot(pocaMom.Unit());
293 } catch (genfit::Exception&) {
294 B2WARNING("recoTrack: " << recoTrack->getArrayIndex() << ": genfit::MeasuredStateOnPlane exception occured");
295 return 0;
296 }
297
298 return std::min(sensor.getWidth(), sensor.getThickness() / std::abs(cosTheta));
299}
300
301
302template <class HitClass> void VXDDedxPIDModule::saveSiHits(VXDDedxTrack* track, const std::vector<HitClass*>& hits,
303 const RecoTrack* recoTrack) const
304{
305 const int numHits = hits.size();
306 if (numHits == 0)
307 return;
308
310
311 //figure out which detector to assign hits to
312 const int currentDetector = geo.getSensorInfo(hits.front()->getSensorID()).getType();
313 assert(currentDetector == VXD::SensorInfoBase::PXD or currentDetector == VXD::SensorInfoBase::SVD);
314 // assert(currentDetector <= 1); //used as array index
315 assert(currentDetector == 0 or currentDetector == 1); //used as array index
316
317 std::vector<double> siliconDedx; //used for averages
318 siliconDedx.reserve(numHits);
319
320 VxdID prevSensor;
321 double psum = 0;
322 double p = 0;
323 for (int i = 0; i < numHits; i++) {
324 const HitClass* hit = hits[i];
325 if (!hit) {
326 B2ERROR("Added hit is a null pointer!");
327 continue;
328 }
329 const VxdID& currentSensor = hit->getSensorID();
330 int layer = -currentSensor.getLayerNumber();
331 assert(layer >= -6 && layer < 0);
332
333 //active medium traversed, in cm (can traverse one sensor at most)
334 double totalDistance = getTraversedLength(hit, recoTrack, p);
335 psum += p;
336 const double charge = hit->getCharge();
337 const double dedx = charge / totalDistance;
338 if (dedx <= 0) {
339 B2WARNING("dE/dx is " << dedx << " in layer " << layer);
340 } else if (i == 0 or prevSensor != currentSensor) { //only save once per sensor (u and v hits share charge!)
341 prevSensor = currentSensor;
342 //store data
343 if (totalDistance > 0) {
344 siliconDedx.push_back(dedx);
345 track->m_dedxAvg[currentDetector] += dedx;
346 track->addDedx(layer, totalDistance, dedx);
347 }
348 }
349
350 track->addHit(currentSensor, layer, charge, totalDistance, dedx);
351 track->m_p = psum / numHits;
352 }
353
354 //save averages
355 calculateMeans(track->m_dedxAvg[currentDetector],
356 track->m_dedxAvgTruncated[currentDetector],
357 track->m_dedxAvgTruncatedErr[currentDetector],
358 siliconDedx);
359}
360
static const ChargedStable pion
charged pion particle
Definition Const.h:661
A Class to store the Monte Carlo particle information.
Definition MCParticle.h:32
@ c_PrimaryParticle
bit 0: Particle is primary particle.
Definition MCParticle.h:47
bool hasStatus(unsigned short int bitmask) const
Return if specific status bit is set.
Definition MCParticle.h:118
int getPDG() const
Return PDG code of particle.
Definition MCParticle.h:101
ROOT::Math::XYZVector getMomentum() const
Return momentum.
Definition MCParticle.h:187
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
Definition Module.cc:208
Module()
Constructor.
Definition Module.cc:30
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
Definition Module.h:80
This is the Reconstruction Event-Data Model Track.
Definition RecoTrack.h:79
std::vector< Belle2::RecoTrack::UsedPXDHit * > getPXDHitList() const
Return an unsorted list of pxd hits.
Definition RecoTrack.h:449
std::vector< Belle2::RecoTrack::UsedSVDHit * > getSVDHitList() const
Return an unsorted list of svd hits.
Definition RecoTrack.h:452
bool hasTrackFitStatus(const genfit::AbsTrackRep *representation=nullptr) const
Check, if there is a fit status for the given representation or for the cardinal one.
Definition RecoTrack.cc:543
RecoHitInformation * getRecoHitInformation(HitType *hit) const
Return the reco hit information for a generic hit from the storeArray.
Definition RecoTrack.h:312
const genfit::TrackPoint * getCreatedTrackPoint(const RecoHitInformation *recoHitInformation) const
Get a pointer to the TrackPoint that was created from this hit.
Definition RecoTrack.cc:230
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
int getArrayIndex() const
Returns this object's array index (in StoreArray), or -1 if not found.
Values of the result of a track fit with a given particle hypothesis.
short getChargeSign() const
Return track charge (1 or -1).
ROOT::Math::XYZVector getMomentum() const
Getter for vector of momentum at closest approach of track in r/phi projection.
Container for likelihoods obtained by the VXD dE/dx PID (VXDDedxPIDModule).
DBObjPtr< SVDdEdxPDFs > m_SVDDedxPDFs
SVD DB object for dedx vs.
void calculateMeans(double &mean, double &truncatedMean, double &truncatedMeanErr, const std::vector< double > &dedx) const
Save arithmetic and truncated mean for the 'dedx' values.
int m_trackID
counter for tracks in this event
virtual void initialize() override
Initialize the module.
StoreArray< SVDCluster > m_svdClusters
Optional array of SVDClusters.
bool m_usePXD
use PXD data for likelihood
virtual void event() override
This method is called for each event.
StoreArray< MCParticle > m_mcparticles
Optional array of MCParticles.
StoreArray< VXDDedxLikelihood > m_dedxLikelihoods
Output array of VXDDedxLikelihoods.
void checkPDFs()
Check the pdfs for consistency every time they change in the database.
virtual void terminate() override
End of the event processing.
DBObjPtr< PXDdEdxPDFs > m_PXDDedxPDFs
PXD DB object for dedx vs.
StoreArray< VXDDedxTrack > m_dedxTracks
Output array of VXDDedxTracks.
virtual ~VXDDedxPIDModule()
Destructor.
int m_eventID
counter for events
bool m_useSVD
use SVD data for likelihood
void saveSiHits(VXDDedxTrack *track, const std::vector< HitClass * > &hits, const RecoTrack *recoTrack) const
save energy loss and hit information from SVD/PXDHits to track
StoreArray< Track > m_tracks
Required array of Tracks.
bool m_onlyPrimaryParticles
For MC only: if true, only save data for primary particles (as determined by MC truth)
static double getTraversedLength(const HitClass *hit, const RecoTrack *recoTrack, double &p)
returns traversed length through active medium of given hit
StoreArray< PXDCluster > m_pxdClusters
Optional array of PXDClusters.
StoreArray< RecoTrack > m_recoTracks
Required array of RecoTracks.
bool m_useIndividualHits
use individual hits (true) or truncated mean (false) to determine likelihoods
VXDDedxPIDModule()
Default constructor.
Debug output for VXDDedxPID module.
Class to facilitate easy access to sensor information of the VXD like coordinate transformations or p...
Definition GeoCache.h:39
const SensorInfoBase & getSensorInfo(Belle2::VxdID id) const
Return a reference to the SensorInfo of a given SensorID.
Definition GeoCache.cc:67
static GeoCache & getInstance()
Return a reference to the singleton instance.
Definition GeoCache.cc:214
Base class to provide Sensor Information for PXD and SVD.
SensorType getType() const
Return the Type of the Sensor.
Class to uniquely identify a any structure of the PXD and SVD.
Definition VxdID.h:33
baseType getLayerNumber() const
Get the layer id.
Definition VxdID.h:96
void addParam(const std::string &name, T &paramVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
Definition Module.h:559
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
Definition Module.h:649
double sqrt(double a)
sqrt for double
Definition beamHelpers.h:28
MCParticle * getMother() const
Returns a pointer to the mother particle.
Definition MCParticle.h:590
Abstract base class for different kinds of events.