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