Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/standard/src/G4WentzelVIModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 //
 27 // -------------------------------------------------------------------
 28 //
 29 // GEANT4 Class file
 30 //
 31 //
 32 // File name:   G4WentzelVIModel
 33 //
 34 // Author:      V.Ivanchenko 
 35 //
 36 // Creation date: 09.04.2008 from G4MuMscModel
 37 //
 38 // Modifications:
 39 // 27-05-2010 V.Ivanchenko added G4WentzelOKandVIxSection class to
 40 //              compute cross sections and sample scattering angle
 41 //
 42 //
 43 // Class Description:
 44 //
 45 // Implementation of the model of multiple scattering based on
 46 // G.Wentzel, Z. Phys. 40 (1927) 590.
 47 // H.W.Lewis, Phys Rev 78 (1950) 526.
 48 // J.M. Fernandez-Varea et al., NIM B73 (1993) 447.
 49 // L.Urban, CERN-OPEN-2006-077.
 50 
 51 // -------------------------------------------------------------------
 52 //
 53 
 54 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 56 
 57 #include "G4WentzelVIModel.hh"
 58 #include "G4PhysicalConstants.hh"
 59 #include "G4SystemOfUnits.hh"
 60 #include "Randomize.hh"
 61 #include "G4ParticleChangeForMSC.hh"
 62 #include "G4PhysicsTableHelper.hh"
 63 #include "G4ElementVector.hh"
 64 #include "G4ProductionCutsTable.hh"
 65 #include "G4EmParameters.hh"
 66 #include "G4Log.hh"
 67 #include "G4Exp.hh"
 68 
 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 70 
 71 using namespace std;
 72 
 73 const G4double invsqrt12 = 1./std::sqrt(12.);
 74 const G4double numlimit = 0.1;
 75 const G4int minNCollisions = 10;
 76 
 77 G4WentzelVIModel::G4WentzelVIModel(G4bool comb, const G4String& nam) 
 78   : G4VMscModel(nam),
 79     singleScatteringMode(false),
 80     isCombined(comb),
 81     useSecondMoment(false)
 82 {
 83   tlimitminfix = 1.e-6*CLHEP::mm;
 84   lowEnergyLimit = 1.0*CLHEP::eV;
 85   SetSingleScatteringFactor(1.25);
 86   wokvi = new G4WentzelOKandVIxSection(isCombined);
 87 }
 88 
 89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 90 
 91 G4WentzelVIModel::~G4WentzelVIModel()
 92 {
 93   delete wokvi;
 94   if(IsMaster()) {
 95     delete fSecondMoments;
 96     fSecondMoments = nullptr;
 97   }
 98 }
 99 
