Producers/src/HitDropper.cc

// $Id: HitDropper.cc,v 1.9 2009/10/12 21:41:09 loizides Exp $

#include "MitEdm/Producers/interface/HitDropper.h"
#include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
#include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
#include "DataFormats/SiStripDetId/interface/TIBDetId.h"
#include "DataFormats/SiStripDetId/interface/TIDDetId.h"
#include "DataFormats/SiStripDetId/interface/TOBDetId.h"
#include "DataFormats/SiStripDetId/interface/TECDetId.h"
#include "DataFormats/MuonDetId/interface/DTLayerId.h"
#include "DataFormats/MuonDetId/interface/CSCDetId.h"
#include "DataFormats/MuonDetId/interface/RPCDetId.h"
#include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"

using namespace edm;
using namespace mitedm;

//--------------------------------------------------------------------------------------------------
reco::HitPattern HitDropper::CorrectedHits(const reco::TransientTrack *tTrack,
                                           const ThreeVector &vtxPos) const
{
  // Return reco::HitPattern structure for the given track with all hits occuring before vtxPos
  // on the track (relative to the primary vertex if given) removed.
  // This version of the function uses an iterative helix-plane intersector and does not (yet)
  // take into account the uncertainties in vertex position.

  const GlobalPoint vtxPoint(vtxPos.x(),vtxPos.y(),vtxPos.z());
  const TrajectoryStateClosestToPoint vtxTSCP = tTrack->trajectoryStateClosestToPoint(vtxPoint);

  reco::HitPattern hitPattern;
  int nHits = 0;
  for (uint hi=0; hi<tTrack->recHitsSize(); ++hi) {
    const TrackingRecHit *hit = tTrack->recHit(hi).get();
    DetId geoId = hit->geographicalId();
    if(geoId == uint32_t(0)) continue;
    const GeomDet *det = trackerGeo_->idToDet(geoId);
  
    HelixArbitraryPlaneCrossing c(HelixPlaneCrossing::PositionType(vtxTSCP.theState().position()),
                                  HelixPlaneCrossing::DirectionType(vtxTSCP.theState().momentum()),
                                  vtxTSCP.theState().transverseCurvature(),
                                  anyDirection);                                        
                                            
    std::pair<bool,double> crossResult = c.pathLength(det->surface());
    if ( crossResult.first && crossResult.second >= 0 ) {
      hitPattern.set(*hit,nHits);
      nHits++;
    }
  }

  return hitPattern;
}

//--------------------------------------------------------------------------------------------------
reco::HitPattern HitDropper::CorrectedHits(const reco::Track *track,
                                           const ThreeVector &vtxPos) const
{
  // Build the transient track and then return the corrected HitPattern.

  reco::TransientTrack tTrack = builder_->build(track);
  return CorrectedHits(&tTrack, vtxPos);
}

//--------------------------------------------------------------------------------------------------
reco::HitPattern HitDropper::CorrectedHits(const reco::Track *track,
                                           const ThreeVector &vtxPos,
                                           const ThreeVector &trkMom,
                                           Double_t lxyError,
                                           Double_t lzError,
                                           Double_t sigmaTolerance) const
{
  // Return reco::HitPattern structure for the given track with all hits occuring before vtxPos
  // on the track (relative to the primary vertex if given) removed.
  // This version of the function determines this completely analytically, and taking the
  // vertex position uncertainty into account, which might be important for particles which decay
  // within a tracker layer.

  const LocalPoint nullLocal(0,0,0);
  
  reco::HitPattern hitPattern;
  int nHits = 0;
  for (uint hi=0; hi<track->recHitsSize(); ++hi) {
    const TrackingRecHit *hit = track->recHit(hi).get();
    DetId geoId = hit->geographicalId();
    if(geoId == uint32_t(0)) continue;
    const GeomDet *det = trackerGeo_->idToDet(geoId);
  
    const GlobalPoint cPos = det->toGlobal(nullLocal);
    
    const ThreeVector crossPos(cPos.x(), cPos.y(), cPos.z());
    const ThreeVector delta = crossPos - vtxPos;
       
    Double_t lengthOverSigma = 0;
    
    //calculate distance between vertex and approximate hit position, projected into
    //the appropriate plane/axis and compared to the uncertainty in vertex position
    //in that plane/axis
    if (IsBarrel(det)) {
      const ThreeVector trkMomXY(trkMom.x(), trkMom.y(), 0.0);
      Double_t deltaXY = delta.Dot(trkMomXY)/trkMomXY.R();
      lengthOverSigma = deltaXY/lxyError;
    }
    else if (IsDisk(det)) {
      Double_t deltaZ = delta.z()*trkMom.z()/fabs(trkMom.z());
      lengthOverSigma = deltaZ/lzError;
    }
    else
      throw edm::Exception(edm::errors::Configuration, "HitDropper::CorrectedHits\n")
         << "Error! Detector element not in a valid barrel or disk layer." << std::endl; 
    
    //add the hit only if it is after the vertex, 
    //allowing for some uncertainty in the vertex position
    if ( lengthOverSigma>(-sigmaTolerance) ) {
      hitPattern.set(*hit,nHits);
      nHits++;
    } 
  }
  return hitPattern;
}


