|
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
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
|
#include "CalibMuon/DTCalibration/interface/DTSegmentSelector.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/Framework/interface/EventSetup.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/Framework/interface/ESHandle.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "DataFormats/DTRecHit/interface/DTRecSegment2D.h"
#include "DataFormats/DTRecHit/interface/DTRecSegment4D.h"
#include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
#include "DataFormats/DTRecHit/interface/DTRecHit1D.h"
#include "CondFormats/DataRecord/interface/DTStatusFlagRcd.h"
#include "CondFormats/DTObjects/interface/DTStatusFlag.h"
DTSegmentSelector::DTSegmentSelector(edm::ParameterSet const& pset, edm::ConsumesCollector& iC)
: muonTags_(pset.getParameter<edm::InputTag>("Muons")),
checkNoisyChannels_(pset.getParameter<bool>("checkNoisyChannels")),
minHitsPhi_(pset.getParameter<int>("minHitsPhi")),
minHitsZ_(pset.getParameter<int>("minHitsZ")),
maxChi2_(pset.getParameter<double>("maxChi2")),
maxAnglePhi_(pset.getParameter<double>("maxAnglePhi")),
maxAngleZ_(pset.getParameter<double>("maxAngleZ")) {
muonToken_ = iC.consumes<reco::MuonCollection>(muonTags_);
theStatusMapToken_ = iC.esConsumes();
}
bool DTSegmentSelector::operator()(DTRecSegment4D const& segment,
edm::Event const& event,
edm::EventSetup const& setup) {
bool result = true;
/* get the muon collection if one is specified
(not specifying a muon collection switches off muon matching */
if (!muonTags_.label().empty()) {
edm::Handle<reco::MuonCollection> muonH;
event.getByToken(muonToken_, muonH);
const std::vector<reco::Muon>* muons = muonH.product();
//std::cout << " Muon collection size: " << muons->size() << std::endl;
if (muons->empty())
return false;
DTChamberId segId(segment.rawId());
bool matched = false;
for (auto& imuon : *muons)
for (const auto& ch : imuon.matches()) {
DetId chId(ch.id.rawId());
if (chId != segId)
continue;
if (imuon.pt() < 15 || !imuon.isGlobalMuon())
continue;
int nsegs = ch.segmentMatches.size();
if (!nsegs)
continue;
LocalPoint posHit = segment.localPosition();
float dx = (posHit.x() ? posHit.x() - ch.x : 0);
float dy = (posHit.y() ? posHit.y() - ch.y : 0);
float dr = sqrt(dx * dx + dy * dy);
if (dr < 5)
matched = true;
}
if (!matched)
result = false;
}
edm::ESHandle<DTStatusFlag> statusMap;
if (checkNoisyChannels_)
statusMap = setup.getHandle(theStatusMapToken_);
// Get the Phi 2D segment
int nPhiHits = -1;
bool segmentNoisyPhi = false;
if (segment.hasPhi()) {
const DTChamberRecSegment2D* phiSeg = segment.phiSegment(); // phiSeg lives in the chamber RF
std::vector<DTRecHit1D> phiRecHits = phiSeg->specificRecHits();
nPhiHits = phiRecHits.size();
if (checkNoisyChannels_)
segmentNoisyPhi = checkNoisySegment(statusMap, phiRecHits);
//} else result = false;
}
// Get the Theta 2D segment
int nZHits = -1;
bool segmentNoisyZ = false;
if (segment.hasZed()) {
const DTSLRecSegment2D* zSeg = segment.zSegment(); // zSeg lives in the SL RF
std::vector<DTRecHit1D> zRecHits = zSeg->specificRecHits();
nZHits = zRecHits.size();
if (checkNoisyChannels_)
segmentNoisyZ = checkNoisySegment(statusMap, zRecHits);
}
// Segment selection
// Discard segment if it has a noisy cell
if (segmentNoisyPhi || segmentNoisyZ)
result = false;
// 2D-segment number of hits
if (segment.hasPhi() && nPhiHits < minHitsPhi_)
result = false;
if (segment.hasZed() && nZHits < minHitsZ_)
result = false;
// Segment chi2
double chiSquare = segment.chi2() / segment.degreesOfFreedom();
if (chiSquare > maxChi2_)
result = false;
// Segment angle
LocalVector segment4DLocalDir = segment.localDirection();
double angleZ = fabs(atan(segment4DLocalDir.y() / segment4DLocalDir.z()) * 180. / Geom::pi());
if (angleZ > maxAngleZ_)
result = false;
double anglePhi = fabs(atan(segment4DLocalDir.x() / segment4DLocalDir.z()) * 180. / Geom::pi());
if (anglePhi > maxAnglePhi_)
result = false;
return result;
}
bool DTSegmentSelector::checkNoisySegment(edm::ESHandle<DTStatusFlag> const& statusMap,
std::vector<DTRecHit1D> const& dtHits) {
bool segmentNoisy = false;
std::vector<DTRecHit1D>::const_iterator dtHit = dtHits.begin();
std::vector<DTRecHit1D>::const_iterator dtHits_end = dtHits.end();
for (; dtHit != dtHits_end; ++dtHit) {
DTWireId wireId = dtHit->wireId();
// Check for noisy channels to skip them
bool isNoisy = false, isFEMasked = false, isTDCMasked = false, isTrigMask = false, isDead = false, isNohv = false;
statusMap->cellStatus(wireId, isNoisy, isFEMasked, isTDCMasked, isTrigMask, isDead, isNohv);
if (isNoisy) {
LogTrace("Calibration") << "Wire: " << wireId << " is noisy, skipping!";
segmentNoisy = true;
break;
}
}
return segmentNoisy;
}
|