Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:28:45 

0001 #include "RecoTracker/SpecialSeedGenerators/interface/SeedFromGenericPairOrTriplet.h"
0002 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0003 #include "DataFormats/GeometryVector/interface/GlobalPoint.h"
0004 #include "DataFormats/SiStripDetId/interface/StripSubdetector.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "Geometry/CommonDetUnit/interface/GeomDet.h"
0007 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateOnSurface.h"
0008 #include "TrackingTools/TransientTrackingRecHit/interface/TransientTrackingRecHit.h"
0009 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0010 
0011 SeedFromGenericPairOrTriplet::SeedFromGenericPairOrTriplet(const MagneticField* mf,
0012                                                            const TrackerGeometry* geom,
0013                                                            const TransientTrackingRecHitBuilder* builder,
0014                                                            const Propagator* propagatorAlong,
0015                                                            const Propagator* propagatorOpposite,
0016                                                            const std::vector<int>& charges,
0017                                                            bool momFromPSet,
0018                                                            double errorRescaling)
0019     : theMagfield(mf),
0020       theTracker(geom),
0021       theBuilder(builder),
0022       thePropagatorAlong(propagatorAlong),
0023       thePropagatorOpposite(propagatorOpposite),
0024       theSetMomentum(momFromPSet),
0025       theCharges(charges),
0026       theErrorRescaling(errorRescaling) {}
0027 
0028 std::vector<TrajectorySeed*> SeedFromGenericPairOrTriplet::seed(const SeedingHitSet& hits,
0029                                                                 const PropagationDirection& dir,
0030                                                                 const NavigationDirection& seedDir,
0031                                                                 const edm::EventSetup& iSetup) {
0032   std::vector<TrajectorySeed*> seeds;
0033   if (hits.size() == 3) {
0034     if (!theSetMomentum) {
0035       TrajectorySeed* seed = seedFromTriplet(hits, dir, seedDir, iSetup);
0036       if (seed)
0037         seeds.push_back(seed);
0038     } else {
0039       for (std::vector<int>::const_iterator iCh = theCharges.begin(); iCh != theCharges.end(); iCh++) {
0040         TrajectorySeed* seed = seedFromTriplet(hits, dir, seedDir, iSetup, *iCh);
0041         if (seed)
0042           seeds.push_back(seed);
0043       }
0044     }
0045 
0046   } else if (hits.size() == 2) {
0047     if (!theSetMomentum) {
0048       TrajectorySeed* seed = seedFromPair(hits, dir, seedDir);
0049       if (seed)
0050         seeds.push_back(seed);
0051     } else {
0052       for (std::vector<int>::const_iterator iCh = theCharges.begin(); iCh != theCharges.end(); iCh++) {
0053         TrajectorySeed* seed = seedFromPair(hits, dir, seedDir, *iCh);
0054         if (seed)
0055           seeds.push_back(seed);
0056       }
0057     }
0058   } else {
0059     throw cms::Exception("CombinatorialSeedGeneratorForCosmics")
0060         << " Wrong number of hits in Set: " << hits.size() << ", should be 2 or 3 ";
0061   }
0062   return seeds;
0063 }
0064 
0065 TrajectorySeed* SeedFromGenericPairOrTriplet::seedFromTriplet(const SeedingHitSet& hits,
0066                                                               const PropagationDirection& dir,
0067                                                               const NavigationDirection& seedDir,
0068                                                               const edm::EventSetup& iSetup,
0069                                                               int charge) const {
0070   if (hits.size() != 3) {
0071     throw cms::Exception("CombinatorialSeedGeneratorForCosmics")
0072         << "call to SeedFromGenericPairOrTriplet::seedFromTriplet with " << hits.size() << " hits ";
0073   }
0074 
0075   auto innerHit = hits[0];
0076   auto middleHit = hits[1];
0077   auto outerHit = hits[2];
0078   GlobalPoint inner = theTracker->idToDet(innerHit->geographicalId())->surface().toGlobal(innerHit->localPosition());
0079   GlobalPoint middle = theTracker->idToDet(middleHit->geographicalId())->surface().toGlobal(middleHit->localPosition());
0080   GlobalPoint outer = theTracker->idToDet(outerHit->geographicalId())->surface().toGlobal(outerHit->localPosition());
0081   if (theSetMomentum) {
0082     LogDebug("SeedFromGenericPairOrTriplet")
0083         << "Using the following hits: outer(r, phi, theta) " << outerHit->geographicalId().subdetId() << " ("
0084         << outer.perp() << ", " << outer.phi() << "," << outer.theta() << ")    "
0085         << middleHit->geographicalId().subdetId() << " middle (" << middle.perp() << ", " << middle.phi() << ","
0086         << middle.theta() << ")    " << innerHit->geographicalId().subdetId() << " inner (" << inner.perp() << ", "
0087         << inner.phi() << "," << inner.theta() << ")   "
0088         << "   (x, y, z)   outer (" << inner.x() << ", " << inner.y() << ", " << inner.z() << ")    middle ("
0089         << middle.x() << ", " << middle.y() << ", " << middle.z() << ")";
0090     if (!qualityFilter(hits))
0091       return nullptr;
0092     bool sdir = (seedDir == outsideIn);
0093     SeedingHitSet newSet(sdir ? hits[1] : hits[0], sdir ? hits[2] : hits[1]);
0094     TrajectorySeed* seed = seedFromPair(newSet, dir, seedDir, charge);
0095     return seed;
0096   }
0097   GlobalPoint* firstPoint = nullptr;
0098   GlobalPoint* secondPoint = nullptr;
0099   GlobalPoint* thirdPoint = nullptr;
0100   int momentumSign = 1;
0101   //const TrackingRecHit* firstHit  = 0;
0102   //const TrackingRecHit* secondHit = 0;
0103   //choose the prop dir and hit order accordingly to where the seed is made
0104   std::vector<const BaseTrackerRecHit*> trHits;
0105   if (seedDir == outsideIn) {
0106     LogDebug("SeedFromGenericPairOrTriplet") << "Seed from outsideIn alongMomentum OR insideOut oppositeToMomentum";
0107     firstPoint = &outer;
0108     secondPoint = &middle;
0109     thirdPoint = &inner;
0110     trHits.push_back(outerHit);
0111     trHits.push_back(middleHit);
0112     //firstHit  = outerHit;
0113     //secondHit = middleHit;
0114   } else {
0115     LogDebug("SeedFromGenericPairOrTriplet") << "Seed from outsideIn oppositeToMomentum OR insideOut alongMomentum";
0116     firstPoint = &inner;
0117     secondPoint = &middle;
0118     thirdPoint = &outer;
0119     trHits.push_back(innerHit);
0120     trHits.push_back(middleHit);
0121     //firstHit  = innerHit;
0122     //secondHit = middleHit;
0123   }
0124   if (dir == oppositeToMomentum)
0125     momentumSign = -1;
0126   FastHelix helix(*thirdPoint, *secondPoint, *firstPoint, theMagfield->nominalValue(), theMagfield);
0127   FreeTrajectoryState originalStartingState = helix.stateAtVertex();
0128   LogDebug("SeedFromGenericPairOrTriplet") << "originalStartingState " << originalStartingState;
0129   /*GlobalTrajectoryParameters originalPar = originalStartingState.parameters();
0130         GlobalTrajectoryParameters newPar = GlobalTrajectoryParameters(originalPar.position(), 
0131                                                                        momentumSign*originalPar.momentum(),
0132                                                                        originalPar.charge(),
0133                                                                        &originalPar.magneticField());
0134     */
0135 
0136   /*FastCircle helix(*thirdPoint, *secondPoint, *firstPoint);
0137     GlobalTrajectoryParameters newPar = GlobalTrajectoryParameters(*secondPoint,
0138                                                                        momentumSign*originalPar.momentum(),
0139                                                                        originalPar.charge(),
0140                                                                        &originalPar.magneticField());*/
0141   //FreeTrajectoryState* startingState = new FreeTrajectoryState(newPar, initialError(trHits[0]));//trHits[1]));
0142   if (!qualityFilter(momentumSign * originalStartingState.momentum()))
0143     return nullptr;
0144   //TrajectorySeed* seed = buildSeed(startingState, trHits, dir);
0145   TrajectorySeed* seed =
0146       buildSeed(momentumSign * originalStartingState.momentum(), originalStartingState.charge(), trHits, dir);
0147   return seed;
0148 }
0149 
0150 TrajectorySeed* SeedFromGenericPairOrTriplet::seedFromPair(const SeedingHitSet& hits,
0151                                                            const PropagationDirection& dir,
0152                                                            const NavigationDirection& seedDir,
0153                                                            int charge) const {
0154   if (hits.size() != 2) {
0155     throw cms::Exception("CombinatorialSeedGeneratorForCosmics")
0156         << "call to SeedFromGenericPairOrTriplet::seedFromPair with " << hits.size() << " hits ";
0157   }
0158   auto innerHit = hits[0];
0159   auto outerHit = hits[1];
0160   GlobalPoint inner = theTracker->idToDet(innerHit->geographicalId())->surface().toGlobal(innerHit->localPosition());
0161   GlobalPoint outer = theTracker->idToDet(outerHit->geographicalId())->surface().toGlobal(outerHit->localPosition());
0162   LogDebug("SeedFromGenericPairOrTriplet")
0163       << "Using the following hits: outer(r, phi, theta) (" << outer.perp() << ", " << outer.phi() << ","
0164       << outer.theta() << ")    inner (" << inner.perp() << ", " << inner.phi() << "," << inner.theta() << ")";
0165   GlobalPoint* firstPoint = nullptr;
0166   GlobalPoint* secondPoint = nullptr;
0167   int momentumSign = 1;
0168   //const TrackingRecHit* firstHit  = 0;
0169   //const TrackingRecHit* secondHit = 0;
0170   std::vector<const BaseTrackerRecHit*> trHits;
0171   //choose the prop dir and hit order accordingly to where the seed is made
0172   if (seedDir == outsideIn) {
0173     LogDebug("SeedFromGenericPairOrTriplet") << "Seed from outsideIn alongMomentum OR insideOut oppositeToMomentum";
0174     firstPoint = &outer;
0175     secondPoint = &inner;
0176     //firstHit  = outerHit;
0177     //secondHit = innerHit;
0178     trHits.push_back(outerHit);
0179     trHits.push_back(innerHit);
0180   } else {
0181     LogDebug("SeedFromGenericPairOrTriplet") << "Seed from outsideIn oppositeToMomentum OR insideOut alongMomentum";
0182     firstPoint = &inner;
0183     secondPoint = &outer;
0184     momentumSign = -1;
0185     //firstHit  = innerHit;
0186     //secondHit = outerHit;
0187     trHits.push_back(innerHit);
0188     trHits.push_back(outerHit);
0189   }
0190   if (dir == oppositeToMomentum)
0191     momentumSign = -1;
0192   GlobalVector momentum = momentumSign * theP * (*secondPoint - *firstPoint).unit();
0193   /*GlobalTrajectoryParameters gtp(*firstPoint, 
0194                                        momentum,
0195                                        charge,
0196                                        &(*theMagfield));*/
0197   //FreeTrajectoryState* startingState = new FreeTrajectoryState(gtp,initialError(trHits[0]));//trHits[1]));
0198   //if (!qualityFilter(startingState, hits)) return 0;
0199   //TrajectorySeed* seed = buildSeed(startingState, firstHit, dir);
0200   //TrajectorySeed* seed = buildSeed(startingState, trHits, dir);
0201   TrajectorySeed* seed = buildSeed(momentum, charge, trHits, dir);
0202   return seed;
0203 }
0204 
0205 TrajectorySeed* SeedFromGenericPairOrTriplet::buildSeed(const GlobalVector& momentum,
0206                                                         int charge,
0207                                                         //const TrackingRecHit* firsthit,
0208                                                         std::vector<const BaseTrackerRecHit*>& trHits,
0209                                                         const PropagationDirection& dir) const {
0210   const Propagator* propagator = thePropagatorAlong;
0211   if (dir == oppositeToMomentum)
0212     propagator = thePropagatorOpposite;
0213 
0214   //debug
0215   GlobalPoint first = theTracker->idToDet(trHits[0]->geographicalId())->surface().toGlobal(trHits[0]->localPosition());
0216   GlobalPoint second = theTracker->idToDet(trHits[1]->geographicalId())->surface().toGlobal(trHits[1]->localPosition());
0217   LogDebug("SeedFromGenericPairOrTriplet")
0218       << "Using the following hits: first(r, phi, theta) (" << first.perp() << ", " << first.phi() << ","
0219       << first.theta() << ")    second (" << second.perp() << ", " << second.phi() << "," << second.theta() << ")";
0220   //build an initial trajectory state with big errors,
0221   //then update it with the first hit
0222   auto transHit = trHits[0];
0223   LocalVector lmom = transHit->surface()->toLocal(momentum);
0224   LocalTrajectoryParameters ltp(charge / momentum.mag(),
0225                                 lmom.x() / lmom.z(),
0226                                 lmom.y() / lmom.z(),
0227                                 transHit->localPosition().x(),
0228                                 transHit->localPosition().y(),
0229                                 lmom.z() > 0. ? 1. : -1.);
0230 
0231   LocalTrajectoryError lte(100.,
0232                            100.,
0233                            (1. + (lmom.x() / lmom.z()) * (lmom.x() / lmom.z())),
0234                            (1. + (lmom.y() / lmom.z()) * (lmom.y() / lmom.z())),
0235                            1. / momentum.mag());
0236 
0237   TrajectoryStateOnSurface initialState(ltp, lte, *transHit->surface(), &(*theMagfield));
0238   KFUpdator updator;
0239   TrajectoryStateOnSurface startingState = updator.update(initialState, *transHit);
0240 
0241   auto transHit2 = trHits[1];
0242 
0243   //TrajectoryStateOnSurface seedTSOS(*startingState, *plane);
0244   const TrajectoryStateOnSurface propTSOS = propagator->propagate(startingState, *(transHit2->surface()));
0245   if (!propTSOS.isValid()) {
0246     LogDebug("SeedFromGenericPairOrTriplet") << "first propagation failed";
0247     return nullptr;
0248   }
0249   TrajectoryStateOnSurface seedTSOS = updator.update(propTSOS, *transHit2);
0250   if (!seedTSOS.isValid()) {
0251     LogDebug("SeedFromGenericPairOrTriplet") << "first update failed";
0252     return nullptr;
0253   }
0254   LogDebug("SeedFromGenericPairOrTriplet") << "starting TSOS " << seedTSOS;
0255   seedTSOS.rescaleError(theErrorRescaling);
0256   PTrajectoryStateOnDet PTraj =
0257       trajectoryStateTransform::persistentState(seedTSOS, trHits[1]->geographicalId().rawId());
0258   edm::OwnVector<TrackingRecHit> seed_hits;
0259   //build the transientTrackingRecHit for the starting hit, then call the method clone to rematch if needed
0260   for (auto ihits = trHits.begin(); ihits != trHits.end(); ihits++) {
0261     seed_hits.push_back((*ihits)->clone());
0262   }
0263   TrajectorySeed* seed = new TrajectorySeed(PTraj, seed_hits, dir);
0264   return seed;
0265 }
0266 /*
0267 CurvilinearTrajectoryError SeedFromGenericPairOrTriplet::initialError(const TrackingRecHit* rechit) {
0268         SeedingHitSet::RecHitPointer transHit =  theBuilder->build(rechit);
0269         //AlgebraicSymMatrix C(5,1);
0270     AlgebraicSymMatrix55 C = AlgebraicMatrixID();
0271         //C*=0.1;
0272         if (theSetMomentum){
0273         C(0,0)=1/theP;
0274                 C(3,3)=transHit->globalPositionError().cxx();
0275                 C(4,4)=transHit->globalPositionError().czz();
0276         
0277                 //C[0][0]=1/theP;
0278                 //C[3][3]=transHit->globalPositionError().cxx();
0279                 //C[4][4]=transHit->globalPositionError().czz();
0280         } else {
0281                 float zErr = transHit->globalPositionError().czz();
0282                 float transverseErr = transHit->globalPositionError().cxx(); // assume equal cxx cyy
0283         C(3,3) = transverseErr;
0284                 C(4,4) = zErr;
0285                 //C[3][3] = transverseErr;
0286                 //C[4][4] = zErr;
0287         }
0288 
0289         return CurvilinearTrajectoryError(C);
0290 }
0291 */
0292 bool SeedFromGenericPairOrTriplet::qualityFilter(const SeedingHitSet& hits) const {
0293   //if we are setting the momentum with the PSet we look for aligned hits
0294   if (theSetMomentum) {
0295     if (hits.size() == 3) {
0296       std::vector<GlobalPoint> gPoints;
0297       unsigned int nHits = hits.size();
0298       for (unsigned int iHit = 0; iHit < nHits; ++iHit)
0299         gPoints.push_back(hits[iHit]->globalPosition());
0300       unsigned int subid = (*hits[0]).geographicalId().subdetId();
0301       if (subid == StripSubdetector::TEC || subid == StripSubdetector::TID) {
0302         LogDebug("SeedFromGenericPairOrTriplet")
0303             << "In the endcaps we cannot decide if hits are aligned with only phi and z";
0304         return true;
0305       }
0306       FastCircle circle(gPoints[0], gPoints[1], gPoints[2]);
0307       if (circle.rho() < 500 && circle.rho() != 0) {
0308         LogDebug("SeedFromGenericPairOrTriplet") << "Seed qualityFilter rejected because rho = " << circle.rho();
0309         return false;
0310       }
0311     }
0312   }
0313   return true;
0314 }
0315 
0316 bool SeedFromGenericPairOrTriplet::qualityFilter(const GlobalVector& momentum) const {
0317   if (!theSetMomentum) {
0318     if (momentum.perp() < theP)
0319       return false;
0320   }
0321   return true;
0322 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!