MacroParticle Class Reference

Detailed Description

class MacroParticle represents a spherical distribution of charge with uniform density and specified radius, following a specified trajectory.

APPROXIMATION: the nonzero radius is used only to prevent infinities for field computations inside the macro-particle. The speed-of-light delay within the macro-particle is neglected -- IOW the entire macro particle is treated as if it were located at its center, but the charge of the macro-particle is scaled by r/radius when r<radius.

tmin and tmax default to -infinity and +infinity; they are used to limit the field of the macroparticle when it is created or destroyed (i.e. secondaries should set tmin, and particles that interact or decay should set tmax when they are killed).

Public Member Functions

MacroParticle (G4double _charge, G4double _radius, G4int _K=1)

Constructor -- no trajectory point.

MacroParticle (const G4Track *track, G4double _charge, G4double _radius, G4int _K=1)

constructor from track, includes trajectory point.

void setTmin (G4double v)

setTmin() sets tmin (time of creation).

void setTmax (G4double v)

setTmax() sets tmax (time of destruction).

void addPoint (const G4Track *track)

addPoint() adds a trajectory point to the MacroParticle. Points must be added in order of increasing time (e.g. in UserSteppingAction()).

void addPoint (const G4double pos[4], const G4double v[3])

addPoint() adds a trajectory point to the MacroParticle. Points must be added in order of increasing time (e.g. in UserSteppingAction()). pos={x,y,z,t} v={Vx,Vy,Vz}

G4double getMostRecentTime ()

getMostRecentTime() returns the time of the most recent point in trajectory[].

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

computeField() computes the E and B fields due to this macro particle at the specified point. point[] = {x,y,z,t}, field[]={Bx,By,Bz,Ex,Ey,Ez}; global coordinates. The MacroParticle must have at least one TrajectoryPoint.

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

addField() is just like computeField(), except the value for this MacroParticle is added to the current value of field[].

void pruneTrajectory ()

pruneTrajectory() will remove unneeded entries in trajectory[] It requires 100 calls to addfield() between prunes, leaves 10 extra entries (just in case of error in minIndex), and always leaves at least 20 entries in trajectory[].

Static Public Member Functions

static void testUnits ()

testUnits will test the E and B field units. User validates by looking at stdout (correct and computed values are printed).

Static Public Attributes

static int verbose = 0

Private Attributes

G4double charge

G4double radius

G4int K

G4double tmin

G4double tmax

G4int minIndex

G4int nCalls

std::vector< TrajectoryPoint > trajectory

Constructor & Destructor Documentation

MacroParticle::MacroParticle	(	G4double	_charge,
		G4double	_radius,
		G4int	_K = `1`
	)			`[inline]`

Constructor -- no trajectory point.

References charge, K, minIndex, nCalls, radius, tmax, and tmin.

Referenced by testUnits().

00080                                                                       :
00081                                                                 trajectory()
00082                 { charge=_charge; radius=_radius; K=_K;  tmin=-DBL_MAX;
00083                   tmax=DBL_MAX; minIndex=0; nCalls=0; }

MacroParticle::MacroParticle	(	const G4Track *	track,
		G4double	_charge,
		G4double	_radius,
		G4int	_K = `1`
	)			`[inline]`

constructor from track, includes trajectory point.

References addPoint(), charge, K, minIndex, nCalls, radius, tmax, and tmin.

00087                                                             : trajectory() {
00088                 charge = track->GetDefinition()->GetPDGCharge() * _charge;
00089                 radius=_radius;
00090                 K=_K;
00091                 tmin=-DBL_MAX;
00092                 tmax=DBL_MAX;
00093                 minIndex=0;
00094                 nCalls=0;
00095                 addPoint(track);
00096         }

Member Function Documentation

void MacroParticle::setTmin ( G4double v ) [inline]

setTmin() sets tmin (time of creation).

References tmin.

00099 { tmin = v; }

void MacroParticle::setTmax ( G4double v ) [inline]

setTmax() sets tmax (time of destruction).

References tmax.

00102 { tmax = v; }

void MacroParticle::addPoint ( const G4Track * track ) [inline]

addPoint() adds a trajectory point to the MacroParticle. Points must be added in order of increasing time (e.g. in UserSteppingAction()).

Referenced by MacroParticle(), and testUnits().

