File indexing completed on 2021-05-18 09:15:26
0001 #include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
0002 #include "DataFormats/GeometrySurface/interface/Surface.h"
0003 #include "DataFormats/GeometrySurface/interface/TangentPlane.h"
0004 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0005 #include "DataFormats/ParticleFlowReco/interface/ConvBremSeed.h"
0006 #include "DataFormats/ParticleFlowReco/interface/ConvBremSeedFwd.h"
0007 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrack.h"
0008 #include "DataFormats/ParticleFlowReco/interface/GsfPFRecTrackFwd.h"
0009 #include "DataFormats/ParticleFlowReco/interface/PFBrem.h"
0010 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0011 #include "DataFormats/ParticleFlowReco/interface/PFClusterFwd.h"
0012 #include "DataFormats/TrackReco/interface/Track.h"
0013 #include "DataFormats/TrackerRecHit2D/interface/SiPixelRecHitCollection.h"
0014 #include "DataFormats/TrackerRecHit2D/interface/SiStripMatchedRecHit2DCollection.h"
0015 #include "DataFormats/TrackerRecHit2D/interface/SiStripRecHit2DCollection.h"
0016 #include "DataFormats/TrackingRecHit/interface/TrackingRecHit.h"
0017 #include "DataFormats/TrackingRecHit/interface/TrackingRecHitFwd.h"
0018 #include "DataFormats/TrajectorySeed/interface/TrajectorySeed.h"
0019 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0020 #include "FWCore/Framework/interface/stream/EDProducer.h"
0021 #include "FWCore/Framework/interface/ESHandle.h"
0022 #include "FWCore/Framework/interface/Event.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0025 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0026 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0027 #include "FastSimulation/ParticlePropagator/interface/MagneticFieldMapRecord.h"
0028 #include "FastSimulation/ParticlePropagator/interface/ParticlePropagator.h"
0029 #include "FastSimulation/TrackerSetup/interface/TrackerInteractionGeometry.h"
0030 #include "FastSimulation/TrackerSetup/interface/TrackerInteractionGeometryRecord.h"
0031 #include "FastSimulation/TrajectoryManager/interface/InsideBoundsMeasurementEstimator.h"
0032 #include "FastSimulation/TrajectoryManager/interface/LocalMagneticField.h"
0033 #include "Geometry/Records/interface/TrackerDigiGeometryRecord.h"
0034 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0035 #include "MagneticField/Engine/interface/MagneticField.h"
0036 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0037 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0038 #include "RecoTracker/TkDetLayers/interface/GeometricSearchTracker.h"
0039 #include "RecoTracker/TkSeedGenerator/interface/FastHelix.h"
0040 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0041 #include "TrackingTools/DetLayers/interface/BarrelDetLayer.h"
0042 #include "TrackingTools/DetLayers/interface/DetLayer.h"
0043 #include "TrackingTools/DetLayers/interface/ForwardDetLayer.h"
0044 #include "TrackingTools/GeomPropagators/interface/AnalyticalPropagator.h"
0045 #include "TrackingTools/KalmanUpdators/interface/KFUpdator.h"
0046 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0047 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0048 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0049
0050 #include "TMath.h"
0051
0052 #include <memory>
0053
0054 class ConvBremSeedProducer : public edm::stream::EDProducer<> {
0055 public:
0056 explicit ConvBremSeedProducer(const edm::ParameterSet&);
0057
0058 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0059
0060 private:
0061 void beginRun(const edm::Run&, const edm::EventSetup&) override;
0062 void produce(edm::Event&, const edm::EventSetup&) override;
0063 void endRun(const edm::Run&, const edm::EventSetup&) override;
0064 void initializeLayerMap();
0065 std::vector<const DetLayer*> theLayerMap;
0066 TrajectoryStateOnSurface makeTrajectoryState(const DetLayer* layer,
0067 const ParticlePropagator& pp,
0068 const MagneticField* field) const;
0069 const DetLayer* detLayer(const TrackerLayer& layer, float zpos) const;
0070
0071 bool isGsfTrack(const reco::Track&, const TrackingRecHit*);
0072
0073 int GoodCluster(const BaseParticlePropagator& bpg,
0074 const reco::PFClusterCollection& pfc,
0075 float minep,
0076 bool sec = false);
0077
0078 std::vector<bool> sharedHits(const std::vector<std::pair<TrajectorySeed, std::pair<GlobalVector, float> > >&);
0079
0080 edm::ParameterSet conf_;
0081 const GeometricSearchTracker* geomSearchTracker_;
0082 const TrackerInteractionGeometry* geometry_;
0083 const TrackerGeometry* tracker_;
0084 const MagneticField* magfield_;
0085 const MagneticFieldMap* fieldMap_;
0086 const PropagatorWithMaterial* propagator_;
0087 const KFUpdator* kfUpdator_;
0088 const TransientTrackingRecHitBuilder* hitBuilder_;
0089 std::vector<const DetLayer*> layerMap_;
0090 int negLayerOffset_;
0091
0092 math::XYZVector B_;
0093
0094
0095 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_;
0096 const edm::ESGetToken<GeometricSearchTracker, TrackerRecoGeometryRecord> geomSearchTrackerToken_;
0097 const edm::ESGetToken<TrackerInteractionGeometry, TrackerInteractionGeometryRecord> geometryToken_;
0098 const edm::ESGetToken<TrackerGeometry, TrackerDigiGeometryRecord> trackerToken_;
0099 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magFieldToken_beginRun_;
0100 const edm::ESGetToken<MagneticFieldMap, MagneticFieldMapRecord> magFieldMapToken_;
0101 const edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> hitBuilderToken_;
0102 };
0103
0104 #include "FWCore/Framework/interface/MakerMacros.h"
0105 DEFINE_FWK_MODULE(ConvBremSeedProducer);
0106
0107 void ConvBremSeedProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0108
0109 edm::ParameterSetDescription desc;
0110 desc.add<edm::InputTag>("pixelRecHits", edm::InputTag("gsPixelRecHits"));
0111 desc.add<edm::InputTag>("matchedrecHits", edm::InputTag("gsStripRecHits", "matchedRecHit"));
0112 desc.add<std::string>("TTRHBuilder", "WithTrackAngle");
0113 desc.add<edm::InputTag>("rphirecHits", edm::InputTag("gsStripRecHits", "rphiRecHit"));
0114 desc.add<edm::InputTag>("PFClusters", edm::InputTag("particleFlowClusterECAL"));
0115 desc.add<edm::InputTag>("PFRecTrackLabel", edm::InputTag("pfTrackElec"));
0116 }
0117
0118 using namespace edm;
0119 using namespace std;
0120 using namespace reco;
0121
0122 ConvBremSeedProducer::ConvBremSeedProducer(const ParameterSet& iConfig)
0123 : conf_(iConfig),
0124 fieldMap_(nullptr),
0125 layerMap_(56, static_cast<const DetLayer*>(nullptr)),
0126 negLayerOffset_(27),
0127 magFieldToken_(esConsumes()),
0128 geomSearchTrackerToken_(esConsumes<edm::Transition::BeginRun>()),
0129 geometryToken_(esConsumes<edm::Transition::BeginRun>()),
0130 trackerToken_(esConsumes<edm::Transition::BeginRun>()),
0131 magFieldToken_beginRun_(esConsumes<edm::Transition::BeginRun>()),
0132 magFieldMapToken_(esConsumes<edm::Transition::BeginRun>()),
0133 hitBuilderToken_(
0134 esConsumes<edm::Transition::BeginRun>(edm::ESInputTag("", conf_.getParameter<string>("TTRHBuilder")))) {
0135 produces<ConvBremSeedCollection>();
0136 }
0137
0138 void ConvBremSeedProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0139 LogDebug("ConvBremSeedProducerProducer") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run();
0140
0141 constexpr float pfmass = 0.0005;
0142
0143
0144
0145
0146 Handle<PFClusterCollection> PfC;
0147 iEvent.getByLabel(conf_.getParameter<InputTag>("PFClusters"), PfC);
0148 const PFClusterCollection& PPP = *(PfC.product());
0149
0150
0151 Handle<SiPixelRecHitCollection> pixelHits;
0152 iEvent.getByLabel(conf_.getParameter<InputTag>("pixelRecHits"), pixelHits);
0153
0154
0155 Handle<SiStripRecHit2DCollection> rphirecHits;
0156 iEvent.getByLabel(conf_.getParameter<InputTag>("rphirecHits"), rphirecHits);
0157 Handle<SiStripMatchedRecHit2DCollection> matchedrecHits;
0158 iEvent.getByLabel(conf_.getParameter<InputTag>("matchedrecHits"), matchedrecHits);
0159
0160
0161 Handle<GsfPFRecTrackCollection> thePfRecTrackCollection;
0162 iEvent.getByLabel(conf_.getParameter<InputTag>("PFRecTrackLabel"), thePfRecTrackCollection);
0163 const GsfPFRecTrackCollection PfRTkColl = *(thePfRecTrackCollection.product());
0164
0165
0166 auto output = std::make_unique<ConvBremSeedCollection>();
0167
0168
0169 vector<pair<TrajectorySeed, pair<GlobalVector, float> > > unclean;
0170
0171 vector<vector<long int> > tripl;
0172
0173 initializeLayerMap();
0174
0175
0176
0177 for (unsigned int ipft = 0; ipft < PfRTkColl.size(); ipft++) {
0178 GsfPFRecTrackRef pft(thePfRecTrackCollection, ipft);
0179 LogDebug("ConvBremSeedProducerProducer") << "NEW GsfPFRecTRACK ";
0180 float eta_br = 0;
0181 unclean.clear();
0182 tripl.clear();
0183 vector<int> gc;
0184 auto const& gsfRecHits = *pft->gsfTrackRef();
0185 float pfoutenergy = sqrt((pfmass * pfmass) + pft->gsfTrackRef()->outerMomentum().Mag2());
0186 XYZTLorentzVector mom = XYZTLorentzVector(pft->gsfTrackRef()->outerMomentum().x(),
0187 pft->gsfTrackRef()->outerMomentum().y(),
0188 pft->gsfTrackRef()->outerMomentum().z(),
0189 pfoutenergy);
0190 XYZTLorentzVector pos = XYZTLorentzVector(pft->gsfTrackRef()->outerPosition().x(),
0191 pft->gsfTrackRef()->outerPosition().y(),
0192 pft->gsfTrackRef()->outerPosition().z(),
0193 0.);
0194 BaseParticlePropagator theOutParticle =
0195 BaseParticlePropagator(RawParticle(mom, pos, pft->gsfTrackRef()->charge()), 0, 0, B_.z());
0196
0197
0198 gc.push_back(GoodCluster(theOutParticle, PPP, 0.5));
0199
0200 vector<PFBrem> brem = (*pft).PFRecBrem();
0201 vector<PFBrem>::iterator ib = brem.begin();
0202 vector<PFBrem>::iterator ib_end = brem.end();
0203 LogDebug("ConvBremSeedProducerProducer") << "NUMBER OF BREMS " << brem.size();
0204
0205
0206 for (; ib != ib_end; ++ib) {
0207 XYZTLorentzVector mom = pft->trajectoryPoint(ib->indTrajPoint()).momentum();
0208 XYZTLorentzVector pos = XYZTLorentzVector(pft->trajectoryPoint(ib->indTrajPoint()).position().x(),
0209 pft->trajectoryPoint(ib->indTrajPoint()).position().y(),
0210 pft->trajectoryPoint(ib->indTrajPoint()).position().z(),
0211 0);
0212
0213
0214 if (pos.Rho() > 80)
0215 continue;
0216 if ((pos.Rho() > 5) && (fabs(ib->SigmaDeltaP() / ib->DeltaP()) > 3))
0217 continue;
0218 if (fabs(ib->DeltaP()) < 3)
0219 continue;
0220 eta_br = mom.eta();
0221 vector<vector<long int> > Idd;
0222
0223 BaseParticlePropagator p(RawParticle(mom, pos, 0.), 0, 0, B_.z());
0224 gc.push_back(GoodCluster(p, PPP, 0.2));
0225
0226 ParticlePropagator PP(p, fieldMap_, nullptr);
0227
0228
0229 list<TrackerLayer>::const_iterator cyliter = geometry_->cylinderBegin();
0230 for (; cyliter != geometry_->cylinderEnd(); ++cyliter) {
0231
0232 if (!(cyliter->sensitive()))
0233 continue;
0234 PP.setPropagationConditions(*cyliter);
0235 PP.propagate();
0236 if (PP.getSuccess() == 0)
0237 continue;
0238
0239
0240 LocalMagneticField mf(PP.getMagneticField());
0241 AnalyticalPropagator alongProp(&mf, anyDirection);
0242 InsideBoundsMeasurementEstimator est;
0243 const DetLayer* tkLayer = detLayer(*cyliter, PP.particle().Z());
0244 if (&(*tkLayer) == nullptr)
0245 continue;
0246 TrajectoryStateOnSurface trajState = makeTrajectoryState(tkLayer, PP, &mf);
0247
0248 auto compat = tkLayer->compatibleDets(trajState, alongProp, est);
0249 vector<long int> temp;
0250 if (compat.empty())
0251 continue;
0252
0253 for (auto i = compat.begin(); i != compat.end(); i++) {
0254 long int detid = i->first->geographicalId().rawId();
0255
0256 if (!GeomDetEnumerators::isTrackerPixel(tkLayer->subDetector())) {
0257 auto DetMatch = (rphirecHits.product())->find((detid));
0258 auto MDetMatch = (matchedrecHits.product())->find((detid));
0259
0260 long int DetID = (DetMatch != rphirecHits->end()) ? detid : 0;
0261
0262 if ((MDetMatch != matchedrecHits->end()) && !MDetMatch->empty()) {
0263 long int pii = MDetMatch->begin()->monoId();
0264 auto CDetMatch = (rphirecHits.product())->find((pii));
0265 DetID = (CDetMatch != rphirecHits->end()) ? pii : 0;
0266 }
0267
0268 temp.push_back(DetID);
0269
0270 } else {
0271 auto DetMatch = (pixelHits.product())->find((detid));
0272 long int DetID = (DetMatch != pixelHits->end()) ? detid : 0;
0273 temp.push_back(DetID);
0274 }
0275 }
0276
0277 Idd.push_back(temp);
0278
0279 }
0280 if (Idd.size() < 2)
0281 continue;
0282
0283
0284 for (unsigned int i = 0; i < Idd.size() - 2; i++) {
0285 for (unsigned int i1 = 0; i1 < Idd[i].size(); i1++) {
0286 for (unsigned int i2 = 0; i2 < Idd[i + 1].size(); i2++) {
0287 for (unsigned int i3 = 0; i3 < Idd[i + 2].size(); i3++) {
0288 if ((Idd[i][i1] != 0) && (Idd[i + 1][i2] != 0) && (Idd[i + 2][i3] != 0)) {
0289 vector<long int> tmp;
0290 tmp.push_back(Idd[i][i1]);
0291 tmp.push_back(Idd[i + 1][i2]);
0292 tmp.push_back(Idd[i + 2][i3]);
0293
0294 bool newTrip = true;
0295 for (unsigned int iv = 0; iv < tripl.size(); iv++) {
0296 if ((tripl[iv][0] == tmp[0]) && (tripl[iv][1] == tmp[1]) && (tripl[iv][2] == tmp[2]))
0297 newTrip = false;
0298 }
0299 if (newTrip) {
0300 tripl.push_back(tmp);
0301 }
0302 }
0303 }
0304 }
0305 }
0306 }
0307 }
0308
0309 float sineta_brem = sinh(eta_br);
0310
0311
0312 auto bfield = iSetup.getHandle(magFieldToken_);
0313 float nomField = bfield->nominalValue();
0314
0315 TransientTrackingRecHit::ConstRecHitContainer glob_hits;
0316 OwnVector<TrackingRecHit> loc_hits;
0317 for (unsigned int i = 0; i < tripl.size(); i++) {
0318 auto DetMatch1 = (rphirecHits.product())->find(tripl[i][0]);
0319 auto DetMatch2 = (rphirecHits.product())->find(tripl[i][1]);
0320 auto DetMatch3 = (rphirecHits.product())->find(tripl[i][2]);
0321 if ((DetMatch1 == rphirecHits->end()) || (DetMatch2 == rphirecHits->end()) || (DetMatch3 == rphirecHits->end()))
0322 continue;
0323 auto DetSet1 = *DetMatch1;
0324 auto DetSet2 = *DetMatch2;
0325 auto DetSet3 = *DetMatch3;
0326
0327 for (auto it1 = DetSet1.begin(); it1 != DetSet1.end(); ++it1) {
0328 GlobalPoint gp1 = tracker_->idToDet(tripl[i][0])->surface().toGlobal(it1->localPosition());
0329
0330 bool tak1 = isGsfTrack(gsfRecHits, &(*it1));
0331
0332 for (auto it2 = DetSet2.begin(); it2 != DetSet2.end(); ++it2) {
0333 GlobalPoint gp2 = tracker_->idToDet(tripl[i][1])->surface().toGlobal(it2->localPosition());
0334 bool tak2 = isGsfTrack(gsfRecHits, &(*it2));
0335
0336 for (auto it3 = DetSet3.begin(); it3 != DetSet3.end(); ++it3) {
0337
0338 GlobalPoint gp3 = tracker_->idToDet(tripl[i][2])->surface().toGlobal(it3->localPosition());
0339 bool tak3 = isGsfTrack(gsfRecHits, &(*it3));
0340
0341 FastHelix helix(gp3, gp2, gp1, nomField, &*bfield);
0342 GlobalVector gv = helix.stateAtVertex().momentum();
0343 GlobalVector gv_corr(gv.x(), gv.y(), gv.perp() * sineta_brem);
0344 float ene = sqrt(gv_corr.mag2() + (pfmass * pfmass));
0345
0346 GlobalPoint gp = helix.stateAtVertex().position();
0347 float ch = helix.stateAtVertex().charge();
0348
0349 XYZTLorentzVector mom = XYZTLorentzVector(gv.x(), gv.y(), gv_corr.z(), ene);
0350 XYZTLorentzVector pos = XYZTLorentzVector(gp.x(), gp.y(), gp.z(), 0.);
0351 BaseParticlePropagator theOutParticle(RawParticle(mom, pos, ch), 0, 0, B_.z());
0352 int bgc = GoodCluster(theOutParticle, PPP, 0.3, true);
0353
0354 if (gv.perp() < 0.5)
0355 continue;
0356
0357 if (tak1 + tak2 + tak3 > 2)
0358 continue;
0359
0360 if (bgc == -1)
0361 continue;
0362 bool clTak = false;
0363 for (unsigned int igcc = 0; igcc < gc.size(); igcc++) {
0364 if (clTak)
0365 continue;
0366 if (bgc == gc[igcc])
0367 clTak = true;
0368 }
0369 if (clTak)
0370 continue;
0371
0372 GlobalTrajectoryParameters Gtp(gp1, gv, int(ch), &(*magfield_));
0373 glob_hits.clear();
0374 loc_hits.clear();
0375 glob_hits.push_back(hitBuilder_->build(it1->clone()));
0376 glob_hits.push_back(hitBuilder_->build(it2->clone()));
0377 glob_hits.push_back(hitBuilder_->build(it3->clone()));
0378
0379
0380
0381 FreeTrajectoryState CSeed(Gtp, CurvilinearTrajectoryError(AlgebraicSymMatrix55(AlgebraicMatrixID())));
0382 TrajectoryStateOnSurface updatedState;
0383
0384 for (int ih = 0; ih < 3; ih++) {
0385 TrajectoryStateOnSurface state =
0386 (ih == 0) ? propagator_->propagate(CSeed, tracker_->idToDet(tripl[i][ih])->surface())
0387 : propagator_->propagate(updatedState, tracker_->idToDet(tripl[i][ih])->surface());
0388
0389 if (!state.isValid()) {
0390 ih = 3;
0391 continue;
0392 }
0393
0394 updatedState = kfUpdator_->update(state, *glob_hits[ih]);
0395 loc_hits.push_back(glob_hits[ih]->hit()->clone());
0396 if (ih == 2) {
0397 PTrajectoryStateOnDet const& PTraj =
0398 trajectoryStateTransform::persistentState(updatedState, tripl[i][2]);
0399
0400 unclean.push_back(make_pair(TrajectorySeed(PTraj, loc_hits, alongMomentum), make_pair(gv_corr, ch)));
0401 }
0402
0403 }
0404 }
0405 }
0406 }
0407 }
0408 vector<bool> inPhot = sharedHits(unclean);
0409 for (unsigned int iu = 0; iu < unclean.size(); iu++) {
0410 if (inPhot[iu])
0411 output->push_back(ConvBremSeed(unclean[iu].first, pft));
0412 }
0413
0414 }
0415 LogDebug("ConvBremSeedProducerProducer") << "END";
0416 iEvent.put(std::move(output));
0417 }
0418
0419 void ConvBremSeedProducer::beginRun(const edm::Run& run, const EventSetup& iSetup) {
0420 geomSearchTracker_ = &iSetup.getData(geomSearchTrackerToken_);
0421
0422 geometry_ = &iSetup.getData(geometryToken_);
0423
0424 tracker_ = &iSetup.getData(trackerToken_);
0425
0426 magfield_ = &iSetup.getData(magFieldToken_beginRun_);
0427 B_ = magfield_->inTesla(GlobalPoint(0, 0, 0));
0428
0429 fieldMap_ = &iSetup.getData(magFieldMapToken_);
0430
0431 hitBuilder_ = &iSetup.getData(hitBuilderToken_);
0432
0433 propagator_ = new PropagatorWithMaterial(alongMomentum, 0.0005, magfield_);
0434 kfUpdator_ = new KFUpdator();
0435 }
0436
0437 void ConvBremSeedProducer::endRun(const edm::Run& run, const EventSetup& iSetup) {
0438 delete propagator_;
0439 delete kfUpdator_;
0440 }
0441
0442 void ConvBremSeedProducer::initializeLayerMap() {
0443
0444
0445
0446
0447
0448
0449
0450
0451
0452
0453
0454 const std::vector<const BarrelDetLayer*>& barrelLayers = geomSearchTracker_->barrelLayers();
0455 LogDebug("FastTracker") << "Barrel DetLayer dump: ";
0456 for (auto bl = barrelLayers.begin(); bl != barrelLayers.end(); ++bl) {
0457 LogDebug("FastTracker") << "radius " << (**bl).specificSurface().radius();
0458 }
0459
0460 const std::vector<const ForwardDetLayer*>& posForwardLayers = geomSearchTracker_->posForwardLayers();
0461 LogDebug("FastTracker") << "Positive Forward DetLayer dump: ";
0462 for (auto fl = posForwardLayers.begin(); fl != posForwardLayers.end(); ++fl) {
0463 LogDebug("FastTracker") << "Z pos " << (**fl).surface().position().z() << " radii "
0464 << (**fl).specificSurface().innerRadius() << ", " << (**fl).specificSurface().outerRadius();
0465 }
0466
0467 const float rTolerance = 1.5;
0468 const float zTolerance = 3.;
0469
0470 LogDebug("FastTracker") << "Dump of TrackerInteractionGeometry cylinders:";
0471 for (std::list<TrackerLayer>::const_iterator i = geometry_->cylinderBegin(); i != geometry_->cylinderEnd(); ++i) {
0472 const BoundCylinder* cyl = i->cylinder();
0473 const BoundDisk* disk = i->disk();
0474
0475 LogDebug("FastTracker") << "Famos Layer no " << i->layerNumber() << " is sensitive? " << i->sensitive() << " pos "
0476 << i->surface().position();
0477 if (!i->sensitive())
0478 continue;
0479
0480 if (cyl != nullptr) {
0481 LogDebug("FastTracker") << " cylinder radius " << cyl->radius();
0482 bool found = false;
0483
0484 for (auto bl = barrelLayers.begin(); bl != barrelLayers.end(); ++bl) {
0485 if (fabs(cyl->radius() - (**bl).specificSurface().radius()) < rTolerance) {
0486 layerMap_[i->layerNumber()] = *bl;
0487 found = true;
0488 LogDebug("FastTracker") << "Corresponding DetLayer found with radius " << (**bl).specificSurface().radius();
0489
0490 break;
0491 }
0492 }
0493 if (!found) {
0494 LogError("FastTracker") << "FAILED to find a corresponding DetLayer!";
0495 }
0496 } else {
0497 LogDebug("FastTracker") << " disk radii " << disk->innerRadius() << ", " << disk->outerRadius();
0498
0499 bool found = false;
0500
0501 for (auto fl = posForwardLayers.begin(); fl != posForwardLayers.end(); ++fl) {
0502 if (fabs(disk->position().z() - (**fl).surface().position().z()) < zTolerance) {
0503 layerMap_[i->layerNumber()] = *fl;
0504 found = true;
0505 LogDebug("FastTracker") << "Corresponding DetLayer found with Z pos " << (**fl).surface().position().z()
0506 << " and radii " << (**fl).specificSurface().innerRadius() << ", "
0507 << (**fl).specificSurface().outerRadius();
0508 break;
0509 }
0510 }
0511 if (!found) {
0512 LogError("FastTracker") << "FAILED to find a corresponding DetLayer!";
0513 }
0514 }
0515 }
0516 }
0517 const DetLayer* ConvBremSeedProducer::detLayer(const TrackerLayer& layer, float zpos) const {
0518 if (zpos > 0 || !layer.forward())
0519 return layerMap_[layer.layerNumber()];
0520 else
0521 return layerMap_[layer.layerNumber() + negLayerOffset_];
0522 }
0523
0524 TrajectoryStateOnSurface ConvBremSeedProducer::makeTrajectoryState(const DetLayer* layer,
0525 const ParticlePropagator& pp,
0526 const MagneticField* field) const {
0527 GlobalPoint pos(pp.particle().X(), pp.particle().Y(), pp.particle().Z());
0528 GlobalVector mom(pp.particle().Px(), pp.particle().Py(), pp.particle().Pz());
0529
0530 auto plane = layer->surface().tangentPlane(pos);
0531
0532 return TrajectoryStateOnSurface(GlobalTrajectoryParameters(pos, mom, TrackCharge(pp.particle().charge()), field),
0533 *plane);
0534 }
0535 bool ConvBremSeedProducer::isGsfTrack(const reco::Track& tkv, const TrackingRecHit* h) {
0536 bool istaken = false;
0537 for (auto const& hit : tkv.recHits()) {
0538 if (istaken || !hit->isValid())
0539 continue;
0540 istaken = hit->sharesInput(h, TrackingRecHit::all);
0541 }
0542 return istaken;
0543 }
0544 vector<bool> ConvBremSeedProducer::sharedHits(const vector<pair<TrajectorySeed, pair<GlobalVector, float> > >& unclean) {
0545 vector<bool> goodseed;
0546 goodseed.clear();
0547 if (unclean.size() < 2) {
0548 for (unsigned int i = 0; i < unclean.size(); i++)
0549 goodseed.push_back(true);
0550 } else {
0551 for (unsigned int i = 0; i < unclean.size(); i++)
0552 goodseed.push_back(true);
0553
0554 for (unsigned int iu = 0; iu < unclean.size() - 1; iu++) {
0555 if (!goodseed[iu])
0556 continue;
0557 for (unsigned int iu2 = iu + 1; iu2 < unclean.size(); iu2++) {
0558 if (!goodseed[iu])
0559 continue;
0560 if (!goodseed[iu2])
0561 continue;
0562
0563
0564 unsigned int shar = 0;
0565 for (auto const& sh : unclean[iu].first.recHits()) {
0566 for (auto const& sh2 : unclean[iu2].first.recHits()) {
0567 if (sh.sharesInput(&sh2, TrackingRecHit::all))
0568
0569 shar++;
0570 }
0571 }
0572 if (shar >= 2) {
0573 if (unclean[iu].second.first.perp() < unclean[iu2].second.first.perp())
0574 goodseed[iu] = false;
0575 else
0576 goodseed[iu2] = false;
0577 }
0578 }
0579 }
0580 }
0581 return goodseed;
0582 }
0583
0584 int ConvBremSeedProducer::GoodCluster(const BaseParticlePropagator& ubpg,
0585 const PFClusterCollection& pfc,
0586 float minep,
0587 bool sec) {
0588 BaseParticlePropagator bpg = ubpg;
0589 bpg.propagateToEcalEntrance(false);
0590 float dr = 1000;
0591 float de = 1000;
0592 float df = 1000;
0593 int ibest = -1;
0594
0595 if (bpg.getSuccess() != 0) {
0596 for (unsigned int i = 0; i < pfc.size(); i++) {
0597 float tmp_ep = pfc[i].energy() / bpg.particle().momentum().e();
0598 float tmp_phi = fabs(pfc[i].position().phi() - bpg.particle().vertex().phi());
0599 if (tmp_phi > TMath::TwoPi())
0600 tmp_phi -= TMath::TwoPi();
0601 float tmp_eta = fabs(pfc[i].position().eta() - bpg.particle().vertex().eta());
0602 float tmp_dr = sqrt(pow(tmp_phi, 2) + pow(tmp_eta, 2));
0603 bool isBet = (tmp_dr < dr);
0604 if (sec)
0605 isBet = (tmp_phi < df);
0606 if ((isBet) && (tmp_ep > minep) && (tmp_ep < 10)) {
0607 dr = tmp_dr;
0608 de = tmp_eta;
0609 df = tmp_phi;
0610 ibest = i;
0611 }
0612 }
0613 bool isBad = (dr > 0.1);
0614 if (sec)
0615 isBad = ((df > 0.25) || (de > 0.5));
0616
0617 if (isBad)
0618 ibest = -1;
0619 }
0620 return ibest;
0621 }
