ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LivermorePolarizedComptonModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LivermorePolarizedComptonModel.cc
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 //
27 // Authors: G.Depaola & F.Longo
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 
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4Electron.hh"
45 #include "G4LossTableManager.hh"
46 #include "G4VAtomDeexcitation.hh"
47 #include "G4AtomicShell.hh"
48 #include "G4Gamma.hh"
49 #include "G4ShellData.hh"
50 #include "G4DopplerProfile.hh"
51 #include "G4Log.hh"
52 #include "G4Exp.hh"
53 #include "G4LogLogInterpolation.hh"
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
56 
57 using namespace std;
58 
64 
65 //static const G4double ln10 = G4Log(10.);
66 
67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
68 
70  :G4VEmModel(nam),isInitialised(false)
71 {
72  verboseLevel= 1;
73  // Verbosity scale:
74  // 0 = nothing
75  // 1 = warning for energy non-conservation
76  // 2 = details of energy budget
77  // 3 = calculation of cross sections, file openings, sampling of atoms
78  // 4 = entering in methods
79 
80  if( verboseLevel>1 )
81  G4cout << "Livermore Polarized Compton is constructed " << G4endl;
82 
83  //Mark this model as "applicable" for atomic deexcitation
84  SetDeexcitationFlag(true);
85 
86  fParticleChange = 0;
88 }
89 
90 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91 
93 {
94  if(IsMaster()) {
95  delete shellData;
96  shellData = 0;
97  delete profileData;
98  profileData = 0;
99  delete scatterFunctionData;
101  for(G4int i=0; i<maxZ; ++i) {
102  if(data[i]) {
103  delete data[i];
104  data[i] = 0;
105  }
106  }
107  }
108 }
109 
110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
111 
113  const G4DataVector& cuts)
114 {
115  if (verboseLevel > 1)
116  G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl;
117 
118  // Initialise element selector
119 
120  if(IsMaster()) {
121 
122  // Access to elements
123 
124  char* path = std::getenv("G4LEDATA");
125 
126  G4ProductionCutsTable* theCoupleTable =
128 
129  G4int numOfCouples = theCoupleTable->GetTableSize();
130 
131  for(G4int i=0; i<numOfCouples; ++i) {
132  const G4Material* material =
133  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
134  const G4ElementVector* theElementVector = material->GetElementVector();
135  G4int nelm = material->GetNumberOfElements();
136 
137  for (G4int j=0; j<nelm; ++j) {
138  G4int Z = G4lrint((*theElementVector)[j]->GetZ());
139  if(Z < 1) { Z = 1; }
140  else if(Z > maxZ){ Z = maxZ; }
141 
142  if( (!data[Z]) ) { ReadData(Z, path); }
143  }
144  }
145 
146  // For Doppler broadening
147  if(!shellData) {
148  shellData = new G4ShellData();
150  G4String file = "/doppler/shell-doppler";
151  shellData->LoadData(file);
152  }
153  if(!profileData) { profileData = new G4DopplerProfile(); }
154 
155  // Scattering Function
156 
158  {
159 
160  G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation;
161  G4String scatterFile = "comp/ce-sf-";
162  scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.);
163  scatterFunctionData->LoadData(scatterFile);
164  }
165 
166  InitialiseElementSelectors(particle, cuts);
167  }
168 
169  if (verboseLevel > 2) {
170  G4cout << "Loaded cross section files" << G4endl;
171  }
172 
173  if( verboseLevel>1 ) {
174  G4cout << "G4LivermoreComptonModel is initialized " << G4endl
175  << "Energy range: "
176  << LowEnergyLimit() / eV << " eV - "
177  << HighEnergyLimit() / GeV << " GeV"
178  << G4endl;
179  }
180  //
181  if(isInitialised) { return; }
182 
185  isInitialised = true;
186 }
187 
188 
190  G4VEmModel* masterModel)
191 {
193 }
194 
195 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
196 
197 void G4LivermorePolarizedComptonModel::ReadData(size_t Z, const char* path)
198 {
199  if (verboseLevel > 1)
200  {
201  G4cout << "G4LivermorePolarizedComptonModel::ReadData()"
202  << G4endl;
203  }
204  if(data[Z]) { return; }
205  const char* datadir = path;
206  if(!datadir)
207  {
208  datadir = std::getenv("G4LEDATA");
209  if(!datadir)
210  {
211  G4Exception("G4LivermorePolarizedComptonModel::ReadData()",
212  "em0006",FatalException,
213  "Environment variable G4LEDATA not defined");
214  return;
215  }
216  }
217 
218  data[Z] = new G4LPhysicsFreeVector();
219 
220  // Activation of spline interpolation
221  data[Z]->SetSpline(false);
222 
223  std::ostringstream ost;
224  ost << datadir << "/livermore/comp/ce-cs-" << Z <<".dat";
225  std::ifstream fin(ost.str().c_str());
226 
227  if( !fin.is_open())
228  {
230  ed << "G4LivermorePolarizedComptonModel data file <" << ost.str().c_str()
231  << "> is not opened!" << G4endl;
232  G4Exception("G4LivermoreComptonModel::ReadData()",
233  "em0003",FatalException,
234  ed,"G4LEDATA version should be G4EMLOW6.34 or later");
235  return;
236  } else {
237  if(verboseLevel > 3) {
238  G4cout << "File " << ost.str()
239  << " is opened by G4LivermorePolarizedComptonModel" << G4endl;
240  }
241  data[Z]->Retrieve(fin, true);
242  data[Z]->ScaleVector(MeV, MeV*barn);
243  }
244  fin.close();
245 }
246 
247 
248 
249 
250 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
251 
253  const G4ParticleDefinition*,
254  G4double GammaEnergy,
257 {
258  if (verboseLevel > 3)
259  G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl;
260 
261  G4double cs = 0.0;
262 
263  if (GammaEnergy < LowEnergyLimit())
264  return 0.0;
265 
266  G4int intZ = G4lrint(Z);
267  if(intZ < 1 || intZ > maxZ) { return cs; }
268 
269  G4LPhysicsFreeVector* pv = data[intZ];
270 
271  // if element was not initialised
272  // do initialisation safely for MT mode
273  if(!pv)
274  {
275  InitialiseForElement(0, intZ);
276  pv = data[intZ];
277  if(!pv) { return cs; }
278  }
279 
280  G4int n = pv->GetVectorLength() - 1;
281  G4double e1 = pv->Energy(0);
282  G4double e2 = pv->Energy(n);
283 
284  if(GammaEnergy <= e1) { cs = GammaEnergy/(e1*e1)*pv->Value(e1); }
285  else if(GammaEnergy <= e2) { cs = pv->Value(GammaEnergy)/GammaEnergy; }
286  else if(GammaEnergy > e2) { cs = pv->Value(e2)/GammaEnergy; }
287 
288  return cs;
289 
290 }
291 
292 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
293 
294 void G4LivermorePolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect,
295  const G4MaterialCutsCouple* couple,
296  const G4DynamicParticle* aDynamicGamma,
297  G4double,
298  G4double)
299 {
300  // The scattered gamma energy is sampled according to Klein - Nishina formula.
301  // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15).
302  // GEANT4 internal units
303  //
304  // Note : Effects due to binding of atomic electrons are negliged.
305 
306  if (verboseLevel > 3)
307  G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl;
308 
309  G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy();
310 
311  // do nothing below the threshold
312  // should never get here because the XS is zero below the limit
313  if (gammaEnergy0 < LowEnergyLimit())
314  return ;
315 
316 
317  G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization();
318 
319  // Protection: a polarisation parallel to the
320  // direction causes problems;
321  // in that case find a random polarization
322 
323  G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection();
324 
325  // Make sure that the polarization vector is perpendicular to the
326  // gamma direction. If not
327 
328  if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0))
329  { // only for testing now
330  gammaPolarization0 = GetRandomPolarization(gammaDirection0);
331  }
332  else
333  {
334  if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0)
335  {
336  gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0);
337  }
338  }
339 
340  // End of Protection
341 
342  G4double E0_m = gammaEnergy0 / electron_mass_c2 ;
343 
344  // Select randomly one element in the current material
345  //G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0);
346  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
347  const G4Element* elm = SelectRandomAtom(couple,particle,gammaEnergy0);
348  G4int Z = (G4int)elm->GetZ();
349 
350  // Sample the energy and the polarization of the scattered photon
351 
352  G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ;
353 
354  G4double epsilon0Local = 1./(1. + 2*E0_m);
355  G4double epsilon0Sq = epsilon0Local*epsilon0Local;
356  G4double alpha1 = - std::log(epsilon0Local);
357  G4double alpha2 = 0.5*(1.- epsilon0Sq);
358 
359  G4double wlGamma = h_Planck*c_light/gammaEnergy0;
360  G4double gammaEnergy1;
361  G4ThreeVector gammaDirection1;
362 
363  do {
364  if ( alpha1/(alpha1+alpha2) > G4UniformRand() )
365  {
366  epsilon = G4Exp(-alpha1*G4UniformRand());
367  epsilonSq = epsilon*epsilon;
368  }
369  else
370  {
371  epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand();
372  epsilon = std::sqrt(epsilonSq);
373  }
374 
375  onecost = (1.- epsilon)/(epsilon*E0_m);
376  sinThetaSqr = onecost*(2.-onecost);
377 
378  // Protection
379  if (sinThetaSqr > 1.)
380  {
381  G4cout
382  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
383  << "sin(theta)**2 = "
384  << sinThetaSqr
385  << "; set to 1"
386  << G4endl;
387  sinThetaSqr = 1.;
388  }
389  if (sinThetaSqr < 0.)
390  {
391  G4cout
392  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
393  << "sin(theta)**2 = "
394  << sinThetaSqr
395  << "; set to 0"
396  << G4endl;
397  sinThetaSqr = 0.;
398  }
399  // End protection
400 
401  G4double x = std::sqrt(onecost/2.) / (wlGamma/cm);;
402  G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1);
403  greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction;
404 
405  } while(greject < G4UniformRand()*Z);
406 
407 
408  // ****************************************************
409  // Phi determination
410  // ****************************************************
411 
412  G4double phi = SetPhi(epsilon,sinThetaSqr);
413 
414  //
415  // scattered gamma angles. ( Z - axis along the parent gamma)
416  //
417 
418  G4double cosTheta = 1. - onecost;
419 
420  // Protection
421 
422  if (cosTheta > 1.)
423  {
424  G4cout
425  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
426  << "cosTheta = "
427  << cosTheta
428  << "; set to 1"
429  << G4endl;
430  cosTheta = 1.;
431  }
432  if (cosTheta < -1.)
433  {
434  G4cout
435  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
436  << "cosTheta = "
437  << cosTheta
438  << "; set to -1"
439  << G4endl;
440  cosTheta = -1.;
441  }
442  // End protection
443 
444 
445  G4double sinTheta = std::sqrt (sinThetaSqr);
446 
447  // Protection
448  if (sinTheta > 1.)
449  {
450  G4cout
451  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
452  << "sinTheta = "
453  << sinTheta
454  << "; set to 1"
455  << G4endl;
456  sinTheta = 1.;
457  }
458  if (sinTheta < -1.)
459  {
460  G4cout
461  << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries "
462  << "sinTheta = "
463  << sinTheta
464  << "; set to -1"
465  << G4endl;
466  sinTheta = -1.;
467  }
468  // End protection
469 
470 
471  G4double dirx = sinTheta*std::cos(phi);
472  G4double diry = sinTheta*std::sin(phi);
473  G4double dirz = cosTheta ;
474 
475 
476  // oneCosT , eom
477 
478  // Doppler broadening - Method based on:
479  // Y. Namito, S. Ban and H. Hirayama,
480  // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code"
481  // NIM A 349, pp. 489-494, 1994
482 
483  // Maximum number of sampling iterations
484 
485  static G4int maxDopplerIterations = 1000;
486  G4double bindingE = 0.;
487  G4double photonEoriginal = epsilon * gammaEnergy0;
488  G4double photonE = -1.;
489  G4int iteration = 0;
490  G4double eMax = gammaEnergy0;
491 
492  G4int shellIdx = 0;
493 
494  if (verboseLevel > 3) {
495  G4cout << "Started loop to sample broading" << G4endl;
496  }
497 
498  do
499  {
500  iteration++;
501  // Select shell based on shell occupancy
502  shellIdx = shellData->SelectRandomShell(Z);
503  bindingE = shellData->BindingEnergy(Z,shellIdx);
504 
505  if (verboseLevel > 3) {
506  G4cout << "Shell ID= " << shellIdx
507  << " Ebind(keV)= " << bindingE/keV << G4endl;
508  }
509  eMax = gammaEnergy0 - bindingE;
510 
511  // Randomly sample bound electron momentum (memento: the data set is in Atomic Units)
512  G4double pSample = profileData->RandomSelectMomentum(Z,shellIdx);
513 
514  if (verboseLevel > 3) {
515  G4cout << "pSample= " << pSample << G4endl;
516  }
517  // Rescale from atomic units
518  G4double pDoppler = pSample * fine_structure_const;
519  G4double pDoppler2 = pDoppler * pDoppler;
520  G4double var2 = 1. + onecost * E0_m;
521  G4double var3 = var2*var2 - pDoppler2;
522  G4double var4 = var2 - pDoppler2 * cosTheta;
523  G4double var = var4*var4 - var3 + pDoppler2 * var3;
524  if (var > 0.)
525  {
526  G4double varSqrt = std::sqrt(var);
527  G4double scale = gammaEnergy0 / var3;
528  // Random select either root
529  if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale;
530  else photonE = (var4 + varSqrt) * scale;
531  }
532  else
533  {
534  photonE = -1.;
535  }
536  } while ( iteration <= maxDopplerIterations &&
537  (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) );
538  //while (iteration <= maxDopplerIterations && photonE > eMax); ???
539 
540 
541  // End of recalculation of photon energy with Doppler broadening
542  // Revert to original if maximum number of iterations threshold has been reached
543  if (iteration >= maxDopplerIterations)
544  {
545  photonE = photonEoriginal;
546  bindingE = 0.;
547  }
548 
549  gammaEnergy1 = photonE;
550 
551  //
552  // update G4VParticleChange for the scattered photon
553  //
554 
555 
556 
557  // gammaEnergy1 = epsilon*gammaEnergy0;
558 
559 
560  // New polarization
561 
562  G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon,
563  sinThetaSqr,
564  phi,
565  cosTheta);
566 
567  // Set new direction
568  G4ThreeVector tmpDirection1( dirx,diry,dirz );
569  gammaDirection1 = tmpDirection1;
570 
571  // Change reference frame.
572 
573  SystemOfRefChange(gammaDirection0,gammaDirection1,
574  gammaPolarization0,gammaPolarization1);
575 
576  if (gammaEnergy1 > 0.)
577  {
578  fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ;
579  fParticleChange->ProposeMomentumDirection( gammaDirection1 );
580  fParticleChange->ProposePolarization( gammaPolarization1 );
581  }
582  else
583  {
584  gammaEnergy1 = 0.;
587  }
588 
589  //
590  // kinematic of the scattered electron
591  //
592 
593  G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE;
594 
595  // SI -protection against negative final energy: no e- is created
596  // like in G4LivermoreComptonModel.cc
597  if(ElecKineEnergy < 0.0) {
598  fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0 - gammaEnergy1);
599  return;
600  }
601 
602  // SI - Removed range test
603 
604  G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2));
605 
606  G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 -
607  gammaEnergy1 * gammaDirection1) * (1./ElecMomentum));
608 
609  G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ;
610  fvect->push_back(dp);
611 
612  // sample deexcitation
613  //
614 
615  if (verboseLevel > 3) {
616  G4cout << "Started atomic de-excitation " << fAtomDeexcitation << G4endl;
617  }
618 
619  if(fAtomDeexcitation && iteration < maxDopplerIterations) {
620  G4int index = couple->GetIndex();
622  size_t nbefore = fvect->size();
624  const G4AtomicShell* shell = fAtomDeexcitation->GetAtomicShell(Z, as);
625  fAtomDeexcitation->GenerateParticles(fvect, shell, Z, index);
626  size_t nafter = fvect->size();
627  if(nafter > nbefore) {
628  for (size_t i=nbefore; i<nafter; ++i) {
629  //Check if there is enough residual energy
630  if (bindingE >= ((*fvect)[i])->GetKineticEnergy())
631  {
632  //Ok, this is a valid secondary: keep it
633  bindingE -= ((*fvect)[i])->GetKineticEnergy();
634  }
635  else
636  {
637  //Invalid secondary: not enough energy to create it!
638  //Keep its energy in the local deposit
639  delete (*fvect)[i];
640  (*fvect)[i]=0;
641  }
642  }
643  }
644  }
645  }
646  //This should never happen
647  if(bindingE < 0.0)
648  G4Exception("G4LivermoreComptonModel::SampleSecondaries()",
649  "em2050",FatalException,"Negative local energy deposit");
650 
652 
653 }
654 
655 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
656 
658  G4double sinSqrTh)
659 {
660  G4double rand1;
661  G4double rand2;
662  G4double phiProbability;
663  G4double phi;
664  G4double a, b;
665 
666  do
667  {
668  rand1 = G4UniformRand();
669  rand2 = G4UniformRand();
670  phiProbability=0.;
671  phi = twopi*rand1;
672 
673  a = 2*sinSqrTh;
674  b = energyRate + 1/energyRate;
675 
676  phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi));
677 
678 
679 
680  }
681  while ( rand2 > phiProbability );
682  return phi;
683 }
684 
685 
686 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
687 
689 {
690  G4double dx = a.x();
691  G4double dy = a.y();
692  G4double dz = a.z();
693  G4double x = dx < 0.0 ? -dx : dx;
694  G4double y = dy < 0.0 ? -dy : dy;
695  G4double z = dz < 0.0 ? -dz : dz;
696  if (x < y) {
697  return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy);
698  }else{
699  return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0);
700  }
701 }
702 
703 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
704 
706 {
707  G4ThreeVector d0 = direction0.unit();
708  G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal
709  G4ThreeVector a0 = a1.unit(); // unit vector
710 
712 
713  G4double angle = twopi*rand1; // random polar angle
714  G4ThreeVector b0 = d0.cross(a0); // cross product
715 
717 
718  c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x());
719  c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y());
720  c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z());
721 
722  G4ThreeVector c0 = c.unit();
723 
724  return c0;
725 
726 }
727 
728 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
729 
731 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const
732 {
733 
734  //
735  // The polarization of a photon is always perpendicular to its momentum direction.
736  // Therefore this function removes those vector component of gammaPolarization, which
737  // points in direction of gammaDirection
738  //
739  // Mathematically we search the projection of the vector a on the plane E, where n is the
740  // plains normal vector.
741  // The basic equation can be found in each geometry book (e.g. Bronstein):
742  // p = a - (a o n)/(n o n)*n
743 
744  return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection;
745 }
746 
747 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
748 
750  G4double sinSqrTh,
751  G4double phi,
752  G4double costheta)
753 {
754  G4double rand1;
755  G4double rand2;
756  G4double cosPhi = std::cos(phi);
757  G4double sinPhi = std::sin(phi);
758  G4double sinTheta = std::sqrt(sinSqrTh);
759  G4double cosSqrPhi = cosPhi*cosPhi;
760  // G4double cossqrth = 1.-sinSqrTh;
761  // G4double sinsqrphi = sinPhi*sinPhi;
762  G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh);
763 
764 
765  // Determination of Theta
766 
767  // ---- MGP ---- Commented out the following 3 lines to avoid compilation
768  // warnings (unused variables)
769  // G4double thetaProbability;
770  G4double theta;
771  // G4double a, b;
772  // G4double cosTheta;
773 
774  /*
775 
776  depaola method
777 
778  do
779  {
780  rand1 = G4UniformRand();
781  rand2 = G4UniformRand();
782  thetaProbability=0.;
783  theta = twopi*rand1;
784  a = 4*normalisation*normalisation;
785  b = (epsilon + 1/epsilon) - 2;
786  thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b);
787  cosTheta = std::cos(theta);
788  }
789  while ( rand2 > thetaProbability );
790 
791  G4double cosBeta = cosTheta;
792 
793  */
794 
795 
796  // Dan Xu method (IEEE TNS, 52, 1160 (2005))
797 
798  rand1 = G4UniformRand();
799  rand2 = G4UniformRand();
800 
801  if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi))
802  {
803  if (rand2<0.5)
804  theta = pi/2.0;
805  else
806  theta = 3.0*pi/2.0;
807  }
808  else
809  {
810  if (rand2<0.5)
811  theta = 0;
812  else
813  theta = pi;
814  }
815  G4double cosBeta = std::cos(theta);
816  G4double sinBeta = std::sqrt(1-cosBeta*cosBeta);
817 
818  G4ThreeVector gammaPolarization1;
819 
820  G4double xParallel = normalisation*cosBeta;
821  G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation;
822  G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation;
823  G4double xPerpendicular = 0.;
824  G4double yPerpendicular = (costheta)*sinBeta/normalisation;
825  G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation;
826 
827  G4double xTotal = (xParallel + xPerpendicular);
828  G4double yTotal = (yParallel + yPerpendicular);
829  G4double zTotal = (zParallel + zPerpendicular);
830 
831  gammaPolarization1.setX(xTotal);
832  gammaPolarization1.setY(yTotal);
833  gammaPolarization1.setZ(zTotal);
834 
835  return gammaPolarization1;
836 
837 }
838 
839 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
840 
842  G4ThreeVector& direction1,
843  G4ThreeVector& polarization0,
844  G4ThreeVector& polarization1)
845 {
846  // direction0 is the original photon direction ---> z
847  // polarization0 is the original photon polarization ---> x
848  // need to specify y axis in the real reference frame ---> y
849  G4ThreeVector Axis_Z0 = direction0.unit();
850  G4ThreeVector Axis_X0 = polarization0.unit();
851  G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed;
852 
853  G4double direction_x = direction1.getX();
854  G4double direction_y = direction1.getY();
855  G4double direction_z = direction1.getZ();
856 
857  direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit();
858  G4double polarization_x = polarization1.getX();
859  G4double polarization_y = polarization1.getY();
860  G4double polarization_z = polarization1.getZ();
861 
862  polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit();
863 
864 }
865 
866 
867 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
868 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
869 
870 #include "G4AutoLock.hh"
871 namespace { G4Mutex LivermorePolarizedComptonModelMutex = G4MUTEX_INITIALIZER; }
872 
873 void
875  G4int Z)
876 {
877  G4AutoLock l(&LivermorePolarizedComptonModelMutex);
878  // G4cout << "G4LivermoreComptonModel::InitialiseForElement Z= "
879  // << Z << G4endl;
880  if(!data[Z]) { ReadData(Z); }
881  l.unlock();
882 }