BeamSpot/src/BeamSpot.cc

0001 /**_________________________________________________________________
0002    class:   BeamSpot.cc
0003    package: DataFormats/BeamSpot
0004
0005  A reconstructed beam spot providing position, width, slopes,
0006  and errors.
0007
0008  author: Francisco Yumiceva, Fermilab (yumiceva@fnal.gov)
0009
0010
0011  ________________________________________________________________**/
0012
0013 #include "DataFormats/BeamSpot/interface/BeamSpot.h"
0014 #include "DataFormats/CLHEP/interface/AlgebraicObjects.h"
0015 #include "DataFormats/GeometrySurface/interface/TkRotation.h"
0016
0017 #include <iostream>
0018
0019 namespace reco {
0020
0021   using namespace math;
0022
0023   BeamSpot::BeamSpot() {
0024     // initialize
0025     position_ = Point(0., 0., 0.);
0026     sigmaZ_ = 0.;
0027     dxdz_ = 0.;
0028     dydz_ = 0.;
0029     BeamWidthX_ = 0.;
0030     BeamWidthY_ = 0;
0031     for (int j = 0; j < 7; j++) {
0032       for (int k = j; k < 7; k++) {
0033         error_(j, k) = 0.;
0034       }
0035     }
0036     type_ = Unknown;
0037     emittanceX_ = 0;
0038     emittanceY_ = 0;
0039     betaStar_ = 0;
0040   }
0041
0042   const BeamSpot::Point BeamSpot::position(const double z) const {
0043     Point pos(x(z), y(z), z);
0044     return pos;
0045   }
0046
0047   void BeamSpot::print(std::stringstream& ss) const {
0048     ss << "-----------------------------------------------------\n"
0049        << "              Beam Spot Data\n\n"
0050        << " Beam type    = " << type() << "\n"
0051        << "       X0     = " << x0() << " +/- " << x0Error() << " [cm]\n"
0052        << "       Y0     = " << y0() << " +/- " << y0Error() << " [cm]\n"
0053        << "       Z0     = " << z0() << " +/- " << z0Error() << " [cm]\n"
0054        << " Sigma Z0     = " << sigmaZ() << " +/- " << sigmaZ0Error() << " [cm]\n"
0055        << " dxdz         = " << dxdz() << " +/- " << dxdzError() << " [radians]\n"
0056        << " dydz         = " << dydz() << " +/- " << dydzError() << " [radians]\n"
0057        << " Beam width X = " << BeamWidthX() << " +/- " << BeamWidthXError() << " [cm]\n"
0058        << " Beam width Y = " << BeamWidthY() << " +/- " << BeamWidthYError() << " [cm]\n"
0059        << " EmittanceX   = " << emittanceX() << " [cm]\n"
0060        << " EmittanceY   = " << emittanceY() << " [cm]\n"
0061        << " beta-star    = " << betaStar() << " [cm]\n"
0062        << "-----------------------------------------------------\n\n";
0063   }
0064
0065   //
0066   std::ostream& operator<<(std::ostream& os, BeamSpot beam) {
0067     std::stringstream ss;
0068     beam.print(ss);
0069     os << ss.str();
0070     return os;
0071   }
0072
0073   BeamSpot::Covariance3DMatrix BeamSpot::rotatedCovariance3D() const {
0074     AlgebraicVector3 newZ(dxdz(), dydz(), 1.);
0075     AlgebraicVector3 globalZ(0., 0., 1.);
0076     AlgebraicVector3 rotationAxis = ROOT::Math::Cross(globalZ.Unit(), newZ.Unit());
0077     float rotationAngle = -acos(ROOT::Math::Dot(globalZ.Unit(), newZ.Unit()));
0078     Basic3DVector<float> aa(rotationAxis[0], rotationAxis[1], rotationAxis[2]);
0079     TkRotation<float> rotation(aa, rotationAngle);
0080     AlgebraicMatrix33 rotationMatrix;
0081     rotationMatrix(0, 0) = rotation.xx();
0082     rotationMatrix(0, 1) = rotation.xy();
0083     rotationMatrix(0, 2) = rotation.xz();
0084     rotationMatrix(1, 0) = rotation.yx();
0085     rotationMatrix(1, 1) = rotation.yy();
0086     rotationMatrix(1, 2) = rotation.yz();
0087     rotationMatrix(2, 0) = rotation.zx();
0088     rotationMatrix(2, 1) = rotation.zy();
0089     rotationMatrix(2, 2) = rotation.zz();
0090
0091     AlgebraicSymMatrix33 diagError;
0092     diagError(0, 0) = pow(BeamWidthX(), 2);
0093     diagError(1, 1) = pow(BeamWidthY(), 2);
0094     diagError(2, 2) = pow(sigmaZ(), 2);
0095
0096     Covariance3DMatrix matrix;
0097     matrix = ROOT::Math::Similarity(rotationMatrix, diagError) + covariance3D();
0098     return matrix;
0099   }
0100
0101 }  // namespace reco