File indexing completed on 2021-11-05 00:08:12
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019
0020
0021
0022 #include "CommonTools/BaseParticlePropagator/interface/BaseParticlePropagator.h"
0023 #include "CommonTools/MVAUtils/interface/GBRForestTools.h"
0024 #include "CondFormats/GBRForest/interface/GBRForest.h"
0025 #include "DataFormats/Common/interface/ValueMap.h"
0026 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0027 #include "DataFormats/Math/interface/deltaR.h"
0028 #include "DataFormats/Math/interface/LorentzVector.h"
0029 #include "DataFormats/ParticleFlowReco/interface/PFCluster.h"
0030 #include "DataFormats/ParticleFlowReco/interface/PreId.h"
0031 #include "DataFormats/ParticleFlowReco/interface/PreIdFwd.h"
0032 #include "DataFormats/TrackReco/interface/Track.h"
0033 #include "DataFormats/TrajectorySeed/interface/PropagationDirection.h"
0034 #include "DataFormats/TrajectorySeed/interface/TrajectorySeedCollection.h"
0035 #include "FWCore/Framework/interface/ESHandle.h"
0036 #include "FWCore/Framework/interface/Event.h"
0037 #include "FWCore/Framework/interface/MakerMacros.h"
0038 #include "FWCore/Framework/interface/stream/EDProducer.h"
0039 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0040 #include "FWCore/ParameterSet/interface/FileInPath.h"
0041 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0042 #include "MagneticField/Engine/interface/MagneticField.h"
0043 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0044 #include "RecoParticleFlow/PFClusterTools/interface/PFResolutionMap.h"
0045 #include "RecoParticleFlow/PFTracking/interface/PFGeometry.h"
0046 #include "RecoParticleFlow/PFTracking/interface/PFTrackTransformer.h"
0047 #include "RecoTracker/TransientTrackingRecHit/interface/TkTransientTrackingRecHitBuilder.h"
0048 #include "TrackingTools/PatternTools/interface/Trajectory.h"
0049 #include "TrackingTools/PatternTools/interface/TrajectorySmoother.h"
0050 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0051 #include "TrackingTools/Records/interface/TransientRecHitRecord.h"
0052 #include "TrackingTools/TrackFitters/interface/TrajectoryFitter.h"
0053 #include "TrackingTools/TrajectoryState/interface/TrajectoryStateTransform.h"
0054
0055 #include "TMath.h"
0056 #include "Math/VectorUtil.h"
0057
0058 #include <fstream>
0059
0060 namespace goodseedhelpers {
0061 struct HeavyObjectCache {
0062 constexpr static unsigned int kMaxWeights = 9;
0063 HeavyObjectCache(const edm::ParameterSet& conf);
0064 std::array<std::unique_ptr<const GBRForest>, kMaxWeights> gbr;
0065 };
0066 }
0067
0068 class GoodSeedProducer final : public edm::stream::EDProducer<edm::GlobalCache<goodseedhelpers::HeavyObjectCache>> {
0069 typedef TrajectoryStateOnSurface TSOS;
0070
0071 public:
0072 explicit GoodSeedProducer(const edm::ParameterSet&, const goodseedhelpers::HeavyObjectCache*);
0073
0074 static std::unique_ptr<goodseedhelpers::HeavyObjectCache> initializeGlobalCache(const edm::ParameterSet& conf) {
0075 return std::make_unique<goodseedhelpers::HeavyObjectCache>(conf);
0076 }
0077
0078 static void globalEndJob(goodseedhelpers::HeavyObjectCache const*) {}
0079
0080 private:
0081 void beginRun(const edm::Run& run, const edm::EventSetup&) override;
0082 void produce(edm::Event&, const edm::EventSetup&) override;
0083
0084
0085 int getBin(float, float);
0086
0087 void fillPreIdRefValueMap(edm::Handle<reco::TrackCollection> tkhandle,
0088 const edm::OrphanHandle<reco::PreIdCollection>&,
0089 edm::ValueMap<reco::PreIdRef>::Filler& filler);
0090
0091
0092
0093 std::string preidckf_;
0094
0095
0096 std::string preidgsf_;
0097
0098
0099 std::string preidname_;
0100
0101
0102 TkClonerImpl hitCloner;
0103
0104
0105 std::unique_ptr<PFTrackTransformer> pfTransformer_;
0106
0107
0108 int nHitsInSeed_;
0109
0110
0111 double minPt_;
0112 double maxPt_;
0113 double maxEta_;
0114
0115 double HcalIsolWindow_;
0116 double EcalStripSumE_minClusEnergy_;
0117 double EcalStripSumE_deltaEta_;
0118 double EcalStripSumE_deltaPhiOverQ_minValue_;
0119 double EcalStripSumE_deltaPhiOverQ_maxValue_;
0120 double minEoverP_;
0121 double maxHoverP_;
0122
0123
0124 double clusThreshold_;
0125
0126
0127 double minEp_;
0128 double maxEp_;
0129
0130
0131 bool produceCkfseed_;
0132
0133
0134 bool disablePreId_;
0135
0136
0137 bool producePreId_;
0138
0139
0140 double PtThresholdSavePredId_;
0141
0142
0143 float thr[150];
0144
0145
0146 edm::ParameterSet conf_;
0147 edm::EDGetTokenT<reco::PFClusterCollection> pfCLusTagPSLabel_;
0148 edm::EDGetTokenT<reco::PFClusterCollection> pfCLusTagECLabel_;
0149 edm::EDGetTokenT<reco::PFClusterCollection> pfCLusTagHCLabel_;
0150 std::vector<edm::EDGetTokenT<std::vector<Trajectory>>> trajContainers_;
0151 std::vector<edm::EDGetTokenT<reco::TrackCollection>> tracksContainers_;
0152
0153 std::string fitterName_;
0154 std::string smootherName_;
0155 std::string propagatorName_;
0156 std::string trackerRecHitBuilderName_;
0157
0158 std::unique_ptr<PFResolutionMap> resMapEtaECAL_;
0159 std::unique_ptr<PFResolutionMap> resMapPhiECAL_;
0160
0161 const edm::ESGetToken<TrajectoryFitter, TrajectoryFitter::Record> fitterToken_;
0162 const edm::ESGetToken<TrajectorySmoother, TrajectoryFitter::Record> smootherToken_;
0163 const edm::ESGetToken<TransientTrackingRecHitBuilder, TransientRecHitRecord> trackerRecHitBuilderToken_;
0164 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldToken_;
0165 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> magneticFieldTokenBeginRun_;
0166
0167
0168 bool useQuality_;
0169 reco::TrackBase::TrackQuality trackQuality_;
0170
0171
0172 float nhit, dpt, chired, chiRatio;
0173 float chikfred, trk_ecalDeta, trk_ecalDphi;
0174 double Min_dr_;
0175
0176
0177 bool useTmva_;
0178
0179
0180 math::XYZVector B_;
0181
0182
0183 std::map<reco::TrackRef, unsigned> refMap_;
0184 };
0185
0186 using namespace edm;
0187 using namespace std;
0188 using namespace reco;
0189
0190 GoodSeedProducer::GoodSeedProducer(const ParameterSet& iConfig, const goodseedhelpers::HeavyObjectCache*)
0191 : pfTransformer_(nullptr),
0192 conf_(iConfig),
0193 resMapEtaECAL_(nullptr),
0194 resMapPhiECAL_(nullptr),
0195 fitterToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<string>("Fitter")))),
0196 smootherToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<string>("Smoother")))),
0197 trackerRecHitBuilderToken_(esConsumes(edm::ESInputTag("", iConfig.getParameter<std::string>("TTRHBuilder")))),
0198 magneticFieldToken_(esConsumes()),
0199 magneticFieldTokenBeginRun_(esConsumes<edm::Transition::BeginRun>()) {
0200 LogInfo("GoodSeedProducer") << "Electron PreIdentification started ";
0201
0202
0203 std::vector<edm::InputTag> tags = iConfig.getParameter<vector<InputTag>>("TkColList");
0204 for (unsigned int i = 0; i < tags.size(); ++i) {
0205 trajContainers_.push_back(consumes<vector<Trajectory>>(tags[i]));
0206 tracksContainers_.push_back(consumes<reco::TrackCollection>(tags[i]));
0207 }
0208
0209 minPt_ = iConfig.getParameter<double>("MinPt");
0210 maxPt_ = iConfig.getParameter<double>("MaxPt");
0211 maxEta_ = iConfig.getParameter<double>("MaxEta");
0212
0213 HcalIsolWindow_ = iConfig.getParameter<double>("HcalWindow");
0214 EcalStripSumE_minClusEnergy_ = iConfig.getParameter<double>("EcalStripSumE_minClusEnergy");
0215 EcalStripSumE_deltaEta_ = iConfig.getParameter<double>("EcalStripSumE_deltaEta");
0216 EcalStripSumE_deltaPhiOverQ_minValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_minValue");
0217 EcalStripSumE_deltaPhiOverQ_maxValue_ = iConfig.getParameter<double>("EcalStripSumE_deltaPhiOverQ_maxValue");
0218 minEoverP_ = iConfig.getParameter<double>("EOverPLead_minValue");
0219 maxHoverP_ = iConfig.getParameter<double>("HOverPLead_maxValue");
0220
0221 pfCLusTagECLabel_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFEcalClusterLabel"));
0222
0223 pfCLusTagHCLabel_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFHcalClusterLabel"));
0224
0225 pfCLusTagPSLabel_ = consumes<reco::PFClusterCollection>(iConfig.getParameter<InputTag>("PFPSClusterLabel"));
0226
0227 preidgsf_ = iConfig.getParameter<string>("PreGsfLabel");
0228 preidckf_ = iConfig.getParameter<string>("PreCkfLabel");
0229 preidname_ = iConfig.getParameter<string>("PreIdLabel");
0230
0231 fitterName_ = iConfig.getParameter<string>("Fitter");
0232 smootherName_ = iConfig.getParameter<string>("Smoother");
0233
0234 nHitsInSeed_ = iConfig.getParameter<int>("NHitsInSeed");
0235
0236 clusThreshold_ = iConfig.getParameter<double>("ClusterThreshold");
0237
0238 minEp_ = iConfig.getParameter<double>("MinEOverP");
0239 maxEp_ = iConfig.getParameter<double>("MaxEOverP");
0240
0241
0242 produceCkfseed_ = iConfig.getUntrackedParameter<bool>("ProduceCkfSeed", false);
0243
0244
0245 disablePreId_ = iConfig.getUntrackedParameter<bool>("DisablePreId", false);
0246
0247 producePreId_ = iConfig.getUntrackedParameter<bool>("ProducePreId", true);
0248
0249 if (disablePreId_)
0250 producePreId_ = false;
0251 PtThresholdSavePredId_ = iConfig.getUntrackedParameter<double>("PtThresholdSavePreId", 1.);
0252
0253 LogDebug("GoodSeedProducer") << "Seeds for GSF will be produced ";
0254
0255
0256 produces<ElectronSeedCollection>(preidgsf_);
0257
0258 if (produceCkfseed_) {
0259 LogDebug("GoodSeedProducer") << "Seeds for CKF will be produced ";
0260 produces<TrajectorySeedCollection>(preidckf_);
0261 }
0262
0263 if (producePreId_) {
0264 LogDebug("GoodSeedProducer") << "PreId debugging information will be produced ";
0265
0266 produces<PreIdCollection>(preidname_);
0267 if (tracksContainers_.size() == 1)
0268 produces<edm::ValueMap<reco::PreIdRef>>(preidname_);
0269 }
0270
0271 useQuality_ = iConfig.getParameter<bool>("UseQuality");
0272 trackQuality_ = TrackBase::qualityByName(iConfig.getParameter<std::string>("TrackQuality"));
0273
0274 useTmva_ = iConfig.getUntrackedParameter<bool>("UseTMVA", false);
0275
0276 Min_dr_ = iConfig.getParameter<double>("Min_dr");
0277
0278 trackerRecHitBuilderName_ = iConfig.getParameter<std::string>("TTRHBuilder");
0279 }
0280
0281
0282
0283
0284
0285
0286 void GoodSeedProducer::produce(Event& iEvent, const EventSetup& iSetup) {
0287 LogDebug("GoodSeedProducer") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run();
0288
0289 auto output_preid = std::make_unique<ElectronSeedCollection>();
0290 auto output_nopre = std::make_unique<TrajectorySeedCollection>();
0291 auto output_preidinfo = std::make_unique<PreIdCollection>();
0292 auto preIdMap_p = std::make_unique<edm::ValueMap<reco::PreIdRef>>();
0293 edm::ValueMap<reco::PreIdRef>::Filler mapFiller(*preIdMap_p);
0294
0295 std::unique_ptr<TrajectoryFitter> fitter;
0296 std::unique_ptr<TrajectorySmoother> smoother;
0297
0298
0299 if (!disablePreId_) {
0300 auto const& aFitter = &iSetup.getData(fitterToken_);
0301 auto const& aSmoother = &iSetup.getData(smootherToken_);
0302
0303 smoother.reset(aSmoother->clone());
0304 fitter = aFitter->clone();
0305 auto const& theTrackerRecHitBuilder = &iSetup.getData(trackerRecHitBuilderToken_);
0306 hitCloner = static_cast<TkTransientTrackingRecHitBuilder const*>(theTrackerRecHitBuilder)->cloner();
0307 fitter->setHitCloner(&hitCloner);
0308 smoother->setHitCloner(&hitCloner);
0309 }
0310
0311
0312 refMap_.clear();
0313
0314
0315 auto const& magneticField = &iSetup.getData(magneticFieldToken_);
0316
0317
0318
0319 Handle<PFClusterCollection> theECPfClustCollection;
0320 iEvent.getByToken(pfCLusTagECLabel_, theECPfClustCollection);
0321
0322 vector<PFCluster const*> basClus;
0323 for (auto const& klus : *theECPfClustCollection.product()) {
0324 if (klus.correctedEnergy() > clusThreshold_)
0325 basClus.push_back(&klus);
0326 }
0327
0328
0329 Handle<PFClusterCollection> theHCPfClustCollection;
0330 iEvent.getByToken(pfCLusTagHCLabel_, theHCPfClustCollection);
0331
0332
0333 Handle<PFClusterCollection> thePSPfClustCollection;
0334 iEvent.getByToken(pfCLusTagPSLabel_, thePSPfClustCollection);
0335
0336
0337 for (unsigned int istr = 0; istr < tracksContainers_.size(); ++istr) {
0338
0339 Handle<TrackCollection> tkRefCollection;
0340 iEvent.getByToken(tracksContainers_[istr], tkRefCollection);
0341 const TrackCollection& Tk = *(tkRefCollection.product());
0342
0343 LogDebug("GoodSeedProducer") << "Number of tracks in collection "
0344 << "tracksContainers_[" << istr << "] to be analyzed " << Tk.size();
0345
0346
0347 for (unsigned int i = 0; i < Tk.size(); ++i) {
0348 if (useQuality_ && (!(Tk[i].quality(trackQuality_))))
0349 continue;
0350
0351 reco::PreId myPreId;
0352 bool GoodPreId = false;
0353
0354 TrackRef trackRef(tkRefCollection, i);
0355 math::XYZVectorF tkmom(Tk[i].momentum());
0356 auto tketa = tkmom.eta();
0357 auto tkpt = std::sqrt(tkmom.perp2());
0358 auto const& Seed = (*trackRef->seedRef());
0359
0360 if (!disablePreId_) {
0361 int ipteta = getBin(Tk[i].eta(), Tk[i].pt());
0362 int ibin = ipteta * 9;
0363
0364
0365 float oPTOB = 1.f / std::sqrt(Tk[i].innerMomentum().mag2());
0366
0367 float nchi = Tk[i].normalizedChi2();
0368
0369 int nhitpi = Tk[i].found();
0370 float EP = 0;
0371
0372
0373 myPreId.setTrack(trackRef);
0374
0375
0376 auto pfmass = 0.0005;
0377 auto pfoutenergy = sqrt((pfmass * pfmass) + Tk[i].outerMomentum().Mag2());
0378
0379 XYZTLorentzVector mom = XYZTLorentzVector(
0380 Tk[i].outerMomentum().x(), Tk[i].outerMomentum().y(), Tk[i].outerMomentum().z(), pfoutenergy);
0381 XYZTLorentzVector pos =
0382 XYZTLorentzVector(Tk[i].outerPosition().x(), Tk[i].outerPosition().y(), Tk[i].outerPosition().z(), 0.);
0383
0384 BaseParticlePropagator theOutParticle(RawParticle(mom, pos, Tk[i].charge()), 0, 0, B_.z());
0385
0386 theOutParticle.propagateToEcalEntrance(false);
0387
0388 float toteta = 1000.f;
0389 float totphi = 1000.f;
0390 float dr = 1000.f;
0391 float EE = 0.f;
0392 float feta = 0.f;
0393 GlobalPoint ElecTrkEcalPos(0, 0, 0);
0394
0395 PFClusterRef clusterRef;
0396 math::XYZPoint meanShowerSaved;
0397 if (theOutParticle.getSuccess() != 0) {
0398 ElecTrkEcalPos = GlobalPoint(theOutParticle.particle().vertex().x(),
0399 theOutParticle.particle().vertex().y(),
0400 theOutParticle.particle().vertex().z());
0401
0402 constexpr float psLim = 2.50746495928f;
0403 bool isBelowPS = (ElecTrkEcalPos.z() * ElecTrkEcalPos.z()) > (psLim * psLim) * ElecTrkEcalPos.perp2();
0404
0405
0406 unsigned clusCounter = 0;
0407 float max_ee = 0;
0408 for (auto aClus : basClus) {
0409 float tmp_ep = float(aClus->correctedEnergy()) * oPTOB;
0410 if ((tmp_ep < minEp_) | (tmp_ep > maxEp_)) {
0411 ++clusCounter;
0412 continue;
0413 }
0414
0415 double ecalShowerDepth = PFCluster::getDepthCorrection(aClus->correctedEnergy(), isBelowPS, false);
0416 auto mom = theOutParticle.particle().momentum().Vect();
0417 auto meanShower = ElecTrkEcalPos + GlobalVector(mom.x(), mom.y(), mom.z()).unit() * ecalShowerDepth;
0418
0419 float etarec = meanShower.eta();
0420 float phirec = meanShower.phi();
0421
0422 float tmp_phi = std::abs(aClus->positionREP().phi() - phirec);
0423 if (tmp_phi > float(TMath::Pi()))
0424 tmp_phi -= float(TMath::TwoPi());
0425
0426 float tmp_dr = std::sqrt(std::pow(tmp_phi, 2.f) + std::pow(aClus->positionREP().eta() - etarec, 2.f));
0427
0428 if (tmp_dr < dr) {
0429 dr = tmp_dr;
0430 if (dr < Min_dr_) {
0431 if (aClus->correctedEnergy() > max_ee) {
0432 toteta = aClus->positionREP().eta() - etarec;
0433 totphi = tmp_phi;
0434 EP = tmp_ep;
0435 EE = aClus->correctedEnergy();
0436 feta = aClus->positionREP().eta();
0437 clusterRef = PFClusterRef(theECPfClustCollection, clusCounter);
0438 meanShowerSaved = meanShower;
0439 }
0440 }
0441 }
0442 ++clusCounter;
0443 }
0444 }
0445 float trk_ecalDeta_ = fabs(toteta);
0446 float trk_ecalDphi_ = fabs(totphi);
0447
0448
0449 auto ecaletares = resMapEtaECAL_->GetBinContent(resMapEtaECAL_->FindBin(feta, EE));
0450 auto ecalphires = resMapPhiECAL_->GetBinContent(resMapPhiECAL_->FindBin(feta, EE));
0451
0452
0453 float chieta = (toteta != 1000.f) ? toteta / ecaletares : toteta;
0454 float chiphi = (totphi != 1000.f) ? totphi / ecalphires : totphi;
0455 float chichi = sqrt(chieta * chieta + chiphi * chiphi);
0456
0457
0458 float eta_cut = thr[ibin + 0];
0459 float phi_cut = thr[ibin + 1];
0460 float ep_cutmin = thr[ibin + 2];
0461 bool GoodMatching =
0462 ((trk_ecalDeta_ < eta_cut) && (trk_ecalDphi_ < phi_cut) && (EP > ep_cutmin) && (nhitpi > 10));
0463
0464 bool EcalMatching = GoodMatching;
0465
0466 if (tkpt > maxPt_)
0467 GoodMatching = true;
0468 if (tkpt < minPt_)
0469 GoodMatching = false;
0470
0471 math::XYZPoint myPoint(ElecTrkEcalPos.x(), ElecTrkEcalPos.y(), ElecTrkEcalPos.z());
0472 myPreId.setECALMatchingProperties(
0473 clusterRef, myPoint, meanShowerSaved, std::abs(toteta), std::abs(totphi), chieta, chiphi, chichi, EP);
0474 myPreId.setECALMatching(EcalMatching);
0475
0476 bool GoodRange = ((std::abs(tketa) < maxEta_) & (tkpt > minPt_));
0477
0478 int hit1max = int(thr[ibin + 3]);
0479 float chiredmin = thr[ibin + 4];
0480 bool GoodKFFiltering = ((nchi > chiredmin) | (nhitpi < hit1max));
0481
0482 myPreId.setTrackFiltering(GoodKFFiltering);
0483
0484 bool GoodTkId = false;
0485
0486 if ((!GoodMatching) && (GoodKFFiltering) && (GoodRange)) {
0487 chired = 1000;
0488 chiRatio = 1000;
0489 dpt = 0;
0490 nhit = nhitpi;
0491 chikfred = nchi;
0492 trk_ecalDeta = trk_ecalDeta_;
0493 trk_ecalDphi = trk_ecalDphi_;
0494
0495 Trajectory::ConstRecHitContainer tmp;
0496 for (auto const& hit : Tk[i].recHits())
0497 tmp.push_back(hit->cloneSH());
0498 auto const& theTrack = Tk[i];
0499 GlobalVector gv(theTrack.innerMomentum().x(), theTrack.innerMomentum().y(), theTrack.innerMomentum().z());
0500 GlobalPoint gp(theTrack.innerPosition().x(), theTrack.innerPosition().y(), theTrack.innerPosition().z());
0501 GlobalTrajectoryParameters gtps(gp, gv, theTrack.charge(), &*magneticField);
0502 TrajectoryStateOnSurface tsos(gtps, theTrack.innerStateCovariance(), *tmp[0]->surface());
0503 Trajectory&& FitTjs = fitter->fitOne(Seed, tmp, tsos);
0504
0505 if (FitTjs.isValid()) {
0506 Trajectory&& SmooTjs = smoother->trajectory(FitTjs);
0507 if (SmooTjs.isValid()) {
0508
0509
0510 float pt_out = SmooTjs.firstMeasurement().updatedState().globalMomentum().perp();
0511 float pt_in = SmooTjs.lastMeasurement().updatedState().globalMomentum().perp();
0512 dpt = (pt_in > 0) ? fabs(pt_out - pt_in) / pt_in : 0.;
0513
0514 chiRatio = SmooTjs.chiSquared() / Tk[i].chi2();
0515 chired = chiRatio * chikfred;
0516 }
0517 }
0518
0519
0520 if (useTmva_) {
0521 float vars[10] = {nhit, chikfred, dpt, EP, chiRatio, chired, trk_ecalDeta, trk_ecalDphi, tkpt, tketa};
0522
0523 float Ytmva = globalCache()->gbr[ipteta]->GetClassifier(vars);
0524
0525 float BDTcut = thr[ibin + 5];
0526 if (Ytmva > BDTcut)
0527 GoodTkId = true;
0528 myPreId.setMVA(GoodTkId, Ytmva);
0529 myPreId.setTrackProperties(chired, chiRatio, dpt);
0530 } else {
0531 float chiratiocut = thr[ibin + 6];
0532 float gschicut = thr[ibin + 7];
0533 float gsptmin = thr[ibin + 8];
0534
0535 GoodTkId = ((dpt > gsptmin) & (chired < gschicut) & (chiRatio < chiratiocut));
0536 }
0537 }
0538
0539 GoodPreId = GoodTkId | GoodMatching;
0540
0541 myPreId.setFinalDecision(GoodPreId);
0542
0543 #ifdef EDM_ML_DEBUG
0544 if (GoodPreId)
0545 LogDebug("GoodSeedProducer") << "Track (pt= " << Tk[i].pt() << "GeV/c, eta= " << Tk[i].eta()
0546 << ") preidentified for agreement between track and ECAL cluster";
0547 if (GoodPreId && (!GoodMatching))
0548 LogDebug("GoodSeedProducer") << "Track (pt= " << Tk[i].pt() << "GeV/c, eta= " << Tk[i].eta()
0549 << ") preidentified only for track properties";
0550 #endif
0551
0552 }
0553
0554 if (GoodPreId) {
0555
0556 ElectronSeed NewSeed(Seed);
0557 NewSeed.setCtfTrack(trackRef);
0558 output_preid->push_back(NewSeed);
0559 } else {
0560 if (produceCkfseed_) {
0561 output_nopre->push_back(Seed);
0562 }
0563 }
0564 if (producePreId_ && myPreId.pt() > PtThresholdSavePredId_) {
0565
0566 refMap_[trackRef] = output_preidinfo->size();
0567 output_preidinfo->push_back(myPreId);
0568 }
0569 }
0570 }
0571
0572
0573 iEvent.put(std::move(output_preid), preidgsf_);
0574 if (produceCkfseed_)
0575 iEvent.put(std::move(output_nopre), preidckf_);
0576 if (producePreId_) {
0577 const edm::OrphanHandle<reco::PreIdCollection> preIdRefProd = iEvent.put(std::move(output_preidinfo), preidname_);
0578
0579 if (tracksContainers_.size() == 1) {
0580 Handle<TrackCollection> tkRefCollection;
0581 iEvent.getByToken(tracksContainers_[0], tkRefCollection);
0582 fillPreIdRefValueMap(tkRefCollection, preIdRefProd, mapFiller);
0583 mapFiller.fill();
0584 iEvent.put(std::move(preIdMap_p), preidname_);
0585 }
0586 }
0587
0588 refMap_.clear();
0589 }
0590
0591
0592 namespace goodseedhelpers {
0593 HeavyObjectCache::HeavyObjectCache(const edm::ParameterSet& conf) {
0594
0595 const bool useTmva = conf.getUntrackedParameter<bool>("UseTMVA", false);
0596
0597 if (useTmva) {
0598 std::array<edm::FileInPath, kMaxWeights> weights = {{edm::FileInPath(conf.getParameter<string>("Weights1")),
0599 edm::FileInPath(conf.getParameter<string>("Weights2")),
0600 edm::FileInPath(conf.getParameter<string>("Weights3")),
0601 edm::FileInPath(conf.getParameter<string>("Weights4")),
0602 edm::FileInPath(conf.getParameter<string>("Weights5")),
0603 edm::FileInPath(conf.getParameter<string>("Weights6")),
0604 edm::FileInPath(conf.getParameter<string>("Weights7")),
0605 edm::FileInPath(conf.getParameter<string>("Weights8")),
0606 edm::FileInPath(conf.getParameter<string>("Weights9"))}};
0607
0608 for (UInt_t j = 0; j < gbr.size(); ++j) {
0609 gbr[j] = createGBRForest(weights[j]);
0610 }
0611 }
0612 }
0613 }
0614
0615
0616 void GoodSeedProducer::beginRun(const edm::Run& run, const EventSetup& es) {
0617
0618 auto const& magneticField = &es.getData(magneticFieldTokenBeginRun_);
0619 B_ = magneticField->inTesla(GlobalPoint(0, 0, 0));
0620
0621 pfTransformer_ = std::make_unique<PFTrackTransformer>(B_);
0622 pfTransformer_->OnlyProp();
0623
0624
0625 FileInPath ecalEtaMap(conf_.getParameter<string>("EtaMap"));
0626 FileInPath ecalPhiMap(conf_.getParameter<string>("PhiMap"));
0627 resMapEtaECAL_ = std::make_unique<PFResolutionMap>("ECAL_eta", ecalEtaMap.fullPath().c_str());
0628 resMapPhiECAL_ = std::make_unique<PFResolutionMap>("ECAL_phi", ecalPhiMap.fullPath().c_str());
0629
0630
0631 FileInPath parFile(conf_.getParameter<string>("ThresholdFile"));
0632 ifstream ifs(parFile.fullPath().c_str());
0633 for (int iy = 0; iy < 81; ++iy)
0634 ifs >> thr[iy];
0635 }
0636
0637 int GoodSeedProducer::getBin(float eta, float pt) {
0638 int ie = 0;
0639 int ip = 0;
0640 if (fabs(eta) < 0.8)
0641 ie = 0;
0642 else {
0643 if (fabs(eta) < 1.479)
0644 ie = 1;
0645 else
0646 ie = 2;
0647 }
0648 if (pt < 6)
0649 ip = 0;
0650 else {
0651 if (pt < 12)
0652 ip = 1;
0653 else
0654 ip = 2;
0655 }
0656 int iep = ie * 3 + ip;
0657 LogDebug("GoodSeedProducer") << "Track pt =" << pt << " eta=" << eta << " bin=" << iep;
0658 return iep;
0659 }
0660
0661 void GoodSeedProducer::fillPreIdRefValueMap(Handle<TrackCollection> tracks,
0662 const edm::OrphanHandle<reco::PreIdCollection>& preidhandle,
0663 edm::ValueMap<reco::PreIdRef>::Filler& filler) {
0664 std::vector<reco::PreIdRef> values;
0665
0666 unsigned ntracks = tracks->size();
0667 for (unsigned itrack = 0; itrack < ntracks; ++itrack) {
0668 reco::TrackRef theTrackRef(tracks, itrack);
0669 std::map<reco::TrackRef, unsigned>::const_iterator itcheck = refMap_.find(theTrackRef);
0670 if (itcheck == refMap_.end()) {
0671
0672 values.push_back(reco::PreIdRef());
0673 } else {
0674 edm::Ref<reco::PreIdCollection> preIdRef(preidhandle, itcheck->second);
0675 values.push_back(preIdRef);
0676
0677 }
0678 }
0679 filler.insert(tracks, values.begin(), values.end());
0680 }
0681
0682 DEFINE_FWK_MODULE(GoodSeedProducer);