Geant4 microyz modified to neurons geometry parameters

Geometry, Fields and Transportation
1

Please fill out the following information to help in answering your question, and also see tips for posting code snippets. If you don’t provide this information it will take more time to help with your problem!

_Geant4 Version:_11.3
_Operating System:_Sequoia 15.2
Compiler/Version:
CMake Version:

Hello!
I try to modify microyz module for a network of 10 neurons. The simulation volume is 640 x 280 x 25 um. The diameters and lengths of the cylinders, representing the dendrites, are in the ranges of 0.226–6.525 um and 0.200–9.918 um, respectively. The sphere diameters within the soma area are ran-
ged from 2.025 um to 19.35 um. The volumes of spherical and cylindrical compartments simulating different parts of the neuron are in the ranges of 34.765–3791.590 um3 and 0.011–280.911 um3, correspondingly. In total, the neuronal network is reconstructed by 96 spheres and 18,082 cylinders. This information is from Radiation damage to neuronal cells: Simulating the energy deposition and water radiolysis in a small neural network, 2016.

I try to add the parameters of the network in DetectorConstruction.cc

#include "DetectorConstruction.hh"

#include "DetectorMessenger.hh"
#include "TrackerSD.hh"
//#include "NeuronSD.hh"

#include "G4Box.hh"
#include "G4GeometryManager.hh"
#include "G4GeometryTolerance.hh"
#include "G4NistManager.hh"
#include "G4PVPlacement.hh"
#include "G4SDManager.hh"
#include "G4SystemOfUnits.hh"
#include "G4UnitsTable.hh"
#include "G4UserLimits.hh"
#include "G4Orb.hh" //add
#include "G4LogicalVolume.hh"//add
//#include "NeuronSensitiveDetector.hh" //add then removed
#include "G4Sphere.hh" //add
#include "G4Tubs.hh" //add
#include "G4ThreeVector.hh"//add
#include "G4RunManager.hh" //add

#include "DetectorConstruction.hh"

#include "DetectorMessenger.hh"
#include "TrackerSD.hh"
//#include "NeuronSD.hh"

// sarah add
DetectorConstruction::DetectorConstruction() : G4VUserDetectorConstruction(), fpWaterMaterial(nullptr),
      fSomaRadius(1 * um),
      fDendriteRadius(3 * um),
      fDendriteLength(6 * um),
      fSpacing(50),
      fTotalNeurons(10),
      fpTrackingCut(11 * eV),
      fpMaxStepSize(1.0 * um)
{
  // Create commands for interactive definition of the detector

  fDetectorMessenger = std::make_unique<DetectorMessenger>(this);
} 

//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  
DetectorConstruction::~DetectorConstruction() = default;
  
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  
G4VPhysicalVolume* DetectorConstruction::Construct()
{
  // Define materials
  DefineMaterials();                  

  // Define volumes
  return DefineVolumes();
}

  
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  
void DetectorConstruction::DefineMaterials()
{
  // Water is defined from NIST material database
      G4NistManager* nistManager = G4NistManager::Instance();
    G4Material* water = nistManager->FindOrBuildMaterial("G4_WATER");

    fpWaterMaterial = water;
  
  // Print materials
  G4cout << *(G4Material::GetMaterialTable()) << G4endl;
}
  
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  //Define Volumes
G4VPhysicalVolume* DetectorConstruction::DefineVolumes()
{
    
    // Define network size in micrometers (211.82 x 236.12 x 18.9 micrometers)
    G4double worldLengthX = 640 * um;  // World length in X ( micrometers)
    G4double worldLengthY = 280 * um;  // World length in Y ( micrometers)
    G4double worldLengthZ = 25 * um;    // World length in Z ( micrometers)
    

  //G4double worldLength = 10 * m; m ---> 10 cm
//G4double worldLength = 10 * cm;  // World size
/*
    // Define neuron dimensions
    G4double somaRadius = 10 * um;   // Soma radius
    G4double dendriteRadius = 2 * um;  // Dendrite radius
    G4double dendriteLength = 50 * um; // Dendrite length
    G4double spacing = 50 * um;  // Spacing between neurons
*/
  // World

  G4GeometryManager::GetInstance()->SetWorldMaximumExtent(worldLengthX);
  
  G4cout << "Computed tolerance = "
         << G4GeometryTolerance::GetInstance()->GetSurfaceTolerance() / um << " um" << G4endl;
  
// Define world volume (box with specified dimensions)
  auto* worldS = new G4Box("world",  // its name
                           worldLengthX / 2, worldLengthY / 2, worldLengthZ / 2);  // its size
  
  auto* worldLV = new G4LogicalVolume(worldS,  // its solid
                                      fpWaterMaterial,  // its material
                                      "World_LV");  // its name
  G4VPhysicalVolume* worldPV = new G4PVPlacement(nullptr,  // no rotation
                                                 G4ThreeVector(),  // at (0,0,0)
                                                 worldLV,  // its logical volume
                                                 "World",  // its name
                                                 nullptr,  // its mother  volume
                                                 false,  // no boolean operations
                                                 0,  // copy number
                                                 false);  // checking overlaps
  // user limits for world
  worldLV->SetUserLimits(new G4UserLimits(fpMaxStepSize, DBL_MAX, DBL_MAX, fpTrackingCut));   
  
// loop to create 10 neurons
       for (int i = 0; i < fTotalNeurons; ++i) {
        // Soma (shape sphere)
        G4ThreeVector somaPos(i * fSpacing, 0, 0);
        auto* soma = new G4Orb("Soma_" + std::to_string(i), fSomaRadius);
        auto* somaLV = new G4LogicalVolume(soma, fpWaterMaterial, "Soma_LV_" + std::to_string(i));
        new G4PVPlacement(nullptr, somaPos, somaLV, "Soma_PV_" + std::to_string(i), worldLV, false, i);

        // Dendrite (shape cylinder)
        G4ThreeVector dendritePos(i * fSpacing, 0, fDendriteLength / 2);
        auto* dendrite = new G4Tubs("Dendrite_" + std::to_string(i), 0, fDendriteRadius, fDendriteLength / 2, 0, 360 * deg);
        auto* dendriteLV = new G4LogicalVolume(dendrite, fpWaterMaterial, "Dendrite_LV_" + std::to_string(i));
        new G4PVPlacement(nullptr, dendritePos, dendriteLV, "Dendrite_PV_" + std::to_string(i), worldLV, false, i + 10);
    }

// Print parameters used

  PrintParameters();
                           
  return worldPV;
}
                                      
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
  
