BLCMDbeam Class Reference

Inheritance diagram for BLCMDbeam:

Detailed Description

class BLCMDbeam implements the beam command.


Public Member Functions
	BLCMDbeam ()
	Constructor.
virtual G4String	commandName ()
	commandName() returns "beam".
virtual int	command (BLArgumentVector &argv, BLArgumentMap &namedArgs)
	command() implements the beam command.
virtual void	printBeam ()
	printBeam() will print a description.
virtual void	defineNamedArgs ()
	defineNamedArgs() defines the named arguments for this command.
virtual int	getNEvents () const
	getNEvents() returns the # events to process.
virtual void	init ()
	init() will initialize internal variables.
virtual bool	generateReferenceParticle (G4Event *event)
	generateReferenceParticle() generates the reference particle.
virtual bool	nextBeamEvent (G4Event *event)
	nextBeamEvent() generates the next beam event.
virtual void	summary ()
	summary() will print a summary, if necessary.
Private Types
enum	BeamType { NONE, GAUSSIAN, RECTANGULAR, ASCII, ELLIPSE, ROOT }
Private Member Functions
G4double	myrand (G4double mean, G4double sigma)
Private Attributes
BeamType	type
G4String	particle
G4int	eventsGenerated
G4int	nEvents
G4int	eventID
G4int	firstEvent
G4int	lastEvent
G4double	beamX
G4double	beamY
G4double	beamZ
G4double	maxR
G4String	rotation
G4RotationMatrix *	rotationMatrix
G4ThreeVector	position
G4int	renumber
G4double	weight
G4int	secondaryTrackID
G4double	meanMomentum
G4double	sigmaX
G4double	sigmaY
G4double	sigmaZ
G4double	sigmaXp
G4double	sigmaYp
G4double	sigmaP
G4double	sigmaT
G4double	sigmaE
G4double	meanXp
G4double	meanYp
G4double	meanT
G4double	beamHeight
G4double	beamWidth
G4String	filename
G4String	directory
G4int	uid
G4String	name
BLNTuple *	ntuple
double *	data
G4String	format
G4int	index
G4ParticleGun *	particleGun
G4ParticleDefinition *	particleDefinition
BLTrackFile *	trackFile
G4int	nvar
G4int	prevEventID

Member Enumeration Documentation

enum BLCMDbeam::BeamType [private]

Enumerator:

NONE
GAUSSIAN
RECTANGULAR
ASCII
ELLIPSE
ROOT

00046 { NONE, GAUSSIAN, RECTANGULAR, ASCII, ELLIPSE, ROOT };

Constructor & Destructor Documentation

BLCMDbeam::BLCMDbeam ( )

Constructor.

References ALL_EVENTS, beamHeight, beamWidth, beamX, beamY, beamZ, BLCMDTYPE_BEAM, data, directory, eventID, eventsGenerated, filename, firstEvent, format, index, lastEvent, maxR, meanMomentum, meanT, meanXp, meanYp, name, nEvents, NONE, ntuple, nvar, particle, particleDefinition, particleGun, prevEventID, BLCommand::registerCommand(), renumber, rotationMatrix, secondaryTrackID, BLCommand::setDescription(), BLCommand::setSynopsis(), sigmaE, sigmaP, sigmaT, sigmaX, sigmaXp, sigmaY, sigmaYp, sigmaZ, trackFile, type, uid, UNDEFINED, and weight.

Referenced by command().

