Producers/src/HitDropper.cc

// $Id: HitDropper.cc,v 1.6 2009/07/12 20:53:07 bendavid Exp $

#include "MitEdm/Producers/interface/HitDropper.h"
#include "DataFormats/TrackingRecHit/interface/InvalidTrackingRecHit.h"
#include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
#include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
#include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
#include "TrackingTools/GeomPropagators/interface/StraightLineBarrelCylinderCrossing.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 StraightLinePlaneCrossing::PositionType vtxPosition(vtxPos);
  const StraightLinePlaneCrossing::DirectionType trkMomentum(trkMom);
  
  StraightLinePlaneCrossing crossing(vtxPosition,trkMomentum,anyDirection);
  
  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);
    
    //calculate intersection of straight line with plane
    const StraightLinePlaneCrossing::PositionType cPos = crossing.position(det->surface()).second;
    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.

  // *** This function is not working yet ***
  
  return reco::HitPattern();
  
  if (0) {
  
    const StraightLinePlaneCrossing::PositionType vtxPosition(vtxPos);
    const StraightLinePlaneCrossing::DirectionType trkMomentum(trkMom);
    
    const GlobalPoint cVtxPosition(vtxPos.x(), vtxPos.y(), vtxPos.z());
    const GlobalVector cVtxMomentum(trkMom.x(), trkMom.y(), trkMom.z());
    
    StraightLinePlaneCrossing planeCrossing(vtxPosition,trkMomentum,anyDirection);
    StraightLineBarrelCylinderCrossing cylinderCrossing(cVtxPosition,cVtxMomentum,anyDirection);
    
    reco::HitPattern hitPattern;
    int nHits = 0;
    const reco::HitPattern &inhits = track->hitPattern();
    for (Int_t hi=0; hi<inhits.numberOfHits(); hi++) {
      uint32_t hit = inhits.getHitPattern(hi);
      uint32_t layerid = inhits.getLayer(hit);
      //if(layerid == uint32_t(0)) continue;
      //const GeomDet *detg = trackerGeo_->idToDet(layerid);
      //if (!detg) continue;
      //const DetLayer *det = detg->layer();
      const DetLayer *det = trackerGeoSearch_->detLayer(DetId(layerid));
      if (!det) continue;
      
      //calculate intersection of straight line with plane
      //const StraightLinePlaneCrossing::PositionType crossPosition = 
      //  crossing.position(det->surface()).second;
      //const ThreeVector crossPos(crossPosition.x(), crossPosition.y(), crossPosition.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 (det->location()==GeomDetEnumerators::barrel) {
        const BarrelDetLayer *barrelDet = static_cast<const BarrelDetLayer*>(det);
        double pathLength = cylinderCrossing.pathLength(barrelDet->specificSurface()).second;
        const GlobalPoint crossPosition = cylinderCrossing.position(pathLength);
        const ThreeVector crossPos(crossPosition.x(), crossPosition.y(), crossPosition.z());
        const ThreeVector delta = crossPos - vtxPos;
        const ThreeVector trkMomXY(trkMom.x(), trkMom.y(), 0.0);
        Double_t deltaXY = delta.Dot(trkMomXY)/trkMomXY.R();
        lengthOverSigma = deltaXY/lxyError;
      }
      else if (det->location()==GeomDetEnumerators::endcap) {
        const ForwardDetLayer *forwardDet = static_cast<const ForwardDetLayer*>(det);
        const StraightLinePlaneCrossing::PositionType cPos = 
          planeCrossing.position(forwardDet->specificSurface()).second;
        const ThreeVector crossPos(cPos.x(), cPos.y(), cPos.z());
        const ThreeVector delta = crossPos - vtxPos;
        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) ) {
        TrackingRecHit::Type hitType = static_cast<TrackingRecHit::Type>(inhits.getHitType(hit));
        InvalidTrackingRecHit dummyhit(layerid, hitType);
        hitPattern.set(dummyhit,nHits);
        printf("inhit = %i, outhit = %i\n",hit,hitPattern.getHitPattern(nHits));
        nHits++;
      } 
    }
    return hitPattern;
  }
}
Revision:	1.7
Committed:	Wed Jul 15 20:38:24 2009 UTC (15 years, 9 months ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_011, Mit_010a, Mit_010
Changes since 1.6:	+20 -18 lines
Log Message:	Cleanup
#	User	Rev	Content
1	loizides	1.7	// $Id: HitDropper.cc,v 1.6 2009/07/12 20:53:07 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			#include "TrackingTools/GeomPropagators/interface/HelixArbitraryPlaneCrossing.h"
7			#include "TrackingTools/GeomPropagators/interface/StraightLinePlaneCrossing.h"
8	bendavid	1.5	#include "TrackingTools/GeomPropagators/interface/StraightLineBarrelCylinderCrossing.h"
9	bendavid	1.1
10			using namespace edm;
11			using namespace mitedm;
12
13			//--------------------------------------------------------------------------------------------------
14			reco::HitPattern HitDropper::CorrectedHits(const reco::TransientTrack *tTrack,
15	loizides	1.4	const ThreeVector &vtxPos) const
16	bendavid	1.1	{
17	loizides	1.4	// Return reco::HitPattern structure for the given track with all hits occuring before vtxPos
18			// on the track (relative to the primary vertex if given) removed.
19			// This version of the function uses an iterative helix-plane intersector and does not (yet)
20			// take into account the uncertainties in vertex position.
21	bendavid	1.2
22			const GlobalPoint vtxPoint(vtxPos.x(),vtxPos.y(),vtxPos.z());
23			const TrajectoryStateClosestToPoint vtxTSCP = tTrack->trajectoryStateClosestToPoint(vtxPoint);
24
25	bendavid	1.1	reco::HitPattern hitPattern;
26			int nHits = 0;
27			for (uint hi=0; hi<tTrack->recHitsSize(); ++hi) {
28			const TrackingRecHit *hit = tTrack->recHit(hi).get();
29	mrudolph	1.3	DetId geoId = hit->geographicalId();
30			if(geoId == uint32_t(0)) continue;
31			const GeomDet *det = trackerGeo_->idToDet(geoId);
32	bendavid	1.2
33	loizides	1.7	HelixArbitraryPlaneCrossing c(HelixPlaneCrossing::PositionType(vtxTSCP.theState().position()),
34			HelixPlaneCrossing::DirectionType(vtxTSCP.theState().momentum()),
35			vtxTSCP.theState().transverseCurvature(),
36			anyDirection);
37	bendavid	1.2
38	loizides	1.7	std::pair<bool,double> crossResult = c.pathLength(det->surface());
39	bendavid	1.2	if ( crossResult.first && crossResult.second >= 0 ) {
40	bendavid	1.1	hitPattern.set(*hit,nHits);
41			nHits++;
42			}
43			}
44	bendavid	1.2
45	bendavid	1.1	return hitPattern;
46			}
47
48			//--------------------------------------------------------------------------------------------------
49			reco::HitPattern HitDropper::CorrectedHits(const reco::Track *track,
50	loizides	1.4	const ThreeVector &vtxPos) const
51	bendavid	1.1	{
52	loizides	1.4	// Build the transient track and then return the corrected HitPattern.
53	bendavid	1.1
54			reco::TransientTrack tTrack = builder_->build(track);
55	bendavid	1.2	return CorrectedHits(&tTrack, vtxPos);
56			}
57
58			//--------------------------------------------------------------------------------------------------
59			reco::HitPattern HitDropper::CorrectedHits(const reco::Track *track,
60	loizides	1.4	const ThreeVector &vtxPos,
61			const ThreeVector &trkMom,
62			Double_t lxyError,
63			Double_t lzError,
64			Double_t sigmaTolerance) const
65	bendavid	1.2	{
66	loizides	1.4	// Return reco::HitPattern structure for the given track with all hits occuring before vtxPos
67			// on the track (relative to the primary vertex if given) removed.
68			// This version of the function determines this completely analytically, and taking the
69			// vertex position uncertainty into account, which might be important for particles which decay
70			// within a tracker layer.
71	bendavid	1.2
72			const StraightLinePlaneCrossing::PositionType vtxPosition(vtxPos);
73			const StraightLinePlaneCrossing::DirectionType trkMomentum(trkMom);
74
75			StraightLinePlaneCrossing crossing(vtxPosition,trkMomentum,anyDirection);
76
77			reco::HitPattern hitPattern;
78			int nHits = 0;
79			for (uint hi=0; hi<track->recHitsSize(); ++hi) {
80			const TrackingRecHit *hit = track->recHit(hi).get();
81	mrudolph	1.3	DetId geoId = hit->geographicalId();
82			if(geoId == uint32_t(0)) continue;
83			const GeomDet *det = trackerGeo_->idToDet(geoId);
84	bendavid	1.2
85			//calculate intersection of straight line with plane
86	loizides	1.7	const StraightLinePlaneCrossing::PositionType cPos = crossing.position(det->surface()).second;
87			const ThreeVector crossPos(cPos.x(), cPos.y(), cPos.z());
88	bendavid	1.2	const ThreeVector delta = crossPos - vtxPos;
89
90			Double_t lengthOverSigma = 0;
91
92			//calculate distance between vertex and approximate hit position, projected into
93			//the appropriate plane/axis and compared to the uncertainty in vertex position
94			//in that plane/axis
95			if (IsBarrel(det)) {
96			const ThreeVector trkMomXY(trkMom.x(), trkMom.y(), 0.0);
97			Double_t deltaXY = delta.Dot(trkMomXY)/trkMomXY.R();
98			lengthOverSigma = deltaXY/lxyError;
99			}
100			else if (IsDisk(det)) {
101			Double_t deltaZ = delta.z()*trkMom.z()/fabs(trkMom.z());
102			lengthOverSigma = deltaZ/lzError;
103			}
104			else
105			throw edm::Exception(edm::errors::Configuration, "HitDropper::CorrectedHits\n")
106			<< "Error! Detector element not in a valid barrel or disk layer." << std::endl;
107
108	loizides	1.7	//add the hit only if it is after the vertex,
109			//allowing for some uncertainty in the vertex position
110	bendavid	1.2	if ( lengthOverSigma>(-sigmaTolerance) ) {
111			hitPattern.set(*hit,nHits);
112			nHits++;
113			}
114			}
115			return hitPattern;
116	bendavid	1.1	}
117	bendavid	1.5
118
119			//--------------------------------------------------------------------------------------------------
120			reco::HitPattern HitDropper::CorrectedHitsAOD(const reco::Track *track,
121			const ThreeVector &vtxPos,
122			const ThreeVector &trkMom,
123			Double_t lxyError,
124			Double_t lzError,
125			Double_t sigmaTolerance) const
126			{
127			// Return reco::HitPattern structure for the given track with all hits occuring before vtxPos
128			// on the track (relative to the primary vertex if given) removed.
129			// This version of the function determines this completely analytically, and taking the
130			// vertex position uncertainty into account, which might be important for particles which decay
131			// within a tracker layer.
132
133	loizides	1.7	// * This function is not working yet *
134	bendavid	1.5
135			return reco::HitPattern();
136
137			if (0) {
138
139			const StraightLinePlaneCrossing::PositionType vtxPosition(vtxPos);
140			const StraightLinePlaneCrossing::DirectionType trkMomentum(trkMom);
141
142			const GlobalPoint cVtxPosition(vtxPos.x(), vtxPos.y(), vtxPos.z());
143			const GlobalVector cVtxMomentum(trkMom.x(), trkMom.y(), trkMom.z());
144
145			StraightLinePlaneCrossing planeCrossing(vtxPosition,trkMomentum,anyDirection);
146			StraightLineBarrelCylinderCrossing cylinderCrossing(cVtxPosition,cVtxMomentum,anyDirection);
147
148			reco::HitPattern hitPattern;
149			int nHits = 0;
150			const reco::HitPattern &inhits = track->hitPattern();
151			for (Int_t hi=0; hi<inhits.numberOfHits(); hi++) {
152			uint32_t hit = inhits.getHitPattern(hi);
153			uint32_t layerid = inhits.getLayer(hit);
154	bendavid	1.6	//if(layerid == uint32_t(0)) continue;
155			//const GeomDet *detg = trackerGeo_->idToDet(layerid);
156			//if (!detg) continue;
157			//const DetLayer *det = detg->layer();
158			const DetLayer *det = trackerGeoSearch_->detLayer(DetId(layerid));
159			if (!det) continue;
160	bendavid	1.5
161			//calculate intersection of straight line with plane
162	loizides	1.7	//const StraightLinePlaneCrossing::PositionType crossPosition =
163			// crossing.position(det->surface()).second;
164			//const ThreeVector crossPos(crossPosition.x(), crossPosition.y(), crossPosition.z());
165			//const ThreeVector delta = crossPos - vtxPos;
166	bendavid	1.5
167			Double_t lengthOverSigma = 0;
168
169			//calculate distance between vertex and approximate hit position, projected into
170			//the appropriate plane/axis and compared to the uncertainty in vertex position
171			//in that plane/axis
172			if (det->location()==GeomDetEnumerators::barrel) {
173			const BarrelDetLayer barrelDet = static_cast<const BarrelDetLayer>(det);
174			double pathLength = cylinderCrossing.pathLength(barrelDet->specificSurface()).second;
175			const GlobalPoint crossPosition = cylinderCrossing.position(pathLength);
176			const ThreeVector crossPos(crossPosition.x(), crossPosition.y(), crossPosition.z());
177			const ThreeVector delta = crossPos - vtxPos;
178			const ThreeVector trkMomXY(trkMom.x(), trkMom.y(), 0.0);
179			Double_t deltaXY = delta.Dot(trkMomXY)/trkMomXY.R();
180			lengthOverSigma = deltaXY/lxyError;
181			}
182			else if (det->location()==GeomDetEnumerators::endcap) {
183			const ForwardDetLayer forwardDet = static_cast<const ForwardDetLayer>(det);
184	loizides	1.7	const StraightLinePlaneCrossing::PositionType cPos =
185			planeCrossing.position(forwardDet->specificSurface()).second;
186			const ThreeVector crossPos(cPos.x(), cPos.y(), cPos.z());
187	bendavid	1.5	const ThreeVector delta = crossPos - vtxPos;
188			Double_t deltaZ = delta.z()*trkMom.z()/fabs(trkMom.z());
189			lengthOverSigma = deltaZ/lzError;
190			}
191			else
192			throw edm::Exception(edm::errors::Configuration, "HitDropper::CorrectedHits\n")
193			<< "Error! Detector element not in a valid barrel or disk layer." << std::endl;
194
195	loizides	1.7	//add the hit only if it is after the vertex,
196			//allowing for some uncertainty in the vertex position
197	bendavid	1.5	if ( lengthOverSigma>(-sigmaTolerance) ) {
198	bendavid	1.6	TrackingRecHit::Type hitType = static_cast<TrackingRecHit::Type>(inhits.getHitType(hit));
199			InvalidTrackingRecHit dummyhit(layerid, hitType);
200			hitPattern.set(dummyhit,nHits);
201			printf("inhit = %i, outhit = %i\n",hit,hitPattern.getHitPattern(nHits));
202	bendavid	1.5	nHits++;
203			}
204			}
205			return hitPattern;
206			}
207			}