CMS 3D CMS Logo

All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Properties Friends Macros Pages
EgammaHLTGsfTrackVarProducer.cc
Go to the documentation of this file.
1 
9 //this class is designed to calculate dEtaIn,dPhiIn gsf track - supercluster pairs
10 //it can take as input std::vector<Electron> which the gsf track-sc is already done
11 //or it can run over the std::vector<GsfTrack> directly in which case it will pick the smallest dEta,dPhi
12 //the dEta, dPhi do not have to be from the same track
13 //it can optionally set dEta, dPhi to 0 based on the number of tracks found
14 
17 
21 
25 
27 
35 
42 
49 
52 
55 
57 
58 #include <memory>
59 
61 private:
63  unsigned long long cacheIDTDGeom_;
64  unsigned long long cacheIDMagField_;
65 
68 
70 
71  public:
72  TrackExtrapolator() : cacheIDTDGeom_(0), cacheIDMagField_(0), mtsTransform_(nullptr) {}
76 
77  void setup(const edm::EventSetup& iSetup);
78 
79  GlobalPoint extrapolateTrackPosToPoint(const reco::GsfTrack& gsfTrack, const GlobalPoint& pointToExtrapTo);
80  GlobalVector extrapolateTrackMomToPoint(const reco::GsfTrack& gsfTrack, const GlobalPoint& pointToExtrapTo);
81 
84  };
85 
86 public:
89  void produce(edm::Event&, const edm::EventSetup&) override;
90  static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
91 
92 private:
97 
103 };
104 
107  consumes<reco::RecoEcalCandidateCollection>(config.getParameter<edm::InputTag>("recoEcalCandidateProducer"))),
108  inputCollectionTag1_(consumes<reco::ElectronCollection>(config.getParameter<edm::InputTag>("inputCollection"))),
109  inputCollectionTag2_(consumes<reco::GsfTrackCollection>(config.getParameter<edm::InputTag>("inputCollection"))),
110  beamSpotTag_(consumes<reco::BeamSpot>(config.getParameter<edm::InputTag>("beamSpotProducer"))),
111  upperTrackNrToRemoveCut_(config.getParameter<int>("upperTrackNrToRemoveCut")),
112  lowerTrackNrToRemoveCut_(config.getParameter<int>("lowerTrackNrToRemoveCut")),
113  useDefaultValuesForBarrel_(config.getParameter<bool>("useDefaultValuesForBarrel")),
114  useDefaultValuesForEndcap_(config.getParameter<bool>("useDefaultValuesForEndcap")) {
115  //register your products
116  produces<reco::RecoEcalCandidateIsolationMap>("Deta").setBranchAlias("deta");
117  produces<reco::RecoEcalCandidateIsolationMap>("DetaSeed").setBranchAlias("detaseed");
118  produces<reco::RecoEcalCandidateIsolationMap>("Dphi").setBranchAlias("dphi");
119  produces<reco::RecoEcalCandidateIsolationMap>("OneOESuperMinusOneOP");
120  produces<reco::RecoEcalCandidateIsolationMap>("OneOESeedMinusOneOP");
121  produces<reco::RecoEcalCandidateIsolationMap>("MissingHits").setBranchAlias("missinghits");
122  produces<reco::RecoEcalCandidateIsolationMap>("Chi2").setBranchAlias("chi2");
123  produces<reco::RecoEcalCandidateIsolationMap>("ValidHits").setBranchAlias("validhits");
124 }
125 
127 
130  desc.add<edm::InputTag>(("recoEcalCandidateProducer"), edm::InputTag("hltRecoEcalSuperClusterActivityCandidate"));
131  desc.add<edm::InputTag>(("inputCollection"), edm::InputTag("hltActivityElectronGsfTracks"));
132  desc.add<edm::InputTag>(("beamSpotProducer"), edm::InputTag("hltOnlineBeamSpot"));
133  desc.add<int>(("upperTrackNrToRemoveCut"), 9999);
134  desc.add<int>(("lowerTrackNrToRemoveCut"), -1);
135  desc.add<bool>(("useDefaultValuesForBarrel"), false);
136  desc.add<bool>(("useDefaultValuesForEndcap"), false);
137 
138  descriptions.add("hltEgammaHLTGsfTrackVarProducer", desc);
139 }
141  trackExtrapolator_.setup(iSetup);
142 
143  // Get the HLT filtered objects
145  iEvent.getByToken(recoEcalCandTag_, recoEcalCandHandle);
146 
148  iEvent.getByToken(inputCollectionTag1_, electronHandle);
149 
151  if (!electronHandle.isValid())
152  iEvent.getByToken(inputCollectionTag2_, gsfTracksHandle);
153 
154  edm::Handle<reco::BeamSpot> recoBeamSpotHandle;
155  iEvent.getByToken(beamSpotTag_, recoBeamSpotHandle);
156  // gets its position
157  const reco::BeamSpot& beamSpot = *recoBeamSpotHandle;
158 
159  edm::ESHandle<MagneticField> theMagField;
160  iSetup.get<IdealMagneticFieldRecord>().get(theMagField);
161 
162  reco::RecoEcalCandidateIsolationMap dEtaMap(recoEcalCandHandle);
163  reco::RecoEcalCandidateIsolationMap dEtaSeedMap(recoEcalCandHandle);
164  reco::RecoEcalCandidateIsolationMap dPhiMap(recoEcalCandHandle);
165  reco::RecoEcalCandidateIsolationMap oneOverESuperMinusOneOverPMap(recoEcalCandHandle);
166  reco::RecoEcalCandidateIsolationMap oneOverESeedMinusOneOverPMap(recoEcalCandHandle);
167  reco::RecoEcalCandidateIsolationMap missingHitsMap(recoEcalCandHandle);
168  reco::RecoEcalCandidateIsolationMap validHitsMap(recoEcalCandHandle);
169  reco::RecoEcalCandidateIsolationMap chi2Map(recoEcalCandHandle);
170 
171  for (reco::RecoEcalCandidateCollection::const_iterator iRecoEcalCand = recoEcalCandHandle->begin();
172  iRecoEcalCand != recoEcalCandHandle->end();
173  iRecoEcalCand++) {
174  reco::RecoEcalCandidateRef recoEcalCandRef(recoEcalCandHandle, iRecoEcalCand - recoEcalCandHandle->begin());
175 
176  const reco::SuperClusterRef scRef = recoEcalCandRef->superCluster();
177  bool isBarrel = std::abs(recoEcalCandRef->eta()) < 1.479;
178  //the idea is that we can take the tracks from properly associated electrons or just take all gsf tracks with that sc as a seed
179  std::vector<const reco::GsfTrack*> gsfTracks;
180  if (electronHandle.isValid()) {
181  for (reco::ElectronCollection::const_iterator eleIt = electronHandle->begin(); eleIt != electronHandle->end();
182  eleIt++) {
183  if (eleIt->superCluster() == scRef) {
184  gsfTracks.push_back(&*eleIt->gsfTrack());
185  }
186  }
187  } else {
188  for (reco::GsfTrackCollection::const_iterator trkIt = gsfTracksHandle->begin(); trkIt != gsfTracksHandle->end();
189  ++trkIt) {
190  edm::RefToBase<TrajectorySeed> seed = trkIt->extra()->seedRef();
192  edm::RefToBase<reco::CaloCluster> caloCluster = elseed->caloCluster();
193  reco::SuperClusterRef scRefFromTrk = caloCluster.castTo<reco::SuperClusterRef>();
194  if (scRefFromTrk == scRef) {
195  gsfTracks.push_back(&*trkIt);
196  }
197  }
198  }
199 
200  int validHitsValue = 0;
201  float chi2Value = 9999999.;
202  float missingHitsValue = 9999999;
203  float dEtaInValue = 999999;
204  float dEtaSeedInValue = 999999;
205  float dPhiInValue = 999999;
206  float oneOverESuperMinusOneOverPValue = 999999;
207  float oneOverESeedMinusOneOverPValue = 999999;
208 
209  const int nrTracks = gsfTracks.size();
210  const bool rmCutsDueToNrTracks = nrTracks <= lowerTrackNrToRemoveCut_ || nrTracks >= upperTrackNrToRemoveCut_;
211  //to use the default values, we require at least one track...
212  const bool useDefaultValues =
213  isBarrel ? useDefaultValuesForBarrel_ && nrTracks >= 1 : useDefaultValuesForEndcap_ && nrTracks >= 1;
214 
215  if (rmCutsDueToNrTracks || useDefaultValues) {
216  dEtaInValue = 0;
217  dEtaSeedInValue = 0;
218  dPhiInValue = 0;
219  missingHitsValue = 0;
220  validHitsValue = 100;
221  chi2Value = 0;
222  oneOverESuperMinusOneOverPValue = 0;
223  oneOverESeedMinusOneOverPValue = 0;
224  } else {
225  for (size_t trkNr = 0; trkNr < gsfTracks.size(); trkNr++) {
226  GlobalPoint scPos(scRef->x(), scRef->y(), scRef->z());
227  GlobalPoint trackExtrapToSC = trackExtrapolator_.extrapolateTrackPosToPoint(*gsfTracks[trkNr], scPos);
228  EleRelPointPair scAtVtx(scRef->position(), trackExtrapToSC, beamSpot.position());
229 
230  float trkP = gsfTracks[trkNr]->p();
231  if (scRef->energy() != 0 && trkP != 0) {
232  if (fabs(1 / scRef->energy() - 1 / trkP) < oneOverESuperMinusOneOverPValue)
233  oneOverESuperMinusOneOverPValue = fabs(1 / scRef->energy() - 1 / trkP);
234  }
235  if (scRef->seed().isNonnull() && scRef->seed()->energy() != 0 && trkP != 0) {
236  if (fabs(1 / scRef->seed()->energy() - 1 / trkP) < oneOverESeedMinusOneOverPValue)
237  oneOverESeedMinusOneOverPValue = fabs(1 / scRef->seed()->energy() - 1 / trkP);
238  }
239 
240  if (gsfTracks[trkNr]->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS) < missingHitsValue)
241  missingHitsValue = gsfTracks[trkNr]->hitPattern().numberOfLostHits(reco::HitPattern::MISSING_INNER_HITS);
242 
243  if (gsfTracks[trkNr]->numberOfValidHits() < validHitsValue)
244  validHitsValue = gsfTracks[trkNr]->numberOfValidHits();
245 
246  if (gsfTracks[trkNr]->numberOfValidHits() < chi2Value)
247  chi2Value = gsfTracks[trkNr]->normalizedChi2();
248 
249  if (fabs(scAtVtx.dEta()) < dEtaInValue)
250  dEtaInValue = fabs(scAtVtx.dEta()); //we are allowing them to come from different tracks
251 
252  if (fabs(scAtVtx.dEta()) < dEtaSeedInValue)
253  dEtaSeedInValue = fabs(scAtVtx.dEta() - scRef->position().eta() + scRef->seed()->position().eta());
254 
255  if (fabs(scAtVtx.dPhi()) < dPhiInValue)
256  dPhiInValue = fabs(scAtVtx.dPhi()); //we are allowing them to come from different tracks
257  }
258  }
259 
260  dEtaMap.insert(recoEcalCandRef, dEtaInValue);
261  dEtaSeedMap.insert(recoEcalCandRef, dEtaSeedInValue);
262  dPhiMap.insert(recoEcalCandRef, dPhiInValue);
263  oneOverESuperMinusOneOverPMap.insert(recoEcalCandRef, oneOverESuperMinusOneOverPValue);
264  oneOverESeedMinusOneOverPMap.insert(recoEcalCandRef, oneOverESeedMinusOneOverPValue);
265  missingHitsMap.insert(recoEcalCandRef, missingHitsValue);
266  validHitsMap.insert(recoEcalCandRef, validHitsValue);
267  chi2Map.insert(recoEcalCandRef, chi2Value);
268  }
269 
270  auto dEtaMapForEvent = std::make_unique<reco::RecoEcalCandidateIsolationMap>(dEtaMap);
271  auto dEtaSeedMapForEvent = std::make_unique<reco::RecoEcalCandidateIsolationMap>(dEtaSeedMap);
272  auto dPhiMapForEvent = std::make_unique<reco::RecoEcalCandidateIsolationMap>(dPhiMap);
273  auto oneOverESuperMinusOneOverPMapForEvent =
274  std::make_unique<reco::RecoEcalCandidateIsolationMap>(oneOverESuperMinusOneOverPMap);
275  auto oneOverESeedMinusOneOverPMapForEvent =
276  std::make_unique<reco::RecoEcalCandidateIsolationMap>(oneOverESeedMinusOneOverPMap);
277  auto missingHitsForEvent = std::make_unique<reco::RecoEcalCandidateIsolationMap>(missingHitsMap);
278  auto validHitsForEvent = std::make_unique<reco::RecoEcalCandidateIsolationMap>(validHitsMap);
279  auto chi2ForEvent = std::make_unique<reco::RecoEcalCandidateIsolationMap>(chi2Map);
280 
281  iEvent.put(std::move(dEtaMapForEvent), "Deta");
282  iEvent.put(std::move(dEtaSeedMapForEvent), "DetaSeed");
283  iEvent.put(std::move(dPhiMapForEvent), "Dphi");
284  iEvent.put(std::move(oneOverESuperMinusOneOverPMapForEvent), "OneOESuperMinusOneOP");
285  iEvent.put(std::move(oneOverESeedMinusOneOverPMapForEvent), "OneOESeedMinusOneOP");
286  iEvent.put(std::move(missingHitsForEvent), "MissingHits");
287  iEvent.put(std::move(validHitsForEvent), "ValidHits");
288  iEvent.put(std::move(chi2ForEvent), "Chi2");
289 }
290 
293  : cacheIDTDGeom_(rhs.cacheIDTDGeom_),
294  cacheIDMagField_(rhs.cacheIDMagField_),
295  magField_(rhs.magField_),
296  trackerHandle_(rhs.trackerHandle_) {
297  if (rhs.mtsTransform_)
299  else
300  mtsTransform_ = nullptr;
301 }
302 
305  if (this != &rhs) { //just to ensure we're not copying ourselves
308  magField_ = rhs.magField_;
310 
311  delete mtsTransform_;
312  if (rhs.mtsTransform_)
314  else
315  mtsTransform_ = nullptr;
316  }
317  return this;
318 }
319 
321  bool updateField(false);
323  updateField = true;
324  cacheIDMagField_ = iSetup.get<IdealMagneticFieldRecord>().cacheIdentifier();
325  iSetup.get<IdealMagneticFieldRecord>().get(magField_);
326  }
327 
328  bool updateGeometry(false);
330  updateGeometry = true;
331  cacheIDTDGeom_ = iSetup.get<TrackerDigiGeometryRecord>().cacheIdentifier();
333  }
334 
335  if (updateField || updateGeometry || !mtsTransform_) {
336  delete mtsTransform_;
338  }
339 }
340 
342  const reco::GsfTrack& gsfTrack, const GlobalPoint& pointToExtrapTo) {
344  TrajectoryStateOnSurface posTSOS = mtsTransform()->extrapolatedState(innTSOS, pointToExtrapTo);
345  GlobalPoint extrapolatedPos;
346  multiTrajectoryStateMode::positionFromModeCartesian(posTSOS, extrapolatedPos);
347  return extrapolatedPos;
348 }
349 
351  const reco::GsfTrack& gsfTrack, const GlobalPoint& pointToExtrapTo) {
353  TrajectoryStateOnSurface posTSOS = mtsTransform()->extrapolatedState(innTSOS, pointToExtrapTo);
354  GlobalVector extrapolatedMom;
355  multiTrajectoryStateMode::momentumFromModeCartesian(posTSOS, extrapolatedMom);
356  return extrapolatedMom;
357 }
358 
unsigned long long cacheIdentifier() const
OrphanHandle< PROD > put(std::unique_ptr< PROD > product)
Put a new product.
Definition: Event.h:131
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:525
#define nullptr
Definition: config.py:1
const MultiTrajectoryStateTransform * mtsTransform() const
const edm::EDGetTokenT< reco::RecoEcalCandidateCollection > recoEcalCandTag_
GlobalPoint extrapolateTrackPosToPoint(const reco::GsfTrack &gsfTrack, const GlobalPoint &pointToExtrapTo)
const edm::EDGetTokenT< reco::GsfTrackCollection > inputCollectionTag2_
GlobalVector extrapolateTrackMomToPoint(const reco::GsfTrack &gsfTrack, const GlobalPoint &pointToExtrapTo)
int iEvent
Definition: GenABIO.cc:224
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
TrajectoryStateOnSurface extrapolatedState(const TrajectoryStateOnSurface tsos, const GlobalPoint &point) const
bool positionFromModeCartesian(TrajectoryStateOnSurface const &tsos, GlobalPoint &position)
const edm::EDGetTokenT< reco::BeamSpot > beamSpotTag_
std::vector< GsfTrack > GsfTrackCollection
collection of GsfTracks
Definition: GsfTrackFwd.h:9
edm::ESHandle< TrackerGeometry > trackerGeomHandle() const
Abs< T >::type abs(const T &t)
Definition: Abs.h:22
EgammaHLTGsfTrackVarProducer(const edm::ParameterSet &)
ParameterDescriptionBase * add(U const &iLabel, T const &value)
bool isValid() const
Definition: HandleBase.h:70
std::vector< Electron > ElectronCollection
collectin of Electron objects
Definition: ElectronFwd.h:9
void produce(edm::Event &, const edm::EventSetup &) override
const edm::EDGetTokenT< reco::ElectronCollection > inputCollectionTag1_
TrackExtrapolator * operator=(const TrackExtrapolator &rhs)
void insert(const key_type &k, const data_type &v)
insert an association
REF castTo() const
Definition: RefToBase.h:257
static void fillDescriptions(edm::ConfigurationDescriptions &descriptions)
void add(std::string const &label, ParameterSetDescription const &psetDescription)
TrajectoryStateOnSurface innerStateOnSurface(const reco::GsfTrack &tk) const
std::vector< RecoEcalCandidate > RecoEcalCandidateCollection
collectin of RecoEcalCandidate objects
fixed size matrix
HLT enums.
T get() const
Definition: EventSetup.h:73
const Point & position() const
position
Definition: BeamSpot.h:59
bool momentumFromModeCartesian(TrajectoryStateOnSurface const &tsos, GlobalVector &momentum)
T const * product() const
Definition: ESHandle.h:86
def move(src, dest)
Definition: eostools.py:511