Geant4 10.7.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 (std::string)
 
void readDXS ()
 
void BuildCDXS ()
 
void BuildCDXS (G4double, G4double)
 
void NormalizeCDXS ()
 
void InterpolateCDXS ()
 
void PrintDXS (int)
 
G4double SampleAngle (G4double)
 
G4double SampleAngleMT (G4double, G4double)
 
G4double SampleAngleEthylene (G4double, G4double)
 
G4bool IsFileFound () const
 

Detailed Description

Definition at line 32 of file G4LEPTSDiffXS.hh.

Constructor & Destructor Documentation

◆ G4LEPTSDiffXS()

G4LEPTSDiffXS::G4LEPTSDiffXS ( std::string  file)

Definition at line 41 of file G4LEPTSDiffXS.cc.

41 {
42 fileName = file;
43
44 readDXS();
45 BuildCDXS();
46 //BuildCDXS(1.0, 0.5);
47 NormalizeCDXS();
48 InterpolateCDXS();
49}
G4LEPTSDiffXS::InterpolateCDXS
void InterpolateCDXS()
Definition: G4LEPTSDiffXS.cc:172
G4LEPTSDiffXS::NormalizeCDXS
void NormalizeCDXS()
Definition: G4LEPTSDiffXS.cc:155

Member Function Documentation

◆ BuildCDXS() [1/2]

void G4LEPTSDiffXS::BuildCDXS ( )

Definition at line 147 of file G4LEPTSDiffXS.cc.

147 {
148
149 BuildCDXS(1.0, 0.0); // El = 0
150}

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

◆ BuildCDXS() [2/2]

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

Definition at line 124 of file G4LEPTSDiffXS.cc.

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

◆ InterpolateCDXS()

void G4LEPTSDiffXS::InterpolateCDXS ( )

Definition at line 172 of file G4LEPTSDiffXS.cc.

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

Referenced by G4LEPTSDiffXS(), and SampleAngleEthylene().

◆ IsFileFound()

G4bool G4LEPTSDiffXS::IsFileFound ( ) const
inline

Definition at line 48 of file G4LEPTSDiffXS.hh.

48 {
49 return bFileFound;
50 }

◆ NormalizeCDXS()

void G4LEPTSDiffXS::NormalizeCDXS ( )

Definition at line 155 of file G4LEPTSDiffXS.cc.

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

Referenced by G4LEPTSDiffXS(), and SampleAngleEthylene().

◆ PrintDXS()

void G4LEPTSDiffXS::PrintDXS ( int  NE)

Definition at line 330 of file G4LEPTSDiffXS.cc.

330 {
331// Debug
332//#include <string>
333//using namespace std;
334
335
336 G4double dxs;
337
338 G4cout << G4endl<< "DXS & CDXS: " << fileName << G4endl<< G4endl;
339
340 for (G4int aBin=0; aBin<NumAng; aBin++) {
341 if( aBin>0)
342 dxs = (CDXS[NE][aBin] - CDXS[NE][aBin-1])/(CDXS[0][aBin] - CDXS[0][aBin-1]);
343 else
344 dxs = CDXS[NE][aBin];
345
346 G4cout << CDXS[0][aBin] << " " << dxs << " " << CDXS[NE][aBin] << G4endl;
347 }
348
349 G4cout << G4endl<< "IDXS & ICDXS: " << fileName << G4endl<< G4endl;
350
351 for (G4int aBin=0; aBin<INumAng; aBin++) {
352 if( aBin>0)
353 dxs = (ICDXS[NE][aBin] - ICDXS[NE][aBin-1])/(ICDXS[0][aBin] - ICDXS[0][aBin-1]);
354 else
355 dxs = ICDXS[NE][aBin];
356
357 G4cout << ICDXS[0][aBin] << " " << dxs << " " << ICDXS[NE][aBin] << G4endl;
358 }
359
360
361 // if(jmGlobals->VerboseHeaders) {
362 // G4cout << G4endl << "dxskt1" << G4endl;
363 // for (G4int aBin=0;aBin<NumAng;aBin++){
364 // G4cout << DXS[0][aBin] << "\t" << DXS[1][aBin] << "\t" << DXS[2][aBin] << "\t"
365 // << CDXS[1][aBin] << "\t" << KT[12][aBin] << G4endl;
366 // }
367 // }
368
369}
NE
@ NE
Definition: Evaluator.cc:65
G4endl
#define G4endl
Definition: G4ios.hh:57
G4cout
G4GLOB_DLL std::ostream G4cout

◆ readDXS()

void G4LEPTSDiffXS::readDXS ( )

Definition at line 54 of file G4LEPTSDiffXS.cc.

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

Referenced by G4LEPTSDiffXS().

◆ SampleAngle()

G4double G4LEPTSDiffXS::SampleAngle ( G4double  Energy)

Definition at line 219 of file G4LEPTSDiffXS.cc.

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

Referenced by SampleAngleEthylene().

◆ SampleAngleEthylene()

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

Definition at line 251 of file G4LEPTSDiffXS.cc.

251 {
252
253 BuildCDXS(E, El);
254 NormalizeCDXS();
255 InterpolateCDXS();
256
257 return( SampleAngle(E) );
258}
G4LEPTSDiffXS::SampleAngle
G4double SampleAngle(G4double)
Definition: G4LEPTSDiffXS.cc:219

◆ SampleAngleMT()

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

Definition at line 263 of file G4LEPTSDiffXS.cc.

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