00001 #include "RecoTracker/TkDetLayers/interface/TIBLayer.h"
00002
00003 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00004
00005 #include "DataFormats/GeometrySurface/interface/SimpleCylinderBounds.h"
00006
00007 #include "RecoTracker/TkDetLayers/interface/LayerCrossingSide.h"
00008 #include "RecoTracker/TkDetLayers/interface/DetGroupMerger.h"
00009 #include "RecoTracker/TkDetLayers/interface/CompatibleDetToGroupAdder.h"
00010
00011 #include "TrackingTools/DetLayers/interface/DetLayerException.h"
00012 #include "TrackingTools/PatternTools/interface/MeasurementEstimator.h"
00013 #include "TrackingTools/GeomPropagators/interface/HelixBarrelCylinderCrossing.h"
00014 #include "TrackingTools/DetLayers/src/DetLessZ.h"
00015
00016 using namespace std;
00017
00018 typedef GeometricSearchDet::DetWithState DetWithState;
00019
00020 TIBLayer::TIBLayer(vector<const TIBRing*>& innerRings,
00021 vector<const TIBRing*>& outerRings) :
00022 theInnerComps(innerRings.begin(),innerRings.end()),
00023 theOuterComps(outerRings.begin(),outerRings.end())
00024 {
00025 theComps.assign(theInnerComps.begin(),theInnerComps.end());
00026 theComps.insert(theComps.end(),theOuterComps.begin(),theOuterComps.end());
00027
00028 sort(theComps.begin(),theComps.end(),DetLessZ());
00029 sort(theInnerComps.begin(),theInnerComps.end(),DetLessZ());
00030 sort(theOuterComps.begin(),theOuterComps.end(),DetLessZ());
00031
00032 for(vector<const GeometricSearchDet*>::const_iterator it=theComps.begin();
00033 it!=theComps.end();it++){
00034 theBasicComps.insert(theBasicComps.end(),
00035 (**it).basicComponents().begin(),
00036 (**it).basicComponents().end());
00037 }
00038
00039
00040 theInnerCylinder = cylinder( theInnerComps);
00041 theOuterCylinder = cylinder( theOuterComps);
00042 initialize();
00043
00044 LogDebug("TkDetLayers") << "==== DEBUG TIBLayer =====" ;
00045 LogDebug("TkDetLayers") << "innerCyl radius, thickness, lenght: "
00046 << theInnerCylinder->radius() << " , "
00047 << theInnerCylinder->bounds().thickness() << " , "
00048 << theInnerCylinder->bounds().length() ;
00049
00050 LogDebug("TkDetLayers") << "outerCyl radius, thickness, lenght: "
00051 << theOuterCylinder->radius() << " , "
00052 << theOuterCylinder->bounds().thickness() << " , "
00053 << theOuterCylinder->bounds().length() ;
00054
00055 LogDebug("TkDetLayers") << "Cyl radius, thickness, lenght: "
00056 << specificSurface().radius() << " , "
00057 << specificSurface().bounds().thickness() << " , "
00058 << specificSurface().bounds().length() ;
00059
00060 for (vector<const GeometricSearchDet*>::const_iterator i=theInnerComps.begin();
00061 i != theInnerComps.end(); i++){
00062 LogDebug("TkDetLayers") << "inner TIBRing pos z,radius,eta,phi: "
00063 << (**i).position().z() << " , "
00064 << (**i).position().perp() << " , "
00065 << (**i).position().eta() << " , "
00066 << (**i).position().phi() ;
00067 }
00068
00069 for (vector<const GeometricSearchDet*>::const_iterator i=theOuterComps.begin();
00070 i != theOuterComps.end(); i++){
00071 LogDebug("TkDetLayers") << "outer TIBRing pos z,radius,eta,phi: "
00072 << (**i).position().z() << " , "
00073 << (**i).position().perp() << " , "
00074 << (**i).position().eta() << " , "
00075 << (**i).position().phi() ;
00076 }
00077
00078
00079
00080
00081
00082
00083
00084 theInnerBinFinder = GeneralBinFinderInZforGeometricSearchDet<float>(theInnerComps.begin(),
00085 theInnerComps.end());
00086
00087 theOuterBinFinder = GeneralBinFinderInZforGeometricSearchDet<float>(theOuterComps.begin(),
00088 theOuterComps.end());
00089 }
00090
00091 TIBLayer::~TIBLayer(){
00092 vector<const GeometricSearchDet*>::const_iterator i;
00093 for (i=theComps.begin(); i!=theComps.end(); i++) {
00094 delete *i;
00095 }
00096 }
00097
00098
00099
00100
00101 BoundCylinder*
00102 TIBLayer::cylinder( const vector<const GeometricSearchDet*>& rings)
00103 {
00104 float leftPos = rings.front()->surface().position().z();
00105 float rightPos = rings.back()->surface().position().z();
00106
00107 const BoundCylinder & frontRing = static_cast<const BoundCylinder &>(rings.front()->surface());
00108 const BoundCylinder & backRing = static_cast<const BoundCylinder &>(rings.back()->surface());
00109 float r = frontRing.radius();
00110 const Bounds& leftBounds = frontRing.bounds();
00111 const Bounds& rightBounds = backRing.bounds();
00112
00113
00114
00115
00116
00117 float thick = leftBounds.thickness() / 2;
00118 float zmin = leftPos - leftBounds.length() / 2;
00119 float zmax = rightPos + rightBounds.length() / 2;
00120 float rmin = r-thick;
00121 float rmax = r+thick;
00122 float zpos = 0.5*(leftPos+rightPos);
00123
00124 return new BoundCylinder( Surface::PositionType( 0, 0, zpos),
00125 rings.front()->surface().rotation(),
00126 SimpleCylinderBounds( rmin, rmax,
00127 zmin-zpos, zmax-zpos));
00128 }
00129
00130
00131
00132 void
00133 TIBLayer::groupedCompatibleDetsV( const TrajectoryStateOnSurface& tsos,
00134 const Propagator& prop,
00135 const MeasurementEstimator& est,
00136 std::vector<DetGroup> & result) const {
00137 SubLayerCrossings crossings;
00138 crossings = computeCrossings( tsos, prop.propagationDirection());
00139 if(! crossings.isValid()) return;
00140
00141 vector<DetGroup> closestResult;
00142 addClosest( tsos, prop, est, crossings.closest(), closestResult);
00143
00144 if (closestResult.empty()) return;
00145
00146
00147 DetGroupElement closestGel( closestResult.front().front());
00148 float window = computeWindowSize( closestGel.det(), closestGel.trajectoryState(), est);
00149
00150 searchNeighbors( tsos, prop, est, crossings.closest(), window,
00151 closestResult, false);
00152
00153 vector<DetGroup> nextResult;
00154 searchNeighbors( tsos, prop, est, crossings.other(), window,
00155 nextResult, true);
00156
00157 int crossingSide = LayerCrossingSide().barrelSide( closestGel.trajectoryState(), prop);
00158 DetGroupMerger::orderAndMergeTwoLevels( closestResult, nextResult, result,
00159 crossings.closestIndex(), crossingSide);
00160 }
00161
00162 SubLayerCrossings TIBLayer::computeCrossings( const TrajectoryStateOnSurface& startingState,
00163 PropagationDirection propDir) const
00164 {
00165 GlobalPoint startPos( startingState.globalPosition());
00166 GlobalVector startDir( startingState.globalMomentum());
00167 double rho( startingState.transverseCurvature());
00168
00169 HelixBarrelCylinderCrossing innerCrossing( startPos, startDir, rho,
00170 propDir,*theInnerCylinder);
00171 if (!innerCrossing.hasSolution()) return SubLayerCrossings();
00172
00173 GlobalPoint gInnerPoint( innerCrossing.position());
00174 int innerIndex = theInnerBinFinder.binIndex(gInnerPoint.z());
00175 const GeometricSearchDet* innerRing( theInnerComps[innerIndex]);
00176 float innerDist = fabs( innerRing->surface().position().z() - gInnerPoint.z());
00177 SubLayerCrossing innerSLC( 0, innerIndex, gInnerPoint);
00178
00179 HelixBarrelCylinderCrossing outerCrossing( startPos, startDir, rho,
00180 propDir,*theOuterCylinder);
00181 if (!outerCrossing.hasSolution()) return SubLayerCrossings();
00182
00183 GlobalPoint gOuterPoint( outerCrossing.position());
00184 int outerIndex = theOuterBinFinder.binIndex(gOuterPoint.z());
00185 const GeometricSearchDet* outerRing( theOuterComps[outerIndex]);
00186 float outerDist = fabs( outerRing->surface().position().z() - gOuterPoint.z());
00187 SubLayerCrossing outerSLC( 1, outerIndex, gOuterPoint);
00188
00189 if (innerDist < outerDist) {
00190 return SubLayerCrossings( innerSLC, outerSLC, 0);
00191 }
00192 else {
00193 return SubLayerCrossings( outerSLC, innerSLC, 1);
00194 }
00195 }
00196
00197 bool TIBLayer::addClosest( const TrajectoryStateOnSurface& tsos,
00198 const Propagator& prop,
00199 const MeasurementEstimator& est,
00200 const SubLayerCrossing& crossing,
00201 vector<DetGroup>& result) const
00202 {
00203
00204
00205 const vector<const GeometricSearchDet*>& sub( subLayer( crossing.subLayerIndex()));
00206 const Det* det(sub[crossing.closestDetIndex()]);
00207 return CompatibleDetToGroupAdder().add( *det, tsos, prop, est, result);
00208 }
00209
00210 void TIBLayer::searchNeighbors( const TrajectoryStateOnSurface& tsos,
00211 const Propagator& prop,
00212 const MeasurementEstimator& est,
00213 const SubLayerCrossing& crossing,
00214 float window,
00215 vector<DetGroup>& result,
00216 bool checkClosest) const
00217 {
00218 GlobalPoint gCrossingPos = crossing.position();
00219
00220 const vector<const GeometricSearchDet*>& sLayer( subLayer( crossing.subLayerIndex()));
00221
00222 int closestIndex = crossing.closestDetIndex();
00223 int negStartIndex = closestIndex-1;
00224 int posStartIndex = closestIndex+1;
00225
00226 if (checkClosest) {
00227 if (gCrossingPos.z() < sLayer[closestIndex]->surface().position().z()) {
00228 posStartIndex = closestIndex;
00229 }
00230 else {
00231 negStartIndex = closestIndex;
00232 }
00233 }
00234
00235 typedef CompatibleDetToGroupAdder Adder;
00236 for (int idet=negStartIndex; idet >= 0; idet--) {
00237 const GeometricSearchDet* neighborRing = sLayer[idet];
00238 if (!overlap( gCrossingPos, *neighborRing, window)) break;
00239 if (!Adder::add( *neighborRing, tsos, prop, est, result)) break;
00240 }
00241 for (int idet=posStartIndex; idet < static_cast<int>(sLayer.size()); idet++) {
00242 const GeometricSearchDet* neighborRing = sLayer[idet];
00243 if (!overlap( gCrossingPos, *neighborRing, window)) break;
00244 if (!Adder::add( *neighborRing, tsos, prop, est, result)) break;
00245 }
00246 }
00247
00248 bool TIBLayer::overlap( const GlobalPoint& crossPoint,
00249 const GeometricSearchDet& det,
00250 float window) const
00251 {
00252 float halfLength = det.surface().bounds().length()/2.;
00253
00254
00255
00256 if ( fabs( crossPoint.z()-det.position().z()) < (halfLength + window)) {
00257
00258 return true;
00259 } else {
00260
00261 return false;
00262 }
00263 }
00264
00265 float TIBLayer::computeWindowSize( const GeomDet* det,
00266 const TrajectoryStateOnSurface& tsos,
00267 const MeasurementEstimator& est) const
00268 {
00269
00270
00271
00272
00273
00274 MeasurementEstimator::Local2DVector localError( est.maximalLocalDisplacement(tsos, det->surface()));
00275 float yError = localError.y();
00276
00277 float tanTheta = tan( tsos.globalMomentum().theta());
00278 float thickCorrection = det->surface().bounds().thickness() / (2.*fabs( tanTheta));
00279
00280
00281
00282 return yError + thickCorrection;
00283 }
00284
00285
00286