00138                      : rotation() , position()
00139 {
00140         registerCommand(BLCMDTYPE_BEAM);
00141         setSynopsis("Define the Beam.");
00142         setDescription("The beam command is: beam type arg1=v1 ...\n\n"
00143                 "Types are: gaussian, rectangular, ellipse, ascii, and root.\n\n"
00144                 "Gaussian beams are randomly generated to emanate from "
00145                 "beamX,beamY,beamZ with the given "
00146                 "sigmas; negative sigma means flat with |sigma| as "
00147                 "halfwidth.\n\n"
00148                 "Rectangular beams are randomly generated to emanate "
00149                 "from the rectangle beamHeight by beamWidth centered at "
00150                 "beamX,beamY,beamZ. \n\n"
00151                 "Ellipse beams are randomly generated on "
00152                 "the ellipses in (X,Xp), (Y,Yp), (T,E), with meanE determined "
00153                 "from meanP and the sigmas used as half-widths; tracks\n"
00154                 "are generated on the ellipse with uniform density when "
00155                 "plotted with scales such that the ellipse is a circle.\n\n"
00156                 "ASCII beams are read from a file using the format "
00157                 "specified; the formats supported are: BLTrackFile\n\n"
00158                 "Root beams are read from a .root file using the "
00159                 "TNtuple named directory/name in the file. It must have the "
00160                 "same fields as used in BLTrackFile format. Note that EventIDs "
00161                 "greater than 16,777,216 will be trncated.\n\n"
00162                 "When reading a file (ascii or root), "
00163                 "beamX and beamY are added to input tracks; if beamZ is set "
00164                 "it will overwrite the z of the track, but if it is not set "
00165                 "the z of the track in the file is kept.\n\n"
00166                 "All coordinates are centerline coordinates.\n\n"
00167                 "Multiple beam commands can be given, and they will generate\n"
00168                 "events in the order they appear in input.file.\n\n"
00169                 "Events are generated starting at firstEvent, until either\n"
00170                 "nEvents have been generated or lastEvent would be exceeded.\n\n"
00171                 "For gaussian, rectangular, and ellipse beams, the beam "
00172                 "particle can be given as either a particle name or "
00173                 "its integer PDGid. Some common beam particle names are: "
00174                 "proton, anti_proton, pi+, pi-, mu+, mu-, e+, e-, kaon+, "
00175                 "kaon-, kaon0, nu_e, anti_nu_e.\n"
00176                 "See the User's Guide for a complete list of particle names.");
00177         // initial default values:
00178         type = NONE;
00179         particle = "mu+";
00180         eventsGenerated = 0;
00181         nEvents = ALL_EVENTS;
00182         eventID = ALL_EVENTS;
00183         firstEvent = -1;
00184         lastEvent = ALL_EVENTS;
00185         meanMomentum = 200.0*MeV;
00186         beamX = 0.0;
00187         beamY = 0.0;
00188         beamZ = UNDEFINED;
00189         maxR = 1.0*kilometer;
00190         renumber = 0;
00191         weight = 1.0;
00192         secondaryTrackID = 1001;
00193         rotationMatrix = 0;
00194         sigmaX = 0.0;
00195         sigmaY = 0.0;
00196         sigmaZ = 0.0;
00197         sigmaXp = 0.0;
00198         sigmaYp = 0.0;
00199         sigmaP = 0.0;
00200         sigmaT = 0.0;
00201         sigmaE = UNDEFINED;
00202         meanXp = 0.0;
00203         meanYp = 0.0;
00204         meanT = 0.0;
00205         beamHeight = 0.0;
00206         beamWidth = 0.0;
00207         filename = "";
00208         directory = "";
00209         uid = 0;
00210         name = "";
00211         ntuple = 0;
00212         data = 0;
00213         format = "BLTrackFile";
00214         index = 0;
00215         particleDefinition = 0;
00216         particleGun = 0;
00217         trackFile = 0;
00218         nvar = 0;
00219         prevEventID = -9999;
00220 }

Member Function Documentation

G4double BLCMDbeam::myrand	(	G4double	mean,
		G4double	sigma
	)			`[inline, private]`

Referenced by nextBeamEvent().

00098                                                        {
00099                 if(sigma >= 0.0) return sigma*CLHEP::RandGauss::shoot() + mean;
00100                 return mean+sigma-2.0*sigma*G4UniformRand();
00101         }

virtual G4String BLCMDbeam::commandName ( ) [inline, virtual]

commandName() returns "beam".

Implements BLCommand.

00109 { return "beam"; }

int BLCMDbeam::command	(	BLArgumentVector &	argv,
		BLArgumentMap &	namedArgs
	)			`[virtual]`

command() implements the beam command.

Implements BLCommand.

References ALL_EVENTS, ASCII, beamX, beamY, beamZ, BLCMDbeam(), ELLIPSE, firstEvent, GAUSSIAN, BLCoordinates::getCurrentGlobal(), BLCoordinates::getCurrentRotation(), BLManager::getObject(), BLGroup::getWorld(), BLCommand::handleNamedArgs(), lastEvent, nEvents, position, printBeam(), BLCommand::printError(), RECTANGULAR, BLManager::registerBeam(), ROOT, rotation, rotationMatrix, BLGroup::setMinHeight(), BLGroup::setMinLength(), BLGroup::setMinWidth(), BLCommand::stringToRotationMatrix(), type, and UNDEFINED.

