IsolationUtils/src/IntegrandThetaFunction.cc

0001 #include "PhysicsTools/IsolationUtils/interface/IntegrandThetaFunction.h"
0002
0003 // -*- C++ -*-
0004 //
0005 // Package:    IntegrandThetaFunction
0006 // Class:      IntegrandThetaFunction
0007 //
0008 /**\class IntegrandThetaFunction IntegrandThetaFunction.cc PhysicsTools/IsolationUtils/src/IntegrandThetaFunction.cc
0009
0010  Description: auxialiary class for fixed area isolation cone computation
0011               (this class performs the integration over the polar angle)
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 // system 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/IntegralOverPhiFunction.h"
0034 #include "DataFormats/Math/interface/normalizedPhi.h"
0035
0036 //
0037 // constructors and destructor
0038 //
0039
0040 IntegrandThetaFunction::IntegrandThetaFunction()
0041     : ROOT::Math::ParamFunction<ROOT::Math::IParametricGradFunctionOneDim>(3) {
0042   theta0_ = 0.;
0043   phi0_ = 0.;
0044   alpha_ = 0.;
0045
0046   fPhi_ = new IntegralOverPhiFunction();
0047 }
0048
0049 IntegrandThetaFunction::IntegrandThetaFunction(const IntegrandThetaFunction& bluePrint)
0050     : ROOT::Math::ParamFunction<ROOT::Math::IParametricGradFunctionOneDim>(bluePrint) {
0051   theta0_ = bluePrint.theta0_;
0052   phi0_ = bluePrint.phi0_;
0053   alpha_ = bluePrint.alpha_;
0054
0055   fPhi_ = new IntegralOverPhiFunction(*bluePrint.fPhi_);
0056 }
0057
0058 IntegrandThetaFunction::~IntegrandThetaFunction() { delete fPhi_; }
0059
0060 //
0061 // assignment operator
0062 //
0063
0064 IntegrandThetaFunction& IntegrandThetaFunction::operator=(const IntegrandThetaFunction& bluePrint) {
0065   theta0_ = bluePrint.theta0_;
0066   phi0_ = bluePrint.phi0_;
0067   alpha_ = bluePrint.alpha_;
0068
0069   (*fPhi_) = (*bluePrint.fPhi_);
0070
0071   return (*this);
0072 }
0073
0074 //
0075 // member functions
0076 //
0077
0078 void IntegrandThetaFunction::SetParameterTheta0(double theta0) { theta0_ = theta0; }
0079
0080 void IntegrandThetaFunction::SetParameterPhi0(double phi0) {
0081   phi0_ = normalizedPhi(phi0);  // map azimuth angle into interval [-pi,+pi]
0082 }
0083
0084 void IntegrandThetaFunction::SetParameterAlpha(double alpha) { alpha_ = alpha; }
0085
0086 void IntegrandThetaFunction::SetParameters(const double* param) {
0087   theta0_ = param[0];
0088   phi0_ = param[1];
0089   alpha_ = param[2];
0090 }
0091
0092 double IntegrandThetaFunction::DoEvalPar(double x, const double* param) const  //FIXME: constness
0093 {
0094   theta0_ = param[0];
0095   phi0_ = param[1];
0096   alpha_ = param[2];
0097
0098   return DoEval(x);
0099 }
0100
0101 double IntegrandThetaFunction::DoEval(double x) const {
0102   //--- return zero if theta either close  to zero or close to Pi
0103   //    (numerical expressions might become "NaN"s)
0104   const double epsilon = 1.e-3;
0105   if (x < epsilon || x > (TMath::Pi() - epsilon))
0106     return 0.;
0107
0108   //--- calculate trigonometric expressions
0109   //    (dependend on angle theta;
0110   //     polar angle of point within cone)
0111   double sinTheta = TMath::Sin(x);
0112   double cscTheta = 1. / sinTheta;
0113
0114   double detJacobi =
0115       -cscTheta;  // partial derrivative dEta/dTheta (for constant particle density in tau id. isolation cone)
0116   //double detJacobi = 1.; // ordinary solid angle (FOR TESTING ONLY)
0117
0118   //--- evaluate integral over angle phi
0119   //    (azimuth angle of point within cone)
0120   fPhi_->SetParameterTheta0(theta0_);
0121   fPhi_->SetParameterPhi0(phi0_);
0122   fPhi_->SetParameterAlpha(alpha_);
0123
0124   double integralOverPhi = (*fPhi_)(x);
0125
0126   if (debugLevel_ > 0) {
0127     edm::LogVerbatim("") << "integralOverPhi = " << integralOverPhi << std::endl
0128                          << " theta0 = " << theta0_ << std::endl
0129                          << " phi0 = " << phi0_ << std::endl
0130                          << " alpha = " << alpha_ << std::endl
0131                          << " theta = " << x << std::endl
0132                          << std::endl;
0133   }
0134
0135   //--- integrand for integration over theta
0136   //    equals
0137   //      |dEta/dTheta| * integral over phi * sin(theta)
0138   //
0139   //    (o the factor dEta/dTheta represents the particle density as function of theta,
0140   //       assuming that the particle density is flat in eta;
0141   //     o the factor sin(theta) originates from the solid angle surface element
0142   //       expressed in spherical polar coordinates)
0143   //
0144   //return TMath::Abs(detJacobi)*integralOverPhi*sinTheta;
0145   return TMath::Abs(detJacobi) * integralOverPhi;
0146 }
0147
0148 double IntegrandThetaFunction::DoDerivative(double x) const {
0149   //--- virtual function inherited from ROOT::Math::ParamFunction base class;
0150   //    not implemented, because not neccessary, but needs to be defined to make code compile...
0151   edm::LogWarning("") << "Function not implemented yet !" << std::endl;
0152
0153   return 0.;
0154 }
0155
0156 double IntegrandThetaFunction::DoParameterDerivative(double, const double*, unsigned int) const {
0157   //--- virtual function inherited from ROOT::Math::ParamFunction base class;
0158   //    not implemented, because not neccessary, but needs to be defined to make code compile...
0159   edm::LogWarning("") << "Function not implemented yet !" << std::endl;
0160
0161   return 0.;
0162 }
0163
0164 void IntegrandThetaFunction::DoParameterGradient(double x, double* paramGradient) const {
0165   //--- virtual function inherited from ROOT::Math::ParamFunction base class;
0166   //    not implemented, because not neccessary, but needs to be defined to make code compile...
0167   edm::LogWarning("") << "Function not implemented yet !" << std::endl;
0168 }