Bunch Class Reference

Detailed Description

class Bunch represents a bunch.

It maintains the trajectory for a reference particle and computes the field for all particles in the bunch related to that reference.

A Bunch is created in PreUserTrackingAction() when tracking a reference particle. The tracks in the bunch are added in beginCollectiveTracking(). The fields from the Bunch are computed in computeField(), which is called by collectiveStep(), and they are then used for tracking via addFieldValue().

The grid for PoissonConvolve3D is always centered at the reference particle, and moves with it between time steps. No solution is possible if any charge is at or very close to the boundary, so the specified percentile of the charge distribution is used to determine the size of the grid, keeping the percentile of the distribution between MINGRID and MAXGRID of the grid size. Any charge located >SUPPRESSGRID of the grid size is suppressed (in x, y, or z).

The challenge is to provide suitably accurate boundary conditions without taking too much CPU. For a perfectly Gaussian bunch that can be done by applying an approximation grid 4 or 8 times smaller than the Poisson3D grid. For a cylindrical bunch that does not work. So for now it is left up to the user to specify the approximation grid size. Unlike the Poisson3D grid, the approximation grid can be any size. An attempt is made to merge cells with tiny charge (<1% of the total) with its neighbor, but this does not always reduce the number of approximation charges very much.

This version requires the reference particle's momentum to be parallel to the z axis.

Public Member Functions

Bunch (G4ParticleDefinition *def, int nx, int ny, int nz, double dx, double dy, double dz, double _charge, double _maxBeta, bool _fixedGrid, double _percentile, int nxApprox, int nyApprox, int nzApprox, int _minActive)

void addReferencePoint (const G4Track *track)

add a point to the reference trajectory.

void addReferencePoint (double t, G4ThreeVector pos, G4ThreeVector mom)

add a point to the reference trajectory.

int nTrajectoryPoints ()

nTrajectoryPoints() returns the # points in the reference trajectory

int getNParticles () const

getNParticles() returns # particles in bunch

G4ParticleDefinition * getDefinition ()

getParticleDefinition()

bool addTrack (G4Track *t)

Add a track to the Bunch. returns false if this track is not part of this bunch.

int size ()

int getNapprox ()

int getMinActive ()

G4ThreeVector getBound ()

bool getReferencePos4Mom (double time, G4ThreeVector &position, G4LorentzVector &fourMomentum) const

Interpolate in the point[] vector to get position and 4-momentum of the reference for this bunch. Returns true if OK, false if time is after the valid range for this reference. If time is before the reference was created, its first two trajectory points are linearly extrapolated backwards. Because it must search for the enclosing points of the trajectory, this is most efficient when time does not jump very much (it starts the search at the place in the trajectory of the previous call). (Note: the const is a lie, but is required by BLElementField; only index is modified, which is essentially benign.).

bool init ()

init() initializes the computation, sizing the grid for the initial bunch, and computing the field for the initial time step. Called after all tracks are added to the bunch, before stepping. Returns false if failure.

bool isInBunch (G4Track *track)

isInBunch() returns true if this track is in the bunch. Tracks are selected by these criteria: same particle as reference within bound_x,bound_y,bound_z of the bunch's reference traj. v/c in the bunch's beam frame less than maxBeta

bool computeField ()

All tracks in the bunch are added into the PoissonConvolve3D and the field is computed. returns false if current step time is beyond the range of the bunch's reference.

bool resizeGrid (bool force=false)

resizes the grid, based on current tracks. returns true if grid was resized.

void addFieldValue (const G4double point[4], G4double field[6]) const

addFieldValue adds the field for this bunch to field[].

Static Public Member Functions

static void testEfield ()

testEfield() will test the units and value of the E field in the beam frame for a point charge.

static void testBfield ()

testBfield() will test the B field in the lab frame for a bunch that is a line of moving charges, 1000 Amps.

static void setVerbose (int v)

setVerbose() sets the debug print level

static void setUseApproximationOnly (int v)

setUseApproximationOnly() sets it

Private Attributes

std::vector< Point > point

unsigned index

double prevTime

G4ParticleDefinition * definition

double charge

double maxBeta

double bound_x

double bound_y

double bound_z

int nParticles

double percentile

int minActive

PoissonConvolve3D * poisson

bool goodGuess

bool fixedGrid

Approx approx

std::vector< G4Track * > tracks

BLRunManager * runManager

Static Private Attributes

static const double extrapolateFrac = 0.01

static const double minDeltaT = 1.0e-6*ns

static int verbose = 0

static int useApproximationOnly = 0

Friends

class Approx

Classes

class Approx

class Approx does a 3-d histogram of the charge to generate an approximation for PoissonConvolve3D. Within each bin (voxel) it use the mean position and the total charge. More...

struct Point

Constructor & Destructor Documentation

