PFProducer/src/PFAlgo.cc

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   // HF resolution parameters
0037   assert(resolHF_square_.size() == 3);  // make sure that stochastic, constant, noise (i.e. three) terms are specified.
0038
0039   // Muon parameters
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   // Bad Hcal Track Parameters
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   //Set the vertices for muon cleaning
0114   pfmu_->setInputsForCleaning(primaryVertices);
0115
0116   //Now find the primary vertex!
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   //Use vertices if the user wants to but only if it exists a good vertex
0127   useVertices_ = useVertex && primaryVertexFound;
0128 }
0129
0130 void PFAlgo::reconstructParticles(const reco::PFBlockHandle& blockHandle, PFEGammaFilters const* pfegamma) {
0131   auto const& blocks = *blockHandle;
0132
0133   // reset output collection
0134   pfCandidates_->clear();
0135
0136   LogTrace("PFAlgo|reconstructParticles")
0137       << "start of function PFAlgo::reconstructParticles, blocks.size()=" << blocks.size();
0138
0139   // sort elements in three lists:
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         // Single HO elements are likely to be noise. Not considered for now.
0163         hoBlockRefs.push_back(blockref);
0164         singleEcalOrHcal = true;
0165       }
0166     }
0167
0168     if (!singleEcalOrHcal) {
0169       otherBlockRefs.push_back(blockref);
0170     }
0171   }  //loop blocks
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   // loop on blocks that are not single ecal,
0178   // and not single hcal.
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   // process remaining single hcal blocks
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   // process remaining single ecal blocks
0197   for (auto const& ecal : ecalBlockRefs) {
0198     LogTrace("PFAlgo|reconstructParticles") << "processBlock, ECAL block number " << eblcks++;
0199     processBlock(ecal, empty, empty, pfegamma);
0200   }
0201
0202   // Post HF Cleaning
0203   pfCleanedCandidates_.clear();
0204   // Check if the post HF Cleaning was requested - if not, do nothing
0205   if (postHFCleaning_) {
0206     postCleaning();
0207   }
0208
0209   //Muon post cleaning
0210   pfmu_->postClean(pfCandidates_.get());
0211
0212   //Add Missing muons
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   // const edm::ValueMap<reco::GsfElectronRef> & myGedElectronValMap(*valueMapGedElectrons_);
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     //      const reco::PFCandidate & egmcand((*pfEgammaCandidates_)[ieg]);
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     }  // end sameBlock
0269
0270     if (isGoodElectron && isGoodPhoton) {
0271       if (isPrimaryElectron)
0272         isGoodPhoton = false;
0273       else
0274         isGoodElectron = false;
0275     }
0276
0277     // isElectron
0278     if (isGoodElectron) {
0279       reco::GsfElectronRef gedEleRef = (*valueMapGedElectrons_)[pfEgmRef];
0280       reco::PFCandidate myPFElectron = *pfEgmRef;
0281       // run protections
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         // The electron is considered safe for JetMET and the additional tracks pointing to it are locked
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     }  //isGoodElectron
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         // DNN pfid
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         // Lock all the elements
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       }  // end isSafe
0365     }    // end isGoodPhoton
0366   }      // end loop on EGM candidates
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   // vector<unsigned> elementIndices;
0413   // elementIndices.push_back(iTrack);
0414
0415   // The track momentum
0416   double trackMomentum = elements[iTrack].trackRef()->p();
0417   LogTrace("PFAlgo|recoTracksNotHCAL") << elements[iTrack];
0418
0419   // Is it a "tight" muon ? In that case assume no
0420   //Track momentum corresponds to the calorimeter
0421   //energy for now
0422   bool thisIsAMuon = PFMuonAlgo::isMuon(elements[iTrack]);
0423   if (thisIsAMuon)
0424     trackMomentum = 0.;
0425
0426   // If it is not a muon check if Is it a fake track ?
0427   //Michalis: I wonder if we should convert this to dpt/pt?
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   }  // !thisIsaMuon
0469
0470   // Create a track candidate
0471   // unsigned tmpi = reconstructTrack( elements[iTrack] );
0472   //active[iTrack] = false;
0473   std::vector<unsigned> tmpi;
0474   std::vector<unsigned> kTrack;
0475
0476   // Some cleaning : secondary tracks without calo's and large momentum must be fake
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   }  //rejectTracks_Step45_ && ...
0499
0500   tmpi.push_back(reconstructTrack(elements[iTrack]));
0501
0502   kTrack.push_back(iTrack);
0503   active[iTrack] = false;
0504
0505   // No ECAL cluster either ... continue...
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   // Look for closest ECAL cluster
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   // Set ECAL energy for muons
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   // Consider charged particles closest to the same ECAL cluster
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     // GP: only allow the ecal cluster closest to this track to be considered
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     // Don't consider already used tracks
0572     if (!active[jTrack])
0573       continue;
0574
0575     // The loop is on the other tracks !
0576     if (jTrack == iTrack)
0577       continue;
0578
0579     // Check if the ECAL cluster closest to this track is
0580     // indeed the current ECAL cluster. Only tracks not
0581     // linked to an HCAL cluster are considered here
0582     // (GP: this is because if there's a jTrack linked
0583     // to the same Ecal cluster and that also has an Hcal link
0584     // we would have also linked iTrack to that Hcal above)
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     // Check if this track is a muon
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     // Check if this track is not a fake
0596     bool rejectFake = false;
0597     reco::TrackRef trackRef = elements[jTrack].trackRef();
0598     if (!thatIsAMuon && trackRef->ptError() > ptError_) {
0599       // GP: FIXME this selection should be adapted depending on hasDeadHcal
0600       //     but we don't know if jTrack is linked to a dead Hcal or not
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     // Otherwise, add this track momentum to the total track momentum
0622     /* */
0623     // reco::TrackRef trackRef = elements[jTrack].trackRef();
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     // And create a charged particle candidate !
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   // Loop over all ECAL linked clusters ordered by increasing distance.
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     // Just skip clusters already taken
0660     if (!active[index]) {
0661       LogTrace("PFAlgo|recoTracksNotHCAL") << "is not active  - ignore ";
0662       continue;
0663     }
0664
0665     // Just skip this cluster if it's closer to another track
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     // The corresponding PFCluster and energy
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     // Check the presence of preshower clusters in the vicinity
0684     // Preshower cluster closer to another ECAL cluster are ignored.
0685     //JOSH: This should use the association map already produced by the cluster corrector for consistency,
0686     //but should be ok for now
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     // KH: use raw ECAL energy for PF hadron calibration. use calibrated ECAL energy when adding PF photons
0693     const double ecalEnergy = clusterRef->energy();
0694     const double ecalEnergyCalibrated = clusterRef->correctedEnergy();  // calibrated based on the egamma hypothesis
0695     LogTrace("PFAlgo|recoTracksNotHCAL") << "Corrected ECAL(+PS) energy = " << ecalEnergy;
0696
0697     // Since the electrons were found beforehand, this track must be a hadron. Calibrate
0698     // the energy under the hadron hypothesis.
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     // Stop the loop when adding more ECAL clusters ruins the compatibility
0713     if (connectedToEcal && calibEcal - trackMomentum >= 0.) {
0714       // if ( connectedToEcal && calibEcal - trackMomentum >=
0715       //     nSigmaECAL_*neutralHadronEnergyResolution(trackMomentum,clusterRef->positionREP().Eta())  ) {
0716       calibEcal = previousCalibEcal;
0717       slopeEcal = previousSlopeEcal;
0718       totalEcal = calibEcal / slopeEcal;
0719
0720       // Turn this last cluster in a photon
0721       // (The PS clusters are already locked in "associatePSClusters")
0722       active[index] = false;
0723
0724       // Find the associated tracks
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)];  // KH: use the PF ECAL cluster calibrated energy
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       // Check that there is at least one track
0737       if (!assTracks.empty()) {
0738         ecalCand.addElementInBlock(blockref, assTracks.begin()->second);
0739
0740         // Assign the position of the track at the ECAL entrance
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     // Lock used clusters.
0750     connectedToEcal = true;
0751     iEcal = index;
0752     active[index] = false;
0753     for (unsigned ic : tmpi)
0754       (*pfCandidates_)[ic].addElementInBlock(blockref, iEcal);
0755
0756   }  // Loop ecal elements
0757
0758   bool bNeutralProduced = false;
0759
0760   // Create a photon if the ecal energy is too large
0761   if (connectedToEcal) {
0762     reco::PFClusterRef pivotalRef = elements[iEcal].clusterRef();
0763
0764     double neutralEnergy = calibEcal - trackMomentum;
0765
0766     /*
0767     // Check if there are other tracks linked to that ECAL
0768     for(IE ie = sortedTracks.begin(); ie != sortedTracks.end() && neutralEnergy > 0; ++ie ) {
0769       unsigned jTrack = ie->second;
0770
0771       // The loop is on the other tracks !
0772       if ( jTrack == iTrack ) continue;
0773
0774       // Check if the ECAL closest to this track is the current ECAL
0775       // Otherwise ignore this track in the neutral energy determination
0776       std::multimap<double, unsigned> sortedECAL;
0777       block.associatedElements( jTrack,  linkData,
0778                     sortedECAL,
0779                     reco::PFBlockElement::ECAL,
0780                     reco::PFBlock::LINKTEST_ALL );
0781       if ( sortedECAL.begin()->second != iEcal ) continue;
0782
0783       // Check if this track is also linked to an HCAL
0784       // (in which case it goes to the next loop and is ignored here)
0785       std::multimap<double, unsigned> sortedHCAL;
0786       block.associatedElements( jTrack,  linkData,
0787                     sortedHCAL,
0788                     reco::PFBlockElement::HCAL,
0789                     reco::PFBlock::LINKTEST_ALL );
0790       if ( sortedHCAL.size() ) continue;
0791
0792       // Otherwise, subtract this track momentum from the ECAL energy
0793       reco::TrackRef trackRef = elements[jTrack].trackRef();
0794       neutralEnergy -= trackRef->p();
0795
0796     } // End other tracks
0797     */
0798
0799     // Add a photon if the energy excess is large enough
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     }  // End neutral energy
0815
0816     // Set elements in blocks and ECAL energies to all tracks
0817     for (unsigned ic = 0; ic < tmpi.size(); ++ic) {
0818       // Skip muons
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   }  // End connected ECAL
0838
0839   // Fill the element_in_block for tracks that are eventually linked to no ECAL clusters at all.
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 //Check if the track is a primary track of a secondary interaction
0853 //If that is the case reconstruct a charged hadron only using that
0854 //track
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   // there's 3 possible options possible here, in principle:
0877   //    1) flag everything that may be associated to a dead hcal marker
0878   //    2) flag everything whose closest hcal link is a dead hcal marker
0879   //    3) flag only things that are linked only to dead hcal marker
0880   // in our first test we go for (2)
0881   //--- option (1) --
0882   //bool hasDeadHcal = false;
0883   //for (auto it = hcalElems.begin(), ed = hcalElems.end(); it != ed; /*NOTE NO ++it HERE */ ) {
0884   //    if (deadArea[it->second]) { hasDeadHcal = true; it = hcalElems.erase(it); } // std::multimap::erase returns iterator to next
0885   //    else ++it;
0886   //}
0887   //--- option (2) --
0888   bool hasDeadHcal = false;
0889   if (!hcalElems.empty() && deadArea[hcalElems.begin()->second]) {
0890     hasDeadHcal = true;
0891   }
0892   //--- option (3) --
0893   //bool hasDeadHcal = true;
0894   //for (auto it = hcalElems.begin(), ed = hcalElems.end(); it != ed; /*NOTE NO ++it HERE */ ) {
0895   //    if (deadArea[it->second]) { it = hcalElems.erase(it); } // std::multimap::erase returns iterator to next
0896   //    else { hasDeadHcal = false; }
0897   //}
0898   return hasDeadHcal;
0899 }
0900
0901 // for tracks with bad Hcal, check the track quality
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     // now we add an extra block for tracks at high |eta|
0911     if (!goodTrackDeadHcal && std::abs(trackRef->eta()) > goodPixelTrackDeadHcal_minEta_ &&  // high eta
0912         trackRef->hitPattern().pixelLayersWithMeasurement() >= 3 &&                          // pixel track
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_ &&  // tighter cut
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() &&  // sanity
0922         trackRef->pt() < goodPixelTrackDeadHcal_maxPt_) {                           // sanity
0923       goodTrackDeadHcal = true;
0924       // FIXME: may decide to do something to the track pT
0925     }
0926     //if (!goodTrackDeadHcal && trackRef->hitPattern().trackerLayersWithMeasurement() == 4 && trackRef->validFraction() == 1
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   // Loop over all tracks
0956   for (auto const& trk : sortedTracks) {
0957     unsigned jTrack = trk.second;
0958     LogTrace("PFAlgo|elementLoop") << "jTrack=" << jTrack;
0959     // Track must be active
0960     if (!active[jTrack])
0961       continue;
0962     LogTrace("PFAlgo|elementLoop") << "active[jTrack]=" << active[jTrack];
0963
0964     // The loop is on the other tracks !
0965     if (jTrack == iTrack)
0966       continue;
0967     LogTrace("PFAlgo|elementLoop") << "skipping jTrack=" << jTrack << " for same iTrack";
0968
0969     // Check if the ECAL closest to this track is the current ECAL
0970     // Otherwise ignore this track in the neutral energy determination
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     // Check if this track is also linked to an HCAL
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     // In that case establish a link with the first track
0986     block.setLink(iTrack, sortedHCAL.begin()->second, sortedECAL.begin()->first, linkData, PFBlock::LINKTEST_RECHIT);
0987
0988   }  // End other tracks
0989
0990   // Redefine HCAL elements
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     //only process TRACK elements, but fill the ElementIndices vector with indices for all elements.
1011     //Mark the active & deadArea for bad HCAL
1012     auto ret_decideType = decideType(elements, type, active, inds, deadArea, iEle);
1013     if (ret_decideType == 1) {
1014       LogTrace("PFAlgo|elementLoop") << "ret_decideType==1, continuing";
1015       continue;
1016     }
1017     LogTrace("PFAlgo|elementLoop") << "ret_decideType=" << ret_decideType << " type=" << type;
1018
1019     active[iEle] = checkAndReconstructSecondaryInteraction(blockref, elements, active[iEle], iEle);
1020
1021     if (!active[iEle]) {
1022       LogTrace("PFAlgo|elementLoop") << "Already used by electrons, muons, conversions";
1023       continue;
1024     }
1025
1026     reco::TrackRef trackRef = elements[iEle].trackRef();
1027     assert(!trackRef.isNull());
1028
1029     LogTrace("PFAlgo|elementLoop") << "PFAlgo:processBlock"
1030                                    << " trackIs.size()=" << inds.trackIs.size()
1031                                    << " ecalIs.size()=" << inds.ecalIs.size() << " hcalIs.size()=" << inds.hcalIs.size()
1032                                    << " hoIs.size()=" << inds.hoIs.size() << " hfEmIs.size()=" << inds.hfEmIs.size()
1033                                    << " hfHadIs.size()=" << inds.hfHadIs.size();
1034
1035     // look for associated elements of all types
1036     //COLINFEB16
1037     // all types of links are considered.
1038     // the elements are sorted by increasing distance
1039     std::multimap<double, unsigned> ecalElems;
1040     block.associatedElements(iEle, linkData, ecalElems, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
1041
1042     std::multimap<double, unsigned> hcalElems;
1043     block.associatedElements(iEle, linkData, hcalElems, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
1044
1045     std::multimap<double, unsigned> hfEmElems;
1046     std::multimap<double, unsigned> hfHadElems;
1047     block.associatedElements(iEle, linkData, hfEmElems, reco::PFBlockElement::HFEM, reco::PFBlock::LINKTEST_ALL);
1048     block.associatedElements(iEle, linkData, hfHadElems, reco::PFBlockElement::HFHAD, reco::PFBlock::LINKTEST_ALL);
1049
1050     LogTrace("PFAlgo|elementLoop") << "\tTrack " << iEle << " is linked to " << ecalElems.size() << " ecal and "
1051                                    << hcalElems.size() << " hcal and " << hfEmElems.size() << " hfEm and "
1052                                    << hfHadElems.size() << " hfHad elements";
1053
1054 #ifdef EDM_ML_DEBUG
1055     for (const auto& pair : ecalElems) {
1056       LogTrace("PFAlgo|elementLoop") << "ecal: dist " << pair.first << "\t elem " << pair.second;
1057     }
1058     for (const auto& pair : hcalElems) {
1059       LogTrace("PFAlgo|elementLoop") << "hcal: dist " << pair.first << "\t elem " << pair.second
1060                                      << (deadArea[pair.second] ? "  DEAD AREA MARKER" : "");
1061     }
1062 #endif
1063
1064     const bool hasDeadHcal = checkHasDeadHcal(hcalElems, deadArea);
1065     if (hasDeadHcal) {
1066       hcalElems.clear();
1067     }
1068     const bool goodTrackDeadHcal = checkGoodTrackDeadHcal(trackRef, hasDeadHcal);
1069
1070     // When a track has no HCAL cluster linked, but another track is linked to the same
1071     // ECAL cluster and an HCAL cluster, link the track to the HCAL cluster for
1072     // later analysis
1073     if (hcalElems.empty() && !ecalElems.empty() && !hasDeadHcal) {
1074       relinkTrackToHcal(block, ecalElems, hcalElems, active, linkData, iEle);
1075     }
1076
1077     //MICHELE
1078     //TEMPORARY SOLUTION FOR ELECTRON REJECTION IN PFTAU
1079     //COLINFEB16
1080     // in case particle flow electrons are not reconstructed,
1081     // the mva_e_pi of the charged hadron will be set to 1
1082     // if a GSF element is associated to the current TRACK element
1083     // This information will be used in the electron rejection for tau ID.
1084     std::multimap<double, unsigned> gsfElems;
1085     block.associatedElements(iEle, linkData, gsfElems, reco::PFBlockElement::GSF);
1086
1087     if (hcalElems.empty()) {
1088       LogTrace("PFAlgo|elementLoop") << "no hcal element connected to track " << iEle;
1089     }
1090
1091     // will now loop on associated elements ...
1092     bool hcalFound = false;
1093     bool hfhadFound = false;
1094
1095     LogTrace("PFAlgo|elementLoop") << "now looping on elements associated to the track: ecalElems";
1096
1097     // ... first on associated ECAL elements
1098     // Check if there is still a free ECAL for this track
1099     for (auto const& ecal : ecalElems) {
1100       unsigned index = ecal.second;
1101       // Sanity checks and optional printout
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       // This ECAL is not free (taken by an electron?) - just skip it
1112       if (!active[index])
1113         continue;
1114
1115       // Flag ECAL clusters for which the corresponding track is not linked to an HCAL cluster
1116
1117       //reco::PFClusterRef clusterRef = elements[index].clusterRef();
1118       //assert( !clusterRef.isNull() );
1119       if (!hcalElems.empty())
1120         LogTrace("PFAlgo|elementLoop") << "\t\tat least one hcal element connected to the track."
1121                                        << " Sparing Ecal cluster for the hcal loop";
1122
1123     }  //loop print ecal elements
1124
1125     // tracks which are not linked to an HCAL (or HFHAD)
1126     // are reconstructed now.
1127
1128     if (hcalElems.empty() && hfHadElems.empty()) {
1129       auto ret_continue = recoTracksNotHCAL(
1130           block, linkData, elements, blockref, active, goodTrackDeadHcal, hasDeadHcal, iEle, ecalElems, trackRef);
1131       if (ret_continue) {
1132         continue;
1133       }
1134     }  // end if( hcalElems.empty() && hfHadElems.empty() )
1135
1136     // In case several HCAL elements are linked to this track,
1137     // unlinking all of them except the closest.
1138     for (auto const& hcal : hcalElems) {
1139       unsigned index = hcal.second;
1140
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       // all hcal clusters except the closest
1150       // will be unlinked from the track
1151       if (!hcalFound) {  // closest hcal
1152         LogTrace("PFAlgo|elementLoop") << "\t\tclosest hcal cluster, doing nothing";
1153
1154         hcalFound = true;
1155
1156         // active[index] = false;
1157         // hcalUsed.push_back( index );
1158       } else {  // other associated hcal
1159         // unlink from the track
1160         LogTrace("PFAlgo|elementLoop") << "\t\tsecondary hcal cluster. unlinking";
1161         block.setLink(iEle, index, -1., linkData, PFBlock::LINKTEST_RECHIT);
1162       }
1163     }  //loop hcal elements
1164
1165     // ---Same for HFHAD---
1166     // In case several HFHAD elements are linked to this track,
1167     // unlinking all of them except the closest.
1168     for (auto const& hfhad : hfHadElems) {
1169       unsigned index = hfhad.second;
1170
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       // all hfhad clusters except the closest
1180       // will be unlinked from the track
1181       if (!hfhadFound) {  // closest hfhad
1182         LogTrace("PFAlgo|elementLoop") << "\t\tclosest hfhad cluster, doing nothing";
1183
1184         hfhadFound = true;
1185
1186       } else {  // other associated hfhad
1187         // unlink from the track
1188         LogTrace("PFAlgo|elementLoop") << "\t\tsecondary hfhad cluster. unlinking";
1189         block.setLink(iEle, index, -1., linkData, PFBlock::LINKTEST_RECHIT);
1190       }
1191     }  //loop hfhad elements
1192
1193     LogTrace("PFAlgo|elementLoop") << "end of loop over iEle";
1194   }  // end of outer loop on elements iEle of any type
1195   LogTrace("PFAlgo|elementLoop") << "end of function PFAlgo::elementLoop";
1196 }
1197
1198 //Arranges the PFBlock elements according to type into the ElementIndices output vector.
1199 //Also checks for dead HCAL area and updates the active and deadArea vectors.
1200 //Returns 0 for elements of TRACK type, 1 otherwise
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;  //continue
1220     case PFBlockElement::HCAL:
1221       if (active[iEle]) {
1222         if (elements[iEle].clusterRef()->flags() & reco::CaloCluster::badHcalMarker) {
1223           LogTrace("PFAlgo|decideType") << "HCAL DEAD AREA: remember and skip.";
1224           active[iEle] = false;
1225           deadArea[iEle] = true;
1226           return 1;  //continue
1227         }
1228         inds.hcalIs.push_back(iEle);
1229         LogTrace("PFAlgo|decideType") << "HCAL, stored index, continue";
1230       }
1231       return 1;  //continue
1232     case PFBlockElement::HO:
1233       if (useHO_) {
1234         if (active[iEle]) {
1235           inds.hoIs.push_back(iEle);
1236           LogTrace("PFAlgo|decideType") << "HO, stored index, continue";
1237         }
1238       }
1239       return 1;  //continue
1240     case PFBlockElement::HFEM:
1241       if (active[iEle]) {
1242         inds.hfEmIs.push_back(iEle);
1243         LogTrace("PFAlgo|decideType") << "HFEM, stored index, continue";
1244       }
1245       return 1;  //continue
1246     case PFBlockElement::HFHAD:
1247       if (active[iEle]) {
1248         inds.hfHadIs.push_back(iEle);
1249         LogTrace("PFAlgo|decideType") << "HFHAD, stored index, continue";
1250       }
1251       return 1;  //continue
1252     default:
1253       return 1;  //continue
1254   }
1255   LogTrace("PFAlgo|decideType") << "Did not match type to anything, return 0";
1256   return 0;
1257 }
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   // inds.trackIs can be empty, even if there are tracks in this block,
1269   // but what we want to check is if this block has any track including inactive ones
1270   if (!trackInBlock)
1271     for (unsigned iEle = 0; iEle < elements.size(); iEle++) {
1272       PFBlockElement::Type type = elements[iEle].type();
1273       if (type == PFBlockElement::TRACK) {
1274         trackInBlock = true;
1275         break;
1276       }
1277     }
1278   // there is at least one HF element in this block.
1279   // in case of no track, all elements must be HF
1280   if (!trackInBlock)
1281     assert(inds.hfEmIs.size() + inds.hfHadIs.size() == elements.size());
1282
1283   //
1284   // Dealing with a block with at least one track
1285   // Occasionally, there are only inactive tracks and multiple HF clusters. Consider such blocks too.
1286   //
1287   if (trackInBlock) {  // count any tracks (not only active tracks)
1288     // sorted tracks associated with a HfHad cluster
1289     std::multimap<double, unsigned> sortedTracks;
1290     std::multimap<double, unsigned> sortedTracksActive;  // only active ones
1291     // HfEms associated with tracks linked to a HfHad cluster
1292     std::multimap<unsigned, std::pair<double, unsigned>> associatedHfEms;
1293     // Temporary map for HfEm satellite clusters
1294     std::multimap<double, std::pair<unsigned, double>> hfemSatellites;
1295
1296     //
1297     // Loop over active HfHad clusters
1298     //
1299     for (unsigned iHfHad : inds.hfHadIs) {
1300       PFBlockElement::Type type = elements[iHfHad].type();
1301       assert(type == PFBlockElement::HFHAD);
1302
1303       PFClusterRef hclusterRef = elements[iHfHad].clusterRef();
1304       assert(!hclusterRef.isNull());
1305
1306       sortedTracks.clear();
1307       sortedTracksActive.clear();
1308       associatedHfEms.clear();
1309       hfemSatellites.clear();
1310
1311       // Look for associated tracks
1312       block.associatedElements(
1313           iHfHad, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
1314
1315       LogTrace("PFAlgo|createCandidatesHF") << "elements[" << iHfHad << "]=" << elements[iHfHad];
1316
1317       if (sortedTracks.empty()) {
1318         LogTrace("PFAlgo|createCandidatesHCF") << "\tno associated tracks, keep for later";
1319         continue;
1320       }
1321
1322       // Lock the HFHAD cluster
1323       active[iHfHad] = false;
1324
1325       LogTrace("PFAlgo|createCandidatesHF") << sortedTracks.size() << " associated tracks:";
1326
1327       double totalChargedMomentum = 0.;
1328       double sumpError2 = 0.;
1329
1330       //
1331       // Loop over all tracks associated to this HFHAD cluster
1332       //
1333       for (auto const& trk : sortedTracks) {
1334         unsigned iTrack = trk.second;
1335
1336         // Check the track has not already been used
1337         if (!active[iTrack])
1338           continue;
1339         // Sanity check 1
1340         PFBlockElement::Type type = elements[iTrack].type();
1341         assert(type == reco::PFBlockElement::TRACK);
1342         // Sanity check 2
1343         reco::TrackRef trackRef = elements[iTrack].trackRef();
1344         assert(!trackRef.isNull());
1345
1346         // Introduce tracking errors
1347         double trackMomentum = trackRef->p();
1348         totalChargedMomentum += trackMomentum;
1349
1350         // Also keep the total track momentum error for comparison with the calo energy
1351         double dp = trackRef->qoverpError() * trackMomentum * trackMomentum;
1352         sumpError2 += dp * dp;
1353
1354         // Store active tracks for 2nd loop to create charged hadrons
1355         sortedTracksActive.emplace(trk);
1356
1357         // look for HFEM elements associated to iTrack (associated to iHfHad)
1358         std::multimap<double, unsigned> sortedHfEms;
1359         block.associatedElements(
1360             iTrack, linkData, sortedHfEms, reco::PFBlockElement::HFEM, reco::PFBlock::LINKTEST_ALL);
1361
1362         LogTrace("PFAlgo|createCandidatesHF") << "number of HfEm elements linked to this track: " << sortedHfEms.size();
1363
1364         bool connectedToHfEm = false;  // Will become true if there is at least one HFEM cluster connected
1365
1366         //
1367         // Loop over all HFEM clusters connected to iTrack
1368         //
1369         for (auto const& hfem : sortedHfEms) {
1370           unsigned iHfEm = hfem.second;
1371           double dist = hfem.first;
1372
1373           // Ignore HFEM cluters already used
1374           if (!active[iHfEm]) {
1375             LogTrace("PFAlgo|createCandidatesHF") << "cluster locked";
1376             continue;
1377           }
1378
1379           // Sanity checks
1380           PFBlockElement::Type type = elements[iHfEm].type();
1381           assert(type == PFBlockElement::HFEM);
1382           PFClusterRef eclusterRef = elements[iHfEm].clusterRef();
1383           assert(!eclusterRef.isNull());
1384
1385           // Check if this HFEM is not closer to another track - ignore it in that case
1386           std::multimap<double, unsigned> sortedTracksHfEm;
1387           block.associatedElements(
1388               iHfEm, linkData, sortedTracksHfEm, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
1389           unsigned jTrack = sortedTracksHfEm.begin()->second;
1390           if (jTrack != iTrack)
1391             continue;
1392
1393           double distHfEm = block.dist(jTrack, iHfEm, linkData, reco::PFBlock::LINKTEST_ALL);
1394           double hfemEnergy = eclusterRef->energy();
1395
1396           if (!connectedToHfEm) {  // This is the closest HFEM cluster - will add its energy later
1397
1398             LogTrace("PFAlgo|createCandidatesHF") << "closest: " << elements[iHfEm];
1399             connectedToHfEm = true;
1400             std::pair<unsigned, double> satellite(iHfEm, hfemEnergy);
1401             hfemSatellites.emplace(-1., satellite);
1402
1403           } else {  // Keep satellite clusters for later
1404
1405             // KH: same as above.
1406             std::pair<unsigned, double> satellite(iHfEm, hfemEnergy);
1407             hfemSatellites.emplace(dist, satellite);
1408           }
1409
1410           std::pair<double, unsigned> associatedHfEm(distHfEm, iHfEm);
1411           associatedHfEms.emplace(iTrack, associatedHfEm);
1412
1413         }  // End loop hfem associated to iTrack
1414       }    // sortedTracks
1415
1416       // HfHad energy
1417       double uncalibratedenergyHfHad = hclusterRef->energy();
1418       double energyHfHad = uncalibratedenergyHfHad;
1419       if (thepfEnergyCalibrationHF_.getcalibHF_use()) {
1420         energyHfHad = thepfEnergyCalibrationHF_.energyHad(  // HAD only calibration
1421             uncalibratedenergyHfHad,
1422             hclusterRef->positionREP().Eta(),
1423             hclusterRef->positionREP().Phi());
1424       }
1425       double calibFactorHfHad = (uncalibratedenergyHfHad > 0.) ? energyHfHad / uncalibratedenergyHfHad : 1.;
1426
1427       // HfEm energy
1428       double energyHfEmTmp = 0.;
1429       double uncalibratedenergyHfEmTmp = 0.;
1430       double energyHfEm = 0.;
1431       double uncalibratedenergyHfEm = 0.;
1432
1433       // estimated HF resolution and track p error
1434       double caloResolution = hfEnergyResolution(totalChargedMomentum);
1435       caloResolution *= totalChargedMomentum;
1436       double totalError = sqrt(caloResolution * caloResolution + sumpError2);
1437       double nsigmaHFEM = nSigmaHFEM(totalChargedMomentum);
1438       double nsigmaHFHAD = nSigmaHFHAD(totalChargedMomentum);
1439
1440       // Handle case that no active track gets associated to HfHad cluster
1441       if (sortedTracksActive.empty()) {
1442         // look for HFEM elements associated to iHfHad
1443         std::multimap<double, unsigned> sortedHfEms;
1444         std::multimap<double, unsigned> sortedHfEmsActive;
1445         block.associatedElements(
1446             iHfHad, linkData, sortedHfEms, reco::PFBlockElement::HFEM, reco::PFBlock::LINKTEST_ALL);
1447         //
1448         // If iHfHad is connected to HFEM cluster, Loop over all of them
1449         //
1450         if (!sortedHfEms.empty()) {
1451           for (auto const& hfem : sortedHfEms) {
1452             unsigned iHfEm = hfem.second;
1453             // Ignore HFEM cluters already used
1454             if (!active[iHfEm])
1455               continue;
1456             sortedHfEmsActive.emplace(hfem);
1457             PFClusterRef eclusterRef = elements[iHfEm].clusterRef();
1458             assert(!eclusterRef.isNull());
1459             double hfemEnergy = eclusterRef->energy();
1460             uncalibratedenergyHfEm += hfemEnergy;
1461             energyHfEm = uncalibratedenergyHfEm;
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             }  // calib true
1468           }    // loop over sortedHfEm
1469         }      // if !sortedHfEms.empty()
1470         //
1471         // Create HF candidates
1472         unsigned tmpi = reconstructCluster(*hclusterRef, energyHfEm + energyHfHad);
1473         (*pfCandidates_)[tmpi].setHcalEnergy(uncalibratedenergyHfHad, energyHfHad);
1474         (*pfCandidates_)[tmpi].setEcalEnergy(uncalibratedenergyHfEm, energyHfEm);
1475         (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfHad);
1476         for (auto const& hfem : sortedHfEmsActive) {
1477           unsigned iHfEm = hfem.second;
1478           (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfEm);
1479           active[iHfEm] = false;
1480         }
1481
1482       }  // if sortedTracksActive.empty() ends
1483       //
1484       // Active tracks are associated.
1485       // Create HFHAD candidates from excess energy w.r.t. tracks
1486       else if ((energyHfHad - totalChargedMomentum) > nsigmaHFHAD * totalError) {  // HfHad is excessive
1487         assert(energyHfEm == 0.);
1488         // HfHad candidate from excess
1489         double energyHfHadExcess = max(energyHfHad - totalChargedMomentum, 0.);
1490         double uncalibratedenergyHfHadExcess = energyHfHadExcess / calibFactorHfHad;
1491         unsigned tmpi = reconstructCluster(*hclusterRef, energyHfHadExcess);
1492         (*pfCandidates_)[tmpi].setHcalEnergy(uncalibratedenergyHfHadExcess, energyHfHadExcess);
1493         (*pfCandidates_)[tmpi].setEcalEnergy(0., 0.);
1494         (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfHad);
1495         energyHfHad = max(energyHfHad - energyHfHadExcess, 0.);
1496         uncalibratedenergyHfHad = max(uncalibratedenergyHfHad - uncalibratedenergyHfHadExcess, 0.);
1497       }
1498       //
1499       // If there is a room for HFEM satellites to get associated,
1500       // loop over all HFEM satellites, starting for the closest to the various tracks
1501       // and adding other satellites until saturation of the total track momentum
1502       //
1503       else {
1504         for (auto const& hfemSatellite : hfemSatellites) {
1505           //
1506           uncalibratedenergyHfEmTmp += std::get<1>(hfemSatellite.second);  // KH: raw HFEM energy
1507           energyHfEmTmp = uncalibratedenergyHfEmTmp;
1508           double energyHfHadTmp = uncalibratedenergyHfHad;  // now to test hfhad calibration with EM+HAD cases
1509           unsigned iHfEm = std::get<0>(hfemSatellite.second);
1510           PFClusterRef eclusterRef = elements[iHfEm].clusterRef();
1511           assert(!eclusterRef.isNull());
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           // Continue looping until all closest clusters are exhausted and as long as
1523           // the calorimetric energy does not saturate the total momentum.
1524           if (hfemSatellite.first < 0. || caloEnergyTmp < totalChargedMomentum) {
1525             LogTrace("PFAlgo|createCandidatesHF")
1526                 << "\t\t\tactive, adding " << std::get<1>(hfemSatellite.second) << " to HFEM energy, and locking";
1527             active[std::get<0>(hfemSatellite.second)] = false;
1528             // HfEm is excessive (possible for the first hfemSatellite)
1529             if (hfemSatellite.first < 0. && (caloEnergyTmp - totalChargedMomentum) > nsigmaHFEM * totalError) {
1530               // HfEm candidate from excess
1531               double energyHfEmExcess = max(caloEnergyTmp - totalChargedMomentum, 0.);
1532               double uncalibratedenergyHfEmExcess = energyHfEmExcess / calibFactorHfEm;
1533               unsigned tmpi = reconstructCluster(*eclusterRef, energyHfEmExcess);
1534               (*pfCandidates_)[tmpi].setEcalEnergy(uncalibratedenergyHfEmExcess, energyHfEmExcess);
1535               (*pfCandidates_)[tmpi].setHcalEnergy(0, 0.);
1536               (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfEm);
1537               energyHfEmTmp = max(energyHfEmTmp - energyHfEmExcess, 0.);
1538               uncalibratedenergyHfEmTmp = max(uncalibratedenergyHfEmTmp - uncalibratedenergyHfEmExcess, 0.);
1539             }
1540             energyHfEm = energyHfEmTmp;
1541             uncalibratedenergyHfEm = uncalibratedenergyHfEmTmp;
1542             energyHfHad = energyHfHadTmp;
1543             continue;
1544           }
1545           break;
1546         }  // loop over hfemsattelites ends
1547       }    // if HFHAD is excessive or not
1548
1549       //
1550       // Loop over all tracks associated to this HFHAD cluster *again* in order to produce charged hadrons
1551       //
1552       for (auto const& trk : sortedTracksActive) {
1553         unsigned iTrack = trk.second;
1554
1555         // Sanity check
1556         reco::TrackRef trackRef = elements[iTrack].trackRef();
1557         assert(!trackRef.isNull());
1558
1559         //
1560         // Reconstructing charged hadrons
1561         //
1562         unsigned tmpi = reconstructTrack(elements[iTrack]);
1563         active[iTrack] = false;
1564         (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfHad);
1565         auto myHfEms = associatedHfEms.equal_range(iTrack);
1566         for (auto ii = myHfEms.first; ii != myHfEms.second; ++ii) {
1567           unsigned iHfEm = ii->second.second;
1568           if (active[iHfEm])
1569             continue;
1570           (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHfEm);
1571         }
1572         double frac = 0.;
1573         if (totalChargedMomentum)
1574           frac = trackRef->p() / totalChargedMomentum;
1575         (*pfCandidates_)[tmpi].setEcalEnergy(uncalibratedenergyHfEm * frac, energyHfEm * frac);
1576         (*pfCandidates_)[tmpi].setHcalEnergy(uncalibratedenergyHfHad * frac, energyHfHad * frac);
1577
1578       }  // sortedTracks loop ends
1579
1580     }  // iHfHad element loop ends
1581
1582     //
1583     // Loop over remaining active HfEm clusters
1584     //
1585     for (unsigned iHfEm = 0; iHfEm < elements.size(); iHfEm++) {
1586       PFBlockElement::Type type = elements[iHfEm].type();
1587       if (type == PFBlockElement::HFEM && active[iHfEm]) {
1588         reco::PFClusterRef eclusterRef = elements[iHfEm].clusterRef();
1589         double energyHF = 0.;
1590         double uncalibratedenergyHF = 0.;
1591         unsigned tmpi = 0;
1592         // do EM-only calibration here
1593         energyHF = eclusterRef->energy();
1594         uncalibratedenergyHF = energyHF;
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     }  // remaining active HfEm cluster loop ends
1607
1608   }  // if-statement for blocks including tracks ends here
1609   //
1610   // -----------------------------------------------
1611   // From here, traditional PF HF candidate creation
1612   // -----------------------------------------------
1613   //
1614   else if (elements.size() == 1) {
1615     //Auguste: HAD-only calibration here
1616     reco::PFClusterRef clusterRef = elements[0].clusterRef();
1617     double energyHF = 0.;
1618     double uncalibratedenergyHF = 0.;
1619     unsigned tmpi = 0;
1620     switch (clusterRef->layer()) {
1621       case PFLayer::HF_EM:
1622         // do EM-only calibration here
1623         energyHF = clusterRef->energy();
1624         uncalibratedenergyHF = energyHF;
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         // do HAD-only calibration here
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     //Auguste: EM + HAD calibration here
1660     reco::PFClusterRef c0 = elements[0].clusterRef();
1661     reco::PFClusterRef c1 = elements[1].clusterRef();
1662     // 2 HF elements. Must be in each layer.
1663     reco::PFClusterRef cem = (c0->layer() == PFLayer::HF_EM ? c0 : c1);
1664     reco::PFClusterRef chad = (c1->layer() == PFLayer::HF_HAD ? c1 : c0);
1665
1666     if (cem->layer() != PFLayer::HF_EM || chad->layer() != PFLayer::HF_HAD) {
1667       edm::LogError("PFAlgo::createCandidatesHF") << "Error: 2 elements, but not 1 HFEM and 1 HFHAD";
1668       edm::LogError("PFAlgo::createCandidatesHF") << block;
1669       assert(0);
1670       //    assert( c1->layer()== PFLayer::HF_EM &&
1671       //        c0->layer()== PFLayer::HF_HAD );
1672     }
1673     // do EM+HAD calibration here
1674     double energyHfEm = cem->energy();
1675     double energyHfHad = chad->energy();
1676     double uncalibratedenergyHfEm = energyHfEm;
1677     double uncalibratedenergyHfHad = energyHfHad;
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     // Unusual blocks including HF elements, but do not fit any of the above categories
1695     edm::LogWarning("PFAlgo::createCandidatesHF")
1696         << "Warning: HF, but n elem different from 1 or 2 or >=3 or !trackIs.empty()";
1697     edm::LogWarning("PFAlgo::createCandidatesHF") << block;
1698   }
1699   LogTrace("PFAlgo|createCandidateHF") << "end of function PFAlgo::createCandidateHF";
1700 }
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   // --------------- loop hcal ------------------
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     // associated tracks
1722     std::multimap<double, unsigned> sortedTracks;
1723     block.associatedElements(iHcal, linkData, sortedTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
1724
1725     std::multimap<unsigned, std::pair<double, unsigned>> associatedEcals;
1726
1727     std::map<unsigned, std::pair<double, double>> associatedPSs;
1728
1729     std::multimap<double, std::pair<unsigned, bool>> associatedTracks;
1730
1731     // A temporary maps for ECAL satellite clusters
1732     std::multimap<double, std::tuple<unsigned, ::math::XYZVector, double>>
1733         ecalSatellites;  // last element (double) : correction under the egamma hypothesis
1734     std::tuple<unsigned, ::math::XYZVector, double> fakeSatellite(iHcal, ::math::XYZVector(0., 0., 0.), 1.);
1735     ecalSatellites.emplace(-1., fakeSatellite);
1736
1737     std::multimap<unsigned, std::pair<double, unsigned>> associatedHOs;
1738
1739     PFClusterRef hclusterref = elements[iHcal].clusterRef();
1740     assert(!hclusterref.isNull());
1741
1742     //if there is no track attached to that HCAL, then do not
1743     //reconstruct an HCAL alone, check if it can be recovered
1744     //first
1745
1746     // if no associated tracks, create a neutral hadron
1747     //if(sortedTracks.empty() ) {
1748
1749     if (sortedTracks.empty()) {
1750       LogTrace("PFAlgo|createCandidatesHCAL") << "\tno associated tracks, keep for later";
1751       continue;
1752     }
1753
1754     // Lock the HCAL cluster
1755     active[iHcal] = false;
1756
1757     // in the following, tracks are associated to this hcal cluster.
1758     // will look for an excess of energy in the calorimeters w/r to
1759     // the charged energy, and turn this excess into a neutral hadron or
1760     // a photon.
1761
1762     LogTrace("PFAlgo|createCandidatesHCAL") << sortedTracks.size() << " associated tracks:";
1763
1764     double totalChargedMomentum = 0;
1765     double sumpError2 = 0.;
1766     double totalHO = 0.;
1767     double totalEcal = 0.;
1768     double totalEcalEGMCalib = 0.;
1769     double totalHcal = hclusterref->energy();
1770     vector<double> hcalP;
1771     vector<double> hcalDP;
1772     vector<unsigned> tkIs;
1773     double maxDPovP = -9999.;
1774
1775     //Keep track of how much energy is assigned to calorimeter-vs-track energy/momentum excess
1776     vector<unsigned> chargedHadronsIndices;
1777     vector<unsigned> chargedHadronsInBlock;
1778     double mergedNeutralHadronEnergy = 0;
1779     double mergedPhotonEnergy = 0;
1780     double muonHCALEnergy = 0.;
1781     double muonECALEnergy = 0.;
1782     double muonHCALError = 0.;
1783     double muonECALError = 0.;
1784     unsigned nMuons = 0;
1785
1786     ::math::XYZVector photonAtECAL(0., 0., 0.);
1787     std::vector<std::tuple<unsigned, ::math::XYZVector, double>>
1788         ecalClusters;  // last element (double) : correction under the egamma hypothesis
1789     double sumEcalClusters = 0;
1790     ::math::XYZVector hadronDirection(
1791         hclusterref->position().X(), hclusterref->position().Y(), hclusterref->position().Z());
1792     hadronDirection = hadronDirection.Unit();
1793     ::math::XYZVector hadronAtECAL = totalHcal * hadronDirection;
1794
1795     // Loop over all tracks associated to this HCAL cluster
1796     for (auto const& trk : sortedTracks) {
1797       unsigned iTrack = trk.second;
1798
1799       // Check the track has not already been used (e.g., in electrons, conversions...)
1800       if (!active[iTrack])
1801         continue;
1802       // Sanity check 1
1803       PFBlockElement::Type type = elements[iTrack].type();
1804       assert(type == reco::PFBlockElement::TRACK);
1805       // Sanity check 2
1806       reco::TrackRef trackRef = elements[iTrack].trackRef();
1807       assert(!trackRef.isNull());
1808
1809       // The direction at ECAL entrance
1810       const ::math::XYZPointF& chargedPosition =
1811           dynamic_cast<const reco::PFBlockElementTrack*>(&elements[iTrack])->positionAtECALEntrance();
1812       ::math::XYZVector chargedDirection(chargedPosition.X(), chargedPosition.Y(), chargedPosition.Z());
1813       chargedDirection = chargedDirection.Unit();
1814
1815       // look for ECAL elements associated to iTrack (associated to iHcal)
1816       std::multimap<double, unsigned> sortedEcals;
1817       block.associatedElements(iTrack, linkData, sortedEcals, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
1818
1819       LogTrace("PFAlgo|createCandidatesHCAL") << "number of Ecal elements linked to this track: " << sortedEcals.size();
1820
1821       // look for HO elements associated to iTrack (associated to iHcal)
1822       std::multimap<double, unsigned> sortedHOs;
1823       if (useHO_) {
1824         block.associatedElements(iTrack, linkData, sortedHOs, reco::PFBlockElement::HO, reco::PFBlock::LINKTEST_ALL);
1825       }
1826       LogTrace("PFAlgo|createCandidatesHCAL")
1827           << "PFAlgo : number of HO elements linked to this track: " << sortedHOs.size();
1828
1829       // Create a PF Candidate right away if the track is a tight muon
1830       reco::MuonRef muonRef = elements[iTrack].muonRef();
1831
1832       bool thisIsAMuon = PFMuonAlgo::isMuon(elements[iTrack]);
1833       bool thisIsAnIsolatedMuon = PFMuonAlgo::isIsolatedMuon(elements[iTrack]);
1834       bool thisIsALooseMuon = false;
1835
1836       if (!thisIsAMuon) {
1837         thisIsALooseMuon = PFMuonAlgo::isLooseMuon(elements[iTrack]);
1838       }
1839
1840       if (thisIsAMuon) {
1841         LogTrace("PFAlgo|createCandidatesHCAL") << "This track is identified as a muon - remove it from the stack";
1842         LogTrace("PFAlgo|createCandidatesHCAL") << elements[iTrack];
1843
1844         // Create a muon.
1845
1846         unsigned tmpi = reconstructTrack(elements[iTrack]);
1847
1848         (*pfCandidates_)[tmpi].addElementInBlock(blockref, iTrack);
1849         (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHcal);
1850         double muonHcal = std::min(muonHCAL_[0] + muonHCAL_[1], totalHcal);
1851
1852         // if muon is isolated and muon momentum exceeds the calo energy, absorb the calo energy
1853         // rationale : there has been a EM showering by the muon in the calorimeter - or the coil -
1854         // and we don't want to double count the energy
1855         bool letMuonEatCaloEnergy = false;
1856
1857         if (thisIsAnIsolatedMuon) {
1858           // The factor 1.3 is the e/pi factor in HCAL...
1859           double totalCaloEnergy = totalHcal / 1.30;
1860           unsigned iEcal = 0;
1861           if (!sortedEcals.empty()) {
1862             iEcal = sortedEcals.begin()->second;
1863             PFClusterRef eclusterref = elements[iEcal].clusterRef();
1864             totalCaloEnergy += eclusterref->energy();
1865           }
1866
1867           if (useHO_) {
1868             // The factor 1.3 is assumed to be the e/pi factor in HO, too.
1869             unsigned iHO = 0;
1870             if (!sortedHOs.empty()) {
1871               iHO = sortedHOs.begin()->second;
1872               PFClusterRef eclusterref = elements[iHO].clusterRef();
1873               totalCaloEnergy += eclusterref->energy() / 1.30;
1874             }
1875           }
1876
1877           if ((pfCandidates_->back()).p() > totalCaloEnergy)
1878             letMuonEatCaloEnergy = true;
1879         }
1880
1881         if (letMuonEatCaloEnergy)
1882           muonHcal = totalHcal;
1883         double muonEcal = 0.;
1884         unsigned iEcal = 0;
1885         if (!sortedEcals.empty()) {
1886           iEcal = sortedEcals.begin()->second;
1887           PFClusterRef eclusterref = elements[iEcal].clusterRef();
1888           (*pfCandidates_)[tmpi].addElementInBlock(blockref, iEcal);
1889           muonEcal = std::min(muonECAL_[0] + muonECAL_[1], eclusterref->energy());
1890           if (letMuonEatCaloEnergy)
1891             muonEcal = eclusterref->energy();
1892           // If the muon expected energy accounts for the whole ecal cluster energy, lock the ecal cluster
1893           if (eclusterref->energy() - muonEcal < 0.2)
1894             active[iEcal] = false;
1895           (*pfCandidates_)[tmpi].setEcalEnergy(eclusterref->energy(), muonEcal);
1896         }
1897         unsigned iHO = 0;
1898         double muonHO = 0.;
1899         if (useHO_) {
1900           if (!sortedHOs.empty()) {
1901             iHO = sortedHOs.begin()->second;
1902             PFClusterRef hoclusterref = elements[iHO].clusterRef();
1903             (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHO);
1904             muonHO = std::min(muonHO_[0] + muonHO_[1], hoclusterref->energy());
1905             if (letMuonEatCaloEnergy)
1906               muonHO = hoclusterref->energy();
1907             // If the muon expected energy accounts for the whole HO cluster energy, lock the HO cluster
1908             if (hoclusterref->energy() - muonHO < 0.2)
1909               active[iHO] = false;
1910             (*pfCandidates_)[tmpi].setHcalEnergy(totalHcal, muonHcal);
1911             (*pfCandidates_)[tmpi].setHoEnergy(hoclusterref->energy(), muonHO);
1912           }
1913         } else {
1914           (*pfCandidates_)[tmpi].setHcalEnergy(totalHcal, muonHcal);
1915         }
1916         setHcalDepthInfo((*pfCandidates_)[tmpi], *hclusterref);
1917
1918         if (letMuonEatCaloEnergy) {
1919           muonHCALEnergy += totalHcal;
1920           if (useHO_)
1921             muonHCALEnergy += muonHO;
1922           muonHCALError += 0.;
1923           muonECALEnergy += muonEcal;
1924           muonECALError += 0.;
1925           photonAtECAL -= muonEcal * chargedDirection;
1926           hadronAtECAL -= totalHcal * chargedDirection;
1927           if (!sortedEcals.empty())
1928             active[iEcal] = false;
1929           active[iHcal] = false;
1930           if (useHO_ && !sortedHOs.empty())
1931             active[iHO] = false;
1932         } else {
1933           // Estimate of the energy deposit & resolution in the calorimeters
1934           muonHCALEnergy += muonHCAL_[0];
1935           muonHCALError += muonHCAL_[1] * muonHCAL_[1];
1936           if (muonHO > 0.) {
1937             muonHCALEnergy += muonHO_[0];
1938             muonHCALError += muonHO_[1] * muonHO_[1];
1939           }
1940           muonECALEnergy += muonECAL_[0];
1941           muonECALError += muonECAL_[1] * muonECAL_[1];
1942           // ... as well as the equivalent "momentum" at ECAL entrance
1943           photonAtECAL -= muonECAL_[0] * chargedDirection;
1944           hadronAtECAL -= muonHCAL_[0] * chargedDirection;
1945         }
1946
1947         // Remove it from the stack
1948         active[iTrack] = false;
1949         // Go to next track
1950         continue;
1951       }
1952
1953       //
1954
1955       LogTrace("PFAlgo|createCandidatesHCAL") << "elements[iTrack=" << iTrack << "]=" << elements[iTrack];
1956
1957       // introduce tracking errors
1958       double trackMomentum = trackRef->p();
1959       totalChargedMomentum += trackMomentum;
1960
1961       // If the track is not a tight muon, but still resembles a muon
1962       // keep it for later for blocks with too large a charged energy
1963       if (thisIsALooseMuon && !thisIsAMuon)
1964         nMuons += 1;
1965
1966       // ... and keep anyway the pt error for possible fake rejection
1967       // ... blow up errors of 5th and 4th iteration, to reject those
1968       // ... tracks first (in case it's needed)
1969       double dpt = trackRef->ptError();
1970       double blowError = PFTrackAlgoTools::errorScale(trackRef->algo(), factors45_);
1971       // except if it is from an interaction
1972       bool isPrimaryOrSecondary = isFromSecInt(elements[iTrack], "all");
1973
1974       if (isPrimaryOrSecondary)
1975         blowError = 1.;
1976
1977       std::pair<unsigned, bool> tkmuon(iTrack, thisIsALooseMuon);
1978       associatedTracks.emplace(-dpt * blowError, tkmuon);
1979
1980       // Also keep the total track momentum error for comparison with the calo energy
1981       double dp = trackRef->qoverpError() * trackMomentum * trackMomentum;
1982       sumpError2 += dp * dp;
1983
1984       bool connectedToEcal = false;  // Will become true if there is at least one ECAL cluster connected
1985       if (!sortedEcals.empty()) {    // start case: at least one ecal element associated to iTrack
1986
1987         // Loop over all connected ECAL clusters
1988         for (auto const& ecal : sortedEcals) {
1989           unsigned iEcal = ecal.second;
1990           double dist = ecal.first;
1991
1992           // Ignore ECAL cluters already used
1993           if (!active[iEcal]) {
1994             LogTrace("PFAlgo|createCandidatesHCAL") << "cluster locked";
1995             continue;
1996           }
1997
1998           // Sanity checks
1999           PFBlockElement::Type type = elements[iEcal].type();
2000           assert(type == PFBlockElement::ECAL);
2001           PFClusterRef eclusterref = elements[iEcal].clusterRef();
2002           assert(!eclusterref.isNull());
2003
2004           // Check if this ECAL is not closer to another track - ignore it in that case
2005           std::multimap<double, unsigned> sortedTracksEcal;
2006           block.associatedElements(
2007               iEcal, linkData, sortedTracksEcal, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
2008           unsigned jTrack = sortedTracksEcal.begin()->second;
2009           if (jTrack != iTrack)
2010             continue;
2011
2012           double distEcal = block.dist(jTrack, iEcal, linkData, reco::PFBlock::LINKTEST_ALL);
2013
2014           float ecalEnergyCalibrated = eclusterref->correctedEnergy();  // calibrated based on the egamma hypothesis
2015           float ecalEnergy = eclusterref->energy();
2016           ::math::XYZVector photonDirection(
2017               eclusterref->position().X(), eclusterref->position().Y(), eclusterref->position().Z());
2018           photonDirection = photonDirection.Unit();
2019
2020           if (!connectedToEcal) {  // This is the closest ECAL cluster - will add its energy later
2021
2022             LogTrace("PFAlgo|createCandidatesHCAL") << "closest: " << elements[iEcal];
2023
2024             connectedToEcal = true;
2025             // PJ 1st-April-09 : To be done somewhere !!! (Had to comment it, but it is needed)
2026             // currentChargedHadron.addElementInBlock( blockref, iEcal );
2027
2028             // KH: we don't know if this satellite is due to egamma or hadron shower. use raw energy for PF hadron calibration._ store also calibration constant.
2029             double ecalCalibFactor = (ecalEnergy > 1E-9) ? ecalEnergyCalibrated / ecalEnergy : 1.;
2030             std::tuple<unsigned, ::math::XYZVector, double> satellite(
2031                 iEcal, ecalEnergy * photonDirection, ecalCalibFactor);
2032             ecalSatellites.emplace(-1., satellite);
2033
2034           } else {  // Keep satellite clusters for later
2035
2036             // KH: same as above.
2037             double ecalCalibFactor = (ecalEnergy > 1E-9) ? ecalEnergyCalibrated / ecalEnergy : 1.;
2038             std::tuple<unsigned, ::math::XYZVector, double> satellite(
2039                 iEcal, ecalEnergy * photonDirection, ecalCalibFactor);
2040             ecalSatellites.emplace(dist, satellite);
2041           }
2042
2043           std::pair<double, unsigned> associatedEcal(distEcal, iEcal);
2044           associatedEcals.emplace(iTrack, associatedEcal);
2045
2046         }  // End loop ecal associated to iTrack
2047       }    // end case: at least one ecal element associated to iTrack
2048
2049       if (useHO_ && !sortedHOs.empty()) {  // start case: at least one ho element associated to iTrack
2050
2051         // Loop over all connected HO clusters
2052         for (auto const& ho : sortedHOs) {
2053           unsigned iHO = ho.second;
2054           double distHO = ho.first;
2055
2056           // Ignore HO cluters already used
2057           if (!active[iHO]) {
2058             LogTrace("PFAlgo|createCandidatesHCAL") << "cluster locked";
2059             continue;
2060           }
2061
2062           // Sanity checks
2063           PFBlockElement::Type type = elements[iHO].type();
2064           assert(type == PFBlockElement::HO);
2065           PFClusterRef hoclusterref = elements[iHO].clusterRef();
2066           assert(!hoclusterref.isNull());
2067
2068           // Check if this HO is not closer to another track - ignore it in that case
2069           std::multimap<double, unsigned> sortedTracksHO;
2070           block.associatedElements(
2071               iHO, linkData, sortedTracksHO, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
2072           unsigned jTrack = sortedTracksHO.begin()->second;
2073           if (jTrack != iTrack)
2074             continue;
2075
2076           // double chi2HO = block.chi2(jTrack,iHO,linkData,
2077           //                              reco::PFBlock::LINKTEST_ALL);
2078           //double distHO = block.dist(jTrack,iHO,linkData,
2079           //               reco::PFBlock::LINKTEST_ALL);
2080
2081           // Increment the total energy by the energy of this HO cluster
2082           totalHO += hoclusterref->energy();
2083           active[iHO] = false;
2084           // Keep track for later reference in the PFCandidate.
2085           std::pair<double, unsigned> associatedHO(distHO, iHO);
2086           associatedHOs.emplace(iTrack, associatedHO);
2087
2088         }  // End loop ho associated to iTrack
2089       }    // end case: at least one ho element associated to iTrack
2090
2091     }  // end loop on tracks associated to hcal element iHcal
2092
2093     // Include totalHO in totalHCAL for the time being (it will be calibrated as HCAL energy)
2094     totalHcal += totalHO;
2095
2096     // test compatibility between calo and tracker. //////////////
2097
2098     double caloEnergy = 0.;
2099     double slopeEcal = 1.0;
2100     double calibEcal = 0.;
2101     double calibHcal = 0.;
2102     hadronDirection = hadronAtECAL.Unit();
2103
2104     // Determine the expected calo resolution from the total charged momentum
2105     double caloResolution = neutralHadronEnergyResolution(totalChargedMomentum, hclusterref->positionREP().Eta());
2106     caloResolution *= totalChargedMomentum;
2107     // Account for muons
2108     caloResolution = std::sqrt(caloResolution * caloResolution + muonHCALError + muonECALError);
2109     totalEcal -= std::min(totalEcal, muonECALEnergy);
2110     totalEcalEGMCalib -= std::min(totalEcalEGMCalib, muonECALEnergy);
2111     totalHcal -= std::min(totalHcal, muonHCALEnergy);
2112     if (totalEcal < 1E-9)
2113       photonAtECAL = ::math::XYZVector(0., 0., 0.);
2114     if (totalHcal < 1E-9)
2115       hadronAtECAL = ::math::XYZVector(0., 0., 0.);
2116
2117     // Loop over all ECAL satellites, starting for the closest to the various tracks
2118     // and adding other satellites until saturation of the total track momentum
2119     // Note : for code simplicity, the first element of the loop is the HCAL cluster
2120     // with 0 energy in the ECAL
2121     for (auto const& ecalSatellite : ecalSatellites) {
2122       // Add the energy of this ECAL cluster
2123       double previousCalibEcal = calibEcal;
2124       double previousCalibHcal = calibHcal;
2125       double previousCaloEnergy = caloEnergy;
2126       double previousSlopeEcal = slopeEcal;
2127       ::math::XYZVector previousHadronAtECAL = hadronAtECAL;
2128       //
2129       totalEcal +=
2130           sqrt(std::get<1>(ecalSatellite.second).Mag2());  // KH: raw ECAL energy for input to PF hadron calibration
2131       totalEcalEGMCalib += sqrt(std::get<1>(ecalSatellite.second).Mag2()) *
2132                            std::get<2>(ecalSatellite.second);  // KH: calibrated ECAL energy under the egamma hypothesis
2133       photonAtECAL += std::get<1>(ecalSatellite.second) *
2134                       std::get<2>(ecalSatellite.second);  // KH: calibrated ECAL energy under the egamma hypothesis
2135       calibEcal = std::max(0., totalEcal);                // KH: preparing for hadron calibration
2136       calibHcal = std::max(0., totalHcal);
2137       hadronAtECAL = calibHcal * hadronDirection;
2138       // Calibrate ECAL and HCAL energy under the hadron hypothesis.
2139       calibration_.energyEmHad(totalChargedMomentum,
2140                                calibEcal,
2141                                calibHcal,
2142                                hclusterref->positionREP().Eta(),
2143                                hclusterref->positionREP().Phi());
2144       caloEnergy = calibEcal + calibHcal;
2145       if (totalEcal > 0.)
2146         slopeEcal = calibEcal / totalEcal;
2147
2148       hadronAtECAL = calibHcal * hadronDirection;
2149
2150       // Continue looping until all closest clusters are exhausted and as long as
2151       // the calorimetric energy does not saturate the total momentum.
2152       if (ecalSatellite.first < 0. || caloEnergy - totalChargedMomentum <= 0.) {
2153         LogTrace("PFAlgo|createCandidatesHCAL")
2154             << "\t\t\tactive, adding " << std::get<1>(ecalSatellite.second) << " to ECAL energy, and locking";
2155         active[std::get<0>(ecalSatellite.second)] = false;
2156         double clusterEnergy =
2157             sqrt(std::get<1>(ecalSatellite.second).Mag2()) *
2158             std::get<2>(ecalSatellite.second);  // KH: ECAL energy calibrated under the egamma hypothesis
2159         if (clusterEnergy > 50) {               // KH: used to split energetic ecal clusters (E>50 GeV)
2160           ecalClusters.push_back(ecalSatellite.second);
2161           sumEcalClusters += clusterEnergy;
2162         }
2163         continue;
2164       }
2165
2166       // Otherwise, do not consider the last cluster examined and exit.
2167       // active[is->second.first] = true;
2168       totalEcal -= sqrt(std::get<1>(ecalSatellite.second).Mag2());
2169       totalEcalEGMCalib -= sqrt(std::get<1>(ecalSatellite.second).Mag2()) * std::get<2>(ecalSatellite.second);
2170       photonAtECAL -= std::get<1>(ecalSatellite.second) * std::get<2>(ecalSatellite.second);
2171       calibEcal = previousCalibEcal;
2172       calibHcal = previousCalibHcal;
2173       hadronAtECAL = previousHadronAtECAL;
2174       slopeEcal = previousSlopeEcal;
2175       caloEnergy = previousCaloEnergy;
2176
2177       break;
2178     }
2179
2180     // Sanity check !
2181     assert(caloEnergy >= 0);
2182
2183     // And now check for hadronic energy excess...
2184
2185     //colin: resolution should be measured on the ecal+hcal case.
2186     // however, the result will be close.
2187     // double caloResolution = neutralHadronEnergyResolution( caloEnergy );
2188     // caloResolution *= caloEnergy;
2189     // PJ The resolution is on the expected charged calo energy !
2190     //double caloResolution = neutralHadronEnergyResolution( totalChargedMomentum, hclusterref->positionREP().Eta());
2191     //caloResolution *= totalChargedMomentum;
2192     // that of the charged particles linked to the cluster!
2193
2194     ////////////////////// TRACKER MUCH LARGER THAN CALO /////////////////////////
2195     if (totalChargedMomentum - caloEnergy > nSigmaTRACK_ * caloResolution) {
2196       // First consider loose muons
2197       if (nMuons > 0) {
2198         for (auto const& trk : associatedTracks) {
2199           // Only muons
2200           if (!trk.second.second)
2201             continue;
2202
2203           const unsigned int iTrack = trk.second.first;
2204           // Only active tracks
2205           if (!active[iTrack])
2206             continue;
2207
2208           const double trackMomentum = elements[trk.second.first].trackRef()->p();
2209
2210           // look for ECAL elements associated to iTrack (associated to iHcal)
2211           std::multimap<double, unsigned> sortedEcals;
2212           block.associatedElements(
2213               iTrack, linkData, sortedEcals, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
2214           std::multimap<double, unsigned> sortedHOs;
2215           block.associatedElements(iTrack, linkData, sortedHOs, reco::PFBlockElement::HO, reco::PFBlock::LINKTEST_ALL);
2216
2217           //Here allow for loose muons!
2218           auto& muon = (*pfCandidates_)[reconstructTrack(elements[iTrack], true)];
2219
2220           muon.addElementInBlock(blockref, iTrack);
2221           muon.addElementInBlock(blockref, iHcal);
2222           const double muonHcal = std::min(muonHCAL_[0] + muonHCAL_[1], totalHcal - totalHO);
2223           double muonHO = 0.;
2224           muon.setHcalEnergy(totalHcal, muonHcal);
2225           if (!sortedEcals.empty()) {
2226             const unsigned int iEcal = sortedEcals.begin()->second;
2227             PFClusterRef eclusterref = elements[iEcal].clusterRef();
2228             muon.addElementInBlock(blockref, iEcal);
2229             const double muonEcal = std::min(muonECAL_[0] + muonECAL_[1], eclusterref->energy());
2230             muon.setEcalEnergy(eclusterref->energy(), muonEcal);
2231           }
2232           if (useHO_ && !sortedHOs.empty()) {
2233             const unsigned int iHO = sortedHOs.begin()->second;
2234             PFClusterRef hoclusterref = elements[iHO].clusterRef();
2235             muon.addElementInBlock(blockref, iHO);
2236             muonHO = std::min(muonHO_[0] + muonHO_[1], hoclusterref->energy());
2237             muon.setHcalEnergy(max(totalHcal - totalHO, 0.0), muonHcal);
2238             muon.setHoEnergy(hoclusterref->energy(), muonHO);
2239           }
2240           setHcalDepthInfo(muon, *hclusterref);
2241           // Remove it from the block
2242           const ::math::XYZPointF& chargedPosition =
2243               dynamic_cast<const reco::PFBlockElementTrack*>(&elements[trk.second.first])->positionAtECALEntrance();
2244           ::math::XYZVector chargedDirection(chargedPosition.X(), chargedPosition.Y(), chargedPosition.Z());
2245           chargedDirection = chargedDirection.Unit();
2246           totalChargedMomentum -= trackMomentum;
2247           // Update the calo energies
2248           if (totalEcal > 0.)
2249             calibEcal -= std::min(calibEcal, muonECAL_[0] * calibEcal / totalEcal);
2250           if (totalHcal > 0.)
2251             calibHcal -= std::min(calibHcal, muonHCAL_[0] * calibHcal / totalHcal);
2252           totalEcal -= std::min(totalEcal, muonECAL_[0]);
2253           totalHcal -= std::min(totalHcal, muonHCAL_[0]);
2254           if (totalEcal > muonECAL_[0])
2255             photonAtECAL -= muonECAL_[0] * chargedDirection;
2256           if (totalHcal > muonHCAL_[0])
2257             hadronAtECAL -= muonHCAL_[0] * calibHcal / totalHcal * chargedDirection;
2258           caloEnergy = calibEcal + calibHcal;
2259           muonHCALEnergy += muonHCAL_[0];
2260           muonHCALError += muonHCAL_[1] * muonHCAL_[1];
2261           if (muonHO > 0.) {
2262             muonHCALEnergy += muonHO_[0];
2263             muonHCALError += muonHO_[1] * muonHO_[1];
2264             if (totalHcal > 0.) {
2265               calibHcal -= std::min(calibHcal, muonHO_[0] * calibHcal / totalHcal);
2266               totalHcal -= std::min(totalHcal, muonHO_[0]);
2267             }
2268           }
2269           muonECALEnergy += muonECAL_[0];
2270           muonECALError += muonECAL_[1] * muonECAL_[1];
2271           active[iTrack] = false;
2272           // Stop the loop whenever enough muons are removed
2273           //Commented out: Keep looking for muons since they often come in pairs -Matt
2274           //if ( totalChargedMomentum < caloEnergy ) break;
2275         }
2276         // New calo resolution.
2277         caloResolution = neutralHadronEnergyResolution(totalChargedMomentum, hclusterref->positionREP().Eta());
2278         caloResolution *= totalChargedMomentum;
2279         caloResolution = std::sqrt(caloResolution * caloResolution + muonHCALError + muonECALError);
2280       }
2281     }
2282
2283 #ifdef EDM_ML_DEBUG
2284     LogTrace("PFAlgo|createCandidatesHCAL") << "\tBefore Cleaning ";
2285     LogTrace("PFAlgo|createCandidatesHCAL") << "\tCompare Calo Energy to total charged momentum ";
2286     LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tsum p    = " << totalChargedMomentum << " +- " << sqrt(sumpError2);
2287     LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tsum ecal = " << totalEcal;
2288     LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tsum hcal = " << totalHcal;
2289     LogTrace("PFAlgo|createCandidatesHCAL") << "\t\t => Calo Energy = " << caloEnergy << " +- " << caloResolution;
2290     LogTrace("PFAlgo|createCandidatesHCAL")
2291         << "\t\t => Calo Energy- total charged momentum = " << caloEnergy - totalChargedMomentum << " +- "
2292         << sqrt(sumpError2 + caloResolution * caloResolution);
2293 #endif
2294
2295     // Second consider bad tracks (if still needed after muon removal)
2296     unsigned corrTrack = 10000000;
2297     double corrFact = 1.;
2298
2299     if (rejectTracks_Bad_ && totalChargedMomentum - caloEnergy > nSigmaTRACK_ * caloResolution) {
2300       for (auto const& trk : associatedTracks) {
2301         const unsigned iTrack = trk.second.first;
2302         // Only active tracks
2303         if (!active[iTrack])
2304           continue;
2305         const reco::TrackRef& trackref = elements[trk.second.first].trackRef();
2306
2307         const double dptRel = fabs(trk.first) / trackref->pt() * 100;
2308         const bool isSecondary = isFromSecInt(elements[iTrack], "secondary");
2309         const bool isPrimary = isFromSecInt(elements[iTrack], "primary");
2310
2311         if (isSecondary && dptRel < dptRel_DispVtx_)
2312           continue;
2313         // Consider only bad tracks
2314         if (fabs(trk.first) < ptError_)
2315           break;
2316         // What would become the block charged momentum if this track were removed
2317         const double wouldBeTotalChargedMomentum = totalChargedMomentum - trackref->p();
2318         // Reject worst tracks, as long as the total charged momentum
2319         // is larger than the calo energy
2320
2321         if (wouldBeTotalChargedMomentum > caloEnergy) {
2322           if (isSecondary) {
2323             LogTrace("PFAlgo|createCandidatesHCAL")
2324                 << "In bad track rejection step dptRel = " << dptRel << " dptRel_DispVtx_ = " << dptRel_DispVtx_;
2325             LogTrace("PFAlgo|createCandidatesHCAL")
2326                 << "The calo energy would be still smaller even without this track but it is attached to a NI";
2327           }
2328
2329           if (isPrimary || (isSecondary && dptRel < dptRel_DispVtx_))
2330             continue;
2331           active[iTrack] = false;
2332           totalChargedMomentum = wouldBeTotalChargedMomentum;
2333           LogTrace("PFAlgo|createCandidatesHCAL")
2334               << "\tElement  " << elements[iTrack] << " rejected (dpt = " << -trk.first
2335               << " GeV/c, algo = " << trackref->algo() << ")";
2336           // Just rescale the nth worst track momentum to equalize the calo energy
2337         } else {
2338           if (isPrimary)
2339             break;
2340           corrTrack = iTrack;
2341           corrFact = (caloEnergy - wouldBeTotalChargedMomentum) / elements[trk.second.first].trackRef()->p();
2342           if (trackref->p() * corrFact < 0.05) {
2343             corrFact = 0.;
2344             active[iTrack] = false;
2345           }
2346           totalChargedMomentum -= trackref->p() * (1. - corrFact);
2347           LogTrace("PFAlgo|createCandidatesHCAL")
2348               << "\tElement  " << elements[iTrack] << " (dpt = " << -trk.first << " GeV/c, algo = " << trackref->algo()
2349               << ") rescaled by " << corrFact << " Now the total charged momentum is " << totalChargedMomentum;
2350           break;
2351         }
2352       }
2353     }
2354
2355     // New determination of the calo and track resolution avec track deletion/rescaling.
2356     caloResolution = neutralHadronEnergyResolution(totalChargedMomentum, hclusterref->positionREP().Eta());
2357     caloResolution *= totalChargedMomentum;
2358     caloResolution = std::sqrt(caloResolution * caloResolution + muonHCALError + muonECALError);
2359
2360     // Check if the charged momentum is still very inconsistent with the calo measurement.
2361     // In this case, just drop all tracks from 4th and 5th iteration linked to this block
2362
2363     if (rejectTracks_Step45_ && sortedTracks.size() > 1 &&
2364         totalChargedMomentum - caloEnergy > nSigmaTRACK_ * caloResolution) {
2365       for (auto const& trk : associatedTracks) {
2366         unsigned iTrack = trk.second.first;
2367         reco::TrackRef trackref = elements[iTrack].trackRef();
2368         if (!active[iTrack])
2369           continue;
2370
2371         double dptRel = fabs(trk.first) / trackref->pt() * 100;
2372         bool isPrimaryOrSecondary = isFromSecInt(elements[iTrack], "all");
2373
2374         if (isPrimaryOrSecondary && dptRel < dptRel_DispVtx_)
2375           continue;
2376
2377         if (PFTrackAlgoTools::step5(trackref->algo())) {
2378           active[iTrack] = false;
2379           totalChargedMomentum -= trackref->p();
2380
2381           LogTrace("PFAlgo|createCandidatesHCAL")
2382               << "\tElement  " << elements[iTrack] << " rejected (dpt = " << -trk.first
2383               << " GeV/c, algo = " << trackref->algo() << ")";
2384         }
2385       }
2386     }
2387
2388     // New determination of the calo and track resolution avec track deletion/rescaling.
2389     caloResolution = neutralHadronEnergyResolution(totalChargedMomentum, hclusterref->positionREP().Eta());
2390     caloResolution *= totalChargedMomentum;
2391     caloResolution = std::sqrt(caloResolution * caloResolution + muonHCALError + muonECALError);
2392
2393     // Make PF candidates with the remaining tracks in the block
2394     sumpError2 = 0.;
2395     for (auto const& trk : associatedTracks) {
2396       unsigned iTrack = trk.second.first;
2397       if (!active[iTrack])
2398         continue;
2399       reco::TrackRef trackRef = elements[iTrack].trackRef();
2400       double trackMomentum = trackRef->p();
2401       double dp = trackRef->qoverpError() * trackMomentum * trackMomentum;
2402       unsigned tmpi = reconstructTrack(elements[iTrack]);
2403
2404       (*pfCandidates_)[tmpi].addElementInBlock(blockref, iTrack);
2405       (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHcal);
2406       setHcalDepthInfo((*pfCandidates_)[tmpi], *hclusterref);
2407       auto myEcals = associatedEcals.equal_range(iTrack);
2408       for (auto ii = myEcals.first; ii != myEcals.second; ++ii) {
2409         unsigned iEcal = ii->second.second;
2410         if (active[iEcal])
2411           continue;
2412         (*pfCandidates_)[tmpi].addElementInBlock(blockref, iEcal);
2413       }
2414
2415       if (useHO_) {
2416         auto myHOs = associatedHOs.equal_range(iTrack);
2417         for (auto ii = myHOs.first; ii != myHOs.second; ++ii) {
2418           unsigned iHO = ii->second.second;
2419           if (active[iHO])
2420             continue;
2421           (*pfCandidates_)[tmpi].addElementInBlock(blockref, iHO);
2422         }
2423       }
2424
2425       if (iTrack == corrTrack) {
2426         if (corrFact < 0.)
2427           corrFact = 0.;  // protect against negative scaling
2428         (*pfCandidates_)[tmpi].rescaleMomentum(corrFact);
2429         trackMomentum *= corrFact;
2430       }
2431       chargedHadronsIndices.push_back(tmpi);
2432       chargedHadronsInBlock.push_back(iTrack);
2433       active[iTrack] = false;
2434       hcalP.push_back(trackMomentum);
2435       hcalDP.push_back(dp);
2436       if (dp / trackMomentum > maxDPovP)
2437         maxDPovP = dp / trackMomentum;
2438       sumpError2 += dp * dp;
2439     }
2440
2441     // The total uncertainty of the difference Calo-Track
2442     double totalError = sqrt(sumpError2 + caloResolution * caloResolution);
2443
2444 #ifdef EDM_ML_DEBUG
2445     LogTrace("PFAlgo|createCandidatesHCAL")
2446         << "\tCompare Calo Energy to total charged momentum " << endl
2447         << "\t\tsum p    = " << totalChargedMomentum << " +- " << sqrt(sumpError2) << endl
2448         << "\t\tsum ecal = " << totalEcal << endl
2449         << "\t\tsum hcal = " << totalHcal << endl
2450         << "\t\t => Calo Energy = " << caloEnergy << " +- " << caloResolution << endl
2451         << "\t\t => Calo Energy- total charged momentum = " << caloEnergy - totalChargedMomentum << " +- "
2452         << totalError;
2453 #endif
2454
2455     /* */
2456
2457     /////////////// TRACKER AND CALO COMPATIBLE  ////////////////
2458     double nsigma = nSigmaHCAL(totalChargedMomentum, hclusterref->positionREP().Eta());
2459     //double nsigma = nSigmaHCAL(caloEnergy,hclusterref->positionREP().Eta());
2460     if (abs(totalChargedMomentum - caloEnergy) < nsigma * totalError) {
2461       // deposited caloEnergy compatible with total charged momentum
2462       // if tracking errors are large take weighted average
2463
2464 #ifdef EDM_ML_DEBUG
2465       LogTrace("PFAlgo|createCandidatesHCAL")
2466           << "\t\tcase 1: COMPATIBLE "
2467           << "|Calo Energy- total charged momentum| = " << abs(caloEnergy - totalChargedMomentum) << " < " << nsigma
2468           << " x " << totalError;
2469       if (maxDPovP < 0.1)
2470         LogTrace("PFAlgo|createCandidatesHCAL") << "\t\t\tmax DP/P = " << maxDPovP << " less than 0.1: do nothing ";
2471       else
2472         LogTrace("PFAlgo|createCandidatesHCAL")
2473             << "\t\t\tmax DP/P = " << maxDPovP << " >  0.1: take weighted averages ";
2474 #endif
2475
2476       // if max DP/P < 10%  do nothing
2477       if (maxDPovP > 0.1) {
2478         // for each track associated to hcal
2479         //      int nrows = tkIs.size();
2480         int nrows = chargedHadronsIndices.size();
2481         TMatrixTSym<double> a(nrows);
2482         TVectorD b(nrows);
2483         TVectorD check(nrows);
2484         double sigma2E = caloResolution * caloResolution;
2485         for (int i = 0; i < nrows; i++) {
2486           double sigma2i = hcalDP[i] * hcalDP[i];
2487           LogTrace("PFAlgo|createCandidatesHCAL")
2488               << "\t\t\ttrack associated to hcal " << i << " P = " << hcalP[i] << " +- " << hcalDP[i];
2489           a(i, i) = 1. / sigma2i + 1. / sigma2E;
2490           b(i) = hcalP[i] / sigma2i + caloEnergy / sigma2E;
2491           for (int j = 0; j < nrows; j++) {
2492             if (i == j)
2493               continue;
2494             a(i, j) = 1. / sigma2E;
2495           }  // end loop on j
2496         }    // end loop on i
2497
2498         // solve ax = b
2499         TDecompChol decomp(a);
2500         bool ok = false;
2501         TVectorD x = decomp.Solve(b, ok);
2502         // for each track create a PFCandidate track
2503         // with a momentum rescaled to weighted average
2504         if (ok) {
2505           for (int i = 0; i < nrows; i++) {
2506             //      unsigned iTrack = trackInfos[i].index;
2507             unsigned ich = chargedHadronsIndices[i];
2508             double rescaleFactor = x(i) / hcalP[i];
2509             if (rescaleFactor < 0.)
2510               rescaleFactor = 0.;  // protect against negative scaling
2511             (*pfCandidates_)[ich].rescaleMomentum(rescaleFactor);
2512
2513             LogTrace("PFAlgo|createCandidatesHCAL")
2514                 << "\t\t\told p " << hcalP[i] << " new p " << x(i) << " rescale " << rescaleFactor;
2515           }
2516         } else {
2517           edm::LogError("PFAlgo::createCandidatesHCAL") << "TDecompChol.Solve returned ok=false";
2518           assert(0);
2519         }
2520       }
2521     }
2522
2523     /////////////// NEUTRAL DETECTION  ////////////////
2524     else if (caloEnergy > totalChargedMomentum) {
2525       //case 2: caloEnergy > totalChargedMomentum + nsigma*totalError
2526       //there is an excess of energy in the calos
2527       //create a neutral hadron or a photon
2528
2529       double eNeutralHadron = caloEnergy - totalChargedMomentum;
2530       double ePhoton = (caloEnergy - totalChargedMomentum) /
2531                        slopeEcal;  // KH: this slopeEcal is computed based on ECAL energy under the hadron hypothesis,
2532                                    // thought we are creating photons.
2533       // This is a fuzzy case, but it should be better than corrected twice under both egamma and hadron hypotheses.
2534
2535 #ifdef EDM_ML_DEBUG
2536       if (!sortedTracks.empty()) {
2537         LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tcase 2: NEUTRAL DETECTION " << caloEnergy << " > " << nsigma
2538                                                 << "x" << totalError << " + " << totalChargedMomentum;
2539         LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tneutral activity detected: " << endl
2540                                                 << "\t\t\t           photon = " << ePhoton << endl
2541                                                 << "\t\t\tor neutral hadron = " << eNeutralHadron;
2542
2543         LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tphoton or hadron ?";
2544       }
2545
2546       if (sortedTracks.empty())
2547         LogTrace("PFAlgo|createCandidatesHCAL") << "\t\tno track -> hadron ";
2548       else
2549         LogTrace("PFAlgo|createCandidatesHCAL")
2550             << "\t\t" << sortedTracks.size() << " tracks -> check if the excess is photonic or hadronic";
2551 #endif
2552
2553       double ratioMax = 0.;
2554       reco::PFClusterRef maxEcalRef;
2555       unsigned maxiEcal = 9999;
2556
2557       // for each track associated to hcal: iterator IE ie :
2558
2559       LogTrace("PFAlgo|createCandidatesHCAL") << "loop over sortedTracks.size()=" << sortedTracks.size();
2560       for (auto const& trk : sortedTracks) {
2561         unsigned iTrack = trk.second;
2562
2563         PFBlockElement::Type type = elements[iTrack].type();
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         // double distECAL = iae->second.first;
2575         unsigned iEcal = iae->second.second;
2576
2577         PFBlockElement::Type typeEcal = elements[iEcal].type();
2578         assert(typeEcal == reco::PFBlockElement::ECAL);
2579
2580         reco::PFClusterRef clusterRef = elements[iEcal].clusterRef();
2581         assert(!clusterRef.isNull());
2582
2583         double pTrack = trackRef->p();
2584         double eECAL = clusterRef->energy();
2585         double eECALOverpTrack = eECAL / pTrack;
2586
2587         if (eECALOverpTrack > ratioMax) {
2588           ratioMax = eECALOverpTrack;
2589           maxEcalRef = clusterRef;
2590           maxiEcal = iEcal;
2591         }
2592
2593       }  // end loop on tracks associated to hcal: iterator IE ie
2594
2595       std::vector<reco::PFClusterRef> pivotalClusterRef;
2596       std::vector<unsigned> iPivotal;
2597       std::vector<double> particleEnergy, ecalEnergy, hcalEnergy, rawecalEnergy, rawhcalEnergy;
2598       std::vector<::math::XYZVector> particleDirection;
2599
2600       // If the excess is smaller than the ecal energy, assign the whole
2601       // excess to photons
2602       if (ePhoton < totalEcal || eNeutralHadron - calibEcal < 1E-10) {
2603         if (!maxEcalRef.isNull()) {
2604           // So the merged photon energy is,
2605           mergedPhotonEnergy = ePhoton;
2606         }
2607       } else {
2608         // Otherwise assign the whole ECAL energy to the photons
2609         if (!maxEcalRef.isNull()) {
2610           // So the merged photon energy is,
2611           mergedPhotonEnergy = totalEcalEGMCalib;  // KH: use calibrated ECAL energy under the egamma hypothesis
2612         }
2613         // ... and assign the remaining excess to neutral hadrons using the direction of ecal clusters
2614         mergedNeutralHadronEnergy = eNeutralHadron - calibEcal;
2615       }
2616
2617       if (mergedPhotonEnergy > 0) {
2618         // Split merged photon into photons for each energetic ecal cluster (necessary for jet substructure reconstruction)
2619         // make only one merged photon if less than 2 ecal clusters
2620         // KH: this part still needs review, after using non-corrected ECAL energy for PF hadron calibrations
2621         if (ecalClusters.size() <= 1) {
2622           ecalClusters.clear();
2623           ecalClusters.emplace_back(
2624               maxiEcal,
2625               photonAtECAL,
2626               1.);  // KH: calibration factor of 1, which should be ok as long as sumEcalClusters is consistent with photonAtECAL in this case
2627           sumEcalClusters = sqrt(photonAtECAL.Mag2());
2628         }
2629         for (auto const& pae : ecalClusters) {
2630           const double clusterEnergyCalibrated =
2631               sqrt(std::get<1>(pae).Mag2()) *
2632               std::get<2>(
2633                   pae);  // KH: calibrated under the egamma hypothesis. Note: sumEcalClusters is normally calibrated under egamma hypothesis
2634           particleEnergy.push_back(mergedPhotonEnergy * clusterEnergyCalibrated / sumEcalClusters);
2635           particleDirection.push_back(std::get<1>(pae));
2636           ecalEnergy.push_back(mergedPhotonEnergy * clusterEnergyCalibrated / sumEcalClusters);
2637           hcalEnergy.push_back(0.);
2638           rawecalEnergy.push_back(totalEcal);
2639           rawhcalEnergy.push_back(0.);
2640           pivotalClusterRef.push_back(elements[std::get<0>(pae)].clusterRef());
2641           iPivotal.push_back(std::get<0>(pae));
2642         }
2643       }  // mergedPhotonEnergy > 0
2644
2645       if (mergedNeutralHadronEnergy > 1.0) {
2646         // Split merged neutral hadrons according to directions of energetic ecal clusters (necessary for jet substructure reconstruction)
2647         // make only one merged neutral hadron if less than 2 ecal clusters
2648         if (ecalClusters.size() <= 1) {
2649           ecalClusters.clear();
2650           ecalClusters.emplace_back(
2651               iHcal,
2652               hadronAtECAL,
2653               1.);  // KH: calibration factor of 1, which should be ok as long as sumEcalClusters is consistent with photonAtECAL
2654           sumEcalClusters = sqrt(hadronAtECAL.Mag2());
2655         }
2656         for (auto const& pae : ecalClusters) {
2657           const double clusterEnergyCalibrated =
2658               sqrt(std::get<1>(pae).Mag2()) *
2659               std::get<2>(
2660                   pae);  // KH: calibrated under the egamma hypothesis. Note: sumEcalClusters is normally calibrated under egamma hypothesis
2661           particleEnergy.push_back(mergedNeutralHadronEnergy * clusterEnergyCalibrated / sumEcalClusters);
2662           particleDirection.push_back(std::get<1>(pae));
2663           ecalEnergy.push_back(0.);
2664           hcalEnergy.push_back(mergedNeutralHadronEnergy * clusterEnergyCalibrated / sumEcalClusters);
2665           rawecalEnergy.push_back(0.);
2666           rawhcalEnergy.push_back(totalHcal);
2667           pivotalClusterRef.push_back(hclusterref);
2668           iPivotal.push_back(iHcal);
2669         }
2670       }  //mergedNeutralHadronEnergy > 1.0
2671
2672       // reconstructing a merged neutral
2673       // the type of PFCandidate is known from the
2674       // reference to the pivotal cluster.
2675
2676       for (unsigned iPivot = 0; iPivot < iPivotal.size(); ++iPivot) {
2677         if (particleEnergy[iPivot] < 0.)
2678           edm::LogWarning("PFAlgo|createCandidatesHCAL")
2679               << "ALARM = Negative energy for iPivot=" << iPivot << ", " << particleEnergy[iPivot];
2680
2681         const bool useDirection = true;
2682         auto& neutral = (*pfCandidates_)[reconstructCluster(*pivotalClusterRef[iPivot],
2683                                                             particleEnergy[iPivot],
2684                                                             useDirection,
2685                                                             particleDirection[iPivot].X(),
2686                                                             particleDirection[iPivot].Y(),
2687                                                             particleDirection[iPivot].Z())];
2688
2689         neutral.setEcalEnergy(rawecalEnergy[iPivot], ecalEnergy[iPivot]);
2690         if (!useHO_) {
2691           neutral.setHcalEnergy(rawhcalEnergy[iPivot], hcalEnergy[iPivot]);
2692           neutral.setHoEnergy(0., 0.);
2693         } else {                              // useHO_
2694           if (rawhcalEnergy[iPivot] == 0.) {  // photons should be here
2695             neutral.setHcalEnergy(0., 0.);
2696             neutral.setHoEnergy(0., 0.);
2697           } else {
2698             neutral.setHcalEnergy(max(rawhcalEnergy[iPivot] - totalHO, 0.0),
2699                                   hcalEnergy[iPivot] * max(1. - totalHO / rawhcalEnergy[iPivot], 0.));
2700             neutral.setHoEnergy(totalHO, totalHO * hcalEnergy[iPivot] / rawhcalEnergy[iPivot]);
2701           }
2702         }
2703         neutral.setPs1Energy(0.);
2704         neutral.setPs2Energy(0.);
2705         neutral.set_mva_nothing_gamma(-1.);
2706         // neutral.addElement(&elements[iPivotal]);
2707         // neutral.addElementInBlock(blockref, iPivotal[iPivot]);
2708         neutral.addElementInBlock(blockref, iHcal);
2709         for (unsigned iTrack : chargedHadronsInBlock) {
2710           neutral.addElementInBlock(blockref, iTrack);
2711           // Assign the position of the track at the ECAL entrance
2712           const ::math::XYZPointF& chargedPosition =
2713               dynamic_cast<const reco::PFBlockElementTrack*>(&elements[iTrack])->positionAtECALEntrance();
2714           neutral.setPositionAtECALEntrance(chargedPosition);
2715
2716           auto myEcals = associatedEcals.equal_range(iTrack);
2717           for (auto ii = myEcals.first; ii != myEcals.second; ++ii) {
2718             unsigned iEcal = ii->second.second;
2719             if (active[iEcal])
2720               continue;
2721             neutral.addElementInBlock(blockref, iEcal);
2722           }
2723         }
2724       }
2725
2726     }  // excess of energy
2727
2728     // will now share the hcal energy between the various charged hadron
2729     // candidates, taking into account the potential neutral hadrons
2730
2731     //JB: The question is: we've resolved the merged photons cleanly, but how should
2732     //the remaining hadrons be assigned the remaining ecal energy?
2733     //*Temporary solution*: follow HCAL example with fractions...
2734
2735     // remove the energy of the potential neutral hadron
2736     double totalHcalEnergyCalibrated = std::max(calibHcal - mergedNeutralHadronEnergy, 0.);
2737     // similarly for the merged photons
2738     double totalEcalEnergyCalibrated = std::max(calibEcal - mergedPhotonEnergy, 0.);
2739     // share between the charged hadrons
2740
2741     //COLIN can compute this before
2742     // not exactly equal to sum p, this is sum E
2743     double chargedHadronsTotalEnergy = 0;
2744     for (unsigned index : chargedHadronsIndices) {
2745       reco::PFCandidate& chargedHadron = (*pfCandidates_)[index];
2746       chargedHadronsTotalEnergy += chargedHadron.energy();
2747     }
2748
2749     for (unsigned index : chargedHadronsIndices) {
2750       reco::PFCandidate& chargedHadron = (*pfCandidates_)[index];
2751       float fraction = chargedHadron.energy() / chargedHadronsTotalEnergy;
2752
2753       if (!useHO_) {
2754         chargedHadron.setHcalEnergy(fraction * totalHcal, fraction * totalHcalEnergyCalibrated);
2755         chargedHadron.setHoEnergy(0., 0.);
2756       } else {
2757         chargedHadron.setHcalEnergy(fraction * max(totalHcal - totalHO, 0.0),
2758                                     fraction * totalHcalEnergyCalibrated * (1. - totalHO / totalHcal));
2759         chargedHadron.setHoEnergy(fraction * totalHO, fraction * totalHO * totalHcalEnergyCalibrated / totalHcal);
2760       }
2761       //JB: fixing up (previously omitted) setting of ECAL energy gouzevit
2762       chargedHadron.setEcalEnergy(fraction * totalEcal, fraction * totalEcalEnergyCalibrated);
2763     }
2764
2765     // Finally treat unused ecal satellites as photons.
2766     for (auto const& ecalSatellite : ecalSatellites) {
2767       // Ignore satellites already taken
2768       unsigned iEcal = std::get<0>(ecalSatellite.second);
2769       if (!active[iEcal])
2770         continue;
2771
2772       // Sanity checks again (well not useful, this time!)
2773       PFBlockElement::Type type = elements[iEcal].type();
2774       assert(type == PFBlockElement::ECAL);
2775       PFClusterRef eclusterref = elements[iEcal].clusterRef();
2776       assert(!eclusterref.isNull());
2777
2778       // Lock the cluster
2779       active[iEcal] = false;
2780
2781       // Find the associated tracks
2782       std::multimap<double, unsigned> assTracks;
2783       block.associatedElements(iEcal, linkData, assTracks, reco::PFBlockElement::TRACK, reco::PFBlock::LINKTEST_ALL);
2784
2785       // Create a photon
2786       double ecalClusterEnergyCalibrated =
2787           sqrt(std::get<1>(ecalSatellite.second).Mag2()) *
2788           std::get<2>(
2789               ecalSatellite.second);  // KH: calibrated under the egamma hypothesis (rawEcalClusterEnergy * calibration)
2790       auto& cand = (*pfCandidates_)[reconstructCluster(*eclusterref, ecalClusterEnergyCalibrated)];
2791       cand.setEcalEnergy(eclusterref->energy(), ecalClusterEnergyCalibrated);
2792       cand.setHcalEnergy(0., 0.);
2793       cand.setHoEnergy(0., 0.);
2794       cand.setPs1Energy(associatedPSs[iEcal].first);
2795       cand.setPs2Energy(associatedPSs[iEcal].second);
2796       cand.addElementInBlock(blockref, iEcal);
2797       cand.addElementInBlock(blockref, sortedTracks.begin()->second);
2798
2799       if (fabs(eclusterref->energy() - sqrt(std::get<1>(ecalSatellite.second).Mag2())) > 1e-3 ||
2800           fabs(eclusterref->correctedEnergy() - ecalClusterEnergyCalibrated) > 1e-3)
2801         edm::LogWarning("PFAlgo:processBlock")
2802             << "ecalCluster vs ecalSatellites look inconsistent (eCluster E, calibE, ecalSatellite E, calib E): "
2803             << eclusterref->energy() << " " << eclusterref->correctedEnergy() << " "
2804             << sqrt(std::get<1>(ecalSatellite.second).Mag2()) << " " << ecalClusterEnergyCalibrated;
2805
2806     }  // ecalSatellites
2807
2808   }  // hcalIs
2809   // end loop on hcal element iHcal= hcalIs[i]
2810   LogTrace("PFAlgo|createCandidatesHCAL") << "end of function PFAlgo::createCandidatesHCAL";
2811 }
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   // Processing the remaining HCAL clusters
2821   LogTrace("PFAlgo|createCandidatesHCALUnlinked")
2822       << "start of function PFAlgo::createCandidatesHCALUnlinked, hcalIs.size()=" << inds.hcalIs.size();
2823
2824   // --------------- loop remaining hcal ------------------
2825
2826   for (unsigned iHcal : inds.hcalIs) {
2827     // Keep ECAL and HO elements for reference in the PFCandidate
2828     std::vector<unsigned> ecalRefs;
2829     std::vector<unsigned> hoRefs;
2830
2831     LogTrace("PFAlgo|createCandidatesHCALUnlinked") << elements[iHcal] << " ";
2832
2833     if (!active[iHcal]) {
2834       LogTrace("PFAlgo|createCandidatesHCALUnlinked") << "not active " << iHcal;
2835       continue;
2836     }
2837
2838     // Find the ECAL elements linked to it
2839     std::multimap<double, unsigned> ecalElems;
2840     block.associatedElements(iHcal, linkData, ecalElems, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
2841
2842     // Loop on these ECAL elements
2843     float totalEcal = 0.;
2844     float ecalMax = 0.;
2845     reco::PFClusterRef eClusterRef;
2846     for (auto const& ecal : ecalElems) {
2847       unsigned iEcal = ecal.second;
2848       double dist = ecal.first;
2849       PFBlockElement::Type type = elements[iEcal].type();
2850       assert(type == PFBlockElement::ECAL);
2851
2852       // Check if already used
2853       if (!active[iEcal])
2854         continue;
2855
2856       // Check the distance (one HCALPlusECAL tower, roughly)
2857       // if ( dist > 0.15 ) continue;
2858
2859       //COLINFEB16
2860       // what could be done is to
2861       // - link by rechit.
2862       // - take in the neutral hadron all the ECAL clusters
2863       // which are within the same CaloTower, according to the distance,
2864       // except the ones which are closer to another HCAL cluster.
2865       // - all the other ECAL linked to this HCAL are photons.
2866       //
2867       // about the closest HCAL cluster.
2868       // it could maybe be easier to loop on the ECAL clusters first
2869       // to cut the links to all HCAL clusters except the closest, as is
2870       // done in the first track loop. But maybe not!
2871       // or add an helper function to the PFAlgo class to ask
2872       // if a given element is the closest of a given type to another one?
2873
2874       // Check if not closer from another free HCAL
2875       std::multimap<double, unsigned> hcalElems;
2876       block.associatedElements(iEcal, linkData, hcalElems, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
2877
2878       const bool isClosest = std::none_of(hcalElems.begin(), hcalElems.end(), [&](auto const& hcal) {
2879         return active[hcal.second] && hcal.first < dist;
2880       });
2881
2882       if (!isClosest)
2883         continue;
2884
2885 #ifdef EDM_ML_DEBUG
2886       LogTrace("PFAlgo|createCandidatesHCALUnlinked")
2887           << "\telement " << elements[iEcal] << " linked with dist " << dist;
2888       LogTrace("PFAlgo|createCandidatesHCALUnlinked") << "Added to HCAL cluster to form a neutral hadron";
2889 #endif
2890
2891       reco::PFClusterRef eclusterRef = elements[iEcal].clusterRef();
2892       assert(!eclusterRef.isNull());
2893
2894       // KH: use raw ECAL energy for PF hadron calibration_.
2895       double ecalEnergy = eclusterRef->energy();  // ecalEnergy = eclusterRef->correctedEnergy();
2896
2897       totalEcal += ecalEnergy;
2898       if (ecalEnergy > ecalMax) {
2899         ecalMax = ecalEnergy;
2900         eClusterRef = eclusterRef;
2901       }
2902
2903       ecalRefs.push_back(iEcal);
2904       active[iEcal] = false;
2905
2906     }  // End loop ECAL
2907
2908     // Now find the HO clusters linked to the HCAL cluster
2909     double totalHO = 0.;
2910     double hoMax = 0.;
2911     //unsigned jHO = 0;
2912     if (useHO_) {
2913       std::multimap<double, unsigned> hoElems;
2914       block.associatedElements(iHcal, linkData, hoElems, reco::PFBlockElement::HO, reco::PFBlock::LINKTEST_ALL);
2915
2916       // Loop on these HO elements
2917       //      double totalHO = 0.;
2918       //      double hoMax = 0.;
2919       //      unsigned jHO = 0;
2920       reco::PFClusterRef hoClusterRef;
2921       for (auto const& ho : hoElems) {
2922         unsigned iHO = ho.second;
2923         double dist = ho.first;
2924         PFBlockElement::Type type = elements[iHO].type();
2925         assert(type == PFBlockElement::HO);
2926
2927         // Check if already used
2928         if (!active[iHO])
2929           continue;
2930
2931         // Check the distance (one HCALPlusHO tower, roughly)
2932         // if ( dist > 0.15 ) continue;
2933
2934         // Check if not closer from another free HCAL
2935         std::multimap<double, unsigned> hcalElems;
2936         block.associatedElements(iHO, linkData, hcalElems, reco::PFBlockElement::HCAL, reco::PFBlock::LINKTEST_ALL);
2937
2938         const bool isClosest = std::none_of(hcalElems.begin(), hcalElems.end(), [&](auto const& hcal) {
2939           return active[hcal.second] && hcal.first < dist;
2940         });
2941
2942         if (!isClosest)
2943           continue;
2944
2945 #ifdef EDM_ML_DEBUG
2946         if (useHO_) {
2947           LogTrace("PFAlgo|createCandidatesHCALUnlinked")
2948               << "\telement " << elements[iHO] << " linked with dist " << dist;
2949           LogTrace("PFAlgo|createCandidatesHCALUnlinked") << "Added to HCAL cluster to form a neutral hadron";
2950         }
2951 #endif
2952
2953         reco::PFClusterRef hoclusterRef = elements[iHO].clusterRef();
2954         assert(!hoclusterRef.isNull());
2955
2956         double hoEnergy =
2957             hoclusterRef->energy();  // calibration_.energyEm(*hoclusterRef,ps1Ene,ps2Ene,crackCorrection);
2958
2959         totalHO += hoEnergy;
2960         if (hoEnergy > hoMax) {
2961           hoMax = hoEnergy;
2962           hoClusterRef = hoclusterRef;
2963           //jHO = iHO;
2964         }
2965
2966         hoRefs.push_back(iHO);
2967         active[iHO] = false;
2968
2969       }  // End loop HO
2970     }
2971
2972     PFClusterRef hclusterRef = elements[iHcal].clusterRef();
2973     assert(!hclusterRef.isNull());
2974
2975     // HCAL energy
2976     double totalHcal = hclusterRef->energy();
2977     // Include the HO energy
2978     if (useHO_)
2979       totalHcal += totalHO;
2980
2981     // Calibration
2982     double calibEcal = totalEcal > 0. ? totalEcal : 0.;
2983     double calibHcal = std::max(0., totalHcal);
2984     if (hclusterRef->layer() == PFLayer::HF_HAD || hclusterRef->layer() == PFLayer::HF_EM) {
2985       calibEcal = totalEcal;
2986     } else {
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   }  //loop hcal elements
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   // --------------- loop ecal ------------------
3023
3024   // for each ecal element iEcal = ecalIs[i] in turn:
3025
3026   for (unsigned i = 0; i < inds.ecalIs.size(); i++) {
3027     unsigned iEcal = inds.ecalIs[i];
3028
3029     LogTrace("PFAlgo|createCandidatesECAL") << "elements[" << iEcal << "]=" << elements[iEcal];
3030
3031     if (!active[iEcal]) {
3032       LogTrace("PFAlgo|createCandidatesECAL") << "iEcal=" << iEcal << " not active";
3033       continue;
3034     }
3035
3036     PFBlockElement::Type type = elements[iEcal].type();
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     // float ecalEnergy = calibration_.energyEm( clusterref->energy() );
3046     double particleEnergy = ecalEnergy;
3047
3048     auto& cand = (*pfCandidates_)[reconstructCluster(*clusterref, particleEnergy)];
3049
3050     cand.setEcalEnergy(clusterref->energy(), ecalEnergy);
3051     cand.setHcalEnergy(0., 0.);
3052     cand.setHoEnergy(0., 0.);
3053     cand.setPs1Energy(0.);
3054     cand.setPs2Energy(0.);
3055     cand.addElementInBlock(blockref, iEcal);
3056
3057   }  // end loop on ecal elements iEcal = ecalIs[i]
3058   LogTrace("PFAlgo|createCandidatesECAL") << "end of function PFALgo::createCandidatesECAL";
3059 }
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   // make a copy of the link data, which will be edited.
3073   PFBlock::LinkData linkData = block.linkData();
3074
3075   // keep track of the elements which are still active.
3076   vector<bool> active(elements.size(), true);
3077
3078   // //PFElectrons:
3079   // usePFElectrons_ external configurable parameter to set the usage of pf electron
3080   std::vector<reco::PFCandidate> tempElectronCandidates;
3081   tempElectronCandidates.clear();
3082
3083   // New EGamma Reconstruction 10/10/2013
3084   if (useEGammaFilters_) {
3085     egammaFilters(blockref, active, pfegamma);
3086   }  // end if use EGammaFilters
3087
3088   //Lock extra conversion tracks not used by Photon Algo
3089   if (usePFConversions_) {
3090     conversionAlgo(elements, active);
3091   }
3092
3093   // In the following elementLoop() function, the primary goal is to deal with tracks that are:
3094   // - not associated to an HCAL cluster
3095   // - not identified as an electron.
3096   // Those tracks should be predominantly relatively low energy charged
3097   // hadrons which are not detected in the ECAL.
3098
3099   // The secondary goal is to prepare for the next loops
3100   // - The ecal and hcal elements are sorted in separate vectors in `ElementIndices inds`
3101   // which will be used as a base for the corresponding loops.
3102   // - For tracks which are connected to more than one HCAL cluster,
3103   // the links between the track and the cluster are cut for all clusters
3104   // but the closest one.
3105   // - HF only blocks ( HFEM, HFHAD, HFEM+HFAD) are identified
3106
3107   // obsolete comments?
3108   // loop1:
3109   // - sort ecal and hcal elements in separate vectors
3110   // - for tracks:
3111   //       - lock closest ecal cluster
3112   //       - cut link to farthest hcal cluster, if more than 1.
3113
3114   vector<bool> deadArea(elements.size(), false);
3115
3116   // vectors to store element indices to ho, hcal and ecal elements, will be filled by elementLoop()
3117   ElementIndices inds;
3118
3119   elementLoop(block, linkData, elements, active, blockref, inds, deadArea);
3120
3121   // Reconstruct pfCandidate from HF (either EM-only, Had-only or both)
3122   // For phase2, process also pfblocks containing HF clusters and linked tracks
3123   if (!(inds.hfEmIs.empty() && inds.hfHadIs.empty())) {
3124     createCandidatesHF(block, linkData, elements, active, blockref, inds);
3125     if (inds.hcalIs.empty() && inds.ecalIs.empty())
3126       return;
3127     LogDebug("PFAlgo::processBlock")
3128         << "Block contains HF clusters, and also contains ECAL or HCAL clusters. Continue.\n"
3129         << block;
3130   }
3131
3132   createCandidatesHCAL(block, linkData, elements, active, blockref, inds, deadArea);
3133   // COLINFEB16: now dealing with the HCAL elements that are not linked to any track
3134   createCandidatesHCALUnlinked(block, linkData, elements, active, blockref, inds, deadArea);
3135   createCandidatesECAL(block, linkData, elements, active, blockref, inds, deadArea);
3136
3137   LogTrace("PFAlgo|processBlock") << "end of function PFAlgo::processBlock";
3138 }  // end processBlock
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   // Assume this particle is a charged Hadron
3150   double px = track.px();
3151   double py = track.py();
3152   double pz = track.pz();
3153   double energy = sqrt(track.p() * track.p() + 0.13957 * 0.13957);
3154
3155   LogTrace("PFAlgo|reconstructTrack") << "Reconstructing PF candidate from track of pT = " << track.pt()
3156                                       << " eta = " << track.eta() << " phi = " << track.phi() << " px = " << px
3157                                       << " py = " << py << " pz = " << pz << " energy = " << energy;
3158
3159   // Create a PF Candidate
3160   ::math::XYZTLorentzVector momentum(px, py, pz, energy);
3161   reco::PFCandidate::ParticleType particleType = reco::PFCandidate::h;
3162
3163   // Add it to the stack
3164   LogTrace("PFAlgo|reconstructTrack") << "Creating PFCandidate charge=" << charge << ", type=" << particleType
3165                                       << ", pt=" << momentum.pt() << ", eta=" << momentum.eta()
3166                                       << ", phi=" << momentum.phi();
3167   pfCandidates_->push_back(PFCandidate(charge, momentum, particleType));
3168   //Set vertex and stuff like this
3169   pfCandidates_->back().setVertex(trackRef->vertex());
3170   pfCandidates_->back().setTrackRef(trackRef);
3171   pfCandidates_->back().setPositionAtECALEntrance(eltTrack->positionAtECALEntrance());
3172   if (muonRef.isNonnull())
3173     pfCandidates_->back().setMuonRef(muonRef);
3174
3175   //Set time
3176   if (elt.isTimeValid())
3177     pfCandidates_->back().setTime(elt.time(), elt.timeError());
3178
3179   //OK Now try to reconstruct the particle as a muon
3180   bool isMuon = pfmu_->reconstructMuon(pfCandidates_->back(), muonRef, allowLoose);
3181   bool isFromDisp = isFromSecInt(elt, "secondary");
3182
3183   if ((!isMuon) && isFromDisp) {
3184     double dpt = trackRef->ptError();
3185     double dptRel = dpt / trackRef->pt() * 100;
3186     //If the track is ill measured it is better to not refit it, since the track information probably would not be used.
3187     //In the PFAlgo we use the trackref information. If the track error is too big the refitted information might be very different
3188     // from the not refitted one.
3189     if (dptRel < dptRel_DispVtx_) {
3190       LogTrace("PFAlgo|reconstructTrack")
3191           << "Not refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy;
3192       //reco::TrackRef trackRef = eltTrack->trackRef();
3193       reco::PFDisplacedVertexRef vRef =
3194           eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef();
3195       reco::Track trackRefit = vRef->refittedTrack(trackRef);
3196       //change the momentum with the refitted track
3197       ::math::XYZTLorentzVector momentum(
3198           trackRefit.px(), trackRefit.py(), trackRefit.pz(), sqrt(trackRefit.p() * trackRefit.p() + 0.13957 * 0.13957));
3199       LogTrace("PFAlgo|reconstructTrack")
3200           << "Refitted px = " << px << " py = " << py << " pz = " << pz << " energy = " << energy;
3201     }
3202     pfCandidates_->back().setFlag(reco::PFCandidate::T_FROM_DISP, true);
3203     pfCandidates_->back().setDisplacedVertexRef(
3204         eltTrack->displacedVertexRef(reco::PFBlockElement::T_FROM_DISP)->displacedVertexRef(),
3205         reco::PFCandidate::T_FROM_DISP);
3206   }
3207
3208   // do not label as primary a track which would be recognised as a muon. A muon cannot produce NI. It is with high probability a fake
3209   if (isFromSecInt(elt, "primary") && !isMuon) {
3210     pfCandidates_->back().setFlag(reco::PFCandidate::T_TO_DISP, true);
3211     pfCandidates_->back().setDisplacedVertexRef(
3212         eltTrack->displacedVertexRef(reco::PFBlockElement::T_TO_DISP)->displacedVertexRef(),
3213         reco::PFCandidate::T_TO_DISP);
3214   }
3215
3216   // returns index to the newly created PFCandidate
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   // need to convert the ::math::XYZPoint data member of the PFCluster class=
3234   // to a displacement vector:
3235
3236   // Transform particleX,Y,Z to a position at ECAL/HCAL entrance
3237   double factor = 1.;
3238   if (useDirection) {
3239     switch (cluster.layer()) {
3240       case PFLayer::ECAL_BARREL:
3241       case PFLayer::HCAL_BARREL1:
3242         factor = std::sqrt(cluster.position().Perp2() / (particleX * particleX + particleY * particleY));
3243         break;
3244       case PFLayer::ECAL_ENDCAP:
3245       case PFLayer::HCAL_ENDCAP:
3246       case PFLayer::HF_HAD:
3247       case PFLayer::HF_EM:
3248         factor = cluster.position().Z() / particleZ;
3249         break;
3250       default:
3251         assert(0);
3252     }
3253   }
3254   //MIKE First of all let's check if we have vertex.
3255   ::math::XYZPoint vertexPos;
3256   if (useVertices_)
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   //::math::XYZVector clusterPos( cluster.position().X(), cluster.position().Y(),cluster.position().Z() );
3268   //::math::XYZVector particleDirection ( particleX, particleY, particleZ );
3269
3270   clusterPos = useDirection ? particleDirection.Unit() : clusterPos.Unit();
3271   clusterPos *= particleEnergy;
3272
3273   // clusterPos is now a vector along the cluster direction,
3274   // with a magnitude equal to the cluster energy.
3275
3276   double mass = 0;
3277   ROOT::Math::LorentzVector<ROOT::Math::PxPyPzM4D<double>> momentum(
3278       clusterPos.X(), clusterPos.Y(), clusterPos.Z(), mass);
3279   // mathcore is a piece of #$%
3280   ::math::XYZTLorentzVector tmp;
3281   // implicit constructor not allowed
3282   tmp = momentum;
3283
3284   // Charge
3285   int charge = 0;
3286
3287   // Type
3288   switch (cluster.layer()) {
3289     case PFLayer::ECAL_BARREL:
3290     case PFLayer::ECAL_ENDCAP:
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   // The pf candidate
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   // The position at ECAL entrance (well: watch out, it is not true
3313   // for HCAL clusters... to be fixed)
3314   pfCandidates_->back().setPositionAtECALEntrance(
3315       ::math::XYZPointF(cluster.position().X(), cluster.position().Y(), cluster.position().Z()));
3316
3317   //Set the cnadidate Vertex
3318   pfCandidates_->back().setVertex(vertexPos);
3319
3320   // depth info
3321   setHcalDepthInfo(pfCandidates_->back(), cluster);
3322
3323   //*TODO* cluster time is not reliable at the moment, so only use track timing
3324
3325   LogTrace("PFAlgo|reconstructCluster") << "** candidate: " << pfCandidates_->back();
3326
3327   // returns index to the newly created PFCandidate
3328   return pfCandidates_->size() - 1;
3329 }
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 //GMA need the followign two for HO also
3360
3361 double PFAlgo::neutralHadronEnergyResolution(double clusterEnergyHCAL, double eta) const {
3362   // Add a protection
3363   clusterEnergyHCAL = std::max(clusterEnergyHCAL, 1.);
3364
3365   double resol = fabs(eta) < 1.48 ? sqrt(1.02 * 1.02 / clusterEnergyHCAL + 0.065 * 0.065)
3366                                   : sqrt(1.20 * 1.20 / clusterEnergyHCAL + 0.028 * 0.028);
3367
3368   return resol;
3369 }
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   // Add a protection
3380   clusterEnergyHF = std::max(clusterEnergyHF, 1.);
3381
3382   double resol =
3383       sqrt(resolHF_square_[0] / clusterEnergyHF + resolHF_square_[1] + resolHF_square_[2] / pow(clusterEnergyHF, 2));
3384   // 0: stochastic term, 1: constant term, 2: noise term
3385   // Note: resolHF_square_[0,1,2] should be already squared
3386
3387   return resol;
3388 }
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   // Find all PS clusters with type psElement associated to ECAL cluster iEcal,
3433   // within all PFBlockElement "elements" of a given PFBlock "block"
3434   // psElement can be reco::PFBlockElement::PS1 or reco::PFBlockElement::PS2
3435   // Returns a vector of PS cluster energies, and updates the "active" vector.
3436
3437   // Find all PS clusters linked to the iEcal cluster
3438   std::multimap<double, unsigned> sortedPS;
3439   block.associatedElements(iEcal, linkData, sortedPS, psElementType, reco::PFBlock::LINKTEST_ALL);
3440
3441   // Loop over these PS clusters
3442   double totalPS = 0.;
3443   for (auto const& ps : sortedPS) {
3444     // CLuster index and distance to iEcal
3445     unsigned iPS = ps.second;
3446     // double distPS = ps.first;
3447
3448     // Ignore clusters already in use
3449     if (!active[iPS])
3450       continue;
3451
3452     // Check that this cluster is not closer to another ECAL cluster
3453     std::multimap<double, unsigned> sortedECAL;
3454     block.associatedElements(iPS, linkData, sortedECAL, reco::PFBlockElement::ECAL, reco::PFBlock::LINKTEST_ALL);
3455     unsigned jEcal = sortedECAL.begin()->second;
3456     if (jEcal != iEcal)
3457       continue;
3458
3459     // Update PS energy
3460     PFBlockElement::Type pstype = elements[iPS].type();
3461     assert(pstype == psElementType);
3462     PFClusterRef psclusterref = elements[iPS].clusterRef();
3463     assert(!psclusterref.isNull());
3464     totalPS += psclusterref->energy();
3465     psEne[0] += psclusterref->energy();
3466     active[iPS] = false;
3467   }
3468 }
3469
3470 bool PFAlgo::isFromSecInt(const reco::PFBlockElement& eTrack, string order) const {
3471   reco::PFBlockElement::TrackType T_TO_DISP = reco::PFBlockElement::T_TO_DISP;
3472   reco::PFBlockElement::TrackType T_FROM_DISP = reco::PFBlockElement::T_FROM_DISP;
3473   //  reco::PFBlockElement::TrackType T_FROM_GAMMACONV = reco::PFBlockElement::T_FROM_GAMMACONV;
3474   reco::PFBlockElement::TrackType T_FROM_V0 = reco::PFBlockElement::T_FROM_V0;
3475
3476   bool bPrimary = (order.find("primary") != string::npos);
3477   bool bSecondary = (order.find("secondary") != string::npos);
3478   bool bAll = (order.find("all") != string::npos);
3479
3480   bool isToDisp = usePFNuclearInteractions_ && eTrack.trackType(T_TO_DISP);
3481   bool isFromDisp = usePFNuclearInteractions_ && eTrack.trackType(T_FROM_DISP);
3482
3483   if (bPrimary && isToDisp)
3484     return true;
3485   if (bSecondary && isFromDisp)
3486     return true;
3487   if (bAll && (isToDisp || isFromDisp))
3488     return true;
3489
3490   //   bool isFromConv = usePFConversions_ && eTrack.trackType(T_FROM_GAMMACONV);
3491
3492   //   if ((bAll || bSecondary)&& isFromConv) return true;
3493
3494   bool isFromDecay = (bAll || bSecondary) && usePFDecays_ && eTrack.trackType(T_FROM_V0);
3495
3496   return isFromDecay;
3497 }
3498
3499 void PFAlgo::postCleaning() {
3500   //Compute met and met significance (met/sqrt(SumEt))
3501   double metX = 0.;
3502   double metY = 0.;
3503   double sumet = 0;
3504   std::vector<unsigned int> pfCandidatesToBeRemoved;
3505   for (auto const& pfc : *pfCandidates_) {
3506     metX += pfc.px();
3507     metY += pfc.py();
3508     sumet += pfc.pt();
3509   }
3510   double met2 = metX * metX + metY * metY;
3511   // Select events with large MET significance.
3512   double significance = std::sqrt(met2 / sumet);
3513   double significanceCor = significance;
3514   if (significance > minSignificance_) {
3515     double metXCor = metX;
3516     double metYCor = metY;
3517     double sumetCor = sumet;
3518     double met2Cor = met2;
3519     double deltaPhi = 3.14159;
3520     double deltaPhiPt = 100.;
3521     bool next = true;
3522     unsigned iCor = 1E9;
3523
3524     // Find the HF candidate with the largest effect on the MET
3525     while (next) {
3526       double metReduc = -1.;
3527       // Loop on the candidates
3528       for (unsigned i = 0; i < pfCandidates_->size(); ++i) {
3529         const PFCandidate& pfc = (*pfCandidates_)[i];
3530
3531         // Check that the pfCandidate is in the HF
3532         if (pfc.particleId() != reco::PFCandidate::h_HF && pfc.particleId() != reco::PFCandidate::egamma_HF)
3533           continue;
3534
3535         // Check if has meaningful pt
3536         if (pfc.pt() < minHFCleaningPt_)
3537           continue;
3538
3539         // Check that it is  not already scheculed to be cleaned
3540         const bool skip = std::any_of(
3541             pfCandidatesToBeRemoved.begin(), pfCandidatesToBeRemoved.end(), [&](unsigned int j) { return i == j; });
3542         if (skip)
3543           continue;
3544
3545         // Check that the pt and the MET are aligned
3546         deltaPhi = std::acos((metX * pfc.px() + metY * pfc.py()) / (pfc.pt() * std::sqrt(met2)));
3547         deltaPhiPt = deltaPhi * pfc.pt();
3548         if (deltaPhiPt > maxDeltaPhiPt_)
3549           continue;
3550
3551         // Now remove the candidate from the MET
3552         double metXInt = metX - pfc.px();
3553         double metYInt = metY - pfc.py();
3554         double sumetInt = sumet - pfc.pt();
3555         double met2Int = metXInt * metXInt + metYInt * metYInt;
3556         if (met2Int < met2Cor) {
3557           metXCor = metXInt;
3558           metYCor = metYInt;
3559           metReduc = (met2 - met2Int) / met2Int;
3560           met2Cor = met2Int;
3561           sumetCor = sumetInt;
3562           significanceCor = std::sqrt(met2Cor / sumetCor);
3563           iCor = i;
3564         }
3565       }
3566       //
3567       // If the MET must be significanly reduced, schedule the candidate to be cleaned
3568       if (metReduc > minDeltaMet_) {
3569         pfCandidatesToBeRemoved.push_back(iCor);
3570         metX = metXCor;
3571         metY = metYCor;
3572         sumet = sumetCor;
3573         met2 = met2Cor;
3574       } else {
3575         // Otherwise just stop the loop
3576         next = false;
3577       }
3578     }
3579     //
3580     // The significance must be significantly reduced to indeed clean the candidates
3581     if (significance - significanceCor > minSignificanceReduction_ && significanceCor < maxSignificance_) {
3582       edm::LogInfo("PFAlgo|postCleaning") << "Significance reduction = " << significance << " -> " << significanceCor
3583                                           << " = " << significanceCor - significance;
3584       for (unsigned int toRemove : pfCandidatesToBeRemoved) {
3585         edm::LogInfo("PFAlgo|postCleaning") << "Removed : " << (*pfCandidates_)[toRemove];
3586         pfCleanedCandidates_.push_back((*pfCandidates_)[toRemove]);
3587         (*pfCandidates_)[toRemove].rescaleMomentum(1E-6);
3588         //reco::PFCandidate::ParticleType unknown = reco::PFCandidate::X;
3589         //(*pfCandidates_)[toRemove].setParticleType(unknown);
3590       }
3591     }
3592   }  //significance
3593 }  //postCleaning
3594
3595 void PFAlgo::checkCleaning(const reco::PFRecHitCollection& cleanedHits) {
3596   // No hits to recover, leave.
3597   if (cleanedHits.empty())
3598     return;
3599
3600   //Compute met and met significance (met/sqrt(SumEt))
3601   double metX = 0.;
3602   double metY = 0.;
3603   double sumet = 0;
3604   std::vector<unsigned int> hitsToBeAdded;
3605   for (auto const& pfc : *pfCandidates_) {
3606     metX += pfc.px();
3607     metY += pfc.py();
3608     sumet += pfc.pt();
3609   }
3610   double met2 = metX * metX + metY * metY;
3611   double met2_Original = met2;
3612   // Select events with large MET significance.
3613   // double significance = std::sqrt(met2/sumet);
3614   // double significanceCor = significance;
3615   double metXCor = metX;
3616   double metYCor = metY;
3617   double sumetCor = sumet;
3618   double met2Cor = met2;
3619   bool next = true;
3620   unsigned iCor = 1E9;
3621
3622   // Find the cleaned hit with the largest effect on the MET
3623   while (next) {
3624     double metReduc = -1.;
3625     // Loop on the candidates
3626     for (unsigned i = 0; i < cleanedHits.size(); ++i) {
3627       const PFRecHit& hit = cleanedHits[i];
3628       double length = std::sqrt(hit.position().mag2());
3629       double px = hit.energy() * hit.position().x() / length;
3630       double py = hit.energy() * hit.position().y() / length;
3631       double pt = std::sqrt(px * px + py * py);
3632
3633       // Check that it is  not already scheculed to be cleaned
3634       bool skip = false;
3635       for (unsigned int hitIdx : hitsToBeAdded) {
3636         if (i == hitIdx)
3637           skip = true;
3638         if (skip)
3639           break;
3640       }
3641       if (skip)
3642         continue;
3643
3644       // Now add the candidate to the MET
3645       double metXInt = metX + px;
3646       double metYInt = metY + py;
3647       double sumetInt = sumet + pt;
3648       double met2Int = metXInt * metXInt + metYInt * metYInt;
3649
3650       // And check if it could contribute to a MET reduction
3651       if (met2Int < met2Cor) {
3652         metXCor = metXInt;
3653         metYCor = metYInt;
3654         metReduc = (met2 - met2Int) / met2Int;
3655         met2Cor = met2Int;
3656         sumetCor = sumetInt;
3657         // significanceCor = std::sqrt(met2Cor/sumetCor);
3658         iCor = i;
3659       }
3660     }
3661     //
3662     // If the MET must be significanly reduced, schedule the candidate to be added
3663     //
3664     if (metReduc > minDeltaMet_) {
3665       hitsToBeAdded.push_back(iCor);
3666       metX = metXCor;
3667       metY = metYCor;
3668       sumet = sumetCor;
3669       met2 = met2Cor;
3670     } else {
3671       // Otherwise just stop the loop
3672       next = false;
3673     }
3674   }
3675   //
3676   // At least 10 GeV MET reduction
3677   if (std::sqrt(met2_Original) - std::sqrt(met2) > 5.) {
3678     LogTrace("PFAlgo|checkCleaning") << hitsToBeAdded.size() << " hits were re-added ";
3679     LogTrace("PFAlgo|checkCleaning") << "MET reduction = " << std::sqrt(met2_Original) << " -> " << std::sqrt(met2Cor)
3680                                      << " = " << std::sqrt(met2Cor) - std::sqrt(met2_Original);
3681     LogTrace("PFAlgo|checkCleaning") << "Added after cleaning check : ";
3682     for (unsigned int hitIdx : hitsToBeAdded) {
3683       const PFRecHit& hit = cleanedHits[hitIdx];
3684       PFCluster cluster(hit.layer(), hit.energy(), hit.position().x(), hit.position().y(), hit.position().z());
3685       reconstructCluster(cluster, hit.energy());
3686       LogTrace("PFAlgo|checkCleaning") << pfCandidates_->back() << ". time = " << hit.time();
3687     }
3688   }
3689 }  //PFAlgo::checkCleaning