VHbb/interface/VHbbNameSpace.h

#include "TLorentzVector.h"
#include "TVector3.h"
#include "TVector2.h"
#include "TMath.h"
/*#if !defined(__CINT__) && !defined(__MAKECINT__)
#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
#endif*/

namespace VHbb {
  

  double deltaPhi(double phi1,double phi2)
  {
    double result = phi1 - phi2;
    while (result > TMath::Pi()) result -= 2*TMath::Pi();
    while (result <= -TMath::Pi()) result += 2*TMath::Pi();
    return result;
  }

  inline double deltaR(double eta1,double phi1,double eta2,double phi2)
    {
      double deta = eta1 - eta2;
      double dphi = deltaPhi(phi1, phi2);
      return TMath::Sqrt(deta*deta + dphi*dphi);
    }


  double Hmass( double V_eta,double V_phi,double V_pt, 
                double hJet1_eta,double hJet1_phi,double hJet1_pt, 
                double hJet2_eta,double hJet2_phi,double hJet2_pt ){
    
    TVector3 V(1,1,1);
    V.SetPtEtaPhi(V_pt,V_eta,V_phi);
    
    TVector3 H1(1,1,1);
    H1.SetPtEtaPhi(hJet1_pt,hJet1_eta,hJet1_phi);
    H1.SetMag(1/sin(H1.Theta()));
    
    TVector3 H2(1,1,1);
    H2.SetPtEtaPhi(hJet2_pt,hJet2_eta,hJet2_phi);
    H2.SetMag(1/sin(H2.Theta()));
    
    TVector3 n1(H1);
    TVector3 n2(H2);
    
    double det= n1.Px() * n2.Py() - n2.Px() * n1.Py();
    
    H1.SetMag( (  - n2.Py() * V.Px() + n2.Px() * V.Py() )  / (sin(n1.Theta()) *det ) );
    H2.SetMag( ( + n1.Py() * V.Px() - n1.Px() * V.Py() )  / (sin(n2.Theta())  *det ) );
    
    double mass=TMath::Sqrt( TMath::Power( (H1.Mag()+H2.Mag()),2 ) - TMath::Power(( ( H1+H2 ).Mag()),2) );
    
    return mass;
    
  }
  
  double Hmass_comb(double hJet1_eta,double hJet1_phi,double hJet1_pt, double hJet1_mass,
                    double hJet2_eta,double hJet2_phi,double hJet2_pt, double hJet2_mass){

    TLorentzVector H1, H2;
    H1.SetPtEtaPhiM(hJet1_pt,hJet1_eta,hJet1_phi, hJet1_mass);;
    H2.SetPtEtaPhiM(hJet2_pt,hJet2_eta,hJet2_phi, hJet2_mass);

    return (H1 + H2).M();

  }

  double Hmass_3j(double h_eta,double h_phi,double h_pt, double h_mass,
                  double aJet_eta,double aJet_phi,double aJet_pt, double aJet_mass){

    TLorentzVector H, H3;
    H.SetPtEtaPhiM( h_pt,h_eta,h_phi, h_mass);;
    H3.SetPtEtaPhiM(aJet_pt,aJet_eta,aJet_phi, aJet_mass);

    return (H + H3).M();


  }
  
  double ANGLELZ(double pt, double eta, double phi, double mass, double pt2, double eta2, double phi2, double mass2){
  TLorentzVector m1, m2, msum;
  m1.SetPtEtaPhiM(pt, eta, phi, mass);
  m2.SetPtEtaPhiM(pt2, eta2, phi2, mass2);
  msum = m1 + m2;

  TVector3 bZ =  msum.BoostVector();

  m1.Boost(-bZ);  
  m2.Boost(-bZ);  

  TVector3 b1;


  if((int) (pt) % 2 == 0)
    b1 =  m1.BoostVector();
  else
    b1 =  m2.BoostVector();

 double cosTheta = b1.Dot(msum.BoostVector()) / (b1.Mag()*msum.BoostVector().Mag());
 return(cosTheta);
  }


   double ANGLEHB(double pt, double eta, double phi, double e, double pt2, double eta2, double phi2, double e2){
 TLorentzVector m1, m2, msum;
 m1.SetPtEtaPhiE(pt, eta, phi, e);
 m2.SetPtEtaPhiE(pt2, eta2, phi2, e2);
 msum = m1 + m2;

 TVector3 bZ =  msum.BoostVector();

 m1.Boost(-bZ);
 m2.Boost(-bZ);  

 TVector3 b1;

  if((int) (pt) % 2 == 0)
    b1 =  m1.BoostVector();
  else
    b1 =  m2.BoostVector();

 double cosTheta = b1.Dot(msum.BoostVector()) / (b1.Mag()*msum.BoostVector().Mag());
 return(cosTheta);
   }

