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.

                                        {
  fileName = file;
 
  readDXS();
  BuildCDXS();
  //BuildCDXS(1.0, 0.5);
  NormalizeCDXS();
  InterpolateCDXS();
}

Member Function Documentation

BuildCDXS() [1/2]

void G4LEPTSDiffXS::BuildCDXS

(

)

Definition at line 147 of file G4LEPTSDiffXS.cc.

                              {
 
  BuildCDXS(1.0, 0.0); // El = 0
}

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

BuildCDXS() [2/2]

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

Definition at line 124 of file G4LEPTSDiffXS.cc.

                                                     {
 
  for(G4int aBin=0;aBin<NumAng;aBin++) {
    for(G4int eBin=0;eBin<=NumEn;eBin++){
      CDXS[eBin][aBin]=0.0;
    }
  }
 
  for(G4int aBin=0;aBin<NumAng;aBin++)
    CDXS[0][aBin] = DXS[0][aBin];
 
  for (G4int eBin=1;eBin<=NumEn;eBin++){
    G4double sum=0.0;
    for (G4int aBin=0;aBin<NumAng;aBin++){
      sum += pow(DXS[eBin][aBin], (1.0-El/E) );
      CDXS[eBin][aBin]=sum;
    }
  }
}

InterpolateCDXS()

void G4LEPTSDiffXS::InterpolateCDXS

(

)

Definition at line 172 of file G4LEPTSDiffXS.cc.

                                    {  // *10 angles, linear
 
  G4double eps = 1e-5;
  G4int ia = 0;
 
  for( G4int aBin=0; aBin<NumAng-1; aBin++) {
    G4double x1 = CDXS[0][aBin] + eps;
    G4double x2 = CDXS[0][aBin+1] + eps;
    G4double dx = (x2-x1)/100;
 
    //if( x1<10 || x1) dx = (x2-x1)/100;
 
    for( G4double x=x1; x < (x2-dx/10); x += dx) {
      for( G4int eBin=0; eBin<=NumEn; eBin++) {
    G4double y1 = CDXS[eBin][aBin];
    G4double y2 = CDXS[eBin][aBin+1];
    G4double z1 = KT[eBin][aBin];
    G4double z2 = KT[eBin][aBin+1];
 
    if( aBin==0 ) {
      y1 /=100;
      z1 /=100;
    }
 
    if( eBin==0 ) {   //linear abscisa
      ICDXS[eBin][ia] = (y1*(x2-x) + y2*(x-x1))/(x2-x1);
    }
    else {           //log-log ordenada
      ICDXS[eBin][ia] = G4Exp( (log(y1)*log(x2/x)+log(y2)*log(x/x1))/log(x2/x1) );
    }
 
    IKT[eBin][ia] = (z1*(x2-x) + z2*(x-x1))/(x2-x1);
    //IKT[eBin][ia] = exp( (log(z1)*log(x2/x)+log(z2)*log(x/x1))/log(x2/x1) );
      }
 
      ia++;
    }
 
  }
 
  INumAng = ia;
}

Referenced by G4LEPTSDiffXS(), and SampleAngleEthylene().

IsFileFound()

G4bool G4LEPTSDiffXS::IsFileFound ( ) const

inline

Definition at line 48 of file G4LEPTSDiffXS.hh.

                             {
    return bFileFound;
  }

NormalizeCDXS()

void G4LEPTSDiffXS::NormalizeCDXS

(

)

Definition at line 155 of file G4LEPTSDiffXS.cc.

                                  {
 
  // Normalize:  1/area
  for (G4int eBin=1; eBin<=NumEn; eBin++){
    G4double area = CDXS[eBin][NumAng-1];
    //G4cout << eBin << " area = " << area << G4endl;
 
    for (G4int aBin=0; aBin<NumAng; aBin++) {
      CDXS[eBin][aBin] /= area;
      //DXS[eBin][aBin]  /= area;
    }
  }
}

Referenced by G4LEPTSDiffXS(), and SampleAngleEthylene().

PrintDXS()

void G4LEPTSDiffXS::PrintDXS

(

int

NE

)

