TagAndProbe/interface/ColinsSoperVariables.h

0001 #ifndef PhysicsTools_TagAndProbe_ColinsSuperVariables
0002 #define PhysicsTools_TagAndProbe_ColinsSuperVariables
0003
0004 #define CM_ENERGY 7000.0
0005 #include "TLorentzVector.h"
0006 #include "TVector3.h"
0007
0008 // calculate the Colins-Soper variables;
0009 // everything is in the lab frame
0010 inline void calCSVariables(TLorentzVector mu, TLorentzVector mubar, double *res, bool swap) {
0011   // convention. beam direction is on the positive Z direction.
0012   // beam contains quark flux.
0013   TLorentzVector Pbeam(0, 0, CM_ENERGY / 2.0, CM_ENERGY / 2.0);
0014   TLorentzVector Ptarget(0, 0, -CM_ENERGY / 2.0, CM_ENERGY / 2.0);
0015
0016   TLorentzVector Q(mu + mubar);
0017   /************************************************************************
0018    *
0019    * 1) cos(theta) = 2 Q^-1 (Q^2+Qt^2)^-1 (mu^+ mubar^- - mu^- mubar^+)
0020    *
0021    *
0022    ************************************************************************/
0023   double muplus = 1.0 / sqrt(2.0) * (mu.E() + mu.Z());
0024   double muminus = 1.0 / sqrt(2.0) * (mu.E() - mu.Z());
0025
0026   double mubarplus = 1.0 / sqrt(2.0) * (mubar.E() + mubar.Z());
0027   double mubarminus = 1.0 / sqrt(2.0) * (mubar.E() - mubar.Z());
0028
0029   double costheta =
0030       2.0 / Q.Mag() / sqrt(pow(Q.Mag(), 2) + pow(Q.Pt(), 2)) * (muplus * mubarminus - muminus * mubarplus);
0031   if (swap)
0032     costheta = -costheta;
0033
0034   /************************************************************************
0035    *
0036    * 2) sin2(theta) = Q^-2 Dt^2 - Q^-2 (Q^2 + Qt^2)^-1 * (Dt dot Qt)^2
0037    *
0038    ************************************************************************/
0039   TLorentzVector D(mu - mubar);
0040   double dt_qt = D.X() * Q.X() + D.Y() * Q.Y();
0041   double sin2theta =
0042       pow(D.Pt() / Q.Mag(), 2) - 1.0 / pow(Q.Mag(), 2) / (pow(Q.Mag(), 2) + pow(Q.Pt(), 2)) * pow(dt_qt, 2);
0043
0044   /************************************************************************
0045    *
0046    * 3) tanphi = (Q^2 + Qt^2)^1/2 / Q (Dt dot R unit) /(Dt dot Qt unit)
0047    *
0048    ************************************************************************/
0049   // unit vector on R direction
0050   TVector3 R = Pbeam.Vect().Cross(Q.Vect());
0051   TVector3 Runit = R.Unit();
0052
0053   // unit vector on Qt
0054   TVector3 Qt = Q.Vect();
0055   Qt.SetZ(0);
0056   TVector3 Qtunit = Qt.Unit();
0057
0058   TVector3 Dt = D.Vect();
0059   Dt.SetZ(0);
0060   double tanphi = sqrt(pow(Q.Mag(), 2) + pow(Q.Pt(), 2)) / Q.Mag() * Dt.Dot(Runit) / Dt.Dot(Qtunit);
0061   if (swap)
0062     tanphi = -tanphi;
0063
0064   res[0] = costheta;
0065   res[1] = sin2theta;
0066   res[2] = tanphi;
0067 }
0068
0069 #endif