00001 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00002 #include "DataFormats/SiPixelDetId/interface/PXBDetId.h"
00003 #include "DataFormats/SiPixelDetId/interface/PXFDetId.h"
00004 #include "DataFormats/SiStripDetId/interface/TIBDetId.h"
00005 #include "DataFormats/SiStripDetId/interface/TIDDetId.h"
00006 #include "DataFormats/SiStripDetId/interface/TOBDetId.h"
00007 #include "DataFormats/SiStripDetId/interface/TECDetId.h"
00008
00009 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
00010 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00011 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
00012
00013 #include "RecoTracker/TrackProducer/interface/TrackMerger.h"
00014
00015 #include <boost/foreach.hpp>
00016 #define foreach BOOST_FOREACH
00017
00018 #define TRACK_SORT 1 // just use all hits from inner track, then append from outer outside it
00019 #define DET_SORT 0 // sort hits using global position of detector
00020 #define HIT_SORT 0 // sort hits using global position of hit
00021
00022
00023 TrackMerger::TrackMerger(const edm::ParameterSet &iConfig) :
00024 useInnermostState_(iConfig.getParameter<bool>("useInnermostState")),
00025 debug_(iConfig.getUntrackedParameter<bool>("debug",false)),
00026 theBuilderName(iConfig.getParameter<std::string>("ttrhBuilderName"))
00027 {
00028 }
00029
00030 TrackMerger::~TrackMerger()
00031 {
00032 }
00033
00034 void TrackMerger::init(const edm::EventSetup &iSetup)
00035 {
00036 iSetup.get<TrackerDigiGeometryRecord>().get(theGeometry);
00037 iSetup.get<IdealMagneticFieldRecord>().get(theMagField);
00038 iSetup.get<TransientRecHitRecord>().get(theBuilderName,theBuilder);
00039 }
00040
00041 TrackCandidate TrackMerger::merge(const reco::Track &inner, const reco::Track &outer) const
00042 {
00043 std::vector<const TrackingRecHit *> hits;
00044 hits.reserve(inner.recHitsSize() + outer.recHitsSize());
00045 if (debug_) std::cout << "Inner track hits: " << std::endl;
00046 for (trackingRecHit_iterator it = inner.recHitsBegin(), ed = inner.recHitsEnd(); it != ed; ++it) {
00047 hits.push_back(&**it);
00048 if (debug_) {
00049 DetId id(hits.back()->geographicalId());
00050 std::cout << " subdet " << id.subdetId() << " layer " << layer(id) << " valid " << hits.back()->isValid() << " detid: " << id() << std::endl;
00051 }
00052 }
00053 if (debug_) std::cout << "Outer track hits: " << std::endl;
00054
00055 #if TRACK_SORT
00056 DetId lastId(hits.back()->geographicalId());
00057 int lastSubdet = lastId.subdetId(), lastLayer = layer(lastId);
00058 for (trackingRecHit_iterator it = outer.recHitsBegin(), ed = outer.recHitsEnd(); it != ed; ++it) {
00059 const TrackingRecHit *hit = &**it;
00060 DetId id(hit->geographicalId());
00061 int thisSubdet = id.subdetId();
00062 if (thisSubdet > lastSubdet || (thisSubdet == lastSubdet && layer(id) > lastLayer)) {
00063 hits.push_back(hit);
00064 if (debug_) std::cout << " adding subdet " << id.subdetId() << " layer " << layer(id) << " valid " << hit->isValid() << " detid: " << id() << std::endl;
00065 } else {
00066 if (debug_) std::cout << " skipping subdet " << thisSubdet << " layer " << layer(id) << " valid " << hit->isValid() << " detid: " << id() << std::endl;
00067 }
00068 }
00069 #else
00070 size_t nHitsFirstTrack = hits.size();
00071 for (trackingRecHit_iterator it = outer.recHitsBegin(), ed = outer.recHitsEnd(); it != ed; ++it) {
00072 const TrackingRecHit *hit = &**it;
00073 DetId id(hit->geographicalId());
00074 int lay = layer(id);
00075 bool shared = false;
00076 bool valid = hit->isValid();
00077 if (debug_) std::cout << " subdet " << id.subdetId() << " layer " << layer(id) << " valid " << valid << " detid: " << id() << std::endl;
00078 size_t iHit = 0;
00079 foreach(const TrackingRecHit *& hit2, hits) {
00080 ++iHit; if (iHit > nHitsFirstTrack) break;
00081 DetId id2 = hit2->geographicalId();
00082 if (id.subdetId() != id2.subdetId()) continue;
00083 if (layer(id2) != lay) continue;
00084 if (hit->sharesInput(hit2, TrackingRecHit::all)) {
00085 if (debug_) std::cout << " discared as duplicate of other hit" << id() << std::endl;
00086 shared = true; break;
00087 }
00088 if (hit2->isValid() && !valid) {
00089 if (debug_) std::cout << " replacing old invalid hit on detid " << id2() << std::endl;
00090 hit2 = hit; shared = true; break;
00091 }
00092 if (debug_) std::cout << " discared as additional hit on layer that already contains hit with detid " << id() << std::endl;
00093 shared = true; break;
00094 }
00095 if (shared) continue;
00096 hits.push_back(hit);
00097 }
00098 #endif
00099
00100 math::XYZVector p = (inner.innerMomentum() + outer.outerMomentum());
00101 GlobalVector v(p.x(), p.y(), p.z());
00102 if (!useInnermostState_) v *= -1;
00103
00104 TrackCandidate::RecHitContainer ownHits;
00105 unsigned int nhits = hits.size();
00106 ownHits.reserve(nhits);
00107
00108 #if TRACK_SORT
00109 if (!useInnermostState_) std::reverse(hits.begin(), hits.end());
00110 foreach(const TrackingRecHit * hit, hits) {
00111 ownHits.push_back(*hit);
00112 }
00113 #elif DET_SORT
00114
00115 std::sort(hits.begin(), hits.end(), MomentumSort(v, &*theGeometry));
00116
00117 foreach(const TrackingRecHit * hit, hits) {
00118 ownHits.push_back(*hit);
00119 }
00120 #else
00121
00122 std::vector<TransientTrackingRecHit::RecHitPointer> ttrh(nhits);
00123 for (unsigned int i = 0; i < nhits; ++i) ttrh[i] = theBuilder->build(hits[i]);
00124 std::sort(ttrh.begin(), ttrh.end(), GlobalMomentumSort(v));
00125
00126 foreach(const TransientTrackingRecHit::RecHitPointer & hit, ttrh) {
00127 ownHits.push_back(*hit->hit());
00128 }
00129 #endif
00130
00131 PTrajectoryStateOnDet state;
00132 PropagationDirection pdir;
00133 if (useInnermostState_) {
00134 pdir = alongMomentum;
00135 if ((inner.outerPosition()-inner.innerPosition()).Dot(inner.momentum()) >= 0) {
00136
00137 TrajectoryStateOnSurface originalTsosIn(trajectoryStateTransform::innerStateOnSurface(inner, *theGeometry, &*theMagField));
00138 state = trajectoryStateTransform::persistentState( originalTsosIn, DetId(inner.innerDetId()) );
00139 } else {
00140
00141 TrajectoryStateOnSurface originalTsosOut(trajectoryStateTransform::outerStateOnSurface(inner, *theGeometry, &*theMagField));
00142 state = trajectoryStateTransform::persistentState( originalTsosOut, DetId(inner.outerDetId()) );
00143 }
00144 } else {
00145 pdir = oppositeToMomentum;
00146 if ((outer.outerPosition()-inner.innerPosition()).Dot(inner.momentum()) >= 0) {
00147
00148 TrajectoryStateOnSurface originalTsosOut(trajectoryStateTransform::outerStateOnSurface(outer, *theGeometry, &*theMagField));
00149 state = trajectoryStateTransform::persistentState( originalTsosOut, DetId(outer.outerDetId()) );
00150 } else {
00151
00152 TrajectoryStateOnSurface originalTsosIn(trajectoryStateTransform::innerStateOnSurface(outer, *theGeometry, &*theMagField));
00153 state = trajectoryStateTransform::persistentState( originalTsosIn, DetId(outer.innerDetId()) );
00154 }
00155
00156 }
00157 TrajectorySeed seed(state, TrackCandidate::RecHitContainer(), pdir);
00158 return TrackCandidate(ownHits, seed, state, (useInnermostState_ ? inner : outer).seedRef());
00159 }
00160
00161 int TrackMerger::layer(DetId detid) const {
00162 switch (detid.subdetId()) {
00163 case PixelSubdetector::PixelBarrel: return PXBDetId(detid).layer();
00164 case PixelSubdetector::PixelEndcap: return PXFDetId(detid).disk();
00165 case StripSubdetector::TIB: return TIBDetId(detid).layer();
00166 case StripSubdetector::TID: return TIDDetId(detid).wheel();
00167 case StripSubdetector::TOB: return TOBDetId(detid).layer();
00168 case StripSubdetector::TEC: return TECDetId(detid).wheel();
00169 }
00170 return -1;
00171 }
00172
00173 bool TrackMerger::MomentumSort::operator()(const TrackingRecHit *hit1, const TrackingRecHit *hit2) const
00174 {
00175 const GeomDet *det1 = geom_->idToDet(hit1->geographicalId());
00176 const GeomDet *det2 = geom_->idToDet(hit2->geographicalId());
00177 GlobalPoint p1 = det1->toGlobal(LocalPoint(0,0,0));
00178 GlobalPoint p2 = det2->toGlobal(LocalPoint(0,0,0));
00179 return (p2 - p1).dot(dir_) > 0;
00180 }
00181
00182 bool TrackMerger::GlobalMomentumSort::operator()(const TransientTrackingRecHit::RecHitPointer &hit1, const TransientTrackingRecHit::RecHitPointer &hit2) const
00183 {
00184 GlobalPoint p1 = hit1->isValid() ? hit1->globalPosition() : hit1->det()->position();
00185 GlobalPoint p2 = hit2->isValid() ? hit2->globalPosition() : hit2->det()->position();
00186 return (p2 - p1).dot(dir_) > 0;
00187 }