Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // Gorad (Geant4 Open-source Radiation Analysis and Design) 27 // 28 // Author : Makoto Asai (SLAC National Accelerator Laboratory) 29 // 30 // Development of Gorad is funded by NASA Johnson Space Center (JSC) 31 // under the contract NNJ15HK11B. 32 // 33 // ******************************************************************** 34 // 35 // GRParallelWorldBiasingProcess.cc 36 // A process that takes care of the geometry importance biasing. 37 // This class assumes the existance of a parallel world dedicated 38 // to the geometry importance biasing and biasing parameters are 39 // associated to the region defined in that parallel world. 40 // 41 // History 42 // September 8th, 2020 : first implementation 43 // 44 // ******************************************************************** 45 46 #include "G4ios.hh" 47 #include "GRParallelWorldBiasingProcess.hh" 48 #include "G4Step.hh" 49 #include "G4Track.hh" 50 #include "G4Navigator.hh" 51 #include "G4PathFinder.hh" 52 #include "G4VTouchable.hh" 53 #include "G4VPhysicalVolume.hh" 54 #include "G4ParticleChange.hh" 55 #include "G4ParticleChangeForNothing.hh" 56 #include "Randomize.hh" 57 #include "G4Region.hh" 58 #include "GRBiasingRegionInfo.hh" 59 60 GRParallelWorldBiasingProcess:: 61 GRParallelWorldBiasingProcess(const G4String& processName,G4ProcessType theType) 62 :G4ParallelWorldProcess(processName,theType) 63 { 64 particleChange = new G4ParticleChange(); 65 emptyParticleChange = new G4ParticleChangeForNothing(); 66 } 67 68 GRParallelWorldBiasingProcess::~GRParallelWorldBiasingProcess() 69 {;} 70 71 G4VParticleChange* GRParallelWorldBiasingProcess::PostStepDoIt( 72 const G4Track& track, 73 const G4Step& step) 74 { 75 fOldGhostTouchable = fGhostPostStepPoint->GetTouchableHandle(); 76 CopyStep(step); 77 78 if(fOnBoundary) 79 { 80 fNewGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID); 81 } 82 else 83 { 84 fNewGhostTouchable = fOldGhostTouchable; 85 } 86 87 fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable); 88 fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable); 89 90 if(fNewGhostTouchable->GetVolume() == nullptr) 91 { 92 // world volume boundary 93 emptyParticleChange->Initialize(track); 94 return emptyParticleChange; 95 } 96 if(fNewGhostTouchable->GetVolume() == fOldGhostTouchable->GetVolume()) 97 { 98 // stayed in the same volume 99 emptyParticleChange->Initialize(track); 100 return emptyParticleChange; 101 } 102 103 G4int preStepCopyNo = fOldGhostTouchable->GetReplicaNumber(); 104 G4int postStepCopyNo = fNewGhostTouchable->GetReplicaNumber(); 105 106 GRBiasingRegionInfo* biasInfo = static_cast<GRBiasingRegionInfo*> 107 (fNewGhostTouchable->GetVolume()->GetLogicalVolume()->GetRegion()->GetUserInformation()); 108 auto probability = biasInfo->GetProbability(); 109 if(probability < 1.) 110 { 111 // skip biasing to avoid over biasing 112 G4double ran = G4UniformRand(); 113 if(ran > probability) 114 { 115 emptyParticleChange->Initialize(track); 116 return emptyParticleChange; 117 } 118 } 119 120 G4double trackWeight = track.GetWeight(); 121 particleChange->Initialize(track); 122 123 auto biasFactor = biasInfo->GetBiasingFactor(); 124 if(biasFactor==1) 125 { 126 // not biasing 127 emptyParticleChange->Initialize(track); 128 return emptyParticleChange; 129 } 130 131 if(preStepCopyNo < postStepCopyNo) 132 { 133 // getting into more important volume, i.e. do splitting 134 trackWeight /= biasFactor; 135 particleChange->ProposeParentWeight(trackWeight); 136 for(G4int iClone=1;iClone<biasFactor;iClone++) 137 { 138 G4Track* clonedTrack = new G4Track(track); 139 clonedTrack->SetWeight(trackWeight); 140 particleChange->AddSecondary(clonedTrack); 141 } 142 particleChange->SetSecondaryWeightByProcess(true); 143 } 144 else if(preStepCopyNo > postStepCopyNo) 145 { 146 // getting into less important volume, i.e. do Russian Rouletting 147 G4double survivalRate = 1.0/biasFactor; 148 G4double ran = G4UniformRand(); 149 if(ran > survivalRate) 150 { 151 // dead 152 particleChange->ProposeTrackStatus(fStopAndKill); 153 } 154 else 155 { 156 // survive 157 trackWeight *= biasFactor; 158 particleChange->ProposeParentWeight(trackWeight); 159 } 160 } 161 else 162 { 163 // This should not happen!! 164 G4ExceptionDescription ed; 165 ed << "pre step point : "<<fOldGhostTouchable->GetVolume()->GetName()<<" - copy no. "<<fOldGhostTouchable->GetReplicaNumber()<<"\n" 166 << "post step point : "<<fNewGhostTouchable->GetVolume()->GetName()<<" - copy no. "<<fNewGhostTouchable->GetReplicaNumber(); 167 G4Exception("GRParallelWorldBiasingProcess::PostStepDoIt()","GoradProc0001",FatalException,ed); 168 } 169 170 return particleChange; 171 } 172 173