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