/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_7/src/Alignment/TwoBodyDecay/src/TwoBodyDecayDerivatives.cc

Go to the documentation of this file.

#include "Alignment/TwoBodyDecay/interface/TwoBodyDecayDerivatives.h"
#include "Alignment/TwoBodyDecay/interface/TwoBodyDecayModel.h"

#include "FWCore/Utilities/interface/Exception.h"

#include <algorithm>

TwoBodyDecayDerivatives::TwoBodyDecayDerivatives( double mPrimary, double mSecondary ) : 
  thePrimaryMass( mPrimary ), theSecondaryMass( mSecondary ) {}


TwoBodyDecayDerivatives::~TwoBodyDecayDerivatives() {}


const std::pair< AlgebraicMatrix, AlgebraicMatrix >
TwoBodyDecayDerivatives::derivatives( const TwoBodyDecay & tbd )  const
{
  return derivatives( tbd.decayParameters() );
}

const std::pair< AlgebraicMatrix, AlgebraicMatrix >
TwoBodyDecayDerivatives::derivatives( const TwoBodyDecayParameters & param )  const
{
  // get the derivatives with respect to all parameters
  std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpx = this->dqsdpx( param );
  std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpy = this->dqsdpy( param );
  std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdpz = this->dqsdpz( param );
  std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdtheta = this->dqsdtheta( param );
  std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdphi = this->dqsdphi( param );
  std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdm = this->dqsdm( param );

  AlgebraicMatrix dqplusdz( 3, dimension );
  dqplusdz.sub( 1, px, dqsdpx.first );
  dqplusdz.sub( 1, py, dqsdpy.first );
  dqplusdz.sub( 1, pz, dqsdpz.first );
  dqplusdz.sub( 1, theta, dqsdtheta.first );
  dqplusdz.sub( 1, phi, dqsdphi.first );
  dqplusdz.sub( 1, mass, dqsdm.first );

  AlgebraicMatrix dqminusdz( 3, dimension );
  dqminusdz.sub( 1, px, dqsdpx.second );
  dqminusdz.sub( 1, py, dqsdpy.second );
  dqminusdz.sub( 1, pz, dqsdpz.second );
  dqminusdz.sub( 1, theta, dqsdtheta.second );
  dqminusdz.sub( 1, phi, dqsdphi.second );
  dqminusdz.sub( 1, mass, dqsdm.second );

  return std::make_pair( dqplusdz, dqminusdz );
}


const std::pair< AlgebraicMatrix, AlgebraicMatrix >
TwoBodyDecayDerivatives::selectedDerivatives( const TwoBodyDecay & tbd, const std::vector< bool > & selector ) const
{
  return selectedDerivatives( tbd.decayParameters(), selector );
}


const std::pair< AlgebraicMatrix, AlgebraicMatrix >
TwoBodyDecayDerivatives::selectedDerivatives( const TwoBodyDecayParameters & param, const std::vector< bool > & selector ) const
{
  if ( selector.size() != dimension )
  {
    throw cms::Exception( "BadConfig" ) << "@SUB=TwoBodyDecayDerivatives::selectedDerivatives"
                                        << "selector has bad dimension (size=" << selector.size() << ").";
  }

  int nSelected = std::count( selector.begin(), selector.end(), true );
  int iSelected = 1;

  AlgebraicMatrix dqplusdz( 3, nSelected );
  AlgebraicMatrix dqminusdz( 3, nSelected );
  std::pair< AlgebraicMatrix, AlgebraicMatrix > dqsdzi;

  for ( unsigned int i = 1; i <= dimension; i++ )
  {
    if ( selector[i] )
    {
      dqsdzi = this->dqsdzi( param, DerivativeParameterName(i) );
      dqplusdz.sub( 1, iSelected, dqsdzi.first );
      dqminusdz.sub( 1, iSelected, dqsdzi.second );
      iSelected++;
    }
  }

  return std::make_pair( dqplusdz, dqminusdz );
}


