| 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
29 #include "G4ChannelingFastSimCrystalData.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4PhysicalConstants.hh"
32
33 G4ChannelingFastSimCrystalData::G4ChannelingFa
34 {
35
36 }
37
38 //....oooOO0OOooo........oooOO0OOooo........oo
39
40 void G4ChannelingFastSimCrystalData::SetMateri
41 const G4M
42 const G4S
43 const G4S
44 {
45 G4String filename=crystal->GetName(); //in
46 filename.erase(0,3);
47
48 if (fVerbosity)
49 {
50 G4cout <<
51 "=========================================
52 << G4endl;
53 G4cout <<
54 "====== Crystal lattice
55 << G4endl;
56 G4cout <<
57 "=========================================
58 << G4endl;
59 G4cout << "Crystal material: " << filename <
60 }
61
62 //choice between planes (1D model) and axe
63 if (lattice.compare(0,1,"(")==0)
64 {
65 iModel=1; //planes
66 filename = filename + "_planes_"; //temp
67 if (fVerbosity)
68 G4cout << "Crystal planes: " << lattice << G
69 }
70 else if (lattice.compare(0,1,"<")==0)
71 {
72 iModel=2; //axes
73 filename = filename + "_axes_"; //tempor
74 if (fVerbosity)
75 G4cout << "Crystal axes: " << lattice << G4e
76 }
77
78 //input file:
79 filename = filename + lattice.substr(1,(la
80
81 if(filePath=="")
82 {
83 //standard file path if another one is n
84 filename = "/" + filename;
85 filename = G4FindDataDir("G4CHANNELINGDA
86 }
87 else
88 {
89 //custom file path
90 filename = filePath + filename;
91 }
92
93 fNelements=(G4int)crystal->GetNumberOfElem
94 for(G4int i=0; i<fNelements; i++)
95 {
96 fZ1.push_back(crystal->GetElement(i)->Ge
97 fAN.push_back(crystal->GetElement(i)->Ge
98 fI0.push_back(crystal->GetElement(i)->Ge
99 }
100
101 G4double var;//just variable
102 G4double unitIF;//unit of interpolation fu
103
104 std::ifstream vfilein;
105 vfilein.open(filename);
106 //check if the input file was found, other
107 if(!vfilein.is_open())
108 {
109 G4String outputMessage="Input file " +
110 filename +
111 " is not found!
112 G4Exception("SetMaterialProperties",
113 "001",
114 FatalException,
115 outputMessage);
116 }
117
118 //read nuclear concentration
119 for(G4int i=0; i<fNelements; i++)
120 {
121 vfilein >> var;
122 fN0.push_back(var/cm3);
123 }
124
125 //read amplitude of thermal oscillations
126 for(G4int i=0; i<fNelements; i++)
127 {
128 vfilein >> var;
129 fU1.push_back(var*cm);
130 }
131
132 if (iModel==1)
133 {
134 // read channel dimensions
135 vfilein >> fDx;
136 fDx*=cm;
137 // read interpolation step size
138 vfilein >> fNpointsx;
139
140 fDy = fDx;
141 fNpointsy = 0;
142 }
143 else if (iModel==2)
144 {
145 // read channel dimensions
146 vfilein >> fDx >> fDy;
147 fDx*=cm;
148 fDy*=cm;
149 // read the number of nodes of interpol
150 vfilein >> fNpointsx >> fNpointsy;
151 }
152
153 //read the height of the potential well, n
154 vfilein >> fVmax;
155 fVmax*=eV;
156 fVmax2=2.*fVmax;
157
158 //read the on-zero minimal potential insid
159 //necessary for angle recalculation for en
160 //the crystal lateral surface
161 vfilein >> fVMinCrystal;
162 fVMinCrystal*=eV;
163
164 // to create the class of interpolation fo
165 fElectricFieldX =
166 new G4ChannelingFastSimInterpolati
167 if(iModel==2) {fElectricFieldY =
168 new G4ChannelingFastSimInterpolati
169 fElectronDensity =
170 new G4ChannelingFastSimInterpolati
171 fMinIonizationEnergy =
172 new G4ChannelingFastSimInterpolati
173
174 // do it for any element of crystal materi
175 for(G4int i=0; i<fNelements; i++)
176 {
177 fNucleiDensity.push_back(
178 new G4ChannelingFastSimInterpolati
179 }
180
181 if (iModel==1)
182 {
183 G4double ai, bi, ci, di;
184 for(G4int i=0; i<fNpointsx; i++)
185 {
186 //reading the coefficients of cubic sp
187 vfilein >> ai >> bi >> ci >> di;
188 //setting spline coefficients for elec
189 unitIF=eV/cm;
190 fElectricFieldX->SetCoefficients1D(ai*
191 ci*
192
193 //reading the coefficients of cubic sp
194 vfilein >> ai >> bi >> ci >> di;
195 //setting spline coefficients for nucl
196 unitIF=1.;
197 fNucleiDensity[0]->SetCoefficients1D(a
198 c
199
200 //reading the coefficients of cubic sp
201 vfilein >> ai >> bi >> ci >> di;
202 //setting spline coefficients for elec
203 unitIF=1./cm3;
204 fElectronDensity->SetCoefficients1D(ai
205 ci
206
207 //reading the coefficients of cubic sp
208 vfilein >> ai >> bi >> ci >> di;
209 //setting spline coefficients for mini
210 unitIF=eV;
211 fMinIonizationEnergy->SetCoefficients1
212
213
214 for(G4int ii=1; ii<fNelements; ii++)
215 {
216 //reading the coefficients of cubi
217 vfilein >> ai >> bi >> ci >> di;
218 //setting spline coefficients for
219 unitIF=1.;
220 fNucleiDensity[ii]->SetCoefficient
221
222 }
223 }
224 }
225 else if (iModel==2)
226 {
227 G4double ai3D, bi3D, ci3D;
228 for(G4int j=0; j<fNpointsy; j++)
229 {
230 for(G4int i=0; i<fNpointsx+1; i++)
231 {
232 for(G4int k=0; k<2; k++)
233 {
234 //reading the coefficients of cubi
235 vfilein >> ai3D >> bi3D >> ci3D;
236 unitIF=eV;
237 //setting spline coefficients for
238 fMinIonizationEnergy->SetCoefficie
239
240
241 //reading the coefficients of cubi
242 vfilein >> ai3D >> bi3D >> ci3D;
243 //setting spline coefficients for
244 unitIF=eV/cm;
245 fElectricFieldX->SetCoefficients2D
246
247
248 //reading the coefficients of cubi
249 vfilein >> ai3D >> bi3D >> ci3D;
250 //setting spline coefficients for
251 unitIF=eV/cm;
252 fElectricFieldY->SetCoefficients2D
253
254
255 //reading the coefficients of cubi
256 vfilein >> ai3D >> bi3D >> ci3D;
257 //setting spline coefficients for
258 unitIF=1.;
259 fNucleiDensity[0]->SetCoefficients
260
261
262 //reading the coefficients of cubi
263 vfilein >> ai3D >> bi3D >> ci3D;
264 //setting spline coefficients for
265 unitIF=1./cm3;
266 fElectronDensity->SetCoefficients2
267
268
269 for(G4int ii=1; ii<fNelements; ii+
270 {
271 //reading the coefficients of
272 vfilein >> ai3D >> bi3D >> ci3
273 //setting spline coefficients
274 unitIF=1.;
275 fNucleiDensity[ii]->SetCoeffic
276
277
278 }
279
280 }
281 }
282 }
283 }
284
285 vfilein.close();
286
287 //set special values and coefficients
288 G4double alphahbarc2=std::pow(CLHEP::fine_
289 fK30=2.*CLHEP::pi*alphahbarc2/CLHEP::elect
290
291 for(G4int i=0; i<fNelements; i++)
292 {
293 fRF.push_back((std::pow(9*CLHEP::pi*CLHE
294 *0.5291772109217*angstrom)
295
296 fTetamax0.push_back(CLHEP::hbarc/(fR0*st
297 fTeta10.push_back(CLHEP::hbarc/fRF[i]);
298 fPu11.push_back(std::pow(fU1[i]/CLHEP::h
299
300 fK20.push_back(alphahbarc2*4*CLHEP::pi*f
301
302 fK40.push_back(3.76*std::pow(CLHEP::fine
303
304 fKD.push_back(fK30*fZ1[i]*fN0[i]);
305
306 fLogPlasmaEdI0.push_back(G4Log((crystal-
307 }
308
309 fBB.resize(fNelements);
310 fE1XBbb.resize(fNelements);
311 fBBDEXP.resize(fNelements);
312 fPzu11.resize(fNelements);
313 fTeta12.resize(fNelements);
314 fTetamax2.resize(fNelements);
315 fTetamax12.resize(fNelements);
316 fK2.resize(fNelements);
317
318 fChangeStep = CLHEP::pi*std::min(fDx,fDy)/
319 //ste
320 // =
321 }
322
323 //....oooOO0OOooo........oooOO0OOooo........oo
324
325 G4ThreeVector G4ChannelingFastSimCrystalData::
326
327 {
328 G4double x0=pos0.x(),y0=pos0.y(),z0=pos0.z(
329 z0+=fHalfDimBoundingBox.z();
330 G4double x,y,z;
331
332 if (fBent)
333 {
334 // for bent crystal
335 G4double rsqrt = std::sqrt(fBendingRsqu
336 fBending2R*(
337 x0*x0 + z0*z
338 //transform to co-rotating reference sy
339 x = fBendingR - rsqrt;
340 y = y0;
341 z = fBendingR*std::asin((z0*fCosMiscutA
342 }
343 else
344 {
345 //for straight crystal
346 x = x0*fCosMiscutAngle - z0*fSinMiscutA
347 y = y0;
348 z = x0*fSinMiscutAngle + z0*fCosMiscutA
349
350 //for crystalline undulator
351 if(fCU){x-=GetCUx(z);}
352 }
353
354 //calculation of coordinates within a chann
355 fNChannelx=std::floor(x/fDx); //remember th
356 //to track th
357 x-=fNChannelx*fDx;
358 fNChannely=std::floor(y/fDy);//remember the
359 //to track the
360 y-=fNChannely*fDy;
361 //correction of the longitudinal coordinate
362 if (fBent) {fCorrectionZ = fBendingR/(fBend
363
364 return G4ThreeVector(x,y,z);
365 }
366
367 //....oooOO0OOooo........oooOO0OOooo........oo
368
369 G4ThreeVector G4ChannelingFastSimCrystalData::
370 const G4ThreeVector &pos)
371 {
372 G4double x=pos.x(),y=pos.y(),z=pos.z();
373
374 //transform to co-rotating reference system
375 x+=fNChannelx*fDx;
376 y+=fNChannely*fDy;
377
378 G4double x0,y0,z0;
379
380 if (fBent)
381 {
382 // for bent crystal
383 G4double rcos = (fBendingR - x)*(1. - s
384 G4double a = x + rcos;
385 G4double b = std::sqrt(x*x + fBending2R
386
387 //transform to Box coordinates
388 x0 = a*fCosMiscutAngle + b*fSinMiscutAn
389 y0 = y;
390 z0 = b*fCosMiscutAngle - a*fSinMiscutAn
391 }
392 else
393 {
394 //for crystalline undulator
395 if(fCU){x+=GetCUx(z);}
396
397 //for straight crystal
398 x0 = x*fCosMiscutAngle + z*fSinMiscutAn
399 y0 = y;
400 z0 =-x*fSinMiscutAngle + z*fCosMiscutAn
401 }
402
403 return G4ThreeVector(x0,y0,z0-fHalfDimBound
404 }
405
406 //....oooOO0OOooo........oooOO0OOooo........oo
407
408 G4ThreeVector G4ChannelingFastSimCrystalData::
409 G4double& y,
410 G4double& z)
411 {
412
413 //test of enter in other channel
414 if (x<0)
415 {
416 fNChannelx-=1;
417 x+=fDx; //enter in other channel
418 //correction of the longitudinal coord
419 if (fBent) {fCorrectionZ = fBendingR/(
420 }
421 else if (x>=fDx)
422 {
423 fNChannelx+=1;
424 x-=fDx; //enter in other channel
425 //correction of the longitudinal coord
426 if (fBent) {fCorrectionZ = fBendingR/(
427 }
428
429 //test of enter in other channel
430 if (y<0)
431 {
432 fNChannely-=1;
433 y+=fDy; //enter in other channel
434 }
435 else if (y>=fDy)
436 {
437 fNChannely+=1;
438 y-=fDy; //enter in other channel
439 }
440
441 return G4ThreeVector(x,y,z);
442 }
443
444 //....oooOO0OOooo........oooOO0OOooo........oo
445