Bunch::Bunch	(	G4ParticleDefinition *	def,
		int	nx,
		int	ny,
		int	nz,
		double	dx,
		double	dy,
		double	dz,
		double	_charge,
		double	_maxBeta,
		bool	_fixedGrid,
		double	_percentile,
		int	nxApprox,
		int	nyApprox,
		int	nzApprox,
		int	_minActive
	)			`[inline]`

References approx, bound_x, bound_y, bound_z, charge, definition, fixedGrid, BLRunManager::getObject(), goodGuess, index, maxBeta, minActive, nParticles, percentile, poisson, prevTime, Bunch::Approx::reset(), and runManager.

00144                                                                           :
00145                 point(), approx(nxApprox,nyApprox,nzApprox), tracks()
00146                 { index=0; prevTime=-DBL_MAX; definition=def; 
00147                   charge=_charge; maxBeta=_maxBeta; fixedGrid=_fixedGrid;
00148                   bound_x=dx; bound_y=dy; bound_z=dz;
00149                   nParticles=0; percentile=_percentile; minActive=_minActive;
00150                   poisson = new PoissonConvolve3D(nx,ny,nz,-dx,dx,-dy,dy,-dz,dz);
00151                   goodGuess=false; runManager=BLRunManager::getObject();
00152                   approx.reset(this);
00153                 }

Member Function Documentation

void Bunch::addReferencePoint ( const G4Track * track ) [inline]

add a point to the reference trajectory.

References minDeltaT, point, and prevTime.

Referenced by BLCMDspacecharge::UserSteppingAction().

00156                                                      {
00157                 double time=track->GetGlobalTime();
00158                 if(time >= prevTime+minDeltaT) { // valid for prevTime=-DBL_MAX
00159                         point.push_back(Point(time, track->GetPosition(),
00160                                 track->GetMomentum()));
00161                         prevTime = time;
00162                 }
00163         }

void Bunch::addReferencePoint	(	double	t,
		G4ThreeVector	pos,
		G4ThreeVector	mom
	)			`[inline]`

add a point to the reference trajectory.

References point.

00167                 { point.push_back(Point(t,pos,mom)); }

int Bunch::nTrajectoryPoints ( ) [inline]

nTrajectoryPoints() returns the # points in the reference trajectory

References point.

Referenced by BLCMDspacecharge::PostUserTrackingAction().

00170 { return point.size(); }

int Bunch::getNParticles ( ) const [inline]

getNParticles() returns # particles in bunch

References nParticles.

00173 { return nParticles; }

G4ParticleDefinition* Bunch::getDefinition ( ) [inline]

getParticleDefinition()

References definition.

00176 { return definition; }

bool Bunch::addTrack ( G4Track * t ) [inline]

Add a track to the Bunch. returns false if this track is not part of this bunch.

References isInBunch(), and tracks.

00180                                   {
00181                 if(!isInBunch(t)) return false;
00182                 tracks.push_back(t);
00183                 return true;
00184         }

int Bunch::size ( ) [inline]

References tracks.

00186 { return tracks.size(); }

int Bunch::getNapprox ( ) [inline]

References approx, and Bunch::Approx::getVector().

00187 { return approx.getVector().size(); }

int Bunch::getMinActive ( ) [inline]

References minActive.

00188 { return minActive; }

G4ThreeVector Bunch::getBound ( ) [inline]

References bound_x, bound_y, and bound_z.

00190                 { return G4ThreeVector(bound_x,bound_y,bound_z); }

bool Bunch::getReferencePos4Mom	(	double	time,
		G4ThreeVector &	position,
		G4LorentzVector &	fourMomentum
	)			const

Interpolate in the point[] vector to get position and 4-momentum of the reference for this bunch. Returns true if OK, false if time is after the valid range for this reference. If time is before the reference was created, its first two trajectory points are linearly extrapolated backwards. Because it must search for the enclosing points of the trajectory, this is most efficient when time does not jump very much (it starts the search at the place in the trajectory of the previous call). (Note: the const is a lie, but is required by BLElementField; only index is modified, which is essentially benign.).

References BLAssert, definition, index, minDeltaT, point, and verbose.

Referenced by addFieldValue(), computeField(), init(), isInBunch(), and resizeGrid().

