Simulating cylindrical HPGe detectors using B1 example

Particles, Track, Event, Run and Biasing
1

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:

image
image864×301 8.26 KB

Here’s my detectorconstruction.cc code:

#include "G4Tubs.hh"


B1DetectorConstruction::B1DetectorConstruction()
: G4VUserDetectorConstruction(),
  fScoringVolume(0)
{ }

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

B1DetectorConstruction::~B1DetectorConstruction()
{ }

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

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:

2

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

B1PrimaryGeneratorAction::B1PrimaryGeneratorAction()
: G4VUserPrimaryGeneratorAction(),
  fParticleGun(0), 
  fEnvelopeBox(0)
{
  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");
  fParticleGun->SetParticleDefinition(particle);
  fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0.,0.,1.*cm));
  fParticleGun->SetParticleEnergy(6.*MeV);
}

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

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

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

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.";
    G4Exception("B1PrimaryGeneratorAction::GeneratePrimaries()",
     "MyCode0002",JustWarning,msg);
  }

  //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;
  
  fParticleGun->SetParticlePosition(G4ThreeVector(x0,y0,z0));
  
  G4cout << "x0, y0, z0 = " << x0 << ",\t" << y0 << ",\t" << z0 << G4endl;

  fParticleGun->GeneratePrimaryVertex(anEvent);
  
  //fParticleSource->GeneratePrimaryVertex(anEvent);
}

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

3

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?

4

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.

5

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:

6

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.

7 
#include "G4Tubs.hh"

?

8

Yes I already did that but didn’t work

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

9

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
10

@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 PrimaryGeneratorAction.cc file. Everything works properly now and I get random positions for the photons. Thank you!!!