Geant4  10.00.p02
G4LivermoreRayleighModel.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 // Author: Sebastien Incerti
27 // 31 March 2012
28 // on base of G4LivermoreRayleighModel
29 //
30 
32 #include "G4SystemOfUnits.hh"
34 
35 using namespace std;
36 
37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 
42 
44  :G4VEmModel("LivermoreRayleigh"),isInitialised(false)
45 {
46  fParticleChange = 0;
47  lowEnergyLimit = 10 * eV;
48 
50 
51  verboseLevel= 0;
52  // Verbosity scale for debugging purposes:
53  // 0 = nothing
54  // 1 = calculation of cross sections, file openings...
55  // 2 = entering in methods
56 
57  if(verboseLevel > 0)
58  {
59  G4cout << "G4LivermoreRayleighModel is constructed " << G4endl;
60  }
61 
62 }
63 
64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
65 
67 {}
68 
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70 
72  const G4DataVector& cuts)
73 {
74  if (verboseLevel > 1)
75  {
76  G4cout << "Calling Initialise() of G4LivermoreRayleighModel." << G4endl
77  << "Energy range: "
78  << LowEnergyLimit() / eV << " eV - "
79  << HighEnergyLimit() / GeV << " GeV"
80  << G4endl;
81  }
82 
83  if(IsMaster()) {
84 
85  // Initialise element selector
86  InitialiseElementSelectors(particle, cuts);
87 
88  // Access to elements
89  char* path = getenv("G4LEDATA");
90  G4ProductionCutsTable* theCoupleTable =
92  G4int numOfCouples = theCoupleTable->GetTableSize();
93 
94  for(G4int i=0; i<numOfCouples; ++i)
95  {
96  const G4MaterialCutsCouple* couple =
97  theCoupleTable->GetMaterialCutsCouple(i);
98  const G4Material* material = couple->GetMaterial();
99  const G4ElementVector* theElementVector = material->GetElementVector();
100  G4int nelm = material->GetNumberOfElements();
101 
102  for (G4int j=0; j<nelm; ++j)
103  {
104  G4int Z = G4lrint((*theElementVector)[j]->GetZ());
105  if(Z < 1) { Z = 1; }
106  else if(Z > maxZ) { Z = maxZ; }
107  if( (!dataCS[Z]) ) { ReadData(Z, path); }
108  }
109  }
110  }
111 
112  if(isInitialised) { return; }
114  isInitialised = true;
115 
116 }
117 
118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
119 
121  G4VEmModel* masterModel)
122 {
124 }
125 
126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
127 
128 void G4LivermoreRayleighModel::ReadData(size_t Z, const char* path)
129 {
130  if (verboseLevel > 1)
131  {
132  G4cout << "Calling ReadData() of G4LivermoreRayleighModel"
133  << G4endl;
134  }
135 
136  if(dataCS[Z]) { return; }
137 
138  const char* datadir = path;
139 
140  if(!datadir)
141  {
142  datadir = getenv("G4LEDATA");
143  if(!datadir)
144  {
145  G4Exception("G4LivermoreRayleighModelModel::ReadData()","em0006",
147  "Environment variable G4LEDATA not defined");
148  return;
149  }
150  }
151 
152  //
153 
154  dataCS[Z] = new G4LPhysicsFreeVector();
155 
156  // Activation of spline interpolation
157  //dataCS[Z] ->SetSpline(true);
158 
159  std::ostringstream ostCS;
160  ostCS << datadir << "/livermore/rayl/re-cs-" << Z <<".dat";
161  std::ifstream finCS(ostCS.str().c_str());
162 
163  if( !finCS .is_open() )
164  {
166  ed << "G4LivermoreRayleighModel data file <" << ostCS.str().c_str()
167  << "> is not opened!" << G4endl;
168  G4Exception("G4LivermoreRayleighModel::ReadData()","em0003",FatalException,
169  ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
170  return;
171  }
172  else
173  {
174  if(verboseLevel > 3) {
175  G4cout << "File " << ostCS.str()
176  << " is opened by G4LivermoreRayleighModel" << G4endl;
177  }
178  dataCS[Z]->Retrieve(finCS, true);
179  }
180 }
181 
182 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
183 
185  const G4ParticleDefinition*,
186  G4double GammaEnergy,
187  G4double Z, G4double,
189 {
190  if (verboseLevel > 1)
191  {
192  G4cout << "G4LivermoreRayleighModel::ComputeCrossSectionPerAtom()"
193  << G4endl;
194  }
195 
196  if(GammaEnergy < lowEnergyLimit) { return 0.0; }
197 
198  G4double xs = 0.0;
199 
200  G4int intZ = G4lrint(Z);
201 
202  if(intZ < 1 || intZ > maxZ) { return xs; }
203 
204  G4LPhysicsFreeVector* pv = dataCS[intZ];
205 
206  // if element was not initialised
207  // do initialisation safely for MT mode
208  if(!pv) {
209  InitialiseForElement(0, intZ);
210  pv = dataCS[intZ];
211  if(!pv) { return xs; }
212  }
213 
214  G4int n = pv->GetVectorLength() - 1;
215  G4double e = GammaEnergy/MeV;
216  if(e >= pv->Energy(n)) {
217  xs = (*pv)[n]/(e*e);
218  } else if(e >= pv->Energy(0)) {
219  xs = pv->Value(e)/(e*e);
220  }
221 
222  if(verboseLevel > 0)
223  {
224  G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
225  << e << G4endl;
226  G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
227  G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0]
228  << G4endl;
229  G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n]
230  << G4endl;
231  G4cout << "*********************************************************"
232  << G4endl;
233  }
234  return xs;
235 }
236 
237 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
238 
240  std::vector<G4DynamicParticle*>*,
241  const G4MaterialCutsCouple* couple,
242  const G4DynamicParticle* aDynamicGamma,
244 {
245  if (verboseLevel > 1) {
246  G4cout << "Calling SampleSecondaries() of G4LivermoreRayleighModel"
247  << G4endl;
248  }
249  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
250 
251  // absorption of low-energy gamma
252  if (photonEnergy0 <= lowEnergyLimit)
253  {
257  return ;
258  }
259 
260  // Select randomly one element in the current material
261  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
262  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
263  G4int Z = G4lrint(elm->GetZ());
264 
265  // Sample the angle of the scattered photon
266 
267  G4ThreeVector photonDirection =
268  GetAngularDistribution()->SampleDirection(aDynamicGamma,
269  photonEnergy0,
270  Z, couple->GetMaterial());
271  fParticleChange->ProposeMomentumDirection(photonDirection);
272 }
273 
274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
275 
276 #include "G4AutoLock.hh"
277 namespace { G4Mutex LivermoreRayleighModelMutex = G4MUTEX_INITIALIZER; }
278 
279 void
281  G4int Z)
282 {
283  G4AutoLock l(&LivermoreRayleighModelMutex);
284  // G4cout << "G4LivermoreRayleighModel::InitialiseForElement Z= "
285  // << Z << G4endl;
286  if(!dataCS[Z]) { ReadData(Z); }
287  l.unlock();
288 }
289 
290 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:599
static const double MeV
Definition: G4SIunits.hh:193
std::vector< G4Element * > G4ElementVector
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4double GetKineticEnergy() const
virtual void InitialiseLocal(const G4ParticleDefinition *, G4VEmModel *masterModel)
CLHEP::Hep3Vector G4ThreeVector
void InitialiseElementSelectors(const G4ParticleDefinition *, const G4DataVector &)
Definition: G4VEmModel.cc:135
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:592
G4double GetZ() const
Definition: G4Element.hh:131
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:578
G4bool IsMaster() const
Definition: G4VEmModel.hh:676
G4ParticleDefinition * GetDefinition() const
size_t GetVectorLength() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:158
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4GLOB_DLL std::ostream G4cout
virtual G4ThreeVector & SampleDirection(const G4DynamicParticle *dp, G4double finalTotalEnergy, G4int Z, const G4Material *)=0
void ReadData(size_t Z, const char *path=0)
static const double GeV
Definition: G4SIunits.hh:196
const G4int n
G4double Energy(size_t index) const
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:760
G4double Value(G4double theEnergy, size_t &lastidx) const
static G4LPhysicsFreeVector * dataCS[101]
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
static G4ProductionCutsTable * GetProductionCutsTable()
static const double eV
Definition: G4SIunits.hh:194
G4int G4Mutex
Definition: G4Threading.hh:156
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double kinEnergy, G4double Z, G4double A=0, G4double cut=0, G4double emax=DBL_MAX)
void SetElementSelectors(std::vector< G4EmElementSelector * > *)
Definition: G4VEmModel.hh:768
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
int G4lrint(double ad)
Definition: templates.hh:163
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:585
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii=false)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void SetProposedKineticEnergy(G4double proposedKinEnergy)
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin, G4double maxEnergy)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
G4ParticleChangeForGamma * fParticleChange
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:510
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:121
const G4Material * GetMaterial() const