ViewVC Help

View File | Revision Log | Show Annotations | Root Listing

root/cvsroot/UserCode/VHbb/interface/VHbbNameSpace.h

Comparing UserCode/VHbb/interface/VHbbNameSpace.h (file contents):
Revision 1.1 by bortigno, Fri Dec 14 15:59:45 2012 UTC vs.
Revision 1.11 by bortigno, Fri Mar 29 14:42:55 2013 UTC

#	Line 1 | Line 1
1	<
1	>	#include "TLorentzVector.h"
2	>	#include "TVector3.h"
3	>	#include "TVector2.h"
4		#include "TMath.h"
5	+	/*#if !defined(__CINT__) && !defined(__MAKECINT__)
6	+	#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
7	+	#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
8	+	#endif*/
9
10		namespace VHbb {
11	+
12
13		double deltaPhi(double phi1,double phi2)
14		{
#	Line 15 | Line 22 | namespace VHbb {
22		{
23		double deta = eta1 - eta2;
24		double dphi = deltaPhi(phi1, phi2);
25	<	return std::sqrt(deta*deta + dphi*dphi);
25	>	return TMath::Sqrt(deta*deta + dphi*dphi);
26	>	}
27	>
28	>
29	>	double Hmass( double V_eta,double V_phi,double V_pt,
30	>	double hJet1_eta,double hJet1_phi,double hJet1_pt,
31	>	double hJet2_eta,double hJet2_phi,double hJet2_pt ){
32	>
33	>	TVector3 V(1,1,1);
34	>	V.SetPtEtaPhi(V_pt,V_eta,V_phi);
35	>
36	>	TVector3 H1(1,1,1);
37	>	H1.SetPtEtaPhi(hJet1_pt,hJet1_eta,hJet1_phi);
38	>	H1.SetMag(1/sin(H1.Theta()));
39	>
40	>	TVector3 H2(1,1,1);
41	>	H2.SetPtEtaPhi(hJet2_pt,hJet2_eta,hJet2_phi);
42	>	H2.SetMag(1/sin(H2.Theta()));
43	>
44	>	TVector3 n1(H1);
45	>	TVector3 n2(H2);
46	>
47	>	double det= n1.Px() * n2.Py() - n2.Px() * n1.Py();
48	>
49	>	H1.SetMag( (  - n2.Py() * V.Px() + n2.Px() * V.Py() )  / (sin(n1.Theta()) *det ) );
50	>	H2.SetMag( ( + n1.Py() * V.Px() - n1.Px() * V.Py() )  / (sin(n2.Theta())  *det ) );
51	>
52	>	double mass=TMath::Sqrt( TMath::Power( (H1.Mag()+H2.Mag()),2 ) - TMath::Power(( ( H1+H2 ).Mag()),2) );
53	>
54	>	return mass;
55	>
56	>	}
57	>
58	>	double Hmass_comb(double hJet1_eta,double hJet1_phi,double hJet1_pt, double hJet1_mass,
59	>	double hJet2_eta,double hJet2_phi,double hJet2_pt, double hJet2_mass){
60	>
61	>	TLorentzVector H1, H2;
62	>	H1.SetPtEtaPhiM(hJet1_pt,hJet1_eta,hJet1_phi, hJet1_mass);;
63	>	H2.SetPtEtaPhiM(hJet2_pt,hJet2_eta,hJet2_phi, hJet2_mass);
64	>
65	>	return (H1 + H2).M();
66	>
67	>	}
68	>
69	>	double Hmass_3j(double h_eta,double h_phi,double h_pt, double h_mass,
70	>	double aJet_eta,double aJet_phi,double aJet_pt, double aJet_mass){
71	>
72	>	TLorentzVector H, H3;
73	>	H.SetPtEtaPhiM( h_pt,h_eta,h_phi, h_mass);;
74	>	H3.SetPtEtaPhiM(aJet_pt,aJet_eta,aJet_phi, aJet_mass);
75	>
76	>	return (H + H3).M();
77	>
78	>
79	>	}
80	>
81	>	double ANGLELZ(double pt, double eta, double phi, double mass, double pt2, double eta2, double phi2, double mass2){
82	>	TLorentzVector m1, m2, msum;
83	>	m1.SetPtEtaPhiM(pt, eta, phi, mass);
84	>	m2.SetPtEtaPhiM(pt2, eta2, phi2, mass2);
85	>	msum = m1 + m2;
86	>
87	>	TVector3 bZ =  msum.BoostVector();
88	>
89	>	m1.Boost(-bZ);
90	>	m2.Boost(-bZ);
91	>
92	>	TVector3 b1;
93	>
94	>
95	>	if((int) (pt) % 2 == 0)
96	>	b1 =  m1.BoostVector();
97	>	else
98	>	b1 =  m2.BoostVector();
99	>
100	>	double cosTheta = b1.Dot(msum.BoostVector()) / (b1.Mag()*msum.BoostVector().Mag());
101	>	return(cosTheta);
102	>	}
103	>
104	>
105	>	double ANGLEHB(double pt, double eta, double phi, double e, double pt2, double eta2, double phi2, double e2){
106	>	TLorentzVector m1, m2, msum;
107	>	m1.SetPtEtaPhiE(pt, eta, phi, e);
108	>	m2.SetPtEtaPhiE(pt2, eta2, phi2, e2);
109	>	msum = m1 + m2;
110	>
111	>	TVector3 bZ =  msum.BoostVector();
112	>
113	>	m1.Boost(-bZ);
114	>	m2.Boost(-bZ);
115	>
116	>	TVector3 b1;
117	>
118	>	if((int) (pt) % 2 == 0)
119	>	b1 =  m1.BoostVector();
120	>	else
121	>	b1 =  m2.BoostVector();
122	>
123	>	double cosTheta = b1.Dot(msum.BoostVector()) / (b1.Mag()*msum.BoostVector().Mag());
124	>	return(cosTheta);
125	>	}
126	>
127	>	double metCorSysShift(double met, double metphi, int Nvtx, int EVENT_run)
128	>	{
129	>	double metx = met * cos(metphi);
130	>	double mety = met * sin(metphi);
131	>	double px = 0.0, py = 0.0;
132	>	if (EVENT_run!=1) {
133	>	//pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_data
134	>	px = +1.68804e-01 + 3.37139e-01*Nvtx;
135	>	py = -1.72555e-01 - 1.79594e-01*Nvtx;
136	>	} else {
137	>	//pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_mc
138	>	px = +2.22335e-02 - 6.59183e-02*Nvtx;
139	>	py = +1.52720e-01 - 1.28052e-01*Nvtx;
140		}
141	+	metx -= px;
142	+	mety -= py;
143	+	return std::sqrt(metx*metx + mety*mety);
144	+	}
145	+
146	+	double metphiCorSysShift(double met, double metphi, int Nvtx, int EVENT_run)
147	+	{
148	+	double metx = met * cos(metphi);
149	+	double mety = met * sin(metphi);
150	+	double px = 0.0, py = 0.0;
151	+	if (EVENT_run!=1) {
152	+
153	+	//pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_data
154	+	px = +1.68804e-01 + 3.37139e-01*Nvtx;
155	+	py = -1.72555e-01 - 1.79594e-01*Nvtx;
156	+	} else {
157	+	//pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_mc
158	+	px = +2.22335e-02 - 6.59183e-02*Nvtx;
159	+	py = +1.52720e-01 - 1.28052e-01*Nvtx;
160	+	}
161	+	metx -= px;
162	+	mety -= py;
163	+	if (metx == 0.0 && mety == 0.0)
164	+	return 0.0;
165	+
166	+	double phi1 = std::atan2(mety,metx);
167	+	double phi2 = std::atan2(mety,metx)-2.0*M_PI;
168	+	if (std::abs(phi1-metphi) < std::abs(phi2-metphi)+0.5*M_PI)
169	+	return phi1;
170	+	else
171	+	return phi2;
172	+	}
173	+
174	+	TVector2 metType1Reg(double met, double metphi, double corr1, double corr2, double pt1, double eta1, double phi1, double e1, double pt2, double eta2, double phi2, double e2)
175	+	{
176	+	double metx = met * cos(metphi);
177	+	double mety = met * sin(metphi);
178	+	TLorentzVector j1;
179	+	TLorentzVector j2;
180	+	j1.SetPtEtaPhiE(pt1,eta1,phi1, e1 );
181	+	j2.SetPtEtaPhiE(pt2,eta2,phi2, e2 );
182	+	metx += j1.Px()*(1-corr1);
183	+	metx += j2.Px()*(1-corr2);
184	+	mety += j1.Py()*(1-corr1);
185	+	mety += j2.Py()*(1-corr2);
186	+	TVector2 corrMET(metx, mety);
187	+	return corrMET;
188	+	}
189	+
190	+	double metType1Phi(double met, double metphi, double corr1, double corr2, double pt1, double eta1, double phi1, double e1, double pt2, double eta2, double phi2, double e2){
191	+	return metType1Reg(met, metphi, corr1, corr2, pt1, eta1, phi1, e1, pt2, eta2, phi2, e2).Phi();
192	+
193	+	}
194	+	double metType1Et(double met, double metphi, double corr1, double corr2, double pt1, double eta1, double phi1, double e1, double pt2, double eta2, double phi2, double e2){
195	+	return metType1Reg(met, metphi, corr1, corr2, pt1, eta1, phi1, e1, pt2, eta2, phi2, e2).Mod();
196	+
197	+	}
198	+
199	+
200	+	double met_MPF(double met, double metphi, double pt, double phi)
201	+	{
202	+	return 1.+met*pt*std::cos( deltaPhi(metphi,phi) ) / (pt*pt);
203	+
204	+	}
205	+
206	+	double resolutionBias(double eta)
207	+	{
208	+	// return 0;//Nominal!
209	+	if(eta< 0.5) return 0.052;
210	+	if(eta< 1.1) return 0.057;
211	+	if(eta< 1.7) return 0.096;
212	+	if(eta< 2.3) return 0.134;
213	+	if(eta< 5) return 0.28;
214	+	return 0;
215	+	}
216	+
217	+	double evalJERBias( double ptreco, double ptgen, double eta1){
218	+	double eta = fabs(eta1);
219	+	double cor =1;
220	+	if ((fabs(ptreco - ptgen)/ ptreco)<0.5) { //Limit the effect to the core
221	+	cor = (ptreco +resolutionBias(eta) *(ptreco-ptgen))/ptreco;
222	+	}
223	+	if (ptgen > 0.) return ptreco*cor;
224	+	else return ptreco;
225	+	}
226	+
227	+	double evalEt( double pt, double eta, double phi, double e){
228	+	TLorentzVector j;
229	+	j.SetPtEtaPhiE(pt,eta,phi, e );
230	+	return j.Et();
231	+
232	+	}
233	+
234	+	double evalMt( double pt, double eta, double phi, double e){
235	+	TLorentzVector j;
236	+	j.SetPtEtaPhiE(pt,eta,phi, e );
237	+	return j.Mt();
238	+
239	+	}
240	+	/*double evalJECUnc( double pt, double eta){
241	+	// Total uncertainty for reference
242	+	JetCorrectionUncertainty *total = new JetCorrectionUncertainty("/shome/nmohr/CMSSW_5_2_6_patch1/src/UserCode/VHbb/data/START53_V15MC_Uncertainty_AK5PFchs.txt");
243	+
244	+	total->setJetPt(pt);
245	+	total->setJetEta(eta);
246	+	double uncert =  total->getUncertainty(true);
247	+	delete total;
248	+	return uncert;
249	+	}*/
250	+
251	+	double ptWeightDY( double lheV_pt, double sign = 1.)
252	+	{
253	+	double SF = 1.;
254	+	if (50. < lheV_pt && lheV_pt < 100.){
255	+	SF = 0.873885+0.00175853*sign*lheV_pt;
256	+	}
257	+	else if (lheV_pt > 100){
258	+	SF = 1.10651-0.000705265*sign*lheV_pt;
259	+	}
260	+	return SF;
261	+	}
262	+
263	+	// weights correction for EWK NLO correction
264	+	double ptWeightZllH( double genHPt){
265	+	double SF = 1.;
266	+	if ( genHPt > 1.) SF = 0.999757-0.000174527*genHPt;
267	+	return SF>0?SF:0;
268	+	}
269	+
270	+
271	+	double minCSVold(int EVENT_run, double hJet_csv0, double hJet_csv1, double hJet_csvOld0, double hJet_csvOld1){
272	+
273	+	if( EVENT_run < 2 )
274	+	return min(hJet_csvOld0,hJet_csvOld1);
275	+	else
276	+	return min(hJet_csv0,hJet_csv1);
277	+	}
278	+
279	+	double maxCSVold(int EVENT_run, double hJet_csv0, double hJet_csv1, double hJet_csvOld0, double hJet_csvOld1){
280	+
281	+	if( EVENT_run < 2 )
282	+	return max(hJet_csvOld0,hJet_csvOld1);
283	+	else
284	+	return max(hJet_csv0,hJet_csv1);
285	+	}
286	+
287	+
288	+	// weights correction for EWK NLO correction from Atlas
289	+	double ewkAtlas8TeVZllH( double genHPt, double genZPt){
290	+	if (genHPt < 1.) return 1.;
291	+	double hll8_contents[95] = {0.000664024, -0.00357095, -0.00767076, -0.00967366, -0.0134844, -0.0157148, -0.0181885, -0.0209647, -0.0232788, -0.0252373, -0.0265634, -0.0275069, -0.0285776, -0.0281683, -0.0294206, -0.0299975, -0.0308047, -0.0311716, -0.030913, -0.0324821, -0.0323192, -0.0324639, -0.0319356, -0.0322621, -0.0331146, -0.0338905, -0.0345189, -0.0358591, -0.0358407, -0.040018, -0.0396389, -0.0407177, -0.0445103, -0.0441406, -0.0471215, -0.0463301, -0.0513777, -0.0536773, -0.0546446, -0.0568508, -0.0590333, -0.0612157, -0.0633981, -0.0655805, -0.067763, -0.0699454, -0.0721278, -0.0743103, -0.0764927, -0.0786751, -0.0808575, -0.08304, -0.0852224, -0.0874048, -0.0895872, -0.0917697, -0.0939521, -0.0961345, -0.098317, -0.100499, -0.102682, -0.104864, -0.107047, -0.109229, -0.111412, -0.113594, -0.115776, -0.117959, -0.120141, -0.122324, -0.124506, -0.126689, -0.128871, -0.131053, -0.133236, -0.135418, -0.137601, -0.139783, -0.141965, -0.144148, -0.14633, -0.148513, -0.150695, -0.152878, -0.15506, -0.157242, -0.159425, -0.161607, -0.16379, -0.165972, -0.168155, -0.170337, -0.172519, -0.174702, -0.176884};
292	+	int corrBin = (int) ( (genZPt - 25.) / 5. );
293	+	if (corrBin < 0) corrBin = 0;
294	+	if (corrBin > 94) corrBin = 94;
295	+	double SF = 1.+hll8_contents[corrBin];
296	+	return SF;
297	+	}
298	+
299		}
300

Diff Legend

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