Back to home page

Project CMSSW displayed by LXR

 
 

    


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   ///B field
0092   math::XYZVector B_;
0093 
0094   // Event setup tokens
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   // convBremSeeds
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   ///INPUT COLLECTIONS
0144 
0145   ///PF CLUSTERS
0146   Handle<PFClusterCollection> PfC;
0147   iEvent.getByLabel(conf_.getParameter<InputTag>("PFClusters"), PfC);
0148   const PFClusterCollection& PPP = *(PfC.product());
0149 
0150   ///PIXEL
0151   Handle<SiPixelRecHitCollection> pixelHits;
0152   iEvent.getByLabel(conf_.getParameter<InputTag>("pixelRecHits"), pixelHits);
0153 
0154   ///STRIP
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   //GSFPFRECTRACKS
0161   Handle<GsfPFRecTrackCollection> thePfRecTrackCollection;
0162   iEvent.getByLabel(conf_.getParameter<InputTag>("PFRecTrackLabel"), thePfRecTrackCollection);
0163   const GsfPFRecTrackCollection PfRTkColl = *(thePfRecTrackCollection.product());
0164 
0165   ///OUTPUT COLLECTION
0166   auto output = std::make_unique<ConvBremSeedCollection>();
0167 
0168   ///INITIALIZE
0169   vector<pair<TrajectorySeed, pair<GlobalVector, float> > > unclean;
0170   //TRIPLET OF MODULES TO BE USED FOR SEEDING
0171   vector<vector<long int> > tripl;
0172   //LAYER MAP
0173   initializeLayerMap();
0174 
0175   ///LOOP OVER GSF TRACK COLLECTION
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     ///FIND THE CLUSTER ASSOCIATED TO THE GSF TRACK
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     ///LOOP OVER BREM PHOTONS
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       ///BREM SELECTION
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       ///LOOP OVER TRACKER LAYER
0229       list<TrackerLayer>::const_iterator cyliter = geometry_->cylinderBegin();
0230       for (; cyliter != geometry_->cylinderEnd(); ++cyliter) {
0231         ///TRACKER LAYER SELECTION
0232         if (!(cyliter->sensitive()))
0233           continue;
0234         PP.setPropagationConditions(*cyliter);
0235         PP.propagate();
0236         if (PP.getSuccess() == 0)
0237           continue;
0238 
0239         ///FIND COMPATIBLE MODULES
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       }  //END TRACKER LAYER LOOP
0280       if (Idd.size() < 2)
0281         continue;
0282 
0283       ///MODULE TRIPLETS SELECTION
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     }  //END BREM LOOP
0308 
0309     float sineta_brem = sinh(eta_br);
0310 
0311     //OUTPUT COLLECTION
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             //  ips++;
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             ///SEED CREATION
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                 //      output->push_back(Trajectoryseed(PTraj,loc_hits,alongMomentum));
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   }  //END GSF TRACK COLLECTION LOOP
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   // These are the BoundSurface&, the BoundDisk* and the BoundCylinder* for that layer
0444   //   const BoundSurface& theSurface = layer.surface();
0445   //   BoundDisk* theDisk = layer.disk();  // non zero for endcaps
0446   //   BoundCylinder* theCylinder = layer.cylinder(); // non zero for barrel
0447   //   int theLayer = layer.layerNumber(); // 1->3 PixB, 4->5 PixD,
0448   //                                       // 6->9 TIB, 10->12 TID,
0449   //                                       // 13->18 TOB, 19->27 TEC
0450 
0451   /// ATTENTION: HARD CODED LOGIC! If Famos layer numbering changes this logic needs to
0452   /// be adapted to the new numbering!
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         //    if (unclean[iu].second.second *unclean[iu2].second.second >0)continue;
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 }