egamma Namespace Reference

Namespaces
	conv

Classes
struct	ElectronMomentum

Enumerations
enum	HoeMode { HoeMode::SingleTower = 0, HoeMode::TowersBehindCluster = 1 }

Functions
float	classBasedElectronEnergy (reco::GsfElectron const &, reco::BeamSpot const &, EcalClusterFunctionBaseClass const &crackCorrectionFunction)
double	classBasedElectronEnergyUncertainty (reco::GsfElectron const &)
reco::GsfElectron::Classification	classifyElectron (reco::GsfElectron const &)
ElectronMomentum	correctElectronMomentum (reco::GsfElectron const &, TrajectoryStateOnSurface const &)
double	depth1HcalESum (std::vector< CaloTowerDetId > const &towers, CaloTowerCollection const &)
double	depth2HcalESum (std::vector< CaloTowerDetId > const &towers, CaloTowerCollection const &)
float	ecalClusterEnergyUncertaintyElectronSpecific (reco::SuperCluster const &superCluster)
double	electronEnergyUncertainty (reco::GsfElectron::Classification c, double eta, double brem, double energy)
std::pair< reco::TrackRef, float >	getClosestCtfToGsf (reco::GsfTrackRef const &, edm::Handle< reco::TrackCollection > const &ctfTracksH, edm::soa::EtaPhiTableView trackEtaPhiTable)
template<typename T >
uint32_t	getRandomSeedFromObj (const edm::Event &iEvent, const T &obj, size_t nrObjs, size_t objNr)
uint32_t	getRandomSeedFromSC (const edm::Event &iEvent, const reco::SuperClusterRef scRef)
bool	hasActiveHcal (std::vector< CaloTowerDetId > const &towers, CaloTowerConstituentsMap const &towerMap, HcalChannelQuality const &hcalQuality, HcalTopology const &hcalTopology)
bool	isBetterElectron (reco::GsfElectron const &, reco::GsfElectron const &)
bool	isInnermostElectron (reco::GsfElectron const &, reco::GsfElectron const &)
int	sharedDets (reco::GsfTrackRef const &, reco::GsfTrackRef const &)
float	sharedEnergy (reco::CaloCluster const &clu1, reco::CaloCluster const &clu2, EcalRecHitCollection const &barrelRecHits, EcalRecHitCollection const &endcapRecHits)
float	sharedEnergy (reco::SuperClusterRef const &sc1, reco::SuperClusterRef const &sc2, EcalRecHitCollection const &barrelRecHits, EcalRecHitCollection const &endcapRecHits)
int	sharedHits (reco::GsfTrackRef const &, reco::GsfTrackRef const &)
double	simpleElectronEnergyUncertainty (reco::GsfElectron const &)
CaloTowerDetId	towerOf (reco::CaloCluster const &cluster, CaloTowerConstituentsMap const &towerMap)
std::vector< CaloTowerDetId >	towersOf (reco::SuperCluster const &sc, CaloTowerConstituentsMap const &towerMap, HoeMode mode=HoeMode::SingleTower)

Enumeration Type Documentation

HoeMode

enum egamma::HoeMode

strong

Enumerator
SingleTower
TowersBehindCluster

Definition at line 16 of file EgammaHadTower.h.

{ SingleTower = 0, TowersBehindCluster = 1 };

Function Documentation

classBasedElectronEnergy()

float egamma::classBasedElectronEnergy	(	reco::GsfElectron const &	electron,
		reco::BeamSpot const &	bs,
		EcalClusterFunctionBaseClass const &	crackCorrectionFunction
	)

Definition at line 296 of file ElectronEnergyCorrector.cc.

References gpuPixelDoublets::cc, alignCSCRings::corr, HPSPFTauProducerPuppi_cfi::electron, hcalRecHitTable_cff::energy, EgHLTOffHistBins_cfi::et, JetChargeProducer_cfi::exp, hcaldqm::quantity::fEnergy, l1ScoutingRun3::demux::fEt(), l1ScoutingRun3::ugmt::fEta(), EcalClusterFunctionBaseClass::getValue(), and funct::sin().

Referenced by GsfElectronAlgo::createElectron().

                                                                                                    {
  auto elClass = electron.classification();

  // new corrections from N. Chanon et al., taken from EcalClusterCorrectionObjectSpecific.cc
  float corr = 1.;
  float corr2 = 1.;
  float energy = electron.superCluster()->energy();

  //int subdet = electron.superCluster()->seed()->hitsAndFractions()[0].first.subdetId();

  if (electron.isEB()) {
    float cetacorr = fEta(electron.superCluster()->rawEnergy(), electron.superCluster()->eta(), 0) /
                     electron.superCluster()->rawEnergy();
    energy = electron.superCluster()->rawEnergy() * cetacorr;  //previously in CMSSW
    //energy = superCluster.rawEnergy()*fEta(e5x5, superCluster.seed()->eta(), 0)/e5x5;
  } else if (electron.isEE()) {
    energy = electron.superCluster()->rawEnergy() + electron.superCluster()->preshowerEnergy();
  } else {
    edm::LogWarning("ElectronEnergyCorrector::classBasedParameterizationEnergy")
        << "nor barrel neither endcap electron !";
  }

  corr = fBremEta(electron.superCluster()->phiWidth() / electron.superCluster()->etaWidth(),
                  electron.superCluster()->eta(),
                  0,
                  elClass);

  float et = energy * std::sin(2 * std::atan(std::exp(-electron.superCluster()->eta()))) / corr;

  if (electron.isEB()) {
    corr2 = corr * fEt(et, 0, elClass);
  } else if (electron.isEE()) {
    corr2 = corr * fEnergy(energy / corr, 1, elClass);
  } else {
    edm::LogWarning("ElectronEnergyCorrector::classBasedParameterizationEnergy")
        << "nor barrel neither endcap electron !";
  }

  float newEnergy = energy / corr2;

  // cracks
  double crackcor = 1.;
  for (reco::CaloCluster_iterator cIt = electron.superCluster()->clustersBegin();
       cIt != electron.superCluster()->clustersEnd();
       ++cIt) {
    const reco::CaloClusterPtr cc = *cIt;
    crackcor *= (electron.superCluster()->rawEnergy() + cc->energy() * (crackCorrectionFunction.getValue(*cc) - 1.)) /
                electron.superCluster()->rawEnergy();
  }
  newEnergy *= crackcor;

  return newEnergy;
}

classBasedElectronEnergyUncertainty()

