File indexing completed on 2023-01-21 00:19:49
0001 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0002 #include "RecoParticleFlow/PFProducer/interface/PFAlgo.h"
0003 #include "RecoParticleFlow/PFProducer/interface/PFMuonAlgo.h"
0004 #include "RecoParticleFlow/PFProducer/interface/PFElectronExtraEqual.h"
0005 #include "RecoParticleFlow/PFTracking/interface/PFTrackAlgoTools.h"
0006
0007 #include "TDecompChol.h"
0008
0009 #include <numeric>
0010 #include <fstream>
0011
0012 using namespace std;
0013 using namespace reco;
0014
0015 PFAlgo::PFAlgo(double nSigmaECAL,
0016 double nSigmaHCAL,
0017 double nSigmaHFEM,
0018 double nSigmaHFHAD,
0019 std::vector<double> resolHF_square,
0020 PFEnergyCalibration& calibration,
0021 PFEnergyCalibrationHF& thepfEnergyCalibrationHF,
0022 const edm::ParameterSet& pset)
0023 : pfCandidates_(new PFCandidateCollection),
0024 nSigmaECAL_(nSigmaECAL),
0025 nSigmaHCAL_(nSigmaHCAL),
0026 nSigmaHFEM_(nSigmaHFEM),
0027 nSigmaHFHAD_(nSigmaHFHAD),
0028 resolHF_square_(resolHF_square),
0029 calibration_(calibration),
0030 thepfEnergyCalibrationHF_(thepfEnergyCalibrationHF),
0031 connector_() {
0032 const edm::ParameterSet pfMuonAlgoParams = pset.getParameter<edm::ParameterSet>("PFMuonAlgoParameters");
0033 bool postMuonCleaning = pset.getParameter<bool>("postMuonCleaning");
0034 pfmu_ = std::make_unique<PFMuonAlgo>(pfMuonAlgoParams, postMuonCleaning);
0035
0036
0037 assert(resolHF_square_.size() == 3);
0038
0039
0040 muonHCAL_ = pset.getParameter<std::vector<double>>("muon_HCAL");
0041 muonECAL_ = pset.getParameter<std::vector<double>>("muon_ECAL");
0042 muonHO_ = pset.getParameter<std::vector<double>>("muon_HO");
0043 assert(muonHCAL_.size() == 2 && muonECAL_.size() == 2 && muonHO_.size() == 2);
0044 nSigmaTRACK_ = pset.getParameter<double>("nsigma_TRACK");
0045 ptError_ = pset.getParameter<double>("pt_Error");
0046 factors45_ = pset.getParameter<std::vector<double>>("factors_45");
0047 assert(factors45_.size() == 2);
0048
0049
0050 goodTrackDeadHcal_ptErrRel_ = pset.getParameter<double>("goodTrackDeadHcal_ptErrRel");
0051 goodTrackDeadHcal_chi2n_ = pset.getParameter<double>("goodTrackDeadHcal_chi2n");
0052 goodTrackDeadHcal_layers_ = pset.getParameter<uint32_t>("goodTrackDeadHcal_layers");
0053 goodTrackDeadHcal_validFr_ = pset.getParameter<double>("goodTrackDeadHcal_validFr");
0054 goodTrackDeadHcal_dxy_ = pset.getParameter<double>("goodTrackDeadHcal_dxy");
0055
0056 goodPixelTrackDeadHcal_minEta_ = pset.getParameter<double>("goodPixelTrackDeadHcal_minEta");
0057 goodPixelTrackDeadHcal_maxPt_ = pset.getParameter<double>("goodPixelTrackDeadHcal_maxPt");
0058 goodPixelTrackDeadHcal_ptErrRel_ = pset.getParameter<double>("goodPixelTrackDeadHcal_ptErrRel");
0059 goodPixelTrackDeadHcal_chi2n_ = pset.getParameter<double>("goodPixelTrackDeadHcal_chi2n");
0060 goodPixelTrackDeadHcal_maxLost3Hit_ = pset.getParameter<int32_t>("goodPixelTrackDeadHcal_maxLost3Hit");
0061 goodPixelTrackDeadHcal_maxLost4Hit_ = pset.getParameter<int32_t>("goodPixelTrackDeadHcal_maxLost4Hit");
0062 goodPixelTrackDeadHcal_dxy_ = pset.getParameter<double>("goodPixelTrackDeadHcal_dxy");
0063 goodPixelTrackDeadHcal_dz_ = pset.getParameter<double>("goodPixelTrackDeadHcal_dz");
0064 }
0065
0066 PFMuonAlgo* PFAlgo::getPFMuonAlgo() { return pfmu_.get(); }
0067
0068 void PFAlgo::setEGammaParameters(bool use_EGammaFilters, bool useProtectionsForJetMET) {
0069 useEGammaFilters_ = use_EGammaFilters;
0070
0071 if (!useEGammaFilters_)
0072 return;
0073
0074 useProtectionsForJetMET_ = useProtectionsForJetMET;
0075 }
0076 void PFAlgo::setEGammaCollections(const edm::View<reco::PFCandidate>& pfEgammaCandidates,
0077 const edm::ValueMap<reco::GsfElectronRef>& valueMapGedElectrons,
0078 const edm::ValueMap<reco::PhotonRef>& valueMapGedPhotons) {
0079 if (useEGammaFilters_) {
0080 pfEgammaCandidates_ = &pfEgammaCandidates;
0081 valueMapGedElectrons_ = &valueMapGedElectrons;
0082 valueMapGedPhotons_ = &valueMapGedPhotons;
0083 }
0084 }
0085
0086 void PFAlgo::setPostHFCleaningParameters(bool postHFCleaning, const edm::ParameterSet& pfHFCleaningParams) {
0087 postHFCleaning_ = postHFCleaning;
0088 minHFCleaningPt_ = pfHFCleaningParams.getParameter<double>("minHFCleaningPt");
0089 minSignificance_ = pfHFCleaningParams.getParameter<double>("minSignificance");
0090 maxSignificance_ = pfHFCleaningParams.getParameter<double>("maxSignificance");
0091 minSignificanceReduction_ = pfHFCleaningParams.getParameter<double>("minSignificanceReduction");
0092 maxDeltaPhiPt_ = pfHFCleaningParams.getParameter<double>("maxDeltaPhiPt");
0093 minDeltaMet_ = pfHFCleaningParams.getParameter<double>("minDeltaMet");
0094 }
0095
0096 void PFAlgo::setDisplacedVerticesParameters(bool rejectTracks_Bad,
0097 bool rejectTracks_Step45,
0098 bool usePFNuclearInteractions,
0099 bool usePFConversions,
0100 bool usePFDecays,
0101 double dptRel_DispVtx) {
0102 rejectTracks_Bad_ = rejectTracks_Bad;
0103 rejectTracks_Step45_ = rejectTracks_Step45;
0104 usePFNuclearInteractions_ = usePFNuclearInteractions;
0105 usePFConversions_ = usePFConversions;
0106 usePFDecays_ = usePFDecays;
0107 dptRel_DispVtx_ = dptRel_DispVtx;
0108 }
0109
0110 void PFAlgo::setPFVertexParameters(bool useVertex, reco::VertexCollection const& primaryVertices) {
0111 useVertices_ = useVertex;
0112
0113
0114 pfmu_->setInputsForCleaning(primaryVertices);
0115
0116
0117 bool primaryVertexFound = false;
0118 nVtx_ = primaryVertices.size();
0119 for (auto const& vertex : primaryVertices) {
0120 if (vertex.isValid() && (!vertex.isFake())) {
0121 primaryVertex_ = vertex;
0122 primaryVertexFound = true;
0123 break;
0124 }
0125 }
0126
0127 useVertices_ = useVertex && primaryVertexFound;
0128 }
0129
0130 void PFAlgo::reconstructParticles(const reco::PFBlockHandle& blockHandle, PFEGammaFilters const* pfegamma) {
0131 auto const& blocks = *blockHandle;
0132
0133
0134 pfCandidates_->clear();
0135
0136 LogTrace("PFAlgo|reconstructParticles")
0137 << "start of function PFAlgo::reconstructParticles, blocks.size()=" << blocks.size();
0138
0139
0140 std::list<reco::PFBlockRef> hcalBlockRefs;
0141 std::list<reco::PFBlockRef> ecalBlockRefs;
0142 std::list<reco::PFBlockRef> hoBlockRefs;
0143 std::list<reco::PFBlockRef> otherBlockRefs;
0144
0145 for (unsigned i = 0; i < blocks.size(); ++i) {
0146 reco::PFBlockRef blockref = reco::PFBlockRef(blockHandle, i);
0147
0148 const reco::PFBlock& block = *blockref;
0149 const edm::OwnVector<reco::PFBlockElement>& elements = block.elements();
0150
0151 bool singleEcalOrHcal = false;
0152 if (elements.size() == 1) {
0153 if (elements[0].type() == reco::PFBlockElement::ECAL) {
0154 ecalBlockRefs.push_back(blockref);
0155 singleEcalOrHcal = true;
0156 }
0157 if (elements[0].type() == reco::PFBlockElement::HCAL) {
0158 hcalBlockRefs.push_back(blockref);
0159 singleEcalOrHcal = true;
0160 }
0161 if (elements[0].type() == reco::PFBlockElement::HO) {
0162
0163 hoBlockRefs.push_back(blockref);
0164 singleEcalOrHcal = true;
0165 }
0166 }
0167
0168 if (!singleEcalOrHcal) {
0169 otherBlockRefs.push_back(blockref);
0170 }
0171 }
0172
0173 LogTrace("PFAlgo|reconstructParticles")
0174 << "# Ecal blocks: " << ecalBlockRefs.size() << ", # Hcal blocks: " << hcalBlockRefs.size()
0175 << ", # HO blocks: " << hoBlockRefs.size() << ", # Other blocks: " << otherBlockRefs.size();
0176
0177
0178
0179
0180 unsigned nblcks = 0;
0181 for (auto const& other : otherBlockRefs) {
0182 LogTrace("PFAlgo|reconstructParticles") << "processBlock, Block number " << nblcks++;
0183 processBlock(other, hcalBlockRefs, ecalBlockRefs, pfegamma);
0184 }
0185
0186 std::list<reco::PFBlockRef> empty;
0187
0188 unsigned hblcks = 0;
0189
0190 for (auto const& hcal : hcalBlockRefs) {
0191 LogTrace("PFAlgo|reconstructParticles") << "processBlock, HCAL block number " << hblcks++;
0192 processBlock(hcal, empty, empty, pfegamma);
0193 }
0194
0195 unsigned eblcks = 0;
0196
0197 for (auto const& ecal : ecalBlockRefs) {
0198 LogTrace("PFAlgo|reconstructParticles") << "processBlock, ECAL block number " << eblcks++;
0199 processBlock(ecal, empty, empty, pfegamma);
0200 }
0201
0202
0203 pfCleanedCandidates_.clear();
0204
0205 if (postHFCleaning_) {
0206 postCleaning();
0207 }
0208
0209
0210 pfmu_->postClean(pfCandidates_.get());
0211
0212
0213 if (muonHandle_.isValid())
0214 pfmu_->addMissingMuons(muonHandle_, pfCandidates_.get());
0215
0216 LogTrace("PFAlgo|reconstructParticles")
0217 << "end of function PFAlgo::reconstructParticles, pfCandidates_->size()=" << pfCandidates_->size();
0218 }
0219
0220 void PFAlgo::egammaFilters(const reco::PFBlockRef& blockref,
0221 std::vector<bool>& active,
0222 PFEGammaFilters const* pfegamma) {
0223
0224
0225 unsigned int negmcandidates = pfEgammaCandidates_->size();
0226 LogTrace("PFAlgo|egammaFilters") << "start of function PFAlgo::egammaFilters(), negmcandidates=" << negmcandidates;
0227
0228 for (unsigned int ieg = 0; ieg < negmcandidates; ++ieg) {
0229
0230 reco::PFCandidateRef pfEgmRef = pfEgammaCandidates_->refAt(ieg).castTo<reco::PFCandidateRef>();
0231
0232 const PFCandidate::ElementsInBlocks& theElements = (*pfEgmRef).elementsInBlocks();
0233 PFCandidate::ElementsInBlocks::const_iterator iegfirst = theElements.begin();
0234 bool sameBlock = false;
0235 bool isGoodElectron = false;
0236 bool isGoodPhoton = false;
0237 bool isPrimaryElectron = false;
0238 if (iegfirst->first == blockref)
0239 sameBlock = true;
0240 if (sameBlock) {
0241 LogTrace("PFAlgo|egammaFilters") << " I am in looping on EGamma Candidates: pt " << (*pfEgmRef).pt()
0242 << " eta,phi " << (*pfEgmRef).eta() << ", " << (*pfEgmRef).phi() << " charge "
0243 << (*pfEgmRef).charge();
0244
0245 if ((*pfEgmRef).gsfTrackRef().isNonnull()) {
0246 reco::GsfElectronRef gedEleRef = (*valueMapGedElectrons_)[pfEgmRef];
0247 if (gedEleRef.isNonnull()) {
0248 isGoodElectron = pfegamma->passElectronSelection(*gedEleRef, *pfEgmRef, nVtx_);
0249 isPrimaryElectron = pfegamma->isElectron(*gedEleRef);
0250 if (isGoodElectron)
0251 LogTrace("PFAlgo|egammaFilters")
0252 << "** Good Electron, pt " << gedEleRef->pt() << " eta, phi " << gedEleRef->eta() << ", "
0253 << gedEleRef->phi() << " charge " << gedEleRef->charge() << " isPrimary " << isPrimaryElectron
0254 << " isoDr03 "
0255 << (gedEleRef->dr03TkSumPt() + gedEleRef->dr03EcalRecHitSumEt() + gedEleRef->dr03HcalTowerSumEt())
0256 << " mva_isolated " << gedEleRef->mva_Isolated() << " mva_e_pi " << gedEleRef->mva_e_pi();
0257 }
0258 }
0259 if ((*pfEgmRef).superClusterRef().isNonnull()) {
0260 reco::PhotonRef gedPhoRef = (*valueMapGedPhotons_)[pfEgmRef];
0261 if (gedPhoRef.isNonnull()) {
0262 isGoodPhoton = pfegamma->passPhotonSelection(*gedPhoRef);
0263 if (isGoodPhoton)
0264 LogTrace("PFAlgo|egammaFilters") << "** Good Photon, pt " << gedPhoRef->pt() << " eta, phi "
0265 << gedPhoRef->eta() << ", " << gedPhoRef->phi() << endl;
0266 }
0267 }
0268 }
0269
0270 if (isGoodElectron && isGoodPhoton) {
0271 if (isPrimaryElectron)
0272 isGoodPhoton = false;
0273 else
0274 isGoodElectron = false;
0275 }
0276
0277
0278 if (isGoodElectron) {
0279 reco::GsfElectronRef gedEleRef = (*valueMapGedElectrons_)[pfEgmRef];
0280 reco::PFCandidate myPFElectron = *pfEgmRef;
0281
0282 bool lockTracks = false;
0283 bool isSafe = true;
0284 if (useProtectionsForJetMET_) {
0285 lockTracks = true;
0286 isSafe = pfegamma->isElectronSafeForJetMET(*gedEleRef, myPFElectron, primaryVertex_, lockTracks);
0287 }
0288
0289 if (isSafe) {
0290 reco::PFCandidate::ParticleType particleType = reco::PFCandidate::e;
0291 myPFElectron.setParticleType(particleType);
0292 myPFElectron.setCharge(gedEleRef->charge());
0293 myPFElectron.setP4(gedEleRef->p4());
0294 myPFElectron.set_mva_e_pi(gedEleRef->mva_e_pi());
0295 myPFElectron.set_mva_Isolated(gedEleRef->mva_Isolated());
0296
0297 myPFElectron.set_dnn_e_sigIsolated(gedEleRef->dnn_signal_Isolated());
0298 myPFElectron.set_dnn_e_sigNonIsolated(gedEleRef->dnn_signal_nonIsolated());
0299 myPFElectron.set_dnn_e_bkgNonIsolated(gedEleRef->dnn_bkg_nonIsolated());
0300 myPFElectron.set_dnn_e_bkgTau(gedEleRef->dnn_bkg_Tau());
0301 myPFElectron.set_dnn_e_bkgPhoton(gedEleRef->dnn_bkg_Photon());
0302
0303 LogTrace("PFAlgo|egammaFilters") << " PFAlgo: found an electron with NEW EGamma code ";
0304 LogTrace("PFAlgo|egammaFilters") << " myPFElectron: pt " << myPFElectron.pt() << " eta,phi "
0305 << myPFElectron.eta() << ", " << myPFElectron.phi() << " mva "
0306 << myPFElectron.mva_e_pi() << " charge " << myPFElectron.charge();
0307
0308 LogTrace("PFAlgo|egammaFilters") << " THE BLOCK " << *blockref;
0309 for (auto const& eb : theElements) {
0310 active[eb.second] = false;
0311 LogTrace("PFAlgo|egammaFilters") << " Elements used " << eb.second;
0312 }
0313
0314
0315 if (lockTracks) {
0316 const PFCandidate::ElementsInBlocks& extraTracks = myPFElectron.egammaExtraRef()->extraNonConvTracks();
0317 for (auto const& trk : extraTracks) {
0318 LogTrace("PFAlgo|egammaFilters") << " Extra locked track " << trk.second;
0319 active[trk.second] = false;
0320 }
0321 }
0322
0323 LogTrace("PFAlgo|egammaFilters") << "Creating PF electron: pt=" << myPFElectron.pt()
0324 << " eta=" << myPFElectron.eta() << " phi=" << myPFElectron.phi();
0325 pfCandidates_->push_back(myPFElectron);
0326
0327 } else {
0328 LogTrace("PFAlgo|egammaFilters") << "PFAlgo: Electron DISCARDED, NOT SAFE FOR JETMET ";
0329 }
0330 }
0331
0332 if (isGoodPhoton) {
0333 reco::PhotonRef gedPhoRef = (*valueMapGedPhotons_)[pfEgmRef];
0334 reco::PFCandidate myPFPhoton = *pfEgmRef;
0335 bool isSafe = true;
0336 if (useProtectionsForJetMET_) {
0337 isSafe = pfegamma->isPhotonSafeForJetMET(*gedPhoRef, myPFPhoton);
0338 }
0339
0340 if (isSafe) {
0341 reco::PFCandidate::ParticleType particleType = reco::PFCandidate::gamma;
0342 myPFPhoton.setParticleType(particleType);
0343 myPFPhoton.setCharge(0);
0344 myPFPhoton.set_mva_nothing_gamma(1.);
0345 ::math::XYZPoint v(primaryVertex_.x(), primaryVertex_.y(), primaryVertex_.z());
0346 myPFPhoton.setVertex(v);
0347 myPFPhoton.setP4(gedPhoRef->p4());
0348
0349 myPFPhoton.set_dnn_gamma(gedPhoRef->pfDNN());
0350 LogTrace("PFAlgo|egammaFilters") << " PFAlgo: found a photon with NEW EGamma code ";
0351 LogTrace("PFAlgo|egammaFilters") << " myPFPhoton: pt " << myPFPhoton.pt() << " eta,phi " << myPFPhoton.eta()
0352 << ", " << myPFPhoton.phi() << " charge " << myPFPhoton.charge();
0353
0354
0355 LogTrace("PFAlgo|egammaFilters") << " THE BLOCK " << *blockref;
0356 for (auto const& eb : theElements) {
0357 active[eb.second] = false;
0358 LogTrace("PFAlgo|egammaFilters") << " Elements used " << eb.second;
0359 }
0360 LogTrace("PFAlgo|egammaFilters") << "Creating PF photon: pt=" << myPFPhoton.pt() << " eta=" << myPFPhoton.eta()
0361 << " phi=" << myPFPhoton.phi();
0362 pfCandidates_->push_back(myPFPhoton);
0363
0364 }
0365 }
0366 }
0367 LogTrace("PFAlgo|egammaFilters") << "end of function PFAlgo::egammaFilters";
0368 }
0369
0370 void PFAlgo::conversionAlgo(const edm::OwnVector<reco::PFBlockElement>& elements, std::vector<bool>& active) {
0371 LogTrace("PFAlgo|conversionAlgo") << "start of function PFAlgo::conversionAlgo(), elements.size()="
0372 << elements.size();
0373 for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
0374 PFBlockElement::Type type = elements[iEle].type();
0375 if (type == PFBlockElement::TRACK) {
0376 LogTrace("PFAlgo|conversionAlgo") << "elements[" << iEle << "].type() == TRACK, active[" << iEle
0377 << "]=" << active[iEle];
0378 if (elements[iEle].trackRef()->algo() == reco::TrackBase::conversionStep) {
0379 active[iEle] = false;
0380 }
0381 if (elements[iEle].trackRef()->quality(reco::TrackBase::highPurity)) {
0382 LogTrace("PFAlgo|conversionAlgo") << "Track is high purity";
0383 continue;
0384 }
0385 const auto* trackRef = dynamic_cast<const reco::PFBlockElementTrack*>((&elements[iEle]));
0386 if (!(trackRef->trackType(reco::PFBlockElement::T_FROM_GAMMACONV))) {
0387 LogTrace("PFAlgo|conversionAlgo") << "!trackType(T_FROM_GAMMACONV)";
0388 continue;
0389 }
0390 if (!elements[iEle].convRefs().empty()) {
0391 active[iEle] = false;
0392 }
0393 LogTrace("PFAlgo|conversionAlgo") << "active[iEle=" << iEle << "]=" << active[iEle];
0394 }
0395 }
0396 LogTrace("PFAlgo|conversionAlgo") << "end of function PFAlgo::conversionAlgo";
0397 }
0398
0399 bool PFAlgo::recoTracksNotHCAL(const reco::PFBlock& block,
0400 reco::PFBlock::LinkData& linkData,
0401 const edm::OwnVector<reco::PFBlockElement>& elements,
0402 const reco::PFBlockRef& blockref,
0403 std::vector<bool>& active,
0404 bool goodTrackDeadHcal,
0405 bool hasDeadHcal,
0406 unsigned int iTrack,
0407 std::multimap<double, unsigned>& ecalElems,
0408 reco::TrackRef& trackRef) {
0409 LogTrace("PFAlgo|recoTracksNotHCAL") << "start of function PFAlgo::recoTracksNotHCAL(), now dealing with tracks "
0410 "linked to no HCAL clusters. Was HCal active? "
0411 << (!hasDeadHcal);
0412
0413
0414
0415
0416 double trackMomentum = elements[iTrack].trackRef()->p();
0417 LogTrace("PFAlgo|recoTracksNotHCAL") << elements[iTrack];
0418
0419
0420
0421
0422 bool thisIsAMuon = PFMuonAlgo::isMuon(elements[iTrack]);
0423 if (thisIsAMuon)
0424 trackMomentum = 0.;
0425
0426
0427
0428 if (!thisIsAMuon && elements[iTrack].trackRef()->ptError() > ptError_) {
0429 double deficit = trackMomentum;
0430 double resol = neutralHadronEnergyResolution(trackMomentum, elements[iTrack].trackRef()->eta());
0431 resol *= trackMomentum;
0432
0433 if (!ecalElems.empty()) {
0434 unsigned thisEcal = ecalElems.begin()->second;
0435 reco::PFClusterRef clusterRef = elements[thisEcal].clusterRef();
0436 bool useCluster = true;
0437 if (hasDeadHcal) {
0438 std::multimap<double, unsigned> sortedTracks;
0439 block.associatedElements(
0440 thisEcal, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
0441 useCluster = (sortedTracks.begin()->second == iTrack);
0442 }
0443 if (useCluster) {
0444 deficit -= clusterRef->energy();
0445 resol = neutralHadronEnergyResolution(trackMomentum, clusterRef->positionREP().Eta());
0446 resol *= trackMomentum;
0447 }
0448 }
0449
0450 bool isPrimary = isFromSecInt(elements[iTrack], "primary");
0451
0452 if (deficit > nSigmaTRACK_ * resol && !isPrimary && !goodTrackDeadHcal) {
0453 active[iTrack] = false;
0454 LogTrace("PFAlgo|recoTracksNotHCAL")
0455 << elements[iTrack] << std::endl
0456 << " deficit " << deficit << " > " << nSigmaTRACK_ << " x " << resol << " track pt " << trackRef->pt()
0457 << " +- " << trackRef->ptError() << " layers valid " << trackRef->hitPattern().trackerLayersWithMeasurement()
0458 << ", lost " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::TRACK_HITS)
0459 << ", lost outer " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_OUTER_HITS)
0460 << ", lost inner " << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_INNER_HITS)
0461 << "(valid fraction " << trackRef->validFraction() << ")"
0462 << " chi2/ndf " << trackRef->normalizedChi2() << " |dxy| "
0463 << std::abs(trackRef->dxy(primaryVertex_.position())) << " +- " << trackRef->dxyError() << " |dz| "
0464 << std::abs(trackRef->dz(primaryVertex_.position())) << " +- " << trackRef->dzError()
0465 << "is probably a fake (1) --> lock the track";
0466 return true;
0467 }
0468 }
0469
0470
0471
0472
0473 std::vector<unsigned> tmpi;
0474 std::vector<unsigned> kTrack;
0475
0476
0477 double dpt = trackRef->ptError();
0478
0479 if (rejectTracks_Step45_ && ecalElems.empty() && trackMomentum > 30. && dpt > 0.5 &&
0480 (PFTrackAlgoTools::step45(trackRef->algo()) && !goodTrackDeadHcal)) {
0481 double dptRel = dpt / trackRef->pt() * 100;
0482 bool isPrimaryOrSecondary = isFromSecInt(elements[iTrack], "all");
0483
0484 if (isPrimaryOrSecondary && dptRel < dptRel_DispVtx_) {
0485 return true;
0486 };
0487 unsigned nHits = elements[iTrack].trackRef()->hitPattern().trackerLayersWithMeasurement();
0488 unsigned int NLostHit = trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::TRACK_HITS);
0489
0490 LogTrace("PFAlgo|recoTracksNotHCAL") << "A track (algo = " << trackRef->algo() << ") with momentum "
0491 << trackMomentum << " / " << elements[iTrack].trackRef()->pt() << " +/- "
0492 << dpt << " / " << elements[iTrack].trackRef()->eta()
0493 << " without any link to ECAL/HCAL and with " << nHits << " (" << NLostHit
0494 << ") hits (lost hits) has been cleaned";
0495
0496 active[iTrack] = false;
0497 return true;
0498 }
0499
0500 tmpi.push_back(reconstructTrack(elements[iTrack]));
0501
0502 kTrack.push_back(iTrack);
0503 active[iTrack] = false;
0504
0505
0506 if (ecalElems.empty()) {
0507 (*pfCandidates_)[tmpi[0]].setEcalEnergy(0., 0.);
0508 (*pfCandidates_)[tmpi[0]].setHcalEnergy(0., 0.);
0509 (*pfCandidates_)[tmpi[0]].setHoEnergy(0., 0.);
0510 (*pfCandidates_)[tmpi[0]].setPs1Energy(0);
0511 (*pfCandidates_)[tmpi[0]].setPs2Energy(0);
0512 (*pfCandidates_)[tmpi[0]].addElementInBlock(blockref, kTrack[0]);
0513 return true;
0514 }
0515
0516
0517 const unsigned int thisEcal = ecalElems.begin()->second;
0518 reco::PFClusterRef clusterRef = elements[thisEcal].clusterRef();
0519 LogTrace("PFAlgo|recoTracksNotHCAL") << " is associated to " << elements[thisEcal];
0520
0521
0522 if (thisIsAMuon) {
0523 (*pfCandidates_)[tmpi[0]].setEcalEnergy(clusterRef->energy(), std::min(clusterRef->energy(), muonECAL_[0]));
0524 (*pfCandidates_)[tmpi[0]].setHcalEnergy(0., 0.);
0525 (*pfCandidates_)[tmpi[0]].setHoEnergy(0., 0.);
0526 (*pfCandidates_)[tmpi[0]].setPs1Energy(0);
0527 (*pfCandidates_)[tmpi[0]].setPs2Energy(0);
0528 (*pfCandidates_)[tmpi[0]].addElementInBlock(blockref, kTrack[0]);
0529 }
0530
0531 double slopeEcal = 1.;
0532 bool connectedToEcal = false;
0533 unsigned iEcal = 99999;
0534 double calibEcal = 0.;
0535 double calibHcal = 0.;
0536 double totalEcal = thisIsAMuon ? -muonECAL_[0] : 0.;
0537
0538
0539 std::multimap<double, unsigned> sortedTracks;
0540 block.associatedElements(thisEcal, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
0541 LogTrace("PFAlgo|recoTracksNotHCAL") << "the closest ECAL cluster, id " << thisEcal << ", has " << sortedTracks.size()
0542 << " associated tracks, now processing them. ";
0543
0544 if (hasDeadHcal && !sortedTracks.empty()) {
0545
0546 LogTrace("PFAlgo|recoTracksNotHCAL") << " the closest track to ECAL " << thisEcal << " is "
0547 << sortedTracks.begin()->second << " (distance " << sortedTracks.begin()->first
0548 << ")";
0549 if (sortedTracks.begin()->second != iTrack) {
0550 LogTrace("PFAlgo|recoTracksNotHCAL")
0551 << " the closest track to ECAL " << thisEcal << " is " << sortedTracks.begin()->second
0552 << " which is not the one being processed. Will skip ECAL linking for this track";
0553 (*pfCandidates_)[tmpi[0]].setEcalEnergy(0., 0.);
0554 (*pfCandidates_)[tmpi[0]].setHcalEnergy(0., 0.);
0555 (*pfCandidates_)[tmpi[0]].setHoEnergy(0., 0.);
0556 (*pfCandidates_)[tmpi[0]].setPs1Energy(0);
0557 (*pfCandidates_)[tmpi[0]].setPs2Energy(0);
0558 (*pfCandidates_)[tmpi[0]].addElementInBlock(blockref, kTrack[0]);
0559 return true;
0560 } else {
0561 LogTrace("PFAlgo|recoTracksNotHCAL")
0562 << " the closest track to ECAL " << thisEcal << " is " << sortedTracks.begin()->second
0563 << " which is the one being processed. Will skip ECAL linking for all other track";
0564 sortedTracks.clear();
0565 }
0566 }
0567
0568 for (auto const& trk : sortedTracks) {
0569 unsigned jTrack = trk.second;
0570
0571
0572 if (!active[jTrack])
0573 continue;
0574
0575
0576 if (jTrack == iTrack)
0577 continue;
0578
0579
0580
0581
0582
0583
0584
0585 std::multimap<double, unsigned> sortedECAL;
0586 block.associatedElements(jTrack, linkData, sortedECAL, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
0587 if (sortedECAL.begin()->second != thisEcal)
0588 continue;
0589
0590
0591 bool thatIsAMuon = PFMuonAlgo::isMuon(elements[jTrack]);
0592 LogTrace("PFAlgo|recoTracksNotHCAL") << " found track " << jTrack << (thatIsAMuon ? " (muon) " : " (non-muon)")
0593 << ", with distance = " << sortedECAL.begin()->first;
0594
0595
0596 bool rejectFake = false;
0597 reco::TrackRef trackRef = elements[jTrack].trackRef();
0598 if (!thatIsAMuon && trackRef->ptError() > ptError_) {
0599
0600
0601 double deficit = trackMomentum + trackRef->p() - clusterRef->energy();
0602 double resol =
0603 nSigmaTRACK_ * neutralHadronEnergyResolution(trackMomentum + trackRef->p(), clusterRef->positionREP().Eta());
0604 resol *= (trackMomentum + trackRef->p());
0605 if (deficit > nSigmaTRACK_ * resol && !goodTrackDeadHcal) {
0606 rejectFake = true;
0607 kTrack.push_back(jTrack);
0608 active[jTrack] = false;
0609 LogTrace("PFAlgo|recoTracksNotHCAL")
0610 << elements[jTrack] << std::endl
0611 << "is probably a fake (2) --> lock the track"
0612 << "(trackMomentum = " << trackMomentum << ", clusterEnergy = " << clusterRef->energy()
0613 << ", deficit = " << deficit << " > " << nSigmaTRACK_ << " x " << resol
0614 << " assuming neutralHadronEnergyResolution "
0615 << neutralHadronEnergyResolution(trackMomentum + trackRef->p(), clusterRef->positionREP().Eta()) << ") ";
0616 }
0617 }
0618 if (rejectFake)
0619 continue;
0620
0621
0622
0623
0624 if (!thatIsAMuon) {
0625 LogTrace("PFAlgo|recoTracksNotHCAL") << "Track momentum increased from " << trackMomentum << " GeV ";
0626 trackMomentum += trackRef->p();
0627 LogTrace("PFAlgo|recoTracksNotHCAL") << "to " << trackMomentum << " GeV.";
0628 LogTrace("PFAlgo|recoTracksNotHCAL") << "with " << elements[jTrack];
0629 } else {
0630 totalEcal -= muonECAL_[0];
0631 totalEcal = std::max(totalEcal, 0.);
0632 }
0633
0634
0635
0636 tmpi.push_back(reconstructTrack(elements[jTrack]));
0637
0638 kTrack.push_back(jTrack);
0639 active[jTrack] = false;
0640
0641 if (thatIsAMuon) {
0642 (*pfCandidates_)[tmpi.back()].setEcalEnergy(clusterRef->energy(), std::min(clusterRef->energy(), muonECAL_[0]));
0643 (*pfCandidates_)[tmpi.back()].setHcalEnergy(0., 0.);
0644 (*pfCandidates_)[tmpi.back()].setHoEnergy(0., 0.);
0645 (*pfCandidates_)[tmpi.back()].setPs1Energy(0);
0646 (*pfCandidates_)[tmpi.back()].setPs2Energy(0);
0647 (*pfCandidates_)[tmpi.back()].addElementInBlock(blockref, kTrack.back());
0648 }
0649 }
0650
0651 LogTrace("PFAlgo|recoTracksNotHCAL") << "Loop over all associated ECAL clusters";
0652
0653 for (auto const& ecal : ecalElems) {
0654 const unsigned index = ecal.second;
0655 const PFBlockElement::Type type = elements[index].type();
0656 assert(type == PFBlockElement::ECAL);
0657 LogTrace("PFAlgo|recoTracksNotHCAL") << elements[index];
0658
0659
0660 if (!active[index]) {
0661 LogTrace("PFAlgo|recoTracksNotHCAL") << "is not active - ignore ";
0662 continue;
0663 }
0664
0665
0666 block.associatedElements(index, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
0667
0668 const bool skip = std::none_of(
0669 kTrack.begin(), kTrack.end(), [&](unsigned iTrack) { return sortedTracks.begin()->second == iTrack; });
0670
0671 if (skip) {
0672 LogTrace("PFAlgo|recoTracksNotHCAL") << "is closer to another track - ignore ";
0673 continue;
0674 }
0675
0676
0677 const reco::PFClusterRef clusterRef = elements[index].clusterRef();
0678 assert(!clusterRef.isNull());
0679
0680 LogTrace("PFAlgo|recoTracksNotHCAL") << "Ecal cluster with raw energy = " << clusterRef->energy()
0681 << " linked with distance = " << ecal.first;
0682
0683
0684
0685
0686
0687 vector<double> ps1Ene{0.};
0688 associatePSClusters(index, reco::PFBlockElement::PS1, block, elements, linkData, active, ps1Ene);
0689 vector<double> ps2Ene{0.};
0690 associatePSClusters(index, reco::PFBlockElement::PS2, block, elements, linkData, active, ps2Ene);
0691
0692
0693 const double ecalEnergy = clusterRef->energy();
0694 const double ecalEnergyCalibrated = clusterRef->correctedEnergy();
0695 LogTrace("PFAlgo|recoTracksNotHCAL") << "Corrected ECAL(+PS) energy = " << ecalEnergy;
0696
0697
0698
0699 totalEcal += ecalEnergy;
0700 const double previousCalibEcal = calibEcal;
0701 const double previousSlopeEcal = slopeEcal;
0702 calibEcal = std::max(totalEcal, 0.);
0703 calibHcal = 0.;
0704 calibration_.energyEmHad(
0705 trackMomentum, calibEcal, calibHcal, clusterRef->positionREP().Eta(), clusterRef->positionREP().Phi());
0706 if (totalEcal > 0.)
0707 slopeEcal = calibEcal / totalEcal;
0708
0709 LogTrace("PFAlgo|recoTracksNotHCAL") << "The total calibrated energy so far amounts to = " << calibEcal
0710 << " (slope = " << slopeEcal << ")";
0711
0712
0713 if (connectedToEcal && calibEcal - trackMomentum >= 0.) {
0714
0715
0716 calibEcal = previousCalibEcal;
0717 slopeEcal = previousSlopeEcal;
0718 totalEcal = calibEcal / slopeEcal;
0719
0720
0721
0722 active[index] = false;
0723
0724
0725 std::multimap<double, unsigned> assTracks;
0726 block.associatedElements(index, linkData, assTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
0727
0728 auto& ecalCand = (*pfCandidates_)[reconstructCluster(
0729 *clusterRef, ecalEnergyCalibrated)];
0730 ecalCand.setEcalEnergy(clusterRef->energy(), ecalEnergyCalibrated);
0731 ecalCand.setHcalEnergy(0., 0.);
0732 ecalCand.setHoEnergy(0., 0.);
0733 ecalCand.setPs1Energy(ps1Ene[0]);
0734 ecalCand.setPs2Energy(ps2Ene[0]);
0735 ecalCand.addElementInBlock(blockref, index);
0736
0737 if (!assTracks.empty()) {
0738 ecalCand.addElementInBlock(blockref, assTracks.begin()->second);
0739
0740
0741 const ::math::XYZPointF& chargedPosition =
0742 dynamic_cast<const reco::PFBlockElementTrack*>(&elements[assTracks.begin()->second])
0743 ->positionAtECALEntrance();
0744 ecalCand.setPositionAtECALEntrance(chargedPosition);
0745 }
0746 break;
0747 }
0748
0749
0750 connectedToEcal = true;
0751 iEcal = index;
0752 active[index] = false;
0753 for (unsigned ic : tmpi)
0754 (*pfCandidates_)[ic].addElementInBlock(blockref, iEcal);
0755
0756 }
0757
0758 bool bNeutralProduced = false;
0759
0760
0761 if (connectedToEcal) {
0762 reco::PFClusterRef pivotalRef = elements[iEcal].clusterRef();
0763
0764 double neutralEnergy = calibEcal - trackMomentum;
0765
0766
0767
0768
0769
0770
0771
0772
0773
0774
0775
0776
0777
0778
0779
0780
0781
0782
0783
0784
0785
0786
0787
0788
0789
0790
0791
0792
0793
0794
0795
0796
0797
0798
0799
0800 double resol = neutralHadronEnergyResolution(trackMomentum, pivotalRef->positionREP().Eta());
0801 resol *= trackMomentum;
0802 if (neutralEnergy > std::max(0.5, nSigmaECAL_ * resol)) {
0803 neutralEnergy /= slopeEcal;
0804 unsigned tmpj = reconstructCluster(*pivotalRef, neutralEnergy);
0805 (*pfCandidates_)[tmpj].setEcalEnergy(pivotalRef->energy(), neutralEnergy);
0806 (*pfCandidates_)[tmpj].setHcalEnergy(0., 0.);
0807 (*pfCandidates_)[tmpj].setHoEnergy(0., 0.);
0808 (*pfCandidates_)[tmpj].setPs1Energy(0.);
0809 (*pfCandidates_)[tmpj].setPs2Energy(0.);
0810 (*pfCandidates_)[tmpj].addElementInBlock(blockref, iEcal);
0811 bNeutralProduced = true;
0812 for (unsigned ic = 0; ic < kTrack.size(); ++ic)
0813 (*pfCandidates_)[tmpj].addElementInBlock(blockref, kTrack[ic]);
0814 }
0815
0816
0817 for (unsigned ic = 0; ic < tmpi.size(); ++ic) {
0818
0819 if ((*pfCandidates_)[tmpi[ic]].particleId() == reco::PFCandidate::mu)
0820 continue;
0821
0822 double fraction = trackMomentum > 0 ? (*pfCandidates_)[tmpi[ic]].trackRef()->p() / trackMomentum : 0;
0823 double ecalCal = bNeutralProduced ? (calibEcal - neutralEnergy * slopeEcal) * fraction : calibEcal * fraction;
0824 double ecalRaw = totalEcal * fraction;
0825
0826 LogTrace("PFAlgo|recoTracksNotHCAL")
0827 << "The fraction after photon supression is " << fraction << " calibrated ecal = " << ecalCal;
0828
0829 (*pfCandidates_)[tmpi[ic]].setEcalEnergy(ecalRaw, ecalCal);
0830 (*pfCandidates_)[tmpi[ic]].setHcalEnergy(0., 0.);
0831 (*pfCandidates_)[tmpi[ic]].setHoEnergy(0., 0.);
0832 (*pfCandidates_)[tmpi[ic]].setPs1Energy(0);
0833 (*pfCandidates_)[tmpi[ic]].setPs2Energy(0);
0834 (*pfCandidates_)[tmpi[ic]].addElementInBlock(blockref, kTrack[ic]);
0835 }
0836
0837 }
0838
0839
0840 for (unsigned ic = 0; ic < tmpi.size(); ++ic) {
0841 const PFCandidate& pfc = (*pfCandidates_)[tmpi[ic]];
0842 const PFCandidate::ElementsInBlocks& eleInBlocks = pfc.elementsInBlocks();
0843 if (eleInBlocks.empty()) {
0844 LogTrace("PFAlgo|recoTracksNotHCAL") << "Single track / Fill element in block! ";
0845 (*pfCandidates_)[tmpi[ic]].addElementInBlock(blockref, kTrack[ic]);
0846 }
0847 }
0848 LogTrace("PFAlgo|recoTracksNotHCAL") << "end of function PFAlgo::recoTracksNotHCAL";
0849 return false;
0850 }
0851
0852
0853
0854
0855 bool PFAlgo::checkAndReconstructSecondaryInteraction(const reco::PFBlockRef& blockref,
0856 const edm::OwnVector<reco::PFBlockElement>& elements,
0857 bool isActive,
0858 int iElement) {
0859 bool ret = isActive;
0860 if (isActive && isFromSecInt(elements[iElement], "primary")) {
0861 bool isPrimaryTrack =
0862 elements[iElement].displacedVertexRef(PFBlockElement::T_TO_DISP)->displacedVertexRef()->isTherePrimaryTracks();
0863 if (isPrimaryTrack) {
0864 LogTrace("PFAlgo|elementLoop") << "Primary Track reconstructed alone";
0865
0866 unsigned tmpi = reconstructTrack(elements[iElement]);
0867 (*pfCandidates_)[tmpi].addElementInBlock(blockref, iElement);
0868 ret = false;
0869 }
0870 }
0871
0872 return ret;
0873 }
0874
0875 bool PFAlgo::checkHasDeadHcal(const std::multimap<double, unsigned>& hcalElems, const std::vector<bool>& deadArea) {
0876
0877
0878
0879
0880
0881
0882
0883
0884
0885
0886
0887
0888 bool hasDeadHcal = false;
0889 if (!hcalElems.empty() && deadArea[hcalElems.begin()->second]) {
0890 hasDeadHcal = true;
0891 }
0892
0893
0894
0895
0896
0897
0898 return hasDeadHcal;
0899 }
0900
0901
0902 bool PFAlgo::checkGoodTrackDeadHcal(const reco::TrackRef& trackRef, bool hasDeadHcal) {
0903 bool goodTrackDeadHcal = false;
0904 if (hasDeadHcal) {
0905 goodTrackDeadHcal = (trackRef->ptError() < goodTrackDeadHcal_ptErrRel_ * trackRef->pt() &&
0906 trackRef->normalizedChi2() < goodTrackDeadHcal_chi2n_ &&
0907 trackRef->hitPattern().trackerLayersWithMeasurement() >= goodTrackDeadHcal_layers_ &&
0908 trackRef->validFraction() > goodTrackDeadHcal_validFr_ &&
0909 std::abs(trackRef->dxy(primaryVertex_.position())) < goodTrackDeadHcal_dxy_);
0910
0911 if (!goodTrackDeadHcal && std::abs(trackRef->eta()) > goodPixelTrackDeadHcal_minEta_ &&
0912 trackRef->hitPattern().pixelLayersWithMeasurement() >= 3 &&
0913 trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::TRACK_HITS) == 0 &&
0914 trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_INNER_HITS) == 0 &&
0915 trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_OUTER_HITS) <=
0916 (trackRef->hitPattern().pixelLayersWithMeasurement() > 3 ? goodPixelTrackDeadHcal_maxLost4Hit_
0917 : goodPixelTrackDeadHcal_maxLost3Hit_) &&
0918 trackRef->normalizedChi2() < goodPixelTrackDeadHcal_chi2n_ &&
0919 std::abs(trackRef->dxy(primaryVertex_.position())) < goodPixelTrackDeadHcal_dxy_ &&
0920 std::abs(trackRef->dz(primaryVertex_.position())) < goodPixelTrackDeadHcal_dz_ &&
0921 trackRef->ptError() < goodPixelTrackDeadHcal_ptErrRel_ * trackRef->pt() &&
0922 trackRef->pt() < goodPixelTrackDeadHcal_maxPt_) {
0923 goodTrackDeadHcal = true;
0924
0925 }
0926
0927 LogTrace("PFAlgo|elementLoop")
0928 << " track pt " << trackRef->pt() << " +- " << trackRef->ptError() << " layers valid "
0929 << trackRef->hitPattern().trackerLayersWithMeasurement() << ", lost "
0930 << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::TRACK_HITS) << ", lost outer "
0931 << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_OUTER_HITS) << ", lost inner "
0932 << trackRef->hitPattern().trackerLayersWithoutMeasurement(HitPattern::MISSING_INNER_HITS) << "(valid fraction "
0933 << trackRef->validFraction() << ")"
0934 << " chi2/ndf " << trackRef->normalizedChi2() << " |dxy| " << std::abs(trackRef->dxy(primaryVertex_.position()))
0935 << " +- " << trackRef->dxyError() << " |dz| " << std::abs(trackRef->dz(primaryVertex_.position())) << " +- "
0936 << trackRef->dzError() << (goodTrackDeadHcal ? " passes " : " fails ") << "quality cuts";
0937 }
0938 return goodTrackDeadHcal;
0939 }
0940
0941 void PFAlgo::relinkTrackToHcal(const reco::PFBlock& block,
0942 std::multimap<double, unsigned>& ecalElems,
0943 std::multimap<double, unsigned>& hcalElems,
0944 const std::vector<bool>& active,
0945 reco::PFBlock::LinkData& linkData,
0946 unsigned int iTrack) {
0947 unsigned ntt = 1;
0948 unsigned index = ecalElems.begin()->second;
0949 std::multimap<double, unsigned> sortedTracks;
0950 block.associatedElements(index, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
0951 LogTrace("PFAlgo|elementLoop") << "The closest ECAL cluster is linked to " << sortedTracks.size()
0952 << " tracks, with distance = " << ecalElems.begin()->first;
0953
0954 LogTrace("PFAlgo|elementLoop") << "Looping over sortedTracks";
0955
0956 for (auto const& trk : sortedTracks) {
0957 unsigned jTrack = trk.second;
0958 LogTrace("PFAlgo|elementLoop") << "jTrack=" << jTrack;
0959
0960 if (!active[jTrack])
0961 continue;
0962 LogTrace("PFAlgo|elementLoop") << "active[jTrack]=" << active[jTrack];
0963
0964
0965 if (jTrack == iTrack)
0966 continue;
0967 LogTrace("PFAlgo|elementLoop") << "skipping jTrack=" << jTrack << " for same iTrack";
0968
0969
0970
0971 std::multimap<double, unsigned> sortedECAL;
0972 block.associatedElements(jTrack, linkData, sortedECAL, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
0973 if (sortedECAL.begin()->second != index)
0974 continue;
0975 LogTrace("PFAlgo|elementLoop") << " track " << jTrack << " with closest ECAL identical ";
0976
0977
0978 std::multimap<double, unsigned> sortedHCAL;
0979 block.associatedElements(jTrack, linkData, sortedHCAL, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
0980 if (sortedHCAL.empty())
0981 continue;
0982 LogTrace("PFAlgo|elementLoop") << " and with an HCAL cluster " << sortedHCAL.begin()->second;
0983 ntt++;
0984
0985
0986 block.setLink(iTrack, sortedHCAL.begin()->second, sortedECAL.begin()->first, linkData, PFBlock::LINKTEST_RECHIT);
0987
0988 }
0989
0990
0991 block.associatedElements(iTrack, linkData, hcalElems, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
0992
0993 if (!hcalElems.empty())
0994 LogTrace("PFAlgo|elementLoop") << "Track linked back to HCAL due to ECAL sharing with other tracks";
0995 }
0996
0997 void PFAlgo::elementLoop(const reco::PFBlock& block,
0998 reco::PFBlock::LinkData& linkData,
0999 const edm::OwnVector<reco::PFBlockElement>& elements,
1000 std::vector<bool>& active,
1001 const reco::PFBlockRef& blockref,
1002 ElementIndices& inds,
1003 std::vector<bool>& deadArea) {
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
1011
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
1036
1037
1038
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
1071
1072
1073 if (hcalElems.empty() && !ecalElems.empty() && !hasDeadHcal) {
1074 relinkTrackToHcal(block, ecalElems, hcalElems, active, linkData, iEle);
1075 }
1076
1077
1078
1079
1080
1081
1082
1083
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
1092 bool hcalFound = false;
1093 bool hfhadFound = false;
1094
1095 LogTrace("PFAlgo|elementLoop") << "now looping on elements associated to the track: ecalElems";
1096
1097
1098
1099 for (auto const& ecal : ecalElems) {
1100 unsigned index = ecal.second;
1101
1102 PFBlockElement::Type type = elements[index].type();
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
1109 assert(type == PFBlockElement::ECAL);
1110
1111
1112 if (!active[index])
1113 continue;
1114
1115
1116
1117
1118
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 }
1124
1125
1126
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 }
1135
1136
1137
1138 for (auto const& hcal : hcalElems) {
1139 unsigned index = hcal.second;
1140
1141 PFBlockElement::Type type = elements[index].type();
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
1147 assert(type == PFBlockElement::HCAL);
1148
1149
1150
1151 if (!hcalFound) {
1152 LogTrace("PFAlgo|elementLoop") << "\t\tclosest hcal cluster, doing nothing";
1153
1154 hcalFound = true;
1155
1156
1157
1158 } else {
1159
1160 LogTrace("PFAlgo|elementLoop") << "\t\tsecondary hcal cluster. unlinking";
1161 block.setLink(iEle, index, -1., linkData, PFBlock::LINKTEST_RECHIT);
1162 }
1163 }
1164
1165
1166
1167
1168 for (auto const& hfhad : hfHadElems) {
1169 unsigned index = hfhad.second;
1170
1171 PFBlockElement::Type type = elements[index].type();
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
1180
1181 if (!hfhadFound) {
1182 LogTrace("PFAlgo|elementLoop") << "\t\tclosest hfhad cluster, doing nothing";
1183
1184 hfhadFound = true;
1185
1186 } else {
1187
1188 LogTrace("PFAlgo|elementLoop") << "\t\tsecondary hfhad cluster. unlinking";
1189 block.setLink(iEle, index, -1., linkData, PFBlock::LINKTEST_RECHIT);
1190 }
1191 }
1192
1193 LogTrace("PFAlgo|elementLoop") << "end of loop over iEle";
1194 }
1195 LogTrace("PFAlgo|elementLoop") << "end of function PFAlgo::elementLoop";
1196 }
1197
1198
1199
1200
1201 int PFAlgo::decideType(const edm::OwnVector<reco::PFBlockElement>& elements,
1202 const reco::PFBlockElement::Type type,
1203 std::vector<bool>& active,
1204 ElementIndices& inds,
1205 std::vector<bool>& deadArea,
1206 unsigned int iEle) {
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;
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;
1227 }
1228 inds.hcalIs.push_back(iEle);
1229 LogTrace("PFAlgo|decideType") << "HCAL, stored index, continue";
1230 }
1231 return 1;
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;
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;
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;
1252 default:
1253 return 1;
1254 }
1255 LogTrace("PFAlgo|decideType") << "Did not match type to anything, return 0";
1256 return 0;
1257 }
1258
1259 void PFAlgo::createCandidatesHF(const reco::PFBlock& block,
1260 reco::PFBlock::LinkData& linkData,
1261 const edm::OwnVector<reco::PFBlockElement>& elements,
1262 std::vector<bool>& active,
1263 const reco::PFBlockRef& blockref,
1264 ElementIndices& inds) {
1265 LogTrace("PFAlgo|createCandidatesHF") << "starting function PFAlgo::createCandidatesHF";
1266
1267 bool trackInBlock = !inds.trackIs.empty();
1268
1269
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
1279
1280 if (!trackInBlock)
1281 assert(inds.hfEmIs.size() + inds.hfHadIs.size() == elements.size());
1282
1283
1284
1285
1286
1287 if (trackInBlock) {
1288
1289 std::multimap<double, unsigned> sortedTracks;
1290 std::multimap<double, unsigned> sortedTracksActive;
1291
1292 std::multimap<unsigned, std::pair<double, unsigned>> associatedHfEms;
1293
1294 std::multimap<double, std::pair<unsigned, double>> hfemSatellites;
1295
1296
1297
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
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
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
1332
1333 for (auto const& trk : sortedTracks) {
1334 unsigned iTrack = trk.second;
1335
1336
1337 if (!active[iTrack])
1338 continue;
1339
1340 PFBlockElement::Type type = elements[iTrack].type();
1341 assert(type == reco::PFBlockElement::TRACK);
1342
1343 reco::TrackRef trackRef = elements[iTrack].trackRef();
1344 assert(!trackRef.isNull());
1345
1346
1347 double trackMomentum = trackRef->p();
1348 totalChargedMomentum += trackMomentum;
1349
1350
1351 double dp = trackRef->qoverpError() * trackMomentum * trackMomentum;
1352 sumpError2 += dp * dp;
1353
1354
1355 sortedTracksActive.emplace(trk);
1356
1357
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;
1365
1366
1367
1368
1369 for (auto const& hfem : sortedHfEms) {
1370 unsigned iHfEm = hfem.second;
1371 double dist = hfem.first;
1372
1373
1374 if (!active[iHfEm]) {
1375 LogTrace("PFAlgo|createCandidatesHF") << "cluster locked";
1376 continue;
1377 }
1378
1379
1380 PFBlockElement::Type type = elements[iHfEm].type();
1381 assert(type == PFBlockElement::HFEM);
1382 PFClusterRef eclusterRef = elements[iHfEm].clusterRef();
1383 assert(!eclusterRef.isNull());
1384
1385
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) {
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 {
1404
1405
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 }
1414 }
1415
1416
1417 double uncalibratedenergyHfHad = hclusterRef->energy();
1418 double energyHfHad = uncalibratedenergyHfHad;
1419 if (thepfEnergyCalibrationHF_.getcalibHF_use()) {
1420 energyHfHad = thepfEnergyCalibrationHF_.energyHad(
1421 uncalibratedenergyHfHad,
1422 hclusterRef->positionREP().Eta(),
1423 hclusterRef->positionREP().Phi());
1424 }
1425 double calibFactorHfHad = (uncalibratedenergyHfHad > 0.) ? energyHfHad / uncalibratedenergyHfHad : 1.;
1426
1427
1428 double energyHfEmTmp = 0.;
1429 double uncalibratedenergyHfEmTmp = 0.;
1430 double energyHfEm = 0.;
1431 double uncalibratedenergyHfEm = 0.;
1432
1433
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
1441 if (sortedTracksActive.empty()) {
1442
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
1449
1450 if (!sortedHfEms.empty()) {
1451 for (auto const& hfem : sortedHfEms) {
1452 unsigned iHfEm = hfem.second;
1453
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;
1462 if (thepfEnergyCalibrationHF_.getcalibHF_use()) {
1463 energyHfEm = thepfEnergyCalibrationHF_.energyEmHad(
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 }
1468 }
1469 }
1470
1471
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 }
1483
1484
1485
1486 else if ((energyHfHad - totalChargedMomentum) > nsigmaHFHAD * totalError) {
1487 assert(energyHfEm == 0.);
1488
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
1500
1501
1502
1503 else {
1504 for (auto const& hfemSatellite : hfemSatellites) {
1505
1506 uncalibratedenergyHfEmTmp += std::get<1>(hfemSatellite.second);
1507 energyHfEmTmp = uncalibratedenergyHfEmTmp;
1508 double energyHfHadTmp = uncalibratedenergyHfHad;
1509 unsigned iHfEm = std::get<0>(hfemSatellite.second);
1510 PFClusterRef eclusterRef = elements[iHfEm].clusterRef();
1511 assert(!eclusterRef.isNull());
1512 if (thepfEnergyCalibrationHF_.getcalibHF_use()) {
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
1523
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
1529 if (hfemSatellite.first < 0. && (caloEnergyTmp - totalChargedMomentum) > nsigmaHFEM * totalError) {
1530
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 }
1547 }
1548
1549
1550
1551
1552 for (auto const& trk : sortedTracksActive) {
1553 unsigned iTrack = trk.second;
1554
1555
1556 reco::TrackRef trackRef = elements[iTrack].trackRef();
1557 assert(!trackRef.isNull());
1558
1559
1560
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 }
1579
1580 }
1581
1582
1583
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
1593 energyHF = eclusterRef->energy();
1594 uncalibratedenergyHF = energyHF;
1595 if (thepfEnergyCalibrationHF_.getcalibHF_use()) {
1596 energyHF = thepfEnergyCalibrationHF_.energyEm(
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 }
1607
1608 }
1609
1610
1611
1612
1613
1614 else if (elements.size() == 1) {
1615
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
1623 energyHF = clusterRef->energy();
1624 uncalibratedenergyHF = energyHF;
1625 if (thepfEnergyCalibrationHF_.getcalibHF_use()) {
1626 energyHF = thepfEnergyCalibrationHF_.energyEm(
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
1640 energyHF = clusterRef->energy();
1641 uncalibratedenergyHF = energyHF;
1642 if (thepfEnergyCalibrationHF_.getcalibHF_use()) {
1643 energyHF = thepfEnergyCalibrationHF_.energyHad(
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
1660 reco::PFClusterRef c0 = elements[0].clusterRef();
1661 reco::PFClusterRef c1 = elements[1].clusterRef();
1662
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
1671
1672 }
1673
1674 double energyHfEm = cem->energy();
1675 double energyHfHad = chad->energy();
1676 double uncalibratedenergyHfEm = energyHfEm;
1677 double uncalibratedenergyHfHad = energyHfHad;
1678 if (thepfEnergyCalibrationHF_.getcalibHF_use()) {
1679 energyHfEm = thepfEnergyCalibrationHF_.energyEmHad(
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
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 }
1701
1702 void PFAlgo::createCandidatesHCAL(const reco::PFBlock& block,
1703 reco::PFBlock::LinkData& linkData,
1704 const edm::OwnVector<reco::PFBlockElement>& elements,
1705 std::vector<bool>& active,
1706 const reco::PFBlockRef& blockref,
1707 ElementIndices& inds,
1708 std::vector<bool>& deadArea) {
1709 LogTrace("PFAlgo|createCandidatesHCAL")
1710 << "start of function PFAlgo::createCandidatesHCAL, inds.hcalIs.size()=" << inds.hcalIs.size();
1711
1712
1713
1714 for (unsigned iHcal : inds.hcalIs) {
1715 PFBlockElement::Type type = elements[iHcal].type();
1716
1717 assert(type == PFBlockElement::HCAL);
1718
1719 LogTrace("PFAlgo|createCandidatesHCAL") << "elements[" << iHcal << "]=" << elements[iHcal];
1720
1721
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
1732 std::multimap<double, std::tuple<unsigned, ::math::XYZVector, double>>
1733 ecalSatellites;
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
1743
1744
1745
1746
1747
1748
1749 if (sortedTracks.empty()) {
1750 LogTrace("PFAlgo|createCandidatesHCAL") << "\tno associated tracks, keep for later";
1751 continue;
1752 }
1753
1754
1755 active[iHcal] = false;
1756
1757
1758
1759
1760
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
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;
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
1796 for (auto const& trk : sortedTracks) {
1797 unsigned iTrack = trk.second;
1798
1799
1800 if (!active[iTrack])
1801 continue;
1802
1803 PFBlockElement::Type type = elements[iTrack].type();
1804 assert(type == reco::PFBlockElement::TRACK);
1805
1806 reco::TrackRef trackRef = elements[iTrack].trackRef();
1807 assert(!trackRef.isNull());
1808
1809
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
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
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
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
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
1853
1854
1855 bool letMuonEatCaloEnergy = false;
1856
1857 if (thisIsAnIsolatedMuon) {
1858
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
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
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
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
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
1943 photonAtECAL -= muonECAL_[0] * chargedDirection;
1944 hadronAtECAL -= muonHCAL_[0] * chargedDirection;
1945 }
1946
1947
1948 active[iTrack] = false;
1949
1950 continue;
1951 }
1952
1953
1954
1955 LogTrace("PFAlgo|createCandidatesHCAL") << "elements[iTrack=" << iTrack << "]=" << elements[iTrack];
1956
1957
1958 double trackMomentum = trackRef->p();
1959 totalChargedMomentum += trackMomentum;
1960
1961
1962
1963 if (thisIsALooseMuon && !thisIsAMuon)
1964 nMuons += 1;
1965
1966
1967
1968
1969 double dpt = trackRef->ptError();
1970 double blowError = PFTrackAlgoTools::errorScale(trackRef->algo(), factors45_);
1971
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
1981 double dp = trackRef->qoverpError() * trackMomentum * trackMomentum;
1982 sumpError2 += dp * dp;
1983
1984 bool connectedToEcal = false;
1985 if (!sortedEcals.empty()) {
1986
1987
1988 for (auto const& ecal : sortedEcals) {
1989 unsigned iEcal = ecal.second;
1990 double dist = ecal.first;
1991
1992
1993 if (!active[iEcal]) {
1994 LogTrace("PFAlgo|createCandidatesHCAL") << "cluster locked";
1995 continue;
1996 }
1997
1998
1999 PFBlockElement::Type type = elements[iEcal].type();
2000 assert(type == PFBlockElement::ECAL);
2001 PFClusterRef eclusterref = elements[iEcal].clusterRef();
2002 assert(!eclusterref.isNull());
2003
2004
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();
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) {
2021
2022 LogTrace("PFAlgo|createCandidatesHCAL") << "closest: " << elements[iEcal];
2023
2024 connectedToEcal = true;
2025
2026
2027
2028
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 {
2035
2036
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 }
2047 }
2048
2049 if (useHO_ && !sortedHOs.empty()) {
2050
2051
2052 for (auto const& ho : sortedHOs) {
2053 unsigned iHO = ho.second;
2054 double distHO = ho.first;
2055
2056
2057 if (!active[iHO]) {
2058 LogTrace("PFAlgo|createCandidatesHCAL") << "cluster locked";
2059 continue;
2060 }
2061
2062
2063 PFBlockElement::Type type = elements[iHO].type();
2064 assert(type == PFBlockElement::HO);
2065 PFClusterRef hoclusterref = elements[iHO].clusterRef();
2066 assert(!hoclusterref.isNull());
2067
2068
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
2077
2078
2079
2080
2081
2082 totalHO += hoclusterref->energy();
2083 active[iHO] = false;
2084
2085 std::pair<double, unsigned> associatedHO(distHO, iHO);
2086 associatedHOs.emplace(iTrack, associatedHO);
2087
2088 }
2089 }
2090
2091 }
2092
2093
2094 totalHcal += totalHO;
2095
2096
2097
2098 double caloEnergy = 0.;
2099 double slopeEcal = 1.0;
2100 double calibEcal = 0.;
2101 double calibHcal = 0.;
2102 hadronDirection = hadronAtECAL.Unit();
2103
2104
2105 double caloResolution = neutralHadronEnergyResolution(totalChargedMomentum, hclusterref->positionREP().Eta());
2106 caloResolution *= totalChargedMomentum;
2107
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
2118
2119
2120
2121 for (auto const& ecalSatellite : ecalSatellites) {
2122
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());
2131 totalEcalEGMCalib += sqrt(std::get<1>(ecalSatellite.second).Mag2()) *
2132 std::get<2>(ecalSatellite.second);
2133 photonAtECAL += std::get<1>(ecalSatellite.second) *
2134 std::get<2>(ecalSatellite.second);
2135 calibEcal = std::max(0., totalEcal);
2136 calibHcal = std::max(0., totalHcal);
2137 hadronAtECAL = calibHcal * hadronDirection;
2138
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
2151
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);
2159 if (clusterEnergy > 50) {
2160 ecalClusters.push_back(ecalSatellite.second);
2161 sumEcalClusters += clusterEnergy;
2162 }
2163 continue;
2164 }
2165
2166
2167
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
2181 assert(caloEnergy >= 0);
2182
2183
2184
2185
2186
2187
2188
2189
2190
2191
2192
2193
2194
2195 if (totalChargedMomentum - caloEnergy > nSigmaTRACK_ * caloResolution) {
2196
2197 if (nMuons > 0) {
2198 for (auto const& trk : associatedTracks) {
2199
2200 if (!trk.second.second)
2201 continue;
2202
2203 const unsigned int iTrack = trk.second.first;
2204
2205 if (!active[iTrack])
2206 continue;
2207
2208 const double trackMomentum = elements[trk.second.first].trackRef()->p();
2209
2210
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
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
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
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
2273
2274
2275 }
2276
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
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
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
2314 if (fabs(trk.first) < ptError_)
2315 break;
2316
2317 const double wouldBeTotalChargedMomentum = totalChargedMomentum - trackref->p();
2318
2319
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
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
2356 caloResolution = neutralHadronEnergyResolution(totalChargedMomentum, hclusterref->positionREP().Eta());
2357 caloResolution *= totalChargedMomentum;
2358 caloResolution = std::sqrt(caloResolution * caloResolution + muonHCALError + muonECALError);
2359
2360
2361
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
2389 caloResolution = neutralHadronEnergyResolution(totalChargedMomentum, hclusterref->positionREP().Eta());
2390 caloResolution *= totalChargedMomentum;
2391 caloResolution = std::sqrt(caloResolution * caloResolution + muonHCALError + muonECALError);
2392
2393
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.;
2428 (*pfCandidates_)[tmpi].rescaleMomentum(corrFact);
2429 trackMomentum *= corrFact;
2430 }
2431 chargedHadronsIndices.push_back(tmpi);
2432 chargedHadronsInBlock.push_back(iTrack);
2433 active[iTrack] = false;
2434 hcalP.push_back(trackMomentum);
2435 hcalDP.push_back(dp);
2436 if (dp / trackMomentum > maxDPovP)
2437 maxDPovP = dp / trackMomentum;
2438 sumpError2 += dp * dp;
2439 }
2440
2441
2442 double totalError = sqrt(sumpError2 + caloResolution * caloResolution);
2443
2444 #ifdef EDM_ML_DEBUG
2445 LogTrace("PFAlgo|createCandidatesHCAL")
2446 << "\tCompare Calo Energy to total charged momentum " << endl
2447 << "\t\tsum p = " << totalChargedMomentum << " +- " << sqrt(sumpError2) << endl
2448 << "\t\tsum ecal = " << totalEcal << endl
2449 << "\t\tsum hcal = " << totalHcal << endl
2450 << "\t\t => Calo Energy = " << caloEnergy << " +- " << caloResolution << endl
2451 << "\t\t => Calo Energy- total charged momentum = " << caloEnergy - totalChargedMomentum << " +- "
2452 << totalError;
2453 #endif
2454
2455
2456
2457
2458 double nsigma = nSigmaHCAL(totalChargedMomentum, hclusterref->positionREP().Eta());
2459
2460 if (abs(totalChargedMomentum - caloEnergy) < nsigma * totalError) {
2461
2462
2463
2464 #ifdef EDM_ML_DEBUG
2465 LogTrace("PFAlgo|createCandidatesHCAL")
2466 << "\t\tcase 1: COMPATIBLE "
2467 << "|Calo Energy- total charged momentum| = " << abs(caloEnergy - totalChargedMomentum) << " < " << nsigma
2468 << " x " << totalError;
2469 if (maxDPovP < 0.1)
2470 LogTrace("PFAlgo|createCandidatesHCAL") << "\t\t\tmax DP/P = " << maxDPovP << " less than 0.1: do nothing ";
2471 else
2472 LogTrace("PFAlgo|createCandidatesHCAL")
2473 << "\t\t\tmax DP/P = " << maxDPovP << " > 0.1: take weighted averages ";
2474 #endif
2475
2476
2477 if (maxDPovP > 0.1) {
2478
2479
2480 int nrows = chargedHadronsIndices.size();
2481 TMatrixTSym<double> a(nrows);
2482 TVectorD b(nrows);
2483 TVectorD check(nrows);
2484 double sigma2E = caloResolution * caloResolution;
2485 for (int i = 0; i < nrows; i++) {
2486 double sigma2i = hcalDP[i] * hcalDP[i];
2487 LogTrace("PFAlgo|createCandidatesHCAL")
2488 << "\t\t\ttrack associated to hcal " << i << " P = " << hcalP[i] << " +- " << hcalDP[i];
2489 a(i, i) = 1. / sigma2i + 1. / sigma2E;
2490 b(i) = hcalP[i] / sigma2i + caloEnergy / sigma2E;
2491 for (int j = 0; j < nrows; j++) {
2492 if (i == j)
2493 continue;
2494 a(i, j) = 1. / sigma2E;
2495 }
2496 }
2497
2498
2499 TDecompChol decomp(a);
2500 bool ok = false;
2501 TVectorD x = decomp.Solve(b, ok);
2502
2503
2504 if (ok) {
2505 for (int i = 0; i < nrows; i++) {
2506
2507 unsigned ich = chargedHadronsIndices[i];
2508 double rescaleFactor = x(i) / hcalP[i];
2509 if (rescaleFactor < 0.)
2510 rescaleFactor = 0.;
2511 (*pfCandidates_)[ich].rescaleMomentum(rescaleFactor);
2512
2513 LogTrace("PFAlgo|createCandidatesHCAL")
2514 << "\t\t\told p " << hcalP[i] << " new p " << x(i) << " rescale " << rescaleFactor;
2515 }
2516 } else {
2517 edm::LogError("PFAlgo::createCandidatesHCAL") << "TDecompChol.Solve returned ok=false";
2518 assert(0);
2519 }
2520 }
2521 }
2522
2523
2524 else if (caloEnergy > totalChargedMomentum) {
2525
2526
2527
2528
2529 double eNeutralHadron = caloEnergy - totalChargedMomentum;
2530 double ePhoton = (caloEnergy - totalChargedMomentum) /
2531 slopeEcal;
2532
2533
2534
2535 #ifdef EDM_ML_DEBUG
2536 if (!sortedTracks.empty()) {
2537 LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tcase 2: NEUTRAL DETECTION " << caloEnergy << " > " << nsigma
2538 << "x" << totalError << " + " << totalChargedMomentum;
2539 LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tneutral activity detected: " << endl
2540 << "\t\t\t photon = " << ePhoton << endl
2541 << "\t\t\tor neutral hadron = " << eNeutralHadron;
2542
2543 LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tphoton or hadron ?";
2544 }
2545
2546 if (sortedTracks.empty())
2547 LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tno track -> hadron ";
2548 else
2549 LogTrace("PFAlgo|createCandidatesHCAL")
2550 << "\t\t" << sortedTracks.size() << " tracks -> check if the excess is photonic or hadronic";
2551 #endif
2552
2553 double ratioMax = 0.;
2554 reco::PFClusterRef maxEcalRef;
2555 unsigned maxiEcal = 9999;
2556
2557
2558
2559 LogTrace("PFAlgo|createCandidatesHCAL") << "loop over sortedTracks.size()=" << sortedTracks.size();
2560 for (auto const& trk : sortedTracks) {
2561 unsigned iTrack = trk.second;
2562
2563 PFBlockElement::Type type = elements[iTrack].type();
2564 assert(type == reco::PFBlockElement::TRACK);
2565
2566 reco::TrackRef trackRef = elements[iTrack].trackRef();
2567 assert(!trackRef.isNull());
2568
2569 auto iae = associatedEcals.find(iTrack);
2570
2571 if (iae == associatedEcals.end())
2572 continue;
2573
2574
2575 unsigned iEcal = iae->second.second;
2576
2577 PFBlockElement::Type typeEcal = elements[iEcal].type();
2578 assert(typeEcal == reco::PFBlockElement::ECAL);
2579
2580 reco::PFClusterRef clusterRef = elements[iEcal].clusterRef();
2581 assert(!clusterRef.isNull());
2582
2583 double pTrack = trackRef->p();
2584 double eECAL = clusterRef->energy();
2585 double eECALOverpTrack = eECAL / pTrack;
2586
2587 if (eECALOverpTrack > ratioMax) {
2588 ratioMax = eECALOverpTrack;
2589 maxEcalRef = clusterRef;
2590 maxiEcal = iEcal;
2591 }
2592
2593 }
2594
2595 std::vector<reco::PFClusterRef> pivotalClusterRef;
2596 std::vector<unsigned> iPivotal;
2597 std::vector<double> particleEnergy, ecalEnergy, hcalEnergy, rawecalEnergy, rawhcalEnergy;
2598 std::vector<::math::XYZVector> particleDirection;
2599
2600
2601
2602 if (ePhoton < totalEcal || eNeutralHadron - calibEcal < 1E-10) {
2603 if (!maxEcalRef.isNull()) {
2604
2605 mergedPhotonEnergy = ePhoton;
2606 }
2607 } else {
2608
2609 if (!maxEcalRef.isNull()) {
2610
2611 mergedPhotonEnergy = totalEcalEGMCalib;
2612 }
2613
2614 mergedNeutralHadronEnergy = eNeutralHadron - calibEcal;
2615 }
2616
2617 if (mergedPhotonEnergy > 0) {
2618
2619
2620
2621 if (ecalClusters.size() <= 1) {
2622 ecalClusters.clear();
2623 ecalClusters.emplace_back(
2624 maxiEcal,
2625 photonAtECAL,
2626 1.);
2627 sumEcalClusters = sqrt(photonAtECAL.Mag2());
2628 }
2629 for (auto const& pae : ecalClusters) {
2630 const double clusterEnergyCalibrated =
2631 sqrt(std::get<1>(pae).Mag2()) *
2632 std::get<2>(
2633 pae);
2634 particleEnergy.push_back(mergedPhotonEnergy * clusterEnergyCalibrated / sumEcalClusters);
2635 particleDirection.push_back(std::get<1>(pae));
2636 ecalEnergy.push_back(mergedPhotonEnergy * clusterEnergyCalibrated / sumEcalClusters);
2637 hcalEnergy.push_back(0.);
2638 rawecalEnergy.push_back(totalEcal);
2639 rawhcalEnergy.push_back(0.);
2640 pivotalClusterRef.push_back(elements[std::get<0>(pae)].clusterRef());
2641 iPivotal.push_back(std::get<0>(pae));
2642 }
2643 }
2644
2645 if (mergedNeutralHadronEnergy > 1.0) {
2646
2647
2648 if (ecalClusters.size() <= 1) {
2649 ecalClusters.clear();
2650 ecalClusters.emplace_back(
2651 iHcal,
2652 hadronAtECAL,
2653 1.);
2654 sumEcalClusters = sqrt(hadronAtECAL.Mag2());
2655 }
2656 for (auto const& pae : ecalClusters) {
2657 const double clusterEnergyCalibrated =
2658 sqrt(std::get<1>(pae).Mag2()) *
2659 std::get<2>(
2660 pae);
2661 particleEnergy.push_back(mergedNeutralHadronEnergy * clusterEnergyCalibrated / sumEcalClusters);
2662 particleDirection.push_back(std::get<1>(pae));
2663 ecalEnergy.push_back(0.);
2664 hcalEnergy.push_back(mergedNeutralHadronEnergy * clusterEnergyCalibrated / sumEcalClusters);
2665 rawecalEnergy.push_back(0.);
2666 rawhcalEnergy.push_back(totalHcal);
2667 pivotalClusterRef.push_back(hclusterref);
2668 iPivotal.push_back(iHcal);
2669 }
2670 }
2671
2672
2673
2674
2675
2676 for (unsigned iPivot = 0; iPivot < iPivotal.size(); ++iPivot) {
2677 if (particleEnergy[iPivot] < 0.)
2678 edm::LogWarning("PFAlgo|createCandidatesHCAL")
2679 << "ALARM = Negative energy for iPivot=" << iPivot << ", " << particleEnergy[iPivot];
2680
2681 const bool useDirection = true;
2682 auto& neutral = (*pfCandidates_)[reconstructCluster(*pivotalClusterRef[iPivot],
2683 particleEnergy[iPivot],
2684 useDirection,
2685 particleDirection[iPivot].X(),
2686 particleDirection[iPivot].Y(),
2687 particleDirection[iPivot].Z())];
2688
2689 neutral.setEcalEnergy(rawecalEnergy[iPivot], ecalEnergy[iPivot]);
2690 if (!useHO_) {
2691 neutral.setHcalEnergy(rawhcalEnergy[iPivot], hcalEnergy[iPivot]);
2692 neutral.setHoEnergy(0., 0.);
2693 } else {
2694 if (rawhcalEnergy[iPivot] == 0.) {
2695 neutral.setHcalEnergy(0., 0.);
2696 neutral.setHoEnergy(0., 0.);
2697 } else {
2698 neutral.setHcalEnergy(max(rawhcalEnergy[iPivot] - totalHO, 0.0),
2699 hcalEnergy[iPivot] * max(1. - totalHO / rawhcalEnergy[iPivot], 0.));
2700 neutral.setHoEnergy(totalHO, totalHO * hcalEnergy[iPivot] / rawhcalEnergy[iPivot]);
2701 }
2702 }
2703 neutral.setPs1Energy(0.);
2704 neutral.setPs2Energy(0.);
2705 neutral.set_mva_nothing_gamma(-1.);
2706
2707
2708 neutral.addElementInBlock(blockref, iHcal);
2709 for (unsigned iTrack : chargedHadronsInBlock) {
2710 neutral.addElementInBlock(blockref, iTrack);
2711
2712 const ::math::XYZPointF& chargedPosition =
2713 dynamic_cast<const reco::PFBlockElementTrack*>(&elements[iTrack])->positionAtECALEntrance();
2714 neutral.setPositionAtECALEntrance(chargedPosition);
2715
2716 auto myEcals = associatedEcals.equal_range(iTrack);
2717 for (auto ii = myEcals.first; ii != myEcals.second; ++ii) {
2718 unsigned iEcal = ii->second.second;
2719 if (active[iEcal])
2720 continue;
2721 neutral.addElementInBlock(blockref, iEcal);
2722 }
2723 }
2724 }
2725
2726 }
2727
2728
2729
2730
2731
2732
2733
2734
2735
2736 double totalHcalEnergyCalibrated = std::max(calibHcal - mergedNeutralHadronEnergy, 0.);
2737
2738 double totalEcalEnergyCalibrated = std::max(calibEcal - mergedPhotonEnergy, 0.);
2739
2740
2741
2742
2743 double chargedHadronsTotalEnergy = 0;
2744 for (unsigned index : chargedHadronsIndices) {
2745 reco::PFCandidate& chargedHadron = (*pfCandidates_)[index];
2746 chargedHadronsTotalEnergy += chargedHadron.energy();
2747 }
2748
2749 for (unsigned index : chargedHadronsIndices) {
2750 reco::PFCandidate& chargedHadron = (*pfCandidates_)[index];
2751 float fraction = chargedHadron.energy() / chargedHadronsTotalEnergy;
2752
2753 if (!useHO_) {
2754 chargedHadron.setHcalEnergy(fraction * totalHcal, fraction * totalHcalEnergyCalibrated);
2755 chargedHadron.setHoEnergy(0., 0.);
2756 } else {
2757 chargedHadron.setHcalEnergy(fraction * max(totalHcal - totalHO, 0.0),
2758 fraction * totalHcalEnergyCalibrated * (1. - totalHO / totalHcal));
2759 chargedHadron.setHoEnergy(fraction * totalHO, fraction * totalHO * totalHcalEnergyCalibrated / totalHcal);
2760 }
2761
2762 chargedHadron.setEcalEnergy(fraction * totalEcal, fraction * totalEcalEnergyCalibrated);
2763 }
2764
2765
2766 for (auto const& ecalSatellite : ecalSatellites) {
2767
2768 unsigned iEcal = std::get<0>(ecalSatellite.second);
2769 if (!active[iEcal])
2770 continue;
2771
2772
2773 PFBlockElement::Type type = elements[iEcal].type();
2774 assert(type == PFBlockElement::ECAL);
2775 PFClusterRef eclusterref = elements[iEcal].clusterRef();
2776 assert(!eclusterref.isNull());
2777
2778
2779 active[iEcal] = false;
2780
2781
2782 std::multimap<double, unsigned> assTracks;
2783 block.associatedElements(iEcal, linkData, assTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
2784
2785
2786 double ecalClusterEnergyCalibrated =
2787 sqrt(std::get<1>(ecalSatellite.second).Mag2()) *
2788 std::get<2>(
2789 ecalSatellite.second);
2790 auto& cand = (*pfCandidates_)[reconstructCluster(*eclusterref, ecalClusterEnergyCalibrated)];
2791 cand.setEcalEnergy(eclusterref->energy(), ecalClusterEnergyCalibrated);
2792 cand.setHcalEnergy(0., 0.);
2793 cand.setHoEnergy(0., 0.);
2794 cand.setPs1Energy(associatedPSs[iEcal].first);
2795 cand.setPs2Energy(associatedPSs[iEcal].second);
2796 cand.addElementInBlock(blockref, iEcal);
2797 cand.addElementInBlock(blockref, sortedTracks.begin()->second);
2798
2799 if (fabs(eclusterref->energy() - sqrt(std::get<1>(ecalSatellite.second).Mag2())) > 1e-3 ||
2800 fabs(eclusterref->correctedEnergy() - ecalClusterEnergyCalibrated) > 1e-3)
2801 edm::LogWarning("PFAlgo:processBlock")
2802 << "ecalCluster vs ecalSatellites look inconsistent (eCluster E, calibE, ecalSatellite E, calib E): "
2803 << eclusterref->energy() << " " << eclusterref->correctedEnergy() << " "
2804 << sqrt(std::get<1>(ecalSatellite.second).Mag2()) << " " << ecalClusterEnergyCalibrated;
2805
2806 }
2807
2808 }
2809
2810 LogTrace("PFAlgo|createCandidatesHCAL") << "end of function PFAlgo::createCandidatesHCAL";
2811 }
2812
2813 void PFAlgo::createCandidatesHCALUnlinked(const reco::PFBlock& block,
2814 reco::PFBlock::LinkData& linkData,
2815 const edm::OwnVector<reco::PFBlockElement>& elements,
2816 std::vector<bool>& active,
2817 const reco::PFBlockRef& blockref,
2818 ElementIndices& inds,
2819 std::vector<bool>& deadArea) {
2820
2821 LogTrace("PFAlgo|createCandidatesHCALUnlinked")
2822 << "start of function PFAlgo::createCandidatesHCALUnlinked, hcalIs.size()=" << inds.hcalIs.size();
2823
2824
2825
2826 for (unsigned iHcal : inds.hcalIs) {
2827
2828 std::vector<unsigned> ecalRefs;
2829 std::vector<unsigned> hoRefs;
2830
2831 LogTrace("PFAlgo|createCandidatesHCALUnlinked") << elements[iHcal] << " ";
2832
2833 if (!active[iHcal]) {
2834 LogTrace("PFAlgo|createCandidatesHCALUnlinked") << "not active " << iHcal;
2835 continue;
2836 }
2837
2838
2839 std::multimap<double, unsigned> ecalElems;
2840 block.associatedElements(iHcal, linkData, ecalElems, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
2841
2842
2843 float totalEcal = 0.;
2844 float ecalMax = 0.;
2845 reco::PFClusterRef eClusterRef;
2846 for (auto const& ecal : ecalElems) {
2847 unsigned iEcal = ecal.second;
2848 double dist = ecal.first;
2849 PFBlockElement::Type type = elements[iEcal].type();
2850 assert(type == PFBlockElement::ECAL);
2851
2852
2853 if (!active[iEcal])
2854 continue;
2855
2856
2857
2858
2859
2860
2861
2862
2863
2864
2865
2866
2867
2868
2869
2870
2871
2872
2873
2874
2875 std::multimap<double, unsigned> hcalElems;
2876 block.associatedElements(iEcal, linkData, hcalElems, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
2877
2878 const bool isClosest = std::none_of(hcalElems.begin(), hcalElems.end(), [&](auto const& hcal) {
2879 return active[hcal.second] && hcal.first < dist;
2880 });
2881
2882 if (!isClosest)
2883 continue;
2884
2885 #ifdef EDM_ML_DEBUG
2886 LogTrace("PFAlgo|createCandidatesHCALUnlinked")
2887 << "\telement " << elements[iEcal] << " linked with dist " << dist;
2888 LogTrace("PFAlgo|createCandidatesHCALUnlinked") << "Added to HCAL cluster to form a neutral hadron";
2889 #endif
2890
2891 reco::PFClusterRef eclusterRef = elements[iEcal].clusterRef();
2892 assert(!eclusterRef.isNull());
2893
2894
2895 double ecalEnergy = eclusterRef->energy();
2896
2897 totalEcal += ecalEnergy;
2898 if (ecalEnergy > ecalMax) {
2899 ecalMax = ecalEnergy;
2900 eClusterRef = eclusterRef;
2901 }
2902
2903 ecalRefs.push_back(iEcal);
2904 active[iEcal] = false;
2905
2906 }
2907
2908
2909 double totalHO = 0.;
2910 double hoMax = 0.;
2911
2912 if (useHO_) {
2913 std::multimap<double, unsigned> hoElems;
2914 block.associatedElements(iHcal, linkData, hoElems, reco::PFBlockElement::HO, reco::PFBlock::LINKTEST_ALL);
2915
2916
2917
2918
2919
2920 reco::PFClusterRef hoClusterRef;
2921 for (auto const& ho : hoElems) {
2922 unsigned iHO = ho.second;
2923 double dist = ho.first;
2924 PFBlockElement::Type type = elements[iHO].type();
2925 assert(type == PFBlockElement::HO);
2926
2927
2928 if (!active[iHO])
2929 continue;
2930
2931
2932
2933
2934
2935 std::multimap<double, unsigned> hcalElems;
2936 block.associatedElements(iHO, linkData, hcalElems, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
2937
2938 const bool isClosest = std::none_of(hcalElems.begin(), hcalElems.end(), [&](auto const& hcal) {
2939 return active[hcal.second] && hcal.first < dist;
2940 });
2941
2942 if (!isClosest)
2943 continue;
2944
2945 #ifdef EDM_ML_DEBUG
2946 if (useHO_) {
2947 LogTrace("PFAlgo|createCandidatesHCALUnlinked")
2948 << "\telement " << elements[iHO] << " linked with dist " << dist;
2949 LogTrace("PFAlgo|createCandidatesHCALUnlinked") << "Added to HCAL cluster to form a neutral hadron";
2950 }
2951 #endif
2952
2953 reco::PFClusterRef hoclusterRef = elements[iHO].clusterRef();
2954 assert(!hoclusterRef.isNull());
2955
2956 double hoEnergy =
2957 hoclusterRef->energy();
2958
2959 totalHO += hoEnergy;
2960 if (hoEnergy > hoMax) {
2961 hoMax = hoEnergy;
2962 hoClusterRef = hoclusterRef;
2963
2964 }
2965
2966 hoRefs.push_back(iHO);
2967 active[iHO] = false;
2968
2969 }
2970 }
2971
2972 PFClusterRef hclusterRef = elements[iHcal].clusterRef();
2973 assert(!hclusterRef.isNull());
2974
2975
2976 double totalHcal = hclusterRef->energy();
2977
2978 if (useHO_)
2979 totalHcal += totalHO;
2980
2981
2982 double calibEcal = totalEcal > 0. ? totalEcal : 0.;
2983 double calibHcal = std::max(0., totalHcal);
2984 if (hclusterRef->layer() == PFLayer::HF_HAD || hclusterRef->layer() == PFLayer::HF_EM) {
2985 calibEcal = totalEcal;
2986 } else {
2987 calibration_.energyEmHad(
2988 -1., calibEcal, calibHcal, hclusterRef->positionREP().Eta(), hclusterRef->positionREP().Phi());
2989 }
2990
2991 auto& cand = (*pfCandidates_)[reconstructCluster(*hclusterRef, calibEcal + calibHcal)];
2992
2993 cand.setEcalEnergy(totalEcal, calibEcal);
2994 if (!useHO_) {
2995 cand.setHcalEnergy(totalHcal, calibHcal);
2996 cand.setHoEnergy(0., 0.);
2997 } else {
2998 cand.setHcalEnergy(max(totalHcal - totalHO, 0.0), calibHcal * (1. - totalHO / totalHcal));
2999 cand.setHoEnergy(totalHO, totalHO * calibHcal / totalHcal);
3000 }
3001 cand.setPs1Energy(0.);
3002 cand.setPs2Energy(0.);
3003 cand.addElementInBlock(blockref, iHcal);
3004 for (auto const& ec : ecalRefs)
3005 cand.addElementInBlock(blockref, ec);
3006 for (auto const& ho : hoRefs)
3007 cand.addElementInBlock(blockref, ho);
3008
3009 }
3010 }
3011
3012 void PFAlgo::createCandidatesECAL(const reco::PFBlock& block,
3013 reco::PFBlock::LinkData& linkData,
3014 const edm::OwnVector<reco::PFBlockElement>& elements,
3015 std::vector<bool>& active,
3016 const reco::PFBlockRef& blockref,
3017 ElementIndices& inds,
3018 std::vector<bool>& deadArea) {
3019 LogTrace("PFAlgo|createCandidatesECAL")
3020 << "start of function PFAlgo::createCandidatesECAL(), ecalIs.size()=" << inds.ecalIs.size();
3021
3022
3023
3024
3025
3026 for (unsigned i = 0; i < inds.ecalIs.size(); i++) {
3027 unsigned iEcal = inds.ecalIs[i];
3028
3029 LogTrace("PFAlgo|createCandidatesECAL") << "elements[" << iEcal << "]=" << elements[iEcal];
3030
3031 if (!active[iEcal]) {
3032 LogTrace("PFAlgo|createCandidatesECAL") << "iEcal=" << iEcal << " not active";
3033 continue;
3034 }
3035
3036 PFBlockElement::Type type = elements[iEcal].type();
3037 assert(type == PFBlockElement::ECAL);
3038
3039 PFClusterRef clusterref = elements[iEcal].clusterRef();
3040 assert(!clusterref.isNull());
3041
3042 active[iEcal] = false;
3043
3044 float ecalEnergy = clusterref->correctedEnergy();
3045
3046 double particleEnergy = ecalEnergy;
3047
3048 auto& cand = (*pfCandidates_)[reconstructCluster(*clusterref, particleEnergy)];
3049
3050 cand.setEcalEnergy(clusterref->energy(), ecalEnergy);
3051 cand.setHcalEnergy(0., 0.);
3052 cand.setHoEnergy(0., 0.);
3053 cand.setPs1Energy(0.);
3054 cand.setPs2Energy(0.);
3055 cand.addElementInBlock(blockref, iEcal);
3056
3057 }
3058 LogTrace("PFAlgo|createCandidatesECAL") << "end of function PFALgo::createCandidatesECAL";
3059 }
3060
3061 void PFAlgo::processBlock(const reco::PFBlockRef& blockref,
3062 std::list<reco::PFBlockRef>& hcalBlockRefs,
3063 std::list<reco::PFBlockRef>& ecalBlockRefs,
3064 PFEGammaFilters const* pfegamma) {
3065 assert(!blockref.isNull());
3066 const reco::PFBlock& block = *blockref;
3067
3068 LogTrace("PFAlgo|processBlock") << "start of function PFAlgo::processBlock, block=" << block;
3069
3070 const edm::OwnVector<reco::PFBlockElement>& elements = block.elements();
3071 LogTrace("PFAlgo|processBlock") << "elements.size()=" << elements.size();
3072
3073 PFBlock::LinkData linkData = block.linkData();
3074
3075
3076 vector<bool> active(elements.size(), true);
3077
3078
3079
3080 std::vector<reco::PFCandidate> tempElectronCandidates;
3081 tempElectronCandidates.clear();
3082
3083
3084 if (useEGammaFilters_) {
3085 egammaFilters(blockref, active, pfegamma);
3086 }
3087
3088
3089 if (usePFConversions_) {
3090 conversionAlgo(elements, active);
3091 }
3092
3093
3094
3095
3096
3097
3098
3099
3100
3101
3102
3103
3104
3105
3106
3107
3108
3109
3110
3111
3112
3113
3114 vector<bool> deadArea(elements.size(), false);
3115
3116
3117 ElementIndices inds;
3118
3119 elementLoop(block, linkData, elements, active, blockref, inds, deadArea);
3120
3121
3122
3123 if (!(inds.hfEmIs.empty() && inds.hfHadIs.empty())) {
3124 createCandidatesHF(block, linkData, elements, active, blockref, inds);
3125 if (inds.hcalIs.empty() && inds.ecalIs.empty())
3126 return;
3127 LogDebug("PFAlgo::processBlock")
3128 << "Block contains HF clusters, and also contains ECAL or HCAL clusters. Continue.\n"
3129 << block;
3130 }
3131
3132 createCandidatesHCAL(block, linkData, elements, active, blockref, inds, deadArea);
3133
3134 createCandidatesHCALUnlinked(block, linkData, elements, active, blockref, inds, deadArea);
3135 createCandidatesECAL(block, linkData, elements, active, blockref, inds, deadArea);
3136
3137 LogTrace("PFAlgo|processBlock") << "end of function PFAlgo::processBlock";
3138 }
3139
3140
3141 unsigned PFAlgo::reconstructTrack(const reco::PFBlockElement& elt, bool allowLoose) {
3142 const auto* eltTrack = dynamic_cast<const reco::PFBlockElementTrack*>(&elt);
3143
3144 const reco::TrackRef& trackRef = eltTrack->trackRef();
3145 const reco::Track& track = *trackRef;
3146 const reco::MuonRef& muonRef = eltTrack->muonRef();
3147 int charge = track.charge() > 0 ? 1 : -1;
3148
3149
3150 double px = track.px();
3151 double py = track.py();
3152 double pz = track.pz();
3153 double energy = sqrt(track.p() * track.p() + 0.13957 * 0.13957);
3154
3155 LogTrace("PFAlgo|reconstructTrack") << "Reconstructing PF candidate from track of pT = " << track.pt()
3156 << " eta = " << track.eta() << " phi = " << track.phi() << " px = " << px
3157 << " py = " << py << " pz = " << pz << " energy = " << energy;
3158
3159
3160 ::math::XYZTLorentzVector momentum(px, py, pz, energy);
3161 reco::PFCandidate::ParticleType particleType = reco::PFCandidate::h;
3162
3163
3164 LogTrace("PFAlgo|reconstructTrack") << "Creating PFCandidate charge=" << charge << ", type=" << particleType
3165 << ", pt=" << momentum.pt() << ", eta=" << momentum.eta()
3166 << ", phi=" << momentum.phi();
3167 pfCandidates_->push_back(PFCandidate(charge, momentum, particleType));
3168
3169 pfCandidates_->back().setVertex(trackRef->vertex());
3170 pfCandidates_->back().setTrackRef(trackRef);
3171 pfCandidates_->back().setPositionAtECALEntrance(eltTrack->positionAtECALEntrance());
3172 if (muonRef.isNonnull())
3173 pfCandidates_->back().setMuonRef(muonRef);
3174
3175
3176 if (elt.isTimeValid())
3177 pfCandidates_->back().setTime(elt.time(), elt.timeError());
3178
3179
3180 bool isMuon = pfmu_->reconstructMuon(pfCandidates_->back(), muonRef, allowLoose);
3181 bool isFromDisp = isFromSecInt(elt, "secondary");
3182
3183 if ((!isMuon) && isFromDisp) {
3184 double dpt = trackRef->ptError();
3185 double dptRel = dpt / trackRef->pt() * 100;
3186
3187
3188
3189 if (dptRel < dptRel_DispVtx_) {
3190 LogTrace("PFAlgo|reconstructTrack")
3191 << "Not refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy;
3192
3193 reco::PFDisplacedVertexRef vRef =
3194 eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef();
3195 reco::Track trackRefit = vRef->refittedTrack(trackRef);
3196
3197 ::math::XYZTLorentzVector momentum(
3198 trackRefit.px(), trackRefit.py(), trackRefit.pz(), sqrt(trackRefit.p() * trackRefit.p() + 0.13957 * 0.13957));
3199 LogTrace("PFAlgo|reconstructTrack")
3200 << "Refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy;
3201 }
3202 pfCandidates_->back().setFlag(reco::PFCandidate::T_FROM_DISP, true);
3203 pfCandidates_->back().setDisplacedVertexRef(
3204 eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef(),
3205 reco::PFCandidate::T_FROM_DISP);
3206 }
3207
3208
3209 if (isFromSecInt(elt, "primary") && !isMuon) {
3210 pfCandidates_->back().setFlag(reco::PFCandidate::T_TO_DISP, true);
3211 pfCandidates_->back().setDisplacedVertexRef(
3212 eltTrack->displacedVertexRef(reco::PFBlockElement::T_TO_DISP)->displacedVertexRef(),
3213 reco::PFCandidate::T_TO_DISP);
3214 }
3215
3216
3217 return pfCandidates_->size() - 1;
3218 }
3219
3220 unsigned PFAlgo::reconstructCluster(const reco::PFCluster& cluster,
3221 double particleEnergy,
3222 bool useDirection,
3223 double particleX,
3224 double particleY,
3225 double particleZ) {
3226 LogTrace("PFAlgo|reconstructCluster") << "start of function PFAlgo::reconstructCluster, cluster=" << cluster
3227 << "particleEnergy=" << particleEnergy << "useDirection=" << useDirection
3228 << "particleX=" << particleX << "particleY=" << particleY
3229 << "particleZ=" << particleZ;
3230
3231 reco::PFCandidate::ParticleType particleType = reco::PFCandidate::X;
3232
3233
3234
3235
3236
3237 double factor = 1.;
3238 if (useDirection) {
3239 switch (cluster.layer()) {
3240 case PFLayer::ECAL_BARREL:
3241 case PFLayer::HCAL_BARREL1:
3242 factor = std::sqrt(cluster.position().Perp2() / (particleX * particleX + particleY * particleY));
3243 break;
3244 case PFLayer::ECAL_ENDCAP:
3245 case PFLayer::HCAL_ENDCAP:
3246 case PFLayer::HF_HAD:
3247 case PFLayer::HF_EM:
3248 factor = cluster.position().Z() / particleZ;
3249 break;
3250 default:
3251 assert(0);
3252 }
3253 }
3254
3255 ::math::XYZPoint vertexPos;
3256 if (useVertices_)
3257 vertexPos = ::math::XYZPoint(primaryVertex_.x(), primaryVertex_.y(), primaryVertex_.z());
3258 else
3259 vertexPos = ::math::XYZPoint(0.0, 0.0, 0.0);
3260
3261 ::math::XYZVector clusterPos(cluster.position().X() - vertexPos.X(),
3262 cluster.position().Y() - vertexPos.Y(),
3263 cluster.position().Z() - vertexPos.Z());
3264 ::math::XYZVector particleDirection(
3265 particleX * factor - vertexPos.X(), particleY * factor - vertexPos.Y(), particleZ * factor - vertexPos.Z());
3266
3267
3268
3269
3270 clusterPos = useDirection ? particleDirection.Unit() : clusterPos.Unit();
3271 clusterPos *= particleEnergy;
3272
3273
3274
3275
3276 double mass = 0;
3277 ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> momentum(
3278 clusterPos.X(), clusterPos.Y(), clusterPos.Z(), mass);
3279
3280 ::math::XYZTLorentzVector tmp;
3281
3282 tmp = momentum;
3283
3284
3285 int charge = 0;
3286
3287
3288 switch (cluster.layer()) {
3289 case PFLayer::ECAL_BARREL:
3290 case PFLayer::ECAL_ENDCAP:
3291 particleType = PFCandidate::gamma;
3292 break;
3293 case PFLayer::HCAL_BARREL1:
3294 case PFLayer::HCAL_ENDCAP:
3295 particleType = PFCandidate::h0;
3296 break;
3297 case PFLayer::HF_HAD:
3298 particleType = PFCandidate::h_HF;
3299 break;
3300 case PFLayer::HF_EM:
3301 particleType = PFCandidate::egamma_HF;
3302 break;
3303 default:
3304 assert(0);
3305 }
3306
3307
3308 LogTrace("PFAlgo|reconstructCluster") << "Creating PFCandidate charge=" << charge << ", type=" << particleType
3309 << ", pt=" << tmp.pt() << ", eta=" << tmp.eta() << ", phi=" << tmp.phi();
3310 pfCandidates_->push_back(PFCandidate(charge, tmp, particleType));
3311
3312
3313
3314 pfCandidates_->back().setPositionAtECALEntrance(
3315 ::math::XYZPointF(cluster.position().X(), cluster.position().Y(), cluster.position().Z()));
3316
3317
3318 pfCandidates_->back().setVertex(vertexPos);
3319
3320
3321 setHcalDepthInfo(pfCandidates_->back(), cluster);
3322
3323
3324
3325 LogTrace("PFAlgo|reconstructCluster") << "** candidate: " << pfCandidates_->back();
3326
3327
3328 return pfCandidates_->size() - 1;
3329 }
3330
3331 void PFAlgo::setHcalDepthInfo(reco::PFCandidate& cand, const reco::PFCluster& cluster) const {
3332 std::array<double, 7> energyPerDepth;
3333 std::fill(energyPerDepth.begin(), energyPerDepth.end(), 0.0);
3334 for (auto& hitRefAndFrac : cluster.recHitFractions()) {
3335 const auto& hit = *hitRefAndFrac.recHitRef();
3336 if (DetId(hit.detId()).det() == DetId::Hcal) {
3337 if (hit.depth() == 0) {
3338 edm::LogWarning("setHcalDepthInfo") << "Depth zero found";
3339 continue;
3340 }
3341 if (hit.depth() < 1 || hit.depth() > 7) {
3342 throw cms::Exception("CorruptData") << "Bogus depth " << hit.depth() << " at detid " << hit.detId() << "\n";
3343 }
3344 energyPerDepth[hit.depth() - 1] += hitRefAndFrac.fraction() * hit.energy();
3345 }
3346 }
3347 double sum = std::accumulate(energyPerDepth.begin(), energyPerDepth.end(), 0.);
3348 std::array<float, 7> depthFractions;
3349 if (sum > 0) {
3350 for (unsigned int i = 0; i < depthFractions.size(); ++i) {
3351 depthFractions[i] = energyPerDepth[i] / sum;
3352 }
3353 } else {
3354 std::fill(depthFractions.begin(), depthFractions.end(), 0.f);
3355 }
3356 cand.setHcalDepthEnergyFractions(depthFractions);
3357 }
3358
3359
3360
3361 double PFAlgo::neutralHadronEnergyResolution(double clusterEnergyHCAL, double eta) const {
3362
3363 clusterEnergyHCAL = std::max(clusterEnergyHCAL, 1.);
3364
3365 double resol = fabs(eta) < 1.48 ? sqrt(1.02 * 1.02 / clusterEnergyHCAL + 0.065 * 0.065)
3366 : sqrt(1.20 * 1.20 / clusterEnergyHCAL + 0.028 * 0.028);
3367
3368 return resol;
3369 }
3370
3371 double PFAlgo::nSigmaHCAL(double clusterEnergyHCAL, double eta) const {
3372 double nS = fabs(eta) < 1.48 ? nSigmaHCAL_ * (1. + exp(-clusterEnergyHCAL / nSigmaEConstHCAL))
3373 : nSigmaHCAL_ * (1. + exp(-clusterEnergyHCAL / nSigmaEConstHCAL));
3374
3375 return nS;
3376 }
3377
3378 double PFAlgo::hfEnergyResolution(double clusterEnergyHF) const {
3379
3380 clusterEnergyHF = std::max(clusterEnergyHF, 1.);
3381
3382 double resol =
3383 sqrt(resolHF_square_[0] / clusterEnergyHF + resolHF_square_[1] + resolHF_square_[2] / pow(clusterEnergyHF, 2));
3384
3385
3386
3387 return resol;
3388 }
3389
3390 double PFAlgo::nSigmaHFEM(double clusterEnergyHF) const {
3391 double nS = nSigmaHFEM_ * (1. + exp(-clusterEnergyHF / nSigmaEConstHFEM));
3392 return nS;
3393 }
3394
3395 double PFAlgo::nSigmaHFHAD(double clusterEnergyHF) const {
3396 double nS = nSigmaHFHAD_ * (1. + exp(-clusterEnergyHF / nSigmaEConstHFHAD));
3397 return nS;
3398 }
3399
3400 ostream& operator<<(ostream& out, const PFAlgo& algo) {
3401 if (!out)
3402 return out;
3403
3404 out << "====== Particle Flow Algorithm ======= ";
3405 out << endl;
3406 out << "nSigmaECAL_ " << algo.nSigmaECAL_ << endl;
3407 out << "nSigmaHCAL_ " << algo.nSigmaHCAL_ << endl;
3408 out << "nSigmaHFEM_ " << algo.nSigmaHFEM_ << endl;
3409 out << "nSigmaHFHAD_ " << algo.nSigmaHFHAD_ << endl;
3410 out << endl;
3411 out << algo.calibration_ << endl;
3412 out << endl;
3413 out << "reconstructed particles: " << endl;
3414
3415 if (!algo.pfCandidates_.get()) {
3416 out << "candidates already transfered" << endl;
3417 return out;
3418 }
3419 for (auto const& c : *algo.pfCandidates_)
3420 out << c << endl;
3421
3422 return out;
3423 }
3424
3425 void PFAlgo::associatePSClusters(unsigned iEcal,
3426 reco::PFBlockElement::Type psElementType,
3427 const reco::PFBlock& block,
3428 const edm::OwnVector<reco::PFBlockElement>& elements,
3429 const reco::PFBlock::LinkData& linkData,
3430 std::vector<bool>& active,
3431 std::vector<double>& psEne) {
3432
3433
3434
3435
3436
3437
3438 std::multimap<double, unsigned> sortedPS;
3439 block.associatedElements(iEcal, linkData, sortedPS, psElementType, reco::PFBlock::LINKTEST_ALL);
3440
3441
3442 for (auto const& ps : sortedPS) {
3443
3444 unsigned iPS = ps.second;
3445
3446
3447
3448 if (!active[iPS])
3449 continue;
3450
3451
3452 std::multimap<double, unsigned> sortedECAL;
3453 block.associatedElements(iPS, linkData, sortedECAL, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
3454 unsigned jEcal = sortedECAL.begin()->second;
3455 if (jEcal != iEcal)
3456 continue;
3457
3458
3459 PFBlockElement::Type pstype = elements[iPS].type();
3460 assert(pstype == psElementType);
3461 PFClusterRef psclusterref = elements[iPS].clusterRef();
3462 assert(!psclusterref.isNull());
3463 psEne[0] += psclusterref->energy();
3464 active[iPS] = false;
3465 }
3466 }
3467
3468 bool PFAlgo::isFromSecInt(const reco::PFBlockElement& eTrack, string order) const {
3469 reco::PFBlockElement::TrackType T_TO_DISP = reco::PFBlockElement::T_TO_DISP;
3470 reco::PFBlockElement::TrackType T_FROM_DISP = reco::PFBlockElement::T_FROM_DISP;
3471
3472 reco::PFBlockElement::TrackType T_FROM_V0 = reco::PFBlockElement::T_FROM_V0;
3473
3474 bool bPrimary = (order.find("primary") != string::npos);
3475 bool bSecondary = (order.find("secondary") != string::npos);
3476 bool bAll = (order.find("all") != string::npos);
3477
3478 bool isToDisp = usePFNuclearInteractions_ && eTrack.trackType(T_TO_DISP);
3479 bool isFromDisp = usePFNuclearInteractions_ && eTrack.trackType(T_FROM_DISP);
3480
3481 if (bPrimary && isToDisp)
3482 return true;
3483 if (bSecondary && isFromDisp)
3484 return true;
3485 if (bAll && (isToDisp || isFromDisp))
3486 return true;
3487
3488
3489
3490
3491
3492 bool isFromDecay = (bAll || bSecondary) && usePFDecays_ && eTrack.trackType(T_FROM_V0);
3493
3494 return isFromDecay;
3495 }
3496
3497 void PFAlgo::postCleaning() {
3498
3499 double metX = 0.;
3500 double metY = 0.;
3501 double sumet = 0;
3502 std::vector<unsigned int> pfCandidatesToBeRemoved;
3503 for (auto const& pfc : *pfCandidates_) {
3504 metX += pfc.px();
3505 metY += pfc.py();
3506 sumet += pfc.pt();
3507 }
3508 double met2 = metX * metX + metY * metY;
3509
3510 double significance = std::sqrt(met2 / sumet);
3511 double significanceCor = significance;
3512 if (significance > minSignificance_) {
3513 double metXCor = metX;
3514 double metYCor = metY;
3515 double sumetCor = sumet;
3516 double met2Cor = met2;
3517 double deltaPhi = 3.14159;
3518 double deltaPhiPt = 100.;
3519 bool next = true;
3520 unsigned iCor = 1E9;
3521
3522
3523 while (next) {
3524 double metReduc = -1.;
3525
3526 for (unsigned i = 0; i < pfCandidates_->size(); ++i) {
3527 const PFCandidate& pfc = (*pfCandidates_)[i];
3528
3529
3530 if (pfc.particleId() != reco::PFCandidate::h_HF && pfc.particleId() != reco::PFCandidate::egamma_HF)
3531 continue;
3532
3533
3534 if (pfc.pt() < minHFCleaningPt_)
3535 continue;
3536
3537
3538 const bool skip = std::any_of(
3539 pfCandidatesToBeRemoved.begin(), pfCandidatesToBeRemoved.end(), [&](unsigned int j) { return i == j; });
3540 if (skip)
3541 continue;
3542
3543
3544 deltaPhi = std::acos((metX * pfc.px() + metY * pfc.py()) / (pfc.pt() * std::sqrt(met2)));
3545 deltaPhiPt = deltaPhi * pfc.pt();
3546 if (deltaPhiPt > maxDeltaPhiPt_)
3547 continue;
3548
3549
3550 double metXInt = metX - pfc.px();
3551 double metYInt = metY - pfc.py();
3552 double sumetInt = sumet - pfc.pt();
3553 double met2Int = metXInt * metXInt + metYInt * metYInt;
3554 if (met2Int < met2Cor) {
3555 metXCor = metXInt;
3556 metYCor = metYInt;
3557 metReduc = (met2 - met2Int) / met2Int;
3558 met2Cor = met2Int;
3559 sumetCor = sumetInt;
3560 significanceCor = std::sqrt(met2Cor / sumetCor);
3561 iCor = i;
3562 }
3563 }
3564
3565
3566 if (metReduc > minDeltaMet_) {
3567 pfCandidatesToBeRemoved.push_back(iCor);
3568 metX = metXCor;
3569 metY = metYCor;
3570 sumet = sumetCor;
3571 met2 = met2Cor;
3572 } else {
3573
3574 next = false;
3575 }
3576 }
3577
3578
3579 if (significance - significanceCor > minSignificanceReduction_ && significanceCor < maxSignificance_) {
3580 edm::LogInfo("PFAlgo|postCleaning") << "Significance reduction = " << significance << " -> " << significanceCor
3581 << " = " << significanceCor - significance;
3582 for (unsigned int toRemove : pfCandidatesToBeRemoved) {
3583 edm::LogInfo("PFAlgo|postCleaning") << "Removed : " << (*pfCandidates_)[toRemove];
3584 pfCleanedCandidates_.push_back((*pfCandidates_)[toRemove]);
3585 (*pfCandidates_)[toRemove].rescaleMomentum(1E-6);
3586
3587
3588 }
3589 }
3590 }
3591 }
3592
3593 void PFAlgo::checkCleaning(const reco::PFRecHitCollection& cleanedHits) {
3594
3595 if (cleanedHits.empty())
3596 return;
3597
3598
3599 double metX = 0.;
3600 double metY = 0.;
3601 double sumet = 0;
3602 std::vector<unsigned int> hitsToBeAdded;
3603 for (auto const& pfc : *pfCandidates_) {
3604 metX += pfc.px();
3605 metY += pfc.py();
3606 sumet += pfc.pt();
3607 }
3608 double met2 = metX * metX + metY * metY;
3609 double met2_Original = met2;
3610
3611
3612
3613 double metXCor = metX;
3614 double metYCor = metY;
3615 double sumetCor = sumet;
3616 double met2Cor = met2;
3617 bool next = true;
3618 unsigned iCor = 1E9;
3619
3620
3621 while (next) {
3622 double metReduc = -1.;
3623
3624 for (unsigned i = 0; i < cleanedHits.size(); ++i) {
3625 const PFRecHit& hit = cleanedHits[i];
3626 double length = std::sqrt(hit.position().mag2());
3627 double px = hit.energy() * hit.position().x() / length;
3628 double py = hit.energy() * hit.position().y() / length;
3629 double pt = std::sqrt(px * px + py * py);
3630
3631
3632 bool skip = false;
3633 for (unsigned int hitIdx : hitsToBeAdded) {
3634 if (i == hitIdx)
3635 skip = true;
3636 if (skip)
3637 break;
3638 }
3639 if (skip)
3640 continue;
3641
3642
3643 double metXInt = metX + px;
3644 double metYInt = metY + py;
3645 double sumetInt = sumet + pt;
3646 double met2Int = metXInt * metXInt + metYInt * metYInt;
3647
3648
3649 if (met2Int < met2Cor) {
3650 metXCor = metXInt;
3651 metYCor = metYInt;
3652 metReduc = (met2 - met2Int) / met2Int;
3653 met2Cor = met2Int;
3654 sumetCor = sumetInt;
3655
3656 iCor = i;
3657 }
3658 }
3659
3660
3661
3662 if (metReduc > minDeltaMet_) {
3663 hitsToBeAdded.push_back(iCor);
3664 metX = metXCor;
3665 metY = metYCor;
3666 sumet = sumetCor;
3667 met2 = met2Cor;
3668 } else {
3669
3670 next = false;
3671 }
3672 }
3673
3674
3675 if (std::sqrt(met2_Original) - std::sqrt(met2) > 5.) {
3676 LogTrace("PFAlgo|checkCleaning") << hitsToBeAdded.size() << " hits were re-added ";
3677 LogTrace("PFAlgo|checkCleaning") << "MET reduction = " << std::sqrt(met2_Original) << " -> " << std::sqrt(met2Cor)
3678 << " = " << std::sqrt(met2Cor) - std::sqrt(met2_Original);
3679 LogTrace("PFAlgo|checkCleaning") << "Added after cleaning check : ";
3680 for (unsigned int hitIdx : hitsToBeAdded) {
3681 const PFRecHit& hit = cleanedHits[hitIdx];
3682 PFCluster cluster(hit.layer(), hit.energy(), hit.position().x(), hit.position().y(), hit.position().z());
3683 reconstructCluster(cluster, hit.energy());
3684 LogTrace("PFAlgo|checkCleaning") << pfCandidates_->back() << ". time = " << hit.time();
3685 }
3686 }
3687 }