Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:22:48 

0001 #include "MuonAnalysis/MuonAssociators/interface/PropagateToMuon.h"
0002 
0003 #include <cmath>
0004 
0005 #include "DataFormats/RecoCandidate/interface/RecoCandidate.h"
0006 #include "DataFormats/GeometrySurface/interface/TrapezoidalPlaneBounds.h"
0007 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0008 #include "RecoMuon/Records/interface/MuonRecoGeometryRecord.h"
0009 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0010 #include "TrackingTools/KalmanUpdators/interface/Chi2MeasurementEstimator.h"
0011 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0012 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0013 
0014 PropagateToMuon::PropagateToMuon(edm::ESHandle<MagneticField> magfield,
0015                                  edm::ESHandle<Propagator> propagator,
0016                                  edm::ESHandle<Propagator> propagatorAny,
0017                                  edm::ESHandle<Propagator> propagatorOpposite,
0018                                  edm::ESHandle<MuonDetLayerGeometry> muonGeometry,
0019                                  bool useSimpleGeometry,
0020                                  bool useMB2,
0021                                  bool fallbackToME1,
0022                                  WhichTrack whichTrack,
0023                                  WhichState whichState,
0024                                  bool cosmicPropagation,
0025                                  bool useMB2InOverlap)
0026     : magfield_(magfield),
0027       propagator_(propagator),
0028       propagatorAny_(propagatorAny),
0029       propagatorOpposite_(propagatorOpposite),
0030       muonGeometry_(muonGeometry),
0031       useSimpleGeometry_(useSimpleGeometry),
0032       useMB2_(useMB2),
0033       fallbackToME1_(fallbackToME1),
0034       whichTrack_(whichTrack),
0035       whichState_(whichState),
0036       cosmicPropagation_(cosmicPropagation),
0037       useMB2InOverlap_(useMB2InOverlap) {
0038   // Get the barrel cylinder
0039   const DetLayer *dtLay = muonGeometry_->allDTLayers()[useMB2_ ? 1 : 0];
0040   barrelCylinder_ = dynamic_cast<const BoundCylinder *>(&dtLay->surface());
0041   barrelHalfLength_ = barrelCylinder_->bounds().length() / 2;
0042   ;
0043   //std::cout << "L1MuonMatcher: barrel radius = " << barrelCylinder_->radius() << ", half length = " << barrelHalfLength_ << std::endl;
0044 
0045   // Get the endcap disks. Note that ME1 has two disks (ME1/1 and ME2/1-ME3/2-ME4/1), so there's one more index
0046   for (size_t i = 0; i <= (useMB2_ ? 2 : 1); ++i) {
0047     endcapDiskPos_[i] = dynamic_cast<const BoundDisk *>(&muonGeometry_->forwardCSCLayers()[i]->surface());
0048     endcapDiskNeg_[i] = dynamic_cast<const BoundDisk *>(&muonGeometry_->backwardCSCLayers()[i]->surface());
0049     endcapRadii_[i] = std::make_pair(endcapDiskPos_[i]->innerRadius(), endcapDiskPos_[i]->outerRadius());
0050     //std::cout << "L1MuonMatcher: endcap " << i << " Z = " << endcapDiskPos_[i]->position().z() << ", radii = " << endcapRadii_[i].first << "," << endcapRadii_[i].second << std::endl;
0051   }
0052 
0053   if (useMB2_ && useMB2InOverlap_)
0054     barrelHalfLength_ = endcapDiskPos_[2]->position().z();
0055 }
0056 
0057 FreeTrajectoryState PropagateToMuon::startingState(const reco::Candidate &reco) const {
0058   FreeTrajectoryState ret;
0059   if (whichTrack_ != None) {
0060     const reco::RecoCandidate *rc = dynamic_cast<const reco::RecoCandidate *>(&reco);
0061     if (rc == nullptr)
0062       throw cms::Exception("Invalid Data") << "Input object is not a RecoCandidate.\n";
0063     reco::TrackRef tk;
0064     switch (whichTrack_) {
0065       case TrackerTk:
0066         tk = rc->track();
0067         break;
0068       case MuonTk:
0069         tk = rc->standAloneMuon();
0070         break;
0071       case GlobalTk:
0072         tk = rc->combinedMuon();
0073         break;
0074       default:
0075         break;  // just to make gcc happy
0076     }
0077     if (tk.isNull()) {
0078       ret = FreeTrajectoryState();
0079     } else {
0080       ret = startingState(*tk);
0081     }
0082   } else {
0083     ret = FreeTrajectoryState(GlobalPoint(reco.vx(), reco.vy(), reco.vz()),
0084                               GlobalVector(reco.px(), reco.py(), reco.pz()),
0085                               reco.charge(),
0086                               magfield_.product());
0087   }
0088   return ret;
0089 }
0090 
0091 FreeTrajectoryState PropagateToMuon::startingState(const reco::Track &tk) const {
0092   if (!magfield_.isValid())
0093     throw cms::Exception("NotInitialized")
0094         << "PropagateToMuon: You must call init(const edm::EventSetup &iSetup) before using this object.\n";
0095   WhichState state = whichState_;
0096   if (cosmicPropagation_) {
0097     if (whichState_ == Innermost) {
0098       state = tk.innerPosition().Mag2() <= tk.outerPosition().Mag2() ? Innermost : Outermost;
0099     } else if (whichState_ == Outermost) {
0100       state = tk.innerPosition().Mag2() <= tk.outerPosition().Mag2() ? Outermost : Innermost;
0101     }
0102   }
0103   switch (state) {
0104     case Innermost:
0105       return trajectoryStateTransform::innerFreeState(tk, magfield_.product());
0106     case Outermost:
0107       return trajectoryStateTransform::outerFreeState(tk, magfield_.product());
0108 
0109     case AtVertex:
0110     default:
0111       return trajectoryStateTransform::initialFreeState(tk, magfield_.product());
0112   }
0113 }
0114 
0115 FreeTrajectoryState PropagateToMuon::startingState(const SimTrack &tk, const edm::SimVertexContainer &vtxs) const {
0116   if (!magfield_.isValid())
0117     throw cms::Exception("NotInitialized")
0118         << "PropagateToMuon: You must call init(const edm::EventSetup &iSetup) before using this object.\n";
0119   if (tk.noVertex())
0120     throw cms::Exception("UnsupportedOperation")
0121         << "I can't propagate simtracks without a vertex, I don't know where to start from.\n";
0122   const math::XYZTLorentzVectorD &vtx = (vtxs)[tk.vertIndex()].position();
0123   return FreeTrajectoryState(GlobalPoint(vtx.X(), vtx.Y(), vtx.Z()),
0124                              GlobalVector(tk.momentum().X(), tk.momentum().Y(), tk.momentum().Z()),
0125                              int(tk.charge()),
0126                              magfield_.product());
0127 }
0128 
0129 TrajectoryStateOnSurface PropagateToMuon::extrapolate(const FreeTrajectoryState &start) const {
0130   if (!magfield_.isValid() || barrelCylinder_ == nullptr) {
0131     throw cms::Exception("NotInitialized")
0132         << "PropagateToMuon: You must call init(const edm::EventSetup &iSetup) before using this object.\n";
0133   }
0134 
0135   TrajectoryStateOnSurface final;
0136   if (start.momentum().mag() == 0)
0137     return final;
0138   double eta = start.momentum().eta();
0139 
0140   const Propagator *propagatorBarrel = &*propagator_;
0141   const Propagator *propagatorEndcaps = &*propagator_;
0142   if (whichState_ != AtVertex) {
0143     if (start.position().perp() > barrelCylinder_->radius())
0144       propagatorBarrel = &*propagatorOpposite_;
0145     if (fabs(start.position().z()) > endcapDiskPos_[useMB2_ ? 2 : 1]->position().z())
0146       propagatorEndcaps = &*propagatorOpposite_;
0147   }
0148   if (cosmicPropagation_) {
0149     if (start.momentum().dot(GlobalVector(start.position().x(), start.position().y(), start.position().z())) < 0) {
0150       // must flip the propagations
0151       propagatorBarrel = (propagatorBarrel == &*propagator_ ? &*propagatorOpposite_ : &*propagator_);
0152       propagatorEndcaps = (propagatorEndcaps == &*propagator_ ? &*propagatorOpposite_ : &*propagator_);
0153     }
0154   }
0155 
0156   TrajectoryStateOnSurface tsos = propagatorBarrel->propagate(start, *barrelCylinder_);
0157   if (tsos.isValid()) {
0158     if (useSimpleGeometry_) {
0159       //std::cout << "  propagated to barrel, z = " << tsos.globalPosition().z() << ", bound = " << barrelHalfLength_ << std::endl;
0160       if (fabs(tsos.globalPosition().z()) <= barrelHalfLength_)
0161         final = tsos;
0162     } else {
0163       final = getBestDet(tsos, muonGeometry_->allDTLayers()[1]);
0164     }
0165   }
0166   if (!final.isValid()) {
0167     for (int ie = (useMB2_ ? 2 : 1); ie >= 0; --ie) {
0168       tsos = propagatorEndcaps->propagate(start, (eta > 0 ? *endcapDiskPos_[ie] : *endcapDiskNeg_[ie]));
0169       if (tsos.isValid()) {
0170         if (useSimpleGeometry_) {
0171           float rho = tsos.globalPosition().perp();
0172           //std::cout << "  propagated to endcap " << ie << ", rho = " << rho << ", bounds [ " << endcapRadii_[ie].first << ", " << endcapRadii_[ie].second << "]" << std::endl;
0173           if ((rho >= endcapRadii_[ie].first) && (rho <= endcapRadii_[ie].second))
0174             final = tsos;
0175         } else {
0176           final = getBestDet(
0177               tsos, (eta > 0 ? muonGeometry_->forwardCSCLayers()[ie] : muonGeometry_->backwardCSCLayers()[ie]));
0178         }
0179       }  //else //std::cout << "  failed to propagated to endcap " << ie  << std::endl;
0180       if (final.isValid())
0181         break;
0182       if (ie == 2 && !fallbackToME1_)
0183         break;
0184     }
0185   }
0186   return final;
0187 }
0188 
0189 TrajectoryStateOnSurface PropagateToMuon::getBestDet(const TrajectoryStateOnSurface &tsos,
0190                                                      const DetLayer *layer) const {
0191   TrajectoryStateOnSurface ret;  // start as null
0192   // require compatibility at 3 sigma
0193   Chi2MeasurementEstimator estimator(1e10, 3.);
0194   std::vector<GeometricSearchDet::DetWithState> dets = layer->compatibleDets(tsos, *propagatorAny_, estimator);
0195   if (!dets.empty()) {
0196     ret = dets.front().second;
0197   }
0198   return ret;
0199 }
[ 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!