AnaTools/src/SFWeight.cc

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

MuonSFWeight::MuonSFWeight (const string &sfFile, const string &dataOverMC)
{

  TFile *fin = TFile::Open (sfFile.c_str ());
  TGraphAsymmErrors *data = (TGraphAsymmErrors *) fin->Get (dataOverMC.c_str ());

  double bins[data->GetN () + 1];
  for (int i = 0; i < data->GetN (); i++)
    {
      double x, y, exl, exh;

      data->GetPoint (i, x, y);
      exl = data->GetErrorXlow (i);
      exh = data->GetErrorXhigh (i);
      if (i == 0)
        bins[i] = x - exl;
      bins[i + 1] = x + exh;
    }

  muonSFWeight_ = new TH1D ("muonSFWeight", "", data->GetN (), bins);
  muonSFWeight_->SetDirectory (0);
  for (int i = 0; i < data->GetN (); i++)
    {
      double dataX, dataY;

      data->GetPoint (i, dataX, dataY);
      muonSFWeight_->SetBinContent (i + 1, dataY);
    }

  fin->Close ();
  delete data;
}

MuonSFWeight::MuonSFWeight (const string &sfFile, const string &dataEff, const string &mcEff)
{
  TFile *fin = TFile::Open (sfFile.c_str ());
  TGraphAsymmErrors *mc = (TGraphAsymmErrors *) fin->Get (mcEff.c_str ());
  TGraphAsymmErrors *data = (TGraphAsymmErrors *) fin->Get (dataEff.c_str ());

  double bins[mc->GetN () + 1];
  for (int i = 0; i < mc->GetN (); i++)
    {
      double x, y, exl, exh;

      mc->GetPoint (i, x, y);
      exl = mc->GetErrorXlow (i);
      exh = mc->GetErrorXhigh (i);
      if (i == 0)
        bins[i] = x - exl;
      bins[i + 1] = x + exh;
    }

  muonSFWeight_ = new TH1D ("muonSFWeight", "", mc->GetN (), bins);
  muonSFWeight_->SetDirectory (0);
  for (int i = 0; i < mc->GetN (); i++)
    {
      double mcX, mcY, dataX, dataY;

      mc->GetPoint (i, mcX, mcY);
      data->GetPoint (i, dataX, dataY);
      muonSFWeight_->SetBinContent (i + 1, dataY / mcY);
    }

  fin->Close ();
  delete mc;
  delete data;
}

MuonSFWeight::~MuonSFWeight ()
{
  delete muonSFWeight_;
}

ElectronSFWeight::ElectronSFWeight (const string &cmsswRelease, const string &id) :
  cmsswRelease_ (cmsswRelease),
  id_ (id)
{
}

