Simulating cylindrical HPGe detectors using B1 example

Hello experts,

I am running a simulation for a HPGe detector with a hole in the middle and I have developed this from the example B1. The issue I am having is I cannot run the simulation without randomizing the position of x0, y0, z0 and I will like to do this in such a way that the particle source fires in a circular manner along the symmetry axis of the detector because it is a cylinder. I have tried using GPS but I don’t think I implemented properly because all my initial position is set to zero.

I get this output messsage:

Here’s my code:

#include "G4Tubs.hh"

: G4VUserDetectorConstruction(),
{ }


{ }


G4VPhysicalVolume* B1DetectorConstruction::Construct()
  // Get nist material manager
  G4NistManager* nist = G4NistManager::Instance();
  // Envelope parameters
  G4double env_sizeXY = 10*cm, env_sizeZ = 10*cm;
  G4Material* env_mat = nist->FindOrBuildMaterial("G4_Ge");

  G4double PI = 3.141592653589793238;
  G4double pRMin = 0.2*cm, pRMax = 2.5*cm, pDz = 2.1*cm, pSPhi = 0., pDPhi = 2*PI;
  // Option to switch on/off checking of volumes overlaps
  G4bool checkOverlaps = true;

  // World
  G4double world_sizeXY = 1.2*env_sizeXY;
  G4double world_sizeZ  = 1.2*env_sizeZ;
  G4Material* world_mat = nist->FindOrBuildMaterial("G4_AIR");
  G4Box* solidWorld =    
    new G4Box("World",                       //its name
       0.5*world_sizeXY, 0.5*world_sizeXY, 0.5*world_sizeZ);     //its size
  G4LogicalVolume* logicWorld =                         
    new G4LogicalVolume(solidWorld,          //its solid
                        world_mat,           //its material
                        "World");            //its name
  G4VPhysicalVolume* physWorld = 
    new G4PVPlacement(0,                     //no rotation
                      G4ThreeVector(),       //at (0,0,0)
                      logicWorld,            //its logical volume
                      "World",               //its name
                      0,                     //its mother  volume
                      false,                 //no boolean operation
                      0,                     //copy number
                      checkOverlaps);        //overlaps checking
  // Envelope

  G4Tubs* solidEnv =    
    new G4Tubs("Envelope",                    //its name
        pRMin, pRMax, pDz, pSPhi, pDPhi);    //its size*/
    G4LogicalVolume* logicEnv =                         
    new G4LogicalVolume(solidEnv,            //its solid
                        env_mat,             //its material
                        "Envelope");         //its name
  new G4PVPlacement(0,                       //no rotation
                    G4ThreeVector(),         //at (0,0,0)
                    logicEnv,                //its logical volume
                    "Envelope",              //its name
                    logicWorld,              //its mother  volume
                    false,                   //no boolean operation
                    0,                       //copy number
                    checkOverlaps);          //overlaps checking

  // Set Envelope as scoring volume
  fScoringVolume = logicEnv;

  //always return the physical World
  return physWorld;

And here is my PrimaryGeneratorAction file:

Perhaps I also have to modify the EnvelopeBox params to EnvelopeTubs in Tried but didn’t succeed. Any help will be appreciated!

: G4VUserPrimaryGeneratorAction(),
  G4int n_particle = 1;
  fParticleGun  = new G4ParticleGun(n_particle);
  //fParticleGun  = new G4GeneralParticleSource();

  // default particle kinematic
  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
  G4String particleName;
  G4ParticleDefinition* particle
    = particleTable->FindParticle(particleName="gamma");


  delete fParticleGun;


void B1PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
  //this function is called at the begining of each event

  // In order to avoid dependence of PrimaryGeneratorAction
  // on DetectorConstruction class we get Envelope volume
  // from G4LogicalVolumeStore.
  G4double envSizeXY = 0;
  G4double envSizeZ = 0;

  if (!fEnvelopeBox)
    G4LogicalVolume* envLV
      = G4LogicalVolumeStore::GetInstance()->GetVolume("Envelope");
    //if ( envLV ) fEnvelopeBox = dynamic_cast<G4Box*>(envLV->GetSolid());
    if ( envLV ) fEnvelopeBox = dynamic_cast<G4Box*>(envLV->GetSolid());

  if ( fEnvelopeBox ) {
    envSizeXY = fEnvelopeBox->GetXHalfLength()*2;
    envSizeZ = fEnvelopeBox->GetZHalfLength()*2;
  else  {
    G4ExceptionDescription msg;
    msg << "Envelope volume of box shape not found.\n"; 
    msg << "Perhaps you have changed geometry.\n";
    msg << "The gun will be place at the center.";

  //G4double size = 1.0; 
  G4double x0 = envSizeXY * (G4UniformRand()-0.5);
  G4double y0 = envSizeXY * (G4UniformRand()-0.5);
  G4double z0 = -5.5*cm;
  //G4double z0 = fParticleGun->GetParticlePosition();
  G4double x0, y0, z0;
  G4cout << "x0, y0, z0 = " << x0 << ",\t" << y0 << ",\t" << z0 << G4endl;


Sorry to bother you again @weller and @mkelsey but please kindly take a look

get rid of this line:

  G4double x0, y0, z0;

you re-declare these variables, after setting values to them, hence they are 0 (or undefined?) Does this not throw compiler warnings with clear text notification?

Thank you for your quick response @weller

I already got rid of it. I only used it to check what happens if I hard code the coordinates but forgot to make the edit here before posting.

you certainly have to adapt the code to your changes. For example,

if ( envLV ) fEnvelopeBox = dynamic_cast<G4Box*>(envLV->GetSolid());

is a G4Box, but you have a G4Tubs. Accordingly, fEnvelopeBox->GetXHalfLength() is not defined anymore, you will need to find the class reference, and look up the commands you need. Hint: GetZHalfLength() and GetOuterRadius() could be prime candidates :wink:

Thank you for the tips @weller but I have tried checking the reference for G4Tubs but haven’t made any progress because of fEnvelopeBox class. How can I change this to something like fEnvelopeTubs?

I tried this:

if ( envLV ) fEnvelopeBox = dynamic_cast<G4Tubs*>(envLV->GetSolid());

but G4Tubs is not recognized for dynamic_cast.

I apologize for the naive question but I am still a Geant4 beginner.

#include "G4Tubs.hh"


Yes I already did that but didn’t work

I even tried replacing EnvelopeBox with EnvelopeTubs, but ran into different errors. @weller

this is now c++ related only, and not geant4 anymore. I would recommend tutorials to get better understanding of basics.

don’t forget to also change fEnvelopeBox as defined in the B1PrimaryGeneratorAction.hh file…

1 Like

@weller I cannot thank you enough! I really appreciate the timely response; you have been so helpful!

I found out that the issue was because I didn’t define EnvelopeTubs in my PrimaryGeneratorAction.hh header file and also deleted the line: #include G4Tubs.hh in my file. Everything works properly now and I get random positions for the photons. Thank you!!!