Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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 ( 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();
49 }
void InterpolateCDXS()
void NormalizeCDXS()

Member Function Documentation

void G4LEPTSDiffXS::BuildCDXS ( )

Definition at line 147 of file G4LEPTSDiffXS.cc.

147  {
148 
149  BuildCDXS(1.0, 0.0); // El = 0
150 }
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 }
int G4int
Definition: G4Types.hh:78
double G4double
Definition: G4Types.hh:76
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 }
static const G4double eps
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4bool G4LEPTSDiffXS::IsFileFound ( ) const
inline

Definition at line 48 of file G4LEPTSDiffXS.hh.

48  {
49  return bFileFound;
50  }
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 }
int G4int
Definition: G4Types.hh:78
double G4double
Definition: G4Types.hh:76
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 }
Definition: Evaluator.cc:66
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
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  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  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  fscanf(fp,"%f ",&data);
91  DXS[0][aBin]=(G4double)data;
92  for (G4int eBin=1;eBin<=NumEn;eBin++){
93  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  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 }
const char * p
Definition: xmltok.h:285
int G4int
Definition: G4Types.hh:78
const XML_Char const XML_Char * data
Definition: expat.h:268
G4GLOB_DLL std::ostream G4cout
double A(double temperature)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76
static constexpr double twopi
Definition: SystemOfUnits.h:55

Here is the call graph for this function:

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 }
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76
static constexpr double twopi
Definition: SystemOfUnits.h:55
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 }
G4double SampleAngle(G4double)
void InterpolateCDXS()
void NormalizeCDXS()
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 }
tuple x
Definition: test.py:50
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76

The documentation for this class was generated from the following files: