| 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 #ifndef G4VChannelingFastSimCrystalData_h
32 #define G4VChannelingFastSimCrystalData_h 1
33
34 #include "G4ios.hh"
35 #include "globals.hh"
36 #include "G4ThreeVector.hh"
37 #include "Randomize.hh"
38 #include "G4LogicalVolume.hh"
39 #include "G4Material.hh"
40 #include "G4VSolid.hh"
41 #include <unordered_map>
42
43 #include "G4ChannelingFastSimInterpolation.hh"
44
45 /** \file G4VChannelingFastSimCrystalData.hh
46 * \brief Definition of the G4VChannelingFastSimCrystalData class
47 * The class contains the data and properties related to the crystal lattice as well as
48 * functions to simulate of important physical processes, i.e. coulomb scattering on
49 * screened atomic potential, on single electrons and ionization energy losses;
50 * functions of electric fields, nuclear and electron densities and minimum energy
51 * of ionization (the corresponding interpolation coefficients are in
52 * G4ChannelingFastSimInterpolation).
53 * The functions related to the crystal geometry (transformation of coordinates and angles
54 * from the reference system of the bounding box of the local volume to
55 * the crystal lattice co-rotating reference system and vice versa) and
56 * initialization function SetMaterialProperties are created as virtual to make
57 * material data input and geometry functions flexible for modification.
58 */
59
60 class G4VChannelingFastSimCrystalData{
61 public:
62
63 G4VChannelingFastSimCrystalData();
64 virtual ~G4VChannelingFastSimCrystalData();
65
66 ///electric fields produced by crystal lattice
67 G4double Ex(G4double x,G4double y) {return (fElectricFieldX->GetIF(x,y))*(-fZ2/fPV);}
68 G4double Ey(G4double x,G4double y) {return (fElectricFieldY->GetIF(x,y))*(-fZ2/fPV);}
69
70 ///electron density function
71 G4double ElectronDensity(G4double x,G4double y)
72 {
73 G4double nel0=fElectronDensity->GetIF(x,y);
74 if(nel0<0.) {nel0=0.;}//exception, errors of interpolation functions
75 return nel0;
76 }
77 ///minimum energy of ionization function
78 G4double MinIonizationEnergy(G4double x,G4double y)
79 {return fMinIonizationEnergy->GetIF(x,y);}
80 ///nuclear density function (normalized to average nuclear density)
81 G4double NuclearDensity(G4double x,G4double y, G4int ielement)
82 {return std::abs(fNucleiDensity[ielement]->GetIF(x,y));}
83 //abs to describe exception, errors of interpolation functions,
84 //don't put it =0, otherwise division on 0 in CoulombAtomicScattering
85
86 ///Calculate the value of the Lindhard angle (!!! the value for a straight crystal)
87 G4double GetLindhardAngle(G4double etotal, G4double mass, G4double charge);
88 ///Calculate the value of the Lindhard angle (!!! the value for a straight crystal)
89 G4double GetLindhardAngle();//return the Lindhard angle value calculated in
90 //SetParticleProperties
91
92 ///Calculate simulation step (standard value for channeling particles and
93 ///reduced value for overbarrier particles)
94 G4double GetSimulationStep(G4double tx,G4double ty);
95 ///Calculate maximal simulation step (standard value for channeling particles)
96 G4double GetMaxSimulationStep(G4double etotal, G4double mass, G4double charge);
97
98 ///get particle velocity/c
99 G4double GetBeta(){return fBeta;}
100
101 G4int GetNelements() {return fNelements;}
102 G4int GetModel() {return iModel;}//=1 for planes, =2 for axes
103
104 ///get bending angle of the crystal planes/axes
105 ///(default BendingAngle=0 => straight crystal);
106 G4double GetBendingAngle(){return fBendingAngle;}
107 ///fBendingAngle MAY BE NOT THE SAME AS THE BENDING ANGLE OF THE CRYSTAL VOLUME:
108 ///THE VOLUME OF A BENT CRYSTAL MAY BE G4Box, while the planes/axes inside may be bent
109
110 G4double GetMiscutAngle(){return fMiscutAngle;}
111
112 ///get crystal curvature
113 ///for crystalline undulator the curvature is a function, otherwise it's a constant
114 G4double GetCurv(G4double z){return fCU ? -fCUK2*GetCUx(z) : fCurv;}
115
116 ///get crystalline undulator wave function
117 G4double GetCUx(G4double z){return fCUAmplitude*std::cos(fCUK*z+fCUPhase);}
118 ///get crystalline undulator wave 1st derivative function
119 G4double GetCUtetax(G4double z){
120 return fCU ? -fCUAmplitudeK*std::sin(fCUK*z+fCUPhase) : 0;}
121
122 ///find and upload crystal lattice input files, calculate all the basic values
123 ///(to do only once)
124 virtual void SetMaterialProperties(const G4Material* crystal,
125 const G4String &lattice,
126 const G4String &filePath) = 0;
127
128 ///set geometry parameters from current logical volume
129 void SetGeometryParameters(const G4LogicalVolume *crystallogic);
130
131 ///set bending angle of the crystal planes/axes
132 ///(default fBendingAngle=0 => straight crystal);
133 ///only non-negative values! crystal is bent in the positive direction of x
134 void SetBendingAngle(G4double tetab, const G4LogicalVolume *crystallogic);
135 ///fBendingAngle MAY BE NOT THE SAME AS THE BENDING ANGLE OF THE CRYSTAL VOLUME
136 ///THE VOLUME OF A BENT CRYSTAL MAY BE G4Box, while the planes/axes inside may be bent
137
138 ///set miscut angle (default fMiscutAngle=0), acceptable range +-1 mrad,
139 ///otherwise geometry routines may be unstable
140 void SetMiscutAngle(G4double tetam, const G4LogicalVolume *crystallogic);
141
142 ///set crystalline undulator parameters: amplitude, period and phase
143 /// (default: all 3 value = 0)
144 /// function to use in Detector Construction
145 void SetCrystallineUndulatorParameters(G4double amplitude,
146 G4double period,
147 G4double phase,
148 const G4LogicalVolume *crystallogic);
149
150 ///set crystalline undulator parameters (internal function of the model)
151 ///for convenience we put amplitude, period and phase in a G4ThreeVector
152 void SetCUParameters(const G4ThreeVector &litudePeriodPhase,
153 const G4LogicalVolume *crystallogic);
154
155 ///recalculate all the important values
156 ///(to do both at the trajectory start and after energy loss)
157 void SetParticleProperties(G4double etotal,
158 G4double mp,
159 G4double charge,
160 const G4String& particleName);
161
162 ///calculate the coordinates in the co-rotating reference system
163 ///within a channel (periodic cell)
164 ///(connected with crystal planes/axes either bent or straight)
165 virtual G4ThreeVector CoordinatesFromBoxToLattice(const G4ThreeVector &pos0) = 0;
166
167 ///calculate the coordinates in the Box reference system
168 ///(connected with the bounding box of the volume)
169 virtual G4ThreeVector CoordinatesFromLatticeToBox(const G4ThreeVector &pos) = 0;
170
171 ///change the channel if necessary, recalculate x o y
172 virtual G4ThreeVector ChannelChange(G4double& x, G4double& y, G4double& z) = 0;
173
174 ///return correction of the longitudinal coordinate
175 /// (along current plane/axis vs "central plane/axis")
176 G4double GetCorrectionZ(){return fCorrectionZ;}
177
178 ///calculate the horizontal angle in the co-rotating reference system
179 ///within a channel (periodic cell)
180 ///(connected with crystal planes/axes either bent or straight)
181 virtual G4double AngleXFromBoxToLattice(G4double tx, G4double z)=0;
182
183 ///calculate the horizontal angle in the Box reference system
184 ///(connected with the bounding box of the volume)
185 virtual G4double AngleXFromLatticeToBox(G4double tx, G4double z)=0;
186
187 ///auxialiary function to transform the horizontal angle
188 virtual G4double AngleXShift(G4double z)=0;
189
190 ///multiple and single scattering on screened potential
191 G4ThreeVector CoulombAtomicScattering(
192 G4double effectiveStep,
193 G4double step,
194 G4int ielement);
195 ///multiple and single scattering on electrons
196 G4ThreeVector CoulombElectronScattering(G4double eMinIonization,
197 G4double electronDensity,
198 G4double step);
199 ///ionization losses
200 G4double IonizationLosses(G4double dz, G4int ielement);
201
202 void SetVerbosity(G4int ver){fVerbosity = ver;}
203
204 protected:
205 ///classes containing interpolation coefficients
206 //horizontal electric field data
207 G4ChannelingFastSimInterpolation* fElectricFieldX{nullptr};
208 //vertical electric field data
209 G4ChannelingFastSimInterpolation* fElectricFieldY{nullptr};
210 //electron density data
211 G4ChannelingFastSimInterpolation* fElectronDensity{nullptr};
212 //minimal energy of ionization data
213 G4ChannelingFastSimInterpolation* fMinIonizationEnergy{nullptr};
214 //nuclear density distributions data
215 std::vector <G4ChannelingFastSimInterpolation*> fNucleiDensity;
216
217 ///values related to the crystal geometry
218 G4ThreeVector fHalfDimBoundingBox;//bounding box half dimensions
219 G4int fBent=0;//flag of bent crystal,
220 //=0 for straight and =1 for bent, by default straight crystal
221
222 G4double fBendingAngle=0.;// angle of bending of the crystal planes/axes
223 //inside the crystal volume
224 //MAY BE NOT THE SAME AS THE BENDING ANGLE OF THE CRYSTAL VOLUME
225 //THE VOLUME OF A BENT CRYSTAL MAY BE G4Box,
226 //while the planes/axes inside may be bent
227 G4double fBendingR = 0.; // bending radius of the crystal planes/axes
228 G4double fBending2R=0.; // =2*fBendingR
229 G4double fBendingRsquare=0.; // =fBendingR**2
230 G4double fCurv=0.; //=1/fBendingR bending curvature of the crystal planes/axes
231
232 G4double fMiscutAngle = 0.;// miscut angle, can be of either sign or 0;
233 //safe values |ThetaMiscut|<0.001
234 G4double fCosMiscutAngle=1.;// = std::cos(fMiscutAngle), to economy operations
235 G4double fSinMiscutAngle=0.;// = std::sin(fMiscutAngle), to economy operations
236
237 G4double fCorrectionZ = 1.;//correction of the longitudinal coordinate
238 //(along current plane/axis vs "central plane/axis"), 1 is default value
239 //(for "central plane/axis" or a straight crystal)
240
241 G4bool fCU = false;//flag of crystalline undulator geometry
242 //(periodically bent crystal)
243 G4double fCUAmplitude=0.; //Amplitude of a crystalline undulator
244 G4double fCUK=0.; //2*pi/period of a crystalline undulator
245 G4double fCUPhase=0.;//Phase of a crystalline undulator
246 G4double fCUAmplitudeK=0.;//fCUAmplitude*fCUK
247 G4double fCUK2=0.; //fCUK^2
248
249 ///values related to the crystal lattice
250 G4int fNelements=1;//number of nuclear elements in a crystal
251 G4int iModel=1;// model type (iModel=1 for interplanar potential,
252 //iModel=2 for the interaxial one)
253
254 G4double fVmax=0; // the height of the potential well
255 G4double fVmax2=0; // =2*fVmax
256 G4double fVMinCrystal=0;// non-zero minimal potential inside the crystal,
257 // necessary for angle recalculation for entrance/exit
258 //through the crystal lateral surface
259
260 G4double fChangeStep=0;// fChannelingStep = fChangeStep/fTetaL
261
262 std::vector <G4double> fI0; //Mean excitation energy
263
264 std::vector <G4double> fRF;//Thomas-Fermi screening radius
265
266 ///angles necessary for multiple and single coulomb scattering
267
268 //minimal scattering angle by coulomb scattering on nuclei
269 //defined by shielding by electrons
270 std::vector <G4double> fTeta10;//(in the Channeling model
271 //teta1=fTeta10/fPz*(1.13+fK40/vz**2)
272 //maximal scattering angle by coulomb scattering on nuclei defined by nucleus radius
273 std::vector <G4double> fTetamax0;//(in the Channeling model tetamax=fTetamax0/fPz)
274 std::vector <G4double> fTetamax2;//=tetamax*tetamax
275 std::vector <G4double> fTetamax12;//=teta1*teta1+tetamax*tetamax
276 std::vector <G4double> fTeta12; //= teta1*teta1
277
278 ///coefficients necessary for multiple and single coulomb scattering
279 std::vector <G4double> fK20; //a useful coefficient, fK2=fK20/fPV/fPV
280 std::vector <G4double> fK2; //a useful coefficient,
281 //fK2=(fZ2*alpha*hdc)**2*4.*pi*fN0*(fZ1/fPV)**2
282 std::vector <G4double> fK40; //a useful coefficient, fK40=3.76D0*(alpha*fZ1)**2
283 G4double fK30=0;//a useful coefficient, fK3=fK30/fPV/fPV
284 G4double fK3=0;//a useful coefficient, fK3=2.*pi*alpha*hdc/electron_mass_c2/(fPV)**2
285
286 std::vector <G4double> fKD; //a useful coefficient for dE/dx
287 std::vector <G4double> fLogPlasmaEdI0; //item of delta-correction of ionization loss
288
289 ///coefficients for multiple scattering suppression
290 std::vector <G4double> fPu11;//a useful coefficient for exponent containing u1
291 std::vector <G4double> fPzu11;//a useful coefficient for exponent containing u1
292 std::vector <G4double> fBB;//a useful coefficient
293 std::vector <G4double> fE1XBbb;//a useful coefficient
294 std::vector <G4double> fBBDEXP;//a useful coefficient
295
296 //Variable to control printout
297 G4int fVerbosity = 1;
298
299
300 private:
301
302 //exponential integral
303 G4double expint(G4double x);
304
305 ///private variables
306
307 std::unordered_map<G4int, G4double> fMapBendingAngle;//the map fBendingAngle
308 //for different logical volumes
309
310 std::unordered_map<G4int, G4double> fMapMiscutAngle;//the map fMiscutAngle
311 //for different logical volumes
312
313 std::unordered_map<G4int, G4ThreeVector> fMapCUAmplitudePeriodPhase;//the map of
314 //AmplitudePeriodPhase
315 //for different logical volumes
316
317 G4double fChannelingStep=0;// simulation step under the channeling conditions =
318 //channeling oscillation length/fNsteps
319 // channeling oscillation length: Biryukov book Eq. (1.24)
320
321 ///energy depended values
322 G4double fPz=0; // particle momentum absolute value
323 G4double fPV=0; // pv
324 G4double fTetaL=0; //Lindhard angle
325 G4double fBeta=0; //particle (velocity/c)
326 G4double fV2=0; // particle (velocity/c)^2
327 G4double fGamma=0; //Lorentz factor
328 G4double fMe2Gamma=0; // me^2*fGamma
329 G4double fTmax=0; // max ionization losses
330
331 ///particle properties flags
332 G4String fParticleName = "";
333 G4double fZ2=0; //particle charge
334
335 };
336
337 #endif
338
339