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