From 7c82e21f49dd16b141ae3ac8f21b0b4b16b6f7b0 Mon Sep 17 00:00:00 2001 From: kurt Date: Tue, 31 Oct 2017 22:56:49 +0100 Subject: [PATCH 1/3] adding 13 TeV JECs and ghost flavor assn for pp --- .../interface/HiInclusiveJetAnalyzer.h | 100 ++- .../python/FullJetSequence_XeXe_mc.py | 4 +- .../python/FullJetSequence_cleanedPbPb.py | 16 +- .../python/FullJetSequence_nominalPP.py | 30 +- .../JetAnalysis/python/bTaggers_cff.py | 72 +- .../python/jets/makeJetSequences.sh | 32 +- .../jets/templateSequence_bTag_cff.py.txt | 36 +- .../JetAnalysis/src/HiInclusiveJetAnalyzer.cc | 622 ++++++++++++++++-- .../test/runForestAOD_XeXe_DATA_92X.py | 17 +- .../test/runForestAOD_XeXe_ppMC_92X.py | 26 +- 10 files changed, 819 insertions(+), 136 deletions(-) diff --git a/HeavyIonsAnalysis/JetAnalysis/interface/HiInclusiveJetAnalyzer.h b/HeavyIonsAnalysis/JetAnalysis/interface/HiInclusiveJetAnalyzer.h index 0c838e533651b..4a1f9590c3ba6 100644 --- a/HeavyIonsAnalysis/JetAnalysis/interface/HiInclusiveJetAnalyzer.h +++ b/HeavyIonsAnalysis/JetAnalysis/interface/HiInclusiveJetAnalyzer.h @@ -33,6 +33,16 @@ #include "DataFormats/L1GlobalTrigger/interface/L1GlobalTriggerObjectMap.h" #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" +#include "RecoBTag/SecondaryVertex/interface/TrackKinematics.h" +#include "DataFormats/Common/interface/View.h" + +#include "DataFormats/PatCandidates/interface/PackedCandidate.h" + +#include "SimDataFormats/JetMatching/interface/JetFlavourInfo.h" +#include "SimDataFormats/JetMatching/interface/JetFlavourInfoMatching.h" +#include "RecoBTag/SecondaryVertex/interface/V0Filter.h" +#include "RecoBTag/SecondaryVertex/interface/TrackSelector.h" + #include "fastjet/contrib/Njettiness.hh" #include "TMVA/Reader.h" // @@ -43,9 +53,6 @@ \date November 2010 */ - - - class HiInclusiveJetAnalyzer : public edm::EDAnalyzer { public: @@ -84,7 +91,7 @@ class HiInclusiveJetAnalyzer : public edm::EDAnalyzer { math::XYZPoint getPosition(const DetId &id, reco::Vertex::Point vtx = reco::Vertex::Point(0,0,0)); int TaggedJet(reco::Jet calojet, edm::Handle jetTags ); float getTau(unsigned num, const reco::GenJet object) const; - void analyzeSubjets(const reco::Jet jet); + void analyzeSubjets(const reco::Jet jet, int idx, edm::Handle, edm::Handle >, edm::Handle jetTags_CombinedSvtxV2, edm::Handle jetTags_negCombinedSvtxV2, edm::Handle jetTags_JP, edm::Handle > subjetTagInfo, edm::Handle > subjetSV, edm::Handle > subjetTagInfoNegSVx); void fillNewJetVarsRecoJet(const reco::Jet jet); void fillNewJetVarsRefJet(const reco::GenJet jet); void fillNewJetVarsGenJet(const reco::GenJet jet); @@ -92,7 +99,14 @@ class HiInclusiveJetAnalyzer : public edm::EDAnalyzer { void analyzeRefSubjets(const reco::GenJet jet); void analyzeGenSubjets(const reco::GenJet jet); float getAboveCharmThresh(reco::TrackRefVector& selTracks, const reco::TrackIPTagInfo& ipData, int sigOrVal); - + + int findMatchedParton(float eta, float phi, float maxDr, edm::Handle genparts, int partonFlavor); + int getFlavorProcess(int index, edm::Handle genparts); + + reco::TrackSelector trackSelector; + reco::TrackSelector trackPseudoSelector; + reco::V0Filter pseudoVertexV0Filter; + reco::V0Filter trackPairV0Filter; std::auto_ptr routine_; class ExtraInfo : public fastjet::PseudoJet::UserInfoBase { @@ -155,6 +169,8 @@ class HiInclusiveJetAnalyzer : public edm::EDAnalyzer { bool useQuality_; std::string trackQuality_; + bool isPythia6_; + bool doSubEvent_; double genPtMin_; bool doLifeTimeTagging_; @@ -210,6 +226,18 @@ class HiInclusiveJetAnalyzer : public edm::EDAnalyzer { edm::EDGetTokenT PositiveCombinedSecondaryVertexV2BJetTags_; edm::EDGetTokenT NegativeSoftPFMuonByPtBJetTags_; edm::EDGetTokenT PositiveSoftPFMuonByPtBJetTags_; + edm::EDGetTokenT SubjetJetProbabilityBJetTags_; + edm::EDGetTokenT > svImpactParameterTagInfos_; + edm::EDGetTokenT > svSubjetTagInfos_; + edm::EDGetTokenT > svSubjetNegTagInfos_; + edm::EDGetTokenT CombinedSubjetSecondaryVertexBJetTags_; + edm::EDGetTokenT CombinedSubjetNegativeSecondaryVertexBJetTags_; + + bool doExtendedFlavorTagging_; + edm::EDGetTokenT jetFlavourInfosToken_; + edm::EDGetTokenT subjetFlavourInfosToken_; + edm::EDGetTokenT > groomedJetsToken_; + bool useSubjets_; static const int MAXJETS = 1000; static const int MAXTRACKS = 5000; @@ -325,10 +353,38 @@ class HiInclusiveJetAnalyzer : public edm::EDAnalyzer { float jttau2[MAXJETS]; float jttau3[MAXJETS]; + float jtHadronFlavor[MAXJETS]; + float jtPartonFlavor[MAXJETS]; + std::vector> jtSubJetPt; std::vector> jtSubJetEta; std::vector> jtSubJetPhi; std::vector> jtSubJetM; + std::vector> jtSubJetHadronFlavor; + std::vector> jtSubJetPartonFlavor; + std::vector> jtSubJetcsvV2; + std::vector> jtSubJetNegCsvV2; + std::vector> jtSubJetJP; + std::vector> jtSubJetVtxType; + std::vector>> jtSubJetSvtxm; + std::vector>> jtSubJetSvtxpt; + std::vector>> jtSubJetSvtxeta; + std::vector>> jtSubJetSvtxphi; + std::vector>> jtSubJetSvtxNtrk; + std::vector>> jtSubJetSvtxdl; + std::vector>> jtSubJetSvtxdls; + + std::vector>> jtSubJetHadronDR; + std::vector>> jtSubJetHadronPt; + std::vector>> jtSubJetHadronEta; + std::vector>> jtSubJetHadronPhi; + std::vector>> jtSubJetHadronPdg; + std::vector>> jtSubJetPartonDR; + std::vector>> jtSubJetPartonPt; + std::vector>> jtSubJetPartonEta; + std::vector>> jtSubJetPartonPhi; + std::vector>> jtSubJetPartonPdg; + std::vector> jtConstituentsId; std::vector> jtConstituentsE; @@ -442,16 +498,17 @@ class HiInclusiveJetAnalyzer : public edm::EDAnalyzer { float pdiscr_csvV2[MAXJETS]; int nsvtx[MAXJETS]; - int svtxntrk[MAXJETS]; - float svtxdl[MAXJETS]; - float svtxdls[MAXJETS]; - float svtxdl2d[MAXJETS]; - float svtxdls2d[MAXJETS]; - float svtxm[MAXJETS]; - float svtxpt[MAXJETS]; + std::vector >svType; + std::vector >svtxntrk; + std::vector >svtxdl; + std::vector >svtxdls; + std::vector >svtxdl2d; + std::vector >svtxdls2d; + std::vector >svtxm; + std::vector >svtxpt; float svtxmcorr[MAXJETS]; float svtxnormchi2[MAXJETS]; - float svJetDeltaR[MAXJETS]; + std::vector >svJetDeltaR; float svtxTrkSumChi2[MAXJETS]; int svtxTrkNetCharge[MAXJETS]; int svtxNtrkInCone[MAXJETS]; @@ -504,6 +561,23 @@ class HiInclusiveJetAnalyzer : public edm::EDAnalyzer { float discr_jetID_cuts[MAXJETS]; float discr_jetID_bdt[MAXJETS]; + int refparton_flavorProcess[MAXJETS]; + float refGSP_gpt [MAXJETS]; + float refGSP_geta [MAXJETS]; + float refGSP_gphi [MAXJETS]; + float refGSP_gidx [MAXJETS]; + float refGSP_b1pt [MAXJETS]; + float refGSP_b1eta [MAXJETS]; + float refGSP_b1phi [MAXJETS]; + float refGSP_b2pt [MAXJETS]; + float refGSP_b2eta [MAXJETS]; + float refGSP_b2phi [MAXJETS]; + float refGSP_b1Match_jtdR [MAXJETS]; + float refGSP_b2Match_jtdR [MAXJETS]; + float refGSP_bbdR [MAXJETS]; + float refGSP_bbzg [MAXJETS]; + int refGSP_SubJtMatched [MAXJETS]; + float refpt[MAXJETS]; float refeta[MAXJETS]; float refphi[MAXJETS]; diff --git a/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_XeXe_mc.py b/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_XeXe_mc.py index c74d882bf4333..778df09c2e86c 100644 --- a/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_XeXe_mc.py +++ b/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_XeXe_mc.py @@ -89,8 +89,8 @@ akSoftDrop4GenJets + highPurityTracks + ak4PFJetSequence + - ak4CaloJetSequence + - akSoftDrop4PFJetSequence + +# ak4CaloJetSequence + +# akSoftDrop4PFJetSequence + akPu4PFJetSequence + akCs4PFJetSequence ) diff --git a/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_cleanedPbPb.py b/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_cleanedPbPb.py index 3c744c2d28e65..828cbb9680a8e 100644 --- a/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_cleanedPbPb.py +++ b/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_cleanedPbPb.py @@ -48,16 +48,16 @@ #akCsSoftDrop4PFJets + #akCsSoftDrop5PFJets + - akPu3CaloJetSequence + - akPu3PFJetSequence + - akCs3PFJetSequence + +# akPu3CaloJetSequence + +# akPu3PFJetSequence + +# akCs3PFJetSequence + - akPu4CaloJetSequence + - akPu4PFJetSequence + - akCs4PFJetSequence + +# akPu4CaloJetSequence + + akPu4PFJetSequence #+ +# akCs4PFJetSequence + - akPu5CaloJetSequence + - akPu5PFJetSequence# + +# akPu5CaloJetSequence + +# akPu5PFJetSequence# + #to be added later #akCsSoftDrop4PFJetSequence + diff --git a/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_nominalPP.py b/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_nominalPP.py index 291c5a9c8e6db..99263201c4632 100644 --- a/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_nominalPP.py +++ b/HeavyIonsAnalysis/JetAnalysis/python/FullJetSequence_nominalPP.py @@ -54,9 +54,9 @@ from HeavyIonsAnalysis.JetAnalysis.jets.ak3PFJetSequence_pp_mc_cff import * from HeavyIonsAnalysis.JetAnalysis.jets.ak4PFJetSequence_pp_mc_cff import * from HeavyIonsAnalysis.JetAnalysis.jets.ak5PFJetSequence_pp_mc_cff import * -from HeavyIonsAnalysis.JetAnalysis.jets.ak4CaloJetSequence_pp_mc_cff import * -from HeavyIonsAnalysis.JetAnalysis.jets.akSoftDrop4PFJetSequence_pp_mc_cff import * -from HeavyIonsAnalysis.JetAnalysis.jets.akSoftDrop5PFJetSequence_pp_mc_cff import * +#from HeavyIonsAnalysis.JetAnalysis.jets.ak4CaloJetSequence_pp_mc_cff import * +#from HeavyIonsAnalysis.JetAnalysis.jets.akSoftDrop4PFJetSequence_pp_mc_cff import * +#from HeavyIonsAnalysis.JetAnalysis.jets.akSoftDrop5PFJetSequence_pp_mc_cff import * highPurityTracks = cms.EDFilter("TrackSelector", src = cms.InputTag("generalTracks"), @@ -72,17 +72,17 @@ ak5GenJets + ak3PFJets + ak5PFJets + - akSoftDrop4PFJets + - akSoftDrop5PFJets + - akFilter4PFJets + - akFilter5PFJets + - akSoftDrop4GenJets + - akSoftDrop5GenJets + + #akSoftDrop4PFJets + + #akSoftDrop5PFJets + + #akFilter4PFJets + + #akFilter5PFJets + + #akSoftDrop4GenJets + + #akSoftDrop5GenJets + highPurityTracks + - ak3PFJetSequence + - ak4PFJetSequence + - ak5PFJetSequence + - ak4CaloJetSequence + - akSoftDrop4PFJetSequence + - akSoftDrop5PFJetSequence + #ak3PFJetSequence + + ak4PFJetSequence + #ak5PFJetSequence + + #ak4CaloJetSequence + + #akSoftDrop4PFJetSequence + + #akSoftDrop5PFJetSequence ) diff --git a/HeavyIonsAnalysis/JetAnalysis/python/bTaggers_cff.py b/HeavyIonsAnalysis/JetAnalysis/python/bTaggers_cff.py index 7e225d8aab209..782ed947e82ac 100644 --- a/HeavyIonsAnalysis/JetAnalysis/python/bTaggers_cff.py +++ b/HeavyIonsAnalysis/JetAnalysis/python/bTaggers_cff.py @@ -22,10 +22,13 @@ from RecoBTag.Configuration.RecoBTag_cff import * class bTaggers: - def __init__(self,jetname,rParam): + def __init__(self,jetname,rParam,ispp,doSubjets): self.JetTracksAssociatorAtVertex = ak5JetTracksAssociatorAtVertex.clone() self.JetTracksAssociatorAtVertex.jets = cms.InputTag(jetname+"Jets") - self.JetTracksAssociatorAtVertex.tracks = cms.InputTag("hiGeneralTracks") + if(ispp): + self.JetTracksAssociatorAtVertex.tracks = cms.InputTag("generalTracks") + else: + self.JetTracksAssociatorAtVertex.tracks = cms.InputTag("hiGeneralTracks") self.ImpactParameterTagInfos = impactParameterTagInfos.clone() self.ImpactParameterTagInfos.jetTracks = cms.InputTag(jetname+"JetTracksAssociatorAtVertex") self.ImpactParameterTagInfos.primaryVertex = cms.InputTag("offlinePrimaryVertices") @@ -38,18 +41,6 @@ def __init__(self,jetname,rParam): self.JetBProbabilityBJetTags = jetBProbabilityBJetTags.clone() self.JetBProbabilityBJetTags.tagInfos = cms.VInputTag(cms.InputTag(jetname+"ImpactParameterTagInfos")) - self.SecondaryVertexTagInfos = secondaryVertexTagInfos.clone() - self.SecondaryVertexTagInfos.trackIPTagInfos = cms.InputTag(jetname+"ImpactParameterTagInfos") - #self.SimpleSecondaryVertexBJetTags = simpleSecondaryVertexBJetTags.clone() - #self.SimpleSecondaryVertexBJetTags.tagInfos = cms.VInputTag(cms.InputTag(jetname+"SecondaryVertexTagInfos")) - self.CombinedSecondaryVertexBJetTags = combinedSecondaryVertexV2BJetTags.clone() - self.CombinedSecondaryVertexBJetTags.tagInfos = cms.VInputTag(cms.InputTag(jetname+"ImpactParameterTagInfos"), - cms.InputTag(jetname+"SecondaryVertexTagInfos")) - self.CombinedSecondaryVertexV2BJetTags = combinedSecondaryVertexV2BJetTags.clone() - self.CombinedSecondaryVertexV2BJetTags.tagInfos = cms.VInputTag(cms.InputTag(jetname+"ImpactParameterTagInfos"), - cms.InputTag(jetname+"SecondaryVertexTagInfos")) - - # secondary vertex b-tag self.SecondaryVertexTagInfos = secondaryVertexTagInfos.clone() self.SecondaryVertexTagInfos.trackIPTagInfos = cms.InputTag(jetname+"ImpactParameterTagInfos") @@ -103,6 +94,51 @@ def __init__(self,jetname,rParam): self.NegativeSoftPFMuonByPtBJetTags = negativeSoftPFMuonByPtBJetTags.clone() self.NegativeSoftPFMuonByPtBJetTags.tagInfos = cms.VInputTag(cms.InputTag(jetname+"SoftPFMuonsTagInfos")) + if doSubjets: + self.SubjetImpactParameterTagInfos = impactParameterTagInfos.clone() + self.SubjetImpactParameterTagInfos.jetTracks = cms.InputTag(jetname+"SubjetJetTracksAssociatorAtVertex") + self.SubjetJetProbabilityBJetTags = jetProbabilityBJetTags.clone() + self.SubjetJetProbabilityBJetTags.tagInfos = cms.VInputTag(cms.InputTag(jetname+"SubjetImpactParameterTagInfos")) + self.SubjetSecondaryVertexTagInfos = secondaryVertexTagInfos.clone() + self.SubjetSecondaryVertexTagInfos.trackIPTagInfos = cms.InputTag(jetname+'SubjetImpactParameterTagInfos') + self.SubjetSecondaryVertexTagInfos.fatJets = cms.InputTag('ak4PFJets') + self.SubjetSecondaryVertexTagInfos.groomedFatJets = cms.InputTag(jetname+'Jets') + self.SubjetSecondaryVertexTagInfos.useSVClustering = cms.bool(True) + self.SubjetSecondaryVertexTagInfos.useExternalSV = cms.bool(True) + self.SubjetSecondaryVertexTagInfos.jetAlgorithm = cms.string('AntiKt') + self.SubjetSecondaryVertexTagInfos.useSVMomentum = cms.bool(True) + self.SubjetSecondaryVertexTagInfos.rParam = cms.double(0.4) + self.SubjetSecondaryVertexTagInfos.extSVCollection = cms.InputTag('inclusiveSecondaryVertices') + self.SubjetSecondaryVertexTagInfos.vertexCuts.maxDeltaRToJetAxis = cms.double(0.1) + + self.SubjetSecondaryVertexNegativeTagInfos = self.SubjetSecondaryVertexTagInfos.clone() + self.SubjetSecondaryVertexNegativeTagInfos.vertexCuts.distVal2dMin = -2.5 + self.SubjetSecondaryVertexNegativeTagInfos.vertexCuts.distVal2dMax = -0.01 + self.SubjetSecondaryVertexNegativeTagInfos.vertexCuts.distSig2dMin = -99999.9 + self.SubjetSecondaryVertexNegativeTagInfos.vertexCuts.distSig2dMax = -3.0 + self.SubjetSecondaryVertexNegativeTagInfos.vertexCuts.maxDeltaRToJetAxis = -0.5 + + self.SubjetJetTracksAssociatorAtVertex = cms.EDProducer("JetTracksAssociatorExplicit", + jets = cms.InputTag(jetname+'Jets','SubJets') + ) + if ispp: + self.SubjetJetTracksAssociatorAtVertex.tracks = cms.InputTag('highPurityTracks') + else: + self.SubjetJetTracksAssociatorAtVertex.tracks = cms.InputTag('highPurityTracks') + + self.CombinedSubjetSecondaryVertexV2BJetTags = combinedSecondaryVertexV2BJetTags.clone( + tagInfos = cms.VInputTag(cms.InputTag(jetname+"SubjetImpactParameterTagInfos"), + cms.InputTag(jetname+"SubjetSecondaryVertexTagInfos")) + ) + self.CombinedSubjetSecondaryVertexBJetTags = combinedSecondaryVertexV2BJetTags.clone( + tagInfos = cms.VInputTag(cms.InputTag(jetname+"SubjetImpactParameterTagInfos"), + cms.InputTag(jetname+"SubjetSecondaryVertexTagInfos")) + ) + self.CombinedSubjetNegativeSecondaryVertexV2BJetTags = negativeCombinedSecondaryVertexV2BJetTags.clone( + tagInfos = cms.VInputTag(cms.InputTag(jetname+"SubjetImpactParameterTagInfos"), + cms.InputTag(jetname+"SubjetSecondaryVertexNegativeTagInfos")) + ) + self.JetTracksAssociator = cms.Sequence(self.JetTracksAssociatorAtVertex) self.JetBtaggingIP = cms.Sequence(self.ImpactParameterTagInfos * ( self.TrackCountingHighEffBJetTags + @@ -163,9 +199,15 @@ def __init__(self,jetname,rParam): rParam = rParam, bHadrons = cms.InputTag(jetname+"PatJetPartons","bHadrons"), cHadrons = cms.InputTag(jetname+"PatJetPartons","cHadrons"), - partons = cms.InputTag(jetname+"PatJetPartons","partons") + partons = cms.InputTag(jetname+"PatJetPartons","algorithmicPartons"), + leptons = cms.InputTag(jetname+"PatJetPartons","leptons"), ) + if doSubjets: + self.PatJetFlavourAssociation.jets="ak4PFJets" + self.PatJetFlavourAssociation.groomedJets = cms.InputTag(jetname+'Jets') + self.PatJetFlavourAssociation.subjets = cms.InputTag(jetname+'Jets', 'SubJets') + self.PatJetFlavourId = cms.Sequence(self.PatJetPartons*self.PatJetFlavourAssociation) #self.match = patJetGenJetMatch.clone( # src = cms.InputTag(jetname+"Jets"), diff --git a/HeavyIonsAnalysis/JetAnalysis/python/jets/makeJetSequences.sh b/HeavyIonsAnalysis/JetAnalysis/python/jets/makeJetSequences.sh index 8d9cfa4ccd99f..8a93a10f29fc6 100755 --- a/HeavyIonsAnalysis/JetAnalysis/python/jets/makeJetSequences.sh +++ b/HeavyIonsAnalysis/JetAnalysis/python/jets/makeJetSequences.sh @@ -1,6 +1,6 @@ #!/bin/sh -for system in PbPb pp pPb +for system in PbPb pp do for sample in data jec mc mb do @@ -18,11 +18,7 @@ do for object in PF Calo do # no Cs Calo or pp jets - - if ( [ $system == "pPb" ] ) && ( ( [ $radius -ne 4 ] && [ $radius -ne 3 ] ) || ( [ $sub != "NONE" ] && [ $sub != "Pu" ] ) || [ $groom != "NONE" ] ) ; then - continue - fi - if ( [ $object == "Calo" ] ) && ( [ $sub == "Cs" ] ) ; then + if ( [ $object == "Calo" ] ) && ( [ $sub == "Cs" ] ) ; then continue fi subt=$sub @@ -46,7 +42,8 @@ do resolveByDist="False" fi genjets="HiGenJets" - ismc="False" + ispp="False" + ismc="False" corrlabel="_offline" domatch="True" tracks="hiGeneralTracks" @@ -61,15 +58,16 @@ do jetcorrectionlevels="\'L2Relative\',\'L3Absolute\'" #echo "" > $algo$subt$radius${object}JetSequence_${system}_${sample}_cff.py - if [ $system == "pp" ] || [ $system == "pPb" ]; then + if [ $system == "pp" ]; then #corrlabel="_generalTracks" - corrlabel="_offline" tracks="generalTracks" vertex="offlinePrimaryVertices" genparticles="genParticles" partons="genParticles" pflow="particleFlow" - if [ $sample == "data" ] && [ $system != "pPb" ] && [ $sub == "NONE" ] && [ $radius == 4 ] && [ $object == "PF" ]; then + doTower="False" + ispp="True" + if [ $sample == "data" ] && [ $sub == "NONE" ] && [ $radius == 4 ] && [ $object == "PF" ]; then jetcorrectionlevels="\'L2Relative\',\'L3Absolute\',\'L2L3Residual\'" fi fi @@ -78,13 +76,11 @@ do ismc="True" fi - if [ $system == "pp" ] || [ $system == "pPb" ]; then + if [ $system == "pp" ]; then genjets="GenJets" matchGenjets="GenJets" fi - if [ $system == "pp" ]; then - doTower="False" - fi + if [ $sub == "Pu" ]; then corrname=`echo ${algo} | sed 's/\(.*\)/\U\1/'`${sub}${radius}${object}${corrlabel} else @@ -108,6 +104,7 @@ do -e "s/CORRNAME_/$corrname/g" \ -e "s/MATCHED_/$match/g" \ -e "s/ISMC/$ismc/g" \ + -e "s/ISPP/$ispp/g" \ -e "s/MATCHGENJETS/$matchGenjets/g" \ -e "s/GENJETS/$genjets/g" \ -e "s/GENPARTICLES/$genparticles/g" \ @@ -124,6 +121,13 @@ do -e "s/RESOLVEBYDIST_/$resolveByDist/g" \ > $algo$subt$groomt$radius${object}JetSequence_${system}_${sample}_cff.py + if [ $doSubJets == "True" ]; then + sed -i 's/\#SUBJETDUMMY_//g' $algo$subt$groomt$radius${object}JetSequence_${system}_${sample}_cff.py + fi + + if [ $ispp == "True" ]; then + sed -i 's/\#ppDummy_//g' $algo$subt$groomt$radius${object}JetSequence_${system}_${sample}_cff.py + fi # skip no sub if [ $sample == "jec" ]; then echo "${algo}${subt}${groomt}${radius}${object}JetAnalyzer.genPtMin = cms.untracked.double(1)" >> $algo$subt$groomt$radius${object}JetSequence_${system}_${sample}_cff.py diff --git a/HeavyIonsAnalysis/JetAnalysis/python/jets/templateSequence_bTag_cff.py.txt b/HeavyIonsAnalysis/JetAnalysis/python/jets/templateSequence_bTag_cff.py.txt index 250ce0f4c67fe..b04c3c8c2dbea 100644 --- a/HeavyIonsAnalysis/JetAnalysis/python/jets/templateSequence_bTag_cff.py.txt +++ b/HeavyIonsAnalysis/JetAnalysis/python/jets/templateSequence_bTag_cff.py.txt @@ -37,7 +37,7 @@ ALGO_SUB_GROOM_RADIUS_OBJECT_JetID= cms.EDProducer('JetIDProducer', JetIDParams, #ALGO_SUB_GROOM_RADIUS_OBJECT_clean = heavyIonCleanedGenJets.clone(src = cms.InputTag('ALGO_RADIUS_MATCHGENJETS')) -ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger = bTaggers("ALGO_SUB_GROOM_RADIUS_OBJECT_",0.RADIUS_) +ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger = bTaggers("ALGO_SUB_GROOM_RADIUS_OBJECT_",0.RADIUS_,ISPP,DOSUBJETS_) #create objects locally since they dont load properly otherwise #ALGO_SUB_GROOM_RADIUS_OBJECT_match = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.match @@ -83,9 +83,19 @@ ALGO_SUB_GROOM_RADIUS_OBJECT_SoftPFMuonByPtBJetTags = ALGO_SUB_GROOM_RADIUS_OBJE ALGO_SUB_GROOM_RADIUS_OBJECT_NegativeSoftPFMuonByPtBJetTags = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.NegativeSoftPFMuonByPtBJetTags ALGO_SUB_GROOM_RADIUS_OBJECT_PositiveSoftPFMuonByPtBJetTags = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.PositiveSoftPFMuonByPtBJetTags ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourIdLegacy = cms.Sequence(ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetPartonAssociationLegacy*ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourAssociationLegacy) -#Not working with our PU sub, but keep it here for reference -#ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourAssociation = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.PatJetFlavourAssociation -#ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourId = cms.Sequence(ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetPartons*ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourAssociation) +#Not working with our PU sub +ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourAssociation = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.PatJetFlavourAssociation +ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourId = cms.Sequence(ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetPartons*ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourAssociation) + +#adding the subjet taggers +#SUBJETDUMMY_ALGO_SUB_GROOM_RADIUS_OBJECT_SubjetImpactParameterTagInfos = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.SubjetImpactParameterTagInfos +#SUBJETDUMMY_ALGO_SUB_GROOM_RADIUS_OBJECT_SubjetJetProbabilityBJetTags = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.SubjetJetProbabilityBJetTags +#SUBJETDUMMY_ALGO_SUB_GROOM_RADIUS_OBJECT_SubjetSecondaryVertexTagInfos = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.SubjetSecondaryVertexTagInfos +#SUBJETDUMMY_ALGO_SUB_GROOM_RADIUS_OBJECT_SubjetSecondaryVertexNegativeTagInfos = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.SubjetSecondaryVertexNegativeTagInfos +#SUBJETDUMMY_ALGO_SUB_GROOM_RADIUS_OBJECT_SubjetJetTracksAssociatorAtVertex = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.SubjetJetTracksAssociatorAtVertex +#SUBJETDUMMY_ALGO_SUB_GROOM_RADIUS_OBJECT_CombinedSubjetSecondaryVertexBJetTags = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.CombinedSubjetSecondaryVertexBJetTags +#SUBJETDUMMY_ALGO_SUB_GROOM_RADIUS_OBJECT_CombinedSubjetSecondaryVertexV2BJetTags = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.CombinedSubjetSecondaryVertexV2BJetTags +#SUBJETDUMMY_ALGO_SUB_GROOM_RADIUS_OBJECT_CombinedSubjetNegativeSecondaryVertexV2BJetTags = ALGO_SUB_GROOM_RADIUS_OBJECT_bTagger.CombinedSubjetNegativeSecondaryVertexV2BJetTags ALGO_SUB_GROOM_RADIUS_OBJECT_JetBtaggingIP = cms.Sequence(ALGO_SUB_GROOM_RADIUS_OBJECT_ImpactParameterTagInfos * (ALGO_SUB_GROOM_RADIUS_OBJECT_TrackCountingHighEffBJetTags + @@ -139,10 +149,11 @@ ALGO_SUB_GROOM_RADIUS_OBJECT_patJetsWithBtagging = patJets.clone(jetSource = cms genJetMatch = cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_match"), genPartonMatch = cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_parton"), jetCorrFactorsSource = cms.VInputTag(cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_corr")), - JetPartonMapSource = cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourAssociationLegacy"), + #JetPartonMapSource = cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourAssociationLegacy"), + JetPartonMapSource = cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourAssociation"), JetFlavourInfoSource = cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourAssociation"), trackAssociationSource = cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_JetTracksAssociatorAtVertex"), - useLegacyJetMCFlavour = True, + useLegacyJetMCFlavour = False, discriminatorSources = cms.VInputTag(cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_SimpleSecondaryVertexHighEffBJetTags"), cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_SimpleSecondaryVertexHighPurBJetTags"), cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_CombinedSecondaryVertexBJetTags"), @@ -199,8 +210,13 @@ ALGO_SUB_GROOM_RADIUS_OBJECT_JetAnalyzer = inclusiveJetAnalyzer.clone(jetTag = c doSubJets = cms.untracked.bool(DOSUBJETS_), doGenSubJets = cms.untracked.bool(DOGENSUBJETS_), subjetGenTag = cms.untracked.InputTag("ALGO_GROOM_RADIUS_GenJets"), + doExtendedFlavorTagging = cms.untracked.bool(ISMC), + jetFlavourInfos = cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourAssociation"), + subjetFlavourInfos = cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourAssociation","SubJets"), + groomedJets = cms.InputTag("ALGO_SUB_GROOM_RADIUS_OBJECT_Jets"), + isPythia6 = cms.untracked.bool(False), doGenTaus = ISMC - ) + ) ALGO_SUB_GROOM_RADIUS_OBJECT_JetSequence_mc = cms.Sequence( #ALGO_SUB_GROOM_RADIUS_OBJECT_clean @@ -215,10 +231,10 @@ ALGO_SUB_GROOM_RADIUS_OBJECT_JetSequence_mc = cms.Sequence( * #ALGO_SUB_GROOM_RADIUS_OBJECT_JetID #* - ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourIdLegacy + #ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourIdLegacy # works for PbPb #* - #ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourId # Use legacy algo till PU implemented - * + #ppDummy_ALGO_SUB_GROOM_RADIUS_OBJECT_PatJetFlavourId # doesn't work for PbPb yet + #ppDummy_* ALGO_SUB_GROOM_RADIUS_OBJECT_JetTracksAssociatorAtVertex * ALGO_SUB_GROOM_RADIUS_OBJECT_JetBtagging diff --git a/HeavyIonsAnalysis/JetAnalysis/src/HiInclusiveJetAnalyzer.cc b/HeavyIonsAnalysis/JetAnalysis/src/HiInclusiveJetAnalyzer.cc index d8084f1e9ef42..1d312f0121c27 100644 --- a/HeavyIonsAnalysis/JetAnalysis/src/HiInclusiveJetAnalyzer.cc +++ b/HeavyIonsAnalysis/JetAnalysis/src/HiInclusiveJetAnalyzer.cc @@ -56,6 +56,7 @@ using namespace std; using namespace edm; using namespace reco; +//template class std::vector > >; // class ExtraInfo : public fastjet::PseudoJet::UserInfoBase { // public: @@ -68,7 +69,11 @@ using namespace reco; HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) : - geo(0) + trackSelector(iConfig.getParameter("trackSelection")), + trackPseudoSelector(iConfig.getParameter("trackSelection")), + pseudoVertexV0Filter(iConfig.getParameter("trackPairV0Filter")), + trackPairV0Filter(iConfig.getParameter("trackPairV0Filter")), + geo(0) { doMatch_ = iConfig.getUntrackedParameter("matchJets",false); @@ -79,7 +84,7 @@ HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) matchTagPat_ = consumes (iConfig.getUntrackedParameter("matchTag")); // vtxTag_ = iConfig.getUntrackedParameter("vtxTag",edm::InputTag("hiSelectedVertex")); - vtxTag_ = consumes > (iConfig.getUntrackedParameter("vtxTag",edm::InputTag("hiSelectedVertex"))); + vtxTag_ = consumes > (iConfig.getUntrackedParameter("vtxTag",edm::InputTag("offlinePrimaryVertices"))); // iConfig.getUntrackedParameter("vtxTag",edm::InputTag("hiSelectedVertex")); trackTag_ = consumes (iConfig.getParameter("trackTag")); useQuality_ = iConfig.getUntrackedParameter("useQuality",1); @@ -103,11 +108,13 @@ HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) // if (iConfig.exists("genTau3")) // tokenGenTau3_ = consumes >(iConfig.getParameter("genTau3")); // else useGenTaus = false; - + + isPythia6_ = iConfig.getUntrackedParameter("isPythia6",false); isMC_ = iConfig.getUntrackedParameter("isMC",false); useHepMC_ = iConfig.getUntrackedParameter ("useHepMC",false); fillGenJets_ = iConfig.getUntrackedParameter("fillGenJets",false); + doExtendedFlavorTagging_ = iConfig.getUntrackedParameter("doExtendedFlavorTagging",false); doTrigger_ = iConfig.getUntrackedParameter("doTrigger",false); doHiJetID_ = iConfig.getUntrackedParameter("doHiJetID",false); if(doHiJetID_) jetIDweightFile_ = iConfig.getUntrackedParameter("jetIDWeight","weights.xml"); @@ -121,11 +128,11 @@ HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) genjetTag_ = consumes > (iConfig.getParameter("genjetTag")); //genjetTag_ = consumes>(iConfig.getParameter("genjetTag")); if(useHepMC_) eventInfoTag_ = consumes (iConfig.getParameter("eventInfoTag")); - eventGenInfoTag_ = consumes (iConfig.getParameter("eventInfoTag")); + //eventGenInfoTag_ = consumes (iConfig.getParameter("eventInfoTag")); } verbose_ = iConfig.getUntrackedParameter("verbose",false); - useVtx_ = iConfig.getUntrackedParameter("useVtx",false); + useVtx_ = iConfig.getUntrackedParameter("useVtx",true); useJEC_ = iConfig.getUntrackedParameter("useJEC",true); usePat_ = iConfig.getUntrackedParameter("usePAT",true); @@ -148,6 +155,13 @@ HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) doExtraCTagging_ = iConfig.getUntrackedParameter("doExtraCTagging",false); genParticleSrc_ = consumes (iConfig.getUntrackedParameter("genParticles",edm::InputTag("hiGenParticles"))); + jetFlavourInfosToken_ = consumes( iConfig.getParameter("jetFlavourInfos") ); + + if(doSubJets_ && doExtendedFlavorTagging_){ + subjetFlavourInfosToken_ = mayConsume( iConfig.exists("subjetFlavourInfos") ? iConfig.getParameter("subjetFlavourInfos") : edm::InputTag() ); + groomedJetsToken_ = mayConsume >( iConfig.exists("groomedJets") ? iConfig.getParameter("groomedJets") : edm::InputTag() ); + useSubjets_ = ( iConfig.exists("subjetFlavourInfos") && iConfig.exists("groomedJets") ); + } if(doTrigger_){ L1gtReadout_ = consumes< L1GlobalTriggerReadoutRecord > (iConfig.getParameter("L1gtReadout")); @@ -185,6 +199,14 @@ HiInclusiveJetAnalyzer::HiInclusiveJetAnalyzer(const edm::ParameterSet& iConfig) PositiveCombinedSecondaryVertexV2BJetTags_ = consumes (iConfig.getUntrackedParameter("PositiveCombinedSecondaryVertexV2BJetTags",(bTagJetName_+"PositiveCombinedSecondaryVertexV2BJetTags"))); NegativeSoftPFMuonByPtBJetTags_ = consumes (iConfig.getUntrackedParameter("NegativeSoftPFMuonByPtBJetTags",(bTagJetName_+"NegativeSoftPFMuonByPtBJetTags"))); PositiveSoftPFMuonByPtBJetTags_ = consumes (iConfig.getUntrackedParameter("PositiveSoftPFMuonByPtBJetTags",(bTagJetName_+"PositiveSoftPFMuonByPtBJetTags"))); + if(doSubJets_){ + CombinedSubjetSecondaryVertexBJetTags_ = mayConsume (iConfig.getUntrackedParameter("CombinedSubjetSecondaryVertexV2BJetTags",(bTagJetName_+"CombinedSubjetSecondaryVertexV2BJetTags"))); + SubjetJetProbabilityBJetTags_ = mayConsume (iConfig.getUntrackedParameter("SubjetJetProbabilityBJetTags",(bTagJetName_+"SubjetJetProbabilityBJetTags"))); + svSubjetTagInfos_ = mayConsume >(iConfig.getUntrackedParameter("SecondaryVertexTagInfos",(bTagJetName_+"SubjetSecondaryVertexTagInfos"))); + svSubjetNegTagInfos_ = mayConsume >(iConfig.getUntrackedParameter("NegativeSecondaryVertexTagInfos",(bTagJetName_+"SubjetNegativeSecondaryVertexTagInfos"))); + svImpactParameterTagInfos_ = mayConsume >(iConfig.getUntrackedParameter("ImpactParameterTagInfos",(bTagJetName_+"SubjetImpactParameterTagInfos"))); + CombinedSubjetNegativeSecondaryVertexBJetTags_ = mayConsume (iConfig.getUntrackedParameter("CombinedSubjetNegativeSecondaryVertexV2BJetTags",(bTagJetName_+"CombinedSubjetNegativeSecondaryVertexV2BJetTags"))); + } } // cout<<" jet collection : "<Branch("jttau1",jets_.jttau1,"jttau1[nref]/F"); t->Branch("jttau2",jets_.jttau2,"jttau2[nref]/F"); t->Branch("jttau3",jets_.jttau3,"jttau3[nref]/F"); - + + if(doExtendedFlavorTagging_){ + t->Branch("jtHadronFlavor",jets_.jtHadronFlavor,"jtHadronFlavor[nref]/F"); + t->Branch("jtPartonFlavor",jets_.jtPartonFlavor,"jtPartonFlavor[nref]/F"); + } + if(doSubJets_) { t->Branch("jtSubJetPt",&jets_.jtSubJetPt); t->Branch("jtSubJetEta",&jets_.jtSubJetEta); t->Branch("jtSubJetPhi",&jets_.jtSubJetPhi); t->Branch("jtSubJetM",&jets_.jtSubJetM); + t->Branch("jtSubJetcsvV2",&jets_.jtSubJetcsvV2); + t->Branch("jtSubJetNegCsvV2",&jets_.jtSubJetNegCsvV2); + t->Branch("jtSubJetJP",&jets_.jtSubJetJP); + t->Branch("jtSubJetVtxType", &jets_.jtSubJetVtxType); + t->Branch("jtSubJetSvtxm",&jets_.jtSubJetSvtxm); + t->Branch("jtSubJetSvtxpt",&jets_.jtSubJetSvtxpt); + t->Branch("jtSubJetSvtxeta",&jets_.jtSubJetSvtxeta); + t->Branch("jtSubJetSvtxphi",&jets_.jtSubJetSvtxphi); + t->Branch("jtSubJetSvtxNtrk",&jets_.jtSubJetSvtxNtrk); + t->Branch("jtSubJetSvtxdl",&jets_.jtSubJetSvtxdl); + t->Branch("jtSubJetSvtxdls",&jets_.jtSubJetSvtxdls); + if(doExtendedFlavorTagging_){ + t->Branch("jtSubJetHadronFlavor",&jets_.jtSubJetHadronFlavor); + t->Branch("jtSubJetPartonFlavor",&jets_.jtSubJetPartonFlavor); + t->Branch("jtSubJetHadronDR",&jets_.jtSubJetHadronDR); + t->Branch("jtSubJetHadronPt",&jets_.jtSubJetHadronPt); + t->Branch("jtSubJetHadronEta",&jets_.jtSubJetHadronEta); + t->Branch("jtSubJetHadronPhi",&jets_.jtSubJetHadronPhi); + t->Branch("jtSubJetHadronPdg",&jets_.jtSubJetHadronPdg); + t->Branch("jtSubJetPartonDR",&jets_.jtSubJetPartonDR); + t->Branch("jtSubJetPartonPt",&jets_.jtSubJetPartonPt); + t->Branch("jtSubJetPartonEta",&jets_.jtSubJetPartonEta); + t->Branch("jtSubJetPartonPhi",&jets_.jtSubJetPartonPhi); + t->Branch("jtSubJetPartonPdg",&jets_.jtSubJetPartonPdg); + } } if(doJetConstituents_){ @@ -473,13 +525,14 @@ HiInclusiveJetAnalyzer::beginJob() { t->Branch("pdiscr_csvV2",jets_.pdiscr_csvV2,"pdiscr_csvV2[nref]/F"); t->Branch("nsvtx", jets_.nsvtx, "nsvtx[nref]/I"); - t->Branch("svtxntrk", jets_.svtxntrk, "svtxntrk[nref]/I"); - t->Branch("svtxdl", jets_.svtxdl, "svtxdl[nref]/F"); - t->Branch("svtxdls", jets_.svtxdls, "svtxdls[nref]/F"); - t->Branch("svtxdl2d", jets_.svtxdl2d, "svtxdl2d[nref]/F"); - t->Branch("svtxdls2d", jets_.svtxdls2d, "svtxdls2d[nref]/F"); - t->Branch("svtxm", jets_.svtxm, "svtxm[nref]/F"); - t->Branch("svtxpt", jets_.svtxpt, "svtxpt[nref]/F"); + t->Branch("svType", &jets_.svType); + t->Branch("svtxntrk", &jets_.svtxntrk); + t->Branch("svtxdl", &jets_.svtxdl); + t->Branch("svtxdls", &jets_.svtxdls); + t->Branch("svtxdl2d", &jets_.svtxdl2d); + t->Branch("svtxdls2d", &jets_.svtxdls2d); + t->Branch("svtxm", &jets_.svtxm); + t->Branch("svtxpt", &jets_.svtxpt); t->Branch("svtxmcorr", jets_.svtxmcorr, "svtxmcorr[nref]/F"); t->Branch("nIPtrk",jets_.nIPtrk,"nIPtrk[nref]/I"); @@ -514,7 +567,7 @@ HiInclusiveJetAnalyzer::beginJob() { t->Branch("svtxTrkNetCharge",jets_.svtxTrkNetCharge,"svtxTrkNetCharge[nref]/I"); t->Branch("svtxNtrkInCone",jets_.svtxNtrkInCone,"svtxNtrkInCone[nref]/I"); - t->Branch("svJetDeltaR", jets_.svJetDeltaR, "svJetDeltaR[nref]/F"); + t->Branch("svJetDeltaR", &jets_.svJetDeltaR); t->Branch("trackSip2dSigAboveCharm",jets_.trackSip2dSigAboveCharm, "trackSip2dSigAboveCharm[nref]/F"); t->Branch("trackSip3dSigAboveCharm",jets_.trackSip3dSigAboveCharm, "trackSip3dSigAboveCharm[nref]/F"); t->Branch("trackSip2dValAboveCharm",jets_.trackSip2dValAboveCharm, "trackSip2dValAboveCharm[nref]/F"); @@ -640,6 +693,37 @@ HiInclusiveJetAnalyzer::beginJob() { t->Branch("refparton_flavor",jets_.refparton_flavor,"refparton_flavor[nref]/I"); t->Branch("refparton_flavorForB",jets_.refparton_flavorForB,"refparton_flavorForB[nref]/I"); + if(isPythia6_){ + t->Branch("refparton_flavorProcess",jets_.refparton_flavorProcess,"refparton_flavorProcess[nref]/I"); + } + + if(isPythia6_){ // for GSP to bb identify for a event, not for a jet + // t->Branch("refGSP_Channel" ,jets_.refGSP_Channel, "refGSP_Channel/I"); + t->Branch("refGSP_gpt" ,jets_.refGSP_gpt, "refGSP_gpt[nref]/F"); + t->Branch("refGSP_geta" ,jets_.refGSP_geta, "refGSP_geta[nref]/F"); + t->Branch("refGSP_gphi" ,jets_.refGSP_gphi, "refGSP_gphi[nref]/F"); + t->Branch("refGSP_gidx" ,jets_.refGSP_gidx, "refGSP_gidx[nref]/F"); + t->Branch("refGSP_b1pt" ,jets_.refGSP_b1pt, "refGSP_b1pt[nref]/F"); + t->Branch("refGSP_b1eta" ,jets_.refGSP_b1eta, "refGSP_b1eta[nref]/F"); + t->Branch("refGSP_b1phi" ,jets_.refGSP_b1phi, "refGSP_b1phi[nref]/F"); + t->Branch("refGSP_b2pt" ,jets_.refGSP_b2pt, "refGSP_b2pt[nref]/F"); + t->Branch("refGSP_b2eta" ,jets_.refGSP_b2eta, "refGSP_b2eta[nref]/F"); + t->Branch("refGSP_b2phi" ,jets_.refGSP_b2phi, "refGSP_b2phi[nref]/F"); + // t->Branch("refGSP_b1Match_jtIdx" ,jets_.refGSP_b1Match_jtIdx, "refGSP_b1Match_jtIdx/I"); + t->Branch("refGSP_b1Match_jtdR" ,jets_.refGSP_b1Match_jtdR, "refGSP_b1Match_jtdR[nref]/F"); + // t->Branch("refGSP_b1Match_Subjt1dR" ,jets_.refGSP_b1Match_Subjt1dR, "refGSP_b1Match_Subjt1dR/F"); + // t->Branch("refGSP_b1Match_Subjt2dR" ,jets_.refGSP_b1Match_Subjt2dR, "refGSP_b1Match_Subjt2dR/F"); + // t->Branch("refGSP_b2Match_jtIdx" ,jets_.refGSP_b2Match_jtIdx, "refGSP_b2Match_jtIdx/I"); + t->Branch("refGSP_b2Match_jtdR" ,jets_.refGSP_b2Match_jtdR, "refGSP_b2Match_jtdR[nref]/F"); + // t->Branch("refGSP_b1Match_Subjt1dR" ,jets_.refGSP_b1Match_Subjt1dR, "refGSP_b1Match_Subjt1dR/F"); + // t->Branch("refGSP_b1Match_Subjt2dR" ,jets_.refGSP_b1Match_Subjt2dR, "refGSP_b1Match_Subjt2dR/F"); + t->Branch("refGSP_bbdR" ,jets_.refGSP_bbdR, "refGSP_bbdR[nref]/F"); + // t->Branch("refGSP_bbdxy" ,jets_.refGSP_bbdxy, "refGSP_bbdxy/F"); + // t->Branch("refGSP_bbdz" ,jets_.refGSP_bbdz, "refGSP_bbdz/F"); + t->Branch("refGSP_bbzg" ,jets_.refGSP_bbzg, "refGSP_bbzg[nref]/F"); + t->Branch("refGSP_SubJtMatched" ,jets_.refGSP_SubJtMatched, "refGSP_SubJtMatched[nref]/I"); + } + if(doGenSubJets_) { t->Branch("refptG",jets_.refptG,"refptG[nref]/F"); t->Branch("refetaG",jets_.refetaG,"refetaG[nref]/F"); @@ -881,6 +965,10 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent, jets_.vz=vertex->begin()->z(); vtx = vertex->begin()->position(); } + //cout << " all vertices: "<< endl; + for(unsigned ivtx=0; ivtxsize(); ivtx++){ + // cout << " vtx x: " << vertex->at(ivtx).x() << " y: "<< vertex->at(ivtx).y() << " z: " << vertex->at(ivtx).z() << endl; + } } edm::Handle patjets; @@ -926,6 +1014,7 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent, Handle > tagInfo; Handle jetTags_JP; Handle jetTags_JB; + Handle jetTags_subjet_JP; //------------------------------------------------------ // Secondary vertex taggers @@ -942,6 +1031,12 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent, Handle jetTags_CombinedSvtxV2; Handle jetTags_negCombinedSvtxV2; Handle jetTags_posCombinedSvtxV2; + + Handle jetTags_subjet_combinedSvtx; + Handle > subjetTagInfo; + Handle > subjetTagInfoSVx; + Handle > subjetTagInfoNegSVx; + Handle jetTags_subjet_negCombinedSvtx; //------------------------------------------------------ // Soft muon tagger @@ -969,11 +1064,18 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent, iEvent.getByToken(CombinedSecondaryVertexV2BJetTags_, jetTags_CombinedSvtxV2); iEvent.getByToken(NegativeCombinedSecondaryVertexV2BJetTags_, jetTags_negCombinedSvtxV2); iEvent.getByToken(PositiveCombinedSecondaryVertexV2BJetTags_, jetTags_posCombinedSvtxV2); + if(doSubJets_){ + iEvent.getByToken(CombinedSubjetSecondaryVertexBJetTags_, jetTags_subjet_combinedSvtx); + iEvent.getByToken(CombinedSubjetNegativeSecondaryVertexBJetTags_, jetTags_subjet_negCombinedSvtx); + iEvent.getByToken(SubjetJetProbabilityBJetTags_, jetTags_subjet_JP); + iEvent.getByToken(svImpactParameterTagInfos_,subjetTagInfo); + iEvent.getByToken(svSubjetTagInfos_,subjetTagInfoSVx); + iEvent.getByToken(svSubjetNegTagInfos_, subjetTagInfoNegSVx); + } //iEvent.getByToken(NegativeSoftPFMuonByPtBJetTags_, jetTags_softMuneg); //iEvent.getByToken(PositiveSoftPFMuonByPtBJetTags_, jetTags_softMu); } - // get tower information edm::Handle towers; if(doTower){ @@ -1046,31 +1148,45 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent, const SecondaryVertexTagInfo &tagInfoSV = (*tagInfoSVx)[ith_tagged]; jets_.nsvtx[jets_.nref] = tagInfoSV.nVertices(); - if ( jets_.nsvtx[jets_.nref] > 0) { - - jets_.svtxntrk[jets_.nref] = tagInfoSV.nVertexTracks(0); - if(doExtraCTagging_) jets_.svJetDeltaR[jets_.nref] = reco::deltaR(tagInfoSV.flightDirection(0),jetDir); + std::vector svtxntrks; + std::vector svjetDr, svtxdl, svtxdls, svtxdl2d, svtxdls2d, svtxm, svtxpt; + std::vector svType; + if(jets_.nsvtx[jets_.nref]==0){ + svtxntrks.push_back(-999); + svjetDr.push_back(-999); + svtxdl.push_back(-999); + svtxdls.push_back(-999); + svtxdl2d.push_back(-999); + svtxm.push_back(-999); + svtxpt.push_back(-999); + svType.push_back(-999); + } + for (int ivtx=0; ivtx1) cout << "ntrks: " << tagInfoSV.nVertexTracks(ivtx) << endl; + svtxntrks.push_back(tagInfoSV.nVertexTracks(ivtx)); + if(doExtraCTagging_) svjetDr.push_back(reco::deltaR(tagInfoSV.flightDirection(ivtx),jetDir)); // this is the 3d flight distance, for 2-D use (0,true) - Measurement1D m1D = tagInfoSV.flightDistance(0); - jets_.svtxdl[jets_.nref] = m1D.value(); - jets_.svtxdls[jets_.nref] = m1D.significance(); - Measurement1D m2D = tagInfoSV.flightDistance(0,true); - jets_.svtxdl2d[jets_.nref] = m2D.value(); - jets_.svtxdls2d[jets_.nref] = m2D.significance(); - const Vertex& svtx = tagInfoSV.secondaryVertex(0); + Measurement1D m1D = tagInfoSV.flightDistance(ivtx); + svtxdl.push_back(m1D.value()); + svtxdls.push_back(m1D.significance()); + Measurement1D m2D = tagInfoSV.flightDistance(ivtx,true); + svtxdl2d.push_back(m2D.value()); + svtxdls2d.push_back(m2D.significance()); + const Vertex& svtx = tagInfoSV.secondaryVertex(ivtx); double svtxM = svtx.p4().mass(); double svtxPt = svtx.p4().pt(); - jets_.svtxm[jets_.nref] = svtxM; - jets_.svtxpt[jets_.nref] = svtxPt; + //if(jets_.nsvtx[jets_.nref]>1) cout << "svtxm: "<< svtxM << " svtxpt: "<< svtxPt << endl; + svtxm.push_back(svtxM); + svtxpt.push_back(svtxPt); if(doExtraCTagging_){ double trkSumChi2=0; int trkNetCharge=0; int nTrkInCone=0; int nsvtxtrks=0; - TrackRefVector svtxTracks = tagInfoSV.vertexTracks(0); + TrackRefVector svtxTracks = tagInfoSV.vertexTracks(ivtx); for(unsigned int itrk=0; itrkchi2(); @@ -1086,10 +1202,33 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent, double sinth = svtx.p4().Vect().Unit().Cross(tagInfoSV.flightDirection(0).unit()).Mag2(); sinth = sqrt(sinth); jets_.svtxmcorr[jets_.nref] = sqrt(pow(svtxM,2)+(pow(svtxPt,2)*pow(sinth,2)))+svtxPt*sinth; - } + } + TrackRefVector svtxTracks = tagInfoSV.vertexTracks(ivtx); + unsigned nVtxTypeTrks=0, nVtxTypeVertices=0, nAssocTrks=0; + for(unsigned int itrk=0; itrk= 0.5) nAssocTrks++; //0.5 is default value + } + if(nAssocTrks>0){ + bool isTrackVertex = (nAssocTrks==1); + ++*(isTrackVertex ? &nVtxTypeTrks : &nVtxTypeVertices); // <-- holy complicated logic, batman! From GhostTrackComputer.cc + } + if(nVtxTypeVertices) svType.push_back(0); + else if(nVtxTypeTrks) svType.push_back(1); + else svType.push_back(2); + if(svtx.ndof()>0)jets_.svtxnormchi2[jets_.nref] = svtx.chi2()/svtx.ndof(); } + jets_.svType.push_back(svType); + jets_.svtxntrk.push_back(svtxntrks); + jets_.svJetDeltaR.push_back(svjetDr); + jets_.svtxdl.push_back(svtxdl); + jets_.svtxdls.push_back(svtxdls); + jets_.svtxdl2d.push_back(svtxdl2d); + jets_.svtxdls2d.push_back(svtxdls2d); + jets_.svtxm.push_back(svtxm); + jets_.svtxpt.push_back(svtxpt); + } ith_tagged = this->TaggedJet(jet,jetTags_negSvtxHighEff); if(ith_tagged >= 0) jets_.ndiscr_ssvHighEff[jets_.nref] = (*jetTags_negSvtxHighEff)[ith_tagged].second; @@ -1139,8 +1278,6 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent, jets_.ipNHit[jets_.nIP + it] = selTracks[it]->numberOfValidHits(); jets_.ipNHitPixel[jets_.nIP + it] = selTracks[it]->hitPattern().numberOfValidPixelHits(); jets_.ipNHitStrip[jets_.nIP + it] = selTracks[it]->hitPattern().numberOfValidStripHits(); - //K.Jung fix for 92X compilation - //jets_.ipIsHitL1[jets_.nIP + it] = selTracks[it]->hitPattern().hasValidHitInFirstPixelBarrel(); jets_.ipIsHitL1[jets_.nIP + it] = selTracks[it]->hitPattern().hasValidHitInPixelLayer(PixelSubdetector::SubDetector(1),1); //(inBarrel, layer) jets_.ipProb0[jets_.nIP + it] = tagInfoIP.probabilities(0)[it]; jets_.ip2d[jets_.nIP + it] = data.ip2d.value(); @@ -1481,7 +1618,28 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent, jets_.jttau2[jets_.nref] = -999.; jets_.jttau3[jets_.nref] = -999.; - if(doSubJets_) analyzeSubjets(jet); + edm::Handle theJetFlavourInfos; + edm::Handle theSubjetFlavourInfos; + edm::Handle > groomedJets; + if(doExtendedFlavorTagging_){ + iEvent.getByToken(jetFlavourInfosToken_, theJetFlavourInfos ); + if(useSubjets_){ + iEvent.getByToken(subjetFlavourInfosToken_, theSubjetFlavourInfos); + iEvent.getByToken(groomedJetsToken_, groomedJets); + } + + int ijet=0; + for(reco::JetFlavourInfoMatchingCollection::const_iterator j= theJetFlavourInfos->begin(); j!=theJetFlavourInfos->end(); j++, ijet++){ + const reco::Jet *flavorMatch = (*j).first.get(); +// if(abs(flavorMatch->pt() - jets_.rawpt[jets_.nref]) < 0.1){ + if( sqrt(reco::deltaR2(flavorMatch->eta(), flavorMatch->phi(), jet.eta(), jet.phi() )) <0.05 ){ + jets_.jtHadronFlavor[jets_.nref] = (*j).second.getHadronFlavour(); + jets_.jtPartonFlavor[jets_.nref] = (*j).second.getPartonFlavour(); + } + } + } + + if(doSubJets_) analyzeSubjets(jet, jets_.nref, theSubjetFlavourInfos, groomedJets, jetTags_subjet_combinedSvtx, jetTags_subjet_negCombinedSvtx, jetTags_subjet_JP, subjetTagInfo, subjetTagInfoSVx, subjetTagInfoNegSVx); if(usePat_){ if( (*patjets)[j].hasUserFloat(jetName_+"Njettiness:tau1") ) @@ -1737,6 +1895,92 @@ HiInclusiveJetAnalyzer::analyze(const Event& iEvent, jets_.reftau3[jets_.nref] = -999.; jets_.refparton_flavorForB[jets_.nref] = (*patjets)[j].partonFlavour(); + + //Matt's flavor production code + if(isPythia6_){ + //int partonFlavor = (*patjets)[j].partonFlavour(); + // jets_.refparton_flavorForB[jets_.nref] = partonFlavor; + + if(jets_.jtHadronFlavor[jets_.nref]==4||jets_.jtHadronFlavor[jets_.nref]==5){ + + int partonMatchIndex = findMatchedParton(jet.eta(), jet.phi(), 0.3, genparts, jets_.jtHadronFlavor[jets_.nref]); + if(partonMatchIndex<0){ + cout<< "jet " << jets_.nref << " declared flavor " << jets_.jtHadronFlavor[jets_.nref] << " couldn't find the parton "<eta() , (float)(*genparts)[partonMatchIndex].mother(0)->phi(),(float)0.001,genparts, 0); + if (momIndex >=0 && (*genparts)[momIndex].pdgId()==21){ // double check momidex && momId = gluon + for( UInt_t i_dau =0; i_dau < (*genparts)[momIndex].numberOfDaughters(); i_dau++){ + if( (*genparts)[momIndex].daughter(i_dau)->pdgId()== -MatchedPartonId ){ + partonPairIndex=findMatchedParton( (*genparts)[momIndex].daughter(i_dau)->eta(), (*genparts)[momIndex].daughter(i_dau)->phi(),0.001,genparts , 0); + break; // terminate for i_dau + } + } // end i_dau < (*genparts)[momIndex].numberOfDaughters() + } // end if momIndex >=0 && (*genparts)[momIndex].pdgId()==21 + + if(partonPairIndex<=0) {cout<<"event : "<maxDr) continue; + if(genCand.pt() > highestPartonPt){ + matchIndex = i; + highestPartonPt = genCand.pt(); + } + } + return matchIndex; +} + +int HiInclusiveJetAnalyzer::getFlavorProcess(int index, Handle genparts){ + + const GenParticle & constParton = (*genparts)[ index ]; + GenParticle *matchedParton = const_cast(&constParton); + if(matchedParton->numberOfMothers()>1) cout<<" too many parents "<numberOfMothers()==0){ /*cout <<" found primary quark of pdg: " << matchedParton->pdgId() << endl;*/ return 1;} + int momID = matchedParton->mother(0)->pdgId(); + int momIndex = findMatchedParton(matchedParton->mother(0)->eta(),matchedParton->mother(0)->phi(),0.001,genparts); + if(momIndex==index){ cout << " WARNING! Particle is its own mother??" << endl; return -1; } + int gmomIdx=-1; + if(matchedParton->mother(0)->numberOfMothers()>0){ + gmomIdx = findMatchedParton(matchedParton->mother(0)->mother(0)->eta(),matchedParton->mother(0)->mother(0)->phi(),0.001,genparts); + if(gmomIdx==index) { cout << "WARNING! Particle is its own grandmother" << endl; return -1; } + } + + while(abs(matchedParton->mother(0)->pdgId())==5 && matchedParton->numberOfMothers()>0){ + if(matchedParton->mother(0)->status()==21) return 1; // primary b-quark + else { + matchedParton = (reco::GenParticle*)(reco::LeafCandidate*)(matchedParton->mother(0)); + //return getFlavorProcess(momIndex,genparts); + } + } + if(abs(momID)==5 && matchedParton->numberOfMothers()==0){ cout << "warning - found orphaned b-quark!" << endl; return -1; } + + if(!isPythia6_){ + if(matchedParton->status()>20 && matchedParton->status()<30) return 2; + if(matchedParton->status()>30 && matchedParton->status()<40) return 3; + if(matchedParton->status()>40 && matchedParton->status()<50) return 4; + if(matchedParton->status()>50){ + if(abs(momID)==21) return 5; + return 6; + } + } + else{ + if(momIndex<2) return 4; // initial state GSP + else if(momIndex<4) return 3; // sometimes GSP, sometimes associated to FEX + else if(momIndex<6) return 2; // primaries + else if(momIndex<8){ + if(momID==21) return 6; // final state hard GSP + else return 5; // final state soft GSP + } + } + //else cout<<" should never get here "<<" momID "< theSubjetFlavourInfos, edm::Handle > groomedJets, Handle jetTags_CombinedSvtx, Handle jetTags_negCombinedSvtx, Handle jetTags_JP, Handle > subjetIP, Handle > subjetSV, Handle > subjetNegSV) { std::vector sjpt; std::vector sjeta; std::vector sjphi; std::vector sjm; + std::vector hadronFlavor; + std::vector partonFlavor; + std::vector csvV2; + std::vector negCsvV2; + std::vector jptag; + std::vector sjVtxType; + std::vector> subjetVtxMass; + std::vector> subjetVtxPt; + std::vector> subjetVtxEta; + std::vector> subjetVtxPhi; + std::vector> subjetVtxNtrk; + std::vector> subjetVtxdl; + std::vector> subjetVtxdls; + std::vector> hadronDR; + std::vector> hadronPt; + std::vector> hadronEta; + std::vector> hadronPhi; + std::vector> hadronPdg; + std::vector> partonDR; + std::vector> partonPt; + std::vector> partonEta; + std::vector> partonPhi; + std::vector> partonPdg; + if(jet.numberOfDaughters()>0) { - for (unsigned k = 0; k < jet.numberOfDaughters(); ++k) { - const reco::Candidate & dp = *jet.daughter(k); - sjpt.push_back(dp.pt()); - sjeta.push_back(dp.eta()); - sjphi.push_back(dp.phi()); - sjm.push_back(dp.mass()); - } + //cout << "nsubjets: "<< jet.numberOfDaughters() << endl; + for (unsigned k = 0; k < jet.numberOfDaughters(); ++k) { + const reco::Candidate & dp = *jet.daughter(k); + sjpt.push_back(dp.pt()); + sjeta.push_back(dp.eta()); + sjphi.push_back(dp.phi()); + sjm.push_back(dp.mass()); + + vector svm, svpt, sveta, svphi, svntrk, svdl, svdls; + const reco::Jet *subjd = dynamic_cast(jet.daughter(k)); + int ith_tagged = TaggedJet(*subjd,jetTags_CombinedSvtx); + if(ith_tagged >= 0) csvV2.push_back((*jetTags_CombinedSvtx)[ith_tagged].second); + else csvV2.push_back(-1.); + ith_tagged = TaggedJet(*subjd,jetTags_JP); + if(ith_tagged >= 0) jptag.push_back((*jetTags_JP)[ith_tagged].second); + else jptag.push_back(-1.); + ith_tagged = TaggedJet(*subjd,jetTags_negCombinedSvtx); + if(ith_tagged >= 0) negCsvV2.push_back((*jetTags_negCombinedSvtx)[ith_tagged].second); + else negCsvV2.push_back(-1.); + + ith_tagged = TaggedJet(*subjd,jetTags_CombinedSvtx); + if(ith_tagged>=0){ + const SecondaryVertexTagInfo &subjetInfo = (*subjetSV)[ith_tagged]; + for(unsigned isv=0; isv 0.9 && subjetInfo.nVertices()>0){ +// cout << "full vtx criteria met for subjet " << k << " (pt=" << subjd->pt() << ", csv="<< (*jetTags_CombinedSvtx)[ith_tagged].second << ") " << endl; +// cout << "full vtx reco tracks: "<< endl; + int itrack=0; + for(reco::Vertex::trackRef_iterator it = subjetInfo.secondaryVertex(0).tracks_begin(); it != subjetInfo.secondaryVertex(0).tracks_end(); it++, itrack++){ +// cout << "sv track " << itrack << " pt: "<< (*it)->pt() << " eta: " << (*it)->eta() << " phi: "<< (*it)->phi() << endl; + } + + } + //if((*jetTags_CombinedSvtx)[ith_tagged].second > 0.9 && subjetInfo.nVertices()==0) cout << "WARNING! No vertex found in a large CSV subjet!!" << endl; + if(subjetInfo.nVertices()>0) sjVtxType.push_back(0); + else{ + const TrackIPTagInfo &ipInfo = (*subjetIP)[ith_tagged]; + GlobalPoint pv(ipInfo.primaryVertex()->position().x(),ipInfo.primaryVertex()->position().y(),ipInfo.primaryVertex()->position().z()); + std::vector indices = ipInfo.sortedIndexes(reco::btag::IP3DSig); + const std::vector &ipData = ipInfo.impactParameterData(); + + const TrackIPTagInfo::input_container &tracks = ipInfo.selectedTracks(); + std::vector pseudoVertexTracks; + int nVtxTrks=0; + + const Track * trackPairV0Test[2]; + // cout << "pseudoVtx criteria met for subjet " << k << " (pt=" << subjd->pt() << ", csv=" << (*jetTags_CombinedSvtx)[ith_tagged].second <<") " << endl; + // if((*jetTags_CombinedSvtx)[ith_tagged].second > 0.9) cout << " track debug info for PSVtx: " << endl; + for(std::size_t i=0; i 0.9) cout << " track " << i << " pt: "<< track.pt() << " eta: "<< track.eta() << " phi: " << track.phi() << " displacement: "<< data.ip3d.value() << " dist to jet axis: "<< data.distanceToJetAxis.value() << endl; + nVtxTrks++; + } + // check against all other tracks for V0 track pairs + trackPairV0Test[0] = reco::btag::toTrack(tracks[idx]); + bool ok = true; + for(std::size_t j=0; j= 2 && pseudoVertexV0Filter(pseudoVertexTracks)) sjVtxType.push_back(1); + else sjVtxType.push_back(2); + } + subjetVtxMass.push_back(svm); + subjetVtxPt.push_back(svpt); + subjetVtxEta.push_back(sveta); + subjetVtxPhi.push_back(svphi); + subjetVtxNtrk.push_back(svntrk); + subjetVtxdl.push_back(svdl); + subjetVtxdls.push_back(svdls); + + if(doExtendedFlavorTagging_){ + vector hdr, hpt, heta, hphi, hpdg, pdr, ppt, peta, pphi, ppdg; + for ( reco::JetFlavourInfoMatchingCollection::const_iterator sj = theSubjetFlavourInfos->begin(); sj != theSubjetFlavourInfos->end(); sj++ ) { + if( sqrt(reco::deltaR2(subjd->eta(), subjd->phi(), (*sj).first.get()->eta(), (*sj).first.get()->phi() )) <0.01 ){ + + reco::Jet *aSubjet = const_cast((*sj).first.get()); + + reco::JetFlavourInfo aInfo = (*sj).second; + hadronFlavor.push_back(aInfo.getHadronFlavour()); + partonFlavor.push_back(aInfo.getPartonFlavour()); + const reco::GenParticleRefVector &bHadrons = aInfo.getbHadrons(); + for(reco::GenParticleRefVector::const_iterator it = bHadrons.begin(); it != bHadrons.end(); ++it){ + hdr.push_back(reco::deltaR(aSubjet->eta(), aSubjet->phi(), (*it)->eta(), (*it)->phi())); + hpt.push_back((*it)->pt()); + heta.push_back((*it)->eta()); + hphi.push_back((*it)->phi()); + hpdg.push_back((*it)->pdgId()); + } + const reco::GenParticleRefVector &cHadrons = aInfo.getcHadrons(); + for(reco::GenParticleRefVector::const_iterator it = cHadrons.begin(); it != cHadrons.end(); ++it){ + hdr.push_back(reco::deltaR(aSubjet->eta(), aSubjet->phi(), (*it)->eta(), (*it)->phi())); + hpt.push_back((*it)->pt()); + heta.push_back((*it)->eta()); + hphi.push_back((*it)->phi()); + hpdg.push_back((*it)->pdgId()); + } + const reco::GenParticleRefVector & partons = aInfo.getPartons(); + for(reco::GenParticleRefVector::const_iterator it = partons.begin(); it != partons.end(); ++it){ + pdr.push_back(reco::deltaR( aSubjet->eta(), aSubjet->phi(), (*it)->eta(), (*it)->phi() )); + ppt.push_back((*it)->pt()); + peta.push_back((*it)->eta()); + pphi.push_back((*it)->phi()); + ppdg.push_back((*it)->pdgId()); + } + break; + } // end if reco::deltaR2 subjd (*sj).first.get() + } // end for const_iterator sj = theSubjetFlavourInfos->begin() + hadronDR.push_back(hdr); + hadronPt.push_back(hpt); + hadronEta.push_back(heta); + hadronPhi.push_back(hphi); + hadronPdg.push_back(hpdg); + partonDR.push_back(pdr); + partonPt.push_back(ppt); + partonEta.push_back(peta); + partonPhi.push_back(pphi); + partonPdg.push_back(ppdg); + } + + } + } } else { - sjpt.push_back(-999.); - sjeta.push_back(-999.); - sjphi.push_back(-999.); - sjm.push_back(-999.); + sjpt.push_back(-999.); + sjeta.push_back(-999.); + sjphi.push_back(-999.); + sjm.push_back(-999.); + sjVtxType.push_back(-999.); + csvV2.push_back(-999.); + negCsvV2.push_back(-999); + jptag.push_back(-999.); + if(doExtendedFlavorTagging_){ + hadronFlavor.push_back(-999.); + partonFlavor.push_back(-999.); + } } + //cout << "subjet pt size: "<< sjpt.size() << endl; + //cout << "csv size: "<< csvV1.size() << endl; + //cout << "hadron flavor size: " << hadronFlavor.size() << endl; + //cout << "subjet SV mass size: "<< subjetVtxMass.size() << endl; jets_.jtSubJetPt.push_back(sjpt); jets_.jtSubJetEta.push_back(sjeta); jets_.jtSubJetPhi.push_back(sjphi); jets_.jtSubJetM.push_back(sjm); - + jets_.jtSubJetVtxType.push_back(sjVtxType); + jets_.jtSubJetcsvV2.push_back(csvV2); + jets_.jtSubJetNegCsvV2.push_back(negCsvV2); + jets_.jtSubJetJP.push_back(jptag); + jets_.jtSubJetSvtxm.push_back(subjetVtxMass); + jets_.jtSubJetSvtxpt.push_back(subjetVtxPt); + jets_.jtSubJetSvtxeta.push_back(subjetVtxEta); + jets_.jtSubJetSvtxphi.push_back(subjetVtxPhi); + jets_.jtSubJetSvtxNtrk.push_back(subjetVtxNtrk); + jets_.jtSubJetSvtxdl.push_back(subjetVtxdl); + jets_.jtSubJetSvtxdls.push_back(subjetVtxdls); + if(doExtendedFlavorTagging_){ + jets_.jtSubJetHadronFlavor.push_back(hadronFlavor); + jets_.jtSubJetPartonFlavor.push_back(partonFlavor); + jets_.jtSubJetHadronDR.push_back(hadronDR); + jets_.jtSubJetHadronPt.push_back(hadronPt); + jets_.jtSubJetHadronEta.push_back(hadronEta); + jets_.jtSubJetHadronPhi.push_back(hadronPhi); + jets_.jtSubJetHadronPdg.push_back(hadronPdg); + jets_.jtSubJetPartonDR.push_back(partonDR); + jets_.jtSubJetPartonPt.push_back(partonPt); + jets_.jtSubJetPartonEta.push_back(partonEta); + jets_.jtSubJetPartonPhi.push_back(partonPhi); + jets_.jtSubJetPartonPdg.push_back(partonPdg); + } } diff --git a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_DATA_92X.py b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_DATA_92X.py index 4e0b66973c03e..49b51414dfebe 100644 --- a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_DATA_92X.py +++ b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_DATA_92X.py @@ -49,8 +49,8 @@ process.GlobalTag = GlobalTag(process.GlobalTag, '92X_dataRun2_Prompt_Candidate_forXeXe_v1', '') process.HiForest.GlobalTagLabel = process.GlobalTag.globaltag -from HeavyIonsAnalysis.Configuration.CommonFunctions_cff import overrideJEC_pp5020 -process = overrideJEC_pp5020(process) +#from HeavyIonsAnalysis.Configuration.CommonFunctions_cff import overrideJEC_pp5020 +#process = overrideJEC_pp5020(process) ##################################################################################### # Define tree output @@ -81,7 +81,10 @@ process.load('HeavyIonsAnalysis.JetAnalysis.jets.akPu4PFJetSequence_pp_data_cff') process.load('HeavyIonsAnalysis.JetAnalysis.jets.akCs4PFJetSequence_pp_data_cff') -process.akPu4PFcorr.payload = "AK4PF_offline" +process.akPu4PFcorr.payload = "AK4PF" +process.akCs4PFcorr.payload = "AK4PF" +process.ak4PFcorr.payload = "AK4PF" +process.ak4Calocorr.payload = "AK4PF" process.load('RecoJets.JetProducers.kt4PFJets_cfi') process.load('RecoHI.HiJetAlgos.hiFJRhoProducer') @@ -258,3 +261,11 @@ process.pAna = cms.EndPath(process.skimanalysis) # Customization +process.akPu4PFJetAnalyzer.trackSelection = process.ak4PFSecondaryVertexTagInfos.trackSelection +process.akPu4PFJetAnalyzer.trackPairV0Filter = process.ak4PFSecondaryVertexTagInfos.vertexCuts.v0Filter +process.akCs4PFJetAnalyzer.trackSelection = process.ak4PFSecondaryVertexTagInfos.trackSelection +process.akCs4PFJetAnalyzer.trackPairV0Filter = process.ak4PFSecondaryVertexTagInfos.vertexCuts.v0Filter +process.ak4CaloJetAnalyzer.trackSelection = process.ak4PFSecondaryVertexTagInfos.trackSelection +process.ak4CaloJetAnalyzer.trackPairV0Filter = process.ak4PFSecondaryVertexTagInfos.vertexCuts.v0Filter +process.ak4PFJetAnalyzer.trackSelection = process.ak4PFSecondaryVertexTagInfos.trackSelection +process.ak4PFJetAnalyzer.trackPairV0Filter = process.ak4PFSecondaryVertexTagInfos.vertexCuts.v0Filter diff --git a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_ppMC_92X.py b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_ppMC_92X.py index 53649b4dac27d..1d57847605d52 100644 --- a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_ppMC_92X.py +++ b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_ppMC_92X.py @@ -26,14 +26,14 @@ process.source = cms.Source("PoolSource", duplicateCheckMode = cms.untracked.string("noDuplicateCheck"), fileNames = cms.untracked.vstring( - #"/store/user/gsfs/Hydjet_Quenched_MinBias_Cymbal5Ev8_XeXe_9_2/ppRECO__201711004/171005_005307/0000/step3_MinBias_ppReco_XeXe_RAW2DIGI_L1Reco_RECO_138.root" - "file:step3_MinBias_ppReco_XeXe_RAW2DIGI_L1Reco_RECO_138.root" + "/store/user/gsfs/Hydjet_Quenched_MinBias_Cymbal5Ev8_XeXe_9_2/ppRECO__201711004/171005_005307/0000/step3_MinBias_ppReco_XeXe_RAW2DIGI_L1Reco_RECO_1.root" + #"file:step3_MinBias_ppReco_XeXe_RAW2DIGI_L1Reco_RECO_138.root" ) ) # Number of events we want to process, -1 = all events process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(10)) + input = cms.untracked.int32(100)) ##################################################################################### @@ -63,8 +63,8 @@ #]) # Customization -from HeavyIonsAnalysis.Configuration.CommonFunctions_cff import overrideJEC_pp5020 -process = overrideJEC_pp5020(process) +#from HeavyIonsAnalysis.Configuration.CommonFunctions_cff import overrideJEC_pp5020 +#process = overrideJEC_pp5020(process) ##################################################################################### # Define tree output @@ -84,7 +84,9 @@ ############################# process.load("HeavyIonsAnalysis.JetAnalysis.FullJetSequence_XeXe_mc") -process.akPu4PFcorr.payload = "AK4PF_offline" +process.akPu4PFcorr.payload = "AK4PF" +process.akCs4PFcorr.payload = "AK4PF" +process.ak4PFcorr.payload = "AK4PF" # Use this version for JEC #process.load("HeavyIonsAnalysis.JetAnalysis.FullJetSequence_JECPP") @@ -220,3 +222,15 @@ process.pAna = cms.EndPath(process.skimanalysis) # Customization +process.akPu4PFJetAnalyzer.doExtendedFlavorTagging = cms.untracked.bool(False) +process.akPu4PFJetAnalyzer.trackSelection = process.ak4PFSecondaryVertexTagInfos.trackSelection +process.akPu4PFJetAnalyzer.trackPairV0Filter = process.ak4PFSecondaryVertexTagInfos.vertexCuts.v0Filter +process.akCs4PFJetAnalyzer.doExtendedFlavorTagging = cms.untracked.bool(False) +process.akCs4PFJetAnalyzer.trackSelection = process.ak4PFSecondaryVertexTagInfos.trackSelection +process.akCs4PFJetAnalyzer.trackPairV0Filter = process.ak4PFSecondaryVertexTagInfos.vertexCuts.v0Filter +process.ak4CaloJetAnalyzer.doExtendedFlavorTagging = cms.untracked.bool(False) +process.ak4CaloJetAnalyzer.trackSelection = process.ak4PFSecondaryVertexTagInfos.trackSelection +process.ak4CaloJetAnalyzer.trackPairV0Filter = process.ak4PFSecondaryVertexTagInfos.vertexCuts.v0Filter +process.ak4PFJetAnalyzer.doExtendedFlavorTagging = cms.untracked.bool(False) +process.ak4PFJetAnalyzer.trackSelection = process.ak4PFSecondaryVertexTagInfos.trackSelection +process.ak4PFJetAnalyzer.trackPairV0Filter = process.ak4PFSecondaryVertexTagInfos.vertexCuts.v0Filter From b28be494ca655d11a84e698ba040f7ac1a4a89fb Mon Sep 17 00:00:00 2001 From: kurt Date: Tue, 31 Oct 2017 23:09:32 +0100 Subject: [PATCH 2/3] adding 92X pp MC file --- .../JetAnalysis/test/runForestAOD_pp_MC_92X.py | 15 +++++++++------ 1 file changed, 9 insertions(+), 6 deletions(-) diff --git a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pp_MC_92X.py b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pp_MC_92X.py index 55c03e3545281..c8436280f4858 100644 --- a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pp_MC_92X.py +++ b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_pp_MC_92X.py @@ -26,7 +26,7 @@ process.source = cms.Source("PoolSource", duplicateCheckMode = cms.untracked.string("noDuplicateCheck"), fileNames = cms.untracked.vstring( - "file:/afs/cern.ch/user/k/kjung/eos/cms/store/cmst3/group/hintt/Pbp/8_0_26_patch2/PYQUEN_TTbar/PPb/Events_1.root" + "/store/user/gsfs/Pythia8_Dijet80_pp_CUETP8M1_5020GeV/RECO__201711004/171004_122752/0000/step3_pp_RAW2DIGI_L1Reco_RECO_102.root" ) ) @@ -62,8 +62,8 @@ ]) # Customization -from HeavyIonsAnalysis.Configuration.CommonFunctions_cff import overrideJEC_pp5020 -process = overrideJEC_pp5020(process) +#from HeavyIonsAnalysis.Configuration.CommonFunctions_cff import overrideJEC_pp5020 +#process = overrideJEC_pp5020(process) ##################################################################################### # Define tree output @@ -83,6 +83,7 @@ ############################# process.load("HeavyIonsAnalysis.JetAnalysis.FullJetSequence_nominalPP") +process.ak4PFcorr.payload = "AK4PF" # Use this version for JEC #process.load("HeavyIonsAnalysis.JetAnalysis.FullJetSequence_JECPP") @@ -164,9 +165,9 @@ process.hiEvtAnalyzer * process.HiGenParticleAna* process.jetSequences + - process.egmGsfElectronIDSequence + #Should be added in the path for VID module - process.ggHiNtuplizer + - process.ggHiNtuplizerGED + + #process.egmGsfElectronIDSequence + #Should be added in the path for VID module + #process.ggHiNtuplizer + + #process.ggHiNtuplizerGED + process.pfcandAnalyzer + process.HiForest + process.trackSequencesPP + @@ -215,3 +216,5 @@ process.pAna = cms.EndPath(process.skimanalysis) # Customization +process.ak4PFJetAnalyzer.trackSelection = process.ak4PFSecondaryVertexTagInfos.trackSelection +process.ak4PFJetAnalyzer.trackPairV0Filter = process.ak4PFSecondaryVertexTagInfos.vertexCuts.v0Filter From e715213e0d709777ed9c50dd7b862c8da08c8ef7 Mon Sep 17 00:00:00 2001 From: kurt Date: Tue, 31 Oct 2017 23:13:11 +0100 Subject: [PATCH 3/3] updating XeXe and PbPb cfgs for flavor tagging --- .../test/runForestAOD_PbPb_MB_92X.py | 24 ++++++++++++++++--- .../test/runForestAOD_XeXe_DATA_92X.py | 2 +- 2 files changed, 22 insertions(+), 4 deletions(-) diff --git a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_PbPb_MB_92X.py b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_PbPb_MB_92X.py index 3725460da2f98..07989e0299a52 100644 --- a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_PbPb_MB_92X.py +++ b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_PbPb_MB_92X.py @@ -26,14 +26,14 @@ process.source = cms.Source("PoolSource", duplicateCheckMode = cms.untracked.string("noDuplicateCheck"), fileNames = cms.untracked.vstring( - "file:/afs/cern.ch/user/m/mverweij/merge/xexe/step3_RAW2DIGI_L1Reco_RECO_1.root", + "file:/afs/cern.ch/user/k/kjung/gitForest/CMSSW_9_3_1/src/HeavyIonsAnalysis/JetAnalysis/test/step3_MinBias_XeXe_RAW2DIGI_L1Reco_RECO_MERGED.root" #samples/PbPb_MC_RECODEBUG.root" ) ) # Number of events we want to process, -1 = all events process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(10) + input = cms.untracked.int32(2) ) @@ -48,12 +48,26 @@ process.load('FWCore.MessageService.MessageLogger_cfi') from Configuration.AlCa.GlobalTag_condDBv2 import GlobalTag -process.GlobalTag = GlobalTag(process.GlobalTag, '91X_mcRun2_HeavyIon_v3', '') #for now track GT manually, since centrality tables updated ex post facto +process.GlobalTag = GlobalTag(process.GlobalTag, '93X_mcRun2_HeavyIon_v1', '') #for now track GT manually, since centrality tables updated ex post facto process.HiForest.GlobalTagLabel = process.GlobalTag.globaltag from HeavyIonsAnalysis.Configuration.CommonFunctions_cff import overrideJEC_PbPb5020 process = overrideJEC_PbPb5020(process) +from CondCore.DBCommon.CondDBSetup_cfi import * +process.akPu4PFJECoverride = cms.ESSource("PoolDBESSource", CondDBSetup, + connect = cms.string('sqlite_file:testakPu4PF_JetCorrections_newFunction.db'), + toGet = cms.VPSet( # overide Global Tag use EcalTBWeights_EBEE_offline + cms.PSet( + record = cms.string('JetCorrectionsRecord'), + tag = cms.string('JetCorrectorParametersCollection_PbPb_92X_5020GeV_v0_AK4PF_offline'), + label = cms.untracked.string('AKPu4PF_offline') + ) + ) + ) +process.es_prefer_jpTagConds = cms.ESPrefer("PoolDBESSource","akPu4PFJECoverride") + + process.load("RecoHI.HiCentralityAlgos.CentralityBin_cfi") ##################################################################################### @@ -85,6 +99,10 @@ #rho analyzer process.load('HeavyIonsAnalysis.JetAnalysis.hiFJRhoAnalyzer_cff') +process.akPu4PFcorr.levels = cms.vstring('L2Relative','L3Absolute') +#process.akCs4PFcorr.levels = cms.vstring('L1FastJet','L2Relative','L3Absolute','L2L3Residual') +#process.akPu4PFcorr.rho = cms.InputTag('hiFJRhoProducer','mapToRho') + #################################################################################### ############################# diff --git a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_DATA_92X.py b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_DATA_92X.py index dbfff347c24a5..8243149136f11 100644 --- a/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_DATA_92X.py +++ b/HeavyIonsAnalysis/JetAnalysis/test/runForestAOD_XeXe_DATA_92X.py @@ -32,7 +32,7 @@ # Number of events we want to process, -1 = all events process.maxEvents = cms.untracked.PSet( - input = cms.untracked.int32(10)) + input = cms.untracked.int32(100)) #####################################################################################