131{
134 G4double etotalToSetParticleProperties;
135
143
145
148
152 G4double lindhardAngleNumberHighLimit0;
153
154
155
156 G4double x1=0.,x2=0.,x3=0.,x4=0.,y1=0.,y2=0.,y3=0.,y4=0.;
157
158 G4double tx1=0.,tx2=0.,tx3=0.,tx4=0.,ty1=0.,ty2=0.,ty3=0.,ty4=0.;
159
160 G4double kvx1=0.,kvx2=0.,kvx3=0.,kvx4=0.,kvy1=0.,kvy2=0.,kvy3=0.,kvy4=0.;
161
163
165
167
168
170
173
174
175 if (fRad)
176 {
178
180 }
181
185
186
189
190 lindhardAngleNumberHighLimit0 =
192 GetParticleDefinition()->GetParticleDefinitionID());
193
194
196
197
199
200
206
208
209
210
211 tx0 = std::atan(momentumDirection.
x()/momentumDirection.
z());
212 ty0 = std::atan(momentumDirection.
y()/momentumDirection.
z());
213
214
216 ty = ty0;
217
218 etotalToSetParticleProperties = etotal*0.999;
220
221 do
222 {
223
224 tGlobalPreStep=tGlobal;
225
226 xyz0PreStep = xyz0;
227
228 txPreStep = tx;
229 tyPreStep = ty;
230 etotalPreStep = etotal;
231
233 dzd3=dz/3;
234 dzd8=dz/8;
235
236
237
238
239
240
241
242 kvx1=fCrystalData->
Ex(x,y);
243 x1=x+tx*dzd3;
246 {
247 kvy1=fCrystalData->
Ey(x,y);
248 y1=y+ty*dzd3;
249 ty1=ty+kvy1*dzd3;
250 }
251
252
253 kvx2=fCrystalData->
Ex(x1,y1);
254 x2=x-tx*dzd3+tx1*dz;
258 {
259 kvy2=fCrystalData->
Ey(x1,y1);
260 y2=y-ty*dzd3+ty1*dz;
261 ty2=ty-kvy1*dzd3+kvy2*dz;
262 }
263
264
265 kvx3=fCrystalData->
Ex(x2,y2);
266 x3=x+(tx-tx1+tx2)*dz;
269 {
270 kvy3=fCrystalData->
Ey(x2,y2);
271 y3=y+(ty-ty1+ty2)*dz;
272 ty3=ty+(kvy1-kvy2+kvy3)*dz;
273 }
274
275
276 kvx4=fCrystalData->
Ex(x3,y3);
277 x4=x+(tx+3.*tx1+3.*tx2+tx3)*dzd8;
278 tx4=tx+(kvx1+3.*kvx2+3.*kvx3+kvx4)*dzd8-
281 {
282 kvy4=fCrystalData->
Ey(x3,y3);
283 y4=y+(ty+3.*ty1+3.*ty2+ty3)*dzd8;
284 ty4=ty+(kvy1+3.*kvy2+3.*kvy3+kvy4)*dzd8;
285 }
286 else
287 {
288 y4 =y+ty*dz;
289 ty4=ty;
290 }
291
292 x=x4;
293 tx=tx4;
294 y=y4;
295 ty=ty4;
296
298
299
304
305
306
307
309
310 momentumDirectionStep=
311 dz*std::sqrt(1+std::pow(std::tan(tx),2)+std::pow(std::tan(ty),2));
312 tGlobal+=momentumDirectionStep/(fCrystalData->
GetBeta())/CLHEP::c_light;
313
314
316
317
319 {
320
321 effectiveStep = momentumDirectionStep*fCrystalData->
NuclearDensity(x,y,i);
322
323 scatteringAnglesAndEnergyLoss += fCrystalData->
324 CoulombAtomicScattering(effectiveStep,momentumDirectionStep,i);
325
326
328 }
329
333 momentumDirectionStep);
334 tx += scatteringAnglesAndEnergyLoss.
x();
335 ty += scatteringAnglesAndEnergyLoss.
y();
336 etotal -= scatteringAnglesAndEnergyLoss.
z();
337
338
339
340 if (etotalToSetParticleProperties>etotal)
341 {
343 etotalToSetParticleProperties = etotal*0.999;
344 }
345
346
347
348
349 if (inside)
350 {
351
353 GetParticleDefinition()->
354 GetParticleDefinitionID()))
355 {inCrystal = false;}
356
357
358 else if (fCrystalData->
GetModel()==1)
359 {
360
361 if (std::abs(tx) >=
363 {inCrystal = false;}
364 }
365 else if (fCrystalData->
GetModel()==2)
366 {
367
368 if (std::sqrt(tx*tx+ty*ty) >= lindhardAngleNumberHighLimit0*
370 {inCrystal = false;}
371 }
372
373
374
375 if (fRad)
376 {
377
379 ty0 = ty;
380
381
382
384 tx0,ty0,
385 scatteringAnglesAndEnergyLoss.
x(),
386 scatteringAnglesAndEnergyLoss.
y(),
387 momentumDirectionStep,tGlobal,xyz0,
388 crystallogic->
389 GetSolid()->
390 Inside(xyz0)!=
kInside&&inCrystal))
391
392
393 {
394
395
401
402
404
405
407
408
413
414
416 ty = ty0;
417 }
418 }
419
420
423 {
424
425
426 tGlobal = tGlobalPreStep;
427 xyz0 = xyz0PreStep;
428 tx = txPreStep;
429 ty = tyPreStep;
430 etotal = etotalPreStep;
432
433
434
435 inCrystal = false;
436 }
437 }
438 else
439 {
440
443 {inside = true;}
444
445
448 {inCrystal = false;}
449 }
450 }
451 while (inCrystal);
452
453
455 ty0 = ty;
456
457
459
461
464 GetParticleDefinition()->GetPDGMass());
465
467 1./std::sqrt(1.+std::pow(std::tan(tx0),2)+std::pow(std::tan(ty0),2));
470 momentumDirectionZ*std::tan(ty0),
471 momentumDirectionZ));
472}
CLHEP::Hep3Vector G4ThreeVector
void GeneratePhoton(G4FastStep &fastStep)
G4double GetParticleNewAngleX()
G4double GetParticleNewTotalEnergy()
G4bool DoRadiation(G4double etotal, G4double mass, G4double angleX, G4double angleY, G4double angleScatteringX, G4double angleScatteringY, G4double step, G4double globalTime, G4ThreeVector coordinateXYZ, G4bool flagEndTrajectory=false)
G4double GetParticleNewAngleY()
G4double GetNewGlobalTime()
const G4ThreeVector & GetParticleNewCoordinateXYZ()
G4ThreeVector CoordinatesFromBoxToLattice(const G4ThreeVector &pos0)
G4double AngleXFromBoxToLattice(G4double tx, G4double z)
G4ThreeVector ChannelChange(G4double &x, G4double &y, G4double &z)
change the channel if necessary, recalculate x o y
G4ThreeVector CoordinatesFromLatticeToBox(const G4ThreeVector &pos)
G4double AngleXFromLatticeToBox(G4double tx, G4double z)
G4double GetLowKineticEnergyLimit(const G4String &particleName)
get cuts
G4double GetLindhardAngleNumberHighLimit(const G4String &particleName)
void SetNumberOfSecondaryTracks(G4int)
void ProposePrimaryTrackFinalPosition(const G4ThreeVector &, G4bool localCoordinates=true)
void ProposePrimaryTrackFinalTime(G4double)
void ProposePrimaryTrackFinalKineticEnergy(G4double)
void ProposePrimaryTrackFinalMomentumDirection(const G4ThreeVector &, G4bool localCoordinates=true)
G4ThreeVector GetPrimaryTrackLocalPosition() const
const G4Track * GetPrimaryTrack() const
G4ThreeVector GetPrimaryTrackLocalDirection() const
G4LogicalVolume * GetEnvelopeLogicalVolume() const
G4VSolid * GetSolid() const
G4double GetPDGMass() const
G4int GetLeptonNumber() const
G4double GetPDGCharge() const
const G4ParticleDefinition * GetParticleDefinition() const
G4double GetGlobalTime() const
G4double GetTotalEnergy() const
G4double IonizationLosses(G4double dz, G4int ielement)
ionization losses
G4double GetLindhardAngle(G4double etotal, G4double mass)
Calculate the value of the Lindhard angle (!!! the value for a straight crystal)
G4double GetBeta()
get particle velocity/c
void SetGeometryParameters(const G4LogicalVolume *crystallogic)
set geometry parameters from current logical volume
G4double ElectronDensity(G4double x, G4double y)
electron density function
G4double GetSimulationStep(G4double tx, G4double ty)
G4double MinIonizationEnergy(G4double x, G4double y)
minimum energy of ionization function
G4double Ex(G4double x, G4double y)
electric fields produced by crystal lattice
G4ThreeVector CoulombElectronScattering(G4double eMinIonization, G4double electronDensity, G4double step)
multiple and single scattering on electrons
void SetParticleProperties(G4double etotal, G4double mp, G4double charge, G4bool ifhadron)
G4double NuclearDensity(G4double x, G4double y, G4int ielement)
nuclear density function (normalized to average nuclear density)
G4double Ey(G4double x, G4double y)
G4double GetCorrectionZ()
G4double GetCurv(G4double z)