double egamma::classBasedElectronEnergyUncertainty

(

reco::GsfElectron const &

electron

)

Definition at line 271 of file ElectronEnergyCorrector.cc.

References electrons_cff::ecalEnergy, HPSPFTauProducerPuppi_cfi::electron, and electronEnergyUncertainty().

Referenced by GsfElectronAlgo::createElectron().

                                                                                  {
  double ecalEnergy = electron.correctedEcalEnergy();
  double eleEta = electron.superCluster()->eta();
  double brem = electron.superCluster()->etaWidth() / electron.superCluster()->phiWidth();
  return egamma::electronEnergyUncertainty(electron.classification(), eleEta, brem, ecalEnergy);
}

classifyElectron()

reco::GsfElectron::Classification egamma::classifyElectron

(

reco::GsfElectron const &

)

Referenced by GsfElectronAlgo::createElectron().

correctElectronMomentum()

egamma::ElectronMomentum egamma::correctElectronMomentum	(	reco::GsfElectron const &	electron,
		TrajectoryStateOnSurface const &	vtxTsos
	)

Definition at line 23 of file ElectronMomentumCorrector.cc.

References reco::GsfElectron::BADTRACK, reco::GsfElectron::BIGBREM, HPSPFTauProducerPuppi_cfi::electron, EgHLTOffHistBins_cfi::eOverP, reco::GsfElectron::GAP, reco::GsfElectron::GOLDEN, GaussianSumUtilities1D::mode(), MultiGaussianStateTransform::multiState1D(), pfClustersFromCombinedCaloHF_cfi::scale, reco::GsfElectron::SHOWERING, mathSSE::sqrt(), reco::btau::trackMomentum, and SingleGaussianState1D::variance().

Referenced by GsfElectronAlgo::createElectron().

                                                                                                  {
  int elClass = electron.classification();

  //=======================================================================================
  // cluster energy
  //=======================================================================================

  float scEnergy = electron.correctedEcalEnergy();
  float errorEnergy = electron.correctedEcalEnergyError();

  //=======================================================================================
  // track  momentum
  //=======================================================================================

  // basic values
  float trackMomentum = electron.trackMomentumAtVtx().R();
  //float errorTrackMomentum = 999. ;

  // tracker momentum scale corrections (Mykhailo Dalchenko)
  double scale = 1.;
  if (electron.isEB()) {
    if (elClass == 0) {
      scale = 1. / (0.00104 * sqrt(trackMomentum) + 1);
    }
    if (elClass == 1) {
      scale = 1. / (0.0017 * sqrt(trackMomentum) + 0.9986);
    }
    if (elClass == 3) {
      scale = 1. / (1.004 - 0.00021 * trackMomentum);
    }
    if (elClass == 4) {
      scale = 0.995;
    }
  } else if (electron.isEE()) {
    if (elClass == 3) {
      scale = 1. / (1.01432 - 0.00201872 * trackMomentum + 0.0000142621 * trackMomentum * trackMomentum);
    }
    if (elClass == 4) {
      scale = 1. / (0.996859 - 0.000345347 * trackMomentum);
    }
  }
  if (scale < 0.)
    scale = 1.;  // CC added protection
  trackMomentum = trackMomentum * scale;

  // error (must be done after trackMomentum rescaling)
  MultiGaussianState1D qpState(MultiGaussianStateTransform::multiState1D(vtxTsos, 0));
  GaussianSumUtilities1D qpUtils(qpState);
  float errorTrackMomentum = trackMomentum * trackMomentum * sqrt(qpUtils.mode().variance());

  //=======================================================================================
  // combination
  //=======================================================================================

  float finalMomentum = electron.p4().t();  // initial
  float finalMomentumError = 999.;

  // first check for large errors

  if (errorTrackMomentum / trackMomentum > 0.5 && errorEnergy / scEnergy <= 0.5) {
    finalMomentum = scEnergy;
    finalMomentumError = errorEnergy;
  } else if (errorTrackMomentum / trackMomentum <= 0.5 && errorEnergy / scEnergy > 0.5) {
    finalMomentum = trackMomentum;
    finalMomentumError = errorTrackMomentum;
  } else if (errorTrackMomentum / trackMomentum > 0.5 && errorEnergy / scEnergy > 0.5) {
    if (errorTrackMomentum / trackMomentum < errorEnergy / scEnergy) {
      finalMomentum = trackMomentum;
      finalMomentumError = errorTrackMomentum;
    } else {
      finalMomentum = scEnergy;
      finalMomentumError = errorEnergy;
    }
  }

  // then apply the combination algorithm
  else {
    // calculate E/p and corresponding error
    float eOverP = scEnergy / trackMomentum;
    float errorEOverP = sqrt((errorEnergy / trackMomentum) * (errorEnergy / trackMomentum) +
                             (scEnergy * errorTrackMomentum / trackMomentum / trackMomentum) *
                                 (scEnergy * errorTrackMomentum / trackMomentum / trackMomentum));

    bool eleIsNotInCombination = false;
    if ((eOverP > 1 + 2.5 * errorEOverP) || (eOverP < 1 - 2.5 * errorEOverP) || (eOverP < 0.8) || (eOverP > 1.3)) {
      eleIsNotInCombination = true;
    }
    if (eleIsNotInCombination) {
      if (eOverP > 1) {
        finalMomentum = scEnergy;
        finalMomentumError = errorEnergy;
      } else {
        if (elClass == reco::GsfElectron::GOLDEN) {
          finalMomentum = scEnergy;
          finalMomentumError = errorEnergy;
        }
        if (elClass == reco::GsfElectron::BIGBREM) {
          if (scEnergy < 36) {
            finalMomentum = trackMomentum;
            finalMomentumError = errorTrackMomentum;
          } else {
            finalMomentum = scEnergy;
            finalMomentumError = errorEnergy;
          }
        }
        if (elClass == reco::GsfElectron::BADTRACK) {
          finalMomentum = scEnergy;
          finalMomentumError = errorEnergy;
        }
        if (elClass == reco::GsfElectron::SHOWERING) {
          if (scEnergy < 30) {
            finalMomentum = trackMomentum;
            finalMomentumError = errorTrackMomentum;
          } else {
            finalMomentum = scEnergy;
            finalMomentumError = errorEnergy;
          }
        }
        if (elClass == reco::GsfElectron::GAP) {
          if (scEnergy < 60) {
            finalMomentum = trackMomentum;
            finalMomentumError = errorTrackMomentum;
          } else {
            finalMomentum = scEnergy;
            finalMomentumError = errorEnergy;
          }
        }
      }
    }

    else {
      // combination
      finalMomentum = (scEnergy / errorEnergy / errorEnergy + trackMomentum / errorTrackMomentum / errorTrackMomentum) /
                      (1 / errorEnergy / errorEnergy + 1 / errorTrackMomentum / errorTrackMomentum);
      float finalMomentumVariance = 1 / (1 / errorEnergy / errorEnergy + 1 / errorTrackMomentum / errorTrackMomentum);
      finalMomentumError = sqrt(finalMomentumVariance);
    }
  }

  //=======================================================================================
  // final set
  //=======================================================================================

  auto const& oldMomentum = electron.p4();

  return {{oldMomentum.x() * finalMomentum / oldMomentum.t(),
           oldMomentum.y() * finalMomentum / oldMomentum.t(),
           oldMomentum.z() * finalMomentum / oldMomentum.t(),
           finalMomentum},
          errorTrackMomentum,
          finalMomentumError};
}

