CMS 3D CMS Logo

List of all members | Public Member Functions | Private Member Functions | Private Attributes | Friends
PFAlgo Class Reference

#include <PFAlgo.h>

Public Member Functions

void checkCleaning (const reco::PFRecHitCollection &cleanedHF)
 Check HF Cleaning. More...
 
reco::PFCandidateCollectiongetCleanedCandidates ()
 
PFMuonAlgogetPFMuonAlgo ()
 
reco::PFCandidateCollection makeConnectedCandidates ()
 
 PFAlgo (double nSigmaECAL, double nSigmaHCAL, double nSigmaHFEM, double nSigmaHFHAD, std::vector< double > resolHF_square, PFEnergyCalibration &calibration, PFEnergyCalibrationHF &thepfEnergyCalibrationHF, const edm::ParameterSet &pset)
 constructor More...
 
void reconstructParticles (const reco::PFBlockHandle &blockHandle, PFEGammaFilters const *pfegamma)
 reconstruct particles More...
 
void setCandConnectorParameters (const edm::ParameterSet &iCfgCandConnector)
 
void setCandConnectorParameters (bool bCorrect, bool bCalibPrimary, double dptRel_PrimaryTrack, double dptRel_MergedTrack, double ptErrorSecondary, const std::vector< double > &nuclCalibFactors)
 
void setDisplacedVerticesParameters (bool rejectTracks_Bad, bool rejectTracks_Step45, bool usePFNuclearInteractions, bool usePFConversions, bool usePFDecays, double dptRel_DispVtx)
 
void setEGammaCollections (const edm::View< reco::PFCandidate > &pfEgammaCandidates, const edm::ValueMap< reco::GsfElectronRef > &valueMapGedElectrons, const edm::ValueMap< reco::PhotonRef > &valueMapGedPhotons)
 
void setEGammaParameters (bool use_EGammaFilters, bool useProtectionsForJetMET)
 
void setEGElectronCollection (const reco::GsfElectronCollection &egelectrons)
 
void setHOTag (bool ho)
 
void setMuonHandle (const edm::Handle< reco::MuonCollection > &muons)
 
void setPFVertexParameters (bool useVertex, reco::VertexCollection const &primaryVertices)
 
void setPostHFCleaningParameters (bool postHFCleaning, const edm::ParameterSet &pfHFCleaningParams)
 

Private Member Functions

void associatePSClusters (unsigned iEcal, reco::PFBlockElement::Type psElementType, const reco::PFBlock &block, const edm::OwnVector< reco::PFBlockElement > &elements, const reco::PFBlock::LinkData &linkData, std::vector< bool > &active, std::vector< double > &psEne)
 Associate PS clusters to a given ECAL cluster, and return their energy. More...
 
bool checkAndReconstructSecondaryInteraction (const reco::PFBlockRef &blockref, const edm::OwnVector< reco::PFBlockElement > &elements, bool isActive, int iElement)
 
bool checkGoodTrackDeadHcal (const reco::TrackRef &trackRef, bool hasDeadHcal)
 
bool checkHasDeadHcal (const std::multimap< double, unsigned > &hcalElems, const std::vector< bool > &deadArea)
 
void conversionAlgo (const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active)
 
void createCandidatesECAL (const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active, const reco::PFBlockRef &blockref, ElementIndices &inds, std::vector< bool > &deadArea)
 
void createCandidatesHCAL (const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active, const reco::PFBlockRef &blockref, ElementIndices &inds, std::vector< bool > &deadArea)
 
void createCandidatesHCALUnlinked (const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active, const reco::PFBlockRef &blockref, ElementIndices &inds, std::vector< bool > &deadArea)
 
void createCandidatesHF (const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active, const reco::PFBlockRef &blockref, ElementIndices &inds)
 
int decideType (const edm::OwnVector< reco::PFBlockElement > &elements, const reco::PFBlockElement::Type type, std::vector< bool > &active, ElementIndices &inds, std::vector< bool > &deadArea, unsigned int iEle)
 
void egammaFilters (const reco::PFBlockRef &blockref, std::vector< bool > &active, PFEGammaFilters const *pfegamma)
 
void elementLoop (const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active, const reco::PFBlockRef &blockref, ElementIndices &inds, std::vector< bool > &deadArea)
 
double hfEnergyResolution (double clusterEnergy) const
 
bool isFromSecInt (const reco::PFBlockElement &eTrack, std::string order) const
 
double neutralHadronEnergyResolution (double clusterEnergy, double clusterEta) const
 todo: use PFClusterTools for this More...
 
double nSigmaHCAL (double clusterEnergy, double clusterEta) const
 
double nSigmaHFEM (double clusterEnergy) const
 
double nSigmaHFHAD (double clusterEnergy) const
 
void postCleaning ()
 
void processBlock (const reco::PFBlockRef &blockref, std::list< reco::PFBlockRef > &hcalBlockRefs, std::list< reco::PFBlockRef > &ecalBlockRefs, PFEGammaFilters const *pfegamma)
 
unsigned reconstructCluster (const reco::PFCluster &cluster, double particleEnergy, bool useDirection=false, double particleX=0., double particleY=0., double particleZ=0.)
 
unsigned reconstructTrack (const reco::PFBlockElement &elt, bool allowLoose=false)
 
bool recoTracksNotHCAL (const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, const reco::PFBlockRef &blockref, std::vector< bool > &active, bool goodTrackDeadHcal, bool hasDeadHcal, unsigned int iTrack, std::multimap< double, unsigned > &ecalElems, reco::TrackRef &trackRef)
 
void relinkTrackToHcal (const reco::PFBlock &block, std::multimap< double, unsigned > &ecalElems, std::multimap< double, unsigned > &hcalElems, const std::vector< bool > &active, reco::PFBlock::LinkData &linkData, unsigned int iTrack)
 
void setHcalDepthInfo (reco::PFCandidate &cand, const reco::PFCluster &cluster) const
 

Private Attributes

PFEnergyCalibrationcalibration_
 
PFCandConnector connector_
 
double dptRel_DispVtx_
 
std::vector< double > factors45_
 
float goodPixelTrackDeadHcal_chi2n_
 
float goodPixelTrackDeadHcal_dxy_
 
float goodPixelTrackDeadHcal_dz_
 
int goodPixelTrackDeadHcal_maxLost3Hit_
 
int goodPixelTrackDeadHcal_maxLost4Hit_
 
float goodPixelTrackDeadHcal_maxPt_
 
float goodPixelTrackDeadHcal_minEta_
 
float goodPixelTrackDeadHcal_ptErrRel_
 
float goodTrackDeadHcal_chi2n_
 
float goodTrackDeadHcal_dxy_
 
int goodTrackDeadHcal_layers_
 
float goodTrackDeadHcal_ptErrRel_
 Variables for track cleaning in bad HCal areas. More...
 
float goodTrackDeadHcal_validFr_
 
double maxDeltaPhiPt_
 
double maxSignificance_
 
double minDeltaMet_
 
double minHFCleaningPt_
 
double minSignificance_
 
double minSignificanceReduction_
 
std::vector< double > muonECAL_
 
edm::Handle< reco::MuonCollectionmuonHandle_
 
std::vector< double > muonHCAL_
 Variables for muons and fakes. More...
 
std::vector< double > muonHO_
 
const double nSigmaECAL_
 number of sigma to judge energy excess in ECAL More...
 
const double nSigmaEConstHCAL = 100.
 
const double nSigmaEConstHFEM = 100.
 
const double nSigmaEConstHFHAD = 100.
 
const double nSigmaHCAL_
 number of sigma to judge energy excess in HCAL More...
 
const double nSigmaHFEM_
 number of sigma to judge energy excess in HF More...
 
const double nSigmaHFHAD_
 
double nSigmaTRACK_
 
int nVtx_
 
std::unique_ptr< reco::PFCandidateCollectionpfCandidates_
 
reco::PFCandidateCollection pfCleanedCandidates_
 
const edm::View< reco::PFCandidate > * pfEgammaCandidates_
 
std::unique_ptr< PFMuonAlgopfmu_
 
bool postHFCleaning_
 
bool postMuonCleaning_
 
reco::Vertex primaryVertex_
 
double ptError_
 
bool rejectTracks_Bad_
 
bool rejectTracks_Step45_
 
const std::vector< double > resolHF_square_
 
PFEnergyCalibrationHFthepfEnergyCalibrationHF_
 
double useBestMuonTrack_
 
bool useEGammaFilters_
 Variables for NEW EGAMMA selection. More...
 
bool useHO_
 
bool usePFConversions_
 
bool usePFDecays_
 
bool usePFMuonMomAssign_
 
bool usePFNuclearInteractions_
 
bool useProtectionsForJetMET_
 
bool useVertices_ = false
 
const edm::ValueMap< reco::GsfElectronRef > * valueMapGedElectrons_
 
const edm::ValueMap< reco::PhotonRef > * valueMapGedPhotons_
 

Friends

std::ostream & operator<< (std::ostream &out, const PFAlgo &algo)
 

Detailed Description

Definition at line 53 of file PFAlgo.h.

Constructor & Destructor Documentation

◆ PFAlgo()

PFAlgo::PFAlgo ( double  nSigmaECAL,
double  nSigmaHCAL,
double  nSigmaHFEM,
double  nSigmaHFHAD,
std::vector< double >  resolHF_square,
PFEnergyCalibration calibration,
PFEnergyCalibrationHF thepfEnergyCalibrationHF,
const edm::ParameterSet pset 
)

constructor

Definition at line 15 of file PFAlgo.cc.

References cms::cuda::assert(), factors45_, goodPixelTrackDeadHcal_chi2n_, goodPixelTrackDeadHcal_dxy_, goodPixelTrackDeadHcal_dz_, goodPixelTrackDeadHcal_maxLost3Hit_, goodPixelTrackDeadHcal_maxLost4Hit_, goodPixelTrackDeadHcal_maxPt_, goodPixelTrackDeadHcal_minEta_, goodPixelTrackDeadHcal_ptErrRel_, goodTrackDeadHcal_chi2n_, goodTrackDeadHcal_dxy_, goodTrackDeadHcal_layers_, goodTrackDeadHcal_ptErrRel_, goodTrackDeadHcal_validFr_, muonECAL_, muonHCAL_, muonHO_, nSigmaTRACK_, pfmu_, HLT_2022v12_cff::postMuonCleaning, muonDTDigis_cfi::pset, ptError_, and resolHF_square_.

24  nSigmaECAL_(nSigmaECAL),
30  thepfEnergyCalibrationHF_(thepfEnergyCalibrationHF),
31  connector_() {
32  const edm::ParameterSet pfMuonAlgoParams = pset.getParameter<edm::ParameterSet>("PFMuonAlgoParameters");
33  bool postMuonCleaning = pset.getParameter<bool>("postMuonCleaning");
34  pfmu_ = std::make_unique<PFMuonAlgo>(pfMuonAlgoParams, postMuonCleaning);
35 
36  // HF resolution parameters
37  assert(resolHF_square_.size() == 3); // make sure that stochastic, constant, noise (i.e. three) terms are specified.
38 
39  // Muon parameters
40  muonHCAL_ = pset.getParameter<std::vector<double>>("muon_HCAL");
41  muonECAL_ = pset.getParameter<std::vector<double>>("muon_ECAL");
42  muonHO_ = pset.getParameter<std::vector<double>>("muon_HO");
43  assert(muonHCAL_.size() == 2 && muonECAL_.size() == 2 && muonHO_.size() == 2);
44  nSigmaTRACK_ = pset.getParameter<double>("nsigma_TRACK");
45  ptError_ = pset.getParameter<double>("pt_Error");
46  factors45_ = pset.getParameter<std::vector<double>>("factors_45");
47  assert(factors45_.size() == 2);
48 
49  // Bad Hcal Track Parameters
50  goodTrackDeadHcal_ptErrRel_ = pset.getParameter<double>("goodTrackDeadHcal_ptErrRel");
51  goodTrackDeadHcal_chi2n_ = pset.getParameter<double>("goodTrackDeadHcal_chi2n");
52  goodTrackDeadHcal_layers_ = pset.getParameter<uint32_t>("goodTrackDeadHcal_layers");
53  goodTrackDeadHcal_validFr_ = pset.getParameter<double>("goodTrackDeadHcal_validFr");
54  goodTrackDeadHcal_dxy_ = pset.getParameter<double>("goodTrackDeadHcal_dxy");
55 
56  goodPixelTrackDeadHcal_minEta_ = pset.getParameter<double>("goodPixelTrackDeadHcal_minEta");
57  goodPixelTrackDeadHcal_maxPt_ = pset.getParameter<double>("goodPixelTrackDeadHcal_maxPt");
58  goodPixelTrackDeadHcal_ptErrRel_ = pset.getParameter<double>("goodPixelTrackDeadHcal_ptErrRel");
59  goodPixelTrackDeadHcal_chi2n_ = pset.getParameter<double>("goodPixelTrackDeadHcal_chi2n");
60  goodPixelTrackDeadHcal_maxLost3Hit_ = pset.getParameter<int32_t>("goodPixelTrackDeadHcal_maxLost3Hit");
61  goodPixelTrackDeadHcal_maxLost4Hit_ = pset.getParameter<int32_t>("goodPixelTrackDeadHcal_maxLost4Hit");
62  goodPixelTrackDeadHcal_dxy_ = pset.getParameter<double>("goodPixelTrackDeadHcal_dxy");
63  goodPixelTrackDeadHcal_dz_ = pset.getParameter<double>("goodPixelTrackDeadHcal_dz");
64 }
double ptError_
Definition: PFAlgo.h:297
float goodTrackDeadHcal_dxy_
Definition: PFAlgo.h:305
std::vector< double > muonHCAL_
Variables for muons and fakes.
Definition: PFAlgo.h:293
const double nSigmaHFEM_
number of sigma to judge energy excess in HF
Definition: PFAlgo.h:251
const double nSigmaHCAL_
number of sigma to judge energy excess in HCAL
Definition: PFAlgo.h:248
std::unique_ptr< reco::PFCandidateCollection > pfCandidates_
Definition: PFAlgo.h:226
const std::vector< double > resolHF_square_
Definition: PFAlgo.h:255
float goodPixelTrackDeadHcal_minEta_
Definition: PFAlgo.h:307
PFCandConnector connector_
Definition: PFAlgo.h:290
assert(be >=bs)
float goodPixelTrackDeadHcal_dxy_
Definition: PFAlgo.h:313
std::vector< double > factors45_
Definition: PFAlgo.h:298
double nSigmaHFEM(double clusterEnergy) const
Definition: PFAlgo.cc:3390
float goodTrackDeadHcal_validFr_
Definition: PFAlgo.h:304
std::vector< double > muonECAL_
Definition: PFAlgo.h:294
PFEnergyCalibration & calibration_
Definition: PFAlgo.h:257
int goodPixelTrackDeadHcal_maxLost3Hit_
Definition: PFAlgo.h:311
float goodPixelTrackDeadHcal_maxPt_
Definition: PFAlgo.h:308
std::vector< reco::PFCandidate > PFCandidateCollection
collection of PFCandidates
double nSigmaHFHAD(double clusterEnergy) const
Definition: PFAlgo.cc:3395
int goodTrackDeadHcal_layers_
Definition: PFAlgo.h:303
std::vector< double > muonHO_
Definition: PFAlgo.h:295
float goodPixelTrackDeadHcal_ptErrRel_
Definition: PFAlgo.h:309
const double nSigmaECAL_
number of sigma to judge energy excess in ECAL
Definition: PFAlgo.h:245
PFEnergyCalibrationHF & thepfEnergyCalibrationHF_
Definition: PFAlgo.h:258
double nSigmaHCAL(double clusterEnergy, double clusterEta) const
Definition: PFAlgo.cc:3371
float goodTrackDeadHcal_ptErrRel_
Variables for track cleaning in bad HCal areas.
Definition: PFAlgo.h:301
int goodPixelTrackDeadHcal_maxLost4Hit_
Definition: PFAlgo.h:312
const double nSigmaHFHAD_
Definition: PFAlgo.h:252
float goodPixelTrackDeadHcal_chi2n_
Definition: PFAlgo.h:310
double nSigmaTRACK_
Definition: PFAlgo.h:296
std::unique_ptr< PFMuonAlgo > pfmu_
Definition: PFAlgo.h:262
float goodTrackDeadHcal_chi2n_
Definition: PFAlgo.h:302
float goodPixelTrackDeadHcal_dz_
Definition: PFAlgo.h:314

Member Function Documentation

◆ associatePSClusters()

void PFAlgo::associatePSClusters ( unsigned  iEcal,
reco::PFBlockElement::Type  psElementType,
const reco::PFBlock block,
const edm::OwnVector< reco::PFBlockElement > &  elements,
const reco::PFBlock::LinkData linkData,
std::vector< bool > &  active,
std::vector< double > &  psEne 
)
private

Associate PS clusters to a given ECAL cluster, and return their energy.

Definition at line 3425 of file PFAlgo.cc.

References cms::cuda::assert(), groupFilesInBlocks::block, reco::PFBlockElement::ECAL, bookConverter::elements, edm::Ref< C, T, F >::isNull(), and reco::PFBlock::LINKTEST_ALL.

Referenced by recoTracksNotHCAL().

3431  {
3432  // Find all PS clusters with type psElement associated to ECAL cluster iEcal,
3433  // within all PFBlockElement "elements" of a given PFBlock "block"
3434  // psElement can be reco::PFBlockElement::PS1 or reco::PFBlockElement::PS2
3435  // Returns a vector of PS cluster energies, and updates the "active" vector.
3436 
3437  // Find all PS clusters linked to the iEcal cluster
3438  std::multimap<double, unsigned> sortedPS;
3439  block.associatedElements(iEcal, linkData, sortedPS, psElementType, reco::PFBlock::LINKTEST_ALL);
3440 
3441  // Loop over these PS clusters
3442  double totalPS = 0.;
3443  for (auto const& ps : sortedPS) {
3444  // CLuster index and distance to iEcal
3445  unsigned iPS = ps.second;
3446  // double distPS = ps.first;
3447 
3448  // Ignore clusters already in use
3449  if (!active[iPS])
3450  continue;
3451 
3452  // Check that this cluster is not closer to another ECAL cluster
3453  std::multimap<double, unsigned> sortedECAL;
3454  block.associatedElements(iPS, linkData, sortedECAL, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
3455  unsigned jEcal = sortedECAL.begin()->second;
3456  if (jEcal != iEcal)
3457  continue;
3458 
3459  // Update PS energy
3460  PFBlockElement::Type pstype = elements[iPS].type();
3461  assert(pstype == psElementType);
3462  PFClusterRef psclusterref = elements[iPS].clusterRef();
3463  assert(!psclusterref.isNull());
3464  totalPS += psclusterref->energy();
3465  psEne[0] += psclusterref->energy();
3466  active[iPS] = false;
3467  }
3468 }
assert(be >=bs)
bool isNull() const
Checks for null.
Definition: Ref.h:235

◆ checkAndReconstructSecondaryInteraction()

bool PFAlgo::checkAndReconstructSecondaryInteraction ( const reco::PFBlockRef blockref,
const edm::OwnVector< reco::PFBlockElement > &  elements,
bool  isActive,
int  iElement 
)
private

Definition at line 855 of file PFAlgo.cc.

References bookConverter::elements, isFromSecInt(), LogTrace, reconstructTrack(), and runTheMatrix::ret.

Referenced by elementLoop().

858  {
859  bool ret = isActive;
860  if (isActive && isFromSecInt(elements[iElement], "primary")) {
861  bool isPrimaryTrack =
862  elements[iElement].displacedVertexRef(PFBlockElement::T_TO_DISP)->displacedVertexRef()->isTherePrimaryTracks();
863  if (isPrimaryTrack) {
864  LogTrace("PFAlgo|elementLoop") << "Primary Track reconstructed alone";
865 
866  unsigned tmpi = reconstructTrack(elements[iElement]);
867  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iElement);
868  ret = false;
869  }
870  }
871 
872  return ret;
873 }
ret
prodAgent to be discontinued
bool isFromSecInt(const reco::PFBlockElement &eTrack, std::string order) const
Definition: PFAlgo.cc:3470
#define LogTrace(id)
unsigned reconstructTrack(const reco::PFBlockElement &elt, bool allowLoose=false)
Definition: PFAlgo.cc:3141

◆ checkCleaning()

void PFAlgo::checkCleaning ( const reco::PFRecHitCollection cleanedHF)

Check HF Cleaning.

Definition at line 3595 of file PFAlgo.cc.

References mps_fire::i, LogTrace, minDeltaMet_, GetRecoTauVFromDQM_MC_cff::next, pfCandidates_, DiDispStaMuonMonitor_cfi::pt, multPhiCorr_741_25nsDY_cfi::px, multPhiCorr_741_25nsDY_cfi::py, reconstructCluster(), optionsL1T::skip, mathSSE::sqrt(), hit::x, hit::y, and hit::z.

Referenced by PFProducer::produce().

3595  {
3596  // No hits to recover, leave.
3597  if (cleanedHits.empty())
3598  return;
3599 
3600  //Compute met and met significance (met/sqrt(SumEt))
3601  double metX = 0.;
3602  double metY = 0.;
3603  double sumet = 0;
3604  std::vector<unsigned int> hitsToBeAdded;
3605  for (auto const& pfc : *pfCandidates_) {
3606  metX += pfc.px();
3607  metY += pfc.py();
3608  sumet += pfc.pt();
3609  }
3610  double met2 = metX * metX + metY * metY;
3611  double met2_Original = met2;
3612  // Select events with large MET significance.
3613  // double significance = std::sqrt(met2/sumet);
3614  // double significanceCor = significance;
3615  double metXCor = metX;
3616  double metYCor = metY;
3617  double sumetCor = sumet;
3618  double met2Cor = met2;
3619  bool next = true;
3620  unsigned iCor = 1E9;
3621 
3622  // Find the cleaned hit with the largest effect on the MET
3623  while (next) {
3624  double metReduc = -1.;
3625  // Loop on the candidates
3626  for (unsigned i = 0; i < cleanedHits.size(); ++i) {
3627  const PFRecHit& hit = cleanedHits[i];
3628  double length = std::sqrt(hit.position().mag2());
3629  double px = hit.energy() * hit.position().x() / length;
3630  double py = hit.energy() * hit.position().y() / length;
3631  double pt = std::sqrt(px * px + py * py);
3632 
3633  // Check that it is not already scheculed to be cleaned
3634  bool skip = false;
3635  for (unsigned int hitIdx : hitsToBeAdded) {
3636  if (i == hitIdx)
3637  skip = true;
3638  if (skip)
3639  break;
3640  }
3641  if (skip)
3642  continue;
3643 
3644  // Now add the candidate to the MET
3645  double metXInt = metX + px;
3646  double metYInt = metY + py;
3647  double sumetInt = sumet + pt;
3648  double met2Int = metXInt * metXInt + metYInt * metYInt;
3649 
3650  // And check if it could contribute to a MET reduction
3651  if (met2Int < met2Cor) {
3652  metXCor = metXInt;
3653  metYCor = metYInt;
3654  metReduc = (met2 - met2Int) / met2Int;
3655  met2Cor = met2Int;
3656  sumetCor = sumetInt;
3657  // significanceCor = std::sqrt(met2Cor/sumetCor);
3658  iCor = i;
3659  }
3660  }
3661  //
3662  // If the MET must be significanly reduced, schedule the candidate to be added
3663  //
3664  if (metReduc > minDeltaMet_) {
3665  hitsToBeAdded.push_back(iCor);
3666  metX = metXCor;
3667  metY = metYCor;
3668  sumet = sumetCor;
3669  met2 = met2Cor;
3670  } else {
3671  // Otherwise just stop the loop
3672  next = false;
3673  }
3674  }
3675  //
3676  // At least 10 GeV MET reduction
3677  if (std::sqrt(met2_Original) - std::sqrt(met2) > 5.) {
3678  LogTrace("PFAlgo|checkCleaning") << hitsToBeAdded.size() << " hits were re-added ";
3679  LogTrace("PFAlgo|checkCleaning") << "MET reduction = " << std::sqrt(met2_Original) << " -> " << std::sqrt(met2Cor)
3680  << " = " << std::sqrt(met2Cor) - std::sqrt(met2_Original);
3681  LogTrace("PFAlgo|checkCleaning") << "Added after cleaning check : ";
3682  for (unsigned int hitIdx : hitsToBeAdded) {
3683  const PFRecHit& hit = cleanedHits[hitIdx];
3684  PFCluster cluster(hit.layer(), hit.energy(), hit.position().x(), hit.position().y(), hit.position().z());
3685  reconstructCluster(cluster, hit.energy());
3686  LogTrace("PFAlgo|checkCleaning") << pfCandidates_->back() << ". time = " << hit.time();
3687  }
3688  }
3689 } //PFAlgo::checkCleaning
Particle flow cluster, see clustering algorithm in PFClusterAlgo.
Definition: PFCluster.h:42
double minDeltaMet_
Definition: PFAlgo.h:324
std::unique_ptr< reco::PFCandidateCollection > pfCandidates_
Definition: PFAlgo.h:226
unsigned reconstructCluster(const reco::PFCluster &cluster, double particleEnergy, bool useDirection=false, double particleX=0., double particleY=0., double particleZ=0.)
Definition: PFAlgo.cc:3220
#define LogTrace(id)
Particle flow rechit (rechit + geometry and topology information). See clustering algorithm in PFClus...
Definition: PFRecHit.h:31
T sqrt(T t)
Definition: SSEVec.h:19

◆ checkGoodTrackDeadHcal()

bool PFAlgo::checkGoodTrackDeadHcal ( const reco::TrackRef trackRef,
bool  hasDeadHcal 
)
private

Definition at line 902 of file PFAlgo.cc.

References funct::abs(), goodPixelTrackDeadHcal_chi2n_, goodPixelTrackDeadHcal_dxy_, goodPixelTrackDeadHcal_dz_, goodPixelTrackDeadHcal_maxLost3Hit_, goodPixelTrackDeadHcal_maxLost4Hit_, goodPixelTrackDeadHcal_maxPt_, goodPixelTrackDeadHcal_minEta_, goodPixelTrackDeadHcal_ptErrRel_, goodTrackDeadHcal_chi2n_, goodTrackDeadHcal_dxy_, goodTrackDeadHcal_layers_, goodTrackDeadHcal_ptErrRel_, goodTrackDeadHcal_validFr_, LogTrace, reco::Vertex::position(), and primaryVertex_.

Referenced by elementLoop().

902  {
903  bool goodTrackDeadHcal = false;
904  if (hasDeadHcal) {
905  goodTrackDeadHcal = (trackRef->ptError() < goodTrackDeadHcal_ptErrRel_ * trackRef->pt() &&
906  trackRef->normalizedChi2() < goodTrackDeadHcal_chi2n_ &&
907  trackRef->hitPattern().trackerLayersWithMeasurement() >= goodTrackDeadHcal_layers_ &&
908  trackRef->validFraction() > goodTrackDeadHcal_validFr_ &&
910  // now we add an extra block for tracks at high |eta|
911  if (!goodTrackDeadHcal && std::abs(trackRef->eta()) > goodPixelTrackDeadHcal_minEta_ && // high eta
912  trackRef->hitPattern().pixelLayersWithMeasurement() >= 3 && // pixel track
913  trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::TRACK_HITS) == 0 &&
914  trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_INNER_HITS) == 0 &&
915  trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_OUTER_HITS) <=
916  (trackRef->hitPattern().pixelLayersWithMeasurement() > 3 ? goodPixelTrackDeadHcal_maxLost4Hit_
918  trackRef->normalizedChi2() < goodPixelTrackDeadHcal_chi2n_ && // tighter cut
921  trackRef->ptError() < goodPixelTrackDeadHcal_ptErrRel_ * trackRef->pt() && // sanity
922  trackRef->pt() < goodPixelTrackDeadHcal_maxPt_) { // sanity
923  goodTrackDeadHcal = true;
924  // FIXME: may decide to do something to the track pT
925  }
926  //if (!goodTrackDeadHcal && trackRef->hitPattern().trackerLayersWithMeasurement() == 4 && trackRef->validFraction() == 1
927  LogTrace("PFAlgo|elementLoop")
928  << " track pt " << trackRef->pt() << " +- " << trackRef->ptError() << " layers valid "
929  << trackRef->hitPattern().trackerLayersWithMeasurement() << ", lost "
930  << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::TRACK_HITS) << ", lost outer "
931  << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_OUTER_HITS) << ", lost inner "
932  << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_INNER_HITS) << "(valid fraction "
933  << trackRef->validFraction() << ")"
934  << " chi2/ndf " << trackRef->normalizedChi2() << " |dxy| " << std::abs(trackRef->dxy(primaryVertex_.position()))
935  << " +- " << trackRef->dxyError() << " |dz| " << std::abs(trackRef->dz(primaryVertex_.position())) << " +- "
936  << trackRef->dzError() << (goodTrackDeadHcal ? " passes " : " fails ") << "quality cuts";
937  }
938  return goodTrackDeadHcal;
939 }
float goodTrackDeadHcal_dxy_
Definition: PFAlgo.h:305
const Point & position() const
position
Definition: Vertex.h:127
float goodPixelTrackDeadHcal_minEta_
Definition: PFAlgo.h:307
float goodPixelTrackDeadHcal_dxy_
Definition: PFAlgo.h:313
#define LogTrace(id)
float goodTrackDeadHcal_validFr_
Definition: PFAlgo.h:304
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
reco::Vertex primaryVertex_
Definition: PFAlgo.h:328
int goodPixelTrackDeadHcal_maxLost3Hit_
Definition: PFAlgo.h:311
float goodPixelTrackDeadHcal_maxPt_
Definition: PFAlgo.h:308
int goodTrackDeadHcal_layers_
Definition: PFAlgo.h:303
float goodPixelTrackDeadHcal_ptErrRel_
Definition: PFAlgo.h:309
float goodTrackDeadHcal_ptErrRel_
Variables for track cleaning in bad HCal areas.
Definition: PFAlgo.h:301
int goodPixelTrackDeadHcal_maxLost4Hit_
Definition: PFAlgo.h:312
float goodPixelTrackDeadHcal_chi2n_
Definition: PFAlgo.h:310
float goodTrackDeadHcal_chi2n_
Definition: PFAlgo.h:302
float goodPixelTrackDeadHcal_dz_
Definition: PFAlgo.h:314

◆ checkHasDeadHcal()

bool PFAlgo::checkHasDeadHcal ( const std::multimap< double, unsigned > &  hcalElems,
const std::vector< bool > &  deadArea 
)
private

Definition at line 875 of file PFAlgo.cc.

Referenced by elementLoop().

