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