Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-05-08 23:19:13 

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   reco::CaloClusterPtrVector scPtrVec;
0338   for (auto const& aClus : scHandle->ptrs()) {
0339     // preselection based in Et and H/E cut
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     //LogDebug("ConvertedPhotonProducer") << "ConvertedPhotonProducer SC energy " << aClus->energy() << " eta " <<  aClus->eta() << " phi " <<  aClus->phi() << "\n";
0356 
0357     //// Set here first quantities for the converted photon
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             //std::cout << " cms::Exception caught in ConvertedPhotonProducer::produce" << "\n" ;
0386             edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
0387                                      << e.explainSelf();
0388           }
0389 
0390           // Old TwoTrackMinimumDistance md;
0391           // Old md.calculate  (  (iPair->first)[0].initialFreeState(),  (iPair->first)[1].initialFreeState() );
0392           // Old minAppDist = md.distance();
0393 
0394           /*
0395     for ( unsigned int i=0; i< matchingBC.size(); ++i) {
0396           if (  matchingBC[i].isNull() )  std::cout << " This ref to BC is null: skipping " <<  "\n";
0397           else 
0398         std::cout << " BC energy " << matchingBC[i]->energy() <<  "\n";
0399     }
0400     */
0401 
0402           //// loop over tracks in the pair  for creating a reference
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             //LogDebug("ConvertedPhotonProducer")  << "  ConvertedPhotonProducer Transient Tracks in the pair  charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";
0412 
0413             auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
0414             reco::TrackRef myTkRef = ttt->persistentTrackRef();
0415 
0416             //LogDebug("ConvertedPhotonProducer")  << " ConvertedPhotonProducer Ref to Rec Tracks in the pair  charge " << myTkRef->charge() << " Num of RecHits " << myTkRef->recHitsSize() << " inner momentum " << myTkRef->innerMomentum() << "\n";
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           //    std::cout << " ConvertedPhotonProducer trackPin size " << trackPin.size() << std::endl;
0426           //LogDebug("ConvertedPhotonProducer")  << " ConvertedPhotonProducer SC energy " <<  aClus->energy() << "\n";
0427           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer photon p4 " << p4  << "\n";
0428           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer vtx " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
0429           if (theConversionVertex.isValid()) {
0430             //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer theConversionVertex " << theConversionVertex.position().x() << " " << theConversionVertex.position().y() << " " << theConversionVertex.position().z() << "\n";
0431           }
0432           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer trackPairRef  " << trackPairRef.size() <<  "\n";
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           //    std::cout << "like = " << like << std::endl;
0450           newCandidate.setMVAout(like);
0451           outputConvPhotonCollection.push_back(newCandidate);
0452 
0453           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer Put the ConvertedPhotonCollection a candidate in the Barrel " << "\n";
0454 
0455         } else {
0456           //      std::cout << "   ConvertedPhotonProducer case with only one track found " <<  "\n";
0457 
0458           //std::cout << "   ConvertedPhotonProducer recovering one track " <<  "\n";
0459           trackPairRef.clear();
0460           trackInnPos.clear();
0461           trackPin.clear();
0462           trackPout.clear();
0463           std::vector<reco::TransientTrack>::const_iterator iTk = (iPair->first).begin();
0464           //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";
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           //std::cout << " Provenance " << myTk->algoName() << std::endl;
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               //    std::cout << "  ConvertedPhotonProducer recovering general transient track charge " << trRef->charge() << " momentum " << trRef->innerMomentum() << " dcotTheta " << fabs(dCotTheta) << std::endl;
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               // std::cout << "  ConvertedPhotonProducer chosen dCotTheta " <<  fabs(dCotTheta) << std::endl;
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                 //      std::cout << " ConvertedPhotonProducer adding opposite charge track from generalTrackCollection charge " <<  goodRef ->charge() << " pt " << sqrt(goodRef->innerMomentum().perp2())  << " trackPairRef size " << trackPairRef.size() << std::endl;
0508                 //std::cout << " Track Provenenance " << goodRef->algoName() << std::endl;
0509 
0510                 try {
0511                   vertexFinder_.run(iPair->first, theConversionVertex);
0512 
0513                 } catch (cms::Exception& e) {
0514                   //std::cout << " cms::Exception caught in ConvertedPhotonProducer::produce" << "\n" ;
0515                   edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
0516                                            << e.explainSelf();
0517                 }
0518               }
0519             }
0520 
0521           }  // bool On/Off one track case recovery using generalTracks
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         }  // case with only on track: looking in general tracks
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     // SC energy preselection
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       }  // solve or not the ambiguity of many conversion candidates
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     //std::cout << " cpRef " << cpRef->nTracks() << " " <<  cpRef ->caloCluster()[0]->energy() << std::endl;
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     //    std::cout << " Like " << like << std::endl;
0599     convMap.emplace(like, cpRef);
0600   }
0601 
0602   //  std::cout << " convMap size " << convMap.size() << std::endl;
0603 
0604   std::vector<reco::ConversionRef> bestRefs;
0605   for (auto iMap = convMap.begin(); iMap != convMap.end(); iMap++) {
0606     //    std::cout << " Like list in the map " <<  iMap->first << " " << (iMap->second)->EoverP() << std::endl;
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();  //inner position and inner momentum need track Extra!
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 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!