ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc

Comparing UserCode/OSUT3Analysis/AnaTools/src/OSUAnalysis.cc (file contents):
Revision 1.1 by ahart, Thu Jun 14 17:42:31 2012 UTC vs.
Revision 1.2 by ahart, Mon Jul 23 09:30:37 2012 UTC

#	Line 1 | Line 1
1		#include "OSUT3Analysis/AnaTools/interface/OSUAnalysis.h"
2
3		OSUAnalysis::OSUAnalysis (const edm::ParameterSet &cfg) :
4	<	// Retrieve parameters from the configuration file
4	>	// Retrieve parameters from the configuration file.
5		muons_ (cfg.getParameter<edm::InputTag> ("muons")),
6		maxEta_ (cfg.getParameter<double> ("maxEta")),
7		minPt_ (cfg.getParameter<double> ("minPt"))
8		{
9	<	// Register cut names
10	<	cutNames_.push_back ("minMuons");
9	>	// Construct CutFlow objects. These store the results of cut decisions and
10	>	// handle filling cut flow histograms.
11	>	cuts_ = new CutFlow (hists1D_, fs_);
12	>	// muonCuts_ = new CutFlow (hists1D_, fs_, "muons");
13
14	<	// Create additional histograms
15	<	hists1D_["muonPT"] = fs_->make<TH1D> ("muonPt", ";p_{T} (GeV/c)", 100, 0.0, 300.0);
14	<	hists1D_["muonEta"] = fs_->make<TH1D> ("muonEta", ";#eta", 100, -3.0, 3.0);
14	>	// Create additional histograms.
15	>	TH1::SetDefaultSumw2 ();
16		hists1D_["muonPhi"] = fs_->make<TH1D> ("muonPhi", ";#phi", 100, -5.0, 5.0);
17	+	hists1D_["muonEta"] = fs_->make<TH1D> ("muonEta", ";#eta", 100, -3.0, 3.0);
18	+	hists1D_["muonPT"] = fs_->make<TH1D> ("muonPt", ";p_{T} (GeV/c)", 100, 0.0, 300.0);
19		hists1D_["dimuonMass"] = fs_->make<TH1D> ("dimuonMass", ";m_{#mu^{+} #mu^{-}} (GeV/c^{2})", 90, 30.0, 120.0);
20	<
21	<	TH1::SetDefaultSumw2 ();
22	<	hists1D_["cutFlow"] = fs_->make<TH1D> ("cutFlow", "", cutNames_.size () + 1, 0.0, cutNames_.size () + 1);
23	<	hists1D_["selection"] = fs_->make<TH1D> ("selection", "", cutNames_.size () + 1, 0.0, cutNames_.size () + 1);
21	<	hists1D_["complementarySelection"] = fs_->make<TH1D> ("complementarySelection", "", cutNames_.size () + 1, 0.0, cutNames_.size () + 1);
22	<
23	<	hists1D_["cutFlow"]->GetXaxis ()->SetBinLabel (1, "Total Events");
24	<	hists1D_["selection"]->GetXaxis ()->SetBinLabel (1, "Total Events");
25	<	hists1D_["complementarySelection"]->GetXaxis ()->SetBinLabel (1, "Total Events");
26	<	for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
27	<	{
28	<	hists1D_["cutFlow"]->GetXaxis ()->SetBinLabel (cut - cutNames_.begin () + 2, cut->c_str ());
29	<	hists1D_["selection"]->GetXaxis ()->SetBinLabel (cut - cutNames_.begin () + 2, cut->c_str ());
30	<	hists1D_["complementarySelection"]->GetXaxis ()->SetBinLabel (cut - cutNames_.begin () + 2, cut->c_str ());
31	<	}
32	<	}
33	<
34	<	void
35	<	OSUAnalysis::beginJob ()
36	<	{
37	<	sw_.Start ();
20	>	hists1D_["backMuonPhi"] = fs_->make<TH1D> ("backMuonPhi", ";#phi", 100, -5.0, 5.0);
21	>	hists1D_["backMuonEta"] = fs_->make<TH1D> ("backMuonEta", ";#eta", 100, -3.0, 3.0);
22	>	hists1D_["backMuonPT"] = fs_->make<TH1D> ("backMuonPt", ";p_{T} (GeV/c)", 100, 0.0, 300.0);
23	>	hists1D_["backDimuonMass"] = fs_->make<TH1D> ("backDimuonMass", ";m_{#mu^{+} #mu^{-}} (GeV/c^{2})", 90, 30.0, 120.0);
24		}
25
26	<	void
41	<	OSUAnalysis::endJob ()
26	>	OSUAnalysis::~OSUAnalysis ()
27		{
28	<	sw_.Stop ();
29	<	outputTime ();
30	<	outputCutFlow ();
28	>	// Destroying the CutFlow objects causes the cut flow numbers to be printed
29	>	// to standard output, as well as time information.
30	>	delete cuts_;
31		}
32
33		void
34		OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
35		{
36	<	// Retrieve necessary collections from the event
36	>	// Retrieve necessary collections from the event.
37		edm::Handle<BNmuonCollection> muons;
38		event.getByLabel (muons_, muons);
39
40	<	// Perform selection on objects within the event
41	<	vector<BNmuonCollection::const_iterator> goodMuons;
42	<	selectGoodMuons (muons, goodMuons);
43	<
44	<	// Perform selection calculations and store the results in the map cuts_. The
45	<	// keys to cuts_ must match the strings in cutNames_
46	<	cuts_["minMuons"] = goodMuons.size () >= 2;
40	>	// Perform selection on objects within the event. The full collection should
41	>	// only be gone through once.
42	>	vector<BNmuonCollection::const_iterator> leadMuons;
43	>	selectLeadMuons (muons, leadMuons);
44	>
45	>	// Calculate cut decisions and store the results in the CutFlow object, just
46	>	// like a map.
47	>	cuts_->at ("minMuons") = leadMuons.size () == 2;
48	>	if (cuts_->at ("minMuons"))
49	>	{
50	>	cuts_->at ("maxEta") = leadMuons.at (0)->eta < maxEta_ && leadMuons.at (1)->eta < maxEta_;
51	>	cuts_->at ("minPT") = leadMuons.at (0)->pt > minPt_ && leadMuons.at (1)->pt < minPt_;
52	>	cuts_->at ("oppositeSign") = leadMuons.at (0)->charge * leadMuons.at (1)->charge < 0;
53	>	}
54	>	else
55	>	cuts_->at ("maxEta") = cuts_->at ("minPT") = cuts_->at ("oppositeSign") = false;
56
57	<	fillCutFlow ();
57	>	// This causes the cut flow histograms to be filled.
58	>	cuts_->fillCutFlow ();
59
60	<	// Fill histograms while applying desired cuts.
61	<	for (vector<BNmuonCollection::const_iterator>::const_iterator muon1 = goodMuons.begin (); muon1 != goodMuons.end (); muon1++)
60	>	// Fill additional histograms while applying desired cuts.
61	>	if (cuts_->at ("minMuons"))
62		{
63	<	hists1D_["muonPT"]->Fill ((*muon1)->pt);
69	<	hists1D_["muonEta"]->Fill ((*muon1)->eta);
70	<	hists1D_["muonPhi"]->Fill ((*muon1)->phi);
71	<	for (vector<BNmuonCollection::const_iterator>::const_iterator muon2 = muon1 + 1; muon2 != goodMuons.end (); muon2++)
63	>	if (cuts_->at ("oppositeSign"))
64		{
65	<	if ((*muon1)->charge * (*muon2)->charge > 0)
66	<	continue;
67	<	hists1D_["dimuonMass"]->Fill (dimuonMass (*muon1, *muon2));
65	>	if (cuts_->at ("minPT"))
66	>	{
67	>	hists1D_["muonEta"]->Fill (leadMuons.at (0)->eta);
68	>	hists1D_["muonEta"]->Fill (leadMuons.at (1)->eta);
69	>	}
70	>	if (cuts_->at ("maxEta"))
71	>	{
72	>	hists1D_["muonPT"]->Fill (leadMuons.at (0)->pt);
73	>	hists1D_["muonPT"]->Fill (leadMuons.at (1)->pt);
74	>	}
75	>	if (cuts_->at ("minPT") && cuts_->at ("maxEta"))
76	>	{
77	>	hists1D_["muonPhi"]->Fill (leadMuons.at (0)->phi);
78	>	hists1D_["muonPhi"]->Fill (leadMuons.at (1)->phi);
79	>	hists1D_["dimuonMass"]->Fill (dimuonMass (leadMuons.at (0), leadMuons.at (1)));
80	>	}
81	>	}
82	>	else
83	>	{
84	>	if (cuts_->at ("minPT"))
85	>	{
86	>	hists1D_["backMuonEta"]->Fill (leadMuons.at (0)->eta);
87	>	hists1D_["backMuonEta"]->Fill (leadMuons.at (1)->eta);
88	>	}
89	>	if (cuts_->at ("maxEta"))
90	>	{
91	>	hists1D_["backMuonPT"]->Fill (leadMuons.at (0)->pt);
92	>	hists1D_["backMuonPT"]->Fill (leadMuons.at (1)->pt);
93	>	}
94	>	if (cuts_->at ("minPT") && cuts_->at ("maxEta"))
95	>	{
96	>	hists1D_["backMuonPhi"]->Fill (leadMuons.at (0)->phi);
97	>	hists1D_["backMuonPhi"]->Fill (leadMuons.at (1)->phi);
98	>	hists1D_["backDimuonMass"]->Fill (dimuonMass (leadMuons.at (0), leadMuons.at (1)));
99	>	}
100		}
101		}
102		}
103
104		void
105	<	OSUAnalysis::selectGoodMuons (const edm::Handle<BNmuonCollection> &muons, vector<BNmuonCollection::const_iterator> &goodMuons)
82	<	{
83	<	for (BNmuonCollection::const_iterator muon = muons->begin (); muon != muons->end (); muon++)
84	<	{
85	<	if (fabs (muon->eta) > maxEta_)
86	<	continue;
87	<	if (muon->pt < minPt_)
88	<	continue;
89	<	goodMuons.push_back (muon);
90	<	}
91	<	}
92	<
93	<	double
94	<	OSUAnalysis::dimuonMass (const BNmuonCollection::const_iterator &muon1, const BNmuonCollection::const_iterator &muon2)
95	<	{
96	<	double dimuonE, dimuonPx, dimuonPy, dimuonPz;
97	<
98	<	dimuonE = muon1->energy + muon2->energy;
99	<	dimuonPx = muon1->px + muon2->px;
100	<	dimuonPy = muon1->py + muon2->py;
101	<	dimuonPz = muon1->pz + muon2->pz;
102	<
103	<	return sqrt (dimuonE * dimuonE - dimuonPx * dimuonPx - dimuonPy * dimuonPy - dimuonPz * dimuonPz);
104	<	}
105	<
106	<	// Everything past this point can be broken off into a separate library.
107	<
108	<	void
109	<	OSUAnalysis::fillCutFlow ()
105	>	OSUAnalysis::selectLeadMuons (const edm::Handle<BNmuonCollection> &muons, vector<BNmuonCollection::const_iterator> &leadMuons)
106		{
107	<
112	<	// Fills the cut flow histograms from the contents of cuts_. Three histograms
113	<	// are filled: cutFlow shows the number of events passing each cut, applied
114	<	// cumulatively; selection shows the number of events passing each cut,
115	<	// applied independently; and complementarySelection shows the number of
116	<	// events passing all but one cut.
117	<
118	<	bool fillCumulative = true, fillComplement = true;
119	<	double complement = -1.0;
120	<
121	<	hists1D_["cutFlow"]->Fill (0.5);
122	<	hists1D_["selection"]->Fill (0.5);
123	<	for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
107	>	for (BNmuonCollection::const_iterator newMuon = muons->begin (); newMuon != muons->end (); newMuon++)
108		{
109	<	if (cuts_[*cut])
109	>	vector<BNmuonCollection::const_iterator>::iterator muon;
110	>	for (muon = leadMuons.begin (); muon != leadMuons.end (); muon++)
111		{
112	<	double binCenter = hists1D_["selection"]->GetBinCenter (cut - cutNames_.begin () + 2);
113	<
129	<	hists1D_["selection"]->Fill (binCenter);
130	<	if (fillCumulative)
131	<	hists1D_["cutFlow"]->Fill (binCenter);
132	<	}
133	<	else
134	<	{
135	<	fillCumulative = false;
136	<	if (complement < 0.0)
137	<	complement = hists1D_["complementarySelection"]->GetBinCenter (cut - cutNames_.begin () + 2);
138	<	else
139	<	fillComplement = false;
140	<	}
141	<	}
142	<	if (fillCumulative)
143	<	{
144	<	hists1D_["complementarySelection"]->Fill (0.5);
145	<	for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
146	<	{
147	<	double binCenter = hists1D_["complementarySelection"]->GetBinCenter (cut - cutNames_.begin () + 2);
148	<
149	<	hists1D_["complementarySelection"]->Fill (binCenter);
112	>	if ((*muon)->pt < newMuon->pt)
113	>	break;
114		}
115	+	leadMuons.insert (muon, newMuon);
116		}
117	<	if (!fillCumulative && fillComplement)
118	<	hists1D_["complementarySelection"]->Fill (complement);
154	<	}
155	<
156	<	void
157	<	OSUAnalysis::outputTime ()
158	<	{
159	<	// Outputs the run time of the job, both CPU time and real time.
160	<
161	<	double cpu, real;
162	<	int days, hours, minutes;
163	<
164	<	cpu = sw_.CpuTime ();
165	<	real = sw_.RealTime ();
166	<
167	<	days = (int) (cpu / (60.0 * 60.0 * 24.0));
168	<	cpu -= days * (60.0 * 60.0 * 24.0);
169	<	hours = (int) (cpu / (60.0 * 60.0));
170	<	cpu -= hours * (60.0 * 60.0);
171	<	minutes = (int) (cpu / 60.0);
172	<	cpu -= minutes * 60.0;
173	<
174	<	cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
175	<	cout << "CPU Time:  ";
176	<	if (days)
177	<	cout << days << " days, ";
178	<	if (days || hours)
179	<	cout << hours << " hours, ";
180	<	if (days || hours || minutes)
181	<	cout << minutes << " minutes, ";
182	<	if (days || hours || minutes || cpu)
183	<	cout << cpu << " seconds" << endl;
184	<
185	<	days = (int) (real / (60.0 * 60.0 * 24.0));
186	<	real -= days * (60.0 * 60.0 * 24.0);
187	<	hours = (int) (real / (60.0 * 60.0));
188	<	real -= hours * (60.0 * 60.0);
189	<	minutes = (int) (real / 60.0);
190	<	real -= minutes * 60.0;
191	<
192	<	cout << "Real Time: ";
193	<	if (days)
194	<	cout << days << " days, ";
195	<	if (days || hours)
196	<	cout << hours << " hours, ";
197	<	if (days || hours || minutes)
198	<	cout << minutes << " minutes, ";
199	<	if (days || hours || minutes || cpu)
200	<	cout << real << " seconds" << endl;
201	<	cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
117	>	if (leadMuons.size () > 2)
118	>	leadMuons.resize (2);
119		}
120
121	<	void
122	<	OSUAnalysis::outputCutFlow ()
121	>	double
122	>	OSUAnalysis::dimuonMass (const BNmuonCollection::const_iterator &muon1, const BNmuonCollection::const_iterator &muon2)
123		{
124	<	// Outputs cut flow numbers from the histograms written to the output file.
125	<
209	<	int totalEvents, totalEventsComplement;
210	<
211	<	cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
212	<	cout << setw (25) << left << "Cut Name" << right << setw (25) << "Cut Flow" << setw (25) << "Efficiency" << endl;
213	<	cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
214	<	totalEvents = hists1D_["cutFlow"]->GetBinContent (1);
215	<	cout << setw (25) << left << "Total Events:" << right << setw (25) << totalEvents << setw (25) << "100%" << endl;
216	<	for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
217	<	{
218	<	double cutFlow = hists1D_["cutFlow"]->GetBinContent (cut - cutNames_.begin () + 2);
219	<
220	<	cout << setw (25) << left << (*cut + ":") << right << setw (25) << cutFlow << setw (24) << 100.0 * (cutFlow / (double) totalEvents) << "%" << endl;
221	<	}
222	<	cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
124	>	TLorentzVector muVec1 (muon1->px, muon1->py, muon1->pz, muon1->energy),
125	>	muVec2 (muon2->px, muon2->py, muon2->pz, muon2->energy);
126
127	<	cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
225	<	cout << setw (25) << left << "Cut Name" << right << setw (25) << "Selection" << setw (25) << "Efficiency" << endl;
226	<	cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
227	<	cout << setw (25) << left << "Total Events:" << right << setw (25) << totalEvents << setw (25) << "100%" << endl;
228	<	for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
229	<	{
230	<	double selection  = hists1D_["selection"]->GetBinContent (cut - cutNames_.begin () + 2);
231	<
232	<	cout << setw (25) << left << (*cut + ":") << right << setw (25) << selection << setw (24) << 100.0 * (selection / (double) totalEvents) << "%" << endl;
233	<	}
234	<	cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
235	<
236	<	cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
237	<	cout << setw (25) << left << "Cut Name" << right << setw (25) << "Complementary Selection" << setw (25) << "Efficiency" << endl;
238	<	cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
239	<	totalEventsComplement = hists1D_["complementarySelection"]->GetBinContent (1);
240	<	cout << setw (25) << left << "Total Events:" << right << setw (25) << totalEventsComplement << setw (24) << 100.0 * (totalEventsComplement / (double) totalEvents) << "%" << endl;
241	<	for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
242	<	{
243	<	double complement = hists1D_["complementarySelection"]->GetBinContent (cut - cutNames_.begin () + 2);
244	<
245	<	cout << setw (25) << left << (*cut + ":") << right << setw (25) << complement << setw (24) << 100.0 * (complement / (double) totalEvents) << "%" << endl;
246	<	}
247	<	cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
127	>	return (muVec1 + muVec2).M ();
128		}
129
130		DEFINE_FWK_MODULE(OSUAnalysis);

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines