Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-01-18 09:58:06

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 // Author: V. Grichine (Vladimir,Grichine@cern.ch)
0029 //
0030 //
0031 // G4 Model: diffuse optical elastic scattering with 4-momentum balance 
0032 //
0033 // Class Description
0034 // Final state production model for hadron nuclear elastic scattering; 
0035 // Class Description - End
0036 //
0037 //
0038 // 24.05.07 V. Grichine, first implementation for hadron (no Coulomb) elastic scattering
0039 // 04.09.07 V. Grichine, implementation for Coulomb elastic scattering
0040 // 12.06.11 V. Grichine, new interface to G4hadronElastic
0041 
0042 #ifndef G4DiffuseElastic_h
0043 #define G4DiffuseElastic_h 1
0044 
0045 #include <CLHEP/Units/PhysicalConstants.h>
0046 #include "globals.hh"
0047 #include "G4HadronElastic.hh"
0048 #include "G4HadProjectile.hh"
0049 #include "G4Nucleus.hh"
0050 
0051 #include "G4Pow.hh"
0052 
0053 
0054 class G4ParticleDefinition;
0055 class G4PhysicsTable;
0056 class G4PhysicsLogVector;
0057 
0058 class G4DiffuseElastic : public G4HadronElastic // G4HadronicInteraction
0059 {
0060 public:
0061 
0062   G4DiffuseElastic();
0063 
0064   // G4DiffuseElastic(const G4ParticleDefinition* aParticle);
0065 
0066 
0067 
0068 
0069 
0070   virtual ~G4DiffuseElastic();
0071 
0072   virtual G4bool IsApplicable(const G4HadProjectile &/*aTrack*/, 
0073                   G4Nucleus & /*targetNucleus*/);
0074 
0075   void Initialise();
0076 
0077   void InitialiseOnFly(G4double Z, G4double A);
0078 
0079   void BuildAngleTable();
0080 
0081  
0082   // G4HadFinalState* ApplyYourself(const G4HadProjectile & aTrack, G4Nucleus & targetNucleus);
0083 
0084   virtual G4double SampleInvariantT(const G4ParticleDefinition* p, 
0085                     G4double plab,
0086                     G4int Z, G4int A);
0087 
0088   G4double NeutronTuniform(G4int Z);
0089 
0090   void SetPlabLowLimit(G4double value);
0091 
0092   void SetHEModelLowLimit(G4double value);
0093 
0094   void SetQModelLowLimit(G4double value);
0095 
0096   void SetLowestEnergyLimit(G4double value);
0097 
0098   void SetRecoilKinEnergyLimit(G4double value);
0099 
0100   G4double SampleT(const G4ParticleDefinition* aParticle, 
0101                          G4double p, G4double A);
0102 
0103   G4double SampleTableT(const G4ParticleDefinition* aParticle, 
0104                          G4double p, G4double Z, G4double A);
0105 
0106   G4double SampleThetaCMS(const G4ParticleDefinition* aParticle, G4double p, G4double A);
0107 
0108   G4double SampleTableThetaCMS(const G4ParticleDefinition* aParticle, G4double p, 
0109                                      G4double Z, G4double A);
0110 
0111   G4double GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position);
0112 
0113   G4double SampleThetaLab(const G4HadProjectile* aParticle, 
0114                                 G4double tmass, G4double A);
0115 
0116   G4double GetDiffuseElasticXsc( const G4ParticleDefinition* particle, 
0117                                  G4double theta, 
0118                      G4double momentum, 
0119                  G4double A         );
0120 
0121   G4double GetInvElasticXsc( const G4ParticleDefinition* particle, 
0122                                  G4double theta, 
0123                      G4double momentum, 
0124                  G4double A, G4double Z );
0125 
0126   G4double GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle, 
0127                                  G4double theta, 
0128                      G4double momentum, 
0129                  G4double A, G4double Z );
0130 
0131   G4double GetInvElasticSumXsc( const G4ParticleDefinition* particle, 
0132                                  G4double tMand, 
0133                      G4double momentum, 
0134                  G4double A, G4double Z );
0135 
0136   G4double IntegralElasticProb( const G4ParticleDefinition* particle, 
0137                                  G4double theta, 
0138                      G4double momentum, 
0139                  G4double A            );
0140   
0141 
0142   G4double GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
0143                                  G4double theta, 
0144                      G4double momentum, 
0145                  G4double Z         );
0146 
0147   G4double GetInvCoulombElasticXsc( const G4ParticleDefinition* particle, 
0148                                  G4double tMand, 
0149                      G4double momentum, 
0150                  G4double A, G4double Z         );
0151 
0152   G4double GetCoulombTotalXsc( const G4ParticleDefinition* particle,  
0153                      G4double momentum, G4double Z       );
0154 
0155   G4double GetCoulombIntegralXsc( const G4ParticleDefinition* particle,  
0156                      G4double momentum, G4double Z, 
0157                                  G4double theta1, G4double theta2         );
0158 
0159 
0160   G4double CalculateParticleBeta( const G4ParticleDefinition* particle, 
0161                                     G4double momentum    );
0162 
0163   G4double CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 );
0164 
0165   G4double CalculateAm( G4double momentum, G4double n, G4double Z);
0166 
0167   G4double CalculateNuclearRad( G4double A);
0168 
0169   G4double ThetaCMStoThetaLab(const G4DynamicParticle* aParticle, 
0170                                 G4double tmass, G4double thetaCMS);
0171 
0172   G4double ThetaLabToThetaCMS(const G4DynamicParticle* aParticle, 
0173                                 G4double tmass, G4double thetaLab);
0174 
0175   void TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
0176               G4double Z, G4double A);
0177 
0178 
0179 
0180   G4double BesselJzero(G4double z);
0181   G4double BesselJone(G4double z);
0182   G4double DampFactor(G4double z);
0183   G4double BesselOneByArg(G4double z);
0184 
0185   G4double GetDiffElasticProb(G4double theta);
0186   G4double GetDiffElasticSumProb(G4double theta);
0187   G4double GetDiffElasticSumProbA(G4double alpha);
0188   G4double GetIntegrandFunction(G4double theta);
0189 
0190 
0191   G4double GetNuclearRadius(){return fNuclearRadius;};
0192 
0193 private:
0194 
0195 
0196   G4ParticleDefinition* theProton;
0197   G4ParticleDefinition* theNeutron;
0198   G4ParticleDefinition* theDeuteron;
0199   G4ParticleDefinition* theAlpha;
0200 
0201   const G4ParticleDefinition* thePionPlus;
0202   const G4ParticleDefinition* thePionMinus;
0203 
0204   G4double lowEnergyRecoilLimit;  
0205   G4double lowEnergyLimitHE;  
0206   G4double lowEnergyLimitQ;  
0207   G4double lowestEnergyLimit;  
0208   G4double plabLowLimit;
0209 
0210   G4int fEnergyBin;
0211   G4int fAngleBin;
0212 
0213   G4PhysicsLogVector*           fEnergyVector;
0214   G4PhysicsTable*               fAngleTable;
0215   std::vector<G4PhysicsTable*>  fAngleBank;
0216 
0217   std::vector<G4double> fElementNumberVector;
0218   std::vector<G4String> fElementNameVector;
0219 
0220   const G4ParticleDefinition* fParticle;
0221   G4double fWaveVector;
0222   G4double fAtomicWeight;
0223   G4double fAtomicNumber;
0224   G4double fNuclearRadius;
0225   G4double fBeta;
0226   G4double fZommerfeld;
0227   G4double fAm;
0228   G4bool fAddCoulomb;
0229 
0230 };
0231 
0232 inline G4bool G4DiffuseElastic::IsApplicable(const G4HadProjectile & projectile, 
0233                   G4Nucleus & nucleus)
0234 {
0235   if( ( projectile.GetDefinition() == G4Proton::Proton() ||
0236         projectile.GetDefinition() == G4Neutron::Neutron() ||
0237         projectile.GetDefinition() == G4PionPlus::PionPlus() ||
0238         projectile.GetDefinition() == G4PionMinus::PionMinus() ||
0239         projectile.GetDefinition() == G4KaonPlus::KaonPlus() ||
0240         projectile.GetDefinition() == G4KaonMinus::KaonMinus() ) &&
0241 
0242         nucleus.GetZ_asInt() >= 2 ) return true;
0243   else                              return false;
0244 }
0245 
0246 inline void G4DiffuseElastic::SetRecoilKinEnergyLimit(G4double value)
0247 {
0248   lowEnergyRecoilLimit = value;
0249 }
0250 
0251 inline void G4DiffuseElastic::SetPlabLowLimit(G4double value)
0252 {
0253   plabLowLimit = value;
0254 }
0255 
0256 inline void G4DiffuseElastic::SetHEModelLowLimit(G4double value)
0257 {
0258   lowEnergyLimitHE = value;
0259 }
0260 
0261 inline void G4DiffuseElastic::SetQModelLowLimit(G4double value)
0262 {
0263   lowEnergyLimitQ = value;
0264 }
0265 
0266 inline void G4DiffuseElastic::SetLowestEnergyLimit(G4double value)
0267 {
0268   lowestEnergyLimit = value;
0269 }
0270 
0271 
0272 /////////////////////////////////////////////////////////////
0273 //
0274 // Bessel J0 function based on rational approximation from 
0275 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 
0276 
0277 inline G4double G4DiffuseElastic::BesselJzero(G4double value)
0278 {
0279   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
0280 
0281   modvalue = std::fabs(value);
0282 
0283   if ( value < 8.0 && value > -8.0 )
0284   {
0285     value2 = value*value;
0286 
0287     fact1  = 57568490574.0 + value2*(-13362590354.0 
0288                            + value2*( 651619640.7 
0289                            + value2*(-11214424.18 
0290                            + value2*( 77392.33017 
0291                            + value2*(-184.9052456   ) ) ) ) );
0292 
0293     fact2  = 57568490411.0 + value2*( 1029532985.0 
0294                            + value2*( 9494680.718
0295                            + value2*(59272.64853
0296                            + value2*(267.8532712 
0297                            + value2*1.0               ) ) ) );
0298 
0299     bessel = fact1/fact2;
0300   } 
0301   else 
0302   {
0303     arg    = 8.0/modvalue;
0304 
0305     value2 = arg*arg;
0306 
0307     shift  = modvalue-0.785398164;
0308 
0309     fact1  = 1.0 + value2*(-0.1098628627e-2 
0310                  + value2*(0.2734510407e-4
0311                  + value2*(-0.2073370639e-5 
0312                  + value2*0.2093887211e-6    ) ) );
0313 
0314     fact2  = -0.1562499995e-1 + value2*(0.1430488765e-3
0315                               + value2*(-0.6911147651e-5 
0316                               + value2*(0.7621095161e-6
0317                               - value2*0.934945152e-7    ) ) );
0318 
0319     bessel = std::sqrt(0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2 );
0320   }
0321   return bessel;
0322 }
0323 
0324 /////////////////////////////////////////////////////////////
0325 //
0326 // Bessel J1 function based on rational approximation from 
0327 // J.F. Hart, Computer Approximations, New York, Willey 1968, p. 141 
0328 
0329 inline G4double G4DiffuseElastic::BesselJone(G4double value)
0330 {
0331   G4double modvalue, value2, fact1, fact2, arg, shift, bessel;
0332 
0333   modvalue = std::fabs(value);
0334 
0335   if ( modvalue < 8.0 ) 
0336   {
0337     value2 = value*value;
0338 
0339     fact1  = value*(72362614232.0 + value2*(-7895059235.0 
0340                                   + value2*( 242396853.1
0341                                   + value2*(-2972611.439 
0342                                   + value2*( 15704.48260 
0343                                   + value2*(-30.16036606  ) ) ) ) ) );
0344 
0345     fact2  = 144725228442.0 + value2*(2300535178.0 
0346                             + value2*(18583304.74
0347                             + value2*(99447.43394 
0348                             + value2*(376.9991397 
0349                             + value2*1.0             ) ) ) );
0350     bessel = fact1/fact2;
0351   } 
0352   else 
0353   {
0354     arg    = 8.0/modvalue;
0355 
0356     value2 = arg*arg;
0357 
0358     shift  = modvalue - 2.356194491;
0359 
0360     fact1  = 1.0 + value2*( 0.183105e-2 
0361                  + value2*(-0.3516396496e-4
0362                  + value2*(0.2457520174e-5 
0363                  + value2*(-0.240337019e-6          ) ) ) );
0364 
0365     fact2 = 0.04687499995 + value2*(-0.2002690873e-3
0366                           + value2*( 0.8449199096e-5
0367                           + value2*(-0.88228987e-6
0368                           + value2*0.105787412e-6       ) ) );
0369 
0370     bessel = std::sqrt( 0.636619772/modvalue)*(std::cos(shift)*fact1 - arg*std::sin(shift)*fact2);
0371 
0372     if (value < 0.0) bessel = -bessel;
0373   }
0374   return bessel;
0375 }
0376 
0377 ////////////////////////////////////////////////////////////////////
0378 //
0379 // damp factor in diffraction x/sh(x), x was already *pi
0380 
0381 inline G4double G4DiffuseElastic::DampFactor(G4double x)
0382 {
0383   G4double df;
0384   G4double f2 = 2., f3 = 6., f4 = 24.; // first factorials
0385 
0386   // x *= pi;
0387 
0388   if( std::fabs(x) < 0.01 )
0389   { 
0390     df = 1./(1. + x/f2 + x*x/f3 + x*x*x/f4);
0391   }
0392   else
0393   {
0394     df = x/std::sinh(x); 
0395   }
0396   return df;
0397 }
0398 
0399 
0400 ////////////////////////////////////////////////////////////////////
0401 //
0402 // return J1(x)/x with special case for small x
0403 
0404 inline G4double G4DiffuseElastic::BesselOneByArg(G4double x)
0405 {
0406   G4double x2, result;
0407   
0408   if( std::fabs(x) < 0.01 )
0409   { 
0410    x     *= 0.5;
0411    x2     = x*x;
0412    result = 2. - x2 + x2*x2/6.;
0413   }
0414   else
0415   {
0416     result = BesselJone(x)/x; 
0417   }
0418   return result;
0419 }
0420 
0421 ////////////////////////////////////////////////////////////////////
0422 //
0423 // return particle beta
0424 
0425 inline  G4double G4DiffuseElastic::CalculateParticleBeta( const G4ParticleDefinition* particle, 
0426                                     G4double momentum    )
0427 {
0428   G4double mass = particle->GetPDGMass();
0429   G4double a    = momentum/mass;
0430   fBeta         = a/std::sqrt(1+a*a);
0431 
0432   return fBeta; 
0433 }
0434 
0435 ////////////////////////////////////////////////////////////////////
0436 //
0437 // return Zommerfeld parameter for Coulomb scattering
0438 
0439 inline  G4double G4DiffuseElastic::CalculateZommerfeld( G4double beta, G4double Z1, G4double Z2 )
0440 {
0441   fZommerfeld = CLHEP::fine_structure_const*Z1*Z2/beta;
0442 
0443   return fZommerfeld; 
0444 }
0445 
0446 ////////////////////////////////////////////////////////////////////
0447 //
0448 // return Wentzel correction for Coulomb scattering
0449 
0450 inline  G4double G4DiffuseElastic::CalculateAm( G4double momentum, G4double n, G4double Z)
0451 {
0452   G4double k   = momentum/CLHEP::hbarc;
0453   G4double ch  = 1.13 + 3.76*n*n;
0454   G4double zn  = 1.77*k*(1.0/G4Pow::GetInstance()->A13(Z))*CLHEP::Bohr_radius;
0455   G4double zn2 = zn*zn;
0456   fAm          = ch/zn2;
0457 
0458   return fAm;
0459 }
0460 
0461 ////////////////////////////////////////////////////////////////////
0462 //
0463 // calculate nuclear radius for different atomic weights using different approximations
0464 
0465 inline  G4double G4DiffuseElastic::CalculateNuclearRad( G4double A)
0466 {
0467   G4double R, r0, a11, a12, a13, a2, a3;
0468 
0469   a11 = 1.26;  // 1.08, 1.16
0470   a12 = 1.;  // 1.08, 1.16
0471   a13 = 1.12;  // 1.08, 1.16
0472   a2 = 1.1;
0473   a3 = 1.;
0474 
0475   // Special rms radii for light nucleii
0476 
0477   if (A < 50.)
0478   {
0479     if     (std::abs(A-1.) < 0.5)                         return 0.89*CLHEP::fermi; // p
0480     else if(std::abs(A-2.) < 0.5)                         return 2.13*CLHEP::fermi; // d
0481     else if(  // std::abs(Z-1.) < 0.5 && 
0482 std::abs(A-3.) < 0.5) return 1.80*CLHEP::fermi; // t
0483 
0484     // else if(std::abs(Z-2.) < 0.5 && std::abs(A-3.) < 0.5) return 1.96CLHEP::fermi; // He3
0485     else if( // std::abs(Z-2.) < 0.5 && 
0486 std::abs(A-4.) < 0.5) return 1.68*CLHEP::fermi; // He4
0487 
0488     else if(  // std::abs(Z-3.) < 0.5
0489         std::abs(A-7.) < 0.5   )                         return 2.40*CLHEP::fermi; // Li7
0490     else if(  // std::abs(Z-4.) < 0.5 
0491 std::abs(A-9.) < 0.5)                         return 2.51*CLHEP::fermi; // Be9
0492 
0493     else if( 10.  < A && A <= 16. ) r0  = a11*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;   // 1.08CLHEP::fermi;
0494     else if( 15.  < A && A <= 20. ) r0  = a12*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
0495     else if( 20.  < A && A <= 30. ) r0  = a13*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;
0496     else                            r0  = a2*CLHEP::fermi;
0497 
0498     R = r0*G4Pow::GetInstance()->A13(A);
0499   }
0500   else
0501   {
0502     r0 = a3*CLHEP::fermi;
0503 
0504     R  = r0*G4Pow::GetInstance()->powA(A, 0.27);
0505   }
0506   fNuclearRadius = R;
0507   return R;
0508   /*
0509   G4double r0;
0510   if( A < 50. )
0511   {
0512     if( A > 10. ) r0  = 1.16*( 1 - (1.0/G4Pow::GetInstance()->A23(A)) )*CLHEP::fermi;   // 1.08CLHEP::fermi;
0513     else          r0  = 1.1*CLHEP::fermi;
0514     fNuclearRadius = r0*G4Pow::GetInstance()->A13(A);
0515   }
0516   else
0517   {
0518     r0 = 1.7*CLHEP::fermi;   // 1.7*CLHEP::fermi;
0519     fNuclearRadius = r0*G4Pow::GetInstance()->powA(A, 0.27); // 0.27);
0520   }
0521   return fNuclearRadius;
0522   */
0523 }
0524 
0525 ////////////////////////////////////////////////////////////////////
0526 //
0527 // return Coulomb scattering differential xsc with Wentzel correction  
0528 
0529 inline  G4double G4DiffuseElastic::GetCoulombElasticXsc( const G4ParticleDefinition* particle, 
0530                                  G4double theta, 
0531                      G4double momentum, 
0532                  G4double Z         )
0533 {
0534   G4double sinHalfTheta  = std::sin(0.5*theta);
0535   G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
0536   G4double beta          = CalculateParticleBeta( particle, momentum);
0537   G4double z             = particle->GetPDGCharge();
0538   G4double n             = CalculateZommerfeld( beta, z, Z );
0539   G4double am            = CalculateAm( momentum, n, Z);
0540   G4double k             = momentum/CLHEP::hbarc;
0541   G4double ch            = 0.5*n/k;
0542   G4double ch2           = ch*ch;
0543   G4double xsc           = ch2/(sinHalfTheta2+am)/(sinHalfTheta2+am);
0544 
0545   return xsc;
0546 }
0547 
0548 
0549 ////////////////////////////////////////////////////////////////////
0550 //
0551 // return Coulomb scattering total xsc with Wentzel correction  
0552 
0553 inline  G4double G4DiffuseElastic::GetCoulombTotalXsc( const G4ParticleDefinition* particle,  
0554                                                  G4double momentum, G4double Z  )
0555 {
0556   G4double beta          = CalculateParticleBeta( particle, momentum);
0557   G4cout<<"beta = "<<beta<<G4endl;
0558   G4double z             = particle->GetPDGCharge();
0559   G4double n             = CalculateZommerfeld( beta, z, Z );
0560   G4cout<<"fZomerfeld = "<<n<<G4endl;
0561   G4double am            = CalculateAm( momentum, n, Z);
0562   G4cout<<"cof Am = "<<am<<G4endl;
0563   G4double k             = momentum/CLHEP::hbarc;
0564   G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
0565   G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
0566   G4double ch            = n/k;
0567   G4double ch2           = ch*ch;
0568   G4double xsc           = ch2*CLHEP::pi/(am +am*am);
0569 
0570   return xsc;
0571 }
0572 
0573 ////////////////////////////////////////////////////////////////////
0574 //
0575 // return Coulomb scattering xsc with Wentzel correction  integrated between
0576 // theta1 and < theta2
0577 
0578 inline  G4double G4DiffuseElastic::GetCoulombIntegralXsc( const G4ParticleDefinition* particle,  
0579                      G4double momentum, G4double Z, 
0580                                  G4double theta1, G4double theta2 )
0581 {
0582   G4double c1 = std::cos(theta1);
0583   G4cout<<"c1 = "<<c1<<G4endl;
0584   G4double c2 = std::cos(theta2);
0585   G4cout<<"c2 = "<<c2<<G4endl;
0586   G4double beta          = CalculateParticleBeta( particle, momentum);
0587   // G4cout<<"beta = "<<beta<<G4endl;
0588   G4double z             = particle->GetPDGCharge();
0589   G4double n             = CalculateZommerfeld( beta, z, Z );
0590   // G4cout<<"fZomerfeld = "<<n<<G4endl;
0591   G4double am            = CalculateAm( momentum, n, Z);
0592   // G4cout<<"cof Am = "<<am<<G4endl;
0593   G4double k             = momentum/CLHEP::hbarc;
0594   // G4cout<<"k = "<<k*CLHEP::fermi<<" 1/fermi"<<G4endl;
0595   // G4cout<<"k*Bohr_radius = "<<k*CLHEP::Bohr_radius<<G4endl;
0596   G4double ch            = n/k;
0597   G4double ch2           = ch*ch;
0598   am *= 2.;
0599   G4double xsc           = ch2*CLHEP::twopi*(c1-c2);
0600            xsc          /= (1 - c1 + am)*(1 - c2 + am);
0601 
0602   return xsc;
0603 }
0604 
0605 #endif