9#include <reconstruction/modules/VXDDedxPID/VXDDedxPIDModule.h>
11#include <framework/gearbox/Const.h>
13#include <vxd/geometry/GeoCache.h>
15#include <genfit/MaterialEffects.h>
17#include <Math/VectorUtil.h>
35 setDescription(
"Extract dE/dx and corresponding log-likelihood from fitted tracks and hits in the SVD and PXD.");
39 "Use individual hits (true) or truncated mean (false) to determine likelihoods",
false);
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);
54 if (not
m_PXDDedxPDFs) B2FATAL(
"No PXD dE/dx PDF's available");
56 if (not ok) B2FATAL(
"Binning or ranges of PXD dE/dx PDF's differ");
59 if (not
m_SVDDedxPDFs) B2FATAL(
"No SVD Dedx PDF's available");
61 if (not ok) B2FATAL(
"Binning or ranges of SVD dE/dx PDF's differ");
101 if (!genfit::MaterialEffects::getInstance()->isInitialized()) {
102 B2FATAL(
"Need to have SetupGenfitExtrapolationModule in path before this one.");
128 for (
const auto& track :
m_tracks) {
131 std::shared_ptr<VXDDedxTrack> dedxTrack = std::make_shared<VXDDedxTrack>();
139 B2WARNING(
"No related fit for track ...");
143 if (numMCParticles != 0) {
153 dedxTrack->m_pdg = mcpart->
getPDG();
155 dedxTrack->m_motherPDG = mother ? mother->
getPDG() : 0;
157 const ROOT::Math::XYZVector trueMomentum = mcpart->
getMomentum();
158 dedxTrack->m_pTrue = trueMomentum.R();
163 const ROOT::Math::XYZVector& trackMom = fitResult->
getMomentum();
164 dedxTrack->m_p = trackMom.R();
165 dedxTrack->m_cosTheta = cos(trackMom.Theta());
171 B2WARNING(
"No related track for this fit...");
179 B2ERROR(
"GFTrack is pruned, please run VXDDedxPID only on unpruned tracks! Skipping this track.");
186 const std::vector<PXDCluster*>& pxdClusters = recoTrack->
getPXDHitList();
187 saveSiHits(dedxTrack.get(), pxdClusters, recoTrack);
191 const std::vector<SVDCluster*>& svdClusters = recoTrack->
getSVDHitList();
192 saveSiHits(dedxTrack.get(), svdClusters, recoTrack);
195 if (dedxTrack->dedx.empty()) {
196 B2DEBUG(50,
"Found track with no hits, ignoring.");
201 const int numDedx = dedxTrack->dedx.size();
202 dedxTrack->m_nHits = numDedx;
204 dedxTrack->m_nHitsUsed = numDedx - 2;
207 dedxTrack->clearLogLikelihoods();
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);
213 if (dedxTrack->areLogLikelihoodsAvailable()) {
215 track.addRelationTo(likelihoodObj);
220 track.addRelationTo(newVXDDedxTrack);
228 B2DEBUG(50,
"VXDDedxPIDModule exiting after processing " <<
m_trackID <<
229 " tracks in " <<
m_eventID + 1 <<
" events.");
235 const std::vector<double>& dedx)
const
238 std::vector<double> sortedDedx = dedx;
239 std::sort(sortedDedx.begin(), sortedDedx.end());
241 double truncatedMeanTmp = 0.0;
242 double meanTmp = 0.0;
243 double sumOfSquares = 0.0;
244 const int numDedx = sortedDedx.size();
247 for (
int i = 0; i < numDedx; i++) {
248 meanTmp += sortedDedx[i];
254 for (
int i = 0; i < numDedx - 2; i++) {
255 truncatedMeanTmp += sortedDedx[i];
256 sumOfSquares += sortedDedx[i] * sortedDedx[i];
258 if (numDedx - 2 != 0) {
259 truncatedMeanTmp /= numDedx - 2;
263 truncatedMean = truncatedMeanTmp;
265 if (numDedx - 2 > 1) {
266 truncatedMeanErr =
sqrt(sumOfSquares /
double(numDedx - 2) - truncatedMeanTmp * truncatedMeanTmp) / double((numDedx - 2) - 1);
268 truncatedMeanErr = 0;
276 const ROOT::Math::XYZVector& sensorNormal = sensor.vectorToGlobal(ROOT::Math::XYZVector(0.0, 0.0, 1.0));
279 if (not hitInfo)
return 0;
282 if (not trackPoint)
return 0;
284 const auto* fitterInfo = trackPoint->getFitterInfo();
285 if (not fitterInfo)
return 0;
289 const genfit::MeasuredStateOnPlane& mop = fitterInfo->getFittedState();
290 auto pocaMom = ROOT::Math::XYZVector(mop.getMom());
292 cosTheta = sensorNormal.Dot(pocaMom.Unit());
293 }
catch (genfit::Exception&) {
294 B2WARNING(
"recoTrack: " << recoTrack->
getArrayIndex() <<
": genfit::MeasuredStateOnPlane exception occured");
298 return std::min(sensor.getWidth(), sensor.getThickness() / std::abs(cosTheta));
305 const int numHits = hits.size();
315 assert(currentDetector == 0 or currentDetector == 1);
317 std::vector<double> siliconDedx;
318 siliconDedx.reserve(numHits);
323 for (
int i = 0; i < numHits; i++) {
324 const HitClass* hit = hits[i];
326 B2ERROR(
"Added hit is a null pointer!");
329 const VxdID& currentSensor = hit->getSensorID();
331 assert(layer >= -6 && layer < 0);
336 const double charge = hit->getCharge();
337 const double dedx = charge / totalDistance;
339 B2WARNING(
"dE/dx is " << dedx <<
" in layer " << layer);
340 }
else if (i == 0 or prevSensor != currentSensor) {
341 prevSensor = currentSensor;
343 if (totalDistance > 0) {
344 siliconDedx.push_back(dedx);
345 track->m_dedxAvg[currentDetector] += dedx;
346 track->addDedx(layer, totalDistance, dedx);
350 track->addHit(currentSensor, layer, charge, totalDistance, dedx);
351 track->m_p = psum / numHits;
356 track->m_dedxAvgTruncated[currentDetector],
357 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.