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_2023v12_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:3404
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:3409
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:3385
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 3439 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().

3445  {
3446  // Find all PS clusters with type psElement associated to ECAL cluster iEcal,
3447  // within all PFBlockElement "elements" of a given PFBlock "block"
3448  // psElement can be reco::PFBlockElement::PS1 or reco::PFBlockElement::PS2
3449  // Returns a vector of PS cluster energies, and updates the "active" vector.
3450 
3451  // Find all PS clusters linked to the iEcal cluster
3452  std::multimap<double, unsigned> sortedPS;
3453  block.associatedElements(iEcal, linkData, sortedPS, psElementType, reco::PFBlock::LINKTEST_ALL);
3454 
3455  // Loop over these PS clusters
3456  for (auto const& ps : sortedPS) {
3457  // CLuster index and distance to iEcal
3458  unsigned iPS = ps.second;
3459  // double distPS = ps.first;
3460 
3461  // Ignore clusters already in use
3462  if (!active[iPS])
3463  continue;
3464 
3465  // Check that this cluster is not closer to another ECAL cluster
3466  std::multimap<double, unsigned> sortedECAL;
3467  block.associatedElements(iPS, linkData, sortedECAL, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
3468  unsigned jEcal = sortedECAL.begin()->second;
3469  if (jEcal != iEcal)
3470  continue;
3471 
3472  // Update PS energy
3473  PFBlockElement::Type pstype = elements[iPS].type();
3474  assert(pstype == psElementType);
3475  PFClusterRef psclusterref = elements[iPS].clusterRef();
3476  assert(!psclusterref.isNull());
3477  psEne[0] += psclusterref->energy();
3478  active[iPS] = false;
3479  }
3480 }
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:3482
#define LogTrace(id)
unsigned reconstructTrack(const reco::PFBlockElement &elt, bool allowLoose=false)
Definition: PFAlgo.cc:3155

◆ checkCleaning()

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

Check HF Cleaning.

Definition at line 3607 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().

3607  {
3608  // No hits to recover, leave.
3609  if (cleanedHits.empty())
3610  return;
3611 
3612  //Compute met and met significance (met/sqrt(SumEt))
3613  double metX = 0.;
3614  double metY = 0.;
3615  double sumet = 0;
3616  std::vector<unsigned int> hitsToBeAdded;
3617  for (auto const& pfc : *pfCandidates_) {
3618  metX += pfc.px();
3619  metY += pfc.py();
3620  sumet += pfc.pt();
3621  }
3622  double met2 = metX * metX + metY * metY;
3623  double met2_Original = met2;
3624  // Select events with large MET significance.
3625  // double significance = std::sqrt(met2/sumet);
3626  // double significanceCor = significance;
3627  double metXCor = metX;
3628  double metYCor = metY;
3629  double sumetCor = sumet;
3630  double met2Cor = met2;
3631  bool next = true;
3632  unsigned iCor = 1E9;
3633 
3634  // Find the cleaned hit with the largest effect on the MET
3635  while (next) {
3636  double metReduc = -1.;
3637  // Loop on the candidates
3638  for (unsigned i = 0; i < cleanedHits.size(); ++i) {
3639  const PFRecHit& hit = cleanedHits[i];
3640  double length = std::sqrt(hit.position().mag2());
3641  double px = hit.energy() * hit.position().x() / length;
3642  double py = hit.energy() * hit.position().y() / length;
3643  double pt = std::sqrt(px * px + py * py);
3644 
3645  // Check that it is not already scheculed to be cleaned
3646  bool skip = false;
3647  for (unsigned int hitIdx : hitsToBeAdded) {
3648  if (i == hitIdx)
3649  skip = true;
3650  if (skip)
3651  break;
3652  }
3653  if (skip)
3654  continue;
3655 
3656  // Now add the candidate to the MET
3657  double metXInt = metX + px;
3658  double metYInt = metY + py;
3659  double sumetInt = sumet + pt;
3660  double met2Int = metXInt * metXInt + metYInt * metYInt;
3661 
3662  // And check if it could contribute to a MET reduction
3663  if (met2Int < met2Cor) {
3664  metXCor = metXInt;
3665  metYCor = metYInt;
3666  metReduc = (met2 - met2Int) / met2Int;
3667  met2Cor = met2Int;
3668  sumetCor = sumetInt;
3669  // significanceCor = std::sqrt(met2Cor/sumetCor);
3670  iCor = i;
3671  }
3672  }
3673  //
3674  // If the MET must be significanly reduced, schedule the candidate to be added
3675  //
3676  if (metReduc > minDeltaMet_) {
3677  hitsToBeAdded.push_back(iCor);
3678  metX = metXCor;
3679  metY = metYCor;
3680  sumet = sumetCor;
3681  met2 = met2Cor;
3682  } else {
3683  // Otherwise just stop the loop
3684  next = false;
3685  }
3686  }
3687  //
3688  // At least 10 GeV MET reduction
3689  if (std::sqrt(met2_Original) - std::sqrt(met2) > 5.) {
3690  LogTrace("PFAlgo|checkCleaning") << hitsToBeAdded.size() << " hits were re-added ";
3691  LogTrace("PFAlgo|checkCleaning") << "MET reduction = " << std::sqrt(met2_Original) << " -> " << std::sqrt(met2Cor)
3692  << " = " << std::sqrt(met2Cor) - std::sqrt(met2_Original);
3693  LogTrace("PFAlgo|checkCleaning") << "Added after cleaning check : ";
3694  for (unsigned int hitIdx : hitsToBeAdded) {
3695  const PFRecHit& hit = cleanedHits[hitIdx];
3696  PFCluster cluster(hit.layer(), hit.energy(), hit.position().x(), hit.position().y(), hit.position().z());
3697  reconstructCluster(cluster, hit.energy());
3698  LogTrace("PFAlgo|checkCleaning") << pfCandidates_->back() << ". time = " << hit.time();
3699  }
3700  }
3701 } //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:3234
#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:128
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 3026 of file PFAlgo.cc.

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

Referenced by processBlock().

3032  {
3033  LogTrace("PFAlgo|createCandidatesECAL")
3034  << "start of function PFAlgo::createCandidatesECAL(), ecalIs.size()=" << inds.ecalIs.size();
3035 
3036  // --------------- loop ecal ------------------
3037 
3038  // for each ecal element iEcal = ecalIs[i] in turn:
3039 
3040  for (unsigned i = 0; i < inds.ecalIs.size(); i++) {
3041  unsigned iEcal = inds.ecalIs[i];
3042 
3043  LogTrace("PFAlgo|createCandidatesECAL") << "elements[" << iEcal << "]=" << elements[iEcal];
3044 
3045  if (!active[iEcal]) {
3046  LogTrace("PFAlgo|createCandidatesECAL") << "iEcal=" << iEcal << " not active";
3047  continue;
3048  }
3049 
3050  PFBlockElement::Type type = elements[iEcal].type();
3052 
3053  PFClusterRef clusterref = elements[iEcal].clusterRef();
3054  assert(!clusterref.isNull());
3055 
3056  active[iEcal] = false;
3057 
3058  float ecalEnergy = clusterref->correctedEnergy();
3059  // float ecalEnergy = calibration_.energyEm( clusterref->energy() );
3060  double particleEnergy = ecalEnergy;
3061 
3062  auto& cand = (*pfCandidates_)[reconstructCluster(*clusterref, particleEnergy)];
3063 
3064  cand.setEcalEnergy(clusterref->energy(), ecalEnergy);
3065  cand.setHcalEnergy(0., 0.);
3066  cand.setHoEnergy(0., 0.);
3067  cand.setPs1Energy(0.);
3068  cand.setPs2Energy(0.);
3069  cand.addElementInBlock(blockref, iEcal);
3070 
3071  } // end loop on ecal elements iEcal = ecalIs[i]
3072  LogTrace("PFAlgo|createCandidatesECAL") << "end of function PFALgo::createCandidatesECAL";
3073 }
unsigned reconstructCluster(const reco::PFCluster &cluster, double particleEnergy, bool useDirection=false, double particleX=0., double particleY=0., double particleZ=0.)
Definition: PFAlgo.cc:3234
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 1700 of file PFAlgo.cc.

References a, funct::abs(), cms::cuda::assert(), b, groupFilesInBlocks::block, calibration_, l1ctLayer1_cff::caloResolution, HPSPFTauProducerPuppi_cfi::chargedHadron, RPCNoise_example::check, Calorimetry_cff::dp, dptRel_DispVtx_, MillePedeFileConverter_cfg::e, ECAL, reco::PFBlockElement::ECAL, RecoEcal_cff::ecalClusters, electrons_cff::ecalEnergy, bookConverter::elements, PFEnergyCalibration::energyEmHad(), PFTrackAlgoTools::errorScale(), factors45_, dqmdumpme::first, HLT_2023v12_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().

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

References cms::cuda::assert(), groupFilesInBlocks::block, calibration_, ECAL, reco::PFBlockElement::ECAL, electrons_cff::ecalEnergy, 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().

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

References cms::cuda::assert(), groupFilesInBlocks::block, fftjetpileupestimator_calo_uncalib_cfi::c0, alignmentValidation::c1, l1ctLayer1_cff::caloResolution, 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().

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

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

Referenced by elementLoop().

1204  {
1205  switch (type) {
1206  case PFBlockElement::TRACK:
1207  if (active[iEle]) {
1208  inds.trackIs.push_back(iEle);
1209  LogTrace("PFAlgo|decideType") << "TRACK, stored index, continue";
1210  }
1211  break;
1212  case PFBlockElement::ECAL:
1213  if (active[iEle]) {
1214  inds.ecalIs.push_back(iEle);
1215  LogTrace("PFAlgo|decideType") << "ECAL, stored index, continue";
1216  }
1217  return 1; //continue
1218  case PFBlockElement::HCAL:
1219  if (active[iEle]) {
1220  if (elements[iEle].clusterRef()->flags() & reco::CaloCluster::badHcalMarker) {
1221  LogTrace("PFAlgo|decideType") << "HCAL DEAD AREA: remember and skip.";
1222  active[iEle] = false;
1223  deadArea[iEle] = true;
1224  return 1; //continue
1225  }
1226  inds.hcalIs.push_back(iEle);
1227  LogTrace("PFAlgo|decideType") << "HCAL, stored index, continue";
1228  }
1229  return 1; //continue
1230  case PFBlockElement::HO:
1231  if (useHO_) {
1232  if (active[iEle]) {
1233  inds.hoIs.push_back(iEle);
1234  LogTrace("PFAlgo|decideType") << "HO, stored index, continue";
1235  }
1236  }
1237  return 1; //continue
1238  case PFBlockElement::HFEM:
1239  if (active[iEle]) {
1240  inds.hfEmIs.push_back(iEle);
1241  LogTrace("PFAlgo|decideType") << "HFEM, stored index, continue";
1242  }
1243  return 1; //continue
1244  case PFBlockElement::HFHAD:
1245  if (active[iEle]) {
1246  inds.hfHadIs.push_back(iEle);
1247  LogTrace("PFAlgo|decideType") << "HFHAD, stored index, continue";
1248  }
1249  return 1; //continue
1250  default:
1251  return 1; //continue
1252  }
1253  LogTrace("PFAlgo|decideType") << "Did not match type to anything, return 0";
1254  return 0;
1255 }
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:134
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:583
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:130
double y() const
y coordinate
Definition: Vertex.h:132
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 995 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().

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

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

Referenced by createCandidatesHF().

3392  {
3393  // Add a protection
3394  clusterEnergyHF = std::max(clusterEnergyHF, 1.);
3395 
3396  double resol =
3397  sqrt(resolHF_square_[0] / clusterEnergyHF + resolHF_square_[1] + resolHF_square_[2] / pow(clusterEnergyHF, 2));
3398  // 0: stochastic term, 1: constant term, 2: noise term
3399  // Note: resolHF_square_[0,1,2] should be already squared
3400 
3401  return resol;
3402 }
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 3482 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().

3482  {
3485  // reco::PFBlockElement::TrackType T_FROM_GAMMACONV = reco::PFBlockElement::T_FROM_GAMMACONV;
3487 
3488  bool bPrimary = (order.find("primary") != string::npos);
3489  bool bSecondary = (order.find("secondary") != string::npos);
3490  bool bAll = (order.find("all") != string::npos);
3491 
3492  bool isToDisp = usePFNuclearInteractions_ && eTrack.trackType(T_TO_DISP);
3493  bool isFromDisp = usePFNuclearInteractions_ && eTrack.trackType(T_FROM_DISP);
3494 
3495  if (bPrimary && isToDisp)
3496  return true;
3497  if (bSecondary && isFromDisp)
3498  return true;
3499  if (bAll && (isToDisp || isFromDisp))
3500  return true;
3501 
3502  // bool isFromConv = usePFConversions_ && eTrack.trackType(T_FROM_GAMMACONV);
3503 
3504  // if ((bAll || bSecondary)&& isFromConv) return true;
3505 
3506  bool isFromDecay = (bAll || bSecondary) && usePFDecays_ && eTrack.trackType(T_FROM_V0);
3507 
3508  return isFromDecay;
3509 }
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 3375 of file PFAlgo.cc.

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

Referenced by createCandidatesHCAL(), and recoTracksNotHCAL().

3375  {
3376  // Add a protection
3377  clusterEnergyHCAL = std::max(clusterEnergyHCAL, 1.);
3378 
3379  double resol = fabs(eta) < 1.48 ? sqrt(1.02 * 1.02 / clusterEnergyHCAL + 0.065 * 0.065)
3380  : sqrt(1.20 * 1.20 / clusterEnergyHCAL + 0.028 * 0.028);
3381 
3382  return resol;
3383 }
T sqrt(T t)
Definition: SSEVec.h:19

◆ nSigmaHCAL()

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

Definition at line 3385 of file PFAlgo.cc.

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

Referenced by createCandidatesHCAL().

3385  {
3386  double nS = fabs(eta) < 1.48 ? nSigmaHCAL_ * (1. + exp(-clusterEnergyHCAL / nSigmaEConstHCAL))
3387  : nSigmaHCAL_ * (1. + exp(-clusterEnergyHCAL / nSigmaEConstHCAL));
3388 
3389  return nS;
3390 }
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 3404 of file PFAlgo.cc.

References JetChargeProducer_cfi::exp, nSigmaEConstHFEM, and nSigmaHFEM_.

Referenced by createCandidatesHF().

3404  {
3405  double nS = nSigmaHFEM_ * (1. + exp(-clusterEnergyHF / nSigmaEConstHFEM));
3406  return nS;
3407 }
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 3409 of file PFAlgo.cc.

References JetChargeProducer_cfi::exp, nSigmaEConstHFHAD, and nSigmaHFHAD_.

Referenced by createCandidatesHF().

3409  {
3410  double nS = nSigmaHFHAD_ * (1. + exp(-clusterEnergyHF / nSigmaEConstHFHAD));
3411  return nS;
3412 }
const double nSigmaEConstHFHAD
Definition: PFAlgo.h:336
const double nSigmaHFHAD_
Definition: PFAlgo.h:252

◆ postCleaning()

void PFAlgo::postCleaning ( )
private

Definition at line 3511 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().

3511  {
3512  //Compute met and met significance (met/sqrt(SumEt))
3513  double metX = 0.;
3514  double metY = 0.;
3515  double sumet = 0;
3516  std::vector<unsigned int> pfCandidatesToBeRemoved;
3517  for (auto const& pfc : *pfCandidates_) {
3518  metX += pfc.px();
3519  metY += pfc.py();
3520  sumet += pfc.pt();
3521  }
3522  double met2 = metX * metX + metY * metY;
3523  // Select events with large MET significance.
3524  double significance = std::sqrt(met2 / sumet);
3525  double significanceCor = significance;
3527  double metXCor = metX;
3528  double metYCor = metY;
3529  double sumetCor = sumet;
3530  double met2Cor = met2;
3531  double deltaPhi = 3.14159;
3532  double deltaPhiPt = 100.;
3533  bool next = true;
3534  unsigned iCor = 1E9;
3535 
3536  // Find the HF candidate with the largest effect on the MET
3537  while (next) {
3538  double metReduc = -1.;
3539  // Loop on the candidates
3540  for (unsigned i = 0; i < pfCandidates_->size(); ++i) {
3541  const PFCandidate& pfc = (*pfCandidates_)[i];
3542 
3543  // Check that the pfCandidate is in the HF
3545  continue;
3546 
3547  // Check if has meaningful pt
3548  if (pfc.pt() < minHFCleaningPt_)
3549  continue;
3550 
3551  // Check that it is not already scheculed to be cleaned
3552  const bool skip = std::any_of(
3553  pfCandidatesToBeRemoved.begin(), pfCandidatesToBeRemoved.end(), [&](unsigned int j) { return i == j; });
3554  if (skip)
3555  continue;
3556 
3557  // Check that the pt and the MET are aligned
3558  deltaPhi = std::acos((metX * pfc.px() + metY * pfc.py()) / (pfc.pt() * std::sqrt(met2)));
3559  deltaPhiPt = deltaPhi * pfc.pt();
3560  if (deltaPhiPt > maxDeltaPhiPt_)
3561  continue;
3562 
3563  // Now remove the candidate from the MET
3564  double metXInt = metX - pfc.px();
3565  double metYInt = metY - pfc.py();
3566  double sumetInt = sumet - pfc.pt();
3567  double met2Int = metXInt * metXInt + metYInt * metYInt;
3568  if (met2Int < met2Cor) {
3569  metXCor = metXInt;
3570  metYCor = metYInt;
3571  metReduc = (met2 - met2Int) / met2Int;
3572  met2Cor = met2Int;
3573  sumetCor = sumetInt;
3574  significanceCor = std::sqrt(met2Cor / sumetCor);
3575  iCor = i;
3576  }
3577  }
3578  //
3579  // If the MET must be significanly reduced, schedule the candidate to be cleaned
3580  if (metReduc > minDeltaMet_) {
3581  pfCandidatesToBeRemoved.push_back(iCor);
3582  metX = metXCor;
3583  metY = metYCor;
3584  sumet = sumetCor;
3585  met2 = met2Cor;
3586  } else {
3587  // Otherwise just stop the loop
3588  next = false;
3589  }
3590  }
3591  //
3592  // The significance must be significantly reduced to indeed clean the candidates
3593  if (significance - significanceCor > minSignificanceReduction_ && significanceCor < maxSignificance_) {
3594  edm::LogInfo("PFAlgo|postCleaning") << "Significance reduction = " << significance << " -> " << significanceCor
3595  << " = " << significanceCor - significance;
3596  for (unsigned int toRemove : pfCandidatesToBeRemoved) {
3597  edm::LogInfo("PFAlgo|postCleaning") << "Removed : " << (*pfCandidates_)[toRemove];
3598  pfCleanedCandidates_.push_back((*pfCandidates_)[toRemove]);
3599  (*pfCandidates_)[toRemove].rescaleMomentum(1E-6);
3600  //reco::PFCandidate::ParticleType unknown = reco::PFCandidate::X;
3601  //(*pfCandidates_)[toRemove].setParticleType(unknown);
3602  }
3603  }
3604  } //significance
3605 } //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:15
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 3075 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().

