File indexing completed on 2021-05-18 09:15:25
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010 #include "FWCore/Framework/interface/ESHandle.h"
0011 #include "FWCore/Framework/interface/EventSetup.h"
0012 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0013 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibration.h"
0014 #include "RecoParticleFlow/PFClusterTools/interface/PFEnergyCalibrationHF.h"
0015 #include "CondFormats/PhysicsToolsObjects/interface/PerformancePayloadFromTFormula.h"
0016 #include "CondFormats/DataRecord/interface/PFCalibrationRcd.h"
0017 #include "CondFormats/DataRecord/interface/GBRWrapperRcd.h"
0018 #include "DataFormats/EgammaCandidates/interface/GsfElectron.h"
0019 #include "DataFormats/ParticleFlowReco/interface/PFRecHit.h"
0020 #include "DataFormats/ParticleFlowReco/interface/PFBlockElementSuperCluster.h"
0021 #include "CondFormats/DataRecord/interface/ESEEIntercalibConstantsRcd.h"
0022 #include "CondFormats/DataRecord/interface/ESChannelStatusRcd.h"
0023 #include "CondFormats/ESObjects/interface/ESEEIntercalibConstants.h"
0024 #include "CondFormats/ESObjects/interface/ESChannelStatus.h"
0025 #include "DataFormats/Common/interface/RefToPtr.h"
0026 #include "FWCore/Framework/interface/global/EDProducer.h"
0027 #include "FWCore/Framework/interface/Event.h"
0028 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0029 #include "DataFormats/VertexReco/interface/Vertex.h"
0030 #include "DataFormats/EgammaReco/interface/PreshowerCluster.h"
0031 #include "DataFormats/EgammaReco/interface/SuperCluster.h"
0032 #include "DataFormats/CaloRecHit/interface/CaloCluster.h"
0033 #include "DataFormats/ParticleFlowReco/interface/PFBlock.h"
0034 #include "RecoParticleFlow/PFProducer/interface/PFEGammaAlgo.h"
0035
0036 class PFEGammaProducer : public edm::global::EDProducer<> {
0037 public:
0038 explicit PFEGammaProducer(const edm::ParameterSet&);
0039
0040 void produce(edm::StreamID, edm::Event&, const edm::EventSetup&) const override;
0041
0042 static void fillDescriptions(edm::ConfigurationDescriptions& descriptions);
0043
0044 private:
0045 reco::ConversionCollection createSingleLegConversions(reco::PFCandidateEGammaExtraCollection& extras,
0046 const edm::RefProd<reco::ConversionCollection>& convProd) const;
0047
0048 const edm::EDGetTokenT<reco::PFBlockCollection> inputTagBlocks_;
0049 const edm::EDGetTokenT<reco::PFCluster::EEtoPSAssociation> eetopsSrc_;
0050 const edm::EDGetTokenT<reco::VertexCollection> vertices_;
0051
0052 const edm::EDPutTokenT<reco::PFCandidateCollection> pfCandidateCollectionPutToken_;
0053 const edm::EDPutTokenT<reco::PFCandidateEGammaExtraCollection> pfCandidateEGammaExtraCollectionPutToken_;
0054 const edm::EDPutTokenT<reco::SuperClusterCollection> superClusterCollectionPutToken_;
0055 const edm::EDPutTokenT<reco::CaloClusterCollection> caloClusterCollectionEBEEPutToken_;
0056 const edm::EDPutTokenT<reco::CaloClusterCollection> caloClusterCollectionESPutToken_;
0057 const edm::EDPutTokenT<reco::ConversionCollection> conversionCollectionPutToken_;
0058
0059
0060 const PFEGammaAlgo::PFEGConfigInfo pfEGConfigInfo_;
0061 const PFEGammaAlgo::GBRForests gbrForests_;
0062
0063 const edm::ESGetToken<ESEEIntercalibConstants, ESEEIntercalibConstantsRcd> esEEInterCalibToken_;
0064 const edm::ESGetToken<ESChannelStatus, ESChannelStatusRcd> esChannelStatusToken_;
0065 };
0066
0067 #include "FWCore/Framework/interface/MakerMacros.h"
0068 DEFINE_FWK_MODULE(PFEGammaProducer);
0069
0070 #ifdef PFLOW_DEBUG
0071 #define LOGDRESSED(x) edm::LogInfo(x)
0072 #else
0073 #define LOGDRESSED(x) LogDebug(x)
0074 #endif
0075
0076 PFEGammaProducer::PFEGammaProducer(const edm::ParameterSet& iConfig)
0077 : inputTagBlocks_(consumes<reco::PFBlockCollection>(iConfig.getParameter<edm::InputTag>("blocks"))),
0078 eetopsSrc_(consumes<reco::PFCluster::EEtoPSAssociation>(iConfig.getParameter<edm::InputTag>("EEtoPS_source"))),
0079 vertices_(consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexCollection"))),
0080 pfCandidateCollectionPutToken_{produces<reco::PFCandidateCollection>()},
0081 pfCandidateEGammaExtraCollectionPutToken_{produces<reco::PFCandidateEGammaExtraCollection>()},
0082 superClusterCollectionPutToken_{produces<reco::SuperClusterCollection>()},
0083 caloClusterCollectionEBEEPutToken_{produces<reco::CaloClusterCollection>("EBEEClusters")},
0084 caloClusterCollectionESPutToken_{produces<reco::CaloClusterCollection>("ESClusters")},
0085 conversionCollectionPutToken_{produces<reco::ConversionCollection>()},
0086 pfEGConfigInfo_{
0087 .mvaEleCut = iConfig.getParameter<double>("pf_electron_mvaCut"),
0088 .applyCrackCorrections = iConfig.getParameter<bool>("pf_electronID_crackCorrection"),
0089 .produceEGCandsWithNoSuperCluster = iConfig.getParameter<bool>("produceEGCandsWithNoSuperCluster"),
0090 .mvaConvCut = iConfig.getParameter<double>("pf_conv_mvaCut"),
0091 },
0092 gbrForests_{
0093 iConfig,
0094 },
0095 esEEInterCalibToken_(esConsumes()),
0096 esChannelStatusToken_(esConsumes()) {}
0097
0098 void PFEGammaProducer::produce(edm::StreamID, edm::Event& iEvent, const edm::EventSetup& iSetup) const {
0099 LOGDRESSED("PFEGammaProducer") << "START event: " << iEvent.id().event() << " in run " << iEvent.id().run()
0100 << std::endl;
0101
0102
0103 reco::PFCandidateCollection egCandidates{};
0104 reco::PFCandidateEGammaExtraCollection egExtra{};
0105 reco::SuperClusterCollection sClusters{};
0106
0107
0108 auto esEEInterCalibHandle_ = iSetup.getHandle(esEEInterCalibToken_);
0109 auto esChannelStatusHandle_ = iSetup.getHandle(esChannelStatusToken_);
0110
0111
0112 auto const& primaryVertices = iEvent.get(vertices_);
0113 auto const* primaryVertex = &primaryVertices.front();
0114 for (auto const& pv : primaryVertices) {
0115 if (pv.isValid() && !pv.isFake()) {
0116 primaryVertex = &pv;
0117 break;
0118 }
0119 }
0120
0121 PFEGammaAlgo pfEGammaAlgo{pfEGConfigInfo_,
0122 gbrForests_,
0123 iEvent.get(eetopsSrc_),
0124 *esEEInterCalibHandle_,
0125 *esChannelStatusHandle_,
0126 *primaryVertex};
0127
0128
0129
0130 LOGDRESSED("PFEGammaProducer") << "getting blocks" << std::endl;
0131 auto blocks = iEvent.getHandle(inputTagBlocks_);
0132
0133 LOGDRESSED("PFEGammaProducer") << "EGPFlow is starting..." << std::endl;
0134
0135 #ifdef PFLOW_DEBUG
0136 assert(blocks.isValid() && "edm::Handle to blocks was null!");
0137 std::ostringstream str;
0138
0139
0140 LOGDRESSED("PFEGammaProducer") << str.str() << std::endl;
0141 #endif
0142
0143
0144 std::list<reco::PFBlockRef> hcalBlockRefs;
0145 std::list<reco::PFBlockRef> ecalBlockRefs;
0146 std::list<reco::PFBlockRef> hoBlockRefs;
0147 std::list<reco::PFBlockRef> otherBlockRefs;
0148
0149 for (unsigned i = 0; i < blocks->size(); ++i) {
0150 reco::PFBlockRef blockref(blocks, i);
0151
0152 const edm::OwnVector<reco::PFBlockElement>& elements = blockref->elements();
0153
0154 LOGDRESSED("PFEGammaProducer") << "Found " << elements.size() << " PFBlockElements in block: " << i << std::endl;
0155
0156 bool singleEcalOrHcal = false;
0157 if (elements.size() == 1) {
0158 switch (elements[0].type()) {
0159 case reco::PFBlockElement::SC:
0160 edm::LogError("PFEGammaProducer") << "PFBLOCKALGO BUG!!!! Found a SuperCluster in a block by itself!";
0161 break;
0162 case reco::PFBlockElement::PS1:
0163 case reco::PFBlockElement::PS2:
0164 case reco::PFBlockElement::ECAL:
0165 ecalBlockRefs.push_back(blockref);
0166 singleEcalOrHcal = true;
0167 break;
0168 case reco::PFBlockElement::HFEM:
0169 case reco::PFBlockElement::HFHAD:
0170 case reco::PFBlockElement::HCAL:
0171 if (elements[0].clusterRef()->flags() & reco::CaloCluster::badHcalMarker)
0172 continue;
0173 hcalBlockRefs.push_back(blockref);
0174 singleEcalOrHcal = true;
0175 break;
0176 case reco::PFBlockElement::HO:
0177
0178 hoBlockRefs.push_back(blockref);
0179 singleEcalOrHcal = true;
0180 break;
0181 default:
0182 break;
0183 }
0184 }
0185
0186 if (!singleEcalOrHcal) {
0187 otherBlockRefs.push_back(blockref);
0188 }
0189 }
0190
0191
0192
0193
0194
0195 unsigned nblcks = 0;
0196
0197
0198 for (const auto& blockref : otherBlockRefs) {
0199 ++nblcks;
0200
0201 const auto& elements = blockref->elements();
0202
0203
0204
0205 auto output = pfEGammaAlgo(blockref);
0206
0207 if (!output.candidates.empty()) {
0208 LOGDRESSED("PFEGammaProducer") << "Block with " << elements.size() << " elements produced "
0209 << output.candidates.size() << " e-g candidates!" << std::endl;
0210 }
0211
0212 const size_t egsize = egCandidates.size();
0213 egCandidates.resize(egsize + output.candidates.size());
0214 std::move(output.candidates.begin(), output.candidates.end(), egCandidates.begin() + egsize);
0215
0216 const size_t egxsize = egExtra.size();
0217 egExtra.resize(egxsize + output.candidateExtras.size());
0218 std::move(output.candidateExtras.begin(), output.candidateExtras.end(), egExtra.begin() + egxsize);
0219
0220 const size_t rscsize = sClusters.size();
0221 sClusters.resize(rscsize + output.refinedSuperClusters.size());
0222 std::move(output.refinedSuperClusters.begin(), output.refinedSuperClusters.end(), sClusters.begin() + rscsize);
0223 }
0224
0225 LOGDRESSED("PFEGammaProducer") << "Running PFEGammaAlgo on all blocks produced = " << egCandidates.size()
0226 << " e-g candidates!" << std::endl;
0227
0228 auto sClusterProd = iEvent.getRefBeforePut<reco::SuperClusterCollection>();
0229 auto egXtraProd = iEvent.getRefBeforePut<reco::PFCandidateEGammaExtraCollection>();
0230
0231
0232 for (unsigned int i = 0; i < egCandidates.size(); ++i) {
0233 reco::PFCandidate& cand = egCandidates.at(i);
0234 reco::PFCandidateEGammaExtra& xtra = egExtra.at(i);
0235
0236 reco::PFCandidateEGammaExtraRef extraref(egXtraProd, i);
0237 reco::SuperClusterRef refinedSCRef(sClusterProd, i);
0238
0239 xtra.setSuperClusterRef(refinedSCRef);
0240 cand.setSuperClusterRef(refinedSCRef);
0241 cand.setPFEGammaExtraRef(extraref);
0242 }
0243
0244
0245 reco::CaloClusterCollection caloClustersEBEE{};
0246 reco::CaloClusterCollection caloClustersES{};
0247
0248 std::map<edm::Ptr<reco::CaloCluster>, unsigned int> pfClusterMapEBEE;
0249 std::map<edm::Ptr<reco::CaloCluster>, unsigned int> pfClusterMapES;
0250
0251 for (const auto& sc : sClusters) {
0252 for (reco::CaloCluster_iterator pfclus = sc.clustersBegin(); pfclus != sc.clustersEnd(); ++pfclus) {
0253 if (!pfClusterMapEBEE.count(*pfclus)) {
0254 reco::CaloCluster caloclus(**pfclus);
0255 caloClustersEBEE.push_back(caloclus);
0256 pfClusterMapEBEE[*pfclus] = caloClustersEBEE.size() - 1;
0257 } else {
0258 throw cms::Exception("PFEgammaProducer::produce")
0259 << "Found an EB/EE pfcluster matched to more than one supercluster!" << std::dec << std::endl;
0260 }
0261 }
0262 for (reco::CaloCluster_iterator pfclus = sc.preshowerClustersBegin(); pfclus != sc.preshowerClustersEnd();
0263 ++pfclus) {
0264 if (!pfClusterMapES.count(*pfclus)) {
0265 reco::CaloCluster caloclus(**pfclus);
0266 caloClustersES.push_back(caloclus);
0267 pfClusterMapES[*pfclus] = caloClustersES.size() - 1;
0268 } else {
0269 throw cms::Exception("PFEgammaProducer::produce")
0270 << "Found an ES pfcluster matched to more than one supercluster!" << std::dec << std::endl;
0271 }
0272 }
0273 }
0274
0275
0276 auto const& caloClusHandleEBEE = iEvent.emplace(caloClusterCollectionEBEEPutToken_, std::move(caloClustersEBEE));
0277 auto const& caloClusHandleES = iEvent.emplace(caloClusterCollectionESPutToken_, std::move(caloClustersES));
0278
0279
0280 for (auto& sc : sClusters) {
0281 edm::Ptr<reco::CaloCluster> seedptr(caloClusHandleEBEE, pfClusterMapEBEE[sc.seed()]);
0282 sc.setSeed(seedptr);
0283
0284 reco::CaloClusterPtrVector clusters;
0285 for (reco::CaloCluster_iterator pfclus = sc.clustersBegin(); pfclus != sc.clustersEnd(); ++pfclus) {
0286 edm::Ptr<reco::CaloCluster> clusptr(caloClusHandleEBEE, pfClusterMapEBEE[*pfclus]);
0287 clusters.push_back(clusptr);
0288 }
0289 sc.setClusters(clusters);
0290
0291 reco::CaloClusterPtrVector psclusters;
0292 for (reco::CaloCluster_iterator pfclus = sc.preshowerClustersBegin(); pfclus != sc.preshowerClustersEnd();
0293 ++pfclus) {
0294 edm::Ptr<reco::CaloCluster> clusptr(caloClusHandleES, pfClusterMapES[*pfclus]);
0295 psclusters.push_back(clusptr);
0296 }
0297 sc.setPreshowerClusters(psclusters);
0298 }
0299
0300
0301 auto singleLegConv = createSingleLegConversions(egExtra, iEvent.getRefBeforePut<reco::ConversionCollection>());
0302
0303
0304 iEvent.emplace(superClusterCollectionPutToken_, std::move(sClusters));
0305 iEvent.emplace(pfCandidateEGammaExtraCollectionPutToken_, std::move(egExtra));
0306 iEvent.emplace(conversionCollectionPutToken_, std::move(singleLegConv));
0307 iEvent.emplace(pfCandidateCollectionPutToken_, std::move(egCandidates));
0308 }
0309
0310 reco::ConversionCollection PFEGammaProducer::createSingleLegConversions(
0311 reco::PFCandidateEGammaExtraCollection& extras, const edm::RefProd<reco::ConversionCollection>& convProd) const {
0312 reco::ConversionCollection oneLegConversions{};
0313 math::Error<3>::type error;
0314 for (auto& extra : extras) {
0315 for (const auto& tkrefmva : extra.singleLegConvTrackRefMva()) {
0316 const reco::Track& trk = *tkrefmva.first;
0317
0318 const reco::Vertex convVtx(trk.innerPosition(), error);
0319 std::vector<reco::TrackRef> OneLegConvVector;
0320 OneLegConvVector.push_back(tkrefmva.first);
0321 std::vector<float> OneLegMvaVector;
0322 OneLegMvaVector.push_back(tkrefmva.second);
0323 std::vector<reco::CaloClusterPtr> dummymatchingBC;
0324 reco::CaloClusterPtrVector scPtrVec;
0325 scPtrVec.push_back(edm::refToPtr(extra.superClusterRef()));
0326
0327 std::vector<math::XYZPointF> trackPositionAtEcalVec;
0328 std::vector<math::XYZPointF> innPointVec;
0329 std::vector<math::XYZVectorF> trackPinVec;
0330 std::vector<math::XYZVectorF> trackPoutVec;
0331 math::XYZPointF trackPositionAtEcal(trk.outerPosition().X(), trk.outerPosition().Y(), trk.outerPosition().Z());
0332 trackPositionAtEcalVec.push_back(trackPositionAtEcal);
0333
0334 math::XYZPointF innPoint(trk.innerPosition().X(), trk.innerPosition().Y(), trk.innerPosition().Z());
0335 innPointVec.push_back(innPoint);
0336
0337 math::XYZVectorF trackPin(trk.innerMomentum().X(), trk.innerMomentum().Y(), trk.innerMomentum().Z());
0338 trackPinVec.push_back(trackPin);
0339
0340 math::XYZVectorF trackPout(trk.outerMomentum().X(), trk.outerMomentum().Y(), trk.outerMomentum().Z());
0341 trackPoutVec.push_back(trackPout);
0342
0343 float DCA = trk.d0();
0344 float mvaval = tkrefmva.second;
0345 reco::Conversion singleLegConvCandidate(scPtrVec,
0346 OneLegConvVector,
0347 trackPositionAtEcalVec,
0348 convVtx,
0349 dummymatchingBC,
0350 DCA,
0351 innPointVec,
0352 trackPinVec,
0353 trackPoutVec,
0354 mvaval,
0355 reco::Conversion::pflow);
0356 singleLegConvCandidate.setOneLegMVA(OneLegMvaVector);
0357 oneLegConversions.push_back(singleLegConvCandidate);
0358
0359 reco::ConversionRef convref(convProd, oneLegConversions.size() - 1);
0360 extra.addSingleLegConversionRef(convref);
0361 }
0362 }
0363 return oneLegConversions;
0364 }
0365
0366 void PFEGammaProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0367 edm::ParameterSetDescription desc;
0368 desc.add<bool>("produceEGCandsWithNoSuperCluster", false)
0369 ->setComment("Allow building of candidates with no input or output supercluster?");
0370 desc.add<double>("pf_electron_mvaCut", -0.1);
0371 desc.add<bool>("pf_electronID_crackCorrection", false);
0372 desc.add<double>("pf_conv_mvaCut", 0.0);
0373 desc.add<edm::InputTag>("blocks", edm::InputTag("particleFlowBlock"))->setComment("PF Blocks label");
0374 desc.add<edm::InputTag>("EEtoPS_source", edm::InputTag("particleFlowClusterECAL"))
0375 ->setComment("EE to PS association");
0376 desc.add<edm::InputTag>("vertexCollection", edm::InputTag("offlinePrimaryVertices"));
0377 desc.add<edm::FileInPath>("pf_electronID_mvaWeightFile",
0378 edm::FileInPath("RecoParticleFlow/PFProducer/data/PfElectrons23Jan_BDT.weights.xml.gz"));
0379 desc.add<edm::FileInPath>("pf_convID_mvaWeightFile",
0380 edm::FileInPath("RecoParticleFlow/PFProducer/data/pfConversionAug0411_BDT.weights.xml.gz"));
0381 descriptions.add("particleFlowEGamma", desc);
0382 }