Geant4 11.2.2
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4VLEPTSModel.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#include "G4VLEPTSModel.hh"
27
29
30//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31G4VLEPTSModel::G4VLEPTSModel(const G4String& modelName) : G4VEmModel(modelName)
32{
34
36
37 verboseLevel = 0;
38
39 theLowestEnergyLimit = 0.5*CLHEP::eV;
40
41 theHighestEnergyLimit = 1.0*CLHEP::MeV;
42
44}
45
46
47//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56
57
58//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60{
61 theLowestEnergyLimit = 0.5*CLHEP::eV;
62 theHighestEnergyLimit = 1.0*CLHEP::MeV;
63 //t theHighestEnergyLimit = 15.0*CLHEP::MeV;
66 theNumbBinTable = 100;
67
68}
69
70
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
75 G4double kineticEnergy )
76{
77 G4double MeanFreePath;
78
79 if( verboseLevel >= 3 ) G4cout << aMaterial->GetIndex() << " G4VLEPTSModel::GetMeanFreePath " << kineticEnergy << " > " << theHighestEnergyLimit << " < " << theLowestEnergyLimit << G4endl;
80 if (kineticEnergy > theHighestEnergyLimit || kineticEnergy < theLowestEnergyLimit)
81 MeanFreePath = DBL_MAX;
82 else
83 MeanFreePath = (*theMeanFreePathTable)(aMaterial->GetIndex())->Value(kineticEnergy);
84
85 return MeanFreePath;
86}
87
88
89//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
91{
92 //CHECK IF PATH VARIABLE IS DEFINED
93 const char* path = G4FindDataDir("G4LEDATA");
94 if( path == nullptr ) {
95 G4Exception("G4VLEPTSModel",
96 "",
98 "variable G4LEDATA not defined");
99 }
100
101 // Build microscopic cross section table and mean free path table
102
103 G4String aParticleName = aParticleType.GetParticleName();
104
105 if (theMeanFreePathTable != nullptr) {
108 }
109
111
112 //LOOP TO MATERIALS IN GEOMETRY
113 const G4MaterialTable * materialTable = G4Material::GetMaterialTable() ;
114 std::vector<G4Material*>::const_iterator matite;
115 for( matite = materialTable->begin(); matite != materialTable->end(); matite++ ) {
116 const G4Material * aMaterial = (*matite);
117 G4String mateName = aMaterial->GetName();
118
119 //READ PARAMETERS FOR THIS MATERIAL
120 std::string dirName = std::string(path) + "/lepts/";
121 std::string fnParam = dirName + mateName + "." + aParticleName + ".param.dat";
122 std::string baseName = std::string(path) + "/lepts/" + mateName + "." + aParticleName;
123 G4bool bData = ReadParam( fnParam, aMaterial );
124 if( !bData ) continue; // MATERIAL NOT EXISTING, DO NOT READ OTHER FILES
125
126 //READ INTEGRAL CROSS SECTION FOR THIS MATERIAL
127 std::string fnIXS = baseName + ".IXS.dat";
128
129 std::map< G4int, std::vector<G4double> > integralXS = ReadIXS(fnIXS, aMaterial);
130 if( verboseLevel >= 2 ) G4cout << GetName() << " : " << theXSType << " " << mateName << " INTEGRALXS " << integralXS.size() << G4endl;
131
132 if( integralXS.empty() ) {
133 G4cerr << " Integral cross sections will be set to 0. for material " << mateName << G4endl;
135 ptrVector->PutValue(0, DBL_MAX);
136 ptrVector->PutValue(1, DBL_MAX);
137
138 std::size_t matIdx = aMaterial->GetIndex();
139 theMeanFreePathTable->insertAt( matIdx , ptrVector ) ;
140
141 } else {
142
143 if( verboseLevel >= 2 ) {
144 std::map< G4int, std::vector<G4double> >::const_iterator itei;
145 for( itei = integralXS.begin(); itei != integralXS.end(); itei++ ){
146 G4cout << GetName() << " : " << (*itei).first << " INTEGRALXS NDATA " << (*itei).second.size() << G4endl;
147 }
148 }
149
150 BuildMeanFreePathTable( aMaterial, integralXS );
151
152 std::string fnDXS = baseName + ".DXS.dat";
153 std::string fnRMT = baseName + ".RMT.dat";
154 std::string fnEloss = baseName + ".Eloss.dat";
155 std::string fnEloss2 = baseName + ".Eloss2.dat";
156
157 theDiffXS[aMaterial] = new G4LEPTSDiffXS(fnDXS);
158 if( !theDiffXS[aMaterial]->IsFileFound() ) {
159 G4Exception("G4VLEPTSModel::BuildPhysicsTable",
160 "",
162 (G4String("File not found :" + fnDXS).c_str()));
163 }
164
165 theRMTDistr[aMaterial] = new G4LEPTSDistribution();
166 theRMTDistr[aMaterial]->ReadFile(fnRMT);
167
168 theElostDistr[aMaterial] = new G4LEPTSElossDistr(fnEloss);
169 if( !theElostDistr[aMaterial]->IsFileFound() ) {
170 G4Exception("G4VLEPTSModel::BuildPhysicsTable",
171 "",
173 (G4String("File not found :" + fnEloss).c_str()));
174 }
175 }
176
177 }
178
179}
180
181void G4VLEPTSModel::BuildMeanFreePathTable( const G4Material* aMaterial, std::map< G4int, std::vector<G4double> >& integralXS )
182{
183 G4double LowEdgeEnergy, fValue;
184
185 //BUILD MEAN FREE PATH TABLE FROM INTEGRAL CROSS SECTION
186 std::size_t matIdx = aMaterial->GetIndex();
188
189 for (G4int ii=0; ii < theNumbBinTable; ++ii) {
190 LowEdgeEnergy = ptrVector->Energy(ii);
191 if( verboseLevel >= 2 ) G4cout << GetName() << " " << ii << " Energy " << LowEdgeEnergy << " > " << theLowestEnergyLimit << " < " << theHighestEnergyLimit << G4endl;
192 //- fValue = ComputeMFP(LowEdgeEnergy, material, aParticleName);
193 fValue = 0.;
194 if( LowEdgeEnergy >= theLowestEnergyLimit &&
195 LowEdgeEnergy <= theHighestEnergyLimit) {
196 G4double NbOfMoleculesPerVolume = aMaterial->GetDensity()/theMolecularMass[aMaterial]*CLHEP::Avogadro;
197
198 G4double SIGMA = 0. ;
199 //- for ( std::size_t elm=0 ; elm < aMaterial->GetNumberOfElements() ; elm++ ) {
200 G4double crossSection = 0.;
201
202 G4double eVEnergy = LowEdgeEnergy/CLHEP::eV;
203
204 //- if( verboseLevel >= 2 ) G4cout << " eVEnergy " << eVEnergy << " LowEdgeE " << LowEdgeEnergy << " " << integralXS[theXSType][1] << G4endl;
205
206 if(eVEnergy < integralXS[0][1] ) {
207 crossSection = 0.;
208 } else {
209 G4int Bin = 0; // locate bin
210 G4double aa, bb;
211 for( G4int jj=1; jj<theNXSdat[aMaterial]; jj++) { // Extrapolate for E > Emax !!!
212 if( verboseLevel >= 3 ) G4cout << " GET BIN " << jj << " "<< eVEnergy << " > " << integralXS[0][jj] << G4endl;
213 if( eVEnergy > integralXS[0][jj]) {
214 Bin = jj;
215 } else {
216 break;
217 }
218 }
219 aa = integralXS[0][Bin];
220 bb = integralXS[0][Bin+1];
221 crossSection = (integralXS[theXSType][Bin] + (integralXS[theXSType][Bin+1]-integralXS[theXSType][Bin])/(bb-aa)*(eVEnergy-aa) ) * 1.e-16*CLHEP::cm2;
222
223 if( verboseLevel >= 3 ) G4cout << " crossSection " << crossSection << " " <<integralXS[theXSType][Bin] << " + " << (integralXS[theXSType][Bin+1]-integralXS[theXSType][Bin]) << " / " << (bb-aa) << " *" << (eVEnergy-aa) << " * " << 1.e-16*CLHEP::cm2 << G4endl;;
224
225 // SIGMA += NbOfAtomsPerVolume[elm] * crossSection;
226 SIGMA = NbOfMoleculesPerVolume * crossSection;
227 if( verboseLevel >= 2 ) G4cout << GetName() << " ADDING SIGMA " << SIGMA << " NAtoms " << NbOfMoleculesPerVolume
228 << " Bin " << Bin << " TOTAL " << aa << " " << bb
229 << " XS " << integralXS[theXSType][Bin] << " " << integralXS[theXSType][Bin+1] << G4endl;
230 }
231
232 //-}
233
234 fValue = SIGMA > DBL_MIN ? 1./SIGMA : DBL_MAX;
235 }
236
237 ptrVector->PutValue(ii, fValue);
238 if( verboseLevel >= 2 ) G4cout << GetName() << " BUILDXS " << ii << " : " << LowEdgeEnergy << " = " << fValue << G4endl;
239 }
240
241 theMeanFreePathTable->insertAt( matIdx , ptrVector ) ;
242}
243
244
245//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
247{
248 G4double x;
249
250 if( e < 10001) {
251 x = theDiffXS[aMaterial]->SampleAngleMT(e, el);
252 }
253 else {
254 G4double Ei = e; //incidente
255 G4double Ed = e -el; //dispersado
256
257 G4double Pi = std::sqrt( std::pow( (Ei/27.2/137),2) +2*Ei/27.2); //incidente
258 G4double Pd = std::sqrt( std::pow( (Ed/27.2/137),2) +2*Ed/27.2); //dispersado
259
260 G4double Kmin = Pi - Pd;
261 G4double Kmax = Pi + Pd;
262
263 G4double KR = theRMTDistr[aMaterial]->Sample(Kmin, Kmax); //sorteo mom. transf.
264
265 G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2*Pi*Pd); //cos ang. disp.
266 if( co > 1. ) co = 1.;
267 x = std::acos(co); //*360/twopi; //ang. dispers.
268 }
269 return(x);
270}
271
272//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
274
275 G4double x = SampleAngle(aMaterial, e, el);
276
277 G4double cosTeta = std::cos(x); //*twopi/360.0);
278 G4double sinTeta = std::sqrt(1.0-cosTeta*cosTeta);
279 G4double Phi = CLHEP::twopi * G4UniformRand() ;
280 G4double dirx = sinTeta*std::cos(Phi) , diry = sinTeta*std::sin(Phi) , dirz = cosTeta ;
281
282 G4ThreeVector P1Dir(dirx, diry, dirz);
283#ifdef DEBUG_LEPTS
284 if( verboseLevel >= 2 ) G4cout << " G4VLEPTSModel::SampleNewDirection " <<P1Dir << G4endl;
285#endif
286 P1Dir.rotateUz(P0Dir);
287#ifdef DEBUG_LEPTS
288 if( verboseLevel >= 2 ) G4cout << " G4VLEPTSModel::SampleNewDirection rotated " <<P1Dir << " " << P0Dir << G4endl;
289#endif
290
291 return(P1Dir);
292}
293
294
295//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
297{
298 G4double cosTeta = std::cos(x); //*twopi/360.0);
299 G4double sinTeta = std::sqrt(1.0-cosTeta*cosTeta);
300 G4double Phi = CLHEP::twopi * G4UniformRand() ;
301 G4double dirx = sinTeta*std::cos(Phi) , diry = sinTeta*std::sin(Phi) , dirz = cosTeta ;
302
303 G4ThreeVector P1Dir( dirx,diry,dirz );
304 P1Dir.rotateUz(P0Dir);
305
306 return(P1Dir);
307}
308
309
310//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
312{
313 G4double el;
314 el = theElostDistr[aMaterial]->Sample(eMin/CLHEP::eV, eMax/CLHEP::eV)*CLHEP::eV;
315
316#ifdef DEBUG_LEPTS
317 if( verboseLevel >= 2 ) G4cout << aMaterial->GetName() <<"SampleEnergyLoss/eV " << el/CLHEP::eV << " eMax/eV " << eMax/CLHEP::eV << " "
318 << " " << GetName() << G4endl;
319#endif
320 return el;
321}
322
323
324//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
326{
327 std::ifstream fin(fnParam);
328 if (!fin.is_open()) {
329 G4Exception("G4VLEPTSModel::ReadParam",
330 "",
332 (G4String("File not found: ")+ fnParam).c_str());
333 return false;
334 }
335
336 G4double IonisPot, IonisPotInt;
337
338 fin >> IonisPot >> IonisPotInt;
339 if( verboseLevel >= 1 ) G4cout << "Read param (" << fnParam << ")\t IonisPot: " << IonisPot
340 << " IonisPotInt: " << IonisPotInt << G4endl;
341
342 theIonisPot[aMaterial] = IonisPot * CLHEP::eV;
343 theIonisPotInt[aMaterial] = IonisPotInt * CLHEP::eV;
344
345 G4double MolecularMass = 0;
346 auto nelem = (G4int)aMaterial->GetNumberOfElements();
347 const G4int* atomsV = aMaterial->GetAtomsVector();
348 for( G4int ii = 0; ii < nelem; ++ii ) {
349 MolecularMass += aMaterial->GetElement(ii)->GetA()*atomsV[ii]/CLHEP::g;
350 // G4cout << " MMASS1 " << mmass/CLHEP::g << " " << aMaterial->GetElement(ii)->GetName() << " " << aMaterial->GetElement(ii)->GetA()/CLHEP::g << G4endl;
351 }
352 // G4cout << " MMASS " << MolecularMass << " " << MolecularMass*CLHEP::g << " ME " << mmass << " " << mmass/CLHEP::g << G4endl;
353 theMolecularMass[aMaterial] = MolecularMass* CLHEP::g/CLHEP::mole;
354 // theMolecularMass[aMaterial] = aMaterial->GetMassOfMolecule()*CLHEP::Avogadro; // Material mixtures do not calculate molecular mass
355
356 if( verboseLevel >= 1) G4cout << " IonisPot: " << IonisPot/CLHEP::eV << " eV "
357 << " IonisPotInt: " << IonisPotInt/CLHEP::eV << " eV"
358 << " MolecularMass " << MolecularMass/(CLHEP::g/CLHEP::mole) << " g/mole" << G4endl;
359
360 return true;
361}
362
363//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
364std::map< G4int, std::vector<G4double> > G4VLEPTSModel::ReadIXS(G4String fnIXS, const G4Material* aMaterial )
365{
366 std::map< G4int, std::vector<G4double> > integralXS; // process type - energy
367 //G4cout << "fnIXS (" << fnIXS << ")" << G4endl;
368
369 std::ifstream fin(fnIXS);
370 if (!fin.is_open()) {
371 G4Exception("G4VLEPTSModel::ReadIXS",
372 "",
374 (G4String("File not found: ")+ fnIXS).c_str());
375 return integralXS;
376 }
377
378 G4int nXSdat, nXSsub;
379 fin >> nXSdat >> nXSsub;
380 if( verboseLevel >= 1 ) G4cout << "Read IXS (" << fnIXS << ")\t nXSdat: " << nXSdat
381 << " nXSsub: " << nXSsub << G4endl;
382 theNXSdat[aMaterial] = nXSdat;
383 theNXSsub[aMaterial] = nXSsub;
384
385 G4double xsdat;
386 for (G4int ip=0; ip<=nXSsub; ip++) {
387 integralXS[ip].push_back(0.);
388 }
389 for (G4int ie=1; ie<=nXSdat; ie++) {
390 for (G4int ip=0; ip<=nXSsub; ip++) {
391 fin >> xsdat;
392 integralXS[ip].push_back(xsdat);
393 if( verboseLevel >= 3 ) G4cout << GetName() << " FILL IXS " << ip << " " << ie << " = " << integralXS[ip][ie] << " " << xsdat << G4endl;
394 // xsdat 1e-16*cm2
395 }
396 }
397 fin.close();
398
399 return integralXS;
400}
401
const char * G4FindDataDir(const char *)
@ JustWarning
@ FatalException
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
std::vector< G4Material * > G4MaterialTable
double G4double
Definition G4Types.hh:83
bool G4bool
Definition G4Types.hh:86
int G4int
Definition G4Types.hh:85
@ XSEnergy
G4GLOB_DLL std::ostream G4cerr
#define G4endl
Definition G4ios.hh:67
G4GLOB_DLL std::ostream G4cout
#define G4UniformRand()
Definition Randomize.hh:52
Hep3Vector & rotateUz(const Hep3Vector &)
G4double GetA() const
Definition G4Element.hh:127
G4double GetDensity() const
const G4Element * GetElement(G4int iel) const
const G4int * GetAtomsVector() const
static std::size_t GetNumberOfMaterials()
std::size_t GetIndex() const
static G4MaterialTable * GetMaterialTable()
std::size_t GetNumberOfElements() const
const G4String & GetName() const
const G4String & GetParticleName() const
void clearAndDestroy()
void insertAt(std::size_t, G4PhysicsVector *)
void SetHighEnergyLimit(G4double)
void SetLowEnergyLimit(G4double)
const G4String & GetName() const
G4bool ReadParam(G4String fileName, const G4Material *aMaterial)
std::map< const G4Material *, G4LEPTSDiffXS * > theDiffXS
G4ThreeVector SampleNewDirection(const G4Material *aMaterial, G4ThreeVector Dir, G4double e, G4double el)
G4double SampleEnergyLoss(const G4Material *aMaterial, G4double eMin, G4double eMax)
std::map< const G4Material *, G4double > theMolecularMass
~G4VLEPTSModel() override
G4double GetMeanFreePath(const G4Material *mate, const G4ParticleDefinition *aParticle, G4double kineticEnergy)
std::map< const G4Material *, G4int > theNXSsub
G4VLEPTSModel(const G4String &processName)
std::map< const G4Material *, G4LEPTSElossDistr * > theElostDistr
std::map< const G4Material *, G4double > theIonisPotInt
std::map< const G4Material *, G4double > theIonisPot
G4PhysicsTable * theMeanFreePathTable
G4double theHighestEnergyLimit
void BuildPhysicsTable(const G4ParticleDefinition &aParticleType)
virtual std::map< G4int, std::vector< G4double > > ReadIXS(G4String fileName, const G4Material *aMaterial)
G4double theLowestEnergyLimit
std::map< const G4Material *, G4LEPTSDistribution * > theRMTDistr
void BuildMeanFreePathTable(const G4Material *aMaterial, std::map< G4int, std::vector< G4double > > &integralXS)
std::map< const G4Material *, G4int > theNXSdat
G4double SampleAngle(const G4Material *aMaterial, G4double e, G4double el)
#define DBL_MIN
Definition templates.hh:54
#define DBL_MAX
Definition templates.hh:62