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

Last change on this file was 1347, checked in by garnier, 14 years ago

geant4 tag 9.4
File size: 20.5 KB
Line 
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: G4LivermorePolarizedComptonModel.cc,v 1.6 2009/05/03 08:29:55 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
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
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)
50:G4VEmModel(nam),isInitialised(false),meanFreePathTable(0),scatterFunctionData(0),crossSectionHandler(0)
51{
52  lowEnergyLimit = 250 * eV; 
53  highEnergyLimit = 100 * GeV;
54  //SetLowEnergyLimit(lowEnergyLimit);
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
65  if( verboseLevel>0 ) { 
66  G4cout << "Livermore Polarized Compton is constructed " << G4endl
67         << "Energy range: "
68         << lowEnergyLimit / eV << " eV - "
69         << highEnergyLimit / GeV << " GeV"
70         << G4endl;
71  }
72}
73
74//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
76G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel()
77{ 
78  if (meanFreePathTable)   delete meanFreePathTable;
79  if (crossSectionHandler) delete crossSectionHandler;
80  if (scatterFunctionData) delete scatterFunctionData;
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
91  if (crossSectionHandler)
92  {
93    crossSectionHandler->Clear();
94    delete crossSectionHandler;
95  }
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
120  InitialiseElementSelectors(particle,cuts);
121
122  if(  verboseLevel>0 ) { 
123    G4cout << "Livermore Polarized Compton model is initialized " << G4endl
124         << "Energy range: "
125         << LowEnergyLimit() / eV << " eV - "
126         << HighEnergyLimit() / GeV << " GeV"
127         << G4endl;
128  }
129 
130  //
131   
132  if(isInitialised) return;
133  fParticleChange = GetParticleChangeForGamma();
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
148  if (GammaEnergy < lowEnergyLimit || GammaEnergy > highEnergyLimit) return 0.0;
149
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
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();
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    {
431      gammaEnergy1 = 0.;
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
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 
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.