00223 {
00224         if(argv.size() != 1) {
00225                 printError("Invalid beam command -- need beam type");
00226                 return -1;
00227         }
00228 
00229         BLCMDbeam *b = new BLCMDbeam(*this);
00230 
00231         if(argv[0] == "gaussian" || argv[0] == "gauss") {
00232                 b->type = GAUSSIAN;
00233                 b->nEvents = 1;
00234         } else if(argv[0] == "rectangular" || argv[0] == "rect") {
00235                 b->type = RECTANGULAR;
00236                 b->nEvents = 1;
00237         } else if(argv[0] == "ascii" || argv[0] == "ASCII") {
00238                 b->type = ASCII;
00239                 b->nEvents = ALL_EVENTS;
00240         } else if(argv[0] == "ellipse" || argv[0] == "ELLIPSE") {
00241                 b->type = ELLIPSE;
00242                 b->nEvents = 1;
00243         } else if(argv[0] == "root") {
00244                 b->type = ROOT;
00245                 b->nEvents = ALL_EVENTS;
00246         } else {
00247                 printError("beam: invalid type '%s'",argv[0].c_str());
00248         }
00249 
00250         b->rotation = "";
00251         int retval = b->handleNamedArgs(namedArgs);
00252 
00253         if(b->firstEvent != -1 && b->lastEvent != ALL_EVENTS)
00254                 b->nEvents = b->lastEvent - b->firstEvent + 1;
00255 
00256         if(b->rotation != "") {
00257                 b->rotationMatrix = stringToRotationMatrix(b->rotation);
00258         } else {
00259                 b->rotationMatrix = new G4RotationMatrix();
00260         }
00261         // as usual, the order seems backward because (C R C^-1) C = C R.
00262         *b->rotationMatrix = *BLCoordinates::getCurrentRotation() * 
00263                                                         *b->rotationMatrix;
00264         G4ThreeVector local(b->beamX,b->beamY,
00265                                         (b->beamZ==UNDEFINED ? 0.0 : b->beamZ));
00266         BLCoordinates::getCurrentGlobal(local,b->position);
00267 
00268         // ensure beam is within the world
00269         BLGroup::getWorld()->setMinWidth(fabs(b->position[0])*2.0);
00270         BLGroup::getWorld()->setMinHeight(fabs(b->position[1])*2.0);
00271         BLGroup::getWorld()->setMinLength(fabs(b->position[2])*2.0);
00272 
00273         BLManager::getObject()->registerBeam(b);
00274 
00275         b->printBeam();
00276 
00277         return retval;
00278 }

void BLCMDbeam::printBeam ( ) [virtual]

printBeam() will print a description.

References ASCII, beamHeight, beamWidth, beamX, beamY, beamZ, directory, ELLIPSE, filename, firstEvent, format, GAUSSIAN, lastEvent, maxR, meanMomentum, meanT, meanXp, meanYp, name, nEvents, NONE, particle, RECTANGULAR, renumber, ROOT, sigmaE, sigmaP, sigmaT, sigmaX, sigmaXp, sigmaY, sigmaYp, sigmaZ, type, and weight.

Referenced by command().