00619 {
00620         // find the current entry in point[] -- return false if beyond the end.
00621         while(time < point[index].time && index > 0)
00622                 --((Bunch*)this)->index; // discard const
00623         while(index+1 < point.size()-1 && time > point[index+1].time) {
00624                 ++((Bunch*)this)->index; // discard const
00625         }
00626         if(index+1 >= point.size()) return false;
00627         if(!(index+1 < point.size()))
00628                 printf("index=%d point.size=%ld time=%.3f last_time=%.3f\n",
00629                         index,(long)point.size(),time,point[point.size()-1].time);
00630         BLAssert(index+1 < point.size());
00631 
00632         // Because minDeltaT was imposed while filling point[], we know the
00633         // points are distinct.
00634         double dt = point[index+1].time - point[index].time;
00635         BLAssert(dt >= minDeltaT);
00636 
00637         // interpolate linearly. Works correctly even if time < point[0].time.
00638         double f=(time-point[index].time)/dt;
00639         position = (1.0-f)*point[index].position + f*point[index+1].position;
00640         G4ThreeVector momentum = (1.0-f)*point[index].momentum +
00641                                                 f*point[index+1].momentum;
00642         fourMomentum.setVectM(momentum,definition->GetPDGMass());
00643 
00644         if(verbose >= 3) printf("Bunch::getReferencePos4Mom(%.4f) "
00645                         "pos=%.3f,%.3f,%.3f  4-mom=%.3f,%.3f,%.3f,%.3f\n",
00646                         time,position[0],position[1],position[2],
00647                         fourMomentum[0],fourMomentum[1],fourMomentum[2],
00648                         fourMomentum[3]);
00649 
00650         return true;
00651 }

bool Bunch::init ( )

init() initializes the computation, sizing the grid for the initial bunch, and computing the field for the initial time step. Called after all tracks are added to the bunch, before stepping. Returns false if failure.

References bound_x, bound_y, bound_z, computeField(), getReferencePos4Mom(), BLRunManager::getStepTime(), minActive, poisson, resizeGrid(), runManager, tracks, PoissonConvolve3D::updatePosition(), and verbose.

00654 {
00655         if(minActive < 0)
00656                 minActive = (tracks.size() * (-minActive))/100;
00657 
00658         resizeGrid(true);
00659 
00660         // get position and gamma of reference
00661         G4double time=runManager->getStepTime();
00662         G4ThreeVector refPos;
00663         G4LorentzVector ref4mom;
00664         if(!getReferencePos4Mom(time,refPos,ref4mom)) {
00665             if(verbose >= 3) printf("Bunch::init returning false\n");
00666             return false;
00667         }
00668         G4double gamma=ref4mom.gamma();
00669 
00670         poisson->updatePosition(-bound_x,bound_x,-bound_y,
00671                                         bound_y,-gamma*bound_z,gamma*bound_z);
00672 
00673         bool retval=computeField();
00674 
00675         return retval;
00676 }

bool Bunch::isInBunch ( G4Track * track )

isInBunch() returns true if this track is in the bunch. Tracks are selected by these criteria: same particle as reference within bound_x,bound_y,bound_z of the bunch's reference traj. v/c in the bunch's beam frame less than maxBeta

same particle as reference

within bound_x,bound_y,bound_z of the bunch's reference traj.

v/c in the bunch's beam frame less than maxBeta

References bound_x, bound_y, bound_z, definition, getReferencePos4Mom(), maxBeta, and verbose.

Referenced by addTrack().

00679 {
00680         if(verbose >= 2) {
00681                 G4ThreeVector pos=track->GetPosition();
00682                 G4ThreeVector mom=track->GetMomentum();
00683                 printf("isInBunch? TrackID=%d pos=%.3f,%.3f,%.3f 3-momentum="
00684                         "%.3f,%.3f,%.3f\n", track->GetTrackID(),
00685                         pos[0],pos[1],pos[2],mom[0],mom[1],mom[2]);
00686         }
00687 
00688         // track must be alive
00689         G4TrackStatus trackStatus = track->GetTrackStatus();
00690         if(trackStatus != fAlive && trackStatus != fStopButAlive) {
00691                 if(verbose >= 3) printf("Not in bunch, not alive\n");
00692                 return false;
00693         }
00694 
00695         ///     same particle as reference
00696         if(track->GetDefinition() != definition) {
00697                 if(verbose >= 3) printf("Not in bunch, wrong particle\n");
00698                 return false;
00699         }
00700 
00701         ///     within bound_x,bound_y,bound_z of the bunch's reference traj.
00702         G4ThreeVector refPos;
00703         G4LorentzVector ref4mom;
00704         if(!getReferencePos4Mom(track->GetGlobalTime(),refPos,ref4mom)) {
00705                 if(verbose >= 3) printf("Not in bunch, time too large\n");
00706                 return false;
00707         }
00708         G4ThreeVector p = track->GetPosition() - refPos;
00709         if(fabs(p.x()) > bound_x || fabs(p.y()) > bound_y || 
00710                                                         fabs(p.z()) > bound_z) {
00711                 if(verbose >= 3) printf("Not in bunch, too far away\n");
00712                 return false;
00713         }
00714 
00715         ///     v/c in the bunch's beam frame less than maxBeta
00716         G4ThreeVector boost=ref4mom.findBoostToCM();
00717         G4LorentzVector v;
00718         v.setVectM(track->GetMomentum(),track->GetDefinition()->GetPDGMass());
00719         if(v.boost(boost).beta() > maxBeta) {
00720                 if(verbose >= 3)
00721                         printf("Not in bunch, beta=%.4f > %.4f=maxBeta\n",
00722                                         v.boost(boost).beta(),maxBeta);
00723                 return false;
00724         }
00725 
00726         if(verbose >= 3) printf("Is in bunch\n");
00727         return true;
00728 }

