gorad/src/GRParallelWorldBiasingProcess.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 //  Gorad (Geant4 Open-source Radiation Analysis and Design)
0027 //
0028 //  Author : Makoto Asai (SLAC National Accelerator Laboratory)
0029 //
0030 //  Development of Gorad is funded by NASA Johnson Space Center (JSC)
0031 //  under the contract NNJ15HK11B.
0032 //
0033 // ********************************************************************
0034 //
0035 // GRParallelWorldBiasingProcess.cc
0036 //   A process that takes care of the geometry importance biasing.
0037 //   This class assumes the existance of a parallel world dedicated
0038 //   to the geometry importance biasing and biasing parameters are
0039 //   associated to the region defined in that parallel world.
0040 //
0041 // History
0042 //   September 8th, 2020 : first implementation
0043 //
0044 // ********************************************************************
0045
0046 #include "G4ios.hh"
0047 #include "GRParallelWorldBiasingProcess.hh"
0048 #include "G4Step.hh"
0049 #include "G4Track.hh"
0050 #include "G4Navigator.hh"
0051 #include "G4PathFinder.hh"
0052 #include "G4VTouchable.hh"
0053 #include "G4VPhysicalVolume.hh"
0054 #include "G4ParticleChange.hh"
0055 #include "G4ParticleChangeForNothing.hh"
0056 #include "Randomize.hh"
0057 #include "G4Region.hh"
0058 #include "GRBiasingRegionInfo.hh"
0059
0060 GRParallelWorldBiasingProcess::
0061 GRParallelWorldBiasingProcess(const G4String& processName,G4ProcessType theType)
0062 :G4ParallelWorldProcess(processName,theType)
0063 {
0064   particleChange = new G4ParticleChange();
0065   emptyParticleChange = new G4ParticleChangeForNothing();
0066 }
0067
0068 GRParallelWorldBiasingProcess::~GRParallelWorldBiasingProcess()
0069 {;}
0070
0071 G4VParticleChange* GRParallelWorldBiasingProcess::PostStepDoIt(
0072      const G4Track& track,
0073      const G4Step& step)
0074 {
0075   fOldGhostTouchable = fGhostPostStepPoint->GetTouchableHandle();
0076   CopyStep(step);
0077
0078   if(fOnBoundary)
0079   {
0080     fNewGhostTouchable = fPathFinder->CreateTouchableHandle(fNavigatorID);
0081   }
0082   else
0083   {
0084     fNewGhostTouchable = fOldGhostTouchable;
0085   }
0086
0087   fGhostPreStepPoint->SetTouchableHandle(fOldGhostTouchable);
0088   fGhostPostStepPoint->SetTouchableHandle(fNewGhostTouchable);
0089
0090   if(fNewGhostTouchable->GetVolume() == nullptr)
0091   {
0092     // world volume boundary
0093     emptyParticleChange->Initialize(track);
0094     return emptyParticleChange;
0095   }
0096   if(fNewGhostTouchable->GetVolume() == fOldGhostTouchable->GetVolume())
0097   {
0098     // stayed in the same volume
0099     emptyParticleChange->Initialize(track);
0100     return emptyParticleChange;
0101   }
0102
0103   G4int preStepCopyNo = fOldGhostTouchable->GetReplicaNumber();
0104   G4int postStepCopyNo = fNewGhostTouchable->GetReplicaNumber();
0105
0106   GRBiasingRegionInfo* biasInfo = static_cast<GRBiasingRegionInfo*>
0107     (fNewGhostTouchable->GetVolume()->GetLogicalVolume()->GetRegion()->GetUserInformation());
0108   auto probability = biasInfo->GetProbability();
0109   if(probability < 1.)
0110   {
0111     // skip biasing to avoid over biasing
0112     G4double ran = G4UniformRand();
0113     if(ran > probability)
0114     {
0115       emptyParticleChange->Initialize(track);
0116       return emptyParticleChange;
0117     }
0118   }
0119
0120   G4double trackWeight = track.GetWeight();
0121   particleChange->Initialize(track);
0122
0123   auto biasFactor = biasInfo->GetBiasingFactor();
0124   if(biasFactor==1)
0125   {
0126     // not biasing
0127     emptyParticleChange->Initialize(track);
0128     return emptyParticleChange;
0129   }
0130
0131   if(preStepCopyNo < postStepCopyNo)
0132   {
0133     // getting into more important volume, i.e. do splitting
0134     trackWeight /= biasFactor;
0135     particleChange->ProposeParentWeight(trackWeight);
0136     for(G4int iClone=1;iClone<biasFactor;iClone++)
0137     {
0138       G4Track* clonedTrack = new G4Track(track);
0139       clonedTrack->SetWeight(trackWeight);
0140       particleChange->AddSecondary(clonedTrack);
0141     }
0142     particleChange->SetSecondaryWeightByProcess(true);
0143   }
0144   else if(preStepCopyNo > postStepCopyNo)
0145   {
0146     // getting into less important volume, i.e. do Russian Rouletting
0147     G4double survivalRate = 1.0/biasFactor;
0148     G4double ran = G4UniformRand();
0149     if(ran > survivalRate)
0150     {
0151       // dead
0152       particleChange->ProposeTrackStatus(fStopAndKill);
0153     }
0154     else
0155     {
0156       // survive
0157       trackWeight *= biasFactor;
0158       particleChange->ProposeParentWeight(trackWeight);
0159     }
0160   }
0161   else
0162   {
0163     // This should not happen!!
0164     G4ExceptionDescription ed;
0165     ed << "pre step point :  "<<fOldGhostTouchable->GetVolume()->GetName()<<" - copy no. "<<fOldGhostTouchable->GetReplicaNumber()<<"\n"
0166        << "post step point : "<<fNewGhostTouchable->GetVolume()->GetName()<<" - copy no. "<<fNewGhostTouchable->GetReplicaNumber();
0167     G4Exception("GRParallelWorldBiasingProcess::PostStepDoIt()","GoradProc0001",FatalException,ed);
0168   }
0169
0170   return particleChange;
0171 }
0172