VertexFinder/src/Stub.cc

0001 #include "L1Trigger/VertexFinder/interface/Stub.h"
0002
0003 namespace l1tVertexFinder {
0004
0005   //=== Store useful info about this stub.
0006
0007   Stub::Stub()
0008       : settings_(nullptr),
0009         phi_(0.),
0010         r_(0.),
0011         z_(0.),
0012         idDet_(0),
0013         moduleMinR_(0.),
0014         moduleMaxR_(0.),
0015         moduleMinPhi_(0.),
0016         moduleMaxPhi_(0.),
0017         moduleMinZ_(0.),
0018         moduleMaxZ_(0.),
0019         psModule_(false),
0020         layerId_(0),
0021         endcapRing_(0),
0022         barrel_(false),
0023         sigmaPerp_(0.),
0024         sigmaPar_(0.),
0025         stripPitch_(0.),
0026         stripLength_(0.),
0027         nStrips_(0),
0028         sensorWidth_(0.),
0029         outerModuleAtSmallerR_(false) {}
0030
0031   Stub::Stub(const TTStubRef& ttStubRef,
0032              const AnalysisSettings& settings,
0033              const TrackerGeometry* trackerGeometry,
0034              const TrackerTopology* trackerTopology)
0035       : TTStubRef(ttStubRef), settings_(&settings) {
0036     auto geoDetId = trackerGeometry->idToDet(ttStubRef->clusterRef(0)->getDetId())->geographicalId();
0037     auto theGeomDet = trackerGeometry->idToDet(geoDetId);
0038     auto measurementPoint = ttStubRef->clusterRef(0)->findAverageLocalCoordinatesCentered();
0039     auto pos = theGeomDet->surface().toGlobal(theGeomDet->topology().localPosition(measurementPoint));
0040
0041     phi_ = pos.phi();
0042     r_ = pos.perp();
0043     z_ = pos.z();
0044
0045     if (r_ < settings_->trackerInnerRadius() || r_ > settings_->trackerOuterRadius() ||
0046         std::abs(z_) > settings_->trackerHalfLength()) {
0047       throw cms::Exception(
0048           "Stub: Stub found outside assumed tracker volume. Please update tracker dimensions specified in Settings.h!")
0049           << " r=" << r_ << " z=" << z_ << " " << ttStubRef->getDetId().subdetId() << std::endl;
0050     }
0051
0052     // Set info about the module this stub is in
0053     this->setModuleInfo(trackerGeometry, trackerTopology, geoDetId);
0054   }
0055
0056   //=== Note which tracking particle(s), if any, produced this stub.
0057   void Stub::fillTruth(edm::Handle<TTStubAssMap> mcTruthTTStubHandle,
0058                        edm::Handle<TTClusterAssMap> mcTruthTTClusterHandle) {
0059     const TTStubRef& ttStubRef(*this);  // Cast to base class
0060
0061     //--- Fill assocTP_ info. If both clusters in this stub were produced by the same single tracking particle, find out which one it was.
0062
0063     // Require same TP contributed to both clusters.
0064     if (mcTruthTTStubHandle->isGenuine(ttStubRef)) {
0065       assocTP_ = mcTruthTTStubHandle->findTrackingParticlePtr(ttStubRef);
0066     }
0067
0068     // Fill assocTPs_ info.
0069     if (settings_->stubMatchStrict()) {
0070       // We consider only stubs in which this TP contributed to both clusters.
0071       if (!assocTP_.isNull())
0072         assocTPs_.insert(assocTP_);
0073     } else {
0074       // We consider stubs in which this TP contributed to either cluster.
0075
0076       for (unsigned int iClus = 0; iClus <= 1; iClus++) {  // Loop over both clusters that make up stub.
0077         const TTClusterRef& ttClusterRef = ttStubRef->clusterRef(iClus);
0078
0079         // Now identify all TP's contributing to either cluster in stub.
0080         std::vector<edm::Ptr<TrackingParticle>> vecTpPtr =
0081             mcTruthTTClusterHandle->findTrackingParticlePtrs(ttClusterRef);
0082
0083         for (const edm::Ptr<TrackingParticle>& tpPtr : vecTpPtr) {
0084           assocTPs_.insert(tpPtr);
0085         }
0086       }
0087     }
0088   }
0089
0090   void Stub::setModuleInfo(const TrackerGeometry* trackerGeometry,
0091                            const TrackerTopology* trackerTopology,
0092                            const DetId& detId) {
0093     idDet_ = detId();
0094
0095     // Get min & max (r,phi,z) coordinates of the centre of the two sensors containing this stub.
0096     const GeomDetUnit* det0 = trackerGeometry->idToDetUnit(detId);
0097     const GeomDetUnit* det1 = trackerGeometry->idToDetUnit(trackerTopology->partnerDetId(detId));
0098
0099     float R0 = det0->position().perp();
0100     float R1 = det1->position().perp();
0101     float PHI0 = det0->position().phi();
0102     float PHI1 = det1->position().phi();
0103     float Z0 = det0->position().z();
0104     float Z1 = det1->position().z();
0105     moduleMinR_ = std::min(R0, R1);
0106     moduleMaxR_ = std::max(R0, R1);
0107     moduleMinPhi_ = std::min(PHI0, PHI1);
0108     moduleMaxPhi_ = std::max(PHI0, PHI1);
0109     moduleMinZ_ = std::min(Z0, Z1);
0110     moduleMaxZ_ = std::max(Z0, Z1);
0111
0112     // Note if module is PS or 2S, and whether in barrel or endcap.
0113     psModule_ =
0114         trackerGeometry->getDetectorType(detId) ==
0115         TrackerGeometry::ModuleType::
0116             Ph2PSP;  // From https://github.com/cms-sw/cmssw/blob/CMSSW_8_1_X/Geometry/TrackerGeometryBuilder/README.md
0117     barrel_ = detId.subdetId() == StripSubdetector::TOB || detId.subdetId() == StripSubdetector::TIB;
0118
0119     // Encode layer ID.
0120     if (barrel_) {
0121       layerId_ = trackerTopology->layer(detId);  // barrel layer 1-6 encoded as 1-6
0122     } else {
0123       // layerId_ = 10*detId.iSide() + detId.iDisk(); // endcap layer 1-5 encoded as 11-15 (endcap A) or 21-25 (endcapB)
0124       // EJC This seems to give the same encoding
0125       layerId_ = 10 * trackerTopology->side(detId) + trackerTopology->tidWheel(detId);
0126     }
0127
0128     // Note module ring in endcap
0129     // endcapRing_ = barrel_  ?  0  :  detId.iRing();
0130     endcapRing_ = barrel_ ? 0 : trackerTopology->tidRing(detId);
0131
0132     // Get sensor strip or pixel pitch using innermost sensor of pair.
0133
0134     const PixelGeomDetUnit* unit = reinterpret_cast<const PixelGeomDetUnit*>(det0);
0135     const PixelTopology& topo = unit->specificTopology();
0136     const Bounds& bounds = det0->surface().bounds();
0137
0138     std::pair<float, float> pitch = topo.pitch();
0139     stripPitch_ = pitch.first;      // Strip pitch (or pixel pitch along shortest axis)
0140     stripLength_ = pitch.second;    //  Strip length (or pixel pitch along longest axis)
0141     nStrips_ = topo.nrows();        // No. of strips in sensor
0142     sensorWidth_ = bounds.width();  // Width of sensitive region of sensor (= stripPitch * nStrips).
0143
0144     outerModuleAtSmallerR_ = false;
0145     if (barrel_ && det0->position().perp() > det1->position().perp()) {
0146       outerModuleAtSmallerR_ = true;
0147     }
0148
0149     sigmaPerp_ = stripPitch_ / sqrt(12.);  // resolution perpendicular to strip (or to longest pixel axis)
0150     sigmaPar_ = stripLength_ / sqrt(12.);  // resolution parallel to strip (or to longest pixel axis)
0151   }
0152
0153 }  // end namespace l1tVertexFinder