00107                                             {
00108                 G4double pos[4], vel[3];
00109                 G4ThreeVector p=track->GetPosition();
00110                 pos[0] = p[0];
00111                 pos[1] = p[1];
00112                 pos[2] = p[2];
00113                 pos[3] = track->GetGlobalTime();
00114                 p = track->GetMomentumDirection();
00115                 G4double speed = track->GetVelocity();
00116                 vel[0] = p[0] * speed;
00117                 vel[1] = p[1] * speed;
00118                 vel[2] = p[2] * speed;
00119                 addPoint(pos,vel);
00120         }

void MacroParticle::addPoint	(	const G4double	pos[4],
		const G4double	v[3]
	)			`[inline]`

addPoint() adds a trajectory point to the MacroParticle. Points must be added in order of increasing time (e.g. in UserSteppingAction()). pos={x,y,z,t} v={Vx,Vy,Vz}

References trajectory.

00127                 { trajectory.push_back(TrajectoryPoint(pos,v)); }

G4double MacroParticle::getMostRecentTime ( ) [inline]

getMostRecentTime() returns the time of the most recent point in trajectory[].

References BLAssert, and trajectory.

00131                                      {
00132                 BLAssert(trajectory.size() >= 1);
00133                 return trajectory[trajectory.size()-1].pos[3];
00134         }

void MacroParticle::computeField	(	const G4double	point[4],
		G4double	field[6]
	)			`[inline]`

computeField() computes the E and B fields due to this macro particle at the specified point. point[] = {x,y,z,t}, field[]={Bx,By,Bz,Ex,Ey,Ez}; global coordinates. The MacroParticle must have at least one TrajectoryPoint.

References addField().

Referenced by testUnits().

00140                                                                       {
00141                 field[0]=field[1]=field[2]=field[3]=field[4]=field[5]=0.0;
00142                 addField(point,field);
00143         }

void MacroParticle::addField	(	const G4double	point[4],
		G4double	field[6]
	)

addField() is just like computeField(), except the value for this MacroParticle is added to the current value of field[].

References BLAssert, charge, E, K, minIndex, nCalls, radius, tmax, tmin, trajectory, and verbose.

Referenced by computeField().

