Ignore:
Timestamp:
Nov 5, 2010, 3:45:55 PM (14 years ago)
Author:
garnier
Message:

update ti head

File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/source/processes/hadronic/cross_sections/src/G4GGNuclNuclCrossSection.cc

    r1228 r1340  
    3333#include "G4IonTable.hh"
    3434#include "G4ParticleDefinition.hh"
    35 
    36 
    37 
    38 ////////////////////////////////////////////////////////////////////////////////
     35#include "G4HadTmpUtil.hh"
     36
     37
     38///////////////////////////////////////////////////////////////////////////////
    3939//
    4040//
    … …  
    4949}
    5050
    51 ///////////////////////////////////////////////////////////////////////////////////////
     51///////////////////////////////////////////////////////////////////////////////
    5252//
    5353//
    5454
    5555G4GGNuclNuclCrossSection::~G4GGNuclNuclCrossSection()
    56 {
    57 }
    58 
    59 
    60 ////////////////////////////////////////////////////////////////////////////////////////
     56{}
     57
     58///////////////////////////////////////////////////////////////////////////////
    6159//
    6260//
    … …  
    6563G4bool
    6664G4GGNuclNuclCrossSection::IsApplicable(const G4DynamicParticle* aDP,
    67                                           const G4Element*  anElement)
    68 {
    69   return IsZAApplicable(aDP, anElement->GetZ(), anElement->GetN());
     65                                       const G4Element*  anElement)
     66{
     67  G4int Z = G4lrint(anElement->GetZ());
     68  G4int N = G4lrint(anElement->GetN());
     69  return IsIsoApplicable(aDP, Z, N);
    7070}
    7171
    72 ////////////////////////////////////////////////////////////////////////////////////////
     72///////////////////////////////////////////////////////////////////////////////
    7373//
    7474//
    7575
    7676G4bool
    77 G4GGNuclNuclCrossSection::IsZAApplicable(const G4DynamicParticle* aDP,
    78                                             G4double Z, G4double)
    79 {
    80   G4bool applicable      = false;
    81   // G4int baryonNumber     = aDP->GetDefinition()->GetBaryonNumber();
     77G4GGNuclNuclCrossSection::IsIsoApplicable(const G4DynamicParticle* aDP,
     78                                          G4int Z, G4int)
     79{
     80  G4bool applicable = false;
    8281  G4double kineticEnergy = aDP->GetKineticEnergy();
    8382
    84   //  const G4ParticleDefinition* theParticle = aDP->GetDefinition();
    85  
    86   if ( kineticEnergy  >= fLowerLimit && Z > 1.5 ) applicable = true;
    87 
     83  if (kineticEnergy >= fLowerLimit && Z > 1) applicable = true;
    8884  return applicable;
    8985}
    9086
    91 ////////////////////////////////////////////////////////////////////////////////////////
    92 //
    93 // Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
    94 // Glauber model with Gribov correction calculated in the dipole approximation on
    95 // light cone. Gaussian density helps to calculate rest integrals of the model.
    96 // [1] B.Z. Kopeliovich, nucl-th/0306044
     87///////////////////////////////////////////////////////////////////////////////
     88//
     89// Calculates total and inelastic Xsc, derives elastic as total - inelastic
     90// accordong to Glauber model with Gribov correction calculated in the dipole
     91// approximation on light cone. Gaussian density helps to calculate rest
     92// integrals of the model. [1] B.Z. Kopeliovich, nucl-th/0306044
    9793
    9894
    9995G4double G4GGNuclNuclCrossSection::
    100 GetCrossSection(const G4DynamicParticle* aParticle, const G4Element* anElement, G4double T)
    101 {
    102   return GetIsoZACrossSection(aParticle, anElement->GetZ(), anElement->GetN(), T);
    103 }
    104 
    105 ////////////////////////////////////////////////////////////////////////////////////////
    106 //
    107 // Calculates total and inelastic Xsc, derives elastic as total - inelastic accordong to
    108 // Glauber model with Gribov correction calculated in the dipole approximation on
    109 // light cone. Gaussian density of point-like nucleons helps to calculate rest integrals of the model.
    110 // [1] B.Z. Kopeliovich, nucl-th/0306044 + simplification above
    111 
     96GetCrossSection(const G4DynamicParticle* aParticle, const G4Element* anElement,
     97                G4double T)
     98{
     99  G4int Z = G4lrint(anElement->GetZ());
     100  G4int N = G4lrint(anElement->GetN());
     101  return GetZandACrossSection(aParticle, Z, N, T);
     102}
     103
     104///////////////////////////////////////////////////////////////////////////////
     105//
     106// Calculates total and inelastic Xsc, derives elastic as total - inelastic
     107// accordong to Glauber model with Gribov correction calculated in the dipole
     108// approximation on light cone. Gaussian density of point-like nucleons helps
     109// to calculate rest integrals of the model. [1] B.Z. Kopeliovich,
     110// nucl-th/0306044 + simplification above
    112111
    113112
    114113G4double G4GGNuclNuclCrossSection::
    115 GetIsoZACrossSection(const G4DynamicParticle* aParticle, G4double tZ, G4double tA, G4double)
    116 {
    117   G4double xsection, sigma, cofInelastic = 2.4, cofTotal = 2.0, nucleusSquare, cB, ratio;
     114GetZandACrossSection(const G4DynamicParticle* aParticle,
     115                     G4int tZ, G4int tA, G4double)
     116{
     117  G4double xsection;
     118  G4double sigma;
     119  G4double cofInelastic = 2.4;
     120  G4double cofTotal = 2.0;
     121  G4double nucleusSquare;
     122  G4double cB;
     123  G4double ratio;
    118124
    119125  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
    … …  
    132138  G4double pR = GetNucleusRadius(pA);
    133139
    134   cB = GetCoulombBarier(aParticle, tZ, tA, pR, tR);
    135 
    136   if(cB > 0.)
    137   {
     140  cB = GetCoulombBarier(aParticle, G4double(tZ), G4double(tA), pR, tR);
     141  if (cB > 0.) {
    138142
    139143    sigma = (pZ*tZ+pN*tN)*GetHadronNucleonXscNS(theProton, pTkin, theProton) +
    … …  
    143147
    144148    ratio = sigma/nucleusSquare;
    145 
    146149    xsection =  nucleusSquare*std::log( 1. + ratio );
    147 
    148150    fTotalXsc = xsection;
    149 
    150151    fTotalXsc *= cB;
    151152
    … …  
    153154
    154155    fInelasticXsc *= cB;
    155 
    156156    fElasticXsc   = fTotalXsc - fInelasticXsc;
    157157
    … …  
    171171 
    172172    ratio = sigma/nucleusSquare;
    173 
    174173    fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
    175174
    … …  
    186185}
    187186
    188 ////////////////////////////////////////////////////////////////////////////////////////
     187///////////////////////////////////////////////////////////////////////////////
    189188//
    190189//
    191190
    192191G4double G4GGNuclNuclCrossSection::
    193 GetCoulombBarier(const G4DynamicParticle* aParticle, G4double tZ, G4double tA, G4double pR, G4double tR)
     192GetCoulombBarier(const G4DynamicParticle* aParticle, G4double tZ, G4double tA,
     193                 G4double pR, G4double tR)
    194194{
    195195  G4double ratio;
    196196  G4double pZ = aParticle->GetDefinition()->GetPDGCharge();
    197 
    198197
    199198  G4double pTkin = aParticle->GetKineticEnergy();
    … …  
    254253
    255254  nucleusSquare = cofTotal*pi*( pR*pR + tR*tR );   // basically 2piRR
    256 
    257255  ratio = sigma/nucleusSquare;
    258 
    259 
    260   fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
    261    
     256  fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
    262257  G4double difratio = ratio/(1.+ratio);
    263258
    … …  
    272267//////////////////////////////////////////////////////////////////////////
    273268//
    274 // Return suasi-elastic/inelastic cross-section ratio
     269// Return quasi-elastic/inelastic cross-section ratio
    275270
    276271G4double G4GGNuclNuclCrossSection::
    … …  
    298293
    299294  nucleusSquare = cofTotal*pi*( pR*pR + tR*tR );   // basically 2piRR
    300 
    301295  ratio = sigma/nucleusSquare;
    302 
    303   fInelasticXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
     296  fInelasticXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
    304297
    305298  //  sigma = GetHNinelasticXsc(aParticle, tA, tZ);
    306299  ratio = sigma/nucleusSquare;
    307 
    308   fProductionXsc = nucleusSquare*std::log( 1. + cofInelastic*ratio )/cofInelastic;
     300  fProductionXsc = nucleusSquare*std::log(1. + cofInelastic*ratio)/cofInelastic;
    309301
    310302  if (fInelasticXsc > fProductionXsc) ratio = (fInelasticXsc-fProductionXsc)/fInelasticXsc;
    … …  
    315307}
    316308
    317 /////////////////////////////////////////////////////////////////////////////////////
     309///////////////////////////////////////////////////////////////////////////////
    318310//
    319311// Returns hadron-nucleon Xsc according to differnt parametrisations:
    … …  
    324316G4double
    325317G4GGNuclNuclCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
    326                                                   const G4Element* anElement          )
    327 {
    328   G4double At = anElement->GetN();  // number of nucleons
    329   G4double Zt = anElement->GetZ();  // number of protons
    330 
    331 
    332   return GetHadronNucleonXsc( aParticle, At, Zt );
    333 }
    334 
    335 
    336 
    337 
    338 /////////////////////////////////////////////////////////////////////////////////////
     318                                              const G4Element* anElement)
     319{
     320  G4int At = G4lrint(anElement->GetN());  // number of nucleons
     321  G4int Zt = G4lrint(anElement->GetZ());  // number of protons
     322  return GetHadronNucleonXsc(aParticle, At, Zt);
     323}
     324
     325
     326
     327
     328///////////////////////////////////////////////////////////////////////////////
    339329//
    340330// Returns hadron-nucleon Xsc according to differnt parametrisations:
    … …  
    345335G4double
    346336G4GGNuclNuclCrossSection::GetHadronNucleonXsc(const G4DynamicParticle* aParticle,
    347                                                    G4double At,  G4double Zt       )
     337                                                   G4int At, G4int Zt)
    348338{
    349339  G4double xsection = 0.;
    350340
    351 
    352341  G4double targ_mass = G4ParticleTable::GetParticleTable()->
    353   GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
    354 
     342  GetIonTable()->GetIonMass(Zt, At);
    355343  targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
    356344
    357   G4double proj_mass     = aParticle->GetMass();
     345  G4double proj_mass = aParticle->GetMass();
    358346  G4double proj_momentum = aParticle->GetMomentum().mag();
    359347  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
    … …  
    361349  sMand /= GeV*GeV;  // in GeV for parametrisation
    362350  proj_momentum /= GeV;
    363 
    364351  const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
    365  
    366352
    367353  if(pParticle == theNeutron) // as proton ???
    368354  {
    369     xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
     355    xsection = G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
    370356  }
    371357  else if(pParticle == theProton)
    372358  {
    373     xsection = At*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
    374     // xsection = At*( 49.51*std::pow(sMand,-0.097) + 0.314*std::log(sMand)*std::log(sMand) );
    375     // xsection = At*( 38.4 + 0.85*std::abs(std::pow(log(sMand),1.47)) );
     359    xsection = G4double(At)*(21.70*std::pow(sMand,0.0808) + 56.08*std::pow(sMand,-0.4525));
    376360  }
    377361 
    … …  
    381365
    382366
    383 /////////////////////////////////////////////////////////////////////////////////////
     367///////////////////////////////////////////////////////////////////////////////
    384368//
    385369// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
    … …  
    388372G4double
    389373G4GGNuclNuclCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
    390                                                   const G4Element* anElement          )
    391 {
    392   G4double At = anElement->GetN();  // number of nucleons
    393   G4double Zt = anElement->GetZ();  // number of protons
    394 
    395 
     374                                                  const G4Element* anElement)
     375{
     376  G4int At = G4lrint(anElement->GetN());  // number of nucleons
     377  G4int Zt = G4lrint(anElement->GetZ());  // number of protons
    396378  return GetHadronNucleonXscPDG( aParticle, At, Zt );
    397379}
    398380
    399381
    400 
    401 
    402 /////////////////////////////////////////////////////////////////////////////////////
     382///////////////////////////////////////////////////////////////////////////////
    403383//
    404384// Returns hadron-nucleon Xsc according to PDG parametrisation (2005):
    … …  
    408388G4double
    409389G4GGNuclNuclCrossSection::GetHadronNucleonXscPDG(const G4DynamicParticle* aParticle,
    410                                                      G4double At,  G4double Zt )
     390                                                 G4int At, G4int Zt)
    411391{
    412392  G4double xsection = 0.;
    … …  
    415395  if (Nt < 0.) Nt = 0.; 
    416396
    417 
    418397  G4double targ_mass = G4ParticleTable::GetParticleTable()->
    419   GetIonTable()->GetIonMass( G4int(Zt+0.5) , G4int(At+0.5) );
    420 
     398  GetIonTable()->GetIonMass(Zt, At);
    421399  targ_mass = 0.939*GeV;  // ~mean neutron and proton ???
    422400
    423401  G4double proj_mass     = aParticle->GetMass();
    424402  G4double proj_momentum = aParticle->GetMomentum().mag();
    425 
    426403  G4double sMand = CalcMandelstamS ( proj_mass , targ_mass , proj_momentum );
    427 
    428404  sMand         /= GeV*GeV;  // in GeV for parametrisation
    429405
    … …  
    435411  G4double B    = 0.308;
    436412
    437 
    438413  const G4ParticleDefinition* pParticle = aParticle->GetDefinition();
    439414 
    440 
    441415  if(pParticle == theNeutron) // proton-neutron fit
    442416  {
    443     xsection = Zt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
    444                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
    445     xsection  += Nt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
    446                       + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
     417    xsection = G4double(Zt)*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
     418                  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
     419    xsection += Nt*(35.45 + B*std::pow(std::log(sMand/s0),2.)
     420             + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2)); // pp for nn
    447421  }
    448422  else if(pParticle == theProton)
    449423  {
    450      
    451       xsection  = Zt*( 35.45 + B*std::pow(std::log(sMand/s0),2.)
    452                           + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
    453 
    454       xsection += Nt*( 35.80 + B*std::pow(std::log(sMand/s0),2.)
    455                           + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
     424    xsection = G4double(Zt)*(35.45 + B*std::pow(std::log(sMand/s0),2.)
     425                + 42.53*std::pow(sMand,-eta1) - 33.34*std::pow(sMand,-eta2));
     426
     427    xsection += Nt*(35.80 + B*std::pow(std::log(sMand/s0),2.)
     428                  + 40.15*std::pow(sMand,-eta1) - 30.*std::pow(sMand,-eta2));
    456429  }
    457430  xsection *= millibarn; // parametrised in mb
    … …  
    460433
    461434
    462 
    463 
    464 /////////////////////////////////////////////////////////////////////////////////////
     435///////////////////////////////////////////////////////////////////////////////
    465436//
    466437// Returns nucleon-nucleon cross-section based on N. Starkov parametrisation of
    … …  
    469440
    470441G4double
    471 G4GGNuclNuclCrossSection::GetHadronNucleonXscNS( G4ParticleDefinition* pParticle,
     442G4GGNuclNuclCrossSection::GetHadronNucleonXscNS(G4ParticleDefinition* pParticle,
    472443                                                 G4double pTkin,
    473444                                                 G4ParticleDefinition* tParticle)
    … …  
    477448  G4double hnXsc(0);
    478449
    479 
    480   G4double targ_mass     = tParticle->GetPDGMass();
    481   G4double proj_mass     = pParticle->GetPDGMass();
     450  G4double targ_mass = tParticle->GetPDGMass();
     451  G4double proj_mass = pParticle->GetPDGMass();
    482452
    483453  G4double proj_energy   = proj_mass + pTkin;
    … …  
    497467  //  G4double eta2 = 0.458;
    498468  //  G4double B    = 0.308;
    499 
    500  
    501 
    502 
    503469 
    504470  if( proj_momentum >= 10. ) // high energy: pp = nn = np
    … …  
    506472  {
    507473    Delta = 1.;
    508 
    509     if( proj_energy < 40. ) Delta = 0.916+0.0021*proj_energy;
    510 
    511     if( proj_momentum >= 10.)
    512     {
    513         B0 = 7.5;
    514         A0 = 100. - B0*std::log(3.0e7);
    515 
    516         xsection = A0 + B0*std::log(proj_energy) - 11
    517                   + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
    518                      0.93827*0.93827,-0.165);        //  mb
     474    if (proj_energy < 40.) Delta = 0.916+0.0021*proj_energy;
     475
     476    if (proj_momentum >= 10.) {
     477      B0 = 7.5;
     478      A0 = 100. - B0*std::log(3.0e7);
     479
     480      xsection = A0 + B0*std::log(proj_energy) - 11
     481                + 103*std::pow(2*0.93827*proj_energy + proj_mass*proj_mass+
     482                  0.93827*0.93827,-0.165);        //  mb
    519483    }
    520484  }
    … …  
    562526}
    563527
    564 /*
    565 /////////////////////////////////////////////////////////////////////////////////////
    566 //
    567 // Returns hadron-nucleon inelastic cross-section based on proper parametrisation
    568 
    569 G4double
    570 G4GGNuclNuclCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
    571                                                   const G4Element* anElement          )
    572 {
    573   G4double At = anElement->GetN();  // number of nucleons
    574   G4double Zt = anElement->GetZ();  // number of protons
    575 
    576 
    577   return GetHNinelasticXsc( aParticle, At, Zt );
    578 }
    579 
    580 /////////////////////////////////////////////////////////////////////////////////////
     528/////////////////////////////////////////////////////////////////////////////////
    581529//
    582530// Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
    583531
    584532G4double
    585 G4GGNuclNuclCrossSection::GetHNinelasticXsc(const G4DynamicParticle* aParticle,
    586                                                      G4double At,  G4double Zt )
    587 {
    588   //  G4ParticleDefinition* hadron = aParticle->GetDefinition();
    589   G4double sumInelastic, Nt = At - Zt;
    590 
    591   if(Nt < 0.) Nt = 0.;
    592  
    593   sumInelastic  = Zt*GetHadronNucleonXscNS(aParticle, theProton);
    594   sumInelastic += Nt*GetHadronNucleonXscNS(aParticle, theNeutron);   
    595  
    596   return sumInelastic;
    597 }
    598 */
    599 
    600 /////////////////////////////////////////////////////////////////////////////////////
    601 //
    602 // Returns hadron-nucleon inelastic cross-section based on FTF-parametrisation
    603 
    604 G4double
    605533G4GGNuclNuclCrossSection::GetHNinelasticXscVU(const G4DynamicParticle* aParticle,
    606                                                      G4double At,  G4double Zt )
    607 {
    608   G4int PDGcode    = aParticle->GetDefinition()->GetPDGEncoding();
     534                                              G4int At, G4int Zt)
     535{
     536  G4int PDGcode = aParticle->GetDefinition()->GetPDGEncoding();
    609537  G4int absPDGcode = std::abs(PDGcode);
    610 
    611538  G4double Elab = aParticle->GetTotalEnergy();             
    612539                          // (s - 2*0.88*GeV*GeV)/(2*0.939*GeV)/GeV;
    … …  
    650577                           NumberOfTargetNeutrons * XelPN   );
    651578  }
     579
    652580  Xinelastic = Xtotal - Xelastic;
    653 
    654581  if(Xinelastic < 0.) Xinelastic = 0.;
    655582
    … …  
    657584}
    658585
    659 ////////////////////////////////////////////////////////////////////////////////////
    660 //
    661 //
    662 
    663 G4double
    664 G4GGNuclNuclCrossSection::GetNucleusRadius( const G4DynamicParticle* ,
    665                                                const G4Element* anElement)
    666 {
    667   G4double At       = anElement->GetN();
     586///////////////////////////////////////////////////////////////////////////////
     587//
     588//
     589
     590G4double
     591G4GGNuclNuclCrossSection::GetNucleusRadius(const G4DynamicParticle* ,
     592                                           const G4Element* anElement)
     593{
     594  G4double At = anElement->GetN();
    668595  G4double oneThird = 1.0/3.0;
    669596  G4double cubicrAt = std::pow (At, oneThird);
    670597
    671 
    672598  G4double R;  // = fRadiusConst*cubicrAt;
    673   /* 
    674   G4double tmp = std::pow( cubicrAt-1., 3.);
    675   tmp         += At;
    676   tmp         *= 0.5;
    677 
    678   if (At > 20.)   // 20.
    679   {
    680     R = fRadiusConst*std::pow (tmp, oneThird);
    681   }
    682   else
    683   {
    684     R = fRadiusConst*cubicrAt;
    685   }
    686   */
    687  
    688599  R = fRadiusConst*cubicrAt;
    689600
    690   // return R;  // !!!!
    691 
    692 
    693  
    694601  G4double meanA  = 21.;
    695 
    696602  G4double tauA1  = 40.;
    697603  G4double tauA2  = 10.;
    … …  
    715621  {
    716622    R *= ( 1.0 + b3*( 1. - std::exp( (At - meanA)/tauA3) ) );
    717   } 
     623  }
     624
    718625  return R;
    719  
    720 }
    721 
    722 ////////////////////////////////////////////////////////////////////////////////////
     626}
     627
     628///////////////////////////////////////////////////////////////////////////////
    723629//
    724630//
    … …  
    728634{
    729635  G4double R;
    730 
    731   // R = GetNucleusRadiusGG(At);
    732 
    733636  R = GetNucleusRadiusDE(At);
    734637
    … …  
    741644G4GGNuclNuclCrossSection::GetNucleusRadiusGG(G4double At)
    742645{
    743  
    744646  G4double oneThird = 1.0/3.0;
    745647  G4double cubicrAt = std::pow (At, oneThird);
    746648
    747 
    748   G4double R;  // = fRadiusConst*cubicrAt;
    749  
    750   /*
    751   G4double tmp = std::pow( cubicrAt-1., 3.);
    752   tmp         += At;
    753   tmp         *= 0.5;
    754 
    755   if (At > 20.)
    756   {
    757     R = fRadiusConst*std::pow (tmp, oneThird);
     649  G4double R;  // = fRadiusConst*cubicrAt; 
     650  R = fRadiusConst*cubicrAt;
     651
     652  G4double meanA = 20.;
     653  G4double tauA  = 20.;
     654
     655  if ( At > 20.)   // 20.
     656  {
     657    R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) );
    758658  }
    759659  else
    760660  {
    761     R = fRadiusConst*cubicrAt;
    762   }
    763   */
    764  
    765   R = fRadiusConst*cubicrAt;
    766 
    767   G4double meanA = 20.;
    768   G4double tauA  = 20.;
    769 
    770   if ( At > 20.)   // 20.
    771   {
    772     R *= ( 0.8 + 0.2*std::exp( -(At - meanA)/tauA) );
    773   }
    774   else
    775   {
    776661    R *= ( 1.0 + 0.1*( 1. - std::exp( (At - meanA)/tauA) ) );
    777662  }
    778663
    779664  return R;
    780 
    781665}
    782666
    … …  
    785669G4GGNuclNuclCrossSection::GetNucleusRadiusDE(G4double A)
    786670{
    787 
    788671  // algorithm from diffuse-elastic
    789672
    … …  
    797680
    798681
    799   if( A < 50. )
     682  if (A < 50.)
    800683  {
    801684    if( 10  < A && A <= 15. )     r0  = a11*( 1 - std::pow(A, -2./3.) )*fermi;   // 1.08*fermi;
    … …  
    816699
    817700
    818 ////////////////////////////////////////////////////////////////////////////////////
    819 //
    820 //
    821 
    822 G4double G4GGNuclNuclCrossSection::CalculateEcmValue( const G4double mp ,
    823                                                          const G4double mt ,
    824                                                          const G4double Plab )
     701///////////////////////////////////////////////////////////////////////////////
     702//
     703//
     704
     705G4double G4GGNuclNuclCrossSection::CalculateEcmValue(const G4double mp,
     706                                                     const G4double mt,
     707                                                     const G4double Plab)
    825708{
    826709  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
    … …  
    833716
    834717
    835 ////////////////////////////////////////////////////////////////////////////////////
    836 //
    837 //
    838 
    839 G4double G4GGNuclNuclCrossSection::CalcMandelstamS( const G4double mp ,
    840                                                        const G4double mt ,
    841                                                        const G4double Plab )
     718///////////////////////////////////////////////////////////////////////////////
     719//
     720//
     721
     722G4double G4GGNuclNuclCrossSection::CalcMandelstamS(const G4double mp,
     723                                                   const G4double mt,
     724                                                   const G4double Plab)
    842725{
    843726  G4double Elab = std::sqrt ( mp * mp + Plab * Plab );
    … …  
    847730}
    848731
    849 
    850 //
    851 //
    852 ///////////////////////////////////////////////////////////////////////////////////////
     732//
     733//
     734///////////////////////////////////////////////////////////////////////////////
Note: See TracChangeset for help on using the changeset viewer.