//--------------------------------------------------------------------------------------------------
reco::HitPattern HitDropper::CorrectedHitsAOD(const reco::Track *track,
                                           const ThreeVector &vtxPos,
                                           const ThreeVector &trkMom,
                                           Double_t lxyError,
                                           Double_t lzError,
                                           Double_t sigmaTolerance) const
{
  // Return reco::HitPattern structure for the given track with all hits occuring before vtxPos
  // on the track (relative to the primary vertex if given) removed.
  // This version of the function determines this completely analytically, and taking the
  // vertex position uncertainty into account, which might be important for particles which decay
  // within a tracker layer.  In order to function in AOD, the tracker geometry is used in an
  // approximate way, at the level of barrel and disk layers.  The layer thickness is taken
  // into account in the uncertainty tolerance.
  

  //Figure out whether we should use the +z or -z side in the case of a disk layer
  int side = 0;
  if (track->pz() < 0)
    side = -1;
  else
    side = 1;
  
  int nHits = 0;
  reco::HitPattern hitPattern;
  const reco::HitPattern* inHitPatterns[3];
  inHitPatterns[0] = &track->hitPattern();
  inHitPatterns[1] = &track->trackerExpectedHitsInner();
  inHitPatterns[2] = &track->trackerExpectedHitsOuter();
  
  for (Int_t hp=0; hp<3; ++hp) {
    const reco::HitPattern &inhits = *inHitPatterns[hp];
    for (Int_t hi=0; hi<inhits.numberOfHits(); hi++) {
      uint32_t hit = inhits.getHitPattern(hi);
      uint32_t subdet = inhits.getSubStructure(hit);
      uint32_t layerid = inhits.getLayer(hit);
      const DetLayer *det = FindLayer(subdet,layerid,side);
      if (!det) continue;
          
      bool goodHit = false;
      
      //calculate distance between vertex and approximate hit position, projected into
      //the appropriate plane/axis and compared to the uncertainty in vertex position
      //in that plane/axis
      if (det->location()==GeomDetEnumerators::barrel) {
        const BarrelDetLayer *barrelDet = static_cast<const BarrelDetLayer*>(det);
        
        const double barrelRho = barrelDet->specificSurface().radius();
        const double barrelHalfThickness = barrelDet->specificSurface().bounds().thickness()/2.0;
        const ThreeVector crossPos(barrelRho*cos(track->phi()),barrelRho*sin(track->phi()),0.0);
        
        const ThreeVector delta = crossPos - vtxPos;
        const ThreeVector trkMomXY(trkMom.x(), trkMom.y(), 0.0);
        const double deltaXY = delta.Dot(trkMomXY)/trkMomXY.R();
        goodHit = deltaXY>(-sigmaTolerance*lxyError - barrelHalfThickness);
      }
      else if (det->location()==GeomDetEnumerators::endcap) {
        const ForwardDetLayer *forwardDet = static_cast<const ForwardDetLayer*>(det);
        const double diskZ = forwardDet->specificSurface().position().z();
        const double diskHalfThickness = forwardDet->specificSurface().bounds().thickness()/2.0;
        const double deltaZ = (diskZ - vtxPos.z())*trkMom.z()/fabs(trkMom.z());;
        goodHit = deltaZ>(-sigmaTolerance*lzError - diskHalfThickness);
      }
      else
        throw edm::Exception(edm::errors::Configuration, "HitDropper::CorrectedHitsAOD\n")
          << "Error! Detector element not in a valid barrel or disk layer." << std::endl; 
      
      //add the hit only if it is after the vertex, 
      //allowing for some uncertainty in the vertex position
      if ( goodHit ) {
        bool isStereo = inhits.getSide(hit);
        DetId dummyid = det->basicComponents().front()->geographicalId();
        if (isStereo) {
          dummyid = StereoDetId(dummyid);
        }
          
        const TrackingRecHit::Type hitType = static_cast<TrackingRecHit::Type>(inhits.getHitType(hit));
        InvalidTrackingRecHit dummyhit(dummyid, hitType);
        hitPattern.set(dummyhit,nHits);
        if ( hit != hitPattern.getHitPattern(nHits) )
          throw edm::Exception(edm::errors::Configuration, "HitDropper::CorrectedHitsAOD\n")
          << "Error! Mismatch in copying hit pattern." << std::endl; 
          
        nHits++;
      } 
    }
  }

  return hitPattern;
}

