| 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 // Geant4 Header : G4NeutronElectronElModel 27 // Geant4 Header : G4NeutronElectronElModel
28 // 28 //
29 // 16.5.17: V.Grichine 29 // 16.5.17: V.Grichine
30 // 30 //
31 31
32 #include "G4NeutronElectronElModel.hh" 32 #include "G4NeutronElectronElModel.hh"
33 #include "G4SystemOfUnits.hh" 33 #include "G4SystemOfUnits.hh"
34 #include "G4ParticleTable.hh" 34 #include "G4ParticleTable.hh"
35 #include "G4ParticleDefinition.hh" 35 #include "G4ParticleDefinition.hh"
36 #include "G4IonTable.hh" 36 #include "G4IonTable.hh"
37 #include "Randomize.hh" 37 #include "Randomize.hh"
38 #include "G4Integrator.hh" 38 #include "G4Integrator.hh"
39 #include "G4Electron.hh" 39 #include "G4Electron.hh"
40 #include "G4PhysicsTable.hh" 40 #include "G4PhysicsTable.hh"
41 #include "G4PhysicsLogVector.hh" 41 #include "G4PhysicsLogVector.hh"
42 #include "G4PhysicsFreeVector.hh" 42 #include "G4PhysicsFreeVector.hh"
43 #include "G4PhysicsModelCatalog.hh" 43 #include "G4PhysicsModelCatalog.hh"
44 44
45 45
46 using namespace std; 46 using namespace std;
47 using namespace CLHEP; 47 using namespace CLHEP;
48 48
49 G4NeutronElectronElModel::G4NeutronElectronElM 49 G4NeutronElectronElModel::G4NeutronElectronElModel(const G4String& name)
50 : G4HadronElastic(name) 50 : G4HadronElastic(name)
51 { 51 {
52 secID = G4PhysicsModelCatalog::GetModelID( " 52 secID = G4PhysicsModelCatalog::GetModelID( "model_" + name );
53 53
54 // neutron magneton squared 54 // neutron magneton squared
55 55
56 fM = neutron_mass_c2; // neutron mass 56 fM = neutron_mass_c2; // neutron mass
57 fM2 = fM*fM; 57 fM2 = fM*fM;
58 fme = electron_mass_c2; 58 fme = electron_mass_c2;
59 fme2 = fme*fme; 59 fme2 = fme*fme;
60 fMv2 = 0.7056*GeV*GeV; 60 fMv2 = 0.7056*GeV*GeV;
61 61
62 SetMinEnergy( 0.001*GeV ); 62 SetMinEnergy( 0.001*GeV );
63 SetMaxEnergy( 10.*TeV ); 63 SetMaxEnergy( 10.*TeV );
64 SetLowestEnergyLimit(1.e-6*eV); 64 SetLowestEnergyLimit(1.e-6*eV);
65 65
66 theElectron = G4Electron::Electron(); 66 theElectron = G4Electron::Electron();
67 // PDG2016: sin^2 theta Weinberg 67 // PDG2016: sin^2 theta Weinberg
68 68
69 fEnergyBin = 200; 69 fEnergyBin = 200;
70 fMinEnergy = 1.*MeV; 70 fMinEnergy = 1.*MeV;
71 fMaxEnergy = 10000.*GeV; 71 fMaxEnergy = 10000.*GeV;
72 fEnergyVector = new G4PhysicsLogVector(fMinE << 72 fEnergyVector = new G4PhysicsLogVector(fMinEnergy, fMaxEnergy, fEnergyBin);
73 73
74 fAngleBin = 500; 74 fAngleBin = 500;
75 fAngleTable = 0; 75 fAngleTable = 0;
76 76
77 fCutEnergy = 0.; // default value 77 fCutEnergy = 0.; // default value
78 78
79 Initialise(); 79 Initialise();
80 } 80 }
81 81
82 ////////////////////////////////////////////// 82 ////////////////////////////////////////////////
83 83
84 G4NeutronElectronElModel::~G4NeutronElectronEl 84 G4NeutronElectronElModel::~G4NeutronElectronElModel()
85 { 85 {
86 if( fEnergyVector ) 86 if( fEnergyVector )
87 { 87 {
88 delete fEnergyVector; 88 delete fEnergyVector;
89 fEnergyVector = nullptr; << 89 fEnergyVector = 0;
90 } 90 }
91 if( fAngleTable ) 91 if( fAngleTable )
92 { 92 {
93 fAngleTable->clearAndDestroy(); 93 fAngleTable->clearAndDestroy();
94 delete fAngleTable; 94 delete fAngleTable;
95 fAngleTable = nullptr; 95 fAngleTable = nullptr;
96 } 96 }
97 } 97 }
98 98
99 ///////////////////////////////////////// 99 /////////////////////////////////////////
100 100
101 void G4NeutronElectronElModel::ModelDescriptio 101 void G4NeutronElectronElModel::ModelDescription(std::ostream& outFile) const
102 { 102 {
103 outFile << "G4NeutronElectronElModel is a ne << 103
104 << "model which uses the standard model \n << 104 outFile << "G4NeutronElectronElModel is a neutrino-electron (neutral current) elastic scattering\n"
105 << "transfer parameterization. The model << 105 << "model which uses the standard model \n"
>> 106 << "transfer parameterization. The model is fully relativistic\n";
>> 107
106 } 108 }
107 109
108 ////////////////////////////////////////////// 110 /////////////////////////////////////////////////////////
109 111
110 G4bool G4NeutronElectronElModel::IsApplicable( << 112 G4bool G4NeutronElectronElModel::IsApplicable(const G4HadProjectile & aTrack,
>> 113 G4Nucleus & targetNucleus)
111 { 114 {
>> 115 G4bool result = false;
112 G4String pName = aTrack.GetDefinition()->Get 116 G4String pName = aTrack.GetDefinition()->GetParticleName();
>> 117 // G4double minEnergy = 0.;
113 G4double energy = aTrack.GetTotalEnergy(); 118 G4double energy = aTrack.GetTotalEnergy();
114 119
115 return (pName == "neutron" && energy >= fMin << 120 if( fCutEnergy > 0. ) // min detected recoil electron energy
>> 121 {
>> 122 // minEnergy = 0.5*(fCutEnergy+sqrt(fCutEnergy*(fCutEnergy+2.*electron_mass_c2)));
>> 123 }
>> 124 if( pName == "neutron" &&
>> 125 energy >= fMinEnergy && energy <= fMaxEnergy )
>> 126 {
>> 127 result = true;
>> 128 }
>> 129 G4int Z = targetNucleus.GetZ_asInt();
>> 130 Z *= 1;
>> 131
>> 132 return result;
116 } 133 }
117 134
118 ////////////////////////////////////////////// 135 ////////////////////////////////////////////////////
119 136
120 void G4NeutronElectronElModel::Initialise() 137 void G4NeutronElectronElModel::Initialise()
121 { 138 {
122 G4double result = 0., sum, Tkin, dt, t1, t2; 139 G4double result = 0., sum, Tkin, dt, t1, t2;
123 G4int iTkin, jTransfer; 140 G4int iTkin, jTransfer;
124 G4Integrator<G4NeutronElectronElModel, G4dou 141 G4Integrator<G4NeutronElectronElModel, G4double(G4NeutronElectronElModel::*)(G4double)> integral;
125 142
126 fAngleTable = new G4PhysicsTable(fEnergyBin) 143 fAngleTable = new G4PhysicsTable(fEnergyBin);
127 144
128 for( iTkin = 0; iTkin < fEnergyBin; iTkin++) 145 for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
129 { 146 {
130 Tkin = fEnergyVector->GetLowEdgeEnergy(iT 147 Tkin = fEnergyVector->GetLowEdgeEnergy(iTkin);
131 fAm = CalculateAm(Tkin); 148 fAm = CalculateAm(Tkin);
132 dt = 1./fAngleBin; 149 dt = 1./fAngleBin;
133 150
134 G4PhysicsFreeVector* vectorT = new G4Physi 151 G4PhysicsFreeVector* vectorT = new G4PhysicsFreeVector(fAngleBin);
135 152
136 sum = 0.; 153 sum = 0.;
137 154
138 for( jTransfer = 0; jTransfer < fAngleBin; 155 for( jTransfer = 0; jTransfer < fAngleBin; jTransfer++)
139 { 156 {
140 t1 = dt*jTransfer; 157 t1 = dt*jTransfer;
141 t2 = t1 + dt; 158 t2 = t1 + dt;
142 159
143 result = integral.Legendre96( this, &G4N 160 result = integral.Legendre96( this, &G4NeutronElectronElModel::XscIntegrand, t1, t2 );
144 161
145 sum += result; 162 sum += result;
146 // G4cout<<sum<<", "; 163 // G4cout<<sum<<", ";
147 vectorT->PutValue(jTransfer, t1, sum); 164 vectorT->PutValue(jTransfer, t1, sum);
148 } 165 }
149 // G4cout<<G4endl; 166 // G4cout<<G4endl;
150 fAngleTable->insertAt(iTkin,vectorT); 167 fAngleTable->insertAt(iTkin,vectorT);
151 } 168 }
152 return; 169 return;
153 } 170 }
154 171
155 ////////////////////////////////////////////// 172 //////////////////////////////////////////////////////
156 // 173 //
157 // sample recoil electron energy in lab frame 174 // sample recoil electron energy in lab frame
158 175
159 G4double G4NeutronElectronElModel::SampleSin2H 176 G4double G4NeutronElectronElModel::SampleSin2HalfTheta(G4double Tkin)
160 { 177 {
161 G4double result = 0., position; 178 G4double result = 0., position;
162 G4int iTkin, iTransfer; 179 G4int iTkin, iTransfer;
163 180
164 for( iTkin = 0; iTkin < fEnergyBin; iTkin++) 181 for( iTkin = 0; iTkin < fEnergyBin; iTkin++)
165 { 182 {
166 if( Tkin < fEnergyVector->Energy(iTkin) ) << 183 if( Tkin < fEnergyVector->GetLowEdgeEnergy(iTkin) ) break;
167 } 184 }
168 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBi 185 if ( iTkin >= fEnergyBin ) iTkin = fEnergyBin-1; // Tkin is more then theMaxEnergy
169 if ( iTkin < 0 ) iTkin = 0; // aga 186 if ( iTkin < 0 ) iTkin = 0; // against negative index, Tkin < theMinEnergy
170 187
171 position = (*(*fAngleTable)(iTkin))(fAngle 188 position = (*(*fAngleTable)(iTkin))(fAngleBin-1)*G4UniformRand();
172 189
173 // G4cout<<"position = "<<position<<G4endl 190 // G4cout<<"position = "<<position<<G4endl;
174 191
175 for( iTransfer = 0; iTransfer < fAngleBin; 192 for( iTransfer = 0; iTransfer < fAngleBin; iTransfer++)
176 { 193 {
177 if( position <= (*(*fAngleTable)(iTkin)) 194 if( position <= (*(*fAngleTable)(iTkin))(iTransfer) ) break;
178 } 195 }
179 if (iTransfer >= fAngleBin-1) iTransfer = 196 if (iTransfer >= fAngleBin-1) iTransfer = fAngleBin-1;
180 197
181 // G4cout<<"iTransfer = "<<iTransfer<<G4en 198 // G4cout<<"iTransfer = "<<iTransfer<<G4endl;
182 199
183 result = GetTransfer(iTkin, iTransfer, pos 200 result = GetTransfer(iTkin, iTransfer, position);
184 201
185 // G4cout<<"t = "<<t<<G4endl; 202 // G4cout<<"t = "<<t<<G4endl;
186 203
187 204
188 return result; 205 return result;
189 } 206 }
190 207
191 ////////////////////////////////////////////// 208 /////////////////////////////////////////////////
192 209
193 G4double 210 G4double
194 G4NeutronElectronElModel:: GetTransfer( G4int 211 G4NeutronElectronElModel:: GetTransfer( G4int iTkin, G4int iTransfer, G4double position )
195 { 212 {
196 G4double x1, x2, y1, y2, randTransfer, delta 213 G4double x1, x2, y1, y2, randTransfer, delta, mean, epsilon = 1.e-6;
197 214
198 if( iTransfer == 0 || iTransfer == fAngleBi 215 if( iTransfer == 0 || iTransfer == fAngleBin-1 )
199 { 216 {
200 randTransfer = (*fAngleTable)(iTkin)->Ener << 217 randTransfer = (*fAngleTable)(iTkin)->GetLowEdgeEnergy(iTransfer);
>> 218 // iTransfer++;
201 } 219 }
202 else 220 else
203 { 221 {
204 if ( iTransfer >= G4int((*fAngleTable)(iTk 222 if ( iTransfer >= G4int((*fAngleTable)(iTkin)->GetVectorLength()) )
205 { 223 {
206 iTransfer = G4int((*fAngleTable)(iTkin)- << 224 iTransfer = (*fAngleTable)(iTkin)->GetVectorLength() - 1;
207 } 225 }
208 y1 = (*(*fAngleTable)(iTkin))(iTransfer-1) 226 y1 = (*(*fAngleTable)(iTkin))(iTransfer-1);
209 y2 = (*(*fAngleTable)(iTkin))(iTransfer); 227 y2 = (*(*fAngleTable)(iTkin))(iTransfer);
210 228
211 x1 = (*fAngleTable)(iTkin)->Energy(iTransf << 229 x1 = (*fAngleTable)(iTkin)->GetLowEdgeEnergy(iTransfer-1);
212 x2 = (*fAngleTable)(iTkin)->Energy(iTransf << 230 x2 = (*fAngleTable)(iTkin)->GetLowEdgeEnergy(iTransfer);
213 231
214 delta = y2 - y1; 232 delta = y2 - y1;
215 mean = y2 + y1; 233 mean = y2 + y1;
216 234
217 if ( x1 == x2 ) randTransfer = x2; 235 if ( x1 == x2 ) randTransfer = x2;
218 else 236 else
219 { 237 {
>> 238 // if ( y1 == y2 )
>> 239
220 if ( delta < epsilon*mean ) 240 if ( delta < epsilon*mean )
221 { 241 {
222 randTransfer = x1 + ( x2 - x1 )*G4Unif 242 randTransfer = x1 + ( x2 - x1 )*G4UniformRand();
223 } 243 }
224 else 244 else
225 { 245 {
226 randTransfer = x1 + ( position - y1 )* 246 randTransfer = x1 + ( position - y1 )*( x2 - x1 )/delta; // ( y2 - y1 );
227 } 247 }
228 } 248 }
229 } 249 }
230 return randTransfer; 250 return randTransfer;
231 } 251 }
232 252
233 ////////////////////////////////////////////// 253 //////////////////////////////////////////////////////////////
234 // 254 //
235 // Rosenbluth relation (ultra-relativistic!) i 255 // Rosenbluth relation (ultra-relativistic!) in the neutron rest frame,
236 // x = sin^2(theta/2), theta is the electron s 256 // x = sin^2(theta/2), theta is the electron scattering angle
237 // Magnetic form factor in the dipole approxim 257 // Magnetic form factor in the dipole approximation.
238 258
239 G4double G4NeutronElectronElModel::XscIntegran 259 G4double G4NeutronElectronElModel::XscIntegrand(G4double x)
240 { 260 {
241 G4double result = 1., q2, znq2, znf, znf2, z 261 G4double result = 1., q2, znq2, znf, znf2, znf4;
242 262
243 znq2 = 1. + 2.*fee*x/fM; 263 znq2 = 1. + 2.*fee*x/fM;
244 264
245 q2 = 4.*fee2*x/znq2; 265 q2 = 4.*fee2*x/znq2;
246 266
247 znf = 1 + q2/fMv2; 267 znf = 1 + q2/fMv2;
248 znf2 = znf*znf; 268 znf2 = znf*znf;
249 znf4 = znf2*znf2; 269 znf4 = znf2*znf2;
250 270
251 result /= ( x + fAm )*znq2*znq2*znf4; 271 result /= ( x + fAm )*znq2*znq2*znf4;
252 272
253 result *= ( 1 - x )/( 1 + q2/4./fM2 ) + 2.*x 273 result *= ( 1 - x )/( 1 + q2/4./fM2 ) + 2.*x;
254 274
255 return result; 275 return result;
256 } 276 }
257 277
258 ////////////////////////////////////////////// 278 ////////////////////////////////////////////////
259 // 279 //
260 // 280 //
261 281
262 G4HadFinalState* G4NeutronElectronElModel::App 282 G4HadFinalState* G4NeutronElectronElModel::ApplyYourself(
263 const G4HadProjectile& aTrack, G4Nucleus& << 283 const G4HadProjectile& aTrack, G4Nucleus& targetNucleus)
264 { 284 {
265 theParticleChange.Clear(); 285 theParticleChange.Clear();
266 286
267 const G4HadProjectile* aParticle = &aTrack; 287 const G4HadProjectile* aParticle = &aTrack;
268 G4double Tkin = aParticle->GetKineticEnergy( 288 G4double Tkin = aParticle->GetKineticEnergy();
269 fAm = CalculateAm( Tkin); 289 fAm = CalculateAm( Tkin);
270 // G4double En = aParticle->GetTotalEnergy 290 // G4double En = aParticle->GetTotalEnergy();
271 291
272 if( Tkin <= LowestEnergyLimit() ) 292 if( Tkin <= LowestEnergyLimit() )
273 { 293 {
274 theParticleChange.SetEnergyChange(Tkin); 294 theParticleChange.SetEnergyChange(Tkin);
275 theParticleChange.SetMomentumChange(aTrack 295 theParticleChange.SetMomentumChange(aTrack.Get4Momentum().vect().unit());
276 return &theParticleChange; 296 return &theParticleChange;
277 } 297 }
278 // sample e-scattering angle and make final 298 // sample e-scattering angle and make final state in lab frame
279 299
280 G4double sin2ht = SampleSin2HalfTheta( Tkin) 300 G4double sin2ht = SampleSin2HalfTheta( Tkin); // in n-rrest frame
281 301
282 // G4cout<<"sin2ht = "<<sin2ht<<G4endl; 302 // G4cout<<"sin2ht = "<<sin2ht<<G4endl;
283 303
284 G4double eTkin = fee; // fM; 304 G4double eTkin = fee; // fM;
285 305
286 eTkin /= 1.+2.*fee*sin2ht/fM; // fme/En + 2* 306 eTkin /= 1.+2.*fee*sin2ht/fM; // fme/En + 2*sin2ht;
287 307
288 eTkin -= fme; 308 eTkin -= fme;
289 309
290 // G4cout<<"eTkin = "<<eTkin<<G4endl; 310 // G4cout<<"eTkin = "<<eTkin<<G4endl;
291 311
292 if( eTkin > fCutEnergy ) 312 if( eTkin > fCutEnergy )
293 { 313 {
294 G4double ePlab = sqrt( eTkin*(eTkin + 2.*f 314 G4double ePlab = sqrt( eTkin*(eTkin + 2.*fme) );
295 315
296 // G4cout<<"ePlab = "<<ePlab<<G4endl; 316 // G4cout<<"ePlab = "<<ePlab<<G4endl;
297 317
298 G4double cost = 1. - 2*sin2ht; 318 G4double cost = 1. - 2*sin2ht;
299 319
300 if( cost > 1. ) cost = 1.; 320 if( cost > 1. ) cost = 1.;
301 if( cost < -1. ) cost = -1.; 321 if( cost < -1. ) cost = -1.;
302 322
303 G4double sint = std::sqrt( (1.0 - cost)*(1 323 G4double sint = std::sqrt( (1.0 - cost)*(1.0 + cost) );
304 G4double phi = G4UniformRand()*CLHEP::two 324 G4double phi = G4UniformRand()*CLHEP::twopi;
305 325
306 G4ThreeVector eP( sint*std::cos(phi), sint 326 G4ThreeVector eP( sint*std::cos(phi), sint*std::sin(phi), cost );
307 eP *= ePlab; 327 eP *= ePlab;
308 G4LorentzVector lvt2( eP, eTkin + electron 328 G4LorentzVector lvt2( eP, eTkin + electron_mass_c2 ); // recoil e- in n-rest frame
309 329
310 G4LorentzVector lvp1 = aParticle->Get4Mome 330 G4LorentzVector lvp1 = aParticle->Get4Momentum();
311 G4LorentzVector lvt1(0.,0.,0.,electron_mas 331 G4LorentzVector lvt1(0.,0.,0.,electron_mass_c2);
312 G4LorentzVector lvsum = lvp1+lvt1; 332 G4LorentzVector lvsum = lvp1+lvt1;
313 333
314 G4ThreeVector bst = lvp1.boostVector(); 334 G4ThreeVector bst = lvp1.boostVector();
315 lvt2.boost(bst); 335 lvt2.boost(bst);
316 336
317 // G4cout<<"lvt2 = "<<lvt2<<G4endl; 337 // G4cout<<"lvt2 = "<<lvt2<<G4endl;
318 338
319 G4DynamicParticle * aSec = new G4DynamicPa 339 G4DynamicParticle * aSec = new G4DynamicParticle( theElectron, lvt2 );
320 theParticleChange.AddSecondary( aSec, secI 340 theParticleChange.AddSecondary( aSec, secID );
321 341
322 G4LorentzVector lvp2 = lvsum-lvt2; 342 G4LorentzVector lvp2 = lvsum-lvt2;
323 343
324 // G4cout<<"lvp2 = "<<lvp2<<G4endl; 344 // G4cout<<"lvp2 = "<<lvp2<<G4endl;
325 345
326 G4double Tkin2 = lvp2.e()-aParticle->GetDe 346 G4double Tkin2 = lvp2.e()-aParticle->GetDefinition()->GetPDGMass();
327 theParticleChange.SetEnergyChange(Tkin2); 347 theParticleChange.SetEnergyChange(Tkin2);
328 theParticleChange.SetMomentumChange(lvp2.v 348 theParticleChange.SetMomentumChange(lvp2.vect().unit());
329 } 349 }
330 else if( eTkin > 0.0 ) 350 else if( eTkin > 0.0 )
331 { 351 {
332 theParticleChange.SetLocalEnergyDeposit( e 352 theParticleChange.SetLocalEnergyDeposit( eTkin );
333 Tkin -= eTkin; 353 Tkin -= eTkin;
334 354
335 if( Tkin > 0. ) 355 if( Tkin > 0. )
336 { 356 {
337 theParticleChange.SetEnergyChange( Tkin 357 theParticleChange.SetEnergyChange( Tkin );
338 theParticleChange.SetMomentumChange( aTr 358 theParticleChange.SetMomentumChange( aTrack.Get4Momentum().vect().unit() );
339 } 359 }
340 } 360 }
341 else 361 else
342 { 362 {
343 theParticleChange.SetEnergyChange( Tkin ); 363 theParticleChange.SetEnergyChange( Tkin );
344 theParticleChange.SetMomentumChange( aTrac 364 theParticleChange.SetMomentumChange( aTrack.Get4Momentum().vect().unit() );
345 } 365 }
>> 366 G4int Z = targetNucleus.GetZ_asInt();
>> 367 Z *= 1;
>> 368
346 return &theParticleChange; 369 return &theParticleChange;
347 } 370 }
348 371
349 // 372 //
350 // 373 //
351 /////////////////////////// 374 ///////////////////////////
352 375