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