875  {
876  // there's 3 possible options possible here, in principle:
877  // 1) flag everything that may be associated to a dead hcal marker
878  // 2) flag everything whose closest hcal link is a dead hcal marker
879  // 3) flag only things that are linked only to dead hcal marker
880  // in our first test we go for (2)
881  //--- option (1) --
882  //bool hasDeadHcal = false;
883  //for (auto it = hcalElems.begin(), ed = hcalElems.end(); it != ed; /*NOTE NO ++it HERE */ ) {
884  // if (deadArea[it->second]) { hasDeadHcal = true; it = hcalElems.erase(it); } // std::multimap::erase returns iterator to next
885  // else ++it;
886  //}
887  //--- option (2) --
888  bool hasDeadHcal = false;
889  if (!hcalElems.empty() && deadArea[hcalElems.begin()->second]) {
890  hasDeadHcal = true;
891  }
892  //--- option (3) --
893  //bool hasDeadHcal = true;
894  //for (auto it = hcalElems.begin(), ed = hcalElems.end(); it != ed; /*NOTE NO ++it HERE */ ) {
895  // if (deadArea[it->second]) { it = hcalElems.erase(it); } // std::multimap::erase returns iterator to next
896  // else { hasDeadHcal = false; }
897  //}
898  return hasDeadHcal;
899 }

◆ conversionAlgo()

void PFAlgo::conversionAlgo ( const edm::OwnVector< reco::PFBlockElement > &  elements,
std::vector< bool > &  active 
)
private

Definition at line 370 of file PFAlgo.cc.

References reco::TrackBase::conversionStep, bookConverter::elements, relativeConstraints::empty, reco::TrackBase::highPurity, LogTrace, quality, and reco::PFBlockElement::T_FROM_GAMMACONV.

Referenced by processBlock().

370  {
371  LogTrace("PFAlgo|conversionAlgo") << "start of function PFAlgo::conversionAlgo(), elements.size()="
372  << elements.size();
373  for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
374  PFBlockElement::Type type = elements[iEle].type();
375  if (type == PFBlockElement::TRACK) {
376  LogTrace("PFAlgo|conversionAlgo") << "elements[" << iEle << "].type() == TRACK, active[" << iEle
377  << "]=" << active[iEle];
378  if (elements[iEle].trackRef()->algo() == reco::TrackBase::conversionStep) {
379  active[iEle] = false;
380  }
381  if (elements[iEle].trackRef()->quality(reco::TrackBase::highPurity)) {
382  LogTrace("PFAlgo|conversionAlgo") << "Track is high purity";
383  continue;
384  }
385  const auto* trackRef = dynamic_cast<const reco::PFBlockElementTrack*>((&elements[iEle]));
386  if (!(trackRef->trackType(reco::PFBlockElement::T_FROM_GAMMACONV))) {
387  LogTrace("PFAlgo|conversionAlgo") << "!trackType(T_FROM_GAMMACONV)";
388  continue;
389  }
390  if (!elements[iEle].convRefs().empty()) {
391  active[iEle] = false;
392  }
393  LogTrace("PFAlgo|conversionAlgo") << "active[iEle=" << iEle << "]=" << active[iEle];
394  }
395  }
396  LogTrace("PFAlgo|conversionAlgo") << "end of function PFAlgo::conversionAlgo";
397 }
#define LogTrace(id)
string quality

◆ createCandidatesECAL()

void PFAlgo::createCandidatesECAL ( const reco::PFBlock block,
reco::PFBlock::LinkData linkData,
const edm::OwnVector< reco::PFBlockElement > &  elements,
std::vector< bool > &  active,
const reco::PFBlockRef blockref,
ElementIndices inds,
std::vector< bool > &  deadArea 
)
private

Definition at line 3012 of file PFAlgo.cc.

References cms::cuda::assert(), ECAL, ElementIndices::ecalIs, bookConverter::elements, mps_fire::i, edm::Ref< C, T, F >::isNull(), LogTrace, and reconstructCluster().

Referenced by processBlock().

3018  {
3019  LogTrace("PFAlgo|createCandidatesECAL")
3020  << "start of function PFAlgo::createCandidatesECAL(), ecalIs.size()=" << inds.ecalIs.size();
3021 
3022  // --------------- loop ecal ------------------
3023 
3024  // for each ecal element iEcal = ecalIs[i] in turn:
3025 
3026  for (unsigned i = 0; i < inds.ecalIs.size(); i++) {
3027  unsigned iEcal = inds.ecalIs[i];
3028 
3029  LogTrace("PFAlgo|createCandidatesECAL") << "elements[" << iEcal << "]=" << elements[iEcal];
3030 
3031  if (!active[iEcal]) {
3032  LogTrace("PFAlgo|createCandidatesECAL") << "iEcal=" << iEcal << " not active";
3033  continue;
3034  }
3035 
3036  PFBlockElement::Type type = elements[iEcal].type();
3038 
3039  PFClusterRef clusterref = elements[iEcal].clusterRef();
3040  assert(!clusterref.isNull());
3041 
3042  active[iEcal] = false;
3043 
3044  float ecalEnergy = clusterref->correctedEnergy();
3045  // float ecalEnergy = calibration_.energyEm( clusterref->energy() );
3046  double particleEnergy = ecalEnergy;
3047 
3048  auto& cand = (*pfCandidates_)[reconstructCluster(*clusterref, particleEnergy)];
3049 
3050  cand.setEcalEnergy(clusterref->energy(), ecalEnergy);
3051  cand.setHcalEnergy(0., 0.);
3052  cand.setHoEnergy(0., 0.);
3053  cand.setPs1Energy(0.);
3054  cand.setPs2Energy(0.);
3055  cand.addElementInBlock(blockref, iEcal);
3056 
3057  } // end loop on ecal elements iEcal = ecalIs[i]
3058  LogTrace("PFAlgo|createCandidatesECAL") << "end of function PFALgo::createCandidatesECAL";
3059 }
unsigned reconstructCluster(const reco::PFCluster &cluster, double particleEnergy, bool useDirection=false, double particleX=0., double particleY=0., double particleZ=0.)
Definition: PFAlgo.cc:3220
assert(be >=bs)
#define LogTrace(id)
bool isNull() const
Checks for null.
Definition: Ref.h:235
std::vector< unsigned > ecalIs
Definition: PFAlgo.h:44

◆ createCandidatesHCAL()

void PFAlgo::createCandidatesHCAL ( const reco::PFBlock block,
reco::PFBlock::LinkData linkData,
const edm::OwnVector< reco::PFBlockElement > &  elements,
std::vector< bool > &  active,
const reco::PFBlockRef blockref,
ElementIndices inds,
std::vector< bool > &  deadArea 
)
private

Definition at line 1702 of file PFAlgo.cc.

References a, funct::abs(), cms::cuda::assert(), b, groupFilesInBlocks::block, calibration_, HPSPFTauProducerPuppi_cfi::chargedHadron, RPCNoise_example::check, Calorimetry_cff::dp, dptRel_DispVtx_, MillePedeFileConverter_cfg::e, ECAL, reco::PFBlockElement::ECAL, RecoEcal_cff::ecalClusters, bookConverter::elements, PFEnergyCalibration::energyEmHad(), PFTrackAlgoTools::errorScale(), factors45_, first, HLT_2022v12_cff::fraction, HCAL, ElementIndices::hcalIs, photonIsolationHIProducer_cfi::ho, DigiToRawDM_cff::HO, reco::PFBlockElement::HO, mps_fire::i, cuy::ii, isFromSecInt(), PFMuonAlgo::isIsolatedMuon(), PFMuonAlgo::isLooseMuon(), PFMuonAlgo::isMuon(), edm::Ref< C, T, F >::isNull(), phase2tkutil::isPrimary(), dqmiolumiharvest::j, reco::PFBlock::LINKTEST_ALL, LogTrace, SiStripPI::max, SiStripPI::min, muonECAL_, muonHCAL_, muonHO_, neutralHadronEnergyResolution(), custom_jme_cff::nMuons, nSigmaHCAL(), nSigmaTRACK_, convertSQLiteXML::ok, AlCaHLTBitMon_ParallelJobs::p, pfCandidates_, ptError_, bTagCommon_cff::ratioMax, reconstructCluster(), reconstructTrack(), rejectTracks_Bad_, rejectTracks_Step45_, edm::second(), setHcalDepthInfo(), reco::PFBlockElementTrack::setPositionAtECALEntrance(), mathSSE::sqrt(), PFTrackAlgoTools::step5(), reco::PFBlockElement::TRACK, reco::btau::trackMomentum, useHO_, x, X, BeamSpotPI::Y, and BeamSpotPI::Z.

Referenced by processBlock().

