File indexing completed on 2023-03-17 11:18:32
0001
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014
0015
0016
0017
0018
0019 #include "RecoJets/JetProducers/plugins/PileupJetIdProducer.h"
0020
0021 #include <memory>
0022
0023 GBRForestsAndConstants::GBRForestsAndConstants(edm::ParameterSet const& iConfig)
0024 : runMvas_(iConfig.getParameter<bool>("runMvas")),
0025 produceJetIds_(iConfig.getParameter<bool>("produceJetIds")),
0026 inputIsCorrected_(iConfig.getParameter<bool>("inputIsCorrected")),
0027 applyJec_(iConfig.getParameter<bool>("applyJec")),
0028 jec_(iConfig.getParameter<std::string>("jec")),
0029 residualsFromTxt_(iConfig.getParameter<bool>("residualsFromTxt")),
0030 applyConstituentWeight_(false) {
0031 if (residualsFromTxt_) {
0032 residualsTxt_ = iConfig.getParameter<edm::FileInPath>("residualsTxt");
0033 }
0034
0035 std::vector<edm::ParameterSet> algos = iConfig.getParameter<std::vector<edm::ParameterSet>>("algos");
0036 for (auto const& algoPset : algos) {
0037 vAlgoGBRForestsAndConstants_.emplace_back(algoPset, runMvas_);
0038 }
0039
0040 if (!runMvas_) {
0041 assert(algos.size() == 1);
0042 }
0043
0044 edm::InputTag srcConstituentWeights = iConfig.getParameter<edm::InputTag>("srcConstituentWeights");
0045 if (!srcConstituentWeights.label().empty()) {
0046 applyConstituentWeight_ = true;
0047 }
0048 }
0049
0050
0051 PileupJetIdProducer::PileupJetIdProducer(const edm::ParameterSet& iConfig, GBRForestsAndConstants const* globalCache) {
0052 if (globalCache->produceJetIds()) {
0053 produces<edm::ValueMap<StoredPileupJetIdentifier>>("");
0054 }
0055 for (auto const& algoGBRForestsAndConstants : globalCache->vAlgoGBRForestsAndConstants()) {
0056 std::string const& label = algoGBRForestsAndConstants.label();
0057 algos_.emplace_back(label, std::make_unique<PileupJetIdAlgo>(&algoGBRForestsAndConstants));
0058 if (globalCache->runMvas()) {
0059 produces<edm::ValueMap<float>>(label + "Discriminant");
0060 produces<edm::ValueMap<int>>(label + "Id");
0061 }
0062 }
0063
0064 input_jet_token_ = consumes<edm::View<reco::Jet>>(iConfig.getParameter<edm::InputTag>("jets"));
0065 input_vertex_token_ = consumes<reco::VertexCollection>(iConfig.getParameter<edm::InputTag>("vertexes"));
0066 input_vm_pujetid_token_ =
0067 consumes<edm::ValueMap<StoredPileupJetIdentifier>>(iConfig.getParameter<edm::InputTag>("jetids"));
0068 input_rho_token_ = consumes<double>(iConfig.getParameter<edm::InputTag>("rho"));
0069 parameters_token_ = esConsumes(edm::ESInputTag("", globalCache->jec()));
0070
0071 edm::InputTag srcConstituentWeights = iConfig.getParameter<edm::InputTag>("srcConstituentWeights");
0072 if (!srcConstituentWeights.label().empty()) {
0073 input_constituent_weights_token_ = consumes<edm::ValueMap<float>>(srcConstituentWeights);
0074 }
0075 }
0076
0077
0078 PileupJetIdProducer::~PileupJetIdProducer() {}
0079
0080
0081 void PileupJetIdProducer::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) {
0082 GBRForestsAndConstants const* gc = globalCache();
0083
0084 using namespace edm;
0085 using namespace std;
0086 using namespace reco;
0087
0088
0089 Handle<View<Jet>> jetHandle;
0090 iEvent.getByToken(input_jet_token_, jetHandle);
0091 const View<Jet>& jets = *jetHandle;
0092
0093
0094 edm::ValueMap<float> constituentWeights;
0095 if (!input_constituent_weights_token_.isUninitialized()) {
0096 constituentWeights = iEvent.get(input_constituent_weights_token_);
0097 }
0098
0099
0100 Handle<ValueMap<StoredPileupJetIdentifier>> vmap;
0101 if (!gc->produceJetIds()) {
0102 iEvent.getByToken(input_vm_pujetid_token_, vmap);
0103 }
0104
0105 edm::Handle<double> rhoH;
0106 double rho = 0.;
0107
0108
0109 vector<StoredPileupJetIdentifier> ids;
0110 map<string, vector<float>> mvas;
0111 map<string, vector<int>> idflags;
0112
0113 const VertexCollection* vertexes = nullptr;
0114 VertexCollection::const_iterator vtx;
0115 if (gc->produceJetIds()) {
0116
0117 Handle<VertexCollection> vertexHandle;
0118 iEvent.getByToken(input_vertex_token_, vertexHandle);
0119
0120 vertexes = vertexHandle.product();
0121
0122
0123 vtx = vertexes->begin();
0124 while (vtx != vertexes->end() && (vtx->isFake() || vtx->ndof() < 4)) {
0125 ++vtx;
0126 }
0127 if (vtx == vertexes->end()) {
0128 vtx = vertexes->begin();
0129 }
0130 }
0131
0132
0133 bool ispat = true;
0134 for (unsigned int i = 0; i < jets.size(); ++i) {
0135
0136 auto algoi = algos_.begin();
0137 PileupJetIdAlgo* ialgo = algoi->second.get();
0138
0139 const Jet& jet = jets.at(i);
0140 const pat::Jet* patjet = nullptr;
0141 if (ispat) {
0142 patjet = dynamic_cast<const pat::Jet*>(&jet);
0143 ispat = patjet != nullptr;
0144 }
0145
0146
0147 float jec = 0.;
0148 if (gc->applyJec()) {
0149
0150 if (rho == 0.) {
0151 iEvent.getByToken(input_rho_token_, rhoH);
0152 rho = *rhoH;
0153 }
0154
0155 if (not jecCor_) {
0156 initJetEnergyCorrector(iSetup, iEvent.isRealData());
0157 }
0158 if (ispat) {
0159 jecCor_->setJetPt(patjet->correctedJet(0).pt());
0160 } else {
0161 jecCor_->setJetPt(jet.pt());
0162 }
0163 jecCor_->setJetEta(jet.eta());
0164 jecCor_->setJetA(jet.jetArea());
0165 jecCor_->setRho(rho);
0166 jec = jecCor_->getCorrection();
0167 }
0168
0169 bool applyJec = gc->applyJec() && (ispat || !gc->inputIsCorrected());
0170 std::unique_ptr<reco::Jet> corrJet;
0171
0172 if (applyJec) {
0173 float scale = jec;
0174 if (ispat) {
0175 corrJet = std::make_unique<pat::Jet>(patjet->correctedJet(0));
0176 } else {
0177 corrJet.reset(dynamic_cast<reco::Jet*>(jet.clone()));
0178 }
0179 corrJet->scaleEnergy(scale);
0180 }
0181 const reco::Jet* theJet = (applyJec ? corrJet.get() : &jet);
0182
0183 PileupJetIdentifier puIdentifier;
0184 if (gc->produceJetIds()) {
0185
0186
0187 puIdentifier = ialgo->computeIdVariables(
0188 theJet, jec, &(*vtx), *vertexes, rho, constituentWeights, gc->applyConstituentWeight());
0189 ids.push_back(puIdentifier);
0190 } else {
0191
0192 puIdentifier = (*vmap)[jets.refAt(i)];
0193 puIdentifier.jetPt(theJet->pt());
0194 puIdentifier.jetEta(theJet->eta());
0195 puIdentifier.jetPhi(theJet->phi());
0196 ialgo->set(puIdentifier);
0197 puIdentifier = ialgo->computeMva();
0198 }
0199
0200 if (gc->runMvas()) {
0201
0202 mvas[algoi->first].push_back(puIdentifier.mva());
0203 idflags[algoi->first].push_back(puIdentifier.idFlag());
0204 for (++algoi; algoi != algos_.end(); ++algoi) {
0205 ialgo = algoi->second.get();
0206 ialgo->set(puIdentifier);
0207 PileupJetIdentifier id = ialgo->computeMva();
0208 mvas[algoi->first].push_back(id.mva());
0209 idflags[algoi->first].push_back(id.idFlag());
0210 }
0211 }
0212 }
0213
0214
0215 if (gc->runMvas()) {
0216 for (const auto& ialgo : algos_) {
0217
0218 vector<float>& mva = mvas[ialgo.first];
0219 auto mvaout = std::make_unique<ValueMap<float>>();
0220 ValueMap<float>::Filler mvafiller(*mvaout);
0221 mvafiller.insert(jetHandle, mva.begin(), mva.end());
0222 mvafiller.fill();
0223 iEvent.put(std::move(mvaout), ialgo.first + "Discriminant");
0224
0225
0226 vector<int>& idflag = idflags[ialgo.first];
0227 auto idflagout = std::make_unique<ValueMap<int>>();
0228 ValueMap<int>::Filler idflagfiller(*idflagout);
0229 idflagfiller.insert(jetHandle, idflag.begin(), idflag.end());
0230 idflagfiller.fill();
0231 iEvent.put(std::move(idflagout), ialgo.first + "Id");
0232 }
0233 }
0234
0235 if (gc->produceJetIds()) {
0236 assert(jetHandle->size() == ids.size());
0237 auto idsout = std::make_unique<ValueMap<StoredPileupJetIdentifier>>();
0238 ValueMap<StoredPileupJetIdentifier>::Filler idsfiller(*idsout);
0239 idsfiller.insert(jetHandle, ids.begin(), ids.end());
0240 idsfiller.fill();
0241 iEvent.put(std::move(idsout));
0242 }
0243 }
0244
0245
0246 void PileupJetIdProducer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) {
0247
0248
0249 edm::ParameterSetDescription desc;
0250 desc.setUnknown();
0251 descriptions.addDefault(desc);
0252 }
0253
0254
0255 void PileupJetIdProducer::initJetEnergyCorrector(const edm::EventSetup& iSetup, bool isData) {
0256 GBRForestsAndConstants const* gc = globalCache();
0257
0258
0259 std::vector<std::string> jecLevels;
0260 jecLevels.push_back("L1FastJet");
0261 jecLevels.push_back("L2Relative");
0262 jecLevels.push_back("L3Absolute");
0263 if (isData && !gc->residualsFromTxt())
0264 jecLevels.push_back("L2L3Residual");
0265
0266
0267 auto const& parameters = iSetup.getData(parameters_token_);
0268 for (std::vector<std::string>::const_iterator ll = jecLevels.begin(); ll != jecLevels.end(); ++ll) {
0269 const JetCorrectorParameters& ip = parameters[*ll];
0270 jetCorPars_.push_back(ip);
0271 }
0272 if (isData && gc->residualsFromTxt()) {
0273 jetCorPars_.push_back(JetCorrectorParameters(gc->residualsTxt().fullPath()));
0274 }
0275
0276
0277 jecCor_ = std::make_unique<FactorizedJetCorrector>(jetCorPars_);
0278 }
0279
0280 DEFINE_FWK_MODULE(PileupJetIdProducer);