Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-04-06 12:09:48

0001 #include "DQMOffline/RecoB/interface/JetTagPlotter.h"
0002 #include "DQMOffline/RecoB/interface/Tools.h"
0003 #include "FWCore/Utilities/interface/isFinite.h"
0004 #include "DQMServices/Core/interface/DQMStore.h"
0005 
0006 #include <iostream>
0007 
0008 using namespace std;
0009 using namespace RecoBTag;
0010 
0011 JetTagPlotter::JetTagPlotter(const std::string& tagName,
0012                              const EtaPtBin& etaPtBin,
0013                              const edm::ParameterSet& pSet,
0014                              unsigned int mc,
0015                              bool wf,
0016                              DQMStore::IBooker& ibook,
0017                              bool doCTagPlots /*=false*/,
0018                              bool doDifferentialPlots /*=false*/,
0019                              double discrCut /*=-999.*/)
0020     : BaseBTagPlotter(tagName, etaPtBin),
0021       discrStart_(pSet.getParameter<double>("discriminatorStart")),
0022       discrEnd_(pSet.getParameter<double>("discriminatorEnd")),
0023       nBinEffPur_(pSet.getParameter<int>("nBinEffPur")),
0024       startEffPur_(pSet.getParameter<double>("startEffPur")),
0025       endEffPur_(pSet.getParameter<double>("endEffPur")),
0026       mcPlots_(mc),
0027       willFinalize_(wf),
0028       doCTagPlots_(doCTagPlots),
0029       doDifferentialPlots_(doDifferentialPlots),
0030       cutValue_(discrCut) {
0031   // to have a shorter name .....
0032   const std::string& es = theExtensionString;
0033   const std::string jetTagDir(es.substr(1));
0034 
0035   if (willFinalize_)
0036     return;
0037 
0038   //added to count the number of jets by event : 0=DATA or NI, 1to5=quarks u,d,s,c,b , 6=gluon
0039   int nFl = 1;
0040   if (mcPlots_)
0041     nFl = 8;
0042   nJets_.resize(nFl, 0);
0043 
0044   if (mcPlots_) {
0045     // jet flavour
0046     dJetFlav_ = std::make_unique<FlavourHistograms<int>>(
0047         "jetFlavour" + es, "Jet Flavour", 22, -0.5, 21.5, false, false, false, "b", jetTagDir, mcPlots_, ibook);
0048   }
0049 
0050   // jet multiplicity
0051   jetMultiplicity_ = std::make_unique<FlavourHistograms<int>>(
0052       "jetMultiplicity" + es, "Jet Multiplicity", 11, -0.5, 10.5, false, true, true, "b", jetTagDir, mcPlots_, ibook);
0053 
0054   // Discriminator: again with reasonable binning
0055   dDiscriminator_ = std::make_unique<FlavourHistograms<double>>(
0056       "discr" + es, "Discriminator", 102, discrStart_, discrEnd_, false, true, true, "b", jetTagDir, mcPlots_, ibook);
0057   dDiscriminator_->settitle("Discriminant");
0058   // reconstructed jet momentum
0059   dJetRecMomentum_ = std::make_unique<FlavourHistograms<double>>(
0060       "jetMomentum" + es, "jet momentum", 350, 0.0, 350.0, false, false, true, "b", jetTagDir, mcPlots_, ibook);
0061 
0062   // reconstructed jet transverse momentum
0063   dJetRecPt_ = std::make_unique<FlavourHistograms<double>>(
0064       "jetPt" + es, "jet pt", 350, 0.0, 350.0, false, false, true, "b", jetTagDir, mcPlots_, ibook);
0065 
0066   // reconstructed jet eta
0067   dJetRecPseudoRapidity_ = std::make_unique<FlavourHistograms<double>>("jetEta" + es,
0068                                                                        "jet eta",
0069                                                                        20,
0070                                                                        -etaPtBin.getEtaMax(),
0071                                                                        etaPtBin.getEtaMax(),
0072                                                                        false,
0073                                                                        false,
0074                                                                        true,
0075                                                                        "b",
0076                                                                        jetTagDir,
0077                                                                        mcPlots_,
0078                                                                        ibook);
0079 
0080   // reconstructed jet phi
0081   dJetRecPhi_ = std::make_unique<FlavourHistograms<double>>(
0082       "jetPhi" + es, "jet phi", 20, -M_PI, M_PI, false, false, true, "b", jetTagDir, mcPlots_, ibook);
0083 
0084   if (doDifferentialPlots_) {
0085     // jet Phi larger than requested discrimnator cut
0086     dJetPhiDiscrCut_ = std::make_unique<FlavourHistograms<double>>("jetPhi_diffEff" + es,
0087                                                                    "Efficiency vs. jet Phi for discriminator above cut",
0088                                                                    20,
0089                                                                    -M_PI,
0090                                                                    M_PI,
0091                                                                    false,
0092                                                                    false,
0093                                                                    true,
0094                                                                    "b",
0095                                                                    jetTagDir,
0096                                                                    mcPlots_,
0097                                                                    ibook);
0098 
0099     // jet Eta larger than requested discrimnator cut
0100     dJetPseudoRapidityDiscrCut_ =
0101         std::make_unique<FlavourHistograms<double>>("jetEta_diffEff" + es,
0102                                                     "Efficiency vs. jet eta for discriminator above cut",
0103                                                     20,
0104                                                     -etaPtBin.getEtaMax(),
0105                                                     etaPtBin.getEtaMax(),
0106                                                     false,
0107                                                     false,
0108                                                     true,
0109                                                     "b",
0110                                                     jetTagDir,
0111                                                     mcPlots_,
0112                                                     ibook);
0113   }
0114 }
0115 
0116 JetTagPlotter::~JetTagPlotter() {}
0117 
0118 void JetTagPlotter::epsPlot(const std::string& name) {
0119   if (!willFinalize_) {
0120     dJetFlav_->epsPlot(name);
0121     jetMultiplicity_->epsPlot(name);
0122     dDiscriminator_->epsPlot(name);
0123     dJetRecMomentum_->epsPlot(name);
0124     dJetRecPt_->epsPlot(name);
0125     dJetRecPseudoRapidity_->epsPlot(name);
0126     dJetRecPhi_->epsPlot(name);
0127   } else {
0128     effPurFromHistos_->epsPlot(name);
0129   }
0130 }
0131 
0132 void JetTagPlotter::psPlot(const std::string& name) {
0133   std::string cName = "JetTagPlots" + theExtensionString;
0134   setTDRStyle()->cd();
0135   TCanvas canvas(cName.c_str(), cName.c_str(), 600, 900);
0136   canvas.UseCurrentStyle();
0137 
0138   canvas.Divide(2, 3);
0139   canvas.Print((name + cName + ".ps[").c_str());
0140   if (!willFinalize_) {
0141     canvas.cd(1);
0142     dJetFlav_->plot();
0143     canvas.cd(2);
0144     canvas.cd(3);
0145     dDiscriminator_->plot();
0146     canvas.cd(4);
0147     dJetRecMomentum_->plot();
0148     canvas.cd(5);
0149     dJetRecPt_->plot();
0150     canvas.cd(6);
0151     dJetRecPseudoRapidity_->plot();
0152     canvas.Print((name + cName + ".ps").c_str());
0153     canvas.Clear();
0154     canvas.Divide(2, 3);
0155 
0156     jetMultiplicity_->plot();
0157     canvas.Print((name + cName + ".ps").c_str());
0158     canvas.Clear();
0159 
0160     canvas.cd(1);
0161     dJetRecPhi_->plot();
0162     canvas.cd(2);
0163     canvas.cd(3);
0164     canvas.cd(4);
0165   } else {
0166     canvas.cd(5);
0167     effPurFromHistos_->discriminatorNoCutEffic().plot();
0168     canvas.cd(6);
0169     effPurFromHistos_->discriminatorCutEfficScan().plot();
0170     canvas.Print((name + cName + ".ps").c_str());
0171     canvas.Clear();
0172     canvas.Divide(2, 3);
0173     canvas.cd(1);
0174     effPurFromHistos_->plot();
0175   }
0176   canvas.Print((name + cName + ".ps").c_str());
0177   canvas.Print((name + cName + ".ps]").c_str());
0178 }
0179 
0180 void JetTagPlotter::analyzeTag()  //here jetFlavour not needed
0181 {
0182   //to use on data
0183   jetMultiplicity_->fill(-1, nJets_[0]);
0184   nJets_[0] = 0;  //reset to 0 before the next event
0185 }
0186 
0187 void JetTagPlotter::analyzeTag(float w) {
0188   if (mcPlots_) {
0189     //to use with MC
0190     int totNJets = 0;
0191     int udsNJets = 0;
0192     int udsgNJets = 0;
0193     for (int i = 0; i < 8; i++) {
0194       totNJets += nJets_[i];
0195       if (i > 0 && i < 4)
0196         udsNJets += nJets_[i];
0197       if ((i > 0 && i < 4) || i == 6)
0198         udsgNJets += nJets_[i];
0199       if (i <= 5 && i >= 1)
0200         jetMultiplicity_->fill(i, nJets_[i], w);
0201       else if (i == 6)
0202         jetMultiplicity_->fill(21, nJets_[i], w);
0203       else if (i == 7)
0204         jetMultiplicity_->fill(20, nJets_[i], w);
0205       else
0206         jetMultiplicity_->fill(0, nJets_[i], w);
0207       nJets_[i] = 0;  //reset to 0 before the next event
0208     }
0209     jetMultiplicity_->fill(-1, totNJets, w);  //total number of jets in the event
0210     jetMultiplicity_->fill(123, udsNJets, w);
0211     jetMultiplicity_->fill(12321, udsgNJets, w);
0212   } else {
0213     int totNJets = 0;
0214     for (int i = 0; i < 8; i++) {
0215       totNJets += nJets_[i];
0216       nJets_[i] = 0;
0217     }
0218     jetMultiplicity_->fill(-1, totNJets, w);
0219   }
0220 }
0221 
0222 void JetTagPlotter::analyzeTag(const reco::Jet& jet, double jec, float discriminator, int jetFlavour, float w /*=1*/) {
0223   if (mcPlots_) {
0224     dJetFlav_->fill(jetFlavour, jetFlavour, w);
0225     if (abs(jetFlavour) > 0 && abs(jetFlavour) < 6)
0226       nJets_[abs(jetFlavour)] += 1;  //quarks 1 to 5
0227     else if (abs(jetFlavour) == 21)
0228       nJets_[6] += 1;  //gluons
0229     else if (jetFlavour == 20)
0230       nJets_[7] += 1;  //PU
0231     else
0232       nJets_[0] += 1;  //NI
0233   } else {
0234     nJets_[0] += 1;
0235   }
0236   if (edm::isNotFinite(discriminator))
0237     dDiscriminator_->fill(jetFlavour, -999.0, w);
0238   else
0239     dDiscriminator_->fill(jetFlavour, discriminator, w);
0240   dJetRecMomentum_->fill(jetFlavour, jet.p() * jec, w);
0241   dJetRecPt_->fill(jetFlavour, jet.pt() * jec, w);
0242   dJetRecPseudoRapidity_->fill(jetFlavour, jet.eta(), w);
0243   dJetRecPhi_->fill(jetFlavour, jet.phi(), w);
0244   if (doDifferentialPlots_) {
0245     if (edm::isFinite(discriminator) && discriminator > cutValue_) {
0246       dJetPhiDiscrCut_->fill(jetFlavour, jet.phi(), w);
0247       dJetPseudoRapidityDiscrCut_->fill(jetFlavour, jet.eta(), w);
0248     }
0249   }
0250 }
0251 
0252 void JetTagPlotter::analyzeTag(const reco::JetTag& jetTag, double jec, int jetFlavour, float w /*=1*/) {
0253   if (mcPlots_) {
0254     dJetFlav_->fill(jetFlavour, jetFlavour, w);
0255     if (abs(jetFlavour) > 0 && abs(jetFlavour) < 6)
0256       nJets_[abs(jetFlavour)] += 1;  //quarks 1 to 5
0257     else if (abs(jetFlavour) == 21)
0258       nJets_[6] += 1;  //gluons
0259     else if (jetFlavour == 20)
0260       nJets_[7] += 1;  //PU
0261     else
0262       nJets_[0] += 1;  //NI
0263   } else {
0264     nJets_[0] += 1;
0265   }
0266   const auto& discriminator = jetTag.second;
0267   if (edm::isNotFinite(discriminator))
0268     dDiscriminator_->fill(jetFlavour, -999.0, w);
0269   else
0270     dDiscriminator_->fill(jetFlavour, discriminator, w);
0271   dJetRecMomentum_->fill(jetFlavour, jetTag.first->p() * jec, w);
0272   dJetRecPt_->fill(jetFlavour, jetTag.first->pt() * jec, w);
0273   dJetRecPseudoRapidity_->fill(jetFlavour, jetTag.first->eta(), w);
0274   dJetRecPhi_->fill(jetFlavour, jetTag.first->phi(), w);
0275   if (doDifferentialPlots_) {
0276     if (edm::isFinite(discriminator) && discriminator > cutValue_) {
0277       dJetPhiDiscrCut_->fill(jetFlavour, jetTag.first->phi(), w);
0278       dJetPseudoRapidityDiscrCut_->fill(jetFlavour, jetTag.first->eta(), w);
0279     }
0280   }
0281 }
0282 
0283 void JetTagPlotter::finalize(DQMStore::IBooker& ibook_, DQMStore::IGetter& igetter_) {
0284   //
0285   // final processing:
0286   // produce the misid. vs. eff histograms
0287   //
0288   const std::string& es = theExtensionString;
0289   const std::string jetTagDir(es.substr(1));
0290   dDiscriminator_ = std::make_unique<FlavourHistograms<double>>(
0291       "discr" + es, "Discriminator", 102, discrStart_, discrEnd_, "b", jetTagDir, mcPlots_, igetter_);
0292 
0293   effPurFromHistos_ = std::make_unique<EffPurFromHistos>(
0294       *dDiscriminator_, jetTagDir, mcPlots_, ibook_, nBinEffPur_, startEffPur_, endEffPur_);
0295   effPurFromHistos_->doCTagPlots(doCTagPlots_);
0296   effPurFromHistos_->compute(ibook_);
0297 
0298   // Produce the differentiel efficiency vs. kinematical variables
0299   if (doDifferentialPlots_) {
0300     dJetRecPhi_ = std::make_unique<FlavourHistograms<double>>(
0301         "jetPhi" + es, "jet phi", 20, -M_PI, M_PI, "b", jetTagDir, mcPlots_, igetter_);
0302     dJetPhiDiscrCut_ = std::make_unique<FlavourHistograms<double>>("jetPhi_diffEff" + es,
0303                                                                    "Efficiency vs. jet Phi for discriminator above cut",
0304                                                                    20,
0305                                                                    -M_PI,
0306                                                                    M_PI,
0307                                                                    "b",
0308                                                                    jetTagDir,
0309                                                                    mcPlots_,
0310                                                                    igetter_);
0311     dJetPhiDiscrCut_->divide(*dJetRecPhi_);
0312     dJetPhiDiscrCut_->setEfficiencyFlag();
0313 
0314     dJetRecPseudoRapidity_ = std::make_unique<FlavourHistograms<double>>(
0315         "jetEta" + es, "jet eta", 20, -etaPtBin_.getEtaMax(), etaPtBin_.getEtaMax(), "b", jetTagDir, mcPlots_, igetter_);
0316     dJetPseudoRapidityDiscrCut_ =
0317         std::make_unique<FlavourHistograms<double>>("jetEta_diffEff" + es,
0318                                                     "Efficiency vs. jet eta for discriminator above cut",
0319                                                     20,
0320                                                     -etaPtBin_.getEtaMax(),
0321                                                     etaPtBin_.getEtaMax(),
0322                                                     "b",
0323                                                     jetTagDir,
0324                                                     mcPlots_,
0325                                                     igetter_);
0326     dJetPseudoRapidityDiscrCut_->divide(*dJetRecPseudoRapidity_);
0327     dJetPseudoRapidityDiscrCut_->setEfficiencyFlag();
0328   }
0329 }