1708  {
1709  LogTrace("PFAlgo|createCandidatesHCAL")
1710  << "start of function PFAlgo::createCandidatesHCAL, inds.hcalIs.size()=" << inds.hcalIs.size();
1711 
1712  // --------------- loop hcal ------------------
1713 
1714  for (unsigned iHcal : inds.hcalIs) {
1715  PFBlockElement::Type type = elements[iHcal].type();
1716 
1718 
1719  LogTrace("PFAlgo|createCandidatesHCAL") << "elements[" << iHcal << "]=" << elements[iHcal];
1720 
1721  // associated tracks
1722  std::multimap<double, unsigned> sortedTracks;
1723  block.associatedElements(iHcal, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
1724 
1725  std::multimap<unsigned, std::pair<double, unsigned>> associatedEcals;
1726 
1727  std::map<unsigned, std::pair<double, double>> associatedPSs;
1728 
1729  std::multimap<double, std::pair<unsigned, bool>> associatedTracks;
1730 
1731  // A temporary maps for ECAL satellite clusters
1732  std::multimap<double, std::tuple<unsigned, ::math::XYZVector, double>>
1733  ecalSatellites; // last element (double) : correction under the egamma hypothesis
1734  std::tuple<unsigned, ::math::XYZVector, double> fakeSatellite(iHcal, ::math::XYZVector(0., 0., 0.), 1.);
1735  ecalSatellites.emplace(-1., fakeSatellite);
1736 
1737  std::multimap<unsigned, std::pair<double, unsigned>> associatedHOs;
1738 
1739  PFClusterRef hclusterref = elements[iHcal].clusterRef();
1740  assert(!hclusterref.isNull());
1741 
1742  //if there is no track attached to that HCAL, then do not
1743  //reconstruct an HCAL alone, check if it can be recovered
1744  //first
1745 
1746  // if no associated tracks, create a neutral hadron
1747  //if(sortedTracks.empty() ) {
1748 
1749  if (sortedTracks.empty()) {
1750  LogTrace("PFAlgo|createCandidatesHCAL") << "\tno associated tracks, keep for later";
1751  continue;
1752  }
1753 
1754  // Lock the HCAL cluster
1755  active[iHcal] = false;
1756 
1757  // in the following, tracks are associated to this hcal cluster.
1758  // will look for an excess of energy in the calorimeters w/r to
1759  // the charged energy, and turn this excess into a neutral hadron or
1760  // a photon.
1761 
1762  LogTrace("PFAlgo|createCandidatesHCAL") << sortedTracks.size() << " associated tracks:";
1763 
1764  double totalChargedMomentum = 0;
1765  double sumpError2 = 0.;
1766  double totalHO = 0.;
1767  double totalEcal = 0.;
1768  double totalEcalEGMCalib = 0.;
1769  double totalHcal = hclusterref->energy();
1770  vector<double> hcalP;
1771  vector<double> hcalDP;
1772  vector<unsigned> tkIs;
1773  double maxDPovP = -9999.;
1774 
1775  //Keep track of how much energy is assigned to calorimeter-vs-track energy/momentum excess
1776  vector<unsigned> chargedHadronsIndices;
1777  vector<unsigned> chargedHadronsInBlock;
1778  double mergedNeutralHadronEnergy = 0;
1779  double mergedPhotonEnergy = 0;
1780  double muonHCALEnergy = 0.;
1781  double muonECALEnergy = 0.;
1782  double muonHCALError = 0.;
1783  double muonECALError = 0.;
1784  unsigned nMuons = 0;
1785 
1786  ::math::XYZVector photonAtECAL(0., 0., 0.);
1787  std::vector<std::tuple<unsigned, ::math::XYZVector, double>>
1788  ecalClusters; // last element (double) : correction under the egamma hypothesis
1789  double sumEcalClusters = 0;
1790  ::math::XYZVector hadronDirection(
1791  hclusterref->position().X(), hclusterref->position().Y(), hclusterref->position().Z());
1792  hadronDirection = hadronDirection.Unit();
1793  ::math::XYZVector hadronAtECAL = totalHcal * hadronDirection;
1794 
1795  // Loop over all tracks associated to this HCAL cluster
1796  for (auto const& trk : sortedTracks) {
1797  unsigned iTrack = trk.second;
1798 
1799  // Check the track has not already been used (e.g., in electrons, conversions...)
1800  if (!active[iTrack])
1801  continue;
1802  // Sanity check 1
1803  PFBlockElement::Type type = elements[iTrack].type();
1805  // Sanity check 2
1806  reco::TrackRef trackRef = elements[iTrack].trackRef();
1807  assert(!trackRef.isNull());
1808 
1809  // The direction at ECAL entrance
1810  const ::math::XYZPointF& chargedPosition =
1811  dynamic_cast<const reco::PFBlockElementTrack*>(&elements[iTrack])->positionAtECALEntrance();
1812  ::math::XYZVector chargedDirection(chargedPosition.X(), chargedPosition.Y(), chargedPosition.Z());
1813  chargedDirection = chargedDirection.Unit();
1814 
1815  // look for ECAL elements associated to iTrack (associated to iHcal)
1816  std::multimap<double, unsigned> sortedEcals;
1817  block.associatedElements(iTrack, linkData, sortedEcals, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
1818 
1819  LogTrace("PFAlgo|createCandidatesHCAL") << "number of Ecal elements linked to this track: " << sortedEcals.size();
1820 
1821  // look for HO elements associated to iTrack (associated to iHcal)
1822  std::multimap<double, unsigned> sortedHOs;
1823  if (useHO_) {
1824  block.associatedElements(iTrack, linkData, sortedHOs, reco::PFBlockElement::HO, reco::PFBlock::LINKTEST_ALL);
1825  }
1826  LogTrace("PFAlgo|createCandidatesHCAL")
1827  << "PFAlgo : number of HO elements linked to this track: " << sortedHOs.size();
1828 
1829  // Create a PF Candidate right away if the track is a tight muon
1830  reco::MuonRef muonRef = elements[iTrack].muonRef();
1831 
1832  bool thisIsAMuon = PFMuonAlgo::isMuon(elements[iTrack]);
1833  bool thisIsAnIsolatedMuon = PFMuonAlgo::isIsolatedMuon(elements[iTrack]);
1834  bool thisIsALooseMuon = false;
1835 
1836  if (!thisIsAMuon) {
1837  thisIsALooseMuon = PFMuonAlgo::isLooseMuon(elements[iTrack]);
1838  }
1839 
1840  if (thisIsAMuon) {
1841  LogTrace("PFAlgo|createCandidatesHCAL") << "This track is identified as a muon - remove it from the stack";
1842  LogTrace("PFAlgo|createCandidatesHCAL") << elements[iTrack];
1843 
1844  // Create a muon.
1845 
1846  unsigned tmpi = reconstructTrack(elements[iTrack]);
1847 
1848  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iTrack);
1849  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHcal);
1850  double muonHcal = std::min(muonHCAL_[0] + muonHCAL_[1], totalHcal);
1851 
1852  // if muon is isolated and muon momentum exceeds the calo energy, absorb the calo energy
1853  // rationale : there has been a EM showering by the muon in the calorimeter - or the coil -
1854  // and we don't want to double count the energy
1855  bool letMuonEatCaloEnergy = false;
1856 
1857  if (thisIsAnIsolatedMuon) {
1858  // The factor 1.3 is the e/pi factor in HCAL...
1859  double totalCaloEnergy = totalHcal / 1.30;
1860  unsigned iEcal = 0;
1861  if (!sortedEcals.empty()) {
1862  iEcal = sortedEcals.begin()->second;
1863  PFClusterRef eclusterref = elements[iEcal].clusterRef();
1864  totalCaloEnergy += eclusterref->energy();
1865  }
1866 
1867  if (useHO_) {
1868  // The factor 1.3 is assumed to be the e/pi factor in HO, too.
1869  unsigned iHO = 0;
1870  if (!sortedHOs.empty()) {
1871  iHO = sortedHOs.begin()->second;
1872  PFClusterRef eclusterref = elements[iHO].clusterRef();
1873  totalCaloEnergy += eclusterref->energy() / 1.30;
1874  }
1875  }
1876 
1877  if ((pfCandidates_->back()).p() > totalCaloEnergy)
1878  letMuonEatCaloEnergy = true;
1879  }
1880 
1881  if (letMuonEatCaloEnergy)
1882  muonHcal = totalHcal;
1883  double muonEcal = 0.;
1884  unsigned iEcal = 0;
1885  if (!sortedEcals.empty()) {
1886  iEcal = sortedEcals.begin()->second;
1887  PFClusterRef eclusterref = elements[iEcal].clusterRef();
1888  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iEcal);
1889  muonEcal = std::min(muonECAL_[0] + muonECAL_[1], eclusterref->energy());
1890  if (letMuonEatCaloEnergy)
1891  muonEcal = eclusterref->energy();
1892  // If the muon expected energy accounts for the whole ecal cluster energy, lock the ecal cluster
1893  if (eclusterref->energy() - muonEcal < 0.2)
1894  active[iEcal] = false;
1895  (*pfCandidates_)[tmpi].setEcalEnergy(eclusterref->energy(), muonEcal);
1896  }
1897  unsigned iHO = 0;
1898  double muonHO = 0.;
1899  if (useHO_) {
1900  if (!sortedHOs.empty()) {
1901  iHO = sortedHOs.begin()->second;
1902  PFClusterRef hoclusterref = elements[iHO].clusterRef();
1903  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHO);
1904  muonHO = std::min(muonHO_[0] + muonHO_[1], hoclusterref->energy());
1905  if (letMuonEatCaloEnergy)
1906  muonHO = hoclusterref->energy();
1907  // If the muon expected energy accounts for the whole HO cluster energy, lock the HO cluster
1908  if (hoclusterref->energy() - muonHO < 0.2)
1909  active[iHO] = false;
1910  (*pfCandidates_)[tmpi].setHcalEnergy(totalHcal, muonHcal);
1911  (*pfCandidates_)[tmpi].setHoEnergy(hoclusterref->energy(), muonHO);
1912  }
1913  } else {
1914  (*pfCandidates_)[tmpi].setHcalEnergy(totalHcal, muonHcal);
1915  }
1916  setHcalDepthInfo((*pfCandidates_)[tmpi], *hclusterref);
1917 
1918  if (letMuonEatCaloEnergy) {
1919  muonHCALEnergy += totalHcal;
1920  if (useHO_)
1921  muonHCALEnergy += muonHO;
1922  muonHCALError += 0.;
1923  muonECALEnergy += muonEcal;
1924  muonECALError += 0.;
1925  photonAtECAL -= muonEcal * chargedDirection;
1926  hadronAtECAL -= totalHcal * chargedDirection;
1927  if (!sortedEcals.empty())
1928  active[iEcal] = false;
1929  active[iHcal] = false;
1930  if (useHO_ && !sortedHOs.empty())
1931  active[iHO] = false;
1932  } else {
1933  // Estimate of the energy deposit & resolution in the calorimeters
1934  muonHCALEnergy += muonHCAL_[0];
1935  muonHCALError += muonHCAL_[1] * muonHCAL_[1];
1936  if (muonHO > 0.) {
1937  muonHCALEnergy += muonHO_[0];
1938  muonHCALError += muonHO_[1] * muonHO_[1];
1939  }
1940  muonECALEnergy += muonECAL_[0];
1941  muonECALError += muonECAL_[1] * muonECAL_[1];
1942  // ... as well as the equivalent "momentum" at ECAL entrance
1943  photonAtECAL -= muonECAL_[0] * chargedDirection;
1944  hadronAtECAL -= muonHCAL_[0] * chargedDirection;
1945  }
1946 
1947  // Remove it from the stack
1948  active[iTrack] = false;
1949  // Go to next track
1950  continue;
1951  }
1952 
1953  //
1954 
1955  LogTrace("PFAlgo|createCandidatesHCAL") << "elements[iTrack=" << iTrack << "]=" << elements[iTrack];
1956 
1957  // introduce tracking errors
1958  double trackMomentum = trackRef->p();
1959  totalChargedMomentum += trackMomentum;
1960 
1961  // If the track is not a tight muon, but still resembles a muon
1962  // keep it for later for blocks with too large a charged energy
1963  if (thisIsALooseMuon && !thisIsAMuon)
1964  nMuons += 1;
1965 
1966  // ... and keep anyway the pt error for possible fake rejection
1967  // ... blow up errors of 5th and 4th iteration, to reject those
1968  // ... tracks first (in case it's needed)
1969  double dpt = trackRef->ptError();
1970  double blowError = PFTrackAlgoTools::errorScale(trackRef->algo(), factors45_);
1971  // except if it is from an interaction
1972  bool isPrimaryOrSecondary = isFromSecInt(elements[iTrack], "all");
1973 
1974  if (isPrimaryOrSecondary)
1975  blowError = 1.;
1976 
1977  std::pair<unsigned, bool> tkmuon(iTrack, thisIsALooseMuon);
1978  associatedTracks.emplace(-dpt * blowError, tkmuon);
1979 
1980  // Also keep the total track momentum error for comparison with the calo energy
1981  double dp = trackRef->qoverpError() * trackMomentum * trackMomentum;
1982  sumpError2 += dp * dp;
1983 
1984  bool connectedToEcal = false; // Will become true if there is at least one ECAL cluster connected
1985  if (!sortedEcals.empty()) { // start case: at least one ecal element associated to iTrack
1986 
1987  // Loop over all connected ECAL clusters
1988  for (auto const& ecal : sortedEcals) {
1989  unsigned iEcal = ecal.second;
1990  double dist = ecal.first;
1991 
1992  // Ignore ECAL cluters already used
1993  if (!active[iEcal]) {
1994  LogTrace("PFAlgo|createCandidatesHCAL") << "cluster locked";
1995  continue;
1996  }
1997 
1998  // Sanity checks
1999  PFBlockElement::Type type = elements[iEcal].type();
2001  PFClusterRef eclusterref = elements[iEcal].clusterRef();
2002  assert(!eclusterref.isNull());
2003 
2004  // Check if this ECAL is not closer to another track - ignore it in that case
2005  std::multimap<double, unsigned> sortedTracksEcal;
2006  block.associatedElements(
2007  iEcal, linkData, sortedTracksEcal, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
2008  unsigned jTrack = sortedTracksEcal.begin()->second;
2009  if (jTrack != iTrack)
2010  continue;
2011 
2012  double distEcal = block.dist(jTrack, iEcal, linkData, reco::PFBlock::LINKTEST_ALL);
2013 
2014  float ecalEnergyCalibrated = eclusterref->correctedEnergy(); // calibrated based on the egamma hypothesis
2015  float ecalEnergy = eclusterref->energy();
2016  ::math::XYZVector photonDirection(
2017  eclusterref->position().X(), eclusterref->position().Y(), eclusterref->position().Z());
2018  photonDirection = photonDirection.Unit();
2019 
2020  if (!connectedToEcal) { // This is the closest ECAL cluster - will add its energy later
2021 
2022  LogTrace("PFAlgo|createCandidatesHCAL") << "closest: " << elements[iEcal];
2023 
2024  connectedToEcal = true;
2025  // PJ 1st-April-09 : To be done somewhere !!! (Had to comment it, but it is needed)
2026  // currentChargedHadron.addElementInBlock( blockref, iEcal );
2027 
2028  // KH: we don't know if this satellite is due to egamma or hadron shower. use raw energy for PF hadron calibration._ store also calibration constant.
2029  double ecalCalibFactor = (ecalEnergy > 1E-9) ? ecalEnergyCalibrated / ecalEnergy : 1.;
2030  std::tuple<unsigned, ::math::XYZVector, double> satellite(
2031  iEcal, ecalEnergy * photonDirection, ecalCalibFactor);
2032  ecalSatellites.emplace(-1., satellite);
2033 
2034  } else { // Keep satellite clusters for later
2035 
2036  // KH: same as above.
2037  double ecalCalibFactor = (ecalEnergy > 1E-9) ? ecalEnergyCalibrated / ecalEnergy : 1.;
2038  std::tuple<unsigned, ::math::XYZVector, double> satellite(
2039  iEcal, ecalEnergy * photonDirection, ecalCalibFactor);
2040  ecalSatellites.emplace(dist, satellite);
2041  }
2042 
2043  std::pair<double, unsigned> associatedEcal(distEcal, iEcal);
2044  associatedEcals.emplace(iTrack, associatedEcal);
2045 
2046  } // End loop ecal associated to iTrack
2047  } // end case: at least one ecal element associated to iTrack
2048 
2049  if (useHO_ && !sortedHOs.empty()) { // start case: at least one ho element associated to iTrack
2050 
2051  // Loop over all connected HO clusters
2052  for (auto const& ho : sortedHOs) {
2053  unsigned iHO = ho.second;
2054  double distHO = ho.first;
2055 
2056  // Ignore HO cluters already used
2057  if (!active[iHO]) {
2058  LogTrace("PFAlgo|createCandidatesHCAL") << "cluster locked";
2059  continue;
2060  }
2061 
2062  // Sanity checks
2063  PFBlockElement::Type type = elements[iHO].type();
2065  PFClusterRef hoclusterref = elements[iHO].clusterRef();
2066  assert(!hoclusterref.isNull());
2067 
2068  // Check if this HO is not closer to another track - ignore it in that case
2069  std::multimap<double, unsigned> sortedTracksHO;
2070  block.associatedElements(
2071  iHO, linkData, sortedTracksHO, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
2072  unsigned jTrack = sortedTracksHO.begin()->second;
2073  if (jTrack != iTrack)
2074  continue;
2075 
2076  // double chi2HO = block.chi2(jTrack,iHO,linkData,
2077  // reco::PFBlock::LINKTEST_ALL);
2078  //double distHO = block.dist(jTrack,iHO,linkData,
2079  // reco::PFBlock::LINKTEST_ALL);
2080 
2081  // Increment the total energy by the energy of this HO cluster
2082  totalHO += hoclusterref->energy();
2083  active[iHO] = false;
2084  // Keep track for later reference in the PFCandidate.
2085  std::pair<double, unsigned> associatedHO(distHO, iHO);
2086  associatedHOs.emplace(iTrack, associatedHO);
2087 
2088  } // End loop ho associated to iTrack
2089  } // end case: at least one ho element associated to iTrack
2090 
2091  } // end loop on tracks associated to hcal element iHcal
2092 
2093  // Include totalHO in totalHCAL for the time being (it will be calibrated as HCAL energy)
2094  totalHcal += totalHO;
2095 
2096  // test compatibility between calo and tracker. //////////////
2097 
2098  double caloEnergy = 0.;
2099  double slopeEcal = 1.0;
2100  double calibEcal = 0.;
2101  double calibHcal = 0.;
2102  hadronDirection = hadronAtECAL.Unit();
2103 
2104  // Determine the expected calo resolution from the total charged momentum
2105  double caloResolution = neutralHadronEnergyResolution(totalChargedMomentum, hclusterref->positionREP().Eta());
2106  caloResolution *= totalChargedMomentum;
2107  // Account for muons
2108  caloResolution = std::sqrt(caloResolution * caloResolution + muonHCALError + muonECALError);
2109  totalEcal -= std::min(totalEcal, muonECALEnergy);
2110  totalEcalEGMCalib -= std::min(totalEcalEGMCalib, muonECALEnergy);
2111  totalHcal -= std::min(totalHcal, muonHCALEnergy);
2112  if (totalEcal < 1E-9)
2113  photonAtECAL = ::math::XYZVector(0., 0., 0.);
2114  if (totalHcal < 1E-9)
2115  hadronAtECAL = ::math::XYZVector(0., 0., 0.);
2116 
2117  // Loop over all ECAL satellites, starting for the closest to the various tracks
2118  // and adding other satellites until saturation of the total track momentum
2119  // Note : for code simplicity, the first element of the loop is the HCAL cluster
2120  // with 0 energy in the ECAL
2121  for (auto const& ecalSatellite : ecalSatellites) {
2122  // Add the energy of this ECAL cluster
2123  double previousCalibEcal = calibEcal;
2124  double previousCalibHcal = calibHcal;
2125  double previousCaloEnergy = caloEnergy;
2126  double previousSlopeEcal = slopeEcal;
2127  ::math::XYZVector previousHadronAtECAL = hadronAtECAL;
2128  //
2129  totalEcal +=
2130  sqrt(std::get<1>(ecalSatellite.second).Mag2()); // KH: raw ECAL energy for input to PF hadron calibration
2131  totalEcalEGMCalib += sqrt(std::get<1>(ecalSatellite.second).Mag2()) *
2132  std::get<2>(ecalSatellite.second); // KH: calibrated ECAL energy under the egamma hypothesis
2133  photonAtECAL += std::get<1>(ecalSatellite.second) *
2134  std::get<2>(ecalSatellite.second); // KH: calibrated ECAL energy under the egamma hypothesis
2135  calibEcal = std::max(0., totalEcal); // KH: preparing for hadron calibration
2136  calibHcal = std::max(0., totalHcal);
2137  hadronAtECAL = calibHcal * hadronDirection;
2138  // Calibrate ECAL and HCAL energy under the hadron hypothesis.
2139  calibration_.energyEmHad(totalChargedMomentum,
2140  calibEcal,
2141  calibHcal,
2142  hclusterref->positionREP().Eta(),
2143  hclusterref->positionREP().Phi());
2144  caloEnergy = calibEcal + calibHcal;
2145  if (totalEcal > 0.)
2146  slopeEcal = calibEcal / totalEcal;
2147 
2148  hadronAtECAL = calibHcal * hadronDirection;
2149 
2150  // Continue looping until all closest clusters are exhausted and as long as
2151  // the calorimetric energy does not saturate the total momentum.
2152  if (ecalSatellite.first < 0. || caloEnergy - totalChargedMomentum <= 0.) {
2153  LogTrace("PFAlgo|createCandidatesHCAL")
2154  << "\t\t\tactive, adding " << std::get<1>(ecalSatellite.second) << " to ECAL energy, and locking";
2155  active[std::get<0>(ecalSatellite.second)] = false;
2156  double clusterEnergy =
2157  sqrt(std::get<1>(ecalSatellite.second).Mag2()) *
2158  std::get<2>(ecalSatellite.second); // KH: ECAL energy calibrated under the egamma hypothesis
2159  if (clusterEnergy > 50) { // KH: used to split energetic ecal clusters (E>50 GeV)
2160  ecalClusters.push_back(ecalSatellite.second);
2161  sumEcalClusters += clusterEnergy;
2162  }
2163  continue;
2164  }
2165 
2166  // Otherwise, do not consider the last cluster examined and exit.
2167  // active[is->second.first] = true;
2168  totalEcal -= sqrt(std::get<1>(ecalSatellite.second).Mag2());
2169  totalEcalEGMCalib -= sqrt(std::get<1>(ecalSatellite.second).Mag2()) * std::get<2>(ecalSatellite.second);
2170  photonAtECAL -= std::get<1>(ecalSatellite.second) * std::get<2>(ecalSatellite.second);
2171  calibEcal = previousCalibEcal;
2172  calibHcal = previousCalibHcal;
2173  hadronAtECAL = previousHadronAtECAL;
2174  slopeEcal = previousSlopeEcal;
2175  caloEnergy = previousCaloEnergy;
2176 
2177  break;
2178  }
2179 
2180  // Sanity check !
2181  assert(caloEnergy >= 0);
2182 
2183  // And now check for hadronic energy excess...
2184 
2185  //colin: resolution should be measured on the ecal+hcal case.
2186  // however, the result will be close.
2187  // double caloResolution = neutralHadronEnergyResolution( caloEnergy );
2188  // caloResolution *= caloEnergy;
2189  // PJ The resolution is on the expected charged calo energy !
2190  //double caloResolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta());
2191  //caloResolution *= totalChargedMomentum;
2192  // that of the charged particles linked to the cluster!
2193 
2195  if (totalChargedMomentum - caloEnergy > nSigmaTRACK_ * caloResolution) {
2196  // First consider loose muons
2197  if (nMuons > 0) {
2198  for (auto const& trk : associatedTracks) {
2199  // Only muons
2200  if (!trk.second.second)
2201  continue;
2202 
2203  const unsigned int iTrack = trk.second.first;
2204  // Only active tracks
2205  if (!active[iTrack])
2206  continue;
2207 
2208  const double trackMomentum = elements[trk.second.first].trackRef()->p();
2209 
2210  // look for ECAL elements associated to iTrack (associated to iHcal)
2211  std::multimap<double, unsigned> sortedEcals;
2212  block.associatedElements(
2213  iTrack, linkData, sortedEcals, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
2214  std::multimap<double, unsigned> sortedHOs;
2215  block.associatedElements(iTrack, linkData, sortedHOs, reco::PFBlockElement::HO, reco::PFBlock::LINKTEST_ALL);
2216 
2217  //Here allow for loose muons!
2218  auto& muon = (*pfCandidates_)[reconstructTrack(elements[iTrack], true)];
2219 
2220  muon.addElementInBlock(blockref, iTrack);
2221  muon.addElementInBlock(blockref, iHcal);
2222  const double muonHcal = std::min(muonHCAL_[0] + muonHCAL_[1], totalHcal - totalHO);
2223  double muonHO = 0.;
2224  muon.setHcalEnergy(totalHcal, muonHcal);
2225  if (!sortedEcals.empty()) {
2226  const unsigned int iEcal = sortedEcals.begin()->second;
2227  PFClusterRef eclusterref = elements[iEcal].clusterRef();
2228  muon.addElementInBlock(blockref, iEcal);
2229  const double muonEcal = std::min(muonECAL_[0] + muonECAL_[1], eclusterref->energy());
2230  muon.setEcalEnergy(eclusterref->energy(), muonEcal);
2231  }
2232  if (useHO_ && !sortedHOs.empty()) {
2233  const unsigned int iHO = sortedHOs.begin()->second;
2234  PFClusterRef hoclusterref = elements[iHO].clusterRef();
2235  muon.addElementInBlock(blockref, iHO);
2236  muonHO = std::min(muonHO_[0] + muonHO_[1], hoclusterref->energy());
2237  muon.setHcalEnergy(max(totalHcal - totalHO, 0.0), muonHcal);
2238  muon.setHoEnergy(hoclusterref->energy(), muonHO);
2239  }
2240  setHcalDepthInfo(muon, *hclusterref);
2241  // Remove it from the block
2242  const ::math::XYZPointF& chargedPosition =
2243  dynamic_cast<const reco::PFBlockElementTrack*>(&elements[trk.second.first])->positionAtECALEntrance();
2244  ::math::XYZVector chargedDirection(chargedPosition.X(), chargedPosition.Y(), chargedPosition.Z());
2245  chargedDirection = chargedDirection.Unit();
2246  totalChargedMomentum -= trackMomentum;
2247  // Update the calo energies
2248  if (totalEcal > 0.)
2249  calibEcal -= std::min(calibEcal, muonECAL_[0] * calibEcal / totalEcal);
2250  if (totalHcal > 0.)
2251  calibHcal -= std::min(calibHcal, muonHCAL_[0] * calibHcal / totalHcal);
2252  totalEcal -= std::min(totalEcal, muonECAL_[0]);
2253  totalHcal -= std::min(totalHcal, muonHCAL_[0]);
2254  if (totalEcal > muonECAL_[0])
2255  photonAtECAL -= muonECAL_[0] * chargedDirection;
2256  if (totalHcal > muonHCAL_[0])
2257  hadronAtECAL -= muonHCAL_[0] * calibHcal / totalHcal * chargedDirection;
2258  caloEnergy = calibEcal + calibHcal;
2259  muonHCALEnergy += muonHCAL_[0];
2260  muonHCALError += muonHCAL_[1] * muonHCAL_[1];
2261  if (muonHO > 0.) {
2262  muonHCALEnergy += muonHO_[0];
2263  muonHCALError += muonHO_[1] * muonHO_[1];
2264  if (totalHcal > 0.) {
2265  calibHcal -= std::min(calibHcal, muonHO_[0] * calibHcal / totalHcal);
2266  totalHcal -= std::min(totalHcal, muonHO_[0]);
2267  }
2268  }
2269  muonECALEnergy += muonECAL_[0];
2270  muonECALError += muonECAL_[1] * muonECAL_[1];
2271  active[iTrack] = false;
2272  // Stop the loop whenever enough muons are removed
2273  //Commented out: Keep looking for muons since they often come in pairs -Matt
2274  //if ( totalChargedMomentum < caloEnergy ) break;
2275  }
2276  // New calo resolution.
2277  caloResolution = neutralHadronEnergyResolution(totalChargedMomentum, hclusterref->positionREP().Eta());
2278  caloResolution *= totalChargedMomentum;
2279  caloResolution = std::sqrt(caloResolution * caloResolution + muonHCALError + muonECALError);
2280  }
2281  }
2282 
2283 #ifdef EDM_ML_DEBUG
2284  LogTrace("PFAlgo|createCandidatesHCAL") << "\tBefore Cleaning ";
2285  LogTrace("PFAlgo|createCandidatesHCAL") << "\tCompare Calo Energy to total charged momentum ";
2286  LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tsum p = " << totalChargedMomentum << " +- " << sqrt(sumpError2);
2287  LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tsum ecal = " << totalEcal;
2288  LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tsum hcal = " << totalHcal;
2289  LogTrace("PFAlgo|createCandidatesHCAL") << "\t\t => Calo Energy = " << caloEnergy << " +- " << caloResolution;
2290  LogTrace("PFAlgo|createCandidatesHCAL")
2291  << "\t\t => Calo Energy- total charged momentum = " << caloEnergy - totalChargedMomentum << " +- "
2292  << sqrt(sumpError2 + caloResolution * caloResolution);
2293 #endif
2294 
2295  // Second consider bad tracks (if still needed after muon removal)
2296  unsigned corrTrack = 10000000;
2297  double corrFact = 1.;
2298 
2299  if (rejectTracks_Bad_ && totalChargedMomentum - caloEnergy > nSigmaTRACK_ * caloResolution) {
2300  for (auto const& trk : associatedTracks) {
2301  const unsigned iTrack = trk.second.first;
2302  // Only active tracks
2303  if (!active[iTrack])
2304  continue;
2305  const reco::TrackRef& trackref = elements[trk.second.first].trackRef();
2306 
2307  const double dptRel = fabs(trk.first) / trackref->pt() * 100;
2308  const bool isSecondary = isFromSecInt(elements[iTrack], "secondary");
2309  const bool isPrimary = isFromSecInt(elements[iTrack], "primary");
2310 
2311  if (isSecondary && dptRel < dptRel_DispVtx_)
2312  continue;
2313  // Consider only bad tracks
2314  if (fabs(trk.first) < ptError_)
2315  break;
2316  // What would become the block charged momentum if this track were removed
2317  const double wouldBeTotalChargedMomentum = totalChargedMomentum - trackref->p();
2318  // Reject worst tracks, as long as the total charged momentum
2319  // is larger than the calo energy
2320 
2321  if (wouldBeTotalChargedMomentum > caloEnergy) {
2322  if (isSecondary) {
2323  LogTrace("PFAlgo|createCandidatesHCAL")
2324  << "In bad track rejection step dptRel = " << dptRel << " dptRel_DispVtx_ = " << dptRel_DispVtx_;
2325  LogTrace("PFAlgo|createCandidatesHCAL")
2326  << "The calo energy would be still smaller even without this track but it is attached to a NI";
2327  }
2328 
2329  if (isPrimary || (isSecondary && dptRel < dptRel_DispVtx_))
2330  continue;
2331  active[iTrack] = false;
2332  totalChargedMomentum = wouldBeTotalChargedMomentum;
2333  LogTrace("PFAlgo|createCandidatesHCAL")
2334  << "\tElement " << elements[iTrack] << " rejected (dpt = " << -trk.first
2335  << " GeV/c, algo = " << trackref->algo() << ")";
2336  // Just rescale the nth worst track momentum to equalize the calo energy
2337  } else {
2338  if (isPrimary)
2339  break;
2340  corrTrack = iTrack;
2341  corrFact = (caloEnergy - wouldBeTotalChargedMomentum) / elements[trk.second.first].trackRef()->p();
2342  if (trackref->p() * corrFact < 0.05) {
2343  corrFact = 0.;
2344  active[iTrack] = false;
2345  }
2346  totalChargedMomentum -= trackref->p() * (1. - corrFact);
2347  LogTrace("PFAlgo|createCandidatesHCAL")
2348  << "\tElement " << elements[iTrack] << " (dpt = " << -trk.first << " GeV/c, algo = " << trackref->algo()
2349  << ") rescaled by " << corrFact << " Now the total charged momentum is " << totalChargedMomentum;
2350  break;
2351  }
2352  }
2353  }
2354 
2355  // New determination of the calo and track resolution avec track deletion/rescaling.
2356  caloResolution = neutralHadronEnergyResolution(totalChargedMomentum, hclusterref->positionREP().Eta());
2357  caloResolution *= totalChargedMomentum;
2358  caloResolution = std::sqrt(caloResolution * caloResolution + muonHCALError + muonECALError);
2359 
2360  // Check if the charged momentum is still very inconsistent with the calo measurement.
2361  // In this case, just drop all tracks from 4th and 5th iteration linked to this block
2362 
2363  if (rejectTracks_Step45_ && sortedTracks.size() > 1 &&
2364  totalChargedMomentum - caloEnergy > nSigmaTRACK_ * caloResolution) {
2365  for (auto const& trk : associatedTracks) {
2366  unsigned iTrack = trk.second.first;
2367  reco::TrackRef trackref = elements[iTrack].trackRef();
2368  if (!active[iTrack])
2369  continue;
2370 
2371  double dptRel = fabs(trk.first) / trackref->pt() * 100;
2372  bool isPrimaryOrSecondary = isFromSecInt(elements[iTrack], "all");
2373 
2374  if (isPrimaryOrSecondary && dptRel < dptRel_DispVtx_)
2375  continue;
2376 
2377  if (PFTrackAlgoTools::step5(trackref->algo())) {
2378  active[iTrack] = false;
2379  totalChargedMomentum -= trackref->p();
2380 
2381  LogTrace("PFAlgo|createCandidatesHCAL")
2382  << "\tElement " << elements[iTrack] << " rejected (dpt = " << -trk.first
2383  << " GeV/c, algo = " << trackref->algo() << ")";
2384  }
2385  }
2386  }
2387 
2388  // New determination of the calo and track resolution avec track deletion/rescaling.
2389  caloResolution = neutralHadronEnergyResolution(totalChargedMomentum, hclusterref->positionREP().Eta());
2390  caloResolution *= totalChargedMomentum;
2391  caloResolution = std::sqrt(caloResolution * caloResolution + muonHCALError + muonECALError);
2392 
2393  // Make PF candidates with the remaining tracks in the block
2394  sumpError2 = 0.;
2395  for (auto const& trk : associatedTracks) {
2396  unsigned iTrack = trk.second.first;
2397  if (!active[iTrack])
2398  continue;
2399  reco::TrackRef trackRef = elements[iTrack].trackRef();
2400  double trackMomentum = trackRef->p();
2401  double dp = trackRef->qoverpError() * trackMomentum * trackMomentum;
2402  unsigned tmpi = reconstructTrack(elements[iTrack]);
2403 
2404  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iTrack);
2405  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHcal);
2406  setHcalDepthInfo((*pfCandidates_)[tmpi], *hclusterref);
2407  auto myEcals = associatedEcals.equal_range(iTrack);
2408  for (auto ii = myEcals.first; ii != myEcals.second; ++ii) {
2409  unsigned iEcal = ii->second.second;
2410  if (active[iEcal])
2411  continue;
2412  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iEcal);
2413  }
2414 
2415  if (useHO_) {
2416  auto myHOs = associatedHOs.equal_range(iTrack);
2417  for (auto ii = myHOs.first; ii != myHOs.second; ++ii) {
2418  unsigned iHO = ii->second.second;
2419  if (active[iHO])
2420  continue;
2421  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHO);
2422  }
2423  }
2424 
2425  if (iTrack == corrTrack) {
2426  if (corrFact < 0.)
2427  corrFact = 0.; // protect against negative scaling
2428  (*pfCandidates_)[tmpi].rescaleMomentum(corrFact);
2429  trackMomentum *= corrFact;
2430  }
2431  chargedHadronsIndices.push_back(tmpi);
2432  chargedHadronsInBlock.push_back(iTrack);
2433  active[iTrack] = false;
2434  hcalP.push_back(trackMomentum);
2435  hcalDP.push_back(dp);
2436  if (dp / trackMomentum > maxDPovP)
2437  maxDPovP = dp / trackMomentum;
2438  sumpError2 += dp * dp;
2439  }
2440 
2441  // The total uncertainty of the difference Calo-Track
2442  double totalError = sqrt(sumpError2 + caloResolution * caloResolution);
2443 
2444 #ifdef EDM_ML_DEBUG
2445  LogTrace("PFAlgo|createCandidatesHCAL")
2446  << "\tCompare Calo Energy to total charged momentum " << endl
2447  << "\t\tsum p = " << totalChargedMomentum << " +- " << sqrt(sumpError2) << endl
2448  << "\t\tsum ecal = " << totalEcal << endl
2449  << "\t\tsum hcal = " << totalHcal << endl
2450  << "\t\t => Calo Energy = " << caloEnergy << " +- " << caloResolution << endl
2451  << "\t\t => Calo Energy- total charged momentum = " << caloEnergy - totalChargedMomentum << " +- "
2452  << totalError;
2453 #endif
2454 
2455  /* */
2456 
2458  double nsigma = nSigmaHCAL(totalChargedMomentum, hclusterref->positionREP().Eta());
2459  //double nsigma = nSigmaHCAL(caloEnergy,hclusterref->positionREP().Eta());
2460  if (abs(totalChargedMomentum - caloEnergy) < nsigma * totalError) {
2461  // deposited caloEnergy compatible with total charged momentum
2462  // if tracking errors are large take weighted average
2463 
2464 #ifdef EDM_ML_DEBUG
2465  LogTrace("PFAlgo|createCandidatesHCAL")
2466  << "\t\tcase 1: COMPATIBLE "
2467  << "|Calo Energy- total charged momentum| = " << abs(caloEnergy - totalChargedMomentum) << " < " << nsigma
2468  << " x " << totalError;
2469  if (maxDPovP < 0.1)
2470  LogTrace("PFAlgo|createCandidatesHCAL") << "\t\t\tmax DP/P = " << maxDPovP << " less than 0.1: do nothing ";
2471  else
2472  LogTrace("PFAlgo|createCandidatesHCAL")
2473  << "\t\t\tmax DP/P = " << maxDPovP << " > 0.1: take weighted averages ";
2474 #endif
2475 
2476  // if max DP/P < 10% do nothing
2477  if (maxDPovP > 0.1) {
2478  // for each track associated to hcal
2479  // int nrows = tkIs.size();
2480  int nrows = chargedHadronsIndices.size();
2481  TMatrixTSym<double> a(nrows);
2482  TVectorD b(nrows);
2483  TVectorD check(nrows);
2484  double sigma2E = caloResolution * caloResolution;
2485  for (int i = 0; i < nrows; i++) {
2486  double sigma2i = hcalDP[i] * hcalDP[i];
2487  LogTrace("PFAlgo|createCandidatesHCAL")
2488  << "\t\t\ttrack associated to hcal " << i << " P = " << hcalP[i] << " +- " << hcalDP[i];
2489  a(i, i) = 1. / sigma2i + 1. / sigma2E;
2490  b(i) = hcalP[i] / sigma2i + caloEnergy / sigma2E;
2491  for (int j = 0; j < nrows; j++) {
2492  if (i == j)
2493  continue;
2494  a(i, j) = 1. / sigma2E;
2495  } // end loop on j
2496  } // end loop on i
2497 
2498  // solve ax = b
2499  TDecompChol decomp(a);
2500  bool ok = false;
2501  TVectorD x = decomp.Solve(b, ok);
2502  // for each track create a PFCandidate track
2503  // with a momentum rescaled to weighted average
2504  if (ok) {
2505  for (int i = 0; i < nrows; i++) {
2506  // unsigned iTrack = trackInfos[i].index;
2507  unsigned ich = chargedHadronsIndices[i];
2508  double rescaleFactor = x(i) / hcalP[i];
2509  if (rescaleFactor < 0.)
2510  rescaleFactor = 0.; // protect against negative scaling
2511  (*pfCandidates_)[ich].rescaleMomentum(rescaleFactor);
2512 
2513  LogTrace("PFAlgo|createCandidatesHCAL")
2514  << "\t\t\told p " << hcalP[i] << " new p " << x(i) << " rescale " << rescaleFactor;
2515  }
2516  } else {
2517  edm::LogError("PFAlgo::createCandidatesHCAL") << "TDecompChol.Solve returned ok=false";
2518  assert(0);
2519  }
2520  }
2521  }
2522 
2524  else if (caloEnergy > totalChargedMomentum) {
2525  //case 2: caloEnergy > totalChargedMomentum + nsigma*totalError
2526  //there is an excess of energy in the calos
2527  //create a neutral hadron or a photon
2528 
2529  double eNeutralHadron = caloEnergy - totalChargedMomentum;
2530  double ePhoton = (caloEnergy - totalChargedMomentum) /
2531  slopeEcal; // KH: this slopeEcal is computed based on ECAL energy under the hadron hypothesis,
2532  // thought we are creating photons.
2533  // This is a fuzzy case, but it should be better than corrected twice under both egamma and hadron hypotheses.
2534 
2535 #ifdef EDM_ML_DEBUG
2536  if (!sortedTracks.empty()) {
2537  LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tcase 2: NEUTRAL DETECTION " << caloEnergy << " > " << nsigma
2538  << "x" << totalError << " + " << totalChargedMomentum;
2539  LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tneutral activity detected: " << endl
2540  << "\t\t\t photon = " << ePhoton << endl
2541  << "\t\t\tor neutral hadron = " << eNeutralHadron;
2542 
2543  LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tphoton or hadron ?";
2544  }
2545 
2546  if (sortedTracks.empty())
2547  LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tno track -> hadron ";
2548  else
2549  LogTrace("PFAlgo|createCandidatesHCAL")
2550  << "\t\t" << sortedTracks.size() << " tracks -> check if the excess is photonic or hadronic";
2551 #endif
2552 
2553  double ratioMax = 0.;
2554  reco::PFClusterRef maxEcalRef;
2555  unsigned maxiEcal = 9999;
2556 
2557  // for each track associated to hcal: iterator IE ie :
2558 
2559  LogTrace("PFAlgo|createCandidatesHCAL") << "loop over sortedTracks.size()=" << sortedTracks.size();
2560  for (auto const& trk : sortedTracks) {
2561  unsigned iTrack = trk.second;
2562 
2563  PFBlockElement::Type type = elements[iTrack].type();
2565 
2566  reco::TrackRef trackRef = elements[iTrack].trackRef();
2567  assert(!trackRef.isNull());
2568 
2569  auto iae = associatedEcals.find(iTrack);
2570 
2571  if (iae == associatedEcals.end())
2572  continue;
2573 
2574  // double distECAL = iae->second.first;
2575  unsigned iEcal = iae->second.second;
2576 
2577  PFBlockElement::Type typeEcal = elements[iEcal].type();
2578  assert(typeEcal == reco::PFBlockElement::ECAL);
2579 
2580  reco::PFClusterRef clusterRef = elements[iEcal].clusterRef();
2581  assert(!clusterRef.isNull());
2582 
2583  double pTrack = trackRef->p();
2584  double eECAL = clusterRef->energy();
2585  double eECALOverpTrack = eECAL / pTrack;
2586 
2587  if (eECALOverpTrack > ratioMax) {
2588  ratioMax = eECALOverpTrack;
2589  maxEcalRef = clusterRef;
2590  maxiEcal = iEcal;
2591  }
2592 
2593  } // end loop on tracks associated to hcal: iterator IE ie
2594 
2595  std::vector<reco::PFClusterRef> pivotalClusterRef;
2596  std::vector<unsigned> iPivotal;
2597  std::vector<double> particleEnergy, ecalEnergy, hcalEnergy, rawecalEnergy, rawhcalEnergy;
2598  std::vector<::math::XYZVector> particleDirection;
2599 
2600  // If the excess is smaller than the ecal energy, assign the whole
2601  // excess to photons
2602  if (ePhoton < totalEcal || eNeutralHadron - calibEcal < 1E-10) {
2603  if (!maxEcalRef.isNull()) {
2604  // So the merged photon energy is,
2605  mergedPhotonEnergy = ePhoton;
2606  }
2607  } else {
2608  // Otherwise assign the whole ECAL energy to the photons
2609  if (!maxEcalRef.isNull()) {
2610  // So the merged photon energy is,
2611  mergedPhotonEnergy = totalEcalEGMCalib; // KH: use calibrated ECAL energy under the egamma hypothesis
2612  }
2613  // ... and assign the remaining excess to neutral hadrons using the direction of ecal clusters
2614  mergedNeutralHadronEnergy = eNeutralHadron - calibEcal;
2615  }
2616 
2617  if (mergedPhotonEnergy > 0) {
2618  // Split merged photon into photons for each energetic ecal cluster (necessary for jet substructure reconstruction)
2619  // make only one merged photon if less than 2 ecal clusters
2620  // KH: this part still needs review, after using non-corrected ECAL energy for PF hadron calibrations
2621  if (ecalClusters.size() <= 1) {
2622  ecalClusters.clear();
2623  ecalClusters.emplace_back(
2624  maxiEcal,
2625  photonAtECAL,
2626  1.); // KH: calibration factor of 1, which should be ok as long as sumEcalClusters is consistent with photonAtECAL in this case
2627  sumEcalClusters = sqrt(photonAtECAL.Mag2());
2628  }
2629  for (auto const& pae : ecalClusters) {
2630  const double clusterEnergyCalibrated =
2631  sqrt(std::get<1>(pae).Mag2()) *
2632  std::get<2>(
2633  pae); // KH: calibrated under the egamma hypothesis. Note: sumEcalClusters is normally calibrated under egamma hypothesis
2634  particleEnergy.push_back(mergedPhotonEnergy * clusterEnergyCalibrated / sumEcalClusters);
2635  particleDirection.push_back(std::get<1>(pae));
2636  ecalEnergy.push_back(mergedPhotonEnergy * clusterEnergyCalibrated / sumEcalClusters);
2637  hcalEnergy.push_back(0.);
2638  rawecalEnergy.push_back(totalEcal);
2639  rawhcalEnergy.push_back(0.);
2640  pivotalClusterRef.push_back(elements[std::get<0>(pae)].clusterRef());
2641  iPivotal.push_back(std::get<0>(pae));
2642  }
2643  } // mergedPhotonEnergy > 0
2644 
2645  if (mergedNeutralHadronEnergy > 1.0) {
2646  // Split merged neutral hadrons according to directions of energetic ecal clusters (necessary for jet substructure reconstruction)
2647  // make only one merged neutral hadron if less than 2 ecal clusters
2648  if (ecalClusters.size() <= 1) {
2649  ecalClusters.clear();
2650  ecalClusters.emplace_back(
2651  iHcal,
2652  hadronAtECAL,
2653  1.); // KH: calibration factor of 1, which should be ok as long as sumEcalClusters is consistent with photonAtECAL
2654  sumEcalClusters = sqrt(hadronAtECAL.Mag2());
2655  }
2656  for (auto const& pae : ecalClusters) {
2657  const double clusterEnergyCalibrated =
2658  sqrt(std::get<1>(pae).Mag2()) *
2659  std::get<2>(
2660  pae); // KH: calibrated under the egamma hypothesis. Note: sumEcalClusters is normally calibrated under egamma hypothesis
2661  particleEnergy.push_back(mergedNeutralHadronEnergy * clusterEnergyCalibrated / sumEcalClusters);
2662  particleDirection.push_back(std::get<1>(pae));
2663  ecalEnergy.push_back(0.);
2664  hcalEnergy.push_back(mergedNeutralHadronEnergy * clusterEnergyCalibrated / sumEcalClusters);
2665  rawecalEnergy.push_back(0.);
2666  rawhcalEnergy.push_back(totalHcal);
2667  pivotalClusterRef.push_back(hclusterref);
2668  iPivotal.push_back(iHcal);
2669  }
2670  } //mergedNeutralHadronEnergy > 1.0
2671 
2672  // reconstructing a merged neutral
2673  // the type of PFCandidate is known from the
2674  // reference to the pivotal cluster.
2675 
2676  for (unsigned iPivot = 0; iPivot < iPivotal.size(); ++iPivot) {
2677  if (particleEnergy[iPivot] < 0.)
2678  edm::LogWarning("PFAlgo|createCandidatesHCAL")
2679  << "ALARM = Negative energy for iPivot=" << iPivot << ", " << particleEnergy[iPivot];
2680 
2681  const bool useDirection = true;
2682  auto& neutral = (*pfCandidates_)[reconstructCluster(*pivotalClusterRef[iPivot],
2683  particleEnergy[iPivot],
2684  useDirection,
2685  particleDirection[iPivot].X(),
2686  particleDirection[iPivot].Y(),
2687  particleDirection[iPivot].Z())];
2688 
2689  neutral.setEcalEnergy(rawecalEnergy[iPivot], ecalEnergy[iPivot]);
2690  if (!useHO_) {
2691  neutral.setHcalEnergy(rawhcalEnergy[iPivot], hcalEnergy[iPivot]);
2692  neutral.setHoEnergy(0., 0.);
2693  } else { // useHO_
2694  if (rawhcalEnergy[iPivot] == 0.) { // photons should be here
2695  neutral.setHcalEnergy(0., 0.);
2696  neutral.setHoEnergy(0., 0.);
2697  } else {
2698  neutral.setHcalEnergy(max(rawhcalEnergy[iPivot] - totalHO, 0.0),
2699  hcalEnergy[iPivot] * max(1. - totalHO / rawhcalEnergy[iPivot], 0.));
2700  neutral.setHoEnergy(totalHO, totalHO * hcalEnergy[iPivot] / rawhcalEnergy[iPivot]);
2701  }
2702  }
2703  neutral.setPs1Energy(0.);
2704  neutral.setPs2Energy(0.);
2705  neutral.set_mva_nothing_gamma(-1.);
2706  // neutral.addElement(&elements[iPivotal]);
2707  // neutral.addElementInBlock(blockref, iPivotal[iPivot]);
2708  neutral.addElementInBlock(blockref, iHcal);
2709  for (unsigned iTrack : chargedHadronsInBlock) {
2710  neutral.addElementInBlock(blockref, iTrack);
2711  // Assign the position of the track at the ECAL entrance
2712  const ::math::XYZPointF& chargedPosition =
2713  dynamic_cast<const reco::PFBlockElementTrack*>(&elements[iTrack])->positionAtECALEntrance();
2714  neutral.setPositionAtECALEntrance(chargedPosition);
2715 
2716  auto myEcals = associatedEcals.equal_range(iTrack);
2717  for (auto ii = myEcals.first; ii != myEcals.second; ++ii) {
2718  unsigned iEcal = ii->second.second;
2719  if (active[iEcal])
2720  continue;
2721  neutral.addElementInBlock(blockref, iEcal);
2722  }
2723  }
2724  }
2725 
2726  } // excess of energy
2727 
2728  // will now share the hcal energy between the various charged hadron
2729  // candidates, taking into account the potential neutral hadrons
2730 
2731  //JB: The question is: we've resolved the merged photons cleanly, but how should
2732  //the remaining hadrons be assigned the remaining ecal energy?
2733  //*Temporary solution*: follow HCAL example with fractions...
2734 
2735  // remove the energy of the potential neutral hadron
2736  double totalHcalEnergyCalibrated = std::max(calibHcal - mergedNeutralHadronEnergy, 0.);
2737  // similarly for the merged photons
2738  double totalEcalEnergyCalibrated = std::max(calibEcal - mergedPhotonEnergy, 0.);
2739  // share between the charged hadrons
2740 
2741  //COLIN can compute this before
2742  // not exactly equal to sum p, this is sum E
2743  double chargedHadronsTotalEnergy = 0;
2744  for (unsigned index : chargedHadronsIndices) {
2745  reco::PFCandidate& chargedHadron = (*pfCandidates_)[index];
2746  chargedHadronsTotalEnergy += chargedHadron.energy();
2747  }
2748 
2749  for (unsigned index : chargedHadronsIndices) {
2750  reco::PFCandidate& chargedHadron = (*pfCandidates_)[index];
2751  float fraction = chargedHadron.energy() / chargedHadronsTotalEnergy;
2752 
2753  if (!useHO_) {
2754  chargedHadron.setHcalEnergy(fraction * totalHcal, fraction * totalHcalEnergyCalibrated);
2755  chargedHadron.setHoEnergy(0., 0.);
2756  } else {
2757  chargedHadron.setHcalEnergy(fraction * max(totalHcal - totalHO, 0.0),
2758  fraction * totalHcalEnergyCalibrated * (1. - totalHO / totalHcal));
2759  chargedHadron.setHoEnergy(fraction * totalHO, fraction * totalHO * totalHcalEnergyCalibrated / totalHcal);
2760  }
2761  //JB: fixing up (previously omitted) setting of ECAL energy gouzevit
2762  chargedHadron.setEcalEnergy(fraction * totalEcal, fraction * totalEcalEnergyCalibrated);
2763  }
2764 
2765  // Finally treat unused ecal satellites as photons.
2766  for (auto const& ecalSatellite : ecalSatellites) {
2767  // Ignore satellites already taken
2768  unsigned iEcal = std::get<0>(ecalSatellite.second);
2769  if (!active[iEcal])
2770  continue;
2771 
2772  // Sanity checks again (well not useful, this time!)
2773  PFBlockElement::Type type = elements[iEcal].type();
2775  PFClusterRef eclusterref = elements[iEcal].clusterRef();
2776  assert(!eclusterref.isNull());
2777 
2778  // Lock the cluster
2779  active[iEcal] = false;
2780 
2781  // Find the associated tracks
2782  std::multimap<double, unsigned> assTracks;
2783  block.associatedElements(iEcal, linkData, assTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
2784 
2785  // Create a photon
2786  double ecalClusterEnergyCalibrated =
2787  sqrt(std::get<1>(ecalSatellite.second).Mag2()) *
2788  std::get<2>(
2789  ecalSatellite.second); // KH: calibrated under the egamma hypothesis (rawEcalClusterEnergy * calibration)
2790  auto& cand = (*pfCandidates_)[reconstructCluster(*eclusterref, ecalClusterEnergyCalibrated)];
2791  cand.setEcalEnergy(eclusterref->energy(), ecalClusterEnergyCalibrated);
2792  cand.setHcalEnergy(0., 0.);
2793  cand.setHoEnergy(0., 0.);
2794  cand.setPs1Energy(associatedPSs[iEcal].first);
2795  cand.setPs2Energy(associatedPSs[iEcal].second);
2796  cand.addElementInBlock(blockref, iEcal);
2797  cand.addElementInBlock(blockref, sortedTracks.begin()->second);
2798 
2799  if (fabs(eclusterref->energy() - sqrt(std::get<1>(ecalSatellite.second).Mag2())) > 1e-3 ||
2800  fabs(eclusterref->correctedEnergy() - ecalClusterEnergyCalibrated) > 1e-3)
2801  edm::LogWarning("PFAlgo:processBlock")
2802  << "ecalCluster vs ecalSatellites look inconsistent (eCluster E, calibE, ecalSatellite E, calib E): "
2803  << eclusterref->energy() << " " << eclusterref->correctedEnergy() << " "
2804  << sqrt(std::get<1>(ecalSatellite.second).Mag2()) << " " << ecalClusterEnergyCalibrated;
2805 
2806  } // ecalSatellites
2807 
2808  } // hcalIs
2809  // end loop on hcal element iHcal= hcalIs[i]
2810  LogTrace("PFAlgo|createCandidatesHCAL") << "end of function PFAlgo::createCandidatesHCAL";
2811 }
static bool isIsolatedMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:94
double ptError_
Definition: PFAlgo.h:297
std::vector< double > muonHCAL_
Variables for muons and fakes.
Definition: PFAlgo.h:293
bool isPrimary(const SimTrack &simTrk, const PSimHit *simHit)
static bool isMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:48
double errorScale(const reco::TrackBase::TrackAlgorithm &, const std::vector< double > &)
std::unique_ptr< reco::PFCandidateCollection > pfCandidates_
Definition: PFAlgo.h:226
unsigned reconstructCluster(const reco::PFCluster &cluster, double particleEnergy, bool useDirection=false, double particleX=0., double particleY=0., double particleZ=0.)
Definition: PFAlgo.cc:3220
#define X(str)
Definition: MuonsGrabber.cc:38
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float > > XYZPointF
point in space with cartesian internal representation
Definition: Point3D.h:10
Log< level::Error, false > LogError
bool isFromSecInt(const reco::PFBlockElement &eTrack, std::string order) const
Definition: PFAlgo.cc:3470
assert(be >=bs)
bool rejectTracks_Step45_
Definition: PFAlgo.h:277
#define LogTrace(id)
std::vector< double > factors45_
Definition: PFAlgo.h:298
bool useHO_
Definition: PFAlgo.h:260
U second(std::pair< T, U > const &p)
std::vector< unsigned > hcalIs
Definition: PFAlgo.h:42
T sqrt(T t)
Definition: SSEVec.h:19
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
std::vector< double > muonECAL_
Definition: PFAlgo.h:294
void setPositionAtECALEntrance(float x, float y, float z)
set position at ECAL entrance
PFEnergyCalibration & calibration_
Definition: PFAlgo.h:257
void energyEmHad(double t, double &e, double &h, double eta, double phi) const
bool isNull() const
Checks for null.
Definition: Ref.h:235
ii
Definition: cuy.py:589
unsigned reconstructTrack(const reco::PFBlockElement &elt, bool allowLoose=false)
Definition: PFAlgo.cc:3141
std::vector< double > muonHO_
Definition: PFAlgo.h:295
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
double dptRel_DispVtx_
Definition: PFAlgo.h:285
static bool isLooseMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:57
double neutralHadronEnergyResolution(double clusterEnergy, double clusterEta) const
todo: use PFClusterTools for this
Definition: PFAlgo.cc:3361
double b
Definition: hdecay.h:118
void setHcalDepthInfo(reco::PFCandidate &cand, const reco::PFCluster &cluster) const
Definition: PFAlgo.cc:3331
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:41
double a
Definition: hdecay.h:119
double nSigmaHCAL(double clusterEnergy, double clusterEta) const
Definition: PFAlgo.cc:3371
bool rejectTracks_Bad_
Definition: PFAlgo.h:276
bool step5(const reco::TrackBase::TrackAlgorithm &)
Log< level::Warning, false > LogWarning
double nSigmaTRACK_
Definition: PFAlgo.h:296
math::XYZVector XYZVector
Definition: RawParticle.h:26

