CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/RecoTracker/SpecialSeedGenerators/src/SeedFromGenericPairOrTriplet.cc

Go to the documentation of this file.
00001 #include "RecoTracker/SpecialSeedGenerators/interface/SeedFromGenericPairOrTriplet.h"
00002 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
00003 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
00004 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
00005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
00006 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
00007 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
00008 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
00009 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
00010 
00011 SeedFromGenericPairOrTriplet::SeedFromGenericPairOrTriplet(const MagneticField* mf,
00012                                                            const TrackerGeometry* geom,
00013                                                            const TransientTrackingRecHitBuilder* builder,
00014                                                            const Propagator* propagatorAlong,
00015                                                            const Propagator* propagatorOpposite,
00016                                                            const std::vector<int>& charges,
00017                                                            bool momFromPSet,
00018                                                            double errorRescaling):theMagfield(mf), theTracker(geom), theBuilder(builder), thePropagatorAlong(propagatorAlong), thePropagatorOpposite(propagatorOpposite), theSetMomentum(momFromPSet), theCharges(charges), theErrorRescaling(errorRescaling){}
00019 
00020 
00021 std::vector<TrajectorySeed*> SeedFromGenericPairOrTriplet::seed(const SeedingHitSet& hits,
00022                                                                 const PropagationDirection& dir,
00023                                                                 const NavigationDirection&  seedDir,
00024                                                                 const edm::EventSetup& iSetup){
00025         std::vector<TrajectorySeed*> seeds;
00026         if (hits.size() == 3) {
00027                 if(!theSetMomentum){
00028                         TrajectorySeed* seed = seedFromTriplet(hits, dir, seedDir, iSetup);
00029                         if (seed) seeds.push_back(seed);
00030                 } else {
00031                         for (std::vector<int>::const_iterator iCh = theCharges.begin(); iCh != theCharges.end(); iCh++){
00032                                 TrajectorySeed* seed = seedFromTriplet(hits, dir, seedDir, iSetup, *iCh);
00033                                 if(seed) seeds.push_back(seed);
00034                         }
00035                 } 
00036                 
00037         } else if (hits.size() == 2){
00038                 if(!theSetMomentum){
00039                         TrajectorySeed* seed = seedFromPair(hits, dir, seedDir);
00040                         if (seed) seeds.push_back(seed);
00041                 } else {
00042                         for (std::vector<int>::const_iterator iCh = theCharges.begin(); iCh != theCharges.end(); iCh++){
00043                                 TrajectorySeed* seed = seedFromPair(hits, dir, seedDir, *iCh);
00044                                 if (seed) seeds.push_back(seed);
00045                         }
00046                 }
00047         } else {
00048                 throw cms::Exception("CombinatorialSeedGeneratorForCosmics") << " Wrong number of hits in Set: "
00049                                                         << hits.size() << ", should be 2 or 3 ";        
00050         }
00051         return seeds;
00052 }
00053 
00054 TrajectorySeed* SeedFromGenericPairOrTriplet::seedFromTriplet(const SeedingHitSet& hits,
00055                                                               const PropagationDirection& dir,
00056                                                               const NavigationDirection&  seedDir,
00057                                                               const edm::EventSetup& iSetup, int charge) const{
00058         if (hits.size() != 3) {
00059                 throw cms::Exception("CombinatorialSeedGeneratorForCosmics") <<
00060                         "call to SeedFromGenericPairOrTriplet::seedFromTriplet with " << hits.size() << " hits ";
00061         }
00062 
00063         const TrackingRecHit* innerHit  = hits[0]->hit();
00064         const TrackingRecHit* middleHit = hits[1]->hit();
00065         const TrackingRecHit* outerHit  = hits[2]->hit();
00066         GlobalPoint inner  = theTracker->idToDet(innerHit->geographicalId() )->surface().toGlobal(innerHit->localPosition() );
00067         GlobalPoint middle = theTracker->idToDet(middleHit->geographicalId())->surface().toGlobal(middleHit->localPosition());
00068         GlobalPoint outer  = theTracker->idToDet(outerHit->geographicalId() )->surface().toGlobal(outerHit->localPosition() );
00069         if (theSetMomentum){
00070                 LogDebug("SeedFromGenericPairOrTriplet") << 
00071                         "Using the following hits: outer(r, phi, theta) " << outerHit->geographicalId().subdetId()<< " (" << outer.perp() 
00072                         << ", " << outer.phi() 
00073                         << "," << outer.theta() 
00074                         << ")    " << middleHit->geographicalId().subdetId() << " middle (" 
00075                         << middle.perp() << ", " 
00076                         << middle.phi() << "," 
00077                         << middle.theta() << ")    " 
00078                         <<innerHit->geographicalId().subdetId() <<" inner (" 
00079                         << inner.perp() 
00080                         << ", " << inner.phi() 
00081                         << "," << inner.theta() <<")   "
00082                         << "   (x, y, z)   outer ("  
00083                         << inner.x() << ", " 
00084                         << inner.y() << ", " 
00085                         << inner.z() << ")    middle ("  
00086                         << middle.x() << ", " 
00087                         << middle.y() << ", " 
00088                         << middle.z() << ")";
00089                 if (!qualityFilter(hits)) return 0;
00090                 bool sdir = (seedDir == outsideIn);
00091                 SeedingHitSet newSet(sdir ? hits[1] : hits[0], sdir ? hits[2] : hits[1]); 
00092                 TrajectorySeed* seed = seedFromPair(newSet, dir, seedDir, charge);
00093                 return seed;
00094                                                 
00095         }
00096         GlobalPoint* firstPoint  = 0;
00097         GlobalPoint* secondPoint = 0;
00098         GlobalPoint* thirdPoint  = 0;
00099         int momentumSign         = 1;
00100         //const TrackingRecHit* firstHit  = 0;
00101         //const TrackingRecHit* secondHit = 0;
00102         //choose the prop dir and hit order accordingly to where the seed is made
00103         std::vector<const TrackingRecHit*> trHits;
00104         if (seedDir == outsideIn){
00105                 LogDebug("SeedFromGenericPairOrTriplet") 
00106                         << "Seed from outsideIn alongMomentum OR insideOut oppositeToMomentum";
00107                 firstPoint = &outer;
00108                 secondPoint = &middle;
00109                 thirdPoint = &inner;
00110                 trHits.push_back(outerHit);
00111                 trHits.push_back(middleHit);
00112                 //firstHit  = outerHit;
00113                 //secondHit = middleHit;
00114         } else {
00115                 LogDebug("SeedFromGenericPairOrTriplet") 
00116                         << "Seed from outsideIn oppositeToMomentum OR insideOut alongMomentum";
00117                 firstPoint = &inner;
00118                 secondPoint = &middle;
00119                 thirdPoint = &outer;
00120                 trHits.push_back(innerHit);
00121                 trHits.push_back(middleHit);
00122                 //firstHit  = innerHit;
00123                 //secondHit = middleHit;
00124         }
00125         if (dir == oppositeToMomentum) momentumSign = -1; 
00126         FastHelix helix(*thirdPoint, *secondPoint, *firstPoint, iSetup);
00127         FreeTrajectoryState originalStartingState = helix.stateAtVertex();
00128         LogDebug("SeedFromGenericPairOrTriplet") << "originalStartingState " << originalStartingState;
00129         /*GlobalTrajectoryParameters originalPar = originalStartingState.parameters();
00130         GlobalTrajectoryParameters newPar = GlobalTrajectoryParameters(originalPar.position(), 
00131                                                                        momentumSign*originalPar.momentum(),
00132                                                                        originalPar.charge(),
00133                                                                        &originalPar.magneticField());
00134         */
00135         
00136         /*FastCircle helix(*thirdPoint, *secondPoint, *firstPoint);
00137         GlobalTrajectoryParameters newPar = GlobalTrajectoryParameters(*secondPoint,
00138                                                                        momentumSign*originalPar.momentum(),
00139                                                                        originalPar.charge(),
00140                                                                        &originalPar.magneticField());*/
00141         //FreeTrajectoryState* startingState = new FreeTrajectoryState(newPar, initialError(trHits[0]));//trHits[1]));
00142         if (!qualityFilter(momentumSign*originalStartingState.momentum())) return 0;
00143         //TrajectorySeed* seed = buildSeed(startingState, trHits, dir);
00144         TrajectorySeed* seed = buildSeed(momentumSign*originalStartingState.momentum(),originalStartingState.charge(), trHits, dir);
00145         return seed;    
00146 }
00147 
00148 TrajectorySeed* SeedFromGenericPairOrTriplet::seedFromPair(const SeedingHitSet& hits,
00149                                                            const PropagationDirection& dir,
00150                                                            const NavigationDirection&  seedDir, int charge) const{
00151         if (hits.size() != 2) {
00152                 throw cms::Exception("CombinatorialSeedGeneratorForCosmics") <<
00153                         "call to SeedFromGenericPairOrTriplet::seedFromPair with " << hits.size() << " hits ";
00154         }
00155         const TrackingRecHit* innerHit = hits[0]->hit();
00156         const TrackingRecHit* outerHit = hits[1]->hit();
00157         GlobalPoint inner  = theTracker->idToDet(innerHit->geographicalId() )->surface().toGlobal(innerHit->localPosition() );
00158         GlobalPoint outer  = theTracker->idToDet(outerHit->geographicalId() )->surface().toGlobal(outerHit->localPosition() );
00159         LogDebug("SeedFromGenericPairOrTriplet") <<
00160                         "Using the following hits: outer(r, phi, theta) (" << outer.perp()
00161                         << ", " << outer.phi()
00162                         << "," << outer.theta()
00163                         <<")    inner ("
00164                         << inner.perp()
00165                         << ", " << inner.phi()
00166                         << "," << inner.theta() <<")";
00167         GlobalPoint* firstPoint  = 0;
00168         GlobalPoint* secondPoint = 0;
00169         int momentumSign         = 1;
00170         //const TrackingRecHit* firstHit  = 0;
00171         //const TrackingRecHit* secondHit = 0;
00172         std::vector<const TrackingRecHit*> trHits;
00173         //choose the prop dir and hit order accordingly to where the seed is made
00174         if (seedDir == outsideIn){
00175                 LogDebug("SeedFromGenericPairOrTriplet")
00176                         << "Seed from outsideIn alongMomentum OR insideOut oppositeToMomentum";
00177                 firstPoint = &outer;
00178                 secondPoint = &inner;
00179                 //firstHit  = outerHit;
00180                 //secondHit = innerHit;
00181                 trHits.push_back(outerHit);
00182                 trHits.push_back(innerHit);
00183         } else {
00184                 LogDebug("SeedFromGenericPairOrTriplet")
00185                         << "Seed from outsideIn oppositeToMomentum OR insideOut alongMomentum";
00186                 firstPoint = &inner;
00187                 secondPoint = &outer;
00188                 momentumSign = -1;
00189                 //firstHit  = innerHit;
00190                 //secondHit = outerHit;
00191                 trHits.push_back(innerHit);
00192                 trHits.push_back(outerHit);
00193         }
00194         if (dir == oppositeToMomentum) momentumSign = -1;
00195         GlobalVector momentum = momentumSign*theP*(*secondPoint-*firstPoint).unit();
00196         /*GlobalTrajectoryParameters gtp(*firstPoint, 
00197                                        momentum,
00198                                        charge,
00199                                        &(*theMagfield));*/
00200         //FreeTrajectoryState* startingState = new FreeTrajectoryState(gtp,initialError(trHits[0]));//trHits[1]));
00201         //if (!qualityFilter(startingState, hits)) return 0;
00202         //TrajectorySeed* seed = buildSeed(startingState, firstHit, dir);
00203         //TrajectorySeed* seed = buildSeed(startingState, trHits, dir);
00204         TrajectorySeed* seed = buildSeed(momentum, charge, trHits, dir);
00205         return seed;
00206 }
00207 
00208 
00209 TrajectorySeed* SeedFromGenericPairOrTriplet::buildSeed(const GlobalVector& momentum,
00210                                                         int charge, 
00211                                                         //const TrackingRecHit* firsthit,
00212                                                         std::vector<const TrackingRecHit*>& trHits,
00213                                                         const PropagationDirection& dir) const {
00214         const Propagator* propagator = thePropagatorAlong;
00215         if (dir == oppositeToMomentum) propagator = thePropagatorOpposite;
00216 
00217         //debug
00218         GlobalPoint first = theTracker->idToDet(trHits[0]->geographicalId() )->surface().toGlobal(trHits[0]->localPosition() );
00219         GlobalPoint second  = theTracker->idToDet(trHits[1]->geographicalId() )->surface().toGlobal(trHits[1]->localPosition() );
00220         LogDebug("SeedFromGenericPairOrTriplet") <<
00221                         "Using the following hits: first(r, phi, theta) (" << first.perp()
00222                         << ", " << first.phi()
00223                         << "," << first.theta()
00224                         <<")    second ("
00225                         << second.perp()
00226                         << ", " << second.phi()
00227                         << "," << second.theta() <<")";
00228         //build an initial trajectory state with big errors, 
00229         //then update it with the first hit
00230         TransientTrackingRecHit::RecHitPointer transHit = theBuilder->build(trHits[0]);
00231         LocalVector lmom = transHit->surface()->toLocal(momentum);
00232         LocalTrajectoryParameters ltp(charge/momentum.mag(),
00233                                       lmom.x()/lmom.z(),lmom.y()/lmom.z(),
00234                                       transHit->localPosition().x(),
00235                                       transHit->localPosition().y(),
00236                                       lmom.z()>0.?1.:-1.); 
00237 
00238         LocalTrajectoryError lte(100.,100.,
00239                                 (1.+(lmom.x()/lmom.z())*(lmom.x()/lmom.z())),
00240                                 (1.+(lmom.y()/lmom.z())*(lmom.y()/lmom.z())),
00241                                 1./momentum.mag());
00242         
00243         TrajectoryStateOnSurface initialState(ltp,lte,*transHit->surface(),&(*theMagfield));
00244         KFUpdator updator;
00245         TrajectoryStateOnSurface startingState = updator.update(initialState,*transHit);
00246 
00247         TransientTrackingRecHit::RecHitPointer transHit2 = theBuilder->build(trHits[1]);
00248 
00249         //TrajectoryStateOnSurface seedTSOS(*startingState, *plane);
00250         const TrajectoryStateOnSurface propTSOS = propagator->propagate(startingState,*(transHit2->surface()));
00251         if (!propTSOS.isValid()){
00252                 LogDebug("SeedFromGenericPairOrTriplet") << "first propagation failed";
00253                 return 0;
00254         }
00255         TrajectoryStateOnSurface seedTSOS = updator.update(propTSOS, *transHit2);
00256         if (!seedTSOS.isValid()){
00257                 LogDebug("SeedFromGenericPairOrTriplet") << "first update failed";
00258                 return 0;
00259         }       
00260         LogDebug("SeedFromGenericPairOrTriplet") << "starting TSOS " << seedTSOS ;
00261         seedTSOS.rescaleError(theErrorRescaling);
00262         PTrajectoryStateOnDet PTraj=
00263                 trajectoryStateTransform::persistentState(seedTSOS, trHits[1]->geographicalId().rawId());
00264         edm::OwnVector<TrackingRecHit> seed_hits;
00265         //build the transientTrackingRecHit for the starting hit, then call the method clone to rematch if needed
00266         std::vector<const TrackingRecHit*>::const_iterator ihits;
00267         for (ihits = trHits.begin(); ihits != trHits.end(); ihits++){
00268                 seed_hits.push_back((*ihits)->clone());
00269         }
00270         TrajectorySeed* seed = new TrajectorySeed(PTraj,seed_hits,dir);
00271         return seed;
00272 }
00273 /*
00274 CurvilinearTrajectoryError SeedFromGenericPairOrTriplet::initialError(const TrackingRecHit* rechit) {
00275         TransientTrackingRecHit::RecHitPointer transHit =  theBuilder->build(rechit);
00276         //AlgebraicSymMatrix C(5,1);
00277         AlgebraicSymMatrix55 C = AlgebraicMatrixID();
00278         //C*=0.1;
00279         if (theSetMomentum){
00280                 C(0,0)=1/theP;
00281                 C(3,3)=transHit->globalPositionError().cxx();
00282                 C(4,4)=transHit->globalPositionError().czz();
00283                 
00284                 //C[0][0]=1/theP;
00285                 //C[3][3]=transHit->globalPositionError().cxx();
00286                 //C[4][4]=transHit->globalPositionError().czz();
00287         } else {
00288                 float zErr = transHit->globalPositionError().czz();
00289                 float transverseErr = transHit->globalPositionError().cxx(); // assume equal cxx cyy
00290                 C(3,3) = transverseErr;
00291                 C(4,4) = zErr;
00292                 //C[3][3] = transverseErr;
00293                 //C[4][4] = zErr;
00294         }
00295 
00296         return CurvilinearTrajectoryError(C);
00297 }
00298 */
00299 bool SeedFromGenericPairOrTriplet::qualityFilter(const SeedingHitSet& hits) const{
00300         //if we are setting the momentum with the PSet we look for aligned hits
00301         if (theSetMomentum){
00302                 if (hits.size()==3){
00303                         std::vector<GlobalPoint> gPoints;
00304                         unsigned int nHits = hits.size();
00305                         for (unsigned int iHit = 0; iHit < nHits; ++iHit) gPoints.push_back(hits[iHit]->globalPosition());
00306                         unsigned int subid=(*hits[0]).geographicalId().subdetId();
00307                         if(subid == StripSubdetector::TEC || subid == StripSubdetector::TID){
00308                                 LogDebug("SeedFromGenericPairOrTriplet") 
00309                                         << "In the endcaps we cannot decide if hits are aligned with only phi and z";
00310                                 return true;
00311                         }
00312                         FastCircle circle(gPoints[0],
00313                                           gPoints[1],
00314                                           gPoints[2]);
00315                         if (circle.rho() < 500 && circle.rho() != 0) {
00316                                 LogDebug("SeedFromGenericPairOrTriplet") <<
00317                                         "Seed qualityFilter rejected because rho = " << circle.rho();
00318                                 return false;
00319                         }
00320 
00321                 }
00322         }
00323         return true;
00324 }
00325 
00326 bool SeedFromGenericPairOrTriplet::qualityFilter(const GlobalVector& momentum) const{
00327         if (!theSetMomentum){
00328                 if (momentum.perp() < theP) return false;
00329         }
00330         return true; 
00331 }