diff --git a/DisplacedAdaptiveVertexFinder/python/__init__.py b/DisplacedAdaptiveVertexFinder/python/__init__.py new file mode 100644 index 0000000..513c5d6 --- /dev/null +++ b/DisplacedAdaptiveVertexFinder/python/__init__.py @@ -0,0 +1,3 @@ +#Automatically created by SCRAM +import os +__path__.append(os.path.dirname(os.path.abspath(__file__).rsplit('/RecoVertex/AdaptiveVertexFinder/',1)[0])+'/cfipython/slc6_amd64_gcc481/RecoVertex/AdaptiveVertexFinder') diff --git a/DisplacedAdaptiveVertexFinder/python/displacedInclusiveVertexFinder_cfi.py b/DisplacedAdaptiveVertexFinder/python/displacedInclusiveVertexFinder_cfi.py new file mode 100644 index 0000000..ad82e28 --- /dev/null +++ b/DisplacedAdaptiveVertexFinder/python/displacedInclusiveVertexFinder_cfi.py @@ -0,0 +1,41 @@ +import FWCore.ParameterSet.Config as cms + +displacedInclusiveVertexFinder = cms.EDProducer("InclusiveVertexFinder", + beamSpot = cms.InputTag("offlineBeamSpot"), + primaryVertices = cms.InputTag("offlineSlimmedPrimaryVertices"), + #tracks = cms.InputTag("displacedAssocToTracks","displacedAssocToTracks","ANA"), + tracks = cms.InputTag("unpackedTracksAndVertices"), + minHits = cms.uint32(6), #old 8 -> 0 AOD produciton has problems with nhits + maximumLongitudinalImpactParameter = cms.double(99999), #old .3 -> infty + minPt = cms.double(0.4), #old .8 -> 1 + maxNTracks = cms.uint32(100), #old 30 -> 100 + + clusterizer = cms.PSet( + seedMax3DIPSignificance = cms.double(9999.), + seedMax3DIPValue = cms.double(9999.), + seedMin3DIPSignificance = cms.double(-9999), + seedMin3DIPValue = cms.double(-9999), + clusterMaxDistance = cms.double(0.4), #500um #old .05 -> 1 + clusterMaxSignificance = cms.double(4.5), #4.5 sigma #old 4.5 ---> infty + distanceRatio = cms.double(20), # was cluster scale = 1 / density factor =0.05 + clusterMinAngleCosine = cms.double(0.5), # only forward decays #old accept backward decays (unboosted topologies) .5 -> -9999 + ), + + vertexMinAngleCosine = cms.double(0.95), # scalar prod direction of tracks and flight dir #old accept backward decays .95 -> .6 me 0 + vertexMinDLen2DSig = cms.double(2.5), #2.5 sigma + vertexMinDLenSig = cms.double(0.5), #0.5 sigma + fitterSigmacut = cms.double(3), + fitterTini = cms.double(256), + fitterRatio = cms.double(0.25), + useDirectVertexFitter = cms.bool(True), + useVertexReco = cms.bool(True), + vertexReco = cms.PSet( + finder = cms.string('avr'), + primcut = cms.double(1.0), + seccut = cms.double(3), + smoothing = cms.bool(True) + ) + +) + + diff --git a/DisplacedAdaptiveVertexFinder/python/displacedInclusiveVertexing_cff.py b/DisplacedAdaptiveVertexFinder/python/displacedInclusiveVertexing_cff.py new file mode 100644 index 0000000..5679a63 --- /dev/null +++ b/DisplacedAdaptiveVertexFinder/python/displacedInclusiveVertexing_cff.py @@ -0,0 +1,17 @@ +import FWCore.ParameterSet.Config as cms + +#from HNL.DisplacedAdaptiveVertexFinder.unpackedTracksAndVertices_cfi import * +from PhysicsTools.PatAlgos.slimming.unpackedTracksAndVertices_cfi import * +from HNL.DisplacedAdaptiveVertexFinder.displacedInclusiveVertexFinder_cfi import * +from HNL.DisplacedAdaptiveVertexFinder.displacedVertexMerger_cfi import * +from HNL.DisplacedAdaptiveVertexFinder.displacedTrackVertexArbitrator_cfi import * +from HNL.DisplacedSVAssociator.displacedSVAssociationIVF_cfi import * + +displacedInclusiveSecondaryVertices = displacedVertexMerger.clone() +displacedInclusiveSecondaryVertices.secondaryVertices = cms.InputTag("displacedTrackVertexArbitrator") +displacedInclusiveSecondaryVertices.maxFraction = 0.05 #0.05 #old .2 -> .05 +displacedInclusiveSecondaryVertices.minSignificance = 10 + + +displacedInclusiveVertexing = cms.Sequence(unpackedTracksAndVertices * displacedInclusiveVertexFinder * displacedVertexMerger * displacedTrackVertexArbitrator * displacedInclusiveSecondaryVertices ) + diff --git a/DisplacedAdaptiveVertexFinder/python/displacedTrackVertexArbitrator_cfi.py b/DisplacedAdaptiveVertexFinder/python/displacedTrackVertexArbitrator_cfi.py new file mode 100644 index 0000000..5b9ed1e --- /dev/null +++ b/DisplacedAdaptiveVertexFinder/python/displacedTrackVertexArbitrator_cfi.py @@ -0,0 +1,21 @@ +import FWCore.ParameterSet.Config as cms + +displacedTrackVertexArbitrator = cms.EDProducer("TrackVertexArbitrator", + beamSpot = cms.InputTag("offlineBeamSpot"), + primaryVertices = cms.InputTag("offlineSlimmedPrimaryVertices"), + #tracks = cms.InputTag("displacedAssocToTracks","displacedAssocToTracks","ANA"), + tracks = cms.InputTag("unpackedTracksAndVertices"), + secondaryVertices = cms.InputTag("displacedVertexMerger"), + dLenFraction = cms.double(0.333), #old .333 -> .2 + dRCut = cms.double(1), # old .4 -> 1 me 3 + distCut = cms.double(0.1), #old .04 -> .1 + sigCut = cms.double(5), #old 5->10 + fitterSigmacut = cms.double(3), + fitterTini = cms.double(256), + fitterRatio = cms.double(0.25), + trackMinLayers = cms.int32(4), #old 4-> 0 + trackMinPt = cms.double(.4), + trackMinPixels = cms.int32(0) #old 1 -> 0 +) + + diff --git a/DisplacedAdaptiveVertexFinder/python/displacedVertexMerger_cfi.py b/DisplacedAdaptiveVertexFinder/python/displacedVertexMerger_cfi.py new file mode 100644 index 0000000..01b7066 --- /dev/null +++ b/DisplacedAdaptiveVertexFinder/python/displacedVertexMerger_cfi.py @@ -0,0 +1,8 @@ +import FWCore.ParameterSet.Config as cms + +displacedVertexMerger = cms.EDProducer("VertexMerger", + secondaryVertices = cms.InputTag("displacedInclusiveVertexFinder"), + maxFraction = cms.double(0.7), #old .7 -> .5 + minSignificance = cms.double(2)) + + diff --git a/DisplacedAdaptiveVertexFinder/python/unpackedTracksAndVertices_cfi.py b/DisplacedAdaptiveVertexFinder/python/unpackedTracksAndVertices_cfi.py new file mode 100644 index 0000000..feeb13c --- /dev/null +++ b/DisplacedAdaptiveVertexFinder/python/unpackedTracksAndVertices_cfi.py @@ -0,0 +1,7 @@ +import FWCore.ParameterSet.Config as cms +unpackedTracksAndVertices = cms.EDProducer('PATTrackAndVertexUnpacker', + slimmedVertices = cms.InputTag("offlineSlimmedPrimaryVertices"), + slimmedSecondaryVertices = cms.InputTag("slimmedSecondaryVertices"), + additionalTracks= cms.InputTag("lostTracks"), + packedCandidates = cms.InputTag("packedPFCandidates") + ) diff --git a/DisplacedSVAssociator/interface/VertexAssociation.h b/DisplacedSVAssociator/interface/VertexAssociation.h new file mode 100644 index 0000000..c332164 --- /dev/null +++ b/DisplacedSVAssociator/interface/VertexAssociation.h @@ -0,0 +1,98 @@ +// Muon +#include "DataFormats/PatCandidates/interface/Muon.h" + +// tracks +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/TrackReco/interface/TrackFwd.h" + +// vertex +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/BTauReco/interface/VertexTypes.h" + +// kinematics +#include "DataFormats/Math/interface/deltaR.h" + +class VertexAssociation { + public: + VertexAssociation(std::string nameString, const reco::Vertex & PV, const int& debug_) { + primaryVertex = PV; + name = nameString; + debug = debug_; + } + + // fillers + //void addMuon( const pat::MuonCollection & muon); + void addVertex( const reco::Vertex & vertex); + void setPrimaryVertex(const reco::Vertex & pv){ primaryVertex = pv; }; + + // accessors + int getNVertices(){ return vertexCollection.size(); } + //int getNMuon() { return MuonCollection.size(); } + float getBestVertexScore() { return bestVertexScore; } + + std::string getName() { return name; } + // score + const std::pair getBestVertex(const pat::Electron & , const std::string&); + float getVertexScore(const pat::Electron & , const reco::Vertex & , const std::string &); + + private: + reco::Vertex primaryVertex; + //pat::MuonCollection muonCollection; + reco::VertexCollection vertexCollection; + reco::Vertex bestVertex; + float bestVertexScore; + std::string name; + int debug; +}; +//void VertexAssociation::addMuon(const pat::MuonCollection & muon) { muonCollection.push_back(muon); } +void VertexAssociation::addVertex(const reco::Vertex & vertex) { vertexCollection.push_back(vertex); } +float VertexAssociation::getVertexScore(const pat::Electron & electron , const reco::Vertex & vertex, const std::string & algo) { + float score = 0; + double ele_pt = electron.pt(); + double ele_eta = electron.eta(); + double ele_phi = electron.phi(); + + if (debug > 3) std::cout << "[DEBUG 3] [SVA] electron pT=" << ele_pt << " electron Eta: " << ele_eta <<" electron Phi "<< ele_phi << std::endl; + if (debug > 3) std::cout << "[DEBUG 3] [SVA] iterating over tracks" << std::endl; + reco::Vertex::trackRef_iterator tt = vertex.tracks_begin(); + for(; tt != vertex.tracks_end(); ++tt) { + float eta = (*tt)->eta(); + float phi = (*tt)->phi(); + if (debug > 4) std::cout << "[DEBUG 4] [SVA] Track Eta=" << eta << " Phi: " << phi << std::endl; + float dR = reco::deltaR(ele_eta, ele_phi, eta, phi); + if (dR > 1) // was 0.4 + continue; + else + score += 1.0 / dR; + } + if (debug > 3) std::cout << "[DEBUG 3] [SVA] Vertex Matching Score: " << score << std::endl; + return score; +} + +const std::pair VertexAssociation::getBestVertex(const pat::Electron & electron, const std::string & algo) { + // keep track of the best vertex + float bestScore = -1; + reco::Vertex bestVertex; + + if (debug > 1) std::cout << "[DEBUG 1] [SVA] iterating over Vertices" << std::endl; + reco::VertexCollection::const_iterator ss = vertexCollection.begin(); + for(; ss != vertexCollection.end(); ++ss) { + float score = getVertexScore(electron, *ss, algo); + + if (score > bestScore) { // was bestScore + bestScore = score; + bestVertex = *ss; + } + } + + if (bestScore > 0) { + const std::pair vertexPair(bestVertex, bestScore); + return vertexPair; // we find at least one vertex with a track within dR = 1.0 + } + + else { + const std::pair vertexPair(primaryVertex, 0); + return vertexPair; + } +} diff --git a/DisplacedSVAssociator/plugins/BuildFile.xml b/DisplacedSVAssociator/plugins/BuildFile.xml new file mode 100644 index 0000000..337f3c4 --- /dev/null +++ b/DisplacedSVAssociator/plugins/BuildFile.xml @@ -0,0 +1,13 @@ + + + + + + + + + + + + + diff --git a/DisplacedSVAssociator/plugins/DisplacedSVAssociator.cc b/DisplacedSVAssociator/plugins/DisplacedSVAssociator.cc new file mode 100644 index 0000000..93a5010 --- /dev/null +++ b/DisplacedSVAssociator/plugins/DisplacedSVAssociator.cc @@ -0,0 +1,291 @@ +// -*- C++ -*- +// +// Package: TEST/DisplacedSVAssociator +// Class: DisplacedSVAssociator +// +/**\class DisplacedSVAssociator DisplacedSVAssociator.cc TEST/DisplacedSVAssociator/plugins/DisplacedSVAssociator.cc + + Description: [one line class summary] + + Implementation: + [Notes on implementation] +*/ +// +// Original Author: Mohamed Darwish +// Created: Sun, 18 Jun 2017 15:14:26 GMT +// +// + + + +// system include files +#include + +// user include files +#include "FWCore/Framework/interface/Frameworkfwd.h" +#include "FWCore/Framework/interface/stream/EDProducer.h" +#include "DataFormats/TrackReco/interface/Track.h" +#include "DataFormats/VertexReco/interface/Vertex.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/Math/interface/deltaR.h" + +#include "FWCore/Framework/interface/Event.h" +#include "FWCore/Framework/interface/MakerMacros.h" + + +// kinematics + +#include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "FWCore/Utilities/interface/StreamID.h" +#include "HNL/DisplacedSVAssociator/interface/VertexAssociation.h" + + +// +// class declaration +// + +class DisplacedSVAssociator : public edm::stream::EDProducer<> { + public: + explicit DisplacedSVAssociator(const edm::ParameterSet&); + ~DisplacedSVAssociator(); + + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + typedef std::vector VertexScores; + + private: + virtual void beginStream(edm::StreamID) override; + virtual void produce(edm::Event&, const edm::EventSetup&) override; + virtual void endStream() override; + edm::EDGetTokenT muonToken_; + edm::EDGetTokenT tag_secondaryVertices_; + edm::EDGetTokenT tag_primaryVertices_; + + std::string algoName_; + std::string outputLabel_; + float muonPtCut_; + int debug_; + + //virtual void beginRun(edm::Run const&, edm::EventSetup const&) override; + //virtual void endRun(edm::Run const&, edm::EventSetup const&) override; + //virtual void beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override; + //virtual void endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) override; + + // ----------member data --------------------------- +}; + +// +// constants, enums and typedefs +// + + +// +// static data member definitions +// + +// +// constructors and destructor +// +DisplacedSVAssociator::DisplacedSVAssociator(const edm::ParameterSet& iConfig): + muonToken_(consumes(iConfig.getUntrackedParameter("muonTag"))), + tag_secondaryVertices_(consumes(iConfig.getUntrackedParameter("secondaryVertices"))), + tag_primaryVertices_ (consumes(iConfig.getUntrackedParameter("primaryVertices"))) + +{ + algoName_ = iConfig.getUntrackedParameter("algoName"); + muonPtCut_ = iConfig.getUntrackedParameter("muonPtCut"); + debug_ = iConfig.getUntrackedParameter("debug"); + outputLabel_ = iConfig.getUntrackedParameter("outputLabel"); + + produces >(outputLabel_); + //register your products +/* Examples + produces(); + + //if do put with a label + produces("label"); + + //if you want to put into the Run + produces(); +*/ + //now do what ever other initialization is needed + +} + + +DisplacedSVAssociator::~DisplacedSVAssociator() +{ + + // do anything here that needs to be done at destruction time + // (e.g. close files, deallocate resources etc.) + +} + + +// +// member functions +// + +// ------------ method called to produce the data ------------ +void +DisplacedSVAssociator::produce(edm::Event& iEvent, const edm::EventSetup& iSetup) +{ + using namespace edm; + //std::auto_ptr > VertexAssociationCollection(new std::vector); + std::unique_ptr > VertexAssociationCollection(new std::vector); + edm::Handle mu; + edm::HandlesecondaryVertices; + edm::HandleprimaryVertices; + + //get the products from the tags + iEvent.getByToken(muonToken_, mu); + iEvent.getByToken(tag_secondaryVertices_, secondaryVertices); + iEvent.getByToken(tag_primaryVertices_, primaryVertices); + + const pat::MuonCollection& muon = *(mu.product()); + const reco::VertexCollection & sv = *(secondaryVertices.product()); + const reco::VertexCollection & pv = *(primaryVertices.product()); + + //use the first primary vertex for now + const reco::Vertex primaryVertex = pv[0]; + if (debug_ > 1) std::cout << "[DEBUG] [SV Associator] Loop over muons " << std::endl; + float Muon_pt= 999; + float Muon_eta = 999; + float Muon_phi = 999; + + // loop over each jet + // reco::PFJetCollection::const_iterator jj = jet.begin(); + //for(; jj != jet.end() ; ++jj){ + pat::MuonCollection::const_iterator Mu = muon.begin(); + for(; Mu != muon.end() ; ++Mu){ + reco::TrackRef tunePTrack = Mu->muonBestTrack(); + if(tunePTrack->pt() < Muon_pt){ + Muon_pt = tunePTrack->pt(); + Muon_eta = tunePTrack->eta(); + Muon_phi = tunePTrack->phi(); + } + + // minimum Muon pt + if (Muon_pt < muonPtCut_) continue; + + //start building the vertex scores + VertexScores vertexScores; + float bestScore = -1; + reco::Vertex bestVertex; + + //loop over each of the vertices to be scored + if (debug_ > 1) std::cout << "[DEBUG] [SV Associator] Loop over Vertices " << std::endl; + reco::VertexCollection::const_iterator ss = sv.begin(); + size_t iss = 0; + for(; ss != sv.end(); ++ss, ++iss) { + float score = 0; + + //loop over each track in the vertex + if (debug_ > 1) std::cout << "[DEBUG] [SV Associator] Loop over Tracks " << std::endl; + std::vector::const_iterator tt = ss->tracks_begin(); + for(; tt != ss->tracks_end(); ++tt) { + float track_outerEta = (*tt)->outerEta(); + float track_outerPhi = (*tt)->outerPhi(); + + float dR = reco::deltaR(Muon_eta, Muon_phi, track_outerEta, track_outerPhi); + + if (dR > 1.0) continue; + else score += 1.0 / dR; + } + + if (debug_ > 1) std::cout << "[DEBUG] [SV Associator] Inserting Scores " << std::endl; + + //parse a vertex reference using the handle and size + //std::pair scorePair(*ss, score); + vertexScores.push_back(score); + + if (score > bestScore) { + bestScore = score; + bestVertex = *ss; + } + } + + if (debug_ > 1) std::cout << "[DEBUG] [SV Associator] Building Association" << std::endl; + // put together the scores with the best vertex + //VertexAssociation association(vertexScores, *jj, bestVertex, algoName_); + (*VertexAssociationCollection).push_back(vertexScores); + //(*VertexAssociationCollection).push_back(bestVertex); + + } + outputLabel_.clear(); + iEvent.put(std::move(VertexAssociationCollection), outputLabel_); +/* This is an event example + //Read 'ExampleData' from the Event + Handle pIn; + iEvent.getByLabel("example",pIn); + + //Use the ExampleData to create an ExampleData2 which + // is put into the Event + std::unique_ptr pOut(new ExampleData2(*pIn)); + iEvent.put(std::move(pOut)); +*/ + +/* this is an EventSetup example + //Read SetupData from the SetupRecord in the EventSetup + ESHandle pSetup; + iSetup.get().get(pSetup); +*/ + +} + +// ------------ method called once each stream before processing any runs, lumis or events ------------ +void +DisplacedSVAssociator::beginStream(edm::StreamID) +{ +} + +// ------------ method called once each stream after processing all runs, lumis and events ------------ +void +DisplacedSVAssociator::endStream() { +} + +// ------------ method called when starting to processes a run ------------ +/* +void +DisplacedSVAssociator::beginRun(edm::Run const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method called when ending the processing of a run ------------ +/* +void +DisplacedSVAssociator::endRun(edm::Run const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method called when starting to processes a luminosity block ------------ +/* +void +DisplacedSVAssociator::beginLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method called when ending the processing of a luminosity block ------------ +/* +void +DisplacedSVAssociator::endLuminosityBlock(edm::LuminosityBlock const&, edm::EventSetup const&) +{ +} +*/ + +// ------------ method fills 'descriptions' with the allowed parameters for the module ------------ +void +DisplacedSVAssociator::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { + //The following says we do not know what parameters are allowed so do no validation + // Please change this to state exactly what you do use, even if it is no parameters + edm::ParameterSetDescription desc; + desc.setUnknown(); + descriptions.addDefault(desc); +} + +//define this as a plug-in +DEFINE_FWK_MODULE(DisplacedSVAssociator); diff --git a/DisplacedSVAssociator/python/CfiFile_cfi.py b/DisplacedSVAssociator/python/CfiFile_cfi.py new file mode 100644 index 0000000..228add3 --- /dev/null +++ b/DisplacedSVAssociator/python/CfiFile_cfi.py @@ -0,0 +1,4 @@ +import FWCore.ParameterSet.Config as cms + +demo = cms.EDProducer('DisplacedSVAssociator' +) diff --git a/DisplacedSVAssociator/python/ConfFile_cfg.py b/DisplacedSVAssociator/python/ConfFile_cfg.py new file mode 100644 index 0000000..5c73091 --- /dev/null +++ b/DisplacedSVAssociator/python/ConfFile_cfg.py @@ -0,0 +1,26 @@ +import FWCore.ParameterSet.Config as cms + +process = cms.Process("OWNPARTICLES") + +process.load("FWCore.MessageService.MessageLogger_cfi") + +process.maxEvents = cms.untracked.PSet( input = cms.untracked.int32(-1) ) + +process.source = cms.Source("PoolSource", + # replace 'myfile.root' with the source file you want to use + fileNames = cms.untracked.vstring( + 'file:myfile.root' + ) +) + +process.myProducerLabel = cms.EDProducer('DisplacedSVAssociator' +) + +process.out = cms.OutputModule("PoolOutputModule", + fileName = cms.untracked.string('myOutputFile.root') +) + + +process.p = cms.Path(process.myProducerLabel) + +process.e = cms.EndPath(process.out) diff --git a/DisplacedSVAssociator/python/__init__.py b/DisplacedSVAssociator/python/__init__.py new file mode 100644 index 0000000..5a7e4a2 --- /dev/null +++ b/DisplacedSVAssociator/python/__init__.py @@ -0,0 +1,3 @@ +#Automatically created by SCRAM +import os +__path__.append(os.path.dirname(os.path.abspath(__file__).rsplit('/TEST/DisplacedJetSVAssociator/',1)[0])+'/cfipython/slc6_amd64_gcc530/TEST/DisplacedJetSVAssociator') diff --git a/DisplacedSVAssociator/python/displacedSVAssociationIVF_cfi.py b/DisplacedSVAssociator/python/displacedSVAssociationIVF_cfi.py new file mode 100644 index 0000000..73b6baa --- /dev/null +++ b/DisplacedSVAssociator/python/displacedSVAssociationIVF_cfi.py @@ -0,0 +1,14 @@ +import FWCore.ParameterSet.Config as cms + +displacedSVAssociationIVF = cms.EDProducer("DisplacedSVAssociator", + #pfJets = cms.untracked.InputTag("ak4PFJetsCHS","",""), + #pfJets = cms.untracked.InputTag("ak8PFJetsCHS","",""), + muonTag = cms.untracked.InputTag("slimmedMuons","",""), + secondaryVertices = cms.untracked.InputTag("displacedInclusiveSecondaryVertices"), + primaryVertices = cms.untracked.InputTag("offlinePrimaryVertices"), + outputLabel = cms.untracked.string('displacedIVFAssoc'), + muonPtCut = cms.untracked.double(5.0), + algoName = cms.untracked.string("oneOverDeltaR"), + debug = cms.untracked.int32(1) + #jets = cms.untracked.InputTag("ak4PFJetsCHSL1FastL2L3","","") +) diff --git a/HeavyNeutralLeptonAnalysis/BuildFile.xml b/HeavyNeutralLeptonAnalysis/BuildFile.xml new file mode 100644 index 0000000..2340719 --- /dev/null +++ b/HeavyNeutralLeptonAnalysis/BuildFile.xml @@ -0,0 +1,20 @@ + + + + + + + + + + + + + + + + + + + + diff --git a/HeavyNeutralLeptonAnalysis/interface/BigNtuple.h b/HeavyNeutralLeptonAnalysis/interface/BigNtuple.h new file mode 100644 index 0000000..1dc03e3 --- /dev/null +++ b/HeavyNeutralLeptonAnalysis/interface/BigNtuple.h @@ -0,0 +1,561 @@ +#ifndef HNL_HeavyNeutralLeptonAnalysis_BigNtuple +#define HNL_HeavyNeutralLeptonAnalysis_BigNtuple +/* + Class: BigNtuple + Simple interface class to hide all the ROOT I/O from the plugin and make it more readable +*/ + +#include + +#include "TTree.h" + +#include "DataFormats/PatCandidates/interface/Muon.h" +#include "DataFormats/PatCandidates/interface/Electron.h" +#include "DataFormats/PatCandidates/interface/Jet.h" +#include "DataFormats/PatCandidates/interface/MET.h" +#include "FWCore/Common/interface/TriggerNames.h" +#include "DataFormats/Common/interface/TriggerResults.h" +#include "DataFormats/Math/interface/deltaR.h" +#include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "FWCore/Framework/interface/Event.h" +#include "DataFormats/EcalDetId/interface/EcalSubdetector.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/Records/interface/CaloTopologyRecord.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "DataFormats/EcalDetId/interface/EcalSubdetector.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHcalIsolation.h" +#include "DataFormats/PatCandidates/interface/GenericParticle.h" +#include "DataFormats/Common/interface/Handle.h" + +class BigNtuple { +public: + BigNtuple(){} //default, empty constructor + + void set_evtInfo(TTree* tree); + void fill_evtInfo(const edm::EventID& id); + + void set_pv_genInfo(TTree* tree); + void fill_pv_genInfo(const reco::GenParticle prt , const reco::Candidate* mom); + + void set_sv_genInfo(TTree* tree); + void fill_sv_genInfo(const reco::GenParticle prt , const reco::Candidate* mom); + + void set_pvInfo(TTree* tree); + void fill_pvInfo(const reco::VertexCollection& pvs); + + void set_trigInfo(TTree* tree); + void fill_trigInfo(const edm::TriggerResults& triggerResults, const edm::TriggerNames& trigNames); + + void set_pileupInfo(TTree* tree); + void fill_pileupInfo( float npT, float npIT, float pu_Weight, float pu_WeightUp, float pu_WeightDown); + + void set_svInfo(TTree* tree); + void fill_sv_mu_Info(const reco::Vertex& bestVertex, const reco::Vertex& pv, double match); + void fill_sv_ele_Info(const reco::Vertex& bestVertex, const reco::Vertex& pv, double match , double score); + + void set_muInfo(TTree* tree); + void fill_muInfo(const pat::Muon& mu, const reco::Vertex& pv, double Rho , double match1 , double match2 ); + + void set_jetInfo(TTree* tree); + void fill_jetInfo(const pat::Jet& jet); + + void set_metInfo(TTree* tree); + void fill_metInfo(const pat::MET& met); + + void set_eleInfo(TTree* tree); + void fill_eleInfo(const pat::Electron& ele_ , const reco::Vertex& pv, double Rho, double match1, double match2 , std::auto_ptr recHitEcal); + + + void set_eleIDInfo(TTree* tree); + void fill_eleIDInfo(float ele_mva , bool ele_veto , bool ele_loose , bool ele_medium , bool ele_tight); + + void set_bjetInfo(TTree* tree); + void fill_bjetInfo(const pat::Jet& jet, const std::string & bDiscr, int flavor); + + void reset() { + BigNtuple dummy; //create a new one + *this = dummy; //use assignment to reset + } + +private: + unsigned int lumi_ = 0; + unsigned int run_ = 0; + unsigned long long evt_ = 0; + + float npT_ = 0; + float npIT_ = 0; + float pu_Weight_ = 0; + float pu_WeightUp_ = 0; + float pu_WeightDown_ = 0; + + // primary vertex infos -- they shouldn't be vector + float pvX_ = -1000; + float pvY_ = -1000; + float pvZ_ = -1000; + float pvXErr_ = -1000; + float pvYErr_ = -1000; + float pvZErr_ = -1000; + float pvMass_ = -1000; + float pvLxy_ = -1000; + float pvLxyz_ = -1000; + float pvLxySig_ = -1000; + float pvLxyzSig_ = -1000; + float pvChi2_ = -1000; + int pvNTrack_ = -1000; + float pvSumPtSq_ = -1000; + int numberPV_ = -1000; + + //gen infos mu @ pv + std::vector mu_gen_PID1_; + std::vector mu_gen_Status1_; + std::vector mu_gen_Charge1_; + std::vector mu_gen_Pt1_; + std::vector mu_gen_Eta1_; + std::vector mu_gen_Phi1_; + std::vector mu_gen_VX1_; + std::vector mu_gen_VY1_; + std::vector mu_gen_VZ1_; + std::vector mu_gen_Lxy1_; + std::vector mu_gen_Lxyz1_; + std::vector mu_gen_MomPID1_; + std::vector mu_gen_MomStatus1_; + std::vector mu_gen_MomMass1_; + std::vector mu_gen_MomCharge1_; + std::vector mu_gen_MomPt1_; + std::vector mu_gen_MomEta1_; + std::vector mu_gen_MomPhi1_; + std::vector mu_gen_MomBeta1_; + std::vector mu_gen_MomGamma1_; + std::vector mu_gen_MomLxyz1_; + std::vector mu_gen_MomLz1_; + std::vector mu_gen_MomLxy1_; + + //gen Info mu @ sv + std::vector mu_gen_PID2_; + std::vector mu_gen_Status2_; + std::vector mu_gen_Charge2_; + std::vector mu_gen_Pt2_; + std::vector mu_gen_Eta2_; + std::vector mu_gen_Phi2_; + std::vector mu_gen_VX2_; + std::vector mu_gen_VY2_; + std::vector mu_gen_VZ2_; + std::vector mu_gen_Lxy2_; + std::vector mu_gen_Lxyz2_; + std::vector mu_gen_MomPID2_; + std::vector mu_gen_MomMass2_; + std::vector mu_gen_MomCharge2_; + std::vector mu_gen_MomStatus2_; + std::vector mu_gen_MomPt2_; + std::vector mu_gen_MomEta2_; + std::vector mu_gen_MomPhi2_; + std::vector mu_gen_MomBeta2_; + std::vector mu_gen_MomGamma2_; + std::vector mu_gen_MomLxyz2_; + std::vector mu_gen_MomLz2_; + std::vector mu_gen_MomLxy2_; + std::vector mu_gen_MomCTau02_; + + //gen infos ele @ pv + std::vector ele_gen_PID1_; + std::vector ele_gen_Status1_; + std::vector ele_gen_Charge1_; + std::vector ele_gen_Pt1_; + std::vector ele_gen_Eta1_; + std::vector ele_gen_Phi1_; + std::vector ele_gen_VX1_; + std::vector ele_gen_VY1_; + std::vector ele_gen_VZ1_; + std::vector ele_gen_Lxy1_; + std::vector ele_gen_Lxyz1_; + std::vector ele_gen_MomPID1_; + std::vector ele_gen_MomStatus1_; + std::vector ele_gen_MomMass1_; + std::vector ele_gen_MomCharge1_; + std::vector ele_gen_MomPt1_; + std::vector ele_gen_MomEta1_; + std::vector ele_gen_MomPhi1_; + std::vector ele_gen_MomBeta1_; + std::vector ele_gen_MomGamma1_; + std::vector ele_gen_MomLxyz1_; + std::vector ele_gen_MomLz1_; + std::vector ele_gen_MomLxy1_; + + //gen Info ele @ sv + std::vector ele_gen_PID2_; + std::vector ele_gen_Status2_; + std::vector ele_gen_Charge2_; + std::vector ele_gen_Pt2_; + std::vector ele_gen_Eta2_; + std::vector ele_gen_Phi2_; + std::vector ele_gen_VX2_; + std::vector ele_gen_VY2_; + std::vector ele_gen_VZ2_; + std::vector ele_gen_Lxy2_; + std::vector ele_gen_Lxyz2_; + std::vector ele_gen_MomPID2_; + std::vector ele_gen_MomMass2_; + std::vector ele_gen_MomCharge2_; + std::vector ele_gen_MomStatus2_; + std::vector ele_gen_MomPt2_; + std::vector ele_gen_MomEta2_; + std::vector ele_gen_MomPhi2_; + std::vector ele_gen_MomBeta2_; + std::vector ele_gen_MomGamma2_; + std::vector ele_gen_MomLxyz2_; + std::vector ele_gen_MomLz2_; + std::vector ele_gen_MomLxy2_; + std::vector ele_gen_MomCTau02_; + + // final state hadrons + std::vector had_gen_PID_; + std::vector had_gen_Status_; + std::vector had_gen_Charge_; + std::vector had_gen_Pt_; + std::vector had_gen_Eta_; + std::vector had_gen_Phi_; + std::vector had_gen_Mass_; + + //quarks @ gen + std::vector quarks_gen_PID_; + std::vector quarks_gen_Status_; + std::vector quarks_gen_Charge_; + std::vector quarks_gen_Pt_; + std::vector quarks_gen_Eta_; + std::vector quarks_gen_Phi_; + std::vector quarks_gen_Mass_; + + //trigger infos + + bool passIsoMuTk18_ = 0; + bool passIsoMuTk20_ = 0; + bool passIsoMuTk22_ = 0; + bool passIsoMuTk24_ = 0; + bool passIsoMuTk27_ = 0; + bool passIsoMuTk17e_ = 0; + bool passIsoMuTk22e_ = 0; + + bool passIsoMu18_ = 0; + bool passIsoMu20_ = 0; + bool passIsoMu22_ = 0; + bool passIsoMu24_ = 0; + bool passIsoMu27_ = 0; + bool passIsoMu17e_ = 0; + bool passIsoMu22e_ = 0; + bool passTkMu17_ = 0; + bool passTkMu20_ = 0; + + bool passIsoMu24All_ = 0; + bool passIsoMu27All_ = 0; + + bool passDoubleMu17TrkIsoMu8_ = 0; + bool passDoubleMu17TrkIsoTkMu8_ = 0; + bool passDoubleTkMu17TrkIsoTkMu8_ = 0; + + bool passIsoEle27_ = 0; + bool passNonIsoEle115_ = 0; + bool passDoubleEle23andEle12DZ_ = 0; + bool passDoubleEle23andEle12_ = 0; + + bool passDoubleEle33TrkMW_ = 0; + bool passDoubleEle33MW_ = 0; + bool passDoubleEle33_ = 0; + + bool passDoubleMu33Ele33_ = 0; + + + //secondary verteces info due to mu + std::vector sv_mu_TrackSize_; + std::vector sv_mu_LXYSig_; + std::vector sv_mu_LXYZSig_; + std::vector sv_mu_LXY_; + std::vector sv_mu_LXYZ_; + std::vector sv_mu_mass_; + std::vector sv_mu_charge_; + std::vector sv_mu_eta_; + std::vector sv_mu_phi_; + std::vector sv_mu_pt_; + std::vector sv_mu_p_; + std::vector sv_mu_Beta_; + std::vector sv_mu_Gamma_; + std::vector sv_mu_CTau0_; + std::vector sv_mu_NDof_; + std::vector sv_mu_Chi2_; + std::vector sv_mu_Angle3D_; + std::vector sv_mu_Angle2D_; + + std::vector sv_mu_tracks_charge_; + std::vector sv_mu_tracks_eta_; + std::vector sv_mu_tracks_phi_; + std::vector sv_mu_tracks_pt_; + std::vector sv_mu_tracks_dxySig_; + std::vector sv_mu_tracks_dxy_; + std::vector sv_mu_tracks_dxyz_; + + std::vector sv_mu_tracks_Sumcharge_; + std::vector sv_mu_tracks_Sumpt_; + std::vector sv_mu_match_; + + //secondary verteces info due to ele + std::vector sv_ele_TrackSize_; + std::vector sv_ele_LXYSig_; + std::vector sv_ele_LXYZSig_; + std::vector sv_ele_LXY_; + std::vector sv_ele_LXYZ_; + std::vector sv_ele_mass_; + std::vector sv_ele_charge_; + std::vector sv_ele_eta_; + std::vector sv_ele_phi_; + std::vector sv_ele_pt_; + std::vector sv_ele_p_; + std::vector sv_ele_Beta_; + std::vector sv_ele_Gamma_; + std::vector sv_ele_CTau0_; + std::vector sv_ele_NDof_; + std::vector sv_ele_Chi2_; + std::vector sv_ele_Angle3D_; + std::vector sv_ele_Angle2D_; + + std::vector sv_ele_tracks_charge_; + std::vector sv_ele_tracks_eta_; + std::vector sv_ele_tracks_phi_; + std::vector sv_ele_tracks_pt_; + std::vector sv_ele_tracks_dxySig_; + std::vector sv_ele_tracks_dxy_; + std::vector sv_ele_tracks_dxyz_; + + std::vector sv_ele_tracks_Sumcharge_; + std::vector sv_ele_tracks_Sumpt_; + std::vector sv_ele_match_; + std::vector sv_ele_score_; + + //muon infos + std::vector mu_en_ ; + std::vector mu_pt_ ; + std::vector mu_eta_ ; + std::vector mu_phi_ ; + std::vector mu_et_ ; + std::vector mu_charge_ ; + std::vector mu_FirstGenMatch_ ; + std::vector mu_SecondGenMatch_ ; + std::vector mu_trackiso_ ; + std::vector mu_rhoIso_; + std::vector mu_pfSumChargedHadronPt_ ; + std::vector mu_pfSumNeutralHadronEt_ ; + std::vector mu_PFSumPhotonEt_ ; + std::vector mu_pfSumPUPt_ ; + std::vector mu_numberOfValidMuonHits_ ; + std::vector mu_emIso_ ; + std::vector mu_hadIso_ ; + std::vector mu_normalizedChi2_ ; + std::vector mu_numberOfMatchedStations_ ; + std::vector mu_numberOfValidPixelHits_ ; + std::vector mu_numberOftrackerLayersWithMeasurement_ ; + std::vector mu_numberOfpixelLayersWithMeasurement_ ; + std::vector mu_TrackQuality_ ; + std::vector mu_InnerTrackQuality_ ; + std::vector mu_pxTunePMuonBestTrack_ ; + std::vector mu_pyTunePMuonBestTrack_ ; + std::vector mu_pzTunePMuonBestTrack_ ; + std::vector mu_pTunePMuonBestTrack_ ; + std::vector mu_etaTunePMuonBestTrack_ ; + std::vector mu_LXYZ_ ; + std::vector mu_LXY_ ; + std::vector mu_ptTunePMuonBestTrack_ ; + std::vector mu_phiTunePMuonBestTrack_ ; + std::vector mu_thetaTunePMuonBestTrack_ ; + std::vector mu_chargeTunePMuonBestTrack_ ; + std::vector mu_dPToverPTTunePMuonBestTrack_ ; + std::vector mu_absdxyTunePMuonBestTrack_ ; + std::vector mu_absdxyErrorTunePMuonBestTrack_ ; + std::vector mu_absdxySigTunePMuonBestTrack_ ; + std::vector mu_absdzTunePMuonBestTrack_ ; + std::vector mu_absdzErrorTunePMuonBestTrack_ ; + std::vector mu_absdzSigTunePMuonBestTrack_ ; + std::vector mu_recoDeltaBeta_ ; + std::vector mu_recoiso_ ; + std::vector mu_isGlobalMuon_ ; + std::vector mu_isStandAloneMuon_ ; + std::vector mu_isPF_ ; + std::vector mu_isRPCMuon_ ; + std::vector mu_isTrackerMuon_ ; + std::vector mu_isGoodMuon_ ; + std::vector mu_isSoftMuon_ ; + std::vector mu_isLoose_ ; + std::vector mu_isTightMuon_ ; + std::vector mu_STAnHits_ ; + std::vector mu_STAnLost_ ; + std::vector mu_STAnStationsWithAnyHits_ ; + std::vector mu_STAnCscChambersWithAnyHits_ ; + std::vector mu_STAnDtChambersWithAnyHits_ ; + std::vector mu_STAnRpcChambersWithAnyHits_ ; + std::vector mu_STAinnermostStationWithAnyHits_ ; + std::vector mu_STAoutermostStationWithAnyHits_ ; + std::vector mu_STAnStationsWithValidHits_ ; + std::vector mu_STAnCscChambersWithValidHits_ ; + std::vector mu_STAnDtChambersWithValidHit_ ; + std::vector mu_STAnRpcChambersWithValidHits_ ; + std::vector mu_STAnValidMuonHits_ ; + std::vector mu_STAnValidCscHits_ ; + std::vector mu_STAnValidDtHits_ ; + std::vector mu_STAnValidRpcHits_ ; + std::vector mu_STAinnermostStationWithValidHits_ ; + std::vector mu_STAoutermostStationWithValidHits_ ; + std::vector mu_STATofDirection_ ; + std::vector mu_STATofNDof_ ; + std::vector mu_STATofTimeAtIpInOut_ ; + std::vector mu_STATofTimeAtIpInOutErr_ ; + std::vector mu_STATofTimeAtIpOutIn_ ; + std::vector mu_STATofTimeAtIpOutInErr_ ; + + //jet info + + std::vector jet_charge_ ; + std::vector jet_et_ ; + std::vector jet_pt_ ; + std::vector jet_eta_ ; + std::vector jet_phi_ ; + std::vector jet_theta_ ; + std::vector jet_en_ ; + std::vector jet_chargedEmEnergy_ ; + std::vector jet_neutralEmEnergyFraction_ ; + std::vector jet_chargedHadronEnergy_ ; + std::vector jet_neutralHadronEnergyFraction_ ; + std::vector jet_chargedMuEnergy_ ; + std::vector jet_chargedMuEnergyFraction_ ; + std::vector jet_chargedMultiplicity_ ; + std::vector jet_numberOfDaughters_ ; + std::vector jet_muonEnergy_ ; + std::vector jet_muonEnergyFraction_ ; + std::vector jet_muonMultiplicity_ ; + std::vector jet_neutralEmEnergy_ ; + std::vector jet_neutralHadronEnergy_ ; + std::vector jet_neutralHadronMultiplicity_ ; + std::vector jet_neutralMultiplicity_ ; + + //electron info + + std::vector ele_Et_; + std::vector ele_EtFromCaloEn_; + std::vector ele_pt_; + std::vector ele_etaSC_; + std::vector ele_phiSC_; + std::vector ele_phiWidth_; + std::vector ele_etaWidth_; + std::vector ele_energySC_; + std::vector ele_thetaSC_; + std::vector ele_preshowerEnergySC_; + std::vector ele_etaTrack_; + std::vector ele_phiTrack_; + std::vector ele_thetaTrack_; + std::vector ele_x_; + std::vector ele_y_; + std::vector ele_z_; + std::vector ele_e2x5Max_; + std::vector ele_e1x5_; + std::vector ele_e5x5_; + std::vector ele_e2x5MaxOver5x5_; + std::vector ele_e1x5Over5x5_; + std::vector ele_sigmaIetaIetaFull5x5_; + std::vector ele_e2x5MaxFull5x5_; + std::vector ele_e1x5Full5x5_; + std::vector ele_e5x5Full5x5_; + std::vector ele_e2x5MaxOver5x5Full5x5_; + std::vector ele_e1x5Over5x5Full5x5_; + std::vector ele_zTrackPositionAtVtx_; + std::vector ele_hadronicOverEm_; + std::vector ele_deltaEtaInSC_; + std::vector ele_deltaPhiInSC_; + std::vector ele_deltaEtaInSeedCluster_; + std::vector ele_deltaPhiInSeedCluster_; + std::vector ele_sigmaIetaIeta_; + std::vector ele_e2x5Right_; + std::vector ele_e2x5Left_; + std::vector ele_e2x5Top_; + std::vector ele_e2x5Bottom_; + std::vector ele_eMax_; + std::vector ele_eRight_; + std::vector ele_eLeft_; + std::vector ele_eTop_; + std::vector ele_eBottom_; + std::vector ele_e3x3_; + std::vector ele_frac51_; + std::vector ele_frac15_; + + std::vector ele_rawId_; + std::vector ele_ieta_; + std::vector ele_nbOfMissingHits_; + std::vector ele_charge_; + std::vector ele_isEcalDrivenSeed_; + std::vector ele_isPassConversionVeto_; + + std::vector ele_dxy_; + std::vector ele_dz_; + std::vector ele_rhoIso_; + std::vector ele_fbrem_; + std::vector ele_EoverP_; + std::vector ele_Xposition_; + std::vector ele_Yposition_; + std::vector ele_dr03TkSumPt_; + std::vector ele_hcalDepth1OverEcal_; + std::vector ele_hcalDepth2OverEcal_; + std::vector ele_dr03HcalDepth2TowerSumEt_; + std::vector ele_hcalDepth2TowerSumEtNoVeto_; + std::vector ele_hcalDepth1TowerSumEtNoVeto_; + std::vector ele_EcalPlusHcald1iso_; + std::vector ele_dr03EcalRecHitSumEt_; + std::vector ele_dr03HcalDepth1TowerSumEt_; + std::vector ele_dr03HcalDepth1TowerSumEtBc_; + std::vector ele_pfSumPhotonEt_; + std::vector ele_pfSumChargedHadronPt_; + std::vector ele_pfSumNeutralHadronEt_; + std::vector ele_pfSumPUPt_; + std::vector ele_pfDeltaBeta_; + std::vector ele_FirstGenMatch_; + std::vector ele_SecondGenMatch_; + + std::vector ele_Mva2016_; + std::vector ele_CutVeto_; + std::vector ele_CutLoose_; + std::vector ele_CutMedium_; + std::vector ele_CutTight_; + + /* + std::vector ele_Mva_; + std::vector ele_MvaFall17Iso_; + std::vector ele_MvaFall17NoIso_; + std::vector ele_CutBasedVeto_; + std::vector ele_CutBasedLoose_; + std::vector ele_CutBasedMedium_; + std::vector ele_CutBasedTight_; + */ + + //MET info + + float pfMet_et_ = -1000; + float pfMet_pt_ = -1000; + float pfMet_phi_ = -1000; + float pfMet_en_ = -1000; + float pfMet_px_ = -1000; + float pfMet_py_ = -1000; + float pfMet_pz_ = -1000; + float pfMet_sumEt_ = -1000; + float caloMet_pt_ = -1000; + float caloMet_phi_ = -1000; + + //bJet info + std::vector jet_btag_flavor_; + std::vector jet_btag_pfCSVv2IVF_discriminator_; + std::vector jet_btag_pt_; + std::vector jet_btag_eta_; + std::vector jet_btag_phi_; + +}; + + +#endif diff --git a/HeavyNeutralLeptonAnalysis/plugins/BuildFile.xml b/HeavyNeutralLeptonAnalysis/plugins/BuildFile.xml index 56e9546..293c673 100644 --- a/HeavyNeutralLeptonAnalysis/plugins/BuildFile.xml +++ b/HeavyNeutralLeptonAnalysis/plugins/BuildFile.xml @@ -1,7 +1,4 @@ - - - - + @@ -17,7 +14,8 @@ - + +# @@ -30,11 +28,22 @@ + + + + + + + + + + + diff --git a/HeavyNeutralLeptonAnalysis/plugins/HeavyNeutralLeptonAnalysis.cc b/HeavyNeutralLeptonAnalysis/plugins/HeavyNeutralLeptonAnalysis.cc index 2e5c51f..d1e7d31 100644 --- a/HeavyNeutralLeptonAnalysis/plugins/HeavyNeutralLeptonAnalysis.cc +++ b/HeavyNeutralLeptonAnalysis/plugins/HeavyNeutralLeptonAnalysis.cc @@ -19,6 +19,7 @@ // system include files #include +#include // user include files #include "FWCore/Framework/interface/Frameworkfwd.h" @@ -32,18 +33,42 @@ // for the plotting #include "HNL/HeavyNeutralLeptonAnalysis/interface/SmartSelectionMonitor.h" +// root includes +#include "TSystem.h" +#include "TFile.h" +#include "TTree.h" +#include "TCanvas.h" +#include "TH1F.h" +#include "TH2F.h" +#include "TProfile.h" +#include "TProfile2D.h" +#include "TEventList.h" +#include "TROOT.h" +#include "TNtuple.h" +#include "TLorentzVector.h" +#include + + //Load here all the dataformat that we will need +#include "DataFormats/FWLite/interface/Handle.h" +#include "DataFormats/FWLite/interface/Event.h" +#include "DataFormats/FWLite/interface/ChainEvent.h" + #include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h" #include "SimDataFormats/GeneratorProducts/interface/LHEEventProduct.h" #include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h" #include "DataFormats/VertexReco/interface/VertexFwd.h" +#include "DataFormats/MuonReco/interface/MuonTime.h" #include "DataFormats/PatCandidates/interface/Muon.h" #include "DataFormats/PatCandidates/interface/Electron.h" #include "DataFormats/PatCandidates/interface/Jet.h" #include "DataFormats/PatCandidates/interface/Photon.h" #include "DataFormats/PatCandidates/interface/MET.h" #include "DataFormats/PatCandidates/interface/Tau.h" +#include "FWCore/Common/interface/TriggerNames.h" +#include "DataFormats/Common/interface/TriggerResults.h" +#include "DataFormats/PatCandidates/interface/TriggerObjectStandAlone.h" #include "DataFormats/PatCandidates/interface/PackedTriggerPrescales.h" #include "DataFormats/PatCandidates/interface/GenericParticle.h" #include "DataFormats/Common/interface/MergeableCounter.h" @@ -61,27 +86,25 @@ #include "FWCore/Framework/interface/EventSetup.h" #include "FWCore/ParameterSet/interface/ParameterSet.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalClusterLazyTools.h" +#include "Geometry/CaloGeometry/interface/CaloGeometry.h" +#include "Geometry/Records/interface/CaloGeometryRecord.h" +#include "Geometry/Records/interface/CaloTopologyRecord.h" +#include "Geometry/CaloGeometry/interface/CaloSubdetectorGeometry.h" +#include "DataFormats/EcalDetId/interface/EcalSubdetector.h" +#include "DataFormats/EcalDetId/interface/EBDetId.h" +#include "RecoEcal/EgammaCoreTools/interface/EcalClusterTools.h" +#include "RecoEgamma/EgammaIsolationAlgos/interface/EgammaHcalIsolation.h" + #include "PhysicsTools/Utilities/interface/LumiReWeighting.h" #include "CommonTools/UtilAlgos/interface/TFileService.h" #include "EgammaAnalysis/ElectronTools/interface/ElectronEnergyCalibratorRun2.h" #include "EgammaAnalysis/ElectronTools/interface/PhotonEnergyCalibratorRun2.h" -// root includes +#include "HNL/DisplacedSVAssociator/interface/VertexAssociation.h" -#include "TSystem.h" -#include "TFile.h" -#include "TTree.h" -#include "TCanvas.h" -#include "TH1F.h" -#include "TH2F.h" -#include "TProfile.h" -#include "TProfile2D.h" -#include "TEventList.h" -#include "TROOT.h" -#include "TNtuple.h" -#include "TLorentzVector.h" -#include +#include "HNL/HeavyNeutralLeptonAnalysis/interface/BigNtuple.h" using namespace std; @@ -96,57 +119,102 @@ using namespace std; // This will improve performance in multithreaded jobs. class HeavyNeutralLeptonAnalysis : public edm::one::EDAnalyzer { - public: - explicit HeavyNeutralLeptonAnalysis(const edm::ParameterSet&); - ~HeavyNeutralLeptonAnalysis(); - - static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); - - - private: - virtual void beginJob() override; - virtual void analyze(const edm::Event&, const edm::EventSetup&) override; - virtual void initialize(const edm::Event&); - virtual void endJob() override; - - // ----------member data --------------------------- - bool isMC; - - edm::EDGetTokenT < reco::VertexCollection > vtxMiniAODToken_; - edm::EDGetTokenT < double > rhoToken_; - edm::EDGetTokenT < pat::MuonCollection > muonsMiniAODToken_; - edm::EDGetTokenT < pat::ElectronCollection > electronsMiniAODToken_; - edm::EDGetTokenT < EcalRecHitCollection > recHitEBToken_; - edm::EDGetTokenT < EcalRecHitCollection > recHitEEToken_; - edm::EDGetTokenT < pat::TauCollection > tausMiniAODToken_; - edm::EDGetTokenT < pat::PackedCandidateCollection > packedCandidateToken_; - edm::EDGetTokenT < pat::JetCollection > jetsMiniAODToken_; - edm::EDGetTokenT < pat::METCollection > pfMETAODToken_; - edm::EDGetTokenT < edm::TriggerResults > triggerResultsToken_; - edm::EDGetTokenT < edm::TriggerResults > metFilterResultsToken_; - edm::EDGetTokenT < reco::GenParticleCollection > genParticleToken_; - edm::EDGetTokenT < GenEventInfoProduct > genEventInfoToken_; - edm::EDGetTokenT < std::vector > PUInfoToken_; - edm::EDGetTokenT < LHEEventProduct > lheEventProductToken_; - - protected: - edm::Handle < reco::VertexCollection > vtxHandle; - edm::Handle < double > rhoHandle; - edm::Handle < pat::MuonCollection > muonsHandle; - edm::Handle < pat::ElectronCollection > electronsHandle; - edm::Handle < EcalRecHitCollection > recHitCollectionEBHandle; - edm::Handle < EcalRecHitCollection > recHitCollectionEEHandle; - edm::Handle < pat::TauCollection > tausHandle; - edm::Handle < pat::PackedCandidateCollection > pfCandidatesHandle; - edm::Handle < pat::JetCollection > jetsHandle; - edm::Handle < pat::METCollection > metsHandle; - edm::Handle < edm::TriggerResults > triggerResultsHandle; - edm::Handle < edm::TriggerResults > metFilterResultsHandle; - /* Only for MC */ - edm::Handle < reco::GenParticleCollection > genHandle; - edm::Handle < GenEventInfoProduct > genEventInfoHandle; - edm::Handle < std::vector > puInfoH; - edm::Handle < LHEEventProduct > lheEPHandle; +public: + explicit HeavyNeutralLeptonAnalysis(const edm::ParameterSet&); + ~HeavyNeutralLeptonAnalysis(); + + reco::VertexCollection getMatchedVertex_Muon(const pat::Muon & mu, const reco::VertexCollection& vertexCollection); + reco::VertexCollection PrimaryVertex( const reco::VertexCollection &vtx); + bool isAncestor(const reco::Candidate* ancestor, const reco::Candidate* particle); + double MatchGenMuon(const edm::Event&, reco::TrackRef BestTrack , int pdgId); + double MatchGenElectron(const edm::Event& iEvent, const pat::Electron& ele_, int pdgId); + double MatchGenVertex(const edm::Event& iEvent, reco::Vertex vertex, int pdgId); + static void fillDescriptions(edm::ConfigurationDescriptions& descriptions); + + +private: + virtual void beginJob() override; + virtual void analyze(const edm::Event&, const edm::EventSetup&) override; + virtual void initialize(const edm::Event&); + virtual void endJob() override; + + edm::Service fs; + TTree * tree_; + BigNtuple ntuple_; + + + //--------------Variables------------------------ + int debug; + float pvCompatibilityScore = .05; + + float npT = 0; + float npIT = 0; + + edm::LumiReWeighting Lumiweights_; + edm::LumiReWeighting LumiweightsUp_; + edm::LumiReWeighting LumiweightsDown_; + + bool isMC; + bool isMCSignal; + + SmartSelectionMonitor mon; + + //--------------template------------------------- + + // ----------member data --------------------------- + + edm::EDGetTokenT < reco::VertexCollection > vtxMiniAODToken_; + edm::EDGetTokenT < double > rhoToken_; + edm::EDGetTokenT < pat::MuonCollection > muonsMiniAODToken_; + edm::EDGetTokenT < pat::ElectronCollection > electronsMiniAODToken_; + edm::EDGetTokenT < EcalRecHitCollection > recHitEBToken_; + edm::EDGetTokenT < EcalRecHitCollection > recHitEEToken_; + edm::EDGetTokenT < pat::TauCollection > tausMiniAODToken_; + edm::EDGetTokenT < pat::PackedCandidateCollection > packedCandidateToken_; + edm::EDGetTokenT < pat::JetCollection > jetsMiniAODToken_; + edm::EDGetTokenT < pat::METCollection > pfMETAODToken_; + edm::EDGetTokenT < edm::TriggerResults > triggerResultsToken_; + edm::EDGetTokenT < edm::TriggerResults > metFilterResultsToken_; + edm::EDGetTokenT < reco::GenParticleCollection > genParticleToken_; + edm::EDGetTokenT < GenEventInfoProduct > genEventInfoToken_; + edm::EDGetTokenT < std::vector > PUInfoToken_; + edm::EDGetTokenT < LHEEventProduct > lheEventProductToken_; + edm::EDGetTokenT < reco::VertexCollection > inclusiveSecondaryVertices_; + + const std::vector bDiscriminators_; + + edm::EDGetTokenT > eleMvaToken_; + edm::EDGetTokenT > eleVetoToken_; + edm::EDGetTokenT > eleLooseToken_; + edm::EDGetTokenT > eleMediumToken_; + edm::EDGetTokenT > eleTightToken_; + +protected: + edm::Handle < reco::VertexCollection > vtxHandle; + edm::Handle < double > rhoHandle; + edm::Handle < pat::MuonCollection > muonsHandle; + edm::Handle > electronsHandle; + edm::Handle < EcalRecHitCollection > recHitCollectionEBHandle; + edm::Handle < EcalRecHitCollection > recHitCollectionEEHandle; + edm::Handle < pat::TauCollection > tausHandle; + edm::Handle < pat::PackedCandidateCollection > pfCandidatesHandle; + edm::Handle < pat::JetCollection > jetsHandle; + edm::Handle < pat::METCollection > metsHandle; + edm::Handle < edm::TriggerResults > triggerResultsHandle; + edm::Handle < edm::TriggerResults > metFilterResultsHandle; + edm::Handle < reco::VertexCollection > secondaryVertexHandle; + + edm::Handle > electronsMva; + edm::Handle > electronsVeto; + edm::Handle > electronsLoose; + edm::Handle > electronsMedium; + edm::Handle > electronsTight; + + /* Only for MC */ + edm::Handle < reco::GenParticleCollection > genHandle; + edm::Handle < GenEventInfoProduct > genEventInfoHandle; + edm::Handle < std::vector > puInfoH; + edm::Handle < LHEEventProduct > lheEPHandle; }; // @@ -161,7 +229,10 @@ class HeavyNeutralLeptonAnalysis : public edm::one::EDAnalyzer("debugLevel")), isMC(iConfig.getParameter("isMC")), + isMCSignal(iConfig.getParameter("isMCSignal")), vtxMiniAODToken_(mayConsume(iConfig.getParameter("vtxSrc"))), rhoToken_(mayConsume(iConfig.getParameter("rho"))), muonsMiniAODToken_(mayConsume(iConfig.getParameter("muonSrc"))), @@ -177,11 +248,35 @@ HeavyNeutralLeptonAnalysis::HeavyNeutralLeptonAnalysis(const edm::ParameterSet& genParticleToken_(mayConsume(iConfig.getParameter("genParticleSrc"))), genEventInfoToken_(mayConsume(iConfig.getParameter("genEventInfoProduct"))), PUInfoToken_(consumes>(iConfig.getParameter("PUInfo"))), - lheEventProductToken_(mayConsume(iConfig.getParameter("lheEventProducts"))) + lheEventProductToken_(mayConsume(iConfig.getParameter("lheEventProducts"))), + inclusiveSecondaryVertices_(consumes(iConfig.getParameter("SecondaryVertices"))), + bDiscriminators_(iConfig.getParameter >("bDiscriminators")), + eleMvaToken_(consumes>(iConfig.getParameter("electronsMva"))), + eleVetoToken_(consumes>(iConfig.getParameter("electronsVeto"))), + eleLooseToken_(consumes>(iConfig.getParameter("electronsLoose"))), + eleMediumToken_(consumes>(iConfig.getParameter("electronsMedium"))), + eleTightToken_(consumes>(iConfig.getParameter("electronsTight"))) { - //now do what ever initialization is needed - usesResource("TFileService"); + //now do what ever initialization is needed + usesResource("TFileService"); + + tree_ = fs->make("tree_", "tree"); + ntuple_.set_evtInfo(tree_); + if(isMC && isMCSignal){ + ntuple_.set_pv_genInfo(tree_); + ntuple_.set_sv_genInfo(tree_); + } + ntuple_.set_pileupInfo(tree_); + ntuple_.set_trigInfo(tree_); + ntuple_.set_pvInfo(tree_); + ntuple_.set_muInfo(tree_); + ntuple_.set_eleInfo(tree_); + ntuple_.set_eleIDInfo(tree_); + ntuple_.set_svInfo(tree_); + ntuple_.set_jetInfo(tree_); + ntuple_.set_metInfo(tree_); + ntuple_.set_bjetInfo(tree_); } @@ -205,6 +300,13 @@ void HeavyNeutralLeptonAnalysis::initialize(const edm::Event& iEvent){ iEvent.getByToken(pfMETAODToken_, metsHandle); iEvent.getByToken(triggerResultsToken_, triggerResultsHandle); iEvent.getByToken(metFilterResultsToken_, metFilterResultsHandle); + iEvent.getByToken(inclusiveSecondaryVertices_, secondaryVertexHandle); + + iEvent.getByToken(eleMvaToken_, electronsMva); + iEvent.getByToken(eleVetoToken_,electronsVeto); + iEvent.getByToken(eleLooseToken_,electronsLoose); + iEvent.getByToken(eleMediumToken_,electronsMedium); + iEvent.getByToken(eleTightToken_,electronsTight); if (isMC){ iEvent.getByToken(genParticleToken_, genHandle); @@ -218,74 +320,411 @@ void HeavyNeutralLeptonAnalysis::initialize(const edm::Event& iEvent){ // // member functions -// -// ------------ method called for each event ------------ -void -HeavyNeutralLeptonAnalysis::analyze(const edm::Event& iEvent, const edm::EventSetup& iSetup) +reco::VertexCollection HeavyNeutralLeptonAnalysis::getMatchedVertex_Muon(const pat::Muon & muon, const reco::VertexCollection& vertexCollection){ + reco::VertexCollection matchedVertices; + //cout<(muon.sourceCandidatePtr(0).get()); + //cout << "Packed: " << cand << endl; + if(!cand) { + cout << "THIS SHOULD NEVER HAPPEN! No packed candidated associated to muon?!" << endl; + } + //cout << "Packed: " << cand->hasTrackDetails() << " " << (cand->charge() != 0) << " " << (cand->numberOfHits() > 0) << endl; + if(!(cand->hasTrackDetails() && cand->charge() != 0 && cand->numberOfHits() > 0)) { + cout << "THIS SHOULD NEVER HAPPEN! Muon without track or matched to neutral?!" << endl; + } + //cout << cand->pseudoTrack().pt() << " " << cand->pseudoTrack().eta() << " " << cand->pseudoTrack().phi() << endl; + + for(reco::VertexCollection::const_iterator ss = vertexCollection.begin(); ss != vertexCollection.end(); ++ss) { + cout <<"new vertex"<tracks_begin(); tt != ss->tracks_end(); ++tt) { + //cout<<"Track " << (*tt)->pt() << " "<< (*tt)->eta()<< " " << (*tt)->phi() <pseudoTrack().pt() - tt->castTo()->pt()) / tt->castTo()->pt(); + cout << "match " << (cand->pseudoTrack().pt() == tt->castTo()->pt()) <<" dpt = " << dpt << endl; + if( (cand->pseudoTrack().pt() == tt->castTo()->pt()) || dpt < 0.001) { //Options here: innerTrack, globalTrack, muonBestTrack, outerTrack, pickyTrack, track + matchedVertices.push_back(*ss); + break; + } + } + } + return matchedVertices; +} +reco::VertexCollection HeavyNeutralLeptonAnalysis::PrimaryVertex( const reco::VertexCollection &vtx) { - using namespace edm; - - initialize(iEvent); + reco::VertexCollection allPVs; + + for(reco::VertexCollection::const_iterator PV = vtx.begin(); PV!=vtx.end();++PV) + { + if(!PV->isFake()) { + if(PV->ndof() > 4 && fabs(PV->position().z()) <= 24 && fabs(PV->position().rho()) <= 2 ) allPVs.push_back(*PV); + } + } + return allPVs; +} - reco::VertexCollection vtx; - if(vtxHandle.isValid()){ vtx = *vtxHandle;} +//======================================================================================// +bool HeavyNeutralLeptonAnalysis::isAncestor(const reco::Candidate* ancestor,const reco::Candidate* particle) +{ + if(ancestor == particle ) return true; + for(size_t i=0;i< particle->numberOfMothers();i++) + { + if(isAncestor(ancestor,particle->mother(i))) return true; + } + return false; +} +//===================================== Muon Gen Match ================================================// +double HeavyNeutralLeptonAnalysis::MatchGenMuon(const edm::Event& iEvent, reco::TrackRef BestTrack, int pdgId) { + double minDeltaR=999; + double recoGenDeltaR_ = 0.2; + int genIndex=-999; + for(size_t i = 0; isize(); i++){ + if(abs((*genHandle)[i].pdgId()) == 13 && abs((*genHandle)[i].mother()->pdgId()) == pdgId) { + const reco::Candidate * WBoson = (*genHandle)[i].mother(); + for(size_t j = 0; jsize(); j++){ + const reco::Candidate * part = (*genHandle)[j].mother(0) ; + if(part != nullptr && isAncestor( WBoson , part) ){ + if( (*genHandle)[j].status() == 1 && abs((*genHandle)[j].pdgId()) == 13 ){ + double dR=deltaR(BestTrack->eta(),BestTrack->phi(),(*genHandle)[j].eta(),(*genHandle)[j].phi()); + if (dRsize(); i++){ + if(abs((*genHandle)[i].pdgId()) == 11 && abs((*genHandle)[i].mother()->pdgId()) == pdgId) { + const reco::Candidate * WBoson = (*genHandle)[i].mother(); + for(size_t j = 0; jsize(); j++){ + const reco::Candidate * part = (*genHandle)[j].mother(0) ; + if(part != nullptr && isAncestor( WBoson , part) ){ + if( (*genHandle)[j].status() == 1 && abs((*genHandle)[j].pdgId()) == 11 ){ + double dR=deltaR(ele_.eta(),ele_.phi(),(*genHandle)[j].eta(),(*genHandle)[j].phi()); + if (dRsize(); i++){ + if(abs((*genHandle)[i].pdgId()) == pdgId && (*genHandle)[i].mother()->pdgId() == 9900012) { + const reco::Candidate * HNL = (*genHandle)[i].mother(); + for(size_t j = 0; jsize(); j++){ + const reco::Candidate * part = (*genHandle)[j].mother(0) ; + if(part != nullptr && isAncestor( HNL , part) ){ + if( (*genHandle)[j].status() == 1){ + float vx = vertex.x(), vy = vertex.y(), vz = vertex.z(); + float gx = (*genHandle)[j].vx(), gy = (*genHandle)[j].vy(), gz = (*genHandle)[j].vz(); + //float metric1 = std::sqrt(((gx-vx)*(gx-vx))/(gx*gx) + ((gy-vy)*(gy-vy))/(gy*gy) + ((gz-vz)*(gz-vz))/(gz*gz)); + float metric = std::sqrt(((gx-vx)*(gx-vx)) + ((gy-vy)*(gy-vy)) + ((gz-vz)*(gz-vz))); + if (metricsize(); i++){ + if(abs((*genHandle)[i].pdgId()) == 9900012 && abs((*genHandle)[i].mother()->pdgId()) == 24) { + const reco::Candidate * WBoson = (*genHandle)[i].mother(); + for(size_t j = 0; jsize(); j++){ + const reco::Candidate * part = (*genHandle)[j].mother(0) ; + if(part != nullptr && isAncestor( WBoson , part) ){ + ntuple_.fill_pv_genInfo( (*genHandle)[j] , WBoson); + } + } + } + } + //mu && ele @ sv + for(size_t i = 0; isize(); i++){ + if(abs((*genHandle)[i].pdgId()) < 5 && abs((*genHandle)[i].mother()->pdgId()) == 9900012) { + const reco::Candidate * HNL = (*genHandle)[i].mother(); + for(size_t j = 0; jsize(); j++){ + const reco::Candidate * part = (*genHandle)[j].mother(0) ; + if(part != nullptr && isAncestor( HNL , part) ){ + ntuple_.fill_sv_genInfo( (*genHandle)[j] , HNL); + } + } + } + } + } + //============================================================= + // + // Primary Vertex + // + //============================================================= + if(!vtxHandle.isValid()) return; + reco::VertexCollection pvs = PrimaryVertex(*vtxHandle); + if(!pvs.size()) return; + ntuple_.fill_pvInfo(pvs); + //============================================================= + // + // Trigger Info + // + //============================================================= + const edm::TriggerResults triggerResults = *triggerResultsHandle.product(); + const edm::TriggerNames& trigNames = iEvent.triggerNames(triggerResults); + ntuple_.fill_trigInfo(triggerResults, trigNames); + //============================================================= + // + // Pile Up Info + // + //============================================================= + if(isMC){ + + std::vector::const_iterator PVI; + for(PVI = puInfoH->begin(); PVI != puInfoH->end(); ++PVI) { + int BX = PVI->getBunchCrossing(); + if(BX == 0) { + npT += PVI->getTrueNumInteractions(); + npIT = PVI->getPU_NumInteractions(); + } + } + // calculate weight using above code + float pu_weight = Lumiweights_.weight(npT); + float pu_weightUp = LumiweightsUp_.weight(npT); + float pu_weightDown = LumiweightsDown_.weight(npT); + + ntuple_.fill_pileupInfo(npT, npIT, pu_weight, pu_weightUp, pu_weightDown); + + } + //============================================================= + // + // LHE Info + // + //============================================================= /* - double rho = 0; - if(rhoHandle.isValid()){ rho = *rhoHandle;} - */ - + if (isMC) + { + if(lheEPHandle.isValid()){ + mon.fillHisto ("nup", "", lheEPHandle->hepeup ().NUP, 1); + //if (lheEPHandle->hepeup().NUP > 5) continue; to be check + mon.fillHisto ("nupfilt", "", lheEPHandle->hepeup ().NUP, 1); + } + else{ + printf("Handle to externalLHEProducer is invalid"); + } + }*/ + //============================================================= + // + // Method for Muon Tree + // + //============================================================= + std::string muon; pat::MuonCollection muons; - if(muonsHandle.isValid()){ muons = *muonsHandle;} - - pat::ElectronCollection electrons; - if(electronsHandle.isValid()){ electrons = *electronsHandle;} + vector goodMuons; + vector looseMuons; + + if(muonsHandle.isValid()){ + muons = *muonsHandle; + for (const pat::Muon mu : muons) { + if( mu.pt() < 0.0 ) continue; + if (!( fabs(mu.eta()) < 2.4 && mu.pt() > 5. )) continue; + goodMuons.push_back(mu); + if (!mu.isLooseMuon()) continue; + looseMuons.push_back(mu); + } + } + + for (const pat::Muon mu : goodMuons){ + double rho = *(rhoHandle.product()); + reco::TrackRef bestTrack = mu.muonBestTrack(); + double matching_1stmu = (isMC && isMCSignal) ? MatchGenMuon(iEvent, bestTrack, 24) : -999; + double matching_2ndmu = (isMC && isMCSignal) ? MatchGenMuon(iEvent, bestTrack, 9900012) : -999; + //added the matching + ntuple_.fill_muInfo(mu, pvs.at(0) , rho ,matching_1stmu , matching_2ndmu); + } + + // lambda function to sort this muons + std::sort(looseMuons.begin(), looseMuons.end(), [](pat::Muon a, pat::Muon b) {return a.pt() < b.pt(); }); + ////////////////////////////////////////////// EcalRecHitCollection recHitCollectionEB; if(recHitCollectionEBHandle.isValid()){ recHitCollectionEB = *recHitCollectionEBHandle;} - + EcalRecHitCollection recHitCollectionEE; if(recHitCollectionEEHandle.isValid()){ recHitCollectionEE = *recHitCollectionEEHandle;} - + //============================================================= + // + // Method for electrons + // + //============================================================= + vector looseElectrons; + // using iterator to get ref for ele + for(auto ele = electronsHandle->begin(); ele != electronsHandle->end(); ++ele){ + if(ele->gsfTrack().isNull() || ele->pt() < 5 || fabs(ele->eta()) > 2.5 ) continue; + + double rho = *(rhoHandle.product()); + + auto eleRef = edm::Ref>(electronsHandle, (ele - electronsHandle->begin())); + + std::auto_ptr recHitEcal; + recHitEcal.reset(new EcalClusterLazyTools( iEvent, iSetup, recHitEBToken_, recHitEEToken_ )); + + double matching_1stele = (isMC && isMCSignal) ? MatchGenElectron(iEvent, *ele , 24 ) : -999; + double matching_2ndele = (isMC && isMCSignal) ? MatchGenElectron(iEvent, *ele , 9900012) : -999; + + ntuple_.fill_eleInfo(*ele, pvs.at(0), rho , matching_1stele, matching_2ndele, recHitEcal); + if(ele->full5x5_sigmaIetaIeta() < 0.036 && ele->passConversionVeto() == 1) looseElectrons.push_back(*ele); + + float ele_Mva_ = ((*electronsMva)[eleRef]); + bool ele_Veto_ = ((*electronsVeto)[eleRef]); + bool ele_Loose_ = ((*electronsLoose)[eleRef]); + bool ele_Medium_ = ((*electronsMedium)[eleRef]); + bool ele_Tight_ = ((*electronsTight)[eleRef]); + + ntuple_.fill_eleIDInfo(ele_Mva_, ele_Veto_, ele_Loose_ , ele_Medium_, ele_Tight_); + + } + // lambda function to sort this electrons + std::sort(looseElectrons.begin(), looseElectrons.end(), [](pat::Electron a, pat::Electron b) {return a.pt() < b.pt(); }); + //============================================================= + // + // Secondary Vertex + // + //============================================================= + if(secondaryVertexHandle.isValid()){ + //sv due to muon + if(looseMuons.size()){ + pat::Muon muonHNL = looseMuons[0]; + reco::VertexCollection bestVertices_mu = getMatchedVertex_Muon(muonHNL, *secondaryVertexHandle); + // check if SV doesn't match with the PV + for (const reco::Vertex& vtx_mu : bestVertices_mu){ + float x = vtx_mu.x(), y = vtx_mu.y(), z = vtx_mu.z(); + float dx = x - pvs.at(0).x() , dy = y - pvs.at(0).y(), dz = z - pvs.at(0).z(); + float selIVFIsPVScore = std::sqrt((dx/x)*(dx/x) + (dy/y)*(dy/y) + (dz/z)*(dz/z)); + if (selIVFIsPVScore < pvCompatibilityScore) continue; + double matching_vtx = (isMC && isMCSignal) ? MatchGenVertex(iEvent, vtx_mu , 13) : -999; + ntuple_.fill_sv_mu_Info(vtx_mu, pvs.at(0), matching_vtx); + } + } + //sv due to electron + if(looseElectrons.size()){ + pat::Electron electronHNL = looseElectrons[0]; + VertexAssociation JVAIVF("IVF", pvs.at(0), debug); + reco::VertexCollection::const_iterator vtxIter = secondaryVertexHandle->begin(); + for(; vtxIter != secondaryVertexHandle->end(); ++vtxIter ) { + JVAIVF.addVertex(*vtxIter); + } + const std::pair bestVertexPair = JVAIVF.getBestVertex(electronHNL, "oneOverR"); + const reco::Vertex vtx_ele = bestVertexPair.first; + const float bestVertexScore_ele = bestVertexPair.second; + float x = vtx_ele.x(), y = vtx_ele.y(), z = vtx_ele.z(); + float dx = x - pvs.at(0).x() , dy = y - pvs.at(0).y(), dz = z - pvs.at(0).z(); + float selIVFIsPVScore = std::sqrt((dx/x)*(dx/x) + (dy/y)*(dy/y) + (dz/z)*(dz/z)); + if (selIVFIsPVScore > pvCompatibilityScore) { + double matching_vtx = (isMC && isMCSignal) ? MatchGenVertex(iEvent, vtx_ele , 11) : -999; + ntuple_.fill_sv_ele_Info(vtx_ele, pvs.at(0), matching_vtx , bestVertexScore_ele); + } + } + } + //============================================================= + // + // Jets + // + //============================================================= + pat::JetCollection jets; + if(jetsHandle.isValid()){ + jets = *jetsHandle; + for (const pat::Jet jet : jets) { + if( jet.pt() < 0.0 ) continue; + if (!( fabs(jet.eta()) < 3 && jet.pt() > 5. )) continue; + else ntuple_.fill_jetInfo(jet); + int flavor = std::abs( jet.partonFlavour() ); + for( const std::string &bDiscr : bDiscriminators_ ) + { + ntuple_.fill_bjetInfo(jet, bDiscr, flavor); + } + } + } + //============================================================= + // + // Missing Energy + // + //============================================================= + pat::METCollection mets; + if(metsHandle.isValid()){ + mets = *metsHandle; + const pat::MET met = mets.front(); + ntuple_.fill_metInfo(met); + } pat::TauCollection taus; if(tausHandle.isValid()){ taus = *tausHandle;} pat::PackedCandidateCollection pfCandidates; if (pfCandidatesHandle.isValid()) { pfCandidates = *pfCandidatesHandle; } - pat::JetCollection jets; - if(jetsHandle.isValid()){ jets = *jetsHandle;} - pat::METCollection mets; - if(metsHandle.isValid()){ mets = *metsHandle;} - pat::MET met = mets[0]; - - + tree_->Fill(); } - - + // ------------ method called once each job just before starting event loop ------------ void HeavyNeutralLeptonAnalysis::beginJob() { + if(isMC){ + Lumiweights_ = edm::LumiReWeighting("puData_2016GH_central.root","MCpileUp_25ns_Recent2016.root", "pileup", "pileup"); + LumiweightsUp_ = edm::LumiReWeighting("puData_2016GH_up.root", "MCpileUp_25ns_Recent2016.root", "pileup", "pileup"); + LumiweightsDown_ = edm::LumiReWeighting("puData_2016GH_down.root", "MCpileUp_25ns_Recent2016.root", "pileup", "pileup"); + } } // ------------ method called once each job just after ending the event loop ------------ void -HeavyNeutralLeptonAnalysis::endJob() -{ +HeavyNeutralLeptonAnalysis::endJob() { } // ------------ method fills 'descriptions' with the allowed parameters for the module ------------ -void -HeavyNeutralLeptonAnalysis::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { +void HeavyNeutralLeptonAnalysis::fillDescriptions(edm::ConfigurationDescriptions& descriptions) { //The following says we do not know what parameters are allowed so do no validation // Please change this to state exactly what you do use, even if it is no parameters edm::ParameterSetDescription desc; desc.setUnknown(); descriptions.addDefault(desc); } - //define this as a plug-in DEFINE_FWK_MODULE(HeavyNeutralLeptonAnalysis); diff --git a/HeavyNeutralLeptonAnalysis/python/ele_Sequence_cff.py b/HeavyNeutralLeptonAnalysis/python/ele_Sequence_cff.py new file mode 100644 index 0000000..ef7ca1e --- /dev/null +++ b/HeavyNeutralLeptonAnalysis/python/ele_Sequence_cff.py @@ -0,0 +1,11 @@ +import FWCore.ParameterSet.Config as cms +from PhysicsTools.SelectorUtils.tools.vid_id_tools import * +def addElectronSequence(process): + + process.load("Configuration.StandardSequences.GeometryRecoDB_cff") + switchOnVIDElectronIdProducer(process, DataFormat.MiniAOD) + electronModules = ['RecoEgamma.ElectronIdentification.Identification.mvaElectronID_Spring16_GeneralPurpose_V1_cff', + 'RecoEgamma.ElectronIdentification.Identification.cutBasedElectronID_Summer16_80X_V1_cff'] + for idmod in electronModules: setupAllVIDIdsInModule(process,idmod,setupVIDElectronSelection) + process.ele_Sequence = cms.Sequence(process.egmGsfElectronIDSequence) + diff --git a/HeavyNeutralLeptonAnalysis/src/.BigNtuple.cc.swp b/HeavyNeutralLeptonAnalysis/src/.BigNtuple.cc.swp new file mode 100644 index 0000000..e69de29 diff --git a/HeavyNeutralLeptonAnalysis/src/BigNtuple.cc b/HeavyNeutralLeptonAnalysis/src/BigNtuple.cc new file mode 100644 index 0000000..48fb897 --- /dev/null +++ b/HeavyNeutralLeptonAnalysis/src/BigNtuple.cc @@ -0,0 +1,1225 @@ +#include "HNL/HeavyNeutralLeptonAnalysis/interface/BigNtuple.h" +#include "TVector3.h" + +void BigNtuple::set_evtInfo(TTree* tree) { + tree->Branch("run" , &run_, "run/i"); + tree->Branch("lumi", &lumi_, "lumi/i"); + tree->Branch("evt" , &evt_, "evt/i"); +} + +void BigNtuple::fill_evtInfo(const edm::EventID& id) { + lumi_ = id.run(); + run_ = id.luminosityBlock(); + evt_ = id.event(); +} + +void BigNtuple::set_pv_genInfo(TTree* tree) { + + tree->Branch("mu_gen_PID1" , &mu_gen_PID1_); + tree->Branch("mu_gen_Status1",&mu_gen_Status1_); + tree->Branch("mu_gen_Charge1",&mu_gen_Charge1_); + tree->Branch("mu_gen_Pt1",&mu_gen_Pt1_); + tree->Branch("mu_gen_Eta1",&mu_gen_Eta1_); + tree->Branch("mu_gen_Phi1",&mu_gen_Phi1_); + tree->Branch("mu_gen_VX1",&mu_gen_VX1_); + tree->Branch("mu_gen_VY1",&mu_gen_VY1_); + tree->Branch("mu_gen_VZ1",&mu_gen_VZ1_); + tree->Branch("mu_gen_PartVLxy1",&mu_gen_Lxy1_); + tree->Branch("mu_gen_PartVLxyz1",&mu_gen_Lxyz1_); + tree->Branch("mu_gen_MomPID1",&mu_gen_MomPID1_); + tree->Branch("mu_gen_MomStatus1",&mu_gen_MomStatus1_); + tree->Branch("mu_gen_MomMass1",&mu_gen_MomMass1_); + tree->Branch("mu_gen_MomCharge1",&mu_gen_MomCharge1_); + tree->Branch("mu_gen_MomPt1",&mu_gen_MomPt1_); + tree->Branch("mu_gen_MomEta1",&mu_gen_MomEta1_); + tree->Branch("mu_gen_MomPhi1",&mu_gen_MomPhi1_); + tree->Branch("mu_gen_MomBeta1",&mu_gen_MomBeta1_); + tree->Branch("mu_gen_MomGamma1",&mu_gen_MomGamma1_); + tree->Branch("mu_gen_MomLxyz1",&mu_gen_MomLxyz1_); + tree->Branch("mu_gen_MomLz1",&mu_gen_MomLz1_); + tree->Branch("mu_gen_MomLxy1",&mu_gen_MomLxy1_); + + tree->Branch("ele_gen_PID1" , &ele_gen_PID1_); + tree->Branch("ele_gen_Status1",&ele_gen_Status1_); + tree->Branch("ele_gen_Charge1",&ele_gen_Charge1_); + tree->Branch("ele_gen_Pt1",&ele_gen_Pt1_); + tree->Branch("ele_gen_Eta1",&ele_gen_Eta1_); + tree->Branch("ele_gen_Phi1",&ele_gen_Phi1_); + tree->Branch("ele_gen_VX1",&ele_gen_VX1_); + tree->Branch("ele_gen_VY1",&ele_gen_VY1_); + tree->Branch("ele_gen_VZ1",&ele_gen_VZ1_); + tree->Branch("ele_gen_PartVLxy1",&ele_gen_Lxy1_); + tree->Branch("ele_gen_PartVLxyz1",&ele_gen_Lxyz1_); + tree->Branch("ele_gen_MomPID1",&ele_gen_MomPID1_); + tree->Branch("ele_gen_MomStatus1",&ele_gen_MomStatus1_); + tree->Branch("ele_gen_MomMass1",&ele_gen_MomMass1_); + tree->Branch("ele_gen_MomCharge1",&ele_gen_MomCharge1_); + tree->Branch("ele_gen_MomPt1",&ele_gen_MomPt1_); + tree->Branch("ele_gen_MomEta1",&ele_gen_MomEta1_); + tree->Branch("ele_gen_MomPhi1",&ele_gen_MomPhi1_); + tree->Branch("ele_gen_MomBeta1",&ele_gen_MomBeta1_); + tree->Branch("ele_gen_MomGamma1",&ele_gen_MomGamma1_); + tree->Branch("ele_gen_MomLxyz1",&ele_gen_MomLxyz1_); + tree->Branch("ele_gen_MomLz1",&ele_gen_MomLz1_); + tree->Branch("ele_gen_MomLxy1",&ele_gen_MomLxy1_); + +} +void BigNtuple::fill_pv_genInfo(const reco::GenParticle prt , const reco::Candidate* mom){ + + float vx = prt.vx(), vy = prt.vy(), vz = prt.vz(); + float mx = mom->vx(), my = mom->vy(), mz = mom->vz(); + float dx = vx - mx, dy = vy - my, dz = vz - mz; + + if(abs(prt.pdgId()) == 13 ){ + mu_gen_PID1_.push_back(prt.pdgId() ); + mu_gen_Status1_.push_back(prt.status()); + // genkinematics + mu_gen_Pt1_.push_back(prt.pt()); + mu_gen_Eta1_.push_back(prt.eta()); + mu_gen_Phi1_.push_back(prt.phi()); + mu_gen_Charge1_.push_back(prt.charge()); + // vertexposition + mu_gen_VX1_.push_back(vx); + mu_gen_VY1_.push_back(vy); + mu_gen_VZ1_.push_back(vz); + // vertexflight + mu_gen_Lxy1_.push_back(std::sqrt( vx * vx + vy * vy )); + mu_gen_Lxyz1_.push_back(std::sqrt( vx * vx + vy * vy + vz * vz)); + // mother beta, gamma,ctau + float beta_mom = mom->p() / mom->energy(); + float gamma_mom = mom->energy() / mom->mass(); + // mother quantities related to thedecay + // gen id andstatus + mu_gen_MomPID1_.push_back(mom->pdgId()); + mu_gen_MomStatus1_.push_back(mom->status()); + mu_gen_MomCharge1_.push_back(mom->charge()); + mu_gen_MomMass1_.push_back(mom->mass()); + mu_gen_MomPt1_.push_back(mom->pt()); + mu_gen_MomEta1_.push_back(mom->eta()); + mu_gen_MomPhi1_.push_back(mom->phi()); + mu_gen_MomBeta1_.push_back(beta_mom); + mu_gen_MomGamma1_.push_back(gamma_mom); + mu_gen_MomLxyz1_.push_back(std::sqrt(dx*dx + dy*dy + dz*dz)); + mu_gen_MomLz1_.push_back(dz); + mu_gen_MomLxy1_.push_back(std::sqrt(dx*dx + dy*dy)); + } + + if(abs(prt.pdgId()) == 11 ){ + ele_gen_PID1_.push_back(prt.pdgId() ); + ele_gen_Status1_.push_back(prt.status()); + // genkinematics + ele_gen_Pt1_.push_back(prt.pt()); + ele_gen_Eta1_.push_back(prt.eta()); + ele_gen_Phi1_.push_back(prt.phi()); + ele_gen_Charge1_.push_back(prt.charge()); + // vertexposition + ele_gen_VX1_.push_back(vx); + ele_gen_VY1_.push_back(vy); + ele_gen_VZ1_.push_back(vz); + // vertexflight + ele_gen_Lxy1_.push_back(std::sqrt( vx * vx + vy * vy )); + ele_gen_Lxyz1_.push_back(std::sqrt( vx * vx + vy * vy + vz * vz)); + // mother beta, gamma,ctau + float beta_mom = mom->p() / mom->energy(); + float gamma_mom = mom->energy() / mom->mass(); + // mother quantities related to thedecay + // gen id andstatus + ele_gen_MomPID1_.push_back(mom->pdgId()); + ele_gen_MomStatus1_.push_back(mom->status()); + ele_gen_MomCharge1_.push_back(mom->charge()); + ele_gen_MomMass1_.push_back(mom->mass()); + ele_gen_MomPt1_.push_back(mom->pt()); + ele_gen_MomEta1_.push_back(mom->eta()); + ele_gen_MomPhi1_.push_back(mom->phi()); + ele_gen_MomBeta1_.push_back(beta_mom); + ele_gen_MomGamma1_.push_back(gamma_mom); + ele_gen_MomLxyz1_.push_back(std::sqrt(dx*dx + dy*dy + dz*dz)); + ele_gen_MomLz1_.push_back(dz); + ele_gen_MomLxy1_.push_back(std::sqrt(dx*dx + dy*dy)); + } +} + +void BigNtuple::set_sv_genInfo(TTree* tree) { + + tree->Branch("mu_gen_PID2" , &mu_gen_PID2_); + tree->Branch("mu_gen_Status2",&mu_gen_Status2_); + tree->Branch("mu_gen_Charge2",&mu_gen_Charge2_); + tree->Branch("mu_gen_Pt2",&mu_gen_Pt2_); + tree->Branch("mu_gen_Eta2",&mu_gen_Eta2_); + tree->Branch("mu_gen_Phi2",&mu_gen_Phi2_); + tree->Branch("mu_gen_VX2",&mu_gen_VX2_); + tree->Branch("mu_gen_VY2",&mu_gen_VY2_); + tree->Branch("mu_gen_VZ2",&mu_gen_VZ2_); + tree->Branch("mu_gen_PartVLxy2",&mu_gen_Lxy2_); + tree->Branch("mu_gen_PartVLxyz2",&mu_gen_Lxyz2_); + tree->Branch("mu_gen_MomPID2",&mu_gen_MomPID2_); + tree->Branch("mu_gen_MomStatus2",&mu_gen_MomStatus2_); + tree->Branch("mu_gen_MomMass2",&mu_gen_MomMass2_); + tree->Branch("mu_gen_MomCharge2",&mu_gen_MomCharge2_); + tree->Branch("mu_gen_MomPt2",&mu_gen_MomPt2_); + tree->Branch("mu_gen_MomEta2",&mu_gen_MomEta2_); + tree->Branch("mu_gen_MomPhi2",&mu_gen_MomPhi2_); + tree->Branch("mu_gen_MomBeta2",&mu_gen_MomBeta2_); + tree->Branch("mu_gen_MomGamma2",&mu_gen_MomGamma2_); + tree->Branch("mu_gen_MomLxyz2",&mu_gen_MomLxyz2_); + tree->Branch("mu_gen_MomLz2",&mu_gen_MomLz2_); + tree->Branch("mu_gen_MomLxy2",&mu_gen_MomLxy2_); + tree->Branch("mu_gen_MomCTau02" ,&mu_gen_MomCTau02_); + + tree->Branch("ele_gen_PID2" , &ele_gen_PID2_); + tree->Branch("ele_gen_Status2",&ele_gen_Status2_); + tree->Branch("ele_gen_Charge2",&ele_gen_Charge2_); + tree->Branch("ele_gen_Pt2",&ele_gen_Pt2_); + tree->Branch("ele_gen_Eta2",&ele_gen_Eta2_); + tree->Branch("ele_gen_Phi2",&ele_gen_Phi2_); + tree->Branch("ele_gen_VX2",&ele_gen_VX2_); + tree->Branch("ele_gen_VY2",&ele_gen_VY2_); + tree->Branch("ele_gen_VZ2",&ele_gen_VZ2_); + tree->Branch("ele_gen_PartVLxy2",&ele_gen_Lxy2_); + tree->Branch("ele_gen_PartVLxyz2",&ele_gen_Lxyz2_); + tree->Branch("ele_gen_MomPID2",&ele_gen_MomPID2_); + tree->Branch("ele_gen_MomStatus2",&ele_gen_MomStatus2_); + tree->Branch("ele_gen_MomMass2",&ele_gen_MomMass2_); + tree->Branch("ele_gen_MomCharge2",&ele_gen_MomCharge2_); + tree->Branch("ele_gen_MomPt2",&ele_gen_MomPt2_); + tree->Branch("ele_gen_MomEta2",&ele_gen_MomEta2_); + tree->Branch("ele_gen_MomPhi2",&ele_gen_MomPhi2_); + tree->Branch("ele_gen_MomBeta2",&ele_gen_MomBeta2_); + tree->Branch("ele_gen_MomGamma2",&ele_gen_MomGamma2_); + tree->Branch("ele_gen_MomLxyz2",&ele_gen_MomLxyz2_); + tree->Branch("ele_gen_MomLz2",&ele_gen_MomLz2_); + tree->Branch("ele_gen_MomLxy2",&ele_gen_MomLxy2_); + tree->Branch("ele_gen_MomCTau02" ,&ele_gen_MomCTau02_); + + tree->Branch("had_gen_PID" , &had_gen_PID_); + tree->Branch("had_gen_Status",&had_gen_Status_); + tree->Branch("had_gen_Charge",&had_gen_Charge_); + tree->Branch("had_gen_Pt",&had_gen_Pt_); + tree->Branch("had_gen_Eta",&had_gen_Eta_); + tree->Branch("had_gen_Phi",&had_gen_Phi_); + tree->Branch("had_gen_Mass",&had_gen_Mass_); + + tree->Branch("quarks_gen_PID" , &quarks_gen_PID_); + tree->Branch("quarks_gen_Status",&quarks_gen_Status_); + tree->Branch("quarks_gen_Charge",&quarks_gen_Charge_); + tree->Branch("quarks_gen_Pt",&quarks_gen_Pt_); + tree->Branch("quarks_gen_Eta",&quarks_gen_Eta_); + tree->Branch("quarks_gen_Mass",&quarks_gen_Mass_); + +} + +void BigNtuple::fill_sv_genInfo(const reco::GenParticle prt , const reco::Candidate* mom){ + + float vx = prt.vx(), vy = prt.vy(), vz = prt.vz(); + float mx = mom->vx(), my = mom->vy(), mz = mom->vz(); + float dx = vx - mx, dy = vy - my, dz = vz - mz; + + if(abs(prt.pdgId()) == 13 ){ + mu_gen_PID2_.push_back(prt.pdgId() ); + mu_gen_Status2_.push_back(prt.status()); + // genkinematics + mu_gen_Pt2_.push_back(prt.pt()); + mu_gen_Eta2_.push_back(prt.eta()); + mu_gen_Phi2_.push_back(prt.phi()); + mu_gen_Charge2_.push_back(prt.charge()); + // vertexposition + mu_gen_VX2_.push_back(vx); + mu_gen_VY2_.push_back(vy); + mu_gen_VZ2_.push_back(vz); + // vertexflight + mu_gen_Lxy2_.push_back(std::sqrt( vx * vx + vy * vy )); + mu_gen_Lxyz2_.push_back(std::sqrt( vx * vx + vy * vy + vz * vz)); + // mother beta, gamma,ctau + float beta_mom = mom->p() / mom->energy(); + float gamma_mom = mom->energy() / mom->mass(); + // mother quantities related to thedecay + // gen id andstatus + mu_gen_MomPID2_.push_back(mom->pdgId()); + mu_gen_MomStatus2_.push_back(mom->status()); + mu_gen_MomCharge2_.push_back(mom->charge()); + mu_gen_MomMass2_.push_back(mom->mass()); + mu_gen_MomPt2_.push_back(mom->pt()); + mu_gen_MomEta2_.push_back(mom->eta()); + mu_gen_MomPhi2_.push_back(mom->phi()); + mu_gen_MomBeta2_.push_back(beta_mom); + mu_gen_MomGamma2_.push_back(gamma_mom); + mu_gen_MomLxyz2_.push_back(std::sqrt(dx*dx + dy*dy + dz*dz)); + mu_gen_MomLz2_.push_back(dz); + mu_gen_MomLxy2_.push_back(std::sqrt(dx*dx + dy*dy)); + mu_gen_MomCTau02_.push_back(std::sqrt(dx*dx + dy*dy + dz*dz) / (beta_mom * gamma_mom)); + } + + if(abs(prt.pdgId()) == 11 ){ + ele_gen_PID2_.push_back(prt.pdgId() ); + ele_gen_Status2_.push_back(prt.status()); + // genkinematics + ele_gen_Pt2_.push_back(prt.pt()); + ele_gen_Eta2_.push_back(prt.eta()); + ele_gen_Phi2_.push_back(prt.phi()); + ele_gen_Charge2_.push_back(prt.charge()); + // vertexposition + ele_gen_VX2_.push_back(vx); + ele_gen_VY2_.push_back(vy); + ele_gen_VZ2_.push_back(vz); + // vertexflight + ele_gen_Lxy2_.push_back(std::sqrt( vx * vx + vy * vy )); + ele_gen_Lxyz2_.push_back(std::sqrt( vx * vx + vy * vy + vz * vz)); + // mother beta, gamma,ctau + float beta_mom = mom->p() / mom->energy(); + float gamma_mom = mom->energy() / mom->mass(); + // mother quantities related to thedecay + // gen id andstatus + ele_gen_MomPID2_.push_back(mom->pdgId()); + ele_gen_MomStatus2_.push_back(mom->status()); + ele_gen_MomCharge2_.push_back(mom->charge()); + ele_gen_MomMass2_.push_back(mom->mass()); + ele_gen_MomPt2_.push_back(mom->pt()); + ele_gen_MomEta2_.push_back(mom->eta()); + ele_gen_MomPhi2_.push_back(mom->phi()); + ele_gen_MomBeta2_.push_back(beta_mom); + ele_gen_MomGamma2_.push_back(gamma_mom); + ele_gen_MomLxyz2_.push_back(std::sqrt(dx*dx + dy*dy + dz*dz)); + ele_gen_MomLz2_.push_back(dz); + ele_gen_MomLxy2_.push_back(std::sqrt(dx*dx + dy*dy)); + ele_gen_MomCTau02_.push_back(std::sqrt(dx*dx + dy*dy + dz*dz) / (beta_mom * gamma_mom)); + } + // final state hadrons + if( prt.status() == 1 && + prt.mother()->pdgId() != 9900014 && + prt.mother()->pdgId() != 13 && + prt.mother()->pdgId() != 11 ){ + + had_gen_PID_.push_back(prt.pdgId() ); + had_gen_Status_.push_back(prt.status()); + had_gen_Pt_.push_back(prt.pt()); + had_gen_Eta_.push_back(prt.eta()); + had_gen_Phi_.push_back(prt.phi()); + had_gen_Charge_.push_back(prt.charge()); + had_gen_Mass_.push_back(prt.charge()); + } + // quarks @ sv + if( prt.status() == 23 && + prt.mother()->pdgId() == 9900014 && + (abs(prt.pdgId()) != 13 || abs(prt.pdgId()) != 11) ){ + quarks_gen_PID_.push_back(prt.pdgId() ); + quarks_gen_Status_.push_back(prt.status()); + quarks_gen_Pt_.push_back(prt.pt()); + quarks_gen_Eta_.push_back(prt.eta()); + quarks_gen_Phi_.push_back(prt.phi()); + quarks_gen_Charge_.push_back(prt.charge()); + quarks_gen_Mass_.push_back(prt.charge()); + } +} + +void BigNtuple::set_pvInfo(TTree* tree){ + tree->Branch("pvX" , &pvX_, "pvX/F"); + tree->Branch("pvY" , &pvY_, "pvY/F"); + tree->Branch("pvZ" , &pvZ_, "pvZ/F"); + tree->Branch("pvXErr" , &pvXErr_, "pvXErr/F"); + tree->Branch("pvYErr" , &pvYErr_, "pvYErr/F"); + tree->Branch("pvZErr" , &pvZErr_, "pvZErr/F"); + tree->Branch("pvMass" , &pvMass_, "pvMass/F"); + tree->Branch("pvLxy" , &pvLxy_, "pvLxy/F"); + tree->Branch("pvLxyz" , &pvLxyz_, "pvLxyz/F"); + tree->Branch("pvLxySig" , &pvLxySig_, "pvLxySig/F"); + tree->Branch("pvLxyzSig" , &pvLxyzSig_, "pvLxyzSig/F"); + tree->Branch("pvChi2" , &pvChi2_, "pvChi2/F"); + tree->Branch("pvNTrack" , &pvNTrack_, "pvNTrack/i"); + tree->Branch("pvSumPtSq" , &pvSumPtSq_, "pvSumPtSq/F"); + tree->Branch("numberPV" , &numberPV_, "numberPV/i"); + +} + +void BigNtuple::fill_pvInfo(const reco::VertexCollection& pvs){ + numberPV_ = pvs.size(); + const reco::Vertex& pv = pvs.front(); + + float x = pv.x(), y = pv.y(), z = pv.z(); + float xE = pv.xError(), yE = pv.yError(), zE = pv.zError(); + + pvX_ = x; + pvY_ = y; + pvZ_ = z; + pvXErr_ = xE; + pvYErr_ = yE; + pvZErr_ = zE; + pvMass_ = 0; + pvLxy_ = std::sqrt( x * x + y * y ); + pvLxyz_ = std::sqrt( x * x + y * y + z * z); + pvLxySig_ = std::sqrt( x * x + y * y ) / std::sqrt(xE * xE + yE * yE); + pvLxyzSig_ = std::sqrt( x * x + y * y + z * z )/ std::sqrt(xE * xE + yE * yE + zE * zE); + pvChi2_ = pv.chi2(); + + reco::Vertex::trackRef_iterator vtxIter = pv.tracks_begin(); + float SumPtSq = 0; + int NTrack = 0; + for(; vtxIter != pv.tracks_end(); ++vtxIter) { + NTrack++; + SumPtSq += (*vtxIter)->pt() * (*vtxIter)->pt(); + } + pvNTrack_ = NTrack; + pvSumPtSq_ = SumPtSq; +} + + +void BigNtuple::set_trigInfo(TTree* tree){ + + tree->Branch("passIsoMuTk18" , &passIsoMuTk18_, "passIsoMuTk18/O"); + tree->Branch("passIsoMuTk20" , &passIsoMuTk20_, "passIsoMuTk20/O"); + tree->Branch("passIsoMuTk22" , &passIsoMuTk22_, "passIsoMuTk22/O"); + tree->Branch("passIsoMuTk24" , &passIsoMuTk24_, "passIsoMuTk24/O"); + tree->Branch("passIsoMuTk27" , &passIsoMuTk27_, "passIsoMuTk27/O"); + tree->Branch("passIsoMuTk17e" , &passIsoMuTk17e_, "passIsoMuTk17e/O"); + tree->Branch("passIsoMuTk22e" , &passIsoMuTk22e_, "passIsoMuTk22e/O"); + tree->Branch("passIsoMu18" , &passIsoMu18_, "passIsoMu18/O"); + tree->Branch("passIsoMu20" , &passIsoMu20_, "passIsoMu20/O"); + tree->Branch("passIsoMu22" , &passIsoMu22_, "passIsoMu22/O"); + tree->Branch("passIsoMu24" , &passIsoMu24_, "passIsoMu24/O"); + tree->Branch("passIsoMu27" , &passIsoMu27_, "passIsoMu27/O"); + tree->Branch("passIsoMu17e" , &passIsoMu17e_, "passIsoMu17e/O"); + tree->Branch("passIsoMu22e" , &passIsoMu22e_, "passIsoMu22e/O"); + tree->Branch("passTkMu17" , &passTkMu17_, "passTkMu17/O"); + tree->Branch("passTkMu20" , &passTkMu20_, "passTkMu20/O"); + tree->Branch("passIsoMu24All" , &passIsoMu24All_, "passIsoMu24All/O"); + tree->Branch("passIsoMu27All" , &passIsoMu27All_, "passIsoMu27All/O"); + tree->Branch("passDoubleMu17TrkIsoMu8" , &passDoubleMu17TrkIsoMu8_, "passDoubleMu17TrkIsoMu8/O"); + tree->Branch("passDoubleMu17TrkIsoTkMu8" , &passDoubleMu17TrkIsoTkMu8_, "passDoubleMu17TrkIsoTkMu8/O"); + tree->Branch("passDoubleTkMu17TrkIsoTkMu8" , &passDoubleTkMu17TrkIsoTkMu8_, "passDoubleTkMu17TrkIsoTkMu8/O"); + tree->Branch("passIsoEle27", &passIsoEle27_,"passIsoEle27/O"); + tree->Branch("passNonIsoEle115", &passNonIsoEle115_,"passNonIsoEle115/O"); + tree->Branch("passDoubleEle23andEle12DZ", &passDoubleEle23andEle12DZ_,"passDoubleEle23andEle12DZ/O"); + tree->Branch("passDoubleEle23andEle12", &passDoubleEle23andEle12_,"passDoubleEle23andEle12/O"); + tree->Branch("passDoubleEle33TrkMW", &passDoubleEle33TrkMW_,"passDoubleEle33TrkMW/O"); + tree->Branch("passDoubleEle33MW", &passDoubleEle33MW_,"passDoubleEle33MW/O"); + tree->Branch("passDoubleEle33", &passDoubleEle33_,"passDoubleEle33/O"); + tree->Branch("passDoubleMu33Ele33", &passDoubleMu33Ele33_,"passDoubleMu33Ele33/O"); + +} + +void BigNtuple::fill_trigInfo(const edm::TriggerResults& triggerResults, const edm::TriggerNames& trigNames){ + + for (size_t i = 0; i < trigNames.size(); ++i) { + const std::string &name = trigNames.triggerName(i); + bool fired = triggerResults.accept(i); + if(!fired) continue; + + passIsoMuTk18_ |= name.find("HLT_IsoTkMu18_v") != std::string::npos; + passIsoMuTk20_ |= name.find("HLT_IsoTkMu20_v") != std::string::npos; + passIsoMuTk22_ |= name.find("HLT_IsoTkMu22_v") != std::string::npos; + passIsoMuTk24_ |= name.find("HLT_IsoTkMu24_v") != std::string::npos; + passIsoMuTk27_ |= name.find("HLT_IsoTkMu27_v") != std::string::npos; + passIsoMuTk17e_ |= name.find("HLT_IsoTkMu17_eta2p1_v") != std::string::npos; + passIsoMuTk22e_ |= name.find("HLT_IsoTkMu22_eta2p1_v") != std::string::npos; + passIsoMu18_ |= name.find("HLT_IsoMu18_v") != std::string::npos; + passIsoMu20_ |= name.find("HLT_IsoMu20_v") != std::string::npos; + passIsoMu22_ |= name.find("HLT_IsoMu22_v") != std::string::npos; + passIsoMu24_ |= name.find("HLT_IsoMu24_v") != std::string::npos; + passIsoMu27_ |= name.find("HLT_IsoMu27_v") != std::string::npos; + + passIsoMu17e_ |= name.find("HLT_IsoTkMu17_eta2p1_v") != std::string::npos; + passIsoMu22e_ |= name.find("HLT_IsoTkMu22_eta2p1_v") != std::string::npos; + + passTkMu17_ |= name.find("HLT_TkMu17_v") != std::string::npos; + passTkMu20_ |= name.find("HLT_TkMu20_v") != std::string::npos; + + passDoubleMu17TrkIsoMu8_ |= name.find("HLT_Mu17_TrkIsoVVL_Mu8_TrkIsoVVL_DZ_v") != std::string::npos; + passDoubleMu17TrkIsoTkMu8_ |= name.find("HLT_Mu17_TrkIsoVVL_TkMu8_TrkIsoVVL_DZ_v") != std::string::npos; + passDoubleTkMu17TrkIsoTkMu8_ |= name.find("HLT_TkMu17_TrkIsoVVL_TkMu8_TrkIsoVVL_DZ_v") != std::string::npos; + + passIsoEle27_ |= name.find("HLT_Ele27_WPTight_Gsf_v") != std::string::npos; + passNonIsoEle115_ |= name.find("HLT_Ele115_CaloIdVT_GsfTrkIdT_v") != std::string::npos; + + passDoubleEle23andEle12DZ_ |= name.find("HLT_Ele23_Ele12_CaloIdL_TrackIdL_IsoVL_DZ_v") != std::string::npos; + passDoubleEle23andEle12_ |= name.find("HLT_Ele23_Ele12_CaloIdL_TrackIdL_IsoVL_v") != std::string::npos; + + passDoubleEle33TrkMW_ |= name.find("HLT_DoubleEle33_CaloIdL_GsfTrkIdVL_MW_v") != std::string::npos; + passDoubleEle33MW_ |= name.find("HLT_DoubleEle33_CaloIdL_MW_v") != std::string::npos; + passDoubleEle33_ |= name.find("HLT_DoubleEle33_CaloIdL_v") != std::string::npos; + + passDoubleMu33Ele33_ |= name.find("HLT_Mu33_Ele33_CaloIdL_GsfTrkIdVL_v") != std::string::npos; + + } + + passIsoMu24All_ = passIsoMu24All_ || passIsoMu24_ || passIsoMuTk24_ ; + passIsoMu27All_ = passIsoMu27All_ || passIsoMu27_ || passIsoMuTk27_ ; + +} + +void BigNtuple::set_pileupInfo(TTree* tree){ + + tree->Branch("npT" , &npT_, "npT/O"); + tree->Branch("npIT" , &npIT_, "npIT/O"); + tree->Branch("PU_Weight" , &pu_Weight_, "PU_Weight/O"); + tree->Branch("PU_WeightUp" , &pu_WeightUp_, "PU_WeightUp/O"); + tree->Branch("PU_WeightDown" , &pu_WeightDown_, "PU_WeightDown/O"); +} + +void BigNtuple::fill_pileupInfo( float npT, float npIT, float PU_Weight, float PU_WeightUp, float PU_WeightDown){ + + npT_ = npT; + npIT_ = npIT; + pu_Weight_ = PU_Weight; + pu_WeightUp_ = PU_WeightUp; + pu_WeightDown_ = PU_WeightDown; + +} + + + void BigNtuple::set_muInfo(TTree* tree){ + tree->Branch("mu_en" , &mu_en_); + tree->Branch("mu_pt" , &mu_pt_); + tree->Branch("mu_eta" , &mu_eta_); + tree->Branch("mu_phi" , &mu_phi_); + tree->Branch("mu_et" , &mu_et_); + tree->Branch("mu_rhoIso",&mu_rhoIso_); + tree->Branch("mu_charge" , &mu_charge_); + tree->Branch("mu_trackiso" , &mu_trackiso_); + tree->Branch("mu_pfSumChargedHadronPt" , &mu_pfSumChargedHadronPt_); + tree->Branch("mu_pfSumNeutralHadronEt" , &mu_pfSumNeutralHadronEt_); + tree->Branch("mu_PFSumPhotonEt" , &mu_PFSumPhotonEt_); + tree->Branch("mu_pfSumPUPt" , &mu_pfSumPUPt_); + tree->Branch("mu_numberOfValidMuonHits" , &mu_numberOfValidMuonHits_); + tree->Branch("mu_emIso" , &mu_emIso_); + tree->Branch("mu_hadIso" , &mu_hadIso_); + tree->Branch("mu_normalizedChi2" , &mu_normalizedChi2_); + tree->Branch("mu_numberOfMatchedStations" , &mu_numberOfMatchedStations_); + tree->Branch("mu_numberOfValidPixelHits" , &mu_numberOfValidPixelHits_); + tree->Branch("mu_numberOftrackerLayersWithMeasurement" , &mu_numberOftrackerLayersWithMeasurement_); + tree->Branch("mu_numberOfpixelLayersWithMeasurement" , &mu_numberOfpixelLayersWithMeasurement_); + tree->Branch("mu_TrackQuality" , &mu_TrackQuality_); + tree->Branch("mu_InnerTrackQuality" , &mu_InnerTrackQuality_); + tree->Branch("mu_pxTunePMuonBestTrack" , &mu_pxTunePMuonBestTrack_); + tree->Branch("mu_pyTunePMuonBestTrack" , &mu_pyTunePMuonBestTrack_); + tree->Branch("mu_pzTunePMuonBestTrack" , &mu_pzTunePMuonBestTrack_); + tree->Branch("mu_pTunePMuonBestTrack" , &mu_pTunePMuonBestTrack_); + tree->Branch("mu_etaTunePMuonBestTrack" , &mu_etaTunePMuonBestTrack_); + tree->Branch("mu_LXYZ" , &mu_LXYZ_); + tree->Branch("mu_LXY" , &mu_LXY_); + tree->Branch("mu_ptTunePMuonBestTrack" , &mu_ptTunePMuonBestTrack_); + tree->Branch("mu_phiTunePMuonBestTrack" , &mu_phiTunePMuonBestTrack_); + tree->Branch("mu_thetaTunePMuonBestTrack" , &mu_thetaTunePMuonBestTrack_); + tree->Branch("mu_chargeTunePMuonBestTrack" , &mu_chargeTunePMuonBestTrack_); + tree->Branch("mu_dPToverPTTunePMuonBestTrack" , &mu_dPToverPTTunePMuonBestTrack_); + tree->Branch("mu_absdxyTunePMuonBestTrack" , &mu_absdxyTunePMuonBestTrack_); + tree->Branch("mu_absdxyErrorTunePMuonBestTrack" , &mu_absdxyErrorTunePMuonBestTrack_); + tree->Branch("mu_absdxySigTunePMuonBestTrack" , &mu_absdxySigTunePMuonBestTrack_); + tree->Branch("mu_absdzTunePMuonBestTrack" , &mu_absdzTunePMuonBestTrack_); + tree->Branch("mu_absdzErrorTunePMuonBestTrack" , &mu_absdzErrorTunePMuonBestTrack_); + tree->Branch("mu_absdzSigTunePMuonBestTrack" , &mu_absdzSigTunePMuonBestTrack_); + tree->Branch("mu_recoDeltaBeta" , &mu_recoDeltaBeta_); + tree->Branch("mu_recoiso" , &mu_recoiso_); + tree->Branch("mu_isGlobalMuon" , &mu_isGlobalMuon_); + tree->Branch("mu_isStandAloneMuon" , &mu_isStandAloneMuon_); + tree->Branch("mu_isPF" , &mu_isPF_); + tree->Branch("mu_isRPCMuon" , &mu_isRPCMuon_); + tree->Branch("mu_isTrackerMuon" , &mu_isTrackerMuon_); + tree->Branch("mu_isGoodMuon" , &mu_isGoodMuon_); + tree->Branch("mu_isSoftMuon" , &mu_isSoftMuon_); + tree->Branch("mu_isLoose" , &mu_isLoose_); + tree->Branch("mu_isTightMuon" , &mu_isTightMuon_); + tree->Branch("mu_STAnHits" , &mu_STAnHits_); + tree->Branch("mu_STAnLost" , &mu_STAnLost_); + tree->Branch("mu_STAnStationsWithAnyHits" , &mu_STAnStationsWithAnyHits_); + tree->Branch("mu_STAnCscChambersWithAnyHits" , &mu_STAnCscChambersWithAnyHits_); + tree->Branch("mu_STAnDtChambersWithAnyHits" , &mu_STAnDtChambersWithAnyHits_); + tree->Branch("mu_STAnRpcChambersWithAnyHits" , &mu_STAnRpcChambersWithAnyHits_); + tree->Branch("mu_STAinnermostStationWithAnyHits" , &mu_STAinnermostStationWithAnyHits_); + tree->Branch("mu_STAoutermostStationWithAnyHits" , &mu_STAoutermostStationWithAnyHits_); + tree->Branch("mu_STAnStationsWithValidHits" , &mu_STAnStationsWithValidHits_); + tree->Branch("mu_STAnCscChambersWithValidHits" , &mu_STAnCscChambersWithValidHits_); + tree->Branch("mu_STAnDtChambersWithValidHit" , &mu_STAnDtChambersWithValidHit_); + tree->Branch("mu_STAnRpcChambersWithValidHits" , &mu_STAnRpcChambersWithValidHits_); + tree->Branch("mu_STAnValidMuonHits" , &mu_STAnValidMuonHits_); + tree->Branch("mu_STAnValidCscHits" , &mu_STAnValidCscHits_); + tree->Branch("mu_STAnValidDtHits" , &mu_STAnValidDtHits_); + tree->Branch("mu_STAnValidRpcHits" , &mu_STAnValidRpcHits_); + tree->Branch("mu_STAinnermostStationWithValidHits" , &mu_STAinnermostStationWithValidHits_); + tree->Branch("mu_STAoutermostStationWithValidHits" , &mu_STAoutermostStationWithValidHits_); + tree->Branch("mu_STATofDirection" , &mu_STATofDirection_); + tree->Branch("mu_STATofNDof" , &mu_STATofNDof_); + tree->Branch("mu_STATofTimeAtIpInOut" , &mu_STATofTimeAtIpInOut_); + tree->Branch("mu_STATofTimeAtIpInOutErr" , &mu_STATofTimeAtIpInOutErr_); + tree->Branch("mu_STATofTimeAtIpOutIn" , &mu_STATofTimeAtIpOutIn_); + tree->Branch("mu_STATofTimeAtIpOutInErr" , &mu_STATofTimeAtIpOutInErr_); + tree->Branch("mu_FirstGenMatch" , &mu_FirstGenMatch_); + tree->Branch("mu_SecondGenMatch" , &mu_SecondGenMatch_); + + } + + + +void BigNtuple::fill_muInfo(const pat::Muon& mu, const reco::Vertex& pv , double Rho ,double match1 , double match2){ + + mu_isGlobalMuon_.push_back(mu.isGlobalMuon()); + mu_isPF_.push_back(mu.isPFMuon()); + mu_isTrackerMuon_.push_back(mu.isTrackerMuon()); + mu_isRPCMuon_.push_back(mu.isRPCMuon()); + mu_isStandAloneMuon_.push_back(mu.isStandAloneMuon()); + mu_isSoftMuon_.push_back(mu.isSoftMuon(pv)); + mu_isLoose_.push_back(mu.isLooseMuon()); + mu_isTightMuon_.push_back(mu.isTightMuon(pv)); + mu_en_.push_back(mu.energy()); + mu_et_.push_back(mu.et()); + mu_pt_.push_back(mu.pt()); + mu_eta_.push_back(mu.eta()); + mu_phi_.push_back(mu.phi()); + mu_charge_.push_back(mu.charge()); + mu_FirstGenMatch_.push_back(match1); + mu_SecondGenMatch_.push_back(match2); + + reco::TrackRef tunePTrack = mu.muonBestTrack(); + + mu_ptTunePMuonBestTrack_.push_back(tunePTrack->pt()); // transverse momentum + mu_dPToverPTTunePMuonBestTrack_.push_back(tunePTrack->ptError()/tunePTrack->pt()); // error calculation of transverse momentum + mu_pxTunePMuonBestTrack_.push_back(tunePTrack->px()); //px component of the track + mu_pyTunePMuonBestTrack_.push_back(tunePTrack->py()); //py component of the track + mu_pzTunePMuonBestTrack_.push_back(tunePTrack->pz()); //pz component of the track + mu_pTunePMuonBestTrack_.push_back(tunePTrack->p()); //magnitude of momentum vector + mu_etaTunePMuonBestTrack_.push_back(tunePTrack->eta()); + mu_phiTunePMuonBestTrack_.push_back(tunePTrack->phi()); + mu_thetaTunePMuonBestTrack_.push_back(tunePTrack->theta()); + mu_chargeTunePMuonBestTrack_.push_back(tunePTrack->charge()); + mu_absdxyTunePMuonBestTrack_.push_back(fabs(tunePTrack->dxy(pv.position()))); //transvers impact parameter w.r.t. the primary vertex + mu_absdxyErrorTunePMuonBestTrack_.push_back(fabs(tunePTrack->dxyError())); //transvers impact parameter w.r.t. the primary vertex + mu_absdxySigTunePMuonBestTrack_.push_back(fabs(tunePTrack->dxy(pv.position()))/fabs(tunePTrack->dxyError())); + mu_absdzTunePMuonBestTrack_.push_back(fabs(tunePTrack->dz(pv.position()))); // longitudinal impact parameter w.r.t. the primary vertex + mu_absdzErrorTunePMuonBestTrack_.push_back(fabs(tunePTrack->dzError())); // longitudinal impact parameter w.r.t. the primary vertex + mu_absdzSigTunePMuonBestTrack_.push_back(fabs(tunePTrack->dz(pv.position()))/fabs(tunePTrack->dzError())); + mu_TrackQuality_.push_back(tunePTrack->quality(reco::TrackBase::highPurity)); + + mu_rhoIso_.push_back(Rho); //transverse momentum per unit area + + if(mu.globalTrack().isNonnull() ) { + mu_normalizedChi2_.push_back(mu.globalTrack()->normalizedChi2()); + mu_numberOfValidPixelHits_.push_back(mu.innerTrack()->hitPattern().numberOfValidPixelHits()); + mu_numberOfValidMuonHits_.push_back(mu.globalTrack()->hitPattern().numberOfValidMuonHits()); + mu_numberOftrackerLayersWithMeasurement_.push_back(mu.innerTrack()->hitPattern().trackerLayersWithMeasurement()); + mu_numberOfMatchedStations_.push_back(mu.numberOfMatchedStations()); + mu_numberOfpixelLayersWithMeasurement_.push_back(mu.innerTrack()->hitPattern().pixelLayersWithMeasurement()); + mu_InnerTrackQuality_.push_back(mu.innerTrack()->quality(reco::TrackBase::highPurity)); + } + else{ + mu_normalizedChi2_.push_back(-999); + mu_numberOfValidPixelHits_.push_back(-999); + mu_numberOfValidMuonHits_.push_back(-999); + mu_numberOftrackerLayersWithMeasurement_.push_back(-999); + mu_numberOfMatchedStations_.push_back(-999); + mu_numberOfpixelLayersWithMeasurement_.push_back(-999); + mu_InnerTrackQuality_.push_back(-999); + } + + if(mu.standAloneMuon().isNonnull() ) { + mu_STAnHits_.push_back(mu.standAloneMuon()->numberOfValidHits()); + mu_STAnLost_.push_back(mu.standAloneMuon()->numberOfLostHits()); + mu_STAnStationsWithAnyHits_.push_back(mu.standAloneMuon()->hitPattern().muonStationsWithAnyHits()); + mu_STAnCscChambersWithAnyHits_.push_back(mu.standAloneMuon()->hitPattern().cscStationsWithAnyHits()); //csc chambers in track fit + mu_STAnDtChambersWithAnyHits_.push_back(mu.standAloneMuon()->hitPattern().dtStationsWithAnyHits()); //dt chambers in track fit + mu_STAnRpcChambersWithAnyHits_.push_back(mu.standAloneMuon()->hitPattern().rpcStationsWithAnyHits()); //rpc chambers in track fit + mu_STAinnermostStationWithAnyHits_.push_back(mu.standAloneMuon()->hitPattern().innermostMuonStationWithAnyHits()); + mu_STAoutermostStationWithAnyHits_.push_back(mu.standAloneMuon()->hitPattern().outermostMuonStationWithAnyHits()); + mu_STAnCscChambersWithValidHits_.push_back(mu.standAloneMuon()->hitPattern().cscStationsWithValidHits()); //csc chambers anywhere near track + mu_STAnDtChambersWithValidHit_.push_back(mu.standAloneMuon()->hitPattern().dtStationsWithValidHits()); //dt chambers anywhere near track + mu_STAnRpcChambersWithValidHits_.push_back(mu.standAloneMuon()->hitPattern().rpcStationsWithValidHits()); //rpc chambers anywhere near track + mu_STAnValidCscHits_.push_back(mu.standAloneMuon()->hitPattern().numberOfValidMuonCSCHits()); //CSC hits anywhere near track + mu_STAnValidDtHits_.push_back(mu.standAloneMuon()->hitPattern().numberOfValidMuonDTHits()); //DT hits anywhere near track + mu_STAnValidRpcHits_.push_back(mu.standAloneMuon()->hitPattern().numberOfValidMuonRPCHits()); //RPC hits anywhere near track + mu_STAnValidMuonHits_.push_back(mu.standAloneMuon()->hitPattern().numberOfValidMuonHits()); //muon hits anywhere near track + mu_STAinnermostStationWithValidHits_.push_back(mu.standAloneMuon()->hitPattern().innermostMuonStationWithValidHits()); + mu_STAoutermostStationWithValidHits_.push_back(mu.standAloneMuon()->hitPattern().outermostMuonStationWithValidHits()); + mu_STAnStationsWithValidHits_.push_back(mu.standAloneMuon()->hitPattern().muonStationsWithValidHits()); + } + else{ + mu_STAnHits_.push_back(-999); + mu_STAnLost_.push_back(-999); + mu_STAnStationsWithAnyHits_.push_back(-999); + mu_STAnCscChambersWithAnyHits_.push_back(-999); + mu_STAnDtChambersWithAnyHits_.push_back(-999); + mu_STAnRpcChambersWithAnyHits_.push_back(-999); + mu_STAinnermostStationWithAnyHits_.push_back(-999); + mu_STAoutermostStationWithAnyHits_.push_back(-999); + mu_STAnCscChambersWithValidHits_.push_back(-999); + mu_STAnDtChambersWithValidHit_.push_back(-999); + mu_STAnRpcChambersWithValidHits_.push_back(-999); + mu_STAnValidCscHits_.push_back(-999); + mu_STAnValidDtHits_.push_back(-999); + mu_STAnValidRpcHits_.push_back(-999); + mu_STAnValidMuonHits_.push_back(-999); + mu_STAinnermostStationWithValidHits_.push_back(-999); + mu_STAoutermostStationWithValidHits_.push_back(-999); + mu_STAnStationsWithValidHits_.push_back(-999); + } + + reco::MuonTime tofAll = mu.time(); + mu_STATofDirection_.push_back(tofAll.direction()); + mu_STATofNDof_.push_back(tofAll.nDof); + mu_STATofTimeAtIpInOut_.push_back(tofAll.timeAtIpInOut); + mu_STATofTimeAtIpInOutErr_.push_back(tofAll.timeAtIpInOutErr); + mu_STATofTimeAtIpOutIn_.push_back(tofAll.timeAtIpOutIn); + mu_STATofTimeAtIpOutInErr_.push_back(tofAll.timeAtIpOutInErr); + //============= Parameters related to detector isolation ===================== + double charged = mu.pfIsolationR04().sumChargedHadronPt; + double neutral = mu.pfIsolationR04().sumNeutralHadronEt; + double pileup = mu.pfIsolationR04().sumPUPt; + double sumPhotonEt = mu.pfIsolationR04().sumPhotonEt; //Sum Et of PF photonds + double Mu_iso = 1.0*(charged + neutral + sumPhotonEt )/mu.pt(); //recommended this be < 0.20 (loose) or < 0.12 (tight) + double deltaBeta = (charged + std::max(0.0, neutral+sumPhotonEt-0.5*pileup))/mu.pt(); + mu_recoDeltaBeta_.push_back(deltaBeta); //Delta Beta + mu_recoiso_.push_back(Mu_iso); + mu_emIso_.push_back(mu.isolationR03().emEt); + mu_hadIso_.push_back(mu.isolationR03().hadEt); + mu_trackiso_.push_back(mu.isolationR03().sumPt); + //============= Parameters related to PF isolation ===================== + mu_pfSumPUPt_.push_back(mu.pfIsolationR03().sumPhotonEt); + mu_PFSumPhotonEt_.push_back(mu.pfIsolationR03().sumPhotonEt); + mu_pfSumChargedHadronPt_.push_back(mu.pfIsolationR03().sumChargedHadronPt); + mu_pfSumNeutralHadronEt_.push_back(mu.pfIsolationR03().sumNeutralHadronEt); + +} + +void BigNtuple::set_svInfo(TTree* tree){ + + tree->Branch("sv_mu_TrackSize" , &sv_mu_TrackSize_); + tree->Branch("sv_mu_LXYSig" , &sv_mu_LXYSig_); + tree->Branch("sv_mu_LXYZSig" , &sv_mu_LXYZSig_); + tree->Branch("sv_mu_LXY" , &sv_mu_LXY_); + tree->Branch("sv_mu_LXYZ" , &sv_mu_LXYZ_); + tree->Branch("sv_mu_mass" , &sv_mu_mass_); + tree->Branch("sv_mu_charge" , &sv_mu_charge_); + tree->Branch("sv_mu_eta" , &sv_mu_eta_); + tree->Branch("sv_mu_phi" , &sv_mu_phi_); + tree->Branch("sv_mu_pt" , &sv_mu_pt_); + tree->Branch("sv_mu_p" , &sv_mu_p_); + tree->Branch("sv_mu_Beta" , &sv_mu_Beta_); + tree->Branch("sv_mu_Gamma" , &sv_mu_Gamma_); + tree->Branch("sv_mu_CTau0" , &sv_mu_CTau0_); + tree->Branch("sv_mu_NDof" , &sv_mu_NDof_); + tree->Branch("sv_mu_Chi2" , &sv_mu_Chi2_); + tree->Branch("sv_mu_Angle3D" , &sv_mu_Angle3D_); + tree->Branch("sv_mu_Angle2D" , &sv_mu_Angle2D_); + tree->Branch("sv_mu_tracks_charge" , &sv_mu_tracks_charge_); + tree->Branch("sv_mu_tracks_eta" , &sv_mu_tracks_eta_); + tree->Branch("sv_mu_tracks_phi" , &sv_mu_tracks_phi_); + tree->Branch("sv_mu_tracks_pt" , &sv_mu_tracks_pt_); + tree->Branch("sv_mu_tracks_dxySig" , &sv_mu_tracks_dxySig_); + tree->Branch("sv_mu_tracks_dxy" , &sv_mu_tracks_dxy_); + tree->Branch("sv_mu_tracks_dxyz" , &sv_mu_tracks_dxyz_); + tree->Branch("sv_mu_tracks_Sumcharge" , &sv_mu_tracks_Sumcharge_); + tree->Branch("sv_mu_tracks_Sumpt" , &sv_mu_tracks_Sumpt_); + tree->Branch("sv_mu_match" , &sv_mu_match_); + + tree->Branch("sv_ele_TrackSize" , &sv_ele_TrackSize_); + tree->Branch("sv_ele_LXYSig" , &sv_ele_LXYSig_); + tree->Branch("sv_ele_LXYZSig" , &sv_ele_LXYZSig_); + tree->Branch("sv_ele_LXY" , &sv_ele_LXY_); + tree->Branch("sv_ele_LXYZ" , &sv_ele_LXYZ_); + tree->Branch("sv_ele_mass" , &sv_ele_mass_); + tree->Branch("sv_ele_charge" , &sv_ele_charge_); + tree->Branch("sv_ele_eta" , &sv_ele_eta_); + tree->Branch("sv_ele_phi" , &sv_ele_phi_); + tree->Branch("sv_ele_pt" , &sv_ele_pt_); + tree->Branch("sv_ele_p" , &sv_ele_p_); + tree->Branch("sv_ele_Beta" , &sv_ele_Beta_); + tree->Branch("sv_ele_Gamma" , &sv_ele_Gamma_); + tree->Branch("sv_ele_CTau0" , &sv_ele_CTau0_); + tree->Branch("sv_ele_NDof" , &sv_ele_NDof_); + tree->Branch("sv_ele_Chi2" , &sv_ele_Chi2_); + tree->Branch("sv_ele_Angle3D" , &sv_ele_Angle3D_); + tree->Branch("sv_ele_Angle2D" , &sv_ele_Angle2D_); + tree->Branch("sv_ele_tracks_charge" , &sv_ele_tracks_charge_); + tree->Branch("sv_ele_tracks_eta" , &sv_ele_tracks_eta_); + tree->Branch("sv_ele_tracks_phi" , &sv_ele_tracks_phi_); + tree->Branch("sv_ele_tracks_pt" , &sv_ele_tracks_pt_); + tree->Branch("sv_ele_tracks_dxySig" , &sv_ele_tracks_dxySig_); + tree->Branch("sv_ele_tracks_dxy" , &sv_ele_tracks_dxy_); + tree->Branch("sv_ele_tracks_dxyz" , &sv_ele_tracks_dxyz_); + tree->Branch("sv_ele_tracks_Sumcharge" , &sv_ele_tracks_Sumcharge_); + tree->Branch("sv_ele_tracks_Sumpt" , &sv_ele_tracks_Sumpt_); + tree->Branch("sv_ele_match" , &sv_ele_match_); + tree->Branch("sv_ele_score" , &sv_ele_score_); + +} + + +void BigNtuple::fill_sv_mu_Info(const reco::Vertex& bestVertex, const reco::Vertex& pv , double match){ + + float svChi2 = bestVertex.chi2(); + float svNDof = bestVertex.ndof(); + + //flight distance from the firstPV + float x = bestVertex.x(), y = bestVertex.y(), z = bestVertex.z(); + float dx = x - pv.x() , dy = y - pv.y(), dz = z - pv.z(); + + //build the total error + float svxE = bestVertex.xError(), svyE = bestVertex.yError(), svzE = bestVertex.zError(); + float pvxE = pv.xError(), pvyE = pv.yError(), pvzE = pv.zError(); + float xE = std::sqrt(svxE * svxE + pvxE * pvxE), yE = std::sqrt(svyE * svyE + pvyE * pvyE), zE = std::sqrt(svzE * svzE + pvzE * pvzE); + + // mother beta, gamma, ctau + float beta_mom = bestVertex.p4().P() / bestVertex.p4().energy(); + float gamma_mom = bestVertex.p4().energy() / bestVertex.p4().mass(); + + TVector3 pvVector3D(pv.x(), pv.y(), pv.z()); + TVector3 pvVector2D(pv.x(), pv.y(), 0); + TVector3 svVector3D(x, y, z); + TVector3 svVector2D(x, y, 0); + + // line pointing form the primary vertex through the sceondary vertex + TVector3 svMom3D( bestVertex.p4().x(), bestVertex.p4().y(), bestVertex.p4().z()); + TVector3 svMom2D( bestVertex.p4().x(), bestVertex.p4().y(), 0); + + // you want the negative when the momentum and sv are in the same + // direction relative to the PV + // this makes sure the angle is not pi when the vertex is fit behind + // the primary vertex + + float sign2D = (svMom2D * (svVector2D - pvVector2D)) > 0 ? -1: 1; + float sign3D = (svMom3D * (svVector3D - pvVector3D)) > 0 ? -1: 1; + + TVector3 pvToVertex3D( sign3D * dx, sign3D * dy, sign3D * dz); + TVector3 pvToVertex2D( sign2D * dx, sign2D * dy, 0); + + float svAngle3D = pvToVertex3D.Angle(svMom3D); + float svAngle2D = pvToVertex2D.Angle(svMom2D); + + sv_mu_TrackSize_.push_back(bestVertex.nTracks()); + sv_mu_LXY_.push_back(std::sqrt( dx * dx + dy * dy )); + sv_mu_LXYZ_.push_back(std::sqrt( dx * dx + dy * dy + dz * dz )); + sv_mu_LXYSig_.push_back(std::sqrt( dx * dx + dy * dy ) / std::sqrt(xE * xE + yE * yE)); + sv_mu_LXYZSig_.push_back(std::sqrt( dx * dx + dy * dy + dz * dz) / std::sqrt(xE * xE + yE * yE + zE * zE)); + sv_mu_mass_.push_back(bestVertex.p4().mass()); + sv_mu_eta_.push_back(bestVertex.p4().eta()); + sv_mu_phi_.push_back(bestVertex.p4().phi()); + sv_mu_pt_.push_back(bestVertex.p4().pt()); + sv_mu_p_.push_back(bestVertex.p4().P()); + sv_mu_Beta_.push_back(beta_mom); + sv_mu_Gamma_.push_back(gamma_mom); + sv_mu_CTau0_.push_back(std::sqrt( dx * dx + dy * dy + dz * dz) / (beta_mom * gamma_mom)); + sv_mu_NDof_.push_back(svNDof); + sv_mu_Chi2_.push_back(svChi2); + sv_mu_Angle3D_.push_back(svAngle3D); + sv_mu_Angle2D_.push_back(svAngle2D); + + int ch = 0; + float pt = 0; + + reco::Vertex::trackRef_iterator tt = bestVertex.tracks_begin(); + for(; tt != bestVertex.tracks_end(); ++tt) { + + sv_mu_tracks_charge_.push_back((*tt)->charge()); + sv_mu_tracks_eta_.push_back((*tt)->eta()); + sv_mu_tracks_phi_.push_back((*tt)->phi()); + sv_mu_tracks_pt_.push_back((*tt)->pt()); + sv_mu_tracks_dxySig_.push_back(fabs((*tt)->dxy(pv.position()))/fabs((*tt)->dxyError())); + sv_mu_tracks_dxy_.push_back((*tt)->dxy(pv.position())); + + ROOT::Math::SVector lxyz1((*tt)->vx()-pv.position().x(), (*tt)->vy()-pv.position().y(), (*tt)->vz()-pv.position().z()); + float dxyz = (float)ROOT::Math::Mag(lxyz1); // magntude of the vector + + sv_mu_tracks_dxyz_.push_back(dxyz); + ch+=(*tt)->charge(); + pt+=(*tt)->pt(); + } + + sv_mu_tracks_Sumcharge_.push_back(ch); + sv_mu_tracks_Sumpt_.push_back(pt); + sv_mu_match_.push_back(match); + +} + +void BigNtuple::fill_sv_ele_Info(const reco::Vertex& bestVertex, const reco::Vertex& pv , double match, double score){ + + float svChi2 = bestVertex.chi2(); + float svNDof = bestVertex.ndof(); + + //flight distance from the firstPV + float x = bestVertex.x(), y = bestVertex.y(), z = bestVertex.z(); + float dx = x - pv.x() , dy = y - pv.y(), dz = z - pv.z(); + + //build the total error + float svxE = bestVertex.xError(), svyE = bestVertex.yError(), svzE = bestVertex.zError(); + float pvxE = pv.xError(), pvyE = pv.yError(), pvzE = pv.zError(); + float xE = std::sqrt(svxE * svxE + pvxE * pvxE), yE = std::sqrt(svyE * svyE + pvyE * pvyE), zE = std::sqrt(svzE * svzE + pvzE * pvzE); + + // mother beta, gamma, ctau + float beta_mom = bestVertex.p4().P() / bestVertex.p4().energy(); + float gamma_mom = bestVertex.p4().energy() / bestVertex.p4().mass(); + + TVector3 pvVector3D(pv.x(), pv.y(), pv.z()); + TVector3 pvVector2D(pv.x(), pv.y(), 0); + TVector3 svVector3D(x, y, z); + TVector3 svVector2D(x, y, 0); + + // line pointing form the primary vertex through the sceondary vertex + TVector3 svMom3D( bestVertex.p4().x(), bestVertex.p4().y(), bestVertex.p4().z()); + TVector3 svMom2D( bestVertex.p4().x(), bestVertex.p4().y(), 0); + + // you want the negative when the momentum and sv are in the same + // direction relative to the PV + // this makes sure the angle is not pi when the vertex is fit behind + // the primary vertex + float sign2D = (svMom2D * (svVector2D - pvVector2D)) > 0 ? -1: 1; + float sign3D = (svMom3D * (svVector3D - pvVector3D)) > 0 ? -1: 1; + + TVector3 pvToVertex3D( sign3D * dx, sign3D * dy, sign3D * dz); + TVector3 pvToVertex2D( sign2D * dx, sign2D * dy, 0); + + float svAngle3D = pvToVertex3D.Angle(svMom3D); + float svAngle2D = pvToVertex2D.Angle(svMom2D); + + sv_ele_TrackSize_.push_back(bestVertex.nTracks()); + sv_ele_LXY_.push_back(std::sqrt( dx * dx + dy * dy )); + sv_ele_LXYZ_.push_back(std::sqrt( dx * dx + dy * dy + dz * dz )); + sv_ele_LXYSig_.push_back(std::sqrt( dx * dx + dy * dy ) / std::sqrt(xE * xE + yE * yE)); + sv_ele_LXYZSig_.push_back(std::sqrt( dx * dx + dy * dy + dz * dz) / std::sqrt(xE * xE + yE * yE + zE * zE)); + sv_ele_mass_.push_back(bestVertex.p4().mass()); + sv_ele_eta_.push_back(bestVertex.p4().eta()); + sv_ele_phi_.push_back(bestVertex.p4().phi()); + sv_ele_pt_.push_back(bestVertex.p4().pt()); + sv_ele_p_.push_back(bestVertex.p4().P()); + sv_ele_Beta_.push_back(beta_mom); + sv_ele_Gamma_.push_back(gamma_mom); + sv_ele_CTau0_.push_back(std::sqrt( dx * dx + dy * dy + dz * dz) / (beta_mom * gamma_mom)); + sv_ele_NDof_.push_back(svNDof); + sv_ele_Chi2_.push_back(svChi2); + sv_ele_Angle3D_.push_back(svAngle3D); + sv_ele_Angle2D_.push_back(svAngle2D); + sv_ele_score_.push_back(score); + + int ch = 0; + float pt = 0; + + reco::Vertex::trackRef_iterator tt = bestVertex.tracks_begin(); + for(; tt != bestVertex.tracks_end(); ++tt) { + + sv_ele_tracks_charge_.push_back((*tt)->charge()); + sv_ele_tracks_eta_.push_back((*tt)->eta()); + sv_ele_tracks_phi_.push_back((*tt)->phi()); + sv_ele_tracks_pt_.push_back((*tt)->pt()); + sv_ele_tracks_dxySig_.push_back(fabs((*tt)->dxy(pv.position()))/fabs((*tt)->dxyError())); + sv_ele_tracks_dxy_.push_back((*tt)->dxy(pv.position())); + + ROOT::Math::SVector lxyz1((*tt)->vx()-pv.position().x(), (*tt)->vy()-pv.position().y(), (*tt)->vz()-pv.position().z()); + float dxyz = (float)ROOT::Math::Mag(lxyz1); // magntude of the vector + + sv_ele_tracks_dxyz_.push_back(dxyz); + ch+=(*tt)->charge(); + pt+=(*tt)->pt(); + } + + sv_ele_tracks_Sumcharge_.push_back(ch); + sv_ele_tracks_Sumpt_.push_back(pt); + sv_ele_match_.push_back(match); + +} + + +void BigNtuple::set_jetInfo(TTree* tree){ + + tree->Branch("jet_charge" , &jet_charge_); + tree->Branch("jet_et" , &jet_et_); + tree->Branch("jet_pt" , &jet_pt_); + tree->Branch("jet_eta" , &jet_eta_); + tree->Branch("jet_phi" , &jet_phi_); + tree->Branch("jet_theta" , &jet_theta_); + tree->Branch("jet_en" , &jet_en_); + tree->Branch("jet_chargedEmEnergy" , &jet_chargedEmEnergy_); + tree->Branch("jet_neutralEmEnergyFraction" , &jet_neutralEmEnergyFraction_); + tree->Branch("jet_chargedHadronEnergy" , &jet_chargedHadronEnergy_); + tree->Branch("jet_neutralHadronEnergyFraction" , &jet_neutralHadronEnergyFraction_); + tree->Branch("jet_chargedMuEnergy" , &jet_chargedMuEnergy_); + tree->Branch("jet_chargedMuEnergyFraction" , &jet_chargedMuEnergyFraction_); + tree->Branch("jet_numberOfDaughters" , &jet_numberOfDaughters_); + tree->Branch("jet_muonEnergy" , &jet_muonEnergy_); + tree->Branch("jet_muonEnergyFraction" , &jet_muonEnergyFraction_); + tree->Branch("jet_muonMultiplicity" , &jet_muonMultiplicity_); + tree->Branch("jet_neutralEmEnergy" , &jet_neutralEmEnergy_); + tree->Branch("jet_neutralHadronEnergy" , &jet_neutralHadronEnergy_); + tree->Branch("jet_neutralHadronMultiplicity" , &jet_neutralHadronMultiplicity_); + tree->Branch("jet_neutralMultiplicity" , &jet_neutralMultiplicity_); + +} + + + +void BigNtuple::fill_jetInfo(const pat::Jet& jet){ + + jet_charge_.push_back(jet.charge()); + jet_et_.push_back(jet.et()); + jet_pt_.push_back(jet.pt()); + jet_eta_.push_back(jet.eta()); + jet_phi_.push_back(jet.phi()); + jet_theta_.push_back(jet.theta()); + jet_en_.push_back(jet.energy()); + jet_chargedEmEnergy_.push_back(jet.chargedEmEnergy()); + jet_neutralEmEnergyFraction_.push_back(jet.neutralEmEnergyFraction()); + jet_chargedHadronEnergy_.push_back(jet.chargedHadronEnergy()); + jet_neutralHadronEnergyFraction_.push_back(jet.neutralHadronEnergyFraction()); + jet_chargedMuEnergy_.push_back(jet.chargedMuEnergy()); + jet_chargedMuEnergyFraction_.push_back(jet.chargedMuEnergyFraction()); + jet_chargedMultiplicity_.push_back(jet.chargedMultiplicity()); + jet_numberOfDaughters_.push_back(jet.numberOfDaughters()); + jet_muonEnergy_.push_back(jet.muonEnergy()); + jet_muonEnergyFraction_.push_back(jet.muonEnergyFraction()); + jet_muonMultiplicity_.push_back(jet.muonMultiplicity()); + jet_neutralEmEnergy_.push_back(jet.neutralEmEnergy()); + jet_neutralHadronEnergy_.push_back(jet.neutralHadronEnergy()); + jet_neutralHadronMultiplicity_.push_back(jet.neutralHadronMultiplicity()); + jet_neutralMultiplicity_.push_back(jet.neutralMultiplicity()); + +} + +void BigNtuple::set_eleInfo(TTree* tree){ + + tree->Branch("ele_Et",&ele_Et_); + tree->Branch("ele_EtFromCaloEn",&ele_EtFromCaloEn_); + tree->Branch("ele_pt",&ele_pt_); + tree->Branch("ele_etaSC",&ele_etaSC_); + tree->Branch("ele_phiSC",&ele_phiSC_); + tree->Branch("ele_phiWidth",&ele_phiWidth_); + tree->Branch("ele_etaWidth",&ele_etaWidth_); + tree->Branch("ele_energySC",&ele_energySC_); + tree->Branch("ele_thetaSC",&ele_thetaSC_); + tree->Branch("ele_preshowerEnergySC",&ele_preshowerEnergySC_); + tree->Branch("ele_etaTrack",&ele_etaTrack_); + tree->Branch("ele_phiTrack",&ele_phiTrack_); + tree->Branch("ele_thetaTrack",&ele_thetaTrack_); + tree->Branch("ele_x",&ele_x_); + tree->Branch("ele_y",&ele_y_); + tree->Branch("ele_z",&ele_z_); + tree->Branch("ele_e2x5Max",&ele_e2x5Max_); + tree->Branch("ele_e1x5",&ele_e1x5_); + tree->Branch("ele_e5x5",&ele_e5x5_); + tree->Branch("ele_e2x5MaxOver5x5",&ele_e2x5MaxOver5x5_); + tree->Branch("ele_e1x5Over5x5",&ele_e1x5Over5x5_); + tree->Branch("ele_sigmaIetaIetaFull5x5",&ele_sigmaIetaIetaFull5x5_); + tree->Branch("ele_e2x5MaxFull5x5",&ele_e2x5MaxFull5x5_); + tree->Branch("ele_e1x5Full5x5",&ele_e1x5Full5x5_); + tree->Branch("ele_e5x5Full5x5",&ele_e5x5Full5x5_); + tree->Branch("ele_e2x5MaxOver5x5Full5x5",&ele_e2x5MaxOver5x5Full5x5_); + tree->Branch("ele_e1x5Over5x5Full5x5",&ele_e1x5Over5x5Full5x5_); + tree->Branch("ele_zTrackPositionAtVtx",&ele_zTrackPositionAtVtx_); + tree->Branch("ele_hadronicOverEm",&ele_hadronicOverEm_); + tree->Branch("ele_deltaEtaInSC",&ele_deltaEtaInSC_); + tree->Branch("ele_deltaPhiInSC",&ele_deltaPhiInSC_); + tree->Branch("ele_deltaEtaInSeedCluster",&ele_deltaEtaInSeedCluster_); + tree->Branch("ele_deltaPhiInSeedCluster",&ele_deltaPhiInSeedCluster_); + tree->Branch("ele_sigmaIetaIeta",&ele_sigmaIetaIeta_); + tree->Branch("ele_rawId",&ele_rawId_); + tree->Branch("ele_ieta",&ele_ieta_); + tree->Branch("ele_e2x5Right",&ele_e2x5Right_); + tree->Branch("ele_e2x5Left",&ele_e2x5Left_); + tree->Branch("ele_e2x5Top",&ele_e2x5Top_); + tree->Branch("ele_e2x5Bottom",&ele_e2x5Bottom_); + tree->Branch("ele_eMax",&ele_eMax_); + tree->Branch("ele_eRight",&ele_eRight_); + tree->Branch("ele_eLeft",&ele_eLeft_); + tree->Branch("ele_eTop",&ele_eTop_); + tree->Branch("ele_eBottom",&ele_eBottom_); + tree->Branch("ele_e3x3",&ele_e3x3_); + tree->Branch("ele_frac51",&ele_frac51_); + tree->Branch("ele_frac15",&ele_frac15_); + tree->Branch("ele_dxy",&ele_dxy_); + tree->Branch("ele_dz",&ele_dz_); + tree->Branch("ele_isEcalDrivenSeed",&ele_isEcalDrivenSeed_); + tree->Branch("ele_isPassConversionVeto",&ele_isPassConversionVeto_); + tree->Branch("ele_charge",&ele_charge_); + tree->Branch("ele_rhoIso",&ele_rhoIso_); + tree->Branch("ele_nbOfMissingHits",&ele_nbOfMissingHits_); + tree->Branch("ele_fbrem",&ele_fbrem_); + tree->Branch("ele_EoverP",&ele_EoverP_); + tree->Branch("ele_Xposition",&ele_Xposition_); + tree->Branch("ele_Yposition",&ele_Yposition_); + tree->Branch("ele_dr03TkSumPt",&ele_dr03TkSumPt_); + tree->Branch("ele_hcalDepth1OverEcal",&ele_hcalDepth1OverEcal_); + tree->Branch("ele_hcalDepth2OverEcal",&ele_hcalDepth2OverEcal_); + tree->Branch("ele_dr03HcalDepth2TowerSumEt",&ele_dr03HcalDepth2TowerSumEt_); + tree->Branch("ele_hcalDepth2TowerSumEtNoVeto",&ele_hcalDepth2TowerSumEtNoVeto_); + tree->Branch("ele_hcalDepth1TowerSumEtNoVeto",&ele_hcalDepth1TowerSumEtNoVeto_); + tree->Branch("ele_EcalPlusHcald1iso",&ele_EcalPlusHcald1iso_); + tree->Branch("ele_dr03EcalRecHitSumEt",&ele_dr03EcalRecHitSumEt_); + tree->Branch("ele_dr03HcalDepth1TowerSumEt",&ele_dr03HcalDepth1TowerSumEt_); + tree->Branch("ele_dr03HcalDepth1TowerSumEtBc",&ele_dr03HcalDepth1TowerSumEtBc_); + tree->Branch("ele_pfSumPhotonEt",&ele_pfSumPhotonEt_); + tree->Branch("ele_pfSumChargedHadronPt",&ele_pfSumChargedHadronPt_); + tree->Branch("ele_pfSumNeutralHadronEt",&ele_pfSumNeutralHadronEt_); + tree->Branch("ele_pfSumPUPt",&ele_pfSumPUPt_); + tree->Branch("ele_pfDeltaBeta",&ele_pfDeltaBeta_); + tree->Branch("ele_FirstGenMatch",&ele_FirstGenMatch_); + tree->Branch("ele_SecondGenMatch", &ele_SecondGenMatch_); +} + +void BigNtuple::fill_eleInfo(const pat::Electron& ele_, const reco::Vertex& pv, double Rho, double match1, double match2, std::auto_ptr recHitEcal){ + + + ele_FirstGenMatch_.push_back(match1); + ele_SecondGenMatch_.push_back(match2); + ele_Et_.push_back(ele_.superCluster()->energy() * sin(ele_.p4().theta())); + ele_EtFromCaloEn_.push_back(ele_.caloEnergy() * sin(ele_.p4().theta())); + + ele_pt_.push_back(ele_.pt()); + ele_etaSC_.push_back(ele_.superCluster()->eta()); //eta SC + ele_phiSC_.push_back(ele_.superCluster()->phi()); //phi SC + ele_phiWidth_.push_back(ele_.superCluster()->phiWidth()); + ele_etaWidth_.push_back(ele_.superCluster()->etaWidth()); + ele_energySC_.push_back(ele_.superCluster()->energy()); //energy SC + ele_thetaSC_.push_back(ele_.caloPosition().theta()); //theta SC + ele_preshowerEnergySC_.push_back(ele_.superCluster()->preshowerEnergy()); + + ele_etaTrack_.push_back(ele_.p4().eta()); //eta track + ele_phiTrack_.push_back(ele_.p4().phi()); //phi track + ele_thetaTrack_.push_back(ele_.p4().theta()); //theta track + + ele_x_.push_back(ele_.p4().x()); + ele_y_.push_back(ele_.p4().y()); + ele_z_.push_back(ele_.p4().z()); + + ele_e2x5Max_.push_back(ele_.e2x5Max()); + ele_e1x5_.push_back(ele_.e1x5()); + ele_e5x5_.push_back(ele_.e5x5()); + ele_e2x5MaxOver5x5_.push_back(ele_.e2x5Max()/ele_.e5x5()); + ele_e1x5Over5x5_.push_back(ele_.e1x5()/ele_.e5x5()); + ele_sigmaIetaIetaFull5x5_.push_back(ele_.full5x5_sigmaIetaIeta()); + ele_e2x5MaxFull5x5_.push_back(ele_.full5x5_e2x5Max()); + ele_e1x5Full5x5_.push_back(ele_.full5x5_e1x5()); + ele_e5x5Full5x5_.push_back(ele_.full5x5_e5x5()); + ele_e2x5MaxOver5x5Full5x5_.push_back(ele_.full5x5_e2x5Max()/ele_.full5x5_e5x5()); + ele_e1x5Over5x5Full5x5_.push_back(ele_.full5x5_e1x5()/ele_.full5x5_e5x5()); + + ele_zTrackPositionAtVtx_.push_back(ele_.TrackPositionAtVtx().Z()); + ele_hadronicOverEm_.push_back(ele_.hadronicOverEm()); + ele_deltaEtaInSC_.push_back(ele_.deltaEtaSuperClusterTrackAtVtx()); + ele_deltaPhiInSC_.push_back(ele_.deltaPhiSuperClusterTrackAtVtx()); + ele_deltaEtaInSeedCluster_.push_back(ele_.deltaEtaSeedClusterTrackAtVtx()); + ele_deltaPhiInSeedCluster_.push_back(ele_.deltaPhiSeedClusterTrackAtCalo()); + ele_sigmaIetaIeta_.push_back(ele_.sigmaIetaIeta()); + + EBDetId BarrelId = ele_.superCluster()->seed()->seed(); + + ele_rawId_.push_back(BarrelId.rawId()); + ele_ieta_.push_back(BarrelId.ieta()); + + ele_e2x5Right_.push_back(recHitEcal->e2x5Right(*(ele_.superCluster()->seed()))); + ele_e2x5Left_.push_back(recHitEcal->e2x5Left(*(ele_.superCluster()->seed()))); + ele_e2x5Top_.push_back(recHitEcal->e2x5Top(*(ele_.superCluster()->seed()))); + ele_e2x5Bottom_.push_back(recHitEcal->e2x5Bottom(*(ele_.superCluster()->seed()))); + ele_eMax_.push_back(recHitEcal->eMax(*(ele_.superCluster()->seed()))); + ele_eRight_.push_back(recHitEcal->eRight(*(ele_.superCluster()->seed()))); + ele_eLeft_.push_back(recHitEcal->eLeft(*(ele_.superCluster()->seed()))); + ele_eTop_.push_back(recHitEcal->eTop(*(ele_.superCluster()->seed()))); + ele_eBottom_.push_back(recHitEcal->eBottom(*(ele_.superCluster()->seed()))); + ele_e3x3_.push_back(recHitEcal->e3x3(*(ele_.superCluster()->seed()))); + ele_frac51_.push_back( recHitEcal->e5x1(*(ele_.superCluster()->seed()))/ele_.full5x5_e5x5() ); + ele_frac15_.push_back( recHitEcal->e1x5(*(ele_.superCluster()->seed()))/ele_.full5x5_e5x5() ); + + ele_dxy_.push_back(ele_.gsfTrack()->dxy(pv.position())); //GSF -> Gaussian Sum Filter + ele_dz_.push_back(ele_.gsfTrack()->dz(pv.position())); + + ele_isEcalDrivenSeed_.push_back(ele_.ecalDrivenSeed()); + ele_isPassConversionVeto_.push_back(ele_.passConversionVeto()); + ele_charge_.push_back(ele_.gsfTrack()->charge()); + ele_rhoIso_.push_back(Rho); //transverse momentum per unit area + ele_nbOfMissingHits_.push_back(ele_.gsfTrack()->numberOfLostHits()); + ele_fbrem_.push_back(ele_.fbrem()); + ele_EoverP_.push_back(ele_.eSeedClusterOverP()); + ele_Xposition_.push_back(ele_.caloPosition().x()); + ele_Yposition_.push_back(ele_.caloPosition().y()); + + //tracker isolation + ele_dr03TkSumPt_.push_back(ele_.dr03TkSumPt()); + + //------------- detector isolation ------------------------- + ele_hcalDepth1OverEcal_.push_back(ele_.hcalDepth1OverEcal()); + ele_hcalDepth2OverEcal_.push_back(ele_.hcalDepth2OverEcal()); + ele_dr03HcalDepth2TowerSumEt_.push_back(ele_.dr03HcalDepth2TowerSumEt()); + ele_hcalDepth2TowerSumEtNoVeto_.push_back(ele_.isolationVariables03().hcalDepth2TowerSumEt);// hcaldepht2 iso deposit with + // electron footprint removed + ele_hcalDepth1TowerSumEtNoVeto_.push_back(ele_.isolationVariables03().hcalDepth1TowerSumEt);// hcaldepht1 iso deposit with + // electron footprint removed + ele_EcalPlusHcald1iso_.push_back(ele_.dr03EcalRecHitSumEt() + ele_.dr03HcalDepth1TowerSumEt()); + ele_dr03EcalRecHitSumEt_.push_back(ele_.dr03EcalRecHitSumEt()); + ele_dr03HcalDepth1TowerSumEt_.push_back(ele_.dr03HcalDepth1TowerSumEt()); + ele_dr03HcalDepth1TowerSumEtBc_.push_back(ele_.dr03HcalDepth1TowerSumEtBc()); + //------------- PF isolation from pat::electron ------------------------- + ele_pfSumPhotonEt_.push_back(ele_.pfIsolationVariables().sumPhotonEt); + ele_pfSumChargedHadronPt_.push_back(ele_.pfIsolationVariables().sumChargedHadronPt); + ele_pfSumNeutralHadronEt_.push_back(ele_.pfIsolationVariables().sumNeutralHadronEt); + ele_pfSumPUPt_.push_back(ele_.pfIsolationVariables().sumPUPt); + // deltaBeta + double charged = ele_.pfIsolationVariables().sumPhotonEt; + double neutral = ele_.pfIsolationVariables().sumNeutralHadronEt; + double pileup = ele_.pfIsolationVariables().sumPUPt; + double deltaBeta = charged + std::max(0.0, neutral-0.5*pileup); + ele_pfDeltaBeta_.push_back(deltaBeta); + +} +void BigNtuple::set_eleIDInfo(TTree* tree){ + + tree->Branch("ele_Mva2016" , &ele_Mva2016_); + tree->Branch("ele_CutVeto" , &ele_CutVeto_); + tree->Branch("ele_CutLoose" , &ele_CutLoose_); + tree->Branch("ele_CutMedium" , &ele_CutMedium_); + tree->Branch("ele_CutTight" , &ele_CutTight_); + +} + +void BigNtuple::fill_eleIDInfo(float ele_mva , bool ele_veto , bool ele_loose , bool ele_medium ,bool ele_tight){ + + ele_Mva2016_.push_back(ele_mva); + ele_CutVeto_.push_back(ele_veto); + ele_CutLoose_.push_back(ele_loose); + ele_CutMedium_.push_back(ele_medium); + ele_CutTight_.push_back(ele_tight); + +} + +void BigNtuple::set_metInfo(TTree* tree){ + + tree->Branch("pfMet_et" , &pfMet_et_, "pfMet_et/F"); + tree->Branch("pfMet_pt" , &pfMet_pt_, "pfMet_pt/F"); + tree->Branch("pfMet_phi" , &pfMet_phi_, "pfMet_phi/F"); + tree->Branch("pfMet_en" , &pfMet_en_, "pfMet_en/F"); + tree->Branch("pfMet_px" , &pfMet_px_, "pfMet_px/F"); + tree->Branch("pfMet_py" , &pfMet_py_, "pfMet_py/F"); + tree->Branch("pfMet_pz" , &pfMet_pz_, "pfMet_pz/F"); + tree->Branch("pfMet_sumEt" , &pfMet_sumEt_, "pfMet_sumEt/F"); + tree->Branch("caloMet_pt" , &caloMet_pt_, "caloMet_pt/F"); + tree->Branch("caloMet_phi" , &caloMet_phi_, "caloMet_phi/F"); + +} + + + +void BigNtuple::fill_metInfo(const pat::MET& met){ + + pfMet_et_ = met.et(); + pfMet_pt_ = met.pt(); + pfMet_phi_ = met.phi(); + pfMet_en_ = met.energy(); + pfMet_px_ = met.px(); + pfMet_py_ = met.py(); + pfMet_pz_ = met.pz(); + pfMet_sumEt_ = met.sumEt(); + + caloMet_pt_ = met.caloMETPt(); + caloMet_phi_ = met.caloMETPhi(); + + +} + +void BigNtuple::set_bjetInfo(TTree* tree){ + tree->Branch("jet_btag_pt",&jet_btag_pt_); + tree->Branch("jet_btag_eta",&jet_btag_eta_); + tree->Branch("jet_btag_phi",&jet_btag_phi_); + tree->Branch("jet_btag_flavor",&jet_btag_flavor_); + tree->Branch("jet_btag_pfCSVv2IVF_discriminator",&jet_btag_pfCSVv2IVF_discriminator_); +} + +void BigNtuple::fill_bjetInfo(const pat::Jet& jet, const std::string & bDiscr, int flavor){ + + jet_btag_pt_.push_back(jet.pt()); + jet_btag_eta_.push_back(jet.eta()); + jet_btag_phi_.push_back(jet.phi()); + jet_btag_flavor_.push_back(flavor); + jet_btag_pfCSVv2IVF_discriminator_.push_back(jet.bDiscriminator(bDiscr)); + +} diff --git a/HeavyNeutralLeptonAnalysis/test/HeavyNeutralLeptonAnalyzer_cfg.py b/HeavyNeutralLeptonAnalysis/test/HeavyNeutralLeptonAnalyzer_cfg.py index 0f37b99..cd7b285 100644 --- a/HeavyNeutralLeptonAnalysis/test/HeavyNeutralLeptonAnalyzer_cfg.py +++ b/HeavyNeutralLeptonAnalysis/test/HeavyNeutralLeptonAnalyzer_cfg.py @@ -1,49 +1,110 @@ from os import path as path import FWCore.ParameterSet.Config as cms +debugLevel = -1 + +isMC_ = True +isMCSignal_ = True + +jec_tag_DATA = 'JetCorrectorParametersCollection_Summer16_23Sep2016AllV4_DATA_AK4PFchs' +jec_tag_MC = 'JetCorrectorParametersCollection_Summer16_23Sep2016V4_MC_AK4PFchs' +jec_file_DATA = 'sqlite_file:Summer16_23Sep2016AllV4_DATA.db' +jec_file_MC = 'sqlite_file:Summer16_23Sep2016V4_MC.db' +jec_tag = jec_tag_MC if isMC_ else jec_tag_DATA +jec_file = jec_file_MC if isMC_ else jec_file_DATA +algorithm = "AK4PFchs" + process = cms.Process("AnalysisProc") process.load("FWCore.MessageService.MessageLogger_cfi") -process.MessageLogger.cerr.FwkReport.reportEvery = 10000 +process.load('Configuration.EventContent.EventContent_cff') +process.load('Configuration.StandardSequences.GeometryRecoDB_cff') +process.load('Configuration.StandardSequences.MagneticField_cff') +process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_condDBv2_cff') +process.load('Configuration.StandardSequences.Services_cff') +process.load("TrackingTools.TransientTrack.TransientTrackBuilder_cfi") + +process.MessageLogger.cerr.FwkReport.reportEvery = 1 -import PhysicsTools.PythonAnalysis.LumiList as LumiList +import FWCore.PythonUtilities.LumiList as LumiList LumiList.LumiList().getVLuminosityBlockRange() #from Configuration.AlCa.GlobalTag import GlobalTag -process.load('Configuration.StandardSequences.FrontierConditions_GlobalTag_condDBv2_cff') from Configuration.AlCa.GlobalTag_condDBv2 import GlobalTag -process.GlobalTag.globaltag = '' - +#process.GlobalTag.globaltag = '' +process.GlobalTag = GlobalTag(process.GlobalTag, '94X_mc2017_realistic_v10') +process.maxEvents = cms.untracked.PSet( + input = cms.untracked.int32(13) +) process.source = cms.Source("PoolSource", - fileNames = cms.untracked.vstring('root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_105.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_106.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_108.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_109.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_110.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_111.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_113.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_115.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_116.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_117.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_118.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_119.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_120.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_121.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_124.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_129.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_131.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_135.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_136.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_137.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_138.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_139.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_140.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_141.root', - 'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoAOD/Moriond17/displaced/HeavyNeutrino_lljj_M-2_V-0.00836660026534_mu_onshell_pre2017_leptonFirst_NLO/heavyNeutrino_142.root') ) - -process.TFileService = cms.Service("TFileService", fileName = "HeavyNeutralLepton.root") + fileNames = cms.untracked.vstring( +#'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoMiniAOD/Fall17/displaced/HeavyNeutrino_lljj_M-5_V-0.00836660026534_mu_massiveAndCKM_LO/heavyNeutrino_1.root' +#'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoMiniAOD/Fall17/displaced/HeavyNeutrino_lljj_M-5_V-0.00836660026534_mu_massiveAndCKM_LO/heavyNeutrino_1.root', +#'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoMiniAOD/Fall17/displaced/HeavyNeutrino_lljj_M-5_V-0.00836660026534_mu_massiveAndCKM_LO/heavyNeutrino_1.root', +#'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoMiniAOD/Fall17/displaced/HeavyNeutrino_lljj_M-5_V-0.00836660026534_mu_massiveAndCKM_LO/heavyNeutrino_1.root', +#'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoMiniAOD/Fall17/displaced/HeavyNeutrino_lljj_M-5_V-0.00836660026534_mu_massiveAndCKM_LO/heavyNeutrino_1.root', +#'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoMiniAOD/Fall17/displaced/HeavyNeutrino_lljj_M-5_V-0.00836660026534_mu_massiveAndCKM_LO/heavyNeutrino_1.root', +#'root://cms-xrd-global.cern.ch//store/user/tomc/heavyNeutrinoMiniAOD/Fall17/displaced/HeavyNeutrino_lljj_M-5_V-0.00836660026534_mu_massiveAndCKM_LO/heavyNeutrino_1.root' +#'root://cms-xrd-global.cern.ch//store/mc/RunIIFall17MiniAOD/DYJetsToLL_M-50_Zpt-150toInf_TuneCP5_13TeV-madgraphMLM-pythia8/MINIAODSIM/RECOSIMstep_94X_mc2017_realistic_v10-v1/50000/EE9CC3E0-0DED-E711-BCAC-00E081CB560C.root' +'root://cms-xrd-global.cern.ch//store/mc/RunIISummer16MiniAODv2/TTJets_DiLept_TuneCUETP8M1_13TeV-madgraphMLM-pythia8/MINIAODSIM/PUMoriond17_80X_mcRun2_asymptotic_2016_TrancheIV_v6_ext1-v1/80000/4E597432-24BE-E611-ACBB-00266CFFBFC0.root' +#'root://cms-xrd-global.cern.ch//store/data/Run2016G/SingleMuon/MINIAOD/23Sep2016-v1/1110000/72446D9C-D89C-E611-9060-002590A3C984.root' +)) +process.TFileService = cms.Service("TFileService", fileName = cms.string("HeavyNeutralLepton.root")) +process.load('HNL.DisplacedAdaptiveVertexFinder.displacedInclusiveVertexing_cff') + +from HNL.HeavyNeutralLeptonAnalysis.ele_Sequence_cff import addElectronSequence + +addElectronSequence(process) + +process.load("CondCore.CondDB.CondDB_cfi") +process.jec = cms.ESSource("PoolDBESSource", + DBParameters = cms.PSet( + messageLevel = cms.untracked.int32(0) + ), + timetype = cms.string('runnumber'), + toGet = cms.VPSet( + cms.PSet( + record = cms.string('JetCorrectionsRecord'), + tag = cms.string(jec_tag), + label = cms.untracked.string('AK4PFchs') + ), + ), + connect = cms.string(jec_file) +) +#make sure we use the db source and not the global tag +process.es_prefer_jec = cms.ESPrefer('PoolDBESSource','jec') +process.get = cms.EDAnalyzer("EventSetupRecordDataGetter", + toGet = cms.VPSet(cms.PSet( + record = cms.string('JetCorrectionsRecord'), + data = cms.vstring('JetCorrectorParametersCollection/AK4PFchs') + ) + ), + verbose = cms.untracked.bool(True) +) + + +process.load("PhysicsTools.PatAlgos.producersLayer1.jetUpdater_cff") + +process.jetCorrFactors = process.updatedPatJetCorrFactors.clone( + src = cms.InputTag("slimmedJets"), + levels = ['L1FastJet', + 'L2Relative', + 'L3Absolute', + 'L2L3Residual'], + payload = 'AK4PFchs' ) + +process.slimmedJetsJEC = process.updatedPatJets.clone( + jetSource = cms.InputTag("slimmedJets"), + jetCorrFactorsSource = cms.VInputTag(cms.InputTag("jetCorrFactors")) + ) + + +#process.ak4PFJetsCHSL1FastL2L3 = process.ak4PFCHSJetsL1.clone(correctors = ['ak4PFL1FastL2L3Corrector']) +#process.correctJets = cms.Sequence(process.ak4PFL1FastL2L3CorrectorChain * process.ak4PFJetsCHSL1FastL2L3) process.HeavyNeutralLepton = cms.EDAnalyzer('HeavyNeutralLeptonAnalysis', - isMC = cms.bool(True), + debugLevel = cms.int32(debugLevel), + isMC = cms.bool(isMC_), + isMCSignal = cms.bool(isMCSignal_), vtxSrc = cms.InputTag("offlineSlimmedPrimaryVertices"), rho = cms.InputTag("fixedGridRhoFastjetAll"), muonSrc = cms.InputTag("slimmedMuons"), @@ -52,7 +113,7 @@ recHitCollectionEESrc = cms.InputTag("reducedEgamma","reducedEERecHits"), tauSrc = cms.InputTag("slimmedTaus"), packCandSrc = cms.InputTag("packedPFCandidates"), - jetSrc = cms.InputTag("slimmedJets"), + jetSrc = cms.InputTag("slimmedJetsJEC"), pfMETSrc = cms.InputTag("slimmedMETs"), triggerResultSrc = cms.InputTag("TriggerResults","","HLT"), metFilterResultSrc = cms.InputTag("TriggerResults","","PAT"), @@ -60,6 +121,23 @@ genEventInfoProduct = cms.InputTag("generator"), PUInfo = cms.InputTag("slimmedAddPileupInfo"), lheEventProducts = cms.InputTag("externalLHEProducer"), + SecondaryVertices = cms.InputTag("displacedInclusiveSecondaryVertices"), + bDiscriminators = cms.vstring("pfCombinedInclusiveSecondaryVertexV2BJetTags"), + electronsMva = cms.InputTag("electronMVAValueMapProducer:ElectronMVAEstimatorRun2Spring16GeneralPurposeV1Values"), + electronsVeto = cms.InputTag("egmGsfElectronIDs:cutBasedElectronID-Summer16-80X-V1-veto"), + electronsLoose = cms.InputTag("egmGsfElectronIDs:cutBasedElectronID-Summer16-80X-V1-loose"), + electronsMedium= cms.InputTag("egmGsfElectronIDs:cutBasedElectronID-Summer16-80X-V1-medium"), + electronsTight = cms.InputTag("egmGsfElectronIDs:cutBasedElectronID-Summer16-80X-V1-tight") ) -process.p = cms.Path(process.HeavyNeutralLepton) +#process.HeavyNeutralLepton.outputFileName = cms.untracked.string('HeavyNeutralLepton.root') +#process.HeavyNeutralLepton.genTreeName = cms.untracked.string('Tree') + + +process.p = cms.Path( + process.displacedInclusiveVertexing + *process.ele_Sequence + *process.jetCorrFactors + *process.slimmedJetsJEC + *process.HeavyNeutralLepton + ) diff --git a/HeavyNeutralLeptonAnalysis/test/Summer16_23Sep2016AllV4_DATA.db b/HeavyNeutralLeptonAnalysis/test/Summer16_23Sep2016AllV4_DATA.db new file mode 100644 index 0000000..0874161 Binary files /dev/null and b/HeavyNeutralLeptonAnalysis/test/Summer16_23Sep2016AllV4_DATA.db differ diff --git a/HeavyNeutralLeptonAnalysis/test/Summer16_23Sep2016V4_MC.db b/HeavyNeutralLeptonAnalysis/test/Summer16_23Sep2016V4_MC.db new file mode 100644 index 0000000..7d27361 Binary files /dev/null and b/HeavyNeutralLeptonAnalysis/test/Summer16_23Sep2016V4_MC.db differ diff --git a/README.md b/README.md index de8d998..9150928 100644 --- a/README.md +++ b/README.md @@ -1,5 +1,5 @@ # HNL -Repository to analyzr Heavy Neutral Leptons in CMS +Repository to analyzer Heavy Neutral Leptons in CMS Instructions: @@ -7,6 +7,6 @@ Instructions: cmsrel CMSSW_9_4_0 cd CMSSW_9_4_0/src cmsenv -git clone git@github.com:mmusich/HNL.git +git clone git@github.com:HNLto2L2Q/HNL.git scramv1 b -j 8 ```