Analyzer/src/MonitorAnalyzer.cpp

/**
 * MonitorAnalyzer
 * s8
 *
 * Created by Samvel Khalatian on Nov 17, 2010
 * Copyright 2010, All rights reserved
 */

#include <cmath>

#include <iostream>
#include <stdexcept>
#include <algorithm>

#include <TDirectory.h>
#include <TFile.h>
#include <TH1F.h>
#include <TH2F.h>
#include <TH3F.h>
#include <TH3D.h>
#include <TH1D.h>
#include <TLorentzVector.h>
#include <TVector3.h>

#include "s8/IO/interface/Event.h"
#include "Tree/System8/interface/S8GenEvent.h"
#include "Tree/System8/interface/S8Fwd.h"
#include "Tree/System8/interface/S8Jet.h"
#include "Tree/System8/interface/S8Lepton.h"
#include "Tree/System8/interface/S8PrimaryVertex.h"
#include "s8/Selector/interface/S8Selector.h"
#include "s8/Utility/interface/MonitorPlots.h"
#include "s8/Utility/interface/MuonInJet.h"
#include "s8/Utility/interface/TaggerOperatingPoint.h"

#include "s8/Analyzer/interface/MonitorAnalyzer.h"

using std::cerr;
using std::clog;
using std::cout;
using std::endl;
using std::runtime_error;

using s8::MonitorAnalyzer;

MonitorAnalyzer::MonitorAnalyzer() throw()
{
    _s8_selector = 0;
    _reweights = 0;
    _event = 0;
}

MonitorAnalyzer::~MonitorAnalyzer() throw()
{
    if (_s8_selector)
        delete _s8_selector;
}


// Analyzer interface
//
void MonitorAnalyzer::init()
{
    _muonInJetTagger.reset(new TaggerOperatingPoint());
    _awayJetTagger.reset(new TaggerOperatingPoint());

    _muonInJet.reset(new MuonInJet());
    _muonInJet->setDelegate(this);
    _muonInJet->setAwayJetTaggerOperatingPoint(_awayJetTagger.get());

    _muons.reset(new TH1D("muons", "N muons",
                            10, 0, 10));
    // Force Manual memory management
    //
    _muons->SetDirectory(0);
    _muons->GetXaxis()->SetTitle("N_{#mu}");
    _muons->Sumw2();

    _pthat.reset(new TH1D("pthat", "pT hat", 100, 0, 200));
    // Force Manual memory management
    //
    _pthat->SetDirectory(0);
    _pthat->GetXaxis()->SetTitle("#hat{p_{T}}");
    _pthat->Sumw2();

    _jets.reset(new TH1D("jets", "N jets",
                            10, 0, 10));
    // Force Manual memory management
    //
    _jets->SetDirectory(0);
    _jets->GetXaxis()->SetTitle("N_{jet}");
    _jets->Sumw2();

    _primary_vertices.reset(new TH1D("pvs", "Primary Vertices",
                                     10, 0, 10));
    // Forse Manual memory management
    //
    _primary_vertices->SetDirectory(0);
    _primary_vertices->GetXaxis()->SetTitle("N_{PV}");
    _primary_vertices->Sumw2();

    // Generic
    //
    _monitorMuons.reset(new MonitorLepton("allmus"));
    _monitorJets.reset(new MonitorJet("alljets"));
    _monitorJets->plot(MonitorBase::PT)->SetBins(230, 0, 230); 
    _monitorLeadingJet.reset(new MonitorJet("leadingjet"));

    // (n)
    //
    _monitorNMuonInJet.reset(new MonitorMuonInJet("n"));
    _n_primary_vertices.reset(new TH1D("npvs", "(n) Primary Vertices",
                                        10, 0, 10));
    // Forse Manual memory management
    //
    _n_primary_vertices->SetDirectory(0);
    _n_primary_vertices->GetXaxis()->SetTitle("(n) N_{PV}");
    _n_primary_vertices->Sumw2();

    _n_primary_vertices_z.reset(new TH1D("npvz", "(n) Primary Vertices Z",
                                        400, -20, 20));
    // Forse Manual memory management
    //
    _n_primary_vertices_z->SetDirectory(0);
    _n_primary_vertices_z->GetXaxis()->SetTitle("cm");
    _n_primary_vertices_z->Sumw2();

    // (p)
    //
    _monitorPMuonInJet.reset(new MonitorMuonInJet("p"));

    // change binning in the Delta plots
    //
    MonitorDelta *deltas[] = { _monitorNMuonInJet->delta(),
                               _monitorPMuonInJet->delta() };
    for(int i = 0; 2 > i; ++i)
    {
        deltas[i]->plot(MonitorDelta::PTREL)->SetBins(70, 0, 7);
    }

    _s8_selector = new S8Selector();
    _are_n_plots_filled = false;
}

