tauImpactParameter::MultiProngTauSolver Class Reference

#include <MultiProngTauSolver.h>

Inheritance diagram for tauImpactParameter::MultiProngTauSolver:

Public Types
enum	Ambiguity { zero, minus, plus, NAmbiguity }

Public Member Functions
	MultiProngTauSolver ()
virtual	~MultiProngTauSolver ()

Static Public Member Functions
static void	analyticESolver (TLorentzVector &nu_plus, TLorentzVector &nu_minus, const TLorentzVector &A1, bool &isReal)
static LorentzVectorParticle	estimateNu (const LorentzVectorParticle &a1, const TVector3 &pv, int ambiguity, TLorentzVector &tau)
static void	numericalESolver (TLorentzVector &nu_plus, TLorentzVector &nu_minus, const TLorentzVector &A1, bool &isReal)
static void	quadratic (double &x_plus, double &x_minus, double a, double b, double c, bool &isReal)
static TVectorT< double >	rotateToTauFrame (const TVectorT< double > &inpar)
static bool	setTauDirectionatThetaGJMax (const TLorentzVector &a1, double &theta, double &phi, double scale=1.0)
static void	solveByRotation (const TVector3 &TauDir, const TLorentzVector &A1, TLorentzVector &Tau_plus, TLorentzVector &Tau_minus, TLorentzVector &nu_plus, TLorentzVector &nu_minus, bool &isReal, bool rotateback=true)
static double	thetaGJMax (const TLorentzVector &a1)

Detailed Description

Definition at line 17 of file MultiProngTauSolver.h.

Member Enumeration Documentation

enum tauImpactParameter::MultiProngTauSolver::Ambiguity

Enumerator
zero
minus
plus
NAmbiguity

Definition at line 19 of file MultiProngTauSolver.h.

{zero,minus,plus,NAmbiguity};

Constructor & Destructor Documentation

tauImpactParameter::MultiProngTauSolver::MultiProngTauSolver ( )

inline

Definition at line 22 of file MultiProngTauSolver.h.

{};

virtual tauImpactParameter::MultiProngTauSolver::~MultiProngTauSolver ( )

inlinevirtual

Definition at line 23 of file MultiProngTauSolver.h.

References a, analyticESolver(), b, EnergyCorrector::c, estimateNu(), numericalESolver(), MetAnalyzer::pv(), quadratic(), rotateToTauFrame(), Scenarios_cff::scale, setTauDirectionatThetaGJMax(), solveByRotation(), metsig::tau, theta(), and thetaGJMax().

{};

Member Function Documentation

void MultiProngTauSolver::analyticESolver ( TLorentzVector &  nu_plus, TLorentzVector &  nu_minus, const TLorentzVector &  A1, bool &  isReal  )

static

Definition at line 16 of file MultiProngTauSolver.cc.

References a, b, EnergyCorrector::c, quadratic(), mathSSE::sqrt(), and tauImpactParameter::PDGInfo::tau_mass().

Referenced by solveByRotation(), and ~MultiProngTauSolver().

                                                                                                                               {
  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()));
}

LorentzVectorParticle MultiProngTauSolver::estimateNu ( const LorentzVectorParticle &  a1, const TVector3 &  pv, int  ambiguity, TLorentzVector &  tau  )

static

Definition at line 101 of file MultiProngTauSolver.cc.

References tauImpactParameter::Particle::bField(), funct::cos(), mps_fire::i, tauImpactParameter::LorentzVectorParticle::m, minus, tauImpactParameter::LorentzVectorParticle::NLorentzandVertexPar, tauImpactParameter::PDGInfo::nu_tau, tauImpactParameter::LorentzVectorParticle::NVertex, tauImpactParameter::LorentzVectorParticle::p4(), tauImpactParameter::LorentzVectorParticle::parameter(), plus, MetAnalyzer::pv(), tauImpactParameter::LorentzVectorParticle::px, tauImpactParameter::LorentzVectorParticle::py, tauImpactParameter::LorentzVectorParticle::pz, setTauDirectionatThetaGJMax(), funct::sin(), solveByRotation(), jets_cff::tau1, jets_cff::tau2, theta(), findQualityFiles::v, tauImpactParameter::LorentzVectorParticle::vertex(), tauImpactParameter::LorentzVectorParticle::vertexCov(), tauImpactParameter::LorentzVectorParticle::vx, tauImpactParameter::LorentzVectorParticle::vy, tauImpactParameter::LorentzVectorParticle::vz, and zero.

Referenced by ~MultiProngTauSolver().

                                                                                                                                            {
  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());
} 

void MultiProngTauSolver::numericalESolver ( TLorentzVector &  nu_plus, TLorentzVector &  nu_minus, const TLorentzVector &  A1, bool &  isReal  )

static

Definition at line 27 of file MultiProngTauSolver.cc.

References mps_fire::i, createfilelist::int, min(), mathSSE::sqrt(), metsig::tau, and tauImpactParameter::PDGInfo::tau_mass().

Referenced by ~MultiProngTauSolver().

                                                                                                                                {
  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::quadratic ( double &  x_plus, double &  x_minus, double  a, double  b, double  c, bool &  isReal  )

static

Definition at line 8 of file MultiProngTauSolver.cc.

References EnergyCorrector::c, dttmaxenums::R, and mathSSE::sqrt().

Referenced by analyticESolver(), and ~MultiProngTauSolver().

                                                                                                            {
  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);
}

TVectorT< double > MultiProngTauSolver::rotateToTauFrame ( const TVectorT< double > &  inpar )

static

Definition at line 143 of file MultiProngTauSolver.cc.

References funct::cos(), and funct::sin().

Referenced by ~MultiProngTauSolver().

                                                                                   {
  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;
}

bool MultiProngTauSolver::setTauDirectionatThetaGJMax ( const TLorentzVector &  a1, double &  theta, double &  phi, double  scale = 1.0  )

static

Definition at line 76 of file MultiProngTauSolver.cc.

References a, b, funct::cos(), funct::sin(), metsig::tau, theta(), and thetaGJMax().

Referenced by estimateNu(), and ~MultiProngTauSolver().

                                                                                                                     {
  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;
}

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 = true  )

static

Definition at line 51 of file MultiProngTauSolver.cc.

References analyticESolver(), and theta().

Referenced by estimateNu(), tauImpactParameter::TauA1NuConstrainedFitter::fit(), tauImpactParameter::TauA1NuConstrainedFitter::SolveAmbiguityAnalytically(), and ~MultiProngTauSolver().

                                                                                                     {
  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;
  }
}

double MultiProngTauSolver::thetaGJMax ( const TLorentzVector &  a1 )

static

Definition at line 97 of file MultiProngTauSolver.cc.

References tauImpactParameter::PDGInfo::tau_mass().

Referenced by setTauDirectionatThetaGJMax(), tauImpactParameter::TauA1NuConstrainedFitter::SolveAmbiguityAnalyticallywithRot(), tauImpactParameter::TauA1NuConstrainedFitter::TauRot(), and ~MultiProngTauSolver().

                                                              {
  return asin(( PDGInfo::tau_mass()*PDGInfo::tau_mass()-a1.M2())/(2.0*PDGInfo::tau_mass()*fabs(a1.P())));
}