ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LivermoreRayleighModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LivermoreRayleighModel.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 // 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 = std::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 = (*theElementVector)[j]->GetZasInt();
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 = std::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,
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  /*
262  if (photonEnergy0 <= lowEnergyLimit)
263  {
264  fParticleChange->ProposeTrackStatus(fStopAndKill);
265  fParticleChange->SetProposedKineticEnergy(0.);
266  fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
267  return ;
268  }
269  */
270  // Select randomly one element in the current material
271  const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
272  const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
273  G4int Z = G4lrint(elm->GetZ());
274 
275  // Sample the angle of the scattered photon
276 
277  G4ThreeVector photonDirection =
278  GetAngularDistribution()->SampleDirection(aDynamicGamma,
279  photonEnergy0,
280  Z, couple->GetMaterial());
281  fParticleChange->ProposeMomentumDirection(photonDirection);
282 }
283 
284 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
285 
286 #include "G4AutoLock.hh"
287 namespace { G4Mutex LivermoreRayleighModelMutex = G4MUTEX_INITIALIZER; }
288 
289 void
291  G4int Z)
292 {
293  G4AutoLock l(&LivermoreRayleighModelMutex);
294  // G4cout << "G4LivermoreRayleighModel::InitialiseForElement Z= "
295  // << Z << G4endl;
296  if(!dataCS[Z]) { ReadData(Z); }
297  l.unlock();
298 }
299 
300 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....