ATLAS Offline Software
LongBeamspotVertexPositioner.cxx
Go to the documentation of this file.
1 /*
2  Copyright (C) 2002-2022 CERN for the benefit of the ATLAS collaboration
3 */
4 
5 // class header include
7 
8 // For the speed of light
9 #include "GaudiKernel/PhysicalConstants.h"
10 
11 // CLHEP includes
12 #include "CLHEP/Vector/LorentzVector.h"
13 #include "CLHEP/Geometry/Point3D.h"
14 #include "CLHEP/Geometry/Transform3D.h"
15 #include <math.h> /* erf */
16 
17 // RandomNumber generator
19 #include "CLHEP/Random/RandGaussZiggurat.h"
20 #include "CLHEP/Random/RandFlat.h"
21 
22 namespace Simulation
23 {
24 
27  const std::string& n,
28  const IInterface* p )
29  : base_class(t,n,p)
30  {
31  }
32 
35  {
36  ATH_MSG_VERBOSE("Initializing ...");
37 
39  // prepare the RandonNumber generation
40  ATH_CHECK(m_rndGenSvc.retrieve());
41  m_randomEngine = m_rndGenSvc->getEngine(this, m_randomEngineName);
42  if (!m_randomEngine) {
43  ATH_MSG_ERROR("Could not get random number engine from RandomNumberService. Abort.");
44  return StatusCode::FAILURE;
45  }
46 
47  // everything set up properly
48  return StatusCode::SUCCESS;
49  }
50 
51 
54  {
55  ATH_MSG_VERBOSE("Finalizing ...");
56  return StatusCode::SUCCESS;
57  }
58 
60  {
61  //double norm(1.0/232.06);
62  double temp(1.0-std::abs(z)/m_L);
63  return erf(2.5*temp)*heaviside(temp); // zero for |z|>150.0
64  }
65 
66  double LongBeamspotVertexPositioner::getZpos(CLHEP::HepRandomEngine* randomEngine) const
67  {
68  size_t ntries(0);
69  double yval(CLHEP::RandFlat::shoot(randomEngine, 0.0, 1.0));
70  double zval(CLHEP::RandFlat::shoot(randomEngine, -150.0, 150.0));
71  while (this->beamspotFunction(zval)<yval) {
72  if(ntries>1000000) return 0.0; //just so we don't sit in this loop forever
73  yval = CLHEP::RandFlat::shoot(randomEngine, 0.0, 1.0);
74  zval = CLHEP::RandFlat::shoot(randomEngine, -150.0, 150.0);
75  ++ntries;
76  }
77  return zval;
78  }
79 
81  CLHEP::HepLorentzVector *LongBeamspotVertexPositioner::generate(const EventContext& ctx) const
82  {
83  // Prepare the random engine
84  m_randomEngine->setSeed( name(), ctx );
85  CLHEP::HepRandomEngine* randomEngine(m_randomEngine->getEngine(ctx));
87  // See jira issue ATLASSIM-497 for an explanation of why calling
88  // shoot outside the CLHEP::HepLorentzVector constructor is
89  // necessary/preferable.
90  float vertexX = CLHEP::RandGaussZiggurat::shoot(randomEngine)*beamSpotHandle->beamSigma(0);
91  float vertexY = CLHEP::RandGaussZiggurat::shoot(randomEngine)*beamSpotHandle->beamSigma(1);
92  float vertexZ = getZpos(randomEngine);
93  // calculate the vertexSmearing
94  CLHEP::HepLorentzVector *vertexSmearing =
95  new CLHEP::HepLorentzVector( vertexX, vertexY, vertexZ, 0. );
96 
97  // (1) code from: Simulation/G4Atlas/G4AtlasUtilities/VertexPositioner.cxx
98  const double tx = tan( beamSpotHandle->beamTilt(1) );
99  const double ty = tan( beamSpotHandle->beamTilt(0) );
100 
101  const double sqrt_abc = sqrt(1. + tx*tx + ty*ty);
102  const double sqrt_fgh = sqrt(1. + ty*ty);
103 
104  const double a = ty/sqrt_abc;
105  const double b = tx/sqrt_abc;
106  const double c = 1./sqrt_abc;
107 
108  HepGeom::Point3D<double> from1(0,0,1);
109  HepGeom::Point3D<double> to1(a,b,c);
110 
111  const double f = 1./sqrt_fgh;
112  const double g = 0.;
113  const double h = -(ty)/sqrt_fgh;
114 
115  HepGeom::Point3D<double> from2(1,0,0);
116  HepGeom::Point3D<double> to2(f,g,h);
117 
118  // first rotation, then translation
120  HepGeom::Rotate3D(from1, from2, to1, to2).getRotation(),
121  CLHEP::Hep3Vector( beamSpotHandle->beamPos().x(),
122  beamSpotHandle->beamPos().y(),
123  beamSpotHandle->beamPos().z() )
124  );
125 
126  // FIXME: don't use endl in MsgStream printouts
127  ATH_MSG_VERBOSE("BeamSpotSvc reported beam position as " << beamSpotHandle->beamPos() << std::endl
128  << "\tWidth is (" << beamSpotHandle->beamSigma(0)
129  << ", " << beamSpotHandle->beamSigma(1) << ", "
130  << m_L << ")" << std::endl
131  << "\tTilts are " << beamSpotHandle->beamTilt(0) << " and " << beamSpotHandle->beamTilt(1) << std::endl
132  << "\tVertex Position before transform: " << *vertexSmearing);
133 
134  // update with the tilt
135  *vertexSmearing = transform * HepGeom::Point3D<double>(*vertexSmearing);
136 
137  // See if we were asked to do time smearing as well
138  if (m_timeSmearing) {
139  /* This is ballpark code courtesy of Brian Amadio. He provided some functions based on beam parameters.
140  He provided a little trick for pulling out the beam bunch width as well. Hard coding the crossing angle
141  parameter for the time being, as the beam spot service doesn't really provide that yet. */
142  // Assumes that to make these funny beam spots we have a normal-looking bunch
143  double bunch_length_z = (std::sqrt(2)*vertexZ)/0.9; // 0.9 is the crossing angle reduction factor
144  // double tLimit = 2.*(bunch_length_z+bunch_length_z)/Gaudi::Units::c_light;
145  // TF1 func = TF1("func","[0]*exp((-([3]-299792458*x)^2*[2]^2-([3]+299792458*x)^2*[1]^2)/(2*[1]^2*[2]^2))",-1*tLimit,tLimit);
146  // func.SetParameter(0,Gaudi::Units::c_light/(M_PI*bunch_length_z*bunch_length_z));
147  // func.SetParameter(1,bunch_length_z);
148  // func.SetParameter(2,bunch_length_z);
149  // func.SetParameter(3,vertexSmearing->z());
150  // double time_offset = func.GetRandom();
151 
152  // Time should be set in units of distance, which is a little funny
153  double time_offset = CLHEP::RandGaussZiggurat::shoot(
154  randomEngine, vertexSmearing->z()/Gaudi::Units::c_light,
155  bunch_length_z/Gaudi::Units::c_light );
156 
157  vertexSmearing->setT( vertexSmearing->t() + time_offset*Gaudi::Units::c_light );
158  }
159 
160  // and return it
161  return vertexSmearing;
162  }
163 
164 } // namespace Simulation
Simulation::LongBeamspotVertexPositioner::m_beamSpotKey
SG::ReadCondHandleKey< InDet::BeamSpotData > m_beamSpotKey
Definition: LongBeamspotVertexPositioner.h:61
Simulation::LongBeamspotVertexPositioner::initialize
StatusCode initialize() override final
Athena algtool's Hooks.
Definition: LongBeamspotVertexPositioner.cxx:34
LongBeamspotVertexPositioner.h
SG::ReadCondHandle
Definition: ReadCondHandle.h:44
Simulation::LongBeamspotVertexPositioner::m_timeSmearing
Gaudi::Property< bool > m_timeSmearing
Do time smearing.
Definition: LongBeamspotVertexPositioner.h:66
read_hist_ntuple.t
t
Definition: read_hist_ntuple.py:5
ATH_MSG_VERBOSE
#define ATH_MSG_VERBOSE(x)
Definition: AthMsgStreamMacros.h:28
python.utils.AtlRunQueryDQUtils.p
p
Definition: AtlRunQueryDQUtils.py:210
ATH_MSG_ERROR
#define ATH_MSG_ERROR(x)
Definition: AthMsgStreamMacros.h:33
z
#define z
python.CaloCondTools.g
g
Definition: CaloCondTools.py:15
beamspotman.n
n
Definition: beamspotman.py:731
EL::StatusCode
::StatusCode StatusCode
StatusCode definition for legacy code.
Definition: PhysicsAnalysis/D3PDTools/EventLoop/EventLoop/StatusCode.h:22
Amg::Transform3D
Eigen::Affine3d Transform3D
Definition: GeoPrimitives.h:46
Amg::transform
Amg::Vector3D transform(Amg::Vector3D &v, Amg::Transform3D &tr)
Transform a point from a Trasformation3D.
Definition: GeoPrimitivesHelpers.h:156
ATH_CHECK
#define ATH_CHECK
Definition: AthCheckMacros.h:40
Simulation::LongBeamspotVertexPositioner::m_L
Gaudi::Property< double > m_L
Parameter in the Z distribution of the beamspot - default 150.0 mm.
Definition: LongBeamspotVertexPositioner.h:60
hist_file_dump.f
f
Definition: hist_file_dump.py:135
drawFromPickle.tan
tan
Definition: drawFromPickle.py:36
name
std::string name
Definition: Control/AthContainers/Root/debug.cxx:228
plotBeamSpotMon.b
b
Definition: plotBeamSpotMon.py:77
SG::CondHandleKey::initialize
StatusCode initialize(bool used=true)
python.PhysicalConstants.c_light
float c_light
Definition: PhysicalConstants.py:63
plotting.yearwise_efficiency_vs_mu.yval
float yval
Definition: yearwise_efficiency_vs_mu.py:36
RNGWrapper.h
a
TList * a
Definition: liststreamerinfos.cxx:10
h
Simulation::LongBeamspotVertexPositioner::finalize
StatusCode finalize() override final
Athena algtool's Hooks.
Definition: LongBeamspotVertexPositioner.cxx:53
Simulation::LongBeamspotVertexPositioner::m_rndGenSvc
ServiceHandle< IAthRNGSvc > m_rndGenSvc
Definition: LongBeamspotVertexPositioner.h:62
Simulation
Definition: BeamEffectsAlg.cxx:21
Simulation::LongBeamspotVertexPositioner::LongBeamspotVertexPositioner
LongBeamspotVertexPositioner(const std::string &t, const std::string &n, const IInterface *p)
Constructor with parameters.
Definition: LongBeamspotVertexPositioner.cxx:26
Simulation::LongBeamspotVertexPositioner::heaviside
double heaviside(double val) const
Definition: LongBeamspotVertexPositioner.h:57
Simulation::LongBeamspotVertexPositioner::m_randomEngineName
Gaudi::Property< std::string > m_randomEngineName
Name of the random number stream.
Definition: LongBeamspotVertexPositioner.h:65
Simulation::LongBeamspotVertexPositioner::beamspotFunction
double beamspotFunction(double z) const
Definition: LongBeamspotVertexPositioner.cxx:59
python.compressB64.c
def c
Definition: compressB64.py:93
TileDCSDataPlotter.tx
tx
Definition: TileDCSDataPlotter.py:878
Simulation::LongBeamspotVertexPositioner::getZpos
double getZpos(CLHEP::HepRandomEngine *) const
Definition: LongBeamspotVertexPositioner.cxx:66
Simulation::LongBeamspotVertexPositioner::generate
CLHEP::HepLorentzVector * generate(const EventContext &ctx) const override final
computes the vertex displacement
Definition: LongBeamspotVertexPositioner.cxx:81