File indexing completed on 2024-04-06 12:27:29
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
0196
0197 for (const auto& blockref : otherBlockRefs) {
0198
0199 const auto& elements = blockref->elements();
0200
0201
0202
0203 auto output = pfEGammaAlgo(blockref);
0204
0205 if (!output.candidates.empty()) {
0206 LOGDRESSED("PFEGammaProducer") << "Block with " << elements.size() << " elements produced "
0207 << output.candidates.size() << " e-g candidates!" << std::endl;
0208 }
0209
0210 const size_t egsize = egCandidates.size();
0211 egCandidates.resize(egsize + output.candidates.size());
0212 std::move(output.candidates.begin(), output.candidates.end(), egCandidates.begin() + egsize);
0213
0214 const size_t egxsize = egExtra.size();
0215 egExtra.resize(egxsize + output.candidateExtras.size());
0216 std::move(output.candidateExtras.begin(), output.candidateExtras.end(), egExtra.begin() + egxsize);
0217
0218 const size_t rscsize = sClusters.size();
0219 sClusters.resize(rscsize + output.refinedSuperClusters.size());
0220 std::move(output.refinedSuperClusters.begin(), output.refinedSuperClusters.end(), sClusters.begin() + rscsize);
0221 }
0222
0223 LOGDRESSED("PFEGammaProducer") << "Running PFEGammaAlgo on all blocks produced = " << egCandidates.size()
0224 << " e-g candidates!" << std::endl;
0225
0226 auto sClusterProd = iEvent.getRefBeforePut<reco::SuperClusterCollection>();
0227 auto egXtraProd = iEvent.getRefBeforePut<reco::PFCandidateEGammaExtraCollection>();
0228
0229
0230 for (unsigned int i = 0; i < egCandidates.size(); ++i) {
0231 reco::PFCandidate& cand = egCandidates.at(i);
0232 reco::PFCandidateEGammaExtra& xtra = egExtra.at(i);
0233
0234 reco::PFCandidateEGammaExtraRef extraref(egXtraProd, i);
0235 reco::SuperClusterRef refinedSCRef(sClusterProd, i);
0236
0237 xtra.setSuperClusterRef(refinedSCRef);
0238 cand.setSuperClusterRef(refinedSCRef);
0239 cand.setPFEGammaExtraRef(extraref);
0240 }
0241
0242
0243 reco::CaloClusterCollection caloClustersEBEE{};
0244 reco::CaloClusterCollection caloClustersES{};
0245
0246 std::map<edm::Ptr<reco::CaloCluster>, unsigned int> pfClusterMapEBEE;
0247 std::map<edm::Ptr<reco::CaloCluster>, unsigned int> pfClusterMapES;
0248
0249 for (const auto& sc : sClusters) {
0250 for (reco::CaloCluster_iterator pfclus = sc.clustersBegin(); pfclus != sc.clustersEnd(); ++pfclus) {
0251 if (!pfClusterMapEBEE.count(*pfclus)) {
0252 reco::CaloCluster caloclus(**pfclus);
0253 caloClustersEBEE.push_back(caloclus);
0254 pfClusterMapEBEE[*pfclus] = caloClustersEBEE.size() - 1;
0255 } else {
0256 throw cms::Exception("PFEgammaProducer::produce")
0257 << "Found an EB/EE pfcluster matched to more than one supercluster!" << std::dec << std::endl;
0258 }
0259 }
0260 for (reco::CaloCluster_iterator pfclus = sc.preshowerClustersBegin(); pfclus != sc.preshowerClustersEnd();
0261 ++pfclus) {
0262 if (!pfClusterMapES.count(*pfclus)) {
0263 reco::CaloCluster caloclus(**pfclus);
0264 caloClustersES.push_back(caloclus);
0265 pfClusterMapES[*pfclus] = caloClustersES.size() - 1;
0266 } else {
0267 throw cms::Exception("PFEgammaProducer::produce")
0268 << "Found an ES pfcluster matched to more than one supercluster!" << std::dec << std::endl;
0269 }
0270 }
0271 }
0272
0273
0274 auto const& caloClusHandleEBEE = iEvent.emplace(caloClusterCollectionEBEEPutToken_, std::move(caloClustersEBEE));
0275 auto const& caloClusHandleES = iEvent.emplace(caloClusterCollectionESPutToken_, std::move(caloClustersES));
0276
0277
0278 for (auto& sc : sClusters) {
0279 edm::Ptr<reco::CaloCluster> seedptr(caloClusHandleEBEE, pfClusterMapEBEE[sc.seed()]);
0280 sc.setSeed(seedptr);
0281
0282 reco::CaloClusterPtrVector clusters;
0283 for (reco::CaloCluster_iterator pfclus = sc.clustersBegin(); pfclus != sc.clustersEnd(); ++pfclus) {
0284 edm::Ptr<reco::CaloCluster> clusptr(caloClusHandleEBEE, pfClusterMapEBEE[*pfclus]);
0285 clusters.push_back(clusptr);
0286 }
0287 sc.setClusters(clusters);
0288
0289 reco::CaloClusterPtrVector psclusters;
0290 for (reco::CaloCluster_iterator pfclus = sc.preshowerClustersBegin(); pfclus != sc.preshowerClustersEnd();
0291 ++pfclus) {
0292 edm::Ptr<reco::CaloCluster> clusptr(caloClusHandleES, pfClusterMapES[*pfclus]);
0293 psclusters.push_back(clusptr);
0294 }
0295 sc.setPreshowerClusters(psclusters);
0296 }
0297
0298
0299 auto singleLegConv = createSingleLegConversions(egExtra, iEvent.getRefBeforePut<reco::ConversionCollection>());
0300
0301
0302 iEvent.emplace(superClusterCollectionPutToken_, std::move(sClusters));
0303 iEvent.emplace(pfCandidateEGammaExtraCollectionPutToken_, std::move(egExtra));
0304 iEvent.emplace(conversionCollectionPutToken_, std::move(singleLegConv));
0305 iEvent.emplace(pfCandidateCollectionPutToken_, std::move(egCandidates));
0306 }
0307
0308 reco::ConversionCollection PFEGammaProducer::createSingleLegConversions(
0309 reco::PFCandidateEGammaExtraCollection& extras, const edm::RefProd<reco::ConversionCollection>& convProd) const {
0310 reco::ConversionCollection oneLegConversions{};
0311 math::Error<3>::type error;
0312 for (auto& extra : extras) {
0313 for (const auto& tkrefmva : extra.singleLegConvTrackRefMva()) {
0314 const reco::Track& trk = *tkrefmva.first;
0315
0316 const reco::Vertex convVtx(trk.innerPosition(), error);
0317 std::vector<reco::TrackRef> OneLegConvVector;
0318 OneLegConvVector.push_back(tkrefmva.first);
0319 std::vector<float> OneLegMvaVector;
0320 OneLegMvaVector.push_back(tkrefmva.second);
0321 std::vector<reco::CaloClusterPtr> dummymatchingBC;
0322 reco::CaloClusterPtrVector scPtrVec;
0323 scPtrVec.push_back(edm::refToPtr(extra.superClusterRef()));
0324
0325 std::vector<math::XYZPointF> trackPositionAtEcalVec;
0326 std::vector<math::XYZPointF> innPointVec;
0327 std::vector<math::XYZVectorF> trackPinVec;
0328 std::vector<math::XYZVectorF> trackPoutVec;
0329 math::XYZPointF trackPositionAtEcal(trk.outerPosition().X(), trk.outerPosition().Y(), trk.outerPosition().Z());
0330 trackPositionAtEcalVec.push_back(trackPositionAtEcal);
0331
0332 math::XYZPointF innPoint(trk.innerPosition().X(), trk.innerPosition().Y(), trk.innerPosition().Z());
0333 innPointVec.push_back(innPoint);
0334
0335 math::XYZVectorF trackPin(trk.innerMomentum().X(), trk.innerMomentum().Y(), trk.innerMomentum().Z());
0336 trackPinVec.push_back(trackPin);
0337
0338 math::XYZVectorF trackPout(trk.outerMomentum().X(), trk.outerMomentum().Y(), trk.outerMomentum().Z());
0339 trackPoutVec.push_back(trackPout);
0340
0341 float DCA = trk.d0();
0342 float mvaval = tkrefmva.second;
0343 reco::Conversion singleLegConvCandidate(scPtrVec,
0344 OneLegConvVector,
0345 trackPositionAtEcalVec,
0346 convVtx,
0347 dummymatchingBC,
0348 DCA,
0349 innPointVec,
0350 trackPinVec,
0351 trackPoutVec,
0352 mvaval,
0353 reco::Conversion::pflow);
0354 singleLegConvCandidate.setOneLegMVA(OneLegMvaVector);
0355 oneLegConversions.push_back(singleLegConvCandidate);
0356
0357 reco::ConversionRef convref(convProd, oneLegConversions.size() - 1);
0358 extra.addSingleLegConversionRef(convref);
0359 }
0360 }
0361 return oneLegConversions;
0362 }
0363
0364 void PFEGammaProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0365 edm::ParameterSetDescription desc;
0366 desc.add<bool>("produceEGCandsWithNoSuperCluster", false)
0367 ->setComment("Allow building of candidates with no input or output supercluster?");
0368 desc.add<double>("pf_electron_mvaCut", -0.1);
0369 desc.add<bool>("pf_electronID_crackCorrection", false);
0370 desc.add<double>("pf_conv_mvaCut", 0.0);
0371 desc.add<edm::InputTag>("blocks", edm::InputTag("particleFlowBlock"))->setComment("PF Blocks label");
0372 desc.add<edm::InputTag>("EEtoPS_source", edm::InputTag("particleFlowClusterECAL"))
0373 ->setComment("EE to PS association");
0374 desc.add<edm::InputTag>("vertexCollection", edm::InputTag("offlinePrimaryVertices"));
0375 desc.add<edm::FileInPath>("pf_electronID_mvaWeightFile",
0376 edm::FileInPath("RecoParticleFlow/PFProducer/data/PfElectrons23Jan_BDT.weights.xml.gz"));
0377 desc.add<edm::FileInPath>("pf_convID_mvaWeightFile",
0378 edm::FileInPath("RecoParticleFlow/PFProducer/data/pfConversionAug0411_BDT.weights.xml.gz"));
0379 descriptions.add("particleFlowEGamma", desc);
0380 }
