Navigation/src/DirectMuonNavigation.cc

0001 #include "RecoMuon/Navigation/interface/DirectMuonNavigation.h"
0002
0003 /** \file DirectMuonNavigation
0004  *
0005  *  \author Chang Liu  -  Purdue University
0006  */
0007
0008 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0009 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0010 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
0011 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
0012 #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
0013 #include "FWCore/Framework/interface/ESHandle.h"
0014
0015 #include <algorithm>
0016
0017 using namespace std;
0018
0019 DirectMuonNavigation::DirectMuonNavigation(const edm::ESHandle<MuonDetLayerGeometry>& muonLayout)
0020     : theMuonDetLayerGeometry(muonLayout), epsilon_(100.), theEndcapFlag(true), theBarrelFlag(true) {}
0021
0022 DirectMuonNavigation::DirectMuonNavigation(const edm::ESHandle<MuonDetLayerGeometry>& muonLayout,
0023                                            const edm::ParameterSet& par)
0024     : theMuonDetLayerGeometry(muonLayout),
0025       epsilon_(100.),
0026       theEndcapFlag(par.getParameter<bool>("Endcap")),
0027       theBarrelFlag(par.getParameter<bool>("Barrel")) {}
0028
0029 /* return compatible layers for given trajectory state  */
0030 vector<const DetLayer*> DirectMuonNavigation::compatibleLayers(const FreeTrajectoryState& fts,
0031                                                                PropagationDirection dir) const {
0032   float z0 = fts.position().z();
0033   float zm = fts.momentum().z();
0034
0035   bool inOut = outward(fts);
0036
0037   vector<const DetLayer*> output;
0038
0039   // check direction and position of FTS to get a correct order of DetLayers
0040
0041   if (inOut) {
0042     if ((zm * z0) >= 0) {
0043       if (theBarrelFlag)
0044         inOutBarrel(fts, output);
0045       if (theEndcapFlag) {
0046         if (z0 >= 0)
0047           inOutForward(fts, output);
0048         else
0049           inOutBackward(fts, output);
0050       }
0051     } else {
0052       if (theEndcapFlag) {
0053         if (z0 >= 0)
0054           outInForward(fts, output);
0055         else
0056           outInBackward(fts, output);
0057       }
0058       if (theBarrelFlag)
0059         inOutBarrel(fts, output);
0060     }
0061   } else {
0062     if ((zm * z0) >= 0) {
0063       if (theBarrelFlag)
0064         outInBarrel(fts, output);
0065       if (theEndcapFlag) {
0066         if (z0 >= 0)
0067           inOutForward(fts, output);
0068         else
0069           inOutBackward(fts, output);
0070       }
0071     } else {
0072       if (theEndcapFlag) {
0073         if (z0 >= 0)
0074           outInForward(fts, output);
0075         else
0076           outInBackward(fts, output);
0077       }
0078       if (theBarrelFlag)
0079         outInBarrel(fts, output);
0080     }
0081   }
0082
0083   if (dir == oppositeToMomentum)
0084     std::reverse(output.begin(), output.end());
0085
0086   return output;
0087 }
0088
0089 /*
0090 return compatible endcap layers on BOTH ends;
0091 used for beam-halo muons
0092 */
0093 vector<const DetLayer*> DirectMuonNavigation::compatibleEndcapLayers(const FreeTrajectoryState& fts,
0094                                                                      PropagationDirection dir) const {
0095   float zm = fts.momentum().z();
0096
0097   vector<const DetLayer*> output;
0098
0099   // collect all endcap layers on 2 sides
0100   outInBackward(fts, output);
0101   inOutForward(fts, output);
0102
0103   // check direction FTS to get a correct order of DetLayers
0104   if ((zm > 0 && dir == oppositeToMomentum) || (zm < 0 && dir == alongMomentum))
0105     std::reverse(output.begin(), output.end());
0106
0107   return output;
0108 }
0109
0110 void DirectMuonNavigation::inOutBarrel(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
0111   bool cont = false;
0112   const vector<const DetLayer*>& barrel = theMuonDetLayerGeometry->allBarrelLayers();
0113
0114   for (vector<const DetLayer*>::const_iterator iter_B = barrel.begin(); iter_B != barrel.end(); iter_B++) {
0115     if (cont)
0116       output.push_back((*iter_B));
0117     else if (checkCompatible(fts, dynamic_cast<const BarrelDetLayer*>(*iter_B))) {
0118       output.push_back((*iter_B));
0119       cont = true;
0120     }
0121   }
0122 }
0123
0124 void DirectMuonNavigation::outInBarrel(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
0125   // default barrel layers are in out, reverse order
0126   const vector<const DetLayer*>& barrel = theMuonDetLayerGeometry->allBarrelLayers();
0127
0128   bool cont = false;
0129   vector<const DetLayer*>::const_iterator rbegin = barrel.end();
0130   rbegin--;
0131   vector<const DetLayer*>::const_iterator rend = barrel.begin();
0132   rend--;
0133
0134   for (vector<const DetLayer*>::const_iterator iter_B = rbegin; iter_B != rend; iter_B--) {
0135     if (cont)
0136       output.push_back((*iter_B));
0137     else if (checkCompatible(fts, dynamic_cast<const BarrelDetLayer*>(*iter_B))) {
0138       output.push_back((*iter_B));
0139       cont = true;
0140     }
0141   }
0142 }
0143
0144 void DirectMuonNavigation::inOutForward(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
0145   const vector<const DetLayer*>& forward = theMuonDetLayerGeometry->allForwardLayers();
0146   bool cont = false;
0147   for (vector<const DetLayer*>::const_iterator iter_E = forward.begin(); iter_E != forward.end(); iter_E++) {
0148     if (cont)
0149       output.push_back((*iter_E));
0150     else if (checkCompatible(fts, dynamic_cast<const ForwardDetLayer*>(*iter_E))) {
0151       output.push_back((*iter_E));
0152       cont = true;
0153     }
0154   }
0155 }
0156
0157 void DirectMuonNavigation::outInForward(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
0158   // default forward layers are in out, reverse order
0159
0160   bool cont = false;
0161   const vector<const DetLayer*>& forward = theMuonDetLayerGeometry->allForwardLayers();
0162   vector<const DetLayer*>::const_iterator rbegin = forward.end();
0163   rbegin--;
0164   vector<const DetLayer*>::const_iterator rend = forward.begin();
0165   rend--;
0166   for (vector<const DetLayer*>::const_iterator iter_E = rbegin; iter_E != rend; iter_E--) {
0167     if (cont)
0168       output.push_back((*iter_E));
0169     else if (checkCompatible(fts, dynamic_cast<const ForwardDetLayer*>(*iter_E))) {
0170       output.push_back((*iter_E));
0171       cont = true;
0172     }
0173   }
0174 }
0175
0176 void DirectMuonNavigation::inOutBackward(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
0177   bool cont = false;
0178   const vector<const DetLayer*>& backward = theMuonDetLayerGeometry->allBackwardLayers();
0179
0180   for (vector<const DetLayer*>::const_iterator iter_E = backward.begin(); iter_E != backward.end(); iter_E++) {
0181     if (cont)
0182       output.push_back((*iter_E));
0183     else if (checkCompatible(fts, dynamic_cast<const ForwardDetLayer*>(*iter_E))) {
0184       output.push_back((*iter_E));
0185       cont = true;
0186     }
0187   }
0188 }
0189
0190 void DirectMuonNavigation::outInBackward(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
0191   bool cont = false;
0192   const vector<const DetLayer*>& backward = theMuonDetLayerGeometry->allBackwardLayers();
0193
0194   vector<const DetLayer*>::const_iterator rbegin = backward.end();
0195   rbegin--;
0196   vector<const DetLayer*>::const_iterator rend = backward.begin();
0197   rend--;
0198   for (vector<const DetLayer*>::const_iterator iter_E = rbegin; iter_E != rend; iter_E--) {
0199     if (cont)
0200       output.push_back((*iter_E));
0201     else if (checkCompatible(fts, dynamic_cast<const ForwardDetLayer*>(*iter_E))) {
0202       output.push_back((*iter_E));
0203       cont = true;
0204     }
0205   }
0206 }
0207
0208 bool DirectMuonNavigation::checkCompatible(const FreeTrajectoryState& fts, const BarrelDetLayer* dl) const {
0209   float z0 = fts.position().z();
0210   float r0 = fts.position().perp();
0211   float zm = fts.momentum().z();
0212   float rm = fts.momentum().perp();
0213   float slope = zm / rm;
0214   if (!outward(fts))
0215     slope = -slope;
0216   const BoundCylinder& bc = dl->specificSurface();
0217   float radius = bc.radius();
0218   float length = bc.bounds().length() / 2.;
0219
0220   float z1 = slope * (radius - r0) + z0;
0221   return (fabs(z1) <= fabs(length) + epsilon_);
0222 }
0223
0224 bool DirectMuonNavigation::checkCompatible(const FreeTrajectoryState& fts, const ForwardDetLayer* dl) const {
0225   float z0 = fts.position().z();
0226   float r0 = fts.position().perp();
0227   float zm = fts.momentum().z();
0228   float rm = fts.momentum().perp();
0229   float slope = rm / zm;
0230
0231   if (!outward(fts))
0232     slope = -slope;
0233
0234   const BoundDisk& bd = dl->specificSurface();
0235
0236   float outRadius = bd.outerRadius();
0237   float inRadius = bd.innerRadius();
0238   float z = bd.position().z();
0239
0240   float r1 = slope * (z - z0) + r0;
0241   return (r1 >= inRadius - epsilon_ && r1 <= outRadius + epsilon_);
0242 }
0243
0244 bool DirectMuonNavigation::outward(const FreeTrajectoryState& fts) const {
0245   //  return (fts.position().basicVector().dot(fts.momentum().basicVector())>0);
0246
0247   float x0 = fts.position().x();
0248   float y0 = fts.position().y();
0249
0250   float xm = fts.momentum().x();
0251   float ym = fts.momentum().y();
0252
0253   return ((x0 * xm + y0 * ym) > 0);
0254 }