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 <TDirectory.h>
#include <TFile.h>
#include <TH1F.h>
#include <TH3F.h>
#include <TH1D.h>
#include <TLorentzVector.h>

#include "IO/interface/Event.h"
#include "S8Tree/interface/S8GenEvent.h"
#include "S8Tree/interface/S8Fwd.h"
#include "S8Tree/interface/S8Jet.h"
#include "S8Tree/interface/S8Lepton.h"
#include "Selector/interface/S8Selector.h"
#include "Utility/interface/MonitorPlots.h"
#include "Utility/interface/MuonInJet.h"
#include "Utility/interface/TaggerOperatingPoint.h"

#include "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;
}

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"));

    // (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();
    _n_primary_vertices = 0;
}

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());
    _n_primary_vertices = 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);

    (*_muonInJet)(event);
}

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)
{
    double weight = 1;
    const double jet_pt = jet->p4()->Pt();

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

        weight = _reweights->GetBinContent(_reweights->FindBin(jet->p4()->Pt(),
                                                               jet->p4()->Eta(),
                                                               _n_primary_vertices));

    _monitorNMuonInJet->fill(weight, muon, jet, _n_primary_vertices);
    _monitorNMuonInJet->jet()->fillDiscriminator(weight,
            jet->btag(_muonInJetTagger->btag()));
}

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

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

        weight = _reweights->GetBinContent(_reweights->FindBin(jet->p4()->Pt(),
                                                               jet->p4()->Eta(),
                                                               _n_primary_vertices));

    _monitorPMuonInJet->fill(weight, muon, jet, _n_primary_vertices);
    _monitorPMuonInJet->jet()->fillDiscriminator(weight,
            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<TH3 *>(_reweight_file->Get("n_pt_eta_pv"));
    if (!_reweights)
        throw runtime_error("Failed to extract reweights 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);

    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.1
Committed:	Fri May 6 14:31:20 2011 UTC (14 years ago) by samvel
Branch:	MAIN
Log Message:	Import Analyzer
#	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
14			#include <TDirectory.h>
15			#include <TFile.h>
16			#include <TH1F.h>
17			#include <TH3F.h>
18			#include <TH1D.h>
19			#include <TLorentzVector.h>
20
21			#include "IO/interface/Event.h"
22			#include "S8Tree/interface/S8GenEvent.h"
23			#include "S8Tree/interface/S8Fwd.h"
24			#include "S8Tree/interface/S8Jet.h"
25			#include "S8Tree/interface/S8Lepton.h"
26			#include "Selector/interface/S8Selector.h"
27			#include "Utility/interface/MonitorPlots.h"
28			#include "Utility/interface/MuonInJet.h"
29			#include "Utility/interface/TaggerOperatingPoint.h"
30
31			#include "Analyzer/interface/MonitorAnalyzer.h"
32
33			using std::cerr;
34			using std::clog;
35			using std::cout;
36			using std::endl;
37			using std::runtime_error;
38
39			using s8::MonitorAnalyzer;
40
41			MonitorAnalyzer::MonitorAnalyzer() throw()
42			{
43			_s8_selector = 0;
44			_reweights = 0;
45			}
46
47			MonitorAnalyzer::~MonitorAnalyzer() throw()
48			{
49			if (_s8_selector)
50			delete _s8_selector;
51			}
52
53
54
55			// Analyzer interface
56			//
57			void MonitorAnalyzer::init()
58			{
59			_muonInJetTagger.reset(new TaggerOperatingPoint());
60			_awayJetTagger.reset(new TaggerOperatingPoint());
61
62			_muonInJet.reset(new MuonInJet());
63			_muonInJet->setDelegate(this);
64			_muonInJet->setAwayJetTaggerOperatingPoint(_awayJetTagger.get());
65
66			_muons.reset(new TH1D("muons", "N muons",
67			10, 0, 10));
68			// Force Manual memory management
69			//
70			_muons->SetDirectory(0);
71			_muons->GetXaxis()->SetTitle("N_{#mu}");
72			_muons->Sumw2();
73
74			_pthat.reset(new TH1D("pthat", "pT hat", 100, 0, 200));
75			// Force Manual memory management
76			//
77			_pthat->SetDirectory(0);
78			_pthat->GetXaxis()->SetTitle("#hat{p_{T}}");
79			_pthat->Sumw2();
80
81			_jets.reset(new TH1D("jets", "N jets",
82			10, 0, 10));
83			// Force Manual memory management
84			//
85			_jets->SetDirectory(0);
86			_jets->GetXaxis()->SetTitle("N_{jet}");
87			_jets->Sumw2();
88
89			_primary_vertices.reset(new TH1D("pvs", "Primary Vertices",
90			10, 0, 10));
91			// Forse Manual memory management
92			//
93			_primary_vertices->SetDirectory(0);
94			_primary_vertices->GetXaxis()->SetTitle("N_{PV}");
95			_primary_vertices->Sumw2();
96
97			// Generic
98			//
99			_monitorMuons.reset(new MonitorLepton("allmus"));
100			_monitorJets.reset(new MonitorJet("alljets"));
101			_monitorJets->plot(MonitorBase::PT)->SetBins(230, 0, 230);
102			_monitorLeadingJet.reset(new MonitorJet("leadingjet"));
103
104			// (n)
105			//
106			_monitorNMuonInJet.reset(new MonitorMuonInJet("n"));
107
108			// (p)
109			//
110			_monitorPMuonInJet.reset(new MonitorMuonInJet("p"));
111
112			// change binning in the Delta plots
113			//
114			MonitorDelta *deltas[] = { _monitorNMuonInJet->delta(),
115			_monitorPMuonInJet->delta() };
116			for(int i = 0; 2 > i; ++i)
117			{
118			deltas[i]->plot(MonitorDelta::PTREL)->SetBins(70, 0, 7);
119			}
120
121			_s8_selector = new S8Selector();
122			_n_primary_vertices = 0;
123			}
124
125			void MonitorAnalyzer::treeDidLoad(const TreeInfo *,
126			const TriggerCenter *trigger_center)
127			{
128			_s8_selector->treeDidLoad(trigger_center);
129			}
130
131			void MonitorAnalyzer::eventDidLoad(const Event *event)
132			{
133			using s8::Jets;
134			using s8::Leptons;
135
136			// check if event passes the selection
137			//
138			const Event modified_event = (_s8_selector)(event);
139			if (!modified_event)
140			return;
141
142			event = modified_event;
143
144			_pthat->Fill(event->gen()->ptHat());
145			_muons->Fill(event->muons()->size());
146			_jets->Fill(event->jets()->size());
147			_primary_vertices->Fill(event->primaryVertices()->size());
148			_n_primary_vertices = event->primaryVertices()->size();
149
150			for(Leptons::const_iterator muon = event->muons()->begin();
151			event->muons()->end() != muon;
152			++muon)
153
154			_monitorMuons->fill(1, *muon);
155
156			const Jet *leadingJet = 0;
157			for(Jets::const_iterator jet = event->jets()->begin();
158			event->jets()->end() != jet;
159			++jet)
160			{
161			_monitorJets->fill(1, *jet);
162
163			if (leadingJet &&
164			leadingJet->p4()->Pt() >= (*jet)->p4()->Pt())
165
166			continue;
167
168			leadingJet = *jet;
169			}
170
171			if (leadingJet)
172			_monitorLeadingJet->fill(1, leadingJet);
173
174			(*_muonInJet)(event);
175			}
176
177			void MonitorAnalyzer::print(std::ostream &out) const
178			{
179			_s8_selector->print(out);
180			}
181
182			void MonitorAnalyzer::save(TDirectory *directory) const
183			{
184			TDirectory *output = directory->mkdir("MonitorAnalyzer");
185			if (!output)
186			{
187			cerr << "MonitorAnalyzer: Failed to create output folder: "
188			<< "no output is saved" << endl;
189
190			return;
191			}
192
193			output->cd();
194
195			saveGenericPlots(output);
196			saveNPlots(output);
197			savePPlots(output);
198			}
199
200			// MuonInJetOptionsDelegate interface
201			//
202			void MonitorAnalyzer::optionTagIsSet(const std::string &tag)
203			{
204			*_muonInJetTagger << tag;
205			}
206
207			void MonitorAnalyzer::optionAwayTagIsSet(const std::string &tag)
208			{
209			*_awayJetTagger << tag;
210			}
211
212			void MonitorAnalyzer::optionMuonPtIsSet(const Range &value)
213			{
214			_n_muon_pt = value;
215			}
216
217			void MonitorAnalyzer::optionJetPtIsSet(const Range &value)
218			{
219			_n_jet_pt = value;
220			}
221
222			void MonitorAnalyzer::optionJetEtaIsSet(const Range &value)
223			{
224			_n_jet_eta = value;
225			}
226
227			// MuonInJetDelegate interface
228			//
229			bool MonitorAnalyzer::shouldSkipMuonInJetPlusAwayJet(const Lepton *muon,
230			const Jet *jet)
231			{
232			return !isValueInRange(muon->p4()->Pt(), _n_muon_pt) \|\|
233			!isValueInRange(jet->p4()->Pt(), _n_jet_pt) \|\|
234			!isValueInRange(fabs(jet->p4()->Eta()), _n_jet_eta);
235			}
236
237			void MonitorAnalyzer::muonIsInJetPlusAwayJet(const Lepton *muon,
238			const Jet *jet)
239			{
240			double weight = 1;
241			const double jet_pt = jet->p4()->Pt();
242
243			if (_reweights &&
244			isValueInRange(jet_pt, _n_jet_pt))
245
246			weight = _reweights->GetBinContent(_reweights->FindBin(jet->p4()->Pt(),
247			jet->p4()->Eta(),
248			_n_primary_vertices));
249
250			_monitorNMuonInJet->fill(weight, muon, jet, _n_primary_vertices);
251			_monitorNMuonInJet->jet()->fillDiscriminator(weight,
252			jet->btag(_muonInJetTagger->btag()));
253			}
254
255			void MonitorAnalyzer::muonIsInJetPlusTaggedAwayJet(const Lepton *muon,
256			const Jet *jet)
257			{
258			double weight = 1;
259			const double jet_pt = jet->p4()->Pt();
260
261			if (_reweights &&
262			isValueInRange(jet_pt, _n_jet_pt))
263
264			weight = _reweights->GetBinContent(_reweights->FindBin(jet->p4()->Pt(),
265			jet->p4()->Eta(),
266			_n_primary_vertices));
267
268			_monitorPMuonInJet->fill(weight, muon, jet, _n_primary_vertices);
269			_monitorPMuonInJet->jet()->fillDiscriminator(weight,
270			jet->btag(_muonInJetTagger->btag()));
271			}
272
273			// PythiaOptionsDelegate interface
274			//
275			void MonitorAnalyzer::optionGluonSplittingIsSet(const GluonSplitting &value)
276			{
277			_s8_selector->optionGluonSplittingIsSet(value);
278			}
279
280			void MonitorAnalyzer::optionPtHatIsSet(const Range &value)
281			{
282			_s8_selector->optionPtHatIsSet(value);
283			}
284
285			// Trigger options
286			//
287			void MonitorAnalyzer::optionTriggerIsSet(const Trigger &trigger)
288			{
289			_s8_selector->optionTriggerIsSet(trigger);
290			}
291
292			void MonitorAnalyzer::optionSimulateTriggerIsSet(const bool &value)
293			{
294			_s8_selector->optionSimulateTriggerIsSet(value);
295			}
296
297			void MonitorAnalyzer::optionReweightTriggerIsSet(const std::string &filename)
298			{
299			_reweight_file.reset(new TFile(filename.c_str(), "read"));
300			if (!_reweight_file->IsOpen())
301			throw runtime_error("Failed to open reweights file: " + filename);
302
303			_reweights = dynamic_cast<TH3 *>(_reweight_file->Get("n_pt_eta_pv"));
304			if (!_reweights)
305			throw runtime_error("Failed to extract reweights from: " + filename);
306			}
307
308			// Misc Options
309			//
310			void MonitorAnalyzer::optionPrimaryVerticesIsSet(const Range &primary_vertices)
311			{
312			_s8_selector->optionPrimaryVerticesIsSet(primary_vertices);
313			}
314
315
316
317			void MonitorAnalyzer::saveGenericPlots(TDirectory *directory) const
318			{
319			TDirectory *subdir = directory->mkdir("generic");
320			if (!subdir)
321			{
322			cerr << "Failed to create 'generic' TDirectory. Plots are not saved"
323			<< endl;
324
325			return;
326			}
327
328			subdir->cd();
329
330			_muons->Write();
331			_jets->Write();
332			_pthat->Write();
333			_primary_vertices->Write();
334
335			// Generic
336			//
337			_monitorMuons->save(subdir);
338			_monitorJets->save(subdir);
339			_monitorLeadingJet->save(subdir);
340
341			directory->cd();
342			}
343
344			void MonitorAnalyzer::saveNPlots(TDirectory *directory) const
345			{
346			TDirectory *subdir = directory->mkdir("n");
347			if (!subdir)
348			{
349			cerr << "Failed to create 'n' TDirectory. Plots are not saved"
350			<< endl;
351
352			return;
353			}
354
355			subdir->cd();
356
357			// (n)
358			//
359			_monitorNMuonInJet->save(subdir);
360
361			directory->cd();
362			}
363
364			void MonitorAnalyzer::savePPlots(TDirectory *directory) const
365			{
366			TDirectory *subdir = directory->mkdir("p");
367			if (!subdir)
368			{
369			cerr << "Failed to create 'p' TDirectory. Plots are not saved"
370			<< endl;
371
372			return;
373			}
374
375			subdir->cd();
376
377			// (p)
378			//
379			_monitorPMuonInJet->save(subdir);
380
381			directory->cd();
382			}