3078  {
3079  assert(!blockref.isNull());
3080  const reco::PFBlock& block = *blockref;
3081 
3082  LogTrace("PFAlgo|processBlock") << "start of function PFAlgo::processBlock, block=" << block;
3083 
3085  LogTrace("PFAlgo|processBlock") << "elements.size()=" << elements.size();
3086  // make a copy of the link data, which will be edited.
3087  PFBlock::LinkData linkData = block.linkData();
3088 
3089  // keep track of the elements which are still active.
3090  vector<bool> active(elements.size(), true);
3091 
3092  // //PFElectrons:
3093  // usePFElectrons_ external configurable parameter to set the usage of pf electron
3094  std::vector<reco::PFCandidate> tempElectronCandidates;
3095  tempElectronCandidates.clear();
3096 
3097  // New EGamma Reconstruction 10/10/2013
3098  if (useEGammaFilters_) {
3099  egammaFilters(blockref, active, pfegamma);
3100  } // end if use EGammaFilters
3101 
3102  //Lock extra conversion tracks not used by Photon Algo
3103  if (usePFConversions_) {
3104  conversionAlgo(elements, active);
3105  }
3106 
3107  // In the following elementLoop() function, the primary goal is to deal with tracks that are:
3108  // - not associated to an HCAL cluster
3109  // - not identified as an electron.
3110  // Those tracks should be predominantly relatively low energy charged
3111  // hadrons which are not detected in the ECAL.
3112 
3113  // The secondary goal is to prepare for the next loops
3114  // - The ecal and hcal elements are sorted in separate vectors in `ElementIndices inds`
3115  // which will be used as a base for the corresponding loops.
3116  // - For tracks which are connected to more than one HCAL cluster,
3117  // the links between the track and the cluster are cut for all clusters
3118  // but the closest one.
3119  // - HF only blocks ( HFEM, HFHAD, HFEM+HFAD) are identified
3120 
3121  // obsolete comments?
3122  // loop1:
3123  // - sort ecal and hcal elements in separate vectors
3124  // - for tracks:
3125  // - lock closest ecal cluster
3126  // - cut link to farthest hcal cluster, if more than 1.
3127 
3128  vector<bool> deadArea(elements.size(), false);
3129 
3130  // vectors to store element indices to ho, hcal and ecal elements, will be filled by elementLoop()
3131  ElementIndices inds;
3132 
3133  elementLoop(block, linkData, elements, active, blockref, inds, deadArea);
3134 
3135  // Reconstruct pfCandidate from HF (either EM-only, Had-only or both)
3136  // For phase2, process also pfblocks containing HF clusters and linked tracks
3137  if (!(inds.hfEmIs.empty() && inds.hfHadIs.empty())) {
3138  createCandidatesHF(block, linkData, elements, active, blockref, inds);
3139  if (inds.hcalIs.empty() && inds.ecalIs.empty())
3140  return;
3141  LogDebug("PFAlgo::processBlock")
3142  << "Block contains HF clusters, and also contains ECAL or HCAL clusters. Continue.\n"
3143  << block;
3144  }
3145 
3146  createCandidatesHCAL(block, linkData, elements, active, blockref, inds, deadArea);
3147  // COLINFEB16: now dealing with the HCAL elements that are not linked to any track
3148  createCandidatesHCALUnlinked(block, linkData, elements, active, blockref, inds, deadArea);
3149  createCandidatesECAL(block, linkData, elements, active, blockref, inds, deadArea);
3150 
3151  LogTrace("PFAlgo|processBlock") << "end of function PFAlgo::processBlock";
3152 } // 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:995
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:3026
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:1700
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:1257
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:2827
#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 3234 of file PFAlgo.cc.