//--------------------------------------------------------------------------------------------------
const DetLayer *HitDropper::FindLayer(int subdet, int layer, int side) const
{
  // Find a given layer.

  switch(subdet) {
    case StripSubdetector::TIB:
      //edm::LogInfo(TkDetLayers) << "TIB layer n: " << TIBDetId(id).layer() ;
      return trackerGeoSearch_->tibLayers().at(layer-1);
      break;

    case StripSubdetector::TOB:
      //edm::LogInfo(TkDetLayers) << "TOB layer n: " << TOBDetId(id).layer() ;
      return trackerGeoSearch_->tobLayers().at(layer-1);
      break;

    case StripSubdetector::TID:
      //edm::LogInfo(TkDetLayers) << "TID wheel n: " << TIDDetId(id).wheel() ;
      if( side == -1 ) {
        return trackerGeoSearch_->negTidLayers().at(layer-1);
      }else if( side == 1 ) {
        return trackerGeoSearch_->posTidLayers().at(layer-1);
      }
      break;

    case StripSubdetector::TEC:
      //edm::LogInfo(TkDetLayers) << "TEC wheel n: " << TECDetId(id).wheel() ;
      if( side == -1 ) {
        return trackerGeoSearch_->negTecLayers().at(layer-1);
      } else if( side == 1 ) {
        return trackerGeoSearch_->posTecLayers().at(layer-1);
      }
      break;

    case PixelSubdetector::PixelBarrel:
      //edm::LogInfo(TkDetLayers) << "PixelBarrel layer n: " << PXBDetId(id).layer() ;
      return trackerGeoSearch_->pixelBarrelLayers().at(layer-1);
      break;

    case PixelSubdetector::PixelEndcap:
      //edm::LogInfo(TkDetLayers) << "PixelEndcap disk n: " << PXFDetId(id).disk() ;
      if(side == -1 ) {
        return trackerGeoSearch_->negPixelForwardLayers().at(layer-1);
      }else if( side == 1 ) {
        return trackerGeoSearch_->posPixelForwardLayers().at(layer-1);
      }
      break;

    default:    
      return 0;
      // throw(something);
  }
  return 0; //just to avoid compile warnings
}

