TauA1NuConstrainedFitter.cc

Go to the documentation of this file.

/* From SimpleFits Package
 * Designed an written by
 * author: Ian M. Nugent
 * Humboldt Foundations
 */
#include <functional>
#include <cmath>
#include "FWCore/MessageLogger/interface/MessageLogger.h"
#include "RecoTauTag/ImpactParameter/interface/TauA1NuConstrainedFitter.h"
#include "RecoTauTag/ImpactParameter/interface/PDGInfo.h"
#include <iostream>

using namespace tauImpactParameter;

TauA1NuConstrainedFitter::TauA1NuConstrainedFitter(unsigned int ambiguity,const LorentzVectorParticle& A1,const TVector3& PVertex, const TMatrixTSym<double>& VertexCov):
  MultiProngTauSolver(),
  ambiguity_(ambiguity)
{
  TLorentzVector Tau(0,0,0,0);
  //dummy substitution not used later
  TVectorT<double> nu_par(LorentzVectorParticle::NLorentzandVertexPar,1);
  TMatrixTSym<double> nu_cov(LorentzVectorParticle::NLorentzandVertexPar); 
  LorentzVectorParticle Nu(nu_par,nu_cov,PDGInfo::nu_tau,0.0,A1.bField()); 
  particles_.push_back(A1);
  particles_.push_back(Nu);
  
  // setup 13 by 13 matrix
  int size=LorentzVectorParticle::NVertex+particles_.size()*LorentzVectorParticle::NLorentzandVertexPar;
  TVectorT<double>    inpar(size);
  TMatrixTSym<double> incov(size);

  // Get primary vertex information
  if(VertexCov.GetNrows()!=LorentzVectorParticle::NVertex)return;
  inpar(LorentzVectorParticle::vx)=PVertex.X();
  inpar(LorentzVectorParticle::vy)=PVertex.Y();
  inpar(LorentzVectorParticle::vz)=PVertex.Z();  
  for(int i=0; i<LorentzVectorParticle::NVertex;i++){
    for(int j=0; j<LorentzVectorParticle::NVertex;j++)incov(i,j)=VertexCov(i,j);
  }
  int A1offset=LorentzVectorParticle::NVertex;
  int Nuoffset=LorentzVectorParticle::NLorentzandVertexPar+LorentzVectorParticle::NVertex;
  for(int i=0; i<LorentzVectorParticle::NLorentzandVertexPar;i++){
    inpar(i+A1offset)=A1.parameter(i);
    inpar(i+Nuoffset)=Nu.parameter(i)+1.0;// offset by 1 GeV to prevent convergence on first iteration
    for(int j=0; j<LorentzVectorParticle::NLorentzandVertexPar;j++){
      incov(i+A1offset,j+A1offset)=A1.covariance(i,j);
      incov(i+Nuoffset,j+Nuoffset)=Nu.covariance(i,j);
    }
  }

  exppar.ResizeTo(nexpandedpar,1);
  exppar=ComputeInitalExpPar(inpar);
  expcov.ResizeTo(nexpandedpar,nexpandedpar);
  expcov=ErrorMatrixPropagator::propagateError(&TauA1NuConstrainedFitter::ComputeInitalExpPar,inpar,incov);

  TVectorT<double> PAR_0(npar);
  par_0.ResizeTo(npar);
  cov_0.ResizeTo(npar,npar);
  PAR_0=ComputeExpParToPar(exppar);
  for(int i=0; i<npar;i++)par_0(i)=PAR_0(i);
  cov_0=ErrorMatrixPropagator::propagateError(&TauA1NuConstrainedFitter::ComputeExpParToPar,exppar,expcov);

  for(int i=0; i<npar;i++){
    for(int j=0;j<npar;j++){cov_0(i,j)=expcov(i,j);}
  }

  par.ResizeTo(npar);
  par=par_0;
  cov.ResizeTo(npar,npar);
  cov=cov_0;
}

