132{
135 G4double etotalToSetParticleProperties;
136
147
149
152
156 G4double lindhardAngleNumberHighLimit0;
157
159
160
161 G4double x1=0.,x2=0.,x3=0.,x4=0.,y1=0.,y2=0.,y3=0.,y4=0.;
162
163 G4double tx1=0.,tx2=0.,tx3=0.,tx4=0.,ty1=0.,ty2=0.,ty3=0.,ty4=0.;
164
165 G4double kvx1=0.,kvx2=0.,kvx3=0.,kvx4=0.,kvy1=0.,kvy2=0.,kvy3=0.,kvy4=0.;
166
168
170
172
173
175
177 fCrystalData->SetGeometryParameters(crystallogic);
178
179
180 if (fRad)
181 {
183
184 fBaierKatkov->ResetRadIntegral();
185 }
186
190
191 G4String particleName =
193
194 lindhardAngleNumberHighLimit0 =
196 GetParticleDefinition()->GetParticleDefinitionID());
198 GetParticleDefinition()->GetParticleDefinitionID());
199
200
201 fCrystalData->SetParticleProperties(etotal, mass, charge, particleName);
202
203
205
206
208 xyz = fCrystalData->CoordinatesFromBoxToLattice(xyz0);
212
214
215
216
217 tx0 = std::atan(momentumDirection.
x()/momentumDirection.
z());
218 ty0 = std::atan(momentumDirection.
y()/momentumDirection.
z());
219
220
221 tx = fCrystalData->AngleXFromBoxToLattice(tx0,z);
222 ty = ty0;
223
224 etotalToSetParticleProperties = etotal*0.999;
226
227 do
228 {
229
230 tGlobalPreStep=tGlobal;
231
232 xyz0PreStep = xyz0;
233
234 txPreStep = tx;
235 tyPreStep = ty;
236 etotalPreStep = etotal;
237
238 dz = fCrystalData->GetSimulationStep(tx,ty);
239 dzd3=dz/3;
240 dzd8=dz/8;
241
242
243
244
245
246
247
248 kvx1=fCrystalData->Ex(x,y);
249 x1=x+tx*dzd3;
250 tx1=tx+(kvx1-fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dzd3;
251 if (fCrystalData->GetModel()==2)
252 {
253 kvy1=fCrystalData->Ey(x,y);
254 y1=y+ty*dzd3;
255 ty1=ty+kvy1*dzd3;
256 }
257
258
259 kvx2=fCrystalData->Ex(x1,y1);
260 x2=x-tx*dzd3+tx1*dz;
261 tx2=tx-(kvx1-fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dzd3+
262 (kvx2-fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dz;
263 if (fCrystalData->GetModel()==2)
264 {
265 kvy2=fCrystalData->Ey(x1,y1);
266 y2=y-ty*dzd3+ty1*dz;
267 ty2=ty-kvy1*dzd3+kvy2*dz;
268 }
269
270
271 kvx3=fCrystalData->Ex(x2,y2);
272 x3=x+(tx-tx1+tx2)*dz;
273 tx3=tx+(kvx1-kvx2+kvx3-fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ())*dz;
274 if (fCrystalData->GetModel()==2)
275 {
276 kvy3=fCrystalData->Ey(x2,y2);
277 y3=y+(ty-ty1+ty2)*dz;
278 ty3=ty+(kvy1-kvy2+kvy3)*dz;
279 }
280
281
282 kvx4=fCrystalData->Ex(x3,y3);
283 x4=x+(tx+3.*tx1+3.*tx2+tx3)*dzd8;
284 tx4=tx+(kvx1+3.*kvx2+3.*kvx3+kvx4)*dzd8-
285 fCrystalData->GetCurv(z)*fCrystalData->GetCorrectionZ()*dz;
286 if (fCrystalData->GetModel()==2)
287 {
288 kvy4=fCrystalData->Ey(x3,y3);
289 y4=y+(ty+3.*ty1+3.*ty2+ty3)*dzd8;
290 ty4=ty+(kvy1+3.*kvy2+3.*kvy3+kvy4)*dzd8;
291 }
292 else
293 {
294 y4 =y+ty*dz;
295 ty4=ty;
296 }
297
298 x=x4;
299 tx=tx4;
300 y=y4;
301 ty=ty4;
302
303 z+=dz*fCrystalData->GetCorrectionZ();
304
305
306 xyz = fCrystalData->ChannelChange(x,y,z);
310
311
312
313
314 xyz0=fCrystalData->CoordinatesFromLatticeToBox(xyz);
315
316 momentumDirectionStep=
317 dz*std::sqrt(1+std::pow(std::tan(tx),2)+std::pow(std::tan(ty),2));
318 tGlobal+=momentumDirectionStep/(fCrystalData->GetBeta())/CLHEP::c_light;
319
320
322
323
324 for (
G4int i = 0; i < fCrystalData->GetNelements(); i++)
325 {
326
327 effectiveStep = momentumDirectionStep*fCrystalData->NuclearDensity(x,y,i);
328
329 scatteringAnglesAndEnergyLoss += fCrystalData->
330 CoulombAtomicScattering(effectiveStep,momentumDirectionStep,i);
331
332
333 elossAccum += fCrystalData->IonizationLosses(momentumDirectionStep, i);
334 }
335
336 scatteringAnglesAndEnergyLoss += fCrystalData->CoulombElectronScattering(
337 fCrystalData->MinIonizationEnergy(x,y),
338 fCrystalData->ElectronDensity(x,y),
339 momentumDirectionStep);
340 tx += scatteringAnglesAndEnergyLoss.
x();
341 ty += scatteringAnglesAndEnergyLoss.
y();
342 elossAccum += scatteringAnglesAndEnergyLoss.
z();
343
344
345
346 if (etotalToSetParticleProperties>etotal)
347 {
348 fCrystalData->SetParticleProperties(etotal, mass, charge, particleName);
349 etotalToSetParticleProperties = etotal*0.999;
350 }
351
352
353
354
355 if (inside)
356 {
357
359 GetParticleDefinition()->
360 GetParticleDefinitionID()))
361 {inCrystal = false;}
362
363
364 else if (fCrystalData->GetModel()==1)
365 {
366
367 if (std::abs(tx) >=
368 std::max(lindhardAngleNumberHighLimit0*
369 fCrystalData->GetLindhardAngle(),
370 highAngleLimit0))
371 {inCrystal = false;}
372 }
373 else if (fCrystalData->GetModel()==2)
374 {
375
376 if (std::sqrt(tx*tx+ty*ty) >=
377 std::max(lindhardAngleNumberHighLimit0*
378 fCrystalData->GetLindhardAngle(),
379 highAngleLimit0))
380 {inCrystal = false;}
381 }
382
383
384
385 if (fRad)
386 {
387
388 tx0 = fCrystalData->AngleXFromLatticeToBox(tx,z);
389 ty0 = ty;
390
391
392
393 if(fBaierKatkov->DoRadiation(etotal,mass,
394 tx0,ty0,
395 scatteringAnglesAndEnergyLoss.
x(),
396 scatteringAnglesAndEnergyLoss.
y(),
397 momentumDirectionStep,tGlobal,xyz0,
398 crystallogic->
399 GetSolid()->
400 Inside(xyz0)!=
kInside&&inCrystal))
401
402
403 {
404
405
406 etotal = fBaierKatkov->GetParticleNewTotalEnergy();
407 tx0 = fBaierKatkov->GetParticleNewAngleX();
408 ty0 = fBaierKatkov->GetParticleNewAngleY();
409 tGlobal = fBaierKatkov->GetNewGlobalTime();
410 xyz0 = fBaierKatkov->GetParticleNewCoordinateXYZ();
411
412
413 fBaierKatkov->GeneratePhoton(fastStep);
414
415
416 fCrystalData->SetParticleProperties(etotal, mass, charge, particleName);
417
418
419 xyz = fCrystalData->CoordinatesFromBoxToLattice(xyz0);
423
424
425 tx = fCrystalData->AngleXFromBoxToLattice(tx0,z);
426 ty = ty0;
427 }
428 }
429 else
430 {
431
432
433 etotal -= elossAccum;
434 eDeposited += elossAccum;
435 elossAccum=0;
436 ekinetic = etotal-mass;
437 if(ekinetic<1*CLHEP::keV)
438 {
439 G4cout <<
"Warning in G4ChannelingFastSimModel: " <<
440 ekinetic <<
"<" << 1*CLHEP::keV <<
" !" <<
G4endl;
441 eDeposited-=(1*CLHEP::keV-ekinetic);
442 ekinetic = 1*CLHEP::keV;
443 G4cout <<
"Setting deposited energy=" <<
444 eDeposited <<
" & ekinetic=" << ekinetic <<
G4endl;
445 etotal = mass+ekinetic;
446 }
447 }
448
449
452 {
453
454
455 tGlobal = tGlobalPreStep;
456 xyz0 = xyz0PreStep;
457 tx = txPreStep;
458 ty = tyPreStep;
459 etotal = etotalPreStep;
460 z-=dz*fCrystalData->GetCorrectionZ();
461
462
463
464 inCrystal = false;
465 }
466 }
467 else
468 {
469
472 {inside = true;}
473
474
477 {inCrystal = false;}
478 }
479 }
480 while (inCrystal);
481
482
483 tx0 = fCrystalData->AngleXFromLatticeToBox(tx,z);
484 ty0 = ty;
485
486
488
490
491
492 etotal -= elossAccum;
493 eDeposited += elossAccum;
494 ekinetic = etotal-mass;
495 if(ekinetic<1*CLHEP::keV)
496 {
497 G4cout <<
"Warning in G4ChannelingFastSimModel: " <<
498 ekinetic <<
"<" << 1*CLHEP::keV <<
" !" <<
G4endl;
499 eDeposited-=(1*CLHEP::keV-ekinetic);
500 ekinetic = 1*CLHEP::keV;
501 G4cout <<
"Setting deposited energy=" <<
502 eDeposited <<
" & ekinetic=" << ekinetic <<
G4endl;
503 }
505
507
508
509
511 1./std::sqrt(1.+std::pow(std::tan(tx0),2)+std::pow(std::tan(ty0),2));
514 momentumDirectionZ*std::tan(ty0),
515 momentumDirectionZ));
516}
CLHEP::Hep3Vector G4ThreeVector
G4GLOB_DLL std::ostream G4cout
G4double GetHighAngleLimit(G4int particleDefinitionID)
G4double GetLindhardAngleNumberHighLimit(G4int particleDefinitionID)
G4double GetLowKineticEnergyLimit(G4int particleDefinitionID)
get cuts
void SetNumberOfSecondaryTracks(G4int)
void ProposeTotalEnergyDeposited(G4double anEnergyPart)
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
G4double GetPDGCharge() const
const G4String & GetParticleName() const
const G4ParticleDefinition * GetParticleDefinition() const
G4double GetGlobalTime() const
G4double GetKineticEnergy() const