Utils/src/IsolationTools.cc

// $Id: IsolationTools.cc,v 1.8 2011/02/21 13:50:20 ceballos Exp $

#include "MitPhysics/Utils/interface/IsolationTools.h"
#include "MitCommon/MathTools/interface/MathUtils.h"

ClassImp(mithep::IsolationTools)

using namespace mithep;

//--------------------------------------------------------------------------------------------------
Double_t IsolationTools::TrackIsolation(const Track *p, Double_t extRadius, Double_t intRadius,
                                        Double_t ptLow, Double_t maxVtxZDist, 
                                        const Collection<Track> *tracks) 
{
  //Computes the Track Isolation: Summed Transverse Momentum of all tracks inside an 
  //annulus around the electron seed track.  

  Double_t ptSum =0.;  
  for (UInt_t i=0; i<tracks->GetEntries();i++) {   
    Double_t tmpPt = tracks->At(i)->Pt();
    Double_t deltaZ = fabs(p->Z0() - tracks->At(i)->Z0());

    //ignore the track if it is below the pt threshold
    if (tmpPt < ptLow) 
      continue;    
    //ingore the track if it is too far away in Z
    if (deltaZ > maxVtxZDist) 
      continue;
           
    Double_t dr = MathUtils::DeltaR(p->Phi(),p->Eta(),tracks->At(i)->Phi(), tracks->At(i)->Eta());
    //add the track pt if it is inside the annulus
    if ( dr < extRadius && 
         dr >= intRadius ) {
      ptSum += tmpPt;
    }
  }
  return ptSum;  
}

//--------------------------------------------------------------------------------------------------
Double_t IsolationTools::EcalIsolation(const SuperCluster *sc, Double_t coneSize, Double_t etLow, 
                                       const Collection<BasicCluster> *basicClusters) 
{
  //Computes the Ecal Isolation: Summed Transverse Energy of all Basic Clusters inside a  
  //cone around the electron, excluding those that are inside the electron super cluster.

  Double_t ecalIsol=0.;
  const BasicCluster *basicCluster= 0;
  for (UInt_t i=0; i<basicClusters->GetEntries();i++) {    
    basicCluster = basicClusters->At(i);    
    Double_t basicClusterEnergy    = basicCluster->Energy();
    Double_t basicClusterEta  = basicCluster->Eta();
    Double_t basicClusterEt   = basicClusterEnergy*sin(2*atan(exp(basicClusterEta)));           

    if (basicClusterEt > etLow) {            
      bool inSuperCluster = false;        
      
      // loop over the basic clusters of the supercluster
      // to make sure that the basic cluster is not inside
      // the super cluster. We exclude those.
      for (UInt_t j=0; j<sc->ClusterSize(); j++) {
        const BasicCluster *tempBasicClusterInSuperCluster = sc->Cluster(j);    
        if (tempBasicClusterInSuperCluster == basicCluster) {
          inSuperCluster = true;            
        }
      }
      
      if (!inSuperCluster) {        
        Double_t dr = MathUtils::DeltaR(sc->Phi(), sc->Eta(),
                                      basicCluster->Phi(),basicCluster->Eta());
        if(dr < coneSize) {
          ecalIsol += basicClusterEt;
        }
      }
    }
  } 
  return ecalIsol;
}

//--------------------------------------------------------------------------------------------------
Double_t IsolationTools::CaloTowerHadIsolation(const ThreeVector *p, Double_t extRadius, 
                                               Double_t intRadius, Double_t etLow, 
                                               const Collection<CaloTower> *caloTowers) 
{
  //Computes the CaloTower Had Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers 
  //inside an annulus around the electron super cluster position.

  Double_t sumEt = 0;
  for (UInt_t i=0; i<caloTowers->GetEntries();i++) {    
    Double_t caloTowerEt = caloTowers->At(i)->HadEt();
    Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
                                  p->Phi(), p->Eta());
    if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
      sumEt += caloTowerEt;
    }
  }
  return sumEt;
}

