JetCharge/test/JetChargeAnalyzer.cc

0001 #include "FWCore/Framework/interface/Frameworkfwd.h"
0002 #include "FWCore/Framework/interface/one/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::one::EDAnalyzer<edm::one::SharedResources> {
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* const 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   usesResource(TFileService::kSharedResource);
0057   edm::Service<TFileService> fs;
0058   TFileDirectory cwd = fs->mkdir(dir_.c_str());
0059   char buff[255], biff[255];
0060   for (int i = 0; i < 12; i++) {
0061     sprintf(biff, "jch_id_%s%d", (pdgIds[i] >= 0 ? "p" : "m"), abs(pdgIds[i]));
0062     sprintf(buff, "Jet charge for '%s' jets (pdgId %d) [ET > %f]", pdgs[i], pdgIds[i], minET_);
0063     charge_[i] = cwd.make<TH1D>(biff, buff, 22, -1.1, 1.1);
0064   }
0065 }
0066
0067 void JetChargeAnalyzer::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) {
0068   using namespace edm;
0069   using namespace reco;
0070
0071   Handle<JetChargeCollection> hJC;
0072   iEvent.getByToken(srcToken_, hJC);
0073   //
0074   // get jet flavour via genparticles
0075   //
0076   edm::Handle<JetFlavourMatchingCollection> jetMC;
0077   FlavourMap flavours;
0078
0079   iEvent.getByToken(jetMCSrcToken_, jetMC);
0080   for (JetFlavourMatchingCollection::const_iterator iter = jetMC->begin(); iter != jetMC->end(); iter++) {
0081     unsigned int fl = abs(iter->second.getFlavour());
0082     flavours.insert(FlavourMap::value_type(iter->first, fl));
0083   }
0084   //          for (JetChargeCollection::const_iterator it = hJC->begin(), ed = hJC->end(); it != ed; ++it) {
0085   for (unsigned int i = 0; i < hJC->size(); ++i) {
0086     const Jet& jet = *(hJC->key(i));
0087     if (jet.et() < minET_)
0088       continue;
0089     //        int id = jf_.identifyBasedOnPartons(jet).mainFlavour();
0090
0091     edm::RefToBase<reco::Jet> jetr = hJC->key(i);
0092     int id = flavours[jetr];
0093     int k;
0094     for (k = 0; k < 12; k++) {
0095       if (id == pdgIds[k])
0096         break;
0097     }
0098     if (k == 12) {
0099       std::cerr << "Error: jet with flavour " << id << ". !??" << std::endl;
0100       continue;
0101     }
0102     charge_[k]->Fill(hJC->value(i));
0103   }
0104 }
0105
0106 DEFINE_FWK_MODULE(JetChargeAnalyzer);