AnaTools/src/OSUAnalysis.cc

#include "OSUT3Analysis/AnaTools/interface/OSUAnalysis.h"

OSUAnalysis::OSUAnalysis (const edm::ParameterSet &cfg) :
  // Retrieve parameters from the configuration file
  muons_ (cfg.getParameter<edm::InputTag> ("muons")),
  maxEta_ (cfg.getParameter<double> ("maxEta")),
  minPt_ (cfg.getParameter<double> ("minPt"))
{
  // Register cut names
  cutNames_.push_back ("minMuons");

  // Create additional histograms
  hists1D_["muonPT"] = fs_->make<TH1D> ("muonPt", ";p_{T} (GeV/c)", 100, 0.0, 300.0);
  hists1D_["muonEta"] = fs_->make<TH1D> ("muonEta", ";#eta", 100, -3.0, 3.0);
  hists1D_["muonPhi"] = fs_->make<TH1D> ("muonPhi", ";#phi", 100, -5.0, 5.0);
  hists1D_["dimuonMass"] = fs_->make<TH1D> ("dimuonMass", ";m_{#mu^{+} #mu^{-}} (GeV/c^{2})", 90, 30.0, 120.0);

  TH1::SetDefaultSumw2 ();
  hists1D_["cutFlow"] = fs_->make<TH1D> ("cutFlow", "", cutNames_.size () + 1, 0.0, cutNames_.size () + 1);
  hists1D_["selection"] = fs_->make<TH1D> ("selection", "", cutNames_.size () + 1, 0.0, cutNames_.size () + 1);
  hists1D_["complementarySelection"] = fs_->make<TH1D> ("complementarySelection", "", cutNames_.size () + 1, 0.0, cutNames_.size () + 1);

  hists1D_["cutFlow"]->GetXaxis ()->SetBinLabel (1, "Total Events");
  hists1D_["selection"]->GetXaxis ()->SetBinLabel (1, "Total Events");
  hists1D_["complementarySelection"]->GetXaxis ()->SetBinLabel (1, "Total Events");
  for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
    {
      hists1D_["cutFlow"]->GetXaxis ()->SetBinLabel (cut - cutNames_.begin () + 2, cut->c_str ());
      hists1D_["selection"]->GetXaxis ()->SetBinLabel (cut - cutNames_.begin () + 2, cut->c_str ());
      hists1D_["complementarySelection"]->GetXaxis ()->SetBinLabel (cut - cutNames_.begin () + 2, cut->c_str ());
    }
}

void
OSUAnalysis::beginJob ()
{
  sw_.Start ();
}

void
OSUAnalysis::endJob ()
{
  sw_.Stop ();
  outputTime ();
  outputCutFlow ();
}

void
OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
{
  // Retrieve necessary collections from the event
  edm::Handle<BNmuonCollection> muons;
  event.getByLabel (muons_, muons);

  // Perform selection on objects within the event
  vector<BNmuonCollection::const_iterator> goodMuons;
  selectGoodMuons (muons, goodMuons);

  // Perform selection calculations and store the results in the map cuts_. The
  // keys to cuts_ must match the strings in cutNames_
  cuts_["minMuons"] = goodMuons.size () >= 2;

  fillCutFlow ();

  // Fill histograms while applying desired cuts.
  for (vector<BNmuonCollection::const_iterator>::const_iterator muon1 = goodMuons.begin (); muon1 != goodMuons.end (); muon1++)
    {
      hists1D_["muonPT"]->Fill ((*muon1)->pt);
      hists1D_["muonEta"]->Fill ((*muon1)->eta);
      hists1D_["muonPhi"]->Fill ((*muon1)->phi);
      for (vector<BNmuonCollection::const_iterator>::const_iterator muon2 = muon1 + 1; muon2 != goodMuons.end (); muon2++)
        {
          if ((*muon1)->charge * (*muon2)->charge > 0)
            continue;
          hists1D_["dimuonMass"]->Fill (dimuonMass (*muon1, *muon2));
        }
    }
}

void
OSUAnalysis::selectGoodMuons (const edm::Handle<BNmuonCollection> &muons, vector<BNmuonCollection::const_iterator> &goodMuons)
{
  for (BNmuonCollection::const_iterator muon = muons->begin (); muon != muons->end (); muon++)
    {
      if (fabs (muon->eta) > maxEta_)
        continue;
      if (muon->pt < minPt_)
        continue;
      goodMuons.push_back (muon);
    }
}

