| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 26
27 #include "G4AdjointForcedInteractionForGamma.h 27 #include "G4AdjointForcedInteractionForGamma.hh"
28 28
29 #include "G4AdjointCSManager.hh" 29 #include "G4AdjointCSManager.hh"
30 #include "G4AdjointGamma.hh" 30 #include "G4AdjointGamma.hh"
31 #include "G4MaterialCutsCouple.hh" 31 #include "G4MaterialCutsCouple.hh"
32 #include "G4ParticleChange.hh" 32 #include "G4ParticleChange.hh"
33 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh"
34 #include "G4VEmAdjointModel.hh" 34 #include "G4VEmAdjointModel.hh"
35 35
36 G4AdjointForcedInteractionForGamma::G4AdjointF 36 G4AdjointForcedInteractionForGamma::G4AdjointForcedInteractionForGamma(
37 const G4String& process_name) << 37 G4String process_name)
38 : G4VContinuousDiscreteProcess(process_name) 38 : G4VContinuousDiscreteProcess(process_name)
39 , fAdjointComptonModel(nullptr) 39 , fAdjointComptonModel(nullptr)
40 , fAdjointBremModel(nullptr) 40 , fAdjointBremModel(nullptr)
41 { 41 {
42 fCSManager = G4AdjointCSManager::GetAdj 42 fCSManager = G4AdjointCSManager::GetAdjointCSManager();
43 fParticleChange = new G4ParticleChange(); 43 fParticleChange = new G4ParticleChange();
44 } 44 }
45 45
46 ////////////////////////////////////////////// 46 //////////////////////////////////////////////////////////////////////////////
47 G4AdjointForcedInteractionForGamma::~G4Adjoint 47 G4AdjointForcedInteractionForGamma::~G4AdjointForcedInteractionForGamma()
48 { 48 {
49 if(fParticleChange) 49 if(fParticleChange)
50 delete fParticleChange; 50 delete fParticleChange;
51 } 51 }
52 52
53 ////////////////////////////////////////////// 53 //////////////////////////////////////////////////////////////////////////////
54 void G4AdjointForcedInteractionForGamma::Proce 54 void G4AdjointForcedInteractionForGamma::ProcessDescription(
55 std::ostream& out) const 55 std::ostream& out) const
56 { 56 {
57 out << "Forced interaction for gamma.\n"; 57 out << "Forced interaction for gamma.\n";
58 } 58 }
59 59
60 ////////////////////////////////////////////// 60 //////////////////////////////////////////////////////////////////////////////
61 void G4AdjointForcedInteractionForGamma::Build 61 void G4AdjointForcedInteractionForGamma::BuildPhysicsTable(
62 const G4ParticleDefinition&) 62 const G4ParticleDefinition&)
63 { 63 {
64 fCSManager->BuildCrossSectionMatrices(); // 64 fCSManager->BuildCrossSectionMatrices(); // it will be done just once
65 fCSManager->BuildTotalSigmaTables(); 65 fCSManager->BuildTotalSigmaTables();
66 } 66 }
67 67
68 // Note on weight correction for forced intera 68 // Note on weight correction for forced interaction.
69 // For the forced interaction applied here we 69 // For the forced interaction applied here we use a truncated exponential law
70 // for the probability of survival over a fixe 70 // for the probability of survival over a fixed total length. This is done by
71 // using a linear transformation of the non-bi 71 // using a linear transformation of the non-biased probability survival. In
72 // math this is written P'(x)=C1P(x)+C2 , with 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 differen 73 // L can cross different volumes with different cross section sigma. For forced
74 // interaction, we get the limit conditions: 74 // interaction, we get the limit conditions:
75 // P'(L)=0 and P'(0)=1 (L can be used 75 // P'(L)=0 and P'(0)=1 (L can be used over different volumes)
76 // From simple solving of linear equations we 76 // From simple solving of linear equations we
77 // get C1=1/(1-P(L)) and C2=-P(L)/(1-P(L)) 77 // get C1=1/(1-P(L)) and C2=-P(L)/(1-P(L))
78 // P'(x)=(P(x)-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 79 // For the probability over a step x1 to x2, P'(x1->x2)=P'(x2)/P'(x1).
80 // The effective cross 80 // The effective cross
81 // section is defined -d(P'(x))/dx/P'(x). 81 // section is defined -d(P'(x))/dx/P'(x).
82 // We get therefore 82 // We get therefore
83 // sigma_eff = C1sigmaP(x)/(C1P(x)+C2) = sigma 83 // sigma_eff = C1sigmaP(x)/(C1P(x)+C2) = sigmaP(x)/(P(x)+C2/C1)
84 // = sigmaP(x)/(P(x)-P(L)) = sigma/( 84 // = sigmaP(x)/(P(x)-P(L)) = sigma/(1-P(L)/P(x))
85 ////////////////////////////////////////////// 85 //////////////////////////////////////////////////////////////////////////////
86 G4VParticleChange* G4AdjointForcedInteractionF 86 G4VParticleChange* G4AdjointForcedInteractionForGamma::PostStepDoIt(
87 const G4Track& track, const G4Step&) 87 const G4Track& track, const G4Step&)
88 { 88 {
89 fParticleChange->Initialize(track); 89 fParticleChange->Initialize(track);
90 // For the free flight gamma no interaction 90 // For the free flight gamma no interaction occurs but a gamma with same
91 // properties is produced for further forced 91 // properties is produced for further forced interaction. It is done at the
92 // very beginning of the track so that the w 92 // very beginning of the track so that the weight can be the same
93 if(fCopyGammaForForced) 93 if(fCopyGammaForForced)
94 { 94 {
95 G4ThreeVector theGammaMomentum = track.Get 95 G4ThreeVector theGammaMomentum = track.GetMomentum();
96 fParticleChange->AddSecondary( 96 fParticleChange->AddSecondary(
97 new G4DynamicParticle(G4AdjointGamma::Ad 97 new G4DynamicParticle(G4AdjointGamma::AdjointGamma(), theGammaMomentum));
98 fParticleChange->SetParentWeightByProcess( 98 fParticleChange->SetParentWeightByProcess(false);
99 fParticleChange->SetSecondaryWeightByProce 99 fParticleChange->SetSecondaryWeightByProcess(false);
100 } 100 }
101 else 101 else
102 { // Occurrence of forced interaction 102 { // Occurrence of forced interaction
103 // Selection of the model to be called 103 // Selection of the model to be called
104 G4VEmAdjointModel* theSelectedModel = null 104 G4VEmAdjointModel* theSelectedModel = nullptr;
105 G4bool is_scat_proj_to_proj_case = fals 105 G4bool is_scat_proj_to_proj_case = false;
106 G4double factor=1.; 106 G4double factor=1.;
107 if(!fAdjointComptonModel && !fAdjointBremM 107 if(!fAdjointComptonModel && !fAdjointBremModel)
108 return fParticleChange; 108 return fParticleChange;
109 if(!fAdjointComptonModel) 109 if(!fAdjointComptonModel)
110 { 110 {
111 theSelectedModel = fAdjointBrem 111 theSelectedModel = fAdjointBremModel;
112 is_scat_proj_to_proj_case = false; 112 is_scat_proj_to_proj_case = false;
113 // This is needed because the results of 113 // This is needed because the results of it will be used in the post step
114 // do it weight correction inside the mo 114 // do it weight correction inside the model
115 fAdjointBremModel->AdjointCrossSection(t 115 fAdjointBremModel->AdjointCrossSection(track.GetMaterialCutsCouple(),
116 t 116 track.GetKineticEnergy(), false);
117 } 117 }
118 else if(!fAdjointBremModel) 118 else if(!fAdjointBremModel)
119 { 119 {
120 theSelectedModel = fAdjointComp 120 theSelectedModel = fAdjointComptonModel;
121 is_scat_proj_to_proj_case = true; 121 is_scat_proj_to_proj_case = true;
122 } 122 }
123 else 123 else
124 { // Choose the model according to a 50-5 124 { // Choose the model according to a 50-50 % probability
125 G4double bremAdjCS = fAdjointBremModel-> 125 G4double bremAdjCS = fAdjointBremModel->AdjointCrossSection(
126 track.GetMaterialCutsCouple(), track.G 126 track.GetMaterialCutsCouple(), track.GetKineticEnergy(), false);
127 if(G4UniformRand() < 0.5) 127 if(G4UniformRand() < 0.5)
128 { 128 {
129 theSelectedModel = fAdjointBr 129 theSelectedModel = fAdjointBremModel;
130 is_scat_proj_to_proj_case = false; 130 is_scat_proj_to_proj_case = false;
131 factor=bremAdjCS/fLastAdjCS/0.5; 131 factor=bremAdjCS/fLastAdjCS/0.5;
132 } 132 }
133 else 133 else
134 { 134 {
135 theSelectedModel = fAdjointCo 135 theSelectedModel = fAdjointComptonModel;
136 is_scat_proj_to_proj_case = true; 136 is_scat_proj_to_proj_case = true;
137 factor=(fLastAdjCS-bremAdjCS)/fLastAdj 137 factor=(fLastAdjCS-bremAdjCS)/fLastAdjCS/0.5;
138 } 138 }
139 } 139 }
140 140
141 // Compute the weight correction factor 141 // Compute the weight correction factor
142 G4double invEffectiveAdjointCS = 142 G4double invEffectiveAdjointCS =
143 (1. - std::exp(fNbAdjIntLength - fTotNbA 143 (1. - std::exp(fNbAdjIntLength - fTotNbAdjIntLength)) / fLastAdjCS/fCSBias;
144 144
145 // Call the selected model without correc 145 // Call the selected model without correction of the weight in the model
146 theSelectedModel->SetCorrectWeightForPostS 146 theSelectedModel->SetCorrectWeightForPostStepInModel(false);
147 theSelectedModel 147 theSelectedModel
148 ->SetAdditionalWeightCorrectionFactorFor 148 ->SetAdditionalWeightCorrectionFactorForPostStepOutsideModel(
149 factor*fLastAdjCS * invEffectiveAdjoin 149 factor*fLastAdjCS * invEffectiveAdjointCS);
150 theSelectedModel->SampleSecondaries(track, 150 theSelectedModel->SampleSecondaries(track, is_scat_proj_to_proj_case,
151 fParti 151 fParticleChange);
152 theSelectedModel->SetCorrectWeightForPostS 152 theSelectedModel->SetCorrectWeightForPostStepInModel(true);
153 153
154 fContinueGammaAsNewFreeFlight = true; 154 fContinueGammaAsNewFreeFlight = true;
155 } 155 }
156 return fParticleChange; 156 return fParticleChange;
157 } 157 }
158 158
159 ////////////////////////////////////////////// 159 //////////////////////////////////////////////////////////////////////////////
160 G4VParticleChange* G4AdjointForcedInteractionF 160 G4VParticleChange* G4AdjointForcedInteractionForGamma::AlongStepDoIt(
161 const G4Track& track, const G4Step&) 161 const G4Track& track, const G4Step&)
162 { 162 {
163 fParticleChange->Initialize(track); 163 fParticleChange->Initialize(track);
164 // Compute nb of interactions length over st 164 // Compute nb of interactions length over step length
165 G4ThreeVector position = track.GetPosition() 165 G4ThreeVector position = track.GetPosition();
166 G4double stepLength = track.GetStep()->Ge 166 G4double stepLength = track.GetStep()->GetStepLength();
167 G4double ekin = track.GetKineticEne 167 G4double ekin = track.GetKineticEnergy();
168 fLastAdjCS = fCSManager->GetTotalAdjointCS(t 168 fLastAdjCS = fCSManager->GetTotalAdjointCS(track.GetDefinition(), ekin,
169 t 169 track.GetMaterialCutsCouple());
170 G4double nb_fwd_interaction_length_over_step 170 G4double nb_fwd_interaction_length_over_step =
171 stepLength * fCSManager->GetTotalForwardCS 171 stepLength * fCSManager->GetTotalForwardCS(G4AdjointGamma::AdjointGamma(),
172 172 ekin,
173 173 track.GetMaterialCutsCouple());
174 174
175 G4double nb_adj_interaction_length_over_step 175 G4double nb_adj_interaction_length_over_step = stepLength * fLastAdjCS;
176 G4double fwd_survival_probability = 176 G4double fwd_survival_probability =
177 std::exp(-nb_fwd_interaction_length_over_s 177 std::exp(-nb_fwd_interaction_length_over_step);
178 G4double mc_induced_survival_probability = 1 178 G4double mc_induced_survival_probability = 1.;
179 179
180 if(fFreeFlightGamma) 180 if(fFreeFlightGamma)
181 { // for free_flight survival probability s 181 { // for free_flight survival probability stays 1
182 // Accumulate the number of interaction le 182 // Accumulate the number of interaction lengths during free flight of gamma
183 fTotNbAdjIntLength += nb_adj_interaction_l 183 fTotNbAdjIntLength += nb_adj_interaction_length_over_step;
184 fAccTrackLength += stepLength; 184 fAccTrackLength += stepLength;
185 } 185 }
186 else 186 else
187 { 187 {
188 G4double previous_acc_nb_adj_interaction_l 188 G4double previous_acc_nb_adj_interaction_length = fNbAdjIntLength;
189 fNbAdjIntLength += fCSBias*nb_adj_interact 189 fNbAdjIntLength += fCSBias*nb_adj_interaction_length_over_step;
190 theNumberOfInteractionLengthLeft -= fCSBia 190 theNumberOfInteractionLengthLeft -= fCSBias*nb_adj_interaction_length_over_step;
191 191
192 // protection against rare race condition 192 // protection against rare race condition
193 if(std::abs(fTotNbAdjIntLength - previous_ 193 if(std::abs(fTotNbAdjIntLength - previous_acc_nb_adj_interaction_length) <=
194 1.e-15) 194 1.e-15)
195 { 195 {
196 mc_induced_survival_probability = 1.e50; 196 mc_induced_survival_probability = 1.e50;
197 } 197 }
198 else 198 else
199 { 199 {
200 mc_induced_survival_probability = 200 mc_induced_survival_probability =
201 std::exp(-fNbAdjIntLength) - std::exp( 201 std::exp(-fNbAdjIntLength) - std::exp(-fTotNbAdjIntLength);
202 mc_induced_survival_probability /= 202 mc_induced_survival_probability /=
203 (std::exp(-previous_acc_nb_adj_interac 203 (std::exp(-previous_acc_nb_adj_interaction_length) -
204 std::exp(-fTotNbAdjIntLength)); 204 std::exp(-fTotNbAdjIntLength));
205 } 205 }
206 } 206 }
207 G4double weight_correction = 207 G4double weight_correction =
208 fwd_survival_probability / mc_induced_surv 208 fwd_survival_probability / mc_induced_survival_probability;
209 209
210 // Caution!!! 210 // Caution!!!
211 // It is important to select the weight of t 211 // It is important to select the weight of the post_step_point as the
212 // current weight and not the weight of 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 a 213 // is changed after having applied all the along_step_do_it.
214 G4double new_weight = 214 G4double new_weight =
215 weight_correction * track.GetStep()->GetPo 215 weight_correction * track.GetStep()->GetPostStepPoint()->GetWeight();
216 216
217 fParticleChange->SetParentWeightByProcess(fa 217 fParticleChange->SetParentWeightByProcess(false);
218 fParticleChange->SetSecondaryWeightByProcess 218 fParticleChange->SetSecondaryWeightByProcess(false);
219 fParticleChange->ProposeParentWeight(new_wei 219 fParticleChange->ProposeParentWeight(new_weight);
220 220
221 return fParticleChange; 221 return fParticleChange;
222 } 222 }
223 223
224 ////////////////////////////////////////////// 224 //////////////////////////////////////////////////////////////////////////////
225 G4double 225 G4double
226 G4AdjointForcedInteractionForGamma::PostStepGe 226 G4AdjointForcedInteractionForGamma::PostStepGetPhysicalInteractionLength(
227 const G4Track& track, G4double, G4ForceCondi 227 const G4Track& track, G4double, G4ForceCondition* condition)
228 { 228 {
>> 229 static G4int lastFreeFlightTrackId = 1000;
229 G4int step_id = track.G 230 G4int step_id = track.GetCurrentStepNumber();
230 *condition = NotForc 231 *condition = NotForced;
231 fCopyGammaForForced = false; 232 fCopyGammaForForced = false;
232 G4int track_id = track.G 233 G4int track_id = track.GetTrackID();
233 fFreeFlightGamma = 234 fFreeFlightGamma =
234 (track_id != fLastFreeFlightTrackId + 1 || << 235 (track_id != lastFreeFlightTrackId + 1 || fContinueGammaAsNewFreeFlight);
235 if(fFreeFlightGamma) 236 if(fFreeFlightGamma)
236 { 237 {
237 if(step_id == 1 || fContinueGammaAsNewFree 238 if(step_id == 1 || fContinueGammaAsNewFreeFlight)
238 { 239 {
239 *condition = Forced; 240 *condition = Forced;
240 // A gamma with same conditions will be 241 // A gamma with same conditions will be generate at next post_step do it
241 // for the forced interaction 242 // for the forced interaction
242 fCopyGammaForForced = true; 243 fCopyGammaForForced = true;
243 fLastFreeFlightTrackId = track_i << 244 lastFreeFlightTrackId = track_id;
244 fAccTrackLength = 0.; 245 fAccTrackLength = 0.;
245 fTotNbAdjIntLength = 0.; 246 fTotNbAdjIntLength = 0.;
246 fContinueGammaAsNewFreeFlight = false; 247 fContinueGammaAsNewFreeFlight = false;
247 return 1.e-90; 248 return 1.e-90;
248 } 249 }
249 else 250 else
250 { 251 {
251 return DBL_MAX; 252 return DBL_MAX;
252 } 253 }
253 } 254 }
254 else 255 else
255 { // compute the interaction length for for 256 { // compute the interaction length for forced interaction
256 if(step_id == 1) 257 if(step_id == 1)
257 { 258 {
258 fCSBias=0.000001/fTotNbAdjIntLength; 259 fCSBias=0.000001/fTotNbAdjIntLength;
259 fTotNbAdjIntLength*=fCSBias; 260 fTotNbAdjIntLength*=fCSBias;
260 G4double min_val = std::exp(-fTotNbAdjIn 261 G4double min_val = std::exp(-fTotNbAdjIntLength);
261 theNumberOfInteractionLengthLeft = 262 theNumberOfInteractionLengthLeft =
262 -std::log(min_val + G4UniformRand() * 263 -std::log(min_val + G4UniformRand() * (1. - min_val));
263 theInitialNumberOfInteractionLength = th 264 theInitialNumberOfInteractionLength = theNumberOfInteractionLengthLeft;
264 fNbAdjIntLength = 0. 265 fNbAdjIntLength = 0.;
265 } 266 }
266 G4VPhysicalVolume* thePostPhysVolume = 267 G4VPhysicalVolume* thePostPhysVolume =
267 track.GetStep()->GetPreStepPoint()->GetP 268 track.GetStep()->GetPreStepPoint()->GetPhysicalVolume();
268 G4double ekin = track.GetKineticEnergy() 269 G4double ekin = track.GetKineticEnergy();
269 G4double postCS = 0.; 270 G4double postCS = 0.;
270 if(thePostPhysVolume) 271 if(thePostPhysVolume)
271 { 272 {
272 postCS = fCSManager->GetTotalAdjointCS( 273 postCS = fCSManager->GetTotalAdjointCS(
273 G4AdjointGamma::AdjointGamma(), ekin, 274 G4AdjointGamma::AdjointGamma(), ekin,
274 thePostPhysVolume->GetLogicalVolume()- 275 thePostPhysVolume->GetLogicalVolume()->GetMaterialCutsCouple());
275 } 276 }
276 if(postCS > 0.) 277 if(postCS > 0.)
277 return theNumberOfInteractionLengthLeft 278 return theNumberOfInteractionLengthLeft / postCS /fCSBias;
278 else 279 else
279 return DBL_MAX; 280 return DBL_MAX;
280 } 281 }
281 } 282 }
282 283
283 ////////////////////////////////////////////// 284 ////////////////////////////////////////////////////////////////////////////////
284 G4double G4AdjointForcedInteractionForGamma::G 285 G4double G4AdjointForcedInteractionForGamma::GetContinuousStepLimit(
285 const G4Track&, G4double, G4double, G4double 286 const G4Track&, G4double, G4double, G4double&)
286 { 287 {
287 return DBL_MAX; 288 return DBL_MAX;
288 } 289 }
289 290
290 ////////////////////////////////////////////// 291 ////////////////////////////////////////////////////////////////////////////////
291 // Not used in this process but should be impl 292 // Not used in this process but should be implemented as virtual method
292 G4double G4AdjointForcedInteractionForGamma::G 293 G4double G4AdjointForcedInteractionForGamma::GetMeanFreePath(const G4Track&,
293 294 G4double,
294 295 G4ForceCondition*)
295 { 296 {
296 return 0.; 297 return 0.;
297 } 298 }
298 299