Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-02-20 23:14:35

0001 #include "RecoPixelVertexing/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("RecoPixelVertexing/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     assert(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   //endcap
0134   for (auto det : theTracker->detsPXF()) {
0135     // better not to fail..
0136     const PixelGeomDetUnit* pixelDet = dynamic_cast<const PixelGeomDetUnit*>(det);
0137     assert(pixelDet);
0138     PixelData& pd = pixelData[pixelDet->geographicalId()];
0139     pd.det = pixelDet;
0140     pd.part = 1;
0141     pd.layer = 0;
0142     pd.cotangent = getCotangent(pixelDet);
0143     pd.drift = getDrift(pixelDet);
0144   }
0145 }
0146 
0147 void ClusterShapeHitFilter::fillStripData() {
0148   // copied from StripCPE (FIXME maybe we should move all this in LocalReco)
0149   auto const& geom_ = *theTracker;
0150   auto const& dus = geom_.detUnits();
0151   auto offset = dus.size();
0152   for (unsigned int i = 1; i < 7; ++i) {
0153     if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) != dus.size() &&
0154         dus[geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i])]->type().isTrackerStrip()) {
0155       if (geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]) < offset)
0156         offset = geom_.offsetDU(GeomDetEnumerators::tkDetEnum[i]);
0157     }
0158   }
0159 
0160   for (auto i = offset; i != dus.size(); ++i) {
0161     const StripGeomDetUnit* stripdet = (const StripGeomDetUnit*)(dus[i]);
0162     assert(stripdet->index() == int(i));
0163     assert(stripdet->type().isTrackerStrip());  // not pixel
0164     auto const& bounds = stripdet->specificSurface().bounds();
0165     auto detid = stripdet->geographicalId();
0166     auto& p = stripData[detid];
0167     p.det = stripdet;
0168     p.topology = (StripTopology*)(&stripdet->topology());
0169     p.drift = getDrift(stripdet);
0170     p.thickness = bounds.thickness();
0171     p.nstrips = p.topology->nstrips();
0172   }
0173 }
0174 
0175 /*****************************************************************************/
0176 pair<float, float> ClusterShapeHitFilter::getCotangent(const PixelGeomDetUnit* pixelDet) const {
0177   pair<float, float> cotangent;
0178 
0179   cotangent.first = pixelDet->surface().bounds().thickness() / pixelDet->specificTopology().pitch().first;
0180   cotangent.second = pixelDet->surface().bounds().thickness() / pixelDet->specificTopology().pitch().second;
0181 
0182   return cotangent;
0183 }
0184 
0185 /*****************************************************************************/
0186 float ClusterShapeHitFilter::getCotangent(const StripData& sd, const LocalPoint& pos) const {
0187   // FIXME may be problematic in case of RadialStriptolopgy
0188   return sd.thickness / sd.topology->localPitch(pos);
0189 }
0190 
0191 /*****************************************************************************/
0192 pair<float, float> ClusterShapeHitFilter::getDrift(const PixelGeomDetUnit* pixelDet) const {
0193   LocalVector lBfield = (pixelDet->surface()).toLocal(theMagneticField->inTesla(pixelDet->surface().position()));
0194 
0195   double theTanLorentzAnglePerTesla = theSiPixelLorentzAngle->getLorentzAngle(pixelDet->geographicalId().rawId());
0196 
0197   pair<float, float> dir;
0198   dir.first = -theTanLorentzAnglePerTesla * lBfield.y();
0199   dir.second = theTanLorentzAnglePerTesla * lBfield.x();
0200 
0201   return dir;
0202 }
0203 
0204 /*****************************************************************************/
0205 float ClusterShapeHitFilter::getDrift(const StripGeomDetUnit* stripDet) const {
0206   LocalVector lBfield = (stripDet->surface()).toLocal(theMagneticField->inTesla(stripDet->surface().position()));
0207 
0208   double theTanLorentzAnglePerTesla = 0.0;
0209   if (theSiStripLorentzAngle) {
0210     theTanLorentzAnglePerTesla = theSiStripLorentzAngle->getLorentzAngle(stripDet->geographicalId().rawId());
0211   }
0212 
0213   float dir = theTanLorentzAnglePerTesla * lBfield.y();
0214 
0215   return dir;
0216 }
0217 
0218 /*****************************************************************************/
0219 bool ClusterShapeHitFilter::isNormalOriented(const GeomDetUnit* geomDet) const {
0220   if (geomDet->type().isBarrel()) {  // barrel
0221     float perp0 = geomDet->toGlobal(Local3DPoint(0., 0., 0.)).perp();
0222     float perp1 = geomDet->toGlobal(Local3DPoint(0., 0., 1.)).perp();
0223     return (perp1 > perp0);
0224   } else {  // endcap
0225     float rot = geomDet->toGlobal(LocalVector(0., 0., 1.)).z();
0226     float pos = geomDet->toGlobal(Local3DPoint(0., 0., 0.)).z();
0227     return (rot * pos > 0);
0228   }
0229 }
0230 
0231 /*****************************************************************************/
0232 /*****************************************************************************/
0233 
0234 bool ClusterShapeHitFilter::getSizes(const SiPixelRecHit& recHit,
0235                                      const LocalVector& ldir,
0236                                      const SiPixelClusterShapeCache& clusterShapeCache,
0237                                      int& part,
0238                                      ClusterData::ArrayType& meas,
0239                                      pair<float, float>& pred,
0240                                      PixelData const* ipd) const {
0241   // Get detector
0242   const PixelData& pd = getpd(recHit, ipd);
0243 
0244   // Get shape information
0245   const SiPixelClusterShapeData& data = clusterShapeCache.get(recHit.cluster(), pd.det);
0246   bool usable = (data.isStraight() && data.isComplete());
0247 
0248   // Usable?
0249   //if(usable)
0250   {
0251     part = pd.part;
0252 
0253     // Predicted size
0254     pred.first = ldir.x() / ldir.z();
0255     pred.second = ldir.y() / ldir.z();
0256 
0257     SiPixelClusterShapeData::Range sizeRange = data.size();
0258     if (sizeRange.first->second < 0)
0259       pred.second = -pred.second;
0260 
0261     meas.clear();
0262     assert(meas.capacity() >= std::distance(sizeRange.first, sizeRange.second));
0263     for (auto s = sizeRange.first; s != sizeRange.second; ++s) {
0264       meas.push_back_unchecked(*s);
0265     }
0266     if (sizeRange.first->second < 0) {
0267       for (auto& s : meas)
0268         s.second = -s.second;
0269     }
0270 
0271     // Take out drift
0272     std::pair<float, float> const& drift = pd.drift;
0273     pred.first += drift.first;
0274     pred.second += drift.second;
0275 
0276     // Apply cotangent
0277     std::pair<float, float> const& cotangent = pd.cotangent;
0278     pred.first *= cotangent.first;
0279     pred.second *= cotangent.second;
0280   }
0281 
0282   // Usable?
0283   return usable;
0284 }
0285 
0286 bool ClusterShapeHitFilter::isCompatible(const SiPixelRecHit& recHit,
0287                                          const LocalVector& ldir,
0288                                          const SiPixelClusterShapeCache& clusterShapeCache,
0289                                          PixelData const* ipd) const {
0290   // Get detector
0291   if (cutOnPixelCharge_ && (!checkClusterCharge(recHit.geographicalId(), *(recHit.cluster()), ldir)))
0292     return false;
0293   if (!cutOnPixelShape_)
0294     return true;
0295 
0296   const PixelData& pd = getpd(recHit, ipd);
0297 
0298   int part;
0299   ClusterData::ArrayType meas;
0300   pair<float, float> pred;
0301 
0302   PixelLimits const* pl = pd.layer == 1 ? pixelLimitsL1 : pixelLimits;
0303   if (getSizes(recHit, ldir, clusterShapeCache, part, meas, pred, &pd)) {
0304     for (const auto& m : meas) {
0305       PixelKeys key(part, m.first, m.second);
0306       if (!key.isValid())
0307         return true;  // FIXME original logic
0308       if (pl[key].isInside(pred))
0309         return true;
0310     }
0311     // none of the choices worked
0312     return false;
0313   }
0314   // not usable
0315   return true;
0316 }
0317 
0318 bool ClusterShapeHitFilter::isCompatible(const SiPixelRecHit& recHit,
0319                                          const GlobalVector& gdir,
0320                                          const SiPixelClusterShapeCache& clusterShapeCache,
0321                                          PixelData const* ipd) const {
0322   // Get detector
0323   const PixelData& pd = getpd(recHit, ipd);
0324 
0325   LocalVector ldir = pd.det->toLocal(gdir);
0326 
0327   return isCompatible(recHit, ldir, clusterShapeCache, &pd);
0328 }
0329 
0330 /*****************************************************************************/
0331 /*****************************************************************************/
0332 bool ClusterShapeHitFilter::getSizes(DetId id,
0333                                      const SiStripCluster& cluster,
0334                                      const LocalPoint& lpos,
0335                                      const LocalVector& ldir,
0336                                      int& meas,
0337                                      float& pred) const {
0338   // Get detector
0339   auto const& p = getsd(id);
0340 
0341   // Measured width
0342   meas = cluster.amplitudes().size();
0343 
0344   // Usable?
0345   int fs = cluster.firstStrip();
0346   int ns = p.nstrips;
0347   // bool usable = (fs > 1 && fs + meas - 1 < ns);
0348   bool usable = (fs >= 1) & (fs + meas <= ns + 1);
0349 
0350   // Usable?
0351   //if(usable)
0352   {
0353     // Predicted width
0354     pred = ldir.x() / ldir.z();
0355 
0356     // Take out drift
0357     float drift = p.drift;
0358     ;
0359     pred += drift;
0360 
0361     // Apply cotangent
0362     pred *= getCotangent(p, lpos);
0363   }
0364 
0365   return usable;
0366 }
0367 
0368 /*****************************************************************************/
0369 bool ClusterShapeHitFilter::isCompatible(DetId detId,
0370                                          const SiStripCluster& cluster,
0371                                          const LocalPoint& lpos,
0372                                          const LocalVector& ldir) const {
0373   int meas;
0374   float pred;
0375 
0376   if (cutOnStripCharge_ && (!checkClusterCharge(detId, cluster, ldir)))
0377     return false;
0378   if (!cutOnStripShape_)
0379     return true;
0380 
0381   if (getSizes(detId, cluster, lpos, ldir, meas, pred)) {
0382     StripKeys key(meas);
0383     if (key.isValid())
0384       return stripLimits[key].isInside(pred);
0385   }
0386 
0387   // Not usable or no limits
0388   return true;
0389 }
0390 
0391 /*****************************************************************************/
0392 bool ClusterShapeHitFilter::isCompatible(DetId detId,
0393                                          const SiStripCluster& cluster,
0394                                          const GlobalPoint& gpos,
0395                                          const GlobalVector& gdir) const {
0396   const GeomDet* det = getsd(detId).det;
0397   LocalVector ldir = det->toLocal(gdir);
0398   LocalPoint lpos = det->toLocal(gpos);
0399   // now here we do the transformation
0400   lpos -= ldir * lpos.z() / ldir.z();
0401   return isCompatible(detId, cluster, lpos, ldir);
0402 }
0403 bool ClusterShapeHitFilter::isCompatible(DetId detId, const SiStripCluster& cluster, const GlobalVector& gdir) const {
0404   return isCompatible(detId, cluster, getsd(detId).det->toLocal(gdir));
0405 }
0406 
0407 #include "DataFormats/SiStripCluster/interface/SiStripClusterTools.h"
0408 
0409 bool ClusterShapeHitFilter::checkClusterCharge(DetId detId,
0410                                                const SiStripCluster& cluster,
0411                                                const LocalVector& ldir) const {
0412   return siStripClusterTools::chargePerCM(detId, cluster, ldir) > minGoodStripCharge_;
0413 }
0414 
0415 bool ClusterShapeHitFilter::checkClusterCharge(DetId detId,
0416                                                const SiPixelCluster& cluster,
0417                                                const LocalVector& ldir) const {
0418   return siStripClusterTools::chargePerCM(detId, cluster, ldir) > minGoodPixelCharge_;
0419 }
0420 
0421 #include "FWCore/PluginManager/interface/ModuleDef.h"
0422 #include "FWCore/Framework/interface/MakerMacros.h"
0423 
0424 #include "FWCore/Utilities/interface/typelookup.h"
0425 TYPELOOKUP_DATA_REG(ClusterShapeHitFilter);