diff --git a/NtupleProducer/interface/GenHelper.h b/NtupleProducer/interface/GenHelper.h index 35c9cf5b..191dd680 100644 --- a/NtupleProducer/interface/GenHelper.h +++ b/NtupleProducer/interface/GenHelper.h @@ -10,6 +10,9 @@ #ifndef GenHelper_h #define GenHelper_h +#include "TVector3.h" +#include "TLorentzVector.h" + #include #include @@ -50,8 +53,18 @@ namespace genhelper { WDecay GetTopDecay (const reco::Candidate* part); // return final state for top (= final state for W) -> see enum for code reco::GenParticle GetTauHad (const reco::Candidate* part); // build had tau by summing sons without nu + reco::GenParticle GetTauHadNeutrals (const reco::Candidate* part); // build neutral component of had tau by summing sons without nu + const reco::Candidate* IsFromID (const reco::Candidate* part, int targetPDGId); // find if is son of a certain particle (select by targetPDGId); if not found, return NULL, else return its pointer int GetIndexInOutput (const reco::Candidate* part, std::vector cands); -} + typedef reco::GenParticleCollection::const_iterator IG; + typedef reco::GenParticleRefVector::const_iterator IGR; + TVector3 ImpactParameter(const TVector3& pv, const TVector3& sv, const TLorentzVector& p4);//Calculate generator level impact parameter + void GetTausDaughters(const reco::GenParticle& tau, reco::GenParticleRefVector& products, bool ignoreNus, bool direct); + void FindDescendents(const reco::GenParticle& base, reco::GenParticleRefVector& descendents, int status, int pdgId=0, bool skipPhotonsPi0AndFSR=false); + const reco::GenParticleRef GetLeadChParticle(const reco::GenParticleRefVector& products); + int getDetailedTauDecayMode(const reco::GenParticleRefVector& products); + +} #endif diff --git a/NtupleProducer/plugins/BuildFile.xml b/NtupleProducer/plugins/BuildFile.xml index 6227e723..fe19abfa 100644 --- a/NtupleProducer/plugins/BuildFile.xml +++ b/NtupleProducer/plugins/BuildFile.xml @@ -15,6 +15,7 @@ + diff --git a/NtupleProducer/plugins/GenFiller.cc b/NtupleProducer/plugins/GenFiller.cc index 3ddeec77..969b682d 100644 --- a/NtupleProducer/plugins/GenFiller.cc +++ b/NtupleProducer/plugins/GenFiller.cc @@ -213,6 +213,29 @@ void GenFiller::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) if (DEBUG) cout << " --> fromTau: 1, indexTau: " << genhelper::GetIndexInOutput(MothPtr, cands_) << endl; } + + if(filtGenP.hasUserInt("tauGenDecayMode") && + (filtGenP.hasUserInt("ZMothIndex") || filtGenP.hasUserInt("HMothIndex") || filtGenP.hasUserInt("MSSMHMothIndex"))){ + + TVector3 aPVGenPoint = TVector3(genPClone->vx(), genPClone->vy(), genPClone->vz()); + reco::GenParticleRefVector tauDaughters; + genhelper::GetTausDaughters(*genPClone,tauDaughters,true,false); + int detailedDecayMode = genhelper::getDetailedTauDecayMode(tauDaughters); + filtGenP.addUserInt("tauGenDetailedDecayMode", detailedDecayMode); + + reco::GenParticleRef leadChParticleRef = genhelper::GetLeadChParticle(tauDaughters); + + TLorentzVector p4LeadingChParticle(leadChParticleRef->px(), + leadChParticleRef->py(), + leadChParticleRef->pz(), + leadChParticleRef->energy()); + TVector3 tauDecayVertex(leadChParticleRef->vx(), leadChParticleRef->vy(), leadChParticleRef->vz()); + TVector3 pca = genhelper::ImpactParameter(aPVGenPoint, tauDecayVertex, p4LeadingChParticle); + filtGenP.addUserFloat("pca_x",pca.X()); + filtGenP.addUserFloat("pca_y",pca.Y()); + filtGenP.addUserFloat("pca_z",pca.Z()); + } + result->push_back (filtGenP); } @@ -222,24 +245,51 @@ void GenFiller::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) for (unsigned int iTauH = 0; iTauH < tauHadcandsMothers_.size(); iTauH++) { int tauMothInd = tauHadcandsMothers_.at(iTauH); - pat::GenericParticle tauH (genhelper::GetTauHad(cands_.at(tauMothInd))); + pat::GenericParticle tauH (genhelper::GetTauHad(cands_.at(tauMothInd))); + pat::GenericParticle tauH_neutral (genhelper::GetTauHadNeutrals(cands_.at(tauMothInd))); tauH.addUserInt ("TauMothIndex", tauMothInd); + tauH_neutral.addUserInt ("TauMothIndex", tauMothInd); // copy all the other flags from original tau pat::GenericParticle& tauMothGenP = result->at(tauMothInd); - if (tauMothGenP.hasUserInt("HMothIndex") ) tauH.addUserInt ("HMothIndex", tauMothGenP.userInt ("HMothIndex")); - if (tauMothGenP.hasUserInt("MSSMHMothIndex") ) tauH.addUserInt ("MSSMHMothIndex", tauMothGenP.userInt ("MSSMHMothIndex")); - if (tauMothGenP.hasUserInt("TopMothIndex") ) tauH.addUserInt ("TopMothIndex", tauMothGenP.userInt ("TopMothIndex")); - if (tauMothGenP.hasUserInt("bMothIndex") ) tauH.addUserInt ("bMothIndex", tauMothGenP.userInt ("bMothIndex")); - if (tauMothGenP.hasUserInt("WMothIndex") ) tauH.addUserInt ("WMothIndex", tauMothGenP.userInt ("WMothIndex")); - if (tauMothGenP.hasUserInt("ZMothIndex") ) tauH.addUserInt ("ZMothIndex", tauMothGenP.userInt ("ZMothIndex")); + if (tauMothGenP.hasUserInt("HMothIndex") ){ + tauH.addUserInt ("HMothIndex", tauMothGenP.userInt ("HMothIndex")); + tauH_neutral.addUserInt ("HMothIndex", tauMothGenP.userInt ("HMothIndex")); + } + if (tauMothGenP.hasUserInt("MSSMHMothIndex") ){ + tauH.addUserInt ("MSSMHMothIndex", tauMothGenP.userInt ("MSSMHMothIndex")); + tauH_neutral.addUserInt ("MSSMHMothIndex", tauMothGenP.userInt ("MSSMHMothIndex")); + } + if (tauMothGenP.hasUserInt("TopMothIndex") ){ + tauH.addUserInt ("TopMothIndex", tauMothGenP.userInt ("TopMothIndex")); + tauH_neutral.addUserInt ("TopMothIndex", tauMothGenP.userInt ("TopMothIndex")); + } + if (tauMothGenP.hasUserInt("bMothIndex") ){ + tauH.addUserInt ("bMothIndex", tauMothGenP.userInt ("bMothIndex")); + tauH_neutral.addUserInt ("bMothIndex", tauMothGenP.userInt ("bMothIndex")); + } + if (tauMothGenP.hasUserInt("WMothIndex") ){ + tauH.addUserInt ("WMothIndex", tauMothGenP.userInt ("WMothIndex")); + tauH_neutral.addUserInt ("WMothIndex", tauMothGenP.userInt ("WMothIndex")); + } + if (tauMothGenP.hasUserInt("ZMothIndex") ){ + tauH.addUserInt ("ZMothIndex", tauMothGenP.userInt ("ZMothIndex")); + tauH_neutral.addUserInt ("ZMothIndex", tauMothGenP.userInt ("ZMothIndex")); + } + + // many flags change of meaning w.r.t. mother tau, put everything to 0 (can be changed in future) int tauhFlags = 0; - tauH.addUserInt ("generalGenFlags", tauhFlags); // remember! TauH inherits ALL the flags from + tauH.addUserInt ("generalGenFlags", tauhFlags); // remember! TauH inherits ALL the flags from + tauH_neutral.addUserInt ("generalGenFlags", tauhFlags); // remember! TauH inherits ALL the flags from - if (DEBUG) cout << " ++ " << iTauH << " id: " << tauH.pdgId() << " | pt: " << tauH.pt() << " | eta: " << tauH.eta() << endl; + if (DEBUG){ + cout << " ++ " << iTauH << " id: " << tauH.pdgId() << " | pt: " << tauH.pt() << " | eta: " << tauH.eta() << endl; + cout << " ++ " << iTauH << " id: " << tauH_neutral.pdgId() << " | pt: " << tauH_neutral.pt() << " | eta: " << tauH_neutral.eta() << endl; + } result->push_back (tauH); + result->push_back (tauH_neutral); } iEvent.put(result); } diff --git a/NtupleProducer/plugins/HTauTauNtuplizer.cc b/NtupleProducer/plugins/HTauTauNtuplizer.cc index 59f400e8..85b14eab 100644 --- a/NtupleProducer/plugins/HTauTauNtuplizer.cc +++ b/NtupleProducer/plugins/HTauTauNtuplizer.cc @@ -103,6 +103,15 @@ #include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h" #include "JetMETCorrections/Objects/interface/JetCorrectionsRecord.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h" +#include "TrackingTools/IPTools/interface/IPTools.h" +#include "RecoVertex/VertexPrimitives/interface/ConvertToFromReco.h" +#include "TrackingTools/GeomPropagators/interface/AnalyticalTrajectoryExtrapolatorToLine.h" +#include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h" + #include "TLorentzVector.h" namespace { @@ -110,6 +119,7 @@ bool writeJets = true; // Write jets in the tree. bool writeFatJets = true; bool writeSoftLep = false; + bool doPVRefit = true; bool DEBUG = false; } @@ -146,7 +156,7 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { //virtual void FillPhoton(const pat::Photon& photon); int FillJet(const edm::View* jet, const edm::Event&, JetCorrectionUncertainty*); void FillFatJet(const edm::View* fatjets, const edm::Event&); - void FillSoftLeptons(const edm::View *dauhandler, const edm::Event& event, bool theFSR, const edm::View* jets); + void FillSoftLeptons(const edm::View *dauhandler, const edm::Event& event, const edm::EventSetup& setup, bool theFSR, const edm::View* jets); //void FillbQuarks(const edm::Event&); void FillGenInfo(const edm::Event&); void FillGenJetInfo(const edm::Event&); @@ -157,6 +167,11 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { static bool ComparePairsbyPt(pat::CompositeCandidate i, pat::CompositeCandidate j); static bool ComparePairsbyIso(pat::CompositeCandidate i, pat::CompositeCandidate j); + bool refitPV(const edm::Event & iEvent, const edm::EventSetup & iSetup); + bool findPrimaryVertices(const edm::Event & iEvent, const edm::EventSetup & iSetup); + TVector3 getPCA(const edm::Event & iEvent, const edm::EventSetup & iSetup, + const reco::Track *aTrack, const GlobalPoint & aPoint); + // ----------member data --------------------------- //std::map genFlagPosMap_; // to convert from input to output enum format for H/Z decays @@ -206,7 +221,7 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { edm::EDGetTokenT theTotTag; edm::EDGetTokenT thePassTag; edm::EDGetTokenT theLHEPTag; - + edm::EDGetTokenT beamSpotTag; //flags //static const int nOutVars =14; @@ -245,6 +260,11 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { Float_t _PUReweight; Float_t _rho; Int_t _nup; + + Float_t _pv_x=0, _pv_y=0, _pv_z=0; + Float_t _pvGen_x=0, _pvGen_y=0, _pvGen_z=0; + Float_t _pvRefit_x=0, _pvRefit_y=0, _pvRefit_z=0; + bool _isRefitPV=false; // pairs //std::vector _mothers; @@ -261,6 +281,17 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { std::vector _daughters_py; std::vector _daughters_pz; std::vector _daughters_e; + + std::vector _daughters_charged_px; + std::vector _daughters_charged_py; + std::vector _daughters_charged_pz; + std::vector _daughters_charged_e; + + std::vector _daughters_neutral_px; + std::vector _daughters_neutral_py; + std::vector _daughters_neutral_pz; + std::vector _daughters_neutral_e; + std::vector _daughters_TauUpExists; std::vector _daughters_px_TauUp; std::vector _daughters_py_TauUp; @@ -287,6 +318,11 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { std::vector _genpart_py; std::vector _genpart_pz; std::vector _genpart_e; + + std::vector _genpart_pca_x; + std::vector _genpart_pca_y; + std::vector _genpart_pca_z; + std::vector _genpart_pdg; std::vector _genpart_status; //std::vector _genpart_mothInd; @@ -301,6 +337,7 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { std::vector _genpart_WDecayMode; std::vector _genpart_TopDecayMode; std::vector _genpart_TauGenDecayMode; + std::vector _genpart_TauGenDetailedDecayMode; std::vector _genpart_flags; // vector of bit flags bout gen info @@ -484,6 +521,19 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { std::vector _daughters_jetBTagCSV; std::vector _daughters_lepMVA_mvaId; + std::vector _daughters_pca_x; + std::vector _daughters_pca_y; + std::vector _daughters_pca_z; + + std::vector _daughters_pcaRefitPV_x; + std::vector _daughters_pcaRefitPV_y; + std::vector _daughters_pcaRefitPV_z; + + std::vector _daughters_pcaGenPV_x; + std::vector _daughters_pcaGenPV_y; + std::vector _daughters_pcaGenPV_z; + + //Jets variables Int_t _numberOfJets; //std::vector _jets; @@ -577,7 +627,8 @@ HTauTauNtuplizer::HTauTauNtuplizer(const edm::ParameterSet& pset) : reweight(), theGenJetTag (consumes> (pset.getParameter("genjetCollection"))), theTotTag (consumes (pset.getParameter("totCollection"))), thePassTag (consumes (pset.getParameter("passCollection"))), - theLHEPTag (consumes (pset.getParameter("lhepCollection"))) + theLHEPTag (consumes (pset.getParameter("lhepCollection"))), + beamSpotTag (consumes (pset.getParameter("beamSpot"))) { @@ -664,6 +715,17 @@ void HTauTauNtuplizer::Initialize(){ _daughters_py.clear(); _daughters_pz.clear(); _daughters_e.clear(); + + _daughters_charged_px.clear(); + _daughters_charged_py.clear(); + _daughters_charged_pz.clear(); + _daughters_charged_e.clear(); + + _daughters_neutral_px.clear(); + _daughters_neutral_py.clear(); + _daughters_neutral_pz.clear(); + _daughters_neutral_e.clear(); + _daughters_TauUpExists.clear(); _daughters_px_TauUp.clear(); _daughters_py_TauUp.clear(); @@ -741,6 +803,18 @@ void HTauTauNtuplizer::Initialize(){ //_daughter2.clear(); _softLeptons.clear(); //_genDaughters.clear(); + + _daughters_pca_x.clear(); + _daughters_pca_y.clear(); + _daughters_pca_z.clear(); + + _daughters_pcaRefitPV_x.clear(); + _daughters_pcaRefitPV_y.clear(); + _daughters_pcaRefitPV_z.clear(); + + _daughters_pcaGenPV_x.clear(); + _daughters_pcaGenPV_y.clear(); + _daughters_pcaGenPV_z.clear(); //_bquarks.clear(); //_bquarks_px.clear(); @@ -754,6 +828,11 @@ void HTauTauNtuplizer::Initialize(){ _genpart_py.clear(); _genpart_pz.clear(); _genpart_e.clear(); + + _genpart_pca_x.clear(); + _genpart_pca_y.clear(); + _genpart_pca_z.clear(); + _genpart_pdg.clear(); _genpart_status.clear(); //_genpart_mothInd.clear(); @@ -768,6 +847,7 @@ void HTauTauNtuplizer::Initialize(){ _genpart_TopDecayMode.clear(); _genpart_WDecayMode.clear(); _genpart_TauGenDecayMode.clear(); + _genpart_TauGenDetailedDecayMode.clear(); _genpart_flags.clear(); _genjet_px.clear(); @@ -976,6 +1056,16 @@ void HTauTauNtuplizer::beginJob(){ myTree->Branch("daughters_e",&_daughters_e); myTree->Branch("daughters_charge",&_daughters_charge); + myTree->Branch("daughters_charged_px",&_daughters_charged_px); + myTree->Branch("daughters_charged_py",&_daughters_charged_py); + myTree->Branch("daughters_charged_pz",&_daughters_charged_pz); + myTree->Branch("daughters_charged_e",&_daughters_charged_e); + + myTree->Branch("daughters_neutral_px",&_daughters_neutral_px); + myTree->Branch("daughters_neutral_py",&_daughters_neutral_py); + myTree->Branch("daughters_neutral_pz",&_daughters_neutral_pz); + myTree->Branch("daughters_neutral_e",&_daughters_neutral_e); + myTree->Branch("daughters_TauUpExists",&_daughters_TauUpExists); myTree->Branch("daughters_px_TauUp",&_daughters_px_TauUp); myTree->Branch("daughters_py_TauUp",&_daughters_py_TauUp); @@ -1012,6 +1102,9 @@ void HTauTauNtuplizer::beginJob(){ myTree->Branch("genpart_py", &_genpart_py); myTree->Branch("genpart_pz", &_genpart_pz); myTree->Branch("genpart_e", &_genpart_e); + myTree->Branch("genpart_pca_x",&_genpart_pca_x); + myTree->Branch("genpart_pca_y",&_genpart_pca_y); + myTree->Branch("genpart_pca_z",&_genpart_pca_z); myTree->Branch("genpart_pdg", &_genpart_pdg); myTree->Branch("genpart_status", &_genpart_status); //myTree->Branch("genpart_mothInd", _genpart_mothInd); @@ -1026,6 +1119,7 @@ void HTauTauNtuplizer::beginJob(){ myTree->Branch("genpart_TopDecayMode", &_genpart_TopDecayMode); myTree->Branch("genpart_WDecayMode", &_genpart_WDecayMode); myTree->Branch("genpart_TauGenDecayMode", &_genpart_TauGenDecayMode); + myTree->Branch("genpart_TauGenDetailedDecayMode", &_genpart_TauGenDetailedDecayMode); myTree->Branch("genpart_flags", &_genpart_flags); myTree->Branch("genjet_px", &_genjet_px); @@ -1169,7 +1263,17 @@ void HTauTauNtuplizer::beginJob(){ myTree->Branch("daughters_jetPtRatio",&_daughters_jetPtRatio); myTree->Branch("daughters_jetBTagCSV",&_daughters_jetBTagCSV); myTree->Branch("daughters_lepMVA_mvaId",&_daughters_lepMVA_mvaId); - + + myTree->Branch("daughters_pca_x",&_daughters_pca_x); + myTree->Branch("daughters_pca_y",&_daughters_pca_y); + myTree->Branch("daughters_pca_z",&_daughters_pca_z); + myTree->Branch("daughters_pcaRefitPV_x",&_daughters_pcaRefitPV_x); + myTree->Branch("daughters_pcaRefitPV_y",&_daughters_pcaRefitPV_y); + myTree->Branch("daughters_pcaRefitPV_z",&_daughters_pcaRefitPV_z); + myTree->Branch("daughters_pcaGenPV_x",&_daughters_pcaGenPV_x); + myTree->Branch("daughters_pcaGenPV_y",&_daughters_pcaGenPV_y); + myTree->Branch("daughters_pcaGenPV_z",&_daughters_pcaGenPV_z); + myTree->Branch("JetsNumber",&_numberOfJets,"JetsNumber/I"); myTree->Branch("jets_px",&_jets_px); myTree->Branch("jets_py",&_jets_py); @@ -1226,6 +1330,20 @@ void HTauTauNtuplizer::beginJob(){ myTree->Branch("subjets_CSV", &_subjets_CSV); myTree->Branch("subjets_ak8MotherIdx", &_subjets_ak8MotherIdx); + + myTree->Branch("pv_x", &_pv_x); + myTree->Branch("pv_y", &_pv_y); + myTree->Branch("pv_z", &_pv_z); + + myTree->Branch("pvRefit_x", &_pvRefit_x); + myTree->Branch("pvRefit_y", &_pvRefit_y); + myTree->Branch("pvRefit_z", &_pvRefit_z); + + myTree->Branch("pvGen_x", &_pvGen_x); + myTree->Branch("pvGen_y", &_pvGen_y); + myTree->Branch("pvGen_z", &_pvGen_z); + + myTree->Branch("isRefitPV", &_isRefitPV); } Int_t HTauTauNtuplizer::FindCandIndex(const reco::Candidate& cand,Int_t iCand=0){ @@ -1245,8 +1363,10 @@ void HTauTauNtuplizer::analyze(const edm::Event& event, const edm::EventSetup& e { Initialize(); + findPrimaryVertices(event, eSetup); + Handle > vertexs; - //event.getByLabel("offlineSlimmedPrimaryVertices",vertexs); + //event.getByLabel("offlineSlimmedPrimaryVertices",vertex); event.getByToken(theVtxTag,vertexs); //---------------------------------------------------------------------- @@ -1416,7 +1536,7 @@ void HTauTauNtuplizer::analyze(const edm::Event& event, const edm::EventSetup& e FillGenJetInfo(event); // gen jets } //Loop of softleptons and fill them - FillSoftLeptons(daus,event,theFSR,jets); + FillSoftLeptons(daus,event,eSetup,theFSR,jets); //Loop on Jets @@ -1838,7 +1958,9 @@ void HTauTauNtuplizer::FillFatJet(const edm::View* fatjets, const edm: //Fill all leptons (we keep them all for veto purposes -void HTauTauNtuplizer::FillSoftLeptons(const edm::View *daus, const edm::Event& event,bool theFSR, const edm::View *jets){ +void HTauTauNtuplizer::FillSoftLeptons(const edm::View *daus, + const edm::Event& event, const edm::EventSetup& setup, + bool theFSR, const edm::View *jets){ edm::Handle triggerBits; edm::Handle triggerObjects; @@ -1862,8 +1984,12 @@ void HTauTauNtuplizer::FillSoftLeptons(const edm::View *daus, c float rho_miniRelIso = *rhoHandle_miniRelIso; for(edm::View::const_iterator daui = daus->begin(); daui!=daus->end();++daui){ + const reco::Candidate* cand = &(*daui); math::XYZTLorentzVector pfour = cand->p4(); + math::XYZTLorentzVector chargedP4 = cand->p4(); + math::XYZTLorentzVector neutralP4; + TLorentzVector pfourTauUp; TLorentzVector pfourTauDown; @@ -2006,8 +2132,16 @@ void HTauTauNtuplizer::FillSoftLeptons(const edm::View *daus, c float lepMVA_mvaId = -1.; // + GlobalPoint aPVPoint(_pv_x, _pv_y, _pv_z); + GlobalPoint aPVRefitPoint(_pvRefit_x, _pvRefit_y, _pvRefit_z); + GlobalPoint aPVGenPoint(_pvGen_x, _pvGen_y, _pvGen_z); + + TVector3 pcaPV = getPCA(event, setup, cand->bestTrack(), aPVPoint); + TVector3 pcaRefitPV = getPCA(event, setup, cand->bestTrack(), aPVRefitPoint); + TVector3 pcaGenPV; + if(theisMC) pcaGenPV = getPCA(event, setup, cand->bestTrack(), aPVGenPoint); - if(type==ParticleType::MUON){ + if(type==ParticleType::MUON){ muIDflag=userdatahelpers::getUserInt(cand,"muonID"); discr = (float) muIDflag; // not really needed, will use the muonID branch in ntuples... if(userdatahelpers::getUserFloat(cand,"isPFMuon"))typeOfMuon |= 1 << 0; @@ -2103,6 +2237,19 @@ void HTauTauNtuplizer::FillSoftLeptons(const edm::View *daus, c leadChargedParticlePt = userdatahelpers::getUserFloat (cand, "leadChargedParticlePt"); trackRefPt = userdatahelpers::getUserFloat (cand, "trackRefPt"); + const pat::Tau *taon = dynamic_cast(cand); + if(taon){ + pcaPV = getPCA(event, setup, taon->leadChargedHadrCand()->bestTrack(), aPVPoint); + pcaRefitPV = getPCA(event, setup, taon->leadChargedHadrCand()->bestTrack(), aPVRefitPoint); + if(theisMC) pcaGenPV = getPCA(event, setup, taon->leadChargedHadrCand()->bestTrack(), aPVGenPoint); + + reco::CandidatePtrVector chCands = taon->signalChargedHadrCands(); + reco::CandidatePtrVector neCands = taon->signalGammaCands(); + chargedP4 = math::XYZTLorentzVector(); + neutralP4 = math::XYZTLorentzVector(); + for(reco::CandidatePtrVector::const_iterator id=chCands.begin();id!=chCands.end(); ++id) chargedP4 += (*id)->p4(); + for(reco::CandidatePtrVector::const_iterator id=neCands.begin();id!=neCands.end(); ++id) neutralP4 += (*id)->p4(); + } } _discriminator.push_back(discr); _daughters_typeOfMuon.push_back(typeOfMuon); @@ -2170,6 +2317,28 @@ void HTauTauNtuplizer::FillSoftLeptons(const edm::View *daus, c _daughters_jetBTagCSV.push_back(jetBTagCSV); _daughters_lepMVA_mvaId.push_back(lepMVA_mvaId); + _daughters_pca_x.push_back(pcaPV.X()); + _daughters_pca_y.push_back(pcaPV.Y()); + _daughters_pca_z.push_back(pcaPV.Z()); + + _daughters_pcaRefitPV_x.push_back(pcaRefitPV.X()); + _daughters_pcaRefitPV_y.push_back(pcaRefitPV.Y()); + _daughters_pcaRefitPV_z.push_back(pcaRefitPV.Z()); + + _daughters_pcaGenPV_x.push_back(pcaGenPV.X()); + _daughters_pcaGenPV_y.push_back(pcaGenPV.Y()); + _daughters_pcaGenPV_z.push_back(pcaGenPV.Z()); + + _daughters_charged_px.push_back(chargedP4.X()); + _daughters_charged_py.push_back(chargedP4.Y()); + _daughters_charged_pz.push_back(chargedP4.Z()); + _daughters_charged_e.push_back(chargedP4.T()); + + _daughters_neutral_px.push_back(neutralP4.X()); + _daughters_neutral_py.push_back(neutralP4.Y()); + _daughters_neutral_pz.push_back(neutralP4.Z()); + _daughters_neutral_e.push_back(neutralP4.T()); + //TRIGGER MATCHING Long64_t LFtriggerbit=0,L3triggerbit=0,filterFired=0; Long64_t trgMatched = 0; @@ -2356,6 +2525,7 @@ void HTauTauNtuplizer::FillGenInfo(const edm::Event& event) //event.getByLabel ("genInfo", candHandle); event.getByToken (theGenericTag, candHandle); const edm::View* gens = candHandle.product(); + for(edm::View::const_iterator igen = gens->begin(); igen!=gens->end(); ++igen) { // fill gen particle branches @@ -2377,6 +2547,8 @@ void HTauTauNtuplizer::FillGenInfo(const edm::Event& event) int TopDecayMode = -1; int WDecayMode = -1; int TauGenDecayMode = -1; + int TauGenDetailedDecayMode = -1; + TVector3 pca(99,99,99); if (igen->hasUserInt("HMothIndex")) HMIndex = igen->userInt("HMothIndex"); if (igen->hasUserInt("MSSMHMothIndex")) MSSMHMIndex = igen->userInt("MSSMHMothIndex"); @@ -2389,6 +2561,9 @@ void HTauTauNtuplizer::FillGenInfo(const edm::Event& event) if (igen->hasUserInt("TopDecayMode")) TopDecayMode = igen->userInt("TopDecayMode"); if (igen->hasUserInt("WDecayMode")) WDecayMode = igen->userInt("WDecayMode"); if (igen->hasUserInt("tauGenDecayMode")) TauGenDecayMode = igen->userInt("tauGenDecayMode"); + if (igen->hasUserInt("tauGenDetailedDecayMode")) TauGenDetailedDecayMode = igen->userInt("tauGenDetailedDecayMode"); + if (igen->hasUserFloat("pca_x")) pca = TVector3(igen->userFloat("pca_x"),igen->userFloat("pca_y"),igen->userFloat("pca_z")); + _genpart_HMothInd.push_back(HMIndex); _genpart_MSSMHMothInd.push_back(MSSMHMIndex); @@ -2401,11 +2576,21 @@ void HTauTauNtuplizer::FillGenInfo(const edm::Event& event) _genpart_TopDecayMode.push_back(TopDecayMode); _genpart_WDecayMode.push_back(WDecayMode); _genpart_TauGenDecayMode.push_back(TauGenDecayMode); + _genpart_TauGenDetailedDecayMode.push_back(TauGenDetailedDecayMode); + _genpart_pca_x.push_back(pca.X()); + _genpart_pca_y.push_back(pca.Y()); + _genpart_pca_z.push_back(pca.Z()); //const pat::GenericParticle* genClone = &(*igen); //int flags = CreateFlagsWord (genClone); int flags = igen -> userInt ("generalGenFlags"); _genpart_flags.push_back(flags); + + if(igen->hasUserInt("HZDecayMode") || igen->hasUserInt("WDecayMode") || igen->hasUserInt("TopDecayMode")){ + _pvGen_x = igen->vx(); + _pvGen_y = igen->vy(); + _pvGen_z = igen->vz(); + } } } @@ -2719,6 +2904,122 @@ float HTauTauNtuplizer::ComputeMT (math::XYZTLorentzVector visP4, float METx, fl return TMath::Sqrt (scalSum*scalSum - vecSumPt*vecSumPt); } +///////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////// +bool HTauTauNtuplizer::refitPV(const edm::Event & iEvent, const edm::EventSetup & iSetup){ + + edm::ESHandle transTrackBuilder; + iSetup.get().get("TransientTrackBuilder",transTrackBuilder); + + edm::Handle >pfCandHandle; + iEvent.getByToken(thePFCandTag,pfCandHandle); + const edm::View* cands = pfCandHandle.product(); + + Handle > vertices; + iEvent.getByToken(theVtxTag,vertices); + + edm::Handle beamSpot; + iEvent.getByToken(beamSpotTag, beamSpot); + + TransientVertex transVtx; + + //Get tracks associated wiht pfPV + reco::TrackCollection pvTracks; + TLorentzVector aTrack; + for(size_t i=0; isize(); ++i){ + if((*cands)[i].charge()==0 || (*cands)[i].vertexRef().isNull()) continue; + if(!(*cands)[i].bestTrack()) continue; + + unsigned int key = (*cands)[i].vertexRef().key(); + int quality = (*cands)[i].pvAssociationQuality(); + + if(key!=0 || + (quality!=pat::PackedCandidate::UsedInFitTight + && quality!=pat::PackedCandidate::UsedInFitLoose)) continue; + + pvTracks.push_back(*((*cands)[i].bestTrack())); + } + ///Built transient tracks from tracks. + std::vector transTracks; + for(auto iter: pvTracks) transTracks.push_back(transTrackBuilder->build(iter)); + + bool fitOk = false; + if(transTracks.size() >= 2 ) { + AdaptiveVertexFitter avf; + avf.setWeightThreshold(0.001); + try { + transVtx = avf.vertex(transTracks, *beamSpot); + fitOk = true; + } catch (...) { + fitOk = false; + std::cout<<"Vtx fit failed!"< > vertices; + iEvent.getByToken(theVtxTag,vertices); + + if(vertices->size()==0) return false; //at least one vertex + + _pv_x = (*vertices)[0].x(); + _pv_y = (*vertices)[0].y(); + _pv_z = (*vertices)[0].z(); + + _isRefitPV = refitPV(iEvent, iSetup); + + return true; +} +////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////// +TVector3 HTauTauNtuplizer::getPCA(const edm::Event & iEvent, const edm::EventSetup & iSetup, + const reco::Track *aTrack, + const GlobalPoint & aPoint){ + TVector3 aPCA; + if(!doPVRefit || !aTrack || _npv==0 || aTrack->pt()<2) return aPCA; + + edm::ESHandle transTrackBuilder; + iSetup.get().get("TransientTrackBuilder",transTrackBuilder); + if(!transTrackBuilder.isValid()){ + std::cout<<"Problem with TransientTrackBuilder"<build(aTrack); + + AnalyticalImpactPointExtrapolator extrapolator(transTrk.field()); + GlobalPoint pos = extrapolator.extrapolate(transTrk.impactPointState(),aPoint).globalPosition(); + + aPCA.SetX(pos.x() - aPoint.x()); + aPCA.SetY(pos.y() - aPoint.y()); + aPCA.SetZ(pos.z() - aPoint.z()); + + return aPCA; +} +///////////////////////////////////////////////////////////////// +///////////////////////////////////////////////////////////////// + + // // ------------ method fills 'descriptions' with the allowed parameters for the module ------------ // void HTauTauNtuplizer::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { // //The following says we do not know what parameters are allowed so do no validation diff --git a/NtupleProducer/plugins/MuFiller.cc b/NtupleProducer/plugins/MuFiller.cc index d16dad26..ae586744 100644 --- a/NtupleProducer/plugins/MuFiller.cc +++ b/NtupleProducer/plugins/MuFiller.cc @@ -20,6 +20,13 @@ #include #include "DataFormats/VertexReco/interface/Vertex.h" +#include "TrackingTools/TransientTrack/interface/TransientTrackBuilder.h" +#include "TrackingTools/Records/interface/TransientTrackRecord.h" +#include "RecoVertex/VertexPrimitives/interface/TransientVertex.h" +#include "RecoVertex/AdaptiveVertexFit/interface/AdaptiveVertexFitter.h" +#include "TrackingTools/GeomPropagators/interface/AnalyticalTrajectoryExtrapolatorToLine.h" +#include "TrackingTools/GeomPropagators/interface/AnalyticalImpactPointExtrapolator.h" + #include #include #include @@ -125,6 +132,7 @@ MuFiller::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) float dxy_innerTrack = 999.; float dz_innerTrack = 999.; const Vertex* vertex = 0; + if (vertexs->size()>0) { vertex = &(vertexs->front()); dxy = (l.muonBestTrack()->dxy(vertex->position())); @@ -261,7 +269,6 @@ MuFiller::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) iEvent.put(result); } - #include DEFINE_FWK_MODULE(MuFiller); diff --git a/NtupleProducer/plugins/TauFiller.cc b/NtupleProducer/plugins/TauFiller.cc index a62e7583..e5492287 100644 --- a/NtupleProducer/plugins/TauFiller.cc +++ b/NtupleProducer/plugins/TauFiller.cc @@ -161,8 +161,6 @@ TauFiller::TauFiller(const edm::ParameterSet& iConfig) : "chargedIsoPtSum", "neutralIsoPtSum", "puCorrPtSum", - "leadChargedParticlePt", - "trackRefPt" }; @@ -359,7 +357,7 @@ TauFiller::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) l.addUserInt("numNeutralHadronsIsoCone",numNeutralHadronsIsoCone); l.addUserInt("numPhotonsIsoCone",numPhotonsIsoCone); l.addUserInt("numParticlesIsoCone",numParticlesIsoCone); - l.addUserFloat("leadChargedParticlePt",leadChargedParticlePt); + l.addUserFloat("leadChargedParticlePt",leadChargedParticlePt); l.addUserFloat("trackRefPt",trackRefPt); // fill all userfloats diff --git a/NtupleProducer/python/HiggsTauTauProducer.py b/NtupleProducer/python/HiggsTauTauProducer.py index 19dfbb89..e07c8e03 100644 --- a/NtupleProducer/python/HiggsTauTauProducer.py +++ b/NtupleProducer/python/HiggsTauTauProducer.py @@ -719,7 +719,8 @@ srcPFMETCov = cms.InputTag("METSignificance", "METCovariance"), srcPFMETSignificance = cms.InputTag("METSignificance", "METSignificance"), l1extraIsoTau = cms.InputTag("l1extraParticles", "IsoTau"), - HT = cms.InputTag("externalLHEProducer") + HT = cms.InputTag("externalLHEProducer"), + beamSpot = cms.InputTag("offlineBeamSpot") ) if USE_NOHFMET: process.HTauTauTree.metCollection = cms.InputTag("slimmedMETsNoHF") diff --git a/NtupleProducer/python/HiggsTauTauProducer_80X.py b/NtupleProducer/python/HiggsTauTauProducer_80X.py index 0a4bf79e..759b0c19 100644 --- a/NtupleProducer/python/HiggsTauTauProducer_80X.py +++ b/NtupleProducer/python/HiggsTauTauProducer_80X.py @@ -717,7 +717,8 @@ srcPFMETCov = cms.InputTag("METSignificance", "METCovariance"), srcPFMETSignificance = cms.InputTag("METSignificance", "METSignificance"), l1extraIsoTau = cms.InputTag("l1extraParticles", "IsoTau"), - HT = cms.InputTag("externalLHEProducer") + HT = cms.InputTag("externalLHEProducer"), + beamSpot = cms.InputTag("offlineBeamSpot") ) if USE_NOHFMET: process.HTauTauTree.metCollection = cms.InputTag("slimmedMETsNoHF") diff --git a/NtupleProducer/src/GenHelper.cc b/NtupleProducer/src/GenHelper.cc index 30426a4d..7f8485c1 100644 --- a/NtupleProducer/src/GenHelper.cc +++ b/NtupleProducer/src/GenHelper.cc @@ -11,6 +11,8 @@ #include "LLRHiggsTauTau/NtupleProducer/interface/GenHelper.h" #include +#include "DataFormats/TauReco/interface/PFTauDecayMode.h" + bool genhelper::IsLastCopy (const reco::GenParticle& part) { bool isLast = true; @@ -264,6 +266,29 @@ reco::GenParticle genhelper::GetTauHad (const reco::Candidate* part) return TauH; } +reco::GenParticle genhelper::GetTauHadNeutrals (const reco::Candidate* part) +{ + if (abs(part->pdgId()) != 15) + { + reco::GenParticle fakeTauH = reco::GenParticle (0, reco::Candidate::LorentzVector(0.,0.,0.,0.), reco::Candidate::Point (0.,0.,0.), -999999, 0, true); + std::cout << "Warning: building had tau from a particle with pdgId != 15 --> dummy entry returned" << std::endl; + return fakeTauH; + } + + reco::Candidate::LorentzVector p4Had (0,0,0,0); + for (unsigned int iDau = 0; iDau < part->numberOfDaughters(); iDau++) + { + const reco::Candidate * Dau = part->daughter( iDau ); + int dauId = abs(Dau->pdgId()); + if (dauId != 12 && dauId != 14 && dauId != 16 && Dau->charge()==0) // no neutrinos + p4Had += Dau->p4(); + } + + int sign = part->pdgId() / abs(part->pdgId()); + reco::GenParticle TauH = reco::GenParticle (part->charge(), p4Had, part->vertex(), sign*77715, part->status(), true); + return TauH; +} + const reco::Candidate* genhelper::IsFromID (const reco::Candidate* part, int targetPDGId) { if (abs(part->pdgId()) == targetPDGId){ @@ -290,3 +315,194 @@ int genhelper::GetIndexInOutput (const reco::Candidate* part, std::vector allNus; + allNus.insert(12); + allNus.insert(14); + allNus.insert(16); + //allNus.insert(18); + reco::GenParticleRefVector tmp; + for(IGR idr=products.begin(); idr !=products.end(); ++idr) + if(allNus.find(std::abs((*idr)->pdgId()))==allNus.end()) + tmp.push_back((*idr)); + products.swap(tmp); + } +} +/////////////////////////////////////////////////////// +/////////////////////////////////////////////////////// +//copy of function from PhysicsTools/HepMCCandAlgos/interface/GenParticlesHelper.h" which allows ignore status +void genhelper::FindDescendents(const reco::GenParticle& base, + reco::GenParticleRefVector& descendents, + int status, int pdgId, + bool skipPhotonsPi0AndFSR) { + + //one form status or pdgId has to be specifed! + if(status<0 && pdgId==0) return; + + const reco::GenParticleRefVector& daughterRefs = base.daughterRefVector(); + + for(IGR idr = daughterRefs.begin(); idr != daughterRefs.end(); ++idr ) { + + ///Skip leptons from pi0 decays + if(skipPhotonsPi0AndFSR && (*idr)->mother(0) && (abs((*idr)->mother(0)->pdgId())==22 || abs((*idr)->mother(0)->pdgId())==111)) continue; + ///Skip electrons from FSR from muons + if(skipPhotonsPi0AndFSR && (*idr)->mother(0) && abs((*idr)->mother(0)->pdgId())==13 && abs((*idr)->pdgId())==11) continue; + ///Skip muons and electrons from FSR + if(skipPhotonsPi0AndFSR && (*idr)->mother(0) && abs((*idr)->mother(0)->pdgId())==13 && abs((*idr)->pdgId())==11) continue; + if(skipPhotonsPi0AndFSR && (*idr)->mother(0) && abs((*idr)->mother(0)->pdgId())==11 && abs((*idr)->pdgId())==13) continue; + + + if( (status<0 || (*idr)->status() == status ) && + (!pdgId || std::abs((*idr)->pdgId()) == std::abs(pdgId) )) { + descendents.push_back(*idr); + } + else FindDescendents( *(*idr), descendents, status, pdgId, skipPhotonsPi0AndFSR); + } +} +/////////////////////////////////////////////////////// +/////////////////////////////////////////////////////// +const reco::GenParticleRef genhelper::GetLeadChParticle(const reco::GenParticleRefVector& products){ + + float maxPt=0; + reco::GenParticleRef part; + std::set charged; + charged.insert(211);//pi + charged.insert(321);//K + charged.insert(11);//e + charged.insert(13);//mu + + for(IGR idr = products.begin(); idr != products.end(); ++idr ){ + if( (*idr)->pt() > maxPt && + //charged.find( std::abs( (*idr)->pdgId() ) )!=charged.end() //MB: Logix used in pure Pythia code when charge not defined + std::abs( (*idr)->charge() )>0.001 //MB: GenParts have defined charge + ){ + maxPt = (*idr)->pt(); + part = (*idr); + } + } + return part; + } +/////////////////////////////////////////////////////// +/////////////////////////////////////////////////////// +int genhelper::getDetailedTauDecayMode(const reco::GenParticleRefVector& products){ + + int tauDecayMode = reco::PFTauDecayMode::tauDecayOther; + + int numElectrons = 0; + int numMuons = 0; + int numChargedPions = 0; + int numNeutralPions = 0; + int numPhotons = 0; + int numNeutrinos = 0; + int numOtherParticles = 0; + + for(IGR idr = products.begin(); idr != products.end(); ++idr ) { + int pdg_id = std::abs((*idr)->pdgId()); + if(pdg_id == 11) numElectrons++; + else if(pdg_id == 13) numMuons++; + else if(pdg_id == 211 || pdg_id == 321 ) numChargedPions++; //Count both pi+ and K+ + else if(pdg_id == 111 || pdg_id == 130 || pdg_id == 310 ) numNeutralPions++; //Count both pi0 and K0_L/S + else if(pdg_id == 12 || + pdg_id == 14 || + pdg_id == 16) { + numNeutrinos++; + } + else if(pdg_id == 22) numPhotons++; + else { + numOtherParticles++; + } + } + if(numElectrons>1){//sometimes there are gamma->ee conversions + numPhotons += numElectrons/2; + numElectrons -= 2*(numElectrons/2); + } + + if( numOtherParticles == 0 ){ + if( numElectrons == 1 ){ + //--- tau decays into electrons + tauDecayMode = reco::PFTauDecayMode::tauDecaysElectron; + } else if( numMuons == 1 ){ + //--- tau decays into muons + tauDecayMode = reco::PFTauDecayMode::tauDecayMuon; + } else { + //--- hadronic tau decays + switch ( numChargedPions ){ + case 1 : + if( numNeutralPions != 0 ){ + tauDecayMode = reco::PFTauDecayMode::tauDecayOther; + break; + } + switch ( numPhotons ){ + case 0: + tauDecayMode = reco::PFTauDecayMode::tauDecay1ChargedPion0PiZero; + break; + case 2: + tauDecayMode = reco::PFTauDecayMode::tauDecay1ChargedPion1PiZero; + break; + case 4: + tauDecayMode = reco::PFTauDecayMode::tauDecay1ChargedPion2PiZero; + break; + case 6: + tauDecayMode = reco::PFTauDecayMode::tauDecay1ChargedPion3PiZero; + break; + case 8: + tauDecayMode = reco::PFTauDecayMode::tauDecay1ChargedPion4PiZero; + break; + default: + tauDecayMode = reco::PFTauDecayMode::tauDecayOther; + break; + } + break; + case 3 : + if( numNeutralPions != 0 ){ + tauDecayMode = reco::PFTauDecayMode::tauDecayOther; + break; + } + switch ( numPhotons ){ + case 0 : + tauDecayMode = reco::PFTauDecayMode::tauDecay3ChargedPion0PiZero; + break; + case 2 : + tauDecayMode = reco::PFTauDecayMode::tauDecay3ChargedPion1PiZero; + break; + case 4 : + tauDecayMode = reco::PFTauDecayMode::tauDecay3ChargedPion2PiZero; + break; + case 6 : + tauDecayMode = reco::PFTauDecayMode::tauDecay3ChargedPion3PiZero; + break; + case 8 : + tauDecayMode = reco::PFTauDecayMode::tauDecay3ChargedPion4PiZero; + break; + default: + tauDecayMode = reco::PFTauDecayMode::tauDecayOther; + break; + } + break; + } + } + } + return tauDecayMode; +} +/////////////////////////////////////////////////////// +///////////////////////////////////////////////////////