Back to home page

Project CMSSW displayed by LXR
 
 

    

File indexing completed on 2024-04-06 12:23:32 

0001 #include "PhysicsTools/IsolationUtils/interface/IntegralOverPhiFunction.h"
0002 
0003 // -*- C++ -*-
0004 //
0005 // Package:    IntegralOverPhiFunction
0006 // Class:      IntegralOverPhiFunction
0007 //
0008 /**\class IntegralOverPhiFunction IntegralOverPhiFunction.cc PhysicsTools/IsolationUtils/src/IntegralOverPhiFunction.cc
0009 
0010  Description: auxialiary class for fixed area isolation cone computation
0011               (this class performs the integration over the azimuthal 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 "DataFormats/Math/interface/normalizedPhi.h"
0034 
0035 //
0036 // declaration of auxiliary functions
0037 //
0038 
0039 void checkSolutions(unsigned int i,
0040                     unsigned int& numSolution1,
0041                     unsigned int& numSolution2,
0042                     unsigned int& numSolution3,
0043                     unsigned int& numSolution4);
0044 
0045 //
0046 // constructors and destructor
0047 //
0048 
0049 IntegralOverPhiFunction::IntegralOverPhiFunction()
0050     : ROOT::Math::ParamFunction<ROOT::Math::IParametricGradFunctionOneDim>(3) {
0051   theta0_ = 0.;
0052   phi0_ = 0.;
0053   alpha_ = 0.;
0054 
0055   // !!! ONLY FOR TESTING
0056   numSolutionMin1_ = 0;
0057   numSolutionMax1_ = 0;
0058   numSolutionMin2_ = 0;
0059   numSolutionMax2_ = 0;
0060   numSolutionMin3_ = 0;
0061   numSolutionMax3_ = 0;
0062   numSolutionMin4_ = 0;
0063   numSolutionMax4_ = 0;
0064   //     FOR TESTING ONLY !!!
0065 }
0066 
0067 IntegralOverPhiFunction::~IntegralOverPhiFunction() {
0068   // !!! ONLY FOR TESTING
0069   if (debugLevel_ > 0) {
0070     edm::LogVerbatim("") << "<IntegralOverPhiFunction::~IntegralOverPhiFunction>:" << std::endl
0071                          << " numSolutionMin1 = " << numSolutionMin1_ << std::endl
0072                          << " numSolutionMax1 = " << numSolutionMax1_ << std::endl
0073                          << " numSolutionMin2 = " << numSolutionMin2_ << std::endl
0074                          << " numSolutionMax2 = " << numSolutionMax2_ << std::endl
0075                          << " numSolutionMin3 = " << numSolutionMin3_ << std::endl
0076                          << " numSolutionMax3 = " << numSolutionMax3_ << std::endl
0077                          << " numSolutionMin4 = " << numSolutionMin4_ << std::endl
0078                          << " numSolutionMax4 = " << numSolutionMax4_ << std::endl;
0079   }
0080   //     FOR TESTING ONLY !!!
0081 }
0082 
0083 //
0084 // member functions
0085 //
0086 
0087 void IntegralOverPhiFunction::SetParameterTheta0(double theta0) { theta0_ = theta0; }
0088 
0089 void IntegralOverPhiFunction::SetParameterPhi0(double phi0) {
0090   phi0_ = normalizedPhi(phi0);  // map azimuth angle into interval [-pi,+pi]
0091 }
0092 
0093 void IntegralOverPhiFunction::SetParameterAlpha(double alpha) { alpha_ = alpha; }
0094 
0095 void IntegralOverPhiFunction::SetParameters(const double* param) {
0096   theta0_ = param[0];
0097   phi0_ = param[1];
0098   alpha_ = param[2];
0099 }
0100 
0101 double IntegralOverPhiFunction::DoEvalPar(
0102     double x, const double* param) const  //FIXME: in the current implementation const is not entirely true
0103 {
0104   theta0_ = param[0];
0105   phi0_ = param[1];
0106   alpha_ = param[2];
0107 
0108   return DoEval(x);
0109 }
0110 
0111 double IntegralOverPhiFunction::DoEval(double x) const {
0112   //--- return zero if theta either close to zero or close to Pi
0113   //    (phi not well-defined in that case;
0114   //     numerical expressions might become "NaN"s)
0115   const double epsilon = 1.e-3;
0116   if (x < epsilon || x > (TMath::Pi() - epsilon))
0117     return 0.;
0118 
0119   //--- calculate trigonometric expressions
0120   //    (dependend on angle theta0;
0121   //     polar angle of cone axis);
0122   //    avoid theta0 exactly zero or exactly Pi
0123   //    (numerical expressions become "NaN"s)
0124   double theta0 = theta0_;
0125   if (theta0 < epsilon)
0126     theta0 = epsilon;
0127   if (theta0 > (TMath::Pi() - epsilon))
0128     theta0 = (TMath::Pi() - epsilon);
0129   double cosTheta0 = TMath::Cos(theta0);
0130   double cos2Theta0 = TMath::Cos(2 * theta0);
0131   double sinTheta0 = TMath::Sin(theta0);
0132   double cotTheta0 = 1. / TMath::Tan(theta0);
0133   double cscTheta0 = 1. / TMath::Sin(theta0);
0134   //    (dependend on angle phi0;
0135   //     azimuth angle of cone axis)
0136   double cosPhi0 = TMath::Cos(phi0_);
0137   double tanPhi0 = TMath::Tan(phi0_);
0138   //    (dependend on angle theta;
0139   //     polar angle of point within cone)
0140   double cosTheta = TMath::Cos(x);
0141   double cos2Theta = TMath::Cos(2 * x);
0142   double sinTheta = TMath::Sin(x);
0143   double cotTheta = 1. / TMath::Tan(x);
0144   double cscTheta = 1. / TMath::Sin(x);
0145   //    (dependent on angle alpha;
0146   //     opening angle of cone, measured from cone axis)
0147   double cosAlpha = TMath::Cos(alpha_);
0148 
0149   double s = -cosPhi0 * cosPhi0 *
0150              (2 * cosAlpha * cosAlpha + cos2Theta - 4 * cosAlpha * cosTheta * cosTheta0 + cos2Theta0) * sinTheta *
0151              sinTheta * sinTheta0 * sinTheta0;
0152 
0153   //--- return either zero or 2 Pi
0154   //    in case argument of square-root is zero or negative
0155   //    (negative values solely arise from rounding errors);
0156   //    check whether to return zero or 2 Pi:
0157   //     o return zero
0158   //       if |theta- theta0| > alpha,
0159   //       (theta outside cone of opening angle alpha, so vanishing integral over phi)
0160   //     o return 2 Pi
0161   //       if |theta- theta0| < alpha
0162   //       (theta within cone of opening angle alpha;
0163   //        actually theta0 < alpha in forward/backward direction,
0164   //        so that integral includes all phi values, hence yielding 2 Pi)
0165   if (s <= 0) {
0166     if (TMath::Abs(x - theta0) > alpha_)
0167       return 0;
0168     if (TMath::Abs(x - theta0) <= alpha_)
0169       return 2 * TMath::Pi();
0170     std::cerr << "Error in <IntegralOverPhiFunction::operator()>: failed to compute return value !" << std::endl;
0171   }
0172 
0173   double r = (1. / TMath::Sqrt(2.)) * (cscTheta * cscTheta * cscTheta0 * cscTheta0 * TMath::Sqrt(s) * tanPhi0);
0174   double t = cosPhi0 * (-cotTheta * cotTheta0 + cosAlpha * cscTheta * cscTheta0);
0175 
0176   double phi[4];
0177   phi[0] = -TMath::ACos(t - r);
0178   phi[1] = TMath::ACos(t - r);
0179   phi[2] = -TMath::ACos(t + r);
0180   phi[3] = TMath::ACos(t + r);
0181 
0182   if (debugLevel_ > 0) {
0183     edm::LogVerbatim("") << "phi0 = " << phi0_ << std::endl
0184                          << "phi[0] = " << phi[0] << " (phi[0] - phi0 = " << (phi[0] - phi0_) * 180 / TMath::Pi() << ")"
0185                          << std::endl
0186                          << "phi[1] = " << phi[1] << " (phi[1] - phi0 = " << (phi[1] - phi0_) * 180 / TMath::Pi() << ")"
0187                          << std::endl
0188                          << "phi[2] = " << phi[2] << " (phi[2] - phi0 = " << (phi[2] - phi0_) * 180 / TMath::Pi() << ")"
0189                          << std::endl
0190                          << "phi[3] = " << phi[3] << " (phi[3] - phi0 = " << (phi[3] - phi0_) * 180 / TMath::Pi() << ")"
0191                          << std::endl;
0192   }
0193 
0194   double phiMin = 0.;
0195   double phiMax = 0.;
0196   for (unsigned int i = 0; i < 4; ++i) {
0197     for (unsigned int j = (i + 1); j < 4; ++j) {
0198       //--- search for the two solutions for phi
0199       //    that have an equal distance to phi0
0200       //    and are on either side
0201       double dPhi_i = phi[i] - phi0_;
0202       double dPhi_j = phi0_ - phi[j];
0203       if (TMath::Abs(normalizedPhi(dPhi_i - dPhi_j)) <
0204           epsilon) {  // map difference in azimuth angle into interval [-pi,+pi] and require it to be negligible
0205         //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
0206         phiMin = TMath::Min(phi[i], phi[j]);
0207         phiMax = TMath::Max(phi[i], phi[j]);
0208 
0209         // !!! ONLY FOR TESTING
0210         if (phi[i] == phiMin)
0211           checkSolutions(i, numSolutionMin1_, numSolutionMin2_, numSolutionMin3_, numSolutionMin4_);
0212         if (phi[i] == phiMax)
0213           checkSolutions(i, numSolutionMax1_, numSolutionMax2_, numSolutionMax3_, numSolutionMax4_);
0214         if (phi[j] == phiMin)
0215           checkSolutions(j, numSolutionMin1_, numSolutionMin2_, numSolutionMin3_, numSolutionMin4_);
0216         if (phi[j] == phiMax)
0217           checkSolutions(j, numSolutionMax1_, numSolutionMax2_, numSolutionMax3_, numSolutionMax4_);
0218         //     FOR TESTING ONLY !!!
0219       }
0220     }
0221   }
0222 
0223   if (phiMin == 0 && phiMax == 0) {
0224     edm::LogError("") << "failed to compute Return Value !" << std::endl;
0225   }
0226 
0227   return TMath::Abs(normalizedPhi(phi0_ - phiMin)) + TMath::Abs(normalizedPhi(phiMax - phi0_));
0228 }
0229 
0230 double IntegralOverPhiFunction::DoDerivative(double x) const {
0231   //--- virtual function inherited from ROOT::Math::ParamFunction base class;
0232   //    not implemented, because not neccessary, but needs to be defined to make code compile...
0233   edm::LogWarning("") << "Function not implemented yet !" << std::endl;
0234 
0235   return 0.;
0236 }
0237 
0238 double IntegralOverPhiFunction::DoParameterDerivative(double, const double*, unsigned int) const {
0239   //--- virtual function inherited from ROOT::Math::ParamFunction base class;
0240   //    not implemented, because not neccessary, but needs to be defined to make code compile...
0241   edm::LogWarning("") << "Function not implemented yet !" << std::endl;
0242 
0243   return 0.;
0244 }
0245 
0246 void IntegralOverPhiFunction::DoParameterGradient(double x, double* paramGradient) const {
0247   //--- virtual function inherited from ROOT::Math::ParamFunction base class;
0248   //    not implemented, because not neccessary, but needs to be defined to make code compile...
0249   edm::LogWarning("") << "Function not implemented yet !" << std::endl;
0250 }
0251 
0252 //
0253 // definition of auxiliary functions
0254 //
0255 
0256 void checkSolutions(unsigned int i,
0257                     unsigned int& numSolution1,
0258                     unsigned int& numSolution2,
0259                     unsigned int& numSolution3,
0260                     unsigned int& numSolution4) {
0261   switch (i) {
0262     case 0:
0263       ++numSolution1;
0264       break;
0265     case 1:
0266       ++numSolution2;
0267       break;
0268     case 2:
0269       ++numSolution3;
0270       break;
0271     case 3:
0272       ++numSolution4;
0273       break;
0274   }
0275 }
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.2.1 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!