Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-02-14 13:34:02

0001 // -*- C++ -*-
0002 //
0003 // Package:    EgammaElectronAlgos
0004 // Class:      ElectronSeedGenerator.
0005 //
0006 /**\class ElectronSeedGenerator EgammaElectronAlgos/ElectronSeedGenerator
0007 
0008  Description: Top algorithm producing ElectronSeeds, ported from ORCA
0009 
0010  Implementation:
0011      future redesign...
0012 */
0013 //
0014 // Original Author:  Ursula Berthon, Claude Charlot
0015 //         Created:  Mon Mar 27 13:22:06 CEST 2006
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 }  // namespace
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       // so that deltaPhi1 = dPhi1Coef1_ + dPhi1Coef2_/clusterEnergyT
0144       dPhi1Coef2_(dynamicPhiRoad_ ? (deltaPhi1Low_ - deltaPhi1High_) / (1. / lowPtThresh_ - 1. / highPtThresh_) : 0.),
0145       dPhi1Coef1_(dynamicPhiRoad_ ? deltaPhi1Low_ - dPhi1Coef2_ / lowPtThresh_ : 0.),
0146       // use of reco vertex
0147       useRecoVertex_(pset.getParameter<bool>("useRecoVertex")),
0148       // new B/F configurables
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   if (magneticFieldWatcher_.check(setup) || trackerGeometryWatcher_.check(setup)) {
0168     matcher_.setES(setup.getData(magFieldToken_), setup.getData(trackerGeometryToken_));
0169   }
0170 }
0171 
0172 void ElectronSeedGenerator::run(edm::Event &e,
0173                                 const reco::SuperClusterRefVector &sclRefs,
0174                                 const std::vector<const TrajectorySeedCollection *> &seedsV,
0175                                 reco::ElectronSeedCollection &out) {
0176   initialSeedCollectionVector_ = &seedsV;
0177 
0178   // get the beamspot from the Event:
0179   auto const &beamSpot = e.get(beamSpotTag_);
0180 
0181   // if required get the vertices
0182   std::vector<reco::Vertex> const *vertices = nullptr;
0183   if (useRecoVertex_)
0184     vertices = &e.get(verticesTag_);
0185 
0186   for (unsigned int i = 0; i < sclRefs.size(); ++i) {
0187     // Find the seeds
0188     LogDebug("ElectronSeedGenerator") << "new cluster, calling seedsFromThisCluster";
0189     seedsFromThisCluster(sclRefs[i], beamSpot, vertices, out);
0190   }
0191 
0192   LogDebug("ElectronSeedGenerator") << ": For event " << e.id();
0193   LogDebug("ElectronSeedGenerator") << "Nr of superclusters after filter: " << sclRefs.size()
0194                                     << ", no. of ElectronSeeds found  = " << out.size();
0195 }
0196 
0197 void ElectronSeedGenerator::seedsFromThisCluster(edm::Ref<reco::SuperClusterCollection> seedCluster,
0198                                                  reco::BeamSpot const &beamSpot,
0199                                                  std::vector<reco::Vertex> const *vertices,
0200                                                  reco::ElectronSeedCollection &out) {
0201   float clusterEnergy = seedCluster->energy();
0202   GlobalPoint clusterPos(seedCluster->position().x(), seedCluster->position().y(), seedCluster->position().z());
0203   reco::ElectronSeed::CaloClusterRef caloCluster(seedCluster);
0204 
0205   if (dynamicPhiRoad_) {
0206     float clusterEnergyT = clusterEnergy / cosh(EleRelPoint(clusterPos, beamSpot.position()).eta());
0207 
0208     float deltaPhi1;
0209     if (clusterEnergyT < lowPtThresh_) {
0210       deltaPhi1 = deltaPhi1Low_;
0211     } else if (clusterEnergyT > highPtThresh_) {
0212       deltaPhi1 = deltaPhi1High_;
0213     } else {
0214       deltaPhi1 = dPhi1Coef1_ + dPhi1Coef2_ / clusterEnergyT;
0215     }
0216 
0217     matcher_.set1stLayer(-deltaPhi1 * sizeWindowENeg_, deltaPhi1 * (1. - sizeWindowENeg_));
0218     matcher_.set2ndLayer(-deltaPhi2B_ / 2., deltaPhi2B_ / 2., -deltaPhi2F_ / 2., deltaPhi2F_ / 2.);
0219   }
0220 
0221   if (!useRecoVertex_)  // here use the beam spot position
0222   {
0223     double sigmaZ = beamSpot.sigmaZ();
0224     double sigmaZ0Error = beamSpot.sigmaZ0Error();
0225     double sq = sqrt(sigmaZ * sigmaZ + sigmaZ0Error * sigmaZ0Error);
0226     double myZmin1 = beamSpot.position().z() - nSigmasDeltaZ1_ * sq;
0227     double myZmax1 = beamSpot.position().z() + nSigmasDeltaZ1_ * sq;
0228 
0229     GlobalPoint vertexPos;
0230     ele_convert(beamSpot.position(), vertexPos);
0231 
0232     matcher_.set1stLayerZRange(myZmin1, myZmax1);
0233 
0234     // try electron
0235     auto elePixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, -1.);
0236     seedsFromTrajectorySeeds(elePixelSeeds, caloCluster, out, false);
0237     // try positron
0238     auto posPixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, 1.);
0239     seedsFromTrajectorySeeds(posPixelSeeds, caloCluster, out, true);
0240 
0241   } else if (vertices)  // here we use the reco vertices
0242   {
0243     for (auto const &vertex : *vertices) {
0244       GlobalPoint vertexPos(vertex.position().x(), vertex.position().y(), vertex.position().z());
0245       double myZmin1, myZmax1;
0246       if (vertexPos.z() == beamSpot.position().z()) {  // in case vetex not found
0247         double sigmaZ = beamSpot.sigmaZ();
0248         double sigmaZ0Error = beamSpot.sigmaZ0Error();
0249         double sq = sqrt(sigmaZ * sigmaZ + sigmaZ0Error * sigmaZ0Error);
0250         myZmin1 = beamSpot.position().z() - nSigmasDeltaZ1_ * sq;
0251         myZmax1 = beamSpot.position().z() + nSigmasDeltaZ1_ * sq;
0252       } else {  // a vertex has been recoed
0253         myZmin1 = vertex.position().z() - deltaZ1WithVertex_;
0254         myZmax1 = vertex.position().z() + deltaZ1WithVertex_;
0255       }
0256 
0257       matcher_.set1stLayerZRange(myZmin1, myZmax1);
0258 
0259       // try electron
0260       auto elePixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, -1.);
0261       seedsFromTrajectorySeeds(elePixelSeeds, caloCluster, out, false);
0262       // try positron
0263       auto posPixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, 1.);
0264       seedsFromTrajectorySeeds(posPixelSeeds, caloCluster, out, true);
0265     }
0266   }
0267 }