00325 {
00326         switch(type) {
00327         case NONE:
00328                 printf("beam    NONE    Invalid beam command\n");
00329                 break;
00330         case GAUSSIAN:
00331                 printf("beam    GAUSSIAN particle=%s ",particle.c_str());
00332                 printf("nEvents=%d ",nEvents);
00333                 printf("firstEvent=%d ",firstEvent);
00334                 printf("lastEvent=%d ",lastEvent);
00335                 printf("beamX=%.1f ",beamX);
00336                 printf("beamY=%.1f ",beamY);
00337                 printf("beamZ=%.1f ",beamZ);
00338                 printf("maxR=%.1f ",maxR);
00339                 printf("\n\t\t");
00340                 printf("meanMomentum=%.1f ",meanMomentum);
00341                 printf("weight=%1f ",weight);
00342                 printf("\n\t\t");
00343                 printf("sigmaX=%.1f ",sigmaX);
00344                 printf("sigmaY=%.1f ",sigmaY);
00345                 printf("sigmaZ=%.1f ",sigmaZ);
00346                 printf("sigmaXp=%.5f ",sigmaXp);
00347                 printf("sigmaYp=%.5f ",sigmaYp);
00348                 printf("\n\t\t");
00349                 printf("sigmaP=%.1f ",sigmaP);
00350                 printf("sigmaT=%.3f ",sigmaT);
00351                 printf("meanXp=%.5f ",meanXp);
00352                 printf("meanYp=%.5f ",meanYp);
00353                 printf("meanT=%.3f ",meanT);
00354                 printf("\n");
00355                 break;
00356         case RECTANGULAR:
00357                 printf("beam    RECTANGULAR particle=%s ",particle.c_str());
00358                 printf("nEvents=%d ",nEvents);
00359                 printf("firstEvent=%d ",firstEvent);
00360                 printf("lastEvent=%d ",lastEvent);
00361                 printf("beamX=%.1f ",beamX);
00362                 printf("beamY=%.1f ",beamY);
00363                 printf("beamZ=%.1f ",beamZ);
00364                 printf("maxR=%.1f ",maxR);
00365                 printf("\n\t\t");
00366                 printf("meanMomentum=%.1f ",meanMomentum);
00367                 printf("weight=%1f ",weight);
00368                 printf("\n\t\t");
00369                 printf("beamHeight=%.1f ",beamHeight);
00370                 printf("beamWidth=%.1f ",beamWidth);
00371                 printf("sigmaXp=%.5f ",sigmaXp);
00372                 printf("sigmaYp=%.5f ",sigmaYp);
00373                 printf("\n\t\t");
00374                 printf("sigmaP=%.1f ",sigmaP);
00375                 printf("sigmaT=%.3f ",sigmaT);
00376                 printf("meanXp=%.5f ",meanXp);
00377                 printf("meanYp=%.5f ",meanYp);
00378                 printf("meanT=%.3f ",meanT);
00379                 printf("\n");
00380                 break;
00381         case ROOT:
00382                 printf("beam    ROOT   ");
00383                 printf("nEvents=%d ",nEvents);
00384                 printf("firstEvent=%d ",firstEvent);
00385                 printf("lastEvent=%d ",lastEvent);
00386                 printf("beamZ=%.1f ",beamZ);
00387                 printf("maxR=%.1f ",maxR);
00388                 printf("renumber=%d ",renumber);
00389                 printf("weight=%1f ",weight);
00390                 printf("\n\t\t");
00391                 printf("filename=%s ",filename.c_str());
00392                 printf("directory=%s ",directory.c_str());
00393                 printf("name=%s ",name.c_str());
00394                 printf("\n\t\t");
00395                 printf("\n");
00396                 break;
00397         case ASCII:
00398                 printf("beam    ASCII   ");
00399                 printf("nEvents=%d ",nEvents);
00400                 printf("firstEvent=%d ",firstEvent);
00401                 printf("lastEvent=%d ",lastEvent);
00402                 printf("beamZ=%.1f ",beamZ);
00403                 printf("maxR=%.1f ",maxR);
00404                 printf("renumber=%d ",renumber);
00405                 printf("weight=%1f ",weight);
00406                 printf("\n\t\t");
00407                 printf("filename=%s ",filename.c_str());
00408                 printf("format=%s ",format.c_str());
00409                 printf("\n\t\t");
00410                 printf("\n");
00411                 break;
00412         case ELLIPSE:
00413                 printf("beam    ELLIPSE particle=%s ",particle.c_str());
00414                 printf("nEvents=%d ",nEvents);
00415                 printf("firstEvent=%d ",firstEvent);
00416                 printf("lastEvent=%d ",lastEvent);
00417                 printf("beamX=%.1f ",beamX);
00418                 printf("beamY=%.1f ",beamY);
00419                 printf("beamZ=%.1f ",beamZ);
00420                 printf("maxR=%.1f ",maxR);
00421                 printf("\n\t\t");
00422                 printf("meanMomentum=%.1f ",meanMomentum);
00423                 printf("weight=%1f ",weight);
00424                 printf("\n\t\t");
00425                 printf("sigmaX=%.1f ",sigmaX);
00426                 printf("sigmaY=%.1f ",sigmaY);
00427                 printf("sigmaXp=%.5f ",sigmaXp);
00428                 printf("sigmaYp=%.5f ",sigmaYp);
00429                 printf("\n\t\t");
00430                 printf("sigmaE=%.1f ",sigmaE);
00431                 printf("sigmaT=%.3f ",sigmaT);
00432                 printf("meanXp=%.5f ",0.0);
00433                 printf("meanYp=%.5f ",0.0);
00434                 printf("meanT=%.3f ",0.0);
00435                 printf("\n");
00436                 break;
00437         }
00438 }