//--------------------------------------------------------------------------------------------------
DetId HitDropper::StereoDetId(const DetId &i) const
{
  // Deal with stereo det id.

  switch (i.det()) {
    case DetId::Tracker:
      switch (i.subdetId()) {
        case PixelSubdetector::PixelBarrel:
        case PixelSubdetector::PixelEndcap:
          return 0;
        case StripSubdetector::TIB:
        {
          TIBDetId id = i;
          if (id.isStereo())
            return id;
          else {
            const std::vector<unsigned int> idString = id.string();
            return TIBDetId(id.layer(),idString[0],idString[1],idString[2],id.module(),1);
          }
        }
        case StripSubdetector::TID:
        {
          TIDDetId id = i;
          if (id.isStereo())
            return id;
          else {
            const std::vector<unsigned int> idModule = id.module();
            return TIDDetId(id.side(),id.wheel(),id.ring(),idModule[0],idModule[1],1);
          }
        }
      case StripSubdetector::TOB:
      {
        TOBDetId id = i;
        if (id.isStereo())
          return id;
        else {
          const std::vector<unsigned int> idRod = id.rod();
          return TOBDetId(id.layer(),idRod[0],idRod[1],id.module(),1);
        }
      }
        case StripSubdetector::TEC:
        {
            TECDetId id = i;
            if (id.isStereo())
              return id;
            else {
              const std::vector<unsigned int> idPetal = id.petal();
              return TECDetId(id.side(),id.wheel(),idPetal[0],idPetal[1],id.ring(),id.module(),1);
            }
        }
        default:
          return 0;
      }
      break;
    default:
      return 0;
  }
}
Revision:	1.10
Committed:	Mon Nov 2 22:55:58 2009 UTC (15 years, 6 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_012f, Mit_012e, Mit_012d, Mit_012c, Mit_012b, Mit_012a, Mit_012
Changes since 1.9:	+59 -52 lines
Log Message:	Use secondary hit patterns
#	User	Rev	Content
1	bendavid	1.10	// $Id: HitDropper.cc,v 1.9 2009/10/12 21:41:09 loizides Exp $
2	bendavid	1.1
3			#include "MitEdm/Producers/interface/HitDropper.h"
4	bendavid	1.2	#include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
5			#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
6	bendavid	1.8	#include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
7			#include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
8			#include "DataFormats/SiStripDetId/interface/TIBDetId.h"
9			#include "DataFormats/SiStripDetId/interface/TIDDetId.h"
10			#include "DataFormats/SiStripDetId/interface/TOBDetId.h"
11			#include "DataFormats/SiStripDetId/interface/TECDetId.h"
12			#include "DataFormats/MuonDetId/interface/DTLayerId.h"
13			#include "DataFormats/MuonDetId/interface/CSCDetId.h"
14			#include "DataFormats/MuonDetId/interface/RPCDetId.h"
15	bendavid	1.2	#include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
16	bendavid	1.1
17			using namespace edm;
18			using namespace mitedm;
19
20			//--------------------------------------------------------------------------------------------------
21			reco::HitPattern HitDropper::CorrectedHits(const reco::TransientTrack *tTrack,
22	loizides	1.4	const ThreeVector &vtxPos) const
23	bendavid	1.1	{
24	loizides	1.4	// Return reco::HitPattern structure for the given track with all hits occuring before vtxPos
25			// on the track (relative to the primary vertex if given) removed.
26			// This version of the function uses an iterative helix-plane intersector and does not (yet)
27			// take into account the uncertainties in vertex position.
28	bendavid	1.2
29			const GlobalPoint vtxPoint(vtxPos.x(),vtxPos.y(),vtxPos.z());
30			const TrajectoryStateClosestToPoint vtxTSCP = tTrack->trajectoryStateClosestToPoint(vtxPoint);
31
32	bendavid	1.1	reco::HitPattern hitPattern;
33			int nHits = 0;
34			for (uint hi=0; hi<tTrack->recHitsSize(); ++hi) {
35			const TrackingRecHit *hit = tTrack->recHit(hi).get();
36	mrudolph	1.3	DetId geoId = hit->geographicalId();
37			if(geoId == uint32_t(0)) continue;
38			const GeomDet *det = trackerGeo_->idToDet(geoId);
39	bendavid	1.2
40	loizides	1.7	HelixArbitraryPlaneCrossing c(HelixPlaneCrossing::PositionType(vtxTSCP.theState().position()),
41			HelixPlaneCrossing::DirectionType(vtxTSCP.theState().momentum()),
42			vtxTSCP.theState().transverseCurvature(),
43			anyDirection);
44	bendavid	1.2
45	loizides	1.7	std::pair<bool,double> crossResult = c.pathLength(det->surface());
46	bendavid	1.2	if ( crossResult.first && crossResult.second >= 0 ) {
47	bendavid	1.1	hitPattern.set(*hit,nHits);
48			nHits++;
49			}
50			}
51	bendavid	1.2
52	bendavid	1.1	return hitPattern;
53			}
54
55			//--------------------------------------------------------------------------------------------------
56			reco::HitPattern HitDropper::CorrectedHits(const reco::Track *track,
57	loizides	1.4	const ThreeVector &vtxPos) const
58	bendavid	1.1	{
59	loizides	1.4	// Build the transient track and then return the corrected HitPattern.
60	bendavid	1.1
61			reco::TransientTrack tTrack = builder_->build(track);
62	bendavid	1.2	return CorrectedHits(&tTrack, vtxPos);
63			}
64
65			//--------------------------------------------------------------------------------------------------
66			reco::HitPattern HitDropper::CorrectedHits(const reco::Track *track,
67	loizides	1.4	const ThreeVector &vtxPos,
68			const ThreeVector &trkMom,
69			Double_t lxyError,
70			Double_t lzError,
71			Double_t sigmaTolerance) const
72	bendavid	1.2	{
73	loizides	1.4	// Return reco::HitPattern structure for the given track with all hits occuring before vtxPos
74			// on the track (relative to the primary vertex if given) removed.
75			// This version of the function determines this completely analytically, and taking the
76			// vertex position uncertainty into account, which might be important for particles which decay
77			// within a tracker layer.
78	bendavid	1.2
79	bendavid	1.8	const LocalPoint nullLocal(0,0,0);
80	bendavid	1.2
81			reco::HitPattern hitPattern;
82			int nHits = 0;
83			for (uint hi=0; hi<track->recHitsSize(); ++hi) {
84			const TrackingRecHit *hit = track->recHit(hi).get();
85	mrudolph	1.3	DetId geoId = hit->geographicalId();
86			if(geoId == uint32_t(0)) continue;
87			const GeomDet *det = trackerGeo_->idToDet(geoId);
88	bendavid	1.8
89			const GlobalPoint cPos = det->toGlobal(nullLocal);
90	bendavid	1.2
91	loizides	1.7	const ThreeVector crossPos(cPos.x(), cPos.y(), cPos.z());
92	bendavid	1.2	const ThreeVector delta = crossPos - vtxPos;
93	bendavid	1.8
94	bendavid	1.2	Double_t lengthOverSigma = 0;
95
96			//calculate distance between vertex and approximate hit position, projected into
97			//the appropriate plane/axis and compared to the uncertainty in vertex position
98			//in that plane/axis
99			if (IsBarrel(det)) {
100			const ThreeVector trkMomXY(trkMom.x(), trkMom.y(), 0.0);
101			Double_t deltaXY = delta.Dot(trkMomXY)/trkMomXY.R();
102			lengthOverSigma = deltaXY/lxyError;
103			}
104			else if (IsDisk(det)) {
105			Double_t deltaZ = delta.z()*trkMom.z()/fabs(trkMom.z());
106			lengthOverSigma = deltaZ/lzError;
107			}
108			else
109			throw edm::Exception(edm::errors::Configuration, "HitDropper::CorrectedHits\n")
110			<< "Error! Detector element not in a valid barrel or disk layer." << std::endl;
111
112	loizides	1.7	//add the hit only if it is after the vertex,
113			//allowing for some uncertainty in the vertex position
114	bendavid	1.2	if ( lengthOverSigma>(-sigmaTolerance) ) {
115			hitPattern.set(*hit,nHits);
116			nHits++;
117			}
118			}
119			return hitPattern;
120	bendavid	1.1	}
121	bendavid	1.5
122
123			//--------------------------------------------------------------------------------------------------
124			reco::HitPattern HitDropper::CorrectedHitsAOD(const reco::Track *track,
125			const ThreeVector &vtxPos,
126			const ThreeVector &trkMom,
127			Double_t lxyError,
128			Double_t lzError,
129			Double_t sigmaTolerance) const
130			{
131			// Return reco::HitPattern structure for the given track with all hits occuring before vtxPos
132			// on the track (relative to the primary vertex if given) removed.
133			// This version of the function determines this completely analytically, and taking the
134			// vertex position uncertainty into account, which might be important for particles which decay
135	bendavid	1.8	// within a tracker layer. In order to function in AOD, the tracker geometry is used in an
136			// approximate way, at the level of barrel and disk layers. The layer thickness is taken
137			// into account in the uncertainty tolerance.
138
139	bendavid	1.5
140	bendavid	1.8	//Figure out whether we should use the +z or -z side in the case of a disk layer
141			int side = 0;
142			if (track->pz() < 0)
143			side = -1;
144			else
145			side = 1;
146	bendavid	1.5
147	bendavid	1.8	int nHits = 0;
148			reco::HitPattern hitPattern;
149	bendavid	1.10	const reco::HitPattern* inHitPatterns[3];
150			inHitPatterns[0] = &track->hitPattern();
151			inHitPatterns[1] = &track->trackerExpectedHitsInner();
152			inHitPatterns[2] = &track->trackerExpectedHitsOuter();
153
154			for (Int_t hp=0; hp<3; ++hp) {
155			const reco::HitPattern &inhits = *inHitPatterns[hp];
156			for (Int_t hi=0; hi<inhits.numberOfHits(); hi++) {
157			uint32_t hit = inhits.getHitPattern(hi);
158			uint32_t subdet = inhits.getSubStructure(hit);
159			uint32_t layerid = inhits.getLayer(hit);
160			const DetLayer *det = FindLayer(subdet,layerid,side);
161			if (!det) continue;
162
163			bool goodHit = false;
164
165			//calculate distance between vertex and approximate hit position, projected into
166			//the appropriate plane/axis and compared to the uncertainty in vertex position
167			//in that plane/axis
168			if (det->location()==GeomDetEnumerators::barrel) {
169			const BarrelDetLayer barrelDet = static_cast<const BarrelDetLayer>(det);
170
171			const double barrelRho = barrelDet->specificSurface().radius();
172			const double barrelHalfThickness = barrelDet->specificSurface().bounds().thickness()/2.0;
173			const ThreeVector crossPos(barrelRhocos(track->phi()),barrelRhosin(track->phi()),0.0);
174	bendavid	1.8
175	bendavid	1.10	const ThreeVector delta = crossPos - vtxPos;
176			const ThreeVector trkMomXY(trkMom.x(), trkMom.y(), 0.0);
177			const double deltaXY = delta.Dot(trkMomXY)/trkMomXY.R();
178			goodHit = deltaXY>(-sigmaTolerance*lxyError - barrelHalfThickness);
179			}
180			else if (det->location()==GeomDetEnumerators::endcap) {
181			const ForwardDetLayer forwardDet = static_cast<const ForwardDetLayer>(det);
182			const double diskZ = forwardDet->specificSurface().position().z();
183			const double diskHalfThickness = forwardDet->specificSurface().bounds().thickness()/2.0;
184			const double deltaZ = (diskZ - vtxPos.z())*trkMom.z()/fabs(trkMom.z());;
185			goodHit = deltaZ>(-sigmaTolerance*lzError - diskHalfThickness);
186			}
187			else
188			throw edm::Exception(edm::errors::Configuration, "HitDropper::CorrectedHitsAOD\n")
189			<< "Error! Detector element not in a valid barrel or disk layer." << std::endl;
190	bendavid	1.5
191	bendavid	1.10	//add the hit only if it is after the vertex,
192			//allowing for some uncertainty in the vertex position
193			if ( goodHit ) {
194			bool isStereo = inhits.getSide(hit);
195			DetId dummyid = det->basicComponents().front()->geographicalId();
196			if (isStereo) {
197			dummyid = StereoDetId(dummyid);
198			}
199
200			const TrackingRecHit::Type hitType = static_cast<TrackingRecHit::Type>(inhits.getHitType(hit));
201			InvalidTrackingRecHit dummyhit(dummyid, hitType);
202			hitPattern.set(dummyhit,nHits);
203			if ( hit != hitPattern.getHitPattern(nHits) )
204			throw edm::Exception(edm::errors::Configuration, "HitDropper::CorrectedHitsAOD\n")
205			<< "Error! Mismatch in copying hit pattern." << std::endl;
206
207			nHits++;
208			}
209	bendavid	1.8	}
210			}
211
212			return hitPattern;
213			}
214
215			//--------------------------------------------------------------------------------------------------
216			const DetLayer *HitDropper::FindLayer(int subdet, int layer, int side) const
217			{
218	loizides	1.9	// Find a given layer.
219
220	bendavid	1.8	switch(subdet) {
221	loizides	1.9	case StripSubdetector::TIB:
222			//edm::LogInfo(TkDetLayers) << "TIB layer n: " << TIBDetId(id).layer() ;
223			return trackerGeoSearch_->tibLayers().at(layer-1);
224			break;
225
226			case StripSubdetector::TOB:
227			//edm::LogInfo(TkDetLayers) << "TOB layer n: " << TOBDetId(id).layer() ;
228			return trackerGeoSearch_->tobLayers().at(layer-1);
229			break;
230
231			case StripSubdetector::TID:
232			//edm::LogInfo(TkDetLayers) << "TID wheel n: " << TIDDetId(id).wheel() ;
233			if( side == -1 ) {
234			return trackerGeoSearch_->negTidLayers().at(layer-1);
235			}else if( side == 1 ) {
236			return trackerGeoSearch_->posTidLayers().at(layer-1);
237			}
238			break;
239
240			case StripSubdetector::TEC:
241			//edm::LogInfo(TkDetLayers) << "TEC wheel n: " << TECDetId(id).wheel() ;
242			if( side == -1 ) {
243			return trackerGeoSearch_->negTecLayers().at(layer-1);
244			} else if( side == 1 ) {
245			return trackerGeoSearch_->posTecLayers().at(layer-1);
246			}
247			break;
248	bendavid	1.8
249	loizides	1.9	case PixelSubdetector::PixelBarrel:
250			//edm::LogInfo(TkDetLayers) << "PixelBarrel layer n: " << PXBDetId(id).layer() ;
251			return trackerGeoSearch_->pixelBarrelLayers().at(layer-1);
252			break;
253	bendavid	1.8
254	loizides	1.9	case PixelSubdetector::PixelEndcap:
255			//edm::LogInfo(TkDetLayers) << "PixelEndcap disk n: " << PXFDetId(id).disk() ;
256			if(side == -1 ) {
257			return trackerGeoSearch_->negPixelForwardLayers().at(layer-1);
258			}else if( side == 1 ) {
259			return trackerGeoSearch_->posPixelForwardLayers().at(layer-1);
260			}
261			break;
262	bendavid	1.8
263	loizides	1.9	default:
264			return 0;
265			// throw(something);
266	bendavid	1.8	}
267			return 0; //just to avoid compile warnings
268			}
269
270	loizides	1.9	//--------------------------------------------------------------------------------------------------
271	bendavid	1.8	DetId HitDropper::StereoDetId(const DetId &i) const
272			{
273	loizides	1.9	// Deal with stereo det id.
274
275	bendavid	1.8	switch (i.det()) {
276	loizides	1.9	case DetId::Tracker:
277	bendavid	1.8	switch (i.subdetId()) {
278	loizides	1.9	case PixelSubdetector::PixelBarrel:
279			case PixelSubdetector::PixelEndcap:
280			return 0;
281			case StripSubdetector::TIB:
282			{
283			TIBDetId id = i;
284			if (id.isStereo())
285			return id;
286			else {
287			const std::vector<unsigned int> idString = id.string();
288			return TIBDetId(id.layer(),idString[0],idString[1],idString[2],id.module(),1);
289			}
290			}
291			case StripSubdetector::TID:
292			{
293			TIDDetId id = i;
294			if (id.isStereo())
295			return id;
296			else {
297			const std::vector<unsigned int> idModule = id.module();
298			return TIDDetId(id.side(),id.wheel(),id.ring(),idModule[0],idModule[1],1);
299			}
300			}
301	bendavid	1.8	case StripSubdetector::TOB:
302			{
303	loizides	1.9	TOBDetId id = i;
304			if (id.isStereo())
305			return id;
306			else {
307			const std::vector<unsigned int> idRod = id.rod();
308			return TOBDetId(id.layer(),idRod[0],idRod[1],id.module(),1);
309			}
310	bendavid	1.8	}
311	loizides	1.9	case StripSubdetector::TEC:
312			{
313	bendavid	1.8	TECDetId id = i;
314			if (id.isStereo())
315			return id;
316			else {
317			const std::vector<unsigned int> idPetal = id.petal();
318			return TECDetId(id.side(),id.wheel(),idPetal[0],idPetal[1],id.ring(),id.module(),1);
319			}
320	loizides	1.9	}
321			default:
322			return 0;
323	bendavid	1.8	}
324			break;
325	loizides	1.9	default:
326	bendavid	1.8	return 0;
327	bendavid	1.5	}
328			}