TrackFindingTMTT/src/Stub.cc

0001 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0002 #include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
0003 #include "Geometry/CommonTopologies/interface/PixelGeomDetUnit.h"
0004 #include "Geometry/CommonTopologies/interface/PixelTopology.h"
0005 #include "DataFormats/Math/interface/deltaPhi.h"
0006 #include "FWCore/Utilities/interface/Exception.h"
0007
0008 #include "L1Trigger/TrackFindingTMTT/interface/Stub.h"
0009 #include "L1Trigger/TrackFindingTMTT/interface/TP.h"
0010 #include "L1Trigger/TrackFindingTMTT/interface/StubKiller.h"
0011 #include "L1Trigger/TrackFindingTMTT/interface/PrintL1trk.h"
0012
0013 #include <iostream>
0014
0015 using namespace std;
0016
0017 namespace tmtt {
0018
0019   //=== Hybrid L1 tracking: stub constructor.
0020
0021   Stub::Stub(const Settings* settings,
0022              unsigned int idStub,
0023              double phi,
0024              double r,
0025              double z,
0026              double bend,
0027              unsigned int iphi,
0028              double alpha,
0029              unsigned int layerId,
0030              unsigned int iPhiSec,
0031              bool psModule,
0032              bool barrel,
0033              bool tiltedBarrel,
0034              float stripPitch,
0035              float stripLength,
0036              unsigned int nStrips)
0037       : index_in_vStubs_(idStub),  // A unique ID to label the stub.
0038         phi_(phi),
0039         r_(r),
0040         z_(z),
0041         bend_(bend),
0042         iphi_(iphi),
0043         alpha_(alpha),
0044         digitalStub_(std::make_unique<DigitalStub>(settings, r, phi, z, iPhiSec)),
0045         layerId_(layerId),
0046         layerIdReduced_(TrackerModule::calcLayerIdReduced(layerId)),
0047         stripPitch_(stripPitch),
0048         stripLength_(stripLength),
0049         nStrips_(nStrips),
0050         psModule_(psModule),
0051         barrel_(barrel),
0052         tiltedBarrel_(tiltedBarrel) {}
0053
0054   //=== TMTT L1 tracking: stub constructor.
0055
0056   Stub::Stub(const TTStubRef& ttStubRef,
0057              unsigned int index_in_vStubs,
0058              const Settings* settings,
0059              const TrackerTopology* trackerTopology,
0060              const TrackerModule* trackerModule,
0061              const DegradeBend* degradeBend,
0062              const StubKiller* stubKiller)
0063       : ttStubRef_(ttStubRef),
0064         settings_(settings),
0065         index_in_vStubs_(index_in_vStubs),
0066         assocTP_(nullptr),  // Initialize in case job is using no MC truth info.
0067         lastDigiStep_(Stub::DigiStage::NONE),
0068         digitizeWarningsOn_(true),
0069         trackerModule_(trackerModule),  // Info about tracker module containing stub
0070         degradeBend_(degradeBend),      // Used to degrade stub bend information.
0071         // Module related variables (need to be stored for Hybrid)
0072         layerId_(trackerModule->layerId()),
0073         layerIdReduced_(trackerModule->layerIdReduced()),
0074         tiltAngle_(trackerModule->tiltAngle()),
0075         stripPitch_(trackerModule->stripPitch()),
0076         stripLength_(trackerModule->stripLength()),
0077         nStrips_(trackerModule->nStrips()),
0078         psModule_(trackerModule->psModule()),
0079         barrel_(trackerModule->barrel()),
0080         tiltedBarrel_(trackerModule->tiltedBarrel()) {
0081     // Get coordinates of stub.
0082     const TTStub<Ref_Phase2TrackerDigi_>* ttStubP = ttStubRef_.get();
0083
0084     const PixelGeomDetUnit* specDet = trackerModule_->specDet();
0085     const PixelTopology* specTopol = trackerModule_->specTopol();
0086     MeasurementPoint measurementPoint = ttStubRef_->clusterRef(0)->findAverageLocalCoordinatesCentered();
0087     LocalPoint clustlp = specTopol->localPosition(measurementPoint);
0088     GlobalPoint pos = specDet->surface().toGlobal(clustlp);
0089
0090     phi_ = pos.phi();
0091     r_ = pos.perp();
0092     z_ = pos.z();
0093
0094     // Get the coordinates of the two clusters that make up this stub, measured in units of strip pitch, and measured
0095     // in the local frame of the sensor. They have a granularity  of 0.5*pitch.
0096     for (unsigned int iClus = 0; iClus <= 1; iClus++) {  // Loop over two clusters in stub.
0097       localU_cluster_[iClus] = ttStubP->clusterRef(iClus)->findAverageLocalCoordinatesCentered().x();
0098       localV_cluster_[iClus] = ttStubP->clusterRef(iClus)->findAverageLocalCoordinatesCentered().y();
0099     }
0100
0101     // Get location of stub in module in units of strip number (or pixel number along finest granularity axis).
0102     // Range from 0 to (nStrips - 1) inclusive.
0103     // N.B. Since iphi is integer, this degrades the granularity by a factor 2. This seems silly, but track fit wants it.
0104     iphi_ = localU_cluster_[0];  // granularity 1*strip (unclear why we want to degrade it ...)
0105
0106     // Determine alpha correction for non-radial strips in endcap 2S modules.
0107     // (If true hit at larger r than stub r by deltaR, then stub phi needs correcting by +alpha*deltaR).
0108     alpha_ = 0.;
0109     if ((not barrel()) && (not psModule())) {
0110       float fracPosInModule = (float(iphi_) - 0.5 * float(nStrips())) / float(nStrips());
0111       float phiRelToModule = trackerModule_->sensorWidth() * fracPosInModule / r_;
0112       if (z_ < 0)
0113         phiRelToModule *= -1;
0114       if (trackerModule_->outerModuleAtSmallerR())
0115         phiRelToModule *= -1;  // Module flipped.
0116       // If true hit at larger r than stub r by deltaR, then stub phi needs correcting by +alpha*deltaR.
0117       alpha_ = -phiRelToModule / r_;
0118     }
0119
0120     // Calculate variables giving ratio of track intercept angle to stub bend.
0121     this->calcDphiOverBend();
0122
0123     // Get stub bend that is available in front-end electronics, where bend is displacement between
0124     // two hits in stubs in units of strip pitch.
0125     bendInFrontend_ = ttStubRef_->bendFE();
0126     if ((not barrel()) && pos.z() > 0)
0127       bendInFrontend_ *= -1;
0128     if (barrel())
0129       bendInFrontend_ *= -1;
0130
0131     // Get stub bend that is available in off-detector electronics, allowing for degredation of
0132     // bend resolution due to bit encoding by FE chip if required.
0133     numMergedBend_ = 1;  // Number of bend values merged into single degraded one.
0134     if (settings->degradeBendRes() == 2) {
0135       float degradedBend;  // degraded bend
0136       // This returns values of degradedBend & numMergedBend_
0137       this->degradeResolution(bendInFrontend_, degradedBend, numMergedBend_);
0138       bend_ = degradedBend;
0139     } else if (settings->degradeBendRes() == 1) {
0140       bend_ = ttStubRef_->bendBE();  // Degraded bend from official CMS recipe.
0141       if ((not barrel()) && pos.z() > 0)
0142         bend_ *= -1;
0143       if (barrel())
0144         bend_ *= -1;
0145     } else {
0146       bend_ = bendInFrontend_;
0147     }
0148
0149     // Fill frontendPass_ flag, indicating if frontend readout electronics will output this stub.
0150     this->setFrontend(stubKiller);
0151
0152     // Calculate bin range along q/Pt axis of r-phi Hough transform array consistent with bend of this stub.
0153     this->calcQoverPtrange();
0154
0155     // Initialize truth info to false in case job is using no MC truth info.
0156     for (unsigned int iClus = 0; iClus <= 1; iClus++) {
0157       assocTPofCluster_[iClus] = nullptr;
0158     }
0159   }
0160
0161   // KF emualtor: stub constructor
0162   Stub::Stub(const TTStubRef& ttStubRef,
0163              double r,
0164              double phi,
0165              double z,
0166              int layerId,
0167              int layerIdReduced,
0168              double stripPitch,
0169              double stripLength,
0170              bool psModule,
0171              bool barrel,
0172              bool tiltedBarrel)
0173       : ttStubRef_(ttStubRef),
0174         settings_(nullptr),
0175         index_in_vStubs_(0),
0176         phi_(phi),
0177         r_(r),
0178         z_(z),
0179         bend_(0.),
0180         dphiOverBend_(0.),
0181         min_qOverPt_bin_(0),
0182         max_qOverPt_bin_(0),
0183         localU_cluster_({{0., 0.}}),
0184         localV_cluster_({{0., 0.}}),
0185         iphi_(0),
0186         alpha_(0.),
0187         frontendPass_(false),
0188         stubFailedDegradeWindow_(false),
0189         bendInFrontend_(0),
0190         numMergedBend_(0),
0191         assocTP_(nullptr),
0192         assocTPs_(set<const TP*>()),
0193         assocTPofCluster_({{nullptr, nullptr}}),
0194         digitalStub_(nullptr),
0195         lastDigiStep_(DigiStage()),
0196         digitizeWarningsOn_(false),
0197         trackerModule_(nullptr),
0198         degradeBend_(nullptr),
0199         layerId_(layerId),
0200         layerIdReduced_(layerIdReduced),
0201         stripPitch_(stripPitch),
0202         stripLength_(stripLength),
0203         nStrips_(0),
0204         psModule_(psModule),
0205         barrel_(barrel),
0206         tiltedBarrel_(tiltedBarrel) {}
0207
0208   //=== Calculate bin range along q/Pt axis of r-phi Hough transform array consistent with bend of this stub.
0209
0210   void Stub::calcQoverPtrange() {
0211     // First determine bin range along q/Pt axis of HT array
0212     // (Use "int" as nasty things happen if multiply "int" and "unsigned int").
0213     const int nbinsPt = (int)settings_->houghNbinsPt();
0214     const int min_array_bin = 0;
0215     const int max_array_bin = nbinsPt - 1;
0216     // Now calculate range of q/Pt bins allowed by bend filter.
0217     float qOverPtMin = this->qOverPtOverBend() * (this->bend() - this->bendCut());
0218     float qOverPtMax = this->qOverPtOverBend() * (this->bend() + this->bendCut());
0219     int houghNbinsPt = settings_->houghNbinsPt();
0220     const float houghMaxInvPt = 1. / settings_->houghMinPt();
0221     float qOverPtBinSize = (2. * houghMaxInvPt) / houghNbinsPt;
0222     if (settings_->shape() == 2 || settings_->shape() == 1 || settings_->shape() == 3)  // Non-square HT cells.
0223       qOverPtBinSize = 2. * houghMaxInvPt / (houghNbinsPt - 1);
0224     // Convert to bin number along q/Pt axis of HT array.
0225     // N.B. For square HT cells, setting "tmp = -0.5" causeas cell to be accepted if q/Pt at its centre is consistent
0226     // with the stub bend. Instead using "tmp = 0.0" accepts cells if q/Pt at any point in cell is consistent with bend.
0227     // So if you use change from -0.5 to 0.0, you have to tighten the bend cut (by ~0.05) to get similar performance.
0228     // Decision to set tmp = 0.0 taken in softare & GP firmware on 9th August 2016.
0229
0230     float tmp = (settings_->shape() == 2 || settings_->shape() == 1 || settings_->shape() == 3) ? 1. : 0.;
0231     int min_bin = std::floor(-tmp + (qOverPtMin + houghMaxInvPt) / qOverPtBinSize);
0232     int max_bin = std::floor(tmp + (qOverPtMax + houghMaxInvPt) / qOverPtBinSize);
0233
0234     // Limit it to range of HT array.
0235     min_bin = max(min_bin, min_array_bin);
0236     max_bin = min(max_bin, max_array_bin);
0237     // If min_bin > max_bin at this stage, it means that the Pt estimated from the bend is below the cutoff for track-finding.
0238     // Keep min_bin > max_bin, so such stubs can be rejected, but set both variables to values inside the HT bin range.
0239     if (min_bin > max_bin) {
0240       min_bin = max_array_bin;
0241       max_bin = min_array_bin;
0242     }
0243     min_qOverPt_bin_ = (unsigned int)min_bin;
0244     max_qOverPt_bin_ = (unsigned int)max_bin;
0245   }
0246
0247   //=== Digitize stub for input to Geographic Processor, with digitized phi coord. measured relative to closest phi sector.
0248   //=== (This approximation is valid if their are an integer number of digitisation bins inside each phi nonant).
0249
0250   void Stub::digitize(unsigned int iPhiSec, Stub::DigiStage digiStep) {
0251     if (settings_->enableDigitize()) {
0252       bool updated = true;
0253       if (not digitalStub_) {
0254         // Digitize stub if not yet done.
0255         digitalStub_ =
0256             std::make_unique<DigitalStub>(settings_, phi_, r_, z_, min_qOverPt_bin_, max_qOverPt_bin_, bend_, iPhiSec);
0257       } else {
0258         // If digitization already done, redo phi digi if phi sector has changed.
0259         updated = digitalStub_->changePhiSec(iPhiSec);
0260       }
0261
0262       // Save CPU by only updating if something has changed.
0263       if (updated || digiStep != lastDigiStep_) {
0264         lastDigiStep_ = digiStep;
0265
0266         // Replace stub coords with those degraded by digitization process.
0267         if (digiStep == DigiStage::GP) {
0268           phi_ = digitalStub_->phi_GP();
0269         } else {
0270           phi_ = digitalStub_->phi_HT_TF();
0271         }
0272         if (digiStep == DigiStage::GP || digiStep == DigiStage::HT) {
0273           r_ = digitalStub_->r_GP_HT();
0274         } else {
0275           r_ = digitalStub_->r_SF_TF();
0276         }
0277         z_ = digitalStub_->z();
0278         bend_ = digitalStub_->bend();
0279
0280         // Update data members that depend on updated coords.
0281         // (Logically part of digitisation, so disable warnings)
0282         digitizeWarningsOn_ = false;
0283         if (digiStep == DigiStage::GP)
0284           this->calcDphiOverBend();
0285         if (digiStep == DigiStage::HT)
0286           this->calcQoverPtrange();
0287         digitizeWarningsOn_ = true;
0288       }
0289     }
0290   }
0291
0292   //=== Degrade assumed stub bend resolution.
0293   //=== And return an integer indicating how many values of bend are merged into this single one.
0294
0295   void Stub::degradeResolution(float bend, float& degradedBend, unsigned int& num) const {
0296     // If TMTT code is tightening official CMS FE stub window cuts, then calculate TMTT stub windows.
0297     float windowFE;
0298     if (settings_->killLowPtStubs()) {
0299       // Window size corresponding to Pt cut used for tracking.
0300       float invPtMax = 1. / (settings_->houghMinPt());
0301       windowFE = invPtMax / std::abs(this->qOverPtOverBend());
0302       // Increase half-indow size to allow for resolution in bend.
0303       windowFE += this->bendCutInFrontend();
0304     } else {
0305       windowFE = rejectedStubBend_;  // TMTT is not tightening windows.
0306     }
0307
0308     degradeBend_->degrade(bend, psModule(), trackerModule_->detId(), windowFE, degradedBend, num);
0309   }
0310
0311   //=== Set flag indicating if stub will be output by front-end readout electronics
0312   //=== (where we can reconfigure the stub window size and rapidity cut).
0313
0314   void Stub::setFrontend(const StubKiller* stubKiller) {
0315     frontendPass_ = true;              // Did stub pass cuts applied in front-end chip
0316     stubFailedDegradeWindow_ = false;  // Did it only fail cuts corresponding to windows encoded in DegradeBend.h?
0317     // Don't use stubs at large eta, since it is impossible to form L1 tracks from them, so they only contribute to combinatorics.
0318     if (std::abs(this->eta()) > settings_->maxStubEta())
0319       frontendPass_ = false;
0320     // Don't use stubs whose Pt is significantly below the Pt cut used in the L1 tracking, allowing for uncertainty in q/Pt due to stub bend resolution.
0321     const float qOverPtCut = 1. / settings_->houghMinPt();
0322     if (settings_->killLowPtStubs()) {
0323       // Apply this cut in the front-end electronics.
0324       if (std::abs(this->bendInFrontend()) - this->bendCutInFrontend() > qOverPtCut / this->qOverPtOverBend())
0325         frontendPass_ = false;
0326     }
0327
0328     if (frontendPass_ && this->bend() == rejectedStubBend_) {
0329       throw cms::Exception(
0330           "BadConfig: FE stub bend window sizes provided in cfg ES source are tighter than those to make the stubs. "
0331           "Please fix them");
0332     }
0333
0334     if (settings_->killLowPtStubs()) {
0335       // Reapply the same cut using the degraded bend information available in the off-detector electronics.
0336       // The reason is  that the bend degredation can move the Pt below the Pt cut, making the stub useless to the off-detector electronics.
0337       if (std::abs(this->bend()) - this->bendCut() > qOverPtCut / this->qOverPtOverBend())
0338         frontendPass_ = false;
0339     }
0340
0341     // Emulate stubs in dead tracker regions..
0342     StubKiller::KillOptions killScenario = static_cast<StubKiller::KillOptions>(settings_->killScenario());
0343     if (killScenario != StubKiller::KillOptions::none) {
0344       bool kill = stubKiller->killStub(ttStubRef_.get());
0345       if (kill)
0346         frontendPass_ = false;
0347     }
0348   }
0349
0350   //=== Function to calculate approximation for dphiOverBendCorrection aka B
0351   double Stub::approxB() {
0352     if (tiltedBarrel()) {
0353       return settings_->bApprox_gradient() * std::abs(z_) / r_ + settings_->bApprox_intercept();
0354     } else {
0355       return barrel() ? 1 : std::abs(z_) / r_;
0356     }
0357   }
0358
0359   //=== Calculate variables giving ratio of track intercept angle to stub bend.
0360
0361   void Stub::calcDphiOverBend() {
0362     // Uses stub (r,z) instead of module (r,z). Logically correct but has negligable effect on results.
0363     if (settings_->useApproxB()) {
0364       float dphiOverBendCorrection_approx_ = approxB();
0365       dphiOverBend_ = trackerModule_->pitchOverSep() * dphiOverBendCorrection_approx_;
0366     } else {
0367       float dphiOverBendCorrection_ = std::abs(cos(this->theta() - trackerModule_->tiltAngle()) / sin(this->theta()));
0368       dphiOverBend_ = trackerModule_->pitchOverSep() * dphiOverBendCorrection_;
0369     }
0370   }
0371
0372   //=== Note which tracking particle(s), if any, produced this stub.
0373   //=== The 1st argument is a map relating TrackingParticles to TP.
0374
0375   void Stub::fillTruth(const map<edm::Ptr<TrackingParticle>, const TP*>& translateTP,
0376                        const edm::Handle<TTStubAssMap>& mcTruthTTStubHandle,
0377                        const edm::Handle<TTClusterAssMap>& mcTruthTTClusterHandle) {
0378     //--- Fill assocTP_ info. If both clusters in this stub were produced by the same single tracking particle, find out which one it was.
0379
0380     bool genuine = mcTruthTTStubHandle->isGenuine(ttStubRef_);  // Same TP contributed to both clusters?
0381     assocTP_ = nullptr;
0382
0383     // Require same TP contributed to both clusters.
0384     if (genuine) {
0385       edm::Ptr<TrackingParticle> tpPtr = mcTruthTTStubHandle->findTrackingParticlePtr(ttStubRef_);
0386       if (translateTP.find(tpPtr) != translateTP.end()) {
0387         assocTP_ = translateTP.at(tpPtr);
0388         // N.B. Since not all tracking particles are stored in InputData::vTPs_, sometimes no match will be found.
0389       }
0390     }
0391
0392     // Fill assocTPs_ info.
0393
0394     if (settings_->stubMatchStrict()) {
0395       // We consider only stubs in which this TP contributed to both clusters.
0396       if (assocTP_ != nullptr)
0397         assocTPs_.insert(assocTP_);
0398
0399     } else {
0400       // We consider stubs in which this TP contributed to either cluster.
0401
0402       for (unsigned int iClus = 0; iClus <= 1; iClus++) {  // Loop over both clusters that make up stub.
0403         const TTClusterRef& ttClusterRef = ttStubRef_->clusterRef(iClus);
0404
0405         // Now identify all TP's contributing to either cluster in stub.
0406         vector<edm::Ptr<TrackingParticle> > vecTpPtr = mcTruthTTClusterHandle->findTrackingParticlePtrs(ttClusterRef);
0407
0408         for (const edm::Ptr<TrackingParticle>& tpPtr : vecTpPtr) {
0409           if (translateTP.find(tpPtr) != translateTP.end()) {
0410             assocTPs_.insert(translateTP.at(tpPtr));
0411             // N.B. Since not all tracking particles are stored in InputData::vTPs_, sometimes no match will be found.
0412           }
0413         }
0414       }
0415     }
0416
0417     //--- Also note which tracking particles produced the two clusters that make up the stub
0418
0419     for (unsigned int iClus = 0; iClus <= 1; iClus++) {  // Loop over both clusters that make up stub.
0420       const TTClusterRef& ttClusterRef = ttStubRef_->clusterRef(iClus);
0421
0422       bool genuineCluster = mcTruthTTClusterHandle->isGenuine(ttClusterRef);  // Only 1 TP made cluster?
0423       assocTPofCluster_[iClus] = nullptr;
0424
0425       // Only consider clusters produced by just one TP.
0426       if (genuineCluster) {
0427         edm::Ptr<TrackingParticle> tpPtr = mcTruthTTClusterHandle->findTrackingParticlePtr(ttClusterRef);
0428
0429         if (translateTP.find(tpPtr) != translateTP.end()) {
0430           assocTPofCluster_[iClus] = translateTP.at(tpPtr);
0431           // N.B. Since not all tracking particles are stored in InputData::vTPs_, sometimes no match will be found.
0432         }
0433       }
0434     }
0435   }
0436 }  // namespace tmtt