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   //=== Calculate bin range along q/Pt axis of r-phi Hough transform array consistent with bend of this stub.
0162
0163   void Stub::calcQoverPtrange() {
0164     // First determine bin range along q/Pt axis of HT array
0165     // (Use "int" as nasty things happen if multiply "int" and "unsigned int").
0166     const int nbinsPt = (int)settings_->houghNbinsPt();
0167     const int min_array_bin = 0;
0168     const int max_array_bin = nbinsPt - 1;
0169     // Now calculate range of q/Pt bins allowed by bend filter.
0170     float qOverPtMin = this->qOverPtOverBend() * (this->bend() - this->bendCut());
0171     float qOverPtMax = this->qOverPtOverBend() * (this->bend() + this->bendCut());
0172     int houghNbinsPt = settings_->houghNbinsPt();
0173     const float houghMaxInvPt = 1. / settings_->houghMinPt();
0174     float qOverPtBinSize = (2. * houghMaxInvPt) / houghNbinsPt;
0175     if (settings_->shape() == 2 || settings_->shape() == 1 || settings_->shape() == 3)  // Non-square HT cells.
0176       qOverPtBinSize = 2. * houghMaxInvPt / (houghNbinsPt - 1);
0177     // Convert to bin number along q/Pt axis of HT array.
0178     // N.B. For square HT cells, setting "tmp = -0.5" causeas cell to be accepted if q/Pt at its centre is consistent
0179     // with the stub bend. Instead using "tmp = 0.0" accepts cells if q/Pt at any point in cell is consistent with bend.
0180     // 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.
0181     // Decision to set tmp = 0.0 taken in softare & GP firmware on 9th August 2016.
0182
0183     float tmp = (settings_->shape() == 2 || settings_->shape() == 1 || settings_->shape() == 3) ? 1. : 0.;
0184     int min_bin = std::floor(-tmp + (qOverPtMin + houghMaxInvPt) / qOverPtBinSize);
0185     int max_bin = std::floor(tmp + (qOverPtMax + houghMaxInvPt) / qOverPtBinSize);
0186
0187     // Limit it to range of HT array.
0188     min_bin = max(min_bin, min_array_bin);
0189     max_bin = min(max_bin, max_array_bin);
0190     // If min_bin > max_bin at this stage, it means that the Pt estimated from the bend is below the cutoff for track-finding.
0191     // Keep min_bin > max_bin, so such stubs can be rejected, but set both variables to values inside the HT bin range.
0192     if (min_bin > max_bin) {
0193       min_bin = max_array_bin;
0194       max_bin = min_array_bin;
0195     }
0196     min_qOverPt_bin_ = (unsigned int)min_bin;
0197     max_qOverPt_bin_ = (unsigned int)max_bin;
0198   }
0199
0200   //=== Digitize stub for input to Geographic Processor, with digitized phi coord. measured relative to closest phi sector.
0201   //=== (This approximation is valid if their are an integer number of digitisation bins inside each phi nonant).
0202
0203   void Stub::digitize(unsigned int iPhiSec, Stub::DigiStage digiStep) {
0204     if (settings_->enableDigitize()) {
0205       bool updated = true;
0206       if (not digitalStub_) {
0207         // Digitize stub if not yet done.
0208         digitalStub_ =
0209             std::make_unique<DigitalStub>(settings_, phi_, r_, z_, min_qOverPt_bin_, max_qOverPt_bin_, bend_, iPhiSec);
0210       } else {
0211         // If digitization already done, redo phi digi if phi sector has changed.
0212         updated = digitalStub_->changePhiSec(iPhiSec);
0213       }
0214
0215       // Save CPU by only updating if something has changed.
0216       if (updated || digiStep != lastDigiStep_) {
0217         lastDigiStep_ = digiStep;
0218
0219         // Replace stub coords with those degraded by digitization process.
0220         if (digiStep == DigiStage::GP) {
0221           phi_ = digitalStub_->phi_GP();
0222         } else {
0223           phi_ = digitalStub_->phi_HT_TF();
0224         }
0225         if (digiStep == DigiStage::GP || digiStep == DigiStage::HT) {
0226           r_ = digitalStub_->r_GP_HT();
0227         } else {
0228           r_ = digitalStub_->r_SF_TF();
0229         }
0230         z_ = digitalStub_->z();
0231         bend_ = digitalStub_->bend();
0232
0233         // Update data members that depend on updated coords.
0234         // (Logically part of digitisation, so disable warnings)
0235         digitizeWarningsOn_ = false;
0236         if (digiStep == DigiStage::GP)
0237           this->calcDphiOverBend();
0238         if (digiStep == DigiStage::HT)
0239           this->calcQoverPtrange();
0240         digitizeWarningsOn_ = true;
0241       }
0242     }
0243   }
0244
0245   //=== Degrade assumed stub bend resolution.
0246   //=== And return an integer indicating how many values of bend are merged into this single one.
0247
0248   void Stub::degradeResolution(float bend, float& degradedBend, unsigned int& num) const {
0249     // If TMTT code is tightening official CMS FE stub window cuts, then calculate TMTT stub windows.
0250     float windowFE;
0251     if (settings_->killLowPtStubs()) {
0252       // Window size corresponding to Pt cut used for tracking.
0253       float invPtMax = 1. / (settings_->houghMinPt());
0254       windowFE = invPtMax / std::abs(this->qOverPtOverBend());
0255       // Increase half-indow size to allow for resolution in bend.
0256       windowFE += this->bendCutInFrontend();
0257     } else {
0258       windowFE = rejectedStubBend_;  // TMTT is not tightening windows.
0259     }
0260
0261     degradeBend_->degrade(bend, psModule(), trackerModule_->detId(), windowFE, degradedBend, num);
0262   }
0263
0264   //=== Set flag indicating if stub will be output by front-end readout electronics
0265   //=== (where we can reconfigure the stub window size and rapidity cut).
0266
0267   void Stub::setFrontend(const StubKiller* stubKiller) {
0268     frontendPass_ = true;              // Did stub pass cuts applied in front-end chip
0269     stubFailedDegradeWindow_ = false;  // Did it only fail cuts corresponding to windows encoded in DegradeBend.h?
0270     // Don't use stubs at large eta, since it is impossible to form L1 tracks from them, so they only contribute to combinatorics.
0271     if (std::abs(this->eta()) > settings_->maxStubEta())
0272       frontendPass_ = false;
0273     // 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.
0274     const float qOverPtCut = 1. / settings_->houghMinPt();
0275     if (settings_->killLowPtStubs()) {
0276       // Apply this cut in the front-end electronics.
0277       if (std::abs(this->bendInFrontend()) - this->bendCutInFrontend() > qOverPtCut / this->qOverPtOverBend())
0278         frontendPass_ = false;
0279     }
0280
0281     if (frontendPass_ && this->bend() == rejectedStubBend_) {
0282       throw cms::Exception(
0283           "BadConfig: FE stub bend window sizes provided in cfg ES source are tighter than those to make the stubs. "
0284           "Please fix them");
0285     }
0286
0287     if (settings_->killLowPtStubs()) {
0288       // Reapply the same cut using the degraded bend information available in the off-detector electronics.
0289       // The reason is  that the bend degredation can move the Pt below the Pt cut, making the stub useless to the off-detector electronics.
0290       if (std::abs(this->bend()) - this->bendCut() > qOverPtCut / this->qOverPtOverBend())
0291         frontendPass_ = false;
0292     }
0293
0294     // Emulate stubs in dead tracker regions..
0295     StubKiller::KillOptions killScenario = static_cast<StubKiller::KillOptions>(settings_->killScenario());
0296     if (killScenario != StubKiller::KillOptions::none) {
0297       bool kill = stubKiller->killStub(ttStubRef_.get());
0298       if (kill)
0299         frontendPass_ = false;
0300     }
0301   }
0302
0303   //=== Function to calculate approximation for dphiOverBendCorrection aka B
0304   double Stub::approxB() {
0305     if (tiltedBarrel()) {
0306       return settings_->bApprox_gradient() * std::abs(z_) / r_ + settings_->bApprox_intercept();
0307     } else {
0308       return barrel() ? 1 : std::abs(z_) / r_;
0309     }
0310   }
0311
0312   //=== Calculate variables giving ratio of track intercept angle to stub bend.
0313
0314   void Stub::calcDphiOverBend() {
0315     // Uses stub (r,z) instead of module (r,z). Logically correct but has negligable effect on results.
0316     if (settings_->useApproxB()) {
0317       float dphiOverBendCorrection_approx_ = approxB();
0318       dphiOverBend_ = trackerModule_->pitchOverSep() * dphiOverBendCorrection_approx_;
0319     } else {
0320       float dphiOverBendCorrection_ = std::abs(cos(this->theta() - trackerModule_->tiltAngle()) / sin(this->theta()));
0321       dphiOverBend_ = trackerModule_->pitchOverSep() * dphiOverBendCorrection_;
0322     }
0323   }
0324
0325   //=== Note which tracking particle(s), if any, produced this stub.
0326   //=== The 1st argument is a map relating TrackingParticles to TP.
0327
0328   void Stub::fillTruth(const map<edm::Ptr<TrackingParticle>, const TP*>& translateTP,
0329                        const edm::Handle<TTStubAssMap>& mcTruthTTStubHandle,
0330                        const edm::Handle<TTClusterAssMap>& mcTruthTTClusterHandle) {
0331     //--- Fill assocTP_ info. If both clusters in this stub were produced by the same single tracking particle, find out which one it was.
0332
0333     bool genuine = mcTruthTTStubHandle->isGenuine(ttStubRef_);  // Same TP contributed to both clusters?
0334     assocTP_ = nullptr;
0335
0336     // Require same TP contributed to both clusters.
0337     if (genuine) {
0338       edm::Ptr<TrackingParticle> tpPtr = mcTruthTTStubHandle->findTrackingParticlePtr(ttStubRef_);
0339       if (translateTP.find(tpPtr) != translateTP.end()) {
0340         assocTP_ = translateTP.at(tpPtr);
0341         // N.B. Since not all tracking particles are stored in InputData::vTPs_, sometimes no match will be found.
0342       }
0343     }
0344
0345     // Fill assocTPs_ info.
0346
0347     if (settings_->stubMatchStrict()) {
0348       // We consider only stubs in which this TP contributed to both clusters.
0349       if (assocTP_ != nullptr)
0350         assocTPs_.insert(assocTP_);
0351
0352     } else {
0353       // We consider stubs in which this TP contributed to either cluster.
0354
0355       for (unsigned int iClus = 0; iClus <= 1; iClus++) {  // Loop over both clusters that make up stub.
0356         const TTClusterRef& ttClusterRef = ttStubRef_->clusterRef(iClus);
0357
0358         // Now identify all TP's contributing to either cluster in stub.
0359         vector<edm::Ptr<TrackingParticle> > vecTpPtr = mcTruthTTClusterHandle->findTrackingParticlePtrs(ttClusterRef);
0360
0361         for (const edm::Ptr<TrackingParticle>& tpPtr : vecTpPtr) {
0362           if (translateTP.find(tpPtr) != translateTP.end()) {
0363             assocTPs_.insert(translateTP.at(tpPtr));
0364             // N.B. Since not all tracking particles are stored in InputData::vTPs_, sometimes no match will be found.
0365           }
0366         }
0367       }
0368     }
0369
0370     //--- Also note which tracking particles produced the two clusters that make up the stub
0371
0372     for (unsigned int iClus = 0; iClus <= 1; iClus++) {  // Loop over both clusters that make up stub.
0373       const TTClusterRef& ttClusterRef = ttStubRef_->clusterRef(iClus);
0374
0375       bool genuineCluster = mcTruthTTClusterHandle->isGenuine(ttClusterRef);  // Only 1 TP made cluster?
0376       assocTPofCluster_[iClus] = nullptr;
0377
0378       // Only consider clusters produced by just one TP.
0379       if (genuineCluster) {
0380         edm::Ptr<TrackingParticle> tpPtr = mcTruthTTClusterHandle->findTrackingParticlePtr(ttClusterRef);
0381
0382         if (translateTP.find(tpPtr) != translateTP.end()) {
0383           assocTPofCluster_[iClus] = translateTP.at(tpPtr);
0384           // N.B. Since not all tracking particles are stored in InputData::vTPs_, sometimes no match will be found.
0385         }
0386       }
0387     }
0388   }
0389 }  // namespace tmtt