   double metCorSysShift(double met, double metphi, int Nvtx, int EVENT_run)
{
    double metx = met * cos(metphi);
    double mety = met * sin(metphi);
    double px = 0.0, py = 0.0;
    if (EVENT_run!=1) {
        //pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_data
        px = +1.68804e-01 + 3.37139e-01*Nvtx;
        py = -1.72555e-01 - 1.79594e-01*Nvtx;
    } else {
        //pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_mc
        px = +2.22335e-02 - 6.59183e-02*Nvtx;
        py = +1.52720e-01 - 1.28052e-01*Nvtx;
    }
    metx -= px;
    mety -= py;
    return std::sqrt(metx*metx + mety*mety);
}

    double metphiCorSysShift(double met, double metphi, int Nvtx, int EVENT_run)
{
    double metx = met * cos(metphi);
    double mety = met * sin(metphi);
    double px = 0.0, py = 0.0;
    if (EVENT_run!=1) {

        //pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_data
        px = +1.68804e-01 + 3.37139e-01*Nvtx;
        py = -1.72555e-01 - 1.79594e-01*Nvtx;
    } else {
        //pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_mc
        px = +2.22335e-02 - 6.59183e-02*Nvtx;
        py = +1.52720e-01 - 1.28052e-01*Nvtx;
    }
    metx -= px;
    mety -= py;
    if (metx == 0.0 && mety == 0.0)
        return 0.0;

    double phi1 = std::atan2(mety,metx);
    double phi2 = std::atan2(mety,metx)-2.0*M_PI;
    if (std::abs(phi1-metphi) < std::abs(phi2-metphi)+0.5*M_PI)
        return phi1;
    else
        return phi2;
}

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)
{
    double metx = met * cos(metphi);
    double mety = met * sin(metphi);
    TLorentzVector j1;
    TLorentzVector j2;
    j1.SetPtEtaPhiE(pt1,eta1,phi1, e1 );
    j2.SetPtEtaPhiE(pt2,eta2,phi2, e2 );
    metx += j1.Px()*(1-corr1);
    metx += j2.Px()*(1-corr2);
    mety += j1.Py()*(1-corr1);
    mety += j2.Py()*(1-corr2);
    TVector2 corrMET(metx, mety);
     return corrMET;
}

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){
    return metType1Reg(met, metphi, corr1, corr2, pt1, eta1, phi1, e1, pt2, eta2, phi2, e2).Phi();

}
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){
    return metType1Reg(met, metphi, corr1, corr2, pt1, eta1, phi1, e1, pt2, eta2, phi2, e2).Mod();

}


   double met_MPF(double met, double metphi, double pt, double phi)
{
    return 1.+met*pt*std::cos( deltaPhi(metphi,phi) ) / (pt*pt);

}

double resolutionBias(double eta)
{
// return 0;//Nominal!
  if(eta< 0.5) return 0.052;
  if(eta< 1.1) return 0.057;
  if(eta< 1.7) return 0.096;
  if(eta< 2.3) return 0.134;
  if(eta< 5) return 0.28;
  return 0;
}

double evalJERBias( double ptreco, double ptgen, double eta1){
  double eta = fabs(eta1);
  double cor =1;   
  if ((fabs(ptreco - ptgen)/ ptreco)<0.5) { //Limit the effect to the core 
     cor = (ptreco +resolutionBias(eta) *(ptreco-ptgen))/ptreco;   
  }
  if (ptgen > 0.) return ptreco*cor;
  else return ptreco;
}

double evalEt( double pt, double eta, double phi, double e){
  TLorentzVector j;
  j.SetPtEtaPhiE(pt,eta,phi, e );
  return j.Et(); 

}

double evalMt( double pt, double eta, double phi, double e){
  TLorentzVector j;
  j.SetPtEtaPhiE(pt,eta,phi, e );
  return j.Mt(); 

}
/*double evalJECUnc( double pt, double eta){
// Total uncertainty for reference
JetCorrectionUncertainty *total = new JetCorrectionUncertainty("/shome/nmohr/CMSSW_5_2_6_patch1/src/UserCode/VHbb/data/START53_V15MC_Uncertainty_AK5PFchs.txt");

total->setJetPt(pt);
total->setJetEta(eta);
double uncert =  total->getUncertainty(true);
delete total;
return uncert;
}*/


}

Revision:	1.6
Committed:	Thu Mar 14 14:21:36 2013 UTC (12 years, 2 months ago) by nmohr
Content type:	text/plain
Branch:	MAIN
Changes since 1.5:	+58 -5 lines
Log Message:	Add new resolution
#	Content
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	{
15	double result = phi1 - phi2;
16	while (result > TMath::Pi()) result -= 2*TMath::Pi();
17	while (result <= -TMath::Pi()) result += 2*TMath::Pi();
18	return result;
19	}
20
21	inline double deltaR(double eta1,double phi1,double eta2,double phi2)
22	{
23	double deta = eta1 - eta2;
24	double dphi = deltaPhi(phi1, phi2);
25	return TMath::Sqrt(detadeta + dphidphi);
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(metxmetx + metymety);
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.+metptstd::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
252
253	}
254