Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2023-03-17 11:16:31

0001 // -*- C++ -*-
0002 //
0003 // Package:    PatBTag
0004 // Class:      PatBTagCommonHistos
0005 //
0006 /**\class PatBTagCommonHistos PatBTagCommonHistos.cc
0007 
0008  Description: <Define and Fill common set of histograms depending on flavor>
0009 
0010  Implementation:
0011  
0012  Create a container of histograms. 
0013 */
0014 //
0015 // Original Author:  J. E. Ramirez
0016 //
0017 //
0018 
0019 // system include files
0020 #include <memory>
0021 
0022 // user include files
0023 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0024 #include "FWCore/Utilities/interface/InputTag.h"
0025 #include "DataFormats/PatCandidates/interface/Jet.h"
0026 #include "PhysicsTools/PatExamples/interface/PatBTagCommonHistos.h"
0027 
0028 #include "FWCore/ServiceRegistry/interface/Service.h"
0029 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0030 
0031 //
0032 // constructors and destructor
0033 //
0034 PatBTagCommonHistos::PatBTagCommonHistos(const edm::ParameterSet& iConfig)
0035     : histocontainer_(),
0036       BTagger(iConfig.getParameter<edm::ParameterSet>("BJetOperatingPoints"))
0037 //now do what ever initialization is needed
0038 {
0039   edm::ParameterSet PatBjet_(iConfig.getParameter<edm::ParameterSet>("BjetTag"));
0040   BTagdisccut_ = PatBjet_.getUntrackedParameter<double>("mindiscriminatorcut", 5.0);
0041   BTagdiscriminator_ = PatBjet_.getParameter<std::string>("discriminator");
0042   BTagpurity_ = PatBjet_.getParameter<std::string>("purity");
0043 }
0044 
0045 PatBTagCommonHistos::~PatBTagCommonHistos() {
0046   // do anything here that needs to be done at desctruction time
0047   // (e.g. close files, deallocate resources etc.)
0048 }
0049 
0050 //
0051 // member functions
0052 //
0053 
0054 // ------------ method called to for each event  ------------
0055 void PatBTagCommonHistos::Fill(edm::View<pat::Jet>::const_iterator& jet_iter, std::string flavor) {
0056   float isb = jet_iter->bDiscriminator(BTagdiscriminator_);
0057   //
0058   // no pt cut histograms
0059   //
0060   if (fabs(jet_iter->eta()) < 2.4 && fabs(jet_iter->eta()) > 0) {
0061     //
0062     //Fill histos for Loose(defined in configuration file) cut
0063     //tagged using auxiliar function (Loose)
0064     //
0065     if (BTagger.IsbTag(*jet_iter, BTagpurity_, BTagdiscriminator_)) {
0066       histocontainer_["jet_pt_uncorr_" + flavor + "_tagged"]->Fill(jet_iter->correctedJet("raw").pt());
0067       histocontainer_["jet_pt_" + flavor + "_tagged"]->Fill(jet_iter->pt());
0068       histocontainer_["jet_eta_" + flavor + "_tagged"]->Fill(jet_iter->eta());
0069     }
0070     //
0071     // Fill histos
0072     //    tagged minimum discriminant cut criteria
0073     //
0074     if (isb > float(BTagdisccut_)) {
0075       histocontainer_["jet_pt_uncorr_" + flavor + "_taggedmincut"]->Fill(jet_iter->correctedJet("raw").pt());
0076       histocontainer_["jet_pt_" + flavor + "_taggedmincut"]->Fill(jet_iter->pt());
0077       histocontainer_["jet_eta_" + flavor + "_taggedmincut"]->Fill(jet_iter->eta());
0078     }
0079     //
0080     //fill taggable jets histograms (tracks in jet > 0,1,2)
0081     std::map<int, std::string> tagbl;
0082     tagbl[0] = "0";
0083     tagbl[1] = "1";
0084     tagbl[2] = "2";
0085     for (size_t i = 0; i < tagbl.size(); i++)
0086       if (jet_iter->associatedTracks().size() > i) {
0087         histocontainer_["jet_pt_uncorr_" + flavor + "_taggable" + tagbl[i]]->Fill(jet_iter->correctedJet("raw").pt());
0088         histocontainer_["jet_pt_" + flavor + "_taggable" + tagbl[i]]->Fill(jet_iter->pt());
0089         histocontainer_["jet_eta_" + flavor + "_taggable" + tagbl[i]]->Fill(jet_iter->eta());
0090         if (jet_iter->pt() < 30) {
0091           histocontainer_["jet_eta_" + flavor + "_taggable" + tagbl[i] + "_030"]->Fill(jet_iter->eta());
0092         }
0093         if (jet_iter->pt() > 30 && jet_iter->pt() < 50) {
0094           histocontainer_["jet_eta_" + flavor + "_taggable" + tagbl[i] + "_3050"]->Fill(jet_iter->eta());
0095         }
0096         if (jet_iter->pt() > 50) {
0097           histocontainer_["jet_eta_" + flavor + "_taggable" + tagbl[i] + "_50"]->Fill(jet_iter->eta());
0098         }
0099       }
0100     //
0101     // Fill histos needed for normalization
0102     // uncorrected pt distributions
0103     // corrected pt distributions
0104     // eta distributions
0105     // scatter plots for pts
0106     histocontainer_["jet_pt_uncorr_" + flavor]->Fill(jet_iter->correctedJet("raw").pt());
0107     histocontainer_["jet_pt_" + flavor]->Fill(jet_iter->pt());
0108     histocontainer_["jet_eta_" + flavor]->Fill(jet_iter->eta());
0109     histocontainer_["tracks_in_jet_" + flavor]->Fill(jet_iter->associatedTracks().size());
0110     h2_["jet_scatter_pt_" + flavor]->Fill(jet_iter->pt(), jet_iter->correctedJet("raw").pt());
0111     for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0112       histocontainer_["pt_tracks_in_jet_" + flavor]->Fill(jet_iter->associatedTracks()[itrack]->pt());
0113     }
0114     if (jet_iter->pt() < 30) {
0115       histocontainer_["jet_eta_" + flavor + "030"]->Fill(jet_iter->eta());
0116       histocontainer_["tracks_in_jet_" + flavor + "030"]->Fill(jet_iter->associatedTracks().size());
0117       for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0118         histocontainer_["pt_tracks_in_jet_" + flavor + "030"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
0119       }
0120       if (fabs(jet_iter->eta()) < 1.5) {
0121         for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0122           histocontainer_["pt_tracks_in_jet_" + flavor + "_center030"]->Fill(
0123               jet_iter->associatedTracks()[itrack]->pt());
0124         }
0125       }
0126       if (fabs(jet_iter->eta()) > 1.5 && fabs(jet_iter->eta()) < 2.4) {
0127         for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0128           histocontainer_["pt_tracks_in_jet_" + flavor + "_sides030"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
0129         }
0130       }
0131     }
0132     if (jet_iter->pt() > 30 && jet_iter->pt() < 50) {
0133       histocontainer_["jet_eta_" + flavor + "3050"]->Fill(jet_iter->eta());
0134       histocontainer_["tracks_in_jet_" + flavor + "3050"]->Fill(jet_iter->associatedTracks().size());
0135       for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0136         histocontainer_["pt_tracks_in_jet_" + flavor + "3050"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
0137       }
0138       if (fabs(jet_iter->eta()) < 1.5) {
0139         for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0140           histocontainer_["pt_tracks_in_jet_" + flavor + "_center3050"]->Fill(
0141               jet_iter->associatedTracks()[itrack]->pt());
0142         }
0143       }
0144       if (fabs(jet_iter->eta()) > 1.5 && fabs(jet_iter->eta()) < 2.4) {
0145         for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0146           histocontainer_["pt_tracks_in_jet_" + flavor + "_sides3050"]->Fill(
0147               jet_iter->associatedTracks()[itrack]->pt());
0148         }
0149       }
0150     }
0151     if (jet_iter->pt() > 50) {
0152       histocontainer_["jet_eta_" + flavor + "50"]->Fill(jet_iter->eta());
0153       histocontainer_["tracks_in_jet_" + flavor + "50"]->Fill(jet_iter->associatedTracks().size());
0154       for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0155         histocontainer_["pt_tracks_in_jet_" + flavor + "50"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
0156       }
0157       if (fabs(jet_iter->eta()) < 1.5) {
0158         for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0159           histocontainer_["pt_tracks_in_jet_" + flavor + "_center50"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
0160         }
0161       }
0162       if (fabs(jet_iter->eta()) > 1.5 && fabs(jet_iter->eta()) < 2.4) {
0163         for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0164           histocontainer_["pt_tracks_in_jet_" + flavor + "_sides50"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
0165         }
0166       }
0167     }
0168     if (fabs(jet_iter->eta()) < 1.5) {
0169       for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0170         histocontainer_["pt_tracks_in_jet_" + flavor + "_center"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
0171       }
0172     }
0173     if (fabs(jet_iter->eta()) > 1.5 && fabs(jet_iter->eta()) < 2.4) {
0174       for (size_t itrack = 0; itrack < jet_iter->associatedTracks().size(); ++itrack) {
0175         histocontainer_["pt_tracks_in_jet_" + flavor + "_sides"]->Fill(jet_iter->associatedTracks()[itrack]->pt());
0176       }
0177     }
0178 
0179   }  //endif
0180 }  //end function
0181 // ------------ method called once each job just before starting event loop  ------------
0182 // ------------  This function is needed to set a group of histogram  -------------------
0183 void PatBTagCommonHistos::Set(std::string flavor) {
0184   const int ntptarray = 23;
0185   const int njptarray = 14;
0186   const int netaarray = 19;
0187   Double_t jetetabins[netaarray] = {
0188       -2.5, -2.0, -1.75, -1.5, -1.25, -1.0, -0.75, -0.5, -0.25, 0.0, 0.25, 0.5, 0.75, 1.0, 1.25, 1.5, 1.75, 2.0, 2.5};
0189   Double_t jetptxbins[njptarray] = {0., 10., 20., 30., 40., 50., 60., 70., 80, 90., 100., 120., 140., 230.};
0190   Double_t jetpttbins[ntptarray] = {0.,  1.,  2.,  4.,  6.,  8., 10., 15.,  20.,  25.,  30., 35.,
0191                                     40., 45., 50., 60., 70., 80, 90., 100., 120., 140., 230.};
0192   edm::Service<TFileService> fs;
0193   std::string histoid;
0194   std::string histotitle;
0195 
0196   //Define histograms
0197   histoid = "jet_pt_" + flavor;
0198   histotitle = flavor + " jet p_{T} [GeV/c]";
0199   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), njptarray - 1, jetptxbins);
0200   histoid = "jet_pt_uncorr_" + flavor;
0201   histotitle = flavor + " jet uncorrected p_{T} [GeV/c]";
0202   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), njptarray - 1, jetptxbins);
0203   histoid = "jet_eta_" + flavor;
0204   histocontainer_["jet_eta_" + flavor] = fs->make<TH1D>(histoid.c_str(), "jet #eta", netaarray - 1, jetetabins);
0205   histoid = "jet_eta_" + flavor + "030";
0206   histocontainer_["jet_eta_" + flavor + "030"] = fs->make<TH1D>(histoid.c_str(), "jet #eta", netaarray - 1, jetetabins);
0207   histoid = "jet_eta_" + flavor + "3050";
0208   histocontainer_["jet_eta_" + flavor + "3050"] =
0209       fs->make<TH1D>(histoid.c_str(), "jet #eta", netaarray - 1, jetetabins);
0210   histoid = "jet_eta_" + flavor + "50";
0211   histocontainer_["jet_eta_" + flavor + "50"] = fs->make<TH1D>(histoid.c_str(), "jet #eta", netaarray - 1, jetetabins);
0212   histoid = "jet_pt_" + flavor + "_tagged";
0213   histocontainer_["jet_pt_" + flavor + "_tagged"] =
0214       fs->make<TH1D>(histoid.c_str(), "jet_tagged p_{T} [GeV/c]", njptarray - 1, jetptxbins);
0215   histoid = "jet_pt_uncorr_" + flavor + "_tagged";
0216   histocontainer_["jet_pt_uncorr_" + flavor + "_tagged"] =
0217       fs->make<TH1D>(histoid.c_str(), "jet_tagged p_{T} [GeV/c]", njptarray - 1, jetptxbins);
0218   histoid = "jet_eta_" + flavor + "_tagged";
0219   histocontainer_["jet_eta_" + flavor + "_tagged"] =
0220       fs->make<TH1D>(histoid.c_str(), "jet_tagged #eta", netaarray - 1, jetetabins);
0221 
0222   //tagged min cut
0223   histoid = "jet_pt_" + flavor + "_taggedmincut";
0224   histocontainer_["jet_pt_" + flavor + "_taggedmincut"] =
0225       fs->make<TH1D>(histoid.c_str(), "jet_tagged p_{T} [GeV/c]", njptarray - 1, jetptxbins);
0226   histoid = "jet_pt_uncorr_" + flavor + "_taggedmincut";
0227   histocontainer_["jet_pt_uncorr_" + flavor + "_taggedmincut"] =
0228       fs->make<TH1D>(histoid.c_str(), "jet_tagged p_{T} [GeV/c]", njptarray - 1, jetptxbins);
0229   histoid = "jet_eta_" + flavor + "_taggedmincut";
0230   histocontainer_["jet_eta_" + flavor + "_taggedmincut"] =
0231       fs->make<TH1D>(histoid.c_str(), "jet_tagged #eta", netaarray - 1, jetetabins);
0232 
0233   //#tracks per jet
0234   histoid = "tracks_in_jet_" + flavor;
0235   histotitle = "traks per jet " + flavor;
0236   histocontainer_["tracks_in_jet_" + flavor] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), 31, -0.5, 30.5);
0237   histoid = "tracks_in_jet_" + flavor + "030";
0238   histotitle = "traks per jet " + flavor + "pt_{T} < 30[GeV/c]";
0239   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), 31, -0.5, 30.5);
0240   histoid = "tracks_in_jet_" + flavor + "3050";
0241   histotitle = "traks per jet " + flavor + "30 < pt_{T} < 50[GeV/c]";
0242   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), 31, -0.5, 30.5);
0243   histoid = "tracks_in_jet_" + flavor + "50";
0244   histotitle = "traks per jet " + flavor + "pt_{T} > 50[GeV/c]";
0245   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), 31, -0.5, 30.5);
0246 
0247   // pt of tracks in bins of jet pt 0-30,30-50,50
0248   histoid = "pt_tracks_in_jet_" + flavor;
0249   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0250   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0251   histoid = "pt_tracks_in_jet_" + flavor + "030";
0252   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0253   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0254   histoid = "pt_tracks_in_jet_" + flavor + "3050";
0255   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0256   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0257   histoid = "pt_tracks_in_jet_" + flavor + "50";
0258   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0259   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0260 
0261   // pt of tracks in bins of jet eta abs(eta)<1.5, 1.5<abs(eta)<2.4
0262   // combined with bins of jet pt
0263   histoid = "pt_tracks_in_jet_" + flavor + "_center";
0264   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0265   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0266   histoid = "pt_tracks_in_jet_" + flavor + "_center030";
0267   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0268   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0269   histoid = "pt_tracks_in_jet_" + flavor + "_center3050";
0270   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0271   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0272   histoid = "pt_tracks_in_jet_" + flavor + "_center50";
0273   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0274   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0275   histoid = "pt_tracks_in_jet_" + flavor + "_sides";
0276   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0277   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0278   histoid = "pt_tracks_in_jet_" + flavor + "_sides030";
0279   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0280   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0281   histoid = "pt_tracks_in_jet_" + flavor + "_sides3050";
0282   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0283   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0284   histoid = "pt_tracks_in_jet_" + flavor + "_sides50";
0285   histotitle = "track p_{T} [GeV/c] " + flavor + " jets";
0286   histocontainer_[histoid] = fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), ntptarray - 1, jetpttbins);
0287 
0288   //taggable 0,1,2
0289   std::map<int, std::string> tagbl;
0290   tagbl[0] = "0";
0291   tagbl[1] = "1";
0292   tagbl[2] = "2";
0293   for (size_t i = 0; i < tagbl.size(); i++) {
0294     histoid = "jet_pt_" + flavor + "_taggable" + tagbl[i];
0295     histotitle = flavor + " jet_taggable p_{T} [GeV/c]";
0296     histocontainer_["jet_pt_" + flavor + "_taggable" + tagbl[i]] =
0297         fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), njptarray - 1, jetptxbins);
0298     histoid = "jet_pt_uncorr_" + flavor + "_taggable" + tagbl[i];
0299     histotitle = flavor + " jet_taggable p_{T} uncorr [GeV/c]";
0300     histocontainer_["jet_pt_uncorr_" + flavor + "_taggable" + tagbl[i]] =
0301         fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), njptarray - 1, jetptxbins);
0302     histoid = "jet_eta_" + flavor + "_taggable" + tagbl[i];
0303     histotitle = flavor + " jet_taggable #eta";
0304     histocontainer_["jet_eta_" + flavor + "_taggable" + tagbl[i]] =
0305         fs->make<TH1D>(histoid.c_str(), histotitle.c_str(), netaarray - 1, jetetabins);
0306     histoid = "jet_eta_" + flavor + "_taggable" + tagbl[i] + "_030";
0307     histocontainer_["jet_eta_" + flavor + "_taggable" + tagbl[i] + "_030"] =
0308         fs->make<TH1D>(histoid.c_str(), "jet #eta", netaarray - 1, jetetabins);
0309     histoid = "jet_eta_" + flavor + "_taggable" + tagbl[i] + "_3050";
0310     histocontainer_["jet_eta_" + flavor + "_taggable" + tagbl[i] + "_3050"] =
0311         fs->make<TH1D>(histoid.c_str(), "jet #eta", netaarray - 1, jetetabins);
0312     histoid = "jet_eta_" + flavor + "_taggable" + tagbl[i] + "_50";
0313     histocontainer_["jet_eta_" + flavor + "_taggable" + tagbl[i] + "_50"] =
0314         fs->make<TH1D>(histoid.c_str(), "jet #eta", netaarray - 1, jetetabins);
0315   }
0316 
0317   histoid = "jet_scatter_pt_" + flavor;
0318   h2_["jet_scatter_pt_" + flavor] = fs->make<TH2D>(
0319       histoid.c_str(), "jet p_{T} uncorr(y) vs corr(x) [GeV/c]", njptarray - 1, jetptxbins, njptarray - 1, jetptxbins);
0320 
0321 }  //end function Set
0322 
0323 // ------------ method called once each job just before starting event loop  ------------
0324 // ------------              after setting histograms                --------------------
0325 // ------------  This function is needed to save histogram errors -----------------------
0326 void PatBTagCommonHistos::Sumw2() {
0327   for (std::map<std::string, TH1D*>::const_iterator ih = histocontainer_.begin(); ih != histocontainer_.end(); ++ih) {
0328     TH1D* htemp = (TH1D*)ih->second;
0329     htemp->Sumw2();
0330   }
0331 }  //end function Sumw2