bool Bunch::computeField ( )

All tracks in the bunch are added into the PoissonConvolve3D and the field is computed. returns false if current step time is beyond the range of the bunch's reference.

References approx, bound_x, bound_y, bound_z, charge, Bunch::Approx::computeApprox(), definition, extrapolateFrac, fixedGrid, BLRunManager::getDeltaT(), getReferencePos4Mom(), BLRunManager::getStepTime(), Bunch::Approx::getVector(), goodGuess, minActive, nParticles, poisson, Bunch::Approx::putCharge(), PoissonConvolve3D::putCharge(), Bunch::Approx::reset(), resizeGrid(), runManager, PoissonConvolve3D::setApproximation(), PoissonConvolve3D::solve(), SUPPRESSGRID, tracks, PoissonConvolve3D::updatePosition(), useApproximationOnly, verbose, and PoissonConvolve3D::zeroRhs().

Referenced by init().

00731 {
00732         if(verbose >= 3) printf("Bunch::computeField entered\n");
00733 
00734         // get position and gamma of reference
00735         G4double time=runManager->getStepTime();
00736         G4ThreeVector refPos;
00737         G4LorentzVector ref4mom;
00738         if(!getReferencePos4Mom(time,refPos,ref4mom)) {
00739             printf("***Bunch: Invalid Reference\n");
00740             return false;
00741         }
00742         G4double gamma=ref4mom.gamma();
00743         if(verbose >= 3) 
00744             printf("Bunch::computeField t=%.3f refPos=%.3f,%.3f,%.3f gamma=%.4f\n",
00745                 time,refPos[0],refPos[1],refPos[2],gamma);
00746 
00747 
00748         // resize the grid, if necessary
00749         if(fixedGrid == 0) {
00750                 if(resizeGrid()) {
00751                         poisson->updatePosition(-bound_x,bound_x,-bound_y,
00752                                         bound_y,-gamma*bound_z,gamma*bound_z);
00753                         goodGuess = false;
00754                 }
00755         }
00756 
00757         poisson->zeroRhs();
00758         nParticles = 0;
00759         approx.reset(this);
00760 
00761         int nSuppress=0;
00762         for(unsigned i=0; i<tracks.size(); ++i) {
00763                 G4Track *track = tracks[i];
00764                 if(!track) continue;
00765                 G4TrackStatus trackStatus = track->GetTrackStatus();
00766                 if(trackStatus != fAlive && trackStatus != fStopButAlive) {
00767                         tracks[i] = 0; // erase is expensive for large # tracks
00768                         continue;
00769                 }
00770 
00771                 G4ThreeVector p=track->GetPosition();
00772                 double dt=track->GetGlobalTime() - runManager->getStepTime();
00773                 G4double dtExtrapolate=runManager->getDeltaT()*extrapolateFrac;
00774                 if(fabs(dt) > dtExtrapolate)
00775                         p -= dt * track->GetVelocity();
00776                 p -= refPos;
00777                 if(fabs(p.x()) > bound_x*SUPPRESSGRID ||
00778                    fabs(p.y()) > bound_y*SUPPRESSGRID ||
00779                    fabs(p.z()) > bound_z*SUPPRESSGRID) {
00780                         if(verbose >= 3) printf("spacecharge: suppressed "
00781                                         "p=%.3f,%.3f,%.3f\n",p[0],p[1],p[2]);
00782                         ++nSuppress;
00783                         continue;
00784                 }
00785                 double c=charge*definition->GetPDGCharge();
00786                 if(!poisson->putCharge(p.x(),p.y(),gamma*p.z(),c/epsilon0)) {
00787                         printf("***Not in Poisson grid; p=%.3f,%.3f,%.3f\n",
00788                                                         p.x(),p.y(),p.z());
00789                         continue;
00790                 }
00791                 approx.putCharge(p.x(),p.y(),gamma*p.z(),c/epsilon0);
00792                 ++nParticles;
00793         }
00794 
00795         if(nParticles < minActive) {
00796                 printf("***Bunch: Too few tracks; have %d, require %d\n",
00797                                         nParticles,minActive);
00798                 return false;
00799         }
00800 
00801         if(nSuppress > 0) {
00802                 if(verbose >= 2) 
00803                         printf("Bunch::computeField suppressed %d tracks\n",
00804                                                                 nSuppress);
00805         }
00806 
00807         // Solve Poisson's equation in the beam frame.
00808         approx.computeApprox();
00809         poisson->setApproximation(approx.getVector());
00810         if(useApproximationOnly == 0) {
00811                 if(verbose >= 3) printf("Bunch::computeField calling solve\n");
00812                 poisson->solve();
00813         }
00814 
00815         goodGuess = true;
00816 
00817         if(verbose >= 3) printf("Bunch::computeField returning true\n");
00818         return true;
00819 }

