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 | // $Id: G4AdjointBremsstrahlungModel.cc,v 1.5 2009/12/16 17:50:01 gunter Exp $ |
---|
27 | // GEANT4 tag $Name: geant4-09-03 $ |
---|
28 | // |
---|
29 | #include "G4AdjointBremsstrahlungModel.hh" |
---|
30 | #include "G4AdjointCSManager.hh" |
---|
31 | #include "G4Integrator.hh" |
---|
32 | #include "G4TrackStatus.hh" |
---|
33 | #include "G4ParticleChange.hh" |
---|
34 | #include "G4AdjointElectron.hh" |
---|
35 | #include "G4AdjointGamma.hh" |
---|
36 | #include "G4Electron.hh" |
---|
37 | |
---|
38 | #include "G4Timer.hh" |
---|
39 | //#include "G4PenelopeBremsstrahlungModel.hh" |
---|
40 | |
---|
41 | |
---|
42 | //////////////////////////////////////////////////////////////////////////////// |
---|
43 | // |
---|
44 | G4AdjointBremsstrahlungModel::G4AdjointBremsstrahlungModel(): |
---|
45 | G4VEmAdjointModel("AdjointeBremModel"), |
---|
46 | MigdalConstant(classic_electr_radius*electron_Compton_length*electron_Compton_length*4.0*pi) |
---|
47 | { |
---|
48 | SetUseMatrix(false); |
---|
49 | SetUseMatrixPerElement(false); |
---|
50 | |
---|
51 | theDirectStdBremModel = new G4eBremsstrahlungModel(G4Electron::Electron(),"TheDirecteBremModel"); |
---|
52 | theDirectEMModel=theDirectStdBremModel; |
---|
53 | // theDirectPenelopeBremModel =0; |
---|
54 | |
---|
55 | SetApplyCutInRange(true); |
---|
56 | highKinEnergy= 100.*TeV; |
---|
57 | lowKinEnergy = 1.0*keV; |
---|
58 | theTimer =new G4Timer(); |
---|
59 | |
---|
60 | theAdjEquivOfDirectPrimPartDef =G4AdjointElectron::AdjointElectron(); |
---|
61 | theAdjEquivOfDirectSecondPartDef=G4AdjointGamma::AdjointGamma(); |
---|
62 | theDirectPrimaryPartDef=G4Electron::Electron(); |
---|
63 | second_part_of_same_type=false; |
---|
64 | |
---|
65 | /*UsePenelopeModel=false; |
---|
66 | if (UsePenelopeModel) { |
---|
67 | G4PenelopeBremsstrahlungModel* thePenelopeModel = new G4PenelopeBremsstrahlungModel(G4Electron::Electron(),"PenelopeBrem"); |
---|
68 | theEmModelManagerForFwdModels = new G4EmModelManager(); |
---|
69 | isPenelopeModelInitialised = false; |
---|
70 | G4VEmFluctuationModel* f=0; |
---|
71 | G4Region* r=0; |
---|
72 | theDirectEMModel=thePenelopeModel; |
---|
73 | theEmModelManagerForFwdModels->AddEmModel(1, thePenelopeModel, f, r); |
---|
74 | } |
---|
75 | */ |
---|
76 | |
---|
77 | |
---|
78 | |
---|
79 | } |
---|
80 | |
---|
81 | //////////////////////////////////////////////////////////////////////////////// |
---|
82 | // |
---|
83 | void G4AdjointBremsstrahlungModel::SampleSecondaries(const G4Track& aTrack, |
---|
84 | G4bool IsScatProjToProjCase, |
---|
85 | G4ParticleChange* fParticleChange) |
---|
86 | { |
---|
87 | if (!UseMatrix) return RapidSampleSecondaries(aTrack,IsScatProjToProjCase,fParticleChange); |
---|
88 | |
---|
89 | const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); |
---|
90 | DefineCurrentMaterial(aTrack.GetMaterialCutsCouple()); |
---|
91 | |
---|
92 | |
---|
93 | G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); |
---|
94 | G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy(); |
---|
95 | |
---|
96 | if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ |
---|
97 | return; |
---|
98 | } |
---|
99 | |
---|
100 | G4double projectileKinEnergy = SampleAdjSecEnergyFromCSMatrix(adjointPrimKinEnergy, |
---|
101 | IsScatProjToProjCase); |
---|
102 | //Weight correction |
---|
103 | //----------------------- |
---|
104 | CorrectPostStepWeight(fParticleChange, |
---|
105 | aTrack.GetWeight(), |
---|
106 | adjointPrimKinEnergy, |
---|
107 | projectileKinEnergy, |
---|
108 | IsScatProjToProjCase); |
---|
109 | |
---|
110 | |
---|
111 | //Kinematic |
---|
112 | //--------- |
---|
113 | G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); |
---|
114 | G4double projectileTotalEnergy = projectileM0+projectileKinEnergy; |
---|
115 | G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; |
---|
116 | G4double projectileP = std::sqrt(projectileP2); |
---|
117 | |
---|
118 | |
---|
119 | //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel |
---|
120 | //------------------------------------------------ |
---|
121 | G4double u; |
---|
122 | const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ; |
---|
123 | |
---|
124 | if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1; |
---|
125 | else u = - std::log(G4UniformRand()*G4UniformRand())/a2; |
---|
126 | |
---|
127 | G4double theta = u*electron_mass_c2/projectileTotalEnergy; |
---|
128 | |
---|
129 | G4double sint = std::sin(theta); |
---|
130 | G4double cost = std::cos(theta); |
---|
131 | |
---|
132 | G4double phi = twopi * G4UniformRand() ; |
---|
133 | |
---|
134 | G4ThreeVector projectileMomentum; |
---|
135 | projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame |
---|
136 | if (IsScatProjToProjCase) {//the adjoint primary is the scattered e- |
---|
137 | G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.); |
---|
138 | G4ThreeVector dirProd=projectileMomentum-gammaMomentum; |
---|
139 | G4double cost1 = std::cos(dirProd.angle(projectileMomentum)); |
---|
140 | G4double sint1 = std::sqrt(1.-cost1*cost1); |
---|
141 | projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP; |
---|
142 | |
---|
143 | } |
---|
144 | |
---|
145 | projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection()); |
---|
146 | |
---|
147 | |
---|
148 | |
---|
149 | if (!IsScatProjToProjCase ){ //kill the primary and add a secondary |
---|
150 | fParticleChange->ProposeTrackStatus(fStopAndKill); |
---|
151 | fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum)); |
---|
152 | } |
---|
153 | else { |
---|
154 | fParticleChange->ProposeEnergy(projectileKinEnergy); |
---|
155 | fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); |
---|
156 | |
---|
157 | } |
---|
158 | } |
---|
159 | //////////////////////////////////////////////////////////////////////////////// |
---|
160 | // |
---|
161 | void G4AdjointBremsstrahlungModel::RapidSampleSecondaries(const G4Track& aTrack, |
---|
162 | G4bool IsScatProjToProjCase, |
---|
163 | G4ParticleChange* fParticleChange) |
---|
164 | { |
---|
165 | |
---|
166 | const G4DynamicParticle* theAdjointPrimary =aTrack.GetDynamicParticle(); |
---|
167 | DefineCurrentMaterial(aTrack.GetMaterialCutsCouple()); |
---|
168 | |
---|
169 | |
---|
170 | G4double adjointPrimKinEnergy = theAdjointPrimary->GetKineticEnergy(); |
---|
171 | G4double adjointPrimTotalEnergy = theAdjointPrimary->GetTotalEnergy(); |
---|
172 | |
---|
173 | if (adjointPrimKinEnergy>HighEnergyLimit*0.999){ |
---|
174 | return; |
---|
175 | } |
---|
176 | |
---|
177 | G4double projectileKinEnergy =0.; |
---|
178 | G4double gammaEnergy=0.; |
---|
179 | G4double diffCSUsed=0.; |
---|
180 | if (!IsScatProjToProjCase){ |
---|
181 | gammaEnergy=adjointPrimKinEnergy; |
---|
182 | G4double Emax = GetSecondAdjEnergyMaxForProdToProjCase(adjointPrimKinEnergy); |
---|
183 | G4double Emin= GetSecondAdjEnergyMinForProdToProjCase(adjointPrimKinEnergy);; |
---|
184 | if (Emin>=Emax) return; |
---|
185 | projectileKinEnergy=Emin*std::pow(Emax/Emin,G4UniformRand()); |
---|
186 | diffCSUsed=lastCZ/projectileKinEnergy; |
---|
187 | |
---|
188 | } |
---|
189 | else { G4double Emax = GetSecondAdjEnergyMaxForScatProjToProjCase(adjointPrimKinEnergy); |
---|
190 | G4double Emin = GetSecondAdjEnergyMinForScatProjToProjCase(adjointPrimKinEnergy,currentTcutForDirectSecond); |
---|
191 | if (Emin>=Emax) return; |
---|
192 | G4double f1=(Emin-adjointPrimKinEnergy)/Emin; |
---|
193 | G4double f2=(Emax-adjointPrimKinEnergy)/Emax/f1; |
---|
194 | //G4cout<<"f1 and f2 "<<f1<<'\t'<<f2<<G4endl; |
---|
195 | projectileKinEnergy=adjointPrimKinEnergy/(1.-f1*std::pow(f2,G4UniformRand())); |
---|
196 | gammaEnergy=projectileKinEnergy-adjointPrimKinEnergy; |
---|
197 | diffCSUsed=lastCZ*adjointPrimKinEnergy/projectileKinEnergy/gammaEnergy; |
---|
198 | |
---|
199 | } |
---|
200 | |
---|
201 | |
---|
202 | |
---|
203 | |
---|
204 | //Weight correction |
---|
205 | //----------------------- |
---|
206 | //First w_corr is set to the ratio between adjoint total CS and fwd total CS |
---|
207 | G4double w_corr=G4AdjointCSManager::GetAdjointCSManager()->GetPostStepWeightCorrection(); |
---|
208 | |
---|
209 | //Then another correction is needed due to the fact that a biaised differential CS has been used rather than the one consistent with the direct model |
---|
210 | //Here we consider the true diffCS as the one obtained by the numericla differentiation over Tcut of the direct CS, corrected by the Migdal term. |
---|
211 | //Basically any other differential CS diffCS could be used here (example Penelope). |
---|
212 | |
---|
213 | G4double diffCS = DiffCrossSectionPerVolumePrimToSecond(currentMaterial, projectileKinEnergy, gammaEnergy); |
---|
214 | w_corr*=diffCS/diffCSUsed; |
---|
215 | |
---|
216 | G4double new_weight = aTrack.GetWeight()*w_corr; |
---|
217 | fParticleChange->SetParentWeightByProcess(false); |
---|
218 | fParticleChange->SetSecondaryWeightByProcess(false); |
---|
219 | fParticleChange->ProposeParentWeight(new_weight); |
---|
220 | |
---|
221 | //Kinematic |
---|
222 | //--------- |
---|
223 | G4double projectileM0 = theAdjEquivOfDirectPrimPartDef->GetPDGMass(); |
---|
224 | G4double projectileTotalEnergy = projectileM0+projectileKinEnergy; |
---|
225 | G4double projectileP2 = projectileTotalEnergy*projectileTotalEnergy - projectileM0*projectileM0; |
---|
226 | G4double projectileP = std::sqrt(projectileP2); |
---|
227 | |
---|
228 | |
---|
229 | //Angle of the gamma direction with the projectile taken from G4eBremsstrahlungModel |
---|
230 | //------------------------------------------------ |
---|
231 | G4double u; |
---|
232 | const G4double a1 = 0.625 , a2 = 3.*a1 , d = 27. ; |
---|
233 | |
---|
234 | if (9./(9.+d) > G4UniformRand()) u = - std::log(G4UniformRand()*G4UniformRand())/a1; |
---|
235 | else u = - std::log(G4UniformRand()*G4UniformRand())/a2; |
---|
236 | |
---|
237 | G4double theta = u*electron_mass_c2/projectileTotalEnergy; |
---|
238 | |
---|
239 | G4double sint = std::sin(theta); |
---|
240 | G4double cost = std::cos(theta); |
---|
241 | |
---|
242 | G4double phi = twopi * G4UniformRand() ; |
---|
243 | |
---|
244 | G4ThreeVector projectileMomentum; |
---|
245 | projectileMomentum=G4ThreeVector(std::cos(phi)*sint,std::sin(phi)*sint,cost)*projectileP; //gamma frame |
---|
246 | if (IsScatProjToProjCase) {//the adjoint primary is the scattered e- |
---|
247 | G4ThreeVector gammaMomentum = (projectileTotalEnergy-adjointPrimTotalEnergy)*G4ThreeVector(0.,0.,1.); |
---|
248 | G4ThreeVector dirProd=projectileMomentum-gammaMomentum; |
---|
249 | G4double cost1 = std::cos(dirProd.angle(projectileMomentum)); |
---|
250 | G4double sint1 = std::sqrt(1.-cost1*cost1); |
---|
251 | projectileMomentum=G4ThreeVector(std::cos(phi)*sint1,std::sin(phi)*sint1,cost1)*projectileP; |
---|
252 | |
---|
253 | } |
---|
254 | |
---|
255 | projectileMomentum.rotateUz(theAdjointPrimary->GetMomentumDirection()); |
---|
256 | |
---|
257 | |
---|
258 | |
---|
259 | if (!IsScatProjToProjCase ){ //kill the primary and add a secondary |
---|
260 | fParticleChange->ProposeTrackStatus(fStopAndKill); |
---|
261 | fParticleChange->AddSecondary(new G4DynamicParticle(theAdjEquivOfDirectPrimPartDef,projectileMomentum)); |
---|
262 | } |
---|
263 | else { |
---|
264 | fParticleChange->ProposeEnergy(projectileKinEnergy); |
---|
265 | fParticleChange->ProposeMomentumDirection(projectileMomentum.unit()); |
---|
266 | |
---|
267 | } |
---|
268 | } |
---|
269 | //////////////////////////////////////////////////////////////////////////////// |
---|
270 | // |
---|
271 | G4AdjointBremsstrahlungModel::~G4AdjointBremsstrahlungModel() |
---|
272 | {;} |
---|
273 | |
---|
274 | //////////////////////////////////////////////////////////////////////////////// |
---|
275 | // |
---|
276 | G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecond(const G4Material* aMaterial, |
---|
277 | G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
278 | G4double kinEnergyProd // kinetic energy of the secondary particle |
---|
279 | ) |
---|
280 | {/*if (UsePenelopeModel && !isPenelopeModelInitialised) { |
---|
281 | theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0); |
---|
282 | isPenelopeModelInitialised =true; |
---|
283 | } |
---|
284 | */ |
---|
285 | return DiffCrossSectionPerVolumePrimToSecondApproximated2(aMaterial, |
---|
286 | kinEnergyProj, |
---|
287 | kinEnergyProd); |
---|
288 | /*return G4VEmAdjointModel::DiffCrossSectionPerVolumePrimToSecond(aMaterial, |
---|
289 | kinEnergyProj, |
---|
290 | kinEnergyProd);*/ |
---|
291 | } |
---|
292 | |
---|
293 | //////////////////////////////////////////////////////////////////////////////// |
---|
294 | // |
---|
295 | G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated1( |
---|
296 | const G4Material* aMaterial, |
---|
297 | G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
298 | G4double kinEnergyProd // kinetic energy of the secondary particle |
---|
299 | ) |
---|
300 | { |
---|
301 | G4double dCrossEprod=0.; |
---|
302 | G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(kinEnergyProd); |
---|
303 | G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(kinEnergyProd); |
---|
304 | |
---|
305 | |
---|
306 | //In this approximation we consider that the secondary gammas are sampled with 1/Egamma energy distribution |
---|
307 | //This is what is applied in the discrete standard model before the rejection test that make a cooerction |
---|
308 | //The application of the same rejection function is not possble here. |
---|
309 | //The differentiation of the CS over Ecut does not produce neither a good differential CS. That is due to the |
---|
310 | // fact that in the discrete model the differential CS and the integrated CS are both fitted but separatly and |
---|
311 | // therefore do not allow a correct numerical differentiation of the integrated CS to get the differential one. |
---|
312 | // In the future we plan to use the brem secondary spectra from the G4Penelope implementation |
---|
313 | |
---|
314 | if (kinEnergyProj>Emin_proj && kinEnergyProj<=Emax_proj){ |
---|
315 | G4double sigma=theDirectEMModel->CrossSectionPerVolume(aMaterial,theDirectPrimaryPartDef,kinEnergyProj,1.*keV); |
---|
316 | dCrossEprod=sigma/kinEnergyProd/std::log(kinEnergyProj/keV); |
---|
317 | } |
---|
318 | return dCrossEprod; |
---|
319 | |
---|
320 | } |
---|
321 | |
---|
322 | //////////////////////////////////////////////////////////////////////////////// |
---|
323 | // |
---|
324 | G4double G4AdjointBremsstrahlungModel::DiffCrossSectionPerVolumePrimToSecondApproximated2( |
---|
325 | const G4Material* material, |
---|
326 | G4double kinEnergyProj, // kinetic energy of the primary particle before the interaction |
---|
327 | G4double kinEnergyProd // kinetic energy of the secondary particle |
---|
328 | ) |
---|
329 | { |
---|
330 | //In this approximation we derive the direct cross section over Tcut=gamma energy, en after apply the Migdla correction factor |
---|
331 | //used in the direct model |
---|
332 | |
---|
333 | G4double dCrossEprod=0.; |
---|
334 | |
---|
335 | const G4ElementVector* theElementVector = material->GetElementVector(); |
---|
336 | const double* theAtomNumDensityVector = material->GetAtomicNumDensityVector(); |
---|
337 | G4double dum=0.; |
---|
338 | G4double E1=kinEnergyProd,E2=kinEnergyProd*1.001; |
---|
339 | G4double dE=E2-E1; |
---|
340 | for (size_t i=0; i<material->GetNumberOfElements(); i++) { |
---|
341 | G4double C1=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum ,E1); |
---|
342 | G4double C2=theDirectEMModel->ComputeCrossSectionPerAtom(theDirectPrimaryPartDef,kinEnergyProj,(*theElementVector)[i]->GetZ(),dum,E2); |
---|
343 | dCrossEprod += theAtomNumDensityVector[i] * (C1-C2)/dE; |
---|
344 | |
---|
345 | } |
---|
346 | |
---|
347 | //Now the Migdal correction |
---|
348 | |
---|
349 | G4double totalEnergy = kinEnergyProj+electron_mass_c2 ; |
---|
350 | G4double kp2 = MigdalConstant*totalEnergy*totalEnergy |
---|
351 | *(material->GetElectronDensity()); |
---|
352 | |
---|
353 | |
---|
354 | G4double MigdalFactor = 1./(1.+kp2/(kinEnergyProd*kinEnergyProd)); // its seems that the factor used in the CS compuation i the direct |
---|
355 | //model is different than the one used in the secondary sampling by a |
---|
356 | //factor (1.+kp2) To be checked! |
---|
357 | |
---|
358 | dCrossEprod*=MigdalFactor; |
---|
359 | return dCrossEprod; |
---|
360 | |
---|
361 | } |
---|
362 | //////////////////////////////////////////////////////////////////////////////// |
---|
363 | // |
---|
364 | G4double G4AdjointBremsstrahlungModel::AdjointCrossSection(const G4MaterialCutsCouple* aCouple, |
---|
365 | G4double primEnergy, |
---|
366 | G4bool IsScatProjToProjCase) |
---|
367 | {/* if (UsePenelopeModel && !isPenelopeModelInitialised) { |
---|
368 | theEmModelManagerForFwdModels->Initialise(G4Electron::Electron(),G4Gamma::Gamma(),1.,0); |
---|
369 | isPenelopeModelInitialised =true; |
---|
370 | } |
---|
371 | */ |
---|
372 | if (UseMatrix) return G4VEmAdjointModel::AdjointCrossSection(aCouple,primEnergy,IsScatProjToProjCase); |
---|
373 | DefineCurrentMaterial(aCouple); |
---|
374 | G4double Cross=0.; |
---|
375 | lastCZ=theDirectEMModel->CrossSectionPerVolume(aCouple->GetMaterial(),theDirectPrimaryPartDef,100.*MeV,100.*MeV/std::exp(1.));//this give the constant above |
---|
376 | |
---|
377 | if (!IsScatProjToProjCase ){ |
---|
378 | G4double Emax_proj = GetSecondAdjEnergyMaxForProdToProjCase(primEnergy); |
---|
379 | G4double Emin_proj = GetSecondAdjEnergyMinForProdToProjCase(primEnergy); |
---|
380 | if (Emax_proj>Emin_proj && primEnergy > currentTcutForDirectSecond) Cross= lastCZ*std::log(Emax_proj/Emin_proj); |
---|
381 | } |
---|
382 | else { |
---|
383 | G4double Emax_proj = GetSecondAdjEnergyMaxForScatProjToProjCase(primEnergy); |
---|
384 | G4double Emin_proj = GetSecondAdjEnergyMinForScatProjToProjCase(primEnergy,currentTcutForDirectSecond); |
---|
385 | if (Emax_proj>Emin_proj) Cross= lastCZ*std::log((Emax_proj-primEnergy)*Emin_proj/Emax_proj/(Emin_proj-primEnergy)); |
---|
386 | |
---|
387 | } |
---|
388 | return Cross; |
---|
389 | } |
---|
390 | |
---|
391 | |
---|
392 | |
---|
393 | |
---|
394 | |
---|