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;
  }*/

  double ptWeightDY( double lheV_pt, double sign = 1.)
  {
    double SF = 1.;
    if (50. < lheV_pt && lheV_pt < 100.){
      SF = 0.873885+0.00175853*sign*lheV_pt;
    }
    else if (lheV_pt > 100){
      SF = 1.10651-0.000705265*sign*lheV_pt;
    }
    return SF;
  }

  // weights correction for EWK NLO correction
  double ptWeightZllH( double genHPt){
    double SF = 1.;
    if ( genHPt > 1.) SF = 0.999757-0.000174527*genHPt;
    return SF>0?SF:0;
  }


  double minCSVold(int EVENT_run, double hJet_csv0, double hJet_csv1, double hJet_csvOld0, double hJet_csvOld1){

    if( EVENT_run < 2 )
      return std::min(hJet_csvOld0,hJet_csvOld1);
    else
      return std::min(hJet_csv0,hJet_csv1);
  }

  double maxCSVold(int EVENT_run, double hJet_csv0, double hJet_csv1, double hJet_csvOld0, double hJet_csvOld1){

    if( EVENT_run < 2 )
      return std::max(hJet_csvOld0,hJet_csvOld1);
    else
      return std::max(hJet_csv0,hJet_csv1);
  }


  // weights correction for EWK NLO correction from Atlas
  double ewkAtlas8TeVZllH( double genHPt, double genZPt){
    if (genHPt < 1.) return 1.;
    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};
    int corrBin = (int) ( (genZPt - 25.) / 5. );
    if (corrBin < 0) corrBin = 0;
    if (corrBin > 94) corrBin = 94;
    double SF = 1.+hll8_contents[corrBin];
    return SF;
  }
  
  double mueEff(int Vtype, double eta0, double eta1, double pt0, double pt1){
      if (Vtype == 0) return 1.;
      if (Vtype == 1) {
            double corr = 1.0;
            if (abs(eta0) < 1.3) corr*=0.95;
            if (abs(eta1) < 1.3) corr*=0.95;
            if (min(pt0,pt1) < 25.) corr*=0.6;
            else if (min(pt0,pt1) < 30.) corr*=0.8;
            // scale to lumi ABC/D
            corr = 1.-(1.-corr)*6.9/12.1;
            return corr;
      }
      return 1.;
  }

}

