Geant4 11.3.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
Public Member Functions | List of all members
G4LEPTSDiffXS Class Reference

#include <G4LEPTSDiffXS.hh>

Public Member Functions

 G4LEPTSDiffXS (const G4String &)
 
void readDXS ()
 
void BuildCDXS ()
 
void BuildCDXS (G4double, G4double)
 
void NormalizeCDXS ()
 
void InterpolateCDXS ()
 
void PrintDXS (G4int)
 
G4double SampleAngle (G4double)
 
G4double SampleAngleMT (G4double, G4double)
 
G4double SampleAngleEthylene (G4double, G4double)
 
G4bool IsFileFound () const
 

Detailed Description

Definition at line 31 of file G4LEPTSDiffXS.hh.

Constructor & Destructor Documentation

◆ G4LEPTSDiffXS()

G4LEPTSDiffXS::G4LEPTSDiffXS ( const G4String & file)

Definition at line 37 of file G4LEPTSDiffXS.cc.

37 {
38 fileName = file;
39
40 readDXS();
41 BuildCDXS();
42 //BuildCDXS(1.0, 0.5);
43 NormalizeCDXS();
44 InterpolateCDXS();
45}

Member Function Documentation

◆ BuildCDXS() [1/2]

void G4LEPTSDiffXS::BuildCDXS ( )

Definition at line 143 of file G4LEPTSDiffXS.cc.

143 {
144
145 BuildCDXS(1.0, 0.0); // El = 0
146}

Referenced by BuildCDXS(), G4LEPTSDiffXS(), and SampleAngleEthylene().

◆ BuildCDXS() [2/2]

void G4LEPTSDiffXS::BuildCDXS ( G4double E,
G4double El )

Definition at line 120 of file G4LEPTSDiffXS.cc.

120 {
121
122 for(G4int aBin=0;aBin<NumAng;aBin++) {
123 for(G4int eBin=0;eBin<=NumEn;eBin++){
124 CDXS[eBin][aBin]=0.0;
125 }
126 }
127
128 for(G4int aBin=0;aBin<NumAng;aBin++)
129 CDXS[0][aBin] = DXS[0][aBin];
130
131 for (G4int eBin=1;eBin<=NumEn;eBin++){
132 G4double sum=0.0;
133 for (G4int aBin=0;aBin<NumAng;aBin++){
134 sum += std::pow(DXS[eBin][aBin], (1.0-El/E) );
135 CDXS[eBin][aBin]=sum;
136 }
137 }
138}
G4double
double G4double
Definition G4Types.hh:83
G4int
int G4int
Definition G4Types.hh:85

◆ InterpolateCDXS()

void G4LEPTSDiffXS::InterpolateCDXS ( )

Definition at line 168 of file G4LEPTSDiffXS.cc.

168 { // *10 angles, linear
169
170 G4double eps = 1e-5;
171 G4int ia = 0;
172
173 for( G4int aBin=0; aBin<NumAng-1; aBin++) {
174 G4double x1 = CDXS[0][aBin] + eps;
175 G4double x2 = CDXS[0][aBin+1] + eps;
176 G4double dx = (x2-x1)/100;
177
178 //if( x1<10 || x1) dx = (x2-x1)/100;
179
180 for( G4double x=x1; x < (x2-dx/10); x += dx) {
181 for( G4int eBin=0; eBin<=NumEn; eBin++) {
182 G4double y1 = CDXS[eBin][aBin];
183 G4double y2 = CDXS[eBin][aBin+1];
184 G4double z1 = KT[eBin][aBin];
185 G4double z2 = KT[eBin][aBin+1];
186
187 if( aBin==0 ) {
188 y1 /=100;
189 z1 /=100;
190 }
191
192 if( eBin==0 ) { //linear abscisa
193 ICDXS[eBin][ia] = (y1*(x2-x) + y2*(x-x1))/(x2-x1);
194 }
195 else { //log-log ordenada
196 ICDXS[eBin][ia] = G4Exp( (std::log(y1)*std::log(x2/x)+std::log(y2)*std::log(x/x1))/std::log(x2/x1) );
197 }
198
199 IKT[eBin][ia] = (z1*(x2-x) + z2*(x-x1))/(x2-x1);
200 //IKT[eBin][ia] = exp( (log(z1)*log(x2/x)+log(z2)*log(x/x1))/log(x2/x1) );
201 }
202
203 ia++;
204 }
205
206 }
207
208 INumAng = ia;
209}
G4Exp
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition G4Exp.hh:180