void BLCMDbeam::defineNamedArgs ( ) [virtual]

defineNamedArgs() defines the named arguments for this command.

Reimplemented from BLCommand.

References BLCommand::argDouble(), BLCommand::argInt(), BLCommand::argString(), beamHeight, beamWidth, beamX, beamY, beamZ, directory, filename, firstEvent, format, lastEvent, maxR, meanMomentum, meanT, meanXp, meanYp, name, nEvents, particle, renumber, rotation, secondaryTrackID, sigmaE, sigmaP, sigmaT, sigmaX, sigmaXp, sigmaY, sigmaYp, sigmaZ, uid, and weight.

00281 {
00282         argString(particle,"particle","Beam particle name");
00283         argInt(nEvents,"nEvents","Number of events to process (default=1 "
00284                 "for generating events, default=ALL for reading files) "
00285                 "set to lastEvent-firstEvent+1 if both are set.");
00286         argInt(firstEvent,"firstEvent","First event # to process (default "
00287                 "is the next sequential eventID, 1 if none)");
00288         argInt(lastEvent,"lastEvent","Last  (highest) event # to process");
00289         argDouble(beamX,"beamX","Beam location in X (mm)");
00290         argDouble(beamY,"beamY","Beam location in Y (mm)");
00291         argDouble(beamZ,"beamZ","Beam location in Z (mm)");
00292         argDouble(maxR,"maxR","Beam maximum radius (mm)");
00293         argString(rotation,"rotation","Rotation of the beam");
00294         argInt(renumber,"renumber","If nonzero, renumber events sequentially."); 
00295         argDouble(weight,"weight","Weight for events, overwritten by value from input file (1.0)."); 
00296         argInt(secondaryTrackID,"secondaryTrackID","The next TrackID for secondaries (1001)."); 
00297         argDouble(meanMomentum,"meanMomentum",
00298                                 "Gaussian Beam mean momentum (MeV/c)");
00299         argDouble(meanMomentum,"meanP","Synonym for meanMomentum.");
00300         argDouble(sigmaX,"sigmaX","Gaussian Beam sigma in X (mm)");
00301         argDouble(sigmaY,"sigmaY","Gaussian Beam sigma in Y (mm)");
00302         argDouble(sigmaZ,"sigmaZ","Gaussian Beam sigma in Z (mm)");
00303         argDouble(sigmaXp,"sigmaXp","Gaussian Beam sigma in dxdz (slope)");
00304         argDouble(sigmaYp,"sigmaYp","Gaussian Beam sigma in dydz (slope)");
00305         argDouble(sigmaP,"sigmaP","Gaussian Beam sigma in P (MeV/c)");
00306         argDouble(sigmaT,"sigmaT","Gaussian Beam sigma in T (ns)");
00307         argDouble(sigmaE,"sigmaE","Elliptical Beam sigma in E (MeV)");
00308         argDouble(meanXp,"meanXp","Gaussian Beam mean in Xp (slope)");
00309         argDouble(meanYp,"meanYp","Gaussian Beam mean in Yp (slope)");
00310         argDouble(meanT,"meanT","Gaussian Beam mean in T (ns)");
00311         argDouble(beamHeight,"beamHeight","Rectangular Beam height (mm)");
00312         argDouble(beamWidth,"beamWidth","Rectangular Beam width (mm)");
00313         argString(filename,"filename","input file name");
00314         argString(filename,"file","synonym for filename.");
00315         argString(directory,"directory","Root-file directory of NTuple");
00316         argString(directory,"category","Deprecated synonym for directory.");
00317         argInt(uid,"uid","HistoScope uid of NTuple");
00318         argString(name,"name","Root name of NTuple.");
00319         argString(format,"format","ASCII file format (Default=BLTrackFile)");
00320         argDouble(meanXp,"beamXp","Synonym for meanXp.");
00321         argDouble(meanYp,"beamYp","Synonym for meanYp.");
00322 }

