QwGeant4
QweakSimWentzelVIModel.cc
Go to the documentation of this file.
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 // $Id: G4WentzelVIModel.cc 81865 2014-06-06 11:32:58Z gcosmo $
27 //
28 // -------------------------------------------------------------------
29 //
30 // GEANT4 Class file
31 //
32 //
33 // File name: G4WentzelVIModel
34 //
35 // Author: V.Ivanchenko
36 //
37 // Creation date: 09.04.2008 from G4MuMscModel
38 //
39 // Modifications:
40 // 27-05-2010 V.Ivanchenko added G4WentzelOKandVIxSection class to
41 // compute cross sections and sample scattering angle
42 //
43 //
44 // Class Description:
45 //
46 // Implementation of the model of multiple scattering based on
47 // G.Wentzel, Z. Phys. 40 (1927) 590.
48 // H.W.Lewis, Phys Rev 78 (1950) 526.
49 // J.M. Fernandez-Varea et al., NIM B73 (1993) 447.
50 // L.Urban, CERN-OPEN-2006-077.
51 
52 // -------------------------------------------------------------------
53 //
54 
55 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57 
59 
60 #include "G4PhysicalConstants.hh"
61 #include "G4SystemOfUnits.hh"
62 #include "Randomize.hh"
63 #include "G4ParticleChangeForMSC.hh"
64 #include "G4PhysicsTableHelper.hh"
65 #include "G4ElementVector.hh"
66 #include "G4ProductionCutsTable.hh"
67 #include "G4LossTableManager.hh"
68 #include "G4Log.hh"
69 #include "G4Exp.hh"
70 
71 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
72 
73 using namespace std;
74 
76  G4VMscModel(nam),
77  numlimit(0.1),
78  currentCouple(0),
79  cosThetaMin(1.0),
80  inside(false),
81  singleScatteringMode(false)
82 {
83  invsqrt12 = 1./sqrt(12.);
84  tlimitminfix = 1.e-6*mm;
85  lowEnergyLimit = 1.0*eV;
86  particle = 0;
87  nelments = 5;
88  xsecn.resize(nelments);
89  prob.resize(nelments);
90  theManager = G4LossTableManager::Instance();
91  wokvi = new G4WentzelOKandVIxSection();
92  fixedCut = -1.0;
93 
95  = xtsec = 0;
97  cosThetaMax = cosTetMaxNuc = 1.0;
98 
99  fParticleChange = 0;
100  currentCuts = 0;
101  currentMaterial = 0;
102 }
103 
104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
105 
107 {
108  delete wokvi;
109 }
110 
111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
112 
113 void QweakSimWentzelVIModel::Initialise(const G4ParticleDefinition* p,
114  const G4DataVector& cuts)
115 {
116  // reset parameters
117  SetupParticle(p);
118  currentRange = 0.0;
119 
120  cosThetaMax = cos(PolarAngleLimit());
121  //G4cout << "QweakSimWentzelVIModel::Initialise " << p->GetParticleName() << G4endl;
122  wokvi->Initialise(p, cosThetaMax);
123  /*
124  G4cout << "QweakSimWentzelVIModel: " << particle->GetParticleName()
125  << " 1-cos(ThetaLimit)= " << 1 - cosThetaMax
126  << G4endl;
127  */
128  currentCuts = &cuts;
129 
130  // set values of some data members
131  fParticleChange = GetParticleChangeForMSC(p);
132 }
133 
134 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
135 
137  const G4ParticleDefinition* p,
138  G4double kinEnergy,
139  G4double Z, G4double,
140  G4double cutEnergy, G4double)
141 {
142  G4double cross = 0.0;
143  if(p != particle) { SetupParticle(p); }
144  if(kinEnergy < lowEnergyLimit) { return cross; }
145  if(!CurrentCouple()) {
146  G4Exception("QweakSimWentzelVIModel::ComputeCrossSectionPerAtom", "em0011",
147  FatalException, " G4MaterialCutsCouple is not defined");
148  return 0.0;
149  }
150  DefineMaterial(CurrentCouple());
151  cosTetMaxNuc = wokvi->SetupKinematic(kinEnergy, currentMaterial);
152  if(cosTetMaxNuc < 1.0) {
153  G4double cut = cutEnergy;
154  if(fixedCut > 0.0) { cut = fixedCut; }
155  cosTetMaxNuc = wokvi->SetupTarget(G4lrint(Z), cut);
156  cross = wokvi->ComputeTransportCrossSectionPerAtom(cosTetMaxNuc);
157  /*
158  if(p->GetParticleName() == "e-")
159  G4cout << "QweakSimWentzelVIModel::CS: Z= " << G4int(Z) << " e(MeV)= " << kinEnergy
160  << " 1-cosN= " << 1 - cosTetMaxNuc << " cross(bn)= " << cross/barn
161  << " " << particle->GetParticleName() << G4endl;
162  */
163  }
164  return cross;
165 }
166 
167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
168 
170 {
171  SetupParticle(track->GetDynamicParticle()->GetDefinition());
172  inside = false;
173 }
174 
175 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
176 
178  const G4Track& track,
179  G4double& currentMinimalStep)
180 {
181  G4double tlimit = currentMinimalStep;
182  const G4DynamicParticle* dp = track.GetDynamicParticle();
183  G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
184  G4StepStatus stepStatus = sp->GetStepStatus();
185  singleScatteringMode = false;
186 
187  // FIXME
188  ePolarized=false;
189  debugPrint=false;
190  if(strcmp(track.GetParticleDefinition()->GetParticleName().data() , "e-") == 0)
191  if(strcmp(track.GetMaterial()->GetName(),"PBA") == 0){
192  if(track.GetPolarization().getR() >= 0.1) debugPrint=true;
193  if(sqrt(pow(track.GetPolarization().getX(),2)+pow(track.GetPolarization().getY(),2))>0.01){
194  ePolarized=true;
195  polarization=track.GetPolarization();
196  eEnergy=track.GetTotalEnergy();
197  }
198  }
199  debugPrint=false;
200  // FIXME
201 
202  //G4cout << "QweakSimWentzelVIModel::ComputeTruePathLengthLimit stepStatus= "
203  // << stepStatus << " " << track.GetDefinition()->GetParticleName()
204  // << G4endl;
205 
206  // initialisation for each step, lambda may be computed from scratch
207  preKinEnergy = dp->GetKineticEnergy();
208  DefineMaterial(track.GetMaterialCutsCouple());
209  lambdaeff = GetTransportMeanFreePath(particle,preKinEnergy);
211  cosTetMaxNuc = wokvi->SetupKinematic(preKinEnergy, currentMaterial);
212 
213  //G4cout << "lambdaeff= " << lambdaeff << " Range= " << currentRange
214  // << " tlimit= " << tlimit << " 1-cost= " << 1 - cosTetMaxNuc << G4endl;
215 
216  // extra check for abnormal situation
217  // this check needed to run MSC with eIoni and eBrem inactivated
218  if(tlimit > currentRange) { tlimit = currentRange; }
219 
220  // stop here if small range particle
221  if(inside || tlimit < tlimitminfix) {
222  return ConvertTrueToGeom(tlimit, currentMinimalStep);
223  }
224 
225  // pre step
226  G4double presafety = sp->GetSafety();
227  // far from geometry boundary
228  if(currentRange < presafety) {
229  inside = true;
230  return ConvertTrueToGeom(tlimit, currentMinimalStep);
231  }
232 
233  // compute presafety again if presafety <= 0 and no boundary
234  // i.e. when it is needed for optimization purposes
235  if(stepStatus != fGeomBoundary && presafety < tlimitminfix) {
236  presafety = ComputeSafety(sp->GetPosition(), tlimit);
237  if(currentRange < presafety) {
238  inside = true;
239  return ConvertTrueToGeom(tlimit, currentMinimalStep);
240  }
241  }
242  /*
243  G4cout << "e(MeV)= " << preKinEnergy/MeV
244  << " " << particle->GetParticleName()
245  << " CurLimit(mm)= " << tlimit/mm <<" safety(mm)= " << presafety/mm
246  << " R(mm)= " <<currentRange/mm
247  << " L0(mm^-1)= " << lambdaeff*mm
248  <<G4endl;
249  */
250 
251  // natural limit for high energy
252  G4double rlimit = std::max(facrange*currentRange,
253  0.7*(1.0 - cosTetMaxNuc)*lambdaeff);
254 
255  // low-energy e-
256  if(cosThetaMax > cosTetMaxNuc) {
257  rlimit = std::min(rlimit, facsafety*presafety);
258  }
259 
260  // cut correction
261  G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
262  //G4cout << "rcut= " << rcut << " rlimit= " << rlimit << " presafety= "
263  // << presafety << " 1-cosThetaMax= " <<1-cosThetaMax
264  //<< " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc << G4endl;
265  if(rcut > rlimit) { rlimit = std::min(rlimit, rcut*sqrt(rlimit/rcut)); }
266 
267  if(rlimit < tlimit) { tlimit = rlimit; }
268 
269  tlimit = std::max(tlimit, tlimitminfix);
270 
271  // step limit in infinite media
272  tlimit = std::min(tlimit, 50*currentMaterial->GetRadlen()/facgeom);
273 
274  //compute geomlimit and force few steps within a volume
275  if (steppingAlgorithm == fUseDistanceToBoundary
276  && stepStatus == fGeomBoundary) {
277 
278  G4double geomlimit = ComputeGeomLimit(track, presafety, currentRange);
279  tlimit = std::min(tlimit, geomlimit/facgeom);
280  }
281 
282  /*
283  G4cout << particle->GetParticleName() << " e= " << preKinEnergy
284  << " L0= " << lambdaeff << " R= " << currentRange
285  << "tlimit= " << tlimit
286  << " currentMinimalStep= " << currentMinimalStep << G4endl;
287  */
288  return ConvertTrueToGeom(tlimit, currentMinimalStep);
289 }
290 
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
292 
294 {
295  tPathLength = truelength;
297 
298  if(lambdaeff > 0.0 && lambdaeff != DBL_MAX) {
299  G4double tau = tPathLength/lambdaeff;
300  //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
301  // << " Leff= " << lambdaeff << " tau= " << tau << G4endl;
302  // small step
303  if(tau < numlimit) {
304  zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
305 
306  // medium step
307  } else {
308  G4double e1 = 0.0;
309  if(currentRange > tPathLength) {
311  }
312  e1 = 0.5*(e1 + preKinEnergy);
313  cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
314  lambdaeff = GetTransportMeanFreePath(particle,e1);
315  zPathLength = lambdaeff*(1.0 - G4Exp(-tPathLength/lambdaeff));
316  }
317  } else { lambdaeff = DBL_MAX; }
318  //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
319  return zPathLength;
320 }
321 
322 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
323 
324 G4double QweakSimWentzelVIModel::ComputeTrueStepLength(G4double geomStepLength)
325 {
326  // initialisation of single scattering x-section
327  xtsec = 0.0;
329  /*
330  G4cout << "ComputeTrueStepLength: Step= " << geomStepLength
331  << " Lambda= " << lambdaeff
332  << " 1-cosThetaMaxNuc= " << 1 - cosTetMaxNuc << G4endl;
333  */
334  // pathalogical case
335  if(lambdaeff == DBL_MAX) {
336  singleScatteringMode = true;
337  zPathLength = geomStepLength;
338  tPathLength = geomStepLength;
339  cosThetaMin = 1.0;
340 
341  // normal case
342  } else {
343 
344  // small step use only single scattering
345  static const G4double singleScatLimit = 1.0e-7;
346  if(geomStepLength < lambdaeff*singleScatLimit*(1.0 - cosTetMaxNuc)) {
347  singleScatteringMode = true;
348  cosThetaMin = 1.0;
349  zPathLength = geomStepLength;
350  tPathLength = geomStepLength;
351 
352  // step defined by transportation
353  } else if(geomStepLength != zPathLength) {
354 
355  // step defined by transportation
356  zPathLength = geomStepLength;
357  G4double tau = geomStepLength/lambdaeff;
358  tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
359 
360  // energy correction for a big step
361  if(tau > numlimit) {
362  G4double e1 = 0.0;
363  if(currentRange > tPathLength) {
365  }
366  e1 = 0.5*(e1 + preKinEnergy);
367  cosTetMaxNuc = wokvi->SetupKinematic(e1, currentMaterial);
368  lambdaeff = GetTransportMeanFreePath(particle,e1);
369  tau = zPathLength/lambdaeff;
370 
371  if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); }
372  else { tPathLength = currentRange; }
373  }
374  }
375  }
376 
377  // check of step length
378  // define threshold angle between single and multiple scattering
380 
381  // recompute transport cross section - do not change energy
382  // anymore - cannot be applied for big steps
383  if(cosThetaMin > cosTetMaxNuc) {
384 
385  // new computation
386  G4double cross = ComputeXSectionPerVolume();
387  //G4cout << "%%%% cross= " << cross << " xtsec= " << xtsec << G4endl;
388  if(cross <= 0.0) {
389  singleScatteringMode = true;
391  lambdaeff = DBL_MAX;
392  cosThetaMin = 1.0;
393  } else if(xtsec > 0.0) {
394 
395  lambdaeff = 1./cross;
396  G4double tau = zPathLength*cross;
397  if(tau < numlimit) {
398  tPathLength = zPathLength*(1.0 + 0.5*tau + tau*tau/3.0);
399  }
400  else if(tau < 0.999999) { tPathLength = -lambdaeff*G4Log(1.0 - tau); }
401  else { tPathLength = currentRange; }
402 
404  }
405  }
406 
407  /*
408  G4cout <<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength
409  <<" Leff(mm)= "<<lambdaeff/mm<<" sig0(1/mm)= " << xtsec <<G4endl;
410  G4cout << particle->GetParticleName() << " 1-cosThetaMin= " << 1-cosThetaMin
411  << " 1-cosTetMaxNuc= " << 1-cosTetMaxNuc
412  << " e(MeV)= " << preKinEnergy/MeV << " "
413  << singleScatteringMode << G4endl;
414  */
415  return tPathLength;
416 }
417 
418 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
419 
420 G4ThreeVector&
421 QweakSimWentzelVIModel::SampleScattering(const G4ThreeVector& oldDirection,
422  G4double safety)
423 {
424  fDisplacement.set(0.0,0.0,0.0);
425  //G4cout << "!##! QweakSimWentzelVIModel::SampleScattering for "
426  // << particle->GetParticleName() << G4endl;
427 
428  // ignore scattering for zero step length and energy below the limit
429  if(preKinEnergy < lowEnergyLimit || tPathLength <= 0.0)
430  { return fDisplacement; }
431 
432  G4double invlambda = 0.0;
433  if(lambdaeff < DBL_MAX) { invlambda = 0.5/lambdaeff; }
434 
435  // use average kinetic energy over the step
436  G4double cut = (*currentCuts)[currentMaterialIndex];
437  if(fixedCut > 0.0) { cut = fixedCut; }
438  /*
439  G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy/MeV
440  << " Leff= " << lambdaeff <<" sig0(1/mm)= " << xtsec
441  << " xmsc= " << tPathLength*invlambda
442  << " safety= " << safety << G4endl;
443  */
444 
445  // step limit due msc
446  G4int nMscSteps = 1;
447  G4double x0 = tPathLength;
448  G4double z0 = x0*invlambda;
449  G4double zzz = 0.0;
450 
451  // large scattering angle case - two step approach
452 
453  if(!singleScatteringMode) {
454  static const G4double zzmin = 0.05;
455  if(z0 > zzmin && safety > tlimitminfix) {
456  x0 *= 0.5;
457  z0 *= 0.5;
458  nMscSteps = 2;
459  }
460  if(z0 > zzmin) { zzz = G4Exp(-1.0/z0); }
461  }
462 
463  // step limit due to single scattering
464  G4double x1 = 2*tPathLength;
465  if(0.0 < xtsec) { x1 = -G4Log(G4UniformRand())/xtsec; }
466 
467  // no scattering case
468  if(singleScatteringMode && x1 > tPathLength)
469  { return fDisplacement; }
470 
471  const G4ElementVector* theElementVector =
472  currentMaterial->GetElementVector();
473  G4int nelm = currentMaterial->GetNumberOfElements();
474 
475  // geometry
476  G4double sint, cost, phi;
477  G4ThreeVector temp(0.0,0.0,1.0);
478 
479  // current position and direction relative to the end point
480  // because of magnetic field geometry is computed relatively to the
481  // end point of the step
482  G4ThreeVector dir(0.0,0.0,1.0);
483  fDisplacement.set(0.0,0.0,-zPathLength);
484 
485  G4double mscfac = zPathLength/tPathLength;
486 
487  // start a loop
488  G4double x2 = x0;
489  G4double step, z;
490  G4bool singleScat;
491  /*
492  G4cout << "Start of the loop x1(mm)= " << x1 << " x2(mm)= " << x2
493  << " 1-cost1= " << 1 - cosThetaMin << " " << singleScatteringMode
494  << " xtsec= " << xtsec << " " << nMscSteps << G4endl;
495  */
496  do {
497 
498  //G4cout << "# x1(mm)= "<< x1<< " x2(mm)= "<< x2 << G4endl;
499  // single scattering case
500  if(singleScatteringMode && x1 > x2) { break; }
501 
502  // what is next single of multiple?
503  if(x1 <= x2) {
504  step = x1;
505  singleScat = true;
506  } else {
507  step = x2;
508  singleScat = false;
509  }
510 
511  //G4cout << "# step(mm)= "<< step<< " singlScat= "<< singleScat << G4endl;
512 
513  // new position
514  fDisplacement += step*mscfac*dir;
515 
516  if(singleScat) {
517 
518  // select element
519  G4int i = 0;
520  if(nelm > 1) {
521  G4double qsec = G4UniformRand()*xtsec;
522  for (; i<nelm; ++i) { if(xsecn[i] >= qsec) { break; } }
523  }
524  G4double cosTetM =
525  wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
526  //G4cout << "!!! " << cosThetaMin << " " << cosTetM << " "
527  // << prob[i] << G4endl;
528  temp = wokvi->SampleSingleScattering(cosThetaMin, cosTetM, prob[i]);
529 
530  //FIXME
531  cost=cos(temp.getTheta());
532  sint = sqrt((1.0 - cost)*(1.0 + cost));
533  phi=temp.getPhi();
534  if(ePolarized){
535  if(debugPrint) G4cout<<" ~~sc~~ "<<nMscSteps<<G4endl;
536  G4double _prob=G4UniformRand();
537  G4double _amplitude=1.0/eEnergy * sint *
538  sqrt(pow(polarization.getX(),2)+pow(polarization.getY(),2));//scale by transvers polarization
539  if(_amplitude > 1 ) _amplitude=1;
540  if( _prob < _amplitude * sin(phi-pi) )
541  phi-=pi;
542  phi+= polarization.getPhi() - oldDirection.getPhi();
543  if(phi<0) phi+=twopi;
544  else if(phi>twopi) phi=fmod(phi,twopi);
545  G4double vx1 = sint*cos(phi);
546  G4double vy1 = sint*sin(phi);
547  temp.set(vx1,vy1,cost);
548  }
549  //FIXME
550 
551  // direction is changed
552  temp.rotateUz(dir);
553  dir = temp;
554 
555  // new proposed step length
556  x2 -= step;
557  x1 = -G4Log(G4UniformRand())/xtsec;
558 
559  // multiple scattering
560  } else {
561  --nMscSteps;
562  x1 -= step;
563  x2 = x0;
564 
565  // sample z in interval 0 - 1
566  do {
567  z = -z0*G4Log(1.0 - (1.0 - zzz)*G4UniformRand());
568  } while(z > 1.0);
569  cost = 1.0 - 2.0*z/*factCM*/;
570  if(cost > 1.0) { cost = 1.0; }
571  else if(cost < -1.0) { cost =-1.0; }
572  sint = sqrt((1.0 - cost)*(1.0 + cost));
573  phi = twopi*G4UniformRand();
574  //FIXME
575  if(ePolarized){
576  if(debugPrint) G4cout<<" ~~msc~~ "<<nMscSteps<<G4endl;
577  G4double _prob=G4UniformRand();
578  G4double _amplitude=1.0/eEnergy * sint*
579  sqrt(pow(polarization.getX(),2)+pow(polarization.getY(),2));//scale by transvers polarization
580  if(_amplitude > 1 ) _amplitude=1;
581  if( _prob < _amplitude * sin(phi-pi) )
582  phi-=pi;
583  phi+= polarization.getPhi() - oldDirection.getPhi();
584  if(phi<0) phi+=twopi;
585  else if(phi>twopi) phi=fmod(phi,twopi);
586  }
587  //FIXME
588  G4double vx1 = sint*cos(phi);
589  G4double vy1 = sint*sin(phi);
590 
591  // lateral displacement
592  if (latDisplasment && safety > tlimitminfix) {
593  G4double rms = invsqrt12*sqrt(2*z0);
594  G4double r = x0*mscfac;
595  G4double dx = r*(0.5*vx1 + rms*G4RandGauss::shoot(0.0,1.0));
596  G4double dy = r*(0.5*vy1 + rms*G4RandGauss::shoot(0.0,1.0));
597  G4double dz;
598  G4double d = r*r - dx*dx - dy*dy;
599 
600  // change position
601  if(d >= 0.0) {
602  dz = sqrt(d) - r;
603  temp.set(dx,dy,dz);
604  temp.rotateUz(dir);
605  fDisplacement += temp;
606  }
607  }
608  // change direction
609  temp.set(vx1,vy1,cost);
610  temp.rotateUz(dir);
611  dir = temp;
612  }
613  } while (0 < nMscSteps);
614 
615  dir.rotateUz(oldDirection);
616 
617  //FIXME
618  if(debugPrint){
619  G4cout<<" WentzelMS cth, th, phi old.angle(new)" << cost << " " << acos(cost) << " " << phi << " " <<oldDirection.angle(dir) << G4endl;
620  G4cout<<" WentzelMS old dir: R th phi "<<oldDirection.getR()<<" "<<oldDirection.getTheta()<<" "<<oldDirection.getPhi()<<G4endl;
621  G4cout<<" WentzelMS new dir: R th phi "<<dir.getR()<<" "<<dir.getTheta()<<" "<<dir.getPhi()<<G4endl;
622  }
623  //FIXME
624 
625  //G4cout<<"QweakSimWentzelVIModel sampling is done 1-cost= "<< 1.-dir.z()<<G4endl;
626  // end of sampling -------------------------------
627 
628  fParticleChange->ProposeMomentumDirection(dir);
629 
630  // lateral displacement
631  fDisplacement.rotateUz(oldDirection);
632 
633  /*
634  G4cout << " r(mm)= " << fDisplacement.mag()
635  << " safety= " << safety
636  << " trueStep(mm)= " << tPathLength
637  << " geomStep(mm)= " << zPathLength
638  << " x= " << fDisplacement.x()
639  << " y= " << fDisplacement.y()
640  << " z= " << fDisplacement.z()
641  << G4endl;
642  */
643 
644  //G4cout<< "QweakSimWentzelVIModel::SampleScattering end NewDir= " << dir<< G4endl;
645  return fDisplacement;
646 }
647 
648 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
649 
651 {
652  // prepare recomputation of x-sections
653  const G4ElementVector* theElementVector = currentMaterial->GetElementVector();
654  const G4double* theAtomNumDensityVector =
655  currentMaterial->GetVecNbOfAtomsPerVolume();
656  G4int nelm = currentMaterial->GetNumberOfElements();
657  if(nelm > nelments) {
658  nelments = nelm;
659  xsecn.resize(nelm);
660  prob.resize(nelm);
661  }
662  G4double cut = (*currentCuts)[currentMaterialIndex];
663  if(fixedCut > 0.0) { cut = fixedCut; }
664  // cosTetMaxNuc = wokvi->GetCosThetaNuc();
665 
666  // check consistency
667  xtsec = 0.0;
668  if(cosTetMaxNuc > cosThetaMin) { return 0.0; }
669 
670  // loop over elements
671  G4double xs = 0.0;
672  for (G4int i=0; i<nelm; ++i) {
673  G4double costm =
674  wokvi->SetupTarget(G4lrint((*theElementVector)[i]->GetZ()), cut);
675  G4double density = theAtomNumDensityVector[i];
676 
677  G4double esec = 0.0;
678  if(costm < cosThetaMin) {
679 
680  // recompute the transport x-section
681  if(1.0 > cosThetaMin) {
682  xs += density*wokvi->ComputeTransportCrossSectionPerAtom(cosThetaMin);
683  }
684  // recompute the total x-section
685  G4double nucsec = wokvi->ComputeNuclearCrossSection(cosThetaMin, costm);
686  esec = wokvi->ComputeElectronCrossSection(cosThetaMin, costm);
687  nucsec += esec;
688  if(nucsec > 0.0) { esec /= nucsec; }
689  xtsec += nucsec*density;
690  }
691  xsecn[i] = xtsec;
692  prob[i] = esec;
693  //G4cout << i << " xs= " << xs << " xtsec= " << xtsec
694  // << " 1-cosThetaMin= " << 1-cosThetaMin
695  // << " 1-cosTetMaxNuc2= " <<1-cosTetMaxNuc2<< G4endl;
696  }
697 
698  //G4cout << "ComputeXS result: xsec(1/mm)= " << xs
699  // << " txsec(1/mm)= " << xtsec <<G4endl;
700  return xs;
701 }
702 
703 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double KineticEnergy, G4double AtomicNumber, G4double AtomicWeight=0., G4double cut=DBL_MAX, G4double emax=DBL_MAX)
std::vector< G4double > prob
virtual G4ThreeVector & SampleScattering(const G4ThreeVector &, G4double safety)
virtual G4double ComputeGeomPathLength(G4double truePathLength)
void SetupParticle(const G4ParticleDefinition *)
virtual G4double ComputeTrueStepLength(G4double geomStepLength)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
std::vector< G4double > xsecn
const G4Material * currentMaterial
const G4DataVector * currentCuts
G4LossTableManager * theManager
const G4MaterialCutsCouple * currentCouple
void DefineMaterial(const G4MaterialCutsCouple *)
const G4ParticleDefinition * particle
G4ParticleChangeForMSC * fParticleChange
G4WentzelOKandVIxSection * wokvi
QweakSimWentzelVIModel(const G4String &nam="WentzelVIUni")
virtual G4double ComputeTruePathLengthLimit(const G4Track &track, G4double &currentMinimalStep)