00671 {
00672         // p[] is the point at which the field is to be evaluated.
00673         // q[] is a point on the trajectory, v[] is its velocity at q.
00674         // pq[] = p - q.
00675         double q[4], v[3], pq[4];
00676         pq[0] = pq[1] = pq[2] = pq[3] = 0.0;
00677 
00678         if(charge == 0.0) return;
00679         BLAssert(trajectory.size() > 0);
00680 
00681         /* find the TrajectoryPoint closest before the lightcone from p
00682            Use first point (i==0) if before trajectory began. */
00683         int i;
00684         for(i=(int)trajectory.size()-1; i>=0; --i) {
00685                 const double *qq = trajectory[i].pos;
00686                 pq[0] = p[0] - qq[0];
00687                 pq[1] = p[1] - qq[1];
00688                 pq[2] = p[2] - qq[2];
00689                 pq[3] = p[3] - qq[3];
00690                 double d = sqrt(pq[0]*pq[0]+pq[1]*pq[1]+pq[2]*pq[2]);
00691                 if(i == 0 || d/c_light < pq[3])
00692                         break;
00693         }
00694         q[0] = trajectory[i].pos[0];
00695         q[1] = trajectory[i].pos[1];
00696         q[2] = trajectory[i].pos[2];
00697         q[3] = trajectory[i].pos[3];
00698         v[0] = trajectory[i].v[0];
00699         v[1] = trajectory[i].v[1];
00700         v[2] = trajectory[i].v[2];
00701 
00702         // update trajectory pruning info
00703         if(i < minIndex) minIndex = i;
00704         ++nCalls;
00705 
00706         /* now solve for the time of the intersection of the lightcone from p 
00707            with the trajectory, assuming v is constant; t is the offset 
00708            from q[3]. The other root is the unphysical advanced field. */
00709         double v2 = v[0]*v[0]+v[1]*v[1]+v[2]*v[2];
00710         double c2 = c_light*c_light;
00711         double a = c2 - v2;
00712         double b = 2.0*(pq[0]*v[0]+pq[1]*v[1]+pq[2]*v[2]-c2*pq[3]);
00713         double c = c2*pq[3]*pq[3]-(pq[0]*pq[0]+pq[1]*pq[1]+pq[2]*pq[2]);
00714         double discriminant=b*b-4.0*a*c;
00715         // if discriminant ~ 0, p is deep inside this MacroParticle -> no field
00716         if(discriminant < 0.001) {
00717                 if(verbose >= 3) printf("addField: discriminant < 0.001\n");
00718                 return;
00719         }
00720         // note that a>0 so -sqrt(discriminant) gives the retarded root
00721         double t = (-b - sqrt(discriminant))/(2.0*a);
00722 
00723         /* update q[] and pq[] to the actual intersection */
00724         q[0] += t*v[0];
00725         q[1] += t*v[1];
00726         q[2] += t*v[2];
00727         q[3] += t;
00728         pq[0] = p[0] - q[0];
00729         pq[1] = p[1] - q[1];
00730         pq[2] = p[2] - q[2];
00731         pq[3] = p[3] - q[3];
00732 
00733         // apply creation/destruction limits.
00734         if(q[3] < tmin || q[3] > tmax) {
00735                 if(verbose >= 3) printf("addField: t=%.4f tmin=%.4f tmax=%.4f\n",
00736                                         q[3],tmin,tmax);
00737                 return;
00738         }
00739 
00740         // r = distance from q to p
00741         double r = sqrt(pq[0]*pq[0]+pq[1]*pq[1]+pq[2]*pq[2]);
00742         // require accuracy in the lightcone: 1 ppm or 0.01 micron
00743         double tolerance=(r>10.0*mm ? r/1000000.0 : 0.00001*mm);
00744         double light_cone_error=fabs(r-c_light*pq[3]);
00745         BLAssert(light_cone_error < tolerance);
00746 
00747         /* effective charge */
00748         double f = 1.0;
00749         if(r < radius) {
00750                 double a=r/radius;
00751                 switch(K) {
00752                 case 0: f = 1.0;        break;
00753                 case 1: f = a;          break;
00754                 case 2: f = a*a;        break;
00755                 case 3: f = a*a*a;      break;
00756                 case 4: f = a*a*a*a;    break;
00757                 }
00758         }
00759         if(f < 1.0e-4) return;
00760         double ec = f * charge;
00761 
00762         /* get E and B fields  - Landau and Lifshitz, _Classical_Theory_of_
00763            _Fields_, eq. 63.8-9 p 162. Ignore radiation term.  
00764            Units done by hand: the usual Coulomb constant is divided by 1000,
00765            and B is divided by c_light.
00766         */
00767         double k = r - (pq[0]*v[0]+pq[1]*v[1]+pq[2]*v[2])/c_light;
00768         double E[3];
00769         k = ec * e_SI * 8.987551787E6 * (1.0-v2/c_light/c_light)/k/k/k;
00770         E[0] = k * (pq[0]-v[0]*r/c_light);
00771         E[1] = k * (pq[1]-v[1]*r/c_light);
00772         E[2] = k * (pq[2]-v[2]*r/c_light);
00773         field[0] += (pq[1]*E[2] - pq[2]*E[1])/r/c_light;
00774         field[1] += (pq[2]*E[0] - pq[0]*E[2])/r/c_light;
00775         field[2] += (pq[0]*E[1] - pq[1]*E[0])/r/c_light;
00776         field[3] += E[0];
00777         field[4] += E[1];
00778         field[5] += E[2];
00779 }

void MacroParticle::pruneTrajectory ( )

pruneTrajectory() will remove unneeded entries in trajectory[] It requires 100 calls to addfield() between prunes, leaves 10 extra entries (just in case of error in minIndex), and always leaves at least 20 entries in trajectory[].

References minIndex, nCalls, trajectory, and verbose.

00782 {
00783         const int safety=10;    // extra entries left in trajectory[]
00784         const int minKeep=20;   // min # entries to keep
00785         const int minCalls=100; // min # calls to addfield() between prunes
00786 
00787         if(nCalls < minCalls || minIndex <= safety+1) return;
00788 
00789         minIndex -= safety;
00790         if(trajectory.size()-(unsigned)minIndex < (unsigned)minKeep) return;
00791 
00792         for(int i=0; i<minIndex; ++i)
00793                 trajectory.erase(trajectory.begin());
00794 
00795         if(verbose >= 3) printf("MacroParticle:prune removed %d newSize=%ld tMin=%.6f\n",
00796                         minIndex,(long)trajectory.size(),trajectory[0].pos[3]);
00797 
00798         minIndex = trajectory.size();
00799         nCalls = 0;
00800 }