void MonitorAnalyzer::treeDidLoad(const TreeInfo *,
                                  const TriggerCenter *trigger_center)
{
    _s8_selector->treeDidLoad(trigger_center);
}

void MonitorAnalyzer::eventDidLoad(const Event *event)
{
    using s8::Jets;
    using s8::Leptons;

    // check if event passes the selection
    //
    const Event *modified_event = (*_s8_selector)(event);
    if (!modified_event)
        return;

    event = modified_event;

    _pthat->Fill(event->gen()->ptHat());
    _muons->Fill(event->muons()->size());
    _jets->Fill(event->jets()->size());
    _primary_vertices->Fill(event->primaryVertices()->size());

    for(Leptons::const_iterator muon = event->muons()->begin();
        event->muons()->end() != muon;
        ++muon)

        _monitorMuons->fill(1, *muon);

    const Jet *leadingJet = 0;
    for(Jets::const_iterator jet = event->jets()->begin();
        event->jets()->end() != jet;
        ++jet)
    {
        _monitorJets->fill(1, *jet);

        if (leadingJet &&
            leadingJet->p4()->Pt() >= (*jet)->p4()->Pt())

            continue;

        leadingJet = *jet;
    }

    if (leadingJet)
        _monitorLeadingJet->fill(1, leadingJet);

    _are_n_plots_filled = false;
    _event = event;

    (*_muonInJet)(event);

    _event = 0;
}

void MonitorAnalyzer::print(std::ostream &out) const
{
    _s8_selector->print(out);
}

void MonitorAnalyzer::save(TDirectory *directory) const
{
    TDirectory *output = directory->mkdir("MonitorAnalyzer");
    if (!output)
    {
        cerr << "MonitorAnalyzer: Failed to create output folder: "
            << "no output is saved" << endl;

        return;
    }

    output->cd();

    saveGenericPlots(output);
    saveNPlots(output);
    savePPlots(output);
}

// MuonInJetOptionsDelegate interface
//
void MonitorAnalyzer::optionTagIsSet(const std::string &tag)
{
    *_muonInJetTagger << tag;
}

void MonitorAnalyzer::optionAwayTagIsSet(const std::string &tag)
{
    *_awayJetTagger << tag;
}

void MonitorAnalyzer::optionMuonPtIsSet(const Range &value)
{
    _n_muon_pt = value;
}

void MonitorAnalyzer::optionJetPtIsSet(const Range &value)
{
    _n_jet_pt = value;
}

void MonitorAnalyzer::optionJetEtaIsSet(const Range &value)
{
    _n_jet_eta = value;
}

// MuonInJetDelegate interface
//
bool MonitorAnalyzer::shouldSkipMuonInJetPlusAwayJet(const Lepton *muon,
                                                     const Jet *jet)
{
    return !isValueInRange(muon->p4()->Pt(), _n_muon_pt) ||
           !isValueInRange(jet->p4()->Pt(), _n_jet_pt) ||
           !isValueInRange(fabs(jet->p4()->Eta()), _n_jet_eta);
}

