Skip to content

Commit 273a64f

Browse files
JanGerritSchulzJanGerritSchulz
committed
Add new class SimDoublets and a corresponding EDProducer
1 parent 11e37dc commit 273a64f

File tree

6 files changed

+664
-0
lines changed

6 files changed

+664
-0
lines changed

‎SimDataFormats/TrackingAnalysis/interface/SimDoublets.h‎

Lines changed: 155 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,155 @@
1+
#ifndef SimDataFormats_TrackingAnalysis_SimDoublets_h
2+
#define SimDataFormats_TrackingAnalysis_SimDoublets_h
3+
4+
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
5+
#include "SimDataFormats/TrackingAnalysis/interface/TrackingParticleFwd.h"
6+
#include "DataFormats/TrackerCommon/interface/TrackerTopology.h"
7+
#include "DataFormats/BeamSpot/interface/BeamSpot.h"
8+
#include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitFwd.h"
9+
10+
#include <cstddef>
11+
#include <cstdint>
12+
#include <utility>
13+
#include <vector>
14+
15+
/** @brief Semi-Monte Carlo truth information used for pixel-tracking opimization.
16+
*
17+
* SimDoublets hold references to all pixel RecHits of a simulated TrackingParticle.
18+
* Ones those RecHits are sorted according to their position relative to the particle vertex
19+
* by the method sortRecHits(), you can create the true doublets of RecHits that the
20+
* TrackingParticle left in the detector. These SimDoublets::Doublet objects can be used
21+
* to optimize the doublet creation in the reconstruction.
22+
*
23+
* The Doublets are generated as the RecHit pairs between two consecutively hit layers.
24+
* I.e., if a TrackingParticle produces
25+
* - 1 hit (A) in 1st layer
26+
* - 2 hits (B, C) in 3rd layer
27+
* - 1 hit (D) in 4th layer
28+
* then, the true Doublets are:
29+
* (A-B), (A-C), (B-D) and (C-D).
30+
* So, neither does it matter that the 2nd layer got "skipped" as there are no hits,
31+
* nor is the Doublet of (A-D) formed since there is a layer with hits in between.
32+
* Doublets are not created between hits within the same layer.
33+
*
34+
* @author Jan Schulz ([email protected])
35+
* @date January 2025
36+
*/
37+
class SimDoublets {
38+
public:
39+
/**
40+
* Sub-class for true doublets of RecHits
41+
* - first hit = inner RecHit
42+
* - second hit = outer RecHit
43+
*/
44+
class Doublet {
45+
public:
46+
// default constructor
47+
Doublet() = default;
48+
49+
// constructor
50+
Doublet(SimDoublets const&, size_t const, size_t const, const TrackerTopology*);
51+
52+
// method to get the layer pair
53+
std::pair<uint8_t, uint8_t> layerIds() const { return layerIds_; }
54+
55+
// method to get the RecHit pair
56+
std::pair<SiPixelRecHitRef, SiPixelRecHitRef> recHits() const { return recHitRefs_; }
57+
58+
// method to get the number of skipped layers
59+
int8_t numSkippedLayers() const { return numSkippedLayers_; }
60+
61+
// method to get the layer pair ID
62+
int16_t layerPairId() const { return layerPairId_; }
63+
64+
// method to get the inner layerId
65+
uint8_t innerLayerId() const { return layerIds_.first; }
66+
67+
// method to get the outer layerId
68+
uint8_t outerLayerId() const { return layerIds_.second; }
69+
70+
// method to get a reference to the inner RecHit
71+
SiPixelRecHitRef innerRecHit() const { return recHitRefs_.first; }
72+
73+
// method to get a reference to the outer RecHit
74+
SiPixelRecHitRef outerRecHit() const { return recHitRefs_.second; }
75+
76+
// method to get the global position of the inner RecHit
77+
GlobalPoint innerGlobalPos() const;
78+
79+
// method to get the global position of the outer RecHit
80+
GlobalPoint outerGlobalPos() const;
81+
82+
private:
83+
TrackingParticleRef trackingParticleRef_; // reference to the TrackingParticle
84+
std::pair<SiPixelRecHitRef, SiPixelRecHitRef> recHitRefs_; // reference pair to RecHits of the Doublet
85+
std::pair<uint8_t, uint8_t> layerIds_; // pair of layer IDs corresponding to the RecHits
86+
int8_t numSkippedLayers_; // number of layers skipped by the Doublet
87+
int16_t layerPairId_; // ID of the layer pair as defined in the reconstruction for the doublets
88+
GlobalVector beamSpotPosition_; // global position of the beam spot (needed to correct the global RecHit position)
89+
};
90+
91+
// default contructor
92+
SimDoublets() = default;
93+
94+
// constructor
95+
SimDoublets(TrackingParticleRef const trackingParticleRef, reco::BeamSpot const& beamSpot)
96+
: trackingParticleRef_(trackingParticleRef), beamSpotPosition_(beamSpot.x0(), beamSpot.y0(), beamSpot.z0()) {}
97+
98+
// method to add a RecHitRef with its layer
99+
void addRecHit(SiPixelRecHitRef const recHitRef, uint8_t const layerId) {
100+
recHitsAreSorted_ = false; // set sorted-bool to false again
101+
102+
// check if the layerId is not present in the layerIdVector yet
103+
if (std::find(layerIdVector_.begin(), layerIdVector_.end(), layerId) == layerIdVector_.end()) {
104+
// if it does not exist, increment number of layers
105+
numLayers_++;
106+
}
107+
108+
// add recHit and layerId to the vectors
109+
recHitRefVector_.push_back(recHitRef);
110+
layerIdVector_.push_back(layerId);
111+
}
112+
113+
// method to get the reference to the TrackingParticle
114+
TrackingParticleRef trackingParticle() const { return trackingParticleRef_; }
115+
116+
// method to get the reference vector to the RecHits
117+
SiPixelRecHitRefVector recHits() const { return recHitRefVector_; }
118+
119+
// method to get a reference to the RecHit at index i
120+
SiPixelRecHitRef recHits(size_t i) const { return recHitRefVector_[i]; }
121+
122+
// method to get the layer id vector
123+
std::vector<uint8_t> layerIds() const { return layerIdVector_; }
124+
125+
// method to get the layer id at index i
126+
uint8_t layerIds(size_t i) const { return layerIdVector_[i]; }
127+
128+
// method to get the beam spot position
129+
GlobalVector beamSpotPosition() const { return beamSpotPosition_; }
130+
131+
// method to get the number of layers
132+
int numLayers() const { return numLayers_; }
133+
134+
// method to get number of RecHits in the SimDoublets
135+
int numRecHits() const { return layerIdVector_.size(); }
136+
137+
// method to sort the RecHits according to the position
138+
void sortRecHits();
139+
140+
// method to produce the SimDoublets from the RecHits
141+
std::vector<Doublet> getSimDoublets(const TrackerTopology* trackerTopology = nullptr) const;
142+
143+
private:
144+
TrackingParticleRef trackingParticleRef_; // reference to the TrackingParticle
145+
SiPixelRecHitRefVector recHitRefVector_; // reference vector to RecHits associated to the TP (sorted afer building)
146+
std::vector<uint8_t> layerIdVector_; // vector of layer IDs corresponding to the RecHits
147+
GlobalVector beamSpotPosition_; // global position of the beam spot (needed to correct the global RecHit position)
148+
bool recHitsAreSorted_{false}; // true if RecHits were sorted
149+
int numLayers_{0}; // number of layers hit by the TrackingParticle
150+
};
151+
152+
// collection of SimDoublets
153+
typedef std::vector<SimDoublets> SimDoubletsCollection;
154+
155+
#endif

