Utils/src/IsolationTools.cc

// $Id: IsolationTools.cc,v 1.15 2011/05/12 12:26:42 mzanetti Exp $

#include "MitPhysics/Utils/interface/IsolationTools.h"
#include "MitPhysics/Utils/interface/PhotonTools.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)
{
  //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);
    
    // pt cut applied to neutrals
    if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;

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

    // ignore the pf candidate if it is too far away in Z
    Double_t deltaZ = 0.0;
    if(pf->HasTrk()) {
      deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
    }
    if (deltaZ >= delta_z) 
      continue;

    // add the pf pt if it is inside the extRadius and outside the intRadius           
    Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
    if ( dr < extRadius && dr >= intRadius ) ptSum += pf->Pt();
  
  }
  return ptSum;
}

//--------------------------------------------------------------------------------------------------
Double_t IsolationTools::PFMuonIsolation(const Muon *p, const Collection<PFCandidate> *PFCands, const Vertex *vertex, 
                                         const MuonCol *goodMuons, const ElectronCol *goodElectrons,
                                         Double_t  delta_z, Double_t ptMin, Double_t extRadius,
                                         Double_t intRadius, int isoType, Double_t beta)
{
  //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;

    // 0.1 pt cut applied to charged
    if( pf->BestTrk() && pf->Pt() <= 0.1)   continue;
    // pt cut applied to neutrals
    if(!pf->HasTrk() && 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) 
{

  //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);
    
    // 0.1 pt cut applied to charged
    if( pf->BestTrk() && pf->Pt() <= 0.1)   continue;

    // pt cut applied to neutrals
    if(!pf->HasTrk() && 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 ) {

      //EtaStrip Veto for Gamma 
      if (pf->PFType() == PFCandidate::eGamma && fabs(p->Eta() - pf->Eta()) < 0.025) continue;

      //InnerCone (One Tower = dR < 0.07) Veto for non-gamma neutrals
      if (!pf->HasTrk() && pf->PFType() == PFCandidate::eNeutralHadron
          && MathUtils::DeltaR(p->Mom(), pf->Mom()) < 0.07 ) continue; 

      ptSum += pf->Pt();            

    }
  }
  return ptSum;
}


//--------------------------------------------------------------------------------------------------
Double_t IsolationTools::PFElectronIsolation(const Electron *p, const PFCandidateCol *PFCands, 
                                             const Vertex *vertex, const MuonCol *goodMuons, 
                                             const ElectronCol *goodElectrons, 
                                             Double_t delta_z, Double_t ptMin,
                                             Double_t extRadius, Double_t intRadius, int isoType,
                                             Double_t beta)
{
  //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;


    // 0.1 pt cut applied to charged
    if( pf->BestTrk() && pf->Pt() <= 0.1)   continue;

    // pt cut applied to neutrals
    if(!pf->HasTrk() && 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 == kTRUE) continue;

      //EtaStrip Veto for Gamma 
      if (pf->PFType() == PFCandidate::eGamma && fabs(p->Eta() - pf->Eta()) < 0.025) continue;

      //InnerCone (One Tower = dR < 0.07) Veto for non-gamma neutrals
      if (!pf->HasTrk() && pf->PFType() == PFCandidate::eNeutralHadron
          && MathUtils::DeltaR(p->Mom(), pf->Mom()) < 0.07 ) continue; 


      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;
}


