Belle II Software development
PXDDQMEfficiencyModule.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/PXDDQMEfficiencyModule.h>
10#include <pxd/dataobjects/PXDCluster.h>
11#include <tracking/dataobjects/ROIid.h>
12#include <tracking/dataobjects/RecoTrack.h>
13#include <tracking/dataobjects/PXDIntercept.h>
14
15#include <pxd/reconstruction/PXDPixelMasker.h>
16#include <mdst/dataobjects/Track.h>
17#include <framework/gearbox/Const.h>
18
19#include <Math/Vector3D.h>
20#include "TDirectory.h"
21#include "TMatrixDSym.h"
22using namespace Belle2;
23
24//-----------------------------------------------------------------
25// Register the Module
26//-----------------------------------------------------------------
27REG_MODULE(PXDDQMEfficiency);
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("PXDInterceptListName", m_PXDInterceptListName, "name of the list of interceptions", std::string(""));
46 addParam("histogramDirectoryName", m_histogramDirectoryName, "Name of the directory where histograms will be placed",
47 std::string("PXDEFF"));
48 addParam("binsU", m_u_bins, "histogram bins in u direction", int(16));
49 addParam("binsV", m_v_bins, "histogram bins in v direction", int(48));
50 addParam("distCut", m_distcut, "max distance in [cm] for cluster to be counted to a track", double(0.0500));
51 addParam("pCut", m_pcut, "Set a cut on the track p-value", double(1e-20));
52 addParam("requireROIs", m_requireROIs, "require tracks to lie inside a ROI", bool(false));
53 addParam("useAlignment", m_useAlignment, "if true the alignment will be used", bool(true));
54 addParam("maskDeadPixels", m_maskDeadPixels, "Do not consider tracks going through known dead or hot pixels for the efficiency",
55 bool(false));
56 addParam("minSVDHits", m_minSVDHits, "Number of SVD hits required in a track to be considered", 0u);
57 addParam("momCut", m_momCut, "Set a cut on the track momentum in GeV/c, 0 disables", double(0));
58 addParam("pTCut", m_pTCut, "Set a cut on the track pT in GeV/c, 0 disables", double(1));
59 addParam("cutBorders", m_cutBorders, "Do not use tracks near the borders of the sensor", bool(true));
60 addParam("maskedDistance", m_maskedDistance, "Distance inside which no masked pixel or sensor border is allowed", int(10));
61 addParam("trackUFactorDistCut", m_uFactor, "Set a cut on u error of track (factor*err<dist), 0 disables", double(2.0));
62 addParam("trackVFactorDistCut", m_vFactor, "Set a cut on v error of track (factor*err<dist), 0 disables", double(2.0));
63 addParam("z0minCut", m_z0minCut, "Set a cut z0 minimum in cm (large negative value eg -9999 disables)", double(-1));
64 addParam("z0maxCut", m_z0maxCut, "Set a cut z0 maximum in cm (large positive value eg 9999 disables)", double(1));
65 addParam("d0Cut", m_d0Cut, "Set a cut abs(d0) in cm (large positive value eg 9999 disables)", double(0.5));
66 addParam("verboseHistos", m_verboseHistos, "Add more verbose histograms for cuts (not for ereoc)", bool(false));
67}
68
69
71{
72 //calls the define histogram function
73 REG_HISTOGRAM;
74
75 //register the required arrays
76 //Register as optional so validation for cases where they are not available still succeeds, but module will not do any meaningful work without them
78 m_tracks.isOptional(m_tracksName);
79 m_ROIs.isOptional(m_ROIsName);
81}
82
84{
85 if (m_h_combined) m_h_combined->Reset();
86 for (auto& h : m_h_track_hits) if (h.second) h.second->Reset();
87 for (auto& h : m_h_matched_cluster) if (h.second) h.second->Reset();
88 for (auto& h : m_h_p) if (h.second) h.second->Reset();
89 for (auto& h : m_h_pt) if (h.second) h.second->Reset();
90 for (auto& h : m_h_su) if (h.second) h.second->Reset();
91 for (auto& h : m_h_sv) if (h.second) h.second->Reset();
92 for (auto& h : m_h_p2) if (h.second) h.second->Reset();
93 for (auto& h : m_h_pt2) if (h.second) h.second->Reset();
94 for (auto& h : m_h_su2) if (h.second) h.second->Reset();
95 for (auto& h : m_h_sv2) if (h.second) h.second->Reset();
96}
97
99{
100 if (!m_pxdclusters.isValid()) {
101 B2INFO("PXDClusters array is missing, no efficiencies");
102 return;
103 }
104 if (!m_tracks.isValid()) {
105 B2INFO("RecoTrack array is missing, no efficiencies");
106 return;
107 }
108 if (!m_ROIs.isValid() && m_requireROIs) {
109 B2INFO("ROI array is missing but required hits in ROIs, aborting");
110 return;
111 }
112 if (!m_intercepts.isValid()) {
113 B2INFO("Intercept array is missing, no efficiencies");
114 return;
115 }
116
117
118 for (auto& a_track : m_tracks) {
119
120 //If fit failed assume position pointed to is useless anyway
121 if (!a_track.wasFitSuccessful()) continue;
122
123 if (a_track.getNumberOfSVDHits() < m_minSVDHits) continue;
124
125 RelationVector<PXDIntercept> interceptList = a_track.getRelationsTo<PXDIntercept>(m_PXDInterceptListName);
126 if (!interceptList.size()) continue;
127
128 const genfit::FitStatus* fitstatus = a_track.getTrackFitStatus();
129 if (fitstatus->getPVal() < m_pcut) continue;
130
131 genfit::MeasuredStateOnPlane trackstate;
132 trackstate = a_track.getMeasuredStateOnPlaneFromFirstHit();
133 if (trackstate.getMom().Mag() < m_momCut) continue;
134 if (trackstate.getMom().Pt() < m_pTCut) continue;
135
136 auto ptr = a_track.getRelated<Track>("Tracks");
137
138 if (!ptr) {
139 B2ERROR("expect a track for fitted recotracks");
140 continue;
141 }
143 if (!ptr2) {
144 B2ERROR("expect a track fit result for mass");
145 continue;
146 }
147
148 // Vertex cut
149 if (ptr2->getZ0() < m_z0minCut || ptr2->getZ0() > m_z0maxCut || fabs(ptr2->getD0()) > m_d0Cut) continue;
150
151 // Loop over all intercepts to sensors
152 for (auto intercept : interceptList) {
153 auto const aVxdID = intercept.getSensorID();
154 auto& info = m_vxdGeometry.getSensorInfo(aVxdID);
155 //Search for intersections of the track with all PXD layers
156 //Traditional (aka the person before did it like this) method
157 //If there is a way to find out sensors crossed by a track directly, that would most likely be faster
158
159 //true = track intersects current sensor
160 double sigu(-9999);
161 double sigv(-9999);
162 {
163 if (m_verboseHistos) {
164 if (m_h_p[aVxdID]) m_h_p[aVxdID]->Fill(trackstate.getMom().Mag());
165 if (m_h_pt[aVxdID]) m_h_pt[aVxdID]->Fill(trackstate.getMom().Pt());
166 if (m_h_su[aVxdID]) m_h_su[aVxdID]->Fill(sigu);
167 if (m_h_sv[aVxdID]) m_h_sv[aVxdID]->Fill(sigv);
168 }
169 if (m_uFactor * sigu > m_distcut) continue; // Error ufak*SigmaU > cut
170 if (m_vFactor * sigv > m_distcut) continue; // Error vfak*SigmaV > cut
171
172 double u_fit = intercept.getCoorU();
173 double v_fit = intercept.getCoorV();
174
175 // cell ID is not limited to sensor, can be <0 and >250,768
176 int ucell_fit = info.getUCellID(u_fit);
177 int vcell_fit = info.getVCellID(v_fit);
178
179 if (m_cutBorders && isCloseToBorder(ucell_fit, vcell_fit, m_maskedDistance)) {
180 continue;
181 }
182
183 if (m_maskDeadPixels && isDeadPixelClose(ucell_fit, vcell_fit, m_maskedDistance, aVxdID)) {
184 continue;
185 }
186
187 if (m_requireROIs) {
188 //Check if the intersection is inside a ROI
189 //If not, even if measured the cluster was thrown away->Not PXD's fault
190 bool fitInsideROI = false;
191 for (auto& roit : m_ROIs) {
192 if (aVxdID != (roit.getSensorID()).getID()) {
193 continue; //ROI on other sensor
194 }
195
196 if (ucell_fit < roit.getMaxUid()
197 && ucell_fit > roit.getMinUid()
198 && vcell_fit < roit.getMaxVid()
199 && vcell_fit > roit.getMinVid()) {
200 fitInsideROI = true;
201 }
202 }
203 if (!fitInsideROI) {
204 continue;//Hit wouldn't have been recorded
205 }
206 }
207
208 //This track should be on the sensor
209 m_h_track_hits[aVxdID]->Fill(ucell_fit, vcell_fit);
210 auto sensor_index = revLUT[aVxdID];
211 m_h_combined->Fill(sensor_index * 2);
212
213 //Now check if the sensor measured a hit here
214
215 int bestcluster = findClosestCluster(aVxdID, ROOT::Math::XYZVector(u_fit, v_fit, 0));
216 if (bestcluster >= 0) {
217 double u_clus = m_pxdclusters[bestcluster]->getU();
218 double v_clus = m_pxdclusters[bestcluster]->getV();
219
220 //is the closest cluster close enough to the track to count as measured?
221 ROOT::Math::XYZVector dist_clus(u_fit - u_clus, v_fit - v_clus, 0);
222 if (dist_clus.R() <= m_distcut) {
223 m_h_matched_cluster[aVxdID]->Fill(ucell_fit, vcell_fit);
224 m_h_combined->Fill(sensor_index * 2 + 1);
225 if (m_verboseHistos) {
226 if (m_h_p2[aVxdID]) m_h_p2[aVxdID]->Fill(trackstate.getMom().Mag());
227 if (m_h_pt2[aVxdID]) m_h_pt2[aVxdID]->Fill(trackstate.getMom().Pt());
228 if (m_h_su2[aVxdID]) m_h_su2[aVxdID]->Fill(sigu);
229 if (m_h_sv2[aVxdID]) m_h_sv2[aVxdID]->Fill(sigv);
230 }
231 }
232 }
233 }
234 }
235 }
236}
237
239{
240 // Create a separate histogram directory and cd into it.
241 TDirectory* oldDir = gDirectory;
242 if (m_histogramDirectoryName != "") {
243 oldDir->mkdir(m_histogramDirectoryName.c_str());
244 oldDir->cd(m_histogramDirectoryName.c_str());
245 }
246
247 std::vector<VxdID> sensors = m_vxdGeometry.getListOfSensors();
248 std::sort(sensors.begin(), sensors.end()); // make sure it is our natural order
249 int sensor_index = 0;
250 for (VxdID& avxdid : sensors) {
251 VXD::SensorInfoBase info = m_vxdGeometry.getSensorInfo(avxdid);
252 if (info.getType() != VXD::SensorInfoBase::PXD) continue;
253 //Only interested in PXD sensors
254 revLUT[avxdid] = sensor_index++;
255
256 TString buff = (std::string)avxdid;
257 buff.ReplaceAll(".", "_");
258
259 int nu = info.getUCells();
260 int nv = info.getVCells();
261
262 //nu + 1,nv + 1 Bin numbers when using one bin per pixel
263 m_h_track_hits[avxdid] = new TH2F("track_hits_" + buff, "tracks through sensor " + buff,
264 m_u_bins, -0.5, nu - 0.5, m_v_bins, -0.5, nv - 0.5);
265 m_h_matched_cluster[avxdid] = new TH2F("matched_cluster_" + buff, "clusters matched to track intersections " + buff,
266 m_u_bins, -0.5, nu - 0.5, m_v_bins, -0.5, nv - 0.5);
267
268 if (m_verboseHistos) {
269 m_h_p[avxdid] = new TH1F("p_" + buff, "p " + buff, 100, 0, 10);
270 m_h_pt[avxdid] = new TH1F("pt_" + buff, "pt " + buff, 100, 0, 10);
271 m_h_su[avxdid] = new TH1F("su_" + buff, "su " + buff, 1000, 0, 1);
272 m_h_sv[avxdid] = new TH1F("sv_" + buff, "sv " + buff, 1000, 0, 1);
273 m_h_p2[avxdid] = new TH1F("p2_" + buff, "p2 " + buff, 100, 0, 10);
274 m_h_pt2[avxdid] = new TH1F("pt2_" + buff, "pt2 " + buff, 100, 0, 10);
275 m_h_su2[avxdid] = new TH1F("su2_" + buff, "su2 " + buff, 1000, 0, 1);
276 m_h_sv2[avxdid] = new TH1F("sv2_" + buff, "sv2 " + buff, 1000, 0, 1);
277 }
278 }
279 m_h_combined = new TH1D("PXD_Eff_combined", "Efficiency combined", revLUT.size() * 2, 0, revLUT.size() * 2);
280 // cd back to root directory
281 oldDir->cd();
282}
283
284
285int
286PXDDQMEfficiencyModule::findClosestCluster(const VxdID& avxdid, ROOT::Math::XYZVector intersection)
287{
288 int closest = -1;
289 double mindist = 999999999999; //definitely outside of the sensor
290
291 VXD::SensorInfoBase info = m_vxdGeometry.getSensorInfo(avxdid);
292
293 //loop the clusters
294 for (int iclus = 0; iclus < m_pxdclusters.getEntries(); iclus++) {
295 //Do not consider as different if only segment differs!
296 //As of this writing segment is never filled for clusters, but just to be sure
297 VxdID clusterID = m_pxdclusters[iclus]->getSensorID();
298 if (avxdid.getLayerNumber() != clusterID.getLayerNumber() ||
299 avxdid.getLadderNumber() != clusterID.getLadderNumber() ||
300 avxdid.getSensorNumber() != clusterID.getSensorNumber()) {
301 continue;
302 }
303 //only cluster on the correct sensor and direction should survive
304
305 double u = m_pxdclusters[iclus]->getU();
306 double v = m_pxdclusters[iclus]->getV();
307 ROOT::Math::XYZVector current(u, v, 0);
308
309 //2D dist squared
310 double dist = (intersection - current).R();
311 if (dist < mindist) {
312 closest = iclus;
313 mindist = dist;
314 }
315 }
316
317 return closest;
318
319}
320
321bool PXDDQMEfficiencyModule::isCloseToBorder(int u, int v, int checkDistance)
322{
323
324 if (u - checkDistance < 0 || u + checkDistance >= 250 ||
325 v - checkDistance < 0 || v + checkDistance >= 768) {
326 return true;
327 }
328 return false;
329}
330
331bool PXDDQMEfficiencyModule::isDeadPixelClose(int u, int v, int checkDistance, const VxdID& moduleID)
332{
333
334 //Iterate over square around the intersection to see if any close pixel is dead
335 for (int u_iter = u - checkDistance; u_iter <= u + checkDistance ; ++u_iter) {
336 for (int v_iter = v - checkDistance; v_iter <= v + checkDistance ; ++v_iter) {
337 if (PXD::PXDPixelMasker::getInstance().pixelDead(moduleID, u_iter, v_iter)
338 || !PXD::PXDPixelMasker::getInstance().pixelOK(moduleID, u_iter, v_iter)) {
339 return true;
340 }
341 }
342 }
343 return false;
344}
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
PXDDQMEfficiencyModule()
Constructor: Sets the description, the properties and the parameters of the module.
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
double m_distcut
distance cut in cm!
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::string m_PXDInterceptListName
intercept list name
double m_momCut
Cut on fitted track momentum.
std::map< VxdID, int > revLUT
reverse lookup sensor id -> index in histogram
std::map< VxdID, TH1F * > m_h_sv2
histrograms of sv2
int findClosestCluster(const VxdID &vxdid, ROOT::Math::XYZVector intersection)
find the closest cluster
void defineHisto() override final
actually defines the trees and histograms
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)
double m_pcut
pValue-Cut for tracks
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.
TH1D * m_h_combined
combined histograms to workaround dqm glitch
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.
double m_pTCut
Cut on fitted track pT.
StoreArray< PXDIntercept > m_intercepts
store array of PXD Intercepts
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
VXD::GeoCache & m_vxdGeometry
the geometry
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
PXDIntercept stores the U,V coordinates and uncertainties of the intersection of a track with an PXD ...
static PXDPixelMasker & getInstance()
Main (and only) way to access the PXDPixelMasker.
Class for type safe access to objects that are referred to in relations.
size_t size() const
Get number of relations.
T * getRelated(const std::string &name="", const std::string &namedRelation="") const
Get the object to or from which this object has a relation.
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.
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
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.