IntegralOverPhiFunction.cc

Go to the documentation of this file.

#include "PhysicsTools/IsolationUtils/interface/IntegralOverPhiFunction.h"

// -*- C++ -*-
//
// Package:    IntegralOverPhiFunction
// Class:      IntegralOverPhiFunction
//
//
// Original Author:  Christian Veelken, UC Davis
//         Created:  Thu Nov  2 13:47:40 CST 2006
//
//

// system include files
#include <iostream>
#include <iomanip>
#include <vector>

// ROOT include files
#include <TMath.h>

// CMSSW include files
#include "FWCore/MessageLogger/interface/MessageLogger.h"

#include "DataFormats/Math/interface/normalizedPhi.h"

//
// declaration of auxiliary functions
//

void checkSolutions(unsigned int i,
                    unsigned int& numSolution1,
                    unsigned int& numSolution2,
                    unsigned int& numSolution3,
                    unsigned int& numSolution4);

//
// constructors and destructor
//

IntegralOverPhiFunction::IntegralOverPhiFunction()
    : ROOT::Math::ParamFunction<ROOT::Math::IParametricGradFunctionOneDim>(3) {
  theta0_ = 0.;
  phi0_ = 0.;
  alpha_ = 0.;

  // !!! ONLY FOR TESTING
  numSolutionMin1_ = 0;
  numSolutionMax1_ = 0;
  numSolutionMin2_ = 0;
  numSolutionMax2_ = 0;
  numSolutionMin3_ = 0;
  numSolutionMax3_ = 0;
  numSolutionMin4_ = 0;
  numSolutionMax4_ = 0;
  //     FOR TESTING ONLY !!!
}

IntegralOverPhiFunction::~IntegralOverPhiFunction() {
  // !!! ONLY FOR TESTING
  if (debugLevel_ > 0) {
    edm::LogVerbatim("") << "<IntegralOverPhiFunction::~IntegralOverPhiFunction>:" << std::endl
                         << " numSolutionMin1 = " << numSolutionMin1_ << std::endl
                         << " numSolutionMax1 = " << numSolutionMax1_ << std::endl
                         << " numSolutionMin2 = " << numSolutionMin2_ << std::endl
                         << " numSolutionMax2 = " << numSolutionMax2_ << std::endl
                         << " numSolutionMin3 = " << numSolutionMin3_ << std::endl
                         << " numSolutionMax3 = " << numSolutionMax3_ << std::endl
                         << " numSolutionMin4 = " << numSolutionMin4_ << std::endl
                         << " numSolutionMax4 = " << numSolutionMax4_ << std::endl;
  }
  //     FOR TESTING ONLY !!!
}

//
// member functions
//

void IntegralOverPhiFunction::SetParameterTheta0(double theta0) { theta0_ = theta0; }

void IntegralOverPhiFunction::SetParameterPhi0(double phi0) {
  phi0_ = normalizedPhi(phi0);  // map azimuth angle into interval [-pi,+pi]
}

void IntegralOverPhiFunction::SetParameterAlpha(double alpha) { alpha_ = alpha; }

void IntegralOverPhiFunction::SetParameters(const double* param) {
  theta0_ = param[0];
  phi0_ = param[1];
  alpha_ = param[2];
}

double IntegralOverPhiFunction::DoEvalPar(
    double x, const double* param) const  //FIXME: in the current implementation const is not entirely true
{
  theta0_ = param[0];
  phi0_ = param[1];
  alpha_ = param[2];

  return DoEval(x);
}

