Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:32:56 

0001 #include "DQMOffline/RecoB/interface/TagInfoPlotterFactory.h"
0002 #include "FWCore/Framework/interface/MakerMacros.h"
0003 #include "FWCore/Framework/interface/makeRefToBaseProdFrom.h"
0004 #include "Validation/RecoB/plugins/BTagPerformanceAnalyzerMC.h"
0005 
0006 using namespace reco;
0007 using namespace edm;
0008 using namespace std;
0009 using namespace RecoBTag;
0010 
0011 typedef std::pair<Jet, reco::JetFlavourInfo> JetWithFlavour;
0012 
0013 BTagPerformanceAnalyzerMC::BTagPerformanceAnalyzerMC(const edm::ParameterSet &pSet)
0014     : jetSelector(pSet.getParameter<double>("etaMin"),
0015                   pSet.getParameter<double>("etaMax"),
0016                   pSet.getParameter<double>("ptRecJetMin"),
0017                   pSet.getParameter<double>("ptRecJetMax"),
0018                   0.0,
0019                   99999.0,
0020                   pSet.getParameter<double>("ratioMin"),
0021                   pSet.getParameter<double>("ratioMax"),
0022                   pSet.getParameter<bool>("doJetID")),
0023       etaRanges(pSet.getParameter<vector<double>>("etaRanges")),
0024       ptRanges(pSet.getParameter<vector<double>>("ptRanges")),
0025       useOldFlavourTool(pSet.getParameter<bool>("useOldFlavourTool")),
0026       doJEC(pSet.getParameter<bool>("doJEC")),
0027       ptHatWeight(pSet.getParameter<bool>("applyPtHatWeight")),
0028       moduleConfig(pSet.getParameter<vector<edm::ParameterSet>>("tagConfig")),
0029       flavPlots_(pSet.getParameter<std::string>("flavPlots")),
0030       jetMatcher(pSet.getParameter<edm::ParameterSet>("recJetMatching")),
0031       doPUid(pSet.getParameter<bool>("doPUid")) {
0032   // mcPlots_ : 1=b+c+l+ni; 2=all+1; 3=1+d+u+s+g; 4=3+all . Default is 2. Don't
0033   // use 0.
0034   if (flavPlots_.find("dusg") < 15) {
0035     if (flavPlots_.find("all") < 15)
0036       mcPlots_ = 4;
0037     else
0038       mcPlots_ = 3;
0039   } else {
0040     if (flavPlots_.find("all") < 15)
0041       mcPlots_ = 2;
0042     else
0043       mcPlots_ = 1;
0044   }
0045   double ptRecJetMin = pSet.getParameter<double>("ptRecJetMin");
0046   jetMatcher.setThreshold(0.25 * ptRecJetMin);
0047   switch (pSet.getParameter<unsigned int>("leptonPlots")) {
0048     case 11:
0049       electronPlots = true;
0050       muonPlots = false;
0051       tauPlots = false;
0052       break;
0053     case 13:
0054       muonPlots = true;
0055       electronPlots = false;
0056       tauPlots = false;
0057       break;
0058     case 15:
0059       tauPlots = true;
0060       electronPlots = false;
0061       tauPlots = false;
0062       break;
0063     default:
0064       electronPlots = false;
0065       muonPlots = false;
0066       tauPlots = false;
0067   }
0068 
0069   if (etaRanges.size() <= 1)
0070     etaRanges = {pSet.getParameter<double>("etaMin"), pSet.getParameter<double>("etaMax")};
0071   if (ptRanges.size() <= 1)
0072     ptRanges = {pSet.getParameter<double>("ptRecJetMin"), pSet.getParameter<double>("ptRecJetMax")};
0073 
0074   genToken = mayConsume<GenEventInfoProduct>(edm::InputTag("generator"));
0075   genJetsMatchedToken =
0076       mayConsume<edm::Association<reco::GenJetCollection>>(pSet.getParameter<InputTag>("genJetsMatched"));
0077   jetToken = consumes<JetFlavourInfoMatchingCollection>(pSet.getParameter<InputTag>("jetMCSrc"));
0078   caloJetToken = mayConsume<JetFlavourMatchingCollection>(pSet.getParameter<InputTag>("caloJetMCSrc"));
0079   slInfoToken = consumes<SoftLeptonTagInfoCollection>(pSet.getParameter<InputTag>("softLeptonInfo"));
0080   jecMCToken = consumes<JetCorrector>(pSet.getParameter<edm::InputTag>("JECsourceMC"));
0081   jecDataToken = mayConsume<JetCorrector>(pSet.getParameter<edm::InputTag>("JECsourceData"));
0082 
0083   for (vector<edm::ParameterSet>::const_iterator iModule = moduleConfig.begin(); iModule != moduleConfig.end();
0084        ++iModule) {
0085     const string &dataFormatType = iModule->exists("type") ? iModule->getParameter<string>("type") : "JetTag";
0086     if (dataFormatType == "JetTag") {
0087       const InputTag &moduleLabel = iModule->getParameter<InputTag>("label");
0088       jetTagInputTags.push_back(moduleLabel);
0089       binJetTagPlotters.push_back(vector<std::unique_ptr<JetTagPlotter>>());
0090       jetTagToken.push_back(consumes<JetTagCollection>(moduleLabel));
0091     } else if (dataFormatType == "TagCorrelation") {
0092       const InputTag &label1 = iModule->getParameter<InputTag>("label1");
0093       const InputTag &label2 = iModule->getParameter<InputTag>("label2");
0094       tagCorrelationInputTags.push_back(std::pair<edm::InputTag, edm::InputTag>(label1, label2));
0095       binTagCorrelationPlotters.push_back(vector<std::unique_ptr<TagCorrelationPlotter>>());
0096       tagCorrelationToken.push_back(
0097           std::pair<edm::EDGetTokenT<reco::JetTagCollection>, edm::EDGetTokenT<reco::JetTagCollection>>(
0098               consumes<JetTagCollection>(label1), consumes<JetTagCollection>(label2)));
0099     } else {
0100       vector<edm::InputTag> vIP;
0101       tiDataFormatType.push_back(dataFormatType);
0102       binTagInfoPlotters.push_back(vector<std::unique_ptr<BaseTagInfoPlotter>>());
0103       std::vector<edm::EDGetTokenT<edm::View<reco::BaseTagInfo>>> tokens;
0104       if (dataFormatType == "GenericMVA") {
0105         const std::vector<InputTag> listInfo = iModule->getParameter<vector<InputTag>>("listTagInfos");
0106         for (unsigned int ITi = 0; ITi < listInfo.size(); ITi++) {
0107           tokens.push_back(consumes<View<BaseTagInfo>>(listInfo[ITi]));
0108           vIP.push_back(listInfo[ITi]);
0109         }
0110       } else {
0111         const InputTag &moduleLabel = iModule->getParameter<InputTag>("label");
0112         tokens.push_back(consumes<View<BaseTagInfo>>(moduleLabel));
0113         vIP.push_back(moduleLabel);
0114       }
0115       tagInfoToken.push_back(tokens);
0116       tagInfoInputTags.push_back(vIP);
0117     }
0118   }
0119 }
0120 
0121 void BTagPerformanceAnalyzerMC::bookHistograms(DQMStore::IBooker &ibook,
0122                                                edm::Run const &run,
0123                                                edm::EventSetup const &es) {
0124   //
0125   // Book all histograms.
0126   //
0127 
0128   // iterate over ranges:
0129   const int iEtaStart = -1;  // this will be the inactive one
0130   const int iEtaEnd = etaRanges.size() > 2 ? etaRanges.size() - 1
0131                                            : 0;  // if there is only one bin defined, leave it as the inactive one
0132   const int iPtStart = -1;                       // this will be the inactive one
0133   const int iPtEnd =
0134       ptRanges.size() > 2 ? ptRanges.size() - 1 : 0;  // if there is only one bin defined, leave it as the inactive one
0135 
0136   setTDRStyle();
0137 
0138   TagInfoPlotterFactory theFactory;
0139   int iTag = -1;
0140   int iTagCorr = -1;
0141   int iInfoTag = -1;
0142   for (vector<edm::ParameterSet>::const_iterator iModule = moduleConfig.begin(); iModule != moduleConfig.end();
0143        ++iModule) {
0144     const string &dataFormatType = iModule->exists("type") ? iModule->getParameter<string>("type") : "JetTag";
0145     if (dataFormatType == "JetTag") {
0146       iTag++;
0147       const string &folderName = iModule->getParameter<string>("folder");
0148 
0149       bool doDifferentialPlots = false;
0150       double discrCut = -999.;
0151       if (iModule->exists("differentialPlots") && iModule->getParameter<bool>("differentialPlots") == true) {
0152         doDifferentialPlots = true;
0153         discrCut = iModule->getParameter<double>("discrCut");
0154       }
0155 
0156       // eta loop
0157       for (int iEta = iEtaStart; iEta < iEtaEnd; iEta++) {
0158         // pt loop
0159         for (int iPt = iPtStart; iPt < iPtEnd; iPt++) {
0160           const EtaPtBin &etaPtBin = getEtaPtBin(iEta, iPt);
0161 
0162           // Instantiate the genertic b tag plotter
0163           binJetTagPlotters.at(iTag).push_back(
0164               std::make_unique<JetTagPlotter>(folderName,
0165                                               etaPtBin,
0166                                               iModule->getParameter<edm::ParameterSet>("parameters"),
0167                                               mcPlots_,
0168                                               false,
0169                                               ibook,
0170                                               false,
0171                                               doDifferentialPlots,
0172                                               discrCut));
0173         }
0174       }
0175     } else if (dataFormatType == "TagCorrelation") {
0176       iTagCorr++;
0177       const InputTag &label1 = iModule->getParameter<InputTag>("label1");
0178       const InputTag &label2 = iModule->getParameter<InputTag>("label2");
0179 
0180       // eta loop
0181       for (int iEta = iEtaStart; iEta != iEtaEnd; ++iEta) {
0182         // pt loop
0183         for (int iPt = iPtStart; iPt != iPtEnd; ++iPt) {
0184           const EtaPtBin &etaPtBin = getEtaPtBin(iEta, iPt);
0185           // Instantiate the generic b tag correlation plotter
0186           binTagCorrelationPlotters.at(iTagCorr).push_back(
0187               std::make_unique<TagCorrelationPlotter>(label1.label(),
0188                                                       label2.label(),
0189                                                       etaPtBin,
0190                                                       iModule->getParameter<edm::ParameterSet>("parameters"),
0191                                                       mcPlots_,
0192                                                       false,
0193                                                       false,
0194                                                       ibook));
0195         }
0196       }
0197     } else {
0198       iInfoTag++;
0199       // tag info retrievel is deferred (needs availability of EventSetup)
0200       const InputTag &moduleLabel = iModule->getParameter<InputTag>("label");
0201       const string &folderName = iModule->getParameter<string>("folder");
0202       // eta loop
0203       for (int iEta = iEtaStart; iEta < iEtaEnd; iEta++) {
0204         // pt loop
0205         for (int iPt = iPtStart; iPt < iPtEnd; iPt++) {
0206           const EtaPtBin &etaPtBin = getEtaPtBin(iEta, iPt);
0207 
0208           // Instantiate the tagInfo plotter
0209           binTagInfoPlotters.at(iInfoTag).push_back(
0210               theFactory.buildPlotter(dataFormatType,
0211                                       moduleLabel.label(),
0212                                       etaPtBin,
0213                                       iModule->getParameter<edm::ParameterSet>("parameters"),
0214                                       folderName,
0215                                       mcPlots_,
0216                                       false,
0217                                       ibook));
0218         }
0219       }
0220     }
0221   }
0222 }
0223 
0224 EtaPtBin BTagPerformanceAnalyzerMC::getEtaPtBin(const int &iEta, const int &iPt) {
0225   // DEFINE BTagBin:
0226   bool etaActive_, ptActive_;
0227   double etaMin_, etaMax_, ptMin_, ptMax_;
0228 
0229   if (iEta != -1) {
0230     etaActive_ = true;
0231     etaMin_ = etaRanges[iEta];
0232     etaMax_ = etaRanges[iEta + 1];
0233   } else {
0234     etaActive_ = false;
0235     etaMin_ = etaRanges[0];
0236     etaMax_ = etaRanges[etaRanges.size() - 1];
0237   }
0238 
0239   if (iPt != -1) {
0240     ptActive_ = true;
0241     ptMin_ = ptRanges[iPt];
0242     ptMax_ = ptRanges[iPt + 1];
0243   } else {
0244     ptActive_ = false;
0245     ptMin_ = ptRanges[0];
0246     ptMax_ = ptRanges[ptRanges.size() - 1];
0247   }
0248   return EtaPtBin(etaActive_, etaMin_, etaMax_, ptActive_, ptMin_, ptMax_);
0249 }
0250 
0251 BTagPerformanceAnalyzerMC::~BTagPerformanceAnalyzerMC() {}
0252 
0253 void BTagPerformanceAnalyzerMC::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0254   eventInitialized = false;
0255 
0256   float weight = 1;  // event weight
0257 
0258   if (ptHatWeight) {
0259     /* APPLY PTHAT EVENT  WEIGHT */
0260 
0261     edm::Handle<GenEventInfoProduct> genInfoHandle;
0262     iEvent.getByToken(genToken, genInfoHandle);
0263 
0264     if (genInfoHandle.isValid()) {
0265       weight = weight * static_cast<float>(genInfoHandle->weight());
0266     }
0267   }
0268 
0269   LogDebug("Info") << "Event weight is: " << weight;
0270 
0271   FlavourMap flavours;
0272   LeptonMap leptons;
0273 
0274   if (!useOldFlavourTool) {
0275     edm::Handle<JetFlavourInfoMatchingCollection> jetMC;
0276     iEvent.getByToken(jetToken, jetMC);
0277     for (JetFlavourInfoMatchingCollection::const_iterator iter = jetMC->begin(); iter != jetMC->end(); ++iter) {
0278       unsigned int fl = std::abs(iter->second.getPartonFlavour());
0279       flavours.insert(std::make_pair(iter->first, fl));
0280       const GenParticleRefVector &lep = iter->second.getLeptons();
0281       reco::JetFlavour::Leptons lepCount;
0282       for (unsigned int i = 0; i < lep.size(); i++) {
0283         if (abs(lep[i]->pdgId()) == 11)
0284           lepCount.electron++;
0285         else if (abs(lep[i]->pdgId()) == 13)
0286           lepCount.muon++;
0287         else if (abs(lep[i]->pdgId()) == 15)
0288           lepCount.tau++;
0289       }
0290       leptons.insert(std::make_pair(iter->first, lepCount));
0291     }
0292   } else {
0293     edm::Handle<JetFlavourMatchingCollection> jetMC;
0294     iEvent.getByToken(caloJetToken, jetMC);
0295     for (JetFlavourMatchingCollection::const_iterator iter = jetMC->begin(); iter != jetMC->end(); ++iter) {
0296       unsigned int fl = std::abs(iter->second.getFlavour());
0297       flavours.insert(std::make_pair(iter->first, fl));
0298       const reco::JetFlavour::Leptons &lep = iter->second.getLeptons();
0299       leptons.insert(std::make_pair(iter->first, lep));
0300     }
0301   }
0302 
0303   edm::Handle<reco::SoftLeptonTagInfoCollection> infoHandle;
0304   iEvent.getByToken(slInfoToken, infoHandle);
0305 
0306   edm::Handle<edm::Association<reco::GenJetCollection>> genJetsMatched;
0307   if (doPUid) {
0308     iEvent.getByToken(genJetsMatchedToken, genJetsMatched);
0309   }
0310 
0311   // Get JEC
0312   const JetCorrector *corrector = nullptr;
0313   if (doJEC) {
0314     edm::Handle<GenEventInfoProduct> genInfoHandle;  // check if data or MC
0315     iEvent.getByToken(genToken, genInfoHandle);
0316     edm::Handle<JetCorrector> corrHandle;
0317     if (!genInfoHandle.isValid())
0318       iEvent.getByToken(jecDataToken, corrHandle);
0319     else
0320       iEvent.getByToken(jecMCToken, corrHandle);
0321     corrector = corrHandle.product();
0322   }
0323   //
0324 
0325   // Look first at the jetTags
0326   for (unsigned int iJetLabel = 0; iJetLabel != jetTagInputTags.size(); ++iJetLabel) {
0327     edm::Handle<reco::JetTagCollection> tagHandle;
0328     iEvent.getByToken(jetTagToken[iJetLabel], tagHandle);
0329     const reco::JetTagCollection &tagColl = *(tagHandle.product());
0330     LogDebug("Info") << "Found " << tagColl.size() << " B candidates in collection " << jetTagInputTags[iJetLabel];
0331 
0332     int plotterSize = binJetTagPlotters[iJetLabel].size();
0333     for (JetTagCollection::const_iterator tagI = tagColl.begin(); tagI != tagColl.end(); ++tagI) {
0334       // Identify parton associated to jet.
0335 
0336       /// needed for lepton specific plots
0337       if (flavours[tagI->first] == 5 &&
0338           ((electronPlots && !leptons[tagI->first].electron) || (muonPlots && !leptons[tagI->first].muon) ||
0339            (tauPlots && !leptons[tagI->first].tau)))
0340         continue;
0341       // JEC
0342       double jec = 1.0;
0343       /*reco::Jet correctedJet = *(tagI->first);
0344       if(doJEC && corrector) {
0345         jec = corrector->correction(*(tagI->first));
0346       }*/
0347 
0348       JetWithFlavour jetWithFlavour;
0349       // also applies JEC to jet
0350       if (!getJetWithFlavour(iEvent, tagI->first, flavours, jetWithFlavour, corrector, genJetsMatched))
0351         continue;
0352       if (!jetSelector(jetWithFlavour.first, std::abs(jetWithFlavour.second.getPartonFlavour()), infoHandle, jec))
0353         continue;
0354 
0355       for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter) {
0356         bool inBin = false;
0357         inBin = binJetTagPlotters[iJetLabel][iPlotter]->etaPtBin().inBin(jetWithFlavour.first, jec);
0358         // Fill histograms if in desired pt/rapidity bin.
0359         if (inBin)
0360           binJetTagPlotters[iJetLabel][iPlotter]->analyzeTag(
0361               jetWithFlavour.first, jec, tagI->second, std::abs(jetWithFlavour.second.getPartonFlavour()), weight);
0362       }
0363     }
0364     for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter) {
0365       binJetTagPlotters[iJetLabel][iPlotter]->analyzeTag(weight);
0366     }
0367   }
0368 
0369   // Now look at Tag Correlations
0370   for (unsigned int iJetLabel = 0; iJetLabel != tagCorrelationInputTags.size(); ++iJetLabel) {
0371     const std::pair<edm::EDGetTokenT<reco::JetTagCollection>, edm::EDGetTokenT<reco::JetTagCollection>> &inputTokens =
0372         tagCorrelationToken[iJetLabel];
0373     edm::Handle<reco::JetTagCollection> tagHandle1;
0374     iEvent.getByToken(inputTokens.first, tagHandle1);
0375     const reco::JetTagCollection &tagColl1 = *(tagHandle1.product());
0376 
0377     edm::Handle<reco::JetTagCollection> tagHandle2;
0378     iEvent.getByToken(inputTokens.second, tagHandle2);
0379     const reco::JetTagCollection &tagColl2 = *(tagHandle2.product());
0380 
0381     int plotterSize = binTagCorrelationPlotters[iJetLabel].size();
0382     for (JetTagCollection::const_iterator tagI = tagColl1.begin(); tagI != tagColl1.end(); ++tagI) {
0383       if (flavours[tagI->first] == 5 &&
0384           ((electronPlots && !leptons[tagI->first].electron) || (muonPlots && !leptons[tagI->first].muon) ||
0385            (tauPlots && !leptons[tagI->first].tau)))
0386         continue;
0387 
0388       // JEC
0389       double jec = 1.0;
0390       /*reco::Jet correctedJet = *(tagI->first);
0391       if(doJEC && corrector) {
0392         jec = corrector->correction(*(tagI->first));
0393       }*/
0394 
0395       JetWithFlavour jetWithFlavour;
0396       // also applies JEC to jet
0397       if (!getJetWithFlavour(iEvent, tagI->first, flavours, jetWithFlavour, corrector, genJetsMatched))
0398         continue;
0399       if (!jetSelector(jetWithFlavour.first, std::abs(jetWithFlavour.second.getPartonFlavour()), infoHandle, jec))
0400         continue;
0401 
0402       for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter) {
0403         bool inBin = false;
0404         inBin = binTagCorrelationPlotters[iJetLabel][iPlotter]->etaPtBin().inBin(jetWithFlavour.first, jec);
0405 
0406         if (inBin) {
0407           double discr2 = tagColl2[tagI->first];
0408           binTagCorrelationPlotters[iJetLabel][iPlotter]->analyzeTags(
0409               tagI->second, discr2, std::abs(jetWithFlavour.second.getPartonFlavour()), weight);
0410         }
0411       }
0412     }
0413   }
0414 
0415   // Now look at the TagInfos
0416 
0417   for (unsigned int iJetLabel = 0; iJetLabel != tiDataFormatType.size(); ++iJetLabel) {
0418     int plotterSize = binTagInfoPlotters[iJetLabel].size();
0419     for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter)
0420       binTagInfoPlotters[iJetLabel][iPlotter]->setEventSetup(iSetup);
0421 
0422     vector<edm::EDGetTokenT<edm::View<reco::BaseTagInfo>>> &tokens = tagInfoToken[iJetLabel];
0423     // check number of tag infos = expected number of tag infos
0424     vector<string> labels = binTagInfoPlotters[iJetLabel][0]->tagInfoRequirements();
0425     if (labels.empty())
0426       labels.push_back("label");
0427     if (labels.size() != tokens.size())
0428       throw cms::Exception("Configuration")
0429           << "Different number of Tag Infos than expected" << labels.size() << tokens.size() << endl;
0430 
0431     unsigned int nInputTags = tokens.size();
0432     vector<edm::Handle<View<BaseTagInfo>>> tagInfoHandles(nInputTags);
0433     edm::ProductID jetProductID;
0434     unsigned int nTagInfos = 0;
0435     for (unsigned int iInputTags = 0; iInputTags < tokens.size(); ++iInputTags) {
0436       edm::Handle<View<BaseTagInfo>> &tagInfoHandle = tagInfoHandles[iInputTags];
0437       iEvent.getByToken(tokens[iInputTags], tagInfoHandle);
0438       if (tagInfoHandle.isValid() == false) {
0439         edm::LogWarning("BTagPerformanceAnalyzerMC")
0440             << " Collection " << tagInfoInputTags[iJetLabel][iInputTags] << " not present. Skipping it for this event.";
0441         continue;
0442       }
0443 
0444       unsigned int size = tagInfoHandle->size();
0445       LogDebug("Info") << "Found " << size << " B candidates in collection " << tagInfoInputTags[iJetLabel][iInputTags];
0446       edm::ProductID thisProductID = (size > 0) ? (*tagInfoHandle)[0].jet().id() : edm::ProductID();
0447       if (iInputTags == 0) {
0448         jetProductID = thisProductID;
0449         nTagInfos = size;
0450       } else if (jetProductID != thisProductID)
0451         throw cms::Exception("Configuration") << "TagInfos are referencing a different jet collection." << endl;
0452       else if (nTagInfos != size)
0453         throw cms::Exception("Configuration") << "TagInfo collections are having a different size." << endl;
0454     }
0455 
0456     for (unsigned int iTagInfos = 0; iTagInfos < nTagInfos; ++iTagInfos) {
0457       vector<const BaseTagInfo *> baseTagInfos(nInputTags);
0458       edm::RefToBase<Jet> jetRef;
0459       for (unsigned int iTagInfo = 0; iTagInfo < nInputTags; iTagInfo++) {
0460         const BaseTagInfo &baseTagInfo = (*tagInfoHandles[iTagInfo])[iTagInfos];
0461         if (iTagInfo == 0)
0462           jetRef = baseTagInfo.jet();
0463         else if (baseTagInfo.jet() != jetRef)
0464           throw cms::Exception("Configuration") << "TagInfos pointing to different jets." << endl;
0465         baseTagInfos[iTagInfo] = &baseTagInfo;
0466       }
0467 
0468       // Identify parton associated to jet.
0469 
0470       /// needed for lepton specific plots
0471       if (flavours[jetRef] == 5 && ((electronPlots && !leptons[jetRef].electron) ||
0472                                     (muonPlots && !leptons[jetRef].muon) || (tauPlots && !leptons[jetRef].tau)))
0473         continue;
0474 
0475       // JEC
0476       double jec = 1.0;
0477       /*reco::Jet correctedJet = *(jetRef);
0478       if(doJEC && corrector) {
0479         jec = corrector->correction(*(jetRef));
0480       }*/
0481 
0482       JetWithFlavour jetWithFlavour;
0483       // also applies JEC to jet
0484       if (!getJetWithFlavour(iEvent, jetRef, flavours, jetWithFlavour, corrector, genJetsMatched))
0485         continue;
0486       if (!jetSelector(jetWithFlavour.first, std::abs(jetWithFlavour.second.getPartonFlavour()), infoHandle, jec))
0487         continue;
0488 
0489       for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter) {
0490         bool inBin = false;
0491         inBin = binTagInfoPlotters[iJetLabel][iPlotter]->etaPtBin().inBin(*jetRef, jec);
0492         // Fill histograms if in desired pt/rapidity bin.
0493         if (inBin)
0494           binTagInfoPlotters[iJetLabel][iPlotter]->analyzeTag(
0495               baseTagInfos, jec, std::abs(jetWithFlavour.second.getPartonFlavour()), weight);
0496       }
0497     }
0498   }
0499 }
0500 
0501 bool BTagPerformanceAnalyzerMC::getJetWithGenJet(edm::RefToBase<Jet> jetRef,
0502                                                  edm::Handle<edm::Association<reco::GenJetCollection>> genJetsMatched) {
0503   if (!doPUid)
0504     return true;
0505   reco::GenJetRef genjet = (*genJetsMatched)[jetRef];
0506   if (genjet.isNonnull() && genjet.isAvailable())
0507     return true;
0508   return false;
0509 }
0510 
0511 bool BTagPerformanceAnalyzerMC::getJetWithFlavour(const edm::Event &iEvent,
0512                                                   edm::RefToBase<Jet> jetRef,
0513                                                   const FlavourMap &flavours,
0514                                                   JetWithFlavour &jetWithFlavour,
0515                                                   const JetCorrector *corrector,
0516                                                   edm::Handle<edm::Association<reco::GenJetCollection>> genJetsMatched) {
0517   edm::ProductID recProdId = jetRef.id();
0518   edm::ProductID refProdId = (flavours.begin() == flavours.end()) ? recProdId : flavours.begin()->first.id();
0519 
0520   if (!eventInitialized) {
0521     jetCorrector.setCorrector(corrector);
0522     if (recProdId != refProdId) {
0523       edm::RefToBaseVector<Jet> refJets;
0524       for (FlavourMap::const_iterator iter = flavours.begin(); iter != flavours.end(); ++iter)
0525         refJets.push_back(iter->first);
0526       const edm::RefToBaseProd<Jet> recJetsProd(edm::makeRefToBaseProdFrom(jetRef, iEvent));
0527       edm::RefToBaseVector<Jet> recJets;
0528       for (unsigned int i = 0; i < recJetsProd->size(); i++)
0529         recJets.push_back(edm::RefToBase<Jet>(recJetsProd, i));
0530       jetMatcher.matchCollections(refJets, recJets, corrector);
0531     }
0532     eventInitialized = true;
0533   }
0534 
0535   if (recProdId != refProdId) {
0536     jetRef = jetMatcher(jetRef);
0537     if (jetRef.isNull())
0538       return false;
0539   }
0540 
0541   // apply JEC
0542   jetWithFlavour.first = jetCorrector(*jetRef);
0543 
0544   auto itFound = flavours.find(jetRef);
0545   unsigned int flavour = itFound != flavours.end() ? itFound->second : 0;
0546 
0547   if (doPUid) {
0548     bool isNotPU = getJetWithGenJet(jetRef, genJetsMatched);
0549     if (!isNotPU)
0550       flavour = 20;
0551   }
0552 
0553   jetWithFlavour.second = reco::JetFlavourInfo(flavour, flavour);
0554 
0555   LogTrace("Info") << "Found jet with flavour " << jetWithFlavour.second.getPartonFlavour() << endl;
0556   LogTrace("Info") << jetWithFlavour.first.p() << " , " << jetWithFlavour.first.pt() << " - " << endl;
0557   //<< jetWithFlavour.second.getLorentzVector().P()<<" , "<<
0558   // jetWithFlavour.second.getLorentzVector().Pt()<<endl;
0559 
0560   return true;
0561 }
0562 
0563 // define this as a plug-in
0564 DEFINE_FWK_MODULE(BTagPerformanceAnalyzerMC);
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!