| 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 Authors:
28 M. Omer and R. Hajima on 15 November 2019
29 contact:
30 omer.mohamed@jaea.go.jp and hajima.ryoichi@q
31 Publication Information:
32 1- M. Omer, R. Hajima, Validating polarizati
33 Carlo simulation, New J. Phys., vol. 21, 201
34 https://doi.org/10.1088/1367-2630/ab4d8a
35 */
36
37 #include "G4JAEAPolarizedElasticScatteringMode
38 #include "G4SystemOfUnits.hh"
39 #include "G4AutoLock.hh"
40 namespace { G4Mutex G4JAEAPolarizedElasticScat
41 using namespace std;
42
43 //....oooOO0OOooo........oooOO0OOooo........oo
44 //....oooOO0OOooo........oooOO0OOooo........oo
45
46 G4PhysicsFreeVector* G4JAEAPolarizedElasticSca
47 G4DataVector* G4JAEAPolarizedElasticScattering
48
49 G4JAEAPolarizedElasticScatteringModel::G4JAEAP
50 :G4VEmModel("G4JAEAPolarizedElasticScatterin
51 {
52 fParticleChange = 0;
53 lowEnergyLimit = 100 * keV; //low energy
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 o
62 // 2 = entering in methods
63
64 if(verboseLevel > 0)
65 {
66 G4cout << "G4JAEAPolarizedElasticScatter
67 }
68
69 }
70
71 //....oooOO0OOooo........oooOO0OOooo........oo
72
73 G4JAEAPolarizedElasticScatteringModel::~G4JAEA
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........oo
90
91 void G4JAEAPolarizedElasticScatteringModel::In
92 const G4DataVector& cuts)
93 {
94 if (verboseLevel > 1)
95 {
96 G4cout << "Calling Initialise() of G4JAE
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::GetProductionCuts
112 G4int numOfCouples = (G4int)theCoupleTable
113
114 for(G4int i=0; i<numOfCouples; ++i)
115 {
116 const G4MaterialCutsCouple* couple =
117 theCoupleTable->GetMaterialCutsCouple(i);
118 const G4Material* material = couple->GetMate
119 const G4ElementVector* theElementVector = ma
120 std::size_t nelm = material->GetNumberOfElem
121
122 for (std::size_t j=0; j<nelm; ++j)
123 {
124 G4int Z = G4lrint((*theElementVector)[j]
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........oo
139
140 void G4JAEAPolarizedElasticScatteringModel::In
141 G4VEmModel* masterModel)
142 {
143 SetElementSelectors(masterModel->GetElementS
144 }
145
146 //....oooOO0OOooo........oooOO0OOooo........oo
147
148 void G4JAEAPolarizedElasticScatteringModel::Re
149 {
150 if (verboseLevel > 1)
151 {
152 G4cout << "Calling ReadData() of G4JAEAP
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("G4JAEAPolarizedElasticScatter
165 FatalException,
166 "Environment variable G4LEDATA not d
167 return;
168 }
169 }
170
171 std::ostringstream ostCS;
172 ostCS << datadir << "/JAEAESData/amp_Z_" <<
173 std::ifstream ES_Data_Buffer(ostCS.str().c_s
174 if( !ES_Data_Buffer.is_open() )
175 {
176 G4ExceptionDescription ed;
177 ed << "G4JAEAPolarizedElasticScattering
178 << "> is not opened!" << G4endl;
179 G4Exception("G4JAEAPolarizedElasticScatt
180 ed,"G4LEDATA version should be G4EMLOW7.
181 return;
182 }
183 else
184 {
185 if(verboseLevel > 3) {
186 G4cout << "File " << ostCS.str()
187 << " is opened by G4JAEAPolarizedElas
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<
197 {
198 Polarized_ES_Data[Z]->push_back(buffer_v
199 }
200
201 dataCS[Z] = new G4PhysicsFreeVector(300,0.01
202
203 for (G4int i=0;i<300;++i)
204 dataCS[Z]->PutValues(i,10.*i*1e-3,Polarize
205
206 // Activation of spline interpolation
207 dataCS[Z] ->FillSecondDerivatives();
208 }
209
210 //....oooOO0OOooo........oooOO0OOooo........oo
211
212 G4double G4JAEAPolarizedElasticScatteringModel
213 const G4ParticleDefinitio
214 G4double GammaEnergy,
215 G4double Z, G4double,
216 G4double, G4double)
217 {
218 //Select the energy-grid point closest to th
219 // G4double *whichenergy = lower_bound(ES
220 // int energyindex = max(0,(int)(whichene
221
222 if (verboseLevel > 1)
223 {
224 G4cout << "G4JAEAPolarizedElasticScatter
225 << G4endl;
226 }
227
228 if(GammaEnergy < lowEnergyLimit) { return 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
258 << e << G4endl;
259 G4cout << " cs (Geant4 internal unit)
260 G4cout << " -> first E*E*cs value i
261 << G4endl;
262 G4cout << " -> last E*E*cs value i
263 << G4endl;
264 G4cout << "***************************
265 << G4endl;
266 }
267
268 return (xs);
269 }
270
271 //....oooOO0OOooo........oooOO0OOooo........oo
272
273 void G4JAEAPolarizedElasticScatteringModel::Sa
274 std::vector<G4DynamicParti
275 const G4MaterialCutsCouple
276 const G4DynamicParticle* a
277 G4double, G4double)
278 {
279 if (verboseLevel > 1) {
280
281 G4cout << "Calling SampleSecondaries() of
282 << G4endl;
283 }
284 G4double photonEnergy0 = aDynamicGamma->GetK
285
286 // absorption of low-energy gamma
287 if (photonEnergy0 <= lowEnergyLimit)
288 {
289 fParticleChange->ProposeTrackStatus(fSto
290 fParticleChange->SetProposedKineticEnerg
291 fParticleChange->ProposeLocalEnergyDepos
292 return ;
293 }
294
295 const G4ParticleDefinition* particle = aDyn
296 const G4Element* elm = SelectRandomAtom(coup
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*
305 a2=Polarized_ES_Data[Z]->at(4*i+1+300+18
306 a3=Polarized_ES_Data[Z]->at(4*i+2+300+18
307 a4=Polarized_ES_Data[Z]->at(4*i+3+300+18
308 distribution[i]=a1*a1+a2*a2+a3*a3+a4*a4;
309 }
310
311 CLHEP::RandGeneral GenThetaDist(distribution
312 //Intial sampling of the scattering angle. T
313 G4double theta = CLHEP::pi*GenThetaDist.shoo
314 //G4double theta =45.*CLHEP::pi/180.;
315 //Theta is in degree to call scattering ampl
316 G4int theta_in_degree =round(theta*180./CLHE
317
318 //theta_in_degree=45;
319
320 G4double am1=0,am2=0,am3=0,am4=0,aparaSquare
321 apara_aper_Asterisk=0,img_apara_aper_Aster
322 am1=Polarized_ES_Data[Z]->at(4*theta_in_degr
323 am2=Polarized_ES_Data[Z]->at(4*theta_in_degr
324 am3=Polarized_ES_Data[Z]->at(4*theta_in_degr
325 am4=Polarized_ES_Data[Z]->at(4*theta_in_degr
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.
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 out
341 G4double Xi1=0, Xi2=0, Xi3=0, Xi1_Prime=0,Xi
342
343 //Getting the Stokes parameters for the inco
344 G4ThreeVector gammaPolarization0 = aDynamicG
345 Xi1=gammaPolarization0.x();
346 Xi2=gammaPolarization0.y();
347 Xi3=gammaPolarization0.z();
348
349 //Polarization vector must be unit vector (5
350 if ((gammaPolarization0.mag())>1.05 || (Xi1*
351 {
352 G4Exception("G4JAEAPolarizedElasticScatt
353 JustWarning,
354 "WARNING: G4JAEAPolarizedElasticScatteri
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(aparaSq
363 else
364 Phi_Unpolarized=CLHEP::twopi*G4UniformRand()
365 Direction_Unpolarized.setX(sin(theta)*co
366 Direction_Unpolarized.setY(sin(theta)*si
367 Direction_Unpolarized.setZ(cos(theta));
368 Direction_Unpolarized.rotateUz(aDynamicG
369 Xi1_Prime=(aparaSquare-aperpSquare)/(apa
370 Polarization_Unpolarized.setX(Xi1_Prime)
371 Polarization_Unpolarized.setY(0.);
372 Polarization_Unpolarized.setZ(0.);
373 fParticleChange->ProposeMomentumDirectio
374 fParticleChange->ProposePolarization(Pol
375 return;
376 }
377
378 //Linear polarization defined by first Stoke
379 G4double InitialAzimuth=aDynamicGamma->GetMo
380 if(InitialAzimuth<0) InitialAzimuth=InitialA
381
382 G4double Phi_Linear1=0.;
383
384 Phi_Linear1 = GeneratePolarizedPhi(aparaSqua
385 aparaSquare+aperpSquare-Xi1*(apar
386
387 Xi1_Prime=((aparaSquare-aperpSquare)+Xi1*(ap
388 ((aparaSquare+aperpSquare)+Xi1*(aparaSquar
389 Xi2_Prime=(-Xi1*apara_aper_Asterisk*sin(2*Ph
390 ((aparaSquare+aperpSquare)+Xi1*(aparaSquar
391 Xi3_Prime=(-Xi1*img_apara_aper_Asterisk*sin(
392 ((aparaSquare+aperpSquare)+Xi1*(aparaSquar
393 //Store momentum direction and po;arization
394 Direction_Linear1.setX(sin(theta)*cos(Phi_Li
395 Direction_Linear1.setY(sin(theta)*sin(Phi_Li
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 sensitiv
402 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensi
403 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensi
404 Xi3_Prime=Xi3_Prime*fCircularPolarizationSen
405
406 G4double dsigmaL1=0.0;
407 if(abs(Xi1)>0.0) dsigmaL1=0.25*((aparaSquare
408 (1+Xi1*Xi1_Prime*cos(2*Phi_Linear1))
409 (aparaSquare-aperpSquare)*(Xi1*cos(2
410 -Xi1*Xi2_Prime*apara_aper_Asterisk*s
411 Xi1*Xi3_Prime*img_apara_aper_Asteris
412
413 //Linear polarization defined by second Stok
414 //G4double IntialAzimuth=aDynamicGamma->GetM
415 G4double Phi_Linear2=0.;
416
417 InitialAzimuth=InitialAzimuth-CLHEP::pi/4.;
418 if(InitialAzimuth<0) InitialAzimuth=InitialA
419
420 Phi_Linear2 = GeneratePolarizedPhi(aparaSqua
421 ,aparaSquare+aperpSquare-Xi1*(apa
422
423 Xi1_Prime=((aparaSquare-aperpSquare)+Xi2*(ap
424 ((aparaSquare+aperpSquare)+Xi2*(aparaSquar
425 Xi2_Prime=(Xi2*apara_aper_Asterisk*cos(2*Phi
426 ((aparaSquare+aperpSquare)+Xi2*(aparaSquar
427 Xi3_Prime=(Xi2*img_apara_aper_Asterisk*cos(2
428 ((aparaSquare+aperpSquare)+Xi2*(aparaSquar
429 //Store momentum direction and polarization
430 Direction_Linear2.setX(sin(theta)*cos(Phi_Li
431 Direction_Linear2.setY(sin(theta)*sin(Phi_Li
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 sensitiv
438 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensi
439 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensi
440 Xi3_Prime=Xi3_Prime*fCircularPolarizationSen
441
442 G4double dsigmaL2=0.0;
443 if(abs(Xi2)>0.0)
444 dsigmaL2=0.25*((aparaSquare+aperpSquare)*(
445 (aparaSquare-aperpSquare)*(Xi2*sin(2*Ph
446 +Xi2*Xi2_Prime*apara_aper_Asterisk*cos(
447 Xi2*Xi3_Prime*img_apara_aper_Asterisk*c
448
449 //Circular polarization
450 G4double Phi_Circular = CLHEP::twopi*G4Unifo
451 G4double Theta_Circular = 0.;
452
453 Xi1_Prime=(aparaSquare-aperpSquare)/(aparaSq
454 Xi2_Prime=(-Xi3*img_apara_aper_Asterisk)/(ap
455 Xi3_Prime=(Xi3*apara_aper_Asterisk)/(aparaSq
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 sensitiv
462 Xi1_Prime=Xi1_Prime*fLinearPolarizationSensi
463 Xi2_Prime=Xi2_Prime*fLinearPolarizationSensi
464 Xi3_Prime=Xi3_Prime*fCircularPolarizationSen
465
466 G4double dsigmaC=0.0;
467 if(abs(Xi3)>0.0)
468 dsigmaC=0.25*(aparaSquare+aperpSquare+Xi1_
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(P
475 Direction_Circular.setY(sin(theta)*sin(P
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*
484 c2=Polarized_ES_Data[Z]->at(4*i+1+300+181*
485 c3=Polarized_ES_Data[Z]->at(4*i+2+300+181*
486 c4=Polarized_ES_Data[Z]->at(4*i+3+300+181*
487 cdistribution[i]=0.25*((c1*c1+c2*c2+c3*c3+
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(cd
493 Theta_Circular=CLHEP::pi*GenTheta_Circ_D
494 Direction_Circular.setX(sin(Theta_Circul
495 Direction_Circular.setY(sin(Theta_Circul
496 Direction_Circular.setZ(cos(Theta_Circul
497 }
498
499 // Sampling scattered photon direction based
500 G4double totalSigma= dsigmaL1+dsigmaL2+dsigm
501 G4double prob1=dsigmaL1/totalSigma;
502 G4double prob2=dsigmaL2/totalSigma;
503 G4double probc=1-(prob1+prob2);
504
505 //Check the Probability of polarization mixi
506 if (abs(probc - dsigmaC/totalSigma)>=0.0001)
507 {
508 G4Exception("G4JAEAPolarizedElasticScatt
509 JustWarning,
510 "WARNING: Polarization mixing might be i
511 }
512
513 // Generate outgoing photon direction
514 G4ThreeVector finaldirection(0.0,0.0,0.0);
515 G4ThreeVector outcomingPhotonPolarization(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(Polariz
525 outcomingPhotonPolarization.setY(Polariz
526 outcomingPhotonPolarization.setZ(Polariz
527 }
528 else if ((polmix>prob1) && (polmix<=prob1+pr
529 {
530 finaldirection.setX(Direction_Linear2.x(
531 finaldirection.setY(Direction_Linear2.y(
532 finaldirection.setZ(Direction_Linear2.z(
533 outcomingPhotonPolarization.setX(Polariz
534 outcomingPhotonPolarization.setY(Polariz
535 outcomingPhotonPolarization.setZ(Polariz
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(Polariz
543 outcomingPhotonPolarization.setY(Polariz
544 outcomingPhotonPolarization.setZ(Polariz
545 }
546
547 //Sampling the Final State
548 finaldirection.rotateUz(aDynamicGamma->GetMo
549 fParticleChange->ProposeMomentumDirection(fi
550 fParticleChange->SetProposedKineticEnergy(ph
551 fParticleChange->ProposePolarization(outcomi
552
553 }
554
555 //....oooOO0OOooo........oooOO0OOooo........oo
556
557 G4double G4JAEAPolarizedElasticScatteringModel
558 G4double Sigma_perp,
559 G4double initial_Pol_Plan
560 {
561 G4double phi;
562 G4double phiProbability;
563 G4double Probability=Sigma_perp/(Sigma_para+
564 if (Probability<=G4UniformRand())
565 {
566 do
567 {
568 phi = CLHEP::twopi * G4UniformRand();
569 phiProbability = 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
580 }
581 while (phiProbability < G4UniformRand())
582 }
583 return phi;
584
585 }
586 //....oooOO0OOooo........oooOO0OOooo........oo
587
588 void
589 G4JAEAPolarizedElasticScatteringModel::Initial
590 G4int Z)
591 {
592 G4AutoLock l(&G4JAEAPolarizedElasticScatteri
593 // G4cout << "G4JAEAPolarizedElasticScatter
594 // << Z << G4endl;
595 if(!dataCS[Z]) { ReadData(Z); }
596 l.unlock();
597 }
598
599 //....oooOO0OOooo........oooOO0OOooo........oo
600