Analyzer/src/S8TreeMaker.cc

/**
 * S8TreeMaker
 * 
 *
 * Created by Samvel Khalatian on Sep 29, 2010
 * Copyright 2010, All rights reserved
 */

#include <iostream>
#include <vector>

#include <boost/algorithm/string.hpp>
#include <boost/lexical_cast.hpp>
#include <boost/logic/tribool.hpp>
#include <boost/regex.hpp>

#include <TTree.h>
#include <TLorentzVector.h>
#include <TVector3.h>

#include "CommonTools/UtilAlgos/interface/TFileService.h"
#include "DataFormats/Common/interface/Handle.h"
#include "DataFormats/Common/interface/TriggerResults.h"
#include "DataFormats/PatCandidates/interface/Jet.h"
#include "DataFormats/PatCandidates/interface/Muon.h"
#include "DataFormats/PatCandidates/interface/Electron.h"
#include "DataFormats/VertexReco/interface/Vertex.h"
#include "FWCore/Common/interface/TriggerNames.h"
#include "FWCore/Framework/interface/MakerMacros.h"
#include "FWCore/Framework/interface/Event.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "FWCore/ParameterSet/interface/ParameterSet.h"
#include "FWCore/ServiceRegistry/interface/Service.h"
#include "FWCore/Utilities/interface/EDMException.h"
#include "FWCore/Utilities/interface/InputTag.h"
#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"

#include "Tree/System8/interface/S8TreeInfo.h"
#include "Tree/System8/interface/S8EventID.h"
#include "Tree/System8/interface/S8GenEvent.h"
#include "Tree/System8/interface/S8GenParticle.h"
#include "Tree/System8/interface/S8Jet.h"
#include "Tree/System8/interface/S8Lepton.h"
#include "Tree/System8/interface/S8PrimaryVertex.h"
#include "Tree/System8/interface/S8TreeInfo.h"
#include "Tree/System8/interface/S8Trigger.h"
#include "Tree/System8/interface/S8TriggerCenter.h"
#include "EDModule/Analyzer/interface/Tools.h"

#include "EDModule/Analyzer/interface/S8TreeMaker.h"

using std::cout;
using std::endl;
using std::string;
using std::vector;

using boost::lexical_cast;
using boost::regex;
using boost::smatch;

using edm::errors::ErrorCodes;
using edm::Handle;
using edm::InputTag;
using edm::LogWarning;
using edm::LogInfo;
using edm::ParameterSet;
using edm::TriggerResults;
using edm::TriggerNames;

using reco::Vertex;

using s8::TreeInfo;
using s8::TriggerCenter;

using namespace top::tools;

void fillSplit(const int &status, boost::tribool &split)
{
    switch(status)
    {
        case 3: split = false;
                break;

        case 2: if (boost::indeterminate(split))
                    split = true;

                break;
    }
}

S8TreeMaker::S8TreeMaker(const edm::ParameterSet &config):
    _didInitializeHltConfigProvider(false)
{
    _primaryVertices = config.getParameter<string>("primaryVertices");
    _jets = config.getParameter<string>("jets");
    _muons = config.getParameter<string>("muons");
    _electrons = config.getParameter<string>("electrons");
    _triggers = config.getParameter<string>("triggers");

    _jetSelector = config.getParameter<ParameterSet>("jetSelector");

    _isPythia = config.getParameter<bool>("isPythia");
    _saveTriggers = config.getParameter<bool>("saveTriggers");
}

S8TreeMaker::~S8TreeMaker()
{
}

void S8TreeMaker::beginJob()
{
    cout << "Start BeginJob" << endl;

    edm::Service<TFileService> fileService;

    _treeInfo.reset(new TreeInfo);
    _triggerCenter.reset(new TriggerCenter);

    _tree = fileService->make<TTree>("s8", "System8 tree.");

    // Prepare Branches
    //
    _eventID.reset(new s8::EventID());
    _tree->Branch("eventID", _eventID.get());

    _genEvent.reset(new s8::GenEvent());
    _tree->Branch("genEvent", _genEvent.get());

    _s8Electrons.reset(new s8::Leptons());
    _tree->Branch("electrons.", _s8Electrons.get());

    _s8Jets.reset(new s8::Jets());
    _tree->Branch("jets", _s8Jets.get());

    _s8Muons.reset(new s8::Leptons());
    _tree->Branch("muons.", _s8Muons.get());

    _s8PrimaryVertices.reset(new s8::PrimaryVertices());
    _tree->Branch("primaryVertices", _s8PrimaryVertices.get());

    _s8Triggers.reset(new s8::Triggers());
    _tree->Branch("triggers", _s8Triggers.get());

    _didInitializeHltConfigProvider = false;

    cout << "BeginJob is over" << endl;
}

void S8TreeMaker::endJob()
{
    if (!_tree)
        return;

    // Tree Info is disabled for the moment until Hadd is fixed
    //
    edm::Service<TFileService> fileService;
//    TDirectory *dir = fileService->cd();
    TDirectory *dir = fileService->getBareDirectory();
    dir->WriteObject(_treeInfo.get(), "s8info");
    dir->WriteObject(_triggerCenter.get(), "s8triggers");
}

