19#include "BDSBunchGaussian.hh"
21#include "BDSException.hh"
22#include "BDSGlobalConstants.hh"
24#include "parser/beam.h"
26#include "Randomize.hh"
27#include "CLHEP/Matrix/SymMatrix.h"
28#include "CLHEP/Matrix/Vector.h"
29#include "CLHEP/RandomObjects/RandMultiGauss.h"
30#include "CLHEP/Units/PhysicalConstants.h"
38 bool isPositiveDefinite(CLHEP::HepSymMatrix& S)
40 CLHEP::HepSymMatrix temp (S);
42 CLHEP::HepMatrix U = diagonalize ( &temp );
43 CLHEP::HepSymMatrix D = S.similarityT(U);
44 for (G4int i = 1; i <= S.num_row(); i++)
54BDSBunchGaussian::BDSBunchGaussian(
const G4String& nameIn):
56 meansGM(CLHEP::HepVector(6)),
57 sigmaGM(CLHEP::HepSymMatrix(6)),
58 gaussMultiGen(nullptr),
59 offsetSampleMean(false),
62 coordinates = {&x0_v, &xp_v, &y0_v, &yp_v, &z0_v, &zp_v, &E_v, &t_v, &weight_v};
65BDSBunchGaussian::~BDSBunchGaussian()
73 G4Transform3D beamlineTransformIn,
74 const G4double beamlineSIn)
90 meansGM[0] =
X0 / CLHEP::m;
92 meansGM[2] =
Y0 / CLHEP::m;
94 meansGM[4] =
T0 / CLHEP::s;
111 const CLHEP::HepVector& mu,
112 CLHEP::HepSymMatrix& sigma)
116 G4cout << __METHOD_NAME__ <<
"checking 6x6 sigma matrix is positive definite" << G4endl;
117 if (!isPositiveDefinite(sigma))
119 G4cout << __METHOD_NAME__ <<
"WARNING bunch generator sigma matrix is not positive definite" << G4endl;
120 G4cout << sigma << G4endl;
121 G4cout << __METHOD_NAME__ <<
"adding a small error to zero diagonal elements" << G4endl;
124 G4double small_error = 1e-50;
126 for (G4int i=0; i<6; i++)
129 {sigma[i][i] += small_error;}
132 if (!isPositiveDefinite(sigma))
134 G4cout << __METHOD_NAME__ <<
"WARNING bunch generator sigma matrix is still not positive definite" << G4endl;
135 G4cout << sigma << G4endl;
136 G4cout << __METHOD_NAME__ <<
"adding a small error to all elements" << G4endl;
137 for (G4int i=0; i<6; i++)
139 for (G4int j=0; j<6; j++)
142 {sigma[i][j] += small_error;}
145 if (!isPositiveDefinite(sigma))
147 G4cout << sigma << G4endl;
148 throw BDSException(__METHOD_NAME__,
"bunch generator sigma matrix is still not positive definite, giving up");
154 G4cout << __METHOD_NAME__ <<
"mean " << G4endl
156 << __METHOD_NAME__ <<
"sigma " << G4endl
159 G4cout << __METHOD_NAME__ <<
"confirmed: positive definite matrix" << G4endl;
160 return new CLHEP::RandMultiGauss(anEngine,mu,sigma);
165 G4cout << __METHOD_NAME__ <<
"Pregenerating " << nGenerate <<
" events." << G4endl;
167 G4double x_a = 0.0, xp_a = 0.0, y_a = 0.0, yp_a = 0.0;
168 G4double z_a = 0.0, zp_a = 0.0, E_a = 0.0, t_a = 0.0;
170 for (G4int iParticle = 0; iParticle < nGenerate; ++iParticle)
174 G4double nT = (G4double)iParticle + 1;
179 xp_a = xp_a + (d/nT);
183 yp_a = yp_a + (d/nT);
187 zp_a = zp_a + (d/nT);
188 d = fire.totalEnergy - E_a;
193 x0_v.push_back(fire.x);
194 xp_v.push_back(fire.xp);
195 y0_v.push_back(fire.y);
196 yp_v.push_back(fire.yp);
197 z0_v.push_back(fire.z);
198 zp_v.push_back(fire.zp);
199 E_v.push_back(fire.totalEnergy);
200 t_v.push_back(fire.T);
215 for (G4int iParticle = 0; iParticle < nGenerate; ++iParticle)
217 x0_v[iParticle] -= x_a;
218 xp_v[iParticle] -= xp_a;
219 y0_v[iParticle] -= y_a;
220 yp_v[iParticle] -= yp_a;
221 z0_v[iParticle] -= z_a;
222 zp_v[iParticle] -= zp_a;
223 E_v[iParticle] -= E_a;
224 t_v[iParticle] -= t_a;
256 G4double x = v[0] * CLHEP::m;
258 G4double y = v[2] * CLHEP::m;
269 dz = t * CLHEP::c_light;
CLHEP::RandMultiGauss * CreateMultiGauss(CLHEP::HepRandomEngine &anEngine, const CLHEP::HepVector &mu, CLHEP::HepSymMatrix &sigma)
std::vector< G4double > z0_v
Holder for pre-calculated coordinates.
std::vector< G4double > y0_v
Holder for pre-calculated coordinates.
std::vector< G4double > xp_v
Holder for pre-calculated coordinates.
virtual BDSParticleCoordsFull GetNextParticleLocal()
std::vector< G4double > t_v
Holder for pre-calculated coordinates.
CLHEP::RandMultiGauss * gaussMultiGen
Randon number generator with sigma matrix and mean.
virtual void SetOptions(const BDSParticleDefinition *beamParticle, const GMAD::Beam &beam, const BDSBunchType &distrType, G4Transform3D beamlineTransformIn=G4Transform3D::Identity, const G4double beamlineS=0)
std::vector< std::vector< G4double > * > coordinates
Convenience vector of vectors for clearing up.
std::vector< G4double > weight_v
Holder for pre-calculated coordinates.
virtual void BeginOfRunAction(G4int numberOfEvents)
std::vector< G4double > zp_v
Holder for pre-calculated coordinates.
G4int iPartIteration
Iterator for reading out pre-calculate coordinates.
std::vector< G4double > yp_v
Holder for pre-calculated coordinates.
virtual BDSParticleCoordsFull GetNextParticleLocalCoords()
Fire random number generator and get coordinates. Can be overloaded if required.
std::vector< G4double > E_v
Holder for pre-calculated coordinates.
void PreGenerateEvents(G4int nGenerate)
Pre-generate all the particle coordinates and subtract the sample mean.
G4bool offsetSampleMean
Whether to offset the sample mean.
std::vector< G4double > x0_v
Holder for pre-calculated coordinates.
The base class for bunch distribution generators.
G4double sigmaE
Centre of distributions.
G4double Yp0
Centre of distributions.
G4double T0
Centre of distributions.
G4bool finiteSigmaE
Flags to ignore random number generator in case of no finite E or T.
G4double S0
Centre of distributions.
G4double sigmaT
Centre of distributions.
G4double Z0
Centre of distributions.
G4double X0
Centre of distributions.
G4double Zp0
Centre of distributions.
G4double Xp0
Centre of distributions.
G4double E0
Centre of distributions.
G4double Y0
Centre of distributions.
virtual void SetOptions(const BDSParticleDefinition *beamParticle, const GMAD::Beam &beam, const BDSBunchType &distrType, G4Transform3D beamlineTransformIn=G4Transform3D::Identity, const G4double beamlineS=0)
G4bool finiteSigmaT
Flags to ignore random number generator in case of no finite E or T.
static G4double CalculateZp(G4double xp, G4double yp, G4double Zp0)
Calculate zp safely based on other components.
General exception with possible name of object and message.
A set of particle coordinates including energy and weight.
Wrapper for particle definition.
Improve type-safety of native enum data type in C++.