void MacroParticle::testUnits ( ) [static]

testUnits will test the E and B field units. User validates by looking at stdout (correct and computed values are printed).

References addPoint(), charge, computeField(), and MacroParticle().

Referenced by BLCMDspacechargelw::command().

00803 {
00804         printf("MacroParticle::testUnits()\n");
00805         G4double pos[4]={0,0,0,0};
00806         G4double v[3]={0,0,0};
00807         G4double field[6]; 
00808 
00809         /* UNITS TEST
00810            E = k*Q/r/r,  k=8.987551787E9 (Newton*meter*meter/Coulomb/Coulomb)
00811            And remember that 1 volt/meter == 1 Newton/Coulomb.
00812            So a 1 Coulomb charge has E=8.99E3 MV/m at a distance of 1 meter.
00813         */
00814         MacroParticle q(coulomb,0.1*mm);
00815         pos[0]=pos[1]=pos[2]=pos[3]=0.0;
00816         v[0]=v[1]=v[2]=0.0;
00817         q.addPoint(pos,v);
00818         pos[0] = 1.0*meter;
00819         q.computeField(pos,field);
00820         printf("E field test: Ex should be 8988 MV/m   pos=%.0f,%.0f,%.0f  E=%.4g,%.4g,%.4g MV/m\n",
00821                 pos[0]/mm,pos[1]/mm,pos[2]/mm,
00822                 field[3]/(megavolt/meter),field[4]/(megavolt/meter),
00823                 field[5]/(megavolt/meter));
00824 
00825         pos[0] = 0.0;
00826         pos[1] = 2.0*meter;
00827         q.computeField(pos,field);
00828         printf("E field test: Ey should be 2247 MV/m   pos=%.0f,%.0f,%.0f  E=%.4g,%.4g,%.4g MV/m\n",
00829                 pos[0]/mm,pos[1]/mm,pos[2]/mm,
00830                 field[3]/(megavolt/meter),field[4]/(megavolt/meter),
00831                 field[5]/(megavolt/meter));
00832 
00833 
00834         pos[1] = 0.0;
00835         pos[2] = -1.0*meter;
00836         q.computeField(pos,field);
00837         printf("E field test: Ez should be -8988 MV/m   pos=%.0f,%.0f,%.0f  E=%.4g,%.4g,%.4g MV/m\n",
00838                 pos[0]/mm,pos[1]/mm,pos[2]/mm,
00839                 field[3]/(megavolt/meter),field[4]/(megavolt/meter),
00840                 field[5]/(megavolt/meter));
00841 
00842 
00843         /* Magnetic field test
00844            A uniform current of 1000 Amps, at radius 0.1 m, has B=0.002 T.
00845            1000 Amps = 0.001C charges 1 cm apart moving at 10000 m/s.
00846         */
00847         int N=8001;             // # charges
00848         double dz=10*mm;        // their interval along z
00849         double charge = 0.001*coulomb;  // their charge
00850         double speed=10000*meter/second; // their speed along z
00851         std::vector<MacroParticle> array;
00852         for(int i=0; i<N; ++i) {
00853                 array.push_back(MacroParticle(charge,0.1*mm));
00854                 pos[0]=pos[1]=pos[2]=pos[3]=0.0;
00855                 pos[2] = (i-N/2) * dz;
00856                 v[0]=v[1]=0.0;
00857                 v[2] = speed;
00858                 array[i].addPoint(pos,v);
00859         }
00860         field[0]=field[1]=field[2]=field[3]=field[4]=field[5]=0.0;
00861         pos[0]=pos[1]=pos[2]=pos[3]=0.0;
00862         pos[0] = 0.1*meter;
00863         for(int i=0; i<N; ++i) {
00864                 double f[6];
00865                 array[i].computeField(pos,f);
00866                 for(int j=0; j<6; ++j) field[j] += f[j];
00867         }
00868         printf("B field test: By should be 0.002 T  pos=%.0f,%.0f,%.0f  B=%.6g,%.6g,%.6g T\n",
00869                 pos[0]/mm,pos[1]/mm,pos[2]/mm,
00870                 field[0]/tesla,field[1]/tesla,field[2]/tesla);
00871         field[0]=field[1]=field[2]=field[3]=field[4]=field[5]=0.0;
00872         pos[0]=pos[1]=pos[2]=pos[3]=0.0;
00873         pos[1] = 0.2*meter; // further away
00874         for(int i=0; i<N; ++i) {
00875                 double f[6];
00876                 array[i].computeField(pos,f);
00877                 for(int j=0; j<6; ++j) field[j] += f[j];
00878         }
00879         printf("B field test: Bx should be -0.001 T  pos=%.0f,%.0f,%.0f  B=%.6g,%.6g,%.6g T\n",
00880                 pos[0]/mm,pos[1]/mm,pos[2]/mm,
00881                 field[0]/tesla,field[1]/tesla,field[2]/tesla);
00882 }

