Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2021-12-13 13:14:13

0001 //
0002 // Example to calculated top mass from genjets corrected on the fly with parton calibration
0003 // Top mass is calculated from:
0004 // 1) Uncorrected jets
0005 // 2) Corrected jets with flavour ttbar correction assuming jet mass = quark mass
0006 // 3) The same as 3 assuming massless jets
0007 // 4) Corrected jets with mixed ttbar correction assuming massless jets
0008 //
0009 // Author: Attilio Santocchia
0010 // Date: 15.06.2009
0011 // Tested on CMSSW_3_1_0
0012 //
0013 // user include files
0014 #include "FWCore/Framework/interface/Frameworkfwd.h"
0015 #include "FWCore/Framework/interface/one/EDAnalyzer.h"
0016 #include "FWCore/Framework/interface/ESHandle.h"
0017 #include "FWCore/Framework/interface/Event.h"
0018 #include "FWCore/Framework/interface/EventSetup.h"
0019 #include "FWCore/Framework/interface/MakerMacros.h"
0020 #include "FWCore/ParameterSet/interface/ParameterSet.h"
0021 #include "FWCore/ServiceRegistry/interface/Service.h"
0022 
0023 #include "DataFormats/Candidate/interface/Candidate.h"
0024 #include "DataFormats/Candidate/interface/CandidateFwd.h"
0025 #include "DataFormats/JetReco/interface/Jet.h"
0026 #include "DataFormats/JetReco/interface/JetCollection.h"
0027 #include "DataFormats/JetReco/interface/JetFloatAssociation.h"
0028 
0029 #include "JetMETCorrections/Objects/interface/JetCorrector.h"
0030 
0031 #include "DataFormats/Common/interface/View.h"
0032 #include "DataFormats/Common/interface/Ref.h"
0033 #include "DataFormats/Common/interface/getRef.h"
0034 
0035 #include "DataFormats/JetMatching/interface/JetFlavour.h"
0036 #include "DataFormats/JetMatching/interface/JetFlavourMatching.h"
0037 #include "DataFormats/JetMatching/interface/MatchedPartons.h"
0038 #include "DataFormats/JetMatching/interface/JetMatchedPartons.h"
0039 
0040 #include "CommonTools/UtilAlgos/interface/TFileService.h"
0041 
0042 #include "TFile.h"
0043 #include "TH1.h"
0044 #include "TF1.h"
0045 #include "TLorentzVector.h"
0046 
0047 // system include files
0048 #include <memory>
0049 #include <string>
0050 #include <iostream>
0051 #include <vector>
0052 #include <Math/VectorUtil.h>
0053 #include <TMath.h>
0054 
0055 class calcTopMass : public edm::one::EDAnalyzer<edm::one::SharedResources> {
0056 public:
0057   explicit calcTopMass(const edm::ParameterSet &);
0058   ~calcTopMass() override{};
0059   void analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) override;
0060 
0061 private:
0062   edm::InputTag sourcePartons_;
0063   edm::InputTag sourceByRefer_;
0064   edm::Handle<reco::JetMatchedPartonsCollection> theJetPartonMatch;
0065 
0066   std::string m_qT_CorrectorName;
0067   std::string m_cT_CorrectorName;
0068   std::string m_bT_CorrectorName;
0069   std::string m_tT_CorrectorName;
0070 
0071   float bMass;
0072   float cMass;
0073   float qMass;
0074 
0075   TH1F *hMassNoCorr;
0076   TH1F *hMassCorFl0;
0077   TH1F *hMassCorFlM;
0078   TH1F *hMassCorMix;
0079 };
0080 
0081 using namespace std;
0082 using namespace reco;
0083 using namespace edm;
0084 using namespace ROOT::Math::VectorUtil;
0085 
0086 calcTopMass::calcTopMass(const edm::ParameterSet &iConfig) {
0087   sourceByRefer_ = iConfig.getParameter<InputTag>("srcByReference");
0088   m_qT_CorrectorName = iConfig.getParameter<std::string>("qTopCorrector");
0089   m_cT_CorrectorName = iConfig.getParameter<std::string>("cTopCorrector");
0090   m_bT_CorrectorName = iConfig.getParameter<std::string>("bTopCorrector");
0091   m_tT_CorrectorName = iConfig.getParameter<std::string>("tTopCorrector");
0092 
0093   bMass = 4.5;
0094   cMass = 1.5;
0095   qMass = 0.3;
0096 
0097   usesResource(TFileService::kSharedResource);
0098   Service<TFileService> fs;
0099   hMassNoCorr = fs->make<TH1F>("hMassNoCorr", "", 100, 100, 300);
0100   hMassCorFl0 = fs->make<TH1F>("hMassCorFl0", "", 100, 100, 300);
0101   hMassCorFlM = fs->make<TH1F>("hMassCorFlM", "", 100, 100, 300);
0102   hMassCorMix = fs->make<TH1F>("hMassCorMix", "", 100, 100, 300);
0103 }
0104 
0105 void calcTopMass::analyze(const edm::Event &iEvent, const edm::EventSetup &iSetup) {
0106   cout << "[calcTopMass] analysing event " << iEvent.id() << endl;
0107   try {
0108     iEvent.getByLabel(sourceByRefer_, theJetPartonMatch);
0109   } catch (std::exception &ce) {
0110     cerr << "[calcTopMass] caught std::exception " << ce.what() << endl;
0111     return;
0112   }
0113 
0114   // get all correctors from top events
0115   const JetCorrector *qTopCorrector = JetCorrector::getJetCorrector(m_qT_CorrectorName, iSetup);
0116   const JetCorrector *cTopCorrector = JetCorrector::getJetCorrector(m_cT_CorrectorName, iSetup);
0117   const JetCorrector *bTopCorrector = JetCorrector::getJetCorrector(m_bT_CorrectorName, iSetup);
0118   const JetCorrector *tTopCorrector = JetCorrector::getJetCorrector(m_tT_CorrectorName, iSetup);
0119 
0120   TLorentzVector jetsNoCorr[6];
0121   TLorentzVector jetsCorFl0[6];
0122   TLorentzVector jetsCorFlM[6];
0123   TLorentzVector jetsCorMix[6];
0124 
0125   bool isQuarkFound[6] = {false};
0126 
0127   for (JetMatchedPartonsCollection::const_iterator j = theJetPartonMatch->begin(); j != theJetPartonMatch->end(); j++) {
0128     const math::XYZTLorentzVector theJet = (*j).first.get()->p4();
0129     const MatchedPartons aMatch = (*j).second;
0130     const GenParticleRef &thePhyDef = aMatch.physicsDefinitionParton();
0131 
0132     if (thePhyDef.isNonnull()) {
0133       int particPDG = thePhyDef.get()->pdgId();
0134 
0135       if (particPDG == 5) {  //b from t
0136         double bTcorr = bTopCorrector->correction(theJet);
0137         double tTcorr = tTopCorrector->correction(theJet);
0138         jetsNoCorr[0].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
0139         jetsCorFl0[0].SetPtEtaPhiM(theJet.pt() * bTcorr, theJet.eta(), theJet.phi(), 0);
0140         jetsCorFlM[0].SetPtEtaPhiM(theJet.pt() * bTcorr, theJet.eta(), theJet.phi(), bMass);
0141         jetsCorMix[0].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
0142         isQuarkFound[0] = true;
0143       } else if (particPDG == -5) {  //bbar from tbar
0144         double bTcorr = bTopCorrector->correction(theJet);
0145         double tTcorr = tTopCorrector->correction(theJet);
0146         jetsNoCorr[3].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
0147         jetsCorFl0[3].SetPtEtaPhiM(theJet.pt() * bTcorr, theJet.eta(), theJet.phi(), 0);
0148         jetsCorFlM[3].SetPtEtaPhiM(theJet.pt() * bTcorr, theJet.eta(), theJet.phi(), bMass);
0149         jetsCorMix[3].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
0150         isQuarkFound[3] = true;
0151       } else if (particPDG == 2) {  //W+ from t
0152         double qTcorr = qTopCorrector->correction(theJet);
0153         double tTcorr = tTopCorrector->correction(theJet);
0154         jetsNoCorr[1].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
0155         jetsCorFl0[1].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), 0);
0156         jetsCorFlM[1].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), qMass);
0157         jetsCorMix[1].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
0158         isQuarkFound[1] = true;
0159       } else if (particPDG == 4) {  //W+ from t
0160         double cTcorr = cTopCorrector->correction(theJet);
0161         double tTcorr = tTopCorrector->correction(theJet);
0162         jetsNoCorr[1].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
0163         jetsCorFl0[1].SetPtEtaPhiM(theJet.pt() * cTcorr, theJet.eta(), theJet.phi(), 0);
0164         jetsCorFlM[1].SetPtEtaPhiM(theJet.pt() * cTcorr, theJet.eta(), theJet.phi(), cMass);
0165         jetsCorMix[1].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
0166         isQuarkFound[1] = true;
0167       } else if (particPDG == -1) {  //W+ from t
0168         double qTcorr = qTopCorrector->correction(theJet);
0169         double tTcorr = tTopCorrector->correction(theJet);
0170         jetsNoCorr[2].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
0171         jetsCorFl0[2].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), 0);
0172         jetsCorFlM[2].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), qMass);
0173         jetsCorMix[2].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
0174         isQuarkFound[2] = true;
0175       } else if (particPDG == -3) {  //W+ from t
0176         double qTcorr = qTopCorrector->correction(theJet);
0177         double tTcorr = tTopCorrector->correction(theJet);
0178         jetsNoCorr[2].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
0179         jetsCorFl0[2].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), 0);
0180         jetsCorFlM[2].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), qMass);
0181         jetsCorMix[2].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
0182         isQuarkFound[2] = true;
0183       } else if (particPDG == -2) {  //W- from tbar
0184         double qTcorr = qTopCorrector->correction(theJet);
0185         double tTcorr = tTopCorrector->correction(theJet);
0186         jetsNoCorr[4].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
0187         jetsCorFl0[4].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), 0);
0188         jetsCorFlM[4].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), qMass);
0189         jetsCorMix[4].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
0190         isQuarkFound[4] = true;
0191       } else if (particPDG == -4) {  //W- from tbar
0192         double qTcorr = qTopCorrector->correction(theJet);
0193         double tTcorr = tTopCorrector->correction(theJet);
0194         jetsNoCorr[4].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
0195         jetsCorFl0[4].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), 0);
0196         jetsCorFlM[4].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), qMass);
0197         jetsCorMix[4].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
0198         isQuarkFound[4] = true;
0199       } else if (particPDG == 1) {  //W- from tbar
0200         double qTcorr = qTopCorrector->correction(theJet);
0201         double tTcorr = tTopCorrector->correction(theJet);
0202         jetsNoCorr[5].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
0203         jetsCorFl0[5].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), 0);
0204         jetsCorFlM[5].SetPtEtaPhiM(theJet.pt() * qTcorr, theJet.eta(), theJet.phi(), qMass);
0205         jetsCorMix[5].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
0206         isQuarkFound[5] = true;
0207       } else if (particPDG == 3) {  //W- from tbar
0208         double cTcorr = cTopCorrector->correction(theJet);
0209         double tTcorr = tTopCorrector->correction(theJet);
0210         jetsNoCorr[5].SetPtEtaPhiM(theJet.pt(), theJet.eta(), theJet.phi(), theJet.mass());
0211         jetsCorFl0[5].SetPtEtaPhiM(theJet.pt() * cTcorr, theJet.eta(), theJet.phi(), 0);
0212         jetsCorFlM[5].SetPtEtaPhiM(theJet.pt() * cTcorr, theJet.eta(), theJet.phi(), cMass);
0213         jetsCorMix[5].SetPtEtaPhiM(theJet.pt() * tTcorr, theJet.eta(), theJet.phi(), 0);
0214         isQuarkFound[5] = true;
0215       }
0216     }
0217   }
0218 
0219   if (isQuarkFound[0] && isQuarkFound[1] && isQuarkFound[2]) {
0220     TLorentzVector *theTopPNoCorr = new TLorentzVector(jetsNoCorr[0] + jetsNoCorr[1] + jetsNoCorr[2]);
0221     TLorentzVector *theTopPCorFl0 = new TLorentzVector(jetsCorFl0[0] + jetsCorFl0[1] + jetsCorFl0[2]);
0222     TLorentzVector *theTopPCorFlM = new TLorentzVector(jetsCorFlM[0] + jetsCorFlM[1] + jetsCorFlM[2]);
0223     TLorentzVector *theTopPCorMix = new TLorentzVector(jetsCorMix[0] + jetsCorMix[1] + jetsCorMix[2]);
0224     hMassNoCorr->Fill(theTopPNoCorr->M());
0225     hMassCorFl0->Fill(theTopPCorFl0->M());
0226     hMassCorFlM->Fill(theTopPCorFlM->M());
0227     hMassCorMix->Fill(theTopPCorMix->M());
0228   }
0229 
0230   if (isQuarkFound[3] && isQuarkFound[4] && isQuarkFound[5]) {
0231     TLorentzVector *theTopMNoCorr = new TLorentzVector(jetsNoCorr[3] + jetsNoCorr[4] + jetsNoCorr[5]);
0232     TLorentzVector *theTopMCorFl0 = new TLorentzVector(jetsCorFl0[3] + jetsCorFl0[4] + jetsCorFl0[5]);
0233     TLorentzVector *theTopMCorFlM = new TLorentzVector(jetsCorFlM[3] + jetsCorFlM[4] + jetsCorFlM[5]);
0234     TLorentzVector *theTopMCorMix = new TLorentzVector(jetsCorMix[3] + jetsCorMix[4] + jetsCorMix[5]);
0235     hMassNoCorr->Fill(theTopMNoCorr->M());
0236     hMassCorFl0->Fill(theTopMCorFl0->M());
0237     hMassCorFlM->Fill(theTopMCorFlM->M());
0238     hMassCorMix->Fill(theTopMCorMix->M());
0239   }
0240 }
0241 
0242 DEFINE_FWK_MODULE(calcTopMass);