Hi G4 experts!
I am completely new to g4, but i need to simulate some 176Lu decays in a LYSO scintillator windowlet embedded in a complex geometry. Once the geometry is built, i need to create the isotope and let it decay, therefore I use the following code as B1PrimaryGeneratorAction.cc
:
#include "B1PrimaryGeneratorAction.hh"
#include <stdio.h>
#include <math.h>
#include "G4LogicalVolumeStore.hh"
#include "G4LogicalVolume.hh"
#include "G4Box.hh"
#include "G4RunManager.hh"
#include "G4ParticleGun.hh"
#include "G4ParticleTable.hh"
#include "G4IonTable.hh"
#include "G4ParticleDefinition.hh"
#include "G4SystemOfUnits.hh"
#include "Randomize.hh"
#include "time.h"
#include "G4GenericMessenger.hh"
#include "G4PhysicalVolumeStore.hh"
#include "G4Event.hh"
#include "G4Geantino.hh"
B1PrimaryGeneratorAction::B1PrimaryGeneratorAction()
: G4VUserPrimaryGeneratorAction(),
fParticleGun(0),
fEnvelopeBox(0)
{
G4int n_particle = 1;
fParticleGun = new G4ParticleGun(n_particle);
fParticleGun->SetParticleEnergy(0*eV);
fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.));
}
B1PrimaryGeneratorAction::~B1PrimaryGeneratorAction()
{
delete fParticleGun;
}
void B1PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
G4double xlo = -0.5*cm;
G4double ylo = -0.5*cm;
G4double zlo = -0.15004*cm;
G4double xhi = 0.5*cm;
G4double yhi = 0.5*cm;
G4double zhi = -0.00004*cm;
const G4String name="Scintillator";
G4bool verbose=true;
G4PhysicalVolumeStore * volumestore;
G4VPhysicalVolume * pvol = volumestore->GetVolume(name,verbose);
G4LogicalVolume * lvol = pvol->GetLogicalVolume();
G4VSolid * solid = lvol->GetSolid();
G4ThreeVector point;
G4int maxtries=10000, itry=1;
do {
point.set(xlo + G4UniformRand()*(xhi-xlo),
ylo + G4UniformRand()*(yhi-ylo),
zlo + G4UniformRand()*(zhi-zlo));
} while (!solid->Inside(point) && ++itry < maxtries);
if (itry == maxtries)
G4cerr << "Unable to find a point inside your volume!" << G4endl;
fParticleGun->SetParticlePosition(point);
if (fParticleGun->GetParticleDefinition() == G4Geantino::Geantino()) {
G4int Z = 71, A = 176;
G4int Level = 0;
G4double ionCharge = 0.*eplus;
G4double excitEnergy = 0.*keV;
G4ParticleDefinition* ion = G4IonTable::GetIonTable()->GetIon(Z,A,excitEnergy);
fParticleGun->SetParticleDefinition(ion);
fParticleGun->SetParticleCharge(ionCharge);
}
//create vertex
//
fParticleGun->GeneratePrimaryVertex(anEvent);
}
And a PhysicsList.cc
file in which i included the following lines:
#include "PhysicsList.hh"
#include "G4EmStandardPhysics.hh"
#include "G4ParticleTypes.hh"
#include "G4ParticleTable.hh"
#include "G4IonConstructor.hh"
#include "G4Gamma.hh"
#include "G4Electron.hh"
#include "G4Positron.hh"
#include "G4PhysicsListHelper.hh"
#include "G4RadioactiveDecay.hh"
#include "G4UnitsTable.hh"
#include "G4LossTableManager.hh"
#include "G4RadioactiveDecayPhysics.hh"
#include "G4NeutronTrackingCut.hh"
#include "G4DecayPhysics.hh"
#include "G4SystemOfUnits.hh"
PhysicsList::PhysicsList() :
G4VUserPhysicsList()
{
//add new units for radioActive decays
//
const G4double minute = 60*second;
const G4double hour = 60*minute;
const G4double day = 24*hour;
const G4double year = 365*day;
new G4UnitDefinition("minute", "min", "Time", minute);
new G4UnitDefinition("hour", "h", "Time", hour);
new G4UnitDefinition("day", "d", "Time", day);
new G4UnitDefinition("year", "y", "Time", year);
//rdecay02constructor
G4LossTableManager::Instance();
defaultCutValue =1.*mm;
SetVerboseLevel(1);
//default physics
fParticleList = new G4DecayPhysics();
//default physics
fRaddecayList = new G4RadioactiveDecayPhysics();
// EM physics
fEmPhysicsList = new G4EmStandardPhysics();
}
PhysicsList::~PhysicsList()
{
delete fRaddecayList;
delete fEmPhysicsList;
}
void PhysicsList::ConstructParticle()
{
fParticleList->ConstructParticle();
}
void PhysicsList::ConstructProcess()
{
AddTransportation();
// em
fEmPhysicsList->ConstructProcess();
// decays
fParticleList->ConstructProcess();
fRaddecayList->ConstructProcess();
//radioactive decay
G4RadioactiveDecay* radioactiveDecay = new G4RadioactiveDecay();
radioactiveDecay->SetHLThreshold(-1.*s);
radioactiveDecay->SetICM(true);
//Internal Conversion
radioactiveDecay->SetARM(false);
//Atomic Rearangement
G4PhysicsListHelper* ph = G4PhysicsListHelper::GetPhysicsListHelper();
ph->RegisterProcess(radioactiveDecay, G4GenericIon::GenericIon());
//Deexcitation (in case of Atomic Rearangement)
G4UAtomicDeexcitation* de = new G4UAtomicDeexcitation();
de->SetFluo(true);
de->SetAuger(true);
de->SetPIXE(false);
G4LossTableManager::Instance()->SetAtomDeexcitation(de);
G4cout << "### PhysicsList::ConstructProcess is done" << G4endl;
}
void PhysicsList::SetCuts()
{
SetCutsWithDefault();
}
The problem is that if I run the simulation it seems that no radioactive decay is produced… I am probably wrongly writing/importing the correct physics to be used, but i really don’t know what i should modify to get it right…
Can you please help me?
Thanks in advance,
Cecilia