EventMix/src/sPHElectronPairv1.cc

0001 #include <TLorentzVector.h>
0002
0003 #include "sPHElectronPairv1.h"
0004
0005 sPHElectronPairv1::sPHElectronPairv1()
0006 {
0007   _id = -99999;
0008   _type = 0;
0009   _e1 = sPHElectronv1();
0010   _e2 = sPHElectronv1();
0011 }
0012
0013 sPHElectronPairv1::sPHElectronPairv1(sPHElectronv1* e1, sPHElectronv1* e2)
0014 {
0015   _id = e1->get_id() + e2->get_id()*100000;
0016   _type = 0;
0017   _e1 = *e1;
0018   _e2 = *e2;
0019 }
0020
0021 double sPHElectronPairv1::get_mass() const
0022 {
0023   double px1 = _e1.get_px();
0024   double py1 = _e1.get_py();
0025   double pz1 = _e1.get_pz();
0026   double ee1 = sqrt(px1*px1+py1*py1+pz1*pz1+0.000511*0.000511);
0027   TLorentzVector tmp1 = TLorentzVector(px1,py1,pz1,ee1);
0028   double px2 = _e2.get_px();
0029   double py2 = _e2.get_py();
0030   double pz2 = _e2.get_pz();
0031   double ee2 = sqrt(px2*px2+py2*py2+pz2*pz2+0.000511*0.000511);
0032   TLorentzVector tmp2 = TLorentzVector(px2,py2,pz2,ee2);
0033   TLorentzVector pair = tmp1 + tmp2;
0034   return pair.M();
0035 }
0036
0037 double sPHElectronPairv1::get_pt() const
0038 {
0039   double px1 = _e1.get_px();
0040   double py1 = _e1.get_py();
0041   double pz1 = _e1.get_pz();
0042   double ee1 = sqrt(px1*px1+py1*py1+pz1*pz1+0.000511*0.000511);
0043   TLorentzVector tmp1 = TLorentzVector(px1,py1,pz1,ee1);
0044   double px2 = _e2.get_px();
0045   double py2 = _e2.get_py();
0046   double pz2 = _e2.get_pz();
0047   double ee2 = sqrt(px2*px2+py2*py2+pz2*pz2+0.000511*0.000511);
0048   TLorentzVector tmp2 = TLorentzVector(px2,py2,pz2,ee2);
0049   TLorentzVector pair = tmp1 + tmp2;
0050   return pair.Pt();
0051 }
0052
0053 double sPHElectronPairv1::get_eta() const
0054 {
0055   double px1 = _e1.get_px();
0056   double py1 = _e1.get_py();
0057   double pz1 = _e1.get_pz();
0058   double ee1 = sqrt(px1*px1+py1*py1+pz1*pz1+0.000511*0.000511);
0059   TLorentzVector tmp1 = TLorentzVector(px1,py1,pz1,ee1);
0060   double px2 = _e2.get_px();
0061   double py2 = _e2.get_py();
0062   double pz2 = _e2.get_pz();
0063   double ee2 = sqrt(px2*px2+py2*py2+pz2*pz2+0.000511*0.000511);
0064   TLorentzVector tmp2 = TLorentzVector(px2,py2,pz2,ee2);
0065   TLorentzVector pair = tmp1 + tmp2;
0066   return pair.Eta();
0067 }
0068
0069 double sPHElectronPairv1::get_phiv() const
0070 {
0071 // assumes that the first element of the pair is the positron, and second the electron
0072   double px1 = _e1.get_px();
0073   double py1 = _e1.get_py();
0074   double pz1 = _e1.get_pz();
0075   double px2 = _e2.get_px();
0076   double py2 = _e2.get_py();
0077   double pz2 = _e2.get_pz();
0078   double px = px1 + px2;
0079   double py = py1 + py2;
0080   double pz = pz1 + pz2;
0081   double pp = sqrt(px*px+py*py+pz*pz);
0082   double ux = px/pp;
0083   double uy = py/pp;
0084   double uz = pz/pp;
0085
0086 // axis defined by (ux,uy,ux)X(0,0,1).
0087 // this is the axis that is perpendicular to the direction of
0088 // pair, and also perpendicular to the Z axis (field direction).
0089 // If the pair is conversion at R!=0, it must have (apparent)
0090 // momentum component in this axis (caused by field intergral from the
0091 // vertex point to the conversion point).
0092 // The sign of the component is opposite for e+ and e-.
0093 //
0094 // (ux,uy,ux)X(0,0,1)=(uy,-ux,0)
0095 //
0096
0097   double ax =  uy/sqrt(ux*ux+uy*uy);
0098   double ay = -ux/sqrt(ux*ux+uy*uy);
0099
0100 //momentum of e+ and e- in (ax,ay,az) axis. Note that az=0 by definition.
0101
0102 //vector product of pep X pem
0103   double vpx  = py1*pz2 - pz1*py2;
0104   double vpy  = pz1*px2 - px1*pz2;
0105   double vpz  = px1*py2 - py1*px2;
0106   double vp   = sqrt(vpx*vpx+vpy*vpy+vpz*vpz);
0107
0108 //unit vector of pep X pem
0109   double vx = vpx/vp;
0110   double vy = vpy/vp;
0111   double vz = vpz/vp;
0112
0113 //The third axis defined by vector product (ux,uy,uz)X(vx,vy,vz)
0114   double wx = uy*vz - uz*vy;
0115   double wy = uz*vx - ux*vz;
0116 //  double wz = ux*vy - uy*vx;
0117 //  double wl = sqrt(wx*wx+wy*wy+wz*wz);
0118 //// by construction, (wx,wy,wz) must be a unit vector.
0119 //  if(fabs(wl - 1.0) > 0.00001) std:;cout << "Calculation error in W vector"<<std::endl;
0120
0121 // measure angle between (wx,wy,wz) and (ax,ay,0). The angle between them
0122 // should be small if the pair is conversion
0123   double cosPhiV = wx*ax + wy*ay;
0124   double phiV = acos(cosPhiV);
0125
0126   return phiV;
0127 }
0128