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

Last change on this file since 1034 was 1007, checked in by garnier, 15 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.