|
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
|
#include "FastSimulation/Muons/plugins/FastTSGFromIOHit.h"
/** \class FastTSGFromIOHit
*
* Emulate TSGFromIOHit in RecoMuon
*
* \author Adam Everett - Purdue University
*/
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/Framework/interface/Event.h"
#include "SimDataFormats/Track/interface/SimTrackContainer.h"
#include "DataFormats/TrackerRecHit2D/interface/FastTrackerRecHit.h"
#include "DataFormats/MuonSeed/interface/L3MuonTrajectorySeed.h"
#include "RecoTracker/TkTrackingRegions/interface/RectangularEtaPhiTrackingRegion.h"
#include "DataFormats/Math/interface/deltaPhi.h"
FastTSGFromIOHit::FastTSGFromIOHit(const edm::ParameterSet& iConfig, edm::ConsumesCollector& iC) {
theCategory = "FastSimulation|Muons||FastTSGFromIOHit";
thePtCut = iConfig.getParameter<double>("PtCut");
simTracksTk = iC.consumes<edm::SimTrackContainer>(iConfig.getParameter<edm::InputTag>("SimTrackCollectionLabel"));
const auto& seedLabels = iConfig.getParameter<std::vector<edm::InputTag> >("SeedCollectionLabels");
for (const auto& seedLabel : seedLabels) {
seedsTks.push_back(iC.consumes<TrajectorySeedCollection>(seedLabel));
}
}
FastTSGFromIOHit::~FastTSGFromIOHit() { LogTrace(theCategory) << " FastTSGFromIOHit dtor called "; }
void FastTSGFromIOHit::trackerSeeds(const TrackCand& staMuon,
const TrackingRegion& region,
const TrackerTopology* tTopo,
std::vector<TrajectorySeed>& result) {
// Make a ref to l2 muon
reco::TrackRef muRef(staMuon.second);
// Cut on muons with low momenta
if (muRef->pt() < thePtCut || muRef->innerMomentum().Rho() < thePtCut || muRef->innerMomentum().R() < 2.5) {
return;
}
// Retrieve the Monte Carlo truth (SimTracks)
edm::Handle<edm::SimTrackContainer> simTracks;
getEvent()->getByToken(simTracksTk, simTracks);
// Retrieve Seed collection
std::vector<edm::Handle<TrajectorySeedCollection> > seedCollections;
seedCollections.resize(seedsTks.size());
for (unsigned iSeed = 0; iSeed < seedsTks.size(); iSeed++) {
getEvent()->getByToken(seedsTks[iSeed], seedCollections[iSeed]);
}
// cast the tracking region
const RectangularEtaPhiTrackingRegion& regionRef = dynamic_cast<const RectangularEtaPhiTrackingRegion&>(region);
// select and store seeds
std::set<unsigned> simTrackIds;
for (const auto& seeds : seedCollections) {
for (const auto& seed : *seeds) {
// Find the simtrack corresponding to the seed
int simTrackId = static_cast<FastTrackerRecHit const&>(*seed.recHits().begin()).simTrackId(0);
const SimTrack& simTrack = (*simTracks)[simTrackId];
// skip if simTrack already associated to a seed
if (simTrackIds.find(simTrackId) != simTrackIds.end()) {
continue;
}
simTrackIds.insert(simTrackId);
// filter seed
if (!clean(muRef, regionRef, &seed, simTrack)) {
continue;
}
// store results
result.push_back(L3MuonTrajectorySeed(seed, muRef));
}
}
}
bool FastTSGFromIOHit::clean(reco::TrackRef muRef,
const RectangularEtaPhiTrackingRegion& region,
const TrajectorySeed* aSeed,
const SimTrack& theSimTrack) {
// Eta cleaner
const PixelRecoRange<float>& etaRange = region.etaRange();
double etaSeed = theSimTrack.momentum().Eta();
double etaLimit = (fabs(fabs(etaRange.max()) - fabs(etaRange.mean())) < 0.05)
? 0.05
: fabs(fabs(etaRange.max()) - fabs(etaRange.mean()));
bool inEtaRange = etaSeed >= (etaRange.mean() - etaLimit) && etaSeed <= (etaRange.mean() + etaLimit);
if (!inEtaRange) {
return false;
}
// Phi cleaner
const TkTrackingRegionsMargin<float>& phiMargin = region.phiMargin();
double phiSeed = theSimTrack.momentum().Phi();
double phiLimit = (phiMargin.right() < 0.05) ? 0.05 : phiMargin.right();
bool inPhiRange = (fabs(deltaPhi(phiSeed, double(region.direction().phi()))) < phiLimit);
if (!inPhiRange) {
return false;
}
// pt cleaner
double ptSeed = std::sqrt(theSimTrack.momentum().Perp2());
double ptMin = (region.ptMin() > 3.5) ? 3.5 : region.ptMin();
bool inPtRange = ptSeed >= ptMin && ptSeed <= 2 * (muRef->pt());
if (!inPtRange) {
return false;
}
return true;
}
|