source: trunk/source/processes/hadronic/models/coherent_elastic/src/G4NuclNuclDiffuseElastic.cc @ 1228

Last change on this file since 1228 was 1228, checked in by garnier, 14 years ago

update geant4.9.3 tag
File size: 44.2 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: G4NuclNuclDiffuseElastic.cc,v 1.2 2009/04/10 13:22:25 grichine Exp $
27// GEANT4 tag $Name: geant4-09-03 $
28//
29//
30// Physics model class G4NuclNuclDiffuseElastic
31//
32//
33// G4 Model: optical diffuse elastic scattering with 4-momentum balance
34//                         
35// 24-May-07 V. Grichine
36//
37
38#include "G4NuclNuclDiffuseElastic.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
65G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseElastic() 
66  : G4HadronicInteraction(), fParticle(0)
67{
68  SetMinEnergy( 0.01*GeV );
69  SetMaxEnergy( 1.*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  = 200;
86
87  fEnergyVector =  new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
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  fProfileDelta  = 1.;
101  fProfileAlpha   = 0.5;
102
103}
104
105//////////////////////////////////////////////////////////////////////////
106//
107// Constructor with initialisation
108
109G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseElastic(const G4ParticleDefinition* aParticle) 
110  : G4HadronicInteraction(), fParticle(aParticle)
111{
112  SetMinEnergy( 0.01*GeV );
113  SetMaxEnergy( 1.*TeV );
114  verboseLevel = 0;
115  lowEnergyRecoilLimit = 100.*keV; 
116  lowEnergyLimitQ  = 0.0*GeV; 
117  lowEnergyLimitHE = 0.0*GeV; 
118  lowestEnergyLimit= 0.0*keV; 
119  plabLowLimit     = 20.0*MeV;
120
121  theProton   = G4Proton::Proton();
122  theNeutron  = G4Neutron::Neutron();
123  theDeuteron = G4Deuteron::Deuteron();
124  theAlpha    = G4Alpha::Alpha();
125  thePionPlus = G4PionPlus::PionPlus();
126  thePionMinus= G4PionMinus::PionMinus();
127
128  fEnergyBin = 200; // 200; // 100;
129  fAngleBin  = 400; // 200; // 100;
130
131  // fEnergyVector = 0;
132  fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
133  fAngleTable = 0;
134
135  fParticle = aParticle;
136  fWaveVector = 0.;
137  fAtomicWeight = 0.;
138  fAtomicNumber = 0.;
139  fNuclearRadius = 0.;
140  fBeta = 0.;
141  fZommerfeld = 0.;
142  fAm = 0.;
143  fAddCoulomb = false;
144  // Initialise();
145}
146
147//////////////////////////////////////////////////////////////////////////////
148//
149// Destructor
150
151G4NuclNuclDiffuseElastic::~G4NuclNuclDiffuseElastic()
152{
153  if(fEnergyVector) delete fEnergyVector;
154
155  if( fAngleTable )
156  {
157      fAngleTable->clearAndDestroy();
158      delete fAngleTable ;
159  }
160}
161
162//////////////////////////////////////////////////////////////////////////////
163//
164// Initialisation for given particle using element table of application
165
166void G4NuclNuclDiffuseElastic::Initialise() 
167{
168
169  // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
170
171  const G4ElementTable* theElementTable = G4Element::GetElementTable();
172
173  size_t jEl, numOfEl = G4Element::GetNumberOfElements();
174
175  for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
176  {
177    fAtomicNumber = (*theElementTable)[jEl]->GetZ();     // atomic number
178    fAtomicWeight = (*theElementTable)[jEl]->GetN();     // number of nucleons
179    fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
180
181    if(verboseLevel > 0) 
182    {   
183      G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element: "
184            <<(*theElementTable)[jEl]->GetName()<<G4endl;
185    }
186    fElementNumberVector.push_back(fAtomicNumber);
187    fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
188
189    BuildAngleTable();
190    fAngleBank.push_back(fAngleTable);
191  } 
192  return;
193}
194
195////////////////////////////////////////////////////////////////////////////////
196//
197// Model analog of DoIt function
198
199G4HadFinalState* 
200G4NuclNuclDiffuseElastic::ApplyYourself( const G4HadProjectile& aTrack, 
201                                       G4Nucleus& targetNucleus )
202{
203  theParticleChange.Clear();
204
205  const G4HadProjectile* aParticle = &aTrack;
206
207  G4double ekin = aParticle->GetKineticEnergy();
208
209  if(ekin <= lowestEnergyLimit) 
210  {
211    theParticleChange.SetEnergyChange(ekin);
212    theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
213    return &theParticleChange;
214  }
215
216  G4double aTarget = targetNucleus.GetN();
217  G4double zTarget = targetNucleus.GetZ();
218
219  G4double plab = aParticle->GetTotalMomentum();
220
221  if (verboseLevel >1)
222  { 
223    G4cout << "G4NuclNuclDiffuseElastic::DoIt: Incident particle plab=" 
224           << plab/GeV << " GeV/c " 
225           << " ekin(MeV) = " << ekin/MeV << "  " 
226           << aParticle->GetDefinition()->GetParticleName() << G4endl;
227  }
228  // Scattered particle referred to axis of incident particle
229
230  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
231  G4double m1 = theParticle->GetPDGMass();
232
233  G4int Z = static_cast<G4int>(zTarget+0.5);
234  G4int A = static_cast<G4int>(aTarget+0.5);
235  G4int N = A - Z;
236
237  G4int projPDG = theParticle->GetPDGEncoding();
238
239  if (verboseLevel>1) 
240  {
241    G4cout << "G4NuclNuclDiffuseElastic for " << theParticle->GetParticleName()
242           << " PDGcode= " << projPDG << " on nucleus Z= " << Z
243           << " A= " << A << " N= " << N
244           << G4endl;
245  }
246  G4ParticleDefinition * theDef = 0;
247
248  if(Z == 1 && A == 1)       theDef = theProton;
249  else if (Z == 1 && A == 2) theDef = theDeuteron;
250  else if (Z == 1 && A == 3) theDef = G4Triton::Triton();
251  else if (Z == 2 && A == 3) theDef = G4He3::He3();
252  else if (Z == 2 && A == 4) theDef = theAlpha;
253  else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
254 
255  G4double m2 = theDef->GetPDGMass();
256  G4LorentzVector lv1 = aParticle->Get4Momentum();
257  G4LorentzVector lv(0.0,0.0,0.0,m2);   
258  lv += lv1;
259
260  G4ThreeVector bst = lv.boostVector();
261  lv1.boost(-bst);
262
263  G4ThreeVector p1 = lv1.vect();
264  G4double ptot = p1.mag();
265  G4double tmax = 4.0*ptot*ptot;
266  G4double t = 0.0;
267
268
269  //
270  // Sample t
271  //
272 
273  // t = SampleT( theParticle, ptot, A);
274
275  t = SampleTableT( theParticle, ptot, Z, A); // use initialised table
276
277  // NaN finder
278  if(!(t < 0.0 || t >= 0.0)) 
279  {
280    if (verboseLevel > 0) 
281    {
282      G4cout << "G4NuclNuclDiffuseElastic:WARNING: Z= " << Z << " N= " 
283             << N << " pdg= " <<  projPDG
284             << " mom(GeV)= " << plab/GeV
285              << " S-wave will be sampled" 
286             << G4endl; 
287    }
288    t = G4UniformRand()*tmax; 
289  }
290  if(verboseLevel>1)
291  {
292    G4cout <<" t= " << t << " tmax= " << tmax
293           << " ptot= " << ptot << G4endl;
294  }
295  // Sampling of angles in CM system
296
297  G4double phi  = G4UniformRand()*twopi;
298  G4double cost = 1. - 2.0*t/tmax;
299  G4double sint;
300
301  if( cost >= 1.0 ) 
302  {
303    cost = 1.0;
304    sint = 0.0;
305  }
306  else if( cost <= -1.0) 
307  {
308    cost = -1.0;
309    sint =  0.0;
310  }
311  else 
312  {
313    sint = std::sqrt((1.0-cost)*(1.0+cost));
314  }   
315  if (verboseLevel>1) 
316    G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
317
318  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
319  v1 *= ptot;
320  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
321
322  nlv1.boost(bst); 
323
324  G4double eFinal = nlv1.e() - m1;
325
326  if (verboseLevel > 1)
327  { 
328    G4cout << "Scattered: "
329           << nlv1<<" m= " << m1 << " ekin(MeV)= " << eFinal
330           << " Proj: 4-mom " << lv1
331           <<G4endl;
332  }
333  if(eFinal < 0.0) 
334  {
335    G4cout << "G4NuclNuclDiffuseElastic WARNING ekin= " << eFinal
336           << " after scattering of " 
337           << aParticle->GetDefinition()->GetParticleName()
338           << " p(GeV/c)= " << plab
339           << " on " << theDef->GetParticleName()
340           << G4endl;
341    eFinal = 0.0;
342    nlv1.setE(m1);
343  }
344
345  theParticleChange.SetMomentumChange(nlv1.vect().unit());
346  theParticleChange.SetEnergyChange(eFinal);
347 
348  G4LorentzVector nlv0 = lv - nlv1;
349  G4double erec =  nlv0.e() - m2;
350
351  if (verboseLevel > 1) 
352  {
353    G4cout << "Recoil: "
354           << nlv0<<" m= " << m2 << " ekin(MeV)= " << erec
355           <<G4endl;
356  }
357  if(erec > lowEnergyRecoilLimit) 
358  {
359    G4DynamicParticle * aSec = new G4DynamicParticle(theDef, nlv0);
360    theParticleChange.AddSecondary(aSec);
361  } else {
362    if(erec < 0.0) erec = 0.0;
363    theParticleChange.SetLocalEnergyDeposit(erec);
364  }
365
366  return &theParticleChange;
367}
368
369
370////////////////////////////////////////////////////////////////////////////
371//
372// return differential elastic cross section d(sigma)/d(omega)
373
374G4double
375G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc( const G4ParticleDefinition* particle, 
376                                        G4double theta, 
377                                        G4double momentum, 
378                                        G4double A         )
379{
380  fParticle      = particle;
381  fWaveVector    = momentum/hbarc;
382  fAtomicWeight  = A;
383  fAddCoulomb    = false;
384  fNuclearRadius = CalculateNuclearRad(A);
385
386  G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
387
388  return sigma;
389}
390
391////////////////////////////////////////////////////////////////////////////
392//
393// return invariant differential elastic cross section d(sigma)/d(tMand)
394
395G4double
396G4NuclNuclDiffuseElastic::GetInvElasticXsc( const G4ParticleDefinition* particle, 
397                                        G4double tMand, 
398                                        G4double plab, 
399                                        G4double A, G4double Z         )
400{
401  G4double m1 = particle->GetPDGMass();
402  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
403
404  G4int iZ = static_cast<G4int>(Z+0.5);
405  G4int iA = static_cast<G4int>(A+0.5);
406  G4ParticleDefinition * theDef = 0;
407
408  if      (iZ == 1 && iA == 1) theDef = theProton;
409  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
410  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
411  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
412  else if (iZ == 2 && iA == 4) theDef = theAlpha;
413  else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
414 
415  G4double tmass = theDef->GetPDGMass();
416
417  G4LorentzVector lv(0.0,0.0,0.0,tmass);   
418  lv += lv1;
419
420  G4ThreeVector bst = lv.boostVector();
421  lv1.boost(-bst);
422
423  G4ThreeVector p1 = lv1.vect();
424  G4double ptot    = p1.mag();
425  G4double ptot2 = ptot*ptot;
426  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
427
428  if( cost >= 1.0 )      cost = 1.0; 
429  else if( cost <= -1.0) cost = -1.0;
430 
431  G4double thetaCMS = std::acos(cost);
432
433  G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
434
435  sigma *= pi/ptot2;
436
437  return sigma;
438}
439
440////////////////////////////////////////////////////////////////////////////
441//
442// return differential elastic cross section d(sigma)/d(omega) with Coulomb
443// correction
444
445G4double
446G4NuclNuclDiffuseElastic::GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle, 
447                                        G4double theta, 
448                                        G4double momentum, 
449                                        G4double A, G4double Z         )
450{
451  fParticle      = particle;
452  fWaveVector    = momentum/hbarc;
453  fAtomicWeight  = A;
454  fAtomicNumber  = Z;
455  fNuclearRadius = CalculateNuclearRad(A);
456  fAddCoulomb    = false;
457
458  G4double z     = particle->GetPDGCharge();
459
460  G4double kRt   = fWaveVector*fNuclearRadius*theta;
461  G4double kRtC  = 1.9;
462
463  if( z && (kRt > kRtC) )
464  {
465    fAddCoulomb = true;
466    fBeta       = CalculateParticleBeta( particle, momentum);
467    fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
468    fAm         = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
469  }
470  G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
471
472  return sigma;
473}
474
475////////////////////////////////////////////////////////////////////////////
476//
477// return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
478// correction
479
480G4double
481G4NuclNuclDiffuseElastic::GetInvElasticSumXsc( const G4ParticleDefinition* particle, 
482                                        G4double tMand, 
483                                        G4double plab, 
484                                        G4double A, G4double Z         )
485{
486  G4double m1 = particle->GetPDGMass();
487
488  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
489
490  G4int iZ = static_cast<G4int>(Z+0.5);
491  G4int iA = static_cast<G4int>(A+0.5);
492
493  G4ParticleDefinition* theDef = 0;
494
495  if      (iZ == 1 && iA == 1) theDef = theProton;
496  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
497  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
498  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
499  else if (iZ == 2 && iA == 4) theDef = theAlpha;
500  else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
501 
502  G4double tmass = theDef->GetPDGMass();
503
504  G4LorentzVector lv(0.0,0.0,0.0,tmass);   
505  lv += lv1;
506
507  G4ThreeVector bst = lv.boostVector();
508  lv1.boost(-bst);
509
510  G4ThreeVector p1 = lv1.vect();
511  G4double ptot    = p1.mag();
512  G4double ptot2   = ptot*ptot;
513  G4double cost    = 1 - 0.5*std::fabs(tMand)/ptot2;
514
515  if( cost >= 1.0 )      cost = 1.0; 
516  else if( cost <= -1.0) cost = -1.0;
517 
518  G4double thetaCMS = std::acos(cost);
519
520  G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
521
522  sigma *= pi/ptot2;
523
524  return sigma;
525}
526
527////////////////////////////////////////////////////////////////////////////
528//
529// return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
530// correction
531
532G4double
533G4NuclNuclDiffuseElastic::GetInvCoulombElasticXsc( const G4ParticleDefinition* particle, 
534                                        G4double tMand, 
535                                        G4double plab, 
536                                        G4double A, G4double Z         )
537{
538  G4double m1 = particle->GetPDGMass();
539  G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
540
541  G4int iZ = static_cast<G4int>(Z+0.5);
542  G4int iA = static_cast<G4int>(A+0.5);
543  G4ParticleDefinition * theDef = 0;
544
545  if      (iZ == 1 && iA == 1) theDef = theProton;
546  else if (iZ == 1 && iA == 2) theDef = theDeuteron;
547  else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
548  else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
549  else if (iZ == 2 && iA == 4) theDef = theAlpha;
550  else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
551 
552  G4double tmass = theDef->GetPDGMass();
553
554  G4LorentzVector lv(0.0,0.0,0.0,tmass);   
555  lv += lv1;
556
557  G4ThreeVector bst = lv.boostVector();
558  lv1.boost(-bst);
559
560  G4ThreeVector p1 = lv1.vect();
561  G4double ptot    = p1.mag();
562  G4double ptot2 = ptot*ptot;
563  G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
564
565  if( cost >= 1.0 )      cost = 1.0; 
566  else if( cost <= -1.0) cost = -1.0;
567 
568  G4double thetaCMS = std::acos(cost);
569
570  G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
571
572  sigma *= pi/ptot2;
573
574  return sigma;
575}
576
577////////////////////////////////////////////////////////////////////////////
578//
579// return differential elastic probability d(probability)/d(omega)
580
581G4double
582G4NuclNuclDiffuseElastic::GetDiffElasticProb( // G4ParticleDefinition* particle,
583                                        G4double theta
584                                        //  G4double momentum,
585                                        // G4double A         
586                                     )
587{
588  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
589  G4double delta, diffuse, gamma;
590  G4double e1, e2, bone, bone2;
591
592  // G4double wavek = momentum/hbarc;  // wave vector
593  // G4double r0    = 1.08*fermi;
594  // G4double rad   = r0*std::pow(A, 1./3.);
595
596  G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
597  G4double kr2   = kr*kr;
598  G4double krt   = kr*theta;
599
600  bzero      = BesselJzero(krt);
601  bzero2     = bzero*bzero;   
602  bone       = BesselJone(krt);
603  bone2      = bone*bone;
604  bonebyarg  = BesselOneByArg(krt);
605  bonebyarg2 = bonebyarg*bonebyarg; 
606
607  if (fParticle == theProton)
608  {
609    diffuse = 0.63*fermi;
610    gamma   = 0.3*fermi;
611    delta   = 0.1*fermi*fermi;
612    e1      = 0.3*fermi;
613    e2      = 0.35*fermi;
614  }
615  else // as proton, if were not defined
616  {
617    diffuse = 0.63*fermi;
618    gamma   = 0.3*fermi;
619    delta   = 0.1*fermi*fermi;
620    e1      = 0.3*fermi;
621    e2      = 0.35*fermi;
622  }
623  G4double lambda = 15.; // 15 ok
624
625  //  G4double kg    = fWaveVector*gamma;   // wavek*delta;
626
627  G4double kg    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
628  G4double kg2   = kg*kg;
629
630  // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
631  // G4double dk2t2 = dk2t*dk2t;
632  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
633
634  G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
635
636  damp           = DampFactor(pikdt);
637  damp2          = damp*damp;
638
639  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector; 
640  G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
641
642
643  sigma  = kg2;
644  // sigma  += dk2t2;
645  sigma *= bzero2;
646  sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
647  sigma += kr2*bonebyarg2;
648  sigma *= damp2;          // *rad*rad;
649
650  return sigma;
651}
652
653////////////////////////////////////////////////////////////////////////////
654//
655// return differential elastic probability d(probability)/d(omega) with
656// Coulomb correction
657
658G4double
659G4NuclNuclDiffuseElastic::GetDiffElasticSumProb( // G4ParticleDefinition* particle,
660                                        G4double theta
661                                        //  G4double momentum,
662                                        // G4double A         
663                                     )
664{
665  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
666  G4double delta, diffuse, gamma;
667  G4double e1, e2, bone, bone2;
668
669  // G4double wavek = momentum/hbarc;  // wave vector
670  // G4double r0    = 1.08*fermi;
671  // G4double rad   = r0*std::pow(A, 1./3.);
672
673  G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
674  G4double kr2   = kr*kr;
675  G4double krt   = kr*theta;
676
677  bzero      = BesselJzero(krt);
678  bzero2     = bzero*bzero;   
679  bone       = BesselJone(krt);
680  bone2      = bone*bone;
681  bonebyarg  = BesselOneByArg(krt);
682  bonebyarg2 = bonebyarg*bonebyarg; 
683
684  if (fParticle == theProton)
685  {
686    diffuse = 0.63*fermi;
687    // diffuse = 0.6*fermi;
688    gamma   = 0.3*fermi;
689    delta   = 0.1*fermi*fermi;
690    e1      = 0.3*fermi;
691    e2      = 0.35*fermi;
692  }
693  else // as proton, if were not defined
694  {
695    diffuse = 0.63*fermi;
696    gamma   = 0.3*fermi;
697    delta   = 0.1*fermi*fermi;
698    e1      = 0.3*fermi;
699    e2      = 0.35*fermi;
700  }
701  G4double lambda = 15.; // 15 ok
702  // G4double kg    = fWaveVector*gamma;   // wavek*delta;
703  G4double kg    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
704
705  // G4cout<<"kg = "<<kg<<G4endl;
706
707  if(fAddCoulomb)  // add Coulomb correction
708  {
709    G4double sinHalfTheta  = std::sin(0.5*theta);
710    G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
711
712    kg += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
713  // kg += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
714  }
715
716  G4double kg2   = kg*kg;
717
718  // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
719  //   G4cout<<"dk2t = "<<dk2t<<G4endl;
720  // G4double dk2t2 = dk2t*dk2t;
721  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
722
723  G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
724
725  // G4cout<<"pikdt = "<<pikdt<<G4endl;
726
727  damp           = DampFactor(pikdt);
728  damp2          = damp*damp;
729
730  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector; 
731  G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
732
733  sigma  = kg2;
734  // sigma += dk2t2;
735  sigma *= bzero2;
736  sigma += mode2k2*bone2; 
737  sigma += e2dk3t*bzero*bone;
738
739  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2;  // correction at J1()/()
740  sigma += kr2*bonebyarg2;  // correction at J1()/()
741
742  sigma *= damp2;          // *rad*rad;
743
744  return sigma;
745}
746
747
748////////////////////////////////////////////////////////////////////////////
749//
750// return differential elastic probability d(probability)/d(t) with
751// Coulomb correction
752
753G4double
754G4NuclNuclDiffuseElastic::GetDiffElasticSumProbA( G4double alpha )
755{
756  G4double theta; 
757
758  theta = std::sqrt(alpha);
759
760  // theta = std::acos( 1 - alpha/2. );
761
762  G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
763  G4double delta, diffuse, gamma;
764  G4double e1, e2, bone, bone2;
765
766  // G4double wavek = momentum/hbarc;  // wave vector
767  // G4double r0    = 1.08*fermi;
768  // G4double rad   = r0*std::pow(A, 1./3.);
769
770  G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
771  G4double kr2   = kr*kr;
772  G4double krt   = kr*theta;
773
774  bzero      = BesselJzero(krt);
775  bzero2     = bzero*bzero;   
776  bone       = BesselJone(krt);
777  bone2      = bone*bone;
778  bonebyarg  = BesselOneByArg(krt);
779  bonebyarg2 = bonebyarg*bonebyarg; 
780
781  if (fParticle == theProton)
782  {
783    diffuse = 0.63*fermi;
784    // diffuse = 0.6*fermi;
785    gamma   = 0.3*fermi;
786    delta   = 0.1*fermi*fermi;
787    e1      = 0.3*fermi;
788    e2      = 0.35*fermi;
789  }
790  else // as proton, if were not defined
791  {
792    diffuse = 0.63*fermi;
793    gamma   = 0.3*fermi;
794    delta   = 0.1*fermi*fermi;
795    e1      = 0.3*fermi;
796    e2      = 0.35*fermi;
797  }
798  G4double lambda = 15.; // 15 ok
799  // G4double kg    = fWaveVector*gamma;   // wavek*delta;
800  G4double kg    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
801
802  // G4cout<<"kg = "<<kg<<G4endl;
803
804  if(fAddCoulomb)  // add Coulomb correction
805  {
806    G4double sinHalfTheta  = theta*0.5; // std::sin(0.5*theta);
807    G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
808
809    kg += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
810  // kg += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
811  }
812
813  G4double kg2   = kg*kg;
814
815  // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
816  //   G4cout<<"dk2t = "<<dk2t<<G4endl;
817  // G4double dk2t2 = dk2t*dk2t;
818  // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
819
820  G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
821
822  // G4cout<<"pikdt = "<<pikdt<<G4endl;
823
824  damp           = DampFactor(pikdt);
825  damp2          = damp*damp;
826
827  G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector; 
828  G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
829
830  sigma  = kg2;
831  // sigma += dk2t2;
832  sigma *= bzero2;
833  sigma += mode2k2*bone2; 
834  sigma += e2dk3t*bzero*bone;
835
836  // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2;  // correction at J1()/()
837  sigma += kr2*bonebyarg2;  // correction at J1()/()
838
839  sigma *= damp2;          // *rad*rad;
840
841  return sigma;
842}
843
844
845////////////////////////////////////////////////////////////////////////////
846//
847// return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega)
848
849G4double
850G4NuclNuclDiffuseElastic::GetIntegrandFunction( G4double alpha )
851{
852  G4double result;
853
854  result  = GetDiffElasticSumProbA(alpha);
855
856  // result *= 2*pi*std::sin(theta);
857
858  return result;
859}
860
861////////////////////////////////////////////////////////////////////////////
862//
863// return integral elastic cross section d(sigma)/d(omega) integrated 0 - theta
864
865G4double
866G4NuclNuclDiffuseElastic::IntegralElasticProb(  const G4ParticleDefinition* particle, 
867                                        G4double theta, 
868                                        G4double momentum, 
869                                        G4double A         )
870{
871  G4double result;
872  fParticle      = particle;
873  fWaveVector    = momentum/hbarc;
874  fAtomicWeight  = A;
875
876  fNuclearRadius = CalculateNuclearRad(A);
877
878
879  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
880
881  // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
882  result = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta ); 
883
884  return result;
885}
886
887////////////////////////////////////////////////////////////////////////////
888//
889// Return inv momentum transfer -t > 0
890
891G4double G4NuclNuclDiffuseElastic::SampleT( const G4ParticleDefinition* aParticle, G4double p, G4double A)
892{
893  G4double theta = SampleThetaCMS( aParticle,  p, A); // sample theta in cms
894  G4double t     = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
895  return t;
896}
897
898////////////////////////////////////////////////////////////////////////////
899//
900// Return scattering angle sampled in cms
901
902
903G4double
904G4NuclNuclDiffuseElastic::SampleThetaCMS(const G4ParticleDefinition* particle, 
905                                       G4double momentum, G4double A)
906{
907  G4int i, iMax = 100; 
908  G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
909
910  fParticle      = particle;
911  fWaveVector    = momentum/hbarc;
912  fAtomicWeight  = A;
913
914  fNuclearRadius = CalculateNuclearRad(A);
915 
916  thetaMax = 10.174/fWaveVector/fNuclearRadius;
917
918  if (thetaMax > pi) thetaMax = pi;
919
920  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
921
922  // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
923  norm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., thetaMax );
924
925  norm *= G4UniformRand();
926
927  for(i = 1; i <= iMax; i++)
928  {
929    theta1 = (i-1)*thetaMax/iMax; 
930    theta2 = i*thetaMax/iMax;
931    sum   += integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, theta1, theta2);
932
933    if ( sum >= norm ) 
934    {
935      result = 0.5*(theta1 + theta2);
936      break;
937    }
938  }
939  if (i > iMax ) result = 0.5*(theta1 + theta2);
940
941  G4double sigma = pi*thetaMax/iMax;
942
943  result += G4RandGauss::shoot(0.,sigma);
944
945  if(result < 0.) result = 0.;
946  if(result > thetaMax) result = thetaMax;
947
948  return result;
949}
950
951/////////////////////////////////////////////////////////////////////////////
952/////////////////////  Table preparation and reading ////////////////////////
953////////////////////////////////////////////////////////////////////////////
954//
955// Return inv momentum transfer -t > 0 from initialisation table
956
957G4double G4NuclNuclDiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p, 
958                                               G4double Z, G4double A)
959{
960  G4double alpha = SampleTableThetaCMS( aParticle,  p, Z, A); // sample theta2 in cms
961  // G4double t     = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) );             // -t !!!
962  G4double t     = p*p*alpha;             // -t !!!
963  return t;
964}
965
966////////////////////////////////////////////////////////////////////////////
967//
968// Return scattering angle2 sampled in cms according to precalculated table.
969
970
971G4double
972G4NuclNuclDiffuseElastic::SampleTableThetaCMS(const G4ParticleDefinition* particle, 
973                                       G4double momentum, G4double Z, G4double A)
974{
975  size_t iElement;
976  G4int iMomentum, iAngle; 
977  G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W; 
978  G4double m1 = particle->GetPDGMass();
979
980  for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
981  {
982    if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
983  }
984  if ( iElement == fElementNumberVector.size() ) 
985  {
986    InitialiseOnFly(Z,A); // table preparation, if needed
987
988    // iElement--;
989
990    // G4cout << "G4NuclNuclDiffuseElastic: Element with atomic number " << Z
991    // << " is not found, return zero angle" << G4endl;
992    // return 0.; // no table for this element
993  }
994  // G4cout<<"iElement = "<<iElement<<G4endl;
995
996  fAngleTable = fAngleBank[iElement];
997
998  G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
999
1000  for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
1001  {
1002    if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
1003  }
1004  if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;   // kinE is more then theMaxEnergy
1005  if ( iMomentum < 0 )           iMomentum = 0; // against negative index, kinE < theMinEnergy
1006
1007  // G4cout<<"iMomentum = "<<iMomentum<<G4endl;
1008
1009  if (iMomentum == fEnergyBin -1 || iMomentum == 0 )   // the table edges
1010  {
1011    position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
1012
1013    // G4cout<<"position = "<<position<<G4endl;
1014
1015    for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
1016    {
1017      if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
1018    }
1019    if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
1020
1021    // G4cout<<"iAngle = "<<iAngle<<G4endl;
1022
1023    randAngle = GetScatteringAngle(iMomentum, iAngle, position);
1024
1025    // G4cout<<"randAngle = "<<randAngle<<G4endl;
1026  }
1027  else  // kinE inside between energy table edges
1028  {
1029    // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
1030    position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
1031
1032    // G4cout<<"position = "<<position<<G4endl;
1033
1034    for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
1035    {
1036      // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
1037      if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
1038    }
1039    if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
1040
1041    // G4cout<<"iAngle = "<<iAngle<<G4endl;
1042
1043    theta2  = GetScatteringAngle(iMomentum, iAngle, position);
1044
1045    // G4cout<<"theta2 = "<<theta2<<G4endl;
1046    E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
1047
1048    // G4cout<<"E2 = "<<E2<<G4endl;
1049   
1050    iMomentum--;
1051   
1052    // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
1053
1054    // G4cout<<"position = "<<position<<G4endl;
1055
1056    for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
1057    {
1058      // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
1059      if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
1060    }
1061    if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
1062   
1063    theta1  = GetScatteringAngle(iMomentum, iAngle, position);
1064
1065    // G4cout<<"theta1 = "<<theta1<<G4endl;
1066
1067    E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
1068
1069    // G4cout<<"E1 = "<<E1<<G4endl;
1070
1071    W  = 1.0/(E2 - E1);
1072    W1 = (E2 - kinE)*W;
1073    W2 = (kinE - E1)*W;
1074
1075    randAngle = W1*theta1 + W2*theta2;
1076   
1077    // randAngle = theta2;
1078    // G4cout<<"randAngle = "<<randAngle<<G4endl;
1079  }
1080  // G4double angle = randAngle;
1081  // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
1082
1083  return randAngle;
1084}
1085
1086//////////////////////////////////////////////////////////////////////////////
1087//
1088// Initialisation for given particle on fly using new element number
1089
1090void G4NuclNuclDiffuseElastic::InitialiseOnFly(G4double Z, G4double A) 
1091{
1092  fAtomicNumber  = Z;     // atomic number
1093  fAtomicWeight  = A;     // number of nucleons
1094
1095  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1096 
1097  if( verboseLevel > 0 )   
1098  {
1099    G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element with Z = "
1100          <<Z<<"; and A = "<<A<<G4endl;
1101  }
1102  fElementNumberVector.push_back(fAtomicNumber);
1103
1104  BuildAngleTable();
1105
1106  fAngleBank.push_back(fAngleTable);
1107
1108  return;
1109}
1110
1111///////////////////////////////////////////////////////////////////////////////
1112//
1113// Build for given particle and element table of momentum, angle probability.
1114// For the moment in lab system.
1115
1116void G4NuclNuclDiffuseElastic::BuildAngleTable() 
1117{
1118  G4int i, j;
1119  G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1120  G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1121
1122  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1123 
1124  fAngleTable = new G4PhysicsTable(fEnergyBin);
1125
1126  for( i = 0; i < fEnergyBin; i++)
1127  {
1128    kinE        = fEnergyVector->GetLowEdgeEnergy(i);
1129    partMom     = std::sqrt( kinE*(kinE + 2*m1) );
1130
1131    fWaveVector = partMom/hbarc;
1132
1133    G4double kR     = fWaveVector*fNuclearRadius;
1134    G4double kR2    = kR*kR;
1135    G4double kRmax  = 18.6; // 10.6; 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1136    G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; // on the first slope of J1
1137    // G4double kRlim  = 1.2;
1138    // G4double kRlim2 = kRlim*kRlim/kR2;
1139
1140    alphaMax = kRmax*kRmax/kR2;
1141
1142    if (alphaMax > 4.) alphaMax = 4.;  // vmg05-02-09: was pi2
1143
1144    alphaCoulomb = kRcoul*kRcoul/kR2;
1145
1146    if( z )
1147    {
1148      a           = partMom/m1;         // beta*gamma for m1
1149      fBeta       = a/std::sqrt(1+a*a);
1150      fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1151      fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1152    }
1153    G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1154
1155    // G4PhysicsLogVector*  angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1156
1157    G4double delth = alphaMax/fAngleBin;
1158       
1159    sum = 0.;
1160
1161    // fAddCoulomb = false;
1162    fAddCoulomb = true;
1163
1164    // for(j = 1; j < fAngleBin; j++)
1165    for(j = fAngleBin-1; j >= 1; j--)
1166    {
1167      // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1168      // alpha2 = angleBins->GetLowEdgeEnergy(j);
1169
1170      // alpha1 = alphaMax*(j-1)/fAngleBin;
1171      // alpha2 = alphaMax*( j )/fAngleBin;
1172
1173      alpha1 = delth*(j-1);
1174      // if(alpha1 < kRlim2) alpha1 = kRlim2;
1175      alpha2 = alpha1 + delth;
1176
1177      // if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1178      if( ( alpha1 < alphaCoulomb ) && z ) fAddCoulomb = false;
1179
1180      delta = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1181      // delta = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1182
1183      sum += delta;
1184     
1185      angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
1186      // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2<<"; sum = "<<sum<<G4endl;
1187    }
1188    fAngleTable->insertAt(i,angleVector);
1189
1190    // delete[] angleVector;
1191    // delete[] angleBins;
1192  }
1193  return;
1194}
1195
1196/////////////////////////////////////////////////////////////////////////////////
1197//
1198//
1199
1200G4double
1201G4NuclNuclDiffuseElastic:: GetScatteringAngle( G4int iMomentum, G4int iAngle, G4double position )
1202{
1203 G4double x1, x2, y1, y2, randAngle;
1204
1205  if( iAngle == 0 )
1206  {
1207    randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1208    // iAngle++;
1209  }
1210  else
1211  {
1212    if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1213    {
1214      iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1215    }
1216    y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1217    y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1218
1219    x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1220    x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1221
1222    if ( x1 == x2 )   randAngle = x2;
1223    else
1224    {
1225      if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
1226      else
1227      {
1228        randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1229      }
1230    }
1231  }
1232  return randAngle;
1233}
1234
1235
1236
1237////////////////////////////////////////////////////////////////////////////
1238//
1239// Return scattering angle sampled in lab system (target at rest)
1240
1241
1242
1243G4double
1244G4NuclNuclDiffuseElastic::SampleThetaLab( const G4HadProjectile* aParticle, 
1245                                        G4double tmass, G4double A)
1246{
1247  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1248  G4double m1 = theParticle->GetPDGMass();
1249  G4double plab = aParticle->GetTotalMomentum();
1250  G4LorentzVector lv1 = aParticle->Get4Momentum();
1251  G4LorentzVector lv(0.0,0.0,0.0,tmass);   
1252  lv += lv1;
1253
1254  G4ThreeVector bst = lv.boostVector();
1255  lv1.boost(-bst);
1256
1257  G4ThreeVector p1 = lv1.vect();
1258  G4double ptot    = p1.mag();
1259  G4double tmax    = 4.0*ptot*ptot;
1260  G4double t       = 0.0;
1261
1262
1263  //
1264  // Sample t
1265  //
1266 
1267  t = SampleT( theParticle, ptot, A);
1268
1269  // NaN finder
1270  if(!(t < 0.0 || t >= 0.0)) 
1271  {
1272    if (verboseLevel > 0) 
1273    {
1274      G4cout << "G4NuclNuclDiffuseElastic:WARNING: A = " << A
1275             << " mom(GeV)= " << plab/GeV
1276             << " S-wave will be sampled" 
1277             << G4endl; 
1278    }
1279    t = G4UniformRand()*tmax; 
1280  }
1281  if(verboseLevel>1)
1282  {
1283    G4cout <<" t= " << t << " tmax= " << tmax
1284           << " ptot= " << ptot << G4endl;
1285  }
1286  // Sampling of angles in CM system
1287
1288  G4double phi  = G4UniformRand()*twopi;
1289  G4double cost = 1. - 2.0*t/tmax;
1290  G4double sint;
1291
1292  if( cost >= 1.0 ) 
1293  {
1294    cost = 1.0;
1295    sint = 0.0;
1296  }
1297  else if( cost <= -1.0) 
1298  {
1299    cost = -1.0;
1300    sint =  0.0;
1301  }
1302  else 
1303  {
1304    sint = std::sqrt((1.0-cost)*(1.0+cost));
1305  }   
1306  if (verboseLevel>1) 
1307  {
1308    G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
1309  }
1310  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1311  v1 *= ptot;
1312  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1313
1314  nlv1.boost(bst); 
1315
1316  G4ThreeVector np1 = nlv1.vect();
1317
1318    // G4double theta = std::acos( np1.z()/np1.mag() );  // degree;
1319
1320  G4double theta = np1.theta();
1321
1322  return theta;
1323}
1324
1325////////////////////////////////////////////////////////////////////////////
1326//
1327// Return scattering angle in lab system (target at rest) knowing theta in CMS
1328
1329
1330
1331G4double
1332G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab( const G4DynamicParticle* aParticle, 
1333                                        G4double tmass, G4double thetaCMS)
1334{
1335  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1336  G4double m1 = theParticle->GetPDGMass();
1337  // G4double plab = aParticle->GetTotalMomentum();
1338  G4LorentzVector lv1 = aParticle->Get4Momentum();
1339  G4LorentzVector lv(0.0,0.0,0.0,tmass);   
1340
1341  lv += lv1;
1342
1343  G4ThreeVector bst = lv.boostVector();
1344
1345  lv1.boost(-bst);
1346
1347  G4ThreeVector p1 = lv1.vect();
1348  G4double ptot    = p1.mag();
1349
1350  G4double phi  = G4UniformRand()*twopi;
1351  G4double cost = std::cos(thetaCMS);
1352  G4double sint;
1353
1354  if( cost >= 1.0 ) 
1355  {
1356    cost = 1.0;
1357    sint = 0.0;
1358  }
1359  else if( cost <= -1.0) 
1360  {
1361    cost = -1.0;
1362    sint =  0.0;
1363  }
1364  else 
1365  {
1366    sint = std::sqrt((1.0-cost)*(1.0+cost));
1367  }   
1368  if (verboseLevel>1) 
1369  {
1370    G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
1371  }
1372  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1373  v1 *= ptot;
1374  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1375
1376  nlv1.boost(bst); 
1377
1378  G4ThreeVector np1 = nlv1.vect();
1379
1380
1381  G4double thetaLab = np1.theta();
1382
1383  return thetaLab;
1384}
1385
1386////////////////////////////////////////////////////////////////////////////
1387//
1388// Return scattering angle in CMS system (target at rest) knowing theta in Lab
1389
1390
1391
1392G4double
1393G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS( const G4DynamicParticle* aParticle, 
1394                                        G4double tmass, G4double thetaLab)
1395{
1396  const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1397  G4double m1 = theParticle->GetPDGMass();
1398  G4double plab = aParticle->GetTotalMomentum();
1399  G4LorentzVector lv1 = aParticle->Get4Momentum();
1400  G4LorentzVector lv(0.0,0.0,0.0,tmass);   
1401
1402  lv += lv1;
1403
1404  G4ThreeVector bst = lv.boostVector();
1405
1406  // lv1.boost(-bst);
1407
1408  // G4ThreeVector p1 = lv1.vect();
1409  // G4double ptot    = p1.mag();
1410
1411  G4double phi  = G4UniformRand()*twopi;
1412  G4double cost = std::cos(thetaLab);
1413  G4double sint;
1414
1415  if( cost >= 1.0 ) 
1416  {
1417    cost = 1.0;
1418    sint = 0.0;
1419  }
1420  else if( cost <= -1.0) 
1421  {
1422    cost = -1.0;
1423    sint =  0.0;
1424  }
1425  else 
1426  {
1427    sint = std::sqrt((1.0-cost)*(1.0+cost));
1428  }   
1429  if (verboseLevel>1) 
1430  {
1431    G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
1432  }
1433  G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1434  v1 *= plab;
1435  G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1436
1437  nlv1.boost(-bst); 
1438
1439  G4ThreeVector np1 = nlv1.vect();
1440
1441
1442  G4double thetaCMS = np1.theta();
1443
1444  return thetaCMS;
1445}
1446
1447///////////////////////////////////////////////////////////////////////////////
1448//
1449// Test for given particle and element table of momentum, angle probability.
1450// For the moment in lab system.
1451
1452void G4NuclNuclDiffuseElastic::TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
1453                                            G4double Z, G4double A) 
1454{
1455  fAtomicNumber  = Z;     // atomic number
1456  fAtomicWeight  = A;     // number of nucleons
1457  fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1458 
1459     
1460 
1461  G4cout<<"G4NuclNuclDiffuseElastic::TestAngleTable() init the element with Z = "
1462          <<Z<<"; and A = "<<A<<G4endl;
1463 
1464  fElementNumberVector.push_back(fAtomicNumber);
1465
1466 
1467
1468
1469  G4int i=0, j;
1470  G4double a = 0., z = theParticle->GetPDGCharge(),  m1 = fParticle->GetPDGMass();
1471  G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1472  G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1473  G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1474  G4double epsilon = 0.001;
1475
1476  G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1477 
1478  fAngleTable = new G4PhysicsTable(fEnergyBin);
1479
1480  fWaveVector = partMom/hbarc;
1481
1482  G4double kR     = fWaveVector*fNuclearRadius;
1483  G4double kR2    = kR*kR;
1484  G4double kRmax  = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1485  G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
1486
1487  alphaMax = kRmax*kRmax/kR2;
1488
1489  if (alphaMax > 4.) alphaMax = 4.;  // vmg05-02-09: was pi2
1490
1491  alphaCoulomb = kRcoul*kRcoul/kR2;
1492
1493  if( z )
1494  {
1495      a           = partMom/m1; // beta*gamma for m1
1496      fBeta       = a/std::sqrt(1+a*a);
1497      fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1498      fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1499  }
1500  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1501
1502  // G4PhysicsLogVector*  angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1503   
1504 
1505  fAddCoulomb = false;
1506
1507  for(j = 1; j < fAngleBin; j++)
1508  {
1509      // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1510      // alpha2 = angleBins->GetLowEdgeEnergy(j);
1511
1512    alpha1 = alphaMax*(j-1)/fAngleBin;
1513    alpha2 = alphaMax*( j )/fAngleBin;
1514
1515    if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1516
1517    deltaL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1518    deltaL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1519    deltaAG  = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 
1520                                       alpha1, alpha2,epsilon);
1521
1522      // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1523      //     <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1524
1525    sumL10 += deltaL10;
1526    sumL96 += deltaL96;
1527    sumAG  += deltaAG;
1528
1529    G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1530            <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1531     
1532    angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
1533  }
1534  fAngleTable->insertAt(i,angleVector);
1535  fAngleBank.push_back(fAngleTable);
1536
1537  /*
1538  // Integral over all angle range - Bad accuracy !!!
1539
1540  sumL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1541  sumL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1542  sumAG  = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction,
1543                                       0., alpha2,epsilon);
1544  G4cout<<G4endl;
1545  G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1546            <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1547  */
1548  return;
1549}
1550
1551//
1552//
1553/////////////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.