Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:10:19 

0001 #include <TFile.h>
0002 #include "EgammaAnalysis/ElectronTools/interface/PFIsolationEstimator.h"
0003 #include <cmath>
0004 #include "DataFormats/Math/interface/deltaR.h"
0005 
0006 #ifndef STANDALONE
0007 #include "DataFormats/TrackReco/interface/Track.h"
0008 #include "DataFormats/GsfTrackReco/interface/GsfTrack.h"
0009 #include "DataFormats/GsfTrackReco/interface/GsfTrackFwd.h"
0010 #include "DataFormats/TrackReco/interface/TrackFwd.h"
0011 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0012 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0013 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidate.h"
0014 #include "DataFormats/ParticleFlowCandidate/interface/PFCandidateFwd.h"
0015 #include "DataFormats/Common/interface/RefToPtr.h"
0016 #include "DataFormats/VertexReco/interface/Vertex.h"
0017 #include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h"
0018 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0019 #include "TrackingTools/IPTools/interface/IPTools.h"
0020 
0021 #endif
0022 
0023 //--------------------------------------------------------------------------------------------------
0024 PFIsolationEstimator::PFIsolationEstimator() : fisInitialized(kFALSE) {
0025   // Constructor.
0026 }
0027 
0028 //--------------------------------------------------------------------------------------------------
0029 PFIsolationEstimator::~PFIsolationEstimator() {}
0030 
0031 //--------------------------------------------------------------------------------------------------
0032 void PFIsolationEstimator::initialize(Bool_t bApplyVeto, int iParticleType) {
0033   setParticleType(iParticleType);
0034 
0035   //By default check for an option vertex association
0036   checkClosestZVertex = kTRUE;
0037 
0038   //Apply vetoes
0039   setApplyVeto(bApplyVeto);
0040 
0041   setDeltaRVetoBarrelPhotons();
0042   setDeltaRVetoBarrelNeutrals();
0043   setDeltaRVetoBarrelCharged();
0044   setDeltaRVetoEndcapPhotons();
0045   setDeltaRVetoEndcapNeutrals();
0046   setDeltaRVetoEndcapCharged();
0047 
0048   setRectangleDeltaPhiVetoBarrelPhotons();
0049   setRectangleDeltaPhiVetoBarrelNeutrals();
0050   setRectangleDeltaPhiVetoBarrelCharged();
0051   setRectangleDeltaPhiVetoEndcapPhotons();
0052   setRectangleDeltaPhiVetoEndcapNeutrals();
0053   setRectangleDeltaPhiVetoEndcapCharged();
0054 
0055   setRectangleDeltaEtaVetoBarrelPhotons();
0056   setRectangleDeltaEtaVetoBarrelNeutrals();
0057   setRectangleDeltaEtaVetoBarrelCharged();
0058   setRectangleDeltaEtaVetoEndcapPhotons();
0059   setRectangleDeltaEtaVetoEndcapNeutrals();
0060   setRectangleDeltaEtaVetoEndcapCharged();
0061 
0062   if (bApplyVeto && iParticleType == kElectron) {
0063     //Setup veto conditions for electrons
0064     setDeltaRVetoBarrel(kTRUE);
0065     setDeltaRVetoEndcap(kTRUE);
0066     setRectangleVetoBarrel(kFALSE);
0067     setRectangleVetoEndcap(kFALSE);
0068     setApplyDzDxyVeto(kFALSE);
0069     setApplyPFPUVeto(kTRUE);
0070     setApplyMissHitPhVeto(kTRUE);  //NOTE: decided to go for this on the 26May 2012
0071     //Current recommended default value for the electrons
0072     setUseCrystalSize(kFALSE);
0073 
0074     // setDeltaRVetoBarrelPhotons(1E-5);   //NOTE: just to be in synch with the isoDep: fixed isoDep in 26May
0075     // setDeltaRVetoBarrelCharged(1E-5);    //NOTE: just to be in synch with the isoDep: fixed isoDep in 26May
0076     // setDeltaRVetoBarrelNeutrals(1E-5);   //NOTE: just to be in synch with the isoDep: fixed isoDep in 26May
0077     setDeltaRVetoEndcapPhotons(0.08);
0078     setDeltaRVetoEndcapCharged(0.015);
0079     // setDeltaRVetoEndcapNeutrals(1E-5);  //NOTE: just to be in synch with the isoDep: fixed isoDep in 26May
0080 
0081     setConeSize(0.4);
0082 
0083   } else {
0084     //Setup veto conditions for photons
0085     setApplyDzDxyVeto(kTRUE);
0086     setApplyPFPUVeto(kTRUE);
0087     setApplyMissHitPhVeto(kFALSE);
0088     setDeltaRVetoBarrel(kTRUE);
0089     setDeltaRVetoEndcap(kTRUE);
0090     setRectangleVetoBarrel(kTRUE);
0091     setRectangleVetoEndcap(kFALSE);
0092     setUseCrystalSize(kTRUE);
0093     setConeSize(0.3);
0094 
0095     setDeltaRVetoBarrelPhotons(-1);
0096     setDeltaRVetoBarrelNeutrals(-1);
0097     setDeltaRVetoBarrelCharged(0.02);
0098     setRectangleDeltaPhiVetoBarrelPhotons(1.);
0099     setRectangleDeltaPhiVetoBarrelNeutrals(-1);
0100     setRectangleDeltaPhiVetoBarrelCharged(-1);
0101     setRectangleDeltaEtaVetoBarrelPhotons(0.015);
0102     setRectangleDeltaEtaVetoBarrelNeutrals(-1);
0103     setRectangleDeltaEtaVetoBarrelCharged(-1);
0104 
0105     setDeltaRVetoEndcapPhotons(0.07);
0106     setDeltaRVetoEndcapNeutrals(-1);
0107     setDeltaRVetoEndcapCharged(0.02);
0108     setRectangleDeltaPhiVetoEndcapPhotons(-1);
0109     setRectangleDeltaPhiVetoEndcapNeutrals(-1);
0110     setRectangleDeltaPhiVetoEndcapCharged(-1);
0111     setRectangleDeltaEtaVetoEndcapPhotons(-1);
0112     setRectangleDeltaEtaVetoEndcapNeutrals(-1);
0113     setRectangleDeltaEtaVetoEndcapCharged(-1);
0114     setNumberOfCrystalEndcapPhotons(4.);
0115   }
0116 }
0117 
0118 //--------------------------------------------------------------------------------------------------
0119 void PFIsolationEstimator::initializeElectronIsolation(Bool_t bApplyVeto) {
0120   initialize(bApplyVeto, kElectron);
0121   initializeRings(1, fConeSize);
0122 
0123   //   std::cout << " ********* Init Entering in kElectron setup "
0124   //        << " bApplyVeto " << bApplyVeto
0125   //        << " bDeltaRVetoBarrel " << bDeltaRVetoBarrel
0126   //        << " bDeltaRVetoEndcap " << bDeltaRVetoEndcap
0127   //        << " cone size " << fConeSize
0128   //        << " fDeltaRVetoEndcapPhotons " << fDeltaRVetoEndcapPhotons
0129   //        << " fDeltaRVetoEndcapNeutrals " << fDeltaRVetoEndcapNeutrals
0130   //        << " fDeltaRVetoEndcapCharged " << fDeltaRVetoEndcapCharged << std::endl;
0131 }
0132 
0133 //--------------------------------------------------------------------------------------------------
0134 void PFIsolationEstimator::initializePhotonIsolation(Bool_t bApplyVeto) {
0135   initialize(bApplyVeto, kPhoton);
0136   initializeRings(1, fConeSize);
0137 }
0138 
0139 //--------------------------------------------------------------------------------------------------
0140 void PFIsolationEstimator::initializeElectronIsolationInRings(Bool_t bApplyVeto, int iNumberOfRings, float fRingSize) {
0141   initialize(bApplyVeto, kElectron);
0142   initializeRings(iNumberOfRings, fRingSize);
0143 }
0144 
0145 //--------------------------------------------------------------------------------------------------
0146 void PFIsolationEstimator::initializePhotonIsolationInRings(Bool_t bApplyVeto, int iNumberOfRings, float fRingSize) {
0147   initialize(bApplyVeto, kPhoton);
0148   initializeRings(iNumberOfRings, fRingSize);
0149 }
0150 
0151 //--------------------------------------------------------------------------------------------------
0152 void PFIsolationEstimator::initializeRings(int iNumberOfRings, float fRingSize) {
0153   setRingSize(fRingSize);
0154   setNumbersOfRings(iNumberOfRings);
0155 
0156   fIsolationInRings.clear();
0157   for (int isoBin = 0; isoBin < iNumberOfRings; isoBin++) {
0158     float fTemp = 0.0;
0159     fIsolationInRings.push_back(fTemp);
0160 
0161     float fTempPhoton = 0.0;
0162     fIsolationInRingsPhoton.push_back(fTempPhoton);
0163 
0164     float fTempNeutral = 0.0;
0165     fIsolationInRingsNeutral.push_back(fTempNeutral);
0166 
0167     float fTempCharged = 0.0;
0168     fIsolationInRingsCharged.push_back(fTempCharged);
0169 
0170     float fTempChargedAll = 0.0;
0171     fIsolationInRingsChargedAll.push_back(fTempChargedAll);
0172   }
0173 
0174   fConeSize = fRingSize * (float)iNumberOfRings;
0175 }
0176 
0177 //--------------------------------------------------------------------------------------------------
0178 float PFIsolationEstimator::fGetIsolation(const reco::PFCandidate* pfCandidate,
0179                                           const reco::PFCandidateCollection* pfParticlesColl,
0180                                           reco::VertexRef vtx,
0181                                           edm::Handle<reco::VertexCollection> vertices) {
0182   fGetIsolationInRings(pfCandidate, pfParticlesColl, vtx, vertices);
0183   refSC = reco::SuperClusterRef();
0184   fIsolation = fIsolationInRings[0];
0185 
0186   return fIsolation;
0187 }
0188 
0189 //--------------------------------------------------------------------------------------------------
0190 std::vector<float> PFIsolationEstimator::fGetIsolationInRings(const reco::PFCandidate* pfCandidate,
0191                                                               const reco::PFCandidateCollection* pfParticlesColl,
0192                                                               reco::VertexRef vtx,
0193                                                               edm::Handle<reco::VertexCollection> vertices) {
0194   int isoBin;
0195 
0196   for (isoBin = 0; isoBin < iNumberOfRings; isoBin++) {
0197     fIsolationInRings[isoBin] = 0.;
0198     fIsolationInRingsPhoton[isoBin] = 0.;
0199     fIsolationInRingsNeutral[isoBin] = 0.;
0200     fIsolationInRingsCharged[isoBin] = 0.;
0201     fIsolationInRingsChargedAll[isoBin] = 0.;
0202   }
0203 
0204   fEta = pfCandidate->eta();
0205   fPhi = pfCandidate->phi();
0206   fPt = pfCandidate->pt();
0207   fVx = pfCandidate->vx();
0208   fVy = pfCandidate->vy();
0209   fVz = pfCandidate->vz();
0210 
0211   pivotInBarrel = std::abs(pfCandidate->positionAtECALEntrance().eta()) < 1.479;
0212 
0213   for (unsigned iPF = 0; iPF < pfParticlesColl->size(); iPF++) {
0214     const reco::PFCandidate& pfParticle = (*pfParticlesColl)[iPF];
0215 
0216     if (&pfParticle == (pfCandidate))
0217       continue;
0218 
0219     if (pfParticle.pdgId() == 22) {
0220       if (isPhotonParticleVetoed(&pfParticle) >= 0.) {
0221         isoBin = (int)(fDeltaR / fRingSize);
0222         fIsolationInRingsPhoton[isoBin] = fIsolationInRingsPhoton[isoBin] + pfParticle.pt();
0223       }
0224 
0225     } else if (std::abs(pfParticle.pdgId()) == 130) {
0226       if (isNeutralParticleVetoed(&pfParticle) >= 0.) {
0227         isoBin = (int)(fDeltaR / fRingSize);
0228         fIsolationInRingsNeutral[isoBin] = fIsolationInRingsNeutral[isoBin] + pfParticle.pt();
0229       }
0230 
0231       //}else if(std::abs(pfParticle.pdgId()) == 11 ||abs(pfParticle.pdgId()) == 13 || std::abs(pfParticle.pdgId()) == 211){
0232     } else if (std::abs(pfParticle.pdgId()) == 211) {
0233       if (isChargedParticleVetoed(&pfParticle, vtx, vertices) >= 0.) {
0234         isoBin = (int)(fDeltaR / fRingSize);
0235         fIsolationInRingsCharged[isoBin] = fIsolationInRingsCharged[isoBin] + pfParticle.pt();
0236       }
0237     }
0238   }
0239 
0240   for (int isoBin = 0; isoBin < iNumberOfRings; isoBin++) {
0241     fIsolationInRings[isoBin] =
0242         fIsolationInRingsPhoton[isoBin] + fIsolationInRingsNeutral[isoBin] + fIsolationInRingsCharged[isoBin];
0243   }
0244 
0245   return fIsolationInRings;
0246 }
0247 
0248 //--------------------------------------------------------------------------------------------------
0249 float PFIsolationEstimator::fGetIsolation(const reco::Photon* photon,
0250                                           const reco::PFCandidateCollection* pfParticlesColl,
0251                                           reco::VertexRef vtx,
0252                                           edm::Handle<reco::VertexCollection> vertices) {
0253   fGetIsolationInRings(photon, pfParticlesColl, vtx, vertices);
0254   fIsolation = fIsolationInRings[0];
0255 
0256   return fIsolation;
0257 }
0258 
0259 //--------------------------------------------------------------------------------------------------
0260 std::vector<float> PFIsolationEstimator::fGetIsolationInRings(const reco::Photon* photon,
0261                                                               const reco::PFCandidateCollection* pfParticlesColl,
0262                                                               reco::VertexRef vtx,
0263                                                               edm::Handle<reco::VertexCollection> vertices) {
0264   int isoBin;
0265 
0266   for (isoBin = 0; isoBin < iNumberOfRings; isoBin++) {
0267     fIsolationInRings[isoBin] = 0.;
0268     fIsolationInRingsPhoton[isoBin] = 0.;
0269     fIsolationInRingsNeutral[isoBin] = 0.;
0270     fIsolationInRingsCharged[isoBin] = 0.;
0271     fIsolationInRingsChargedAll[isoBin] = 0.;
0272   }
0273 
0274   iMissHits = 0;
0275 
0276   refSC = photon->superCluster();
0277   pivotInBarrel = std::abs((refSC->position().eta())) < 1.479;
0278 
0279   for (unsigned iPF = 0; iPF < pfParticlesColl->size(); iPF++) {
0280     const reco::PFCandidate& pfParticle = (*pfParticlesColl)[iPF];
0281 
0282     if (pfParticle.superClusterRef().isNonnull() && photon->superCluster().isNonnull() &&
0283         pfParticle.superClusterRef() == photon->superCluster())
0284       continue;
0285 
0286     if (pfParticle.pdgId() == 22) {
0287       // Set the vertex of reco::Photon to the first PV
0288       math::XYZVector direction = math::XYZVector(photon->superCluster()->x() - pfParticle.vx(),
0289                                                   photon->superCluster()->y() - pfParticle.vy(),
0290                                                   photon->superCluster()->z() - pfParticle.vz());
0291 
0292       fEta = direction.Eta();
0293       fPhi = direction.Phi();
0294       fVx = pfParticle.vx();
0295       fVy = pfParticle.vy();
0296       fVz = pfParticle.vz();
0297 
0298       if (isPhotonParticleVetoed(&pfParticle) >= 0.) {
0299         isoBin = (int)(fDeltaR / fRingSize);
0300         fIsolationInRingsPhoton[isoBin] = fIsolationInRingsPhoton[isoBin] + pfParticle.pt();
0301       }
0302 
0303     } else if (std::abs(pfParticle.pdgId()) == 130) {
0304       // Set the vertex of reco::Photon to the first PV
0305       math::XYZVector direction = math::XYZVector(photon->superCluster()->x() - pfParticle.vx(),
0306                                                   photon->superCluster()->y() - pfParticle.vy(),
0307                                                   photon->superCluster()->z() - pfParticle.vz());
0308 
0309       fEta = direction.Eta();
0310       fPhi = direction.Phi();
0311       fVx = pfParticle.vx();
0312       fVy = pfParticle.vy();
0313       fVz = pfParticle.vz();
0314 
0315       if (isNeutralParticleVetoed(&pfParticle) >= 0.) {
0316         isoBin = (int)(fDeltaR / fRingSize);
0317         fIsolationInRingsNeutral[isoBin] = fIsolationInRingsNeutral[isoBin] + pfParticle.pt();
0318       }
0319 
0320       //}else if(std::abs(pfParticle.pdgId()) == 11 ||abs(pfParticle.pdgId()) == 13 || std::abs(pfParticle.pdgId()) == 211){
0321     } else if (std::abs(pfParticle.pdgId()) == 211) {
0322       // Set the vertex of reco::Photon to the first PV
0323       math::XYZVector direction = math::XYZVector(photon->superCluster()->x() - (*vtx).x(),
0324                                                   photon->superCluster()->y() - (*vtx).y(),
0325                                                   photon->superCluster()->z() - (*vtx).z());
0326 
0327       fEta = direction.Eta();
0328       fPhi = direction.Phi();
0329       fVx = (*vtx).x();
0330       fVy = (*vtx).y();
0331       fVz = (*vtx).z();
0332 
0333       if (isChargedParticleVetoed(&pfParticle, vtx, vertices) >= 0.) {
0334         isoBin = (int)(fDeltaR / fRingSize);
0335         fIsolationInRingsCharged[isoBin] = fIsolationInRingsCharged[isoBin] + pfParticle.pt();
0336       }
0337     }
0338   }
0339 
0340   for (int isoBin = 0; isoBin < iNumberOfRings; isoBin++) {
0341     fIsolationInRings[isoBin] =
0342         fIsolationInRingsPhoton[isoBin] + fIsolationInRingsNeutral[isoBin] + fIsolationInRingsCharged[isoBin];
0343   }
0344 
0345   return fIsolationInRings;
0346 }
0347 
0348 //--------------------------------------------------------------------------------------------------
0349 float PFIsolationEstimator::fGetIsolation(const reco::GsfElectron* electron,
0350                                           const reco::PFCandidateCollection* pfParticlesColl,
0351                                           reco::VertexRef vtx,
0352                                           edm::Handle<reco::VertexCollection> vertices) {
0353   fGetIsolationInRings(electron, pfParticlesColl, vtx, vertices);
0354   fIsolation = fIsolationInRings[0];
0355 
0356   return fIsolation;
0357 }
0358 
0359 //--------------------------------------------------------------------------------------------------
0360 std::vector<float> PFIsolationEstimator::fGetIsolationInRings(const reco::GsfElectron* electron,
0361                                                               const reco::PFCandidateCollection* pfParticlesColl,
0362                                                               reco::VertexRef vtx,
0363                                                               edm::Handle<reco::VertexCollection> vertices) {
0364   int isoBin;
0365 
0366   for (isoBin = 0; isoBin < iNumberOfRings; isoBin++) {
0367     fIsolationInRings[isoBin] = 0.;
0368     fIsolationInRingsPhoton[isoBin] = 0.;
0369     fIsolationInRingsNeutral[isoBin] = 0.;
0370     fIsolationInRingsCharged[isoBin] = 0.;
0371     fIsolationInRingsChargedAll[isoBin] = 0.;
0372   }
0373 
0374   //  int iMatch =  matchPFObject(electron,pfParticlesColl);
0375 
0376   fEta = electron->eta();
0377   fPhi = electron->phi();
0378   fPt = electron->pt();
0379   fVx = electron->vx();
0380   fVy = electron->vy();
0381   fVz = electron->vz();
0382   iMissHits = electron->gsfTrack()->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
0383 
0384   //  if(electron->ecalDrivenSeed())
0385   refSC = electron->superCluster();
0386   pivotInBarrel = std::abs((refSC->position().eta())) < 1.479;
0387 
0388   for (unsigned iPF = 0; iPF < pfParticlesColl->size(); iPF++) {
0389     const reco::PFCandidate& pfParticle = (*pfParticlesColl)[iPF];
0390 
0391     if (pfParticle.pdgId() == 22) {
0392       if (isPhotonParticleVetoed(&pfParticle) >= 0.) {
0393         isoBin = (int)(fDeltaR / fRingSize);
0394         fIsolationInRingsPhoton[isoBin] = fIsolationInRingsPhoton[isoBin] + pfParticle.pt();
0395       }
0396 
0397     } else if (std::abs(pfParticle.pdgId()) == 130) {
0398       if (isNeutralParticleVetoed(&pfParticle) >= 0.) {
0399         isoBin = (int)(fDeltaR / fRingSize);
0400         fIsolationInRingsNeutral[isoBin] = fIsolationInRingsNeutral[isoBin] + pfParticle.pt();
0401       }
0402 
0403       //}else if(std::abs(pfParticle.pdgId()) == 11 ||abs(pfParticle.pdgId()) == 13 || std::abs(pfParticle.pdgId()) == 211){
0404     } else if (std::abs(pfParticle.pdgId()) == 211) {
0405       if (isChargedParticleVetoed(&pfParticle, vtx, vertices) >= 0.) {
0406         isoBin = (int)(fDeltaR / fRingSize);
0407 
0408         fIsolationInRingsCharged[isoBin] = fIsolationInRingsCharged[isoBin] + pfParticle.pt();
0409       }
0410     }
0411   }
0412 
0413   for (int isoBin = 0; isoBin < iNumberOfRings; isoBin++) {
0414     fIsolationInRings[isoBin] =
0415         fIsolationInRingsPhoton[isoBin] + fIsolationInRingsNeutral[isoBin] + fIsolationInRingsCharged[isoBin];
0416   }
0417 
0418   return fIsolationInRings;
0419 }
0420 
0421 //--------------------------------------------------------------------------------------------------
0422 float PFIsolationEstimator::isPhotonParticleVetoed(const reco::PFCandidate* pfIsoCand) {
0423   fDeltaR = deltaR(fEta, fPhi, pfIsoCand->eta(), pfIsoCand->phi());
0424 
0425   if (fDeltaR > fConeSize)
0426     return -999.;
0427 
0428   fDeltaPhi = deltaPhi(fPhi, pfIsoCand->phi());
0429   fDeltaEta = fEta - pfIsoCand->eta();
0430 
0431   if (!bApplyVeto)
0432     return fDeltaR;
0433 
0434   //NOTE: get the direction for the EB/EE transition region from the deposit just to be in synch with the isoDep
0435   //      this will be changed in the future
0436 
0437   if (bApplyMissHitPhVeto) {
0438     if (iMissHits > 0)
0439       if (pfIsoCand->mva_nothing_gamma() > 0.99) {
0440         if (pfIsoCand->superClusterRef().isNonnull() && refSC.isNonnull()) {
0441           if (pfIsoCand->superClusterRef() == refSC)
0442             return -999.;
0443         }
0444       }
0445   }
0446 
0447   if (pivotInBarrel) {
0448     if (bDeltaRVetoBarrel) {
0449       if (fDeltaR < fDeltaRVetoBarrelPhotons)
0450         return -999.;
0451     }
0452 
0453     if (bRectangleVetoBarrel) {
0454       if (std::abs(fDeltaEta) < fRectangleDeltaEtaVetoBarrelPhotons &&
0455           std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoBarrelPhotons) {
0456         return -999.;
0457       }
0458     }
0459   } else {
0460     if (bUseCrystalSize == true) {
0461       fDeltaRVetoEndcapPhotons = 0.00864 * std::abs(sinh(refSC->position().eta())) * fNumberOfCrystalEndcapPhotons;
0462     }
0463 
0464     if (bDeltaRVetoEndcap) {
0465       if (fDeltaR < fDeltaRVetoEndcapPhotons)
0466         return -999.;
0467     }
0468     if (bRectangleVetoEndcap) {
0469       if (std::abs(fDeltaEta) < fRectangleDeltaEtaVetoEndcapPhotons &&
0470           std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoEndcapPhotons) {
0471         return -999.;
0472       }
0473     }
0474   }
0475 
0476   return fDeltaR;
0477 }
0478 
0479 //--------------------------------------------------------------------------------------------------
0480 float PFIsolationEstimator::isNeutralParticleVetoed(const reco::PFCandidate* pfIsoCand) {
0481   fDeltaR = deltaR(fEta, fPhi, pfIsoCand->eta(), pfIsoCand->phi());
0482 
0483   if (fDeltaR > fConeSize)
0484     return -999;
0485 
0486   fDeltaPhi = deltaPhi(fPhi, pfIsoCand->phi());
0487   fDeltaEta = fEta - pfIsoCand->eta();
0488 
0489   if (!bApplyVeto)
0490     return fDeltaR;
0491 
0492   //NOTE: get the direction for the EB/EE transition region from the deposit just to be in synch with the isoDep
0493   //      this will be changed in the future
0494   if (pivotInBarrel) {
0495     if (!bDeltaRVetoBarrel && !bRectangleVetoBarrel) {
0496       return fDeltaR;
0497     }
0498 
0499     if (bDeltaRVetoBarrel) {
0500       if (fDeltaR < fDeltaRVetoBarrelNeutrals)
0501         return -999.;
0502     }
0503     if (bRectangleVetoBarrel) {
0504       if (std::abs(fDeltaEta) < fRectangleDeltaEtaVetoBarrelNeutrals &&
0505           std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoBarrelNeutrals) {
0506         return -999.;
0507       }
0508     }
0509 
0510   } else {
0511     if (!bDeltaRVetoEndcap && !bRectangleVetoEndcap) {
0512       return fDeltaR;
0513     }
0514     if (bDeltaRVetoEndcap) {
0515       if (fDeltaR < fDeltaRVetoEndcapNeutrals)
0516         return -999.;
0517     }
0518     if (bRectangleVetoEndcap) {
0519       if (std::abs(fDeltaEta) < fRectangleDeltaEtaVetoEndcapNeutrals &&
0520           std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoEndcapNeutrals) {
0521         return -999.;
0522       }
0523     }
0524   }
0525 
0526   return fDeltaR;
0527 }
0528 
0529 //----------------------------------------------------------------------------------------------------
0530 float PFIsolationEstimator::isChargedParticleVetoed(const reco::PFCandidate* pfIsoCand,
0531                                                     edm::Handle<reco::VertexCollection> vertices) {
0532   //need code to handle special conditions
0533 
0534   return -999;
0535 }
0536 
0537 //-----------------------------------------------------------------------------------------------------
0538 float PFIsolationEstimator::isChargedParticleVetoed(const reco::PFCandidate* pfIsoCand,
0539                                                     reco::VertexRef vtxMain,
0540                                                     edm::Handle<reco::VertexCollection> vertices) {
0541   reco::VertexRef vtx = chargedHadronVertex(vertices, *pfIsoCand);
0542   if (vtx.isNull())
0543     return -999.;
0544 
0545   //   float fVtxMainX = (*vtxMain).x();
0546   //   float fVtxMainY = (*vtxMain).y();
0547   float fVtxMainZ = (*vtxMain).z();
0548 
0549   if (bApplyPFPUVeto) {
0550     if (vtx != vtxMain)
0551       return -999.;
0552   }
0553 
0554   if (bApplyDzDxyVeto) {
0555     if (iParticleType == kPhoton) {
0556       float dz = std::abs(pfIsoCand->trackRef()->dz((*vtxMain).position()));
0557       if (dz > 0.2)
0558         return -999.;
0559 
0560       double dxy = pfIsoCand->trackRef()->dxy((*vtxMain).position());
0561       if (std::abs(dxy) > 0.1)
0562         return -999.;
0563 
0564       /*
0565       float dz = std::abs(vtx->z() - fVtxMainZ);
0566       if (dz > 1.)
0567     return -999.;
0568       
0569       
0570       double dxy = ( -(vtx->x() - fVtxMainX)*pfIsoCand->py() + (vtx->y() - fVtxMainY)*pfIsoCand->px()) / pfIsoCand->pt();
0571       
0572       if(std::abs(dxy) > 0.2)
0573     return -999.;
0574       */
0575     } else {
0576       float dz = std::abs(vtx->z() - fVtxMainZ);
0577       if (dz > 1.)
0578         return -999.;
0579 
0580       double dxy = (-(vtx->x() - fVx) * pfIsoCand->py() + (vtx->y() - fVy) * pfIsoCand->px()) / pfIsoCand->pt();
0581       if (std::abs(dxy) > 0.1)
0582         return -999.;
0583     }
0584   }
0585 
0586   fDeltaR = deltaR(pfIsoCand->eta(), pfIsoCand->phi(), fEta, fPhi);
0587 
0588   if (fDeltaR > fConeSize)
0589     return -999.;
0590 
0591   fDeltaPhi = deltaPhi(fPhi, pfIsoCand->phi());
0592   fDeltaEta = fEta - pfIsoCand->eta();
0593 
0594   //   std::cout << " charged hadron: DR " <<  fDeltaR
0595   //        << " pt " <<  pfIsoCand->pt() << " eta,phi " << pfIsoCand->eta() << ", " << pfIsoCand->phi()
0596   //        << " fVtxMainZ " << (*vtxMain).z() << " cand z " << vtx->z() << std::endl;
0597 
0598   if (!bApplyVeto)
0599     return fDeltaR;
0600 
0601   //NOTE: get the direction for the EB/EE transition region from the deposit just to be in synch with the isoDep
0602   //      this will be changed in the future
0603   if (pivotInBarrel) {
0604     if (!bDeltaRVetoBarrel && !bRectangleVetoBarrel) {
0605       return fDeltaR;
0606     }
0607 
0608     if (bDeltaRVetoBarrel) {
0609       if (fDeltaR < fDeltaRVetoBarrelCharged)
0610         return -999.;
0611     }
0612     if (bRectangleVetoBarrel) {
0613       if (std::abs(fDeltaEta) < fRectangleDeltaEtaVetoBarrelCharged &&
0614           std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoBarrelCharged) {
0615         return -999.;
0616       }
0617     }
0618 
0619   } else {
0620     if (!bDeltaRVetoEndcap && !bRectangleVetoEndcap) {
0621       return fDeltaR;
0622     }
0623     if (bDeltaRVetoEndcap) {
0624       if (fDeltaR < fDeltaRVetoEndcapCharged)
0625         return -999.;
0626     }
0627     if (bRectangleVetoEndcap) {
0628       if (std::abs(fDeltaEta) < fRectangleDeltaEtaVetoEndcapCharged &&
0629           std::abs(fDeltaPhi) < fRectangleDeltaPhiVetoEndcapCharged) {
0630         return -999.;
0631       }
0632     }
0633   }
0634 
0635   return fDeltaR;
0636 }
0637 
0638 //--------------------------------------------------------------------------------------------------
0639 reco::VertexRef PFIsolationEstimator::chargedHadronVertex(edm::Handle<reco::VertexCollection> verticesColl,
0640                                                           const reco::PFCandidate& pfcand) {
0641   //code copied from Florian's PFNoPU class
0642 
0643   reco::TrackBaseRef trackBaseRef(pfcand.trackRef());
0644 
0645   size_t iVertex = 0;
0646   unsigned index = 0;
0647   unsigned nFoundVertex = 0;
0648 
0649   float bestweight = 0;
0650 
0651   const reco::VertexCollection& vertices = *(verticesColl.product());
0652 
0653   for (reco::VertexCollection::const_iterator iv = vertices.begin(); iv != vertices.end(); ++iv, ++index) {
0654     const reco::Vertex& vtx = *iv;
0655 
0656     // loop on tracks in vertices
0657     for (reco::Vertex::trackRef_iterator iTrack = vtx.tracks_begin(); iTrack != vtx.tracks_end(); ++iTrack) {
0658       const reco::TrackBaseRef& baseRef = *iTrack;
0659 
0660       // one of the tracks in the vertex is the same as
0661       // the track considered in the function
0662       if (baseRef == trackBaseRef) {
0663         float w = vtx.trackWeight(baseRef);
0664         //select the vertex for which the track has the highest weight
0665         if (w > bestweight) {
0666           bestweight = w;
0667           iVertex = index;
0668           nFoundVertex++;
0669         }
0670       }
0671     }
0672   }
0673 
0674   if (nFoundVertex > 0) {
0675     if (nFoundVertex != 1)
0676       edm::LogWarning("TrackOnTwoVertex") << "a track is shared by at least two verteces. Used to be an assert";
0677     return reco::VertexRef(verticesColl, iVertex);
0678   }
0679   // no vertex found with this track.
0680 
0681   // optional: as a secondary solution, associate the closest vertex in z
0682   if (checkClosestZVertex) {
0683     double dzmin = 10000.;
0684     double ztrack = pfcand.vertex().z();
0685     bool foundVertex = false;
0686     index = 0;
0687     for (reco::VertexCollection::const_iterator iv = vertices.begin(); iv != vertices.end(); ++iv, ++index) {
0688       double dz = std::abs(ztrack - iv->z());
0689       if (dz < dzmin) {
0690         dzmin = dz;
0691         iVertex = index;
0692         foundVertex = true;
0693       }
0694     }
0695 
0696     if (foundVertex)
0697       return reco::VertexRef(verticesColl, iVertex);
0698   }
0699 
0700   return reco::VertexRef();
0701 }
0702 
0703 int PFIsolationEstimator::matchPFObject(const reco::Photon* photon, const reco::PFCandidateCollection* Candidates) {
0704   Int_t iMatch = -1;
0705 
0706   int i = 0;
0707   for (reco::PFCandidateCollection::const_iterator iPF = Candidates->begin(); iPF != Candidates->end(); iPF++) {
0708     const reco::PFCandidate& pfParticle = (*iPF);
0709     //    if((((pfParticle.pdgId()==22 && pfParticle.mva_nothing_gamma()>0.01) || TMath::Abs(pfParticle.pdgId())==11) )){
0710     if ((((pfParticle.pdgId() == 22) || TMath::Abs(pfParticle.pdgId()) == 11))) {
0711       if (pfParticle.superClusterRef() == photon->superCluster())
0712         iMatch = i;
0713     }
0714 
0715     i++;
0716   }
0717 
0718   /*
0719   if(iMatch == -1){
0720     i=0;
0721     float fPt = -1;
0722     for(reco::PFCandidateCollection::const_iterator iPF=Candidates->begin();iPF !=Candidates->end();iPF++){
0723       const reco::PFCandidate& pfParticle = (*iPF);
0724       if((((pfParticle.pdgId()==22 ) || TMath::Abs(pfParticle.pdgId())==11) )){
0725     if(pfParticle.pt()>fPt){
0726       fDeltaR = deltaR(pfParticle.eta(),pfParticle.phi(),photon->eta(),photon->phi());
0727       if(fDeltaR<0.1){
0728         iMatch = i;
0729         fPt = pfParticle.pt();
0730       }
0731     }
0732       }
0733       i++;
0734     }
0735   }
0736 */
0737 
0738   return iMatch;
0739 }
0740 
0741 int PFIsolationEstimator::matchPFObject(const reco::GsfElectron* electron,
0742                                         const reco::PFCandidateCollection* Candidates) {
0743   Int_t iMatch = -1;
0744 
0745   int i = 0;
0746   for (reco::PFCandidateCollection::const_iterator iPF = Candidates->begin(); iPF != Candidates->end(); iPF++) {
0747     const reco::PFCandidate& pfParticle = (*iPF);
0748     //    if((((pfParticle.pdgId()==22 && pfParticle.mva_nothing_gamma()>0.01) || TMath::Abs(pfParticle.pdgId())==11) )){
0749     if ((((pfParticle.pdgId() == 22) || TMath::Abs(pfParticle.pdgId()) == 11))) {
0750       if (pfParticle.superClusterRef() == electron->superCluster())
0751         iMatch = i;
0752     }
0753 
0754     i++;
0755   }
0756 
0757   if (iMatch == -1) {
0758     i = 0;
0759     float fPt = -1;
0760     for (reco::PFCandidateCollection::const_iterator iPF = Candidates->begin(); iPF != Candidates->end(); iPF++) {
0761       const reco::PFCandidate& pfParticle = (*iPF);
0762       if ((((pfParticle.pdgId() == 22) || TMath::Abs(pfParticle.pdgId()) == 11))) {
0763         if (pfParticle.pt() > fPt) {
0764           fDeltaR = deltaR(pfParticle.eta(), pfParticle.phi(), electron->eta(), electron->phi());
0765           if (fDeltaR < 0.1) {
0766             iMatch = i;
0767             fPt = pfParticle.pt();
0768           }
0769         }
0770       }
0771       i++;
0772     }
0773   }
0774 
0775   return iMatch;
0776 }
[ 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!