Using Particle Gun to Sample Particle Energy from Spectra?

Particles, Track, Event, Run and Biasing
1

Hello All,

I understand how to sample particle’s energy using a given spectra in a separate file (energy with counts). I use the following (as in here):

# energy distribution
/gps/ene/type User
/gps/hist/type energy
/control/execute spec.dat

My question is, what is the equivalent of the above when particle gun is used instead of gps?

1 Like
2

The closest equivalent would be to store the vector in your PrimaryGeneratorAction, sample from it in GeneratePrimaries, then set the gun’s energy before calling GeneratePrimaryVertex. In general you can use the gps exactly like the particle gun, so if you want this functionality there’s little reason to stick tot he gun.

Also note that there is a /gps/hist/file command, which lets you store your spectrum in a more natural format, rather than a bunch of /gps/hist/point lines.

3

Thanks @bloer, are there any relevant examples/resources where the energy is sampled from a vector using the ParticleGun?

4

I use the following:

In BeginOfRunAction, fEnergyOverride is a G4PhysicsVector, which has a handy method FindLinearEnergy to sample the distribution if it’s a CDF:

if(overrideEnergyFile != ""){
    std::ifstream fin(overrideEnergyFile);
    if(!fEnergyOverride.Retrieve(fin, /*ascii=*/true)){
      G4ExceptionDescription msg;
      msg<<"Unable to read energy histogram from "<<overrideEnergyFile;
      G4Exception("BGMPrimaryGeneratorAction::BeginOfRunAction","1",
                  FatalException, msg);
    }
    // convert to a CDF
    double sum = 0;
    for(size_t i=0; i<fEnergyOverride.GetVectorLength(); ++i){
      sum += fEnergyOverride[i];
      fEnergyOverride.PutValue(i, sum);
    }
    if(sum <= 0){
      G4ExceptionDescription msg;
      msg<<"No entries in histogram from "<<overrideEnergyFile;
      G4Exception("BGMPrimaryGeneratorAction::BeginOfRunAction","1",
                  FatalException, msg);
    }
    fEnergyOverride.ScaleVector(1, 1./sum);
  }

Then to sample in `GeneratePrimaryEvent:

fGun->GeneratePrimaryVertex(anEvent);
if(fEnergyOverride.IsFilledVectorExist()){
    G4PrimaryVertex* vertex = anEvent->GetPrimaryVertex();
    while(vertex){
      for(int i=0; i<vertex->GetNumberOfParticle(); ++i){
        G4PrimaryParticle* primary = vertex->GetPrimary(i);
        G4double energy = fEnergyOverride.FindLinearEnergy(G4UniformRand());
        primary->SetKineticEnergy(energy);
      }
      vertex = vertex->GetNext();
    }
  }

Note here I let the gun do it’s thing and then override the energy, rather than set the gun’s energy before calling GeneratePrimaryVertex. The latter would be more efficient, and mine is only done this way for bad historical reasons…

5

See also examples/extended/eventgenerator/particleGun : case 2