00001 #include "PhysicsTools/IsolationUtils/interface/ConeAreaFunction.h" 00002 00003 // -*- C++ -*- 00004 // 00005 // Package: ConeAreaFunction 00006 // Class: ConeAreaFunction 00007 // 00016 // 00017 // Original Author: Christian Veelken, UC Davis 00018 // Created: Thu Nov 2 13:47:40 CST 2006 00019 // $Id: ConeAreaFunction.cc,v 1.3 2009/01/14 10:53:15 hegner Exp $ 00020 // 00021 // 00022 00023 // C++ standard library include files 00024 #include <iostream> 00025 #include <iomanip> 00026 #include <vector> 00027 00028 // ROOT include files 00029 #include <TMath.h> 00030 00031 // CMSSW include files 00032 #include "FWCore/MessageLogger/interface/MessageLogger.h" 00033 00034 #include "PhysicsTools/IsolationUtils/interface/IntegrandThetaFunction.h" 00035 #include "DataFormats/Math/interface/normalizedPhi.h" 00036 00037 // 00038 // constructors and destructor 00039 // 00040 00041 ConeAreaFunction::ConeAreaFunction() 00042 : ROOT::Math::ParamFunction<ROOT::Math::IParametricGradFunctionOneDim>(2) 00043 { 00044 theta0_ = 0.; 00045 phi0_ = 0.; 00046 00047 etaMax_ = -1; 00048 00049 fTheta_ = new IntegrandThetaFunction(); 00050 integrator_ = new ROOT::Math::Integrator(*fTheta_); 00051 } 00052 00053 ConeAreaFunction::ConeAreaFunction(const ConeAreaFunction& bluePrint) 00054 { 00055 theta0_ = bluePrint.theta0_; 00056 phi0_ = bluePrint.phi0_; 00057 00058 etaMax_ = bluePrint.etaMax_; 00059 00060 fTheta_ = new IntegrandThetaFunction(*bluePrint.fTheta_); 00061 integrator_ = new ROOT::Math::Integrator(*fTheta_); 00062 } 00063 00064 ConeAreaFunction::~ConeAreaFunction() 00065 { 00066 delete fTheta_; // function gets deleted automatically by Integrator ? 00067 delete integrator_; 00068 } 00069 00070 // 00071 // assignment operator 00072 // 00073 00074 ConeAreaFunction& ConeAreaFunction::operator=(const ConeAreaFunction& bluePrint) 00075 { 00076 theta0_ = bluePrint.theta0_; 00077 phi0_ = bluePrint.phi0_; 00078 00079 etaMax_ = bluePrint.etaMax_; 00080 00081 (*fTheta_) = (*bluePrint.fTheta_); 00082 integrator_->SetFunction(*fTheta_); 00083 00084 return (*this); 00085 } 00086 00087 // 00088 // member functions 00089 // 00090 00091 void ConeAreaFunction::SetParameterTheta0(double theta0) 00092 { 00093 theta0_ = theta0; 00094 } 00095 00096 void ConeAreaFunction::SetParameterPhi0(double phi0) 00097 { 00098 phi0_ = normalizedPhi(phi0); // map azimuth angle into interval [-pi,+pi] 00099 } 00100 00101 void ConeAreaFunction::SetParameters(double* param) 00102 { 00103 if ( debugLevel_ > 0 ) { 00104 edm::LogVerbatim("") << "<ConeAreaFunction::SetParameters>:" << std::endl 00105 << " theta0 = " << param[0] << std::endl 00106 << " phi0 = " << param[1] << std::endl; 00107 } 00108 00109 theta0_ = param[0]; 00110 phi0_ = param[1]; 00111 } 00112 00113 void ConeAreaFunction::SetAcceptanceLimit(double etaMax) 00114 { 00115 //--- check that pseudo-rapidity given as function argument is positive 00116 // (assume equal acceptance for positive and negative pseudo-rapidities) 00117 00118 if ( etaMax > 0 ) { 00119 etaMax_ = etaMax; 00120 } else { 00121 edm::LogError("") << "etaMax cannot be negative !" << std::endl; 00122 } 00123 } 00124 00125 double ConeAreaFunction::DoEvalPar(double x, const double* param) const 00126 { 00127 //--- calculate area covered by cone of opening angle alpha 00128 // (measured from cone axis); 00129 // evaluate integral over angle theta 00130 // (polar angle of point within cone) 00131 // FIXME: the const above is actually not true as it is implemented now. 00132 00133 theta0_ = param[0]; 00134 phi0_ = param[1]; 00135 00136 return DoEval(x); 00137 } 00138 00139 00140 double ConeAreaFunction::DoEval(double x) const 00141 { 00142 //--- calculate area covered by cone of opening angle alpha 00143 // (measured from cone axis); 00144 // evaluate integral over angle theta 00145 // (polar angle of point within cone) 00146 00147 fTheta_->SetParameterTheta0(theta0_); 00148 fTheta_->SetParameterPhi0(phi0_); 00149 fTheta_->SetParameterAlpha(x); 00150 00151 integrator_->SetFunction(*fTheta_); // set updated parameter values in Integrator 00152 00153 double thetaMin = (etaMax_ > 0) ? 2*TMath::ATan(TMath::Exp(-etaMax_)) : 0.; 00154 double thetaMax = TMath::Pi() - thetaMin; 00155 00156 double integralOverTheta = integrator_->Integral(thetaMin, thetaMax); 00157 00158 return integralOverTheta; 00159 } 00160 00161 double ConeAreaFunction::DoDerivative(double x) const 00162 { 00163 //--- virtual function inherited from ROOT::Math::ParamFunction base class; 00164 // not implemented, because not neccessary, but needs to be defined to make code compile... 00165 edm::LogWarning("") << "Function not implemented yet !" << std::endl; 00166 00167 return 0.; 00168 } 00169 00170 double ConeAreaFunction::DoParameterDerivative(double, const double*, unsigned int) const 00171 { 00172 //--- virtual function inherited from ROOT::Math::ParamFunction base class; 00173 // not implemented, because not neccessary, but needs to be defined to make code compile... 00174 edm::LogWarning("") << "Function not implemented yet !" << std::endl; 00175 00176 return 0.; 00177 } 00178 00179 00180 00181 void ConeAreaFunction::DoParameterGradient(double x, double* paramGradient) const 00182 { 00183 //--- virtual function inherited from ROOT::Math::ParamFunction base class; 00184 // not implemented, because not neccessary, but needs to be defined to make code compile... 00185 edm::LogWarning("") << "Function not implemented yet !" << std::endl; 00186 } 00187