depth1HcalESum()

double egamma::depth1HcalESum	(	std::vector< CaloTowerDetId > const &	towers,
		CaloTowerCollection const &	towerCollection
	)

Definition at line 73 of file EgammaHadTower.cc.

References spr::find(), l1tHGCalTowerProducer_cfi::tower, and HLT_2024v13_cff::towers.

Referenced by EgammaHLTBcHcalIsolationProducersRegional::produce().

                                                                                                                 {
  double esum = 0.;
  for (auto const& tower : towerCollection) {
    if (std::find(towers.begin(), towers.end(), tower.id()) != towers.end()) {
      esum += tower.ietaAbs() < 18 || tower.ietaAbs() > 29 ? tower.hadEnergy() : tower.hadEnergyHeInnerLayer();
    }
  }
  return esum;
}

depth2HcalESum()

double egamma::depth2HcalESum	(	std::vector< CaloTowerDetId > const &	towers,
		CaloTowerCollection const &	towerCollection
	)

Definition at line 83 of file EgammaHadTower.cc.

References spr::find(), l1tHGCalTowerProducer_cfi::tower, and HLT_2024v13_cff::towers.

Referenced by EgammaHLTBcHcalIsolationProducersRegional::produce().

                                                                                                                 {
  double esum = 0.;
  for (auto const& tower : towerCollection) {
    if (std::find(towers.begin(), towers.end(), tower.id()) != towers.end()) {
      esum += tower.hadEnergyHeOuterLayer();
    }
  }
  return esum;
}

ecalClusterEnergyUncertaintyElectronSpecific()

float egamma::ecalClusterEnergyUncertaintyElectronSpecific

(

reco::SuperCluster const &

superCluster

)

Definition at line 4 of file ecalClusterEnergyUncertaintyElectronSpecific.cc.

References reco::CaloCluster::energy(), EgHLTOffHistBins_cfi::et, PVValHelper::eta, reco::CaloCluster::eta(), MuonErrorMatrixAnalyzer_cfi::EtaBins, reco::SuperCluster::etaWidth(), l1tPhase2CaloJetEmulator_cfi::iEta, l1tHGCalTowerMapProducer_cfi::nBinsEta, pfMETCorrectionType0_cfi::par0, pfMETCorrectionType0_cfi::par1, pfMETCorrectionType0_cfi::par2, pfMETCorrectionType0_cfi::par3, and reco::SuperCluster::phiWidth().

