How to create multiple sources using ParticleGun

Particles, Track, Event, Run and Biasing
1

Hi everyone,
I want to use ParticleGun to create multiple sources with the following parameters:

  • Shape: Cylinder, r=1.35 mm, halfz=14.5 mm
  • Particle: gamma
  • Energy: 1.173 and 1.332 MeV (cobalt 60)
  • Number of sources: 30. Each source will be placed on a sphere. They focus on the center of the sphere. The position and rotation of each source will depend on spherical coordinates, including the phi and theta angles.

Here is my code to implement a single source. Can anyone help me please to create multiple sources using the ParticleGun? Thank you in advance.

PrimaryGeneratorAction::PrimaryGeneratorAction()
{
    fParticleGun  = new G4ParticleGun(1);

    // Particle
    //
    G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
    G4String particleName;
    G4ParticleDefinition* particle = particleTable->FindParticle(particleName="gamma");
    fParticleGun->SetParticleDefinition(particle);

    InitFunction();
}

PrimaryGeneratorAction::~PrimaryGeneratorAction()
{
  	delete fParticleGun;
}

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
    // Shape and Position
    //
    G4double x0 = 0.*cm, y0 = 0.*cm, z0 = 384.2*mm;
    G4double radius = 1.35*mm;
    G4double halfz = 14.5*mm;
    G4ThreeVector pos;
    G4double rho, alpha;

    rho = G4UniformRand()*radius;
    alpha = G4UniformRand()*twopi;
	pos.setX(rho*std::sin(alpha));
	pos.setY(rho*std::cos(alpha));

    G4double z = z0 + halfz*(2*G4UniformRand() - 1);  //z uniform in (-halfz, +halfz)
    pos.setZ(z);

    fParticleGun->SetParticlePosition(pos);

    // Direction uniform in solid angle
    //
    G4double alphaMin =  0*deg;               //alpha in [0,pi]
    G4double alphaMax =  3*deg;
    G4double fCosAlphaMin, fCosAlphaMax;      //solid angle
	G4double fPsiMin, fPsiMax;
    fCosAlphaMin = std::cos(alphaMin);
    fCosAlphaMax = std::cos(alphaMax);
    
    fPsiMin = 0*deg;       //psi in [0, 2*pi]
    fPsiMax = 360*deg;

    G4double cosAlpha = fCosAlphaMin-G4UniformRand()*(fCosAlphaMin-fCosAlphaMax);
    G4double sinAlpha = std::sqrt(1. - cosAlpha*cosAlpha);
    G4double psi = fPsiMin + G4UniformRand()*(fPsiMax - fPsiMin);

    G4double ux = sinAlpha*std::cos(psi),
             uy = sinAlpha*std::sin(psi),
             uz = cosAlpha;

    G4double theta = 180*deg, phi = 0*deg;
    G4ThreeVector rotationAxis = G4ThreeVector(std::sin(theta)*std::cos(phi),
                         std::sin(theta)*std::sin(phi),
                         std::cos(theta));
    G4ThreeVector dir(ux, uy, uz);
    dir.rotateUz(rotationAxis);

    fParticleGun->SetParticleMomentumDirection(dir);

    // Energy
    //
    G4double energy = InverseCumul();  
    fParticleGun->SetParticleEnergy(energy);    


    fParticleGun->GeneratePrimaryVertex(anEvent);
}

G4double PrimaryGeneratorAction::InverseCumul()
{
    // tabulated function
    // Y is assumed positive, linear per segment, continuous 
    // --> cumulative function is second order polynomial
    
    //choose y randomly
    G4double Yrndm = G4UniformRand()*fYC[fNPoints-1];
    //find bin
    G4int j = fNPoints-2;
    while ((fYC[j] > Yrndm) && (j > 0)) j--;
    //y_rndm --> x_rndm :  fYC(x) is second order polynomial
    G4double Xrndm = fX[j];
    G4double a = fSlp[j];
    if (a != 0.) 
    {
        G4double b = fY[j]/a, c = 2*(Yrndm - fYC[j])/a;
        G4double delta = b*b + c;
        G4int sign = 1; if (a < 0.) sign = -1;
        Xrndm += sign*std::sqrt(delta) - b;    
    } 
    else if (fY[j] > 0.)
    {
        Xrndm += (Yrndm - fYC[j])/fY[j];
    };
    return Xrndm;
}

void PrimaryGeneratorAction::InitFunction()
{
    // tabulated function 
    // Y is assumed positive, linear per segment, continuous
    //
    fNPoints = 2;
    const G4double xx[] = {1.173*MeV, 1.332*MeV}; 
        
    const G4double yy[] = {0.5, 0.5};
    
    //copy arrays in std::vector and compute fMax
    //
    fX.resize(fNPoints); fY.resize(fNPoints);
    fYmax = 0.;
    for (G4int j=0; j<fNPoints; j++) {
        fX[j] = xx[j]; fY[j] = yy[j];
        if (fYmax < fY[j]) fYmax = fY[j];
    };
        
    //compute slopes
    //
    fSlp.resize(fNPoints);
    for (G4int j=0; j<fNPoints-1; j++) { 
        fSlp[j] = (fY[j+1] - fY[j])/(fX[j+1] - fX[j]);
    };
    
    //compute cumulative function
    //
    fYC.resize(fNPoints);  
    fYC[0] = 0.;
    for (G4int j=1; j<fNPoints; j++) {
        fYC[j] = fYC[j-1] + 0.5*(fY[j] + fY[j-1])*(fX[j] - fX[j-1]);
    };     
}
2

Easy peasy! Each time fParticleGun->GeneratePrimaryVertex() is called, it adds a single primary particle (with a separate primary vertex) to the event.

In your GeneratePrimaries() function, I’d factor out a lot of that into a local “AddOnePrimary(anEvent)” function. Then just modify GeneratePrimaries() to have a loop, and call your AddOnePrimary(anEvent) multiple times, within the loop.

If you want to get fancy, you can skip using G4ParticleGun entirely and do the instantiation of G4PrimaryVertex and G4PrimaryParticle yourself. Especially if you want a bunch of particles all coming from the same position (so a single G4PrimaryVertex), that’s more efficient for processing.

1 Like
3

Gents!
I revive the thread as I am stuck with the same problem. What I don’t understand is how fParticleGun->GeneratePrimaryVertex(anEvent) does actually work in the context of several sources: If I use a for-loop according to your recommendations, @mkelsey , I can clearly see the different sources in my GUI but the energy spectrum collected in my scintillator shows me a kind of Gaussian instead of a typical gamma-ray spectrum. (the proper function of the scoring volume is verified by using one source).
To me, it appears as if each single event is “accumulating” the individual source gamma energies?
Thanks!
Jurij

4

Well, yes. If you generate multiple primaries in a single event, then obviously your detector response is going to be the sum of all of those primaries combined. That’s kind of the point.

If you want to record your detector’s response to a single primary, then you should only generate a single primary per event.

5

Hi! I am stuck with a problem very similar to this. I want to generate a single primary per event but multiple events in total. I have my own primary generator class where I am looping over different particle types and positions to generate multiple vertex but they are all inside one single event. I want them to have different events as I would like to make some event based analysis of the secondaries produced by them. I cannot use gps or other types of particle generator because I am drawing the number of primaries from a distribution and their position are inside a specific volume (using fSolid → Inside() logic). Is there a way I can do this?

6

Also, I would just like to clarify, I have different particle types and their position orientation is different. They also go through different Physics process depending on the particle type.