| 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 Authors:
28 M. Omer and R. Hajima on 15 November 2019
29 contact:
30 omer.mohamed@jaea.go.jp and hajima.ryoichi@qst.go.jp
31 Publication Information:
32 1- M. Omer, R. Hajima, Validating polarization effects in gamma-rays elastic scattering by Monte
33 Carlo simulation, New J. Phys., vol. 21, 2019, pp. 113006 (1-10),
34 https://doi.org/10.1088/1367-2630/ab4d8a
35 */
36
37 #include "G4JAEAPolarizedElasticScatteringModel.hh"
38 #include "G4SystemOfUnits.hh"
39 #include "G4AutoLock.hh"
40 namespace { G4Mutex G4JAEAPolarizedElasticScatteringModelMutex = G4MUTEX_INITIALIZER; }
41 using namespace std;
42
43 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
45
46 G4PhysicsFreeVector* G4JAEAPolarizedElasticScatteringModel::dataCS[] = {nullptr};
47 G4DataVector* G4JAEAPolarizedElasticScatteringModel::Polarized_ES_Data[] = {nullptr};
48
49 G4JAEAPolarizedElasticScatteringModel::G4JAEAPolarizedElasticScatteringModel()
50 :G4VEmModel("G4JAEAPolarizedElasticScatteringModel"),isInitialised(false)
51 {
52 fParticleChange = 0;
53 lowEnergyLimit = 100 * keV; //low energy limit for JAEAElasticScattering cross section data
54 fLinearPolarizationSensitvity1=1;
55 fLinearPolarizationSensitvity2=1;
56 fCircularPolarizationSensitvity=1;
57
58 verboseLevel= 0;
59 // Verbosity scale for debugging purposes:
60 // 0 = nothing
61 // 1 = calculation of cross sections, file openings...
62 // 2 = entering in methods
63
64 if(verboseLevel > 0)
65 {
66 G4cout << "G4JAEAPolarizedElasticScatteringModel is constructed " << G4endl;
67 }
68
69 }
70
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
72
73 G4JAEAPolarizedElasticScatteringModel::~G4JAEAPolarizedElasticScatteringModel()
74 {
75 if(IsMaster()) {
76 for(G4int i=0; i<=maxZ; ++i) {
77 if(dataCS[i]) {
78 delete dataCS[i];
79 dataCS[i] = nullptr;
80 }
81 if (Polarized_ES_Data[i]){
82 delete Polarized_ES_Data[i];
83 Polarized_ES_Data[i] = nullptr;
84 }
85 }
86 }
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90
91 void G4JAEAPolarizedElasticScatteringModel::Initialise(const G4ParticleDefinition* particle,
92 const G4DataVector& cuts)
93 {
94 if (verboseLevel > 1)
95 {
96 G4cout << "Calling Initialise() of G4JAEAPolarizedElasticScatteringModel." << G4endl
97 << "Energy range: "
98 << LowEnergyLimit() / eV << " eV - "
99 << HighEnergyLimit() / GeV << " GeV"
100 << G4endl;
101 }
102
103 if(IsMaster()) {
104
105 // Initialise element selector
106 InitialiseElementSelectors(particle, cuts);
107
108 // Access to elements
109 const char* path = G4FindDataDir("G4LEDATA");
110 G4ProductionCutsTable* theCoupleTable =
111 G4ProductionCutsTable::GetProductionCutsTable();
112 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
113
114 for(G4int i=0; i<numOfCouples; ++i)
115 {
116 const G4MaterialCutsCouple* couple =
117 theCoupleTable->GetMaterialCutsCouple(i);
118 const G4Material* material = couple->GetMaterial();
119 const G4ElementVector* theElementVector = material->GetElementVector();
120 std::size_t nelm = material->GetNumberOfElements();
121
122 for (std::size_t j=0; j<nelm; ++j)
123 {
124 G4int Z = G4lrint((*theElementVector)[j]->GetZ());
125 if(Z < 1) { Z = 1; }
126 else if(Z > maxZ) { Z = maxZ; }
127 if( (!dataCS[Z]) ) { ReadData(Z, path); }
128 }
129 }
130 }
131
132 if(isInitialised) { return; }
133 fParticleChange = GetParticleChangeForGamma();
134 isInitialised = true;
135
136 }
137
138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139
140 void G4JAEAPolarizedElasticScatteringModel::InitialiseLocal(const G4ParticleDefinition*,
141 G4VEmModel* masterModel)
142 {
143 SetElementSelectors(masterModel->GetElementSelectors());
144 }
145
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147
148 void G4JAEAPolarizedElasticScatteringModel::ReadData(std::size_t Z, const char* path)
149 {
150 if (verboseLevel > 1)
151 {
152 G4cout << "Calling ReadData() of G4JAEAPolarizedElasticScatteringModel"
153 << G4endl;
154 }
155
156 if(dataCS[Z]) { return; }
157
158 const char* datadir = path;
159 if(!datadir)
160 {
161 datadir = G4FindDataDir("G4LEDATA");
162 if(!datadir)
163 {
164 G4Exception("G4JAEAPolarizedElasticScatteringModel::ReadData()","em0006",
165 FatalException,
166 "Environment variable G4LEDATA not defined");
167 return;
168 }
169 }
170
171 std::ostringstream ostCS;
172 ostCS << datadir << "/JAEAESData/amp_Z_" << Z ;
173 std::ifstream ES_Data_Buffer(ostCS.str().c_str(),ios::binary);
174 if( !ES_Data_Buffer.is_open() )
175 {
176 G4ExceptionDescription ed;
177 ed << "G4JAEAPolarizedElasticScattering Model data file <" << ostCS.str().c_str()
178 << "> is not opened!" << G4endl;
179 G4Exception("G4JAEAPolarizedElasticScatteringModel::ReadData()","em0003",FatalException,
180 ed,"G4LEDATA version should be G4EMLOW7.11 or later. Polarized Elastic Scattering Data are not loaded");
181 return;
182 }
183 else
184 {
185 if(verboseLevel > 3) {
186 G4cout << "File " << ostCS.str()
187 << " is opened by G4JAEAPolarizedElasticScatteringModel" << G4endl;
188 }
189 }
190
191
192 if (!Polarized_ES_Data[Z])
193 Polarized_ES_Data[Z] = new G4DataVector();
194
195 G4float buffer_var;
196 while (ES_Data_Buffer.read(reinterpret_cast<char*>(&buffer_var),sizeof(float)))
197 {
198 Polarized_ES_Data[Z]->push_back(buffer_var);
199 }
200
201 dataCS[Z] = new G4PhysicsFreeVector(300,0.01,3.,/*spline=*/true);
202
203 for (G4int i=0;i<300;++i)
204 dataCS[Z]->PutValues(i,10.*i*1e-3,Polarized_ES_Data[Z]->at(i)*1e-22);
205
206 // Activation of spline interpolation
207 dataCS[Z] ->FillSecondDerivatives();
208 }
209
210 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
211
212 G4double G4JAEAPolarizedElasticScatteringModel::ComputeCrossSectionPerAtom(
213 const G4ParticleDefinition*,
214 G4double GammaEnergy,
215 G4double Z, G4double,
216 G4double, G4double)
217 {
218 //Select the energy-grid point closest to the photon energy
219 // G4double *whichenergy = lower_bound(ESdata[0],ESdata[0]+300,GammaEnergy);
220 // int energyindex = max(0,(int)(whichenergy-ESdata[0]-1));
221
222 if (verboseLevel > 1)
223 {
224 G4cout << "G4JAEAPolarizedElasticScatteringModel::ComputeCrossSectionPerAtom()"
225 << G4endl;
226 }
227
228 if(GammaEnergy < lowEnergyLimit) { return 0.0; }
229
230 G4double xs = 0.0;
231
232 G4int intZ = G4lrint(Z);
233
234 if(intZ < 1 || intZ > maxZ) { return xs; }
235
236 G4PhysicsFreeVector* pv = dataCS[intZ];
237
238 // if element was not initialised
239 // do initialisation safely for MT mode
240 if(!pv) {
241 InitialiseForElement(0, intZ);
242 pv = dataCS[intZ];
243 if(!pv) { return xs; }
244 }
245
246 std::size_t n = pv->GetVectorLength() - 1;
247
248 G4double e = GammaEnergy;
249 if(e >= pv->Energy(n)) {
250 xs = (*pv)[n];
251 } else if(e >= pv->Energy(0)) {
252 xs = pv->Value(e);
253 }
254
255 if(verboseLevel > 0)
256 {
257 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
258 << e << G4endl;
259 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
260 G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0]
261 << G4endl;
262 G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n]
263 << G4endl;
264 G4cout << "*********************************************************"
265 << G4endl;
266 }
267
268 return (xs);
269 }
270
271 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
272
273 void G4JAEAPolarizedElasticScatteringModel::SampleSecondaries(
274 std::vector<G4DynamicParticle*>*,
275 const G4MaterialCutsCouple* couple,
276 const G4DynamicParticle* aDynamicGamma,
277 G4double, G4double)
278 {
279 if (verboseLevel > 1) {
280
281 G4cout << "Calling SampleSecondaries() of G4JAEAPolarizedElasticScatteringModel."
282 << G4endl;
283 }
284 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
285
286 // absorption of low-energy gamma
287 if (photonEnergy0 <= lowEnergyLimit)
288 {
289 fParticleChange->ProposeTrackStatus(fStopAndKill);
290 fParticleChange->SetProposedKineticEnergy(0.);
291 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
292 return ;
293 }
294
295 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
296 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
297 G4int Z = G4lrint(elm->GetZ());
298
299 //Getting the corresponding distrbution
300 G4int energyindex=round(100*photonEnergy0)-1;
301 G4double a1=0, a2=0, a3=0,a4=0;
302 for (G4int i=0;i<=180;++i)
303 {
304 a1=Polarized_ES_Data[Z]->at(4*i+300+181*4*(energyindex));
305 a2=Polarized_ES_Data[Z]->at(4*i+1+300+181*4*(energyindex));
306 a3=Polarized_ES_Data[Z]->at(4*i+2+300+181*4*(energyindex));
307 a4=Polarized_ES_Data[Z]->at(4*i+3+300+181*4*(energyindex));
308 distribution[i]=a1*a1+a2*a2+a3*a3+a4*a4;
309 }
310
311 CLHEP::RandGeneral GenThetaDist(distribution,180);
312 //Intial sampling of the scattering angle. To be updated for the circular polarization
313 G4double theta = CLHEP::pi*GenThetaDist.shoot();
314 //G4double theta =45.*CLHEP::pi/180.;
315 //Theta is in degree to call scattering amplitudes
316 G4int theta_in_degree =round(theta*180./CLHEP::pi);
317
318 //theta_in_degree=45;
319
320 G4double am1=0,am2=0,am3=0,am4=0,aparaSquare=0,aperpSquare=0,
321 apara_aper_Asterisk=0,img_apara_aper_Asterisk=0;
322 am1=Polarized_ES_Data[Z]->at(4*theta_in_degree+300+181*4*(energyindex));
323 am2=Polarized_ES_Data[Z]->at(4*theta_in_degree+1+300+181*4*(energyindex));
324 am3=Polarized_ES_Data[Z]->at(4*theta_in_degree+2+300+181*4*(energyindex));
325 am4=Polarized_ES_Data[Z]->at(4*theta_in_degree+3+300+181*4*(energyindex));
326 aparaSquare=am1*am1+am2*am2;
327 aperpSquare=am3*am3+am4*am4;
328 apara_aper_Asterisk=2*a1*a3+2*a2*a4;
329 img_apara_aper_Asterisk=2*a1*a4-2*a2*a3;
330
331 G4ThreeVector Direction_Unpolarized(0.,0.,0.);
332 G4ThreeVector Direction_Linear1(0.,0.,0.);
333 G4ThreeVector Direction_Linear2(0.,0.,0.);
334 G4ThreeVector Direction_Circular(0.,0.,0.);
335 G4ThreeVector Polarization_Unpolarized(0.,0.,0.);
336 G4ThreeVector Polarization_Linear1(0.,0.,0.);
337 G4ThreeVector Polarization_Linear2(0.,0.,0.);
338 G4ThreeVector Polarization_Circular(0.,0.,0.);
339
340 //Stokes parameters for the incoming and outgoing photon
341 G4double Xi1=0, Xi2=0, Xi3=0, Xi1_Prime=0,Xi2_Prime=0,Xi3_Prime=0;
342
343 //Getting the Stokes parameters for the incoming photon
344 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
345 Xi1=gammaPolarization0.x();
346 Xi2=gammaPolarization0.y();
347 Xi3=gammaPolarization0.z();
348
349 //Polarization vector must be unit vector (5% tolerance)
350 if ((gammaPolarization0.mag())>1.05 || (Xi1*Xi1>1.05) || (Xi2*Xi2>1.05) || (Xi3*Xi3>1.05))
351 {
352 G4Exception("G4JAEAPolarizedElasticScatteringModel::SampleSecondaries()","em1006",
353 JustWarning,
354 "WARNING: G4JAEAPolarizedElasticScatteringModel is only compatible with a unit polarization vector.");
355 return;
356 }
357 //Unpolarized gamma rays
358 if (Xi1==0 && Xi2==0 && Xi3==0)
359 {
360 G4double Phi_Unpolarized=0;
361 if (fLinearPolarizationSensitvity1)
362 Phi_Unpolarized=GeneratePolarizedPhi(aparaSquare,aperpSquare,0.);
363 else
364 Phi_Unpolarized=CLHEP::twopi*G4UniformRand();
365 Direction_Unpolarized.setX(sin(theta)*cos(Phi_Unpolarized));
366 Direction_Unpolarized.setY(sin(theta)*sin(Phi_Unpolarized));
367 Direction_Unpolarized.setZ(cos(theta));
368 Direction_Unpolarized.rotateUz(aDynamicGamma->GetMomentumDirection());
369 Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSquare+aperpSquare);
370 Polarization_Unpolarized.setX(Xi1_Prime);
371 Polarization_Unpolarized.setY(0.);
372 Polarization_Unpolarized.setZ(0.);
373 fParticleChange->ProposeMomentumDirection(Direction_Unpolarized);
374 fParticleChange->ProposePolarization(Polarization_Unpolarized);
375 return;
376 }
377
378 //Linear polarization defined by first Stokes parameter
379 G4double InitialAzimuth=aDynamicGamma->GetMomentumDirection().phi();
380 if(InitialAzimuth<0) InitialAzimuth=InitialAzimuth+CLHEP::twopi;
381
382 G4double Phi_Linear1=0.;
383
384 Phi_Linear1 = GeneratePolarizedPhi(aparaSquare+aperpSquare+Xi1*(aparaSquare-aperpSquare),
385 aparaSquare+aperpSquare-Xi1*(aparaSquare-aperpSquare),InitialAzimuth);
386
387 Xi1_Prime=((aparaSquare-aperpSquare)+Xi1*(aparaSquare+aperpSquare)*cos(2*Phi_Linear1))/
388 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
389 Xi2_Prime=(-Xi1*apara_aper_Asterisk*sin(2*Phi_Linear1))/
390 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
391 Xi3_Prime=(-Xi1*img_apara_aper_Asterisk*sin(2*Phi_Linear1))/
392 ((aparaSquare+aperpSquare)+Xi1*(aparaSquare-aperpSquare)*cos(2*Phi_Linear1));
393 //Store momentum direction and po;arization
394 Direction_Linear1.setX(sin(theta)*cos(Phi_Linear1));
395 Direction_Linear1.setY(sin(theta)*sin(Phi_Linear1));
396 Direction_Linear1.setZ(cos(theta));
397 Polarization_Linear1.setX(Xi1_Prime);
398 Polarization_Linear1.setY(Xi2_Prime);
399 Polarization_Linear1.setZ(Xi3_Prime);
400
401 //Set scattered photon polarization sensitivity
402 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1;
403 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2;
404 Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity;
405
406 G4double dsigmaL1=0.0;
407 if(abs(Xi1)>0.0) dsigmaL1=0.25*((aparaSquare+aperpSquare)*
408 (1+Xi1*Xi1_Prime*cos(2*Phi_Linear1))+
409 (aparaSquare-aperpSquare)*(Xi1*cos(2*Phi_Linear1)+Xi1_Prime)
410 -Xi1*Xi2_Prime*apara_aper_Asterisk*sin(2*Phi_Linear1)-
411 Xi1*Xi3_Prime*img_apara_aper_Asterisk*sin(2*Phi_Linear1));
412
413 //Linear polarization defined by second Stokes parameter
414 //G4double IntialAzimuth=aDynamicGamma->GetMomentumDirection().phi();
415 G4double Phi_Linear2=0.;
416
417 InitialAzimuth=InitialAzimuth-CLHEP::pi/4.;
418 if(InitialAzimuth<0) InitialAzimuth=InitialAzimuth+CLHEP::twopi;
419
420 Phi_Linear2 = GeneratePolarizedPhi(aparaSquare+aperpSquare+Xi1*(aparaSquare-aperpSquare)
421 ,aparaSquare+aperpSquare-Xi1*(aparaSquare-aperpSquare),InitialAzimuth);
422
423 Xi1_Prime=((aparaSquare-aperpSquare)+Xi2*(aparaSquare+aperpSquare)*sin(2*Phi_Linear2))/
424 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
425 Xi2_Prime=(Xi2*apara_aper_Asterisk*cos(2*Phi_Linear2))/
426 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
427 Xi3_Prime=(Xi2*img_apara_aper_Asterisk*cos(2*Phi_Linear2))/
428 ((aparaSquare+aperpSquare)+Xi2*(aparaSquare-aperpSquare)*sin(2*Phi_Linear2));
429 //Store momentum direction and polarization
430 Direction_Linear2.setX(sin(theta)*cos(Phi_Linear2));
431 Direction_Linear2.setY(sin(theta)*sin(Phi_Linear2));
432 Direction_Linear2.setZ(cos(theta));
433 Polarization_Linear2.setX(Xi1_Prime);
434 Polarization_Linear2.setY(Xi2_Prime);
435 Polarization_Linear2.setZ(Xi3_Prime);
436
437 //Set scattered photon polarization sensitivity
438 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1;
439 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2;
440 Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity;
441
442 G4double dsigmaL2=0.0;
443 if(abs(Xi2)>0.0)
444 dsigmaL2=0.25*((aparaSquare+aperpSquare)*(1+Xi2*Xi1_Prime*sin(2*Phi_Linear2))+
445 (aparaSquare-aperpSquare)*(Xi2*sin(2*Phi_Linear2)+Xi1_Prime)
446 +Xi2*Xi2_Prime*apara_aper_Asterisk*cos(2*Phi_Linear2)-
447 Xi2*Xi3_Prime*img_apara_aper_Asterisk*cos(2*Phi_Linear2));
448
449 //Circular polarization
450 G4double Phi_Circular = CLHEP::twopi*G4UniformRand();
451 G4double Theta_Circular = 0.;
452
453 Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSquare+aperpSquare);
454 Xi2_Prime=(-Xi3*img_apara_aper_Asterisk)/(aparaSquare+aperpSquare);
455 Xi3_Prime=(Xi3*apara_aper_Asterisk)/(aparaSquare+aperpSquare);
456
457 Polarization_Circular.setX(Xi1_Prime);
458 Polarization_Circular.setY(Xi2_Prime);
459 Polarization_Circular.setZ(Xi3_Prime);
460
461 //Set scattered photon polarization sensitivity
462 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensitvity1;
463 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensitvity2;
464 Xi3_Prime=Xi3_Prime*fCircularPolarizationSensitvity;
465
466 G4double dsigmaC=0.0;
467 if(abs(Xi3)>0.0)
468 dsigmaC=0.25*(aparaSquare+aperpSquare+Xi1_Prime*(aparaSquare-aperpSquare)-
469 Xi3*Xi2_Prime*img_apara_aper_Asterisk
470 +Xi3*Xi3_Prime*apara_aper_Asterisk);
471
472 if (abs(Xi3)==0.0 && abs(Xi1_Prime)==0.0)
473 {
474 Direction_Circular.setX(sin(theta)*cos(Phi_Circular));
475 Direction_Circular.setY(sin(theta)*sin(Phi_Circular));
476 Direction_Circular.setZ(cos(theta));
477 }
478 else
479 {
480 G4double c1=0, c2=0, c3=0,c4=0;
481 for (G4int i=0;i<=180;++i)
482 {
483 c1=Polarized_ES_Data[Z]->at(4*i+300+181*4*(energyindex));
484 c2=Polarized_ES_Data[Z]->at(4*i+1+300+181*4*(energyindex));
485 c3=Polarized_ES_Data[Z]->at(4*i+2+300+181*4*(energyindex));
486 c4=Polarized_ES_Data[Z]->at(4*i+3+300+181*4*(energyindex));
487 cdistribution[i]=0.25*((c1*c1+c2*c2+c3*c3+c4*c4)+
488 Xi1_Prime*(c1*c1+c2*c2-c3*c3-c4*c4)-
489 Xi3*Xi2_Prime*(2*c1*c4-2*c2*c3)
490 +Xi3*Xi3_Prime*(2*c1*c4-2*c2*c3));
491 }
492 CLHEP::RandGeneral GenTheta_Circ_Dist(cdistribution,180);
493 Theta_Circular=CLHEP::pi*GenTheta_Circ_Dist.shoot();
494 Direction_Circular.setX(sin(Theta_Circular)*cos(Phi_Circular));
495 Direction_Circular.setY(sin(Theta_Circular)*sin(Phi_Circular));
496 Direction_Circular.setZ(cos(Theta_Circular));
497 }
498
499 // Sampling scattered photon direction based on asymmetry arising from polarization mixing
500 G4double totalSigma= dsigmaL1+dsigmaL2+dsigmaC;
501 G4double prob1=dsigmaL1/totalSigma;
502 G4double prob2=dsigmaL2/totalSigma;
503 G4double probc=1-(prob1+prob2);
504
505 //Check the Probability of polarization mixing
506 if (abs(probc - dsigmaC/totalSigma)>=0.0001)
507 {
508 G4Exception("G4JAEAPolarizedElasticScatteringModel::SampleSecondaries()","em1007",
509 JustWarning,
510 "WARNING: Polarization mixing might be incorrect.");
511 }
512
513 // Generate outgoing photon direction
514 G4ThreeVector finaldirection(0.0,0.0,0.0);
515 G4ThreeVector outcomingPhotonPolarization(0.0,0.0,0.0);
516
517 //Polarization mixing
518 G4double polmix=G4UniformRand();
519 if (polmix<=prob1)
520 {
521 finaldirection.setX(Direction_Linear1.x());
522 finaldirection.setY(Direction_Linear1.y());
523 finaldirection.setZ(Direction_Linear1.z());
524 outcomingPhotonPolarization.setX(Polarization_Linear1.x());
525 outcomingPhotonPolarization.setY(Polarization_Linear1.y());
526 outcomingPhotonPolarization.setZ(Polarization_Linear1.z());
527 }
528 else if ((polmix>prob1) && (polmix<=prob1+prob2))
529 {
530 finaldirection.setX(Direction_Linear2.x());
531 finaldirection.setY(Direction_Linear2.y());
532 finaldirection.setZ(Direction_Linear2.z());
533 outcomingPhotonPolarization.setX(Polarization_Linear2.x());
534 outcomingPhotonPolarization.setY(Polarization_Linear2.y());
535 outcomingPhotonPolarization.setZ(Polarization_Linear2.z());
536 }
537 else if (polmix>prob1+prob2)
538 {
539 finaldirection.setX(Direction_Circular.x());
540 finaldirection.setY(Direction_Circular.y());
541 finaldirection.setZ(Direction_Circular.z());
542 outcomingPhotonPolarization.setX(Polarization_Circular.x());
543 outcomingPhotonPolarization.setY(Polarization_Circular.y());
544 outcomingPhotonPolarization.setZ(Polarization_Circular.z());
545 }
546
547 //Sampling the Final State
548 finaldirection.rotateUz(aDynamicGamma->GetMomentumDirection());
549 fParticleChange->ProposeMomentumDirection(finaldirection);
550 fParticleChange->SetProposedKineticEnergy(photonEnergy0);
551 fParticleChange->ProposePolarization(outcomingPhotonPolarization);
552
553 }
554
555 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
556
557 G4double G4JAEAPolarizedElasticScatteringModel::GeneratePolarizedPhi(G4double Sigma_para,
558 G4double Sigma_perp,
559 G4double initial_Pol_Plane)
560 {
561 G4double phi;
562 G4double phiProbability;
563 G4double Probability=Sigma_perp/(Sigma_para+Sigma_perp);
564 if (Probability<=G4UniformRand())
565 {
566 do
567 {
568 phi = CLHEP::twopi * G4UniformRand();
569 phiProbability = cos(phi+initial_Pol_Plane)*cos(phi+initial_Pol_Plane);
570 }
571 while (phiProbability < G4UniformRand());
572
573 }
574 else
575 {
576 do
577 {
578 phi = CLHEP::twopi * G4UniformRand();
579 phiProbability = sin(phi+initial_Pol_Plane)*sin(phi+initial_Pol_Plane);
580 }
581 while (phiProbability < G4UniformRand());
582 }
583 return phi;
584
585 }
586 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
587
588 void
589 G4JAEAPolarizedElasticScatteringModel::InitialiseForElement(const G4ParticleDefinition*,
590 G4int Z)
591 {
592 G4AutoLock l(&G4JAEAPolarizedElasticScatteringModelMutex);
593 // G4cout << "G4JAEAPolarizedElasticScatteringModel::InitialiseForElement Z= "
594 // << Z << G4endl;
595 if(!dataCS[Z]) { ReadData(Z); }
596 l.unlock();
597 }
598
599 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
600