Referenced by GsfElectronAlgo::createElectron().

                                                                                               {
  float en = superCluster.energy();
  float eta = fabs(superCluster.eta());
  float et = en / cosh(eta);
  float brem = superCluster.etaWidth() != 0 ? superCluster.phiWidth() / superCluster.etaWidth() : 0;

  const int nBinsEta = 6;
  const float EtaBins[nBinsEta + 1] = {0.0, 0.7, 1.15, 1.44, 1.56, 2.0, 2.5};

  const int nBinsBrem = 6;
  const float BremBins[nBinsBrem + 1] = {0.8, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0};

  float par0[nBinsEta][nBinsBrem];
  float par1[nBinsEta][nBinsBrem];
  float par2[nBinsEta][nBinsBrem];
  float par3[nBinsEta][nBinsBrem];

  par0[0][0] = 0.00640519;
  par1[0][0] = 0.257578;
  par2[0][0] = 1.72437;
  par3[0][0] = 4.04686e-06;

  par0[0][1] = 0.00709569;
  par1[0][1] = 0.279844;
  par2[0][1] = 1.13789;
  par3[0][1] = 1.16239e-05;

  par0[0][2] = 0.0075544;
  par1[0][2] = 0.341346;
  par2[0][2] = 0.513396;
  par3[0][2] = 2.90054e-06;

  par0[0][3] = 0.00659365;
  par1[0][3] = 0.517649;
  par2[0][3] = -3.1847;
  par3[0][3] = 7.37152e-07;

  par0[0][4] = 0.00771696;
  par1[0][4] = 0.492897;
  par2[0][4] = -1.42222;
  par3[0][4] = 0.000358677;

  par0[0][5] = 0.00561532;
  par1[0][5] = 0.655138;
  par2[0][5] = -3.29839;
  par3[0][5] = 6.25898e-07;

  par0[1][0] = 0.00273646;
  par1[1][0] = 0.714568;
  par2[1][0] = -4.82956;
  par3[1][0] = 4.45878e-07;

  par0[1][1] = 0.00679797;
  par1[1][1] = 0.472856;
  par2[1][1] = -0.281699;
  par3[1][1] = 5.46479e-05;

  par0[1][2] = 0.00845532;
  par1[1][2] = 0.611624;
  par2[1][2] = -1.10104;
  par3[1][2] = 1.16803e-05;

  par0[1][3] = 0.00831068;
  par1[1][3] = 0.853653;
  par2[1][3] = -4.23761;
  par3[1][3] = 2.61247e-05;

  par0[1][4] = 0.00845457;
  par1[1][4] = 0.984985;
  par2[1][4] = -5.19548;
  par3[1][4] = 2.05044e-07;

  par0[1][5] = 0.0110227;
  par1[1][5] = 1.00356;
  par2[1][5] = -4.31936;
  par3[1][5] = 0.14384;

  par0[2][0] = -0.00192618;
  par1[2][0] = 1.69986;
  par2[2][0] = -16.4355;
  par3[2][0] = 1.94946e-06;

  par0[2][1] = 0.0067622;
  par1[2][1] = 0.792209;
  par2[2][1] = -1.18521;
  par3[2][1] = 0.066577;

  par0[2][2] = 0.00761595;
  par1[2][2] = 1.03058;
  par2[2][2] = -4.17237;
  par3[2][2] = 0.168543;

  par0[2][3] = 0.0119179;
  par1[2][3] = 0.910145;
  par2[2][3] = -2.14122;
  par3[2][3] = 0.00342264;

  par0[2][4] = 0.0139921;
  par1[2][4] = 1.01488;
  par2[2][4] = -2.46637;
  par3[2][4] = 0.0458434;

  par0[2][5] = 0.013724;
  par1[2][5] = 1.49078;
  par2[2][5] = -6.60661;
  par3[2][5] = 0.297821;

  par0[3][0] = -0.00197909;
  par1[3][0] = 4.40696;
  par2[3][0] = -4.88737;
  par3[3][0] = 4.99999;

  par0[3][1] = 0.0340196;
  par1[3][1] = 3.86278;
  par2[3][1] = -10.899;
  par3[3][1] = 0.130098;

  par0[3][2] = 0.0102397;
  par1[3][2] = 8.99643;
  par2[3][2] = -31.5122;
  par3[3][2] = 0.00118335;

  par0[3][3] = 0.0110891;
  par1[3][3] = 8.01794;
  par2[3][3] = -21.9038;
  par3[3][3] = 0.000245975;

  par0[3][4] = 0.0328931;
  par1[3][4] = 4.73441;
  par2[3][4] = -12.1148;
  par3[3][4] = 3.01721e-05;

  par0[3][5] = 0.0395614;
  par1[3][5] = 3.54327;
  par2[3][5] = -12.6514;
  par3[3][5] = 0.119761;

  par0[4][0] = 0.0121809;
  par1[4][0] = 0.965608;
  par2[4][0] = -4.19667;
  par3[4][0] = 0.129896;

  par0[4][1] = 0.0168951;
  par1[4][1] = 1.0218;
  par2[4][1] = -4.03078;
  par3[4][1] = 0.374291;

  par0[4][2] = 0.0213549;
  par1[4][2] = 1.29613;
  par2[4][2] = -4.89024;
  par3[4][2] = 0.0297165;

  par0[4][3] = 0.0262602;
  par1[4][3] = 1.41674;
  par2[4][3] = -5.94928;
  par3[4][3] = 0.19298;

  par0[4][4] = 0.0334892;
  par1[4][4] = 1.48572;
  par2[4][4] = -5.3175;
  par3[4][4] = 0.0157013;

  par0[4][5] = 0.0347093;
  par1[4][5] = 1.63127;
  par2[4][5] = -7.27426;
  par3[4][5] = 0.201164;

  par0[5][0] = 0.0185321;
  par1[5][0] = 0.255205;
  par2[5][0] = 1.56798;
  par3[5][0] = 5.07655e-11;

  par0[5][1] = 0.0182718;
  par1[5][1] = 0.459086;
  par2[5][1] = -0.48198;
  par3[5][1] = 0.00114946;

  par0[5][2] = 0.0175505;
  par1[5][2] = 0.92848;
  par2[5][2] = -4.52737;
  par3[5][2] = 0.154827;

  par0[5][3] = 0.0233833;
  par1[5][3] = 0.804105;
  par2[5][3] = -3.75131;
  par3[5][3] = 2.84172;

  par0[5][4] = 0.0334892;
  par1[5][4] = 1.48572;
  par2[5][4] = -5.3175;
  par3[5][4] = 0.0157013;

  par0[5][5] = 0.0347093;
  par1[5][5] = 1.63127;
  par2[5][5] = -7.27426;
  par3[5][5] = 0.201164;

  int iEtaSl = -1;
  for (int iEta = 0; iEta < nBinsEta; ++iEta) {
    if (EtaBins[iEta] <= eta && eta < EtaBins[iEta + 1]) {
      iEtaSl = iEta;
    }
  }

  int iBremSl = -1;
  for (int iBrem = 0; iBrem < nBinsBrem; ++iBrem) {
    if (BremBins[iBrem] <= brem && brem < BremBins[iBrem + 1]) {
      iBremSl = iBrem;
    }
  }

  //this code is confusing as it has no concept of under and overflow bins
  //we will use Et as an example but also applies to eta
  //underflow is 1st bin (naively currently labeled as 0 to 0.7, its really <0.7)
  //overflow is the final bin (naviely currently labeled as 5 to 10, its really >=5)
  //logic: if brem is 0<=brem <0.7 it will be already set to the 1st bin, this checks if its <0
  //logic: if brem is 5<=brem<10 it will be set to the last bin so this then checks if its >5 at which point
  //it also assigns it to the last bin. The value of 5 will have already been assigned in the for
  //loop above to the last bin so its okay that its a >5 test
  if (eta > EtaBins[nBinsEta - 1])
    iEtaSl = nBinsEta - 1;
  if (brem < BremBins[0])
    iBremSl = 0;
  if (brem > BremBins[nBinsBrem - 1])
    iBremSl = nBinsBrem - 1;

  float uncertainty = 0;
  if (et <= 5)
    uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (5 - par2[iEtaSl][iBremSl]) +
                  par3[iEtaSl][iBremSl] / ((5 - par2[iEtaSl][iBremSl]) * (5 - par2[iEtaSl][iBremSl]));
  if (et >= 200)
    uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (200 - par2[iEtaSl][iBremSl]) +
                  par3[iEtaSl][iBremSl] / ((200 - par2[iEtaSl][iBremSl]) * (200 - par2[iEtaSl][iBremSl]));

  if (et > 5 && et < 200)
    uncertainty = par0[iEtaSl][iBremSl] + par1[iEtaSl][iBremSl] / (et - par2[iEtaSl][iBremSl]) +
                  par3[iEtaSl][iBremSl] / ((et - par2[iEtaSl][iBremSl]) * (et - par2[iEtaSl][iBremSl]));

  return (uncertainty * en);
}

electronEnergyUncertainty()

double egamma::electronEnergyUncertainty	(	reco::GsfElectron::Classification	c,
		double	eta,
		double	brem,
		double	energy
	)

Definition at line 403 of file EnergyUncertaintyElectronSpecific.cc.

References HltBtagPostValidation_cff::c, hcalRecHitTable_cff::energy, PVValHelper::eta, and Exception.

