MuonIdentification/src/CSCTimingExtractor.cc

0001 // -*- C++ -*-
0002 //
0003 // Package:    MuonIdentification
0004 // Class:      CSCTimingExtractor
0005 //
0006 /**\class CSCTimingExtractor CSCTimingExtractor.cc RecoMuon/MuonIdentification/src/CSCTimingExtractor.cc
0007  *
0008  * Description: Produce timing information for a muon track using CSC hits from segments used to build the track
0009  *
0010  */
0011 //
0012 // Original Author:  Traczyk Piotr
0013 //         Created:  Thu Oct 11 15:01:28 CEST 2007
0014 //
0015 //
0016
0017 #include "RecoMuon/MuonIdentification/interface/CSCTimingExtractor.h"
0018
0019 // user include files
0020 #include "FWCore/Framework/interface/ConsumesCollector.h"
0021 #include "FWCore/Framework/interface/Frameworkfwd.h"
0022
0023 #include "FWCore/Framework/interface/Event.h"
0024
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026
0027 #include "RecoMuon/TrackingTools/interface/MuonServiceProxy.h"
0028
0029 #include "Geometry/Records/interface/MuonGeometryRecord.h"
0030 #include "DataFormats/MuonDetId/interface/MuonSubdetId.h"
0031
0032 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0033 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0034 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0035
0036 #include "DataFormats/MuonReco/interface/Muon.h"
0037 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0038 #include "DataFormats/CSCRecHit/interface/CSCSegment.h"
0039 #include "DataFormats/MuonDetId/interface/CSCDetId.h"
0040 #include "DataFormats/CSCRecHit/interface/CSCRecHit2D.h"
0041 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0042 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0043
0044 #include "DataFormats/TrackReco/interface/Track.h"
0045 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0046 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0047 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
0048
0049 // system include files
0050 #include <memory>
0051 #include <vector>
0052 #include <string>
0053 #include <iostream>
0054
0055 namespace edm {
0056   class ParameterSet;
0057   class EventSetup;
0058   class InputTag;
0059 }  // namespace edm
0060
0061 class MuonServiceProxy;
0062
0063 //
0064 // constructors and destructor
0065 //
0066 CSCTimingExtractor::CSCTimingExtractor(const edm::ParameterSet& iConfig,
0067                                        MuonSegmentMatcher* segMatcher,
0068                                        edm::ConsumesCollector& iC)
0069     : thePruneCut_(iConfig.getParameter<double>("PruneCut")),
0070       theStripTimeOffset_(iConfig.getParameter<double>("CSCStripTimeOffset")),
0071       theWireTimeOffset_(iConfig.getParameter<double>("CSCWireTimeOffset")),
0072       theStripError_(iConfig.getParameter<double>("CSCStripError")),
0073       theWireError_(iConfig.getParameter<double>("CSCWireError")),
0074       UseWireTime(iConfig.getParameter<bool>("UseWireTime")),
0075       UseStripTime(iConfig.getParameter<bool>("UseStripTime")),
0076       debug(iConfig.getParameter<bool>("debug")),
0077       thePropagatorToken(iC.esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))) {
0078   edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
0079   theService = std::make_unique<MuonServiceProxy>(serviceParameters, edm::ConsumesCollector(iC));
0080   theMatcher = segMatcher;
0081 }
0082
0083 CSCTimingExtractor::~CSCTimingExtractor() {}
0084
0085 //
0086 // member functions
0087 //
0088
0089 void CSCTimingExtractor::fillTiming(TimeMeasurementSequence& tmSequence,
0090                                     const std::vector<const CSCSegment*>& segments,
0091                                     reco::TrackRef muonTrack,
0092                                     const edm::Event& iEvent,
0093                                     const edm::EventSetup& iSetup) {
0094   theService->update(iSetup);
0095
0096   const GlobalTrackingGeometry* theTrackingGeometry = &*theService->trackingGeometry();
0097
0098   // get the propagator
0099   const Propagator* propag = &iSetup.getData(thePropagatorToken);
0100
0101   math::XYZPoint pos = muonTrack->innerPosition();
0102   math::XYZVector mom = muonTrack->innerMomentum();
0103   if (sqrt(muonTrack->innerPosition().mag2()) > sqrt(muonTrack->outerPosition().mag2())) {
0104     pos = muonTrack->outerPosition();
0105     mom = -1 * muonTrack->outerMomentum();
0106   }
0107   GlobalPoint posp(pos.x(), pos.y(), pos.z());
0108   GlobalVector momv(mom.x(), mom.y(), mom.z());
0109   FreeTrajectoryState muonFTS(posp, momv, (TrackCharge)muonTrack->charge(), theService->magneticField().product());
0110
0111   // create a collection on TimeMeasurements for the track
0112   std::vector<TimeMeasurement> tms;
0113   for (const auto& rechit : segments) {
0114     // Create the ChamberId
0115     DetId id = rechit->geographicalId();
0116     CSCDetId chamberId(id.rawId());
0117     //    int station = chamberId.station();
0118
0119     if (rechit->specificRecHits().empty())
0120       continue;
0121
0122     const std::vector<CSCRecHit2D>& hits2d(rechit->specificRecHits());
0123
0124     // store all the hits from the segment
0125     for (const auto& hiti : hits2d) {
0126       const GeomDet* cscDet = theTrackingGeometry->idToDet(hiti.geographicalId());
0127       TimeMeasurement thisHit;
0128
0129       std::pair<TrajectoryStateOnSurface, double> tsos;
0130       tsos = propag->propagateWithPath(muonFTS, cscDet->surface());
0131
0132       double dist;
0133       if (tsos.first.isValid())
0134         dist = tsos.second + posp.mag();
0135       else
0136         dist = cscDet->toGlobal(hiti.localPosition()).mag();
0137
0138       thisHit.distIP = dist;
0139       if (UseStripTime) {
0140         thisHit.weightInvbeta = dist * dist / (theStripError_ * theStripError_ * 30. * 30.);
0141         thisHit.weightTimeVtx = 1. / (theStripError_ * theStripError_);
0142         thisHit.timeCorr = hiti.tpeak() - theStripTimeOffset_;
0143         tms.push_back(thisHit);
0144       }
0145
0146       if (UseWireTime) {
0147         thisHit.weightInvbeta = dist * dist / (theWireError_ * theWireError_ * 30. * 30.);
0148         thisHit.weightTimeVtx = 1. / (theWireError_ * theWireError_);
0149         thisHit.timeCorr = hiti.wireTime() - theWireTimeOffset_;
0150         tms.push_back(thisHit);
0151       }
0152
0153       //      std::cout << " CSC Hit. Dist= " << dist << "    Time= " << thisHit.timeCorr
0154       //           << "   invBeta= " << (1.+thisHit.timeCorr/dist*30.) << std::endl;
0155     }
0156
0157   }  // rechit
0158
0159   bool modified = false;
0160   std::vector<double> dstnc, local_t0, hitWeightInvbeta, hitWeightTimeVtx;
0161   double totalWeightInvbeta = 0;
0162   double totalWeightTimeVtx = 0;
0163
0164   // Now loop over the measurements, calculate 1/beta and cut away outliers
0165   do {
0166     modified = false;
0167     dstnc.clear();
0168     local_t0.clear();
0169     hitWeightInvbeta.clear();
0170     hitWeightTimeVtx.clear();
0171
0172     totalWeightInvbeta = 0;
0173     totalWeightTimeVtx = 0;
0174
0175     for (auto& tm : tms) {
0176       dstnc.push_back(tm.distIP);
0177       local_t0.push_back(tm.timeCorr);
0178       hitWeightInvbeta.push_back(tm.weightInvbeta);
0179       hitWeightTimeVtx.push_back(tm.weightTimeVtx);
0180       totalWeightInvbeta += tm.weightInvbeta;
0181       totalWeightTimeVtx += tm.weightTimeVtx;
0182     }
0183
0184     if (totalWeightInvbeta == 0)
0185       break;
0186
0187     // calculate the value and error of 1/beta and timeVtx from the complete set of 1D hits
0188     if (debug)
0189       std::cout << " Points for global fit: " << dstnc.size() << std::endl;
0190
0191     double invbeta = 0, invbetaErr = 0;
0192     double timeVtx = 0, timeVtxErr = 0;
0193
0194     for (unsigned int i = 0; i < dstnc.size(); i++) {
0195       invbeta += (1. + local_t0.at(i) / dstnc.at(i) * 30.) * hitWeightInvbeta.at(i) / totalWeightInvbeta;
0196       timeVtx += local_t0.at(i) * hitWeightTimeVtx.at(i) / totalWeightTimeVtx;
0197     }
0198
0199     double chimax = 0.;
0200     std::vector<TimeMeasurement>::iterator tmmax;
0201
0202     // Calculate the inv beta and time at vertex dispersion
0203     double diff_ibeta, diff_tvtx;
0204     for (unsigned int i = 0; i < dstnc.size(); i++) {
0205       diff_ibeta = (1. + local_t0.at(i) / dstnc.at(i) * 30.) - invbeta;
0206       diff_ibeta = diff_ibeta * diff_ibeta * hitWeightInvbeta.at(i);
0207       diff_tvtx = local_t0.at(i) - timeVtx;
0208       diff_tvtx = diff_tvtx * diff_tvtx * hitWeightTimeVtx.at(i);
0209       invbetaErr += diff_ibeta;
0210       timeVtxErr += diff_tvtx;
0211
0212       // decide if we cut away time at vertex outliers or inverse beta outliers
0213       // currently not configurable.
0214       if (diff_tvtx > chimax) {
0215         tmmax = tms.begin() + i;
0216         chimax = diff_tvtx;
0217       }
0218     }
0219
0220     // cut away the outliers
0221     if (chimax > thePruneCut_) {
0222       tms.erase(tmmax);
0223       modified = true;
0224     }
0225
0226     if (debug) {
0227       double cf = 1. / (dstnc.size() - 1);
0228       invbetaErr = sqrt(invbetaErr / totalWeightInvbeta * cf);
0229       timeVtxErr = sqrt(timeVtxErr / totalWeightTimeVtx * cf);
0230       std::cout << " Measured 1/beta: " << invbeta << " +/- " << invbetaErr << std::endl;
0231       std::cout << " Measured time: " << timeVtx << " +/- " << timeVtxErr << std::endl;
0232     }
0233
0234   } while (modified);
0235
0236   for (unsigned int i = 0; i < dstnc.size(); i++) {
0237     tmSequence.dstnc.push_back(dstnc.at(i));
0238     tmSequence.local_t0.push_back(local_t0.at(i));
0239     tmSequence.weightInvbeta.push_back(hitWeightInvbeta.at(i));
0240     tmSequence.weightTimeVtx.push_back(hitWeightTimeVtx.at(i));
0241   }
0242
0243   tmSequence.totalWeightInvbeta = totalWeightInvbeta;
0244   tmSequence.totalWeightTimeVtx = totalWeightTimeVtx;
0245 }
0246
0247 // ------------ method called to produce the data  ------------
0248 void CSCTimingExtractor::fillTiming(TimeMeasurementSequence& tmSequence,
0249                                     reco::TrackRef muonTrack,
0250                                     const edm::Event& iEvent,
0251                                     const edm::EventSetup& iSetup) {
0252   if (debug)
0253     std::cout << " *** CSC Timimng Extractor ***" << std::endl;
0254
0255   // get the CSC segments that were used to construct the muon
0256   std::vector<const CSCSegment*> range = theMatcher->matchCSC(*muonTrack, iEvent);
0257
0258   fillTiming(tmSequence, range, muonTrack, iEvent, iSetup);
0259 }
0260
0261 //define this as a plug-in
0262 //DEFINE_FWK_MODULE(CSCTimingExtractor);