Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/hadronic/models/coherent_elastic/src/G4DiffuseElastic.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/hadronic/models/coherent_elastic/src/G4DiffuseElastic.cc (Version 11.3.0) and /processes/hadronic/models/coherent_elastic/src/G4DiffuseElastic.cc (Version 9.5.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: G4DiffuseElastic.cc,v 1.25 2009-09-22 16:21:46 vnivanch Exp $
                                                   >>  27 // GEANT4 tag $Name: not supported by cvs2svn $
 26 //                                                 28 //
 27 //                                                 29 //
 28 // Physics model class G4DiffuseElastic            30 // Physics model class G4DiffuseElastic 
 29 //                                                 31 //
 30 //                                                 32 //
 31 // G4 Model: optical diffuse elastic scatterin     33 // G4 Model: optical diffuse elastic scattering with 4-momentum balance
 32 //                                                 34 //                         
 33 // 24-May-07 V. Grichine                           35 // 24-May-07 V. Grichine
 34 //                                                 36 //
 35 // 21.10.15 V. Grichine                        << 
 36 //             Bug fixed in BuildAngleTable, i << 
 37 //             angle bins at high energies > 5 << 
 38 //                                             << 
 39                                                    37 
 40 #include "G4DiffuseElastic.hh"                     38 #include "G4DiffuseElastic.hh"
 41 #include "G4ParticleTable.hh"                      39 #include "G4ParticleTable.hh"
 42 #include "G4ParticleDefinition.hh"                 40 #include "G4ParticleDefinition.hh"
 43 #include "G4IonTable.hh"                           41 #include "G4IonTable.hh"
 44 #include "G4NucleiProperties.hh"                   42 #include "G4NucleiProperties.hh"
 45                                                    43 
 46 #include "Randomize.hh"                            44 #include "Randomize.hh"
 47 #include "G4Integrator.hh"                         45 #include "G4Integrator.hh"
 48 #include "globals.hh"                              46 #include "globals.hh"
 49 #include "G4PhysicalConstants.hh"              << 
 50 #include "G4SystemOfUnits.hh"                  << 
 51                                                    47 
 52 #include "G4Proton.hh"                             48 #include "G4Proton.hh"
 53 #include "G4Neutron.hh"                            49 #include "G4Neutron.hh"
 54 #include "G4Deuteron.hh"                           50 #include "G4Deuteron.hh"
 55 #include "G4Alpha.hh"                              51 #include "G4Alpha.hh"
 56 #include "G4PionPlus.hh"                           52 #include "G4PionPlus.hh"
 57 #include "G4PionMinus.hh"                          53 #include "G4PionMinus.hh"
 58                                                    54 
 59 #include "G4Element.hh"                            55 #include "G4Element.hh"
 60 #include "G4ElementTable.hh"                       56 #include "G4ElementTable.hh"
 61 #include "G4NistManager.hh"                    << 
 62 #include "G4PhysicsTable.hh"                       57 #include "G4PhysicsTable.hh"
 63 #include "G4PhysicsLogVector.hh"                   58 #include "G4PhysicsLogVector.hh"
 64 #include "G4PhysicsFreeVector.hh"                  59 #include "G4PhysicsFreeVector.hh"
 65                                                    60 
 66 #include "G4Exp.hh"                            << 
 67                                                << 
 68 #include "G4HadronicParameters.hh"             << 
 69                                                << 
 70 //////////////////////////////////////////////     61 /////////////////////////////////////////////////////////////////////////
 71 //                                                 62 //
 72 // Test Constructor. Just to check xsc             63 // Test Constructor. Just to check xsc
 73                                                    64 
 74                                                    65 
 75 G4DiffuseElastic::G4DiffuseElastic()               66 G4DiffuseElastic::G4DiffuseElastic() 
 76   : G4HadronElastic("DiffuseElastic"), fPartic     67   : G4HadronElastic("DiffuseElastic"), fParticle(0)
 77 {                                                  68 {
 78   SetMinEnergy( 0.01*MeV );                    <<  69   SetMinEnergy( 0.01*GeV );
 79   SetMaxEnergy( G4HadronicParameters::Instance <<  70   SetMaxEnergy( 1.*TeV );
 80                                                <<  71   verboseLevel = 0;
 81   verboseLevel         = 0;                    << 
 82   lowEnergyRecoilLimit = 100.*keV;                 72   lowEnergyRecoilLimit = 100.*keV;  
 83   lowEnergyLimitQ      = 0.0*GeV;              <<  73   lowEnergyLimitQ  = 0.0*GeV;  
 84   lowEnergyLimitHE     = 0.0*GeV;              <<  74   lowEnergyLimitHE = 0.0*GeV;  
 85   lowestEnergyLimit    = 0.0*keV;              <<  75   lowestEnergyLimit= 0.0*keV;  
 86   plabLowLimit         = 20.0*MeV;             <<  76   plabLowLimit     = 20.0*MeV;
 87                                                <<  77 
 88   theProton    = G4Proton::Proton();           <<  78   theProton   = G4Proton::Proton();
 89   theNeutron   = G4Neutron::Neutron();         <<  79   theNeutron  = G4Neutron::Neutron();
 90   theDeuteron  = G4Deuteron::Deuteron();       <<  80   theDeuteron = G4Deuteron::Deuteron();
 91   theAlpha     = G4Alpha::Alpha();             <<  81   theAlpha    = G4Alpha::Alpha();
 92   thePionPlus  = G4PionPlus::PionPlus();       <<  82   thePionPlus = G4PionPlus::PionPlus();
 93   thePionMinus = G4PionMinus::PionMinus();     <<  83   thePionMinus= G4PionMinus::PionMinus();
 94                                                    84 
 95   fEnergyBin = 300;  // Increased from the ori <<  85   fEnergyBin = 200;
 96   fAngleBin  = 200;                                86   fAngleBin  = 200;
 97                                                    87 
 98   fEnergyVector =  new G4PhysicsLogVector( the     88   fEnergyVector =  new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
 99                                                << 
100   fAngleTable = 0;                                 89   fAngleTable = 0;
101                                                    90 
102   fParticle      = 0;                          <<  91   fParticle = 0;
103   fWaveVector    = 0.;                         <<  92   fWaveVector = 0.;
104   fAtomicWeight  = 0.;                         <<  93   fAtomicWeight = 0.;
105   fAtomicNumber  = 0.;                         <<  94   fAtomicNumber = 0.;
106   fNuclearRadius = 0.;                             95   fNuclearRadius = 0.;
107   fBeta          = 0.;                         <<  96   fBeta = 0.;
108   fZommerfeld    = 0.;                         <<  97   fZommerfeld = 0.;
109   fAm = 0.;                                        98   fAm = 0.;
110   fAddCoulomb = false;                             99   fAddCoulomb = false;
111 }                                                 100 }
112                                                   101 
113 //////////////////////////////////////////////    102 //////////////////////////////////////////////////////////////////////////////
114 //                                                103 //
115 // Destructor                                     104 // Destructor
116                                                   105 
117 G4DiffuseElastic::~G4DiffuseElastic()             106 G4DiffuseElastic::~G4DiffuseElastic()
118 {                                                 107 {
119   if ( fEnergyVector )                         << 108   if(fEnergyVector) delete fEnergyVector;
120   {                                            << 
121     delete fEnergyVector;                      << 
122     fEnergyVector = 0;                         << 
123   }                                            << 
124   for ( std::vector<G4PhysicsTable*>::iterator << 
125         it != fAngleBank.end(); ++it )         << 
126   {                                            << 
127     if ( (*it) ) (*it)->clearAndDestroy();     << 
128                                                   109 
129     delete *it;                                << 110   if( fAngleTable )
130     *it = 0;                                   << 111   {
                                                   >> 112       fAngleTable->clearAndDestroy();
                                                   >> 113       delete fAngleTable ;
131   }                                               114   }
132   fAngleTable = 0;                             << 
133 }                                                 115 }
134                                                   116 
135 //////////////////////////////////////////////    117 //////////////////////////////////////////////////////////////////////////////
136 //                                                118 //
137 // Initialisation for given particle using ele    119 // Initialisation for given particle using element table of application
138                                                   120 
139 void G4DiffuseElastic::Initialise()               121 void G4DiffuseElastic::Initialise() 
140 {                                                 122 {
141                                                   123 
142   // fEnergyVector = new G4PhysicsLogVector( t    124   // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
143                                                   125 
144   const G4ElementTable* theElementTable = G4El    126   const G4ElementTable* theElementTable = G4Element::GetElementTable();
145                                                   127 
146   std::size_t jEl, numOfEl = G4Element::GetNum << 128   size_t jEl, numOfEl = G4Element::GetNumberOfElements();
147                                                   129 
148   for( jEl = 0; jEl < numOfEl; ++jEl) // appli << 130   for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
149   {                                               131   {
150     fAtomicNumber = (*theElementTable)[jEl]->G    132     fAtomicNumber = (*theElementTable)[jEl]->GetZ();     // atomic number
151     fAtomicWeight = G4NistManager::Instance()- << 133     fAtomicWeight = (*theElementTable)[jEl]->GetN();     // number of nucleons
152     fNuclearRadius = CalculateNuclearRad(fAtom    134     fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
153                                                   135 
154     if( verboseLevel > 0 )                     << 136     if(verboseLevel > 0) 
155     {                                             137     {   
156       G4cout<<"G4DiffuseElastic::Initialise()     138       G4cout<<"G4DiffuseElastic::Initialise() the element: "
157       <<(*theElementTable)[jEl]->GetName()<<G4    139       <<(*theElementTable)[jEl]->GetName()<<G4endl;
158     }                                             140     }
159     fElementNumberVector.push_back(fAtomicNumb    141     fElementNumberVector.push_back(fAtomicNumber);
160     fElementNameVector.push_back((*theElementT    142     fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
161                                                   143 
162     BuildAngleTable();                            144     BuildAngleTable();
163     fAngleBank.push_back(fAngleTable);            145     fAngleBank.push_back(fAngleTable);
164   }                                               146   }  
165   return;                                         147   return;
166 }                                                 148 }
167                                                   149 
168 //////////////////////////////////////////////    150 ////////////////////////////////////////////////////////////////////////////
169 //                                                151 //
170 // return differential elastic cross section d    152 // return differential elastic cross section d(sigma)/d(omega) 
171                                                   153 
172 G4double                                          154 G4double 
173 G4DiffuseElastic::GetDiffuseElasticXsc( const     155 G4DiffuseElastic::GetDiffuseElasticXsc( const G4ParticleDefinition* particle, 
174                                         G4doub    156                                         G4double theta, 
175                       G4double momentum,          157                       G4double momentum, 
176                                         G4doub    158                                         G4double A         )
177 {                                                 159 {
178   fParticle      = particle;                      160   fParticle      = particle;
179   fWaveVector    = momentum/hbarc;                161   fWaveVector    = momentum/hbarc;
180   fAtomicWeight  = A;                             162   fAtomicWeight  = A;
181   fAddCoulomb    = false;                         163   fAddCoulomb    = false;
182   fNuclearRadius = CalculateNuclearRad(A);        164   fNuclearRadius = CalculateNuclearRad(A);
183                                                   165 
184   G4double sigma = fNuclearRadius*fNuclearRadi    166   G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
185                                                   167 
186   return sigma;                                   168   return sigma;
187 }                                                 169 }
188                                                   170 
189 //////////////////////////////////////////////    171 ////////////////////////////////////////////////////////////////////////////
190 //                                                172 //
191 // return invariant differential elastic cross    173 // return invariant differential elastic cross section d(sigma)/d(tMand) 
192                                                   174 
193 G4double                                          175 G4double 
194 G4DiffuseElastic::GetInvElasticXsc( const G4Pa    176 G4DiffuseElastic::GetInvElasticXsc( const G4ParticleDefinition* particle, 
195                                         G4doub    177                                         G4double tMand, 
196                       G4double plab,              178                       G4double plab, 
197                                         G4doub    179                                         G4double A, G4double Z         )
198 {                                                 180 {
199   G4double m1 = particle->GetPDGMass();           181   G4double m1 = particle->GetPDGMass();
200   G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla    182   G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
201                                                   183 
202   G4int iZ = static_cast<G4int>(Z+0.5);           184   G4int iZ = static_cast<G4int>(Z+0.5);
203   G4int iA = static_cast<G4int>(A+0.5);           185   G4int iA = static_cast<G4int>(A+0.5);
204   G4ParticleDefinition * theDef = 0;              186   G4ParticleDefinition * theDef = 0;
205                                                   187 
206   if      (iZ == 1 && iA == 1) theDef = thePro    188   if      (iZ == 1 && iA == 1) theDef = theProton;
207   else if (iZ == 1 && iA == 2) theDef = theDeu    189   else if (iZ == 1 && iA == 2) theDef = theDeuteron;
208   else if (iZ == 1 && iA == 3) theDef = G4Trit    190   else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
209   else if (iZ == 2 && iA == 3) theDef = G4He3:    191   else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
210   else if (iZ == 2 && iA == 4) theDef = theAlp    192   else if (iZ == 2 && iA == 4) theDef = theAlpha;
211   else theDef = G4ParticleTable::GetParticleTa << 193   else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
212                                                   194  
213   G4double tmass = theDef->GetPDGMass();          195   G4double tmass = theDef->GetPDGMass();
214                                                   196 
215   G4LorentzVector lv(0.0,0.0,0.0,tmass);          197   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
216   lv += lv1;                                      198   lv += lv1;
217                                                   199 
218   G4ThreeVector bst = lv.boostVector();           200   G4ThreeVector bst = lv.boostVector();
219   lv1.boost(-bst);                                201   lv1.boost(-bst);
220                                                   202 
221   G4ThreeVector p1 = lv1.vect();                  203   G4ThreeVector p1 = lv1.vect();
222   G4double ptot    = p1.mag();                    204   G4double ptot    = p1.mag();
223   G4double ptot2 = ptot*ptot;                     205   G4double ptot2 = ptot*ptot;
224   G4double cost = 1 - 0.5*std::fabs(tMand)/pto    206   G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
225                                                   207 
226   if( cost >= 1.0 )      cost = 1.0;              208   if( cost >= 1.0 )      cost = 1.0;  
227   else if( cost <= -1.0) cost = -1.0;             209   else if( cost <= -1.0) cost = -1.0;
228                                                   210   
229   G4double thetaCMS = std::acos(cost);            211   G4double thetaCMS = std::acos(cost);
230                                                   212 
231   G4double sigma = GetDiffuseElasticXsc( parti    213   G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
232                                                   214 
233   sigma *= pi/ptot2;                              215   sigma *= pi/ptot2;
234                                                   216 
235   return sigma;                                   217   return sigma;
236 }                                                 218 }
237                                                   219 
238 //////////////////////////////////////////////    220 ////////////////////////////////////////////////////////////////////////////
239 //                                                221 //
240 // return differential elastic cross section d    222 // return differential elastic cross section d(sigma)/d(omega) with Coulomb
241 // correction                                     223 // correction
242                                                   224 
243 G4double                                          225 G4double 
244 G4DiffuseElastic::GetDiffuseElasticSumXsc( con    226 G4DiffuseElastic::GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle, 
245                                         G4doub    227                                         G4double theta, 
246                       G4double momentum,          228                       G4double momentum, 
247                                         G4doub    229                                         G4double A, G4double Z         )
248 {                                                 230 {
249   fParticle      = particle;                      231   fParticle      = particle;
250   fWaveVector    = momentum/hbarc;                232   fWaveVector    = momentum/hbarc;
251   fAtomicWeight  = A;                             233   fAtomicWeight  = A;
252   fAtomicNumber  = Z;                             234   fAtomicNumber  = Z;
253   fNuclearRadius = CalculateNuclearRad(A);        235   fNuclearRadius = CalculateNuclearRad(A);
254   fAddCoulomb    = false;                         236   fAddCoulomb    = false;
255                                                   237 
256   G4double z     = particle->GetPDGCharge();      238   G4double z     = particle->GetPDGCharge();
257                                                   239 
258   G4double kRt   = fWaveVector*fNuclearRadius*    240   G4double kRt   = fWaveVector*fNuclearRadius*theta;
259   G4double kRtC  = 1.9;                           241   G4double kRtC  = 1.9;
260                                                   242 
261   if( z && (kRt > kRtC) )                         243   if( z && (kRt > kRtC) )
262   {                                               244   {
263     fAddCoulomb = true;                           245     fAddCoulomb = true;
264     fBeta       = CalculateParticleBeta( parti    246     fBeta       = CalculateParticleBeta( particle, momentum);
265     fZommerfeld = CalculateZommerfeld( fBeta,     247     fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
266     fAm         = CalculateAm( momentum, fZomm    248     fAm         = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
267   }                                               249   }
268   G4double sigma = fNuclearRadius*fNuclearRadi    250   G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
269                                                   251 
270   return sigma;                                   252   return sigma;
271 }                                                 253 }
272                                                   254 
273 //////////////////////////////////////////////    255 ////////////////////////////////////////////////////////////////////////////
274 //                                                256 //
275 // return invariant differential elastic cross    257 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
276 // correction                                     258 // correction
277                                                   259 
278 G4double                                          260 G4double 
279 G4DiffuseElastic::GetInvElasticSumXsc( const G    261 G4DiffuseElastic::GetInvElasticSumXsc( const G4ParticleDefinition* particle, 
280                                         G4doub    262                                         G4double tMand, 
281                       G4double plab,              263                       G4double plab, 
282                                         G4doub    264                                         G4double A, G4double Z         )
283 {                                                 265 {
284   G4double m1 = particle->GetPDGMass();           266   G4double m1 = particle->GetPDGMass();
285                                                   267 
286   G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla    268   G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
287                                                   269 
288   G4int iZ = static_cast<G4int>(Z+0.5);           270   G4int iZ = static_cast<G4int>(Z+0.5);
289   G4int iA = static_cast<G4int>(A+0.5);           271   G4int iA = static_cast<G4int>(A+0.5);
290                                                   272 
291   G4ParticleDefinition* theDef = 0;               273   G4ParticleDefinition* theDef = 0;
292                                                   274 
293   if      (iZ == 1 && iA == 1) theDef = thePro    275   if      (iZ == 1 && iA == 1) theDef = theProton;
294   else if (iZ == 1 && iA == 2) theDef = theDeu    276   else if (iZ == 1 && iA == 2) theDef = theDeuteron;
295   else if (iZ == 1 && iA == 3) theDef = G4Trit    277   else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
296   else if (iZ == 2 && iA == 3) theDef = G4He3:    278   else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
297   else if (iZ == 2 && iA == 4) theDef = theAlp    279   else if (iZ == 2 && iA == 4) theDef = theAlpha;
298   else theDef = G4ParticleTable::GetParticleTa << 280   else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
299                                                   281  
300   G4double tmass = theDef->GetPDGMass();          282   G4double tmass = theDef->GetPDGMass();
301                                                   283 
302   G4LorentzVector lv(0.0,0.0,0.0,tmass);          284   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
303   lv += lv1;                                      285   lv += lv1;
304                                                   286 
305   G4ThreeVector bst = lv.boostVector();           287   G4ThreeVector bst = lv.boostVector();
306   lv1.boost(-bst);                                288   lv1.boost(-bst);
307                                                   289 
308   G4ThreeVector p1 = lv1.vect();                  290   G4ThreeVector p1 = lv1.vect();
309   G4double ptot    = p1.mag();                    291   G4double ptot    = p1.mag();
310   G4double ptot2   = ptot*ptot;                   292   G4double ptot2   = ptot*ptot;
311   G4double cost    = 1 - 0.5*std::fabs(tMand)/    293   G4double cost    = 1 - 0.5*std::fabs(tMand)/ptot2;
312                                                   294 
313   if( cost >= 1.0 )      cost = 1.0;              295   if( cost >= 1.0 )      cost = 1.0;  
314   else if( cost <= -1.0) cost = -1.0;             296   else if( cost <= -1.0) cost = -1.0;
315                                                   297   
316   G4double thetaCMS = std::acos(cost);            298   G4double thetaCMS = std::acos(cost);
317                                                   299 
318   G4double sigma = GetDiffuseElasticSumXsc( pa    300   G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
319                                                   301 
320   sigma *= pi/ptot2;                              302   sigma *= pi/ptot2;
321                                                   303 
322   return sigma;                                   304   return sigma;
323 }                                                 305 }
324                                                   306 
325 //////////////////////////////////////////////    307 ////////////////////////////////////////////////////////////////////////////
326 //                                                308 //
327 // return invariant differential elastic cross    309 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
328 // correction                                     310 // correction
329                                                   311 
330 G4double                                          312 G4double 
331 G4DiffuseElastic::GetInvCoulombElasticXsc( con    313 G4DiffuseElastic::GetInvCoulombElasticXsc( const G4ParticleDefinition* particle, 
332                                         G4doub    314                                         G4double tMand, 
333                       G4double plab,              315                       G4double plab, 
334                                         G4doub    316                                         G4double A, G4double Z         )
335 {                                                 317 {
336   G4double m1 = particle->GetPDGMass();           318   G4double m1 = particle->GetPDGMass();
337   G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla    319   G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
338                                                   320 
339   G4int iZ = static_cast<G4int>(Z+0.5);           321   G4int iZ = static_cast<G4int>(Z+0.5);
340   G4int iA = static_cast<G4int>(A+0.5);           322   G4int iA = static_cast<G4int>(A+0.5);
341   G4ParticleDefinition * theDef = 0;              323   G4ParticleDefinition * theDef = 0;
342                                                   324 
343   if      (iZ == 1 && iA == 1) theDef = thePro    325   if      (iZ == 1 && iA == 1) theDef = theProton;
344   else if (iZ == 1 && iA == 2) theDef = theDeu    326   else if (iZ == 1 && iA == 2) theDef = theDeuteron;
345   else if (iZ == 1 && iA == 3) theDef = G4Trit    327   else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
346   else if (iZ == 2 && iA == 3) theDef = G4He3:    328   else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
347   else if (iZ == 2 && iA == 4) theDef = theAlp    329   else if (iZ == 2 && iA == 4) theDef = theAlpha;
348   else theDef = G4ParticleTable::GetParticleTa << 330   else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
349                                                   331  
350   G4double tmass = theDef->GetPDGMass();          332   G4double tmass = theDef->GetPDGMass();
351                                                   333 
352   G4LorentzVector lv(0.0,0.0,0.0,tmass);          334   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
353   lv += lv1;                                      335   lv += lv1;
354                                                   336 
355   G4ThreeVector bst = lv.boostVector();           337   G4ThreeVector bst = lv.boostVector();
356   lv1.boost(-bst);                                338   lv1.boost(-bst);
357                                                   339 
358   G4ThreeVector p1 = lv1.vect();                  340   G4ThreeVector p1 = lv1.vect();
359   G4double ptot    = p1.mag();                    341   G4double ptot    = p1.mag();
360   G4double ptot2 = ptot*ptot;                     342   G4double ptot2 = ptot*ptot;
361   G4double cost = 1 - 0.5*std::fabs(tMand)/pto    343   G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
362                                                   344 
363   if( cost >= 1.0 )      cost = 1.0;              345   if( cost >= 1.0 )      cost = 1.0;  
364   else if( cost <= -1.0) cost = -1.0;             346   else if( cost <= -1.0) cost = -1.0;
365                                                   347   
366   G4double thetaCMS = std::acos(cost);            348   G4double thetaCMS = std::acos(cost);
367                                                   349 
368   G4double sigma = GetCoulombElasticXsc( parti    350   G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
369                                                   351 
370   sigma *= pi/ptot2;                              352   sigma *= pi/ptot2;
371                                                   353 
372   return sigma;                                   354   return sigma;
373 }                                                 355 }
374                                                   356 
375 //////////////////////////////////////////////    357 ////////////////////////////////////////////////////////////////////////////
376 //                                                358 //
377 // return differential elastic probability d(p    359 // return differential elastic probability d(probability)/d(omega) 
378                                                   360 
379 G4double                                          361 G4double 
380 G4DiffuseElastic::GetDiffElasticProb( // G4Par    362 G4DiffuseElastic::GetDiffElasticProb( // G4ParticleDefinition* particle, 
381                                         G4doub    363                                         G4double theta 
382           //  G4double momentum,                  364           //  G4double momentum, 
383           // G4double A                           365           // G4double A         
384                                      )            366                                      )
385 {                                                 367 {
386   G4double sigma, bzero, bzero2, bonebyarg, bo    368   G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
387   G4double delta, diffuse, gamma;                 369   G4double delta, diffuse, gamma;
388   G4double e1, e2, bone, bone2;                   370   G4double e1, e2, bone, bone2;
389                                                   371 
390   // G4double wavek = momentum/hbarc;  // wave    372   // G4double wavek = momentum/hbarc;  // wave vector
391   // G4double r0    = 1.08*fermi;                 373   // G4double r0    = 1.08*fermi;
392   // G4double rad   = r0*G4Pow::GetInstance()- << 374   // G4double rad   = r0*std::pow(A, 1./3.);
                                                   >> 375 
                                                   >> 376   G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
                                                   >> 377   G4double kr2   = kr*kr;
                                                   >> 378   G4double krt   = kr*theta;
                                                   >> 379 
                                                   >> 380   bzero      = BesselJzero(krt);
                                                   >> 381   bzero2     = bzero*bzero;    
                                                   >> 382   bone       = BesselJone(krt);
                                                   >> 383   bone2      = bone*bone;
                                                   >> 384   bonebyarg  = BesselOneByArg(krt);
                                                   >> 385   bonebyarg2 = bonebyarg*bonebyarg;  
393                                                   386 
394   if (fParticle == theProton)                     387   if (fParticle == theProton)
395   {                                               388   {
396     diffuse = 0.63*fermi;                         389     diffuse = 0.63*fermi;
397     gamma   = 0.3*fermi;                          390     gamma   = 0.3*fermi;
398     delta   = 0.1*fermi*fermi;                    391     delta   = 0.1*fermi*fermi;
399     e1      = 0.3*fermi;                          392     e1      = 0.3*fermi;
400     e2      = 0.35*fermi;                         393     e2      = 0.35*fermi;
401   }                                               394   }
402   else if (fParticle == theNeutron)            << 
403   {                                            << 
404     diffuse =  0.63*fermi; // 1.63*fermi; //   << 
405     G4double k0 = 1*GeV/hbarc;                 << 
406     diffuse *= k0/fWaveVector;                 << 
407                                                << 
408     gamma   = 0.3*fermi;                       << 
409     delta   = 0.1*fermi*fermi;                 << 
410     e1      = 0.3*fermi;                       << 
411     e2      = 0.35*fermi;                      << 
412   }                                            << 
413   else // as proton, if were not defined          395   else // as proton, if were not defined 
414   {                                               396   {
415     diffuse = 0.63*fermi;                         397     diffuse = 0.63*fermi;
416     gamma   = 0.3*fermi;                          398     gamma   = 0.3*fermi;
417     delta   = 0.1*fermi*fermi;                    399     delta   = 0.1*fermi*fermi;
418     e1      = 0.3*fermi;                          400     e1      = 0.3*fermi;
419     e2      = 0.35*fermi;                         401     e2      = 0.35*fermi;
420   }                                               402   }
421   G4double kr    = fWaveVector*fNuclearRadius; << 
422   G4double kr2   = kr*kr;                      << 
423   G4double krt   = kr*theta;                   << 
424                                                << 
425   bzero      = BesselJzero(krt);               << 
426   bzero2     = bzero*bzero;                    << 
427   bone       = BesselJone(krt);                << 
428   bone2      = bone*bone;                      << 
429   bonebyarg  = BesselOneByArg(krt);            << 
430   bonebyarg2 = bonebyarg*bonebyarg;            << 
431                                                << 
432   G4double lambda = 15.; // 15 ok                 403   G4double lambda = 15.; // 15 ok
433                                                   404 
434   //  G4double kgamma    = fWaveVector*gamma;  << 405   //  G4double kg    = fWaveVector*gamma;   // wavek*delta;
435                                                   406 
436   G4double kgamma    = lambda*(1.-G4Exp(-fWave << 407   G4double kg    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
437   G4double kgamma2   = kgamma*kgamma;          << 408   G4double kg2   = kg*kg;
438                                                   409 
439   // G4double dk2t  = delta*fWaveVector*fWaveV    410   // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
440   // G4double dk2t2 = dk2t*dk2t;                  411   // G4double dk2t2 = dk2t*dk2t;
441   // G4double pikdt = pi*fWaveVector*diffuse*t    412   // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
442                                                   413 
443   G4double pikdt    = lambda*(1.-G4Exp(-pi*fWa << 414   G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
444                                                   415 
445   damp           = DampFactor(pikdt);             416   damp           = DampFactor(pikdt);
446   damp2          = damp*damp;                     417   damp2          = damp*damp;
447                                                   418 
448   G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector    419   G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
449   G4double e2dk3t  = -2.*e2*delta*fWaveVector*    420   G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
450                                                   421 
451                                                   422 
452   sigma  = kgamma2;                            << 423   sigma  = kg2;
453   // sigma  += dk2t2;                             424   // sigma  += dk2t2;
454   sigma *= bzero2;                                425   sigma *= bzero2;
455   sigma += mode2k2*bone2 + e2dk3t*bzero*bone;     426   sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
456   sigma += kr2*bonebyarg2;                        427   sigma += kr2*bonebyarg2;
457   sigma *= damp2;          // *rad*rad;           428   sigma *= damp2;          // *rad*rad;
458                                                   429 
459   return sigma;                                   430   return sigma;
460 }                                                 431 }
461                                                   432 
462 //////////////////////////////////////////////    433 ////////////////////////////////////////////////////////////////////////////
463 //                                                434 //
464 // return differential elastic probability d(p    435 // return differential elastic probability d(probability)/d(omega) with 
465 // Coulomb correction                             436 // Coulomb correction
466                                                   437 
467 G4double                                          438 G4double 
468 G4DiffuseElastic::GetDiffElasticSumProb( // G4    439 G4DiffuseElastic::GetDiffElasticSumProb( // G4ParticleDefinition* particle, 
469                                         G4doub    440                                         G4double theta 
470           //  G4double momentum,                  441           //  G4double momentum, 
471           // G4double A                           442           // G4double A         
472                                      )            443                                      )
473 {                                                 444 {
474   G4double sigma, bzero, bzero2, bonebyarg, bo    445   G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
475   G4double delta, diffuse, gamma;                 446   G4double delta, diffuse, gamma;
476   G4double e1, e2, bone, bone2;                   447   G4double e1, e2, bone, bone2;
477                                                   448 
478   // G4double wavek = momentum/hbarc;  // wave    449   // G4double wavek = momentum/hbarc;  // wave vector
479   // G4double r0    = 1.08*fermi;                 450   // G4double r0    = 1.08*fermi;
480   // G4double rad   = r0*G4Pow::GetInstance()- << 451   // G4double rad   = r0*std::pow(A, 1./3.);
481                                                   452 
482   G4double kr    = fWaveVector*fNuclearRadius;    453   G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
483   G4double kr2   = kr*kr;                         454   G4double kr2   = kr*kr;
484   G4double krt   = kr*theta;                      455   G4double krt   = kr*theta;
485                                                   456 
486   bzero      = BesselJzero(krt);                  457   bzero      = BesselJzero(krt);
487   bzero2     = bzero*bzero;                       458   bzero2     = bzero*bzero;    
488   bone       = BesselJone(krt);                   459   bone       = BesselJone(krt);
489   bone2      = bone*bone;                         460   bone2      = bone*bone;
490   bonebyarg  = BesselOneByArg(krt);               461   bonebyarg  = BesselOneByArg(krt);
491   bonebyarg2 = bonebyarg*bonebyarg;               462   bonebyarg2 = bonebyarg*bonebyarg;  
492                                                   463 
493   if (fParticle == theProton)                     464   if (fParticle == theProton)
494   {                                               465   {
495     diffuse = 0.63*fermi;                         466     diffuse = 0.63*fermi;
496     // diffuse = 0.6*fermi;                       467     // diffuse = 0.6*fermi;
497     gamma   = 0.3*fermi;                          468     gamma   = 0.3*fermi;
498     delta   = 0.1*fermi*fermi;                    469     delta   = 0.1*fermi*fermi;
499     e1      = 0.3*fermi;                          470     e1      = 0.3*fermi;
500     e2      = 0.35*fermi;                         471     e2      = 0.35*fermi;
501   }                                               472   }
502   else if (fParticle == theNeutron)            << 
503   {                                            << 
504     diffuse = 0.63*fermi;                      << 
505     // diffuse = 0.6*fermi;                    << 
506     G4double k0 = 1*GeV/hbarc;                 << 
507     diffuse *= k0/fWaveVector;                 << 
508     gamma   = 0.3*fermi;                       << 
509     delta   = 0.1*fermi*fermi;                 << 
510     e1      = 0.3*fermi;                       << 
511     e2      = 0.35*fermi;                      << 
512   }                                            << 
513   else // as proton, if were not defined          473   else // as proton, if were not defined 
514   {                                               474   {
515     diffuse = 0.63*fermi;                         475     diffuse = 0.63*fermi;
516     gamma   = 0.3*fermi;                          476     gamma   = 0.3*fermi;
517     delta   = 0.1*fermi*fermi;                    477     delta   = 0.1*fermi*fermi;
518     e1      = 0.3*fermi;                          478     e1      = 0.3*fermi;
519     e2      = 0.35*fermi;                         479     e2      = 0.35*fermi;
520   }                                               480   }
521   G4double lambda = 15.; // 15 ok                 481   G4double lambda = 15.; // 15 ok
522   // G4double kgamma    = fWaveVector*gamma;   << 482   // G4double kg    = fWaveVector*gamma;   // wavek*delta;
523   G4double kgamma    = lambda*(1.-G4Exp(-fWave << 483   G4double kg    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
524                                                   484 
525   // G4cout<<"kgamma = "<<kgamma<<G4endl;      << 485   // G4cout<<"kg = "<<kg<<G4endl;
526                                                   486 
527   if(fAddCoulomb)  // add Coulomb correction      487   if(fAddCoulomb)  // add Coulomb correction
528   {                                               488   {
529     G4double sinHalfTheta  = std::sin(0.5*thet    489     G4double sinHalfTheta  = std::sin(0.5*theta);
530     G4double sinHalfTheta2 = sinHalfTheta*sinH    490     G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
531                                                   491 
532     kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta << 492     kg += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
533   // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe << 493   // kg += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
534   }                                               494   }
535                                                   495 
536   G4double kgamma2   = kgamma*kgamma;          << 496   G4double kg2   = kg*kg;
537                                                   497 
538   // G4double dk2t  = delta*fWaveVector*fWaveV    498   // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
539   //   G4cout<<"dk2t = "<<dk2t<<G4endl;           499   //   G4cout<<"dk2t = "<<dk2t<<G4endl;
540   // G4double dk2t2 = dk2t*dk2t;                  500   // G4double dk2t2 = dk2t*dk2t;
541   // G4double pikdt = pi*fWaveVector*diffuse*t    501   // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
542                                                   502 
543   G4double pikdt    = lambda*(1.-G4Exp(-pi*fWa << 503   G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
544                                                   504 
545   // G4cout<<"pikdt = "<<pikdt<<G4endl;           505   // G4cout<<"pikdt = "<<pikdt<<G4endl;
546                                                   506 
547   damp           = DampFactor(pikdt);             507   damp           = DampFactor(pikdt);
548   damp2          = damp*damp;                     508   damp2          = damp*damp;
549                                                   509 
550   G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector    510   G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
551   G4double e2dk3t  = -2.*e2*delta*fWaveVector*    511   G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
552                                                   512 
553   sigma  = kgamma2;                            << 513   sigma  = kg2;
554   // sigma += dk2t2;                              514   // sigma += dk2t2;
555   sigma *= bzero2;                                515   sigma *= bzero2;
556   sigma += mode2k2*bone2;                         516   sigma += mode2k2*bone2; 
557   sigma += e2dk3t*bzero*bone;                     517   sigma += e2dk3t*bzero*bone;
558                                                   518 
559   // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf    519   // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2;  // correction at J1()/()
560   sigma += kr2*bonebyarg2;  // correction at J    520   sigma += kr2*bonebyarg2;  // correction at J1()/()
561                                                   521 
562   sigma *= damp2;          // *rad*rad;           522   sigma *= damp2;          // *rad*rad;
563                                                   523 
564   return sigma;                                   524   return sigma;
565 }                                                 525 }
566                                                   526 
567                                                   527 
568 //////////////////////////////////////////////    528 ////////////////////////////////////////////////////////////////////////////
569 //                                                529 //
570 // return differential elastic probability d(p    530 // return differential elastic probability d(probability)/d(t) with 
571 // Coulomb correction. It is called from Build << 531 // Coulomb correction
572                                                   532 
573 G4double                                          533 G4double 
574 G4DiffuseElastic::GetDiffElasticSumProbA( G4do    534 G4DiffuseElastic::GetDiffElasticSumProbA( G4double alpha )
575 {                                                 535 {
576   G4double theta;                                 536   G4double theta; 
577                                                   537 
578   theta = std::sqrt(alpha);                       538   theta = std::sqrt(alpha);
579                                                   539 
580   // theta = std::acos( 1 - alpha/2. );           540   // theta = std::acos( 1 - alpha/2. );
581                                                   541 
582   G4double sigma, bzero, bzero2, bonebyarg, bo    542   G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
583   G4double delta, diffuse, gamma;                 543   G4double delta, diffuse, gamma;
584   G4double e1, e2, bone, bone2;                   544   G4double e1, e2, bone, bone2;
585                                                   545 
586   // G4double wavek = momentum/hbarc;  // wave    546   // G4double wavek = momentum/hbarc;  // wave vector
587   // G4double r0    = 1.08*fermi;                 547   // G4double r0    = 1.08*fermi;
588   // G4double rad   = r0*G4Pow::GetInstance()- << 548   // G4double rad   = r0*std::pow(A, 1./3.);
589                                                   549 
590   G4double kr    = fWaveVector*fNuclearRadius;    550   G4double kr    = fWaveVector*fNuclearRadius; // wavek*rad;
591   G4double kr2   = kr*kr;                         551   G4double kr2   = kr*kr;
592   G4double krt   = kr*theta;                      552   G4double krt   = kr*theta;
593                                                   553 
594   bzero      = BesselJzero(krt);                  554   bzero      = BesselJzero(krt);
595   bzero2     = bzero*bzero;                       555   bzero2     = bzero*bzero;    
596   bone       = BesselJone(krt);                   556   bone       = BesselJone(krt);
597   bone2      = bone*bone;                         557   bone2      = bone*bone;
598   bonebyarg  = BesselOneByArg(krt);               558   bonebyarg  = BesselOneByArg(krt);
599   bonebyarg2 = bonebyarg*bonebyarg;               559   bonebyarg2 = bonebyarg*bonebyarg;  
600                                                   560 
601   if ( fParticle == theProton )                << 561   if (fParticle == theProton)
602   {                                            << 
603     diffuse = 0.63*fermi;                      << 
604     // diffuse = 0.6*fermi;                    << 
605     gamma   = 0.3*fermi;                       << 
606     delta   = 0.1*fermi*fermi;                 << 
607     e1      = 0.3*fermi;                       << 
608     e2      = 0.35*fermi;                      << 
609   }                                            << 
610   else if ( fParticle == theNeutron )          << 
611   {                                               562   {
612     diffuse = 0.63*fermi;                         563     diffuse = 0.63*fermi;
613     // diffuse = 0.6*fermi;                       564     // diffuse = 0.6*fermi;
614     // G4double k0 = 0.8*GeV/hbarc;            << 
615     // diffuse *= k0/fWaveVector;              << 
616     gamma   = 0.3*fermi;                          565     gamma   = 0.3*fermi;
617     delta   = 0.1*fermi*fermi;                    566     delta   = 0.1*fermi*fermi;
618     e1      = 0.3*fermi;                          567     e1      = 0.3*fermi;
619     e2      = 0.35*fermi;                         568     e2      = 0.35*fermi;
620   }                                               569   }
621   else // as proton, if were not defined          570   else // as proton, if were not defined 
622   {                                               571   {
623     diffuse = 0.63*fermi;                         572     diffuse = 0.63*fermi;
624     gamma   = 0.3*fermi;                          573     gamma   = 0.3*fermi;
625     delta   = 0.1*fermi*fermi;                    574     delta   = 0.1*fermi*fermi;
626     e1      = 0.3*fermi;                          575     e1      = 0.3*fermi;
627     e2      = 0.35*fermi;                         576     e2      = 0.35*fermi;
628   }                                               577   }
629   G4double lambda = 15.; // 15 ok                 578   G4double lambda = 15.; // 15 ok
630   // G4double kgamma    = fWaveVector*gamma;   << 579   // G4double kg    = fWaveVector*gamma;   // wavek*delta;
631   G4double kgamma    = lambda*(1.-G4Exp(-fWave << 580   G4double kg    = lambda*(1.-std::exp(-fWaveVector*gamma/lambda));   // wavek*delta;
632                                                   581 
633   // G4cout<<"kgamma = "<<kgamma<<G4endl;      << 582   // G4cout<<"kg = "<<kg<<G4endl;
634                                                   583 
635   if( fAddCoulomb )  // add Coulomb correction << 584   if(fAddCoulomb)  // add Coulomb correction
636   {                                               585   {
637     G4double sinHalfTheta  = theta*0.5; // std    586     G4double sinHalfTheta  = theta*0.5; // std::sin(0.5*theta);
638     G4double sinHalfTheta2 = sinHalfTheta*sinH    587     G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
639                                                   588 
640     kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta << 589     kg += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
641   // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe << 590   // kg += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
642   }                                               591   }
643   G4double kgamma2   = kgamma*kgamma;          << 592 
                                                   >> 593   G4double kg2   = kg*kg;
644                                                   594 
645   // G4double dk2t  = delta*fWaveVector*fWaveV    595   // G4double dk2t  = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
646   //   G4cout<<"dk2t = "<<dk2t<<G4endl;           596   //   G4cout<<"dk2t = "<<dk2t<<G4endl;
647   // G4double dk2t2 = dk2t*dk2t;                  597   // G4double dk2t2 = dk2t*dk2t;
648   // G4double pikdt = pi*fWaveVector*diffuse*t    598   // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
649                                                   599 
650   G4double pikdt    = lambda*(1. - G4Exp( -pi* << 600   G4double pikdt    = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda));   // wavek*delta;
651                                                   601 
652   // G4cout<<"pikdt = "<<pikdt<<G4endl;           602   // G4cout<<"pikdt = "<<pikdt<<G4endl;
653                                                   603 
654   damp           = DampFactor( pikdt );        << 604   damp           = DampFactor(pikdt);
655   damp2          = damp*damp;                     605   damp2          = damp*damp;
656                                                   606 
657   G4double mode2k2 = ( e1*e1 + e2*e2 )*fWaveVe << 607   G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;  
658   G4double e2dk3t  = -2.*e2*delta*fWaveVector*    608   G4double e2dk3t  = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
659                                                   609 
660   sigma  = kgamma2;                            << 610   sigma  = kg2;
661   // sigma += dk2t2;                              611   // sigma += dk2t2;
662   sigma *= bzero2;                                612   sigma *= bzero2;
663   sigma += mode2k2*bone2;                         613   sigma += mode2k2*bone2; 
664   sigma += e2dk3t*bzero*bone;                     614   sigma += e2dk3t*bzero*bone;
665                                                   615 
666   // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf    616   // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2;  // correction at J1()/()
667   sigma += kr2*bonebyarg2;  // correction at J    617   sigma += kr2*bonebyarg2;  // correction at J1()/()
668                                                   618 
669   sigma *= damp2;          // *rad*rad;           619   sigma *= damp2;          // *rad*rad;
670                                                   620 
671   return sigma;                                   621   return sigma;
672 }                                                 622 }
673                                                   623 
674                                                   624 
675 //////////////////////////////////////////////    625 ////////////////////////////////////////////////////////////////////////////
676 //                                                626 //
677 // return differential elastic probability 2*p    627 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega) 
678                                                   628 
679 G4double                                          629 G4double 
680 G4DiffuseElastic::GetIntegrandFunction( G4doub    630 G4DiffuseElastic::GetIntegrandFunction( G4double alpha )
681 {                                                 631 {
682   G4double result;                                632   G4double result;
683                                                   633 
684   result  = GetDiffElasticSumProbA(alpha);        634   result  = GetDiffElasticSumProbA(alpha);
685                                                   635 
686   // result *= 2*pi*std::sin(theta);              636   // result *= 2*pi*std::sin(theta);
687                                                   637 
688   return result;                                  638   return result;
689 }                                                 639 }
690                                                   640 
691 //////////////////////////////////////////////    641 ////////////////////////////////////////////////////////////////////////////
692 //                                                642 //
693 // return integral elastic cross section d(sig    643 // return integral elastic cross section d(sigma)/d(omega) integrated 0 - theta 
694                                                   644 
695 G4double                                          645 G4double 
696 G4DiffuseElastic::IntegralElasticProb(  const     646 G4DiffuseElastic::IntegralElasticProb(  const G4ParticleDefinition* particle, 
697                                         G4doub    647                                         G4double theta, 
698                       G4double momentum,          648                       G4double momentum, 
699                                         G4doub    649                                         G4double A         )
700 {                                                 650 {
701   G4double result;                                651   G4double result;
702   fParticle      = particle;                      652   fParticle      = particle;
703   fWaveVector    = momentum/hbarc;                653   fWaveVector    = momentum/hbarc;
704   fAtomicWeight  = A;                             654   fAtomicWeight  = A;
705                                                   655 
706   fNuclearRadius = CalculateNuclearRad(A);        656   fNuclearRadius = CalculateNuclearRad(A);
707                                                   657 
708                                                   658 
709   G4Integrator<G4DiffuseElastic,G4double(G4Dif    659   G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
710                                                   660 
711   // result = integral.Legendre10(this,&G4Diff    661   // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 
712   result = integral.Legendre96(this,&G4Diffuse    662   result = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 
713                                                   663 
714   return result;                                  664   return result;
715 }                                                 665 }
716                                                   666 
717 //////////////////////////////////////////////    667 ////////////////////////////////////////////////////////////////////////////
718 //                                                668 //
719 // Return inv momentum transfer -t > 0            669 // Return inv momentum transfer -t > 0
720                                                   670 
721 G4double G4DiffuseElastic::SampleT( const G4Pa    671 G4double G4DiffuseElastic::SampleT( const G4ParticleDefinition* aParticle, G4double p, G4double A)
722 {                                                 672 {
723   G4double theta = SampleThetaCMS( aParticle,     673   G4double theta = SampleThetaCMS( aParticle,  p, A); // sample theta in cms
724   G4double t     = 2*p*p*( 1 - std::cos(theta)    674   G4double t     = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
725   return t;                                       675   return t;
726 }                                                 676 }
727                                                   677 
728 //////////////////////////////////////////////    678 ////////////////////////////////////////////////////////////////////////////
729 //                                                679 //
730 // Return scattering angle sampled in cms         680 // Return scattering angle sampled in cms
731                                                   681 
732                                                   682 
733 G4double                                          683 G4double 
734 G4DiffuseElastic::SampleThetaCMS(const G4Parti    684 G4DiffuseElastic::SampleThetaCMS(const G4ParticleDefinition* particle, 
735                                        G4doubl    685                                        G4double momentum, G4double A)
736 {                                                 686 {
737   G4int i, iMax = 100;                            687   G4int i, iMax = 100;  
738   G4double norm, theta1, theta2, thetaMax;     << 688   G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
739   G4double result = 0., sum = 0.;              << 
740                                                   689 
741   fParticle      = particle;                      690   fParticle      = particle;
742   fWaveVector    = momentum/hbarc;                691   fWaveVector    = momentum/hbarc;
743   fAtomicWeight  = A;                             692   fAtomicWeight  = A;
744                                                   693 
745   fNuclearRadius = CalculateNuclearRad(A);        694   fNuclearRadius = CalculateNuclearRad(A);
746                                                   695   
747   thetaMax = 10.174/fWaveVector/fNuclearRadius    696   thetaMax = 10.174/fWaveVector/fNuclearRadius;
748                                                   697 
749   if (thetaMax > pi) thetaMax = pi;               698   if (thetaMax > pi) thetaMax = pi;
750                                                   699 
751   G4Integrator<G4DiffuseElastic,G4double(G4Dif    700   G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
752                                                   701 
753   // result = integral.Legendre10(this,&G4Diff    702   // result = integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, 0., theta ); 
754   norm = integral.Legendre96(this,&G4DiffuseEl    703   norm = integral.Legendre96(this,&G4DiffuseElastic::GetIntegrandFunction, 0., thetaMax );
755                                                   704 
756   norm *= G4UniformRand();                        705   norm *= G4UniformRand();
757                                                   706 
758   for(i = 1; i <= iMax; i++)                      707   for(i = 1; i <= iMax; i++)
759   {                                               708   {
760     theta1 = (i-1)*thetaMax/iMax;                 709     theta1 = (i-1)*thetaMax/iMax; 
761     theta2 = i*thetaMax/iMax;                     710     theta2 = i*thetaMax/iMax;
762     sum   += integral.Legendre10(this,&G4Diffu    711     sum   += integral.Legendre10(this,&G4DiffuseElastic::GetIntegrandFunction, theta1, theta2);
763                                                   712 
764     if ( sum >= norm )                            713     if ( sum >= norm ) 
765     {                                             714     {
766       result = 0.5*(theta1 + theta2);             715       result = 0.5*(theta1 + theta2);
767       break;                                      716       break;
768     }                                             717     }
769   }                                               718   }
770   if (i > iMax ) result = 0.5*(theta1 + theta2    719   if (i > iMax ) result = 0.5*(theta1 + theta2);
771                                                   720 
772   G4double sigma = pi*thetaMax/iMax;              721   G4double sigma = pi*thetaMax/iMax;
773                                                   722 
774   result += G4RandGauss::shoot(0.,sigma);         723   result += G4RandGauss::shoot(0.,sigma);
775                                                   724 
776   if(result < 0.) result = 0.;                    725   if(result < 0.) result = 0.;
777   if(result > thetaMax) result = thetaMax;        726   if(result > thetaMax) result = thetaMax;
778                                                   727 
779   return result;                                  728   return result;
780 }                                                 729 }
781                                                   730 
782 //////////////////////////////////////////////    731 /////////////////////////////////////////////////////////////////////////////
783 /////////////////////  Table preparation and r    732 /////////////////////  Table preparation and reading ////////////////////////
784 //////////////////////////////////////////////    733 ////////////////////////////////////////////////////////////////////////////
785 //                                                734 //
786 // Return inv momentum transfer -t > 0 from in    735 // Return inv momentum transfer -t > 0 from initialisation table
787                                                   736 
788 G4double G4DiffuseElastic::SampleInvariantT( c    737 G4double G4DiffuseElastic::SampleInvariantT( const G4ParticleDefinition* aParticle, G4double p, 
789                                                   738                                                G4int Z, G4int A)
790 {                                                 739 {
791   fParticle = aParticle;                          740   fParticle = aParticle;
792   G4double m1 = fParticle->GetPDGMass(), t;    << 741   G4double m1 = fParticle->GetPDGMass();
793   G4double totElab = std::sqrt(m1*m1+p*p);        742   G4double totElab = std::sqrt(m1*m1+p*p);
794   G4double mass2 = G4NucleiProperties::GetNucl << 743   G4double m2 = G4NucleiProperties::GetNuclearMass(A, Z);
795   G4LorentzVector lv1(p,0.0,0.0,totElab);         744   G4LorentzVector lv1(p,0.0,0.0,totElab);
796   G4LorentzVector  lv(0.0,0.0,0.0,mass2);      << 745   G4LorentzVector  lv(0.0,0.0,0.0,m2);   
797   lv += lv1;                                      746   lv += lv1;
798                                                   747 
799   G4ThreeVector bst = lv.boostVector();           748   G4ThreeVector bst = lv.boostVector();
800   lv1.boost(-bst);                                749   lv1.boost(-bst);
801                                                   750 
802   G4ThreeVector p1 = lv1.vect();                  751   G4ThreeVector p1 = lv1.vect();
803   G4double momentumCMS = p1.mag();                752   G4double momentumCMS = p1.mag();
804                                                << 
805   if( aParticle == theNeutron)                 << 
806   {                                            << 
807     G4double Tmax = NeutronTuniform( Z );      << 
808     G4double pCMS2 = momentumCMS*momentumCMS;  << 
809     G4double Tkin = std::sqrt(pCMS2+m1*m1)-m1; << 
810                                                << 
811     if( Tkin <= Tmax )                         << 
812     {                                          << 
813       t = 4.*pCMS2*G4UniformRand();            << 
814       // G4cout<<Tkin<<", "<<Tmax<<", "<<std:: << 
815       return t;                                << 
816     }                                          << 
817   }                                            << 
818                                                << 
819   t = SampleTableT( aParticle,  momentumCMS, G << 
820                                                   753 
                                                   >> 754   G4double t = SampleTableT( aParticle,  momentumCMS, G4double(Z), G4double(A) ); // sample theta2 in cms
821   return t;                                       755   return t;
822 }                                                 756 }
823                                                   757 
824 ////////////////////////////////////////////// << 
825                                                << 
826 G4double G4DiffuseElastic::NeutronTuniform(G4i << 
827 {                                              << 
828   G4double elZ  = G4double(Z);                 << 
829   elZ -= 1.;                                   << 
830   // G4double Tkin = 20.*G4Exp(-elZ/10.) + 1.; << 
831   G4double Tkin = 12.*G4Exp(-elZ/10.) + 1.;    << 
832   return Tkin;                                 << 
833 }                                              << 
834                                                << 
835                                                << 
836 //////////////////////////////////////////////    758 ////////////////////////////////////////////////////////////////////////////
837 //                                                759 //
838 // Return inv momentum transfer -t > 0 from in    760 // Return inv momentum transfer -t > 0 from initialisation table
839                                                   761 
840 G4double G4DiffuseElastic::SampleTableT( const    762 G4double G4DiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p, 
841                                                   763                                                G4double Z, G4double A)
842 {                                                 764 {
843   G4double alpha = SampleTableThetaCMS( aParti    765   G4double alpha = SampleTableThetaCMS( aParticle,  p, Z, A); // sample theta2 in cms
844   G4double t     = 2*p*p*( 1 - std::cos(std::s << 766   // G4double t     = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) );             // -t !!!
845   // G4double t     = p*p*alpha;             / << 767   G4double t     = p*p*alpha;             // -t !!!
846   return t;                                       768   return t;
847 }                                                 769 }
848                                                   770 
849 //////////////////////////////////////////////    771 ////////////////////////////////////////////////////////////////////////////
850 //                                                772 //
851 // Return scattering angle2 sampled in cms acc    773 // Return scattering angle2 sampled in cms according to precalculated table.
852                                                   774 
853                                                   775 
854 G4double                                          776 G4double 
855 G4DiffuseElastic::SampleTableThetaCMS(const G4    777 G4DiffuseElastic::SampleTableThetaCMS(const G4ParticleDefinition* particle, 
856                                        G4doubl    778                                        G4double momentum, G4double Z, G4double A)
857 {                                                 779 {
858   std::size_t iElement;                        << 780   size_t iElement;
859   G4int iMomentum, iAngle;                        781   G4int iMomentum, iAngle;  
860   G4double randAngle, position, theta1, theta2    782   G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;  
861   G4double m1 = particle->GetPDGMass();           783   G4double m1 = particle->GetPDGMass();
862                                                   784 
863   for(iElement = 0; iElement < fElementNumberV    785   for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
864   {                                               786   {
865     if( std::fabs(Z - fElementNumberVector[iEl    787     if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
866   }                                               788   }
867   if ( iElement == fElementNumberVector.size()    789   if ( iElement == fElementNumberVector.size() ) 
868   {                                               790   {
869     InitialiseOnFly(Z,A); // table preparation    791     InitialiseOnFly(Z,A); // table preparation, if needed
870                                                   792 
871     // iElement--;                                793     // iElement--;
872                                                   794 
873     // G4cout << "G4DiffuseElastic: Element wi    795     // G4cout << "G4DiffuseElastic: Element with atomic number " << Z
874     // << " is not found, return zero angle" <    796     // << " is not found, return zero angle" << G4endl;
875     // return 0.; // no table for this element    797     // return 0.; // no table for this element
876   }                                               798   }
877   // G4cout<<"iElement = "<<iElement<<G4endl;     799   // G4cout<<"iElement = "<<iElement<<G4endl;
878                                                   800 
879   fAngleTable = fAngleBank[iElement];             801   fAngleTable = fAngleBank[iElement];
880                                                   802 
881   G4double kinE = std::sqrt(momentum*momentum     803   G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
882                                                   804 
883   for( iMomentum = 0; iMomentum < fEnergyBin;     805   for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
884   {                                               806   {
885     if( kinE < fEnergyVector->GetLowEdgeEnergy    807     if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
886   }                                               808   }
887   if ( iMomentum >= fEnergyBin ) iMomentum = f    809   if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1;   // kinE is more then theMaxEnergy
888   if ( iMomentum < 0 )           iMomentum = 0    810   if ( iMomentum < 0 )           iMomentum = 0; // against negative index, kinE < theMinEnergy
889                                                   811 
890   // G4cout<<"iMomentum = "<<iMomentum<<G4endl    812   // G4cout<<"iMomentum = "<<iMomentum<<G4endl;
891                                                   813 
892   if (iMomentum == fEnergyBin -1 || iMomentum     814   if (iMomentum == fEnergyBin -1 || iMomentum == 0 )   // the table edges
893   {                                               815   {
894     position = (*(*fAngleTable)(iMomentum))(fA    816     position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
895                                                   817 
896     // G4cout<<"position = "<<position<<G4endl    818     // G4cout<<"position = "<<position<<G4endl;
897                                                   819 
898     for(iAngle = 0; iAngle < fAngleBin-1; iAng    820     for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
899     {                                             821     {
900       if( position > (*(*fAngleTable)(iMomentu << 822       if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
901     }                                             823     }
902     if (iAngle >= fAngleBin-1) iAngle = fAngle    824     if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
903                                                   825 
904     // G4cout<<"iAngle = "<<iAngle<<G4endl;       826     // G4cout<<"iAngle = "<<iAngle<<G4endl;
905                                                   827 
906     randAngle = GetScatteringAngle(iMomentum,     828     randAngle = GetScatteringAngle(iMomentum, iAngle, position);
907                                                   829 
908     // G4cout<<"randAngle = "<<randAngle<<G4en    830     // G4cout<<"randAngle = "<<randAngle<<G4endl;
909   }                                               831   }
910   else  // kinE inside between energy table ed    832   else  // kinE inside between energy table edges
911   {                                               833   {
912     // position = (*(*fAngleTable)(iMomentum))    834     // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
913     position = (*(*fAngleTable)(iMomentum))(0)    835     position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
914                                                   836 
915     // G4cout<<"position = "<<position<<G4endl    837     // G4cout<<"position = "<<position<<G4endl;
916                                                   838 
917     for(iAngle = 0; iAngle < fAngleBin-1; iAng    839     for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
918     {                                             840     {
919       // if( position < (*(*fAngleTable)(iMome    841       // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
920       if( position > (*(*fAngleTable)(iMomentu    842       if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
921     }                                             843     }
922     if (iAngle >= fAngleBin-1) iAngle = fAngle    844     if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
923                                                   845 
924     // G4cout<<"iAngle = "<<iAngle<<G4endl;       846     // G4cout<<"iAngle = "<<iAngle<<G4endl;
925                                                   847 
926     theta2  = GetScatteringAngle(iMomentum, iA    848     theta2  = GetScatteringAngle(iMomentum, iAngle, position);
927                                                   849 
928     // G4cout<<"theta2 = "<<theta2<<G4endl;       850     // G4cout<<"theta2 = "<<theta2<<G4endl;
929     E2 = fEnergyVector->GetLowEdgeEnergy(iMome    851     E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
930                                                   852 
931     // G4cout<<"E2 = "<<E2<<G4endl;               853     // G4cout<<"E2 = "<<E2<<G4endl;
932                                                   854     
933     iMomentum--;                                  855     iMomentum--;
934                                                   856     
935     // position = (*(*fAngleTable)(iMomentum))    857     // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
936                                                   858 
937     // G4cout<<"position = "<<position<<G4endl    859     // G4cout<<"position = "<<position<<G4endl;
938                                                   860 
939     for(iAngle = 0; iAngle < fAngleBin-1; iAng    861     for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
940     {                                             862     {
941       // if( position < (*(*fAngleTable)(iMome    863       // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
942       if( position > (*(*fAngleTable)(iMomentu    864       if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
943     }                                             865     }
944     if (iAngle >= fAngleBin-1) iAngle = fAngle    866     if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
945                                                   867     
946     theta1  = GetScatteringAngle(iMomentum, iA    868     theta1  = GetScatteringAngle(iMomentum, iAngle, position);
947                                                   869 
948     // G4cout<<"theta1 = "<<theta1<<G4endl;       870     // G4cout<<"theta1 = "<<theta1<<G4endl;
949                                                   871 
950     E1 = fEnergyVector->GetLowEdgeEnergy(iMome    872     E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
951                                                   873 
952     // G4cout<<"E1 = "<<E1<<G4endl;               874     // G4cout<<"E1 = "<<E1<<G4endl;
953                                                   875 
954     W  = 1.0/(E2 - E1);                           876     W  = 1.0/(E2 - E1);
955     W1 = (E2 - kinE)*W;                           877     W1 = (E2 - kinE)*W;
956     W2 = (kinE - E1)*W;                           878     W2 = (kinE - E1)*W;
957                                                   879 
958     randAngle = W1*theta1 + W2*theta2;            880     randAngle = W1*theta1 + W2*theta2;
959                                                   881     
960     // randAngle = theta2;                        882     // randAngle = theta2;
961     // G4cout<<"randAngle = "<<randAngle<<G4en    883     // G4cout<<"randAngle = "<<randAngle<<G4endl;
962   }                                               884   }
963   // G4double angle = randAngle;                  885   // G4double angle = randAngle;
964   // if (randAngle > 0.) randAngle /= 2*pi*std    886   // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
965                                                   887 
966   if(randAngle < 0.) randAngle = 0.;           << 
967                                                << 
968   return randAngle;                               888   return randAngle;
969 }                                                 889 }
970                                                   890 
971 //////////////////////////////////////////////    891 //////////////////////////////////////////////////////////////////////////////
972 //                                                892 //
973 // Initialisation for given particle on fly us    893 // Initialisation for given particle on fly using new element number
974                                                   894 
975 void G4DiffuseElastic::InitialiseOnFly(G4doubl    895 void G4DiffuseElastic::InitialiseOnFly(G4double Z, G4double A) 
976 {                                                 896 {
977   fAtomicNumber  = Z;     // atomic number        897   fAtomicNumber  = Z;     // atomic number
978   fAtomicWeight  = G4NistManager::Instance()-> << 898   fAtomicWeight  = A;     // number of nucleons
979                                                   899 
980   fNuclearRadius = CalculateNuclearRad(fAtomic    900   fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
981                                                << 901   
982   if( verboseLevel > 0 )                          902   if( verboseLevel > 0 )    
983   {                                               903   {
984     G4cout<<"G4DiffuseElastic::InitialiseOnFly << 904     G4cout<<"G4DiffuseElastic::Initialise() the element with Z = "
985     <<Z<<"; and A = "<<A<<G4endl;                 905     <<Z<<"; and A = "<<A<<G4endl;
986   }                                               906   }
987   fElementNumberVector.push_back(fAtomicNumber    907   fElementNumberVector.push_back(fAtomicNumber);
988                                                   908 
989   BuildAngleTable();                              909   BuildAngleTable();
990                                                   910 
991   fAngleBank.push_back(fAngleTable);              911   fAngleBank.push_back(fAngleTable);
992                                                   912 
993   return;                                         913   return;
994 }                                                 914 }
995                                                   915 
996 //////////////////////////////////////////////    916 ///////////////////////////////////////////////////////////////////////////////
997 //                                                917 //
998 // Build for given particle and element table     918 // Build for given particle and element table of momentum, angle probability.
999 // For the moment in lab system.                  919 // For the moment in lab system. 
1000                                                  920 
1001 void G4DiffuseElastic::BuildAngleTable()         921 void G4DiffuseElastic::BuildAngleTable() 
1002 {                                                922 {
1003   G4int i, j;                                    923   G4int i, j;
1004   G4double partMom, kinE, a = 0., z = fPartic    924   G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1005   G4double alpha1, alpha2, alphaMax, alphaCou    925   G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1006                                                  926 
1007   G4Integrator<G4DiffuseElastic,G4double(G4Di    927   G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
1008                                                  928   
1009   fAngleTable = new G4PhysicsTable( fEnergyBi << 929   fAngleTable = new G4PhysicsTable(fEnergyBin);
1010                                                  930 
1011   for( i = 0; i < fEnergyBin; i++)               931   for( i = 0; i < fEnergyBin; i++)
1012   {                                              932   {
1013     kinE        = fEnergyVector->GetLowEdgeEn    933     kinE        = fEnergyVector->GetLowEdgeEnergy(i);
1014     partMom     = std::sqrt( kinE*(kinE + 2*m    934     partMom     = std::sqrt( kinE*(kinE + 2*m1) );
1015                                                  935 
1016     fWaveVector = partMom/hbarc;                 936     fWaveVector = partMom/hbarc;
1017                                                  937 
1018     G4double kR     = fWaveVector*fNuclearRad    938     G4double kR     = fWaveVector*fNuclearRadius;
1019     G4double kR2    = kR*kR;                     939     G4double kR2    = kR*kR;
1020     G4double kRmax  = 18.6; // 10.6; 10.6, 18    940     G4double kRmax  = 18.6; // 10.6; 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1021     G4double kRcoul = 1.9; // 1.2; 1.4, 2.5;     941     G4double kRcoul = 1.9; // 1.2; 1.4, 2.5; // on the first slope of J1
1022     // G4double kRlim  = 1.2;                    942     // G4double kRlim  = 1.2;
1023     // G4double kRlim2 = kRlim*kRlim/kR2;        943     // G4double kRlim2 = kRlim*kRlim/kR2;
1024                                                  944 
1025     alphaMax = kRmax*kRmax/kR2;                  945     alphaMax = kRmax*kRmax/kR2;
1026                                                  946 
1027                                               << 947     if (alphaMax > 4.) alphaMax = 4.;  // vmg05-02-09: was pi2 
1028     // if (alphaMax > 4.) alphaMax = 4.;  //  << 
1029     // if ( alphaMax > 4. || alphaMax < 1. )  << 
1030                                               << 
1031     // if ( alphaMax > 4. || alphaMax < 1. )  << 
1032                                               << 
1033     // G4cout<<"alphaMax = "<<alphaMax<<", "; << 
1034                                               << 
1035     if ( alphaMax >= CLHEP::pi*CLHEP::pi ) al << 
1036                                                  948 
1037     alphaCoulomb = kRcoul*kRcoul/kR2;            949     alphaCoulomb = kRcoul*kRcoul/kR2;
1038                                                  950 
1039     if( z )                                      951     if( z )
1040     {                                            952     {
1041       a           = partMom/m1;         // be    953       a           = partMom/m1;         // beta*gamma for m1
1042       fBeta       = a/std::sqrt(1+a*a);          954       fBeta       = a/std::sqrt(1+a*a);
1043       fZommerfeld = CalculateZommerfeld( fBet    955       fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1044       fAm         = CalculateAm( partMom, fZo    956       fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1045     }                                            957     }
1046     G4PhysicsFreeVector* angleVector = new G4    958     G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1047                                                  959 
1048     // G4PhysicsLogVector*  angleBins = new G    960     // G4PhysicsLogVector*  angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1049                                                  961 
1050     G4double delth = alphaMax/fAngleBin;         962     G4double delth = alphaMax/fAngleBin;
1051                                                  963         
1052     sum = 0.;                                    964     sum = 0.;
1053                                                  965 
1054     // fAddCoulomb = false;                      966     // fAddCoulomb = false;
1055     fAddCoulomb = true;                          967     fAddCoulomb = true;
1056                                                  968 
1057     // for(j = 1; j < fAngleBin; j++)            969     // for(j = 1; j < fAngleBin; j++)
1058     for(j = fAngleBin-1; j >= 1; j--)            970     for(j = fAngleBin-1; j >= 1; j--)
1059     {                                            971     {
1060       // alpha1 = angleBins->GetLowEdgeEnergy    972       // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1061       // alpha2 = angleBins->GetLowEdgeEnergy    973       // alpha2 = angleBins->GetLowEdgeEnergy(j);
1062                                                  974 
1063       // alpha1 = alphaMax*(j-1)/fAngleBin;      975       // alpha1 = alphaMax*(j-1)/fAngleBin;
1064       // alpha2 = alphaMax*( j )/fAngleBin;      976       // alpha2 = alphaMax*( j )/fAngleBin;
1065                                                  977 
1066       alpha1 = delth*(j-1);                      978       alpha1 = delth*(j-1);
1067       // if(alpha1 < kRlim2) alpha1 = kRlim2;    979       // if(alpha1 < kRlim2) alpha1 = kRlim2;
1068       alpha2 = alpha1 + delth;                   980       alpha2 = alpha1 + delth;
1069                                                  981 
1070       // if( ( alpha2 > alphaCoulomb ) && z )    982       // if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1071       if( ( alpha1 < alphaCoulomb ) && z ) fA    983       if( ( alpha1 < alphaCoulomb ) && z ) fAddCoulomb = false;
1072                                                  984 
1073       delta = integral.Legendre10(this, &G4Di    985       delta = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1074       // delta = integral.Legendre96(this, &G    986       // delta = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1075                                                  987 
1076       sum += delta;                              988       sum += delta;
1077                                                  989       
1078       angleVector->PutValue( j-1 , alpha1, su    990       angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
1079       //      G4cout<<"j-1 = "<<j-1<<"; alpha << 991       // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2<<"; sum = "<<sum<<G4endl;
1080     }                                            992     }
1081     fAngleTable->insertAt(i, angleVector);    << 993     fAngleTable->insertAt(i,angleVector);
1082                                                  994 
1083     // delete[] angleVector;                     995     // delete[] angleVector; 
1084     // delete[] angleBins;                       996     // delete[] angleBins; 
1085   }                                              997   }
1086   return;                                        998   return;
1087 }                                                999 }
1088                                                  1000 
1089 /////////////////////////////////////////////    1001 /////////////////////////////////////////////////////////////////////////////////
1090 //                                               1002 //
1091 //                                               1003 //
1092                                                  1004 
1093 G4double                                         1005 G4double 
1094 G4DiffuseElastic:: GetScatteringAngle( G4int     1006 G4DiffuseElastic:: GetScatteringAngle( G4int iMomentum, G4int iAngle, G4double position )
1095 {                                                1007 {
1096  G4double x1, x2, y1, y2, randAngle;             1008  G4double x1, x2, y1, y2, randAngle;
1097                                                  1009 
1098   if( iAngle == 0 )                              1010   if( iAngle == 0 )
1099   {                                              1011   {
1100     randAngle = (*fAngleTable)(iMomentum)->Ge    1012     randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1101     // iAngle++;                                 1013     // iAngle++;
1102   }                                              1014   }
1103   else                                           1015   else
1104   {                                              1016   {
1105     if ( iAngle >= G4int((*fAngleTable)(iMome    1017     if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1106     {                                            1018     {
1107       iAngle = G4int((*fAngleTable)(iMomentum << 1019       iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1108     }                                            1020     }
1109     y1 = (*(*fAngleTable)(iMomentum))(iAngle-    1021     y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1110     y2 = (*(*fAngleTable)(iMomentum))(iAngle)    1022     y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1111                                                  1023 
1112     x1 = (*fAngleTable)(iMomentum)->GetLowEdg    1024     x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1113     x2 = (*fAngleTable)(iMomentum)->GetLowEdg    1025     x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1114                                                  1026 
1115     if ( x1 == x2 )   randAngle = x2;            1027     if ( x1 == x2 )   randAngle = x2;
1116     else                                         1028     else
1117     {                                            1029     {
1118       if ( y1 == y2 ) randAngle = x1 + ( x2 -    1030       if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
1119       else                                       1031       else
1120       {                                          1032       {
1121         randAngle = x1 + ( position - y1 )*(     1033         randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1122       }                                          1034       }
1123     }                                            1035     }
1124   }                                              1036   }
1125   return randAngle;                              1037   return randAngle;
1126 }                                                1038 }
1127                                                  1039 
1128                                                  1040 
1129                                                  1041 
1130 /////////////////////////////////////////////    1042 ////////////////////////////////////////////////////////////////////////////
1131 //                                               1043 //
1132 // Return scattering angle sampled in lab sys    1044 // Return scattering angle sampled in lab system (target at rest)
1133                                                  1045 
1134                                                  1046 
1135                                                  1047 
1136 G4double                                         1048 G4double 
1137 G4DiffuseElastic::SampleThetaLab( const G4Had    1049 G4DiffuseElastic::SampleThetaLab( const G4HadProjectile* aParticle, 
1138                                         G4dou    1050                                         G4double tmass, G4double A)
1139 {                                                1051 {
1140   const G4ParticleDefinition* theParticle = a    1052   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1141   G4double m1 = theParticle->GetPDGMass();       1053   G4double m1 = theParticle->GetPDGMass();
1142   G4double plab = aParticle->GetTotalMomentum    1054   G4double plab = aParticle->GetTotalMomentum();
1143   G4LorentzVector lv1 = aParticle->Get4Moment    1055   G4LorentzVector lv1 = aParticle->Get4Momentum();
1144   G4LorentzVector lv(0.0,0.0,0.0,tmass);         1056   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
1145   lv += lv1;                                     1057   lv += lv1;
1146                                                  1058 
1147   G4ThreeVector bst = lv.boostVector();          1059   G4ThreeVector bst = lv.boostVector();
1148   lv1.boost(-bst);                               1060   lv1.boost(-bst);
1149                                                  1061 
1150   G4ThreeVector p1 = lv1.vect();                 1062   G4ThreeVector p1 = lv1.vect();
1151   G4double ptot    = p1.mag();                   1063   G4double ptot    = p1.mag();
1152   G4double tmax    = 4.0*ptot*ptot;              1064   G4double tmax    = 4.0*ptot*ptot;
1153   G4double t       = 0.0;                        1065   G4double t       = 0.0;
1154                                                  1066 
1155                                                  1067 
1156   //                                             1068   //
1157   // Sample t                                    1069   // Sample t
1158   //                                             1070   //
1159                                                  1071   
1160   t = SampleT( theParticle, ptot, A);            1072   t = SampleT( theParticle, ptot, A);
1161                                                  1073 
1162   // NaN finder                                  1074   // NaN finder
1163   if(!(t < 0.0 || t >= 0.0))                     1075   if(!(t < 0.0 || t >= 0.0)) 
1164   {                                              1076   {
1165     if (verboseLevel > 0)                        1077     if (verboseLevel > 0) 
1166     {                                            1078     {
1167       G4cout << "G4DiffuseElastic:WARNING: A     1079       G4cout << "G4DiffuseElastic:WARNING: A = " << A 
1168        << " mom(GeV)= " << plab/GeV              1080        << " mom(GeV)= " << plab/GeV 
1169              << " S-wave will be sampled"        1081              << " S-wave will be sampled" 
1170        << G4endl;                                1082        << G4endl; 
1171     }                                            1083     }
1172     t = G4UniformRand()*tmax;                    1084     t = G4UniformRand()*tmax; 
1173   }                                              1085   }
1174   if(verboseLevel>1)                             1086   if(verboseLevel>1)
1175   {                                              1087   {
1176     G4cout <<" t= " << t << " tmax= " << tmax    1088     G4cout <<" t= " << t << " tmax= " << tmax 
1177      << " ptot= " << ptot << G4endl;             1089      << " ptot= " << ptot << G4endl;
1178   }                                              1090   }
1179   // Sampling of angles in CM system             1091   // Sampling of angles in CM system
1180                                                  1092 
1181   G4double phi  = G4UniformRand()*twopi;         1093   G4double phi  = G4UniformRand()*twopi;
1182   G4double cost = 1. - 2.0*t/tmax;               1094   G4double cost = 1. - 2.0*t/tmax;
1183   G4double sint;                                 1095   G4double sint;
1184                                                  1096 
1185   if( cost >= 1.0 )                              1097   if( cost >= 1.0 ) 
1186   {                                              1098   {
1187     cost = 1.0;                                  1099     cost = 1.0;
1188     sint = 0.0;                                  1100     sint = 0.0;
1189   }                                              1101   }
1190   else if( cost <= -1.0)                         1102   else if( cost <= -1.0) 
1191   {                                              1103   {
1192     cost = -1.0;                                 1104     cost = -1.0;
1193     sint =  0.0;                                 1105     sint =  0.0;
1194   }                                              1106   }
1195   else                                           1107   else  
1196   {                                              1108   {
1197     sint = std::sqrt((1.0-cost)*(1.0+cost));     1109     sint = std::sqrt((1.0-cost)*(1.0+cost));
1198   }                                              1110   }    
1199   if (verboseLevel>1)                            1111   if (verboseLevel>1) 
1200   {                                              1112   {
1201     G4cout << "cos(t)=" << cost << " std::sin    1113     G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
1202   }                                              1114   }
1203   G4ThreeVector v1(sint*std::cos(phi),sint*st    1115   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1204   v1 *= ptot;                                    1116   v1 *= ptot;
1205   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s    1117   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1206                                                  1118 
1207   nlv1.boost(bst);                               1119   nlv1.boost(bst); 
1208                                                  1120 
1209   G4ThreeVector np1 = nlv1.vect();               1121   G4ThreeVector np1 = nlv1.vect();
1210                                                  1122 
1211     // G4double theta = std::acos( np1.z()/np    1123     // G4double theta = std::acos( np1.z()/np1.mag() );  // degree;
1212                                                  1124 
1213   G4double theta = np1.theta();                  1125   G4double theta = np1.theta();
1214                                                  1126 
1215   return theta;                                  1127   return theta;
1216 }                                                1128 }
1217                                                  1129 
1218 /////////////////////////////////////////////    1130 ////////////////////////////////////////////////////////////////////////////
1219 //                                               1131 //
1220 // Return scattering angle in lab system (tar    1132 // Return scattering angle in lab system (target at rest) knowing theta in CMS
1221                                                  1133 
1222                                                  1134 
1223                                                  1135 
1224 G4double                                         1136 G4double 
1225 G4DiffuseElastic::ThetaCMStoThetaLab( const G    1137 G4DiffuseElastic::ThetaCMStoThetaLab( const G4DynamicParticle* aParticle, 
1226                                         G4dou    1138                                         G4double tmass, G4double thetaCMS)
1227 {                                                1139 {
1228   const G4ParticleDefinition* theParticle = a    1140   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1229   G4double m1 = theParticle->GetPDGMass();       1141   G4double m1 = theParticle->GetPDGMass();
1230   // G4double plab = aParticle->GetTotalMomen    1142   // G4double plab = aParticle->GetTotalMomentum();
1231   G4LorentzVector lv1 = aParticle->Get4Moment    1143   G4LorentzVector lv1 = aParticle->Get4Momentum();
1232   G4LorentzVector lv(0.0,0.0,0.0,tmass);         1144   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
1233                                                  1145 
1234   lv += lv1;                                     1146   lv += lv1;
1235                                                  1147 
1236   G4ThreeVector bst = lv.boostVector();          1148   G4ThreeVector bst = lv.boostVector();
1237                                                  1149 
1238   lv1.boost(-bst);                               1150   lv1.boost(-bst);
1239                                                  1151 
1240   G4ThreeVector p1 = lv1.vect();                 1152   G4ThreeVector p1 = lv1.vect();
1241   G4double ptot    = p1.mag();                   1153   G4double ptot    = p1.mag();
1242                                                  1154 
1243   G4double phi  = G4UniformRand()*twopi;         1155   G4double phi  = G4UniformRand()*twopi;
1244   G4double cost = std::cos(thetaCMS);            1156   G4double cost = std::cos(thetaCMS);
1245   G4double sint;                                 1157   G4double sint;
1246                                                  1158 
1247   if( cost >= 1.0 )                              1159   if( cost >= 1.0 ) 
1248   {                                              1160   {
1249     cost = 1.0;                                  1161     cost = 1.0;
1250     sint = 0.0;                                  1162     sint = 0.0;
1251   }                                              1163   }
1252   else if( cost <= -1.0)                         1164   else if( cost <= -1.0) 
1253   {                                              1165   {
1254     cost = -1.0;                                 1166     cost = -1.0;
1255     sint =  0.0;                                 1167     sint =  0.0;
1256   }                                              1168   }
1257   else                                           1169   else  
1258   {                                              1170   {
1259     sint = std::sqrt((1.0-cost)*(1.0+cost));     1171     sint = std::sqrt((1.0-cost)*(1.0+cost));
1260   }                                              1172   }    
1261   if (verboseLevel>1)                            1173   if (verboseLevel>1) 
1262   {                                              1174   {
1263     G4cout << "cos(tcms)=" << cost << " std::    1175     G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
1264   }                                              1176   }
1265   G4ThreeVector v1(sint*std::cos(phi),sint*st    1177   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1266   v1 *= ptot;                                    1178   v1 *= ptot;
1267   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s    1179   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1268                                                  1180 
1269   nlv1.boost(bst);                               1181   nlv1.boost(bst); 
1270                                                  1182 
1271   G4ThreeVector np1 = nlv1.vect();               1183   G4ThreeVector np1 = nlv1.vect();
1272                                                  1184 
1273                                                  1185 
1274   G4double thetaLab = np1.theta();               1186   G4double thetaLab = np1.theta();
1275                                                  1187 
1276   return thetaLab;                               1188   return thetaLab;
1277 }                                                1189 }
1278 /////////////////////////////////////////////    1190 ////////////////////////////////////////////////////////////////////////////
1279 //                                               1191 //
1280 // Return scattering angle in CMS system (tar    1192 // Return scattering angle in CMS system (target at rest) knowing theta in Lab
1281                                                  1193 
1282                                                  1194 
1283                                                  1195 
1284 G4double                                         1196 G4double 
1285 G4DiffuseElastic::ThetaLabToThetaCMS( const G    1197 G4DiffuseElastic::ThetaLabToThetaCMS( const G4DynamicParticle* aParticle, 
1286                                         G4dou    1198                                         G4double tmass, G4double thetaLab)
1287 {                                                1199 {
1288   const G4ParticleDefinition* theParticle = a    1200   const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1289   G4double m1 = theParticle->GetPDGMass();       1201   G4double m1 = theParticle->GetPDGMass();
1290   G4double plab = aParticle->GetTotalMomentum    1202   G4double plab = aParticle->GetTotalMomentum();
1291   G4LorentzVector lv1 = aParticle->Get4Moment    1203   G4LorentzVector lv1 = aParticle->Get4Momentum();
1292   G4LorentzVector lv(0.0,0.0,0.0,tmass);         1204   G4LorentzVector lv(0.0,0.0,0.0,tmass);   
1293                                                  1205 
1294   lv += lv1;                                     1206   lv += lv1;
1295                                                  1207 
1296   G4ThreeVector bst = lv.boostVector();          1208   G4ThreeVector bst = lv.boostVector();
1297                                                  1209 
1298   // lv1.boost(-bst);                            1210   // lv1.boost(-bst);
1299                                                  1211 
1300   // G4ThreeVector p1 = lv1.vect();              1212   // G4ThreeVector p1 = lv1.vect();
1301   // G4double ptot    = p1.mag();                1213   // G4double ptot    = p1.mag();
1302                                                  1214 
1303   G4double phi  = G4UniformRand()*twopi;         1215   G4double phi  = G4UniformRand()*twopi;
1304   G4double cost = std::cos(thetaLab);            1216   G4double cost = std::cos(thetaLab);
1305   G4double sint;                                 1217   G4double sint;
1306                                                  1218 
1307   if( cost >= 1.0 )                              1219   if( cost >= 1.0 ) 
1308   {                                              1220   {
1309     cost = 1.0;                                  1221     cost = 1.0;
1310     sint = 0.0;                                  1222     sint = 0.0;
1311   }                                              1223   }
1312   else if( cost <= -1.0)                         1224   else if( cost <= -1.0) 
1313   {                                              1225   {
1314     cost = -1.0;                                 1226     cost = -1.0;
1315     sint =  0.0;                                 1227     sint =  0.0;
1316   }                                              1228   }
1317   else                                           1229   else  
1318   {                                              1230   {
1319     sint = std::sqrt((1.0-cost)*(1.0+cost));     1231     sint = std::sqrt((1.0-cost)*(1.0+cost));
1320   }                                              1232   }    
1321   if (verboseLevel>1)                            1233   if (verboseLevel>1) 
1322   {                                              1234   {
1323     G4cout << "cos(tlab)=" << cost << " std::    1235     G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
1324   }                                              1236   }
1325   G4ThreeVector v1(sint*std::cos(phi),sint*st    1237   G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1326   v1 *= plab;                                    1238   v1 *= plab;
1327   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s    1239   G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1328                                                  1240 
1329   nlv1.boost(-bst);                              1241   nlv1.boost(-bst); 
1330                                                  1242 
1331   G4ThreeVector np1 = nlv1.vect();               1243   G4ThreeVector np1 = nlv1.vect();
1332                                                  1244 
1333                                                  1245 
1334   G4double thetaCMS = np1.theta();               1246   G4double thetaCMS = np1.theta();
1335                                                  1247 
1336   return thetaCMS;                               1248   return thetaCMS;
1337 }                                                1249 }
1338                                                  1250 
1339 /////////////////////////////////////////////    1251 ///////////////////////////////////////////////////////////////////////////////
1340 //                                               1252 //
1341 // Test for given particle and element table     1253 // Test for given particle and element table of momentum, angle probability.
1342 // For the moment in lab system.                 1254 // For the moment in lab system. 
1343                                                  1255 
1344 void G4DiffuseElastic::TestAngleTable(const G    1256 void G4DiffuseElastic::TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
1345                                             G    1257                                             G4double Z, G4double A) 
1346 {                                                1258 {
1347   fAtomicNumber  = Z;     // atomic number       1259   fAtomicNumber  = Z;     // atomic number
1348   fAtomicWeight  = A;     // number of nucleo    1260   fAtomicWeight  = A;     // number of nucleons
1349   fNuclearRadius = CalculateNuclearRad(fAtomi    1261   fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1350                                                  1262   
1351                                                  1263      
1352                                                  1264   
1353   G4cout<<"G4DiffuseElastic::TestAngleTable()    1265   G4cout<<"G4DiffuseElastic::TestAngleTable() init the element with Z = "
1354     <<Z<<"; and A = "<<A<<G4endl;                1266     <<Z<<"; and A = "<<A<<G4endl;
1355                                                  1267  
1356   fElementNumberVector.push_back(fAtomicNumbe    1268   fElementNumberVector.push_back(fAtomicNumber);
1357                                                  1269 
1358                                                  1270  
1359                                                  1271 
1360                                                  1272 
1361   G4int i=0, j;                                  1273   G4int i=0, j;
1362   G4double a = 0., z = theParticle->GetPDGCha    1274   G4double a = 0., z = theParticle->GetPDGCharge(),  m1 = fParticle->GetPDGMass();
1363   G4double alpha1=0., alpha2=0., alphaMax=0.,    1275   G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1364   G4double deltaL10 = 0., deltaL96 = 0., delt    1276   G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1365   G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.    1277   G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1366   G4double epsilon = 0.001;                      1278   G4double epsilon = 0.001;
1367                                                  1279 
1368   G4Integrator<G4DiffuseElastic,G4double(G4Di    1280   G4Integrator<G4DiffuseElastic,G4double(G4DiffuseElastic::*)(G4double)> integral;
1369                                                  1281   
1370   fAngleTable = new G4PhysicsTable(fEnergyBin    1282   fAngleTable = new G4PhysicsTable(fEnergyBin);
1371                                                  1283 
1372   fWaveVector = partMom/hbarc;                   1284   fWaveVector = partMom/hbarc;
1373                                                  1285 
1374   G4double kR     = fWaveVector*fNuclearRadiu    1286   G4double kR     = fWaveVector*fNuclearRadius;
1375   G4double kR2    = kR*kR;                       1287   G4double kR2    = kR*kR;
1376   G4double kRmax  = 10.6; // 10.6, 18, 10.174    1288   G4double kRmax  = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1377   G4double kRcoul = 1.2; // 1.4, 2.5; // on t    1289   G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
1378                                                  1290 
1379   alphaMax = kRmax*kRmax/kR2;                    1291   alphaMax = kRmax*kRmax/kR2;
1380                                                  1292 
1381   if (alphaMax > 4.) alphaMax = 4.;  // vmg05    1293   if (alphaMax > 4.) alphaMax = 4.;  // vmg05-02-09: was pi2 
1382                                                  1294 
1383   alphaCoulomb = kRcoul*kRcoul/kR2;              1295   alphaCoulomb = kRcoul*kRcoul/kR2;
1384                                                  1296 
1385   if( z )                                        1297   if( z )
1386   {                                              1298   {
1387       a           = partMom/m1; // beta*gamma    1299       a           = partMom/m1; // beta*gamma for m1
1388       fBeta       = a/std::sqrt(1+a*a);          1300       fBeta       = a/std::sqrt(1+a*a);
1389       fZommerfeld = CalculateZommerfeld( fBet    1301       fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1390       fAm         = CalculateAm( partMom, fZo    1302       fAm         = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1391   }                                              1303   }
1392   G4PhysicsFreeVector* angleVector = new G4Ph    1304   G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1393                                                  1305 
1394   // G4PhysicsLogVector*  angleBins = new G4P    1306   // G4PhysicsLogVector*  angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1395                                                  1307     
1396                                                  1308   
1397   fAddCoulomb = false;                           1309   fAddCoulomb = false;
1398                                                  1310 
1399   for(j = 1; j < fAngleBin; j++)                 1311   for(j = 1; j < fAngleBin; j++)
1400   {                                              1312   {
1401       // alpha1 = angleBins->GetLowEdgeEnergy    1313       // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1402       // alpha2 = angleBins->GetLowEdgeEnergy    1314       // alpha2 = angleBins->GetLowEdgeEnergy(j);
1403                                                  1315 
1404     alpha1 = alphaMax*(j-1)/fAngleBin;           1316     alpha1 = alphaMax*(j-1)/fAngleBin;
1405     alpha2 = alphaMax*( j )/fAngleBin;           1317     alpha2 = alphaMax*( j )/fAngleBin;
1406                                                  1318 
1407     if( ( alpha2 > alphaCoulomb ) && z ) fAdd    1319     if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1408                                                  1320 
1409     deltaL10 = integral.Legendre10(this, &G4D    1321     deltaL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1410     deltaL96 = integral.Legendre96(this, &G4D    1322     deltaL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1411     deltaAG  = integral.AdaptiveGauss(this, &    1323     deltaAG  = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction, 
1412                                        alpha1    1324                                        alpha1, alpha2,epsilon);
1413                                                  1325 
1414       // G4cout<<alpha1<<"\t"<<std::sqrt(alph    1326       // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1415       //     <<deltaL10<<"\t"<<deltaL96<<"\t"    1327       //     <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1416                                                  1328 
1417     sumL10 += deltaL10;                          1329     sumL10 += deltaL10;
1418     sumL96 += deltaL96;                          1330     sumL96 += deltaL96;
1419     sumAG  += deltaAG;                           1331     sumAG  += deltaAG;
1420                                                  1332 
1421     G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/d    1333     G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1422             <<sumL10<<"\t"<<sumL96<<"\t"<<sum    1334             <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1423                                                  1335       
1424     angleVector->PutValue( j-1 , alpha1, sumL    1336     angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
1425   }                                              1337   }
1426   fAngleTable->insertAt(i,angleVector);          1338   fAngleTable->insertAt(i,angleVector);
1427   fAngleBank.push_back(fAngleTable);             1339   fAngleBank.push_back(fAngleTable);
1428                                                  1340 
1429   /*                                             1341   /*
1430   // Integral over all angle range - Bad accu    1342   // Integral over all angle range - Bad accuracy !!!
1431                                                  1343 
1432   sumL10 = integral.Legendre10(this, &G4Diffu    1344   sumL10 = integral.Legendre10(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
1433   sumL96 = integral.Legendre96(this, &G4Diffu    1345   sumL96 = integral.Legendre96(this, &G4DiffuseElastic::GetIntegrandFunction, 0., alpha2);
1434   sumAG  = integral.AdaptiveGauss(this, &G4Di    1346   sumAG  = integral.AdaptiveGauss(this, &G4DiffuseElastic::GetIntegrandFunction, 
1435                                        0., al    1347                                        0., alpha2,epsilon);
1436   G4cout<<G4endl;                                1348   G4cout<<G4endl;
1437   G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/deg    1349   G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1438             <<sumL10<<"\t"<<sumL96<<"\t"<<sum    1350             <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1439   */                                             1351   */
1440   return;                                        1352   return;
1441 }                                                1353 }
1442                                                  1354 
1443 //                                               1355 //
1444 //                                               1356 //
1445 /////////////////////////////////////////////    1357 /////////////////////////////////////////////////////////////////////////////////
1446                                                  1358