GB06/src/GB06BOptnSplitAndKillByImportance.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 /// \file GB06BOptnSplitAndKillByImportance.cc
0028 /// \brief Implementation of the GB06BOptnSplitAndKillByImportance class
0029
0030 #include "GB06BOptnSplitAndKillByImportance.hh"
0031
0032 #include "G4BiasingProcessInterface.hh"
0033 #include "G4BiasingProcessSharedData.hh"
0034 #include "G4ParallelGeometriesLimiterProcess.hh"
0035 #include "G4ProcessManager.hh"
0036 #include "Randomize.hh"
0037
0038 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0039
0040 GB06BOptnSplitAndKillByImportance::GB06BOptnSplitAndKillByImportance(G4String name)
0041   : G4VBiasingOperation(name),
0042     fParallelWorldIndex(-1),
0043     fBiasingSharedData(nullptr),
0044     fParticleChange(),
0045     fDummyParticleChange()
0046 {}
0047
0048 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0049
0050 GB06BOptnSplitAndKillByImportance::~GB06BOptnSplitAndKillByImportance() {}
0051
0052 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0053
0054 G4double GB06BOptnSplitAndKillByImportance::DistanceToApplyOperation(const G4Track*, G4double,
0055                                                                      G4ForceCondition* condition)
0056 {
0057   // -- Remember the touchable history (ie geometry state) at the beginning of the step:
0058   // -- Start by getting the process handling the step limitation in parallel
0059   // -- geometries for the generic biasing:
0060   auto biasingLimiterProcess = fBiasingSharedData->GetParallelGeometriesLimiterProcess();
0061   fPreStepTouchableHistory =
0062     biasingLimiterProcess
0063       ->GetNavigator(fParallelWorldIndex)  // -- get the navigator of the geometry...
0064       ->CreateTouchableHistoryHandle();  // -- ...and collect the geometry state.
0065
0066   // -- We request to be "forced" : meaning GenerateBiasingFinalState(...) will be called
0067   // -- anyway at the PostStepDoIt(...) stage (ie, when the track will have been moved to
0068   // -- its end point position) and, there, we will have to handle the decision to
0069   // -- split/kill if we are on a volume boundary, or do nothing, if we're not:
0070   *condition = Forced;
0071
0072   // -- As we're forced, we make this returned length void:
0073   return DBL_MAX;
0074 }
0075
0076 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0077
0078 G4VParticleChange*
0079 GB06BOptnSplitAndKillByImportance::GenerateBiasingFinalState(const G4Track* track, const G4Step*)
0080 {
0081   // -- Given we used the "Forced" condition, this method is always called.
0082   // -- We check the status of the step in the parallel geometry, and apply
0083   // -- splitting/killing if the step has been limited by the geometry.
0084
0085   // -- We first check if we limit the step in the parallel geometry:
0086   G4bool isLimiting =
0087     fBiasingSharedData->GetParallelGeometriesLimiterProcess()->GetIsLimiting(fParallelWorldIndex);
0088
0089   // -- if not limiting, we leave the track unchanged:
0090   if (!isLimiting) {
0091     fDummyParticleChange.Initialize(*track);
0092     return &fDummyParticleChange;
0093   }
0094
0095   // -- We collect the geometry state at the end point step:
0096   // -- Note this is the same call than in the DistanceToApplyOperation(...) for the
0097   // -- fPreStepTouchableHistory, but the navigator has pushed the track in the next
0098   // -- volume since then (even if the track is still on the boundary), and hence the
0099   // -- geometry state has changed.
0100   auto biasingLimiterProcess = fBiasingSharedData->GetParallelGeometriesLimiterProcess();
0101   fPostStepTouchableHistory =
0102     biasingLimiterProcess->GetNavigator(fParallelWorldIndex)->CreateTouchableHistoryHandle();
0103
0104   // -- We verify we are still in the "same" physical volume:
0105   // -- remember : using a replica, we have "one" physical volume
0106   // -- but representing many different placements, distinguished
0107   // -- by replica number. By checking we are in the same physical
0108   // -- volume, we verify the end step point has not reached the
0109   // -- world volume of the parallel geometry:
0110   if (fPreStepTouchableHistory->GetVolume() != fPostStepTouchableHistory->GetVolume()) {
0111     // -- the track is leaving the volumes in which biasing is applied,
0112     // -- we leave this track unchanged:
0113     fDummyParticleChange.Initialize(*track);
0114     return &fDummyParticleChange;
0115   }
0116
0117   // -------------------------------------------------------------------------------------
0118   // -- At this stage, we know we have a particle crossing a boundary between two slices,
0119   // -- each of this slice has a well defined importance : we apply the biasing.
0120   // -- We will split if the track is entering a volume with higher importance, and
0121   // -- kill (applying Russian roulette) in the other case.
0122   // -------------------------------------------------------------------------------------
0123
0124   // -- We start by getting the replica numbers:
0125   G4int preReplicaNumber = fPreStepTouchableHistory->GetReplicaNumber();
0126   G4int postReplicaNumber = fPostStepTouchableHistory->GetReplicaNumber();
0127
0128   // -- and get the related importance we defined in the importance map:
0129   G4int preImportance = (*fImportanceMap)[preReplicaNumber];
0130   G4int postImportance = (*fImportanceMap)[postReplicaNumber];
0131
0132   // -- Get track weight:
0133   G4double initialWeight = track->GetWeight();
0134
0135   // -- Initialize the "particle change" (it will communicate the new track state to
0136   // -- the tracking):
0137   fParticleChange.Initialize(*track);
0138
0139   if (postImportance > preImportance) {
0140     // -- We split :
0141     G4int splittingFactor = postImportance / preImportance;
0142
0143     // Define the tracks weight:
0144     G4double weightOfTrack = initialWeight / splittingFactor;
0145
0146     // Ask currect track weight to be changed to the new value:
0147     fParticleChange.ProposeParentWeight(weightOfTrack);
0148     // Now we clone this track (this is the actual splitting):
0149     // we will then have the primary and clone of it, hence the
0150     // splitting by a factor 2:
0151     G4Track* clone = new G4Track(*track);
0152     clone->SetWeight(weightOfTrack);
0153     fParticleChange.AddSecondary(clone);
0154     // -- Below's call added for safety & illustration : inform particle change to not
0155     // -- modify the clone (ie : daughter) weight to make it that of the primary.
0156     // -- Here this call is not mandatory as both tracks have same weights.
0157     fParticleChange.SetSecondaryWeightByProcess(true);
0158   }
0159   else {
0160     // -- We apply Russian roulette:
0161     G4double survivingProbability = G4double(postImportance) / G4double(preImportance);
0162
0163     // Shoot a random number (in ]0,1[ segment):
0164     G4double random = G4UniformRand();
0165     if (random > survivingProbability) {
0166       // We ask for the the track to be killed:
0167       fParticleChange.ProposeTrackStatus(fStopAndKill);
0168     }
0169     else {
0170       // In this case, the track survives. We change its weight
0171       // to conserve weight among killed and survival tracks:
0172       fParticleChange.ProposeParentWeight(initialWeight / survivingProbability);
0173     }
0174   }
0175
0176   return &fParticleChange;
0177 }
0178
0179 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......