EventShapeVariables.cc

Go to the documentation of this file.

#include "PhysicsTools/CandUtils/interface/EventShapeVariables.h"

#include "TMath.h"

EventShapeVariables::EventShapeVariables(const edm::View<reco::Candidate>& inputVectors) : eigenVectors_(3, 3) {
  inputVectors_.reserve(inputVectors.size());
  for (const auto& vec : inputVectors) {
    inputVectors_.push_back(math::XYZVector(vec.px(), vec.py(), vec.pz()));
  }
  //default values
  set_r(2.);
  setFWmax(10);
}

EventShapeVariables::EventShapeVariables(const std::vector<math::XYZVector>& inputVectors)
    : inputVectors_(inputVectors), eigenVectors_(3, 3) {
  //default values
  set_r(2.);
  setFWmax(10);
}

EventShapeVariables::EventShapeVariables(const std::vector<math::RhoEtaPhiVector>& inputVectors) : eigenVectors_(3, 3) {
  inputVectors_.reserve(inputVectors.size());
  for (const auto& vec : inputVectors) {
    inputVectors_.push_back(math::XYZVector(vec.x(), vec.y(), vec.z()));
  }
  //default values
  set_r(2.);
  setFWmax(10);
}

EventShapeVariables::EventShapeVariables(const std::vector<math::RThetaPhiVector>& inputVectors) : eigenVectors_(3, 3) {
  inputVectors_.reserve(inputVectors.size());
  for (const auto& vec : inputVectors) {
    inputVectors_.push_back(math::XYZVector(vec.x(), vec.y(), vec.z()));
  }
  //default values
  set_r(2.);
  setFWmax(10);
}

double EventShapeVariables::isotropy(const unsigned int& numberOfSteps) const {
  const double deltaPhi = 2 * TMath::Pi() / numberOfSteps;
  double phi = 0, eIn = -1., eOut = -1.;
  for (unsigned int i = 0; i < numberOfSteps; ++i) {
    phi += deltaPhi;
    double sum = 0;
    double cosphi = TMath::Cos(phi);
    double sinphi = TMath::Sin(phi);
    for (const auto& vec : inputVectors_) {
      // sum over inner product of unit vectors and momenta
      sum += TMath::Abs(cosphi * vec.x() + sinphi * vec.y());
    }
    if (eOut < 0. || sum < eOut)
      eOut = sum;
    if (eIn < 0. || sum > eIn)
      eIn = sum;
  }
  return (eIn - eOut) / eIn;
}

double EventShapeVariables::circularity(const unsigned int& numberOfSteps) const {
  const double deltaPhi = 2 * TMath::Pi() / numberOfSteps;
  double circularity = -1, phi = 0, area = 0;
  for (const auto& vec : inputVectors_) {
    area += TMath::Sqrt(vec.x() * vec.x() + vec.y() * vec.y());
  }
  for (unsigned int i = 0; i < numberOfSteps; ++i) {
    phi += deltaPhi;
    double sum = 0, tmp = 0.;
    double cosphi = TMath::Cos(phi);
    double sinphi = TMath::Sin(phi);
    for (const auto& vec : inputVectors_) {
      sum += TMath::Abs(cosphi * vec.x() + sinphi * vec.y());
    }
    tmp = TMath::Pi() / 2 * sum / area;
    if (circularity < 0 || tmp < circularity) {
      circularity = tmp;
    }
  }
  return circularity;
}

void EventShapeVariables::set_r(double r) {
  r_ = r;
  tensors_computed_ = false;
  eigenValues_ = std::vector<double>(3, 0);
  eigenValuesNoNorm_ = std::vector<double>(3, 0);
}

