TrackingAnalysis/src/SimDoublets.cc

0001 #include "SimDataFormats/TrackingAnalysis/interface/SimDoublets.h"
0002
0003 #include "SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
0004 #include "DataFormats/GeometryVector/interface/GlobalVector.h"
0005 #include "DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
0006 #include "DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
0007
0008 namespace simdoublets {
0009
0010   // Function that gets the global position of a RecHit with respect to a reference point.
0011   GlobalPoint getGlobalHitPosition(SiPixelRecHitRef const& recHit, GlobalVector const& referencePosition) {
0012     return (recHit->globalPosition() - referencePosition);
0013   }
0014
0015   // Function that determines the number of skipped layers for a given pair of RecHits.
0016   int getNumSkippedLayers(std::pair<uint8_t, uint8_t> const& layerIds,
0017                           std::pair<SiPixelRecHitRef, SiPixelRecHitRef> const& recHitRefs,
0018                           const TrackerTopology* trackerTopology) {
0019     // Possibility 0: invalid case (outer layer is not the outer one), set to -1 immediately
0020     if (layerIds.first >= layerIds.second) {
0021       return -1;
0022     }
0023
0024     // get the detector Ids of the two RecHits
0025     DetId innerDetId(recHitRefs.first->geographicalId());
0026     DetId outerDetId(recHitRefs.second->geographicalId());
0027
0028     // determine where the RecHits are
0029     bool innerInBarrel = (innerDetId.subdetId() == PixelSubdetector::PixelBarrel);
0030     bool outerInBarrel = (outerDetId.subdetId() == PixelSubdetector::PixelBarrel);
0031     bool innerInBackward = (!innerInBarrel) && (trackerTopology->pxfSide(innerDetId) == 1);
0032     bool outerInBackward = (!outerInBarrel) && (trackerTopology->pxfSide(outerDetId) == 1);
0033     bool innerInForward = (!innerInBarrel) && (!innerInBackward);
0034     bool outerInForward = (!outerInBarrel) && (!outerInBackward);
0035
0036     // Possibility 1: both RecHits lie in the same detector part (barrel, forward or backward)
0037     if ((innerInBarrel && outerInBarrel) || (innerInForward && outerInForward) ||
0038         (innerInBackward && outerInBackward)) {
0039       return (layerIds.second - layerIds.first - 1);
0040     }
0041     // Possibility 2: the inner RecHit is in the barrel while the outer is in either forward or backward
0042     else if (innerInBarrel) {
0043       return (trackerTopology->pxfDisk(outerDetId) - 1);
0044     }
0045     // Possibility 3: invalid case (one is forward and the other in backward), set to -1
0046     else {
0047       return -1;
0048     }
0049   }
0050
0051   // Function that, for a pair of two layers, gives a unique pair Id (innerLayerId * 100 + outerLayerId).
0052   int getLayerPairId(std::pair<uint8_t, uint8_t> const& layerIds) {
0053     // calculate the unique layer pair Id as (innerLayerId * 100 + outerLayerId)
0054     return (layerIds.first * 100 + layerIds.second);
0055   }
0056 }  // end namespace simdoublets
0057
0058 // ------------------------------------------------------------------------------------------------------
0059 // SimDoublets::Doublet class member functions
0060 // ------------------------------------------------------------------------------------------------------
0061
0062 // constructor
0063 SimDoublets::Doublet::Doublet(SimDoublets const& simDoublets,
0064                               size_t const innerIndex,
0065                               size_t const outerIndex,
0066                               const TrackerTopology* trackerTopology)
0067     : trackingParticleRef_(simDoublets.trackingParticle()), beamSpotPosition_(simDoublets.beamSpotPosition()) {
0068   // fill recHits and layers
0069   recHitRefs_ = std::make_pair(simDoublets.recHits(innerIndex), simDoublets.recHits(outerIndex));
0070   layerIds_ = std::make_pair(simDoublets.layerIds(innerIndex), simDoublets.layerIds(outerIndex));
0071
0072   // determine number of skipped layers
0073   numSkippedLayers_ = simdoublets::getNumSkippedLayers(layerIds_, recHitRefs_, trackerTopology);
0074
0075   // determine Id of the layer pair
0076   layerPairId_ = simdoublets::getLayerPairId(layerIds_);
0077 }
0078
0079 GlobalPoint SimDoublets::Doublet::innerGlobalPos() const {
0080   // get the inner RecHit's global position
0081   return simdoublets::getGlobalHitPosition(recHitRefs_.first, beamSpotPosition_);
0082 }
0083
0084 GlobalPoint SimDoublets::Doublet::outerGlobalPos() const {
0085   // get the outer RecHit's global position
0086   return simdoublets::getGlobalHitPosition(recHitRefs_.second, beamSpotPosition_);
0087 }
0088
0089 // ------------------------------------------------------------------------------------------------------
0090 // SimDoublets class member functions
0091 // ------------------------------------------------------------------------------------------------------
0092
0093 // method to sort the RecHits according to the position
0094 void SimDoublets::sortRecHits() {
0095   // get the production vertex of the TrackingParticle
0096   const GlobalVector vertex(trackingParticleRef_->vx(), trackingParticleRef_->vy(), trackingParticleRef_->vz());
0097
0098   // get the vector of squared magnitudes of the global RecHit positions
0099   std::vector<double> recHitMag2;
0100   recHitMag2.reserve(layerIdVector_.size());
0101   for (const auto& recHit : recHitRefVector_) {
0102     // global RecHit position with respect to the production vertex
0103     Global3DPoint globalPosition = simdoublets::getGlobalHitPosition(recHit, vertex);
0104     recHitMag2.push_back(globalPosition.mag2());
0105   }
0106
0107   // find the permutation vector that sort the magnitudes
0108   std::vector<std::size_t> sortedPerm(recHitMag2.size());
0109   std::iota(sortedPerm.begin(), sortedPerm.end(), 0);
0110   std::sort(sortedPerm.begin(), sortedPerm.end(), [&](std::size_t i, std::size_t j) {
0111     return (recHitMag2[i] < recHitMag2[j]);
0112   });
0113
0114   // create the sorted recHitRefVector and the sorted layerIdVector accordingly
0115   SiPixelRecHitRefVector sorted_recHitRefVector;
0116   std::vector<uint8_t> sorted_layerIdVector;
0117   sorted_recHitRefVector.reserve(sortedPerm.size());
0118   sorted_layerIdVector.reserve(sortedPerm.size());
0119   for (size_t i : sortedPerm) {
0120     sorted_recHitRefVector.push_back(recHitRefVector_[i]);
0121     sorted_layerIdVector.push_back(layerIdVector_[i]);
0122   }
0123
0124   // swap them with the class member
0125   recHitRefVector_.swap(sorted_recHitRefVector);
0126   layerIdVector_.swap(sorted_layerIdVector);
0127
0128   // set sorted bool to true
0129   recHitsAreSorted_ = true;
0130 }
0131
0132 // method to produce the true doublets on the fly
0133 std::vector<SimDoublets::Doublet> SimDoublets::getSimDoublets(const TrackerTopology* trackerTopology) const {
0134   // create output vector for the doublets
0135   std::vector<SimDoublets::Doublet> doubletVector;
0136
0137   // confirm that the RecHits are sorted
0138   assert(recHitsAreSorted_);
0139
0140   // check if there are at least two hits
0141   if (numRecHits() < 2) {
0142     return doubletVector;
0143   }
0144
0145   // loop over the RecHits/layer Ids
0146   for (size_t i = 0; i < layerIdVector_.size(); i++) {
0147     uint8_t innerLayerId = layerIdVector_[i];
0148     uint8_t outerLayerId{};
0149     size_t outerLayerStart{layerIdVector_.size()};
0150
0151     // find the next layer Id + at which hit this layer starts
0152     for (size_t j = i + 1; j < layerIdVector_.size(); j++) {
0153       if (innerLayerId != layerIdVector_[j]) {
0154         outerLayerId = layerIdVector_[j];
0155         outerLayerStart = j;
0156         break;
0157       }
0158     }
0159
0160     // build the doublets of the inner hit i with all outer hits in the layer outerLayerId
0161     for (size_t j = outerLayerStart; j < layerIdVector_.size(); j++) {
0162       // break if the hit doesn't belong to the outer layer anymore
0163       if (outerLayerId != layerIdVector_[j]) {
0164         break;
0165       }
0166
0167       doubletVector.push_back(SimDoublets::Doublet(*this, i, j, trackerTopology));
0168     }
0169   }  // end loop over the RecHits/layer Ids
0170
0171   return doubletVector;
0172 }