Kinematics.cxx

Go to the documentation of this file.

#include "NAGASH/Kinematics.h"
 
using namespace NAGASH;
using namespace std;
 
void Kinematics::GetCSFAngles(const TLorentzVector &lep1, const int &charge1, const TLorentzVector &lep2, double ebeam, double &costh, double &phi)
{
 
    TLorentzVector boson = lep1 + lep2;
    double Lplus = (charge1 < 0) ? lep1.E() + lep1.Pz() : lep2.E() + lep2.Pz();
    double Lminus = (charge1 < 0) ? lep1.E() - lep1.Pz() : lep2.E() - lep2.Pz();
    double Pplus = (charge1 < 0) ? lep2.E() + lep2.Pz() : lep1.E() + lep1.Pz();
    double Pminus = (charge1 < 0) ? lep2.E() - lep2.Pz() : lep1.E() - lep1.Pz();
 
    costh = (Lplus * Pminus - Lminus * Pplus);
    costh *= TMath::Abs(boson.Pz());
    costh /= (boson.Mag() * boson.Pz());
    costh /= TMath::Sqrt(boson.Mag2() + boson.Pt() * boson.Pt());
 
    TVector3 boostV = -boson.BoostVector();
    TLorentzVector lep1_boosted = (charge1 < 0) ? lep1 : lep2;
    lep1_boosted.Boost(boostV);
 
    TVector3 CSAxis, xAxis, yAxis;
    TLorentzVector p1, p2;
    double sign = +1.;
    if (boson.Z() < 0)
        sign = -1.;
    p1.SetXYZM(0., 0., sign * ebeam, 0.938);  // quark (?)
    p2.SetXYZM(0., 0., -sign * ebeam, 0.938); // antiquark (?)
 
    p1.Boost(boostV);
    p2.Boost(boostV);
    CSAxis = (p1.Vect().Unit() - p2.Vect().Unit()).Unit();
    yAxis = (p1.Vect().Unit()).Cross(p2.Vect().Unit());
    yAxis = yAxis.Unit();
    xAxis = yAxis.Cross(CSAxis);
    xAxis = xAxis.Unit();
 
    phi = TMath::ATan2(lep1_boosted.Vect() * yAxis, lep1_boosted.Vect() * xAxis);
}
 
void Kinematics::GetBDFAngles(const TLorentzVector &lep1, const int &charge1, const TLorentzVector &lep2, double &costh, double &phi)
{
    const TLorentzVector dilep = lep1 + lep2;
 
    double cosThetap = (charge1 > 0) ? TMath::TanH((lep2.Eta() - lep1.Eta()) / 2.) : TMath::TanH((lep1.Eta() - lep2.Eta()) / 2.);
    costh = (dilep.Rapidity() < 0) ? -cosThetap : cosThetap;
 
    double sintheta = TMath::Sqrt(1. - costh * costh);
    double phiacop = TMath::Pi() - TVector2::Phi_mpi_pi(lep1.Phi() - lep2.Phi());
    phi = TMath::Tan(phiacop / 2) * sintheta;
}
 
void Kinematics::GetVHFAngles(const TLorentzVector &lep1, const int &charge1, const TLorentzVector &lep2, double &costh, double &phi)
{
 
    TLorentzVector boson = lep1 + lep2;
    TVector3 boostV = boson.BoostVector();
    TLorentzVector lep1_boosted = (charge1 > 0) ? lep2 : lep1;
    lep1_boosted.Boost(-boostV);
    phi = TVector2::Phi_mpi_pi(lep1_boosted.Phi() - boson.Phi());
    double theta = lep1_boosted.Angle(boson.Vect());
 
    costh = TMath::Cos(theta);
}
 
void Kinematics::GetAiPolynominals(double costh, double phi, std::vector<double> &aipols)
{
 
    aipols.resize(8);
 
    double theta = TMath::ACos(costh);
    double sintheta = TMath::Sin(theta);
    double sin2theta = TMath::Sin(2.0 * theta);
 
    double cosphi = TMath::Cos(phi);
    double cos2phi = TMath::Cos(2.0 * phi);
    double sinphi = TMath::Sin(phi);
    double sin2phi = TMath::Sin(2.0 * phi);
 
    aipols[0] = 0.5 - 1.5 * costh * costh;
    aipols[1] = sin2theta * cosphi;
    aipols[2] = sintheta * sintheta * cos2phi;
    aipols[3] = sintheta * cosphi;
    aipols[4] = costh;
    aipols[5] = sintheta * sintheta * sin2phi;
    aipols[6] = sin2theta * sinphi;
    aipols[7] = sintheta * sinphi;
}
 
std::vector<double> Kinematics::GetAiPolynominals(double costh, double phi)
{
    std::vector<double> aipols(8, 0);
 
    double theta = TMath::ACos(costh);
    double sintheta = TMath::Sin(theta);
    double sin2theta = TMath::Sin(2.0 * theta);
 
    double cosphi = TMath::Cos(phi);
    double cos2phi = TMath::Cos(2.0 * phi);
    double sinphi = TMath::Sin(phi);
    double sin2phi = TMath::Sin(2.0 * phi);
 
    aipols[0] = 0.5 - 1.5 * costh * costh;
    aipols[1] = sin2theta * cosphi;
    aipols[2] = sintheta * sintheta * cos2phi;
    aipols[3] = sintheta * cosphi;
    aipols[4] = costh;
    aipols[5] = sintheta * sintheta * sin2phi;
    aipols[6] = sin2theta * sinphi;
    aipols[7] = sintheta * sinphi;
 
    return aipols;
}