void EventShapeVariables::compTensorsAndVectors() {
  if (tensors_computed_)
    return;

  if (inputVectors_.size() < 2) {
    tensors_computed_ = true;
    return;
  }

  TMatrixDSym momentumTensor(3);
  momentumTensor.Zero();

  // fill momentumTensor from inputVectors
  double norm = 0.;
  for (const auto& vec : inputVectors_) {
    double p2 = vec.Dot(vec);
    double pR = (r_ == 2.) ? p2 : TMath::Power(p2, 0.5 * r_);
    norm += pR;
    double pRminus2 = (r_ == 2.) ? 1. : TMath::Power(p2, 0.5 * r_ - 1.);
    momentumTensor(0, 0) += pRminus2 * vec.x() * vec.x();
    momentumTensor(0, 1) += pRminus2 * vec.x() * vec.y();
    momentumTensor(0, 2) += pRminus2 * vec.x() * vec.z();
    momentumTensor(1, 0) += pRminus2 * vec.y() * vec.x();
    momentumTensor(1, 1) += pRminus2 * vec.y() * vec.y();
    momentumTensor(1, 2) += pRminus2 * vec.y() * vec.z();
    momentumTensor(2, 0) += pRminus2 * vec.z() * vec.x();
    momentumTensor(2, 1) += pRminus2 * vec.z() * vec.y();
    momentumTensor(2, 2) += pRminus2 * vec.z() * vec.z();
  }

  if (momentumTensor.IsSymmetric() && (momentumTensor.NonZeros() != 0)) {
    momentumTensor.EigenVectors(eigenValuesNoNormTmp_);
  }
  eigenValuesNoNorm_[0] = eigenValuesNoNormTmp_(0);
  eigenValuesNoNorm_[1] = eigenValuesNoNormTmp_(1);
  eigenValuesNoNorm_[2] = eigenValuesNoNormTmp_(2);

  // momentumTensor normalized to determinant 1
  momentumTensor *= (1. / norm);

  // now get eigens
  if (momentumTensor.IsSymmetric() && (momentumTensor.NonZeros() != 0)) {
    eigenVectors_ = momentumTensor.EigenVectors(eigenValuesTmp_);
  }
  eigenValues_[0] = eigenValuesTmp_(0);
  eigenValues_[1] = eigenValuesTmp_(1);
  eigenValues_[2] = eigenValuesTmp_(2);

  tensors_computed_ = true;
}

double EventShapeVariables::sphericity() {
  if (!tensors_computed_)
    compTensorsAndVectors();
  return 1.5 * (eigenValues_[1] + eigenValues_[2]);
}

double EventShapeVariables::aplanarity() {
  if (!tensors_computed_)
    compTensorsAndVectors();
  return 1.5 * eigenValues_[2];
}

double EventShapeVariables::C() {
  if (!tensors_computed_)
    compTensorsAndVectors();
  return 3. *
         (eigenValues_[0] * eigenValues_[1] + eigenValues_[0] * eigenValues_[2] + eigenValues_[1] * eigenValues_[2]);
}

double EventShapeVariables::D() {
  if (!tensors_computed_)
    compTensorsAndVectors();
  return 27. * eigenValues_[0] * eigenValues_[1] * eigenValues_[2];
}

//========================================================================================================

void EventShapeVariables::setFWmax(unsigned m) {
  fwmom_maxl_ = m;
  fwmom_computed_ = false;
  fwmom_ = std::vector<double>(fwmom_maxl_, 0.);
}

double EventShapeVariables::getFWmoment(unsigned l) {
  if (l > fwmom_maxl_)
    return 0.;

  if (!fwmom_computed_)
    computeFWmoments();

  return fwmom_[l];

}  // getFWmoment

const std::vector<double>& EventShapeVariables::getFWmoments() {
  if (!fwmom_computed_)
    computeFWmoments();

  return fwmom_;
}

void EventShapeVariables::computeFWmoments() {
  if (fwmom_computed_)
    return;

  double esum_total(0.);
  for (unsigned int i = 0; i < inputVectors_.size(); i++) {
    esum_total += inputVectors_[i].R();
  }  // i
  double esum_total_sq = esum_total * esum_total;

  for (unsigned int i = 0; i < inputVectors_.size(); i++) {
    double p_i = inputVectors_[i].R();
    if (p_i <= 0)
      continue;

    for (unsigned int j = 0; j <= i; j++) {
      double p_j = inputVectors_[j].R();
      if (p_j <= 0)
        continue;

      int symmetry_factor = 2;
      if (j == i)
        symmetry_factor = 1;
      double p_ij = p_i * p_j;
      double cosTheta = inputVectors_[i].Dot(inputVectors_[j]) / (p_ij);
      double pi_pj_over_etot2 = p_ij / esum_total_sq;

      double Pn1 = 0;
      double Pn2 = 0;
      for (unsigned n = 0; n < fwmom_maxl_; n++) {
        if (n == 0) {
          Pn2 = pi_pj_over_etot2;
          fwmom_[0] += Pn2 * symmetry_factor;
        } else if (n == 1) {
          Pn1 = pi_pj_over_etot2 * cosTheta;
          fwmom_[1] += Pn1 * symmetry_factor;
        } else {
          double Pn = ((2 * n - 1) * cosTheta * Pn1 - (n - 1) * Pn2) / n;
          fwmom_[n] += Pn * symmetry_factor;
          Pn2 = Pn1;
          Pn1 = Pn;
        }
      }

    }  // j
  }    // i

  fwmom_computed_ = true;

}  // computeFWmoments

//========================================================================================================