Referenced by G4LEPTSDiffXS(), and SampleAngleEthylene().

◆ IsFileFound()

G4bool G4LEPTSDiffXS::IsFileFound ( ) const
inline

Definition at line 47 of file G4LEPTSDiffXS.hh.

47{ return bFileFound; }

◆ NormalizeCDXS()

void G4LEPTSDiffXS::NormalizeCDXS ( )

Definition at line 151 of file G4LEPTSDiffXS.cc.

151 {
152
153 // Normalize: 1/area
154 for (G4int eBin=1; eBin<=NumEn; eBin++){
155 G4double area = CDXS[eBin][NumAng-1];
156 //G4cout << eBin << " area = " << area << G4endl;
157
158 for (G4int aBin=0; aBin<NumAng; aBin++) {
159 CDXS[eBin][aBin] /= area;
160 //DXS[eBin][aBin] /= area;
161 }
162 }
163}

Referenced by G4LEPTSDiffXS(), and SampleAngleEthylene().

◆ PrintDXS()

void G4LEPTSDiffXS::PrintDXS ( G4int NE)

Definition at line 325 of file G4LEPTSDiffXS.cc.

325 {
326
327 G4double dxs;
328
329 G4cout << G4endl<< "DXS & CDXS: " << fileName << G4endl<< G4endl;
330
331 for (G4int aBin=0; aBin<NumAng; aBin++) {
332 if( aBin>0)
333 dxs = (CDXS[NE][aBin] - CDXS[NE][aBin-1])/(CDXS[0][aBin] - CDXS[0][aBin-1]);
334 else
335 dxs = CDXS[NE][aBin];
336
337 G4cout << CDXS[0][aBin] << " " << dxs << " " << CDXS[NE][aBin] << G4endl;
338 }
339
340 G4cout << G4endl<< "IDXS & ICDXS: " << fileName << G4endl<< G4endl;
341
342 for (G4int aBin=0; aBin<INumAng; aBin++) {
343 if( aBin>0)
344 dxs = (ICDXS[NE][aBin] - ICDXS[NE][aBin-1])/(ICDXS[0][aBin] - ICDXS[0][aBin-1]);
345 else
346 dxs = ICDXS[NE][aBin];
347
348 G4cout << ICDXS[0][aBin] << " " << dxs << " " << ICDXS[NE][aBin] << G4endl;
349 }
350
351
352 // if(jmGlobals->VerboseHeaders) {
353 // G4cout << G4endl << "dxskt1" << G4endl;
354 // for (G4int aBin=0;aBin<NumAng;aBin++){
355 // G4cout << DXS[0][aBin] << "\t" << DXS[1][aBin] << "\t" << DXS[2][aBin] << "\t"
356 // << CDXS[1][aBin] << "\t" << KT[12][aBin] << G4endl;
357 // }
358 // }
359
360}
NE
@ NE
Definition Evaluator.cc:68
G4endl
#define G4endl
Definition G4ios.hh:67
G4cout
G4GLOB_DLL std::ostream G4cout

◆ readDXS()

void G4LEPTSDiffXS::readDXS ( )

Definition at line 50 of file G4LEPTSDiffXS.cc.

