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 <reconstruction/dataobjects/VXDDedxTrack.h>
16#include <reconstruction/dataobjects/VXDDedxLikelihood.h>
17#include <mdst/dataobjects/Track.h>
18#include <mdst/dataobjects/MCParticle.h>
19#include <tracking/dataobjects/RecoTrack.h>
20#include <svd/dataobjects/SVDCluster.h>
21#include <pxd/dataobjects/PXDCluster.h>
22#include <svd/dbobjects/SVDdEdxPDFs.h>
23#include <pxd/dbobjects/PXDdEdxPDFs.h>
24
25#include <genfit/MaterialEffects.h>
26
27#include <Math/VectorUtil.h>
28#include <TH2F.h>
29
30#include <memory>
31#include <cassert>
32#include <cmath>
33#include <algorithm>
34
35using namespace Belle2;
36using namespace Dedx;
37
38REG_MODULE(VXDDedxPID);
39
41{
43
44 //Set module properties
45 setDescription("Extract dE/dx and corresponding log-likelihood from fitted tracks and hits in the SVD and PXD.");
46
47 //Parameter definitions
48 addParam("useIndividualHits", m_useIndividualHits,
49 "Use individual hits (true) or truncated mean (false) to determine likelihoods", false);
50 addParam("onlyPrimaryParticles", m_onlyPrimaryParticles,
51 "For MC only: if true, only save data for primary particles (as determined by MC truth)", false);
52 addParam("usePXD", m_usePXD, "Use PXDClusters for dE/dx calculation", false);
53 addParam("useSVD", m_useSVD, "Use SVDClusters for dE/dx calculation", true);
54
55 m_eventID = -1;
56 m_trackID = 0;
57}
58
60
62{
63 if (m_usePXD) {
64 if (not m_PXDDedxPDFs) B2FATAL("No PXD dE/dx PDF's available");
65 bool ok = m_PXDDedxPDFs->checkPDFs(not m_useIndividualHits);
66 if (not ok) B2FATAL("Binning or ranges of PXD dE/dx PDF's differ");
67 }
68 if (m_useSVD) {
69 if (not m_SVDDedxPDFs) B2FATAL("No SVD Dedx PDF's available");
70 bool ok = m_SVDDedxPDFs->checkPDFs(not m_useIndividualHits);
71 if (not ok) B2FATAL("Binning or ranges of SVD dE/dx PDF's differ");
72 }
73}
74
76{
77
78 // required inputs
79 m_tracks.isRequired();
80 m_recoTracks.isRequired();
81
82 //optional inputs
83 m_mcparticles.isOptional();
84 m_tracks.optionalRelationTo(m_mcparticles);
85
86 if (m_useSVD)
87 m_svdClusters.isRequired();
88 else
89 m_svdClusters.isOptional();
90 if (m_usePXD)
91 m_pxdClusters.isRequired();
92 else
93 m_pxdClusters.isOptional();
94
95 // register dE/dx data points
96 m_dedxTracks.registerInDataStore();
97 m_tracks.registerRelationTo(m_dedxTracks);
98
99 // register likelihoods
100 m_dedxLikelihoods.registerInDataStore();
101 m_tracks.registerRelationTo(m_dedxLikelihoods);
102
103 m_SVDDedxPDFs.addCallback([this]() {checkPDFs();});
104 m_PXDDedxPDFs.addCallback([this]() {checkPDFs();});
105 checkPDFs();
106
107
108 // create instances here to not confuse profiling
110
111 if (!genfit::MaterialEffects::getInstance()->isInitialized()) {
112 B2FATAL("Need to have SetupGenfitExtrapolationModule in path before this one.");
113 }
114}
115
117{
118 m_dedxTracks.clear();
119 m_dedxLikelihoods.clear();
120
121 // go through Tracks
122 // get fitresult and gftrack and do extrapolations, save corresponding dE/dx and likelihood values
123 // get genfit::TrackCand through genfit::Track::getCand()
124 // get hit indices through genfit::TrackCand::getHit(...)
125 // create one VXDDedxTrack per fitresult/gftrack
126 //create one DedkLikelihood per Track (plus rel)
127 m_eventID++;
128
129 // inputs
130 const int numMCParticles = m_mcparticles.getEntries();
131
132 // **************************************************
133 //
134 // LOOP OVER TRACKS
135 //
136 // **************************************************
137
138 for (const auto& track : m_tracks) {
139 m_trackID++;
140
141 std::shared_ptr<VXDDedxTrack> dedxTrack = std::make_shared<VXDDedxTrack>();
142 dedxTrack->m_eventID = m_eventID;
143 dedxTrack->m_trackID = m_trackID;
144
145 // load the pion fit hypothesis or the hypothesis which is the closest in mass to a pion
146 // the tracking will not always successfully fit with a pion hypothesis
147 const TrackFitResult* fitResult = track.getTrackFitResultWithClosestMass(Const::pion);
148 if (!fitResult) {
149 B2WARNING("No related fit for track ...");
150 continue;
151 }
152
153 if (numMCParticles != 0) {
154 // find MCParticle corresponding to this track
155 const MCParticle* mcpart = track.getRelatedTo<MCParticle>();
156
157 if (mcpart) {
159 continue; //not a primary particle, ignore
160 }
161
162 //add some MC truths to VXDDedxTrack object
163 dedxTrack->m_pdg = mcpart->getPDG();
164 const MCParticle* mother = mcpart->getMother();
165 dedxTrack->m_motherPDG = mother ? mother->getPDG() : 0;
166
167 const ROOT::Math::XYZVector trueMomentum = mcpart->getMomentum();
168 dedxTrack->m_pTrue = trueMomentum.R();
169 }
170 }
171
172 // get momentum (at origin) from fit result
173 const ROOT::Math::XYZVector& trackMom = fitResult->getMomentum();
174 dedxTrack->m_p = trackMom.R();
175 dedxTrack->m_cosTheta = cos(trackMom.Theta());
176 dedxTrack->m_charge = fitResult->getChargeSign();
177
178 // dE/dx values will be calculated using associated RecoTrack
179 const RecoTrack* recoTrack = track.getRelatedTo<RecoTrack>();
180 if (!recoTrack) {
181 B2WARNING("No related track for this fit...");
182 continue;
183 }
184
185 // Check to see if the track is pruned. We use the cardinal representation for this, as we do not expect that one
186 // representation is pruned whereas the other is not.
187 if (recoTrack->hasTrackFitStatus()) {
188 if (recoTrack->getTrackFitStatus()->isTrackPruned()) {
189 B2ERROR("GFTrack is pruned, please run VXDDedxPID only on unpruned tracks! Skipping this track.");
190 continue;
191 }
192 } else // The RecoTrack might not have a fit status for reasons: let's skip it
193 continue;
194
195 if (m_usePXD) {
196 const std::vector<PXDCluster*>& pxdClusters = recoTrack->getPXDHitList();
197 saveSiHits(dedxTrack.get(), pxdClusters, recoTrack);
198 }
199
200 if (m_useSVD) {
201 const std::vector<SVDCluster*>& svdClusters = recoTrack->getSVDHitList();
202 saveSiHits(dedxTrack.get(), svdClusters, recoTrack);
203 }
204
205 if (dedxTrack->dedx.empty()) {
206 B2DEBUG(50, "Found track with no hits, ignoring.");
207 continue;
208 }
209
210 // add a few last things to the VXDDedxTrack
211 const int numDedx = dedxTrack->dedx.size();
212 dedxTrack->m_nHits = numDedx;
213 // no need to define lowedgetruncated and highedgetruncated as we always remove the highest 2 dE/dx values from 8 dE/dx value
214 dedxTrack->m_nHitsUsed = numDedx - 2;
215
216 // calculate log likelihoods
217 dedxTrack->clearLogLikelihoods();
218 bool truncated = not m_useIndividualHits;
219 if (m_usePXD) dedxTrack->addLogLikelihoods(m_PXDDedxPDFs->getPDFs(truncated), Dedx::c_PXD, truncated);
220 if (m_useSVD) dedxTrack->addLogLikelihoods(m_SVDDedxPDFs->getPDFs(truncated), Dedx::c_SVD, truncated);
221
222 // save log likelihoods
223 if (dedxTrack->areLogLikelihoodsAvailable()) {
224 VXDDedxLikelihood* likelihoodObj = m_dedxLikelihoods.appendNew(dedxTrack->m_vxdLogl);
225 track.addRelationTo(likelihoodObj);
226 }
227
228 // save the VXDDedxTrack
229 VXDDedxTrack* newVXDDedxTrack = m_dedxTracks.appendNew(*dedxTrack);
230 track.addRelationTo(newVXDDedxTrack);
231
232 } // end of loop over tracks
233}
234
236{
237
238 B2DEBUG(50, "VXDDedxPIDModule exiting after processing " << m_trackID <<
239 " tracks in " << m_eventID + 1 << " events.");
240}
241
242
243// calculateMeans need some change as we always remove highest 2 dE/dx values
244void VXDDedxPIDModule::calculateMeans(double& mean, double& truncatedMean, double& truncatedMeanErr,
245 const std::vector<double>& dedx) const
246{
247 // Calculate the truncated average by skipping only highest two value
248 std::vector<double> sortedDedx = dedx;
249 std::sort(sortedDedx.begin(), sortedDedx.end());
250
251 double truncatedMeanTmp = 0.0;
252 double meanTmp = 0.0;
253 double sumOfSquares = 0.0;
254 const int numDedx = sortedDedx.size();
255
256
257 for (int i = 0; i < numDedx; i++) {
258 meanTmp += sortedDedx[i];
259 }
260 if (numDedx != 0) {
261 meanTmp /= numDedx;
262 }
263
264 for (int i = 0; i < numDedx - 2; i++) {
265 truncatedMeanTmp += sortedDedx[i];
266 sumOfSquares += sortedDedx[i] * sortedDedx[i];
267 }
268 if (numDedx - 2 != 0) {
269 truncatedMeanTmp /= numDedx - 2;
270 }
271
272 mean = meanTmp;
273 truncatedMean = truncatedMeanTmp;
274
275 if (numDedx - 2 > 1) {
276 truncatedMeanErr = sqrt(sumOfSquares / double(numDedx - 2) - truncatedMeanTmp * truncatedMeanTmp) / double((numDedx - 2) - 1);
277 } else {
278 truncatedMeanErr = 0;
279 }
280}
281
282template <class HitClass> double VXDDedxPIDModule::getTraversedLength(const HitClass* hit, const RecoTrack* recoTrack, double& p)
283{
285 const VXD::SensorInfoBase& sensor = geo.getSensorInfo(hit->getSensorID());
286 const ROOT::Math::XYZVector& sensorNormal = sensor.vectorToGlobal(ROOT::Math::XYZVector(0.0, 0.0, 1.0));
287
288 const auto* hitInfo = recoTrack->getRecoHitInformation(hit);
289 if (not hitInfo) return 0;
290
291 const auto* trackPoint = recoTrack->getCreatedTrackPoint(hitInfo);
292 if (not trackPoint) return 0;
293
294 const auto* fitterInfo = trackPoint->getFitterInfo();
295 if (not fitterInfo) return 0;
296
297 double cosTheta = 0;
298 try {
299 const genfit::MeasuredStateOnPlane& mop = fitterInfo->getFittedState();
300 auto pocaMom = ROOT::Math::XYZVector(mop.getMom());
301 p = pocaMom.R();
302 cosTheta = sensorNormal.Dot(pocaMom.Unit());
303 } catch (genfit::Exception&) {
304 B2WARNING("recoTrack: " << recoTrack->getArrayIndex() << ": genfit::MeasuredStateOnPlane exception occured");
305 return 0;
306 }
307
308 return std::min(sensor.getWidth(), sensor.getThickness() / std::abs(cosTheta));
309}
310
311
312template <class HitClass> void VXDDedxPIDModule::saveSiHits(VXDDedxTrack* track, const std::vector<HitClass*>& hits,
313 const RecoTrack* recoTrack) const
314{
315 const int numHits = hits.size();
316 if (numHits == 0)
317 return;
318
320
321 //figure out which detector to assign hits to
322 const int currentDetector = geo.getSensorInfo(hits.front()->getSensorID()).getType();
323 assert(currentDetector == VXD::SensorInfoBase::PXD or currentDetector == VXD::SensorInfoBase::SVD);
324 // assert(currentDetector <= 1); //used as array index
325 assert(currentDetector == 0 or currentDetector == 1); //used as array index
326
327 std::vector<double> siliconDedx; //used for averages
328 siliconDedx.reserve(numHits);
329
330 VxdID prevSensor;
331 double psum = 0;
332 double p = 0;
333 for (int i = 0; i < numHits; i++) {
334 const HitClass* hit = hits[i];
335 if (!hit) {
336 B2ERROR("Added hit is a null pointer!");
337 continue;
338 }
339 const VxdID& currentSensor = hit->getSensorID();
340 int layer = -currentSensor.getLayerNumber();
341 assert(layer >= -6 && layer < 0);
342
343 //active medium traversed, in cm (can traverse one sensor at most)
344 double totalDistance = getTraversedLength(hit, recoTrack, p);
345 psum += p;
346 const double charge = hit->getCharge();
347 const double dedx = charge / totalDistance;
348 if (dedx <= 0) {
349 B2WARNING("dE/dx is " << dedx << " in layer " << layer);
350 } else if (i == 0 or prevSensor != currentSensor) { //only save once per sensor (u and v hits share charge!)
351 prevSensor = currentSensor;
352 //store data
353 if (totalDistance > 0) {
354 siliconDedx.push_back(dedx);
355 track->m_dedxAvg[currentDetector] += dedx;
356 track->addDedx(layer, totalDistance, dedx);
357 }
358 }
359
360 track->addHit(currentSensor, layer, charge, totalDistance, dedx);
361 track->m_p = psum / numHits;
362 }
363
364 //save averages
365 calculateMeans(track->m_dedxAvg[currentDetector],
366 track->m_dedxAvgTruncated[currentDetector],
367 track->m_dedxAvgTruncatedErr[currentDetector],
368 siliconDedx);
369}
370
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:38
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:32
baseType getLayerNumber() const
Get the layer id.
Definition VxdID.h:95
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.