Ignore:
Timestamp:
May 28, 2009, 4:26:57 PM (15 years ago)
Author:
garnier
Message:

maj sur la beta de geant 4.9.3

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r1007 r1055  
    2424// ********************************************************************
    2525//
    26 // $Id: G4WentzelVIModel.cc,v 1.16 2008/11/19 11:47:50 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-02 $
     26// $Id: G4WentzelVIModel.cc,v 1.32 2009/05/10 16:09:29 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    5252//
    5353
    54 
    5554//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    5655//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    6059#include "G4LossTableManager.hh"
    6160#include "G4ParticleChangeForMSC.hh"
    62 #include "G4TransportationManager.hh"
    63 #include "G4SafetyHelper.hh"
    6461#include "G4PhysicsTableHelper.hh"
    6562#include "G4ElementVector.hh"
    … …  
    7168
    7269//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
     70
     71G4double G4WentzelVIModel::ScreenRSquare[] = {0.0};
     72G4double G4WentzelVIModel::FormFactor[]    = {0.0};
    7373
    7474using namespace std;
    … …  
    9696  thePositron = G4Positron::Positron();
    9797  theProton   = G4Proton::Proton();
    98   a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885);
     98  lowEnergyLimit = 0.1*keV;
    9999  G4double p0 = electron_mass_c2*classic_electr_radius;
    100100  coeff  = twopi*p0*p0;
    101   constn = 6.937e-6/(MeV*MeV);
    102101  tkin = targetZ = mom2 = DBL_MIN;
    103102  ecut = etag = DBL_MAX;
    … …  
    106105  xsecn.resize(nelments);
    107106  prob.resize(nelments);
    108   for(size_t j=0; j<100; j++) {
    109     FF[j]    = 0.0;
    110   }
     107
     108  // Thomas-Fermi screening radii
     109  // Formfactors from A.V. Butkevich et al., NIM A 488 (2002) 282
     110
     111  if(0.0 == ScreenRSquare[0]) {
     112    G4double a0 = electron_mass_c2/0.88534;
     113    G4double constn = 6.937e-6/(MeV*MeV);
     114
     115    ScreenRSquare[0] = alpha2*a0*a0;
     116    for(G4int j=1; j<100; j++) {
     117      G4double x = a0*fNistManager->GetZ13(j);
     118      ScreenRSquare[j] = alpha2*x*x;
     119      x = fNistManager->GetA27(j);
     120      FormFactor[j] = constn*x*x;
     121    }
     122  }
    111123}
    112124
    … …  
    132144  if(!isInitialized) {
    133145    isInitialized = true;
    134 
    135     if (pParticleChange)
    136       fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange);
    137     else
    138       fParticleChange = new G4ParticleChangeForMSC();
    139 
    140     safetyHelper = G4TransportationManager::GetTransportationManager()
    141       ->GetSafetyHelper();
    142     safetyHelper->InitialiseHelper();
     146    fParticleChange = GetParticleChangeForMSC();
     147    InitialiseSafetyHelper();
    143148  }
    144149}
    … …  
    156161  SetupKinematic(ekin, cutEnergy);
    157162  SetupTarget(Z, ekin);
    158   G4double xsec = ComputeTransportXSectionPerVolume();
     163  G4double xsec = ComputeTransportXSectionPerAtom();
    159164  /* 
    160165  G4cout << "CS: e= " << tkin << " cosEl= " << cosTetMaxElec2
    … …  
    167172//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    168173
    169 G4double G4WentzelVIModel::ComputeTransportXSectionPerVolume()
     174G4double G4WentzelVIModel::ComputeTransportXSectionPerAtom()
    170175{
    171176  G4double xSection = 0.0;
    … …  
    189194      y = 0.0;
    190195    }
    191     xSection += y/targetZ;
    192   }
    193   /*
    194   G4cout << "G4WentzelVIModel:XS per A " << " Z= " << Z << " e(MeV)= " << kinEnergy/MeV
     196    xSection = y;
     197  }
     198  /* 
     199  G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ
     200         << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
    195201         << " cut(MeV)= " << ecut/MeV 
    196202         << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
    197          << " zmaxN= " << (1.0 - cosTetMsxNuc)/screenZ << G4endl;
     203         << " zmaxN= " << (1.0 - cosTetMaxNuc2)/screenZ
     204         << " costm= " << cosTetMaxNuc2 << G4endl;
    198205  */
    199 
    200206  // scattering off nucleus
    201207  if(cosTetMaxNuc2 < 1.0) {
    202208    x  = 1.0 - cosTetMaxNuc2;
    203209    x1 = screenZ*formfactA;
    204     x2 = 1.0/(1.0 - x1);
     210    x2 = 1.0 - x1;
    205211    x3 = x/screenZ;
    206212    x4 = formfactA*x;
    207213    // low-energy limit
    208214    if(x3 < numlimit && x1 < numlimit) {
    209       y = 0.5*x3*x3*x2*x2*x2*(1.0 - 1.333333*x3 + 1.5*x3*x3 - 1.5*x1
    210                               + 3.0*x1*x1 + 2.666666*x3*x1);
     215      y = 0.5*x3*x3*(1.0 - 1.3333333*x3 + 1.5*x3*x3 - 1.5*x1
     216                     + 3.0*x1*x1 + 2.666666*x3*x1)/(x2*x2*x2);
    211217      // high energy limit
    212     } else if(1.0 < x1) {
     218    } else if(x2 <= 0.0) {
    213219      x4 = x1*(1.0 + x3);
    214220      y  = x3*(1.0 + 0.5*x3 - (2.0 - x1)*(1.0 + x3 + x3*x3/3.0)/x4)/(x4*x4);
    … …  
    216222    } else {
    217223      y = ((1.0 + x1)*x2*log((1. + x3)/(1. + x4))
    218            - x3/(1. + x3) - x4/(1. + x4))*x2*x2;
    219     }
     224           - x3/(1. + x3) - x4/(1. + x4))/(x2*x2);
     225    }
     226    //G4cout << "y= " << y << " x1= " <<x1<<"  x2= " <<x2
     227    //     <<"  x3= "<<x3<<"  x4= " << x4<<G4endl;
    220228    if(y < 0.0) {
    221229      nwarnings++;
    … …  
    230238      y = 0.0;
    231239    }
    232     xSection += y;
    233   }
    234   xSection *= (coeff*targetZ*targetZ*chargeSquare*invbeta2/mom2);
    235   //  G4cout << "   XStotal= " << xSection/barn << " screenZ= " << screenZ
    236   //     << " formF= " << formfactA << " for " << p->GetParticleName() << G4endl;
     240    xSection += y*targetZ;
     241  }
     242  xSection *= kinFactor;
     243  /*
     244  G4cout << "Z= " << targetZ << " XStot= " << xSection/barn
     245         << " screenZ= " << screenZ << " formF= " << formfactA
     246         << " for " << particle->GetParticleName()
     247         << " m= " << mass << " 1/v= " << sqrt(invbeta2) << " p= " << sqrt(mom2)
     248         << " x= " << x
     249         << G4endl;
     250  */
    237251  return xSection;
    238252}
    … …  
    277291  // i.e. when it is needed for optimization purposes
    278292  if(stepStatus != fGeomBoundary && presafety < tlimitminfix)
    279     presafety = safetyHelper->ComputeSafety(sp->GetPosition());
     293    presafety = ComputeSafety(sp->GetPosition(), tlimit);
    280294  /*
    281295  G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit tlimit= "
    … …  
    291305    G4double rlimit = facrange*lambda0;
    292306    G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
    293     if(rcut > rlimit) rlimit = std::pow(2.0*rcut*rcut*lambda0,0.33333333);
     307    if(rcut > rlimit) rlimit = std::pow(rcut*rcut*rlimit,0.33333333);
     308    rlimit = std::min(rlimit, facgeom*currentMaterial->GetRadlen());
    294309    if(rlimit < tlimit) tlimit = rlimit;
    295310  }
    … …  
    404419  } else {
    405420
    406     // define threshold angle as 2 sigma of central value
     421    // define threshold angle between single and multiple scattering
    407422    cosThetaMin = 1.0 - 3.0*x1;
    408423
    … …  
    570585
    571586    if(r > tlimitminfix) {
    572       G4ThreeVector Position = *(fParticleChange->GetProposedPosition());
    573       G4double fac= 1.;
    574       if(r >= safety) {
    575         //  ******* so safety is computed at boundary too ************
    576         G4double newsafety =
    577           safetyHelper->ComputeSafety(Position) - tlimitminfix;
    578         if(newsafety <= 0.0) fac = 0.0;
    579         else if(r > newsafety) fac = newsafety/r ;
    580         //G4cout << "NewSafety= " << newsafety << " fac= " << fac
    581         // << " r= " << r << " sint= " << sint << " pos " << Position << G4endl;
    582       } 
    583 
    584       if(fac > 0.) {
    585         // compute new endpoint of the Step
    586         G4ThreeVector newPosition = Position + fac*pos;
    587 
    588         // check safety after displacement
    589         G4double postsafety = safetyHelper->ComputeSafety(newPosition);
    590 
    591         // displacement to boundary
    592         if(postsafety <= 0.0) {
    593           safetyHelper->Locate(newPosition, newDirection);
    594 
    595           // not on the boundary
    596         } else {
    597           safetyHelper->ReLocateWithinVolume(newPosition);
    598           // if(fac < 1.0) G4cout << "NewPosition " << newPosition << G4endl;
    599         }
    600      
    601         fParticleChange->ProposePosition(newPosition);
    602       }
     587      pos /= r;
     588      ComputeDisplacement(fParticleChange, pos, r, safety);
    603589    }
    604590  }
    … …  
    624610  G4double xs = 0.0;
    625611
    626   G4double fac = coeff*chargeSquare*invbeta2/mom2;
    627 
    628612  for (G4int i=0; i<nelm; i++) {
    629613    SetupTarget((*theElementVector)[i]->GetZ(), tkin);
    … …  
    635619    cosTetMaxNuc2  = std::max(cosnm,cosThetaMin);
    636620    cosTetMaxElec2 = std::max(cosem,cosThetaMin);
    637     xtsec += ComputeTransportXSectionPerVolume()*density;
     621    xtsec += ComputeTransportXSectionPerAtom()*density;
    638622    // return limit back
    639623    cosTetMaxElec2 = cosem;
    … …  
    643627    G4double nsec = 0.0;
    644628    G4double x1 = 1.0 - cosThetaMin + screenZ;
    645     G4double f  = fac*targetZ*density;
     629    G4double f  = kinFactor*density;
    646630
    647631    // scattering off electrons
    … …  
    656640      G4double s  = screenZ*formfactA;
    657641      G4double z1 = 1.0 - cosnm + screenZ;
    658       G4double d  = (1.0 - s)/formfactA;
     642      G4double s1 = 1.0 - s;
     643      G4double d  = s1/formfactA;
    659644
    660645      // check numerical limit
    … …  
    667652        G4double x2 = x1 + d;
    668653        G4double z2 = z1 + d;
    669         nsec = (1.0 + 2.0*s)*((cosThetaMin - cosnm)*(1.0/(x1*z1) + 1.0/(x2*z2)) -
    670                            2.0*log(z1*x2/(z2*x1))/d);
     654        nsec = (1.0/x1 - 1.0/z1 + 1.0/x2 - 1.0/z2 - 2.0*log(z1*x2/(z2*x1))/d)/(s1*s1);
    671655      }
    672656      nsec *= f*targetZ;
    … …  
    801785      G4double mom22 = t1*(t1 + 2.0*mass);
    802786      G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
    803       if(ctm < 1.0) cosTetMaxElec = ctm;
     787      if(ctm <  1.0) cosTetMaxElec = ctm;
     788      if(ctm < -1.0) cosTetMaxElec = -1.0;
    804789    }
    805790  }
    … …  
    808793
    809794//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    810 
    811 void G4WentzelVIModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
    812                                      const G4MaterialCutsCouple*,
    813                                      const G4DynamicParticle*,
    814                                      G4double,
    815                                      G4double)
    816 {}
    817 
    818 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracChangeset for help on using the changeset viewer.