File indexing completed on 2023-03-17 11:17:49
0001
0002
0003
0004
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
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
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
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
0089
0090 reco::PhotonCoreCollection outputPhotonCoreCollection;
0091 auto outputPhotonCoreCollection_p = std::make_unique<reco::PhotonCoreCollection>();
0092
0093
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
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
0116 validConversions_ = true;
0117 edm::Handle<reco::ConversionCollection> conversionHandle;
0118 theEvent.getByToken(conversionProducer_, conversionHandle);
0119 if (!conversionHandle.isValid()) {
0120
0121 validConversions_ = false;
0122 }
0123
0124
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
0133
0134 int iSC = 0;
0135
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
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
0158 reco::SuperClusterRef scRef(reco::SuperClusterRef(scHandle, lSC));
0159 iSC++;
0160
0161
0162
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 }
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
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
0231
0232 float ep = 0;
0233 if (max_lh < 0) {
0234
0235
0236 float epMin = 999;
0237
0238 for (iMap = convMap.begin(); iMap != convMap.end(); iMap++) {
0239 reco::ConversionRef convRef = iMap->first;
0240
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
0248 if (ep < epMin) {
0249 epMin = ep;
0250 bestRef = iMap->first;
0251 }
0252 }
0253
0254 }
0255
0256 return bestRef;
0257 }