void S8TreeMaker::beginRun(const edm::Run &run,
                           const edm::EventSetup &eventSetup)
{
    using s8::Trigger;

    // Initialize HLT Config Provider for new Run
    //
    bool didChange = true;
    if (!_hltConfigProvider.init(run,
                                 eventSetup,
                                 InputTag(_triggers).process(),
                                 didChange))

        throw cms::Exception("S8TreeMaker")
            << "Failed to initialize HLTConfig for: " << _triggers;

    // Test if Trigger Menu has changed
    //
    if (!didChange &&
        _didInitializeHltConfigProvider)

        return;

    _hlts.clear();

    // Triggers are found and changed
    //
    typedef std::vector<std::string> Triggers;

    const Triggers &triggerNames = _hltConfigProvider.triggerNames();

    using s8::tools::make_hash;

    TriggerCenter::TriggerMap &triggers = _triggerCenter->triggers();

    for(Triggers::const_iterator trigger = triggerNames.begin();
        triggerNames.end() != trigger;
        ++trigger)
    {
        smatch matches;
        if (!regex_match(*trigger, matches,
                         regex("^(\\w+?)(?:_[vV](\\d+))?$")))
        {
            cout << "Do not understand Trigger Name: " << *trigger
                << endl;

            continue;
        }

        string trigger_name = matches[1];
        boost::to_lower(trigger_name);
        triggers.insert(make_pair(make_hash(trigger_name), matches[1]));

        // Found Trigger of the interest
        //
        HLT foundHLT;

        foundHLT.hash = make_hash(trigger_name);
        foundHLT.id = distance(triggerNames.begin(), trigger);
        foundHLT.version = matches[2].matched
            ? lexical_cast<int>(matches[2])
            : 0;

        _hlts[*trigger] = foundHLT;
    }

    if (_hlts.empty())
        LogWarning("S8TreeMaker")
            << "None of the searched HLT Triggers is found" << endl;

    _didInitializeHltConfigProvider = true;
}

void S8TreeMaker::analyze(const edm::Event &event,
                          const edm::EventSetup &eventSetup)
{
    if (!_tree)
        throw cms::Exception("S8TreeMaker")
            << "Tree does not exist";

    processEventID(event);
    processGenEvent(event);
    processElectrons(event);
    processJets(event);
    processMuons(event);
    processPrimaryVertices(event);

    if (_saveTriggers)
        processTriggers(event, eventSetup);

    // Write Tree entry
    //
    _tree->Fill();

    // Reset event
    //
    _eventID->reset();
    _genEvent->reset();

    for(s8::Jets::iterator jet = _s8Jets->begin();
        _s8Jets->end() != jet;
        ++jet)
    {
        delete *jet;
    }
    _s8Jets->clear();

    for(s8::Leptons::iterator electron = _s8Electrons->begin();
        _s8Electrons->end() != electron;
        ++electron)
    {
        delete *electron;
    }
    _s8Electrons->clear();

    for(s8::Leptons::iterator muon = _s8Muons->begin();
        _s8Muons->end() != muon;
        ++muon)
    {
        delete *muon;
    }
    _s8Muons->clear();

    for(s8::PrimaryVertices::iterator primaryVertex = _s8PrimaryVertices->begin();
        _s8PrimaryVertices->end() != primaryVertex;
        ++primaryVertex)
    {
        delete *primaryVertex;
    }
    _s8PrimaryVertices->clear();

    for(s8::Triggers::iterator trigger = _s8Triggers->begin();
        _s8Triggers->end() != trigger;
        ++trigger)
    {
        delete *trigger;
    }
    _s8Triggers->clear();
}

void S8TreeMaker::processEventID(const edm::Event &event)
{
    _eventID->setRun(event.id().run());
    _eventID->setLumiBlock(event.id().luminosityBlock());
    _eventID->setEvent(event.id().event());
}

void S8TreeMaker::processGenEvent(const edm::Event &event)
{
    using reco::GenParticleCollection;

    if (event.isRealData())
        return;

    // check if Event is Pythia
    //
    if (_isPythia)
    {
        Handle<GenEventInfoProduct> generator;
        event.getByLabel(InputTag("generator"), generator);

        if (!generator.isValid())
        {
            LogWarning("S8TreeMaker")
                << "failed to extract Generator";

            return;
        }

        _genEvent->setPtHat(generator->qScale());
    }

    Handle<GenParticleCollection> genParticles;
    event.getByLabel(InputTag("genParticles"), genParticles);

    boost::tribool bsplit = boost::indeterminate;
    boost::tribool csplit = boost::indeterminate;
    for(GenParticleCollection::const_iterator particle = genParticles->begin();
        genParticles->end() != particle;
        ++particle)
    {
        const int pdgID = abs(particle->pdgId());

        switch(pdgID)
        {
            case 5: fillSplit(particle->status(), bsplit);
                    break;

            case 4: fillSplit(particle->status(), csplit);
                    break;

            default: break;
        }
    }

    if (bsplit)
        _genEvent->setGluonSplitting(s8::GenEvent::BB, true);

    if (csplit)
        _genEvent->setGluonSplitting(s8::GenEvent::CC, true);

    // add beam crossing information
    Handle<std::vector< PileupSummaryInfo > >  PupInfo;
    if ( event.getByLabel("addPileupInfo", PupInfo) )
    {
       std::vector<PileupSummaryInfo>::const_iterator PVInfo;

       for(PVInfo = PupInfo->begin(); PVInfo != PupInfo->end(); ++PVInfo) {

         int BX = PVInfo->getBunchCrossing();

         if(BX == -1) {
           _genEvent->setBx_minus1( PVInfo->getPU_NumInteractions() );
         }
         if(BX == 0) {
           _genEvent->setBx_zero( PVInfo->getPU_NumInteractions() );
         }
         if(BX == 1) {
           _genEvent->setBx_plus1( PVInfo->getPU_NumInteractions() );
         }

       }
    }


}

