HaarigerHarald / geant4_pybind

Alternative Python bindings for Geant4
The Unlicense
35 stars 7 forks source link

I could not start simulation #22

Closed saimegemenyucel closed 9 months ago

saimegemenyucel commented 9 months ago

Hi everyone, I am a newbie in geant4 and stuck in somewhere again. I do not know what did I wrong. When I execute my geant4_pybind python code the UI session starts but do not visualize anything. UI give me my defined materials list and suddenly crushes. Screenshot 2024-02-09 021523 and then when I close interactive mode in visual studio it says: Screenshot 2024-02-09 021530

HaarigerHarald commented 9 months ago

Could you post the simulation code?

saimegemenyucel commented 9 months ago

`from calendar import c import cmd from filecmp import cmp from os import name from re import S from geant4_pybind import *

import sys import cmath import math

G4AnalysisManager = G4RootAnalysisManager

class DetectorConstruction(G4VUserDetectorConstruction):

def __init__(self):
    super().__init__()
    self.fCheckOverLaps = True

def DefineMaterials(self):
    nist = G4NistManager.Instance()

    #vacuum
    G4Material("Galactic", z = 1, a = 1.01*g/mole, density= universe_mean_density, state= kStateGas, temp= 2.73*kelvin, pressure= 3e-18*pascal)

    #PSU Polysulfone
    isotope = False
    C = nist.FindOrBuildElement("C",isotope)
    H = nist.FindOrBuildElement("H",isotope)
    O = nist.FindOrBuildElement("O",isotope)
    S = nist.FindOrBuildElement("S",isotope)
    self.PSU = G4Material(name= "Polysulfone", density= 1.183*g/cm3, nComponents= 4)
    self.PSU.AddElement(C, 6)
    self.PSU.AddElement(H, 4)
    self.PSU.AddElement(O, 2)
    self.PSU.AddElement(S, 1)

    print(G4Material.GetMaterialTable())

def Construct(self):
    self.DefineMaterials()
    nist = G4NistManager.Instance()

    checkOverlaps = True

    world_sizeXY = 80.0*cm
    world_sizeZ = 80.0*cm

    world_mat = G4Material.GetMaterial("Galactic")
    target_mat = G4Material.GetMaterial("Polysulfone")

    if world_mat == None or target_mat == None:
        msg = "Cannot retrieve materials already defined."
        G4Exception("DetectorConstruction::DefineVolumes()",
                    "MyCode0001", FatalException, msg)

    #WORLD
    solidWorld = G4Box("World", world_sizeXY, world_sizeXY, world_sizeZ)

    logicWorld = G4LogicalVolume(solidWorld, world_mat, "World")

    physWorld = G4PVPlacement(None,
                              G4ThreeVector(0, 0, 0),
                              logicWorld,
                              "World",
                              None,
                              False,
                              0,
                              checkOverlaps)

    #TARGET
    solidTarget = G4Box("Target", 1.0*cm, 5.0*cm, 5.0*cm)

    logicTarget = G4LogicalVolume(solidTarget, target_mat, "Target")

    physTarget = G4PVPlacement(None,
                               G4ThreeVector(0, 0, 5.0*cm),
                               logicTarget,
                               "Target",
                               logicWorld,
                               False,
                               0,
                               checkOverlaps)

    #DETECTOR
    caloriMeterSolid = G4Tubs("Detector", 0, 20.0*cm, 50.0*cm, 0*deg, 360*deg)
    caloriMeterLogical = G4LogicalVolume(caloriMeterSolid, world_mat, "Detector")
    caloriMeterPhys = G4PVPlacement(None,
                                    G4ThreeVector(15.0*cm, 0, 0),
                                    None,
                                    "Detector",
                                    caloriMeterLogical,
                                    False,
                                    0,
                                    checkOverlaps)

    caloriMeterVisAtt = G4VisAttributes(G4Colour(1, 0, 0))
    caloriMeterLogical.SetVisibility(True)
    caloriMeterLogical.SetVisAttributes(caloriMeterVisAtt)

    logicWorld.SetVisAttributes(G4VisAttributes.GetInvisible())

    targetVisAtt = G4VisAttributes(G4Colour(1, 0 ,0))
    targetVisAtt.SetVisibility(True)
    logicTarget.SetVisAttributes(targetVisAtt)

    self.fScoringVolume = logicTarget

    return physWorld

    def ConstructSDandField(self):

        fieldValue = G4ThreeVector()
        self.fMagFieldMessenger = G4GlobalMagFieldMessenger(fieldValue)
        self.fMagFieldMessenger.SetVerboseLevel(1)

class PrimaryGeneratorAction(G4VUserPrimaryGeneratorAction):

def __init__(self):
    super().__init__()
    nofParticles = 1
    self.fParticleGun = G4ParticleGun(nofParticles)

    particleDefinition = G4ParticleTable.GetParticleTable().FindParticle("gamma")
    self.fParticleGun.SetParticleDefinition(particleDefinition)
    self.fParticleGun.SetParticleMomentumDirection(G4ThreeVector(1, 0, 0))
    self.fParticleGun.SetParticleEnergy(0.1*MeV)

def GeneratePrimaries(self, anEvent):
    # This function is called at the begining of an event

    # In order to avoid dependence of PrimaryGeneratorAction
    # on DetectorConstruction class we get world volume
    # from G4LogicalVolumeStore
    worldZHalfLength = 0
    logicWorld = G4LogicalVolumeStore.GetInstance().GetVolume("World")

    worldBox = None
    if logicWorld != None:
        worldBox = logicWorld.GetSolid()

    if worldBox != None:
        worldZHalfLength = worldBox.GetZHalfLength()
    else:
        msg = "World volume of box shape not found."
        msg += "Perhaps you have changed geometry."
        msg += "The gun will be place in the center."
        G4Exception("PrimaryGeneratorAction::GeneratePrimaries()",
                    "MyCode0002", JustWarning, msg)

    self.fParticleGun.SetParticlePosition(G4ThreeVector(0, 0, -worldZHalfLength))
    self.fParticleGun.GeneratePrimaryVertex(anEvent)

class EventAction(G4UserEventAction): def BeginOfEventAction(self, event):

initialisation per event

    self.fEnergyAbs = 0
    self.fEnergyGap = 0
    self.fTrackLAbs = 0
    self.fTrackLGap = 0

def EndOfEventAction(self, event):
    # Accumulate statistics

    # get analysis manager
    analysisManager = G4AnalysisManager.Instance()

    # fill histograms
    analysisManager.FillH1(0, self.fEnergyAbs)
    analysisManager.FillH1(1, self.fEnergyGap)
    analysisManager.FillH1(2, self.fTrackLAbs)
    analysisManager.FillH1(3, self.fTrackLGap)

    # fill ntuple
    analysisManager.FillNtupleDColumn(0, self.fEnergyAbs)
    analysisManager.FillNtupleDColumn(1, self.fEnergyGap)
    analysisManager.FillNtupleDColumn(2, self.fTrackLAbs)
    analysisManager.FillNtupleDColumn(3, self.fTrackLGap)
    analysisManager.AddNtupleRow()

    # Print per event (modulo n)
    eventID = event.GetEventID()
    printModulo = G4RunManager.GetRunManager().GetPrintProgress()
    if printModulo > 0 and eventID % printModulo == 0:
        print("---> End of event:", eventID)
        print("   Absorber: total energy:", G4BestUnit(self.fEnergyAbs, "Energy"), end="")
        print("       total track length:", G4BestUnit(self.fTrackLAbs, "Length"))
        print("        Gap: total energy:", G4BestUnit(self.fEnergyGap, "Energy"), end="")
        print("       total track length:", G4BestUnit(self.fTrackLGap, "Length"))

def AddAbs(self, de,  dl):
    self.fEnergyAbs += de
    self.fTrackLAbs += dl

def AddGap(self, de, dl):
    self.fEnergyGap += de
    self.fTrackLGap += dl

class SteppingAction(G4UserSteppingAction): def init(self, detectorConstruction, eventAction): super().init() self.fDetConstruction = detectorConstruction self.fEventAction = eventAction

def UserSteppingAction(self, step):
    # Collect energy and track length step by step

    # get volume of the current step
    volume = step.GetPreStepPoint().GetTouchable().GetVolume()

    # energy deposit
    edep = step.GetTotalEnergyDeposit()

    # step length
    stepLength = 0
    if step.GetTrack().GetDefinition().GetPDGCharge() != 0:
        stepLength = step.GetStepLength()

    if volume == self.fDetConstruction.fAbsorberPV:
        self.fEventAction.AddAbs(edep, stepLength)

    if volume == self.fDetConstruction.fGapPV:
        self.fEventAction.AddGap(edep, stepLength)

class RunAction(G4UserRunAction): def init(self): super().init()

    # set printing event number per each event
    G4RunManager.GetRunManager().SetPrintProgress(1)

    # Create analysis manager
    analysisManager = G4AnalysisManager.Instance()
    print("Using", analysisManager.GetType())

    analysisManager.SetVerboseLevel(1)
    analysisManager.SetNtupleMerging(True)

analysisManager.CreateH1("Eabs", "Edep in absorber", 100, 0, 800MeV) analysisManager.CreateH1("Egap", "Edep in gap", 100, 0, 100MeV) analysisManager.CreateH1("Labs", "trackL in absorber", 100, 0, 1m) analysisManager.CreateH1("Lgap", "trackL in gap", 100, 0, 50cm)

    # Creating ntuple
    analysisManager.CreateNtuple("B4", "Edep and TrackL")
    analysisManager.CreateNtupleDColumn("Eabs")
    analysisManager.CreateNtupleDColumn("Egap")
    analysisManager.CreateNtupleDColumn("Labs")
    analysisManager.CreateNtupleDColumn("Lgap")
    analysisManager.FinishNtuple()

def BeginOfRunAction(self, aRun):

analysisManager = G4AnalysisManager.Instance()

    # Open an output file
    fileName = "B4"
    analysisManager.OpenFile(fileName)

def EndOfRunAction(self, run):

analysisManager = G4AnalysisManager.Instance() if analysisManager.GetH1(1) != None: print("\n ----> print histograms statistic ", end="")

        if self.IsMaster():
            print("for the entire run \n")
        else:
            print("for the local thread \n")

        print(" EAbs : mean =", G4BestUnit(analysisManager.GetH1(0).mean(), "Energy"), end="")
        print(" rms =", G4BestUnit(analysisManager.GetH1(0).rms(),  "Energy"))

        print(" EGap : mean =", G4BestUnit(analysisManager.GetH1(1).mean(), "Energy"), end="")
        print(" rms =", G4BestUnit(analysisManager.GetH1(1).rms(),  "Energy"))

        print(" LAbs : mean = ", G4BestUnit(analysisManager.GetH1(2).mean(), "Length"), end="")
        print(" rms =", G4BestUnit(analysisManager.GetH1(2).rms(),  "Length"))

        print(" LGap : mean =", G4BestUnit(analysisManager.GetH1(3).mean(), "Length"), end="")
        print(" rms =", G4BestUnit(analysisManager.GetH1(3).rms(),  "Length"))

    # save histograms & ntuple
    analysisManager.Write()
    analysisManager.CloseFile()

class ActionInitialization(G4VUserActionInitialization): """ def BuildForMaster(self): self.SetUserAction(RunAction())

def Build(self):
    self.SetUserAction(PrimaryGeneratorAction())

    runAction = RunAction()
    self.SetUserAction(runAction)

    eventAction = EventAction(runAction)
    self.SetUserAction(eventAction)

    self.SetUserAction(SteppingAction(eventAction))

def BeginOfRunAction(self, run):
    # inform the runManager to save random number seed
    # G4RunManager.GetRunManager().SetRandomNumberStore(True)

    # Get analysis manager
    analysisManager = G4AnalysisManager.Instance()

    # Open an output file
    fileName = "deneme4"
    analysisManager.OpenFile(fileName)

def EndOfRunAction(self, run):
    analysisManager = G4AnalysisManager.Instance()

    analysisManager.Write()
    analysisManager.CloseFile()
"""
def __init__(self, detConstruction):
    super().__init__()
    self.fDetConstruction = detConstruction

def BuildForMaster(self):
    self.SetUserAction(RunAction())

def Build(self):
    self.SetUserAction(PrimaryGeneratorAction())
    self.SetUserAction(RunAction())
    eventAction = EventAction()
    self.SetUserAction(eventAction)
    self.SetUserAction(SteppingAction(self.fDetConstruction, eventAction))

def PrintUsage(): print("Usage: ", file=sys.stderr) print(" python exampleB4a.py [-m macro ] [-u UIsession] [-t nThreads]", file=sys.stderr) print(" note: -t option is available only for multi-threaded mode.", file=sys.stderr)

Evaluate arguments

if len(sys.argv) > 7: PrintUsage() sys.exit(1)

macro = "" session = ""

for i in range(1, len(sys.argv), 2): if sys.argv[i] == "-m": macro = sys.argv[i+1] elif sys.argv[i] == "-u": session = sys.argv[i+1] else: PrintUsage() sys.exit(1)

ui = None if len(macro) == 0: ui = G4UIExecutive(len(sys.argv), sys.argv, session)

runManager = G4RunManagerFactory.CreateRunManager(G4RunManagerType.Serial)

detConstruction = DetectorConstruction() runManager.SetUserInitialization(detConstruction)

physicsList = QBBC()

physicsList = FTFP_BERT() physicsList.SetVerboseLevel(1) runManager.SetUserInitialization(physicsList)

User action initialization

actionInitialization = ActionInitialization(detConstruction) runManager.SetUserInitialization(actionInitialization)

runManager.SetUserInitialization(ActionInitialization(detConstruction))

visManager = G4VisExecutive()

G4VisExecutive can take a verbosity argument - see /vis/verbose guidance.

visManager = G4VisExecutive("Quiet")

visManager.Initialize()

Get the User Interface manager

UImanager = G4UImanager.GetUIpointer()

Process macro or start UI session

if ui == None:

batch mode

command = "/control/execute "
fileName = sys.argv[1]
UImanager.ApplyCommand(command+fileName)

else:

interactive mode

UImanager.ApplyCommand("/control/execute init_vis.mac")
ui.SessionStart()`
saimegemenyucel commented 9 months ago

Sorry for my bad posting, it is my first time

HaarigerHarald commented 9 months ago

You could also attach a .py, .ipynb or similar.

saimegemenyucel commented 9 months ago

my code.zip I do not know what are those .prim files

HaarigerHarald commented 9 months ago

The placement caloriMeterPhys doesn't have a logical volume attached to it. You are passing a None. I think this is what you want to do:

caloriMeterPhys = G4PVPlacement(None,
                                G4ThreeVector(15.0*cm, 0, 0),
                                caloriMeterLogical,
                                "Detector",
                                logicTarget,
                                False,
                                0,
                                checkOverlaps)
saimegemenyucel commented 9 months ago

image Actually I had, but I also tried to add logicTarget. Result did not change unfortunately.