Definition at line 330 of file G4LEPTSDiffXS.cc.

                                     {
// Debug
//#include <string>
//using namespace std;
 
 
  G4double dxs;
 
  G4cout << G4endl<< "DXS & CDXS: " << fileName << G4endl<< G4endl;
 
  for (G4int aBin=0; aBin<NumAng; aBin++) {
    if( aBin>0)
      dxs = (CDXS[NE][aBin] - CDXS[NE][aBin-1])/(CDXS[0][aBin] - CDXS[0][aBin-1]);
    else
      dxs = CDXS[NE][aBin];
 
    G4cout << CDXS[0][aBin] << " " << dxs << " " << CDXS[NE][aBin] << G4endl;
  }
 
  G4cout << G4endl<< "IDXS & ICDXS: " << fileName << G4endl<< G4endl;
 
  for (G4int aBin=0; aBin<INumAng; aBin++) {
    if( aBin>0)
      dxs = (ICDXS[NE][aBin] - ICDXS[NE][aBin-1])/(ICDXS[0][aBin] - ICDXS[0][aBin-1]);
    else
      dxs = ICDXS[NE][aBin];
 
    G4cout << ICDXS[0][aBin] << " " << dxs << " " << ICDXS[NE][aBin] << G4endl;
  }
 
 
  // if(jmGlobals->VerboseHeaders) {
  //   G4cout << G4endl << "dxskt1" << G4endl;
  //   for (G4int aBin=0;aBin<NumAng;aBin++){
  //     G4cout << DXS[0][aBin] << "\t" << DXS[1][aBin] << "\t" << DXS[2][aBin] << "\t"
  //         << CDXS[1][aBin] << "\t" << KT[12][aBin] << G4endl;
  //   }
  // }
 
}

readDXS()

void G4LEPTSDiffXS::readDXS

(

)

Definition at line 54 of file G4LEPTSDiffXS.cc.

                             {
 
  FILE   *fp;
  float data, data2;
 
  if ((fp=fopen(fileName.c_str(), "r"))==NULL){
    //G4cout << "Error reading " << fileName << G4endl;
    NumEn = 0;
    bFileFound = false;
    return;
  }
 
  bFileFound = true;
 
  //G4cout << "Reading2 " << fileName << G4endl;
 
  //NumAng = 181;
  (void) fscanf(fp, "%d %d %s", &NumAng, &NumEn, DXSTypeName);
  if( !strcmp(DXSTypeName, "KTC") )     DXSType = 2;  // read DXS & calculate KT
  else if( !strcmp(DXSTypeName, "KT") ) DXSType = 1;  // read DXS & KT
  else                                  DXSType = 0;
 
  //  if( verboseLevel >= 1 ) G4cout << "Read DXS   (" << fileName  <<  ")\t NEg " << NumEn << " NAng " << NumAng
  //       << "DXSType " << DXSTypeName << " " << DXSType << G4endl;
 
  for (G4int eBin=1; eBin<=NumEn; eBin++){
    (void) fscanf(fp,"%f ",&data);
    Eb[eBin] = (G4double)data;
  }
 
 
  //for (aBin=1;aBin<NumAng;aBin++){
 
  if(DXSType==1) {
    G4cout << "DXSTYpe 1" << G4endl;
    for (G4int aBin=0;aBin<NumAng;aBin++){
      (void) fscanf(fp,"%f ",&data);
      DXS[0][aBin]=(G4double)data;
      for (G4int eBin=1;eBin<=NumEn;eBin++){
    (void) fscanf(fp,"%f %f ",&data2, &data);
    DXS[eBin][aBin]=(G4double)data;
    KT[eBin][aBin]=(G4double)data2;
      }
    }
  }
  else {
    for(G4int aBin=0; aBin<NumAng; aBin++){
      for(G4int eBin=0; eBin<=NumEn; eBin++){
    (void) fscanf(fp,"%f ",&data);
    DXS[eBin][aBin] = (G4double)data;
      }
    }
    for(G4int aBin=0; aBin<NumAng; aBin++){
      for(G4int eBin=1; eBin<=NumEn; eBin++){
    G4double A = DXS[0][aBin];                         // Angle
    G4double E = Eb[eBin];                             // Energy
    G4double p = sqrt(pow( (E/27.2/137),2) +2*E/27.2); // Momentum
    KT[eBin][aBin] = p *sqrt(2.-2.*cos(A*CLHEP::twopi/360.)); // Mom. Transfer
    //G4cout << "aEpKt " << aBin << " " << A << " E " << E << " p " << p << " KT "
    //   << KT[eBin][aBin] << " DXS " << DXS[eBin][aBin] << G4endl;
      }
    }
  }
 
  fclose(fp);
}

