File indexing completed on 2021-02-14 14:26:41
0001 #include "RecoParticleFlow/PFTracking/interface/ConvBremPFTrackFinder.h"
0002 #include "FWCore/Framework/interface/ESHandle.h"
0003
0004 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0005 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0006 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0007 #include "DataFormats/EgammaReco/interface/ElectronSeedFwd.h"
0008 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0009 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0010 #include "TrackingTools/IPTools/interface/IPTools.h"
0011 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0012 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
0013 #include "TMath.h"
0014 #include "RecoParticleFlow/PFClusterTools/interface/LinkByRecHit.h"
0015 #include "TMVA/MethodBDT.h"
0016
0017 using namespace edm;
0018 using namespace std;
0019 using namespace reco;
0020
0021 ConvBremPFTrackFinder::ConvBremPFTrackFinder(const TransientTrackBuilder& builder,
0022 double mvaBremConvCutBarrelLowPt,
0023 double mvaBremConvCutBarrelHighPt,
0024 double mvaBremConvCutEndcapsLowPt,
0025 double mvaBremConvCutEndcapsHighPt)
0026 : builder_(builder),
0027 mvaBremConvCutBarrelLowPt_(mvaBremConvCutBarrelLowPt),
0028 mvaBremConvCutBarrelHighPt_(mvaBremConvCutBarrelHighPt),
0029 mvaBremConvCutEndcapsLowPt_(mvaBremConvCutEndcapsLowPt),
0030 mvaBremConvCutEndcapsHighPt_(mvaBremConvCutEndcapsHighPt) {}
0031 ConvBremPFTrackFinder::~ConvBremPFTrackFinder() {}
0032
0033 void ConvBremPFTrackFinder::runConvBremFinder(const Handle<PFRecTrackCollection>& thePfRecTrackCol,
0034 const Handle<VertexCollection>& primaryVertex,
0035 const edm::Handle<reco::PFDisplacedTrackerVertexCollection>& pfNuclears,
0036 const edm::Handle<reco::PFConversionCollection>& pfConversions,
0037 const edm::Handle<reco::PFV0Collection>& pfV0,
0038 const convbremhelpers::HeavyObjectCache* cache,
0039 bool useNuclear,
0040 bool useConversions,
0041 bool useV0,
0042 const reco::PFClusterCollection& theEClus,
0043 const reco::GsfPFRecTrack& gsfpfrectk) {
0044 found_ = false;
0045 bool debug = false;
0046 bool debugRef = false;
0047
0048 if (debug)
0049 cout << "runConvBremFinder:: Entering " << endl;
0050
0051 const reco::GsfTrackRef& refGsf = gsfpfrectk.gsfTrackRef();
0052 const reco::PFRecTrackRef& pfTrackRef = gsfpfrectk.kfPFRecTrackRef();
0053 vector<PFBrem> primPFBrem = gsfpfrectk.PFRecBrem();
0054
0055 const PFRecTrackCollection& PfRTkColl = *(thePfRecTrackCol.product());
0056 reco::PFRecTrackCollection::const_iterator pft = PfRTkColl.begin();
0057 reco::PFRecTrackCollection::const_iterator pftend = PfRTkColl.end();
0058
0059 vector<PFRecTrackRef> AllPFRecTracks;
0060 AllPFRecTracks.clear();
0061 unsigned int ipft = 0;
0062
0063 for (; pft != pftend; ++pft, ipft++) {
0064
0065 if (pfTrackRef.isNonnull())
0066 if (pfTrackRef->trackRef() == pft->trackRef())
0067 continue;
0068
0069 PFRecTrackRef pfRecTrRef(thePfRecTrackCol, ipft);
0070 TrackRef trackRef = pfRecTrRef->trackRef();
0071 reco::TrackBaseRef selTrackBaseRef(trackRef);
0072
0073 if (debug)
0074 cout << "runConvBremFinder:: pushing_back High Purity " << pft->trackRef()->pt() << " eta,phi "
0075 << pft->trackRef()->eta() << ", " << pft->trackRef()->phi() << " Memory Address Ref " << &*trackRef
0076 << " Memory Address BaseRef " << &*selTrackBaseRef << endl;
0077 AllPFRecTracks.push_back(pfRecTrRef);
0078 }
0079
0080 if (useConversions) {
0081 const PFConversionCollection& PfConvColl = *(pfConversions.product());
0082 for (unsigned i = 0; i < PfConvColl.size(); i++) {
0083 reco::PFConversionRef convRef(pfConversions, i);
0084
0085 unsigned int trackSize = (convRef->pfTracks()).size();
0086 if (convRef->pfTracks().size() < 2)
0087 continue;
0088 for (unsigned iTk = 0; iTk < trackSize; iTk++) {
0089 PFRecTrackRef compPFTkRef = convRef->pfTracks()[iTk];
0090 reco::TrackBaseRef newTrackBaseRef(compPFTkRef->trackRef());
0091
0092 if (pfTrackRef.isNonnull()) {
0093 reco::TrackBaseRef primaryTrackBaseRef(pfTrackRef->trackRef());
0094 if (primaryTrackBaseRef == newTrackBaseRef)
0095 continue;
0096 }
0097 bool notFound = true;
0098 for (unsigned iPF = 0; iPF < AllPFRecTracks.size(); iPF++) {
0099 reco::TrackBaseRef selTrackBaseRef(AllPFRecTracks[iPF]->trackRef());
0100
0101 if (debugRef)
0102 cout << "## Track 1 HP pt " << AllPFRecTracks[iPF]->trackRef()->pt() << " eta, phi "
0103 << AllPFRecTracks[iPF]->trackRef()->eta() << ", " << AllPFRecTracks[iPF]->trackRef()->phi()
0104 << " Memory Address Ref " << &(*AllPFRecTracks[iPF]->trackRef()) << " Memory Address BaseRef "
0105 << &*selTrackBaseRef << endl;
0106 if (debugRef)
0107 cout << "** Track 2 CONV pt " << compPFTkRef->trackRef()->pt() << " eta, phi "
0108 << compPFTkRef->trackRef()->eta() << ", " << compPFTkRef->trackRef()->phi() << " Memory Address Ref "
0109 << &*compPFTkRef->trackRef() << " Memory Address BaseRef " << &*newTrackBaseRef << endl;
0110
0111 if (AllPFRecTracks[iPF]->trackRef() == compPFTkRef->trackRef()) {
0112 if (debugRef)
0113 cout << " SAME BREM REF " << endl;
0114 notFound = false;
0115 }
0116 }
0117 if (notFound) {
0118 if (debug)
0119 cout << "runConvBremFinder:: pushing_back Conversions " << compPFTkRef->trackRef()->pt() << " eta,phi "
0120 << compPFTkRef->trackRef()->eta() << " phi " << compPFTkRef->trackRef()->phi() << endl;
0121 AllPFRecTracks.push_back(compPFTkRef);
0122 }
0123 }
0124 }
0125 }
0126
0127 if (useNuclear) {
0128 const PFDisplacedTrackerVertexCollection& PfNuclColl = *(pfNuclears.product());
0129 for (unsigned i = 0; i < PfNuclColl.size(); i++) {
0130 const reco::PFDisplacedTrackerVertexRef dispacedVertexRef(pfNuclears, i);
0131 unsigned int trackSize = dispacedVertexRef->pfRecTracks().size();
0132 for (unsigned iTk = 0; iTk < trackSize; iTk++) {
0133 reco::PFRecTrackRef newPFRecTrackRef = dispacedVertexRef->pfRecTracks()[iTk];
0134 reco::TrackBaseRef newTrackBaseRef(newPFRecTrackRef->trackRef());
0135
0136 if (pfTrackRef.isNonnull()) {
0137 reco::TrackBaseRef primaryTrackBaseRef(pfTrackRef->trackRef());
0138 if (primaryTrackBaseRef == newTrackBaseRef)
0139 continue;
0140 }
0141 bool notFound = true;
0142 for (unsigned iPF = 0; iPF < AllPFRecTracks.size(); iPF++) {
0143 reco::TrackBaseRef selTrackBaseRef(AllPFRecTracks[iPF]->trackRef());
0144 if (selTrackBaseRef == newTrackBaseRef)
0145 notFound = false;
0146 }
0147 if (notFound) {
0148 if (debug)
0149 cout << "runConvBremFinder:: pushing_back displaced Vertex pt " << newPFRecTrackRef->trackRef()->pt()
0150 << " eta,phi " << newPFRecTrackRef->trackRef()->eta() << ", " << newPFRecTrackRef->trackRef()->phi()
0151 << endl;
0152 AllPFRecTracks.push_back(newPFRecTrackRef);
0153 }
0154 }
0155 }
0156 }
0157
0158 if (useV0) {
0159 const PFV0Collection& PfV0Coll = *(pfV0.product());
0160 for (unsigned i = 0; i < PfV0Coll.size(); i++) {
0161 reco::PFV0Ref v0Ref(pfV0, i);
0162 unsigned int trackSize = (v0Ref->pfTracks()).size();
0163 for (unsigned iTk = 0; iTk < trackSize; iTk++) {
0164 reco::PFRecTrackRef newPFRecTrackRef = (v0Ref->pfTracks())[iTk];
0165 reco::TrackBaseRef newTrackBaseRef(newPFRecTrackRef->trackRef());
0166
0167 if (pfTrackRef.isNonnull()) {
0168 reco::TrackBaseRef primaryTrackBaseRef(pfTrackRef->trackRef());
0169 if (primaryTrackBaseRef == newTrackBaseRef)
0170 continue;
0171 }
0172 bool notFound = true;
0173 for (unsigned iPF = 0; iPF < AllPFRecTracks.size(); iPF++) {
0174 reco::TrackBaseRef selTrackBaseRef(AllPFRecTracks[iPF]->trackRef());
0175 if (selTrackBaseRef == newTrackBaseRef)
0176 notFound = false;
0177 }
0178 if (notFound) {
0179 if (debug)
0180 cout << "runConvBremFinder:: pushing_back V0 " << newPFRecTrackRef->trackRef()->pt() << " eta,phi "
0181 << newPFRecTrackRef->trackRef()->eta() << ", " << newPFRecTrackRef->trackRef()->phi() << endl;
0182 AllPFRecTracks.push_back(newPFRecTrackRef);
0183 }
0184 }
0185 }
0186 }
0187
0188 pfRecTrRef_vec_.clear();
0189
0190 for (unsigned iPF = 0; iPF < AllPFRecTracks.size(); iPF++) {
0191 double dphi = fabs(AllPFRecTracks[iPF]->trackRef()->phi() - refGsf->phi());
0192 if (dphi > TMath::Pi())
0193 dphi -= TMath::TwoPi();
0194 double deta = fabs(AllPFRecTracks[iPF]->trackRef()->eta() - refGsf->eta());
0195
0196
0197 if (fabs(dphi) > 1.0 || fabs(deta) > 0.4)
0198 continue;
0199
0200 double minDEtaBremKF = 1000.;
0201 double minDPhiBremKF = 1000.;
0202 double minDRBremKF = 1000.;
0203 double minDEtaBremKFPos = 1000.;
0204 double minDPhiBremKFPos = 1000.;
0205 double minDRBremKFPos = 1000.;
0206 reco::TrackRef trkRef = AllPFRecTracks[iPF]->trackRef();
0207
0208 double secEta = trkRef->innerMomentum().eta();
0209 double secPhi = trkRef->innerMomentum().phi();
0210
0211 for (unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
0212 if (primPFBrem[ipbrem].indTrajPoint() == 99)
0213 continue;
0214 const reco::PFTrajectoryPoint& atPrimECAL =
0215 primPFBrem[ipbrem].extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance);
0216 if (!atPrimECAL.isValid())
0217 continue;
0218 double bremEta = atPrimECAL.momentum().Eta();
0219 double bremPhi = atPrimECAL.momentum().Phi();
0220
0221 double deta = fabs(bremEta - secEta);
0222 double dphi = fabs(bremPhi - secPhi);
0223 if (dphi > TMath::Pi())
0224 dphi -= TMath::TwoPi();
0225 double DR = sqrt(deta * deta + dphi * dphi);
0226
0227 double detaPos = fabs(bremEta - trkRef->innerPosition().eta());
0228 double dphiPos = fabs(bremPhi - trkRef->innerPosition().phi());
0229 if (dphiPos > TMath::Pi())
0230 dphiPos -= TMath::TwoPi();
0231 double DRPos = sqrt(detaPos * detaPos + dphiPos * dphiPos);
0232
0233
0234 if (DR < minDRBremKF) {
0235 minDRBremKF = DR;
0236 minDEtaBremKF = deta;
0237 minDPhiBremKF = fabs(dphi);
0238 }
0239
0240 if (DRPos < minDRBremKFPos) {
0241 minDRBremKFPos = DR;
0242 minDEtaBremKFPos = detaPos;
0243 minDPhiBremKFPos = fabs(dphiPos);
0244 }
0245 }
0246
0247
0248 float gsfR = sqrt(refGsf->innerPosition().x() * refGsf->innerPosition().x() +
0249 refGsf->innerPosition().y() * refGsf->innerPosition().y());
0250
0251
0252 secR = sqrt(trkRef->innerPosition().x() * trkRef->innerPosition().x() +
0253 trkRef->innerPosition().y() * trkRef->innerPosition().y());
0254
0255
0256
0257 if ((minDPhiBremKF < 0.1 || minDPhiBremKFPos < 0.1) && (minDEtaBremKF < 0.02 || minDEtaBremKFPos < 0.02) &&
0258 secR > (gsfR - 8)) {
0259 if (debug)
0260 cout << "runConvBremFinder:: OK Find track and BREM close "
0261 << " MinDphi " << minDPhiBremKF << " MinDeta " << minDEtaBremKF << endl;
0262
0263 float MinDist = 100000.;
0264 float EE_calib = 0.;
0265 PFRecTrack pfrectrack = *AllPFRecTracks[iPF];
0266
0267 for (PFClusterCollection::const_iterator clus = theEClus.begin(); clus != theEClus.end(); clus++) {
0268
0269 clus->position();
0270 double dist = -1.;
0271 PFCluster clust = *clus;
0272 clust.calculatePositionREP();
0273 dist = pfrectrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax).isValid()
0274 ? LinkByRecHit::testTrackAndClusterByRecHit(pfrectrack, clust)
0275 : -1.;
0276
0277 if (dist > 0. && dist < MinDist) {
0278 MinDist = dist;
0279 EE_calib = cache->pfcalib_->energyEm(*clus, 0.0, 0.0, false);
0280 }
0281 }
0282 if (MinDist > 0. && MinDist < 100000.) {
0283
0284
0285 secPout = sqrt(trkRef->outerMomentum().x() * trkRef->outerMomentum().x() +
0286 trkRef->outerMomentum().y() * trkRef->outerMomentum().y() +
0287 trkRef->outerMomentum().z() * trkRef->outerMomentum().z());
0288
0289 secPin = sqrt(trkRef->innerMomentum().x() * trkRef->innerMomentum().x() +
0290 trkRef->innerMomentum().y() * trkRef->innerMomentum().y() +
0291 trkRef->innerMomentum().z() * trkRef->innerMomentum().z());
0292
0293
0294 ptRatioGsfKF = trkRef->pt() / (refGsf->ptMode());
0295
0296 Vertex dummy;
0297 const Vertex* pv = &dummy;
0298 edm::Ref<VertexCollection> pvRef;
0299 if (!primaryVertex->empty()) {
0300 pv = &*primaryVertex->begin();
0301
0302 pvRef = edm::Ref<VertexCollection>(primaryVertex, 0);
0303 } else {
0304 Vertex::Error e;
0305 e(0, 0) = 0.0015 * 0.0015;
0306 e(1, 1) = 0.0015 * 0.0015;
0307 e(2, 2) = 15. * 15.;
0308 Vertex::Point p(0, 0, 0);
0309 dummy = Vertex(p, e, 0, 0, 0);
0310 }
0311
0312
0313 GlobalVector direction(refGsf->innerMomentum().x(), refGsf->innerMomentum().y(), refGsf->innerMomentum().z());
0314
0315 TransientTrack transientTrack = builder_.build(*trkRef);
0316 sTIP = IPTools::signedTransverseImpactParameter(transientTrack, direction, *pv).second.significance();
0317
0318 Epout = EE_calib / secPout;
0319
0320
0321 detaBremKF = minDEtaBremKF;
0322
0323
0324 unsigned int tmpSh = 0;
0325
0326 int kfhitcounter = 0;
0327 for (auto const& nhit : refGsf->recHits())
0328 if (nhit->isValid()) {
0329 kfhitcounter = 0;
0330 for (auto const& ihit : trkRef->recHits())
0331 if (ihit->isValid()) {
0332
0333 if (nhit->sharesInput(ihit, TrackingRecHit::all))
0334 tmpSh++;
0335 kfhitcounter++;
0336
0337
0338
0339 }
0340 }
0341
0342 nHITS1 = tmpSh;
0343
0344 TString weightfilepath = "";
0345 double mvaValue = -10;
0346 double cutvalue = -10;
0347
0348 float vars[6] = {secR, sTIP, nHITS1, Epout, detaBremKF, ptRatioGsfKF};
0349 if (refGsf->pt() < 20 && fabs(refGsf->eta()) < 1.479) {
0350 mvaValue = cache->gbrBarrelLowPt_->GetClassifier(vars);
0351 cutvalue = mvaBremConvCutBarrelLowPt_;
0352 }
0353 if (refGsf->pt() > 20 && fabs(refGsf->eta()) < 1.479) {
0354 mvaValue = cache->gbrBarrelHighPt_->GetClassifier(vars);
0355 cutvalue = mvaBremConvCutBarrelHighPt_;
0356 }
0357 if (refGsf->pt() < 20 && fabs(refGsf->eta()) > 1.479) {
0358 mvaValue = cache->gbrEndcapsLowPt_->GetClassifier(vars);
0359 cutvalue = mvaBremConvCutEndcapsLowPt_;
0360 }
0361 if (refGsf->pt() > 20 && fabs(refGsf->eta()) > 1.479) {
0362 mvaValue = cache->gbrEndcapsHighPt_->GetClassifier(vars);
0363 cutvalue = mvaBremConvCutEndcapsHighPt_;
0364 }
0365
0366 if (debug)
0367 cout << "Gsf track Pt, Eta " << refGsf->pt() << " " << refGsf->eta() << endl;
0368 if (debug)
0369 cout << "Cutvalue " << cutvalue << endl;
0370
0371 if ((kfhitcounter - nHITS1) <= 3 && nHITS1 > 3)
0372 mvaValue = 2;
0373
0374 if (debug)
0375 cout << " The input variables for conv brem tracks identification " << endl
0376 << " secR " << secR << " gsfR " << gsfR << endl
0377 << " N shared hits " << nHITS1 << endl
0378 << " sTIP " << sTIP << endl
0379 << " detaBremKF " << detaBremKF << endl
0380 << " E/pout " << Epout << endl
0381 << " ptRatioKFGsf " << ptRatioGsfKF << endl
0382 << " ***** MVA ***** " << mvaValue << endl;
0383
0384 if (mvaValue > cutvalue) {
0385 found_ = true;
0386 pfRecTrRef_vec_.push_back(AllPFRecTracks[iPF]);
0387 }
0388 }
0389 }
0390 }
0391 }