IsolationUtils/src/ConeAreaFunction.cc

0001 #include "PhysicsTools/IsolationUtils/interface/ConeAreaFunction.h"
0002
0003 // -*- C++ -*-
0004 //
0005 // Package:    ConeAreaFunction
0006 // Class:      ConeAreaFunction
0007 //
0008 /**\class ConeAreaFunction ConeAreaFunction.cc PhysicsTools/IsolationUtils//src/ConeAreaFunction.cc
0009
0010  Description: low level class to compute area of signal cone
0011               corresponding to three-dimensional opening angle alpha given as function argument
0012
0013  Implementation:
0014      imported into CMSSW on 05/18/2007
0015 */
0016 //
0017 // Original Author:  Christian Veelken, UC Davis
0018 //         Created:  Thu Nov  2 13:47:40 CST 2006
0019 //
0020 //
0021
0022 // C++ standard library include files
0023 #include <iostream>
0024 #include <iomanip>
0025 #include <vector>
0026
0027 // ROOT include files
0028 #include <TMath.h>
0029
0030 // CMSSW include files
0031 #include "FWCore/MessageLogger/interface/MessageLogger.h"
0032
0033 #include "PhysicsTools/IsolationUtils/interface/IntegrandThetaFunction.h"
0034 #include "DataFormats/Math/interface/normalizedPhi.h"
0035
0036 //
0037 // constructors and destructor
0038 //
0039
0040 ConeAreaFunction::ConeAreaFunction() : ROOT::Math::ParamFunction<ROOT::Math::IParametricGradFunctionOneDim>(2) {
0041   theta0_ = 0.;
0042   phi0_ = 0.;
0043
0044   etaMax_ = -1;
0045
0046   fTheta_ = new IntegrandThetaFunction();
0047   integrator_ = new ROOT::Math::Integrator(*fTheta_);
0048 }
0049
0050 ConeAreaFunction::ConeAreaFunction(const ConeAreaFunction& bluePrint)
0051     : ROOT::Math::ParamFunction<ROOT::Math::IParametricGradFunctionOneDim>(bluePrint) {
0052   theta0_ = bluePrint.theta0_;
0053   phi0_ = bluePrint.phi0_;
0054
0055   etaMax_ = bluePrint.etaMax_;
0056
0057   fTheta_ = new IntegrandThetaFunction(*bluePrint.fTheta_);
0058   integrator_ = new ROOT::Math::Integrator(*fTheta_);
0059 }
0060
0061 ConeAreaFunction::~ConeAreaFunction() {
0062   delete fTheta_;  // function gets deleted automatically by Integrator ?
0063   delete integrator_;
0064 }
0065
0066 //
0067 // assignment operator
0068 //
0069
0070 ConeAreaFunction& ConeAreaFunction::operator=(const ConeAreaFunction& bluePrint) {
0071   theta0_ = bluePrint.theta0_;
0072   phi0_ = bluePrint.phi0_;
0073
0074   etaMax_ = bluePrint.etaMax_;
0075
0076   (*fTheta_) = (*bluePrint.fTheta_);
0077   integrator_->SetFunction(*fTheta_);
0078
0079   return (*this);
0080 }
0081
0082 //
0083 // member functions
0084 //
0085
0086 void ConeAreaFunction::SetParameterTheta0(double theta0) { theta0_ = theta0; }
0087
0088 void ConeAreaFunction::SetParameterPhi0(double phi0) {
0089   phi0_ = normalizedPhi(phi0);  // map azimuth angle into interval [-pi,+pi]
0090 }
0091
0092 void ConeAreaFunction::SetParameters(const double* param) {
0093   if (debugLevel_ > 0) {
0094     edm::LogVerbatim("") << "<ConeAreaFunction::SetParameters>:" << std::endl
0095                          << " theta0 = " << param[0] << std::endl
0096                          << " phi0 = " << param[1] << std::endl;
0097   }
0098
0099   theta0_ = param[0];
0100   phi0_ = param[1];
0101 }
0102
0103 void ConeAreaFunction::SetAcceptanceLimit(double etaMax) {
0104   //--- check that pseudo-rapidity given as function argument is positive
0105   //    (assume equal acceptance for positive and negative pseudo-rapidities)
0106
0107   if (etaMax > 0) {
0108     etaMax_ = etaMax;
0109   } else {
0110     edm::LogError("") << "etaMax cannot be negative !" << std::endl;
0111   }
0112 }
0113
0114 double ConeAreaFunction::DoEvalPar(double x, const double* param) const {
0115   //--- calculate area covered by cone of opening angle alpha
0116   //    (measured from cone axis);
0117   //    evaluate integral over angle theta
0118   //    (polar angle of point within cone)
0119   // FIXME: the const above is actually not true as it is implemented now.
0120
0121   theta0_ = param[0];
0122   phi0_ = param[1];
0123
0124   return DoEval(x);
0125 }
0126
0127 double ConeAreaFunction::DoEval(double x) const {
0128   //--- calculate area covered by cone of opening angle alpha
0129   //    (measured from cone axis);
0130   //    evaluate integral over angle theta
0131   //    (polar angle of point within cone)
0132
0133   fTheta_->SetParameterTheta0(theta0_);
0134   fTheta_->SetParameterPhi0(phi0_);
0135   fTheta_->SetParameterAlpha(x);
0136
0137   integrator_->SetFunction(*fTheta_);  // set updated parameter values in Integrator
0138
0139   double thetaMin = (etaMax_ > 0) ? 2 * TMath::ATan(TMath::Exp(-etaMax_)) : 0.;
0140   double thetaMax = TMath::Pi() - thetaMin;
0141
0142   double integralOverTheta = integrator_->Integral(thetaMin, thetaMax);
0143
0144   return integralOverTheta;
0145 }
0146
0147 double ConeAreaFunction::DoDerivative(double x) const {
0148   //--- virtual function inherited from ROOT::Math::ParamFunction base class;
0149   //    not implemented, because not neccessary, but needs to be defined to make code compile...
0150   edm::LogWarning("") << "Function not implemented yet !" << std::endl;
0151
0152   return 0.;
0153 }
0154
0155 double ConeAreaFunction::DoParameterDerivative(double, const double*, unsigned int) const {
0156   //--- virtual function inherited from ROOT::Math::ParamFunction base class;
0157   //    not implemented, because not neccessary, but needs to be defined to make code compile...
0158   edm::LogWarning("") << "Function not implemented yet !" << std::endl;
0159
0160   return 0.;
0161 }
0162
0163 void ConeAreaFunction::DoParameterGradient(double x, double* paramGradient) const {
0164   //--- virtual function inherited from ROOT::Math::ParamFunction base class;
0165   //    not implemented, because not neccessary, but needs to be defined to make code compile...
0166   edm::LogWarning("") << "Function not implemented yet !" << std::endl;
0167 }