diff --git a/examples/13.IAXO/GammaBiasing.rml b/examples/13.IAXO/GammaBiasing.rml new file mode 100644 index 00000000..0300b94c --- /dev/null +++ b/examples/13.IAXO/GammaBiasing.rml @@ -0,0 +1,72 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + 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) - + diff --git a/examples/13.IAXO/geometry/setup.gdml b/examples/13.IAXO/geometry/setup.gdml index 4279c32c..7572fe61 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 @@ - + - + diff --git a/examples/15.GammaBiasing/ShieldingBiasing.rml b/examples/15.GammaBiasing/ShieldingBiasing.rml new file mode 100644 index 00000000..4b5329bb --- /dev/null +++ b/examples/15.GammaBiasing/ShieldingBiasing.rml @@ -0,0 +1,69 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/examples/15.GammaBiasing/geometry/setup.gdml b/examples/15.GammaBiasing/geometry/setup.gdml new file mode 100644 index 00000000..f234dd89 --- /dev/null +++ b/examples/15.GammaBiasing/geometry/setup.gdml @@ -0,0 +1,77 @@ + + + + + ]> + + + + + + + + + + + + + + + + + + + + + + + &materials; + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + diff --git a/include/Biasing.h b/include/Biasing.h new file mode 100644 index 00000000..3b5f3aa4 --- /dev/null +++ b/include/Biasing.h @@ -0,0 +1,80 @@ + +#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/PhysicsList.h b/include/PhysicsList.h index ee1cdb8b..0c275159 100644 --- a/include/PhysicsList.h +++ b/include/PhysicsList.h @@ -10,10 +10,14 @@ #include "SimulationManager.h" +class G4GenericBiasingPhysics; + class PhysicsList : public G4VModularPhysicsList { public: PhysicsList() = delete; + explicit PhysicsList(SimulationManager* simulationManager, TRestGeant4PhysicsLists* restPhysicsLists); + ~PhysicsList() override; protected: @@ -21,7 +25,9 @@ class PhysicsList : public G4VModularPhysicsList { virtual void InitializePhysicsLists(); void ConstructParticle() override; + void ConstructProcess() override; + void SetCuts() override; private: @@ -36,6 +42,8 @@ class PhysicsList : public G4VModularPhysicsList { G4VPhysicsConstructor* fRadDecPhysicsList = nullptr; std::vector fHadronPhys; + G4GenericBiasingPhysics* fBiasingPhysicsList = nullptr; + TRestGeant4PhysicsLists* fRestPhysicsLists = nullptr; }; 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/Biasing.cxx b/src/Biasing.cxx new file mode 100644 index 00000000..2f70c74f --- /dev/null +++ b/src/Biasing.cxx @@ -0,0 +1,117 @@ + +#include "Biasing.h" + +#include + +// Operation + +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); + + if (fSplittingFactor == 1) return processFinalState; + + // special case: no secondaries + if (processFinalState->GetNumberOfSecondaries() == 0) return 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.SetSecondaryWeightByProcess(true); + + 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) { + // 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 + weightFactor = fSplittingFactor; // increase weight + if (rand < 1.0 / fSplittingFactor) { + nSecondaries = 1; + } else { + nSecondaries = 0; + } + } + + G4double gammaWeight = track->GetWeight() * weightFactor; + // G4cout << "Gamma weight: " << gammaWeight << " nSecondaries: " << nSecondaries << G4endl; + gammaTrack->SetWeight(gammaWeight); + + fParticleChange.SetNumberOfSecondaries(nSecondaries); + if (nSecondaries > 0) { + fParticleChange.AddSecondary(gammaTrack); + } + + actualParticleChange->Clear(); + + G4int nCalls = 1; + while (nCalls < nSecondaries) { + 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(); + } + + 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/DataModel.cxx b/src/DataModel.cxx index da441ee4..d1f9f577 100644 --- a/src/DataModel.cxx +++ b/src/DataModel.cxx @@ -170,21 +170,38 @@ 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) { processName = "REST-for-physics-kill"; processTypeName = "REST-for-physics"; processID = 1000000; // use id out of range! + const bool computeEnergy = metadata->GetKillVolumesComputeEnergy(); 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/DetectorConstruction.cxx b/src/DetectorConstruction.cxx index e540e933..3591361d 100644 --- a/src/DetectorConstruction.cxx +++ b/src/DetectorConstruction.cxx @@ -17,6 +17,7 @@ #include #include +#include "Biasing.h" #include "SimulationManager.h" using namespace std; @@ -293,10 +294,62 @@ void DetectorConstruction::ConstructSDandField() { auto region = new G4Region(name); logicalVolume->SetRegion(region); } + + // Biasing + 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); + } + } else { + biasingVolumes.insert(name); + } + } + + 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(biasingInfo.GetSplittingFactor(), true, + biasingInfo.GetBiasingCenter()); + 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) { - 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 +358,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/PhysicsList.cxx b/src/PhysicsList.cxx index 55748595..f8051f9a 100644 --- a/src/PhysicsList.cxx +++ b/src/PhysicsList.cxx @@ -12,6 +12,7 @@ #include #include #include +#include #include #include #include @@ -128,6 +129,14 @@ void PhysicsList::InitializePhysicsLists() { emCounter++; } + 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) { RESTWarning << "PhysicsList: No EM physics list has been enabled" << RESTendl; @@ -179,6 +188,10 @@ void PhysicsList::ConstructParticle() { fRadDecPhysicsList->ConstructParticle(); } + if (fBiasingPhysicsList) { + fBiasingPhysicsList->ConstructParticle(); + } + for (auto& hadronicPhysicsList : fHadronPhys) { hadronicPhysicsList->ConstructParticle(); } @@ -229,6 +242,10 @@ void PhysicsList::ConstructProcess() { fRadDecPhysicsList->ConstructProcess(); } + if (fBiasingPhysicsList) { + fBiasingPhysicsList->ConstructProcess(); + } + // Hadronic physics lists for (auto& hadronicPhysicsList : fHadronPhys) { hadronicPhysicsList->ConstructProcess(); diff --git a/src/SensitiveDetector.cxx b/src/SensitiveDetector.cxx index 027ca9c3..9c436d8e 100644 --- a/src/SensitiveDetector.cxx +++ b/src/SensitiveDetector.cxx @@ -4,7 +4,6 @@ #include #include #include -#include #include #include "SimulationManager.h" diff --git a/src/TrackingAction.cxx b/src/TrackingAction.cxx index 93c5d963..47a7020a 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,9 +12,8 @@ using namespace std; TrackingAction::TrackingAction(SimulationManager* simulationManager) : G4UserTrackingAction(), fSimulationManager(simulationManager) {} -TrackingAction::~TrackingAction() {} - void TrackingAction::PreUserTrackingAction(const G4Track* track) { + // G4cout << track->GetVolume()->GetLogicalVolume()->GetName() << G4endl; fSimulationManager->GetOutputManager()->RecordTrack(track); }