PFTracking/src/ConvBremPFTrackFinder.cc

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 #include "DataFormats/Math/interface/deltaPhi.h"
0017
0018 using namespace edm;
0019 using namespace std;
0020 using namespace reco;
0021
0022 ConvBremPFTrackFinder::ConvBremPFTrackFinder(const TransientTrackBuilder& builder,
0023                                              double mvaBremConvCutBarrelLowPt,
0024                                              double mvaBremConvCutBarrelHighPt,
0025                                              double mvaBremConvCutEndcapsLowPt,
0026                                              double mvaBremConvCutEndcapsHighPt)
0027     : builder_(builder),
0028       mvaBremConvCutBarrelLowPt_(mvaBremConvCutBarrelLowPt),
0029       mvaBremConvCutBarrelHighPt_(mvaBremConvCutBarrelHighPt),
0030       mvaBremConvCutEndcapsLowPt_(mvaBremConvCutEndcapsLowPt),
0031       mvaBremConvCutEndcapsHighPt_(mvaBremConvCutEndcapsHighPt) {}
0032 ConvBremPFTrackFinder::~ConvBremPFTrackFinder() {}
0033
0034 void ConvBremPFTrackFinder::runConvBremFinder(const Handle<PFRecTrackCollection>& thePfRecTrackCol,
0035                                               const Handle<VertexCollection>& primaryVertex,
0036                                               const edm::Handle<reco::PFDisplacedTrackerVertexCollection>& pfNuclears,
0037                                               const edm::Handle<reco::PFConversionCollection>& pfConversions,
0038                                               const edm::Handle<reco::PFV0Collection>& pfV0,
0039                                               const convbremhelpers::HeavyObjectCache* cache,
0040                                               bool useNuclear,
0041                                               bool useConversions,
0042                                               bool useV0,
0043                                               const reco::PFClusterCollection& theEClus,
0044                                               const reco::GsfPFRecTrack& gsfpfrectk) {
0045   found_ = false;
0046
0047   LogDebug("ConvBremPFTrackFinder") << "runConvBremFinder:: Entering ";
0048
0049   const reco::GsfTrackRef& refGsf = gsfpfrectk.gsfTrackRef();
0050   float refGsfEta = refGsf->eta();
0051   float refGsfPt = refGsf->pt();
0052   float refGsfPhi = refGsf->phi();
0053   float gsfR = sqrt(refGsf->innerPosition().x() * refGsf->innerPosition().x() +
0054                     refGsf->innerPosition().y() * refGsf->innerPosition().y());
0055   // direction of the Gsf track
0056   GlobalVector direction(refGsf->innerMomentum().x(), refGsf->innerMomentum().y(), refGsf->innerMomentum().z());
0057   float refGsfPtMode = refGsf->ptMode();
0058
0059   const reco::PFRecTrackRef& pfTrackRef = gsfpfrectk.kfPFRecTrackRef();
0060   vector<PFBrem> primPFBrem = gsfpfrectk.PFRecBrem();
0061
0062   const PFRecTrackCollection& PfRTkColl = *(thePfRecTrackCol.product());
0063   reco::PFRecTrackCollection::const_iterator pft = PfRTkColl.begin();
0064   reco::PFRecTrackCollection::const_iterator pftend = PfRTkColl.end();
0065
0066   vector<PFRecTrackRef> AllPFRecTracks;
0067   AllPFRecTracks.clear();
0068   unsigned int ipft = 0;
0069
0070   for (; pft != pftend; ++pft, ipft++) {
0071     // do not consider the kf track already associated to the seed
0072     if (pfTrackRef.isNonnull())
0073       if (pfTrackRef->trackRef() == pft->trackRef())
0074         continue;
0075
0076     PFRecTrackRef pfRecTrRef(thePfRecTrackCol, ipft);
0077     TrackRef trackRef = pfRecTrRef->trackRef();
0078     reco::TrackBaseRef selTrackBaseRef(trackRef);
0079
0080     LogDebug("ConvBremPFTrackFinder") << "runConvBremFinder:: pushing_back High Purity " << pft->trackRef()->pt()
0081                                       << " eta,phi " << pft->trackRef()->eta() << ", " << pft->trackRef()->phi()
0082                                       << " Memory Address Ref  " << &*trackRef << " Memory Address BaseRef  "
0083                                       << &*selTrackBaseRef;
0084     AllPFRecTracks.push_back(pfRecTrRef);
0085   }
0086
0087   if (useConversions) {
0088     const PFConversionCollection& PfConvColl = *(pfConversions.product());
0089     for (unsigned i = 0; i < PfConvColl.size(); i++) {
0090       reco::PFConversionRef convRef(pfConversions, i);
0091
0092       unsigned int trackSize = (convRef->pfTracks()).size();
0093       if (convRef->pfTracks().size() < 2)
0094         continue;
0095       for (unsigned iTk = 0; iTk < trackSize; iTk++) {
0096         PFRecTrackRef compPFTkRef = convRef->pfTracks()[iTk];
0097         reco::TrackBaseRef newTrackBaseRef(compPFTkRef->trackRef());
0098         // do not consider the kf track already associated to the seed
0099         if (pfTrackRef.isNonnull()) {
0100           reco::TrackBaseRef primaryTrackBaseRef(pfTrackRef->trackRef());
0101           if (primaryTrackBaseRef == newTrackBaseRef)
0102             continue;
0103         }
0104         bool notFound = true;
0105         for (unsigned iPF = 0; iPF < AllPFRecTracks.size(); iPF++) {
0106           reco::TrackBaseRef selTrackBaseRef(AllPFRecTracks[iPF]->trackRef());
0107
0108           LogDebug("ConvBremPFTrackFinder")
0109               << "## Track 1 HP pt " << AllPFRecTracks[iPF]->trackRef()->pt() << " eta, phi "
0110               << AllPFRecTracks[iPF]->trackRef()->eta() << ", " << AllPFRecTracks[iPF]->trackRef()->phi()
0111               << " Memory Address Ref  " << &(*AllPFRecTracks[iPF]->trackRef()) << " Memory Address BaseRef  "
0112               << &*selTrackBaseRef;
0113           LogDebug("ConvBremPFTrackFinder")
0114               << "** Track 2 CONV pt " << compPFTkRef->trackRef()->pt() << " eta, phi "
0115               << compPFTkRef->trackRef()->eta() << ", " << compPFTkRef->trackRef()->phi() << " Memory Address Ref "
0116               << &*compPFTkRef->trackRef() << " Memory Address BaseRef " << &*newTrackBaseRef;
0117           //if(selTrackBaseRef == newTrackBaseRef ||  AllPFRecTracks[iPF]->trackRef()== compPFTkRef->trackRef()) {
0118           if (AllPFRecTracks[iPF]->trackRef() == compPFTkRef->trackRef()) {
0119             LogDebug("ConvBremPFTrackFinder") << "  SAME BREM REF " << endl;
0120             notFound = false;
0121           }
0122         }
0123         if (notFound) {
0124           LogDebug("ConvBremPFTrackFinder")
0125               << "runConvBremFinder:: pushing_back Conversions " << compPFTkRef->trackRef()->pt() << " eta,phi "
0126               << compPFTkRef->trackRef()->eta() << " phi " << compPFTkRef->trackRef()->phi();
0127           AllPFRecTracks.push_back(compPFTkRef);
0128         }
0129       }
0130     }
0131   }
0132
0133   if (useNuclear) {
0134     const PFDisplacedTrackerVertexCollection& PfNuclColl = *(pfNuclears.product());
0135     for (unsigned i = 0; i < PfNuclColl.size(); i++) {
0136       const reco::PFDisplacedTrackerVertexRef dispacedVertexRef(pfNuclears, i);
0137       unsigned int trackSize = dispacedVertexRef->pfRecTracks().size();
0138       for (unsigned iTk = 0; iTk < trackSize; iTk++) {
0139         reco::PFRecTrackRef newPFRecTrackRef = dispacedVertexRef->pfRecTracks()[iTk];
0140         reco::TrackBaseRef newTrackBaseRef(newPFRecTrackRef->trackRef());
0141         // do not consider the kf track already associated to the seed
0142         if (pfTrackRef.isNonnull()) {
0143           reco::TrackBaseRef primaryTrackBaseRef(pfTrackRef->trackRef());
0144           if (primaryTrackBaseRef == newTrackBaseRef)
0145             continue;
0146         }
0147         bool notFound = true;
0148         for (unsigned iPF = 0; iPF < AllPFRecTracks.size(); iPF++) {
0149           reco::TrackBaseRef selTrackBaseRef(AllPFRecTracks[iPF]->trackRef());
0150           if (selTrackBaseRef == newTrackBaseRef)
0151             notFound = false;
0152         }
0153         if (notFound) {
0154           LogDebug("ConvBremPFTrackFinder")
0155               << "runConvBremFinder:: pushing_back displaced Vertex pt " << newPFRecTrackRef->trackRef()->pt()
0156               << " eta,phi " << newPFRecTrackRef->trackRef()->eta() << ", " << newPFRecTrackRef->trackRef()->phi();
0157           AllPFRecTracks.push_back(newPFRecTrackRef);
0158         }
0159       }
0160     }
0161   }
0162
0163   if (useV0) {
0164     const PFV0Collection& PfV0Coll = *(pfV0.product());
0165     for (unsigned i = 0; i < PfV0Coll.size(); i++) {
0166       reco::PFV0Ref v0Ref(pfV0, i);
0167       unsigned int trackSize = (v0Ref->pfTracks()).size();
0168       for (unsigned iTk = 0; iTk < trackSize; iTk++) {
0169         reco::PFRecTrackRef newPFRecTrackRef = (v0Ref->pfTracks())[iTk];
0170         reco::TrackBaseRef newTrackBaseRef(newPFRecTrackRef->trackRef());
0171         // do not consider the kf track already associated to the seed
0172         if (pfTrackRef.isNonnull()) {
0173           reco::TrackBaseRef primaryTrackBaseRef(pfTrackRef->trackRef());
0174           if (primaryTrackBaseRef == newTrackBaseRef)
0175             continue;
0176         }
0177         bool notFound = true;
0178         for (unsigned iPF = 0; iPF < AllPFRecTracks.size(); iPF++) {
0179           reco::TrackBaseRef selTrackBaseRef(AllPFRecTracks[iPF]->trackRef());
0180           if (selTrackBaseRef == newTrackBaseRef)
0181             notFound = false;
0182         }
0183         if (notFound) {
0184           LogDebug("ConvBremPFTrackFinder")
0185               << "runConvBremFinder:: pushing_back V0 " << newPFRecTrackRef->trackRef()->pt() << " eta,phi "
0186               << newPFRecTrackRef->trackRef()->eta() << ", " << newPFRecTrackRef->trackRef()->phi();
0187           AllPFRecTracks.push_back(newPFRecTrackRef);
0188         }
0189       }
0190     }
0191   }
0192
0193   pfRecTrRef_vec_.clear();
0194
0195   for (unsigned iPF = 0; iPF < AllPFRecTracks.size(); iPF++) {
0196     double dphi = std::abs(::deltaPhi(AllPFRecTracks[iPF]->trackRef()->phi(), refGsfPhi));
0197     double deta = std::abs(AllPFRecTracks[iPF]->trackRef()->eta() - refGsfEta);
0198
0199     // limiting the phase space (just for saving cpu-time)
0200     if (std::abs(dphi) > 1.0 || std::abs(deta) > 0.4)
0201       continue;
0202
0203     double minDEtaBremKF = std::numeric_limits<double>::max();
0204     double minDPhiBremKF = std::numeric_limits<double>::max();
0205     double minDRBremKF2 = std::numeric_limits<double>::max();
0206     double minDEtaBremKFPos = std::numeric_limits<double>::max();
0207     double minDPhiBremKFPos = std::numeric_limits<double>::max();
0208     double minDRBremKFPos2 = std::numeric_limits<double>::max();
0209     reco::TrackRef trkRef = AllPFRecTracks[iPF]->trackRef();
0210
0211     double secEta = trkRef->innerMomentum().eta();
0212     double secPhi = trkRef->innerMomentum().phi();
0213
0214     double posEta = trkRef->innerPosition().eta();
0215     double posPhi = trkRef->innerPosition().phi();
0216
0217     for (unsigned ipbrem = 0; ipbrem < primPFBrem.size(); ipbrem++) {
0218       if (primPFBrem[ipbrem].indTrajPoint() == 99)
0219         continue;
0220       const reco::PFTrajectoryPoint& atPrimECAL =
0221           primPFBrem[ipbrem].extrapolatedPoint(reco::PFTrajectoryPoint::ECALEntrance);
0222       if (!atPrimECAL.isValid())
0223         continue;
0224       double bremEta = atPrimECAL.momentum().Eta();
0225       double bremPhi = atPrimECAL.momentum().Phi();
0226
0227       double deta = std::abs(bremEta - secEta);
0228       double dphi = std::abs(::deltaPhi(bremPhi, secPhi));
0229       double DR2 = deta * deta + dphi * dphi;
0230
0231       double detaPos = std::abs(bremEta - posEta);
0232       double dphiPos = std::abs(::deltaPhi(bremPhi, posPhi));
0233       double DRPos2 = detaPos * detaPos + dphiPos * dphiPos;
0234
0235       // find the closest track tangent
0236       if (DR2 < minDRBremKF2) {
0237         minDRBremKF2 = DR2;
0238         minDEtaBremKF = deta;
0239         minDPhiBremKF = std::abs(dphi);
0240       }
0241
0242       if (DRPos2 < minDRBremKFPos2) {
0243         minDRBremKFPos2 = DRPos2;
0244         minDEtaBremKFPos = detaPos;
0245         minDPhiBremKFPos = std::abs(dphiPos);
0246       }
0247     }
0248
0249     // secR
0250     secR = sqrt(trkRef->innerPosition().x() * trkRef->innerPosition().x() +
0251                 trkRef->innerPosition().y() * trkRef->innerPosition().y());
0252
0253     // apply loose selection (to be parallel) between the secondary track and brem-tangents.
0254     // Moreover if the secR is internal with respect to the GSF track by two pixel layers discard it.
0255     if ((minDPhiBremKF < 0.1 || minDPhiBremKFPos < 0.1) && (minDEtaBremKF < 0.02 || minDEtaBremKFPos < 0.02) &&
0256         secR > (gsfR - 8)) {
0257       LogDebug("ConvBremPFTrackFinder") << "runConvBremFinder:: OK Find track and BREM close "
0258                                         << " MinDphi " << minDPhiBremKF << " MinDeta " << minDEtaBremKF;
0259
0260       float MinDist = std::numeric_limits<float>::max();
0261       float EE_calib = 0.;
0262       PFRecTrack pfrectrack = *AllPFRecTracks[iPF];
0263       // Find and ECAL associated cluster
0264       for (PFClusterCollection::const_iterator clus = theEClus.begin(); clus != theEClus.end(); clus++) {
0265         // Removed unusd variable, left this in case it has side effects
0266         clus->position();
0267         double dist = -1.;
0268         PFCluster clust = *clus;
0269         clust.calculatePositionREP();
0270         dist = pfrectrack.extrapolatedPoint(reco::PFTrajectoryPoint::ECALShowerMax).isValid()
0271                    ? LinkByRecHit::testTrackAndClusterByRecHit(pfrectrack, clust)
0272                    : -1.;
0273
0274         if (dist > 0. && dist < MinDist) {
0275           MinDist = dist;
0276           EE_calib = cache->pfcalib_->energyEm(*clus, 0.0, 0.0, false);
0277         }
0278       }
0279       if (MinDist > 0. && MinDist < 100000.) {
0280         // compute all the input variables for conv brem selection
0281
0282         secPout = sqrt(trkRef->outerMomentum().x() * trkRef->outerMomentum().x() +
0283                        trkRef->outerMomentum().y() * trkRef->outerMomentum().y() +
0284                        trkRef->outerMomentum().z() * trkRef->outerMomentum().z());
0285
0286         secPin = sqrt(trkRef->innerMomentum().x() * trkRef->innerMomentum().x() +
0287                       trkRef->innerMomentum().y() * trkRef->innerMomentum().y() +
0288                       trkRef->innerMomentum().z() * trkRef->innerMomentum().z());
0289
0290         // maybe put innter momentum pt?
0291         ptRatioGsfKF = trkRef->pt() / refGsfPtMode;
0292
0293         Vertex dummy;
0294         const Vertex* pv = &dummy;
0295         edm::Ref<VertexCollection> pvRef;
0296         if (!primaryVertex->empty()) {
0297           pv = &*primaryVertex->begin();
0298           // we always use the first vertex (at the moment)
0299           pvRef = edm::Ref<VertexCollection>(primaryVertex, 0);
0300         } else {  // create a dummy PV
0301           Vertex::Error e;
0302           e(0, 0) = 0.0015 * 0.0015;
0303           e(1, 1) = 0.0015 * 0.0015;
0304           e(2, 2) = 15. * 15.;
0305           Vertex::Point p(0, 0, 0);
0306           dummy = Vertex(p, e, 0, 0, 0);
0307         }
0308
0309         TransientTrack transientTrack = builder_.build(*trkRef);
0310         sTIP = IPTools::signedTransverseImpactParameter(transientTrack, direction, *pv).second.significance();
0311
0312         Epout = EE_calib / secPout;
0313
0314         // eta distance brem-secondary kf track
0315         detaBremKF = minDEtaBremKF;
0316
0317         // Number of commont hits
0318         unsigned int tmpSh = 0;
0319         //uint ish=0;
0320         int kfhitcounter = 0;
0321         for (auto const& nhit : refGsf->recHits())
0322           if (nhit->isValid()) {
0323             kfhitcounter = 0;
0324             for (auto const& ihit : trkRef->recHits())
0325               if (ihit->isValid()) {
0326                 // method 1
0327                 if (nhit->sharesInput(ihit, TrackingRecHit::all))
0328                   tmpSh++;
0329                 kfhitcounter++;
0330                 // method 2 to switch in case of problem with rechit collections
0331                 // if((ihit->geographicalId()==nhit->geographicalId())&&
0332                 //     ((nhit->localPosition()-ihit->localPosition()).mag()<0.01)) ish++;
0333               }
0334           }
0335
0336         nHITS1 = tmpSh;
0337
0338         double mvaValue = -10;
0339         double cutvalue = -10;
0340
0341         float vars[6] = {secR, sTIP, nHITS1, Epout, detaBremKF, ptRatioGsfKF};
0342         if (refGsfPt < 20 && std::abs(refGsfEta) < 1.479) {
0343           mvaValue = cache->gbrBarrelLowPt_->GetClassifier(vars);
0344           cutvalue = mvaBremConvCutBarrelLowPt_;
0345         }
0346         if (refGsfPt > 20 && std::abs(refGsfEta) < 1.479) {
0347           mvaValue = cache->gbrBarrelHighPt_->GetClassifier(vars);
0348           cutvalue = mvaBremConvCutBarrelHighPt_;
0349         }
0350         if (refGsfPt < 20 && std::abs(refGsfEta) > 1.479) {
0351           mvaValue = cache->gbrEndcapsLowPt_->GetClassifier(vars);
0352           cutvalue = mvaBremConvCutEndcapsLowPt_;
0353         }
0354         if (refGsfPt > 20 && std::abs(refGsfEta) > 1.479) {
0355           mvaValue = cache->gbrEndcapsHighPt_->GetClassifier(vars);
0356           cutvalue = mvaBremConvCutEndcapsHighPt_;
0357         }
0358
0359         LogDebug("ConvBremPFTrackFinder") << "Gsf track Pt, Eta " << refGsfPt << " " << refGsfEta;
0360         LogDebug("ConvBremPFTrackFinder") << "Cutvalue " << cutvalue;
0361
0362         if ((kfhitcounter - nHITS1) <= 3 && nHITS1 > 3)
0363           mvaValue = 2;  // this means duplicates tracks, very likely not physical
0364
0365         LogDebug("ConvBremPFTrackFinder")
0366             << " The input variables for conv brem tracks identification "
0367             << " secR          " << secR << " gsfR " << gsfR << " N shared hits " << nHITS1 << " sTIP          " << sTIP
0368             << " detaBremKF    " << detaBremKF << " E/pout        " << Epout << " ptRatioKFGsf  " << ptRatioGsfKF
0369             << " ***** MVA ***** " << mvaValue;
0370
0371         if (mvaValue > cutvalue) {
0372           found_ = true;
0373           pfRecTrRef_vec_.push_back(AllPFRecTracks[iPF]);
0374         }
0375       }  // end MinDIST
0376     }    // end selection kf - brem tangents
0377   }      // loop on the kf tracks
0378 }