double
OSUAnalysis::dimuonMass (const BNmuonCollection::const_iterator &muon1, const BNmuonCollection::const_iterator &muon2)
{
  double dimuonE, dimuonPx, dimuonPy, dimuonPz;

  dimuonE = muon1->energy + muon2->energy;
  dimuonPx = muon1->px + muon2->px;
  dimuonPy = muon1->py + muon2->py;
  dimuonPz = muon1->pz + muon2->pz;

  return sqrt (dimuonE * dimuonE - dimuonPx * dimuonPx - dimuonPy * dimuonPy - dimuonPz * dimuonPz);
}

// Everything past this point can be broken off into a separate library.

void
OSUAnalysis::fillCutFlow ()
{

  // Fills the cut flow histograms from the contents of cuts_. Three histograms
  // are filled: cutFlow shows the number of events passing each cut, applied
  // cumulatively; selection shows the number of events passing each cut,
  // applied independently; and complementarySelection shows the number of
  // events passing all but one cut.

  bool fillCumulative = true, fillComplement = true;
  double complement = -1.0;

  hists1D_["cutFlow"]->Fill (0.5);
  hists1D_["selection"]->Fill (0.5);
  for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
    {
      if (cuts_[*cut])
        {
          double binCenter = hists1D_["selection"]->GetBinCenter (cut - cutNames_.begin () + 2);

          hists1D_["selection"]->Fill (binCenter);
          if (fillCumulative)
            hists1D_["cutFlow"]->Fill (binCenter);
        }
      else
        {
          fillCumulative = false;
          if (complement < 0.0)
            complement = hists1D_["complementarySelection"]->GetBinCenter (cut - cutNames_.begin () + 2);
          else
            fillComplement = false;
        }
    }
  if (fillCumulative)
    {
      hists1D_["complementarySelection"]->Fill (0.5);
      for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
        {
          double binCenter = hists1D_["complementarySelection"]->GetBinCenter (cut - cutNames_.begin () + 2);

          hists1D_["complementarySelection"]->Fill (binCenter);
        }
    }
  if (!fillCumulative && fillComplement)
    hists1D_["complementarySelection"]->Fill (complement);
}

void
OSUAnalysis::outputTime ()
{
  // Outputs the run time of the job, both CPU time and real time.

  double cpu, real;
  int days, hours, minutes;

  cpu = sw_.CpuTime ();
  real = sw_.RealTime ();

  days = (int) (cpu / (60.0 * 60.0 * 24.0));
  cpu -= days * (60.0 * 60.0 * 24.0);
  hours = (int) (cpu / (60.0 * 60.0));
  cpu -= hours * (60.0 * 60.0);
  minutes = (int) (cpu / 60.0);
  cpu -= minutes * 60.0;

  cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
  cout << "CPU Time:  ";
  if (days)
    cout << days << " days, ";
  if (days || hours)
    cout << hours << " hours, ";
  if (days || hours || minutes)
    cout << minutes << " minutes, ";
  if (days || hours || minutes || cpu)
    cout << cpu << " seconds" << endl;

  days = (int) (real / (60.0 * 60.0 * 24.0));
  real -= days * (60.0 * 60.0 * 24.0);
  hours = (int) (real / (60.0 * 60.0));
  real -= hours * (60.0 * 60.0);
  minutes = (int) (real / 60.0);
  real -= minutes * 60.0;

  cout << "Real Time: ";
  if (days)
    cout << days << " days, ";
  if (days || hours)
    cout << hours << " hours, ";
  if (days || hours || minutes)
    cout << minutes << " minutes, ";
  if (days || hours || minutes || cpu)
    cout << real << " seconds" << endl;
  cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
}