const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdpx( const TwoBodyDecayParameters & param ) const
{
  double px = param[TwoBodyDecayParameters::px];
  double py = param[TwoBodyDecayParameters::py];
  double pz = param[TwoBodyDecayParameters::pz];
  double theta = param[TwoBodyDecayParameters::theta];
  double phi = param[TwoBodyDecayParameters::phi];

  // compute transverse and absolute momentum
  double pT2 = px*px + py*py;
  double p2 = pT2 + pz*pz;
  double pT = sqrt( pT2 );
  double p = sqrt( p2 );

  double sphi = sin( phi );
  double cphi = cos( phi );
  double stheta = sin( theta );
  double ctheta = cos( theta );

  // some constants from kinematics
  double c1 = 0.5*thePrimaryMass/theSecondaryMass;
  double c2 = sqrt( c1*c1 - 1. );
  double c3 = 0.5*c2*ctheta/c1;
  double c4 = sqrt( p2 + thePrimaryMass*thePrimaryMass );

  // momentum of decay particle 1 in the primary's boosted frame
  AlgebraicMatrix pplus( 3, 1 );
  pplus[0][0] = theSecondaryMass*c2*stheta*cphi;
  pplus[1][0] = theSecondaryMass*c2*stheta*sphi;
  pplus[2][0] = 0.5*p + c3*c4;

  // momentum of decay particle 2 in the primary's boosted frame
  AlgebraicMatrix pminus( 3, 1 );
  pminus[0][0] = -pplus[0][0];
  pminus[1][0] = -pplus[1][0];
  pminus[2][0] = 0.5*p - c3*c4;

  // derivative of rotation matrix w.r.t. px
  AlgebraicMatrix dRotMatdpx( 3, 3 );

  dRotMatdpx[0][0] = pz/(pT*p)*(1.-px*px*(1./pT2+1./p2));
  dRotMatdpx[0][1] = px*py/(pT*pT2);
  dRotMatdpx[0][2] = (1.-px*px/p2)/p;

  dRotMatdpx[1][0] = -px*py*pz/(pT*p)*(1./pT2+1./p2);
  dRotMatdpx[1][1] = (1.-px*px/pT2)/pT;
  dRotMatdpx[1][2] = -px*py/(p*p2);

  dRotMatdpx[2][0] = -(1./pT-pT/p2)*px/p;
  dRotMatdpx[2][1] = 0.;
  dRotMatdpx[2][2] = -px*pz/(p*p2);

  // derivative of the momentum of particle 1 in the lab frame w.r.t. px
  double dpplusdpx = px*( 0.5/p + c3/c4 );

  AlgebraicMatrix dqplusdpx = dRotMatdpx*pplus;
  dqplusdpx[0][0] += px*dpplusdpx/p;
  dqplusdpx[1][0] += py*dpplusdpx/p;
  dqplusdpx[2][0] += pz*dpplusdpx/p;

  // derivative of the momentum of particle 2 in the lab frame w.r.t. px
  double dpminusdpx = px*( 0.5/p - c3/c4 );

  AlgebraicMatrix dqminusdpx = dRotMatdpx*pminus; 
  dqminusdpx[0][0] += px*dpminusdpx/p;
  dqminusdpx[1][0] += py*dpminusdpx/p;
  dqminusdpx[2][0] += pz*dpminusdpx/p;

  // return result
  return std::make_pair( dqplusdpx, dqminusdpx );
}


