Launch Different Particles in Same Run

Particles, Track, Event, Run and Biasing
1

Hello, I was wondering if anyone could help me with this.

I am trying to simulate a tau- and tau+ decay. I successfully coded the tau- decay, but then I was told to introduce “with equal probability” a tau+ in the simulation, so I tried to do the following:

/// \file AmBePrimaryGeneratorAction.cc
/// \brief Implementation of the AmBePrimaryGeneratorAction class; adapted from
///        Geant4 example B5 B5PrimaryGeneratorAction.cc

#include "AmBePrimaryGeneratorAction.hh"

#include "G4Event.hh"
#include "G4ParticleTable.hh"
#include "G4IonTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4Geantino.hh"
#include "G4SystemOfUnits.hh"
#include "G4RandomDirection.hh"
#include "Randomize.hh"

#include <iostream>
#include <cstdlib>
#include <ctime>

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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

	G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
	G4ParticleDefinition *particle = particleTable->FindParticle("geantino");

	G4ThreeVector pos(0., 0., 0.);
	G4ThreeVector mom(0., 0., 1.);

	fParticleGun->SetParticlePosition(pos);
	fParticleGun->SetParticleMomentumDirection(mom);
	fParticleGun->SetParticleMomentum(0.*GeV);
	fParticleGun->SetParticleDefinition(particle);
}

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

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

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

void AmBePrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
  //create vertex
  G4ParticleDefinition *particle = fParticleGun->GetParticleDefinition();
  G4double r = ((double) rand() / (RAND_MAX));
  if(particle == G4Geantino::Geantino())
  {
	  G4double energy_tot = 45000 - 1776.8600;
	  G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
      if(r >= 0.5){
		  G4ParticleDefinition *particle = particleTable->FindParticle("tau-");
	  }
	  else{
		  G4ParticleDefinition *particle = particleTable->FindParticle("tau+");
	  }
	  fParticleGun->SetParticleEnergy(energy_tot*MeV);
	  fParticleGun->SetParticlePosition(G4ThreeVector(0.,0.,0.));
	  fParticleGun->SetParticleMomentumDirection(G4RandomDirection());
	  fParticleGun->SetParticleDefinition(particle);
   }

   fParticleGun->GeneratePrimaryVertex(anEvent);
}
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

The only thing that I changed from the first tau- simulation is that I introduced that new variable “r” and the if else condition, but I does not work.

Anyone could help me understand why it does not work and how to make it do? Thank you so much in advance, Diego.

2

G4ParticleDefinition* particle = fParticleGun->GetParticleDefinition();

At the first event, current particle is a geantino
at the next event, current particle is no more a geantino …

3

Yeah thanks Marie!

I already fixed that, the only thing that is causing trouble is that I cannot make r to be actualised for every event. As you see, my intention is that a random value named r defines which tau I am launching, but for some reason it gets stuck in the first value it computes, which is 0.84 or something, and then it fires only tau-.

I am no big expert in C++ (and neither Geant4) so this may be a silly problem, but I cannot fix it.

Thanks in advance, Diego

4

the code inside
if (particle == G4Geantino::Geantino()) { … } is executed only at first event

5

Why? And how can I make it run for all events?

Thanks!

6

The key part here is

After the first event, and thus first call to GeneratePrimaries, fParticleGun’s particle definition will be tau or antitau depending on the random selection. On the second event, and second call to GeneratePrimaries, that definition is extracted and thus will never match G4Geantino::Geantino on that or any subsequent event.

As this use case is potentially swapping the particle definition every event, I’d just drop the conditional on the particle definition.

7

Okay, thanks for the explanation.

However, I still do not know how I could do this. How I could launch tau- or tau+ randomly.

Thanks again :slight_smile:

8

Throw a uniform random number between 0 and 1 (you may look in the documentation to find out how). If the value is < 0.5, set the particle gun to tau+, otherwise, set the particle gun to tau-. You may look in $G4INCLUDE to find out the class names for those two particles (ls $G4INCLUDE/G4Tau* is a good start).

9

It is enough to suppress the lines :

if (particle == G4Geantino::Geantino())
{

}

10

Thank you very much!

Diego

11

This topic was automatically closed 7 days after the last reply. New replies are no longer allowed.