virtual int BLCMDbeam::getNEvents ( ) const [inline, virtual]

getNEvents() returns the # events to process.

Reimplemented from BLBeam.

References nEvents.

00121 { return nEvents; }

void BLCMDbeam::init ( ) [virtual]

init() will initialize internal variables.

Implements BLBeam.

References ASCII, data, directory, ELLIPSE, eventsGenerated, filename, format, GAUSSIAN, BLNTuple::getNData(), index, name, NONE, ntuple, nvar, particle, particleDefinition, particleGun, BLCommand::printError(), BLNTuple::read(), RECTANGULAR, ROOT, sigmaE, sigmaP, TrackFields, trackFile, type, and UNDEFINED.

00441 {
00442         eventsGenerated = 0;
00443         index = 0;
00444         if(particleDefinition != 0) return;
00445 
00446         if(isdigit(particle(0)) || particle(0) == '-')
00447                 particleDefinition = G4ParticleTable::GetParticleTable()->
00448                                         FindParticle(atoi(particle.c_str()));
00449         else
00450                 particleDefinition = G4ParticleTable::GetParticleTable()->
00451                                         FindParticle(particle);
00452         if(!particleDefinition)
00453                 G4Exception("beam command","UnknownParticle",FatalException,
00454                                                 "Unknown particle type");
00455         particleGun = new G4ParticleGun(1);
00456         particleGun->SetParticleDefinition(particleDefinition);
00457 
00458         switch(type) {
00459         case NONE:
00460                 break;
00461         case GAUSSIAN:
00462         case RECTANGULAR:
00463                 break;
00464         case ELLIPSE:
00465                 if(sigmaE == UNDEFINED) sigmaE = sigmaP;
00466                 break;
00467         case ROOT:
00468                 { ntuple = BLNTuple::read("Root",directory,name,TrackFields,
00469                                                                 filename);
00470                   nvar = ntuple->getNData();
00471                   data = new double[nvar];
00472                 }
00473                 break;
00474         case ASCII:
00475                 if(format == "BLTrackFile") {
00476                         trackFile = new BLTrackFile(filename,"","r");
00477                 } else {
00478                         printError("beam: Unknown format '%s', known formats: "
00479                                         "BLTrackFile.",format.c_str());
00480                 }
00481                 break;
00482         }
00483 }

bool BLCMDbeam::generateReferenceParticle ( G4Event * event ) [virtual]

generateReferenceParticle() generates the reference particle.

Implements BLBeam.

00486 {
00487         return false;
00488 }

bool BLCMDbeam::nextBeamEvent ( G4Event * event ) [virtual]

nextBeamEvent() generates the next beam event.

Implements BLBeam.

References ASCII, beamHeight, beamWidth, beamZ, BLAssert, BLTF_ERROR, BLTF_OK, BLManager::clearTrackIDMap(), data, E, ELLIPSE, eventID, eventsGenerated, firstEvent, format, GAUSSIAN, BLManager::getEventID(), BLManager::getObject(), lastEvent, M_PI, maxR, meanMomentum, meanT, meanXp, meanYp, myrand(), nEvents, NONE, ntuple, nvar, particleDefinition, particleGun, position, prevEventID, BLTrackFile::read(), BLNTuple::readRow(), RECTANGULAR, renumber, ROOT, rotationMatrix, secondaryTrackID, BLManager::setNextSecondaryTrackID(), BLManager::setPrimaryTrackID(), BLBeam::setRandomSeedToGenerate(), BLBeam::setRandomSeedToTrack(), sigmaE, sigmaP, sigmaT, sigmaX, sigmaXp, sigmaY, sigmaYp, sigmaZ, trackFile, type, UNDEFINED, and weight.

