| 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 26
27 // G4 Low energy model: n-n or p-p scattering 27 // G4 Low energy model: n-n or p-p scattering
28 // F.W. Jones, L.G. Greeniaus, H.P. Wellisch 28 // F.W. Jones, L.G. Greeniaus, H.P. Wellisch
29 29
30 // FWJ 27-AUG-2010: extended Coulomb-suppresse <<
31 30
32 #include "G4LEpp.hh" 31 #include "G4LEpp.hh"
33 #include "G4PhysicalConstants.hh" <<
34 #include "G4SystemOfUnits.hh" <<
35 #include "Randomize.hh" 32 #include "Randomize.hh"
36 #include "G4ios.hh" 33 #include "G4ios.hh"
37 34
38 // Initialization of static data arrays: 35 // Initialization of static data arrays:
39 #include "G4LEppData.hh" 36 #include "G4LEppData.hh"
40 37
41 #include "G4PhysicsModelCatalog.hh" << 38 G4LEpp::G4LEpp():G4HadronicInteraction("G4LEpp")
>> 39 {
>> 40 // theParticleChange.SetNumberOfSecondaries(1);
>> 41 // SetMinEnergy(10.*MeV);
>> 42 // SetMaxEnergy(1200.*MeV);
42 43
>> 44 SetCoulombEffects(0);
43 45
44 G4LEpp::G4LEpp():G4HadronElastic("G4LEpp") <<
45 { <<
46 secID = G4PhysicsModelCatalog::GetModelID( " <<
47 SetMinEnergy(0.); 46 SetMinEnergy(0.);
48 SetMaxEnergy(5.*GeV); << 47 SetMaxEnergy(1200.*GeV);
49 } 48 }
50 49
51 G4LEpp::~G4LEpp() 50 G4LEpp::~G4LEpp()
52 {} << 51 {
>> 52 // theParticleChange.Clear();
>> 53 }
>> 54
>> 55
>> 56 void
>> 57 G4LEpp::SetCoulombEffects(G4int State)
>> 58 {
>> 59 if (State) {
>> 60 for(G4int i=0; i<NANGLE; i++)
>> 61 {
>> 62 sig[i] = SigCoul[i];
>> 63 }
>> 64 elab = ElabCoul;
>> 65 }
>> 66 else {
>> 67 for(G4int i=0; i<NANGLE; i++)
>> 68 {
>> 69 sig[i] = Sig[i];
>> 70 }
>> 71 elab = Elab;
>> 72 }
>> 73 }
>> 74
53 75
54 G4HadFinalState* 76 G4HadFinalState*
55 G4LEpp::ApplyYourself(const G4HadProjectile& a 77 G4LEpp::ApplyYourself(const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
56 { 78 {
57 theParticleChange.Clear(); 79 theParticleChange.Clear();
58 const G4HadProjectile* aParticle = &aTrack 80 const G4HadProjectile* aParticle = &aTrack;
59 81
60 G4double P = aParticle->GetTotalMomentum() 82 G4double P = aParticle->GetTotalMomentum();
61 G4double Px = aParticle->Get4Momentum().x( 83 G4double Px = aParticle->Get4Momentum().x();
62 G4double Py = aParticle->Get4Momentum().y( 84 G4double Py = aParticle->Get4Momentum().y();
63 G4double Pz = aParticle->Get4Momentum().z( 85 G4double Pz = aParticle->Get4Momentum().z();
64 G4double E = aParticle->GetTotalEnergy(); << 86 G4double ek = aParticle->GetKineticEnergy();
65 G4ThreeVector theInitial = aParticle->Get4 << 87 G4ThreeVector theInitial = aParticle->Get4Momentum().vect();
66 88
67 if (verboseLevel > 1) { << 89 // if (verboseLevel > 1)
68 G4double ek = aParticle->GetKineticEnerg << 90 {
>> 91 G4double E = aParticle->GetTotalEnergy();
69 G4double E0 = aParticle->GetDefinition() 92 G4double E0 = aParticle->GetDefinition()->GetPDGMass();
70 G4double Q = aParticle->GetDefinition()- 93 G4double Q = aParticle->GetDefinition()->GetPDGCharge();
71 G4int A = targetNucleus.GetA_asInt(); << 94 G4double N = targetNucleus.GetN();
72 G4int Z = targetNucleus.GetZ_asInt(); << 95 G4double Z = targetNucleus.GetZ();
73 G4cout << "G4LEpp:ApplyYourself: inciden 96 G4cout << "G4LEpp:ApplyYourself: incident particle: "
74 << aParticle->GetDefinition()->Ge 97 << aParticle->GetDefinition()->GetParticleName() << G4endl;
75 G4cout << "P = " << P/GeV << " GeV/c" 98 G4cout << "P = " << P/GeV << " GeV/c"
76 << ", Px = " << Px/GeV << " GeV/c 99 << ", Px = " << Px/GeV << " GeV/c"
77 << ", Py = " << Py/GeV << " GeV/c 100 << ", Py = " << Py/GeV << " GeV/c"
78 << ", Pz = " << Pz/GeV << " GeV/c 101 << ", Pz = " << Pz/GeV << " GeV/c" << G4endl;
79 G4cout << "E = " << E/GeV << " GeV" 102 G4cout << "E = " << E/GeV << " GeV"
80 << ", kinetic energy = " << ek/Ge 103 << ", kinetic energy = " << ek/GeV << " GeV"
81 << ", mass = " << E0/GeV << " GeV 104 << ", mass = " << E0/GeV << " GeV"
82 << ", charge = " << Q << G4endl; 105 << ", charge = " << Q << G4endl;
83 G4cout << "G4LEpp:ApplyYourself: materia 106 G4cout << "G4LEpp:ApplyYourself: material:" << G4endl;
84 G4cout << "A = " << A << 107 G4cout << "A = " << N
85 << ", Z = " << Z 108 << ", Z = " << Z
86 << ", atomic mass " 109 << ", atomic mass "
87 << G4Proton::Proton()->GetPDGMas 110 << G4Proton::Proton()->GetPDGMass()/GeV << "GeV"
88 << G4endl; 111 << G4endl;
89 // 112 //
90 // GHEISHA ADD operation to get total en 113 // GHEISHA ADD operation to get total energy, mass, charge
91 // 114 //
92 E += proton_mass_c2; << 115 E += G4Proton::Proton()->GetPDGMass();
93 G4double E02 = E*E - P*P; 116 G4double E02 = E*E - P*P;
94 E0 = std::sqrt(std::fabs(E02)); << 117 E0 = std::sqrt(std::abs(E02));
95 if (E02 < 0)E0 *= -1; 118 if (E02 < 0)E0 *= -1;
96 Q += Z; 119 Q += Z;
97 G4cout << "G4LEpp:ApplyYourself: total:" 120 G4cout << "G4LEpp:ApplyYourself: total:" << G4endl;
98 G4cout << "E = " << E/GeV << " GeV" 121 G4cout << "E = " << E/GeV << " GeV"
99 << ", mass = " << E0/GeV << " GeV 122 << ", mass = " << E0/GeV << " GeV"
100 << ", charge = " << Q << G4endl; 123 << ", charge = " << Q << G4endl;
101 } 124 }
102 G4double t = SampleInvariantT(aParticle->G << 125
103 G4double cost = 1.0 - 2*t/(P*P); << 126 // Find energy bin
104 if(cost > 1.0) { cost = 1.0; } << 127
105 if(cost <-1.0) { cost =-1.0; } << 128 G4int je1 = 0;
106 G4double sint = std::sqrt((1.0 - cost)*(1. << 129 G4int je2 = NENERGY - 1;
107 G4double phi = twopi*G4UniformRand(); << 130 ek = ek/GeV;
>> 131 do {
>> 132 G4int midBin = (je1 + je2)/2;
>> 133 if (ek < elab[midBin])
>> 134 je2 = midBin;
>> 135 else
>> 136 je1 = midBin;
>> 137 } while (je2 - je1 > 1);
>> 138 // G4int j;
>> 139 //std::abs(ek-elab[je1]) < std::abs(ek-elab[je2]) ? j = je1 : j = je2;
>> 140 G4double delab = elab[je2] - elab[je1];
>> 141
>> 142 // Sample the angle
>> 143
>> 144 G4float sample = G4UniformRand();
>> 145 G4int ke1 = 0;
>> 146 G4int ke2 = NANGLE - 1;
>> 147 G4double dsig = sig[je2][0] - sig[je1][0];
>> 148 G4double rc = dsig/delab;
>> 149 G4double b = sig[je1][0] - rc*elab[je1];
>> 150 G4double sigint1 = rc*ek + b;
>> 151 G4double sigint2 = 0.;
>> 152
>> 153 if (verboseLevel > 1) G4cout << "sample=" << sample << G4endl
>> 154 << ke1 << " " << ke2 << " "
>> 155 << sigint1 << " " << sigint2 << G4endl;
>> 156
>> 157 do {
>> 158 G4int midBin = (ke1 + ke2)/2;
>> 159 dsig = sig[je2][midBin] - sig[je1][midBin];
>> 160 rc = dsig/delab;
>> 161 b = sig[je1][midBin] - rc*elab[je1];
>> 162 G4double sigint = rc*ek + b;
>> 163 if (sample < sigint) {
>> 164 ke2 = midBin;
>> 165 sigint2 = sigint;
>> 166 }
>> 167 else {
>> 168 ke1 = midBin;
>> 169 sigint1 = sigint;
>> 170 }
>> 171 if (verboseLevel > 1)G4cout << ke1 << " " << ke2 << " "
>> 172 << sigint1 << " " << sigint2 << G4endl;
>> 173 } while (ke2 - ke1 > 1);
>> 174
>> 175 // sigint1 and sigint2 should be recoverable from above loop
>> 176
>> 177 // G4double dsig = sig[je2][ke1] - sig[je1][ke1];
>> 178 // G4double rc = dsig/delab;
>> 179 // G4double b = sig[je1][ke1] - rc*elab[je1];
>> 180 // G4double sigint1 = rc*ek + b;
>> 181
>> 182 // G4double dsig = sig[je2][ke2] - sig[je1][ke2];
>> 183 // G4double rc = dsig/delab;
>> 184 // G4double b = sig[je1][ke2] - rc*elab[je1];
>> 185 // G4double sigint2 = rc*ek + b;
>> 186
>> 187 dsig = sigint2 - sigint1;
>> 188 rc = 1./dsig;
>> 189 b = ke1 - rc*sigint1;
>> 190 G4double kint = rc*sample + b;
>> 191 G4double theta = (0.5 + kint)*pi/180.;
>> 192 if (theta < 0.) theta = 0.;
>> 193
>> 194 // G4int k;
>> 195 //std::abs(sample-sig[j][ke1]) < std::abs(sample-sig[j][ke2]) ? k = ke1 : k = ke2;
>> 196 // G4double theta = (0.5 + k)*pi/180.;
>> 197
>> 198 if (verboseLevel > 1) {
>> 199 G4cout << " energy bin " << je1 << " energy=" << elab[je1] << G4endl;
>> 200 G4cout << " angle bin " << kint << " angle=" << theta/degree << G4endl;
>> 201 }
>> 202
>> 203
108 // Get the target particle 204 // Get the target particle
>> 205
109 G4DynamicParticle* targetParticle = target 206 G4DynamicParticle* targetParticle = targetNucleus.ReturnTargetParticle();
110 207
111 G4double E1 = aParticle->GetTotalEnergy(); 208 G4double E1 = aParticle->GetTotalEnergy();
112 G4double M1 = aParticle->GetDefinition()-> 209 G4double M1 = aParticle->GetDefinition()->GetPDGMass();
113 G4double E2 = targetParticle->GetTotalEner 210 G4double E2 = targetParticle->GetTotalEnergy();
114 G4double M2 = targetParticle->GetDefinitio 211 G4double M2 = targetParticle->GetDefinition()->GetPDGMass();
115 G4double totalEnergy = E1 + E2; 212 G4double totalEnergy = E1 + E2;
116 G4double pseudoMass = std::sqrt(totalEnerg 213 G4double pseudoMass = std::sqrt(totalEnergy*totalEnergy - P*P);
>> 214 // pseudoMass also = std::sqrt(M1*M1 + M2*M2 + 2*M2*E1)
117 215
118 // Transform into centre of mass system 216 // Transform into centre of mass system
119 217
120 G4double px = (M2/pseudoMass)*Px; 218 G4double px = (M2/pseudoMass)*Px;
121 G4double py = (M2/pseudoMass)*Py; 219 G4double py = (M2/pseudoMass)*Py;
122 G4double pz = (M2/pseudoMass)*Pz; 220 G4double pz = (M2/pseudoMass)*Pz;
123 G4double p = std::sqrt(px*px + py*py + pz* 221 G4double p = std::sqrt(px*px + py*py + pz*pz);
124 222
125 if (verboseLevel > 1) { 223 if (verboseLevel > 1) {
126 G4cout << " E1, M1 (GeV) " << E1/GeV << 224 G4cout << " E1, M1 (GeV) " << E1/GeV << " " << M1/GeV << G4endl;
127 G4cout << " E2, M2 (GeV) " << E2/GeV << 225 G4cout << " E2, M2 (GeV) " << E2/GeV << " " << M2/GeV << G4endl;
128 G4cout << " particle 1 momentum in CM << 226 G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
129 << " " << py/GeV << " " << 227 << pz/GeV << " " << p/GeV << G4endl;
130 << pz/GeV << " " << p/GeV << G4endl; <<
131 } 228 }
132 229
133 // First scatter w.r.t. Z axis 230 // First scatter w.r.t. Z axis
134 G4double pxnew = p*sint*std::cos(phi); << 231 G4double phi = G4UniformRand()*twopi;
135 G4double pynew = p*sint*std::sin(phi); << 232 G4double pxnew = p*std::sin(theta)*std::cos(phi);
136 G4double pznew = p*cost; << 233 G4double pynew = p*std::sin(theta)*std::sin(phi);
>> 234 G4double pznew = p*std::cos(theta);
137 235
138 // Rotate according to the direction of th 236 // Rotate according to the direction of the incident particle
139 if (px*px + py*py > 0) { 237 if (px*px + py*py > 0) {
140 G4double ph, cosp, sinp; << 238 G4double cost, sint, ph, cosp, sinp;
141 cost = pz/p; 239 cost = pz/p;
142 sint = (std::sqrt((1-cost)*(1+cost)) + s << 240 sint = (std::sqrt(std::abs((1-cost)*(1+cost))) + std::sqrt(px*px+py*py)/p)/2;
143 py < 0 ? ph = 3*halfpi : ph = halfpi; 241 py < 0 ? ph = 3*halfpi : ph = halfpi;
144 if (std::fabs(px) > 0.000001*GeV) ph = s << 242 if (std::abs(px) > 0.000001*GeV) ph = std::atan2(py,px);
145 cosp = std::cos(ph); 243 cosp = std::cos(ph);
146 sinp = std::sin(ph); 244 sinp = std::sin(ph);
147 px = (cost*cosp*pxnew - sinp*pynew + sin 245 px = (cost*cosp*pxnew - sinp*pynew + sint*cosp*pznew);
148 py = (cost*sinp*pxnew + cosp*pynew + sin 246 py = (cost*sinp*pxnew + cosp*pynew + sint*sinp*pznew);
149 pz = (-sint*pxnew + cos 247 pz = (-sint*pxnew + cost*pznew);
>> 248 // G4ThreeVector it(a,b,c);
>> 249 // p0->SetMomentum(it);
>> 250 // G4ThreeVector aTargetMom = theInitial - it;
>> 251 // targetParticle->SetMomentum(aTargetMom);
150 } 252 }
151 else { 253 else {
152 px = pxnew; 254 px = pxnew;
153 py = pynew; 255 py = pynew;
154 pz = pznew; 256 pz = pznew;
155 } 257 }
156 258
157 if (verboseLevel > 1) { 259 if (verboseLevel > 1) {
158 G4cout << " AFTER SCATTER..." << G4endl 260 G4cout << " AFTER SCATTER..." << G4endl;
159 G4cout << " particle 1 momentum in CM " 261 G4cout << " particle 1 momentum in CM " << px/GeV << " " << py/GeV << " "
160 << pz/GeV << " " << p/GeV << G4endl 262 << pz/GeV << " " << p/GeV << G4endl;
161 } 263 }
162 264
163 // Transform to lab system 265 // Transform to lab system
164 266
165 G4double E1pM2 = E1 + M2; 267 G4double E1pM2 = E1 + M2;
166 G4double betaCM = P/E1pM2; 268 G4double betaCM = P/E1pM2;
167 G4double betaCMx = Px/E1pM2; 269 G4double betaCMx = Px/E1pM2;
168 G4double betaCMy = Py/E1pM2; 270 G4double betaCMy = Py/E1pM2;
169 G4double betaCMz = Pz/E1pM2; 271 G4double betaCMz = Pz/E1pM2;
170 G4double gammaCM = E1pM2/std::sqrt(E1pM2*E 272 G4double gammaCM = E1pM2/std::sqrt(E1pM2*E1pM2 - P*P);
171 273
172 if (verboseLevel > 1) { 274 if (verboseLevel > 1) {
173 G4cout << " betaCM " << betaCMx << " " 275 G4cout << " betaCM " << betaCMx << " " << betaCMy << " "
174 << betaCMz << " " << betaCM << G4 276 << betaCMz << " " << betaCM << G4endl;
175 G4cout << " gammaCM " << gammaCM << G4e 277 G4cout << " gammaCM " << gammaCM << G4endl;
176 } 278 }
177 279
178 // Now following GLOREN... 280 // Now following GLOREN...
179 281
180 G4double BETA[5], PA[5], PB[5]; 282 G4double BETA[5], PA[5], PB[5];
181 BETA[1] = -betaCMx; 283 BETA[1] = -betaCMx;
182 BETA[2] = -betaCMy; 284 BETA[2] = -betaCMy;
183 BETA[3] = -betaCMz; 285 BETA[3] = -betaCMz;
184 BETA[4] = gammaCM; 286 BETA[4] = gammaCM;
185 287
186 //The incident particle... 288 //The incident particle...
187 289
188 PA[1] = px; 290 PA[1] = px;
189 PA[2] = py; 291 PA[2] = py;
190 PA[3] = pz; 292 PA[3] = pz;
191 PA[4] = std::sqrt(M1*M1 + p*p); 293 PA[4] = std::sqrt(M1*M1 + p*p);
192 294
193 G4double BETPA = BETA[1]*PA[1] + BETA[2]* 295 G4double BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
194 G4double BPGAM = (BETPA * BETA[4]/(BETA[4 296 G4double BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
195 297
196 PB[1] = PA[1] + BPGAM * BETA[1]; 298 PB[1] = PA[1] + BPGAM * BETA[1];
197 PB[2] = PA[2] + BPGAM * BETA[2]; 299 PB[2] = PA[2] + BPGAM * BETA[2];
198 PB[3] = PA[3] + BPGAM * BETA[3]; 300 PB[3] = PA[3] + BPGAM * BETA[3];
199 PB[4] = (PA[4] - BETPA) * BETA[4]; 301 PB[4] = (PA[4] - BETPA) * BETA[4];
200 302
201 G4DynamicParticle* newP = new G4DynamicPar 303 G4DynamicParticle* newP = new G4DynamicParticle;
202 newP->SetDefinition(aParticle->GetDefiniti << 304 newP->SetDefinition(const_cast<G4ParticleDefinition *>(aParticle->GetDefinition()) );
203 newP->SetMomentum(G4ThreeVector(PB[1], PB[ 305 newP->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
204 306
>> 307
205 //The target particle... 308 //The target particle...
206 309
207 PA[1] = -px; 310 PA[1] = -px;
208 PA[2] = -py; 311 PA[2] = -py;
209 PA[3] = -pz; 312 PA[3] = -pz;
210 PA[4] = std::sqrt(M2*M2 + p*p); 313 PA[4] = std::sqrt(M2*M2 + p*p);
211 314
212 BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + B 315 BETPA = BETA[1]*PA[1] + BETA[2]*PA[2] + BETA[3]*PA[3];
213 BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - 316 BPGAM = (BETPA * BETA[4]/(BETA[4] + 1.) - PA[4]) * BETA[4];
214 317
215 PB[1] = PA[1] + BPGAM * BETA[1]; 318 PB[1] = PA[1] + BPGAM * BETA[1];
216 PB[2] = PA[2] + BPGAM * BETA[2]; 319 PB[2] = PA[2] + BPGAM * BETA[2];
217 PB[3] = PA[3] + BPGAM * BETA[3]; 320 PB[3] = PA[3] + BPGAM * BETA[3];
218 PB[4] = (PA[4] - BETPA) * BETA[4]; 321 PB[4] = (PA[4] - BETPA) * BETA[4];
219 322
220 targetParticle->SetMomentum(G4ThreeVector( 323 targetParticle->SetMomentum(G4ThreeVector(PB[1], PB[2], PB[3]));
221 324
>> 325 // G4double ektotal = newP->GetKineticEnergy() +
>> 326 // targetParticle->GetKineticEnergy();
>> 327
222 if (verboseLevel > 1) { 328 if (verboseLevel > 1) {
223 G4cout << " particle 1 momentum in LAB 329 G4cout << " particle 1 momentum in LAB "
224 << newP->GetMomentum()/GeV << 330 << newP->GetMomentum()*(1./GeV)
225 << " " << newP->GetTotalMomentum()/ 331 << " " << newP->GetTotalMomentum()/GeV << G4endl;
226 G4cout << " particle 2 momentum in LAB 332 G4cout << " particle 2 momentum in LAB "
227 << targetParticle->GetMomentum()/Ge << 333 << targetParticle->GetMomentum()*(1./GeV)
228 << " " << targetParticle->GetTotalM 334 << " " << targetParticle->GetTotalMomentum()/GeV << G4endl;
229 G4cout << " TOTAL momentum in LAB " 335 G4cout << " TOTAL momentum in LAB "
230 << (newP->GetMomentum()+targetParti << 336 << (newP->GetMomentum()+targetParticle->GetMomentum())*(1./GeV)
231 << " " 337 << " "
232 << (newP->GetMomentum()+targetParti 338 << (newP->GetMomentum()+targetParticle->GetMomentum()).mag()/GeV
233 << G4endl; 339 << G4endl;
234 } 340 }
235 341
236 theParticleChange.SetMomentumChange( newP- << 342 // if (theta < pi/2.) {
237 theParticleChange.SetEnergyChange(newP->Ge << 343 // G4double p = newP->GetMomentum().mag();
238 delete newP; << 344 // G4ThreeVector m = newP->GetMomentum();
>> 345 // if (p > DBL_MIN)
>> 346 // theParticleChange.SetMomentumChange(m.x()/p, m.y()/p, m.z()/p);
>> 347 // else
>> 348 // theParticleChange.SetMomentumChange(0., 0., 0.);
>> 349
>> 350 theParticleChange.SetMomentumChange( newP->GetMomentumDirection());
>> 351 theParticleChange.SetEnergyChange(newP->GetKineticEnergy());
>> 352 delete newP;
>> 353
>> 354 // }
>> 355 // else {
>> 356 // // charge exchange
>> 357 // theParticleChange.SetNumberOfSecondaries(2);
>> 358 // theParticleChange.AddSecondary(newP);
>> 359 // theParticleChange.SetStatusChange(fStopAndKill);
>> 360 // // theParticleChange.SetEnergyChange(0.0);
>> 361 // }
239 362
240 // Recoil particle 363 // Recoil particle
241 theParticleChange.AddSecondary(targetParti << 364 G4DynamicParticle* p1 = new G4DynamicParticle;
>> 365 p1->SetDefinition(targetParticle->GetDefinition());
>> 366 p1->SetMomentum(targetParticle->GetMomentum());
>> 367 theParticleChange.AddSecondary(p1);
>> 368
>> 369
242 return &theParticleChange; 370 return &theParticleChange;
243 } 371 }
244 372
245 ////////////////////////////////////////////// << 373 // end of file
246 // <<
247 // sample momentum transfer using Lab. momentu <<
248 <<
249 G4double G4LEpp::SampleInvariantT(const G4Part <<
250 G4double plab, G4int , G4int ) <<
251 { <<
252 G4double nMass = p->GetPDGMass(); // 939.565 <<
253 G4double ek = std::sqrt(plab*plab+nMass*nMas <<
254 <<
255 // Find energy bin <<
256 <<
257 G4int je1 = 0; <<
258 G4int je2 = NENERGY - 1; <<
259 ek /= GeV; <<
260 <<
261 do <<
262 { <<
263 G4int midBin = (je1 + je2)/2; <<
264 <<
265 if (ek < elab[midBin]) je2 = midBin; <<
266 else je1 = midBin; <<
267 } <<
268 while (je2 - je1 > 1); /* Loop checking, 10 <<
269 <<
270 G4double delab = elab[je2] - elab[je1]; <<
271 <<
272 // Sample the angle <<
273 <<
274 G4double sample = G4UniformRand(); <<
275 G4int ke1 = 0; <<
276 G4int ke2 = NANGLE - 1; <<
277 G4double dsig, b, rc; <<
278 <<
279 dsig = Sig[je2][0] - Sig[je1][0]; <<
280 rc = dsig/delab; <<
281 b = Sig[je1][0] - rc*elab[je1]; <<
282 <<
283 G4double sigint1 = rc*ek + b; <<
284 G4double sigint2 = 0.; <<
285 <<
286 do <<
287 { <<
288 G4int midBin = (ke1 + ke2)/2; <<
289 dsig = Sig[je2][midBin] - Sig[je1][midBi <<
290 rc = dsig/delab; <<
291 b = Sig[je1][midBin] - rc*elab[je1]; <<
292 G4double sigint = rc*ek + b; <<
293 <<
294 if (sample < sigint) <<
295 { <<
296 ke2 = midBin; <<
297 sigint2 = sigint; <<
298 } <<
299 else <<
300 { <<
301 ke1 = midBin; <<
302 sigint1 = sigint; <<
303 } <<
304 } <<
305 while (ke2 - ke1 > 1); /* Loop checking, 10 <<
306 <<
307 dsig = sigint2 - sigint1; <<
308 rc = 1./dsig; <<
309 b = ke1 - rc*sigint1; <<
310 <<
311 G4double kint = rc*sample + b; <<
312 G4double theta = (0.5 + kint)*pi/180.; <<
313 G4double t = 0.5*plab*plab*(1 - std::cos(the <<
314 <<
315 return t; <<
316 } <<
317 // end of file <<
318 374