//--------------------------------------------------------------------------------------------------
Double_t IsolationTools::CaloTowerEmIsolation(const ThreeVector *p, Double_t extRadius, 
                                              Double_t intRadius, Double_t etLow, 
                                              const Collection<CaloTower> *caloTowers) 
{
  //Computes the CaloTower Em Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers 
  //inside an annulus around the electron super cluster position.

  Double_t sumEt = 0;
  for (UInt_t i=0; i<caloTowers->GetEntries();i++) {    
    Double_t caloTowerEt = caloTowers->At(i)->EmEt();
    Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
                                   p->Phi(), p->Eta());
    if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
      sumEt += caloTowerEt;
    }
  }
  return sumEt;
}

//--------------------------------------------------------------------------------------------------
Double_t IsolationTools::PFMuonIsolation(const Muon *p, const Collection<PFCandidate> *PFCands, 
                                         const Vertex *vertex, Double_t  delta_z, Double_t ptMin,
                                         Double_t extRadius, Double_t intRadius, int isoType,
                                         Double_t beta, const MuonCol *goodMuons, 
                                         const ElectronCol *goodElectrons)
{
  //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an 
  //annulus around the particle seed track.  

  Double_t zLepton = 0.0;
  if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);

  Double_t ptSum =0.;  
  for (UInt_t i=0; i<PFCands->GetEntries();i++) {   
    const PFCandidate *pf = PFCands->At(i);
    
    Bool_t isGoodType = kFALSE;
    // all particles
    if     (isoType == 0)                                          isGoodType = kTRUE;
    // charged particles only
    else if(isoType == 1 && pf->BestTrk())                         isGoodType = kTRUE;
    // charged particles and gammas only
    else if(isoType == 2 && 
           (pf->BestTrk() || pf->PFType() == PFCandidate::eGamma)) isGoodType = kTRUE;
     // all particles, rejecting good leptons
    else if(isoType == 3)                                          isGoodType = kTRUE;

    if(isGoodType == kFALSE) continue;

    if(pf->Pt() <= ptMin) continue;

    if(pf->TrackerTrk() && p->TrackerTrk() &&
       pf->TrackerTrk() == p->TrackerTrk()) continue;

    Double_t deltaZ = 0.0;
    if(pf->BestTrk()) {
      deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
    }

    // ignore the pf candidate if it is too far away in Z
    if (deltaZ >= delta_z) 
      continue;
           
    Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
    // add the pf pt if it is inside the extRadius and outside the intRadius
    if ( dr < extRadius && 
         dr >= intRadius ) {
      Bool_t isLepton = kFALSE;
      if(goodMuons && isoType == 3){
        for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
          const Muon *m = goodMuons->At(nl);
          if(pf->TrackerTrk() && m->TrackerTrk() &&
             pf->TrackerTrk() == m->TrackerTrk()) {
            isLepton = kTRUE;
            break;
          }
        }
      }
      if(goodElectrons && isLepton == kFALSE && isoType == 3){
        for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
          const Electron *e = goodElectrons->At(nl);
          if(pf->TrackerTrk() && e->TrackerTrk() &&
             pf->TrackerTrk() == e->TrackerTrk()) {
            isLepton = kTRUE;
            break;
          }
          if(pf->GsfTrk() && e->GsfTrk() &&
             pf->GsfTrk() == e->GsfTrk()) {
            isLepton = kTRUE;
            break;
          }
        }
      }
      if(isLepton == kFALSE){
        if(pf->BestTrk()) ptSum += pf->Pt();
        else              ptSum += pf->Pt()*beta;
      }
    }
  }
  return ptSum;
}
//--------------------------------------------------------------------------------------------------
Double_t IsolationTools::PFElectronIsolation(const Electron *p, const PFCandidateCol *PFCands, 
                                             const Vertex *vertex, Double_t delta_z, Double_t ptMin,
                                             Double_t extRadius, Double_t intRadius, int isoType,
                                             Double_t beta, const MuonCol *goodMuons, 
                                             const ElectronCol *goodElectrons)
{
  //Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an 
  //annulus around the particle seed track.  

  Double_t zLepton = 0.0;
  if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);

  Double_t ptSum =0.;  
  for (UInt_t i=0; i<PFCands->GetEntries();i++) {   
    const PFCandidate *pf = PFCands->At(i);
    
    Bool_t isGoodType = kFALSE;
    // all particles
    if     (isoType == 0)                                          isGoodType = kTRUE;
    // charged particles only
    else if(isoType == 1 && pf->BestTrk())                         isGoodType = kTRUE;
    // charged particles and gammas only
    else if(isoType == 2 && 
           (pf->BestTrk() || pf->PFType() == PFCandidate::eGamma)) isGoodType = kTRUE;
    // all particles, rejecting good leptons
    else if(isoType == 3)                                          isGoodType = kTRUE;

    if(isGoodType == kFALSE) continue;

    if(pf->Pt() <= ptMin) continue;

    if(pf->TrackerTrk() && p->TrackerTrk() &&
       pf->TrackerTrk() == p->TrackerTrk()) continue;

    if(pf->GsfTrk() && p->GsfTrk() &&
       pf->GsfTrk() == p->GsfTrk()) continue;

    Double_t deltaZ = 0.0;
    if(pf->BestTrk()) {
      deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
    }

    // ignore the pf candidate if it is too far away in Z
    if (deltaZ >= delta_z) 
      continue;
           
    Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
    // add the pf pt if it is inside the extRadius and outside the intRadius
    if ( dr < extRadius && 
         dr >= intRadius ) {
      Bool_t isLepton = kFALSE;
      if(goodMuons && isoType == 3){
        for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
          const Muon *m = goodMuons->At(nl);
          if(pf->TrackerTrk() && m->TrackerTrk() &&
             pf->TrackerTrk() == m->TrackerTrk()) {
            isLepton = kTRUE;
            break;
          }
        }
      }
      if(goodElectrons && isLepton == kFALSE && isoType == 3){
        for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
          const Electron *e = goodElectrons->At(nl);
          if(pf->TrackerTrk() && e->TrackerTrk() &&
             pf->TrackerTrk() == e->TrackerTrk()) {
            isLepton = kTRUE;
            break;
          }
          if(pf->GsfTrk() && e->GsfTrk() &&
             pf->GsfTrk() == e->GsfTrk()) {
            isLepton = kTRUE;
            break;
          }
        }
      }
      if(isLepton == kFALSE){
        if(pf->BestTrk()) ptSum += pf->Pt();
        else              ptSum += pf->Pt()*beta;
      }
    }
  }
  return ptSum;
}
//--------------------------------------------------------------------------------------------------
Double_t IsolationTools::BetaM(const TrackCol *tracks, const Muon *p, const Vertex *vertex, 
                               Double_t ptMin, Double_t  delta_z, Double_t extRadius,
                               Double_t intRadius){

  if(!tracks) return 1.0;
  if(tracks->GetEntries() <= 0) return 1.0;

  double Pt_jets_X = 0. ;
  double Pt_jets_Y = 0. ;
  double Pt_jets_X_tot = 0. ;
  double Pt_jets_Y_tot = 0. ;

  for(int i=0;i<int(tracks->GetEntries());i++){
    const Track *pTrack = tracks->At(i);

    if(pTrack && p->TrackerTrk() &&
       pTrack == p->TrackerTrk()) continue;

    if(pTrack->Pt() <= ptMin) continue;

    Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
    if ( dr < extRadius && dr >= intRadius ) {
      Pt_jets_X_tot += pTrack->Px();
      Pt_jets_Y_tot += pTrack->Py();  
      double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
      if(pDz < delta_z){
        Pt_jets_X += pTrack->Px();
        Pt_jets_Y += pTrack->Py();
      }
    }
  }

  if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
    return sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot);

  return 1.0;
}

