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/G4WentzelVIModel.cc

    r1228 r1315  
    2424// ********************************************************************
    2525//
    26 // $Id: G4WentzelVIModel.cc,v 1.37 2009/10/28 10:14:13 vnivanch Exp $
    27 // GEANT4 tag $Name: geant4-09-03 $
     26// $Id: G4WentzelVIModel.cc,v 1.60 2010/06/01 11:13:31 vnivanch Exp $
     27// GEANT4 tag $Name: geant4-09-04-beta-cand-01 $
    2828//
    2929// -------------------------------------------------------------------
    … …  
    3939//
    4040// Modifications:
     41// 27-05-2010 V.Ivanchenko added G4WentzelOKandVIxSection class to
     42//              compute cross sections and sample scattering angle
    4143//
    4244//
    … …  
    5759#include "G4WentzelVIModel.hh"
    5860#include "Randomize.hh"
    59 #include "G4LossTableManager.hh"
    6061#include "G4ParticleChangeForMSC.hh"
    6162#include "G4PhysicsTableHelper.hh"
    6263#include "G4ElementVector.hh"
    6364#include "G4ProductionCutsTable.hh"
    64 #include "G4PhysicsLogVector.hh"
    65 #include "G4Electron.hh"
    66 #include "G4Positron.hh"
    67 #include "G4Proton.hh"
    68 
    69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    70 
    71 G4double G4WentzelVIModel::ScreenRSquare[] = {0.0};
    72 G4double G4WentzelVIModel::FormFactor[]    = {0.0};
     65#include "G4LossTableManager.hh"
     66#include "G4Pow.hh"
     67
     68//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    7369
    7470using namespace std;
    … …  
    7773  G4VMscModel(nam),
    7874  theLambdaTable(0),
    79   theLambda2Table(0),
    80   numlimit(0.2),
    81   nbins(60),
    82   nwarnings(0),
    83   nwarnlimit(50),
     75  numlimit(0.1),
    8476  currentCouple(0),
    8577  cosThetaMin(1.0),
    86   q2Limit(TeV*TeV),
    87   alpha2(fine_structure_const*fine_structure_const),
    8878  isInitialized(false),
    8979  inside(false)
    … …  
    9181  invsqrt12 = 1./sqrt(12.);
    9282  tlimitminfix = 1.e-6*mm;
    93   theManager = G4LossTableManager::Instance();
    94   fNistManager = G4NistManager::Instance();
    95   theElectron = G4Electron::Electron();
    96   thePositron = G4Positron::Positron();
    97   theProton   = G4Proton::Proton();
    98   lowEnergyLimit = 0.1*keV;
    99   G4double p0 = electron_mass_c2*classic_electr_radius;
    100   coeff  = twopi*p0*p0;
    101   tkin = targetZ = mom2 = DBL_MIN;
    102   ecut = etag = DBL_MAX;
     83  lowEnergyLimit = 1.0*eV;
    10384  particle = 0;
    10485  nelments = 5;
    10586  xsecn.resize(nelments);
    10687  prob.resize(nelments);
    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   }
     88  theManager = G4LossTableManager::Instance();
     89  fG4pow = G4Pow::GetInstance();
     90  wokvi = new G4WentzelOKandVIxSection();
    12391}
    12492
    … …  
    12694
    12795G4WentzelVIModel::~G4WentzelVIModel()
    128 {}
     96{
     97  delete wokvi;
     98}
    12999
    130100//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    … …  
    135105  // reset parameters
    136106  SetupParticle(p);
    137   tkin = targetZ = mom2 = 0.0;
    138   ecut = etag = DBL_MAX;
    139107  currentRange = 0.0;
    140108  cosThetaMax = cos(PolarAngleLimit());
     109  wokvi->Initialise(p, cosThetaMax);
     110  /*
     111  G4cout << "G4WentzelVIModel: factorA2(GeV^2) = " << factorA2/(GeV*GeV)
     112         << "  1-cos(ThetaLimit)= " << 1 - cosThetaMax
     113         << G4endl;
     114  */
    141115  currentCuts = &cuts;
    142116
    … …  
    157131                             G4double cutEnergy, G4double)
    158132{
    159   SetupParticle(p);
    160   if(kinEnergy < lowEnergyLimit) return 0.0;
    161   SetupKinematic(kinEnergy, cutEnergy);
    162   SetupTarget(Z, kinEnergy);
    163   G4double xsec = ComputeTransportXSectionPerAtom();
    164   /*   
    165   G4cout << "CS: e= " << tkin << " cosEl= " << cosTetMaxElec2
    166          << " cosN= " << cosTetMaxNuc2 << " xsec(bn)= " << xsec/barn
    167          << " " << particle->GetParticleName() << G4endl;
     133  G4double xsec = 0.0;
     134  if(p != particle) { SetupParticle(p); }
     135  if(kinEnergy < lowEnergyLimit) { return xsec; }
     136  DefineMaterial(CurrentCouple());
     137  cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
     138  if(cosTetMaxNuc < 1.0) {
     139    cosTetMaxNuc = wokvi->SetupTarget(G4int(Z), cutEnergy);
     140    xsec = wokvi->ComputeTransportCrossSectionPerAtom(cosTetMaxNuc);
     141  /*     
     142    G4cout << "G4WentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= " << kinEnergy
     143           << " 1-cosN= " << 1 - costm << " xsec(bn)= " << xsec/barn
     144           << " " << particle->GetParticleName() << G4endl;
    168145  */
     146  }
    169147  return xsec;
    170 }
    171 
    172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    173 
    174 G4double G4WentzelVIModel::ComputeTransportXSectionPerAtom()
    175 {
    176   G4double xSection = 0.0;
    177   G4double x, y, x1, x2, x3, x4;
    178 
    179   // scattering off electrons
    180   if(cosTetMaxElec2 < 1.0) {
    181     x = (1.0 - cosTetMaxElec2)/screenZ;
    182     if(x < numlimit) y = 0.5*x*x*(1.0 - 1.3333333*x + 1.5*x*x);
    183     else             y = log(1.0 + x) - x/(1.0 + x);
    184     if(y < 0.0) {
    185       nwarnings++;
    186       if(nwarnings < nwarnlimit /*&& y < -1.e-10*/) {
    187         G4cout << "Electron scattering <0 for L1 " << y
    188                << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2)
    189                << " Z= " << targetZ << "  "
    190                << particle->GetParticleName() << G4endl;
    191         G4cout << " z= " << 1.0-cosTetMaxElec2 << " screenZ= " << screenZ
    192                << " x= " << x << G4endl;
    193       }
    194       y = 0.0;
    195     }
    196     xSection = y;
    197   }
    198   /* 
    199   G4cout << "G4WentzelVI:XS per A " << " Z= " << targetZ
    200          << " e(MeV)= " << tkin/MeV << " XSel= " << xSection
    201          << " cut(MeV)= " << ecut/MeV 
    202          << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
    203          << " zmaxN= " << (1.0 - cosTetMaxNuc2)/screenZ
    204          << " costm= " << cosTetMaxNuc2 << G4endl;
    205   */
    206   // scattering off nucleus
    207   if(cosTetMaxNuc2 < 1.0) {
    208     x  = 1.0 - cosTetMaxNuc2;
    209     x1 = screenZ*formfactA;
    210     x2 = 1.0 - x1;
    211     x3 = x/screenZ;
    212     x4 = formfactA*x;
    213     // low-energy limit
    214     if(x3 < numlimit && x1 < numlimit) {
    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);
    217       // high energy limit
    218     } else if(x2 <= 0.0) {
    219       x4 = x1*(1.0 + x3);
    220       y  = x3*(1.0 + 0.5*x3 - (2.0 - x1)*(1.0 + x3 + x3*x3/3.0)/x4)/(x4*x4);
    221       // middle energy
    222     } else {
    223       y = ((1.0 + x1)*x2*log((1. + x3)/(1. + x4))
    224            - x3/(1. + x3) - x4/(1. + x4))/(x2*x2);
    225     }
    226     //G4cout << "y= " << y << " x1= " <<x1<<"  x2= " <<x2
    227     //     <<"  x3= "<<x3<<"  x4= " << x4<<G4endl;
    228     if(y < 0.0) {
    229       nwarnings++;
    230       if(nwarnings < nwarnlimit /*&& y < -1.e-10*/) {
    231         G4cout << "Nuclear scattering <0 for L1 " << y
    232                << " e(MeV)= " << tkin << " Z= " << targetZ << "  "
    233                << particle->GetParticleName() << G4endl;
    234         G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
    235                << " x= " << " x1= " << x1 << " x2= " << x2
    236                << " x3= " << x3 << " x4= " << x4 <<G4endl;
    237       }
    238       y = 0.0;
    239     }
    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   */
    251   return xSection;
    252148}
    253149
    … …  
    263159  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
    264160  G4StepStatus stepStatus = sp->GetStepStatus();
     161  //G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit stepStatus= "
     162  //     << stepStatus << G4endl;
    265163
    266164  // initialisation for 1st step 
    … …  
    268166    inside = false;
    269167    SetupParticle(dp->GetDefinition());
    270     theLambdaTable = theTable;
    271168  }
    272169
    … …  
    274171  preKinEnergy  = dp->GetKineticEnergy();
    275172  DefineMaterial(track.GetMaterialCutsCouple());
    276   lambda0 = GetLambda(preKinEnergy);
     173  theLambdaTable = theTable;
     174  lambdaeff = GetLambda(preKinEnergy);
    277175  currentRange =
    278176    theManager->GetRangeFromRestricteDEDX(particle,preKinEnergy,currentCouple);
     177  cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial);
    279178
    280179  // extra check for abnormal situation
    281180  // this check needed to run MSC with eIoni and eBrem inactivated
    282   if(tlimit > currentRange) tlimit = currentRange;
     181  if(tlimit > currentRange) { tlimit = currentRange; }
    283182
    284183  // stop here if small range particle
    285   if(inside) return tlimit;   
     184  if(inside) { return tlimit; }
    286185
    287186  // pre step
    … …  
    290189  // compute presafety again if presafety <= 0 and no boundary
    291190  // i.e. when it is needed for optimization purposes
    292   if(stepStatus != fGeomBoundary && presafety < tlimitminfix)
     191  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
    293192    presafety = ComputeSafety(sp->GetPosition(), tlimit);
     193  }
    294194  /*
    295   G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit tlimit= "
    296          <<tlimit<<" safety= " << presafety
    297          << " range= " <<currentRange<<G4endl;
     195  G4cout << "e(MeV)= " << preKinEnergy/MeV
     196         << "  " << particle->GetParticleName()
     197         << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
     198         << " R(mm)= " <<currentRange/mm
     199         << " L0(mm^-1)= " << lambdaeff*mm
     200         <<G4endl;
    298201  */
    299202  // far from geometry boundary
    300203  if(currentRange < presafety) {
    301204    inside = true; 
    302    
    303     // limit mean scattering angle
    304   } else {
    305     G4double rlimit = facrange*lambda0;
    306     G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
    307     if(rcut > rlimit) rlimit = std::pow(rcut*rcut*rlimit,0.33333333);
    308     rlimit = std::min(rlimit, facgeom*currentMaterial->GetRadlen());
    309     if(rlimit < tlimit) tlimit = rlimit;
    310   }
    311   /*
     205    return tlimit;
     206  }
     207
     208  // natural limit for high energy
     209  G4double rlimit = std::max(facrange*currentRange,
     210                             0.7*(1.0 - cosTetMaxNuc)*lambdaeff);
     211
     212  // low-energy e-
     213  if(cosThetaMax > cosTetMaxNuc) {
     214    rlimit = std::min(rlimit, facsafety*presafety);
     215  }
     216   
     217  // cut correction
     218  G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
     219  //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= " << presafety
     220  // << " 1-cosThetaMax= " <<1-cosThetaMax << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
     221  // << G4endl;
     222  if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
     223
     224  if(rlimit < tlimit) { tlimit = rlimit; }
     225
     226  tlimit = std::max(tlimit, tlimitminfix);
     227
     228  // step limit in infinite media
     229  tlimit = std::min(tlimit, 20*currentMaterial->GetRadlen());
     230  /* 
    312231  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
    313          << " L0= " << lambda0 << " R= " << currentRange
     232         << " L0= " << lambdaeff << " R= " << currentRange
    314233         << "tlimit= " << tlimit 
    315234         << " currentMinimalStep= " << currentMinimalStep << G4endl;
    … …  
    324243  tPathLength  = truelength;
    325244  zPathLength  = tPathLength;
    326   lambdaeff    = lambda0;
    327 
    328   if(lambda0 > 0.0) {
    329     G4double tau = tPathLength/lambda0;
     245
     246  if(lambdaeff > 0.0) {
     247    G4double tau = tPathLength/lambdaeff;
    330248    //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
    331     //   << " lambda0= " << lambda0 << " tau= " << tau << G4endl;
     249    //   << " Leff= " << lambdaeff << " tau= " << tau << G4endl;
    332250    // small step
    333251    if(tau < numlimit) {
    … …  
    336254      // medium step
    337255    } else {
    338       //      zPathLength = lambda0*(1.0 - exp(-tPathLength/lambda0));
    339256      G4double e1 = 0.0;
    340257      if(currentRange > tPathLength) {
    … …  
    343260                                   currentCouple);
    344261      }
    345       lambdaeff = GetLambda(0.5*(e1 + preKinEnergy));
     262      e1 = 0.5*(e1 + preKinEnergy);
     263      cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
     264      lambdaeff = GetLambda(e1);
    346265      zPathLength = lambdaeff*(1.0 - exp(-tPathLength/lambdaeff));
    347266    }
    … …  
    355274G4double G4WentzelVIModel::ComputeTrueStepLength(G4double geomStepLength)
    356275{
    357   // step defined other than transportation
    358   if(geomStepLength == zPathLength) return tPathLength;
     276  // initialisation of single scattering x-section
     277  xtsec = 0.0;
     278
     279  // pathalogical case
     280  if(lambdaeff <= 0.0) {
     281    zPathLength = geomStepLength;
     282    tPathLength = geomStepLength;
     283    return tPathLength;
     284  }
     285
     286  G4double tau = geomStepLength/lambdaeff;
    359287
    360288  // step defined by transportation
    361   tPathLength  = geomStepLength;
    362   zPathLength  = geomStepLength;
    363   G4double tau = zPathLength/lambdaeff;
    364   tPathLength *= (1.0 + 0.5*tau + tau*tau/3.0);
    365 
    366   if(tau > numlimit) {
    367     G4double e1 = 0.0;
    368     if(currentRange > tPathLength) {
    369       e1 = theManager->GetEnergy(particle,
    370                                  currentRange-tPathLength,
    371                                  currentCouple);
    372     }
    373     lambdaeff = GetLambda(0.5*(e1 + preKinEnergy));
    374     tau = zPathLength/lambdaeff;
    375 
    376     if(tau < 0.999999) tPathLength = -lambdaeff*log(1.0 - tau);
    377     else               tPathLength = currentRange;
    378 
    379     if(tPathLength < zPathLength) tPathLength = zPathLength;
    380   }
    381   if(tPathLength > currentRange) tPathLength = currentRange;
    382   //G4cout<<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
     289  if(geomStepLength != zPathLength) {
     290
     291    // step defined by transportation
     292    zPathLength = geomStepLength;
     293    tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
     294
     295    // energy correction for a big step
     296    if(tau > numlimit) {
     297      G4double e1 = 0.0;
     298      if(currentRange > tPathLength) {
     299        e1 = theManager->GetEnergy(particle,
     300                                   currentRange-tPathLength,
     301                                   currentCouple);
     302      }
     303      e1 = 0.5*(e1 + preKinEnergy);
     304      cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
     305      lambdaeff = GetLambda(e1);
     306      tau = zPathLength/lambdaeff;
     307     
     308      if(tau < 0.999999) { tPathLength = -lambdaeff*log(1.0 - tau); }
     309      else               { tPathLength = currentRange; }
     310    }
     311  }
     312
     313  // check of step length
     314  // define threshold angle between single and multiple scattering
     315  cosThetaMin = 1.0 - 1.5*tPathLength/lambdaeff;
     316
     317  // recompute transport cross section - do not change energy
     318  // anymore - cannot be applied for big steps
     319  if(cosThetaMin > cosTetMaxNuc) {
     320
     321    // new computation
     322    G4double xsec = ComputeXSectionPerVolume();
     323    //G4cout << "%%%% xsec= " << xsec << "  xtsec= " << xtsec << G4endl;
     324    if(xtsec > 0.0) {
     325      if(xsec > 0.0) { lambdaeff = 1./xsec; }
     326      else           { lambdaeff = DBL_MAX; }
     327
     328      tau = zPathLength*xsec;
     329      if(tau < numlimit) { tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0); }
     330      else if(tau < 0.999999) { tPathLength = -lambdaeff*log(1.0 - tau); }
     331      else               { tPathLength = currentRange; }
     332    }
     333  }
     334
     335  if(tPathLength > currentRange) { tPathLength = currentRange; }
     336  if(tPathLength < zPathLength)  { tPathLength = zPathLength; }
     337  /*   
     338  G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
     339         <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
     340  G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
     341         << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
     342         << " e(MeV)= " << preKinEnergy/MeV << G4endl;
     343  */
    383344  return tPathLength;
    384345}
    … …  
    391352  //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
    392353  //     << particle->GetParticleName() << G4endl;
    393   G4double kinEnergy = dynParticle->GetKineticEnergy();
    394 
    395   // ignore scattering for zero step length and enegy below the limit
    396   if(kinEnergy < lowEnergyLimit || tPathLength <= DBL_MIN) return;
    397 
    398   G4double ekin = preKinEnergy;
    399   if(ekin - kinEnergy > ekin*dtrl) {
    400     ekin = 0.5*(preKinEnergy + kinEnergy);
    401     lambdaeff = GetLambda(ekin);
    402   } 
     354
     355  // ignore scattering for zero step length and energy below the limit
     356  if(dynParticle->GetKineticEnergy() < lowEnergyLimit ||
     357     tPathLength <= DBL_MIN || lambdaeff <= DBL_MIN)
     358    { return; }
    403359 
    404   G4double x1 = 0.5*tPathLength/lambdaeff;
    405   G4double cut= (*currentCuts)[currentMaterialIndex];
    406   /* 
    407   G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy<<" Eeff(MeV)= "<<ekin/MeV
    408          << " L0= " << lambda0 << " Leff= " << lambdaeff
    409          << " x1= " << x1 << " safety= " << safety << G4endl;
     360  G4double invlambda = 0.0;
     361  if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
     362
     363  // use average kinetic energy over the step
     364  G4double cut = (*currentCuts)[currentMaterialIndex];
     365  /*     
     366  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
     367         << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
     368         << " x1= " <<  tPathLength*invlambda << " safety= " << safety << G4endl;
    410369  */
    411370
    412   G4double xsec = 0.0;
    413   G4bool largeAng = false;
    414 
    415   // large scattering angle case
    416   if(x1 > 0.5) {
    417     x1 *= 0.5;
    418     largeAng = true;
    419 
    420     // normal case
    421   } else {
    422 
    423     // define threshold angle between single and multiple scattering
    424     cosThetaMin = 1.0 - 3.0*x1;
    425 
    426     // for low-energy e-,e+ no limit
    427     SetupKinematic(ekin, cut);
    428  
    429     // recompute transport cross section
    430     if(cosThetaMin > cosTetMaxNuc) {
    431 
    432       xsec = ComputeXSectionPerVolume();
    433 
    434       if(xtsec > 0.0) x1 = 0.5*tPathLength*xtsec;
    435       else            x1 = 0.0;
    436 
    437       /*     
    438         G4cout << "cosTetMaxNuc= " << cosTetMaxNuc
    439         << " cosThetaMin= " << cosThetaMin
    440         << " cosThetaMax= " << cosThetaMax
    441         << " cosTetMaxElec2= " << cosTetMaxElec2 << G4endl;
    442         G4cout << "Recomputed xsec(1/mm)= " << xsec << " x1= " << x1 << G4endl;
    443       */
    444     }
    445   }
    446 
    447   // result of central part sampling
    448   G4double z;
     371  G4double length = tPathLength;
     372  G4double lengthlim = tPathLength*1.e-6;
     373
     374  // step limit due msc
     375  G4double x0 = length;
     376  // large scattering angle case - two step approach
     377  if(tPathLength*invlambda > 0.5 && length > tlimitminfix) { x0 *= 0.5; }
     378
     379  // step limit due single scattering
     380  G4double x1 = length;
     381  if(xtsec > 0.0) { x1 = -log(G4UniformRand())/xtsec; }
     382
     383  const G4ElementVector* theElementVector =
     384    currentMaterial->GetElementVector();
     385  G4int nelm = currentMaterial->GetNumberOfElements();
     386
     387  // geometry
     388  G4double sint, cost, phi;
     389  G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
     390  G4ThreeVector temp(0.0,0.0,1.0);
     391
     392  // current position and direction relative to the end point
     393  // because of magnetic field geometry is computed relatively to the
     394  // end point of the step
     395  G4ThreeVector dir(0.0,0.0,1.0);
     396  G4ThreeVector pos(0.0,0.0,-zPathLength);
     397  G4double mscfac = zPathLength/tPathLength;
     398
     399  // start a loop
    449400  do {
    450     z = -x1*log(G4UniformRand());
    451   } while (z > 1.0);
    452 
    453   // cost is sampled ------------------------------
    454   G4double cost = 1.0 - 2.0*z;
    455   //  if(cost < -1.0) cost = -1.0;
    456   // else if(cost > 1.0) cost = 1.0;
    457   G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
    458 
    459   G4double phi  = twopi*G4UniformRand();
    460 
    461   G4double dirx = sint*cos(phi);
    462   G4double diry = sint*sin(phi);
    463  
    464   //G4cout << "G4WentzelVIModel: step(mm)= " << tPathLength/mm
    465   //     << " sint= " << sint << " cost= " << cost<< G4endl;
    466  
    467   G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
    468   G4ThreeVector newDirection(dirx,diry,cost);
    469   G4ThreeVector temp(0.0,0.0,1.0);
    470   G4ThreeVector pos(0.0,0.0,-zPathLength);
    471   G4ThreeVector dir(0.0,0.0,1.0);
    472   G4bool isscat = false;
    473 
    474   // sample MSC scattering for large angle
    475   // extra central scattering for half step
    476   if(largeAng) {
    477     isscat = true;
    478     pos.setZ(-0.5*zPathLength);
    479     do {
    480       z = -x1*log(G4UniformRand());
    481     } while (z > 1.0);
     401    G4double step = x0;
     402    G4bool singleScat = false;
     403
     404    // single scattering case
     405    if(x1 < x0) {
     406      step = x1;
     407      singleScat = true;
     408    }
     409
     410    // new position
     411    pos += step*mscfac*dir;
     412
     413    // added multiple scattering
     414    G4double z;
     415    G4double tet2 = step*invlambda; 
     416    do { z = -tet2*log(G4UniformRand()); } while (z >= 1.0);
     417
    482418    cost = 1.0 - 2.0*z;
    483     //if(std::fabs(cost) > 1.0) cost = 1.0;
    484 
    485419    sint = sqrt((1.0 - cost)*(1.0 + cost));
    486420    phi  = twopi*G4UniformRand();
    487 
    488     // position and direction for secondary scattering
    489     dir.set(sint*cos(phi),sint*sin(phi),cost);
    490     pos += 0.5*dir*zPathLength;
    491     x1 *= 2.0;
    492   }
    493 
    494   // sample Rutherford scattering for large angle
    495   if(xsec > DBL_MIN) {
    496     G4double t = tPathLength;
    497     G4int nelm = currentMaterial->GetNumberOfElements();
    498     const G4ElementVector* theElementVector =
    499       currentMaterial->GetElementVector();
    500     do{
    501       G4double x  = -log(G4UniformRand())/xsec;     
    502       pos += dir*(zPathLength*std::min(x,t)/tPathLength);
    503       t -= x;
    504       if(t > 0.0) {
    505         G4double zz1 = 1.0;
    506         G4double qsec = G4UniformRand()*xsec;
    507 
    508         // scattering off nucleus
    509         G4int i = 0;
    510         if(nelm > 1) {
    511           for (; i<nelm; i++) {if(xsecn[i] >= qsec) break;}
    512           if(i >= nelm) i = nelm - 1;
    513         }
    514         SetupTarget((*theElementVector)[i]->GetZ(), tkin);
    515         G4double formf = formfactA;
    516         G4double costm = cosTetMaxNuc2;
    517         if(prob[i] > 0.0) {
    518           if(G4UniformRand() <= prob[i]) {
    519             formf = 0.0;
    520             costm = cosTetMaxElec2;
    521           }
    522         }
    523         if(cosThetaMin > costm) {
    524 
    525           G4double w1 = 1. - cosThetaMin + screenZ;
    526           G4double w2 = 1. - costm + screenZ;
    527           G4double w3 = cosThetaMin - costm;
    528           G4double grej, zz;
    529           do {
    530             zz = w1*w2/(w1 + G4UniformRand()*w3) - screenZ;
    531             grej = 1.0/(1.0 + formf*zz);
    532           } while ( G4UniformRand() > grej*grej ); 
    533           if(zz < 0.0) zz = 0.0;
    534           else if(zz > 2.0) zz = 2.0;
    535           zz1 = 1.0 - zz;
    536         }
    537         if(zz1 < 1.0) {
    538           isscat = true;
    539           //G4cout << "Rutherford zz1= " << zz1 << " t= " << t << G4endl;
    540           sint = sqrt((1.0 - zz1)*(1.0 + zz1));
    541           //G4cout << "sint= " << sint << G4endl;
    542           phi  = twopi*G4UniformRand();
    543           G4double vx1 = sint*cos(phi);
    544           G4double vy1 = sint*sin(phi);
    545           temp.set(vx1,vy1,zz1);
    546           temp.rotateUz(dir);
    547           dir = temp;
    548         }
     421    G4double vx1 = sint*cos(phi);
     422    G4double vy1 = sint*sin(phi);
     423
     424    // lateral displacement 
     425    if (latDisplasment && safety > tlimitminfix) {
     426      G4double rms = invsqrt12*sqrt(2.0*tet2);
     427      G4double dx = step*(0.5*vx1 + rms*G4RandGauss::shoot(0.0,1.0));
     428      G4double dy = step*(0.5*vy1 + rms*G4RandGauss::shoot(0.0,1.0));
     429      G4double dz;
     430      G4double d = (dx*dx + dy*dy)/(step*step);
     431      if(d < numlimit)  { dz = -0.5*step*d*(1.0 + 0.25*d); }
     432      else if(d < 1.0)  { dz = -step*(1.0 - sqrt(1.0 - d));}
     433      else              { dx = dy = dz = 0.0; }
     434
     435      // change position
     436      temp.set(dx,dy,dz);
     437      temp.rotateUz(dir);
     438      pos += temp;
     439    }
     440
     441    // direction is changed
     442    temp.set(vx1,vy1,cost);
     443    temp.rotateUz(dir);
     444    dir = temp;
     445
     446    if(singleScat) {
     447
     448      // select element
     449      G4int i = 0;
     450      if(nelm > 1) {
     451        G4double qsec = G4UniformRand()*xtsec;
     452        for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
     453        if(i >= nelm) { i = nelm - 1; }
    549454      }
    550     } while (t > 0.0);
    551   }
    552   if(isscat) newDirection.rotateUz(dir);
    553   newDirection.rotateUz(oldDirection);
     455      G4double cosTetM =
     456        wokvi->SetupTarget(G4int((*theElementVector)[i]->GetZ()), cut);
     457      temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
     458      temp.rotateUz(dir);
     459
     460      // renew direction
     461      dir = temp;
     462
     463      // new single scatetring
     464      x1 = -log(G4UniformRand())/xtsec;
     465    }
     466
     467    // update step
     468    length -= step;
     469
     470  } while (length > lengthlim);
     471   
     472  dir.rotateUz(oldDirection);
     473  pos.rotateUz(oldDirection);
    554474
    555475  //G4cout << "G4WentzelVIModel sampling of scattering is done" << G4endl;
    556476  // end of sampling -------------------------------
    557477
    558   fParticleChange->ProposeMomentumDirection(newDirection);
    559 
    560   if (latDisplasment && safety > tlimitminfix) {
    561     G4double rms = invsqrt12*sqrt(2.0*x1);
    562     G4double dx = zPathLength*(0.5*dirx + rms*G4RandGauss::shoot(0.0,1.0));
    563     G4double dy = zPathLength*(0.5*diry + rms*G4RandGauss::shoot(0.0,1.0));
    564     G4double dz;
    565     G4double d = (dx*dx + dy*dy)/(zPathLength*zPathLength);
    566     if(d < numlimit)  dz = -0.5*zPathLength*d*(1.0 + 0.25*d);
    567     else if(d < 1.0)  dz = -zPathLength*(1.0 - sqrt(1.0 - d));
    568     else {
    569       dx = dy = dz = 0.0;
    570     }
    571 
    572     temp.set(dx,dy,dz);
    573     if(isscat) temp.rotateUz(dir);
    574     pos += temp;
    575    
    576     pos.rotateUz(oldDirection);
    577 
     478  fParticleChange->ProposeMomentumDirection(dir);
     479
     480  // lateral displacement 
     481  if (latDisplasment) {
    578482    G4double r = pos.mag();
    579483
    … …  
    597501G4double G4WentzelVIModel::ComputeXSectionPerVolume()
    598502{
    599   const G4ElementVector* theElementVector =
    600     currentMaterial->GetElementVector();
     503  // prepare recomputation of x-sections
     504  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
    601505  const G4double* theAtomNumDensityVector =
    602506    currentMaterial->GetVecNbOfAtomsPerVolume();
    … …  
    604508  if(nelm > nelments) {
    605509    nelments = nelm;
    606     xsecn.resize(nelments);
    607     prob.resize(nelments);
    608   }
    609 
     510    xsecn.resize(nelm);
     511    prob.resize(nelm);
     512  }
     513  G4double cut = (*currentCuts)[currentMaterialIndex];
     514  cosTetMaxNuc = wokvi->GetCosThetaNuc();
     515
     516  // check consistency
    610517  xtsec = 0.0;
     518  if(cosTetMaxNuc > cosThetaMin) { return 0.0; }
     519
     520  // loop over elements
    611521  G4double xs = 0.0;
    612 
    613   for (G4int i=0; i<nelm; i++) {
    614     SetupTarget((*theElementVector)[i]->GetZ(), tkin);
     522  for (G4int i=0; i<nelm; ++i) {
     523    G4double costm =
     524      wokvi->SetupTarget(G4int((*theElementVector)[i]->GetZ()), cut);
    615525    G4double density = theAtomNumDensityVector[i];
    616     G4double cosnm = cosTetMaxNuc2;
    617     G4double cosem = cosTetMaxElec2;
    618 
    619     // recompute the angular limit
    620     cosTetMaxNuc2  = std::max(cosnm,cosThetaMin);
    621     cosTetMaxElec2 = std::max(cosem,cosThetaMin);
    622     xtsec += ComputeTransportXSectionPerAtom()*density;
    623     // return limit back
    624     cosTetMaxElec2 = cosem;
    625     cosTetMaxNuc2  = cosnm;
    626526
    627527    G4double esec = 0.0;
    628     G4double nsec = 0.0;
    629     G4double x1 = 1.0 - cosThetaMin + screenZ;
    630     G4double f  = kinFactor*density;
    631 
    632     // scattering off electrons
    633     if(cosThetaMin > cosem) {
    634       esec = f*(cosThetaMin - cosem)/(x1*(1.0 - cosem + screenZ));
    635     }
    636 
    637     // scattering off nucleaus
    638     if(cosThetaMin > cosnm) {
    639 
    640       // Rutherford part
    641       G4double s  = screenZ*formfactA;
    642       G4double z1 = 1.0 - cosnm + screenZ;
    643       G4double s1 = 1.0 - s;
    644       G4double d  = s1/formfactA;
    645 
    646       // check numerical limit
    647       if(d < numlimit*x1) {
    648         G4double x2 = x1*x1;
    649         G4double z2 = z1*z1;
    650         nsec = (1.0/(x1*x2) - 1.0/(z1*z2) - d*1.5*(1.0/(x2*x2) - 1.0/(z2*z2)))/
    651           (3.0*formfactA*formfactA);
    652       } else {
    653         G4double x2 = x1 + d;
    654         G4double z2 = z1 + d;
    655         nsec = (1.0/x1 - 1.0/z1 + 1.0/x2 - 1.0/z2 - 2.0*log(z1*x2/(z2*x1))/d)/(s1*s1);
    656       }
    657       nsec *= f*targetZ;
    658     }
    659     nsec += esec;
    660     if(nsec > 0.0) esec /= nsec;
    661     xs += nsec;
    662     xsecn[i] = xs;
     528    if(costm < cosThetaMin) { 
     529
     530      // recompute the transport x-section
     531      xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosThetaMin);
     532
     533      // recompute the total x-section
     534      G4double nsec = wokvi->ComputeNuclearCrossSection(cosThetaMin, costm);
     535      esec = wokvi->ComputeElectronCrossSection(cosThetaMin, costm);
     536      nsec += esec;
     537      if(nsec > 0.0) { esec /= nsec; }
     538      xtsec += nsec*density;
     539    }
     540    xsecn[i] = xtsec;
    663541    prob[i]  = esec;
    664     //G4cout << i << "  xs= " << xs << " cosThetaMin= " << cosThetaMin
    665     //     << " costm= " << costm << G4endl;
     542    //G4cout << i << "  xs= " << xs << " xtsec= " << xtsec << " 1-cosThetaMin= " << 1-cosThetaMin
     543    //     << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl;
    666544  }
    667545 
    668546  //G4cout << "ComputeXS result:  xsec(1/mm)= " << xs
    669   //<< " txsec(1/mm)= " << xtsec <<G4endl;
     547  //     << " txsec(1/mm)= " << xtsec <<G4endl;
    670548  return xs;
    671549}
    672550
    673551//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    674 
    675 /*
    676 G4double G4MuMscModel::ComputeXSectionPerVolume()
    677 {
    678   const G4ElementVector* theElementVector =
    679     currentMaterial->GetElementVector();
    680   const G4double* theAtomNumDensityVector =
    681     currentMaterial->GetVecNbOfAtomsPerVolume();
    682   size_t nelm = currentMaterial->GetNumberOfElements();
    683 
    684   xsece1 = 0.0;
    685   xsece2 = 0.0;
    686   xsecn2 = 0.0;
    687   zcorr  = 0.0;
    688 
    689   G4double fac = coeff*chargeSquare*invbeta2/mom2;
    690 
    691   for (size_t i=0; i<nelm; i++) {
    692     const G4Element* elm = (*theElementVector)[i];
    693     G4double Z = elm->GetZ();
    694     SetupTarget(Z, tkin);
    695     G4double den = fac*theAtomNumDensityVector[i]*Z;
    696 
    697     G4double x  = 1.0 - cosThetaMin;
    698     G4double x1 = x + screenZ;
    699     G4double x2 = 1.0/(x1*x1);
    700     G4double x3 = 1.0 + x*formfactA;
    701    
    702     //G4cout << "x= " << x << " den= " << den << " cosE= " << cosTetMaxElec << G4endl;
    703     //G4cout << "cosThtaMin= " << cosThetaMin << G4endl;
    704     //G4cout << "cosTetMaxNuc= " << cosTetMaxNuc << " q2Limit= " << q2Limit << G4endl;
    705    
    706     // scattering off electrons
    707     if(cosTetMaxElec < cosThetaMin) {
    708 
    709       // flat part
    710       G4double s = den*x2*x;
    711       xsece1 += s;
    712       zcorr  += 0.5*x*s;
    713 
    714       // Rutherford part
    715       G4double z1 = 1.0 - cosTetMaxElec + screenZ;
    716       G4double z2 = (cosThetaMin - cosTetMaxElec)/x1;
    717       if(z2 < 0.2) s = z2*(x - 0.5*z2*(x - screenZ))/x1;
    718       else         s = log(1.0 + z2)  - screenZ*z2/z1;
    719       xsece2  += den*z2/z1;
    720       zcorr   += den*s;
    721     }
    722     den *= Z;
    723 
    724     //G4cout << "Z= " << Z<< " cosL= " << cosTetMaxNuc << " cosMin= " << cosThetaMin << G4endl;
    725     // scattering off nucleaus
    726     if(cosTetMaxNuc < cosThetaMin) {
    727 
    728       // flat part
    729       G4double s = den*x2*x/(x3*x3);
    730       xsece1 += s;
    731       zcorr  += 0.5*x*s;
    732 
    733       // Rutherford part
    734       s  = screenZ*formfactA;
    735       G4double w  = 1.0 + 2.0*s;
    736       G4double z1 = 1.0 - cosTetMaxNuc + screenZ;
    737       G4double d  = (1.0 - s)/formfactA;
    738       G4double x4 = x1 + d;
    739       G4double z4 = z1 + d;
    740       G4double t1 = 1.0/(x1*z1);
    741       G4double t4 = 1.0/(x4*z4);
    742       G4double w1 = cosThetaMin - cosTetMaxNuc;
    743       G4double w2 = log(z1*x4/(x1*z4));
    744 
    745       den *= w;     
    746       xsecn2  += den*(w1*(t1 + t4) - 2.0*w2/d);
    747       zcorr   += den*(w*w2 - w1*(screenZ*t1 + t4/formfactA));
    748     }
    749     xsece[i] = xsece2;
    750     xsecn[i] = xsecn2;
    751     //    G4cout << i << "  xsece2= " << xsece2 << "  xsecn2= " << xsecn2 << G4endl;
    752   }
    753   G4double xsec = xsece1 + xsece2 + xsecn2;
    754  
    755     //G4cout << "xsece1= " << xsece1 << "  xsece2= " << xsece2
    756     //<< "  xsecn2= " << xsecn2
    757         // << " zsec= " << zcorr*0.5*tPathLength << G4endl;
    758   zcorr *= 0.5*tPathLength;
    759 
    760   return xsec;
    761 }
    762 */
    763 
    764 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
    765 
    766 void G4WentzelVIModel::ComputeMaxElectronScattering(G4double cutEnergy)
    767 {
    768   ecut = cutEnergy;
    769   G4double tmax = tkin;
    770   cosTetMaxElec = 1.0;
    771   if(mass > MeV) {
    772     G4double ratio = electron_mass_c2/mass;
    773     G4double tau = tkin/mass;
    774     tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
    775       (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
    776     cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
    777   } else {
    778 
    779     if(particle == theElectron) tmax *= 0.5;
    780     G4double t = std::min(cutEnergy, tmax);
    781     G4double mom21 = t*(t + 2.0*electron_mass_c2);
    782     G4double t1 = tkin - t;
    783     //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
    784     //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
    785     if(t1 > 0.0) {
    786       G4double mom22 = t1*(t1 + 2.0*mass);
    787       G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
    788       if(ctm <  1.0) cosTetMaxElec = ctm;
    789       if(ctm < -1.0) cosTetMaxElec = -1.0;
    790     }
    791   }
    792   if(cosTetMaxElec < cosTetMaxNuc) cosTetMaxElec = cosTetMaxNuc;
    793 }
    794 
    795 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracChangeset for help on using the changeset viewer.