Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-07-15 08:47:55

0001 /** \class ConvertedPhotonProducer
0002  **  
0003  **
0004  **  \author Nancy Marinelli, U. of Notre Dame, US
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   // Register the product
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   // instantiate the Track Pair Finder algorithm
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   // Transform Track into TransientTrack (needed by the Vertex fitter)
0183   transientTrackBuilder_ = &theEventSetup.getData(transientTrackToken_);
0184 }
0185 
0186 void ConvertedPhotonProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
0187   //
0188   // create empty output collections
0189   //
0190   // Converted photon candidates
0191   reco::ConversionCollection outputConvPhotonCollection;
0192   // Converted photon candidates
0193   reco::ConversionCollection cleanedConversionCollection;
0194 
0195   // Get the Super Cluster collection in the Barrel
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   // Get the Super Cluster collection in the Endcap
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   //// Get the Out In CKF tracks from conversions
0213   bool validTrackInputs = true;
0214   auto outInTrkHandle = theEvent.getHandle(conversionOITrackProducer_);
0215   if (!outInTrkHandle.isValid()) {
0216     //std::cout << "Error! Can't get the conversionOITrack " << "\n";
0217     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the conversionOITrack "
0218                                              << "\n";
0219     validTrackInputs = false;
0220   }
0221   //  LogDebug("ConvertedPhotonProducer")<< "ConvertedPhotonProducer  outInTrack collection size " << (*outInTrkHandle).size() << "\n";
0222 
0223   //// Get the association map between CKF Out In tracks and the SC where they originated
0224   auto outInTrkSCAssocHandle = theEvent.getHandle(outInTrackSCAssociationCollection_);
0225   if (!outInTrkSCAssocHandle.isValid()) {
0226     //  std::cout << "Error! Can't get the product " <<  outInTrackSCAssociationCollection_.c_str() <<"\n";
0227     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the outInTrackSCAssociationCollection)";
0228     validTrackInputs = false;
0229   }
0230 
0231   //// Get the In Out  CKF tracks from conversions
0232   auto inOutTrkHandle = theEvent.getHandle(conversionIOTrackProducer_);
0233   if (!inOutTrkHandle.isValid()) {
0234     // std::cout << "Error! Can't get the conversionIOTrack " << "\n";
0235     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the conversionIOTrack "
0236                                              << "\n";
0237     validTrackInputs = false;
0238   }
0239   //  LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer inOutTrack collection size " << (*inOutTrkHandle).size() << "\n";
0240 
0241   //// Get the generalTracks if the recovery of one track cases is switched on
0242 
0243   edm::Handle<reco::TrackCollection> generalTrkHandle;
0244   if (recoverOneTrackCase_) {
0245     theEvent.getByToken(generalTrackProducer_, generalTrkHandle);
0246     if (!generalTrkHandle.isValid()) {
0247       //std::cout << "Error! Can't get the genralTracks " << "\n";
0248       edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the genralTracks "
0249                                                << "\n";
0250     }
0251   }
0252 
0253   //// Get the association map between CKF in out tracks and the SC  where they originated
0254   auto inOutTrkSCAssocHandle = theEvent.getHandle(inOutTrackSCAssociationCollection_);
0255   if (!inOutTrkSCAssocHandle.isValid()) {
0256     //std::cout << "Error! Can't get the product " <<  inOutTrackSCAssociationCollection_.c_str() <<"\n";
0257     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the inOutTrackSCAssociationCollection_.c_str()";
0258     validTrackInputs = false;
0259   }
0260 
0261   // Get the basic cluster collection in the Barrel
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   // Get the basic cluster collection in the Endcap
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     //do the conversion:
0277     std::vector<reco::TransientTrack> t_outInTrk = transientTrackBuilder_->build(outInTrkHandle);
0278     std::vector<reco::TransientTrack> t_inOutTrk = transientTrackBuilder_->build(inOutTrkHandle);
0279 
0280     ///// Find the +/- pairs
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     //LogDebug("ConvertedPhotonProducer")  << "ConvertedPhotonProducer  allPairs.size " << allPairs.size() << "\n";
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   // put the product in the event
0305   auto const conversionHandle =
0306       theEvent.emplace(convertedPhotonCollectionPutToken_, std::move(outputConvPhotonCollection));
0307 
0308   // Loop over barrel and endcap SC collections and fill the  photon collection
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   // instantiate the algorithm for finding the position of the track extrapolation at the Ecal front face
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   //  Loop over SC in the barrel and reconstruct converted photons
0337   int myCands = 0;
0338   reco::CaloClusterPtrVector scPtrVec;
0339   for (auto const& aClus : scHandle->ptrs()) {
0340     // preselection based in Et and H/E cut
0341     if (aClus->energy() / cosh(aClus->eta()) <= minSCEt_)
0342       continue;
0343     const reco::CaloCluster* pClus = &(*aClus);
0344     auto const* sc = dynamic_cast<const reco::SuperCluster*>(pClus);
0345     double HoE = hcalHelper.hcalESum(*sc, 0) / sc->energy();
0346     if (HoE >= maxHOverE_)
0347       continue;
0348     /////
0349 
0350     std::vector<edm::Ref<reco::TrackCollection>> trackPairRef;
0351     std::vector<math::XYZPointF> trackInnPos;
0352     std::vector<math::XYZVectorF> trackPin;
0353     std::vector<math::XYZVectorF> trackPout;
0354     float minAppDist = -99;
0355 
0356     //LogDebug("ConvertedPhotonProducer") << "ConvertedPhotonProducer SC energy " << aClus->energy() << " eta " <<  aClus->eta() << " phi " <<  aClus->phi() << "\n";
0357 
0358     //// Set here first quantities for the converted photon
0359     const reco::Particle::Point vtx(0, 0, 0);
0360 
0361     math::XYZVector direction = aClus->position() - vtx;
0362     math::XYZVector momentum = direction.unit() * aClus->energy();
0363     const reco::Particle::LorentzVector p4(momentum.x(), momentum.y(), momentum.z(), aClus->energy());
0364 
0365     int nFound = 0;
0366     if (!allPairs.empty()) {
0367       nFound = 0;
0368 
0369       for (auto iPair = allPairs.begin(); iPair != allPairs.end(); ++iPair) {
0370         scPtrVec.clear();
0371 
0372         reco::Vertex theConversionVertex;
0373         reco::CaloClusterPtr caloPtr = iPair->second;
0374         if (!(aClus == caloPtr))
0375           continue;
0376 
0377         scPtrVec.push_back(aClus);
0378         nFound++;
0379 
0380         std::vector<math::XYZPointF> trkPositionAtEcal = theEcalImpactPositionFinder.find(iPair->first, bcHandle);
0381         std::vector<reco::CaloClusterPtr> matchingBC = theEcalImpactPositionFinder.matchingBC();
0382 
0383         minAppDist = -99;
0384         const std::string metname = "ConvertedPhotons|ConvertedPhotonProducer";
0385         if ((iPair->first).size() > 1) {
0386           try {
0387             vertexFinder_.run(iPair->first, theConversionVertex);
0388 
0389           } catch (cms::Exception& e) {
0390             //std::cout << " cms::Exception caught in ConvertedPhotonProducer::produce" << "\n" ;
0391             edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
0392                                      << e.explainSelf();
0393           }
0394 
0395           // Old TwoTrackMinimumDistance md;
0396           // Old md.calculate  (  (iPair->first)[0].initialFreeState(),  (iPair->first)[1].initialFreeState() );
0397           // Old minAppDist = md.distance();
0398 
0399           /*
0400     for ( unsigned int i=0; i< matchingBC.size(); ++i) {
0401           if (  matchingBC[i].isNull() )  std::cout << " This ref to BC is null: skipping " <<  "\n";
0402           else 
0403         std::cout << " BC energy " << matchingBC[i]->energy() <<  "\n";
0404     }
0405     */
0406 
0407           //// loop over tracks in the pair  for creating a reference
0408           trackPairRef.clear();
0409           trackInnPos.clear();
0410           trackPin.clear();
0411           trackPout.clear();
0412 
0413           for (std::vector<reco::TransientTrack>::const_iterator iTk = (iPair->first).begin();
0414                iTk != (iPair->first).end();
0415                ++iTk) {
0416             //LogDebug("ConvertedPhotonProducer")  << "  ConvertedPhotonProducer Transient Tracks in the pair  charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";
0417 
0418             auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
0419             reco::TrackRef myTkRef = ttt->persistentTrackRef();
0420 
0421             //LogDebug("ConvertedPhotonProducer")  << " ConvertedPhotonProducer Ref to Rec Tracks in the pair  charge " << myTkRef->charge() << " Num of RecHits " << myTkRef->recHitsSize() << " inner momentum " << myTkRef->innerMomentum() << "\n";
0422             if (myTkRef->extra().isNonnull()) {
0423               trackInnPos.push_back(toFConverterP(myTkRef->innerPosition()));
0424               trackPin.push_back(toFConverterV(myTkRef->innerMomentum()));
0425               trackPout.push_back(toFConverterV(myTkRef->outerMomentum()));
0426             }
0427             trackPairRef.push_back(myTkRef);
0428           }
0429 
0430           //    std::cout << " ConvertedPhotonProducer trackPin size " << trackPin.size() << std::endl;
0431           //LogDebug("ConvertedPhotonProducer")  << " ConvertedPhotonProducer SC energy " <<  aClus->energy() << "\n";
0432           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer photon p4 " << p4  << "\n";
0433           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer vtx " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
0434           if (theConversionVertex.isValid()) {
0435             //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer theConversionVertex " << theConversionVertex.position().x() << " " << theConversionVertex.position().y() << " " << theConversionVertex.position().z() << "\n";
0436           }
0437           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer trackPairRef  " << trackPairRef.size() <<  "\n";
0438 
0439           minAppDist = calculateMinApproachDistance(trackPairRef[0], trackPairRef[1]);
0440 
0441           double like = -999.;
0442           reco::Conversion newCandidate(scPtrVec,
0443                                         trackPairRef,
0444                                         trkPositionAtEcal,
0445                                         theConversionVertex,
0446                                         matchingBC,
0447                                         minAppDist,
0448                                         trackInnPos,
0449                                         trackPin,
0450                                         trackPout,
0451                                         like,
0452                                         algo);
0453           like = likelihoodCalc_.calculateLikelihood(newCandidate);
0454           //    std::cout << "like = " << like << std::endl;
0455           newCandidate.setMVAout(like);
0456           outputConvPhotonCollection.push_back(newCandidate);
0457 
0458           myCands++;
0459           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer Put the ConvertedPhotonCollection a candidate in the Barrel " << "\n";
0460 
0461         } else {
0462           //      std::cout << "   ConvertedPhotonProducer case with only one track found " <<  "\n";
0463 
0464           //std::cout << "   ConvertedPhotonProducer recovering one track " <<  "\n";
0465           trackPairRef.clear();
0466           trackInnPos.clear();
0467           trackPin.clear();
0468           trackPout.clear();
0469           std::vector<reco::TransientTrack>::const_iterator iTk = (iPair->first).begin();
0470           //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";
0471           auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
0472           reco::TrackRef myTk = ttt->persistentTrackRef();
0473           if (myTk->extra().isNonnull()) {
0474             trackInnPos.push_back(toFConverterP(myTk->innerPosition()));
0475             trackPin.push_back(toFConverterV(myTk->innerMomentum()));
0476             trackPout.push_back(toFConverterV(myTk->outerMomentum()));
0477           }
0478           trackPairRef.push_back(myTk);
0479           //std::cout << " Provenance " << myTk->algoName() << std::endl;
0480 
0481           if (recoverOneTrackCase_) {
0482             float theta1 = myTk->innerMomentum().Theta();
0483             float dCot = 999.;
0484             float dCotTheta = -999.;
0485             reco::TrackRef goodRef;
0486             for (auto const& tran : t_generalTrk) {
0487               auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(tran.basicTransientTrack());
0488               reco::TrackRef trRef = ttt->persistentTrackRef();
0489               if (trRef->charge() * myTk->charge() > 0)
0490                 continue;
0491               float dEta = trRef->eta() - myTk->eta();
0492               float dPhi = trRef->phi() - myTk->phi();
0493               if (sqrt(dEta * dEta + dPhi * dPhi) > dRForConversionRecovery_)
0494                 continue;
0495               float theta2 = trRef->innerMomentum().Theta();
0496               dCotTheta = 1. / tan(theta1) - 1. / tan(theta2);
0497               //    std::cout << "  ConvertedPhotonProducer recovering general transient track charge " << trRef->charge() << " momentum " << trRef->innerMomentum() << " dcotTheta " << fabs(dCotTheta) << std::endl;
0498               if (fabs(dCotTheta) < dCot) {
0499                 dCot = fabs(dCotTheta);
0500                 goodRef = trRef;
0501               }
0502             }
0503 
0504             if (goodRef.isNonnull()) {
0505               minAppDist = calculateMinApproachDistance(myTk, goodRef);
0506 
0507               // std::cout << "  ConvertedPhotonProducer chosen dCotTheta " <<  fabs(dCotTheta) << std::endl;
0508               if (fabs(dCotTheta) < deltaCotCut_ && minAppDist > minApproachDisCut_) {
0509                 trackInnPos.push_back(toFConverterP(goodRef->innerPosition()));
0510                 trackPin.push_back(toFConverterV(goodRef->innerMomentum()));
0511                 trackPout.push_back(toFConverterV(goodRef->outerMomentum()));
0512                 trackPairRef.push_back(goodRef);
0513                 //      std::cout << " ConvertedPhotonProducer adding opposite charge track from generalTrackCollection charge " <<  goodRef ->charge() << " pt " << sqrt(goodRef->innerMomentum().perp2())  << " trackPairRef size " << trackPairRef.size() << std::endl;
0514                 //std::cout << " Track Provenenance " << goodRef->algoName() << std::endl;
0515 
0516                 try {
0517                   vertexFinder_.run(iPair->first, theConversionVertex);
0518 
0519                 } catch (cms::Exception& e) {
0520                   //std::cout << " cms::Exception caught in ConvertedPhotonProducer::produce" << "\n" ;
0521                   edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
0522                                            << e.explainSelf();
0523                 }
0524               }
0525             }
0526 
0527           }  // bool On/Off one track case recovery using generalTracks
0528           const double like = -999.;
0529           outputConvPhotonCollection.emplace_back(scPtrVec,
0530                                                   trackPairRef,
0531                                                   trkPositionAtEcal,
0532                                                   theConversionVertex,
0533                                                   matchingBC,
0534                                                   minAppDist,
0535                                                   trackInnPos,
0536                                                   trackPin,
0537                                                   trackPout,
0538                                                   like,
0539                                                   algo);
0540           auto& newCandidate = outputConvPhotonCollection.back();
0541           newCandidate.setMVAout(likelihoodCalc_.calculateLikelihood(newCandidate));
0542 
0543         }  // case with only on track: looking in general tracks
0544       }
0545     }
0546   }
0547 }
0548 
0549 void ConvertedPhotonProducer::cleanCollections(const edm::Handle<edm::View<reco::CaloCluster>>& scHandle,
0550                                                const edm::OrphanHandle<reco::ConversionCollection>& conversionHandle,
0551                                                reco::ConversionCollection& outputConversionCollection) {
0552   reco::Conversion* newCandidate = nullptr;
0553   for (auto const& aClus : scHandle->ptrs()) {
0554     // SC energy preselection
0555     if (aClus->energy() / cosh(aClus->eta()) <= minSCEt_)
0556       continue;
0557 
0558     if (conversionHandle.isValid()) {
0559       if (risolveAmbiguity_) {
0560         std::vector<reco::ConversionRef> bestRef = solveAmbiguity(conversionHandle, aClus);
0561 
0562         for (std::vector<reco::ConversionRef>::iterator iRef = bestRef.begin(); iRef != bestRef.end(); iRef++) {
0563           if (iRef->isNonnull()) {
0564             newCandidate = (*iRef)->clone();
0565             outputConversionCollection.push_back(*newCandidate);
0566             delete newCandidate;
0567           }
0568         }
0569 
0570       } else {
0571         for (unsigned int icp = 0; icp < conversionHandle->size(); icp++) {
0572           reco::ConversionRef cpRef(reco::ConversionRef(conversionHandle, icp));
0573           if (!(aClus.id() == cpRef->caloCluster()[0].id() && aClus.key() == cpRef->caloCluster()[0].key()))
0574             continue;
0575           if (!cpRef->isConverted())
0576             continue;
0577           if (cpRef->nTracks() < 2)
0578             continue;
0579           newCandidate = (&(*cpRef))->clone();
0580           outputConversionCollection.push_back(*newCandidate);
0581           delete newCandidate;
0582         }
0583 
0584       }  // solve or not the ambiguity of many conversion candidates
0585     }
0586   }
0587 }
0588 
0589 std::vector<reco::ConversionRef> ConvertedPhotonProducer::solveAmbiguity(
0590     const edm::OrphanHandle<reco::ConversionCollection>& conversionHandle, reco::CaloClusterPtr const& scRef) {
0591   std::multimap<double, reco::ConversionRef, std::greater<double>> convMap;
0592 
0593   for (unsigned int icp = 0; icp < conversionHandle->size(); icp++) {
0594     reco::ConversionRef cpRef{conversionHandle, icp};
0595 
0596     //std::cout << " cpRef " << cpRef->nTracks() << " " <<  cpRef ->caloCluster()[0]->energy() << std::endl;
0597     if (!(scRef.id() == cpRef->caloCluster()[0].id() && scRef.key() == cpRef->caloCluster()[0].key()))
0598       continue;
0599     if (!cpRef->isConverted())
0600       continue;
0601     double like = cpRef->MVAout();
0602     if (cpRef->nTracks() < 2)
0603       continue;
0604     //    std::cout << " Like " << like << std::endl;
0605     convMap.emplace(like, cpRef);
0606   }
0607 
0608   //  std::cout << " convMap size " << convMap.size() << std::endl;
0609 
0610   std::vector<reco::ConversionRef> bestRefs;
0611   for (auto iMap = convMap.begin(); iMap != convMap.end(); iMap++) {
0612     //    std::cout << " Like list in the map " <<  iMap->first << " " << (iMap->second)->EoverP() << std::endl;
0613     bestRefs.push_back(iMap->second);
0614     if (int(bestRefs.size()) == maxNumOfCandidates_)
0615       break;
0616   }
0617 
0618   return bestRefs;
0619 }
0620 
0621 float ConvertedPhotonProducer::calculateMinApproachDistance(const reco::TrackRef& track1,
0622                                                             const reco::TrackRef& track2) {
0623   double x1, x2, y1, y2;
0624   double xx_1 = track1->innerPosition().x(), yy_1 = track1->innerPosition().y(), zz_1 = track1->innerPosition().z();
0625   double xx_2 = track2->innerPosition().x(), yy_2 = track2->innerPosition().y(), zz_2 = track2->innerPosition().z();
0626   double radius1 =
0627       track1->innerMomentum().Rho() / (.3 * (magneticField_->inTesla(GlobalPoint(xx_1, yy_1, zz_1)).z())) * 100;
0628   double radius2 =
0629       track2->innerMomentum().Rho() / (.3 * (magneticField_->inTesla(GlobalPoint(xx_2, yy_2, zz_2)).z())) * 100;
0630   getCircleCenter(track1, radius1, x1, y1);
0631   getCircleCenter(track2, radius2, x2, y2);
0632 
0633   return std::sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) - radius1 - radius2;
0634 }
0635 
0636 void ConvertedPhotonProducer::getCircleCenter(const reco::TrackRef& tk, double r, double& x0, double& y0) {
0637   double x1, y1, phi;
0638   x1 = tk->innerPosition().x();  //inner position and inner momentum need track Extra!
0639   y1 = tk->innerPosition().y();
0640   phi = tk->innerMomentum().phi();
0641   const int charge = tk->charge();
0642   x0 = x1 + r * sin(phi) * charge;
0643   y0 = y1 - r * cos(phi) * charge;
0644 }