Producers/src/HitDropper.cc

// $Id: HitDropper.cc,v 1.3 2009/02/06 16:08:44 mrudolph 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"

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 crossing(HelixPlaneCrossing::PositionType(vtxTSCP.theState().position()),
                                         HelixPlaneCrossing::DirectionType(vtxTSCP.theState().momentum()),
                                         vtxTSCP.theState().transverseCurvature(),
                                         anyDirection);                                        
                                            
    std::pair<bool,double> crossResult = crossing.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 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 (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;
}
Revision:	1.4
Committed:	Fri Mar 20 18:01:48 2009 UTC (16 years, 1 month ago) by loizides
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_009c, Mit_009b, Mit_009a, Mit_009, Mit_008
Changes since 1.3:	+18 -21 lines
Log Message:	Cleanup
#	User	Rev	Content
1	loizides	1.4	// $Id: HitDropper.cc,v 1.3 2009/02/06 16:08:44 mrudolph 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.1
9			using namespace edm;
10			using namespace mitedm;
11
12			//--------------------------------------------------------------------------------------------------
13			reco::HitPattern HitDropper::CorrectedHits(const reco::TransientTrack *tTrack,
14	loizides	1.4	const ThreeVector &vtxPos) const
15	bendavid	1.1	{
16	loizides	1.4	// Return reco::HitPattern structure for the given track with all hits occuring before vtxPos
17			// on the track (relative to the primary vertex if given) removed.
18			// This version of the function uses an iterative helix-plane intersector and does not (yet)
19			// take into account the uncertainties in vertex position.
20	bendavid	1.2
21			const GlobalPoint vtxPoint(vtxPos.x(),vtxPos.y(),vtxPos.z());
22			const TrajectoryStateClosestToPoint vtxTSCP = tTrack->trajectoryStateClosestToPoint(vtxPoint);
23
24	bendavid	1.1	reco::HitPattern hitPattern;
25			int nHits = 0;
26			for (uint hi=0; hi<tTrack->recHitsSize(); ++hi) {
27			const TrackingRecHit *hit = tTrack->recHit(hi).get();
28	mrudolph	1.3	DetId geoId = hit->geographicalId();
29			if(geoId == uint32_t(0)) continue;
30			const GeomDet *det = trackerGeo_->idToDet(geoId);
31	bendavid	1.2
32			HelixArbitraryPlaneCrossing crossing(HelixPlaneCrossing::PositionType(vtxTSCP.theState().position()),
33			HelixPlaneCrossing::DirectionType(vtxTSCP.theState().momentum()),
34			vtxTSCP.theState().transverseCurvature(),
35			anyDirection);
36
37			std::pair<bool,double> crossResult = crossing.pathLength(det->surface());
38
39			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			const StraightLinePlaneCrossing::PositionType crossPosition = crossing.position(det->surface()).second;
87			const ThreeVector crossPos(crossPosition.x(), crossPosition.y(), crossPosition.z());
88			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			//add the hit only if it is after the vertex, allowing for some uncertainty in the vertex position
109			if ( lengthOverSigma>(-sigmaTolerance) ) {
110			hitPattern.set(*hit,nHits);
111			nHits++;
112			}
113			}
114			return hitPattern;
115	bendavid	1.1	}