PixelLowPtUtilities/src/ClusterShapeHitFilter.cc

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   // Load pixel limits
0048   loadPixelLimits(pixelShapeFile_, pixelLimits);
0049   loadPixelLimits(pixelShapeFileL1_, pixelLimitsL1);
0050   fillPixelData();
0051
0052   // Load strip limits
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;  // 0or 1
0071     inFile >> dx;    // 0 to 10
0072     inFile >> dy;    // 0 to 15 ...
0073
0074     const PixelKeys key(part, dx, dy);
0075     auto& pl = plim[key];
0076
0077     for (int b = 0; b < 2; b++)      // branch
0078       for (int d = 0; d < 2; d++)    // direction
0079         for (int k = 0; k < 2; k++)  // lower and upper
0080           inFile >> pl.data[b][d][k];
0081
0082     double f;
0083     int d;
0084
0085     inFile >> f;  // density
0086     inFile >> d;  // points
0087     inFile >> f;  // density
0088     inFile >> d;  // points
0089   }
0090
0091   inFile.close();
0092
0093   LogTrace("MinBiasTracking|ClusterShapeHitFilter") << " [ClusterShapeHitFilter] pixel-cluster-shape filter loaded";
0094 }
0095
0096 /*****************************************************************************/
0097 void ClusterShapeHitFilter::loadStripLimits() {
0098   // Load strip
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++)    // branch
0110       for (int k = 0; k < 2; k++)  // lower and upper
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   //barrel
0121   for (auto det : theTracker->detsPXB()) {
0122     // better not to fail..
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   //endcap
0135   for (auto det : theTracker->detsPXF()) {
0136     // better not to fail..
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   // copied from StripCPE (FIXME maybe we should move all this in LocalReco)
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());  // not pixel
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   // FIXME may be problematic in case of RadialStriptolopgy
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()) {  // barrel
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 {  // endcap
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   // Get detector
0243   const PixelData& pd = getpd(recHit, ipd);
0244
0245   // Get shape information
0246   const SiPixelClusterShapeData& data = clusterShapeCache.get(recHit.cluster(), pd.det);
0247   bool usable = (data.isStraight() && data.isComplete());
0248
0249   // Usable?
0250   //if(usable)
0251   {
0252     part = pd.part;
0253
0254     // Predicted size
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     // Take out drift
0273     std::pair<float, float> const& drift = pd.drift;
0274     pred.first += drift.first;
0275     pred.second += drift.second;
0276
0277     // Apply cotangent
0278     std::pair<float, float> const& cotangent = pd.cotangent;
0279     pred.first *= cotangent.first;
0280     pred.second *= cotangent.second;
0281   }
0282
0283   // Usable?
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   // Get detector
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;  // FIXME original logic
0309       if (pl[key].isInside(pred))
0310         return true;
0311     }
0312     // none of the choices worked
0313     return false;
0314   }
0315   // not usable
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   // Get detector
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   // Get detector
0340   auto const& p = getsd(id);
0341
0342   // Measured width
0343   meas = cluster.amplitudes().size();
0344
0345   // Usable?
0346   int fs = cluster.firstStrip();
0347   int ns = p.nstrips;
0348   // bool usable = (fs > 1 && fs + meas - 1 < ns);
0349   bool usable = (fs >= 1) & (fs + meas <= ns + 1);
0350
0351   // Usable?
0352   //if(usable)
0353   {
0354     // Predicted width
0355     pred = ldir.x() / ldir.z();
0356
0357     // Take out drift
0358     float drift = p.drift;
0359     ;
0360     pred += drift;
0361
0362     // Apply cotangent
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   // Not usable or no limits
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   // now here we do the transformation
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);