Referenced by classBasedElectronEnergyUncertainty().

                                                                                                            {
  if (c == GsfElectron::GOLDEN)
    return computeEnergyUncertaintyGolden(eta, brem, energy);
  if (c == GsfElectron::BIGBREM)
    return computeEnergyUncertaintyBigbrem(eta, brem, energy);
  if (c == GsfElectron::SHOWERING)
    return computeEnergyUncertaintyShowering(eta, brem, energy);
  if (c == GsfElectron::BADTRACK)
    return computeEnergyUncertaintyBadTrack(eta, brem, energy);
  if (c == GsfElectron::GAP)
    return computeEnergyUncertaintyCracks(eta, brem, energy);
  throw cms::Exception("GsfElectronAlgo|InternalError") << "unknown classification";
}

getClosestCtfToGsf()

std::pair< TrackRef, float > egamma::getClosestCtfToGsf	(	reco::GsfTrackRef const &	gsfTrackRef,
		edm::Handle< reco::TrackCollection > const &	ctfTracksH,
		edm::soa::EtaPhiTableView	trackEtaPhiTable
	)

Definition at line 16 of file GsfElectronTools.cc.

References TrackingRecHit::all, ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), reco::deltaR2(), ALPAKA_ACCELERATOR_NAMESPACE::dR2(), reco::HitPattern::getHitPattern(), reco::HitPattern::hitPattern, SiStripPI::min, reco::HitPattern::pixelHitFilter(), edm::Handle< T >::product(), reco::HitPattern::stripTIBHitFilter(), reco::HitPattern::stripTIDHitFilter(), and reco::HitPattern::TRACK_HITS.

Referenced by LowPtGsfElectronCoreProducer::produce(), GsfElectronCoreEcalDrivenProducer::produce(), and GEDGsfElectronCoreProducer::produceElectronCore().

                                                                                  {
    float maxFracShared = 0;
    TrackRef ctfTrackRef = TrackRef();
    const TrackCollection* ctfTrackCollection = ctfTracksH.product();

    float gsfEta = gsfTrackRef->eta();
    float gsfPhi = gsfTrackRef->phi();
    const HitPattern& gsfHitPattern = gsfTrackRef->hitPattern();

    constexpr float dR2 = 0.3 * 0.3;

    unsigned int counter = 0;
    for (auto ctfTkIter = ctfTrackCollection->begin(); ctfTkIter != ctfTrackCollection->end(); ctfTkIter++, counter++) {
      // dont want to look at every single track in the event!
      using namespace edm::soa::col;
      if (reco::deltaR2(gsfEta, gsfPhi, trackTable.get<Eta>(counter), trackTable.get<Phi>(counter)) > dR2)
        continue;

      unsigned int shared = 0;
      int gsfHitCounter = 0;
      int numGsfInnerHits = 0;
      int numCtfInnerHits = 0;
      // get the CTF Track Hit Pattern
      const HitPattern& ctfHitPattern = ctfTkIter->hitPattern();

      for (auto elHitsIt = gsfTrackRef->recHitsBegin(); elHitsIt != gsfTrackRef->recHitsEnd();
           elHitsIt++, gsfHitCounter++) {
        if (!((**elHitsIt).isValid()))  //count only valid Hits
        {
          continue;
        }

        // look only in the pixels/TIB/TID
        uint32_t gsfHit = gsfHitPattern.getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter);
        if (!(HitPattern::pixelHitFilter(gsfHit) || HitPattern::stripTIBHitFilter(gsfHit) ||
              HitPattern::stripTIDHitFilter(gsfHit))) {
          continue;
        }

        numGsfInnerHits++;

        int ctfHitsCounter = 0;
        numCtfInnerHits = 0;
        for (auto ctfHitsIt = ctfTkIter->recHitsBegin(); ctfHitsIt != ctfTkIter->recHitsEnd();
             ctfHitsIt++, ctfHitsCounter++) {
          if (!((**ctfHitsIt).isValid()))  //count only valid Hits!
          {
            continue;
          }

          uint32_t ctfHit = ctfHitPattern.getHitPattern(HitPattern::TRACK_HITS, ctfHitsCounter);
          if (!(HitPattern::pixelHitFilter(ctfHit) || HitPattern::stripTIBHitFilter(ctfHit) ||
                HitPattern::stripTIDHitFilter(ctfHit))) {
            continue;
          }

          numCtfInnerHits++;

          if ((**elHitsIt).sharesInput(&(**ctfHitsIt), TrackingRecHit::all)) {
            shared++;
            break;
          }

        }  //ctfHits iterator

      }  //gsfHits iterator

      if ((numGsfInnerHits == 0) || (numCtfInnerHits == 0)) {
        continue;
      }

      if (static_cast<float>(shared) / std::min(numGsfInnerHits, numCtfInnerHits) > maxFracShared) {
        maxFracShared = static_cast<float>(shared) / std::min(numGsfInnerHits, numCtfInnerHits);
        ctfTrackRef = TrackRef(ctfTracksH, counter);
      }

    }  //ctfTrack iterator

    return {ctfTrackRef, maxFracShared};
  }

getRandomSeedFromObj()

template<typename T >

uint32_t egamma::getRandomSeedFromObj	(	const edm::Event &	iEvent,
		const T &	obj,
		size_t	nrObjs,
		size_t	objNr
	)

Definition at line 21 of file EgammaRandomSeeds.h.

References iEvent, createfilelist::int, M_PI, SiStripPI::max, getGTfromDQMFile::obj, and fileCollector::seed.

Referenced by CalibratedPhotonProducerT< T >::setSemiDetRandomSeed(), and CalibratedElectronProducerT< T >::setSemiDetRandomSeed().

                                                                                                   {
    std::seed_seq seeder = {int(iEvent.id().event()),
                            int(iEvent.id().luminosityBlock()),
                            int(iEvent.id().run()),
                            int(nrObjs),
                            int(std::numeric_limits<int>::max() * obj.phi() / M_PI) & 0xFFF,
                            int(objNr)};
    uint32_t seed = 0, tries = 10;
    do {
      seeder.generate(&seed, &seed + 1);
      tries++;
    } while (seed == 0 && tries < 10);
    return seed ? seed : iEvent.id().event() + 10000 * objNr;
  }

getRandomSeedFromSC()

uint32_t egamma::getRandomSeedFromSC	(	const edm::Event &	iEvent,
		const reco::SuperClusterRef	scRef
	)

Definition at line 5 of file EgammaRandomSeeds.cc.

References iEvent, createfilelist::int, edm::Ref< C, T, F >::key(), hltrates_dqm_sourceclient-live_cfg::offset, and fileCollector::seed.

