MuonIdentification/src/MuonTimingFiller.cc

0001 //
0002 // Package:    MuonTimingFiller
0003 // Class:      MuonTimingFiller
0004 //
0005 /**\class MuonTimingFiller MuonTimingFiller.cc RecoMuon/MuonIdentification/src/MuonTimingFiller.cc
0006
0007  Description: <one line class summary>
0008
0009  Implementation:
0010      <Notes on implementation>
0011 */
0012 //
0013 // Original Author:  Piotr Traczyk, CERN
0014 //         Created:  Mon Mar 16 12:27:22 CET 2009
0015 //
0016 //
0017
0018 // system include files
0019 #include <memory>
0020
0021 // user include files
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 // constructors and destructor
0040 //
0041 MuonTimingFiller::MuonTimingFiller(const edm::ParameterSet& iConfig, edm::ConsumesCollector&& iC) {
0042   // Load parameters for the DTTimingExtractor
0043   edm::ParameterSet dtTimingParameters = iConfig.getParameter<edm::ParameterSet>("DTTimingParameters");
0044
0045   // Load parameters for the CSCTimingExtractor
0046   edm::ParameterSet cscTimingParameters = iConfig.getParameter<edm::ParameterSet>("CSCTimingParameters");
0047
0048   // Fallback mechanism for old configs (there the segment matcher was built inside the timing extractors)
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 // member functions
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             // Use only the segments that passed arbitration to avoid mixing
0097             // segments from in-time and out-of-time muons that may bias the result
0098             // SegmentAndTrackArbitration
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   // Fill DT-specific timing information block
0115   fillTimeFromMeasurements(dtTmSeq, dtTime);
0116
0117   // Fill CSC-specific timing information block
0118   fillTimeFromMeasurements(cscTmSeq, cscTime);
0119
0120   // Fill RPC-specific timing information block
0121   fillRPCTime(muon, rpcTime, iEvent);
0122
0123   // Combine the TimeMeasurementSequences from DT/CSC subdetectors
0124   TimeMeasurementSequence combinedTmSeq;
0125   combineTMSequences(muon, dtTmSeq, cscTmSeq, combinedTmSeq);
0126   // add ECAL info
0127   if (useECAL_)
0128     addEcalTime(muon, combinedTmSeq);
0129
0130   // Fill the master timing block
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   // currently unused
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   // Cut on the crystal energy and restrict to the ECAL barrel for now
0257   //  if (muonE.emMax<ecalEcut_ || fabs(muon.eta())>1.5) return;
0258   if (muonE.emMax < ecalEcut_)
0259     return;
0260
0261   // A simple parametrization of the error on the ECAL time measurement
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 }