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