Geant4  10.02.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  if(IsMaster()) {
69  for(G4int i=0; i<maxZ; ++i) {
70  if(dataCS[i]) {
71  delete dataCS[i];
72  dataCS[i] = 0;
73  }
74  }
75  }
76 }
77 
78 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79 
81  const G4DataVector& cuts)
82 {
83  if (verboseLevel > 1)
84  {
85  G4cout << "Calling Initialise() of G4LivermoreRayleighModel." << G4endl
86  << "Energy range: "
87  << LowEnergyLimit() / eV << " eV - "
88  << HighEnergyLimit() / GeV << " GeV"
89  << G4endl;
90  }
91 
92  if(IsMaster()) {
93 
94  // Initialise element selector
95  InitialiseElementSelectors(particle, cuts);
96 
97  // Access to elements
98  char* path = getenv("G4LEDATA");
99  G4ProductionCutsTable* theCoupleTable =
101  G4int numOfCouples = theCoupleTable->GetTableSize();
102 
103  for(G4int i=0; i<numOfCouples; ++i)
104  {
105  const G4MaterialCutsCouple* couple =
106  theCoupleTable->GetMaterialCutsCouple(i);
107  const G4Material* material = couple->GetMaterial();
108  const G4ElementVector* theElementVector = material->GetElementVector();
109  G4int nelm = material->GetNumberOfElements();
110 
111  for (G4int j=0; j<nelm; ++j)
112  {
113  G4int Z = G4lrint((*theElementVector)[j]->GetZ());
114  if(Z < 1) { Z = 1; }
115  else if(Z > maxZ) { Z = maxZ; }
116  if( (!dataCS[Z]) ) { ReadData(Z, path); }
117  }
118  }
119  }
120 
121  if(isInitialised) { return; }
123  isInitialised = true;
124 
125 }
126 
127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
128 
130  G4VEmModel* masterModel)
131 {
133 }
134 
135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
136 
137 void G4LivermoreRayleighModel::ReadData(size_t Z, const char* path)
138 {
139  if (verboseLevel > 1)
140  {
141  G4cout << "Calling ReadData() of G4LivermoreRayleighModel"
142  << G4endl;
143  }
144 
145  if(dataCS[Z]) { return; }
146 
147  const char* datadir = path;
148 
149  if(!datadir)
150  {
151  datadir = getenv("G4LEDATA");
152  if(!datadir)
153  {
154  G4Exception("G4LivermoreRayleighModelModel::ReadData()","em0006",
156  "Environment variable G4LEDATA not defined");
157  return;
158  }
159  }
160 
161  //
162 
163  dataCS[Z] = new G4LPhysicsFreeVector();
164 
165  // Activation of spline interpolation
166  //dataCS[Z] ->SetSpline(true);
167 
168  std::ostringstream ostCS;
169  ostCS << datadir << "/livermore/rayl/re-cs-" << Z <<".dat";
170  std::ifstream finCS(ostCS.str().c_str());
171 
172  if( !finCS .is_open() )
173  {
175  ed << "G4LivermoreRayleighModel data file <" << ostCS.str().c_str()
176  << "> is not opened!" << G4endl;
177  G4Exception("G4LivermoreRayleighModel::ReadData()","em0003",FatalException,
178  ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
179  return;
180  }
181  else
182  {
183  if(verboseLevel > 3) {
184  G4cout << "File " << ostCS.str()
185  << " is opened by G4LivermoreRayleighModel" << G4endl;
186  }
187  dataCS[Z]->Retrieve(finCS, true);
188  }
189 }
190 
191 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
192 
194  const G4ParticleDefinition*,
195  G4double GammaEnergy,
196  G4double Z, G4double,
198 {
199  if (verboseLevel > 1)
200  {
201  G4cout << "G4LivermoreRayleighModel::ComputeCrossSectionPerAtom()"
202  << G4endl;
203  }
204 
205  if(GammaEnergy < lowEnergyLimit) { return 0.0; }
206 
207  G4double xs = 0.0;
208 
209  G4int intZ = G4lrint(Z);
210 
211  if(intZ < 1 || intZ > maxZ) { return xs; }
212 
213  G4LPhysicsFreeVector* pv = dataCS[intZ];
214 
215  // if element was not initialised
216  // do initialisation safely for MT mode
217  if(!pv) {
218  InitialiseForElement(0, intZ);
219  pv = dataCS[intZ];
220  if(!pv) { return xs; }
221  }
222 
223  G4int n = pv->GetVectorLength() - 1;
224  G4double e = GammaEnergy/MeV;
225  if(e >= pv->Energy(n)) {
226  xs = (*pv)[n]/(e*e);
227  } else if(e >= pv->Energy(0)) {
228  xs = pv->Value(e)/(e*e);
229  }
230 
231  if(verboseLevel > 0)
232  {
233  G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)="
234  << e << G4endl;
235  G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
236  G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0]
237  << G4endl;
238  G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n]
239  << G4endl;
240  G4cout << "*********************************************************"
241  << G4endl;
242  }
243  return xs;
244 }
245 
246 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
247 
249  std::vector<G4DynamicParticle*>*,
250  const G4MaterialCutsCouple* couple,
251  const G4DynamicParticle* aDynamicGamma,
253 {
254  if (verboseLevel > 1) {
255  G4cout << "Calling SampleSecondaries() of G4LivermoreRayleighModel"
256  << G4endl;
257  }
258  G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
259 
260  // absorption of low-energy gamma
261  if (photonEnergy0 <= lowEnergyLimit)
262  {
266  return ;
267  }
268 
269  // Select randomly one element in the current material
270  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
271  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
272  G4int Z = G4lrint(elm->GetZ());
273 
274  // Sample the angle of the scattered photon
275 
276  G4ThreeVector photonDirection =
277  GetAngularDistribution()->SampleDirection(aDynamicGamma,
278  photonEnergy0,
279  Z, couple->GetMaterial());
280  fParticleChange->ProposeMomentumDirection(photonDirection);
281 }
282 
283 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
284 
285 #include "G4AutoLock.hh"
286 namespace { G4Mutex LivermoreRayleighModelMutex = G4MUTEX_INITIALIZER; }
287 
288 void
290  G4int Z)
291 {
292  G4AutoLock l(&LivermoreRayleighModelMutex);
293  // G4cout << "G4LivermoreRayleighModel::InitialiseForElement Z= "
294  // << Z << G4endl;
295  if(!dataCS[Z]) { ReadData(Z); }
296  l.unlock();
297 }
298 
299 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
G4double LowEnergyLimit() const
Definition: G4VEmModel.hh:641
static const double MeV
Definition: G4SIunits.hh:211
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:148
G4double HighEnergyLimit() const
Definition: G4VEmModel.hh:634
G4double GetZ() const
Definition: G4Element.hh:131
virtual void InitialiseForElement(const G4ParticleDefinition *, G4int Z)
G4VEmAngularDistribution * GetAngularDistribution()
Definition: G4VEmModel.hh:617
G4bool IsMaster() const
Definition: G4VEmModel.hh:718
G4ParticleDefinition * GetDefinition() const
size_t GetVectorLength() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
#define G4MUTEX_INITIALIZER
Definition: G4Threading.hh:175
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:214
const G4int n
G4double Energy(size_t index) const
std::vector< G4EmElementSelector * > * GetElementSelectors()
Definition: G4VEmModel.hh:802
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:212
G4int G4Mutex
Definition: G4Threading.hh:173
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:810
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
int G4lrint(double ad)
Definition: templates.hh:163
void SetAngularDistribution(G4VEmAngularDistribution *)
Definition: G4VEmModel.hh:624
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:186
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:544
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
const G4Material * GetMaterial() const