TVectorT<double> TauA1NuConstrainedFitter::ComputeInitalExpPar(const TVectorT<double>& inpar){
  TVectorT<double> outpar(nexpandedpar);
  int offset=LorentzVectorParticle::NVertex;// for A1
  TVector3 pv(inpar(LorentzVectorParticle::vx),inpar(LorentzVectorParticle::vy),inpar(LorentzVectorParticle::vz));
  TVector3 sv(inpar(LorentzVectorParticle::vx+offset),inpar(LorentzVectorParticle::vy+offset),inpar(LorentzVectorParticle::vz+offset));
  TVector3 TauDir=sv-pv;
  outpar(tau_phi)=TauDir.Phi();
  outpar(tau_theta)=TauDir.Theta();
  outpar(a1_px)=inpar(LorentzVectorParticle::px+offset);
  outpar(a1_py)=inpar(LorentzVectorParticle::py+offset);
  outpar(a1_pz)=inpar(LorentzVectorParticle::pz+offset);
  outpar(a1_m)=inpar(LorentzVectorParticle::m+offset);
  outpar(a1_vx)=inpar(LorentzVectorParticle::vx+offset);
  outpar(a1_vy)=inpar(LorentzVectorParticle::vy+offset);
  outpar(a1_vz)=inpar(LorentzVectorParticle::vz+offset);
  offset+=LorentzVectorParticle::NLorentzandVertexPar; // for Nu
  outpar(nu_px)=inpar(LorentzVectorParticle::px+offset);
  outpar(nu_py)=inpar(LorentzVectorParticle::py+offset);
  outpar(nu_pz)=inpar(LorentzVectorParticle::pz+offset);
  return outpar;
}

TVectorT<double> TauA1NuConstrainedFitter::ComputeExpParToPar(const TVectorT<double>& inpar){
  TVectorT<double> outpar(npar);
  for(int i=0;i<npar;i++){outpar(i)=inpar(i);}
  return outpar;
}

TVectorT<double> TauA1NuConstrainedFitter::ComputeNuLorentzVectorPar(const TVectorT<double>& inpar){
  TVectorT<double> outpar(LorentzVectorParticle::NLorentzandVertexPar);
  outpar(LorentzVectorParticle::vx)=inpar(a1_vx);
  outpar(LorentzVectorParticle::vy)=inpar(a1_vy);
  outpar(LorentzVectorParticle::vz)=inpar(a1_vz);
  outpar(LorentzVectorParticle::px)=inpar(nu_px);
  outpar(LorentzVectorParticle::py)=inpar(nu_py);
  outpar(LorentzVectorParticle::pz)=inpar(nu_pz);
  outpar(LorentzVectorParticle::m)=0;
  return outpar;
}

TVectorT<double> TauA1NuConstrainedFitter::ComputeA1LorentzVectorPar(const TVectorT<double>& inpar){
  TVectorT<double> outpar(LorentzVectorParticle::NLorentzandVertexPar);
  outpar(LorentzVectorParticle::vx)=inpar(a1_vx);
  outpar(LorentzVectorParticle::vy)=inpar(a1_vy);
  outpar(LorentzVectorParticle::vz)=inpar(a1_vz);
  outpar(LorentzVectorParticle::px)=inpar(a1_px);
  outpar(LorentzVectorParticle::py)=inpar(a1_py);
  outpar(LorentzVectorParticle::pz)=inpar(a1_pz);
  outpar(LorentzVectorParticle::m)=inpar(a1_m);
  return outpar;
}

TVectorT<double> TauA1NuConstrainedFitter::ComputeMotherLorentzVectorPar(const TVectorT<double>& inpar){
  TVectorT<double> outpar(LorentzVectorParticle::NLorentzandVertexPar);
  TVectorT<double> nupar=ComputeNuLorentzVectorPar(inpar);
  TVectorT<double> a1par=ComputeA1LorentzVectorPar(inpar);
  for(int i=0;i<LorentzVectorParticle::NLorentzandVertexPar;i++){
    if(i<LorentzVectorParticle::NVertex){outpar(i)=a1par(i);}
    else{outpar(i)=nupar(i)+a1par(i);}
    //if(i==LorentzVectorParticle::m) outpar(i,0)=PDGInfo::tau_mass();
  }
  double nu_px = nupar(LorentzVectorParticle::px);
  double nu_py = nupar(LorentzVectorParticle::py);
  double nu_pz = nupar(LorentzVectorParticle::pz); 
  double Enu2  = nu_px*nu_px + nu_py*nu_py + nu_pz*nu_pz;
  double a1_px = a1par(LorentzVectorParticle::px);
  double a1_py = a1par(LorentzVectorParticle::py);
  double a1_pz = a1par(LorentzVectorParticle::pz);
  double a1_m =  a1par(LorentzVectorParticle::m);
  double Ea12  = a1_px*a1_px + a1_py*a1_py + a1_pz*a1_pz + a1_m*a1_m;
  double outpar_px = outpar(LorentzVectorParticle::px);
  double outpar_py = outpar(LorentzVectorParticle::py);
  double outpar_pz = outpar(LorentzVectorParticle::pz);
  double P2=outpar_px*outpar_px + outpar_py*outpar_py + outpar_pz*outpar_pz;
   outpar(LorentzVectorParticle::m)=sqrt(std::fabs(Enu2 + Ea12 + 2*sqrt(Enu2*Ea12)-P2));
  return outpar;
}

