File indexing completed on 2024-04-06 12:28:45
0001 #include "RecoTracker/SpecialSeedGenerators/interface/SeedGeneratorForCosmics.h"
0002 #include "RecoTracker/TkHitPairs/interface/CosmicLayerPairs.h"
0003 #include "RecoTracker/PixelSeeding/interface/CosmicLayerTriplets.h"
0004 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0005 #include "FWCore/Utilities/interface/isFinite.h"
0006 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0007 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0008 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0009 void SeedGeneratorForCosmics::init(const SiStripRecHit2DCollection& collstereo,
0010 const SiStripRecHit2DCollection& collrphi,
0011 const SiStripMatchedRecHit2DCollection& collmatched,
0012 const edm::EventSetup& iSetup) {
0013 magfield = iSetup.getHandle(theMagfieldToken);
0014 tracker = iSetup.getHandle(theTrackerToken);
0015 thePropagatorAl = new PropagatorWithMaterial(alongMomentum, 0.1057, &(*magfield));
0016 thePropagatorOp = new PropagatorWithMaterial(oppositeToMomentum, 0.1057, &(*magfield));
0017 theUpdator = new KFUpdator();
0018
0019 LogDebug("CosmicSeedFinder") << " Hits built with " << hitsforseeds << " hits";
0020
0021 GeometricSearchTracker const& track = iSetup.getData(theSearchTrackerToken);
0022 TrackerTopology const& ttopo = iSetup.getData(theTTopoToken);
0023
0024 CosmicLayerPairs cosmiclayers(geometry, collrphi, collmatched, track, ttopo);
0025
0026 thePairGenerator = new CosmicHitPairGenerator(cosmiclayers, *tracker);
0027 HitPairs.clear();
0028 if ((hitsforseeds == "pairs") || (hitsforseeds == "pairsandtriplets")) {
0029 thePairGenerator->hitPairs(region, HitPairs);
0030 }
0031
0032 CosmicLayerTriplets cosmiclayers2(geometry, collrphi, track, ttopo);
0033 theTripletGenerator = new CosmicHitTripletGenerator(cosmiclayers2, *tracker);
0034 HitTriplets.clear();
0035 if ((hitsforseeds == "triplets") || (hitsforseeds == "pairsandtriplets")) {
0036 theTripletGenerator->hitTriplets(region, HitTriplets);
0037 }
0038 }
0039
0040 SeedGeneratorForCosmics::SeedGeneratorForCosmics(edm::ParameterSet const& conf, edm::ConsumesCollector iCC)
0041 : maxSeeds_(conf.getParameter<int32_t>("maxSeeds")),
0042 theMagfieldToken(iCC.esConsumes()),
0043 theTrackerToken(iCC.esConsumes()),
0044 theSearchTrackerToken(iCC.esConsumes()),
0045 theTTopoToken(iCC.esConsumes()) {
0046 float ptmin = conf.getParameter<double>("ptMin");
0047 float originradius = conf.getParameter<double>("originRadius");
0048 float halflength = conf.getParameter<double>("originHalfLength");
0049 float originz = conf.getParameter<double>("originZPosition");
0050 seedpt = conf.getParameter<double>("SeedPt");
0051
0052 geometry = conf.getUntrackedParameter<std::string>("GeometricStructure", "STANDARD");
0053 region = GlobalTrackingRegion(ptmin, originradius, halflength, originz);
0054 hitsforseeds = conf.getUntrackedParameter<std::string>("HitsForSeeds", "pairs");
0055 edm::LogInfo("SeedGeneratorForCosmics")
0056 << " PtMin of track is " << ptmin << " The Radius of the cylinder for seeds is " << originradius << "cm"
0057 << " The set Seed Momentum" << seedpt;
0058
0059
0060 positiveYOnly = conf.getParameter<bool>("PositiveYOnly");
0061 negativeYOnly = conf.getParameter<bool>("NegativeYOnly");
0062
0063 }
0064
0065 void SeedGeneratorForCosmics::run(const SiStripRecHit2DCollection& collstereo,
0066 const SiStripRecHit2DCollection& collrphi,
0067 const SiStripMatchedRecHit2DCollection& collmatched,
0068 const edm::EventSetup& c,
0069 TrajectorySeedCollection& output) {
0070 init(collstereo, collrphi, collmatched, c);
0071 seeds(output, region);
0072 delete thePairGenerator;
0073 delete theTripletGenerator;
0074 delete thePropagatorAl;
0075 delete thePropagatorOp;
0076 delete theUpdator;
0077 }
0078 bool SeedGeneratorForCosmics::seeds(TrajectorySeedCollection& output, const TrackingRegion& region) {
0079 LogDebug("CosmicSeedFinder") << "Number of triplets " << HitTriplets.size();
0080 LogDebug("CosmicSeedFinder") << "Number of pairs " << HitPairs.size();
0081
0082 for (unsigned int it = 0; it < HitTriplets.size(); it++) {
0083
0084
0085
0086
0087
0088
0089
0090
0091
0092 GlobalPoint inner = tracker->idToDet((*(HitTriplets[it].inner())).geographicalId())
0093 ->surface()
0094 .toGlobal((*(HitTriplets[it].inner())).localPosition());
0095
0096 GlobalPoint middle = tracker->idToDet((*(HitTriplets[it].middle())).geographicalId())
0097 ->surface()
0098 .toGlobal((*(HitTriplets[it].middle())).localPosition());
0099
0100 GlobalPoint outer = tracker->idToDet((*(HitTriplets[it].outer())).geographicalId())
0101 ->surface()
0102 .toGlobal((*(HitTriplets[it].outer())).localPosition());
0103
0104 SeedingHitSet::ConstRecHitPointer outrhit = HitTriplets[it].outer();
0105
0106 SeedingHitSet::ConstRecHitPointer innrhit = HitTriplets[it].inner();
0107 if (positiveYOnly && (outrhit->globalPosition().y() < 0 || innrhit->globalPosition().y() < 0 ||
0108 outrhit->globalPosition().y() < innrhit->globalPosition().y()))
0109 continue;
0110 if (negativeYOnly && (outrhit->globalPosition().y() > 0 || innrhit->globalPosition().y() > 0 ||
0111 outrhit->globalPosition().y() > innrhit->globalPosition().y()))
0112 continue;
0113
0114
0115 edm::OwnVector<TrackingRecHit> hits;
0116 hits.push_back(HitTriplets[it].outer()->hit()->clone());
0117 FastHelix helix(inner, middle, outer, magfield->nominalValue(), &(*magfield));
0118 GlobalVector gv = helix.stateAtVertex().momentum();
0119 float ch = helix.stateAtVertex().charge();
0120 float mom2 = gv.x() * gv.x() + gv.y() * gv.y() + gv.z() * gv.z();
0121 if (mom2 > 1e06f || edm::isNotFinite(mom2))
0122 continue;
0123
0124 if (gv.y() > 0) {
0125 gv = -1. * gv;
0126 ch = -1. * ch;
0127 }
0128
0129 GlobalTrajectoryParameters Gtp(outer, gv, int(ch), &(*magfield));
0130 FreeTrajectoryState CosmicSeed(Gtp, CurvilinearTrajectoryError(AlgebraicSymMatrix55(AlgebraicMatrixID())));
0131 if ((outer.y() - inner.y()) > 0) {
0132 const TSOS outerState = thePropagatorAl->propagate(
0133 CosmicSeed, tracker->idToDet((*(HitTriplets[it].outer())).geographicalId())->surface());
0134 if (outerState.isValid()) {
0135 LogDebug("CosmicSeedFinder") << "outerState " << outerState;
0136 const TSOS outerUpdated = theUpdator->update(outerState, *outrhit);
0137 if (outerUpdated.isValid()) {
0138 LogDebug("CosmicSeedFinder") << "outerUpdated " << outerUpdated;
0139
0140 output.push_back(TrajectorySeed(trajectoryStateTransform::persistentState(
0141 outerUpdated, (*(HitTriplets[it].outer())).geographicalId().rawId()),
0142 hits,
0143 alongMomentum));
0144
0145 if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
0146 edm::LogError("TooManySeeds") << "Found too many seeds, bailing out.\n";
0147 output.clear();
0148 return false;
0149 }
0150 }
0151 }
0152 } else {
0153 const TSOS outerState = thePropagatorOp->propagate(
0154 CosmicSeed, tracker->idToDet((*(HitTriplets[it].outer())).geographicalId())->surface());
0155 if (outerState.isValid()) {
0156 LogDebug("CosmicSeedFinder") << "outerState " << outerState;
0157 const TSOS outerUpdated = theUpdator->update(outerState, *outrhit);
0158 if (outerUpdated.isValid()) {
0159 LogDebug("CosmicSeedFinder") << "outerUpdated " << outerUpdated;
0160
0161 output.push_back(TrajectorySeed(trajectoryStateTransform::persistentState(
0162 outerUpdated, (*(HitTriplets[it].outer())).geographicalId().rawId()),
0163 hits,
0164 oppositeToMomentum));
0165
0166 if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
0167 edm::LogError("TooManySeeds") << "Found too many seeds, bailing out.\n";
0168 output.clear();
0169 return false;
0170 }
0171 }
0172 }
0173 }
0174 }
0175
0176 for (unsigned int is = 0; is < HitPairs.size(); is++) {
0177 GlobalPoint inner = tracker->idToDet((*(HitPairs[is].inner())).geographicalId())
0178 ->surface()
0179 .toGlobal((*(HitPairs[is].inner())).localPosition());
0180 GlobalPoint outer = tracker->idToDet((*(HitPairs[is].outer())).geographicalId())
0181 ->surface()
0182 .toGlobal((*(HitPairs[is].outer())).localPosition());
0183
0184 LogDebug("CosmicSeedFinder") << "inner point of the seed " << inner << " outer point of the seed " << outer;
0185 SeedingHitSet::ConstRecHitPointer outrhit = HitPairs[is].outer();
0186
0187 SeedingHitSet::ConstRecHitPointer innrhit = HitPairs[is].inner();
0188 if (positiveYOnly && (outrhit->globalPosition().y() < 0 || innrhit->globalPosition().y() < 0 ||
0189 outrhit->globalPosition().y() < innrhit->globalPosition().y()))
0190 continue;
0191 if (negativeYOnly && (outrhit->globalPosition().y() > 0 || innrhit->globalPosition().y() > 0 ||
0192 outrhit->globalPosition().y() > innrhit->globalPosition().y()))
0193 continue;
0194
0195
0196 edm::OwnVector<TrackingRecHit> hits;
0197 hits.push_back(HitPairs[is].outer()->hit()->clone());
0198
0199
0200 for (int i = 0; i < 2; i++) {
0201
0202 int predsign = (2 * i) - 1;
0203 if ((outer.y() - inner.y()) > 0) {
0204 GlobalTrajectoryParameters Gtp(
0205 outer, (inner - outer) * (seedpt / (inner - outer).mag()), predsign, &(*magfield));
0206
0207 FreeTrajectoryState CosmicSeed(Gtp, CurvilinearTrajectoryError(AlgebraicSymMatrix55(AlgebraicMatrixID())));
0208
0209 LogDebug("CosmicSeedFinder") << " FirstTSOS " << CosmicSeed;
0210
0211 const TSOS outerState = thePropagatorAl->propagate(
0212 CosmicSeed, tracker->idToDet((*(HitPairs[is].outer())).geographicalId())->surface());
0213 if (outerState.isValid()) {
0214 LogDebug("CosmicSeedFinder") << "outerState " << outerState;
0215 const TSOS outerUpdated = theUpdator->update(outerState, *outrhit);
0216 if (outerUpdated.isValid()) {
0217 LogDebug("CosmicSeedFinder") << "outerUpdated " << outerUpdated;
0218
0219 PTrajectoryStateOnDet const& PTraj = trajectoryStateTransform::persistentState(
0220 outerUpdated, (*(HitPairs[is].outer())).geographicalId().rawId());
0221
0222 output.push_back(TrajectorySeed(PTraj, hits, alongMomentum));
0223
0224 if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
0225 edm::LogError("TooManySeeds") << "Found too many seeds, bailing out.\n";
0226 output.clear();
0227 return false;
0228 }
0229
0230 } else
0231 edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first update failed ";
0232 } else
0233 edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first propagation failed ";
0234
0235 } else {
0236 GlobalTrajectoryParameters Gtp(
0237 outer, (outer - inner) * (seedpt / (outer - inner).mag()), predsign, &(*magfield));
0238 FreeTrajectoryState CosmicSeed(Gtp, CurvilinearTrajectoryError(AlgebraicSymMatrix55(AlgebraicMatrixID())));
0239 LogDebug("CosmicSeedFinder") << " FirstTSOS " << CosmicSeed;
0240
0241 const TSOS outerState = thePropagatorAl->propagate(
0242 CosmicSeed, tracker->idToDet((*(HitPairs[is].outer())).geographicalId())->surface());
0243 if (outerState.isValid()) {
0244 LogDebug("CosmicSeedFinder") << "outerState " << outerState;
0245 const TSOS outerUpdated = theUpdator->update(outerState, *outrhit);
0246 if (outerUpdated.isValid()) {
0247 LogDebug("CosmicSeedFinder") << "outerUpdated " << outerUpdated;
0248
0249 PTrajectoryStateOnDet const& PTraj = trajectoryStateTransform::persistentState(
0250 outerUpdated, (*(HitPairs[is].outer())).geographicalId().rawId());
0251
0252 output.push_back(TrajectorySeed(PTraj, hits, oppositeToMomentum));
0253
0254 if ((maxSeeds_ > 0) && (output.size() > size_t(maxSeeds_))) {
0255 edm::LogError("TooManySeeds") << "Found too many seeds, bailing out.\n";
0256 output.clear();
0257 return false;
0258 }
0259
0260 } else
0261 edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first update failed ";
0262 } else
0263 edm::LogWarning("CosmicSeedFinder") << " SeedForCosmics first propagation failed ";
0264 }
0265 }
0266 }
0267 return true;
0268 }