source: trunk/source/processes/hadronic/models/coherent_elastic/src/G4DiffuseElastic.cc@ 1036

Last change on this file since 1036 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 35.6 KB
Line 
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: G4DiffuseElastic.cc,v 1.20 2008/01/14 10:39:13 vnivanch Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28//
29//
30// Physics model class G4DiffuseElastic
31//
32//
33// G4 Model: optical diffuse elastic scattering with 4-momentum balance
34//
35// 24-May-07 V. Grichine
36//
37
38#include "G4DiffuseElastic.hh"
39#include "G4ParticleTable.hh"
40#include "G4ParticleDefinition.hh"
41#include "G4IonTable.hh"
42
43#include "Randomize.hh"
44#include "G4Integrator.hh"
45#include "globals.hh"
46
47#include "G4Proton.hh"
48#include "G4Neutron.hh"
49#include "G4Deuteron.hh"
50#include "G4Alpha.hh"
51#include "G4PionPlus.hh"
52#include "G4PionMinus.hh"
53
54#include "G4Element.hh"
55#include "G4ElementTable.hh"
56#include "G4PhysicsTable.hh"
57#include "G4PhysicsLogVector.hh"
58#include "G4PhysicsFreeVector.hh"
59
60/////////////////////////////////////////////////////////////////////////
61//
62// Test Constructor. Just to check xsc
63
64
65G4DiffuseElastic::G4DiffuseElastic()
66 : G4HadronicInteraction(), fParticle(0)
67{
68 SetMinEnergy( 0.01*GeV );
69 SetMaxEnergy( 100.*TeV );
70 verboseLevel = 0;
71 lowEnergyRecoilLimit = 100.*keV;
72 lowEnergyLimitQ = 0.0*GeV;
73 lowEnergyLimitHE = 0.0*GeV;
74 lowestEnergyLimit= 0.0*keV;
75 plabLowLimit = 20.0*MeV;
76
77 theProton = G4Proton::Proton();
78 theNeutron = G4Neutron::Neutron();
79 theDeuteron = G4Deuteron::Deuteron();
80 theAlpha = G4Alpha::Alpha();
81 thePionPlus = G4PionPlus::PionPlus();
82 thePionMinus= G4PionMinus::PionMinus();
83
84 fEnergyBin = 200;
85 fAngleBin = 100;
86
87 fEnergyVector = 0;
88 fAngleTable = 0;
89
90 fParticle = 0;
91 fWaveVector = 0.;
92 fAtomicWeight = 0.;
93 fAtomicNumber = 0.;
94 fNuclearRadius = 0.;
95 fBeta = 0.;
96 fZommerfeld = 0.;
97 fAm = 0.;
98 fAddCoulomb = false;
99}
100
101//////////////////////////////////////////////////////////////////////////
102//
103// Constructor with initialisation
104
105G4DiffuseElastic::G4DiffuseElastic(const G4ParticleDefinition* aParticle)
106 : G4HadronicInteraction(), fParticle(aParticle)
107{
108 SetMinEnergy( 0.01*GeV );
109 SetMaxEnergy( 100.*TeV );
110 verboseLevel = 0;
111 lowEnergyRecoilLimit = 100.*keV;
112 lowEnergyLimitQ = 0.0*GeV;
113 lowEnergyLimitHE = 0.0*GeV;
114 lowestEnergyLimit= 0.0*keV;
115 plabLowLimit = 20.0*MeV;
116
117 theProton = G4Proton::Proton();
118 theNeutron = G4Neutron::Neutron();
119 theDeuteron = G4Deuteron::Deuteron();
120 theAlpha = G4Alpha::Alpha();
121 thePionPlus = G4PionPlus::PionPlus();
122 thePionMinus= G4PionMinus::PionMinus();
123
124 fEnergyBin = 200;
125 fAngleBin = 100;
126
127 // fEnergyVector = 0;
128 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
129 fAngleTable = 0;
130
131 fParticle = aParticle;
132 fWaveVector = 0.;
133 fAtomicWeight = 0.;
134 fAtomicNumber = 0.;
135 fNuclearRadius = 0.;
136 fBeta = 0.;
137 fZommerfeld = 0.;
138 fAm = 0.;
139 fAddCoulomb = false;
140 // Initialise();
141}
142
143//////////////////////////////////////////////////////////////////////////////
144//
145// Destructor
146
147G4DiffuseElastic::~G4DiffuseElastic()
148{
149 if(fEnergyVector) delete fEnergyVector;
150
151 if( fAngleTable )
152 {
153 fAngleTable->clearAndDestroy();
154 delete fAngleTable ;
155 }
156}
157
158//////////////////////////////////////////////////////////////////////////////
159//
160// Initialisation for given particle using element table of application
161
162void G4DiffuseElastic::Initialise()
163{
164
165 // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
166
167 const G4ElementTable* theElementTable = G4Element::GetElementTable();
168
169 size_t jEl, numOfEl = G4Element::GetNumberOfElements();
170
171 for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
172 {
173 fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number
174 fAtomicWeight = (*theElementTable)[jEl]->GetN(); // number of nucleons
175 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
176 if(verboseLevel > 0)
177 G4cout<<"G4DiffuseElastic::Initialise() the element: "
178 <<(*theElementTable)[jEl]->GetName()<<G4endl;
179 fElementNumberVector.push_back(fAtomicNumber);
180 fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
181
182 BuildAngleTable();
183 fAngleBank.push_back(fAngleTable);
184 }
185 return;
186}
187
188//////////////////////////////////////////////////////////////////////////////
189//
190// Initialisation for given particle on fly using new element number
191
192void G4DiffuseElastic::InitialiseOnFly(G4double Z, G4double A)
193{
194 fAtomicNumber = Z; // atomic number
195 fAtomicWeight = A; // number of nucleons
196 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
197 if(verboseLevel > 0)
198 G4cout<<"G4DiffuseElastic::Initialise() the element with Z = "
199 <<Z<<"; and A = "<<A<<G4endl;
200 fElementNumberVector.push_back(fAtomicNumber);
201
202 BuildAngleTable();
203 fAngleBank.push_back(fAngleTable);
204 return;
205}
206
207///////////////////////////////////////////////////////////////////////////////
208//
209// Build for given particle and element table of momentum, angle probability.
210// For the moment in lab system.
211
212void G4DiffuseElastic::BuildAngleTable()
213{
214 G4int i, j;
215 G4double partMom, kinE, a=0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
216 G4double theta1, theta2, thetaMax, thetaCoulomb, sum = 0.;
217
218 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
219
220 fAngleTable = new G4PhysicsTable(fEnergyBin);
221
222 for(i = 0; i < fEnergyBin; i++)
223 {
224 kinE = fEnergyVector->GetLowEdgeEnergy(i);
225 partMom = std::sqrt( kinE*(kinE + 2*m1) );
226 fWaveVector = partMom/hbarc;
227
228 thetaMax = 10.174/fWaveVector/fNuclearRadius;
229
230 if (thetaMax > pi) thetaMax = pi;
231
232 thetaCoulomb = 0.2*thetaMax;
233
234 if(z)
235 {
236 a = partMom/m1;
237 fBeta = a/std::sqrt(1+a*a);
238 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
239 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
240 }
241 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin);
242
243 G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.01*thetaMax, thetaMax, fAngleBin );
244
245 for(j = 1; j < fAngleBin; j++)
246 {
247 theta1 = angleBins->GetLowEdgeEnergy(j-1);
248 theta2 = angleBins->GetLowEdgeEnergy(j);
249
250 if(theta2 > thetaCoulomb && z) fAddCoulomb = true;
251
252 sum += integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, theta1,theta2);
253
254 angleVector->PutValue( j-1 , theta2, sum );
255 // G4cout<<"j-1 = "<<j-1<<"; theta2 = "<<theta2<<"; sum = "<<sum<<G4endl;
256 }
257 fAddCoulomb = false;
258
259 fAngleTable->insertAt(i,angleVector);
260
261 // delete[] angleVector;
262 // delete[] angleBins;
263 }
264 return;
265}
266
267////////////////////////////////////////////////////////////////////////////////
268//
269// Model analog of DoIt function
270
271G4HadFinalState*
272G4DiffuseElastic::ApplyYourself( const G4HadProjectile& aTrack,
273 G4Nucleus& targetNucleus )
274{
275 theParticleChange.Clear();
276
277 const G4HadProjectile* aParticle = &aTrack;
278
279 G4double ekin = aParticle->GetKineticEnergy();
280
281 if(ekin <= lowestEnergyLimit)
282 {
283 theParticleChange.SetEnergyChange(ekin);
284 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
285 return &theParticleChange;
286 }
287
288 G4double aTarget = targetNucleus.GetN();
289 G4double zTarget = targetNucleus.GetZ();
290
291 G4double plab = aParticle->GetTotalMomentum();
292
293 if (verboseLevel >1)
294 {
295 G4cout << "G4DiffuseElastic::DoIt: Incident particle plab="
296 << plab/GeV << " GeV/c "
297 << " ekin(MeV) = " << ekin/MeV << " "
298 << aParticle->GetDefinition()->GetParticleName() << G4endl;
299 }
300 // Scattered particle referred to axis of incident particle
301
302 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
303 G4double m1 = theParticle->GetPDGMass();
304
305 G4int Z = static_cast<G4int>(zTarget+0.5);
306 G4int A = static_cast<G4int>(aTarget+0.5);
307 G4int N = A - Z;
308
309 G4int projPDG = theParticle->GetPDGEncoding();
310
311 if (verboseLevel>1)
312 {
313 G4cout << "G4DiffuseElastic for " << theParticle->GetParticleName()
314 << " PDGcode= " << projPDG << " on nucleus Z= " << Z
315 << " A= " << A << " N= " << N
316 << G4endl;
317 }
318 G4ParticleDefinition * theDef = 0;
319
320 if(Z == 1 && A == 1) theDef = theProton;
321 else if (Z == 1 && A == 2) theDef = theDeuteron;
322 else if (Z == 1 && A == 3) theDef = G4Triton::Triton();
323 else if (Z == 2 && A == 3) theDef = G4He3::He3();
324 else if (Z == 2 && A == 4) theDef = theAlpha;
325 else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
326
327 G4double m2 = theDef->GetPDGMass();
328 G4LorentzVector lv1 = aParticle->Get4Momentum();
329 G4LorentzVector lv(0.0,0.0,0.0,m2);
330 lv += lv1;
331
332 G4ThreeVector bst = lv.boostVector();
333 lv1.boost(-bst);
334
335 G4ThreeVector p1 = lv1.vect();
336 G4double ptot = p1.mag();
337 G4double tmax = 4.0*ptot*ptot;
338 G4double t = 0.0;
339
340
341 //
342 // Sample t
343 //
344
345 // t = SampleT( theParticle, ptot, A);
346
347 t = SampleTableT( theParticle, ptot, Z, A); // use initialised table
348
349 // NaN finder
350 if(!(t < 0.0 || t >= 0.0))
351 {
352 if (verboseLevel > 0)
353 {
354 G4cout << "G4DiffuseElastic:WARNING: Z= " << Z << " N= "
355 << N << " pdg= " << projPDG
356 << " mom(GeV)= " << plab/GeV
357 << " S-wave will be sampled"
358 << G4endl;
359 }
360 t = G4UniformRand()*tmax;
361 }
362 if(verboseLevel>1)
363 {
364 G4cout <<" t= " << t << " tmax= " << tmax
365 << " ptot= " << ptot << G4endl;
366 }
367 // Sampling of angles in CM system
368
369 G4double phi = G4UniformRand()*twopi;
370 G4double cost = 1. - 2.0*t/tmax;
371 G4double sint;
372
373 if( cost >= 1.0 )
374 {
375 cost = 1.0;
376 sint = 0.0;
377 }
378 else if( cost <= -1.0)
379 {
380 cost = -1.0;
381 sint = 0.0;
382 }
383 else
384 {
385 sint = std::sqrt((1.0-cost)*(1.0+cost));
386 }
387 if (verboseLevel>1)
388 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
389
390 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
391 v1 *= ptot;
392 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
393
394 nlv1.boost(bst);
395
396 G4double eFinal = nlv1.e() - m1;
397
398 if (verboseLevel > 1)
399 {
400 G4cout << "Scattered: "
401 << nlv1<<" m= " << m1 << " ekin(MeV)= " << eFinal
402 << " Proj: 4-mom " << lv1
403 <<G4endl;
404 }
405 if(eFinal < 0.0)
406 {
407 G4cout << "G4DiffuseElastic WARNING ekin= " << eFinal
408 << " after scattering of "
409 << aParticle->GetDefinition()->GetParticleName()
410 << " p(GeV/c)= " << plab
411 << " on " << theDef->GetParticleName()
412 << G4endl;
413 eFinal = 0.0;
414 nlv1.setE(m1);
415 }
416
417 theParticleChange.SetMomentumChange(nlv1.vect().unit());
418 theParticleChange.SetEnergyChange(eFinal);
419
420 G4LorentzVector nlv0 = lv - nlv1;
421 G4double erec = nlv0.e() - m2;
422
423 if (verboseLevel > 1)
424 {
425 G4cout << "Recoil: "
426 << nlv0<<" m= " << m2 << " ekin(MeV)= " << erec
427 <<G4endl;
428 }
429 if(erec > lowEnergyRecoilLimit)
430 {
431 G4DynamicParticle * aSec = new G4DynamicParticle(theDef, nlv0);
432 theParticleChange.AddSecondary(aSec);
433 } else {
434 if(erec < 0.0) erec = 0.0;
435 theParticleChange.SetLocalEnergyDeposit(erec);
436 }
437
438 return &theParticleChange;
439}
440
441
442////////////////////////////////////////////////////////////////////////////
443//
444// return differential elastic cross section d(sigma)/d(omega)
445
446G4double
447G4DiffuseElastic::GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
448 G4double theta,
449 G4double momentum,
450 G4double A )
451{
452 fParticle = particle;
453 fWaveVector = momentum/hbarc;
454 fAtomicWeight = A;
455 fAddCoulomb = false;
456 fNuclearRadius = CalculateNuclearRad(A);
457
458 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
459
460 return sigma;
461}
462
463////////////////////////////////////////////////////////////////////////////
464//
465// return invariant differential elastic cross section d(sigma)/d(tMand)
466
467G4double
468G4DiffuseElastic::GetInvElasticXsc( const G4ParticleDefinition* particle,
469 G4double tMand,
470 G4double plab,
471 G4double A, G4double Z )
472{
473 G4double m1 = particle->GetPDGMass();
474 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
475
476 G4int iZ = static_cast<G4int>(Z+0.5);
477 G4int iA = static_cast<G4int>(A+0.5);
478 G4ParticleDefinition * theDef = 0;
479
480 if (iZ == 1 && iA == 1) theDef = theProton;
481 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
482 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
483 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
484 else if (iZ == 2 && iA == 4) theDef = theAlpha;
485 else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
486
487 G4double tmass = theDef->GetPDGMass();
488
489 G4LorentzVector lv(0.0,0.0,0.0,tmass);
490 lv += lv1;
491
492 G4ThreeVector bst = lv.boostVector();
493 lv1.boost(-bst);
494
495 G4ThreeVector p1 = lv1.vect();
496 G4double ptot = p1.mag();
497 G4double ptot2 = ptot*ptot;
498 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
499
500 if( cost >= 1.0 ) cost = 1.0;
501 else if( cost <= -1.0) cost = -1.0;
502
503 G4double thetaCMS = std::acos(cost);
504
505 G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
506
507 sigma *= pi/ptot2;
508
509 return sigma;
510}
511
512////////////////////////////////////////////////////////////////////////////
513//
514// return differential elastic cross section d(sigma)/d(omega) with Coulomb
515// correction
516
517G4double
518G4DiffuseElastic::GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
519 G4double theta,
520 G4double momentum,
521 G4double A, G4double Z )
522{
523 fParticle = particle;
524 fWaveVector = momentum/hbarc;
525 fAtomicWeight = A;
526 fAtomicNumber = Z;
527 G4double z = particle->GetPDGCharge();
528 if(z)
529 {
530 fAddCoulomb = true;
531 fBeta = CalculateParticleBeta( particle, momentum);
532 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
533 fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
534 }
535 fNuclearRadius = CalculateNuclearRad(A);
536
537 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
538
539 return sigma;
540}
541
542////////////////////////////////////////////////////////////////////////////
543//
544// return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
545// correction
546
547G4double
548G4DiffuseElastic::GetInvElasticSumXsc( const G4ParticleDefinition* particle,
549 G4double tMand,
550 G4double plab,
551 G4double A, G4double Z )
552{
553 G4double m1 = particle->GetPDGMass();
554 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
555
556 G4int iZ = static_cast<G4int>(Z+0.5);
557 G4int iA = static_cast<G4int>(A+0.5);
558 G4ParticleDefinition * theDef = 0;
559
560 if (iZ == 1 && iA == 1) theDef = theProton;
561 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
562 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
563 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
564 else if (iZ == 2 && iA == 4) theDef = theAlpha;
565 else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
566
567 G4double tmass = theDef->GetPDGMass();
568
569 G4LorentzVector lv(0.0,0.0,0.0,tmass);
570 lv += lv1;
571
572 G4ThreeVector bst = lv.boostVector();
573 lv1.boost(-bst);
574
575 G4ThreeVector p1 = lv1.vect();
576 G4double ptot = p1.mag();
577 G4double ptot2 = ptot*ptot;
578 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
579
580 if( cost >= 1.0 ) cost = 1.0;
581 else if( cost <= -1.0) cost = -1.0;
582
583 G4double thetaCMS = std::acos(cost);
584
585 G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
586
587 sigma *= pi/ptot2;
588
589 return sigma;
590}
591
592////////////////////////////////////////////////////////////////////////////
593//
594// return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
595// correction
596
597G4double
598G4DiffuseElastic::GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
599 G4double tMand,
600 G4double plab,
601 G4double A, G4double Z )
602{
603 G4double m1 = particle->GetPDGMass();
604 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
605
606 G4int iZ = static_cast<G4int>(Z+0.5);
607 G4int iA = static_cast<G4int>(A+0.5);
608 G4ParticleDefinition * theDef = 0;
609
610 if (iZ == 1 && iA == 1) theDef = theProton;
611 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
612 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
613 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
614 else if (iZ == 2 && iA == 4) theDef = theAlpha;
615 else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
616
617 G4double tmass = theDef->GetPDGMass();
618
619 G4LorentzVector lv(0.0,0.0,0.0,tmass);
620 lv += lv1;
621
622 G4ThreeVector bst = lv.boostVector();
623 lv1.boost(-bst);
624
625 G4ThreeVector p1 = lv1.vect();
626 G4double ptot = p1.mag();
627 G4double ptot2 = ptot*ptot;
628 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
629
630 if( cost >= 1.0 ) cost = 1.0;
631 else if( cost <= -1.0) cost = -1.0;
632
633 G4double thetaCMS = std::acos(cost);
634
635 G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
636
637 sigma *= pi/ptot2;
638
639 return sigma;
640}
641
642////////////////////////////////////////////////////////////////////////////
643//
644// return differential elastic probability d(probability)/d(omega)
645
646G4double
647G4DiffuseElastic::GetDiffElasticProb( // G4ParticleDefinition* particle,
648 G4double theta
649 // G4double momentum,
650 // G4double A
651 )
652{
653 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
654 G4double delta, diffuse, gamma;
655 G4double e1, e2, bone, bone2;
656
657 // G4double wavek = momentum/hbarc; // wave vector
658 // G4double r0 = 1.08*fermi;
659 // G4double rad = r0*std::pow(A, 1./3.);
660 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
661 G4double kr2 = kr*kr;
662 G4double krt = kr*theta;
663
664 bzero = BesselJzero(krt);
665 bzero2 = bzero*bzero;
666 bone = BesselJone(krt);
667 bone2 = bone*bone;
668 bonebyarg = BesselOneByArg(krt);
669 bonebyarg2 = bonebyarg*bonebyarg;
670
671 if (fParticle == theProton)
672 {
673 diffuse = 0.63*fermi;
674 gamma = 0.3*fermi;
675 delta = 0.1*fermi*fermi;
676 e1 = 0.3*fermi;
677 e2 = 0.35*fermi;
678 }
679 else // as proton, if were not defined
680 {
681 diffuse = 0.63*fermi;
682 gamma = 0.3*fermi;
683 delta = 0.1*fermi*fermi;
684 e1 = 0.3*fermi;
685 e2 = 0.35*fermi;
686 }
687 G4double lambda = 15.; // 15 ok
688 // G4double kg = fWaveVector*gamma; // wavek*delta;
689 G4double kg = lambda*(1.-std::exp(-fWaveVector*gamma/lambda)); // wavek*delta;
690 G4double kg2 = kg*kg;
691 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
692 // G4double dk2t2 = dk2t*dk2t;
693 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
694 G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
695
696 damp = DampFactor(pikdt);
697 damp2 = damp*damp;
698
699 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
700 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
701
702
703 sigma = kg2;
704 // sigma += dk2t2;
705 sigma *= bzero2;
706 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
707 sigma += kr2*bonebyarg2;
708 sigma *= damp2; // *rad*rad;
709
710 return sigma;
711}
712
713////////////////////////////////////////////////////////////////////////////
714//
715// return differential elastic probability d(probability)/d(omega) with
716// Coulomb correction
717
718G4double
719G4DiffuseElastic::GetDiffElasticSumProb( // G4ParticleDefinition* particle,
720 G4double theta
721 // G4double momentum,
722 // G4double A
723 )
724{
725 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
726 G4double delta, diffuse, gamma;
727 G4double e1, e2, bone, bone2;
728
729 // G4double wavek = momentum/hbarc; // wave vector
730 // G4double r0 = 1.08*fermi;
731 // G4double rad = r0*std::pow(A, 1./3.);
732 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
733 G4double kr2 = kr*kr;
734 G4double krt = kr*theta;
735
736 bzero = BesselJzero(krt);
737 bzero2 = bzero*bzero;
738 bone = BesselJone(krt);
739 bone2 = bone*bone;
740 bonebyarg = BesselOneByArg(krt);
741 bonebyarg2 = bonebyarg*bonebyarg;
742
743 if (fParticle == theProton)
744 {
745 diffuse = 0.63*fermi;
746 // diffuse = 0.6*fermi;
747 gamma = 0.3*fermi;
748 delta = 0.1*fermi*fermi;
749 e1 = 0.3*fermi;
750 e2 = 0.35*fermi;
751 }
752 else // as proton, if were not defined
753 {
754 diffuse = 0.63*fermi;
755 gamma = 0.3*fermi;
756 delta = 0.1*fermi*fermi;
757 e1 = 0.3*fermi;
758 e2 = 0.35*fermi;
759 }
760 G4double lambda = 15.; // 15 ok
761 // G4double kg = fWaveVector*gamma; // wavek*delta;
762 G4double kg = lambda*(1.-std::exp(-fWaveVector*gamma/lambda)); // wavek*delta;
763
764 // G4cout<<"kg = "<<kg<<G4endl;
765
766 if(fAddCoulomb) // add Coulomb correction
767 {
768 G4double sinHalfTheta = std::sin(0.5*theta);
769 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
770
771 kg += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
772 // kg += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
773 }
774
775 G4double kg2 = kg*kg;
776 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
777
778 // G4cout<<"dk2t = "<<dk2t<<G4endl;
779
780 // G4double dk2t2 = dk2t*dk2t;
781
782 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
783 G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
784
785 // G4cout<<"pikdt = "<<pikdt<<G4endl;
786
787 damp = DampFactor(pikdt);
788 damp2 = damp*damp;
789
790 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
791 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
792
793 sigma = kg2;
794 // sigma += dk2t2;
795 sigma *= bzero2;
796 sigma += mode2k2*bone2;
797 sigma += e2dk3t*bzero*bone;
798
799 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
800 sigma += kr2*bonebyarg2; // correction at J1()/()
801
802 sigma *= damp2; // *rad*rad;
803
804 return sigma;
805}
806
807
808////////////////////////////////////////////////////////////////////////////
809//
810// return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega)
811
812G4double
813G4DiffuseElastic::GetIntegrandFunction( G4double theta )
814{
815 G4double result;
816
817 result = 2*pi*std::sin(theta);
818 result *= GetDiffElasticSumProb(theta);
819 return result;
820}
821
822////////////////////////////////////////////////////////////////////////////
823//
824// return integral elastic cross section d(sigma)/d(omega) integrated 0 - theta
825
826G4double
827G4DiffuseElastic::IntegralElasticProb( const G4ParticleDefinition* particle,
828 G4double theta,
829 G4double momentum,
830 G4double A )
831{
832 G4double result;
833 fParticle = particle;
834 fWaveVector = momentum/hbarc;
835 fAtomicWeight = A;
836
837 fNuclearRadius = CalculateNuclearRad(A);
838
839
840 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
841
842 // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
843 result = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
844
845 return result;
846}
847
848////////////////////////////////////////////////////////////////////////////
849//
850// Return inv momentum transfer -t > 0
851
852G4double G4DiffuseElastic::SampleT( const G4ParticleDefinition* aParticle, G4double p, G4double A)
853{
854 G4double theta = SampleThetaCMS( aParticle, p, A); // sample theta in cms
855 G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
856 return t;
857}
858
859////////////////////////////////////////////////////////////////////////////
860//
861// Return inv momentum transfer -t > 0 from initialisation table
862
863G4double G4DiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p,
864 G4double Z, G4double A)
865{
866 G4double theta = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta in cms
867 G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
868 return t;
869}
870
871////////////////////////////////////////////////////////////////////////////
872//
873// Return scattering angle sampled in cms
874
875
876G4double
877G4DiffuseElastic::SampleThetaCMS(const G4ParticleDefinition* particle,
878 G4double momentum, G4double A)
879{
880 G4int i, iMax = 100;
881 G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
882
883 fParticle = particle;
884 fWaveVector = momentum/hbarc;
885 fAtomicWeight = A;
886
887 fNuclearRadius = CalculateNuclearRad(A);
888
889 thetaMax = 10.174/fWaveVector/fNuclearRadius;
890
891 if (thetaMax > pi) thetaMax = pi;
892
893 G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
894
895 // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta );
896 norm = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., thetaMax );
897
898 norm *= G4UniformRand();
899
900 for(i = 1; i <= iMax; i++)
901 {
902 theta1 = (i-1)*thetaMax/iMax;
903 theta2 = i*thetaMax/iMax;
904 sum += integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, theta1, theta2);
905
906 if ( sum >= norm )
907 {
908 result = 0.5*(theta1 + theta2);
909 break;
910 }
911 }
912 if (i > iMax ) result = 0.5*(theta1 + theta2);
913
914 G4double sigma = pi*thetaMax/iMax;
915
916 result += G4RandGauss::shoot(0.,sigma);
917
918 if(result < 0.) result = 0.;
919 if(result > thetaMax) result = thetaMax;
920
921 return result;
922}
923
924////////////////////////////////////////////////////////////////////////////
925//
926// Return scattering angle sampled in cms according to precalculated table.
927
928
929G4double
930G4DiffuseElastic::SampleTableThetaCMS(const G4ParticleDefinition* particle,
931 G4double momentum, G4double Z, G4double A)
932{
933 size_t iElement;
934 G4int iMomentum, iAngle;
935 G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;
936 G4double m1 = particle->GetPDGMass();
937
938 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
939 {
940 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
941 }
942 if ( iElement == fElementNumberVector.size() )
943 {
944 InitialiseOnFly(Z,A);
945 // iElement--;
946
947 // G4cout << "G4DiffuseElastic: Element with atomic number " << Z
948 // << " is not found, return zero angle" << G4endl;
949 // return 0.; // no table for this element
950 }
951 // G4cout<<"iElement = "<<iElement<<G4endl;
952
953 fAngleTable = fAngleBank[iElement];
954
955 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
956
957 for(iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
958 {
959 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
960 }
961 if ( iMomentum == fEnergyBin ) iMomentum--; // kinE is more then theMaxEnergy
962 if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy
963 // G4cout<<"iMomentum = "<<iMomentum<<G4endl;
964
965 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
966 {
967 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
968 // G4cout<<"position = "<<position<<G4endl;
969
970 for(iAngle = 0; iAngle < fAngleBin; iAngle++)
971 {
972 if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
973 }
974 if (iAngle == fAngleBin) iAngle--;
975 // G4cout<<"iAngle = "<<iAngle<<G4endl;
976
977 randAngle = GetScatteringAngle(iMomentum, iAngle, position);
978 // G4cout<<"randAngle = "<<randAngle<<G4endl;
979 }
980 else
981 {
982 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
983 // G4cout<<"position = "<<position<<G4endl;
984
985 for(iAngle = 0; iAngle < fAngleBin; iAngle++)
986 {
987 if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
988 }
989 if (iAngle == fAngleBin) iAngle--;
990 // G4cout<<"iAngle = "<<iAngle<<G4endl;
991
992 theta2 = GetScatteringAngle(iMomentum, iAngle, position);
993 // G4cout<<"theta2 = "<<theta2<<G4endl;
994 E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
995 // G4cout<<"E2 = "<<E2<<G4endl;
996
997 iMomentum--;
998
999 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
1000 // G4cout<<"position = "<<position<<G4endl;
1001
1002 for(iAngle = 0; iAngle < fAngleBin; iAngle++)
1003 {
1004 if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
1005 }
1006 if (iAngle == fAngleBin) iAngle--;
1007
1008 theta1 = GetScatteringAngle(iMomentum, iAngle, position);
1009 // G4cout<<"theta1 = "<<theta1<<G4endl;
1010 E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
1011 // G4cout<<"E1 = "<<E1<<G4endl;
1012
1013 W = 1.0/(E2 - E1);
1014 W1 = (E2 - kinE)*W;
1015 W2 = (kinE - E1)*W;
1016
1017 randAngle = W1*theta1 + W2*theta2;
1018 // G4cout<<"randAngle = "<<randAngle<<G4endl;
1019 }
1020 return randAngle;
1021}
1022
1023/////////////////////////////////////////////////////////////////////////////////
1024//
1025//
1026
1027G4double
1028G4DiffuseElastic:: GetScatteringAngle(G4int iMomentum, G4int iAngle, G4double position)
1029{
1030 G4double x1, x2, y1, y2, randAngle;
1031
1032 if( iAngle == 0 )
1033 {
1034 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1035 }
1036 else
1037 {
1038 if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1039 {
1040 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1041 }
1042 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1043 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1044
1045 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1046 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1047
1048 if ( x1 == x2 ) randAngle = x2;
1049 else
1050 {
1051 if ( y1 == y2 ) randAngle = x1 + (x2 - x1)*G4UniformRand();
1052 else
1053 {
1054 randAngle = x1 + (position - y1)*(x2 - x1)/(y2 - y1);
1055 }
1056 }
1057 }
1058 return randAngle;
1059}
1060
1061
1062
1063////////////////////////////////////////////////////////////////////////////
1064//
1065// Return scattering angle sampled in lab system (target at rest)
1066
1067
1068
1069G4double
1070G4DiffuseElastic::SampleThetaLab( const G4HadProjectile* aParticle,
1071 G4double tmass, G4double A)
1072{
1073 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1074 G4double m1 = theParticle->GetPDGMass();
1075 G4double plab = aParticle->GetTotalMomentum();
1076 G4LorentzVector lv1 = aParticle->Get4Momentum();
1077 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1078 lv += lv1;
1079
1080 G4ThreeVector bst = lv.boostVector();
1081 lv1.boost(-bst);
1082
1083 G4ThreeVector p1 = lv1.vect();
1084 G4double ptot = p1.mag();
1085 G4double tmax = 4.0*ptot*ptot;
1086 G4double t = 0.0;
1087
1088
1089 //
1090 // Sample t
1091 //
1092
1093 t = SampleT( theParticle, ptot, A);
1094
1095 // NaN finder
1096 if(!(t < 0.0 || t >= 0.0))
1097 {
1098 if (verboseLevel > 0)
1099 {
1100 G4cout << "G4DiffuseElastic:WARNING: A = " << A
1101 << " mom(GeV)= " << plab/GeV
1102 << " S-wave will be sampled"
1103 << G4endl;
1104 }
1105 t = G4UniformRand()*tmax;
1106 }
1107 if(verboseLevel>1)
1108 {
1109 G4cout <<" t= " << t << " tmax= " << tmax
1110 << " ptot= " << ptot << G4endl;
1111 }
1112 // Sampling of angles in CM system
1113
1114 G4double phi = G4UniformRand()*twopi;
1115 G4double cost = 1. - 2.0*t/tmax;
1116 G4double sint;
1117
1118 if( cost >= 1.0 )
1119 {
1120 cost = 1.0;
1121 sint = 0.0;
1122 }
1123 else if( cost <= -1.0)
1124 {
1125 cost = -1.0;
1126 sint = 0.0;
1127 }
1128 else
1129 {
1130 sint = std::sqrt((1.0-cost)*(1.0+cost));
1131 }
1132 if (verboseLevel>1)
1133 {
1134 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
1135 }
1136 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1137 v1 *= ptot;
1138 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1139
1140 nlv1.boost(bst);
1141
1142 G4ThreeVector np1 = nlv1.vect();
1143
1144 // G4double theta = std::acos( np1.z()/np1.mag() ); // degree;
1145
1146 G4double theta = np1.theta();
1147
1148 return theta;
1149}
1150
1151////////////////////////////////////////////////////////////////////////////
1152//
1153// Return scattering angle in lab system (target at rest) knowing theta in CMS
1154
1155
1156
1157G4double
1158G4DiffuseElastic::ThetaCMStoThetaLab( const G4DynamicParticle* aParticle,
1159 G4double tmass, G4double thetaCMS)
1160{
1161 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1162 G4double m1 = theParticle->GetPDGMass();
1163 // G4double plab = aParticle->GetTotalMomentum();
1164 G4LorentzVector lv1 = aParticle->Get4Momentum();
1165 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1166
1167 lv += lv1;
1168
1169 G4ThreeVector bst = lv.boostVector();
1170
1171 lv1.boost(-bst);
1172
1173 G4ThreeVector p1 = lv1.vect();
1174 G4double ptot = p1.mag();
1175
1176 G4double phi = G4UniformRand()*twopi;
1177 G4double cost = std::cos(thetaCMS);
1178 G4double sint;
1179
1180 if( cost >= 1.0 )
1181 {
1182 cost = 1.0;
1183 sint = 0.0;
1184 }
1185 else if( cost <= -1.0)
1186 {
1187 cost = -1.0;
1188 sint = 0.0;
1189 }
1190 else
1191 {
1192 sint = std::sqrt((1.0-cost)*(1.0+cost));
1193 }
1194 if (verboseLevel>1)
1195 {
1196 G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
1197 }
1198 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1199 v1 *= ptot;
1200 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1201
1202 nlv1.boost(bst);
1203
1204 G4ThreeVector np1 = nlv1.vect();
1205
1206
1207 G4double thetaLab = np1.theta();
1208
1209 return thetaLab;
1210}
1211////////////////////////////////////////////////////////////////////////////
1212//
1213// Return scattering angle in CMS system (target at rest) knowing theta in Lab
1214
1215
1216
1217G4double
1218G4DiffuseElastic::ThetaLabToThetaCMS( const G4DynamicParticle* aParticle,
1219 G4double tmass, G4double thetaLab)
1220{
1221 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1222 G4double m1 = theParticle->GetPDGMass();
1223 G4double plab = aParticle->GetTotalMomentum();
1224 G4LorentzVector lv1 = aParticle->Get4Momentum();
1225 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1226
1227 lv += lv1;
1228
1229 G4ThreeVector bst = lv.boostVector();
1230
1231 // lv1.boost(-bst);
1232
1233 // G4ThreeVector p1 = lv1.vect();
1234 // G4double ptot = p1.mag();
1235
1236 G4double phi = G4UniformRand()*twopi;
1237 G4double cost = std::cos(thetaLab);
1238 G4double sint;
1239
1240 if( cost >= 1.0 )
1241 {
1242 cost = 1.0;
1243 sint = 0.0;
1244 }
1245 else if( cost <= -1.0)
1246 {
1247 cost = -1.0;
1248 sint = 0.0;
1249 }
1250 else
1251 {
1252 sint = std::sqrt((1.0-cost)*(1.0+cost));
1253 }
1254 if (verboseLevel>1)
1255 {
1256 G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
1257 }
1258 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1259 v1 *= plab;
1260 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1261
1262 nlv1.boost(-bst);
1263
1264 G4ThreeVector np1 = nlv1.vect();
1265
1266
1267 G4double thetaCMS = np1.theta();
1268
1269 return thetaCMS;
1270}
1271
1272//
1273//
1274/////////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.