Referenced by CalibratedPhotonProducerT< T >::setSemiDetRandomSeed(), and CalibratedElectronProducerT< T >::setSemiDetRandomSeed().

                                                                                            {
  const int offset = 0;  //for future expansion
  std::seed_seq seeder = {int(iEvent.id().event()),
                          int(iEvent.id().luminosityBlock()),
                          int(iEvent.id().run()),
                          int(scRef->seed()->seed().rawId()),
                          int(scRef->seed()->hitsAndFractions().size()),
                          offset};
  uint32_t seed = 0, tries = 10;
  do {
    seeder.generate(&seed, &seed + 1);
    tries++;
  } while (seed == 0 && tries < 10);
  return seed ? seed : iEvent.id().event() + 10000 * scRef.key();
}

hasActiveHcal()

bool egamma::hasActiveHcal	(	std::vector< CaloTowerDetId > const &	towers,
		CaloTowerConstituentsMap const &	towerMap,
		HcalChannelQuality const &	hcalQuality,
		HcalTopology const &	hcalTopology
	)

Definition at line 93 of file EgammaHadTower.cc.

References CaloTowerConstituentsMap::constituentsOf(), gather_cfg::cout, TauDecayModes::dec, HcalDetId::depth(), HcalChannelStatus::getValue(), HcalCondObjectContainer< Item >::getValues(), DetId::Hcal, HcalBarrel, HcalChannelStatus::HcalCellDead, HcalChannelStatus::HcalCellMask, HcalChannelStatus::HcalCellOff, HcalEndcap, HcalTopology::idFront(), HcalDetId::ieta(), HcalDetId::iphi(), DetId::rawId(), mps_update::status, HcalDetId::subdet(), and HLT_2024v13_cff::towers.

Referenced by ElectronHcalHelper::hasActiveHcal().

                                                             {
  bool active = false;
  int statusMask = ((1 << HcalChannelStatus::HcalCellOff) | (1 << HcalChannelStatus::HcalCellMask) |
                    (1 << HcalChannelStatus::HcalCellDead));
#ifdef EDM_ML_DEBUG
  std::cout << "DEBUG: hasActiveHcal called with " << towers.size() << " detids. First tower detid ieta "
            << towers.front().ieta() << " iphi " << towers.front().iphi() << std::endl;
#endif
  for (auto towerid : towers) {
    unsigned int ngood = 0, nbad = 0;
    for (DetId id : towerMap.constituentsOf(towerid)) {
      if (id.det() != DetId::Hcal) {
        continue;
      }
      HcalDetId hid(id);
      if (hid.subdet() != HcalBarrel && hid.subdet() != HcalEndcap)
        continue;
#ifdef EDM_ML_DEBUG
      std::cout << "EgammaHadTower DetId " << std::hex << id.rawId() << "  hid.rawId  " << hid.rawId() << std::dec
                << "   sub " << hid.subdet() << "   ieta " << hid.ieta() << "   iphi " << hid.iphi() << "   depth "
                << hid.depth() << std::endl;
#endif
      // Sunanda's fix for 2017 Plan1
      // and removed protection
      int status =
          hcalQuality.getValues((DetId)(hcalTopology.idFront(HcalDetId(id))), /*throwOnFail=*/true)->getValue();

#ifdef EDM_ML_DEBUG
      std::cout << "channels status = " << std::hex << status << std::dec << "  int value = " << status << std::endl;
#endif

      if (status & statusMask) {
#ifdef EDM_ML_DEBUG
        std::cout << "          BAD!" << std::endl;
#endif
        nbad++;
      } else {
        ngood++;
      }
    }
#ifdef EDM_ML_DEBUG
    std::cout << "    overall ngood " << ngood << " nbad " << nbad << "\n";
#endif
    if (nbad == 0 || (ngood > 0 && nbad < ngood)) {
      active = true;
    }
  }
  return active;
}

isBetterElectron()

bool egamma::isBetterElectron	(	reco::GsfElectron const &	e1,
		reco::GsfElectron const &	e2
	)

Definition at line 16 of file EgAmbiguityTools.cc.

References funct::abs(), StorageManager_cfg::e1, and reco::GsfElectron::eSuperClusterOverP().

Referenced by isInnermostElectron(), and GsfElectronProducer::setAmbiguityData().

                                                                              {
    return (std::abs(e1.eSuperClusterOverP() - 1) < std::abs(e2.eSuperClusterOverP() - 1));
  }

isInnermostElectron()

bool egamma::isInnermostElectron	(	reco::GsfElectron const &	e1,
		reco::GsfElectron const &	e2
	)

Definition at line 20 of file EgAmbiguityTools.cc.

References StorageManager_cfg::e1, pixelTopology::getLayer(), reco::GsfElectron::gsfTrack(), and isBetterElectron().

Referenced by GsfElectronProducer::setAmbiguityData().

                                                                                 {
    // retreive first valid hit
    int gsfHitCounter1 = 0;
    for (auto const& elHit : e1.gsfTrack()->recHits()) {
      if (elHit->isValid())
        break;
      gsfHitCounter1++;
    }

    int gsfHitCounter2 = 0;
    for (auto const& elHit : e2.gsfTrack()->recHits()) {
      if (elHit->isValid())
        break;
      gsfHitCounter2++;
    }

    uint32_t gsfHit1 = e1.gsfTrack()->hitPattern().getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter1);
    uint32_t gsfHit2 = e2.gsfTrack()->hitPattern().getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter2);

    if (HitPattern::getSubStructure(gsfHit1) != HitPattern::getSubStructure(gsfHit2)) {
      return (HitPattern::getSubStructure(gsfHit1) < HitPattern::getSubStructure(gsfHit2));
    } else if (HitPattern::getLayer(gsfHit1) != HitPattern::getLayer(gsfHit2)) {
      return (HitPattern::getLayer(gsfHit1) < HitPattern::getLayer(gsfHit2));
    } else {
      return isBetterElectron(e1, e2);
    }
  }

sharedDets()

int egamma::sharedDets	(	reco::GsfTrackRef const &	gsfTrackRef1,
		reco::GsfTrackRef const &	gsfTrackRef2
	)

Definition at line 95 of file EgAmbiguityTools.cc.

