File indexing completed on 2022-02-28 01:32:08
0001
0002
0003
0004
0005
0006 #include <vector>
0007 #include <memory>
0008 #include <algorithm>
0009 #include <cmath>
0010
0011 #include "DataFormats/Common/interface/Handle.h"
0012 #include "DataFormats/Common/interface/Ref.h"
0013 #include "DataFormats/DetId/interface/DetId.h"
0014 #include "DataFormats/TrackReco/interface/Track.h"
0015 #include "DataFormats/HcalIsolatedTrack/interface/IsolatedPixelTrackCandidate.h"
0016 #include "DataFormats/Common/interface/TriggerResults.h"
0017
0018 #include "DataFormats/L1Trigger/interface/L1EmParticle.h"
0019 #include "DataFormats/L1Trigger/interface/L1JetParticle.h"
0020
0021 #include "DataFormats/L1Trigger/interface/Jet.h"
0022 #include "DataFormats/L1Trigger/interface/Tau.h"
0023 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutSetupFwd.h"
0024 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerReadoutRecord.h"
0025 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapRecord.h"
0026 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMapFwd.h"
0027 #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h"
0028 #include "DataFormats/HLTReco/interface/TriggerFilterObjectWithRefs.h"
0029
0030 #include "DataFormats/VertexReco/interface/Vertex.h"
0031 #include "DataFormats/VertexReco/interface/VertexFwd.h"
0032
0033
0034 #include "FWCore/Framework/interface/Frameworkfwd.h"
0035 #include "FWCore/Framework/interface/stream/EDProducer.h"
0036 #include "FWCore/Framework/interface/Event.h"
0037 #include "FWCore/Framework/interface/EventSetup.h"
0038 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0039 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0040 #include "FWCore/ParameterSet/interface/ConfigurationDescriptions.h"
0041 #include "FWCore/ParameterSet/interface/ParameterSetDescription.h"
0042 #include "FWCore/Utilities/interface/Exception.h"
0043 #include "FWCore/Utilities/interface/transform.h"
0044
0045
0046 #include "Math/GenVector/VectorUtil.h"
0047 #include "Math/GenVector/PxPyPzE4D.h"
0048 #include "DataFormats/Math/interface/deltaR.h"
0049
0050
0051 #include "MagneticField/VolumeBasedEngine/interface/VolumeBasedMagneticField.h"
0052 #include "MagneticField/Engine/interface/MagneticField.h"
0053 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0054
0055
0056 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0057 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0058 #include "Geometry/EcalAlgo/interface/EcalBarrelGeometry.h"
0059 #include "Geometry/EcalAlgo/interface/EcalEndcapGeometry.h"
0060 #include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h"
0061
0062
0063
0064 class IsolatedPixelTrackCandidateL1TProducer : public edm::stream::EDProducer<> {
0065 public:
0066 IsolatedPixelTrackCandidateL1TProducer(const edm::ParameterSet& ps);
0067 ~IsolatedPixelTrackCandidateL1TProducer() override;
0068
0069 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0070
0071 void beginRun(const edm::Run&, const edm::EventSetup&) override;
0072 void produce(edm::Event& evt, const edm::EventSetup& es) override;
0073
0074 double getDistInCM(double eta1, double phi1, double eta2, double phi2);
0075 std::pair<double, double> GetEtaPhiAtEcal(double etaIP, double phiIP, double pT, int charge, double vtxZ);
0076
0077 private:
0078 struct seedAtEC {
0079 seedAtEC(unsigned int i, bool f, double et, double fi) : index(i), ok(f), eta(et), phi(fi) {}
0080 unsigned int index;
0081 bool ok;
0082 double eta, phi;
0083 };
0084
0085 const edm::EDGetTokenT<trigger::TriggerFilterObjectWithRefs> tok_hlt_;
0086 const edm::EDGetTokenT<l1t::TauBxCollection> tok_l1_;
0087 const edm::EDGetTokenT<reco::VertexCollection> tok_vert_;
0088 const std::vector<edm::EDGetTokenT<reco::TrackCollection> > toks_pix_;
0089 const edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> tok_bFieldH_;
0090 const edm::ESGetToken<CaloGeometry, CaloGeometryRecord> tok_geom_;
0091
0092 const std::string bfield_;
0093 const double prelimCone_;
0094 const double pixelIsolationConeSizeAtEC_;
0095 const double vtxCutSeed_;
0096 const double vtxCutIsol_;
0097 const double tauAssocCone_;
0098 const double tauUnbiasCone_;
0099 const double minPTrackValue_;
0100 const double maxPForIsolationValue_;
0101 const double ebEtaBoundary_;
0102
0103
0104 double rEB_;
0105 double zEE_;
0106 double bfVal_;
0107 };
0108
0109 IsolatedPixelTrackCandidateL1TProducer::IsolatedPixelTrackCandidateL1TProducer(const edm::ParameterSet& config)
0110 : tok_hlt_(consumes<trigger::TriggerFilterObjectWithRefs>(config.getParameter<edm::InputTag>("L1GTSeedLabel"))),
0111 tok_l1_(consumes<l1t::TauBxCollection>(config.getParameter<edm::InputTag>("L1eTauJetsSource"))),
0112 tok_vert_(consumes<reco::VertexCollection>(config.getParameter<edm::InputTag>("VertexLabel"))),
0113 toks_pix_(
0114 edm::vector_transform(config.getParameter<std::vector<edm::InputTag> >("PixelTracksSources"),
0115 [this](edm::InputTag const& tag) { return consumes<reco::TrackCollection>(tag); })),
0116 tok_bFieldH_(esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>()),
0117 tok_geom_(esConsumes<CaloGeometry, CaloGeometryRecord, edm::Transition::BeginRun>()),
0118 bfield_(config.getParameter<std::string>("MagFieldRecordName")),
0119 prelimCone_(config.getParameter<double>("ExtrapolationConeSize")),
0120 pixelIsolationConeSizeAtEC_(config.getParameter<double>("PixelIsolationConeSizeAtEC")),
0121 vtxCutSeed_(config.getParameter<double>("MaxVtxDXYSeed")),
0122 vtxCutIsol_(config.getParameter<double>("MaxVtxDXYIsol")),
0123 tauAssocCone_(config.getParameter<double>("tauAssociationCone")),
0124 tauUnbiasCone_(config.getParameter<double>("tauUnbiasCone")),
0125 minPTrackValue_(config.getParameter<double>("minPTrack")),
0126 maxPForIsolationValue_(config.getParameter<double>("maxPTrackForIsolation")),
0127 ebEtaBoundary_(config.getParameter<double>("EBEtaBoundary")),
0128 rEB_(-1),
0129 zEE_(-1),
0130 bfVal_(0) {
0131
0132 produces<reco::IsolatedPixelTrackCandidateCollection>();
0133 }
0134
0135 IsolatedPixelTrackCandidateL1TProducer::~IsolatedPixelTrackCandidateL1TProducer() {}
0136
0137 void IsolatedPixelTrackCandidateL1TProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0138 edm::ParameterSetDescription desc;
0139 std::vector<edm::InputTag> tracksrc = {edm::InputTag("hltPixelTracks")};
0140 desc.add<edm::InputTag>("L1eTauJetsSource", edm::InputTag("hltGtStage2Digis", "Tau"));
0141 desc.add<double>("tauAssociationCone", 0.0);
0142 desc.add<double>("tauUnbiasCone", 1.2);
0143 desc.add<std::vector<edm::InputTag> >("PixelTracksSources", tracksrc);
0144 desc.add<double>("ExtrapolationConeSize", 1.0);
0145 desc.add<double>("PixelIsolationConeSizeAtEC", 40);
0146 desc.add<edm::InputTag>("L1GTSeedLabel", edm::InputTag("hltL1sV0SingleJet60"));
0147 desc.add<double>("MaxVtxDXYSeed", 101.0);
0148 desc.add<double>("MaxVtxDXYIsol", 101.0);
0149 desc.add<edm::InputTag>("VertexLabel", edm::InputTag("hltTrimmedPixelVertices"));
0150 desc.add<std::string>("MagFieldRecordName", "VolumeBasedMagneticField");
0151 desc.add<double>("minPTrack", 5.0);
0152 desc.add<double>("maxPTrackForIsolation", 3.0);
0153 desc.add<double>("EBEtaBoundary", 1.479);
0154 descriptions.add("isolPixelTrackProdL1T", desc);
0155 }
0156
0157 void IsolatedPixelTrackCandidateL1TProducer::beginRun(const edm::Run& run, const edm::EventSetup& theEventSetup) {
0158 const CaloGeometry* pG = &theEventSetup.getData(tok_geom_);
0159
0160 const double rad(dynamic_cast<const EcalBarrelGeometry*>(pG->getSubdetectorGeometry(DetId::Ecal, EcalBarrel))
0161 ->avgRadiusXYFrontFaceCenter());
0162 const double zz(dynamic_cast<const EcalEndcapGeometry*>(pG->getSubdetectorGeometry(DetId::Ecal, EcalEndcap))
0163 ->avgAbsZFrontFaceCenter());
0164
0165 rEB_ = rad;
0166 zEE_ = zz;
0167
0168 const MagneticField* vbfField = &theEventSetup.getData(tok_bFieldH_);
0169 const VolumeBasedMagneticField* vbfCPtr = dynamic_cast<const VolumeBasedMagneticField*>(&(*vbfField));
0170 GlobalVector BField = vbfCPtr->inTesla(GlobalPoint(0, 0, 0));
0171 bfVal_ = BField.mag();
0172 #ifdef EDM_ML_DEBUG
0173 edm::LogVerbatim("IsoTrack") << "rEB " << rEB_ << " zEE " << zEE_ << " B " << bfVal_;
0174 #endif
0175 }
0176
0177 void IsolatedPixelTrackCandidateL1TProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
0178 auto trackCollection = std::make_unique<reco::IsolatedPixelTrackCandidateCollection>();
0179
0180
0181 std::vector<reco::TrackRef> pixelTrackRefs;
0182
0183 for (unsigned int iPix = 0; iPix < toks_pix_.size(); iPix++) {
0184 const edm::Handle<reco::TrackCollection>& iPixCol = theEvent.getHandle(toks_pix_[iPix]);
0185 for (reco::TrackCollection::const_iterator pit = iPixCol->begin(); pit != iPixCol->end(); pit++) {
0186 pixelTrackRefs.push_back(reco::TrackRef(iPixCol, pit - iPixCol->begin()));
0187 }
0188 }
0189
0190 const edm::Handle<l1t::TauBxCollection>& l1eTauJets = theEvent.getHandle(tok_l1_);
0191
0192 const edm::Handle<reco::VertexCollection>& pVert = theEvent.getHandle(tok_vert_);
0193
0194 double ptTriggered = -10;
0195 double etaTriggered = -100;
0196 double phiTriggered = -100;
0197
0198 const edm::Handle<trigger::TriggerFilterObjectWithRefs>& l1trigobj = theEvent.getHandle(tok_hlt_);
0199
0200 std::vector<edm::Ref<l1t::TauBxCollection> > l1tauobjref;
0201 std::vector<edm::Ref<l1t::JetBxCollection> > l1jetobjref;
0202
0203
0204 l1trigobj->getObjects(trigger::TriggerL1Tau, l1tauobjref);
0205
0206 l1trigobj->getObjects(trigger::TriggerL1Jet, l1jetobjref);
0207
0208 for (const auto& p : l1tauobjref) {
0209 if (p->pt() > ptTriggered) {
0210 ptTriggered = p->pt();
0211 phiTriggered = p->phi();
0212 etaTriggered = p->eta();
0213 }
0214 }
0215 for (const auto& p : l1jetobjref) {
0216 if (p->pt() > ptTriggered) {
0217 ptTriggered = p->pt();
0218 phiTriggered = p->phi();
0219 etaTriggered = p->eta();
0220 }
0221 }
0222 #ifdef EDM_ML_DEBUG
0223 edm::LogVerbatim("IsoTrack") << "Sizes " << l1tauobjref.size() << ":" << l1jetobjref.size() << " Trig " << ptTriggered
0224 << ":" << etaTriggered << ":" << phiTriggered;
0225 #endif
0226 double drMaxL1Track_ = tauAssocCone_;
0227 int ntr = 0;
0228 std::vector<seedAtEC> VecSeedsatEC;
0229
0230 for (unsigned iS = 0; iS < pixelTrackRefs.size(); iS++) {
0231 bool vtxMatch = false;
0232
0233 reco::VertexCollection::const_iterator vitSel;
0234 double minDZ = 1000;
0235 bool found(false);
0236 for (reco::VertexCollection::const_iterator vit = pVert->begin(); vit != pVert->end(); vit++) {
0237 if (std::abs(pixelTrackRefs[iS]->dz(vit->position())) < minDZ) {
0238 minDZ = std::abs(pixelTrackRefs[iS]->dz(vit->position()));
0239 vitSel = vit;
0240 found = true;
0241 }
0242 }
0243
0244 if (found) {
0245 if (std::abs(pixelTrackRefs[iS]->dxy(vitSel->position())) < vtxCutSeed_)
0246 vtxMatch = true;
0247 } else {
0248 vtxMatch = true;
0249 }
0250 #ifdef EDM_ML_DEBUG
0251 edm::LogVerbatim("IsoTrack") << "minZD " << minDZ << " Found " << found << ":" << vtxMatch;
0252 #endif
0253
0254 double R = reco::deltaR(etaTriggered, phiTriggered, pixelTrackRefs[iS]->eta(), pixelTrackRefs[iS]->phi());
0255 #ifdef EDM_ML_DEBUG
0256 edm::LogVerbatim("IsoTrack") << "Distance to L1 " << R << ":" << tauUnbiasCone_ << " Result "
0257 << (R < tauUnbiasCone_);
0258 #endif
0259 if (R < tauUnbiasCone_)
0260 continue;
0261
0262
0263 bool tmatch = false;
0264 l1t::TauBxCollection::const_iterator selj;
0265 for (l1t::TauBxCollection::const_iterator tj = l1eTauJets->begin(); tj != l1eTauJets->end(); tj++) {
0266 if (reco::deltaR(pixelTrackRefs[iS]->momentum().eta(),
0267 pixelTrackRefs[iS]->momentum().phi(),
0268 tj->momentum().eta(),
0269 tj->momentum().phi()) > drMaxL1Track_)
0270 continue;
0271 selj = tj;
0272 tmatch = true;
0273 }
0274 #ifdef EDM_ML_DEBUG
0275 edm::LogVerbatim("IsoTrack") << "tMatch " << tmatch;
0276 #endif
0277
0278 std::pair<double, double> seedCooAtEC;
0279
0280 if (found)
0281 seedCooAtEC = GetEtaPhiAtEcal(pixelTrackRefs[iS]->eta(),
0282 pixelTrackRefs[iS]->phi(),
0283 pixelTrackRefs[iS]->pt(),
0284 pixelTrackRefs[iS]->charge(),
0285 vitSel->z());
0286
0287 else
0288 seedCooAtEC = GetEtaPhiAtEcal(pixelTrackRefs[iS]->eta(),
0289 pixelTrackRefs[iS]->phi(),
0290 pixelTrackRefs[iS]->pt(),
0291 pixelTrackRefs[iS]->charge(),
0292 0);
0293 seedAtEC seed(iS, (tmatch || vtxMatch), seedCooAtEC.first, seedCooAtEC.second);
0294 VecSeedsatEC.push_back(seed);
0295 #ifdef EDM_ML_DEBUG
0296 edm::LogVerbatim("IsoTrack") << "Seed " << seedCooAtEC.first << seedCooAtEC.second;
0297 #endif
0298 }
0299 for (unsigned int i = 0; i < VecSeedsatEC.size(); i++) {
0300 unsigned int iSeed = VecSeedsatEC[i].index;
0301 if (!VecSeedsatEC[i].ok)
0302 continue;
0303 if (pixelTrackRefs[iSeed]->p() < minPTrackValue_)
0304 continue;
0305 l1t::TauBxCollection::const_iterator selj;
0306 for (l1t::TauBxCollection::const_iterator tj = l1eTauJets->begin(); tj != l1eTauJets->end(); tj++) {
0307 if (reco::deltaR(pixelTrackRefs[iSeed]->momentum().eta(),
0308 pixelTrackRefs[iSeed]->momentum().phi(),
0309 tj->momentum().eta(),
0310 tj->momentum().phi()) > drMaxL1Track_)
0311 continue;
0312 selj = tj;
0313 }
0314 double maxP = 0;
0315 double sumP = 0;
0316 for (unsigned int j = 0; j < VecSeedsatEC.size(); j++) {
0317 if (i == j)
0318 continue;
0319 unsigned int iSurr = VecSeedsatEC[j].index;
0320
0321 if (reco::deltaR(pixelTrackRefs[iSeed]->eta(),
0322 pixelTrackRefs[iSeed]->phi(),
0323 pixelTrackRefs[iSurr]->eta(),
0324 pixelTrackRefs[iSurr]->phi()) > prelimCone_)
0325 continue;
0326 double minDZ2(1000);
0327 bool found(false);
0328 reco::VertexCollection::const_iterator vitSel2;
0329 for (reco::VertexCollection::const_iterator vit = pVert->begin(); vit != pVert->end(); vit++) {
0330 if (std::abs(pixelTrackRefs[iSurr]->dz(vit->position())) < minDZ2) {
0331 minDZ2 = std::abs(pixelTrackRefs[iSurr]->dz(vit->position()));
0332 vitSel2 = vit;
0333 found = true;
0334 }
0335 }
0336
0337 if (found && std::abs(pixelTrackRefs[iSurr]->dxy(vitSel2->position())) > vtxCutIsol_)
0338 continue;
0339
0340 if (getDistInCM(VecSeedsatEC[i].eta, VecSeedsatEC[i].phi, VecSeedsatEC[j].eta, VecSeedsatEC[j].phi) <
0341 pixelIsolationConeSizeAtEC_) {
0342 sumP += pixelTrackRefs[iSurr]->p();
0343 if (pixelTrackRefs[iSurr]->p() > maxP)
0344 maxP = pixelTrackRefs[iSurr]->p();
0345 }
0346 }
0347 if (maxP < maxPForIsolationValue_) {
0348 reco::IsolatedPixelTrackCandidate newCandidate(
0349 pixelTrackRefs[iSeed], l1t::TauRef(l1eTauJets, selj - l1eTauJets->begin()), maxP, sumP);
0350 newCandidate.setEtaPhiEcal(VecSeedsatEC[i].eta, VecSeedsatEC[i].phi);
0351 trackCollection->push_back(newCandidate);
0352 ntr++;
0353 }
0354 }
0355 #ifdef EDM_ML_DEBUG
0356 edm::LogVerbatim("IsoTrack") << "Number of Isolated Track " << ntr;
0357 #endif
0358
0359 theEvent.put(std::move(trackCollection));
0360 }
0361
0362 double IsolatedPixelTrackCandidateL1TProducer::getDistInCM(double eta1, double phi1, double eta2, double phi2) {
0363 double Rec;
0364 double theta1 = 2 * atan(exp(-eta1));
0365 double theta2 = 2 * atan(exp(-eta2));
0366 if (std::abs(eta1) < 1.479)
0367 Rec = rEB_;
0368 else if (std::abs(eta1) > 1.479 && std::abs(eta1) < 7.0)
0369 Rec = tan(theta1) * zEE_;
0370 else
0371 return 1000;
0372
0373
0374 double angle =
0375 acos((sin(theta1) * sin(theta2) * (sin(phi1) * sin(phi2) + cos(phi1) * cos(phi2)) + cos(theta1) * cos(theta2)));
0376 if (angle < M_PI_2)
0377 return std::abs((Rec / sin(theta1)) * tan(angle));
0378 else
0379 return 1000;
0380 }
0381
0382 std::pair<double, double> IsolatedPixelTrackCandidateL1TProducer::GetEtaPhiAtEcal(
0383 double etaIP, double phiIP, double pT, int charge, double vtxZ) {
0384 double deltaPhi = 0;
0385 double etaEC = 100;
0386 double phiEC = 100;
0387
0388 double Rcurv = 9999999;
0389 if (bfVal_ != 0)
0390 Rcurv = pT * 33.3 * 100 / (bfVal_ * 10);
0391
0392 double ecDist = zEE_;
0393 double ecRad = rEB_;
0394 double theta = 2 * atan(exp(-etaIP));
0395 double zNew = 0;
0396 if (theta > M_PI_2)
0397 theta = M_PI - theta;
0398 if (std::abs(etaIP) < ebEtaBoundary_) {
0399 if ((0.5 * ecRad / Rcurv) > 1) {
0400 etaEC = 10000;
0401 deltaPhi = 0;
0402 } else {
0403 deltaPhi = -charge * asin(0.5 * ecRad / Rcurv);
0404 double alpha1 = 2 * asin(0.5 * ecRad / Rcurv);
0405 double z = ecRad / tan(theta);
0406 if (etaIP > 0)
0407 zNew = z * (Rcurv * alpha1) / ecRad + vtxZ;
0408 else
0409 zNew = -z * (Rcurv * alpha1) / ecRad + vtxZ;
0410 double zAbs = std::abs(zNew);
0411 if (zAbs < ecDist) {
0412 etaEC = -log(tan(0.5 * atan(ecRad / zAbs)));
0413 deltaPhi = -charge * asin(0.5 * ecRad / Rcurv);
0414 }
0415 if (zAbs > ecDist) {
0416 zAbs = (std::abs(etaIP) / etaIP) * ecDist;
0417 double Zflight = std::abs(zAbs - vtxZ);
0418 double alpha = (Zflight * ecRad) / (z * Rcurv);
0419 double Rec = 2 * Rcurv * sin(alpha / 2);
0420 deltaPhi = -charge * alpha / 2;
0421 etaEC = -log(tan(0.5 * atan(Rec / ecDist)));
0422 }
0423 }
0424 } else {
0425 zNew = (std::abs(etaIP) / etaIP) * ecDist;
0426 double Zflight = std::abs(zNew - vtxZ);
0427 double Rvirt = std::abs(Zflight * tan(theta));
0428 double Rec = 2 * Rcurv * sin(Rvirt / (2 * Rcurv));
0429 deltaPhi = -(charge) * (Rvirt / (2 * Rcurv));
0430 etaEC = -log(tan(0.5 * atan(Rec / ecDist)));
0431 }
0432
0433 if (zNew < 0)
0434 etaEC = -etaEC;
0435 phiEC = phiIP + deltaPhi;
0436
0437 if (phiEC < -M_PI)
0438 phiEC += M_2_PI;
0439 if (phiEC > M_PI)
0440 phiEC -= M_2_PI;
0441
0442 std::pair<double, double> retVal(etaEC, phiEC);
0443 return retVal;
0444 }
0445 #include "FWCore/PluginManager/interface/ModuleDef.h"
0446 #include "FWCore/Framework/interface/MakerMacros.h"
0447
0448 DEFINE_FWK_MODULE(IsolatedPixelTrackCandidateL1TProducer);