void DetectorConstruction::ConstructSDandField()
{
  // Sensitive detector
                                                 
   auto* sdManager = G4SDManager::GetSDMpointer();
 // Create a single sensitive detector for the entire neuron network
    auto* trackerSD = new TrackerSD("TrackerNeuronNetworkSD"); // replaced new neuronsensitivedetector with TrackerSD
    sdManager->AddNewDetector(trackerSD);

    // Loop to create and assign sensitive detectors to all soma and dendrite logical volumes
    for (int i = 0; i < fTotalNeurons; ++i) {
//for the entire network
      // Soma logical volume name
        G4String somaLVName = "Soma_LV_" + std::to_string(i);

        // Dendrite logical volume name
        G4String dendriteLVName = "Dendrite_LV_" + std::to_string(i);

        // Assign the same sensitive detector to the soma
        SetSensitiveDetector(somaLVName, trackerSD);

        // Assign the same sensitive detector to the dendrite
        SetSensitiveDetector(dendriteLVName, trackerSD);
    }

    G4cout << "Single sensitive detector 'TrackerNeuronNetworkSD' assigned to the entire neuron network." << G4endl;

    }
      

     // Tracking cut set
               
void DetectorConstruction::SetTrackingCut(const G4double& value)
{
  fpTrackingCut = value;
}
  
//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
         
void DetectorConstruction::SetMaxStepSize(const G4double& value)
{
  fpMaxStepSize = value;
} 

//set raduis of soma and dendrites
void DetectorConstruction::SetRadius(const G4double& somaRadius, const G4double& dendriteRadius)
{
  fSomaRadius = somaRadius;
fDendriteRadius = dendriteRadius;
}

//set length of dendrite
void DetectorConstruction::SetDendriteLength(const G4double& dendriteLength)
{
    fDendriteLength = dendriteLength;
    G4cout << "Dendrite length set to: " << G4BestUnit(fDendriteLength, "Length") << G4endl;
    G4RunManager::GetRunManager()->GeometryHasBeenModified();
}

// set total neurons
void DetectorConstruction::SetTotalNeurons(const G4double& totalNeurons)
{
    fTotalNeurons = static_cast<int>(totalNeurons);
    G4cout << "Total neurons set to: " << fTotalNeurons << G4endl;
    G4RunManager::GetRunManager()->GeometryHasBeenModified();
}

  

//print parameters
void DetectorConstruction::PrintParameters() const
{
  G4cout << "\n---------------------------------------------------------\n";
  G4cout << "---> The tracking cut is set to " << G4BestUnit(fpTrackingCut, "Energy") << G4endl;
  G4cout << "---> The maximum step size is set to " << G4BestUnit(fpMaxStepSize, "Length")
         << G4endl;
  G4cout << "---> The soma radius is set to " << G4BestUnit(fSomaRadius, "Length") << G4endl;
  G4cout << "---> The dendrite  radius is set to " << G4BestUnit(fDendriteRadius, "Length") << G4endl;
  G4cout << "---> The dendrite length is set to " << G4BestUnit(fDendriteLength, "Length") << G4endl;
    G4cout << "  Spacing: " << fSpacing / um << " um" << G4endl;
        G4cout << "  Total neurons: " << fTotalNeurons << G4endl;

  G4cout << "---> The world dimensions are set to: "
           << G4BestUnit(worldLengthX, "Length") << " x "
           << G4BestUnit(worldLengthY, "Length") << " x "
           << G4BestUnit(worldLengthZ, "Length") << G4endl;
  G4cout << "\n---------------------------------------------------------\n";
}