const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdpy( const TwoBodyDecayParameters & param ) const
{
  double px = param[TwoBodyDecayParameters::px];
  double py = param[TwoBodyDecayParameters::py];
  double pz = param[TwoBodyDecayParameters::pz];
  double theta = param[TwoBodyDecayParameters::theta];
  double phi = param[TwoBodyDecayParameters::phi];

  // compute transverse and absolute momentum
  double pT2 = px*px + py*py;
  double p2 = pT2 + pz*pz;
  double pT = sqrt( pT2 );
  double p = sqrt( p2 );

  double sphi = sin( phi );
  double cphi = cos( phi );
  double stheta = sin( theta );
  double ctheta = cos( theta );

  // some constants from kinematics
  double c1 = 0.5*thePrimaryMass/theSecondaryMass;
  double c2 = sqrt( c1*c1 - 1. );
  double c3 = 0.5*c2*ctheta/c1;
  double c4 = sqrt( p2 + thePrimaryMass*thePrimaryMass );

  // momentum of decay particle 1 in the rest frame of the primary
  AlgebraicMatrix pplus( 3, 1 );
  pplus[0][0] = theSecondaryMass*c2*stheta*cphi;
  pplus[1][0] = theSecondaryMass*c2*stheta*sphi;
  pplus[2][0] = 0.5*p + c3*c4;

  // momentum of decay particle 2 in the rest frame of the primary
  AlgebraicMatrix pminus( 3, 1 );
  pminus[0][0] = -pplus[0][0];
  pminus[1][0] = -pplus[1][0];
  pminus[2][0] = 0.5*p - c3*c4;

  // derivative of rotation matrix w.r.t. py
  AlgebraicMatrix dRotMatdpy( 3, 3 );

  dRotMatdpy[0][0] = -px*py*pz/(pT*p)*(1./pT2+1./p2);
  dRotMatdpy[0][1] = (py*py/pT2-1.)/pT;
  dRotMatdpy[0][2] = -px*py/(p*p2);

  dRotMatdpy[1][0] = pz/(pT*p)*(1.-py*py*(1./pT2+1./p2));
  dRotMatdpy[1][1] = -px*py/(pT*pT2);
  dRotMatdpy[1][2] = (1.-py*py/p2)/p;

  dRotMatdpy[2][0] = -(1./pT-pT/p2)*py/p;
  dRotMatdpy[2][1] = 0.;
  dRotMatdpy[2][2] = -py*pz/(p*p2);

  // derivative of the momentum of particle 1 in the lab frame w.r.t. py
  double dpplusdpy = py*( 0.5/p + c3/c4 );

  AlgebraicMatrix dqplusdpy = dRotMatdpy*pplus;
  dqplusdpy[0][0] += px*dpplusdpy/p;
  dqplusdpy[1][0] += py*dpplusdpy/p;
  dqplusdpy[2][0] += pz*dpplusdpy/p;

  // derivative of the momentum of particle 2 in the lab frame w.r.t. py
  double dpminusdpy = py*( 0.5/p - c3/c4 );

  AlgebraicMatrix dqminusdpy = dRotMatdpy*pminus;
  dqminusdpy[0][0] += px*dpminusdpy/p;
  dqminusdpy[1][0] += py*dpminusdpy/p;
  dqminusdpy[2][0] += pz*dpminusdpy/p;

  // return result
  return std::make_pair( dqplusdpy, dqminusdpy );
}


const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdpz( const TwoBodyDecayParameters & param ) const
{
  double px = param[TwoBodyDecayParameters::px];
  double py = param[TwoBodyDecayParameters::py];
  double pz = param[TwoBodyDecayParameters::pz];
  double theta = param[TwoBodyDecayParameters::theta];
  double phi = param[TwoBodyDecayParameters::phi];

  // compute transverse and absolute momentum
  double pT2 = px*px + py*py;
  double p2 = pT2 + pz*pz;
  double pT = sqrt( pT2 );
  double p = sqrt( p2 );

  double sphi = sin( phi );
  double cphi = cos( phi );
  double stheta = sin( theta );
  double ctheta = cos( theta );

  // some constants from kinematics
  double c1 = 0.5*thePrimaryMass/theSecondaryMass;
  double c2 = sqrt( c1*c1 - 1. );
  double c3 = 0.5*c2*ctheta/c1;
  double c4 = sqrt( p2 + thePrimaryMass*thePrimaryMass );

  // momentum of decay particle 1 in the rest frame of the primary
  AlgebraicMatrix pplus( 3, 1 );
  pplus[0][0] = theSecondaryMass*c2*stheta*cphi;
  pplus[1][0] = theSecondaryMass*c2*stheta*sphi;
  pplus[2][0] = 0.5*p + c3*c4;

  // momentum of decay particle 2 in the rest frame of the primary
  AlgebraicMatrix pminus( 3, 1 );
  pminus[0][0] = -pplus[0][0];
  pminus[1][0] = -pplus[1][0];
  pminus[2][0] = 0.5*p - c3*c4;

  // derivative of rotation matrix w.r.t. py
  AlgebraicMatrix dRotMatdpz( 3, 3 );

  dRotMatdpz[0][0] = px/(pT*p)*(1.-pz*pz/p2);
  dRotMatdpz[0][1] = 0.;
  dRotMatdpz[0][2] = -px*pz/(p*p2);

  dRotMatdpz[1][0] = py/(p*pT)*(1.-pz*pz/p2);
  dRotMatdpz[1][1] = 0.;
  dRotMatdpz[1][2] = -py*pz/(p*p2);

  dRotMatdpz[2][0] = pT*pz/(p*p2);
  dRotMatdpz[2][1] = 0.;
  dRotMatdpz[2][2] = (1.-pz*pz/p2)/p;

  // derivative of the momentum of particle 1 in the lab frame w.r.t. pz
  double dpplusdpz = pz*( 0.5/p + c3/c4 );

  AlgebraicMatrix dqplusdpz = dRotMatdpz*pplus;
  dqplusdpz[0][0] += px*dpplusdpz/p;
  dqplusdpz[1][0] += py*dpplusdpz/p;
  dqplusdpz[2][0] += pz*dpplusdpz/p;

  // derivative of the momentum of particle 2 in the lab frame w.r.t. pz
  double dpminusdpz = pz*( 0.5/p - c3/c4 );

  AlgebraicMatrix dqminusdpz = dRotMatdpz*pminus;
  dqminusdpz[0][0] += px*dpminusdpz/p;
  dqminusdpz[1][0] += py*dpminusdpz/p;
  dqminusdpz[2][0] += pz*dpminusdpz/p;

  // return result
  return std::make_pair( dqplusdpz, dqminusdpz );
}


