TrackParam.cxx

Go to the documentation of this file.

// Copyright 2019-2020 CERN and copyright holders of ALICE O2.
// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders.
// All rights not expressly granted are reserved.
//
// This software is distributed under the terms of the GNU General Public
// License v3 (GPL Version 3), copied verbatim in the file "COPYING".
//
// In applying this license CERN does not waive the privileges and immunities
// granted to it by virtue of its status as an Intergovernmental Organization
// or submit itself to any jurisdiction.
 
 
#include "MCHTracking/TrackParam.h"
 
#include <iomanip>
#include <iostream>
 
#include <TMath.h>
 
#include "Framework/Logger.h"
 
#include "DataFormatsMCH/Cluster.h"
 
namespace o2
{
namespace mch
{
 
using namespace std;
 
//_________________________________________________________________________
TrackParam::TrackParam(Double_t z, const Double_t param[5]) : mZ(z)
{
  setParameters(param);
}
 
//_________________________________________________________________________
TrackParam::TrackParam(Double_t z, const Double_t param[5], const Double_t cov[15]) : mZ(z)
{
  setParameters(param);
  setCovariances(cov);
}
 
//_________________________________________________________________________
TrackParam::TrackParam(const TrackParam& tp)
  : mZ(tp.mZ),
    mParameters(tp.mParameters),
    mClusterPtr(tp.mClusterPtr),
    mRemovable(tp.mRemovable),
    mTrackChi2(tp.mTrackChi2),
    mLocalChi2(tp.mLocalChi2)
{
  if (tp.mCovariances) {
    mCovariances = std::make_unique<TMatrixD>(*(tp.mCovariances));
  }
  if (tp.mPropagator) {
    mPropagator = std::make_unique<TMatrixD>(*(tp.mPropagator));
  }
  if (tp.mExtrapParameters) {
    mExtrapParameters = std::make_unique<TMatrixD>(*(tp.mExtrapParameters));
  }
  if (tp.mExtrapCovariances) {
    mExtrapCovariances = std::make_unique<TMatrixD>(*(tp.mExtrapCovariances));
  }
  if (tp.mSmoothParameters) {
    mSmoothParameters = std::make_unique<TMatrixD>(*(tp.mSmoothParameters));
  }
  if (tp.mSmoothCovariances) {
    mSmoothCovariances = std::make_unique<TMatrixD>(*(tp.mSmoothCovariances));
  }
}
 
//_________________________________________________________________________
TrackParam& TrackParam::operator=(const TrackParam& tp)
{
  if (this == &tp) {
    return *this;
  }
 
  mZ = tp.mZ;
 
  mParameters = tp.mParameters;
 
  if (tp.mCovariances) {
    if (mCovariances) {
      *mCovariances = *(tp.mCovariances);
    } else {
      mCovariances = std::make_unique<TMatrixD>(*(tp.mCovariances));
    }
  } else {
    mCovariances.reset();
  }
 
  if (tp.mPropagator) {
    if (mPropagator) {
      *mPropagator = *(tp.mPropagator);
    } else {
      mPropagator = std::make_unique<TMatrixD>(*(tp.mPropagator));
    }
  } else {
    mPropagator.reset();
  }
 
  if (tp.mExtrapParameters) {
    if (mExtrapParameters) {
      *mExtrapParameters = *(tp.mExtrapParameters);
    } else {
      mExtrapParameters = std::make_unique<TMatrixD>(*(tp.mExtrapParameters));
    }
  } else {
    mExtrapParameters.reset();
  }
 
  if (tp.mExtrapCovariances) {
    if (mExtrapCovariances) {
      *mExtrapCovariances = *(tp.mExtrapCovariances);
    } else {
      mExtrapCovariances = std::make_unique<TMatrixD>(*(tp.mExtrapCovariances));
    }
  } else {
    mExtrapCovariances.reset();
  }
 
  if (tp.mSmoothParameters) {
    if (mSmoothParameters) {
      *mSmoothParameters = *(tp.mSmoothParameters);
    } else {
      mSmoothParameters = std::make_unique<TMatrixD>(*(tp.mSmoothParameters));
    }
  } else {
    mSmoothParameters.reset();
  }
 
  if (tp.mSmoothCovariances) {
    if (mSmoothCovariances) {
      *mSmoothCovariances = *(tp.mSmoothCovariances);
    } else {
      mSmoothCovariances = std::make_unique<TMatrixD>(*(tp.mSmoothCovariances));
    }
  } else {
    mSmoothCovariances.reset();
  }
 
  mClusterPtr = tp.mClusterPtr;
 
  mRemovable = tp.mRemovable;
 
  mTrackChi2 = tp.mTrackChi2;
  mLocalChi2 = tp.mLocalChi2;
 
  return *this;
}
 
//__________________________________________________________________________
void TrackParam::clear()
{
  deleteCovariances();
  mPropagator.reset();
  mExtrapParameters.reset();
  mExtrapCovariances.reset();
  mSmoothParameters.reset();
  mSmoothCovariances.reset();
}
 
//__________________________________________________________________________
Double_t TrackParam::px() const
{
  Double_t pZ;
  if (TMath::Abs(mParameters(4, 0)) > 0) {
    Double_t pYZ = TMath::Abs(1.0 / mParameters(4, 0));
    pZ = -pYZ / (TMath::Sqrt(1.0 + mParameters(3, 0) * mParameters(3, 0))); // spectro. (z<0)
  } else {
    pZ = -FLT_MAX / TMath::Sqrt(1.0 + mParameters(3, 0) * mParameters(3, 0) + mParameters(1, 0) * mParameters(1, 0));
  }
  return pZ * mParameters(1, 0);
}
 
//__________________________________________________________________________
Double_t TrackParam::py() const
{
  Double_t pZ;
  if (TMath::Abs(mParameters(4, 0)) > 0) {
    Double_t pYZ = TMath::Abs(1.0 / mParameters(4, 0));
    pZ = -pYZ / (TMath::Sqrt(1.0 + mParameters(3, 0) * mParameters(3, 0))); // spectro. (z<0)
  } else {
    pZ = -FLT_MAX / TMath::Sqrt(1.0 + mParameters(3, 0) * mParameters(3, 0) + mParameters(1, 0) * mParameters(1, 0));
  }
  return pZ * mParameters(3, 0);
}
 
//__________________________________________________________________________
Double_t TrackParam::pz() const
{
  if (TMath::Abs(mParameters(4, 0)) > 0) {
    Double_t pYZ = TMath::Abs(1.0 / mParameters(4, 0));
    return -pYZ / (TMath::Sqrt(1.0 + mParameters(3, 0) * mParameters(3, 0))); // spectro. (z<0)
  } else {
    return -FLT_MAX / TMath::Sqrt(1.0 + mParameters(3, 0) * mParameters(3, 0) + mParameters(1, 0) * mParameters(1, 0));
  }
}
 
//__________________________________________________________________________
Double_t TrackParam::p() const
{
  if (TMath::Abs(mParameters(4, 0)) > 0) {
    Double_t pYZ = TMath::Abs(1.0 / mParameters(4, 0));
    Double_t pZ = -pYZ / (TMath::Sqrt(1.0 + mParameters(3, 0) * mParameters(3, 0))); // spectro. (z<0)
    return -pZ * TMath::Sqrt(1.0 + mParameters(3, 0) * mParameters(3, 0) + mParameters(1, 0) * mParameters(1, 0));
  } else {
    return FLT_MAX;
  }
}
 
//__________________________________________________________________________
const TMatrixD& TrackParam::getCovariances() const
{
  if (!mCovariances) {
    mCovariances = std::make_unique<TMatrixD>(5, 5);
    mCovariances->Zero();
  }
  return *mCovariances;
}
 
//__________________________________________________________________________
void TrackParam::setCovariances(const TMatrixD& covariances)
{
  if (mCovariances) {
    *mCovariances = covariances;
  } else {
    mCovariances = std::make_unique<TMatrixD>(covariances);
  }
}
 
//__________________________________________________________________________
void TrackParam::setCovariances(const Double_t covariances[15])
{
  if (!mCovariances) {
    mCovariances = std::make_unique<TMatrixD>(5, 5);
  }
  for (Int_t i = 0; i < 5; i++) {
    for (Int_t j = 0; j <= i; j++) {
      (*mCovariances)(i, j) = (*mCovariances)(j, i) = covariances[i * (i + 1) / 2 + j];
    }
  }
}
 
//__________________________________________________________________________
void TrackParam::setVariances(const Double_t covariances[15])
{
  static constexpr int varIdx[5] = {0, 2, 5, 9, 14};
  if (!mCovariances) {
    mCovariances = std::make_unique<TMatrixD>(5, 5);
  }
  mCovariances->Zero();
  for (Int_t i = 0; i < 5; i++) {
    (*mCovariances)(i, i) = covariances[varIdx[i]];
  }
}
 
//__________________________________________________________________________
void TrackParam::deleteCovariances()
{
  mCovariances.reset();
}
 
//__________________________________________________________________________
const TMatrixD& TrackParam::getPropagator() const
{
  if (!mPropagator) {
    mPropagator = std::make_unique<TMatrixD>(5, 5);
    mPropagator->UnitMatrix();
  }
  return *mPropagator;
}
 
//__________________________________________________________________________
void TrackParam::resetPropagator()
{
  if (mPropagator) {
    mPropagator->UnitMatrix();
  }
}
 
//__________________________________________________________________________
void TrackParam::updatePropagator(const TMatrixD& propagator)
{
  if (mPropagator) {
    *mPropagator = TMatrixD(propagator, TMatrixD::kMult, *mPropagator);
  } else {
    mPropagator = std::make_unique<TMatrixD>(propagator);
  }
}
 
//__________________________________________________________________________
const TMatrixD& TrackParam::getExtrapParameters() const
{
  if (!mExtrapParameters) {
    mExtrapParameters = std::make_unique<TMatrixD>(5, 1);
    mExtrapParameters->Zero();
  }
  return *mExtrapParameters;
}
 
//__________________________________________________________________________
void TrackParam::setExtrapParameters(const TMatrixD& extrapParameters)
{
  if (mExtrapParameters) {
    *mExtrapParameters = extrapParameters;
  } else {
    mExtrapParameters = std::make_unique<TMatrixD>(extrapParameters);
  }
}
 
//__________________________________________________________________________
const TMatrixD& TrackParam::getExtrapCovariances() const
{
  if (!mExtrapCovariances) {
    mExtrapCovariances = std::make_unique<TMatrixD>(5, 5);
    mExtrapCovariances->Zero();
  }
  return *mExtrapCovariances;
}
 
//__________________________________________________________________________
void TrackParam::setExtrapCovariances(const TMatrixD& extrapCovariances)
{
  if (mExtrapCovariances) {
    *mExtrapCovariances = extrapCovariances;
  } else {
    mExtrapCovariances = std::make_unique<TMatrixD>(extrapCovariances);
  }
}
 
//__________________________________________________________________________
const TMatrixD& TrackParam::getSmoothParameters() const
{
  if (!mSmoothParameters) {
    mSmoothParameters = std::make_unique<TMatrixD>(5, 1);
    mSmoothParameters->Zero();
  }
  return *mSmoothParameters;
}
 
//__________________________________________________________________________
void TrackParam::setSmoothParameters(const TMatrixD& smoothParameters)
{
  if (mSmoothParameters) {
    *mSmoothParameters = smoothParameters;
  } else {
    mSmoothParameters = std::make_unique<TMatrixD>(smoothParameters);
  }
}
 
//__________________________________________________________________________
const TMatrixD& TrackParam::getSmoothCovariances() const
{
  if (!mSmoothCovariances) {
    mSmoothCovariances = std::make_unique<TMatrixD>(5, 5);
    mSmoothCovariances->Zero();
  }
  return *mSmoothCovariances;
}
 
//__________________________________________________________________________
void TrackParam::setSmoothCovariances(const TMatrixD& smoothCovariances)
{
  if (mSmoothCovariances) {
    *mSmoothCovariances = smoothCovariances;
  } else {
    mSmoothCovariances = std::make_unique<TMatrixD>(smoothCovariances);
  }
}
 
//__________________________________________________________________________
Bool_t TrackParam::isCompatibleTrackParam(const TrackParam& trackParam, Double_t sigma2Cut, Double_t& chi2) const
{
 
  // reset chi2 value
  chi2 = 0.;
 
  // ckeck covariance matrices
  if (!mCovariances && !trackParam.mCovariances) {
    LOG(error) << "Covariance matrix must exist for at least one set of parameters";
    return kFALSE;
  }
 
  Double_t maxChi2 = 5. * sigma2Cut * sigma2Cut; // 5 degrees of freedom
 
  // check Z parameters
  if (mZ != trackParam.mZ) {
    LOG(warn) << "Parameters are given at different Z position (" << mZ << " : " << trackParam.mZ
              << "): results are meaningless";
  }
 
  // compute the parameter residuals
  TMatrixD deltaParam(mParameters, TMatrixD::kMinus, trackParam.mParameters);
 
  // build the error matrix
  TMatrixD weight(5, 5);
  if (mCovariances) {
    weight += *mCovariances;
  }
  if (trackParam.mCovariances) {
    weight += *(trackParam.mCovariances);
  }
 
  // invert the error matrix to get the parameter weights if possible
  if (weight.Determinant() == 0) {
    LOG(error) << "Cannot compute the compatibility chi2";
    return kFALSE;
  }
  weight.Invert();
 
  // compute the compatibility chi2
  TMatrixD tmp(deltaParam, TMatrixD::kTransposeMult, weight);
  TMatrixD mChi2(tmp, TMatrixD::kMult, deltaParam);
 
  // set chi2 value
  chi2 = mChi2(0, 0);
 
  // check compatibility
  if (chi2 > maxChi2) {
    return kFALSE;
  }
 
  return kTRUE;
}
 
//__________________________________________________________________________
TrackParamStruct TrackParam::getTrackParamStruct() const
{
 
  TrackParamStruct param{};
 
  param.x = getNonBendingCoor();
  param.y = getBendingCoor();
  param.z = getZ();
  param.px = px();
  param.py = py();
  param.pz = pz();
  param.sign = getCharge();
 
  return param;
}
 
//__________________________________________________________________________
void TrackParam::print() const
{
  cout << "<TrackParam> Bending P=" << setw(5) << setprecision(3) << 1. / mParameters(4, 0)
       << ", NonBendSlope=" << setw(5) << setprecision(3) << mParameters(1, 0) * 180. / TMath::Pi()
       << ", BendSlope=" << setw(5) << setprecision(3) << mParameters(3, 0) * 180. / TMath::Pi() << ", (x,y,z)_IP=("
       << setw(5) << setprecision(3) << mParameters(0, 0) << "," << setw(5) << setprecision(3) << mParameters(2, 0)
       << "," << setw(5) << setprecision(3) << mZ << ") cm, (px,py,pz)=(" << setw(5) << setprecision(3) << px() << ","
       << setw(5) << setprecision(3) << py() << "," << setw(5) << setprecision(3) << pz() << ") GeV/c, "
       << "local chi2=" << getLocalChi2() << endl;
}
 
} // namespace mch
} // namespace o2