◆ createCandidatesHCALUnlinked()

void PFAlgo::createCandidatesHCALUnlinked ( const reco::PFBlock block,
reco::PFBlock::LinkData linkData,
const edm::OwnVector< reco::PFBlockElement > &  elements,
std::vector< bool > &  active,
const reco::PFBlockRef blockref,
ElementIndices inds,
std::vector< bool > &  deadArea 
)
private

Definition at line 2813 of file PFAlgo.cc.

References cms::cuda::assert(), groupFilesInBlocks::block, calibration_, ECAL, reco::PFBlockElement::ECAL, bookConverter::elements, PFEnergyCalibration::energyEmHad(), reco::PFBlockElement::HCAL, hltEgammaHLTExtra_cfi::hcal, ElementIndices::hcalIs, PFLayer::HF_EM, PFLayer::HF_HAD, photonIsolationHIProducer_cfi::ho, DigiToRawDM_cff::HO, reco::PFBlockElement::HO, edm::Ref< C, T, F >::isNull(), reco::PFBlock::LINKTEST_ALL, LogTrace, SiStripPI::max, reconstructCluster(), and useHO_.

Referenced by processBlock().

2819  {
2820  // Processing the remaining HCAL clusters
2821  LogTrace("PFAlgo|createCandidatesHCALUnlinked")
2822  << "start of function PFAlgo::createCandidatesHCALUnlinked, hcalIs.size()=" << inds.hcalIs.size();
2823 
2824  // --------------- loop remaining hcal ------------------
2825 
2826  for (unsigned iHcal : inds.hcalIs) {
2827  // Keep ECAL and HO elements for reference in the PFCandidate
2828  std::vector<unsigned> ecalRefs;
2829  std::vector<unsigned> hoRefs;
2830 
2831  LogTrace("PFAlgo|createCandidatesHCALUnlinked") << elements[iHcal] << " ";
2832 
2833  if (!active[iHcal]) {
2834  LogTrace("PFAlgo|createCandidatesHCALUnlinked") << "not active " << iHcal;
2835  continue;
2836  }
2837 
2838  // Find the ECAL elements linked to it
2839  std::multimap<double, unsigned> ecalElems;
2840  block.associatedElements(iHcal, linkData, ecalElems, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
2841 
2842  // Loop on these ECAL elements
2843  float totalEcal = 0.;
2844  float ecalMax = 0.;
2845  reco::PFClusterRef eClusterRef;
2846  for (auto const& ecal : ecalElems) {
2847  unsigned iEcal = ecal.second;
2848  double dist = ecal.first;
2849  PFBlockElement::Type type = elements[iEcal].type();
2851 
2852  // Check if already used
2853  if (!active[iEcal])
2854  continue;
2855 
2856  // Check the distance (one HCALPlusECAL tower, roughly)
2857  // if ( dist > 0.15 ) continue;
2858 
2859  //COLINFEB16
2860  // what could be done is to
2861  // - link by rechit.
2862  // - take in the neutral hadron all the ECAL clusters
2863  // which are within the same CaloTower, according to the distance,
2864  // except the ones which are closer to another HCAL cluster.
2865  // - all the other ECAL linked to this HCAL are photons.
2866  //
2867  // about the closest HCAL cluster.
2868  // it could maybe be easier to loop on the ECAL clusters first
2869  // to cut the links to all HCAL clusters except the closest, as is
2870  // done in the first track loop. But maybe not!
2871  // or add an helper function to the PFAlgo class to ask
2872  // if a given element is the closest of a given type to another one?
2873 
2874  // Check if not closer from another free HCAL
2875  std::multimap<double, unsigned> hcalElems;
2876  block.associatedElements(iEcal, linkData, hcalElems, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
2877 
2878  const bool isClosest = std::none_of(hcalElems.begin(), hcalElems.end(), [&](auto const& hcal) {
2879  return active[hcal.second] && hcal.first < dist;
2880  });
2881 
2882  if (!isClosest)
2883  continue;
2884 
2885 #ifdef EDM_ML_DEBUG
2886  LogTrace("PFAlgo|createCandidatesHCALUnlinked")
2887  << "\telement " << elements[iEcal] << " linked with dist " << dist;
2888  LogTrace("PFAlgo|createCandidatesHCALUnlinked") << "Added to HCAL cluster to form a neutral hadron";
2889 #endif
2890 
2891  reco::PFClusterRef eclusterRef = elements[iEcal].clusterRef();
2892  assert(!eclusterRef.isNull());
2893 
2894  // KH: use raw ECAL energy for PF hadron calibration_.
2895  double ecalEnergy = eclusterRef->energy(); // ecalEnergy = eclusterRef->correctedEnergy();
2896 
2897  totalEcal += ecalEnergy;
2898  if (ecalEnergy > ecalMax) {
2899  ecalMax = ecalEnergy;
2900  eClusterRef = eclusterRef;
2901  }
2902 
2903  ecalRefs.push_back(iEcal);
2904  active[iEcal] = false;
2905 
2906  } // End loop ECAL
2907 
2908  // Now find the HO clusters linked to the HCAL cluster
2909  double totalHO = 0.;
2910  double hoMax = 0.;
2911  //unsigned jHO = 0;
2912  if (useHO_) {
2913  std::multimap<double, unsigned> hoElems;
2914  block.associatedElements(iHcal, linkData, hoElems, reco::PFBlockElement::HO, reco::PFBlock::LINKTEST_ALL);
2915 
2916  // Loop on these HO elements
2917  // double totalHO = 0.;
2918  // double hoMax = 0.;
2919  // unsigned jHO = 0;
2920  reco::PFClusterRef hoClusterRef;
2921  for (auto const& ho : hoElems) {
2922  unsigned iHO = ho.second;
2923  double dist = ho.first;
2924  PFBlockElement::Type type = elements[iHO].type();
2926 
2927  // Check if already used
2928  if (!active[iHO])
2929  continue;
2930 
2931  // Check the distance (one HCALPlusHO tower, roughly)
2932  // if ( dist > 0.15 ) continue;
2933 
2934  // Check if not closer from another free HCAL
2935  std::multimap<double, unsigned> hcalElems;
2936  block.associatedElements(iHO, linkData, hcalElems, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
2937 
2938  const bool isClosest = std::none_of(hcalElems.begin(), hcalElems.end(), [&](auto const& hcal) {
2939  return active[hcal.second] && hcal.first < dist;
2940  });
2941 
2942  if (!isClosest)
2943  continue;
2944 
2945 #ifdef EDM_ML_DEBUG
2946  if (useHO_) {
2947  LogTrace("PFAlgo|createCandidatesHCALUnlinked")
2948  << "\telement " << elements[iHO] << " linked with dist " << dist;
2949  LogTrace("PFAlgo|createCandidatesHCALUnlinked") << "Added to HCAL cluster to form a neutral hadron";
2950  }
2951 #endif
2952 
2953  reco::PFClusterRef hoclusterRef = elements[iHO].clusterRef();
2954  assert(!hoclusterRef.isNull());
2955 
2956  double hoEnergy =
2957  hoclusterRef->energy(); // calibration_.energyEm(*hoclusterRef,ps1Ene,ps2Ene,crackCorrection);
2958 
2959  totalHO += hoEnergy;
2960  if (hoEnergy > hoMax) {
2961  hoMax = hoEnergy;
2962  hoClusterRef = hoclusterRef;
2963  //jHO = iHO;
2964  }
2965 
2966  hoRefs.push_back(iHO);
2967  active[iHO] = false;
2968 
2969  } // End loop HO
2970  }
2971 
2972  PFClusterRef hclusterRef = elements[iHcal].clusterRef();
2973  assert(!hclusterRef.isNull());
2974 
2975  // HCAL energy
2976  double totalHcal = hclusterRef->energy();
2977  // Include the HO energy
2978  if (useHO_)
2979  totalHcal += totalHO;
2980 
2981  // Calibration
2982  double calibEcal = totalEcal > 0. ? totalEcal : 0.;
2983  double calibHcal = std::max(0., totalHcal);
2984  if (hclusterRef->layer() == PFLayer::HF_HAD || hclusterRef->layer() == PFLayer::HF_EM) {
2985  calibEcal = totalEcal;
2986  } else {
2988  -1., calibEcal, calibHcal, hclusterRef->positionREP().Eta(), hclusterRef->positionREP().Phi());
2989  }
2990 
2991  auto& cand = (*pfCandidates_)[reconstructCluster(*hclusterRef, calibEcal + calibHcal)];
2992 
2993  cand.setEcalEnergy(totalEcal, calibEcal);
2994  if (!useHO_) {
2995  cand.setHcalEnergy(totalHcal, calibHcal);
2996  cand.setHoEnergy(0., 0.);
2997  } else {
2998  cand.setHcalEnergy(max(totalHcal - totalHO, 0.0), calibHcal * (1. - totalHO / totalHcal));
2999  cand.setHoEnergy(totalHO, totalHO * calibHcal / totalHcal);
3000  }
3001  cand.setPs1Energy(0.);
3002  cand.setPs2Energy(0.);
3003  cand.addElementInBlock(blockref, iHcal);
3004  for (auto const& ec : ecalRefs)
3005  cand.addElementInBlock(blockref, ec);
3006  for (auto const& ho : hoRefs)
3007  cand.addElementInBlock(blockref, ho);
3008 
3009  } //loop hcal elements
3010 }
unsigned reconstructCluster(const reco::PFCluster &cluster, double particleEnergy, bool useDirection=false, double particleX=0., double particleY=0., double particleZ=0.)
Definition: PFAlgo.cc:3220
assert(be >=bs)
#define LogTrace(id)
bool useHO_
Definition: PFAlgo.h:260
std::vector< unsigned > hcalIs
Definition: PFAlgo.h:42
PFEnergyCalibration & calibration_
Definition: PFAlgo.h:257
void energyEmHad(double t, double &e, double &h, double eta, double phi) const
bool isNull() const
Checks for null.
Definition: Ref.h:235

◆ createCandidatesHF()

void PFAlgo::createCandidatesHF ( const reco::PFBlock block,
reco::PFBlock::LinkData linkData,
const edm::OwnVector< reco::PFBlockElement > &  elements,
std::vector< bool > &  active,
const reco::PFBlockRef blockref,
ElementIndices inds 
)
private

Definition at line 1259 of file PFAlgo.cc.

References cms::cuda::assert(), groupFilesInBlocks::block, fftjetpileupestimator_calo_uncalib_cfi::c0, alignmentValidation::c1, Calorimetry_cff::dp, bookConverter::elements, PFEnergyCalibrationHF::energyEm(), PFEnergyCalibrationHF::energyEmHad(), PFEnergyCalibrationHF::energyHad(), HcalResponse_cfi::energyHF, DivergingColor::frac, PFEnergyCalibrationHF::getcalibHF_use(), PFLayer::HF_EM, PFLayer::HF_HAD, reco::PFBlockElement::HFEM, ElementIndices::hfEmIs, hfEnergyResolution(), ElementIndices::hfHadIs, cuy::ii, edm::Ref< C, T, F >::isNull(), reco::PFBlock::LINKTEST_ALL, LogTrace, SiStripPI::max, nSigmaHFEM(), nSigmaHFHAD(), reconstructCluster(), reconstructTrack(), mathSSE::sqrt(), thepfEnergyCalibrationHF_, reco::PFBlockElement::TRACK, ElementIndices::trackIs, and reco::btau::trackMomentum.

Referenced by processBlock().

1264  {
1265  LogTrace("PFAlgo|createCandidatesHF") << "starting function PFAlgo::createCandidatesHF";
1266 
1267  bool trackInBlock = !inds.trackIs.empty();
1268  // inds.trackIs can be empty, even if there are tracks in this block,
1269  // but what we want to check is if this block has any track including inactive ones
1270  if (!trackInBlock)
1271  for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
1272  PFBlockElement::Type type = elements[iEle].type();
1273  if (type == PFBlockElement::TRACK) {
1274  trackInBlock = true;
1275  break;
1276  }
1277  }
1278  // there is at least one HF element in this block.
1279  // in case of no track, all elements must be HF
1280  if (!trackInBlock)
1281  assert(inds.hfEmIs.size() + inds.hfHadIs.size() == elements.size());
1282 
1283  //
1284  // Dealing with a block with at least one track
1285  // Occasionally, there are only inactive tracks and multiple HF clusters. Consider such blocks too.
1286  //
1287  if (trackInBlock) { // count any tracks (not only active tracks)
1288  // sorted tracks associated with a HfHad cluster
1289  std::multimap<double, unsigned> sortedTracks;
1290  std::multimap<double, unsigned> sortedTracksActive; // only active ones
1291  // HfEms associated with tracks linked to a HfHad cluster
1292  std::multimap<unsigned, std::pair<double, unsigned>> associatedHfEms;
1293  // Temporary map for HfEm satellite clusters
1294  std::multimap<double, std::pair<unsigned, double>> hfemSatellites;
1295 
1296  //
1297  // Loop over active HfHad clusters
1298  //
1299  for (unsigned iHfHad : inds.hfHadIs) {
1300  PFBlockElement::Type type = elements[iHfHad].type();
1301  assert(type == PFBlockElement::HFHAD);
1302 
1303  PFClusterRef hclusterRef = elements[iHfHad].clusterRef();
1304  assert(!hclusterRef.isNull());
1305 
1306  sortedTracks.clear();
1307  sortedTracksActive.clear();
1308  associatedHfEms.clear();
1309  hfemSatellites.clear();
1310 
1311  // Look for associated tracks
1312  block.associatedElements(
1313  iHfHad, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
1314 
1315  LogTrace("PFAlgo|createCandidatesHF") << "elements[" << iHfHad << "]=" << elements[iHfHad];
1316 
1317  if (sortedTracks.empty()) {
1318  LogTrace("PFAlgo|createCandidatesHCF") << "\tno associated tracks, keep for later";
1319  continue;
1320  }
1321 
1322  // Lock the HFHAD cluster
1323  active[iHfHad] = false;
1324 
1325  LogTrace("PFAlgo|createCandidatesHF") << sortedTracks.size() << " associated tracks:";
1326 
1327  double totalChargedMomentum = 0.;
1328  double sumpError2 = 0.;
1329 
1330  //
1331  // Loop over all tracks associated to this HFHAD cluster
1332  //
1333  for (auto const& trk : sortedTracks) {
1334  unsigned iTrack = trk.second;
1335 
1336  // Check the track has not already been used
1337  if (!active[iTrack])
1338  continue;
1339  // Sanity check 1
1340  PFBlockElement::Type type = elements[iTrack].type();
1342  // Sanity check 2
1343  reco::TrackRef trackRef = elements[iTrack].trackRef();
1344  assert(!trackRef.isNull());
1345 
1346  // Introduce tracking errors
1347  double trackMomentum = trackRef->p();
1348  totalChargedMomentum += trackMomentum;
1349 
1350  // Also keep the total track momentum error for comparison with the calo energy
1351  double dp = trackRef->qoverpError() * trackMomentum * trackMomentum;
1352  sumpError2 += dp * dp;
1353 
1354  // Store active tracks for 2nd loop to create charged hadrons
1355  sortedTracksActive.emplace(trk);
1356 
1357  // look for HFEM elements associated to iTrack (associated to iHfHad)
1358  std::multimap<double, unsigned> sortedHfEms;
1359  block.associatedElements(
1360  iTrack, linkData, sortedHfEms, reco::PFBlockElement::HFEM, reco::PFBlock::LINKTEST_ALL);
1361 
1362  LogTrace("PFAlgo|createCandidatesHF") << "number of HfEm elements linked to this track: " << sortedHfEms.size();
1363 
1364  bool connectedToHfEm = false; // Will become true if there is at least one HFEM cluster connected
1365 
1366  //
1367  // Loop over all HFEM clusters connected to iTrack
1368  //
1369  for (auto const& hfem : sortedHfEms) {
1370  unsigned iHfEm = hfem.second;
1371  double dist = hfem.first;
1372 
1373  // Ignore HFEM cluters already used
1374  if (!active[iHfEm]) {
1375  LogTrace("PFAlgo|createCandidatesHF") << "cluster locked";
1376  continue;
1377  }
1378 
1379  // Sanity checks
1380  PFBlockElement::Type type = elements[iHfEm].type();
1381  assert(type == PFBlockElement::HFEM);
1382  PFClusterRef eclusterRef = elements[iHfEm].clusterRef();
1383  assert(!eclusterRef.isNull());
1384 
1385  // Check if this HFEM is not closer to another track - ignore it in that case
1386  std::multimap<double, unsigned> sortedTracksHfEm;
1387  block.associatedElements(
1388  iHfEm, linkData, sortedTracksHfEm, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
1389  unsigned jTrack = sortedTracksHfEm.begin()->second;
1390  if (jTrack != iTrack)
1391  continue;
1392 
1393  double distHfEm = block.dist(jTrack, iHfEm, linkData, reco::PFBlock::LINKTEST_ALL);
1394  double hfemEnergy = eclusterRef->energy();
1395 
1396  if (!connectedToHfEm) { // This is the closest HFEM cluster - will add its energy later
1397 
1398  LogTrace("PFAlgo|createCandidatesHF") << "closest: " << elements[iHfEm];
1399  connectedToHfEm = true;
1400  std::pair<unsigned, double> satellite(iHfEm, hfemEnergy);
1401  hfemSatellites.emplace(-1., satellite);
1402 
1403  } else { // Keep satellite clusters for later
1404 
1405  // KH: same as above.
1406  std::pair<unsigned, double> satellite(iHfEm, hfemEnergy);
1407  hfemSatellites.emplace(dist, satellite);
1408  }
1409 
1410  std::pair<double, unsigned> associatedHfEm(distHfEm, iHfEm);
1411  associatedHfEms.emplace(iTrack, associatedHfEm);
1412 
1413  } // End loop hfem associated to iTrack
1414  } // sortedTracks
1415 
1416  // HfHad energy
1417  double uncalibratedenergyHfHad = hclusterRef->energy();
1418  double energyHfHad = uncalibratedenergyHfHad;
1420  energyHfHad = thepfEnergyCalibrationHF_.energyHad( // HAD only calibration
1421  uncalibratedenergyHfHad,
1422  hclusterRef->positionREP().Eta(),
1423  hclusterRef->positionREP().Phi());
1424  }
1425  double calibFactorHfHad = (uncalibratedenergyHfHad > 0.) ? energyHfHad / uncalibratedenergyHfHad : 1.;
1426 
1427  // HfEm energy
1428  double energyHfEmTmp = 0.;
1429  double uncalibratedenergyHfEmTmp = 0.;
1430  double energyHfEm = 0.;
1431  double uncalibratedenergyHfEm = 0.;
1432 
1433  // estimated HF resolution and track p error
1434  double caloResolution = hfEnergyResolution(totalChargedMomentum);
1435  caloResolution *= totalChargedMomentum;
1436  double totalError = sqrt(caloResolution * caloResolution + sumpError2);
1437  double nsigmaHFEM = nSigmaHFEM(totalChargedMomentum);
1438  double nsigmaHFHAD = nSigmaHFHAD(totalChargedMomentum);
1439 
1440  // Handle case that no active track gets associated to HfHad cluster
1441  if (sortedTracksActive.empty()) {
1442  // look for HFEM elements associated to iHfHad
1443  std::multimap<double, unsigned> sortedHfEms;
1444  std::multimap<double, unsigned> sortedHfEmsActive;
1445  block.associatedElements(
1446  iHfHad, linkData, sortedHfEms, reco::PFBlockElement::HFEM, reco::PFBlock::LINKTEST_ALL);
1447  //
1448  // If iHfHad is connected to HFEM cluster, Loop over all of them
1449  //
1450  if (!sortedHfEms.empty()) {
1451  for (auto const& hfem : sortedHfEms) {
1452  unsigned iHfEm = hfem.second;
1453  // Ignore HFEM cluters already used
1454  if (!active[iHfEm])
1455  continue;
1456  sortedHfEmsActive.emplace(hfem);
1457  PFClusterRef eclusterRef = elements[iHfEm].clusterRef();
1458  assert(!eclusterRef.isNull());
1459  double hfemEnergy = eclusterRef->energy();
1460  uncalibratedenergyHfEm += hfemEnergy;
1461  energyHfEm = uncalibratedenergyHfEm;
1464  uncalibratedenergyHfEm, 0.0, eclusterRef->positionREP().Eta(), eclusterRef->positionREP().Phi());
1465  energyHfHad = thepfEnergyCalibrationHF_.energyEmHad(
1466  0.0, uncalibratedenergyHfHad, hclusterRef->positionREP().Eta(), hclusterRef->positionREP().Phi());
1467  } // calib true
1468  } // loop over sortedHfEm
1469  } // if !sortedHfEms.empty()
1470  //
1471  // Create HF candidates
1472  unsigned tmpi = reconstructCluster(*hclusterRef, energyHfEm + energyHfHad);
1473  (*pfCandidates_)[tmpi].setHcalEnergy(uncalibratedenergyHfHad, energyHfHad);
1474  (*pfCandidates_)[tmpi].setEcalEnergy(uncalibratedenergyHfEm, energyHfEm);
1475  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfHad);
1476  for (auto const& hfem : sortedHfEmsActive) {
1477  unsigned iHfEm = hfem.second;
1478  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfEm);
1479  active[iHfEm] = false;
1480  }
1481 
1482  } // if sortedTracksActive.empty() ends
1483  //
1484  // Active tracks are associated.
1485  // Create HFHAD candidates from excess energy w.r.t. tracks
1486  else if ((energyHfHad - totalChargedMomentum) > nsigmaHFHAD * totalError) { // HfHad is excessive
1487  assert(energyHfEm == 0.);
1488  // HfHad candidate from excess
1489  double energyHfHadExcess = max(energyHfHad - totalChargedMomentum, 0.);
1490  double uncalibratedenergyHfHadExcess = energyHfHadExcess / calibFactorHfHad;
1491  unsigned tmpi = reconstructCluster(*hclusterRef, energyHfHadExcess);
1492  (*pfCandidates_)[tmpi].setHcalEnergy(uncalibratedenergyHfHadExcess, energyHfHadExcess);
1493  (*pfCandidates_)[tmpi].setEcalEnergy(0., 0.);
1494  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfHad);
1495  energyHfHad = max(energyHfHad - energyHfHadExcess, 0.);
1496  uncalibratedenergyHfHad = max(uncalibratedenergyHfHad - uncalibratedenergyHfHadExcess, 0.);
1497  }
1498  //
1499  // If there is a room for HFEM satellites to get associated,
1500  // loop over all HFEM satellites, starting for the closest to the various tracks
1501  // and adding other satellites until saturation of the total track momentum
1502  //
1503  else {
1504  for (auto const& hfemSatellite : hfemSatellites) {
1505  //
1506  uncalibratedenergyHfEmTmp += std::get<1>(hfemSatellite.second); // KH: raw HFEM energy
1507  energyHfEmTmp = uncalibratedenergyHfEmTmp;
1508  double energyHfHadTmp = uncalibratedenergyHfHad; // now to test hfhad calibration with EM+HAD cases
1509  unsigned iHfEm = std::get<0>(hfemSatellite.second);
1510  PFClusterRef eclusterRef = elements[iHfEm].clusterRef();
1511  assert(!eclusterRef.isNull());
1513  energyHfEmTmp = thepfEnergyCalibrationHF_.energyEmHad(
1514  uncalibratedenergyHfEmTmp, 0.0, eclusterRef->positionREP().Eta(), eclusterRef->positionREP().Phi());
1515  energyHfHadTmp = thepfEnergyCalibrationHF_.energyEmHad(
1516  0.0, uncalibratedenergyHfHad, hclusterRef->positionREP().Eta(), hclusterRef->positionREP().Phi());
1517  }
1518 
1519  double caloEnergyTmp = energyHfEmTmp + energyHfHadTmp;
1520  double calibFactorHfEm = (uncalibratedenergyHfEmTmp > 0.) ? energyHfEmTmp / uncalibratedenergyHfEmTmp : 1.;
1521 
1522  // Continue looping until all closest clusters are exhausted and as long as
1523  // the calorimetric energy does not saturate the total momentum.
1524  if (hfemSatellite.first < 0. || caloEnergyTmp < totalChargedMomentum) {
1525  LogTrace("PFAlgo|createCandidatesHF")
1526  << "\t\t\tactive, adding " << std::get<1>(hfemSatellite.second) << " to HFEM energy, and locking";
1527  active[std::get<0>(hfemSatellite.second)] = false;
1528  // HfEm is excessive (possible for the first hfemSatellite)
1529  if (hfemSatellite.first < 0. && (caloEnergyTmp - totalChargedMomentum) > nsigmaHFEM * totalError) {
1530  // HfEm candidate from excess
1531  double energyHfEmExcess = max(caloEnergyTmp - totalChargedMomentum, 0.);
1532  double uncalibratedenergyHfEmExcess = energyHfEmExcess / calibFactorHfEm;
1533  unsigned tmpi = reconstructCluster(*eclusterRef, energyHfEmExcess);
1534  (*pfCandidates_)[tmpi].setEcalEnergy(uncalibratedenergyHfEmExcess, energyHfEmExcess);
1535  (*pfCandidates_)[tmpi].setHcalEnergy(0, 0.);
1536  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfEm);
1537  energyHfEmTmp = max(energyHfEmTmp - energyHfEmExcess, 0.);
1538  uncalibratedenergyHfEmTmp = max(uncalibratedenergyHfEmTmp - uncalibratedenergyHfEmExcess, 0.);
1539  }
1540  energyHfEm = energyHfEmTmp;
1541  uncalibratedenergyHfEm = uncalibratedenergyHfEmTmp;
1542  energyHfHad = energyHfHadTmp;
1543  continue;
1544  }
1545  break;
1546  } // loop over hfemsattelites ends
1547  } // if HFHAD is excessive or not
1548 
1549  //
1550  // Loop over all tracks associated to this HFHAD cluster *again* in order to produce charged hadrons
1551  //
1552  for (auto const& trk : sortedTracksActive) {
1553  unsigned iTrack = trk.second;
1554 
1555  // Sanity check
1556  reco::TrackRef trackRef = elements[iTrack].trackRef();
1557  assert(!trackRef.isNull());
1558 
1559  //
1560  // Reconstructing charged hadrons
1561  //
1562  unsigned tmpi = reconstructTrack(elements[iTrack]);
1563  active[iTrack] = false;
1564  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfHad);
1565  auto myHfEms = associatedHfEms.equal_range(iTrack);
1566  for (auto ii = myHfEms.first; ii != myHfEms.second; ++ii) {
1567  unsigned iHfEm = ii->second.second;
1568  if (active[iHfEm])
1569  continue;
1570  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfEm);
1571  }
1572  double frac = 0.;
1573  if (totalChargedMomentum)
1574  frac = trackRef->p() / totalChargedMomentum;
1575  (*pfCandidates_)[tmpi].setEcalEnergy(uncalibratedenergyHfEm * frac, energyHfEm * frac);
1576  (*pfCandidates_)[tmpi].setHcalEnergy(uncalibratedenergyHfHad * frac, energyHfHad * frac);
1577 
1578  } // sortedTracks loop ends
1579 
1580  } // iHfHad element loop ends
1581 
1582  //
1583  // Loop over remaining active HfEm clusters
1584  //
1585  for (unsigned iHfEm = 0; iHfEm < elements.size(); iHfEm++) {
1586  PFBlockElement::Type type = elements[iHfEm].type();
1587  if (type == PFBlockElement::HFEM && active[iHfEm]) {
1588  reco::PFClusterRef eclusterRef = elements[iHfEm].clusterRef();
1589  double energyHF = 0.;
1590  double uncalibratedenergyHF = 0.;
1591  unsigned tmpi = 0;
1592  // do EM-only calibration here
1593  energyHF = eclusterRef->energy();
1594  uncalibratedenergyHF = energyHF;
1597  uncalibratedenergyHF, eclusterRef->positionREP().Eta(), eclusterRef->positionREP().Phi());
1598  }
1599  tmpi = reconstructCluster(*eclusterRef, energyHF);
1600  (*pfCandidates_)[tmpi].setEcalEnergy(uncalibratedenergyHF, energyHF);
1601  (*pfCandidates_)[tmpi].setHcalEnergy(0., 0.);
1602  (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfEm);
1603  active[iHfEm] = false;
1604  LogTrace("PFAlgo|createCandidatesHF") << "HF EM alone from blocks with tracks! " << energyHF;
1605  }
1606  } // remaining active HfEm cluster loop ends
1607 
1608  } // if-statement for blocks including tracks ends here
1609  //
1610  // -----------------------------------------------
1611  // From here, traditional PF HF candidate creation
1612  // -----------------------------------------------
1613  //
1614  else if (elements.size() == 1) {
1615  //Auguste: HAD-only calibration here
1616  reco::PFClusterRef clusterRef = elements[0].clusterRef();
1617  double energyHF = 0.;
1618  double uncalibratedenergyHF = 0.;
1619  unsigned tmpi = 0;
1620  switch (clusterRef->layer()) {
1621  case PFLayer::HF_EM:
1622  // do EM-only calibration here
1623  energyHF = clusterRef->energy();
1624  uncalibratedenergyHF = energyHF;
1627  uncalibratedenergyHF, clusterRef->positionREP().Eta(), clusterRef->positionREP().Phi());
1628  }
1629  tmpi = reconstructCluster(*clusterRef, energyHF);
1630  (*pfCandidates_)[tmpi].setEcalEnergy(uncalibratedenergyHF, energyHF);
1631  (*pfCandidates_)[tmpi].setHcalEnergy(0., 0.);
1632  (*pfCandidates_)[tmpi].setHoEnergy(0., 0.);
1633  (*pfCandidates_)[tmpi].setPs1Energy(0.);
1634  (*pfCandidates_)[tmpi].setPs2Energy(0.);
1635  (*pfCandidates_)[tmpi].addElementInBlock(blockref, inds.hfEmIs[0]);
1636  LogTrace("PFAlgo|createCandidatesHF") << "HF EM alone ! " << energyHF;
1637  break;
1638  case PFLayer::HF_HAD:
1639  // do HAD-only calibration here
1640  energyHF = clusterRef->energy();
1641  uncalibratedenergyHF = energyHF;
1644  uncalibratedenergyHF, clusterRef->positionREP().Eta(), clusterRef->positionREP().Phi());
1645  }
1646  tmpi = reconstructCluster(*clusterRef, energyHF);
1647  (*pfCandidates_)[tmpi].setHcalEnergy(uncalibratedenergyHF, energyHF);
1648  (*pfCandidates_)[tmpi].setEcalEnergy(0., 0.);
1649  (*pfCandidates_)[tmpi].setHoEnergy(0., 0.);
1650  (*pfCandidates_)[tmpi].setPs1Energy(0.);
1651  (*pfCandidates_)[tmpi].setPs2Energy(0.);
1652  (*pfCandidates_)[tmpi].addElementInBlock(blockref, inds.hfHadIs[0]);
1653  LogTrace("PFAlgo|createCandidatesHF") << "HF Had alone ! " << energyHF;
1654  break;
1655  default:
1656  assert(0);
1657  }
1658  } else if (elements.size() == 2) {
1659  //Auguste: EM + HAD calibration here
1660  reco::PFClusterRef c0 = elements[0].clusterRef();
1661  reco::PFClusterRef c1 = elements[1].clusterRef();
1662  // 2 HF elements. Must be in each layer.
1663  reco::PFClusterRef cem = (c0->layer() == PFLayer::HF_EM ? c0 : c1);
1664  reco::PFClusterRef chad = (c1->layer() == PFLayer::HF_HAD ? c1 : c0);
1665 
1666  if (cem->layer() != PFLayer::HF_EM || chad->layer() != PFLayer::HF_HAD) {
1667  edm::LogError("PFAlgo::createCandidatesHF") << "Error: 2 elements, but not 1 HFEM and 1 HFHAD";
1668  edm::LogError("PFAlgo::createCandidatesHF") << block;
1669  assert(0);
1670  // assert( c1->layer()== PFLayer::HF_EM &&
1671  // c0->layer()== PFLayer::HF_HAD );
1672  }
1673  // do EM+HAD calibration here
1674  double energyHfEm = cem->energy();
1675  double energyHfHad = chad->energy();
1676  double uncalibratedenergyHfEm = energyHfEm;
1677  double uncalibratedenergyHfHad = energyHfHad;
1680  uncalibratedenergyHfEm, 0.0, c0->positionREP().Eta(), c0->positionREP().Phi());
1681  energyHfHad = thepfEnergyCalibrationHF_.energyEmHad(
1682  0.0, uncalibratedenergyHfHad, c1->positionREP().Eta(), c1->positionREP().Phi());
1683  }
1684  auto& cand = (*pfCandidates_)[reconstructCluster(*chad, energyHfEm + energyHfHad)];
1685  cand.setEcalEnergy(uncalibratedenergyHfEm, energyHfEm);
1686  cand.setHcalEnergy(uncalibratedenergyHfHad, energyHfHad);
1687  cand.setHoEnergy(0., 0.);
1688  cand.setPs1Energy(0.);
1689  cand.setPs2Energy(0.);
1690  cand.addElementInBlock(blockref, inds.hfEmIs[0]);
1691  cand.addElementInBlock(blockref, inds.hfHadIs[0]);
1692  LogTrace("PFAlgo|createCandidatesHF") << "HF EM+HAD found ! " << energyHfEm << " " << energyHfHad;
1693  } else {
1694  // Unusual blocks including HF elements, but do not fit any of the above categories
1695  edm::LogWarning("PFAlgo::createCandidatesHF")
1696  << "Warning: HF, but n elem different from 1 or 2 or >=3 or !trackIs.empty()";
1697  edm::LogWarning("PFAlgo::createCandidatesHF") << block;
1698  }
1699  LogTrace("PFAlgo|createCandidateHF") << "end of function PFAlgo::createCandidateHF";
1700 }
unsigned reconstructCluster(const reco::PFCluster &cluster, double particleEnergy, bool useDirection=false, double particleX=0., double particleY=0., double particleZ=0.)
Definition: PFAlgo.cc:3220
Log< level::Error, false > LogError
assert(be >=bs)
double energyEmHad(double uncalibratedEnergyECAL, double uncalibratedEnergyHCAL, double eta, double phi)
#define LogTrace(id)
double hfEnergyResolution(double clusterEnergy) const
Definition: PFAlgo.cc:3378
std::vector< unsigned > hfHadIs
Definition: PFAlgo.h:50
const bool & getcalibHF_use() const
double energyEm(double uncalibratedEnergyECAL, double eta, double phi)
T sqrt(T t)
Definition: SSEVec.h:19
double nSigmaHFEM(double clusterEnergy) const
Definition: PFAlgo.cc:3390
bool isNull() const
Checks for null.
Definition: Ref.h:235
ii
Definition: cuy.py:589
double nSigmaHFHAD(double clusterEnergy) const
Definition: PFAlgo.cc:3395
unsigned reconstructTrack(const reco::PFBlockElement &elt, bool allowLoose=false)
Definition: PFAlgo.cc:3141
std::vector< unsigned > hfEmIs
Definition: PFAlgo.h:49
double energyHad(double uncalibratedEnergyHCAL, double eta, double phi)
PFEnergyCalibrationHF & thepfEnergyCalibrationHF_
Definition: PFAlgo.h:258
Log< level::Warning, false > LogWarning
std::vector< unsigned > trackIs
Definition: PFAlgo.h:45