void TauA1NuConstrainedFitter::UpdateExpandedPar(){
  // assumes changes to a1 correlation to vertex is small
  if(par.GetNrows()==npar && cov.GetNrows()==npar && exppar.GetNrows()==npar && expcov.GetNrows()==npar) return;
  for(int i=0; i<npar;i++){
    exppar(i)=par(i);
    for(int j=0; j<npar;j++){expcov(i,j)=cov(i,j);}
  }
}

std::vector<LorentzVectorParticle> TauA1NuConstrainedFitter::getRefitDaughters(){
  std::vector<LorentzVectorParticle> refitParticles;
  UpdateExpandedPar();
  double c(0),b(0);
  for(unsigned int i=0;i<particles_.size();i++){c+=particles_[i].charge();b=particles_[i].bField();}
  TVectorT<double> a1=ComputeA1LorentzVectorPar(exppar);
  TMatrixTSym<double> a1cov=ErrorMatrixPropagator::propagateError(&TauA1NuConstrainedFitter::ComputeA1LorentzVectorPar,exppar,expcov);
  refitParticles.push_back(LorentzVectorParticle(a1,a1cov,particles_[0].pdgId(),c,b));
  TVectorT<double> nu=ComputeNuLorentzVectorPar(exppar);
  TMatrixTSym<double> nucov=ErrorMatrixPropagator::propagateError(&TauA1NuConstrainedFitter::ComputeNuLorentzVectorPar,exppar,expcov);
  refitParticles.push_back(LorentzVectorParticle(nu,nucov,PDGInfo::nu_tau,0.0,b));
  return refitParticles;
}

LorentzVectorParticle TauA1NuConstrainedFitter::getMother(){
  UpdateExpandedPar();
  double c(0),b(0);
  for(unsigned int i=0;i<particles_.size();i++){c+=particles_[i].charge();b=particles_[i].bField();}
  TVectorT<double> m=ComputeMotherLorentzVectorPar(exppar);
  TMatrixTSym<double> mcov=ErrorMatrixPropagator::propagateError(&TauA1NuConstrainedFitter::ComputeMotherLorentzVectorPar,exppar,expcov);
  LorentzVectorParticle mymother= LorentzVectorParticle(m,mcov,(int)(-1.0*std::abs(PDGInfo::tau_minus)*c),c,b);
  return mymother;
}

void TauA1NuConstrainedFitter::CovertParToObjects(const TVectorD &v,TLorentzVector &a1,TLorentzVector &nu,double &phi,double &theta,TVector3 &TauDir){
  a1=TLorentzVector(v(a1_px),v(a1_py),v(a1_pz),sqrt(v(a1_m)*v(a1_m)+v(a1_px)*v(a1_px)+v(a1_py)*v(a1_py)+v(a1_pz)*v(a1_pz)));
  nu=TLorentzVector(v(nu_px),v(nu_py),v(nu_pz),sqrt(v(nu_px)*v(nu_px)+v(nu_py)*v(nu_py)+v(nu_pz)*v(nu_pz)));
  phi=v(tau_phi);
  theta=v(tau_theta);
  TauDir.SetMagThetaPhi(1.0,theta,phi);
}

bool TauA1NuConstrainedFitter::fit(){
  // Check if Tau Direction is unphysical and if nessicary set the starting point to Theta_{GJ-Max}
  TLorentzVector a1(par(a1_px),par(a1_py),par(a1_pz),sqrt(par(a1_m)*par(a1_m)+par(a1_px)*par(a1_px)+par(a1_py)*par(a1_py)+par(a1_pz)*par(a1_pz)));
  double phi(par(tau_phi)),theta(par(tau_theta));
  TLorentzVector Tau_plus,Tau_minus,nu_plus,nu_minus;
  TVector3 TauDir(cos(phi)*sin(theta),sin(phi)*sin(theta),cos(theta));
  bool isReal;
  solveByRotation(TauDir,a1,Tau_plus,Tau_minus,nu_plus,nu_minus,isReal);

  //check that the do product of the a1 and tau is positive, otherwise there is no information for tau direction -> use zero solution
  if(TauDir.Dot(a1.Vect())<0){
    isReal=false;
  }

  //case 1: is real then solve analytically
  if(isReal && (ambiguity_==plus || ambiguity_==minus)){
    // popogate errors
    TVectorT<double> par_tmp=TauA1NuConstrainedFitter::SolveAmbiguityAnalytically(par,ambiguity_);
    cov=ErrorMatrixPropagator::propagateError(std::bind(TauA1NuConstrainedFitter::SolveAmbiguityAnalytically,std::placeholders::_1,ambiguity_),par,cov_0);
    for(int i=0; i<npar;i++) par(i)=par_tmp(i);
    return true;
  }
  // case 2 is in unphsyical region - rotate and substitue \theta_{GJ} with \theta_{GJ}^{Max} and then solve analytically
  else if(!isReal && ambiguity_==zero){
    TVectorT<double> par_tmp=TauA1NuConstrainedFitter::SolveAmbiguityAnalyticallywithRot(par,ambiguity_);
    cov=ErrorMatrixPropagator::propagateError(std::bind(TauA1NuConstrainedFitter::SolveAmbiguityAnalyticallywithRot,std::placeholders::_1,ambiguity_),par,cov_0);
    for(int i=0; i<npar;i++) par(i)=par_tmp(i);
    return true;
  }
  return false;
}

