CMS 3D CMS Logo

ConvertedPhotonProducer.cc
Go to the documentation of this file.
1 
44 
45 #include <vector>
46 
48 public:
50 
51  void produce(edm::Event& evt, const edm::EventSetup& es) override;
52 
53 private:
54  void buildCollections(
55  edm::EventSetup const& es,
58  ElectronHcalHelper const& hcalHelper,
59  const edm::Handle<reco::TrackCollection>& trkHandle,
60  std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors>& allPairs,
61  reco::ConversionCollection& outputConvPhotonCollection);
64  reco::ConversionCollection& outputCollection);
65 
67  reco::CaloClusterPtr const& sc);
68 
69  float calculateMinApproachDistance(const reco::TrackRef& track1, const reco::TrackRef& track2);
70  void getCircleCenter(const reco::TrackRef& tk, double r, double& x0, double& y0);
71 
74 
77 
79 
82  HcalPFCuts const* hcalCuts_ = nullptr;
83 
84  // Register the product
87 
93 
96 
100 
104 
106  double maxHOverE_;
107  double minSCEt_;
110  double deltaCotCut_;
114 
115  std::unique_ptr<ElectronHcalHelper> hcalHelper_;
116 
119 
121 
123 };
124 
127 
129  : conversionOITrackProducer_{consumes(config.getParameter<std::string>("conversionOITrackProducer"))},
130  conversionIOTrackProducer_{consumes(config.getParameter<std::string>("conversionIOTrackProducer"))},
131  outInTrackSCAssociationCollection_{consumes({config.getParameter<std::string>("conversionOITrackProducer"),
132  config.getParameter<std::string>("outInTrackSCAssociation")})},
133  inOutTrackSCAssociationCollection_{consumes({config.getParameter<std::string>("conversionIOTrackProducer"),
134  config.getParameter<std::string>("inOutTrackSCAssociation")})},
135 
136  generalTrackProducer_{consumes(config.getParameter<edm::InputTag>("generalTracksSrc"))},
137  convertedPhotonCollectionPutToken_{
138  produces<reco::ConversionCollection>(config.getParameter<std::string>("convertedPhotonCollection"))},
139  cleanedConvertedPhotonCollectionPutToken_{
140  produces<reco::ConversionCollection>(config.getParameter<std::string>("cleanedConvertedPhotonCollection"))},
141 
142  bcBarrelCollection_{consumes(config.getParameter<edm::InputTag>("bcBarrelCollection"))},
143  bcEndcapCollection_{consumes(config.getParameter<edm::InputTag>("bcEndcapCollection"))},
144  scHybridBarrelProducer_{consumes(config.getParameter<edm::InputTag>("scHybridBarrelProducer"))},
145  scIslandEndcapProducer_{consumes(config.getParameter<edm::InputTag>("scIslandEndcapProducer"))},
146  hbheRecHits_{consumes(config.getParameter<edm::InputTag>("hbheRecHits"))},
147  caloGeomToken_{esConsumes()},
148  mFToken_{esConsumes()},
149  transientTrackToken_{esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))},
150  vertexFinder_{config},
151  algoName_{config.getParameter<std::string>("AlgorithmName")},
152 
153  hOverEConeSize_{config.getParameter<double>("hOverEConeSize")},
154  maxHOverE_{config.getParameter<double>("maxHOverE")},
155  minSCEt_{config.getParameter<double>("minSCEt")},
156  recoverOneTrackCase_{config.getParameter<bool>("recoverOneTrackCase")},
157  dRForConversionRecovery_{config.getParameter<double>("dRForConversionRecovery")},
158  deltaCotCut_{config.getParameter<double>("deltaCotCut")},
159  minApproachDisCut_{config.getParameter<double>("minApproachDisCut")},
160 
161  maxNumOfCandidates_{config.getParameter<int>("maxNumOfCandidates")},
162  risolveAmbiguity_{config.getParameter<bool>("risolveConversionAmbiguity")},
163  likelihoodWeights_{config.getParameter<std::string>("MVA_weights_location")} {
164  // instantiate the Track Pair Finder algorithm
165  likelihoodCalc_.setWeightsFile(edm::FileInPath{likelihoodWeights_.c_str()}.fullPath().c_str());
166 
167  cutsFromDB_ = config.getParameter<bool>("usePFThresholdsFromDB");
168  if (cutsFromDB_) {
169  hcalCutsToken_ = esConsumes<HcalPFCuts, HcalPFCutsRcd>(edm::ESInputTag("", "withTopo"));
170  }
171 
173  cfgCone.hOverEConeSize = hOverEConeSize_;
174  if (cfgCone.hOverEConeSize > 0) {
175  cfgCone.onlyBehindCluster = false;
176  cfgCone.checkHcalStatus = false;
177 
178  cfgCone.hbheRecHits = hbheRecHits_;
179 
180  cfgCone.eThresHB = config.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
181  cfgCone.maxSeverityHB = config.getParameter<int>("maxHcalRecHitSeverity");
182  cfgCone.eThresHE = config.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
183  cfgCone.maxSeverityHE = cfgCone.maxSeverityHB;
184  }
185 
186  hcalHelper_ = std::make_unique<ElectronHcalHelper>(cfgCone, consumesCollector());
187 }
188 
189 void ConvertedPhotonProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
190  magneticField_ = &theEventSetup.getData(mFToken_);
191 
192  // Transform Track into TransientTrack (needed by the Vertex fitter)
194 
195  if (cutsFromDB_) {
196  hcalCuts_ = &theEventSetup.getData(hcalCutsToken_);
197  }
198 
199  //
200  // create empty output collections
201  //
202  // Converted photon candidates
203  reco::ConversionCollection outputConvPhotonCollection;
204  // Converted photon candidates
205  reco::ConversionCollection cleanedConversionCollection;
206 
207  // Get the Super Cluster collection in the Barrel
208  bool validBarrelSCHandle = true;
209  auto scBarrelHandle = theEvent.getHandle(scHybridBarrelProducer_);
210  if (!scBarrelHandle.isValid()) {
211  edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the scHybridBarrelProducer";
212  validBarrelSCHandle = false;
213  }
214 
215  // Get the Super Cluster collection in the Endcap
216  bool validEndcapSCHandle = true;
218  theEvent.getByToken(scIslandEndcapProducer_, scEndcapHandle);
219  if (!scEndcapHandle.isValid()) {
220  edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the scIslandEndcapProducer";
221  validEndcapSCHandle = false;
222  }
223 
225  bool validTrackInputs = true;
226  auto outInTrkHandle = theEvent.getHandle(conversionOITrackProducer_);
227  if (!outInTrkHandle.isValid()) {
228  //std::cout << "Error! Can't get the conversionOITrack " << "\n";
229  edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the conversionOITrack "
230  << "\n";
231  validTrackInputs = false;
232  }
233  // LogDebug("ConvertedPhotonProducer")<< "ConvertedPhotonProducer outInTrack collection size " << (*outInTrkHandle).size() << "\n";
234 
236  auto outInTrkSCAssocHandle = theEvent.getHandle(outInTrackSCAssociationCollection_);
237  if (!outInTrkSCAssocHandle.isValid()) {
238  // std::cout << "Error! Can't get the product " << outInTrackSCAssociationCollection_.c_str() <<"\n";
239  edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the outInTrackSCAssociationCollection)";
240  validTrackInputs = false;
241  }
242 
244  auto inOutTrkHandle = theEvent.getHandle(conversionIOTrackProducer_);
245  if (!inOutTrkHandle.isValid()) {
246  // std::cout << "Error! Can't get the conversionIOTrack " << "\n";
247  edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the conversionIOTrack "
248  << "\n";
249  validTrackInputs = false;
250  }
251  // LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer inOutTrack collection size " << (*inOutTrkHandle).size() << "\n";
252 
254 
255  edm::Handle<reco::TrackCollection> generalTrkHandle;
256  if (recoverOneTrackCase_) {
257  theEvent.getByToken(generalTrackProducer_, generalTrkHandle);
258  if (!generalTrkHandle.isValid()) {
259  //std::cout << "Error! Can't get the genralTracks " << "\n";
260  edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the genralTracks "
261  << "\n";
262  }
263  }
264 
266  auto inOutTrkSCAssocHandle = theEvent.getHandle(inOutTrackSCAssociationCollection_);
267  if (!inOutTrkSCAssocHandle.isValid()) {
268  //std::cout << "Error! Can't get the product " << inOutTrackSCAssociationCollection_.c_str() <<"\n";
269  edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the inOutTrackSCAssociationCollection_.c_str()";
270  validTrackInputs = false;
271  }
272 
273  // Get the basic cluster collection in the Barrel
275  theEvent.getByToken(bcBarrelCollection_, bcBarrelHandle);
276  if (!bcBarrelHandle.isValid()) {
277  edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the bcBarrelCollection";
278  }
279 
280  // Get the basic cluster collection in the Endcap
282  theEvent.getByToken(bcEndcapCollection_, bcEndcapHandle);
283  if (!bcEndcapHandle.isValid()) {
284  edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the bcEndcapCollection";
285  }
286 
287  if (validTrackInputs) {
288  //do the conversion:
289  std::vector<reco::TransientTrack> t_outInTrk = transientTrackBuilder_->build(outInTrkHandle);
290  std::vector<reco::TransientTrack> t_inOutTrk = transientTrackBuilder_->build(inOutTrkHandle);
291 
293  std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors> allPairs;
294  allPairs = trackPairFinder_.run(
295  t_outInTrk, outInTrkHandle, outInTrkSCAssocHandle, t_inOutTrk, inOutTrkHandle, inOutTrkSCAssocHandle);
296  //LogDebug("ConvertedPhotonProducer") << "ConvertedPhotonProducer allPairs.size " << allPairs.size() << "\n";
297 
298  hcalHelper_->beginEvent(theEvent, theEventSetup);
299 
300  buildCollections(theEventSetup,
301  scBarrelHandle,
302  bcBarrelHandle,
303  *hcalHelper_,
304  generalTrkHandle,
305  allPairs,
306  outputConvPhotonCollection);
307  buildCollections(theEventSetup,
308  scEndcapHandle,
309  bcEndcapHandle,
310  *hcalHelper_,
311  generalTrkHandle,
312  allPairs,
313  outputConvPhotonCollection);
314  }
315 
316  // put the product in the event
317  auto const conversionHandle =
318  theEvent.emplace(convertedPhotonCollectionPutToken_, std::move(outputConvPhotonCollection));
319 
320  // Loop over barrel and endcap SC collections and fill the photon collection
321  if (validBarrelSCHandle)
322  cleanCollections(scBarrelHandle, conversionHandle, cleanedConversionCollection);
323  if (validEndcapSCHandle)
324  cleanCollections(scEndcapHandle, conversionHandle, cleanedConversionCollection);
325 
326  theEvent.emplace(cleanedConvertedPhotonCollectionPutToken_, std::move(cleanedConversionCollection));
327 }
328 
330  edm::EventSetup const& es,
331  const edm::Handle<edm::View<reco::CaloCluster>>& scHandle,
332  const edm::Handle<edm::View<reco::CaloCluster>>& bcHandle,
333  ElectronHcalHelper const& hcalHelper,
334  const edm::Handle<reco::TrackCollection>& generalTrkHandle,
335  std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors>& allPairs,
336  reco::ConversionCollection& outputConvPhotonCollection)
337 
338 {
339  // instantiate the algorithm for finding the position of the track extrapolation at the Ecal front face
340  ConversionTrackEcalImpactPoint theEcalImpactPositionFinder(magneticField_);
341 
343 
344  std::vector<reco::TransientTrack> t_generalTrk;
346  t_generalTrk = transientTrackBuilder_->build(generalTrkHandle);
347 
348  // Loop over SC in the barrel and reconstruct converted photons
350  for (auto const& aClus : scHandle->ptrs()) {
351  // preselection based in Et and H/E cut
352  if (aClus->energy() / cosh(aClus->eta()) <= minSCEt_)
353  continue;
354  const reco::CaloCluster* pClus = &(*aClus);
355  auto const* sc = dynamic_cast<const reco::SuperCluster*>(pClus);
356  double HoE = hcalHelper.hcalESum(*sc, 0, hcalCuts_) / sc->energy();
357  if (HoE >= maxHOverE_)
358  continue;
360 
361  std::vector<edm::Ref<reco::TrackCollection>> trackPairRef;
362  std::vector<math::XYZPointF> trackInnPos;
363  std::vector<math::XYZVectorF> trackPin;
364  std::vector<math::XYZVectorF> trackPout;
365  float minAppDist = -99;
366 
367  //LogDebug("ConvertedPhotonProducer") << "ConvertedPhotonProducer SC energy " << aClus->energy() << " eta " << aClus->eta() << " phi " << aClus->phi() << "\n";
368 
370  const reco::Particle::Point vtx(0, 0, 0);
371 
372  math::XYZVector direction = aClus->position() - vtx;
373  math::XYZVector momentum = direction.unit() * aClus->energy();
374  const reco::Particle::LorentzVector p4(momentum.x(), momentum.y(), momentum.z(), aClus->energy());
375 
376  if (!allPairs.empty()) {
377  for (auto iPair = allPairs.begin(); iPair != allPairs.end(); ++iPair) {
378  scPtrVec.clear();
379 
380  reco::Vertex theConversionVertex;
381  reco::CaloClusterPtr caloPtr = iPair->second;
382  if (!(aClus == caloPtr))
383  continue;
384 
385  scPtrVec.push_back(aClus);
386 
387  std::vector<math::XYZPointF> trkPositionAtEcal = theEcalImpactPositionFinder.find(iPair->first, bcHandle);
388  std::vector<reco::CaloClusterPtr> matchingBC = theEcalImpactPositionFinder.matchingBC();
389 
390  minAppDist = -99;
391  const std::string metname = "ConvertedPhotons|ConvertedPhotonProducer";
392  if ((iPair->first).size() > 1) {
393  try {
394  vertexFinder_.run(iPair->first, theConversionVertex);
395 
396  } catch (cms::Exception& e) {
397  //std::cout << " cms::Exception caught in ConvertedPhotonProducer::produce" << "\n" ;
398  edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
399  << e.explainSelf();
400  }
401 
402  // Old TwoTrackMinimumDistance md;
403  // Old md.calculate ( (iPair->first)[0].initialFreeState(), (iPair->first)[1].initialFreeState() );
404  // Old minAppDist = md.distance();
405 
406  /*
407  for ( unsigned int i=0; i< matchingBC.size(); ++i) {
408  if ( matchingBC[i].isNull() ) std::cout << " This ref to BC is null: skipping " << "\n";
409  else
410  std::cout << " BC energy " << matchingBC[i]->energy() << "\n";
411  }
412  */
413 
415  trackPairRef.clear();
416  trackInnPos.clear();
417  trackPin.clear();
418  trackPout.clear();
419 
420  for (std::vector<reco::TransientTrack>::const_iterator iTk = (iPair->first).begin();
421  iTk != (iPair->first).end();
422  ++iTk) {
423  //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer Transient Tracks in the pair charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";
424 
425  auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
426  reco::TrackRef myTkRef = ttt->persistentTrackRef();
427 
428  //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer Ref to Rec Tracks in the pair charge " << myTkRef->charge() << " Num of RecHits " << myTkRef->recHitsSize() << " inner momentum " << myTkRef->innerMomentum() << "\n";
429  if (myTkRef->extra().isNonnull()) {
430  trackInnPos.push_back(toFConverterP(myTkRef->innerPosition()));
431  trackPin.push_back(toFConverterV(myTkRef->innerMomentum()));
432  trackPout.push_back(toFConverterV(myTkRef->outerMomentum()));
433  }
434  trackPairRef.push_back(myTkRef);
435  }
436 
437  // std::cout << " ConvertedPhotonProducer trackPin size " << trackPin.size() << std::endl;
438  //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer SC energy " << aClus->energy() << "\n";
439  //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer photon p4 " << p4 << "\n";
440  //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer vtx " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
441  if (theConversionVertex.isValid()) {
442  //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer theConversionVertex " << theConversionVertex.position().x() << " " << theConversionVertex.position().y() << " " << theConversionVertex.position().z() << "\n";
443  }
444  //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer trackPairRef " << trackPairRef.size() << "\n";
445 
446  minAppDist = calculateMinApproachDistance(trackPairRef[0], trackPairRef[1]);
447 
448  double like = -999.;
449  reco::Conversion newCandidate(scPtrVec,
450  trackPairRef,
451  trkPositionAtEcal,
452  theConversionVertex,
453  matchingBC,
454  minAppDist,
455  trackInnPos,
456  trackPin,
457  trackPout,
458  like,
459  algo);
460  like = likelihoodCalc_.calculateLikelihood(newCandidate);
461  // std::cout << "like = " << like << std::endl;
462  newCandidate.setMVAout(like);
463  outputConvPhotonCollection.push_back(newCandidate);
464 
465  //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer Put the ConvertedPhotonCollection a candidate in the Barrel " << "\n";
466 
467  } else {
468  // std::cout << " ConvertedPhotonProducer case with only one track found " << "\n";
469 
470  //std::cout << " ConvertedPhotonProducer recovering one track " << "\n";
471  trackPairRef.clear();
472  trackInnPos.clear();
473  trackPin.clear();
474  trackPout.clear();
475  std::vector<reco::TransientTrack>::const_iterator iTk = (iPair->first).begin();
476  //std::cout << " ConvertedPhotonProducer Transient Tracks in the pair charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << " pt " << sqrt(iTk->track().innerMomentum().perp2()) << "\n";
477  auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
478  reco::TrackRef myTk = ttt->persistentTrackRef();
479  if (myTk->extra().isNonnull()) {
480  trackInnPos.push_back(toFConverterP(myTk->innerPosition()));
481  trackPin.push_back(toFConverterV(myTk->innerMomentum()));
482  trackPout.push_back(toFConverterV(myTk->outerMomentum()));
483  }
484  trackPairRef.push_back(myTk);
485  //std::cout << " Provenance " << myTk->algoName() << std::endl;
486 
487  if (recoverOneTrackCase_) {
488  float theta1 = myTk->innerMomentum().Theta();
489  float dCot = 999.;
490  float dCotTheta = -999.;
491  reco::TrackRef goodRef;
492  for (auto const& tran : t_generalTrk) {
493  auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(tran.basicTransientTrack());
494  reco::TrackRef trRef = ttt->persistentTrackRef();
495  if (trRef->charge() * myTk->charge() > 0)
496  continue;
497  float dEta = trRef->eta() - myTk->eta();
498  float dPhi = trRef->phi() - myTk->phi();
500  continue;
501  float theta2 = trRef->innerMomentum().Theta();
502  dCotTheta = 1. / tan(theta1) - 1. / tan(theta2);
503  // std::cout << " ConvertedPhotonProducer recovering general transient track charge " << trRef->charge() << " momentum " << trRef->innerMomentum() << " dcotTheta " << fabs(dCotTheta) << std::endl;
504  if (fabs(dCotTheta) < dCot) {
505  dCot = fabs(dCotTheta);
506  goodRef = trRef;
507  }
508  }
509 
510  if (goodRef.isNonnull()) {
511  minAppDist = calculateMinApproachDistance(myTk, goodRef);
512 
513  // std::cout << " ConvertedPhotonProducer chosen dCotTheta " << fabs(dCotTheta) << std::endl;
514  if (fabs(dCotTheta) < deltaCotCut_ && minAppDist > minApproachDisCut_) {
515  trackInnPos.push_back(toFConverterP(goodRef->innerPosition()));
516  trackPin.push_back(toFConverterV(goodRef->innerMomentum()));
517  trackPout.push_back(toFConverterV(goodRef->outerMomentum()));
518  trackPairRef.push_back(goodRef);
519  // std::cout << " ConvertedPhotonProducer adding opposite charge track from generalTrackCollection charge " << goodRef ->charge() << " pt " << sqrt(goodRef->innerMomentum().perp2()) << " trackPairRef size " << trackPairRef.size() << std::endl;
520  //std::cout << " Track Provenenance " << goodRef->algoName() << std::endl;
521 
522  try {
523  vertexFinder_.run(iPair->first, theConversionVertex);
524 
525  } catch (cms::Exception& e) {
526  //std::cout << " cms::Exception caught in ConvertedPhotonProducer::produce" << "\n" ;
527  edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
528  << e.explainSelf();
529  }
530  }
531  }
532 
533  } // bool On/Off one track case recovery using generalTracks
534  const double like = -999.;
535  outputConvPhotonCollection.emplace_back(scPtrVec,
536  trackPairRef,
537  trkPositionAtEcal,
538  theConversionVertex,
539  matchingBC,
540  minAppDist,
541  trackInnPos,
542  trackPin,
543  trackPout,
544  like,
545  algo);
546  auto& newCandidate = outputConvPhotonCollection.back();
547  newCandidate.setMVAout(likelihoodCalc_.calculateLikelihood(newCandidate));
548 
549  } // case with only on track: looking in general tracks
550  }
551  }
552  }
553 }
554 
557  reco::ConversionCollection& outputConversionCollection) {
558  reco::Conversion* newCandidate = nullptr;
559  for (auto const& aClus : scHandle->ptrs()) {
560  // SC energy preselection
561  if (aClus->energy() / cosh(aClus->eta()) <= minSCEt_)
562  continue;
563 
564  if (conversionHandle.isValid()) {
565  if (risolveAmbiguity_) {
566  std::vector<reco::ConversionRef> bestRef = solveAmbiguity(conversionHandle, aClus);
567 
568  for (std::vector<reco::ConversionRef>::iterator iRef = bestRef.begin(); iRef != bestRef.end(); iRef++) {
569  if (iRef->isNonnull()) {
570  newCandidate = (*iRef)->clone();
571  outputConversionCollection.push_back(*newCandidate);
572  delete newCandidate;
573  }
574  }
575 
576  } else {
577  for (unsigned int icp = 0; icp < conversionHandle->size(); icp++) {
579  if (!(aClus.id() == cpRef->caloCluster()[0].id() && aClus.key() == cpRef->caloCluster()[0].key()))
580  continue;
581  if (!cpRef->isConverted())
582  continue;
583  if (cpRef->nTracks() < 2)
584  continue;
585  newCandidate = (&(*cpRef))->clone();
586  outputConversionCollection.push_back(*newCandidate);
587  delete newCandidate;
588  }
589 
590  } // solve or not the ambiguity of many conversion candidates
591  }
592  }
593 }
594 
595 std::vector<reco::ConversionRef> ConvertedPhotonProducer::solveAmbiguity(
597  std::multimap<double, reco::ConversionRef, std::greater<double>> convMap;
598 
599  for (unsigned int icp = 0; icp < conversionHandle->size(); icp++) {
601 
602  //std::cout << " cpRef " << cpRef->nTracks() << " " << cpRef ->caloCluster()[0]->energy() << std::endl;
603  if (!(scRef.id() == cpRef->caloCluster()[0].id() && scRef.key() == cpRef->caloCluster()[0].key()))
604  continue;
605  if (!cpRef->isConverted())
606  continue;
607  double like = cpRef->MVAout();
608  if (cpRef->nTracks() < 2)
609  continue;
610  // std::cout << " Like " << like << std::endl;
611  convMap.emplace(like, cpRef);
612  }
613 
614  // std::cout << " convMap size " << convMap.size() << std::endl;
615 
616  std::vector<reco::ConversionRef> bestRefs;
617  for (auto iMap = convMap.begin(); iMap != convMap.end(); iMap++) {
618  // std::cout << " Like list in the map " << iMap->first << " " << (iMap->second)->EoverP() << std::endl;
619  bestRefs.push_back(iMap->second);
620  if (int(bestRefs.size()) == maxNumOfCandidates_)
621  break;
622  }
623 
624  return bestRefs;
625 }
626 
628  const reco::TrackRef& track2) {
629  double x1, x2, y1, y2;
630  double xx_1 = track1->innerPosition().x(), yy_1 = track1->innerPosition().y(), zz_1 = track1->innerPosition().z();
631  double xx_2 = track2->innerPosition().x(), yy_2 = track2->innerPosition().y(), zz_2 = track2->innerPosition().z();
632  double radius1 =
633  track1->innerMomentum().Rho() / (.3 * (magneticField_->inTesla(GlobalPoint(xx_1, yy_1, zz_1)).z())) * 100;
634  double radius2 =
635  track2->innerMomentum().Rho() / (.3 * (magneticField_->inTesla(GlobalPoint(xx_2, yy_2, zz_2)).z())) * 100;
636  getCircleCenter(track1, radius1, x1, y1);
637  getCircleCenter(track2, radius2, x2, y2);
638 
639  return std::sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) - radius1 - radius2;
640 }
641 
642 void ConvertedPhotonProducer::getCircleCenter(const reco::TrackRef& tk, double r, double& x0, double& y0) {
643  double x1, y1, phi;
644  x1 = tk->innerPosition().x(); //inner position and inner momentum need track Extra!
645  y1 = tk->innerPosition().y();
646  phi = tk->innerMomentum().phi();
647  const int charge = tk->charge();
648  x0 = x1 + r * sin(phi) * charge;
649  y0 = y1 - r * cos(phi) * charge;
650 }
edm::EDPutTokenT< reco::ConversionCollection > cleanedConvertedPhotonCollectionPutToken_
edm::ESGetToken< HcalPFCuts, HcalPFCutsRcd > hcalCutsToken_
edm::EDGetTokenT< reco::TrackCaloClusterPtrAssociation > outInTrackSCAssociationCollection_
ESGetTokenH3DDVariant esConsumes(std::string const &Record, edm::ConsumesCollector &)
Definition: DeDxTools.cc:283
edm::EDGetTokenT< edm::View< reco::CaloCluster > > bcEndcapCollection_
TransientVertex run(const std::vector< reco::TransientTrack > &pair)
ProductID id() const
Accessor for product ID.
Definition: Ref.h:244
T const & getData(const ESGetToken< T, R > &iToken) const noexcept(false)
Definition: EventSetup.h:119
ProductID id() const
Accessor for product ID.
Definition: Ptr.h:160
std::string fullPath() const
Definition: FileInPath.cc:161
edm::EDPutTokenT< reco::ConversionCollection > convertedPhotonCollectionPutToken_
T z() const
Definition: PV3DBase.h:61
const std::string metname
virtual GlobalVector inTesla(const GlobalPoint &gp) const =0
Field value ad specified global point, in Tesla.
void push_back(Ptr< T > const &iPtr)
Definition: PtrVector.h:152
edm::EDGetTokenT< reco::TrackCollection > conversionIOTrackProducer_
Sin< T >::type sin(const T &t)
Definition: Sin.h:22
Global3DPoint GlobalPoint
Definition: GlobalPoint.h:10
edm::EDGetTokenT< edm::View< reco::CaloCluster > > scIslandEndcapProducer_
static ConversionAlgorithm algoByName(const std::string &name)
Definition: Conversion.cc:139
bool isNonnull() const
Checks for non-null.
Definition: Ref.h:238
edm::ESGetToken< CaloGeometry, CaloGeometryRecord > caloGeomToken_
bool getByToken(EDGetToken token, Handle< PROD > &result) const
Definition: Event.h:528
Definition: config.py:1
std::unique_ptr< ElectronHcalHelper > hcalHelper_
edm::Ptr< CaloCluster > CaloClusterPtr
ROOT::Math::PositionVector3D< ROOT::Math::Cartesian3D< float > > XYZPointF
point in space with cartesian internal representation
Definition: Point3D.h:10
edm::EDGetTokenT< HBHERecHitCollection > hbheRecHits_
Log< level::Error, false > LogError
MagneticField const * magneticField_
edm::EDGetTokenT< HBHERecHitCollection > hbheRecHits
key_type key() const
Accessor for product key.
Definition: Ref.h:250
std::vector< Conversion > ConversionCollection
collectin of Conversion objects
Definition: ConversionFwd.h:9
std::vector< reco::CaloClusterPtr > matchingBC() const
reco::TransientTrack build(const reco::Track *p) const
edm::EDGetTokenT< reco::TrackCollection > generalTrackProducer_
math::XYZVectorF toFConverterV(const math::XYZVector &val)
void buildCollections(edm::EventSetup const &es, const edm::Handle< edm::View< reco::CaloCluster >> &scHandle, const edm::Handle< edm::View< reco::CaloCluster >> &bcHandle, ElectronHcalHelper const &hcalHelper, const edm::Handle< reco::TrackCollection > &trkHandle, std::map< std::vector< reco::TransientTrack >, reco::CaloClusterPtr, CompareTwoTracksVectors > &allPairs, reco::ConversionCollection &outputConvPhotonCollection)
T sqrt(T t)
Definition: SSEVec.h:19
Cos< T >::type cos(const T &t)
Definition: Cos.h:22
Conversion * clone() const
returns a clone of the candidate
Definition: Conversion.cc:148
Tan< T >::type tan(const T &t)
Definition: Tan.h:22
math::XYZPoint Point
point in the space
Definition: Particle.h:25
std::vector< reco::ConversionRef > solveAmbiguity(const edm::OrphanHandle< reco::ConversionCollection > &conversionHandle, reco::CaloClusterPtr const &sc)
edm::ESGetToken< TransientTrackBuilder, TransientTrackRecord > transientTrackToken_
float calculateMinApproachDistance(const reco::TrackRef &track1, const reco::TrackRef &track2)
edm::EDGetTokenT< reco::TrackCollection > conversionOITrackProducer_
#define DEFINE_FWK_MODULE(type)
Definition: MakerMacros.h:16
edm::EDGetTokenT< reco::TrackCaloClusterPtrAssociation > inOutTrackSCAssociationCollection_
std::vector< math::XYZPointF > find(const std::vector< reco::TransientTrack > &tracks, const edm::Handle< edm::View< reco::CaloCluster > > &bcHandle)
ROOT::Math::DisplacementVector3D< ROOT::Math::Cartesian3D< float > > XYZVectorF
spatial vector with cartesian internal representation
Definition: Vector3D.h:16
ConversionLikelihoodCalculator likelihoodCalc_
double calculateLikelihood(reco::ConversionRef conversion)
ProductIndex id() const
Definition: ProductID.h:35
edm::EDGetTokenT< edm::View< reco::CaloCluster > > bcBarrelCollection_
edm::EDGetTokenT< edm::View< reco::CaloCluster > > scHybridBarrelProducer_
OrphanHandle< PROD > emplace(EDPutTokenT< PROD > token, Args &&... args)
puts a new product
Definition: Event.h:431
XYZVectorD XYZVector
spatial vector with cartesian internal representation
Definition: Vector3D.h:31
XYZPointD XYZPoint
point in space with cartesian internal representation
Definition: Point3D.h:12
void produce(edm::Event &evt, const edm::EventSetup &es) override
ConversionVertexFinder vertexFinder_
void cleanCollections(const edm::Handle< edm::View< reco::CaloCluster >> &scHandle, const edm::OrphanHandle< reco::ConversionCollection > &conversionHandle, reco::ConversionCollection &outputCollection)
TEveGeoShape * clone(const TEveElement *element, TEveElement *parent)
Definition: eve_macros.cc:135
void setMVAout(const float &mva)
set the value of the TMVA output
Definition: Conversion.h:161
bool isValid() const
Definition: HandleBase.h:70
ConvertedPhotonProducer(const edm::ParameterSet &ps)
ConversionTrackPairFinder trackPairFinder_
TransientTrackBuilder const * transientTrackBuilder_
void clear()
Clear the PtrVector.
Definition: PtrVectorBase.h:81
Handle< PROD > getHandle(EDGetTokenT< PROD > token) const
Definition: Event.h:552
key_type key() const
Definition: Ptr.h:165
EgammaHcalIsolation::arrayHB eThresHB
math::XYZPointF toFConverterP(const math::XYZPoint &val)
EgammaHcalIsolation::arrayHE eThresHE
Log< level::Warning, false > LogWarning
double hcalESum(const reco::SuperCluster &, int depth, const HcalPFCuts *hcalCuts) const
std::array< double, 4 > arrayHB
std::map< std::vector< reco::TransientTrack >, reco::CaloClusterPtr, CompareTwoTracksVectors > run(const std::vector< reco::TransientTrack > &outIn, const edm::Handle< reco::TrackCollection > &outInTrkHandle, const edm::Handle< reco::TrackCaloClusterPtrAssociation > &outInTrackSCAssH, const std::vector< reco::TransientTrack > &inOut, const edm::Handle< reco::TrackCollection > &inOutTrkHandle, const edm::Handle< reco::TrackCaloClusterPtrAssociation > &inOutTrackSCAssH)
math::XYZTLorentzVector LorentzVector
Lorentz vector.
Definition: Particle.h:21
edm::ESGetToken< MagneticField, IdealMagneticFieldRecord > mFToken_
def move(src, dest)
Definition: eostools.py:511
bool isValid() const
Tells whether the vertex is valid.
Definition: Vertex.h:73
std::array< double, 7 > arrayHE
void getCircleCenter(const reco::TrackRef &tk, double r, double &x0, double &y0)