◆ decideType()

int PFAlgo::decideType ( const edm::OwnVector< reco::PFBlockElement > &  elements,
const reco::PFBlockElement::Type  type,
std::vector< bool > &  active,
ElementIndices inds,
std::vector< bool > &  deadArea,
unsigned int  iEle 
)
private

Definition at line 1201 of file PFAlgo.cc.

References reco::CaloCluster::badHcalMarker, ECAL, ElementIndices::ecalIs, bookConverter::elements, HLT_2022v12_cff::flags, HCAL, ElementIndices::hcalIs, ElementIndices::hfEmIs, ElementIndices::hfHadIs, DigiToRawDM_cff::HO, ElementIndices::hoIs, LogTrace, ElementIndices::trackIs, and useHO_.

Referenced by elementLoop().

1206  {
1207  switch (type) {
1208  case PFBlockElement::TRACK:
1209  if (active[iEle]) {
1210  inds.trackIs.push_back(iEle);
1211  LogTrace("PFAlgo|decideType") << "TRACK, stored index, continue";
1212  }
1213  break;
1214  case PFBlockElement::ECAL:
1215  if (active[iEle]) {
1216  inds.ecalIs.push_back(iEle);
1217  LogTrace("PFAlgo|decideType") << "ECAL, stored index, continue";
1218  }
1219  return 1; //continue
1220  case PFBlockElement::HCAL:
1221  if (active[iEle]) {
1222  if (elements[iEle].clusterRef()->flags() & reco::CaloCluster::badHcalMarker) {
1223  LogTrace("PFAlgo|decideType") << "HCAL DEAD AREA: remember and skip.";
1224  active[iEle] = false;
1225  deadArea[iEle] = true;
1226  return 1; //continue
1227  }
1228  inds.hcalIs.push_back(iEle);
1229  LogTrace("PFAlgo|decideType") << "HCAL, stored index, continue";
1230  }
1231  return 1; //continue
1232  case PFBlockElement::HO:
1233  if (useHO_) {
1234  if (active[iEle]) {
1235  inds.hoIs.push_back(iEle);
1236  LogTrace("PFAlgo|decideType") << "HO, stored index, continue";
1237  }
1238  }
1239  return 1; //continue
1240  case PFBlockElement::HFEM:
1241  if (active[iEle]) {
1242  inds.hfEmIs.push_back(iEle);
1243  LogTrace("PFAlgo|decideType") << "HFEM, stored index, continue";
1244  }
1245  return 1; //continue
1246  case PFBlockElement::HFHAD:
1247  if (active[iEle]) {
1248  inds.hfHadIs.push_back(iEle);
1249  LogTrace("PFAlgo|decideType") << "HFHAD, stored index, continue";
1250  }
1251  return 1; //continue
1252  default:
1253  return 1; //continue
1254  }
1255  LogTrace("PFAlgo|decideType") << "Did not match type to anything, return 0";
1256  return 0;
1257 }
std::vector< unsigned > hoIs
Definition: PFAlgo.h:43
#define LogTrace(id)
bool useHO_
Definition: PFAlgo.h:260
std::vector< unsigned > hfHadIs
Definition: PFAlgo.h:50
std::vector< unsigned > hcalIs
Definition: PFAlgo.h:42
std::vector< unsigned > hfEmIs
Definition: PFAlgo.h:49
std::vector< unsigned > ecalIs
Definition: PFAlgo.h:44
std::vector< unsigned > trackIs
Definition: PFAlgo.h:45

◆ egammaFilters()

void PFAlgo::egammaFilters ( const reco::PFBlockRef blockref,
std::vector< bool > &  active,
PFEGammaFilters const *  pfegamma 
)
private

Definition at line 220 of file PFAlgo.cc.

References reco::LeafCandidate::charge(), reco::PFCandidate::e, reco::PFCandidate::egammaExtraRef(), reco::LeafCandidate::eta(), reco::PFCandidate::gamma, PFEGammaFilters::isElectron(), PFEGammaFilters::isElectronSafeForJetMET(), edm::Ref< C, T, F >::isNonnull(), PFEGammaFilters::isPhotonSafeForJetMET(), LogTrace, reco::PFCandidate::mva_e_pi(), nVtx_, PbPb_ZMuSkimMuonDPG_cff::particleType, PFEGammaFilters::passElectronSelection(), PFEGammaFilters::passPhotonSelection(), pfCandidates_, pfEgammaCandidates_, reco::LeafCandidate::phi(), primaryVertex_, reco::LeafCandidate::pt(), edm::View< T >::refAt(), reco::PFCandidate::set_dnn_e_bkgNonIsolated(), reco::PFCandidate::set_dnn_e_bkgPhoton(), reco::PFCandidate::set_dnn_e_bkgTau(), reco::PFCandidate::set_dnn_e_sigIsolated(), reco::PFCandidate::set_dnn_e_sigNonIsolated(), reco::PFCandidate::set_dnn_gamma(), reco::PFCandidate::set_mva_e_pi(), reco::PFCandidate::set_mva_Isolated(), reco::PFCandidate::set_mva_nothing_gamma(), reco::LeafCandidate::setCharge(), reco::LeafCandidate::setP4(), reco::PFCandidate::setParticleType(), reco::LeafCandidate::setVertex(), edm::View< T >::size(), useProtectionsForJetMET_, findQualityFiles::v, reco::Vertex::x(), reco::Vertex::y(), and reco::Vertex::z().

Referenced by processBlock().

222  {
223  // const edm::ValueMap<reco::GsfElectronRef> & myGedElectronValMap(*valueMapGedElectrons_);
224 
225  unsigned int negmcandidates = pfEgammaCandidates_->size();
226  LogTrace("PFAlgo|egammaFilters") << "start of function PFAlgo::egammaFilters(), negmcandidates=" << negmcandidates;
227 
228  for (unsigned int ieg = 0; ieg < negmcandidates; ++ieg) {
229  // const reco::PFCandidate & egmcand((*pfEgammaCandidates_)[ieg]);
231 
232  const PFCandidate::ElementsInBlocks& theElements = (*pfEgmRef).elementsInBlocks();
233  PFCandidate::ElementsInBlocks::const_iterator iegfirst = theElements.begin();
234  bool sameBlock = false;
235  bool isGoodElectron = false;
236  bool isGoodPhoton = false;
237  bool isPrimaryElectron = false;
238  if (iegfirst->first == blockref)
239  sameBlock = true;
240  if (sameBlock) {
241  LogTrace("PFAlgo|egammaFilters") << " I am in looping on EGamma Candidates: pt " << (*pfEgmRef).pt()
242  << " eta,phi " << (*pfEgmRef).eta() << ", " << (*pfEgmRef).phi() << " charge "
243  << (*pfEgmRef).charge();
244 
245  if ((*pfEgmRef).gsfTrackRef().isNonnull()) {
246  reco::GsfElectronRef gedEleRef = (*valueMapGedElectrons_)[pfEgmRef];
247  if (gedEleRef.isNonnull()) {
248  isGoodElectron = pfegamma->passElectronSelection(*gedEleRef, *pfEgmRef, nVtx_);
249  isPrimaryElectron = pfegamma->isElectron(*gedEleRef);
250  if (isGoodElectron)
251  LogTrace("PFAlgo|egammaFilters")
252  << "** Good Electron, pt " << gedEleRef->pt() << " eta, phi " << gedEleRef->eta() << ", "
253  << gedEleRef->phi() << " charge " << gedEleRef->charge() << " isPrimary " << isPrimaryElectron
254  << " isoDr03 "
255  << (gedEleRef->dr03TkSumPt() + gedEleRef->dr03EcalRecHitSumEt() + gedEleRef->dr03HcalTowerSumEt())
256  << " mva_isolated " << gedEleRef->mva_Isolated() << " mva_e_pi " << gedEleRef->mva_e_pi();
257  }
258  }
259  if ((*pfEgmRef).superClusterRef().isNonnull()) {
260  reco::PhotonRef gedPhoRef = (*valueMapGedPhotons_)[pfEgmRef];
261  if (gedPhoRef.isNonnull()) {
262  isGoodPhoton = pfegamma->passPhotonSelection(*gedPhoRef);
263  if (isGoodPhoton)
264  LogTrace("PFAlgo|egammaFilters") << "** Good Photon, pt " << gedPhoRef->pt() << " eta, phi "
265  << gedPhoRef->eta() << ", " << gedPhoRef->phi() << endl;
266  }
267  }
268  } // end sameBlock
269 
270  if (isGoodElectron && isGoodPhoton) {
271  if (isPrimaryElectron)
272  isGoodPhoton = false;
273  else
274  isGoodElectron = false;
275  }
276 
277  // isElectron
278  if (isGoodElectron) {
279  reco::GsfElectronRef gedEleRef = (*valueMapGedElectrons_)[pfEgmRef];
280  reco::PFCandidate myPFElectron = *pfEgmRef;
281  // run protections
282  bool lockTracks = false;
283  bool isSafe = true;
285  lockTracks = true;
286  isSafe = pfegamma->isElectronSafeForJetMET(*gedEleRef, myPFElectron, primaryVertex_, lockTracks);
287  }
288 
289  if (isSafe) {
291  myPFElectron.setParticleType(particleType);
292  myPFElectron.setCharge(gedEleRef->charge());
293  myPFElectron.setP4(gedEleRef->p4());
294  myPFElectron.set_mva_e_pi(gedEleRef->mva_e_pi());
295  myPFElectron.set_mva_Isolated(gedEleRef->mva_Isolated());
296 
297  myPFElectron.set_dnn_e_sigIsolated(gedEleRef->dnn_signal_Isolated());
298  myPFElectron.set_dnn_e_sigNonIsolated(gedEleRef->dnn_signal_nonIsolated());
299  myPFElectron.set_dnn_e_bkgNonIsolated(gedEleRef->dnn_bkg_nonIsolated());
300  myPFElectron.set_dnn_e_bkgTau(gedEleRef->dnn_bkg_Tau());
301  myPFElectron.set_dnn_e_bkgPhoton(gedEleRef->dnn_bkg_Photon());
302 
303  LogTrace("PFAlgo|egammaFilters") << " PFAlgo: found an electron with NEW EGamma code ";
304  LogTrace("PFAlgo|egammaFilters") << " myPFElectron: pt " << myPFElectron.pt() << " eta,phi "
305  << myPFElectron.eta() << ", " << myPFElectron.phi() << " mva "
306  << myPFElectron.mva_e_pi() << " charge " << myPFElectron.charge();
307 
308  LogTrace("PFAlgo|egammaFilters") << " THE BLOCK " << *blockref;
309  for (auto const& eb : theElements) {
310  active[eb.second] = false;
311  LogTrace("PFAlgo|egammaFilters") << " Elements used " << eb.second;
312  }
313 
314  // The electron is considered safe for JetMET and the additional tracks pointing to it are locked
315  if (lockTracks) {
316  const PFCandidate::ElementsInBlocks& extraTracks = myPFElectron.egammaExtraRef()->extraNonConvTracks();
317  for (auto const& trk : extraTracks) {
318  LogTrace("PFAlgo|egammaFilters") << " Extra locked track " << trk.second;
319  active[trk.second] = false;
320  }
321  }
322 
323  LogTrace("PFAlgo|egammaFilters") << "Creating PF electron: pt=" << myPFElectron.pt()
324  << " eta=" << myPFElectron.eta() << " phi=" << myPFElectron.phi();
325  pfCandidates_->push_back(myPFElectron);
326 
327  } else {
328  LogTrace("PFAlgo|egammaFilters") << "PFAlgo: Electron DISCARDED, NOT SAFE FOR JETMET ";
329  }
330  } //isGoodElectron
331 
332  if (isGoodPhoton) {
333  reco::PhotonRef gedPhoRef = (*valueMapGedPhotons_)[pfEgmRef];
334  reco::PFCandidate myPFPhoton = *pfEgmRef;
335  bool isSafe = true;
337  isSafe = pfegamma->isPhotonSafeForJetMET(*gedPhoRef, myPFPhoton);
338  }
339 
340  if (isSafe) {
342  myPFPhoton.setParticleType(particleType);
343  myPFPhoton.setCharge(0);
344  myPFPhoton.set_mva_nothing_gamma(1.);
346  myPFPhoton.setVertex(v);
347  myPFPhoton.setP4(gedPhoRef->p4());
348  // DNN pfid
349  myPFPhoton.set_dnn_gamma(gedPhoRef->pfDNN());
350  LogTrace("PFAlgo|egammaFilters") << " PFAlgo: found a photon with NEW EGamma code ";
351  LogTrace("PFAlgo|egammaFilters") << " myPFPhoton: pt " << myPFPhoton.pt() << " eta,phi " << myPFPhoton.eta()
352  << ", " << myPFPhoton.phi() << " charge " << myPFPhoton.charge();
353 
354  // Lock all the elements
355  LogTrace("PFAlgo|egammaFilters") << " THE BLOCK " << *blockref;
356  for (auto const& eb : theElements) {
357  active[eb.second] = false;
358  LogTrace("PFAlgo|egammaFilters") << " Elements used " << eb.second;
359  }
360  LogTrace("PFAlgo|egammaFilters") << "Creating PF photon: pt=" << myPFPhoton.pt() << " eta=" << myPFPhoton.eta()
361  << " phi=" << myPFPhoton.phi();
362  pfCandidates_->push_back(myPFPhoton);
363 
364  } // end isSafe
365  } // end isGoodPhoton
366  } // end loop on EGM candidates
367  LogTrace("PFAlgo|egammaFilters") << "end of function PFAlgo::egammaFilters";
368 }
void set_dnn_e_sigIsolated(float mva)
Definition: PFCandidate.h:349
ParticleType
particle types
Definition: PFCandidate.h:44
RefToBase< value_type > refAt(size_type i) const
double pt() const final
transverse momentum
double z() const
z coordinate
Definition: Vertex.h:133
void set_mva_nothing_gamma(float mva)
set mva for gamma detection
Definition: PFCandidate.h:332
std::unique_ptr< reco::PFCandidateCollection > pfCandidates_
Definition: PFAlgo.h:226
bool useProtectionsForJetMET_
Definition: PFAlgo.h:266
void set_dnn_gamma(float mva)
Definition: PFCandidate.h:369
void set_dnn_e_bkgNonIsolated(float mva)
Definition: PFCandidate.h:357
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
int nVtx_
Definition: PFAlgo.h:286
#define LogTrace(id)
void setVertex(const Point &vertex) override
set vertex
size_type size() const
std::vector< ElementInBlock > ElementsInBlocks
Definition: PFCandidate.h:404
reco::PFCandidateEGammaExtraRef egammaExtraRef() const
return a reference to the EGamma extra
Definition: PFCandidate.cc:576
void setCharge(Charge q) final
set electric charge
void set_mva_e_pi(float mvaNI)
Definition: PFCandidate.h:315
void set_mva_Isolated(float mvaI)
Definition: PFCandidate.h:311
const edm::View< reco::PFCandidate > * pfEgammaCandidates_
Definition: PFAlgo.h:267
void setParticleType(ParticleType type)
set Particle Type
Definition: PFCandidate.cc:278
reco::Vertex primaryVertex_
Definition: PFAlgo.h:328
double x() const
x coordinate
Definition: Vertex.h:129
double y() const
y coordinate
Definition: Vertex.h:131
void set_dnn_e_sigNonIsolated(float mva)
Definition: PFCandidate.h:353
void set_dnn_e_bkgPhoton(float mva)
Definition: PFCandidate.h:365
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
void set_dnn_e_bkgTau(float mva)
Definition: PFCandidate.h:361
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:41
double phi() const final
momentum azimuthal angle
void setP4(const LorentzVector &p4) final
set 4-momentum
float mva_e_pi() const
mva for electron-pion discrimination
Definition: PFCandidate.h:317
int charge() const final
electric charge
double eta() const final
momentum pseudorapidity

◆ elementLoop()

void PFAlgo::elementLoop ( const reco::PFBlock block,
reco::PFBlock::LinkData linkData,
const edm::OwnVector< reco::PFBlockElement > &  elements,
std::vector< bool > &  active,
const reco::PFBlockRef blockref,
ElementIndices inds,
std::vector< bool > &  deadArea 
)
private

Definition at line 997 of file PFAlgo.cc.

References cms::cuda::assert(), groupFilesInBlocks::block, checkAndReconstructSecondaryInteraction(), checkGoodTrackDeadHcal(), checkHasDeadHcal(), decideType(), ECAL, reco::PFBlockElement::ECAL, ElementIndices::ecalIs, bookConverter::elements, reco::PFBlockElement::GSF, HCAL, reco::PFBlockElement::HCAL, ElementIndices::hcalIs, reco::PFBlockElement::HFEM, ElementIndices::hfEmIs, reco::PFBlockElement::HFHAD, ElementIndices::hfHadIs, ElementIndices::hoIs, edm::Ref< C, T, F >::isNull(), reco::PFBlock::LINKTEST_ALL, LogTrace, recoTracksNotHCAL(), relinkTrackToHcal(), and ElementIndices::trackIs.

Referenced by processBlock().

