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

Last change on this file since 968 was 968, checked in by garnier, 15 years ago

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