CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/PhysicsTools/IsolationUtils/src/IntegrandThetaFunction.cc

Go to the documentation of this file.
00001 #include "PhysicsTools/IsolationUtils/interface/IntegrandThetaFunction.h"
00002 
00003 // -*- C++ -*-
00004 //
00005 // Package:    IntegrandThetaFunction
00006 // Class:      IntegrandThetaFunction
00007 // 
00016 //
00017 // Original Author:  Christian Veelken, UC Davis
00018 //         Created:  Thu Nov  2 13:47:40 CST 2006
00019 // $Id: IntegrandThetaFunction.cc,v 1.3 2009/01/14 10:53:15 hegner Exp $
00020 //
00021 //
00022 
00023 // system 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/IntegralOverPhiFunction.h"
00035 #include "DataFormats/Math/interface/normalizedPhi.h"
00036 
00037 //
00038 // constructors and destructor
00039 //
00040 
00041 IntegrandThetaFunction::IntegrandThetaFunction()
00042   : ROOT::Math::ParamFunction<ROOT::Math::IParametricGradFunctionOneDim>(3)
00043 {
00044   theta0_ = 0.; 
00045   phi0_ = 0.; 
00046   alpha_ = 0.; 
00047 
00048   fPhi_ = new IntegralOverPhiFunction();
00049 }
00050 
00051 IntegrandThetaFunction::IntegrandThetaFunction(const IntegrandThetaFunction& bluePrint)
00052 {
00053   theta0_ = bluePrint.theta0_;
00054   phi0_ = bluePrint.phi0_;
00055   alpha_ = bluePrint.alpha_;
00056 
00057   fPhi_ = new IntegralOverPhiFunction(*bluePrint.fPhi_);
00058 }
00059 
00060 IntegrandThetaFunction::~IntegrandThetaFunction()
00061 {
00062   delete fPhi_;
00063 }
00064 
00065 //
00066 // assignment operator
00067 //
00068 
00069 IntegrandThetaFunction& IntegrandThetaFunction::operator=(const IntegrandThetaFunction& bluePrint)
00070 {
00071   theta0_ = bluePrint.theta0_;
00072   phi0_ = bluePrint.phi0_;
00073   alpha_ = bluePrint.alpha_;
00074   
00075   (*fPhi_) = (*bluePrint.fPhi_);
00076 
00077   return (*this);
00078 }
00079 
00080 //
00081 // member functions
00082 //
00083 
00084 void IntegrandThetaFunction::SetParameterTheta0(double theta0)
00085 {
00086   theta0_ = theta0;
00087 }
00088 
00089 void IntegrandThetaFunction::SetParameterPhi0(double phi0)
00090 {
00091   phi0_ = normalizedPhi(phi0); // map azimuth angle into interval [-pi,+pi]
00092 }
00093 
00094 void IntegrandThetaFunction::SetParameterAlpha(double alpha)
00095 {
00096   alpha_ = alpha;
00097 }
00098 
00099 void IntegrandThetaFunction::SetParameters(double* param)
00100 {
00101   theta0_ = param[0];
00102   phi0_ = param[1];
00103   alpha_ = param[2];
00104 }
00105 
00106 double IntegrandThetaFunction::DoEvalPar(double x, const double* param) const  //FIXME: constness
00107 {
00108   theta0_ = param[0];
00109   phi0_ = param[1];
00110   alpha_ = param[2];
00111 
00112   return DoEval(x);
00113 }
00114 
00115 double IntegrandThetaFunction::DoEval(double x) const
00116 {
00117 //--- return zero if theta either close  to zero or close to Pi
00118 //    (numerical expressions might become "NaN"s)
00119   const double epsilon = 1.e-3;
00120   if ( x < epsilon || x > (TMath::Pi() - epsilon) ) return 0.;
00121 
00122 //--- calculate trigonometric expressions
00123 //    (dependend on angle theta;
00124 //     polar angle of point within cone)
00125   double sinTheta = TMath::Sin(x);
00126   double cscTheta = 1./sinTheta;
00127 
00128   double detJacobi = -cscTheta; // partial derrivative dEta/dTheta (for constant particle density in tau id. isolation cone)
00129   //double detJacobi = 1.; // ordinary solid angle (FOR TESTING ONLY)
00130 
00131 //--- evaluate integral over angle phi
00132 //    (azimuth angle of point within cone)
00133   fPhi_->SetParameterTheta0(theta0_);
00134   fPhi_->SetParameterPhi0(phi0_);
00135   fPhi_->SetParameterAlpha(alpha_);
00136 
00137   double integralOverPhi = (*fPhi_)(x);
00138 
00139   if ( debugLevel_ > 0 ) {
00140     edm::LogVerbatim("") << "integralOverPhi = " << integralOverPhi << std::endl
00141                          << " theta0 = " << theta0_ << std::endl
00142                          << " phi0 = " << phi0_ << std::endl
00143                          << " alpha = " << alpha_ << std::endl
00144                          << " theta = " << x << std::endl
00145                          << std::endl;
00146   }
00147   
00148 //--- integrand for integration over theta
00149 //    equals 
00150 //      |dEta/dTheta| * integral over phi * sin(theta)
00151 //
00152 //    (o the factor dEta/dTheta represents the particle density as function of theta, 
00153 //       assuming that the particle density is flat in eta;
00154 //     o the factor sin(theta) originates from the solid angle surface element 
00155 //       expressed in spherical polar coordinates)
00156 //
00157   //return TMath::Abs(detJacobi)*integralOverPhi*sinTheta;
00158   return TMath::Abs(detJacobi)*integralOverPhi;
00159 }  
00160 
00161 double IntegrandThetaFunction::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 IntegrandThetaFunction::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 IntegrandThetaFunction::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