TVectorT<double> TauA1NuConstrainedFitter::SolveAmbiguityAnalytically(const TVectorT<double>& inpar, unsigned int amb){
  // Solve equation quadratic equation
  TVectorT<double> outpar(inpar.GetNrows());
  TLorentzVector a1,nu;
  double phi(0),theta(0);
  TVector3 TauDir;
  CovertParToObjects(inpar,a1,nu,phi,theta,TauDir);
  TLorentzVector a1_d=a1;
  TLorentzVector nu_d=nu;
  TLorentzVector Tau_plus,Tau_minus,nu_plus,nu_minus;
  bool isReal;
  solveByRotation(TauDir,a1_d,Tau_plus,Tau_minus,nu_plus,nu_minus,isReal,true);
  if(amb==plus)nu=nu_plus;
  else nu=nu_minus;
  for(int i=0; i<outpar.GetNrows();i++){ outpar(i)=inpar(i);}
  outpar(nu_px)=nu.Px(); 
  outpar(nu_py)=nu.Py(); 
  outpar(nu_pz)=nu.Pz(); 
  return outpar;
}

TVectorT<double> TauA1NuConstrainedFitter::SolveAmbiguityAnalyticallywithRot(const TVectorT<double>& inpar, unsigned int ambiguity){
  // Rotate and subsitute \theta_{GJ} with \theta_{GJ}^{Max} - assumes uncertianty on thata and phi of the a1 or small compared to the tau direction. 
  TVectorT<double> outpar(inpar.GetNrows());
  TLorentzVector a1,nu;
  double phi(0),theta(0);
  TVector3 TauDir;
  CovertParToObjects(inpar,a1,nu,phi,theta,TauDir);
  double theta_a1(a1.Theta()),phi_a1(a1.Phi()),theta_GJMax(thetaGJMax(a1));
  TauDir.RotateZ(-phi_a1);
  TauDir.RotateY(-theta_a1);
  double phiprime(TauDir.Phi());
  TauDir=TVector3(sin(theta_GJMax)*cos(phiprime),sin(theta_GJMax)*sin(phiprime),cos(theta_GJMax));
  TauDir.RotateY(theta_a1);
  TauDir.RotateZ(phi_a1);
  for(int i=0; i<outpar.GetNrows();i++) outpar(i)=inpar(i);
  outpar(tau_phi)=TauDir.Phi();
  outpar(tau_theta)=TauDir.Theta();
  return SolveAmbiguityAnalytically(outpar,ambiguity);
}

// Return the significance of the rotation when the tau direction is in the unphysical region
double TauA1NuConstrainedFitter::getTauRotationSignificance(){
  TVectorT<double> par_tmp=TauA1NuConstrainedFitter::TauRot(par);
  TMatrixTSym<double> cov_tmp=ErrorMatrixPropagator::propagateError(&TauA1NuConstrainedFitter::TauRot,par,cov_0);
  if(!(cov_tmp(0,0)>0)) return -999; // return invalid value if the covariance is unphysical (safety flag)
  if(par_tmp(0)>0)    return par_tmp(0)/sqrt(cov_tmp(0,0)); // return the significance if the value is in the unphysical region
  return 0; // reutrn 0 for the rotation significance if the tau is in the physical region
}


TVectorT<double> TauA1NuConstrainedFitter::TauRot(const TVectorT<double>& inpar){
  TVectorT<double> outpar(1);
  TLorentzVector a1,nu;
  double phi(0),theta(0);
  TVector3 TauDir;
  CovertParToObjects(inpar,a1,nu,phi,theta,TauDir);
  double theta_a1(a1.Theta()),phi_a1(a1.Phi()),theta_GJMax(thetaGJMax(a1));
  TauDir.RotateZ(-phi_a1);
  TauDir.RotateY(-theta_a1);
  outpar(0)=(TauDir.Theta()-theta_GJMax);
  return outpar;
}