Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:27:02 

0001 // -*- C++ -*-
0002 //
0003 // Package:    MuonIdentification
0004 // Class:      DTTimingExtractor
0005 //
0006 /**\class DTTimingExtractor DTTimingExtractor.cc RecoMuon/MuonIdentification/src/DTTimingExtractor.cc
0007  *
0008  * Description: Produce timing information for a muon track using 1D DT 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/DTTimingExtractor.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 #include "Geometry/DTGeometry/interface/DTGeometry.h"
0032 #include "Geometry/DTGeometry/interface/DTLayer.h"
0033 #include "Geometry/DTGeometry/interface/DTChamber.h"
0034 #include "Geometry/DTGeometry/interface/DTSuperLayer.h"
0035 #include "DataFormats/DTRecHit/interface/DTSLRecSegment2D.h"
0036 #include "DataFormats/DTRecHit/interface/DTRecSegment2D.h"
0037 
0038 #include "Geometry/Records/interface/GlobalTrackingGeometryRecord.h"
0039 #include "Geometry/CommonDetUnit/interface/GlobalTrackingGeometry.h"
0040 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0041 
0042 #include "DataFormats/MuonReco/interface/Muon.h"
0043 #include "DataFormats/MuonReco/interface/MuonFwd.h"
0044 #include "DataFormats/DTRecHit/interface/DTRecSegment2D.h"
0045 #include "DataFormats/DTRecHit/interface/DTRecSegment4DCollection.h"
0046 #include "DataFormats/DTRecHit/interface/DTRecSegment2DCollection.h"
0047 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0048 #include "RecoMuon/TransientTrackingRecHit/interface/MuonTransientTrackingRecHit.h"
0049 
0050 #include "DataFormats/TrackReco/interface/Track.h"
0051 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0052 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0053 #include "DataFormats/TrackReco/interface/TrackExtraFwd.h"
0054 
0055 // system include files
0056 #include <memory>
0057 #include <vector>
0058 #include <string>
0059 #include <iostream>
0060 
0061 namespace edm {
0062   class ParameterSet;
0063   class EventSetup;
0064   class InputTag;
0065 }  // namespace edm
0066 
0067 class MuonServiceProxy;
0068 
0069 //
0070 // constructors and destructor
0071 //
0072 DTTimingExtractor::DTTimingExtractor(const edm::ParameterSet& iConfig,
0073                                      MuonSegmentMatcher* segMatcher,
0074                                      edm::ConsumesCollector& iC)
0075     : theHitsMin_(iConfig.getParameter<int>("HitsMin")),
0076       thePruneCut_(iConfig.getParameter<double>("PruneCut")),
0077       theTimeOffset_(iConfig.getParameter<double>("DTTimeOffset")),
0078       theError_(iConfig.getParameter<double>("HitError")),
0079       useSegmentT0_(iConfig.getParameter<bool>("UseSegmentT0")),
0080       doWireCorr_(iConfig.getParameter<bool>("DoWireCorr")),
0081       dropTheta_(iConfig.getParameter<bool>("DropTheta")),
0082       requireBothProjections_(iConfig.getParameter<bool>("RequireBothProjections")),
0083       debug(iConfig.getParameter<bool>("debug")),
0084       theDTGeomToken(iC.esConsumes()),
0085       thePropagatorToken(iC.esConsumes(edm::ESInputTag("", "SteppingHelixPropagatorAny"))) {
0086   edm::ParameterSet serviceParameters = iConfig.getParameter<edm::ParameterSet>("ServiceParameters");
0087   theService = std::make_unique<MuonServiceProxy>(serviceParameters, edm::ConsumesCollector(iC));
0088   theMatcher = segMatcher;
0089 }
0090 
0091 DTTimingExtractor::~DTTimingExtractor() {}
0092 
0093 //
0094 // member functions
0095 //
0096 void DTTimingExtractor::fillTiming(TimeMeasurementSequence& tmSequence,
0097                                    const std::vector<const DTRecSegment4D*>& segments,
0098                                    reco::TrackRef muonTrack,
0099                                    const edm::Event& iEvent,
0100                                    const edm::EventSetup& iSetup) {
0101   if (debug)
0102     std::cout << " *** DT Timimng Extractor ***" << std::endl;
0103 
0104   theService->update(iSetup);
0105 
0106   const GlobalTrackingGeometry* theTrackingGeometry = &*theService->trackingGeometry();
0107 
0108   // get the DT geometry
0109   DTGeometry const& theDTGeom = iSetup.getData(theDTGeomToken);
0110 
0111   // get the propagator
0112   const Propagator* propag = &iSetup.getData(thePropagatorToken);
0113 
0114   math::XYZPoint pos = muonTrack->innerPosition();
0115   math::XYZVector mom = muonTrack->innerMomentum();
0116   GlobalPoint posp(pos.x(), pos.y(), pos.z());
0117   GlobalVector momv(mom.x(), mom.y(), mom.z());
0118   FreeTrajectoryState muonFTS(posp, momv, (TrackCharge)muonTrack->charge(), theService->magneticField().product());
0119 
0120   // create a collection on TimeMeasurements for the track
0121   std::vector<TimeMeasurement> tms;
0122   for (const auto& rechit : segments) {
0123     // Create the ChamberId
0124     DetId id = rechit->geographicalId();
0125     DTChamberId chamberId(id.rawId());
0126     int station = chamberId.station();
0127     if (debug)
0128       std::cout << "Matched DT segment in station " << station << std::endl;
0129 
0130     // use only segments with both phi and theta projections present (optional)
0131     bool bothProjections = ((rechit->hasPhi()) && (rechit->hasZed()));
0132     if (requireBothProjections_ && !bothProjections)
0133       continue;
0134 
0135     // loop over (theta, phi) segments
0136     for (int phi = 0; phi < 2; phi++) {
0137       if (dropTheta_ && !phi)
0138         continue;
0139 
0140       const DTRecSegment2D* segm;
0141       if (phi) {
0142         segm = dynamic_cast<const DTRecSegment2D*>(rechit->phiSegment());
0143         if (debug)
0144           std::cout << " *** Segment t0: " << segm->t0() << std::endl;
0145       } else
0146         segm = dynamic_cast<const DTRecSegment2D*>(rechit->zSegment());
0147 
0148       if (segm == nullptr)
0149         continue;
0150       if (segm->specificRecHits().empty())
0151         continue;
0152 
0153       const GeomDet* geomDet = theTrackingGeometry->idToDet(segm->geographicalId());
0154       const std::vector<DTRecHit1D>& hits1d(segm->specificRecHits());
0155 
0156       // store all the hits from the segment
0157       for (const auto& hiti : hits1d) {
0158         const GeomDet* dtcell = theTrackingGeometry->idToDet(hiti.geographicalId());
0159         TimeMeasurement thisHit;
0160 
0161         std::pair<TrajectoryStateOnSurface, double> tsos;
0162         tsos = propag->propagateWithPath(muonFTS, dtcell->surface());
0163 
0164         double dist;
0165         double dist_straight = dtcell->toGlobal(hiti.localPosition()).mag();
0166         if (tsos.first.isValid()) {
0167           dist = tsos.second + posp.mag();
0168           //      std::cout << "Propagate distance: " << dist << " ( innermost: " << posp.mag() << ")" << std::endl;
0169         } else {
0170           dist = dist_straight;
0171           //      std::cout << "Geom distance: " << dist << std::endl;
0172         }
0173 
0174         thisHit.driftCell = hiti.geographicalId();
0175         if (hiti.lrSide() == DTEnums::Left)
0176           thisHit.isLeft = true;
0177         else
0178           thisHit.isLeft = false;
0179         thisHit.isPhi = phi;
0180         thisHit.posInLayer = geomDet->toLocal(dtcell->toGlobal(hiti.localPosition())).x();
0181         thisHit.distIP = dist;
0182         thisHit.station = station;
0183         if (useSegmentT0_ && segm->ist0Valid())
0184           thisHit.timeCorr = segm->t0();
0185         else
0186           thisHit.timeCorr = 0.;
0187         thisHit.timeCorr += theTimeOffset_;
0188 
0189         // signal propagation along the wire correction for unmached theta or phi segment hits
0190         if (doWireCorr_ && !bothProjections && tsos.first.isValid()) {
0191           const DTLayer* layer = theDTGeom.layer(hiti.wireId());
0192           float propgL = layer->toLocal(tsos.first.globalPosition()).y();
0193           float wirePropCorr = propgL / 24.4 * 0.00543;  // signal propagation speed along the wire
0194           if (thisHit.isLeft)
0195             wirePropCorr = -wirePropCorr;
0196           thisHit.posInLayer += wirePropCorr;
0197           const DTSuperLayer* sl = layer->superLayer();
0198           float tofCorr = sl->surface().position().mag() - tsos.first.globalPosition().mag();
0199           tofCorr = (tofCorr / 29.979) * 0.00543;
0200           if (thisHit.isLeft)
0201             tofCorr = -tofCorr;
0202           thisHit.posInLayer += tofCorr;
0203         } else {
0204           // non straight-line path correction
0205           float slCorr = (dist_straight - dist) / 29.979 * 0.00543;
0206           if (thisHit.isLeft)
0207             slCorr = -slCorr;
0208           thisHit.posInLayer += slCorr;
0209         }
0210 
0211         if (debug)
0212           std::cout << " dist: " << dist << "  t0: " << thisHit.posInLayer << std::endl;
0213 
0214         tms.push_back(thisHit);
0215       }
0216     }  // phi = (0,1)
0217   }    // rechit
0218 
0219   bool modified = false;
0220   std::vector<double> dstnc, local_t0, hitWeightTimeVtx, hitWeightInvbeta, left;
0221   double totalWeightInvbeta = 0;
0222   double totalWeightTimeVtx = 0;
0223 
0224   // Now loop over the measurements, calculate 1/beta and time at vertex and cut away outliers
0225   do {
0226     modified = false;
0227     dstnc.clear();
0228     local_t0.clear();
0229     hitWeightTimeVtx.clear();
0230     hitWeightInvbeta.clear();
0231     left.clear();
0232 
0233     std::vector<int> hit_idx;
0234     totalWeightInvbeta = 0;
0235     totalWeightTimeVtx = 0;
0236 
0237     // Rebuild segments
0238     for (int sta = 1; sta < 5; sta++)
0239       for (int phi = 0; phi < 2; phi++) {
0240         std::vector<TimeMeasurement> seg;
0241         std::vector<int> seg_idx;
0242         int tmpos = 0;
0243         for (auto& tm : tms) {
0244           if ((tm.station == sta) && (tm.isPhi == phi)) {
0245             seg.push_back(tm);
0246             seg_idx.push_back(tmpos);
0247           }
0248           tmpos++;
0249         }
0250 
0251         unsigned int segsize = seg.size();
0252         if (segsize < theHitsMin_)
0253           continue;
0254 
0255         double a = 0, b = 0;
0256         std::vector<double> hitxl, hitxr, hityl, hityr;
0257 
0258         for (auto& tm : seg) {
0259           DetId id = tm.driftCell;
0260           const GeomDet* dtcell = theTrackingGeometry->idToDet(id);
0261           DTChamberId chamberId(id.rawId());
0262           const GeomDet* dtcham = theTrackingGeometry->idToDet(chamberId);
0263 
0264           double celly = dtcham->toLocal(dtcell->position()).z();
0265 
0266           if (tm.isLeft) {
0267             hitxl.push_back(celly);
0268             hityl.push_back(tm.posInLayer);
0269           } else {
0270             hitxr.push_back(celly);
0271             hityr.push_back(tm.posInLayer);
0272           }
0273         }
0274 
0275         if (!fitT0(a, b, hitxl, hityl, hitxr, hityr)) {
0276           if (debug)
0277             std::cout << "     t0 = zero, Left hits: " << hitxl.size() << " Right hits: " << hitxr.size() << std::endl;
0278           continue;
0279         }
0280 
0281         // a segment must have at least one left and one right hit
0282         if ((hitxl.empty()) || (hityl.empty()))
0283           continue;
0284 
0285         int segidx = 0;
0286         for (const auto& tm : seg) {
0287           DetId id = tm.driftCell;
0288           const GeomDet* dtcell = theTrackingGeometry->idToDet(id);
0289           DTChamberId chamberId(id.rawId());
0290           const GeomDet* dtcham = theTrackingGeometry->idToDet(chamberId);
0291 
0292           double layerZ = dtcham->toLocal(dtcell->position()).z();
0293           double segmLocalPos = b + layerZ * a;
0294           double hitLocalPos = tm.posInLayer;
0295           int hitSide = -tm.isLeft * 2 + 1;
0296           double t0_segm = (-(hitSide * segmLocalPos) + (hitSide * hitLocalPos)) / 0.00543 + tm.timeCorr;
0297 
0298           if (debug)
0299             std::cout << "   Segm hit.  dstnc: " << tm.distIP << "   t0: " << t0_segm << std::endl;
0300 
0301           dstnc.push_back(tm.distIP);
0302           local_t0.push_back(t0_segm);
0303           left.push_back(hitSide);
0304           hitWeightInvbeta.push_back(((double)seg.size() - 2.) * tm.distIP * tm.distIP /
0305                                      ((double)seg.size() * 30. * 30. * theError_ * theError_));
0306           hitWeightTimeVtx.push_back(((double)seg.size() - 2.) / ((double)seg.size() * theError_ * theError_));
0307           hit_idx.push_back(seg_idx.at(segidx));
0308           segidx++;
0309           totalWeightInvbeta += ((double)seg.size() - 2.) * tm.distIP * tm.distIP /
0310                                 ((double)seg.size() * 30. * 30. * theError_ * theError_);
0311           totalWeightTimeVtx += ((double)seg.size() - 2.) / ((double)seg.size() * theError_ * theError_);
0312         }
0313       }
0314 
0315     if (totalWeightInvbeta == 0)
0316       break;
0317 
0318     // calculate the value and error of 1/beta and timeVtx from the complete set of 1D hits
0319     if (debug)
0320       std::cout << " Points for global fit: " << dstnc.size() << std::endl;
0321 
0322     double invbeta = 0, invbetaErr = 0;
0323     double timeVtx = 0, timeVtxErr = 0;
0324 
0325     for (unsigned int i = 0; i < dstnc.size(); i++) {
0326       invbeta += (1. + local_t0.at(i) / dstnc.at(i) * 30.) * hitWeightInvbeta.at(i) / totalWeightInvbeta;
0327       timeVtx += local_t0.at(i) * hitWeightTimeVtx.at(i) / totalWeightTimeVtx;
0328     }
0329 
0330     double chimax = 0.;
0331     std::vector<TimeMeasurement>::iterator tmmax;
0332 
0333     // Calculate the inv beta and time at vertex dispersion
0334     double diff_ibeta, diff_tvtx;
0335     for (unsigned int i = 0; i < dstnc.size(); i++) {
0336       diff_ibeta = (1. + local_t0.at(i) / dstnc.at(i) * 30.) - invbeta;
0337       diff_ibeta = diff_ibeta * diff_ibeta * hitWeightInvbeta.at(i);
0338       diff_tvtx = local_t0.at(i) - timeVtx;
0339       diff_tvtx = diff_tvtx * diff_tvtx * hitWeightTimeVtx.at(i);
0340       invbetaErr += diff_ibeta;
0341       timeVtxErr += diff_tvtx;
0342 
0343       // decide if we cut away time at vertex outliers or inverse beta outliers
0344       // currently not configurable.
0345       if (diff_tvtx > chimax) {
0346         tmmax = tms.begin() + hit_idx.at(i);
0347         chimax = diff_tvtx;
0348       }
0349     }
0350 
0351     // cut away the outliers
0352     if (chimax > thePruneCut_) {
0353       tms.erase(tmmax);
0354       modified = true;
0355     }
0356 
0357     if (debug) {
0358       double cf = 1. / (dstnc.size() - 1);
0359       invbetaErr = sqrt(invbetaErr / totalWeightInvbeta * cf);
0360       timeVtxErr = sqrt(timeVtxErr / totalWeightTimeVtx * cf);
0361       std::cout << " Measured 1/beta: " << invbeta << " +/- " << invbetaErr << std::endl;
0362       std::cout << " Measured time: " << timeVtx << " +/- " << timeVtxErr << std::endl;
0363     }
0364 
0365   } while (modified);
0366 
0367   for (unsigned int i = 0; i < dstnc.size(); i++) {
0368     tmSequence.dstnc.push_back(dstnc.at(i));
0369     tmSequence.local_t0.push_back(local_t0.at(i));
0370     tmSequence.weightInvbeta.push_back(hitWeightInvbeta.at(i));
0371     tmSequence.weightTimeVtx.push_back(hitWeightTimeVtx.at(i));
0372   }
0373 
0374   tmSequence.totalWeightInvbeta = totalWeightInvbeta;
0375   tmSequence.totalWeightTimeVtx = totalWeightTimeVtx;
0376 }
0377 
0378 // ------------ method called to produce the data  ------------
0379 void DTTimingExtractor::fillTiming(TimeMeasurementSequence& tmSequence,
0380                                    reco::TrackRef muonTrack,
0381                                    const edm::Event& iEvent,
0382                                    const edm::EventSetup& iSetup) {
0383   // get the DT segments that were used to construct the muon
0384   std::vector<const DTRecSegment4D*> range = theMatcher->matchDT(*muonTrack, iEvent);
0385 
0386   if (debug)
0387     std::cout << " The muon track matches " << range.size() << " segments." << std::endl;
0388   fillTiming(tmSequence, range, muonTrack, iEvent, iSetup);
0389 }
0390 
0391 double DTTimingExtractor::fitT0(double& a,
0392                                 double& b,
0393                                 const std::vector<double>& xl,
0394                                 const std::vector<double>& yl,
0395                                 const std::vector<double>& xr,
0396                                 const std::vector<double>& yr) {
0397   double sx = 0, sy = 0, sxy = 0, sxx = 0, ssx = 0, ssy = 0, s = 0, ss = 0;
0398 
0399   for (unsigned int i = 0; i < xl.size(); i++) {
0400     sx += xl[i];
0401     sy += yl[i];
0402     sxy += xl[i] * yl[i];
0403     sxx += xl[i] * xl[i];
0404     s++;
0405     ssx += xl[i];
0406     ssy += yl[i];
0407     ss++;
0408   }
0409 
0410   for (unsigned int i = 0; i < xr.size(); i++) {
0411     sx += xr[i];
0412     sy += yr[i];
0413     sxy += xr[i] * yr[i];
0414     sxx += xr[i] * xr[i];
0415     s++;
0416     ssx -= xr[i];
0417     ssy -= yr[i];
0418     ss--;
0419   }
0420 
0421   double delta = ss * ss * sxx + s * sx * sx + s * ssx * ssx - s * s * sxx - 2 * ss * sx * ssx;
0422 
0423   double t0_corr = 0.;
0424 
0425   if (delta) {
0426     a = (ssy * s * ssx + sxy * ss * ss + sy * sx * s - sy * ss * ssx - ssy * sx * ss - sxy * s * s) / delta;
0427     b = (ssx * sy * ssx + sxx * ssy * ss + sx * sxy * s - sxx * sy * s - ssx * sxy * ss - sx * ssy * ssx) / delta;
0428     t0_corr = (ssx * s * sxy + sxx * ss * sy + sx * sx * ssy - sxx * s * ssy - sx * ss * sxy - ssx * sx * sy) / delta;
0429   }
0430 
0431   // convert drift distance to time
0432   t0_corr /= -0.00543;
0433 
0434   return t0_corr;
0435 }
0436 
0437 //define this as a plug-in
0438 //DEFINE_FWK_MODULE(DTTimingExtractor);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!