void MonitorAnalyzer::muonIsInJetPlusAwayJet(const Lepton *muon,
                                             const Jet *jet,
                                             const Event* event)
{
    double weight    = 1;
    double weight_pv = 1;
    const double jet_pt = jet->p4()->Pt();
    using std::min;

    if (_reweights &&
        isValueInRange(jet_pt, _n_jet_pt))

        weight = _reweights->GetBinContent(_reweights->FindBin(jet->p4()->Pt(),
                                                               jet->p4()->Eta()));
    if (_reweightsPV &&
        isValueInRange(jet_pt, _n_jet_pt)) {
        int mc_npvminus1 = min(event->gen()->bx_minus1(),34);
        int mc_npv0      = min(event->gen()->bx_zero(),34);
        int mc_npvplus1  = min(event->gen()->bx_plus1(),34);
        if (mc_npvminus1 || mc_npv0 || mc_npvplus1 )
            weight_pv = _reweightsPV->GetBinContent(mc_npvminus1,mc_npv0,mc_npvplus1);
    }

    _monitorNMuonInJet->fill(weight * weight_pv, muon, jet, _event->primaryVertices()->size());
    _monitorNMuonInJet->jet()->fillDiscriminator(weight * weight_pv,
            jet->btag(_muonInJetTagger->btag()));

    if (!_are_n_plots_filled)
    {
        _n_primary_vertices->Fill(_event->primaryVertices()->size(), weight);

        for(PrimaryVertices::const_iterator pv = _event->primaryVertices()->begin();
            _event->primaryVertices()->end() != pv;
            ++pv)
        {
            _n_primary_vertices_z->Fill((*pv)->vertex()->Z(), weight);
        }

        _are_n_plots_filled = true;
    }
}

void MonitorAnalyzer::muonIsInJetPlusTaggedAwayJet(const Lepton *muon,
                                                   const Jet *jet,
                                                   const Event * event)
{
    double weight    = 1;
    double weight_pv = 1;
    const double jet_pt = jet->p4()->Pt();
    using std::min;

    if (_reweights &&
        isValueInRange(jet_pt, _n_jet_pt))

        weight = _reweights->GetBinContent(_reweights->FindBin(jet->p4()->Pt(),
                                                               jet->p4()->Eta()));
    if (_reweightsPV &&
        isValueInRange(jet_pt, _n_jet_pt)) {
        int mc_npvminus1 = min(event->gen()->bx_minus1(),34);
        int mc_npv0      = min(event->gen()->bx_zero(),34);
        int mc_npvplus1  = min(event->gen()->bx_plus1(),34);
        if (mc_npvminus1 || mc_npv0 || mc_npvplus1 )
            weight_pv = _reweightsPV->GetBinContent(mc_npvminus1,mc_npv0,mc_npvplus1);
    }

    _monitorPMuonInJet->fill(weight * weight_pv, muon, jet, _event->primaryVertices()->size());
    _monitorPMuonInJet->jet()->fillDiscriminator(weight * weight_pv,
            jet->btag(_muonInJetTagger->btag()));
}

// PythiaOptionsDelegate interface
//
void MonitorAnalyzer::optionGluonSplittingIsSet(const GluonSplitting &value)
{
    _s8_selector->optionGluonSplittingIsSet(value);
}

void MonitorAnalyzer::optionPtHatIsSet(const Range &value)
{
    _s8_selector->optionPtHatIsSet(value);
}

// Trigger options
//
void MonitorAnalyzer::optionTriggerIsSet(const Trigger &trigger)
{
    _s8_selector->optionTriggerIsSet(trigger);
}

void MonitorAnalyzer::optionSimulateTriggerIsSet(const bool &value)
{
    _s8_selector->optionSimulateTriggerIsSet(value);
}

