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