| 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
29 #include "G4ChannelingFastSimCrystalData.hh"
30 #include "G4SystemOfUnits.hh"
31 #include "G4PhysicalConstants.hh"
32
33 G4ChannelingFastSimCrystalData::G4ChannelingFastSimCrystalData()
34 {
35
36 }
37
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39
40 void G4ChannelingFastSimCrystalData::SetMaterialProperties(
41 const G4Material *crystal,
42 const G4String &lattice,
43 const G4String &filePath)
44 {
45 G4String filename=crystal->GetName(); //input file
46 filename.erase(0,3);
47
48 if (fVerbosity)
49 {
50 G4cout <<
51 "======================================================================="
52 << G4endl;
53 G4cout <<
54 "====== Crystal lattice data ========"
55 << G4endl;
56 G4cout <<
57 "======================================================================="
58 << G4endl;
59 G4cout << "Crystal material: " << filename << G4endl;
60 }
61
62 //choice between planes (1D model) and axes (2D model)
63 if (lattice.compare(0,1,"(")==0)
64 {
65 iModel=1; //planes
66 filename = filename + "_planes_"; //temporary name
67 if (fVerbosity)
68 G4cout << "Crystal planes: " << lattice << G4endl;
69 }
70 else if (lattice.compare(0,1,"<")==0)
71 {
72 iModel=2; //axes
73 filename = filename + "_axes_"; //temporary name
74 if (fVerbosity)
75 G4cout << "Crystal axes: " << lattice << G4endl;
76 }
77
78 //input file:
79 filename = filename + lattice.substr(1,(lattice.length())-2) + ".dat";
80
81 if(filePath=="")
82 {
83 //standard file path if another one is not set
84 filename = "/" + filename;
85 filename = G4FindDataDir("G4CHANNELINGDATA") + filename;
86 }
87 else
88 {
89 //custom file path
90 filename = filePath + filename;
91 }
92
93 fNelements=(G4int)crystal->GetNumberOfElements();
94 for(G4int i=0; i<fNelements; i++)
95 {
96 fZ1.push_back(crystal->GetElement(i)->GetZ());
97 fAN.push_back(crystal->GetElement(i)->GetAtomicMassAmu());
98 fI0.push_back(crystal->GetElement(i)->GetIonisation()->GetMeanExcitationEnergy());
99 }
100
101 G4double var;//just variable
102 G4double unitIF;//unit of interpolation function
103
104 std::ifstream vfilein;
105 vfilein.open(filename);
106 //check if the input file was found, otherwise return an exception
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 interpolation
150 vfilein >> fNpointsx >> fNpointsy;
151 }
152
153 //read the height of the potential well, necessary only for step length calculation
154 vfilein >> fVmax;
155 fVmax*=eV;
156 fVmax2=2.*fVmax;
157
158 //read the on-zero minimal potential inside the crystal,
159 //necessary for angle recalculation for entrance/exit through
160 //the crystal lateral surface
161 vfilein >> fVMinCrystal;
162 fVMinCrystal*=eV;
163
164 // to create the class of interpolation for any function
165 fElectricFieldX =
166 new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);
167 if(iModel==2) {fElectricFieldY =
168 new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);}
169 fElectronDensity =
170 new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);
171 fMinIonizationEnergy =
172 new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);
173
174 // do it for any element of crystal material
175 for(G4int i=0; i<fNelements; i++)
176 {
177 fNucleiDensity.push_back(
178 new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel));
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 spline
187 vfilein >> ai >> bi >> ci >> di;
188 //setting spline coefficients for electric field
189 unitIF=eV/cm;
190 fElectricFieldX->SetCoefficients1D(ai*unitIF, bi*unitIF,
191 ci*unitIF, di*unitIF, i);
192
193 //reading the coefficients of cubic spline
194 vfilein >> ai >> bi >> ci >> di;
195 //setting spline coefficients for nuclear density (first element)
196 unitIF=1.;
197 fNucleiDensity[0]->SetCoefficients1D(ai*unitIF, bi*unitIF,
198 ci*unitIF, di*unitIF, i);
199
200 //reading the coefficients of cubic spline
201 vfilein >> ai >> bi >> ci >> di;
202 //setting spline coefficients for electron density
203 unitIF=1./cm3;
204 fElectronDensity->SetCoefficients1D(ai*unitIF, bi*unitIF,
205 ci*unitIF, di*unitIF, i);
206
207 //reading the coefficients of cubic spline
208 vfilein >> ai >> bi >> ci >> di;
209 //setting spline coefficients for minimal ionization energy
210 unitIF=eV;
211 fMinIonizationEnergy->SetCoefficients1D(ai*unitIF, bi*unitIF,
212 ci*unitIF, di*unitIF, i);
213
214 for(G4int ii=1; ii<fNelements; ii++)
215 {
216 //reading the coefficients of cubic spline
217 vfilein >> ai >> bi >> ci >> di;
218 //setting spline coefficients for nuclear density (other elements if any)
219 unitIF=1.;
220 fNucleiDensity[ii]->SetCoefficients1D(ai*unitIF, bi*unitIF,
221 ci*unitIF, di*unitIF, i);
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 cubic spline
235 vfilein >> ai3D >> bi3D >> ci3D;
236 unitIF=eV;
237 //setting spline coefficients for minimal ionization energy
238 fMinIonizationEnergy->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
239 i, j, k);
240
241 //reading the coefficients of cubic spline
242 vfilein >> ai3D >> bi3D >> ci3D;
243 //setting spline coefficients for horizontal electric field
244 unitIF=eV/cm;
245 fElectricFieldX->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
246 i, j, k);
247
248 //reading the coefficients of cubic spline
249 vfilein >> ai3D >> bi3D >> ci3D;
250 //setting spline coefficients for vertical electric field
251 unitIF=eV/cm;
252 fElectricFieldY->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
253 i, j, k);
254
255 //reading the coefficients of cubic spline
256 vfilein >> ai3D >> bi3D >> ci3D;
257 //setting spline coefficients for nuclear density (first element)
258 unitIF=1.;
259 fNucleiDensity[0]->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
260 i, j, k);
261
262 //reading the coefficients of cubic spline
263 vfilein >> ai3D >> bi3D >> ci3D;
264 //setting spline coefficients for electron density
265 unitIF=1./cm3;
266 fElectronDensity->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
267 i, j, k);
268
269 for(G4int ii=1; ii<fNelements; ii++)
270 {
271 //reading the coefficients of cubic spline
272 vfilein >> ai3D >> bi3D >> ci3D;
273 //setting spline coefficients for nuclear density (other elements if any)
274 unitIF=1.;
275 fNucleiDensity[ii]->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF,
276 ci3D*unitIF,
277 i, j, k);
278 }
279
280 }
281 }
282 }
283 }
284
285 vfilein.close();
286
287 //set special values and coefficients
288 G4double alphahbarc2=std::pow(CLHEP::fine_structure_const*CLHEP::hbarc ,2.);
289 fK30=2.*CLHEP::pi*alphahbarc2/CLHEP::electron_mass_c2;
290
291 for(G4int i=0; i<fNelements; i++)
292 {
293 fRF.push_back((std::pow(9*CLHEP::pi*CLHEP::pi/128/fZ1[i],1/3.))
294 *0.5291772109217*angstrom);//Thomas-Fermi screening radius
295
296 fTetamax0.push_back(CLHEP::hbarc/(fR0*std::pow(fAN[i],1./3.)));
297 fTeta10.push_back(CLHEP::hbarc/fRF[i]);
298 fPu11.push_back(std::pow(fU1[i]/CLHEP::hbarc,2.));
299
300 fK20.push_back(alphahbarc2*4*CLHEP::pi*fN0[i]*fZ1[i]*fZ1[i]);
301
302 fK40.push_back(3.76*std::pow(CLHEP::fine_structure_const*fZ1[i],2.));
303
304 fKD.push_back(fK30*fZ1[i]*fN0[i]);
305
306 fLogPlasmaEdI0.push_back(G4Log((crystal->GetIonisation()->GetPlasmaEnergy())/fI0[i]));
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)/fNsteps;//necessary to define simulation
319 //step = fChannelingStep =
320 // = fChannelingStep0*sqrt(pv)
321 }
322
323 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
324
325 G4ThreeVector G4ChannelingFastSimCrystalData::CoordinatesFromBoxToLattice
326 (const G4ThreeVector &pos0)
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(fBendingRsquare -
336 fBending2R*(x0*fCosMiscutAngle - z0*fSinMiscutAngle) +
337 x0*x0 + z0*z0);
338 //transform to co-rotating reference system connected with "central plane/axis"
339 x = fBendingR - rsqrt;
340 y = y0;
341 z = fBendingR*std::asin((z0*fCosMiscutAngle + x0*fSinMiscutAngle)/rsqrt);
342 }
343 else
344 {
345 //for straight crystal
346 x = x0*fCosMiscutAngle - z0*fSinMiscutAngle;
347 y = y0;
348 z = x0*fSinMiscutAngle + z0*fCosMiscutAngle;
349
350 //for crystalline undulator
351 if(fCU){x-=GetCUx(z);}
352 }
353
354 //calculation of coordinates within a channel (periodic cell)
355 fNChannelx=std::floor(x/fDx); //remember the horizontal channel number
356 //to track the particle
357 x-=fNChannelx*fDx;
358 fNChannely=std::floor(y/fDy);//remember the vertical channel number
359 //to track the particle (=0 for planar case)
360 y-=fNChannely*fDy;
361 //correction of the longitudinal coordinate
362 if (fBent) {fCorrectionZ = fBendingR/(fBendingR-fNChannelx*fDx);}
363
364 return G4ThreeVector(x,y,z);
365 }
366
367 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
368
369 G4ThreeVector G4ChannelingFastSimCrystalData::CoordinatesFromLatticeToBox(
370 const G4ThreeVector &pos)
371 {
372 G4double x=pos.x(),y=pos.y(),z=pos.z();
373
374 //transform to co-rotating reference system connected with "central plane/axis"
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. - std::cos(z/fBendingR));
384 G4double a = x + rcos;
385 G4double b = std::sqrt(x*x + fBending2R*rcos - a*a);
386
387 //transform to Box coordinates
388 x0 = a*fCosMiscutAngle + b*fSinMiscutAngle;
389 y0 = y;
390 z0 = b*fCosMiscutAngle - a*fSinMiscutAngle;
391 }
392 else
393 {
394 //for crystalline undulator
395 if(fCU){x+=GetCUx(z);}
396
397 //for straight crystal
398 x0 = x*fCosMiscutAngle + z*fSinMiscutAngle;
399 y0 = y;
400 z0 =-x*fSinMiscutAngle + z*fCosMiscutAngle;
401 }
402
403 return G4ThreeVector(x0,y0,z0-fHalfDimBoundingBox.z());
404 }
405
406 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
407
408 G4ThreeVector G4ChannelingFastSimCrystalData::ChannelChange(G4double& x,
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 coordinate
419 if (fBent) {fCorrectionZ = fBendingR/(fBendingR-fNChannelx*fDx);}
420 }
421 else if (x>=fDx)
422 {
423 fNChannelx+=1;
424 x-=fDx; //enter in other channel
425 //correction of the longitudinal coordinate
426 if (fBent) {fCorrectionZ = fBendingR/(fBendingR-fNChannelx*fDx);}
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........oooOO0OOooo........oooOO0OOooo....
445