100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
101 
102 void G4WentzelVIModel::Initialise(const G4ParticleDefinition* p,
103                                   const G4DataVector& cuts)
104 {
105   // reset parameters
106   SetupParticle(p);
107   InitialiseParameters(p);
108   currentRange = 0.0;
109 
110   if(isCombined) {
111     G4double tet = PolarAngleLimit();
112     if(tet <= 0.0)           { cosThetaMax = 1.0; }
113     else if(tet < CLHEP::pi) { cosThetaMax = cos(tet); }
114   }
115   //G4cout << "G4WentzelVIModel::Initialise " << p->GetParticleName() 
116   //   << " " << this << " " << wokvi << G4endl;
117 
118   wokvi->Initialise(p, cosThetaMax);
119   /* 
120   G4cout << "G4WentzelVIModel: " << particle->GetParticleName()
121          << "  1-cos(ThetaLimit)= " << 1 - cosThetaMax 
122          << " SingScatFactor= " << ssFactor
123          << G4endl;
124   */
125   currentCuts = &cuts;
126 
127   // set values of some data members
128   fParticleChange = GetParticleChangeForMSC(p);
129 
130   // Access to materials
131   const G4ProductionCutsTable* theCoupleTable =
132     G4ProductionCutsTable::GetProductionCutsTable();
133   G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
134   nelments = 0;
135   for(G4int i=0; i<numOfCouples; ++i) {
136     G4int nelm = (G4int)theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetNumberOfElements();
137     nelments = std::max(nelments, nelm);
138   }
139   xsecn.resize(nelments);
140   prob.resize(nelments);
141 
142   // build second moment table only if transport table is build
143   G4PhysicsTable* table = GetCrossSectionTable();
144   if(useSecondMoment && IsMaster() && nullptr != table) {
145 
146     //G4cout << "### G4WentzelVIModel::Initialise: build 2nd moment table "
147     //           << table << G4endl;
148     fSecondMoments =  
149       G4PhysicsTableHelper::PreparePhysicsTable(fSecondMoments);
150 
151     G4bool splineFlag = true;
152     G4PhysicsVector* aVector = nullptr;
153     G4PhysicsVector* bVector = nullptr;
154     G4double emin = std::max(LowEnergyLimit(), LowEnergyActivationLimit());
155     G4double emax = std::min(HighEnergyLimit(), HighEnergyActivationLimit());
156     if(emin < emax) {
157       std::size_t n = G4EmParameters::Instance()->NumberOfBinsPerDecade()
158         *G4lrint(std::log10(emax/emin));
159       if(n < 3) { n = 3; }
160 
161       for(G4int i=0; i<numOfCouples; ++i) {
162 
163         //G4cout<< "i= " << i << " Flag=  " << fSecondMoments->GetFlag(i) 
164         //      << G4endl;
165         if(fSecondMoments->GetFlag(i)) {
166           DefineMaterial(theCoupleTable->GetMaterialCutsCouple(i));
167        
168           delete (*fSecondMoments)[i];
169           if(nullptr == aVector) { 
170             aVector = new G4PhysicsLogVector(emin, emax, n, splineFlag);
171             bVector = aVector;
172           } else {
173             bVector = new G4PhysicsVector(*aVector);
174           }
175           for(std::size_t j=0; j<n; ++j) {
176             G4double e = bVector->Energy(j); 
177             bVector->PutValue(j, ComputeSecondMoment(p, e)*e*e);
178           }
179           if(splineFlag) { bVector->FillSecondDerivatives(); }
180           (*fSecondMoments)[i] = bVector;  
181         }
182       }
183     } 
184     //G4cout << *fSecondMoments << G4endl;
185   }
186 }
187 
188 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
189 
190 void G4WentzelVIModel::InitialiseLocal(const G4ParticleDefinition*, 
191                                        G4VEmModel* masterModel)
192 {
193   fSecondMoments = static_cast<G4WentzelVIModel*>(masterModel)
194     ->GetSecondMomentTable(); 
195 }
196 
197 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198 
199 void G4WentzelVIModel::DefineMaterial(const G4MaterialCutsCouple* cup) 
200 { 
201   if(cup != currentCouple) {
202     currentCouple = cup;
203     SetCurrentCouple(cup); 
204     currentMaterial = cup->GetMaterial();
205     currentMaterialIndex = currentCouple->GetIndex(); 
206   }
207 }
208 
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
210 
211 G4double G4WentzelVIModel::ComputeCrossSectionPerAtom( 
212                              const G4ParticleDefinition* p,
213                              G4double kinEnergy,
214                              G4double Z, G4double,
215                              G4double cutEnergy, G4double)
216 {
217   G4double cross = 0.0;
218   SetupParticle(p); 
219   if(kinEnergy < lowEnergyLimit) { return cross; }
220   if(nullptr == CurrentCouple()) {
221     G4Exception("G4WentzelVIModel::ComputeCrossSectionPerAtom", "em0011",
222                 FatalException, " G4MaterialCutsCouple is not defined");
223     return 0.0;
224   }
225   DefineMaterial(CurrentCouple());
226   cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
227   if(cosTetMaxNuc < 1.0) {
228     G4double cut  = (0.0 < fixedCut) ? fixedCut : cutEnergy;
229     G4double cost = wokvi->SetupTarget(G4lrint(Z), cut);
230     cross = wokvi->ComputeTransportCrossSectionPerAtom(cost);
231     /*
232     if(p->GetParticleName() == "e-")      
233     G4cout << "G4WentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= "<<kinEnergy 
234            << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
235            << " " << particle->GetParticleName() << G4endl;
236     */
237   }
238   return cross;
239 }
240 
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242 
243 void G4WentzelVIModel::StartTracking(G4Track* track)
244 {
245   /*
246   G4cout << "G4WentzelVIModel::StartTracking " << track << "  " << this << "  "
247    << track->GetParticleDefinition()->GetParticleName() 
248    << "   workvi: " << wokvi << G4endl; 
249   */
250   SetupParticle(track->GetParticleDefinition());
251 }
252 
253 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
254 
255 G4double G4WentzelVIModel::ComputeTruePathLengthLimit(
256                              const G4Track& track,
257                              G4double& currentMinimalStep)
258 {
259   G4double tlimit = currentMinimalStep;
260   const G4DynamicParticle* dp = track.GetDynamicParticle();
261   const G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
262   G4StepStatus stepStatus = sp->GetStepStatus();
263   singleScatteringMode = false;
264 
265   //G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit stepStatus= " 
266   //         << stepStatus << "  " << track.GetDefinition()->GetParticleName() 
267   //         << G4endl;
268 
269   // initialisation for each step, lambda may be computed from scratch
270   preKinEnergy = dp->GetKineticEnergy();
271   effKinEnergy = preKinEnergy;
272   DefineMaterial(track.GetMaterialCutsCouple());
273   const G4double logPreKinEnergy = dp->GetLogKineticEnergy();
274   lambdaeff = GetTransportMeanFreePath(particle,preKinEnergy,logPreKinEnergy);
275   currentRange = GetRange(particle,preKinEnergy,currentCouple,logPreKinEnergy);
276   cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial);
277   
278   //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange
279   // << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl;
280   
281   // extra check for abnormal situation
282   // this check needed to run MSC with eIoni and eBrem inactivated
283   if(tlimit > currentRange) { tlimit = currentRange; }
284 
285   // stop here if small range particle
286   if(tlimit < tlimitminfix) { 
287     return ConvertTrueToGeom(tlimit, currentMinimalStep); 
288   }
289 
290   // pre step
291   G4double presafety = sp->GetSafety();
292   // far from geometry boundary
293   if(currentRange < presafety) {
294     return ConvertTrueToGeom(tlimit, currentMinimalStep);
295   }
296 
297   // compute presafety again if presafety <= 0 and no boundary
298   // i.e. when it is needed for optimization purposes
299   if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
300     presafety = ComputeSafety(sp->GetPosition(), tlimit); 
301     if(currentRange < presafety) {
302       return ConvertTrueToGeom(tlimit, currentMinimalStep);
303     }
304   }
305   /*       
306   G4cout << "e(MeV)= " << preKinEnergy/MeV
307          << "  " << particle->GetParticleName() 
308          << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
309          << " R(mm)= " <<currentRange/mm
310          << " L0(mm^-1)= " << lambdaeff*mm 
311          << G4endl;
312   */
313   // natural limit for high energy
314   G4double rlimit = std::max(facrange*currentRange, 
315                              (1.0 - cosTetMaxNuc)*lambdaeff*invssFactor);
316 
317   // low-energy e-
318   if(cosThetaMax > cosTetMaxNuc) {
319     rlimit = std::min(rlimit, facsafety*presafety);
320   }
321    
322   // cut correction
323   G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
324   //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= " 
325   // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax 
326   //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl;
327   if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*std::sqrt(rlimit/rcut)); }
328  
329   tlimit = std::min(tlimit, rlimit);
330   tlimit = std::max(tlimit, tlimitminfix);
331 
332   // step limit in infinite media
333   tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
334 
335   //compute geomlimit and force few steps within a volume
336   if (steppingAlgorithm == fUseDistanceToBoundary 
337       && stepStatus == fGeomBoundary) {
338 
339     G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
340     tlimit = std::min(tlimit, geomlimit/facgeom);
341   } 
342   /*         
343   G4cout << particle->GetParticleName() << " e= " << preKinEnergy
344          << " L0= " << lambdaeff << " R= " << currentRange
345          << " tlimit= " << tlimit  
346            << " currentMinimalStep= " << currentMinimalStep << G4endl;
347   */
348   return ConvertTrueToGeom(tlimit, currentMinimalStep);
349 }
350 
351 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
352 
353 G4double G4WentzelVIModel::ComputeGeomPathLength(G4double truelength)
354 {
355   zPathLength = tPathLength = truelength;
356 
357   // small step use only single scattering
358   cosThetaMin = 1.0;
359   ComputeTransportXSectionPerVolume(cosThetaMin);
360   //G4cout << "xtsec= " << xtsec << "  Nav= " 
361   //         << zPathLength*xtsec << G4endl;
362   if(0.0 >= lambdaeff || G4int(zPathLength*xtsec) < minNCollisions) {
363     singleScatteringMode = true;
364     lambdaeff = DBL_MAX;
365 
366   } else {
367     //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
368     //           << " Leff= " << lambdaeff << G4endl; 
369     // small step
370     if(tPathLength < numlimit*lambdaeff) {
371       G4double tau = tPathLength/lambdaeff;
372       zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
373 
374       // medium step
375     } else {
376       G4double e1 = 0.0;
377       if(currentRange > tPathLength) {
378         e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
379       }
380       effKinEnergy = 0.5*(e1 + preKinEnergy);
381       cosTetMaxNuc = wokvi->SetupKinematic(effKinEnergy, currentMaterial);
382       lambdaeff = GetTransportMeanFreePath(particle, effKinEnergy);
383       //G4cout << " tLength= "<< tPathLength<< " Leff= " << lambdaeff << G4endl;
384       zPathLength = lambdaeff;
385       if(tPathLength*numlimit < lambdaeff) {
386         zPathLength *= (1.0 - G4Exp(-tPathLength/lambdaeff));
387       }
388     }
389   }
390   //G4cout << "Comp.geom: zLength= "<<zPathLength<<" tLength= "
391   //         << tPathLength<< " Leff= " << lambdaeff << G4endl;
392   return zPathLength;
393 }
394 
395 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
396 
397 G4double G4WentzelVIModel::ComputeTrueStepLength(G4double geomStepLength)
398 {
399   // initialisation of single scattering x-section
400   /*
401   G4cout << "ComputeTrueStepLength: Step= " << geomStepLength 
402          << "  geomL= " << zPathLength
403          << "  Lambda= " <<  lambdaeff 
404            << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl;
405   */
406   if(singleScatteringMode) {
407     zPathLength = tPathLength = geomStepLength;
408 
409   } else {
410 
411     // step defined by transportation
412     // change both geom and true step lengths 
413     if(geomStepLength < zPathLength) { 
414 
415       // single scattering
416       if(G4int(geomStepLength*xtsec) < minNCollisions) {
417         zPathLength = tPathLength = geomStepLength;
418         lambdaeff = DBL_MAX;
419         singleScatteringMode = true;
420 
421         // multiple scattering
422       } else {
423         // small step
424         if(geomStepLength < numlimit*lambdaeff) {
425           G4double tau = geomStepLength/lambdaeff;
426           tPathLength = geomStepLength*(1.0 + 0.5*tau + tau*tau/3.0); 
427 
428           // energy correction for a big step
429         } else {
430           tPathLength *= geomStepLength/zPathLength;
431           G4double e1 = 0.0;
432           if(currentRange > tPathLength) {
433             e1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
434           }
435           effKinEnergy = 0.5*(e1 + preKinEnergy);
436           cosTetMaxNuc = wokvi->SetupKinematic(effKinEnergy, currentMaterial);
437           lambdaeff = GetTransportMeanFreePath(particle, effKinEnergy);
438           G4double tau = geomStepLength/lambdaeff;
439 
440           if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); } 
441           else               { tPathLength = currentRange; }
442         }
443         zPathLength = geomStepLength;
444       }
445     }
446   }
447   // check of step length
448   // define threshold angle between single and multiple scattering 
449   if(!singleScatteringMode) {
450     cosThetaMin -= ssFactor*tPathLength/lambdaeff; 
451     xtsec = 0.0;
452 
453     // recompute transport cross section - do not change energy
454     // anymore - cannot be applied for big steps
455     if(cosThetaMin > cosTetMaxNuc) {
456       // new computation
457       G4double cross = ComputeTransportXSectionPerVolume(cosThetaMin);
458       //G4cout << "%%%% cross= " << cross << "  xtsec= " << xtsec 
459       //           << " 1-cosTMin= " << 1.0 - cosThetaMin << G4endl;
460       if(cross <= 0.0) {
461         singleScatteringMode = true;
462         tPathLength = zPathLength; 
463         lambdaeff = DBL_MAX;
464         cosThetaMin = 1.0;
465       } else if(xtsec > 0.0) {
466         
467         lambdaeff = 1./cross; 
468         G4double tau = zPathLength*cross;
469         if(tau < numlimit) { 
470           tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0); 
471         } else if(tau < 0.999999) { 
472           tPathLength = -lambdaeff*G4Log(1.0 - tau); 
473         } else { 
474           tPathLength = currentRange;
475         }
476       }
477     } 
478   }
479   tPathLength = std::min(tPathLength, currentRange);
480   /*      
481   G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
482          <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
483   G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
484          << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc 
485          << " e(MeV)= " << preKinEnergy/MeV << "  "  
486          << " SSmode= " << singleScatteringMode << G4endl;
487   */
488   return tPathLength;
489 }
490 
491 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
492 
493 G4ThreeVector& 
494 G4WentzelVIModel::SampleScattering(const G4ThreeVector& oldDirection,
495                                    G4double /*safety*/)
496 {
497   fDisplacement.set(0.0,0.0,0.0);
498   //G4cout << "!##! G4WentzelVIModel::SampleScattering for " 
499   //         << particle->GetParticleName() << G4endl;
500 
501   // ignore scattering for zero step length and energy below the limit
502   if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0) 
503     { return fDisplacement; }
504   
505   G4double invlambda = 0.0;
506   if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
507 
508   // use average kinetic energy over the step
509   G4double cut = (*currentCuts)[currentMaterialIndex];
510   if(fixedCut > 0.0) { cut = fixedCut; }
511   /*  
512   G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
513            << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec 
514           << " xmsc= " <<  tPathLength*invlambda 
515          << " safety= " << safety << G4endl;
516   */
517   // step limit due msc
518   G4int nMscSteps = 1;
519   G4double x0 = tPathLength;
520   G4double z0 = x0*invlambda;
521   G4double prob2 = 0.0;
522 
523   CLHEP::HepRandomEngine* rndmEngine = G4Random::getTheEngine();
524 
525   // large scattering angle case - two step approach
526   if(!singleScatteringMode) {
527     if(useSecondMoment) { 
528       G4double z1 = invlambda*invlambda;
529       G4double z2 = SecondMoment(particle, currentCouple, effKinEnergy);
530       prob2 = (z2 - z1)/(1.5*z1 - z2);
531     }
532     // numerical limitation on step length for 2-step mode
533     if (z0 > 0.05 && z0 < 30.) { 
534       x0 *= 0.5; 
535       z0 *= 0.5;
536       nMscSteps = 2;
537       G4double zzz = G4Exp(-1.0/z0); 
538       z0 += zzz; 
539       prob2 *= (1.0 + zzz);
540     }
541     prob2 /= (1.0 + prob2);
542   } 
543 
544   // step limit due to single scattering
545   G4double x1 = 2*tPathLength;
546   if(0.0 < xtsec) { x1 = -G4Log(rndmEngine->flat())/xtsec; }
547 
548   // no scattering case
549   if(singleScatteringMode && x1 > tPathLength)  
550     { return fDisplacement; }
551 
552   const G4ElementVector* theElementVector = 
553     currentMaterial->GetElementVector();
554   std::size_t nelm = currentMaterial->GetNumberOfElements();
555 
556   // geometry
557   G4double sint, cost, phi;
558   G4ThreeVector temp(0.0,0.0,1.0);
559 
560   // current position and direction relative to the end point
561   // because of magnetic field geometry is computed relatively to the 
562   // end point of the step 
563   G4ThreeVector dir(0.0,0.0,1.0);
564   fDisplacement.set(0.0,0.0,-zPathLength);
565 
566   G4double mscfac = zPathLength/tPathLength;
567 
568   // start a loop 
569   G4double x2 = x0;
570   G4double step, z;
571   G4bool singleScat;
572   /*   
573     G4cout << "Start of the loop x1(mm)= " << x1 << "  x2(mm)= " << x2 
574     << " 1-cost1= " << 1 - cosThetaMin << " SSmode= " << singleScatteringMode 
575            << " xtsec= " << xtsec << " Nst= "  << nMscSteps << G4endl;
576   */
577   do {
578 
579     //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= "<< x2 << G4endl;
580     // single scattering case
581     if(singleScatteringMode && x1 > x2) { 
582       fDisplacement += x2*mscfac*dir;
583       break; 
584     }
585 
586     // what is next single of multiple?
587     if(x1 <= x2) { 
588       step = x1;
589       singleScat = true;
590     } else {
591       step = x2;
592       singleScat = false;
593     }
594 
595     //G4cout << "# step(mm)= "<< step<< "  singlScat= "<< singleScat << G4endl;
596 
597     // new position
598     fDisplacement += step*mscfac*dir;
599 
600     if(singleScat) {
601 
602       // select element
603       std::size_t i = 0;
604       if(nelm > 1) {
605         G4double qsec = rndmEngine->flat()*xtsec;
606         for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
607       }
608       G4double cosTetM = 
609         wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut);
610       //G4cout << "!!! " << cosThetaMin << "  " << cosTetM << "  " 
611       //     << prob[i] << G4endl;
612       temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
613 
614       // direction is changed
615       temp.rotateUz(dir);
616       dir = temp;
617       //G4cout << dir << G4endl;
618 
619       // new proposed step length
620       x2 -= step; 
621       x1  = -G4Log(rndmEngine->flat())/xtsec; 
622 
623     // multiple scattering
624     } else { 
625       --nMscSteps;
626       x1 -= step;
627       x2  = x0;
628 
629       // sample z in interval 0 - 1
630       G4bool isFirst = true;
631       if(prob2 > 0.0 && rndmEngine->flat() < prob2) { isFirst = false; } 
632       do {
633         if(isFirst) { z = -G4Log(rndmEngine->flat()); }
634         else        { z = G4RandGamma::shoot(rndmEngine, 2.0, 2.0); }
635         z *= z0;
636         // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
637       } while(z > 1.0);
638 
639       cost = 1.0 - 2.0*z/*factCM*/;
640       if(cost > 1.0)       { cost = 1.0; }
641       else if(cost < -1.0) { cost =-1.0; }
642       sint = sqrt((1.0 - cost)*(1.0 + cost));
643       phi  = twopi*rndmEngine->flat();
644       G4double vx1 = sint*std::cos(phi);
645       G4double vy1 = sint*std::sin(phi);
646 
647       // lateral displacement  
648       if (latDisplasment) {
649         G4double rms = z0 > 0.0 ? invsqrt12*std::sqrt(2*z0) : 0.0;
650         G4double r  = x0*mscfac;
651         G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(rndmEngine,0.0,1.0));
652         G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(rndmEngine,0.0,1.0));
653         G4double d  = r*r - dx*dx - dy*dy;
654 
655         // change position
656         if(d >= 0.0)  { 
657           temp.set(dx, dy, std::sqrt(d) - r);
658           temp.rotateUz(dir); 
659           fDisplacement += temp;
660         }
661       }
662       // change direction
663       temp.set(vx1,vy1,cost);
664       temp.rotateUz(dir);
665       dir = temp;
666     }
667     // Loop checking, 03-Aug-2015, Vladimir Ivanchenko
668   } while (0 < nMscSteps);
669     
670   dir.rotateUz(oldDirection);
671 
672   //G4cout<<"G4WentzelVIModel sampling is done 1-cost= "<< 1.-dir.z()<<G4endl;
673   // end of sampling -------------------------------
674 
675   fParticleChange->ProposeMomentumDirection(dir);
676 
677   // lateral displacement  
678   fDisplacement.rotateUz(oldDirection);
679 
680   /*
681          G4cout << " r(mm)= " << fDisplacement.mag() 
682                 << " safety= " << safety
683                 << " trueStep(mm)= " << tPathLength
684                 << " geomStep(mm)= " << zPathLength
685                 << " x= " << fDisplacement.x() 
686                 << " y= " << fDisplacement.y() 
687                 << " z= " << fDisplacement.z()
688                 << G4endl;
689   */
690 
691   //G4cout<< "G4WentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl;
692   return fDisplacement;
693 }
694 
695 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
696 
697 G4double G4WentzelVIModel::ComputeTransportXSectionPerVolume(G4double cosTheta)
698 {
699   // prepare recomputation of x-sections
700   const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
701   const G4double* theAtomNumDensityVector = 
702     currentMaterial->GetVecNbOfAtomsPerVolume();
703   G4int nelm = (G4int)currentMaterial->GetNumberOfElements();
704   if(nelm > nelments) {
705     nelments = nelm;
706     xsecn.resize(nelm);
707     prob.resize(nelm);
708   }
709 
710   // check consistency
711   xtsec = 0.0;
712   if(cosTetMaxNuc >= cosTheta) { return 0.0; }
713 
714   G4double cut = (*currentCuts)[currentMaterialIndex];
715   if(fixedCut > 0.0) { cut = fixedCut; }
716 
717   // loop over elements
718   G4double xs = 0.0;
719   for (G4int i=0; i<nelm; ++i) {
720     G4double costm = 
721       wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut);
722     G4double density = theAtomNumDensityVector[i];
723 
724     G4double esec = 0.0;
725     if(costm < cosTheta) {  
726 
727       // recompute the transport x-section
728       if(1.0 > cosTheta) {
729         xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosTheta);
730       }
731       // recompute the total x-section
732       G4double nucsec = wokvi->ComputeNuclearCrossSection(cosTheta, costm);
733       esec = wokvi->ComputeElectronCrossSection(cosTheta, costm);
734       nucsec += esec;
735       if(nucsec > 0.0) { esec /= nucsec; }
736       xtsec += nucsec*density;
737     }
738     xsecn[i] = xtsec;
739     prob[i]  = esec;
740     //G4cout << i << "  xs= " << xs << " xtsec= " << xtsec 
741     //       << " 1-cosTheta= " << 1-cosTheta 
742     //           << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl;
743   }
744   
745   //G4cout << "ComputeXS result:  xsec(1/mm)= " << xs 
746   //         << " txsec(1/mm)= " << xtsec <<G4endl; 
747   return xs;
748 }
749 
750 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
751 
752 G4double G4WentzelVIModel::ComputeSecondMoment(const G4ParticleDefinition* p,
753                  G4double kinEnergy)
754 {
755   G4double xs = 0.0;
756 
757   SetupParticle(p);
758   cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
759 
760   if(cosTetMaxNuc >= 1.0) { return xs; }
761 
762   const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
763   const G4double* theAtomNumDensityVector = 
764     currentMaterial->GetVecNbOfAtomsPerVolume();
765   std::size_t nelm = currentMaterial->GetNumberOfElements();
766 
767   G4double cut = (*currentCuts)[currentMaterialIndex];
768   if(fixedCut > 0.0) { cut = fixedCut; }
769 
770   // loop over elements
771   for (std::size_t i=0; i<nelm; ++i) {
772     G4double costm = 
773       wokvi->SetupTarget((*theElementVector)[i]->GetZasInt(), cut);
774     xs += theAtomNumDensityVector[i]
775       *wokvi->ComputeSecondTransportMoment(costm);
776   }
777   return xs;
778 }
779 
780 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
781 
782 void G4WentzelVIModel::SetSingleScatteringFactor(G4double val)
783 {
784   if(val > 0.05) {
785     ssFactor = val;
786     invssFactor = 1.0/(val - 0.05);
787   }
788 }
789 
790 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
791