bool Bunch::resizeGrid ( bool force = false )

resizes the grid, based on current tracks. returns true if grid was resized.

References bound_x, bound_y, bound_z, extrapolateFrac, fixedGrid, BLRunManager::getDeltaT(), getReferencePos4Mom(), BLRunManager::getStepTime(), MAXGRID, MINGRID, percentile, runManager, tracks, and verbose.

Referenced by computeField(), and init().

00822 {
00823         if(fixedGrid != 0) return false;
00824 
00825         // get position of reference
00826         G4double time=runManager->getStepTime();
00827         G4ThreeVector refPos;
00828         G4LorentzVector ref4mom;
00829         getReferencePos4Mom(time,refPos,ref4mom);
00830 
00831         // histogram x,y,z for tracks in Bunch
00832         // (histograms extend to twice the bounds.)
00833         const int NHIST=200;
00834         int histX[NHIST], histY[NHIST], histZ[NHIST];
00835         for(int i=0; i<NHIST; ++i)
00836                 histX[i] = histY[i] = histZ[i] = 0;
00837         double dx=2.0*bound_x/NHIST;
00838         double dy=2.0*bound_y/NHIST;
00839         double dz=2.0*bound_z/NHIST;
00840         int n=0;
00841         for(unsigned i=0; i<tracks.size(); ++i) {
00842                 G4Track *track = tracks[i];
00843                 if(track == 0) continue;
00844                 G4TrackStatus trackStatus = track->GetTrackStatus();
00845                 if(trackStatus != fAlive && trackStatus != fStopButAlive) {
00846                         tracks[i] = 0; // erase is expensive for large # tracks
00847                         continue;
00848                 }
00849                 G4ThreeVector p=track->GetPosition();
00850                 double dt=track->GetGlobalTime() - runManager->getStepTime();
00851                 G4double dtExtrapolate=runManager->getDeltaT()*extrapolateFrac;
00852                 if(fabs(dt) > dtExtrapolate)
00853                         p -= dt * track->GetVelocity();
00854                 p -= refPos;
00855                 ++n;
00856                 int ix=(int)(fabs(p.x())/dx);
00857                 if(ix > NHIST-1) ix=NHIST-1;
00858                 ++histX[ix];
00859                 int iy=(int)(fabs(p.y())/dy);
00860                 if(iy > NHIST-1) iy=NHIST-1;
00861                 ++histY[iy];
00862                 int iz=(int)(fabs(p.z())/dz);
00863                 if(iz > NHIST-1) iz=NHIST-1;
00864                 ++histZ[iz];
00865         }
00866         if(n == 0) return false;
00867 
00868         int need=(int)(percentile*n/100.0);
00869         double avg = (MINGRID+MAXGRID)/2.0;
00870 
00871         // find percentile in x, and set bound_x so it's in [MINGRID,MAXGRID]
00872         int t=0, ix=0, iy=0, iz=0;
00873         for(ix=0; ix<NHIST; ++ix) {
00874                 t += histX[ix];
00875                 if(t >= need) break;
00876         };
00877         double new_x = ix*dx;
00878         if(new_x < MINGRID*bound_x || new_x > MAXGRID*bound_x) {
00879                 bound_x = new_x/avg;
00880                 force = true;
00881         }
00882 
00883         // find percentile in y, and set bound_y so it's in [MINGRID,MAXGRID]
00884         t = 0;
00885         for(iy=0; iy<NHIST; ++iy) {
00886                 t += histY[iy];
00887                 if(t >= need) break;
00888         };
00889         double new_y = iy*dy;
00890         if(new_y < MINGRID*bound_y || new_y > MAXGRID*bound_y) {
00891                 bound_y = new_y/avg;
00892                 force = true;
00893         }
00894 
00895         // find percentile in z, and set bound_z so it's in [MINGRID,MAXGRID]
00896         t = 0;
00897         for(iz=0; iz<NHIST; ++iz) {
00898                 t += histZ[iz];
00899                 if(t >= need) break;
00900         };
00901         double new_z = iz*dz;
00902         if(new_z < MINGRID*bound_z || new_z > MAXGRID*bound_z) {
00903                 bound_z = new_z/avg;
00904                 force = true;
00905         }
00906 
00907         if(!force) return false; 
00908 
00909         if(verbose >= 2) printf("Bunch new grid size:%.3f,%.3f,%.3f\n",
00910                                                 bound_x,bound_y,bound_z);
00911         return true;
00912 }

void Bunch::addFieldValue	(	const G4double	point[4],
		G4double	field[6]
	)			const