But I get this error although try to change spacing between neurons or change the width of world fron 640 to 699 um.

### ===  Deexcitation model UAtomDeexcitation is activated for 1 region:
          DefaultRegionForTheWorld  1  1  0

### === G4UAtomicDeexcitation::InitialiseForNewRun()
### ===  Auger flag: 1
### ===  Ignore cuts flag:   1
### ===  PIXE model for hadrons: Empirical
### ===  PIXE model for e+-:     Livermore

-------- EEEE ------- G4Exception-START -------- EEEE -------
*** G4Exception : GeomMgt0002
      issued by : G4SmartVoxelHeader::BuildNodes()
PANIC! - Overlapping daughter with mother volume.
         Daughter physical volume Soma_PV_1
         is entirely outside mother logical volume World_LV !!
*** Fatal Exception *** core dump ***
 **** Track information is not available at this moment
 **** Step information is not available at this moment

-------- EEEE -------- G4Exception-END --------- EEEE -------


*** G4Exception: Aborting execution ***

Could you help me, please? I changed also in src and include files for microyz but the error refers to the geometery of this network…

2

It looks that you forgotten to specify units for fSpacing.

3

Yes. I have a question. Why in the origin code, using world length of box of water to be 10 m! and the fradius of sensitive detector to collect the hits to be 5 nm. Why using a huge place? Is it normal to replace this dimension of the box with the dimension 640 x 280 x 25 um, with the same radius of target?

4

If the code of PrimaryGeneratorAction.cc

PrimaryGeneratorAction::PrimaryGeneratorAction() : G4VUserPrimaryGeneratorAction()
{
  G4int nofParticles = 1;
  fParticleGun = std::make_unique<G4ParticleGun>(nofParticles);

  G4ParticleDefinition* particleDefinition =
    G4ParticleTable::GetParticleTable()->FindParticle("e-");

  fParticleGun->SetParticleDefinition(particleDefinition);
  fParticleGun->SetParticleMomentumDirection(G4RandomDirection());
  fParticleGun->SetParticleEnergy(1 * keV);
}

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

PrimaryGeneratorAction::~PrimaryGeneratorAction() = default;

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

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
  fParticleGun->SetParticlePosition(G4ThreeVector(0., 0., 1.));

  fParticleGun->SetParticleMomentumDirection(G4RandomDirection());

  fParticleGun->GeneratePrimaryVertex(anEvent);
}

To shoot particles to XY plane of the full 640 x 280 x 25 um simulation volume of box of water . The trajectories of
particle beams were generated inside the box in the direction par-
allel to Z axis. The distribution of the particle tracks within the XY
plane was set randomly.
Should I change anything in this code, please?

5

Hi Sarah

Not sure what you mean,“The trajectories of particle beams were generated inside the box in the direction parallel to Z axis. The distribution of the particle tracks within the XY plane was set randomly.” How can they be random in the xy-plane and simultaneously parallel to the z-axis? Do you mean perpendicular to the z-axis? That’s strange. G4RandomDirection() should generate directions random in 3D.

You don’t need

fParticleGun->SetParticleMomentumDirection(G4RandomDirection());

in the constructor.

Otherwise the code looks OK and looks like it will generate tracks in random directions.

John

6

Thanks!
Happy New Year !
I think Particles travel in the Z direction with random distributions of positions in XY plane.
if the full dimension of the world is 640 x 280 x 25 um fill with just water. So, is it OK to modify the code to be as this description, please?

// sarah modify here
PrimaryGeneratorAction::PrimaryGeneratorAction()
    : G4VUserPrimaryGeneratorAction()
{
    G4int nofParticles = 1;
    fParticleGun = std::make_unique<G4ParticleGun>(nofParticles);

    // Define the particle type (proton in this case)
    G4ParticleDefinition* particleDefinition =
        G4ParticleTable::GetParticleTable()->FindParticle("proton");
    fParticleGun->SetParticleDefinition(particleDefinition);

    // Set default properties
    fParticleGun->SetParticleEnergy(1 * MeV);
}

PrimaryGeneratorAction::~PrimaryGeneratorAction() = default;

void PrimaryGeneratorAction::GeneratePrimaries(G4Event* anEvent)
{
    // Simulation box dimensions
    G4double boundingXHalfLength = 320 * um;
    G4double boundingYHalfLength = 140 * um;
    G4double zPosition = -12.5 * um; // Mid-plane of the 25-um box

    // Set random position within the box (uniform distribution)
    G4double x = CLHEP::RandFlat::shoot(-boundingXHalfLength, boundingXHalfLength);
    G4double y = CLHEP::RandFlat::shoot(-boundingYHalfLength, boundingYHalfLength);

    fParticleGun->SetParticlePosition(G4ThreeVector(x, y, zPosition));

    // Set momentum direction along +Z axis
    fParticleGun->SetParticleMomentumDirection(G4ThreeVector(0., 0., 1.));

    // Generate the particle
    fParticleGun->GeneratePrimaryVertex(anEvent);
}
7

Thanks!
Happy New Year !