00001 #include "BDSBunchTwiss.hh"
00002 #include "BDSDebug.hh"
00003
00004 BDSBunchTwiss::BDSBunchTwiss() :
00005 BDSBunchInterface(), betaX(0.0), betaY(0.0), alphaX(0.0), alphaY(0.0), emitX(0.0), emitY(0.0), gammaX(0.0), gammaY(0.0)
00006 {
00007 #ifdef BDSDEBUG
00008 G4cout << __METHOD_NAME__ << G4endl;
00009 #endif
00010
00011 GaussMultiGen = NULL;
00012 }
00013
00014 BDSBunchTwiss::BDSBunchTwiss(G4double betaXIn, G4double betaYIn,
00015 G4double alphaXIn, G4double alphaYIn,
00016 G4double emitXIn, G4double emitYIn,
00017 G4double X0In, G4double Y0In, G4double Z0In, G4double T0In,
00018 G4double Xp0In, G4double Yp0In, G4double Zp0In,
00019 G4double sigmaTIn, G4double sigmaEIn) :
00020 BDSBunchInterface(X0In,Y0In,Z0In,T0In,Xp0In,Yp0In,Zp0In,sigmaTIn,sigmaEIn), betaX(betaXIn), betaY(betaYIn), alphaX(alphaXIn), alphaY(alphaYIn), emitX(emitXIn), emitY(emitYIn)
00021 {
00022 #ifdef BDSDEBUG
00023 G4cout << __METHOD_NAME__ << G4endl;
00024 #endif
00025
00026 GaussMultiGen = NULL;
00027
00028 sigmaT = sigmaTIn;
00029 sigmaE = sigmaEIn;
00030 gammaX = (1.0+alphaX*alphaX)/betaX;
00031 gammaY = (1.0+alphaY*alphaY)/betaY;
00032
00033 CommonConstruction();
00034 }
00035
00036 BDSBunchTwiss::~BDSBunchTwiss() {
00037 #ifdef BDSDEBUG
00038 G4cout << __METHOD_NAME__ << G4endl;
00039 #endif
00040
00041 delete GaussMultiGen;
00042 }
00043
00044 void BDSBunchTwiss::SetOptions(struct Options& opt) {
00045 #ifdef BDSDEBUG
00046 G4cout << __METHOD_NAME__ << G4endl;
00047 #endif
00048
00049 BDSBunchInterface::SetOptions(opt);
00050 SetBetaX(opt.betx);
00051 SetBetaY(opt.bety);
00052 SetAlphaX(opt.alfx);
00053 SetAlphaY(opt.alfy);
00054 SetEmitX(opt.emitx);
00055 SetEmitY(opt.emity);
00056 gammaX = (1.0+alphaX*alphaX)/betaX;
00057 gammaY = (1.0+alphaY*alphaY)/betaY;
00058
00059
00060 CommonConstruction();
00061
00062 return;
00063 }
00064
00065 void BDSBunchTwiss::CommonConstruction() {
00066 #ifdef BDSDEBUG
00067 G4cout << __METHOD_NAME__ << G4endl;
00068 #endif
00069
00070 meansGM = CLHEP::HepVector(6);
00071
00072
00073 meansGM[0] = X0;
00074 meansGM[1] = Xp0;
00075 meansGM[2] = Y0;
00076 meansGM[3] = Yp0;
00077 meansGM[4] = T0;
00078 meansGM[5] = 1;
00079
00080 sigmaGM = CLHEP::HepSymMatrix(6);
00081
00082
00083 sigmaGM[0][0] = emitX*betaX;
00084 sigmaGM[0][1] = -emitX*alphaX;
00085 sigmaGM[1][0] = -emitX*alphaX;
00086 sigmaGM[1][1] = emitX*gammaX;
00087 sigmaGM[2][2] = emitY*betaY;
00088 sigmaGM[2][3] = -emitY*alphaY;
00089 sigmaGM[3][2] = -emitY*alphaY;
00090 sigmaGM[3][3] = emitY*gammaY;
00091 sigmaGM[4][4] = pow(sigmaT,2);
00092 sigmaGM[5][5] = pow(sigmaE,2);
00093
00094 delete GaussMultiGen;
00095 GaussMultiGen = CreateMultiGauss(*CLHEP::HepRandom::getTheEngine(),meansGM,sigmaGM);
00096
00097 return;
00098 }
00099
00100 void BDSBunchTwiss::GetNextParticle(G4double& x0, G4double& y0, G4double& z0,
00101 G4double& xp, G4double& yp, G4double& zp,
00102 G4double& t , G4double& E, G4double& weight) {
00103 #ifdef BDSDEBUG
00104 G4cout << __METHOD_NAME__ << G4endl;
00105 #endif
00106
00107 CLHEP::HepVector v = GaussMultiGen->fire();
00108 x0 = v[0] * CLHEP::m;
00109 xp = v[1] * CLHEP::rad;
00110 y0 = v[2] * CLHEP::m;
00111 yp = v[3] * CLHEP::rad;
00112 t = v[4] * CLHEP::s;
00113 zp = 0.0 * CLHEP::rad;
00114 z0 = Z0*CLHEP::m + t*CLHEP::c_light;
00115 E = BDSGlobalConstants::Instance()->GetParticleKineticEnergy() * v[5];
00116
00117 zp = CalculateZp(xp,yp,Zp0);
00118
00119 weight = 1.0;
00120 return;
00121 }
