RecoB/plugins/BTagPerformanceAnalyzerOnData.cc

0001 #include "DQMOffline/RecoB/plugins/BTagPerformanceAnalyzerOnData.h"
0002 #include "FWCore/Framework/interface/MakerMacros.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "DQMOffline/RecoB/interface/TagInfoPlotterFactory.h"
0005
0006 using namespace reco;
0007 using namespace edm;
0008 using namespace std;
0009 using namespace RecoBTag;
0010
0011 BTagPerformanceAnalyzerOnData::BTagPerformanceAnalyzerOnData(const edm::ParameterSet& pSet)
0012     : jetSelector(pSet.getParameter<double>("etaMin"),
0013                   pSet.getParameter<double>("etaMax"),
0014                   pSet.getParameter<double>("ptRecJetMin"),
0015                   pSet.getParameter<double>("ptRecJetMax"),
0016                   0.0,
0017                   99999.0,
0018                   pSet.getParameter<double>("ratioMin"),
0019                   pSet.getParameter<double>("ratioMax"),
0020                   pSet.getParameter<bool>("doJetID")),
0021       etaRanges(pSet.getParameter<vector<double>>("etaRanges")),
0022       ptRanges(pSet.getParameter<vector<double>>("ptRanges")),
0023       doJEC(pSet.getParameter<bool>("doJEC")),
0024       moduleConfig(pSet.getParameter<vector<edm::ParameterSet>>("tagConfig")) {
0025   genToken = mayConsume<GenEventInfoProduct>(edm::InputTag("generator"));
0026   slInfoToken = consumes<SoftLeptonTagInfoCollection>(pSet.getParameter<InputTag>("softLeptonInfo"));
0027   jecMCToken = mayConsume<JetCorrector>(pSet.getParameter<edm::InputTag>("JECsourceMC"));
0028   jecDataToken = consumes<JetCorrector>(pSet.getParameter<edm::InputTag>("JECsourceData"));
0029
0030   if (etaRanges.size() <= 1)
0031     etaRanges = {pSet.getParameter<double>("etaMin"), pSet.getParameter<double>("etaMax")};
0032   if (ptRanges.size() <= 1)
0033     ptRanges = {pSet.getParameter<double>("ptRecJetMin"), pSet.getParameter<double>("ptRecJetMax")};
0034
0035   for (vector<edm::ParameterSet>::const_iterator iModule = moduleConfig.begin(); iModule != moduleConfig.end();
0036        ++iModule) {
0037     const string& dataFormatType = iModule->exists("type") ? iModule->getParameter<string>("type") : "JetTag";
0038     if (dataFormatType == "JetTag") {
0039       const InputTag& moduleLabel = iModule->getParameter<InputTag>("label");
0040       jetTagInputTags.push_back(moduleLabel);
0041       binJetTagPlotters.push_back(vector<std::unique_ptr<JetTagPlotter>>());
0042       jetTagToken.push_back(consumes<JetTagCollection>(moduleLabel));
0043     } else if (dataFormatType == "TagCorrelation") {
0044       const InputTag& label1 = iModule->getParameter<InputTag>("label1");
0045       const InputTag& label2 = iModule->getParameter<InputTag>("label2");
0046       tagCorrelationInputTags.push_back(std::pair<edm::InputTag, edm::InputTag>(label1, label2));
0047       binTagCorrelationPlotters.push_back(vector<std::unique_ptr<TagCorrelationPlotter>>());
0048       tagCorrelationToken.push_back(
0049           std::pair<edm::EDGetTokenT<reco::JetTagCollection>, edm::EDGetTokenT<reco::JetTagCollection>>(
0050               consumes<JetTagCollection>(label1), consumes<JetTagCollection>(label2)));
0051     } else {
0052       vector<edm::InputTag> vIP;
0053       tiDataFormatType.push_back(dataFormatType);
0054       binTagInfoPlotters.push_back(vector<std::unique_ptr<BaseTagInfoPlotter>>());
0055       std::vector<edm::EDGetTokenT<edm::View<reco::BaseTagInfo>>> tokens;
0056       if (dataFormatType == "GenericMVA") {
0057         const std::vector<InputTag> listInfo = iModule->getParameter<vector<InputTag>>("listTagInfos");
0058         for (unsigned int ITi = 0; ITi < listInfo.size(); ITi++) {
0059           tokens.push_back(consumes<View<BaseTagInfo>>(listInfo[ITi]));
0060           vIP.push_back(listInfo[ITi]);
0061         }
0062       } else {
0063         const InputTag& moduleLabel = iModule->getParameter<InputTag>("label");
0064         tokens.push_back(consumes<View<BaseTagInfo>>(moduleLabel));
0065         vIP.push_back(moduleLabel);
0066       }
0067       tagInfoToken.push_back(tokens);
0068       tagInfoInputTags.push_back(vIP);
0069     }
0070   }
0071 }
0072
0073 void BTagPerformanceAnalyzerOnData::bookHistograms(DQMStore::IBooker& ibook,
0074                                                    edm::Run const& run,
0075                                                    edm::EventSetup const& es) {
0076   // Book all histograms.
0077
0078   // iterate over ranges:
0079   const int iEtaStart = -1;  // this will be the inactive one
0080   const int iEtaEnd = etaRanges.size() > 2 ? etaRanges.size() - 1
0081                                            : 0;  // if there is only one bin defined, leave it as the inactive one
0082   const int iPtStart = -1;                       // this will be the inactive one
0083   const int iPtEnd =
0084       ptRanges.size() > 2 ? ptRanges.size() - 1 : 0;  // if there is only one bin defined, leave it as the inactive one
0085   setTDRStyle();
0086
0087   TagInfoPlotterFactory theFactory;
0088   int iTag = -1;
0089   int iTagCorr = -1;
0090   int iInfoTag = -1;
0091   for (vector<edm::ParameterSet>::const_iterator iModule = moduleConfig.begin(); iModule != moduleConfig.end();
0092        ++iModule) {
0093     const string& dataFormatType = iModule->exists("type") ? iModule->getParameter<string>("type") : "JetTag";
0094     if (dataFormatType == "JetTag") {
0095       iTag++;
0096       const string& folderName = iModule->getParameter<string>("folder");
0097
0098       bool doDifferentialPlots = false;
0099       double discrCut = -999.;
0100       if (iModule->exists("differentialPlots") && iModule->getParameter<bool>("differentialPlots") == true) {
0101         doDifferentialPlots = true;
0102         discrCut = iModule->getParameter<double>("discrCut");
0103       }
0104
0105       // eta loop
0106       for (int iEta = iEtaStart; iEta < iEtaEnd; ++iEta) {
0107         // pt loop
0108         for (int iPt = iPtStart; iPt < iPtEnd; ++iPt) {
0109           const EtaPtBin& etaPtBin = getEtaPtBin(iEta, iPt);
0110
0111           // Instantiate the generic b tag plotter
0112           binJetTagPlotters.at(iTag).push_back(
0113               std::make_unique<JetTagPlotter>(folderName,
0114                                               etaPtBin,
0115                                               iModule->getParameter<edm::ParameterSet>("parameters"),
0116                                               0,
0117                                               false,
0118                                               ibook,
0119                                               false,
0120                                               doDifferentialPlots,
0121                                               discrCut));
0122         }
0123       }
0124
0125     } else if (dataFormatType == "TagCorrelation") {
0126       iTagCorr++;
0127       const InputTag& label1 = iModule->getParameter<InputTag>("label1");
0128       const InputTag& label2 = iModule->getParameter<InputTag>("label2");
0129
0130       // eta loop
0131       for (int iEta = iEtaStart; iEta != iEtaEnd; ++iEta) {
0132         // pt loop
0133         for (int iPt = iPtStart; iPt != iPtEnd; ++iPt) {
0134           const EtaPtBin& etaPtBin = getEtaPtBin(iEta, iPt);
0135           // Instantiate the generic b tag correlation plotter
0136           binTagCorrelationPlotters.at(iTagCorr).push_back(
0137               std::make_unique<TagCorrelationPlotter>(label1.label(),
0138                                                       label2.label(),
0139                                                       etaPtBin,
0140                                                       iModule->getParameter<edm::ParameterSet>("parameters"),
0141                                                       0,
0142                                                       false,
0143                                                       false,
0144                                                       ibook));
0145         }
0146       }
0147     } else {
0148       iInfoTag++;
0149       // tag info retrievel is deferred (needs availability of EventSetup)
0150       const InputTag& moduleLabel = iModule->getParameter<InputTag>("label");
0151       const string& folderName = iModule->getParameter<string>("folder");
0152
0153       // eta loop
0154       for (int iEta = iEtaStart; iEta < iEtaEnd; ++iEta) {
0155         // pt loop
0156         for (int iPt = iPtStart; iPt < iPtEnd; ++iPt) {
0157           const EtaPtBin& etaPtBin = getEtaPtBin(iEta, iPt);
0158
0159           // Instantiate the tagInfo plotter
0160           binTagInfoPlotters.at(iInfoTag).push_back(
0161               theFactory.buildPlotter(dataFormatType,
0162                                       moduleLabel.label(),
0163                                       etaPtBin,
0164                                       iModule->getParameter<edm::ParameterSet>("parameters"),
0165                                       folderName,
0166                                       0,
0167                                       false,
0168                                       ibook));
0169         }
0170       }
0171     }
0172   }
0173 }
0174
0175 EtaPtBin BTagPerformanceAnalyzerOnData::getEtaPtBin(const int& iEta, const int& iPt) {
0176   // DEFINE BTagBin:
0177   bool etaActive_, ptActive_;
0178   double etaMin_, etaMax_, ptMin_, ptMax_;
0179
0180   if (iEta != -1) {
0181     etaActive_ = true;
0182     etaMin_ = etaRanges[iEta];
0183     etaMax_ = etaRanges[iEta + 1];
0184   } else {
0185     etaActive_ = false;
0186     etaMin_ = etaRanges[0];
0187     etaMax_ = etaRanges[etaRanges.size() - 1];
0188   }
0189
0190   if (iPt != -1) {
0191     ptActive_ = true;
0192     ptMin_ = ptRanges[iPt];
0193     ptMax_ = ptRanges[iPt + 1];
0194   } else {
0195     ptActive_ = false;
0196     ptMin_ = ptRanges[0];
0197     ptMax_ = ptRanges[ptRanges.size() - 1];
0198   }
0199   return EtaPtBin(etaActive_, etaMin_, etaMax_, ptActive_, ptMin_, ptMax_);
0200 }
0201
0202 BTagPerformanceAnalyzerOnData::~BTagPerformanceAnalyzerOnData() {
0203   /*for (unsigned int iJetLabel = 0; iJetLabel != binJetTagPlotters.size(); ++iJetLabel) {
0204     int plotterSize =  binJetTagPlotters[iJetLabel].size();
0205     for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter) {
0206       delete binJetTagPlotters[iJetLabel][iPlotter];
0207     }
0208   }
0209
0210   for (vector<vector<TagCorrelationPlotter*> >::iterator iJetLabel = binTagCorrelationPlotters.begin();
0211        iJetLabel != binTagCorrelationPlotters.end(); ++iJetLabel)
0212     for(vector<TagCorrelationPlotter*>::iterator iPlotter = iJetLabel->begin(); iPlotter != iJetLabel->end(); ++iPlotter)
0213       delete *iPlotter;
0214
0215   for (unsigned int iJetLabel = 0; iJetLabel != binTagInfoPlotters.size(); ++iJetLabel) {
0216     int plotterSize =  binTagInfoPlotters[iJetLabel].size();
0217     for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter) {
0218       delete binTagInfoPlotters[iJetLabel][iPlotter];
0219     }
0220   }*/
0221 }
0222
0223 void BTagPerformanceAnalyzerOnData::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0224   edm::Handle<reco::SoftLeptonTagInfoCollection> infoHandle;
0225   iEvent.getByToken(slInfoToken, infoHandle);
0226
0227   //Get JEC
0228   const JetCorrector* corrector = nullptr;
0229   if (doJEC) {
0230     edm::Handle<GenEventInfoProduct> genInfoHandle;  //check if data or MC
0231     iEvent.getByToken(genToken, genInfoHandle);
0232     edm::Handle<JetCorrector> corrHandle;
0233     if (!genInfoHandle.isValid())
0234       iEvent.getByToken(jecDataToken, corrHandle);
0235     else
0236       iEvent.getByToken(jecMCToken, corrHandle);
0237     corrector = corrHandle.product();
0238   }
0239
0240   // Look first at the jetTags
0241
0242   for (unsigned int iJetLabel = 0; iJetLabel != jetTagInputTags.size(); ++iJetLabel) {
0243     edm::Handle<reco::JetTagCollection> tagHandle;
0244     iEvent.getByToken(jetTagToken[iJetLabel], tagHandle);
0245     //
0246     // insert check on the presence of the collections
0247     //
0248
0249     if (!tagHandle.isValid()) {
0250       edm::LogWarning("BTagPerformanceAnalyzerOnData")
0251           << " Collection " << jetTagInputTags[iJetLabel] << " not present. Skipping it for this event.";
0252       continue;
0253     }
0254
0255     const reco::JetTagCollection& tagColl = *(tagHandle.product());
0256     LogDebug("Info") << "Found " << tagColl.size() << " B candidates in collection " << jetTagInputTags[iJetLabel];
0257
0258     int plotterSize = binJetTagPlotters[iJetLabel].size();
0259     for (JetTagCollection::const_iterator tagI = tagColl.begin(); tagI != tagColl.end(); ++tagI) {
0260       //JEC
0261       double jec = 1.0;
0262       if (doJEC && corrector) {
0263         jec = corrector->correction(*(tagI->first));
0264       }
0265
0266       if (!jetSelector(*(tagI->first), -1, infoHandle, jec))
0267         continue;
0268
0269       for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter) {
0270         bool inBin = binJetTagPlotters[iJetLabel][iPlotter]->etaPtBin().inBin(*tagI->first, jec);
0271         // Fill histograms if in desired pt/rapidity bin.
0272         if (inBin)
0273           binJetTagPlotters[iJetLabel][iPlotter]->analyzeTag(*tagI, jec, -1);
0274       }
0275     }
0276     for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter) {
0277       binJetTagPlotters[iJetLabel][iPlotter]->analyzeTag();
0278     }
0279   }
0280
0281   // Now look at Tag Correlations
0282   for (unsigned int iJetLabel = 0; iJetLabel != tagCorrelationInputTags.size(); ++iJetLabel) {
0283     const std::pair<edm::EDGetTokenT<reco::JetTagCollection>, edm::EDGetTokenT<reco::JetTagCollection>>& inputTokens =
0284         tagCorrelationToken[iJetLabel];
0285     edm::Handle<reco::JetTagCollection> tagHandle1;
0286     iEvent.getByToken(inputTokens.first, tagHandle1);
0287     const reco::JetTagCollection& tagColl1 = *(tagHandle1.product());
0288
0289     edm::Handle<reco::JetTagCollection> tagHandle2;
0290     iEvent.getByToken(inputTokens.second, tagHandle2);
0291     const reco::JetTagCollection& tagColl2 = *(tagHandle2.product());
0292
0293     int plotterSize = binTagCorrelationPlotters[iJetLabel].size();
0294     for (JetTagCollection::const_iterator tagI = tagColl1.begin(); tagI != tagColl1.end(); ++tagI) {
0295       //JEC
0296       double jec = 1.0;
0297       if (doJEC && corrector) {
0298         jec = corrector->correction(*(tagI->first));
0299       }
0300
0301       if (!jetSelector(*(tagI->first), -1, infoHandle, jec))
0302         continue;
0303
0304       for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter) {
0305         bool inBin = binTagCorrelationPlotters[iJetLabel][iPlotter]->etaPtBin().inBin(*(tagI->first), jec);
0306
0307         if (inBin) {
0308           double discr2 = tagColl2[tagI->first];
0309           binTagCorrelationPlotters[iJetLabel][iPlotter]->analyzeTags(tagI->second, discr2, -1);
0310         }
0311       }
0312     }
0313   }
0314
0315   // Now look at the TagInfos
0316   for (unsigned int iJetLabel = 0; iJetLabel != tiDataFormatType.size(); ++iJetLabel) {
0317     int plotterSize = binTagInfoPlotters[iJetLabel].size();
0318     for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter)
0319       binTagInfoPlotters[iJetLabel][iPlotter]->setEventSetup(iSetup);
0320
0321     vector<edm::EDGetTokenT<edm::View<reco::BaseTagInfo>>>& tokens = tagInfoToken[iJetLabel];
0322     //check number of tag infos = expected number of tag infos
0323     vector<string> labels = binTagInfoPlotters[iJetLabel][0]->tagInfoRequirements();
0324     if (labels.empty())
0325       labels.push_back("label");
0326     if (labels.size() != tokens.size())
0327       throw cms::Exception("Configuration")
0328           << "Different number of Tag Infos than expected" << labels.size() << tokens.size() << endl;
0329
0330     unsigned int nInputTags = tokens.size();
0331     vector<edm::Handle<View<BaseTagInfo>>> tagInfoHandles(nInputTags);
0332     edm::ProductID jetProductID;
0333     unsigned int nTagInfos = 0;
0334     for (unsigned int iInputTags = 0; iInputTags < tokens.size(); ++iInputTags) {
0335       edm::Handle<View<BaseTagInfo>>& tagInfoHandle = tagInfoHandles[iInputTags];
0336       iEvent.getByToken(tokens[iInputTags], tagInfoHandle);
0337       //
0338       // protect against missing products
0339       //
0340       if (tagInfoHandle.isValid() == false) {
0341         edm::LogWarning("BTagPerformanceAnalyzerOnData")
0342             << " Collection " << tagInfoInputTags[iJetLabel][iInputTags] << " not present. Skipping it for this event.";
0343         continue;
0344       }
0345
0346       unsigned int size = tagInfoHandle->size();
0347       LogDebug("Info") << "Found " << size << " B candidates in collection " << tagInfoInputTags[iJetLabel][iInputTags];
0348       edm::ProductID thisProductID = (size > 0) ? (*tagInfoHandle)[0].jet().id() : edm::ProductID();
0349       if (iInputTags == 0) {
0350         jetProductID = thisProductID;
0351         nTagInfos = size;
0352       } else if (jetProductID != thisProductID)
0353         throw cms::Exception("Configuration") << "TagInfos are referencing a different jet collection." << endl;
0354       else if (nTagInfos != size)
0355         throw cms::Exception("Configuration") << "TagInfo collections are having a different size." << endl;
0356     }
0357
0358     for (unsigned int iTagInfos = 0; iTagInfos < nTagInfos; ++iTagInfos) {
0359       vector<const BaseTagInfo*> baseTagInfos(nInputTags);
0360       edm::RefToBase<Jet> jetRef;
0361       for (unsigned int iTagInfo = 0; iTagInfo < nInputTags; iTagInfo++) {
0362         const BaseTagInfo& baseTagInfo = (*tagInfoHandles[iTagInfo])[iTagInfos];
0363         if (iTagInfo == 0)
0364           jetRef = baseTagInfo.jet();
0365         else if (baseTagInfo.jet() != jetRef)
0366           throw cms::Exception("Configuration") << "TagInfos pointing to different jets." << endl;
0367         baseTagInfos[iTagInfo] = &baseTagInfo;
0368       }
0369
0370       //JEC
0371       double jec = 1.0;
0372       if (doJEC && corrector) {
0373         jec = corrector->correction(*(jetRef));
0374       }
0375
0376       if (!jetSelector(*jetRef, -1, infoHandle, jec))
0377         continue;
0378
0379       for (int iPlotter = 0; iPlotter != plotterSize; ++iPlotter) {
0380         bool inBin = binTagInfoPlotters[iJetLabel][iPlotter]->etaPtBin().inBin(*jetRef, jec);
0381         // Fill histograms if in desired pt/rapidity bin.
0382         if (inBin)
0383           binTagInfoPlotters[iJetLabel][iPlotter]->analyzeTag(baseTagInfos, jec, -1);
0384       }
0385     }
0386   }
0387 }
0388
0389 //define this as a plug-in
0390 DEFINE_FWK_MODULE(BTagPerformanceAnalyzerOnData);