// method added by F.Stoeckli: computes the track isolation with NO constrint on the OV-track compatibility
Double_t IsolationTools::TrackIsolationNoPV(const mithep::Particle* p, const BaseVertex* bsp, 
                                            Double_t extRadius, 
                                            Double_t intRadius, 
                                            Double_t ptLow, 
                                            Double_t etaStrip,
                                            Double_t maxD0,
                                            mithep::TrackQuality::EQuality quality,
                                            const mithep::Collection<mithep::Track> *tracks,
                                            UInt_t maxNExpectedHitsInner,
                                            const mithep::DecayParticleCol *conversions) {
  
  // loop over all tracks
  Double_t tPt = 0.;
  for(UInt_t i=0; i<tracks->GetEntries(); ++i) {
    const Track* t = tracks->At(i);
    if ( t->Pt() < ptLow ) continue;
    if ( ! t->Quality().Quality(quality) ) continue;
    // only check for beamspot if available, otherwise ignore cut
    if ( bsp && fabs(t->D0Corrected( *bsp) ) > maxD0) continue;
    if (t->NExpectedHitsInner()>maxNExpectedHitsInner) continue;
    if (conversions && PhotonTools::MatchedConversion(t,conversions,bsp)) continue;
    Double_t dR   = MathUtils::DeltaR(t->Mom(),p->Mom());
    Double_t dEta = fabs(t->Eta()-p->Eta());
    if(dR < extRadius && dR > intRadius && dEta > etaStrip) tPt += t->Pt();
  }
  return tPt;
}