//--------------------------------------------------------------------------------------------------
Double_t IsolationTools::BetaE(const TrackCol *tracks, const Electron *p, const Vertex *vertex, 
                               Double_t ptMin, Double_t  delta_z, Double_t extRadius,
                               Double_t intRadius){

  if(!tracks) return 1.0;
  if(tracks->GetEntries() <= 0) return 1.0;

  double Pt_jets_X = 0. ;
  double Pt_jets_Y = 0. ;
  double Pt_jets_X_tot = 0. ;
  double Pt_jets_Y_tot = 0. ;

  for(int i=0;i<int(tracks->GetEntries());i++){
    const Track *pTrack = tracks->At(i);

    if(pTrack && p->TrackerTrk() &&
       pTrack == p->TrackerTrk()) continue;

    if(pTrack && p->GsfTrk() &&
       pTrack == p->GsfTrk()) continue;

    if(pTrack->Pt() <= ptMin) continue;

    Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
    if ( dr < extRadius && dr >= intRadius ) {
      Pt_jets_X_tot += pTrack->Px();
      Pt_jets_Y_tot += pTrack->Py();  
      double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
      if(pDz < delta_z){
        Pt_jets_X += pTrack->Px();
        Pt_jets_Y += pTrack->Py();
      }
    }
  }

  if(sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot) > 0)
    return sqrt(Pt_jets_X*Pt_jets_X + Pt_jets_Y*Pt_jets_Y) / sqrt(Pt_jets_X_tot*Pt_jets_X_tot + Pt_jets_Y_tot*Pt_jets_Y_tot);

  return 1.0;
}
Revision:	1.9
Committed:	Mon Mar 7 12:46:02 2011 UTC (14 years, 1 month ago) by ceballos
Content type:	text/plain
Branch:	MAIN
CVS Tags:	Mit_020a, Mit_020, Mit_020pre1
Changes since 1.8:	+2 -1 lines
Log Message:	new stuff
#	User	Rev	Content
1	ceballos	1.9	// $Id: IsolationTools.cc,v 1.8 2011/02/21 13:50:20 ceballos Exp $
2	loizides	1.1
3			#include "MitPhysics/Utils/interface/IsolationTools.h"
4			#include "MitCommon/MathTools/interface/MathUtils.h"
5
6	loizides	1.4	ClassImp(mithep::IsolationTools)
7
8	loizides	1.1	using namespace mithep;
9
10			//--------------------------------------------------------------------------------------------------
11			Double_t IsolationTools::TrackIsolation(const Track *p, Double_t extRadius, Double_t intRadius,
12			Double_t ptLow, Double_t maxVtxZDist,
13	phedex	1.2	const Collection<Track> *tracks)
14	loizides	1.1	{
15			//Computes the Track Isolation: Summed Transverse Momentum of all tracks inside an
16			//annulus around the electron seed track.
17
18			Double_t ptSum =0.;
19			for (UInt_t i=0; i<tracks->GetEntries();i++) {
20			Double_t tmpPt = tracks->At(i)->Pt();
21			Double_t deltaZ = fabs(p->Z0() - tracks->At(i)->Z0());
22
23			//ignore the track if it is below the pt threshold
24			if (tmpPt < ptLow)
25			continue;
26			//ingore the track if it is too far away in Z
27			if (deltaZ > maxVtxZDist)
28			continue;
29
30	phedex	1.2	Double_t dr = MathUtils::DeltaR(p->Phi(),p->Eta(),tracks->At(i)->Phi(), tracks->At(i)->Eta());
31	loizides	1.1	//add the track pt if it is inside the annulus
32			if ( dr < extRadius &&
33			dr >= intRadius ) {
34			ptSum += tmpPt;
35			}
36			}
37			return ptSum;
38			}
39
40			//--------------------------------------------------------------------------------------------------
41			Double_t IsolationTools::EcalIsolation(const SuperCluster *sc, Double_t coneSize, Double_t etLow,
42	phedex	1.2	const Collection<BasicCluster> *basicClusters)
43	loizides	1.1	{
44			//Computes the Ecal Isolation: Summed Transverse Energy of all Basic Clusters inside a
45			//cone around the electron, excluding those that are inside the electron super cluster.
46
47			Double_t ecalIsol=0.;
48			const BasicCluster *basicCluster= 0;
49			for (UInt_t i=0; i<basicClusters->GetEntries();i++) {
50			basicCluster = basicClusters->At(i);
51			Double_t basicClusterEnergy = basicCluster->Energy();
52			Double_t basicClusterEta = basicCluster->Eta();
53			Double_t basicClusterEt = basicClusterEnergysin(2atan(exp(basicClusterEta)));
54
55	bendavid	1.3	if (basicClusterEt > etLow) {
56	loizides	1.1	bool inSuperCluster = false;
57
58			// loop over the basic clusters of the supercluster
59			// to make sure that the basic cluster is not inside
60			// the super cluster. We exclude those.
61			for (UInt_t j=0; j<sc->ClusterSize(); j++) {
62			const BasicCluster *tempBasicClusterInSuperCluster = sc->Cluster(j);
63			if (tempBasicClusterInSuperCluster == basicCluster) {
64			inSuperCluster = true;
65			}
66			}
67
68			if (!inSuperCluster) {
69	phedex	1.2	Double_t dr = MathUtils::DeltaR(sc->Phi(), sc->Eta(),
70			basicCluster->Phi(),basicCluster->Eta());
71	loizides	1.1	if(dr < coneSize) {
72			ecalIsol += basicClusterEt;
73			}
74			}
75			}
76			}
77			return ecalIsol;
78			}
79
80			//--------------------------------------------------------------------------------------------------
81			Double_t IsolationTools::CaloTowerHadIsolation(const ThreeVector *p, Double_t extRadius,
82			Double_t intRadius, Double_t etLow,
83	phedex	1.2	const Collection<CaloTower> *caloTowers)
84	loizides	1.1	{
85			//Computes the CaloTower Had Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers
86			//inside an annulus around the electron super cluster position.
87
88			Double_t sumEt = 0;
89			for (UInt_t i=0; i<caloTowers->GetEntries();i++) {
90			Double_t caloTowerEt = caloTowers->At(i)->HadEt();
91			Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
92			p->Phi(), p->Eta());
93			if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
94			sumEt += caloTowerEt;
95			}
96			}
97			return sumEt;
98			}
99
100			//--------------------------------------------------------------------------------------------------
101			Double_t IsolationTools::CaloTowerEmIsolation(const ThreeVector *p, Double_t extRadius,
102			Double_t intRadius, Double_t etLow,
103	phedex	1.2	const Collection<CaloTower> *caloTowers)
104	loizides	1.1	{
105			//Computes the CaloTower Em Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers
106			//inside an annulus around the electron super cluster position.
107
108			Double_t sumEt = 0;
109			for (UInt_t i=0; i<caloTowers->GetEntries();i++) {
110			Double_t caloTowerEt = caloTowers->At(i)->EmEt();
111			Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
112			p->Phi(), p->Eta());
113			if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
114			sumEt += caloTowerEt;
115			}
116			}
117			return sumEt;
118			}
119	ceballos	1.5
120			//--------------------------------------------------------------------------------------------------
121			Double_t IsolationTools::PFMuonIsolation(const Muon p, const Collection<PFCandidate> PFCands,
122	ceballos	1.8	const Vertex *vertex, Double_t delta_z, Double_t ptMin,
123			Double_t extRadius, Double_t intRadius, int isoType,
124			Double_t beta, const MuonCol *goodMuons,
125			const ElectronCol *goodElectrons)
126	ceballos	1.5	{
127			//Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
128			//annulus around the particle seed track.
129
130			Double_t zLepton = 0.0;
131	ceballos	1.8	if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
132	ceballos	1.5
133			Double_t ptSum =0.;
134			for (UInt_t i=0; i<PFCands->GetEntries();i++) {
135			const PFCandidate *pf = PFCands->At(i);
136
137			Bool_t isGoodType = kFALSE;
138			// all particles
139	ceballos	1.6	if (isoType == 0) isGoodType = kTRUE;
140	ceballos	1.5	// charged particles only
141	ceballos	1.6	else if(isoType == 1 && pf->BestTrk()) isGoodType = kTRUE;
142	ceballos	1.5	// charged particles and gammas only
143	ceballos	1.6	else if(isoType == 2 &&
144			(pf->BestTrk() \|\| pf->PFType() == PFCandidate::eGamma)) isGoodType = kTRUE;
145	ceballos	1.8	// all particles, rejecting good leptons
146			else if(isoType == 3) isGoodType = kTRUE;
147	ceballos	1.5
148			if(isGoodType == kFALSE) continue;
149
150	ceballos	1.6	if(pf->Pt() <= ptMin) continue;
151
152	ceballos	1.5	if(pf->TrackerTrk() && p->TrackerTrk() &&
153			pf->TrackerTrk() == p->TrackerTrk()) continue;
154
155			Double_t deltaZ = 0.0;
156			if(pf->BestTrk()) {
157	ceballos	1.8	deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
158	ceballos	1.5	}
159
160			// ignore the pf candidate if it is too far away in Z
161	ceballos	1.6	if (deltaZ >= delta_z)
162	ceballos	1.5	continue;
163
164			Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
165			// add the pf pt if it is inside the extRadius and outside the intRadius
166			if ( dr < extRadius &&
167			dr >= intRadius ) {
168	ceballos	1.8	Bool_t isLepton = kFALSE;
169			if(goodMuons && isoType == 3){
170			for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
171			const Muon *m = goodMuons->At(nl);
172			if(pf->TrackerTrk() && m->TrackerTrk() &&
173			pf->TrackerTrk() == m->TrackerTrk()) {
174			isLepton = kTRUE;
175			break;
176			}
177			}
178			}
179			if(goodElectrons && isLepton == kFALSE && isoType == 3){
180			for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
181			const Electron *e = goodElectrons->At(nl);
182			if(pf->TrackerTrk() && e->TrackerTrk() &&
183			pf->TrackerTrk() == e->TrackerTrk()) {
184			isLepton = kTRUE;
185			break;
186			}
187			if(pf->GsfTrk() && e->GsfTrk() &&
188			pf->GsfTrk() == e->GsfTrk()) {
189			isLepton = kTRUE;
190			break;
191			}
192			}
193			}
194			if(isLepton == kFALSE){
195			if(pf->BestTrk()) ptSum += pf->Pt();
196			else ptSum += pf->Pt()*beta;
197			}
198	ceballos	1.5	}
199			}
200			return ptSum;
201			}
202			//--------------------------------------------------------------------------------------------------
203			Double_t IsolationTools::PFElectronIsolation(const Electron p, const PFCandidateCol PFCands,
204	ceballos	1.8	const Vertex *vertex, Double_t delta_z, Double_t ptMin,
205			Double_t extRadius, Double_t intRadius, int isoType,
206			Double_t beta, const MuonCol *goodMuons,
207			const ElectronCol *goodElectrons)
208	ceballos	1.5	{
209			//Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
210			//annulus around the particle seed track.
211
212			Double_t zLepton = 0.0;
213	ceballos	1.8	if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
214	ceballos	1.5
215			Double_t ptSum =0.;
216			for (UInt_t i=0; i<PFCands->GetEntries();i++) {
217			const PFCandidate *pf = PFCands->At(i);
218
219			Bool_t isGoodType = kFALSE;
220			// all particles
221	ceballos	1.6	if (isoType == 0) isGoodType = kTRUE;
222	ceballos	1.5	// charged particles only
223	ceballos	1.6	else if(isoType == 1 && pf->BestTrk()) isGoodType = kTRUE;
224	ceballos	1.5	// charged particles and gammas only
225	ceballos	1.6	else if(isoType == 2 &&
226			(pf->BestTrk() \|\| pf->PFType() == PFCandidate::eGamma)) isGoodType = kTRUE;
227	ceballos	1.8	// all particles, rejecting good leptons
228			else if(isoType == 3) isGoodType = kTRUE;
229	ceballos	1.5
230			if(isGoodType == kFALSE) continue;
231
232	ceballos	1.6	if(pf->Pt() <= ptMin) continue;
233
234	ceballos	1.5	if(pf->TrackerTrk() && p->TrackerTrk() &&
235			pf->TrackerTrk() == p->TrackerTrk()) continue;
236
237			if(pf->GsfTrk() && p->GsfTrk() &&
238			pf->GsfTrk() == p->GsfTrk()) continue;
239
240			Double_t deltaZ = 0.0;
241			if(pf->BestTrk()) {
242	ceballos	1.8	deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
243	ceballos	1.5	}
244
245			// ignore the pf candidate if it is too far away in Z
246	ceballos	1.6	if (deltaZ >= delta_z)
247	ceballos	1.5	continue;
248
249			Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
250			// add the pf pt if it is inside the extRadius and outside the intRadius
251			if ( dr < extRadius &&
252			dr >= intRadius ) {
253	ceballos	1.8	Bool_t isLepton = kFALSE;
254			if(goodMuons && isoType == 3){
255			for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
256			const Muon *m = goodMuons->At(nl);
257			if(pf->TrackerTrk() && m->TrackerTrk() &&
258			pf->TrackerTrk() == m->TrackerTrk()) {
259			isLepton = kTRUE;
260			break;
261			}
262			}
263			}
264			if(goodElectrons && isLepton == kFALSE && isoType == 3){
265			for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
266			const Electron *e = goodElectrons->At(nl);
267			if(pf->TrackerTrk() && e->TrackerTrk() &&
268			pf->TrackerTrk() == e->TrackerTrk()) {
269			isLepton = kTRUE;
270			break;
271			}
272			if(pf->GsfTrk() && e->GsfTrk() &&
273			pf->GsfTrk() == e->GsfTrk()) {
274			isLepton = kTRUE;
275			break;
276			}
277			}
278			}
279			if(isLepton == kFALSE){
280			if(pf->BestTrk()) ptSum += pf->Pt();
281			else ptSum += pf->Pt()*beta;
282			}
283	ceballos	1.5	}
284			}
285			return ptSum;
286			}
287	ceballos	1.7	//--------------------------------------------------------------------------------------------------
288	ceballos	1.8	Double_t IsolationTools::BetaM(const TrackCol tracks, const Muon p, const Vertex *vertex,
289	ceballos	1.7	Double_t ptMin, Double_t delta_z, Double_t extRadius,
290			Double_t intRadius){
291
292	ceballos	1.9	if(!tracks) return 1.0;
293	ceballos	1.7	if(tracks->GetEntries() <= 0) return 1.0;
294
295			double Pt_jets_X = 0. ;
296			double Pt_jets_Y = 0. ;
297			double Pt_jets_X_tot = 0. ;
298			double Pt_jets_Y_tot = 0. ;
299
300			for(int i=0;i<int(tracks->GetEntries());i++){
301			const Track *pTrack = tracks->At(i);
302
303			if(pTrack && p->TrackerTrk() &&
304			pTrack == p->TrackerTrk()) continue;
305
306			if(pTrack->Pt() <= ptMin) continue;
307
308			Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
309			if ( dr < extRadius && dr >= intRadius ) {
310			Pt_jets_X_tot += pTrack->Px();
311			Pt_jets_Y_tot += pTrack->Py();
312	ceballos	1.8	double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
313	ceballos	1.7	if(pDz < delta_z){
314			Pt_jets_X += pTrack->Px();
315			Pt_jets_Y += pTrack->Py();
316			}
317			}
318			}
319
320			if(sqrt(Pt_jets_X_totPt_jets_X_tot + Pt_jets_Y_totPt_jets_Y_tot) > 0)
321			return sqrt(Pt_jets_XPt_jets_X + Pt_jets_YPt_jets_Y) / sqrt(Pt_jets_X_totPt_jets_X_tot + Pt_jets_Y_totPt_jets_Y_tot);
322
323			return 1.0;
324			}
325
326			//--------------------------------------------------------------------------------------------------
327	ceballos	1.8	Double_t IsolationTools::BetaE(const TrackCol tracks, const Electron p, const Vertex *vertex,
328	ceballos	1.7	Double_t ptMin, Double_t delta_z, Double_t extRadius,
329			Double_t intRadius){
330
331			if(!tracks) return 1.0;
332			if(tracks->GetEntries() <= 0) return 1.0;
333
334			double Pt_jets_X = 0. ;
335			double Pt_jets_Y = 0. ;
336			double Pt_jets_X_tot = 0. ;
337			double Pt_jets_Y_tot = 0. ;
338
339			for(int i=0;i<int(tracks->GetEntries());i++){
340			const Track *pTrack = tracks->At(i);
341
342			if(pTrack && p->TrackerTrk() &&
343			pTrack == p->TrackerTrk()) continue;
344
345			if(pTrack && p->GsfTrk() &&
346			pTrack == p->GsfTrk()) continue;
347
348			if(pTrack->Pt() <= ptMin) continue;
349
350			Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
351			if ( dr < extRadius && dr >= intRadius ) {
352			Pt_jets_X_tot += pTrack->Px();
353			Pt_jets_Y_tot += pTrack->Py();
354	ceballos	1.8	double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
355	ceballos	1.7	if(pDz < delta_z){
356			Pt_jets_X += pTrack->Px();
357			Pt_jets_Y += pTrack->Py();
358			}
359			}
360			}
361
362			if(sqrt(Pt_jets_X_totPt_jets_X_tot + Pt_jets_Y_totPt_jets_Y_tot) > 0)
363			return sqrt(Pt_jets_XPt_jets_X + Pt_jets_YPt_jets_Y) / sqrt(Pt_jets_X_totPt_jets_X_tot + Pt_jets_Y_totPt_jets_Y_tot);
364
365			return 1.0;
366			}