void S8TreeMaker::processElectrons(const edm::Event &event)
{
    using pat::ElectronCollection;

    // Extract Electrons
    //
    Handle<ElectronCollection> electrons;
    event.getByLabel(InputTag(_electrons), electrons);

    if (!electrons.isValid())
    {
        LogWarning("S8TreeMaker")
            << "failed to extract Electrons.";

        return;
    }

    // Process all Electrons
    //
    for(ElectronCollection::const_iterator electron = electrons->begin();
        electrons->end() != electron;
        ++electron)
    {
        s8::Lepton *s8Electron = new s8::Lepton();;

        setP4(s8Electron->p4(), electron->p4());
        setVertex(s8Electron->vertex(), electron->vertex());

        s8Electron->impactParameter()->first = electron->dB();
        s8Electron->impactParameter()->second = electron->edB();

        // Extract GenParticle information
        //
        if (electron->genLepton())
        {
            s8::GenParticle *s8GenParticle = s8Electron->genParticle();

            setP4(s8GenParticle->p4(), electron->genLepton()->p4());
            setVertex(s8GenParticle->vertex(), electron->genLepton()->vertex());

            s8GenParticle->setId(electron->genLepton()->pdgId());
            s8GenParticle->setStatus(electron->genLepton()->status());
            if (electron->genLepton()->mother())
                s8GenParticle->setParentId(electron->genLepton()->mother()->pdgId());
        }

        _s8Electrons->push_back(s8Electron);
    } // End loop over electrons
}

void S8TreeMaker::processJets(const edm::Event &event)
{
    using pat::JetCollection;

    // Extract Jets
    //
    Handle<JetCollection> jets;
    event.getByLabel(InputTag(_jets), jets);

    if (!jets.isValid())
    {
        LogWarning("S8TreeMaker")
            << "failed to extract Jets.";

        return;
    }

    // Process All Jets
    //
    //pat::strbitset jetBitset = _jetSelector.getBitTemplate();
    for(JetCollection::const_iterator jet = jets->begin();
        jets->end() != jet;
        ++jet)
    {
      //if (!_jetSelector(*jet, jetBitset))
      //    continue;

        using s8::Jet;
        Jet *s8Jet = new Jet();

        setP4(s8Jet->p4(), jet->p4());

        s8Jet->setFlavour(jet->partonFlavour());
        s8Jet->setTracks(jet->associatedTracks().size());

        if (jet->genParton())
        {
            s8::GenParticle *s8GenParticle = s8Jet->genParticle();

            setP4(s8GenParticle->p4(), jet->genParton()->p4());
            setVertex(s8GenParticle->vertex(), jet->genParton()->vertex());

            s8GenParticle->setId(jet->genParton()->pdgId());
            s8GenParticle->setStatus(jet->genParton()->status());
            if (jet->genParton()->mother())
                s8GenParticle->setParentId(jet->genParton()->mother()->pdgId());
        }

        // Save b-taggers
        //
        s8Jet->setBTag(Jet::TCHE,
                       jet->bDiscriminator("trackCountingHighEffBJetTags"));

        s8Jet->setBTag(Jet::TCHP,
                       jet->bDiscriminator("trackCountingHighPurBJetTags"));

        s8Jet->setBTag(Jet::JP,
                       jet->bDiscriminator("jetProbabilityBJetTags"));

        s8Jet->setBTag(Jet::SSV,
                       jet->bDiscriminator("simpleSecondaryVertexBJetTags"));

        s8Jet->setBTag(Jet::SSVHE,
                       jet->bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));

        s8Jet->setBTag(Jet::SSVHP,
                       jet->bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));

        s8Jet->setBTag(Jet::CSV,
                       jet->bDiscriminator("combinedSecondaryVertexBJetTags"));

        s8Jet->setBTag(Jet::CSV_MVA,
                       jet->bDiscriminator("combinedSecondaryVertexMVABJetTags"));

        s8Jet->setBTag(Jet::JBP,
                       jet->bDiscriminator("jetBProbabilityBJetTags"));

        _s8Jets->push_back(s8Jet);
    }
}