1003  {
1004  LogTrace("PFAlgo|elementLoop") << "start of function PFAlgo::elementLoop, elements.size()" << elements.size();
1005 
1006  for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
1007  PFBlockElement::Type type = elements[iEle].type();
1008 
1009  LogTrace("PFAlgo|elementLoop") << "elements[iEle=" << iEle << "]=" << elements[iEle];
1010  //only process TRACK elements, but fill the ElementIndices vector with indices for all elements.
1011  //Mark the active & deadArea for bad HCAL
1012  auto ret_decideType = decideType(elements, type, active, inds, deadArea, iEle);
1013  if (ret_decideType == 1) {
1014  LogTrace("PFAlgo|elementLoop") << "ret_decideType==1, continuing";
1015  continue;
1016  }
1017  LogTrace("PFAlgo|elementLoop") << "ret_decideType=" << ret_decideType << " type=" << type;
1018 
1019  active[iEle] = checkAndReconstructSecondaryInteraction(blockref, elements, active[iEle], iEle);
1020 
1021  if (!active[iEle]) {
1022  LogTrace("PFAlgo|elementLoop") << "Already used by electrons, muons, conversions";
1023  continue;
1024  }
1025 
1026  reco::TrackRef trackRef = elements[iEle].trackRef();
1027  assert(!trackRef.isNull());
1028 
1029  LogTrace("PFAlgo|elementLoop") << "PFAlgo:processBlock"
1030  << " trackIs.size()=" << inds.trackIs.size()
1031  << " ecalIs.size()=" << inds.ecalIs.size() << " hcalIs.size()=" << inds.hcalIs.size()
1032  << " hoIs.size()=" << inds.hoIs.size() << " hfEmIs.size()=" << inds.hfEmIs.size()
1033  << " hfHadIs.size()=" << inds.hfHadIs.size();
1034 
1035  // look for associated elements of all types
1036  //COLINFEB16
1037  // all types of links are considered.
1038  // the elements are sorted by increasing distance
1039  std::multimap<double, unsigned> ecalElems;
1040  block.associatedElements(iEle, linkData, ecalElems, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
1041 
1042  std::multimap<double, unsigned> hcalElems;
1043  block.associatedElements(iEle, linkData, hcalElems, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
1044 
1045  std::multimap<double, unsigned> hfEmElems;
1046  std::multimap<double, unsigned> hfHadElems;
1047  block.associatedElements(iEle, linkData, hfEmElems, reco::PFBlockElement::HFEM, reco::PFBlock::LINKTEST_ALL);
1048  block.associatedElements(iEle, linkData, hfHadElems, reco::PFBlockElement::HFHAD, reco::PFBlock::LINKTEST_ALL);
1049 
1050  LogTrace("PFAlgo|elementLoop") << "\tTrack " << iEle << " is linked to " << ecalElems.size() << " ecal and "
1051  << hcalElems.size() << " hcal and " << hfEmElems.size() << " hfEm and "
1052  << hfHadElems.size() << " hfHad elements";
1053 
1054 #ifdef EDM_ML_DEBUG
1055  for (const auto& pair : ecalElems) {
1056  LogTrace("PFAlgo|elementLoop") << "ecal: dist " << pair.first << "\t elem " << pair.second;
1057  }
1058  for (const auto& pair : hcalElems) {
1059  LogTrace("PFAlgo|elementLoop") << "hcal: dist " << pair.first << "\t elem " << pair.second
1060  << (deadArea[pair.second] ? " DEAD AREA MARKER" : "");
1061  }
1062 #endif
1063 
1064  const bool hasDeadHcal = checkHasDeadHcal(hcalElems, deadArea);
1065  if (hasDeadHcal) {
1066  hcalElems.clear();
1067  }
1068  const bool goodTrackDeadHcal = checkGoodTrackDeadHcal(trackRef, hasDeadHcal);
1069 
1070  // When a track has no HCAL cluster linked, but another track is linked to the same
1071  // ECAL cluster and an HCAL cluster, link the track to the HCAL cluster for
1072  // later analysis
1073  if (hcalElems.empty() && !ecalElems.empty() && !hasDeadHcal) {
1074  relinkTrackToHcal(block, ecalElems, hcalElems, active, linkData, iEle);
1075  }
1076 
1077  //MICHELE
1078  //TEMPORARY SOLUTION FOR ELECTRON REJECTION IN PFTAU
1079  //COLINFEB16
1080  // in case particle flow electrons are not reconstructed,
1081  // the mva_e_pi of the charged hadron will be set to 1
1082  // if a GSF element is associated to the current TRACK element
1083  // This information will be used in the electron rejection for tau ID.
1084  std::multimap<double, unsigned> gsfElems;
1085  block.associatedElements(iEle, linkData, gsfElems, reco::PFBlockElement::GSF);
1086 
1087  if (hcalElems.empty()) {
1088  LogTrace("PFAlgo|elementLoop") << "no hcal element connected to track " << iEle;
1089  }
1090 
1091  // will now loop on associated elements ...
1092  bool hcalFound = false;
1093  bool hfhadFound = false;
1094 
1095  LogTrace("PFAlgo|elementLoop") << "now looping on elements associated to the track: ecalElems";
1096 
1097  // ... first on associated ECAL elements
1098  // Check if there is still a free ECAL for this track
1099  for (auto const& ecal : ecalElems) {
1100  unsigned index = ecal.second;
1101  // Sanity checks and optional printout
1103 #ifdef EDM_ML_DEBUG
1104  double dist = ecal.first;
1105  LogTrace("PFAlgo|elementLoop") << "\telement " << elements[index] << " linked with distance = " << dist;
1106  if (!active[index])
1107  LogTrace("PFAlgo|elementLoop") << "This ECAL is already used - skip it";
1108 #endif
1110 
1111  // This ECAL is not free (taken by an electron?) - just skip it
1112  if (!active[index])
1113  continue;
1114 
1115  // Flag ECAL clusters for which the corresponding track is not linked to an HCAL cluster
1116 
1117  //reco::PFClusterRef clusterRef = elements[index].clusterRef();
1118  //assert( !clusterRef.isNull() );
1119  if (!hcalElems.empty())
1120  LogTrace("PFAlgo|elementLoop") << "\t\tat least one hcal element connected to the track."
1121  << " Sparing Ecal cluster for the hcal loop";
1122 
1123  } //loop print ecal elements
1124 
1125  // tracks which are not linked to an HCAL (or HFHAD)
1126  // are reconstructed now.
1127 
1128  if (hcalElems.empty() && hfHadElems.empty()) {
1129  auto ret_continue = recoTracksNotHCAL(
1130  block, linkData, elements, blockref, active, goodTrackDeadHcal, hasDeadHcal, iEle, ecalElems, trackRef);
1131  if (ret_continue) {
1132  continue;
1133  }
1134  } // end if( hcalElems.empty() && hfHadElems.empty() )
1135 
1136  // In case several HCAL elements are linked to this track,
1137  // unlinking all of them except the closest.
1138  for (auto const& hcal : hcalElems) {
1139  unsigned index = hcal.second;
1140 
1142 
1143 #ifdef EDM_ML_DEBUG
1144  double dist = block.dist(iEle, index, linkData, reco::PFBlock::LINKTEST_ALL);
1145  LogTrace("PFAlgo|elementLoop") << "\telement " << elements[index] << " linked with distance " << dist;
1146 #endif
1148 
1149  // all hcal clusters except the closest
1150  // will be unlinked from the track
1151  if (!hcalFound) { // closest hcal
1152  LogTrace("PFAlgo|elementLoop") << "\t\tclosest hcal cluster, doing nothing";
1153 
1154  hcalFound = true;
1155 
1156  // active[index] = false;
1157  // hcalUsed.push_back( index );
1158  } else { // other associated hcal
1159  // unlink from the track
1160  LogTrace("PFAlgo|elementLoop") << "\t\tsecondary hcal cluster. unlinking";
1161  block.setLink(iEle, index, -1., linkData, PFBlock::LINKTEST_RECHIT);
1162  }
1163  } //loop hcal elements
1164 
1165  // ---Same for HFHAD---
1166  // In case several HFHAD elements are linked to this track,
1167  // unlinking all of them except the closest.
1168  for (auto const& hfhad : hfHadElems) {
1169  unsigned index = hfhad.second;
1170 
1172 
1173 #ifdef EDM_ML_DEBUG
1174  double dist = block.dist(iEle, index, linkData, reco::PFBlock::LINKTEST_ALL);
1175  LogTrace("PFAlgo|elementLoop") << "\telement " << elements[index] << " linked with distance " << dist;
1176 #endif
1177  assert(type == PFBlockElement::HFHAD);
1178 
1179  // all hfhad clusters except the closest
1180  // will be unlinked from the track
1181  if (!hfhadFound) { // closest hfhad
1182  LogTrace("PFAlgo|elementLoop") << "\t\tclosest hfhad cluster, doing nothing";
1183 
1184  hfhadFound = true;
1185 
1186  } else { // other associated hfhad
1187  // unlink from the track
1188  LogTrace("PFAlgo|elementLoop") << "\t\tsecondary hfhad cluster. unlinking";
1189  block.setLink(iEle, index, -1., linkData, PFBlock::LINKTEST_RECHIT);
1190  }
1191  } //loop hfhad elements
1192 
1193  LogTrace("PFAlgo|elementLoop") << "end of loop over iEle";
1194  } // end of outer loop on elements iEle of any type
1195  LogTrace("PFAlgo|elementLoop") << "end of function PFAlgo::elementLoop";
1196 }
bool checkHasDeadHcal(const std::multimap< double, unsigned > &hcalElems, const std::vector< bool > &deadArea)
Definition: PFAlgo.cc:875
std::vector< unsigned > hoIs
Definition: PFAlgo.h:43
bool checkGoodTrackDeadHcal(const reco::TrackRef &trackRef, bool hasDeadHcal)
Definition: PFAlgo.cc:902
assert(be >=bs)
#define LogTrace(id)
std::vector< unsigned > hfHadIs
Definition: PFAlgo.h:50
std::vector< unsigned > hcalIs
Definition: PFAlgo.h:42
bool recoTracksNotHCAL(const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, const reco::PFBlockRef &blockref, std::vector< bool > &active, bool goodTrackDeadHcal, bool hasDeadHcal, unsigned int iTrack, std::multimap< double, unsigned > &ecalElems, reco::TrackRef &trackRef)
Definition: PFAlgo.cc:399
bool isNull() const
Checks for null.
Definition: Ref.h:235
int decideType(const edm::OwnVector< reco::PFBlockElement > &elements, const reco::PFBlockElement::Type type, std::vector< bool > &active, ElementIndices &inds, std::vector< bool > &deadArea, unsigned int iEle)
Definition: PFAlgo.cc:1201
std::vector< unsigned > hfEmIs
Definition: PFAlgo.h:49
std::vector< unsigned > ecalIs
Definition: PFAlgo.h:44
void relinkTrackToHcal(const reco::PFBlock &block, std::multimap< double, unsigned > &ecalElems, std::multimap< double, unsigned > &hcalElems, const std::vector< bool > &active, reco::PFBlock::LinkData &linkData, unsigned int iTrack)
Definition: PFAlgo.cc:941
std::vector< unsigned > trackIs
Definition: PFAlgo.h:45
bool checkAndReconstructSecondaryInteraction(const reco::PFBlockRef &blockref, const edm::OwnVector< reco::PFBlockElement > &elements, bool isActive, int iElement)
Definition: PFAlgo.cc:855

◆ getCleanedCandidates()

reco::PFCandidateCollection& PFAlgo::getCleanedCandidates ( )
inline
Returns
collection of cleaned HF candidates

Definition at line 112 of file PFAlgo.h.

References pfCleanedCandidates_.

Referenced by PFProducer::produce().

112 { return pfCleanedCandidates_; }
reco::PFCandidateCollection pfCleanedCandidates_
Definition: PFAlgo.h:228

◆ getPFMuonAlgo()

PFMuonAlgo * PFAlgo::getPFMuonAlgo ( )

Definition at line 66 of file PFAlgo.cc.

References pfmu_.

Referenced by PFProducer::produce().

66 { return pfmu_.get(); }
std::unique_ptr< PFMuonAlgo > pfmu_
Definition: PFAlgo.h:262

◆ hfEnergyResolution()

double PFAlgo::hfEnergyResolution ( double  clusterEnergy) const
private

Definition at line 3378 of file PFAlgo.cc.

References SiStripPI::max, funct::pow(), pfClustersFromCombinedCaloHCal_cfi::resol, resolHF_square_, and mathSSE::sqrt().

Referenced by createCandidatesHF().

3378  {
3379  // Add a protection
3380  clusterEnergyHF = std::max(clusterEnergyHF, 1.);
3381 
3382  double resol =
3383  sqrt(resolHF_square_[0] / clusterEnergyHF + resolHF_square_[1] + resolHF_square_[2] / pow(clusterEnergyHF, 2));
3384  // 0: stochastic term, 1: constant term, 2: noise term
3385  // Note: resolHF_square_[0,1,2] should be already squared
3386 
3387  return resol;
3388 }
const std::vector< double > resolHF_square_
Definition: PFAlgo.h:255
T sqrt(T t)
Definition: SSEVec.h:19
Power< A, B >::type pow(const A &a, const B &b)
Definition: Power.h:29

◆ isFromSecInt()

bool PFAlgo::isFromSecInt ( const reco::PFBlockElement eTrack,
std::string  order 
) const
private

Definition at line 3470 of file PFAlgo.cc.

References eventshapeDQM_cfi::order, reco::PFBlockElement::T_FROM_DISP, reco::PFBlockElement::T_FROM_V0, reco::PFBlockElement::T_TO_DISP, reco::PFBlockElement::trackType(), usePFDecays_, and usePFNuclearInteractions_.

Referenced by checkAndReconstructSecondaryInteraction(), createCandidatesHCAL(), reconstructTrack(), and recoTracksNotHCAL().

3470  {
3473  // reco::PFBlockElement::TrackType T_FROM_GAMMACONV = reco::PFBlockElement::T_FROM_GAMMACONV;
3475 
3476  bool bPrimary = (order.find("primary") != string::npos);
3477  bool bSecondary = (order.find("secondary") != string::npos);
3478  bool bAll = (order.find("all") != string::npos);
3479 
3480  bool isToDisp = usePFNuclearInteractions_ && eTrack.trackType(T_TO_DISP);
3481  bool isFromDisp = usePFNuclearInteractions_ && eTrack.trackType(T_FROM_DISP);
3482 
3483  if (bPrimary && isToDisp)
3484  return true;
3485  if (bSecondary && isFromDisp)
3486  return true;
3487  if (bAll && (isToDisp || isFromDisp))
3488  return true;
3489 
3490  // bool isFromConv = usePFConversions_ && eTrack.trackType(T_FROM_GAMMACONV);
3491 
3492  // if ((bAll || bSecondary)&& isFromConv) return true;
3493 
3494  bool isFromDecay = (bAll || bSecondary) && usePFDecays_ && eTrack.trackType(T_FROM_V0);
3495 
3496  return isFromDecay;
3497 }
bool usePFNuclearInteractions_
Definition: PFAlgo.h:279
virtual bool trackType(TrackType trType) const
bool usePFDecays_
Definition: PFAlgo.h:281

◆ makeConnectedCandidates()

reco::PFCandidateCollection PFAlgo::makeConnectedCandidates ( )
inline
Returns
the collection of candidates

Definition at line 115 of file PFAlgo.h.

References PFCandConnector::connect(), connector_, and pfCandidates_.

Referenced by PFProducer::produce().

115 { return connector_.connect(*pfCandidates_); }
std::unique_ptr< reco::PFCandidateCollection > pfCandidates_
Definition: PFAlgo.h:226
PFCandConnector connector_
Definition: PFAlgo.h:290
reco::PFCandidateCollection connect(reco::PFCandidateCollection &pfCand) const

◆ neutralHadronEnergyResolution()

double PFAlgo::neutralHadronEnergyResolution ( double  clusterEnergy,
double  clusterEta 
) const
private

todo: use PFClusterTools for this

Definition at line 3361 of file PFAlgo.cc.

References PVValHelper::eta, SiStripPI::max, pfClustersFromCombinedCaloHCal_cfi::resol, and mathSSE::sqrt().

Referenced by createCandidatesHCAL(), and recoTracksNotHCAL().

3361  {
3362  // Add a protection
3363  clusterEnergyHCAL = std::max(clusterEnergyHCAL, 1.);
3364 
3365  double resol = fabs(eta) < 1.48 ? sqrt(1.02 * 1.02 / clusterEnergyHCAL + 0.065 * 0.065)
3366  : sqrt(1.20 * 1.20 / clusterEnergyHCAL + 0.028 * 0.028);
3367 
3368  return resol;
3369 }
T sqrt(T t)
Definition: SSEVec.h:19

◆ nSigmaHCAL()

double PFAlgo::nSigmaHCAL ( double  clusterEnergy,
double  clusterEta 
) const
private

Definition at line 3371 of file PFAlgo.cc.

References PVValHelper::eta, JetChargeProducer_cfi::exp, nSigmaEConstHCAL, and nSigmaHCAL_.

Referenced by createCandidatesHCAL().

3371  {
3372  double nS = fabs(eta) < 1.48 ? nSigmaHCAL_ * (1. + exp(-clusterEnergyHCAL / nSigmaEConstHCAL))
3373  : nSigmaHCAL_ * (1. + exp(-clusterEnergyHCAL / nSigmaEConstHCAL));
3374 
3375  return nS;
3376 }
const double nSigmaHCAL_
number of sigma to judge energy excess in HCAL
Definition: PFAlgo.h:248
const double nSigmaEConstHCAL
Definition: PFAlgo.h:334

◆ nSigmaHFEM()

double PFAlgo::nSigmaHFEM ( double  clusterEnergy) const
private

Definition at line 3390 of file PFAlgo.cc.

References JetChargeProducer_cfi::exp, nSigmaEConstHFEM, and nSigmaHFEM_.

Referenced by createCandidatesHF().

3390  {
3391  double nS = nSigmaHFEM_ * (1. + exp(-clusterEnergyHF / nSigmaEConstHFEM));
3392  return nS;
3393 }
const double nSigmaHFEM_
number of sigma to judge energy excess in HF
Definition: PFAlgo.h:251
const double nSigmaEConstHFEM
Definition: PFAlgo.h:335

◆ nSigmaHFHAD()

double PFAlgo::nSigmaHFHAD ( double  clusterEnergy) const
private

Definition at line 3395 of file PFAlgo.cc.

References JetChargeProducer_cfi::exp, nSigmaEConstHFHAD, and nSigmaHFHAD_.

Referenced by createCandidatesHF().

3395  {
3396  double nS = nSigmaHFHAD_ * (1. + exp(-clusterEnergyHF / nSigmaEConstHFHAD));
3397  return nS;
3398 }
const double nSigmaEConstHFHAD
Definition: PFAlgo.h:336
const double nSigmaHFHAD_
Definition: PFAlgo.h:252

◆ postCleaning()

void PFAlgo::postCleaning ( )
private

Definition at line 3499 of file PFAlgo.cc.

References SiPixelRawToDigiRegional_cfi::deltaPhi, reco::PFCandidate::egamma_HF, reco::PFCandidate::h_HF, mps_fire::i, dqmiolumiharvest::j, maxDeltaPhiPt_, maxSignificance_, minDeltaMet_, minHFCleaningPt_, minSignificance_, minSignificanceReduction_, GetRecoTauVFromDQM_MC_cff::next, reco::PFCandidate::particleId(), pfCandidates_, pfCleanedCandidates_, reco::LeafCandidate::pt(), reco::LeafCandidate::px(), reco::LeafCandidate::py(), met_cff::significance, optionsL1T::skip, and mathSSE::sqrt().

Referenced by reconstructParticles().

3499  {
3500  //Compute met and met significance (met/sqrt(SumEt))
3501  double metX = 0.;
3502  double metY = 0.;
3503  double sumet = 0;
3504  std::vector<unsigned int> pfCandidatesToBeRemoved;
3505  for (auto const& pfc : *pfCandidates_) {
3506  metX += pfc.px();
3507  metY += pfc.py();
3508  sumet += pfc.pt();
3509  }
3510  double met2 = metX * metX + metY * metY;
3511  // Select events with large MET significance.
3512  double significance = std::sqrt(met2 / sumet);
3513  double significanceCor = significance;
3515  double metXCor = metX;
3516  double metYCor = metY;
3517  double sumetCor = sumet;
3518  double met2Cor = met2;
3519  double deltaPhi = 3.14159;
3520  double deltaPhiPt = 100.;
3521  bool next = true;
3522  unsigned iCor = 1E9;
3523 
3524  // Find the HF candidate with the largest effect on the MET
3525  while (next) {
3526  double metReduc = -1.;
3527  // Loop on the candidates
3528  for (unsigned i = 0; i < pfCandidates_->size(); ++i) {
3529  const PFCandidate& pfc = (*pfCandidates_)[i];
3530 
3531  // Check that the pfCandidate is in the HF
3533  continue;
3534 
3535  // Check if has meaningful pt
3536  if (pfc.pt() < minHFCleaningPt_)
3537  continue;
3538 
3539  // Check that it is not already scheculed to be cleaned
3540  const bool skip = std::any_of(
3541  pfCandidatesToBeRemoved.begin(), pfCandidatesToBeRemoved.end(), [&](unsigned int j) { return i == j; });
3542  if (skip)
3543  continue;
3544 
3545  // Check that the pt and the MET are aligned
3546  deltaPhi = std::acos((metX * pfc.px() + metY * pfc.py()) / (pfc.pt() * std::sqrt(met2)));
3547  deltaPhiPt = deltaPhi * pfc.pt();
3548  if (deltaPhiPt > maxDeltaPhiPt_)
3549  continue;
3550 
3551  // Now remove the candidate from the MET
3552  double metXInt = metX - pfc.px();
3553  double metYInt = metY - pfc.py();
3554  double sumetInt = sumet - pfc.pt();
3555  double met2Int = metXInt * metXInt + metYInt * metYInt;
3556  if (met2Int < met2Cor) {
3557  metXCor = metXInt;
3558  metYCor = metYInt;
3559  metReduc = (met2 - met2Int) / met2Int;
3560  met2Cor = met2Int;
3561  sumetCor = sumetInt;
3562  significanceCor = std::sqrt(met2Cor / sumetCor);
3563  iCor = i;
3564  }
3565  }
3566  //
3567  // If the MET must be significanly reduced, schedule the candidate to be cleaned
3568  if (metReduc > minDeltaMet_) {
3569  pfCandidatesToBeRemoved.push_back(iCor);
3570  metX = metXCor;
3571  metY = metYCor;
3572  sumet = sumetCor;
3573  met2 = met2Cor;
3574  } else {
3575  // Otherwise just stop the loop
3576  next = false;
3577  }
3578  }
3579  //
3580  // The significance must be significantly reduced to indeed clean the candidates
3581  if (significance - significanceCor > minSignificanceReduction_ && significanceCor < maxSignificance_) {
3582  edm::LogInfo("PFAlgo|postCleaning") << "Significance reduction = " << significance << " -> " << significanceCor
3583  << " = " << significanceCor - significance;
3584  for (unsigned int toRemove : pfCandidatesToBeRemoved) {
3585  edm::LogInfo("PFAlgo|postCleaning") << "Removed : " << (*pfCandidates_)[toRemove];
3586  pfCleanedCandidates_.push_back((*pfCandidates_)[toRemove]);
3587  (*pfCandidates_)[toRemove].rescaleMomentum(1E-6);
3588  //reco::PFCandidate::ParticleType unknown = reco::PFCandidate::X;
3589  //(*pfCandidates_)[toRemove].setParticleType(unknown);
3590  }
3591  }
3592  } //significance
3593 } //postCleaning
double maxDeltaPhiPt_
Definition: PFAlgo.h:323
double maxSignificance_
Definition: PFAlgo.h:321
double pt() const final
transverse momentum
double minHFCleaningPt_
Definition: PFAlgo.h:319
double minDeltaMet_
Definition: PFAlgo.h:324
double minSignificance_
Definition: PFAlgo.h:320
std::unique_ptr< reco::PFCandidateCollection > pfCandidates_
Definition: PFAlgo.h:226
reco::PFCandidateCollection pfCleanedCandidates_
Definition: PFAlgo.h:228
double minSignificanceReduction_
Definition: PFAlgo.h:322
double px() const final
x coordinate of momentum vector
T sqrt(T t)
Definition: SSEVec.h:19
double py() const final
y coordinate of momentum vector
Log< level::Info, false > LogInfo
significance
Definition: met_cff.py:18
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:41
virtual ParticleType particleId() const
Definition: PFCandidate.h:392

◆ processBlock()

void PFAlgo::processBlock ( const reco::PFBlockRef blockref,
std::list< reco::PFBlockRef > &  hcalBlockRefs,
std::list< reco::PFBlockRef > &  ecalBlockRefs,
PFEGammaFilters const *  pfegamma 
)
private

process one block. can be reimplemented in more sophisticated algorithms

Definition at line 3061 of file PFAlgo.cc.

References cms::cuda::assert(), groupFilesInBlocks::block, conversionAlgo(), createCandidatesECAL(), createCandidatesHCAL(), createCandidatesHCALUnlinked(), createCandidatesHF(), egammaFilters(), elementLoop(), bookConverter::elements, edm::Ref< C, T, F >::isNull(), LogDebug, LogTrace, useEGammaFilters_, and usePFConversions_.

Referenced by reconstructParticles().

3064  {
3065  assert(!blockref.isNull());
3066  const reco::PFBlock& block = *blockref;
3067 
3068  LogTrace("PFAlgo|processBlock") << "start of function PFAlgo::processBlock, block=" << block;
3069 
3071  LogTrace("PFAlgo|processBlock") << "elements.size()=" << elements.size();
3072  // make a copy of the link data, which will be edited.
3073  PFBlock::LinkData linkData = block.linkData();
3074 
3075  // keep track of the elements which are still active.
3076  vector<bool> active(elements.size(), true);
3077 
3078  // //PFElectrons:
3079  // usePFElectrons_ external configurable parameter to set the usage of pf electron
3080  std::vector<reco::PFCandidate> tempElectronCandidates;
3081  tempElectronCandidates.clear();
3082 
3083  // New EGamma Reconstruction 10/10/2013
3084  if (useEGammaFilters_) {
3085  egammaFilters(blockref, active, pfegamma);
3086  } // end if use EGammaFilters
3087 
3088  //Lock extra conversion tracks not used by Photon Algo
3089  if (usePFConversions_) {
3090  conversionAlgo(elements, active);
3091  }
3092 
3093  // In the following elementLoop() function, the primary goal is to deal with tracks that are:
3094  // - not associated to an HCAL cluster
3095  // - not identified as an electron.
3096  // Those tracks should be predominantly relatively low energy charged
3097  // hadrons which are not detected in the ECAL.
3098 
3099  // The secondary goal is to prepare for the next loops
3100  // - The ecal and hcal elements are sorted in separate vectors in `ElementIndices inds`
3101  // which will be used as a base for the corresponding loops.
3102  // - For tracks which are connected to more than one HCAL cluster,
3103  // the links between the track and the cluster are cut for all clusters
3104  // but the closest one.
3105  // - HF only blocks ( HFEM, HFHAD, HFEM+HFAD) are identified
3106 
3107  // obsolete comments?
3108  // loop1:
3109  // - sort ecal and hcal elements in separate vectors
3110  // - for tracks:
3111  // - lock closest ecal cluster
3112  // - cut link to farthest hcal cluster, if more than 1.
3113 
3114  vector<bool> deadArea(elements.size(), false);
3115 
3116  // vectors to store element indices to ho, hcal and ecal elements, will be filled by elementLoop()
3117  ElementIndices inds;
3118 
3119  elementLoop(block, linkData, elements, active, blockref, inds, deadArea);
3120 
3121  // Reconstruct pfCandidate from HF (either EM-only, Had-only or both)
3122  // For phase2, process also pfblocks containing HF clusters and linked tracks
3123  if (!(inds.hfEmIs.empty() && inds.hfHadIs.empty())) {
3124  createCandidatesHF(block, linkData, elements, active, blockref, inds);
3125  if (inds.hcalIs.empty() && inds.ecalIs.empty())
3126  return;
3127  LogDebug("PFAlgo::processBlock")
3128  << "Block contains HF clusters, and also contains ECAL or HCAL clusters. Continue.\n"
3129  << block;
3130  }
3131 
3132  createCandidatesHCAL(block, linkData, elements, active, blockref, inds, deadArea);
3133  // COLINFEB16: now dealing with the HCAL elements that are not linked to any track
3134  createCandidatesHCALUnlinked(block, linkData, elements, active, blockref, inds, deadArea);
3135  createCandidatesECAL(block, linkData, elements, active, blockref, inds, deadArea);
3136 
3137  LogTrace("PFAlgo|processBlock") << "end of function PFAlgo::processBlock";
3138 } // end processBlock
bool usePFConversions_
Definition: PFAlgo.h:280
void elementLoop(const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active, const reco::PFBlockRef &blockref, ElementIndices &inds, std::vector< bool > &deadArea)
Definition: PFAlgo.cc:997
std::map< unsigned int, Link > LinkData
Definition: PFBlock.h:39
void egammaFilters(const reco::PFBlockRef &blockref, std::vector< bool > &active, PFEGammaFilters const *pfegamma)
Definition: PFAlgo.cc:220
assert(be >=bs)
#define LogTrace(id)
void createCandidatesECAL(const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active, const reco::PFBlockRef &blockref, ElementIndices &inds, std::vector< bool > &deadArea)
Definition: PFAlgo.cc:3012
bool isNull() const
Checks for null.
Definition: Ref.h:235
void createCandidatesHCAL(const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active, const reco::PFBlockRef &blockref, ElementIndices &inds, std::vector< bool > &deadArea)
Definition: PFAlgo.cc:1702
void conversionAlgo(const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active)
Definition: PFAlgo.cc:370
void createCandidatesHF(const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active, const reco::PFBlockRef &blockref, ElementIndices &inds)
Definition: PFAlgo.cc:1259
bool useEGammaFilters_
Variables for NEW EGAMMA selection.
Definition: PFAlgo.h:265
void createCandidatesHCALUnlinked(const reco::PFBlock &block, reco::PFBlock::LinkData &linkData, const edm::OwnVector< reco::PFBlockElement > &elements, std::vector< bool > &active, const reco::PFBlockRef &blockref, ElementIndices &inds, std::vector< bool > &deadArea)
Definition: PFAlgo.cc:2813
#define LogDebug(id)
Block of elements.
Definition: PFBlock.h:26

◆ reconstructCluster()

unsigned PFAlgo::reconstructCluster ( const reco::PFCluster cluster,
double  particleEnergy,
bool  useDirection = false,
double  particleX = 0.,
double  particleY = 0.,
double  particleZ = 0. 
)
private

Reconstruct a neutral particle from a cluster. If chargedEnergy is specified, the neutral particle is created only if the cluster energy is significantly larger than the chargedEnergy. In this case, the energy of the neutral particle is cluster energy - chargedEnergy

Definition at line 3220 of file PFAlgo.cc.

References cms::cuda::assert(), ALCARECOTkAlJpsiMuMu_cff::charge, PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, DQMScaleToClient_cfi::factor, CustomPhysics_cfi::gamma, PFLayer::HCAL_BARREL1, PFLayer::HCAL_ENDCAP, PFLayer::HF_EM, PFLayer::HF_HAD, reco::PFCluster::layer(), LogTrace, EgHLTOffHistBins_cfi::mass, PbPb_ZMuSkimMuonDPG_cff::particleType, pfElectronTranslator_cfi::PFCandidate, pfCandidates_, reco::CaloCluster::position(), primaryVertex_, setHcalDepthInfo(), mathSSE::sqrt(), createJobs::tmp, useVertices_, reco::PFCandidate::X, reco::Vertex::x(), reco::Vertex::y(), and reco::Vertex::z().

Referenced by checkCleaning(), createCandidatesECAL(), createCandidatesHCAL(), createCandidatesHCALUnlinked(), createCandidatesHF(), and recoTracksNotHCAL().

3225  {
3226  LogTrace("PFAlgo|reconstructCluster") << "start of function PFAlgo::reconstructCluster, cluster=" << cluster
3227  << "particleEnergy=" << particleEnergy << "useDirection=" << useDirection
3228  << "particleX=" << particleX << "particleY=" << particleY
3229  << "particleZ=" << particleZ;
3230 
3232 
3233  // need to convert the ::math::XYZPoint data member of the PFCluster class=
3234  // to a displacement vector:
3235 
3236  // Transform particleX,Y,Z to a position at ECAL/HCAL entrance
3237  double factor = 1.;
3238  if (useDirection) {
3239  switch (cluster.layer()) {
3240  case PFLayer::ECAL_BARREL:
3241  case PFLayer::HCAL_BARREL1:
3242  factor = std::sqrt(cluster.position().Perp2() / (particleX * particleX + particleY * particleY));
3243  break;
3244  case PFLayer::ECAL_ENDCAP:
3245  case PFLayer::HCAL_ENDCAP:
3246  case PFLayer::HF_HAD:
3247  case PFLayer::HF_EM:
3248  factor = cluster.position().Z() / particleZ;
3249  break;
3250  default:
3251  assert(0);
3252  }
3253  }
3254  //MIKE First of all let's check if we have vertex.
3255  ::math::XYZPoint vertexPos;
3256  if (useVertices_)
3258  else
3259  vertexPos = ::math::XYZPoint(0.0, 0.0, 0.0);
3260 
3261  ::math::XYZVector clusterPos(cluster.position().X() - vertexPos.X(),
3262  cluster.position().Y() - vertexPos.Y(),
3263  cluster.position().Z() - vertexPos.Z());
3264  ::math::XYZVector particleDirection(
3265  particleX * factor - vertexPos.X(), particleY * factor - vertexPos.Y(), particleZ * factor - vertexPos.Z());
3266 
3267  //::math::XYZVector clusterPos( cluster.position().X(), cluster.position().Y(),cluster.position().Z() );
3268  //::math::XYZVector particleDirection ( particleX, particleY, particleZ );
3269 
3270  clusterPos = useDirection ? particleDirection.Unit() : clusterPos.Unit();
3271  clusterPos *= particleEnergy;
3272 
3273  // clusterPos is now a vector along the cluster direction,
3274  // with a magnitude equal to the cluster energy.
3275 
3276  double mass = 0;
3277  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> momentum(
3278  clusterPos.X(), clusterPos.Y(), clusterPos.Z(), mass);
3279  // mathcore is a piece of #$%
3281  // implicit constructor not allowed
3282  tmp = momentum;
3283 
3284  // Charge
3285  int charge = 0;
3286 
3287  // Type
3288  switch (cluster.layer()) {
3289  case PFLayer::ECAL_BARREL:
3290  case PFLayer::ECAL_ENDCAP:
3292  break;
3293  case PFLayer::HCAL_BARREL1:
3294  case PFLayer::HCAL_ENDCAP:
3295  particleType = PFCandidate::h0;
3296  break;
3297  case PFLayer::HF_HAD:
3298  particleType = PFCandidate::h_HF;
3299  break;
3300  case PFLayer::HF_EM:
3301  particleType = PFCandidate::egamma_HF;
3302  break;
3303  default:
3304  assert(0);
3305  }
3306 
3307  // The pf candidate
3308  LogTrace("PFAlgo|reconstructCluster") << "Creating PFCandidate charge=" << charge << ", type=" << particleType
3309  << ", pt=" << tmp.pt() << ", eta=" << tmp.eta() << ", phi=" << tmp.phi();
3311 
3312  // The position at ECAL entrance (well: watch out, it is not true
3313  // for HCAL clusters... to be fixed)
3314  pfCandidates_->back().setPositionAtECALEntrance(
3315  ::math::XYZPointF(cluster.position().X(), cluster.position().Y(), cluster.position().Z()));
3316 
3317  //Set the cnadidate Vertex
3318  pfCandidates_->back().setVertex(vertexPos);
3319 
3320  // depth info
3321  setHcalDepthInfo(pfCandidates_->back(), cluster);
3322 
3323  //*TODO* cluster time is not reliable at the moment, so only use track timing
3324 
3325  LogTrace("PFAlgo|reconstructCluster") << "** candidate: " << pfCandidates_->back();
3326 
3327  // returns index to the newly created PFCandidate
3328  return pfCandidates_->size() - 1;
3329 }
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:154
ParticleType
particle types
Definition: PFCandidate.h:44
double z() const
z coordinate
Definition: Vertex.h:133
std::unique_ptr< reco::PFCandidateCollection > pfCandidates_
Definition: PFAlgo.h:226
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float > > XYZPointF
point in space with cartesian internal representation
Definition: Point3D.h:10
PFLayer::Layer layer() const
cluster layer, see PFLayer.h in this directory
Definition: PFCluster.cc:56
assert(be >=bs)
#define LogTrace(id)
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
T sqrt(T t)
Definition: SSEVec.h:19
reco::Vertex primaryVertex_
Definition: PFAlgo.h:328
double x() const
x coordinate
Definition: Vertex.h:129
double y() const
y coordinate
Definition: Vertex.h:131
bool useVertices_
Definition: PFAlgo.h:329
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
math::XYZVector XYZPoint
void setHcalDepthInfo(reco::PFCandidate &cand, const reco::PFCluster &cluster) const
Definition: PFAlgo.cc:3331
tmp
align.sh
Definition: createJobs.py:716

