| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 // Author: Alexei Sytov
27 // Co-author: Gianfranco Paternò (modificat
28 // On the base of the CRYSTALRAD realization o
29 // A. I. Sytov, V. V. Tikhomirov, and L. Bandi
30
31 /// \file G4ChannelingFastSimModel.cc
32 /// \brief Implementation of the G4ChannelingF
33 //
34 //
35 //
36 #include "G4ChannelingFastSimModel.hh"
37
38 #include "Randomize.hh"
39
40 #include "G4TransportationManager.hh"
41 #include "G4SystemOfUnits.hh"
42
43 //....oooOO0OOooo........oooOO0OOooo........oo
44
45 G4ChannelingFastSimModel::G4ChannelingFastSimM
46
47 : G4VFastSimulationModel(modelName, envelope)
48 {
49 }
50
51 //....oooOO0OOooo........oooOO0OOooo........oo
52
53 G4ChannelingFastSimModel::G4ChannelingFastSimM
54 : G4VFastSimulationModel(modelName)
55 {
56
57 }
58
59 //....oooOO0OOooo........oooOO0OOooo........oo
60
61 G4ChannelingFastSimModel::~G4ChannelingFastSim
62 {
63 }
64
65 //....oooOO0OOooo........oooOO0OOooo........oo
66
67 G4bool G4ChannelingFastSimModel::IsApplicable(
68 {
69 return std::abs(particleType.GetPDGCharge())
70 }
71
72 //....oooOO0OOooo........oooOO0OOooo........oo
73
74 G4bool G4ChannelingFastSimModel::ModelTrigger(
75 {
76 //default output
77 G4bool modelTrigger = false;
78
79 G4int particleDefinitionID =
80 fastTrack.GetPrimaryTrack()->GetPart
81 //kinetic energy
82 G4double ekinetic = fastTrack.GetPrimaryTrac
83
84 //energy cut, at the beginning, to not check
85 if(ekinetic > GetLowKineticEnergyLimit(parti
86 {
87 //current logical volume
88 G4LogicalVolume* crystallogic = fastTrac
89 fCrystalData->SetGeometryParameters(crys
90
91 G4ThreeVector momentumDirection = fastTr
92 // the particle angle vs crystal plane o
93 G4double angle = std::atan(momentumDirec
94 //recalculate angle into the lattice ref
95 angle = fCrystalData->
96 AngleXFromBoxToLattice(angle,
97 (fCrystal
98 fast
99 if (fCrystalData->GetModel()==2)
100 {
101 angle = std::sqrt(angle*angle+
102 std::pow(std::atan
103
104 }
105
106 //particle mass
107 G4double mass = fastTrack.GetPrimaryTrac
108 //particle total energy
109 G4double etotal = mass + ekinetic;
110 //particle charge
111 G4double charge = fastTrack.GetPrimaryTr
112 GetParticleDefinition(
113
114 //Particle position
115 G4ThreeVector xyz0 = fastTrack.GetPrimar
116 //Step estimate
117 G4double dz0 = fCrystalData->GetMaxSimul
118 xyz0 += 2*dz0*momentumDirection;//overes
119 //in cha
120
121 //Applies the parameterisation not at th
122 //above low energy limit and below angul
123
124 modelTrigger = (crystallogic->GetSolid()
125 Inside(xyz0)==kInside) &
126 momentumDirection.z()>0.
127 std::abs(angle) <
128 std::max(
129 GetLindhardAngleNumb
130 fCrystalData->GetLin
131
132
133 GetHighAngleLimit(pa
134 }
135
136 return modelTrigger;
137 }
138
139 //....oooOO0OOooo........oooOO0OOooo........oo
140
141 void G4ChannelingFastSimModel::DoIt(const G4Fa
142 G4FastStep& fastStep)
143 {
144 G4double etotal;//particle total energy
145 G4double etotalPreStep;//etotal at the previ
146 G4double etotalToSetParticleProperties;//eto
147 //Set
148 G4double ekinetic = 0;//kinetic energy
149 G4double eDeposited = 0.;//deposited energy
150 G4double elossAccum = 0;// accumulate local
151 G4double mass; //particle mass
152 G4double charge;//particle charge
153 G4double tGlobal; //global time
154 G4double tGlobalPreStep; //global time at th
155 G4ThreeVector xyz0;// the coordinates in the
156 G4ThreeVector xyz0PreStep;// xyz at the prev
157 G4ThreeVector xyz;// the coordinates in the
158 //a channel (elementary pe
159 G4double x,y,z; // the coordinates in the
160 //a channel (elementary pe
161 G4double tx0,ty0; // the angles in the local
162 G4double tx,ty; // the angles in the co-ro
163 //a channel (elementary pe
164 G4double txPreStep,tyPreStep;// tx,ty at the
165 G4ThreeVector momentumDirection;
166 G4ThreeVector scatteringAnglesAndEnergyLoss;
167 G4double lindhardAngleNumberHighLimit0; //cu
168 //[L
169 G4double highAngleLimit0; //current absolute
170
171 //coordinates in Runge-Kutta calculations
172 G4double x1=0.,x2=0.,x3=0.,x4=0.,y1=0.,y2=0.
173 //angles in Runge-Kutta calculations
174 G4double tx1=0.,tx2=0.,tx3=0.,tx4=0.,ty1=0.,
175 //variables in Runge-Kutta calculations
176 G4double kvx1=0.,kvx2=0.,kvx3=0.,kvx4=0.,kvy
177 //simulation step along z (internal step of
178 G4double dz,dzd3,dzd8;//dzd3 = dz/3; dzd8 =
179 //simulation step along the momentum directi
180 G4double momentumDirectionStep;
181 //effective simulation step (taking into acc
182 G4double effectiveStep;
183
184 // flag, if Inside(xyz0) switches to kInside
185 G4bool inside = false;
186
187 G4LogicalVolume* crystallogic = fastTrack.Ge
188 fCrystalData->SetGeometryParameters(crystall
189
190 //set the max number of secondaries (photons
191 if (fRad)
192 {
193 fastStep.SetNumberOfSecondaryTracks(fMax
194 //reseting the BaierKatkov integral to s
195 fBaierKatkov->ResetRadIntegral();//to av
196 }
197
198 mass = fastTrack.GetPrimaryTrack()->GetParti
199 etotal = mass + fastTrack.GetPrimaryTrack()-
200 charge = fastTrack.GetPrimaryTrack()->GetPar
201
202 G4String particleName =
203 fastTrack.GetPrimaryTrack()->GetParticle
204
205 lindhardAngleNumberHighLimit0 =
206 GetLindhardAngleNumberHighLimit(fastTrac
207 GetParti
208 highAngleLimit0 = GetHighAngleLimit(fastTrac
209 GetParti
210
211 //set fCrystalData parameters depending on t
212 fCrystalData->SetParticleProperties(etotal,
213
214 //global time
215 tGlobal = fastTrack.GetPrimaryTrack()->GetGl
216
217 //coordinates in the co-rotating reference s
218 xyz0= fastTrack.GetPrimaryTrackLocalPosition
219 xyz = fCrystalData->CoordinatesFromBoxToLatt
220 x=xyz.x();
221 y=xyz.y();
222 z=xyz.z();
223
224 momentumDirection=fastTrack.GetPrimaryTrackL
225 //angle in the co-rotating reference system
226 //(!!! ONLY FORWARD DIRECTION, momentumDirec
227 //valid for high energies defined by the sta
228 tx0 = std::atan(momentumDirection.x()/moment
229 ty0 = std::atan(momentumDirection.y()/moment
230
231 //angles in the co-rotating reference system
232 tx = fCrystalData->AngleXFromBoxToLattice(tx
233 ty = ty0;
234
235 etotalToSetParticleProperties = etotal*0.999
236 G4bool inCrystal=true;//flag necessary to es
237 //do calculations until the particle is insi
238 do
239 {
240 //remember the global time before the ne
241 tGlobalPreStep=tGlobal;
242 //remember the coordinates before the ne
243 xyz0PreStep = xyz0;
244 //remember the angles and the total ener
245 txPreStep = tx;
246 tyPreStep = ty;
247 etotalPreStep = etotal;
248
249 dz = fCrystalData->GetSimulationStep(tx,
250 dzd3=dz/3;
251 dzd8=dz/8;
252
253 //trajectory calculation:
254 //Runge-Cutt "3/8"
255 //fCrystalData->GetCurv(z)*fCrystalData-
256 //the radius on x; GetCurv gets 1/R for
257
258 //first step
259 kvx1=fCrystalData->Ex(x,y);
260 x1=x+tx*dzd3;
261 tx1=tx+(kvx1-fCrystalData->GetCurv(z)*fC
262 if (fCrystalData->GetModel()==2)
263 {
264 kvy1=fCrystalData->Ey(x,y);
265 y1=y+ty*dzd3;
266 ty1=ty+kvy1*dzd3;
267 }
268
269 //second step
270 kvx2=fCrystalData->Ex(x1,y1);
271 x2=x-tx*dzd3+tx1*dz;
272 tx2=tx-(kvx1-fCrystalData->GetCurv(z)*fC
273 (kvx2-fCrystalData->GetCurv(z)*f
274 if (fCrystalData->GetModel()==2)
275 {
276 kvy2=fCrystalData->Ey(x1,y1);
277 y2=y-ty*dzd3+ty1*dz;
278 ty2=ty-kvy1*dzd3+kvy2*dz;
279 }
280
281 //third step
282 kvx3=fCrystalData->Ex(x2,y2);
283 x3=x+(tx-tx1+tx2)*dz;
284 tx3=tx+(kvx1-kvx2+kvx3-fCrystalData->Get
285 if (fCrystalData->GetModel()==2)
286 {
287 kvy3=fCrystalData->Ey(x2,y2);
288 y3=y+(ty-ty1+ty2)*dz;
289 ty3=ty+(kvy1-kvy2+kvy3)*dz;
290 }
291
292 //fourth step
293 kvx4=fCrystalData->Ex(x3,y3);
294 x4=x+(tx+3.*tx1+3.*tx2+tx3)*dzd8;
295 tx4=tx+(kvx1+3.*kvx2+3.*kvx3+kvx4)*dzd8-
296 fCrystalData->GetCurv(z)*fCrysta
297 if (fCrystalData->GetModel()==2)
298 {
299 kvy4=fCrystalData->Ey(x3,y3);
300 y4=y+(ty+3.*ty1+3.*ty2+ty3)*dzd8;
301 ty4=ty+(kvy1+3.*kvy2+3.*kvy3+kvy4)*d
302 }
303 else
304 {
305 y4 =y+ty*dz;
306 ty4=ty;
307 }
308
309 x=x4;
310 tx=tx4;
311 y=y4;
312 ty=ty4;
313
314 z+=dz*fCrystalData->GetCorrectionZ();//m
315 //("
316
317 xyz = fCrystalData->ChannelChange(x,y,z)
318 x=xyz.x();
319 y=xyz.y();
320 z=xyz.z();
321
322 //the coordinates in the local reference
323 //this vector will be used in the cycle
324 //in the radiation model (if activated)
325 xyz0=fCrystalData->CoordinatesFromLattic
326
327 momentumDirectionStep=
328 dz*std::sqrt(1+std::pow(std::tan
329 tGlobal+=momentumDirectionStep/(fCrystal
330
331 //default scattering and energy loss 0
332 scatteringAnglesAndEnergyLoss = G4ThreeV
333
334 //calculate separately for each element
335 for (G4int i = 0; i < fCrystalData->GetN
336 {
337 //effective step taking into account
338 effectiveStep = momentumDirectionSte
339 //Coulomb scattering on screened ato
340 scatteringAnglesAndEnergyLoss += fCr
341 CoulombAtomicScatteri
342
343 //Amorphous part of ionization energ
344 elossAccum += fCrystalData->Ionizati
345 }
346 //electron scattering and coherent part
347 scatteringAnglesAndEnergyLoss += fCrysta
348
349
350
351 tx += scatteringAnglesAndEnergyLoss.x();
352 ty += scatteringAnglesAndEnergyLoss.y();
353 elossAccum += scatteringAnglesAndEnergyL
354
355 // recalculate the energy depended param
356 //(only if the energy decreased enough,
357 if (etotalToSetParticleProperties>etotal
358 {
359 fCrystalData->SetParticleProperties(
360 etotalToSetParticleProperties = etot
361 }
362
363 //chain of conditions to escape the cycl
364 // if Inside(xyz0)==kInside has been alr
365 //(a particle has been inside the crysta
366 if (inside)
367 {
368 // if low energy
369 if (etotal-mass<=GetLowKineticEnergyL
370
371
372 {inCrystal = false;}//escape the
373 //check if the angle w.r.t. the axes
374 //return to standard Geant4:
375 else if (fCrystalData->GetModel()==1)
376 {
377 //if the angle w.r.t. the planes
378 if (std::abs(tx) >=
379 std::max(lindhardAngleNumberH
380 fCrystalData->Ge
381 highAngleLimit0))
382 {inCrystal = false;}//escape t
383 }
384 else if (fCrystalData->GetModel()==2)
385 {
386 //if the angle w.r.t. the axes is
387 if (std::sqrt(tx*tx+ty*ty) >=
388 std::max(lindhardAngleNumberH
389 fCrystalData->Ge
390 highAngleLimit0))
391 {inCrystal = false;}//escape t
392 }
393
394 //radiation production & radiation
395 //works only if the radiation model
396 if (fRad)
397 {
398 //back to the local reference s
399 tx0 = fCrystalData->AngleXFromL
400 ty0 = ty;
401 //xyz0 was calculated above
402
403 //running the radiation model a
404 if(fBaierKatkov->DoRadiation(et
405 tx0
406 sca
407 sca
408 mom
409 cry
410 Get
411 Ins
412 // also it was checked if the p
413 // calculate the radiation inte
414 {
415 //a photon has been emitted
416 //shift the particle back i
417 etotal = fBaierKatkov->GetP
418 tx0 = fBaierKatkov->GetPart
419 ty0 = fBaierKatkov->GetPart
420 tGlobal = fBaierKatkov->Get
421 xyz0 = fBaierKatkov->GetPar
422
423 //add secondary photon
424 fBaierKatkov->GeneratePhoto
425
426 //particle energy was chang
427 fCrystalData->SetParticlePr
428
429 //coordinates in the co-rot
430 xyz = fCrystalData->Coordin
431 x=xyz.x();
432 y=xyz.y();
433 z=xyz.z();
434
435 //angles in the co-rotating
436 tx = fCrystalData->AngleXFr
437 ty = ty0;
438 }
439 }
440 else
441 {
442 //we calculate deposited energy
443 //of radiation otherwise we do
444 etotal -= elossAccum;
445 eDeposited += elossAccum;
446 elossAccum=0;
447 ekinetic = etotal-mass;
448 if(ekinetic<1*keV)
449 {
450 G4cout << "Warning in G4Cha
451 ekinetic << "<" << 1*keV <<
452 eDeposited-=(1*keV-ekinetic
453 ekinetic = 1*keV;
454 G4cout << "Setting deposite
455 eDeposited << " & ekinetic=
456 etotal = mass+ekinetic;
457 }
458 }
459
460 //precise check if the particle is
461 if (crystallogic->GetSolid()->
462 Inside(xyz0)!=kInside)
463 {
464 //one step back to remain insid
465 //after the escape of the volum
466 tGlobal = tGlobalPreStep;
467 xyz0 = xyz0PreStep;
468 tx = txPreStep;
469 ty = tyPreStep;
470 etotal = etotalPreStep;
471 z-=dz*fCrystalData->GetCorrecti
472 // change the flag => this part
473 // the model before escape this
474
475 inCrystal = false; //escape the
476 }
477 }
478 else
479 {
480 // if Inside(xyz0)==kInside we can e
481 if (crystallogic->GetSolid()->
482 Inside(xyz0)==kInsi
483 {inside = true;}
484 // a very rare case, if a particle r
485 // on the boundary and escapes the c
486 else if (crystallogic->GetSolid()->
487 Inside(xyz0)==kOuts
488 {inCrystal = false;}//escape
489 }
490 }
491 while (inCrystal);
492
493 //the angles in the local reference system o
494 tx0 = fCrystalData->AngleXFromLatticeToBox(t
495 ty0 = ty;
496
497 //set global time
498 fastStep.ProposePrimaryTrackFinalTime(tGloba
499 //set final position
500 fastStep.ProposePrimaryTrackFinalPosition(xy
501
502 //set deposited energy (due to ionization)
503 etotal -= elossAccum;
504 eDeposited += elossAccum;
505 ekinetic = etotal-mass;
506 if(ekinetic<1*keV)
507 {
508 G4cout << "Warning in G4ChannelingFastSi
509 ekinetic << "<" << 1*keV << " !" << G4en
510 eDeposited-=(1*keV-ekinetic);
511 ekinetic = 1*keV;
512 G4cout << "Setting deposited energy=" <<
513 eDeposited << " & ekinetic=" << ekinetic
514 }
515 fastStep.ProposeTotalEnergyDeposited(eDeposi
516 //set final kinetic energy
517 fastStep.ProposePrimaryTrackFinalKineticEner
518
519
520 //set final momentum direction
521 G4double momentumDirectionZ =
522 1./std::sqrt(1.+std::pow(std::tan(tx
523 fastStep.ProposePrimaryTrackFinalMomentumDir
524 G4ThreeVector(momentumDirectionZ
525 momentumDirectionZ
526 momentumDirectionZ
527 }
528
529 //....oooOO0OOooo........oooOO0OOooo........oo
530
531 void G4ChannelingFastSimModel::Input(const G4M
532 const G4S
533 const G4S
534 {
535 //initializing the class with containing al
536 //the crystal material and crystal lattice
537 //Channeling scattering and ionization proc
538 fCrystalData = new G4ChannelingFastSimCryst
539 //setting all the crystal material and latt
540 fCrystalData->SetMaterialProperties(crystal
541
542 //setting default low energy cuts for kinet
543 SetLowKineticEnergyLimit(1*GeV,"proton");
544 SetLowKineticEnergyLimit(1*GeV,"anti_proton
545 SetLowKineticEnergyLimit(200*MeV,"e-");
546 SetLowKineticEnergyLimit(200*MeV,"e+");
547
548 //set the model high limit of the angle exp
549 SetLindhardAngleNumberHighLimit(100.,"proto
550 SetLindhardAngleNumberHighLimit(100.,"anti_
551 SetLindhardAngleNumberHighLimit(100.,"e-");
552 SetLindhardAngleNumberHighLimit(100.,"e+");
553 }
554
555 //....oooOO0OOooo........oooOO0OOooo........oo
556
557 void G4ChannelingFastSimModel::RadiationModelA
558 {
559 fRad = true;
560 //activate the Baier-Katkov radiation mode
561 fBaierKatkov = new G4BaierKatkov();
562 }
563
564 //....oooOO0OOooo........oooOO0OOooo........oo
565