File indexing completed on 2023-03-17 11:17:30
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronSeedGenerator.h"
0019 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronUtilities.h"
0020 #include "RecoTracker/Record/interface/TrackerRecoGeometryRecord.h"
0021 #include "RecoTracker/Record/interface/CkfComponentsRecord.h"
0022 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0023 #include "Geometry/Records/interface/TrackerTopologyRcd.h"
0024 #include "TrackingTools/Records/interface/TrackingComponentsRecord.h"
0025 #include "TrackingTools/MaterialEffects/interface/PropagatorWithMaterial.h"
0026 #include "RecoTracker/MeasurementDet/interface/MeasurementTrackerEvent.h"
0027 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0028
0029 #include <vector>
0030 #include <utility>
0031
0032 namespace {
0033
0034 bool equivalent(const TrajectorySeed &s1, const TrajectorySeed &s2) {
0035 if (s1.nHits() != s2.nHits())
0036 return false;
0037
0038 const TrajectorySeed::RecHitRange r1 = s1.recHits();
0039 const TrajectorySeed::RecHitRange r2 = s2.recHits();
0040 for (auto i1 = r1.begin(), i2 = r2.begin(); i1 != r1.end(); ++i1, ++i2) {
0041 if (!i1->isValid() || !i2->isValid())
0042 return false;
0043 if (i1->geographicalId() != i2->geographicalId())
0044 return false;
0045 if (!(i1->localPosition() == i2->localPosition()))
0046 return false;
0047 }
0048
0049 return true;
0050 }
0051
0052 void addSeed(reco::ElectronSeed &seed, const SeedWithInfo *info, bool positron, reco::ElectronSeedCollection &out) {
0053 if (!info) {
0054 out.emplace_back(seed);
0055 return;
0056 }
0057
0058 if (positron) {
0059 seed.setPosAttributes(info->dRz2, info->dPhi2, info->dRz1, info->dPhi1);
0060 } else {
0061 seed.setNegAttributes(info->dRz2, info->dPhi2, info->dRz1, info->dPhi1);
0062 }
0063 for (auto &res : out) {
0064 if ((seed.caloCluster().key() == res.caloCluster().key()) && (seed.hitsMask() == res.hitsMask()) &&
0065 equivalent(seed, res)) {
0066 if (positron) {
0067 if (res.dRZPos(1) == std::numeric_limits<float>::infinity() &&
0068 res.dRZNeg(1) != std::numeric_limits<float>::infinity()) {
0069 res.setPosAttributes(info->dRz2, info->dPhi2, info->dRz1, info->dPhi1);
0070 seed.setNegAttributes(res.dRZNeg(1), res.dPhiNeg(1), res.dRZNeg(0), res.dPhiNeg(0));
0071 break;
0072 } else {
0073 if (res.dRZPos(1) != std::numeric_limits<float>::infinity()) {
0074 if (res.dRZPos(1) != seed.dRZPos(1)) {
0075 edm::LogWarning("ElectronSeedGenerator|BadValue")
0076 << "this similar old seed already has another dRz2Pos"
0077 << "\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: " << (unsigned int)res.hitsMask() << "/"
0078 << res.dRZNeg(1) << "/" << res.dPhiNeg(1) << "/" << res.dRZPos(1) << "/" << res.dPhiPos(1)
0079 << "\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: " << (unsigned int)seed.hitsMask() << "/"
0080 << seed.dRZNeg(1) << "/" << seed.dPhiNeg(1) << "/" << seed.dRZPos(1) << "/" << seed.dPhiPos(1);
0081 }
0082 }
0083 }
0084 } else {
0085 if (res.dRZNeg(1) == std::numeric_limits<float>::infinity() &&
0086 res.dRZPos(1) != std::numeric_limits<float>::infinity()) {
0087 res.setNegAttributes(info->dRz2, info->dPhi2, info->dRz1, info->dPhi1);
0088 seed.setPosAttributes(res.dRZPos(1), res.dPhiPos(1), res.dRZPos(0), res.dPhiPos(0));
0089 break;
0090 } else {
0091 if (res.dRZNeg(1) != std::numeric_limits<float>::infinity()) {
0092 if (res.dRZNeg(1) != seed.dRZNeg(1)) {
0093 edm::LogWarning("ElectronSeedGenerator|BadValue")
0094 << "this old seed already has another dRz2"
0095 << "\nold seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: " << (unsigned int)res.hitsMask() << "/"
0096 << res.dRZNeg(1) << "/" << res.dPhiNeg(1) << "/" << res.dRZPos(1) << "/" << res.dPhiPos(1)
0097 << "\nnew seed mask/dRz2/dPhi2/dRz2Pos/dPhi2Pos: " << (unsigned int)seed.hitsMask() << "/"
0098 << seed.dRZNeg(1) << "/" << seed.dPhiNeg(1) << "/" << seed.dRZPos(1) << "/" << seed.dPhiPos(1);
0099 }
0100 }
0101 }
0102 }
0103 }
0104 }
0105
0106 out.emplace_back(seed);
0107 }
0108
0109 void seedsFromTrajectorySeeds(const std::vector<SeedWithInfo> &pixelSeeds,
0110 const reco::ElectronSeed::CaloClusterRef &cluster,
0111 reco::ElectronSeedCollection &out,
0112 bool positron) {
0113 if (!pixelSeeds.empty()) {
0114 LogDebug("ElectronSeedGenerator") << "Compatible " << (positron ? "positron" : "electron") << " seeds found.";
0115 }
0116
0117 std::vector<SeedWithInfo>::const_iterator s;
0118 for (s = pixelSeeds.begin(); s != pixelSeeds.end(); s++) {
0119 reco::ElectronSeed seed(s->seed);
0120 seed.setCaloCluster(cluster);
0121 seed.initTwoHitSeed(s->hitsMask);
0122 addSeed(seed, &*s, positron, out);
0123 }
0124 }
0125
0126 }
0127
0128 ElectronSeedGenerator::ElectronSeedGenerator(const edm::ParameterSet &pset,
0129 const ElectronSeedGenerator::Tokens &ts,
0130 edm::ConsumesCollector &&cc)
0131 : dynamicPhiRoad_(pset.getParameter<bool>("dynamicPhiRoad")),
0132 verticesTag_(ts.token_vtx),
0133 beamSpotTag_(ts.token_bs),
0134 magFieldToken_{cc.esConsumes()},
0135 trackerGeometryToken_{cc.esConsumes()},
0136 lowPtThresh_(pset.getParameter<double>("LowPtThreshold")),
0137 highPtThresh_(pset.getParameter<double>("HighPtThreshold")),
0138 nSigmasDeltaZ1_(pset.getParameter<double>("nSigmasDeltaZ1")),
0139 deltaZ1WithVertex_(pset.getParameter<double>("deltaZ1WithVertex")),
0140 sizeWindowENeg_(pset.getParameter<double>("SizeWindowENeg")),
0141 deltaPhi1Low_(pset.getParameter<double>("DeltaPhi1Low")),
0142 deltaPhi1High_(pset.getParameter<double>("DeltaPhi1High")),
0143
0144 dPhi1Coef2_(dynamicPhiRoad_ ? (deltaPhi1Low_ - deltaPhi1High_) / (1. / lowPtThresh_ - 1. / highPtThresh_) : 0.),
0145 dPhi1Coef1_(dynamicPhiRoad_ ? deltaPhi1Low_ - dPhi1Coef2_ / lowPtThresh_ : 0.),
0146
0147 useRecoVertex_(pset.getParameter<bool>("useRecoVertex")),
0148
0149 deltaPhi2B_(pset.getParameter<double>("DeltaPhi2B")),
0150 deltaPhi2F_(pset.getParameter<double>("DeltaPhi2F")),
0151 phiMin2B_(-pset.getParameter<double>("PhiMax2B")),
0152 phiMin2F_(-pset.getParameter<double>("PhiMax2F")),
0153 phiMax2B_(pset.getParameter<double>("PhiMax2B")),
0154 phiMax2F_(pset.getParameter<double>("PhiMax2F")),
0155 matcher_(pset.getParameter<double>("ePhiMin1"),
0156 pset.getParameter<double>("ePhiMax1"),
0157 phiMin2B_,
0158 phiMax2B_,
0159 phiMin2F_,
0160 phiMax2F_,
0161 pset.getParameter<double>("z2MaxB"),
0162 pset.getParameter<double>("r2MaxF"),
0163 pset.getParameter<double>("rMaxI"),
0164 useRecoVertex_) {}
0165
0166 void ElectronSeedGenerator::setupES(const edm::EventSetup &setup) {
0167 auto newMagField = magneticFieldWatcher_.check(setup);
0168 auto newTrackerGeom = trackerGeometryWatcher_.check(setup);
0169 if (newMagField || newTrackerGeom) {
0170 matcher_.setES(setup.getData(magFieldToken_), setup.getData(trackerGeometryToken_));
0171 }
0172 }
0173
0174 void ElectronSeedGenerator::run(edm::Event &e,
0175 const reco::SuperClusterRefVector &sclRefs,
0176 const std::vector<const TrajectorySeedCollection *> &seedsV,
0177 reco::ElectronSeedCollection &out) {
0178 initialSeedCollectionVector_ = &seedsV;
0179
0180
0181 auto const &beamSpot = e.get(beamSpotTag_);
0182
0183
0184 std::vector<reco::Vertex> const *vertices = nullptr;
0185 if (useRecoVertex_)
0186 vertices = &e.get(verticesTag_);
0187
0188 for (unsigned int i = 0; i < sclRefs.size(); ++i) {
0189
0190 LogDebug("ElectronSeedGenerator") << "new cluster, calling seedsFromThisCluster";
0191 seedsFromThisCluster(sclRefs[i], beamSpot, vertices, out);
0192 }
0193
0194 LogDebug("ElectronSeedGenerator") << ": For event " << e.id();
0195 LogDebug("ElectronSeedGenerator") << "Nr of superclusters after filter: " << sclRefs.size()
0196 << ", no. of ElectronSeeds found = " << out.size();
0197 }
0198
0199 void ElectronSeedGenerator::seedsFromThisCluster(edm::Ref<reco::SuperClusterCollection> seedCluster,
0200 reco::BeamSpot const &beamSpot,
0201 std::vector<reco::Vertex> const *vertices,
0202 reco::ElectronSeedCollection &out) {
0203 float clusterEnergy = seedCluster->energy();
0204 GlobalPoint clusterPos(seedCluster->position().x(), seedCluster->position().y(), seedCluster->position().z());
0205 reco::ElectronSeed::CaloClusterRef caloCluster(seedCluster);
0206
0207 if (dynamicPhiRoad_) {
0208 float clusterEnergyT = clusterEnergy / cosh(EleRelPoint(clusterPos, beamSpot.position()).eta());
0209
0210 float deltaPhi1;
0211 if (clusterEnergyT < lowPtThresh_) {
0212 deltaPhi1 = deltaPhi1Low_;
0213 } else if (clusterEnergyT > highPtThresh_) {
0214 deltaPhi1 = deltaPhi1High_;
0215 } else {
0216 deltaPhi1 = dPhi1Coef1_ + dPhi1Coef2_ / clusterEnergyT;
0217 }
0218
0219 matcher_.set1stLayer(-deltaPhi1 * sizeWindowENeg_, deltaPhi1 * (1. - sizeWindowENeg_));
0220 matcher_.set2ndLayer(-deltaPhi2B_ / 2., deltaPhi2B_ / 2., -deltaPhi2F_ / 2., deltaPhi2F_ / 2.);
0221 }
0222
0223 if (!useRecoVertex_)
0224 {
0225 double sigmaZ = beamSpot.sigmaZ();
0226 double sigmaZ0Error = beamSpot.sigmaZ0Error();
0227 double sq = sqrt(sigmaZ * sigmaZ + sigmaZ0Error * sigmaZ0Error);
0228 double myZmin1 = beamSpot.position().z() - nSigmasDeltaZ1_ * sq;
0229 double myZmax1 = beamSpot.position().z() + nSigmasDeltaZ1_ * sq;
0230
0231 GlobalPoint vertexPos;
0232 ele_convert(beamSpot.position(), vertexPos);
0233
0234 matcher_.set1stLayerZRange(myZmin1, myZmax1);
0235
0236
0237 auto elePixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, -1.);
0238 seedsFromTrajectorySeeds(elePixelSeeds, caloCluster, out, false);
0239
0240 auto posPixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, 1.);
0241 seedsFromTrajectorySeeds(posPixelSeeds, caloCluster, out, true);
0242
0243 } else if (vertices)
0244 {
0245 for (auto const &vertex : *vertices) {
0246 GlobalPoint vertexPos(vertex.position().x(), vertex.position().y(), vertex.position().z());
0247 double myZmin1, myZmax1;
0248 if (vertexPos.z() == beamSpot.position().z()) {
0249 double sigmaZ = beamSpot.sigmaZ();
0250 double sigmaZ0Error = beamSpot.sigmaZ0Error();
0251 double sq = sqrt(sigmaZ * sigmaZ + sigmaZ0Error * sigmaZ0Error);
0252 myZmin1 = beamSpot.position().z() - nSigmasDeltaZ1_ * sq;
0253 myZmax1 = beamSpot.position().z() + nSigmasDeltaZ1_ * sq;
0254 } else {
0255 myZmin1 = vertex.position().z() - deltaZ1WithVertex_;
0256 myZmax1 = vertex.position().z() + deltaZ1WithVertex_;
0257 }
0258
0259 matcher_.set1stLayerZRange(myZmin1, myZmax1);
0260
0261
0262 auto elePixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, -1.);
0263 seedsFromTrajectorySeeds(elePixelSeeds, caloCluster, out, false);
0264
0265 auto posPixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, 1.);
0266 seedsFromTrajectorySeeds(posPixelSeeds, caloCluster, out, true);
0267 }
0268 }
0269 }