8#include <tracking/vxdHoughTracking/findlets/ROIFinder.h>
9#include <tracking/vxdHoughTracking/findlets/RawTrackCandCleaner.icc.h>
10#include <tracking/trackFindingCDC/utilities/StringManipulation.h>
11#include <tracking/trackFindingVXD/trackQualityEstimators/QualityEstimatorCircleFit.h>
12#include <tracking/trackFindingVXD/trackQualityEstimators/QualityEstimatorRiemannHelixFit.h>
13#include <tracking/trackFindingVXD/trackQualityEstimators/QualityEstimatorTripletFit.h>
14#include <framework/logging/Logger.h>
15#include <framework/core/ModuleParamList.h>
16#include <framework/geometry/BFieldManager.h>
17#include <framework/geometry/B2Vector3.h>
18#include <framework/database/DBObjPtr.h>
19#include <mdst/dbobjects/BeamSpot.h>
20#include <pxd/geometry/SensorInfo.h>
21#include <vxd/geometry/GeoCache.h>
24using namespace TrackFindingCDC;
25using namespace vxdHoughTracking;
38 "Calculate PXDIntercepts and ROIs in this findlet based on a simple circle extrapolation (r-phi) and straight line extrapolation (z, theta)?",
42 "Name of the PXDIntercepts StoreArray produced by SVDHoughTracking using a simple circle extrapolation in r-phi and a straight line extrapolation in theta.",
46 "Name of the ROIs StoreArray produced by SVDHoughTracking using a simple circle extrapolation in r-phi and a straight line extrapolation in theta.",
50 "Allowed tolerance in phi (in radians) (for ROI calculation in u direction). If the intercept is within this range of the active region, an intercept is created.",
53 "Allowed tolerance in z (in cm) (for ROI calculation in v direction). If the intercept is within this range of the active region of a sensor, an intercept is created.",
57 "Correct charge-dependent bias of the extrapolated hit on L1 with radiusCorrectionFactor * trackCharge / trackRadius.",
60 "Correct charge-dependent bias of the extrapolated hit on L2 with radiusCorrectionFactor * trackCharge / trackRadius.",
63 "Correct sin(trackPhi) dependent bias with sinPhiCorrectionFactor * sin(trackPhi).",
66 "Correct cos(trackPhi) dependent bias with cosPhiCorrectionFactor * cos(trackPhi).",
70 "Correction factor for the extrapolated z position.",
83 "Multiplier term for ROI size estimation for L1 u direction. Usage: multiplierUL1 * 1/R + minimumROISizeUL1, with R in cm.",
86 "Multiplier term for ROI size estimation for L2 u direction. Usage: multiplierUL2 * 1/R + minimumROISizeUL2, with R in cm.",
89 "Multiplier term for ROI size estimation for L1 v direction. Usage: (1 + abs(tan(lambda)) * multiplierVL1) + minimumROISizeVL1.",
92 "Multiplier term for ROI size estimation for L2 v direction. Usage: (1 + abs(tan(lambda)) * multiplierVL2) + minimumROISizeVL2.",
102 "Refit the track with trackQualityEstimationMethod?",
m_refit);
104 "Add virtual IP to fit the tracks again for ROI creation?",
m_addVirtualIP);
106 "Identifier which fit method to use. Valid identifiers are: [circleFit, tripletFit, helixFit]",
129 B2ERROR(
"Please check the ROI size parameters. "
130 "None of minimumROISizeUL1, minimumROISizeUL2, minimumROISizeVL1, minimumROISizeVL2, "
131 "multiplierUL1, multiplierUL2, multiplierVL1, multiplierVL2, overlapU, overlapV must be < 0!");
135 m_estimator = std::make_unique<QualityEstimatorTripletFit>();
137 m_estimator = std::make_unique<QualityEstimatorCircleFit>();
139 m_estimator = std::make_unique<QualityEstimatorRiemannHelixFit>();
160 const TMatrixDSym posErr = (*beamSpotDB).getIPPositionCovMatrix();
186 for (
auto& track : finalTracks) {
194 auto sortedHits = track.getSortedHits();
199 sortedHits.push_back(&virtualIPSpacePoint);
202 std::sort(sortedHits.begin(), sortedHits.end(),
203 [](
const SpacePoint * a,
const SpacePoint * b) {return a->getPosition().Perp() < b->getPosition().Perp(); });
205 refittedTrackEstimate =
m_estimator->estimateQualityAndProperties(sortedHits);
208 std::vector<PXDIntercept> thisTracksIntercepts;
209 thisTracksIntercepts.reserve(8);
211 B2Vector3D momentumEstimate = track.getMomSeed();
213 momentumEstimate = *refittedTrackEstimate.
p;
216 const double trackPhi = momentumEstimate.
Phi();
217 const double trackTheta = momentumEstimate.
Theta();
218 const double tanTrackLambda =
tan(M_PI_2 - trackTheta);
219 const double trackRadius = momentumEstimate.
Perp() / (0.00299792458 *
m_bFieldZ) ;
220 const double trackCharge = track.getChargeSeed();
223 for (
int layer = 1; layer <= 2; layer++) {
227 for (
int ladder = 1; ladder <= (layer == 1 ? 8 : 12); ladder++) {
233 const double sensorLocalX = sensorPerpRadius - rotatedBeamSpotX;
235 double phiDiff = trackPhi - sensorPhi;
236 if (phiDiff < -M_PI) {
237 phiDiff += 2. * M_PI;
239 if (phiDiff > M_PI) {
240 phiDiff -= 2. * M_PI;
246 if (fabs(phiDiff) > 0.25 * M_PI and fabs(phiDiff) < 0.75 * M_PI) {
249 if (fabs(phiDiff) > 0.75 * M_PI and fabs(tanTrackLambda) > 0.2 ) {
257 double relTrackCenterPhi = 0;
258 if (trackCharge < 0) {
259 relTrackCenterPhi = trackPhi - M_PI_2 - sensorPhi;
260 }
else if (trackCharge > 0) {
261 relTrackCenterPhi = trackPhi + M_PI_2 - sensorPhi;
263 const double xCenter = trackRadius * cos(relTrackCenterPhi);
264 const double yCenter = trackRadius * sin(relTrackCenterPhi);
267 if ((trackRadius * trackRadius - (sensorLocalX - xCenter) * (sensorLocalX - xCenter)) < 0) {
270 const double ytmp =
sqrt(trackRadius * trackRadius - (sensorLocalX - xCenter) * (sensorLocalX - xCenter));
271 const double yplus = yCenter + ytmp + rotatedBeamSpotY;
272 const double yminus = yCenter - ytmp + rotatedBeamSpotY;
279 const double localUPositionPlus = yplus -
c_shiftY + correctionterm;
280 const double localUPositionMinus = yminus -
c_shiftY + correctionterm;
297 for (
int sensor = 1; sensor <= 2; sensor++) {
301 double localVPosition = z - shiftZ;
302 double localVPositionPlus = zPlus - shiftZ;
303 double localVPositionMinus = zMinus - shiftZ;
308 const VxdID sensorID =
VxdID(layer, ladder, sensor);
314 intercept.
setCoorU(localUPositionPlus);
315 intercept.
setCoorV(localVPositionPlus);
318 thisTracksIntercepts.push_back(intercept);
324 intercept.
setCoorU(localUPositionMinus);
325 intercept.
setCoorV(localVPositionMinus);
328 thisTracksIntercepts.push_back(intercept);
335 const double omega = 1. / trackRadius;
340 for (
auto& intercept : thisTracksIntercepts) {
341 const VxdID& interceptSensorID = intercept.getSensorID();
343 double uCoordinate = intercept.getCoorU();
344 double vCoordinate = intercept.getCoorV();
349 int interceptUCell = currentSensor->
getUCellID(uCoordinate, vCoordinate,
false);
350 int interceptVCell = currentSensor->
getVCellID(vCoordinate,
false);
351 int nUCells = currentSensor->
getUCells();
352 int nVCells = currentSensor->
getVCells();
354 unsigned short uSize = (interceptSensorID.
getLayerNumber() == 1 ? uSizeL1 : uSizeL2);
355 unsigned short vSize = (interceptSensorID.
getLayerNumber() == 1 ? vSizeL1 : vSizeL2);
364 short uCellDownLeft = interceptUCell - uSize / 2;
365 short vCellDownLeft = interceptVCell - vSize / 2;
368 short uCellUpRight = interceptUCell + uSize / 2;
369 short vCellUpRight = interceptVCell + vSize / 2;
371 if (uCellDownLeft >= nUCells or vCellDownLeft >= nVCells or uCellUpRight < 0 or vCellUpRight < 0) {
375 if (uCellDownLeft < 0) {
378 if (vCellDownLeft < 0) {
383 if (uCellUpRight >= nUCells) {
384 uCellUpRight = nUCells - 1;
386 if (vCellUpRight >= nVCells) {
387 vCellUpRight = nVCells - 1;
390 m_storeROIs.appendNew(
ROIid(uCellDownLeft, uCellUpRight, vCellDownLeft, vCellUpRight, interceptSensorID));
DataType Phi() const
The azimuth angle.
DataType Theta() const
The polar angle.
DataType Perp() const
The transverse component (R in cylindrical coordinate system).
bool isValid() const
Check whether a valid object was obtained from the database.
Class for accessing objects in the database.
The Module parameter list class.
PXDIntercept stores the U,V coordinates and uncertainties of the intersection of a track with an PXD ...
Specific implementation of SensorInfo for PXD Sensors which provides additional pixel specific inform...
ROIid stores the U and V ids and the sensor id of the Region Of Interest.
@ c_isActive
bit 11: SPTC is active (i.e.
SpacePoint typically is build from 1 PXDCluster or 1-2 SVDClusters.
void initialize() override
void beginEvent() override
virtual void exposeParameters(ModuleParamList *moduleParamList, const std::string &prefix)
static const double T
[tesla]
void setVxdID(VxdID::baseType user_vxdID)
set the sensor ID
void setCoorU(double user_coorU)
set the U coordinate of the intercept
void setCoorV(double user_coorV)
set the V coordinate of the intercept
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.
int getVCells() const
Return number of pixel/strips in v direction.
@ VXD
Any type of VXD Sensor.
int getUCells() const
Return number of pixel/strips in u direction.
int getVCellID(double v, bool clamp=false) const
Return the corresponding pixel/strip ID of a given v coordinate.
int getUCellID(double u, double v=0, bool clamp=false) const
Return the corresponding pixel/strip ID of a given u coordinate.
Class to uniquely identify a any structure of the PXD and SVD.
baseType getLayerNumber() const
Get the layer id.
double m_radiusCorrectionFactorL2
Correction factor for radial bias for L2: factor * charge / radius.
unsigned short m_maximumROISizeV
maximum ROI size in v in pixel
double m_minimumROISizeUL1
Minimum size of ROI in u-direction on L1 in pixel.
const double c_ladderPhiL2[12]
Phi values of the ladders of L2.
B2Vector3D m_BeamSpotPosition
B2Vector3D actually containing the BeamSpot position. This will be used as the starting point of the ...
bool m_calculateROI
Calculate ROI in this findlet?
~ROIFinder()
Default destructor.
double m_multiplierVL1
Multiplier term for v-direction on L1.
double m_radiusCorrectionFactorL1
Correction factor for radial bias for L1: factor * charge / radius.
const double c_centerZShiftLayer2[2]
Shift of the center of the active area of each sensor in a ladder of layer 2 For use of mhp_z > (leng...
double m_minimumROISizeVL2
Minimum size of ROI in v-direction on L2 in pixel.
double m_multiplierUL2
Multiplier term for u-direction on L2.
const double c_centerZShiftLayer1[2]
Shift of the center of the active area of each sensor in a ladder of layer 1 For use of mhp_z > (leng...
void initialize() override
Initialize the StoreArrays.
double m_multiplierVL2
Multiplier term for v-direction on L2.
StoreArray< ROIid > m_storeROIs
ROIs StoreArray.
double m_sinPhiCorrectionFactor
Correction factor for the sin(phi) modulation.
const double c_activeSensorLength[2]
Length of the active region for L1 and L2.
const double c_sensorMaxY
Maximum y coordinate if x-axis is perpendicular to the sensor:
std::string m_ROIFitMethod
Refit with this estimator, options are circleFit, tripletFit, helixFit.
const double c_ladderPhiL1[8]
Phi values of the ladders of L1.
double m_tolerancePhi
Allowed tolerance (in radians) phi to create intercepts per sensor.
std::string m_storeROIsName
Name of the ROIs StoreArray.
std::string m_storePXDInterceptsName
Name of the PXDIntercepts StoreArray.
const double c_layerRadius[2]
Radius of the two layers.
void apply(const std::vector< SpacePointTrackCand > &finalTracks) override
Function to call all the sub-findlets.
bool m_addVirtualIP
Add a virtual IP for the refit?
StoreArray< PXDIntercept > m_storePXDIntercepts
PXDIntercepts StoreArray.
void beginRun() override
Initialize the BField.
B2Vector3D m_BeamSpotPositionError
B2Vector3D actually containing the BeamSpot position error.
bool m_refit
Refit the tracks with m_ROIFitMethod.
double m_zPositionCorrectionFactor
Correction factor for the z position.
void beginEvent() override
Clear the object pools.
double m_cosPhiCorrectionFactor
Correction factor for the cos(phi) modulation.
std::unique_ptr< QualityEstimatorBase > m_estimator
pointer to the selected QualityEstimator
double m_bFieldZ
BField in Tesla.
double m_minimumROISizeVL1
Minimum size of ROI in v-direction on L1 in pixel.
const double c_shiftY
Shift of the center position in y if the x-axis is perpendicular to the sensor:
const double c_sensorMinY
PXD is shifted to create the windmill structure.
double m_minimumROISizeUL2
Minimum size of ROI in u-direction on L2 in pixel.
unsigned short m_maximumROISizeU
maximum ROI size in u in pixel
void exposeParameters(ModuleParamList *moduleParamList, const std::string &prefix) override
Expose the parameters of the sub findlets.
ROIFinder()
Constructor for adding the subfindlets.
double m_multiplierUL1
Multiplier term in ROI size estimation For u: size = multiplier * 1/R + minimumROISize For v: size = ...
double m_toleranceZ
Allowed tolerance (in cm) in z to create intercepts per sensor.
const double c_activeSensorWidth
PXD sensors in L1 and L2 have the same size in u direction (=width):
void addParameter(const std::string &name, T ¶mVariable, const std::string &description, const T &defaultValue)
Adds a new parameter to the module list.
B2Vector3< double > B2Vector3D
typedef for common usage with double
static void getField(const double *pos, double *field)
return the magnetic field at a given position.
double tan(double a)
tan for double
double sqrt(double a)
sqrt for double
Abstract base class for different kinds of events.
Container for complete fit/estimation results.
std::optional< B2Vector3D > p
momentum vector estimate from the QE