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

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