void MonitorAnalyzer::optionReweightTriggerIsSet(const std::string &filename)
{
    _reweight_file.reset(new TFile(filename.c_str(), "read"));
    if (!_reweight_file->IsOpen())
        throw runtime_error("Failed to open reweights file: " + filename);

    _reweights = dynamic_cast<TH2 *>(_reweight_file->Get("n_pt_eta"));
    if (!_reweights)
        throw runtime_error("Failed to extract reweights from: " + filename);
}

void MonitorAnalyzer::optionReweightTriggerPVIsSet(const std::string &filename)
{

  _reweight_pv_file.reset(new TFile(filename.c_str(), "read"));
  if (!_reweight_pv_file->IsOpen())
    throw runtime_error("Failed to open PVreweights file: " + filename);
  
  _reweightsPV = dynamic_cast<TH3D *>(_reweight_pv_file->Get("WHist"));
  
  if (!_reweightsPV)
    throw runtime_error("Failed to extract PVreweights from: " + filename);

}
// Misc Options
//
void MonitorAnalyzer::optionPrimaryVerticesIsSet(const Range &primary_vertices)
{
    _s8_selector->optionPrimaryVerticesIsSet(primary_vertices);
}


void MonitorAnalyzer::saveGenericPlots(TDirectory *directory) const
{
    TDirectory *subdir = directory->mkdir("generic");
    if (!subdir)
    {
        cerr << "Failed to create 'generic' TDirectory. Plots are not saved"
            << endl;

        return;
    }

    subdir->cd();

    _muons->Write();
    _jets->Write();
    _pthat->Write();
    _primary_vertices->Write();

    // Generic
    //
    _monitorMuons->save(subdir);
    _monitorJets->save(subdir);
    _monitorLeadingJet->save(subdir);

    directory->cd();
}

void MonitorAnalyzer::saveNPlots(TDirectory *directory) const
{
    TDirectory *subdir = directory->mkdir("n");
    if (!subdir)
    {
        cerr << "Failed to create 'n' TDirectory. Plots are not saved"
            << endl;

        return;
    }

    subdir->cd();

    // (n)
    //
    _monitorNMuonInJet->save(subdir);
    _n_primary_vertices->Write();
    _n_primary_vertices_z->Write();

    directory->cd();
}

