MultiProngTauSolver.cc

Go to the documentation of this file.

#include "RecoTauTag/ImpactParameter/interface/MultiProngTauSolver.h"
#include "TMatrixTSym.h"
#include "TVectorT.h"
#include "FWCore/MessageLogger/interface/MessageLogger.h"

using namespace tauImpactParameter;

void MultiProngTauSolver::quadratic(double& x_plus, double& x_minus, double a, double b, double c, bool& isReal) {
  double R = b * b - 4 * a * c;
  isReal = true;
  if (R < 0) {
    isReal = false;
  }                                            // flag cases when R<0 but compute quadratic equation with |R|
  x_minus = (-b + sqrt(fabs(R))) / (2.0 * a);  // opposite sign is smaller
  x_plus = (-b - sqrt(fabs(R))) / (2.0 * a);
}

void MultiProngTauSolver::analyticESolver(TLorentzVector& nu_plus,
                                          TLorentzVector& nu_minus,
                                          const TLorentzVector& A1,
                                          bool& isReal) {
  double a = (A1.Pz() * A1.Pz()) / (A1.E() * A1.E()) - 1.0;
  double K = (PDGInfo::tau_mass() * PDGInfo::tau_mass() - A1.M2() - 2.0 * A1.Pt() * A1.Pt()) / (2.0 * A1.E());
  double b = 2.0 * K * A1.Pz() / A1.E();
  double c = K * K - A1.Pt() * A1.Pt();
  double z_plus(0), z_minus(0);
  quadratic(z_plus, z_minus, a, b, c, isReal);
  nu_plus = TLorentzVector(-A1.Px(), -A1.Py(), z_plus, sqrt(z_plus * z_plus + A1.Pt() * A1.Pt()));
  nu_minus = TLorentzVector(-A1.Px(), -A1.Py(), z_minus, sqrt(z_minus * z_minus + A1.Pt() * A1.Pt()));
}

void MultiProngTauSolver::numericalESolver(TLorentzVector& nu_plus,
                                           TLorentzVector& nu_minus,
                                           const TLorentzVector& A1,
                                           bool& isReal) {
  double rmin(-100), rmax(100), step(0.01), mtau2(PDGInfo::tau_mass() * PDGInfo::tau_mass()), z1(-999), z2(-999),
      zmin(-999), min(9999), prev(9999);
  double z = rmin;
  TLorentzVector nu, tau;
  for (int i = 0; i <= (int)(rmax - rmin) / step; i++) {
    nu.SetPxPyPzE(-A1.Px(), -A1.Py(), z, sqrt(z * z + A1.Pt() * A1.Pt()));
    tau = A1 + nu;
    double m2 = tau.M2();
    if (m2 - mtau2 < 0 && prev - mtau2 >= 0)
      z1 = z;
    if (m2 - mtau2 > 0 && prev - mtau2 <= 0)
      z2 = z;
    if (min > m2) {
      zmin = z;
      min = m2;
    }
    prev = m2;
    z += step;
  }
  if (z1 != -999 && z2 != -999) {
    nu_plus = TLorentzVector(-A1.Px(), -A1.Py(), z1, sqrt(z1 * z1 + A1.Pt() * A1.Pt()));
    nu_minus = TLorentzVector(-A1.Px(), -A1.Py(), z2, sqrt(z2 * z2 + A1.Pt() * A1.Pt()));
  } else {
    nu_plus = TLorentzVector(-A1.Px(), -A1.Py(), zmin, sqrt(zmin * zmin + A1.Pt() * A1.Pt()));
    nu_minus = TLorentzVector(-A1.Px(), -A1.Py(), zmin, sqrt(zmin * zmin + A1.Pt() * A1.Pt()));
  }
}

void MultiProngTauSolver::solveByRotation(const TVector3& TauDir,
                                          const TLorentzVector& A1,
                                          TLorentzVector& Tau_plus,
                                          TLorentzVector& Tau_minus,
                                          TLorentzVector& nu_plus,
                                          TLorentzVector& nu_minus,
                                          bool& isReal,
                                          bool rotateback) {
  TLorentzVector A1rot = A1;
  double phi(TauDir.Phi()), theta(TauDir.Theta());
  A1rot.RotateZ(-phi);
  A1rot.RotateY(-theta);
  //  NumericalESolver(nu_plus,nu_minus,A1rot); // for debugging analyticESolver (slow)
  analyticESolver(nu_plus, nu_minus, A1rot, isReal);
  if (rotateback) {
    nu_plus.RotateY(theta);
    nu_plus.RotateZ(phi);
    Tau_plus = A1 + nu_plus;
    //
    nu_minus.RotateY(theta);
    nu_minus.RotateZ(phi);
    Tau_minus = A1 + nu_minus;
  } else {
    Tau_plus = A1rot + nu_plus;
    Tau_minus = A1rot + nu_minus;
  }
}