00491 {
00492         G4double mass = particleDefinition->GetPDGMass();
00493         G4ThreeVector pos;
00494         G4ThreeVector direction;
00495         G4double time = 0.0;
00496         G4double momentum = 0.0;
00497         G4double ke = 0.0;
00498         G4int PDGid, trackID, parentID;
00499 
00500         BLManager *manager = BLManager::getObject();
00501 
00502         // default eventID -- changed when reading a file unless renumber!=0
00503         eventID = manager->getEventID();
00504         if(eventsGenerated == 0 && firstEvent != -1)
00505                 eventID = firstEvent;
00506 
00507         if(++eventsGenerated > nEvents) return false;
00508 
00509         // loop - must not generate any event with r > maxR
00510         // beam pos = (0,0,0), beam nominal direction = (0,0,1)
00511         // rotation and offset come later
00512         bool again=false;
00513         trackID = -1;
00514         parentID = 0;
00515         do {
00516             again = false;
00517             switch(type) {
00518             case NONE:
00519                 return false;
00520             case GAUSSIAN:
00521                 setRandomSeedToGenerate(eventID);
00522                 // meanX=beamX, meanY=beamY (implicitly)
00523                 pos[0] = myrand(0.0,sigmaX);
00524                 pos[1] = myrand(0.0,sigmaY);
00525                 pos[2] = myrand(0.0,sigmaZ);
00526                 direction[0] = myrand(meanXp,sigmaXp);
00527                 direction[1] = myrand(meanYp,sigmaYp);
00528                 direction[2] = 1.0/sqrt(1.0 + direction[0]*direction[0] +
00529                                                 direction[1]*direction[1]);
00530                 direction[0] *= direction[2];
00531                 direction[1] *= direction[2];
00532                 momentum = myrand(meanMomentum,sigmaP);
00533                 time = myrand(meanT,sigmaT);
00534                 trackID = 1;
00535                 break;
00536             case RECTANGULAR:
00537                 setRandomSeedToGenerate(eventID);
00538                 pos[0] = beamWidth*G4UniformRand()-beamWidth/2.0;
00539                 pos[1] = beamHeight*G4UniformRand()-beamHeight/2.0;
00540                 pos[2] = myrand(0.0,sigmaZ);
00541                 direction[0] = myrand(meanXp,sigmaXp);
00542                 direction[1] = myrand(meanYp,sigmaYp);
00543                 direction[2] = 1.0/sqrt(1.0 + direction[0]*direction[0] +
00544                                                 direction[1]*direction[1]);
00545                 direction[0] *= direction[2];
00546                 direction[1] *= direction[2];
00547                 momentum = myrand(meanMomentum,sigmaP);
00548                 time = myrand(meanT,sigmaT);
00549                 trackID = 1;
00550                 break;
00551             case ROOT:
00552                 { do {
00553                         if(!ntuple->readRow(data,nvar)) return false;
00554                         // data[] = x, y, z, Px, Py, Pz, t, PDGid, EventID, 
00555                         // TrackID, ParentID, Weight
00556                         if(renumber == 0) eventID = (G4int)data[8];
00557                   } while(renumber == 0 && eventID < firstEvent);
00558                   particleDefinition = G4ParticleTable::GetParticleTable()->
00559                                                 FindParticle((G4int)data[7]);
00560                   if(!particleDefinition) {
00561                         char tmp[64];
00562                         sprintf(tmp,"Invalid PDGid=%d -- Track abandoned",
00563                                                                 (int)data[7]);
00564                         G4Exception("beam command","Invalid PDGid",
00565                                                         JustWarning, tmp);
00566                         again = true;
00567                         continue;
00568                   }
00569                   mass = particleDefinition->GetPDGMass();
00570                   particleGun->SetParticleDefinition(particleDefinition);
00571                   pos[0] = data[0]*mm;
00572                   pos[1] = data[1]*mm;
00573                   pos[2] = data[2]*mm;
00574                   data[3] *= MeV;
00575                   data[4] *= MeV;
00576                   data[5] *= MeV;
00577                   momentum = sqrt(data[3]*data[3]+data[4]*data[4]+data[5]*data[5]);
00578                   direction[0] = data[3]/momentum;
00579                   direction[1] = data[4]/momentum;
00580                   direction[2] = data[5]/momentum;
00581                   time = data[6]*ns;
00582                   trackID = (G4int)data[9];
00583                   parentID = (G4int)data[10];
00584                   if(nvar >= 12) weight = data[11];
00585                 }
00586                 if(beamZ != UNDEFINED) pos[2] = 0.0;
00587                 break;
00588             case ASCII:
00589                 if(format == "BLTrackFile") {
00590                         do {
00591                                 int tmpID=0;
00592                                 BLTrackFileStatus stat=trackFile->read(pos,time,
00593                                         direction,PDGid,tmpID,trackID,parentID,
00594                                         weight);
00595                                 if(stat == BLTF_ERROR)
00596                                         G4Exception("beam command",
00597                                                 "Invalid BLTrackFile input",
00598                                                 JustWarning, "File abandoned");
00599                                 if(stat != BLTF_OK)
00600                                         return false;
00601                                 if(renumber == 0) eventID = tmpID;
00602                         } while(renumber == 0 && eventID < firstEvent);
00603                         particleDefinition = G4ParticleTable::GetParticleTable()
00604                                                         ->FindParticle(PDGid);
00605                         if(!particleDefinition) {
00606                             char tmp[64];
00607                             sprintf(tmp,"Invalid PDGid=%d -- Track abandoned",
00608                                                                         PDGid);
00609                             G4Exception("beam command","Invalid PDGid",
00610                                                         JustWarning, tmp);
00611                             again = true;
00612                             continue;
00613                         }
00614                         mass = particleDefinition->GetPDGMass();
00615                         particleGun->SetParticleDefinition(particleDefinition);
00616                         momentum = sqrt(direction[0]*direction[0]+
00617                                         direction[1]*direction[1]+
00618                                         direction[2]*direction[2]);
00619                         direction[0] /= momentum;
00620                         direction[1] /= momentum;
00621                         direction[2] /= momentum;
00622                 } else {
00623                         G4Exception("beam command","Invalid Format",
00624                                         FatalException, "");
00625                 }
00626                 if(beamZ != UNDEFINED) pos[2] = 0.0;
00627                 break;
00628             case ELLIPSE:
00629                 setRandomSeedToGenerate(eventID);
00630                 // meanX=beamX, meanY=beamY (implicitly)
00631                 double a = 2.0*M_PI*G4UniformRand();
00632                 double b = 2.0*M_PI*G4UniformRand();
00633                 double c = 2.0*M_PI*G4UniformRand();
00634                 pos[0] = sigmaX*cos(a);
00635                 pos[1] = sigmaY*cos(b);
00636                 pos[2] = myrand(0.0,sigmaZ);
00637                 direction[0] = sigmaXp*sin(a);
00638                 direction[1] = sigmaYp*sin(b);
00639                 direction[2] = 1.0/sqrt(1.0 + direction[0]*direction[0] +
00640                                                 direction[1]*direction[1]);
00641                 direction[0] *= direction[2];
00642                 direction[1] *= direction[2];
00643                 double E = sqrt(meanMomentum*meanMomentum + mass*mass) +
00644                                                                 sigmaE*sin(c);
00645                 momentum = sqrt(E*E - mass*mass);
00646                 time = sigmaT*cos(c);
00647                 trackID = 1;
00648                 break;
00649             }
00650         } while(again || sqrt(pos[0]*pos[0]+pos[1]*pos[1]) > maxR);
00651 
00652         BLAssert(trackID >= 0);
00653 
00654         if(eventID > lastEvent) return false;
00655 
00656         // apply rotation and offset
00657         if(rotationMatrix) {
00658                 direction = *rotationMatrix * direction;
00659                 pos = *rotationMatrix * pos;
00660         }
00661         pos += position;
00662 
00663         ke = sqrt(momentum*momentum + mass*mass) - mass;
00664         particleGun->SetParticleTime(time);
00665         particleGun->SetParticlePosition(pos);
00666         particleGun->SetParticleEnergy(ke);
00667         particleGun->SetParticleMomentumDirection(direction);
00668         particleGun->GeneratePrimaryVertex(event);
00669         event->SetEventID(eventID);
00670         event->GetPrimaryVertex()->SetWeight(weight);
00671         if(eventID != prevEventID) {
00672                 setRandomSeedToTrack(eventID);
00673                 manager->clearTrackIDMap();
00674                 manager->setNextSecondaryTrackID(secondaryTrackID);
00675                 prevEventID = eventID;
00676         }
00677         manager->setPrimaryTrackID(trackID,parentID);
00678         if(trackID >= secondaryTrackID) {
00679                 G4Exception("beam command","Large Primary TrackID",JustWarning,
00680                         "Confusion with secondary tracks is likely");
00681         }
00682         return true;
00683 }