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