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:3406
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:3411
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:3387
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 3441 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().

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

◆ checkCleaning()

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

Check HF Cleaning.

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

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

◆ checkGoodTrackDeadHcal()

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

Definition at line 902 of file PFAlgo.cc.

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

Referenced by elementLoop().

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

◆ checkHasDeadHcal()

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

Definition at line 875 of file PFAlgo.cc.

Referenced by elementLoop().

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

◆ conversionAlgo()

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

Definition at line 370 of file PFAlgo.cc.

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

Referenced by processBlock().

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

◆ createCandidatesECAL()

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

Definition at line 3028 of file PFAlgo.cc.

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

Referenced by processBlock().

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

◆ createCandidatesHCAL()

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

Definition at line 1702 of file PFAlgo.cc.

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

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

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

Referenced by processBlock().

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

◆ createCandidatesHF()

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

Definition at line 1259 of file PFAlgo.cc.

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

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

◆ decideType()

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

Definition at line 1201 of file PFAlgo.cc.

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

Referenced by elementLoop().

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

◆ egammaFilters()

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

Definition at line 220 of file PFAlgo.cc.

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

Referenced by processBlock().

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

◆ elementLoop()

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

Definition at line 997 of file PFAlgo.cc.

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

Referenced by processBlock().

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

◆ getCleanedCandidates()

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

Definition at line 112 of file PFAlgo.h.

References pfCleanedCandidates_.

Referenced by PFProducer::produce().

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

◆ getPFMuonAlgo()

PFMuonAlgo * PFAlgo::getPFMuonAlgo ( )

Definition at line 66 of file PFAlgo.cc.

References pfmu_.

Referenced by PFProducer::produce().

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

◆ hfEnergyResolution()

double PFAlgo::hfEnergyResolution ( double  clusterEnergy) const
private

Definition at line 3394 of file PFAlgo.cc.

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

Referenced by createCandidatesHF().

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

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

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

Referenced by createCandidatesHCAL(), and recoTracksNotHCAL().

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

◆ nSigmaHCAL()

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

Definition at line 3387 of file PFAlgo.cc.

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

Referenced by createCandidatesHCAL().

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

References JetChargeProducer_cfi::exp, nSigmaEConstHFEM, and nSigmaHFEM_.

Referenced by createCandidatesHF().

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

References JetChargeProducer_cfi::exp, nSigmaEConstHFHAD, and nSigmaHFHAD_.

Referenced by createCandidatesHF().

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

◆ postCleaning()

void PFAlgo::postCleaning ( )
private

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

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

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

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

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

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