void S8TreeMaker::processMuons(const edm::Event &event)
{
    using pat::MuonCollection;

    // Extract Muons
    //
    Handle<MuonCollection> muons;
    event.getByLabel(InputTag(_muons), muons);

    if (!muons.isValid())
    {
        LogWarning("S8TreeMaker")
            << "failed to extract Muons.";

        return;
    }

    // Process all Muons
    //
    for(MuonCollection::const_iterator muon = muons->begin();
        muons->end() != muon;
        ++muon)
    {
        // Only Muons with basic cuts are saved:
        //
        if (1 >= muon->numberOfMatches())
            continue;

        s8::Lepton *s8Muon = new s8::Lepton();

        setP4(s8Muon->p4(), muon->p4());
        setVertex(s8Muon->vertex(), muon->vertex());

        s8Muon->impactParameter()->first = muon->dB();
        s8Muon->impactParameter()->second = muon->edB();

        // Extract GenParticle information
        //
        if (muon->genLepton())
        {
            s8::GenParticle *s8GenParticle = s8Muon->genParticle();

            setP4(s8GenParticle->p4(), muon->genLepton()->p4());
            setVertex(s8GenParticle->vertex(), muon->genLepton()->vertex());

            s8GenParticle->setId(muon->genLepton()->pdgId());
            s8GenParticle->setStatus(muon->genLepton()->status());
            if (muon->genLepton()->mother())
                s8GenParticle->setParentId(muon->genLepton()->mother()->pdgId());
        }

        _s8Muons->push_back(s8Muon);
    } // End loop over muons
}

void S8TreeMaker::processPrimaryVertices(const edm::Event &event)
{
    // Primary Vertices
    //
    typedef vector<Vertex> PVCollection;

    Handle<PVCollection> primaryVertices;
    event.getByLabel(InputTag(_primaryVertices), primaryVertices);

    if (!primaryVertices.isValid())
    {
        LogWarning("S8TreeMaker")
            << "failed to extract Primary Vertices.";

        return;
    }

    if (primaryVertices->empty())
    {
        LogWarning("S8TreeMaker")
            << "primary vertices collection is empty.";

        return;
    }
    // Process Primary Vertices
    //
    for(PVCollection::const_iterator vertex = primaryVertices->begin();
        primaryVertices->end() != vertex;
        ++vertex)
    {
        if (!isGoodPrimaryVertex(*vertex, event.isRealData()))
            continue;

        s8::PrimaryVertex *s8Vertex = new s8::PrimaryVertex();

        setVertex(s8Vertex->vertex(), vertex->position());
        s8Vertex->setNdof(vertex->ndof());
        s8Vertex->setRho(vertex->position().Rho());

        _s8PrimaryVertices->push_back(s8Vertex);
    }
}

void S8TreeMaker::processTriggers(const edm::Event &event,
                                  const edm::EventSetup &eventSetup)
{
    if (_hlts.empty())
        return;

    // Triggers
    //
    Handle<TriggerResults> triggers;
    event.getByLabel(InputTag(_triggers), triggers);

    if (!triggers.isValid())
    {
        LogWarning("S8TreeMaker")
            << "failed to extract Triggers";

        return;
    }

    // Process only found HLTs
    //
    for(HLTs::const_iterator hlt = _hlts.begin();
        _hlts.end() != hlt;
        ++hlt)
    {
        s8::Trigger *s8Trigger = new s8::Trigger();
        s8Trigger->setHash(hlt->second.hash);
        if (hlt->second.version)
            s8Trigger->setVersion(hlt->second.version);

        s8Trigger->setIsPass(triggers->accept(hlt->second.id));
        s8Trigger->setPrescale(_hltConfigProvider.prescaleValue(event,
            eventSetup, hlt->first));


        _s8Triggers->push_back(s8Trigger);
    }
}

bool S8TreeMaker::isGoodPrimaryVertex(const Vertex &vertex,
                                    const bool &isData)
{
    return !vertex.isFake() &&
            4 <= vertex.ndof() &&
            (isData ? 24 : 15) >= fabs(vertex.z()) &&
            2 >= fabs(vertex.position().Rho());
}

