Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2023-03-17 11:17:49 

0001 /** \class PhotonCoreProducer
0002  **  
0003  **
0004  **  \author Nancy Marinelli, U. of Notre Dame, US
0005  **
0006  ***/
0007 
0008 #include "DataFormats/Common/interface/Handle.h"
0009 #include "DataFormats/EgammaCandidates/interface/Conversion.h"
0010 #include "DataFormats/EgammaCandidates/interface/Photon.h"
0011 #include "DataFormats/EgammaCandidates/interface/PhotonCore.h"
0012 #include "DataFormats/EgammaReco/interface/BasicCluster.h"
0013 #include "DataFormats/EgammaReco/interface/BasicClusterShapeAssociation.h"
0014 #include "DataFormats/EgammaReco/interface/ClusterShape.h"
0015 #include "DataFormats/EgammaReco/interface/ElectronSeed.h"
0016 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0017 #include "FWCore/Framework/interface/ESHandle.h"
0018 #include "FWCore/Framework/interface/Event.h"
0019 #include "FWCore/Framework/interface/EventSetup.h"
0020 #include "FWCore/Framework/interface/stream/EDProducer.h"
0021 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0022 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0023 #include "FWCore/Utilities/interface/Exception.h"
0024 
0025 #include <vector>
0026 
0027 // PhotonCoreProducer inherits from EDProducer, so it can be a module:
0028 class PhotonCoreProducer : public edm::stream::EDProducer<> {
0029 public:
0030   PhotonCoreProducer(const edm::ParameterSet& ps);
0031   ~PhotonCoreProducer() override;
0032 
0033   void produce(edm::Event& evt, const edm::EventSetup& es) override;
0034 
0035 private:
0036   void fillPhotonCollection(edm::Event& evt,
0037                             edm::EventSetup const& es,
0038                             const edm::Handle<reco::SuperClusterCollection>& scHandle,
0039                             const edm::Handle<reco::ConversionCollection>& conversionHandle,
0040                             const edm::Handle<reco::ElectronSeedCollection>& pixelSeeds,
0041                             reco::PhotonCoreCollection& outputCollection,
0042                             int& iSC);
0043 
0044   reco::ConversionRef solveAmbiguity(const edm::Handle<reco::ConversionCollection>& conversionHandle,
0045                                      reco::SuperClusterRef& sc);
0046 
0047   std::string PhotonCoreCollection_;
0048   edm::EDGetTokenT<reco::SuperClusterCollection> scHybridBarrelProducer_;
0049   edm::EDGetTokenT<reco::SuperClusterCollection> scIslandEndcapProducer_;
0050   edm::EDGetTokenT<reco::ConversionCollection> conversionProducer_;
0051   edm::EDGetTokenT<reco::ElectronSeedCollection> pixelSeedProducer_;
0052 
0053   double minSCEt_;
0054   bool validConversions_;
0055   edm::ParameterSet conf_;
0056   bool validPixelSeeds_;
0057   bool risolveAmbiguity_;
0058   bool endcapOnly_;
0059 };
0060 
0061 #include "FWCore/Framework/interface/MakerMacros.h"
0062 DEFINE_FWK_MODULE(PhotonCoreProducer);
0063 
0064 PhotonCoreProducer::PhotonCoreProducer(const edm::ParameterSet& config)
0065     : conf_(config)
0066 
0067 {
0068   // use onfiguration file to setup input/output collection names
0069   scHybridBarrelProducer_ =
0070       consumes<reco::SuperClusterCollection>(conf_.getParameter<edm::InputTag>("scHybridBarrelProducer"));
0071   scIslandEndcapProducer_ =
0072       consumes<reco::SuperClusterCollection>(conf_.getParameter<edm::InputTag>("scIslandEndcapProducer"));
0073   conversionProducer_ = consumes<reco::ConversionCollection>(conf_.getParameter<edm::InputTag>("conversionProducer"));
0074   PhotonCoreCollection_ = conf_.getParameter<std::string>("photonCoreCollection");
0075   pixelSeedProducer_ = consumes<reco::ElectronSeedCollection>(conf_.getParameter<edm::InputTag>("pixelSeedProducer"));
0076   minSCEt_ = conf_.getParameter<double>("minSCEt");
0077   risolveAmbiguity_ = conf_.getParameter<bool>("risolveConversionAmbiguity");
0078   endcapOnly_ = conf_.getParameter<bool>("endcapOnly");
0079 
0080   // Register the product
0081   produces<reco::PhotonCoreCollection>(PhotonCoreCollection_);
0082 }
0083 
0084 PhotonCoreProducer::~PhotonCoreProducer() {}
0085 
0086 void PhotonCoreProducer::produce(edm::Event& theEvent, const edm::EventSetup& theEventSetup) {
0087   using namespace edm;
0088   //  nEvt_++;
0089 
0090   reco::PhotonCoreCollection outputPhotonCoreCollection;
0091   auto outputPhotonCoreCollection_p = std::make_unique<reco::PhotonCoreCollection>();
0092 
0093   // Get the  Barrel Super Cluster collection
0094   bool validBarrelSCHandle = true;
0095   if (endcapOnly_) {
0096     validBarrelSCHandle = false;
0097   }
0098 
0099   Handle<reco::SuperClusterCollection> scBarrelHandle;
0100   theEvent.getByToken(scHybridBarrelProducer_, scBarrelHandle);
0101   if (!endcapOnly_ && !scBarrelHandle.isValid()) {
0102     edm::LogError("PhotonCoreProducer") << "Error! Can't get the scHybridBarrelProducer";
0103     validBarrelSCHandle = false;
0104   }
0105 
0106   // Get the  Endcap Super Cluster collection
0107   bool validEndcapSCHandle = true;
0108   Handle<reco::SuperClusterCollection> scEndcapHandle;
0109   theEvent.getByToken(scIslandEndcapProducer_, scEndcapHandle);
0110   if (!scEndcapHandle.isValid()) {
0111     edm::LogError("PhotonCoreProducer") << "Error! Can't get the scIslandEndcapProducer";
0112     validEndcapSCHandle = false;
0113   }
0114 
0115   ///// Get the conversion collection
0116   validConversions_ = true;
0117   edm::Handle<reco::ConversionCollection> conversionHandle;
0118   theEvent.getByToken(conversionProducer_, conversionHandle);
0119   if (!conversionHandle.isValid()) {
0120     //edm::LogError("PhotonCoreProducer") << "Error! Can't get the product "<< conversionProducer_.label() << "\n" ;
0121     validConversions_ = false;
0122   }
0123 
0124   // Get ElectronPixelSeeds
0125   validPixelSeeds_ = true;
0126   Handle<reco::ElectronSeedCollection> pixelSeedHandle;
0127   reco::ElectronSeedCollection pixelSeeds;
0128   theEvent.getByToken(pixelSeedProducer_, pixelSeedHandle);
0129   if (!pixelSeedHandle.isValid()) {
0130     validPixelSeeds_ = false;
0131   }
0132   //  if ( validPixelSeeds_) pixelSeeds = *(pixelSeedHandle.product());
0133 
0134   int iSC = 0;  // index in photon collection
0135   // Loop over barrel and endcap SC collections and fill the  photon collection
0136   if (validBarrelSCHandle)
0137     fillPhotonCollection(
0138         theEvent, theEventSetup, scBarrelHandle, conversionHandle, pixelSeedHandle, outputPhotonCoreCollection, iSC);
0139   if (validEndcapSCHandle)
0140     fillPhotonCollection(
0141         theEvent, theEventSetup, scEndcapHandle, conversionHandle, pixelSeedHandle, outputPhotonCoreCollection, iSC);
0142 
0143   // put the product in the event
0144   edm::LogInfo("PhotonCoreProducer") << " Put in the event " << iSC << " Photon Candidates \n";
0145   outputPhotonCoreCollection_p->assign(outputPhotonCoreCollection.begin(), outputPhotonCoreCollection.end());
0146   theEvent.put(std::move(outputPhotonCoreCollection_p), PhotonCoreCollection_);
0147 }
0148 
0149 void PhotonCoreProducer::fillPhotonCollection(edm::Event& evt,
0150                                               edm::EventSetup const& es,
0151                                               const edm::Handle<reco::SuperClusterCollection>& scHandle,
0152                                               const edm::Handle<reco::ConversionCollection>& conversionHandle,
0153                                               const edm::Handle<reco::ElectronSeedCollection>& pixelSeedHandle,
0154                                               reco::PhotonCoreCollection& outputPhotonCoreCollection,
0155                                               int& iSC) {
0156   for (unsigned int lSC = 0; lSC < scHandle->size(); lSC++) {
0157     // get SuperClusterRef
0158     reco::SuperClusterRef scRef(reco::SuperClusterRef(scHandle, lSC));
0159     iSC++;
0160     //const reco::SuperCluster* pClus=&(*scRef);
0161 
0162     // SC energy preselection
0163     if (scRef->energy() / cosh(scRef->eta()) <= minSCEt_)
0164       continue;
0165 
0166     reco::PhotonCore newCandidate(scRef);
0167     newCandidate.setParentSuperCluster(scRef);
0168     if (validConversions_) {
0169       if (risolveAmbiguity_) {
0170         reco::ConversionRef bestRef = solveAmbiguity(conversionHandle, scRef);
0171         if (bestRef.isNonnull())
0172           newCandidate.addConversion(bestRef);
0173 
0174       } else {
0175         for (unsigned int icp = 0; icp < conversionHandle->size(); icp++) {
0176           reco::ConversionRef cpRef(reco::ConversionRef(conversionHandle, icp));
0177           if (cpRef->caloCluster().empty())
0178             continue;
0179           if (!(scRef.id() == cpRef->caloCluster()[0].id() && scRef.key() == cpRef->caloCluster()[0].key()))
0180             continue;
0181           if (!cpRef->isConverted())
0182             continue;
0183           newCandidate.addConversion(cpRef);
0184         }
0185 
0186       }  // solve or not the ambiguity of many conversion candidates
0187     }
0188 
0189     if (validPixelSeeds_) {
0190       for (unsigned int icp = 0; icp < pixelSeedHandle->size(); icp++) {
0191         reco::ElectronSeedRef cpRef(reco::ElectronSeedRef(pixelSeedHandle, icp));
0192         if (!cpRef->isEcalDriven())
0193           continue;
0194         if (!(scRef.id() == cpRef->caloCluster().id() && scRef.key() == cpRef->caloCluster().key()))
0195           continue;
0196         newCandidate.addElectronPixelSeed(cpRef);
0197       }
0198     }
0199 
0200     outputPhotonCoreCollection.push_back(newCandidate);
0201   }
0202 }
0203 
0204 reco::ConversionRef PhotonCoreProducer::solveAmbiguity(const edm::Handle<reco::ConversionCollection>& conversionHandle,
0205                                                        reco::SuperClusterRef& scRef) {
0206   std::multimap<reco::ConversionRef, double> convMap;
0207   for (unsigned int icp = 0; icp < conversionHandle->size(); icp++) {
0208     reco::ConversionRef cpRef(reco::ConversionRef(conversionHandle, icp));
0209 
0210     if (!(scRef.id() == cpRef->caloCluster()[0].id() && scRef.key() == cpRef->caloCluster()[0].key()))
0211       continue;
0212     if (!cpRef->isConverted())
0213       continue;
0214     double like = cpRef->MVAout();
0215     convMap.insert(std::make_pair(cpRef, like));
0216   }
0217 
0218   std::multimap<reco::ConversionRef, double>::iterator iMap;
0219   double max_lh = -1.;
0220   reco::ConversionRef bestRef;
0221   //  std::cout << " Pick up the best conv " << std::endl;
0222   for (iMap = convMap.begin(); iMap != convMap.end(); iMap++) {
0223     double like = iMap->second;
0224     if (like > max_lh) {
0225       max_lh = like;
0226       bestRef = iMap->first;
0227     }
0228   }
0229 
0230   //std::cout << " Best conv like " << max_lh << std::endl;
0231 
0232   float ep = 0;
0233   if (max_lh < 0) {
0234     //    std::cout << " Candidates with only one track " << std::endl;
0235     /// only one track reconstructed. Pick the one with best E/P
0236     float epMin = 999;
0237 
0238     for (iMap = convMap.begin(); iMap != convMap.end(); iMap++) {
0239       reco::ConversionRef convRef = iMap->first;
0240       // std::vector<reco::TrackRef> tracks = convRef->tracks();
0241       const std::vector<edm::RefToBase<reco::Track> > tracks = convRef->tracks();
0242       float px = tracks[0]->innerMomentum().x();
0243       float py = tracks[0]->innerMomentum().y();
0244       float pz = tracks[0]->innerMomentum().z();
0245       float p = sqrt(px * px + py * py + pz * pz);
0246       ep = fabs(1. - convRef->caloCluster()[0]->energy() / p);
0247       //    std::cout << " 1-E/P = " << ep << std::endl;
0248       if (ep < epMin) {
0249         epMin = ep;
0250         bestRef = iMap->first;
0251       }
0252     }
0253     //  std::cout << " Best conv 1-E/P " << ep << std::endl;
0254   }
0255 
0256   return bestRef;
0257 }
[ 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!