Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2022-06-10 01:54:01

0001 #include "L1Trigger/Phase2L1ParticleFlow/interface/RegionMapper.h"
0002 #include "FWCore/Framework/interface/Event.h"
0003 #include "DataFormats/Common/interface/RefToPtr.h"
0004 
0005 #include "DataFormats/L1TCorrelator/interface/TkElectron.h"
0006 #include "DataFormats/L1TCorrelator/interface/TkElectronFwd.h"
0007 #include "DataFormats/L1Trigger/interface/EGamma.h"
0008 #include "DataFormats/L1TCorrelator/interface/TkEm.h"
0009 #include "DataFormats/L1TCorrelator/interface/TkEmFwd.h"
0010 
0011 using namespace l1tpf_impl;
0012 
0013 RegionMapper::RegionMapper(const edm::ParameterSet &iConfig) : useRelativeRegionalCoordinates_(false) {
0014   if (iConfig.existsAs<std::vector<edm::ParameterSet>>("regions")) {
0015     useRelativeRegionalCoordinates_ = iConfig.getParameter<bool>("useRelativeRegionalCoordinates");
0016     for (const edm::ParameterSet &preg : iConfig.getParameter<std::vector<edm::ParameterSet>>("regions")) {
0017       std::vector<double> etaBoundaries = preg.getParameter<std::vector<double>>("etaBoundaries");
0018       unsigned int phiSlices = preg.getParameter<uint32_t>("phiSlices");
0019       float etaExtra = preg.getParameter<double>("etaExtra");
0020       float phiExtra = preg.getParameter<double>("phiExtra");
0021       float phiWidth = 2 * M_PI / phiSlices;
0022       unsigned int ncalomax = 0, nemcalomax = 0, ntrackmax = 0, nmuonmax = 0, npfmax = 0, npuppimax = 0;
0023       if (preg.existsAs<uint32_t>("caloNMax"))
0024         ncalomax = preg.getParameter<uint32_t>("caloNMax");
0025       if (preg.existsAs<uint32_t>("emcaloNMax"))
0026         nemcalomax = preg.getParameter<uint32_t>("emcaloNMax");
0027       if (preg.existsAs<uint32_t>("trackNMax"))
0028         ntrackmax = preg.getParameter<uint32_t>("trackNMax");
0029       if (preg.existsAs<uint32_t>("muonNMax"))
0030         nmuonmax = preg.getParameter<uint32_t>("muonNMax");
0031       if (preg.existsAs<uint32_t>("pfNMax"))
0032         npfmax = preg.getParameter<uint32_t>("pfNMax");
0033       if (preg.existsAs<uint32_t>("puppiNMax"))
0034         npuppimax = preg.getParameter<uint32_t>("puppiNMax");
0035       for (unsigned int ieta = 0, neta = etaBoundaries.size() - 1; ieta < neta; ++ieta) {
0036         for (unsigned int iphi = 0; iphi < phiSlices; ++iphi) {
0037           //float phiCenter = (iphi + 0.5) * phiWidth - M_PI;
0038           float phiCenter = reco::reduceRange(iphi * phiWidth);  //align with L1 TrackFinder phi sector indexing for now
0039           regions_.push_back(Region(etaBoundaries[ieta],
0040                                     etaBoundaries[ieta + 1],
0041                                     phiCenter,
0042                                     phiWidth,
0043                                     etaExtra,
0044                                     phiExtra,
0045                                     useRelativeRegionalCoordinates_,
0046                                     ncalomax,
0047                                     nemcalomax,
0048                                     ntrackmax,
0049                                     nmuonmax,
0050                                     npfmax,
0051                                     npuppimax));
0052         }
0053       }
0054     }
0055     std::string trackRegionMode = "TrackAssoMode::any";
0056     if (iConfig.existsAs<std::string>("trackRegionMode"))
0057       trackRegionMode = iConfig.getParameter<std::string>("trackRegionMode");
0058     if (trackRegionMode == "atVertex")
0059       trackRegionMode_ = TrackAssoMode::atVertex;
0060     else if (trackRegionMode == "atCalo")
0061       trackRegionMode_ = TrackAssoMode::atCalo;
0062     else if (trackRegionMode == "any")
0063       trackRegionMode_ = TrackAssoMode::any;
0064     else
0065       throw cms::Exception(
0066           "Configuration",
0067           "Unsupported value for trackRegionMode: " + trackRegionMode + " (allowed are 'atVertex', 'atCalo', 'any')");
0068   } else {
0069     // start off with a dummy region
0070     unsigned int ncalomax = 0, nemcalomax = 0, ntrackmax = 0, nmuonmax = 0, npfmax = 0, npuppimax = 0;
0071     regions_.emplace_back(-5.5,
0072                           5.5,
0073                           0,
0074                           2 * M_PI,
0075                           0.5,
0076                           0.5,
0077                           useRelativeRegionalCoordinates_,
0078                           ncalomax,
0079                           nemcalomax,
0080                           ntrackmax,
0081                           nmuonmax,
0082                           npfmax,
0083                           npuppimax);
0084   }
0085 }
0086 
0087 void RegionMapper::clear() {
0088   for (Region &r : regions_)
0089     r.zero();
0090   clusterRefMap_.clear();
0091   trackRefMap_.clear();
0092 }
0093 
0094 void RegionMapper::addTrack(const l1t::PFTrack &t) {
0095   // now let's be optimistic and make things very simple
0096   // we propagate in floating point the track to the calo
0097   // we add the track to the region corresponding to its vertex (eta,phi) coordinates AND its (eta,phi) calo coordinates
0098   for (Region &r : regions_) {
0099     bool inside = true;
0100     switch (trackRegionMode_) {
0101       case TrackAssoMode::atVertex:
0102         inside = r.contains(t.eta(), t.phi());
0103         break;
0104       case TrackAssoMode::atCalo:
0105         inside = r.contains(t.caloEta(), t.caloPhi());
0106         break;
0107       case TrackAssoMode::any:
0108         inside = r.contains(t.eta(), t.phi()) || r.contains(t.caloEta(), t.caloPhi());
0109         break;
0110     }
0111     if (inside) {
0112       PropagatedTrack prop;
0113       prop.fillInput(t.pt(), r.localEta(t.eta()), r.localPhi(t.phi()), t.charge(), t.vertex().Z(), t.quality(), &t);
0114       prop.fillPropagated(t.pt(),
0115                           t.trkPtError(),
0116                           t.caloPtError(),
0117                           r.localEta(t.caloEta()),
0118                           r.localPhi(t.caloPhi()),
0119                           t.quality(),
0120                           t.isMuon());
0121       prop.hwStubs = t.nStubs();
0122       prop.hwChi2 = round(t.chi2() * 10);
0123       r.track.push_back(prop);
0124     }
0125   }
0126 }
0127 void RegionMapper::addTrack(const l1t::PFTrack &t, l1t::PFTrackRef ref) {
0128   addTrack(t);
0129   trackRefMap_[&t] = ref;
0130 }
0131 
0132 void RegionMapper::addMuon(const l1t::Muon &mu) {
0133   // now let's be optimistic and make things very simple
0134   // we don't propagate anything
0135   for (Region &r : regions_) {
0136     if (r.contains(mu.eta(), mu.phi())) {
0137       Muon prop;
0138       prop.fill(mu.pt(), r.localEta(mu.eta()), r.localPhi(mu.phi()), mu.charge(), mu.hwQual(), &mu);
0139       r.muon.push_back(prop);
0140     }
0141   }
0142 }
0143 
0144 void RegionMapper::addMuon(const l1t::TkMuon &mu) {
0145   // now let's be optimistic and make things very simple
0146   // we don't propagate anything
0147   for (Region &r : regions_) {
0148     if (r.contains(mu.eta(), mu.phi())) {
0149       Muon prop;
0150       prop.fill(mu.pt(), r.localEta(mu.eta()), r.localPhi(mu.phi()), mu.charge(), mu.hwQual());
0151       r.muon.push_back(prop);
0152     }
0153   }
0154 }
0155 
0156 void RegionMapper::addCalo(const l1t::PFCluster &p) {
0157   if (p.pt() == 0)
0158     return;
0159   for (Region &r : regions_) {
0160     if (r.contains(p.eta(), p.phi())) {
0161       CaloCluster calo;
0162       calo.fill(p.pt(), p.emEt(), p.ptError(), r.localEta(p.eta()), r.localPhi(p.phi()), p.isEM(), 0, &p);
0163       r.calo.push_back(calo);
0164     }
0165   }
0166 }
0167 void RegionMapper::addCalo(const l1t::PFCluster &p, l1t::PFClusterRef ref) {
0168   addCalo(p);
0169   clusterRefMap_[&p] = ref;
0170 }
0171 
0172 void RegionMapper::addEmCalo(const l1t::PFCluster &p) {
0173   if (p.pt() == 0)
0174     return;
0175   for (Region &r : regions_) {
0176     if (r.contains(p.eta(), p.phi())) {
0177       CaloCluster calo;
0178       calo.fill(p.pt(), p.emEt(), p.ptError(), r.localEta(p.eta()), r.localPhi(p.phi()), p.isEM(), p.hwQual(), &p);
0179       r.emcalo.push_back(calo);
0180     }
0181   }
0182 }
0183 void RegionMapper::addEmCalo(const l1t::PFCluster &p, l1t::PFClusterRef ref) {
0184   addEmCalo(p);
0185   clusterRefMap_[&p] = ref;
0186 }
0187 
0188 std::unique_ptr<l1t::PFCandidateCollection> RegionMapper::fetch(bool puppi, float ptMin) const {
0189   auto ret = std::make_unique<l1t::PFCandidateCollection>();
0190   for (const Region &r : regions_) {
0191     for (const PFParticle &p : (puppi ? r.puppi : r.pf)) {
0192       bool inside = true;
0193       switch (trackRegionMode_) {
0194         case TrackAssoMode::atVertex:
0195           inside = r.fiducialLocal(p.floatVtxEta(), p.floatVtxPhi());
0196           break;
0197         case TrackAssoMode::atCalo:
0198           inside = r.fiducialLocal(p.floatEta(), p.floatPhi());
0199           break;
0200         case TrackAssoMode::any:
0201           inside = r.fiducialLocal(p.floatVtxEta(), p.floatVtxPhi());
0202           break;  // WARNING: this may not be the best choice
0203       }
0204       if (!inside)
0205         continue;
0206       if (p.floatPt() > ptMin) {
0207         reco::Particle::PolarLorentzVector p4(
0208             p.floatPt(), r.globalEta(p.floatVtxEta()), r.globalPhi(p.floatVtxPhi()), 0.13f);
0209         ret->emplace_back(
0210             l1t::PFCandidate::ParticleType(p.hwId), p.intCharge(), p4, p.floatPuppiW(), p.hwPt, p.hwEta, p.hwPhi);
0211         ret->back().setZ0(p.floatDZ());
0212         ret->back().setHwZ0(p.track.hwZ0);
0213         ret->back().setStatus(p.hwStatus);
0214         if (p.cluster.src) {
0215           auto match = clusterRefMap_.find(p.cluster.src);
0216           if (match == clusterRefMap_.end()) {
0217             throw cms::Exception("CorruptData") << "Invalid cluster pointer in PF candidate id " << p.hwId << " pt "
0218                                                 << p4.pt() << " eta " << p4.eta() << " phi " << p4.phi();
0219           }
0220           ret->back().setPFCluster(match->second);
0221         }
0222         if (p.track.src) {
0223           auto match = trackRefMap_.find(p.track.src);
0224           if (match == trackRefMap_.end()) {
0225             throw cms::Exception("CorruptData") << "Invalid track pointer in PF candidate id " << p.hwId << " pt "
0226                                                 << p4.pt() << " eta " << p4.eta() << " phi " << p4.phi();
0227           }
0228           ret->back().setPFTrack(match->second);
0229         }
0230       }
0231     }
0232   }
0233   return ret;
0234 }
0235 
0236 std::unique_ptr<l1t::PFCandidateCollection> RegionMapper::fetchCalo(float ptMin, bool emcalo) const {
0237   auto ret = std::make_unique<l1t::PFCandidateCollection>();
0238   for (const Region &r : regions_) {
0239     for (const CaloCluster &p : (emcalo ? r.emcalo : r.calo)) {
0240       if (!r.fiducialLocal(p.floatEta(), p.floatPhi()))
0241         continue;
0242       if (p.floatPt() > ptMin) {
0243         reco::Particle::PolarLorentzVector p4(p.floatPt(), r.globalEta(p.floatEta()), r.globalPhi(p.floatPhi()), 0.13f);
0244         l1t::PFCandidate::ParticleType kind =
0245             (p.isEM || emcalo) ? l1t::PFCandidate::Photon : l1t::PFCandidate::NeutralHadron;
0246         ret->emplace_back(kind, /*charge=*/0, p4, /*puppiW=*/1, p.hwPt, p.hwEta, p.hwPhi);
0247         if (p.src) {
0248           auto match = clusterRefMap_.find(p.src);
0249           if (match == clusterRefMap_.end()) {
0250             throw cms::Exception("CorruptData")
0251                 << "Invalid cluster pointer in cluster pt " << p4.pt() << " eta " << p4.eta() << " phi " << p4.phi();
0252           }
0253           ret->back().setPFCluster(match->second);
0254         }
0255       }
0256     }
0257   }
0258   return ret;
0259 }
0260 
0261 std::unique_ptr<l1t::PFCandidateCollection> RegionMapper::fetchTracks(float ptMin, bool fromPV) const {
0262   auto ret = std::make_unique<l1t::PFCandidateCollection>();
0263   for (const Region &r : regions_) {
0264     for (const PropagatedTrack &p : r.track) {
0265       if (fromPV && !p.fromPV)
0266         continue;
0267       bool inside = true;
0268       switch (trackRegionMode_) {
0269         case TrackAssoMode::atVertex:
0270           inside = r.fiducialLocal(p.floatVtxEta(), p.floatVtxPhi());
0271           break;
0272         case TrackAssoMode::atCalo:
0273           inside = r.fiducialLocal(p.floatEta(), p.floatPhi());
0274           break;
0275         case TrackAssoMode::any:
0276           inside = r.fiducialLocal(p.floatVtxEta(), p.floatVtxPhi());
0277           break;  // WARNING: this may not be the best choice
0278       }
0279       if (!inside)
0280         continue;
0281       if (p.floatPt() > ptMin) {
0282         reco::Particle::PolarLorentzVector p4(
0283             p.floatVtxPt(), r.globalEta(p.floatVtxEta()), r.globalPhi(p.floatVtxPhi()), 0.13f);
0284         l1t::PFCandidate::ParticleType kind = p.muonLink ? l1t::PFCandidate::Muon : l1t::PFCandidate::ChargedHadron;
0285         ret->emplace_back(kind, p.intCharge(), p4, /*puppiW=*/float(p.fromPV), p.hwPt, p.hwEta, p.hwPhi);
0286         ret->back().setZ0(p.floatDZ());
0287         ret->back().setHwZ0(p.hwZ0);
0288         if (p.src) {
0289           auto match = trackRefMap_.find(p.src);
0290           if (match == trackRefMap_.end()) {
0291             throw cms::Exception("CorruptData")
0292                 << "Invalid track pointer in PF track  pt " << p4.pt() << " eta " << p4.eta() << " phi " << p4.phi();
0293           }
0294           ret->back().setPFTrack(match->second);
0295         }
0296       }
0297     }
0298   }
0299   return ret;
0300 }
0301 
0302 void RegionMapper::putEgObjects(edm::Event &iEvent,
0303                                 const bool writeEgSta,
0304                                 const std::string &egLablel,
0305                                 const std::string &tkEmLabel,
0306                                 const std::string &tkEleLabel,
0307                                 const float ptMin) const {
0308   auto egs = std::make_unique<BXVector<l1t::EGamma>>();
0309   auto tkems = std::make_unique<l1t::TkEmCollection>();
0310   auto tkeles = std::make_unique<l1t::TkElectronCollection>();
0311 
0312   edm::RefProd<BXVector<l1t::EGamma>> ref_egs;
0313   if (writeEgSta)
0314     ref_egs = iEvent.getRefBeforePut<BXVector<l1t::EGamma>>(egLablel);
0315 
0316   edm::Ref<BXVector<l1t::EGamma>>::key_type idx = 0;
0317 
0318   for (const Region &r : regions_) {
0319     for (const auto &egphoton : r.egphotons) {
0320       if (egphoton.floatPt() < ptMin)
0321         continue;
0322 
0323       if (!r.fiducialLocal(egphoton.floatEta(), egphoton.floatPhi()))
0324         continue;
0325 
0326       edm::Ref<BXVector<l1t::EGamma>> reg;
0327       auto mom = reco::Candidate::PolarLorentzVector(
0328           egphoton.floatPt(), r.globalEta(egphoton.floatEta()), r.globalPhi(egphoton.floatPhi()), 0.);
0329       if (writeEgSta) {
0330         l1t::EGamma eg(mom);
0331         eg.setHwQual(egphoton.hwQual);
0332         egs->push_back(0, eg);
0333         reg = edm::Ref<BXVector<l1t::EGamma>>(ref_egs, idx++);
0334       } else {
0335         auto egptr = egphoton.cluster.src->constituentsAndFractions()[0].first;
0336         reg = edm::Ref<BXVector<l1t::EGamma>>(egptr.id(), dynamic_cast<const l1t::EGamma *>(egptr.get()), egptr.key());
0337       }
0338 
0339       l1t::TkEm tkem(reco::Candidate::LorentzVector(mom), reg, egphoton.floatIso(), egphoton.floatIsoPV());
0340       tkem.setHwQual(egphoton.hwQual);
0341       tkem.setPFIsol(egphoton.floatPFIso());
0342       tkem.setPFIsolPV(egphoton.floatPFIsoPV());
0343       tkems->push_back(tkem);
0344 
0345       if (egphoton.ele_idx == -1)
0346         continue;
0347 
0348       const auto &egele = r.egeles[egphoton.ele_idx];
0349 
0350       if (!r.fiducialLocal(egele.floatEta(), egele.floatPhi()))
0351         continue;
0352 
0353       auto mom_ele = reco::Candidate::PolarLorentzVector(
0354           egele.floatPt(), r.globalEta(egele.floatEta()), r.globalPhi(egele.floatPhi()), 0.);
0355 
0356       l1t::TkElectron tkele(
0357           reco::Candidate::LorentzVector(mom_ele), reg, edm::refToPtr(egele.track.src->track()), egele.floatIso());
0358       tkele.setHwQual(egele.hwQual);
0359       tkele.setPFIsol(egele.floatPFIso());
0360       tkeles->push_back(tkele);
0361     }
0362   }
0363   if (writeEgSta)
0364     iEvent.put(std::move(egs), egLablel);
0365   iEvent.put(std::move(tkems), tkEmLabel);
0366   iEvent.put(std::move(tkeles), tkEleLabel);
0367 }
0368 
0369 std::pair<unsigned, unsigned> RegionMapper::totAndMaxInput(int type) const {
0370   unsigned ntot = 0, nmax = 0;
0371   for (const auto &r : regions_) {
0372     unsigned int ni = r.nInput(Region::InputType(type));
0373     ntot += ni;
0374     nmax = std::max(nmax, ni);
0375   }
0376   return std::make_pair(ntot, nmax);
0377 }
0378 
0379 std::unique_ptr<std::vector<unsigned>> RegionMapper::vecInput(int type) const {
0380   auto v = std::make_unique<std::vector<unsigned>>();
0381   for (const auto &r : regions_) {
0382     unsigned ni = r.nInput(Region::InputType(type));
0383     v->push_back(ni);
0384   }
0385   return v;
0386 }
0387 
0388 std::pair<unsigned, unsigned> RegionMapper::totAndMaxOutput(int type, bool puppi) const {
0389   unsigned ntot = 0, nmax = 0;
0390   for (const auto &r : regions_) {
0391     unsigned int ni = r.nOutput(Region::OutputType(type), puppi);
0392     ntot += ni;
0393     nmax = std::max(nmax, ni);
0394   }
0395   return std::make_pair(ntot, nmax);
0396 }
0397 
0398 std::unique_ptr<std::vector<unsigned>> RegionMapper::vecOutput(int type, bool puppi) const {
0399   auto v = std::make_unique<std::vector<unsigned>>();
0400   for (const auto &r : regions_) {
0401     unsigned ni = r.nOutput(Region::OutputType(type), puppi);
0402     v->push_back(ni);
0403   }
0404   return v;
0405 }