void
OSUAnalysis::outputCutFlow ()
{
  // Outputs cut flow numbers from the histograms written to the output file.

  int totalEvents, totalEventsComplement;

  cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
  cout << setw (25) << left << "Cut Name" << right << setw (25) << "Cut Flow" << setw (25) << "Efficiency" << endl;
  cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
  totalEvents = hists1D_["cutFlow"]->GetBinContent (1);
  cout << setw (25) << left << "Total Events:" << right << setw (25) << totalEvents << setw (25) << "100%" << endl;
  for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
    {
      double cutFlow = hists1D_["cutFlow"]->GetBinContent (cut - cutNames_.begin () + 2);

      cout << setw (25) << left << (*cut + ":") << right << setw (25) << cutFlow << setw (24) << 100.0 * (cutFlow / (double) totalEvents) << "%" << endl;
    }
  cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;

  cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
  cout << setw (25) << left << "Cut Name" << right << setw (25) << "Selection" << setw (25) << "Efficiency" << endl;
  cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
  cout << setw (25) << left << "Total Events:" << right << setw (25) << totalEvents << setw (25) << "100%" << endl;
  for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
    {
      double selection  = hists1D_["selection"]->GetBinContent (cut - cutNames_.begin () + 2);

      cout << setw (25) << left << (*cut + ":") << right << setw (25) << selection << setw (24) << 100.0 * (selection / (double) totalEvents) << "%" << endl;
    }
  cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;

  cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
  cout << setw (25) << left << "Cut Name" << right << setw (25) << "Complementary Selection" << setw (25) << "Efficiency" << endl;
  cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
  totalEventsComplement = hists1D_["complementarySelection"]->GetBinContent (1);
  cout << setw (25) << left << "Total Events:" << right << setw (25) << totalEventsComplement << setw (24) << 100.0 * (totalEventsComplement / (double) totalEvents) << "%" << endl;
  for (vector<string>::const_iterator cut = cutNames_.begin (); cut != cutNames_.end (); cut++)
    {
      double complement = hists1D_["complementarySelection"]->GetBinContent (cut - cutNames_.begin () + 2);

      cout << setw (25) << left << (*cut + ":") << right << setw (25) << complement << setw (24) << 100.0 * (complement / (double) totalEvents) << "%" << endl;
    }
  cout << setw (75) << setfill ('-') << '-' << setfill (' ') << endl;
}

DEFINE_FWK_MODULE(OSUAnalysis);
Revision:	1.1
Committed:	Thu Jun 14 17:42:31 2012 UTC (12 years, 11 months ago) by ahart
Content type:	text/plain
Branch:	MAIN
Log Message:	Allow fuzzy dataset names for osudb and osusub. First commit of AnaTools directory.
#	User	Rev	Content
1	ahart	1.1	#include "OSUT3Analysis/AnaTools/interface/OSUAnalysis.h"
2
3			OSUAnalysis::OSUAnalysis (const edm::ParameterSet &cfg) :
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");
11
12			// Create additional histograms
13			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);
15			hists1D_["muonPhi"] = fs_->make<TH1D> ("muonPhi", ";#phi", 100, -5.0, 5.0);
16			hists1D_["dimuonMass"] = fs_->make<TH1D> ("dimuonMass", ";m_{#mu^{+} #mu^{-}} (GeV/c^{2})", 90, 30.0, 120.0);
17
18			TH1::SetDefaultSumw2 ();
19			hists1D_["cutFlow"] = fs_->make<TH1D> ("cutFlow", "", cutNames_.size () + 1, 0.0, cutNames_.size () + 1);
20			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 ();
38			}
39
40			void
41			OSUAnalysis::endJob ()
42			{
43			sw_.Stop ();
44			outputTime ();
45			outputCutFlow ();
46			}
47
48			void
49			OSUAnalysis::analyze (const edm::Event &event, const edm::EventSetup &setup)
50			{
51			// Retrieve necessary collections from the event
52			edm::Handle<BNmuonCollection> muons;
53			event.getByLabel (muons_, muons);
54
55			// Perform selection on objects within the event
56			vector<BNmuonCollection::const_iterator> goodMuons;
57			selectGoodMuons (muons, goodMuons);
58
59			// Perform selection calculations and store the results in the map cuts_. The
60			// keys to cuts_ must match the strings in cutNames_
61			cuts_["minMuons"] = goodMuons.size () >= 2;
62
63			fillCutFlow ();
64
65			// Fill histograms while applying desired cuts.
66			for (vector<BNmuonCollection::const_iterator>::const_iterator muon1 = goodMuons.begin (); muon1 != goodMuons.end (); muon1++)
67			{
68			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++)
72			{
73			if ((muon1)->charge (*muon2)->charge > 0)
74			continue;
75			hists1D_["dimuonMass"]->Fill (dimuonMass (muon1, muon2));
76			}
77			}
78			}
79
80			void
81			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 ()
110			{
111
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++)
124			{
125			if (cuts_[*cut])
126			{
127			double binCenter = hists1D_["selection"]->GetBinCenter (cut - cutNames_.begin () + 2);
128
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);
150			}
151			}
152			if (!fillCumulative && fillComplement)
153			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;
202			}
203
204			void
205			OSUAnalysis::outputCutFlow ()
206			{
207			// Outputs cut flow numbers from the histograms written to the output file.
208
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;
223
224			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;
248			}
249
250			DEFINE_FWK_MODULE(OSUAnalysis);