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