addFieldValue adds the field for this bunch to field[].

References PoissonConvolve3D::approxE(), BLRunManager::getDeltaT(), PoissonConvolve3D::getE(), getReferencePos4Mom(), BLRunManager::getStepTime(), nParticles, poisson, runManager, useApproximationOnly, and verbose.

00915 {
00916         // return if no particles in bunch
00917         if(nParticles == 0) return;
00918 
00919         // check extrapolation of reference position is OK
00920         G4double stepTime = runManager->getStepTime();
00921         G4double deltaT = runManager->getDeltaT();
00922 
00923         // get position of reference
00924         G4ThreeVector pos;
00925         G4LorentzVector ref4mom;
00926         if(!getReferencePos4Mom(point[3],pos,ref4mom))
00927                 return;
00928         G4ThreeVector boost=ref4mom.findBoostToCM();
00929         double gamma=ref4mom.gamma();
00930         if(verbose >= 3) printf("Bunch::addFieldValue: t=%.3f refPos=%.3f,%.3f,%.3f point=%.3f,%.3f,%.3f\n",point[3],pos[0],pos[1],pos[2],point[0],point[1],point[2]);
00931 
00932         G4ThreeVector Ebeam;
00933         if(useApproximationOnly)
00934                 Ebeam = poisson->approxE(point[0]-pos[0],
00935                                 point[1]-pos[1],gamma*(point[2]-pos[2]));
00936         else
00937                 Ebeam = poisson->getE(point[0]-pos[0],
00938                                 point[1]-pos[1],gamma*(point[2]-pos[2]));
00939 
00940         // boost of EM fields from Jefimenko, AJP _64_(5), p618 (1996), 
00941         // eq12-17 with xyz->zxy and B'=0; boost[] is beta[]=v[]/c_light;
00942         // There is an additional minus sign in boost[] from findBoostToCM().
00943         // Reference is parallel to the z axis.
00944         double f[6];
00945         f[0] = gamma*boost[2]*Ebeam[1]/c_light;         // Bx
00946         f[1] = -gamma*boost[2]*Ebeam[0]/c_light;        // By
00947         f[2] = 0.0;                                     // Bz
00948         f[3] = gamma*Ebeam[0];          // Ex
00949         f[4] = gamma*Ebeam[1];          // Ey
00950         f[5] = Ebeam[2];                // Ez
00951 
00952         if(verbose >= 3) 
00953                 printf("Bunch::addFieldValue(%.3f,%.3f,%.3f,%.3f)\n"
00954                         "    pos=%.3f,%.3f,%.3f boost=%.3f,%.3f,%.3f "
00955                         "gamma=%.3f  Ebeam=%.3f,%.3f,%.3f\n"
00956                         "    B=%.3f,%.3f,%.3f  E=%.3f,%.3f,%.3f\n",
00957                         point[0]/mm,point[1]/mm,point[2]/mm,point[3]/ns,
00958                         pos[0]/mm,pos[1]/mm,pos[2]/mm,
00959                         boost[0],boost[1],boost[2],gamma,
00960                         Ebeam[0]/(megavolt/meter),
00961                         Ebeam[1]/(megavolt/meter),
00962                         Ebeam[2]/(megavolt/meter),
00963                         f[0]/tesla,f[1]/tesla,f[2]/tesla,
00964                         f[3]/(megavolt/meter),f[4]/(megavolt/meter),
00965                         f[5]/(megavolt/meter));
00966 
00967         field[0] += f[0];
00968         field[1] += f[1];
00969         field[2] += f[2];
00970         field[3] += f[3];
00971         field[4] += f[4];
00972         field[5] += f[5];
00973 }

void Bunch::testEfield ( ) [static]

testEfield() will test the units and value of the E field in the beam frame for a point charge.

References BLAssert, PoissonConvolve3D::defaultApproximation(), E, PoissonConvolve3D::getE(), PoissonConvolve3D::putCharge(), PoissonConvolve3D::solve(), and PoissonConvolve3D::zeroRhs().

Referenced by BLCMDspacecharge::command().

