PFProducer/plugins/PFEGammaProducer.cc

0001 /**\class PFEGammaProducer
0002 \brief Producer for particle flow reconstructed particles (PFCandidates)
0003
0004 This producer makes use of PFAlgo, the particle flow algorithm.
0005
0006 \author Colin Bernet
0007 \date   July 2006
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   /// particle flow algorithm configuration
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   // output collections
0103   reco::PFCandidateCollection egCandidates{};
0104   reco::PFCandidateEGammaExtraCollection egExtra{};
0105   reco::SuperClusterCollection sClusters{};
0106
0107   // preshower conditions
0108   auto esEEInterCalibHandle_ = iSetup.getHandle(esEEInterCalibToken_);
0109   auto esChannelStatusHandle_ = iSetup.getHandle(esChannelStatusToken_);
0110
0111   //Assign the PFAlgo Parameters
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   // get the collection of blocks
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   //str<<(*pfAlgo_)<<std::endl;
0139   //    cout << (*pfAlgo_) << std::endl;
0140   LOGDRESSED("PFEGammaProducer") << str.str() << std::endl;
0141 #endif
0142
0143   // sort elements in three lists:
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           // Single HO elements are likely to be noise. Not considered for now.
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   }  //loop blocks
0190
0191   // loop on blocks that are not single ecal, single ps1, single ps2 , or
0192   // single hcal and produce unbiased collection of EGamma Candidates
0193
0194   //printf("loop over blocks\n");
0195   unsigned nblcks = 0;
0196
0197   // this auto is a const reco::PFBlockRef&
0198   for (const auto& blockref : otherBlockRefs) {
0199     ++nblcks;
0200     // this auto is a: const edm::OwnVector< reco::PFBlockElement >&
0201     const auto& elements = blockref->elements();
0202     // make a copy of the link data, which will be edited.
0203     //PFBlock::LinkData linkData =  block.linkData();
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   //set the correct references to refined SC and EG extra using the refprods
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   //build collections of output CaloClusters from the used PFClusters
0245   reco::CaloClusterCollection caloClustersEBEE{};
0246   reco::CaloClusterCollection caloClustersES{};
0247
0248   std::map<edm::Ptr<reco::CaloCluster>, unsigned int> pfClusterMapEBEE;  //maps of pfclusters to caloclusters
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   //put calocluster output collections in event and get orphan handles to create ptrs
0276   auto const& caloClusHandleEBEE = iEvent.emplace(caloClusterCollectionEBEEPutToken_, std::move(caloClustersEBEE));
0277   auto const& caloClusHandleES = iEvent.emplace(caloClusterCollectionESPutToken_, std::move(caloClustersES));
0278
0279   //relink superclusters to output caloclusters
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   //create and fill references to single leg conversions
0301   auto singleLegConv = createSingleLegConversions(egExtra, iEvent.getRefBeforePut<reco::ConversionCollection>());
0302
0303   // release our demonspawn into the wild to cause havoc
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 }