G4WentzelVIModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4WentzelVIModel.cc 66596 2012-12-23 14:57:45Z vnivanch $
00027 //
00028 // -------------------------------------------------------------------
00029 //
00030 // GEANT4 Class file
00031 //
00032 //
00033 // File name:   G4WentzelVIModel
00034 //
00035 // Author:      V.Ivanchenko 
00036 //
00037 // Creation date: 09.04.2008 from G4MuMscModel
00038 //
00039 // Modifications:
00040 // 27-05-2010 V.Ivanchenko added G4WentzelOKandVIxSection class to
00041 //              compute cross sections and sample scattering angle
00042 //
00043 //
00044 // Class Description:
00045 //
00046 // Implementation of the model of multiple scattering based on
00047 // G.Wentzel, Z. Phys. 40 (1927) 590.
00048 // H.W.Lewis, Phys Rev 78 (1950) 526.
00049 // J.M. Fernandez-Varea et al., NIM B73 (1993) 447.
00050 // L.Urban, CERN-OPEN-2006-077.
00051 
00052 // -------------------------------------------------------------------
00053 //
00054 
00055 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00056 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00057 
00058 #include "G4WentzelVIModel.hh"
00059 #include "G4PhysicalConstants.hh"
00060 #include "G4SystemOfUnits.hh"
00061 #include "Randomize.hh"
00062 #include "G4ParticleChangeForMSC.hh"
00063 #include "G4PhysicsTableHelper.hh"
00064 #include "G4ElementVector.hh"
00065 #include "G4ProductionCutsTable.hh"
00066 #include "G4LossTableManager.hh"
00067 #include "G4Pow.hh"
00068 #include "G4NistManager.hh"
00069 
00070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00071 
00072 using namespace std;
00073 
00074 G4WentzelVIModel::G4WentzelVIModel(const G4String& nam) :
00075   G4VMscModel(nam),
00076   numlimit(0.1),
00077   currentCouple(0),
00078   cosThetaMin(1.0),
00079   inside(false),
00080   singleScatteringMode(false)
00081 {
00082   invsqrt12 = 1./sqrt(12.);
00083   tlimitminfix = 1.e-6*mm;
00084   lowEnergyLimit = 1.0*eV;
00085   particle = 0;
00086   nelments = 5;
00087   xsecn.resize(nelments);
00088   prob.resize(nelments);
00089   theManager = G4LossTableManager::Instance();
00090   fNistManager = G4NistManager::Instance();
00091   fG4pow = G4Pow::GetInstance();
00092   wokvi = new G4WentzelOKandVIxSection();
00093 
00094   preKinEnergy = tPathLength = zPathLength = lambdaeff = currentRange = xtsec = 0;
00095   currentMaterialIndex = 0;
00096   cosThetaMax = cosTetMaxNuc = 1.0;
00097 
00098   fParticleChange = 0;
00099   currentCuts = 0;
00100   currentMaterial = 0;
00101 }
00102 
00103 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00104 
00105 G4WentzelVIModel::~G4WentzelVIModel()
00106 {
00107   delete wokvi;
00108 }
00109 
00110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00111 
00112 void G4WentzelVIModel::Initialise(const G4ParticleDefinition* p,
00113                                   const G4DataVector& cuts)
00114 {
00115   // reset parameters
00116   SetupParticle(p);
00117   currentRange = 0.0;
00118 
00119   cosThetaMax = cos(PolarAngleLimit());
00120   wokvi->Initialise(p, cosThetaMax);
00121   /*
00122   G4cout << "G4WentzelVIModel: " << particle->GetParticleName()
00123          << "  1-cos(ThetaLimit)= " << 1 - cosThetaMax 
00124          << G4endl;
00125   */
00126   currentCuts = &cuts;
00127 
00128   // set values of some data members
00129   fParticleChange = GetParticleChangeForMSC(p);
00130 }
00131 
00132 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00133 
00134 G4double G4WentzelVIModel::ComputeCrossSectionPerAtom( 
00135                              const G4ParticleDefinition* p,
00136                              G4double kinEnergy,
00137                              G4double Z, G4double,
00138                              G4double cutEnergy, G4double)
00139 {
00140   G4double cross = 0.0;
00141   if(p != particle) { SetupParticle(p); }
00142   if(kinEnergy < lowEnergyLimit) { return cross; }
00143   if(!CurrentCouple()) {
00144     G4Exception("G4WentzelVIModel::ComputeCrossSectionPerAtom", "em0011",
00145                 FatalException, " G4MaterialCutsCouple is not defined");
00146     return 0.0;
00147   }
00148   DefineMaterial(CurrentCouple());
00149   cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
00150   if(cosTetMaxNuc < 1.0) {
00151     cosTetMaxNuc = wokvi->SetupTarget(G4lrint(Z), cutEnergy);
00152     cross = wokvi->ComputeTransportCrossSectionPerAtom(cosTetMaxNuc);
00153     /*
00154     if(p->GetParticleName() == "e-")      
00155     G4cout << "G4WentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= " << kinEnergy 
00156            << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
00157            << " " << particle->GetParticleName() << G4endl;
00158     */
00159   }
00160   return cross;
00161 }
00162 
00163 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00164 
00165 void G4WentzelVIModel::StartTracking(G4Track* track)
00166 {
00167   SetupParticle(track->GetDynamicParticle()->GetDefinition());
00168   inside = false;
00169 }
00170 
00171 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00172 
00173 G4double G4WentzelVIModel::ComputeTruePathLengthLimit(
00174                              const G4Track& track,
00175                              G4double& currentMinimalStep)
00176 {
00177   G4double tlimit = currentMinimalStep;
00178   const G4DynamicParticle* dp = track.GetDynamicParticle();
00179   G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
00180   G4StepStatus stepStatus = sp->GetStepStatus();
00181   singleScatteringMode = false;
00182   //G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit stepStatus= " 
00183   //     << stepStatus << G4endl;
00184 
00185   // initialisation for each step, lambda may be computed from scratch
00186   preKinEnergy  = dp->GetKineticEnergy();
00187   DefineMaterial(track.GetMaterialCutsCouple());
00188   lambdaeff = GetTransportMeanFreePath(particle,preKinEnergy);
00189   currentRange = GetRange(particle,preKinEnergy,currentCouple);
00190   cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial);
00191 
00192   //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange
00193   //     << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl;
00194 
00195   // extra check for abnormal situation
00196   // this check needed to run MSC with eIoni and eBrem inactivated
00197   if(tlimit > currentRange) { tlimit = currentRange; }
00198 
00199   // stop here if small range particle
00200   if(inside || tlimit < tlimitminfix) { 
00201     return ConvertTrueToGeom(tlimit, currentMinimalStep); 
00202   }
00203 
00204   // pre step
00205   G4double presafety = sp->GetSafety();
00206   // far from geometry boundary
00207   if(currentRange < presafety) {
00208     inside = true;  
00209     return ConvertTrueToGeom(tlimit, currentMinimalStep);
00210   }
00211 
00212   // compute presafety again if presafety <= 0 and no boundary
00213   // i.e. when it is needed for optimization purposes
00214   if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
00215     presafety = ComputeSafety(sp->GetPosition(), tlimit); 
00216     if(currentRange < presafety) {
00217       inside = true;  
00218       return ConvertTrueToGeom(tlimit, currentMinimalStep);
00219     }
00220   }
00221   /*
00222   G4cout << "e(MeV)= " << preKinEnergy/MeV
00223          << "  " << particle->GetParticleName() 
00224          << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
00225          << " R(mm)= " <<currentRange/mm
00226          << " L0(mm^-1)= " << lambdaeff*mm 
00227          <<G4endl;
00228   */
00229 
00230   // natural limit for high energy
00231   G4double rlimit = std::max(facrange*currentRange, 
00232                              0.7*(1.0 - cosTetMaxNuc)*lambdaeff);
00233 
00234   // low-energy e-
00235   if(cosThetaMax > cosTetMaxNuc) {
00236     rlimit = std::min(rlimit, facsafety*presafety);
00237   }
00238    
00239   // cut correction
00240   G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
00241   //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= " 
00242   // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax 
00243   //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl;
00244   if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
00245  
00246   if(rlimit < tlimit) { tlimit = rlimit; }
00247 
00248   tlimit = std::max(tlimit, tlimitminfix);
00249 
00250   // step limit in infinite media
00251   tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
00252 
00253   //compute geomlimit and force few steps within a volume
00254   if (steppingAlgorithm == fUseDistanceToBoundary 
00255       && stepStatus == fGeomBoundary) {
00256 
00257     G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
00258     tlimit = std::min(tlimit, geomlimit/facgeom);
00259   } 
00260 
00261   /*  
00262   G4cout << particle->GetParticleName() << " e= " << preKinEnergy
00263          << " L0= " << lambdaeff << " R= " << currentRange
00264          << "tlimit= " << tlimit  
00265          << " currentMinimalStep= " << currentMinimalStep << G4endl;
00266   */
00267   return ConvertTrueToGeom(tlimit, currentMinimalStep);
00268 }
00269 
00270 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00271 
00272 G4double G4WentzelVIModel::ComputeGeomPathLength(G4double truelength)
00273 {
00274   tPathLength  = truelength;
00275   zPathLength  = tPathLength;
00276 
00277   if(lambdaeff > 0.0 && lambdaeff != DBL_MAX) {
00278     G4double tau = tPathLength/lambdaeff;
00279     //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
00280     //   << " Leff= " << lambdaeff << " tau= " << tau << G4endl; 
00281     // small step
00282     if(tau < numlimit) {
00283       zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
00284 
00285       // medium step
00286     } else {
00287       G4double e1 = 0.0;
00288       if(currentRange > tPathLength) {
00289         e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
00290       }
00291       e1 = 0.5*(e1 + preKinEnergy);
00292       cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
00293       lambdaeff = GetTransportMeanFreePath(particle,e1);
00294       zPathLength = lambdaeff*(1.0 - exp(-tPathLength/lambdaeff));
00295     }
00296   } else { lambdaeff = DBL_MAX; }
00297   //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
00298   return zPathLength;
00299 }
00300 
00301 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00302 
00303 G4double G4WentzelVIModel::ComputeTrueStepLength(G4double geomStepLength)
00304 {
00305   // initialisation of single scattering x-section
00306   xtsec = 0.0;
00307   cosThetaMin = cosTetMaxNuc;
00308   /*
00309   G4cout << "ComputeTrueStepLength: Step= " << geomStepLength 
00310          << "  Lambda= " <<  lambdaeff 
00311          << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl;
00312   */
00313   // pathalogical case
00314   if(lambdaeff == DBL_MAX) { 
00315     singleScatteringMode = true;
00316     zPathLength  = geomStepLength;
00317     tPathLength  = geomStepLength;
00318     cosThetaMin = 1.0;
00319 
00320     // normal case
00321   } else {
00322 
00323     // small step use only single scattering
00324     const G4double singleScatLimit = 1.0e-7;
00325     if(geomStepLength < lambdaeff*singleScatLimit*(1.0 - cosTetMaxNuc)) {
00326       singleScatteringMode = true;
00327       cosThetaMin = 1.0;
00328       zPathLength  = geomStepLength;
00329       tPathLength  = geomStepLength;
00330 
00331       // step defined by transportation 
00332     } else if(geomStepLength != zPathLength) { 
00333 
00334       // step defined by transportation 
00335       zPathLength  = geomStepLength;
00336       G4double tau = geomStepLength/lambdaeff;
00337       tPathLength  = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0); 
00338 
00339       // energy correction for a big step
00340       if(tau > numlimit) {
00341         G4double e1 = 0.0;
00342         if(currentRange > tPathLength) {
00343           e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
00344         }
00345         e1 = 0.5*(e1 + preKinEnergy);
00346         cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
00347         lambdaeff = GetTransportMeanFreePath(particle,e1);
00348         tau = zPathLength/lambdaeff;
00349       
00350         if(tau < 0.999999) { tPathLength = -lambdaeff*log(1.0 - tau); } 
00351         else               { tPathLength = currentRange; }
00352       }
00353     }
00354   }
00355 
00356   // check of step length
00357   // define threshold angle between single and multiple scattering 
00358   if(!singleScatteringMode) { cosThetaMin = 1.0 - 1.5*tPathLength/lambdaeff; }
00359 
00360   // recompute transport cross section - do not change energy
00361   // anymore - cannot be applied for big steps
00362   if(cosThetaMin > cosTetMaxNuc) {
00363 
00364     // new computation
00365     G4double cross = ComputeXSectionPerVolume();
00366     //G4cout << "%%%% cross= " << cross << "  xtsec= " << xtsec << G4endl;
00367     if(cross <= 0.0) {
00368       singleScatteringMode = true;
00369       tPathLength = zPathLength; 
00370       lambdaeff = DBL_MAX;
00371       cosThetaMin = 1.0;
00372     } else if(xtsec > 0.0) {
00373 
00374       lambdaeff = 1./cross; 
00375       G4double tau = zPathLength*cross;
00376       if(tau < numlimit) { tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0); } 
00377       else if(tau < 0.999999) { tPathLength = -lambdaeff*log(1.0 - tau); } 
00378       else                    { tPathLength = currentRange; }
00379 
00380       if(tPathLength > currentRange) { tPathLength = currentRange; }
00381     } 
00382   }
00383 
00384   /*   
00385   G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
00386          <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
00387   G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
00388          << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc 
00389          << " e(MeV)= " << preKinEnergy/MeV << "  "  
00390          << singleScatteringMode << G4endl;
00391   */
00392   return tPathLength;
00393 }
00394 
00395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00396 
00397 G4ThreeVector& 
00398 G4WentzelVIModel::SampleScattering(const G4ThreeVector& oldDirection,
00399                                    G4double safety)
00400 {
00401   fDisplacement.set(0.0,0.0,0.0);
00402   //G4cout << "!##! G4WentzelVIModel::SampleScattering for " 
00403   //     << particle->GetParticleName() << G4endl;
00404 
00405   // ignore scattering for zero step length and energy below the limit
00406   if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0) 
00407     { return fDisplacement; }
00408   
00409   G4double invlambda = 0.0;
00410   if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
00411 
00412   // use average kinetic energy over the step
00413   G4double cut = (*currentCuts)[currentMaterialIndex];
00414   /*  
00415   G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
00416          << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec 
00417          << " xmsc= " <<  tPathLength*invlambda 
00418          << " safety= " << safety << G4endl;
00419   */
00420 
00421   // step limit due msc
00422   G4double x0 = tPathLength;
00423   //  const G4double thinlimit = 0.5; 
00424   const G4double thinlimit = 0.1; 
00425   G4int nMscSteps = 1;
00426   // large scattering angle case - two step approach
00427   if(!singleScatteringMode && tPathLength*invlambda > thinlimit 
00428      && safety > tlimitminfix) { 
00429     x0 *= 0.5; 
00430     nMscSteps = 2; 
00431   } 
00432 
00433   // step limit due to single scattering
00434   G4double x1 = 2*tPathLength;
00435   if(0.0 < xtsec) { x1 = -log(G4UniformRand())/xtsec; }
00436 
00437   // no scattering case
00438   if(singleScatteringMode && x1 > tPathLength)  
00439     { return fDisplacement; }
00440 
00441   const G4ElementVector* theElementVector = 
00442     currentMaterial->GetElementVector();
00443   G4int nelm = currentMaterial->GetNumberOfElements();
00444 
00445   // geometry
00446   G4double sint, cost, phi;
00447   G4ThreeVector temp(0.0,0.0,1.0);
00448 
00449   // current position and direction relative to the end point
00450   // because of magnetic field geometry is computed relatively to the 
00451   // end point of the step 
00452   G4ThreeVector dir(0.0,0.0,1.0);
00453   fDisplacement.set(0.0,0.0,-zPathLength);
00454 
00455   G4double mscfac = zPathLength/tPathLength;
00456 
00457   // start a loop 
00458   G4double x2 = x0;
00459   G4double step, z;
00460   G4bool singleScat;
00461   /*  
00462       G4cout << "Start of the loop x1(mm)= " << x1 << "  x2(mm)= " << x2 
00463       << " 1-cost1= " << 1 - cosThetaMin << "  " << singleScatteringMode 
00464       << " xtsec= " << xtsec << G4endl;
00465   */
00466   do {
00467 
00468     // single scattering case
00469     if(singleScatteringMode || x1 <= x2) { 
00470       step = x1;
00471       singleScat = true;
00472     } else {
00473       step = x2;
00474       singleScat = false;
00475     }
00476 
00477     // new position
00478     fDisplacement += step*mscfac*dir;
00479 
00480     if(singleScat) {
00481 
00482       // select element
00483       G4int i = 0;
00484       if(nelm > 1) {
00485         G4double qsec = G4UniformRand()*xtsec;
00486         for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
00487       }
00488       G4double cosTetM = 
00489         wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
00490       //G4cout << "!!! " << cosThetaMin << "  " << cosTetM << "  " << prob[i] << G4endl;
00491       temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
00492 
00493       // direction is changed
00494       temp.rotateUz(dir);
00495       dir = temp;
00496 
00497       // new proposed step length
00498       x2 -= step; 
00499       x1  = -log(G4UniformRand())/xtsec; 
00500       if(singleScatteringMode && x1 > x2) { break; }
00501 
00502     // multiple scattering
00503     } else { 
00504       --nMscSteps;
00505       x1 -= step;
00506       x2  = x0;
00507 
00508       // for multiple scattering x0 should be used as a step size
00509       G4double z0 = x0*invlambda;
00510 
00511       // sample z in interval 0 - 1
00512       do {
00513         G4double zzz = 0.0;
00514         if(z0 > 0.1) { zzz = exp(-1.0/z0); }
00515         z = -z0*log(1.0 - (1.0 - zzz)*G4UniformRand());
00516       } while(z > 1.0);
00517       cost = 1.0 - 2.0*z/*factCM*/;
00518       if(cost > 1.0)       { cost = 1.0; }
00519       else if(cost < -1.0) { cost =-1.0; }
00520       sint = sqrt((1.0 - cost)*(1.0 + cost));
00521       phi  = twopi*G4UniformRand();
00522       G4double vx1 = sint*cos(phi);
00523       G4double vy1 = sint*sin(phi);
00524 
00525       // change direction
00526       temp.set(vx1,vy1,cost);
00527       temp.rotateUz(dir);
00528       dir = temp;
00529 
00530       // lateral displacement  
00531       if (latDisplasment && safety > tlimitminfix) {
00532         G4double rms = invsqrt12*sqrt(2*z0);
00533         G4double r   = x0*mscfac;
00534         G4double dx  = r*(0.5*vx1 + rms*G4RandGauss::shoot(0.0,1.0));
00535         G4double dy  = r*(0.5*vy1 + rms*G4RandGauss::shoot(0.0,1.0));
00536         G4double dz;
00537         G4double d   = r*r - dx*dx - dy*dy;
00538         if(d >= 0.0)  { dz = sqrt(d) - r; }
00539         else          { dx = dy = dz = 0.0; }
00540 
00541         // change position
00542         temp.set(dx,dy,dz);
00543         temp.rotateUz(dir); 
00544         fDisplacement += temp;
00545       }
00546     }
00547   } while (0 < nMscSteps);
00548     
00549   dir.rotateUz(oldDirection);
00550 
00551   //G4cout << "G4WentzelVIModel sampling of scattering is done" << G4endl;
00552   // end of sampling -------------------------------
00553 
00554   fParticleChange->ProposeMomentumDirection(dir);
00555 
00556   // lateral displacement  
00557   fDisplacement.rotateUz(oldDirection);
00558 
00559   /*            
00560          G4cout << " r(mm)= " << fDisplacement.mag() 
00561                 << " safety= " << safety
00562                 << " trueStep(mm)= " << tPathLength
00563                 << " geomStep(mm)= " << zPathLength
00564                 << " x= " << fDisplacement.x() 
00565                 << " y= " << fDisplacement.y() 
00566                 << " z= " << fDisplacement.z()
00567                 << G4endl;
00568   */
00569 
00570   //G4cout<< "G4WentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl;
00571   return fDisplacement;
00572 }
00573 
00574 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
00575 
00576 G4double G4WentzelVIModel::ComputeXSectionPerVolume()
00577 {
00578   // prepare recomputation of x-sections
00579   const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
00580   const G4double* theAtomNumDensityVector = 
00581     currentMaterial->GetVecNbOfAtomsPerVolume();
00582   G4int nelm = currentMaterial->GetNumberOfElements();
00583   if(nelm > nelments) {
00584     nelments = nelm;
00585     xsecn.resize(nelm);
00586     prob.resize(nelm);
00587   }
00588   G4double cut = (*currentCuts)[currentMaterialIndex];
00589   //  cosTetMaxNuc = wokvi->GetCosThetaNuc();
00590 
00591   // check consistency
00592   xtsec = 0.0;
00593   if(cosTetMaxNuc > cosThetaMin) { return 0.0; }
00594 
00595   // loop over elements
00596   G4double xs = 0.0;
00597   for (G4int i=0; i<nelm; ++i) {
00598     G4double costm = 
00599       wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
00600     G4double density = theAtomNumDensityVector[i];
00601 
00602     G4double esec = 0.0;
00603     if(costm < cosThetaMin) {  
00604 
00605       // recompute the transport x-section
00606       if(1.0 > cosThetaMin) {
00607         xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosThetaMin);
00608       }
00609       // recompute the total x-section
00610       G4double nucsec = wokvi->ComputeNuclearCrossSection(cosThetaMin, costm);
00611       esec = wokvi->ComputeElectronCrossSection(cosThetaMin, costm);
00612       nucsec += esec;
00613       if(nucsec > 0.0) { esec /= nucsec; }
00614       xtsec += nucsec*density;
00615     }
00616     xsecn[i] = xtsec;
00617     prob[i]  = esec;
00618     //G4cout << i << "  xs= " << xs << " xtsec= " << xtsec 
00619     //       << " 1-cosThetaMin= " << 1-cosThetaMin 
00620     //     << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl;
00621   }
00622   
00623   //G4cout << "ComputeXS result:  xsec(1/mm)= " << xs 
00624   //     << " txsec(1/mm)= " << xtsec <<G4endl; 
00625   return xs;
00626 }
00627 
00628 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......

Generated on Mon May 27 17:50:25 2013 for Geant4 by  doxygen 1.4.7