| Geant4 Cross Reference |
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: G4NuclNuclDiffuseElastic.cc,v 1.5 2010-11-09 09:04:29 grichine Exp $
>> 27 // GEANT4 tag $Name: geant4-09-04-patch-01 $
26 // 28 //
27 // 29 //
28 // Physics model class G4NuclNuclDiffuseElasti 30 // Physics model class G4NuclNuclDiffuseElastic
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 37
36 #include "G4NuclNuclDiffuseElastic.hh" 38 #include "G4NuclNuclDiffuseElastic.hh"
37 #include "G4ParticleTable.hh" 39 #include "G4ParticleTable.hh"
38 #include "G4ParticleDefinition.hh" 40 #include "G4ParticleDefinition.hh"
39 #include "G4IonTable.hh" 41 #include "G4IonTable.hh"
40 #include "G4NucleiProperties.hh" <<
41 42
42 #include "Randomize.hh" 43 #include "Randomize.hh"
43 #include "G4Integrator.hh" 44 #include "G4Integrator.hh"
44 #include "globals.hh" 45 #include "globals.hh"
45 #include "G4PhysicalConstants.hh" <<
46 #include "G4SystemOfUnits.hh" <<
47 46
48 #include "G4Proton.hh" 47 #include "G4Proton.hh"
49 #include "G4Neutron.hh" 48 #include "G4Neutron.hh"
50 #include "G4Deuteron.hh" 49 #include "G4Deuteron.hh"
51 #include "G4Alpha.hh" 50 #include "G4Alpha.hh"
52 #include "G4PionPlus.hh" 51 #include "G4PionPlus.hh"
53 #include "G4PionMinus.hh" 52 #include "G4PionMinus.hh"
54 53
55 #include "G4Element.hh" 54 #include "G4Element.hh"
56 #include "G4ElementTable.hh" 55 #include "G4ElementTable.hh"
57 #include "G4NistManager.hh" <<
58 #include "G4PhysicsTable.hh" 56 #include "G4PhysicsTable.hh"
59 #include "G4PhysicsLogVector.hh" 57 #include "G4PhysicsLogVector.hh"
60 #include "G4PhysicsFreeVector.hh" 58 #include "G4PhysicsFreeVector.hh"
61 #include "G4HadronicParameters.hh" <<
62 59
63 ////////////////////////////////////////////// 60 /////////////////////////////////////////////////////////////////////////
64 // 61 //
65 // Test Constructor. Just to check xsc 62 // Test Constructor. Just to check xsc
66 63
67 64
68 G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseEla 65 G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseElastic()
69 : G4HadronElastic("NNDiffuseElastic"), fPart << 66 : G4HadronicInteraction(), fParticle(0)
70 { 67 {
71 SetMinEnergy( 50*MeV ); 68 SetMinEnergy( 50*MeV );
72 SetMaxEnergy( G4HadronicParameters::Instance << 69 SetMaxEnergy( 1.*TeV );
73 verboseLevel = 0; 70 verboseLevel = 0;
74 lowEnergyRecoilLimit = 100.*keV; 71 lowEnergyRecoilLimit = 100.*keV;
75 lowEnergyLimitQ = 0.0*GeV; 72 lowEnergyLimitQ = 0.0*GeV;
76 lowEnergyLimitHE = 0.0*GeV; 73 lowEnergyLimitHE = 0.0*GeV;
77 lowestEnergyLimit= 0.0*keV; 74 lowestEnergyLimit= 0.0*keV;
78 plabLowLimit = 20.0*MeV; 75 plabLowLimit = 20.0*MeV;
79 76
80 theProton = G4Proton::Proton(); 77 theProton = G4Proton::Proton();
81 theNeutron = G4Neutron::Neutron(); 78 theNeutron = G4Neutron::Neutron();
82 theDeuteron = G4Deuteron::Deuteron(); 79 theDeuteron = G4Deuteron::Deuteron();
83 theAlpha = G4Alpha::Alpha(); 80 theAlpha = G4Alpha::Alpha();
84 thePionPlus = G4PionPlus::PionPlus(); 81 thePionPlus = G4PionPlus::PionPlus();
85 thePionMinus= G4PionMinus::PionMinus(); 82 thePionMinus= G4PionMinus::PionMinus();
86 83
87 fEnergyBin = 300; // Increased from the ori << 84 fEnergyBin = 200;
88 fAngleBin = 200; 85 fAngleBin = 200;
89 86
90 fEnergyVector = new G4PhysicsLogVector( the 87 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
91 fAngleTable = 0; 88 fAngleTable = 0;
92 89
93 fParticle = 0; 90 fParticle = 0;
94 fWaveVector = 0.; 91 fWaveVector = 0.;
95 fAtomicWeight = 0.; 92 fAtomicWeight = 0.;
96 fAtomicNumber = 0.; 93 fAtomicNumber = 0.;
97 fNuclearRadius = 0.; 94 fNuclearRadius = 0.;
98 fBeta = 0.; 95 fBeta = 0.;
99 fZommerfeld = 0.; 96 fZommerfeld = 0.;
100 fAm = 0.; 97 fAm = 0.;
101 fAddCoulomb = false; 98 fAddCoulomb = false;
102 // Ranges of angle table relative to current 99 // Ranges of angle table relative to current Rutherford (Coulomb grazing) angle
103 100
104 // Empirical parameters <<
105 <<
106 fCofAlphaMax = 1.5; 101 fCofAlphaMax = 1.5;
107 fCofAlphaCoulomb = 0.5; 102 fCofAlphaCoulomb = 0.5;
108 103
109 fProfileDelta = 1.; 104 fProfileDelta = 1.;
110 fProfileAlpha = 0.5; << 105 fProfileAlpha = 0.5;
>> 106
>> 107 }
>> 108
>> 109 //////////////////////////////////////////////////////////////////////////
>> 110 //
>> 111 // Constructor with initialisation
>> 112
>> 113 G4NuclNuclDiffuseElastic::G4NuclNuclDiffuseElastic(const G4ParticleDefinition* aParticle)
>> 114 : G4HadronicInteraction(), fParticle(aParticle)
>> 115 {
>> 116 SetMinEnergy( 50.*MeV); // 0.01*GeV );
>> 117 SetMaxEnergy( 1.*TeV); // 1.*TeV );
>> 118 verboseLevel = 0;
>> 119 lowEnergyRecoilLimit = 100.*keV;
>> 120 lowEnergyLimitQ = 0.0*GeV;
>> 121 lowEnergyLimitHE = 0.0*GeV;
>> 122 lowestEnergyLimit= 0.0*keV;
>> 123 plabLowLimit = 20.0*MeV;
>> 124
>> 125 theProton = G4Proton::Proton();
>> 126 theNeutron = G4Neutron::Neutron();
>> 127 theDeuteron = G4Deuteron::Deuteron();
>> 128 theAlpha = G4Alpha::Alpha();
>> 129 thePionPlus = G4PionPlus::PionPlus();
>> 130 thePionMinus= G4PionMinus::PionMinus();
>> 131
>> 132 fEnergyBin = 200; //200; // 200; // 100;
>> 133 fAngleBin = 400; // 400; // 200; // 100;
>> 134
>> 135 // fEnergyVector = 0;
>> 136 fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
>> 137 fAngleTable = 0;
>> 138
>> 139 fParticle = aParticle;
>> 140 fWaveVector = 0.;
>> 141 fAtomicWeight = 0.;
>> 142 fAtomicNumber = 0.;
>> 143 fNuclearRadius = 0.;
>> 144 fBeta = 0.;
>> 145 fZommerfeld = 0.;
>> 146 fAm = 0.;
>> 147 fAddCoulomb = false;
>> 148
>> 149 // Ranges of angle table relative to current Rutherford (Coulomb grazing) angle
111 150
112 fCofLambda = 1.0; << 151 fCofAlphaMax = 1.5;
113 fCofDelta = 0.04; << 152 fCofAlphaCoulomb = 0.5;
114 fCofAlpha = 0.095; <<
115 <<
116 fNuclearRadius1 = fNuclearRadius2 = fNuclear <<
117 = fRutherfordRatio = fCoulombPhase0 = fHal <<
118 = fRutherfordTheta = fProfileLambda = fCof <<
119 = fSumSigma = fEtaRatio = fReZ = 0.0; <<
120 fMaxL = 0; <<
121 153
122 fNuclearRadiusCof = 1.0; << 154 fProfileDelta = 1.;
123 fCoulombMuC = 0.0; << 155 fProfileAlpha = 0.5;
>> 156
>> 157 // Initialise();
124 } 158 }
125 159
126 ////////////////////////////////////////////// 160 //////////////////////////////////////////////////////////////////////////////
127 // 161 //
128 // Destructor 162 // Destructor
129 163
130 G4NuclNuclDiffuseElastic::~G4NuclNuclDiffuseEl 164 G4NuclNuclDiffuseElastic::~G4NuclNuclDiffuseElastic()
131 { 165 {
132 if ( fEnergyVector ) { << 166 if(fEnergyVector) delete fEnergyVector;
133 delete fEnergyVector; <<
134 fEnergyVector = 0; <<
135 } <<
136 167
137 for ( std::vector<G4PhysicsTable*>::iterator << 168 if( fAngleTable )
138 it != fAngleBank.end(); ++it ) { << 169 {
139 if ( (*it) ) (*it)->clearAndDestroy(); << 170 fAngleTable->clearAndDestroy();
140 delete *it; << 171 delete fAngleTable ;
141 *it = 0; <<
142 } 172 }
143 fAngleTable = 0; <<
144 } 173 }
145 174
146 ////////////////////////////////////////////// 175 //////////////////////////////////////////////////////////////////////////////
147 // 176 //
148 // Initialisation for given particle using ele 177 // Initialisation for given particle using element table of application
149 178
150 void G4NuclNuclDiffuseElastic::Initialise() 179 void G4NuclNuclDiffuseElastic::Initialise()
151 { 180 {
152 181
153 // fEnergyVector = new G4PhysicsLogVector( t 182 // fEnergyVector = new G4PhysicsLogVector( theMinEnergy, theMaxEnergy, fEnergyBin );
154 183
155 const G4ElementTable* theElementTable = G4El 184 const G4ElementTable* theElementTable = G4Element::GetElementTable();
156 std::size_t jEl, numOfEl = G4Element::GetNum << 185 size_t jEl, numOfEl = G4Element::GetNumberOfElements();
157 186
158 // projectile radius 187 // projectile radius
159 188
160 G4double A1 = G4double( fParticle->GetBaryon 189 G4double A1 = G4double( fParticle->GetBaryonNumber() );
161 G4double R1 = CalculateNuclearRad(A1); 190 G4double R1 = CalculateNuclearRad(A1);
162 191
163 for(jEl = 0 ; jEl < numOfEl; ++jEl) // appli 192 for(jEl = 0 ; jEl < numOfEl; ++jEl) // application element loop
164 { 193 {
165 fAtomicNumber = (*theElementTable)[jEl]->G 194 fAtomicNumber = (*theElementTable)[jEl]->GetZ(); // atomic number
166 fAtomicWeight = G4NistManager::Instance()- << 195 fAtomicWeight = (*theElementTable)[jEl]->GetN(); // number of nucleons
167 196
168 fNuclearRadius = CalculateNuclearRad(fAtom 197 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
169 fNuclearRadius += R1; 198 fNuclearRadius += R1;
170 199
171 if(verboseLevel > 0) 200 if(verboseLevel > 0)
172 { 201 {
173 G4cout<<"G4NuclNuclDiffuseElastic::Initi 202 G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element: "
174 <<(*theElementTable)[jEl]->GetName()<<G4 203 <<(*theElementTable)[jEl]->GetName()<<G4endl;
175 } 204 }
176 fElementNumberVector.push_back(fAtomicNumb 205 fElementNumberVector.push_back(fAtomicNumber);
177 fElementNameVector.push_back((*theElementT 206 fElementNameVector.push_back((*theElementTable)[jEl]->GetName());
178 207
179 BuildAngleTable(); 208 BuildAngleTable();
180 fAngleBank.push_back(fAngleTable); 209 fAngleBank.push_back(fAngleTable);
181 } 210 }
>> 211 return;
>> 212 }
>> 213
>> 214 ////////////////////////////////////////////////////////////////////////////////
>> 215 //
>> 216 // Model analog of DoIt function
>> 217
>> 218 G4HadFinalState*
>> 219 G4NuclNuclDiffuseElastic::ApplyYourself( const G4HadProjectile& aTrack,
>> 220 G4Nucleus& targetNucleus )
>> 221 {
>> 222 theParticleChange.Clear();
>> 223
>> 224 const G4HadProjectile* aParticle = &aTrack;
>> 225
>> 226 G4double ekin = aParticle->GetKineticEnergy();
>> 227
>> 228 if(ekin <= lowestEnergyLimit)
>> 229 {
>> 230 theParticleChange.SetEnergyChange(ekin);
>> 231 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
>> 232 return &theParticleChange;
>> 233 }
>> 234
>> 235 G4double aTarget = targetNucleus.GetN();
>> 236 G4double zTarget = targetNucleus.GetZ();
>> 237
>> 238 G4double plab = aParticle->GetTotalMomentum();
>> 239
>> 240 if (verboseLevel >1)
>> 241 {
>> 242 G4cout << "G4NuclNuclDiffuseElastic::DoIt: Incident particle plab="
>> 243 << plab/GeV << " GeV/c "
>> 244 << " ekin(MeV) = " << ekin/MeV << " "
>> 245 << aParticle->GetDefinition()->GetParticleName() << G4endl;
>> 246 }
>> 247 // Scattered particle referred to axis of incident particle
>> 248
>> 249 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
>> 250 G4double m1 = theParticle->GetPDGMass();
>> 251
>> 252 G4int Z = static_cast<G4int>(zTarget+0.5);
>> 253 G4int A = static_cast<G4int>(aTarget+0.5);
>> 254 G4int N = A - Z;
>> 255
>> 256 G4int projPDG = theParticle->GetPDGEncoding();
>> 257
>> 258 if (verboseLevel>1)
>> 259 {
>> 260 G4cout << "G4NuclNuclDiffuseElastic for " << theParticle->GetParticleName()
>> 261 << " PDGcode= " << projPDG << " on nucleus Z= " << Z
>> 262 << " A= " << A << " N= " << N
>> 263 << G4endl;
>> 264 }
>> 265 G4ParticleDefinition * theDef = 0;
>> 266
>> 267 if(Z == 1 && A == 1) theDef = theProton;
>> 268 else if (Z == 1 && A == 2) theDef = theDeuteron;
>> 269 else if (Z == 1 && A == 3) theDef = G4Triton::Triton();
>> 270 else if (Z == 2 && A == 3) theDef = G4He3::He3();
>> 271 else if (Z == 2 && A == 4) theDef = theAlpha;
>> 272 else theDef = G4ParticleTable::GetParticleTable()->FindIon(Z,A,0,Z);
>> 273
>> 274 G4double m2 = theDef->GetPDGMass();
>> 275 G4LorentzVector lv1 = aParticle->Get4Momentum();
>> 276 G4LorentzVector lv(0.0,0.0,0.0,m2);
>> 277 lv += lv1;
>> 278
>> 279 G4ThreeVector bst = lv.boostVector();
>> 280 lv1.boost(-bst);
>> 281
>> 282 G4ThreeVector p1 = lv1.vect();
>> 283 G4double ptot = p1.mag();
>> 284 G4double tmax = 4.0*ptot*ptot;
>> 285 G4double t = 0.0;
>> 286
>> 287
>> 288 //
>> 289 // Sample t
>> 290 //
>> 291
>> 292 // t = SampleT( theParticle, ptot, A);
>> 293
>> 294 t = SampleTableT( theParticle, ptot, Z, A); // use initialised table
>> 295
>> 296 // NaN finder
>> 297 if(!(t < 0.0 || t >= 0.0))
>> 298 {
>> 299 if (verboseLevel > 0)
>> 300 {
>> 301 G4cout << "G4NuclNuclDiffuseElastic:WARNING: Z= " << Z << " N= "
>> 302 << N << " pdg= " << projPDG
>> 303 << " mom(GeV)= " << plab/GeV
>> 304 << " S-wave will be sampled"
>> 305 << G4endl;
>> 306 }
>> 307 t = G4UniformRand()*tmax;
>> 308 }
>> 309 if(verboseLevel>1)
>> 310 {
>> 311 G4cout <<" t= " << t << " tmax= " << tmax
>> 312 << " ptot= " << ptot << G4endl;
>> 313 }
>> 314 // Sampling of angles in CM system
>> 315
>> 316 G4double phi = G4UniformRand()*twopi;
>> 317 G4double cost = 1. - 2.0*t/tmax;
>> 318 G4double sint;
>> 319
>> 320 if( cost >= 1.0 )
>> 321 {
>> 322 cost = 1.0;
>> 323 sint = 0.0;
>> 324 }
>> 325 else if( cost <= -1.0)
>> 326 {
>> 327 cost = -1.0;
>> 328 sint = 0.0;
>> 329 }
>> 330 else
>> 331 {
>> 332 sint = std::sqrt((1.0-cost)*(1.0+cost));
>> 333 }
>> 334 if (verboseLevel>1)
>> 335 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
>> 336
>> 337 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
>> 338 v1 *= ptot;
>> 339 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
>> 340
>> 341 nlv1.boost(bst);
>> 342
>> 343 G4double eFinal = nlv1.e() - m1;
>> 344
>> 345 if (verboseLevel > 1)
>> 346 {
>> 347 G4cout << "Scattered: "
>> 348 << nlv1<<" m= " << m1 << " ekin(MeV)= " << eFinal
>> 349 << " Proj: 4-mom " << lv1
>> 350 <<G4endl;
>> 351 }
>> 352 if(eFinal < 0.0)
>> 353 {
>> 354 G4cout << "G4NuclNuclDiffuseElastic WARNING ekin= " << eFinal
>> 355 << " after scattering of "
>> 356 << aParticle->GetDefinition()->GetParticleName()
>> 357 << " p(GeV/c)= " << plab
>> 358 << " on " << theDef->GetParticleName()
>> 359 << G4endl;
>> 360 eFinal = 0.0;
>> 361 nlv1.setE(m1);
>> 362 }
>> 363
>> 364 theParticleChange.SetMomentumChange(nlv1.vect().unit());
>> 365 theParticleChange.SetEnergyChange(eFinal);
>> 366
>> 367 G4LorentzVector nlv0 = lv - nlv1;
>> 368 G4double erec = nlv0.e() - m2;
>> 369
>> 370 if (verboseLevel > 1)
>> 371 {
>> 372 G4cout << "Recoil: "
>> 373 << nlv0<<" m= " << m2 << " ekin(MeV)= " << erec
>> 374 <<G4endl;
>> 375 }
>> 376 if(erec > lowEnergyRecoilLimit)
>> 377 {
>> 378 G4DynamicParticle * aSec = new G4DynamicParticle(theDef, nlv0);
>> 379 theParticleChange.AddSecondary(aSec);
>> 380 } else {
>> 381 if(erec < 0.0) erec = 0.0;
>> 382 theParticleChange.SetLocalEnergyDeposit(erec);
>> 383 }
>> 384
>> 385 return &theParticleChange;
182 } 386 }
183 387
184 388
185 ////////////////////////////////////////////// 389 ////////////////////////////////////////////////////////////////////////////
186 // 390 //
187 // return differential elastic cross section d 391 // return differential elastic cross section d(sigma)/d(omega)
188 392
189 G4double 393 G4double
190 G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc 394 G4NuclNuclDiffuseElastic::GetDiffuseElasticXsc( const G4ParticleDefinition* particle,
191 G4doub 395 G4double theta,
192 G4double momentum, 396 G4double momentum,
193 G4doub 397 G4double A )
194 { 398 {
195 fParticle = particle; 399 fParticle = particle;
196 fWaveVector = momentum/hbarc; 400 fWaveVector = momentum/hbarc;
197 fAtomicWeight = A; 401 fAtomicWeight = A;
198 fAddCoulomb = false; 402 fAddCoulomb = false;
199 fNuclearRadius = CalculateNuclearRad(A); 403 fNuclearRadius = CalculateNuclearRad(A);
200 404
201 G4double sigma = fNuclearRadius*fNuclearRadi 405 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticProb(theta);
202 406
203 return sigma; 407 return sigma;
204 } 408 }
205 409
206 ////////////////////////////////////////////// 410 ////////////////////////////////////////////////////////////////////////////
207 // 411 //
208 // return invariant differential elastic cross 412 // return invariant differential elastic cross section d(sigma)/d(tMand)
209 413
210 G4double 414 G4double
211 G4NuclNuclDiffuseElastic::GetInvElasticXsc( co 415 G4NuclNuclDiffuseElastic::GetInvElasticXsc( const G4ParticleDefinition* particle,
212 G4doub 416 G4double tMand,
213 G4double plab, 417 G4double plab,
214 G4doub 418 G4double A, G4double Z )
215 { 419 {
216 G4double m1 = particle->GetPDGMass(); 420 G4double m1 = particle->GetPDGMass();
217 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 421 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
218 422
219 G4int iZ = static_cast<G4int>(Z+0.5); 423 G4int iZ = static_cast<G4int>(Z+0.5);
220 G4int iA = static_cast<G4int>(A+0.5); 424 G4int iA = static_cast<G4int>(A+0.5);
221 G4ParticleDefinition * theDef = 0; 425 G4ParticleDefinition * theDef = 0;
222 426
223 if (iZ == 1 && iA == 1) theDef = thePro 427 if (iZ == 1 && iA == 1) theDef = theProton;
224 else if (iZ == 1 && iA == 2) theDef = theDeu 428 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
225 else if (iZ == 1 && iA == 3) theDef = G4Trit 429 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
226 else if (iZ == 2 && iA == 3) theDef = G4He3: 430 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
227 else if (iZ == 2 && iA == 4) theDef = theAlp 431 else if (iZ == 2 && iA == 4) theDef = theAlpha;
228 else theDef = G4ParticleTable::GetParticleTa << 432 else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
229 433
230 G4double tmass = theDef->GetPDGMass(); 434 G4double tmass = theDef->GetPDGMass();
231 435
232 G4LorentzVector lv(0.0,0.0,0.0,tmass); 436 G4LorentzVector lv(0.0,0.0,0.0,tmass);
233 lv += lv1; 437 lv += lv1;
234 438
235 G4ThreeVector bst = lv.boostVector(); 439 G4ThreeVector bst = lv.boostVector();
236 lv1.boost(-bst); 440 lv1.boost(-bst);
237 441
238 G4ThreeVector p1 = lv1.vect(); 442 G4ThreeVector p1 = lv1.vect();
239 G4double ptot = p1.mag(); 443 G4double ptot = p1.mag();
240 G4double ptot2 = ptot*ptot; 444 G4double ptot2 = ptot*ptot;
241 G4double cost = 1 - 0.5*std::fabs(tMand)/pto 445 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
242 446
243 if( cost >= 1.0 ) cost = 1.0; 447 if( cost >= 1.0 ) cost = 1.0;
244 else if( cost <= -1.0) cost = -1.0; 448 else if( cost <= -1.0) cost = -1.0;
245 449
246 G4double thetaCMS = std::acos(cost); 450 G4double thetaCMS = std::acos(cost);
247 451
248 G4double sigma = GetDiffuseElasticXsc( parti 452 G4double sigma = GetDiffuseElasticXsc( particle, thetaCMS, ptot, A);
249 453
250 sigma *= pi/ptot2; 454 sigma *= pi/ptot2;
251 455
252 return sigma; 456 return sigma;
253 } 457 }
254 458
255 ////////////////////////////////////////////// 459 ////////////////////////////////////////////////////////////////////////////
256 // 460 //
257 // return differential elastic cross section d 461 // return differential elastic cross section d(sigma)/d(omega) with Coulomb
258 // correction 462 // correction
259 463
260 G4double 464 G4double
261 G4NuclNuclDiffuseElastic::GetDiffuseElasticSum 465 G4NuclNuclDiffuseElastic::GetDiffuseElasticSumXsc( const G4ParticleDefinition* particle,
262 G4doub 466 G4double theta,
263 G4double momentum, 467 G4double momentum,
264 G4doub 468 G4double A, G4double Z )
265 { 469 {
266 fParticle = particle; 470 fParticle = particle;
267 fWaveVector = momentum/hbarc; 471 fWaveVector = momentum/hbarc;
268 fAtomicWeight = A; 472 fAtomicWeight = A;
269 fAtomicNumber = Z; 473 fAtomicNumber = Z;
270 fNuclearRadius = CalculateNuclearRad(A); 474 fNuclearRadius = CalculateNuclearRad(A);
271 fAddCoulomb = false; 475 fAddCoulomb = false;
272 476
273 G4double z = particle->GetPDGCharge(); 477 G4double z = particle->GetPDGCharge();
274 478
275 G4double kRt = fWaveVector*fNuclearRadius* 479 G4double kRt = fWaveVector*fNuclearRadius*theta;
276 G4double kRtC = 1.9; 480 G4double kRtC = 1.9;
277 481
278 if( z && (kRt > kRtC) ) 482 if( z && (kRt > kRtC) )
279 { 483 {
280 fAddCoulomb = true; 484 fAddCoulomb = true;
281 fBeta = CalculateParticleBeta( parti 485 fBeta = CalculateParticleBeta( particle, momentum);
282 fZommerfeld = CalculateZommerfeld( fBeta, 486 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
283 fAm = CalculateAm( momentum, fZomm 487 fAm = CalculateAm( momentum, fZommerfeld, fAtomicNumber);
284 } 488 }
285 G4double sigma = fNuclearRadius*fNuclearRadi 489 G4double sigma = fNuclearRadius*fNuclearRadius*GetDiffElasticSumProb(theta);
286 490
287 return sigma; 491 return sigma;
288 } 492 }
289 493
290 ////////////////////////////////////////////// 494 ////////////////////////////////////////////////////////////////////////////
291 // 495 //
292 // return invariant differential elastic cross 496 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
293 // correction 497 // correction
294 498
295 G4double 499 G4double
296 G4NuclNuclDiffuseElastic::GetInvElasticSumXsc( 500 G4NuclNuclDiffuseElastic::GetInvElasticSumXsc( const G4ParticleDefinition* particle,
297 G4doub 501 G4double tMand,
298 G4double plab, 502 G4double plab,
299 G4doub 503 G4double A, G4double Z )
300 { 504 {
301 G4double m1 = particle->GetPDGMass(); 505 G4double m1 = particle->GetPDGMass();
302 506
303 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 507 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
304 508
305 G4int iZ = static_cast<G4int>(Z+0.5); 509 G4int iZ = static_cast<G4int>(Z+0.5);
306 G4int iA = static_cast<G4int>(A+0.5); 510 G4int iA = static_cast<G4int>(A+0.5);
307 511
308 G4ParticleDefinition* theDef = 0; 512 G4ParticleDefinition* theDef = 0;
309 513
310 if (iZ == 1 && iA == 1) theDef = thePro 514 if (iZ == 1 && iA == 1) theDef = theProton;
311 else if (iZ == 1 && iA == 2) theDef = theDeu 515 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
312 else if (iZ == 1 && iA == 3) theDef = G4Trit 516 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
313 else if (iZ == 2 && iA == 3) theDef = G4He3: 517 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
314 else if (iZ == 2 && iA == 4) theDef = theAlp 518 else if (iZ == 2 && iA == 4) theDef = theAlpha;
315 else theDef = G4ParticleTable::GetParticleTa << 519 else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
316 520
317 G4double tmass = theDef->GetPDGMass(); 521 G4double tmass = theDef->GetPDGMass();
318 522
319 G4LorentzVector lv(0.0,0.0,0.0,tmass); 523 G4LorentzVector lv(0.0,0.0,0.0,tmass);
320 lv += lv1; 524 lv += lv1;
321 525
322 G4ThreeVector bst = lv.boostVector(); 526 G4ThreeVector bst = lv.boostVector();
323 lv1.boost(-bst); 527 lv1.boost(-bst);
324 528
325 G4ThreeVector p1 = lv1.vect(); 529 G4ThreeVector p1 = lv1.vect();
326 G4double ptot = p1.mag(); 530 G4double ptot = p1.mag();
327 G4double ptot2 = ptot*ptot; 531 G4double ptot2 = ptot*ptot;
328 G4double cost = 1 - 0.5*std::fabs(tMand)/ 532 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
329 533
330 if( cost >= 1.0 ) cost = 1.0; 534 if( cost >= 1.0 ) cost = 1.0;
331 else if( cost <= -1.0) cost = -1.0; 535 else if( cost <= -1.0) cost = -1.0;
332 536
333 G4double thetaCMS = std::acos(cost); 537 G4double thetaCMS = std::acos(cost);
334 538
335 G4double sigma = GetDiffuseElasticSumXsc( pa 539 G4double sigma = GetDiffuseElasticSumXsc( particle, thetaCMS, ptot, A, Z );
336 540
337 sigma *= pi/ptot2; 541 sigma *= pi/ptot2;
338 542
339 return sigma; 543 return sigma;
340 } 544 }
341 545
342 ////////////////////////////////////////////// 546 ////////////////////////////////////////////////////////////////////////////
343 // 547 //
344 // return invariant differential elastic cross 548 // return invariant differential elastic cross section d(sigma)/d(tMand) with Coulomb
345 // correction 549 // correction
346 550
347 G4double 551 G4double
348 G4NuclNuclDiffuseElastic::GetInvCoulombElastic 552 G4NuclNuclDiffuseElastic::GetInvCoulombElasticXsc( const G4ParticleDefinition* particle,
349 G4doub 553 G4double tMand,
350 G4double plab, 554 G4double plab,
351 G4doub 555 G4double A, G4double Z )
352 { 556 {
353 G4double m1 = particle->GetPDGMass(); 557 G4double m1 = particle->GetPDGMass();
354 G4LorentzVector lv1(0.,0.,plab,std::sqrt(pla 558 G4LorentzVector lv1(0.,0.,plab,std::sqrt(plab*plab+m1*m1));
355 559
356 G4int iZ = static_cast<G4int>(Z+0.5); 560 G4int iZ = static_cast<G4int>(Z+0.5);
357 G4int iA = static_cast<G4int>(A+0.5); 561 G4int iA = static_cast<G4int>(A+0.5);
358 G4ParticleDefinition * theDef = 0; 562 G4ParticleDefinition * theDef = 0;
359 563
360 if (iZ == 1 && iA == 1) theDef = thePro 564 if (iZ == 1 && iA == 1) theDef = theProton;
361 else if (iZ == 1 && iA == 2) theDef = theDeu 565 else if (iZ == 1 && iA == 2) theDef = theDeuteron;
362 else if (iZ == 1 && iA == 3) theDef = G4Trit 566 else if (iZ == 1 && iA == 3) theDef = G4Triton::Triton();
363 else if (iZ == 2 && iA == 3) theDef = G4He3: 567 else if (iZ == 2 && iA == 3) theDef = G4He3::He3();
364 else if (iZ == 2 && iA == 4) theDef = theAlp 568 else if (iZ == 2 && iA == 4) theDef = theAlpha;
365 else theDef = G4ParticleTable::GetParticleTa << 569 else theDef = G4ParticleTable::GetParticleTable()->FindIon(iZ,iA,0,iZ);
366 570
367 G4double tmass = theDef->GetPDGMass(); 571 G4double tmass = theDef->GetPDGMass();
368 572
369 G4LorentzVector lv(0.0,0.0,0.0,tmass); 573 G4LorentzVector lv(0.0,0.0,0.0,tmass);
370 lv += lv1; 574 lv += lv1;
371 575
372 G4ThreeVector bst = lv.boostVector(); 576 G4ThreeVector bst = lv.boostVector();
373 lv1.boost(-bst); 577 lv1.boost(-bst);
374 578
375 G4ThreeVector p1 = lv1.vect(); 579 G4ThreeVector p1 = lv1.vect();
376 G4double ptot = p1.mag(); 580 G4double ptot = p1.mag();
377 G4double ptot2 = ptot*ptot; 581 G4double ptot2 = ptot*ptot;
378 G4double cost = 1 - 0.5*std::fabs(tMand)/pto 582 G4double cost = 1 - 0.5*std::fabs(tMand)/ptot2;
379 583
380 if( cost >= 1.0 ) cost = 1.0; 584 if( cost >= 1.0 ) cost = 1.0;
381 else if( cost <= -1.0) cost = -1.0; 585 else if( cost <= -1.0) cost = -1.0;
382 586
383 G4double thetaCMS = std::acos(cost); 587 G4double thetaCMS = std::acos(cost);
384 588
385 G4double sigma = GetCoulombElasticXsc( parti 589 G4double sigma = GetCoulombElasticXsc( particle, thetaCMS, ptot, Z );
386 590
387 sigma *= pi/ptot2; 591 sigma *= pi/ptot2;
388 592
389 return sigma; 593 return sigma;
390 } 594 }
391 595
392 ////////////////////////////////////////////// 596 ////////////////////////////////////////////////////////////////////////////
393 // 597 //
394 // return differential elastic probability d(p 598 // return differential elastic probability d(probability)/d(omega)
395 599
396 G4double 600 G4double
397 G4NuclNuclDiffuseElastic::GetDiffElasticProb( 601 G4NuclNuclDiffuseElastic::GetDiffElasticProb( // G4ParticleDefinition* particle,
398 G4doub 602 G4double theta
399 // G4double momentum, 603 // G4double momentum,
400 // G4double A 604 // G4double A
401 ) 605 )
402 { 606 {
403 G4double sigma, bzero, bzero2, bonebyarg, bo 607 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
404 G4double delta, diffuse, gamma; 608 G4double delta, diffuse, gamma;
405 G4double e1, e2, bone, bone2; 609 G4double e1, e2, bone, bone2;
406 610
407 // G4double wavek = momentum/hbarc; // wave 611 // G4double wavek = momentum/hbarc; // wave vector
408 // G4double r0 = 1.08*fermi; 612 // G4double r0 = 1.08*fermi;
409 // G4double rad = r0*G4Pow::GetInstance()- << 613 // G4double rad = r0*std::pow(A, 1./3.);
410 614
411 G4double kr = fWaveVector*fNuclearRadius; 615 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
412 G4double kr2 = kr*kr; 616 G4double kr2 = kr*kr;
413 G4double krt = kr*theta; 617 G4double krt = kr*theta;
414 618
415 bzero = BesselJzero(krt); 619 bzero = BesselJzero(krt);
416 bzero2 = bzero*bzero; 620 bzero2 = bzero*bzero;
417 bone = BesselJone(krt); 621 bone = BesselJone(krt);
418 bone2 = bone*bone; 622 bone2 = bone*bone;
419 bonebyarg = BesselOneByArg(krt); 623 bonebyarg = BesselOneByArg(krt);
420 bonebyarg2 = bonebyarg*bonebyarg; 624 bonebyarg2 = bonebyarg*bonebyarg;
421 625
422 // VI - Coverity complains <<
423 /* <<
424 if (fParticle == theProton) 626 if (fParticle == theProton)
425 { 627 {
426 diffuse = 0.63*fermi; 628 diffuse = 0.63*fermi;
427 gamma = 0.3*fermi; 629 gamma = 0.3*fermi;
428 delta = 0.1*fermi*fermi; 630 delta = 0.1*fermi*fermi;
429 e1 = 0.3*fermi; 631 e1 = 0.3*fermi;
430 e2 = 0.35*fermi; 632 e2 = 0.35*fermi;
431 } 633 }
432 else // as proton, if were not defined 634 else // as proton, if were not defined
433 { 635 {
434 */ <<
435 diffuse = 0.63*fermi; 636 diffuse = 0.63*fermi;
436 gamma = 0.3*fermi; 637 gamma = 0.3*fermi;
437 delta = 0.1*fermi*fermi; 638 delta = 0.1*fermi*fermi;
438 e1 = 0.3*fermi; 639 e1 = 0.3*fermi;
439 e2 = 0.35*fermi; 640 e2 = 0.35*fermi;
440 //} << 641 }
441 G4double lambda = 15.; // 15 ok 642 G4double lambda = 15.; // 15 ok
442 643
443 // G4double kgamma = fWaveVector*gamma; << 644 // G4double kg = fWaveVector*gamma; // wavek*delta;
444 645
445 G4double kgamma = lambda*(1.-G4Exp(-fWave << 646 G4double kg = lambda*(1.-std::exp(-fWaveVector*gamma/lambda)); // wavek*delta;
446 G4double kgamma2 = kgamma*kgamma; << 647 G4double kg2 = kg*kg;
447 648
448 // G4double dk2t = delta*fWaveVector*fWaveV 649 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
449 // G4double dk2t2 = dk2t*dk2t; 650 // G4double dk2t2 = dk2t*dk2t;
450 // G4double pikdt = pi*fWaveVector*diffuse*t 651 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
451 652
452 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa << 653 G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
453 654
454 damp = DampFactor(pikdt); 655 damp = DampFactor(pikdt);
455 damp2 = damp*damp; 656 damp2 = damp*damp;
456 657
457 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 658 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
458 G4double e2dk3t = -2.*e2*delta*fWaveVector* 659 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
459 660
460 661
461 sigma = kgamma2; << 662 sigma = kg2;
462 // sigma += dk2t2; 663 // sigma += dk2t2;
463 sigma *= bzero2; 664 sigma *= bzero2;
464 sigma += mode2k2*bone2 + e2dk3t*bzero*bone; 665 sigma += mode2k2*bone2 + e2dk3t*bzero*bone;
465 sigma += kr2*bonebyarg2; 666 sigma += kr2*bonebyarg2;
466 sigma *= damp2; // *rad*rad; 667 sigma *= damp2; // *rad*rad;
467 668
468 return sigma; 669 return sigma;
469 } 670 }
470 671
471 ////////////////////////////////////////////// 672 ////////////////////////////////////////////////////////////////////////////
472 // 673 //
473 // return differential elastic probability d(p 674 // return differential elastic probability d(probability)/d(omega) with
474 // Coulomb correction 675 // Coulomb correction
475 676
476 G4double 677 G4double
477 G4NuclNuclDiffuseElastic::GetDiffElasticSumPro 678 G4NuclNuclDiffuseElastic::GetDiffElasticSumProb( // G4ParticleDefinition* particle,
478 G4doub 679 G4double theta
479 // G4double momentum, 680 // G4double momentum,
480 // G4double A 681 // G4double A
481 ) 682 )
482 { 683 {
483 G4double sigma, bzero, bzero2, bonebyarg, bo 684 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
484 G4double delta, diffuse, gamma; 685 G4double delta, diffuse, gamma;
485 G4double e1, e2, bone, bone2; 686 G4double e1, e2, bone, bone2;
486 687
487 // G4double wavek = momentum/hbarc; // wave 688 // G4double wavek = momentum/hbarc; // wave vector
488 // G4double r0 = 1.08*fermi; 689 // G4double r0 = 1.08*fermi;
489 // G4double rad = r0*G4Pow::GetInstance()- << 690 // G4double rad = r0*std::pow(A, 1./3.);
490 691
491 G4double kr = fWaveVector*fNuclearRadius; 692 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
492 G4double kr2 = kr*kr; 693 G4double kr2 = kr*kr;
493 G4double krt = kr*theta; 694 G4double krt = kr*theta;
494 695
495 bzero = BesselJzero(krt); 696 bzero = BesselJzero(krt);
496 bzero2 = bzero*bzero; 697 bzero2 = bzero*bzero;
497 bone = BesselJone(krt); 698 bone = BesselJone(krt);
498 bone2 = bone*bone; 699 bone2 = bone*bone;
499 bonebyarg = BesselOneByArg(krt); 700 bonebyarg = BesselOneByArg(krt);
500 bonebyarg2 = bonebyarg*bonebyarg; 701 bonebyarg2 = bonebyarg*bonebyarg;
501 702
502 if (fParticle == theProton) 703 if (fParticle == theProton)
503 { 704 {
504 diffuse = 0.63*fermi; 705 diffuse = 0.63*fermi;
505 // diffuse = 0.6*fermi; 706 // diffuse = 0.6*fermi;
506 gamma = 0.3*fermi; 707 gamma = 0.3*fermi;
507 delta = 0.1*fermi*fermi; 708 delta = 0.1*fermi*fermi;
508 e1 = 0.3*fermi; 709 e1 = 0.3*fermi;
509 e2 = 0.35*fermi; 710 e2 = 0.35*fermi;
510 } 711 }
511 else // as proton, if were not defined 712 else // as proton, if were not defined
512 { 713 {
513 diffuse = 0.63*fermi; 714 diffuse = 0.63*fermi;
514 gamma = 0.3*fermi; 715 gamma = 0.3*fermi;
515 delta = 0.1*fermi*fermi; 716 delta = 0.1*fermi*fermi;
516 e1 = 0.3*fermi; 717 e1 = 0.3*fermi;
517 e2 = 0.35*fermi; 718 e2 = 0.35*fermi;
518 } 719 }
519 G4double lambda = 15.; // 15 ok 720 G4double lambda = 15.; // 15 ok
520 // G4double kgamma = fWaveVector*gamma; << 721 // G4double kg = fWaveVector*gamma; // wavek*delta;
521 G4double kgamma = lambda*(1.-G4Exp(-fWave << 722 G4double kg = lambda*(1.-std::exp(-fWaveVector*gamma/lambda)); // wavek*delta;
522 723
523 // G4cout<<"kgamma = "<<kgamma<<G4endl; << 724 // G4cout<<"kg = "<<kg<<G4endl;
524 725
525 if(fAddCoulomb) // add Coulomb correction 726 if(fAddCoulomb) // add Coulomb correction
526 { 727 {
527 G4double sinHalfTheta = std::sin(0.5*thet 728 G4double sinHalfTheta = std::sin(0.5*theta);
528 G4double sinHalfTheta2 = sinHalfTheta*sinH 729 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
529 730
530 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta << 731 kg += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
531 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe << 732 // kg += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
532 } 733 }
533 734
534 G4double kgamma2 = kgamma*kgamma; << 735 G4double kg2 = kg*kg;
535 736
536 // G4double dk2t = delta*fWaveVector*fWaveV 737 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
537 // G4cout<<"dk2t = "<<dk2t<<G4endl; 738 // G4cout<<"dk2t = "<<dk2t<<G4endl;
538 // G4double dk2t2 = dk2t*dk2t; 739 // G4double dk2t2 = dk2t*dk2t;
539 // G4double pikdt = pi*fWaveVector*diffuse*t 740 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
540 741
541 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa << 742 G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
542 743
543 // G4cout<<"pikdt = "<<pikdt<<G4endl; 744 // G4cout<<"pikdt = "<<pikdt<<G4endl;
544 745
545 damp = DampFactor(pikdt); 746 damp = DampFactor(pikdt);
546 damp2 = damp*damp; 747 damp2 = damp*damp;
547 748
548 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 749 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
549 G4double e2dk3t = -2.*e2*delta*fWaveVector* 750 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
550 751
551 sigma = kgamma2; << 752 sigma = kg2;
552 // sigma += dk2t2; 753 // sigma += dk2t2;
553 sigma *= bzero2; 754 sigma *= bzero2;
554 sigma += mode2k2*bone2; 755 sigma += mode2k2*bone2;
555 sigma += e2dk3t*bzero*bone; 756 sigma += e2dk3t*bzero*bone;
556 757
557 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf 758 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
558 sigma += kr2*bonebyarg2; // correction at J 759 sigma += kr2*bonebyarg2; // correction at J1()/()
559 760
560 sigma *= damp2; // *rad*rad; 761 sigma *= damp2; // *rad*rad;
561 762
562 return sigma; 763 return sigma;
563 } 764 }
564 765
565 766
566 ////////////////////////////////////////////// 767 ////////////////////////////////////////////////////////////////////////////
567 // 768 //
568 // return differential elastic probability d(p 769 // return differential elastic probability d(probability)/d(t) with
569 // Coulomb correction 770 // Coulomb correction
570 771
571 G4double 772 G4double
572 G4NuclNuclDiffuseElastic::GetDiffElasticSumPro 773 G4NuclNuclDiffuseElastic::GetDiffElasticSumProbA( G4double alpha )
573 { 774 {
574 G4double theta; 775 G4double theta;
575 776
576 theta = std::sqrt(alpha); 777 theta = std::sqrt(alpha);
577 778
578 // theta = std::acos( 1 - alpha/2. ); 779 // theta = std::acos( 1 - alpha/2. );
579 780
580 G4double sigma, bzero, bzero2, bonebyarg, bo 781 G4double sigma, bzero, bzero2, bonebyarg, bonebyarg2, damp, damp2;
581 G4double delta, diffuse, gamma; 782 G4double delta, diffuse, gamma;
582 G4double e1, e2, bone, bone2; 783 G4double e1, e2, bone, bone2;
583 784
584 // G4double wavek = momentum/hbarc; // wave 785 // G4double wavek = momentum/hbarc; // wave vector
585 // G4double r0 = 1.08*fermi; 786 // G4double r0 = 1.08*fermi;
586 // G4double rad = r0*G4Pow::GetInstance()- << 787 // G4double rad = r0*std::pow(A, 1./3.);
587 788
588 G4double kr = fWaveVector*fNuclearRadius; 789 G4double kr = fWaveVector*fNuclearRadius; // wavek*rad;
589 G4double kr2 = kr*kr; 790 G4double kr2 = kr*kr;
590 G4double krt = kr*theta; 791 G4double krt = kr*theta;
591 792
592 bzero = BesselJzero(krt); 793 bzero = BesselJzero(krt);
593 bzero2 = bzero*bzero; 794 bzero2 = bzero*bzero;
594 bone = BesselJone(krt); 795 bone = BesselJone(krt);
595 bone2 = bone*bone; 796 bone2 = bone*bone;
596 bonebyarg = BesselOneByArg(krt); 797 bonebyarg = BesselOneByArg(krt);
597 bonebyarg2 = bonebyarg*bonebyarg; 798 bonebyarg2 = bonebyarg*bonebyarg;
598 799
599 if (fParticle == theProton) 800 if (fParticle == theProton)
600 { 801 {
601 diffuse = 0.63*fermi; 802 diffuse = 0.63*fermi;
602 // diffuse = 0.6*fermi; 803 // diffuse = 0.6*fermi;
603 gamma = 0.3*fermi; 804 gamma = 0.3*fermi;
604 delta = 0.1*fermi*fermi; 805 delta = 0.1*fermi*fermi;
605 e1 = 0.3*fermi; 806 e1 = 0.3*fermi;
606 e2 = 0.35*fermi; 807 e2 = 0.35*fermi;
607 } 808 }
608 else // as proton, if were not defined 809 else // as proton, if were not defined
609 { 810 {
610 diffuse = 0.63*fermi; 811 diffuse = 0.63*fermi;
611 gamma = 0.3*fermi; 812 gamma = 0.3*fermi;
612 delta = 0.1*fermi*fermi; 813 delta = 0.1*fermi*fermi;
613 e1 = 0.3*fermi; 814 e1 = 0.3*fermi;
614 e2 = 0.35*fermi; 815 e2 = 0.35*fermi;
615 } 816 }
616 G4double lambda = 15.; // 15 ok 817 G4double lambda = 15.; // 15 ok
617 // G4double kgamma = fWaveVector*gamma; << 818 // G4double kg = fWaveVector*gamma; // wavek*delta;
618 G4double kgamma = lambda*(1.-G4Exp(-fWave << 819 G4double kg = lambda*(1.-std::exp(-fWaveVector*gamma/lambda)); // wavek*delta;
619 820
620 // G4cout<<"kgamma = "<<kgamma<<G4endl; << 821 // G4cout<<"kg = "<<kg<<G4endl;
621 822
622 if(fAddCoulomb) // add Coulomb correction 823 if(fAddCoulomb) // add Coulomb correction
623 { 824 {
624 G4double sinHalfTheta = theta*0.5; // std 825 G4double sinHalfTheta = theta*0.5; // std::sin(0.5*theta);
625 G4double sinHalfTheta2 = sinHalfTheta*sinH 826 G4double sinHalfTheta2 = sinHalfTheta*sinHalfTheta;
626 827
627 kgamma += 0.5*fZommerfeld/kr/(sinHalfTheta << 828 kg += 0.5*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
628 // kgamma += 0.65*fZommerfeld/kr/(sinHalfThe << 829 // kg += 0.65*fZommerfeld/kr/(sinHalfTheta2+fAm); // correction at J0()
629 } 830 }
630 831
631 G4double kgamma2 = kgamma*kgamma; << 832 G4double kg2 = kg*kg;
632 833
633 // G4double dk2t = delta*fWaveVector*fWaveV 834 // G4double dk2t = delta*fWaveVector*fWaveVector*theta; // delta*wavek*wavek*theta;
634 // G4cout<<"dk2t = "<<dk2t<<G4endl; 835 // G4cout<<"dk2t = "<<dk2t<<G4endl;
635 // G4double dk2t2 = dk2t*dk2t; 836 // G4double dk2t2 = dk2t*dk2t;
636 // G4double pikdt = pi*fWaveVector*diffuse*t 837 // G4double pikdt = pi*fWaveVector*diffuse*theta;// pi*wavek*diffuse*theta;
637 838
638 G4double pikdt = lambda*(1.-G4Exp(-pi*fWa << 839 G4double pikdt = lambda*(1.-std::exp(-pi*fWaveVector*diffuse*theta/lambda)); // wavek*delta;
639 840
640 // G4cout<<"pikdt = "<<pikdt<<G4endl; 841 // G4cout<<"pikdt = "<<pikdt<<G4endl;
641 842
642 damp = DampFactor(pikdt); 843 damp = DampFactor(pikdt);
643 damp2 = damp*damp; 844 damp2 = damp*damp;
644 845
645 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector 846 G4double mode2k2 = (e1*e1+e2*e2)*fWaveVector*fWaveVector;
646 G4double e2dk3t = -2.*e2*delta*fWaveVector* 847 G4double e2dk3t = -2.*e2*delta*fWaveVector*fWaveVector*fWaveVector*theta;
647 848
648 sigma = kgamma2; << 849 sigma = kg2;
649 // sigma += dk2t2; 850 // sigma += dk2t2;
650 sigma *= bzero2; 851 sigma *= bzero2;
651 sigma += mode2k2*bone2; 852 sigma += mode2k2*bone2;
652 sigma += e2dk3t*bzero*bone; 853 sigma += e2dk3t*bzero*bone;
653 854
654 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerf 855 // sigma += kr2*(1 + 8.*fZommerfeld*fZommerfeld/kr2)*bonebyarg2; // correction at J1()/()
655 sigma += kr2*bonebyarg2; // correction at J 856 sigma += kr2*bonebyarg2; // correction at J1()/()
656 857
657 sigma *= damp2; // *rad*rad; 858 sigma *= damp2; // *rad*rad;
658 859
659 return sigma; 860 return sigma;
660 } 861 }
661 862
662 863
663 ////////////////////////////////////////////// 864 ////////////////////////////////////////////////////////////////////////////
664 // 865 //
665 // return differential elastic probability 2*p 866 // return differential elastic probability 2*pi*sin(theta)*d(probability)/d(omega)
666 867
667 G4double 868 G4double
668 G4NuclNuclDiffuseElastic::GetIntegrandFunction 869 G4NuclNuclDiffuseElastic::GetIntegrandFunction( G4double alpha )
669 { 870 {
670 G4double result; 871 G4double result;
671 872
672 result = GetDiffElasticSumProbA(alpha); 873 result = GetDiffElasticSumProbA(alpha);
673 874
674 // result *= 2*pi*std::sin(theta); 875 // result *= 2*pi*std::sin(theta);
675 876
676 return result; 877 return result;
677 } 878 }
678 879
679 ////////////////////////////////////////////// 880 ////////////////////////////////////////////////////////////////////////////
680 // 881 //
681 // return integral elastic cross section d(sig 882 // return integral elastic cross section d(sigma)/d(omega) integrated 0 - theta
682 883
683 G4double 884 G4double
684 G4NuclNuclDiffuseElastic::IntegralElasticProb( 885 G4NuclNuclDiffuseElastic::IntegralElasticProb( const G4ParticleDefinition* particle,
685 G4doub 886 G4double theta,
686 G4double momentum, 887 G4double momentum,
687 G4doub 888 G4double A )
688 { 889 {
689 G4double result; 890 G4double result;
690 fParticle = particle; 891 fParticle = particle;
691 fWaveVector = momentum/hbarc; 892 fWaveVector = momentum/hbarc;
692 fAtomicWeight = A; 893 fAtomicWeight = A;
693 894
694 fNuclearRadius = CalculateNuclearRad(A); 895 fNuclearRadius = CalculateNuclearRad(A);
695 896
696 897
697 G4Integrator<G4NuclNuclDiffuseElastic,G4doub 898 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
698 899
699 // result = integral.Legendre10(this,&G4Nucl 900 // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
700 result = integral.Legendre96(this,&G4NuclNuc 901 result = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
701 902
702 return result; 903 return result;
703 } 904 }
704 905
705 ////////////////////////////////////////////// 906 ////////////////////////////////////////////////////////////////////////////
706 // 907 //
707 // Return inv momentum transfer -t > 0 908 // Return inv momentum transfer -t > 0
708 909
709 G4double G4NuclNuclDiffuseElastic::SampleT( co << 910 G4double G4NuclNuclDiffuseElastic::SampleT( const G4ParticleDefinition* aParticle, G4double p, G4double A)
710 G4double p, G4double A) <<
711 { 911 {
712 G4double theta = SampleThetaCMS( aParticle, 912 G4double theta = SampleThetaCMS( aParticle, p, A); // sample theta in cms
713 G4double t = 2*p*p*( 1 - std::cos(theta) 913 G4double t = 2*p*p*( 1 - std::cos(theta) ); // -t !!!
714 return t; 914 return t;
715 } 915 }
716 916
717 ////////////////////////////////////////////// 917 ////////////////////////////////////////////////////////////////////////////
718 // 918 //
719 // Return scattering angle sampled in cms 919 // Return scattering angle sampled in cms
720 920
721 921
722 G4double 922 G4double
723 G4NuclNuclDiffuseElastic::SampleThetaCMS(const 923 G4NuclNuclDiffuseElastic::SampleThetaCMS(const G4ParticleDefinition* particle,
724 G4doubl 924 G4double momentum, G4double A)
725 { 925 {
726 G4int i, iMax = 100; 926 G4int i, iMax = 100;
727 G4double norm, theta1, theta2, thetaMax; << 927 G4double norm, result, theta1, theta2, thetaMax, sum = 0.;
728 G4double result = 0., sum = 0.; <<
729 928
730 fParticle = particle; 929 fParticle = particle;
731 fWaveVector = momentum/hbarc; 930 fWaveVector = momentum/hbarc;
732 fAtomicWeight = A; 931 fAtomicWeight = A;
733 932
734 fNuclearRadius = CalculateNuclearRad(A); 933 fNuclearRadius = CalculateNuclearRad(A);
735 934
736 thetaMax = 10.174/fWaveVector/fNuclearRadius 935 thetaMax = 10.174/fWaveVector/fNuclearRadius;
737 936
738 if (thetaMax > pi) thetaMax = pi; 937 if (thetaMax > pi) thetaMax = pi;
739 938
740 G4Integrator<G4NuclNuclDiffuseElastic,G4doub 939 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
741 940
742 // result = integral.Legendre10(this,&G4Nucl 941 // result = integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., theta );
743 norm = integral.Legendre96(this,&G4NuclNuclD 942 norm = integral.Legendre96(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., thetaMax );
744 943
745 norm *= G4UniformRand(); 944 norm *= G4UniformRand();
746 945
747 for(i = 1; i <= iMax; i++) 946 for(i = 1; i <= iMax; i++)
748 { 947 {
749 theta1 = (i-1)*thetaMax/iMax; 948 theta1 = (i-1)*thetaMax/iMax;
750 theta2 = i*thetaMax/iMax; 949 theta2 = i*thetaMax/iMax;
751 sum += integral.Legendre10(this,&G4NuclN 950 sum += integral.Legendre10(this,&G4NuclNuclDiffuseElastic::GetIntegrandFunction, theta1, theta2);
752 951
753 if ( sum >= norm ) 952 if ( sum >= norm )
754 { 953 {
755 result = 0.5*(theta1 + theta2); 954 result = 0.5*(theta1 + theta2);
756 break; 955 break;
757 } 956 }
758 } 957 }
759 if (i > iMax ) result = 0.5*(theta1 + theta2 958 if (i > iMax ) result = 0.5*(theta1 + theta2);
760 959
761 G4double sigma = pi*thetaMax/iMax; 960 G4double sigma = pi*thetaMax/iMax;
762 961
763 result += G4RandGauss::shoot(0.,sigma); 962 result += G4RandGauss::shoot(0.,sigma);
764 963
765 if(result < 0.) result = 0.; 964 if(result < 0.) result = 0.;
766 if(result > thetaMax) result = thetaMax; 965 if(result > thetaMax) result = thetaMax;
767 966
768 return result; 967 return result;
769 } 968 }
770 969
771 ////////////////////////////////////////////// 970 /////////////////////////////////////////////////////////////////////////////
772 ///////////////////// Table preparation and r 971 ///////////////////// Table preparation and reading ////////////////////////
773 <<
774 ////////////////////////////////////////////// <<
775 // <<
776 // Interface function. Return inv momentum tra <<
777 <<
778 G4double G4NuclNuclDiffuseElastic::SampleInvar <<
779 <<
780 { <<
781 fParticle = aParticle; <<
782 fAtomicNumber = G4double(Z); <<
783 fAtomicWeight = G4double(A); <<
784 <<
785 G4double t(0.), m1 = fParticle->GetPDGMass() <<
786 G4double totElab = std::sqrt(m1*m1+p*p); <<
787 G4double mass2 = G4NucleiProperties::GetNucl <<
788 G4LorentzVector lv1(p,0.0,0.0,totElab); <<
789 G4LorentzVector lv(0.0,0.0,0.0,mass2); <<
790 lv += lv1; <<
791 <<
792 G4ThreeVector bst = lv.boostVector(); <<
793 lv1.boost(-bst); <<
794 <<
795 G4ThreeVector p1 = lv1.vect(); <<
796 G4double momentumCMS = p1.mag(); <<
797 <<
798 // t = SampleTableT( aParticle, momentumCMS <<
799 <<
800 t = SampleCoulombMuCMS( aParticle, momentum <<
801 <<
802 <<
803 return t; <<
804 } <<
805 <<
806 ////////////////////////////////////////////// <<
807 // <<
808 // Return inv momentum transfer -t > 0 as Coul <<
809 <<
810 G4double G4NuclNuclDiffuseElastic::SampleCoulo <<
811 { <<
812 G4double t(0.), rand(0.), mu(0.); <<
813 <<
814 G4double A1 = G4double( aParticle->GetBaryon <<
815 G4double R1 = CalculateNuclearRad(A1); <<
816 <<
817 fNuclearRadius = CalculateNuclearRad(fAtomi <<
818 fNuclearRadius += R1; <<
819 <<
820 InitDynParameters(fParticle, p); <<
821 <<
822 fCoulombMuC = fHalfRutThetaTg2/(1.+fHalfRutT <<
823 <<
824 rand = G4UniformRand(); <<
825 <<
826 // sample (1-cosTheta) in 0 < Theta < ThetaC <<
827 <<
828 mu = fCoulombMuC*rand*fAm; <<
829 mu /= fAm + fCoulombMuC*( 1. - rand ); <<
830 <<
831 t = 4.*p*p*mu; <<
832 <<
833 return t; <<
834 } <<
835 <<
836 <<
837 ////////////////////////////////////////////// 972 ////////////////////////////////////////////////////////////////////////////
838 // 973 //
839 // Return inv momentum transfer -t > 0 from in 974 // Return inv momentum transfer -t > 0 from initialisation table
840 975
841 G4double G4NuclNuclDiffuseElastic::SampleTable 976 G4double G4NuclNuclDiffuseElastic::SampleTableT( const G4ParticleDefinition* aParticle, G4double p,
842 977 G4double Z, G4double A)
843 { 978 {
844 G4double alpha = SampleTableThetaCMS( aParti 979 G4double alpha = SampleTableThetaCMS( aParticle, p, Z, A); // sample theta2 in cms
845 // G4double t = 2*p*p*( 1 - std::cos(std 980 // G4double t = 2*p*p*( 1 - std::cos(std::sqrt(alpha)) ); // -t !!!
846 G4double t = p*p*alpha; // - 981 G4double t = p*p*alpha; // -t !!!
847 return t; 982 return t;
848 } 983 }
849 984
850 ////////////////////////////////////////////// 985 ////////////////////////////////////////////////////////////////////////////
851 // 986 //
852 // Return scattering angle2 sampled in cms acc 987 // Return scattering angle2 sampled in cms according to precalculated table.
853 988
854 989
855 G4double 990 G4double
856 G4NuclNuclDiffuseElastic::SampleTableThetaCMS( 991 G4NuclNuclDiffuseElastic::SampleTableThetaCMS(const G4ParticleDefinition* particle,
857 G4doubl 992 G4double momentum, G4double Z, G4double A)
858 { 993 {
859 std::size_t iElement; << 994 size_t iElement;
860 G4int iMomentum, iAngle; 995 G4int iMomentum, iAngle;
861 G4double randAngle, position, theta1, theta2 996 G4double randAngle, position, theta1, theta2, E1, E2, W1, W2, W;
862 G4double m1 = particle->GetPDGMass(); 997 G4double m1 = particle->GetPDGMass();
863 998
864 for(iElement = 0; iElement < fElementNumberV 999 for(iElement = 0; iElement < fElementNumberVector.size(); iElement++)
865 { 1000 {
866 if( std::fabs(Z - fElementNumberVector[iEl 1001 if( std::fabs(Z - fElementNumberVector[iElement]) < 0.5) break;
867 } 1002 }
868 if ( iElement == fElementNumberVector.size() 1003 if ( iElement == fElementNumberVector.size() )
869 { 1004 {
870 InitialiseOnFly(Z,A); // table preparation 1005 InitialiseOnFly(Z,A); // table preparation, if needed
871 1006
872 // iElement--; 1007 // iElement--;
873 1008
874 // G4cout << "G4NuclNuclDiffuseElastic: El 1009 // G4cout << "G4NuclNuclDiffuseElastic: Element with atomic number " << Z
875 // << " is not found, return zero angle" < 1010 // << " is not found, return zero angle" << G4endl;
876 // return 0.; // no table for this element 1011 // return 0.; // no table for this element
877 } 1012 }
878 // G4cout<<"iElement = "<<iElement<<G4endl; 1013 // G4cout<<"iElement = "<<iElement<<G4endl;
879 1014
880 fAngleTable = fAngleBank[iElement]; 1015 fAngleTable = fAngleBank[iElement];
881 1016
882 G4double kinE = std::sqrt(momentum*momentum 1017 G4double kinE = std::sqrt(momentum*momentum + m1*m1) - m1;
883 1018
884 for( iMomentum = 0; iMomentum < fEnergyBin; 1019 for( iMomentum = 0; iMomentum < fEnergyBin; iMomentum++)
885 { 1020 {
886 // G4cout<<iMomentum<<", kinE = "<<kinE/Me 1021 // G4cout<<iMomentum<<", kinE = "<<kinE/MeV<<", vectorE = "<<fEnergyVector->GetLowEdgeEnergy(iMomentum)/MeV<<G4endl;
887 1022
888 if( kinE < fEnergyVector->GetLowEdgeEnergy 1023 if( kinE < fEnergyVector->GetLowEdgeEnergy(iMomentum) ) break;
889 } 1024 }
890 // G4cout<<"iMomentum = "<<iMomentum<<", fEn 1025 // G4cout<<"iMomentum = "<<iMomentum<<", fEnergyBin -1 = "<<fEnergyBin -1<<G4endl;
891 1026
892 1027
893 if ( iMomentum >= fEnergyBin ) iMomentum = f 1028 if ( iMomentum >= fEnergyBin ) iMomentum = fEnergyBin-1; // kinE is more then theMaxEnergy
894 if ( iMomentum < 0 ) iMomentum = 0 1029 if ( iMomentum < 0 ) iMomentum = 0; // against negative index, kinE < theMinEnergy
895 1030
896 1031
897 if (iMomentum == fEnergyBin -1 || iMomentum 1032 if (iMomentum == fEnergyBin -1 || iMomentum == 0 ) // the table edges
898 { 1033 {
899 position = (*(*fAngleTable)(iMomentum))(fA 1034 position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
900 1035
901 // G4cout<<"position = "<<position<<G4endl 1036 // G4cout<<"position = "<<position<<G4endl;
902 1037
903 for(iAngle = 0; iAngle < fAngleBin-1; iAng 1038 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
904 { 1039 {
905 if( position < (*(*fAngleTable)(iMomentu 1040 if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
906 } 1041 }
907 if (iAngle >= fAngleBin-1) iAngle = fAngle 1042 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
908 1043
909 // G4cout<<"iAngle = "<<iAngle<<G4endl; 1044 // G4cout<<"iAngle = "<<iAngle<<G4endl;
910 1045
911 randAngle = GetScatteringAngle(iMomentum, 1046 randAngle = GetScatteringAngle(iMomentum, iAngle, position);
912 1047
913 // G4cout<<"randAngle = "<<randAngle<<G4en 1048 // G4cout<<"randAngle = "<<randAngle<<G4endl;
914 } 1049 }
915 else // kinE inside between energy table ed 1050 else // kinE inside between energy table edges
916 { 1051 {
917 // position = (*(*fAngleTable)(iMomentum)) 1052 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
918 position = (*(*fAngleTable)(iMomentum))(0) 1053 position = (*(*fAngleTable)(iMomentum))(0)*G4UniformRand();
919 1054
920 // G4cout<<"position = "<<position<<G4endl 1055 // G4cout<<"position = "<<position<<G4endl;
921 1056
922 for(iAngle = 0; iAngle < fAngleBin-1; iAng 1057 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
923 { 1058 {
924 // if( position < (*(*fAngleTable)(iMome 1059 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
925 if( position > (*(*fAngleTable)(iMomentu 1060 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
926 } 1061 }
927 if (iAngle >= fAngleBin-1) iAngle = fAngle 1062 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
928 1063
929 // G4cout<<"iAngle = "<<iAngle<<G4endl; 1064 // G4cout<<"iAngle = "<<iAngle<<G4endl;
930 1065
931 theta2 = GetScatteringAngle(iMomentum, iA 1066 theta2 = GetScatteringAngle(iMomentum, iAngle, position);
932 1067
933 // G4cout<<"theta2 = "<<theta2<<G4endl; 1068 // G4cout<<"theta2 = "<<theta2<<G4endl;
934 1069
935 E2 = fEnergyVector->GetLowEdgeEnergy(iMome 1070 E2 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
936 1071
937 // G4cout<<"E2 = "<<E2<<G4endl; 1072 // G4cout<<"E2 = "<<E2<<G4endl;
938 1073
939 iMomentum--; 1074 iMomentum--;
940 1075
941 // position = (*(*fAngleTable)(iMomentum)) 1076 // position = (*(*fAngleTable)(iMomentum))(fAngleBin-2)*G4UniformRand();
942 1077
943 // G4cout<<"position = "<<position<<G4endl 1078 // G4cout<<"position = "<<position<<G4endl;
944 1079
945 for(iAngle = 0; iAngle < fAngleBin-1; iAng 1080 for(iAngle = 0; iAngle < fAngleBin-1; iAngle++)
946 { 1081 {
947 // if( position < (*(*fAngleTable)(iMome 1082 // if( position < (*(*fAngleTable)(iMomentum))(iAngle) ) break;
948 if( position > (*(*fAngleTable)(iMomentu 1083 if( position > (*(*fAngleTable)(iMomentum))(iAngle) ) break;
949 } 1084 }
950 if (iAngle >= fAngleBin-1) iAngle = fAngle 1085 if (iAngle >= fAngleBin-1) iAngle = fAngleBin-2;
951 1086
952 theta1 = GetScatteringAngle(iMomentum, iA 1087 theta1 = GetScatteringAngle(iMomentum, iAngle, position);
953 1088
954 // G4cout<<"theta1 = "<<theta1<<G4endl; 1089 // G4cout<<"theta1 = "<<theta1<<G4endl;
955 1090
956 E1 = fEnergyVector->GetLowEdgeEnergy(iMome 1091 E1 = fEnergyVector->GetLowEdgeEnergy(iMomentum);
957 1092
958 // G4cout<<"E1 = "<<E1<<G4endl; 1093 // G4cout<<"E1 = "<<E1<<G4endl;
959 1094
960 W = 1.0/(E2 - E1); 1095 W = 1.0/(E2 - E1);
961 W1 = (E2 - kinE)*W; 1096 W1 = (E2 - kinE)*W;
962 W2 = (kinE - E1)*W; 1097 W2 = (kinE - E1)*W;
963 1098
964 randAngle = W1*theta1 + W2*theta2; 1099 randAngle = W1*theta1 + W2*theta2;
965 1100
966 // randAngle = theta2; 1101 // randAngle = theta2;
967 } 1102 }
968 // G4double angle = randAngle; 1103 // G4double angle = randAngle;
969 // if (randAngle > 0.) randAngle /= 2*pi*std 1104 // if (randAngle > 0.) randAngle /= 2*pi*std::sin(angle);
970 // G4cout<<"randAngle = "<<randAngle/degree< 1105 // G4cout<<"randAngle = "<<randAngle/degree<<G4endl;
971 1106
972 return randAngle; 1107 return randAngle;
973 } 1108 }
974 1109
975 ////////////////////////////////////////////// 1110 //////////////////////////////////////////////////////////////////////////////
976 // 1111 //
977 // Initialisation for given particle on fly us 1112 // Initialisation for given particle on fly using new element number
978 1113
979 void G4NuclNuclDiffuseElastic::InitialiseOnFly 1114 void G4NuclNuclDiffuseElastic::InitialiseOnFly(G4double Z, G4double A)
980 { 1115 {
981 fAtomicNumber = Z; // atomic number 1116 fAtomicNumber = Z; // atomic number
982 fAtomicWeight = G4NistManager::Instance()-> << 1117 fAtomicWeight = A; // number of nucleons
983 1118
984 G4double A1 = G4double( fParticle->GetBaryon 1119 G4double A1 = G4double( fParticle->GetBaryonNumber() );
985 G4double R1 = CalculateNuclearRad(A1); 1120 G4double R1 = CalculateNuclearRad(A1);
986 1121
987 fNuclearRadius = CalculateNuclearRad(fAtomic 1122 fNuclearRadius = CalculateNuclearRad(fAtomicWeight) + R1;
988 1123
989 if( verboseLevel > 0 ) 1124 if( verboseLevel > 0 )
990 { 1125 {
991 G4cout<<"G4NuclNuclDiffuseElastic::Initial 1126 G4cout<<"G4NuclNuclDiffuseElastic::Initialise() the element with Z = "
992 <<Z<<"; and A = "<<A<<G4endl; 1127 <<Z<<"; and A = "<<A<<G4endl;
993 } 1128 }
994 fElementNumberVector.push_back(fAtomicNumber 1129 fElementNumberVector.push_back(fAtomicNumber);
995 1130
996 BuildAngleTable(); 1131 BuildAngleTable();
997 1132
998 fAngleBank.push_back(fAngleTable); 1133 fAngleBank.push_back(fAngleTable);
999 1134
1000 return; 1135 return;
1001 } 1136 }
1002 1137
1003 ///////////////////////////////////////////// 1138 ///////////////////////////////////////////////////////////////////////////////
1004 // 1139 //
1005 // Build for given particle and element table 1140 // Build for given particle and element table of momentum, angle probability.
1006 // For the moment in lab system. 1141 // For the moment in lab system.
1007 1142
1008 void G4NuclNuclDiffuseElastic::BuildAngleTabl 1143 void G4NuclNuclDiffuseElastic::BuildAngleTable()
1009 { 1144 {
1010 G4int i, j; 1145 G4int i, j;
1011 G4double partMom, kinE, m1 = fParticle->Get << 1146 G4double partMom, kinE, a = 0., z = fParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1012 G4double alpha1, alpha2, alphaMax, alphaCou 1147 G4double alpha1, alpha2, alphaMax, alphaCoulomb, delta = 0., sum = 0.;
1013 1148
1014 // G4cout<<"particle z = "<<z<<"; particle 1149 // G4cout<<"particle z = "<<z<<"; particle m1 = "<<m1/GeV<<" GeV"<<G4endl;
1015 1150
1016 G4Integrator<G4NuclNuclDiffuseElastic,G4dou 1151 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1017 1152
1018 fAngleTable = new G4PhysicsTable(fEnergyBin 1153 fAngleTable = new G4PhysicsTable(fEnergyBin);
1019 1154
1020 for( i = 0; i < fEnergyBin; i++) 1155 for( i = 0; i < fEnergyBin; i++)
1021 { 1156 {
1022 kinE = fEnergyVector->GetLowEdgeEn 1157 kinE = fEnergyVector->GetLowEdgeEnergy(i);
1023 1158
1024 // G4cout<<G4endl; 1159 // G4cout<<G4endl;
1025 // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G 1160 // G4cout<<"kinE = "<<kinE/MeV<<" MeV"<<G4endl;
1026 1161
1027 partMom = std::sqrt( kinE*(kinE + 2*m 1162 partMom = std::sqrt( kinE*(kinE + 2*m1) );
1028 1163
1029 InitDynParameters(fParticle, partMom); << 1164 fWaveVector = partMom/hbarc;
>> 1165
>> 1166 G4double kR = fWaveVector*fNuclearRadius;
>> 1167
>> 1168 if( z )
>> 1169 {
>> 1170 a = partMom/m1; // beta*gamma for m1
>> 1171 fBeta = a/std::sqrt(1+a*a);
>> 1172 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
>> 1173 fRutherfordRatio = fZommerfeld/fWaveVector;
>> 1174 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
>> 1175 }
>> 1176 // G4cout<<"fZommerfeld = "<<fZommerfeld<<G4endl;
>> 1177
>> 1178 fProfileLambda = kR; // *std::sqrt(1.-2*fZommerfeld/kR);
>> 1179
>> 1180 // G4cout<<"fProfileLambda = "<<fProfileLambda<<G4endl;
>> 1181
>> 1182 fProfileDelta = fCofDelta*fProfileLambda;
>> 1183 fProfileAlpha = fCofAlpha*fProfileLambda;
>> 1184
>> 1185 // CalculateCoulombPhaseZero();
>> 1186
>> 1187 CalculateRutherfordAnglePar();
1030 1188
1031 alphaMax = fRutherfordTheta*fCofAlphaMax; 1189 alphaMax = fRutherfordTheta*fCofAlphaMax;
1032 1190
1033 if(alphaMax > pi) alphaMax = pi; 1191 if(alphaMax > pi) alphaMax = pi;
1034 1192
1035 // VI: Coverity complain <<
1036 //alphaMax = pi2; <<
1037 1193
1038 alphaCoulomb = fRutherfordTheta*fCofAlpha 1194 alphaCoulomb = fRutherfordTheta*fCofAlphaCoulomb;
1039 1195
1040 // G4cout<<"alphaCoulomb = "<<alphaCoulom 1196 // G4cout<<"alphaCoulomb = "<<alphaCoulomb/degree<<"; alphaMax = "<<alphaMax/degree<<G4endl;
1041 1197
1042 G4PhysicsFreeVector* angleVector = new G4 1198 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1043 1199
1044 // G4PhysicsLogVector* angleBins = new G 1200 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1045 1201
1046 G4double delth = (alphaMax-alphaCoulomb)/ 1202 G4double delth = (alphaMax-alphaCoulomb)/fAngleBin;
1047 1203
1048 sum = 0.; 1204 sum = 0.;
1049 1205
1050 // fAddCoulomb = false; 1206 // fAddCoulomb = false;
1051 fAddCoulomb = true; 1207 fAddCoulomb = true;
1052 1208
1053 // for(j = 1; j < fAngleBin; j++) 1209 // for(j = 1; j < fAngleBin; j++)
1054 for(j = fAngleBin-1; j >= 1; j--) 1210 for(j = fAngleBin-1; j >= 1; j--)
1055 { 1211 {
1056 // alpha1 = angleBins->GetLowEdgeEnergy 1212 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1057 // alpha2 = angleBins->GetLowEdgeEnergy 1213 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1058 1214
1059 // alpha1 = alphaMax*(j-1)/fAngleBin; 1215 // alpha1 = alphaMax*(j-1)/fAngleBin;
1060 // alpha2 = alphaMax*( j )/fAngleBin; 1216 // alpha2 = alphaMax*( j )/fAngleBin;
1061 1217
1062 alpha1 = alphaCoulomb + delth*(j-1); 1218 alpha1 = alphaCoulomb + delth*(j-1);
1063 // if(alpha1 < kRlim2) alpha1 = kRlim2; 1219 // if(alpha1 < kRlim2) alpha1 = kRlim2;
1064 alpha2 = alpha1 + delth; 1220 alpha2 = alpha1 + delth;
1065 1221
1066 delta = integral.Legendre10(this, &G4Nu << 1222 delta = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetFresnelDiffuseXsc, alpha1, alpha2);
1067 // delta = integral.Legendre96(this, &G 1223 // delta = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1068 1224
1069 sum += delta; 1225 sum += delta;
1070 1226
1071 angleVector->PutValue( j-1 , alpha1, su 1227 angleVector->PutValue( j-1 , alpha1, sum ); // alpha2
1072 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = " 1228 // G4cout<<"j-1 = "<<j-1<<"; alpha2 = "<<alpha2/degree<<"; sum = "<<sum<<G4endl;
1073 } 1229 }
1074 fAngleTable->insertAt(i,angleVector); 1230 fAngleTable->insertAt(i,angleVector);
1075 1231
1076 // delete[] angleVector; 1232 // delete[] angleVector;
1077 // delete[] angleBins; 1233 // delete[] angleBins;
1078 } 1234 }
1079 return; 1235 return;
1080 } 1236 }
1081 1237
1082 ///////////////////////////////////////////// 1238 /////////////////////////////////////////////////////////////////////////////////
1083 // 1239 //
1084 // 1240 //
1085 1241
1086 G4double 1242 G4double
1087 G4NuclNuclDiffuseElastic:: GetScatteringAngle 1243 G4NuclNuclDiffuseElastic:: GetScatteringAngle( G4int iMomentum, G4int iAngle, G4double position )
1088 { 1244 {
1089 G4double x1, x2, y1, y2, randAngle; 1245 G4double x1, x2, y1, y2, randAngle;
1090 1246
1091 if( iAngle == 0 ) 1247 if( iAngle == 0 )
1092 { 1248 {
1093 randAngle = (*fAngleTable)(iMomentum)->Ge 1249 randAngle = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1094 // iAngle++; 1250 // iAngle++;
1095 } 1251 }
1096 else 1252 else
1097 { 1253 {
1098 if ( iAngle >= G4int((*fAngleTable)(iMome 1254 if ( iAngle >= G4int((*fAngleTable)(iMomentum)->GetVectorLength()) )
1099 { 1255 {
1100 iAngle = G4int((*fAngleTable)(iMomentum << 1256 iAngle = (*fAngleTable)(iMomentum)->GetVectorLength() - 1;
1101 } 1257 }
1102 y1 = (*(*fAngleTable)(iMomentum))(iAngle- 1258 y1 = (*(*fAngleTable)(iMomentum))(iAngle-1);
1103 y2 = (*(*fAngleTable)(iMomentum))(iAngle) 1259 y2 = (*(*fAngleTable)(iMomentum))(iAngle);
1104 1260
1105 x1 = (*fAngleTable)(iMomentum)->GetLowEdg 1261 x1 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle-1);
1106 x2 = (*fAngleTable)(iMomentum)->GetLowEdg 1262 x2 = (*fAngleTable)(iMomentum)->GetLowEdgeEnergy(iAngle);
1107 1263
1108 if ( x1 == x2 ) randAngle = x2; 1264 if ( x1 == x2 ) randAngle = x2;
1109 else 1265 else
1110 { 1266 {
1111 if ( y1 == y2 ) randAngle = x1 + ( x2 - 1267 if ( y1 == y2 ) randAngle = x1 + ( x2 - x1 )*G4UniformRand();
1112 else 1268 else
1113 { 1269 {
1114 randAngle = x1 + ( position - y1 )*( 1270 randAngle = x1 + ( position - y1 )*( x2 - x1 )/( y2 - y1 );
1115 } 1271 }
1116 } 1272 }
1117 } 1273 }
1118 return randAngle; 1274 return randAngle;
1119 } 1275 }
1120 1276
1121 1277
1122 1278
1123 ///////////////////////////////////////////// 1279 ////////////////////////////////////////////////////////////////////////////
1124 // 1280 //
1125 // Return scattering angle sampled in lab sys 1281 // Return scattering angle sampled in lab system (target at rest)
1126 1282
1127 1283
1128 1284
1129 G4double 1285 G4double
1130 G4NuclNuclDiffuseElastic::SampleThetaLab( con 1286 G4NuclNuclDiffuseElastic::SampleThetaLab( const G4HadProjectile* aParticle,
1131 G4dou 1287 G4double tmass, G4double A)
1132 { 1288 {
1133 const G4ParticleDefinition* theParticle = a 1289 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1134 G4double m1 = theParticle->GetPDGMass(); 1290 G4double m1 = theParticle->GetPDGMass();
1135 G4double plab = aParticle->GetTotalMomentum 1291 G4double plab = aParticle->GetTotalMomentum();
1136 G4LorentzVector lv1 = aParticle->Get4Moment 1292 G4LorentzVector lv1 = aParticle->Get4Momentum();
1137 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1293 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1138 lv += lv1; 1294 lv += lv1;
1139 1295
1140 G4ThreeVector bst = lv.boostVector(); 1296 G4ThreeVector bst = lv.boostVector();
1141 lv1.boost(-bst); 1297 lv1.boost(-bst);
1142 1298
1143 G4ThreeVector p1 = lv1.vect(); 1299 G4ThreeVector p1 = lv1.vect();
1144 G4double ptot = p1.mag(); 1300 G4double ptot = p1.mag();
1145 G4double tmax = 4.0*ptot*ptot; 1301 G4double tmax = 4.0*ptot*ptot;
1146 G4double t = 0.0; 1302 G4double t = 0.0;
1147 1303
1148 1304
1149 // 1305 //
1150 // Sample t 1306 // Sample t
1151 // 1307 //
1152 1308
1153 t = SampleT( theParticle, ptot, A); 1309 t = SampleT( theParticle, ptot, A);
1154 1310
1155 // NaN finder 1311 // NaN finder
1156 if(!(t < 0.0 || t >= 0.0)) 1312 if(!(t < 0.0 || t >= 0.0))
1157 { 1313 {
1158 if (verboseLevel > 0) 1314 if (verboseLevel > 0)
1159 { 1315 {
1160 G4cout << "G4NuclNuclDiffuseElastic:WAR 1316 G4cout << "G4NuclNuclDiffuseElastic:WARNING: A = " << A
1161 << " mom(GeV)= " << plab/GeV 1317 << " mom(GeV)= " << plab/GeV
1162 << " S-wave will be sampled" 1318 << " S-wave will be sampled"
1163 << G4endl; 1319 << G4endl;
1164 } 1320 }
1165 t = G4UniformRand()*tmax; 1321 t = G4UniformRand()*tmax;
1166 } 1322 }
1167 if(verboseLevel>1) 1323 if(verboseLevel>1)
1168 { 1324 {
1169 G4cout <<" t= " << t << " tmax= " << tmax 1325 G4cout <<" t= " << t << " tmax= " << tmax
1170 << " ptot= " << ptot << G4endl; 1326 << " ptot= " << ptot << G4endl;
1171 } 1327 }
1172 // Sampling of angles in CM system 1328 // Sampling of angles in CM system
1173 1329
1174 G4double phi = G4UniformRand()*twopi; 1330 G4double phi = G4UniformRand()*twopi;
1175 G4double cost = 1. - 2.0*t/tmax; 1331 G4double cost = 1. - 2.0*t/tmax;
1176 G4double sint; 1332 G4double sint;
1177 1333
1178 if( cost >= 1.0 ) 1334 if( cost >= 1.0 )
1179 { 1335 {
1180 cost = 1.0; 1336 cost = 1.0;
1181 sint = 0.0; 1337 sint = 0.0;
1182 } 1338 }
1183 else if( cost <= -1.0) 1339 else if( cost <= -1.0)
1184 { 1340 {
1185 cost = -1.0; 1341 cost = -1.0;
1186 sint = 0.0; 1342 sint = 0.0;
1187 } 1343 }
1188 else 1344 else
1189 { 1345 {
1190 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1346 sint = std::sqrt((1.0-cost)*(1.0+cost));
1191 } 1347 }
1192 if (verboseLevel>1) 1348 if (verboseLevel>1)
1193 { 1349 {
1194 G4cout << "cos(t)=" << cost << " std::sin 1350 G4cout << "cos(t)=" << cost << " std::sin(t)=" << sint << G4endl;
1195 } 1351 }
1196 G4ThreeVector v1(sint*std::cos(phi),sint*st 1352 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1197 v1 *= ptot; 1353 v1 *= ptot;
1198 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1354 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1199 1355
1200 nlv1.boost(bst); 1356 nlv1.boost(bst);
1201 1357
1202 G4ThreeVector np1 = nlv1.vect(); 1358 G4ThreeVector np1 = nlv1.vect();
1203 1359
1204 // G4double theta = std::acos( np1.z()/np 1360 // G4double theta = std::acos( np1.z()/np1.mag() ); // degree;
1205 1361
1206 G4double theta = np1.theta(); 1362 G4double theta = np1.theta();
1207 1363
1208 return theta; 1364 return theta;
1209 } 1365 }
1210 1366
1211 ///////////////////////////////////////////// 1367 ////////////////////////////////////////////////////////////////////////////
1212 // 1368 //
1213 // Return scattering angle in lab system (tar 1369 // Return scattering angle in lab system (target at rest) knowing theta in CMS
1214 1370
1215 1371
1216 1372
1217 G4double 1373 G4double
1218 G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab( 1374 G4NuclNuclDiffuseElastic::ThetaCMStoThetaLab( const G4DynamicParticle* aParticle,
1219 G4dou 1375 G4double tmass, G4double thetaCMS)
1220 { 1376 {
1221 const G4ParticleDefinition* theParticle = a 1377 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1222 G4double m1 = theParticle->GetPDGMass(); 1378 G4double m1 = theParticle->GetPDGMass();
1223 // G4double plab = aParticle->GetTotalMomen 1379 // G4double plab = aParticle->GetTotalMomentum();
1224 G4LorentzVector lv1 = aParticle->Get4Moment 1380 G4LorentzVector lv1 = aParticle->Get4Momentum();
1225 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1381 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1226 1382
1227 lv += lv1; 1383 lv += lv1;
1228 1384
1229 G4ThreeVector bst = lv.boostVector(); 1385 G4ThreeVector bst = lv.boostVector();
1230 1386
1231 lv1.boost(-bst); 1387 lv1.boost(-bst);
1232 1388
1233 G4ThreeVector p1 = lv1.vect(); 1389 G4ThreeVector p1 = lv1.vect();
1234 G4double ptot = p1.mag(); 1390 G4double ptot = p1.mag();
1235 1391
1236 G4double phi = G4UniformRand()*twopi; 1392 G4double phi = G4UniformRand()*twopi;
1237 G4double cost = std::cos(thetaCMS); 1393 G4double cost = std::cos(thetaCMS);
1238 G4double sint; 1394 G4double sint;
1239 1395
1240 if( cost >= 1.0 ) 1396 if( cost >= 1.0 )
1241 { 1397 {
1242 cost = 1.0; 1398 cost = 1.0;
1243 sint = 0.0; 1399 sint = 0.0;
1244 } 1400 }
1245 else if( cost <= -1.0) 1401 else if( cost <= -1.0)
1246 { 1402 {
1247 cost = -1.0; 1403 cost = -1.0;
1248 sint = 0.0; 1404 sint = 0.0;
1249 } 1405 }
1250 else 1406 else
1251 { 1407 {
1252 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1408 sint = std::sqrt((1.0-cost)*(1.0+cost));
1253 } 1409 }
1254 if (verboseLevel>1) 1410 if (verboseLevel>1)
1255 { 1411 {
1256 G4cout << "cos(tcms)=" << cost << " std:: 1412 G4cout << "cos(tcms)=" << cost << " std::sin(tcms)=" << sint << G4endl;
1257 } 1413 }
1258 G4ThreeVector v1(sint*std::cos(phi),sint*st 1414 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1259 v1 *= ptot; 1415 v1 *= ptot;
1260 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1416 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(ptot*ptot + m1*m1));
1261 1417
1262 nlv1.boost(bst); 1418 nlv1.boost(bst);
1263 1419
1264 G4ThreeVector np1 = nlv1.vect(); 1420 G4ThreeVector np1 = nlv1.vect();
1265 1421
1266 1422
1267 G4double thetaLab = np1.theta(); 1423 G4double thetaLab = np1.theta();
1268 1424
1269 return thetaLab; 1425 return thetaLab;
1270 } 1426 }
1271 1427
1272 ///////////////////////////////////////////// 1428 ////////////////////////////////////////////////////////////////////////////
1273 // 1429 //
1274 // Return scattering angle in CMS system (tar 1430 // Return scattering angle in CMS system (target at rest) knowing theta in Lab
1275 1431
1276 1432
1277 1433
1278 G4double 1434 G4double
1279 G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS( 1435 G4NuclNuclDiffuseElastic::ThetaLabToThetaCMS( const G4DynamicParticle* aParticle,
1280 G4dou 1436 G4double tmass, G4double thetaLab)
1281 { 1437 {
1282 const G4ParticleDefinition* theParticle = a 1438 const G4ParticleDefinition* theParticle = aParticle->GetDefinition();
1283 G4double m1 = theParticle->GetPDGMass(); 1439 G4double m1 = theParticle->GetPDGMass();
1284 G4double plab = aParticle->GetTotalMomentum 1440 G4double plab = aParticle->GetTotalMomentum();
1285 G4LorentzVector lv1 = aParticle->Get4Moment 1441 G4LorentzVector lv1 = aParticle->Get4Momentum();
1286 G4LorentzVector lv(0.0,0.0,0.0,tmass); 1442 G4LorentzVector lv(0.0,0.0,0.0,tmass);
1287 1443
1288 lv += lv1; 1444 lv += lv1;
1289 1445
1290 G4ThreeVector bst = lv.boostVector(); 1446 G4ThreeVector bst = lv.boostVector();
1291 1447
1292 // lv1.boost(-bst); 1448 // lv1.boost(-bst);
1293 1449
1294 // G4ThreeVector p1 = lv1.vect(); 1450 // G4ThreeVector p1 = lv1.vect();
1295 // G4double ptot = p1.mag(); 1451 // G4double ptot = p1.mag();
1296 1452
1297 G4double phi = G4UniformRand()*twopi; 1453 G4double phi = G4UniformRand()*twopi;
1298 G4double cost = std::cos(thetaLab); 1454 G4double cost = std::cos(thetaLab);
1299 G4double sint; 1455 G4double sint;
1300 1456
1301 if( cost >= 1.0 ) 1457 if( cost >= 1.0 )
1302 { 1458 {
1303 cost = 1.0; 1459 cost = 1.0;
1304 sint = 0.0; 1460 sint = 0.0;
1305 } 1461 }
1306 else if( cost <= -1.0) 1462 else if( cost <= -1.0)
1307 { 1463 {
1308 cost = -1.0; 1464 cost = -1.0;
1309 sint = 0.0; 1465 sint = 0.0;
1310 } 1466 }
1311 else 1467 else
1312 { 1468 {
1313 sint = std::sqrt((1.0-cost)*(1.0+cost)); 1469 sint = std::sqrt((1.0-cost)*(1.0+cost));
1314 } 1470 }
1315 if (verboseLevel>1) 1471 if (verboseLevel>1)
1316 { 1472 {
1317 G4cout << "cos(tlab)=" << cost << " std:: 1473 G4cout << "cos(tlab)=" << cost << " std::sin(tlab)=" << sint << G4endl;
1318 } 1474 }
1319 G4ThreeVector v1(sint*std::cos(phi),sint*st 1475 G4ThreeVector v1(sint*std::cos(phi),sint*std::sin(phi),cost);
1320 v1 *= plab; 1476 v1 *= plab;
1321 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),s 1477 G4LorentzVector nlv1(v1.x(),v1.y(),v1.z(),std::sqrt(plab*plab + m1*m1));
1322 1478
1323 nlv1.boost(-bst); 1479 nlv1.boost(-bst);
1324 1480
1325 G4ThreeVector np1 = nlv1.vect(); 1481 G4ThreeVector np1 = nlv1.vect();
1326 1482
1327 1483
1328 G4double thetaCMS = np1.theta(); 1484 G4double thetaCMS = np1.theta();
1329 1485
1330 return thetaCMS; 1486 return thetaCMS;
1331 } 1487 }
1332 1488
1333 ///////////////////////////////////////////// 1489 ///////////////////////////////////////////////////////////////////////////////
1334 // 1490 //
1335 // Test for given particle and element table 1491 // Test for given particle and element table of momentum, angle probability.
1336 // For the moment in lab system. 1492 // For the moment in lab system.
1337 1493
1338 void G4NuclNuclDiffuseElastic::TestAngleTable 1494 void G4NuclNuclDiffuseElastic::TestAngleTable(const G4ParticleDefinition* theParticle, G4double partMom,
1339 G 1495 G4double Z, G4double A)
1340 { 1496 {
1341 fAtomicNumber = Z; // atomic number 1497 fAtomicNumber = Z; // atomic number
1342 fAtomicWeight = A; // number of nucleo 1498 fAtomicWeight = A; // number of nucleons
1343 fNuclearRadius = CalculateNuclearRad(fAtomi 1499 fNuclearRadius = CalculateNuclearRad(fAtomicWeight);
1344 1500
1345 1501
1346 1502
1347 G4cout<<"G4NuclNuclDiffuseElastic::TestAngl 1503 G4cout<<"G4NuclNuclDiffuseElastic::TestAngleTable() init the element with Z = "
1348 <<Z<<"; and A = "<<A<<G4endl; 1504 <<Z<<"; and A = "<<A<<G4endl;
1349 1505
1350 fElementNumberVector.push_back(fAtomicNumbe 1506 fElementNumberVector.push_back(fAtomicNumber);
1351 1507
1352 1508
1353 1509
1354 1510
1355 G4int i=0, j; 1511 G4int i=0, j;
1356 G4double a = 0., z = theParticle->GetPDGCha 1512 G4double a = 0., z = theParticle->GetPDGCharge(), m1 = fParticle->GetPDGMass();
1357 G4double alpha1=0., alpha2=0., alphaMax=0., 1513 G4double alpha1=0., alpha2=0., alphaMax=0., alphaCoulomb=0.;
1358 G4double deltaL10 = 0., deltaL96 = 0., delt 1514 G4double deltaL10 = 0., deltaL96 = 0., deltaAG = 0.;
1359 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0. 1515 G4double sumL10 = 0.,sumL96 = 0.,sumAG = 0.;
1360 G4double epsilon = 0.001; 1516 G4double epsilon = 0.001;
1361 1517
1362 G4Integrator<G4NuclNuclDiffuseElastic,G4dou 1518 G4Integrator<G4NuclNuclDiffuseElastic,G4double(G4NuclNuclDiffuseElastic::*)(G4double)> integral;
1363 1519
1364 fAngleTable = new G4PhysicsTable(fEnergyBin 1520 fAngleTable = new G4PhysicsTable(fEnergyBin);
1365 1521
1366 fWaveVector = partMom/hbarc; 1522 fWaveVector = partMom/hbarc;
1367 1523
1368 G4double kR = fWaveVector*fNuclearRadiu 1524 G4double kR = fWaveVector*fNuclearRadius;
1369 G4double kR2 = kR*kR; 1525 G4double kR2 = kR*kR;
1370 G4double kRmax = 10.6; // 10.6, 18, 10.174 1526 G4double kRmax = 10.6; // 10.6, 18, 10.174; ~ 3 maxima of J1 or 15., 25.
1371 G4double kRcoul = 1.2; // 1.4, 2.5; // on t 1527 G4double kRcoul = 1.2; // 1.4, 2.5; // on the first slope of J1
1372 1528
1373 alphaMax = kRmax*kRmax/kR2; 1529 alphaMax = kRmax*kRmax/kR2;
1374 1530
1375 if (alphaMax > 4.) alphaMax = 4.; // vmg05 1531 if (alphaMax > 4.) alphaMax = 4.; // vmg05-02-09: was pi2
1376 1532
1377 alphaCoulomb = kRcoul*kRcoul/kR2; 1533 alphaCoulomb = kRcoul*kRcoul/kR2;
1378 1534
1379 if( z ) 1535 if( z )
1380 { 1536 {
1381 a = partMom/m1; // beta*gamma 1537 a = partMom/m1; // beta*gamma for m1
1382 fBeta = a/std::sqrt(1+a*a); 1538 fBeta = a/std::sqrt(1+a*a);
1383 fZommerfeld = CalculateZommerfeld( fBet 1539 fZommerfeld = CalculateZommerfeld( fBeta, z, fAtomicNumber);
1384 fAm = CalculateAm( partMom, fZo 1540 fAm = CalculateAm( partMom, fZommerfeld, fAtomicNumber);
1385 } 1541 }
1386 G4PhysicsFreeVector* angleVector = new G4Ph 1542 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fAngleBin-1);
1387 1543
1388 // G4PhysicsLogVector* angleBins = new G4P 1544 // G4PhysicsLogVector* angleBins = new G4PhysicsLogVector( 0.001*alphaMax, alphaMax, fAngleBin );
1389 1545
1390 1546
1391 fAddCoulomb = false; 1547 fAddCoulomb = false;
1392 1548
1393 for(j = 1; j < fAngleBin; j++) 1549 for(j = 1; j < fAngleBin; j++)
1394 { 1550 {
1395 // alpha1 = angleBins->GetLowEdgeEnergy 1551 // alpha1 = angleBins->GetLowEdgeEnergy(j-1);
1396 // alpha2 = angleBins->GetLowEdgeEnergy 1552 // alpha2 = angleBins->GetLowEdgeEnergy(j);
1397 1553
1398 alpha1 = alphaMax*(j-1)/fAngleBin; 1554 alpha1 = alphaMax*(j-1)/fAngleBin;
1399 alpha2 = alphaMax*( j )/fAngleBin; 1555 alpha2 = alphaMax*( j )/fAngleBin;
1400 1556
1401 if( ( alpha2 > alphaCoulomb ) && z ) fAdd 1557 if( ( alpha2 > alphaCoulomb ) && z ) fAddCoulomb = true;
1402 1558
1403 deltaL10 = integral.Legendre10(this, &G4N 1559 deltaL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1404 deltaL96 = integral.Legendre96(this, &G4N 1560 deltaL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, alpha1, alpha2);
1405 deltaAG = integral.AdaptiveGauss(this, & 1561 deltaAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction,
1406 alpha1 1562 alpha1, alpha2,epsilon);
1407 1563
1408 // G4cout<<alpha1<<"\t"<<std::sqrt(alph 1564 // G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1409 // <<deltaL10<<"\t"<<deltaL96<<"\t" 1565 // <<deltaL10<<"\t"<<deltaL96<<"\t"<<deltaAG<<G4endl;
1410 1566
1411 sumL10 += deltaL10; 1567 sumL10 += deltaL10;
1412 sumL96 += deltaL96; 1568 sumL96 += deltaL96;
1413 sumAG += deltaAG; 1569 sumAG += deltaAG;
1414 1570
1415 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/d 1571 G4cout<<alpha1<<"\t"<<std::sqrt(alpha1)/degree<<"\t"
1416 <<sumL10<<"\t"<<sumL96<<"\t"<<sum 1572 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1417 1573
1418 angleVector->PutValue( j-1 , alpha1, sumL 1574 angleVector->PutValue( j-1 , alpha1, sumL10 ); // alpha2
1419 } 1575 }
1420 fAngleTable->insertAt(i,angleVector); 1576 fAngleTable->insertAt(i,angleVector);
1421 fAngleBank.push_back(fAngleTable); 1577 fAngleBank.push_back(fAngleTable);
1422 1578
1423 /* 1579 /*
1424 // Integral over all angle range - Bad accu 1580 // Integral over all angle range - Bad accuracy !!!
1425 1581
1426 sumL10 = integral.Legendre10(this, &G4NuclN 1582 sumL10 = integral.Legendre10(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1427 sumL96 = integral.Legendre96(this, &G4NuclN 1583 sumL96 = integral.Legendre96(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction, 0., alpha2);
1428 sumAG = integral.AdaptiveGauss(this, &G4Nu 1584 sumAG = integral.AdaptiveGauss(this, &G4NuclNuclDiffuseElastic::GetIntegrandFunction,
1429 0., al 1585 0., alpha2,epsilon);
1430 G4cout<<G4endl; 1586 G4cout<<G4endl;
1431 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/deg 1587 G4cout<<alpha2<<"\t"<<std::sqrt(alpha2)/degree<<"\t"
1432 <<sumL10<<"\t"<<sumL96<<"\t"<<sum 1588 <<sumL10<<"\t"<<sumL96<<"\t"<<sumAG<<G4endl;
1433 */ 1589 */
1434 return; 1590 return;
1435 } <<
1436 <<
1437 ///////////////////////////////////////////// <<
1438 // <<
1439 // <<
1440 <<
1441 G4double G4NuclNuclDiffuseElastic::GetLegend <<
1442 { <<
1443 G4double legPol, epsilon = 1.e-6; <<
1444 G4double x = std::cos(theta); <<
1445 <<
1446 if ( n < 0 ) legPol = 0.; <<
1447 else if( n == 0 ) legPol = 1.; <<
1448 else if( n == 1 ) legPol = x; <<
1449 else if( n == 2 ) legPol = (3.*x*x-1.)/2.; <<
1450 else if( n == 3 ) legPol = (5.*x*x*x-3.*x)/ <<
1451 else if( n == 4 ) legPol = (35.*x*x*x*x-30. <<
1452 else if( n == 5 ) legPol = (63.*x*x*x*x*x-7 <<
1453 else if( n == 6 ) legPol = (231.*x*x*x*x*x* <<
1454 else <<
1455 { <<
1456 // legPol = ( (2*n-1)*x*GetLegendrePol(n- <<
1457 <<
1458 legPol = std::sqrt( 2./(n*CLHEP::pi*std:: <<
1459 } <<
1460 return legPol; <<
1461 } <<
1462 <<
1463 ///////////////////////////////////////////// <<
1464 // <<
1465 // <<
1466 <<
1467 G4complex G4NuclNuclDiffuseElastic::GetErfCom <<
1468 { <<
1469 G4int n; <<
1470 G4double n2, cofn, shny, chny, fn, gn; <<
1471 <<
1472 G4double x = z.real(); <<
1473 G4double y = z.imag(); <<
1474 <<
1475 G4double outRe = 0., outIm = 0.; <<
1476 <<
1477 G4double twox = 2.*x; <<
1478 G4double twoxy = twox*y; <<
1479 G4double twox2 = twox*twox; <<
1480 <<
1481 G4double cof1 = G4Exp(-x*x)/CLHEP::pi; <<
1482 <<
1483 G4double cos2xy = std::cos(twoxy); <<
1484 G4double sin2xy = std::sin(twoxy); <<
1485 <<
1486 G4double twoxcos2xy = twox*cos2xy; <<
1487 G4double twoxsin2xy = twox*sin2xy; <<
1488 <<
1489 for( n = 1; n <= nMax; n++) <<
1490 { <<
1491 n2 = n*n; <<
1492 <<
1493 cofn = G4Exp(-0.5*n2)/(n2+twox2); // /(n <<
1494 <<
1495 chny = std::cosh(n*y); <<
1496 shny = std::sinh(n*y); <<
1497 <<
1498 fn = twox - twoxcos2xy*chny + n*sin2xy*s <<
1499 gn = twoxsin2xy*chny + n*cos2xy*s <<
1500 <<
1501 fn *= cofn; <<
1502 gn *= cofn; <<
1503 <<
1504 outRe += fn; <<
1505 outIm += gn; <<
1506 } <<
1507 outRe *= 2*cof1; <<
1508 outIm *= 2*cof1; <<
1509 <<
1510 if(std::abs(x) < 0.0001) <<
1511 { <<
1512 outRe += GetErf(x); <<
1513 outIm += cof1*y; <<
1514 } <<
1515 else <<
1516 { <<
1517 outRe += GetErf(x) + cof1*(1-cos2xy)/twox <<
1518 outIm += cof1*sin2xy/twox; <<
1519 } <<
1520 return G4complex(outRe, outIm); <<
1521 } <<
1522 <<
1523 <<
1524 ///////////////////////////////////////////// <<
1525 // <<
1526 // <<
1527 <<
1528 G4complex G4NuclNuclDiffuseElastic::GetErfInt <<
1529 { <<
1530 G4double outRe, outIm; <<
1531 <<
1532 G4double x = z.real(); <<
1533 G4double y = z.imag(); <<
1534 fReZ = x; <<
1535 <<
1536 G4Integrator<G4NuclNuclDiffuseElastic,G4dou <<
1537 <<
1538 outRe = integral.Legendre96(this,&G4NuclNuc <<
1539 outIm = integral.Legendre96(this,&G4NuclNuc <<
1540 <<
1541 outRe *= 2./std::sqrt(CLHEP::pi); <<
1542 outIm *= 2./std::sqrt(CLHEP::pi); <<
1543 <<
1544 outRe += GetErf(x); <<
1545 <<
1546 return G4complex(outRe, outIm); <<
1547 } <<
1548 <<
1549 ///////////////////////////////////////////// <<
1550 // <<
1551 // <<
1552 <<
1553 <<
1554 G4complex G4NuclNuclDiffuseElastic::GammaLess <<
1555 { <<
1556 G4double sinThetaR = 2.*fHalfRutThetaT <<
1557 G4double cosHalfThetaR2 = 1./(1. + fHalfRut <<
1558 <<
1559 G4double u = std::sqrt(0.5*fPr <<
1560 G4double kappa = u/std::sqrt(CLHEP <<
1561 G4double dTheta = theta - fRutherfo <<
1562 u *= dTheta; <<
1563 G4double u2 = u*u; <<
1564 G4double u2m2p3 = u2*2./3.; <<
1565 <<
1566 G4complex im = G4complex(0.,1.); <<
1567 G4complex order = G4complex(u,u); <<
1568 order /= std::sqrt(2.); <<
1569 <<
1570 G4complex gamma = CLHEP::pi*kappa*G <<
1571 G4complex a0 = 0.5*(1. + 4.*(1.+ <<
1572 G4complex a1 = 0.5*(1. + 2.*(1.+ <<
1573 G4complex out = gamma*(1. - a1*dT <<
1574 <<
1575 return out; <<
1576 } <<
1577 <<
1578 ///////////////////////////////////////////// <<
1579 // <<
1580 // <<
1581 <<
1582 G4complex G4NuclNuclDiffuseElastic::GammaMore <<
1583 { <<
1584 G4double sinThetaR = 2.*fHalfRutThetaT <<
1585 G4double cosHalfThetaR2 = 1./(1. + fHalfRut <<
1586 <<
1587 G4double u = std::sqrt(0.5*fPr <<
1588 G4double kappa = u/std::sqrt(CLHEP <<
1589 G4double dTheta = theta - fRutherfo <<
1590 u *= dTheta; <<
1591 G4double u2 = u*u; <<
1592 G4double u2m2p3 = u2*2./3.; <<
1593 <<
1594 G4complex im = G4complex(0.,1.); <<
1595 G4complex order = G4complex(u,u); <<
1596 order /= std::sqrt(2.); <<
1597 G4complex gamma = CLHEP::pi*kappa*G <<
1598 G4complex a0 = 0.5*(1. + 4.*(1.+ <<
1599 G4complex a1 = 0.5*(1. + 2.*(1.+ <<
1600 G4complex out = -gamma*(1. - a1*d <<
1601 <<
1602 return out; <<
1603 } <<
1604 <<
1605 ///////////////////////////////////////////// <<
1606 // <<
1607 // <<
1608 <<
1609 G4complex G4NuclNuclDiffuseElastic::Amplitud <<
1610 { <<
1611 G4double kappa = std::sqrt(0.5*fProfileLamb <<
1612 G4complex out = G4complex(kappa/fWaveVector <<
1613 <<
1614 out *= PhaseNear(theta); <<
1615 <<
1616 if( theta <= fRutherfordTheta ) <<
1617 { <<
1618 out *= GammaLess(theta) + ProfileNear(the <<
1619 // out *= GammaMore(theta) + ProfileNear( <<
1620 out += CoulombAmplitude(theta); <<
1621 } <<
1622 else <<
1623 { <<
1624 out *= GammaMore(theta) + ProfileNear(the <<
1625 // out *= GammaLess(theta) + ProfileNear( <<
1626 } <<
1627 return out; <<
1628 } <<
1629 <<
1630 ///////////////////////////////////////////// <<
1631 // <<
1632 // <<
1633 <<
1634 G4complex G4NuclNuclDiffuseElastic::Amplitude <<
1635 { <<
1636 G4double sinThetaR = 2.*fHalfRutThetaTg/(1 <<
1637 G4double dTheta = 0.5*(theta - fRutherf <<
1638 G4double sindTheta = std::sin(dTheta); <<
1639 G4double persqrt2 = std::sqrt(0.5); <<
1640 <<
1641 G4complex order = G4complex(persqrt2,pe <<
1642 order *= std::sqrt(0.5*fProfil <<
1643 // order *= std::sqrt(0.5*fPro <<
1644 <<
1645 G4complex out; <<
1646 <<
1647 if ( theta <= fRutherfordTheta ) <<
1648 { <<
1649 out = 1. - 0.5*GetErfcInt(-order)*Profile <<
1650 } <<
1651 else <<
1652 { <<
1653 out = 0.5*GetErfcInt(order)*ProfileNear(t <<
1654 } <<
1655 <<
1656 out *= CoulombAmplitude(theta); <<
1657 <<
1658 return out; <<
1659 } <<
1660 <<
1661 ///////////////////////////////////////////// <<
1662 // <<
1663 // <<
1664 <<
1665 G4complex G4NuclNuclDiffuseElastic::Amplitud <<
1666 { <<
1667 G4int n; <<
1668 G4double T12b, b, b2; // cosTheta = std::co <<
1669 G4complex out = G4complex(0.,0.), shiftC, s <<
1670 G4complex im = G4complex(0.,1.); <<
1671 <<
1672 for( n = 0; n < fMaxL; n++) <<
1673 { <<
1674 shiftC = std::exp( im*2.*CalculateCoulomb <<
1675 // b = ( fZommerfeld + std::sqrt( fZommer <<
1676 b = ( std::sqrt( G4double(n*(n+1)) ) )/fW <<
1677 b2 = b*b; <<
1678 T12b = fSumSigma*G4Exp(-b2/fNuclearRadius <<
1679 shiftN = std::exp( -0.5*(1.-im*fEtaRatio) <<
1680 out += (2.*n+1.)*shiftC*shiftN*GetLegend <<
1681 } <<
1682 out /= 2.*im*fWaveVector; <<
1683 out += CoulombAmplitude(theta); <<
1684 return out; <<
1685 } <<
1686 <<
1687 <<
1688 ///////////////////////////////////////////// <<
1689 // <<
1690 // <<
1691 <<
1692 G4complex G4NuclNuclDiffuseElastic::Amplitud <<
1693 { <<
1694 G4int n; <<
1695 G4double T12b, a, aTemp, b2, sinThetaH = st <<
1696 G4double sinThetaH2 = sinThetaH*sinThetaH; <<
1697 G4complex out = G4complex(0.,0.); <<
1698 G4complex im = G4complex(0.,1.); <<
1699 <<
1700 a = -fSumSigma/CLHEP::twopi/fNuclearRadius <<
1701 b2 = fWaveVector*fWaveVector*fNuclearRadius <<
1702 <<
1703 aTemp = a; <<
1704 <<
1705 for( n = 1; n < fMaxL; n++) <<
1706 { <<
1707 T12b = aTemp*G4Exp(-b2/n)/n; <<
1708 aTemp *= a; <<
1709 out += T12b; <<
1710 G4cout<<"out = "<<out<<G4endl; <<
1711 } <<
1712 out *= -4.*im*fWaveVector/CLHEP::pi; <<
1713 out += CoulombAmplitude(theta); <<
1714 return out; <<
1715 } <<
1716 <<
1717 <<
1718 ///////////////////////////////////////////// <<
1719 // <<
1720 // Test for given particle and element table <<
1721 // For the partMom in CMS. <<
1722 <<
1723 void G4NuclNuclDiffuseElastic::InitParameters <<
1724 G4d <<
1725 { <<
1726 fAtomicNumber = Z; // atomic number <<
1727 fAtomicWeight = A; // number of nucleo <<
1728 <<
1729 fNuclearRadius2 = CalculateNuclearRad(fAtom <<
1730 G4double A1 = G4double( theParticle->Ge <<
1731 fNuclearRadius1 = CalculateNuclearRad(A1); <<
1732 // fNuclearRadius = std::sqrt(fNuclearRadiu <<
1733 fNuclearRadius = fNuclearRadius1 + fNuclear <<
1734 <<
1735 G4double a = 0.; <<
1736 G4double z = theParticle->GetPDGCharge(); <<
1737 G4double m1 = theParticle->GetPDGMass(); <<
1738 <<
1739 fWaveVector = partMom/CLHEP::hbarc; <<
1740 <<
1741 G4double lambda = fCofLambda*fWaveVector*fN <<
1742 G4cout<<"kR = "<<lambda<<G4endl; <<
1743 <<
1744 if( z ) <<
1745 { <<
1746 a = partMom/m1; // beta*gamma f <<
1747 fBeta = a/std::sqrt(1+a*a); <<
1748 fZommerfeld = CalculateZommerfeld( fBeta, <<
1749 fRutherfordRatio = fZommerfeld/fWaveVecto <<
1750 fAm = CalculateAm( partMom, fZomm <<
1751 } <<
1752 G4cout<<"fZommerfeld = "<<fZommerfeld<<G4en <<
1753 fProfileLambda = lambda; // *std::sqrt(1.-2 <<
1754 G4cout<<"fProfileLambda = "<<fProfileLambda <<
1755 fProfileDelta = fCofDelta*fProfileLambda; <<
1756 fProfileAlpha = fCofAlpha*fProfileLambda; <<
1757 <<
1758 CalculateCoulombPhaseZero(); <<
1759 CalculateRutherfordAnglePar(); <<
1760 <<
1761 return; <<
1762 } <<
1763 ///////////////////////////////////////////// <<
1764 // <<
1765 // Test for given particle and element table <<
1766 // For the partMom in CMS. <<
1767 <<
1768 void G4NuclNuclDiffuseElastic::InitDynParamet <<
1769 G4d <<
1770 { <<
1771 G4double a = 0.; <<
1772 G4double z = theParticle->GetPDGCharge(); <<
1773 G4double m1 = theParticle->GetPDGMass(); <<
1774 <<
1775 fWaveVector = partMom/CLHEP::hbarc; <<
1776 <<
1777 G4double lambda = fCofLambda*fWaveVector*fN <<
1778 <<
1779 if( z ) <<
1780 { <<
1781 a = partMom/m1; // beta*gamma f <<
1782 fBeta = a/std::sqrt(1+a*a); <<
1783 fZommerfeld = CalculateZommerfeld( fBeta, <<
1784 fRutherfordRatio = fZommerfeld/fWaveVecto <<
1785 fAm = CalculateAm( partMom, fZomm <<
1786 } <<
1787 fProfileLambda = lambda; // *std::sqrt(1.-2 <<
1788 fProfileDelta = fCofDelta*fProfileLambda; <<
1789 fProfileAlpha = fCofAlpha*fProfileLambda; <<
1790 <<
1791 CalculateCoulombPhaseZero(); <<
1792 CalculateRutherfordAnglePar(); <<
1793 <<
1794 return; <<
1795 } <<
1796 <<
1797 <<
1798 ///////////////////////////////////////////// <<
1799 // <<
1800 // Test for given particle and element table <<
1801 // For the partMom in CMS. <<
1802 <<
1803 void G4NuclNuclDiffuseElastic::InitParameters <<
1804 G4d <<
1805 { <<
1806 fAtomicNumber = Z; // target atomic nu <<
1807 fAtomicWeight = A; // target number of <<
1808 <<
1809 fNuclearRadius2 = CalculateNuclearRad(fAtom <<
1810 G4double A1 = G4double( aParticle->GetD <<
1811 fNuclearRadius1 = CalculateNuclearRad(A1); <<
1812 fNuclearRadiusSquare = fNuclearRadius1*fNuc <<
1813 <<
1814 <<
1815 G4double a = 0., kR12; <<
1816 G4double z = aParticle->GetDefinition()->G <<
1817 G4double m1 = aParticle->GetDefinition()->G <<
1818 <<
1819 fWaveVector = partMom/CLHEP::hbarc; <<
1820 <<
1821 G4double pN = A1 - z; <<
1822 if( pN < 0. ) pN = 0.; <<
1823 <<
1824 G4double tN = A - Z; <<
1825 if( tN < 0. ) tN = 0.; <<
1826 <<
1827 G4double pTkin = aParticle->GetKineticEnerg <<
1828 pTkin /= A1; <<
1829 <<
1830 <<
1831 fSumSigma = (Z*z+pN*tN)*GetHadronNucleonXsc <<
1832 (z*tN+pN*Z)*GetHadronNucleonXsc <<
1833 <<
1834 G4cout<<"fSumSigma = "<<fSumSigma/CLHEP::mi <<
1835 G4cout<<"pi*R2 = "<<CLHEP::pi*fNuclearRadiu <<
1836 kR12 = fWaveVector*std::sqrt(fNuclearRadius <<
1837 G4cout<<"k*sqrt(R2) = "<<kR12<<" "<<G4endl; <<
1838 fMaxL = (G4int(kR12)+1)*4; <<
1839 G4cout<<"fMaxL = "<<fMaxL<<" "<<G4endl; <<
1840 <<
1841 if( z ) <<
1842 { <<
1843 a = partMom/m1; // beta*gamma <<
1844 fBeta = a/std::sqrt(1+a*a); <<
1845 fZommerfeld = CalculateZommerfeld( fBet <<
1846 fAm = CalculateAm( partMom, fZo <<
1847 } <<
1848 <<
1849 CalculateCoulombPhaseZero(); <<
1850 <<
1851 <<
1852 return; <<
1853 } <<
1854 <<
1855 <<
1856 ///////////////////////////////////////////// <<
1857 // <<
1858 // Returns nucleon-nucleon cross-section base <<
1859 // data from mainly http://wwwppds.ihep.su:80 <<
1860 // projectile nucleon is pParticle with pTkin <<
1861 <<
1862 G4double <<
1863 G4NuclNuclDiffuseElastic::GetHadronNucleonXsc <<
1864 <<
1865 <<
1866 { <<
1867 G4double xsection(0), /*Delta,*/ A0, B0; <<
1868 G4double hpXsc(0); <<
1869 G4double hnXsc(0); <<
1870 <<
1871 <<
1872 G4double targ_mass = tParticle->GetPDGM <<
1873 G4double proj_mass = pParticle->GetPDGM <<
1874 <<
1875 G4double proj_energy = proj_mass + pTkin; <<
1876 G4double proj_momentum = std::sqrt(pTkin*(p <<
1877 <<
1878 G4double sMand = CalcMandelstamS ( proj_mas <<
1879 <<
1880 sMand /= CLHEP::GeV*CLHEP::GeV; // <<
1881 proj_momentum /= CLHEP::GeV; <<
1882 proj_energy /= CLHEP::GeV; <<
1883 proj_mass /= CLHEP::GeV; <<
1884 G4double logS = G4Log(sMand); <<
1885 <<
1886 // General PDG fit constants <<
1887 <<
1888 <<
1889 // fEtaRatio=Re[f(0)]/Im[f(0)] <<
1890 <<
1891 if( proj_momentum >= 1.2 ) <<
1892 { <<
1893 fEtaRatio = 0.13*(logS - 5.8579332)*G4Po <<
1894 } <<
1895 else if( proj_momentum >= 0.6 ) <<
1896 { <<
1897 fEtaRatio = -75.5*(G4Pow::GetInstance()-> <<
1898 (G4Pow::GetInstance()->powA(3*proj_moment <<
1899 } <<
1900 else <<
1901 { <<
1902 fEtaRatio = 15.5*proj_momentum/(27*proj_m <<
1903 } <<
1904 G4cout<<"fEtaRatio = "<<fEtaRatio<<G4endl; <<
1905 <<
1906 // xsc <<
1907 <<
1908 if( proj_momentum >= 10. ) // high energy: <<
1909 // if( proj_momentum >= 2.) <<
1910 { <<
1911 //Delta = 1.; <<
1912 <<
1913 //if( proj_energy < 40. ) Delta = 0.916+0 <<
1914 <<
1915 //AR-12Aug2016 if( proj_momentum >= 10.) <<
1916 { <<
1917 B0 = 7.5; <<
1918 A0 = 100. - B0*G4Log(3.0e7); <<
1919 <<
1920 xsection = A0 + B0*G4Log(proj_energy) <<
1921 + 103*G4Pow::GetInstance()- <<
1922 0.93827*0.93827,-0.165); <<
1923 } <<
1924 } <<
1925 else // low energy pp = nn != np <<
1926 { <<
1927 if(pParticle == tParticle) // pp or nn <<
1928 { <<
1929 if( proj_momentum < 0.73 ) <<
1930 { <<
1931 hnXsc = 23 + 50*( G4Pow::GetInstanc <<
1932 } <<
1933 else if( proj_momentum < 1.05 ) <<
1934 { <<
1935 hnXsc = 23 + 40*(G4Log(proj_momentu <<
1936 (G4Log(proj_momentum <<
1937 } <<
1938 else // if( proj_momentum < 10. ) <<
1939 { <<
1940 hnXsc = 39.0 + <<
1941 75*(proj_momentum - 1.2)/(G4Pow <<
1942 } <<
1943 xsection = hnXsc; <<
1944 } <<
1945 else // pn to be np <<
1946 { <<
1947 if( proj_momentum < 0.8 ) <<
1948 { <<
1949 hpXsc = 33+30*G4Pow::GetInstance()- <<
1950 } <<
1951 else if( proj_momentum < 1.4 ) <<
1952 { <<
1953 hpXsc = 33+30*G4Pow::GetInstance()- <<
1954 } <<
1955 else // if( proj_momentum < 10. ) <<
1956 { <<
1957 hpXsc = 33.3+ <<
1958 20.8*(G4Pow::GetInstance()->pow <<
1959 (G4Pow::GetInstance()->powA( <<
1960 } <<
1961 xsection = hpXsc; <<
1962 } <<
1963 } <<
1964 xsection *= CLHEP::millibarn; // parametris <<
1965 G4cout<<"xsection = "<<xsection/CLHEP::mill <<
1966 return xsection; <<
1967 } <<
1968 <<
1969 ///////////////////////////////////////////// <<
1970 // <<
1971 // The ratio el/ruth for Fresnel smooth nucle <<
1972 <<
1973 G4double G4NuclNuclDiffuseElastic::GetRatioGe <<
1974 { <<
1975 G4double sinThetaR = 2.*fHalfRutThetaTg/(1 <<
1976 G4double dTheta = 0.5*(theta - fRutherf <<
1977 G4double sindTheta = std::sin(dTheta); <<
1978 <<
1979 G4double prof = Profile(theta); <<
1980 G4double prof2 = prof*prof; <<
1981 // G4double profmod = std::abs(prof); <<
1982 G4double order = std::sqrt(fProfileLam <<
1983 <<
1984 order = std::abs(order); // since sin chang <<
1985 // G4cout<<"order = "<<order<<G4endl; <<
1986 <<
1987 G4double cosFresnel = GetCint(order); <<
1988 G4double sinFresnel = GetSint(order); <<
1989 <<
1990 G4double out; <<
1991 <<
1992 if ( theta <= fRutherfordTheta ) <<
1993 { <<
1994 out = 1. + 0.5*( (0.5-cosFresnel)*(0.5-c <<
1995 out += ( cosFresnel + sinFresnel - 1. )*p <<
1996 } <<
1997 else <<
1998 { <<
1999 out = 0.5*( (0.5-cosFresnel)*(0.5-cosFres <<
2000 } <<
2001 <<
2002 return out; <<
2003 } <<
2004 <<
2005 ///////////////////////////////////////////// <<
2006 // <<
2007 // For the calculation of arg Gamma(z) one ne <<
2008 // of ln(Gamma(z)) <<
2009 <<
2010 G4complex G4NuclNuclDiffuseElastic::GammaLoga <<
2011 { <<
2012 const G4double cof[6] = { 76.18009172947146 <<
2013 24.0140982408309 <<
2014 0.1208650973866 <<
2015 G4int j; <<
2016 G4complex z = zz - 1.0; <<
2017 G4complex tmp = z + 5.5; <<
2018 tmp -= (z + 0.5) * std::log(tmp); <<
2019 G4complex ser = G4complex(1.000000000190015 <<
2020 <<
2021 for ( j = 0; j <= 5; j++ ) <<
2022 { <<
2023 z += 1.0; <<
2024 ser += cof[j]/z; <<
2025 } <<
2026 return -tmp + std::log(2.5066282746310005*s <<
2027 } <<
2028 <<
2029 ///////////////////////////////////////////// <<
2030 // <<
2031 // Bessel J0 function based on rational appro <<
2032 // J.F. Hart, Computer Approximations, New Yo <<
2033 <<
2034 G4double G4NuclNuclDiffuseElastic::BesselJzer <<
2035 { <<
2036 G4double modvalue, value2, fact1, fact2, ar <<
2037 <<
2038 modvalue = std::fabs(value); <<
2039 <<
2040 if ( value < 8.0 && value > -8.0 ) <<
2041 { <<
2042 value2 = value*value; <<
2043 <<
2044 fact1 = 57568490574.0 + value2*(-1336259 <<
2045 + value2*( 6516196 <<
2046 + value2*(-1121442 <<
2047 + value2*( 77392.3 <<
2048 + value2*(-184.905 <<
2049 <<
2050 fact2 = 57568490411.0 + value2*( 1029532 <<
2051 + value2*( 9494680 <<
2052 + value2*(59272.64 <<
2053 + value2*(267.8532 <<
2054 + value2*1.0 <<
2055 <<
2056 bessel = fact1/fact2; <<
2057 } <<
2058 else <<
2059 { <<
2060 arg = 8.0/modvalue; <<
2061 <<
2062 value2 = arg*arg; <<
2063 <<
2064 shift = modvalue-0.785398164; <<
2065 <<
2066 fact1 = 1.0 + value2*(-0.1098628627e-2 <<
2067 + value2*(0.2734510407e-4 <<
2068 + value2*(-0.2073370639e-5 <<
2069 + value2*0.2093887211e-6 <<
2070 <<
2071 fact2 = -0.1562499995e-1 + value2*(0.143 <<
2072 + value2*(-0.69 <<
2073 + value2*(0.762 <<
2074 - value2*0.9349 <<
2075 <<
2076 bessel = std::sqrt(0.636619772/modvalue)* <<
2077 } <<
2078 return bessel; <<
2079 } <<
2080 <<
2081 ///////////////////////////////////////////// <<
2082 // <<
2083 // Bessel J1 function based on rational appro <<
2084 // J.F. Hart, Computer Approximations, New Yo <<
2085 <<
2086 G4double G4NuclNuclDiffuseElastic::BesselJone <<
2087 { <<
2088 G4double modvalue, value2, fact1, fact2, ar <<
2089 <<
2090 modvalue = std::fabs(value); <<
2091 <<
2092 if ( modvalue < 8.0 ) <<
2093 { <<
2094 value2 = value*value; <<
2095 <<
2096 fact1 = value*(72362614232.0 + value2*(- <<
2097 + value2*( <<
2098 + value2*(- <<
2099 + value2*( <<
2100 + value2*(- <<
2101 <<
2102 fact2 = 144725228442.0 + value2*(2300535 <<
2103 + value2*(1858330 <<
2104 + value2*(99447.4 <<
2105 + value2*(376.999 <<
2106 + value2*1.0 <<
2107 bessel = fact1/fact2; <<
2108 } <<
2109 else <<
2110 { <<
2111 arg = 8.0/modvalue; <<
2112 <<
2113 value2 = arg*arg; <<
2114 <<
2115 shift = modvalue - 2.356194491; <<
2116 <<
2117 fact1 = 1.0 + value2*( 0.183105e-2 <<
2118 + value2*(-0.3516396496e-4 <<
2119 + value2*(0.2457520174e-5 <<
2120 + value2*(-0.240337019e-6 <<
2121 <<
2122 fact2 = 0.04687499995 + value2*(-0.200269 <<
2123 + value2*( 0.844919 <<
2124 + value2*(-0.882289 <<
2125 + value2*0.10578741 <<
2126 <<
2127 bessel = std::sqrt( 0.636619772/modvalue) <<
2128 <<
2129 if (value < 0.0) bessel = -bessel; <<
2130 } <<
2131 return bessel; <<
2132 } 1591 }
2133 1592
2134 // 1593 //
2135 // 1594 //
2136 ///////////////////////////////////////////// 1595 /////////////////////////////////////////////////////////////////////////////////
2137 1596