| 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 #include "G4VChannelingFastSimCrystalData.hh"
32 #include "G4SystemOfUnits.hh"
33 #include "G4PhysicalConstants.hh"
34 #include "G4Log.hh"
35
36 G4VChannelingFastSimCrystalData::G4VChanneling
37 {
38
39 }
40
41 //....oooOO0OOooo........oooOO0OOooo........oo
42
43 G4VChannelingFastSimCrystalData::~G4VChannelin
44
45 //....oooOO0OOooo........oooOO0OOooo........oo
46
47 void G4VChannelingFastSimCrystalData::SetGeome
48 (const G4
49 {
50 G4int crystalID = crystallogic->GetInstanc
51
52 //set bending angle if the it exists in th
53 (fMapBendingAngle.count(crystalID) > 0)
54 ? SetBendingAngle(fMapBendingAngle[crystal
55 : SetBendingAngle(0.,crystallogic);
56
57 //set miscut angle if the it exists in the
58 (fMapMiscutAngle.count(crystalID) > 0)
59 ? SetMiscutAngle(fMapMiscutAngle[crystalID
60 : SetMiscutAngle(0.,crystallogic);
61
62 //set crystalline undulator parameters if
63 //otherwise default = G4ThreeVector(0,0,0)
64 (fMapCUAmplitudePeriodPhase.count(crystalI
65 ? SetCUParameters(fMapCUAmplitudePeriodPha
66 : SetCUParameters(G4ThreeVector(0.,0.,0.),
67 }
68
69 //....oooOO0OOooo........oooOO0OOooo........oo
70
71 void G4VChannelingFastSimCrystalData::SetBendi
72 const G4LogicalVolume* cryst
73 {
74 G4int crystalID = crystallogic->GetInstanc
75
76 //set the bending angle for this logical v
77 fMapBendingAngle[crystalID]=tetab;
78
79 G4ThreeVector limboxmin;//minimal limits o
80 G4ThreeVector limboxmax;//maximal limits o
81 //save the limits of the box bounding the
82 crystallogic->GetSolid()->BoundingLimits(l
83
84 //bounding box half dimensions
85 fHalfDimBoundingBox = (limboxmax-limboxmin
86
87 G4double lcr = limboxmax.getZ()-limboxmin.
88
89 fBendingAngle=std::abs(tetab);
90 if (fBendingAngle<0.000001)//no bending le
91 {
92 if(fBendingAngle>DBL_EPSILON)
93 {
94 G4cout << "Channeling model: volum
95 G4cout << "Warning: bending angle
96 }
97
98 fBent=0;
99 fBendingAngle=0.;
100 fBendingR=0.;//just for convenience (in
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
118 G4cout << "Warning: bending angle i
119 }
120 }
121
122 }
123
124 //....oooOO0OOooo........oooOO0OOooo........oo
125
126 void G4VChannelingFastSimCrystalData::SetMiscu
127 const G4LogicalVolume *crysta
128 {
129 G4int crystalID = crystallogic->GetInstanc
130
131 //set the bending angle for this logical v
132 fMapMiscutAngle[crystalID]=tetam;
133
134 // fMiscutAngle>0: rotation of xz coordina
135 fMiscutAngle=tetam;
136 if (std::abs(tetam)>1.*mrad)
137 {
138 G4cout << "Channeling model: volume "
139 G4cout << "Warning: miscut angle is hi
140 G4cout << "coordinate transformation r
141 }
142 fCosMiscutAngle=std::cos(fMiscutAngle);
143 fSinMiscutAngle=std::sin(fMiscutAngle);
144 }
145
146 //....oooOO0OOooo........oooOO0OOooo........oo
147
148 void G4VChannelingFastSimCrystalData::SetCryst
149 G4doubl
150 G4doubl
151 G4doubl
152 const G
153 {
154 if (amplitude<DBL_EPSILON||period<DBL_EPSI
155 {
156 amplitude = 0.;
157 period=0.;
158 phase=0.;
159 G4cout << "Channeling model: volume "
160 G4cout << "Warning: The crystalline un
161 "=> the crystalline undulato
162 }
163
164 SetCUParameters(G4ThreeVector(amplitude,pe
165 }
166
167 //....oooOO0OOooo........oooOO0OOooo........oo
168
169 void G4VChannelingFastSimCrystalData::SetCUPar
170 const G4
171 const G4
172 {
173 G4int crystalID = crystallogic->GetInstanc
174
175 //set the crystalline undulator parameters
176 fMapCUAmplitudePeriodPhase[crystalID]=ampl
177 fCUAmplitude=amplitudePeriodPhase.x();
178 G4double period = amplitudePeriodPhase.y()
179 fCUPhase = amplitudePeriodPhase.z();
180
181 //if the amplidude of the crystalline undu
182 if(fCUAmplitude>DBL_EPSILON&&period>DBL_EP
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 cryst
192 SetBendingAngle(0,crystallogic);
193
194 G4cout << "Channeling model: volum
195 G4cout << "Warning: crystalline un
196 "a bent crystal mode =>
197 }
198 }
199 else
200 {
201 fCU = false;
202 fCUAmplitude = 0.;
203 fCUK = 0.;
204 fCUPhase = 0.;
205 fMapCUAmplitudePeriodPhase[crystalID]
206 }
207
208 fCUK2 = fCUK*fCUK;
209 fCUAmplitudeK = fCUAmplitude*fCUK;
210 }
211
212 //....oooOO0OOooo........oooOO0OOooo........oo
213
214 void G4VChannelingFastSimCrystalData::SetParti
215
216
217
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; //
226 fPz=std::sqrt(t); // momentum of parti
227 fPV=t/etotal; // pv
228 fBeta=fPz/etotal; // velocity/c
229 fTetaL = std::sqrt(std::abs(fZ2)*fVmax2
230 fChannelingStep = fChangeStep/fTetaL; /
231
232 // Energy losses
233 fV2 = fBeta*fBeta; // particle (velocit
234 fGamma = etotal/mass; // Lorentz facto
235 fMe2Gamma = 2*CLHEP::electron_mass_c2*f
236 // max ionization losses
237 fTmax = fMe2Gamma*fGamma*fV2/
238 (CLHEP::electron_mass_c2/mass*C
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 s
247 // defining by shielding by electrons
248 // teta1=hdc/(fPz*fRF)*DSQRT(1.13D0+3.76
249 teta1=fTeta10[i]*std::sqrt(1.13+fK40[
250
251
252 // the coefficient for multiple scatteri
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 incohere
257 // by the atomic correlations in crystal
258 // (screened atomic potential): EXP(-(fP
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 s
265 // defining by nucleus radius
266 // tetamax=hc/(fPz*1.D-6*fR0*fAN**(1.D0/
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 scatte
272 // fK2=(fZ2)**2*alphahbarc2*4.*pi*fN0*(f
273 fK2[i]=fK20[i]*zz22/fPV/fPV;
274 }
275
276 // fK3=(fZ2)**2*alphahbarc2*pi/electron_ma
277 fK3=fK30*zz22/fV2;
278 }
279
280 //....oooOO0OOooo........oooOO0OOooo........oo
281
282 G4double G4VChannelingFastSimCrystalData::GetL
283
284
285 {
286 G4double pv0 = etotal-mass*mass/etotal;
287 return std::sqrt(2*std::abs(charge)*fVm
288 //(!!! the
289 }
290
291 //....oooOO0OOooo........oooOO0OOooo........oo
292
293 G4double G4VChannelingFastSimCrystalData::GetL
294 {
295 return fTetaL; //return the Lindhard angle
296 }
297
298 //....oooOO0OOooo........oooOO0OOooo........oo
299
300 G4double G4VChannelingFastSimCrystalData::GetS
301 {
302 G4double simulationstep;
303 //find angle of particle w.r.t. the plane
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 ang
315 if (angle<fTetaL)
316 {
317 simulationstep = fChannelingStep;
318 }
319 else
320 {
321 simulationstep = fChangeStep;
322 if (angle > 0.0) { simulationstep /= ang
323 }
324
325 return simulationstep;
326 }
327
328 //....oooOO0OOooo........oooOO0OOooo........oo
329
330 G4double G4VChannelingFastSimCrystalData::GetM
331
332
333 {
334 //standard value of step for channeling pa
335 return fChangeStep/GetLindhardAngle(etotal
336 }
337
338 //....oooOO0OOooo........oooOO0OOooo........oo
339
340 G4ThreeVector G4VChannelingFastSimCrystalData:
341
342
343
344 {
345 G4double tx = 0.;//horizontal scatteri
346 G4double ty = 0.;//vertical scattering
347
348 G4double ksi=0.1;
349
350 // calculation of the teta2-minimal possi
351 G4double e1=fK2[ielement]*effectiveSte
352 // (real formula is (4*pi*fN0*wpl(x)*dz*f
353 G4double teta122=fTetamax12[ielement]/
354 // teta122=fTeta12+teta22=teta1^2+teta
355
356 G4double teta22;
357 G4double t;
358 // if the angle of a single scattering is
359 // angle of coulomb scattering defining b
360 // multiple scattering by both nuclei and
361 // occur => minimal possible angle of a s
362 if (teta122<=fTeta12[ielement]*1.00012
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*fB
376 fBBDEXP[ielement]*
377 (expint(fBB[ielemen
378
379 // sumilation of multiple coulomb scatt
380 // for high speed of a program, real fo
381 // 4*pi*fN0*wpl(x)*dz*(fZ1*zz2*alpha*hd
382 // *(ln(1+a)+(1-exp(-a*b))/(1+a)+(1+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 scatterin
394 G4double zss=0.;
395 G4double dzss=step;
396
397 // (calculation of a distance, at which a
398 ksi=G4UniformRand();
399
400 zss=-G4Log(ksi)*step/(e1*(1./teta122-1
401 G4double tt;
402
403 // At some step several single scattering
404 // So,if the distance of the next scatter
405 // another scattering can occur. If the d
406 // is less than the difference between th
407 // the previous scattering, another scatt
408 // In the cycle we simulate each of them.
409 // the remaining part of step is less tha
410 //********************************************
411 // if at a step a single scattering occur
412 while (zss<dzss)
413 {
414
415 // simulation by Monte-Carlo of angles
416 ksi=G4UniformRand();
417
418 tt=fTetamax12[ielement]/(1.+ksi*(fTe
419 fTeta12[ielement];
420
421 ksi=G4UniformRand();
422
423 // suppression of incoherent scattering
424 t=fPzu11[ielement]*tt;
425 t=std::exp(-t);
426
427 if (t<ksi) //if scattering takes pla
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
439 ksi=G4UniformRand();
440
441 zss=-G4Log(ksi)*step/(e1*(1./teta122
442 }
443 //********************************************
444 return G4ThreeVector(tx,ty,0.);
445 }
446
447 //....oooOO0OOooo........oooOO0OOooo........oo
448
449 G4ThreeVector G4VChannelingFastSimCrystalData:
450
451
452
453 {
454
455 G4double zss=0.;
456 G4double dzss=step;
457 G4double ksi = 0.;
458
459 G4double tx = 0.;//horizontal scattering a
460 G4double ty = 0.;//vertical scattering ang
461 G4double eloss = 0.;//energy loss
462
463 // eMinIonization - minimal energy transfe
464 // a cut to reduce the number of calls of
465 // is needed only at low density regions,
466 if (eMinIonization<0.5*eV){eMinIonization=
467
468 // single scattering on electrons routine
469 if ((eMinIonization<fTmax)&&(electronDensi
470 {
471
472 // (calculation of a distance, at which ano
473 // simulation of scattering length (by the
474 ksi=G4UniformRand();
475
476 zss=-1.0*G4Log(ksi)/(fK3*electronDensity
477
478 //********************************************
479 // if at a step a single scattering occur
480 while (zss<dzss)
481 {
482 // simulation by Monte-Carlo of angles of
483 ksi=G4UniformRand();
484
485 // energy transfered to electron
486 G4double e1=eMinIonization/(1.-ksi*(1.
487
488 // scattering angle
489 G4double t=0;
490 if(fTmax-e1>DBL_EPSILON) //to be sure
491 {
492 t=std::sqrt(2.*CLHEP::electron_mas
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 a
504 // simulation of scattering length
505 // (by the same way single scattering by
506 ksi=G4UniformRand();
507
508 zss=-1.0*G4Log(ksi)/(fK3*electronDensi
509 }
510 //********************************************
511 }
512 return G4ThreeVector(tx,ty,eloss);
513 }
514
515 //....oooOO0OOooo........oooOO0OOooo........oo
516
517 G4double G4VChannelingFastSimCrystalData::Ioni
518
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
526 G4double delta= 2*(G4Log(fBeta*fGamma)+fLo
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*G4
533 1/8*((fGamma - 1)/fGamma)*((fGa
534 }
535 else if(fParticleName=="e+")
536 {
537 loge+=(-fV2/12*(11 + 14/(fGamma + 1) +
538 4/(fGamma + 1)/(fGamma +
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........oo
549
550 G4double G4VChannelingFastSimCrystalData::expi
551 {
552 // ======================================
553 // Purpose: Compute exponential integral
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)
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