3157  {
3158  const auto* eltTrack = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
3159 
3160  const reco::TrackRef& trackRef = eltTrack->trackRef();
3161  const reco::Track& track = *trackRef;
3162  const reco::MuonRef& muonRef = eltTrack->muonRef();
3163  int charge = track.charge() > 0 ? 1 : -1;
3164 
3165  // Assume this particle is a charged Hadron
3166  double px = track.px();
3167  double py = track.py();
3168  double pz = track.pz();
3169  double energy = sqrt(track.p() * track.p() + 0.13957 * 0.13957);
3170 
3171  LogTrace("PFAlgo|reconstructTrack") << "Reconstructing PF candidate from track of pT = " << track.pt()
3172  << " eta = " << track.eta() << " phi = " << track.phi() << " px = " << px
3173  << " py = " << py << " pz = " << pz << " energy = " << energy;
3174 
3175  // Create a PF Candidate
3176  ::math::XYZTLorentzVector momentum(px, py, pz, energy);
3178 
3179  // Add it to the stack
3180  LogTrace("PFAlgo|reconstructTrack") << "Creating PFCandidate charge=" << charge << ", type=" << particleType
3181  << ", pt=" << momentum.pt() << ", eta=" << momentum.eta()
3182  << ", phi=" << momentum.phi();
3183  pfCandidates_->push_back(PFCandidate(charge, momentum, particleType));
3184  //Set vertex and stuff like this
3185  pfCandidates_->back().setVertex(trackRef->vertex());
3186  pfCandidates_->back().setTrackRef(trackRef);
3187  pfCandidates_->back().setPositionAtECALEntrance(eltTrack->positionAtECALEntrance());
3188  if (muonRef.isNonnull())
3189  pfCandidates_->back().setMuonRef(muonRef);
3190 
3191  //Set time
3192  if (elt.isTimeValid())
3193  pfCandidates_->back().setTime(elt.time(), elt.timeError());
3194 
3195  //OK Now try to reconstruct the particle as a muon
3196  bool isMuon = pfmu_->reconstructMuon(pfCandidates_->back(), muonRef, allowLoose);
3197  bool isFromDisp = isFromSecInt(elt, "secondary");
3198 
3199  if ((!isMuon) && isFromDisp) {
3200  double dpt = trackRef->ptError();
3201  double dptRel = dpt / trackRef->pt() * 100;
3202  //If the track is ill measured it is better to not refit it, since the track information probably would not be used.
3203  //In the PFAlgo we use the trackref information. If the track error is too big the refitted information might be very different
3204  // from the not refitted one.
3205  if (dptRel < dptRel_DispVtx_) {
3206  LogTrace("PFAlgo|reconstructTrack")
3207  << "Not refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy;
3208  //reco::TrackRef trackRef = eltTrack->trackRef();
3210  eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef();
3211  reco::Track trackRefit = vRef->refittedTrack(trackRef);
3212  //change the momentum with the refitted track
3213  ::math::XYZTLorentzVector momentum(
3214  trackRefit.px(), trackRefit.py(), trackRefit.pz(), sqrt(trackRefit.p() * trackRefit.p() + 0.13957 * 0.13957));
3215  LogTrace("PFAlgo|reconstructTrack")
3216  << "Refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy;
3217  }
3218  pfCandidates_->back().setFlag(reco::PFCandidate::T_FROM_DISP, true);
3219  pfCandidates_->back().setDisplacedVertexRef(
3220  eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef(),
3222  }
3223 
3224  // 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
3225  if (isFromSecInt(elt, "primary") && !isMuon) {
3226  pfCandidates_->back().setFlag(reco::PFCandidate::T_TO_DISP, true);
3227  pfCandidates_->back().setDisplacedVertexRef(
3228  eltTrack->displacedVertexRef(reco::PFBlockElement::T_TO_DISP)->displacedVertexRef(),
3230  }
3231 
3232  // returns index to the newly created PFCandidate
3233  return pfCandidates_->size() - 1;
3234 }
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:3484
#define LogTrace(id)
XYZTLorentzVectorD XYZTLorentzVector
Lorentz vector with cylindrical internal representation using pseudorapidity.
Definition: LorentzVector.h:29
T sqrt(T t)
Definition: SSEVec.h:19
bool isTimeValid() const
do we have a valid time information
double dptRel_DispVtx_
Definition: PFAlgo.h:285
double pz() const
z coordinate of momentum vector
Definition: TrackBase.h:646
std::unique_ptr< PFMuonAlgo > pfmu_
Definition: PFAlgo.h:262

◆ recoTracksNotHCAL()

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

Definition at line 399 of file PFAlgo.cc.

References funct::abs(), cms::cuda::assert(), associatePSClusters(), groupFilesInBlocks::block, calibration_, dptRel_DispVtx_, ECAL, reco::PFBlockElement::ECAL, bookConverter::elements, reco::PFCandidate::elementsInBlocks(), PFEnergyCalibration::energyEmHad(), PVValHelper::eta, HLT_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:3236
const Point & position() const
position
Definition: Vertex.h:127
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float > > XYZPointF
point in space with cartesian internal representation
Definition: Point3D.h:10
bool isFromSecInt(const reco::PFBlockElement &eTrack, std::string order) const
Definition: PFAlgo.cc:3484
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:3441
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:3157
double dptRel_DispVtx_
Definition: PFAlgo.h:285
double neutralHadronEnergyResolution(double clusterEnergy, double clusterEta) const
todo: use PFClusterTools for this
Definition: PFAlgo.cc:3377
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 ntt = 1;
948  unsigned index = ecalElems.begin()->second;
949  std::multimap<double, unsigned> sortedTracks;
950  block.associatedElements(index, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
951  LogTrace("PFAlgo|elementLoop") << "The closest ECAL cluster is linked to " << sortedTracks.size()
952  << " tracks, with distance = " << ecalElems.begin()->first;
953 
954  LogTrace("PFAlgo|elementLoop") << "Looping over sortedTracks";
955  // Loop over all tracks
956  for (auto const& trk : sortedTracks) {
957  unsigned jTrack = trk.second;
958  LogTrace("PFAlgo|elementLoop") << "jTrack=" << jTrack;
959  // Track must be active
960  if (!active[jTrack])
961  continue;
962  LogTrace("PFAlgo|elementLoop") << "active[jTrack]=" << active[jTrack];
963 
964  // The loop is on the other tracks !
965  if (jTrack == iTrack)
966  continue;
967  LogTrace("PFAlgo|elementLoop") << "skipping jTrack=" << jTrack << " for same iTrack";
968 
969  // Check if the ECAL closest to this track is the current ECAL
970  // Otherwise ignore this track in the neutral energy determination
971  std::multimap<double, unsigned> sortedECAL;
972  block.associatedElements(jTrack, linkData, sortedECAL, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
973  if (sortedECAL.begin()->second != index)
974  continue;
975  LogTrace("PFAlgo|elementLoop") << " track " << jTrack << " with closest ECAL identical ";
976 
977  // Check if this track is also linked to an HCAL
978  std::multimap<double, unsigned> sortedHCAL;
979  block.associatedElements(jTrack, linkData, sortedHCAL, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
980  if (sortedHCAL.empty())
981  continue;
982  LogTrace("PFAlgo|elementLoop") << " and with an HCAL cluster " << sortedHCAL.begin()->second;
983  ntt++;
984 
985  // In that case establish a link with the first track
986  block.setLink(iTrack, sortedHCAL.begin()->second, sortedECAL.begin()->first, linkData, PFBlock::LINKTEST_RECHIT);
987 
988  } // End other tracks
989 
990  // Redefine HCAL elements
991  block.associatedElements(iTrack, linkData, hcalElems, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
992 
993  if (!hcalElems.empty())
994  LogTrace("PFAlgo|elementLoop") << "Track linked back to HCAL due to ECAL sharing with other tracks";
995 }
#define LogTrace(id)

◆ setCandConnectorParameters() [1/2]

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

Definition at line 68 of file PFAlgo.h.

References connector_, HLT_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 3347 of file PFAlgo.cc.

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

Referenced by createCandidatesHCAL(), and reconstructCluster().

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

◆ setHOTag()

void PFAlgo::setHOTag ( bool  ho)
inline

Definition at line 65 of file PFAlgo.h.

References photonIsolationHIProducer_cfi::ho, and useHO_.

◆ setMuonHandle()

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

Definition at line 66 of file PFAlgo.h.

References muonHandle_, and PDWG_BPHSkim_cff::muons.

Referenced by PFProducer::produce().

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

◆ setPFVertexParameters()

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

Definition at line 110 of file PFAlgo.cc.

References nVtx_, pfmu_, primaryVertex_, HLT_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:303
double maxDeltaPhiPt_
Definition: PFAlgo.h:323
double maxSignificance_
Definition: PFAlgo.h:321
double minHFCleaningPt_
Definition: PFAlgo.h:319
double minDeltaMet_
Definition: PFAlgo.h:324
double minSignificance_
Definition: PFAlgo.h:320
double minSignificanceReduction_
Definition: PFAlgo.h:322
bool postHFCleaning_
Definition: PFAlgo.h:317

Friends And Related Function Documentation

◆ operator<<

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

Member Data Documentation

◆ calibration_

PFEnergyCalibration& PFAlgo::calibration_
private

Definition at line 257 of file PFAlgo.h.

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

◆ connector_

PFCandConnector PFAlgo::connector_
private

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

Definition at line 290 of file PFAlgo.h.

Referenced by makeConnectedCandidates(), and setCandConnectorParameters().

◆ dptRel_DispVtx_

double PFAlgo::dptRel_DispVtx_
private

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

Definition at line 285 of file PFAlgo.h.

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

◆ factors45_

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

Definition at line 298 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), and PFAlgo().

◆ goodPixelTrackDeadHcal_chi2n_

float PFAlgo::goodPixelTrackDeadHcal_chi2n_
private

Definition at line 310 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_dxy_

float PFAlgo::goodPixelTrackDeadHcal_dxy_
private

Definition at line 313 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_dz_

float PFAlgo::goodPixelTrackDeadHcal_dz_
private

Definition at line 314 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_maxLost3Hit_

int PFAlgo::goodPixelTrackDeadHcal_maxLost3Hit_
private

Definition at line 311 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_maxLost4Hit_

int PFAlgo::goodPixelTrackDeadHcal_maxLost4Hit_
private

Definition at line 312 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_maxPt_

float PFAlgo::goodPixelTrackDeadHcal_maxPt_
private

Definition at line 308 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_minEta_

float PFAlgo::goodPixelTrackDeadHcal_minEta_
private

Definition at line 307 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodPixelTrackDeadHcal_ptErrRel_

float PFAlgo::goodPixelTrackDeadHcal_ptErrRel_
private

Definition at line 309 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodTrackDeadHcal_chi2n_

float PFAlgo::goodTrackDeadHcal_chi2n_
private

Definition at line 302 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodTrackDeadHcal_dxy_

float PFAlgo::goodTrackDeadHcal_dxy_
private

Definition at line 305 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodTrackDeadHcal_layers_

int PFAlgo::goodTrackDeadHcal_layers_
private

Definition at line 303 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodTrackDeadHcal_ptErrRel_

float PFAlgo::goodTrackDeadHcal_ptErrRel_
private

Variables for track cleaning in bad HCal areas.

Definition at line 301 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ goodTrackDeadHcal_validFr_

float PFAlgo::goodTrackDeadHcal_validFr_
private

Definition at line 304 of file PFAlgo.h.

Referenced by checkGoodTrackDeadHcal(), and PFAlgo().

◆ maxDeltaPhiPt_

double PFAlgo::maxDeltaPhiPt_
private

Definition at line 323 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

◆ maxSignificance_

double PFAlgo::maxSignificance_
private

Definition at line 321 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

◆ minDeltaMet_

double PFAlgo::minDeltaMet_
private

Definition at line 324 of file PFAlgo.h.

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

◆ minHFCleaningPt_

double PFAlgo::minHFCleaningPt_
private

Definition at line 319 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

◆ minSignificance_

double PFAlgo::minSignificance_
private

Definition at line 320 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

◆ minSignificanceReduction_

double PFAlgo::minSignificanceReduction_
private

Definition at line 322 of file PFAlgo.h.

Referenced by postCleaning(), and setPostHFCleaningParameters().

◆ muonECAL_

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

Definition at line 294 of file PFAlgo.h.

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

◆ muonHandle_

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

Definition at line 331 of file PFAlgo.h.

Referenced by reconstructParticles(), and setMuonHandle().

◆ muonHCAL_

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

Variables for muons and fakes.

Definition at line 293 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), and PFAlgo().

◆ muonHO_

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

Definition at line 295 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), and PFAlgo().

◆ nSigmaECAL_

const double PFAlgo::nSigmaECAL_
private

number of sigma to judge energy excess in ECAL

Definition at line 245 of file PFAlgo.h.

Referenced by recoTracksNotHCAL().

◆ nSigmaEConstHCAL

const double PFAlgo::nSigmaEConstHCAL = 100.
private

Definition at line 334 of file PFAlgo.h.

Referenced by nSigmaHCAL().

◆ nSigmaEConstHFEM

const double PFAlgo::nSigmaEConstHFEM = 100.
private

Definition at line 335 of file PFAlgo.h.

Referenced by nSigmaHFEM().

◆ nSigmaEConstHFHAD

const double PFAlgo::nSigmaEConstHFHAD = 100.
private

Definition at line 336 of file PFAlgo.h.

Referenced by nSigmaHFHAD().

◆ nSigmaHCAL_

const double PFAlgo::nSigmaHCAL_
private

number of sigma to judge energy excess in HCAL

Definition at line 248 of file PFAlgo.h.

Referenced by nSigmaHCAL().

◆ nSigmaHFEM_

const double PFAlgo::nSigmaHFEM_
private

number of sigma to judge energy excess in HF

Definition at line 251 of file PFAlgo.h.

Referenced by nSigmaHFEM().

◆ nSigmaHFHAD_

const double PFAlgo::nSigmaHFHAD_
private

Definition at line 252 of file PFAlgo.h.

Referenced by nSigmaHFHAD().

◆ nSigmaTRACK_

double PFAlgo::nSigmaTRACK_
private

Definition at line 296 of file PFAlgo.h.

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

◆ nVtx_

int PFAlgo::nVtx_
private

Definition at line 286 of file PFAlgo.h.

Referenced by egammaFilters(), and setPFVertexParameters().

◆ pfCandidates_

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

◆ pfCleanedCandidates_

reco::PFCandidateCollection PFAlgo::pfCleanedCandidates_
private

Definition at line 228 of file PFAlgo.h.

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

◆ pfEgammaCandidates_

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

Definition at line 267 of file PFAlgo.h.

Referenced by egammaFilters(), and setEGammaCollections().

◆ pfmu_

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

◆ postHFCleaning_

bool PFAlgo::postHFCleaning_
private

Definition at line 317 of file PFAlgo.h.

Referenced by reconstructParticles(), and setPostHFCleaningParameters().

◆ postMuonCleaning_

bool PFAlgo::postMuonCleaning_
private

Definition at line 318 of file PFAlgo.h.

◆ primaryVertex_

reco::Vertex PFAlgo::primaryVertex_
private

◆ ptError_

double PFAlgo::ptError_
private

Definition at line 297 of file PFAlgo.h.

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

◆ rejectTracks_Bad_

bool PFAlgo::rejectTracks_Bad_
private

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

Definition at line 276 of file PFAlgo.h.

Referenced by createCandidatesHCAL(), and setDisplacedVerticesParameters().

◆ rejectTracks_Step45_

bool PFAlgo::rejectTracks_Step45_
private

◆ resolHF_square_

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

Definition at line 255 of file PFAlgo.h.

Referenced by hfEnergyResolution(), and PFAlgo().

◆ thepfEnergyCalibrationHF_

PFEnergyCalibrationHF& PFAlgo::thepfEnergyCalibrationHF_
private

Definition at line 258 of file PFAlgo.h.

Referenced by createCandidatesHF().

◆ useBestMuonTrack_

double PFAlgo::useBestMuonTrack_
private

Definition at line 325 of file PFAlgo.h.

◆ useEGammaFilters_

bool PFAlgo::useEGammaFilters_
private

Variables for NEW EGAMMA selection.

Definition at line 265 of file PFAlgo.h.

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

◆ useHO_

bool PFAlgo::useHO_
private

Definition at line 260 of file PFAlgo.h.

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

◆ usePFConversions_

bool PFAlgo::usePFConversions_
private

Definition at line 280 of file PFAlgo.h.

Referenced by processBlock(), and setDisplacedVerticesParameters().

◆ usePFDecays_

bool PFAlgo::usePFDecays_
private

Definition at line 281 of file PFAlgo.h.

Referenced by isFromSecInt(), and setDisplacedVerticesParameters().

◆ usePFMuonMomAssign_

bool PFAlgo::usePFMuonMomAssign_
private

Definition at line 272 of file PFAlgo.h.

◆ usePFNuclearInteractions_

bool PFAlgo::usePFNuclearInteractions_
private

Definition at line 279 of file PFAlgo.h.

Referenced by isFromSecInt(), and setDisplacedVerticesParameters().

◆ useProtectionsForJetMET_

bool PFAlgo::useProtectionsForJetMET_
private

Definition at line 266 of file PFAlgo.h.

Referenced by egammaFilters(), and setEGammaParameters().

◆ useVertices_

bool PFAlgo::useVertices_ = false
private

Definition at line 329 of file PFAlgo.h.

Referenced by reconstructCluster(), and setPFVertexParameters().

◆ valueMapGedElectrons_

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

Definition at line 268 of file PFAlgo.h.

Referenced by setEGammaCollections().

◆ valueMapGedPhotons_

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

Definition at line 269 of file PFAlgo.h.

Referenced by setEGammaCollections().