Referenced by G4LEPTSDiffXS().

SampleAngle()

G4double G4LEPTSDiffXS::SampleAngle

(

G4double

Energy

)

Definition at line 219 of file G4LEPTSDiffXS.cc.

                                                   {
  G4int  ii,jj,kk=0, Ebin;
 
  Ebin=1;
  for(ii=2; ii<=NumEn; ii++)
    if(Energy >= Eb[ii])
      Ebin=ii;
  if(Energy > Eb[NumEn]) Ebin=NumEn;
  else if(Energy > (Eb[Ebin]+Eb[Ebin+1])*0.5 ) Ebin++;
 
  //G4cout << "SampleAngle E " << Energy << " Ebin " << Ebin << " E[] " << Eb[Ebin] << G4endl;
 
  ii=0;
  jj=INumAng-1;
  G4double rnd=G4UniformRand();
 
  while ((jj-ii)>1){
    kk=(ii+jj)/2;
    G4double dxs = ICDXS[Ebin][kk];
    if (dxs < rnd) ii=kk;
    else           jj=kk;
  }
 
 
  //G4double x = ICDXS[0][jj];
  G4double x = ICDXS[0][kk] *CLHEP::twopi/360.;
 
  return(x);
}

Referenced by SampleAngleEthylene().

SampleAngleEthylene()

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

Definition at line 251 of file G4LEPTSDiffXS.cc.

                                                                   {
 
  BuildCDXS(E, El);
  NormalizeCDXS();
  InterpolateCDXS();
 
  return( SampleAngle(E) );
}

SampleAngleMT()

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

Definition at line 263 of file G4LEPTSDiffXS.cc.

                                                                     {
  G4int  ii, jj, kk=0, Ebin, iMin, iMax;
 
  G4double Ei = Energy;
  G4double Ed = Energy - Elost;
  G4double Pi = sqrt( pow( (Ei/27.2/137),2) +2*Ei/27.2); //incidente
  G4double Pd = sqrt( pow( (Ed/27.2/137),2) +2*Ed/27.2); //dispersado
  G4double Kmin = Pi - Pd;
  G4double Kmax = Pi + Pd;
 
  if(Pd <= 1e-9 ) return (0.0);
 
 
  // locate Energy bin
  Ebin=1;
  for(ii=2; ii<=NumEn; ii++)
    if(Energy > Eb[ii]) Ebin=ii;
  if(Energy > Eb[NumEn]) Ebin=NumEn;
  else if(Energy > (Eb[Ebin]+Eb[Ebin+1])*0.5 ) Ebin++;
 
  //G4cout << "SampleAngle2 E " << Energy << " Ebin " << Ebin << " E[] " << Eb[Ebin] << G4endl;
 
  ii=0; jj=INumAng-1;
  while ((jj-ii)>1) {
    kk=(ii+jj)/2;
    if( IKT[Ebin][kk] < Kmin ) ii=kk;
    else                      jj=kk;
  }
  iMin = ii;
 
  ii=0; jj=INumAng-1;
  while ((jj-ii)>1) {
    kk=(ii+jj)/2;
    if( IKT[Ebin][kk] < Kmax ) ii=kk;
    else                      jj=kk;
  }
  iMax = ii;
 
 
  // r -> a + (b-a)*r = a*(1-r) + b*r
  G4double rnd = G4UniformRand();
  rnd = (1-rnd)*ICDXS[Ebin][iMin] + rnd*ICDXS[Ebin][iMax];
  //G4double rnd = (ICDXS[Ebin][iMax] - ICDXS[Ebin][iMin]) * G4UniformRand()
  //  +  ICDXS[Ebin][iMin];
 
  ii=0; jj=INumAng-1;
  while ((jj-ii)>1){
    kk=(ii+jj)/2;
    if( ICDXS[Ebin][kk] < rnd) ii=kk;
    else                      jj=kk;
  }
 
  //Sampled
  G4double KR = IKT[Ebin][kk];
 
  G4double co = (Pi*Pi + Pd*Pd - KR*KR) / (2*Pi*Pd); //cos ang. disp.
  if(co > 1) co =1;
  G4double x = acos(co); //*360/twopi;            //ang. dispers.
 
  // Elastic aprox.
  //x = 2*asin(KR/Pi/2)*360/twopi;
 
  return(x);
}

---

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

• G4LEPTSDiffXS.hh
• G4LEPTSDiffXS.cc