Compare commits

...

4 commits

Author SHA1 Message Date
Ava
324f479434 stand 2025-12-01, threads fehlerhaft 2025-12-01 17:50:50 +01:00
Ava
8809c2e595 define outputcolumns via macro 2025-12-01 14:16:36 +01:00
Ava
02f57fd562 particle type via macro 2025-12-01 13:46:34 +01:00
Ava
a25c6ee795 particle energy via macro 2025-12-01 13:27:21 +01:00
11 changed files with 342 additions and 162 deletions

View file

@ -1,16 +1,24 @@
#ifndef AACTIONINITIALIZATION_H
#define AACTIONINITIALIZATION_H
#pragma once
#include "G4VUserActionInitialization.hh"
class ARunAction;
class APrimaryGenerator;
class ASensitiveDetector;
class AActionInitialization : public G4VUserActionInitialization
{
public:
AActionInitialization();
AActionInitialization(ASensitiveDetector* sd);
virtual ~AActionInitialization();
virtual void Build() const override;
virtual void BuildForMaster() const override;
private:
ASensitiveDetector* fSD; // Pointer auf den Sensitive Detector
};
#endif

View file

@ -32,6 +32,8 @@ public:
virtual G4VPhysicalVolume* Construct();
void ConstructSDandField();
ASensitiveDetector* GetSensitiveDetector() const { return fSD; }
private:
std::vector<G4LogicalVolume*> detectorVolumes;

View file

@ -1,22 +1,25 @@
#ifndef AMESSENGER_HH
#define AMESSENGER_HH
#pragma once
#include "G4UImessenger.hh"
#include "G4UIcmdWithAString.hh"
#include "G4UIcmdWithADoubleAndUnit.hh"
class ASensitiveDetector;
class APrimaryGenerator;
class AMessenger : public G4UImessenger
{
public:
AMessenger(ASensitiveDetector* sd);
AMessenger(ASensitiveDetector* sd, APrimaryGenerator* generator = nullptr);
~AMessenger();
virtual void SetNewValue(G4UIcommand* cmd, G4String value) override;
void SetNewValue(G4UIcommand* cmd, G4String value) override;
private:
ASensitiveDetector* fSD;
G4UIcmdWithAString* fOutputFileCmd;
};
APrimaryGenerator* fGenerator; // Referenz zum Primärgenerator
#endif
G4UIcmdWithAString* fOutputFileCmd;
G4UIcmdWithAString* fParticleTypeCmd;
G4UIcmdWithADoubleAndUnit* fEnergyCmd;
G4UIcmdWithAString* fOutputColumnsCmd;
};

View file

@ -13,10 +13,15 @@ public:
APrimaryGenerator();
~APrimaryGenerator();
virtual void GeneratePrimaries(G4Event *);
virtual void GeneratePrimaries(G4Event * anEvent) override;
void SetParticleEnergy(G4double energy) { fParticleEnergy = energy; }
void SetParticleName(const G4String& name) { fParticleName = name; }
private:
G4ParticleGun *fParticleGun;
G4double fParticleEnergy;
G4String fParticleName;
};
#endif

View file