Revision:	1.16
Committed:	Thu May 12 16:08:36 2011 UTC (13 years, 11 months ago) by sixie
Content type:	text/plain
Branch:	MAIN
Changes since 1.15:	+61 -3 lines
Log Message:	modify pfiso for ele
#	User	Rev	Content
1	sixie	1.16	// $Id: IsolationTools.cc,v 1.15 2011/05/12 12:26:42 mzanetti Exp $
2	loizides	1.1
3			#include "MitPhysics/Utils/interface/IsolationTools.h"
4	bendavid	1.11	#include "MitPhysics/Utils/interface/PhotonTools.h"
5	loizides	1.1	#include "MitCommon/MathTools/interface/MathUtils.h"
6
7	loizides	1.4	ClassImp(mithep::IsolationTools)
8
9	loizides	1.1	using namespace mithep;
10
11			//--------------------------------------------------------------------------------------------------
12			Double_t IsolationTools::TrackIsolation(const Track *p, Double_t extRadius, Double_t intRadius,
13			Double_t ptLow, Double_t maxVtxZDist,
14	phedex	1.2	const Collection<Track> *tracks)
15	loizides	1.1	{
16			//Computes the Track Isolation: Summed Transverse Momentum of all tracks inside an
17			//annulus around the electron seed track.
18
19			Double_t ptSum =0.;
20			for (UInt_t i=0; i<tracks->GetEntries();i++) {
21			Double_t tmpPt = tracks->At(i)->Pt();
22			Double_t deltaZ = fabs(p->Z0() - tracks->At(i)->Z0());
23
24			//ignore the track if it is below the pt threshold
25			if (tmpPt < ptLow)
26			continue;
27			//ingore the track if it is too far away in Z
28			if (deltaZ > maxVtxZDist)
29			continue;
30
31	phedex	1.2	Double_t dr = MathUtils::DeltaR(p->Phi(),p->Eta(),tracks->At(i)->Phi(), tracks->At(i)->Eta());
32	loizides	1.1	//add the track pt if it is inside the annulus
33			if ( dr < extRadius &&
34			dr >= intRadius ) {
35			ptSum += tmpPt;
36			}
37			}
38			return ptSum;
39			}
40
41			//--------------------------------------------------------------------------------------------------
42			Double_t IsolationTools::EcalIsolation(const SuperCluster *sc, Double_t coneSize, Double_t etLow,
43	phedex	1.2	const Collection<BasicCluster> *basicClusters)
44	loizides	1.1	{
45			//Computes the Ecal Isolation: Summed Transverse Energy of all Basic Clusters inside a
46			//cone around the electron, excluding those that are inside the electron super cluster.
47
48			Double_t ecalIsol=0.;
49			const BasicCluster *basicCluster= 0;
50			for (UInt_t i=0; i<basicClusters->GetEntries();i++) {
51			basicCluster = basicClusters->At(i);
52			Double_t basicClusterEnergy = basicCluster->Energy();
53			Double_t basicClusterEta = basicCluster->Eta();
54			Double_t basicClusterEt = basicClusterEnergysin(2atan(exp(basicClusterEta)));
55
56	bendavid	1.3	if (basicClusterEt > etLow) {
57	loizides	1.1	bool inSuperCluster = false;
58
59			// loop over the basic clusters of the supercluster
60			// to make sure that the basic cluster is not inside
61			// the super cluster. We exclude those.
62			for (UInt_t j=0; j<sc->ClusterSize(); j++) {
63			const BasicCluster *tempBasicClusterInSuperCluster = sc->Cluster(j);
64			if (tempBasicClusterInSuperCluster == basicCluster) {
65			inSuperCluster = true;
66			}
67			}
68
69			if (!inSuperCluster) {
70	phedex	1.2	Double_t dr = MathUtils::DeltaR(sc->Phi(), sc->Eta(),
71			basicCluster->Phi(),basicCluster->Eta());
72	loizides	1.1	if(dr < coneSize) {
73			ecalIsol += basicClusterEt;
74			}
75			}
76			}
77			}
78			return ecalIsol;
79			}
80
81			//--------------------------------------------------------------------------------------------------
82			Double_t IsolationTools::CaloTowerHadIsolation(const ThreeVector *p, Double_t extRadius,
83			Double_t intRadius, Double_t etLow,
84	phedex	1.2	const Collection<CaloTower> *caloTowers)
85	loizides	1.1	{
86			//Computes the CaloTower Had Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers
87			//inside an annulus around the electron super cluster position.
88
89			Double_t sumEt = 0;
90			for (UInt_t i=0; i<caloTowers->GetEntries();i++) {
91			Double_t caloTowerEt = caloTowers->At(i)->HadEt();
92			Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
93			p->Phi(), p->Eta());
94			if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
95			sumEt += caloTowerEt;
96			}
97			}
98			return sumEt;
99			}
100
101			//--------------------------------------------------------------------------------------------------
102			Double_t IsolationTools::CaloTowerEmIsolation(const ThreeVector *p, Double_t extRadius,
103			Double_t intRadius, Double_t etLow,
104	phedex	1.2	const Collection<CaloTower> *caloTowers)
105	loizides	1.1	{
106			//Computes the CaloTower Em Et Isolation: Summed Hadronic Transverse Energy of all Calo Towers
107			//inside an annulus around the electron super cluster position.
108
109			Double_t sumEt = 0;
110			for (UInt_t i=0; i<caloTowers->GetEntries();i++) {
111			Double_t caloTowerEt = caloTowers->At(i)->EmEt();
112			Double_t dr = MathUtils::DeltaR(caloTowers->At(i)->Phi(), caloTowers->At(i)->Eta(),
113			p->Phi(), p->Eta());
114			if (dr < extRadius && dr > intRadius && caloTowerEt > etLow) {
115			sumEt += caloTowerEt;
116			}
117			}
118			return sumEt;
119			}
120	ceballos	1.5
121			//--------------------------------------------------------------------------------------------------
122			Double_t IsolationTools::PFMuonIsolation(const Muon p, const Collection<PFCandidate> PFCands,
123	ceballos	1.8	const Vertex *vertex, Double_t delta_z, Double_t ptMin,
124	mzanetti	1.15	Double_t extRadius, Double_t intRadius)
125			{
126			//Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
127			//annulus around the particle seed track.
128
129			Double_t zLepton = 0.0;
130			if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
131
132			Double_t ptSum =0.;
133			for (UInt_t i=0; i<PFCands->GetEntries();i++) {
134			const PFCandidate *pf = PFCands->At(i);
135
136			// pt cut applied to neutrals
137			if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
138
139			// exclude muon
140			if(pf->TrackerTrk() && p->TrackerTrk() &&
141			pf->TrackerTrk() == p->TrackerTrk()) continue;
142
143			// ignore the pf candidate if it is too far away in Z
144			Double_t deltaZ = 0.0;
145			if(pf->HasTrk()) {
146			deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
147			}
148			if (deltaZ >= delta_z)
149			continue;
150
151			// add the pf pt if it is inside the extRadius and outside the intRadius
152			Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
153			if ( dr < extRadius && dr >= intRadius ) ptSum += pf->Pt();
154
155			}
156			return ptSum;
157			}
158
159			//--------------------------------------------------------------------------------------------------
160			Double_t IsolationTools::PFMuonIsolation(const Muon p, const Collection<PFCandidate> PFCands, const Vertex *vertex,
161			const MuonCol goodMuons, const ElectronCol goodElectrons,
162			Double_t delta_z, Double_t ptMin, Double_t extRadius,
163			Double_t intRadius, int isoType, Double_t beta)
164	ceballos	1.5	{
165			//Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
166			//annulus around the particle seed track.
167
168			Double_t zLepton = 0.0;
169	ceballos	1.8	if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
170	ceballos	1.5
171			Double_t ptSum =0.;
172			for (UInt_t i=0; i<PFCands->GetEntries();i++) {
173			const PFCandidate *pf = PFCands->At(i);
174
175			Bool_t isGoodType = kFALSE;
176			// all particles
177	ceballos	1.6	if (isoType == 0) isGoodType = kTRUE;
178	ceballos	1.5	// charged particles only
179	ceballos	1.6	else if(isoType == 1 && pf->BestTrk()) isGoodType = kTRUE;
180	ceballos	1.5	// charged particles and gammas only
181	ceballos	1.6	else if(isoType == 2 &&
182			(pf->BestTrk() \|\| pf->PFType() == PFCandidate::eGamma)) isGoodType = kTRUE;
183	ceballos	1.8	// all particles, rejecting good leptons
184			else if(isoType == 3) isGoodType = kTRUE;
185	ceballos	1.5
186			if(isGoodType == kFALSE) continue;
187
188	mzanetti	1.13	// 0.1 pt cut applied to charged
189	ceballos	1.12	if( pf->BestTrk() && pf->Pt() <= 0.1) continue;
190	mzanetti	1.13	// pt cut applied to neutrals
191			if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
192	ceballos	1.6
193	ceballos	1.5	if(pf->TrackerTrk() && p->TrackerTrk() &&
194			pf->TrackerTrk() == p->TrackerTrk()) continue;
195
196			Double_t deltaZ = 0.0;
197			if(pf->BestTrk()) {
198	ceballos	1.8	deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
199	ceballos	1.5	}
200
201			// ignore the pf candidate if it is too far away in Z
202	ceballos	1.6	if (deltaZ >= delta_z)
203	ceballos	1.5	continue;
204
205			Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
206			// add the pf pt if it is inside the extRadius and outside the intRadius
207			if ( dr < extRadius &&
208			dr >= intRadius ) {
209	ceballos	1.8	Bool_t isLepton = kFALSE;
210			if(goodMuons && isoType == 3){
211			for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
212			const Muon *m = goodMuons->At(nl);
213			if(pf->TrackerTrk() && m->TrackerTrk() &&
214			pf->TrackerTrk() == m->TrackerTrk()) {
215			isLepton = kTRUE;
216			break;
217			}
218			}
219			}
220			if(goodElectrons && isLepton == kFALSE && isoType == 3){
221			for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
222			const Electron *e = goodElectrons->At(nl);
223			if(pf->TrackerTrk() && e->TrackerTrk() &&
224			pf->TrackerTrk() == e->TrackerTrk()) {
225			isLepton = kTRUE;
226			break;
227			}
228			if(pf->GsfTrk() && e->GsfTrk() &&
229			pf->GsfTrk() == e->GsfTrk()) {
230			isLepton = kTRUE;
231			break;
232			}
233			}
234			}
235			if(isLepton == kFALSE){
236			if(pf->BestTrk()) ptSum += pf->Pt();
237			else ptSum += pf->Pt()*beta;
238			}
239	ceballos	1.5	}
240			}
241			return ptSum;
242			}
243	mzanetti	1.15
244	ceballos	1.5	//--------------------------------------------------------------------------------------------------
245			Double_t IsolationTools::PFElectronIsolation(const Electron p, const PFCandidateCol PFCands,
246	ceballos	1.8	const Vertex *vertex, Double_t delta_z, Double_t ptMin,
247	sixie	1.16	Double_t extRadius, Double_t intRadius)
248			{
249
250			//Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
251			//annulus around the particle seed track.
252
253			Double_t zLepton = 0.0;
254			if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
255
256			Double_t ptSum =0.;
257			for (UInt_t i=0; i<PFCands->GetEntries();i++) {
258			const PFCandidate *pf = PFCands->At(i);
259
260			// 0.1 pt cut applied to charged
261			if( pf->BestTrk() && pf->Pt() <= 0.1) continue;
262
263			// pt cut applied to neutrals
264			if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
265
266			if(pf->TrackerTrk() && p->TrackerTrk() &&
267			pf->TrackerTrk() == p->TrackerTrk()) continue;
268
269			if(pf->GsfTrk() && p->GsfTrk() &&
270			pf->GsfTrk() == p->GsfTrk()) continue;
271
272			Double_t deltaZ = 0.0;
273			if(pf->BestTrk()) {
274			deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
275			}
276
277			// ignore the pf candidate if it is too far away in Z
278			if (deltaZ >= delta_z)
279			continue;
280
281			Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
282			// add the pf pt if it is inside the extRadius and outside the intRadius
283			if ( dr < extRadius &&
284			dr >= intRadius ) {
285
286			//EtaStrip Veto for Gamma
287			if (pf->PFType() == PFCandidate::eGamma && fabs(p->Eta() - pf->Eta()) < 0.025) continue;
288
289			//InnerCone (One Tower = dR < 0.07) Veto for non-gamma neutrals
290			if (!pf->HasTrk() && pf->PFType() == PFCandidate::eNeutralHadron
291			&& MathUtils::DeltaR(p->Mom(), pf->Mom()) < 0.07 ) continue;
292
293			ptSum += pf->Pt();
294
295			}
296			}
297			return ptSum;
298			}
299
300
301			//--------------------------------------------------------------------------------------------------
302			Double_t IsolationTools::PFElectronIsolation(const Electron p, const PFCandidateCol PFCands,
303			const Vertex vertex, const MuonCol goodMuons,
304			const ElectronCol *goodElectrons,
305			Double_t delta_z, Double_t ptMin,
306	ceballos	1.8	Double_t extRadius, Double_t intRadius, int isoType,
307	sixie	1.16	Double_t beta)
308	ceballos	1.5	{
309			//Computes the PF Isolation: Summed Transverse Momentum of all PF candidates inside an
310			//annulus around the particle seed track.
311
312			Double_t zLepton = 0.0;
313	ceballos	1.8	if(p->BestTrk()) zLepton = p->BestTrk()->DzCorrected(*vertex);
314	ceballos	1.5
315			Double_t ptSum =0.;
316			for (UInt_t i=0; i<PFCands->GetEntries();i++) {
317			const PFCandidate *pf = PFCands->At(i);
318
319			Bool_t isGoodType = kFALSE;
320			// all particles
321	ceballos	1.6	if (isoType == 0) isGoodType = kTRUE;
322	ceballos	1.5	// charged particles only
323	ceballos	1.6	else if(isoType == 1 && pf->BestTrk()) isGoodType = kTRUE;
324	ceballos	1.5	// charged particles and gammas only
325	ceballos	1.6	else if(isoType == 2 &&
326			(pf->BestTrk() \|\| pf->PFType() == PFCandidate::eGamma)) isGoodType = kTRUE;
327	ceballos	1.8	// all particles, rejecting good leptons
328			else if(isoType == 3) isGoodType = kTRUE;
329	ceballos	1.5
330			if(isGoodType == kFALSE) continue;
331
332	sixie	1.14
333			// 0.1 pt cut applied to charged
334	ceballos	1.12	if( pf->BestTrk() && pf->Pt() <= 0.1) continue;
335
336	sixie	1.14	// pt cut applied to neutrals
337			if(!pf->HasTrk() && pf->Pt() <= ptMin) continue;
338	ceballos	1.6
339	ceballos	1.5	if(pf->TrackerTrk() && p->TrackerTrk() &&
340			pf->TrackerTrk() == p->TrackerTrk()) continue;
341
342			if(pf->GsfTrk() && p->GsfTrk() &&
343			pf->GsfTrk() == p->GsfTrk()) continue;
344
345			Double_t deltaZ = 0.0;
346			if(pf->BestTrk()) {
347	ceballos	1.8	deltaZ = TMath::Abs(pf->BestTrk()->DzCorrected(*vertex) - zLepton);
348	ceballos	1.5	}
349
350			// ignore the pf candidate if it is too far away in Z
351	ceballos	1.6	if (deltaZ >= delta_z)
352	ceballos	1.5	continue;
353
354			Double_t dr = MathUtils::DeltaR(p->Mom(), pf->Mom());
355			// add the pf pt if it is inside the extRadius and outside the intRadius
356			if ( dr < extRadius &&
357			dr >= intRadius ) {
358	ceballos	1.8	Bool_t isLepton = kFALSE;
359			if(goodMuons && isoType == 3){
360			for (UInt_t nl=0; nl<goodMuons->GetEntries();nl++) {
361			const Muon *m = goodMuons->At(nl);
362			if(pf->TrackerTrk() && m->TrackerTrk() &&
363			pf->TrackerTrk() == m->TrackerTrk()) {
364			isLepton = kTRUE;
365			break;
366			}
367			}
368			}
369			if(goodElectrons && isLepton == kFALSE && isoType == 3){
370			for (UInt_t nl=0; nl<goodElectrons->GetEntries();nl++) {
371			const Electron *e = goodElectrons->At(nl);
372			if(pf->TrackerTrk() && e->TrackerTrk() &&
373			pf->TrackerTrk() == e->TrackerTrk()) {
374			isLepton = kTRUE;
375			break;
376			}
377			if(pf->GsfTrk() && e->GsfTrk() &&
378			pf->GsfTrk() == e->GsfTrk()) {
379			isLepton = kTRUE;
380			break;
381			}
382			}
383			}
384	sixie	1.14
385			if (isLepton == kTRUE) continue;
386
387			//EtaStrip Veto for Gamma
388			if (pf->PFType() == PFCandidate::eGamma && fabs(p->Eta() - pf->Eta()) < 0.025) continue;
389
390			//InnerCone (One Tower = dR < 0.07) Veto for non-gamma neutrals
391			if (!pf->HasTrk() && pf->PFType() == PFCandidate::eNeutralHadron
392			&& MathUtils::DeltaR(p->Mom(), pf->Mom()) < 0.07 ) continue;
393
394
395			if(pf->BestTrk()) ptSum += pf->Pt();
396			else ptSum += pf->Pt()*beta;
397
398
399	ceballos	1.5	}
400			}
401			return ptSum;
402			}
403	ceballos	1.7	//--------------------------------------------------------------------------------------------------
404	ceballos	1.8	Double_t IsolationTools::BetaM(const TrackCol tracks, const Muon p, const Vertex *vertex,
405	ceballos	1.7	Double_t ptMin, Double_t delta_z, Double_t extRadius,
406			Double_t intRadius){
407
408	ceballos	1.9	if(!tracks) return 1.0;
409	ceballos	1.7	if(tracks->GetEntries() <= 0) return 1.0;
410
411			double Pt_jets_X = 0. ;
412			double Pt_jets_Y = 0. ;
413			double Pt_jets_X_tot = 0. ;
414			double Pt_jets_Y_tot = 0. ;
415
416			for(int i=0;i<int(tracks->GetEntries());i++){
417			const Track *pTrack = tracks->At(i);
418
419			if(pTrack && p->TrackerTrk() &&
420			pTrack == p->TrackerTrk()) continue;
421
422			if(pTrack->Pt() <= ptMin) continue;
423
424			Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
425			if ( dr < extRadius && dr >= intRadius ) {
426			Pt_jets_X_tot += pTrack->Px();
427			Pt_jets_Y_tot += pTrack->Py();
428	ceballos	1.8	double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
429	ceballos	1.7	if(pDz < delta_z){
430			Pt_jets_X += pTrack->Px();
431			Pt_jets_Y += pTrack->Py();
432			}
433			}
434			}
435
436			if(sqrt(Pt_jets_X_totPt_jets_X_tot + Pt_jets_Y_totPt_jets_Y_tot) > 0)
437			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);
438
439			return 1.0;
440			}
441
442			//--------------------------------------------------------------------------------------------------
443	ceballos	1.8	Double_t IsolationTools::BetaE(const TrackCol tracks, const Electron p, const Vertex *vertex,
444	ceballos	1.7	Double_t ptMin, Double_t delta_z, Double_t extRadius,
445			Double_t intRadius){
446
447			if(!tracks) return 1.0;
448			if(tracks->GetEntries() <= 0) return 1.0;
449
450			double Pt_jets_X = 0. ;
451			double Pt_jets_Y = 0. ;
452			double Pt_jets_X_tot = 0. ;
453			double Pt_jets_Y_tot = 0. ;
454
455			for(int i=0;i<int(tracks->GetEntries());i++){
456			const Track *pTrack = tracks->At(i);
457
458			if(pTrack && p->TrackerTrk() &&
459			pTrack == p->TrackerTrk()) continue;
460
461			if(pTrack && p->GsfTrk() &&
462			pTrack == p->GsfTrk()) continue;
463
464			if(pTrack->Pt() <= ptMin) continue;
465
466			Double_t dr = MathUtils::DeltaR(pTrack->Mom(),p->Mom());
467			if ( dr < extRadius && dr >= intRadius ) {
468			Pt_jets_X_tot += pTrack->Px();
469			Pt_jets_Y_tot += pTrack->Py();
470	ceballos	1.8	double pDz = TMath::Abs(pTrack->DzCorrected(*vertex));
471	ceballos	1.7	if(pDz < delta_z){
472			Pt_jets_X += pTrack->Px();
473			Pt_jets_Y += pTrack->Py();
474			}
475			}
476			}
477
478			if(sqrt(Pt_jets_X_totPt_jets_X_tot + Pt_jets_Y_totPt_jets_Y_tot) > 0)
479			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);
480
481			return 1.0;
482			}
483	fabstoec	1.10
484
485			// method added by F.Stoeckli: computes the track isolation with NO constrint on the OV-track compatibility
486			Double_t IsolationTools::TrackIsolationNoPV(const mithep::Particle* p, const BaseVertex* bsp,
487			Double_t extRadius,
488			Double_t intRadius,
489			Double_t ptLow,
490			Double_t etaStrip,
491			Double_t maxD0,
492			mithep::TrackQuality::EQuality quality,
493	bendavid	1.11	const mithep::Collection<mithep::Track> *tracks,
494			UInt_t maxNExpectedHitsInner,
495			const mithep::DecayParticleCol *conversions) {
496	fabstoec	1.10
497			// loop over all tracks
498			Double_t tPt = 0.;
499			for(UInt_t i=0; i<tracks->GetEntries(); ++i) {
500			const Track* t = tracks->At(i);
501			if ( t->Pt() < ptLow ) continue;
502			if ( ! t->Quality().Quality(quality) ) continue;
503			// only check for beamspot if available, otherwise ignore cut
504			if ( bsp && fabs(t->D0Corrected( *bsp) ) > maxD0) continue;
505	bendavid	1.11	if (t->NExpectedHitsInner()>maxNExpectedHitsInner) continue;
506			if (conversions && PhotonTools::MatchedConversion(t,conversions,bsp)) continue;
507	fabstoec	1.10	Double_t dR = MathUtils::DeltaR(t->Mom(),p->Mom());
508			Double_t dEta = fabs(t->Eta()-p->Eta());
509			if(dR < extRadius && dR > intRadius && dEta > etaStrip) tPt += t->Pt();
510			}
511			return tPt;
512			}
513