From 085a7cbcc7023da7e0f733e8cac2c0ca83657c47 Mon Sep 17 00:00:00 2001 From: Luis Obis Date: Thu, 15 Dec 2022 11:53:53 +0100 Subject: [PATCH 01/24] fail when process is not understood --- src/DataModel.cxx | 19 +++++++++++++++---- 1 file changed, 15 insertions(+), 4 deletions(-) diff --git a/src/DataModel.cxx b/src/DataModel.cxx index da441ee4..7a75fc22 100644 --- a/src/DataModel.cxx +++ b/src/DataModel.cxx @@ -170,14 +170,25 @@ void TRestGeant4Hits::InsertStep(const G4Step* step) { metadata->fGeant4PhysicsInfo.InsertParticleName(particleID, particleName); const auto process = step->GetPostStepPoint()->GetProcessDefinedStep(); - G4String processName = "Init"; - G4String processTypeName = "Init"; - Int_t processID = 0; - if (track->GetCurrentStepNumber() != 0) { + G4String processName; + G4String processTypeName; + Int_t processID; + if (process != nullptr) { // 0 = Init step (G4SteppingVerbose) process is not defined for this step processName = process->GetProcessName(); processTypeName = G4VProcess::GetProcessTypeName(process->GetProcessType()); processID = TRestGeant4PhysicsInfo::GetProcessIDFromGeant4Process(process); + } else { + if (track->GetCurrentStepNumber() == 0) { + processName = "Init"; + processTypeName = "Init"; + processID = 0; + } else { + cerr << "Cannot identify the process of track: " << track->GetTrackID() << " (" + << track->GetParticleDefinition()->GetParticleName() << ") " + << " at step number " << track->GetCurrentStepNumber() << endl; + // exit(1); + } } if (kill) { From 0201f6bf4601ecc8be06e2458df8f1292d8c8a51 Mon Sep 17 00:00:00 2001 From: Luis Obis Date: Thu, 15 Dec 2022 11:54:13 +0100 Subject: [PATCH 02/24] minor cleanup --- include/TrackingAction.h | 1 - src/DetectorConstruction.cxx | 4 ++-- src/TrackingAction.cxx | 5 ----- 3 files changed, 2 insertions(+), 8 deletions(-) diff --git a/include/TrackingAction.h b/include/TrackingAction.h index ae2d2f3f..0b56b5bd 100644 --- a/include/TrackingAction.h +++ b/include/TrackingAction.h @@ -17,7 +17,6 @@ class SimulationManager; class TrackingAction : public G4UserTrackingAction { public: TrackingAction(SimulationManager*); - ~TrackingAction(); virtual void PreUserTrackingAction(const G4Track*); virtual void PostUserTrackingAction(const G4Track*); diff --git a/src/DetectorConstruction.cxx b/src/DetectorConstruction.cxx index 2fd554ed..14063a27 100644 --- a/src/DetectorConstruction.cxx +++ b/src/DetectorConstruction.cxx @@ -296,7 +296,7 @@ void DetectorConstruction::ConstructSDandField() { } void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume* world) { - auto detector = (DetectorConstruction*)G4RunManager::GetRunManager()->GetUserDetectorConstruction(); + const auto detector = (DetectorConstruction*)G4RunManager::GetRunManager()->GetUserDetectorConstruction(); TRestGeant4Metadata* restG4Metadata = detector->fSimulationManager->GetRestMetadata(); const size_t n = int(world->GetLogicalVolume()->GetNoDaughters()); @@ -305,7 +305,7 @@ void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume* w if (i == n) { volume = const_cast(world); } else { - volume = world->GetLogicalVolume()->GetDaughter(i); + volume = world->GetLogicalVolume()->GetDaughter(G4int(i)); } TString namePhysical = (TString)volume->GetName(); if (fGdmlNewPhysicalNames.size() > i) { diff --git a/src/TrackingAction.cxx b/src/TrackingAction.cxx index 93c5d963..7799557b 100644 --- a/src/TrackingAction.cxx +++ b/src/TrackingAction.cxx @@ -2,11 +2,8 @@ #include "TrackingAction.h" -#include #include -#include #include -#include #include "SimulationManager.h" @@ -15,8 +12,6 @@ using namespace std; TrackingAction::TrackingAction(SimulationManager* simulationManager) : G4UserTrackingAction(), fSimulationManager(simulationManager) {} -TrackingAction::~TrackingAction() {} - void TrackingAction::PreUserTrackingAction(const G4Track* track) { fSimulationManager->GetOutputManager()->RecordTrack(track); } From 20efca0b5909a5b702ee00d00bf9f176060a8de2 Mon Sep 17 00:00:00 2001 From: Luis Obis Date: Thu, 15 Dec 2022 11:54:25 +0100 Subject: [PATCH 03/24] add prototype for biasing example --- examples/13.IAXO/GammaBiasing.rml | 65 +++++++++++++++++++++++++++++++ 1 file changed, 65 insertions(+) create mode 100644 examples/13.IAXO/GammaBiasing.rml diff --git a/examples/13.IAXO/GammaBiasing.rml b/examples/13.IAXO/GammaBiasing.rml new file mode 100644 index 00000000..d6f2faa9 --- /dev/null +++ b/examples/13.IAXO/GammaBiasing.rml @@ -0,0 +1,65 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + From aef3b204032afea1cf5b49cbf0608cac541933da Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Wed, 22 Feb 2023 11:04:17 +0100 Subject: [PATCH 04/24] biasing templated based on example 4 (Geant4) compiles --- include/GammaBiasingOperation.h | 47 ++++++++++++ include/GammaBiasingOperator.h | 48 ++++++++++++ src/DetectorConstruction.cxx | 129 ++++++++++++++++++-------------- src/GammaBiasingOperation.cxx | 104 +++++++++++++++++++++++++ src/GammaBiasingOperator.cxx | 52 +++++++++++++ 5 files changed, 324 insertions(+), 56 deletions(-) create mode 100644 include/GammaBiasingOperation.h create mode 100644 include/GammaBiasingOperator.h create mode 100644 src/GammaBiasingOperation.cxx create mode 100644 src/GammaBiasingOperator.cxx diff --git a/include/GammaBiasingOperation.h b/include/GammaBiasingOperation.h new file mode 100644 index 00000000..3579a03d --- /dev/null +++ b/include/GammaBiasingOperation.h @@ -0,0 +1,47 @@ +// +// Created by lobis on 2/22/2023. +// + +#ifndef REST_GAMMABIASINGOPERATION_H +#define REST_GAMMABIASINGOPERATION_H + +#include "G4VBiasingOperation.hh" +#include "G4ParticleChange.hh" +#include "G4BiasingProcessInterface.hh" + +class GammaBiasingOperation : public G4VBiasingOperation { +public: + GammaBiasingOperation(G4String name); + + virtual ~GammaBiasingOperation(); + +public: + virtual const G4VBiasingInteractionLaw * + ProvideOccurenceBiasingInteractionLaw(const G4BiasingProcessInterface *, + G4ForceCondition &) { return 0; } + + virtual G4VParticleChange *ApplyFinalStateBiasing(const G4BiasingProcessInterface *, + const G4Track *, + const G4Step *, + G4bool &); + + virtual G4double DistanceToApplyOperation(const G4Track *, + G4double, + G4ForceCondition *) { return DBL_MAX; } + + virtual G4VParticleChange *GenerateBiasingFinalState(const G4Track *, + const G4Step *) { return 0; } + + +public: + void SetSplittingFactor(G4int splittingFactor) { fSplittingFactor = splittingFactor; } + + G4int GetSplittingFactor() const { return fSplittingFactor; } + +private: + G4int fSplittingFactor; + G4ParticleChange fParticleChange; +}; + + +#endif //REST_GAMMABIASINGOPERATION_H diff --git a/include/GammaBiasingOperator.h b/include/GammaBiasingOperator.h new file mode 100644 index 00000000..c0ce5253 --- /dev/null +++ b/include/GammaBiasingOperator.h @@ -0,0 +1,48 @@ +// +// Created by lobis on 2/22/2023. +// + +#ifndef REST_GAMMABIASINGOPERATOR_H +#define REST_GAMMABIASINGOPERATOR_H + +#include "G4VBiasingOperator.hh" + +class GammaBiasingOperation; + +class GammaBiasingOperator : public G4VBiasingOperator { +public: + GammaBiasingOperator(); + + virtual ~GammaBiasingOperator() {} + +public: + virtual void StartRun(); + + virtual void StartTracking(const G4Track *track); + +private: + virtual G4VBiasingOperation * + ProposeNonPhysicsBiasingOperation(const G4Track *, + const G4BiasingProcessInterface *) { return 0; } + + virtual G4VBiasingOperation * + ProposeOccurenceBiasingOperation(const G4Track *, + const G4BiasingProcessInterface *) { return 0; } + + virtual G4VBiasingOperation * + ProposeFinalStateBiasingOperation(const G4Track *track, + const G4BiasingProcessInterface *callingProcess); + +private: + using G4VBiasingOperator::OperationApplied; + +private: + GammaBiasingOperation *fBremSplittingOperation; + G4int fSplittingFactor; + G4bool fBiasPrimaryOnly; + G4bool fBiasOnlyOnce; + G4int fNInteractions; +}; + + +#endif //REST_GAMMABIASINGOPERATOR_H diff --git a/src/DetectorConstruction.cxx b/src/DetectorConstruction.cxx index 99c6044a..193550f9 100644 --- a/src/DetectorConstruction.cxx +++ b/src/DetectorConstruction.cxx @@ -18,20 +18,21 @@ #include #include "SimulationManager.h" +#include "GammaBiasingOperator.h" using namespace std; using namespace TRestGeant4PrimaryGeneratorTypes; -DetectorConstruction::DetectorConstruction(SimulationManager* simulationManager) - : fSimulationManager(simulationManager) { +DetectorConstruction::DetectorConstruction(SimulationManager *simulationManager) + : fSimulationManager(simulationManager) { G4cout << "Detector Construction" << G4endl; fGdmlParser = new G4GDMLParser(); } DetectorConstruction::~DetectorConstruction() { delete fGdmlParser; } -G4VPhysicalVolume* DetectorConstruction::Construct() { - TRestGeant4Metadata* restG4Metadata = fSimulationManager->GetRestMetadata(); +G4VPhysicalVolume *DetectorConstruction::Construct() { + TRestGeant4Metadata *restG4Metadata = fSimulationManager->GetRestMetadata(); cout << "Isotope table " << endl; cout << *(G4Isotope::GetIsotopeTable()) << endl; @@ -44,31 +45,31 @@ G4VPhysicalVolume* DetectorConstruction::Construct() { const auto startingPath = filesystem::current_path(); const auto [gdmlPath, gdmlToRead] = - TRestTools::SeparatePathAndName((string)restG4Metadata->GetGdmlFilename()); + TRestTools::SeparatePathAndName((string) restG4Metadata->GetGdmlFilename()); filesystem::current_path(gdmlPath); G4cout << "gdmlToRead: " << gdmlToRead << G4endl; fGdmlParser->Read(gdmlToRead, false); - G4VPhysicalVolume* worldVolume = fGdmlParser->GetWorldVolume(); + G4VPhysicalVolume *worldVolume = fGdmlParser->GetWorldVolume(); - const auto worldSolid = dynamic_cast(worldVolume->GetLogicalVolume()->GetSolid()); + const auto worldSolid = dynamic_cast(worldVolume->GetLogicalVolume()->GetSolid()); restG4Metadata->fGeant4PrimaryGeneratorInfo.fSpatialGeneratorWorldSize = { - worldSolid->GetXHalfLength(), worldSolid->GetYHalfLength(), worldSolid->GetZHalfLength()}; + worldSolid->GetXHalfLength(), worldSolid->GetYHalfLength(), worldSolid->GetZHalfLength()}; restG4Metadata->fGeant4GeometryInfo.InitializeOnDetectorConstruction(gdmlToRead, worldVolume); restG4Metadata->ReadDetector(); restG4Metadata->PrintMetadata(); // now we have detector info - const auto& geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); + const auto &geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); geometryInfo.Print(); // do some checks { // Check all physical volume names are unique - G4PhysicalVolumeStore* physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); + G4PhysicalVolumeStore *physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); set physicalVolumeNames; - vector::const_iterator physicalVolume; + vector::const_iterator physicalVolume; for (physicalVolume = physicalVolumeStore->begin(); physicalVolume != physicalVolumeStore->end(); physicalVolume++) { auto name = (*physicalVolume)->GetName(); @@ -83,9 +84,9 @@ G4VPhysicalVolume* DetectorConstruction::Construct() { } // Check all logical volume names are unique - G4LogicalVolumeStore* logicalVolumeStore = G4LogicalVolumeStore::GetInstance(); + G4LogicalVolumeStore *logicalVolumeStore = G4LogicalVolumeStore::GetInstance(); set logicalVolumeNames; - vector::const_iterator logicalVolume; + vector::const_iterator logicalVolume; for (logicalVolume = logicalVolumeStore->begin(); logicalVolume != logicalVolumeStore->end(); logicalVolume++) { auto name = (*logicalVolume)->GetName(); @@ -101,15 +102,15 @@ G4VPhysicalVolume* DetectorConstruction::Construct() { } filesystem::current_path(startingPath); - auto sensitiveVolume = (string)restG4Metadata->GetSensitiveVolume(); - G4VPhysicalVolume* physicalVolume = GetPhysicalVolume(sensitiveVolume); + auto sensitiveVolume = (string) restG4Metadata->GetSensitiveVolume(); + G4VPhysicalVolume *physicalVolume = GetPhysicalVolume(sensitiveVolume); if (physicalVolume == nullptr) { // sensitive volume was not found, perhaps the user specified a logical volume auto physicalVolumes = geometryInfo.GetAllPhysicalVolumesFromLogical(sensitiveVolume); if (physicalVolumes.size() == 1) { restG4Metadata->InsertSensitiveVolume( - geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolumes[0])); - sensitiveVolume = (string)restG4Metadata->GetSensitiveVolume(); + geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolumes[0])); + sensitiveVolume = (string) restG4Metadata->GetSensitiveVolume(); physicalVolume = GetPhysicalVolume(sensitiveVolume); } } @@ -123,20 +124,20 @@ G4VPhysicalVolume* DetectorConstruction::Construct() { Double_t my = restG4Metadata->GetMagneticField().Y() * tesla; Double_t mz = restG4Metadata->GetMagneticField().Z() * tesla; - G4MagneticField* magField = new G4UniformMagField(G4ThreeVector(mx, my, mz)); + G4MagneticField *magField = new G4UniformMagField(G4ThreeVector(mx, my, mz)); // G4FieldManager* localFieldMgr = new G4FieldManager(magField); - G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager(); + G4FieldManager *fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager(); fieldMgr->SetDetectorField(magField); fieldMgr->CreateChordFinder(magField); - G4LogicalVolume* volume = physicalVolume->GetLogicalVolume(); - G4Material* material = volume->GetMaterial(); + G4LogicalVolume *volume = physicalVolume->GetLogicalVolume(); + G4Material *material = volume->GetMaterial(); G4cout << "Sensitive volume properties:" << G4endl; G4cout << "\t- Material: " << material->GetName() << G4endl; G4cout << "\t- Temperature: " << material->GetTemperature() << " K" << G4endl; G4cout << "\t- Density: " << material->GetDensity() / (g / cm3) << " g/cm3" << G4endl; - const auto& primaryGeneratorInfo = restG4Metadata->GetGeant4PrimaryGeneratorInfo(); + const auto &primaryGeneratorInfo = restG4Metadata->GetGeant4PrimaryGeneratorInfo(); // Getting generation volume const auto fromVolume = primaryGeneratorInfo.GetSpatialGeneratorFrom(); if (fromVolume != "NO_SUCH_PARA") { @@ -145,11 +146,11 @@ G4VPhysicalVolume* DetectorConstruction::Construct() { cout << "Generator type: " << primaryGeneratorInfo.GetSpatialGeneratorType() << endl; const auto spatialGeneratorTypeEnum = - StringToSpatialGeneratorTypes(primaryGeneratorInfo.GetSpatialGeneratorType().Data()); + StringToSpatialGeneratorTypes(primaryGeneratorInfo.GetSpatialGeneratorType().Data()); if (spatialGeneratorTypeEnum == TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::VOLUME && primaryGeneratorInfo.GetSpatialGeneratorFrom() != "Not defined") { - G4VPhysicalVolume* pVol = GetPhysicalVolume(primaryGeneratorInfo.GetSpatialGeneratorFrom().Data()); + G4VPhysicalVolume *pVol = GetPhysicalVolume(primaryGeneratorInfo.GetSpatialGeneratorFrom().Data()); if (pVol == nullptr) { cout << "ERROR: The generator volume '" << primaryGeneratorInfo.GetSpatialGeneratorFrom() << "'was not found in the geometry" << endl; @@ -161,7 +162,7 @@ G4VPhysicalVolume* DetectorConstruction::Construct() { if (spatialGeneratorTypeEnum == TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::SURFACE || spatialGeneratorTypeEnum == TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::VOLUME) { restG4Metadata->fGeant4PrimaryGeneratorInfo.fSpatialGeneratorPosition = { - fGeneratorTranslation.x(), fGeneratorTranslation.y(), fGeneratorTranslation.z()}; + fGeneratorTranslation.x(), fGeneratorTranslation.y(), fGeneratorTranslation.z()}; } fGeneratorSolid = pVol->GetLogicalVolume()->GetSolid(); @@ -202,9 +203,9 @@ G4VPhysicalVolume* DetectorConstruction::Construct() { for (unsigned int id = 0; id < restG4Metadata->GetNumberOfActiveVolumes(); id++) { TString activeVolumeName = restG4Metadata->GetActiveVolumeName(id); - G4VPhysicalVolume* pVol = GetPhysicalVolume((G4String)activeVolumeName); + G4VPhysicalVolume *pVol = GetPhysicalVolume((G4String) activeVolumeName); if (pVol != nullptr) { - G4LogicalVolume* lVol = pVol->GetLogicalVolume(); + G4LogicalVolume *lVol = pVol->GetLogicalVolume(); if (restG4Metadata->GetMaxStepSize(activeVolumeName) > 0) { RESTInfo << "Setting maxStepSize of " << restG4Metadata->GetMaxStepSize(activeVolumeName) * mm << "mm for volume '" << activeVolumeName << "'" << RESTendl; @@ -220,15 +221,15 @@ G4VPhysicalVolume* DetectorConstruction::Construct() { return worldVolume; } -G4VPhysicalVolume* DetectorConstruction::GetPhysicalVolume(const G4String& physicalVolumeName) const { - G4PhysicalVolumeStore* physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); - TRestGeant4Metadata* restG4Metadata = fSimulationManager->GetRestMetadata(); - const auto& geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); - vector::const_iterator physicalVolume; +G4VPhysicalVolume *DetectorConstruction::GetPhysicalVolume(const G4String &physicalVolumeName) const { + G4PhysicalVolumeStore *physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); + TRestGeant4Metadata *restG4Metadata = fSimulationManager->GetRestMetadata(); + const auto &geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); + vector::const_iterator physicalVolume; for (physicalVolume = physicalVolumeStore->begin(); physicalVolume != physicalVolumeStore->end(); physicalVolume++) { auto name = (*physicalVolume)->GetName(); - auto alternativeName = (G4String)geometryInfo.GetAlternativeNameFromGeant4PhysicalName(name); + auto alternativeName = (G4String) geometryInfo.GetAlternativeNameFromGeant4PhysicalName(name); if (name == physicalVolumeName || alternativeName == physicalVolumeName) { return *physicalVolume; } @@ -238,20 +239,20 @@ G4VPhysicalVolume* DetectorConstruction::GetPhysicalVolume(const G4String& physi } void DetectorConstruction::ConstructSDandField() { - const TRestGeant4Metadata& metadata = *fSimulationManager->GetRestMetadata(); + const TRestGeant4Metadata &metadata = *fSimulationManager->GetRestMetadata(); - set logicalVolumesSelected; - for (const auto& userSensitiveVolume : metadata.GetSensitiveVolumes()) { + set logicalVolumesSelected; + for (const auto &userSensitiveVolume: metadata.GetSensitiveVolumes()) { // Each sensitive detector declaration in the RML should correspond to at least one logical volume - G4LogicalVolume* logicalVolume; + G4LogicalVolume *logicalVolume; // Check if user selected a Geant4 physical volume by name - G4VPhysicalVolume* physicalVolume = - G4PhysicalVolumeStore::GetInstance()->GetVolume(userSensitiveVolume.Data(), false); + G4VPhysicalVolume *physicalVolume = + G4PhysicalVolumeStore::GetInstance()->GetVolume(userSensitiveVolume.Data(), false); if (physicalVolume == nullptr) { const G4String geant4VolumeName = - metadata.GetGeant4GeometryInfo() - .GetGeant4PhysicalNameFromAlternativeName(userSensitiveVolume.Data()) - .Data(); + metadata.GetGeant4GeometryInfo() + .GetGeant4PhysicalNameFromAlternativeName(userSensitiveVolume.Data()) + .Data(); // Check if user selected a Geant4 physical volume by REST name (only relevant for assemblies) physicalVolume = G4PhysicalVolumeStore::GetInstance()->GetVolume(geant4VolumeName, false); } @@ -268,53 +269,69 @@ void DetectorConstruction::ConstructSDandField() { // Check if the user string matches a logical volume by expanding the input as a regex // This can have multiple hits const auto logicalVolumesMatchingExpression = - metadata.GetGeant4GeometryInfo().GetAllLogicalVolumesMatchingExpression(userSensitiveVolume); + metadata.GetGeant4GeometryInfo().GetAllLogicalVolumesMatchingExpression(userSensitiveVolume); if (logicalVolumesMatchingExpression.empty()) { cerr << "Detector construction error: could not find matching logical volume(s) for '" << userSensitiveVolume << "'" << endl; exit(1); } else { - for (const auto& logicalVolumeName : logicalVolumesMatchingExpression) { + for (const auto &logicalVolumeName: logicalVolumesMatchingExpression) { logicalVolumesSelected.insert( - G4LogicalVolumeStore::GetInstance()->GetVolume(logicalVolumeName.Data(), false)); + G4LogicalVolumeStore::GetInstance()->GetVolume(logicalVolumeName.Data(), false)); } } } } - G4SDManager* SDManager = G4SDManager::GetSDMpointer(); + G4SDManager *SDManager = G4SDManager::GetSDMpointer(); - for (G4LogicalVolume* logicalVolume : logicalVolumesSelected) { + for (G4LogicalVolume *logicalVolume: logicalVolumesSelected) { auto name = logicalVolume->GetName(); - G4VSensitiveDetector* sensitiveDetector = new SensitiveDetector(fSimulationManager, name); + G4VSensitiveDetector *sensitiveDetector = new SensitiveDetector(fSimulationManager, name); SDManager->AddNewDetector(sensitiveDetector); logicalVolume->SetSensitiveDetector(sensitiveDetector); auto region = new G4Region(name); logicalVolume->SetRegion(region); } + + // Biasing + const auto biasingLogicalVolumeName = "ShieldingVolume"; + G4LogicalVolume *biasingLogicalVolume = + G4LogicalVolumeStore::GetInstance()->GetVolume(biasingLogicalVolumeName, false); + if (biasingLogicalVolume == nullptr) { + cerr << "Detector construction error: could not find logical volume '" << biasingLogicalVolumeName << "'" + << endl; + exit(1); + } + + GammaBiasingOperator *bremSplittingOperator = new GammaBiasingOperator(); + bremSplittingOperator->AttachTo(biasingLogicalVolume); + G4cout << " Attaching biasing operator " << bremSplittingOperator->GetName() + << " to logical volume " << biasingLogicalVolume->GetName() + << G4endl; } -void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume* world) { - const auto detector = (DetectorConstruction*)G4RunManager::GetRunManager()->GetUserDetectorConstruction(); - TRestGeant4Metadata* restG4Metadata = detector->fSimulationManager->GetRestMetadata(); +void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume *world) { + const auto detector = (DetectorConstruction *) G4RunManager::GetRunManager()->GetUserDetectorConstruction(); + TRestGeant4Metadata *restG4Metadata = detector->fSimulationManager->GetRestMetadata(); const size_t n = int(world->GetLogicalVolume()->GetNoDaughters()); for (size_t i = 0; i < n + 1; i++) { // world is the + 1 - G4VPhysicalVolume* volume; + G4VPhysicalVolume *volume; if (i == n) { - volume = const_cast(world); + volume = const_cast(world); } else { volume = world->GetLogicalVolume()->GetDaughter(G4int(i)); } - TString namePhysical = (TString)volume->GetName(); + TString namePhysical = (TString) volume->GetName(); if (fGdmlNewPhysicalNames.size() > i) { // it has been filled fGeant4PhysicalNameToNewPhysicalNameMap[namePhysical] = fGdmlNewPhysicalNames[i]; } TString physicalNewName = GetAlternativeNameFromGeant4PhysicalName(namePhysical); - TString nameLogical = (TString)volume->GetLogicalVolume()->GetName(); - TString nameMaterial = (TString)volume->GetLogicalVolume()->GetMaterial()->GetName(); + TString nameLogical = (TString) volume->GetLogicalVolume()->GetName(); + TString nameMaterial = (TString) volume->GetLogicalVolume()->GetMaterial()->GetName(); auto position = volume->GetTranslation(); fPhysicalToLogicalVolumeMap[physicalNewName] = nameLogical; diff --git a/src/GammaBiasingOperation.cxx b/src/GammaBiasingOperation.cxx new file mode 100644 index 00000000..8bfd5fc4 --- /dev/null +++ b/src/GammaBiasingOperation.cxx @@ -0,0 +1,104 @@ + +#include "GammaBiasingOperation.h" +#include "GammaBiasingOperator.h" + +#include "G4ParticleChangeForLoss.hh" + +GammaBiasingOperation::GammaBiasingOperation(G4String name) + : G4VBiasingOperation(name), + fSplittingFactor(1), + fParticleChange() { +} + + +GammaBiasingOperation::~GammaBiasingOperation() {} + +G4VParticleChange * +GammaBiasingOperation:: +ApplyFinalStateBiasing(const G4BiasingProcessInterface *callingProcess, + const G4Track *track, + const G4Step *step, + G4bool &) { + + // -- Collect brem. process (wrapped process) final state: + G4VParticleChange *processFinalState = + callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); + + // -- if no splitting requested, let the brem. process to return directly its + // -- generated final state: + if (fSplittingFactor == 1) return processFinalState; + + // -- a special case here: the brem. process corrects for cross-section change + // -- over the step due to energy loss by sometimes "abandoning" the interaction, + // -- returning an unchanged incoming electron/positron. + // -- We respect this correction, and if no secondary is produced, its means this + // -- case is happening: + if (processFinalState->GetNumberOfSecondaries() == 0) return processFinalState; + + // -- Now start the biasing: + // -- - the electron state will be taken as the first one produced by the brem. + // -- process, hence the one stored in above processFinalState particle change. + // -- This state will be stored in our fParticleChange object. + // -- - the photon accompagnying the electron will be stored also this way. + // -- - we will then do fSplittingFactor - 1 call to the brem. process to collect + // -- fSplittingFactor - 1 additionnal gammas. All these will be stored in our + // -- fParticleChange object. + + // -- We called the brem. process above. Its concrete particle change is indeed + // -- a "G4ParticleChangeForLoss" object. We cast this particle change to access + // -- methods of the concrete G4ParticleChangeForLoss type: + G4ParticleChangeForLoss *actualParticleChange = + (G4ParticleChangeForLoss *) processFinalState; + + fParticleChange.Initialize(*track); + + // -- Store electron final state: + fParticleChange. + ProposeTrackStatus(actualParticleChange->GetTrackStatus()); + fParticleChange. + ProposeEnergy(actualParticleChange->GetProposedKineticEnergy()); + fParticleChange. + ProposeMomentumDirection(actualParticleChange->GetProposedMomentumDirection()); + + // -- Now deal with the gamma's: + // -- their common weight: + G4double gammaWeight = track->GetWeight() / fSplittingFactor; + + // -- inform we will have fSplittingFactor gamma's: + fParticleChange.SetNumberOfSecondaries(fSplittingFactor); + + // -- inform we take care of secondaries weight (otherwise these + // -- secondaries are by default given the primary weight). + fParticleChange.SetSecondaryWeightByProcess(true); + + // -- Store first gamma: + G4Track *gammaTrack = actualParticleChange->GetSecondary(0); + gammaTrack->SetWeight(gammaWeight); + fParticleChange.AddSecondary(gammaTrack); + // -- and clean-up the brem. process particle change: + actualParticleChange->Clear(); + + // -- now start the fSplittingFactor-1 calls to the brem. process to store each + // -- related gamma: + G4int nCalls = 1; + while (nCalls < fSplittingFactor) { + // ( note: we don't need to cast to actual type here, as methods for accessing + // secondary particles are from base class G4VParticleChange ) + processFinalState = callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); + if (processFinalState->GetNumberOfSecondaries() == 1) { + gammaTrack = processFinalState->GetSecondary(0); + gammaTrack->SetWeight(gammaWeight); + fParticleChange.AddSecondary(gammaTrack); + nCalls++; + } + // -- very rare special case: we ignore for now. + else if (processFinalState->GetNumberOfSecondaries() > 1) { + for (G4int i = 0; i < processFinalState->GetNumberOfSecondaries(); i++) + delete processFinalState->GetSecondary(i); + } + processFinalState->Clear(); + } + + // -- we are done: + return &fParticleChange; +} diff --git a/src/GammaBiasingOperator.cxx b/src/GammaBiasingOperator.cxx new file mode 100644 index 00000000..4cefb4bd --- /dev/null +++ b/src/GammaBiasingOperator.cxx @@ -0,0 +1,52 @@ +// +// Created by lobis on 2/22/2023. +// + +#include "GammaBiasingOperation.h" +#include "GammaBiasingOperator.h" + +#include "G4BiasingProcessInterface.hh" +#include "G4GenericMessenger.hh" + +GammaBiasingOperator::GammaBiasingOperator() + : G4VBiasingOperator("BremSplittingOperator"), + fSplittingFactor(1), + fBiasPrimaryOnly(true), + fBiasOnlyOnce(true) { + fBremSplittingOperation = new GammaBiasingOperation("BremSplittingOperation"); +} + + +void GammaBiasingOperator::StartRun() { + fBremSplittingOperation->SetSplittingFactor(fSplittingFactor); + G4cout << GetName() << " : starting run with brem. splitting factor = " + << fSplittingFactor; + if (fBiasPrimaryOnly) G4cout << ", biasing only primaries "; + else + G4cout << ", biasing primary and secondary tracks "; + if (fBiasOnlyOnce) G4cout << ", biasing only once per track "; + else + G4cout << ", biasing several times per track "; + G4cout << " . " << G4endl; +} + +void GammaBiasingOperator::StartTracking(const G4Track * /* track */ ) { + // -- reset the number of times the brem. splitting was applied: + fNInteractions = 0; +} + +G4VBiasingOperation * +GammaBiasingOperator::ProposeFinalStateBiasingOperation(const G4Track *track, + const G4BiasingProcessInterface * /* callingProcess */) { + // -- Check if biasing of primary particle only is requested. If so, and + // -- if particle is not a primary one, don't ask for biasing: + if (fBiasPrimaryOnly && (track->GetParentID() != 0)) return nullptr; + // -- Check if brem. splitting should be applied only once to the track, + // -- and if so, and if brem. splitting already occured, don't ask for biasing: + if (fBiasOnlyOnce && (fNInteractions > 0)) return nullptr; + + // -- Count the number of times the brem. splitting is applied: + fNInteractions++; + // -- Return the brem. splitting operation: + return fBremSplittingOperation; +} \ No newline at end of file From 4d886500c6920ca3c3a8060d0d9640bcc053475d Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Wed, 22 Feb 2023 11:49:18 +0100 Subject: [PATCH 05/24] working splitting in physics list --- include/PhysicsList.h | 30 ++++++++++++------- src/GammaBiasingOperation.cxx | 18 +++++------ src/GammaBiasingOperator.cxx | 5 ++-- src/PhysicsList.cxx | 56 ++++++++++++++++++++++------------- src/TrackingAction.cxx | 9 +++--- 5 files changed, 72 insertions(+), 46 deletions(-) diff --git a/include/PhysicsList.h b/include/PhysicsList.h index ee1cdb8b..1f56feef 100644 --- a/include/PhysicsList.h +++ b/include/PhysicsList.h @@ -10,33 +10,43 @@ #include "SimulationManager.h" + +class G4GenericBiasingPhysics; + class PhysicsList : public G4VModularPhysicsList { - public: +public: PhysicsList() = delete; - explicit PhysicsList(SimulationManager* simulationManager, TRestGeant4PhysicsLists* restPhysicsLists); + + explicit PhysicsList(SimulationManager *simulationManager, TRestGeant4PhysicsLists *restPhysicsLists); + ~PhysicsList() override; - protected: +protected: // Construct particle and physics virtual void InitializePhysicsLists(); void ConstructParticle() override; + void ConstructProcess() override; + void SetCuts() override; - private: - SimulationManager* fSimulationManager = nullptr; +private: + SimulationManager *fSimulationManager = nullptr; G4EmConfigurator fEmConfig; - G4VPhysicsConstructor* fEmPhysicsList = nullptr; + G4VPhysicsConstructor *fEmPhysicsList = nullptr; std::string fEmPhysicsListName; // Can be different from the output of GetPhysicsName - G4VPhysicsConstructor* fDecPhysicsList = nullptr; - G4VPhysicsConstructor* fRadDecPhysicsList = nullptr; - std::vector fHadronPhys; + G4VPhysicsConstructor *fDecPhysicsList = nullptr; + G4VPhysicsConstructor *fRadDecPhysicsList = nullptr; + std::vector fHadronPhys; + + G4GenericBiasingPhysics *fBiasingPhysicsList = nullptr; + + TRestGeant4PhysicsLists *fRestPhysicsLists = nullptr; - TRestGeant4PhysicsLists* fRestPhysicsLists = nullptr; }; #endif diff --git a/src/GammaBiasingOperation.cxx b/src/GammaBiasingOperation.cxx index 8bfd5fc4..893fe7ea 100644 --- a/src/GammaBiasingOperation.cxx +++ b/src/GammaBiasingOperation.cxx @@ -1,17 +1,19 @@ +#include + #include "GammaBiasingOperation.h" #include "GammaBiasingOperator.h" #include "G4ParticleChangeForLoss.hh" GammaBiasingOperation::GammaBiasingOperation(G4String name) - : G4VBiasingOperation(name), + : G4VBiasingOperation(std::move(name)), fSplittingFactor(1), fParticleChange() { } -GammaBiasingOperation::~GammaBiasingOperation() {} +GammaBiasingOperation::~GammaBiasingOperation() = default; G4VParticleChange * GammaBiasingOperation:: @@ -20,12 +22,9 @@ ApplyFinalStateBiasing(const G4BiasingProcessInterface *callingProcess, const G4Step *step, G4bool &) { - // -- Collect brem. process (wrapped process) final state: G4VParticleChange *processFinalState = callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); - // -- if no splitting requested, let the brem. process to return directly its - // -- generated final state: if (fSplittingFactor == 1) return processFinalState; // -- a special case here: the brem. process corrects for cross-section change @@ -39,16 +38,17 @@ ApplyFinalStateBiasing(const G4BiasingProcessInterface *callingProcess, // -- - the electron state will be taken as the first one produced by the brem. // -- process, hence the one stored in above processFinalState particle change. // -- This state will be stored in our fParticleChange object. - // -- - the photon accompagnying the electron will be stored also this way. + // -- - the photon accompanying the electron will be stored also this way. // -- - we will then do fSplittingFactor - 1 call to the brem. process to collect - // -- fSplittingFactor - 1 additionnal gammas. All these will be stored in our + // -- fSplittingFactor - 1 additional gammas. All these will be stored in our // -- fParticleChange object. // -- We called the brem. process above. Its concrete particle change is indeed // -- a "G4ParticleChangeForLoss" object. We cast this particle change to access // -- methods of the concrete G4ParticleChangeForLoss type: - G4ParticleChangeForLoss *actualParticleChange = - (G4ParticleChangeForLoss *) processFinalState; + + + auto actualParticleChange = (G4ParticleChangeForLoss *) processFinalState; fParticleChange.Initialize(*track); diff --git a/src/GammaBiasingOperator.cxx b/src/GammaBiasingOperator.cxx index 4cefb4bd..9f780a2e 100644 --- a/src/GammaBiasingOperator.cxx +++ b/src/GammaBiasingOperator.cxx @@ -10,8 +10,8 @@ GammaBiasingOperator::GammaBiasingOperator() : G4VBiasingOperator("BremSplittingOperator"), - fSplittingFactor(1), - fBiasPrimaryOnly(true), + fSplittingFactor(10), + fBiasPrimaryOnly(false), fBiasOnlyOnce(true) { fBremSplittingOperation = new GammaBiasingOperation("BremSplittingOperation"); } @@ -38,6 +38,7 @@ void GammaBiasingOperator::StartTracking(const G4Track * /* track */ ) { G4VBiasingOperation * GammaBiasingOperator::ProposeFinalStateBiasingOperation(const G4Track *track, const G4BiasingProcessInterface * /* callingProcess */) { + // -- Check if biasing of primary particle only is requested. If so, and // -- if particle is not a primary one, don't ask for biasing: if (fBiasPrimaryOnly && (track->GetParentID() != 0)) return nullptr; diff --git a/src/PhysicsList.cxx b/src/PhysicsList.cxx index 55748595..592e3801 100644 --- a/src/PhysicsList.cxx +++ b/src/PhysicsList.cxx @@ -3,6 +3,7 @@ #include +#include #include #include #include @@ -42,8 +43,8 @@ using namespace std; -PhysicsList::PhysicsList(SimulationManager* simulationManager, TRestGeant4PhysicsLists* physicsLists) - : G4VModularPhysicsList(), fSimulationManager(simulationManager) { +PhysicsList::PhysicsList(SimulationManager *simulationManager, TRestGeant4PhysicsLists *physicsLists) + : G4VModularPhysicsList(), fSimulationManager(simulationManager) { // add new units for radioActive decays const G4double G4minute = 60 * second; const G4double G4hour = 60 * G4minute; @@ -60,8 +61,8 @@ PhysicsList::PhysicsList(SimulationManager* simulationManager, TRestGeant4Physic G4LossTableManager::Instance(); // fix lower limit for cut G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange( - fRestPhysicsLists->GetMinimumEnergyProductionCuts() * keV, - fRestPhysicsLists->GetMaximumEnergyProductionCuts() * keV); + fRestPhysicsLists->GetMinimumEnergyProductionCuts() * keV, + fRestPhysicsLists->GetMaximumEnergyProductionCuts() * keV); InitializePhysicsLists(); } @@ -70,7 +71,7 @@ PhysicsList::~PhysicsList() { delete fEmPhysicsList; delete fDecPhysicsList; delete fRadDecPhysicsList; - for (auto& hadronicPhysicsList : fHadronPhys) { + for (auto &hadronicPhysicsList: fHadronPhys) { delete hadronicPhysicsList; } } @@ -128,6 +129,11 @@ void PhysicsList::InitializePhysicsLists() { emCounter++; } + fBiasingPhysicsList = new G4GenericBiasingPhysics(); + std::vector processToBias = {"eBrem"}; + fBiasingPhysicsList->PhysicsBias("e-", processToBias); + fBiasingPhysicsList->PhysicsBias("e+", processToBias); + if (fRestPhysicsLists->GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Essential && emCounter == 0) { RESTWarning << "PhysicsList: No EM physics list has been enabled" << RESTendl; @@ -179,7 +185,11 @@ void PhysicsList::ConstructParticle() { fRadDecPhysicsList->ConstructParticle(); } - for (auto& hadronicPhysicsList : fHadronPhys) { + if (fBiasingPhysicsList) { + fBiasingPhysicsList->ConstructParticle(); + } + + for (auto &hadronicPhysicsList: fHadronPhys) { hadronicPhysicsList->ConstructParticle(); } } @@ -192,27 +202,27 @@ void PhysicsList::ConstructProcess() { fEmPhysicsList->ConstructProcess(); fEmConfig.AddModels(); - G4UImanager* UI = G4UImanager::GetUIpointer(); + G4UImanager *UI = G4UImanager::GetUIpointer(); UI->ApplyCommand("/process/em/fluo true"); UI->ApplyCommand("/process/em/auger true"); UI->ApplyCommand("/process/em/pixe true"); bool boolEmOptionPixe = StringToBool( - fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "pixe", "false").Data()); + fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "pixe", "false").Data()); string stringEmOptionPixe = (boolEmOptionPixe ? "true" : "false"); G4cout << "Setting EM option '/process/em/pixe' to '" << stringEmOptionPixe << "' for physics list '" << fEmPhysicsListName << "'" << endl; UI->ApplyCommand(string("/process/em/pixe ") + stringEmOptionPixe); bool boolEmOptionFluo = StringToBool( - fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "fluo", "true").Data()); + fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "fluo", "true").Data()); string stringEmOptionFluo = (boolEmOptionFluo ? "true" : "false"); G4cout << "Setting EM option '/process/em/fluo' to '" << stringEmOptionFluo << "' for physics list '" << fEmPhysicsListName << "'" << endl; UI->ApplyCommand(string("/process/em/fluo ") + stringEmOptionFluo); bool boolEmOptionAuger = StringToBool( - fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "auger", "true").Data()); + fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "auger", "true").Data()); string stringEmOptionAuger = (boolEmOptionAuger ? "true" : "false"); G4cout << "Setting EM option '/process/em/auger' to '" << stringEmOptionAuger << "' for physics list '" << fEmPhysicsListName << "'" << endl; @@ -229,8 +239,12 @@ void PhysicsList::ConstructProcess() { fRadDecPhysicsList->ConstructProcess(); } + if (fBiasingPhysicsList) { + fBiasingPhysicsList->ConstructProcess(); + } + // Hadronic physics lists - for (auto& hadronicPhysicsList : fHadronPhys) { + for (auto &hadronicPhysicsList: fHadronPhys) { hadronicPhysicsList->ConstructProcess(); } @@ -282,15 +296,15 @@ void PhysicsList::ConstructProcess() { if (!(fRestPhysicsLists->GetPhysicsListOptionValue("G4RadioactiveDecay", "TritiumDecay") == "false") && (fRestPhysicsLists->GetPhysicsListOptionValue("G4RadioactiveDecay", "TritiumDecay", "false") == - "true" || + "true" || fSimulationManager->GetRestMetadata()->GetParticleSource()->GetParticleName() == "H3")) { // Tritium (H3) fix | https://geant4-forum.web.cern.ch/t/triton-decay-with-rdecay01-example/2616 - G4ParticleDefinition* tritium = G4Triton::Definition(); + G4ParticleDefinition *tritium = G4Triton::Definition(); tritium->SetPDGStable(false); - G4VProcess* decay = nullptr; - G4ProcessManager* tritiumProcessManager = tritium->GetProcessManager(); - G4ProcessVector* tritiumProcessVector = tritiumProcessManager->GetAtRestProcessVector(); + G4VProcess *decay = nullptr; + G4ProcessManager *tritiumProcessManager = tritium->GetProcessManager(); + G4ProcessVector *tritiumProcessVector = tritiumProcessManager->GetAtRestProcessVector(); for (unsigned int i = 0; i < tritiumProcessVector->size() && decay == nullptr; i++) { if ((*tritiumProcessVector)[i]->GetProcessName() == "Decay") decay = (*tritiumProcessVector)[i]; @@ -311,9 +325,9 @@ void PhysicsList::ConstructProcess() { // To implement UserLimits to StepSize inside the gas theParticleIterator->reset(); while ((*theParticleIterator)()) { - G4ParticleDefinition* particle = theParticleIterator->value(); - const auto& particleName = particle->GetParticleName(); - G4ProcessManager* processManager = particle->GetProcessManager(); + G4ParticleDefinition *particle = theParticleIterator->value(); + const auto &particleName = particle->GetParticleName(); + G4ProcessManager *processManager = particle->GetProcessManager(); if (particleName == "e-") { processManager->AddDiscreteProcess(new G4StepLimiter("e-Step")); @@ -333,10 +347,10 @@ void PhysicsList::ConstructProcess() { for (int A = 2 * Z; A <= 3 * Z; A++) { for (unsigned int n = 0; n < fRestPhysicsLists->GetIonStepList().size(); n++) { if (fRestPhysicsLists->GetIonStepList()[n] == G4IonTable::GetIonTable()->GetIonName(Z, A)) { - G4ParticleDefinition* particle = G4IonTable::GetIonTable()->GetIon(Z, A, 0); + G4ParticleDefinition *particle = G4IonTable::GetIonTable()->GetIon(Z, A, 0); G4String particle_name = G4IonTable::GetIonTable()->GetIonName(Z, A, 0); cout << "Found ion: " << particle_name << " Z " << Z << " A " << A << endl; - G4ProcessManager* processManager = particle->GetProcessManager(); + G4ProcessManager *processManager = particle->GetProcessManager(); processManager->AddDiscreteProcess(new G4StepLimiter("ionStep")); } } diff --git a/src/TrackingAction.cxx b/src/TrackingAction.cxx index 7799557b..095a7b9d 100644 --- a/src/TrackingAction.cxx +++ b/src/TrackingAction.cxx @@ -9,13 +9,14 @@ using namespace std; -TrackingAction::TrackingAction(SimulationManager* simulationManager) - : G4UserTrackingAction(), fSimulationManager(simulationManager) {} +TrackingAction::TrackingAction(SimulationManager *simulationManager) + : G4UserTrackingAction(), fSimulationManager(simulationManager) {} -void TrackingAction::PreUserTrackingAction(const G4Track* track) { +void TrackingAction::PreUserTrackingAction(const G4Track *track) { + // G4cout << track->GetVolume()->GetLogicalVolume()->GetName() << G4endl; fSimulationManager->GetOutputManager()->RecordTrack(track); } -void TrackingAction::PostUserTrackingAction(const G4Track* track) { +void TrackingAction::PostUserTrackingAction(const G4Track *track) { fSimulationManager->GetOutputManager()->UpdateTrack(track); } From df8707b70da84a523526c10e63740ce679329f96 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Wed, 22 Feb 2023 12:28:28 +0100 Subject: [PATCH 06/24] split only when pointing towards (0,0,0) --- src/GammaBiasingOperation.cxx | 28 +++++++++++++++++++++++++--- src/GammaBiasingOperator.cxx | 2 +- 2 files changed, 26 insertions(+), 4 deletions(-) diff --git a/src/GammaBiasingOperation.cxx b/src/GammaBiasingOperation.cxx index 893fe7ea..e56794bc 100644 --- a/src/GammaBiasingOperation.cxx +++ b/src/GammaBiasingOperation.cxx @@ -62,10 +62,8 @@ ApplyFinalStateBiasing(const G4BiasingProcessInterface *callingProcess, // -- Now deal with the gamma's: // -- their common weight: - G4double gammaWeight = track->GetWeight() / fSplittingFactor; // -- inform we will have fSplittingFactor gamma's: - fParticleChange.SetNumberOfSecondaries(fSplittingFactor); // -- inform we take care of secondaries weight (otherwise these // -- secondaries are by default given the primary weight). @@ -73,15 +71,39 @@ ApplyFinalStateBiasing(const G4BiasingProcessInterface *callingProcess, // -- Store first gamma: G4Track *gammaTrack = actualParticleChange->GetSecondary(0); + + + // print gamma info (energy, direction) + /* + G4cout << "Gamma info. Weight: " << gammaTrack->GetWeight() + << " Energy (keV): " << gammaTrack->GetKineticEnergy() / CLHEP::keV << " Direction: " + << gammaTrack->GetMomentumDirection() << G4endl; + */ + // if direction points towards (0,0,0) + bool split = false; + const auto diff = gammaTrack->GetPosition() - G4ThreeVector(0, 0, 0); + if (gammaTrack->GetMomentumDirection().dot(diff) < 0) { + // G4cout << "Gamma points towards (0,0,0)" << G4endl; + split = true; + } + + const int nSecondaries = split ? fSplittingFactor : 1; + + G4double gammaWeight = track->GetWeight() / nSecondaries; + // G4cout << "Gamma weight: " << gammaWeight << " nSecondaries: " << nSecondaries << G4endl; gammaTrack->SetWeight(gammaWeight); + + + fParticleChange.SetNumberOfSecondaries(nSecondaries); fParticleChange.AddSecondary(gammaTrack); // -- and clean-up the brem. process particle change: + actualParticleChange->Clear(); // -- now start the fSplittingFactor-1 calls to the brem. process to store each // -- related gamma: G4int nCalls = 1; - while (nCalls < fSplittingFactor) { + while (nCalls < nSecondaries) { // ( note: we don't need to cast to actual type here, as methods for accessing // secondary particles are from base class G4VParticleChange ) processFinalState = callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); diff --git a/src/GammaBiasingOperator.cxx b/src/GammaBiasingOperator.cxx index 9f780a2e..4cae07e0 100644 --- a/src/GammaBiasingOperator.cxx +++ b/src/GammaBiasingOperator.cxx @@ -10,7 +10,7 @@ GammaBiasingOperator::GammaBiasingOperator() : G4VBiasingOperator("BremSplittingOperator"), - fSplittingFactor(10), + fSplittingFactor(2), fBiasPrimaryOnly(false), fBiasOnlyOnce(true) { fBremSplittingOperation = new GammaBiasingOperation("BremSplittingOperation"); From 04522541354e5c74ccdb0dd63db80779462ba5ec Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Wed, 22 Feb 2023 13:26:03 +0100 Subject: [PATCH 07/24] attempt to integrate with TRestGeant4Metadata --- examples/13.IAXO/GammaBiasing.rml | 6 ++++ include/GammaBiasingOperator.h | 1 - src/DetectorConstruction.cxx | 60 ++++++++++++++++++++++++------- src/GammaBiasingOperation.cxx | 44 +++-------------------- src/GammaBiasingOperator.cxx | 28 +++------------ 5 files changed, 63 insertions(+), 76 deletions(-) diff --git a/examples/13.IAXO/GammaBiasing.rml b/examples/13.IAXO/GammaBiasing.rml index d6f2faa9..0a6478b4 100644 --- a/examples/13.IAXO/GammaBiasing.rml +++ b/examples/13.IAXO/GammaBiasing.rml @@ -22,6 +22,12 @@ Author: Luis Obis (lobis@unizar.es) + + + + + + diff --git a/include/GammaBiasingOperator.h b/include/GammaBiasingOperator.h index c0ce5253..5f6729e9 100644 --- a/include/GammaBiasingOperator.h +++ b/include/GammaBiasingOperator.h @@ -39,7 +39,6 @@ class GammaBiasingOperator : public G4VBiasingOperator { private: GammaBiasingOperation *fBremSplittingOperation; G4int fSplittingFactor; - G4bool fBiasPrimaryOnly; G4bool fBiasOnlyOnce; G4int fNInteractions; }; diff --git a/src/DetectorConstruction.cxx b/src/DetectorConstruction.cxx index 193550f9..a32f83e0 100644 --- a/src/DetectorConstruction.cxx +++ b/src/DetectorConstruction.cxx @@ -296,20 +296,54 @@ void DetectorConstruction::ConstructSDandField() { } // Biasing - const auto biasingLogicalVolumeName = "ShieldingVolume"; - G4LogicalVolume *biasingLogicalVolume = - G4LogicalVolumeStore::GetInstance()->GetVolume(biasingLogicalVolumeName, false); - if (biasingLogicalVolume == nullptr) { - cerr << "Detector construction error: could not find logical volume '" << biasingLogicalVolumeName << "'" - << endl; - exit(1); - } + const auto &biasingInfo = metadata.GetGeant4BiasingInfo(); + if (biasingInfo.GetEnabled()) { + const auto &geometryInfo = metadata.GetGeant4GeometryInfo(); + + set biasingVolumes; + for (const auto &biasingVolume: biasingInfo.GetBiasingVolumes()) { + auto name = biasingVolume; + if (!geometryInfo.IsValidLogicalVolume(name)) { + if (geometryInfo.IsValidPhysicalVolume(name)) { + name = geometryInfo.GetLogicalVolumeFromPhysical(name); + biasingVolumes.insert(name); + } else { + RESTError << "volume name '" << name + << "' inside biasing section is invalid. Please check it belongs to a logical volume" + << RESTendl; + exit(1); + } + } + } - GammaBiasingOperator *bremSplittingOperator = new GammaBiasingOperator(); - bremSplittingOperator->AttachTo(biasingLogicalVolume); - G4cout << " Attaching biasing operator " << bremSplittingOperator->GetName() - << " to logical volume " << biasingLogicalVolume->GetName() - << G4endl; + for (const auto &biasingLogicalVolumeName: biasingVolumes) { + G4LogicalVolume *biasingLogicalVolume = + G4LogicalVolumeStore::GetInstance()->GetVolume(biasingLogicalVolumeName, false); + if (biasingLogicalVolume == nullptr) { + cerr << "Detector construction error: could not find logical volume '" << biasingLogicalVolumeName + << "'" + << endl; + exit(1); + } + + // TODO: better memory management + auto gammaBiasing = new GammaBiasingOperator(); + gammaBiasing->AttachTo(biasingLogicalVolume); + G4cout << " Attaching biasing operator " << gammaBiasing->GetName() + << " to logical volume " << biasingLogicalVolume->GetName() + << G4endl; + } + + // Print info + cout << "Biasing is enabled " << endl; + cout << "Splitting factor: " << biasingInfo.GetSplittingFactor() << endl; + cout << "Biasing volumes: " << endl; + for (const auto &volume: biasingVolumes) { + cout << "\t" << volume << endl; + } + const auto center = biasingInfo.GetBiasingCenter(); + cout << "Biasing center: " << center.x() << ", " << center.y() << ", " << center.z() << endl; + } } void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume *world) { diff --git a/src/GammaBiasingOperation.cxx b/src/GammaBiasingOperation.cxx index e56794bc..9f5d8f6c 100644 --- a/src/GammaBiasingOperation.cxx +++ b/src/GammaBiasingOperation.cxx @@ -15,39 +15,19 @@ GammaBiasingOperation::GammaBiasingOperation(G4String name) GammaBiasingOperation::~GammaBiasingOperation() = default; -G4VParticleChange * -GammaBiasingOperation:: -ApplyFinalStateBiasing(const G4BiasingProcessInterface *callingProcess, - const G4Track *track, - const G4Step *step, - G4bool &) { +G4VParticleChange *GammaBiasingOperation::ApplyFinalStateBiasing(const G4BiasingProcessInterface *callingProcess, + const G4Track *track, + const G4Step *step, + G4bool &) { G4VParticleChange *processFinalState = callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); if (fSplittingFactor == 1) return processFinalState; - // -- a special case here: the brem. process corrects for cross-section change - // -- over the step due to energy loss by sometimes "abandoning" the interaction, - // -- returning an unchanged incoming electron/positron. - // -- We respect this correction, and if no secondary is produced, its means this - // -- case is happening: + // special case: no secondaries if (processFinalState->GetNumberOfSecondaries() == 0) return processFinalState; - // -- Now start the biasing: - // -- - the electron state will be taken as the first one produced by the brem. - // -- process, hence the one stored in above processFinalState particle change. - // -- This state will be stored in our fParticleChange object. - // -- - the photon accompanying the electron will be stored also this way. - // -- - we will then do fSplittingFactor - 1 call to the brem. process to collect - // -- fSplittingFactor - 1 additional gammas. All these will be stored in our - // -- fParticleChange object. - - // -- We called the brem. process above. Its concrete particle change is indeed - // -- a "G4ParticleChangeForLoss" object. We cast this particle change to access - // -- methods of the concrete G4ParticleChangeForLoss type: - - auto actualParticleChange = (G4ParticleChangeForLoss *) processFinalState; fParticleChange.Initialize(*track); @@ -60,16 +40,8 @@ ApplyFinalStateBiasing(const G4BiasingProcessInterface *callingProcess, fParticleChange. ProposeMomentumDirection(actualParticleChange->GetProposedMomentumDirection()); - // -- Now deal with the gamma's: - // -- their common weight: - - // -- inform we will have fSplittingFactor gamma's: - - // -- inform we take care of secondaries weight (otherwise these - // -- secondaries are by default given the primary weight). fParticleChange.SetSecondaryWeightByProcess(true); - // -- Store first gamma: G4Track *gammaTrack = actualParticleChange->GetSecondary(0); @@ -96,16 +68,11 @@ ApplyFinalStateBiasing(const G4BiasingProcessInterface *callingProcess, fParticleChange.SetNumberOfSecondaries(nSecondaries); fParticleChange.AddSecondary(gammaTrack); - // -- and clean-up the brem. process particle change: actualParticleChange->Clear(); - // -- now start the fSplittingFactor-1 calls to the brem. process to store each - // -- related gamma: G4int nCalls = 1; while (nCalls < nSecondaries) { - // ( note: we don't need to cast to actual type here, as methods for accessing - // secondary particles are from base class G4VParticleChange ) processFinalState = callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); if (processFinalState->GetNumberOfSecondaries() == 1) { gammaTrack = processFinalState->GetSecondary(0); @@ -121,6 +88,5 @@ ApplyFinalStateBiasing(const G4BiasingProcessInterface *callingProcess, processFinalState->Clear(); } - // -- we are done: return &fParticleChange; } diff --git a/src/GammaBiasingOperator.cxx b/src/GammaBiasingOperator.cxx index 4cae07e0..92bc3330 100644 --- a/src/GammaBiasingOperator.cxx +++ b/src/GammaBiasingOperator.cxx @@ -10,8 +10,7 @@ GammaBiasingOperator::GammaBiasingOperator() : G4VBiasingOperator("BremSplittingOperator"), - fSplittingFactor(2), - fBiasPrimaryOnly(false), + fSplittingFactor(1), fBiasOnlyOnce(true) { fBremSplittingOperation = new GammaBiasingOperation("BremSplittingOperation"); } @@ -19,35 +18,18 @@ GammaBiasingOperator::GammaBiasingOperator() void GammaBiasingOperator::StartRun() { fBremSplittingOperation->SetSplittingFactor(fSplittingFactor); - G4cout << GetName() << " : starting run with brem. splitting factor = " - << fSplittingFactor; - if (fBiasPrimaryOnly) G4cout << ", biasing only primaries "; - else - G4cout << ", biasing primary and secondary tracks "; - if (fBiasOnlyOnce) G4cout << ", biasing only once per track "; - else - G4cout << ", biasing several times per track "; - G4cout << " . " << G4endl; } -void GammaBiasingOperator::StartTracking(const G4Track * /* track */ ) { - // -- reset the number of times the brem. splitting was applied: +void GammaBiasingOperator::StartTracking(const G4Track *) { fNInteractions = 0; } G4VBiasingOperation * -GammaBiasingOperator::ProposeFinalStateBiasingOperation(const G4Track *track, - const G4BiasingProcessInterface * /* callingProcess */) { - - // -- Check if biasing of primary particle only is requested. If so, and - // -- if particle is not a primary one, don't ask for biasing: - if (fBiasPrimaryOnly && (track->GetParentID() != 0)) return nullptr; - // -- Check if brem. splitting should be applied only once to the track, - // -- and if so, and if brem. splitting already occured, don't ask for biasing: +GammaBiasingOperator::ProposeFinalStateBiasingOperation(const G4Track *, + const G4BiasingProcessInterface *) { if (fBiasOnlyOnce && (fNInteractions > 0)) return nullptr; - // -- Count the number of times the brem. splitting is applied: fNInteractions++; - // -- Return the brem. splitting operation: + return fBremSplittingOperation; } \ No newline at end of file From 520aed1d33d978b46e488a63bfc6f19cae3c33e7 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Wed, 22 Feb 2023 12:26:36 +0000 Subject: [PATCH 08/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- include/GammaBiasingOperation.h | 36 ++++---- include/GammaBiasingOperator.h | 36 ++++---- include/PhysicsList.h | 24 +++--- src/DetectorConstruction.cxx | 145 ++++++++++++++++---------------- src/GammaBiasingOperation.cxx | 39 +++------ src/GammaBiasingOperator.cxx | 22 ++--- src/PhysicsList.cxx | 44 +++++----- src/TrackingAction.cxx | 8 +- 8 files changed, 162 insertions(+), 192 deletions(-) diff --git a/include/GammaBiasingOperation.h b/include/GammaBiasingOperation.h index 3579a03d..75dddba7 100644 --- a/include/GammaBiasingOperation.h +++ b/include/GammaBiasingOperation.h @@ -5,43 +5,37 @@ #ifndef REST_GAMMABIASINGOPERATION_H #define REST_GAMMABIASINGOPERATION_H -#include "G4VBiasingOperation.hh" -#include "G4ParticleChange.hh" #include "G4BiasingProcessInterface.hh" +#include "G4ParticleChange.hh" +#include "G4VBiasingOperation.hh" class GammaBiasingOperation : public G4VBiasingOperation { -public: + public: GammaBiasingOperation(G4String name); virtual ~GammaBiasingOperation(); -public: - virtual const G4VBiasingInteractionLaw * - ProvideOccurenceBiasingInteractionLaw(const G4BiasingProcessInterface *, - G4ForceCondition &) { return 0; } - - virtual G4VParticleChange *ApplyFinalStateBiasing(const G4BiasingProcessInterface *, - const G4Track *, - const G4Step *, - G4bool &); + public: + virtual const G4VBiasingInteractionLaw* ProvideOccurenceBiasingInteractionLaw( + const G4BiasingProcessInterface*, G4ForceCondition&) { + return 0; + } - virtual G4double DistanceToApplyOperation(const G4Track *, - G4double, - G4ForceCondition *) { return DBL_MAX; } + virtual G4VParticleChange* ApplyFinalStateBiasing(const G4BiasingProcessInterface*, const G4Track*, + const G4Step*, G4bool&); - virtual G4VParticleChange *GenerateBiasingFinalState(const G4Track *, - const G4Step *) { return 0; } + virtual G4double DistanceToApplyOperation(const G4Track*, G4double, G4ForceCondition*) { return DBL_MAX; } + virtual G4VParticleChange* GenerateBiasingFinalState(const G4Track*, const G4Step*) { return 0; } -public: + public: void SetSplittingFactor(G4int splittingFactor) { fSplittingFactor = splittingFactor; } G4int GetSplittingFactor() const { return fSplittingFactor; } -private: + private: G4int fSplittingFactor; G4ParticleChange fParticleChange; }; - -#endif //REST_GAMMABIASINGOPERATION_H +#endif // REST_GAMMABIASINGOPERATION_H diff --git a/include/GammaBiasingOperator.h b/include/GammaBiasingOperator.h index 5f6729e9..7a845db2 100644 --- a/include/GammaBiasingOperator.h +++ b/include/GammaBiasingOperator.h @@ -10,38 +10,38 @@ class GammaBiasingOperation; class GammaBiasingOperator : public G4VBiasingOperator { -public: + public: GammaBiasingOperator(); virtual ~GammaBiasingOperator() {} -public: + public: virtual void StartRun(); - virtual void StartTracking(const G4Track *track); + virtual void StartTracking(const G4Track* track); -private: - virtual G4VBiasingOperation * - ProposeNonPhysicsBiasingOperation(const G4Track *, - const G4BiasingProcessInterface *) { return 0; } + private: + virtual G4VBiasingOperation* ProposeNonPhysicsBiasingOperation(const G4Track*, + const G4BiasingProcessInterface*) { + return 0; + } - virtual G4VBiasingOperation * - ProposeOccurenceBiasingOperation(const G4Track *, - const G4BiasingProcessInterface *) { return 0; } + virtual G4VBiasingOperation* ProposeOccurenceBiasingOperation(const G4Track*, + const G4BiasingProcessInterface*) { + return 0; + } - virtual G4VBiasingOperation * - ProposeFinalStateBiasingOperation(const G4Track *track, - const G4BiasingProcessInterface *callingProcess); + virtual G4VBiasingOperation* ProposeFinalStateBiasingOperation( + const G4Track* track, const G4BiasingProcessInterface* callingProcess); -private: + private: using G4VBiasingOperator::OperationApplied; -private: - GammaBiasingOperation *fBremSplittingOperation; + private: + GammaBiasingOperation* fBremSplittingOperation; G4int fSplittingFactor; G4bool fBiasOnlyOnce; G4int fNInteractions; }; - -#endif //REST_GAMMABIASINGOPERATOR_H +#endif // REST_GAMMABIASINGOPERATOR_H diff --git a/include/PhysicsList.h b/include/PhysicsList.h index 1f56feef..0c275159 100644 --- a/include/PhysicsList.h +++ b/include/PhysicsList.h @@ -10,18 +10,17 @@ #include "SimulationManager.h" - class G4GenericBiasingPhysics; class PhysicsList : public G4VModularPhysicsList { -public: + public: PhysicsList() = delete; - explicit PhysicsList(SimulationManager *simulationManager, TRestGeant4PhysicsLists *restPhysicsLists); + explicit PhysicsList(SimulationManager* simulationManager, TRestGeant4PhysicsLists* restPhysicsLists); ~PhysicsList() override; -protected: + protected: // Construct particle and physics virtual void InitializePhysicsLists(); @@ -31,22 +30,21 @@ class PhysicsList : public G4VModularPhysicsList { void SetCuts() override; -private: - SimulationManager *fSimulationManager = nullptr; + private: + SimulationManager* fSimulationManager = nullptr; G4EmConfigurator fEmConfig; - G4VPhysicsConstructor *fEmPhysicsList = nullptr; + G4VPhysicsConstructor* fEmPhysicsList = nullptr; std::string fEmPhysicsListName; // Can be different from the output of GetPhysicsName - G4VPhysicsConstructor *fDecPhysicsList = nullptr; - G4VPhysicsConstructor *fRadDecPhysicsList = nullptr; - std::vector fHadronPhys; - - G4GenericBiasingPhysics *fBiasingPhysicsList = nullptr; + G4VPhysicsConstructor* fDecPhysicsList = nullptr; + G4VPhysicsConstructor* fRadDecPhysicsList = nullptr; + std::vector fHadronPhys; - TRestGeant4PhysicsLists *fRestPhysicsLists = nullptr; + G4GenericBiasingPhysics* fBiasingPhysicsList = nullptr; + TRestGeant4PhysicsLists* fRestPhysicsLists = nullptr; }; #endif diff --git a/src/DetectorConstruction.cxx b/src/DetectorConstruction.cxx index a32f83e0..3dc6ad49 100644 --- a/src/DetectorConstruction.cxx +++ b/src/DetectorConstruction.cxx @@ -17,22 +17,22 @@ #include #include -#include "SimulationManager.h" #include "GammaBiasingOperator.h" +#include "SimulationManager.h" using namespace std; using namespace TRestGeant4PrimaryGeneratorTypes; -DetectorConstruction::DetectorConstruction(SimulationManager *simulationManager) - : fSimulationManager(simulationManager) { +DetectorConstruction::DetectorConstruction(SimulationManager* simulationManager) + : fSimulationManager(simulationManager) { G4cout << "Detector Construction" << G4endl; fGdmlParser = new G4GDMLParser(); } DetectorConstruction::~DetectorConstruction() { delete fGdmlParser; } -G4VPhysicalVolume *DetectorConstruction::Construct() { - TRestGeant4Metadata *restG4Metadata = fSimulationManager->GetRestMetadata(); +G4VPhysicalVolume* DetectorConstruction::Construct() { + TRestGeant4Metadata* restG4Metadata = fSimulationManager->GetRestMetadata(); cout << "Isotope table " << endl; cout << *(G4Isotope::GetIsotopeTable()) << endl; @@ -45,31 +45,31 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { const auto startingPath = filesystem::current_path(); const auto [gdmlPath, gdmlToRead] = - TRestTools::SeparatePathAndName((string) restG4Metadata->GetGdmlFilename()); + TRestTools::SeparatePathAndName((string)restG4Metadata->GetGdmlFilename()); filesystem::current_path(gdmlPath); G4cout << "gdmlToRead: " << gdmlToRead << G4endl; fGdmlParser->Read(gdmlToRead, false); - G4VPhysicalVolume *worldVolume = fGdmlParser->GetWorldVolume(); + G4VPhysicalVolume* worldVolume = fGdmlParser->GetWorldVolume(); - const auto worldSolid = dynamic_cast(worldVolume->GetLogicalVolume()->GetSolid()); + const auto worldSolid = dynamic_cast(worldVolume->GetLogicalVolume()->GetSolid()); restG4Metadata->fGeant4PrimaryGeneratorInfo.fSpatialGeneratorWorldSize = { - worldSolid->GetXHalfLength(), worldSolid->GetYHalfLength(), worldSolid->GetZHalfLength()}; + worldSolid->GetXHalfLength(), worldSolid->GetYHalfLength(), worldSolid->GetZHalfLength()}; restG4Metadata->fGeant4GeometryInfo.InitializeOnDetectorConstruction(gdmlToRead, worldVolume); restG4Metadata->ReadDetector(); restG4Metadata->PrintMetadata(); // now we have detector info - const auto &geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); + const auto& geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); geometryInfo.Print(); // do some checks { // Check all physical volume names are unique - G4PhysicalVolumeStore *physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); + G4PhysicalVolumeStore* physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); set physicalVolumeNames; - vector::const_iterator physicalVolume; + vector::const_iterator physicalVolume; for (physicalVolume = physicalVolumeStore->begin(); physicalVolume != physicalVolumeStore->end(); physicalVolume++) { auto name = (*physicalVolume)->GetName(); @@ -84,9 +84,9 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { } // Check all logical volume names are unique - G4LogicalVolumeStore *logicalVolumeStore = G4LogicalVolumeStore::GetInstance(); + G4LogicalVolumeStore* logicalVolumeStore = G4LogicalVolumeStore::GetInstance(); set logicalVolumeNames; - vector::const_iterator logicalVolume; + vector::const_iterator logicalVolume; for (logicalVolume = logicalVolumeStore->begin(); logicalVolume != logicalVolumeStore->end(); logicalVolume++) { auto name = (*logicalVolume)->GetName(); @@ -102,15 +102,15 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { } filesystem::current_path(startingPath); - auto sensitiveVolume = (string) restG4Metadata->GetSensitiveVolume(); - G4VPhysicalVolume *physicalVolume = GetPhysicalVolume(sensitiveVolume); + auto sensitiveVolume = (string)restG4Metadata->GetSensitiveVolume(); + G4VPhysicalVolume* physicalVolume = GetPhysicalVolume(sensitiveVolume); if (physicalVolume == nullptr) { // sensitive volume was not found, perhaps the user specified a logical volume auto physicalVolumes = geometryInfo.GetAllPhysicalVolumesFromLogical(sensitiveVolume); if (physicalVolumes.size() == 1) { restG4Metadata->InsertSensitiveVolume( - geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolumes[0])); - sensitiveVolume = (string) restG4Metadata->GetSensitiveVolume(); + geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolumes[0])); + sensitiveVolume = (string)restG4Metadata->GetSensitiveVolume(); physicalVolume = GetPhysicalVolume(sensitiveVolume); } } @@ -124,20 +124,20 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { Double_t my = restG4Metadata->GetMagneticField().Y() * tesla; Double_t mz = restG4Metadata->GetMagneticField().Z() * tesla; - G4MagneticField *magField = new G4UniformMagField(G4ThreeVector(mx, my, mz)); + G4MagneticField* magField = new G4UniformMagField(G4ThreeVector(mx, my, mz)); // G4FieldManager* localFieldMgr = new G4FieldManager(magField); - G4FieldManager *fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager(); + G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager(); fieldMgr->SetDetectorField(magField); fieldMgr->CreateChordFinder(magField); - G4LogicalVolume *volume = physicalVolume->GetLogicalVolume(); - G4Material *material = volume->GetMaterial(); + G4LogicalVolume* volume = physicalVolume->GetLogicalVolume(); + G4Material* material = volume->GetMaterial(); G4cout << "Sensitive volume properties:" << G4endl; G4cout << "\t- Material: " << material->GetName() << G4endl; G4cout << "\t- Temperature: " << material->GetTemperature() << " K" << G4endl; G4cout << "\t- Density: " << material->GetDensity() / (g / cm3) << " g/cm3" << G4endl; - const auto &primaryGeneratorInfo = restG4Metadata->GetGeant4PrimaryGeneratorInfo(); + const auto& primaryGeneratorInfo = restG4Metadata->GetGeant4PrimaryGeneratorInfo(); // Getting generation volume const auto fromVolume = primaryGeneratorInfo.GetSpatialGeneratorFrom(); if (fromVolume != "NO_SUCH_PARA") { @@ -146,11 +146,11 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { cout << "Generator type: " << primaryGeneratorInfo.GetSpatialGeneratorType() << endl; const auto spatialGeneratorTypeEnum = - StringToSpatialGeneratorTypes(primaryGeneratorInfo.GetSpatialGeneratorType().Data()); + StringToSpatialGeneratorTypes(primaryGeneratorInfo.GetSpatialGeneratorType().Data()); if (spatialGeneratorTypeEnum == TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::VOLUME && primaryGeneratorInfo.GetSpatialGeneratorFrom() != "Not defined") { - G4VPhysicalVolume *pVol = GetPhysicalVolume(primaryGeneratorInfo.GetSpatialGeneratorFrom().Data()); + G4VPhysicalVolume* pVol = GetPhysicalVolume(primaryGeneratorInfo.GetSpatialGeneratorFrom().Data()); if (pVol == nullptr) { cout << "ERROR: The generator volume '" << primaryGeneratorInfo.GetSpatialGeneratorFrom() << "'was not found in the geometry" << endl; @@ -162,7 +162,7 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { if (spatialGeneratorTypeEnum == TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::SURFACE || spatialGeneratorTypeEnum == TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::VOLUME) { restG4Metadata->fGeant4PrimaryGeneratorInfo.fSpatialGeneratorPosition = { - fGeneratorTranslation.x(), fGeneratorTranslation.y(), fGeneratorTranslation.z()}; + fGeneratorTranslation.x(), fGeneratorTranslation.y(), fGeneratorTranslation.z()}; } fGeneratorSolid = pVol->GetLogicalVolume()->GetSolid(); @@ -203,9 +203,9 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { for (unsigned int id = 0; id < restG4Metadata->GetNumberOfActiveVolumes(); id++) { TString activeVolumeName = restG4Metadata->GetActiveVolumeName(id); - G4VPhysicalVolume *pVol = GetPhysicalVolume((G4String) activeVolumeName); + G4VPhysicalVolume* pVol = GetPhysicalVolume((G4String)activeVolumeName); if (pVol != nullptr) { - G4LogicalVolume *lVol = pVol->GetLogicalVolume(); + G4LogicalVolume* lVol = pVol->GetLogicalVolume(); if (restG4Metadata->GetMaxStepSize(activeVolumeName) > 0) { RESTInfo << "Setting maxStepSize of " << restG4Metadata->GetMaxStepSize(activeVolumeName) * mm << "mm for volume '" << activeVolumeName << "'" << RESTendl; @@ -221,15 +221,15 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { return worldVolume; } -G4VPhysicalVolume *DetectorConstruction::GetPhysicalVolume(const G4String &physicalVolumeName) const { - G4PhysicalVolumeStore *physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); - TRestGeant4Metadata *restG4Metadata = fSimulationManager->GetRestMetadata(); - const auto &geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); - vector::const_iterator physicalVolume; +G4VPhysicalVolume* DetectorConstruction::GetPhysicalVolume(const G4String& physicalVolumeName) const { + G4PhysicalVolumeStore* physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); + TRestGeant4Metadata* restG4Metadata = fSimulationManager->GetRestMetadata(); + const auto& geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); + vector::const_iterator physicalVolume; for (physicalVolume = physicalVolumeStore->begin(); physicalVolume != physicalVolumeStore->end(); physicalVolume++) { auto name = (*physicalVolume)->GetName(); - auto alternativeName = (G4String) geometryInfo.GetAlternativeNameFromGeant4PhysicalName(name); + auto alternativeName = (G4String)geometryInfo.GetAlternativeNameFromGeant4PhysicalName(name); if (name == physicalVolumeName || alternativeName == physicalVolumeName) { return *physicalVolume; } @@ -239,20 +239,20 @@ G4VPhysicalVolume *DetectorConstruction::GetPhysicalVolume(const G4String &physi } void DetectorConstruction::ConstructSDandField() { - const TRestGeant4Metadata &metadata = *fSimulationManager->GetRestMetadata(); + const TRestGeant4Metadata& metadata = *fSimulationManager->GetRestMetadata(); - set logicalVolumesSelected; - for (const auto &userSensitiveVolume: metadata.GetSensitiveVolumes()) { + set logicalVolumesSelected; + for (const auto& userSensitiveVolume : metadata.GetSensitiveVolumes()) { // Each sensitive detector declaration in the RML should correspond to at least one logical volume - G4LogicalVolume *logicalVolume; + G4LogicalVolume* logicalVolume; // Check if user selected a Geant4 physical volume by name - G4VPhysicalVolume *physicalVolume = - G4PhysicalVolumeStore::GetInstance()->GetVolume(userSensitiveVolume.Data(), false); + G4VPhysicalVolume* physicalVolume = + G4PhysicalVolumeStore::GetInstance()->GetVolume(userSensitiveVolume.Data(), false); if (physicalVolume == nullptr) { const G4String geant4VolumeName = - metadata.GetGeant4GeometryInfo() - .GetGeant4PhysicalNameFromAlternativeName(userSensitiveVolume.Data()) - .Data(); + metadata.GetGeant4GeometryInfo() + .GetGeant4PhysicalNameFromAlternativeName(userSensitiveVolume.Data()) + .Data(); // Check if user selected a Geant4 physical volume by REST name (only relevant for assemblies) physicalVolume = G4PhysicalVolumeStore::GetInstance()->GetVolume(geant4VolumeName, false); } @@ -269,25 +269,25 @@ void DetectorConstruction::ConstructSDandField() { // Check if the user string matches a logical volume by expanding the input as a regex // This can have multiple hits const auto logicalVolumesMatchingExpression = - metadata.GetGeant4GeometryInfo().GetAllLogicalVolumesMatchingExpression(userSensitiveVolume); + metadata.GetGeant4GeometryInfo().GetAllLogicalVolumesMatchingExpression(userSensitiveVolume); if (logicalVolumesMatchingExpression.empty()) { cerr << "Detector construction error: could not find matching logical volume(s) for '" << userSensitiveVolume << "'" << endl; exit(1); } else { - for (const auto &logicalVolumeName: logicalVolumesMatchingExpression) { + for (const auto& logicalVolumeName : logicalVolumesMatchingExpression) { logicalVolumesSelected.insert( - G4LogicalVolumeStore::GetInstance()->GetVolume(logicalVolumeName.Data(), false)); + G4LogicalVolumeStore::GetInstance()->GetVolume(logicalVolumeName.Data(), false)); } } } } - G4SDManager *SDManager = G4SDManager::GetSDMpointer(); + G4SDManager* SDManager = G4SDManager::GetSDMpointer(); - for (G4LogicalVolume *logicalVolume: logicalVolumesSelected) { + for (G4LogicalVolume* logicalVolume : logicalVolumesSelected) { auto name = logicalVolume->GetName(); - G4VSensitiveDetector *sensitiveDetector = new SensitiveDetector(fSimulationManager, name); + G4VSensitiveDetector* sensitiveDetector = new SensitiveDetector(fSimulationManager, name); SDManager->AddNewDetector(sensitiveDetector); logicalVolume->SetSensitiveDetector(sensitiveDetector); @@ -296,49 +296,48 @@ void DetectorConstruction::ConstructSDandField() { } // Biasing - const auto &biasingInfo = metadata.GetGeant4BiasingInfo(); + const auto& biasingInfo = metadata.GetGeant4BiasingInfo(); if (biasingInfo.GetEnabled()) { - const auto &geometryInfo = metadata.GetGeant4GeometryInfo(); + const auto& geometryInfo = metadata.GetGeant4GeometryInfo(); set biasingVolumes; - for (const auto &biasingVolume: biasingInfo.GetBiasingVolumes()) { + for (const auto& biasingVolume : biasingInfo.GetBiasingVolumes()) { auto name = biasingVolume; if (!geometryInfo.IsValidLogicalVolume(name)) { if (geometryInfo.IsValidPhysicalVolume(name)) { name = geometryInfo.GetLogicalVolumeFromPhysical(name); biasingVolumes.insert(name); } else { - RESTError << "volume name '" << name - << "' inside biasing section is invalid. Please check it belongs to a logical volume" - << RESTendl; + RESTError + << "volume name '" << name + << "' inside biasing section is invalid. Please check it belongs to a logical volume" + << RESTendl; exit(1); } } } - for (const auto &biasingLogicalVolumeName: biasingVolumes) { - G4LogicalVolume *biasingLogicalVolume = - G4LogicalVolumeStore::GetInstance()->GetVolume(biasingLogicalVolumeName, false); + for (const auto& biasingLogicalVolumeName : biasingVolumes) { + G4LogicalVolume* biasingLogicalVolume = + G4LogicalVolumeStore::GetInstance()->GetVolume(biasingLogicalVolumeName, false); if (biasingLogicalVolume == nullptr) { - cerr << "Detector construction error: could not find logical volume '" << biasingLogicalVolumeName - << "'" - << endl; + cerr << "Detector construction error: could not find logical volume '" + << biasingLogicalVolumeName << "'" << endl; exit(1); } // TODO: better memory management auto gammaBiasing = new GammaBiasingOperator(); gammaBiasing->AttachTo(biasingLogicalVolume); - G4cout << " Attaching biasing operator " << gammaBiasing->GetName() - << " to logical volume " << biasingLogicalVolume->GetName() - << G4endl; + G4cout << " Attaching biasing operator " << gammaBiasing->GetName() << " to logical volume " + << biasingLogicalVolume->GetName() << G4endl; } // Print info cout << "Biasing is enabled " << endl; cout << "Splitting factor: " << biasingInfo.GetSplittingFactor() << endl; cout << "Biasing volumes: " << endl; - for (const auto &volume: biasingVolumes) { + for (const auto& volume : biasingVolumes) { cout << "\t" << volume << endl; } const auto center = biasingInfo.GetBiasingCenter(); @@ -346,26 +345,26 @@ void DetectorConstruction::ConstructSDandField() { } } -void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume *world) { - const auto detector = (DetectorConstruction *) G4RunManager::GetRunManager()->GetUserDetectorConstruction(); - TRestGeant4Metadata *restG4Metadata = detector->fSimulationManager->GetRestMetadata(); +void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume* world) { + const auto detector = (DetectorConstruction*)G4RunManager::GetRunManager()->GetUserDetectorConstruction(); + TRestGeant4Metadata* restG4Metadata = detector->fSimulationManager->GetRestMetadata(); const size_t n = int(world->GetLogicalVolume()->GetNoDaughters()); for (size_t i = 0; i < n + 1; i++) { // world is the + 1 - G4VPhysicalVolume *volume; + G4VPhysicalVolume* volume; if (i == n) { - volume = const_cast(world); + volume = const_cast(world); } else { volume = world->GetLogicalVolume()->GetDaughter(G4int(i)); } - TString namePhysical = (TString) volume->GetName(); + TString namePhysical = (TString)volume->GetName(); if (fGdmlNewPhysicalNames.size() > i) { // it has been filled fGeant4PhysicalNameToNewPhysicalNameMap[namePhysical] = fGdmlNewPhysicalNames[i]; } TString physicalNewName = GetAlternativeNameFromGeant4PhysicalName(namePhysical); - TString nameLogical = (TString) volume->GetLogicalVolume()->GetName(); - TString nameMaterial = (TString) volume->GetLogicalVolume()->GetMaterial()->GetName(); + TString nameLogical = (TString)volume->GetLogicalVolume()->GetName(); + TString nameMaterial = (TString)volume->GetLogicalVolume()->GetMaterial()->GetName(); auto position = volume->GetTranslation(); fPhysicalToLogicalVolumeMap[physicalNewName] = nameLogical; diff --git a/src/GammaBiasingOperation.cxx b/src/GammaBiasingOperation.cxx index 9f5d8f6c..7ae6a336 100644 --- a/src/GammaBiasingOperation.cxx +++ b/src/GammaBiasingOperation.cxx @@ -1,49 +1,37 @@ -#include - #include "GammaBiasingOperation.h" -#include "GammaBiasingOperator.h" + +#include #include "G4ParticleChangeForLoss.hh" +#include "GammaBiasingOperator.h" GammaBiasingOperation::GammaBiasingOperation(G4String name) - : G4VBiasingOperation(std::move(name)), - fSplittingFactor(1), - fParticleChange() { -} - + : G4VBiasingOperation(std::move(name)), fSplittingFactor(1), fParticleChange() {} GammaBiasingOperation::~GammaBiasingOperation() = default; -G4VParticleChange *GammaBiasingOperation::ApplyFinalStateBiasing(const G4BiasingProcessInterface *callingProcess, - const G4Track *track, - const G4Step *step, - G4bool &) { - - G4VParticleChange *processFinalState = - callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); +G4VParticleChange* GammaBiasingOperation::ApplyFinalStateBiasing( + const G4BiasingProcessInterface* callingProcess, const G4Track* track, const G4Step* step, G4bool&) { + G4VParticleChange* processFinalState = callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); if (fSplittingFactor == 1) return processFinalState; // special case: no secondaries if (processFinalState->GetNumberOfSecondaries() == 0) return processFinalState; - auto actualParticleChange = (G4ParticleChangeForLoss *) processFinalState; + auto actualParticleChange = (G4ParticleChangeForLoss*)processFinalState; fParticleChange.Initialize(*track); // -- Store electron final state: - fParticleChange. - ProposeTrackStatus(actualParticleChange->GetTrackStatus()); - fParticleChange. - ProposeEnergy(actualParticleChange->GetProposedKineticEnergy()); - fParticleChange. - ProposeMomentumDirection(actualParticleChange->GetProposedMomentumDirection()); + fParticleChange.ProposeTrackStatus(actualParticleChange->GetTrackStatus()); + fParticleChange.ProposeEnergy(actualParticleChange->GetProposedKineticEnergy()); + fParticleChange.ProposeMomentumDirection(actualParticleChange->GetProposedMomentumDirection()); fParticleChange.SetSecondaryWeightByProcess(true); - G4Track *gammaTrack = actualParticleChange->GetSecondary(0); - + G4Track* gammaTrack = actualParticleChange->GetSecondary(0); // print gamma info (energy, direction) /* @@ -65,7 +53,6 @@ G4VParticleChange *GammaBiasingOperation::ApplyFinalStateBiasing(const G4Biasing // G4cout << "Gamma weight: " << gammaWeight << " nSecondaries: " << nSecondaries << G4endl; gammaTrack->SetWeight(gammaWeight); - fParticleChange.SetNumberOfSecondaries(nSecondaries); fParticleChange.AddSecondary(gammaTrack); @@ -80,7 +67,7 @@ G4VParticleChange *GammaBiasingOperation::ApplyFinalStateBiasing(const G4Biasing fParticleChange.AddSecondary(gammaTrack); nCalls++; } - // -- very rare special case: we ignore for now. + // -- very rare special case: we ignore for now. else if (processFinalState->GetNumberOfSecondaries() > 1) { for (G4int i = 0; i < processFinalState->GetNumberOfSecondaries(); i++) delete processFinalState->GetSecondary(i); diff --git a/src/GammaBiasingOperator.cxx b/src/GammaBiasingOperator.cxx index 92bc3330..a184d6fc 100644 --- a/src/GammaBiasingOperator.cxx +++ b/src/GammaBiasingOperator.cxx @@ -2,34 +2,26 @@ // Created by lobis on 2/22/2023. // -#include "GammaBiasingOperation.h" #include "GammaBiasingOperator.h" #include "G4BiasingProcessInterface.hh" #include "G4GenericMessenger.hh" +#include "GammaBiasingOperation.h" GammaBiasingOperator::GammaBiasingOperator() - : G4VBiasingOperator("BremSplittingOperator"), - fSplittingFactor(1), - fBiasOnlyOnce(true) { + : G4VBiasingOperator("BremSplittingOperator"), fSplittingFactor(1), fBiasOnlyOnce(true) { fBremSplittingOperation = new GammaBiasingOperation("BremSplittingOperation"); } +void GammaBiasingOperator::StartRun() { fBremSplittingOperation->SetSplittingFactor(fSplittingFactor); } -void GammaBiasingOperator::StartRun() { - fBremSplittingOperation->SetSplittingFactor(fSplittingFactor); -} - -void GammaBiasingOperator::StartTracking(const G4Track *) { - fNInteractions = 0; -} +void GammaBiasingOperator::StartTracking(const G4Track*) { fNInteractions = 0; } -G4VBiasingOperation * -GammaBiasingOperator::ProposeFinalStateBiasingOperation(const G4Track *, - const G4BiasingProcessInterface *) { +G4VBiasingOperation* GammaBiasingOperator::ProposeFinalStateBiasingOperation( + const G4Track*, const G4BiasingProcessInterface*) { if (fBiasOnlyOnce && (fNInteractions > 0)) return nullptr; fNInteractions++; return fBremSplittingOperation; -} \ No newline at end of file +} diff --git a/src/PhysicsList.cxx b/src/PhysicsList.cxx index 592e3801..2dd8b5ec 100644 --- a/src/PhysicsList.cxx +++ b/src/PhysicsList.cxx @@ -3,7 +3,6 @@ #include -#include #include #include #include @@ -13,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -43,8 +43,8 @@ using namespace std; -PhysicsList::PhysicsList(SimulationManager *simulationManager, TRestGeant4PhysicsLists *physicsLists) - : G4VModularPhysicsList(), fSimulationManager(simulationManager) { +PhysicsList::PhysicsList(SimulationManager* simulationManager, TRestGeant4PhysicsLists* physicsLists) + : G4VModularPhysicsList(), fSimulationManager(simulationManager) { // add new units for radioActive decays const G4double G4minute = 60 * second; const G4double G4hour = 60 * G4minute; @@ -61,8 +61,8 @@ PhysicsList::PhysicsList(SimulationManager *simulationManager, TRestGeant4Physic G4LossTableManager::Instance(); // fix lower limit for cut G4ProductionCutsTable::GetProductionCutsTable()->SetEnergyRange( - fRestPhysicsLists->GetMinimumEnergyProductionCuts() * keV, - fRestPhysicsLists->GetMaximumEnergyProductionCuts() * keV); + fRestPhysicsLists->GetMinimumEnergyProductionCuts() * keV, + fRestPhysicsLists->GetMaximumEnergyProductionCuts() * keV); InitializePhysicsLists(); } @@ -71,7 +71,7 @@ PhysicsList::~PhysicsList() { delete fEmPhysicsList; delete fDecPhysicsList; delete fRadDecPhysicsList; - for (auto &hadronicPhysicsList: fHadronPhys) { + for (auto& hadronicPhysicsList : fHadronPhys) { delete hadronicPhysicsList; } } @@ -189,7 +189,7 @@ void PhysicsList::ConstructParticle() { fBiasingPhysicsList->ConstructParticle(); } - for (auto &hadronicPhysicsList: fHadronPhys) { + for (auto& hadronicPhysicsList : fHadronPhys) { hadronicPhysicsList->ConstructParticle(); } } @@ -202,27 +202,27 @@ void PhysicsList::ConstructProcess() { fEmPhysicsList->ConstructProcess(); fEmConfig.AddModels(); - G4UImanager *UI = G4UImanager::GetUIpointer(); + G4UImanager* UI = G4UImanager::GetUIpointer(); UI->ApplyCommand("/process/em/fluo true"); UI->ApplyCommand("/process/em/auger true"); UI->ApplyCommand("/process/em/pixe true"); bool boolEmOptionPixe = StringToBool( - fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "pixe", "false").Data()); + fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "pixe", "false").Data()); string stringEmOptionPixe = (boolEmOptionPixe ? "true" : "false"); G4cout << "Setting EM option '/process/em/pixe' to '" << stringEmOptionPixe << "' for physics list '" << fEmPhysicsListName << "'" << endl; UI->ApplyCommand(string("/process/em/pixe ") + stringEmOptionPixe); bool boolEmOptionFluo = StringToBool( - fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "fluo", "true").Data()); + fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "fluo", "true").Data()); string stringEmOptionFluo = (boolEmOptionFluo ? "true" : "false"); G4cout << "Setting EM option '/process/em/fluo' to '" << stringEmOptionFluo << "' for physics list '" << fEmPhysicsListName << "'" << endl; UI->ApplyCommand(string("/process/em/fluo ") + stringEmOptionFluo); bool boolEmOptionAuger = StringToBool( - fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "auger", "true").Data()); + fRestPhysicsLists->GetPhysicsListOptionValue(fEmPhysicsListName.c_str(), "auger", "true").Data()); string stringEmOptionAuger = (boolEmOptionAuger ? "true" : "false"); G4cout << "Setting EM option '/process/em/auger' to '" << stringEmOptionAuger << "' for physics list '" << fEmPhysicsListName << "'" << endl; @@ -244,7 +244,7 @@ void PhysicsList::ConstructProcess() { } // Hadronic physics lists - for (auto &hadronicPhysicsList: fHadronPhys) { + for (auto& hadronicPhysicsList : fHadronPhys) { hadronicPhysicsList->ConstructProcess(); } @@ -296,15 +296,15 @@ void PhysicsList::ConstructProcess() { if (!(fRestPhysicsLists->GetPhysicsListOptionValue("G4RadioactiveDecay", "TritiumDecay") == "false") && (fRestPhysicsLists->GetPhysicsListOptionValue("G4RadioactiveDecay", "TritiumDecay", "false") == - "true" || + "true" || fSimulationManager->GetRestMetadata()->GetParticleSource()->GetParticleName() == "H3")) { // Tritium (H3) fix | https://geant4-forum.web.cern.ch/t/triton-decay-with-rdecay01-example/2616 - G4ParticleDefinition *tritium = G4Triton::Definition(); + G4ParticleDefinition* tritium = G4Triton::Definition(); tritium->SetPDGStable(false); - G4VProcess *decay = nullptr; - G4ProcessManager *tritiumProcessManager = tritium->GetProcessManager(); - G4ProcessVector *tritiumProcessVector = tritiumProcessManager->GetAtRestProcessVector(); + G4VProcess* decay = nullptr; + G4ProcessManager* tritiumProcessManager = tritium->GetProcessManager(); + G4ProcessVector* tritiumProcessVector = tritiumProcessManager->GetAtRestProcessVector(); for (unsigned int i = 0; i < tritiumProcessVector->size() && decay == nullptr; i++) { if ((*tritiumProcessVector)[i]->GetProcessName() == "Decay") decay = (*tritiumProcessVector)[i]; @@ -325,9 +325,9 @@ void PhysicsList::ConstructProcess() { // To implement UserLimits to StepSize inside the gas theParticleIterator->reset(); while ((*theParticleIterator)()) { - G4ParticleDefinition *particle = theParticleIterator->value(); - const auto &particleName = particle->GetParticleName(); - G4ProcessManager *processManager = particle->GetProcessManager(); + G4ParticleDefinition* particle = theParticleIterator->value(); + const auto& particleName = particle->GetParticleName(); + G4ProcessManager* processManager = particle->GetProcessManager(); if (particleName == "e-") { processManager->AddDiscreteProcess(new G4StepLimiter("e-Step")); @@ -347,10 +347,10 @@ void PhysicsList::ConstructProcess() { for (int A = 2 * Z; A <= 3 * Z; A++) { for (unsigned int n = 0; n < fRestPhysicsLists->GetIonStepList().size(); n++) { if (fRestPhysicsLists->GetIonStepList()[n] == G4IonTable::GetIonTable()->GetIonName(Z, A)) { - G4ParticleDefinition *particle = G4IonTable::GetIonTable()->GetIon(Z, A, 0); + G4ParticleDefinition* particle = G4IonTable::GetIonTable()->GetIon(Z, A, 0); G4String particle_name = G4IonTable::GetIonTable()->GetIonName(Z, A, 0); cout << "Found ion: " << particle_name << " Z " << Z << " A " << A << endl; - G4ProcessManager *processManager = particle->GetProcessManager(); + G4ProcessManager* processManager = particle->GetProcessManager(); processManager->AddDiscreteProcess(new G4StepLimiter("ionStep")); } } diff --git a/src/TrackingAction.cxx b/src/TrackingAction.cxx index 095a7b9d..47a7020a 100644 --- a/src/TrackingAction.cxx +++ b/src/TrackingAction.cxx @@ -9,14 +9,14 @@ using namespace std; -TrackingAction::TrackingAction(SimulationManager *simulationManager) - : G4UserTrackingAction(), fSimulationManager(simulationManager) {} +TrackingAction::TrackingAction(SimulationManager* simulationManager) + : G4UserTrackingAction(), fSimulationManager(simulationManager) {} -void TrackingAction::PreUserTrackingAction(const G4Track *track) { +void TrackingAction::PreUserTrackingAction(const G4Track* track) { // G4cout << track->GetVolume()->GetLogicalVolume()->GetName() << G4endl; fSimulationManager->GetOutputManager()->RecordTrack(track); } -void TrackingAction::PostUserTrackingAction(const G4Track *track) { +void TrackingAction::PostUserTrackingAction(const G4Track* track) { fSimulationManager->GetOutputManager()->UpdateTrack(track); } From a5482dfb5b8b7ddc554e9cd1f8e83b0c0dfba510 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Wed, 22 Feb 2023 18:53:20 +0100 Subject: [PATCH 09/24] attempt to fix pipeline failure --- examples/13.IAXO/GammaBiasing.rml | 6 ------ include/GammaBiasingOperation.h | 5 ++++- include/GammaBiasingOperator.h | 7 ++++++- src/DetectorConstruction.cxx | 5 ++++- src/GammaBiasingOperation.cxx | 8 ++++++-- src/GammaBiasingOperator.cxx | 11 ++++++++--- src/PhysicsList.cxx | 11 +++++++---- 7 files changed, 35 insertions(+), 18 deletions(-) diff --git a/examples/13.IAXO/GammaBiasing.rml b/examples/13.IAXO/GammaBiasing.rml index 0a6478b4..d6f2faa9 100644 --- a/examples/13.IAXO/GammaBiasing.rml +++ b/examples/13.IAXO/GammaBiasing.rml @@ -22,12 +22,6 @@ Author: Luis Obis (lobis@unizar.es) - - - - - - diff --git a/include/GammaBiasingOperation.h b/include/GammaBiasingOperation.h index 75dddba7..3db6c990 100644 --- a/include/GammaBiasingOperation.h +++ b/include/GammaBiasingOperation.h @@ -5,13 +5,15 @@ #ifndef REST_GAMMABIASINGOPERATION_H #define REST_GAMMABIASINGOPERATION_H +#include + #include "G4BiasingProcessInterface.hh" #include "G4ParticleChange.hh" #include "G4VBiasingOperation.hh" class GammaBiasingOperation : public G4VBiasingOperation { public: - GammaBiasingOperation(G4String name); + GammaBiasingOperation(const G4String& name, G4int splittingFactor, const TVector3& biasingCenter); virtual ~GammaBiasingOperation(); @@ -36,6 +38,7 @@ class GammaBiasingOperation : public G4VBiasingOperation { private: G4int fSplittingFactor; G4ParticleChange fParticleChange; + TVector3 fBiasingCenter; }; #endif // REST_GAMMABIASINGOPERATION_H diff --git a/include/GammaBiasingOperator.h b/include/GammaBiasingOperator.h index 7a845db2..d4fa40fd 100644 --- a/include/GammaBiasingOperator.h +++ b/include/GammaBiasingOperator.h @@ -5,13 +5,15 @@ #ifndef REST_GAMMABIASINGOPERATOR_H #define REST_GAMMABIASINGOPERATOR_H +#include + #include "G4VBiasingOperator.hh" class GammaBiasingOperation; class GammaBiasingOperator : public G4VBiasingOperator { public: - GammaBiasingOperator(); + GammaBiasingOperator(int splittingFactor, bool biasOnlyOnce, const TVector3& biasingCenter); virtual ~GammaBiasingOperator() {} @@ -20,6 +22,8 @@ class GammaBiasingOperator : public G4VBiasingOperator { virtual void StartTracking(const G4Track* track); + TVector3 GetBiasingCenter() const { return fBiasingCenter; } + private: virtual G4VBiasingOperation* ProposeNonPhysicsBiasingOperation(const G4Track*, const G4BiasingProcessInterface*) { @@ -42,6 +46,7 @@ class GammaBiasingOperator : public G4VBiasingOperator { G4int fSplittingFactor; G4bool fBiasOnlyOnce; G4int fNInteractions; + TVector3 fBiasingCenter; }; #endif // REST_GAMMABIASINGOPERATOR_H diff --git a/src/DetectorConstruction.cxx b/src/DetectorConstruction.cxx index 3dc6ad49..8e918b06 100644 --- a/src/DetectorConstruction.cxx +++ b/src/DetectorConstruction.cxx @@ -314,6 +314,8 @@ void DetectorConstruction::ConstructSDandField() { << RESTendl; exit(1); } + } else { + biasingVolumes.insert(name); } } @@ -327,7 +329,8 @@ void DetectorConstruction::ConstructSDandField() { } // TODO: better memory management - auto gammaBiasing = new GammaBiasingOperator(); + auto gammaBiasing = new GammaBiasingOperator(biasingInfo.GetSplittingFactor(), true, + biasingInfo.GetBiasingCenter()); gammaBiasing->AttachTo(biasingLogicalVolume); G4cout << " Attaching biasing operator " << gammaBiasing->GetName() << " to logical volume " << biasingLogicalVolume->GetName() << G4endl; diff --git a/src/GammaBiasingOperation.cxx b/src/GammaBiasingOperation.cxx index 7ae6a336..f89e064c 100644 --- a/src/GammaBiasingOperation.cxx +++ b/src/GammaBiasingOperation.cxx @@ -6,8 +6,12 @@ #include "G4ParticleChangeForLoss.hh" #include "GammaBiasingOperator.h" -GammaBiasingOperation::GammaBiasingOperation(G4String name) - : G4VBiasingOperation(std::move(name)), fSplittingFactor(1), fParticleChange() {} +GammaBiasingOperation::GammaBiasingOperation(const G4String& name, int splittingFactor, + const TVector3& biasingCenter) + : G4VBiasingOperation(name), + fSplittingFactor(splittingFactor), + fBiasingCenter(biasingCenter), + fParticleChange() {} GammaBiasingOperation::~GammaBiasingOperation() = default; diff --git a/src/GammaBiasingOperator.cxx b/src/GammaBiasingOperator.cxx index a184d6fc..14e37a79 100644 --- a/src/GammaBiasingOperator.cxx +++ b/src/GammaBiasingOperator.cxx @@ -8,9 +8,14 @@ #include "G4GenericMessenger.hh" #include "GammaBiasingOperation.h" -GammaBiasingOperator::GammaBiasingOperator() - : G4VBiasingOperator("BremSplittingOperator"), fSplittingFactor(1), fBiasOnlyOnce(true) { - fBremSplittingOperation = new GammaBiasingOperation("BremSplittingOperation"); +GammaBiasingOperator::GammaBiasingOperator(int splittingFactor, bool biasOnlyOnce, + const TVector3& biasingCenter) + : G4VBiasingOperator("BremSplittingOperator"), + fSplittingFactor(splittingFactor), + fBiasOnlyOnce(biasOnlyOnce), + fBiasingCenter(biasingCenter) { + fBremSplittingOperation = + new GammaBiasingOperation("BremSplittingOperation", splittingFactor, biasingCenter); } void GammaBiasingOperator::StartRun() { fBremSplittingOperation->SetSplittingFactor(fSplittingFactor); } diff --git a/src/PhysicsList.cxx b/src/PhysicsList.cxx index 2dd8b5ec..f8051f9a 100644 --- a/src/PhysicsList.cxx +++ b/src/PhysicsList.cxx @@ -129,10 +129,13 @@ void PhysicsList::InitializePhysicsLists() { emCounter++; } - fBiasingPhysicsList = new G4GenericBiasingPhysics(); - std::vector processToBias = {"eBrem"}; - fBiasingPhysicsList->PhysicsBias("e-", processToBias); - fBiasingPhysicsList->PhysicsBias("e+", processToBias); + const auto& biasingInfo = fSimulationManager->GetRestMetadata()->GetGeant4BiasingInfo(); + if (biasingInfo.GetEnabled()) { + fBiasingPhysicsList = new G4GenericBiasingPhysics(); + std::vector processToBias = {"eBrem"}; + fBiasingPhysicsList->PhysicsBias("e-", processToBias); + fBiasingPhysicsList->PhysicsBias("e+", processToBias); + } if (fRestPhysicsLists->GetVerboseLevel() >= TRestStringOutput::REST_Verbose_Level::REST_Essential && emCounter == 0) { From 2c660564cceefc13110e700753f35bf06a5513d4 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Thu, 23 Feb 2023 13:27:29 +0100 Subject: [PATCH 10/24] support for old biasing --- examples/13.IAXO/GammaBiasing.rml | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/examples/13.IAXO/GammaBiasing.rml b/examples/13.IAXO/GammaBiasing.rml index d6f2faa9..b9f3274f 100644 --- a/examples/13.IAXO/GammaBiasing.rml +++ b/examples/13.IAXO/GammaBiasing.rml @@ -31,6 +31,13 @@ Author: Luis Obis (lobis@unizar.es) + + + + + + + From 5165af18105bb2dbaf0619fe6d265160cf9a00b6 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Mon, 27 Feb 2023 09:47:40 +0100 Subject: [PATCH 11/24] update constructor order --- src/GammaBiasingOperation.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/GammaBiasingOperation.cxx b/src/GammaBiasingOperation.cxx index f89e064c..9467b8a2 100644 --- a/src/GammaBiasingOperation.cxx +++ b/src/GammaBiasingOperation.cxx @@ -10,8 +10,8 @@ GammaBiasingOperation::GammaBiasingOperation(const G4String& name, int splitting const TVector3& biasingCenter) : G4VBiasingOperation(name), fSplittingFactor(splittingFactor), - fBiasingCenter(biasingCenter), - fParticleChange() {} + fParticleChange(), + fBiasingCenter(biasingCenter) {} GammaBiasingOperation::~GammaBiasingOperation() = default; From e61475a3972357f5ceedc66185bd664b43c5304d Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Mon, 27 Feb 2023 10:47:56 +0100 Subject: [PATCH 12/24] update volume names --- examples/13.IAXO/geometry/setup.gdml | 14 +++++++------- 1 file changed, 7 insertions(+), 7 deletions(-) diff --git a/examples/13.IAXO/geometry/setup.gdml b/examples/13.IAXO/geometry/setup.gdml index 4279c32c..acf7c7f4 100644 --- a/examples/13.IAXO/geometry/setup.gdml +++ b/examples/13.IAXO/geometry/setup.gdml @@ -465,7 +465,7 @@ - + @@ -474,13 +474,13 @@ - + - + - + @@ -895,11 +895,11 @@ - + - + @@ -912,4 +912,4 @@ - + \ No newline at end of file From 5558cccaacb2486eb97091ddac3650c139cca811 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 27 Feb 2023 09:48:12 +0000 Subject: [PATCH 13/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- examples/13.IAXO/geometry/setup.gdml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/13.IAXO/geometry/setup.gdml b/examples/13.IAXO/geometry/setup.gdml index acf7c7f4..7572fe61 100644 --- a/examples/13.IAXO/geometry/setup.gdml +++ b/examples/13.IAXO/geometry/setup.gdml @@ -912,4 +912,4 @@ - \ No newline at end of file + From b6ac1aaba67b895510ca5cec7e2c65a9ede8e0f9 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Mon, 27 Feb 2023 11:01:48 +0100 Subject: [PATCH 14/24] reduce number of backwards tracks --- examples/13.IAXO/GammaBiasing.rml | 4 +- src/GammaBiasingOperation.cxx | 65 ++++++++++++++++++------------- 2 files changed, 39 insertions(+), 30 deletions(-) diff --git a/examples/13.IAXO/GammaBiasing.rml b/examples/13.IAXO/GammaBiasing.rml index b9f3274f..d30c012d 100644 --- a/examples/13.IAXO/GammaBiasing.rml +++ b/examples/13.IAXO/GammaBiasing.rml @@ -33,10 +33,10 @@ Author: Luis Obis (lobis@unizar.es) - + - + l diff --git a/src/GammaBiasingOperation.cxx b/src/GammaBiasingOperation.cxx index 9467b8a2..41d74c9a 100644 --- a/src/GammaBiasingOperation.cxx +++ b/src/GammaBiasingOperation.cxx @@ -6,25 +6,25 @@ #include "G4ParticleChangeForLoss.hh" #include "GammaBiasingOperator.h" -GammaBiasingOperation::GammaBiasingOperation(const G4String& name, int splittingFactor, - const TVector3& biasingCenter) - : G4VBiasingOperation(name), - fSplittingFactor(splittingFactor), - fParticleChange(), - fBiasingCenter(biasingCenter) {} +GammaBiasingOperation::GammaBiasingOperation(const G4String &name, int splittingFactor, + const TVector3 &biasingCenter) + : G4VBiasingOperation(name), + fSplittingFactor(splittingFactor), + fParticleChange(), + fBiasingCenter(biasingCenter) {} GammaBiasingOperation::~GammaBiasingOperation() = default; -G4VParticleChange* GammaBiasingOperation::ApplyFinalStateBiasing( - const G4BiasingProcessInterface* callingProcess, const G4Track* track, const G4Step* step, G4bool&) { - G4VParticleChange* processFinalState = callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); +G4VParticleChange *GammaBiasingOperation::ApplyFinalStateBiasing( + const G4BiasingProcessInterface *callingProcess, const G4Track *track, const G4Step *step, G4bool &) { + G4VParticleChange *processFinalState = callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); if (fSplittingFactor == 1) return processFinalState; // special case: no secondaries if (processFinalState->GetNumberOfSecondaries() == 0) return processFinalState; - auto actualParticleChange = (G4ParticleChangeForLoss*)processFinalState; + auto actualParticleChange = (G4ParticleChangeForLoss *) processFinalState; fParticleChange.Initialize(*track); @@ -35,30 +35,39 @@ G4VParticleChange* GammaBiasingOperation::ApplyFinalStateBiasing( fParticleChange.SetSecondaryWeightByProcess(true); - G4Track* gammaTrack = actualParticleChange->GetSecondary(0); - - // print gamma info (energy, direction) - /* - G4cout << "Gamma info. Weight: " << gammaTrack->GetWeight() - << " Energy (keV): " << gammaTrack->GetKineticEnergy() / CLHEP::keV << " Direction: " - << gammaTrack->GetMomentumDirection() << G4endl; - */ - // if direction points towards (0,0,0) - bool split = false; - const auto diff = gammaTrack->GetPosition() - G4ThreeVector(0, 0, 0); + G4Track *gammaTrack = actualParticleChange->GetSecondary(0); + + int nSecondaries = 1; + double weightFactor = 1.0; + + const auto diff = + gammaTrack->GetPosition() - G4ThreeVector(fBiasingCenter.X(), fBiasingCenter.Y(), fBiasingCenter.Z()); if (gammaTrack->GetMomentumDirection().dot(diff) < 0) { - // G4cout << "Gamma points towards (0,0,0)" << G4endl; - split = true; + // pointing towards point of interest + nSecondaries = fSplittingFactor; + weightFactor = 1.0 / fSplittingFactor; + } else { + // pointing away from point of interest + // random number between 0 and 1 + const double rand = G4UniformRand(); + // if random number is less than 1 / fSplittingFactor, keep alive + if (rand < 1.0 / fSplittingFactor) { + nSecondaries = 1; + weightFactor = fSplittingFactor; // increase weight + } else { + nSecondaries = 0; + weightFactor = 0; + } } - const int nSecondaries = split ? fSplittingFactor : 1; - - G4double gammaWeight = track->GetWeight() / nSecondaries; + G4double gammaWeight = track->GetWeight() * weightFactor; // G4cout << "Gamma weight: " << gammaWeight << " nSecondaries: " << nSecondaries << G4endl; gammaTrack->SetWeight(gammaWeight); fParticleChange.SetNumberOfSecondaries(nSecondaries); - fParticleChange.AddSecondary(gammaTrack); + if (nSecondaries > 0) { + fParticleChange.AddSecondary(gammaTrack); + } actualParticleChange->Clear(); @@ -71,7 +80,7 @@ G4VParticleChange* GammaBiasingOperation::ApplyFinalStateBiasing( fParticleChange.AddSecondary(gammaTrack); nCalls++; } - // -- very rare special case: we ignore for now. + // -- very rare special case: we ignore for now. else if (processFinalState->GetNumberOfSecondaries() > 1) { for (G4int i = 0; i < processFinalState->GetNumberOfSecondaries(); i++) delete processFinalState->GetSecondary(i); From fb596cb2189b916d296a9af250ed64fcfb57c8e1 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Mon, 27 Feb 2023 11:02:26 +0100 Subject: [PATCH 15/24] fix bad volume name --- examples/13.IAXO/GammaBiasing.rml | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/examples/13.IAXO/GammaBiasing.rml b/examples/13.IAXO/GammaBiasing.rml index d30c012d..347b1b31 100644 --- a/examples/13.IAXO/GammaBiasing.rml +++ b/examples/13.IAXO/GammaBiasing.rml @@ -29,14 +29,15 @@ Author: Luis Obis (lobis@unizar.es) - + - l + + l From fd1e9f80b418564c94b6d796318d94e29dd5ab8a Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Mon, 27 Feb 2023 11:02:58 +0100 Subject: [PATCH 16/24] fix bad volume name --- examples/13.IAXO/Neutrons.rml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/examples/13.IAXO/Neutrons.rml b/examples/13.IAXO/Neutrons.rml index c4dde63b..a35b32e4 100644 --- a/examples/13.IAXO/Neutrons.rml +++ b/examples/13.IAXO/Neutrons.rml @@ -29,7 +29,7 @@ Author: Luis Obis (lobis@unizar.es) - + From 84a86054331a29a4d002e0d64d1c44224722dee7 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 27 Feb 2023 10:04:28 +0000 Subject: [PATCH 17/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/GammaBiasingOperation.cxx | 28 ++++++++++++++-------------- 1 file changed, 14 insertions(+), 14 deletions(-) diff --git a/src/GammaBiasingOperation.cxx b/src/GammaBiasingOperation.cxx index 41d74c9a..ce8d8d5f 100644 --- a/src/GammaBiasingOperation.cxx +++ b/src/GammaBiasingOperation.cxx @@ -6,25 +6,25 @@ #include "G4ParticleChangeForLoss.hh" #include "GammaBiasingOperator.h" -GammaBiasingOperation::GammaBiasingOperation(const G4String &name, int splittingFactor, - const TVector3 &biasingCenter) - : G4VBiasingOperation(name), - fSplittingFactor(splittingFactor), - fParticleChange(), - fBiasingCenter(biasingCenter) {} +GammaBiasingOperation::GammaBiasingOperation(const G4String& name, int splittingFactor, + const TVector3& biasingCenter) + : G4VBiasingOperation(name), + fSplittingFactor(splittingFactor), + fParticleChange(), + fBiasingCenter(biasingCenter) {} GammaBiasingOperation::~GammaBiasingOperation() = default; -G4VParticleChange *GammaBiasingOperation::ApplyFinalStateBiasing( - const G4BiasingProcessInterface *callingProcess, const G4Track *track, const G4Step *step, G4bool &) { - G4VParticleChange *processFinalState = callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); +G4VParticleChange* GammaBiasingOperation::ApplyFinalStateBiasing( + const G4BiasingProcessInterface* callingProcess, const G4Track* track, const G4Step* step, G4bool&) { + G4VParticleChange* processFinalState = callingProcess->GetWrappedProcess()->PostStepDoIt(*track, *step); if (fSplittingFactor == 1) return processFinalState; // special case: no secondaries if (processFinalState->GetNumberOfSecondaries() == 0) return processFinalState; - auto actualParticleChange = (G4ParticleChangeForLoss *) processFinalState; + auto actualParticleChange = (G4ParticleChangeForLoss*)processFinalState; fParticleChange.Initialize(*track); @@ -35,13 +35,13 @@ G4VParticleChange *GammaBiasingOperation::ApplyFinalStateBiasing( fParticleChange.SetSecondaryWeightByProcess(true); - G4Track *gammaTrack = actualParticleChange->GetSecondary(0); + G4Track* gammaTrack = actualParticleChange->GetSecondary(0); int nSecondaries = 1; double weightFactor = 1.0; const auto diff = - gammaTrack->GetPosition() - G4ThreeVector(fBiasingCenter.X(), fBiasingCenter.Y(), fBiasingCenter.Z()); + gammaTrack->GetPosition() - G4ThreeVector(fBiasingCenter.X(), fBiasingCenter.Y(), fBiasingCenter.Z()); if (gammaTrack->GetMomentumDirection().dot(diff) < 0) { // pointing towards point of interest nSecondaries = fSplittingFactor; @@ -53,7 +53,7 @@ G4VParticleChange *GammaBiasingOperation::ApplyFinalStateBiasing( // if random number is less than 1 / fSplittingFactor, keep alive if (rand < 1.0 / fSplittingFactor) { nSecondaries = 1; - weightFactor = fSplittingFactor; // increase weight + weightFactor = fSplittingFactor; // increase weight } else { nSecondaries = 0; weightFactor = 0; @@ -80,7 +80,7 @@ G4VParticleChange *GammaBiasingOperation::ApplyFinalStateBiasing( fParticleChange.AddSecondary(gammaTrack); nCalls++; } - // -- very rare special case: we ignore for now. + // -- very rare special case: we ignore for now. else if (processFinalState->GetNumberOfSecondaries() > 1) { for (G4int i = 0; i < processFinalState->GetNumberOfSecondaries(); i++) delete processFinalState->GetSecondary(i); From 7a63d7cdd86155599529a1ab28a76973207aad34 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Mon, 27 Feb 2023 11:26:45 +0100 Subject: [PATCH 18/24] grouped biasing in a single header/src file --- include/Biasing.h | 81 +++++++++++++++++++ include/GammaBiasingOperation.h | 44 ---------- include/GammaBiasingOperator.h | 52 ------------ ...{GammaBiasingOperation.cxx => Biasing.cxx} | 37 +++++++-- src/GammaBiasingOperator.cxx | 32 -------- 5 files changed, 112 insertions(+), 134 deletions(-) create mode 100644 include/Biasing.h delete mode 100644 include/GammaBiasingOperation.h delete mode 100644 include/GammaBiasingOperator.h rename src/{GammaBiasingOperation.cxx => Biasing.cxx} (74%) delete mode 100644 src/GammaBiasingOperator.cxx diff --git a/include/Biasing.h b/include/Biasing.h new file mode 100644 index 00000000..600f2d55 --- /dev/null +++ b/include/Biasing.h @@ -0,0 +1,81 @@ + +#ifndef REST_BIASING_H +#define REST_BIASING_H + +#include + +#include +#include +#include +#include + +class GammaBiasingOperation : public G4VBiasingOperation { +public: + GammaBiasingOperation(const G4String &name, G4int splittingFactor, const TVector3 &biasingCenter); + + virtual ~GammaBiasingOperation(); + +public: + virtual const G4VBiasingInteractionLaw *ProvideOccurenceBiasingInteractionLaw( + const G4BiasingProcessInterface *, G4ForceCondition &) { + return 0; + } + + virtual G4VParticleChange *ApplyFinalStateBiasing(const G4BiasingProcessInterface *, const G4Track *, + const G4Step *, G4bool &); + + virtual G4double DistanceToApplyOperation(const G4Track *, G4double, G4ForceCondition *) { return DBL_MAX; } + + virtual G4VParticleChange *GenerateBiasingFinalState(const G4Track *, const G4Step *) { return 0; } + +public: + void SetSplittingFactor(G4int splittingFactor) { fSplittingFactor = splittingFactor; } + + G4int GetSplittingFactor() const { return fSplittingFactor; } + +private: + G4int fSplittingFactor; + G4ParticleChange fParticleChange; + TVector3 fBiasingCenter; +}; + + +class GammaBiasingOperator : public G4VBiasingOperator { +public: + GammaBiasingOperator(int splittingFactor, bool biasOnlyOnce, const TVector3 &biasingCenter); + + virtual ~GammaBiasingOperator() {} + +public: + virtual void StartRun(); + + virtual void StartTracking(const G4Track *track); + + TVector3 GetBiasingCenter() const { return fBiasingCenter; } + +private: + virtual G4VBiasingOperation *ProposeNonPhysicsBiasingOperation(const G4Track *, + const G4BiasingProcessInterface *) { + return 0; + } + + virtual G4VBiasingOperation *ProposeOccurenceBiasingOperation(const G4Track *, + const G4BiasingProcessInterface *) { + return 0; + } + + virtual G4VBiasingOperation *ProposeFinalStateBiasingOperation( + const G4Track *track, const G4BiasingProcessInterface *callingProcess); + +private: + using G4VBiasingOperator::OperationApplied; + +private: + GammaBiasingOperation *fBremSplittingOperation; + G4int fSplittingFactor; + G4bool fBiasOnlyOnce; + G4int fNInteractions; + TVector3 fBiasingCenter; +}; + +#endif //REST_BIASING_H diff --git a/include/GammaBiasingOperation.h b/include/GammaBiasingOperation.h deleted file mode 100644 index 3db6c990..00000000 --- a/include/GammaBiasingOperation.h +++ /dev/null @@ -1,44 +0,0 @@ -// -// Created by lobis on 2/22/2023. -// - -#ifndef REST_GAMMABIASINGOPERATION_H -#define REST_GAMMABIASINGOPERATION_H - -#include - -#include "G4BiasingProcessInterface.hh" -#include "G4ParticleChange.hh" -#include "G4VBiasingOperation.hh" - -class GammaBiasingOperation : public G4VBiasingOperation { - public: - GammaBiasingOperation(const G4String& name, G4int splittingFactor, const TVector3& biasingCenter); - - virtual ~GammaBiasingOperation(); - - public: - virtual const G4VBiasingInteractionLaw* ProvideOccurenceBiasingInteractionLaw( - const G4BiasingProcessInterface*, G4ForceCondition&) { - return 0; - } - - virtual G4VParticleChange* ApplyFinalStateBiasing(const G4BiasingProcessInterface*, const G4Track*, - const G4Step*, G4bool&); - - virtual G4double DistanceToApplyOperation(const G4Track*, G4double, G4ForceCondition*) { return DBL_MAX; } - - virtual G4VParticleChange* GenerateBiasingFinalState(const G4Track*, const G4Step*) { return 0; } - - public: - void SetSplittingFactor(G4int splittingFactor) { fSplittingFactor = splittingFactor; } - - G4int GetSplittingFactor() const { return fSplittingFactor; } - - private: - G4int fSplittingFactor; - G4ParticleChange fParticleChange; - TVector3 fBiasingCenter; -}; - -#endif // REST_GAMMABIASINGOPERATION_H diff --git a/include/GammaBiasingOperator.h b/include/GammaBiasingOperator.h deleted file mode 100644 index d4fa40fd..00000000 --- a/include/GammaBiasingOperator.h +++ /dev/null @@ -1,52 +0,0 @@ -// -// Created by lobis on 2/22/2023. -// - -#ifndef REST_GAMMABIASINGOPERATOR_H -#define REST_GAMMABIASINGOPERATOR_H - -#include - -#include "G4VBiasingOperator.hh" - -class GammaBiasingOperation; - -class GammaBiasingOperator : public G4VBiasingOperator { - public: - GammaBiasingOperator(int splittingFactor, bool biasOnlyOnce, const TVector3& biasingCenter); - - virtual ~GammaBiasingOperator() {} - - public: - virtual void StartRun(); - - virtual void StartTracking(const G4Track* track); - - TVector3 GetBiasingCenter() const { return fBiasingCenter; } - - private: - virtual G4VBiasingOperation* ProposeNonPhysicsBiasingOperation(const G4Track*, - const G4BiasingProcessInterface*) { - return 0; - } - - virtual G4VBiasingOperation* ProposeOccurenceBiasingOperation(const G4Track*, - const G4BiasingProcessInterface*) { - return 0; - } - - virtual G4VBiasingOperation* ProposeFinalStateBiasingOperation( - const G4Track* track, const G4BiasingProcessInterface* callingProcess); - - private: - using G4VBiasingOperator::OperationApplied; - - private: - GammaBiasingOperation* fBremSplittingOperation; - G4int fSplittingFactor; - G4bool fBiasOnlyOnce; - G4int fNInteractions; - TVector3 fBiasingCenter; -}; - -#endif // REST_GAMMABIASINGOPERATOR_H diff --git a/src/GammaBiasingOperation.cxx b/src/Biasing.cxx similarity index 74% rename from src/GammaBiasingOperation.cxx rename to src/Biasing.cxx index 41d74c9a..ea64e86c 100644 --- a/src/GammaBiasingOperation.cxx +++ b/src/Biasing.cxx @@ -1,10 +1,9 @@ -#include "GammaBiasingOperation.h" +#include "Biasing.h" -#include +#include -#include "G4ParticleChangeForLoss.hh" -#include "GammaBiasingOperator.h" +// Operation GammaBiasingOperation::GammaBiasingOperation(const G4String &name, int splittingFactor, const TVector3 &biasingCenter) @@ -51,12 +50,11 @@ G4VParticleChange *GammaBiasingOperation::ApplyFinalStateBiasing( // random number between 0 and 1 const double rand = G4UniformRand(); // if random number is less than 1 / fSplittingFactor, keep alive + weightFactor = fSplittingFactor; // increase weight if (rand < 1.0 / fSplittingFactor) { nSecondaries = 1; - weightFactor = fSplittingFactor; // increase weight } else { nSecondaries = 0; - weightFactor = 0; } } @@ -90,3 +88,30 @@ G4VParticleChange *GammaBiasingOperation::ApplyFinalStateBiasing( return &fParticleChange; } + +// Operator + +#include "G4BiasingProcessInterface.hh" + +GammaBiasingOperator::GammaBiasingOperator(int splittingFactor, bool biasOnlyOnce, + const TVector3 &biasingCenter) + : G4VBiasingOperator("BremSplittingOperator"), + fSplittingFactor(splittingFactor), + fBiasOnlyOnce(biasOnlyOnce), + fBiasingCenter(biasingCenter) { + fBremSplittingOperation = + new GammaBiasingOperation("BremSplittingOperation", splittingFactor, biasingCenter); +} + +void GammaBiasingOperator::StartRun() { fBremSplittingOperation->SetSplittingFactor(fSplittingFactor); } + +void GammaBiasingOperator::StartTracking(const G4Track *) { fNInteractions = 0; } + +G4VBiasingOperation *GammaBiasingOperator::ProposeFinalStateBiasingOperation( + const G4Track *, const G4BiasingProcessInterface *) { + if (fBiasOnlyOnce && (fNInteractions > 0)) return nullptr; + + fNInteractions++; + + return fBremSplittingOperation; +} diff --git a/src/GammaBiasingOperator.cxx b/src/GammaBiasingOperator.cxx deleted file mode 100644 index 14e37a79..00000000 --- a/src/GammaBiasingOperator.cxx +++ /dev/null @@ -1,32 +0,0 @@ -// -// Created by lobis on 2/22/2023. -// - -#include "GammaBiasingOperator.h" - -#include "G4BiasingProcessInterface.hh" -#include "G4GenericMessenger.hh" -#include "GammaBiasingOperation.h" - -GammaBiasingOperator::GammaBiasingOperator(int splittingFactor, bool biasOnlyOnce, - const TVector3& biasingCenter) - : G4VBiasingOperator("BremSplittingOperator"), - fSplittingFactor(splittingFactor), - fBiasOnlyOnce(biasOnlyOnce), - fBiasingCenter(biasingCenter) { - fBremSplittingOperation = - new GammaBiasingOperation("BremSplittingOperation", splittingFactor, biasingCenter); -} - -void GammaBiasingOperator::StartRun() { fBremSplittingOperation->SetSplittingFactor(fSplittingFactor); } - -void GammaBiasingOperator::StartTracking(const G4Track*) { fNInteractions = 0; } - -G4VBiasingOperation* GammaBiasingOperator::ProposeFinalStateBiasingOperation( - const G4Track*, const G4BiasingProcessInterface*) { - if (fBiasOnlyOnce && (fNInteractions > 0)) return nullptr; - - fNInteractions++; - - return fBremSplittingOperation; -} From 1a9cc40e3b289971989fb3e1fcd7eef930c8fc0d Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 27 Feb 2023 10:30:50 +0000 Subject: [PATCH 19/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- include/Biasing.h | 53 +++++++------- src/Biasing.cxx | 20 +++--- src/DetectorConstruction.cxx | 132 +++++++++++++++++------------------ 3 files changed, 102 insertions(+), 103 deletions(-) diff --git a/include/Biasing.h b/include/Biasing.h index 600f2d55..3b5f3aa4 100644 --- a/include/Biasing.h +++ b/include/Biasing.h @@ -10,72 +10,71 @@ #include class GammaBiasingOperation : public G4VBiasingOperation { -public: - GammaBiasingOperation(const G4String &name, G4int splittingFactor, const TVector3 &biasingCenter); + public: + GammaBiasingOperation(const G4String& name, G4int splittingFactor, const TVector3& biasingCenter); virtual ~GammaBiasingOperation(); -public: - virtual const G4VBiasingInteractionLaw *ProvideOccurenceBiasingInteractionLaw( - const G4BiasingProcessInterface *, G4ForceCondition &) { + public: + virtual const G4VBiasingInteractionLaw* ProvideOccurenceBiasingInteractionLaw( + const G4BiasingProcessInterface*, G4ForceCondition&) { return 0; } - virtual G4VParticleChange *ApplyFinalStateBiasing(const G4BiasingProcessInterface *, const G4Track *, - const G4Step *, G4bool &); + virtual G4VParticleChange* ApplyFinalStateBiasing(const G4BiasingProcessInterface*, const G4Track*, + const G4Step*, G4bool&); - virtual G4double DistanceToApplyOperation(const G4Track *, G4double, G4ForceCondition *) { return DBL_MAX; } + virtual G4double DistanceToApplyOperation(const G4Track*, G4double, G4ForceCondition*) { return DBL_MAX; } - virtual G4VParticleChange *GenerateBiasingFinalState(const G4Track *, const G4Step *) { return 0; } + virtual G4VParticleChange* GenerateBiasingFinalState(const G4Track*, const G4Step*) { return 0; } -public: + public: void SetSplittingFactor(G4int splittingFactor) { fSplittingFactor = splittingFactor; } G4int GetSplittingFactor() const { return fSplittingFactor; } -private: + private: G4int fSplittingFactor; G4ParticleChange fParticleChange; TVector3 fBiasingCenter; }; - class GammaBiasingOperator : public G4VBiasingOperator { -public: - GammaBiasingOperator(int splittingFactor, bool biasOnlyOnce, const TVector3 &biasingCenter); + public: + GammaBiasingOperator(int splittingFactor, bool biasOnlyOnce, const TVector3& biasingCenter); virtual ~GammaBiasingOperator() {} -public: + public: virtual void StartRun(); - virtual void StartTracking(const G4Track *track); + virtual void StartTracking(const G4Track* track); TVector3 GetBiasingCenter() const { return fBiasingCenter; } -private: - virtual G4VBiasingOperation *ProposeNonPhysicsBiasingOperation(const G4Track *, - const G4BiasingProcessInterface *) { + private: + virtual G4VBiasingOperation* ProposeNonPhysicsBiasingOperation(const G4Track*, + const G4BiasingProcessInterface*) { return 0; } - virtual G4VBiasingOperation *ProposeOccurenceBiasingOperation(const G4Track *, - const G4BiasingProcessInterface *) { + virtual G4VBiasingOperation* ProposeOccurenceBiasingOperation(const G4Track*, + const G4BiasingProcessInterface*) { return 0; } - virtual G4VBiasingOperation *ProposeFinalStateBiasingOperation( - const G4Track *track, const G4BiasingProcessInterface *callingProcess); + virtual G4VBiasingOperation* ProposeFinalStateBiasingOperation( + const G4Track* track, const G4BiasingProcessInterface* callingProcess); -private: + private: using G4VBiasingOperator::OperationApplied; -private: - GammaBiasingOperation *fBremSplittingOperation; + private: + GammaBiasingOperation* fBremSplittingOperation; G4int fSplittingFactor; G4bool fBiasOnlyOnce; G4int fNInteractions; TVector3 fBiasingCenter; }; -#endif //REST_BIASING_H +#endif // REST_BIASING_H diff --git a/src/Biasing.cxx b/src/Biasing.cxx index 911bad1a..2f70c74f 100644 --- a/src/Biasing.cxx +++ b/src/Biasing.cxx @@ -50,7 +50,7 @@ G4VParticleChange* GammaBiasingOperation::ApplyFinalStateBiasing( // random number between 0 and 1 const double rand = G4UniformRand(); // if random number is less than 1 / fSplittingFactor, keep alive - weightFactor = fSplittingFactor; // increase weight + weightFactor = fSplittingFactor; // increase weight if (rand < 1.0 / fSplittingFactor) { nSecondaries = 1; } else { @@ -94,21 +94,21 @@ G4VParticleChange* GammaBiasingOperation::ApplyFinalStateBiasing( #include "G4BiasingProcessInterface.hh" GammaBiasingOperator::GammaBiasingOperator(int splittingFactor, bool biasOnlyOnce, - const TVector3 &biasingCenter) - : G4VBiasingOperator("BremSplittingOperator"), - fSplittingFactor(splittingFactor), - fBiasOnlyOnce(biasOnlyOnce), - fBiasingCenter(biasingCenter) { + const TVector3& biasingCenter) + : G4VBiasingOperator("BremSplittingOperator"), + fSplittingFactor(splittingFactor), + fBiasOnlyOnce(biasOnlyOnce), + fBiasingCenter(biasingCenter) { fBremSplittingOperation = - new GammaBiasingOperation("BremSplittingOperation", splittingFactor, biasingCenter); + new GammaBiasingOperation("BremSplittingOperation", splittingFactor, biasingCenter); } void GammaBiasingOperator::StartRun() { fBremSplittingOperation->SetSplittingFactor(fSplittingFactor); } -void GammaBiasingOperator::StartTracking(const G4Track *) { fNInteractions = 0; } +void GammaBiasingOperator::StartTracking(const G4Track*) { fNInteractions = 0; } -G4VBiasingOperation *GammaBiasingOperator::ProposeFinalStateBiasingOperation( - const G4Track *, const G4BiasingProcessInterface *) { +G4VBiasingOperation* GammaBiasingOperator::ProposeFinalStateBiasingOperation( + const G4Track*, const G4BiasingProcessInterface*) { if (fBiasOnlyOnce && (fNInteractions > 0)) return nullptr; fNInteractions++; diff --git a/src/DetectorConstruction.cxx b/src/DetectorConstruction.cxx index 7d8bcf99..3591361d 100644 --- a/src/DetectorConstruction.cxx +++ b/src/DetectorConstruction.cxx @@ -23,16 +23,16 @@ using namespace std; using namespace TRestGeant4PrimaryGeneratorTypes; -DetectorConstruction::DetectorConstruction(SimulationManager *simulationManager) - : fSimulationManager(simulationManager) { +DetectorConstruction::DetectorConstruction(SimulationManager* simulationManager) + : fSimulationManager(simulationManager) { G4cout << "Detector Construction" << G4endl; fGdmlParser = new G4GDMLParser(); } DetectorConstruction::~DetectorConstruction() { delete fGdmlParser; } -G4VPhysicalVolume *DetectorConstruction::Construct() { - TRestGeant4Metadata *restG4Metadata = fSimulationManager->GetRestMetadata(); +G4VPhysicalVolume* DetectorConstruction::Construct() { + TRestGeant4Metadata* restG4Metadata = fSimulationManager->GetRestMetadata(); cout << "Isotope table " << endl; cout << *(G4Isotope::GetIsotopeTable()) << endl; @@ -45,31 +45,31 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { const auto startingPath = filesystem::current_path(); const auto [gdmlPath, gdmlToRead] = - TRestTools::SeparatePathAndName((string) restG4Metadata->GetGdmlFilename()); + TRestTools::SeparatePathAndName((string)restG4Metadata->GetGdmlFilename()); filesystem::current_path(gdmlPath); G4cout << "gdmlToRead: " << gdmlToRead << G4endl; fGdmlParser->Read(gdmlToRead, false); - G4VPhysicalVolume *worldVolume = fGdmlParser->GetWorldVolume(); + G4VPhysicalVolume* worldVolume = fGdmlParser->GetWorldVolume(); - const auto worldSolid = dynamic_cast(worldVolume->GetLogicalVolume()->GetSolid()); + const auto worldSolid = dynamic_cast(worldVolume->GetLogicalVolume()->GetSolid()); restG4Metadata->fGeant4PrimaryGeneratorInfo.fSpatialGeneratorWorldSize = { - worldSolid->GetXHalfLength(), worldSolid->GetYHalfLength(), worldSolid->GetZHalfLength()}; + worldSolid->GetXHalfLength(), worldSolid->GetYHalfLength(), worldSolid->GetZHalfLength()}; restG4Metadata->fGeant4GeometryInfo.InitializeOnDetectorConstruction(gdmlToRead, worldVolume); restG4Metadata->ReadDetector(); restG4Metadata->PrintMetadata(); // now we have detector info - const auto &geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); + const auto& geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); geometryInfo.Print(); // do some checks { // Check all physical volume names are unique - G4PhysicalVolumeStore *physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); + G4PhysicalVolumeStore* physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); set physicalVolumeNames; - vector::const_iterator physicalVolume; + vector::const_iterator physicalVolume; for (physicalVolume = physicalVolumeStore->begin(); physicalVolume != physicalVolumeStore->end(); physicalVolume++) { auto name = (*physicalVolume)->GetName(); @@ -84,9 +84,9 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { } // Check all logical volume names are unique - G4LogicalVolumeStore *logicalVolumeStore = G4LogicalVolumeStore::GetInstance(); + G4LogicalVolumeStore* logicalVolumeStore = G4LogicalVolumeStore::GetInstance(); set logicalVolumeNames; - vector::const_iterator logicalVolume; + vector::const_iterator logicalVolume; for (logicalVolume = logicalVolumeStore->begin(); logicalVolume != logicalVolumeStore->end(); logicalVolume++) { auto name = (*logicalVolume)->GetName(); @@ -102,15 +102,15 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { } filesystem::current_path(startingPath); - auto sensitiveVolume = (string) restG4Metadata->GetSensitiveVolume(); - G4VPhysicalVolume *physicalVolume = GetPhysicalVolume(sensitiveVolume); + auto sensitiveVolume = (string)restG4Metadata->GetSensitiveVolume(); + G4VPhysicalVolume* physicalVolume = GetPhysicalVolume(sensitiveVolume); if (physicalVolume == nullptr) { // sensitive volume was not found, perhaps the user specified a logical volume auto physicalVolumes = geometryInfo.GetAllPhysicalVolumesFromLogical(sensitiveVolume); if (physicalVolumes.size() == 1) { restG4Metadata->InsertSensitiveVolume( - geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolumes[0])); - sensitiveVolume = (string) restG4Metadata->GetSensitiveVolume(); + geometryInfo.GetAlternativeNameFromGeant4PhysicalName(physicalVolumes[0])); + sensitiveVolume = (string)restG4Metadata->GetSensitiveVolume(); physicalVolume = GetPhysicalVolume(sensitiveVolume); } } @@ -124,20 +124,20 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { Double_t my = restG4Metadata->GetMagneticField().Y() * tesla; Double_t mz = restG4Metadata->GetMagneticField().Z() * tesla; - G4MagneticField *magField = new G4UniformMagField(G4ThreeVector(mx, my, mz)); + G4MagneticField* magField = new G4UniformMagField(G4ThreeVector(mx, my, mz)); // G4FieldManager* localFieldMgr = new G4FieldManager(magField); - G4FieldManager *fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager(); + G4FieldManager* fieldMgr = G4TransportationManager::GetTransportationManager()->GetFieldManager(); fieldMgr->SetDetectorField(magField); fieldMgr->CreateChordFinder(magField); - G4LogicalVolume *volume = physicalVolume->GetLogicalVolume(); - G4Material *material = volume->GetMaterial(); + G4LogicalVolume* volume = physicalVolume->GetLogicalVolume(); + G4Material* material = volume->GetMaterial(); G4cout << "Sensitive volume properties:" << G4endl; G4cout << "\t- Material: " << material->GetName() << G4endl; G4cout << "\t- Temperature: " << material->GetTemperature() << " K" << G4endl; G4cout << "\t- Density: " << material->GetDensity() / (g / cm3) << " g/cm3" << G4endl; - const auto &primaryGeneratorInfo = restG4Metadata->GetGeant4PrimaryGeneratorInfo(); + const auto& primaryGeneratorInfo = restG4Metadata->GetGeant4PrimaryGeneratorInfo(); // Getting generation volume const auto fromVolume = primaryGeneratorInfo.GetSpatialGeneratorFrom(); if (fromVolume != "NO_SUCH_PARA") { @@ -146,11 +146,11 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { cout << "Generator type: " << primaryGeneratorInfo.GetSpatialGeneratorType() << endl; const auto spatialGeneratorTypeEnum = - StringToSpatialGeneratorTypes(primaryGeneratorInfo.GetSpatialGeneratorType().Data()); + StringToSpatialGeneratorTypes(primaryGeneratorInfo.GetSpatialGeneratorType().Data()); if (spatialGeneratorTypeEnum == TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::VOLUME && primaryGeneratorInfo.GetSpatialGeneratorFrom() != "Not defined") { - G4VPhysicalVolume *pVol = GetPhysicalVolume(primaryGeneratorInfo.GetSpatialGeneratorFrom().Data()); + G4VPhysicalVolume* pVol = GetPhysicalVolume(primaryGeneratorInfo.GetSpatialGeneratorFrom().Data()); if (pVol == nullptr) { cout << "ERROR: The generator volume '" << primaryGeneratorInfo.GetSpatialGeneratorFrom() << "'was not found in the geometry" << endl; @@ -162,7 +162,7 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { if (spatialGeneratorTypeEnum == TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::SURFACE || spatialGeneratorTypeEnum == TRestGeant4PrimaryGeneratorTypes::SpatialGeneratorTypes::VOLUME) { restG4Metadata->fGeant4PrimaryGeneratorInfo.fSpatialGeneratorPosition = { - fGeneratorTranslation.x(), fGeneratorTranslation.y(), fGeneratorTranslation.z()}; + fGeneratorTranslation.x(), fGeneratorTranslation.y(), fGeneratorTranslation.z()}; } fGeneratorSolid = pVol->GetLogicalVolume()->GetSolid(); @@ -203,9 +203,9 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { for (unsigned int id = 0; id < restG4Metadata->GetNumberOfActiveVolumes(); id++) { TString activeVolumeName = restG4Metadata->GetActiveVolumeName(id); - G4VPhysicalVolume *pVol = GetPhysicalVolume((G4String) activeVolumeName); + G4VPhysicalVolume* pVol = GetPhysicalVolume((G4String)activeVolumeName); if (pVol != nullptr) { - G4LogicalVolume *lVol = pVol->GetLogicalVolume(); + G4LogicalVolume* lVol = pVol->GetLogicalVolume(); if (restG4Metadata->GetMaxStepSize(activeVolumeName) > 0) { RESTInfo << "Setting maxStepSize of " << restG4Metadata->GetMaxStepSize(activeVolumeName) * mm << "mm for volume '" << activeVolumeName << "'" << RESTendl; @@ -221,15 +221,15 @@ G4VPhysicalVolume *DetectorConstruction::Construct() { return worldVolume; } -G4VPhysicalVolume *DetectorConstruction::GetPhysicalVolume(const G4String &physicalVolumeName) const { - G4PhysicalVolumeStore *physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); - TRestGeant4Metadata *restG4Metadata = fSimulationManager->GetRestMetadata(); - const auto &geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); - vector::const_iterator physicalVolume; +G4VPhysicalVolume* DetectorConstruction::GetPhysicalVolume(const G4String& physicalVolumeName) const { + G4PhysicalVolumeStore* physicalVolumeStore = G4PhysicalVolumeStore::GetInstance(); + TRestGeant4Metadata* restG4Metadata = fSimulationManager->GetRestMetadata(); + const auto& geometryInfo = restG4Metadata->GetGeant4GeometryInfo(); + vector::const_iterator physicalVolume; for (physicalVolume = physicalVolumeStore->begin(); physicalVolume != physicalVolumeStore->end(); physicalVolume++) { auto name = (*physicalVolume)->GetName(); - auto alternativeName = (G4String) geometryInfo.GetAlternativeNameFromGeant4PhysicalName(name); + auto alternativeName = (G4String)geometryInfo.GetAlternativeNameFromGeant4PhysicalName(name); if (name == physicalVolumeName || alternativeName == physicalVolumeName) { return *physicalVolume; } @@ -239,20 +239,20 @@ G4VPhysicalVolume *DetectorConstruction::GetPhysicalVolume(const G4String &physi } void DetectorConstruction::ConstructSDandField() { - const TRestGeant4Metadata &metadata = *fSimulationManager->GetRestMetadata(); + const TRestGeant4Metadata& metadata = *fSimulationManager->GetRestMetadata(); - set logicalVolumesSelected; - for (const auto &userSensitiveVolume: metadata.GetSensitiveVolumes()) { + set logicalVolumesSelected; + for (const auto& userSensitiveVolume : metadata.GetSensitiveVolumes()) { // Each sensitive detector declaration in the RML should correspond to at least one logical volume - G4LogicalVolume *logicalVolume; + G4LogicalVolume* logicalVolume; // Check if user selected a Geant4 physical volume by name - G4VPhysicalVolume *physicalVolume = - G4PhysicalVolumeStore::GetInstance()->GetVolume(userSensitiveVolume.Data(), false); + G4VPhysicalVolume* physicalVolume = + G4PhysicalVolumeStore::GetInstance()->GetVolume(userSensitiveVolume.Data(), false); if (physicalVolume == nullptr) { const G4String geant4VolumeName = - metadata.GetGeant4GeometryInfo() - .GetGeant4PhysicalNameFromAlternativeName(userSensitiveVolume.Data()) - .Data(); + metadata.GetGeant4GeometryInfo() + .GetGeant4PhysicalNameFromAlternativeName(userSensitiveVolume.Data()) + .Data(); // Check if user selected a Geant4 physical volume by REST name (only relevant for assemblies) physicalVolume = G4PhysicalVolumeStore::GetInstance()->GetVolume(geant4VolumeName, false); } @@ -269,25 +269,25 @@ void DetectorConstruction::ConstructSDandField() { // Check if the user string matches a logical volume by expanding the input as a regex // This can have multiple hits const auto logicalVolumesMatchingExpression = - metadata.GetGeant4GeometryInfo().GetAllLogicalVolumesMatchingExpression(userSensitiveVolume); + metadata.GetGeant4GeometryInfo().GetAllLogicalVolumesMatchingExpression(userSensitiveVolume); if (logicalVolumesMatchingExpression.empty()) { cerr << "Detector construction error: could not find matching logical volume(s) for '" << userSensitiveVolume << "'" << endl; exit(1); } else { - for (const auto &logicalVolumeName: logicalVolumesMatchingExpression) { + for (const auto& logicalVolumeName : logicalVolumesMatchingExpression) { logicalVolumesSelected.insert( - G4LogicalVolumeStore::GetInstance()->GetVolume(logicalVolumeName.Data(), false)); + G4LogicalVolumeStore::GetInstance()->GetVolume(logicalVolumeName.Data(), false)); } } } } - G4SDManager *SDManager = G4SDManager::GetSDMpointer(); + G4SDManager* SDManager = G4SDManager::GetSDMpointer(); - for (G4LogicalVolume *logicalVolume: logicalVolumesSelected) { + for (G4LogicalVolume* logicalVolume : logicalVolumesSelected) { auto name = logicalVolume->GetName(); - G4VSensitiveDetector *sensitiveDetector = new SensitiveDetector(fSimulationManager, name); + G4VSensitiveDetector* sensitiveDetector = new SensitiveDetector(fSimulationManager, name); SDManager->AddNewDetector(sensitiveDetector); logicalVolume->SetSensitiveDetector(sensitiveDetector); @@ -296,12 +296,12 @@ void DetectorConstruction::ConstructSDandField() { } // Biasing - const auto &biasingInfo = metadata.GetGeant4BiasingInfo(); + const auto& biasingInfo = metadata.GetGeant4BiasingInfo(); if (biasingInfo.GetEnabled()) { - const auto &geometryInfo = metadata.GetGeant4GeometryInfo(); + const auto& geometryInfo = metadata.GetGeant4GeometryInfo(); set biasingVolumes; - for (const auto &biasingVolume: biasingInfo.GetBiasingVolumes()) { + for (const auto& biasingVolume : biasingInfo.GetBiasingVolumes()) { auto name = biasingVolume; if (!geometryInfo.IsValidLogicalVolume(name)) { if (geometryInfo.IsValidPhysicalVolume(name)) { @@ -309,9 +309,9 @@ void DetectorConstruction::ConstructSDandField() { biasingVolumes.insert(name); } else { RESTError - << "volume name '" << name - << "' inside biasing section is invalid. Please check it belongs to a logical volume" - << RESTendl; + << "volume name '" << name + << "' inside biasing section is invalid. Please check it belongs to a logical volume" + << RESTendl; exit(1); } } else { @@ -319,9 +319,9 @@ void DetectorConstruction::ConstructSDandField() { } } - for (const auto &biasingLogicalVolumeName: biasingVolumes) { - G4LogicalVolume *biasingLogicalVolume = - G4LogicalVolumeStore::GetInstance()->GetVolume(biasingLogicalVolumeName, false); + for (const auto& biasingLogicalVolumeName : biasingVolumes) { + G4LogicalVolume* biasingLogicalVolume = + G4LogicalVolumeStore::GetInstance()->GetVolume(biasingLogicalVolumeName, false); if (biasingLogicalVolume == nullptr) { cerr << "Detector construction error: could not find logical volume '" << biasingLogicalVolumeName << "'" << endl; @@ -340,7 +340,7 @@ void DetectorConstruction::ConstructSDandField() { cout << "Biasing is enabled " << endl; cout << "Splitting factor: " << biasingInfo.GetSplittingFactor() << endl; cout << "Biasing volumes: " << endl; - for (const auto &volume: biasingVolumes) { + for (const auto& volume : biasingVolumes) { cout << "\t" << volume << endl; } const auto center = biasingInfo.GetBiasingCenter(); @@ -348,26 +348,26 @@ void DetectorConstruction::ConstructSDandField() { } } -void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume *world) { - const auto detector = (DetectorConstruction *) G4RunManager::GetRunManager()->GetUserDetectorConstruction(); - TRestGeant4Metadata *restG4Metadata = detector->fSimulationManager->GetRestMetadata(); +void TRestGeant4GeometryInfo::PopulateFromGeant4World(const G4VPhysicalVolume* world) { + const auto detector = (DetectorConstruction*)G4RunManager::GetRunManager()->GetUserDetectorConstruction(); + TRestGeant4Metadata* restG4Metadata = detector->fSimulationManager->GetRestMetadata(); const size_t n = int(world->GetLogicalVolume()->GetNoDaughters()); for (size_t i = 0; i < n + 1; i++) { // world is the + 1 - G4VPhysicalVolume *volume; + G4VPhysicalVolume* volume; if (i == n) { - volume = const_cast(world); + volume = const_cast(world); } else { volume = world->GetLogicalVolume()->GetDaughter(G4int(i)); } - TString namePhysical = (TString) volume->GetName(); + TString namePhysical = (TString)volume->GetName(); if (fGdmlNewPhysicalNames.size() > i) { // it has been filled fGeant4PhysicalNameToNewPhysicalNameMap[namePhysical] = fGdmlNewPhysicalNames[i]; } TString physicalNewName = GetAlternativeNameFromGeant4PhysicalName(namePhysical); - TString nameLogical = (TString) volume->GetLogicalVolume()->GetName(); - TString nameMaterial = (TString) volume->GetLogicalVolume()->GetMaterial()->GetName(); + TString nameLogical = (TString)volume->GetLogicalVolume()->GetName(); + TString nameMaterial = (TString)volume->GetLogicalVolume()->GetMaterial()->GetName(); auto position = volume->GetTranslation(); fPhysicalToLogicalVolumeMap[physicalNewName] = nameLogical; From ae791225e2cb5926ce658e0baf4ceaea392a2448 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Mon, 27 Feb 2023 12:16:21 +0100 Subject: [PATCH 20/24] add generic "cosmic" distribution (power -2.7) --- examples/13.IAXO/GammaBiasing.rml | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/examples/13.IAXO/GammaBiasing.rml b/examples/13.IAXO/GammaBiasing.rml index 347b1b31..4e839490 100644 --- a/examples/13.IAXO/GammaBiasing.rml +++ b/examples/13.IAXO/GammaBiasing.rml @@ -15,10 +15,10 @@ Author: Luis Obis (lobis@unizar.es) - + - - + + @@ -27,8 +27,8 @@ Author: Luis Obis (lobis@unizar.es) - - + + From f687f5832906f08359cac37fb57470cd1aaf1f92 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Mon, 27 Feb 2023 17:12:55 +0100 Subject: [PATCH 21/24] add option to save kill volume energy --- examples/15.GammaBiasing/ShieldingBiasing.rml | 41 +++++++++++ examples/15.GammaBiasing/geometry/setup.gdml | 70 ++++++++++++++++++ src/DataModel.cxx | 71 ++++++++++--------- 3 files changed, 147 insertions(+), 35 deletions(-) create mode 100644 examples/15.GammaBiasing/ShieldingBiasing.rml create mode 100644 examples/15.GammaBiasing/geometry/setup.gdml diff --git a/examples/15.GammaBiasing/ShieldingBiasing.rml b/examples/15.GammaBiasing/ShieldingBiasing.rml new file mode 100644 index 00000000..872129d4 --- /dev/null +++ b/examples/15.GammaBiasing/ShieldingBiasing.rml @@ -0,0 +1,41 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/15.GammaBiasing/geometry/setup.gdml b/examples/15.GammaBiasing/geometry/setup.gdml new file mode 100644 index 00000000..475af4cb --- /dev/null +++ b/examples/15.GammaBiasing/geometry/setup.gdml @@ -0,0 +1,70 @@ + + + + + ]> + + + + + + + + + + + + + + + + &materials; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/src/DataModel.cxx b/src/DataModel.cxx index 7a75fc22..b96a9902 100644 --- a/src/DataModel.cxx +++ b/src/DataModel.cxx @@ -10,21 +10,21 @@ using namespace std; -TRestGeant4Event::TRestGeant4Event(const G4Event* event) : TRestGeant4Event() { +TRestGeant4Event::TRestGeant4Event(const G4Event *event) : TRestGeant4Event() { SetID(event->GetEventID()); SetOK(true); time_t system_time = time(nullptr); - SetTime((Double_t)system_time); + SetTime((Double_t) system_time); auto primaryVertex = event->GetPrimaryVertex(); - const auto& position = primaryVertex->GetPosition(); + const auto &position = primaryVertex->GetPosition(); fPrimaryPosition = {position.x() / CLHEP::mm, position.y() / CLHEP::mm, position.z() / CLHEP::mm}; for (int i = 0; i < primaryVertex->GetNumberOfParticle(); i++) { - const auto& primaryParticle = primaryVertex->GetPrimary(i); + const auto &primaryParticle = primaryVertex->GetPrimary(i); fPrimaryParticleNames.emplace_back(primaryParticle->GetParticleDefinition()->GetParticleName()); fPrimaryEnergies.emplace_back(primaryParticle->GetKineticEnergy() / CLHEP::keV); - const auto& momentum = primaryParticle->GetMomentumDirection(); + const auto &momentum = primaryParticle->GetMomentumDirection(); fPrimaryDirections.emplace_back(momentum.x(), momentum.y(), momentum.z()); } @@ -49,7 +49,7 @@ TRestGeant4Event::TRestGeant4Event(const G4Event* event) : TRestGeant4Event() { */ } -bool TRestGeant4Event::InsertTrack(const G4Track* track) { +bool TRestGeant4Event::InsertTrack(const G4Track *track) { if (fInitialStep.GetNumberOfHits() != 1) { cout << "fInitialStep does not have exactly one step! Problem with stepping verbose" << endl; exit(1); @@ -60,10 +60,10 @@ bool TRestGeant4Event::InsertTrack(const G4Track* track) { // First track of sub-event (primary) fSubEventPrimaryParticleName = track->GetParticleDefinition()->GetParticleName(); fSubEventPrimaryEnergy = track->GetKineticEnergy() / CLHEP::keV; - const auto& position = track->GetPosition(); + const auto &position = track->GetPosition(); fSubEventPrimaryPosition = {position.x() / CLHEP::mm, position.y() / CLHEP::mm, position.z() / CLHEP::mm}; - const auto& momentum = track->GetMomentumDirection(); + const auto &momentum = track->GetMomentumDirection(); fSubEventPrimaryDirection = {momentum.x(), momentum.y(), momentum.z()}; } @@ -71,12 +71,12 @@ bool TRestGeant4Event::InsertTrack(const G4Track* track) { fTracks.emplace_back(track); - auto& insertedTrack = fTracks.back(); + auto &insertedTrack = fTracks.back(); insertedTrack.SetHits(fInitialStep); insertedTrack.SetEvent(this); - TRestGeant4Track* parentTrack = GetTrackByID(track->GetParentID()); + TRestGeant4Track *parentTrack = GetTrackByID(track->GetParentID()); if (parentTrack) { parentTrack->AddSecondaryTrackID(track->GetTrackID()); } @@ -84,9 +84,9 @@ bool TRestGeant4Event::InsertTrack(const G4Track* track) { return true; } -void TRestGeant4Event::UpdateTrack(const G4Track* track) { fTracks.back().UpdateTrack(track); } +void TRestGeant4Event::UpdateTrack(const G4Track *track) { fTracks.back().UpdateTrack(track); } -void TRestGeant4Event::InsertStep(const G4Step* step) { +void TRestGeant4Event::InsertStep(const G4Step *step) { if (step->GetTrack()->GetCurrentStepNumber() == 0) { // initial step (from SteppingVerbose) is generated before TrackingAction can insert the first track fInitialStep = TRestGeant4Hits(); @@ -97,11 +97,11 @@ void TRestGeant4Event::InsertStep(const G4Step* step) { } } -bool OutputManager::IsValidTrack(const G4Track*) const { return true; } +bool OutputManager::IsValidTrack(const G4Track *) const { return true; } -bool OutputManager::IsValidStep(const G4Step*) const { return true; } +bool OutputManager::IsValidStep(const G4Step *) const { return true; } -TRestGeant4Track::TRestGeant4Track(const G4Track* track) : TRestGeant4Track() { +TRestGeant4Track::TRestGeant4Track(const G4Track *track) : TRestGeant4Track() { fTrackID = track->GetTrackID(); fParentID = track->GetParentID(); @@ -124,13 +124,13 @@ TRestGeant4Track::TRestGeant4Track(const G4Track* track) : TRestGeant4Track() { fGlobalTimestamp = track->GetGlobalTime() / CLHEP::microsecond; - const G4ThreeVector& trackOrigin = track->GetPosition(); + const G4ThreeVector &trackOrigin = track->GetPosition(); fInitialPosition = {trackOrigin.x(), trackOrigin.y(), trackOrigin.z()}; } -void TRestGeant4Track::InsertStep(const G4Step* step) { fHits.InsertStep(step); } +void TRestGeant4Track::InsertStep(const G4Step *step) { fHits.InsertStep(step); } -void TRestGeant4Track::UpdateTrack(const G4Track* track) { +void TRestGeant4Track::UpdateTrack(const G4Track *track) { if (track->GetTrackID() != fTrackID) { G4cout << "Geant4Track::UpdateTrack - mismatch of trackID!" << endl; exit(1); @@ -140,19 +140,19 @@ void TRestGeant4Track::UpdateTrack(const G4Track* track) { fTimeLength = track->GetGlobalTime() / CLHEP::microsecond - fGlobalTimestamp; } -Int_t TRestGeant4PhysicsInfo::GetProcessIDFromGeant4Process(const G4VProcess* process) { +Int_t TRestGeant4PhysicsInfo::GetProcessIDFromGeant4Process(const G4VProcess *process) { return process->GetProcessType() * 1000 + process->GetProcessSubType(); } -void TRestGeant4Hits::InsertStep(const G4Step* step) { - const G4Track* track = step->GetTrack(); +void TRestGeant4Hits::InsertStep(const G4Step *step) { + const G4Track *track = step->GetTrack(); - TRestGeant4Metadata* metadata = GetGeant4Metadata(); + TRestGeant4Metadata *metadata = GetGeant4Metadata(); - const auto& geometryInfo = metadata->GetGeant4GeometryInfo(); + const auto &geometryInfo = metadata->GetGeant4GeometryInfo(); - const auto& volumeNameGeant4 = step->GetPreStepPoint()->GetPhysicalVolume()->GetName(); - const auto& volumeName = geometryInfo.GetAlternativeNameFromGeant4PhysicalName(volumeNameGeant4); + const auto &volumeNameGeant4 = step->GetPreStepPoint()->GetPhysicalVolume()->GetName(); + const auto &volumeName = geometryInfo.GetAlternativeNameFromGeant4PhysicalName(volumeNameGeant4); if (!metadata->IsActiveVolume(volumeName) && step->GetTrack()->GetCurrentStepNumber() != 0) { // we always store the first step @@ -161,9 +161,9 @@ void TRestGeant4Hits::InsertStep(const G4Step* step) { const bool kill = metadata->IsKillVolume(volumeName); - const auto& particle = step->GetTrack()->GetDefinition(); - const auto& particleID = particle->GetPDGEncoding(); - const auto& particleName = particle->GetParticleName(); + const auto &particle = step->GetTrack()->GetDefinition(); + const auto &particleID = particle->GetPDGEncoding(); + const auto &particleName = particle->GetParticleName(); auto energy = step->GetTotalEnergyDeposit() / CLHEP::keV; @@ -195,7 +195,8 @@ void TRestGeant4Hits::InsertStep(const G4Step* step) { processName = "REST-for-physics-kill"; processTypeName = "REST-for-physics"; processID = 1000000; // use id out of range! - energy = 0; + const bool computeEnergy = metadata->GetKillVolumesComputeEnergy(); + energy = computeEnergy ? step->GetTrack()->GetKineticEnergy() / CLHEP::keV : 0; step->GetTrack()->SetTrackStatus(fStopAndKill); } @@ -203,9 +204,9 @@ void TRestGeant4Hits::InsertStep(const G4Step* step) { metadata->fGeant4PhysicsInfo.InsertProcessName(processID, processName, processTypeName); auto sensitiveVolumeName = - geometryInfo.GetAlternativeNameFromGeant4PhysicalName(metadata->GetSensitiveVolume()); + geometryInfo.GetAlternativeNameFromGeant4PhysicalName(metadata->GetSensitiveVolume()); - G4Track* aTrack = step->GetTrack(); + G4Track *aTrack = step->GetTrack(); Double_t x = aTrack->GetPosition().x() / CLHEP::mm; Double_t y = aTrack->GetPosition().y() / CLHEP::mm; @@ -213,7 +214,7 @@ void TRestGeant4Hits::InsertStep(const G4Step* step) { const TVector3 hitPosition(x, y, z); const Double_t hitGlobalTime = track->GetGlobalTime() / CLHEP::microsecond; - const G4ThreeVector& momentum = track->GetMomentumDirection(); + const G4ThreeVector &momentum = track->GetMomentumDirection(); AddHit(hitPosition, energy, hitGlobalTime); // this increases fNHits @@ -227,9 +228,9 @@ void TRestGeant4Hits::InsertStep(const G4Step* step) { } void OutputManager::RemoveUnwantedTracks() { - const auto& metadata = fSimulationManager->GetRestMetadata(); + const auto &metadata = fSimulationManager->GetRestMetadata(); set trackIDsToKeep; // We populate this container with the tracks we want to keep - for (const auto& track : fEvent->fTracks) { + for (const auto &track: fEvent->fTracks) { // If one children track is kept, we keep all the parents if (trackIDsToKeep.count(track.GetTrackID()) > 0) { continue; @@ -255,7 +256,7 @@ void OutputManager::RemoveUnwantedTracks() { // const size_t numberOfTracksBefore = fEvent->fTracks.size(); vector tracksAfterRemoval; - for (const auto& track : fEvent->fTracks) { + for (const auto &track: fEvent->fTracks) { // we do this to preserve original order if (trackIDsToKeep.count(track.GetTrackID()) > 0) { tracksAfterRemoval.push_back(*(fEvent->GetTrackByID(track.GetTrackID()))); From 3b897b11d7cb27e8c62e11b9f13b158380e13572 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Mon, 27 Feb 2023 16:13:12 +0000 Subject: [PATCH 22/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- examples/15.GammaBiasing/geometry/setup.gdml | 2 +- src/DataModel.cxx | 68 ++++++++++---------- 2 files changed, 35 insertions(+), 35 deletions(-) diff --git a/examples/15.GammaBiasing/geometry/setup.gdml b/examples/15.GammaBiasing/geometry/setup.gdml index 475af4cb..05f05a19 100644 --- a/examples/15.GammaBiasing/geometry/setup.gdml +++ b/examples/15.GammaBiasing/geometry/setup.gdml @@ -62,7 +62,7 @@ - + diff --git a/src/DataModel.cxx b/src/DataModel.cxx index b96a9902..3c1a585c 100644 --- a/src/DataModel.cxx +++ b/src/DataModel.cxx @@ -10,21 +10,21 @@ using namespace std; -TRestGeant4Event::TRestGeant4Event(const G4Event *event) : TRestGeant4Event() { +TRestGeant4Event::TRestGeant4Event(const G4Event* event) : TRestGeant4Event() { SetID(event->GetEventID()); SetOK(true); time_t system_time = time(nullptr); - SetTime((Double_t) system_time); + SetTime((Double_t)system_time); auto primaryVertex = event->GetPrimaryVertex(); - const auto &position = primaryVertex->GetPosition(); + const auto& position = primaryVertex->GetPosition(); fPrimaryPosition = {position.x() / CLHEP::mm, position.y() / CLHEP::mm, position.z() / CLHEP::mm}; for (int i = 0; i < primaryVertex->GetNumberOfParticle(); i++) { - const auto &primaryParticle = primaryVertex->GetPrimary(i); + const auto& primaryParticle = primaryVertex->GetPrimary(i); fPrimaryParticleNames.emplace_back(primaryParticle->GetParticleDefinition()->GetParticleName()); fPrimaryEnergies.emplace_back(primaryParticle->GetKineticEnergy() / CLHEP::keV); - const auto &momentum = primaryParticle->GetMomentumDirection(); + const auto& momentum = primaryParticle->GetMomentumDirection(); fPrimaryDirections.emplace_back(momentum.x(), momentum.y(), momentum.z()); } @@ -49,7 +49,7 @@ TRestGeant4Event::TRestGeant4Event(const G4Event *event) : TRestGeant4Event() { */ } -bool TRestGeant4Event::InsertTrack(const G4Track *track) { +bool TRestGeant4Event::InsertTrack(const G4Track* track) { if (fInitialStep.GetNumberOfHits() != 1) { cout << "fInitialStep does not have exactly one step! Problem with stepping verbose" << endl; exit(1); @@ -60,10 +60,10 @@ bool TRestGeant4Event::InsertTrack(const G4Track *track) { // First track of sub-event (primary) fSubEventPrimaryParticleName = track->GetParticleDefinition()->GetParticleName(); fSubEventPrimaryEnergy = track->GetKineticEnergy() / CLHEP::keV; - const auto &position = track->GetPosition(); + const auto& position = track->GetPosition(); fSubEventPrimaryPosition = {position.x() / CLHEP::mm, position.y() / CLHEP::mm, position.z() / CLHEP::mm}; - const auto &momentum = track->GetMomentumDirection(); + const auto& momentum = track->GetMomentumDirection(); fSubEventPrimaryDirection = {momentum.x(), momentum.y(), momentum.z()}; } @@ -71,12 +71,12 @@ bool TRestGeant4Event::InsertTrack(const G4Track *track) { fTracks.emplace_back(track); - auto &insertedTrack = fTracks.back(); + auto& insertedTrack = fTracks.back(); insertedTrack.SetHits(fInitialStep); insertedTrack.SetEvent(this); - TRestGeant4Track *parentTrack = GetTrackByID(track->GetParentID()); + TRestGeant4Track* parentTrack = GetTrackByID(track->GetParentID()); if (parentTrack) { parentTrack->AddSecondaryTrackID(track->GetTrackID()); } @@ -84,9 +84,9 @@ bool TRestGeant4Event::InsertTrack(const G4Track *track) { return true; } -void TRestGeant4Event::UpdateTrack(const G4Track *track) { fTracks.back().UpdateTrack(track); } +void TRestGeant4Event::UpdateTrack(const G4Track* track) { fTracks.back().UpdateTrack(track); } -void TRestGeant4Event::InsertStep(const G4Step *step) { +void TRestGeant4Event::InsertStep(const G4Step* step) { if (step->GetTrack()->GetCurrentStepNumber() == 0) { // initial step (from SteppingVerbose) is generated before TrackingAction can insert the first track fInitialStep = TRestGeant4Hits(); @@ -97,11 +97,11 @@ void TRestGeant4Event::InsertStep(const G4Step *step) { } } -bool OutputManager::IsValidTrack(const G4Track *) const { return true; } +bool OutputManager::IsValidTrack(const G4Track*) const { return true; } -bool OutputManager::IsValidStep(const G4Step *) const { return true; } +bool OutputManager::IsValidStep(const G4Step*) const { return true; } -TRestGeant4Track::TRestGeant4Track(const G4Track *track) : TRestGeant4Track() { +TRestGeant4Track::TRestGeant4Track(const G4Track* track) : TRestGeant4Track() { fTrackID = track->GetTrackID(); fParentID = track->GetParentID(); @@ -124,13 +124,13 @@ TRestGeant4Track::TRestGeant4Track(const G4Track *track) : TRestGeant4Track() { fGlobalTimestamp = track->GetGlobalTime() / CLHEP::microsecond; - const G4ThreeVector &trackOrigin = track->GetPosition(); + const G4ThreeVector& trackOrigin = track->GetPosition(); fInitialPosition = {trackOrigin.x(), trackOrigin.y(), trackOrigin.z()}; } -void TRestGeant4Track::InsertStep(const G4Step *step) { fHits.InsertStep(step); } +void TRestGeant4Track::InsertStep(const G4Step* step) { fHits.InsertStep(step); } -void TRestGeant4Track::UpdateTrack(const G4Track *track) { +void TRestGeant4Track::UpdateTrack(const G4Track* track) { if (track->GetTrackID() != fTrackID) { G4cout << "Geant4Track::UpdateTrack - mismatch of trackID!" << endl; exit(1); @@ -140,19 +140,19 @@ void TRestGeant4Track::UpdateTrack(const G4Track *track) { fTimeLength = track->GetGlobalTime() / CLHEP::microsecond - fGlobalTimestamp; } -Int_t TRestGeant4PhysicsInfo::GetProcessIDFromGeant4Process(const G4VProcess *process) { +Int_t TRestGeant4PhysicsInfo::GetProcessIDFromGeant4Process(const G4VProcess* process) { return process->GetProcessType() * 1000 + process->GetProcessSubType(); } -void TRestGeant4Hits::InsertStep(const G4Step *step) { - const G4Track *track = step->GetTrack(); +void TRestGeant4Hits::InsertStep(const G4Step* step) { + const G4Track* track = step->GetTrack(); - TRestGeant4Metadata *metadata = GetGeant4Metadata(); + TRestGeant4Metadata* metadata = GetGeant4Metadata(); - const auto &geometryInfo = metadata->GetGeant4GeometryInfo(); + const auto& geometryInfo = metadata->GetGeant4GeometryInfo(); - const auto &volumeNameGeant4 = step->GetPreStepPoint()->GetPhysicalVolume()->GetName(); - const auto &volumeName = geometryInfo.GetAlternativeNameFromGeant4PhysicalName(volumeNameGeant4); + const auto& volumeNameGeant4 = step->GetPreStepPoint()->GetPhysicalVolume()->GetName(); + const auto& volumeName = geometryInfo.GetAlternativeNameFromGeant4PhysicalName(volumeNameGeant4); if (!metadata->IsActiveVolume(volumeName) && step->GetTrack()->GetCurrentStepNumber() != 0) { // we always store the first step @@ -161,9 +161,9 @@ void TRestGeant4Hits::InsertStep(const G4Step *step) { const bool kill = metadata->IsKillVolume(volumeName); - const auto &particle = step->GetTrack()->GetDefinition(); - const auto &particleID = particle->GetPDGEncoding(); - const auto &particleName = particle->GetParticleName(); + const auto& particle = step->GetTrack()->GetDefinition(); + const auto& particleID = particle->GetPDGEncoding(); + const auto& particleName = particle->GetParticleName(); auto energy = step->GetTotalEnergyDeposit() / CLHEP::keV; @@ -204,9 +204,9 @@ void TRestGeant4Hits::InsertStep(const G4Step *step) { metadata->fGeant4PhysicsInfo.InsertProcessName(processID, processName, processTypeName); auto sensitiveVolumeName = - geometryInfo.GetAlternativeNameFromGeant4PhysicalName(metadata->GetSensitiveVolume()); + geometryInfo.GetAlternativeNameFromGeant4PhysicalName(metadata->GetSensitiveVolume()); - G4Track *aTrack = step->GetTrack(); + G4Track* aTrack = step->GetTrack(); Double_t x = aTrack->GetPosition().x() / CLHEP::mm; Double_t y = aTrack->GetPosition().y() / CLHEP::mm; @@ -214,7 +214,7 @@ void TRestGeant4Hits::InsertStep(const G4Step *step) { const TVector3 hitPosition(x, y, z); const Double_t hitGlobalTime = track->GetGlobalTime() / CLHEP::microsecond; - const G4ThreeVector &momentum = track->GetMomentumDirection(); + const G4ThreeVector& momentum = track->GetMomentumDirection(); AddHit(hitPosition, energy, hitGlobalTime); // this increases fNHits @@ -228,9 +228,9 @@ void TRestGeant4Hits::InsertStep(const G4Step *step) { } void OutputManager::RemoveUnwantedTracks() { - const auto &metadata = fSimulationManager->GetRestMetadata(); + const auto& metadata = fSimulationManager->GetRestMetadata(); set trackIDsToKeep; // We populate this container with the tracks we want to keep - for (const auto &track: fEvent->fTracks) { + for (const auto& track : fEvent->fTracks) { // If one children track is kept, we keep all the parents if (trackIDsToKeep.count(track.GetTrackID()) > 0) { continue; @@ -256,7 +256,7 @@ void OutputManager::RemoveUnwantedTracks() { // const size_t numberOfTracksBefore = fEvent->fTracks.size(); vector tracksAfterRemoval; - for (const auto &track: fEvent->fTracks) { + for (const auto& track : fEvent->fTracks) { // we do this to preserve original order if (trackIDsToKeep.count(track.GetTrackID()) > 0) { tracksAfterRemoval.push_back(*(fEvent->GetTrackByID(track.GetTrackID()))); From 7245b172c4db2c0c8cc67604a7ba6a41fafb2627 Mon Sep 17 00:00:00 2001 From: Luis Antonio Obis Aparicio Date: Tue, 28 Feb 2023 10:35:57 +0100 Subject: [PATCH 23/24] working biasing example --- examples/13.IAXO/GammaBiasing.rml | 1 - examples/15.GammaBiasing/ShieldingBiasing.rml | 68 +++++++++++++------ examples/15.GammaBiasing/geometry/setup.gdml | 11 ++- src/DataModel.cxx | 7 +- src/SensitiveDetector.cxx | 13 ++-- 5 files changed, 69 insertions(+), 31 deletions(-) diff --git a/examples/13.IAXO/GammaBiasing.rml b/examples/13.IAXO/GammaBiasing.rml index 4e839490..0300b94c 100644 --- a/examples/13.IAXO/GammaBiasing.rml +++ b/examples/13.IAXO/GammaBiasing.rml @@ -37,7 +37,6 @@ Author: Luis Obis (lobis@unizar.es) - l diff --git a/examples/15.GammaBiasing/ShieldingBiasing.rml b/examples/15.GammaBiasing/ShieldingBiasing.rml index 872129d4..4b5329bb 100644 --- a/examples/15.GammaBiasing/ShieldingBiasing.rml +++ b/examples/15.GammaBiasing/ShieldingBiasing.rml @@ -1,41 +1,69 @@ + + - - - - - - - - - - - - - - - + + - - + + + + + + + - + - + + + + + - + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/15.GammaBiasing/geometry/setup.gdml b/examples/15.GammaBiasing/geometry/setup.gdml index 475af4cb..f234dd89 100644 --- a/examples/15.GammaBiasing/geometry/setup.gdml +++ b/examples/15.GammaBiasing/geometry/setup.gdml @@ -20,9 +20,16 @@ - + + + + + + + + &materials; @@ -62,7 +69,7 @@ - + diff --git a/src/DataModel.cxx b/src/DataModel.cxx index b96a9902..e8ac2b3d 100644 --- a/src/DataModel.cxx +++ b/src/DataModel.cxx @@ -196,7 +196,12 @@ void TRestGeant4Hits::InsertStep(const G4Step *step) { processTypeName = "REST-for-physics"; processID = 1000000; // use id out of range! const bool computeEnergy = metadata->GetKillVolumesComputeEnergy(); - energy = computeEnergy ? step->GetTrack()->GetKineticEnergy() / CLHEP::keV : 0; + energy = 0; + if (computeEnergy) { + // needs to be here because track is killed before entering sensitive detector + energy = step->GetTrack()->GetKineticEnergy() / CLHEP::keV; + SimulationManager::GetOutputManager()->AddSensitiveEnergy(energy); + } step->GetTrack()->SetTrackStatus(fStopAndKill); } diff --git a/src/SensitiveDetector.cxx b/src/SensitiveDetector.cxx index 027ca9c3..20c0733b 100644 --- a/src/SensitiveDetector.cxx +++ b/src/SensitiveDetector.cxx @@ -4,22 +4,21 @@ #include #include #include -#include #include #include "SimulationManager.h" using namespace std; -SensitiveDetector::SensitiveDetector(SimulationManager* simulationManager, const G4String& name) - : G4VSensitiveDetector(name), fSimulationManager(simulationManager) {} +SensitiveDetector::SensitiveDetector(SimulationManager *simulationManager, const G4String &name) + : G4VSensitiveDetector(name), fSimulationManager(simulationManager) {} -G4bool SensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) { +G4bool SensitiveDetector::ProcessHits(G4Step *step, G4TouchableHistory *) { // return value will always be ignored, its present for backwards compatibility (I guess) const auto volumeName = fSimulationManager->GetRestMetadata() - ->GetGeant4GeometryInfo() - .GetAlternativeNameFromGeant4PhysicalName( - step->GetPreStepPoint()->GetPhysicalVolume()->GetName()); + ->GetGeant4GeometryInfo() + .GetAlternativeNameFromGeant4PhysicalName( + step->GetPreStepPoint()->GetPhysicalVolume()->GetName()); const bool isGeantino = step->GetTrack()->GetParticleDefinition() == G4Geantino::Definition(); From 9acbf0f53ddc426e0bcfae0f90e81a37be6e2084 Mon Sep 17 00:00:00 2001 From: "pre-commit-ci[bot]" <66853113+pre-commit-ci[bot]@users.noreply.github.com> Date: Tue, 28 Feb 2023 09:40:23 +0000 Subject: [PATCH 24/24] [pre-commit.ci] auto fixes from pre-commit.com hooks for more information, see https://pre-commit.ci --- src/SensitiveDetector.cxx | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/SensitiveDetector.cxx b/src/SensitiveDetector.cxx index 20c0733b..9c436d8e 100644 --- a/src/SensitiveDetector.cxx +++ b/src/SensitiveDetector.cxx @@ -10,15 +10,15 @@ using namespace std; -SensitiveDetector::SensitiveDetector(SimulationManager *simulationManager, const G4String &name) - : G4VSensitiveDetector(name), fSimulationManager(simulationManager) {} +SensitiveDetector::SensitiveDetector(SimulationManager* simulationManager, const G4String& name) + : G4VSensitiveDetector(name), fSimulationManager(simulationManager) {} -G4bool SensitiveDetector::ProcessHits(G4Step *step, G4TouchableHistory *) { +G4bool SensitiveDetector::ProcessHits(G4Step* step, G4TouchableHistory*) { // return value will always be ignored, its present for backwards compatibility (I guess) const auto volumeName = fSimulationManager->GetRestMetadata() - ->GetGeant4GeometryInfo() - .GetAlternativeNameFromGeant4PhysicalName( - step->GetPreStepPoint()->GetPhysicalVolume()->GetName()); + ->GetGeant4GeometryInfo() + .GetAlternativeNameFromGeant4PhysicalName( + step->GetPreStepPoint()->GetPhysicalVolume()->GetName()); const bool isGeantino = step->GetTrack()->GetParticleDefinition() == G4Geantino::Definition();