50 {
51
52 FILE *fp;
53 G4float data, data2;
54
55 if ((fp=fopen(fileName.c_str(), "r"))==nullptr){
56 //G4cout << "Error reading " << fileName << G4endl;
57 NumEn = 0;
58 bFileFound = false;
59 return;
60 }
61
62 bFileFound = true;
63
64 //G4cout << "Reading2 " << fileName << G4endl;
65
66 //NumAng = 181;
67 (void) fscanf(fp, "%d %d %s", &NumAng, &NumEn, DXSTypeName);
68 if( strcmp(DXSTypeName, "KTC") == 0 ) DXSType = 2; // read DXS & calculate KT
69 else if( strcmp(DXSTypeName, "KT") == 0 ) DXSType = 1; // read DXS & KT
70 else DXSType = 0;
71
72 // if( verboseLevel >= 1 ) G4cout << "Read DXS (" << fileName << ")\t NEg " << NumEn << " NAng " << NumAng
73 // << "DXSType " << DXSTypeName << " " << DXSType << G4endl;
74
75 for (G4int eBin=1; eBin<=NumEn; eBin++){
76 (void) fscanf(fp,"%f ",&data);
77 Eb[eBin] = (G4double)data;
78 }
79
80
81 //for (aBin=1;aBin<NumAng;aBin++){
82
83 if(DXSType==1) {
84 G4cout << "DXSTYpe 1" << G4endl;
85 for (G4int aBin=0;aBin<NumAng;aBin++){
86 (void) fscanf(fp,"%f ",&data);
87 DXS[0][aBin]=(G4double)data;
88 for (G4int eBin=1;eBin<=NumEn;eBin++){
89 (void) fscanf(fp,"%f %f ",&data2, &data);
90 DXS[eBin][aBin]=(G4double)data;
91 KT[eBin][aBin]=(G4double)data2;
92 }
93 }
94 }
95 else {
96 for(G4int aBin=0; aBin<NumAng; aBin++){
97 for(G4int eBin=0; eBin<=NumEn; eBin++){
98 (void) fscanf(fp,"%f ",&data);
99 DXS[eBin][aBin] = (G4double)data;
100 }
101 }
102 for(G4int aBin=0; aBin<NumAng; aBin++){
103 for(G4int eBin=1; eBin<=NumEn; eBin++){
104 G4double A = DXS[0][aBin]; // Angle
105 G4double E = Eb[eBin]; // Energy
106 G4double p = std::sqrt(std::pow( (E/27.2/137),2) +2*E/27.2); // Momentum
107 KT[eBin][aBin] = p *std::sqrt(2.-2.*std::cos(A*CLHEP::twopi/360.)); // Mom. Transfer
108 //G4cout << "aEpKt " << aBin << " " << A << " E " << E << " p " << p << " KT "
109 // << KT[eBin][aBin] << " DXS " << DXS[eBin][aBin] << G4endl;
110 }
111 }
112 }
113
114 fclose(fp);
115}
G4float
float G4float
Definition G4Types.hh:84
A
const G4double A[17]
Definition G4WaterStopping.cc:53

Referenced by G4LEPTSDiffXS().

◆ SampleAngle()

G4double G4LEPTSDiffXS::SampleAngle ( G4double Energy)

Definition at line 214 of file G4LEPTSDiffXS.cc.

214 {
215 G4int ii,jj,kk=0, Ebin;
216
217 Ebin=1;
218 for(ii=2; ii<=NumEn; ii++)
219 if(Energy >= Eb[ii])
220 Ebin=ii;
221 if(Energy > Eb[NumEn]) Ebin=NumEn;
222 else if(Energy > (Eb[Ebin]+Eb[Ebin+1])*0.5 ) Ebin++;
223
224 //G4cout << "SampleAngle E " << Energy << " Ebin " << Ebin << " E[] " << Eb[Ebin] << G4endl;
225
226 ii=0;
227 jj=INumAng-1;
228 G4double rnd=G4UniformRand();
229
230 while ((jj-ii)>1){
231 kk=(ii+jj)/2;
232 G4double dxs = ICDXS[Ebin][kk];
233 if (dxs < rnd) ii=kk;
234 else jj=kk;
235 }
236
237
238 //G4double x = ICDXS[0][jj];
239 G4double x = ICDXS[0][kk] *CLHEP::twopi/360.;
240
241 return(x);
242}
G4UniformRand
#define G4UniformRand()
Definition Randomize.hh:52

