CMS 3D CMS Logo

PropagateToMuon.cc
Go to the documentation of this file.
2 
3 #include <cmath>
4 
6 
14 
16 
18  useSimpleGeometry_(iConfig.getParameter<bool>("useSimpleGeometry")),
19  useMB2_(iConfig.existsAs<bool>("useStation2") ? iConfig.getParameter<bool>("useStation2") : true),
20  fallbackToME1_(iConfig.existsAs<bool>("fallbackToME1") ? iConfig.getParameter<bool>("fallbackToME1") : false),
21  whichTrack_(None), whichState_(AtVertex),
22  cosmicPropagation_(iConfig.existsAs<bool>("cosmicPropagationHypothesis") ? iConfig.getParameter<bool>("cosmicPropagationHypothesis") : false),
23  useMB2InOverlap_(iConfig.existsAs<bool>("useMB2InOverlap") ? iConfig.getParameter<bool>("useMB2InOverlap") : false)
24 
25 {
26  std::string whichTrack = iConfig.getParameter<std::string>("useTrack");
27  if (whichTrack == "none") { whichTrack_ = None; }
28  else if (whichTrack == "tracker") { whichTrack_ = TrackerTk; }
29  else if (whichTrack == "muon") { whichTrack_ = MuonTk; }
30  else if (whichTrack == "global") { whichTrack_ = GlobalTk; }
31  else throw cms::Exception("Configuration") << "Parameter 'useTrack' must be 'none', 'tracker', 'muon', 'global'\n";
32  if (whichTrack_ != None) {
33  std::string whichState = iConfig.getParameter<std::string>("useState");
34  if (whichState == "atVertex") { whichState_ = AtVertex; }
35  else if (whichState == "innermost") { whichState_ = Innermost; }
36  else if (whichState == "outermost") { whichState_ = Outermost; }
37  else throw cms::Exception("Configuration") << "Parameter 'useState' must be 'atVertex', 'innermost', 'outermost'\n";
38  }
40  throw cms::Exception("Configuration") << "When using 'cosmicPropagationHypothesis' useTrack must not be 'none', and the state must not be 'atVertex'\n";
41  }
42 }
43 
45 
46 void
47 PropagateToMuon::init(const edm::EventSetup & iSetup) {
48  iSetup.get<IdealMagneticFieldRecord>().get(magfield_);
49  iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAlong", propagator_);
50  iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorOpposite", propagatorOpposite_);
51  iSetup.get<TrackingComponentsRecord>().get("SteppingHelixPropagatorAny", propagatorAny_);
53 
54  // Get the barrel cylinder
55  const DetLayer * dtLay = muonGeometry_->allDTLayers()[useMB2_ ? 1 : 0];
56  barrelCylinder_ = dynamic_cast<const BoundCylinder *>(&dtLay->surface());
57  barrelHalfLength_ = barrelCylinder_->bounds().length()/2;;
58  //std::cout << "L1MuonMatcher: barrel radius = " << barrelCylinder_->radius() << ", half length = " << barrelHalfLength_ << std::endl;
59 
60  // 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
61  for (size_t i = 0; i <= (useMB2_ ? 2 : 1); ++i) {
62  endcapDiskPos_[i] = dynamic_cast<const BoundDisk *>(& muonGeometry_->forwardCSCLayers()[i]->surface());
63  endcapDiskNeg_[i] = dynamic_cast<const BoundDisk *>(& muonGeometry_->backwardCSCLayers()[i]->surface());
65  //std::cout << "L1MuonMatcher: endcap " << i << " Z = " << endcapDiskPos_[i]->position().z() << ", radii = " << endcapRadii_[i].first << "," << endcapRadii_[i].second << std::endl;
66  }
67 
69  barrelHalfLength_ = endcapDiskPos_[2]->position().z();
70 
71 }
72 
76  if (whichTrack_ != None) {
77  const reco::RecoCandidate *rc = dynamic_cast<const reco::RecoCandidate *>(&reco);
78  if (rc == nullptr) throw cms::Exception("Invalid Data") << "Input object is not a RecoCandidate.\n";
79  reco::TrackRef tk;
80  switch (whichTrack_) {
81  case TrackerTk: tk = rc->track(); break;
82  case MuonTk : tk = rc->standAloneMuon(); break;
83  case GlobalTk : tk = rc->combinedMuon(); break;
84  default: break; // just to make gcc happy
85  }
86  if (tk.isNull()) {
87  ret = FreeTrajectoryState();
88  } else {
89  ret = startingState(*tk);
90  }
91  } else {
92  ret = FreeTrajectoryState( GlobalPoint( reco.vx(), reco.vy(), reco.vz()),
93  GlobalVector(reco.px(), reco.py(), reco.pz()),
94  reco.charge(),
95  magfield_.product());
96  }
97  return ret;
98 }
99 
102  if (!magfield_.isValid()) throw cms::Exception("NotInitialized") << "PropagateToMuon: You must call init(const edm::EventSetup &iSetup) before using this object.\n";
103  WhichState state = whichState_;
104  if (cosmicPropagation_) {
105  if (whichState_ == Innermost) {
106  state = tk.innerPosition().Mag2() <= tk.outerPosition().Mag2() ? Innermost : Outermost;
107  } else if (whichState_ == Outermost) {
108  state = tk.innerPosition().Mag2() <= tk.outerPosition().Mag2() ? Outermost : Innermost;
109  }
110  }
111  switch (state) {
114 
115  case AtVertex:
116  default:
118  }
119 
120 }
121 
124  if (!magfield_.isValid()) throw cms::Exception("NotInitialized") << "PropagateToMuon: You must call init(const edm::EventSetup &iSetup) before using this object.\n";
125  if (tk.noVertex()) throw cms::Exception("UnsupportedOperation") << "I can't propagate simtracks without a vertex, I don't know where to start from.\n";
126  const math::XYZTLorentzVectorD & vtx = (vtxs)[tk.vertIndex()].position();
127  return FreeTrajectoryState( GlobalPoint(vtx.X(), vtx.Y(), vtx.Z()),
128  GlobalVector(tk.momentum().X(), tk.momentum().Y(), tk.momentum().Z()),
129  int(tk.charge()),
130  magfield_.product());
131 }
132 
133 
134 
137  if (!magfield_.isValid() || barrelCylinder_ == nullptr) {
138  throw cms::Exception("NotInitialized") << "PropagateToMuon: You must call init(const edm::EventSetup &iSetup) before using this object.\n";
139  }
140 
142  if (start.momentum().mag() == 0) return final;
143  double eta = start.momentum().eta();
144 
145  const Propagator * propagatorBarrel = &*propagator_;
146  const Propagator * propagatorEndcaps = &*propagator_;
147  if (whichState_ != AtVertex) {
148  if (start.position().perp() > barrelCylinder_->radius()) propagatorBarrel = &*propagatorOpposite_;
149  if (fabs(start.position().z()) > endcapDiskPos_[useMB2_?2:1]->position().z()) propagatorEndcaps = &*propagatorOpposite_;
150  }
151  if (cosmicPropagation_) {
152  if (start.momentum().dot(GlobalVector(start.position().x(), start.position().y(), start.position().z())) < 0) {
153  // must flip the propagations
154  propagatorBarrel = (propagatorBarrel == &*propagator_ ? &*propagatorOpposite_ : &*propagator_);
155  propagatorEndcaps = (propagatorEndcaps == &*propagator_ ? &*propagatorOpposite_ : &*propagator_);
156  }
157  }
158 
159  TrajectoryStateOnSurface tsos = propagatorBarrel->propagate(start, *barrelCylinder_);
160  if (tsos.isValid()) {
161  if (useSimpleGeometry_) {
162  //std::cout << " propagated to barrel, z = " << tsos.globalPosition().z() << ", bound = " << barrelHalfLength_ << std::endl;
163  if (fabs(tsos.globalPosition().z()) <= barrelHalfLength_) final = tsos;
164  } else {
165  final = getBestDet(tsos, muonGeometry_->allDTLayers()[1]);
166  }
167  }
168  if (!final.isValid()) {
169  for (int ie = (useMB2_ ? 2 : 1); ie >= 0; --ie) {
170  tsos = propagatorEndcaps->propagate(start, (eta > 0 ? *endcapDiskPos_[ie] : *endcapDiskNeg_[ie]));
171  if (tsos.isValid()) {
172  if (useSimpleGeometry_) {
173  float rho = tsos.globalPosition().perp();
174  //std::cout << " propagated to endcap " << ie << ", rho = " << rho << ", bounds [ " << endcapRadii_[ie].first << ", " << endcapRadii_[ie].second << "]" << std::endl;
175  if ((rho >= endcapRadii_[ie].first) && (rho <= endcapRadii_[ie].second)) final = tsos;
176  } else {
177  final = getBestDet(tsos, (eta > 0 ? muonGeometry_->forwardCSCLayers()[ie] : muonGeometry_->backwardCSCLayers()[ie]));
178  }
179  } //else //std::cout << " failed to propagated to endcap " << ie << std::endl;
180  if (final.isValid()) break;
181  if (ie == 2 && !fallbackToME1_) break;
182  }
183  }
184  return final;
185 }
186 
188 PropagateToMuon::getBestDet(const TrajectoryStateOnSurface &tsos, const DetLayer *layer) const {
189  TrajectoryStateOnSurface ret; // start as null
190  Chi2MeasurementEstimator estimator(1e10, 3.); // require compatibility at 3 sigma
191  std::vector<GeometricSearchDet::DetWithState> dets = layer->compatibleDets(tsos, *propagatorAny_, estimator);
192  if (!dets.empty()) {
193  ret = dets.front().second;
194  }
195  return ret;
196 }
Definition: start.py:1
WhichTrack whichTrack_
Labels for input collections.
T getParameter(std::string const &) const
PropagateToMuon(const edm::ParameterSet &iConfig)
TrajectoryStateOnSurface extrapolate(const reco::Track &tk) const
edm::ESHandle< Propagator > propagatorAny_
virtual double pz() const =0
z coordinate of momentum vector
const BoundDisk * endcapDiskNeg_
T perp() const
Definition: PV3DBase.h:72
FreeTrajectoryState startingState(const reco::Candidate &reco) const
Starting state for the propagation.
ROOT::Math::LorentzVector< ROOT::Math::PxPyPzE4D< double > > XYZTLorentzVectorD
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:14
virtual double vx() const =0
x coordinate of vertex position
WhichState whichState_
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
T y() const
Definition: PV3DBase.h:63
GlobalPoint globalPosition() const
const BoundDisk * endcapDiskPos_
PreciseFloatType< T, U >::Type dot(const Vector3DBase< U, FrameTag > &v) const
Definition: Vector3DBase.h:107
bool useMB2_
Propagate to MB2 (default) instead of MB1.
virtual double vy() const =0
y coordinate of vertex position
virtual reco::TrackRef standAloneMuon() const
reference to a stand-alone muon Track
TrajectoryStateOnSurface getBestDet(const TrajectoryStateOnSurface &tsos, const DetLayer *station) const
const math::XYZPoint & outerPosition() const
position of the outermost hit
Definition: Track.h:67
float charge() const
charge
Definition: CoreSimTrack.cc:17
virtual std::vector< DetWithState > compatibleDets(const TrajectoryStateOnSurface &startingState, const Propagator &prop, const MeasurementEstimator &est) const
virtual reco::TrackRef track() const
reference to a Track
edm::ESHandle< MagneticField > magfield_
const BoundCylinder * barrelCylinder_
const math::XYZPoint & innerPosition() const
position of the innermost hit
Definition: Track.h:57
U second(std::pair< T, U > const &p)
virtual double py() const =0
y coordinate of momentum vector
T mag() const
Definition: PV3DBase.h:67
T z() const
Definition: PV3DBase.h:64
bool noVertex() const
Definition: SimTrack.h:31
virtual const BoundSurface & surface() const =0
The surface of the GeometricSearchDet.
edm::ESHandle< Propagator > propagator_
bool isNull() const
Checks for null.
Definition: Ref.h:248
GlobalVector momentum() const
int vertIndex() const
index of the vertex in the Event container (-1 if no vertex)
Definition: SimTrack.h:30
bool cosmicPropagation_
for cosmics, some things change: the along-opposite is not in-out, nor the innermost/outermost states...
GlobalPoint position() const
const std::vector< const DetLayer * > & allDTLayers() const
return the DT DetLayers (barrel), inside-out
const std::vector< const DetLayer * > & forwardCSCLayers() const
return the forward (+Z) CSC DetLayers, inside-out
std::vector< SimVertex > SimVertexContainer
void init(const edm::EventSetup &iSetup)
TrajectoryStateOnSurface propagate(STA const &state, SUR const &surface) const
Definition: Propagator.h:53
T eta() const
Definition: PV3DBase.h:76
virtual int charge() const =0
electric charge
fixed size matrix
bool useSimpleGeometry_
Labels for input collections.
const math::XYZTLorentzVectorD & momentum() const
Definition: CoreSimTrack.h:19
static int position[264][3]
Definition: ReadPGInfo.cc:509
T get() const
Definition: EventSetup.h:71
edm::ESHandle< MuonDetLayerGeometry > muonGeometry_
bool fallbackToME1_
Fallback to ME1 if propagation to ME2 fails.
double barrelHalfLength_
const std::vector< const DetLayer * > & backwardCSCLayers() const
return the backward (-Z) CSC DetLayers, inside-out
virtual double vz() const =0
z coordinate of vertex position
FreeTrajectoryState innerFreeState(const reco::Track &tk, const MagneticField *field, bool withErr=true)
virtual double px() const =0
x coordinate of momentum vector
bool isValid() const
Definition: ESHandle.h:44
T x() const
Definition: PV3DBase.h:62
FreeTrajectoryState outerFreeState(const reco::Track &tk, const MagneticField *field, bool withErr=true)
FreeTrajectoryState initialFreeState(const reco::Track &tk, const MagneticField *field, bool withErr=true)
edm::ESHandle< Propagator > propagatorOpposite_
T const * product() const
Definition: ESHandle.h:86
Global3DVector GlobalVector
Definition: GlobalVector.h:10
std::pair< float, float > endcapRadii_
virtual reco::TrackRef combinedMuon() const
reference to a stand-alone muon Track