References cms::cuda::assert(), ALCARECOTkAlJpsiMuMu_cff::charge, PFLayer::ECAL_BARREL, PFLayer::ECAL_ENDCAP, 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().

3239  {
3240  LogTrace("PFAlgo|reconstructCluster") << "start of function PFAlgo::reconstructCluster, cluster=" << cluster
3241  << "particleEnergy=" << particleEnergy << "useDirection=" << useDirection
3242  << "particleX=" << particleX << "particleY=" << particleY
3243  << "particleZ=" << particleZ;
3244 
3246 
3247  // need to convert the ::math::XYZPoint data member of the PFCluster class=
3248  // to a displacement vector:
3249 
3250  // Transform particleX,Y,Z to a position at ECAL/HCAL entrance
3251  double factor = 1.;
3252  if (useDirection) {
3253  switch (cluster.layer()) {
3254  case PFLayer::ECAL_BARREL:
3255  case PFLayer::HCAL_BARREL1:
3256  factor = std::sqrt(cluster.position().Perp2() / (particleX * particleX + particleY * particleY));
3257  break;
3258  case PFLayer::ECAL_ENDCAP:
3259  case PFLayer::HCAL_ENDCAP:
3260  case PFLayer::HF_HAD:
3261  case PFLayer::HF_EM:
3262  factor = cluster.position().Z() / particleZ;
3263  break;
3264  default:
3265  assert(0);
3266  }
3267  }
3268  //MIKE First of all let's check if we have vertex.
3269  ::math::XYZPoint vertexPos;
3270  if (useVertices_)
3272  else
3273  vertexPos = ::math::XYZPoint(0.0, 0.0, 0.0);
3274 
3275  ::math::XYZVector clusterPos(cluster.position().X() - vertexPos.X(),
3276  cluster.position().Y() - vertexPos.Y(),
3277  cluster.position().Z() - vertexPos.Z());
3278  ::math::XYZVector particleDirection(
3279  particleX * factor - vertexPos.X(), particleY * factor - vertexPos.Y(), particleZ * factor - vertexPos.Z());
3280 
3281  //::math::XYZVector clusterPos( cluster.position().X(), cluster.position().Y(),cluster.position().Z() );
3282  //::math::XYZVector particleDirection ( particleX, particleY, particleZ );
3283 
3284  clusterPos = useDirection ? particleDirection.Unit() : clusterPos.Unit();
3285  clusterPos *= particleEnergy;
3286 
3287  // clusterPos is now a vector along the cluster direction,
3288  // with a magnitude equal to the cluster energy.
3289 
3290  double mass = 0;
3291  ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> momentum(
3292  clusterPos.X(), clusterPos.Y(), clusterPos.Z(), mass);
3293  // mathcore is a piece of #$%
3295  // implicit constructor not allowed
3296  tmp = momentum;
3297 
3298  // Charge
3299  int charge = 0;
3300 
3301  // Type
3302  switch (cluster.layer()) {
3303  case PFLayer::ECAL_BARREL:
3304  case PFLayer::ECAL_ENDCAP:
3306  break;
3307  case PFLayer::HCAL_BARREL1:
3308  case PFLayer::HCAL_ENDCAP:
3309  particleType = PFCandidate::h0;
3310  break;
3311  case PFLayer::HF_HAD:
3312  particleType = PFCandidate::h_HF;
3313  break;
3314  case PFLayer::HF_EM:
3315  particleType = PFCandidate::egamma_HF;
3316  break;
3317  default:
3318  assert(0);
3319  }
3320 
3321  // The pf candidate
3322  LogTrace("PFAlgo|reconstructCluster") << "Creating PFCandidate charge=" << charge << ", type=" << particleType
3323  << ", pt=" << tmp.pt() << ", eta=" << tmp.eta() << ", phi=" << tmp.phi();
3325 
3326  // The position at ECAL entrance (well: watch out, it is not true
3327  // for HCAL clusters... to be fixed)
3328  pfCandidates_->back().setPositionAtECALEntrance(
3329  ::math::XYZPointF(cluster.position().X(), cluster.position().Y(), cluster.position().Z()));
3330 
3331  //Set the cnadidate Vertex
3332  pfCandidates_->back().setVertex(vertexPos);
3333 
3334  // depth info
3335  setHcalDepthInfo(pfCandidates_->back(), cluster);
3336 
3337  //*TODO* cluster time is not reliable at the moment, so only use track timing
3338 
3339  LogTrace("PFAlgo|reconstructCluster") << "** candidate: " << pfCandidates_->back();
3340 
3341  // returns index to the newly created PFCandidate
3342  return pfCandidates_->size() - 1;
3343 }
const math::XYZPoint & position() const
cluster centroid position
Definition: CaloCluster.h:153
ParticleType
particle types
Definition: PFCandidate.h:44
double z() const
z coordinate
Definition: Vertex.h:134
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:130
double y() const
y coordinate
Definition: Vertex.h:132
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:3345
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:3075
bool isValid() const
Definition: HandleBase.h:70
void postCleaning()
Definition: PFAlgo.cc:3511
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 3155 of file PFAlgo.cc.

