CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_6_1_2_SLHC2/src/RecoTracker/TrackProducer/src/TrackMerger.cc

Go to the documentation of this file.
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     // OLD METHOD, sometimes fails
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     // NEW sort, more accurate
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             // use inner state
00137             TrajectoryStateOnSurface originalTsosIn(trajectoryStateTransform::innerStateOnSurface(inner, *theGeometry, &*theMagField));
00138             state = trajectoryStateTransform::persistentState( originalTsosIn, DetId(inner.innerDetId()) );
00139         } else { 
00140             // use outer state
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             // use outer state
00148             TrajectoryStateOnSurface originalTsosOut(trajectoryStateTransform::outerStateOnSurface(outer, *theGeometry, &*theMagField));
00149             state = trajectoryStateTransform::persistentState( originalTsosOut, DetId(outer.outerDetId()) );
00150         } else { 
00151             // use inner state
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; // never match
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 }