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