Revision:	1.13
Committed:	Mon Apr 15 11:45:39 2013 UTC (12 years ago) by peller
Content type:	text/plain
Branch:	MAIN
Changes since 1.12:	+12 -3 lines
Log Message:	correction for electrons
#	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	bortigno	1.8	#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	bortigno	1.8	TLorentzVector m1, m2, msum;
83			m1.SetPtEtaPhiM(pt, eta, phi, mass);
84			m2.SetPtEtaPhiM(pt2, eta2, phi2, mass2);
85			msum = m1 + m2;
86	nmohr	1.4
87	bortigno	1.8	TVector3 bZ = msum.BoostVector();
88	nmohr	1.4
89	bortigno	1.8	m1.Boost(-bZ);
90			m2.Boost(-bZ);
91	nmohr	1.4
92	bortigno	1.8	TVector3 b1;
93	nmohr	1.4
94
95	bortigno	1.8	if((int) (pt) % 2 == 0)
96			b1 = m1.BoostVector();
97			else
98			b1 = m2.BoostVector();
99	nmohr	1.4
100	bortigno	1.8	double cosTheta = b1.Dot(msum.BoostVector()) / (b1.Mag()*msum.BoostVector().Mag());
101			return(cosTheta);
102	nmohr	1.4	}
103
104
105	bortigno	1.8	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	nmohr	1.4
111	bortigno	1.8	TVector3 bZ = msum.BoostVector();
112	nmohr	1.4
113	bortigno	1.8	m1.Boost(-bZ);
114			m2.Boost(-bZ);
115	nmohr	1.4
116	bortigno	1.8	TVector3 b1;
117	nmohr	1.4
118	bortigno	1.8	if((int) (pt) % 2 == 0)
119			b1 = m1.BoostVector();
120			else
121			b1 = m2.BoostVector();
122	nmohr	1.4
123	bortigno	1.8	double cosTheta = b1.Dot(msum.BoostVector()) / (b1.Mag()*msum.BoostVector().Mag());
124			return(cosTheta);
125			}
126	nmohr	1.4
127	bortigno	1.8	double metCorSysShift(double met, double metphi, int Nvtx, int EVENT_run)
128			{
129	nmohr	1.5	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	bortigno	1.8	//pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_data
134			px = +1.68804e-01 + 3.37139e-01*Nvtx;
135			py = -1.72555e-01 - 1.79594e-01*Nvtx;
136	nmohr	1.5	} else {
137	bortigno	1.8	//pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_mc
138			px = +2.22335e-02 - 6.59183e-02*Nvtx;
139			py = +1.52720e-01 - 1.28052e-01*Nvtx;
140	nmohr	1.5	}
141			metx -= px;
142			mety -= py;
143			return std::sqrt(metxmetx + metymety);
144	bortigno	1.8	}
145	nmohr	1.5
146	bortigno	1.8	double metphiCorSysShift(double met, double metphi, int Nvtx, int EVENT_run)
147			{
148	nmohr	1.5	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	bortigno	1.8	//pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_data
154			px = +1.68804e-01 + 3.37139e-01*Nvtx;
155			py = -1.72555e-01 - 1.79594e-01*Nvtx;
156	nmohr	1.5	} else {
157	bortigno	1.8	//pfMEtSysShiftCorrParameters_2012runAplusBvsNvtx_mc
158			px = +2.22335e-02 - 6.59183e-02*Nvtx;
159			py = +1.52720e-01 - 1.28052e-01*Nvtx;
160	nmohr	1.5	}
161			metx -= px;
162			mety -= py;
163			if (metx == 0.0 && mety == 0.0)
164	bortigno	1.8	return 0.0;
165	nmohr	1.5
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	bortigno	1.8	return phi1;
170	nmohr	1.5	else
171	bortigno	1.8	return phi2;
172			}
173	bortigno	1.2
174	bortigno	1.8	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	nmohr	1.6	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	bortigno	1.8	return corrMET;
188			}
189	nmohr	1.6
190	bortigno	1.8	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	nmohr	1.6	return metType1Reg(met, metphi, corr1, corr2, pt1, eta1, phi1, e1, pt2, eta2, phi2, e2).Phi();
192
193	bortigno	1.8	}
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	nmohr	1.6	return metType1Reg(met, metphi, corr1, corr2, pt1, eta1, phi1, e1, pt2, eta2, phi2, e2).Mod();
196
197	bortigno	1.8	}
198	nmohr	1.6
199
200	bortigno	1.8	double met_MPF(double met, double metphi, double pt, double phi)
201			{
202	nmohr	1.6	return 1.+metptstd::cos( deltaPhi(metphi,phi) ) / (pt*pt);
203
204	bortigno	1.8	}
205	nmohr	1.6
206	bortigno	1.8	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	nmohr	1.5
217	bortigno	1.8	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	nmohr	1.5	}
226
227	bortigno	1.8	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	nmohr	1.5
232	bortigno	1.8	}
233	nmohr	1.5
234	bortigno	1.8	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	nmohr	1.5
239	bortigno	1.8	}
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	nmohr	1.6
251	nmohr	1.9	double ptWeightDY( double lheV_pt, double sign = 1.)
252	bortigno	1.8	{
253	nmohr	1.7	double SF = 1.;
254			if (50. < lheV_pt && lheV_pt < 100.){
255	nmohr	1.9	SF = 0.873885+0.00175853signlheV_pt;
256	nmohr	1.7	}
257			else if (lheV_pt > 100){
258	nmohr	1.9	SF = 1.10651-0.000705265signlheV_pt;
259	nmohr	1.7	}
260			return SF;
261	bortigno	1.8	}
262	nmohr	1.6
263	bortigno	1.8	// weights correction for EWK NLO correction
264	nmohr	1.9	double ptWeightZllH( double genHPt){
265	bortigno	1.8	double SF = 1.;
266	nmohr	1.9	if ( genHPt > 1.) SF = 0.999757-0.000174527*genHPt;
267	bortigno	1.8	return SF>0?SF:0;
268			}
269	bortigno	1.2
270	bortigno	1.11
271	bortigno	1.10	double minCSVold(int EVENT_run, double hJet_csv0, double hJet_csv1, double hJet_csvOld0, double hJet_csvOld1){
272
273			if( EVENT_run < 2 )
274	nmohr	1.12	return std::min(hJet_csvOld0,hJet_csvOld1);
275	bortigno	1.10	else
276	nmohr	1.12	return std::min(hJet_csv0,hJet_csv1);
277	bortigno	1.10	}
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	nmohr	1.12	return std::max(hJet_csvOld0,hJet_csvOld1);
283	bortigno	1.10	else
284	nmohr	1.12	return std::max(hJet_csv0,hJet_csv1);
285	bortigno	1.10	}
286
287
288	nmohr	1.9	// weights correction for EWK NLO correction from Atlas
289			double ewkAtlas8TeVZllH( double genHPt, double genZPt){
290	bortigno	1.10	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	nmohr	1.9	}
298	nmohr	1.12
299	peller	1.13	double mueEff(int Vtype, double eta0, double eta1, double pt0, double pt1){
300			if (Vtype == 0) return 1.;
301			if (Vtype == 1) {
302			double corr = 1.0;
303			if (abs(eta0) < 1.3) corr*=0.95;
304			if (abs(eta1) < 1.3) corr*=0.95;
305			if (min(pt0,pt1) < 25.) corr*=0.6;
306			else if (min(pt0,pt1) < 30.) corr*=0.8;
307			// scale to lumi ABC/D
308			corr = 1.-(1.-corr)*6.9/12.1;
309			return corr;
310			}
311	nmohr	1.12	return 1.;
312			}
313	nmohr	1.9
314	bortigno	1.1	}
315