Producers/src/HitDropper.cc

// $Id: HitDropper.cc,v 1.10 2009/11/02 22:55:58 bendavid 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;
}

reco::HitPattern HitDropper::SharedHits(const reco::Track *t1, const reco::Track *t2) const
{
  //Return a hit pattern corresponding to the hits on the two tracks which share clusters
  
  int nHits = 0;
  reco::HitPattern sharedHits;
  
  if (!t1->extra().isAvailable() || !t2->extra().isAvailable())
    return sharedHits;
  
  for (trackingRecHit_iterator iHit1 = t1->recHitsBegin();  iHit1 != t1->recHitsEnd(); ++iHit1) { 
    const TrackingRecHit *hit1 = iHit1->get();
    if (hit1->isValid()) {
      for (trackingRecHit_iterator iHit2 = t2->recHitsBegin();  iHit2 != t2->recHitsEnd(); ++iHit2) { 
        const TrackingRecHit *hit2 = iHit2->get();
        if (hit2->isValid() && hit1->sharesInput(hit2,TrackingRecHit::some)) {
          sharedHits.set(*hit1,nHits);
          nHits++;
        }
      }    
    }
  }
  
  return sharedHits;
  
}

//--------------------------------------------------------------------------------------------------
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.11
Committed:	Tue Dec 15 23:27:34 2009 UTC (15 years, 4 months ago) by bendavid
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_014a, Mit_014, Mit_014pre3, Mit_014pre2, Mit_014pre1, Mit_013d, Mit_013c, Mit_013b, Mit_013a, Mit_013, Mit_013pre1, Mit_012i, Mit_012h, Mit_012g
Changes since 1.10:	+28 -1 lines
Log Message:	Add shared layer mask
#	User	Rev	Content
1	bendavid	1.11	// $Id: HitDropper.cc,v 1.10 2009/11/02 22:55:58 bendavid 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	bendavid	1.11	reco::HitPattern HitDropper::SharedHits(const reco::Track t1, const reco::Track t2) const
216			{
217			//Return a hit pattern corresponding to the hits on the two tracks which share clusters
218
219			int nHits = 0;
220			reco::HitPattern sharedHits;
221
222			if (!t1->extra().isAvailable() \|\| !t2->extra().isAvailable())
223			return sharedHits;
224
225			for (trackingRecHit_iterator iHit1 = t1->recHitsBegin(); iHit1 != t1->recHitsEnd(); ++iHit1) {
226			const TrackingRecHit *hit1 = iHit1->get();
227			if (hit1->isValid()) {
228			for (trackingRecHit_iterator iHit2 = t2->recHitsBegin(); iHit2 != t2->recHitsEnd(); ++iHit2) {
229			const TrackingRecHit *hit2 = iHit2->get();
230			if (hit2->isValid() && hit1->sharesInput(hit2,TrackingRecHit::some)) {
231			sharedHits.set(*hit1,nHits);
232			nHits++;
233			}
234			}
235			}
236			}
237
238			return sharedHits;
239
240			}
241
242	bendavid	1.8	//--------------------------------------------------------------------------------------------------
243			const DetLayer *HitDropper::FindLayer(int subdet, int layer, int side) const
244			{
245	loizides	1.9	// Find a given layer.
246
247	bendavid	1.8	switch(subdet) {
248	loizides	1.9	case StripSubdetector::TIB:
249			//edm::LogInfo(TkDetLayers) << "TIB layer n: " << TIBDetId(id).layer() ;
250			return trackerGeoSearch_->tibLayers().at(layer-1);
251			break;
252
253			case StripSubdetector::TOB:
254			//edm::LogInfo(TkDetLayers) << "TOB layer n: " << TOBDetId(id).layer() ;
255			return trackerGeoSearch_->tobLayers().at(layer-1);
256			break;
257
258			case StripSubdetector::TID:
259			//edm::LogInfo(TkDetLayers) << "TID wheel n: " << TIDDetId(id).wheel() ;
260			if( side == -1 ) {
261			return trackerGeoSearch_->negTidLayers().at(layer-1);
262			}else if( side == 1 ) {
263			return trackerGeoSearch_->posTidLayers().at(layer-1);
264			}
265			break;
266
267			case StripSubdetector::TEC:
268			//edm::LogInfo(TkDetLayers) << "TEC wheel n: " << TECDetId(id).wheel() ;
269			if( side == -1 ) {
270			return trackerGeoSearch_->negTecLayers().at(layer-1);
271			} else if( side == 1 ) {
272			return trackerGeoSearch_->posTecLayers().at(layer-1);
273			}
274			break;
275	bendavid	1.8
276	loizides	1.9	case PixelSubdetector::PixelBarrel:
277			//edm::LogInfo(TkDetLayers) << "PixelBarrel layer n: " << PXBDetId(id).layer() ;
278			return trackerGeoSearch_->pixelBarrelLayers().at(layer-1);
279			break;
280	bendavid	1.8
281	loizides	1.9	case PixelSubdetector::PixelEndcap:
282			//edm::LogInfo(TkDetLayers) << "PixelEndcap disk n: " << PXFDetId(id).disk() ;
283			if(side == -1 ) {
284			return trackerGeoSearch_->negPixelForwardLayers().at(layer-1);
285			}else if( side == 1 ) {
286			return trackerGeoSearch_->posPixelForwardLayers().at(layer-1);
287			}
288			break;
289	bendavid	1.8
290	loizides	1.9	default:
291			return 0;
292			// throw(something);
293	bendavid	1.8	}
294			return 0; //just to avoid compile warnings
295			}
296
297	loizides	1.9	//--------------------------------------------------------------------------------------------------
298	bendavid	1.8	DetId HitDropper::StereoDetId(const DetId &i) const
299			{
300	loizides	1.9	// Deal with stereo det id.
301
302	bendavid	1.8	switch (i.det()) {
303	loizides	1.9	case DetId::Tracker:
304	bendavid	1.8	switch (i.subdetId()) {
305	loizides	1.9	case PixelSubdetector::PixelBarrel:
306			case PixelSubdetector::PixelEndcap:
307			return 0;
308			case StripSubdetector::TIB:
309			{
310			TIBDetId id = i;
311			if (id.isStereo())
312			return id;
313			else {
314			const std::vector<unsigned int> idString = id.string();
315			return TIBDetId(id.layer(),idString[0],idString[1],idString[2],id.module(),1);
316			}
317			}
318			case StripSubdetector::TID:
319			{
320			TIDDetId id = i;
321			if (id.isStereo())
322			return id;
323			else {
324			const std::vector<unsigned int> idModule = id.module();
325			return TIDDetId(id.side(),id.wheel(),id.ring(),idModule[0],idModule[1],1);
326			}
327			}
328	bendavid	1.8	case StripSubdetector::TOB:
329			{
330	loizides	1.9	TOBDetId id = i;
331			if (id.isStereo())
332			return id;
333			else {
334			const std::vector<unsigned int> idRod = id.rod();
335			return TOBDetId(id.layer(),idRod[0],idRod[1],id.module(),1);
336			}
337	bendavid	1.8	}
338	loizides	1.9	case StripSubdetector::TEC:
339			{
340	bendavid	1.8	TECDetId id = i;
341			if (id.isStereo())
342			return id;
343			else {
344			const std::vector<unsigned int> idPetal = id.petal();
345			return TECDetId(id.side(),id.wheel(),idPetal[0],idPetal[1],id.ring(),id.module(),1);
346			}
347	loizides	1.9	}
348			default:
349			return 0;
350	bendavid	1.8	}
351			break;
352	loizides	1.9	default:
353	bendavid	1.8	return 0;
354	bendavid	1.5	}
355			}