00976 {
00977         /* UNITS TEST
00978            E = k*Q/r/r,  k=8.987551787E9 (Newton*meter*meter/Coulomb/Coulomb)
00979            And remember that 1 volt/meter == 1 Newton/Coulomb.
00980            So a 1 Coulomb charge has E=8.99E3 MV/m at a distance of 1 meter.
00981            1 coulomb = 6.24150E+18 e+.
00982         */
00983         const double HW=100.0*mm; // half-width of grid
00984         PoissonConvolve3D *p = new PoissonConvolve3D(33,33,33,-HW,HW,-HW,HW,-HW,HW);
00985         p->zeroRhs();
00986         p->putCharge(0.0,0.0,0.0,6.24150E+18/epsilon0);
00987         p->defaultApproximation();
00988         p->solve();
00989         G4ThreeVector E = p->getE(1.0*meter,0.0,0.0);
00990         printf("spacecharge: testEfield: E field should be %.2f MV/m: %.2f\n",
00991                         8987.55, E.x()/(megavolt/meter));
00992         BLAssert(fabs(E.x()/(megavolt/meter)-8987.55) < 0.1 && \
00993                 fabs(E.y()/(megavolt/meter)) < 0.1 && \
00994                 fabs(E.z()/(megavolt/meter)) < 0.1);
00995 /***
00996         // check for continuity at the edge of the grid.
00997         printf("#x Phi Ex\n");
00998         for(double x=0.9*HW; x<1.1*HW; x+=0.01*HW) {
00999                 E = p->getE(x,0.0,0.0);
01000                 printf("%.2f %.5g %.1f\n",x,p->phi(x,0.0,0.0),
01001                                 E.x()/(megavolt/meter));
01002         }
01003 ***/
01004         delete p;
01005 }

void Bunch::testBfield ( ) [static]

testBfield() will test the B field in the lab frame for a bunch that is a line of moving charges, 1000 Amps.

References BLAssert, E, PoissonConvolve3D::getE(), PoissonConvolve3D::putCharge(), PoissonConvolve3D::setApproximation(), PoissonConvolve3D::solve(), and PoissonConvolve3D::zeroRhs().

Referenced by BLCMDspacecharge::command().

01008 {
01009         /* Magnetic field test
01010            A uniform current of 1000 Amps, at radius 0.1 m, has B=0.002 T.
01011            1000 Amps = 0.001C charges 1 cm apart moving at 10000 m/s.
01012         */
01013         double c=0.001*coulomb;
01014         double beta=10000*meter/second/c_light;
01015         double gamma=1.0/sqrt(1.0-beta*beta);
01016         double d=1.0*cm;
01017         PoissonConvolve3D *p = new PoissonConvolve3D(129,129,129,-1000.0,1000.0,
01018                                 -1000.0,1000.0,-gamma*1000.0,gamma*1000.0);
01019         p->zeroRhs();
01020         std::vector<Charge> v;
01021         for(double z=-50.0*d; z<50.01*d; z+=d) {
01022                 v.push_back(Charge(0.0,0.0,gamma*z,c/epsilon0));
01023                 p->putCharge(0.0,0.0,gamma*z,c/epsilon0);
01024         }
01025         p->setApproximation(v);
01026         p->solve();
01027         G4ThreeVector E = p->getE(0.1*meter,0.0,0.0);
01028         double By = gamma*beta*E[0]/c_light;    // By
01029         printf("spacecharge: testBfield: By should be 0.002 tesla: %.6f\n",
01030                                                                 By/tesla);
01031         BLAssert(By-0.002*tesla < 0.00005*tesla);
01032 }

static void Bunch::setVerbose ( int v ) [inline, static]

setVerbose() sets the debug print level

References verbose.

Referenced by BLCMDspacecharge::command().

00239 { verbose=v; }

static void Bunch::setUseApproximationOnly ( int v ) [inline, static]

setUseApproximationOnly() sets it

References useApproximationOnly.

Referenced by BLCMDspacecharge::command().

00242 { useApproximationOnly=v; }

Friends And Related Function Documentation

friend class Approx [friend]

Member Data Documentation

const double Bunch::extrapolateFrac = 0.01 [static, private]

Referenced by computeField(), and resizeGrid().

const double Bunch::minDeltaT = 1.0e-6*ns [static, private]

Referenced by addReferencePoint(), and getReferencePos4Mom().

std::vector<Point> Bunch::point [private]

Referenced by addReferencePoint(), getReferencePos4Mom(), and nTrajectoryPoints().

unsigned Bunch::index [private]

Referenced by Bunch(), and getReferencePos4Mom().

double Bunch::prevTime [private]

Referenced by addReferencePoint(), and Bunch().

G4ParticleDefinition* Bunch::definition [private]

Referenced by Bunch(), computeField(), getDefinition(), getReferencePos4Mom(), and isInBunch().

double Bunch::charge [private]

Referenced by Bunch(), and computeField().

double Bunch::maxBeta [private]

Referenced by Bunch(), and isInBunch().

double Bunch::bound_x [private]

Referenced by Bunch(), computeField(), getBound(), init(), isInBunch(), Bunch::Approx::reset(), and resizeGrid().

double Bunch::bound_y [private]

Referenced by Bunch(), computeField(), getBound(), init(), isInBunch(), Bunch::Approx::reset(), and resizeGrid().

double Bunch::bound_z [private]

Referenced by Bunch(), computeField(), getBound(), init(), isInBunch(), Bunch::Approx::reset(), and resizeGrid().

int Bunch::nParticles [private]

Referenced by addFieldValue(), Bunch(), computeField(), and getNParticles().

