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) {
            LogWarning("S8TreeMaker") << "Got bx_0";
           _genEvent->setBx_zero( PVInfo->getPU_NumInteractions() );
         }
         if(BX == 1) {
            LogWarning("S8TreeMaker") << "Got 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.4
Committed:	Tue Nov 22 23:24:21 2011 UTC (13 years, 5 months ago) by meloam
Content type:	text/plain
Branch:	MAIN
CVS Tags:	V00-00-03
Changes since 1.3:	+28 -1 lines
Log Message:	adding mctrue bx_plus1 bx_zero bx_minus1
#	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			LogWarning("S8TreeMaker") << "Got bx_0";
380			_genEvent->setBx_zero( PVInfo->getPU_NumInteractions() );
381			}
382			if(BX == 1) {
383			LogWarning("S8TreeMaker") << "Got bx_1";
384			_genEvent->setBx_plus1( PVInfo->getPU_NumInteractions() );
385			}
386
387			}
388			}
389
390
391	samvel	1.1	}
392
393			void S8TreeMaker::processElectrons(const edm::Event &event)
394			{
395			using pat::ElectronCollection;
396
397			// Extract Electrons
398			//
399			Handle<ElectronCollection> electrons;
400			event.getByLabel(InputTag(_electrons), electrons);
401
402			if (!electrons.isValid())
403			{
404			LogWarning("S8TreeMaker")
405			<< "failed to extract Electrons.";
406
407			return;
408			}
409
410			// Process all Electrons
411			//
412			for(ElectronCollection::const_iterator electron = electrons->begin();
413			electrons->end() != electron;
414			++electron)
415			{
416			s8::Lepton *s8Electron = new s8::Lepton();;
417
418			setP4(s8Electron->p4(), electron->p4());
419			setVertex(s8Electron->vertex(), electron->vertex());
420
421			s8Electron->impactParameter()->first = electron->dB();
422			s8Electron->impactParameter()->second = electron->edB();
423
424			// Extract GenParticle information
425			//
426			if (electron->genLepton())
427			{
428			s8::GenParticle *s8GenParticle = s8Electron->genParticle();
429
430			setP4(s8GenParticle->p4(), electron->genLepton()->p4());
431			setVertex(s8GenParticle->vertex(), electron->genLepton()->vertex());
432
433			s8GenParticle->setId(electron->genLepton()->pdgId());
434			s8GenParticle->setStatus(electron->genLepton()->status());
435			if (electron->genLepton()->mother())
436			s8GenParticle->setParentId(electron->genLepton()->mother()->pdgId());
437			}
438
439			_s8Electrons->push_back(s8Electron);
440			} // End loop over electrons
441			}
442
443			void S8TreeMaker::processJets(const edm::Event &event)
444			{
445			using pat::JetCollection;
446
447			// Extract Jets
448			//
449			Handle<JetCollection> jets;
450			event.getByLabel(InputTag(_jets), jets);
451
452			if (!jets.isValid())
453			{
454			LogWarning("S8TreeMaker")
455			<< "failed to extract Jets.";
456
457			return;
458			}
459
460			// Process All Jets
461			//
462	yumiceva	1.3	//pat::strbitset jetBitset = _jetSelector.getBitTemplate();
463	samvel	1.1	for(JetCollection::const_iterator jet = jets->begin();
464			jets->end() != jet;
465			++jet)
466			{
467	yumiceva	1.3	//if (!_jetSelector(*jet, jetBitset))
468			// continue;
469	samvel	1.1
470			using s8::Jet;
471			Jet *s8Jet = new Jet();
472
473			setP4(s8Jet->p4(), jet->p4());
474
475			s8Jet->setFlavour(jet->partonFlavour());
476			s8Jet->setTracks(jet->associatedTracks().size());
477
478			if (jet->genParton())
479			{
480			s8::GenParticle *s8GenParticle = s8Jet->genParticle();
481
482			setP4(s8GenParticle->p4(), jet->genParton()->p4());
483			setVertex(s8GenParticle->vertex(), jet->genParton()->vertex());
484
485			s8GenParticle->setId(jet->genParton()->pdgId());
486			s8GenParticle->setStatus(jet->genParton()->status());
487			if (jet->genParton()->mother())
488			s8GenParticle->setParentId(jet->genParton()->mother()->pdgId());
489			}
490
491			// Save b-taggers
492			//
493			s8Jet->setBTag(Jet::TCHE,
494			jet->bDiscriminator("trackCountingHighEffBJetTags"));
495
496			s8Jet->setBTag(Jet::TCHP,
497			jet->bDiscriminator("trackCountingHighPurBJetTags"));
498
499			s8Jet->setBTag(Jet::JP,
500			jet->bDiscriminator("jetProbabilityBJetTags"));
501
502			s8Jet->setBTag(Jet::SSV,
503			jet->bDiscriminator("simpleSecondaryVertexBJetTags"));
504
505			s8Jet->setBTag(Jet::SSVHE,
506			jet->bDiscriminator("simpleSecondaryVertexHighEffBJetTags"));
507
508			s8Jet->setBTag(Jet::SSVHP,
509			jet->bDiscriminator("simpleSecondaryVertexHighPurBJetTags"));
510
511			s8Jet->setBTag(Jet::CSV,
512			jet->bDiscriminator("combinedSecondaryVertexBJetTags"));
513
514			s8Jet->setBTag(Jet::CSV_MVA,
515			jet->bDiscriminator("combinedSecondaryVertexMVABJetTags"));
516
517			s8Jet->setBTag(Jet::JBP,
518			jet->bDiscriminator("jetBProbabilityBJetTags"));
519
520			_s8Jets->push_back(s8Jet);
521			}
522			}
523
524			void S8TreeMaker::processMuons(const edm::Event &event)
525			{
526			using pat::MuonCollection;
527
528			// Extract Muons
529			//
530			Handle<MuonCollection> muons;
531			event.getByLabel(InputTag(_muons), muons);
532
533			if (!muons.isValid())
534			{
535			LogWarning("S8TreeMaker")
536			<< "failed to extract Muons.";
537
538			return;
539			}
540
541			// Process all Muons
542			//
543			for(MuonCollection::const_iterator muon = muons->begin();
544			muons->end() != muon;
545			++muon)
546			{
547			// Only Muons with basic cuts are saved:
548			//
549			if (1 >= muon->numberOfMatches())
550			continue;
551
552			s8::Lepton *s8Muon = new s8::Lepton();
553
554			setP4(s8Muon->p4(), muon->p4());
555			setVertex(s8Muon->vertex(), muon->vertex());
556
557			s8Muon->impactParameter()->first = muon->dB();
558			s8Muon->impactParameter()->second = muon->edB();
559
560			// Extract GenParticle information
561			//
562			if (muon->genLepton())
563			{
564			s8::GenParticle *s8GenParticle = s8Muon->genParticle();
565
566			setP4(s8GenParticle->p4(), muon->genLepton()->p4());
567			setVertex(s8GenParticle->vertex(), muon->genLepton()->vertex());
568
569			s8GenParticle->setId(muon->genLepton()->pdgId());
570			s8GenParticle->setStatus(muon->genLepton()->status());
571			if (muon->genLepton()->mother())
572			s8GenParticle->setParentId(muon->genLepton()->mother()->pdgId());
573			}
574
575			_s8Muons->push_back(s8Muon);
576			} // End loop over muons
577			}
578
579			void S8TreeMaker::processPrimaryVertices(const edm::Event &event)
580			{
581			// Primary Vertices
582			//
583			typedef vector<Vertex> PVCollection;
584
585			Handle<PVCollection> primaryVertices;
586			event.getByLabel(InputTag(_primaryVertices), primaryVertices);
587
588			if (!primaryVertices.isValid())
589			{
590			LogWarning("S8TreeMaker")
591			<< "failed to extract Primary Vertices.";
592
593			return;
594			}
595
596			if (primaryVertices->empty())
597			{
598			LogWarning("S8TreeMaker")
599			<< "primary vertices collection is empty.";
600
601			return;
602			}
603			// Process Primary Vertices
604			//
605			for(PVCollection::const_iterator vertex = primaryVertices->begin();
606			primaryVertices->end() != vertex;
607			++vertex)
608			{
609			if (!isGoodPrimaryVertex(*vertex, event.isRealData()))
610			continue;
611
612			s8::PrimaryVertex *s8Vertex = new s8::PrimaryVertex();
613
614			setVertex(s8Vertex->vertex(), vertex->position());
615			s8Vertex->setNdof(vertex->ndof());
616			s8Vertex->setRho(vertex->position().Rho());
617
618			_s8PrimaryVertices->push_back(s8Vertex);
619			}
620			}
621
622			void S8TreeMaker::processTriggers(const edm::Event &event,
623			const edm::EventSetup &eventSetup)
624			{
625			if (_hlts.empty())
626			return;
627
628			// Triggers
629			//
630			Handle<TriggerResults> triggers;
631			event.getByLabel(InputTag(_triggers), triggers);
632
633			if (!triggers.isValid())
634			{
635			LogWarning("S8TreeMaker")
636			<< "failed to extract Triggers";
637
638			return;
639			}
640
641			// Process only found HLTs
642			//
643			for(HLTs::const_iterator hlt = _hlts.begin();
644			_hlts.end() != hlt;
645			++hlt)
646			{
647			s8::Trigger *s8Trigger = new s8::Trigger();
648			s8Trigger->setHash(hlt->second.hash);
649			if (hlt->second.version)
650			s8Trigger->setVersion(hlt->second.version);
651
652			s8Trigger->setIsPass(triggers->accept(hlt->second.id));
653			s8Trigger->setPrescale(_hltConfigProvider.prescaleValue(event,
654			eventSetup, hlt->first));
655
656
657			_s8Triggers->push_back(s8Trigger);
658			}
659			}
660
661			bool S8TreeMaker::isGoodPrimaryVertex(const Vertex &vertex,
662			const bool &isData)
663			{
664			return !vertex.isFake() &&
665			4 <= vertex.ndof() &&
666			(isData ? 24 : 15) >= fabs(vertex.z()) &&
667			2 >= fabs(vertex.position().Rho());
668			}
669
670			DEFINE_FWK_MODULE(S8TreeMaker);