ViewVC Help
View File | Revision Log | Show Annotations | Root Listing
root/cvsroot/UserCode/MitEdm/Producers/src/HitDropper.cc
Revision: 1.12
Committed: Tue Jun 8 20:17:30 2010 UTC (14 years, 11 months ago) by bendavid
Content type: text/plain
Branch: MAIN
CVS Tags: Mit_032, Mit_031, Mit_025c_branch2, Mit_025c_branch1, Mit_030, Mit_029c, Mit_029b, Mit_030_pre1, Mit_029a, Mit_029, Mit_029_pre1, Mit_028a, Mit_025c_branch0, Mit_028, Mit_027a, Mit_027, Mit_026, Mit_025e, Mit_025d, Mit_025c, Mit_025b, Mit_025a, Mit_025, Mit_025pre2, Mit_024b, Mit_025pre1, Mit_024a, Mit_024, Mit_023, Mit_022a, Mit_022, Mit_020d, TMit_020d, Mit_020c, Mit_021, Mit_021pre2, Mit_021pre1, Mit_020b, Mit_020a, Mit_020, Mit_020pre1, Mit_018, Mit_017, Mit_017pre3, Mit_017pre2, Mit_017pre1, V07-05-00, Mit_016, Mit_015b, Mit_015a, Mit_015, Mit_014e, Mit_014d, Mit_014c, Mit_014b, ConvRejection-10-06-09, HEAD
Branch point for: Mit_025c_branch
Changes since 1.11: +13 -4 lines
Log Message:
Add simple tool to do conversion selection and matching

File Contents

