File indexing completed on 2023-06-16 23:22:36
0001 #include "RecoTracker/PixelLowPtUtilities/interface/ClusterShapeHitFilter.h"
0002
0003 #include "FWCore/Framework/interface/EventSetup.h"
0004 #include "FWCore/Framework/interface/ESHandle.h"
0005 #include "FWCore/ParameterSet/interface/FileInPath.h"
0006 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0007
0008 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHit.h"
0009 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2D.h"
0010
0011 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0012
0013 #include "Geometry/CommonDetUnit/interface/PixelGeomDetUnit.h"
0014 #include "Geometry/TrackerGeometryBuilder/interface/StripGeomDetUnit.h"
0015 #include "Geometry/CommonTopologies/interface/StripTopology.h"
0016 #include "Geometry/CommonDetUnit/interface/GeomDetType.h"
0017
0018 #include "Geometry/TrackerGeometryBuilder/interface/RectangularPixelTopology.h"
0019 #include "Geometry/CommonTopologies/interface/RectangularStripTopology.h"
0020
0021 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0022
0023 #include "MagneticField/Engine/interface/MagneticField.h"
0024
0025 #include "CondFormats/SiPixelObjects/interface/SiPixelLorentzAngle.h"
0026
0027 #include "CondFormats/SiStripObjects/interface/SiStripLorentzAngle.h"
0028
0029 #include <fstream>
0030 #include <cassert>
0031
0032 using namespace std;
0033
0034
0035 ClusterShapeHitFilter::ClusterShapeHitFilter(const TrackerGeometry* theTracker_,
0036 const TrackerTopology* theTkTopol_,
0037 const MagneticField* theMagneticField_,
0038 const SiPixelLorentzAngle* theSiPixelLorentzAngle_,
0039 const SiStripLorentzAngle* theSiStripLorentzAngle_,
0040 const std::string& pixelShapeFile_,
0041 const std::string& pixelShapeFileL1_)
0042 : theTracker(theTracker_),
0043 theTkTopol(theTkTopol_),
0044 theMagneticField(theMagneticField_),
0045 theSiPixelLorentzAngle(theSiPixelLorentzAngle_),
0046 theSiStripLorentzAngle(theSiStripLorentzAngle_) {
0047
0048 loadPixelLimits(pixelShapeFile_, pixelLimits);
0049 loadPixelLimits(pixelShapeFileL1_, pixelLimitsL1);
0050 fillPixelData();
0051
0052
0053 loadStripLimits();
0054 fillStripData();
0055 cutOnPixelCharge_ = cutOnStripCharge_ = false;
0056 cutOnPixelShape_ = cutOnStripShape_ = true;
0057 }
0058
0059
0060 ClusterShapeHitFilter::~ClusterShapeHitFilter() {}
0061
0062
0063 void ClusterShapeHitFilter::loadPixelLimits(std::string const& file, PixelLimits* plim) {
0064 edm::FileInPath fileInPath(file.c_str());
0065 ifstream inFile(fileInPath.fullPath().c_str());
0066
0067 while (inFile.eof() == false) {
0068 int part, dx, dy;
0069
0070 inFile >> part;
0071 inFile >> dx;
0072 inFile >> dy;
0073
0074 const PixelKeys key(part, dx, dy);
0075 auto& pl = plim[key];
0076
0077 for (int b = 0; b < 2; b++)
0078 for (int d = 0; d < 2; d++)
0079 for (int k = 0; k < 2; k++)
0080 inFile >> pl.data[b][d][k];
0081
0082 double f;
0083 int d;
0084
0085 inFile >> f;
0086 inFile >> d;
0087 inFile >> f;
0088 inFile >> d;
0089 }
0090
0091 inFile.close();
0092
0093 LogTrace("MinBiasTracking|ClusterShapeHitFilter") << " [ClusterShapeHitFilter] pixel-cluster-shape filter loaded";
0094 }
0095
0096
0097 void ClusterShapeHitFilter::loadStripLimits() {
0098
0099 edm::FileInPath fileInPath("RecoTracker/PixelLowPtUtilities/data/stripShape.par");
0100 ifstream inFile(fileInPath.fullPath().c_str());
0101
0102 while (inFile.eof() == false) {
0103 int dx;
0104 inFile >> dx;
0105
0106 StripKeys key(dx);
0107 auto& sl = stripLimits[key];
0108
0109 for (int b = 0; b < 2; b++)
0110 for (int k = 0; k < 2; k++)
0111 inFile >> sl.data[b][k];
0112 }
0113
0114 inFile.close();
0115
0116 LogTrace("MinBiasTracking|ClusterShapeHitFilter") << " [ClusterShapeHitFilter] strip-cluster-width filter loaded";
0117 }
0118
0119 void ClusterShapeHitFilter::fillPixelData() {
0120
0121 for (auto det : theTracker->detsPXB()) {
0122
0123 const PixelGeomDetUnit* pixelDet = dynamic_cast<const PixelGeomDetUnit*>(det);
0124 if (pixelDet) {
0125 PixelData& pd = pixelData[pixelDet->geographicalId()];
0126 pd.det = pixelDet;
0127 pd.part = 0;
0128 pd.layer = theTkTopol->pxbLayer(pixelDet->geographicalId());
0129 pd.cotangent = getCotangent(pixelDet);
0130 pd.drift = getDrift(pixelDet);
0131 }
0132 }
0133
0134
0135 for (auto det : theTracker->detsPXF()) {
0136
0137 const PixelGeomDetUnit* pixelDet = dynamic_cast<const PixelGeomDetUnit*>(det);
0138 assert(pixelDet);
0139 PixelData& pd = pixelData[pixelDet->geographicalId()];
0140 pd.det = pixelDet;
0141 pd.part = 1;
0142 pd.layer = 0;
0143 pd.cotangent = getCotangent(pixelDet);
0144 pd.drift = getDrift(pixelDet);
0145 }
0146 }
0147
0148 void ClusterShapeHitFilter::fillStripData() {
0149
0150 auto const& geom_ = *theTracker;
0151 auto const& dus = geom_.detUnits();
0152 auto offset = dus.size();
0153 for (unsigned int i = 1; i < 7; ++i) {
0154 if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) != dus.size() &&
0155 dus[geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i])]->type().isTrackerStrip()) {
0156 if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) < offset)
0157 offset = geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]);
0158 }
0159 }
0160
0161 for (auto i = offset; i != dus.size(); ++i) {
0162 const StripGeomDetUnit* stripdet = (const StripGeomDetUnit*)(dus[i]);
0163 assert(stripdet->index() == int(i));
0164 assert(stripdet->type().isTrackerStrip());
0165 auto const& bounds = stripdet->specificSurface().bounds();
0166 auto detid = stripdet->geographicalId();
0167 auto& p = stripData[detid];
0168 p.det = stripdet;
0169 p.topology = (StripTopology*)(&stripdet->topology());
0170 p.drift = getDrift(stripdet);
0171 p.thickness = bounds.thickness();
0172 p.nstrips = p.topology->nstrips();
0173 }
0174 }
0175
0176
0177 pair<float, float> ClusterShapeHitFilter::getCotangent(const PixelGeomDetUnit* pixelDet) const {
0178 pair<float, float> cotangent;
0179
0180 cotangent.first = pixelDet->surface().bounds().thickness() / pixelDet->specificTopology().pitch().first;
0181 cotangent.second = pixelDet->surface().bounds().thickness() / pixelDet->specificTopology().pitch().second;
0182
0183 return cotangent;
0184 }
0185
0186
0187 float ClusterShapeHitFilter::getCotangent(const StripData& sd, const LocalPoint& pos) const {
0188
0189 return sd.thickness / sd.topology->localPitch(pos);
0190 }
0191
0192
0193 pair<float, float> ClusterShapeHitFilter::getDrift(const PixelGeomDetUnit* pixelDet) const {
0194 LocalVector lBfield = (pixelDet->surface()).toLocal(theMagneticField->inTesla(pixelDet->surface().position()));
0195
0196 double theTanLorentzAnglePerTesla = theSiPixelLorentzAngle->getLorentzAngle(pixelDet->geographicalId().rawId());
0197
0198 pair<float, float> dir;
0199 dir.first = -theTanLorentzAnglePerTesla * lBfield.y();
0200 dir.second = theTanLorentzAnglePerTesla * lBfield.x();
0201
0202 return dir;
0203 }
0204
0205
0206 float ClusterShapeHitFilter::getDrift(const StripGeomDetUnit* stripDet) const {
0207 LocalVector lBfield = (stripDet->surface()).toLocal(theMagneticField->inTesla(stripDet->surface().position()));
0208
0209 double theTanLorentzAnglePerTesla = 0.0;
0210 if (theSiStripLorentzAngle) {
0211 theTanLorentzAnglePerTesla = theSiStripLorentzAngle->getLorentzAngle(stripDet->geographicalId().rawId());
0212 }
0213
0214 float dir = theTanLorentzAnglePerTesla * lBfield.y();
0215
0216 return dir;
0217 }
0218
0219
0220 bool ClusterShapeHitFilter::isNormalOriented(const GeomDetUnit* geomDet) const {
0221 if (geomDet->type().isBarrel()) {
0222 float perp0 = geomDet->toGlobal(Local3DPoint(0., 0., 0.)).perp();
0223 float perp1 = geomDet->toGlobal(Local3DPoint(0., 0., 1.)).perp();
0224 return (perp1 > perp0);
0225 } else {
0226 float rot = geomDet->toGlobal(LocalVector(0., 0., 1.)).z();
0227 float pos = geomDet->toGlobal(Local3DPoint(0., 0., 0.)).z();
0228 return (rot * pos > 0);
0229 }
0230 }
0231
0232
0233
0234
0235 bool ClusterShapeHitFilter::getSizes(const SiPixelRecHit& recHit,
0236 const LocalVector& ldir,
0237 const SiPixelClusterShapeCache& clusterShapeCache,
0238 int& part,
0239 ClusterData::ArrayType& meas,
0240 pair<float, float>& pred,
0241 PixelData const* ipd) const {
0242
0243 const PixelData& pd = getpd(recHit, ipd);
0244
0245
0246 const SiPixelClusterShapeData& data = clusterShapeCache.get(recHit.cluster(), pd.det);
0247 bool usable = (data.isStraight() && data.isComplete());
0248
0249
0250
0251 {
0252 part = pd.part;
0253
0254
0255 pred.first = ldir.x() / ldir.z();
0256 pred.second = ldir.y() / ldir.z();
0257
0258 SiPixelClusterShapeData::Range sizeRange = data.size();
0259 if (sizeRange.first->second < 0)
0260 pred.second = -pred.second;
0261
0262 meas.clear();
0263 assert(meas.capacity() >= std::distance(sizeRange.first, sizeRange.second));
0264 for (auto s = sizeRange.first; s != sizeRange.second; ++s) {
0265 meas.push_back_unchecked(*s);
0266 }
0267 if (sizeRange.first->second < 0) {
0268 for (auto& s : meas)
0269 s.second = -s.second;
0270 }
0271
0272
0273 std::pair<float, float> const& drift = pd.drift;
0274 pred.first += drift.first;
0275 pred.second += drift.second;
0276
0277
0278 std::pair<float, float> const& cotangent = pd.cotangent;
0279 pred.first *= cotangent.first;
0280 pred.second *= cotangent.second;
0281 }
0282
0283
0284 return usable;
0285 }
0286
0287 bool ClusterShapeHitFilter::isCompatible(const SiPixelRecHit& recHit,
0288 const LocalVector& ldir,
0289 const SiPixelClusterShapeCache& clusterShapeCache,
0290 PixelData const* ipd) const {
0291
0292 if (cutOnPixelCharge_ && (!checkClusterCharge(recHit.geographicalId(), *(recHit.cluster()), ldir)))
0293 return false;
0294 if (!cutOnPixelShape_)
0295 return true;
0296
0297 const PixelData& pd = getpd(recHit, ipd);
0298
0299 int part;
0300 ClusterData::ArrayType meas;
0301 pair<float, float> pred;
0302
0303 PixelLimits const* pl = pd.layer == 1 ? pixelLimitsL1 : pixelLimits;
0304 if (getSizes(recHit, ldir, clusterShapeCache, part, meas, pred, &pd)) {
0305 for (const auto& m : meas) {
0306 PixelKeys key(part, m.first, m.second);
0307 if (!key.isValid())
0308 return true;
0309 if (pl[key].isInside(pred))
0310 return true;
0311 }
0312
0313 return false;
0314 }
0315
0316 return true;
0317 }
0318
0319 bool ClusterShapeHitFilter::isCompatible(const SiPixelRecHit& recHit,
0320 const GlobalVector& gdir,
0321 const SiPixelClusterShapeCache& clusterShapeCache,
0322 PixelData const* ipd) const {
0323
0324 const PixelData& pd = getpd(recHit, ipd);
0325
0326 LocalVector ldir = pd.det->toLocal(gdir);
0327
0328 return isCompatible(recHit, ldir, clusterShapeCache, &pd);
0329 }
0330
0331
0332
0333 bool ClusterShapeHitFilter::getSizes(DetId id,
0334 const SiStripCluster& cluster,
0335 const LocalPoint& lpos,
0336 const LocalVector& ldir,
0337 int& meas,
0338 float& pred) const {
0339
0340 auto const& p = getsd(id);
0341
0342
0343 meas = cluster.amplitudes().size();
0344
0345
0346 int fs = cluster.firstStrip();
0347 int ns = p.nstrips;
0348
0349 bool usable = (fs >= 1) & (fs + meas <= ns + 1);
0350
0351
0352
0353 {
0354
0355 pred = ldir.x() / ldir.z();
0356
0357
0358 float drift = p.drift;
0359 ;
0360 pred += drift;
0361
0362
0363 pred *= getCotangent(p, lpos);
0364 }
0365
0366 return usable;
0367 }
0368
0369
0370 bool ClusterShapeHitFilter::isCompatible(DetId detId,
0371 const SiStripCluster& cluster,
0372 const LocalPoint& lpos,
0373 const LocalVector& ldir) const {
0374 int meas;
0375 float pred;
0376
0377 if (cutOnStripCharge_ && (!checkClusterCharge(detId, cluster, ldir)))
0378 return false;
0379 if (!cutOnStripShape_)
0380 return true;
0381
0382 if (getSizes(detId, cluster, lpos, ldir, meas, pred)) {
0383 StripKeys key(meas);
0384 if (key.isValid())
0385 return stripLimits[key].isInside(pred);
0386 }
0387
0388
0389 return true;
0390 }
0391
0392
0393 bool ClusterShapeHitFilter::isCompatible(DetId detId,
0394 const SiStripCluster& cluster,
0395 const GlobalPoint& gpos,
0396 const GlobalVector& gdir) const {
0397 const GeomDet* det = getsd(detId).det;
0398 LocalVector ldir = det->toLocal(gdir);
0399 LocalPoint lpos = det->toLocal(gpos);
0400
0401 lpos -= ldir * lpos.z() / ldir.z();
0402 return isCompatible(detId, cluster, lpos, ldir);
0403 }
0404 bool ClusterShapeHitFilter::isCompatible(DetId detId, const SiStripCluster& cluster, const GlobalVector& gdir) const {
0405 return isCompatible(detId, cluster, getsd(detId).det->toLocal(gdir));
0406 }
0407
0408 #include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h"
0409
0410 bool ClusterShapeHitFilter::checkClusterCharge(DetId detId,
0411 const SiStripCluster& cluster,
0412 const LocalVector& ldir) const {
0413 return siStripClusterTools::chargePerCM(detId, cluster, ldir) > minGoodStripCharge_;
0414 }
0415
0416 bool ClusterShapeHitFilter::checkClusterCharge(DetId detId,
0417 const SiPixelCluster& cluster,
0418 const LocalVector& ldir) const {
0419 return siStripClusterTools::chargePerCM(detId, cluster, ldir) > minGoodPixelCharge_;
0420 }
0421
0422 #include "FWCore/PluginManager/interface/ModuleDef.h"
0423 #include "FWCore/Framework/interface/MakerMacros.h"
0424
0425 #include "FWCore/Utilities/interface/typelookup.h"
0426 TYPELOOKUP_DATA_REG(ClusterShapeHitFilter);