File indexing completed on 2023-05-08 23:19:13
0001
0002
0003
0004
0005
0006
0007
0008 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0009 #include "DataFormats/CaloTowers/interface/CaloTowerCollection.h"
0010 #include "DataFormats/Common/interface/Handle.h"
0011 #include "DataFormats/Common/interface/View.h"
0012 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
0013 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0014 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0015 #include "DataFormats/EgammaTrackReco/interface/TrackCaloClusterAssociation.h"
0016 #include "DataFormats/TrackReco/interface/Track.h"
0017 #include "DataFormats/TrackReco/interface/TrackExtra.h"
0018 #include "DataFormats/VertexReco/interface/Vertex.h"
0019 #include "FWCore/Framework/interface/ConsumesCollector.h"
0020 #include "FWCore/Framework/interface/Event.h"
0021 #include "FWCore/Framework/interface/EventSetup.h"
0022 #include "FWCore/Framework/interface/stream/EDProducer.h"
0023 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0024 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0025 #include "FWCore/Utilities/interface/ESGetToken.h"
0026 #include "FWCore/Utilities/interface/Exception.h"
0027 #include "Geometry/CaloGeometry/interface/CaloGeometry.h"
0028 #include "Geometry/Records/interface/CaloGeometryRecord.h"
0029 #include "Geometry/TrackerGeometryBuilder/interface/TrackerGeometry.h"
0030 #include "MagneticField/Engine/interface/MagneticField.h"
0031 #include "MagneticField/Records/interface/IdealMagneticFieldRecord.h"
0032 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionTrackEcalImpactPoint.h"
0033 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionTrackPairFinder.h"
0034 #include "RecoEgamma/EgammaPhotonAlgos/interface/ConversionVertexFinder.h"
0035 #include "RecoEgamma/EgammaTools/interface/ConversionLikelihoodCalculator.h"
0036 #include "TrackingTools/Records/interface/TransientTrackRecord.h"
0037 #include "TrackingTools/TransientTrack/interface/TrackTransientTrack.h"
0038 #include "TrackingTools/TransientTrack/interface/TransientTrack.h"
0039 #include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h"
0040 #include "RecoEgamma/EgammaElectronAlgos/interface/ElectronHcalHelper.h"
0041
0042 #include <vector>
0043
0044 class ConvertedPhotonProducer : public edm::stream::EDProducer<> {
0045 public:
0046 ConvertedPhotonProducer(const edm::ParameterSet& ps);
0047
0048 void beginRun(edm::Run const&, const edm::EventSetup& es) final;
0049 void produce(edm::Event& evt, const edm::EventSetup& es) override;
0050
0051 private:
0052 void buildCollections(
0053 edm::EventSetup const& es,
0054 const edm::Handle<edm::View<reco::CaloCluster>>& scHandle,
0055 const edm::Handle<edm::View<reco::CaloCluster>>& bcHandle,
0056 ElectronHcalHelper const& hcalHelper,
0057 const edm::Handle<reco::TrackCollection>& trkHandle,
0058 std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors>& allPairs,
0059 reco::ConversionCollection& outputConvPhotonCollection);
0060 void cleanCollections(const edm::Handle<edm::View<reco::CaloCluster>>& scHandle,
0061 const edm::OrphanHandle<reco::ConversionCollection>& conversionHandle,
0062 reco::ConversionCollection& outputCollection);
0063
0064 std::vector<reco::ConversionRef> solveAmbiguity(const edm::OrphanHandle<reco::ConversionCollection>& conversionHandle,
0065 reco::CaloClusterPtr const& sc);
0066
0067 float calculateMinApproachDistance(const reco::TrackRef& track1, const reco::TrackRef& track2);
0068 void getCircleCenter(const reco::TrackRef& tk, double r, double& x0, double& y0);
0069
0070 edm::EDGetTokenT<reco::TrackCollection> conversionOITrackProducer_;
0071 edm::EDGetTokenT<reco::TrackCollection> conversionIOTrackProducer_;
0072
0073 edm::EDGetTokenT<reco::TrackCaloClusterPtrAssociation> outInTrackSCAssociationCollection_;
0074 edm::EDGetTokenT<reco::TrackCaloClusterPtrAssociation> inOutTrackSCAssociationCollection_;
0075
0076 edm::EDGetTokenT<reco::TrackCollection> generalTrackProducer_;
0077
0078
0079 edm::EDPutTokenT<reco::ConversionCollection> convertedPhotonCollectionPutToken_;
0080 edm::EDPutTokenT<reco::ConversionCollection> cleanedConvertedPhotonCollectionPutToken_;
0081
0082 edm::EDGetTokenT<edm::View<reco::CaloCluster>> bcBarrelCollection_;
0083 edm::EDGetTokenT<edm::View<reco::CaloCluster>> bcEndcapCollection_;
0084 edm::EDGetTokenT<edm::View<reco::CaloCluster>> scHybridBarrelProducer_;
0085 edm::EDGetTokenT<edm::View<reco::CaloCluster>> scIslandEndcapProducer_;
0086 edm::EDGetTokenT<HBHERecHitCollection> hbheRecHits_;
0087
0088 MagneticField const* magneticField_;
0089 TransientTrackBuilder const* transientTrackBuilder_;
0090
0091 edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0092 edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> mFToken_;
0093 edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackToken_;
0094
0095 ConversionTrackPairFinder trackPairFinder_;
0096 ConversionVertexFinder vertexFinder_;
0097 std::string algoName_;
0098
0099 double hOverEConeSize_;
0100 double maxHOverE_;
0101 double minSCEt_;
0102 bool recoverOneTrackCase_;
0103 double dRForConversionRecovery_;
0104 double deltaCotCut_;
0105 double minApproachDisCut_;
0106 int maxNumOfCandidates_;
0107 bool risolveAmbiguity_;
0108
0109 std::unique_ptr<ElectronHcalHelper> hcalHelper_;
0110
0111 ConversionLikelihoodCalculator likelihoodCalc_;
0112 std::string likelihoodWeights_;
0113
0114 math::XYZPointF toFConverterP(const math::XYZPoint& val) { return math::XYZPointF(val.x(), val.y(), val.z()); }
0115
0116 math::XYZVectorF toFConverterV(const math::XYZVector& val) { return math::XYZVectorF(val.x(), val.y(), val.z()); }
0117 };
0118
0119 #include "FWCore/Framework/interface/MakerMacros.h"
0120 DEFINE_FWK_MODULE(ConvertedPhotonProducer);
0121
0122 ConvertedPhotonProducer::ConvertedPhotonProducer(const edm::ParameterSet& config)
0123 : conversionOITrackProducer_{consumes(config.getParameter<std::string>("conversionOITrackProducer"))},
0124 conversionIOTrackProducer_{consumes(config.getParameter<std::string>("conversionIOTrackProducer"))},
0125 outInTrackSCAssociationCollection_{consumes({config.getParameter<std::string>("conversionOITrackProducer"),
0126 config.getParameter<std::string>("outInTrackSCAssociation")})},
0127 inOutTrackSCAssociationCollection_{consumes({config.getParameter<std::string>("conversionIOTrackProducer"),
0128 config.getParameter<std::string>("inOutTrackSCAssociation")})},
0129
0130 generalTrackProducer_{consumes(config.getParameter<edm::InputTag>("generalTracksSrc"))},
0131 convertedPhotonCollectionPutToken_{
0132 produces<reco::ConversionCollection>(config.getParameter<std::string>("convertedPhotonCollection"))},
0133 cleanedConvertedPhotonCollectionPutToken_{
0134 produces<reco::ConversionCollection>(config.getParameter<std::string>("cleanedConvertedPhotonCollection"))},
0135
0136 bcBarrelCollection_{consumes(config.getParameter<edm::InputTag>("bcBarrelCollection"))},
0137 bcEndcapCollection_{consumes(config.getParameter<edm::InputTag>("bcEndcapCollection"))},
0138 scHybridBarrelProducer_{consumes(config.getParameter<edm::InputTag>("scHybridBarrelProducer"))},
0139 scIslandEndcapProducer_{consumes(config.getParameter<edm::InputTag>("scIslandEndcapProducer"))},
0140 hbheRecHits_{consumes(config.getParameter<edm::InputTag>("hbheRecHits"))},
0141 caloGeomToken_{esConsumes()},
0142 mFToken_{esConsumes<MagneticField, IdealMagneticFieldRecord, edm::Transition::BeginRun>()},
0143 transientTrackToken_{esConsumes<TransientTrackBuilder, TransientTrackRecord, edm::Transition::BeginRun>(
0144 edm::ESInputTag("", "TransientTrackBuilder"))},
0145 vertexFinder_{config},
0146 algoName_{config.getParameter<std::string>("AlgorithmName")},
0147
0148 hOverEConeSize_{config.getParameter<double>("hOverEConeSize")},
0149 maxHOverE_{config.getParameter<double>("maxHOverE")},
0150 minSCEt_{config.getParameter<double>("minSCEt")},
0151 recoverOneTrackCase_{config.getParameter<bool>("recoverOneTrackCase")},
0152 dRForConversionRecovery_{config.getParameter<double>("dRForConversionRecovery")},
0153 deltaCotCut_{config.getParameter<double>("deltaCotCut")},
0154 minApproachDisCut_{config.getParameter<double>("minApproachDisCut")},
0155
0156 maxNumOfCandidates_{config.getParameter<int>("maxNumOfCandidates")},
0157 risolveAmbiguity_{config.getParameter<bool>("risolveConversionAmbiguity")},
0158 likelihoodWeights_{config.getParameter<std::string>("MVA_weights_location")} {
0159
0160 likelihoodCalc_.setWeightsFile(edm::FileInPath{likelihoodWeights_.c_str()}.fullPath().c_str());
0161
0162 ElectronHcalHelper::Configuration cfgCone;
0163 cfgCone.hOverEConeSize = hOverEConeSize_;
0164 if (cfgCone.hOverEConeSize > 0) {
0165 cfgCone.onlyBehindCluster = false;
0166 cfgCone.checkHcalStatus = false;
0167
0168 cfgCone.hbheRecHits = hbheRecHits_;
0169
0170 cfgCone.eThresHB = config.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0171 cfgCone.maxSeverityHB = config.getParameter<int>("maxHcalRecHitSeverity");
0172 cfgCone.eThresHE = config.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0173 cfgCone.maxSeverityHE = cfgCone.maxSeverityHB;
0174 }
0175
0176 hcalHelper_ = std::make_unique<ElectronHcalHelper>(cfgCone, consumesCollector());
0177 }
0178
0179 void ConvertedPhotonProducer::beginRun(edm::Run const& r, edm::EventSetup const& theEventSetup) {
0180 magneticField_ = &theEventSetup.getData(mFToken_);
0181
0182
0183 transientTrackBuilder_ = &theEventSetup.getData(transientTrackToken_);
0184 }
0185
0186 void ConvertedPhotonProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
0187
0188
0189
0190
0191 reco::ConversionCollection outputConvPhotonCollection;
0192
0193 reco::ConversionCollection cleanedConversionCollection;
0194
0195
0196 bool validBarrelSCHandle = true;
0197 auto scBarrelHandle = theEvent.getHandle(scHybridBarrelProducer_);
0198 if (!scBarrelHandle.isValid()) {
0199 edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the scHybridBarrelProducer";
0200 validBarrelSCHandle = false;
0201 }
0202
0203
0204 bool validEndcapSCHandle = true;
0205 edm::Handle<edm::View<reco::CaloCluster>> scEndcapHandle;
0206 theEvent.getByToken(scIslandEndcapProducer_, scEndcapHandle);
0207 if (!scEndcapHandle.isValid()) {
0208 edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the scIslandEndcapProducer";
0209 validEndcapSCHandle = false;
0210 }
0211
0212
0213 bool validTrackInputs = true;
0214 auto outInTrkHandle = theEvent.getHandle(conversionOITrackProducer_);
0215 if (!outInTrkHandle.isValid()) {
0216
0217 edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the conversionOITrack "
0218 << "\n";
0219 validTrackInputs = false;
0220 }
0221
0222
0223
0224 auto outInTrkSCAssocHandle = theEvent.getHandle(outInTrackSCAssociationCollection_);
0225 if (!outInTrkSCAssocHandle.isValid()) {
0226
0227 edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the outInTrackSCAssociationCollection)";
0228 validTrackInputs = false;
0229 }
0230
0231
0232 auto inOutTrkHandle = theEvent.getHandle(conversionIOTrackProducer_);
0233 if (!inOutTrkHandle.isValid()) {
0234
0235 edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the conversionIOTrack "
0236 << "\n";
0237 validTrackInputs = false;
0238 }
0239
0240
0241
0242
0243 edm::Handle<reco::TrackCollection> generalTrkHandle;
0244 if (recoverOneTrackCase_) {
0245 theEvent.getByToken(generalTrackProducer_, generalTrkHandle);
0246 if (!generalTrkHandle.isValid()) {
0247
0248 edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the genralTracks "
0249 << "\n";
0250 }
0251 }
0252
0253
0254 auto inOutTrkSCAssocHandle = theEvent.getHandle(inOutTrackSCAssociationCollection_);
0255 if (!inOutTrkSCAssocHandle.isValid()) {
0256
0257 edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the inOutTrackSCAssociationCollection_.c_str()";
0258 validTrackInputs = false;
0259 }
0260
0261
0262 edm::Handle<edm::View<reco::CaloCluster>> bcBarrelHandle;
0263 theEvent.getByToken(bcBarrelCollection_, bcBarrelHandle);
0264 if (!bcBarrelHandle.isValid()) {
0265 edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the bcBarrelCollection";
0266 }
0267
0268
0269 edm::Handle<edm::View<reco::CaloCluster>> bcEndcapHandle;
0270 theEvent.getByToken(bcEndcapCollection_, bcEndcapHandle);
0271 if (!bcEndcapHandle.isValid()) {
0272 edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the bcEndcapCollection";
0273 }
0274
0275 if (validTrackInputs) {
0276
0277 std::vector<reco::TransientTrack> t_outInTrk = transientTrackBuilder_->build(outInTrkHandle);
0278 std::vector<reco::TransientTrack> t_inOutTrk = transientTrackBuilder_->build(inOutTrkHandle);
0279
0280
0281 std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors> allPairs;
0282 allPairs = trackPairFinder_.run(
0283 t_outInTrk, outInTrkHandle, outInTrkSCAssocHandle, t_inOutTrk, inOutTrkHandle, inOutTrkSCAssocHandle);
0284
0285
0286 hcalHelper_->beginEvent(theEvent, theEventSetup);
0287
0288 buildCollections(theEventSetup,
0289 scBarrelHandle,
0290 bcBarrelHandle,
0291 *hcalHelper_,
0292 generalTrkHandle,
0293 allPairs,
0294 outputConvPhotonCollection);
0295 buildCollections(theEventSetup,
0296 scEndcapHandle,
0297 bcEndcapHandle,
0298 *hcalHelper_,
0299 generalTrkHandle,
0300 allPairs,
0301 outputConvPhotonCollection);
0302 }
0303
0304
0305 auto const conversionHandle =
0306 theEvent.emplace(convertedPhotonCollectionPutToken_, std::move(outputConvPhotonCollection));
0307
0308
0309 if (validBarrelSCHandle)
0310 cleanCollections(scBarrelHandle, conversionHandle, cleanedConversionCollection);
0311 if (validEndcapSCHandle)
0312 cleanCollections(scEndcapHandle, conversionHandle, cleanedConversionCollection);
0313
0314 theEvent.emplace(cleanedConvertedPhotonCollectionPutToken_, std::move(cleanedConversionCollection));
0315 }
0316
0317 void ConvertedPhotonProducer::buildCollections(
0318 edm::EventSetup const& es,
0319 const edm::Handle<edm::View<reco::CaloCluster>>& scHandle,
0320 const edm::Handle<edm::View<reco::CaloCluster>>& bcHandle,
0321 ElectronHcalHelper const& hcalHelper,
0322 const edm::Handle<reco::TrackCollection>& generalTrkHandle,
0323 std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors>& allPairs,
0324 reco::ConversionCollection& outputConvPhotonCollection)
0325
0326 {
0327
0328 ConversionTrackEcalImpactPoint theEcalImpactPositionFinder(magneticField_);
0329
0330 reco::Conversion::ConversionAlgorithm algo = reco::Conversion::algoByName(algoName_);
0331
0332 std::vector<reco::TransientTrack> t_generalTrk;
0333 if (recoverOneTrackCase_)
0334 t_generalTrk = transientTrackBuilder_->build(generalTrkHandle);
0335
0336
0337 reco::CaloClusterPtrVector scPtrVec;
0338 for (auto const& aClus : scHandle->ptrs()) {
0339
0340 if (aClus->energy() / cosh(aClus->eta()) <= minSCEt_)
0341 continue;
0342 const reco::CaloCluster* pClus = &(*aClus);
0343 auto const* sc = dynamic_cast<const reco::SuperCluster*>(pClus);
0344 double HoE = hcalHelper.hcalESum(*sc, 0) / sc->energy();
0345 if (HoE >= maxHOverE_)
0346 continue;
0347
0348
0349 std::vector<edm::Ref<reco::TrackCollection>> trackPairRef;
0350 std::vector<math::XYZPointF> trackInnPos;
0351 std::vector<math::XYZVectorF> trackPin;
0352 std::vector<math::XYZVectorF> trackPout;
0353 float minAppDist = -99;
0354
0355
0356
0357
0358 const reco::Particle::Point vtx(0, 0, 0);
0359
0360 math::XYZVector direction = aClus->position() - vtx;
0361 math::XYZVector momentum = direction.unit() * aClus->energy();
0362 const reco::Particle::LorentzVector p4(momentum.x(), momentum.y(), momentum.z(), aClus->energy());
0363
0364 if (!allPairs.empty()) {
0365 for (auto iPair = allPairs.begin(); iPair != allPairs.end(); ++iPair) {
0366 scPtrVec.clear();
0367
0368 reco::Vertex theConversionVertex;
0369 reco::CaloClusterPtr caloPtr = iPair->second;
0370 if (!(aClus == caloPtr))
0371 continue;
0372
0373 scPtrVec.push_back(aClus);
0374
0375 std::vector<math::XYZPointF> trkPositionAtEcal = theEcalImpactPositionFinder.find(iPair->first, bcHandle);
0376 std::vector<reco::CaloClusterPtr> matchingBC = theEcalImpactPositionFinder.matchingBC();
0377
0378 minAppDist = -99;
0379 const std::string metname = "ConvertedPhotons|ConvertedPhotonProducer";
0380 if ((iPair->first).size() > 1) {
0381 try {
0382 vertexFinder_.run(iPair->first, theConversionVertex);
0383
0384 } catch (cms::Exception& e) {
0385
0386 edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
0387 << e.explainSelf();
0388 }
0389
0390
0391
0392
0393
0394
0395
0396
0397
0398
0399
0400
0401
0402
0403 trackPairRef.clear();
0404 trackInnPos.clear();
0405 trackPin.clear();
0406 trackPout.clear();
0407
0408 for (std::vector<reco::TransientTrack>::const_iterator iTk = (iPair->first).begin();
0409 iTk != (iPair->first).end();
0410 ++iTk) {
0411
0412
0413 auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
0414 reco::TrackRef myTkRef = ttt->persistentTrackRef();
0415
0416
0417 if (myTkRef->extra().isNonnull()) {
0418 trackInnPos.push_back(toFConverterP(myTkRef->innerPosition()));
0419 trackPin.push_back(toFConverterV(myTkRef->innerMomentum()));
0420 trackPout.push_back(toFConverterV(myTkRef->outerMomentum()));
0421 }
0422 trackPairRef.push_back(myTkRef);
0423 }
0424
0425
0426
0427
0428
0429 if (theConversionVertex.isValid()) {
0430
0431 }
0432
0433
0434 minAppDist = calculateMinApproachDistance(trackPairRef[0], trackPairRef[1]);
0435
0436 double like = -999.;
0437 reco::Conversion newCandidate(scPtrVec,
0438 trackPairRef,
0439 trkPositionAtEcal,
0440 theConversionVertex,
0441 matchingBC,
0442 minAppDist,
0443 trackInnPos,
0444 trackPin,
0445 trackPout,
0446 like,
0447 algo);
0448 like = likelihoodCalc_.calculateLikelihood(newCandidate);
0449
0450 newCandidate.setMVAout(like);
0451 outputConvPhotonCollection.push_back(newCandidate);
0452
0453
0454
0455 } else {
0456
0457
0458
0459 trackPairRef.clear();
0460 trackInnPos.clear();
0461 trackPin.clear();
0462 trackPout.clear();
0463 std::vector<reco::TransientTrack>::const_iterator iTk = (iPair->first).begin();
0464
0465 auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
0466 reco::TrackRef myTk = ttt->persistentTrackRef();
0467 if (myTk->extra().isNonnull()) {
0468 trackInnPos.push_back(toFConverterP(myTk->innerPosition()));
0469 trackPin.push_back(toFConverterV(myTk->innerMomentum()));
0470 trackPout.push_back(toFConverterV(myTk->outerMomentum()));
0471 }
0472 trackPairRef.push_back(myTk);
0473
0474
0475 if (recoverOneTrackCase_) {
0476 float theta1 = myTk->innerMomentum().Theta();
0477 float dCot = 999.;
0478 float dCotTheta = -999.;
0479 reco::TrackRef goodRef;
0480 for (auto const& tran : t_generalTrk) {
0481 auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(tran.basicTransientTrack());
0482 reco::TrackRef trRef = ttt->persistentTrackRef();
0483 if (trRef->charge() * myTk->charge() > 0)
0484 continue;
0485 float dEta = trRef->eta() - myTk->eta();
0486 float dPhi = trRef->phi() - myTk->phi();
0487 if (sqrt(dEta * dEta + dPhi * dPhi) > dRForConversionRecovery_)
0488 continue;
0489 float theta2 = trRef->innerMomentum().Theta();
0490 dCotTheta = 1. / tan(theta1) - 1. / tan(theta2);
0491
0492 if (fabs(dCotTheta) < dCot) {
0493 dCot = fabs(dCotTheta);
0494 goodRef = trRef;
0495 }
0496 }
0497
0498 if (goodRef.isNonnull()) {
0499 minAppDist = calculateMinApproachDistance(myTk, goodRef);
0500
0501
0502 if (fabs(dCotTheta) < deltaCotCut_ && minAppDist > minApproachDisCut_) {
0503 trackInnPos.push_back(toFConverterP(goodRef->innerPosition()));
0504 trackPin.push_back(toFConverterV(goodRef->innerMomentum()));
0505 trackPout.push_back(toFConverterV(goodRef->outerMomentum()));
0506 trackPairRef.push_back(goodRef);
0507
0508
0509
0510 try {
0511 vertexFinder_.run(iPair->first, theConversionVertex);
0512
0513 } catch (cms::Exception& e) {
0514
0515 edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
0516 << e.explainSelf();
0517 }
0518 }
0519 }
0520
0521 }
0522 const double like = -999.;
0523 outputConvPhotonCollection.emplace_back(scPtrVec,
0524 trackPairRef,
0525 trkPositionAtEcal,
0526 theConversionVertex,
0527 matchingBC,
0528 minAppDist,
0529 trackInnPos,
0530 trackPin,
0531 trackPout,
0532 like,
0533 algo);
0534 auto& newCandidate = outputConvPhotonCollection.back();
0535 newCandidate.setMVAout(likelihoodCalc_.calculateLikelihood(newCandidate));
0536
0537 }
0538 }
0539 }
0540 }
0541 }
0542
0543 void ConvertedPhotonProducer::cleanCollections(const edm::Handle<edm::View<reco::CaloCluster>>& scHandle,
0544 const edm::OrphanHandle<reco::ConversionCollection>& conversionHandle,
0545 reco::ConversionCollection& outputConversionCollection) {
0546 reco::Conversion* newCandidate = nullptr;
0547 for (auto const& aClus : scHandle->ptrs()) {
0548
0549 if (aClus->energy() / cosh(aClus->eta()) <= minSCEt_)
0550 continue;
0551
0552 if (conversionHandle.isValid()) {
0553 if (risolveAmbiguity_) {
0554 std::vector<reco::ConversionRef> bestRef = solveAmbiguity(conversionHandle, aClus);
0555
0556 for (std::vector<reco::ConversionRef>::iterator iRef = bestRef.begin(); iRef != bestRef.end(); iRef++) {
0557 if (iRef->isNonnull()) {
0558 newCandidate = (*iRef)->clone();
0559 outputConversionCollection.push_back(*newCandidate);
0560 delete newCandidate;
0561 }
0562 }
0563
0564 } else {
0565 for (unsigned int icp = 0; icp < conversionHandle->size(); icp++) {
0566 reco::ConversionRef cpRef(reco::ConversionRef(conversionHandle, icp));
0567 if (!(aClus.id() == cpRef->caloCluster()[0].id() && aClus.key() == cpRef->caloCluster()[0].key()))
0568 continue;
0569 if (!cpRef->isConverted())
0570 continue;
0571 if (cpRef->nTracks() < 2)
0572 continue;
0573 newCandidate = (&(*cpRef))->clone();
0574 outputConversionCollection.push_back(*newCandidate);
0575 delete newCandidate;
0576 }
0577
0578 }
0579 }
0580 }
0581 }
0582
0583 std::vector<reco::ConversionRef> ConvertedPhotonProducer::solveAmbiguity(
0584 const edm::OrphanHandle<reco::ConversionCollection>& conversionHandle, reco::CaloClusterPtr const& scRef) {
0585 std::multimap<double, reco::ConversionRef, std::greater<double>> convMap;
0586
0587 for (unsigned int icp = 0; icp < conversionHandle->size(); icp++) {
0588 reco::ConversionRef cpRef{conversionHandle, icp};
0589
0590
0591 if (!(scRef.id() == cpRef->caloCluster()[0].id() && scRef.key() == cpRef->caloCluster()[0].key()))
0592 continue;
0593 if (!cpRef->isConverted())
0594 continue;
0595 double like = cpRef->MVAout();
0596 if (cpRef->nTracks() < 2)
0597 continue;
0598
0599 convMap.emplace(like, cpRef);
0600 }
0601
0602
0603
0604 std::vector<reco::ConversionRef> bestRefs;
0605 for (auto iMap = convMap.begin(); iMap != convMap.end(); iMap++) {
0606
0607 bestRefs.push_back(iMap->second);
0608 if (int(bestRefs.size()) == maxNumOfCandidates_)
0609 break;
0610 }
0611
0612 return bestRefs;
0613 }
0614
0615 float ConvertedPhotonProducer::calculateMinApproachDistance(const reco::TrackRef& track1,
0616 const reco::TrackRef& track2) {
0617 double x1, x2, y1, y2;
0618 double xx_1 = track1->innerPosition().x(), yy_1 = track1->innerPosition().y(), zz_1 = track1->innerPosition().z();
0619 double xx_2 = track2->innerPosition().x(), yy_2 = track2->innerPosition().y(), zz_2 = track2->innerPosition().z();
0620 double radius1 =
0621 track1->innerMomentum().Rho() / (.3 * (magneticField_->inTesla(GlobalPoint(xx_1, yy_1, zz_1)).z())) * 100;
0622 double radius2 =
0623 track2->innerMomentum().Rho() / (.3 * (magneticField_->inTesla(GlobalPoint(xx_2, yy_2, zz_2)).z())) * 100;
0624 getCircleCenter(track1, radius1, x1, y1);
0625 getCircleCenter(track2, radius2, x2, y2);
0626
0627 return std::sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) - radius1 - radius2;
0628 }
0629
0630 void ConvertedPhotonProducer::getCircleCenter(const reco::TrackRef& tk, double r, double& x0, double& y0) {
0631 double x1, y1, phi;
0632 x1 = tk->innerPosition().x();
0633 y1 = tk->innerPosition().y();
0634 phi = tk->innerMomentum().phi();
0635 const int charge = tk->charge();
0636 x0 = x1 + r * sin(phi) * charge;
0637 y0 = y1 - r * cos(phi) * charge;
0638 }