ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LivermorePolarizedRayleighModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LivermorePolarizedRayleighModel.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 // Author: Sebastien Incerti
28 // 30 October 2008
29 // on base of G4LowEnergyPolarizedRayleigh developed by R. Capra
30 //
31 // History:
32 // --------
33 // 02 May 2009 S Incerti as V. Ivanchenko proposed in G4LivermoreRayleighModel.cc
34 //
35 // Cleanup initialisation and generation of secondaries:
36 // - apply internal high-energy limit only in constructor
37 // - do not apply low-energy limit (default is 0)
38 // - remove GetMeanFreePath method and table
39 // - remove initialisation of element selector
40 // - use G4ElementSelector
41 
43 #include "G4PhysicalConstants.hh"
44 #include "G4SystemOfUnits.hh"
45 #include "G4LogLogInterpolation.hh"
46 #include "G4CompositeEMDataSet.hh"
47 
48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
49 
50 using namespace std;
51 
52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53 
57 
59  const G4String& nam)
60  :G4VEmModel(nam),fParticleChange(0),isInitialised(false)
61 {
62  fParticleChange =0;
63  lowEnergyLimit = 250 * eV;
64  //SetLowEnergyLimit(lowEnergyLimit);
65  //SetHighEnergyLimit(highEnergyLimit);
66  //
67  verboseLevel= 0;
68  // Verbosity scale:
69  // 0 = nothing
70  // 1 = warning for energy non-conservation
71  // 2 = details of energy budget
72  // 3 = calculation of cross sections, file openings, sampling of atoms
73  // 4 = entering in methods
74 
75  if(verboseLevel > 0) {
76  G4cout << "Livermore Polarized Rayleigh is constructed " << G4endl
77  << "Energy range: "
78  << LowEnergyLimit() / eV << " eV - "
79  << HighEnergyLimit() / GeV << " GeV"
80  << G4endl;
81  }
82 }
83 
84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 
87 {
88  if(IsMaster()) {
89  for(G4int i=0; i<maxZ; ++i) {
90  if(dataCS[i]) {
91  delete dataCS[i];
92  dataCS[i] = 0;
93  }
94  }
95  delete formFactorData;
96  formFactorData = 0;
97 
98  }
99 }
100 
101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
102 
104  const G4DataVector& cuts)
105 {
106 // Rayleigh process: The Quantum Theory of Radiation
107 // W. Heitler, Oxford at the Clarendon Press, Oxford (1954)
108 // Scattering function: A simple model of photon transport
109 // D.E. Cullen, Nucl. Instr. Meth. in Phys. Res. B 101 (1995) 499-510
110 // Polarization of the outcoming photon: Beam test of a prototype detector array for the PoGO astronomical hard X-ray/soft gamma-ray polarimeter
111 // T. Mizuno et al., Nucl. Instr. Meth. in Phys. Res. A 540 (2005) 158-168
112 
113  if (verboseLevel > 3)
114  G4cout << "Calling G4LivermorePolarizedRayleighModel::Initialise()" << G4endl;
115 
116 
117  if(IsMaster()) {
118 
119  // Form Factor
120 
121  G4VDataSetAlgorithm* ffInterpolation = new G4LogLogInterpolation;
122  G4String formFactorFile = "rayl/re-ff-";
123  formFactorData = new G4CompositeEMDataSet(ffInterpolation,1.,1.);
124  formFactorData->LoadData(formFactorFile);
125 
126  // Initialise element selector
127  InitialiseElementSelectors(particle, cuts);
128 
129  // Access to elements
130  char* path = std::getenv("G4LEDATA");
131  G4ProductionCutsTable* theCoupleTable =
133  G4int numOfCouples = theCoupleTable->GetTableSize();
134 
135  for(G4int i=0; i<numOfCouples; ++i)
136  {
137  const G4MaterialCutsCouple* couple =
138  theCoupleTable->GetMaterialCutsCouple(i);
139  const G4Material* material = couple->GetMaterial();
140  const G4ElementVector* theElementVector = material->GetElementVector();
141  G4int nelm = material->GetNumberOfElements();
142 
143  for (G4int j=0; j<nelm; ++j)
144  {
145  G4int Z = G4lrint((*theElementVector)[j]->GetZ());
146  if(Z < 1) { Z = 1; }
147  else if(Z > maxZ) { Z = maxZ; }
148  if( (!dataCS[Z]) ) { ReadData(Z, path); }
149  }
150  }
151  }
152 
153  if(isInitialised) { return; }
155  isInitialised = true;
156 
157 }
158 
159 
160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
161 
163  G4VEmModel* masterModel)
164 {
166 }
167 
168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
169 
170 void G4LivermorePolarizedRayleighModel::ReadData(size_t Z, const char* path)
171 {
172  if (verboseLevel > 1)
173  {
174  G4cout << "Calling ReadData() of G4LivermoreRayleighModel"
175  << G4endl;
176  }
177 
178  if(dataCS[Z]) { return; }
179 
180  const char* datadir = path;
181 
182  if(!datadir)
183  {
184  datadir = std::getenv("G4LEDATA");
185  if(!datadir)
186  {
187  G4Exception("G4LivermoreRayleighModelModel::ReadData()","em0006",
189  "Environment variable G4LEDATA not defined");
190  return;
191  }
192  }
193 
194  //
195 
196  dataCS[Z] = new G4LPhysicsFreeVector();
197 
198  // Activation of spline interpolation
199  //dataCS[Z] ->SetSpline(true);
200 
201  std::ostringstream ostCS;
202  ostCS << datadir << "/livermore/rayl/re-cs-" << Z <<".dat";
203  std::ifstream finCS(ostCS.str().c_str());
204 
205  if( !finCS .is_open() )
206  {
208  ed << "G4LivermorePolarizedRayleighModel data file <" << ostCS.str().c_str()
209  << "> is not opened!" << G4endl;
210  G4Exception("G4LivermorePolarizedRayleighModel::ReadData()","em0003",FatalException,
211  ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
212  return;
213  }
214  else
215  {
216  if(verboseLevel > 3) {
217  G4cout << "File " << ostCS.str()
218  << " is opened by G4LivermoreRayleighModel" << G4endl;
219  }
220  dataCS[Z]->Retrieve(finCS, true);
221  }
222  }
223 
224  //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
225 
227  const G4ParticleDefinition*,
228  G4double GammaEnergy,
231  {
232  if (verboseLevel > 1)
233  {
234  G4cout << "G4LivermoreRayleighModel::ComputeCrossSectionPerAtom()"
235  << G4endl;
236  }
237 
238  if(GammaEnergy < lowEnergyLimit) { return 0.0; }
239 
240  G4double xs = 0.0;
241 
242  G4int intZ = G4lrint(Z);
243 
244  if(intZ < 1 || intZ > maxZ) { return xs; }
245 
246  G4LPhysicsFreeVector* pv = dataCS[intZ];
247 
248  // if element was not initialised
249  // do initialisation safely for MT mode
250  if(!pv) {
251  InitialiseForElement(0, intZ);
252  pv = dataCS[intZ];
253  if(!pv) { return xs; }
254  }
255 
256  G4int n = pv->GetVectorLength() - 1;
257  G4double e = GammaEnergy/MeV;
258  if(e >= pv->Energy(n)) {
259  xs = (*pv)[n]/(e*e);
260  } else if(e >= pv->Energy(0)) {
261  xs = pv->Value(e)/(e*e);
262  }
263 
264  /* if(verboseLevel > 0)
265  {
266  G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
267  << e << G4endl;
268  G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
269  G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0]
270  << G4endl;
271  G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n]
272  << G4endl;
273  G4cout << "*********************************************************"
274  << G4endl;
275  }
276  */
277 
278  return xs;
279  }
280 
281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
282 
283 void G4LivermorePolarizedRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
284  const G4MaterialCutsCouple* couple,
285  const G4DynamicParticle* aDynamicGamma,
286  G4double,
287  G4double)
288 {
289  if (verboseLevel > 3)
290  G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedRayleighModel" << G4endl;
291 
292  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
293 
294  if (photonEnergy0 <= lowEnergyLimit)
295  {
299  return ;
300  }
301 
302  G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection();
303 
304  // Select randomly one element in the current material
305  // G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0);
306  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
307  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
308  G4int Z = (G4int)elm->GetZ();
309 
310  G4double outcomingPhotonCosTheta = GenerateCosTheta(photonEnergy0, Z);
311  G4double outcomingPhotonPhi = GeneratePhi(outcomingPhotonCosTheta);
313 
314  // incomingPhoton reference frame:
315  // z = versor parallel to the incomingPhotonDirection
316  // x = versor parallel to the incomingPhotonPolarization
317  // y = defined as z^x
318 
319  // outgoingPhoton reference frame:
320  // z' = versor parallel to the outgoingPhotonDirection
321  // x' = defined as x-x*z'z' normalized
322  // y' = defined as z'^x'
323 
324  G4ThreeVector z(aDynamicGamma->GetMomentumDirection().unit());
325  G4ThreeVector x(GetPhotonPolarization(*aDynamicGamma));
326  G4ThreeVector y(z.cross(x));
327 
328  // z' = std::cos(phi)*std::sin(theta) x + std::sin(phi)*std::sin(theta) y + std::cos(theta) z
329  G4double xDir;
330  G4double yDir;
331  G4double zDir;
332  zDir=outcomingPhotonCosTheta;
333  xDir=std::sqrt(1-outcomingPhotonCosTheta*outcomingPhotonCosTheta);
334  yDir=xDir;
335  xDir*=std::cos(outcomingPhotonPhi);
336  yDir*=std::sin(outcomingPhotonPhi);
337 
338  G4ThreeVector zPrime((xDir*x + yDir*y + zDir*z).unit());
339  G4ThreeVector xPrime(x.perpPart(zPrime).unit());
340  G4ThreeVector yPrime(zPrime.cross(xPrime));
341 
342  // outgoingPhotonPolarization is directed as x' std::cos(beta) + y' std::sin(beta)
343  G4ThreeVector outcomingPhotonPolarization(xPrime*std::cos(beta) + yPrime*std::sin(beta));
344 
346  fParticleChange->ProposePolarization(outcomingPhotonPolarization);
347  fParticleChange->SetProposedKineticEnergy(photonEnergy0);
348 
349 }
350 
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
352 
354 {
355  // d sigma k0
356  // --------- = r0^2 * pi * F^2(x, Z) * ( 2 - sin^2 theta) * std::sin (theta), x = ---- std::sin(theta/2)
357  // d theta hc
358 
359  // d sigma k0 1 - y
360  // --------- = r0^2 * pi * F^2(x, Z) * ( 1 + y^2), x = ---- std::sqrt ( ------- ), y = std::cos(theta)
361  // d y hc 2
362 
363  // Z
364  // F(x, Z) ~ --------
365  // a + bx
366  //
367  // The time to exit from the outer loop grows as ~ k0
368  // On pcgeant2 the time is ~ 1 s for k0 ~ 1 MeV on the oxygen element. A 100 GeV
369  // event will take ~ 10 hours.
370  //
371  // On the avarage the inner loop does 1.5 iterations before exiting
372 
373  const G4double xFactor = (incomingPhotonEnergy*cm)/(h_Planck*c_light);
374  //const G4VEMDataSet * formFactorData = GetScatterFunctionData();
375 
376  G4double cosTheta;
377  G4double fCosTheta;
378  G4double x;
379  G4double fValue;
380 
381  if (incomingPhotonEnergy > 5.*MeV)
382  {
383  cosTheta = 1.;
384  }
385  else
386  {
387  do
388  {
389  do
390  {
391  cosTheta = 2.*G4UniformRand()-1.;
392  fCosTheta = (1.+cosTheta*cosTheta)/2.;
393  }
394  while (fCosTheta < G4UniformRand());
395 
396  x = xFactor*std::sqrt((1.-cosTheta)/2.);
397 
398  if (x > 1.e+005)
399  fValue = formFactorData->FindValue(x, zAtom-1);
400  else
401  fValue = formFactorData->FindValue(0., zAtom-1);
402 
403  fValue/=zAtom;
404  fValue*=fValue;
405  }
406  while(fValue < G4UniformRand());
407  }
408 
409  return cosTheta;
410 }
411 
412 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
413 
415 {
416  // d sigma
417  // --------- = alpha * ( 1 - sin^2 (theta) * cos^2 (phi) )
418  // d phi
419 
420  // On the average the loop takes no more than 2 iterations before exiting
421 
422  G4double phi;
423  G4double cosPhi;
424  G4double phiProbability;
425  G4double sin2Theta;
426 
427  sin2Theta=1.-cosTheta*cosTheta;
428 
429  do
430  {
431  phi = twopi * G4UniformRand();
432  cosPhi = std::cos(phi);
433  phiProbability= 1. - sin2Theta*cosPhi*cosPhi;
434  }
435  while (phiProbability < G4UniformRand());
436 
437  return phi;
438 }
439 
440 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
441 
443 {
444  // Rayleigh polarization is always on the x' direction
445 
446  return 0;
447 }
448 
449 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
450 
452 {
453 
454 // SI - From G4VLowEnergyDiscretePhotonProcess.cc
455 
456  G4ThreeVector photonMomentumDirection;
457  G4ThreeVector photonPolarization;
458 
459  photonPolarization = photon.GetPolarization();
460  photonMomentumDirection = photon.GetMomentumDirection();
461 
462  if ((!photonPolarization.isOrthogonal(photonMomentumDirection, 1e-6)) || photonPolarization.mag()==0.)
463  {
464  // if |photonPolarization|==0. or |photonPolarization * photonDirection0| > 1e-6 * |photonPolarization ^ photonDirection0|
465  // then polarization is choosen randomly.
466 
467  G4ThreeVector e1(photonMomentumDirection.orthogonal().unit());
468  G4ThreeVector e2(photonMomentumDirection.cross(e1).unit());
469 
471 
472  e1*=std::cos(angle);
473  e2*=std::sin(angle);
474 
475  photonPolarization=e1+e2;
476  }
477  else if (photonPolarization.howOrthogonal(photonMomentumDirection) != 0.)
478  {
479  // if |photonPolarization * photonDirection0| != 0.
480  // then polarization is made orthonormal;
481 
482  photonPolarization=photonPolarization.perpPart(photonMomentumDirection);
483  }
484 
485  return photonPolarization.unit();
486 }
487 
488 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
489 
490  #include "G4AutoLock.hh"
491  namespace { G4Mutex LivermorePolarizedRayleighModelMutex = G4MUTEX_INITIALIZER; }
492 
494  G4int Z)
495  {
496  G4AutoLock l(&LivermorePolarizedRayleighModelMutex);
497  // G4cout << "G4LivermoreRayleighModel::InitialiseForElement Z= "
498  // << Z << G4endl;
499  if(!dataCS[Z]) { ReadData(Z); }
500  l.unlock();
501  }