| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 // Author: Alexei Sytov
27 // Co-author: Gianfranco Paternò (modifications & testing)
28 // On the base of the CRYSTALRAD realization of scattering model:
29 // A. I. Sytov, V. V. Tikhomirov, and L. Bandiera PRAB 22, 064601 (2019)
30
31 #include "G4VChannelingFastSimCrystalData.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4PhysicalConstants.hh"
34 #include "G4Log.hh"
35
36 G4VChannelingFastSimCrystalData::G4VChannelingFastSimCrystalData()
37 {
38
39 }
40
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42
43 G4VChannelingFastSimCrystalData::~G4VChannelingFastSimCrystalData(){;}
44
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46
47 void G4VChannelingFastSimCrystalData::SetGeometryParameters
48 (const G4LogicalVolume *crystallogic)
49 {
50 G4int crystalID = crystallogic->GetInstanceID();
51
52 //set bending angle if the it exists in the list, otherwise default = 0
53 (fMapBendingAngle.count(crystalID) > 0)
54 ? SetBendingAngle(fMapBendingAngle[crystalID],crystallogic)
55 : SetBendingAngle(0.,crystallogic);
56
57 //set miscut angle if the it exists in the list, otherwise default = 0
58 (fMapMiscutAngle.count(crystalID) > 0)
59 ? SetMiscutAngle(fMapMiscutAngle[crystalID],crystallogic)
60 : SetMiscutAngle(0.,crystallogic);
61
62 //set crystalline undulator parameters if they exist in the list,
63 //otherwise default = G4ThreeVector(0,0,0).
64 (fMapCUAmplitudePeriodPhase.count(crystalID) > 0)
65 ? SetCUParameters(fMapCUAmplitudePeriodPhase[crystalID],crystallogic)
66 : SetCUParameters(G4ThreeVector(0.,0.,0.),crystallogic);
67 }
68
69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70
71 void G4VChannelingFastSimCrystalData::SetBendingAngle(G4double tetab,
72 const G4LogicalVolume* crystallogic)
73 {
74 G4int crystalID = crystallogic->GetInstanceID();
75
76 //set the bending angle for this logical volume
77 fMapBendingAngle[crystalID]=tetab;
78
79 G4ThreeVector limboxmin;//minimal limits of the box bounding the logical volume
80 G4ThreeVector limboxmax;//maximal limits of the box bounding the logical volume
81 //save the limits of the box bounding the logical volume
82 crystallogic->GetSolid()->BoundingLimits(limboxmin,limboxmax);
83
84 //bounding box half dimensions
85 fHalfDimBoundingBox = (limboxmax-limboxmin)/2.;
86
87 G4double lcr = limboxmax.getZ()-limboxmin.getZ();//crystal thickness
88
89 fBendingAngle=std::abs(tetab);
90 if (fBendingAngle<0.000001)//no bending less then 1 urad
91 {
92 if(fBendingAngle>DBL_EPSILON)
93 {
94 G4cout << "Channeling model: volume " << crystallogic->GetName() << G4endl;
95 G4cout << "Warning: bending angle is lower than 1 urad => set to 0" << G4endl;
96 }
97
98 fBent=0;
99 fBendingAngle=0.;
100 fBendingR=0.;//just for convenience (infinity in reality)
101 fBending2R=0.;
102 fBendingRsquare=0.;
103 fCurv=0.;
104
105 fCorrectionZ = 1.;
106 }
107 else
108 {
109 fBent=1;
110 fBendingR=lcr/fBendingAngle;
111 fBending2R=2.*fBendingR;
112 fBendingRsquare=fBendingR*fBendingR;
113 fCurv=1./fBendingR;
114
115 if (tetab<0.)
116 {
117 G4cout << "Channeling model: volume " << crystallogic->GetName() << G4endl;
118 G4cout << "Warning: bending angle is negative => set to be positive" << G4endl;
119 }
120 }
121
122 }
123
124 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
125
126 void G4VChannelingFastSimCrystalData::SetMiscutAngle(G4double tetam,
127 const G4LogicalVolume *crystallogic)
128 {
129 G4int crystalID = crystallogic->GetInstanceID();
130
131 //set the bending angle for this logical volume
132 fMapMiscutAngle[crystalID]=tetam;
133
134 // fMiscutAngle>0: rotation of xz coordinate planes clockwise in the xz plane
135 fMiscutAngle=tetam;
136 if (std::abs(tetam)>1.*mrad)
137 {
138 G4cout << "Channeling model: volume " << crystallogic->GetName() << G4endl;
139 G4cout << "Warning: miscut angle is higher than 1 mrad => " << G4endl;
140 G4cout << "coordinate transformation routines may be unstable" << G4endl;
141 }
142 fCosMiscutAngle=std::cos(fMiscutAngle);
143 fSinMiscutAngle=std::sin(fMiscutAngle);
144 }
145
146 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147
148 void G4VChannelingFastSimCrystalData::SetCrystallineUndulatorParameters(
149 G4double amplitude,
150 G4double period,
151 G4double phase,
152 const G4LogicalVolume *crystallogic)
153 {
154 if (amplitude<DBL_EPSILON||period<DBL_EPSILON)
155 {
156 amplitude = 0.;
157 period=0.;
158 phase=0.;
159 G4cout << "Channeling model: volume " << crystallogic->GetName() << G4endl;
160 G4cout << "Warning: The crystalline undulator parameters are out of range "
161 "=> the crystalline undulator mode switched off" << G4endl;
162 }
163
164 SetCUParameters(G4ThreeVector(amplitude,period,phase),crystallogic);
165 }
166
167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168
169 void G4VChannelingFastSimCrystalData::SetCUParameters(
170 const G4ThreeVector &litudePeriodPhase,
171 const G4LogicalVolume *crystallogic)
172 {
173 G4int crystalID = crystallogic->GetInstanceID();
174
175 //set the crystalline undulator parameters for this logical volume
176 fMapCUAmplitudePeriodPhase[crystalID]=amplitudePeriodPhase;
177 fCUAmplitude=amplitudePeriodPhase.x();
178 G4double period = amplitudePeriodPhase.y();
179 fCUPhase = amplitudePeriodPhase.z();
180
181 //if the amplidude of the crystalline undulator is 0 => no undulator
182 if(fCUAmplitude>DBL_EPSILON&&period>DBL_EPSILON)
183 {
184 //crystalline undulator flag
185 fCU = true;
186
187 fCUK = CLHEP::twopi/period;
188
189 if(fBendingAngle>DBL_EPSILON)
190 {
191 //bent and periodically bent crystal are not compatible
192 SetBendingAngle(0,crystallogic);
193
194 G4cout << "Channeling model: volume " << crystallogic->GetName() << G4endl;
195 G4cout << "Warning: crystalline undulator is not compatible with "
196 "a bent crystal mode => setting bending angle to 0." << G4endl;
197 }
198 }
199 else
200 {
201 fCU = false;
202 fCUAmplitude = 0.;
203 fCUK = 0.;
204 fCUPhase = 0.;
205 fMapCUAmplitudePeriodPhase[crystalID] = G4ThreeVector(0.,0.,0.);
206 }
207
208 fCUK2 = fCUK*fCUK;
209 fCUAmplitudeK = fCUAmplitude*fCUK;
210 }
211
212 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
213
214 void G4VChannelingFastSimCrystalData::SetParticleProperties(G4double etotal,
215 G4double mass,
216 G4double charge,
217 const G4String& particleName)
218 {
219 G4double teta1;
220 fZ2=charge;
221 G4double zz22=fZ2*fZ2;
222 fParticleName=particleName;
223
224 // particle momentum and energy
225 G4double t=etotal*etotal-mass*mass; // economy of operations
226 fPz=std::sqrt(t); // momentum of particle
227 fPV=t/etotal; // pv
228 fBeta=fPz/etotal; // velocity/c
229 fTetaL = std::sqrt(std::abs(fZ2)*fVmax2/fPV); //Lindhard angle
230 fChannelingStep = fChangeStep/fTetaL; //standard simulation step
231
232 // Energy losses
233 fV2 = fBeta*fBeta; // particle (velocity/c)^2
234 fGamma = etotal/mass; // Lorentz factor
235 fMe2Gamma = 2*CLHEP::electron_mass_c2*fGamma;
236 // max ionization losses
237 fTmax = fMe2Gamma*fGamma*fV2/
238 (CLHEP::electron_mass_c2/mass*CLHEP::electron_mass_c2/mass +
239 1. + fMe2Gamma/mass);
240 // max ionization losses for electrons
241 if(fParticleName=="e-"){fTmax/=2;}
242
243 for(G4int i=0; i<fNelements; i++)
244 {
245
246 // minimal scattering angle by coulomb scattering on nuclei
247 // defining by shielding by electrons
248 // teta1=hdc/(fPz*fRF)*DSQRT(1.13D0+3.76D0*(alpha*fZ1*fZ2/fBeta)**2){ev*cm/(eV*cm)}
249 teta1=fTeta10[i]*std::sqrt(1.13+fK40[i]*zz22/fV2); // /fPz later to speed up
250 // the calculations
251
252 // the coefficient for multiple scattering
253 fBB[i]=teta1*teta1*fPu11[i];
254 fE1XBbb[i]=expint(fBB[i]);
255 fBBDEXP[i]=(1.+fBB[i])*std::exp(fBB[i]);
256 // necessary for suppression of incoherent scattering
257 // by the atomic correlations in crystals for single scattering on nucleus
258 // (screened atomic potential): EXP(-(fPz*teta*fU1)**2)=EXP(-fPzu11*teta**2).GE.ksi
259 // =>no scattering
260 fPzu11[i]=fPu11[i]*fPz*fPz;
261
262 teta1=teta1/fPz; //
263 fTeta12[i]=teta1*teta1;
264 // maximal scattering angle by coulomb scattering on nuclei
265 // defining by nucleus radius
266 // tetamax=hc/(fPz*1.D-6*fR0*fAN**(1.D0/3.D0))// {Mev*fermi/(MeV*fermi)}
267 G4double tetamax=fTetamax0[i]/fPz;
268 fTetamax2[i]=tetamax*tetamax;
269 fTetamax12[i]=fTeta12[i]+fTetamax2[i];
270
271 // a coefficient in a formula for scattering (for high speed of simulation)
272 // fK2=(fZ2)**2*alphahbarc2*4.*pi*fN0*(fZ1/fPV)**2
273 fK2[i]=fK20[i]*zz22/fPV/fPV;
274 }
275
276 // fK3=(fZ2)**2*alphahbarc2*pi/electron_mass_c2/(fV2)**2
277 fK3=fK30*zz22/fV2;
278 }
279
280 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
281
282 G4double G4VChannelingFastSimCrystalData::GetLindhardAngle(G4double etotal,
283 G4double mass,
284 G4double charge)
285 {
286 G4double pv0 = etotal-mass*mass/etotal;
287 return std::sqrt(2*std::abs(charge)*fVmax/pv0); //Calculate the value of the Lindhard angle
288 //(!!! the value for a straight crystal)
289 }
290
291 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
292
293 G4double G4VChannelingFastSimCrystalData::GetLindhardAngle()
294 {
295 return fTetaL; //return the Lindhard angle value calculated in SetParticleProperties
296 }
297
298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
299
300 G4double G4VChannelingFastSimCrystalData::GetSimulationStep(G4double tx,G4double ty)
301 {
302 G4double simulationstep;
303 //find angle of particle w.r.t. the plane or axis
304 G4double angle=0.;
305 if (iModel==1)//1D model
306 {
307 angle = std::abs(tx);
308 }
309 else if (iModel==2)//2D model
310 {
311 angle = std::sqrt(tx*tx+ty*ty);
312 }
313
314 //compare this angle with the Lindhard angle
315 if (angle<fTetaL)
316 {
317 simulationstep = fChannelingStep;
318 }
319 else
320 {
321 simulationstep = fChangeStep;
322 if (angle > 0.0) { simulationstep /= angle; }
323 }
324
325 return simulationstep;
326 }
327
328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
329
330 G4double G4VChannelingFastSimCrystalData::GetMaxSimulationStep(G4double etotal,
331 G4double mass,
332 G4double charge)
333 {
334 //standard value of step for channeling particles which is the maximal possible step
335 return fChangeStep/GetLindhardAngle(etotal, mass, charge);
336 }
337
338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
339
340 G4ThreeVector G4VChannelingFastSimCrystalData::CoulombAtomicScattering(
341 G4double effectiveStep,
342 G4double step,
343 G4int ielement)
344 {
345 G4double tx = 0.;//horizontal scattering angle
346 G4double ty = 0.;//vertical scattering angle
347
348 G4double ksi=0.1;
349
350 // calculation of the teta2-minimal possible angle of a single scattering
351 G4double e1=fK2[ielement]*effectiveStep; //for high speed of a program
352 // (real formula is (4*pi*fN0*wpl(x)*dz*fZ1*zz2*alpha*hdc/fPV)**2)
353 G4double teta122=fTetamax12[ielement]/(ksi*fTetamax12[ielement]/e1+1.);
354 // teta122=fTeta12+teta22=teta1^2+teta2^2
355
356 G4double teta22;
357 G4double t;
358 // if the angle of a single scattering is less teta1 - minimal possible
359 // angle of coulomb scattering defining by the electron shielding than
360 // multiple scattering by both nuclei and electrons and electrons will not
361 // occur => minimal possible angle of a single scattering is equal to teta1
362 if (teta122<=fTeta12[ielement]*1.000125)
363 {
364 teta22=0.;
365 teta122=fTeta12[ielement];
366 }
367 else
368 {
369 teta22=teta122-fTeta12[ielement];
370 G4double aa=teta22/fTeta12[ielement];
371 G4double aa1=1.+aa;
372
373 // crystal, with scattering suppression
374 G4double tetamsi=e1*(G4Log(aa1)+
375 (1.-std::exp(-aa*fBB[ielement]))/aa1+
376 fBBDEXP[ielement]*
377 (expint(fBB[ielement]*aa1)-fE1XBbb[ielement]));
378
379 // sumilation of multiple coulomb scattering by nuclei and electrons
380 // for high speed of a program, real formula is
381 // 4*pi*fN0*wpl(x)*dz*(fZ1*zz2*alpha*hdc/fPV)**2*
382 // *(ln(1+a)+(1-exp(-a*b))/(1+a)+(1+b)*exp(b)*(E1XB(b*(1+a))-E1XB(b)))
383
384 ksi=G4UniformRand();
385 t=std::sqrt(-tetamsi*G4Log(ksi));
386
387 ksi=G4UniformRand();
388
389 tx+=t*std::cos(CLHEP::twopi*ksi);
390 ty+=t*std::sin(CLHEP::twopi*ksi);
391
392 }
393 // simulation of single coulomb scattering by nuclei (with screened potential)
394 G4double zss=0.;
395 G4double dzss=step;
396
397 // (calculation of a distance, at which another single scattering can happen)
398 ksi=G4UniformRand();
399
400 zss=-G4Log(ksi)*step/(e1*(1./teta122-1./fTetamax12[ielement]));
401 G4double tt;
402
403 // At some step several single scattering can occur.
404 // So,if the distance of the next scattering is less than the step,
405 // another scattering can occur. If the distance of the next scattering
406 // is less than the difference between the step and the distance of
407 // the previous scattering, another scattering can occur. And so on, and so on.
408 // In the cycle we simulate each of them. The cycle is finished, when
409 // the remaining part of step is less than a distance of the next single scattering.
410 //********************************************
411 // if at a step a single scattering occurs
412 while (zss<dzss)
413 {
414
415 // simulation by Monte-Carlo of angles of single scattering
416 ksi=G4UniformRand();
417
418 tt=fTetamax12[ielement]/(1.+ksi*(fTetamax2[ielement]-teta22)/teta122)-
419 fTeta12[ielement];
420
421 ksi=G4UniformRand();
422
423 // suppression of incoherent scattering by the atomic correlations in crystals
424 t=fPzu11[ielement]*tt;
425 t=std::exp(-t);
426
427 if (t<ksi) //if scattering takes place
428 {
429 //scattering angle
430 t=std::sqrt(tt);
431 ksi=G4UniformRand();
432
433 tx+=t*std::cos(CLHEP::twopi*ksi);
434 ty+=t*std::sin(CLHEP::twopi*ksi);
435 }
436
437 dzss-=zss;
438 // (calculation of a distance, at which another single scattering can happen)
439 ksi=G4UniformRand();
440
441 zss=-G4Log(ksi)*step/(e1*(1./teta122-1./fTetamax12[ielement]));
442 }
443 //********************************************
444 return G4ThreeVector(tx,ty,0.);
445 }
446
447 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
448
449 G4ThreeVector G4VChannelingFastSimCrystalData::CoulombElectronScattering(
450 G4double eMinIonization,
451 G4double electronDensity,
452 G4double step)
453 {
454
455 G4double zss=0.;
456 G4double dzss=step;
457 G4double ksi = 0.;
458
459 G4double tx = 0.;//horizontal scattering angle
460 G4double ty = 0.;//vertical scattering angle
461 G4double eloss = 0.;//energy loss
462
463 // eMinIonization - minimal energy transfered to electron
464 // a cut to reduce the number of calls of electron scattering
465 // is needed only at low density regions, in many cases does not do anything at all
466 if (eMinIonization<0.5*eV){eMinIonization=0.5*eV;}
467
468 // single scattering on electrons routine
469 if ((eMinIonization<fTmax)&&(electronDensity>DBL_EPSILON))
470 {
471
472 // (calculation of a distance, at which another single scattering can happen)
473 // simulation of scattering length (by the same way single scattering by nucleus
474 ksi=G4UniformRand();
475
476 zss=-1.0*G4Log(ksi)/(fK3*electronDensity)/(1./eMinIonization-1./fTmax);
477
478 //********************************************
479 // if at a step a single scattering occur
480 while (zss<dzss)
481 {
482 // simulation by Monte-Carlo of angles of single scattering
483 ksi=G4UniformRand();
484
485 // energy transfered to electron
486 G4double e1=eMinIonization/(1.-ksi*(1.-eMinIonization/fTmax));
487
488 // scattering angle
489 G4double t=0;
490 if(fTmax-e1>DBL_EPSILON) //to be sure e1<fTmax
491 {
492 t=std::sqrt(2.*CLHEP::electron_mass_c2*e1*(1-e1/fTmax))/fPz;
493 }
494
495 // energy losses
496 eloss=e1;
497 ksi=G4UniformRand();
498
499 tx+=t*std::cos(CLHEP::twopi*ksi);
500 ty+=t*std::sin(CLHEP::twopi*ksi);
501
502 dzss-=zss;
503 // (calculation of a distance, at which another single scattering can happen)
504 // simulation of scattering length
505 // (by the same way single scattering by nucleus
506 ksi=G4UniformRand();
507
508 zss=-1.0*G4Log(ksi)/(fK3*electronDensity)/(1./eMinIonization-1./fTmax);
509 }
510 //********************************************
511 }
512 return G4ThreeVector(tx,ty,eloss);
513 }
514
515 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
516
517 G4double G4VChannelingFastSimCrystalData::IonizationLosses(G4double dz,
518 G4int ielement)
519 {
520 //amorphous part of ionization losses
521
522 G4double elosses = 0.;
523 // 1/2 already taken into account in fKD
524
525 G4double loge = G4Log(fMe2Gamma*fGamma*fV2/fI0[ielement]);
526 G4double delta= 2*(G4Log(fBeta*fGamma)+fLogPlasmaEdI0[ielement]-0.5);
527 if(delta<0){delta=0;}
528 loge-=delta;
529 if(fParticleName=="e-")
530 {
531 loge+=(-G4Log(2.) + 1
532 -(2*fGamma - 1)/fGamma/fGamma*G4Log(2.) +
533 1/8*((fGamma - 1)/fGamma)*((fGamma - 1)/fGamma));
534 }
535 else if(fParticleName=="e+")
536 {
537 loge+=(-fV2/12*(11 + 14/(fGamma + 1) + 10/(fGamma + 1)/(fGamma + 1) +
538 4/(fGamma + 1)/(fGamma + 1)/(fGamma + 1)));
539 }
540 else
541 {
542 loge-=fV2;
543 }
544 elosses=fZ2*fZ2*fKD[ielement]/fV2*loge*dz;
545
546 return elosses;}
547
548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
549
550 G4double G4VChannelingFastSimCrystalData::expint(G4double X)
551 {
552 // ============================================
553 // Purpose: Compute exponential integral E1(x)
554 // Input : x --- Argument of E1(x)
555 // Output: E1 --- E1(x)
556 // ============================================
557
558 G4double E1, R, T, T0;
559 G4int M;
560
561 if (X==0)
562 {
563 E1=1.e300;
564 }
565 else if (X<=1.)
566 {
567 E1=1.;
568 R=1.;
569
570
571 for(int K=1; K<=25; K++)
572 {
573 R=-R*K*X/std::pow(K+1.,2.);
574 E1=E1+R;
575 if (std::abs(R)<=std::abs(E1)*1.0e-15) {break;}
576 }
577
578 E1=-0.5772156649015328-G4Log(X)+X*E1;
579 }
580 else
581 {
582 M=20+std::trunc(80.0/X);
583 T0=0.;
584
585 for(int K=M; K>=1; K--)
586 {
587 T0=K/(1.0+K/(X+T0));
588 }
589
590 T=1.0/(X+T0);
591 E1=std::exp(-X)*T;
592 }
593
594 return E1;
595 }
596