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

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

update
File size: 20.7 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.1 2008/10/30 14:16:35 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-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)
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  delete meanFreePathTable;
68  delete crossSectionHandler;
69  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  InitialiseElementSelectors(particle,cuts);
81
82  // Energy limits
83 
84  if (LowEnergyLimit() < lowEnergyLimit)
85    {
86      G4cout << "G4LivermorePolarizedComptonModel: low energy limit increased from " << 
87        LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl;
88      SetLowEnergyLimit(lowEnergyLimit);
89    }
90
91  if (HighEnergyLimit() > highEnergyLimit)
92    {
93      G4cout << "G4LivermorePolarizedComptonModel: high energy limit decreased from " << 
94        HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl;
95      SetHighEnergyLimit(highEnergyLimit);
96    }
97
98  // Reading of data files - all materials are read
99 
100  crossSectionHandler = new G4CrossSectionHandler;
101  crossSectionHandler->Clear();
102  G4String crossSectionFile = "comp/ce-cs-";
103  crossSectionHandler->LoadData(crossSectionFile);
104
105  meanFreePathTable = 0;
106  meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials();
107
108  G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
109  G4String scatterFile = "comp/ce-sf-";
110  scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
111  scatterFunctionData->LoadData(scatterFile);
112
113  // For Doppler broadening
114  shellData.SetOccupancyData();
115  G4String file = "/doppler/shell-doppler";
116  shellData.LoadData(file);
117
118  //
119  if (verboseLevel > 2) 
120    G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl;
121
122  G4cout << "Livermore Polarized Compton model is initialized " << G4endl
123         << "Energy range: "
124         << LowEnergyLimit() / keV << " keV - "
125         << HighEnergyLimit() / GeV << " GeV"
126         << G4endl;
127
128  //
129   
130  if(isInitialised) return;
131
132  if(pParticleChange)
133    fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
134  else
135    fParticleChange = new G4ParticleChangeForGamma();
136   
137  isInitialised = true;
138}
139
140//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
141
142G4double G4LivermorePolarizedComptonModel::ComputeCrossSectionPerAtom(
143                                       const G4ParticleDefinition*,
144                                             G4double GammaEnergy,
145                                             G4double Z, G4double,
146                                             G4double, G4double)
147{
148  if (verboseLevel > 3)
149    G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
150
151  G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy);
152  return cs;
153}
154
155//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
156
157void G4LivermorePolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
158                                              const G4MaterialCutsCouple* couple,
159                                              const G4DynamicParticle* aDynamicGamma,
160                                              G4double,
161                                              G4double)
162{
163  // The scattered gamma energy is sampled according to Klein - Nishina formula.
164  // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
165  // GEANT4 internal units
166  //
167  // Note : Effects due to binding of atomic electrons are negliged.
168
169  if (verboseLevel > 3)
170    G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl;
171
172  G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy();
173  G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
174
175  // Protection: a polarisation parallel to the
176  // direction causes problems;
177  // in that case find a random polarization
178
179  G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
180
181  // Make sure that the polarization vector is perpendicular to the
182  // gamma direction. If not
183
184  if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
185    { // only for testing now
186      gammaPolarization0 = GetRandomPolarization(gammaDirection0);
187    }
188  else
189    {
190      if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
191        {
192          gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
193        }
194    }
195
196  // End of Protection
197
198  // Within energy limit?
199
200  if(gammaEnergy0 <= lowEnergyLimit)
201    {
202      fParticleChange->ProposeTrackStatus(fStopAndKill);
203      fParticleChange->SetProposedKineticEnergy(0.);
204      fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0);
205      // SI - IS THE FOLLOWING RETURN NECESSARY ?
206      return;
207    }
208
209  G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
210
211  // Select randomly one element in the current material
212
213  G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
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      fParticleChange->SetProposedKineticEnergy(0.) ;
432      fParticleChange->ProposeTrackStatus(fStopAndKill);
433    }
434
435  //
436  // kinematic of the scattered electron
437  //
438
439  G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
440
441  // SI - Removed range test
442 
443  G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
444
445  G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
446                                   gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
447
448  fParticleChange->ProposeLocalEnergyDeposit(bindingE);
449 
450  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
451  fvect->push_back(dp);
452
453}
454
455//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
456
457G4double G4LivermorePolarizedComptonModel::SetPhi(G4double energyRate,
458                                             G4double sinSqrTh)
459{
460  G4double rand1;
461  G4double rand2;
462  G4double phiProbability;
463  G4double phi;
464  G4double a, b;
465
466  do
467    {
468      rand1 = G4UniformRand();
469      rand2 = G4UniformRand();
470      phiProbability=0.;
471      phi = twopi*rand1;
472     
473      a = 2*sinSqrTh;
474      b = energyRate + 1/energyRate;
475     
476      phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
477
478     
479 
480    }
481  while ( rand2 > phiProbability );
482  return phi;
483}
484
485
486//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
487
488G4ThreeVector G4LivermorePolarizedComptonModel::SetPerpendicularVector(G4ThreeVector& a)
489{
490  G4double dx = a.x();
491  G4double dy = a.y();
492  G4double dz = a.z();
493  G4double x = dx < 0.0 ? -dx : dx;
494  G4double y = dy < 0.0 ? -dy : dy;
495  G4double z = dz < 0.0 ? -dz : dz;
496  if (x < y) {
497    return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
498  }else{
499    return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
500  }
501}
502
503//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
504
505G4ThreeVector G4LivermorePolarizedComptonModel::GetRandomPolarization(G4ThreeVector& direction0)
506{
507  G4ThreeVector d0 = direction0.unit();
508  G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
509  G4ThreeVector a0 = a1.unit(); // unit vector
510
511  G4double rand1 = G4UniformRand();
512 
513  G4double angle = twopi*rand1; // random polar angle
514  G4ThreeVector b0 = d0.cross(a0); // cross product
515 
516  G4ThreeVector c;
517 
518  c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
519  c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
520  c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
521 
522  G4ThreeVector c0 = c.unit();
523
524  return c0;
525 
526}
527
528//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
529
530G4ThreeVector G4LivermorePolarizedComptonModel::GetPerpendicularPolarization
531(const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
532{
533
534  //
535  // The polarization of a photon is always perpendicular to its momentum direction.
536  // Therefore this function removes those vector component of gammaPolarization, which
537  // points in direction of gammaDirection
538  //
539  // Mathematically we search the projection of the vector a on the plane E, where n is the
540  // plains normal vector.
541  // The basic equation can be found in each geometry book (e.g. Bronstein):
542  // p = a - (a o n)/(n o n)*n
543 
544  return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
545}
546
547//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
548
549G4ThreeVector G4LivermorePolarizedComptonModel::SetNewPolarization(G4double epsilon,
550                                                              G4double sinSqrTh, 
551                                                              G4double phi,
552                                                              G4double costheta) 
553{
554  G4double rand1;
555  G4double rand2;
556  G4double cosPhi = std::cos(phi);
557  G4double sinPhi = std::sin(phi);
558  G4double sinTheta = std::sqrt(sinSqrTh);
559  G4double cosSqrPhi = cosPhi*cosPhi;
560  //  G4double cossqrth = 1.-sinSqrTh;
561  //  G4double sinsqrphi = sinPhi*sinPhi;
562  G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
563 
564
565  // Determination of Theta
566 
567  // ---- MGP ---- Commented out the following 3 lines to avoid compilation
568  // warnings (unused variables)
569  // G4double thetaProbability;
570  G4double theta;
571  // G4double a, b;
572  // G4double cosTheta;
573
574  /*
575
576  depaola method
577 
578  do
579  {
580      rand1 = G4UniformRand();
581      rand2 = G4UniformRand();
582      thetaProbability=0.;
583      theta = twopi*rand1;
584      a = 4*normalisation*normalisation;
585      b = (epsilon + 1/epsilon) - 2;
586      thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
587      cosTheta = std::cos(theta);
588    }
589  while ( rand2 > thetaProbability );
590 
591  G4double cosBeta = cosTheta;
592
593  */
594
595
596  // Dan Xu method (IEEE TNS, 52, 1160 (2005))
597
598  rand1 = G4UniformRand();
599  rand2 = G4UniformRand();
600
601  if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
602    {
603      if (rand2<0.5)
604        theta = pi/2.0;
605      else
606        theta = 3.0*pi/2.0;
607    }
608  else
609    {
610      if (rand2<0.5)
611        theta = 0;
612      else
613        theta = pi;
614    }
615  G4double cosBeta = std::cos(theta);
616  G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
617 
618  G4ThreeVector gammaPolarization1;
619
620  G4double xParallel = normalisation*cosBeta;
621  G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
622  G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
623  G4double xPerpendicular = 0.;
624  G4double yPerpendicular = (costheta)*sinBeta/normalisation;
625  G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
626
627  G4double xTotal = (xParallel + xPerpendicular);
628  G4double yTotal = (yParallel + yPerpendicular);
629  G4double zTotal = (zParallel + zPerpendicular);
630 
631  gammaPolarization1.setX(xTotal);
632  gammaPolarization1.setY(yTotal);
633  gammaPolarization1.setZ(zTotal);
634 
635  return gammaPolarization1;
636
637}
638
639//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
640
641void G4LivermorePolarizedComptonModel::SystemOfRefChange(G4ThreeVector& direction0,
642                                                    G4ThreeVector& direction1,
643                                                    G4ThreeVector& polarization0,
644                                                    G4ThreeVector& polarization1)
645{
646  // direction0 is the original photon direction ---> z
647  // polarization0 is the original photon polarization ---> x
648  // need to specify y axis in the real reference frame ---> y
649  G4ThreeVector Axis_Z0 = direction0.unit();
650  G4ThreeVector Axis_X0 = polarization0.unit();
651  G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
652
653  G4double direction_x = direction1.getX();
654  G4double direction_y = direction1.getY();
655  G4double direction_z = direction1.getZ();
656 
657  direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
658  G4double polarization_x = polarization1.getX();
659  G4double polarization_y = polarization1.getY();
660  G4double polarization_z = polarization1.getZ();
661
662  polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
663
664}
665
666//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
667
668G4double G4LivermorePolarizedComptonModel::GetMeanFreePath(const G4Track& track,
669                                                      G4double,
670                                                      G4ForceCondition*)
671{
672  const G4DynamicParticle* photon = track.GetDynamicParticle();
673  G4double energy = photon->GetKineticEnergy();
674  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
675  size_t materialIndex = couple->GetIndex();
676  G4double meanFreePath;
677  if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
678  else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
679  else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
680  return meanFreePath;
681}
682
683
Note: See TracBrowser for help on using the repository browser.