void MonitorAnalyzer::savePPlots(TDirectory *directory) const
{
    TDirectory *subdir = directory->mkdir("p");
    if (!subdir)
    {
        cerr << "Failed to create 'p' TDirectory. Plots are not saved"
            << endl;

        return;
    }

    subdir->cd();

    // (p)
    //
    _monitorPMuonInJet->save(subdir);

    directory->cd();
}
Revision:	1.4
Committed:	Fri Dec 16 00:38:31 2011 UTC (13 years, 4 months ago) by meloam
Branch:	MAIN
CVS Tags:	V00-00-04, HEAD
Changes since 1.3:	+49 -13 lines
Log Message:	adding pu reweighting, new OPs
#	User	Rev	Content
1	samvel	1.1	/**
2			* MonitorAnalyzer
3			* s8
4			*
5			* Created by Samvel Khalatian on Nov 17, 2010
6			* Copyright 2010, All rights reserved
7			*/
8
9			#include <cmath>
10
11			#include <iostream>
12			#include <stdexcept>
13	meloam	1.4	#include <algorithm>
14	samvel	1.1
15			#include <TDirectory.h>
16			#include <TFile.h>
17			#include <TH1F.h>
18	meloam	1.4	#include <TH2F.h>
19	samvel	1.1	#include <TH3F.h>
20	meloam	1.4	#include <TH3D.h>
21	samvel	1.1	#include <TH1D.h>
22			#include <TLorentzVector.h>
23	samvel	1.2	#include <TVector3.h>
24	samvel	1.1
25	samvel	1.3	#include "s8/IO/interface/Event.h"
26			#include "Tree/System8/interface/S8GenEvent.h"
27			#include "Tree/System8/interface/S8Fwd.h"
28			#include "Tree/System8/interface/S8Jet.h"
29			#include "Tree/System8/interface/S8Lepton.h"
30			#include "Tree/System8/interface/S8PrimaryVertex.h"
31			#include "s8/Selector/interface/S8Selector.h"
32			#include "s8/Utility/interface/MonitorPlots.h"
33			#include "s8/Utility/interface/MuonInJet.h"
34			#include "s8/Utility/interface/TaggerOperatingPoint.h"
35	samvel	1.1
36	samvel	1.3	#include "s8/Analyzer/interface/MonitorAnalyzer.h"
37	samvel	1.1
38			using std::cerr;
39			using std::clog;
40			using std::cout;
41			using std::endl;
42			using std::runtime_error;
43
44			using s8::MonitorAnalyzer;
45
46			MonitorAnalyzer::MonitorAnalyzer() throw()
47			{
48			_s8_selector = 0;
49			_reweights = 0;
50	samvel	1.2	_event = 0;
51	samvel	1.1	}
52
53			MonitorAnalyzer::~MonitorAnalyzer() throw()
54			{
55			if (_s8_selector)
56			delete _s8_selector;
57			}
58
59
60
61			// Analyzer interface
62			//
63			void MonitorAnalyzer::init()
64			{
65			_muonInJetTagger.reset(new TaggerOperatingPoint());
66			_awayJetTagger.reset(new TaggerOperatingPoint());
67
68			_muonInJet.reset(new MuonInJet());
69			_muonInJet->setDelegate(this);
70			_muonInJet->setAwayJetTaggerOperatingPoint(_awayJetTagger.get());
71
72			_muons.reset(new TH1D("muons", "N muons",
73			10, 0, 10));
74			// Force Manual memory management
75			//
76			_muons->SetDirectory(0);
77			_muons->GetXaxis()->SetTitle("N_{#mu}");
78			_muons->Sumw2();
79
80			_pthat.reset(new TH1D("pthat", "pT hat", 100, 0, 200));
81			// Force Manual memory management
82			//
83			_pthat->SetDirectory(0);
84			_pthat->GetXaxis()->SetTitle("#hat{p_{T}}");
85			_pthat->Sumw2();
86
87			_jets.reset(new TH1D("jets", "N jets",
88			10, 0, 10));
89			// Force Manual memory management
90			//
91			_jets->SetDirectory(0);
92			_jets->GetXaxis()->SetTitle("N_{jet}");
93			_jets->Sumw2();
94
95			_primary_vertices.reset(new TH1D("pvs", "Primary Vertices",
96			10, 0, 10));
97			// Forse Manual memory management
98			//
99			_primary_vertices->SetDirectory(0);
100			_primary_vertices->GetXaxis()->SetTitle("N_{PV}");
101			_primary_vertices->Sumw2();
102
103			// Generic
104			//
105			_monitorMuons.reset(new MonitorLepton("allmus"));
106			_monitorJets.reset(new MonitorJet("alljets"));
107			_monitorJets->plot(MonitorBase::PT)->SetBins(230, 0, 230);
108			_monitorLeadingJet.reset(new MonitorJet("leadingjet"));
109
110			// (n)
111			//
112			_monitorNMuonInJet.reset(new MonitorMuonInJet("n"));
113	samvel	1.2	_n_primary_vertices.reset(new TH1D("npvs", "(n) Primary Vertices",
114			10, 0, 10));
115			// Forse Manual memory management
116			//
117			_n_primary_vertices->SetDirectory(0);
118			_n_primary_vertices->GetXaxis()->SetTitle("(n) N_{PV}");
119			_n_primary_vertices->Sumw2();
120
121			_n_primary_vertices_z.reset(new TH1D("npvz", "(n) Primary Vertices Z",
122			400, -20, 20));
123			// Forse Manual memory management
124			//
125			_n_primary_vertices_z->SetDirectory(0);
126			_n_primary_vertices_z->GetXaxis()->SetTitle("cm");
127			_n_primary_vertices_z->Sumw2();
128	samvel	1.1
129			// (p)
130			//
131			_monitorPMuonInJet.reset(new MonitorMuonInJet("p"));
132
133			// change binning in the Delta plots
134			//
135			MonitorDelta *deltas[] = { _monitorNMuonInJet->delta(),
136			_monitorPMuonInJet->delta() };
137			for(int i = 0; 2 > i; ++i)
138			{
139			deltas[i]->plot(MonitorDelta::PTREL)->SetBins(70, 0, 7);
140			}
141
142			_s8_selector = new S8Selector();
143	samvel	1.2	_are_n_plots_filled = false;
144	samvel	1.1	}
145
146			void MonitorAnalyzer::treeDidLoad(const TreeInfo *,
147			const TriggerCenter *trigger_center)
148			{
149			_s8_selector->treeDidLoad(trigger_center);
150			}
151
152			void MonitorAnalyzer::eventDidLoad(const Event *event)
153			{
154			using s8::Jets;
155			using s8::Leptons;
156
157			// check if event passes the selection
158			//
159			const Event modified_event = (_s8_selector)(event);
160			if (!modified_event)
161			return;
162
163			event = modified_event;
164
165			_pthat->Fill(event->gen()->ptHat());
166			_muons->Fill(event->muons()->size());
167			_jets->Fill(event->jets()->size());
168			_primary_vertices->Fill(event->primaryVertices()->size());
169
170			for(Leptons::const_iterator muon = event->muons()->begin();
171			event->muons()->end() != muon;
172			++muon)
173
174			_monitorMuons->fill(1, *muon);
175
176			const Jet *leadingJet = 0;
177			for(Jets::const_iterator jet = event->jets()->begin();
178			event->jets()->end() != jet;
179			++jet)
180			{
181			_monitorJets->fill(1, *jet);
182
183			if (leadingJet &&
184			leadingJet->p4()->Pt() >= (*jet)->p4()->Pt())
185
186			continue;
187
188			leadingJet = *jet;
189			}
190
191			if (leadingJet)
192			_monitorLeadingJet->fill(1, leadingJet);
193
194	samvel	1.2	_are_n_plots_filled = false;
195			_event = event;
196
197	samvel	1.1	(*_muonInJet)(event);
198	samvel	1.2
199			_event = 0;
200	samvel	1.1	}
201
202			void MonitorAnalyzer::print(std::ostream &out) const
203			{
204			_s8_selector->print(out);
205			}
206
207			void MonitorAnalyzer::save(TDirectory *directory) const
208			{
209			TDirectory *output = directory->mkdir("MonitorAnalyzer");
210			if (!output)
211			{
212			cerr << "MonitorAnalyzer: Failed to create output folder: "
213			<< "no output is saved" << endl;
214
215			return;
216			}
217
218			output->cd();
219
220			saveGenericPlots(output);
221			saveNPlots(output);
222			savePPlots(output);
223			}
224
225			// MuonInJetOptionsDelegate interface
226			//
227			void MonitorAnalyzer::optionTagIsSet(const std::string &tag)
228			{
229			*_muonInJetTagger << tag;
230			}
231
232			void MonitorAnalyzer::optionAwayTagIsSet(const std::string &tag)
233			{
234			*_awayJetTagger << tag;
235			}
236
237			void MonitorAnalyzer::optionMuonPtIsSet(const Range &value)
238			{
239			_n_muon_pt = value;
240			}
241
242			void MonitorAnalyzer::optionJetPtIsSet(const Range &value)
243			{
244			_n_jet_pt = value;
245			}
246
247			void MonitorAnalyzer::optionJetEtaIsSet(const Range &value)
248			{
249			_n_jet_eta = value;
250			}
251
252			// MuonInJetDelegate interface
253			//
254			bool MonitorAnalyzer::shouldSkipMuonInJetPlusAwayJet(const Lepton *muon,
255			const Jet *jet)
256			{
257			return !isValueInRange(muon->p4()->Pt(), _n_muon_pt) \|\|
258			!isValueInRange(jet->p4()->Pt(), _n_jet_pt) \|\|
259			!isValueInRange(fabs(jet->p4()->Eta()), _n_jet_eta);
260			}
261
262			void MonitorAnalyzer::muonIsInJetPlusAwayJet(const Lepton *muon,
263	meloam	1.4	const Jet *jet,
264			const Event* event)
265	samvel	1.1	{
266	meloam	1.4	double weight = 1;
267			double weight_pv = 1;
268	samvel	1.1	const double jet_pt = jet->p4()->Pt();
269	meloam	1.4	using std::min;
270	samvel	1.1
271			if (_reweights &&
272			isValueInRange(jet_pt, _n_jet_pt))
273
274			weight = _reweights->GetBinContent(_reweights->FindBin(jet->p4()->Pt(),
275	meloam	1.4	jet->p4()->Eta()));
276			if (_reweightsPV &&
277			isValueInRange(jet_pt, _n_jet_pt)) {
278			int mc_npvminus1 = min(event->gen()->bx_minus1(),34);
279			int mc_npv0 = min(event->gen()->bx_zero(),34);
280			int mc_npvplus1 = min(event->gen()->bx_plus1(),34);
281			if (mc_npvminus1 \|\| mc_npv0 \|\| mc_npvplus1 )
282			weight_pv = _reweightsPV->GetBinContent(mc_npvminus1,mc_npv0,mc_npvplus1);
283			}
284	samvel	1.1
285	meloam	1.4	_monitorNMuonInJet->fill(weight * weight_pv, muon, jet, _event->primaryVertices()->size());
286			_monitorNMuonInJet->jet()->fillDiscriminator(weight * weight_pv,
287	samvel	1.1	jet->btag(_muonInJetTagger->btag()));
288	samvel	1.2
289			if (!_are_n_plots_filled)
290			{
291			_n_primary_vertices->Fill(_event->primaryVertices()->size(), weight);
292
293			for(PrimaryVertices::const_iterator pv = _event->primaryVertices()->begin();
294			_event->primaryVertices()->end() != pv;
295			++pv)
296			{
297			_n_primary_vertices_z->Fill((*pv)->vertex()->Z(), weight);
298			}
299
300			_are_n_plots_filled = true;
301			}
302	samvel	1.1	}
303
304			void MonitorAnalyzer::muonIsInJetPlusTaggedAwayJet(const Lepton *muon,
305	meloam	1.4	const Jet *jet,
306			const Event * event)
307	samvel	1.1	{
308	meloam	1.4	double weight = 1;
309			double weight_pv = 1;
310	samvel	1.1	const double jet_pt = jet->p4()->Pt();
311	meloam	1.4	using std::min;
312	samvel	1.1
313			if (_reweights &&
314			isValueInRange(jet_pt, _n_jet_pt))
315
316			weight = _reweights->GetBinContent(_reweights->FindBin(jet->p4()->Pt(),
317	meloam	1.4	jet->p4()->Eta()));
318			if (_reweightsPV &&
319			isValueInRange(jet_pt, _n_jet_pt)) {
320			int mc_npvminus1 = min(event->gen()->bx_minus1(),34);
321			int mc_npv0 = min(event->gen()->bx_zero(),34);
322			int mc_npvplus1 = min(event->gen()->bx_plus1(),34);
323			if (mc_npvminus1 \|\| mc_npv0 \|\| mc_npvplus1 )
324			weight_pv = _reweightsPV->GetBinContent(mc_npvminus1,mc_npv0,mc_npvplus1);
325			}
326	samvel	1.1
327	meloam	1.4	_monitorPMuonInJet->fill(weight * weight_pv, muon, jet, _event->primaryVertices()->size());
328			_monitorPMuonInJet->jet()->fillDiscriminator(weight * weight_pv,
329	samvel	1.1	jet->btag(_muonInJetTagger->btag()));
330			}
331
332			// PythiaOptionsDelegate interface
333			//
334			void MonitorAnalyzer::optionGluonSplittingIsSet(const GluonSplitting &value)
335			{
336			_s8_selector->optionGluonSplittingIsSet(value);
337			}
338
339			void MonitorAnalyzer::optionPtHatIsSet(const Range &value)
340			{
341			_s8_selector->optionPtHatIsSet(value);
342			}
343
344			// Trigger options
345			//
346			void MonitorAnalyzer::optionTriggerIsSet(const Trigger &trigger)
347			{
348			_s8_selector->optionTriggerIsSet(trigger);
349			}
350
351			void MonitorAnalyzer::optionSimulateTriggerIsSet(const bool &value)
352			{
353			_s8_selector->optionSimulateTriggerIsSet(value);
354			}
355
356			void MonitorAnalyzer::optionReweightTriggerIsSet(const std::string &filename)
357			{
358			_reweight_file.reset(new TFile(filename.c_str(), "read"));
359			if (!_reweight_file->IsOpen())
360			throw runtime_error("Failed to open reweights file: " + filename);
361
362	meloam	1.4	_reweights = dynamic_cast<TH2 *>(_reweight_file->Get("n_pt_eta"));
363	samvel	1.1	if (!_reweights)
364			throw runtime_error("Failed to extract reweights from: " + filename);
365			}
366
367	meloam	1.4	void MonitorAnalyzer::optionReweightTriggerPVIsSet(const std::string &filename)
368			{
369
370			_reweight_pv_file.reset(new TFile(filename.c_str(), "read"));
371			if (!_reweight_pv_file->IsOpen())
372			throw runtime_error("Failed to open PVreweights file: " + filename);
373
374			_reweightsPV = dynamic_cast<TH3D *>(_reweight_pv_file->Get("WHist"));
375
376			if (!_reweightsPV)
377			throw runtime_error("Failed to extract PVreweights from: " + filename);
378
379			}
380	samvel	1.1	// Misc Options
381			//
382			void MonitorAnalyzer::optionPrimaryVerticesIsSet(const Range &primary_vertices)
383			{
384			_s8_selector->optionPrimaryVerticesIsSet(primary_vertices);
385			}
386
387
388
389			void MonitorAnalyzer::saveGenericPlots(TDirectory *directory) const
390			{
391			TDirectory *subdir = directory->mkdir("generic");
392			if (!subdir)
393			{
394			cerr << "Failed to create 'generic' TDirectory. Plots are not saved"
395			<< endl;
396
397			return;
398			}
399
400			subdir->cd();
401
402			_muons->Write();
403			_jets->Write();
404			_pthat->Write();
405			_primary_vertices->Write();
406
407			// Generic
408			//
409			_monitorMuons->save(subdir);
410			_monitorJets->save(subdir);
411			_monitorLeadingJet->save(subdir);
412
413			directory->cd();
414			}
415
416			void MonitorAnalyzer::saveNPlots(TDirectory *directory) const
417			{
418			TDirectory *subdir = directory->mkdir("n");
419			if (!subdir)
420			{
421			cerr << "Failed to create 'n' TDirectory. Plots are not saved"
422			<< endl;
423
424			return;
425			}
426
427			subdir->cd();
428
429			// (n)
430			//
431			_monitorNMuonInJet->save(subdir);
432	samvel	1.2	_n_primary_vertices->Write();
433			_n_primary_vertices_z->Write();
434	samvel	1.1
435			directory->cd();
436			}
437
438			void MonitorAnalyzer::savePPlots(TDirectory *directory) const
439			{
440			TDirectory *subdir = directory->mkdir("p");
441			if (!subdir)
442			{
443			cerr << "Failed to create 'p' TDirectory. Plots are not saved"
444			<< endl;
445
446			return;
447			}
448
449			subdir->cd();
450
451			// (p)
452			//
453			_monitorPMuonInJet->save(subdir);
454
455			directory->cd();
456			}