DEFINE_FWK_MODULE(S8TreeMaker);
Revision:	1.5
Committed:	Fri Dec 16 00:35:40 2011 UTC (13 years, 4 months ago) by meloam
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V00-00-04, HEAD
Changes since 1.4:	+0 -2 lines
Log Message:	updating EDModule for dec 2011 run
#	User	Rev	Content
1	samvel	1.1	/**
2			* S8TreeMaker
3			*
4			*
5			* Created by Samvel Khalatian on Sep 29, 2010
6			* Copyright 2010, All rights reserved
7			*/
8
9			#include <iostream>
10			#include <vector>
11
12			#include <boost/algorithm/string.hpp>
13			#include <boost/lexical_cast.hpp>
14			#include <boost/logic/tribool.hpp>
15			#include <boost/regex.hpp>
16
17			#include <TTree.h>
18			#include <TLorentzVector.h>
19			#include <TVector3.h>
20
21			#include "CommonTools/UtilAlgos/interface/TFileService.h"
22			#include "DataFormats/Common/interface/Handle.h"
23			#include "DataFormats/Common/interface/TriggerResults.h"
24			#include "DataFormats/PatCandidates/interface/Jet.h"
25			#include "DataFormats/PatCandidates/interface/Muon.h"
26			#include "DataFormats/PatCandidates/interface/Electron.h"
27			#include "DataFormats/VertexReco/interface/Vertex.h"
28			#include "FWCore/Common/interface/TriggerNames.h"
29			#include "FWCore/Framework/interface/MakerMacros.h"
30			#include "FWCore/Framework/interface/Event.h"
31			#include "FWCore/MessageLogger/interface/MessageLogger.h"
32			#include "FWCore/ParameterSet/interface/ParameterSet.h"
33			#include "FWCore/ServiceRegistry/interface/Service.h"
34			#include "FWCore/Utilities/interface/EDMException.h"
35			#include "FWCore/Utilities/interface/InputTag.h"
36			#include "SimDataFormats/GeneratorProducts/interface/GenEventInfoProduct.h"
37	meloam	1.4	#include "SimDataFormats/PileupSummaryInfo/interface/PileupSummaryInfo.h"
38	samvel	1.1
39			#include "Tree/System8/interface/S8TreeInfo.h"
40			#include "Tree/System8/interface/S8EventID.h"
41			#include "Tree/System8/interface/S8GenEvent.h"
42			#include "Tree/System8/interface/S8GenParticle.h"
43			#include "Tree/System8/interface/S8Jet.h"
44			#include "Tree/System8/interface/S8Lepton.h"
45			#include "Tree/System8/interface/S8PrimaryVertex.h"
46			#include "Tree/System8/interface/S8TreeInfo.h"
47			#include "Tree/System8/interface/S8Trigger.h"
48			#include "Tree/System8/interface/S8TriggerCenter.h"
49			#include "EDModule/Analyzer/interface/Tools.h"
50
51			#include "EDModule/Analyzer/interface/S8TreeMaker.h"
52
53			using std::cout;
54			using std::endl;
55			using std::string;
56			using std::vector;
57
58			using boost::lexical_cast;
59			using boost::regex;
60			using boost::smatch;
61
62			using edm::errors::ErrorCodes;
63			using edm::Handle;
64			using edm::InputTag;
65			using edm::LogWarning;
66			using edm::LogInfo;
67			using edm::ParameterSet;
68			using edm::TriggerResults;
69			using edm::TriggerNames;
70
71			using reco::Vertex;
72
73			using s8::TreeInfo;
74			using s8::TriggerCenter;
75
76			using namespace top::tools;
77
78			void fillSplit(const int &status, boost::tribool &split)
79			{
80			switch(status)
81			{
82			case 3: split = false;
83			break;
84
85			case 2: if (boost::indeterminate(split))
86			split = true;
87
88			break;
89			}
90			}
91
92			S8TreeMaker::S8TreeMaker(const edm::ParameterSet &config):
93			_didInitializeHltConfigProvider(false)
94			{
95			_primaryVertices = config.getParameter<string>("primaryVertices");
96			_jets = config.getParameter<string>("jets");
97			_muons = config.getParameter<string>("muons");
98			_electrons = config.getParameter<string>("electrons");
99			_triggers = config.getParameter<string>("triggers");
100
101			_jetSelector = config.getParameter<ParameterSet>("jetSelector");
102
103			_isPythia = config.getParameter<bool>("isPythia");
104			_saveTriggers = config.getParameter<bool>("saveTriggers");
105			}
106
107			S8TreeMaker::~S8TreeMaker()
108			{
109			}
110
111			void S8TreeMaker::beginJob()
112			{
113			cout << "Start BeginJob" << endl;
114
115			edm::Service<TFileService> fileService;
116
117			_treeInfo.reset(new TreeInfo);
118			_triggerCenter.reset(new TriggerCenter);
119
120			_tree = fileService->make<TTree>("s8", "System8 tree.");
121
122			// Prepare Branches
123			//
124			_eventID.reset(new s8::EventID());
125			_tree->Branch("eventID", _eventID.get());
126
127			_genEvent.reset(new s8::GenEvent());
128			_tree->Branch("genEvent", _genEvent.get());
129
130			_s8Electrons.reset(new s8::Leptons());
131			_tree->Branch("electrons.", _s8Electrons.get());
132
133			_s8Jets.reset(new s8::Jets());
134			_tree->Branch("jets", _s8Jets.get());
135
136			_s8Muons.reset(new s8::Leptons());
137			_tree->Branch("muons.", _s8Muons.get());
138
139			_s8PrimaryVertices.reset(new s8::PrimaryVertices());
140			_tree->Branch("primaryVertices", _s8PrimaryVertices.get());
141
142			_s8Triggers.reset(new s8::Triggers());
143			_tree->Branch("triggers", _s8Triggers.get());
144
145			_didInitializeHltConfigProvider = false;
146
147			cout << "BeginJob is over" << endl;
148			}
149
150			void S8TreeMaker::endJob()
151			{
152			if (!_tree)
153			return;
154
155			// Tree Info is disabled for the moment until Hadd is fixed
156			//
157			edm::Service<TFileService> fileService;
158	yumiceva	1.2	// TDirectory *dir = fileService->cd();
159			TDirectory *dir = fileService->getBareDirectory();
160	samvel	1.1	dir->WriteObject(_treeInfo.get(), "s8info");
161			dir->WriteObject(_triggerCenter.get(), "s8triggers");
162			}
163
164			void S8TreeMaker::beginRun(const edm::Run &run,
165			const edm::EventSetup &eventSetup)
166			{
167			using s8::Trigger;
168
169			// Initialize HLT Config Provider for new Run
170			//
171			bool didChange = true;
172			if (!_hltConfigProvider.init(run,
173			eventSetup,
174			InputTag(_triggers).process(),
175			didChange))
176
177			throw cms::Exception("S8TreeMaker")
178			<< "Failed to initialize HLTConfig for: " << _triggers;
179
180			// Test if Trigger Menu has changed
181			//
182			if (!didChange &&
183			_didInitializeHltConfigProvider)
184
185			return;
186
187			_hlts.clear();
188
189			// Triggers are found and changed
190			//
191			typedef std::vector<std::string> Triggers;
192
193			const Triggers &triggerNames = _hltConfigProvider.triggerNames();
194
195			using s8::tools::make_hash;
196
197			TriggerCenter::TriggerMap &triggers = _triggerCenter->triggers();
198
199			for(Triggers::const_iterator trigger = triggerNames.begin();
200			triggerNames.end() != trigger;
201			++trigger)
202			{
203			smatch matches;
204			if (!regex_match(*trigger, matches,
205			regex("^(\\w+?)(?:_[vV](\\d+))?$")))
206			{
207			cout << "Do not understand Trigger Name: " << *trigger
208			<< endl;
209
210			continue;
211			}
212
213			string trigger_name = matches[1];
214			boost::to_lower(trigger_name);
215			triggers.insert(make_pair(make_hash(trigger_name), matches[1]));
216
217			// Found Trigger of the interest
218			//
219			HLT foundHLT;
220
221			foundHLT.hash = make_hash(trigger_name);
222			foundHLT.id = distance(triggerNames.begin(), trigger);
223			foundHLT.version = matches[2].matched
224			? lexical_cast<int>(matches[2])
225			: 0;
226
227			_hlts[*trigger] = foundHLT;
228			}
229
230			if (_hlts.empty())
231			LogWarning("S8TreeMaker")
232			<< "None of the searched HLT Triggers is found" << endl;
233
234			_didInitializeHltConfigProvider = true;
235			}
236
237			void S8TreeMaker::analyze(const edm::Event &event,
238			const edm::EventSetup &eventSetup)
239			{
240			if (!_tree)
241			throw cms::Exception("S8TreeMaker")
242			<< "Tree does not exist";
243
244			processEventID(event);
245			processGenEvent(event);
246			processElectrons(event);
247			processJets(event);
248			processMuons(event);
249			processPrimaryVertices(event);
250
251			if (_saveTriggers)
252			processTriggers(event, eventSetup);
253
254			// Write Tree entry
255			//
256			_tree->Fill();
257
258			// Reset event
259			//
260			_eventID->reset();
261			_genEvent->reset();
262
263			for(s8::Jets::iterator jet = _s8Jets->begin();
264			_s8Jets->end() != jet;
265			++jet)
266			{
267			delete *jet;
268			}
269			_s8Jets->clear();
270
271			for(s8::Leptons::iterator electron = _s8Electrons->begin();
272			_s8Electrons->end() != electron;
273			++electron)
274			{
275			delete *electron;
276			}
277			_s8Electrons->clear();
278
279			for(s8::Leptons::iterator muon = _s8Muons->begin();
280			_s8Muons->end() != muon;
281			++muon)
282			{
283			delete *muon;
284			}
285			_s8Muons->clear();
286
287			for(s8::PrimaryVertices::iterator primaryVertex = _s8PrimaryVertices->begin();
288			_s8PrimaryVertices->end() != primaryVertex;
289			++primaryVertex)
290			{
291			delete *primaryVertex;
292			}
293			_s8PrimaryVertices->clear();
294
295			for(s8::Triggers::iterator trigger = _s8Triggers->begin();
296			_s8Triggers->end() != trigger;
297			++trigger)
298			{
299			delete *trigger;
300			}
301			_s8Triggers->clear();
302			}
303
304			void S8TreeMaker::processEventID(const edm::Event &event)
305			{
306			_eventID->setRun(event.id().run());
307			_eventID->setLumiBlock(event.id().luminosityBlock());
308			_eventID->setEvent(event.id().event());
309			}
310
311			void S8TreeMaker::processGenEvent(const edm::Event &event)
312			{
313			using reco::GenParticleCollection;
314
315			if (event.isRealData())
316			return;
317
318			// check if Event is Pythia
319			//
320			if (_isPythia)
321			{
322			Handle<GenEventInfoProduct> generator;
323			event.getByLabel(InputTag("generator"), generator);
324
325			if (!generator.isValid())
326			{
327			LogWarning("S8TreeMaker")
328			<< "failed to extract Generator";
329
330			return;
331			}
332
333			_genEvent->setPtHat(generator->qScale());
334			}
335
336			Handle<GenParticleCollection> genParticles;
337			event.getByLabel(InputTag("genParticles"), genParticles);
338
339			boost::tribool bsplit = boost::indeterminate;
340			boost::tribool csplit = boost::indeterminate;
341			for(GenParticleCollection::const_iterator particle = genParticles->begin();
342			genParticles->end() != particle;
343			++particle)
344			{
345			const int pdgID = abs(particle->pdgId());
346
347			switch(pdgID)
348			{
349			case 5: fillSplit(particle->status(), bsplit);
350			break;
351
352			case 4: fillSplit(particle->status(), csplit);
353			break;
354
355			default: break;
356			}
357			}
358
359			if (bsplit)
360			_genEvent->setGluonSplitting(s8::GenEvent::BB, true);
361
362			if (csplit)
363			_genEvent->setGluonSplitting(s8::GenEvent::CC, true);
364	meloam	1.4
365			// add beam crossing information
366			Handle<std::vector< PileupSummaryInfo > > PupInfo;
367			if ( event.getByLabel("addPileupInfo", PupInfo) )
368			{
369			std::vector<PileupSummaryInfo>::const_iterator PVInfo;
370
371			for(PVInfo = PupInfo->begin(); PVInfo != PupInfo->end(); ++PVInfo) {
372
373			int BX = PVInfo->getBunchCrossing();
374
375			if(BX == -1) {
376			_genEvent->setBx_minus1( PVInfo->getPU_NumInteractions() );
377			}
378			if(BX == 0) {
379			_genEvent->setBx_zero( PVInfo->getPU_NumInteractions() );
380			}
381			if(BX == 1) {
382			_genEvent->setBx_plus1( PVInfo->getPU_NumInteractions() );
383			}
384
385			}
386			}
387
388
389	samvel	1.1	}
390
391			void S8TreeMaker::processElectrons(const edm::Event &event)
392			{
393			using pat::ElectronCollection;
394
395			// Extract Electrons
396			//
397			Handle<ElectronCollection> electrons;
398			event.getByLabel(InputTag(_electrons), electrons);
399
400			if (!electrons.isValid())
401			{
402			LogWarning("S8TreeMaker")
403			<< "failed to extract Electrons.";
404
405			return;
406			}
407
408			// Process all Electrons
409			//
410			for(ElectronCollection::const_iterator electron = electrons->begin();
411			electrons->end() != electron;
412			++electron)
413			{
414			s8::Lepton *s8Electron = new s8::Lepton();;
415
416			setP4(s8Electron->p4(), electron->p4());
417			setVertex(s8Electron->vertex(), electron->vertex());
418
419			s8Electron->impactParameter()->first = electron->dB();
420			s8Electron->impactParameter()->second = electron->edB();
421
422			// Extract GenParticle information
423			//
424			if (electron->genLepton())
425			{
426			s8::GenParticle *s8GenParticle = s8Electron->genParticle();
427
428			setP4(s8GenParticle->p4(), electron->genLepton()->p4());
429			setVertex(s8GenParticle->vertex(), electron->genLepton()->vertex());
430
431			s8GenParticle->setId(electron->genLepton()->pdgId());
432			s8GenParticle->setStatus(electron->genLepton()->status());
433			if (electron->genLepton()->mother())
434			s8GenParticle->setParentId(electron->genLepton()->mother()->pdgId());
435			}
436
437			_s8Electrons->push_back(s8Electron);
438			} // End loop over electrons
439			}
440
441			void S8TreeMaker::processJets(const edm::Event &event)
442			{
443			using pat::JetCollection;
444
445			// Extract Jets
446			//
447			Handle<JetCollection> jets;
448			event.getByLabel(InputTag(_jets), jets);
449
450			if (!jets.isValid())
451			{
452			LogWarning("S8TreeMaker")
453			<< "failed to extract Jets.";
454
455			return;
456			}
457
458			// Process All Jets
459			//
460	yumiceva	1.3	//pat::strbitset jetBitset = _jetSelector.getBitTemplate();
461	samvel	1.1	for(JetCollection::const_iterator jet = jets->begin();
462			jets->end() != jet;
463			++jet)
464			{
465	yumiceva	1.3	//if (!_jetSelector(*jet, jetBitset))
466			// continue;
467	samvel	1.1
468			using s8::Jet;
469			Jet *s8Jet = new Jet();
470
471			setP4(s8Jet->p4(), jet->p4());
472
473			s8Jet->setFlavour(jet->partonFlavour());
474			s8Jet->setTracks(jet->associatedTracks().size());
475
476			if (jet->genParton())
477			{
478			s8::GenParticle *s8GenParticle = s8Jet->genParticle();
479
480			setP4(s8GenParticle->p4(), jet->genParton()->p4());
481			setVertex(s8GenParticle->vertex(), jet->genParton()->vertex());
482
483			s8GenParticle->setId(jet->genParton()->pdgId());
484			s8GenParticle->setStatus(jet->genParton()->status());
485			if (jet->genParton()->mother())
486			s8GenParticle->setParentId(jet->genParton()->mother()->pdgId());
487			}
488
489			// Save b-taggers
490			//
491			s8Jet->setBTag(Jet::TCHE,
492			jet->bDiscriminator("trackCountingHighEffBJetTags"));
493
494			s8Jet->setBTag(Jet::TCHP,
495			jet->bDiscriminator("trackCountingHighPurBJetTags"));
496
497			s8Jet->setBTag(Jet::JP,
498			jet->bDiscriminator("jetProbabilityBJetTags"));
499
500			s8Jet->setBTag(Jet::SSV,
501			jet->bDiscriminator("simpleSecondaryVertexBJetTags"));
502
503			s8Jet->setBTag(Jet::SSVHE,
504			jet->bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
505
506			s8Jet->setBTag(Jet::SSVHP,
507			jet->bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
508
509			s8Jet->setBTag(Jet::CSV,
510			jet->bDiscriminator("combinedSecondaryVertexBJetTags"));
511
512			s8Jet->setBTag(Jet::CSV_MVA,
513			jet->bDiscriminator("combinedSecondaryVertexMVABJetTags"));
514
515			s8Jet->setBTag(Jet::JBP,
516			jet->bDiscriminator("jetBProbabilityBJetTags"));
517
518			_s8Jets->push_back(s8Jet);
519			}
520			}
521
522			void S8TreeMaker::processMuons(const edm::Event &event)
523			{
524			using pat::MuonCollection;
525
526			// Extract Muons
527			//
528			Handle<MuonCollection> muons;
529			event.getByLabel(InputTag(_muons), muons);
530
531			if (!muons.isValid())
532			{
533			LogWarning("S8TreeMaker")
534			<< "failed to extract Muons.";
535
536			return;
537			}
538
539			// Process all Muons
540			//
541			for(MuonCollection::const_iterator muon = muons->begin();
542			muons->end() != muon;
543			++muon)
544			{
545			// Only Muons with basic cuts are saved:
546			//
547			if (1 >= muon->numberOfMatches())
548			continue;
549
550			s8::Lepton *s8Muon = new s8::Lepton();
551
552			setP4(s8Muon->p4(), muon->p4());
553			setVertex(s8Muon->vertex(), muon->vertex());
554
555			s8Muon->impactParameter()->first = muon->dB();
556			s8Muon->impactParameter()->second = muon->edB();
557
558			// Extract GenParticle information
559			//
560			if (muon->genLepton())
561			{
562			s8::GenParticle *s8GenParticle = s8Muon->genParticle();
563
564			setP4(s8GenParticle->p4(), muon->genLepton()->p4());
565			setVertex(s8GenParticle->vertex(), muon->genLepton()->vertex());
566
567			s8GenParticle->setId(muon->genLepton()->pdgId());
568			s8GenParticle->setStatus(muon->genLepton()->status());
569			if (muon->genLepton()->mother())
570			s8GenParticle->setParentId(muon->genLepton()->mother()->pdgId());
571			}
572
573			_s8Muons->push_back(s8Muon);
574			} // End loop over muons
575			}
576
577			void S8TreeMaker::processPrimaryVertices(const edm::Event &event)
578			{
579			// Primary Vertices
580			//
581			typedef vector<Vertex> PVCollection;
582
583			Handle<PVCollection> primaryVertices;
584			event.getByLabel(InputTag(_primaryVertices), primaryVertices);
585
586			if (!primaryVertices.isValid())
587			{
588			LogWarning("S8TreeMaker")
589			<< "failed to extract Primary Vertices.";
590
591			return;
592			}
593
594			if (primaryVertices->empty())
595			{
596			LogWarning("S8TreeMaker")
597			<< "primary vertices collection is empty.";
598
599			return;
600			}
601			// Process Primary Vertices
602			//
603			for(PVCollection::const_iterator vertex = primaryVertices->begin();
604			primaryVertices->end() != vertex;
605			++vertex)
606			{
607			if (!isGoodPrimaryVertex(*vertex, event.isRealData()))
608			continue;
609
610			s8::PrimaryVertex *s8Vertex = new s8::PrimaryVertex();
611
612			setVertex(s8Vertex->vertex(), vertex->position());
613			s8Vertex->setNdof(vertex->ndof());
614			s8Vertex->setRho(vertex->position().Rho());
615
616			_s8PrimaryVertices->push_back(s8Vertex);
617			}
618			}
619
620			void S8TreeMaker::processTriggers(const edm::Event &event,
621			const edm::EventSetup &eventSetup)
622			{
623			if (_hlts.empty())
624			return;
625
626			// Triggers
627			//
628			Handle<TriggerResults> triggers;
629			event.getByLabel(InputTag(_triggers), triggers);
630
631			if (!triggers.isValid())
632			{
633			LogWarning("S8TreeMaker")
634			<< "failed to extract Triggers";
635
636			return;
637			}
638
639			// Process only found HLTs
640			//
641			for(HLTs::const_iterator hlt = _hlts.begin();
642			_hlts.end() != hlt;
643			++hlt)
644			{
645			s8::Trigger *s8Trigger = new s8::Trigger();
646			s8Trigger->setHash(hlt->second.hash);
647			if (hlt->second.version)
648			s8Trigger->setVersion(hlt->second.version);
649
650			s8Trigger->setIsPass(triggers->accept(hlt->second.id));
651			s8Trigger->setPrescale(_hltConfigProvider.prescaleValue(event,
652			eventSetup, hlt->first));
653
654
655			_s8Triggers->push_back(s8Trigger);
656			}
657			}
658
659			bool S8TreeMaker::isGoodPrimaryVertex(const Vertex &vertex,
660			const bool &isData)
661			{
662			return !vertex.isFake() &&
663			4 <= vertex.ndof() &&
664			(isData ? 24 : 15) >= fabs(vertex.z()) &&
665			2 >= fabs(vertex.position().Rho());
666			}
667
668			DEFINE_FWK_MODULE(S8TreeMaker);