File indexing completed on 2023-03-17 11:20:36
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include <memory>
0020
0021
0022 #include "FWCore/Framework/interface/Frameworkfwd.h"
0023
0024 #include "FWCore/Framework/interface/Event.h"
0025
0026 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0027
0028 #include "DataFormats/MuonReco/interface/Muon.h"
0029 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0030 #include "DataFormats/MuonReco/interface/MuonTimeExtra.h"
0031 #include "DataFormats/MuonReco/interface/MuonTimeExtraMap.h"
0032 #include "DataFormats/RPCRecHit/interface/RPCRecHit.h"
0033
0034 #include "RecoMuon/MuonIdentification/interface/MuonTimingFiller.h"
0035 #include "RecoMuon/MuonIdentification/interface/TimeMeasurementSequence.h"
0036 #include "DataFormats/EcalDetId/interface/EcalSubdetector.h"
0037
0038
0039
0040
0041 MuonTimingFiller::MuonTimingFiller(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC) {
0042
0043 edm::ParameterSet dtTimingParameters = iConfig.getParameter<edm::ParameterSet>("DTTimingParameters");
0044
0045
0046 edm::ParameterSet cscTimingParameters = iConfig.getParameter<edm::ParameterSet>("CSCTimingParameters");
0047
0048
0049 edm::ParameterSet matchParameters;
0050 if (iConfig.existsAs<edm::ParameterSet>("MatchParameters"))
0051 matchParameters = iConfig.getParameter<edm::ParameterSet>("MatchParameters");
0052 else
0053 matchParameters = dtTimingParameters.getParameter<edm::ParameterSet>("MatchParameters");
0054
0055 theMatcher_ = std::make_unique<MuonSegmentMatcher>(matchParameters, iC);
0056 theDTTimingExtractor_ = std::make_unique<DTTimingExtractor>(dtTimingParameters, theMatcher_.get(), iC);
0057 theCSCTimingExtractor_ = std::make_unique<CSCTimingExtractor>(cscTimingParameters, theMatcher_.get(), iC);
0058
0059 errorEB_ = iConfig.getParameter<double>("ErrorEB");
0060 errorEE_ = iConfig.getParameter<double>("ErrorEE");
0061 ecalEcut_ = iConfig.getParameter<double>("EcalEnergyCut");
0062
0063 useDT_ = iConfig.getParameter<bool>("UseDT");
0064 useCSC_ = iConfig.getParameter<bool>("UseCSC");
0065 useECAL_ = iConfig.getParameter<bool>("UseECAL");
0066 }
0067
0068 MuonTimingFiller::~MuonTimingFiller() {}
0069
0070
0071
0072
0073
0074 void MuonTimingFiller::fillTiming(const reco::Muon& muon,
0075 reco::MuonTimeExtra& dtTime,
0076 reco::MuonTimeExtra& cscTime,
0077 reco::MuonTime& rpcTime,
0078 reco::MuonTimeExtra& combinedTime,
0079 edm::Event& iEvent,
0080 const edm::EventSetup& iSetup) {
0081 TimeMeasurementSequence dtTmSeq, cscTmSeq;
0082
0083 if (!(muon.combinedMuon().isNull())) {
0084 theDTTimingExtractor_->fillTiming(dtTmSeq, muon.combinedMuon(), iEvent, iSetup);
0085 theCSCTimingExtractor_->fillTiming(cscTmSeq, muon.combinedMuon(), iEvent, iSetup);
0086 } else {
0087 if (!(muon.standAloneMuon().isNull())) {
0088 theDTTimingExtractor_->fillTiming(dtTmSeq, muon.standAloneMuon(), iEvent, iSetup);
0089 theCSCTimingExtractor_->fillTiming(cscTmSeq, muon.standAloneMuon(), iEvent, iSetup);
0090 } else {
0091 if (muon.isTrackerMuon()) {
0092 std::vector<const DTRecSegment4D*> dtSegments;
0093 std::vector<const CSCSegment*> cscSegments;
0094 for (auto& chamber : muon.matches()) {
0095 for (auto& segment : chamber.segmentMatches) {
0096
0097
0098
0099 if (segment.isMask(reco::MuonSegmentMatch::BestInStationByDR) &&
0100 segment.isMask(reco::MuonSegmentMatch::BelongsToTrackByDR)) {
0101 if (!(segment.dtSegmentRef.isNull()))
0102 dtSegments.push_back(segment.dtSegmentRef.get());
0103 if (!(segment.cscSegmentRef.isNull()))
0104 cscSegments.push_back(segment.cscSegmentRef.get());
0105 }
0106 }
0107 }
0108 theDTTimingExtractor_->fillTiming(dtTmSeq, dtSegments, muon.innerTrack(), iEvent, iSetup);
0109 theCSCTimingExtractor_->fillTiming(cscTmSeq, cscSegments, muon.innerTrack(), iEvent, iSetup);
0110 }
0111 }
0112 }
0113
0114
0115 fillTimeFromMeasurements(dtTmSeq, dtTime);
0116
0117
0118 fillTimeFromMeasurements(cscTmSeq, cscTime);
0119
0120
0121 fillRPCTime(muon, rpcTime, iEvent);
0122
0123
0124 TimeMeasurementSequence combinedTmSeq;
0125 combineTMSequences(muon, dtTmSeq, cscTmSeq, combinedTmSeq);
0126
0127 if (useECAL_)
0128 addEcalTime(muon, combinedTmSeq);
0129
0130
0131 fillTimeFromMeasurements(combinedTmSeq, combinedTime);
0132
0133 LogTrace("MuonTime") << "Global 1/beta: " << combinedTime.inverseBeta() << " +/- " << combinedTime.inverseBetaErr()
0134 << std::endl;
0135 LogTrace("MuonTime") << " Free 1/beta: " << combinedTime.freeInverseBeta() << " +/- "
0136 << combinedTime.freeInverseBetaErr() << std::endl;
0137 LogTrace("MuonTime") << " Vertex time (in-out): " << combinedTime.timeAtIpInOut() << " +/- "
0138 << combinedTime.timeAtIpInOutErr() << " # of points: " << combinedTime.nDof() << std::endl;
0139 LogTrace("MuonTime") << " Vertex time (out-in): " << combinedTime.timeAtIpOutIn() << " +/- "
0140 << combinedTime.timeAtIpOutInErr() << std::endl;
0141 LogTrace("MuonTime") << " direction: " << combinedTime.direction() << std::endl;
0142 }
0143
0144 void MuonTimingFiller::fillTimeFromMeasurements(const TimeMeasurementSequence& tmSeq, reco::MuonTimeExtra& muTime) {
0145 std::vector<double> x, y;
0146 double invbeta(0), invbetaerr(0);
0147 double vertexTime(0), vertexTimeErr(0), vertexTimeR(0), vertexTimeRErr(0);
0148 double freeBeta(0), freeBetaErr(0), freeTime(0), freeTimeErr(0);
0149
0150 if (tmSeq.dstnc.size() <= 1)
0151 return;
0152
0153 for (unsigned int i = 0; i < tmSeq.dstnc.size(); i++) {
0154 invbeta +=
0155 (1. + tmSeq.local_t0.at(i) / tmSeq.dstnc.at(i) * 30.) * tmSeq.weightInvbeta.at(i) / tmSeq.totalWeightInvbeta;
0156 x.push_back(tmSeq.dstnc.at(i) / 30.);
0157 y.push_back(tmSeq.local_t0.at(i) + tmSeq.dstnc.at(i) / 30.);
0158 vertexTime += tmSeq.local_t0.at(i) * tmSeq.weightTimeVtx.at(i) / tmSeq.totalWeightTimeVtx;
0159 vertexTimeR +=
0160 (tmSeq.local_t0.at(i) + 2 * tmSeq.dstnc.at(i) / 30.) * tmSeq.weightTimeVtx.at(i) / tmSeq.totalWeightTimeVtx;
0161 }
0162
0163 double diff;
0164 for (unsigned int i = 0; i < tmSeq.dstnc.size(); i++) {
0165 diff = (1. + tmSeq.local_t0.at(i) / tmSeq.dstnc.at(i) * 30.) - invbeta;
0166 invbetaerr += diff * diff * tmSeq.weightInvbeta.at(i);
0167 diff = tmSeq.local_t0.at(i) - vertexTime;
0168 vertexTimeErr += diff * diff * tmSeq.weightTimeVtx.at(i);
0169 diff = tmSeq.local_t0.at(i) + 2 * tmSeq.dstnc.at(i) / 30. - vertexTimeR;
0170 vertexTimeRErr += diff * diff * tmSeq.weightTimeVtx.at(i);
0171 }
0172
0173 double cf = 1. / (tmSeq.dstnc.size() - 1);
0174 invbetaerr = sqrt(invbetaerr / tmSeq.totalWeightInvbeta * cf);
0175 vertexTimeErr = sqrt(vertexTimeErr / tmSeq.totalWeightTimeVtx * cf);
0176 vertexTimeRErr = sqrt(vertexTimeRErr / tmSeq.totalWeightTimeVtx * cf);
0177
0178 muTime.setInverseBeta(invbeta);
0179 muTime.setInverseBetaErr(invbetaerr);
0180 muTime.setTimeAtIpInOut(vertexTime);
0181 muTime.setTimeAtIpInOutErr(vertexTimeErr);
0182 muTime.setTimeAtIpOutIn(vertexTimeR);
0183 muTime.setTimeAtIpOutInErr(vertexTimeRErr);
0184
0185 rawFit(freeBeta, freeBetaErr, freeTime, freeTimeErr, x, y);
0186
0187 muTime.setFreeInverseBeta(freeBeta);
0188 muTime.setFreeInverseBetaErr(freeBetaErr);
0189
0190 muTime.setNDof(tmSeq.dstnc.size());
0191 }
0192
0193 void MuonTimingFiller::fillRPCTime(const reco::Muon& muon, reco::MuonTime& rpcTime, edm::Event& iEvent) {
0194 double trpc = 0, trpc2 = 0;
0195
0196 reco::TrackRef staTrack = muon.standAloneMuon();
0197 if (staTrack.isNull())
0198 return;
0199
0200 const std::vector<const RPCRecHit*> rpcHits = theMatcher_->matchRPC(*staTrack, iEvent);
0201 const int nrpc = rpcHits.size();
0202 for (const auto& hitRPC : rpcHits) {
0203 const double time = hitRPC->timeError() < 0 ? hitRPC->BunchX() * 25. : hitRPC->time();
0204 trpc += time;
0205 trpc2 += time * time;
0206 }
0207
0208 if (nrpc == 0)
0209 return;
0210
0211 trpc2 = trpc2 / nrpc;
0212 trpc = trpc / nrpc;
0213 const double trpcerr = sqrt(std::max(0., trpc2 - trpc * trpc));
0214
0215 rpcTime.timeAtIpInOut = trpc;
0216 rpcTime.timeAtIpInOutErr = trpcerr;
0217 rpcTime.nDof = nrpc;
0218
0219
0220 rpcTime.timeAtIpOutIn = 0.;
0221 rpcTime.timeAtIpOutInErr = 0.;
0222 }
0223
0224 void MuonTimingFiller::combineTMSequences(const reco::Muon& muon,
0225 const TimeMeasurementSequence& dtSeq,
0226 const TimeMeasurementSequence& cscSeq,
0227 TimeMeasurementSequence& cmbSeq) {
0228 if (useDT_)
0229 for (unsigned int i = 0; i < dtSeq.dstnc.size(); i++) {
0230 cmbSeq.dstnc.push_back(dtSeq.dstnc.at(i));
0231 cmbSeq.local_t0.push_back(dtSeq.local_t0.at(i));
0232 cmbSeq.weightTimeVtx.push_back(dtSeq.weightTimeVtx.at(i));
0233 cmbSeq.weightInvbeta.push_back(dtSeq.weightInvbeta.at(i));
0234
0235 cmbSeq.totalWeightTimeVtx += dtSeq.weightTimeVtx.at(i);
0236 cmbSeq.totalWeightInvbeta += dtSeq.weightInvbeta.at(i);
0237 }
0238
0239 if (useCSC_)
0240 for (unsigned int i = 0; i < cscSeq.dstnc.size(); i++) {
0241 cmbSeq.dstnc.push_back(cscSeq.dstnc.at(i));
0242 cmbSeq.local_t0.push_back(cscSeq.local_t0.at(i));
0243 cmbSeq.weightTimeVtx.push_back(cscSeq.weightTimeVtx.at(i));
0244 cmbSeq.weightInvbeta.push_back(cscSeq.weightInvbeta.at(i));
0245
0246 cmbSeq.totalWeightTimeVtx += cscSeq.weightTimeVtx.at(i);
0247 cmbSeq.totalWeightInvbeta += cscSeq.weightInvbeta.at(i);
0248 }
0249 }
0250
0251 void MuonTimingFiller::addEcalTime(const reco::Muon& muon, TimeMeasurementSequence& cmbSeq) {
0252 reco::MuonEnergy muonE;
0253 if (muon.isEnergyValid())
0254 muonE = muon.calEnergy();
0255
0256
0257
0258 if (muonE.emMax < ecalEcut_)
0259 return;
0260
0261
0262 double emErr;
0263 if (muonE.ecal_id.subdetId() == EcalBarrel)
0264 emErr = errorEB_ / muonE.emMax;
0265 else
0266 emErr = errorEE_ / muonE.emMax;
0267 double hitWeight = 1 / (emErr * emErr);
0268 double hitDist = muonE.ecal_position.r();
0269
0270 cmbSeq.local_t0.push_back(muonE.ecal_time);
0271 cmbSeq.weightTimeVtx.push_back(hitWeight);
0272 cmbSeq.weightInvbeta.push_back(hitDist * hitDist * hitWeight / (30. * 30.));
0273
0274 cmbSeq.dstnc.push_back(hitDist);
0275
0276 cmbSeq.totalWeightTimeVtx += hitWeight;
0277 cmbSeq.totalWeightInvbeta += hitDist * hitDist * hitWeight / (30. * 30.);
0278 }
0279
0280 void MuonTimingFiller::rawFit(double& freeBeta,
0281 double& freeBetaErr,
0282 double& freeTime,
0283 double& freeTimeErr,
0284 const std::vector<double>& hitsx,
0285 const std::vector<double>& hitsy) {
0286 double s = 0, sx = 0, sy = 0, x, y;
0287 double sxx = 0, sxy = 0;
0288
0289 freeBeta = 0;
0290 freeBetaErr = 0;
0291 freeTime = 0;
0292 freeTimeErr = 0;
0293
0294 if (hitsx.empty())
0295 return;
0296 if (hitsx.size() == 1) {
0297 freeTime = hitsy[0];
0298 } else {
0299 for (unsigned int i = 0; i != hitsx.size(); i++) {
0300 x = hitsx[i];
0301 y = hitsy[i];
0302 sy += y;
0303 sxy += x * y;
0304 s += 1.;
0305 sx += x;
0306 sxx += x * x;
0307 }
0308
0309 double d = s * sxx - sx * sx;
0310 freeTime = (sxx * sy - sx * sxy) / d;
0311 freeBeta = (s * sxy - sx * sy) / d;
0312 freeBetaErr = sqrt(sxx / d);
0313 freeTimeErr = sqrt(s / d);
0314 }
0315 }