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
#	User	Rev	Content
1	bortigno	1.2	#include "TLorentzVector.h"
2			#include "TVector3.h"
3	nmohr	1.6	#include "TVector2.h"
4	bortigno	1.1	#include "TMath.h"
5	nmohr	1.6	/*#if !defined(__CINT__) && !defined(__MAKECINT__)
6			#include "CondFormats/JetMETObjects/interface/JetCorrectorParameters.h"
7			#include "CondFormats/JetMETObjects/interface/JetCorrectionUncertainty.h"
8			#endif*/
9	bortigno	1.1
10			namespace VHbb {
11	nmohr	1.6
12	bortigno	1.1
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	bortigno	1.3	return TMath::Sqrt(detadeta + dphidphi);
26	bortigno	1.1	}
27	bortigno	1.2
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	nmohr	1.5	double det= n1.Px() * n2.Py() - n2.Px() * n1.Py();
48	bortigno	1.2
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	nmohr	1.5	double mass=TMath::Sqrt( TMath::Power( (H1.Mag()+H2.Mag()),2 ) - TMath::Power(( ( H1+H2 ).Mag()),2) );
53	bortigno	1.2
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	nmohr	1.4
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	nmohr	1.5	double cosTheta = b1.Dot(msum.BoostVector()) / (b1.Mag()*msum.BoostVector().Mag());
101	nmohr	1.4	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	nmohr	1.5	double cosTheta = b1.Dot(msum.BoostVector()) / (b1.Mag()*msum.BoostVector().Mag());
124	nmohr	1.4	return(cosTheta);
125			}
126
127	nmohr	1.5	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	bortigno	1.2
174	nmohr	1.6	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	nmohr	1.5	double resolutionBias(double eta)
207			{
208			// return 0;//Nominal!
209	nmohr	1.6	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	nmohr	1.5	}
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	nmohr	1.6	if (ptgen > 0.) return ptreco*cor;
224			else return ptreco;
225	nmohr	1.5	}
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	nmohr	1.6	/*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	bortigno	1.2
253	bortigno	1.1	}
254