References ALCARECOTkAlJpsiMuMu_cff::charge, dptRel_DispVtx_, hcalRecHitTable_cff::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_2023v12_cff::track.

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

3155  {
3156  const auto* eltTrack = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
3157 
3158  const reco::TrackRef& trackRef = eltTrack->trackRef();
3159  const reco::Track& track = *trackRef;
3160  const reco::MuonRef& muonRef = eltTrack->muonRef();
3161  int charge = track.charge() > 0 ? 1 : -1;
3162 
3163  // Assume this particle is a charged Hadron
3164  double px = track.px();
3165  double py = track.py();
3166  double pz = track.pz();
3167  double energy = sqrt(track.p() * track.p() + 0.13957 * 0.13957);
3168 
3169  LogTrace("PFAlgo|reconstructTrack") << "Reconstructing PF candidate from track of pT = " << track.pt()
3170  << " eta = " << track.eta() << " phi = " << track.phi() << " px = " << px
3171  << " py = " << py << " pz = " << pz << " energy = " << energy;
3172 
3173  // Create a PF Candidate
3174  ::math::XYZTLorentzVector momentum(px, py, pz, energy);
3176 
3177  // Add it to the stack
3178  LogTrace("PFAlgo|reconstructTrack") << "Creating PFCandidate charge=" << charge << ", type=" << particleType
3179  << ", pt=" << momentum.pt() << ", eta=" << momentum.eta()
3180  << ", phi=" << momentum.phi();
3181  pfCandidates_->push_back(PFCandidate(charge, momentum, particleType));
3182  //Set vertex and stuff like this
3183  pfCandidates_->back().setVertex(trackRef->vertex());
3184  pfCandidates_->back().setTrackRef(trackRef);
3185  pfCandidates_->back().setPositionAtECALEntrance(eltTrack->positionAtECALEntrance());
3186  if (muonRef.isNonnull())
3187  pfCandidates_->back().setMuonRef(muonRef);
3188 
3189  //Set time
3190  if (elt.isTimeValid())
3191  pfCandidates_->back().setTime(elt.time(), elt.timeError());
3192 
3193  //OK Now try to reconstruct the particle as a muon
3194  bool isMuon = pfmu_->reconstructMuon(pfCandidates_->back(), muonRef, allowLoose);
3195  bool isFromDisp = isFromSecInt(elt, "secondary");
3196 
3197  if ((!isMuon) && isFromDisp) {
3198  double dpt = trackRef->ptError();
3199  double dptRel = dpt / trackRef->pt() * 100;
3200  //If the track is ill measured it is better to not refit it, since the track information probably would not be used.
3201  //In the PFAlgo we use the trackref information. If the track error is too big the refitted information might be very different
3202  // from the not refitted one.
3203  if (dptRel < dptRel_DispVtx_) {
3204  LogTrace("PFAlgo|reconstructTrack")
3205  << "Not refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy;
3206  //reco::TrackRef trackRef = eltTrack->trackRef();
3208  eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef();
3209  reco::Track trackRefit = vRef->refittedTrack(trackRef);
3210  //change the momentum with the refitted track
3211  ::math::XYZTLorentzVector momentum(
3212  trackRefit.px(), trackRefit.py(), trackRefit.pz(), sqrt(trackRefit.p() * trackRefit.p() + 0.13957 * 0.13957));
3213  LogTrace("PFAlgo|reconstructTrack")
3214  << "Refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy;
3215  }
3216  pfCandidates_->back().setFlag(reco::PFCandidate::T_FROM_DISP, true);
3217  pfCandidates_->back().setDisplacedVertexRef(
3218  eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef(),
3220  }
3221 
3222  // 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
3223  if (isFromSecInt(elt, "primary") && !isMuon) {
3224  pfCandidates_->back().setFlag(reco::PFCandidate::T_TO_DISP, true);
3225  pfCandidates_->back().setDisplacedVertexRef(
3226  eltTrack->displacedVertexRef(reco::PFBlockElement::T_TO_DISP)->displacedVertexRef(),
3228  }
3229 
3230  // returns index to the newly created PFCandidate
3231  return pfCandidates_->size() - 1;
3232 }
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:3482
#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, electrons_cff::ecalEnergy, bookConverter::elements, reco::PFCandidate::elementsInBlocks(), PFEnergyCalibration::energyEmHad(), PVValHelper::eta, HLT_2023v12_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_, pfClustersFromCombinedCaloHF_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:3234
const Point & position() const
position
Definition: Vertex.h:128
bool isFromSecInt(const reco::PFBlockElement &eTrack, std::string order) const
Definition: PFAlgo.cc:3482
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:665
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:3439
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:3155
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float > > XYZPointF
point in space with cartesian internal representation
Definition: RntStructs.h:13
double dptRel_DispVtx_
Definition: PFAlgo.h:285
double neutralHadronEnergyResolution(double clusterEnergy, double clusterEta) const
todo: use PFClusterTools for this
Definition: PFAlgo.cc:3375
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
TupleMultiplicity< TrackerTraits > const *__restrict__ uint32_t nHits

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

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

Referenced by createCandidatesHCAL(), and reconstructCluster().

3345  {
3346  std::array<double, 7> energyPerDepth;
3347  std::fill(energyPerDepth.begin(), energyPerDepth.end(), 0.0);
3348  for (auto& hitRefAndFrac : cluster.recHitFractions()) {
3349  const auto& hit = *hitRefAndFrac.recHitRef();
3350  if (DetId(hit.detId()).det() == DetId::Hcal) {
3351  if (hit.depth() == 0) {
3352  edm::LogWarning("setHcalDepthInfo") << "Depth zero found";
3353  continue;
3354  }
3355  if (hit.depth() < 1 || hit.depth() > 7) {
3356  throw cms::Exception("CorruptData") << "Bogus depth " << hit.depth() << " at detid " << hit.detId() << "\n";
3357  }
3358  energyPerDepth[hit.depth() - 1] += hitRefAndFrac.fraction() * hit.energy();
3359  }
3360  }
3361  double sum = std::accumulate(energyPerDepth.begin(), energyPerDepth.end(), 0.);
3362  std::array<float, 7> depthFractions;
3363  if (sum > 0) {
3364  for (unsigned int i = 0; i < depthFractions.size(); ++i) {
3365  depthFractions[i] = energyPerDepth[i] / sum;
3366  }
3367  } else {
3368  std::fill(depthFractions.begin(), depthFractions.end(), 0.f);
3369  }
3370  cand.setHcalDepthEnergyFractions(depthFractions);
3371 }
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 DiMuonV_cfg::muons.

Referenced by PFProducer::produce().

66 { muonHandle_ = muons; }
muons
the two sets of parameters below are mutually exclusive, depending if RECO or ALCARECO is used the us...
Definition: DiMuonV_cfg.py:212
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_2023v12_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_2023v12_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:307
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().