Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:25:00 

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 #include "CondFormats/DataRecord/interface/HcalPFCutsRcd.h"
0042 #include "CondTools/Hcal/interface/HcalPFCutsHandler.h"
0043 #include "Geometry/CaloTopology/interface/HcalTopology.h"
0044 
0045 #include <vector>
0046 
0047 class ConvertedPhotonProducer : public edm::stream::EDProducer<> {
0048 public:
0049   ConvertedPhotonProducer(const edm::ParameterSet& ps);
0050 
0051   void produce(edm::Event& evt, const edm::EventSetup& es) override;
0052 
0053 private:
0054   void buildCollections(
0055       edm::EventSetup const& es,
0056       const edm::Handle<edm::View<reco::CaloCluster>>& scHandle,
0057       const edm::Handle<edm::View<reco::CaloCluster>>& bcHandle,
0058       ElectronHcalHelper const& hcalHelper,
0059       const edm::Handle<reco::TrackCollection>& trkHandle,
0060       std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors>& allPairs,
0061       reco::ConversionCollection& outputConvPhotonCollection);
0062   void cleanCollections(const edm::Handle<edm::View<reco::CaloCluster>>& scHandle,
0063                         const edm::OrphanHandle<reco::ConversionCollection>& conversionHandle,
0064                         reco::ConversionCollection& outputCollection);
0065 
0066   std::vector<reco::ConversionRef> solveAmbiguity(const edm::OrphanHandle<reco::ConversionCollection>& conversionHandle,
0067                                                   reco::CaloClusterPtr const& sc);
0068 
0069   float calculateMinApproachDistance(const reco::TrackRef& track1, const reco::TrackRef& track2);
0070   void getCircleCenter(const reco::TrackRef& tk, double r, double& x0, double& y0);
0071 
0072   edm::EDGetTokenT<reco::TrackCollection> conversionOITrackProducer_;
0073   edm::EDGetTokenT<reco::TrackCollection> conversionIOTrackProducer_;
0074 
0075   edm::EDGetTokenT<reco::TrackCaloClusterPtrAssociation> outInTrackSCAssociationCollection_;
0076   edm::EDGetTokenT<reco::TrackCaloClusterPtrAssociation> inOutTrackSCAssociationCollection_;
0077 
0078   edm::EDGetTokenT<reco::TrackCollection> generalTrackProducer_;
0079 
0080   edm::ESGetToken<HcalPFCuts, HcalPFCutsRcd> hcalCutsToken_;
0081   bool cutsFromDB_;
0082   HcalPFCuts const* hcalCuts_ = nullptr;
0083 
0084   // Register the product
0085   edm::EDPutTokenT<reco::ConversionCollection> convertedPhotonCollectionPutToken_;
0086   edm::EDPutTokenT<reco::ConversionCollection> cleanedConvertedPhotonCollectionPutToken_;
0087 
0088   edm::EDGetTokenT<edm::View<reco::CaloCluster>> bcBarrelCollection_;
0089   edm::EDGetTokenT<edm::View<reco::CaloCluster>> bcEndcapCollection_;
0090   edm::EDGetTokenT<edm::View<reco::CaloCluster>> scHybridBarrelProducer_;
0091   edm::EDGetTokenT<edm::View<reco::CaloCluster>> scIslandEndcapProducer_;
0092   edm::EDGetTokenT<HBHERecHitCollection> hbheRecHits_;
0093 
0094   MagneticField const* magneticField_;
0095   TransientTrackBuilder const* transientTrackBuilder_;
0096 
0097   edm::ESGetToken<CaloGeometry, CaloGeometryRecord> caloGeomToken_;
0098   edm::ESGetToken<MagneticField, IdealMagneticFieldRecord> mFToken_;
0099   edm::ESGetToken<TransientTrackBuilder, TransientTrackRecord> transientTrackToken_;
0100 
0101   ConversionTrackPairFinder trackPairFinder_;
0102   ConversionVertexFinder vertexFinder_;
0103   std::string algoName_;
0104 
0105   double hOverEConeSize_;
0106   double maxHOverE_;
0107   double minSCEt_;
0108   bool recoverOneTrackCase_;
0109   double dRForConversionRecovery_;
0110   double deltaCotCut_;
0111   double minApproachDisCut_;
0112   int maxNumOfCandidates_;
0113   bool risolveAmbiguity_;
0114 
0115   std::unique_ptr<ElectronHcalHelper> hcalHelper_;
0116 
0117   ConversionLikelihoodCalculator likelihoodCalc_;
0118   std::string likelihoodWeights_;
0119 
0120   math::XYZPointF toFConverterP(const math::XYZPoint& val) { return math::XYZPointF(val.x(), val.y(), val.z()); }
0121 
0122   math::XYZVectorF toFConverterV(const math::XYZVector& val) { return math::XYZVectorF(val.x(), val.y(), val.z()); }
0123 };
0124 
0125 #include "FWCore/Framework/interface/MakerMacros.h"
0126 DEFINE_FWK_MODULE(ConvertedPhotonProducer);
0127 
0128 ConvertedPhotonProducer::ConvertedPhotonProducer(const edm::ParameterSet& config)
0129     : conversionOITrackProducer_{consumes(config.getParameter<std::string>("conversionOITrackProducer"))},
0130       conversionIOTrackProducer_{consumes(config.getParameter<std::string>("conversionIOTrackProducer"))},
0131       outInTrackSCAssociationCollection_{consumes({config.getParameter<std::string>("conversionOITrackProducer"),
0132                                                    config.getParameter<std::string>("outInTrackSCAssociation")})},
0133       inOutTrackSCAssociationCollection_{consumes({config.getParameter<std::string>("conversionIOTrackProducer"),
0134                                                    config.getParameter<std::string>("inOutTrackSCAssociation")})},
0135 
0136       generalTrackProducer_{consumes(config.getParameter<edm::InputTag>("generalTracksSrc"))},
0137       convertedPhotonCollectionPutToken_{
0138           produces<reco::ConversionCollection>(config.getParameter<std::string>("convertedPhotonCollection"))},
0139       cleanedConvertedPhotonCollectionPutToken_{
0140           produces<reco::ConversionCollection>(config.getParameter<std::string>("cleanedConvertedPhotonCollection"))},
0141 
0142       bcBarrelCollection_{consumes(config.getParameter<edm::InputTag>("bcBarrelCollection"))},
0143       bcEndcapCollection_{consumes(config.getParameter<edm::InputTag>("bcEndcapCollection"))},
0144       scHybridBarrelProducer_{consumes(config.getParameter<edm::InputTag>("scHybridBarrelProducer"))},
0145       scIslandEndcapProducer_{consumes(config.getParameter<edm::InputTag>("scIslandEndcapProducer"))},
0146       hbheRecHits_{consumes(config.getParameter<edm::InputTag>("hbheRecHits"))},
0147       caloGeomToken_{esConsumes()},
0148       mFToken_{esConsumes()},
0149       transientTrackToken_{esConsumes(edm::ESInputTag("", "TransientTrackBuilder"))},
0150       vertexFinder_{config},
0151       algoName_{config.getParameter<std::string>("AlgorithmName")},
0152 
0153       hOverEConeSize_{config.getParameter<double>("hOverEConeSize")},
0154       maxHOverE_{config.getParameter<double>("maxHOverE")},
0155       minSCEt_{config.getParameter<double>("minSCEt")},
0156       recoverOneTrackCase_{config.getParameter<bool>("recoverOneTrackCase")},
0157       dRForConversionRecovery_{config.getParameter<double>("dRForConversionRecovery")},
0158       deltaCotCut_{config.getParameter<double>("deltaCotCut")},
0159       minApproachDisCut_{config.getParameter<double>("minApproachDisCut")},
0160 
0161       maxNumOfCandidates_{config.getParameter<int>("maxNumOfCandidates")},
0162       risolveAmbiguity_{config.getParameter<bool>("risolveConversionAmbiguity")},
0163       likelihoodWeights_{config.getParameter<std::string>("MVA_weights_location")} {
0164   // instantiate the Track Pair Finder algorithm
0165   likelihoodCalc_.setWeightsFile(edm::FileInPath{likelihoodWeights_.c_str()}.fullPath().c_str());
0166 
0167   cutsFromDB_ = config.getParameter<bool>("usePFThresholdsFromDB");
0168   if (cutsFromDB_) {
0169     hcalCutsToken_ = esConsumes<HcalPFCuts, HcalPFCutsRcd>(edm::ESInputTag("", "withTopo"));
0170   }
0171 
0172   ElectronHcalHelper::Configuration cfgCone;
0173   cfgCone.hOverEConeSize = hOverEConeSize_;
0174   if (cfgCone.hOverEConeSize > 0) {
0175     cfgCone.onlyBehindCluster = false;
0176     cfgCone.checkHcalStatus = false;
0177 
0178     cfgCone.hbheRecHits = hbheRecHits_;
0179 
0180     cfgCone.eThresHB = config.getParameter<EgammaHcalIsolation::arrayHB>("recHitEThresholdHB");
0181     cfgCone.maxSeverityHB = config.getParameter<int>("maxHcalRecHitSeverity");
0182     cfgCone.eThresHE = config.getParameter<EgammaHcalIsolation::arrayHE>("recHitEThresholdHE");
0183     cfgCone.maxSeverityHE = cfgCone.maxSeverityHB;
0184   }
0185 
0186   hcalHelper_ = std::make_unique<ElectronHcalHelper>(cfgCone, consumesCollector());
0187 }
0188 
0189 void ConvertedPhotonProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
0190   magneticField_ = &theEventSetup.getData(mFToken_);
0191 
0192   // Transform Track into TransientTrack (needed by the Vertex fitter)
0193   transientTrackBuilder_ = &theEventSetup.getData(transientTrackToken_);
0194 
0195   if (cutsFromDB_) {
0196     hcalCuts_ = &theEventSetup.getData(hcalCutsToken_);
0197   }
0198 
0199   //
0200   // create empty output collections
0201   //
0202   // Converted photon candidates
0203   reco::ConversionCollection outputConvPhotonCollection;
0204   // Converted photon candidates
0205   reco::ConversionCollection cleanedConversionCollection;
0206 
0207   // Get the Super Cluster collection in the Barrel
0208   bool validBarrelSCHandle = true;
0209   auto scBarrelHandle = theEvent.getHandle(scHybridBarrelProducer_);
0210   if (!scBarrelHandle.isValid()) {
0211     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the scHybridBarrelProducer";
0212     validBarrelSCHandle = false;
0213   }
0214 
0215   // Get the Super Cluster collection in the Endcap
0216   bool validEndcapSCHandle = true;
0217   edm::Handle<edm::View<reco::CaloCluster>> scEndcapHandle;
0218   theEvent.getByToken(scIslandEndcapProducer_, scEndcapHandle);
0219   if (!scEndcapHandle.isValid()) {
0220     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the scIslandEndcapProducer";
0221     validEndcapSCHandle = false;
0222   }
0223 
0224   //// Get the Out In CKF tracks from conversions
0225   bool validTrackInputs = true;
0226   auto outInTrkHandle = theEvent.getHandle(conversionOITrackProducer_);
0227   if (!outInTrkHandle.isValid()) {
0228     //std::cout << "Error! Can't get the conversionOITrack " << "\n";
0229     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the conversionOITrack "
0230                                              << "\n";
0231     validTrackInputs = false;
0232   }
0233   //  LogDebug("ConvertedPhotonProducer")<< "ConvertedPhotonProducer  outInTrack collection size " << (*outInTrkHandle).size() << "\n";
0234 
0235   //// Get the association map between CKF Out In tracks and the SC where they originated
0236   auto outInTrkSCAssocHandle = theEvent.getHandle(outInTrackSCAssociationCollection_);
0237   if (!outInTrkSCAssocHandle.isValid()) {
0238     //  std::cout << "Error! Can't get the product " <<  outInTrackSCAssociationCollection_.c_str() <<"\n";
0239     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the outInTrackSCAssociationCollection)";
0240     validTrackInputs = false;
0241   }
0242 
0243   //// Get the In Out  CKF tracks from conversions
0244   auto inOutTrkHandle = theEvent.getHandle(conversionIOTrackProducer_);
0245   if (!inOutTrkHandle.isValid()) {
0246     // std::cout << "Error! Can't get the conversionIOTrack " << "\n";
0247     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the conversionIOTrack "
0248                                              << "\n";
0249     validTrackInputs = false;
0250   }
0251   //  LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer inOutTrack collection size " << (*inOutTrkHandle).size() << "\n";
0252 
0253   //// Get the generalTracks if the recovery of one track cases is switched on
0254 
0255   edm::Handle<reco::TrackCollection> generalTrkHandle;
0256   if (recoverOneTrackCase_) {
0257     theEvent.getByToken(generalTrackProducer_, generalTrkHandle);
0258     if (!generalTrkHandle.isValid()) {
0259       //std::cout << "Error! Can't get the genralTracks " << "\n";
0260       edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the genralTracks "
0261                                                << "\n";
0262     }
0263   }
0264 
0265   //// Get the association map between CKF in out tracks and the SC  where they originated
0266   auto inOutTrkSCAssocHandle = theEvent.getHandle(inOutTrackSCAssociationCollection_);
0267   if (!inOutTrkSCAssocHandle.isValid()) {
0268     //std::cout << "Error! Can't get the product " <<  inOutTrackSCAssociationCollection_.c_str() <<"\n";
0269     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the inOutTrackSCAssociationCollection_.c_str()";
0270     validTrackInputs = false;
0271   }
0272 
0273   // Get the basic cluster collection in the Barrel
0274   edm::Handle<edm::View<reco::CaloCluster>> bcBarrelHandle;
0275   theEvent.getByToken(bcBarrelCollection_, bcBarrelHandle);
0276   if (!bcBarrelHandle.isValid()) {
0277     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the bcBarrelCollection";
0278   }
0279 
0280   // Get the basic cluster collection in the Endcap
0281   edm::Handle<edm::View<reco::CaloCluster>> bcEndcapHandle;
0282   theEvent.getByToken(bcEndcapCollection_, bcEndcapHandle);
0283   if (!bcEndcapHandle.isValid()) {
0284     edm::LogError("ConvertedPhotonProducer") << "Error! Can't get the bcEndcapCollection";
0285   }
0286 
0287   if (validTrackInputs) {
0288     //do the conversion:
0289     std::vector<reco::TransientTrack> t_outInTrk = transientTrackBuilder_->build(outInTrkHandle);
0290     std::vector<reco::TransientTrack> t_inOutTrk = transientTrackBuilder_->build(inOutTrkHandle);
0291 
0292     ///// Find the +/- pairs
0293     std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors> allPairs;
0294     allPairs = trackPairFinder_.run(
0295         t_outInTrk, outInTrkHandle, outInTrkSCAssocHandle, t_inOutTrk, inOutTrkHandle, inOutTrkSCAssocHandle);
0296     //LogDebug("ConvertedPhotonProducer")  << "ConvertedPhotonProducer  allPairs.size " << allPairs.size() << "\n";
0297 
0298     hcalHelper_->beginEvent(theEvent, theEventSetup);
0299 
0300     buildCollections(theEventSetup,
0301                      scBarrelHandle,
0302                      bcBarrelHandle,
0303                      *hcalHelper_,
0304                      generalTrkHandle,
0305                      allPairs,
0306                      outputConvPhotonCollection);
0307     buildCollections(theEventSetup,
0308                      scEndcapHandle,
0309                      bcEndcapHandle,
0310                      *hcalHelper_,
0311                      generalTrkHandle,
0312                      allPairs,
0313                      outputConvPhotonCollection);
0314   }
0315 
0316   // put the product in the event
0317   auto const conversionHandle =
0318       theEvent.emplace(convertedPhotonCollectionPutToken_, std::move(outputConvPhotonCollection));
0319 
0320   // Loop over barrel and endcap SC collections and fill the  photon collection
0321   if (validBarrelSCHandle)
0322     cleanCollections(scBarrelHandle, conversionHandle, cleanedConversionCollection);
0323   if (validEndcapSCHandle)
0324     cleanCollections(scEndcapHandle, conversionHandle, cleanedConversionCollection);
0325 
0326   theEvent.emplace(cleanedConvertedPhotonCollectionPutToken_, std::move(cleanedConversionCollection));
0327 }
0328 
0329 void ConvertedPhotonProducer::buildCollections(
0330     edm::EventSetup const& es,
0331     const edm::Handle<edm::View<reco::CaloCluster>>& scHandle,
0332     const edm::Handle<edm::View<reco::CaloCluster>>& bcHandle,
0333     ElectronHcalHelper const& hcalHelper,
0334     const edm::Handle<reco::TrackCollection>& generalTrkHandle,
0335     std::map<std::vector<reco::TransientTrack>, reco::CaloClusterPtr, CompareTwoTracksVectors>& allPairs,
0336     reco::ConversionCollection& outputConvPhotonCollection)
0337 
0338 {
0339   // instantiate the algorithm for finding the position of the track extrapolation at the Ecal front face
0340   ConversionTrackEcalImpactPoint theEcalImpactPositionFinder(magneticField_);
0341 
0342   reco::Conversion::ConversionAlgorithm algo = reco::Conversion::algoByName(algoName_);
0343 
0344   std::vector<reco::TransientTrack> t_generalTrk;
0345   if (recoverOneTrackCase_)
0346     t_generalTrk = transientTrackBuilder_->build(generalTrkHandle);
0347 
0348   //  Loop over SC in the barrel and reconstruct converted photons
0349   reco::CaloClusterPtrVector scPtrVec;
0350   for (auto const& aClus : scHandle->ptrs()) {
0351     // preselection based in Et and H/E cut
0352     if (aClus->energy() / cosh(aClus->eta()) <= minSCEt_)
0353       continue;
0354     const reco::CaloCluster* pClus = &(*aClus);
0355     auto const* sc = dynamic_cast<const reco::SuperCluster*>(pClus);
0356     double HoE = hcalHelper.hcalESum(*sc, 0, hcalCuts_) / sc->energy();
0357     if (HoE >= maxHOverE_)
0358       continue;
0359     /////
0360 
0361     std::vector<edm::Ref<reco::TrackCollection>> trackPairRef;
0362     std::vector<math::XYZPointF> trackInnPos;
0363     std::vector<math::XYZVectorF> trackPin;
0364     std::vector<math::XYZVectorF> trackPout;
0365     float minAppDist = -99;
0366 
0367     //LogDebug("ConvertedPhotonProducer") << "ConvertedPhotonProducer SC energy " << aClus->energy() << " eta " <<  aClus->eta() << " phi " <<  aClus->phi() << "\n";
0368 
0369     //// Set here first quantities for the converted photon
0370     const reco::Particle::Point vtx(0, 0, 0);
0371 
0372     math::XYZVector direction = aClus->position() - vtx;
0373     math::XYZVector momentum = direction.unit() * aClus->energy();
0374     const reco::Particle::LorentzVector p4(momentum.x(), momentum.y(), momentum.z(), aClus->energy());
0375 
0376     if (!allPairs.empty()) {
0377       for (auto iPair = allPairs.begin(); iPair != allPairs.end(); ++iPair) {
0378         scPtrVec.clear();
0379 
0380         reco::Vertex theConversionVertex;
0381         reco::CaloClusterPtr caloPtr = iPair->second;
0382         if (!(aClus == caloPtr))
0383           continue;
0384 
0385         scPtrVec.push_back(aClus);
0386 
0387         std::vector<math::XYZPointF> trkPositionAtEcal = theEcalImpactPositionFinder.find(iPair->first, bcHandle);
0388         std::vector<reco::CaloClusterPtr> matchingBC = theEcalImpactPositionFinder.matchingBC();
0389 
0390         minAppDist = -99;
0391         const std::string metname = "ConvertedPhotons|ConvertedPhotonProducer";
0392         if ((iPair->first).size() > 1) {
0393           try {
0394             vertexFinder_.run(iPair->first, theConversionVertex);
0395 
0396           } catch (cms::Exception& e) {
0397             //std::cout << " cms::Exception caught in ConvertedPhotonProducer::produce" << "\n" ;
0398             edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
0399                                      << e.explainSelf();
0400           }
0401 
0402           // Old TwoTrackMinimumDistance md;
0403           // Old md.calculate  (  (iPair->first)[0].initialFreeState(),  (iPair->first)[1].initialFreeState() );
0404           // Old minAppDist = md.distance();
0405 
0406           /*
0407     for ( unsigned int i=0; i< matchingBC.size(); ++i) {
0408           if (  matchingBC[i].isNull() )  std::cout << " This ref to BC is null: skipping " <<  "\n";
0409           else 
0410         std::cout << " BC energy " << matchingBC[i]->energy() <<  "\n";
0411     }
0412     */
0413 
0414           //// loop over tracks in the pair  for creating a reference
0415           trackPairRef.clear();
0416           trackInnPos.clear();
0417           trackPin.clear();
0418           trackPout.clear();
0419 
0420           for (std::vector<reco::TransientTrack>::const_iterator iTk = (iPair->first).begin();
0421                iTk != (iPair->first).end();
0422                ++iTk) {
0423             //LogDebug("ConvertedPhotonProducer")  << "  ConvertedPhotonProducer Transient Tracks in the pair  charge " << iTk->charge() << " Num of RecHits " << iTk->recHitsSize() << " inner momentum " << iTk->track().innerMomentum() << "\n";
0424 
0425             auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
0426             reco::TrackRef myTkRef = ttt->persistentTrackRef();
0427 
0428             //LogDebug("ConvertedPhotonProducer")  << " ConvertedPhotonProducer Ref to Rec Tracks in the pair  charge " << myTkRef->charge() << " Num of RecHits " << myTkRef->recHitsSize() << " inner momentum " << myTkRef->innerMomentum() << "\n";
0429             if (myTkRef->extra().isNonnull()) {
0430               trackInnPos.push_back(toFConverterP(myTkRef->innerPosition()));
0431               trackPin.push_back(toFConverterV(myTkRef->innerMomentum()));
0432               trackPout.push_back(toFConverterV(myTkRef->outerMomentum()));
0433             }
0434             trackPairRef.push_back(myTkRef);
0435           }
0436 
0437           //    std::cout << " ConvertedPhotonProducer trackPin size " << trackPin.size() << std::endl;
0438           //LogDebug("ConvertedPhotonProducer")  << " ConvertedPhotonProducer SC energy " <<  aClus->energy() << "\n";
0439           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer photon p4 " << p4  << "\n";
0440           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer vtx " << vtx.x() << " " << vtx.y() << " " << vtx.z() << "\n";
0441           if (theConversionVertex.isValid()) {
0442             //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer theConversionVertex " << theConversionVertex.position().x() << " " << theConversionVertex.position().y() << " " << theConversionVertex.position().z() << "\n";
0443           }
0444           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer trackPairRef  " << trackPairRef.size() <<  "\n";
0445 
0446           minAppDist = calculateMinApproachDistance(trackPairRef[0], trackPairRef[1]);
0447 
0448           double like = -999.;
0449           reco::Conversion newCandidate(scPtrVec,
0450                                         trackPairRef,
0451                                         trkPositionAtEcal,
0452                                         theConversionVertex,
0453                                         matchingBC,
0454                                         minAppDist,
0455                                         trackInnPos,
0456                                         trackPin,
0457                                         trackPout,
0458                                         like,
0459                                         algo);
0460           like = likelihoodCalc_.calculateLikelihood(newCandidate);
0461           //    std::cout << "like = " << like << std::endl;
0462           newCandidate.setMVAout(like);
0463           outputConvPhotonCollection.push_back(newCandidate);
0464 
0465           //LogDebug("ConvertedPhotonProducer") << " ConvertedPhotonProducer Put the ConvertedPhotonCollection a candidate in the Barrel " << "\n";
0466 
0467         } else {
0468           //      std::cout << "   ConvertedPhotonProducer case with only one track found " <<  "\n";
0469 
0470           //std::cout << "   ConvertedPhotonProducer recovering one track " <<  "\n";
0471           trackPairRef.clear();
0472           trackInnPos.clear();
0473           trackPin.clear();
0474           trackPout.clear();
0475           std::vector<reco::TransientTrack>::const_iterator iTk = (iPair->first).begin();
0476           //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";
0477           auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(iTk->basicTransientTrack());
0478           reco::TrackRef myTk = ttt->persistentTrackRef();
0479           if (myTk->extra().isNonnull()) {
0480             trackInnPos.push_back(toFConverterP(myTk->innerPosition()));
0481             trackPin.push_back(toFConverterV(myTk->innerMomentum()));
0482             trackPout.push_back(toFConverterV(myTk->outerMomentum()));
0483           }
0484           trackPairRef.push_back(myTk);
0485           //std::cout << " Provenance " << myTk->algoName() << std::endl;
0486 
0487           if (recoverOneTrackCase_) {
0488             float theta1 = myTk->innerMomentum().Theta();
0489             float dCot = 999.;
0490             float dCotTheta = -999.;
0491             reco::TrackRef goodRef;
0492             for (auto const& tran : t_generalTrk) {
0493               auto const* ttt = dynamic_cast<const reco::TrackTransientTrack*>(tran.basicTransientTrack());
0494               reco::TrackRef trRef = ttt->persistentTrackRef();
0495               if (trRef->charge() * myTk->charge() > 0)
0496                 continue;
0497               float dEta = trRef->eta() - myTk->eta();
0498               float dPhi = trRef->phi() - myTk->phi();
0499               if (sqrt(dEta * dEta + dPhi * dPhi) > dRForConversionRecovery_)
0500                 continue;
0501               float theta2 = trRef->innerMomentum().Theta();
0502               dCotTheta = 1. / tan(theta1) - 1. / tan(theta2);
0503               //    std::cout << "  ConvertedPhotonProducer recovering general transient track charge " << trRef->charge() << " momentum " << trRef->innerMomentum() << " dcotTheta " << fabs(dCotTheta) << std::endl;
0504               if (fabs(dCotTheta) < dCot) {
0505                 dCot = fabs(dCotTheta);
0506                 goodRef = trRef;
0507               }
0508             }
0509 
0510             if (goodRef.isNonnull()) {
0511               minAppDist = calculateMinApproachDistance(myTk, goodRef);
0512 
0513               // std::cout << "  ConvertedPhotonProducer chosen dCotTheta " <<  fabs(dCotTheta) << std::endl;
0514               if (fabs(dCotTheta) < deltaCotCut_ && minAppDist > minApproachDisCut_) {
0515                 trackInnPos.push_back(toFConverterP(goodRef->innerPosition()));
0516                 trackPin.push_back(toFConverterV(goodRef->innerMomentum()));
0517                 trackPout.push_back(toFConverterV(goodRef->outerMomentum()));
0518                 trackPairRef.push_back(goodRef);
0519                 //      std::cout << " ConvertedPhotonProducer adding opposite charge track from generalTrackCollection charge " <<  goodRef ->charge() << " pt " << sqrt(goodRef->innerMomentum().perp2())  << " trackPairRef size " << trackPairRef.size() << std::endl;
0520                 //std::cout << " Track Provenenance " << goodRef->algoName() << std::endl;
0521 
0522                 try {
0523                   vertexFinder_.run(iPair->first, theConversionVertex);
0524 
0525                 } catch (cms::Exception& e) {
0526                   //std::cout << " cms::Exception caught in ConvertedPhotonProducer::produce" << "\n" ;
0527                   edm::LogWarning(metname) << "cms::Exception caught in ConvertedPhotonProducer::produce\n"
0528                                            << e.explainSelf();
0529                 }
0530               }
0531             }
0532 
0533           }  // bool On/Off one track case recovery using generalTracks
0534           const double like = -999.;
0535           outputConvPhotonCollection.emplace_back(scPtrVec,
0536                                                   trackPairRef,
0537                                                   trkPositionAtEcal,
0538                                                   theConversionVertex,
0539                                                   matchingBC,
0540                                                   minAppDist,
0541                                                   trackInnPos,
0542                                                   trackPin,
0543                                                   trackPout,
0544                                                   like,
0545                                                   algo);
0546           auto& newCandidate = outputConvPhotonCollection.back();
0547           newCandidate.setMVAout(likelihoodCalc_.calculateLikelihood(newCandidate));
0548 
0549         }  // case with only on track: looking in general tracks
0550       }
0551     }
0552   }
0553 }
0554 
0555 void ConvertedPhotonProducer::cleanCollections(const edm::Handle<edm::View<reco::CaloCluster>>& scHandle,
0556                                                const edm::OrphanHandle<reco::ConversionCollection>& conversionHandle,
0557                                                reco::ConversionCollection& outputConversionCollection) {
0558   reco::Conversion* newCandidate = nullptr;
0559   for (auto const& aClus : scHandle->ptrs()) {
0560     // SC energy preselection
0561     if (aClus->energy() / cosh(aClus->eta()) <= minSCEt_)
0562       continue;
0563 
0564     if (conversionHandle.isValid()) {
0565       if (risolveAmbiguity_) {
0566         std::vector<reco::ConversionRef> bestRef = solveAmbiguity(conversionHandle, aClus);
0567 
0568         for (std::vector<reco::ConversionRef>::iterator iRef = bestRef.begin(); iRef != bestRef.end(); iRef++) {
0569           if (iRef->isNonnull()) {
0570             newCandidate = (*iRef)->clone();
0571             outputConversionCollection.push_back(*newCandidate);
0572             delete newCandidate;
0573           }
0574         }
0575 
0576       } else {
0577         for (unsigned int icp = 0; icp < conversionHandle->size(); icp++) {
0578           reco::ConversionRef cpRef(reco::ConversionRef(conversionHandle, icp));
0579           if (!(aClus.id() == cpRef->caloCluster()[0].id() && aClus.key() == cpRef->caloCluster()[0].key()))
0580             continue;
0581           if (!cpRef->isConverted())
0582             continue;
0583           if (cpRef->nTracks() < 2)
0584             continue;
0585           newCandidate = (&(*cpRef))->clone();
0586           outputConversionCollection.push_back(*newCandidate);
0587           delete newCandidate;
0588         }
0589 
0590       }  // solve or not the ambiguity of many conversion candidates
0591     }
0592   }
0593 }
0594 
0595 std::vector<reco::ConversionRef> ConvertedPhotonProducer::solveAmbiguity(
0596     const edm::OrphanHandle<reco::ConversionCollection>& conversionHandle, reco::CaloClusterPtr const& scRef) {
0597   std::multimap<double, reco::ConversionRef, std::greater<double>> convMap;
0598 
0599   for (unsigned int icp = 0; icp < conversionHandle->size(); icp++) {
0600     reco::ConversionRef cpRef{conversionHandle, icp};
0601 
0602     //std::cout << " cpRef " << cpRef->nTracks() << " " <<  cpRef ->caloCluster()[0]->energy() << std::endl;
0603     if (!(scRef.id() == cpRef->caloCluster()[0].id() && scRef.key() == cpRef->caloCluster()[0].key()))
0604       continue;
0605     if (!cpRef->isConverted())
0606       continue;
0607     double like = cpRef->MVAout();
0608     if (cpRef->nTracks() < 2)
0609       continue;
0610     //    std::cout << " Like " << like << std::endl;
0611     convMap.emplace(like, cpRef);
0612   }
0613 
0614   //  std::cout << " convMap size " << convMap.size() << std::endl;
0615 
0616   std::vector<reco::ConversionRef> bestRefs;
0617   for (auto iMap = convMap.begin(); iMap != convMap.end(); iMap++) {
0618     //    std::cout << " Like list in the map " <<  iMap->first << " " << (iMap->second)->EoverP() << std::endl;
0619     bestRefs.push_back(iMap->second);
0620     if (int(bestRefs.size()) == maxNumOfCandidates_)
0621       break;
0622   }
0623 
0624   return bestRefs;
0625 }
0626 
0627 float ConvertedPhotonProducer::calculateMinApproachDistance(const reco::TrackRef& track1,
0628                                                             const reco::TrackRef& track2) {
0629   double x1, x2, y1, y2;
0630   double xx_1 = track1->innerPosition().x(), yy_1 = track1->innerPosition().y(), zz_1 = track1->innerPosition().z();
0631   double xx_2 = track2->innerPosition().x(), yy_2 = track2->innerPosition().y(), zz_2 = track2->innerPosition().z();
0632   double radius1 =
0633       track1->innerMomentum().Rho() / (.3 * (magneticField_->inTesla(GlobalPoint(xx_1, yy_1, zz_1)).z())) * 100;
0634   double radius2 =
0635       track2->innerMomentum().Rho() / (.3 * (magneticField_->inTesla(GlobalPoint(xx_2, yy_2, zz_2)).z())) * 100;
0636   getCircleCenter(track1, radius1, x1, y1);
0637   getCircleCenter(track2, radius2, x2, y2);
0638 
0639   return std::sqrt((x1 - x2) * (x1 - x2) + (y1 - y2) * (y1 - y2)) - radius1 - radius2;
0640 }
0641 
0642 void ConvertedPhotonProducer::getCircleCenter(const reco::TrackRef& tk, double r, double& x0, double& y0) {
0643   double x1, y1, phi;
0644   x1 = tk->innerPosition().x();  //inner position and inner momentum need track Extra!
0645   y1 = tk->innerPosition().y();
0646   phi = tk->innerMomentum().phi();
0647   const int charge = tk->charge();
0648   x0 = x1 + r * sin(phi) * charge;
0649   y0 = y1 - r * cos(phi) * charge;
0650 }
[ 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!