Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:17:30 

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   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   // get the beamspot from the Event:
0181   auto const &beamSpot = e.get(beamSpotTag_);
0182 
0183   // if required get the vertices
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     // Find the seeds
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_)  // here use the beam spot position
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     // try electron
0237     auto elePixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, -1.);
0238     seedsFromTrajectorySeeds(elePixelSeeds, caloCluster, out, false);
0239     // try positron
0240     auto posPixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, 1.);
0241     seedsFromTrajectorySeeds(posPixelSeeds, caloCluster, out, true);
0242 
0243   } else if (vertices)  // here we use the reco 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()) {  // in case vetex not found
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 {  // a vertex has been recoed
0255         myZmin1 = vertex.position().z() - deltaZ1WithVertex_;
0256         myZmax1 = vertex.position().z() + deltaZ1WithVertex_;
0257       }
0258 
0259       matcher_.set1stLayerZRange(myZmin1, myZmax1);
0260 
0261       // try electron
0262       auto elePixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, -1.);
0263       seedsFromTrajectorySeeds(elePixelSeeds, caloCluster, out, false);
0264       // try positron
0265       auto posPixelSeeds = matcher_(*initialSeedCollectionVector_, clusterPos, vertexPos, clusterEnergy, 1.);
0266       seedsFromTrajectorySeeds(posPixelSeeds, caloCluster, out, true);
0267     }
0268   }
0269 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!