Back to home page

Project CMSSW displayed by LXR
 
 

    

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     // do not consider the kf track already associated to the seed
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         // do not consider the kf track already associated to the seed
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           //if(selTrackBaseRef == newTrackBaseRef ||  AllPFRecTracks[iPF]->trackRef()== compPFTkRef->trackRef()) {
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         // do not consider the kf track already associated to the seed
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         // do not consider the kf track already associated to the seed
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     // limiting the phase space (just for saving cpu-time)
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       // find the closest track tangent
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     //gsfR
0248     float gsfR = sqrt(refGsf->innerPosition().x() * refGsf->innerPosition().x() +
0249                       refGsf->innerPosition().y() * refGsf->innerPosition().y());
0250 
0251     // secR
0252     secR = sqrt(trkRef->innerPosition().x() * trkRef->innerPosition().x() +
0253                 trkRef->innerPosition().y() * trkRef->innerPosition().y());
0254 
0255     // apply loose selection (to be parallel) between the secondary track and brem-tangents.
0256     // Moreover if the secR is internal with respect to the GSF track by two pixel layers discard it.
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       // Find and ECAL associated cluster
0267       for (PFClusterCollection::const_iterator clus = theEClus.begin(); clus != theEClus.end(); clus++) {
0268         // Removed unusd variable, left this in case it has side effects
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         // compute all the input variables for conv brem selection
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         // maybe put innter momentum pt?
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           // we always use the first vertex (at the moment)
0302           pvRef = edm::Ref<VertexCollection>(primaryVertex, 0);
0303         } else {  // create a dummy PV
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         // direction of the Gsf track
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         // eta distance brem-secondary kf track
0321         detaBremKF = minDEtaBremKF;
0322 
0323         // Number of commont hits
0324         unsigned int tmpSh = 0;
0325         //uint ish=0;
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                 // method 1
0333                 if (nhit->sharesInput(ihit, TrackingRecHit::all))
0334                   tmpSh++;
0335                 kfhitcounter++;
0336                 // method 2 to switch in case of problem with rechit collections
0337                 // if((ihit->geographicalId()==nhit->geographicalId())&&
0338                 //     ((nhit->localPosition()-ihit->localPosition()).mag()<0.01)) ish++;
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;  // this means duplicates tracks, very likely not physical
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       }  // end MinDIST
0389     }    // end selection kf - brem tangents
0390   }      // loop on the kf tracks
0391 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!