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
#	Content
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	#include <algorithm>
14
15	#include <TDirectory.h>
16	#include <TFile.h>
17	#include <TH1F.h>
18	#include <TH2F.h>
19	#include <TH3F.h>
20	#include <TH3D.h>
21	#include <TH1D.h>
22	#include <TLorentzVector.h>
23	#include <TVector3.h>
24
25	#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
36	#include "s8/Analyzer/interface/MonitorAnalyzer.h"
37
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	_event = 0;
51	}
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	_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
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	_are_n_plots_filled = false;
144	}
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	_are_n_plots_filled = false;
195	_event = event;
196
197	(*_muonInJet)(event);
198
199	_event = 0;
200	}
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	const Jet *jet,
264	const Event* event)
265	{
266	double weight = 1;
267	double weight_pv = 1;
268	const double jet_pt = jet->p4()->Pt();
269	using std::min;
270
271	if (_reweights &&
272	isValueInRange(jet_pt, _n_jet_pt))
273
274	weight = _reweights->GetBinContent(_reweights->FindBin(jet->p4()->Pt(),
275	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
285	_monitorNMuonInJet->fill(weight * weight_pv, muon, jet, _event->primaryVertices()->size());
286	_monitorNMuonInJet->jet()->fillDiscriminator(weight * weight_pv,
287	jet->btag(_muonInJetTagger->btag()));
288
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	}
303
304	void MonitorAnalyzer::muonIsInJetPlusTaggedAwayJet(const Lepton *muon,
305	const Jet *jet,
306	const Event * event)
307	{
308	double weight = 1;
309	double weight_pv = 1;
310	const double jet_pt = jet->p4()->Pt();
311	using std::min;
312
313	if (_reweights &&
314	isValueInRange(jet_pt, _n_jet_pt))
315
316	weight = _reweights->GetBinContent(_reweights->FindBin(jet->p4()->Pt(),
317	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
327	_monitorPMuonInJet->fill(weight * weight_pv, muon, jet, _event->primaryVertices()->size());
328	_monitorPMuonInJet->jet()->fillDiscriminator(weight * weight_pv,
329	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	_reweights = dynamic_cast<TH2 *>(_reweight_file->Get("n_pt_eta"));
363	if (!_reweights)
364	throw runtime_error("Failed to extract reweights from: " + filename);
365	}
366
367	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	// 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	_n_primary_vertices->Write();
433	_n_primary_vertices_z->Write();
434
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	}