double IntegralOverPhiFunction::DoEval(double x) const {
  //--- return zero if theta either close to zero or close to Pi
  //    (phi not well-defined in that case;
  //     numerical expressions might become "NaN"s)
  const double epsilon = 1.e-3;
  if (x < epsilon || x > (TMath::Pi() - epsilon))
    return 0.;

  //--- calculate trigonometric expressions
  //    (dependend on angle theta0;
  //     polar angle of cone axis);
  //    avoid theta0 exactly zero or exactly Pi
  //    (numerical expressions become "NaN"s)
  double theta0 = theta0_;
  if (theta0 < epsilon)
    theta0 = epsilon;
  if (theta0 > (TMath::Pi() - epsilon))
    theta0 = (TMath::Pi() - epsilon);
  double cosTheta0 = TMath::Cos(theta0);
  double cos2Theta0 = TMath::Cos(2 * theta0);
  double sinTheta0 = TMath::Sin(theta0);
  double cotTheta0 = 1. / TMath::Tan(theta0);
  double cscTheta0 = 1. / TMath::Sin(theta0);
  //    (dependend on angle phi0;
  //     azimuth angle of cone axis)
  double cosPhi0 = TMath::Cos(phi0_);
  double tanPhi0 = TMath::Tan(phi0_);
  //    (dependend on angle theta;
  //     polar angle of point within cone)
  double cosTheta = TMath::Cos(x);
  double cos2Theta = TMath::Cos(2 * x);
  double sinTheta = TMath::Sin(x);
  double cotTheta = 1. / TMath::Tan(x);
  double cscTheta = 1. / TMath::Sin(x);
  //    (dependent on angle alpha;
  //     opening angle of cone, measured from cone axis)
  double cosAlpha = TMath::Cos(alpha_);

  double s = -cosPhi0 * cosPhi0 *
             (2 * cosAlpha * cosAlpha + cos2Theta - 4 * cosAlpha * cosTheta * cosTheta0 + cos2Theta0) * sinTheta *
             sinTheta * sinTheta0 * sinTheta0;

  //--- return either zero or 2 Pi
  //    in case argument of square-root is zero or negative
  //    (negative values solely arise from rounding errors);
  //    check whether to return zero or 2 Pi:
  //     o return zero
  //       if |theta- theta0| > alpha,
  //       (theta outside cone of opening angle alpha, so vanishing integral over phi)
  //     o return 2 Pi
  //       if |theta- theta0| < alpha
  //       (theta within cone of opening angle alpha;
  //        actually theta0 < alpha in forward/backward direction,
  //        so that integral includes all phi values, hence yielding 2 Pi)
  if (s <= 0) {
    if (TMath::Abs(x - theta0) > alpha_)
      return 0;
    if (TMath::Abs(x - theta0) <= alpha_)
      return 2 * TMath::Pi();
    std::cerr << "Error in <IntegralOverPhiFunction::operator()>: failed to compute return value !" << std::endl;
  }

  double r = (1. / TMath::Sqrt(2.)) * (cscTheta * cscTheta * cscTheta0 * cscTheta0 * TMath::Sqrt(s) * tanPhi0);
  double t = cosPhi0 * (-cotTheta * cotTheta0 + cosAlpha * cscTheta * cscTheta0);

  double phi[4];
  phi[0] = -TMath::ACos(t - r);
  phi[1] = TMath::ACos(t - r);
  phi[2] = -TMath::ACos(t + r);
  phi[3] = TMath::ACos(t + r);

  if (debugLevel_ > 0) {
    edm::LogVerbatim("") << "phi0 = " << phi0_ << std::endl
                         << "phi[0] = " << phi[0] << " (phi[0] - phi0 = " << (phi[0] - phi0_) * 180 / TMath::Pi() << ")"
                         << std::endl
                         << "phi[1] = " << phi[1] << " (phi[1] - phi0 = " << (phi[1] - phi0_) * 180 / TMath::Pi() << ")"
                         << std::endl
                         << "phi[2] = " << phi[2] << " (phi[2] - phi0 = " << (phi[2] - phi0_) * 180 / TMath::Pi() << ")"
                         << std::endl
                         << "phi[3] = " << phi[3] << " (phi[3] - phi0 = " << (phi[3] - phi0_) * 180 / TMath::Pi() << ")"
                         << std::endl;
  }

  double phiMin = 0.;
  double phiMax = 0.;
  for (unsigned int i = 0; i < 4; ++i) {
    for (unsigned int j = (i + 1); j < 4; ++j) {
      //--- search for the two solutions for phi
      //    that have an equal distance to phi0
      //    and are on either side
      double dPhi_i = phi[i] - phi0_;
      double dPhi_j = phi0_ - phi[j];
      if (TMath::Abs(normalizedPhi(dPhi_i - dPhi_j)) <
          epsilon) {  // map difference in azimuth angle into interval [-pi,+pi] and require it to be negligible
        //if ( TMath::Abs((phi[i] - phi0_) - (phi0_ - phi[j])) < epsilon ) { // map difference in azimuth angle into interval [-pi,+pi] and require it to be negligible
        phiMin = TMath::Min(phi[i], phi[j]);
        phiMax = TMath::Max(phi[i], phi[j]);

        // !!! ONLY FOR TESTING
        if (phi[i] == phiMin)
          checkSolutions(i, numSolutionMin1_, numSolutionMin2_, numSolutionMin3_, numSolutionMin4_);
        if (phi[i] == phiMax)
          checkSolutions(i, numSolutionMax1_, numSolutionMax2_, numSolutionMax3_, numSolutionMax4_);
        if (phi[j] == phiMin)
          checkSolutions(j, numSolutionMin1_, numSolutionMin2_, numSolutionMin3_, numSolutionMin4_);
        if (phi[j] == phiMax)
          checkSolutions(j, numSolutionMax1_, numSolutionMax2_, numSolutionMax3_, numSolutionMax4_);
        //     FOR TESTING ONLY !!!
      }
    }
  }

  if (phiMin == 0 && phiMax == 0) {
    edm::LogError("") << "failed to compute Return Value !" << std::endl;
  }

  return TMath::Abs(normalizedPhi(phi0_ - phiMin)) + TMath::Abs(normalizedPhi(phiMax - phi0_));
}

double IntegralOverPhiFunction::DoDerivative(double x) const {
  //--- virtual function inherited from ROOT::Math::ParamFunction base class;
  //    not implemented, because not neccessary, but needs to be defined to make code compile...
  edm::LogWarning("") << "Function not implemented yet !" << std::endl;

  return 0.;
}

double IntegralOverPhiFunction::DoParameterDerivative(double, const double*, unsigned int) const {
  //--- virtual function inherited from ROOT::Math::ParamFunction base class;
  //    not implemented, because not neccessary, but needs to be defined to make code compile...
  edm::LogWarning("") << "Function not implemented yet !" << std::endl;

  return 0.;
}

void IntegralOverPhiFunction::DoParameterGradient(double x, double* paramGradient) const {
  //--- virtual function inherited from ROOT::Math::ParamFunction base class;
  //    not implemented, because not neccessary, but needs to be defined to make code compile...
  edm::LogWarning("") << "Function not implemented yet !" << std::endl;
}

//
// definition of auxiliary functions
//

void checkSolutions(unsigned int i,
                    unsigned int& numSolution1,
                    unsigned int& numSolution2,
                    unsigned int& numSolution3,
                    unsigned int& numSolution4) {
  switch (i) {
    case 0:
      ++numSolution1;
      break;
    case 1:
      ++numSolution2;
      break;
    case 2:
      ++numSolution3;
      break;
    case 3:
      ++numSolution4;
      break;
  }
}