const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdtheta( const TwoBodyDecayParameters & param ) const
{
  double px = param[TwoBodyDecayParameters::px];
  double py = param[TwoBodyDecayParameters::py];
  double pz = param[TwoBodyDecayParameters::pz];
  double theta = param[TwoBodyDecayParameters::theta];
  double phi = param[TwoBodyDecayParameters::phi];

  // compute transverse and absolute momentum
  double pT2 = px*px + py*py;
  double p2 = pT2 + pz*pz;

  double sphi = sin( phi );
  double cphi = cos( phi );
  double stheta = sin( theta );
  double ctheta = cos( theta );

  // some constants from kinematics
  double c1 = 0.5*thePrimaryMass/theSecondaryMass;
  double c2 = sqrt( c1*c1 - 1. );
  double c3 = -0.5*c2*stheta/c1;
  double c4 = sqrt( p2 + thePrimaryMass*thePrimaryMass );

  // derivative of the momentum of particle 1 in the primary's rest frame w.r.t. theta
  AlgebraicMatrix dpplusdtheta( 3, 1 );
  dpplusdtheta[0][0] = theSecondaryMass*c2*ctheta*cphi;
  dpplusdtheta[1][0] = theSecondaryMass*c2*ctheta*sphi;
  dpplusdtheta[2][0] = c3*c4;

  // derivative of the momentum of particle 2 in the primary's rest frame w.r.t. theta
  AlgebraicMatrix dpminusdtheta( 3, 1 );
  dpminusdtheta[0][0] = -theSecondaryMass*c2*ctheta*cphi;
  dpminusdtheta[1][0] = -theSecondaryMass*c2*ctheta*sphi;
  dpminusdtheta[2][0] = -c3*c4;

  TwoBodyDecayModel decayModel;
  AlgebraicMatrix rotMat = decayModel.rotationMatrix( px, py, pz );

  AlgebraicMatrix dqplusdtheta = rotMat*dpplusdtheta;
  AlgebraicMatrix dqminusdtheta = rotMat*dpminusdtheta;

  return std::make_pair( dqplusdtheta, dqminusdtheta );
}