Member Data Documentation

G4double MacroParticle::charge [private]

Referenced by addField(), MacroParticle(), and testUnits().

G4double MacroParticle::radius [private]

Referenced by addField(), and MacroParticle().

G4int MacroParticle::K [private]

Referenced by addField(), and MacroParticle().

G4double MacroParticle::tmin [private]

Referenced by addField(), MacroParticle(), and setTmin().

G4double MacroParticle::tmax [private]

Referenced by addField(), MacroParticle(), and setTmax().

G4int MacroParticle::minIndex [private]

Referenced by addField(), MacroParticle(), and pruneTrajectory().

G4int MacroParticle::nCalls [private]

Referenced by addField(), MacroParticle(), and pruneTrajectory().

std::vector<TrajectoryPoint> MacroParticle::trajectory [private]

Referenced by addField(), addPoint(), getMostRecentTime(), and pruneTrajectory().

int MacroParticle::verbose = 0 [static]

Referenced by addField(), BLCMDspacechargelw::beginCollectiveTracking(), BLCMDspacechargelw::command(), and pruneTrajectory().

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

BLCMDspacechargelw.cc


Public Member Functions
	MacroParticle (G4double _charge, G4double _radius, G4int _K=1)
	Constructor -- no trajectory point.
	MacroParticle (const G4Track *track, G4double _charge, G4double _radius, G4int _K=1)
	constructor from track, includes trajectory point.
void	setTmin (G4double v)
	setTmin() sets tmin (time of creation).
void	setTmax (G4double v)
	setTmax() sets tmax (time of destruction).
void	addPoint (const G4Track *track)
	addPoint() adds a trajectory point to the MacroParticle. Points must be added in order of increasing time (e.g. in UserSteppingAction()).
void	addPoint (const G4double pos[4], const G4double v[3])
	addPoint() adds a trajectory point to the MacroParticle. Points must be added in order of increasing time (e.g. in UserSteppingAction()). pos={x,y,z,t} v={Vx,Vy,Vz}
G4double	getMostRecentTime ()
	getMostRecentTime() returns the time of the most recent point in trajectory[].
void	computeField (const G4double point[4], G4double field[6])
	computeField() computes the E and B fields due to this macro particle at the specified point. point[] = {x,y,z,t}, field[]={Bx,By,Bz,Ex,Ey,Ez}; global coordinates. The MacroParticle must have at least one TrajectoryPoint.
void	addField (const G4double point[4], G4double field[6])
	addField() is just like computeField(), except the value for this MacroParticle is added to the current value of field[].
void	pruneTrajectory ()
	pruneTrajectory() will remove unneeded entries in trajectory[] It requires 100 calls to addfield() between prunes, leaves 10 extra entries (just in case of error in minIndex), and always leaves at least 20 entries in trajectory[].
Static Public Member Functions
static void	testUnits ()
	testUnits will test the E and B field units. User validates by looking at stdout (correct and computed values are printed).
Static Public Attributes
static int	verbose = 0
Private Attributes
G4double	charge
G4double	radius
G4int	K
G4double	tmin
G4double	tmax
G4int	minIndex
G4int	nCalls
std::vector< TrajectoryPoint >	trajectory