Radioactive isotope in scintillator - no photon production

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

1 Like