00001 #include "RecoMuon/Navigation/interface/DirectMuonNavigation.h"
00002
00010 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
00011 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
00012 #include "DataFormats/GeometrySurface/interface/BoundCylinder.h"
00013 #include "DataFormats/GeometrySurface/interface/BoundDisk.h"
00014 #include "RecoMuon/DetLayers/interface/MuonDetLayerGeometry.h"
00015 #include "FWCore/Framework/interface/ESHandle.h"
00016
00017 #include <algorithm>
00018
00019 using namespace std;
00020
00021 DirectMuonNavigation::DirectMuonNavigation(const edm::ESHandle<MuonDetLayerGeometry>& muonLayout) : theMuonDetLayerGeometry(muonLayout), epsilon_(100.), theEndcapFlag(true), theBarrelFlag(true) {
00022 }
00023
00024 DirectMuonNavigation::DirectMuonNavigation(const edm::ESHandle<MuonDetLayerGeometry>& muonLayout, const edm::ParameterSet& par) : theMuonDetLayerGeometry(muonLayout), epsilon_(100.), theEndcapFlag(par.getUntrackedParameter<bool>("Endcap", true)), theBarrelFlag(par.getUntrackedParameter<bool>("Barrel",true)) {
00025
00026 }
00027
00028
00029
00030 vector<const DetLayer*>
00031 DirectMuonNavigation::compatibleLayers( const FreeTrajectoryState& fts,
00032 PropagationDirection dir ) const {
00033
00034 float z0 = fts.position().z();
00035 float zm = fts.momentum().z();
00036
00037 bool inOut = outward(fts);
00038
00039 vector<const DetLayer*> output;
00040
00041
00042
00043 if (inOut) {
00044 if ((zm * z0) >= 0) {
00045 if (theBarrelFlag) inOutBarrel(fts,output);
00046 if (theEndcapFlag) {
00047 if ( z0 >= 0 ) inOutForward(fts,output);
00048 else inOutBackward(fts,output);
00049 }
00050 } else {
00051 if (theEndcapFlag) {
00052 if ( z0 >= 0 ) outInForward(fts,output);
00053 else outInBackward(fts,output);
00054 }
00055 if (theBarrelFlag) inOutBarrel(fts,output);
00056 }
00057 } else {
00058 if ((zm * z0) >= 0) {
00059 if (theBarrelFlag) outInBarrel(fts,output);
00060 if (theEndcapFlag) {
00061 if ( z0 >= 0 ) inOutForward(fts,output);
00062 else inOutBackward(fts,output);
00063 }
00064 } else {
00065 if (theEndcapFlag) {
00066 if ( z0 >= 0 ) outInForward(fts,output);
00067 else outInBackward(fts,output);
00068 }
00069 if (theBarrelFlag) outInBarrel(fts,output);
00070 }
00071 }
00072
00073 if ( dir == oppositeToMomentum ) std::reverse(output.begin(),output.end());
00074
00075 return output;
00076 }
00077
00078
00079
00080
00081
00082 vector<const DetLayer*>
00083 DirectMuonNavigation::compatibleEndcapLayers( const FreeTrajectoryState& fts,
00084 PropagationDirection dir ) const {
00085
00086 float zm = fts.momentum().z();
00087
00088 vector<const DetLayer*> output;
00089
00090
00091 outInBackward(fts,output);
00092 inOutForward(fts,output);
00093
00094
00095 if ( ( zm > 0 && dir == oppositeToMomentum ) ||
00096 ( zm < 0 && dir == alongMomentum ) )
00097 std::reverse(output.begin(),output.end());
00098
00099 return output;
00100 }
00101
00102 void DirectMuonNavigation::inOutBarrel(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
00103
00104 bool cont = false;
00105 const vector<DetLayer*>& barrel = theMuonDetLayerGeometry->allBarrelLayers();
00106
00107 for (vector<DetLayer*>::const_iterator iter_B = barrel.begin(); iter_B != barrel.end(); iter_B++){
00108
00109 if( cont ) output.push_back((*iter_B));
00110 else if ( checkCompatible(fts,dynamic_cast<const BarrelDetLayer*>(*iter_B))) {
00111 output.push_back((*iter_B));
00112 cont = true;
00113 }
00114 }
00115 }
00116
00117
00118 void DirectMuonNavigation::outInBarrel(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
00119
00120
00121 const vector<DetLayer*>& barrel = theMuonDetLayerGeometry->allBarrelLayers();
00122
00123 bool cont = false;
00124 vector<DetLayer*>::const_iterator rbegin = barrel.end();
00125 rbegin--;
00126 vector<DetLayer*>::const_iterator rend = barrel.begin();
00127 rend--;
00128
00129 for (vector<DetLayer*>::const_iterator iter_B = rbegin; iter_B != rend; iter_B--){
00130 if( cont ) output.push_back((*iter_B));
00131 else if ( checkCompatible(fts,dynamic_cast<BarrelDetLayer*>(*iter_B))) {
00132 output.push_back((*iter_B));
00133 cont = true;
00134 }
00135 }
00136 }
00137
00138 void DirectMuonNavigation::inOutForward(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
00139
00140 const vector<DetLayer*>& forward = theMuonDetLayerGeometry->allForwardLayers();
00141 bool cont = false;
00142 for (vector<DetLayer*>::const_iterator iter_E = forward.begin(); iter_E != forward.end();
00143 iter_E++){
00144 if( cont ) output.push_back((*iter_E));
00145 else if ( checkCompatible(fts,dynamic_cast<ForwardDetLayer*>(*iter_E))) {
00146 output.push_back((*iter_E));
00147 cont = true;
00148 }
00149 }
00150 }
00151
00152 void DirectMuonNavigation::outInForward(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
00153
00154
00155 bool cont = false;
00156 const vector<DetLayer*>& forward = theMuonDetLayerGeometry->allForwardLayers();
00157 vector<DetLayer*>::const_iterator rbegin = forward.end();
00158 rbegin--;
00159 vector<DetLayer*>::const_iterator rend = forward.begin();
00160 rend--;
00161 for (vector<DetLayer*>::const_iterator iter_E = rbegin; iter_E != rend;
00162 iter_E--){
00163 if( cont ) output.push_back((*iter_E));
00164 else if ( checkCompatible(fts,dynamic_cast<ForwardDetLayer*>(*iter_E))) {
00165 output.push_back((*iter_E));
00166 cont = true;
00167 }
00168 }
00169 }
00170
00171 void DirectMuonNavigation::inOutBackward(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
00172 bool cont = false;
00173 const vector<DetLayer*>& backward = theMuonDetLayerGeometry->allBackwardLayers();
00174
00175 for (vector<DetLayer*>::const_iterator iter_E = backward.begin(); iter_E != backward.end();
00176 iter_E++){
00177 if( cont ) output.push_back((*iter_E));
00178 else if ( checkCompatible(fts,dynamic_cast<ForwardDetLayer*>(*iter_E))) {
00179 output.push_back((*iter_E));
00180 cont = true;
00181 }
00182 }
00183 }
00184
00185 void DirectMuonNavigation::outInBackward(const FreeTrajectoryState& fts, vector<const DetLayer*>& output) const {
00186
00187 bool cont = false;
00188 const vector<DetLayer*>& backward = theMuonDetLayerGeometry->allBackwardLayers();
00189
00190 vector<DetLayer*>::const_iterator rbegin = backward.end();
00191 rbegin--;
00192 vector<DetLayer*>::const_iterator rend = backward.begin();
00193 rend--;
00194 for (vector<DetLayer*>::const_iterator iter_E = rbegin; iter_E != rend;
00195 iter_E--){
00196 if( cont ) output.push_back((*iter_E));
00197 else if ( checkCompatible(fts,dynamic_cast<ForwardDetLayer*>(*iter_E))) {
00198 output.push_back((*iter_E));
00199 cont = true;
00200 }
00201 }
00202
00203 }
00204
00205
00206 bool DirectMuonNavigation::checkCompatible(const FreeTrajectoryState& fts,const BarrelDetLayer* dl) const {
00207
00208 float z0 = fts.position().z();
00209 float r0 = fts.position().perp();
00210 float zm = fts.momentum().z();
00211 float rm = fts.momentum().perp();
00212 float slope = zm/rm;
00213 if (!outward(fts) ) slope = -slope;
00214 const BoundCylinder bc = dl->specificSurface();
00215 float radius = bc.radius();
00216 float length = bc.bounds().length()/2.;
00217
00218 float z1 = slope*(radius - r0) + z0;
00219 return ( fabs(z1) <= fabs(length)+epsilon_ );
00220
00221 }
00222
00223 bool DirectMuonNavigation::checkCompatible(const FreeTrajectoryState& fts,const ForwardDetLayer* dl) const {
00224
00225 float z0 = fts.position().z();
00226 float r0 = fts.position().perp();
00227 float zm = fts.momentum().z();
00228 float rm = fts.momentum().perp();
00229 float slope = rm/zm;
00230
00231 if (!outward(fts) ) slope = -slope;
00232
00233 const BoundDisk bd = dl->specificSurface();
00234
00235 float outRadius = bd.outerRadius();
00236 float inRadius = bd.innerRadius();
00237 float z = bd.position().z();
00238
00239 float r1 = slope*(z - z0) + r0;
00240 return (r1 >= inRadius-epsilon_ && r1 <= outRadius+epsilon_);
00241
00242 }
00243
00244 bool DirectMuonNavigation::outward(const FreeTrajectoryState& fts) const {
00245
00246
00247
00248 float x0 = fts.position().x();
00249 float y0 = fts.position().y();
00250
00251 float xm = fts.momentum().x();
00252 float ym = fts.momentum().y();
00253
00254 return ((x0 * xm + y0 * ym ) > 0);
00255
00256
00257 }