double Bunch::percentile [private]

Referenced by Bunch(), and resizeGrid().

int Bunch::minActive [private]

Referenced by Bunch(), computeField(), getMinActive(), and init().

PoissonConvolve3D* Bunch::poisson [private]

Referenced by addFieldValue(), Bunch(), computeField(), and init().

bool Bunch::goodGuess [private]

Referenced by Bunch(), and computeField().

bool Bunch::fixedGrid [private]

Referenced by Bunch(), computeField(), and resizeGrid().

Approx Bunch::approx [private]

Referenced by Bunch(), computeField(), and getNapprox().

std::vector<G4Track*> Bunch::tracks [private]

Referenced by addTrack(), computeField(), init(), resizeGrid(), and size().

BLRunManager* Bunch::runManager [private]

Referenced by addFieldValue(), Bunch(), computeField(), init(), and resizeGrid().

int Bunch::verbose = 0 [static, private]

Referenced by addFieldValue(), computeField(), getReferencePos4Mom(), init(), isInBunch(), resizeGrid(), and setVerbose().

int Bunch::useApproximationOnly = 0 [static, private]

Referenced by addFieldValue(), computeField(), and setUseApproximationOnly().

The documentation for this class was generated from the following file:

BLCMDspacecharge.cc


Public Member Functions
	Bunch (G4ParticleDefinition *def, int nx, int ny, int nz, double dx, double dy, double dz, double _charge, double _maxBeta, bool _fixedGrid, double _percentile, int nxApprox, int nyApprox, int nzApprox, int _minActive)
void	addReferencePoint (const G4Track *track)
	add a point to the reference trajectory.
void	addReferencePoint (double t, G4ThreeVector pos, G4ThreeVector mom)
	add a point to the reference trajectory.
int	nTrajectoryPoints ()
	nTrajectoryPoints() returns the # points in the reference trajectory
int	getNParticles () const
	getNParticles() returns # particles in bunch
G4ParticleDefinition *	getDefinition ()
	getParticleDefinition()
bool	addTrack (G4Track *t)
	Add a track to the Bunch. returns false if this track is not part of this bunch.
int	size ()
int	getNapprox ()
int	getMinActive ()
G4ThreeVector	getBound ()
bool	getReferencePos4Mom (double time, G4ThreeVector &position, G4LorentzVector &fourMomentum) const
	Interpolate in the point[] vector to get position and 4-momentum of the reference for this bunch. Returns true if OK, false if time is after the valid range for this reference. If time is before the reference was created, its first two trajectory points are linearly extrapolated backwards. Because it must search for the enclosing points of the trajectory, this is most efficient when time does not jump very much (it starts the search at the place in the trajectory of the previous call). (Note: the const is a lie, but is required by BLElementField; only index is modified, which is essentially benign.).
bool	init ()
	init() initializes the computation, sizing the grid for the initial bunch, and computing the field for the initial time step. Called after all tracks are added to the bunch, before stepping. Returns false if failure.
bool	isInBunch (G4Track *track)
	isInBunch() returns true if this track is in the bunch. Tracks are selected by these criteria: same particle as reference within bound_x,bound_y,bound_z of the bunch's reference traj. v/c in the bunch's beam frame less than maxBeta
bool	computeField ()
	All tracks in the bunch are added into the PoissonConvolve3D and the field is computed. returns false if current step time is beyond the range of the bunch's reference.
bool	resizeGrid (bool force=false)
	resizes the grid, based on current tracks. returns true if grid was resized.
void	addFieldValue (const G4double point[4], G4double field[6]) const
	addFieldValue adds the field for this bunch to field[].
Static Public Member Functions
static void	testEfield ()
	testEfield() will test the units and value of the E field in the beam frame for a point charge.
static void	testBfield ()
	testBfield() will test the B field in the lab frame for a bunch that is a line of moving charges, 1000 Amps.
static void	setVerbose (int v)
	setVerbose() sets the debug print level
static void	setUseApproximationOnly (int v)
	setUseApproximationOnly() sets it
Private Attributes
std::vector< Point >	point
unsigned	index
double	prevTime
G4ParticleDefinition *	definition
double	charge
double	maxBeta
double	bound_x
double	bound_y
double	bound_z
int	nParticles
double	percentile
int	minActive
PoissonConvolve3D *	poisson
bool	goodGuess
bool	fixedGrid
Approx	approx
std::vector< G4Track * >	tracks
BLRunManager *	runManager
Static Private Attributes
static const double	extrapolateFrac = 0.01
static const double	minDeltaT = 1.0e-6*ns
static int	verbose = 0
static int	useApproximationOnly = 0
Friends
class	Approx
Classes
class	Approx
	class Approx does a 3-d histogram of the charge to generate an approximation for PoissonConvolve3D. Within each bin (voxel) it use the mean position and the total charge. More...
struct	Point