◆ reconstructParticles()

void PFAlgo::reconstructParticles ( const reco::PFBlockHandle blockHandle,
PFEGammaFilters const *  pfegamma 
)

reconstruct particles

Definition at line 130 of file PFAlgo.cc.

References groupFilesInBlocks::block, gather_cfg::blocks, reco::PFBlockElement::ECAL, bookConverter::elements, relativeConstraints::empty, reco::PFBlockElement::HCAL, reco::PFBlockElement::HO, mps_fire::i, edm::HandleBase::isValid(), LogTrace, muonHandle_, trackingPlots::other, pfCandidates_, pfCleanedCandidates_, pfmu_, postCleaning(), postHFCleaning_, and processBlock().

Referenced by PFProducer::produce().

130  {
131  auto const& blocks = *blockHandle;
132 
133  // reset output collection
134  pfCandidates_->clear();
135 
136  LogTrace("PFAlgo|reconstructParticles")
137  << "start of function PFAlgo::reconstructParticles, blocks.size()=" << blocks.size();
138 
139  // sort elements in three lists:
140  std::list<reco::PFBlockRef> hcalBlockRefs;
141  std::list<reco::PFBlockRef> ecalBlockRefs;
142  std::list<reco::PFBlockRef> hoBlockRefs;
143  std::list<reco::PFBlockRef> otherBlockRefs;
144 
145  for (unsigned i = 0; i < blocks.size(); ++i) {
146  reco::PFBlockRef blockref = reco::PFBlockRef(blockHandle, i);
147 
148  const reco::PFBlock& block = *blockref;
150 
151  bool singleEcalOrHcal = false;
152  if (elements.size() == 1) {
154  ecalBlockRefs.push_back(blockref);
155  singleEcalOrHcal = true;
156  }
158  hcalBlockRefs.push_back(blockref);
159  singleEcalOrHcal = true;
160  }
161  if (elements[0].type() == reco::PFBlockElement::HO) {
162  // Single HO elements are likely to be noise. Not considered for now.
163  hoBlockRefs.push_back(blockref);
164  singleEcalOrHcal = true;
165  }
166  }
167 
168  if (!singleEcalOrHcal) {
169  otherBlockRefs.push_back(blockref);
170  }
171  } //loop blocks
172 
173  LogTrace("PFAlgo|reconstructParticles")
174  << "# Ecal blocks: " << ecalBlockRefs.size() << ", # Hcal blocks: " << hcalBlockRefs.size()
175  << ", # HO blocks: " << hoBlockRefs.size() << ", # Other blocks: " << otherBlockRefs.size();
176 
177  // loop on blocks that are not single ecal,
178  // and not single hcal.
179 
180  unsigned nblcks = 0;
181  for (auto const& other : otherBlockRefs) {
182  LogTrace("PFAlgo|reconstructParticles") << "processBlock, Block number " << nblcks++;
183  processBlock(other, hcalBlockRefs, ecalBlockRefs, pfegamma);
184  }
185 
186  std::list<reco::PFBlockRef> empty;
187 
188  unsigned hblcks = 0;
189  // process remaining single hcal blocks
190  for (auto const& hcal : hcalBlockRefs) {
191  LogTrace("PFAlgo|reconstructParticles") << "processBlock, HCAL block number " << hblcks++;
192  processBlock(hcal, empty, empty, pfegamma);
193  }
194 
195  unsigned eblcks = 0;
196  // process remaining single ecal blocks
197  for (auto const& ecal : ecalBlockRefs) {
198  LogTrace("PFAlgo|reconstructParticles") << "processBlock, ECAL block number " << eblcks++;
199  processBlock(ecal, empty, empty, pfegamma);
200  }
201 
202  // Post HF Cleaning
203  pfCleanedCandidates_.clear();
204  // Check if the post HF Cleaning was requested - if not, do nothing
205  if (postHFCleaning_) {
206  postCleaning();
207  }
208 
209  //Muon post cleaning
210  pfmu_->postClean(pfCandidates_.get());
211 
212  //Add Missing muons
213  if (muonHandle_.isValid())
214  pfmu_->addMissingMuons(muonHandle_, pfCandidates_.get());
215 
216  LogTrace("PFAlgo|reconstructParticles")
217  << "end of function PFAlgo::reconstructParticles, pfCandidates_->size()=" << pfCandidates_->size();
218 }
edm::Ref< PFBlockCollection > PFBlockRef
persistent reference to PFCluster objects
Definition: PFBlockFwd.h:19
std::unique_ptr< reco::PFCandidateCollection > pfCandidates_
Definition: PFAlgo.h:226
reco::PFCandidateCollection pfCleanedCandidates_
Definition: PFAlgo.h:228
#define LogTrace(id)
bool postHFCleaning_
Definition: PFAlgo.h:317
edm::Handle< reco::MuonCollection > muonHandle_
Definition: PFAlgo.h:331
void processBlock(const reco::PFBlockRef &blockref, std::list< reco::PFBlockRef > &hcalBlockRefs, std::list< reco::PFBlockRef > &ecalBlockRefs, PFEGammaFilters const *pfegamma)
Definition: PFAlgo.cc:3061
bool isValid() const
Definition: HandleBase.h:70
void postCleaning()
Definition: PFAlgo.cc:3499
std::unique_ptr< PFMuonAlgo > pfmu_
Definition: PFAlgo.h:262
Block of elements.
Definition: PFBlock.h:26

◆ reconstructTrack()

unsigned PFAlgo::reconstructTrack ( const reco::PFBlockElement elt,
bool  allowLoose = false 
)
private

Reconstruct a charged particle from a track Returns the index of the newly created candidate in pfCandidates_ Michalis added a flag here to treat muons inside jets

Definition at line 3141 of file PFAlgo.cc.

References ALCARECOTkAlJpsiMuMu_cff::charge, dptRel_DispVtx_, HCALHighEnergyHPDFilter_cfi::energy, reco::PFCandidate::h, isFromSecInt(), reco::isMuon(), edm::Ref< C, T, F >::isNonnull(), reco::PFBlockElement::isTimeValid(), LogTrace, reco::TrackBase::p(), PbPb_ZMuSkimMuonDPG_cff::particleType, pfElectronTranslator_cfi::PFCandidate, pfCandidates_, pfmu_, multPhiCorr_741_25nsDY_cfi::px, reco::TrackBase::px(), multPhiCorr_741_25nsDY_cfi::py, reco::TrackBase::py(), reco::TrackBase::pz(), mathSSE::sqrt(), reco::PFBlockElement::T_FROM_DISP, reco::PFCandidate::T_FROM_DISP, reco::PFBlockElement::T_TO_DISP, reco::PFCandidate::T_TO_DISP, reco::PFBlockElement::time(), reco::PFBlockElement::timeError(), and HLT_2022v12_cff::track.

Referenced by checkAndReconstructSecondaryInteraction(), createCandidatesHCAL(), createCandidatesHF(), and recoTracksNotHCAL().