‎SimDataFormats/TrackingAnalysis/src/SimDoublets.cc‎

Lines changed: 172 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,172 @@
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+
}

‎SimDataFormats/TrackingAnalysis/src/classes.h‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -3,6 +3,7 @@
33
#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertex.h"
44
#include "SimDataFormats/TrackingAnalysis/interface/TrackingVertexContainer.h"
55
#include "SimDataFormats/TrackingAnalysis/interface/UniqueSimTrackId.h"
6+
#include "SimDataFormats/TrackingAnalysis/interface/SimDoublets.h"
67
#include "DataFormats/TrackReco/interface/Track.h"
78
#include "DataFormats/Common/interface/AssociationMapHelpers.h"
89
#include "DataFormats/Common/interface/Wrapper.h"

‎SimDataFormats/TrackingAnalysis/src/classes_def.xml‎

100755100644
Lines changed: 6 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -40,4 +40,10 @@
4040
<class name="edm::Wrapper<std::unordered_map<UniqueSimTrackId, TrackingParticleRef, UniqueSimTrackIdHash>>"/>
4141
<class name="SimTrackToTPMap"/>
4242
<class name="edm::Wrapper<SimTrackToTPMap>"/>
43+
44+
<class name="SimDoublets" ClassVersion="3">
45+
<version ClassVersion="3" checksum="2009082523"/>
46+
</class>
47+
<class name="std::vector<SimDoublets>"/>
48+
<class name="edm::Wrapper< std::vector<SimDoublets> >"/>
4349
</lcgdict>

‎SimTracker/TrackerHitAssociation/plugins/BuildFile.xml‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,6 @@
11
<use name="SimTracker/TrackerHitAssociation"/>
22
<use name="SimDataFormats/TrackingAnalysis"/>
3+
<use name="Geometry/Records"/>
34
<library file="*.cc" name="SimTrackerTrackerHitAssociationPlugins">
45
<flags EDM_PLUGIN="1"/>
56
</library>

0 commit comments

Comments
 (0)