Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-10-25 10:01:45

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 
0196   // this auto is a const reco::PFBlockRef&
0197   for (const auto& blockref : otherBlockRefs) {
0198     // this auto is a: const edm::OwnVector< reco::PFBlockElement >&
0199     const auto& elements = blockref->elements();
0200     // make a copy of the link data, which will be edited.
0201     //PFBlock::LinkData linkData =  block.linkData();
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   //set the correct references to refined SC and EG extra using the refprods
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   //build collections of output CaloClusters from the used PFClusters
0243   reco::CaloClusterCollection caloClustersEBEE{};
0244   reco::CaloClusterCollection caloClustersES{};
0245 
0246   std::map<edm::Ptr<reco::CaloCluster>, unsigned int> pfClusterMapEBEE;  //maps of pfclusters to caloclusters
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   //put calocluster output collections in event and get orphan handles to create ptrs
0274   auto const& caloClusHandleEBEE = iEvent.emplace(caloClusterCollectionEBEEPutToken_, std::move(caloClustersEBEE));
0275   auto const& caloClusHandleES = iEvent.emplace(caloClusterCollectionESPutToken_, std::move(caloClustersES));
0276 
0277   //relink superclusters to output caloclusters
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   //create and fill references to single leg conversions
0299   auto singleLegConv = createSingleLegConversions(egExtra, iEvent.getRefBeforePut<reco::ConversionCollection>());
0300 
0301   // release our demonspawn into the wild to cause havoc
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 }