double
ElectronSFWeight::at (const double &eta, const double &pt)
{
  //cout << eta << " " << pt << endl;
  //cout << "doing it\n";
  if (cmsswRelease_ == "53X")
    {
      if (id_ == "loose")
        {
          if (fabs (eta) < 0.8)
            {
              if (pt < 15.0 && pt >= 10.0)
                return 0.855;
              if (pt < 20.0 && pt >= 15.0)
                return 0.962;
              if (pt < 30.0 && pt >= 20.0)
                return 1.005;
              if (pt < 40.0 && pt >= 30.0)
                return 1.004;
              if (pt < 50.0 && pt >= 40.0)
                return 1.008;
              if (pt >= 50.0)
                return 1.008;
            }
          if (fabs (eta) < 1.442 && fabs (eta) >= 0.8)
            {
              if (pt < 15.0 && pt >= 10.0)
                return 0.858;
              if (pt < 20.0 && pt >= 15.0)
                return 0.962;
              if (pt < 30.0 && pt >= 20.0)
                return 0.981;
              if (pt < 40.0 && pt >= 30.0)
                return 0.991;
              if (pt < 50.0 && pt >= 40.0)
                return 0.994;
              if (pt >= 50.0)
                return 0.999;
            }
          if (fabs (eta) < 1.556 && fabs (eta) >= 1.442)
            {
              if (pt < 15.0 && pt >= 10.0)
                return 1.109;
              if (pt < 20.0 && pt >= 15.0)
                return 0.903;
              if (pt < 30.0 && pt >= 20.0)
                return 1.044;
              if (pt < 40.0 && pt >= 30.0)
                return 0.998;
              if (pt < 50.0 && pt >= 40.0)
                return 0.989;
              if (pt >= 50.0)
                return 0.994;
            }
          if (fabs (eta) < 2.0 && fabs (eta) >= 1.556)
            {
              if (pt < 15.0 && pt >= 10.0)
                return 0.838;
              if (pt < 20.0 && pt >= 15.0)
                return 0.939;
              if (pt < 30.0 && pt >= 20.0)
                return 0.980;
              if (pt < 40.0 && pt >= 30.0)
                return 0.992;
              if (pt < 50.0 && pt >= 40.0)
                return 1.004;
              if (pt >= 50.0)
                return 1.006;
            }
          if (fabs (eta) < 2.5 && fabs (eta) >= 2.0)
            {
              if (pt < 15.0 && pt >= 10.0)
                return 1.034;
              if (pt < 20.0 && pt >= 15.0)
                return 0.970;
              if (pt < 30.0 && pt >= 20.0)
                return 1.017;
              if (pt < 40.0 && pt >= 30.0)
                return 1.019;
              if (pt < 50.0 && pt >= 40.0)
                return 1.005;
              if (pt >= 50.0)
                return 1.009;
            }
        }
      if (id_ == "tight")
        {
          if (fabs (eta) < 0.8)
            {
              if (pt < 15.0 && pt >= 10.0)
                return 0.818;
              if (pt < 20.0 && pt >= 15.0)
                return 0.928;
              if (pt < 30.0 && pt >= 20.0)
                return 0.973;
              if (pt < 40.0 && pt >= 30.0)
                return 0.979;
              if (pt < 50.0 && pt >= 40.0)
                return 0.984;
              if (pt >= 50.0)
                return 0.983;
            }
          if (fabs (eta) < 1.442 && fabs (eta) >= 0.8)
            {
              if (pt < 15.0 && pt >= 10.0)
                return 0.840;
              if (pt < 20.0 && pt >= 15.0)
                return 0.914;
              if (pt < 30.0 && pt >= 20.0)
                return 0.948;
              if (pt < 40.0 && pt >= 30.0)
                return 0.961;
              if (pt < 50.0 && pt >= 40.0)
                return 0.972;
              if (pt >= 50.0)
                return 0.977;
            }
          if (fabs (eta) < 1.556 && fabs (eta) >= 1.442)
            {
              if (pt < 15.0 && pt >= 10.0)
                return 1.008;
              if (pt < 20.0 && pt >= 15.0)
                return 0.877;
              if (pt < 30.0 && pt >= 20.0)
                return 0.983;
              if (pt < 40.0 && pt >= 30.0)
                return 0.983;
              if (pt < 50.0 && pt >= 40.0)
                return 0.957;
              if (pt >= 50.0)
                return 0.978;
            }
          if (fabs (eta) < 2.0 && fabs (eta) >= 1.556)
            {
              if (pt < 15.0 && pt >= 10.0)
                return 0.906;
              if (pt < 20.0 && pt >= 15.0)
                return 0.907;
              if (pt < 30.0 && pt >= 20.0)
                return 0.957;
              if (pt < 40.0 && pt >= 30.0)
                return 0.962;
              if (pt < 50.0 && pt >= 40.0)
                return 0.985;
              if (pt >= 50.0)
                return 0.986;
            }
          if (fabs (eta) < 2.5 && fabs (eta) >= 2.0)
            {
              if (pt < 15.0 && pt >= 10.0)
                return 0.991;
              if (pt < 20.0 && pt >= 15.0)
                return 0.939;
              if (pt < 30.0 && pt >= 20.0)
                return 1.001;
              if (pt < 40.0 && pt >= 30.0)
                return 1.002;
              if (pt < 50.0 && pt >= 40.0)
                return 0.999;
              if (pt >= 50.0)
                return 0.995;
            }
        }
      if (id_ == "mvaNonTrig")
        {
          //cout << "mvaNonTrig\n";
          if (fabs (eta) < 0.8)
            {
              if (pt < 20.0)
                return 1.0;
              if (pt < 30.0 && pt >= 20.0)
                return 0.994;
              if (pt < 40.0 && pt >= 30.0)
                return 0.974;
              if (pt < 50.0 && pt >= 40.0)
                return 0.988;
              if (pt >= 50.0)
                return 0.999;
            }
          if (fabs (eta) < 1.478 && fabs (eta) >= 0.8)
            {
              if (pt < 20.0)
                return 1.0;
              if (pt < 30.0 && pt >= 20.0)
                return 0.971;
              if (pt < 40.0 && pt >= 30.0)
                return 0.967;
              if (pt < 50.0 && pt >= 40.0)
                return 0.984;
              if (pt >= 50.0)
                return 0.998;
            }
          if (fabs (eta) <= 2.5 && fabs (eta) >= 1.478)
            {
              if (pt < 20.0)
                return 1.0;
              if (pt < 30.0 && pt >= 20.0)
                return 0.921;
              if (pt < 40.0 && pt >= 30.0)
                return 0.932;
              if (pt < 50.0 && pt >= 40.0)
                return 0.968;
              if (pt >= 50.0)
                return 0.983;
            }
        }
    }

  return 0.0;
}