# User Rev Content
1 bendavid 1.12 // $Id: HitDropper.cc,v 1.11 2009/12/15 23:27:34 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 bendavid 1.12 std::pair<reco::HitPattern,uint> HitDropper::CorrectedHitsAOD(const reco::Track *track,
125 bendavid 1.5 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 bendavid 1.12 uint nWrongHits = 0;
155    
156 bendavid 1.10 for (Int_t hp=0; hp<3; ++hp) {
157     const reco::HitPattern &inhits = *inHitPatterns[hp];
158     for (Int_t hi=0; hi<inhits.numberOfHits(); hi++) {
159     uint32_t hit = inhits.getHitPattern(hi);
160     uint32_t subdet = inhits.getSubStructure(hit);
161     uint32_t layerid = inhits.getLayer(hit);
162     const DetLayer *det = FindLayer(subdet,layerid,side);
163     if (!det) continue;
164    
165     bool goodHit = false;
166    
167     //calculate distance between vertex and approximate hit position, projected into
168     //the appropriate plane/axis and compared to the uncertainty in vertex position
169     //in that plane/axis
170     if (det->location()==GeomDetEnumerators::barrel) {
171     const BarrelDetLayer *barrelDet = static_cast<const BarrelDetLayer*>(det);
172    
173     const double barrelRho = barrelDet->specificSurface().radius();
174     const double barrelHalfThickness = barrelDet->specificSurface().bounds().thickness()/2.0;
175     const ThreeVector crossPos(barrelRho*cos(track->phi()),barrelRho*sin(track->phi()),0.0);
176 bendavid 1.8
177 bendavid 1.10 const ThreeVector delta = crossPos - vtxPos;
178     const ThreeVector trkMomXY(trkMom.x(), trkMom.y(), 0.0);
179     const double deltaXY = delta.Dot(trkMomXY)/trkMomXY.R();
180     goodHit = deltaXY>(-sigmaTolerance*lxyError - barrelHalfThickness);
181     }
182     else if (det->location()==GeomDetEnumerators::endcap) {
183     const ForwardDetLayer *forwardDet = static_cast<const ForwardDetLayer*>(det);
184     const double diskZ = forwardDet->specificSurface().position().z();
185     const double diskHalfThickness = forwardDet->specificSurface().bounds().thickness()/2.0;
186     const double deltaZ = (diskZ - vtxPos.z())*trkMom.z()/fabs(trkMom.z());;
187     goodHit = deltaZ>(-sigmaTolerance*lzError - diskHalfThickness);
188     }
189     else
190     throw edm::Exception(edm::errors::Configuration, "HitDropper::CorrectedHitsAOD\n")
191     << "Error! Detector element not in a valid barrel or disk layer." << std::endl;
192 bendavid 1.5
193 bendavid 1.10 //add the hit only if it is after the vertex,
194     //allowing for some uncertainty in the vertex position
195     if ( goodHit ) {
196     bool isStereo = inhits.getSide(hit);
197     DetId dummyid = det->basicComponents().front()->geographicalId();
198     if (isStereo) {
199     dummyid = StereoDetId(dummyid);
200     }
201    
202     const TrackingRecHit::Type hitType = static_cast<TrackingRecHit::Type>(inhits.getHitType(hit));
203     InvalidTrackingRecHit dummyhit(dummyid, hitType);
204     hitPattern.set(dummyhit,nHits);
205     if ( hit != hitPattern.getHitPattern(nHits) )
206     throw edm::Exception(edm::errors::Configuration, "HitDropper::CorrectedHitsAOD\n")
207     << "Error! Mismatch in copying hit pattern." << std::endl;
208    
209     nHits++;
210 bendavid 1.12 }
211     else {
212     if (inhits.validHitFilter(hit)) {
213     if (inhits.trackerHitFilter(hit)) {
214     nWrongHits++;
215     }
216     }
217     }
218 bendavid 1.8 }
219     }
220    
221 bendavid 1.12 return std::pair<reco::HitPattern,uint>(hitPattern,nWrongHits);
222 bendavid 1.8 }
223    
224 bendavid 1.11 reco::HitPattern HitDropper::SharedHits(const reco::Track *t1, const reco::Track *t2) const
225     {
226     //Return a hit pattern corresponding to the hits on the two tracks which share clusters
227    
228     int nHits = 0;
229     reco::HitPattern sharedHits;
230    
231     if (!t1->extra().isAvailable() || !t2->extra().isAvailable())
232     return sharedHits;
233    
234     for (trackingRecHit_iterator iHit1 = t1->recHitsBegin(); iHit1 != t1->recHitsEnd(); ++iHit1) {
235     const TrackingRecHit *hit1 = iHit1->get();
236     if (hit1->isValid()) {
237     for (trackingRecHit_iterator iHit2 = t2->recHitsBegin(); iHit2 != t2->recHitsEnd(); ++iHit2) {
238     const TrackingRecHit *hit2 = iHit2->get();
239     if (hit2->isValid() && hit1->sharesInput(hit2,TrackingRecHit::some)) {
240     sharedHits.set(*hit1,nHits);
241     nHits++;
242     }
243     }
244     }
245     }
246    
247     return sharedHits;
248    
249     }
250    
251 bendavid 1.8 //--------------------------------------------------------------------------------------------------
252     const DetLayer *HitDropper::FindLayer(int subdet, int layer, int side) const
253     {
254 loizides 1.9 // Find a given layer.
255    
256 bendavid 1.8 switch(subdet) {
257 loizides 1.9 case StripSubdetector::TIB:
258     //edm::LogInfo(TkDetLayers) << "TIB layer n: " << TIBDetId(id).layer() ;
259     return trackerGeoSearch_->tibLayers().at(layer-1);
260     break;
261    
262     case StripSubdetector::TOB:
263     //edm::LogInfo(TkDetLayers) << "TOB layer n: " << TOBDetId(id).layer() ;
264     return trackerGeoSearch_->tobLayers().at(layer-1);
265     break;
266    
267     case StripSubdetector::TID:
268     //edm::LogInfo(TkDetLayers) << "TID wheel n: " << TIDDetId(id).wheel() ;
269     if( side == -1 ) {
270     return trackerGeoSearch_->negTidLayers().at(layer-1);
271     }else if( side == 1 ) {
272     return trackerGeoSearch_->posTidLayers().at(layer-1);
273     }
274     break;
275    
276     case StripSubdetector::TEC:
277     //edm::LogInfo(TkDetLayers) << "TEC wheel n: " << TECDetId(id).wheel() ;
278     if( side == -1 ) {
279     return trackerGeoSearch_->negTecLayers().at(layer-1);
280     } else if( side == 1 ) {
281     return trackerGeoSearch_->posTecLayers().at(layer-1);
282     }
283     break;
284 bendavid 1.8
285 loizides 1.9 case PixelSubdetector::PixelBarrel:
286     //edm::LogInfo(TkDetLayers) << "PixelBarrel layer n: " << PXBDetId(id).layer() ;
287     return trackerGeoSearch_->pixelBarrelLayers().at(layer-1);
288     break;
289 bendavid 1.8
290 loizides 1.9 case PixelSubdetector::PixelEndcap:
291     //edm::LogInfo(TkDetLayers) << "PixelEndcap disk n: " << PXFDetId(id).disk() ;
292     if(side == -1 ) {
293     return trackerGeoSearch_->negPixelForwardLayers().at(layer-1);
294     }else if( side == 1 ) {
295     return trackerGeoSearch_->posPixelForwardLayers().at(layer-1);
296     }
297     break;
298 bendavid 1.8
299 loizides 1.9 default:
300     return 0;
301     // throw(something);
302 bendavid 1.8 }
303     return 0; //just to avoid compile warnings
304     }
305    
306 loizides 1.9 //--------------------------------------------------------------------------------------------------
307 bendavid 1.8 DetId HitDropper::StereoDetId(const DetId &i) const
308     {
309 loizides 1.9 // Deal with stereo det id.
310    
311 bendavid 1.8 switch (i.det()) {
312 loizides 1.9 case DetId::Tracker:
313 bendavid 1.8 switch (i.subdetId()) {
314 loizides 1.9 case PixelSubdetector::PixelBarrel:
315     case PixelSubdetector::PixelEndcap:
316     return 0;
317     case StripSubdetector::TIB:
318     {
319     TIBDetId id = i;
320     if (id.isStereo())
321     return id;
322     else {
323     const std::vector<unsigned int> idString = id.string();
324     return TIBDetId(id.layer(),idString[0],idString[1],idString[2],id.module(),1);
325     }
326     }
327     case StripSubdetector::TID:
328     {
329     TIDDetId id = i;
330     if (id.isStereo())
331     return id;
332     else {
333     const std::vector<unsigned int> idModule = id.module();
334     return TIDDetId(id.side(),id.wheel(),id.ring(),idModule[0],idModule[1],1);
335     }
336     }
337 bendavid 1.8 case StripSubdetector::TOB:
338     {
339 loizides 1.9 TOBDetId id = i;
340     if (id.isStereo())
341     return id;
342     else {
343     const std::vector<unsigned int> idRod = id.rod();
344     return TOBDetId(id.layer(),idRod[0],idRod[1],id.module(),1);
345     }
346 bendavid 1.8 }
347 loizides 1.9 case StripSubdetector::TEC:
348     {
349 bendavid 1.8 TECDetId id = i;
350     if (id.isStereo())
351     return id;
352     else {
353     const std::vector<unsigned int> idPetal = id.petal();
354     return TECDetId(id.side(),id.wheel(),idPetal[0],idPetal[1],id.ring(),id.module(),1);
355     }
356 loizides 1.9 }
357     default:
358     return 0;
359 bendavid 1.8 }
360     break;
361 loizides 1.9 default:
362 bendavid 1.8 return 0;
363 bendavid 1.5 }
364     }