1+ #include " SimDataFormats/TrackingAnalysis/interface/SimDoublets.h"
2+
3+ #include " SimDataFormats/TrackingAnalysis/interface/TrackingParticle.h"
4+ #include " DataFormats/GeometryVector/interface/GlobalVector.h"
5+ #include " DataFormats/GeometryCommonDetAlgo/interface/MeasurementPoint.h"
6+ #include " DataFormats/SiPixelDetId/interface/PixelSubdetector.h"
7+
8+ namespace simdoublets {
9+
10+ // Function that gets the global position of a RecHit with respect to a reference point.
11+ GlobalPoint getGlobalHitPosition (SiPixelRecHitRef const & recHit, GlobalVector const & referencePosition) {
12+ return (recHit->globalPosition () - referencePosition);
13+ }
14+
15+ // Function that determines the number of skipped layers for a given pair of RecHits.
16+ int getNumSkippedLayers (std::pair<uint8_t , uint8_t > const & layerIds,
17+ std::pair<SiPixelRecHitRef, SiPixelRecHitRef> const & recHitRefs,
18+ const TrackerTopology* trackerTopology) {
19+ // Possibility 0: invalid case (outer layer is not the outer one), set to -1 immediately
20+ if (layerIds.first >= layerIds.second ) {
21+ return -1 ;
22+ }
23+
24+ // get the detector Ids of the two RecHits
25+ DetId innerDetId (recHitRefs.first ->geographicalId ());
26+ DetId outerDetId (recHitRefs.second ->geographicalId ());
27+
28+ // determine where the RecHits are
29+ bool innerInBarrel = (innerDetId.subdetId () == PixelSubdetector::PixelBarrel);
30+ bool outerInBarrel = (outerDetId.subdetId () == PixelSubdetector::PixelBarrel);
31+ bool innerInBackward = (!innerInBarrel) && (trackerTopology->pxfSide (innerDetId) == 1 );
32+ bool outerInBackward = (!outerInBarrel) && (trackerTopology->pxfSide (outerDetId) == 1 );
33+ bool innerInForward = (!innerInBarrel) && (!innerInBackward);
34+ bool outerInForward = (!outerInBarrel) && (!outerInBackward);
35+
36+ // Possibility 1: both RecHits lie in the same detector part (barrel, forward or backward)
37+ if ((innerInBarrel && outerInBarrel) || (innerInForward && outerInForward) ||
38+ (innerInBackward && outerInBackward)) {
39+ return (layerIds.second - layerIds.first - 1 );
40+ }
41+ // Possibility 2: the inner RecHit is in the barrel while the outer is in either forward or backward
42+ else if (innerInBarrel) {
43+ return (trackerTopology->pxfDisk (outerDetId) - 1 );
44+ }
45+ // Possibility 3: invalid case (one is forward and the other in backward), set to -1
46+ else {
47+ return -1 ;
48+ }
49+ }
50+
51+ // Function that, for a pair of two layers, gives a unique pair Id (innerLayerId * 100 + outerLayerId).
52+ int getLayerPairId (std::pair<uint8_t , uint8_t > const & layerIds) {
53+ // calculate the unique layer pair Id as (innerLayerId * 100 + outerLayerId)
54+ return (layerIds.first * 100 + layerIds.second );
55+ }
56+ } // end namespace simdoublets
57+
58+ // ------------------------------------------------------------------------------------------------------
59+ // SimDoublets::Doublet class member functions
60+ // ------------------------------------------------------------------------------------------------------
61+
62+ // constructor
63+ SimDoublets::Doublet::Doublet (SimDoublets const & simDoublets,
64+ size_t const innerIndex,
65+ size_t const outerIndex,
66+ const TrackerTopology* trackerTopology)
67+ : trackingParticleRef_(simDoublets.trackingParticle()), beamSpotPosition_(simDoublets.beamSpotPosition()) {
68+ // fill recHits and layers
69+ recHitRefs_ = std::make_pair (simDoublets.recHits (innerIndex), simDoublets.recHits (outerIndex));
70+ layerIds_ = std::make_pair (simDoublets.layerIds (innerIndex), simDoublets.layerIds (outerIndex));
71+
72+ // determine number of skipped layers
73+ numSkippedLayers_ = simdoublets::getNumSkippedLayers (layerIds_, recHitRefs_, trackerTopology);
74+
75+ // determine Id of the layer pair
76+ layerPairId_ = simdoublets::getLayerPairId (layerIds_);
77+ }
78+
79+ GlobalPoint SimDoublets::Doublet::innerGlobalPos () const {
80+ // get the inner RecHit's global position
81+ return simdoublets::getGlobalHitPosition (recHitRefs_.first , beamSpotPosition_);
82+ }
83+
84+ GlobalPoint SimDoublets::Doublet::outerGlobalPos () const {
85+ // get the outer RecHit's global position
86+ return simdoublets::getGlobalHitPosition (recHitRefs_.second , beamSpotPosition_);
87+ }
88+
89+ // ------------------------------------------------------------------------------------------------------
90+ // SimDoublets class member functions
91+ // ------------------------------------------------------------------------------------------------------
92+
93+ // method to sort the RecHits according to the position
94+ void SimDoublets::sortRecHits () {
95+ // get the production vertex of the TrackingParticle
96+ const GlobalVector vertex (trackingParticleRef_->vx (), trackingParticleRef_->vy (), trackingParticleRef_->vz ());
97+
98+ // get the vector of squared magnitudes of the global RecHit positions
99+ std::vector<double > recHitMag2;
100+ recHitMag2.reserve (layerIdVector_.size ());
101+ for (const auto & recHit : recHitRefVector_) {
102+ // global RecHit position with respect to the production vertex
103+ Global3DPoint globalPosition = simdoublets::getGlobalHitPosition (recHit, vertex);
104+ recHitMag2.push_back (globalPosition.mag2 ());
105+ }
106+
107+ // find the permutation vector that sort the magnitudes
108+ std::vector<std::size_t > sortedPerm (recHitMag2.size ());
109+ std::iota (sortedPerm.begin (), sortedPerm.end (), 0 );
110+ std::sort (sortedPerm.begin (), sortedPerm.end (), [&](std::size_t i, std::size_t j) {
111+ return (recHitMag2[i] < recHitMag2[j]);
112+ });
113+
114+ // create the sorted recHitRefVector and the sorted layerIdVector accordingly
115+ SiPixelRecHitRefVector sorted_recHitRefVector;
116+ std::vector<uint8_t > sorted_layerIdVector;
117+ sorted_recHitRefVector.reserve (sortedPerm.size ());
118+ sorted_layerIdVector.reserve (sortedPerm.size ());
119+ for (size_t i : sortedPerm) {
120+ sorted_recHitRefVector.push_back (recHitRefVector_[i]);
121+ sorted_layerIdVector.push_back (layerIdVector_[i]);
122+ }
123+
124+ // swap them with the class member
125+ recHitRefVector_.swap (sorted_recHitRefVector);
126+ layerIdVector_.swap (sorted_layerIdVector);
127+
128+ // set sorted bool to true
129+ recHitsAreSorted_ = true ;
130+ }
131+
132+ // method to produce the true doublets on the fly
133+ std::vector<SimDoublets::Doublet> SimDoublets::getSimDoublets (const TrackerTopology* trackerTopology) const {
134+ // create output vector for the doublets
135+ std::vector<SimDoublets::Doublet> doubletVector;
136+
137+ // confirm that the RecHits are sorted
138+ assert (recHitsAreSorted_);
139+
140+ // check if there are at least two hits
141+ if (numRecHits () < 2 ) {
142+ return doubletVector;
143+ }
144+
145+ // loop over the RecHits/layer Ids
146+ for (size_t i = 0 ; i < layerIdVector_.size (); i++) {
147+ uint8_t innerLayerId = layerIdVector_[i];
148+ uint8_t outerLayerId{};
149+ size_t outerLayerStart{layerIdVector_.size ()};
150+
151+ // find the next layer Id + at which hit this layer starts
152+ for (size_t j = i + 1 ; j < layerIdVector_.size (); j++) {
153+ if (innerLayerId != layerIdVector_[j]) {
154+ outerLayerId = layerIdVector_[j];
155+ outerLayerStart = j;
156+ break ;
157+ }
158+ }
159+
160+ // build the doublets of the inner hit i with all outer hits in the layer outerLayerId
161+ for (size_t j = outerLayerStart; j < layerIdVector_.size (); j++) {
162+ // break if the hit doesn't belong to the outer layer anymore
163+ if (outerLayerId != layerIdVector_[j]) {
164+ break ;
165+ }
166+
167+ doubletVector.push_back (SimDoublets::Doublet (*this , i, j, trackerTopology));
168+ }
169+ } // end loop over the RecHits/layer Ids
170+
171+ return doubletVector;
172+ }
0 commit comments