CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_10_patch1/src/RecoVertex/BeamSpotProducer/src/BSpdfsFcn.cc

Go to the documentation of this file.
00001 
00014 #include "RecoVertex/BeamSpotProducer/interface/BSpdfsFcn.h"
00015 #include "TMath.h"
00016 
00017 #include <cmath>
00018 #include <vector>
00019 
00020 //______________________________________________________________________
00021 double BSpdfsFcn::PDFGauss_d(double z, double d, double sigmad,
00022                                                          double phi, const std::vector<double>& parms) const {
00023 
00024   //---------------------------------------------------------------------------
00025   //  PDF for d0 distribution. This PDF is a simple gaussian in the
00026   //  beam reference frame.
00027   //---------------------------------------------------------------------------
00028         double fsqrt2pi = sqrt(2.* TMath::Pi());
00029         
00030         double sig = sqrt(parms[fPar_SigmaBeam]*parms[fPar_SigmaBeam] +
00031                                           sigmad*sigmad);
00032 
00033         double dprime = d - ( ( parms[fPar_X0] + z*parms[fPar_dxdz] )*sin(phi)
00034                 - ( parms[fPar_Y0] + z*parms[fPar_dydz] )*cos(phi) );
00035         
00036         double result = (exp(-(dprime*dprime)/(2.0*sig*sig)))/(sig*fsqrt2pi);
00037                 
00038         return result;
00039 
00040 }
00041 
00042 //______________________________________________________________________
00043 double BSpdfsFcn::PDFGauss_d_resolution(double z, double d, double phi, double pt, const std::vector<double>& parms) const {
00044 
00045   //---------------------------------------------------------------------------
00046   //  PDF for d0 distribution. This PDF is a simple gaussian in the
00047   //  beam reference frame. The IP resolution is parametrize by a linear
00048   //  function as a function of 1/pt.   
00049   //---------------------------------------------------------------------------
00050         double fsqrt2pi = sqrt(2.* TMath::Pi());
00051 
00052         double sigmad = parms[fPar_c0] + parms[fPar_c1]/pt;
00053         
00054         double sig = sqrt(parms[fPar_SigmaBeam]*parms[fPar_SigmaBeam] +
00055                                           sigmad*sigmad);
00056 
00057         double dprime = d - ( ( parms[fPar_X0] + z*parms[fPar_dxdz] )*sin(phi)
00058                 - ( parms[fPar_Y0] + z*parms[fPar_dydz] )*cos(phi) );
00059         
00060         double result = (exp(-(dprime*dprime)/(2.0*sig*sig)))/(sig*fsqrt2pi);
00061                 
00062         return result;
00063 
00064 }
00065 
00066 //______________________________________________________________________
00067 double BSpdfsFcn::PDFGauss_z(double z, double sigmaz, const std::vector<double>& parms) const {
00068 
00069   //---------------------------------------------------------------------------
00070   //  PDF for z-vertex distribution. This distribution
00071   // is parametrized by a simple normalized gaussian distribution.
00072   //---------------------------------------------------------------------------
00073         double fsqrt2pi = sqrt(2.* TMath::Pi());
00074 
00075         double sig = sqrt(sigmaz*sigmaz+parms[fPar_SigmaZ]*parms[fPar_SigmaZ]);
00076         //double sig = sqrt(sigmaz*sigmaz+parms[1]*parms[1]);
00077         double result = (exp(-((z-parms[fPar_Z0])*(z-parms[fPar_Z0]))/(2.0*sig*sig)))/(sig*fsqrt2pi);
00078         //double result = (exp(-((z-parms[0])*(z-parms[0]))/(2.0*sig*sig)))/(sig*fsqrt2pi);
00079         
00080         return result;
00081 
00082 }
00083 
00084 
00085 
00086 
00087 //______________________________________________________________________
00088 double BSpdfsFcn::operator() (const std::vector<double>& params) const {
00089                 
00090         double f = 0.0;
00091 
00092         //std::cout << "fusepdfs=" << fusepdfs << " params.size="<<params.size() << std::endl;
00093         
00094         std::vector<BSTrkParameters>::const_iterator iparam = fBSvector.begin();
00095 
00096         double pdf = 0;
00097         
00098         for( iparam = fBSvector.begin(); iparam != fBSvector.end(); ++iparam) {
00099 
00100                 
00101                 if (fusepdfs == "PDFGauss_z") {
00102                         pdf = PDFGauss_z( iparam->z0(), iparam->sigz0(),params);
00103                 }
00104                 else if (fusepdfs == "PDFGauss_d") {
00105                         pdf = PDFGauss_d( iparam->z0(), iparam->d0(),
00106                                                           iparam->sigd0(), iparam->phi0(),params);
00107                 }
00108                 else if (fusepdfs == "PDFGauss_d_resolution") {
00109                         pdf = PDFGauss_d_resolution( iparam->z0(), iparam->d0(),
00110                                                           iparam->phi0(), iparam->pt(),params);
00111                 }
00112                 else if (fusepdfs == "PDFGauss_d*PDFGauss_z") {
00113                         //std::cout << "pdf= " << pdf << std::endl;
00114                         pdf = PDFGauss_d( iparam->z0(), iparam->d0(),
00115                                                           iparam->sigd0(), iparam->phi0(),params)*
00116                                 PDFGauss_z( iparam->z0(), iparam->sigz0(),params);
00117                 }
00118                 else if (fusepdfs == "PDFGauss_d_resolution*PDFGauss_z") {
00119                         pdf = PDFGauss_d_resolution( iparam->z0(), iparam->d0(),
00120                                                                                  iparam->phi0(), iparam->pt(),params)*
00121                                 PDFGauss_z( iparam->z0(), iparam->sigz0(),params);
00122                 }
00123                 
00124                 f = log(pdf) + f;
00125     }
00126 
00127         f= -2.0*f;
00128         return f;
00129 
00130 }