Belle II Software development
PXDDQMEfficiencySelftrackModule.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 <pxd/modules/pxdDQM/PXDDQMEfficiencySelftrackModule.h>
10#include <pxd/dataobjects/PXDCluster.h>
11#include <tracking/dataobjects/ROIid.h>
12#include <tracking/dataobjects/RecoTrack.h>
13
14#include <pxd/reconstruction/PXDPixelMasker.h>
15#include <mdst/dataobjects/Track.h>
16#include <framework/gearbox/Const.h>
17#include <framework/geometry/VectorUtil.h>
18
19#include "TDirectory.h"
20#include "TMatrixDSym.h"
21
22using namespace Belle2;
23
24//-----------------------------------------------------------------
25// Register the Module
26//-----------------------------------------------------------------
27REG_MODULE(PXDDQMEfficiencySelftrack);
28
29//-----------------------------------------------------------------
30// Implementation
31//-----------------------------------------------------------------
32
34{
35 // Set module properties
36 setDescription("Create basic histograms for PXD efficiency");
37
38 // What exactly is needed for this to be true?
40
41 // Parameter definitions
42 addParam("pxdClustersName", m_pxdClustersName, "name of StoreArray with PXD cluster", std::string(""));
43 addParam("tracksName", m_tracksName, "name of StoreArray with RecoTracks", std::string(""));
44 addParam("ROIsName", m_ROIsName, "name of the list of HLT ROIs, if available in output", std::string(""));
45 addParam("histogramDirectoryName", m_histogramDirectoryName, "Name of the directory where histograms will be placed",
46 std::string("PXDEFF"));
47 addParam("binsU", m_u_bins, "histogram bins in u direction", int(4));
48 addParam("binsV", m_v_bins, "histogram bins in v direction", int(6));
49 addParam("distCut", m_distcut, "max distance in [cm] for cluster to be counted to a track", double(0.0500));
50 addParam("pCut", m_pcut, "Set a cut on the track p-value", double(1e-20));
51 addParam("requireROIs", m_requireROIs, "require tracks to lie inside a ROI", bool(false));
52 addParam("useAlignment", m_useAlignment, "if true the alignment will be used", bool(true));
53 addParam("maskDeadPixels", m_maskDeadPixels, "Do not consider tracks going through known dead or hot pixels for the efficiency",
54 bool(false));
55 addParam("minSVDHits", m_minSVDHits, "Number of SVD hits required in a track to be considered", 0u);
56 addParam("momCut", m_momCut, "Set a cut on the track momentum in GeV/c, 0 disables", double(0));
57 addParam("pTCut", m_pTCut, "Set a cut on the track pT in GeV/c, 0 disables", double(1));
58 addParam("cutBorders", m_cutBorders, "Do not use tracks near the borders of the sensor", bool(true));
59 addParam("maskedDistance", m_maskedDistance, "Distance inside which no masked pixel or sensor border is allowed", int(10));
60 addParam("trackUFactorDistCut", m_uFactor, "Set a cut on u error of track (factor*err<dist), 0 disables", double(2.0));
61 addParam("trackVFactorDistCut", m_vFactor, "Set a cut on v error of track (factor*err<dist), 0 disables", double(2.0));
62 addParam("z0minCut", m_z0minCut, "Set a cut z0 minimum in cm (large negative value eg -9999 disables)", double(-1));
63 addParam("z0maxCut", m_z0maxCut, "Set a cut z0 maximum in cm (large positive value eg 9999 disables)", double(1));
64 addParam("d0Cut", m_d0Cut, "Set a cut abs(d0) in cm (large positive value eg 9999 disables)", double(0.5));
65 addParam("verboseHistos", m_verboseHistos, "Add more verbose histograms for cuts (not for ereoc)", bool(false));
66}
67
68
70{
71 //calls the define histogram function
72 REG_HISTOGRAM;
73
74 //register the required arrays
75 //Register as optional so validation for cases where they are not available still succeeds, but module will not do any meaningful work without them
77 m_tracks.isOptional(m_tracksName);
78 m_ROIs.isOptional(m_ROIsName);
79}
80
82{
83 for (auto& h : m_h_track_hits) if (h.second) h.second->Reset();
84 for (auto& h : m_h_matched_cluster) if (h.second) h.second->Reset();
85 for (auto& h : m_h_p) if (h.second) h.second->Reset();
86 for (auto& h : m_h_pt) if (h.second) h.second->Reset();
87 for (auto& h : m_h_su) if (h.second) h.second->Reset();
88 for (auto& h : m_h_sv) if (h.second) h.second->Reset();
89 for (auto& h : m_h_p2) if (h.second) h.second->Reset();
90 for (auto& h : m_h_pt2) if (h.second) h.second->Reset();
91 for (auto& h : m_h_su2) if (h.second) h.second->Reset();
92 for (auto& h : m_h_sv2) if (h.second) h.second->Reset();
93}
94
96{
97 if (!m_pxdclusters.isValid()) {
98 B2INFO("PXDClusters array is missing, no efficiencies");
99 return;
100 }
101 if (!m_tracks.isValid()) {
102 B2INFO("RecoTrack array is missing, no efficiencies");
103 return;
104 }
105 if (!m_ROIs.isValid() && m_requireROIs) {
106 B2INFO("ROI array is missing but required hits in ROIs, aborting");
107 return;
108 }
109
110
111 for (auto& a_track : m_tracks) {
112
113 //If fit failed assume position pointed to is useless anyway
114 if (!a_track.wasFitSuccessful()) continue;
115
116 if (a_track.getNumberOfSVDHits() < m_minSVDHits) continue;
117
118 const genfit::FitStatus* fitstatus = a_track.getTrackFitStatus();
119 if (fitstatus->getPVal() < m_pcut) continue;
120
121 genfit::MeasuredStateOnPlane trackstate;
122 trackstate = a_track.getMeasuredStateOnPlaneFromFirstHit();
123 if (trackstate.getMom().Mag() < m_momCut) continue;
124 if (trackstate.getMom().Pt() < m_pTCut) continue;
125
126 auto ptr = a_track.getRelated<Track>("Tracks");
127
128 if (!ptr) {
129 B2ERROR("expect a track for fitted recotracks");
130 continue;
131 }
133 if (!ptr2) {
134 B2ERROR("expect a track fit result for mass");
135 continue;
136 }
137
138 // Vertex cut
139 if (ptr2->getZ0() < m_z0minCut || ptr2->getZ0() > m_z0maxCut || fabs(ptr2->getD0()) > m_d0Cut) continue;
140
141 //loop over all PXD sensors to get the intersections
142 std::vector<VxdID> sensors = m_vxdGeometry.getListOfSensors();
143 for (VxdID& aVxdID : sensors) {
144 VXD::SensorInfoBase info = m_vxdGeometry.getSensorInfo(aVxdID);
145 if (info.getType() != VXD::SensorInfoBase::PXD) continue;
146 //Search for intersections of the track with all PXD layers
147 //Traditional (aka the person before did it like this) method
148 //If there is a way to find out sensors crossed by a track directly, that would most likely be faster
149
150 bool isgood = false;
151 //true = track intersects current sensor
152 double sigu(-9999);
153 double sigv(-9999);
154 ROOT::Math::XYZVector intersec_buff = getTrackInterSec(info, a_track, isgood, sigu, sigv);
155
156 if (!isgood) {
157 continue;//track does not go through this sensor-> nothing to measure anyway
158 } else {
159 if (m_verboseHistos) {
160 if (m_h_p[aVxdID]) m_h_p[aVxdID]->Fill(trackstate.getMom().Mag());
161 if (m_h_pt[aVxdID]) m_h_pt[aVxdID]->Fill(trackstate.getMom().Pt());
162 if (m_h_su[aVxdID]) m_h_su[aVxdID]->Fill(sigu);
163 if (m_h_sv[aVxdID]) m_h_sv[aVxdID]->Fill(sigv);
164 }
165 if (m_uFactor * sigu > m_distcut) continue; // Error ufak*SigmaU > cut
166 if (m_vFactor * sigv > m_distcut) continue; // Error vfak*SigmaV > cut
167
168 double u_fit = intersec_buff.X();
169 double v_fit = intersec_buff.Y();
170
171 int ucell_fit = info.getUCellID(intersec_buff.X());
172 int vcell_fit = info.getVCellID(intersec_buff.Y());
173
174 if (m_cutBorders && isCloseToBorder(ucell_fit, vcell_fit, m_maskedDistance)) {
175 continue;
176 }
177
178 if (m_maskDeadPixels && isDeadPixelClose(ucell_fit, vcell_fit, m_maskedDistance, aVxdID)) {
179 continue;
180 }
181
182 if (m_requireROIs) {
183 //Check if the intersection is inside a ROI
184 //If not, even if measured the cluster was thrown away->Not PXD's fault
185 bool fitInsideROI = false;
186 for (auto& roit : m_ROIs) {
187 if (aVxdID != roit.getSensorID()) {
188 continue; //ROI on other sensor
189 }
190
191 if (ucell_fit < roit.getMaxUid()
192 && ucell_fit > roit.getMinUid()
193 && vcell_fit < roit.getMaxVid()
194 && vcell_fit > roit.getMinVid()) {
195 fitInsideROI = true;
196 }
197 }
198 if (!fitInsideROI) {
199 continue;//Hit wouldn't have been recorded
200 }
201 }
202
203 //This track should be on the sensor
204 m_h_track_hits[aVxdID]->Fill(ucell_fit, vcell_fit);
205
206 //Now check if the sensor measured a hit here
207
208 int bestcluster = findClosestCluster(aVxdID, intersec_buff);
209 if (bestcluster >= 0) {
210 double u_clus = m_pxdclusters[bestcluster]->getU();
211 double v_clus = m_pxdclusters[bestcluster]->getV();
212
213 //is the closest cluster close enough to the track to count as measured?
214 ROOT::Math::XYZVector dist_clus(u_fit - u_clus, v_fit - v_clus, 0);
215 if (dist_clus.R() <= m_distcut) {
216 m_h_matched_cluster[aVxdID]->Fill(ucell_fit, vcell_fit);
217 if (m_verboseHistos) {
218 if (m_h_p2[aVxdID]) m_h_p2[aVxdID]->Fill(trackstate.getMom().Mag());
219 if (m_h_pt2[aVxdID]) m_h_pt2[aVxdID]->Fill(trackstate.getMom().Pt());
220 if (m_h_su2[aVxdID]) m_h_su2[aVxdID]->Fill(sigu);
221 if (m_h_sv2[aVxdID]) m_h_sv2[aVxdID]->Fill(sigv);
222 }
223 }
224 }
225 }
226 }
227 }
228}
229
230
231
232ROOT::Math::XYZVector
234 bool& isgood, double& du, double& dv)
235{
236 //will be set true if the intersect was found
237 isgood = false;
238
239 ROOT::Math::XYZVector intersec(99999999, 9999999, 0); //point outside the sensor
240
241 genfit::MeasuredStateOnPlane gfTrackState = aTrack.getMeasuredStateOnPlaneFromFirstHit();
242
243 //adopted (aka stolen) from tracking/modules/pxdClusterRescue/PXDClusterRescueROIModule
244 try {
245 // get sensor plane
246 ROOT::Math::XYZVector zeroVec(0, 0, 0);
247 ROOT::Math::XYZVector uVec(1, 0, 0);
248 ROOT::Math::XYZVector vVec(0, 1, 0);
249
250 genfit::DetPlane* sensorPlane = new genfit::DetPlane();
251 sensorPlane->setO(XYZToTVector(pxdSensorInfo.pointToGlobal(zeroVec, m_useAlignment)));
252 sensorPlane->setUV(XYZToTVector(pxdSensorInfo.vectorToGlobal(uVec, m_useAlignment)),
253 XYZToTVector(pxdSensorInfo.vectorToGlobal(vVec, m_useAlignment)));
254
255 //boost pointer (will be deleted automatically ?!?!?)
256 genfit::SharedPlanePtr sensorPlaneSptr(sensorPlane);
257
258 // do extrapolation
259 gfTrackState.extrapolateToPlane(sensorPlaneSptr);
260 } catch (genfit::Exception& gfException) {
261 B2WARNING("Fitting failed: " << gfException.getExcString());
262 isgood = false;
263 return intersec;
264 } catch (...) {
265 B2WARNING("Fitting failed: for some reason");
266 isgood = false;
267 return intersec;
268 }
269
270 //local position
271 intersec = pxdSensorInfo.pointToLocal(ROOT::Math::XYZVector(gfTrackState.getPos()), m_useAlignment);
272
273 //try to get the momentum
274 B2DEBUG(1, "Fitted momentum on the plane p = " << gfTrackState.getMom().Mag());
275
276 // no tolerance currently! Maybe one should be added!
277 double tolerance = 0.0;
278 bool inside = pxdSensorInfo.inside(intersec.X(), intersec.Y(), tolerance, tolerance);
279
280 // get intersection point in local coordinates with covariance matrix
281 TMatrixDSym covMatrix = gfTrackState.getCov(); // 5D with elements q/p,u',v',u,v in plane system
282
283 // get ROI by covariance matrix and local intersection point
284 du = std::sqrt(covMatrix(3, 3));
285 dv = std::sqrt(covMatrix(4, 4));
286
287 if (inside) isgood = true;
288
289 return intersec;
290}
291
292
294{
295 // Create a separate histogram directory and cd into it.
296 TDirectory* oldDir = gDirectory;
297 if (m_histogramDirectoryName != "") {
298 oldDir->mkdir(m_histogramDirectoryName.c_str());
299 oldDir->cd(m_histogramDirectoryName.c_str());
300 }
301
302 std::vector<VxdID> sensors = m_vxdGeometry.getListOfSensors();
303 for (VxdID& avxdid : sensors) {
304 VXD::SensorInfoBase info = m_vxdGeometry.getSensorInfo(avxdid);
305 if (info.getType() != VXD::SensorInfoBase::PXD) continue;
306 //Only interested in PXD sensors
307
308 TString buff = (std::string)avxdid;
309 buff.ReplaceAll(".", "_");
310
311 int nu = info.getUCells();
312 int nv = info.getVCells();
313
314 //nu + 1,nv + 1 Bin numbers when using one bin per pixel
315 m_h_track_hits[avxdid] = new TH2F("track_hits_" + buff, "tracks through sensor " + buff,
316 m_u_bins, -0.5, nu - 0.5, m_v_bins, -0.5, nv - 0.5);
317 m_h_matched_cluster[avxdid] = new TH2F("matched_cluster_" + buff, "clusters matched to track intersections " + buff,
318 m_u_bins, -0.5, nu - 0.5, m_v_bins, -0.5, nv - 0.5);
319
320 if (m_verboseHistos) {
321 m_h_p[avxdid] = new TH1F("p_" + buff, "p " + buff, 100, 0, 10);
322 m_h_pt[avxdid] = new TH1F("pt_" + buff, "pt " + buff, 100, 0, 10);
323 m_h_su[avxdid] = new TH1F("su_" + buff, "su " + buff, 1000, 0, 1);
324 m_h_sv[avxdid] = new TH1F("sv_" + buff, "sv " + buff, 1000, 0, 1);
325 m_h_p2[avxdid] = new TH1F("p2_" + buff, "p2 " + buff, 100, 0, 10);
326 m_h_pt2[avxdid] = new TH1F("pt2_" + buff, "pt2 " + buff, 100, 0, 10);
327 m_h_su2[avxdid] = new TH1F("su2_" + buff, "su2 " + buff, 1000, 0, 1);
328 m_h_sv2[avxdid] = new TH1F("sv2_" + buff, "sv2 " + buff, 1000, 0, 1);
329 }
330 }
331 // cd back to root directory
332 oldDir->cd();
333}
334
335
336int
337PXDDQMEfficiencySelftrackModule::findClosestCluster(const VxdID& avxdid, ROOT::Math::XYZVector intersection)
338{
339 int closest = -1;
340 double mindist = 999999999999; //definitely outside of the sensor
341
342 VXD::SensorInfoBase info = m_vxdGeometry.getSensorInfo(avxdid);
343
344 //loop the clusters
345 for (int iclus = 0; iclus < m_pxdclusters.getEntries(); iclus++) {
346 //Do not consider as different if only segment differs!
347 //As of this writing segment is never filled for clusters, but just to be sure
348 VxdID clusterID = m_pxdclusters[iclus]->getSensorID();
349 if (avxdid.getLayerNumber() != clusterID.getLayerNumber() ||
350 avxdid.getLadderNumber() != clusterID.getLadderNumber() ||
351 avxdid.getSensorNumber() != clusterID.getSensorNumber()) {
352 continue;
353 }
354 //only cluster on the correct sensor and direction should survive
355
356 double u = m_pxdclusters[iclus]->getU();
357 double v = m_pxdclusters[iclus]->getV();
358 ROOT::Math::XYZVector current(u, v, 0);
359
360 //2D dist squared
361 double dist = (intersection - current).R();
362 if (dist < mindist) {
363 closest = iclus;
364 mindist = dist;
365 }
366 }
367
368 return closest;
369
370}
371
372bool PXDDQMEfficiencySelftrackModule::isCloseToBorder(int u, int v, int checkDistance)
373{
374
375 if (u - checkDistance < 0 || u + checkDistance >= 250 ||
376 v - checkDistance < 0 || v + checkDistance >= 768) {
377 return true;
378 }
379 return false;
380}
381
382bool PXDDQMEfficiencySelftrackModule::isDeadPixelClose(int u, int v, int checkDistance, const VxdID& moduleID)
383{
384
385 //Iterate over square around the intersection to see if any close pixel is dead
386 for (int u_iter = u - checkDistance; u_iter <= u + checkDistance ; ++u_iter) {
387 for (int v_iter = v - checkDistance; v_iter <= v + checkDistance ; ++v_iter) {
388 if (PXD::PXDPixelMasker::getInstance().pixelDead(moduleID, u_iter, v_iter)
389 || !PXD::PXDPixelMasker::getInstance().pixelOK(moduleID, u_iter, v_iter)) {
390 return true;
391 }
392 }
393 }
394 return false;
395}
double R
typedef autogenerated by FFTW
static const ChargedStable pion
charged pion particle
Definition Const.h:661
HistoModule()
Constructor.
Definition HistoModule.h:32
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
@ 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
std::map< VxdID, TH1F * > m_h_su
histograms of su
void initialize() override final
initializes the need store arrays, trees and histograms
std::map< VxdID, TH2F * > m_h_track_hits
histograms of track hits
double m_d0Cut
cut abs(d0) in cm (large positive value eg 9999 disables)
bool isCloseToBorder(int u, int v, int checkDistance)
is it close to the border
std::map< VxdID, TH1F * > m_h_su2
histrograms of su2
double m_vFactor
factor for track-error on distcut comparison
unsigned int m_minSVDHits
Required hits in SVD strips for tracks.
StoreArray< PXDCluster > m_pxdclusters
store array of pxd clusters
bool m_verboseHistos
add some verbose histograms for cuts
std::map< VxdID, TH1F * > m_h_sv2
histrograms of sv2
int findClosestCluster(const VxdID &vxdid, ROOT::Math::XYZVector intersection)
find the closest cluster
PXDDQMEfficiencySelftrackModule()
Constructor: Sets the description, the properties and the parameters of the module.
void defineHisto() override final
actually defines the trees and histograms
ROOT::Math::XYZVector getTrackInterSec(const VXD::SensorInfoBase &pxdSensorInfo, const RecoTrack &aTrack, bool &isgood, double &du, double &dv)
helper functions to do some of the calculations
bool isDeadPixelClose(int u, int v, int checkDistance, const VxdID &moduleID)
is a dead pixel close
double m_z0maxCut
cut z0 maximum in cm (large positive value eg 9999 disables)
std::map< VxdID, TH1F * > m_h_pt2
histograms of pt2
bool m_requireROIs
Require tracks going through ROIs.
void event() override final
main function which fills trees and histograms
std::map< VxdID, TH1F * > m_h_sv
histograms of sv
std::string m_pxdClustersName
name of the store array of pxd clusters
std::string m_histogramDirectoryName
Where to save the histograms too.
double m_z0minCut
cut z0 minimum in cm (large negative value eg -9999 disables)
int m_maskedDistance
Distance inside which no dead pixel or module border is allowed.
std::string m_ROIsName
name of the store array of ROIs
std::map< VxdID, TH2F * > m_h_matched_cluster
histograms of matched clusters
std::map< VxdID, TH1F * > m_h_p2
histograms of p2
bool m_useAlignment
if true alignment will be used!
void beginRun() override final
begin run function which resets histograms
StoreArray< RecoTrack > m_tracks
store array of tracks
double m_uFactor
factor for track-error on distcut comparison
std::map< VxdID, TH1F * > m_h_p
histograms of momenta
std::map< VxdID, TH1F * > m_h_pt
histograms of transverse momenta
std::string m_tracksName
name of the store array of tracks
StoreArray< ROIid > m_ROIs
store array of ROIs
static PXDPixelMasker & getInstance()
Main (and only) way to access the PXDPixelMasker.
This is the Reconstruction Event-Data Model Track.
Definition RecoTrack.h:79
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
Class that bundles various TrackFitResults.
Definition Track.h:25
const TrackFitResult * getTrackFitResultWithClosestMass(const Const::ChargedStable &requestedType) const
Return the track fit for the fit hypothesis with the closest mass.
Definition Track.cc:104
Base class to provide Sensor Information for PXD and SVD.
ROOT::Math::XYZVector pointToLocal(const ROOT::Math::XYZVector &global, bool reco=false) const
Convert a point from global to local coordinates.
ROOT::Math::XYZVector pointToGlobal(const ROOT::Math::XYZVector &local, bool reco=false) const
Convert a point from local to global coordinates.
bool inside(double u, double v, double uTolerance=DBL_EPSILON, double vTolerance=DBL_EPSILON) const
Check whether a given point is inside the active area.
ROOT::Math::XYZVector vectorToGlobal(const ROOT::Math::XYZVector &local, bool reco=false) const
Convert a vector from local to global coordinates.
Class to uniquely identify a any structure of the PXD and SVD.
Definition VxdID.h:32
baseType getSensorNumber() const
Get the sensor id.
Definition VxdID.h:99
baseType getLadderNumber() const
Get the ladder id.
Definition VxdID.h:97
baseType getLayerNumber() const
Get the layer id.
Definition VxdID.h:95
static constexpr auto XYZToTVector
Helper function to convert XYZVector to TVector3.
Definition VectorUtil.h:24
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
Namespace to provide code needed by both Vertex Detectors, PXD and SVD, and also testbeam telescopes.
Definition GeoCache.h:33
Abstract base class for different kinds of events.