9#include <reconstruction/modules/VXDDedxPID/VXDDedxPIDModule.h>
11#include <framework/gearbox/Const.h>
13#include <vxd/geometry/GeoCache.h>
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>
25#include <genfit/MaterialEffects.h>
27#include <Math/VectorUtil.h>
45 setDescription(
"Extract dE/dx and corresponding log-likelihood from fitted tracks and hits in the SVD and PXD.");
49 "Use individual hits (true) or truncated mean (false) to determine likelihoods",
false);
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);
64 if (not
m_PXDDedxPDFs) B2FATAL(
"No PXD dE/dx PDF's available");
66 if (not ok) B2FATAL(
"Binning or ranges of PXD dE/dx PDF's differ");
69 if (not
m_SVDDedxPDFs) B2FATAL(
"No SVD Dedx PDF's available");
71 if (not ok) B2FATAL(
"Binning or ranges of SVD dE/dx PDF's differ");
111 if (!genfit::MaterialEffects::getInstance()->isInitialized()) {
112 B2FATAL(
"Need to have SetupGenfitExtrapolationModule in path before this one.");
138 for (
const auto& track :
m_tracks) {
141 std::shared_ptr<VXDDedxTrack> dedxTrack = std::make_shared<VXDDedxTrack>();
149 B2WARNING(
"No related fit for track ...");
153 if (numMCParticles != 0) {
163 dedxTrack->m_pdg = mcpart->
getPDG();
165 dedxTrack->m_motherPDG = mother ? mother->
getPDG() : 0;
167 const ROOT::Math::XYZVector trueMomentum = mcpart->
getMomentum();
168 dedxTrack->m_pTrue = trueMomentum.R();
173 const ROOT::Math::XYZVector& trackMom = fitResult->
getMomentum();
174 dedxTrack->m_p = trackMom.R();
175 dedxTrack->m_cosTheta = cos(trackMom.Theta());
181 B2WARNING(
"No related track for this fit...");
189 B2ERROR(
"GFTrack is pruned, please run VXDDedxPID only on unpruned tracks! Skipping this track.");
196 const std::vector<PXDCluster*>& pxdClusters = recoTrack->
getPXDHitList();
197 saveSiHits(dedxTrack.get(), pxdClusters, recoTrack);
201 const std::vector<SVDCluster*>& svdClusters = recoTrack->
getSVDHitList();
202 saveSiHits(dedxTrack.get(), svdClusters, recoTrack);
205 if (dedxTrack->dedx.empty()) {
206 B2DEBUG(50,
"Found track with no hits, ignoring.");
211 const int numDedx = dedxTrack->dedx.size();
212 dedxTrack->m_nHits = numDedx;
214 dedxTrack->m_nHitsUsed = numDedx - 2;
217 dedxTrack->clearLogLikelihoods();
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);
223 if (dedxTrack->areLogLikelihoodsAvailable()) {
225 track.addRelationTo(likelihoodObj);
230 track.addRelationTo(newVXDDedxTrack);
238 B2DEBUG(50,
"VXDDedxPIDModule exiting after processing " <<
m_trackID <<
239 " tracks in " <<
m_eventID + 1 <<
" events.");
245 const std::vector<double>& dedx)
const
248 std::vector<double> sortedDedx = dedx;
249 std::sort(sortedDedx.begin(), sortedDedx.end());
251 double truncatedMeanTmp = 0.0;
252 double meanTmp = 0.0;
253 double sumOfSquares = 0.0;
254 const int numDedx = sortedDedx.size();
257 for (
int i = 0; i < numDedx; i++) {
258 meanTmp += sortedDedx[i];
264 for (
int i = 0; i < numDedx - 2; i++) {
265 truncatedMeanTmp += sortedDedx[i];
266 sumOfSquares += sortedDedx[i] * sortedDedx[i];
268 if (numDedx - 2 != 0) {
269 truncatedMeanTmp /= numDedx - 2;
273 truncatedMean = truncatedMeanTmp;
275 if (numDedx - 2 > 1) {
276 truncatedMeanErr =
sqrt(sumOfSquares /
double(numDedx - 2) - truncatedMeanTmp * truncatedMeanTmp) / double((numDedx - 2) - 1);
278 truncatedMeanErr = 0;
286 const ROOT::Math::XYZVector& sensorNormal = sensor.vectorToGlobal(ROOT::Math::XYZVector(0.0, 0.0, 1.0));
289 if (not hitInfo)
return 0;
292 if (not trackPoint)
return 0;
294 const auto* fitterInfo = trackPoint->getFitterInfo();
295 if (not fitterInfo)
return 0;
299 const genfit::MeasuredStateOnPlane& mop = fitterInfo->getFittedState();
300 auto pocaMom = ROOT::Math::XYZVector(mop.getMom());
302 cosTheta = sensorNormal.Dot(pocaMom.Unit());
303 }
catch (genfit::Exception&) {
304 B2WARNING(
"recoTrack: " << recoTrack->
getArrayIndex() <<
": genfit::MeasuredStateOnPlane exception occured");
308 return std::min(sensor.getWidth(), sensor.getThickness() / std::abs(cosTheta));
315 const int numHits = hits.size();
325 assert(currentDetector == 0 or currentDetector == 1);
327 std::vector<double> siliconDedx;
328 siliconDedx.reserve(numHits);
333 for (
int i = 0; i < numHits; i++) {
334 const HitClass* hit = hits[i];
336 B2ERROR(
"Added hit is a null pointer!");
339 const VxdID& currentSensor = hit->getSensorID();
341 assert(layer >= -6 && layer < 0);
346 const double charge = hit->getCharge();
347 const double dedx = charge / totalDistance;
349 B2WARNING(
"dE/dx is " << dedx <<
" in layer " << layer);
350 }
else if (i == 0 or prevSensor != currentSensor) {
351 prevSensor = currentSensor;
353 if (totalDistance > 0) {
354 siliconDedx.push_back(dedx);
355 track->m_dedxAvg[currentDetector] += dedx;
356 track->addDedx(layer, totalDistance, dedx);
360 track->addHit(currentSensor, layer, charge, totalDistance, dedx);
361 track->m_p = psum / numHits;
366 track->m_dedxAvgTruncated[currentDetector],
367 track->m_dedxAvgTruncatedErr[currentDetector],
static const ChargedStable pion
charged pion particle
A Class to store the Monte Carlo particle information.
@ c_PrimaryParticle
bit 0: Particle is primary particle.
bool hasStatus(unsigned short int bitmask) const
Return if specific status bit is set.
int getPDG() const
Return PDG code of particle.
ROOT::Math::XYZVector getMomentum() const
Return momentum.
void setDescription(const std::string &description)
Sets the description of the module.
void setPropertyFlags(unsigned int propertyFlags)
Sets the flags for the module properties.
@ c_ParallelProcessingCertified
This module can be run in parallel processing mode safely (All I/O must be done through the data stor...
This is the Reconstruction Event-Data Model Track.
std::vector< Belle2::RecoTrack::UsedPXDHit * > getPXDHitList() const
Return an unsorted list of pxd hits.
std::vector< Belle2::RecoTrack::UsedSVDHit * > getSVDHitList() const
Return an unsorted list of svd hits.
bool hasTrackFitStatus(const genfit::AbsTrackRep *representation=nullptr) const
Check, if there is a fit status for the given representation or for the cardinal one.
RecoHitInformation * getRecoHitInformation(HitType *hit) const
Return the reco hit information for a generic hit from the storeArray.
const genfit::TrackPoint * getCreatedTrackPoint(const RecoHitInformation *recoHitInformation) const
Get a pointer to the TrackPoint that was created from this hit.
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...
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...
const SensorInfoBase & getSensorInfo(Belle2::VxdID id) const
Return a reference to the SensorInfo of a given SensorID.
static GeoCache & getInstance()
Return a reference to the singleton instance.
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.
baseType getLayerNumber() const
Get the layer id.
void addParam(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module.
#define REG_MODULE(moduleName)
Register the given module (without 'Module' suffix) with the framework.
double sqrt(double a)
sqrt for double
MCParticle * getMother() const
Returns a pointer to the mother particle.
Abstract base class for different kinds of events.