Belle II Software development
PXDDQMEfficiencyNtupleSelftrackModule.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/PXDDQMEfficiencyNtupleSelftrackModule.h>
10
11#include <pxd/dataobjects/PXDCluster.h>
12#include <tracking/dataobjects/ROIid.h>
13#include <mdst/dataobjects/Track.h>
14#include <tracking/dataobjects/RecoTrack.h>
15
16#include <pxd/reconstruction/PXDPixelMasker.h>
17#include <framework/gearbox/Const.h>
18#include <framework/geometry/VectorUtil.h>
19
20#include "TMatrixDSym.h"
21#include <Math/Vector3D.h>
22
23using namespace Belle2;
24
25//-----------------------------------------------------------------
26// Register the Module
27//-----------------------------------------------------------------
28REG_MODULE(PXDDQMEfficiencyNtupleSelftrack);
29
30//-----------------------------------------------------------------
31// Implementation
32//-----------------------------------------------------------------
33
35 m_vxdGeometry(VXD::GeoCache::getInstance())
36{
37 // Set module properties
38 setDescription("Create basic histograms for PXD efficiency");
39
40 // setPropertyFlags(c_ParallelProcessingCertified);// for ntuple not certified!!!
41
42 // Parameter definitions
43 addParam("pxdClustersName", m_pxdClustersName, "name of StoreArray with PXD cluster", std::string(""));
44 addParam("recoTracksName", m_recoTracksName, "name of StoreArray with RecoTracks", std::string(""));
45 addParam("tracksName", m_tracksName, "name of StoreArray with Tracks", std::string(""));
46 addParam("ROIsName", m_ROIsName, "name of the list of HLT ROIs, if available in output", std::string(""));
47 addParam("useAlignment", m_useAlignment, "if true the alignment will be used", true);
48 addParam("pCut", m_pcut, "Set a cut on the track fit p-value (0=no cut)", double(0));
49 addParam("minSVDHits", m_minSVDHits, "Number of SVD hits required in a track to be considered", 0u);
50 addParam("momCut", m_momCut, "Set a cut on the track momentum in GeV/c, 0 disables", double(0));
51 addParam("pTCut", m_pTCut, "Set a cut on the track pT in GeV/c, 0 disables", double(0));
52 addParam("maskedDistance", m_maskedDistance, "Distance inside which no masked pixel or sensor border is allowed", int(10));
53}
54
55
57{
58 auto dir = gDirectory;
59 if (m_tuple) {
60 if (m_file) { // no file -> no write to file
61 m_file->cd();
62 }
63 m_tuple->Write();
64 delete m_tuple;
65 m_tuple = nullptr;
66 }
67 if (m_file) {
68 m_file->Write();
69 m_file->Close();
70 delete m_file;
71 m_file = nullptr;
72 }
73 dir->cd();
74}
75
76
78{
79 m_file = new TFile("test.root", "recreate");
80 if (m_file) m_file->cd();
81 m_tuple = new TNtuple("effcontrol", "effcontrol",
82 "vxdid:u:v:p:pt:distu:distv:sigu:sigv:dist:inroi:clborder:cldead:matched:z0:d0:svdhits:charge:phi:costheta");
83
84 //register the required arrays
85 //Register as optional so validation for cases where they are not available still succeeds, but module will not do any meaningful work without them
88 m_tracks.isOptional(m_tracksName);
89 m_ROIs.isOptional(m_ROIsName);
90}
91
92
94{
95 if (!m_pxdclusters.isValid()) {
96 B2INFO("PXDClusters array is missing, no efficiencies");
97 return;
98 }
99 if (!m_recoTracks.isValid()) {
100 B2INFO("RecoTrack array is missing, no efficiencies");
101 return;
102 }
103 if (!m_tracks.isValid()) {
104 B2INFO("Track array is missing, no efficiencies");
105 return;
106 }
107
108 for (auto& track : m_tracks) {
109 RelationVector<RecoTrack> recoTrack = track.getRelationsTo<RecoTrack>(m_recoTracksName);
110 if (!recoTrack.size()) continue;
111
112 auto a_track = recoTrack[0];
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 const TrackFitResult* ptr2 = track.getTrackFitResultWithClosestMass(Const::pion);
127 if (!ptr2) {
128 B2ERROR("expect a track fit result for mass");
129 continue;
130 }
131
132 //loop over all PXD sensors to get the intersections
133 std::vector<VxdID> sensors = m_vxdGeometry.getListOfSensors();
134 for (VxdID& aVxdID : sensors) {
135 VXD::SensorInfoBase info = m_vxdGeometry.getSensorInfo(aVxdID);
136 if (info.getType() != VXD::SensorInfoBase::PXD) continue;
137 //Search for intersections of the track with all PXD layers
138 //Traditional (aka the person before did it like this) method
139 //If there is a way to find out sensors crossed by a track directly, that would most likely be faster
140
141 bool isgood = false;
142 //true = track intersects current sensor
143 double sigu(-9999);
144 double sigv(-9999);
145 ROOT::Math::XYZVector intersec_buff = getTrackInterSec(info, *a_track, isgood, sigu, sigv);
146
147 if (!isgood) {
148 continue;//track does not go through this sensor-> nothing to measure anyway
149 } else {
150
151 double u_fit = intersec_buff.X();
152 double v_fit = intersec_buff.Y();
153
154 int ucell_fit = info.getUCellID(intersec_buff.X());
155 int vcell_fit = info.getVCellID(intersec_buff.Y());
156
157 bool closeToBoarder = false;
158 if (isCloseToBorder(ucell_fit, vcell_fit, m_maskedDistance)) {
159 closeToBoarder = true;
160 }
161
162 bool closeToDead = false;
163 if (isDeadPixelClose(ucell_fit, vcell_fit, m_maskedDistance, aVxdID)) {
164 closeToDead = true;
165 }
166
167 bool fitInsideROI = false;
168 if (m_ROIs.isValid()) {
169 //Check if the intersection is inside a ROI
170 //If not, even if measured the cluster was thrown away->Not PXD's fault
171 for (auto& roit : m_ROIs) {
172 if (aVxdID != roit.getSensorID()) {
173 continue; //ROI on other sensor
174 }
175
176 if (ucell_fit < roit.getMaxUid()
177 && ucell_fit > roit.getMinUid()
178 && vcell_fit < roit.getMaxVid()
179 && vcell_fit > roit.getMinVid()) {
180 fitInsideROI = true;
181 }
182 }
183 }
184
185 //Now check if the sensor measured a hit here
186
187 int bestcluster = findClosestCluster(aVxdID, intersec_buff);
188 double du_clus = 0;
189 double dv_clus = 0;
190 double d_clus = 0;
191 float charge = 0;
192 bool matched = false;
193 if (bestcluster >= 0) {
194 double u_clus = m_pxdclusters[bestcluster]->getU();
195 double v_clus = m_pxdclusters[bestcluster]->getV();
196
197 //is the closest cluster close enough to the track to count as measured?
198 ROOT::Math::XYZVector dist_clus(u_fit - u_clus, v_fit - v_clus, 0);
199 du_clus = u_fit - u_clus;
200 dv_clus = v_fit - v_clus;
201 d_clus = dist_clus.R();
202 charge = m_pxdclusters[bestcluster]->getCharge();
203 matched = true;
204 }
205 float fill[22] = {float((int)aVxdID), float(u_fit), float(v_fit), float(trackstate.getMom().Mag()), float(trackstate.getMom().Pt()),
206 float(du_clus), float(dv_clus), float(sigu), float(sigv), float(d_clus),
207 float(fitInsideROI), float(closeToBoarder), float(closeToDead), float(matched),
208 float(ptr2->getZ0()), float(ptr2->getD0()), float(a_track->getNumberOfSVDHits()),
209 charge, float(trackstate.getMom().Phi()), float(trackstate.getMom().CosTheta())
210 };
211 m_tuple->Fill(fill);
212 }
213 }
214 }
215}
216
217
218
219ROOT::Math::XYZVector
221 bool& isgood, double& du, double& dv)
222{
223 //will be set true if the intersect was found
224 isgood = false;
225
226 ROOT::Math::XYZVector intersec(99999999, 9999999, 0); //point outside the sensor
227
228 genfit::MeasuredStateOnPlane gfTrackState = aTrack.getMeasuredStateOnPlaneFromFirstHit();
229
230 //adopted (aka stolen) from tracking/modules/pxdClusterRescue/PXDClusterRescueROIModule
231 try {
232 // get sensor plane
233 ROOT::Math::XYZVector zeroVec(0, 0, 0);
234 ROOT::Math::XYZVector uVec(1, 0, 0);
235 ROOT::Math::XYZVector vVec(0, 1, 0);
236
237 genfit::DetPlane* sensorPlane = new genfit::DetPlane();
238 sensorPlane->setO(XYZToTVector(pxdSensorInfo.pointToGlobal(zeroVec, m_useAlignment)));
239 sensorPlane->setUV(XYZToTVector(pxdSensorInfo.vectorToGlobal(uVec, m_useAlignment)),
240 XYZToTVector(pxdSensorInfo.vectorToGlobal(vVec, m_useAlignment)));
241
242 //boost pointer (will be deleted automatically ?!?!?)
243 genfit::SharedPlanePtr sensorPlaneSptr(sensorPlane);
244
245 // do extrapolation
246 gfTrackState.extrapolateToPlane(sensorPlaneSptr);
247 } catch (genfit::Exception& gfException) {
248 B2WARNING("Fitting failed: " << gfException.getExcString());
249 isgood = false;
250 return intersec;
251 } catch (...) {
252 B2WARNING("Fitting failed: for some reason");
253 isgood = false;
254 return intersec;
255 }
256
257 //local position
258 intersec = pxdSensorInfo.pointToLocal(ROOT::Math::XYZVector(gfTrackState.getPos()), m_useAlignment);
259
260 //try to get the momentum
261 B2DEBUG(1, "Fitted momentum on the plane p = " << gfTrackState.getMom().Mag());
262
263 // no tolerance currently! Maybe one should be added!
264 double tolerance = 0.0;
265 bool inside = pxdSensorInfo.inside(intersec.X(), intersec.Y(), tolerance, tolerance);
266
267 // get intersection point in local coordinates with covariance matrix
268 TMatrixDSym covMatrix = gfTrackState.getCov(); // 5D with elements q/p,u',v',u,v in plane system
269
270 // get ROI by covariance matrix and local intersection point
271 du = std::sqrt(covMatrix(3, 3));
272 dv = std::sqrt(covMatrix(4, 4));
273
274 if (inside) isgood = true;
275
276 return intersec;
277}
278
279
280int
281PXDDQMEfficiencyNtupleSelftrackModule::findClosestCluster(const VxdID& avxdid, ROOT::Math::XYZVector intersection)
282{
283 int closest = -1;
284 double mindist = 999999999999; //definitely outside of the sensor
285
286 VXD::SensorInfoBase info = m_vxdGeometry.getSensorInfo(avxdid);
287
288 //loop the clusters
289 for (int iclus = 0; iclus < m_pxdclusters.getEntries(); iclus++) {
290 //Do not consider as different if only segment differs!
291 //As of this writing segment is never filled for clusters, but just to be sure
292 VxdID clusterID = m_pxdclusters[iclus]->getSensorID();
293 if (avxdid.getLayerNumber() != clusterID.getLayerNumber() ||
294 avxdid.getLadderNumber() != clusterID.getLadderNumber() ||
295 avxdid.getSensorNumber() != clusterID.getSensorNumber()) {
296 continue;
297 }
298 //only cluster on the correct sensor and direction should survive
299
300 double u = m_pxdclusters[iclus]->getU();
301 double v = m_pxdclusters[iclus]->getV();
302 ROOT::Math::XYZVector current(u, v, 0);
303
304 //2D dist squared
305 double dist = (intersection - current).R();
306 if (dist < mindist) {
307 closest = iclus;
308 mindist = dist;
309 }
310 }
311
312 return closest;
313
314}
315
317{
318
319 if (u - checkDistance < 0 || u + checkDistance >= 250 ||
320 v - checkDistance < 0 || v + checkDistance >= 768) {
321 return true;
322 }
323 return false;
324}
325
326bool PXDDQMEfficiencyNtupleSelftrackModule::isDeadPixelClose(int u, int v, int checkDistance, const VxdID& moduleID)
327{
328
329 //Iterate over square around the intersection to see if any close pixel is dead
330 for (int u_iter = u - checkDistance; u_iter <= u + checkDistance ; ++u_iter) {
331 for (int v_iter = v - checkDistance; v_iter <= v + checkDistance ; ++v_iter) {
332 if (PXD::PXDPixelMasker::getInstance().pixelDead(moduleID, u_iter, v_iter)
333 || !PXD::PXDPixelMasker::getInstance().pixelOK(moduleID, u_iter, v_iter)) {
334 return true;
335 }
336 }
337 }
338 return false;
339}
double R
typedef autogenerated by FFTW
static const ChargedStable pion
charged pion particle
Definition Const.h:661
void setDescription(const std::string &description)
Sets the description of the module.
Definition Module.cc:214
Module()
Constructor.
Definition Module.cc:30
void initialize() override final
initializes the need store arrays, trees and histograms
bool isCloseToBorder(int u, int v, int checkDistance)
is it close to the border
unsigned int m_minSVDHits
Required hits in SVD strips for tracks.
StoreArray< PXDCluster > m_pxdclusters
store array of pxd clusters
int findClosestCluster(const VxdID &vxdid, ROOT::Math::XYZVector intersection)
find the closest cluster
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
void terminate() override final
terminate , save tuple to file if needed
void event() override final
main function which fills trees and histograms
std::string m_pxdClustersName
name of the store array of pxd clusters
std::string m_recoTracksName
name of the store array of recotracks
int m_maskedDistance
Distance inside which no dead pixel or module border is allowed.
StoreArray< RecoTrack > m_recoTracks
store array of reco tracks
PXDDQMEfficiencyNtupleSelftrackModule()
Constructor: Sets the description, the properties and the parameters of the module.
std::string m_tracksName
name of the store array of tracks
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 for type safe access to objects that are referred to in relations.
size_t size() const
Get number of relations.
Values of the result of a track fit with a given particle hypothesis.
double getD0() const
Getter for d0.
double getZ0() const
Getter for z0.
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.