Referenced by SampleAngleEthylene().

◆ SampleAngleEthylene()

G4double G4LEPTSDiffXS::SampleAngleEthylene ( G4double E,
G4double El )

Definition at line 246 of file G4LEPTSDiffXS.cc.

246 {
247
248 BuildCDXS(E, El);
249 NormalizeCDXS();
250 InterpolateCDXS();
251
252 return( SampleAngle(E) );
253}
G4LEPTSDiffXS::SampleAngle
G4double SampleAngle(G4double)
Definition G4LEPTSDiffXS.cc:214

◆ SampleAngleMT()

G4double G4LEPTSDiffXS::SampleAngleMT ( G4double Energy,
G4double Elost )

Definition at line 258 of file G4LEPTSDiffXS.cc.

258 {
259 G4int ii, jj, kk=0, Ebin, iMin, iMax;
260
261 G4double Ei = Energy;
262 G4double Ed = Energy - Elost;
263 G4double Pi = std::sqrt( std::pow( (Ei/27.2/137),2) +2*Ei/27.2); //incidente
264 G4double Pd = std::sqrt( std::pow( (Ed/27.2/137),2) +2*Ed/27.2); //dispersado
265 G4double Kmin = Pi - Pd;
266 G4double Kmax = Pi + Pd;
267
268 if(Pd <= 1e-9 ) return (0.0);
269
270
271 // locate Energy bin
272 Ebin=1;
273 for(ii=2; ii<=NumEn; ii++)
274 if(Energy > Eb[ii]) Ebin=ii;
275 if(Energy > Eb[NumEn]) Ebin=NumEn;
276 else if(Energy > (Eb[Ebin]+Eb[Ebin+1])*0.5 ) Ebin++;
277
278 //G4cout << "SampleAngle2 E " << Energy << " Ebin " << Ebin << " E[] " << Eb[Ebin] << G4endl;
279
280 ii=0; jj=INumAng-1;
281 while ((jj-ii)>1) {
282 kk=(ii+jj)/2;
283 if( IKT[Ebin][kk] < Kmin ) ii=kk;
284 else jj=kk;
285 }
286 iMin = ii;
287
288 ii=0; jj=INumAng-1;
289 while ((jj-ii)>1) {
290 kk=(ii+jj)/2;
291 if( IKT[Ebin][kk] < Kmax ) ii=kk;
292 else jj=kk;
293 }
294 iMax = ii;
295
296
297 // r -> a + (b-a)*r = a*(1-r) + b*r
298 G4double rnd = G4UniformRand();
299 rnd = (1-rnd)*ICDXS[Ebin][iMin] + rnd*ICDXS[Ebin][iMax];
300 //G4double rnd = (ICDXS[Ebin][iMax] - ICDXS[Ebin][iMin]) * G4UniformRand()
301 // + ICDXS[Ebin][iMin];
302
303 ii=0; jj=INumAng-1;
304 while ((jj-ii)>1){
305 kk=(ii+jj)/2;
306 if( ICDXS[Ebin][kk] < rnd) ii=kk;
307 else jj=kk;
308 }
309
310 //Sampled
311 G4double KR = IKT[Ebin][kk];
312
313 G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2*Pi*Pd); //cos ang. disp.
314 if(co > 1) co =1;
315 G4double x = std::acos(co); //*360/twopi; //ang. dispers.
316
317 // Elastic aprox.
318 //x = 2*asin(KR/Pi/2)*360/twopi;
319
320 return(x);
321}
The documentation for this class was generated from the following files: