| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 #include "G4VLEPTSModel.hh"
27
28 #include "CLHEP/Units/SystemOfUnits.h"
29
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 G4VLEPTSModel::G4VLEPTSModel(const G4String& m
32 {
33 theMeanFreePathTable=nullptr;
34
35 theNumbBinTable=100;
36
37 verboseLevel = 0;
38
39 theLowestEnergyLimit = 0.5*CLHEP::eV;
40
41 theHighestEnergyLimit = 1.0*CLHEP::MeV;
42
43 theXSType = XSEnergy;
44 }
45
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48 G4VLEPTSModel::~G4VLEPTSModel()
49 {
50
51 if(theMeanFreePathTable != nullptr) {
52 theMeanFreePathTable->clearAndDestroy();
53 delete theMeanFreePathTable;
54 }
55 }
56
57
58 //....oooOO0OOooo........oooOO0OOooo........oo
59 void G4VLEPTSModel::Init()
60 {
61 theLowestEnergyLimit = 0.5*CLHEP::eV;
62 theHighestEnergyLimit = 1.0*CLHEP::MeV;
63 //t theHighestEnergyLimit = 15.0*CLHEP::M
64 SetLowEnergyLimit(theLowestEnergyLimit);
65 SetHighEnergyLimit(theHighestEnergyLimit);
66 theNumbBinTable = 100;
67
68 }
69
70
71
72 //....oooOO0OOooo........oooOO0OOooo........oo
73 G4double G4VLEPTSModel::GetMeanFreePath(const
74 const G4ParticleDefinition* ,
75 G4double kineticEnergy )
76 {
77 G4double MeanFreePath;
78
79 if( verboseLevel >= 3 ) G4cout << aMaterial-
80 if (kineticEnergy > theHighestEnergyLimit ||
81 MeanFreePath = DBL_MAX;
82 else
83 MeanFreePath = (*theMeanFreePathTable)(aMa
84
85 return MeanFreePath;
86 }
87
88
89 //....oooOO0OOooo........oooOO0OOooo........oo
90 void G4VLEPTSModel::BuildPhysicsTable(const G4
91 {
92 //CHECK IF PATH VARIABLE IS DEFINED
93 const char* path = G4FindDataDir("G4LEDATA")
94 if( path == nullptr ) {
95 G4Exception("G4VLEPTSModel",
96 "",
97 FatalException,
98 "variable G4LEDATA not defined");
99 }
100
101 // Build microscopic cross section table and
102
103 G4String aParticleName = aParticleType.GetPa
104
105 if (theMeanFreePathTable != nullptr) {
106 theMeanFreePathTable->clearAndDestroy();
107 delete theMeanFreePathTable;
108 }
109
110 theMeanFreePathTable = new G4PhysicsTable(G4
111
112 //LOOP TO MATERIALS IN GEOMETRY
113 const G4MaterialTable * materialTable = G4Ma
114 std::vector<G4Material*>::const_iterator mat
115 for( matite = materialTable->cbegin(); matit
116 const G4Material * aMaterial = (*matite);
117 G4String mateName = aMaterial->GetName();
118
119 //READ PARAMETERS FOR THIS MATERIAL
120 std::string dirName = std::string(path) +
121 std::string fnParam = dirName + mateName
122 std::string baseName = std::string(path) +
123 G4bool bData = ReadParam( fnParam, aMateri
124 if( !bData ) continue; // MATERIAL NOT EX
125
126 //READ INTEGRAL CROSS SECTION FOR THIS MAT
127 std::string fnIXS = baseName + ".IXS.dat"
128
129 std::map< G4int, std::vector<G4double> > i
130 if( verboseLevel >= 2 ) G4cout << GetName(
131
132 if( integralXS.empty() ) {
133 G4cerr << " Integral cross sections will
134 auto ptrVector = new G4PhysicsLogVector
135 ptrVector->PutValue(0, DBL_MAX);
136 ptrVector->PutValue(1, DBL_MAX);
137
138 std::size_t matIdx = aMaterial->GetIndex
139 theMeanFreePathTable->insertAt( matIdx ,
140
141 } else {
142
143 if( verboseLevel >= 2 ) {
144 std::map< G4int, std::vector<G4double> >::co
145 for( itei = integralXS.begin(); itei != inte
146 G4cout << GetName() << " : " << (*itei).fi
147 }
148 }
149
150 BuildMeanFreePathTable( aMaterial, integ
151
152 std::string fnDXS = baseName + ".DXS.dat
153 std::string fnRMT = baseName + ".RMT.dat
154 std::string fnEloss = baseName + ".Eloss
155 std::string fnEloss2 = baseName + ".Elos
156
157 theDiffXS[aMaterial] = new G4LEPTSDiffXS
158 if( !theDiffXS[aMaterial]->IsFileFound()
159 G4Exception("G4VLEPTSModel::BuildPhysicsTabl
160 "",
161 FatalException,
162 (G4String("File not found :" + fnDXS).
163 }
164
165 theRMTDistr[aMaterial] = new G4LEPTSDist
166 theRMTDistr[aMaterial]->ReadFile(fnRMT);
167
168 theElostDistr[aMaterial] = new G4LEPTSEl
169 if( !theElostDistr[aMaterial]->IsFileFou
170 G4Exception("G4VLEPTSModel::BuildPhysicsTabl
171 "",
172 FatalException,
173 (G4String("File not found :" + fnEloss
174 }
175 }
176
177 }
178
179 }
180
181 void G4VLEPTSModel::BuildMeanFreePathTable( co
182 {
183 G4double LowEdgeEnergy, fValue;
184
185 //BUILD MEAN FREE PATH TABLE FROM INTEGRAL C
186 std::size_t matIdx = aMaterial->GetIndex();
187 auto ptrVector = new G4PhysicsLogVector(the
188
189 for (G4int ii=0; ii < theNumbBinTable; ++ii)
190 LowEdgeEnergy = ptrVector->Energy(ii);
191 if( verboseLevel >= 2 ) G4cout << GetName(
192 //- fValue = ComputeMFP(LowEdgeEnergy
193 fValue = 0.;
194 if( LowEdgeEnergy >= theLowestEnergyLimit
195 LowEdgeEnergy <= theHighestEnergyLimit) {
196 G4double NbOfMoleculesPerVolume = aMater
197
198 G4double SIGMA = 0. ;
199 //- for ( std::size_t elm=0 ; elm <
200 G4double crossSection = 0.;
201
202 G4double eVEnergy = LowEdgeEnergy/CLHEP::eV;
203
204 //- if( verboseLevel >= 2 ) G4cout << "
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];
212 if( verboseLevel >= 3 ) G4cout << " GET
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]
222
223 if( verboseLevel >= 3 ) G4cout << " crossS
224
225 // SIGMA += NbOfAtomsPerVolume[elm] * c
226 SIGMA = NbOfMoleculesPerVolume * crossSect
227 if( verboseLevel >= 2 ) G4cout << GetName(
228 << " Bin " << Bin << " TOTAL " << a
229 << " XS " << integralXS[theXSType][
230 }
231
232 //-}
233
234 fValue = SIGMA > DBL_MIN ? 1./SIGMA : DB
235 }
236
237 ptrVector->PutValue(ii, fValue);
238 if( verboseLevel >= 2 ) G4cout << GetName(
239 }
240
241 theMeanFreePathTable->insertAt( matIdx , ptr
242 }
243
244
245 //....oooOO0OOooo........oooOO0OOooo........oo
246 G4double G4VLEPTSModel::SampleAngle(const G4Ma
247 {
248 G4double x;
249
250 if( e < 10001) {
251 x = theDiffXS[aMaterial]->SampleAngleMT(e,
252 }
253 else {
254 G4double Ei = e;
255 G4double Ed = e -el;
256
257 G4double Pi = std::sqrt( std::pow( (Ei/27.
258 G4double Pd = std::sqrt( std::pow( (Ed/27.
259
260 G4double Kmin = Pi - Pd;
261 G4double Kmax = Pi + Pd;
262
263 G4double KR = theRMTDistr[aMaterial]->Samp
264
265 G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2
266 if( co > 1. ) co = 1.;
267 x = std::acos(co); //*360/twopi;
268 }
269 return(x);
270 }
271
272 //....oooOO0OOooo........oooOO0OOooo........oo
273 G4ThreeVector G4VLEPTSModel::SampleNewDirectio
274
275 G4double x = SampleAngle(aMaterial, e, el);
276
277 G4double cosTeta = std::cos(x); //*twopi/360
278 G4double sinTeta = std::sqrt(1.0-cosTeta*cos
279 G4double Phi = CLHEP::twopi * G4UniformR
280 G4double dirx = sinTeta*std::cos(Phi) , d
281
282 G4ThreeVector P1Dir(dirx, diry, dirz);
283 #ifdef DEBUG_LEPTS
284 if( verboseLevel >= 2 ) G4cout << " G4VLEPTS
285 #endif
286 P1Dir.rotateUz(P0Dir);
287 #ifdef DEBUG_LEPTS
288 if( verboseLevel >= 2 ) G4cout << " G4VLEPTS
289 #endif
290
291 return(P1Dir);
292 }
293
294
295 //....oooOO0OOooo........oooOO0OOooo........oo
296 G4ThreeVector G4VLEPTSModel::SampleNewDirectio
297 {
298 G4double cosTeta = std::cos(x); //*twopi/360
299 G4double sinTeta = std::sqrt(1.0-cosTeta*cos
300 G4double Phi = CLHEP::twopi * G4UniformR
301 G4double dirx = sinTeta*std::cos(Phi) , d
302
303 G4ThreeVector P1Dir( dirx,diry,dirz );
304 P1Dir.rotateUz(P0Dir);
305
306 return(P1Dir);
307 }
308
309
310 //....oooOO0OOooo........oooOO0OOooo........oo
311 G4double G4VLEPTSModel::SampleEnergyLoss(const
312 {
313 G4double el;
314 el = theElostDistr[aMaterial]->Sample(eMin/C
315
316 #ifdef DEBUG_LEPTS
317 if( verboseLevel >= 2 ) G4cout << aMaterial-
318 << " " << GetName() << G4endl;
319 #endif
320 return el;
321 }
322
323
324 //....oooOO0OOooo........oooOO0OOooo........oo
325 G4bool G4VLEPTSModel::ReadParam(const G4String
326 {
327 std::ifstream fin(fnParam);
328 if (!fin.is_open()) {
329 G4Exception("G4VLEPTSModel::ReadParam",
330 "",
331 JustWarning,
332 (G4String("File not found: ")+ fnParam).c_
333 return false;
334 }
335
336 G4double IonisPot, IonisPotInt;
337
338 fin >> IonisPot >> IonisPotInt;
339 if( verboseLevel >= 1 ) G4cout << "Read para
340 << " IonisPotInt: " << IonisPotInt << G4en
341
342 theIonisPot[aMaterial] = IonisPot * CLHEP::e
343 theIonisPotInt[aMaterial] = IonisPotInt * CL
344
345 G4double MolecularMass = 0;
346 auto nelem = (G4int)aMaterial->GetNumberOfE
347 const G4int* atomsV = aMaterial->GetAtomsVe
348 for( G4int ii = 0; ii < nelem; ++ii ) {
349 MolecularMass += aMaterial->GetElement(ii)
350 // G4cout << " MMASS1 " << mmass/CLHEP:
351 }
352 // G4cout << " MMASS " << MolecularMass <<
353 theMolecularMass[aMaterial] = MolecularMass*
354 // theMolecularMass[aMaterial] = aMaterial-
355
356 if( verboseLevel >= 1) G4cout << " IonisPot:
357 << " IonisPotInt: " << IonisPotInt/CLH
358 << " MolecularMass " << MolecularMass/
359
360 return true;
361 }
362
363 //....oooOO0OOooo........oooOO0OOooo........oo
364 std::map< G4int, std::vector<G4double> > G4VLE
365 {
366 std::map< G4int, std::vector<G4double> > int
367 //G4cout << "fnIXS (" << fnIXS << ")" << G4e
368
369 std::ifstream fin(fnIXS);
370 if (!fin.is_open()) {
371 G4Exception("G4VLEPTSModel::ReadIXS",
372 "",
373 JustWarning,
374 (G4String("File not found: ")+ fnIXS).c_st
375 return integralXS;
376 }
377
378 G4int nXSdat, nXSsub;
379 fin >> nXSdat >> nXSsub;
380 if( verboseLevel >= 1 ) G4cout << "Read IXS
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 << GetNa
394 // xsdat 1e-16*cm2
395 }
396 }
397 fin.close();
398
399 return integralXS;
400 }
401
402