const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdphi( const TwoBodyDecayParameters & param ) const
{
  double px = param[TwoBodyDecayParameters::px];
  double py = param[TwoBodyDecayParameters::py];
  double pz = param[TwoBodyDecayParameters::pz];
  double theta = param[TwoBodyDecayParameters::theta];
  double phi = param[TwoBodyDecayParameters::phi];

  double sphi = sin( phi );
  double cphi = cos( phi );
  double stheta = sin( theta );

  // some constants from kinematics
  double c1 = 0.5*thePrimaryMass/theSecondaryMass;
  double c2 = sqrt( c1*c1 - 1. );

  // derivative of the momentum of particle 1 in the primary's rest frame w.r.t. phi
  AlgebraicMatrix dpplusdphi( 3, 1 );
  dpplusdphi[0][0] = -theSecondaryMass*c2*stheta*sphi;
  dpplusdphi[1][0] = theSecondaryMass*c2*stheta*cphi;
  dpplusdphi[2][0] = 0.;

  // derivative of the momentum of particle 2 in the primary's rest frame w.r.t. phi
  AlgebraicMatrix dpminusdphi( 3, 1 );
  dpminusdphi[0][0] = theSecondaryMass*c2*stheta*sphi;
  dpminusdphi[1][0] = -theSecondaryMass*c2*stheta*cphi;
  dpminusdphi[2][0] = 0.;

  TwoBodyDecayModel decayModel;
  AlgebraicMatrix rotMat = decayModel.rotationMatrix( px, py, pz );

  AlgebraicMatrix dqplusdphi = rotMat*dpplusdphi;
  AlgebraicMatrix dqminusdphi = rotMat*dpminusdphi;

  return std::make_pair( dqplusdphi, dqminusdphi );
}


const std::pair< AlgebraicMatrix, AlgebraicMatrix > TwoBodyDecayDerivatives::dqsdm( const TwoBodyDecayParameters & param ) const
{
  double px = param[TwoBodyDecayParameters::px];
  double py = param[TwoBodyDecayParameters::py];
  double pz = param[TwoBodyDecayParameters::pz];
  double theta = param[TwoBodyDecayParameters::theta];
  double phi = param[TwoBodyDecayParameters::phi];

  double pT2 = px*px + py*py;
  double p2 = pT2 + pz*pz;

  double sphi = sin( phi );
  double cphi = cos( phi );
  double ctheta = cos( theta );
  double stheta = sin( theta );

  // some constants from kinematics
  double c1 = 0.5*thePrimaryMass/theSecondaryMass;
  double c2 = 1./sqrt( c1*c1 - 1. );
  double m2 = thePrimaryMass*thePrimaryMass;

  // derivative of the momentum of particle 1 in the primary's rest frame w.r.t. the primary's mass
  AlgebraicMatrix dpplusdm( 3, 1 );
  dpplusdm[0][0] = c2*0.5*c1*stheta*cphi;
  dpplusdm[1][0] = c2*0.5*c1*stheta*sphi;
  dpplusdm[2][0] = c2*theSecondaryMass*( c1*c1 + p2/m2 )/sqrt( p2 + m2 )*ctheta;

  // derivative of the momentum of particle 2 in the primary's rest frame w.r.t. the primary's mass
  AlgebraicMatrix dpminusdm( 3, 1 );
  dpminusdm[0][0] = -dpplusdm[0][0];
  dpminusdm[1][0] = -dpplusdm[1][0];
  dpminusdm[2][0] = -dpplusdm[2][0];

  TwoBodyDecayModel decayModel;
  AlgebraicMatrix rotMat = decayModel.rotationMatrix( px, py, pz );

  AlgebraicMatrix dqplusdm = rotMat*dpplusdm;
  AlgebraicMatrix dqminusdm = rotMat*dpminusdm;

  return std::make_pair( dqplusdm, dqminusdm );
}


const std::pair< AlgebraicMatrix, AlgebraicMatrix >
TwoBodyDecayDerivatives::dqsdzi( const TwoBodyDecayParameters & param, const DerivativeParameterName & i ) const
{
  switch ( i )
  {
  case TwoBodyDecayDerivatives::px :
    return dqsdpx( param );
    break;
  case TwoBodyDecayDerivatives::py :
    return dqsdpy( param );
    break;
  case TwoBodyDecayDerivatives::pz :
    return dqsdpz( param );
    break;
  case TwoBodyDecayDerivatives::theta :
    return dqsdtheta( param );
    break;
  case TwoBodyDecayDerivatives::phi :
    return dqsdphi( param );
    break;
  case TwoBodyDecayDerivatives::mass :
    return dqsdm( param );
    break;
  default:
    throw cms::Exception( "BadConfig" ) << "@SUB=TwoBodyDecayDerivatives::dqsdzi"
                                        << "no decay parameter related to selector (" << i << ").";
  };

  return std::make_pair( AlgebraicMatrix( 3, 1, 0 ), AlgebraicMatrix( 3, 1, 0 ) );
}