source: trunk/source/processes/electromagnetic/lowenergy/src/G4LivermorePolarizedComptonModel.cc@ 1229

Last change on this file since 1229 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
File size: 20.5 KB
RevLine 
[968]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//
[1055]26// $Id: G4LivermorePolarizedComptonModel.cc,v 1.6 2009/05/03 08:29:55 sincerti Exp $
[1228]27// GEANT4 tag $Name: geant4-09-03 $
[968]28//
[1055]29// History:
30// --------
31// 02 May 2009 S Incerti as V. Ivanchenko proposed in G4LivermoreComptonModel.cc
32//
33// Cleanup initialisation and generation of secondaries:
34// - apply internal high-energy limit only in constructor
35// - do not apply low-energy limit (default is 0)
36// - remove GetMeanFreePath method and table
37// - added protection against numerical problem in energy sampling
38// - use G4ElementSelector
[968]39
40#include "G4LivermorePolarizedComptonModel.hh"
41
42//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
43
44using namespace std;
45
46//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
47
48G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel(const G4ParticleDefinition*,
49 const G4String& nam)
[1055]50:G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),scatterFunctionData(0),crossSectionHandler(0)
[968]51{
[1055]52 lowEnergyLimit = 250 * eV;
[968]53 highEnergyLimit = 100 * GeV;
[1055]54 //SetLowEnergyLimit(lowEnergyLimit);
[968]55 SetHighEnergyLimit(highEnergyLimit);
56
57 verboseLevel= 0;
58 // Verbosity scale:
59 // 0 = nothing
60 // 1 = warning for energy non-conservation
61 // 2 = details of energy budget
62 // 3 = calculation of cross sections, file openings, sampling of atoms
63 // 4 = entering in methods
64
[1055]65 if( verboseLevel>0 ) {
[968]66 G4cout << "Livermore Polarized Compton is constructed " << G4endl
67 << "Energy range: "
[1055]68 << lowEnergyLimit / eV << " eV - "
[968]69 << highEnergyLimit / GeV << " GeV"
70 << G4endl;
[1055]71 }
[968]72}
73
74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
76G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel()
77{
[1055]78 if (meanFreePathTable) delete meanFreePathTable;
79 if (crossSectionHandler) delete crossSectionHandler;
80 if (scatterFunctionData) delete scatterFunctionData;
[968]81}
82
83//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
84
85void G4LivermorePolarizedComptonModel::Initialise(const G4ParticleDefinition* particle,
86 const G4DataVector& cuts)
87{
88 if (verboseLevel > 3)
89 G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl;
90
[1055]91 if (crossSectionHandler)
92 {
93 crossSectionHandler->Clear();
94 delete crossSectionHandler;
95 }
[968]96
97 // Reading of data files - all materials are read
98
99 crossSectionHandler = new G4CrossSectionHandler;
100 crossSectionHandler->Clear();
101 G4String crossSectionFile = "comp/ce-cs-";
102 crossSectionHandler->LoadData(crossSectionFile);
103
104 meanFreePathTable = 0;
105 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
106
107 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
108 G4String scatterFile = "comp/ce-sf-";
109 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
110 scatterFunctionData->LoadData(scatterFile);
111
112 // For Doppler broadening
113 shellData.SetOccupancyData();
114 G4String file = "/doppler/shell-doppler";
115 shellData.LoadData(file);
116
117 if (verboseLevel > 2)
118 G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl;
119
[1055]120 InitialiseElementSelectors(particle,cuts);
121
122 if( verboseLevel>0 ) {
123 G4cout << "Livermore Polarized Compton model is initialized " << G4endl
[968]124 << "Energy range: "
[1055]125 << LowEnergyLimit() / eV << " eV - "
[968]126 << HighEnergyLimit() / GeV << " GeV"
127 << G4endl;
[1055]128 }
129
[968]130 //
131
132 if(isInitialised) return;
[1055]133 fParticleChange = GetParticleChangeForGamma();
[968]134 isInitialised = true;
135}
136
137//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
138
139G4double G4LivermorePolarizedComptonModel::ComputeCrossSectionPerAtom(
140 const G4ParticleDefinition*,
141 G4double GammaEnergy,
142 G4double Z, G4double,
143 G4double, G4double)
144{
145 if (verboseLevel > 3)
146 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
147
[1055]148 if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
149
[968]150 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
151 return cs;
152}
153
154//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
155
156void G4LivermorePolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
157 const G4MaterialCutsCouple* couple,
158 const G4DynamicParticle* aDynamicGamma,
159 G4double,
160 G4double)
161{
162 // The scattered gamma energy is sampled according to Klein - Nishina formula.
163 // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
164 // GEANT4 internal units
165 //
166 // Note : Effects due to binding of atomic electrons are negliged.
167
168 if (verboseLevel > 3)
169 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl;
170
171 G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy();
172 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
173
174 // Protection: a polarisation parallel to the
175 // direction causes problems;
176 // in that case find a random polarization
177
178 G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
179
180 // Make sure that the polarization vector is perpendicular to the
181 // gamma direction. If not
182
183 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
184 { // only for testing now
185 gammaPolarization0 = GetRandomPolarization(gammaDirection0);
186 }
187 else
188 {
189 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
190 {
191 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
192 }
193 }
194
195 // End of Protection
196
197 // Within energy limit?
198
199 if(gammaEnergy0 <= lowEnergyLimit)
200 {
201 fParticleChange->ProposeTrackStatus(fStopAndKill);
202 fParticleChange->SetProposedKineticEnergy(0.);
203 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0);
204 return;
205 }
206
207 G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
208
209 // Select randomly one element in the current material
[1055]210 //G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
211 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
212 const G4Element* elm = SelectRandomAtom(couple,particle,gammaEnergy0);
213 G4int Z = (G4int)elm->GetZ();
[968]214
215 // Sample the energy and the polarization of the scattered photon
216
217 G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
218
219 G4double epsilon0 = 1./(1. + 2*E0_m);
220 G4double epsilon0Sq = epsilon0*epsilon0;
221 G4double alpha1 = - std::log(epsilon0);
222 G4double alpha2 = 0.5*(1.- epsilon0Sq);
223
224 G4double wlGamma = h_Planck*c_light/gammaEnergy0;
225 G4double gammaEnergy1;
226 G4ThreeVector gammaDirection1;
227
228 do {
229 if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
230 {
231 epsilon = std::exp(-alpha1*G4UniformRand());
232 epsilonSq = epsilon*epsilon;
233 }
234 else
235 {
236 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
237 epsilon = std::sqrt(epsilonSq);
238 }
239
240 onecost = (1.- epsilon)/(epsilon*E0_m);
241 sinThetaSqr = onecost*(2.-onecost);
242
243 // Protection
244 if (sinThetaSqr > 1.)
245 {
246 G4cout
247 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
248 << "sin(theta)**2 = "
249 << sinThetaSqr
250 << "; set to 1"
251 << G4endl;
252 sinThetaSqr = 1.;
253 }
254 if (sinThetaSqr < 0.)
255 {
256 G4cout
257 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
258 << "sin(theta)**2 = "
259 << sinThetaSqr
260 << "; set to 0"
261 << G4endl;
262 sinThetaSqr = 0.;
263 }
264 // End protection
265
266 G4double x = std::sqrt(onecost/2.) / (wlGamma/cm);;
267 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
268 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
269
270 } while(greject < G4UniformRand()*Z);
271
272
273 // ****************************************************
274 // Phi determination
275 // ****************************************************
276
277 G4double phi = SetPhi(epsilon,sinThetaSqr);
278
279 //
280 // scattered gamma angles. ( Z - axis along the parent gamma)
281 //
282
283 G4double cosTheta = 1. - onecost;
284
285 // Protection
286
287 if (cosTheta > 1.)
288 {
289 G4cout
290 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
291 << "cosTheta = "
292 << cosTheta
293 << "; set to 1"
294 << G4endl;
295 cosTheta = 1.;
296 }
297 if (cosTheta < -1.)
298 {
299 G4cout
300 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
301 << "cosTheta = "
302 << cosTheta
303 << "; set to -1"
304 << G4endl;
305 cosTheta = -1.;
306 }
307 // End protection
308
309
310 G4double sinTheta = std::sqrt (sinThetaSqr);
311
312 // Protection
313 if (sinTheta > 1.)
314 {
315 G4cout
316 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
317 << "sinTheta = "
318 << sinTheta
319 << "; set to 1"
320 << G4endl;
321 sinTheta = 1.;
322 }
323 if (sinTheta < -1.)
324 {
325 G4cout
326 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
327 << "sinTheta = "
328 << sinTheta
329 << "; set to -1"
330 << G4endl;
331 sinTheta = -1.;
332 }
333 // End protection
334
335
336 G4double dirx = sinTheta*std::cos(phi);
337 G4double diry = sinTheta*std::sin(phi);
338 G4double dirz = cosTheta ;
339
340
341 // oneCosT , eom
342
343 // Doppler broadening - Method based on:
344 // Y. Namito, S. Ban and H. Hirayama,
345 // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
346 // NIM A 349, pp. 489-494, 1994
347
348 // Maximum number of sampling iterations
349
350 G4int maxDopplerIterations = 1000;
351 G4double bindingE = 0.;
352 G4double photonEoriginal = epsilon * gammaEnergy0;
353 G4double photonE = -1.;
354 G4int iteration = 0;
355 G4double eMax = gammaEnergy0;
356
357 do
358 {
359 iteration++;
360 // Select shell based on shell occupancy
361 G4int shell = shellData.SelectRandomShell(Z);
362 bindingE = shellData.BindingEnergy(Z,shell);
363
364 eMax = gammaEnergy0 - bindingE;
365
366 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
367 G4double pSample = profileData.RandomSelectMomentum(Z,shell);
368 // Rescale from atomic units
369 G4double pDoppler = pSample * fine_structure_const;
370 G4double pDoppler2 = pDoppler * pDoppler;
371 G4double var2 = 1. + onecost * E0_m;
372 G4double var3 = var2*var2 - pDoppler2;
373 G4double var4 = var2 - pDoppler2 * cosTheta;
374 G4double var = var4*var4 - var3 + pDoppler2 * var3;
375 if (var > 0.)
376 {
377 G4double varSqrt = std::sqrt(var);
378 G4double scale = gammaEnergy0 / var3;
379 // Random select either root
380 if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
381 else photonE = (var4 + varSqrt) * scale;
382 }
383 else
384 {
385 photonE = -1.;
386 }
387 } while ( iteration <= maxDopplerIterations &&
388 (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
389
390 // End of recalculation of photon energy with Doppler broadening
391 // Revert to original if maximum number of iterations threshold has been reached
392 if (iteration >= maxDopplerIterations)
393 {
394 photonE = photonEoriginal;
395 bindingE = 0.;
396 }
397
398 gammaEnergy1 = photonE;
399
400 //
401 // update G4VParticleChange for the scattered photon
402 //
403
404 // gammaEnergy1 = epsilon*gammaEnergy0;
405
406
407 // New polarization
408
409 G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
410 sinThetaSqr,
411 phi,
412 cosTheta);
413
414 // Set new direction
415 G4ThreeVector tmpDirection1( dirx,diry,dirz );
416 gammaDirection1 = tmpDirection1;
417
418 // Change reference frame.
419
420 SystemOfRefChange(gammaDirection0,gammaDirection1,
421 gammaPolarization0,gammaPolarization1);
422
423 if (gammaEnergy1 > 0.)
424 {
425 fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ;
426 fParticleChange->ProposeMomentumDirection( gammaDirection1 );
427 fParticleChange->ProposePolarization( gammaPolarization1 );
428 }
429 else
430 {
[1055]431 gammaEnergy1 = 0.;
[968]432 fParticleChange->SetProposedKineticEnergy(0.) ;
433 fParticleChange->ProposeTrackStatus(fStopAndKill);
434 }
435
436 //
437 // kinematic of the scattered electron
438 //
439
440 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
441
[1055]442 // SI -protection against negative final energy: no e- is created
443 // like in G4LivermoreComptonModel.cc
444 if(ElecKineEnergy < 0.0) {
445 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0 - gammaEnergy1);
446 return;
447 }
448
[968]449 // SI - Removed range test
450
451 G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
452
453 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
454 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
455
456 fParticleChange->ProposeLocalEnergyDeposit(bindingE);
457
458 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
459 fvect->push_back(dp);
460
461}
462
463//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
464
465G4double G4LivermorePolarizedComptonModel::SetPhi(G4double energyRate,
466 G4double sinSqrTh)
467{
468 G4double rand1;
469 G4double rand2;
470 G4double phiProbability;
471 G4double phi;
472 G4double a, b;
473
474 do
475 {
476 rand1 = G4UniformRand();
477 rand2 = G4UniformRand();
478 phiProbability=0.;
479 phi = twopi*rand1;
480
481 a = 2*sinSqrTh;
482 b = energyRate + 1/energyRate;
483
484 phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
485
486
487
488 }
489 while ( rand2 > phiProbability );
490 return phi;
491}
492
493
494//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
495
496G4ThreeVector G4LivermorePolarizedComptonModel::SetPerpendicularVector(G4ThreeVector& a)
497{
498 G4double dx = a.x();
499 G4double dy = a.y();
500 G4double dz = a.z();
501 G4double x = dx < 0.0 ? -dx : dx;
502 G4double y = dy < 0.0 ? -dy : dy;
503 G4double z = dz < 0.0 ? -dz : dz;
504 if (x < y) {
505 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
506 }else{
507 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
508 }
509}
510
511//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
512
513G4ThreeVector G4LivermorePolarizedComptonModel::GetRandomPolarization(G4ThreeVector& direction0)
514{
515 G4ThreeVector d0 = direction0.unit();
516 G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
517 G4ThreeVector a0 = a1.unit(); // unit vector
518
519 G4double rand1 = G4UniformRand();
520
521 G4double angle = twopi*rand1; // random polar angle
522 G4ThreeVector b0 = d0.cross(a0); // cross product
523
524 G4ThreeVector c;
525
526 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
527 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
528 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
529
530 G4ThreeVector c0 = c.unit();
531
532 return c0;
533
534}
535
536//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
537
538G4ThreeVector G4LivermorePolarizedComptonModel::GetPerpendicularPolarization
539(const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
540{
541
542 //
543 // The polarization of a photon is always perpendicular to its momentum direction.
544 // Therefore this function removes those vector component of gammaPolarization, which
545 // points in direction of gammaDirection
546 //
547 // Mathematically we search the projection of the vector a on the plane E, where n is the
548 // plains normal vector.
549 // The basic equation can be found in each geometry book (e.g. Bronstein):
550 // p = a - (a o n)/(n o n)*n
551
552 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
553}
554
555//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
556
557G4ThreeVector G4LivermorePolarizedComptonModel::SetNewPolarization(G4double epsilon,
558 G4double sinSqrTh,
559 G4double phi,
560 G4double costheta)
561{
562 G4double rand1;
563 G4double rand2;
564 G4double cosPhi = std::cos(phi);
565 G4double sinPhi = std::sin(phi);
566 G4double sinTheta = std::sqrt(sinSqrTh);
567 G4double cosSqrPhi = cosPhi*cosPhi;
568 // G4double cossqrth = 1.-sinSqrTh;
569 // G4double sinsqrphi = sinPhi*sinPhi;
570 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
571
572
573 // Determination of Theta
574
575 // ---- MGP ---- Commented out the following 3 lines to avoid compilation
576 // warnings (unused variables)
577 // G4double thetaProbability;
578 G4double theta;
579 // G4double a, b;
580 // G4double cosTheta;
581
582 /*
583
584 depaola method
585
586 do
587 {
588 rand1 = G4UniformRand();
589 rand2 = G4UniformRand();
590 thetaProbability=0.;
591 theta = twopi*rand1;
592 a = 4*normalisation*normalisation;
593 b = (epsilon + 1/epsilon) - 2;
594 thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
595 cosTheta = std::cos(theta);
596 }
597 while ( rand2 > thetaProbability );
598
599 G4double cosBeta = cosTheta;
600
601 */
602
603
604 // Dan Xu method (IEEE TNS, 52, 1160 (2005))
605
606 rand1 = G4UniformRand();
607 rand2 = G4UniformRand();
608
609 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
610 {
611 if (rand2<0.5)
612 theta = pi/2.0;
613 else
614 theta = 3.0*pi/2.0;
615 }
616 else
617 {
618 if (rand2<0.5)
619 theta = 0;
620 else
621 theta = pi;
622 }
623 G4double cosBeta = std::cos(theta);
624 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
625
626 G4ThreeVector gammaPolarization1;
627
628 G4double xParallel = normalisation*cosBeta;
629 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
630 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
631 G4double xPerpendicular = 0.;
632 G4double yPerpendicular = (costheta)*sinBeta/normalisation;
633 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
634
635 G4double xTotal = (xParallel + xPerpendicular);
636 G4double yTotal = (yParallel + yPerpendicular);
637 G4double zTotal = (zParallel + zPerpendicular);
638
639 gammaPolarization1.setX(xTotal);
640 gammaPolarization1.setY(yTotal);
641 gammaPolarization1.setZ(zTotal);
642
643 return gammaPolarization1;
644
645}
646
647//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
648
649void G4LivermorePolarizedComptonModel::SystemOfRefChange(G4ThreeVector& direction0,
650 G4ThreeVector& direction1,
651 G4ThreeVector& polarization0,
652 G4ThreeVector& polarization1)
653{
654 // direction0 is the original photon direction ---> z
655 // polarization0 is the original photon polarization ---> x
656 // need to specify y axis in the real reference frame ---> y
657 G4ThreeVector Axis_Z0 = direction0.unit();
658 G4ThreeVector Axis_X0 = polarization0.unit();
659 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
660
661 G4double direction_x = direction1.getX();
662 G4double direction_y = direction1.getY();
663 G4double direction_z = direction1.getZ();
664
665 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
666 G4double polarization_x = polarization1.getX();
667 G4double polarization_y = polarization1.getY();
668 G4double polarization_z = polarization1.getZ();
669
670 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
671
672}
673
674
Note: See TracBrowser for help on using the repository browser.