source: trunk/source/processes/electromagnetic/standard/src/G4WentzelVIModel.cc@ 973

Last change on this file since 973 was 968, checked in by garnier, 17 years ago

fichier ajoutes
File size: 24.1 KB
RevLine 
[968]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,v 1.17 2009/02/19 19:17:15 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
28//
29// -------------------------------------------------------------------
30//
31// GEANT4 Class file
32//
33//
34// File name: G4WentzelVIModel
35//
36// Author: V.Ivanchenko
37//
38// Creation date: 09.04.2008 from G4MuMscModel
39//
40// Modifications:
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
55//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
56//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
57
58#include "G4WentzelVIModel.hh"
59#include "Randomize.hh"
60#include "G4LossTableManager.hh"
61#include "G4ParticleChangeForMSC.hh"
62//#include "G4TransportationManager.hh"
63//#include "G4SafetyHelper.hh"
64#include "G4PhysicsTableHelper.hh"
65#include "G4ElementVector.hh"
66#include "G4ProductionCutsTable.hh"
67#include "G4PhysicsLogVector.hh"
68#include "G4Electron.hh"
69#include "G4Positron.hh"
70#include "G4Proton.hh"
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
73
74using namespace std;
75
76G4WentzelVIModel::G4WentzelVIModel(const G4String& nam) :
77 G4VMscModel(nam),
78 theLambdaTable(0),
79 theLambda2Table(0),
80 numlimit(0.2),
81 nbins(60),
82 nwarnings(0),
83 nwarnlimit(50),
84 currentCouple(0),
85 cosThetaMin(1.0),
86 q2Limit(TeV*TeV),
87 alpha2(fine_structure_const*fine_structure_const),
88 isInitialized(false),
89 inside(false)
90{
91 invsqrt12 = 1./sqrt(12.);
92 tlimitminfix = 1.e-6*mm;
93 theManager = G4LossTableManager::Instance();
94 fNistManager = G4NistManager::Instance();
95 theElectron = G4Electron::Electron();
96 thePositron = G4Positron::Positron();
97 theProton = G4Proton::Proton();
98 a0 = alpha2*electron_mass_c2*electron_mass_c2/(0.885*0.885);
99 G4double p0 = electron_mass_c2*classic_electr_radius;
100 coeff = twopi*p0*p0;
101 constn = 6.937e-6/(MeV*MeV);
102 tkin = targetZ = mom2 = DBL_MIN;
103 ecut = etag = DBL_MAX;
104 particle = 0;
105 nelments = 5;
106 xsecn.resize(nelments);
107 prob.resize(nelments);
108 for(size_t j=0; j<100; j++) {
109 FF[j] = 0.0;
110 }
111}
112
113//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
114
115G4WentzelVIModel::~G4WentzelVIModel()
116{}
117
118//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
119
120void G4WentzelVIModel::Initialise(const G4ParticleDefinition* p,
121 const G4DataVector& cuts)
122{
123 // reset parameters
124 SetupParticle(p);
125 tkin = targetZ = mom2 = DBL_MIN;
126 ecut = etag = DBL_MAX;
127 currentRange = 0.0;
128 cosThetaMax = cos(PolarAngleLimit());
129 currentCuts = &cuts;
130
131 // set values of some data members
132 if(!isInitialized) {
133 isInitialized = true;
134
135 if (pParticleChange)
136 fParticleChange = reinterpret_cast<G4ParticleChangeForMSC*>(pParticleChange);
137 else
138 fParticleChange = new G4ParticleChangeForMSC();
139
140 InitialiseSafetyHelper();
141 }
142}
143
144//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
145
146G4double G4WentzelVIModel::ComputeCrossSectionPerAtom(
147 const G4ParticleDefinition* p,
148 G4double kinEnergy,
149 G4double Z, G4double,
150 G4double cutEnergy, G4double)
151{
152 SetupParticle(p);
153 G4double ekin = std::max(lowEnergyLimit, kinEnergy);
154 SetupKinematic(ekin, cutEnergy);
155 SetupTarget(Z, ekin);
156 G4double xsec = ComputeTransportXSectionPerVolume();
157 /*
158 G4cout << "CS: e= " << tkin << " cosEl= " << cosTetMaxElec2
159 << " cosN= " << cosTetMaxNuc2 << " xsec(bn)= " << xsec/barn
160 << " " << particle->GetParticleName() << G4endl;
161 */
162 return xsec;
163}
164
165//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
166
167G4double G4WentzelVIModel::ComputeTransportXSectionPerVolume()
168{
169 G4double xSection = 0.0;
170 G4double x, y, x1, x2, x3, x4;
171
172 // scattering off electrons
173 if(cosTetMaxElec2 < 1.0) {
174 x = (1.0 - cosTetMaxElec2)/screenZ;
175 if(x < numlimit) y = 0.5*x*x*(1.0 - 1.3333333*x + 1.5*x*x);
176 else y = log(1.0 + x) - x/(1.0 + x);
177 if(y < 0.0) {
178 nwarnings++;
179 if(nwarnings < nwarnlimit /*&& y < -1.e-10*/) {
180 G4cout << "Electron scattering <0 for L1 " << y
181 << " e(MeV)= " << tkin << " p(MeV/c)= " << sqrt(mom2)
182 << " Z= " << targetZ << " "
183 << particle->GetParticleName() << G4endl;
184 G4cout << " z= " << 1.0-cosTetMaxElec2 << " screenZ= " << screenZ
185 << " x= " << x << G4endl;
186 }
187 y = 0.0;
188 }
189 xSection += y/targetZ;
190 }
191 /*
192 G4cout << "G4WentzelVIModel:XS per A " << " Z= " << Z << " e(MeV)= " << kinEnergy/MeV
193 << " cut(MeV)= " << ecut/MeV
194 << " zmaxE= " << (1.0 - cosTetMaxElec)/screenZ
195 << " zmaxN= " << (1.0 - cosTetMsxNuc)/screenZ << G4endl;
196 */
197
198 // scattering off nucleus
199 if(cosTetMaxNuc2 < 1.0) {
200 x = 1.0 - cosTetMaxNuc2;
201 x1 = screenZ*formfactA;
202 x2 = 1.0/(1.0 - x1);
203 x3 = x/screenZ;
204 x4 = formfactA*x;
205 // low-energy limit
206 if(x3 < numlimit && x1 < numlimit) {
207 y = 0.5*x3*x3*x2*x2*x2*(1.0 - 1.333333*x3 + 1.5*x3*x3 - 1.5*x1
208 + 3.0*x1*x1 + 2.666666*x3*x1);
209 // high energy limit
210 } else if(1.0 < x1) {
211 x4 = x1*(1.0 + x3);
212 y = x3*(1.0 + 0.5*x3 - (2.0 - x1)*(1.0 + x3 + x3*x3/3.0)/x4)/(x4*x4);
213 // middle energy
214 } else {
215 y = ((1.0 + x1)*x2*log((1. + x3)/(1. + x4))
216 - x3/(1. + x3) - x4/(1. + x4))*x2*x2;
217 }
218 if(y < 0.0) {
219 nwarnings++;
220 if(nwarnings < nwarnlimit /*&& y < -1.e-10*/) {
221 G4cout << "Nuclear scattering <0 for L1 " << y
222 << " e(MeV)= " << tkin << " Z= " << targetZ << " "
223 << particle->GetParticleName() << G4endl;
224 G4cout << " formfactA= " << formfactA << " screenZ= " << screenZ
225 << " x= " << " x1= " << x1 << " x2= " << x2
226 << " x3= " << x3 << " x4= " << x4 <<G4endl;
227 }
228 y = 0.0;
229 }
230 xSection += y;
231 }
232 xSection *= (coeff*targetZ*targetZ*chargeSquare*invbeta2/mom2);
233 // G4cout << " XStotal= " << xSection/barn << " screenZ= " << screenZ
234 // << " formF= " << formfactA << " for " << p->GetParticleName() << G4endl;
235 return xSection;
236}
237
238//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
239
240G4double G4WentzelVIModel::ComputeTruePathLengthLimit(
241 const G4Track& track,
242 G4PhysicsTable* theTable,
243 G4double currentMinimalStep)
244{
245 G4double tlimit = currentMinimalStep;
246 const G4DynamicParticle* dp = track.GetDynamicParticle();
247 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
248 G4StepStatus stepStatus = sp->GetStepStatus();
249
250 // initialisation for 1st step
251 if(stepStatus == fUndefined) {
252 inside = false;
253 SetupParticle(dp->GetDefinition());
254 theLambdaTable = theTable;
255 }
256
257 // initialisation for each step, lambda may be computed from scratch
258 preKinEnergy = dp->GetKineticEnergy();
259 DefineMaterial(track.GetMaterialCutsCouple());
260 lambda0 = GetLambda(preKinEnergy);
261 currentRange =
262 theManager->GetRangeFromRestricteDEDX(particle,preKinEnergy,currentCouple);
263
264 // extra check for abnormal situation
265 // this check needed to run MSC with eIoni and eBrem inactivated
266 if(tlimit > currentRange) tlimit = currentRange;
267
268 // stop here if small range particle
269 if(inside) return tlimit;
270
271 // pre step
272 G4double presafety = sp->GetSafety();
273
274 // compute presafety again if presafety <= 0 and no boundary
275 // i.e. when it is needed for optimization purposes
276 if(stepStatus != fGeomBoundary && presafety < tlimitminfix)
277 presafety = ComputeSafety(sp->GetPosition(), tlimit);
278 /*
279 G4cout << "G4WentzelVIModel::ComputeTruePathLengthLimit tlimit= "
280 <<tlimit<<" safety= " << presafety
281 << " range= " <<currentRange<<G4endl;
282 */
283 // far from geometry boundary
284 if(currentRange < presafety) {
285 inside = true;
286
287 // limit mean scattering angle
288 } else {
289 G4double rlimit = facrange*lambda0;
290 G4double rcut = currentCouple->GetProductionCuts()->GetProductionCut(1);
291 if(rcut > rlimit) rlimit = std::pow(2.0*rcut*rcut*lambda0,0.33333333);
292 if(rlimit < tlimit) tlimit = rlimit;
293 }
294 /*
295 G4cout << particle->GetParticleName() << " e= " << preKinEnergy
296 << " L0= " << lambda0 << " R= " << currentRange
297 << "tlimit= " << tlimit
298 << " currentMinimalStep= " << currentMinimalStep << G4endl;
299 */
300 return tlimit;
301}
302
303//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
304
305G4double G4WentzelVIModel::ComputeGeomPathLength(G4double truelength)
306{
307 tPathLength = truelength;
308 zPathLength = tPathLength;
309 lambdaeff = lambda0;
310
311 if(lambda0 > 0.0) {
312 G4double tau = tPathLength/lambda0;
313 //G4cout << "ComputeGeomPathLength: tLength= " << tPathLength
314 // << " lambda0= " << lambda0 << " tau= " << tau << G4endl;
315 // small step
316 if(tau < numlimit) {
317 zPathLength *= (1.0 - 0.5*tau + tau*tau/6.0);
318
319 // medium step
320 } else {
321 // zPathLength = lambda0*(1.0 - exp(-tPathLength/lambda0));
322 G4double e1 = 0.0;
323 if(currentRange > tPathLength) {
324 e1 = theManager->GetEnergy(particle,
325 currentRange-tPathLength,
326 currentCouple);
327 }
328 lambdaeff = GetLambda(0.5*(e1 + preKinEnergy));
329 zPathLength = lambdaeff*(1.0 - exp(-tPathLength/lambdaeff));
330 }
331 }
332 //G4cout<<"Comp.geom: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
333 return zPathLength;
334}
335
336//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
337
338G4double G4WentzelVIModel::ComputeTrueStepLength(G4double geomStepLength)
339{
340 // step defined other than transportation
341 if(geomStepLength == zPathLength) return tPathLength;
342
343 // step defined by transportation
344 tPathLength = geomStepLength;
345 zPathLength = geomStepLength;
346 G4double tau = zPathLength/lambdaeff;
347 tPathLength *= (1.0 + 0.5*tau + tau*tau/3.0);
348
349 if(tau > numlimit) {
350 G4double e1 = 0.0;
351 if(currentRange > tPathLength) {
352 e1 = theManager->GetEnergy(particle,
353 currentRange-tPathLength,
354 currentCouple);
355 }
356 lambdaeff = GetLambda(0.5*(e1 + preKinEnergy));
357 tau = zPathLength/lambdaeff;
358
359 if(tau < 0.999999) tPathLength = -lambdaeff*log(1.0 - tau);
360 else tPathLength = currentRange;
361
362 if(tPathLength < zPathLength) tPathLength = zPathLength;
363 }
364 if(tPathLength > currentRange) tPathLength = currentRange;
365 //G4cout<<"Comp.true: zLength= "<<zPathLength<<" tLength= "<<tPathLength<<G4endl;
366 return tPathLength;
367}
368
369//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
370
371void G4WentzelVIModel::SampleScattering(const G4DynamicParticle* dynParticle,
372 G4double safety)
373{
374 //G4cout << "!##! G4WentzelVIModel::SampleScattering for "
375 // << particle->GetParticleName() << G4endl;
376 G4double kinEnergy = dynParticle->GetKineticEnergy();
377 if(kinEnergy <= DBL_MIN || tPathLength <= DBL_MIN) return;
378
379 G4double ekin = preKinEnergy;
380 if(ekin - kinEnergy > ekin*dtrl) {
381 ekin = 0.5*(preKinEnergy + kinEnergy);
382 lambdaeff = GetLambda(ekin);
383 }
384
385 G4double x1 = 0.5*tPathLength/lambdaeff;
386 G4double cut= (*currentCuts)[currentMaterialIndex];
387 /*
388 G4cout <<"SampleScat: E0(MeV)= "<< preKinEnergy<<" Eeff(MeV)= "<<ekin/MeV
389 << " L0= " << lambda0 << " Leff= " << lambdaeff
390 << " x1= " << x1 << " safety= " << safety << G4endl;
391 */
392
393 G4double xsec = 0.0;
394 G4bool largeAng = false;
395
396 // large scattering angle case
397 if(x1 > 0.5) {
398 x1 *= 0.5;
399 largeAng = true;
400
401 // normal case
402 } else {
403
404 // define threshold angle as 2 sigma of central value
405 cosThetaMin = 1.0 - 3.0*x1;
406
407 // for low-energy e-,e+ no limit
408 ekin = std::max(ekin, lowEnergyLimit);
409 SetupKinematic(ekin, cut);
410
411 // recompute transport cross section
412 if(cosThetaMin > cosTetMaxNuc) {
413
414 xsec = ComputeXSectionPerVolume();
415
416 if(xtsec > DBL_MIN) x1 = 0.5*tPathLength*xtsec;
417 else x1 = 0.0;
418
419 /*
420 G4cout << "cosTetMaxNuc= " << cosTetMaxNuc
421 << " cosThetaMin= " << cosThetaMin
422 << " cosThetaMax= " << cosThetaMax
423 << " cosTetMaxElec2= " << cosTetMaxElec2 << G4endl;
424 G4cout << "Recomputed xsec(1/mm)= " << xsec << " x1= " << x1 << G4endl;
425 */
426 }
427 }
428
429 // result of central part sampling
430 G4double z;
431 do {
432 z = -x1*log(G4UniformRand());
433 } while (z > 1.0);
434
435 // cost is sampled ------------------------------
436 G4double cost = 1.0 - 2.0*z;
437 if(cost < -1.0) cost = -1.0;
438 else if(cost > 1.0) cost = 1.0;
439 G4double sint = sqrt((1.0 - cost)*(1.0 + cost));
440
441 G4double phi = twopi*G4UniformRand();
442
443 G4double dirx = sint*cos(phi);
444 G4double diry = sint*sin(phi);
445
446 //G4cout << "G4WentzelVIModel: step(mm)= " << tPathLength/mm
447 // << " sint= " << sint << " cost= " << cost<< G4endl;
448
449 G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
450 G4ThreeVector newDirection(dirx,diry,cost);
451 G4ThreeVector temp(0.0,0.0,1.0);
452 G4ThreeVector pos(0.0,0.0,-zPathLength);
453 G4ThreeVector dir(0.0,0.0,1.0);
454 G4bool isscat = false;
455
456 // sample MSC scattering for large angle
457 // extra central scattering for holf step
458 if(largeAng) {
459 isscat = true;
460 pos.setZ(-0.5*zPathLength);
461 do {
462 z = -x1*log(G4UniformRand());
463 } while (z > 1.0);
464 cost = 1.0 - 2.0*z;
465 if(std::abs(cost) > 1.0) cost = 1.0;
466
467 sint = sqrt((1.0 - cost)*(1.0 + cost));
468 phi = twopi*G4UniformRand();
469
470 // position and direction for secondary scattering
471 dir.set(sint*cos(phi),sint*sin(phi),cost);
472 pos += 0.5*dir*zPathLength;
473 x1 *= 2.0;
474 }
475
476 // sample Reserford scattering for large angle
477 if(xsec > DBL_MIN) {
478 G4double t = tPathLength;
479 G4int nelm = currentMaterial->GetNumberOfElements();
480 const G4ElementVector* theElementVector =
481 currentMaterial->GetElementVector();
482 do{
483 G4double x = -log(G4UniformRand())/xsec;
484 pos += dir*(zPathLength*std::min(x,t)/tPathLength);
485 t -= x;
486 if(t > 0.0) {
487 G4double zz1 = 1.0;
488 G4double qsec = G4UniformRand()*xsec;
489
490 // scattering off nucleus
491 G4int i = 0;
492 if(nelm > 1) {
493 for (; i<nelm; i++) {if(xsecn[i] >= qsec) break;}
494 if(i >= nelm) i = nelm - 1;
495 }
496 SetupTarget((*theElementVector)[i]->GetZ(), tkin);
497 G4double formf = formfactA;
498 G4double costm = cosTetMaxNuc2;
499 if(prob[i] > 0.0) {
500 if(G4UniformRand() <= prob[i]) {
501 formf = 0.0;
502 costm = cosTetMaxElec2;
503 }
504 }
505 if(cosThetaMin > costm) {
506
507 G4double w1 = 1. - cosThetaMin + screenZ;
508 G4double w2 = 1. - costm + screenZ;
509 G4double w3 = cosThetaMin - costm;
510 G4double grej, zz;
511 do {
512 zz = w1*w2/(w1 + G4UniformRand()*w3) - screenZ;
513 grej = 1.0/(1.0 + formf*zz);
514 } while ( G4UniformRand() > grej*grej );
515 if(zz < 0.0) zz = 0.0;
516 else if(zz > 2.0) zz = 2.0;
517 zz1 = 1.0 - zz;
518 }
519 if(zz1 < 1.0) {
520 isscat = true;
521 //G4cout << "Reserford zz1= " << zz1 << " t= " << t << G4endl;
522 sint = sqrt((1.0 - zz1)*(1.0 + zz1));
523 //G4cout << "sint= " << sint << G4endl;
524 phi = twopi*G4UniformRand();
525 G4double vx1 = sint*cos(phi);
526 G4double vy1 = sint*sin(phi);
527 temp.set(vx1,vy1,zz1);
528 temp.rotateUz(dir);
529 dir = temp;
530 }
531 }
532 } while (t > 0.0);
533 }
534 if(isscat) newDirection.rotateUz(dir);
535 newDirection.rotateUz(oldDirection);
536
537 //G4cout << "G4WentzelVIModel sampling of scattering is done" << G4endl;
538 // end of sampling -------------------------------
539
540 fParticleChange->ProposeMomentumDirection(newDirection);
541
542 if (latDisplasment && safety > tlimitminfix) {
543 G4double rms = invsqrt12*sqrt(2.0*x1);
544 G4double dx = zPathLength*(0.5*dirx + rms*G4RandGauss::shoot(0.0,1.0));
545 G4double dy = zPathLength*(0.5*diry + rms*G4RandGauss::shoot(0.0,1.0));
546 G4double dz;
547 G4double d = (dx*dx + dy*dy)/(zPathLength*zPathLength);
548 if(d < numlimit) dz = -0.5*zPathLength*d*(1.0 + 0.25*d);
549 else if(d < 1.0) dz = -zPathLength*(1.0 - sqrt(1.0 - d));
550 else {
551 dx = dy = dz = 0.0;
552 }
553
554 temp.set(dx,dy,dz);
555 if(isscat) temp.rotateUz(dir);
556 pos += temp;
557
558 pos.rotateUz(oldDirection);
559
560 G4double r = pos.mag();
561
562 /*
563 G4cout << " r(mm)= " << r << " safety= " << safety
564 << " trueStep(mm)= " << tPathLength
565 << " geomStep(mm)= " << zPathLength
566 << G4endl;
567 */
568
569 if(r > tlimitminfix) {
570 pos /= r;
571 ComputeDisplacement(fParticleChange, pos, r, safety);
572 }
573 }
574 //G4cout << "G4WentzelVIModel::SampleScattering end" << G4endl;
575}
576
577//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
578
579G4double G4WentzelVIModel::ComputeXSectionPerVolume()
580{
581 const G4ElementVector* theElementVector =
582 currentMaterial->GetElementVector();
583 const G4double* theAtomNumDensityVector =
584 currentMaterial->GetVecNbOfAtomsPerVolume();
585 G4int nelm = currentMaterial->GetNumberOfElements();
586 if(nelm > nelments) {
587 nelments = nelm;
588 xsecn.resize(nelments);
589 prob.resize(nelments);
590 }
591
592 xtsec = 0.0;
593 G4double xs = 0.0;
594
595 G4double fac = coeff*chargeSquare*invbeta2/mom2;
596
597 for (G4int i=0; i<nelm; i++) {
598 SetupTarget((*theElementVector)[i]->GetZ(), tkin);
599 G4double density = theAtomNumDensityVector[i];
600 G4double cosnm = cosTetMaxNuc2;
601 G4double cosem = cosTetMaxElec2;
602
603 // recompute the angular limit
604 cosTetMaxNuc2 = std::max(cosnm,cosThetaMin);
605 cosTetMaxElec2 = std::max(cosem,cosThetaMin);
606 xtsec += ComputeTransportXSectionPerVolume()*density;
607 // return limit back
608 cosTetMaxElec2 = cosem;
609 cosTetMaxNuc2 = cosnm;
610
611 G4double esec = 0.0;
612 G4double nsec = 0.0;
613 G4double x1 = 1.0 - cosThetaMin + screenZ;
614 G4double f = fac*targetZ*density;
615
616 // scattering off electrons
617 if(cosThetaMin > cosem) {
618 esec = f*(cosThetaMin - cosem)/(x1*(1.0 - cosem + screenZ));
619 }
620
621 // scattering off nucleaus
622 if(cosThetaMin > cosnm) {
623
624 // Reserford part
625 G4double s = screenZ*formfactA;
626 G4double z1 = 1.0 - cosnm + screenZ;
627 G4double d = (1.0 - s)/formfactA;
628
629 // check numerical limit
630 if(d < numlimit*x1) {
631 G4double x2 = x1*x1;
632 G4double z2 = z1*z1;
633 nsec = (1.0/(x1*x2) - 1.0/(z1*z2) - d*1.5*(1.0/(x2*x2) - 1.0/(z2*z2)))/
634 (3.0*formfactA*formfactA);
635 } else {
636 G4double x2 = x1 + d;
637 G4double z2 = z1 + d;
638 nsec = (1.0 + 2.0*s)*((cosThetaMin - cosnm)*(1.0/(x1*z1) + 1.0/(x2*z2)) -
639 2.0*log(z1*x2/(z2*x1))/d);
640 }
641 nsec *= f*targetZ;
642 }
643 nsec += esec;
644 if(nsec > 0.0) esec /= nsec;
645 xs += nsec;
646 xsecn[i] = xs;
647 prob[i] = esec;
648 //G4cout << i << " xs= " << xs << " cosThetaMin= " << cosThetaMin
649 // << " costm= " << costm << G4endl;
650 }
651
652 //G4cout << "ComputeXS result: xsec(1/mm)= " << xs
653 //<< " txsec(1/mm)= " << xtsec <<G4endl;
654 return xs;
655}
656
657//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
658
659/*
660G4double G4MuMscModel::ComputeXSectionPerVolume()
661{
662 const G4ElementVector* theElementVector =
663 currentMaterial->GetElementVector();
664 const G4double* theAtomNumDensityVector =
665 currentMaterial->GetVecNbOfAtomsPerVolume();
666 size_t nelm = currentMaterial->GetNumberOfElements();
667
668 xsece1 = 0.0;
669 xsece2 = 0.0;
670 xsecn2 = 0.0;
671 zcorr = 0.0;
672
673 G4double fac = coeff*chargeSquare*invbeta2/mom2;
674
675 for (size_t i=0; i<nelm; i++) {
676 const G4Element* elm = (*theElementVector)[i];
677 G4double Z = elm->GetZ();
678 SetupTarget(Z, tkin);
679 G4double den = fac*theAtomNumDensityVector[i]*Z;
680
681 G4double x = 1.0 - cosThetaMin;
682 G4double x1 = x + screenZ;
683 G4double x2 = 1.0/(x1*x1);
684 G4double x3 = 1.0 + x*formfactA;
685
686 //G4cout << "x= " << x << " den= " << den << " cosE= " << cosTetMaxElec << G4endl;
687 //G4cout << "cosThtaMin= " << cosThetaMin << G4endl;
688 //G4cout << "cosTetMaxNuc= " << cosTetMaxNuc << " q2Limit= " << q2Limit << G4endl;
689
690 // scattering off electrons
691 if(cosTetMaxElec < cosThetaMin) {
692
693 // flat part
694 G4double s = den*x2*x;
695 xsece1 += s;
696 zcorr += 0.5*x*s;
697
698 // Reserford part
699 G4double z1 = 1.0 - cosTetMaxElec + screenZ;
700 G4double z2 = (cosThetaMin - cosTetMaxElec)/x1;
701 if(z2 < 0.2) s = z2*(x - 0.5*z2*(x - screenZ))/x1;
702 else s = log(1.0 + z2) - screenZ*z2/z1;
703 xsece2 += den*z2/z1;
704 zcorr += den*s;
705 }
706 den *= Z;
707
708 //G4cout << "Z= " << Z<< " cosL= " << cosTetMaxNuc << " cosMin= " << cosThetaMin << G4endl;
709 // scattering off nucleaus
710 if(cosTetMaxNuc < cosThetaMin) {
711
712 // flat part
713 G4double s = den*x2*x/(x3*x3);
714 xsece1 += s;
715 zcorr += 0.5*x*s;
716
717 // Reserford part
718 s = screenZ*formfactA;
719 G4double w = 1.0 + 2.0*s;
720 G4double z1 = 1.0 - cosTetMaxNuc + screenZ;
721 G4double d = (1.0 - s)/formfactA;
722 G4double x4 = x1 + d;
723 G4double z4 = z1 + d;
724 G4double t1 = 1.0/(x1*z1);
725 G4double t4 = 1.0/(x4*z4);
726 G4double w1 = cosThetaMin - cosTetMaxNuc;
727 G4double w2 = log(z1*x4/(x1*z4));
728
729 den *= w;
730 xsecn2 += den*(w1*(t1 + t4) - 2.0*w2/d);
731 zcorr += den*(w*w2 - w1*(screenZ*t1 + t4/formfactA));
732 }
733 xsece[i] = xsece2;
734 xsecn[i] = xsecn2;
735 // G4cout << i << " xsece2= " << xsece2 << " xsecn2= " << xsecn2 << G4endl;
736 }
737 G4double xsec = xsece1 + xsece2 + xsecn2;
738
739 //G4cout << "xsece1= " << xsece1 << " xsece2= " << xsece2
740 //<< " xsecn2= " << xsecn2
741 // << " zsec= " << zcorr*0.5*tPathLength << G4endl;
742 zcorr *= 0.5*tPathLength;
743
744 return xsec;
745}
746*/
747
748//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
749
750void G4WentzelVIModel::ComputeMaxElectronScattering(G4double cutEnergy)
751{
752 ecut = cutEnergy;
753 G4double tmax = tkin;
754 cosTetMaxElec = 1.0;
755 if(mass > MeV) {
756 G4double ratio = electron_mass_c2/mass;
757 G4double tau = tkin/mass;
758 tmax = 2.0*electron_mass_c2*tau*(tau + 2.)/
759 (1.0 + 2.0*ratio*(tau + 1.0) + ratio*ratio);
760 cosTetMaxElec = 1.0 - std::min(cutEnergy, tmax)*electron_mass_c2/mom2;
761 } else {
762
763 if(particle == theElectron) tmax *= 0.5;
764 G4double t = std::min(cutEnergy, tmax);
765 G4double mom21 = t*(t + 2.0*electron_mass_c2);
766 G4double t1 = tkin - t;
767 //G4cout <<"tkin=" <<tkin<<" tmax= "<<tmax<<" t= "
768 //<<t<< " t1= "<<t1<<" cut= "<<ecut<<G4endl;
769 if(t1 > 0.0) {
770 G4double mom22 = t1*(t1 + 2.0*mass);
771 G4double ctm = (mom2 + mom22 - mom21)*0.5/sqrt(mom2*mom22);
772 if(ctm < 1.0) cosTetMaxElec = ctm;
773 }
774 }
775 if(cosTetMaxElec < cosTetMaxNuc) cosTetMaxElec = cosTetMaxNuc;
776}
777
778//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
779
780void G4WentzelVIModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
781 const G4MaterialCutsCouple*,
782 const G4DynamicParticle*,
783 G4double,
784 G4double)
785{}
786
787//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.