Ignore:
Timestamp:
Jun 18, 2010, 11:42:07 AM (14 years ago)
Author:
garnier
Message:

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/electromagnetic/standard/src/G4GoudsmitSaundersonMscModel.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4GoudsmitSaundersonMscModel.cc,v 1.20 2009/12/16 17:50:11 gunter Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4GoudsmitSaundersonMscModel.cc,v 1.24 2010/05/17 15:11:30 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5151//                     assuming the case of lambdan<1 as single scattering regime
    5252//                     tuning theta sampling for theta below the screening angle
     53// 08.02.2010 O.Kadri: bugfix in compound xsection calculation and small angle computation
     54//                     adding a rejection condition to hard collision angular sampling
     55//                     ComputeTruePathLengthLimit was taken from G4WentzelVIModel
     56// 26.03.2010 O.Kadri: direct xsection calculation not inverse of the inverse
     57//                     angular sampling without large angle rejection method
     58//                     longitudinal displacement is computed exactly from <z>
     59// 12.05.2010 O.Kadri: exchange between target and projectile has as a condition the particle type (e-/e-)
     60//                     some cleanup to minimize time consuming (adding lamdan12 & Qn12, changing the error to 1.0e-12 for scrA)
    5361//
    5462//
    … …  
    6068//Ref.4:Bielajew et al.,".....", NIMB 173 (2001) 332-343;
    6169//Ref.5:F. Salvat et al.,"ELSEPA--Dirac partial ...molecules", Comp.Phys.Comm.165 (2005) pp 157-190;
    62 //Ref.6:G4UrbanMscModel G4_v9.1Ref09;
    63 //Ref.7:G4eCoulombScatteringModel G4_v9.1Ref09.
    64 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    65 
     70//Ref.6:G4UrbanMscModel G4 9.2;
     71//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    6672#include "G4GoudsmitSaundersonMscModel.hh"
    6773#include "G4GoudsmitSaundersonTable.hh"
    … …  
    7783#include "G4PhysicsTable.hh"
    7884#include "Randomize.hh"
    79 #include "G4Poisson.hh"
    8085
    8186using namespace std;
    … …  
    128133                       G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
    129134{ 
    130   //Build cross section table : Taken from Ref.7
    131135  G4double cs=0.0;
    132136  G4double kinEnergy = kineticEnergy;
    … …  
    134138  if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
    135139
    136   G4double value0,value1;
    137   CalculateIntegrals(p,Z,kinEnergy,value0,value1);
     140  G4double cs0;
     141  CalculateIntegrals(p,Z,kinEnergy,cs0,cs);
    138142 
    139   if(value1 > 0.0) cs = 1./value1;
    140 
    141143  return cs;
    142144} 
    … …  
    148150{   
    149151  G4double kineticEnergy = dynParticle->GetKineticEnergy();
    150   if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)) return ;
    151 
    152   G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
    153   G4double phi1,phi2,cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
    154   G4double q1,Gamma,Eta,delta,nu,nu0,nu1,nu2;
     152  if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)||
     153     (tPathLength/tausmall < lambda1)) return ;
    155154
    156155  ///////////////////////////////////////////
    157   // Effective energy and path-length from Eq. 4.7.15+16 of Ref.4
     156  // Effective energy
    158157  G4double  eloss = theManager->GetEnergy(particle,tPathLength,currentCouple);
    159   if(eloss>0.5*kineticEnergy)eloss=kineticEnergy-eloss;//exchange possibility between target atomic e- and incident particle
     158  if(eloss>0.5*kineticEnergy)
     159   {if((dynParticle->GetCharge())==-eplus)eloss=kineticEnergy-eloss;//exchange between target and projectile if they are electrons
     160    else eloss=0.5*kineticEnergy;
     161   }
    160162  G4double ee       = kineticEnergy - 0.5*eloss;
    161163  G4double ttau     = ee/electron_mass_c2;
    … …  
    166168  kineticEnergy *= (1 - cst1);
    167169  ///////////////////////////////////////////
    168   // additivity rule for mixture and compound xsection calculation
     170  // additivity rule for mixture and compound xsection's
    169171  const G4Material* mat = currentCouple->GetMaterial();
     172  const G4ElementVector* theElementVector = mat->GetElementVector();
     173  const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
    170174  G4int nelm = mat->GetNumberOfElements();
    171   const G4ElementVector* theElementVector = mat->GetElementVector();
    172   const G4double* theFraction = mat->GetFractionVector();
    173   G4double atomPerVolume = mat->GetTotNbOfAtomsPerVolume();
    174   G4double llambda0=0.,llambda1=0.;
     175  G4double s0,s1;
     176  lambda0=0.;
    175177  for(G4int i=0;i<nelm;i++)
    176     {
    177       G4double l0,l1;
    178       CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,l0,l1);
    179       llambda0 += (theFraction[i]/l0);
    180       llambda1 += (theFraction[i]/l1);
     178    {
     179      CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,s0,s1);
     180      lambda0 += (theAtomNumDensityVector[i]*s0);
    181181    }
    182   if(llambda0>DBL_MIN) llambda0 =1./llambda0;
    183   if(llambda1>DBL_MIN) llambda1 =1./llambda1;
     182  if(lambda0>DBL_MIN) lambda0 =1./lambda0;
     183
     184// Newton-Raphson root's finding method of scrA from:
     185// Sig1(PWA)/Sig0(PWA)=g1=2*scrA*((1+scrA)*log(1+1/scrA)-1)
    184186  G4double g1=0.0;
    185   if(llambda1>DBL_MIN) g1 = llambda0/llambda1;
    186 
    187   G4double x1,x0;
    188   x0=g1/2.;
     187  if(lambda1>DBL_MIN) g1 = lambda0/lambda1;
     188
     189  G4double logx0,x1,delta;
     190  G4double x0=g1/2.;
    189191  do
    190192    { 
    191       x1 = x0-(x0*((1.+x0)*log(1.+1./x0)-1.0)-g1/2.)/( (1.+2.*x0)*log(1.+1./x0)-2.0);// x1=x0-f(x0)/f'(x0)
     193      logx0=std::log(1.+1./x0);
     194      x1 = x0-(x0*((1.+x0)*logx0-1.0)-g1/2.)/( (1.+2.*x0)*logx0-2.0);
    192195      delta = std::abs( x1 - x0 );   
    193       x0 = x1;  // new approximation becomes the old approximation for the next iteration
    194     } while (delta > 1e-10);
     196      x0 = x1;
     197    } while (delta > 1.0e-12);
    195198  G4double scrA = x1;
    196199
    197   G4double us=0.0,vs=0.0,ws=1.0,x_coord=0.0,y_coord=0.0,z_coord=1.0;
    198200  G4double lambdan=0.;
    199   G4bool mscatt=false,noscatt=false;
    200 
    201   if(llambda0>0.)lambdan=atomPerVolume*tPathLength/llambda0;
    202   if((lambdan<=1.0e-12))return;
    203 
     201
     202  if(lambda0>0.)lambdan=tPathLength/lambda0;
     203  if(lambdan<=1.0e-12)return;
     204  G4double lambdan12=0.5*lambdan;
     205  Qn1 = lambdan *g1;//2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
     206  Qn12 = 0.5*Qn1;
     207 
     208  G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
     209  G4double cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
     210  G4double us=0.0,vs=0.0,ws=1.0,wss=0.,x_coord=0.0,y_coord=0.0,z_coord=1.0;
     211 
    204212  G4double epsilon1=G4UniformRand();
    205   G4double expn = exp(-lambdan);
    206   if((epsilon1<expn)||insideskin)// no scattering
    207     {noscatt=true;}
    208   else if((epsilon1<((1.+lambdan)*expn)||(lambdan<1.)))
    209     {
    210       mscatt=false;
    211       ws=G4UniformRand();
    212       ws= 1.-2.*scrA*ws/(1.-ws + scrA);
    213       if(acos(ws)<sqrt(scrA))//small angle approximation for theta less than screening angle
    214       {G4int i=0;
    215       do{i++;
    216       ws=1.+0.5*atomPerVolume*tPathLength*log(G4UniformRand())/llambda1;
    217       }while((fabs(ws)>1.)&&(i<20));//i<20 to avoid time consuming during the run
    218       if(i==19)ws=cos(sqrt(scrA));
    219       }
    220       G4double phi0=twopi*G4UniformRand();
    221       us=sqrt(1.-ws*ws)*cos(phi0);
    222       vs=sqrt(1.-ws*ws)*sin(phi0);
    223       G4double rr=G4UniformRand();
    224       x_coord=(rr*us);
    225       y_coord=(rr*vs);
    226       z_coord=((1.-rr)+rr*ws);
    227     }
    228   else
    229     {
    230       mscatt=true;
     213  G4double expn = std::exp(-lambdan);
     214  if(epsilon1<expn)// no scattering
     215    {return;}
     216  else if((epsilon1<((1.+lambdan)*expn))||(lambdan<1.))//single scattering (Rutherford DCS's)
     217    {
     218
     219      G4double xi=G4UniformRand();
     220      xi= 2.*scrA*xi/(1.-xi + scrA);   
     221      if(xi<0.)xi=0.;
     222      else if(xi>2.)xi=2.;
     223      ws=1.-xi;
     224      wss=std::sqrt(xi*(2.-xi));     
     225      G4double phi0=CLHEP::twopi*G4UniformRand();
     226      us=wss*cos(phi0);
     227      vs=wss*sin(phi0);
     228    }
     229  else // multiple scattering
     230    {
    231231      // Ref.2 subsection 4.4 "The best solution found"
    232232      // Sample first substep scattering angle
    233       SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
    234       phi1  = twopi*G4UniformRand();
     233      SampleCosineTheta(lambdan12,scrA,cosTheta1,sinTheta1);
     234      G4double phi1  = CLHEP::twopi*G4UniformRand();
    235235      cosPhi1 = cos(phi1);
    236236      sinPhi1 = sin(phi1);
    237237
    238238      // Sample second substep scattering angle
    239       SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
    240       phi2  = twopi*G4UniformRand();
     239      SampleCosineTheta(lambdan12,scrA,cosTheta2,sinTheta2);
     240      G4double phi2  = CLHEP::twopi*G4UniformRand();
    241241      cosPhi2 = cos(phi2);
    242242      sinPhi2 = sin(phi2);
    243243
    244       // Scattering direction
     244      // Overall scattering direction
    245245      us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
    246246      vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
    247247      ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2;
     248      G4double sqrtA=sqrt(scrA);
     249      if(acos(ws)<sqrtA)//small angle approximation for theta less than screening angle
     250      {
     251       G4int i=0;
     252       do{i++;
     253       ws=1.+Qn12*log(G4UniformRand());
     254       }while((fabs(ws)>1.)&&(i<20));//i<20 to avoid time consuming during the run
     255       if(i>=19)ws=cos(sqrtA);
     256
     257       wss=std::sqrt((1.-ws)*(1.0+ws));     
     258       us=wss*cos(phi1);
     259       vs=wss*sin(phi1);
     260     }
    248261    }
    249262   
    … …  
    253266  fParticleChange->ProposeMomentumDirection(newDirection);
    254267 
    255   if((safety > tlimitminfix)&&(latDisplasment))
    256     { 
    257 
    258       if(mscatt)
    259         {
    260           if(scrA<DBL_MIN)scrA=DBL_MIN;
    261           if(llambda1<DBL_MIN)llambda1=DBL_MIN;
    262           q1 = 2.*scrA*((1. + scrA)*log(1. + 1./scrA) - 1.);
    263           if(q1<DBL_MIN)q1=DBL_MIN;
    264           Gamma  = 6.*scrA*(1. + scrA)*((1. + 2.*scrA)*log(1. + 1./scrA) - 2.)/q1;
    265           Eta    = atomPerVolume*tPathLength/llambda1;     
    266           delta  = 0.90824829 - Eta*(0.102062073-Gamma*0.026374715);
    267 
    268           nu = G4UniformRand();
    269           nu = std::sqrt(nu);
    270           nu0 = (1.0 - nu)/2.;
    271           nu1 = nu*delta;
    272           nu2 = nu*(1.0-delta);
    273           x_coord=(nu1*sinTheta1*cosPhi1+nu2*sinTheta2*(cosPhi1*cosPhi2-cosTheta1*sinPhi1*sinPhi2)+nu0*us);
    274           y_coord=(nu1*sinTheta1*sinPhi1+nu2*sinTheta2*(sinPhi1*cosPhi2+cosTheta1*cosPhi1*sinPhi2)+nu0*vs);
    275           z_coord=(nu0+nu1*cosTheta1+nu2*cosTheta2+ nu0*ws)  ;
    276         }
    277 
     268  if((safety > tlimitminfix)&&latDisplasment)
     269    {
     270      if(Qn1<0.02)// corresponding to error less than 1% in the exact formula of <z>
     271      z_coord = 1.0 - Qn1*(0.5 - Qn1/6.);
     272      else z_coord = (1.-std::exp(-Qn1))/Qn1;
     273
     274      G4double rr=std::sqrt((1.- z_coord*z_coord)/(1.-ws*ws));
     275      x_coord = rr*us;
     276      y_coord = rr*vs;
    278277      // displacement is computed relatively to the end point
    279       if(!noscatt)z_coord -= 1.0;//for noscatt zcoord z_coord !=0.
    280       G4double rr = sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
     278      z_coord -= 1.0;
     279      rr = std::sqrt(x_coord*x_coord+y_coord*y_coord+z_coord*z_coord);
    281280      G4double r  = rr*zPathLength;
    282281      /*
    … …  
    287286             << G4endl;
    288287      */
     288      if(tPathLength<=zPathLength)return;
    289289      if(r > tlimitminfix) {
    290290
    291         G4ThreeVector latDirection(x_coord/rr,y_coord/rr,z_coord/rr);
    292         latDirection.rotateUz(oldDirection);
    293 
    294         ComputeDisplacement(fParticleChange, latDirection, r, safety);
     291        G4ThreeVector Direction(x_coord/rr,y_coord/rr,z_coord/rr);
     292        Direction.rotateUz(oldDirection);
     293
     294        ComputeDisplacement(fParticleChange, Direction, r, safety);
    295295      }     
    296296    }
    … …  
    303303                                                G4double &cost, G4double &sint)
    304304{
    305   G4double u,Qn1,r1,tet;
    306305  G4double xi=0.;
    307   Qn1=2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
    308 
    309 if (Qn1<0.001)xi=-0.5*Qn1*log(G4UniformRand());
    310 else if(Qn1>0.5)xi=2.*G4UniformRand();//isotropic distribution
    311 else
    312 {
    313       // procedure described by Benedito in Ref.1
     306 
     307  if (Qn12<0.001) 
     308  {G4double r1,tet;
    314309      do{
    315         r1=G4UniformRand();
    316         u=GSTable->SampleTheta(lambdan,scrA,G4UniformRand());
    317         xi = 2.*u;
    318         tet=acos(1.-xi);
     310        r1=G4UniformRand();
     311        xi=-Qn12*log(G4UniformRand());
     312        tet=acos(1.-xi);
    319313      }while(tet*r1*r1>sin(tet));
    320 }
     314  }
     315  else if(Qn12>0.5)xi=2.*G4UniformRand();
     316  else xi=2.*(GSTable->SampleTheta(lambdan,scrA,G4UniformRand()));
     317
    321318
    322319  if(xi<0.)xi=0.;
    323   if(xi>2.)xi=2.;
     320  else if(xi>2.)xi=2.;
     321
    324322  cost=(1. - xi);
    325323  sint=sqrt(xi*(2.-xi));
    … …  
    333331void
    334332G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z,
    335                                                  G4double kinEnergy,G4double &Lam0,
    336                                                  G4double &Lam1)
     333                                                 G4double kinEnergy,G4double &Sig0,
     334                                                 G4double &Sig1)
    337335{
    338   G4double summ00=0.0;
    339   G4double summ10=0.0;
    340336  G4double x1,x2,y1,y2,acoeff,bcoeff;
    341337  G4double kineticE = kinEnergy;
    … …  
    343339  if(kineticE>highKEnergy)kineticE=highKEnergy;
    344340  kineticE /= eV;
    345   G4double logE=log(kineticE);
     341  G4double logE=std::log(kineticE);
    346342 
    347343  G4int  iZ = G4int(Z);
    … …  
    349345
    350346  G4int enerInd=0;
    351   for(G4int i=1;i<106;i++)
     347  for(G4int i=0;i<105;i++)
    352348  {
    353349  if((logE>=ener[i])&&(logE<ener[i+1])){enerInd=i;break;}
    … …  
    364360        acoeff=(y2-y1)/(x2*x2-x1*x1);
    365361        bcoeff=y2-acoeff*x2*x2;
    366         summ00=acoeff*logE*logE+bcoeff;
    367         summ00 =exp(summ00);
     362        Sig0=acoeff*logE*logE+bcoeff;
     363        Sig0 =std::exp(Sig0);
    368364        y1=FTCSE[iZ-1][enerInd];
    369365        y2=FTCSE[iZ-1][enerInd+1];
    370366        acoeff=(y2-y1)/(x2*x2-x1*x1);
    371367        bcoeff=y2-acoeff*x2*x2;
    372         summ10=acoeff*logE*logE+bcoeff;
    373         summ10 =exp(summ10);
     368        Sig1=acoeff*logE*logE+bcoeff;
     369        Sig1=std::exp(Sig1);
    374370      }
    375371    else  //Interpolation of the form y=ax+b
    … …  
    379375        y1=TCSE[iZ-1][104];
    380376        y2=TCSE[iZ-1][105];
    381         summ00=(y2-y1)*(logE-x1)/(x2-x1)+y1;
    382         summ00 =exp(summ00);
     377        Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
     378        Sig0=std::exp(Sig0);
    383379        y1=FTCSE[iZ-1][104];
    384380        y2=FTCSE[iZ-1][105];
    385         summ10=(y2-y1)*(logE-x1)/(x2-x1)+y1;
    386         summ10 =exp(summ10);
     381        Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
     382        Sig1=std::exp(Sig1);
    387383      }
    388384    }
    … …  
    397393        acoeff=(y2-y1)/(x2*x2-x1*x1);
    398394        bcoeff=y2-acoeff*x2*x2;
    399         summ00=acoeff*logE*logE+bcoeff;
    400         summ00 =exp(summ00);
     395        Sig0=acoeff*logE*logE+bcoeff;
     396        Sig0 =std::exp(Sig0);
    401397        y1=FTCSP[iZ-1][enerInd];
    402398        y2=FTCSP[iZ-1][enerInd+1];
    403399        acoeff=(y2-y1)/(x2*x2-x1*x1);
    404400        bcoeff=y2-acoeff*x2*x2;
    405         summ10=acoeff*logE*logE+bcoeff;
    406         summ10 =exp(summ10);
     401        Sig1=acoeff*logE*logE+bcoeff;
     402        Sig1=std::exp(Sig1);
    407403      }
    408404    else
    … …  
    412408        y1=TCSP[iZ-1][104];
    413409        y2=TCSP[iZ-1][105];
    414         summ00=(y2-y1)*(logE-x1)/(x2-x1)+y1;
    415         summ00 =exp(summ00);
     410        Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
     411        Sig0 =std::exp(Sig0);
    416412        y1=FTCSP[iZ-1][104];
    417413        y2=FTCSP[iZ-1][105];
    418         summ10=(y2-y1)*(logE-x1)/(x2-x1)+y1;
    419         summ10 =exp(summ10);
     414        Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
     415        Sig1=std::exp(Sig1);
    420416      }
    421417    }
    422418   
    423   summ00 *=barn;
    424   summ10 *=barn;
    425 
    426   Lam0=1./((1.+1./Z)*summ00);
    427   Lam1=1./((1.+1./Z)*summ10);
    428 
    429 }
    430 
    431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    432 //t->g->t step transformations taken from Ref.6
     419  Sig0 *= barn;
     420  Sig1 *= barn;
     421
     422}
     423
     424//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     425//t->g->t step transformations taken from Ref.6
    433426
    434427G4double
    … …  
    635628
    636629//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    637 
     630// taken from Ref.6
    638631G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double)
    639632{
    … …  
    688681  if(samplez) {
    689682
    690     const G4double  ztmax = 0.99, onethird = 1./3. ;
     683    const G4double  ztmax = 0.99;
    691684    G4double zt = zmean/tPathLength ;
    692685
    … …  
    694687
    695688      G4double u,cz1;
    696       if(zt >= onethird) {
     689      if(zt >= 0.333333333) {
    697690
    698691        G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
    … …  
    719712
    720713//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    721 
     714// taken from Ref.6
    722715G4double
    723716G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength)
Note: See TracChangeset for help on using the changeset viewer.