From 9e17afd43f4ede1049cee672b88a906814e15f81 Mon Sep 17 00:00:00 2001 From: Artur Kalinowski Date: Thu, 14 Jul 2016 16:44:27 +0200 Subject: [PATCH 01/10] * add variables necessary for the CP studies * fix bug in TauFiller where leadChargedParticlePt and trackRefPt were added twice. This made the trackRefPt to contain a wrong (~random) value. --- NtupleProducer/interface/GenHelper.h | 12 +- NtupleProducer/plugins/BuildFile.xml | 1 + NtupleProducer/plugins/GenFiller.cc | 22 ++ NtupleProducer/plugins/HTauTauNtuplizer.cc | 301 ++++++++++++++++++- NtupleProducer/plugins/MuFiller.cc | 38 +++ NtupleProducer/plugins/TauFiller.cc | 4 +- NtupleProducer/python/HiggsTauTauProducer.py | 3 +- NtupleProducer/src/GenHelper.cc | 92 ++++++ 8 files changed, 461 insertions(+), 12 deletions(-) diff --git a/NtupleProducer/interface/GenHelper.h b/NtupleProducer/interface/GenHelper.h index 35c9cf5b..f98a0f63 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 @@ -52,6 +55,13 @@ namespace genhelper { reco::GenParticle GetTauHad (const reco::Candidate* part); // build 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); + +} #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..61225722 100644 --- a/NtupleProducer/plugins/GenFiller.cc +++ b/NtupleProducer/plugins/GenFiller.cc @@ -213,6 +213,28 @@ 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); + reco::GenParticleRef leadChParticleRef = genhelper::GetLeadChParticle(tauDaughters); + + std::cout<<"lead ch part pdg: "<pdgId()<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); } diff --git a/NtupleProducer/plugins/HTauTauNtuplizer.cc b/NtupleProducer/plugins/HTauTauNtuplizer.cc index dd5dca91..3ef3cc2c 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 { @@ -146,7 +155,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 +166,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 @@ -203,6 +217,7 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { edm::EDGetTokenT theTotTag; edm::EDGetTokenT thePassTag; edm::EDGetTokenT theLHEPTag; + edm::EDGetTokenT beamSpotTag; //flags @@ -236,6 +251,10 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { Float_t _PUReweight; Float_t _rho; Int_t _nup; + + Float_t _pv_x, _pv_y, _pv_z; + Float_t _pvGen_x, _pvGen_y, _pvGen_z; + Float_t _pvRefit_x, _pvRefit_y, _pvRefit_z; // pairs //std::vector _mothers; @@ -252,6 +271,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; @@ -278,6 +308,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; @@ -473,6 +508,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; @@ -561,7 +609,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"))) { theFileName = pset.getUntrackedParameter("fileName"); @@ -644,6 +693,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(); @@ -719,6 +779,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(); @@ -732,6 +804,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(); @@ -942,6 +1019,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); @@ -977,6 +1064,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); @@ -1132,7 +1222,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); @@ -1188,6 +1288,19 @@ void HTauTauNtuplizer::beginJob(){ myTree->Branch("subjets_CSV", &_subjets_CSV); myTree->Branch("subjets_ak8MotherIdx", &_subjets_ak8MotherIdx); + + myTree->Branch("px_x", &_pv_x); + myTree->Branch("px_y", &_pv_y); + myTree->Branch("px_z", &_pv_z); + + myTree->Branch("pxRefit_x", &_pvRefit_x); + myTree->Branch("pxRefit_y", &_pvRefit_y); + myTree->Branch("pxRefit_z", &_pvRefit_z); + + myTree->Branch("pxGen_x", &_pvGen_x); + myTree->Branch("pxGen_y", &_pvGen_y); + myTree->Branch("pxGen_z", &_pvGen_z); + } Int_t HTauTauNtuplizer::FindCandIndex(const reco::Candidate& cand,Int_t iCand=0){ @@ -1207,8 +1320,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); //---------------------------------------------------------------------- @@ -1365,7 +1480,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 @@ -1785,7 +1900,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; @@ -1809,8 +1926,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; @@ -1948,8 +2069,15 @@ 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 = 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; @@ -2044,6 +2172,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); + 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); @@ -2110,6 +2251,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 int LFtriggerbit=0,L3triggerbit=0,filterFired=0; int triggertypeIsGood = 0; @@ -2261,6 +2424,9 @@ void HTauTauNtuplizer::FillGenInfo(const edm::Event& event) //event.getByLabel ("genInfo", candHandle); event.getByToken (theGenericTag, candHandle); const edm::View* gens = candHandle.product(); + + std::vector tauDaughters; + for(edm::View::const_iterator igen = gens->begin(); igen!=gens->end(); ++igen) { // fill gen particle branches @@ -2282,6 +2448,7 @@ void HTauTauNtuplizer::FillGenInfo(const edm::Event& event) int TopDecayMode = -1; int WDecayMode = -1; int TauGenDecayMode = -1; + TVector3 pca(99,99,99); if (igen->hasUserInt("HMothIndex")) HMIndex = igen->userInt("HMothIndex"); if (igen->hasUserInt("MSSMHMothIndex")) MSSMHMIndex = igen->userInt("MSSMHMothIndex"); @@ -2294,6 +2461,8 @@ 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->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); @@ -2306,11 +2475,20 @@ void HTauTauNtuplizer::FillGenInfo(const edm::Event& event) _genpart_TopDecayMode.push_back(TopDecayMode); _genpart_WDecayMode.push_back(WDecayMode); _genpart_TauGenDecayMode.push_back(TauGenDecayMode); + _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(); + } } } @@ -2624,6 +2802,115 @@ 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(); + + 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(!aTrack) return aPCA; + + edm::ESHandle transTrackBuilder; + iSetup.get().get("TransientTrackBuilder",transTrackBuilder); + reco::TransientTrack transTrk=transTrackBuilder->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..0404d61c 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 @@ -46,6 +53,10 @@ virtual void produce(edm::Event&, const edm::EventSetup&); virtual void endJob(){}; + TVector3 getPCA(const edm::Event & iEvent, const edm::EventSetup & iSetup, + const reco::Track *aTrack, + const GlobalPoint & aPoint); + //const edm::InputTag theCandidateTag; //edm::EDGetTokenT > theCandidateTag; edm::EDGetTokenT theCandidateTag; @@ -125,12 +136,18 @@ MuFiller::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) float dxy_innerTrack = 999.; float dz_innerTrack = 999.; const Vertex* vertex = 0; + + TVector3 pcaAOD, pcaRefit; + if (vertexs->size()>0) { vertex = &(vertexs->front()); dxy = (l.muonBestTrack()->dxy(vertex->position())); dz = (l.muonBestTrack()->dz(vertex->position())); dxy_innerTrack = (l.innerTrack()->dxy(vertex->position())); dz_innerTrack = (l.innerTrack()->dz(vertex->position())); + + GlobalPoint aPoint(vertex->position().x(), vertex->position().y(), vertex->position().z()); + pcaAOD = getPCA(iEvent, iSetup, l.muonBestTrack().get(), aPoint); } float rel_error_trackpt = l.muonBestTrack()->ptError()/l.muonBestTrack()->pt(); @@ -260,7 +277,28 @@ MuFiller::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) } iEvent.put(result); } +////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////// +TVector3 MuFiller::getPCA(const edm::Event & iEvent, const edm::EventSetup & iSetup, + const reco::Track *aTrack, + const GlobalPoint & aPoint){ + + edm::ESHandle transTrackBuilder; + iSetup.get().get("TransientTrackBuilder",transTrackBuilder); + reco::TransientTrack transTrk=transTrackBuilder->build(aTrack); + + AnalyticalImpactPointExtrapolator extrapolator(transTrk.field()); + GlobalPoint pos = extrapolator.extrapolate(transTrk.impactPointState(),aPoint).globalPosition(); + TVector3 aPCA; + aPCA.SetX(pos.x() - aPoint.x()); + aPCA.SetY(pos.y() - aPoint.y()); + aPCA.SetZ(pos.z() - aPoint.z()); + + return aPCA; +} +////////////////////////////////////////////////////////////////// +////////////////////////////////////////////////////////////////// #include DEFINE_FWK_MODULE(MuFiller); diff --git a/NtupleProducer/plugins/TauFiller.cc b/NtupleProducer/plugins/TauFiller.cc index 3aca677a..5d7600cf 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 b7038d3d..f263eddc 100644 --- a/NtupleProducer/python/HiggsTauTauProducer.py +++ b/NtupleProducer/python/HiggsTauTauProducer.py @@ -710,7 +710,8 @@ metFilters = cms.InputTag ("TriggerResults","",METfiltersProcess), PUPPImetCollection = cms.InputTag("slimmedMETsPuppi"), 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..79e62278 100644 --- a/NtupleProducer/src/GenHelper.cc +++ b/NtupleProducer/src/GenHelper.cc @@ -290,3 +290,95 @@ 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; + } +/////////////////////////////////////////////////////// +/////////////////////////////////////////////////////// From 07a9c454f93bb601e5a7e04099c2ce5b7d1f5db5 Mon Sep 17 00:00:00 2001 From: Artur Kalinowski Date: Thu, 14 Jul 2016 16:57:59 +0200 Subject: [PATCH 02/10] * remove test printout. --- NtupleProducer/plugins/GenFiller.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/NtupleProducer/plugins/GenFiller.cc b/NtupleProducer/plugins/GenFiller.cc index 61225722..d18e0df7 100644 --- a/NtupleProducer/plugins/GenFiller.cc +++ b/NtupleProducer/plugins/GenFiller.cc @@ -222,8 +222,6 @@ void GenFiller::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) genhelper::GetTausDaughters(*genPClone,tauDaughters,true,false); reco::GenParticleRef leadChParticleRef = genhelper::GetLeadChParticle(tauDaughters); - std::cout<<"lead ch part pdg: "<pdgId()<px(), leadChParticleRef->py(), leadChParticleRef->pz(), From aba1cda052ef6a4673545133b2d27a97c7f42bff Mon Sep 17 00:00:00 2001 From: Artur Kalinowski Date: Fri, 15 Jul 2016 12:24:17 +0200 Subject: [PATCH 03/10] * save neutral component of the hadronic generator level tau --- NtupleProducer/interface/GenHelper.h | 2 ++ NtupleProducer/plugins/GenFiller.cc | 45 ++++++++++++++++++++++------ NtupleProducer/src/GenHelper.cc | 23 ++++++++++++++ 3 files changed, 61 insertions(+), 9 deletions(-) diff --git a/NtupleProducer/interface/GenHelper.h b/NtupleProducer/interface/GenHelper.h index f98a0f63..cd17e3eb 100644 --- a/NtupleProducer/interface/GenHelper.h +++ b/NtupleProducer/interface/GenHelper.h @@ -53,6 +53,8 @@ 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); diff --git a/NtupleProducer/plugins/GenFiller.cc b/NtupleProducer/plugins/GenFiller.cc index d18e0df7..a771bb64 100644 --- a/NtupleProducer/plugins/GenFiller.cc +++ b/NtupleProducer/plugins/GenFiller.cc @@ -242,24 +242,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/src/GenHelper.cc b/NtupleProducer/src/GenHelper.cc index 79e62278..09aa458a 100644 --- a/NtupleProducer/src/GenHelper.cc +++ b/NtupleProducer/src/GenHelper.cc @@ -264,6 +264,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){ From 5e511ec1b0b132b1a4834422b00dfacb2d6c13dc Mon Sep 17 00:00:00 2001 From: Artur Kalinowski Date: Fri, 15 Jul 2016 12:51:14 +0200 Subject: [PATCH 04/10] * fix PV components branch names --- NtupleProducer/plugins/HTauTauNtuplizer.cc | 18 +++++++++--------- 1 file changed, 9 insertions(+), 9 deletions(-) diff --git a/NtupleProducer/plugins/HTauTauNtuplizer.cc b/NtupleProducer/plugins/HTauTauNtuplizer.cc index 3ef3cc2c..ece7f847 100644 --- a/NtupleProducer/plugins/HTauTauNtuplizer.cc +++ b/NtupleProducer/plugins/HTauTauNtuplizer.cc @@ -1289,17 +1289,17 @@ void HTauTauNtuplizer::beginJob(){ myTree->Branch("subjets_ak8MotherIdx", &_subjets_ak8MotherIdx); - myTree->Branch("px_x", &_pv_x); - myTree->Branch("px_y", &_pv_y); - myTree->Branch("px_z", &_pv_z); + myTree->Branch("pv_x", &_pv_x); + myTree->Branch("pv_y", &_pv_y); + myTree->Branch("pv_z", &_pv_z); - myTree->Branch("pxRefit_x", &_pvRefit_x); - myTree->Branch("pxRefit_y", &_pvRefit_y); - myTree->Branch("pxRefit_z", &_pvRefit_z); + myTree->Branch("pvRefit_x", &_pvRefit_x); + myTree->Branch("pvRefit_y", &_pvRefit_y); + myTree->Branch("pvRefit_z", &_pvRefit_z); - myTree->Branch("pxGen_x", &_pvGen_x); - myTree->Branch("pxGen_y", &_pvGen_y); - myTree->Branch("pxGen_z", &_pvGen_z); + myTree->Branch("pvGen_x", &_pvGen_x); + myTree->Branch("pvGen_y", &_pvGen_y); + myTree->Branch("pvGen_z", &_pvGen_z); } From cd17d0a5455bfa608f79d29d8e7b586f87cc0b78 Mon Sep 17 00:00:00 2001 From: Artur Kalinowski Date: Tue, 19 Jul 2016 15:45:27 +0200 Subject: [PATCH 05/10] * fix wrong refitted PV position when refit fails. --- NtupleProducer/interface/GenHelper.h | 1 + NtupleProducer/plugins/GenFiller.cc | 3 + NtupleProducer/plugins/HTauTauNtuplizer.cc | 28 ++++-- NtupleProducer/src/GenHelper.cc | 104 +++++++++++++++++++++ 4 files changed, 128 insertions(+), 8 deletions(-) diff --git a/NtupleProducer/interface/GenHelper.h b/NtupleProducer/interface/GenHelper.h index cd17e3eb..191dd680 100644 --- a/NtupleProducer/interface/GenHelper.h +++ b/NtupleProducer/interface/GenHelper.h @@ -64,6 +64,7 @@ namespace genhelper { 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/GenFiller.cc b/NtupleProducer/plugins/GenFiller.cc index a771bb64..969b682d 100644 --- a/NtupleProducer/plugins/GenFiller.cc +++ b/NtupleProducer/plugins/GenFiller.cc @@ -220,6 +220,9 @@ void GenFiller::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) 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(), diff --git a/NtupleProducer/plugins/HTauTauNtuplizer.cc b/NtupleProducer/plugins/HTauTauNtuplizer.cc index ece7f847..b214b9ec 100644 --- a/NtupleProducer/plugins/HTauTauNtuplizer.cc +++ b/NtupleProducer/plugins/HTauTauNtuplizer.cc @@ -255,6 +255,7 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { Float_t _pv_x, _pv_y, _pv_z; Float_t _pvGen_x, _pvGen_y, _pvGen_z; Float_t _pvRefit_x, _pvRefit_y, _pvRefit_z; + bool _isRefitPV; // pairs //std::vector _mothers; @@ -327,6 +328,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 @@ -823,6 +825,7 @@ void HTauTauNtuplizer::Initialize(){ _genpart_TopDecayMode.clear(); _genpart_WDecayMode.clear(); _genpart_TauGenDecayMode.clear(); + _genpart_TauGenDetailedDecayMode.clear(); _genpart_flags.clear(); _genjet_px.clear(); @@ -1081,6 +1084,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); @@ -1301,6 +1305,7 @@ void HTauTauNtuplizer::beginJob(){ 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){ @@ -2448,6 +2453,7 @@ 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"); @@ -2461,6 +2467,7 @@ 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")); @@ -2475,6 +2482,7 @@ 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()); @@ -2862,12 +2870,12 @@ bool HTauTauNtuplizer::refitPV(const edm::Event & iEvent, const edm::EventSetup _pvRefit_z = (*vertices)[0].z(); } else { - _pvRefit_x = -999; - _pvRefit_y = -999; - _pvRefit_z = -999; + _pvRefit_x = (*vertices)[0].x(); + _pvRefit_y = (*vertices)[0].y(); + _pvRefit_z = (*vertices)[0].z(); } - return true; + return fitOk && transVtx.isValid(); } ////////////////////////////////////////////////////////////////// ////////////////////////////////////////////////////////////////// @@ -2882,7 +2890,7 @@ bool HTauTauNtuplizer::findPrimaryVertices(const edm::Event & iEvent, const edm: _pv_y = (*vertices)[0].y(); _pv_z = (*vertices)[0].z(); - refitPV(iEvent, iSetup); + _isRefitPV = refitPV(iEvent, iSetup); return true; } @@ -2893,12 +2901,16 @@ TVector3 HTauTauNtuplizer::getPCA(const edm::Event & iEvent, const edm::EventSet const GlobalPoint & aPoint){ TVector3 aPCA; if(!aTrack) return aPCA; - + edm::ESHandle transTrackBuilder; - iSetup.get().get("TransientTrackBuilder",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()); diff --git a/NtupleProducer/src/GenHelper.cc b/NtupleProducer/src/GenHelper.cc index 09aa458a..d5eb9293 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; @@ -405,3 +407,105 @@ const reco::GenParticleRef genhelper::GetLeadChParticle(const reco::GenParticleR } /////////////////////////////////////////////////////// /////////////////////////////////////////////////////// +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; +} +/////////////////////////////////////////////////////// +/////////////////////////////////////////////////////// From f7ff253ee1e5306fe7942e2e489b0b9a35b62eb7 Mon Sep 17 00:00:00 2001 From: Artur Kalinowski Date: Fri, 22 Jul 2016 11:13:46 +0200 Subject: [PATCH 06/10] * add ap rotection against problems in calculating PCA for very low pT tracks. (found for electron candidates) Tracks with pt<1 do not have PCA calculated --- NtupleProducer/plugins/HTauTauNtuplizer.cc | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/NtupleProducer/plugins/HTauTauNtuplizer.cc b/NtupleProducer/plugins/HTauTauNtuplizer.cc index b214b9ec..c00bdc82 100644 --- a/NtupleProducer/plugins/HTauTauNtuplizer.cc +++ b/NtupleProducer/plugins/HTauTauNtuplizer.cc @@ -2078,6 +2078,11 @@ void HTauTauNtuplizer::FillSoftLeptons(const edm::View *daus, GlobalPoint aPVRefitPoint(_pvRefit_x, _pvRefit_y, _pvRefit_z); GlobalPoint aPVGenPoint(_pvGen_x, _pvGen_y, _pvGen_z); + std::cout<<"cand->p4().pt(): "<p4().pt()<<" bestTrack: "<bestTrack()<bestTrack(), aPVPoint); TVector3 pcaRefitPV = getPCA(event, setup, cand->bestTrack(), aPVRefitPoint); TVector3 pcaGenPV = getPCA(event, setup, cand->bestTrack(), aPVGenPoint); @@ -2900,7 +2905,9 @@ TVector3 HTauTauNtuplizer::getPCA(const edm::Event & iEvent, const edm::EventSet const reco::Track *aTrack, const GlobalPoint & aPoint){ TVector3 aPCA; - if(!aTrack) return aPCA; + if(!aTrack || _npv==0 || aTrack->pt()<1.0) return aPCA; + + std::cout< transTrackBuilder; iSetup.get().get("TransientTrackBuilder",transTrackBuilder); @@ -2908,11 +2915,21 @@ TVector3 HTauTauNtuplizer::getPCA(const edm::Event & iEvent, const edm::EventSet std::cout<<"Problem with TransientTrackBuilder"<build(aTrack); + + std::cout<<"Here 2"<pt()<<" impactPointState: "< Date: Thu, 28 Jul 2016 08:34:43 +0200 Subject: [PATCH 07/10] * remove debug printout * add a flag enabling calculation of PCA --- NtupleProducer/plugins/HTauTauNtuplizer.cc | 18 ++---------------- 1 file changed, 2 insertions(+), 16 deletions(-) diff --git a/NtupleProducer/plugins/HTauTauNtuplizer.cc b/NtupleProducer/plugins/HTauTauNtuplizer.cc index c00bdc82..de89fb6a 100644 --- a/NtupleProducer/plugins/HTauTauNtuplizer.cc +++ b/NtupleProducer/plugins/HTauTauNtuplizer.cc @@ -119,6 +119,7 @@ bool writeJets = true; // Write jets in the tree. bool writeFatJets = true; bool writeSoftLep = false; + bool doPVRefit = true; bool DEBUG = false; } @@ -2078,11 +2079,6 @@ void HTauTauNtuplizer::FillSoftLeptons(const edm::View *daus, GlobalPoint aPVRefitPoint(_pvRefit_x, _pvRefit_y, _pvRefit_z); GlobalPoint aPVGenPoint(_pvGen_x, _pvGen_y, _pvGen_z); - std::cout<<"cand->p4().pt(): "<p4().pt()<<" bestTrack: "<bestTrack()<bestTrack(), aPVPoint); TVector3 pcaRefitPV = getPCA(event, setup, cand->bestTrack(), aPVRefitPoint); TVector3 pcaGenPV = getPCA(event, setup, cand->bestTrack(), aPVGenPoint); @@ -2905,9 +2901,7 @@ TVector3 HTauTauNtuplizer::getPCA(const edm::Event & iEvent, const edm::EventSet const reco::Track *aTrack, const GlobalPoint & aPoint){ TVector3 aPCA; - if(!aTrack || _npv==0 || aTrack->pt()<1.0) return aPCA; - - std::cout<pt()<2) return aPCA; edm::ESHandle transTrackBuilder; iSetup.get().get("TransientTrackBuilder",transTrackBuilder); @@ -2916,20 +2910,12 @@ TVector3 HTauTauNtuplizer::getPCA(const edm::Event & iEvent, const edm::EventSet return aPCA; } - std::cout<<"Here 1"<build(aTrack); - std::cout<<"Here 2"<pt()<<" impactPointState: "< Date: Fri, 29 Jul 2016 17:06:05 +0200 Subject: [PATCH 08/10] * add procetion agains issued with wrong PV refit. --- NtupleProducer/plugins/HTauTauNtuplizer.cc | 12 +++++++----- 1 file changed, 7 insertions(+), 5 deletions(-) diff --git a/NtupleProducer/plugins/HTauTauNtuplizer.cc b/NtupleProducer/plugins/HTauTauNtuplizer.cc index 8e2a60cf..5e943242 100644 --- a/NtupleProducer/plugins/HTauTauNtuplizer.cc +++ b/NtupleProducer/plugins/HTauTauNtuplizer.cc @@ -2138,7 +2138,8 @@ void HTauTauNtuplizer::FillSoftLeptons(const edm::View *daus, TVector3 pcaPV = getPCA(event, setup, cand->bestTrack(), aPVPoint); TVector3 pcaRefitPV = getPCA(event, setup, cand->bestTrack(), aPVRefitPoint); - TVector3 pcaGenPV = getPCA(event, setup, cand->bestTrack(), aPVGenPoint); + TVector3 pcaGenPV; + if(theisMC) pcaGenPV = getPCA(event, setup, cand->bestTrack(), aPVGenPoint); if(type==ParticleType::MUON){ muIDflag=userdatahelpers::getUserInt(cand,"muonID"); @@ -2240,7 +2241,7 @@ void HTauTauNtuplizer::FillSoftLeptons(const edm::View *daus, if(taon){ pcaPV = getPCA(event, setup, taon->leadChargedHadrCand()->bestTrack(), aPVPoint); pcaRefitPV = getPCA(event, setup, taon->leadChargedHadrCand()->bestTrack(), aPVRefitPoint); - pcaGenPV = getPCA(event, setup, taon->leadChargedHadrCand()->bestTrack(), aPVGenPoint); + if(theisMC) pcaGenPV = getPCA(event, setup, taon->leadChargedHadrCand()->bestTrack(), aPVGenPoint); reco::CandidatePtrVector chCands = taon->signalChargedHadrCands(); reco::CandidatePtrVector neCands = taon->signalGammaCands(); @@ -2954,8 +2955,10 @@ bool HTauTauNtuplizer::refitPV(const edm::Event & iEvent, const edm::EventSetup std::cout<<"Vtx fit failed!"<build(aTrack); - AnalyticalImpactPointExtrapolator extrapolator(transTrk.field()); GlobalPoint pos = extrapolator.extrapolate(transTrk.impactPointState(),aPoint).globalPosition(); From 47ea60db8fe91e4e0a7513b17938635927d8bd07 Mon Sep 17 00:00:00 2001 From: Artur Kalinowski Date: Sat, 30 Jul 2016 07:11:00 +0200 Subject: [PATCH 09/10] * add misisng input tag for the 80X HiggsTauTauProducer * remove debug printout --- NtupleProducer/python/HiggsTauTauProducer_80X.py | 3 ++- NtupleProducer/src/GenHelper.cc | 3 --- 2 files changed, 2 insertions(+), 4 deletions(-) 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 d5eb9293..7f8485c1 100644 --- a/NtupleProducer/src/GenHelper.cc +++ b/NtupleProducer/src/GenHelper.cc @@ -332,9 +332,6 @@ void genhelper::GetTausDaughters(const reco::GenParticle& tau, FindDescendents(tau, products, 1, 0); else{ const reco::GenParticleRefVector& daughterRefs = tau.daughterRefVector(); - - std::cout<<" tau.daughterRefVector().size(): "<< tau.daughterRefVector().size()< Date: Sun, 31 Jul 2016 21:34:42 +0200 Subject: [PATCH 10/10] * initialize vertex position variables to avoid rare crash due events without initila W/Z boson ot t quark (DY m<50) --- NtupleProducer/plugins/HTauTauNtuplizer.cc | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/NtupleProducer/plugins/HTauTauNtuplizer.cc b/NtupleProducer/plugins/HTauTauNtuplizer.cc index 5e943242..85b14eab 100644 --- a/NtupleProducer/plugins/HTauTauNtuplizer.cc +++ b/NtupleProducer/plugins/HTauTauNtuplizer.cc @@ -261,10 +261,10 @@ class HTauTauNtuplizer : public edm::EDAnalyzer { Float_t _rho; Int_t _nup; - Float_t _pv_x, _pv_y, _pv_z; - Float_t _pvGen_x, _pvGen_y, _pvGen_z; - Float_t _pvRefit_x, _pvRefit_y, _pvRefit_z; - bool _isRefitPV; + 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;