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