Geant4  10.02.p03
Public Member Functions | Private Attributes | List of all members
G4LEPTSDiffXS Class Reference

#include <G4LEPTSDiffXS.hh>

Collaboration diagram for G4LEPTSDiffXS:

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
 

Private Attributes

std::string fileName
 
int NumAng
 
int INumAng
 
int NumEn
 
char DXSTypeName [8]
 
int DXSType
 
G4double Eb [100]
 
G4double DXS [100][190]
 
G4double CDXS [100][190]
 
G4double ICDXS [100][19000]
 
G4double KT [100][190]
 
G4double IKT [100][19000]
 
G4bool bFileFound
 

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 }
file
TFile * file
Definition: comparison_ascii.C:5
G4LEPTSDiffXS::fileName
std::string fileName
Definition: G4LEPTSDiffXS.hh:53
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 }

◆ 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 }
G4int
int G4int
Definition: G4Types.hh:78
G4LEPTSDiffXS::DXS
G4double DXS[100][190]
Definition: G4LEPTSDiffXS.hh:61
G4LEPTSDiffXS::CDXS
G4double CDXS[100][190]
Definition: G4LEPTSDiffXS.hh:61
G4double
double G4double
Definition: G4Types.hh:76

◆ 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] = exp( (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 }
y2
Double_t y2[nxs]
Definition: Hadr00/scripts/Style.C:21
y1
Double_t y1[nxs]
Definition: Hadr00/scripts/Style.C:20
x2
Double_t x2[nxs]
Definition: Hadr00/scripts/Style.C:19
G4LEPTSDiffXS::KT
G4double KT[100][190]
Definition: G4LEPTSDiffXS.hh:63
eps
static const G4double eps
Definition: G4ElectroNuclearCrossSection.cc:99
G4int
int G4int
Definition: G4Types.hh:78
G4LEPTSDiffXS::IKT
G4double IKT[100][19000]
Definition: G4LEPTSDiffXS.hh:63
x1
Double_t x1[nxs]
Definition: Hadr00/scripts/Style.C:18
G4LEPTSDiffXS::CDXS
G4double CDXS[100][190]
Definition: G4LEPTSDiffXS.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4LEPTSDiffXS::ICDXS
G4double ICDXS[100][19000]
Definition: G4LEPTSDiffXS.hh:61

◆ 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 }
G4int
int G4int
Definition: G4Types.hh:78
G4LEPTSDiffXS::CDXS
G4double CDXS[100][190]
Definition: G4LEPTSDiffXS.hh:61
G4double
double G4double
Definition: G4Types.hh:76

◆ 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
Definition: Evaluator.cc:66
G4int
int G4int
Definition: G4Types.hh:78
G4LEPTSDiffXS::fileName
std::string fileName
Definition: G4LEPTSDiffXS.hh:53
G4cout
G4GLOB_DLL std::ostream G4cout
G4LEPTSDiffXS::CDXS
G4double CDXS[100][190]
Definition: G4LEPTSDiffXS.hh:61
G4endl
#define G4endl
Definition: G4ios.hh:61
G4double
double G4double
Definition: G4Types.hh:76
G4LEPTSDiffXS::ICDXS
G4double ICDXS[100][19000]
Definition: G4LEPTSDiffXS.hh:61

◆ 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  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 }
G4LEPTSDiffXS::DXSTypeName
char DXSTypeName[8]
Definition: G4LEPTSDiffXS.hh:57
G4LEPTSDiffXS::KT
G4double KT[100][190]
Definition: G4LEPTSDiffXS.hh:63
G4int
int G4int
Definition: G4Types.hh:78
G4LEPTSDiffXS::fileName
std::string fileName
Definition: G4LEPTSDiffXS.hh:53
G4cout
G4GLOB_DLL std::ostream G4cout
A
double A(double temperature)
Definition: G4DNAElectronHoleRecombination.cc:59
G4LEPTSDiffXS::DXS
G4double DXS[100][190]
Definition: G4LEPTSDiffXS.hh:61
G4endl
#define G4endl
Definition: G4ios.hh:61
fclose
fclose(fp)
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::twopi
static const double twopi
Definition: SystemOfUnits.h:54
G4LEPTSDiffXS::Eb
G4double Eb[100]
Definition: G4LEPTSDiffXS.hh:59
Here is the call graph for this function:

◆ 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 }
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4double
double G4double
Definition: G4Types.hh:76
CLHEP::twopi
static const double twopi
Definition: SystemOfUnits.h:54
G4LEPTSDiffXS::ICDXS
G4double ICDXS[100][19000]
Definition: G4LEPTSDiffXS.hh:61
G4LEPTSDiffXS::Eb
G4double Eb[100]
Definition: G4LEPTSDiffXS.hh:59

◆ 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
G4LEPTSDiffXS::InterpolateCDXS
void InterpolateCDXS()
Definition: G4LEPTSDiffXS.cc:172
G4LEPTSDiffXS::NormalizeCDXS
void NormalizeCDXS()
Definition: G4LEPTSDiffXS.cc:155

◆ 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 }
G4int
int G4int
Definition: G4Types.hh:78
G4UniformRand
#define G4UniformRand()
Definition: Randomize.hh:97
G4LEPTSDiffXS::IKT
G4double IKT[100][19000]
Definition: G4LEPTSDiffXS.hh:63
G4double
double G4double
Definition: G4Types.hh:76
G4LEPTSDiffXS::ICDXS
G4double ICDXS[100][19000]
Definition: G4LEPTSDiffXS.hh:61
G4LEPTSDiffXS::Eb
G4double Eb[100]
Definition: G4LEPTSDiffXS.hh:59

Member Data Documentation

◆ bFileFound

G4bool G4LEPTSDiffXS::bFileFound
private

Definition at line 65 of file G4LEPTSDiffXS.hh.

◆ CDXS

G4double G4LEPTSDiffXS::CDXS[100][190]
private

Definition at line 61 of file G4LEPTSDiffXS.hh.

◆ DXS

G4double G4LEPTSDiffXS::DXS[100][190]
private

Definition at line 61 of file G4LEPTSDiffXS.hh.

◆ DXSType

int G4LEPTSDiffXS::DXSType
private

Definition at line 58 of file G4LEPTSDiffXS.hh.

◆ DXSTypeName

char G4LEPTSDiffXS::DXSTypeName[8]
private

Definition at line 57 of file G4LEPTSDiffXS.hh.

◆ Eb

G4double G4LEPTSDiffXS::Eb[100]
private

Definition at line 59 of file G4LEPTSDiffXS.hh.

◆ fileName

std::string G4LEPTSDiffXS::fileName
private

Definition at line 53 of file G4LEPTSDiffXS.hh.

◆ ICDXS

G4double G4LEPTSDiffXS::ICDXS[100][19000]
private

Definition at line 61 of file G4LEPTSDiffXS.hh.

◆ IKT

G4double G4LEPTSDiffXS::IKT[100][19000]
private

Definition at line 63 of file G4LEPTSDiffXS.hh.

◆ INumAng

int G4LEPTSDiffXS::INumAng
private

Definition at line 55 of file G4LEPTSDiffXS.hh.

◆ KT

G4double G4LEPTSDiffXS::KT[100][190]
private

Definition at line 63 of file G4LEPTSDiffXS.hh.

◆ NumAng

int G4LEPTSDiffXS::NumAng
private

Definition at line 54 of file G4LEPTSDiffXS.hh.

◆ NumEn

int G4LEPTSDiffXS::NumEn
private

Definition at line 56 of file G4LEPTSDiffXS.hh.

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