@ -8,14 +8,20 @@
#include "G4Track.hh"
#include "G4SystemOfUnits.hh"
#include "G4UnitsTable.hh"
#include "G4UImessenger.hh"
#include "G4UIcmdWithAString.hh"
#include <unordered_map>
#include <fstream>
#include <functional>
#include <map>
#include <vector>
class AMessenger;
//
// Makro für die OutputMap gleiche Logik wie beim Kollegen
//
#define OUTMAP [](const HitInfo& hit, std::ostream& out)
class ASensitiveDetector : public G4VSensitiveDetector
{
public:
@ -23,39 +29,89 @@ public:
~ASensitiveDetector();
virtual void Initialize(G4HCofThisEvent* hce) override;
virtual G4bool ProcessHits(G4Step* step, G4TouchableHistory* touch) override;
virtual void EndOfEvent(G4HCofThisEvent* hce) override;
virtual G4bool ProcessHits(G4Step* step, G4TouchableHistory* touchHist) override;
//Filename
void SetOutputFilename(const G4String& name) { fOutputFilename = name; }
void SetMessenger(AMessenger* messenger) { fMessenger = messenger; }
void RegisterVolumeName(G4LogicalVolume* lv, const std::string& cleanName) {
fCleanNames[lv] = cleanName;}
void SetOutputFilename(const G4String& f) { fOutputFilename = f; }
void SetOutputColumns(const std::vector<G4String>& cols) { fOutputColumns = cols; }
void SetMessenger(AMessenger* m) { fMessenger = m; }
private:
AMessenger* fMessenger = nullptr; // jetzt optional
G4String fOutputFilename = "outfiles/shower_setup.txt";
G4int fCurrentEventID = -1; // speichert die aktuelle Event-ID
AMessenger* fMessenger = nullptr;
G4String fOutputFilename;
G4int fCurrentEventID = -1;
struct HitInfo {
G4int eventID;
G4int trackID;
G4int parentID;
G4String particleName;
G4double particleEnergy;
G4String primaryParticleName;
G4double primaryEnergy;
G4double energyDeposited;
G4String detectorName;
G4String cleanDetectorName;
public:
//
// HitInfo enthält ALLE üblichen Parameter
//
struct HitInfo
{
// --- Basis-identifikation ---
G4int eventID; // Nummer des Events
G4int trackID; // Track ID
G4int parentID; // Parent Track ID
HitInfo()
: trackID(-1), parentID(-1), particleEnergy(0.), primaryEnergy(0.), energyDeposited(0.) {}
// --- Teilchen ---
G4String particleName; // Name des Teilchens
G4int pdg; // PDG-Code
G4double particleEnergy;// kinetische Energie am Vertex
// --- Primärinformation ---
G4String primaryName; // Name des Primärteilchens
G4double primaryEnergy; // Primärenergie
// --- Energien ---
G4double edep; // Deposited Energy
G4double ekin; // aktuelle kinetische Energie
// --- Positionen ---
G4ThreeVector prePos; // Position vor dem Step
G4ThreeVector postPos; // Position nach dem Step
// --- Richtung ---
G4ThreeVector preDir; // Richtung vor dem Step
G4ThreeVector postDir; // Richtung nach dem Step
// --- Zeiten ---
G4double globalTime; // absolute Zeit
G4double localTime; // lokale Zeit
// --- Step-Information ---
G4double stepLength; // Step-Länge
G4String processName; // Prozess, der diesen Step erzeugt hat
// --- Detektor ---
G4String detectorName; // voller Name
G4String cleanName; // aufgeräumter Name ohne _phys
HitInfo() :
eventID(0), trackID(0), parentID(0),
pdg(0),
particleEnergy(0.), primaryEnergy(0.),
edep(0.), ekin(0.),
globalTime(0.), localTime(0.),
stepLength(0.)
{}
};
G4String fDetectorName;
std::unordered_map<G4LogicalVolume*, std::string> fCleanNames;
std::unordered_map<G4int, HitInfo> fTrackHitMap; // TrackID → HitInfo
private:
//
// OutputMap: string → Lambda(hit, out)
//
using OutputLambda = std::function<void(const HitInfo&, std::ostream&)>;
std::map<std::string, OutputLambda> outputMap;
//
// interne Daten
//
std::vector<G4String> fOutputColumns;
std::unordered_map<G4int, HitInfo> fTrackHitMap;
//
// interne Hilfsfunktionen
//
void InitializeOutputMap();
};
#endif
#endif

17
macros/defrun.mac Normal file
View file

@ -0,0 +1,17 @@
/run/numberOfThreads 2
# --- Output Datei setzen ---
/my/outputFilename outfiles/test1000_100MeV_e
# --- Particle Type setzen ---
/my/particleType e-
# --- Particle Energy setzen ---
/my/particleEnergy 100 MeV
# --- Spalten setzen (ohne Kommata) ---
/my/outputColumns event track parent part pdg ekin edep x y z dir step proc det
/run/initialize
/run/beamOn 10

6
sim.cc
View file

@ -26,10 +26,12 @@ int main(int argc, char** argv)
runManager->SetUserInitialization(new APhysicsList());
// Detector
runManager->SetUserInitialization(new ADetectorConstruction());
ADetectorConstruction* detector = new ADetectorConstruction();
runManager->SetUserInitialization(detector);
// Actions
runManager->SetUserInitialization(new AActionInitialization());
runManager->SetUserInitialization(new AActionInitialization(detector->GetSensitiveDetector()));
// Visualization
G4VisManager* visManager = new G4VisExecutive();

View file

@ -1,9 +1,10 @@
#include "AActionInitialization.hh"
#include "AMessenger.hh"
#include "APrimaryGenerator.hh"
#include "ARunAction.hh"
AActionInitialization::AActionInitialization()
: G4VUserActionInitialization()
AActionInitialization::AActionInitialization(ASensitiveDetector* sd)
: G4VUserActionInitialization(), fSD(sd)
{
}
@ -12,7 +13,10 @@ AActionInitialization::~AActionInitialization() {}
void AActionInitialization::Build() const
{
// Primärgenerator
SetUserAction(new APrimaryGenerator());
APrimaryGenerator* primaryGen = new APrimaryGenerator();
SetUserAction(primaryGen);
AMessenger* messenger = new AMessenger(fSD, primaryGen);
// RunAction
SetUserAction(new ARunAction());

View file

@ -1,18 +1,42 @@
#include "AMessenger.hh"
#include "ASensitiveDetector.hh"
#include "APrimaryGenerator.hh"
#include "G4SystemOfUnits.hh"
AMessenger::AMessenger(ASensitiveDetector* sd)
: fSD(sd)
AMessenger::AMessenger(ASensitiveDetector* sd, APrimaryGenerator* generator)
: fSD(sd), fGenerator(generator)
{
// eigenes UI Kommando
fOutputFileCmd = new G4UIcmdWithAString("/my/outputFilename", this);
fOutputFileCmd->SetGuidance("Set output filename for SD text output.");
fOutputFileCmd->SetParameterName("filename", false);
// Kommando für Teilchentyp
fParticleTypeCmd = new G4UIcmdWithAString("/my/particleType", this);
fParticleTypeCmd->SetGuidance("Set particle type (e-, e+, gamma, proton, neutron).");
fParticleTypeCmd->SetParameterName("particle", false);
fParticleTypeCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
// Kommando für Teilchenenergie
fEnergyCmd = new G4UIcmdWithADoubleAndUnit("/my/particleEnergy", this);
fEnergyCmd->SetGuidance("Set particle energy.");
fEnergyCmd->SetParameterName("energy", false);
fEnergyCmd->SetUnitCategory("Energy");
fEnergyCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
//Kommando für Parameterreihenfolge
fOutputColumnsCmd = new G4UIcmdWithAString("/my/outputColumns", this);
fOutputColumnsCmd->SetGuidance("Set output columns and order (comma separated). E.g.: particle,Epart,Dname,Edep");
fOutputColumnsCmd->SetParameterName("columns", false);
fOutputColumnsCmd->AvailableForStates(G4State_PreInit, G4State_Idle);
}
AMessenger::~AMessenger()
{
delete fOutputFileCmd;
delete fParticleTypeCmd;
delete fEnergyCmd;
delete fOutputColumnsCmd;
}
void AMessenger::SetNewValue(G4UIcommand* cmd, G4String value)
@ -21,4 +45,24 @@ void AMessenger::SetNewValue(G4UIcommand* cmd, G4String value)
{
fSD->SetOutputFilename(value);
}
else if(cmd == fEnergyCmd && fGenerator)
{
fGenerator->SetParticleEnergy(fEnergyCmd->GetNewDoubleValue(value));
}
else if (cmd == fParticleTypeCmd && fGenerator)
{
fGenerator->SetParticleName(value);
}
else if(cmd == fOutputColumnsCmd && fSD)
{
std::vector<G4String> cols;
std::istringstream iss(value);
G4String token;
while(iss >> token)
{
cols.push_back(token);
}
fSD->SetOutputColumns(cols);
}
}

View file

@ -3,29 +3,16 @@
APrimaryGenerator::APrimaryGenerator()
{
fParticleGun = new G4ParticleGun(1); // 1 particles per event
fParticleEnergy = 100.*MeV; // Default-Energie
fParticleName = "e-"; //default
//Particle position
G4double x = 0. * m;
G4double y = 0. * m;
G4double z = -1. * cm;
G4ThreeVector pos(x, y, z);
//Particle direction
G4double px = 0.;
G4double py = 0.;
G4double pz = 1.;
G4ThreeVector mom(px,py,pz);
//Particle type
G4ParticleTable *particleTable = G4ParticleTable::GetParticleTable();
// possibilities: "e+","e-","gamma", ,"proton", "neutron"
G4ParticleDefinition *particle = particleTable->FindParticle("e-");
G4ThreeVector pos(0.,0.,-1.*cm);
G4ThreeVector mom(0.,0.,1.);
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* particle = particleTable->FindParticle("e-");
fParticleGun->SetParticlePosition(pos);
fParticleGun->SetParticleMomentumDirection(mom);
fParticleGun->SetParticleEnergy(100 * MeV); //or GeV
fParticleGun->SetParticleDefinition(particle);
}
@ -36,6 +23,15 @@ APrimaryGenerator::~APrimaryGenerator()
void APrimaryGenerator::GeneratePrimaries(G4Event *anEvent)
{
//Create vertex
G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
G4ParticleDefinition* particle = particleTable->FindParticle(fParticleName);
if (!particle)
{
G4cerr << "Unknown particle: " << fParticleName << ". Using e- as default." << G4endl;
particle = particleTable->FindParticle("e-");
}
fParticleGun->SetParticleDefinition(particle);
fParticleGun->SetParticleEnergy(fParticleEnergy); // dynamisch vom Macro
fParticleGun->GeneratePrimaryVertex(anEvent);
}

View file

@ -1,131 +1,174 @@
#include "ASensitiveDetector.hh"
#include "AMessenger.hh"
#include "G4RunManager.hh"
#include "G4Event.hh"
#include "G4PrimaryVertex.hh"
#include "G4PrimaryParticle.hh"
#include "G4SystemOfUnits.hh"
#include "G4RunManager.hh"
#include "G4Threading.hh"
#include <filesystem>
#include <iostream>
ASensitiveDetector::ASensitiveDetector(G4String name)
: G4VSensitiveDetector(name), fDetectorName(name)
: G4VSensitiveDetector(name), fOutputFilename("outfiles/shower_setup")
{
fOutputFilename = "outfiles/shower_setup.txt";
std::filesystem::create_directories("outfiles");
InitializeOutputMap();
if(fOutputColumns.empty()) {
fOutputColumns = {"event", "track", "parent", "part", "edep", "det"};
}
}
ASensitiveDetector::~ASensitiveDetector()
{
ASensitiveDetector::~ASensitiveDetector() {}
}
void ASensitiveDetector::Initialize(G4HCofThisEvent* /*hce*/)
void ASensitiveDetector::Initialize(G4HCofThisEvent*)
{
fTrackHitMap.clear();
}
G4bool ASensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory* /*touchHist*/)
void ASensitiveDetector::InitializeOutputMap()
{
outputMap = {
{"event", OUTMAP { out << hit.eventID; }},
{"track", OUTMAP { out << hit.trackID; }},
{"parent", OUTMAP { out << hit.parentID; }},
{"part", OUTMAP { out << hit.particleName; }},
{"pdg", OUTMAP { out << hit.pdg; }},
{"ekin", OUTMAP { out << hit.ekin / MeV; }},
{"edep", OUTMAP { out << hit.edep / MeV; }},
{"x", OUTMAP { out << hit.postPos.x() / mm; }},
{"y", OUTMAP { out << hit.postPos.y() / mm; }},
{"z", OUTMAP { out << hit.postPos.z() / mm; }},
{"pos", OUTMAP { out << hit.postPos.x()/mm<<" "<<hit.postPos.y()/mm<<" "<<hit.postPos.z()/mm; }},
{"dx", OUTMAP { out << hit.postDir.x(); }},
{"dy", OUTMAP { out << hit.postDir.y(); }},
{"dz", OUTMAP { out << hit.postDir.z(); }},
{"dir", OUTMAP { out << hit.postDir.x()<<" "<<hit.postDir.y()<<" "<<hit.postDir.z(); }},
{"global_t", OUTMAP { out << hit.globalTime; }},
{"local_t", OUTMAP { out << hit.localTime; }},
{"step", OUTMAP { out << hit.stepLength / mm; }},
{"proc", OUTMAP { out << hit.processName; }},
{"det", OUTMAP { out << hit.cleanName; }}
};
}
G4bool ASensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*)
{
if(!step) return false;
G4Track* trk = step->GetTrack();
if(!trk) return false;
const G4Event* event = G4RunManager::GetRunManager()->GetCurrentEvent();
if(!event) return false;
if (!event) return false;
// Event-ID direkt hier setzen
fCurrentEventID = event->GetEventID();
G4int id = trk->GetTrackID();
G4Track* track = step->GetTrack();
G4int trackID = track->GetTrackID();
HitInfo& hit = fTrackHitMap[id];
auto it = fTrackHitMap.find(trackID);
if (it == fTrackHitMap.end())
{
HitInfo hit;
hit.eventID = fCurrentEventID;
hit.trackID = trackID;
hit.parentID = track->GetParentID();
hit.particleName = track->GetDefinition()->GetParticleName();
hit.particleEnergy = track->GetVertexKineticEnergy();
hit.detectorName = step->GetPreStepPoint()->GetTouchable()->GetVolume()->GetName();
std::string detNameStr = hit.detectorName;
size_t pos = detNameStr.find("_phys");
if (pos != std::string::npos)
detNameStr = detNameStr.substr(0, pos); // nur Basisname
hit.cleanDetectorName = detNameStr;
// Primary particle info
G4PrimaryVertex* pv = event->GetPrimaryVertex();
if (pv)
{
G4PrimaryParticle* pp = pv->GetPrimary();
if (pp && pp->GetG4code())
hit.primaryParticleName = pp->GetG4code()->GetParticleName();
if (pp)
hit.primaryEnergy = pp->GetKineticEnergy();
}
hit.energyDeposited = step->GetTotalEnergyDeposit();
fTrackHitMap[trackID] = hit;
// Basisdaten
hit.eventID = fCurrentEventID;
hit.trackID = id;
hit.parentID = trk->GetParentID();
if(trk->GetDefinition()) {
hit.particleName = trk->GetDefinition()->GetParticleName();
hit.pdg = trk->GetDefinition()->GetPDGEncoding();
}
else
{
it->second.energyDeposited += step->GetTotalEnergyDeposit();
hit.ekin = trk->GetKineticEnergy();
// Primärteilchen (sicher prüfen)
if(auto vertex = event->GetPrimaryVertex()) {
if(auto prim = vertex->GetPrimary()) {
if(auto code = prim->GetG4code()) {
hit.primaryName = code->GetParticleName();
hit.primaryEnergy = prim->GetKineticEnergy();
}
}
}
// Energie aufsummieren
hit.edep += step->GetTotalEnergyDeposit();
// PreStep/PostStep prüfen
auto pre = step->GetPreStepPoint();
auto post = step->GetPostStepPoint();
if(pre && post) {
hit.prePos = pre->GetPosition();
hit.postPos = post->GetPosition();
hit.preDir = pre->GetMomentumDirection();
hit.postDir = post->GetMomentumDirection();
hit.globalTime = post->GetGlobalTime();
hit.localTime = post->GetLocalTime();
hit.stepLength = step->GetStepLength();
hit.processName = post->GetProcessDefinedStep() ? post->GetProcessDefinedStep()->GetProcessName() : "none";
// Detektorname prüfen
if(pre->GetTouchable()) {
auto vol = pre->GetTouchable()->GetVolume();
if(vol) {
auto full = vol->GetName();
hit.detectorName = full;
auto pos = full.find("_phys");
hit.cleanName = (pos == std::string::npos ? full : full.substr(0, pos));
} else {
hit.detectorName = "none";
hit.cleanName = "none";
}
} else {
hit.detectorName = "none";
hit.cleanName = "none";
}
}
return true;
}
void ASensitiveDetector::EndOfEvent(G4HCofThisEvent* /*hce*/)
void ASensitiveDetector::EndOfEvent(G4HCofThisEvent*)
{
if (fOutputFilename.empty()) fOutputFilename = "outfiles/shower_setup.txt";
std::filesystem::create_directories("outfiles");
// Thread-ID für Datei
G4int tid = G4Threading::G4GetThreadId();
std::ofstream out(fOutputFilename, std::ios::app);
if (!out.is_open()) {
G4cerr << "Cannot open file " << fOutputFilename << "!" << G4endl;
return;
// Thread-local Datei
thread_local std::ofstream out;
// Datei öffnen, wenn noch nicht offen
if(!out.is_open()) {
std::filesystem::create_directories("outfiles");
std::string threadFile = fOutputFilename + "_" + std::to_string(tid) + ".hits";
out.open(threadFile, std::ios::app);
if(!out.is_open()) {
G4cerr << "Cannot open output file: " << threadFile << G4endl;
return;
}
// Header einmal pro Datei
for(size_t i=0; i<fOutputColumns.size(); ++i) {
out << fOutputColumns[i];
if(i+1 < fOutputColumns.size()) out << "\t";
}
out << "\n";
}
static bool headerWritten = false;
if (!headerWritten) {
out << "eventID\ttrackID\tparentID\tparticle\tEpart\tPrimparticle\tEprim\tDname\tEdep\n";
headerWritten = true;
// Hits schreiben
for(auto& [id, hit] : fTrackHitMap) {
for(size_t i=0; i<fOutputColumns.size(); ++i) {
const auto& col = fOutputColumns[i];
if(outputMap.count(col))
outputMap[col](hit, out);
else
out << "NA";
if(i+1 < fOutputColumns.size()) out << "\t";
}
out << "\n";
}
// Hits nur für das aktuelle Event
std::vector<HitInfo> hits;
for (const auto& [trackID, hit] : fTrackHitMap) {
hits.push_back(hit);
}
// Sortieren: parentID -> trackID
std::sort(hits.begin(), hits.end(), [](const HitInfo& a, const HitInfo& b) {
if (a.parentID != b.parentID)
return a.parentID < b.parentID;
return a.trackID < b.trackID;
});
// Schreiben
for (const auto& hit : hits)
{
out << hit.eventID
<< "\t" << hit.trackID
<< "\t" << hit.parentID
<< "\t" << hit.particleName
<< "\t" << hit.particleEnergy / MeV
<< "\t" << hit.primaryParticleName
<< "\t" << hit.primaryEnergy / MeV
<< "\t" << hit.cleanDetectorName
<< "\t" << hit.energyDeposited / MeV
<< "\n";
}
out.close();
// Map leeren für nächstes Event
// Map zurücksetzen
fTrackHitMap.clear();
}