JetCharge/test/JetChargeAnalyzer.cc

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/EDAnalyzer.h"
0003 #include "FWCore/Framework/interface/Event.h"
0004 #include "FWCore/Framework/interface/MakerMacros.h"
0005
0006 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0007 #include "FWCore/Utilities/interface/InputTag.h"
0008
0009 #include "DataFormats/JetReco/interface/Jet.h"
0010 #include "DataFormats/JetReco/interface/JetFloatAssociation.h"
0011
0012 #include "DataFormats/JetMatching/interface/JetFlavour.h"
0013 #include "DataFormats/JetMatching/interface/JetFlavourMatching.h"
0014
0015 //#include "RecoBTag/MCTools/interface/JetFlavourIdentifier.h"
0016
0017 #include "FWCore/ServiceRegistry/interface/Service.h"
0018 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0019 #include <TH1.h>
0020 #include <string>
0021
0022 class JetChargeAnalyzer : public edm::EDAnalyzer {
0023 public:
0024   struct JetRefCompare {
0025     inline bool operator()(const edm::RefToBase<reco::Jet>& j1, const edm::RefToBase<reco::Jet>& j2) const {
0026       return j1.id() < j2.id() || (j1.id() == j2.id() && j1.key() < j2.key());
0027     }
0028   };
0029
0030   typedef std::map<edm::RefToBase<reco::Jet>, unsigned int, JetRefCompare> FlavourMap;
0031   typedef reco::JetFloatAssociation::Container JetChargeCollection;
0032
0033   explicit JetChargeAnalyzer(const edm::ParameterSet&);
0034   ~JetChargeAnalyzer() {}
0035
0036   virtual void analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) override;
0037
0038 private:
0039   // physics stuff
0040   edm::EDGetTokenT<JetChargeCollection> srcToken_;
0041   edm::EDGetTokenT<reco::JetFlavourMatchingCollection> jetMCSrcToken_;
0042   double minET_;
0043   // plot stuff
0044   std::string dir_;
0045   TH1D* charge_[12];
0046 };
0047
0048 const int pdgIds[12] = {0, 1, -1, 2, -2, 3, -3, 4, -4, 5, -5, 21};
0049 const char* pdgs[12] = {"?", "u", "-u", "d", "-d", "s", "-s", "c", "-c", "b", "-b", "g"};
0050
0051 JetChargeAnalyzer::JetChargeAnalyzer(const edm::ParameterSet& iConfig)
0052     : srcToken_(consumes<JetChargeCollection>(iConfig.getParameter<edm::InputTag>("src"))),
0053       jetMCSrcToken_(consumes<reco::JetFlavourMatchingCollection>(iConfig.getParameter<edm::InputTag>("jetFlavour"))),
0054       minET_(iConfig.getParameter<double>("minET")),
0055       dir_(iConfig.getParameter<std::string>("dir")) {
0056   edm::Service<TFileService> fs;
0057   TFileDirectory cwd = fs->mkdir(dir_.c_str());
0058   char buff[255], biff[255];
0059   for (int i = 0; i < 12; i++) {
0060     sprintf(biff, "jch_id_%s%d", (pdgIds[i] >= 0 ? "p" : "m"), abs(pdgIds[i]));
0061     sprintf(buff, "Jet charge for '%s' jets (pdgId %d) [ET > %f]", pdgs[i], pdgIds[i], minET_);
0062     charge_[i] = cwd.make<TH1D>(biff, buff, 22, -1.1, 1.1);
0063   }
0064 }
0065
0066 void JetChargeAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0067   using namespace edm;
0068   using namespace reco;
0069
0070   Handle<JetChargeCollection> hJC;
0071   iEvent.getByToken(srcToken_, hJC);
0072   //
0073   // get jet flavour via genparticles
0074   //
0075   edm::Handle<JetFlavourMatchingCollection> jetMC;
0076   FlavourMap flavours;
0077
0078   iEvent.getByToken(jetMCSrcToken_, jetMC);
0079   for (JetFlavourMatchingCollection::const_iterator iter = jetMC->begin(); iter != jetMC->end(); iter++) {
0080     unsigned int fl = abs(iter->second.getFlavour());
0081     flavours.insert(FlavourMap::value_type(iter->first, fl));
0082   }
0083   //          for (JetChargeCollection::const_iterator it = hJC->begin(), ed = hJC->end(); it != ed; ++it) {
0084   for (unsigned int i = 0; i < hJC->size(); ++i) {
0085     const Jet& jet = *(hJC->key(i));
0086     if (jet.et() < minET_)
0087       continue;
0088     //        int id = jf_.identifyBasedOnPartons(jet).mainFlavour();
0089
0090     edm::RefToBase<reco::Jet> jetr = hJC->key(i);
0091     int id = flavours[jetr];
0092     int k;
0093     for (k = 0; k < 12; k++) {
0094       if (id == pdgIds[k])
0095         break;
0096     }
0097     if (k == 12) {
0098       std::cerr << "Error: jet with flavour " << id << ". !??" << std::endl;
0099       continue;
0100     }
0101     charge_[k]->Fill(hJC->value(i));
0102   }
0103 }
0104
0105 DEFINE_FWK_MODULE(JetChargeAnalyzer);