References reco::HitPattern::getHitPattern().

                                                                                   {
    //get the Hit Pattern for the gsfTracks
    const HitPattern& gsfHitPattern1 = gsfTrackRef1->hitPattern();
    const HitPattern& gsfHitPattern2 = gsfTrackRef2->hitPattern();

    unsigned int shared = 0;

    int gsfHitCounter1 = 0;
    for (trackingRecHit_iterator elHitsIt1 = gsfTrackRef1->recHitsBegin(); elHitsIt1 != gsfTrackRef1->recHitsEnd();
         elHitsIt1++, gsfHitCounter1++) {
      if (!((**elHitsIt1).isValid())) {
        //count only valid Hits
        continue;
      }
      //if (gsfHitCounter1>1) continue; // to test only the first hit of the track 1
      uint32_t gsfHit = gsfHitPattern1.getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter1);
      if (!(HitPattern::pixelHitFilter(gsfHit) || HitPattern::stripTIBHitFilter(gsfHit) ||
            HitPattern::stripTOBHitFilter(gsfHit) || HitPattern::stripTECHitFilter(gsfHit) ||
            HitPattern::stripTIDHitFilter(gsfHit))) {
        continue;
      }

      int gsfHitsCounter2 = 0;
      for (trackingRecHit_iterator gsfHitsIt2 = gsfTrackRef2->recHitsBegin(); gsfHitsIt2 != gsfTrackRef2->recHitsEnd();
           gsfHitsIt2++, gsfHitsCounter2++) {
        if (!((**gsfHitsIt2).isValid())) {
          //count only valid Hits!
          continue;
        }

        uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(HitPattern::TRACK_HITS, gsfHitsCounter2);
        if (!(HitPattern::pixelHitFilter(gsfHit2) || HitPattern::stripTIBHitFilter(gsfHit2) ||
              HitPattern::stripTOBHitFilter(gsfHit2) || HitPattern::stripTECHitFilter(gsfHit2) ||
              HitPattern::stripTIDHitFilter(gsfHit2)))
          continue;
        if ((**elHitsIt1).geographicalId() == (**gsfHitsIt2).geographicalId())
          shared++;
      }  //gsfHits2 iterator
    }    //gsfHits1 iterator

    //std::cout << "[sharedHits] number of shared dets " << shared << std::endl;
    //return shared/min(gsfTrackRef1->numberOfValidHits(),gsfTrackRef2->numberOfValidHits());
    return shared;
  }

sharedEnergy() [1/2]

float egamma::sharedEnergy	(	reco::CaloCluster const &	clu1,
		reco::CaloCluster const &	clu2,
		EcalRecHitCollection const &	barrelRecHits,
		EcalRecHitCollection const &	endcapRecHits
	)

Definition at line 140 of file EgAmbiguityTools.cc.

References HLT_2024v13_cff::barrelRecHits, EcalBarrel, EcalEndcap, HLT_2024v13_cff::endcapRecHits, HLT_2024v13_cff::fractionShared, and reco::CaloCluster::hitsAndFractions().

Referenced by RealisticHitToClusterAssociator::findAndMergeInvisibleClusters(), TSToSimTSHitLCAssociatorByEnergyScoreImpl::makeConnections(), GsfElectronProducer::setAmbiguityData(), and sharedEnergy().

                                                                {
    double fractionShared = 0;

    for (auto const& h1 : clu1.hitsAndFractions()) {
      for (auto const& h2 : clu2.hitsAndFractions()) {
        if (h1.first != h2.first)
          continue;

        // here we have common Xtal id
        EcalRecHitCollection::const_iterator itt;
        if (h1.first.subdetId() == EcalBarrel) {
          if ((itt = barrelRecHits.find(h1.first)) != barrelRecHits.end())
            fractionShared += itt->energy();
        } else if (h1.first.subdetId() == EcalEndcap) {
          if ((itt = endcapRecHits.find(h1.first)) != endcapRecHits.end())
            fractionShared += itt->energy();
        }
      }
    }

    //std::cout << "[sharedEnergy] shared energy /min(energy1,energy2) " << fractionShared << std::endl;
    return fractionShared;
  }

sharedEnergy() [2/2]

float egamma::sharedEnergy	(	reco::SuperClusterRef const &	sc1,
		reco::SuperClusterRef const &	sc2,
		EcalRecHitCollection const &	barrelRecHits,
		EcalRecHitCollection const &	endcapRecHits
	)

Definition at line 167 of file EgAmbiguityTools.cc.

References HLT_2024v13_cff::barrelRecHits, HLT_2024v13_cff::endcapRecHits, and sharedEnergy().

                                                                {
    double energyShared = 0;
    for (CaloCluster_iterator icl1 = sc1->clustersBegin(); icl1 != sc1->clustersEnd(); icl1++) {
      for (CaloCluster_iterator icl2 = sc2->clustersBegin(); icl2 != sc2->clustersEnd(); icl2++) {
        energyShared += sharedEnergy(**icl1, **icl2, barrelRecHits, endcapRecHits);
      }
    }
    return energyShared;
  }

sharedHits()

int egamma::sharedHits	(	reco::GsfTrackRef const &	gsfTrackRef1,
		reco::GsfTrackRef const &	gsfTrackRef2
	)

Definition at line 48 of file EgAmbiguityTools.cc.

References reco::HitPattern::getHitPattern(), and TrackingRecHit::some.

                                                                                   {
    //get the Hit Pattern for the gsfTracks
    HitPattern const& gsfHitPattern1 = gsfTrackRef1->hitPattern();
    HitPattern const& gsfHitPattern2 = gsfTrackRef2->hitPattern();

    unsigned int shared = 0;

    int gsfHitCounter1 = 0;
    for (trackingRecHit_iterator elHitsIt1 = gsfTrackRef1->recHitsBegin(); elHitsIt1 != gsfTrackRef1->recHitsEnd();
         elHitsIt1++, gsfHitCounter1++) {
      if (!(*elHitsIt1)->isValid()) {
        //count only valid Hits
        continue;
      }
      //if (gsfHitCounter1>1) continue; // test only the first hit of the track 1
      uint32_t gsfHit = gsfHitPattern1.getHitPattern(HitPattern::TRACK_HITS, gsfHitCounter1);
      if (!(HitPattern::pixelHitFilter(gsfHit) || HitPattern::stripTIBHitFilter(gsfHit) ||
            HitPattern::stripTOBHitFilter(gsfHit) || HitPattern::stripTECHitFilter(gsfHit) ||
            HitPattern::stripTIDHitFilter(gsfHit))) {
        continue;
      }

      int gsfHitsCounter2 = 0;
      for (trackingRecHit_iterator gsfHitsIt2 = gsfTrackRef2->recHitsBegin(); gsfHitsIt2 != gsfTrackRef2->recHitsEnd();
           gsfHitsIt2++, gsfHitsCounter2++) {
        if (!(**gsfHitsIt2).isValid()) {
          //count only valid Hits
          continue;
        }
        uint32_t gsfHit2 = gsfHitPattern2.getHitPattern(HitPattern::TRACK_HITS, gsfHitsCounter2);
        if (!(HitPattern::pixelHitFilter(gsfHit2) || HitPattern::stripTIBHitFilter(gsfHit2) ||
              HitPattern::stripTOBHitFilter(gsfHit2) || HitPattern::stripTECHitFilter(gsfHit2) ||
              HitPattern::stripTIDHitFilter(gsfHit2))) {
          continue;
        }
        if ((*elHitsIt1)->sharesInput(&(**gsfHitsIt2), TrackingRecHit::some)) {
          //if (comp.equals(&(**elHitsIt1),&(**gsfHitsIt2))) {
          shared++;
        }
      }  //gsfHits2 iterator
    }    //gsfHits1 iterator

    //std::cout << "[sharedHits] number of shared hits " << shared << std::endl;
    return shared;
  }

