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 27 #include "G4AdjointForcedInteractionForGamma.hh" 28 29 #include "G4AdjointCSManager.hh" 30 #include "G4AdjointGamma.hh" 31 #include "G4MaterialCutsCouple.hh" 32 #include "G4ParticleChange.hh" 33 #include "G4SystemOfUnits.hh" 34 #include "G4VEmAdjointModel.hh" 35 36 G4AdjointForcedInteractionForGamma::G4AdjointForcedInteractionForGamma( 37 const G4String& process_name) 38 : G4VContinuousDiscreteProcess(process_name) 39 , fAdjointComptonModel(nullptr) 40 , fAdjointBremModel(nullptr) 41 { 42 fCSManager = G4AdjointCSManager::GetAdjointCSManager(); 43 fParticleChange = new G4ParticleChange(); 44 } 45 46 ////////////////////////////////////////////////////////////////////////////// 47 G4AdjointForcedInteractionForGamma::~G4AdjointForcedInteractionForGamma() 48 { 49 if(fParticleChange) 50 delete fParticleChange; 51 } 52 53 ////////////////////////////////////////////////////////////////////////////// 54 void G4AdjointForcedInteractionForGamma::ProcessDescription( 55 std::ostream& out) const 56 { 57 out << "Forced interaction for gamma.\n"; 58 } 59 60 ////////////////////////////////////////////////////////////////////////////// 61 void G4AdjointForcedInteractionForGamma::BuildPhysicsTable( 62 const G4ParticleDefinition&) 63 { 64 fCSManager->BuildCrossSectionMatrices(); // it will be done just once 65 fCSManager->BuildTotalSigmaTables(); 66 } 67 68 // Note on weight correction for forced interaction. 69 // For the forced interaction applied here we use a truncated exponential law 70 // for the probability of survival over a fixed total length. This is done by 71 // using a linear transformation of the non-biased probability survival. In 72 // math this is written P'(x)=C1P(x)+C2 , with P(x)=exp(-sum(sigma_ixi)) . x and 73 // L can cross different volumes with different cross section sigma. For forced 74 // interaction, we get the limit conditions: 75 // P'(L)=0 and P'(0)=1 (L can be used over different volumes) 76 // From simple solving of linear equations we 77 // get C1=1/(1-P(L)) and C2=-P(L)/(1-P(L)) 78 // P'(x)=(P(x)-P(L))/(1-P(L)) 79 // For the probability over a step x1 to x2, P'(x1->x2)=P'(x2)/P'(x1). 80 // The effective cross 81 // section is defined -d(P'(x))/dx/P'(x). 82 // We get therefore 83 // sigma_eff = C1sigmaP(x)/(C1P(x)+C2) = sigmaP(x)/(P(x)+C2/C1) 84 // = sigmaP(x)/(P(x)-P(L)) = sigma/(1-P(L)/P(x)) 85 ////////////////////////////////////////////////////////////////////////////// 86 G4VParticleChange* G4AdjointForcedInteractionForGamma::PostStepDoIt( 87 const G4Track& track, const G4Step&) 88 { 89 fParticleChange->Initialize(track); 90 // For the free flight gamma no interaction occurs but a gamma with same 91 // properties is produced for further forced interaction. It is done at the 92 // very beginning of the track so that the weight can be the same 93 if(fCopyGammaForForced) 94 { 95 G4ThreeVector theGammaMomentum = track.GetMomentum(); 96 fParticleChange->AddSecondary( 97 new G4DynamicParticle(G4AdjointGamma::AdjointGamma(), theGammaMomentum)); 98 fParticleChange->SetParentWeightByProcess(false); 99 fParticleChange->SetSecondaryWeightByProcess(false); 100 } 101 else 102 { // Occurrence of forced interaction 103 // Selection of the model to be called 104 G4VEmAdjointModel* theSelectedModel = nullptr; 105 G4bool is_scat_proj_to_proj_case = false; 106 G4double factor=1.; 107 if(!fAdjointComptonModel && !fAdjointBremModel) 108 return fParticleChange; 109 if(!fAdjointComptonModel) 110 { 111 theSelectedModel = fAdjointBremModel; 112 is_scat_proj_to_proj_case = false; 113 // This is needed because the results of it will be used in the post step 114 // do it weight correction inside the model 115 fAdjointBremModel->AdjointCrossSection(track.GetMaterialCutsCouple(), 116 track.GetKineticEnergy(), false); 117 } 118 else if(!fAdjointBremModel) 119 { 120 theSelectedModel = fAdjointComptonModel; 121 is_scat_proj_to_proj_case = true; 122 } 123 else 124 { // Choose the model according to a 50-50 % probability 125 G4double bremAdjCS = fAdjointBremModel->AdjointCrossSection( 126 track.GetMaterialCutsCouple(), track.GetKineticEnergy(), false); 127 if(G4UniformRand() < 0.5) 128 { 129 theSelectedModel = fAdjointBremModel; 130 is_scat_proj_to_proj_case = false; 131 factor=bremAdjCS/fLastAdjCS/0.5; 132 } 133 else 134 { 135 theSelectedModel = fAdjointComptonModel; 136 is_scat_proj_to_proj_case = true; 137 factor=(fLastAdjCS-bremAdjCS)/fLastAdjCS/0.5; 138 } 139 } 140 141 // Compute the weight correction factor 142 G4double invEffectiveAdjointCS = 143 (1. - std::exp(fNbAdjIntLength - fTotNbAdjIntLength)) / fLastAdjCS/fCSBias; 144 145 // Call the selected model without correction of the weight in the model 146 theSelectedModel->SetCorrectWeightForPostStepInModel(false); 147 theSelectedModel 148 ->SetAdditionalWeightCorrectionFactorForPostStepOutsideModel( 149 factor*fLastAdjCS * invEffectiveAdjointCS); 150 theSelectedModel->SampleSecondaries(track, is_scat_proj_to_proj_case, 151 fParticleChange); 152 theSelectedModel->SetCorrectWeightForPostStepInModel(true); 153 154 fContinueGammaAsNewFreeFlight = true; 155 } 156 return fParticleChange; 157 } 158 159 ////////////////////////////////////////////////////////////////////////////// 160 G4VParticleChange* G4AdjointForcedInteractionForGamma::AlongStepDoIt( 161 const G4Track& track, const G4Step&) 162 { 163 fParticleChange->Initialize(track); 164 // Compute nb of interactions length over step length 165 G4ThreeVector position = track.GetPosition(); 166 G4double stepLength = track.GetStep()->GetStepLength(); 167 G4double ekin = track.GetKineticEnergy(); 168 fLastAdjCS = fCSManager->GetTotalAdjointCS(track.GetDefinition(), ekin, 169 track.GetMaterialCutsCouple()); 170 G4double nb_fwd_interaction_length_over_step = 171 stepLength * fCSManager->GetTotalForwardCS(G4AdjointGamma::AdjointGamma(), 172 ekin, 173 track.GetMaterialCutsCouple()); 174 175 G4double nb_adj_interaction_length_over_step = stepLength * fLastAdjCS; 176 G4double fwd_survival_probability = 177 std::exp(-nb_fwd_interaction_length_over_step); 178 G4double mc_induced_survival_probability = 1.; 179 180 if(fFreeFlightGamma) 181 { // for free_flight survival probability stays 1 182 // Accumulate the number of interaction lengths during free flight of gamma 183 fTotNbAdjIntLength += nb_adj_interaction_length_over_step; 184 fAccTrackLength += stepLength; 185 } 186 else 187 { 188 G4double previous_acc_nb_adj_interaction_length = fNbAdjIntLength; 189 fNbAdjIntLength += fCSBias*nb_adj_interaction_length_over_step; 190 theNumberOfInteractionLengthLeft -= fCSBias*nb_adj_interaction_length_over_step; 191 192 // protection against rare race condition 193 if(std::abs(fTotNbAdjIntLength - previous_acc_nb_adj_interaction_length) <= 194 1.e-15) 195 { 196 mc_induced_survival_probability = 1.e50; 197 } 198 else 199 { 200 mc_induced_survival_probability = 201 std::exp(-fNbAdjIntLength) - std::exp(-fTotNbAdjIntLength); 202 mc_induced_survival_probability /= 203 (std::exp(-previous_acc_nb_adj_interaction_length) - 204 std::exp(-fTotNbAdjIntLength)); 205 } 206 } 207 G4double weight_correction = 208 fwd_survival_probability / mc_induced_survival_probability; 209 210 // Caution!!! 211 // It is important to select the weight of the post_step_point as the 212 // current weight and not the weight of the track, as the weight of the track 213 // is changed after having applied all the along_step_do_it. 214 G4double new_weight = 215 weight_correction * track.GetStep()->GetPostStepPoint()->GetWeight(); 216 217 fParticleChange->SetParentWeightByProcess(false); 218 fParticleChange->SetSecondaryWeightByProcess(false); 219 fParticleChange->ProposeParentWeight(new_weight); 220 221 return fParticleChange; 222 } 223 224 ////////////////////////////////////////////////////////////////////////////// 225 G4double 226 G4AdjointForcedInteractionForGamma::PostStepGetPhysicalInteractionLength( 227 const G4Track& track, G4double, G4ForceCondition* condition) 228 { 229 G4int step_id = track.GetCurrentStepNumber(); 230 *condition = NotForced; 231 fCopyGammaForForced = false; 232 G4int track_id = track.GetTrackID(); 233 fFreeFlightGamma = 234 (track_id != fLastFreeFlightTrackId + 1 || fContinueGammaAsNewFreeFlight); 235 if(fFreeFlightGamma) 236 { 237 if(step_id == 1 || fContinueGammaAsNewFreeFlight) 238 { 239 *condition = Forced; 240 // A gamma with same conditions will be generate at next post_step do it 241 // for the forced interaction 242 fCopyGammaForForced = true; 243 fLastFreeFlightTrackId = track_id; 244 fAccTrackLength = 0.; 245 fTotNbAdjIntLength = 0.; 246 fContinueGammaAsNewFreeFlight = false; 247 return 1.e-90; 248 } 249 else 250 { 251 return DBL_MAX; 252 } 253 } 254 else 255 { // compute the interaction length for forced interaction 256 if(step_id == 1) 257 { 258 fCSBias=0.000001/fTotNbAdjIntLength; 259 fTotNbAdjIntLength*=fCSBias; 260 G4double min_val = std::exp(-fTotNbAdjIntLength); 261 theNumberOfInteractionLengthLeft = 262 -std::log(min_val + G4UniformRand() * (1. - min_val)); 263 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft; 264 fNbAdjIntLength = 0.; 265 } 266 G4VPhysicalVolume* thePostPhysVolume = 267 track.GetStep()->GetPreStepPoint()->GetPhysicalVolume(); 268 G4double ekin = track.GetKineticEnergy(); 269 G4double postCS = 0.; 270 if(thePostPhysVolume) 271 { 272 postCS = fCSManager->GetTotalAdjointCS( 273 G4AdjointGamma::AdjointGamma(), ekin, 274 thePostPhysVolume->GetLogicalVolume()->GetMaterialCutsCouple()); 275 } 276 if(postCS > 0.) 277 return theNumberOfInteractionLengthLeft / postCS /fCSBias; 278 else 279 return DBL_MAX; 280 } 281 } 282 283 //////////////////////////////////////////////////////////////////////////////// 284 G4double G4AdjointForcedInteractionForGamma::GetContinuousStepLimit( 285 const G4Track&, G4double, G4double, G4double&) 286 { 287 return DBL_MAX; 288 } 289 290 //////////////////////////////////////////////////////////////////////////////// 291 // Not used in this process but should be implemented as virtual method 292 G4double G4AdjointForcedInteractionForGamma::GetMeanFreePath(const G4Track&, 293 G4double, 294 G4ForceCondition*) 295 { 296 return 0.; 297 } 298