diff --git a/TopAnalysis/plugins/BoostedTTbarFlatTreeProducer.cc b/TopAnalysis/plugins/BoostedTTbarFlatTreeProducer.cc index 8079d104139..a2940799ac2 100644 --- a/TopAnalysis/plugins/BoostedTTbarFlatTreeProducer.cc +++ b/TopAnalysis/plugins/BoostedTTbarFlatTreeProducer.cc @@ -8,15 +8,22 @@ #include "TVectorD.h" #include "TMatrixDSym.h" #include "TMatrixDSymEigen.h" +#include "fastjet/contrib/Njettiness.hh" +#include "fastjet/tools/MassDropTagger.hh" +#include "fastjet/contrib/SoftDrop.hh" +#include "DataFormats/JetReco/interface/GenJet.h" +#include "DataFormats/JetReco/interface/GenJetCollection.h" -#include "KKousour/TopAnalysis/plugins/BoostedTTbarFlatTreeProducer.h" +#include "UserCode/TopAnalysis/plugins/BoostedTTbarFlatTreeProducer.h" using namespace std; using namespace reco; +using namespace fastjet; BoostedTTbarFlatTreeProducer::BoostedTTbarFlatTreeProducer(edm::ParameterSet const& cfg) { jetsToken = consumes(cfg.getParameter("jets")); + genjetsToken = consumes(cfg.getUntrackedParameter("genjets",edm::InputTag(""))); muonsToken = consumes(cfg.getParameter("muons")); electronsToken = consumes(cfg.getParameter("electrons")); metToken = consumes(cfg.getParameter("met")); @@ -32,6 +39,7 @@ BoostedTTbarFlatTreeProducer::BoostedTTbarFlatTreeProducer(edm::ParameterSet con runInfoToken = consumes(edm::InputTag("externalLHEProducer")); srcBtag_ = cfg.getParameter("btagger"); xmlFile_ = cfg.getParameter("xmlFile"); + xmlFileGen_ = cfg.getParameter("xmlFileGen"); triggerNames_ = cfg.getParameter >("triggerNames"); etaMax_ = cfg.getParameter("etaMax"); ptMin_ = cfg.getParameter("ptMin"); @@ -41,8 +49,22 @@ BoostedTTbarFlatTreeProducer::BoostedTTbarFlatTreeProducer(edm::ParameterSet con minMuPt_ = cfg.getParameter("minMuPt"); minElPt_ = cfg.getParameter("minElPt"); isMC_ = cfg.getUntrackedParameter("isMC",false); + isPrint_ = cfg.getUntrackedParameter("isPrint",false); saveWeights_ = cfg.getUntrackedParameter("saveWeights",false); debug_ = cfg.getUntrackedParameter("debug",false); + GenptMin_ = cfg.getParameter("GenptMin"); + GenetaMax_ = cfg.getParameter("GenetaMax"); + + jetFlavourInfosToken_ = consumes( cfg.getParameter("jetFlavourInfos")); + + //Gen Jet information + fAKJetDef = new fastjet::JetDefinition(fastjet::antikt_algorithm, 0.8); + int activeAreaRepeats = 1; + double ghostArea = 0.01; + double ghostEtaMax = 7.0; + fActiveArea = new fastjet::ActiveAreaSpec (ghostEtaMax,activeAreaRepeats,ghostArea); + fAreaDefinition = new fastjet::AreaDefinition (fastjet::active_area_explicit_ghosts, *fActiveArea ); + sd = new fastjet::contrib::SoftDrop(0.0,0.1,0.8);//beta_, zCut_, R0 ); } ////////////////////////////////////////////////////////////////////////////////////////// void BoostedTTbarFlatTreeProducer::beginJob() @@ -66,7 +88,9 @@ void BoostedTTbarFlatTreeProducer::beginJob() outTree_->Branch("lumi" ,&lumi_ ,"lumi_/I"); outTree_->Branch("nvtx" ,&nVtx_ ,"nVtx_/I"); outTree_->Branch("nJets" ,&nJets_ ,"nJets_/I"); + outTree_->Branch("nGenJets" ,&nGenJets_ ,"nGenJets_/I"); outTree_->Branch("nLeptons" ,&nLeptons_ ,"nLeptons_/I"); + outTree_->Branch("nGenLeptons" ,&nGenLeptons_ ,"nGenLeptons_/I"); outTree_->Branch("nBJets" ,&nBJets_ ,"nBJets_/I"); outTree_->Branch("pvRho" ,&pvRho_ ,"pvRho_/F"); outTree_->Branch("pvz" ,&pvz_ ,"pvz_/F"); @@ -75,17 +99,28 @@ void BoostedTTbarFlatTreeProducer::beginJob() outTree_->Branch("rho" ,&rho_ ,"rho_/F"); outTree_->Branch("ht" ,&ht_ ,"ht_/F"); outTree_->Branch("mva" ,&mva_ ,"mva_/F"); + outTree_->Branch("mvaGen" ,&mvaGen_ ,"mvaGen_/F"); outTree_->Branch("met" ,&met_ ,"met_/F"); + outTree_->Branch("metGen" ,&metGen_ ,"metGen_/F"); outTree_->Branch("metSig" ,&metSig_ ,"metSig_/F"); + outTree_->Branch("metGenSig" ,&metGenSig_ ,"metGenSig_/F"); outTree_->Branch("mJJ" ,&mJJ_ ,"mJJ_/F"); outTree_->Branch("yJJ" ,&yJJ_ ,"yJJ_/F"); outTree_->Branch("ptJJ" ,&ptJJ_ ,"ptJJ_/F"); outTree_->Branch("dRJJ" ,&dRJJ_ ,"dRJJ_/F"); outTree_->Branch("dPhiJJ" ,&dPhiJJ_ ,"dPhiJJ_/F"); outTree_->Branch("dPhiLJ" ,&dPhiLJ_ ,"dPhiLJ_/F"); + //Dijet GEN observables + outTree_->Branch("mGenJJ" ,&mGenJJ_ ,"mGenJJ_/F"); + outTree_->Branch("yGenJJ" ,&yGenJJ_ ,"yGenJJ_/F"); + outTree_->Branch("ptGenJJ" ,&ptGenJJ_ ,"ptGenJJ_/F"); + outTree_->Branch("dRGenJJ" ,&dRGenJJ_ ,"dRGenJJ_/F"); + outTree_->Branch("dPhiGenJJ" ,&dPhiGenJJ_ ,"dPhiGenJJ_/F"); + outTree_->Branch("dPhiGenLJ" ,&dPhiGenLJ_ ,"dPhiGenLJ_/F"); //------------------------------------------------------------------ isBtag_ = new std::vector; flavor_ = new std::vector; + flavorHadron_ = new std::vector; nSubJets_ = new std::vector; nBSubJets_ = new std::vector; pt_ = new std::vector; @@ -111,16 +146,28 @@ void BoostedTTbarFlatTreeProducer::beginJob() ptSub1_ = new std::vector; etaSub0_ = new std::vector; etaSub1_ = new std::vector; + phiSub0_ = new std::vector; + phiSub1_ = new std::vector; flavorSub0_ = new std::vector; flavorSub1_ = new std::vector; + flavorHadronSub0_ = new std::vector; + flavorHadronSub1_ = new std::vector; lId_ = new std::vector; lPt_ = new std::vector; lEta_ = new std::vector; lPhi_ = new std::vector; lE_ = new std::vector; lIso_ = new std::vector; + + lGenId_ = new std::vector; + lGenPt_ = new std::vector; + lGenEta_ = new std::vector; + lGenPhi_ = new std::vector; + lGenE_ = new std::vector; + outTree_->Branch("jetIsBtag" ,"vector" ,&isBtag_); outTree_->Branch("jetFlavor" ,"vector" ,&flavor_); + outTree_->Branch("jetFlavorHadron" ,"vector" ,&flavorHadron_); outTree_->Branch("jetNSub" ,"vector" ,&nSubJets_); outTree_->Branch("jetNBSub" ,"vector" ,&nBSubJets_); outTree_->Branch("jetPt" ,"vector" ,&pt_); @@ -146,14 +193,20 @@ void BoostedTTbarFlatTreeProducer::beginJob() outTree_->Branch("jetPtSub1" ,"vector" ,&ptSub1_); outTree_->Branch("jetEtaSub0" ,"vector" ,&etaSub0_); outTree_->Branch("jetEtaSub1" ,"vector" ,&etaSub1_); + outTree_->Branch("jetPhiSub0" ,"vector" ,&phiSub0_); + outTree_->Branch("jetPhiSub1" ,"vector" ,&phiSub1_); outTree_->Branch("jetFlavorSub0" ,"vector" ,&flavorSub0_); outTree_->Branch("jetFlavorSub1" ,"vector" ,&flavorSub1_); - outTree_->Branch("lepId" ,"vector" ,&lId_); - outTree_->Branch("lepPt" ,"vector" ,&lPt_); - outTree_->Branch("lepEta" ,"vector" ,&lEta_); - outTree_->Branch("lepPhi" ,"vector" ,&lPhi_); - outTree_->Branch("lepEnergy" ,"vector" ,&lE_); + outTree_->Branch("jetFlavorHadronSub0" ,"vector" ,&flavorHadronSub0_); + outTree_->Branch("jetFlavorHadronSub1" ,"vector" ,&flavorHadronSub1_); outTree_->Branch("lepIso" ,"vector" ,&lIso_); + + outTree_->Branch("lepGenId" ,"vector" ,&lGenId_); + outTree_->Branch("lepGenPt" ,"vector" ,&lGenPt_); + outTree_->Branch("lepGenEta" ,"vector" ,&lGenEta_); + outTree_->Branch("lepGenPhi" ,"vector" ,&lGenPhi_); + outTree_->Branch("lepGenEnergy" ,"vector" ,&lGenE_); + //------------------------------------------------------------------ triggerBit_ = new std::vector; triggerPre_ = new std::vector; @@ -161,6 +214,7 @@ void BoostedTTbarFlatTreeProducer::beginJob() outTree_->Branch("triggerPre" ,"vector" ,&triggerPre_); //------------------------------------------------------------------ discr_ = new BoostedDiscriminatorMVA("KKousour/TopAnalysis/data/"+xmlFile_); + discrGen_ = new BoostedDiscriminatorMVA("KKousour/TopAnalysis/data/"+xmlFileGen_); //------------------- MC --------------------------------- if (isMC_) { outTree_->Branch("decay" ,&decay_ ,"decay_/I"); @@ -194,6 +248,34 @@ void BoostedTTbarFlatTreeProducer::beginJob() outTree_->Branch("partonEta" ,"vector" ,&partonEta_); outTree_->Branch("partonPhi" ,"vector" ,&partonPhi_); outTree_->Branch("partonE" ,"vector" ,&partonE_); + + //gen jets + GenSoftDropMass_ = new std::vector; + GenSoftDropTau_ = new std::vector; + GenJetpt_ = new std::vector; + GenJetphi_ = new std::vector; + GenJeteta_ = new std::vector; + GenJetenergy_ = new std::vector; + GenJetmass_ = new std::vector; + GenJettau1_ = new std::vector; + GenJettau2_ = new std::vector; + GenJettau3_ = new std::vector; + GenSDSimmetry_ = new std::vector; + isBJetGen_ = new std::vector; + + outTree_->Branch("GenSoftDropMass" ,"vector" ,&GenSoftDropMass_); + outTree_->Branch("GenSoftDropTau" ,"vector" ,&GenSoftDropTau_); + outTree_->Branch("GenJetpt" ,"vector" ,&GenJetpt_); + outTree_->Branch("GenJeteta" ,"vector" ,&GenJeteta_); + outTree_->Branch("GenJettau1" ,"vector" ,&GenJettau1_); + outTree_->Branch("GenJettau2" ,"vector" ,&GenJettau2_); + outTree_->Branch("GenJettau3" ,"vector" ,&GenJettau3_); + outTree_->Branch("GenJetphi" ,"vector" ,&GenJetphi_); + outTree_->Branch("GenJetenergy" ,"vector" ,&GenJetenergy_); + outTree_->Branch("GenJetmass" ,"vector" ,&GenJetmass_); + outTree_->Branch("GenSDSimmetry" ,"vector" ,&GenSDSimmetry_); + outTree_->Branch("isBJetGen" ,"vector" ,&isBJetGen_); + } cout<<"Begin job finished"< vP4; @@ -498,9 +609,11 @@ void BoostedTTbarFlatTreeProducer::analyze(edm::Event const& iEvent, edm::EventS bool isBtag = (btag >= btagMin_); bool isLeptonMatched = false; float DRmax = 0.4; + if(isPrint_) {if (isBtag==true) cout<<"RECO "<pt()<<" "<eta()<<" "<phi()<eta(),lep->phi(),ijet->eta(),ijet->phi()) < DRmax ) isLeptonMatched = true; if (!isLeptonMatched) { flavor_ ->push_back(ijet->partonFlavour()); + flavorHadron_ ->push_back(ijet->hadronFlavour()); chf_ ->push_back(ijet->chargedHadronEnergyFraction()); nhf_ ->push_back(ijet->neutralHadronEnergyFraction()); phf_ ->push_back(ijet->photonEnergyFraction()); @@ -531,7 +644,9 @@ void BoostedTTbarFlatTreeProducer::analyze(edm::Event const& iEvent, edm::EventS massSub0_->push_back((ijet->subjets("SoftDrop"))[0]->mass()); ptSub0_->push_back((ijet->subjets("SoftDrop"))[0]->pt()); etaSub0_->push_back((ijet->subjets("SoftDrop"))[0]->eta()); + phiSub0_->push_back((ijet->subjets("SoftDrop"))[0]->phi()); flavorSub0_->push_back((ijet->subjets("SoftDrop"))[0]->partonFlavour()); + flavorHadronSub0_->push_back((ijet->subjets("SoftDrop"))[0]->hadronFlavour()); if ((ijet->subjets("SoftDrop"))[0]->bDiscriminator(srcBtag_.c_str()) >= btagMin_) { nBSub++; } @@ -540,7 +655,9 @@ void BoostedTTbarFlatTreeProducer::analyze(edm::Event const& iEvent, edm::EventS massSub1_->push_back((ijet->subjets("SoftDrop"))[1]->mass()); ptSub1_->push_back((ijet->subjets("SoftDrop"))[1]->pt()); etaSub1_->push_back((ijet->subjets("SoftDrop"))[1]->eta()); + phiSub1_->push_back((ijet->subjets("SoftDrop"))[1]->phi()); flavorSub1_->push_back((ijet->subjets("SoftDrop"))[1]->partonFlavour()); + flavorHadronSub1_->push_back((ijet->subjets("SoftDrop"))[1]->hadronFlavour()); if ((ijet->subjets("SoftDrop"))[1]->bDiscriminator(srcBtag_.c_str()) >= btagMin_) { nBSub++; } @@ -556,6 +673,16 @@ void BoostedTTbarFlatTreeProducer::analyze(edm::Event const& iEvent, edm::EventS if ((*met)[0].sumEt() > 0) { metSig_ = (*met)[0].et()/(*met)[0].sumEt(); } + + if(isMC_){ + metGen_ = (*met)[0].genMET()->et(); + if ((*met)[0].genMET()->sumEt() > 0) { + metGenSig_ = (*met)[0].genMET()->et()/(*met)[0].genMET()->sumEt(); + } + } + + if(isPrint_) cout<<"MET "<size(); run_ = iEvent.id().run(); evt_ = iEvent.id().event(); @@ -580,6 +707,7 @@ void BoostedTTbarFlatTreeProducer::analyze(edm::Event const& iEvent, edm::EventS iEvent.getByToken(lheEvtInfoToken,lheEvtInfo); iEvent.getByToken(genParticlesToken,genParticles); iEvent.getByToken(pupInfoToken,pupInfo); + iEvent.getByToken(genjetsToken,genjets); genEvtWeight_ = genEvtInfo->weight(); lheOriginalXWGTUP_ = lheEvtInfo->originalXWGTUP(); @@ -603,6 +731,9 @@ void BoostedTTbarFlatTreeProducer::analyze(edm::Event const& iEvent, edm::EventS LorentzVector p4T(0,0,0,0),p4Tbar(0,0,0,0); bool WPlusLep(false),WMinusLep(false); + + vector myGenLeptons; + for(unsigned ip = 0; ip < genParticles->size(); ++ ip) { const GenParticle &p = (*genParticles)[ip]; if (p.pdgId() == 24) { @@ -621,6 +752,19 @@ void BoostedTTbarFlatTreeProducer::analyze(edm::Event const& iEvent, edm::EventS } } } + //fill leptons for isolation criteria + if (abs(p.pdgId()) == 11 || abs(p.pdgId())==13) { + if(p.status()==3) + myGenLeptons.push_back(p.p4()); + } + + lGenPt_->push_back(p.pt()); + lGenEta_->push_back(p.eta()); + lGenPhi_->push_back(p.phi()); + lGenE_->push_back(p.energy()); + lGenId_->push_back(p.pdgId()); + //lIso_->push_back(getPFMiniIsolation(cands,myLeptons[ii])/myLeptons[ii]->pt()); + if (fabs(p.pdgId()) == 6 && p.isLastCopy()) { if (p.pdgId() == 6) { p4T = p.p4(); @@ -649,6 +793,9 @@ void BoostedTTbarFlatTreeProducer::analyze(edm::Event const& iEvent, edm::EventS } }// end of particle loop + + nGenLeptons_ = (int)myGenLeptons.size(); + if (p4T.pt() > p4Tbar.pt()) { ptTopParton_[0] = p4T.pt(); ptTopParton_[1] = p4Tbar.pt(); @@ -676,6 +823,91 @@ void BoostedTTbarFlatTreeProducer::analyze(edm::Event const& iEvent, edm::EventS npu_ = PUI->getTrueNumInteractions(); } } + + vector vP4Gen; + + for(GenJetCollection::const_iterator i_gen = genjets->begin(); i_gen != genjets->end(); i_gen++) { + if (i_gen->pt() > GenptMin_ && fabs(i_gen->y()) < GenetaMax_ ) { + + nGenJets_++; + + GenJetpt_->push_back(i_gen->pt()); + GenJetphi_->push_back(i_gen->phi()); + GenJeteta_->push_back(i_gen->y()); + GenJetenergy_->push_back(i_gen->energy()); + GenJetmass_->push_back(i_gen->mass()); + + vP4Gen.push_back(i_gen->p4()); + + std::vector lClusterParticles; + for(unsigned int ic=0; icnumberOfDaughters(); ic++) { + const reco::Candidate* gencand = i_gen->daughter(ic); + fastjet::PseudoJet pPart(gencand->px(),gencand->py(),gencand->pz(),gencand->energy()); + lClusterParticles.push_back(pPart); + } + + fClustering = new fastjet::ClusterSequenceArea(lClusterParticles, *fAKJetDef, *fAreaDefinition); + std::vector lOutJets = fClustering->inclusive_jets(20.0); + + if(lOutJets.size() == 0) { + delete fClustering; + return; + } + + fastjet::PseudoJet pT1JetSD = (*sd)( lOutJets[0]); + float iMSoftDrop = pT1JetSD.m(); + float iSimmetry = pT1JetSD.structure_of().symmetry(); + + fastjet::contrib::NormalizedMeasure normalizedMeasure (1.0,0.8); + fastjet::contrib::Njettiness routine(fastjet::contrib::Njettiness::onepass_kt_axes,normalizedMeasure); + float iTau1 = routine.getTau(1.,lClusterParticles); + float iTau2 = routine.getTau(2.,lClusterParticles); + float iTau3 = routine.getTau(3.,lClusterParticles); + + GenJettau1_->push_back(iTau1); + GenJettau2_->push_back(iTau2); + GenJettau3_->push_back(iTau3); + + GenSDSimmetry_->push_back(iSimmetry); + + float Tau32=iTau3/iTau2; + GenSoftDropTau_->push_back(Tau32); + GenSoftDropMass_->push_back(iMSoftDrop); + + lOutJets.clear(); + delete fClustering; + lClusterParticles.clear(); + + } + } + + edm::Handle theJetFlavourInfos; + iEvent.getByToken(jetFlavourInfosToken_, theJetFlavourInfos ); + + for ( reco::JetFlavourInfoMatchingCollection::const_iterator j = theJetFlavourInfos->begin();j != theJetFlavourInfos->end();++j ) { + const reco::Jet *aJet = (*j).first.get(); + reco::JetFlavourInfo aInfo = (*j).second; + + if (aJet->pt() > GenptMin_ && fabs(aJet->y()) < GenetaMax_ ) { + + int FlavourGenHadron = aInfo.getHadronFlavour(); + if(abs(FlavourGenHadron)==5){ + isBJetGen_->push_back(true); + if(isPrint_) cout<<"GEN "<pt()<<" "<y()<<" "<phi()<push_back(false); + } + } + + if (nGenJets_ > 1) { + dPhiGenJJ_ = fabs(deltaPhi(vP4Gen[0].phi(),vP4Gen[1].phi())); + if(nGenLeptons_>0) dPhiGenLJ_ = fabs(deltaPhi(vP4Gen[0].phi(),myGenLeptons[0].phi())); + dRGenJJ_ = deltaR(vP4Gen[0],vP4Gen[1]); + mGenJJ_ = (vP4Gen[0]+vP4Gen[1]).mass(); + yGenJJ_ = (vP4Gen[0]+vP4Gen[1]).Rapidity(); + ptGenJJ_ = (vP4Gen[0]+vP4Gen[1]).pt(); + mvaGen_ = discrGen_->eval((*GenJettau3_)[0]/(*GenJettau2_)[0],(*GenJettau3_)[0]/(*GenJettau1_)[0],(*GenJettau3_)[1]/(*GenJettau2_)[1],(*GenJettau3_)[1]/(*GenJettau1_)[1]); + } }//--- end of MC ------- cutFlowHisto_->Fill("All",1); @@ -699,23 +931,34 @@ void BoostedTTbarFlatTreeProducer::initialize() mJJ_ = -1; yJJ_ = -1; ptJJ_ = -1; + dPhiGenJJ_ = -1; + dPhiGenLJ_ = -1; + dRGenJJ_ = -1; + mGenJJ_ = -1; + yGenJJ_ = -1; + ptGenJJ_ = -1; run_ = -1; evt_ = -1; lumi_ = -1; nVtx_ = -1; nJets_ = -1; + nGenJets_ = -1; nLeptons_ = -1; nBJets_ = -1; rho_ = -1; met_ = -1; + metGen_ = -1; metSig_ = -1; + metGenSig_ = -1; ht_ = -1; mva_ = -999; + mvaGen_ = -999; pvRho_ = -999; pvz_ = -999; pvndof_ = -999; pvchi2_ = -999; flavor_ ->clear(); + flavorHadron_ ->clear(); nSubJets_ ->clear(); nBSubJets_ ->clear(); pt_ ->clear(); @@ -742,6 +985,8 @@ void BoostedTTbarFlatTreeProducer::initialize() etaSub1_ ->clear(); flavorSub0_ ->clear(); flavorSub1_ ->clear(); + flavorHadronSub0_ ->clear(); + flavorHadronSub1_ ->clear(); btag_ ->clear(); isBtag_ ->clear(); triggerBit_ ->clear(); @@ -752,6 +997,12 @@ void BoostedTTbarFlatTreeProducer::initialize() lEta_ ->clear(); lPhi_ ->clear(); lE_ ->clear(); + lGenId_ ->clear(); + lGenPt_ ->clear(); + lGenEta_ ->clear(); + lGenPhi_ ->clear(); + lGenE_ ->clear(); + //----- MC ------- if (isMC_) { decay_ = -1; @@ -777,6 +1028,19 @@ void BoostedTTbarFlatTreeProducer::initialize() partonEta_->clear(); partonPhi_->clear(); partonE_->clear(); + + isBJetGen_->clear(); + GenSoftDropMass_->clear(); + GenSoftDropTau_->clear(); + GenJetpt_->clear(); + GenJetphi_->clear(); + GenJeteta_->clear(); + GenJetenergy_->clear(); + GenJetmass_->clear(); + GenJettau1_->clear(); + GenJettau2_->clear(); + GenJettau3_->clear(); + } } ////////////////////////////////////////////////////////////////////////////////////////// @@ -785,19 +1049,3 @@ BoostedTTbarFlatTreeProducer::~BoostedTTbarFlatTreeProducer() } DEFINE_FWK_MODULE(BoostedTTbarFlatTreeProducer); - - - - - - - - - - - - - - - - diff --git a/TopAnalysis/plugins/BoostedTTbarFlatTreeProducer.h b/TopAnalysis/plugins/BoostedTTbarFlatTreeProducer.h index f231402737b..c1f0d011c7c 100644 --- a/TopAnalysis/plugins/BoostedTTbarFlatTreeProducer.h +++ b/TopAnalysis/plugins/BoostedTTbarFlatTreeProducer.h @@ -36,11 +36,25 @@ #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" #include "SimDataFormats/GeneratorProducts/interface/LHERunInfoProduct.h" #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h" -#include "KKousour/TopAnalysis/plugins/BoostedDiscriminatorMVA.h" +#include "UserCode/TopAnalysis/plugins/BoostedDiscriminatorMVA.h" #include "TTree.h" #include "TH1F.h" #include "TLorentzVector.h" +#include "SimDataFormats/JetMatching/interface/JetFlavourInfo.h"//add +#include "SimDataFormats/JetMatching/interface/JetFlavourInfoMatching.h"//add + +#include "fastjet/GhostedAreaSpec.hh" +#include "fastjet/ClusterSequenceArea.hh" +#include "fastjet/tools/Filter.hh" +#include "fastjet/tools/Pruner.hh" +#include "fastjet/GhostedAreaSpec.hh" +#include "fastjet/PseudoJet.hh" +#include "fastjet/tools/MassDropTagger.hh" +#include "fastjet/contrib/SoftDrop.hh" + +using namespace reco; + class BoostedTTbarFlatTreeProducer : public edm::EDAnalyzer { public: @@ -61,6 +75,7 @@ class BoostedTTbarFlatTreeProducer : public edm::EDAnalyzer void initialize(); //---- configurable parameters -------- edm::EDGetTokenT jetsToken; + edm::EDGetTokenT genjetsToken; edm::EDGetTokenT muonsToken; edm::EDGetTokenT electronsToken; edm::EDGetTokenT metToken; @@ -75,10 +90,10 @@ class BoostedTTbarFlatTreeProducer : public edm::EDAnalyzer edm::EDGetTokenT lheEvtInfoToken; edm::EDGetTokenT runInfoToken; - std::string srcBtag_,xmlFile_; + std::string srcBtag_,xmlFile_,xmlFileGen_; std::vector triggerNames_; - double etaMax_,ptMin_,ptMinLeading_,massMin_,btagMin_,minMuPt_,minElPt_; - bool isMC_,saveWeights_,debug_; + double etaMax_,ptMin_,ptMinLeading_,massMin_,btagMin_,minMuPt_,minElPt_,GenetaMax_,GenptMin_; + bool isMC_,isPrint_,saveWeights_,debug_; //--------------------------------------------------- edm::Service fs_; TTree *outTree_; @@ -87,24 +102,29 @@ class BoostedTTbarFlatTreeProducer : public edm::EDAnalyzer TH1F *triggerPassHisto_,*triggerNamesHisto_; //---- output TREE variables ------ //---- global event variables ----- - int run_,evt_,nVtx_,lumi_,nJets_,nBJets_,nLeptons_; - float rho_,met_,metSig_,ht_,mva_,pvRho_,pvz_,pvndof_,pvchi2_; + int run_,evt_,nVtx_,lumi_,nJets_,nBJets_,nLeptons_,nGenJets_,nGenLeptons_; + float rho_,met_,metSig_,ht_,mva_,pvRho_,pvz_,pvndof_,pvchi2_,mvaGen_,metGenSig_; std::vector *triggerBit_; std::vector *triggerPre_; //---- top variables -------------- float dRJJ_,dPhiJJ_,mJJ_,yJJ_,ptJJ_; - float dPhiLJ_; + float dRGenJJ_,dPhiGenJJ_,mGenJJ_,yGenJJ_,ptGenJJ_,metGen_; + float dPhiLJ_,dPhiGenLJ_; //---- jet variables -------------- std::vector *isBtag_; - std::vector *flavor_,*nSubJets_,*nBSubJets_; + std::vector *flavor_,*nSubJets_,*nBSubJets_,*flavorHadron_; std::vector *pt_,*eta_,*phi_,*mass_,*massSoftDrop_,*energy_,*chf_,*nhf_,*phf_,*elf_,*muf_,*btag_,*tau1_,*tau2_,*tau3_; - std::vector *btagSub0_,*btagSub1_,*massSub0_,*massSub1_,*ptSub0_,*ptSub1_,*etaSub0_,*etaSub1_; - std::vector *flavorSub0_,*flavorSub1_; + std::vector *btagSub0_,*btagSub1_,*massSub0_,*massSub1_,*ptSub0_,*ptSub1_,*etaSub0_,*etaSub1_,*phiSub0_,*phiSub1_; + std::vector *flavorSub0_,*flavorSub1_,*flavorHadronSub0_,*flavorHadronSub1_; //---- lepton variables ----------- std::vector *lId_; std::vector *lPt_,*lEta_,*lPhi_,*lE_,*lIso_; + + std::vector *lGenId_; + std::vector *lGenPt_,*lGenEta_,*lGenPhi_,*lGenE_; //---- MVA discriminator ---------- BoostedDiscriminatorMVA *discr_; + BoostedDiscriminatorMVA *discrGen_; //---- MC variables --------------- int npu_,decay_; float genEvtWeight_,lheOriginalXWGTUP_; @@ -114,7 +134,23 @@ class BoostedTTbarFlatTreeProducer : public edm::EDAnalyzer std::vector *partonId_,*partonSt_,*partonMatchIdx_; std::vector *partonPt_,*partonEta_,*partonPhi_,*partonE_,*partonMatchDR_; + edm::EDGetTokenT jetFlavourInfosToken_; + + //gen jets + std::vector *GenSoftDropMass_,*GenSoftDropTau_; + std::vector *GenJetpt_; + std::vector *GenJetphi_; + std::vector *GenJeteta_; + std::vector *GenJetenergy_; + std::vector *GenJetmass_; + std::vector *GenJettau1_; + std::vector *GenJettau2_; + std::vector *GenJettau3_; + std::vector *GenSDSimmetry_; + std::vector *isBJetGen_; + edm::Handle jets; + edm::Handle genjets; edm::Handle muons; edm::Handle electrons; edm::Handle met; @@ -128,6 +164,14 @@ class BoostedTTbarFlatTreeProducer : public edm::EDAnalyzer edm::Handle genEvtInfo; edm::Handle lheEvtInfo; edm::Handle runInfo; + + //fastjet + fastjet::Filter* fTrimmer1; + fastjet::JetDefinition* fAKJetDef; + fastjet::ActiveAreaSpec* fActiveArea; + fastjet::AreaDefinition* fAreaDefinition; + fastjet::ClusterSequenceArea* fClustering; + fastjet::contrib::SoftDrop* sd; }; diff --git a/TopAnalysis/test/TTbarMC_run_cfg.py b/TopAnalysis/test/TTbarMC_run_cfg.py new file mode 100644 index 00000000000..32b5bdbab63 --- /dev/null +++ b/TopAnalysis/test/TTbarMC_run_cfg.py @@ -0,0 +1,136 @@ +import FWCore.ParameterSet.Config as cms +process = cms.Process('myprocess') +process.TFileService=cms.Service("TFileService",fileName=cms.string('flatTreeFile3.root')) +#process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_cff') +#process.GlobalTag.globaltag = 'PHYS14_25_V1' +##-------------------- Define the source ---------------------------- +process.maxEvents = cms.untracked.PSet(input = cms.untracked.int32(-1)) +process.source = cms.Source("PoolSource", + fileNames = cms.untracked.vstring( + #"file:10B2302E-DA08-E611-884E-001E67DDD348.root") + #"file:004A0552-3929-E611-BD44-0025905A48F0.root"), + "file:004E1145-C828-E611-9097-0CC47A4D76AA.root"), +) +############# Format MessageLogger ################# +process.load('FWCore.MessageService.MessageLogger_cfi') +process.MessageLogger.cerr.FwkReport.reportEvery = 1000 + +from PhysicsTools.PatAlgos.selectionLayer1.jetSelector_cfi import selectedPatJets +process.goodJets = selectedPatJets.clone(src='slimmedJets',cut='pt>30 & abs(eta)<2.4') + +process.load('RecoJets.JetProducers.QGTagger_cfi') +process.QGTagger.srcJets = cms.InputTag('slimmedJets') +process.QGTagger.jetsLabel = cms.string('QGL_AK4PFchs') + +from RecoJets.Configuration.RecoGenJets_cff import ak4GenJets, ak8GenJets + +genParticleCollection = 'prunedGenParticles' + +from RecoJets.JetProducers.ak4GenJets_cfi import ak4GenJets +from RecoJets.JetProducers.GenJetParameters_cfi import * + +process.ak8GenJetsCustom = ak4GenJets.clone( + src = genParticleCollection, + rParam = cms.double(0.8), + jetAlgorithm = cms.string("AntiKt") +) + +genParticleCollection = 'prunedGenParticles' +genJetCollection = 'slimmedGenJetsAK8' + +from PhysicsTools.JetMCAlgos.HadronAndPartonSelector_cfi import selectedHadronsAndPartons +process.selectedHadronsAndPartons = selectedHadronsAndPartons.clone( + particles = genParticleCollection +) + +from PhysicsTools.JetMCAlgos.AK4PFJetsMCFlavourInfos_cfi import ak4JetFlavourInfos +process.ak8genJetFlavourInfos = ak4JetFlavourInfos.clone( + jets = genJetCollection, + rParam = cms.double(0.8), +) + +#only needed if information of the associated b hadrons are required +from PhysicsTools.JetMCAlgos.GenHFHadronMatcher_cff import matchGenBHadron +process.matchGenBHadron = matchGenBHadron.clone( + genParticles = genParticleCollection, + jetFlavourInfos = cms.InputTag("ak8genJetFlavourInfos"), + flavour = cms.int32(5), + onlyJetClusteredHadrons = cms.bool(True), + noBBbarResonances = cms.bool(False), +) + +##-------------------- User analyzer -------------------------------- +process.boosted = cms.EDAnalyzer('BoostedTTbarFlatTreeProducer', + jets = cms.InputTag('slimmedJetsAK8'), + muons = cms.InputTag('slimmedMuons'), + electrons = cms.InputTag('slimmedElectrons'), + met = cms.InputTag('slimmedMETs'), + vertices = cms.InputTag('offlineSlimmedPrimaryVertices'), + rho = cms.InputTag('fixedGridRhoFastjetAll'), + candidates = cms.InputTag('packedPFCandidates'), + triggerPrescales = cms.InputTag('patTrigger'), + xmlFile = cms.string('boosted_mva_Fisher.weights.xml'), + xmlFileGen = cms.string('boosted_mva_Fisher.weights.xml'), + nJetsMin = cms.int32(6), + nBJetsMin = cms.int32(2), + ptMin = cms.double(30), + minMuPt = cms.double(30), + minElPt = cms.double(30), + ptMinLeading = cms.double(300), + massMin = cms.double(100), + htMin = cms.double(400), + etaMax = cms.double(2.4), + kinfit = cms.string('kinFitTtFullHadEvent'), + btagMinThreshold = cms.double(0.814), + btagMin = cms.double(0.1), + btagMaxThreshold = cms.double(1.1), + btagger = cms.string('pfCombinedInclusiveSecondaryVertexV2BJetTags'), + qgtagger = cms.InputTag('QGTagger','qgLikelihood'), + pu = cms.untracked.string("addPileupInfo"), + genparticles = cms.untracked.InputTag('prunedGenParticles'), + triggerNames = cms.vstring('HLT_PFJet260_v','HLT_PFHT900_v'), + triggerResults = cms.InputTag('TriggerResults','','HLT'), + isMC = cms.untracked.bool(True), + genjets = cms.untracked.InputTag('slimmedGenJetsAK8'), + GenptMin = cms.double(30), + GenetaMax = cms.double(2.5), + jetFlavourInfos = cms.InputTag("ak8genJetFlavourInfos"), + isPrint = cms.untracked.bool(False), +) + +process.load('TopQuarkAnalysis.TopKinFitter.TtFullHadKinFitProducer_cfi') + +process.kinFitTtFullHadEvent.jets = 'goodJets' +process.kinFitTtFullHadEvent.bTagAlgo = 'pfCombinedInclusiveSecondaryVertexV2BJetTags' +process.kinFitTtFullHadEvent.minBTagValueBJet = 0.814 +process.kinFitTtFullHadEvent.maxBTagValueNonBJet = 0.814 +process.kinFitTtFullHadEvent.bTags = 2 +process.kinFitTtFullHadEvent.maxNJets = 8 + +#process.ttFilter = cms.EDFilter('AllHadronicPartonFilter', +# genparticles = cms.InputTag('prunedGenParticles'), +# forceTopDecay = cms.bool(True), +# forceHiggsDecay = cms.bool(False) +#) + +process.kinFitTtFullHadEventNoW = process.kinFitTtFullHadEvent.clone(constraints = cms.vuint32(5)) +process.boostedNoW = process.boosted.clone(kinfit = 'kinFitTtFullHadEventNoW') + +process.p = cms.Path( + #process.ttFilter * + process.selectedHadronsAndPartons*process.ak8genJetFlavourInfos*#*process.matchGenBHadron* + process.goodJets * + #process.QGTagger * + #process.kinFitTtFullHadEventNoW * + #process.boostedNoW * + process.kinFitTtFullHadEvent * + process.boosted +) + + + + + + + +