simpleElectronEnergyUncertainty()

double egamma::simpleElectronEnergyUncertainty

(

reco::GsfElectron const &

electron

)

Definition at line 278 of file ElectronEnergyCorrector.cc.

References electrons_cff::ecalEnergy, HPSPFTauProducerPuppi_cfi::electron, and relativeConstraints::error.

Referenced by GsfElectronAlgo::createElectron().

                                                                              {
  double error = 999.;
  double ecalEnergy = electron.correctedEcalEnergy();

  if (electron.isEB()) {
    float parEB[3] = {5.24e-02, 2.01e-01, 1.00e-02};
    error = ecalEnergy * energyError(ecalEnergy, parEB);
  } else if (electron.isEE()) {
    float parEE[3] = {1.46e-01, 9.21e-01, 1.94e-03};
    error = ecalEnergy * energyError(ecalEnergy, parEE);
  } else {
    edm::LogWarning("ElectronEnergyCorrector::simpleParameterizationUncertainty")
        << "nor barrel neither endcap electron !";
  }

  return error;
}

towerOf()

CaloTowerDetId egamma::towerOf	(	reco::CaloCluster const &	cluster,
		CaloTowerConstituentsMap const &	towerMap
	)

Definition at line 10 of file EgammaHadTower.cc.

References DetId::det(), DetId::Ecal, reco::CaloCluster::hitsAndFractions(), l1ctLayer2EG_cff::id, reco::CaloCluster::seed(), l1tHGCalTowerProducer_cfi::tower, and CaloTowerConstituentsMap::towerOf().

Referenced by EgammaHcalIsolation::getHcalESumBc(), EgammaHcalIsolation::getHcalEtSumBc(), ecaldqm::TrigPrimTask::runOnDigis(), and towersOf().

                                                                                                       {
  DetId detid = cluster.seed();
  if (detid.det() != DetId::Ecal) {
    // Basic clusters of hybrid super-cluster do not have the seed set; take the first DetId instead
    // Should be checked . The single Tower Mode should be favoured until fixed
    detid = cluster.hitsAndFractions()[0].first;
    if (detid.det() != DetId::Ecal) {
      CaloTowerDetId tower;
      return tower;
    }
  }
  CaloTowerDetId id(towerMap.towerOf(detid));
  return id;
}

towersOf()

std::vector< CaloTowerDetId > egamma::towersOf	(	reco::SuperCluster const &	sc,
		CaloTowerConstituentsMap const &	towerMap,
		HoeMode	mode = HoeMode::SingleTower
	)

Definition at line 25 of file EgammaHadTower.cc.

References alignmentValidation::c1, reco::SuperCluster::clustersBegin(), reco::SuperCluster::clustersEnd(), ALPAKA_ACCELERATOR_NAMESPACE::brokenline::constexpr(), gather_cfg::cout, spr::find(), l1ctLayer2EG_cff::id, ALPAKA_ACCELERATOR_NAMESPACE::vertexFinder::it, ALCARECOPromptCalibProdSiPixelAli0T_cff::mode, reco::SuperCluster::seed(), jetUpdater_cfi::sort, towerOf(), and HLT_2024v13_cff::towers.

Referenced by ElectronHcalHelper::hcalTowersBehindClusters(), and EgammaHLTBcHcalIsolationProducersRegional::produce().

                                                           {
  constexpr unsigned int nMaxClusters = 4;

  std::vector<CaloTowerDetId> towers;
  std::vector<reco::CaloClusterPtr> orderedClusters;

  // in this mode, check only the tower behind the seed
  if (mode == HoeMode::SingleTower) {
    towers.push_back(towerOf(*sc.seed(), towerMap));
  }

  // in this mode check the towers behind each basic cluster
  if (mode == HoeMode::TowersBehindCluster) {
    // Loop on the basic clusters
    for (auto it = sc.clustersBegin(); it != sc.clustersEnd(); ++it) {
      orderedClusters.push_back(*it);
    }
    std::sort(orderedClusters.begin(), orderedClusters.end(), [](auto& c1, auto& c2) { return (*c1 > *c2); });
    unsigned nclusters = orderedClusters.size();
    for (unsigned iclus = 0; iclus < nclusters && iclus < nMaxClusters; ++iclus) {
      // Get the tower
      CaloTowerDetId id = towerOf(*(orderedClusters[iclus]), towerMap);
#ifdef EDM_ML_DEBUG
      std::cout << "CaloTowerId " << id << std::endl;
#endif
      if (std::find(towers.begin(), towers.end(), id) == towers.end()) {
        towers.push_back(id);
      }
    }
  }
  //  if(towers.size() > 4) {
  //    std::cout << " NTOWERS " << towers.size() << " ";
  //    for(unsigned i=0; i<towers.size() ; ++i) {
  //      std::cout << towers[i] << " ";
  //    }
  //    std::cout <<  std::endl;
  //    for ( unsigned iclus=0 ; iclus < orderedClusters.size(); ++iclus) {
  //      // Get the tower
  //      CaloTowerDetId id = towerOf(*(orderedClusters[iclus]));
  //      std::cout << " Pos " << orderedClusters[iclus]->position() << " " << orderedClusters[iclus]->energy() << " " << id ;
  //    }
  //    std::cout << std::endl;
  //  }
  return towers;
}