ElectronSFWeight::~ElectronSFWeight ()
{
}
Revision:	1.6
Committed:	Mon Aug 19 15:51:14 2013 UTC (11 years, 8 months ago) by lantonel
Content type:	text/plain
Branch:	MAIN
CVS Tags:	HEAD
Changes since 1.5:	+13 -13 lines
Error occurred while calculating annotation data.
Log Message:	don't veto electrons above 200 GeV
#	Content
1	#include "OSUT3Analysis/AnaTools/interface/SFWeight.h"
2
3	MuonSFWeight::MuonSFWeight (const string &sfFile, const string &dataOverMC)
4	{
5
6	TFile *fin = TFile::Open (sfFile.c_str ());
7	TGraphAsymmErrors data = (TGraphAsymmErrors ) fin->Get (dataOverMC.c_str ());
8
9	double bins[data->GetN () + 1];
10	for (int i = 0; i < data->GetN (); i++)
11	{
12	double x, y, exl, exh;
13
14	data->GetPoint (i, x, y);
15	exl = data->GetErrorXlow (i);
16	exh = data->GetErrorXhigh (i);
17	if (i == 0)
18	bins[i] = x - exl;
19	bins[i + 1] = x + exh;
20	}
21
22	muonSFWeight_ = new TH1D ("muonSFWeight", "", data->GetN (), bins);
23	muonSFWeight_->SetDirectory (0);
24	for (int i = 0; i < data->GetN (); i++)
25	{
26	double dataX, dataY;
27
28	data->GetPoint (i, dataX, dataY);
29	muonSFWeight_->SetBinContent (i + 1, dataY);
30	}
31
32	fin->Close ();
33	delete data;
34	}
35
36	MuonSFWeight::MuonSFWeight (const string &sfFile, const string &dataEff, const string &mcEff)
37	{
38	TFile *fin = TFile::Open (sfFile.c_str ());
39	TGraphAsymmErrors mc = (TGraphAsymmErrors ) fin->Get (mcEff.c_str ());
40	TGraphAsymmErrors data = (TGraphAsymmErrors ) fin->Get (dataEff.c_str ());
41
42	double bins[mc->GetN () + 1];
43	for (int i = 0; i < mc->GetN (); i++)
44	{
45	double x, y, exl, exh;
46
47	mc->GetPoint (i, x, y);
48	exl = mc->GetErrorXlow (i);
49	exh = mc->GetErrorXhigh (i);
50	if (i == 0)
51	bins[i] = x - exl;
52	bins[i + 1] = x + exh;
53	}
54
55	muonSFWeight_ = new TH1D ("muonSFWeight", "", mc->GetN (), bins);
56	muonSFWeight_->SetDirectory (0);
57	for (int i = 0; i < mc->GetN (); i++)
58	{
59	double mcX, mcY, dataX, dataY;
60
61	mc->GetPoint (i, mcX, mcY);
62	data->GetPoint (i, dataX, dataY);
63	muonSFWeight_->SetBinContent (i + 1, dataY / mcY);
64	}
65
66	fin->Close ();
67	delete mc;
68	delete data;
69	}
70
71	MuonSFWeight::~MuonSFWeight ()
72	{
73	delete muonSFWeight_;
74	}
75
76	ElectronSFWeight::ElectronSFWeight (const string &cmsswRelease, const string &id) :
77	cmsswRelease_ (cmsswRelease),
78	id_ (id)
79	{
80	}
81
82	double
83	ElectronSFWeight::at (const double &eta, const double &pt)
84	{
85	//cout << eta << " " << pt << endl;
86	//cout << "doing it\n";
87	if (cmsswRelease_ == "53X")
88	{
89	if (id_ == "loose")
90	{
91	if (fabs (eta) < 0.8)
92	{
93	if (pt < 15.0 && pt >= 10.0)
94	return 0.855;
95	if (pt < 20.0 && pt >= 15.0)
96	return 0.962;
97	if (pt < 30.0 && pt >= 20.0)
98	return 1.005;
99	if (pt < 40.0 && pt >= 30.0)
100	return 1.004;
101	if (pt < 50.0 && pt >= 40.0)
102	return 1.008;
103	if (pt >= 50.0)
104	return 1.008;
105	}
106	if (fabs (eta) < 1.442 && fabs (eta) >= 0.8)
107	{
108	if (pt < 15.0 && pt >= 10.0)
109	return 0.858;
110	if (pt < 20.0 && pt >= 15.0)
111	return 0.962;
112	if (pt < 30.0 && pt >= 20.0)
113	return 0.981;
114	if (pt < 40.0 && pt >= 30.0)
115	return 0.991;
116	if (pt < 50.0 && pt >= 40.0)
117	return 0.994;
118	if (pt >= 50.0)
119	return 0.999;
120	}
121	if (fabs (eta) < 1.556 && fabs (eta) >= 1.442)
122	{
123	if (pt < 15.0 && pt >= 10.0)
124	return 1.109;
125	if (pt < 20.0 && pt >= 15.0)
126	return 0.903;
127	if (pt < 30.0 && pt >= 20.0)
128	return 1.044;
129	if (pt < 40.0 && pt >= 30.0)
130	return 0.998;
131	if (pt < 50.0 && pt >= 40.0)
132	return 0.989;
133	if (pt >= 50.0)
134	return 0.994;
135	}
136	if (fabs (eta) < 2.0 && fabs (eta) >= 1.556)
137	{
138	if (pt < 15.0 && pt >= 10.0)
139	return 0.838;
140	if (pt < 20.0 && pt >= 15.0)
141	return 0.939;
142	if (pt < 30.0 && pt >= 20.0)
143	return 0.980;
144	if (pt < 40.0 && pt >= 30.0)
145	return 0.992;
146	if (pt < 50.0 && pt >= 40.0)
147	return 1.004;
148	if (pt >= 50.0)
149	return 1.006;
150	}
151	if (fabs (eta) < 2.5 && fabs (eta) >= 2.0)
152	{
153	if (pt < 15.0 && pt >= 10.0)
154	return 1.034;
155	if (pt < 20.0 && pt >= 15.0)
156	return 0.970;
157	if (pt < 30.0 && pt >= 20.0)
158	return 1.017;
159	if (pt < 40.0 && pt >= 30.0)
160	return 1.019;
161	if (pt < 50.0 && pt >= 40.0)
162	return 1.005;
163	if (pt >= 50.0)
164	return 1.009;
165	}
166	}
167	if (id_ == "tight")
168	{
169	if (fabs (eta) < 0.8)
170	{
171	if (pt < 15.0 && pt >= 10.0)
172	return 0.818;
173	if (pt < 20.0 && pt >= 15.0)
174	return 0.928;
175	if (pt < 30.0 && pt >= 20.0)
176	return 0.973;
177	if (pt < 40.0 && pt >= 30.0)
178	return 0.979;
179	if (pt < 50.0 && pt >= 40.0)
180	return 0.984;
181	if (pt >= 50.0)
182	return 0.983;
183	}
184	if (fabs (eta) < 1.442 && fabs (eta) >= 0.8)
185	{
186	if (pt < 15.0 && pt >= 10.0)
187	return 0.840;
188	if (pt < 20.0 && pt >= 15.0)
189	return 0.914;
190	if (pt < 30.0 && pt >= 20.0)
191	return 0.948;
192	if (pt < 40.0 && pt >= 30.0)
193	return 0.961;
194	if (pt < 50.0 && pt >= 40.0)
195	return 0.972;
196	if (pt >= 50.0)
197	return 0.977;
198	}
199	if (fabs (eta) < 1.556 && fabs (eta) >= 1.442)
200	{
201	if (pt < 15.0 && pt >= 10.0)
202	return 1.008;
203	if (pt < 20.0 && pt >= 15.0)
204	return 0.877;
205	if (pt < 30.0 && pt >= 20.0)
206	return 0.983;
207	if (pt < 40.0 && pt >= 30.0)
208	return 0.983;
209	if (pt < 50.0 && pt >= 40.0)
210	return 0.957;
211	if (pt >= 50.0)
212	return 0.978;
213	}
214	if (fabs (eta) < 2.0 && fabs (eta) >= 1.556)
215	{
216	if (pt < 15.0 && pt >= 10.0)
217	return 0.906;
218	if (pt < 20.0 && pt >= 15.0)
219	return 0.907;
220	if (pt < 30.0 && pt >= 20.0)
221	return 0.957;
222	if (pt < 40.0 && pt >= 30.0)
223	return 0.962;
224	if (pt < 50.0 && pt >= 40.0)
225	return 0.985;
226	if (pt >= 50.0)
227	return 0.986;
228	}
229	if (fabs (eta) < 2.5 && fabs (eta) >= 2.0)
230	{
231	if (pt < 15.0 && pt >= 10.0)
232	return 0.991;
233	if (pt < 20.0 && pt >= 15.0)
234	return 0.939;
235	if (pt < 30.0 && pt >= 20.0)
236	return 1.001;
237	if (pt < 40.0 && pt >= 30.0)
238	return 1.002;
239	if (pt < 50.0 && pt >= 40.0)
240	return 0.999;
241	if (pt >= 50.0)
242	return 0.995;
243	}
244	}
245	if (id_ == "mvaNonTrig")
246	{
247	//cout << "mvaNonTrig\n";
248	if (fabs (eta) < 0.8)
249	{
250	if (pt < 20.0)
251	return 1.0;
252	if (pt < 30.0 && pt >= 20.0)
253	return 0.994;
254	if (pt < 40.0 && pt >= 30.0)
255	return 0.974;
256	if (pt < 50.0 && pt >= 40.0)
257	return 0.988;
258	if (pt >= 50.0)
259	return 0.999;
260	}
261	if (fabs (eta) < 1.478 && fabs (eta) >= 0.8)
262	{
263	if (pt < 20.0)
264	return 1.0;
265	if (pt < 30.0 && pt >= 20.0)
266	return 0.971;
267	if (pt < 40.0 && pt >= 30.0)
268	return 0.967;
269	if (pt < 50.0 && pt >= 40.0)
270	return 0.984;
271	if (pt >= 50.0)
272	return 0.998;
273	}
274	if (fabs (eta) <= 2.5 && fabs (eta) >= 1.478)
275	{
276	if (pt < 20.0)
277	return 1.0;
278	if (pt < 30.0 && pt >= 20.0)
279	return 0.921;
280	if (pt < 40.0 && pt >= 30.0)
281	return 0.932;
282	if (pt < 50.0 && pt >= 40.0)
283	return 0.968;
284	if (pt >= 50.0)
285	return 0.983;
286	}
287	}
288	}
289
290	return 0.0;
291	}
292
293	ElectronSFWeight::~ElectronSFWeight ()
294	{
295	}