3141  {
3142  const auto* eltTrack = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
3143 
3144  const reco::TrackRef& trackRef = eltTrack->trackRef();
3145  const reco::Track& track = *trackRef;
3146  const reco::MuonRef& muonRef = eltTrack->muonRef();
3147  int charge = track.charge() > 0 ? 1 : -1;
3148 
3149  // Assume this particle is a charged Hadron
3150  double px = track.px();
3151  double py = track.py();
3152  double pz = track.pz();
3153  double energy = sqrt(track.p() * track.p() + 0.13957 * 0.13957);
3154 
3155  LogTrace("PFAlgo|reconstructTrack") << "Reconstructing PF candidate from track of pT = " << track.pt()
3156  << " eta = " << track.eta() << " phi = " << track.phi() << " px = " << px
3157  << " py = " << py << " pz = " << pz << " energy = " << energy;
3158 
3159  // Create a PF Candidate
3160  ::math::XYZTLorentzVector momentum(px, py, pz, energy);
3162 
3163  // Add it to the stack
3164  LogTrace("PFAlgo|reconstructTrack") << "Creating PFCandidate charge=" << charge << ", type=" << particleType
3165  << ", pt=" << momentum.pt() << ", eta=" << momentum.eta()
3166  << ", phi=" << momentum.phi();
3167  pfCandidates_->push_back(PFCandidate(charge, momentum, particleType));
3168  //Set vertex and stuff like this
3169  pfCandidates_->back().setVertex(trackRef->vertex());
3170  pfCandidates_->back().setTrackRef(trackRef);
3171  pfCandidates_->back().setPositionAtECALEntrance(eltTrack->positionAtECALEntrance());
3172  if (muonRef.isNonnull())
3173  pfCandidates_->back().setMuonRef(muonRef);
3174 
3175  //Set time
3176  if (elt.isTimeValid())
3177  pfCandidates_->back().setTime(elt.time(), elt.timeError());
3178 
3179  //OK Now try to reconstruct the particle as a muon
3180  bool isMuon = pfmu_->reconstructMuon(pfCandidates_->back(), muonRef, allowLoose);
3181  bool isFromDisp = isFromSecInt(elt, "secondary");
3182 
3183  if ((!isMuon) && isFromDisp) {
3184  double dpt = trackRef->ptError();
3185  double dptRel = dpt / trackRef->pt() * 100;
3186  //If the track is ill measured it is better to not refit it, since the track information probably would not be used.
3187  //In the PFAlgo we use the trackref information. If the track error is too big the refitted information might be very different
3188  // from the not refitted one.
3189  if (dptRel < dptRel_DispVtx_) {
3190  LogTrace("PFAlgo|reconstructTrack")
3191  << "Not refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy;
3192  //reco::TrackRef trackRef = eltTrack->trackRef();
3194  eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef();
3195  reco::Track trackRefit = vRef->refittedTrack(trackRef);
3196  //change the momentum with the refitted track
3197  ::math::XYZTLorentzVector momentum(
3198  trackRefit.px(), trackRefit.py(), trackRefit.pz(), sqrt(trackRefit.p() * trackRefit.p() + 0.13957 * 0.13957));
3199  LogTrace("PFAlgo|reconstructTrack")
3200  << "Refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy;
3201  }
3202  pfCandidates_->back().setFlag(reco::PFCandidate::T_FROM_DISP, true);
3203  pfCandidates_->back().setDisplacedVertexRef(
3204  eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef(),
3206  }
3207 
3208  // do not label as primary a track which would be recognised as a muon. A muon cannot produce NI. It is with high probability a fake
3209  if (isFromSecInt(elt, "primary") && !isMuon) {
3210  pfCandidates_->back().setFlag(reco::PFCandidate::T_TO_DISP, true);
3211  pfCandidates_->back().setDisplacedVertexRef(
3212  eltTrack->displacedVertexRef(reco::PFBlockElement::T_TO_DISP)->displacedVertexRef(),
3214  }
3215 
3216  // returns index to the newly created PFCandidate
3217  return pfCandidates_->size() - 1;
3218 }
ParticleType
particle types
Definition: PFCandidate.h:44
bool isMuon(const Candidate &part)
Definition: pdgIdUtils.h:9
std::unique_ptr< reco::PFCandidateCollection > pfCandidates_
Definition: PFAlgo.h:226
double px() const
x coordinate of momentum vector
Definition: TrackBase.h:640
float timeError() const
double p() const
momentum vector magnitude
Definition: TrackBase.h:631
double py() const
y coordinate of momentum vector
Definition: TrackBase.h:643
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
bool isFromSecInt(const reco::PFBlockElement &eTrack, std::string order) const
Definition: PFAlgo.cc:3470
#define LogTrace(id)
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
T sqrt(T t)
Definition: SSEVec.h:19
bool isTimeValid() const
do we have a valid time information
double dptRel_DispVtx_
Definition: PFAlgo.h:285
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:646
std::unique_ptr< PFMuonAlgo > pfmu_
Definition: PFAlgo.h:262

◆ recoTracksNotHCAL()

bool PFAlgo::recoTracksNotHCAL ( const reco::PFBlock block,
reco::PFBlock::LinkData linkData,
const edm::OwnVector< reco::PFBlockElement > &  elements,
const reco::PFBlockRef blockref,
std::vector< bool > &  active,
bool  goodTrackDeadHcal,
bool  hasDeadHcal,
unsigned int  iTrack,
std::multimap< double, unsigned > &  ecalElems,
reco::TrackRef trackRef 
)
private

Definition at line 399 of file PFAlgo.cc.

References funct::abs(), cms::cuda::assert(), associatePSClusters(), groupFilesInBlocks::block, calibration_, dptRel_DispVtx_, ECAL, reco::PFBlockElement::ECAL, bookConverter::elements, reco::PFCandidate::elementsInBlocks(), PFEnergyCalibration::energyEmHad(), PVValHelper::eta, HLT_2022v12_cff::fraction, isFromSecInt(), PFMuonAlgo::isMuon(), edm::Ref< C, T, F >::isNull(), phase2tkutil::isPrimary(), reco::PFBlock::LINKTEST_ALL, LogTrace, SiStripPI::max, SiStripPI::min, reco::PFCandidate::mu, muonECAL_, neutralHadronEnergyResolution(), nHits, nSigmaECAL_, nSigmaTRACK_, pfCandidates_, reco::Vertex::position(), primaryVertex_, reco::PFBlockElement::PS1, reco::PFBlockElement::PS2, ptError_, reconstructCluster(), reconstructTrack(), rejectTracks_Step45_, pfClustersFromCombinedCaloHCal_cfi::resol, reco::PFBlockElementTrack::setPositionAtECALEntrance(), optionsL1T::skip, PFTrackAlgoTools::step45(), reco::PFBlockElement::TRACK, and reco::btau::trackMomentum.

Referenced by elementLoop().

408  {
409  LogTrace("PFAlgo|recoTracksNotHCAL") << "start of function PFAlgo::recoTracksNotHCAL(), now dealing with tracks "
410  "linked to no HCAL clusters. Was HCal active? "
411  << (!hasDeadHcal);
412  // vector<unsigned> elementIndices;
413  // elementIndices.push_back(iTrack);
414 
415  // The track momentum
416  double trackMomentum = elements[iTrack].trackRef()->p();
417  LogTrace("PFAlgo|recoTracksNotHCAL") << elements[iTrack];
418 
419  // Is it a "tight" muon ? In that case assume no
420  //Track momentum corresponds to the calorimeter
421  //energy for now
422  bool thisIsAMuon = PFMuonAlgo::isMuon(elements[iTrack]);
423  if (thisIsAMuon)
424  trackMomentum = 0.;
425 
426  // If it is not a muon check if Is it a fake track ?
427  //Michalis: I wonder if we should convert this to dpt/pt?
428  if (!thisIsAMuon && elements[iTrack].trackRef()->ptError() > ptError_) {
429  double deficit = trackMomentum;
430  double resol = neutralHadronEnergyResolution(trackMomentum, elements[iTrack].trackRef()->eta());
431  resol *= trackMomentum;
432 
433  if (!ecalElems.empty()) {
434  unsigned thisEcal = ecalElems.begin()->second;
435  reco::PFClusterRef clusterRef = elements[thisEcal].clusterRef();
436  bool useCluster = true;
437  if (hasDeadHcal) {
438  std::multimap<double, unsigned> sortedTracks;
439  block.associatedElements(
440  thisEcal, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
441  useCluster = (sortedTracks.begin()->second == iTrack);
442  }
443  if (useCluster) {
444  deficit -= clusterRef->energy();
445  resol = neutralHadronEnergyResolution(trackMomentum, clusterRef->positionREP().Eta());
446  resol *= trackMomentum;
447  }
448  }
449 
450  bool isPrimary = isFromSecInt(elements[iTrack], "primary");
451 
452  if (deficit > nSigmaTRACK_ * resol && !isPrimary && !goodTrackDeadHcal) {
453  active[iTrack] = false;
454  LogTrace("PFAlgo|recoTracksNotHCAL")
455  << elements[iTrack] << std::endl
456  << " deficit " << deficit << " > " << nSigmaTRACK_ << " x " << resol << " track pt " << trackRef->pt()
457  << " +- " << trackRef->ptError() << " layers valid " << trackRef->hitPattern().trackerLayersWithMeasurement()
458  << ", lost " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::TRACK_HITS)
459  << ", lost outer " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_OUTER_HITS)
460  << ", lost inner " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_INNER_HITS)
461  << "(valid fraction " << trackRef->validFraction() << ")"
462  << " chi2/ndf " << trackRef->normalizedChi2() << " |dxy| "
463  << std::abs(trackRef->dxy(primaryVertex_.position())) << " +- " << trackRef->dxyError() << " |dz| "
464  << std::abs(trackRef->dz(primaryVertex_.position())) << " +- " << trackRef->dzError()
465  << "is probably a fake (1) --> lock the track";
466  return true;
467  }
468  } // !thisIsaMuon
469 
470  // Create a track candidate
471  // unsigned tmpi = reconstructTrack( elements[iTrack] );
472  //active[iTrack] = false;
473  std::vector<unsigned> tmpi;
474  std::vector<unsigned> kTrack;
475 
476  // Some cleaning : secondary tracks without calo's and large momentum must be fake
477  double dpt = trackRef->ptError();
478 
479  if (rejectTracks_Step45_ && ecalElems.empty() && trackMomentum > 30. && dpt > 0.5 &&
480  (PFTrackAlgoTools::step45(trackRef->algo()) && !goodTrackDeadHcal)) {
481  double dptRel = dpt / trackRef->pt() * 100;
482  bool isPrimaryOrSecondary = isFromSecInt(elements[iTrack], "all");
483 
484  if (isPrimaryOrSecondary && dptRel < dptRel_DispVtx_) {
485  return true;
486  };
487  unsigned nHits = elements[iTrack].trackRef()->hitPattern().trackerLayersWithMeasurement();
488  unsigned int NLostHit = trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::TRACK_HITS);
489 
490  LogTrace("PFAlgo|recoTracksNotHCAL") << "A track (algo = " << trackRef->algo() << ") with momentum "
491  << trackMomentum << " / " << elements[iTrack].trackRef()->pt() << " +/- "
492  << dpt << " / " << elements[iTrack].trackRef()->eta()
493  << " without any link to ECAL/HCAL and with " << nHits << " (" << NLostHit
494  << ") hits (lost hits) has been cleaned";
495 
496  active[iTrack] = false;
497  return true;
498  } //rejectTracks_Step45_ && ...
499 
500  tmpi.push_back(reconstructTrack(elements[iTrack]));
501 
502  kTrack.push_back(iTrack);
503  active[iTrack] = false;
504 
505  // No ECAL cluster either ... continue...
506  if (ecalElems.empty()) {
507  (*pfCandidates_)[tmpi[0]].setEcalEnergy(0., 0.);
508  (*pfCandidates_)[tmpi[0]].setHcalEnergy(0., 0.);
509  (*pfCandidates_)[tmpi[0]].setHoEnergy(0., 0.);
510  (*pfCandidates_)[tmpi[0]].setPs1Energy(0);
511  (*pfCandidates_)[tmpi[0]].setPs2Energy(0);
512  (*pfCandidates_)[tmpi[0]].addElementInBlock(blockref, kTrack[0]);
513  return true;
514  }
515 
516  // Look for closest ECAL cluster
517  const unsigned int thisEcal = ecalElems.begin()->second;
518  reco::PFClusterRef clusterRef = elements[thisEcal].clusterRef();
519  LogTrace("PFAlgo|recoTracksNotHCAL") << " is associated to " << elements[thisEcal];
520 
521  // Set ECAL energy for muons
522  if (thisIsAMuon) {
523  (*pfCandidates_)[tmpi[0]].setEcalEnergy(clusterRef->energy(), std::min(clusterRef->energy(), muonECAL_[0]));
524  (*pfCandidates_)[tmpi[0]].setHcalEnergy(0., 0.);
525  (*pfCandidates_)[tmpi[0]].setHoEnergy(0., 0.);
526  (*pfCandidates_)[tmpi[0]].setPs1Energy(0);
527  (*pfCandidates_)[tmpi[0]].setPs2Energy(0);
528  (*pfCandidates_)[tmpi[0]].addElementInBlock(blockref, kTrack[0]);
529  }
530 
531  double slopeEcal = 1.;
532  bool connectedToEcal = false;
533  unsigned iEcal = 99999;
534  double calibEcal = 0.;
535  double calibHcal = 0.;
536  double totalEcal = thisIsAMuon ? -muonECAL_[0] : 0.;
537 
538  // Consider charged particles closest to the same ECAL cluster
539  std::multimap<double, unsigned> sortedTracks;
540  block.associatedElements(thisEcal, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
541  LogTrace("PFAlgo|recoTracksNotHCAL") << "the closest ECAL cluster, id " << thisEcal << ", has " << sortedTracks.size()
542  << " associated tracks, now processing them. ";
543 
544  if (hasDeadHcal && !sortedTracks.empty()) {
545  // GP: only allow the ecal cluster closest to this track to be considered
546  LogTrace("PFAlgo|recoTracksNotHCAL") << " the closest track to ECAL " << thisEcal << " is "
547  << sortedTracks.begin()->second << " (distance " << sortedTracks.begin()->first
548  << ")";
549  if (sortedTracks.begin()->second != iTrack) {
550  LogTrace("PFAlgo|recoTracksNotHCAL")
551  << " the closest track to ECAL " << thisEcal << " is " << sortedTracks.begin()->second
552  << " which is not the one being processed. Will skip ECAL linking for this track";
553  (*pfCandidates_)[tmpi[0]].setEcalEnergy(0., 0.);
554  (*pfCandidates_)[tmpi[0]].setHcalEnergy(0., 0.);
555  (*pfCandidates_)[tmpi[0]].setHoEnergy(0., 0.);
556  (*pfCandidates_)[tmpi[0]].setPs1Energy(0);
557  (*pfCandidates_)[tmpi[0]].setPs2Energy(0);
558  (*pfCandidates_)[tmpi[0]].addElementInBlock(blockref, kTrack[0]);
559  return true;
560  } else {
561  LogTrace("PFAlgo|recoTracksNotHCAL")
562  << " the closest track to ECAL " << thisEcal << " is " << sortedTracks.begin()->second
563  << " which is the one being processed. Will skip ECAL linking for all other track";
564  sortedTracks.clear();
565  }
566  }
567 
568  for (auto const& trk : sortedTracks) {
569  unsigned jTrack = trk.second;
570 
571  // Don't consider already used tracks
572  if (!active[jTrack])
573  continue;
574 
575  // The loop is on the other tracks !
576  if (jTrack == iTrack)
577  continue;
578 
579  // Check if the ECAL cluster closest to this track is
580  // indeed the current ECAL cluster. Only tracks not
581  // linked to an HCAL cluster are considered here
582  // (GP: this is because if there's a jTrack linked
583  // to the same Ecal cluster and that also has an Hcal link
584  // we would have also linked iTrack to that Hcal above)
585  std::multimap<double, unsigned> sortedECAL;
586  block.associatedElements(jTrack, linkData, sortedECAL, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
587  if (sortedECAL.begin()->second != thisEcal)
588  continue;
589 
590  // Check if this track is a muon
591  bool thatIsAMuon = PFMuonAlgo::isMuon(elements[jTrack]);
592  LogTrace("PFAlgo|recoTracksNotHCAL") << " found track " << jTrack << (thatIsAMuon ? " (muon) " : " (non-muon)")
593  << ", with distance = " << sortedECAL.begin()->first;
594 
595  // Check if this track is not a fake
596  bool rejectFake = false;
597  reco::TrackRef trackRef = elements[jTrack].trackRef();
598  if (!thatIsAMuon && trackRef->ptError() > ptError_) {
599  // GP: FIXME this selection should be adapted depending on hasDeadHcal
600  // but we don't know if jTrack is linked to a dead Hcal or not
601  double deficit = trackMomentum + trackRef->p() - clusterRef->energy();
602  double resol =
603  nSigmaTRACK_ * neutralHadronEnergyResolution(trackMomentum + trackRef->p(), clusterRef->positionREP().Eta());
604  resol *= (trackMomentum + trackRef->p());
605  if (deficit > nSigmaTRACK_ * resol && !goodTrackDeadHcal) {
606  rejectFake = true;
607  kTrack.push_back(jTrack);
608  active[jTrack] = false;
609  LogTrace("PFAlgo|recoTracksNotHCAL")
610  << elements[jTrack] << std::endl
611  << "is probably a fake (2) --> lock the track"
612  << "(trackMomentum = " << trackMomentum << ", clusterEnergy = " << clusterRef->energy()
613  << ", deficit = " << deficit << " > " << nSigmaTRACK_ << " x " << resol
614  << " assuming neutralHadronEnergyResolution "
615  << neutralHadronEnergyResolution(trackMomentum + trackRef->p(), clusterRef->positionREP().Eta()) << ") ";
616  }
617  }
618  if (rejectFake)
619  continue;
620 
621  // Otherwise, add this track momentum to the total track momentum
622  /* */
623  // reco::TrackRef trackRef = elements[jTrack].trackRef();
624  if (!thatIsAMuon) {
625  LogTrace("PFAlgo|recoTracksNotHCAL") << "Track momentum increased from " << trackMomentum << " GeV ";
626  trackMomentum += trackRef->p();
627  LogTrace("PFAlgo|recoTracksNotHCAL") << "to " << trackMomentum << " GeV.";
628  LogTrace("PFAlgo|recoTracksNotHCAL") << "with " << elements[jTrack];
629  } else {
630  totalEcal -= muonECAL_[0];
631  totalEcal = std::max(totalEcal, 0.);
632  }
633 
634  // And create a charged particle candidate !
635 
636  tmpi.push_back(reconstructTrack(elements[jTrack]));
637 
638  kTrack.push_back(jTrack);
639  active[jTrack] = false;
640 
641  if (thatIsAMuon) {
642  (*pfCandidates_)[tmpi.back()].setEcalEnergy(clusterRef->energy(), std::min(clusterRef->energy(), muonECAL_[0]));
643  (*pfCandidates_)[tmpi.back()].setHcalEnergy(0., 0.);
644  (*pfCandidates_)[tmpi.back()].setHoEnergy(0., 0.);
645  (*pfCandidates_)[tmpi.back()].setPs1Energy(0);
646  (*pfCandidates_)[tmpi.back()].setPs2Energy(0);
647  (*pfCandidates_)[tmpi.back()].addElementInBlock(blockref, kTrack.back());
648  }
649  }
650 
651  LogTrace("PFAlgo|recoTracksNotHCAL") << "Loop over all associated ECAL clusters";
652  // Loop over all ECAL linked clusters ordered by increasing distance.
653  for (auto const& ecal : ecalElems) {
654  const unsigned index = ecal.second;
655  const PFBlockElement::Type type = elements[index].type();
657  LogTrace("PFAlgo|recoTracksNotHCAL") << elements[index];
658 
659  // Just skip clusters already taken
660  if (!active[index]) {
661  LogTrace("PFAlgo|recoTracksNotHCAL") << "is not active - ignore ";
662  continue;
663  }
664 
665  // Just skip this cluster if it's closer to another track
666  block.associatedElements(index, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
667 
668  const bool skip = std::none_of(
669  kTrack.begin(), kTrack.end(), [&](unsigned iTrack) { return sortedTracks.begin()->second == iTrack; });
670 
671  if (skip) {
672  LogTrace("PFAlgo|recoTracksNotHCAL") << "is closer to another track - ignore ";
673  continue;
674  }
675 
676  // The corresponding PFCluster and energy
677  const reco::PFClusterRef clusterRef = elements[index].clusterRef();
678  assert(!clusterRef.isNull());
679 
680  LogTrace("PFAlgo|recoTracksNotHCAL") << "Ecal cluster with raw energy = " << clusterRef->energy()
681  << " linked with distance = " << ecal.first;
682 
683  // Check the presence of preshower clusters in the vicinity
684  // Preshower cluster closer to another ECAL cluster are ignored.
685  //JOSH: This should use the association map already produced by the cluster corrector for consistency,
686  //but should be ok for now
687  vector<double> ps1Ene{0.};
688  associatePSClusters(index, reco::PFBlockElement::PS1, block, elements, linkData, active, ps1Ene);
689  vector<double> ps2Ene{0.};
690  associatePSClusters(index, reco::PFBlockElement::PS2, block, elements, linkData, active, ps2Ene);
691 
692  // KH: use raw ECAL energy for PF hadron calibration. use calibrated ECAL energy when adding PF photons
693  const double ecalEnergy = clusterRef->energy();
694  const double ecalEnergyCalibrated = clusterRef->correctedEnergy(); // calibrated based on the egamma hypothesis
695  LogTrace("PFAlgo|recoTracksNotHCAL") << "Corrected ECAL(+PS) energy = " << ecalEnergy;
696 
697  // Since the electrons were found beforehand, this track must be a hadron. Calibrate
698  // the energy under the hadron hypothesis.
699  totalEcal += ecalEnergy;
700  const double previousCalibEcal = calibEcal;
701  const double previousSlopeEcal = slopeEcal;
702  calibEcal = std::max(totalEcal, 0.);
703  calibHcal = 0.;
705  trackMomentum, calibEcal, calibHcal, clusterRef->positionREP().Eta(), clusterRef->positionREP().Phi());
706  if (totalEcal > 0.)
707  slopeEcal = calibEcal / totalEcal;
708 
709  LogTrace("PFAlgo|recoTracksNotHCAL") << "The total calibrated energy so far amounts to = " << calibEcal
710  << " (slope = " << slopeEcal << ")";
711 
712  // Stop the loop when adding more ECAL clusters ruins the compatibility
713  if (connectedToEcal && calibEcal - trackMomentum >= 0.) {
714  // if ( connectedToEcal && calibEcal - trackMomentum >=
715  // nSigmaECAL_*neutralHadronEnergyResolution(trackMomentum,clusterRef->positionREP().Eta()) ) {
716  calibEcal = previousCalibEcal;
717  slopeEcal = previousSlopeEcal;
718  totalEcal = calibEcal / slopeEcal;
719 
720  // Turn this last cluster in a photon
721  // (The PS clusters are already locked in "associatePSClusters")
722  active[index] = false;
723 
724  // Find the associated tracks
725  std::multimap<double, unsigned> assTracks;
726  block.associatedElements(index, linkData, assTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
727 
728  auto& ecalCand = (*pfCandidates_)[reconstructCluster(
729  *clusterRef, ecalEnergyCalibrated)]; // KH: use the PF ECAL cluster calibrated energy
730  ecalCand.setEcalEnergy(clusterRef->energy(), ecalEnergyCalibrated);
731  ecalCand.setHcalEnergy(0., 0.);
732  ecalCand.setHoEnergy(0., 0.);
733  ecalCand.setPs1Energy(ps1Ene[0]);
734  ecalCand.setPs2Energy(ps2Ene[0]);
735  ecalCand.addElementInBlock(blockref, index);
736  // Check that there is at least one track
737  if (!assTracks.empty()) {
738  ecalCand.addElementInBlock(blockref, assTracks.begin()->second);
739 
740  // Assign the position of the track at the ECAL entrance
741  const ::math::XYZPointF& chargedPosition =
742  dynamic_cast<const reco::PFBlockElementTrack*>(&elements[assTracks.begin()->second])
743  ->positionAtECALEntrance();
744  ecalCand.setPositionAtECALEntrance(chargedPosition);
745  }
746  break;
747  }
748 
749  // Lock used clusters.
750  connectedToEcal = true;
751  iEcal = index;
752  active[index] = false;
753  for (unsigned ic : tmpi)
754  (*pfCandidates_)[ic].addElementInBlock(blockref, iEcal);
755 
756  } // Loop ecal elements
757 
758  bool bNeutralProduced = false;
759 
760  // Create a photon if the ecal energy is too large
761  if (connectedToEcal) {
762  reco::PFClusterRef pivotalRef = elements[iEcal].clusterRef();
763 
764  double neutralEnergy = calibEcal - trackMomentum;
765 
766  /*
767  // Check if there are other tracks linked to that ECAL
768  for(IE ie = sortedTracks.begin(); ie != sortedTracks.end() && neutralEnergy > 0; ++ie ) {
769  unsigned jTrack = ie->second;
770 
771  // The loop is on the other tracks !
772  if ( jTrack == iTrack ) continue;
773 
774  // Check if the ECAL closest to this track is the current ECAL
775  // Otherwise ignore this track in the neutral energy determination
776  std::multimap<double, unsigned> sortedECAL;
777  block.associatedElements( jTrack, linkData,
778  sortedECAL,
779  reco::PFBlockElement::ECAL,
780  reco::PFBlock::LINKTEST_ALL );
781  if ( sortedECAL.begin()->second != iEcal ) continue;
782 
783  // Check if this track is also linked to an HCAL
784  // (in which case it goes to the next loop and is ignored here)
785  std::multimap<double, unsigned> sortedHCAL;
786  block.associatedElements( jTrack, linkData,
787  sortedHCAL,
788  reco::PFBlockElement::HCAL,
789  reco::PFBlock::LINKTEST_ALL );
790  if ( sortedHCAL.size() ) continue;
791 
792  // Otherwise, subtract this track momentum from the ECAL energy
793  reco::TrackRef trackRef = elements[jTrack].trackRef();
794  neutralEnergy -= trackRef->p();
795 
796  } // End other tracks
797  */
798 
799  // Add a photon if the energy excess is large enough
800  double resol = neutralHadronEnergyResolution(trackMomentum, pivotalRef->positionREP().Eta());
801  resol *= trackMomentum;
802  if (neutralEnergy > std::max(0.5, nSigmaECAL_ * resol)) {
803  neutralEnergy /= slopeEcal;
804  unsigned tmpj = reconstructCluster(*pivotalRef, neutralEnergy);
805  (*pfCandidates_)[tmpj].setEcalEnergy(pivotalRef->energy(), neutralEnergy);
806  (*pfCandidates_)[tmpj].setHcalEnergy(0., 0.);
807  (*pfCandidates_)[tmpj].setHoEnergy(0., 0.);
808  (*pfCandidates_)[tmpj].setPs1Energy(0.);
809  (*pfCandidates_)[tmpj].setPs2Energy(0.);
810  (*pfCandidates_)[tmpj].addElementInBlock(blockref, iEcal);
811  bNeutralProduced = true;
812  for (unsigned ic = 0; ic < kTrack.size(); ++ic)
813  (*pfCandidates_)[tmpj].addElementInBlock(blockref, kTrack[ic]);
814  } // End neutral energy
815 
816  // Set elements in blocks and ECAL energies to all tracks
817  for (unsigned ic = 0; ic < tmpi.size(); ++ic) {
818  // Skip muons
819  if ((*pfCandidates_)[tmpi[ic]].particleId() == reco::PFCandidate::mu)
820  continue;
821 
822  double fraction = trackMomentum > 0 ? (*pfCandidates_)[tmpi[ic]].trackRef()->p() / trackMomentum : 0;
823  double ecalCal = bNeutralProduced ? (calibEcal - neutralEnergy * slopeEcal) * fraction : calibEcal * fraction;
824  double ecalRaw = totalEcal * fraction;
825 
826  LogTrace("PFAlgo|recoTracksNotHCAL")
827  << "The fraction after photon supression is " << fraction << " calibrated ecal = " << ecalCal;
828 
829  (*pfCandidates_)[tmpi[ic]].setEcalEnergy(ecalRaw, ecalCal);
830  (*pfCandidates_)[tmpi[ic]].setHcalEnergy(0., 0.);
831  (*pfCandidates_)[tmpi[ic]].setHoEnergy(0., 0.);
832  (*pfCandidates_)[tmpi[ic]].setPs1Energy(0);
833  (*pfCandidates_)[tmpi[ic]].setPs2Energy(0);
834  (*pfCandidates_)[tmpi[ic]].addElementInBlock(blockref, kTrack[ic]);
835  }
836 
837  } // End connected ECAL
838 
839  // Fill the element_in_block for tracks that are eventually linked to no ECAL clusters at all.
840  for (unsigned ic = 0; ic < tmpi.size(); ++ic) {
841  const PFCandidate& pfc = (*pfCandidates_)[tmpi[ic]];
842  const PFCandidate::ElementsInBlocks& eleInBlocks = pfc.elementsInBlocks();
843  if (eleInBlocks.empty()) {
844  LogTrace("PFAlgo|recoTracksNotHCAL") << "Single track / Fill element in block! ";
845  (*pfCandidates_)[tmpi[ic]].addElementInBlock(blockref, kTrack[ic]);
846  }
847  }
848  LogTrace("PFAlgo|recoTracksNotHCAL") << "end of function PFAlgo::recoTracksNotHCAL";
849  return false;
850 }
double ptError_
Definition: PFAlgo.h:297
bool isPrimary(const SimTrack &simTrk, const PSimHit *simHit)
static bool isMuon(const reco::PFBlockElement &elt)
Definition: PFMuonAlgo.cc:48
std::unique_ptr< reco::PFCandidateCollection > pfCandidates_
Definition: PFAlgo.h:226
unsigned reconstructCluster(const reco::PFCluster &cluster, double particleEnergy, bool useDirection=false, double particleX=0., double particleY=0., double particleZ=0.)
Definition: PFAlgo.cc:3220
const Point & position() const
position
Definition: Vertex.h:127
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float > > XYZPointF
point in space with cartesian internal representation
Definition: Point3D.h:10
bool isFromSecInt(const reco::PFBlockElement &eTrack, std::string order) const
Definition: PFAlgo.cc:3470
assert(be >=bs)
bool rejectTracks_Step45_
Definition: PFAlgo.h:277
#define LogTrace(id)
std::vector< ElementInBlock > ElementsInBlocks
Definition: PFCandidate.h:404
const ElementsInBlocks & elementsInBlocks() const
Definition: PFCandidate.cc:658
void associatePSClusters(unsigned iEcal, reco::PFBlockElement::Type psElementType, const reco::PFBlock &block, const edm::OwnVector< reco::PFBlockElement > &elements, const reco::PFBlock::LinkData &linkData, std::vector< bool > &active, std::vector< double > &psEne)
Associate PS clusters to a given ECAL cluster, and return their energy.
Definition: PFAlgo.cc:3425
bool step45(const reco::TrackBase::TrackAlgorithm &)
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
reco::Vertex primaryVertex_
Definition: PFAlgo.h:328
std::vector< double > muonECAL_
Definition: PFAlgo.h:294
void setPositionAtECALEntrance(float x, float y, float z)
set position at ECAL entrance
PFEnergyCalibration & calibration_
Definition: PFAlgo.h:257
void energyEmHad(double t, double &e, double &h, double eta, double phi) const
bool isNull() const
Checks for null.
Definition: Ref.h:235
unsigned reconstructTrack(const reco::PFBlockElement &elt, bool allowLoose=false)
Definition: PFAlgo.cc:3141
double dptRel_DispVtx_
Definition: PFAlgo.h:285
double neutralHadronEnergyResolution(double clusterEnergy, double clusterEta) const
todo: use PFClusterTools for this
Definition: PFAlgo.cc:3361
caConstants::TupleMultiplicity const CAHitNtupletGeneratorKernelsGPU::HitToTuple const cms::cuda::AtomicPairCounter GPUCACell const *__restrict__ uint32_t const *__restrict__ gpuPixelDoublets::CellNeighborsVector const gpuPixelDoublets::CellTracksVector const GPUCACell::OuterHitOfCell const int32_t nHits
const double nSigmaECAL_
number of sigma to judge energy excess in ECAL
Definition: PFAlgo.h:245
Particle reconstructed by the particle flow algorithm.
Definition: PFCandidate.h:41
double nSigmaTRACK_
Definition: PFAlgo.h:296

◆ relinkTrackToHcal()

void PFAlgo::relinkTrackToHcal ( const reco::PFBlock block,
std::multimap< double, unsigned > &  ecalElems,
std::multimap< double, unsigned > &  hcalElems,
const std::vector< bool > &  active,
reco::PFBlock::LinkData linkData,
unsigned int  iTrack 
)
private

Definition at line 941 of file PFAlgo.cc.

References groupFilesInBlocks::block, reco::PFBlockElement::ECAL, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL, LogTrace, and reco::PFBlockElement::TRACK.

Referenced by elementLoop().

946  {
947  unsigned ntt = 1;
948  unsigned index = ecalElems.begin()->second;
949  std::multimap<double, unsigned> sortedTracks;
950  block.associatedElements(index, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
951  LogTrace("PFAlgo|elementLoop") << "The closest ECAL cluster is linked to " << sortedTracks.size()
952  << " tracks, with distance = " << ecalElems.begin()->first;
953 
954  LogTrace("PFAlgo|elementLoop") << "Looping over sortedTracks";
955  // Loop over all tracks
956  for (auto const& trk : sortedTracks) {
957  unsigned jTrack = trk.second;
958  LogTrace("PFAlgo|elementLoop") << "jTrack=" << jTrack;
959  // Track must be active
960  if (!active[jTrack])
961  continue;
962  LogTrace("PFAlgo|elementLoop") << "active[jTrack]=" << active[jTrack];
963 
964  // The loop is on the other tracks !
965  if (jTrack == iTrack)
966  continue;
967  LogTrace("PFAlgo|elementLoop") << "skipping jTrack=" << jTrack << " for same iTrack";
968 
969  // Check if the ECAL closest to this track is the current ECAL
970  // Otherwise ignore this track in the neutral energy determination
971  std::multimap<double, unsigned> sortedECAL;
972  block.associatedElements(jTrack, linkData, sortedECAL, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
973  if (sortedECAL.begin()->second != index)
974  continue;
975  LogTrace("PFAlgo|elementLoop") << " track " << jTrack << " with closest ECAL identical ";
976 
977  // Check if this track is also linked to an HCAL
978  std::multimap<double, unsigned> sortedHCAL;
979  block.associatedElements(jTrack, linkData, sortedHCAL, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
980  if (sortedHCAL.empty())
981  continue;
982  LogTrace("PFAlgo|elementLoop") << " and with an HCAL cluster " << sortedHCAL.begin()->second;
983  ntt++;
984 
985  // In that case establish a link with the first track
986  block.setLink(iTrack, sortedHCAL.begin()->second, sortedECAL.begin()->first, linkData, PFBlock::LINKTEST_RECHIT);
987 
988  } // End other tracks
989 
990  // Redefine HCAL elements
991  block.associatedElements(iTrack, linkData, hcalElems, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
992 
993  if (!hcalElems.empty())
994  LogTrace("PFAlgo|elementLoop") << "Track linked back to HCAL due to ECAL sharing with other tracks";
995 }
#define LogTrace(id)

◆ setCandConnectorParameters() [1/2]

void PFAlgo::setCandConnectorParameters ( const edm::ParameterSet iCfgCandConnector)
inline

Definition at line 68 of file PFAlgo.h.

References connector_, HLT_2022v12_cff::iCfgCandConnector, and PFCandConnector::setParameters().

68  {
70  }
PFCandConnector connector_
Definition: PFAlgo.h:290
void setParameters(const edm::ParameterSet &iCfgCandConnector)

◆ setCandConnectorParameters() [2/2]

void PFAlgo::setCandConnectorParameters ( bool  bCorrect,
bool  bCalibPrimary,
double  dptRel_PrimaryTrack,
double  dptRel_MergedTrack,
double  ptErrorSecondary,
const std::vector< double > &  nuclCalibFactors 
)
inline

◆ setDisplacedVerticesParameters()

void PFAlgo::setDisplacedVerticesParameters ( bool  rejectTracks_Bad,
bool  rejectTracks_Step45,
bool  usePFNuclearInteractions,
bool  usePFConversions,
bool  usePFDecays,
double  dptRel_DispVtx 
)

Definition at line 96 of file PFAlgo.cc.

References HLT_2022v12_cff::dptRel_DispVtx, dptRel_DispVtx_, HLT_2022v12_cff::rejectTracks_Bad, rejectTracks_Bad_, HLT_2022v12_cff::rejectTracks_Step45, rejectTracks_Step45_, HLT_2022v12_cff::usePFConversions, usePFConversions_, HLT_2022v12_cff::usePFDecays, usePFDecays_, HLT_2022v12_cff::usePFNuclearInteractions, and usePFNuclearInteractions_.

◆ setEGammaCollections()

void PFAlgo::setEGammaCollections ( const edm::View< reco::PFCandidate > &  pfEgammaCandidates,
const edm::ValueMap< reco::GsfElectronRef > &  valueMapGedElectrons,
const edm::ValueMap< reco::PhotonRef > &  valueMapGedPhotons 
)

Definition at line 76 of file PFAlgo.cc.

References gedPhotonCore_cfi::pfEgammaCandidates, pfEgammaCandidates_, useEGammaFilters_, valueMapGedElectrons_, and valueMapGedPhotons_.

Referenced by PFProducer::produce().

78  {
79  if (useEGammaFilters_) {
81  valueMapGedElectrons_ = &valueMapGedElectrons;
82  valueMapGedPhotons_ = &valueMapGedPhotons;
83  }
84 }
const edm::ValueMap< reco::PhotonRef > * valueMapGedPhotons_
Definition: PFAlgo.h:269
const edm::View< reco::PFCandidate > * pfEgammaCandidates_
Definition: PFAlgo.h:267
bool useEGammaFilters_
Variables for NEW EGAMMA selection.
Definition: PFAlgo.h:265
const edm::ValueMap< reco::GsfElectronRef > * valueMapGedElectrons_
Definition: PFAlgo.h:268

◆ setEGammaParameters()

void PFAlgo::setEGammaParameters ( bool  use_EGammaFilters,
bool  useProtectionsForJetMET 
)

Definition at line 68 of file PFAlgo.cc.

References useEGammaFilters_, HLT_2022v12_cff::useProtectionsForJetMET, and useProtectionsForJetMET_.

68  {
69  useEGammaFilters_ = use_EGammaFilters;
70 
71  if (!useEGammaFilters_)
72  return;
73 
75 }
bool useProtectionsForJetMET_
Definition: PFAlgo.h:266
bool useEGammaFilters_
Variables for NEW EGAMMA selection.
Definition: PFAlgo.h:265

◆ setEGElectronCollection()

void PFAlgo::setEGElectronCollection ( const reco::GsfElectronCollection egelectrons)

◆ setHcalDepthInfo()

void PFAlgo::setHcalDepthInfo ( reco::PFCandidate cand,
const reco::PFCluster cluster 
) const
private

Definition at line 3331 of file PFAlgo.cc.

References Exception, ntuplemaker::fill, DetId::Hcal, mps_fire::i, and reco::PFCluster::recHitFractions().

Referenced by createCandidatesHCAL(), and reconstructCluster().

3331  {
3332  std::array<double, 7> energyPerDepth;
3333  std::fill(energyPerDepth.begin(), energyPerDepth.end(), 0.0);
3334  for (auto& hitRefAndFrac : cluster.recHitFractions()) {
3335  const auto& hit = *hitRefAndFrac.recHitRef();
3336  if (DetId(hit.detId()).det() == DetId::Hcal) {
3337  if (hit.depth() == 0) {
3338  edm::LogWarning("setHcalDepthInfo") << "Depth zero found";
3339  continue;
3340  }
3341  if (hit.depth() < 1 || hit.depth() > 7) {
3342  throw cms::Exception("CorruptData") << "Bogus depth " << hit.depth() << " at detid " << hit.detId() << "\n";
3343  }
3344  energyPerDepth[hit.depth() - 1] += hitRefAndFrac.fraction() * hit.energy();
3345  }
3346  }
3347  double sum = std::accumulate(energyPerDepth.begin(), energyPerDepth.end(), 0.);
3348  std::array<float, 7> depthFractions;
3349  if (sum > 0) {
3350  for (unsigned int i = 0; i < depthFractions.size(); ++i) {
3351  depthFractions[i] = energyPerDepth[i] / sum;
3352  }
3353  } else {
3354  std::fill(depthFractions.begin(), depthFractions.end(), 0.f);
3355  }
3356  cand.setHcalDepthEnergyFractions(depthFractions);
3357 }
const std::vector< reco::PFRecHitFraction > & recHitFractions() const
vector of rechit fractions
Definition: PFCluster.h:65
Definition: DetId.h:17
Log< level::Warning, false > LogWarning

◆ setHOTag()

void PFAlgo::setHOTag ( bool  ho)
inline

Definition at line 65 of file PFAlgo.h.

References photonIsolationHIProducer_cfi::ho, and useHO_.

◆ setMuonHandle()

void PFAlgo::setMuonHandle ( const edm::Handle< reco::MuonCollection > &  muons)
inline

Definition at line 66 of file PFAlgo.h.

References muonHandle_, and PDWG_BPHSkim_cff::muons.

Referenced by PFProducer::produce().

66 { muonHandle_ = muons; }
edm::Handle< reco::MuonCollection > muonHandle_
Definition: PFAlgo.h:331

◆ setPFVertexParameters()

void PFAlgo::setPFVertexParameters ( bool  useVertex,
reco::VertexCollection const &  primaryVertices 
)

Definition at line 110 of file PFAlgo.cc.

References nVtx_, pfmu_, primaryVertex_, HLT_2022v12_cff::primaryVertices, VerticesFromLeptons_cfi::useVertex, useVertices_, and bphysicsOniaDQM_cfi::vertex.

Referenced by PFProducer::produce().

110  {
112 
113  //Set the vertices for muon cleaning
114  pfmu_->setInputsForCleaning(primaryVertices);
115 
116  //Now find the primary vertex!
117  bool primaryVertexFound = false;
118  nVtx_ = primaryVertices.size();
119  for (auto const& vertex : primaryVertices) {
120  if (vertex.isValid() && (!vertex.isFake())) {
122  primaryVertexFound = true;
123  break;
124  }
125  }
126  //Use vertices if the user wants to but only if it exists a good vertex
127  useVertices_ = useVertex && primaryVertexFound;
128 }
int nVtx_
Definition: PFAlgo.h:286
reco::Vertex primaryVertex_
Definition: PFAlgo.h:328
bool useVertices_
Definition: PFAlgo.h:329
std::unique_ptr< PFMuonAlgo > pfmu_
Definition: PFAlgo.h:262

◆ setPostHFCleaningParameters()

void PFAlgo::setPostHFCleaningParameters ( bool  postHFCleaning,
const edm::ParameterSet pfHFCleaningParams 
)

Definition at line 86 of file PFAlgo.cc.

References edm::ParameterSet::getParameter(), maxDeltaPhiPt_, maxSignificance_, minDeltaMet_, minHFCleaningPt_, minSignificance_, minSignificanceReduction_, HLT_2022v12_cff::postHFCleaning, and postHFCleaning_.

86  {
88  minHFCleaningPt_ = pfHFCleaningParams.getParameter<double>("minHFCleaningPt");
89  minSignificance_ = pfHFCleaningParams.getParameter<double>("minSignificance");
90  maxSignificance_ = pfHFCleaningParams.getParameter<double>("maxSignificance");
91  minSignificanceReduction_ = pfHFCleaningParams.getParameter<double>("minSignificanceReduction");
92  maxDeltaPhiPt_ = pfHFCleaningParams.getParameter<double>("maxDeltaPhiPt");
93  minDeltaMet_ = pfHFCleaningParams.getParameter<double>("minDeltaMet");
94 }
T getParameter(std::string const &) const
Definition: ParameterSet.h:303
double maxDeltaPhiPt_
Definition: PFAlgo.h:323
double maxSignificance_
Definition: PFAlgo.h:321
double minHFCleaningPt_
Definition: PFAlgo.h:319
double minDeltaMet_
Definition: PFAlgo.h:324
double minSignificance_
Definition: PFAlgo.h:320
double minSignificanceReduction_
Definition: PFAlgo.h:322
bool postHFCleaning_
Definition: PFAlgo.h:317

Friends And Related Function Documentation

◆ operator<<

std::ostream& operator<< ( std::ostream &  out,
const PFAlgo algo 
)
friend

Member Data Documentation

◆ calibration_

PFEnergyCalibration& PFAlgo::calibration_
private

Definition at line 257 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), createCandidatesHCALUnlinked(), and recoTracksNotHCAL().

◆ connector_

PFCandConnector PFAlgo::connector_
private

A tool used for a postprocessing of displaced vertices based on reconstructed PFCandidates

Definition at line 290 of file PFAlgo.h.

Referenced by makeConnectedCandidates(), and setCandConnectorParameters().

◆ dptRel_DispVtx_

double PFAlgo::dptRel_DispVtx_
private

Maximal relative uncertainty on the tracks going to or incoming from the displcaed vertex to be used in the PFAlgo

Definition at line 285 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), reconstructTrack(), recoTracksNotHCAL(), and setDisplacedVerticesParameters().

◆ factors45_

std::vector<double> PFAlgo::factors45_
private

Definition at line 298 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), and PFAlgo().

◆ goodPixelTrackDeadHcal_chi2n_

float PFAlgo::goodPixelTrackDeadHcal_chi2n_
private

Definition at line 310 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_dxy_

float PFAlgo::goodPixelTrackDeadHcal_dxy_
private

Definition at line 313 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_dz_

float PFAlgo::goodPixelTrackDeadHcal_dz_
private

Definition at line 314 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_maxLost3Hit_

int PFAlgo::goodPixelTrackDeadHcal_maxLost3Hit_
private

Definition at line 311 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_maxLost4Hit_

int PFAlgo::goodPixelTrackDeadHcal_maxLost4Hit_
private

Definition at line 312 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_maxPt_

float PFAlgo::goodPixelTrackDeadHcal_maxPt_
private

Definition at line 308 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_minEta_

float PFAlgo::goodPixelTrackDeadHcal_minEta_
private

Definition at line 307 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_ptErrRel_

float PFAlgo::goodPixelTrackDeadHcal_ptErrRel_
private

Definition at line 309 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodTrackDeadHcal_chi2n_

float PFAlgo::goodTrackDeadHcal_chi2n_
private

Definition at line 302 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodTrackDeadHcal_dxy_

float PFAlgo::goodTrackDeadHcal_dxy_
private

Definition at line 305 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodTrackDeadHcal_layers_

int PFAlgo::goodTrackDeadHcal_layers_
private

Definition at line 303 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodTrackDeadHcal_ptErrRel_

float PFAlgo::goodTrackDeadHcal_ptErrRel_
private

Variables for track cleaning in bad HCal areas.

Definition at line 301 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodTrackDeadHcal_validFr_

float PFAlgo::goodTrackDeadHcal_validFr_
private

Definition at line 304 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ maxDeltaPhiPt_

double PFAlgo::maxDeltaPhiPt_
private

Definition at line 323 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

◆ maxSignificance_

double PFAlgo::maxSignificance_
private

Definition at line 321 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

◆ minDeltaMet_

double PFAlgo::minDeltaMet_
private

Definition at line 324 of file PFAlgo.h.

Referenced by checkCleaning(), postCleaning(), and setPostHFCleaningParameters().

◆ minHFCleaningPt_

double PFAlgo::minHFCleaningPt_
private

Definition at line 319 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

◆ minSignificance_

double PFAlgo::minSignificance_
private

Definition at line 320 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

◆ minSignificanceReduction_

double PFAlgo::minSignificanceReduction_
private

Definition at line 322 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

◆ muonECAL_

std::vector<double> PFAlgo::muonECAL_
private

Definition at line 294 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), PFAlgo(), and recoTracksNotHCAL().

◆ muonHandle_

edm::Handle<reco::MuonCollection> PFAlgo::muonHandle_
private

Definition at line 331 of file PFAlgo.h.

Referenced by reconstructParticles(), and setMuonHandle().

◆ muonHCAL_

std::vector<double> PFAlgo::muonHCAL_
private

Variables for muons and fakes.

Definition at line 293 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), and PFAlgo().

◆ muonHO_

std::vector<double> PFAlgo::muonHO_
private

Definition at line 295 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), and PFAlgo().

◆ nSigmaECAL_

const double PFAlgo::nSigmaECAL_
private

number of sigma to judge energy excess in ECAL

Definition at line 245 of file PFAlgo.h.

Referenced by recoTracksNotHCAL().

◆ nSigmaEConstHCAL

const double PFAlgo::nSigmaEConstHCAL = 100.
private

Definition at line 334 of file PFAlgo.h.

Referenced by nSigmaHCAL().

◆ nSigmaEConstHFEM

const double PFAlgo::nSigmaEConstHFEM = 100.
private

Definition at line 335 of file PFAlgo.h.

Referenced by nSigmaHFEM().

◆ nSigmaEConstHFHAD

const double PFAlgo::nSigmaEConstHFHAD = 100.
private

Definition at line 336 of file PFAlgo.h.

Referenced by nSigmaHFHAD().

◆ nSigmaHCAL_

const double PFAlgo::nSigmaHCAL_
private

number of sigma to judge energy excess in HCAL

Definition at line 248 of file PFAlgo.h.

Referenced by nSigmaHCAL().

◆ nSigmaHFEM_

const double PFAlgo::nSigmaHFEM_
private

number of sigma to judge energy excess in HF

Definition at line 251 of file PFAlgo.h.

Referenced by nSigmaHFEM().

◆ nSigmaHFHAD_

const double PFAlgo::nSigmaHFHAD_
private

Definition at line 252 of file PFAlgo.h.

Referenced by nSigmaHFHAD().

◆ nSigmaTRACK_

double PFAlgo::nSigmaTRACK_
private

Definition at line 296 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), PFAlgo(), and recoTracksNotHCAL().

◆ nVtx_

int PFAlgo::nVtx_
private

Definition at line 286 of file PFAlgo.h.

Referenced by egammaFilters(), and setPFVertexParameters().

◆ pfCandidates_

std::unique_ptr<reco::PFCandidateCollection> PFAlgo::pfCandidates_
private

◆ pfCleanedCandidates_

reco::PFCandidateCollection PFAlgo::pfCleanedCandidates_
private

Definition at line 228 of file PFAlgo.h.

Referenced by getCleanedCandidates(), postCleaning(), and reconstructParticles().

◆ pfEgammaCandidates_

const edm::View<reco::PFCandidate>* PFAlgo::pfEgammaCandidates_
private

Definition at line 267 of file PFAlgo.h.

Referenced by egammaFilters(), and setEGammaCollections().

◆ pfmu_

std::unique_ptr<PFMuonAlgo> PFAlgo::pfmu_
private

◆ postHFCleaning_

bool PFAlgo::postHFCleaning_
private

Definition at line 317 of file PFAlgo.h.

Referenced by reconstructParticles(), and setPostHFCleaningParameters().

◆ postMuonCleaning_

bool PFAlgo::postMuonCleaning_
private

Definition at line 318 of file PFAlgo.h.

◆ primaryVertex_

reco::Vertex PFAlgo::primaryVertex_
private

◆ ptError_

double PFAlgo::ptError_
private

Definition at line 297 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), PFAlgo(), and recoTracksNotHCAL().

◆ rejectTracks_Bad_

bool PFAlgo::rejectTracks_Bad_
private

Flags to use the protection against fakes and not reconstructed displaced vertices

Definition at line 276 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), and setDisplacedVerticesParameters().

◆ rejectTracks_Step45_

bool PFAlgo::rejectTracks_Step45_
private

◆ resolHF_square_

const std::vector<double> PFAlgo::resolHF_square_
private

Definition at line 255 of file PFAlgo.h.

Referenced by hfEnergyResolution(), and PFAlgo().

◆ thepfEnergyCalibrationHF_

PFEnergyCalibrationHF& PFAlgo::thepfEnergyCalibrationHF_
private

Definition at line 258 of file PFAlgo.h.

Referenced by createCandidatesHF().

◆ useBestMuonTrack_

double PFAlgo::useBestMuonTrack_
private

Definition at line 325 of file PFAlgo.h.

◆ useEGammaFilters_

bool PFAlgo::useEGammaFilters_
private

Variables for NEW EGAMMA selection.

Definition at line 265 of file PFAlgo.h.

Referenced by processBlock(), setEGammaCollections(), and setEGammaParameters().

◆ useHO_

bool PFAlgo::useHO_
private

Definition at line 260 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), createCandidatesHCALUnlinked(), decideType(), and setHOTag().

◆ usePFConversions_

bool PFAlgo::usePFConversions_
private

Definition at line 280 of file PFAlgo.h.

Referenced by processBlock(), and setDisplacedVerticesParameters().

◆ usePFDecays_

bool PFAlgo::usePFDecays_
private

Definition at line 281 of file PFAlgo.h.

Referenced by isFromSecInt(), and setDisplacedVerticesParameters().

◆ usePFMuonMomAssign_

bool PFAlgo::usePFMuonMomAssign_
private

Definition at line 272 of file PFAlgo.h.

◆ usePFNuclearInteractions_

bool PFAlgo::usePFNuclearInteractions_
private

Definition at line 279 of file PFAlgo.h.

Referenced by isFromSecInt(), and setDisplacedVerticesParameters().

◆ useProtectionsForJetMET_

bool PFAlgo::useProtectionsForJetMET_
private

Definition at line 266 of file PFAlgo.h.

Referenced by egammaFilters(), and setEGammaParameters().

◆ useVertices_

bool PFAlgo::useVertices_ = false
private

Definition at line 329 of file PFAlgo.h.

Referenced by reconstructCluster(), and setPFVertexParameters().

◆ valueMapGedElectrons_

const edm::ValueMap<reco::GsfElectronRef>* PFAlgo::valueMapGedElectrons_
private

Definition at line 268 of file PFAlgo.h.

Referenced by setEGammaCollections().

◆ valueMapGedPhotons_

const edm::ValueMap<reco::PhotonRef>* PFAlgo::valueMapGedPhotons_
private

Definition at line 269 of file PFAlgo.h.

Referenced by setEGammaCollections().