bool MultiProngTauSolver::setTauDirectionatThetaGJMax(const TLorentzVector& a1,
                                                      double& theta,
                                                      double& phi,
                                                      double scale) {
  double thetaGJMaxvar = thetaGJMax(a1);
  TVector3 a1v(a1.Vect());
  if (a1v.Mag() != 0)
    a1v *= 1 / a1v.Mag();
  TVector3 tau(cos(phi) * sin(theta), sin(phi) * sin(theta), cos(theta));
  double dphitheta = acos(a1v.Dot(tau) / (a1v.Mag() * tau.Mag()));
  if (thetaGJMaxvar < dphitheta || scale < 0) {
    if (scale < 0)
      scale = 1.0;
    double a = (thetaGJMaxvar / dphitheta) - (1 - scale);
    double b = 1 - (thetaGJMaxvar / dphitheta) + (1 - scale);
    edm::LogInfo("RecoTauTag/ImpactParameter")
        << "SetTauDirectionatThetaGJMax before GF " << thetaGJMaxvar << " dot "
        << acos(a1v.Dot(tau) / (a1v.Mag() * tau.Mag())) << " a1 phi " << a1v.Phi() << " tau phi " << tau.Phi()
        << " a1 theta " << a1v.Theta() << " tau theta " << tau.Theta();
    tau *= a;
    a1v *= b;
    tau += a1v;
    theta = tau.Theta();
    phi = tau.Phi();
    edm::LogInfo("RecoTauTag/ImpactParameter")
        << "SetTauDirectionatThetaGJMax GF " << thetaGJMaxvar << " dot " << acos(a1v.Dot(tau) / (a1v.Mag() * tau.Mag()))
        << " phi " << phi << " theta " << theta;
    return true;
  }
  return false;
}

double MultiProngTauSolver::thetaGJMax(const TLorentzVector& a1) {
  return asin((PDGInfo::tau_mass() * PDGInfo::tau_mass() - a1.M2()) / (2.0 * PDGInfo::tau_mass() * fabs(a1.P())));
}

LorentzVectorParticle MultiProngTauSolver::estimateNu(const LorentzVectorParticle& a1,
                                                      const TVector3& pv,
                                                      int ambiguity,
                                                      TLorentzVector& tau) {
  TLorentzVector lorentzA1 = a1.p4();
  TVector3 sv = a1.vertex();
  TVector3 tauFlghtDir = sv - pv;
  TLorentzVector nuGuess;

  TVector3 startingtauFlghtDir = tauFlghtDir.Unit();
  if (ambiguity == zero) {
    double theta = tauFlghtDir.Theta();
    double phi = tauFlghtDir.Phi();
    setTauDirectionatThetaGJMax(lorentzA1, theta, phi);
    startingtauFlghtDir = TVector3(sin(theta) * cos(phi), sin(theta) * sin(phi), cos(theta));
  }
  TLorentzVector tau1, tau2, nu1, nu2;
  bool isReal;
  solveByRotation(startingtauFlghtDir, lorentzA1, tau1, tau2, nu1, nu2, isReal);
  if (ambiguity == plus) {
    nuGuess = nu1;
    tau = tau1;
  }
  if (ambiguity == minus) {
    nuGuess = nu2;
    tau = tau1;
  }
  if (ambiguity == zero) {
    nuGuess = nu1;
    tau = tau1;
  }
  TVectorT<double> par(LorentzVectorParticle::NLorentzandVertexPar);
  par(LorentzVectorParticle::vx) = a1.parameter(LorentzVectorParticle::vx);
  par(LorentzVectorParticle::vy) = a1.parameter(LorentzVectorParticle::vy);
  par(LorentzVectorParticle::vz) = a1.parameter(LorentzVectorParticle::vz);
  par(LorentzVectorParticle::px) = nuGuess.Px();
  par(LorentzVectorParticle::py) = nuGuess.Py();
  par(LorentzVectorParticle::pz) = nuGuess.Pz();
  par(LorentzVectorParticle::m) = nuGuess.M();
  TMatrixTSym<double> Cov(LorentzVectorParticle::NLorentzandVertexPar);
  TMatrixTSym<double> pvCov = a1.vertexCov();
  for (int i = 0; i < LorentzVectorParticle::NLorentzandVertexPar; i++) {
    for (int j = 0; j <= i; j++) {
      if (i < LorentzVectorParticle::NVertex)
        Cov(i, j) = pvCov(i, j);
      else
        Cov(i, j) = 0;
    }
    double v = 0;
    if (i == LorentzVectorParticle::px || i == LorentzVectorParticle::py || i == LorentzVectorParticle::pz)
      v = 10 * par(i) * par(i);
    if (v < 1000.0)
      v = 1000.0;  // try lowing to test impact
    Cov(i, i) += v;
  }
  return LorentzVectorParticle(par, Cov, PDGInfo::nu_tau, 0, a1.bField());
}

TVectorT<double> MultiProngTauSolver::rotateToTauFrame(const TVectorT<double>& inpar) {
  TVectorT<double> outpar(3);
  TVector3 res(inpar(0), inpar(1), inpar(2));
  TVector3 Uz(sin(inpar(4)) * cos(inpar(3)), sin(inpar(4)) * sin(inpar(3)), cos(inpar(4)));
  res.RotateUz(Uz);
  /*  double phi=inpar(3,0);
  double theta=inpar(4,0);
  res.RotateZ(-phi);
  TVector3 Y(0,1,0); 
  TVector3 thetadir=res.Cross(Y); 
  thetadir.RotateY(-theta);
  res.RotateY(-theta);
  res.RotateZ(thetadir.Phi());*/
  outpar(0) = res.X();
  outpar(1) = res.Y();
  outpar(2) = res.Z();
  return outpar;
}