59 if ((fp=fopen(fileName.c_str(),
"r"))==NULL){
71 (
void) fscanf(fp,
"%d %d %s", &NumAng, &NumEn, DXSTypeName);
72 if( !strcmp(DXSTypeName,
"KTC") ) DXSType = 2;
73 else if( !strcmp(DXSTypeName,
"KT") ) DXSType = 1;
79 for (
G4int eBin=1; eBin<=NumEn; eBin++){
80 (
void) fscanf(fp,
"%f ",&data);
89 for (
G4int aBin=0;aBin<NumAng;aBin++){
90 (
void) fscanf(fp,
"%f ",&data);
92 for (
G4int eBin=1;eBin<=NumEn;eBin++){
93 (
void) fscanf(fp,
"%f %f ",&data2, &data);
100 for(
G4int aBin=0; aBin<NumAng; aBin++){
101 for(
G4int eBin=0; eBin<=NumEn; eBin++){
102 (
void) fscanf(fp,
"%f ",&data);
106 for(
G4int aBin=0; aBin<NumAng; aBin++){
107 for(
G4int eBin=1; eBin<=NumEn; eBin++){
110 G4double p = sqrt(pow( (E/27.2/137),2) +2*E/27.2);
111 KT[eBin][aBin] = p *sqrt(2.-2.*cos(A*
CLHEP::twopi/360.));
126 for(
G4int aBin=0;aBin<NumAng;aBin++) {
127 for(
G4int eBin=0;eBin<=NumEn;eBin++){
128 CDXS[eBin][aBin]=0.0;
132 for(
G4int aBin=0;aBin<NumAng;aBin++)
133 CDXS[0][aBin] = DXS[0][aBin];
135 for (
G4int eBin=1;eBin<=NumEn;eBin++){
137 for (
G4int aBin=0;aBin<NumAng;aBin++){
138 sum += pow(DXS[eBin][aBin], (1.0-El/E) );
139 CDXS[eBin][aBin]=
sum;
158 for (
G4int eBin=1; eBin<=NumEn; eBin++){
159 G4double area = CDXS[eBin][NumAng-1];
162 for (
G4int aBin=0; aBin<NumAng; aBin++) {
163 CDXS[eBin][aBin] /= area;
177 for(
G4int aBin=0; aBin<NumAng-1; aBin++) {
185 for(
G4int eBin=0; eBin<=NumEn; eBin++) {
197 ICDXS[eBin][ia] = (y1*(x2-
x) + y2*(x-x1))/(x2-x1);
200 ICDXS[eBin][ia] =
G4Exp( (log(y1)*log(x2/x)+log(y2)*log(x/x1))/log(x2/x1) );
203 IKT[eBin][ia] = (z1*(x2-
x) + z2*(x-x1))/(x2-x1);
220 G4int ii,jj,kk=0, Ebin;
223 for(ii=2; ii<=NumEn; ii++)
226 if(Energy > Eb[NumEn]) Ebin=NumEn;
227 else if(Energy > (Eb[Ebin]+Eb[Ebin+1])*0.5 ) Ebin++;
238 if (dxs < rnd) ii=kk;
257 return( SampleAngle(E) );
264 G4int ii, jj, kk=0, Ebin, iMin, iMax;
268 G4double Pi = sqrt( pow( (Ei/27.2/137),2) +2*Ei/27.2);
269 G4double Pd = sqrt( pow( (Ed/27.2/137),2) +2*Ed/27.2);
273 if(Pd <= 1
e-9 )
return (0.0);
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++;
288 if( IKT[Ebin][kk] < Kmin ) ii=kk;
296 if( IKT[Ebin][kk] < Kmax ) ii=kk;
304 rnd = (1-rnd)*ICDXS[Ebin][iMin] + rnd*ICDXS[Ebin][iMax];
311 if( ICDXS[Ebin][kk] < rnd) ii=kk;
318 G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2*Pi*Pd);
340 for (
G4int aBin=0; aBin<NumAng; aBin++) {
342 dxs = (CDXS[
NE][aBin] - CDXS[
NE][aBin-1])/(CDXS[0][aBin] - CDXS[0][aBin-1]);
344 dxs = CDXS[
NE][aBin];
346 G4cout << CDXS[0][aBin] <<
" " << dxs <<
" " << CDXS[
NE][aBin] <<
G4endl;
349 G4cout << G4endl<<
"IDXS & ICDXS: " << fileName << G4endl<<
G4endl;
351 for (
G4int aBin=0; aBin<INumAng; aBin++) {
353 dxs = (ICDXS[
NE][aBin] - ICDXS[
NE][aBin-1])/(ICDXS[0][aBin] - ICDXS[0][aBin-1]);
355 dxs = ICDXS[
NE][aBin];
357 G4cout << ICDXS[0][aBin] <<
" " << dxs <<
" " << ICDXS[
NE][aBin] <<
G4endl;
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Float_t y1[n_points_granero]
G4LEPTSDiffXS(std::string)
Float_t x1[n_points_granero]
Float_t y2[n_points_geant4]
const XML_Char const XML_Char * data
G4double SampleAngleEthylene(G4double, G4double)
double A(double temperature)
static constexpr double twopi
G4GLOB_DLL std::ostream G4cout
G4double SampleAngleMT(G4double, G4double)
G4double SampleAngle(G4double)
static const G4double eps
Float_t x2[n_points_geant4]
typedef void(XMLCALL *XML_ElementDeclHandler)(void *userData