9 #include "GaudiKernel/PhysicalConstants.h"
12 #include "CLHEP/Vector/LorentzVector.h"
13 #include "CLHEP/Geometry/Point3D.h"
14 #include "CLHEP/Geometry/Transform3D.h"
18 #include "CLHEP/Random/RandGaussZiggurat.h"
40 if (!m_randomEngine) {
41 ATH_MSG_ERROR(
"Could not get random number engine from RandomNumberService. Abort.");
42 return StatusCode::FAILURE;
46 return StatusCode::SUCCESS;
53 return StatusCode::SUCCESS;
60 m_randomEngine->setSeed(
name(), ctx );
61 CLHEP::HepRandomEngine* randomEngine(m_randomEngine->getEngine(ctx));
66 float vertexX = CLHEP::RandGaussZiggurat::shoot(randomEngine)*beamSpotHandle->beamSigma(0);
67 float vertexY = CLHEP::RandGaussZiggurat::shoot(randomEngine)*beamSpotHandle->beamSigma(1);
68 float vertexZ = CLHEP::RandGaussZiggurat::shoot(randomEngine)*beamSpotHandle->beamSigma(2);
70 CLHEP::HepLorentzVector *vertexSmearing =
71 new CLHEP::HepLorentzVector( vertexX, vertexY, vertexZ, 0. );
74 const double tx =
tan( beamSpotHandle->beamTilt(1) );
75 const double ty =
tan( beamSpotHandle->beamTilt(0) );
77 const double sqrt_abc = sqrt(1. +
tx*
tx + ty*ty);
78 const double sqrt_fgh = sqrt(1. + ty*ty);
80 const double a = ty/sqrt_abc;
81 const double b =
tx/sqrt_abc;
82 const double c = 1./sqrt_abc;
84 HepGeom::Point3D<double> from1(0,0,1);
85 HepGeom::Point3D<double> to1(
a,
b,
c);
87 const double f = 1./sqrt_fgh;
89 const double h = -(ty)/sqrt_fgh;
91 HepGeom::Point3D<double> from2(1,0,0);
92 HepGeom::Point3D<double> to2(
f,
g,
h);
96 HepGeom::Rotate3D(from1, from2, to1, to2).getRotation(),
97 CLHEP::Hep3Vector( beamSpotHandle->beamPos().x(),
98 beamSpotHandle->beamPos().y(),
99 beamSpotHandle->beamPos().z() )
103 ATH_MSG_VERBOSE(
"BeamSpotSvc reported beam position as " << beamSpotHandle->beamPos() << std::endl
104 <<
"\tWidth is (" << beamSpotHandle->beamSigma(0)
105 <<
", " << beamSpotHandle->beamSigma(1) <<
", "
106 << beamSpotHandle->beamSigma(2) <<
")" << std::endl
108 <<
"\tTilts are " << beamSpotHandle->beamTilt(0) <<
" and " << beamSpotHandle->beamTilt(1) << std::endl
109 <<
"\tVertex Position before transform: " << *vertexSmearing);
112 *vertexSmearing =
transform * HepGeom::Point3D<double>(*vertexSmearing);
115 vertexSmearing->setT(vertexT);
118 return vertexSmearing;