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
0102
0103
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
0113
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
0122
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
0130
0131
0132
0133
0134
0135
0136
0137
0138
0139
0140
0141
0142 if (!qualityFilter(momentumSign * originalStartingState.momentum()))
0143 return nullptr;
0144
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
0169
0170 std::vector<const BaseTrackerRecHit*> trHits;
0171
0172 if (seedDir == outsideIn) {
0173 LogDebug("SeedFromGenericPairOrTriplet") << "Seed from outsideIn alongMomentum OR insideOut oppositeToMomentum";
0174 firstPoint = &outer;
0175 secondPoint = &inner;
0176
0177
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
0186
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
0194
0195
0196
0197
0198
0199
0200
0201 TrajectorySeed* seed = buildSeed(momentum, charge, trHits, dir);
0202 return seed;
0203 }
0204
0205 TrajectorySeed* SeedFromGenericPairOrTriplet::buildSeed(const GlobalVector& momentum,
0206 int charge,
0207
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
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
0221
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
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
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
0268
0269
0270
0271
0272
0273
0274
0275
0276
0277
0278
0279
0280
0281
0282
0283
0284
0285
0286
0287
0288
0289
0290
0291
0292 bool SeedFromGenericPairOrTriplet::qualityFilter(const SeedingHitSet& hits) const {
0293
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 }