| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Author: Alexei Sytov <<
27 // Co-author: Gianfranco Paternò (modificat <<
28 26
29 #include "G4ChannelingFastSimCrystalData.hh" 27 #include "G4ChannelingFastSimCrystalData.hh"
30 #include "G4SystemOfUnits.hh" 28 #include "G4SystemOfUnits.hh"
31 #include "G4PhysicalConstants.hh" 29 #include "G4PhysicalConstants.hh"
32 30
33 G4ChannelingFastSimCrystalData::G4ChannelingFa 31 G4ChannelingFastSimCrystalData::G4ChannelingFastSimCrystalData()
34 { 32 {
35 33
36 } 34 }
37 35
38 //....oooOO0OOooo........oooOO0OOooo........oo 36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39 37
40 void G4ChannelingFastSimCrystalData::SetMateri << 38 void G4ChannelingFastSimCrystalData::SetMaterialProperties(const G4Material *crystal,
41 const G4M << 39 const G4String &lattice)
42 const G4S <<
43 const G4S <<
44 { 40 {
45 G4String filename=crystal->GetName(); //in 41 G4String filename=crystal->GetName(); //input file
46 filename.erase(0,3); 42 filename.erase(0,3);
47 43
48 if (fVerbosity) 44 if (fVerbosity)
49 { 45 {
50 G4cout << 46 G4cout <<
51 "========================================= 47 "======================================================================="
52 << G4endl; 48 << G4endl;
53 G4cout << << 49 G4cout <<
54 "====== Crystal lattice << 50 "====== Crystal lattice data ========"
55 << G4endl; << 51 << G4endl;
56 G4cout << 52 G4cout <<
57 "========================================= 53 "======================================================================="
58 << G4endl; 54 << G4endl;
59 G4cout << "Crystal material: " << filename < 55 G4cout << "Crystal material: " << filename << G4endl;
60 } 56 }
61 57
62 //choice between planes (1D model) and axe 58 //choice between planes (1D model) and axes (2D model)
63 if (lattice.compare(0,1,"(")==0) 59 if (lattice.compare(0,1,"(")==0)
64 { 60 {
65 iModel=1; //planes 61 iModel=1; //planes
66 filename = filename + "_planes_"; //temp 62 filename = filename + "_planes_"; //temporary name
67 if (fVerbosity) 63 if (fVerbosity)
68 G4cout << "Crystal planes: " << lattice << G 64 G4cout << "Crystal planes: " << lattice << G4endl;
69 } 65 }
70 else if (lattice.compare(0,1,"<")==0) 66 else if (lattice.compare(0,1,"<")==0)
71 { 67 {
72 iModel=2; //axes 68 iModel=2; //axes
73 filename = filename + "_axes_"; //tempor 69 filename = filename + "_axes_"; //temporary name
74 if (fVerbosity) 70 if (fVerbosity)
75 G4cout << "Crystal axes: " << lattice << G4e 71 G4cout << "Crystal axes: " << lattice << G4endl;
76 } 72 }
77 73
78 //input file: 74 //input file:
79 filename = filename + lattice.substr(1,(la 75 filename = filename + lattice.substr(1,(lattice.length())-2) + ".dat";
80 76
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 77 fNelements=(G4int)crystal->GetNumberOfElements();
94 for(G4int i=0; i<fNelements; i++) 78 for(G4int i=0; i<fNelements; i++)
95 { 79 {
96 fZ1.push_back(crystal->GetElement(i)->Ge 80 fZ1.push_back(crystal->GetElement(i)->GetZ());
97 fAN.push_back(crystal->GetElement(i)->Ge 81 fAN.push_back(crystal->GetElement(i)->GetAtomicMassAmu());
98 fI0.push_back(crystal->GetElement(i)->Ge 82 fI0.push_back(crystal->GetElement(i)->GetIonisation()->GetMeanExcitationEnergy());
99 } 83 }
100 84
101 G4double var;//just variable 85 G4double var;//just variable
102 G4double unitIF;//unit of interpolation fu 86 G4double unitIF;//unit of interpolation function
103 87
104 std::ifstream vfilein; 88 std::ifstream vfilein;
105 vfilein.open(filename); 89 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 90
118 //read nuclear concentration 91 //read nuclear concentration
119 for(G4int i=0; i<fNelements; i++) 92 for(G4int i=0; i<fNelements; i++)
120 { 93 {
121 vfilein >> var; 94 vfilein >> var;
122 fN0.push_back(var/cm3); 95 fN0.push_back(var/cm3);
123 } 96 }
124 97
125 //read amplitude of thermal oscillations 98 //read amplitude of thermal oscillations
126 for(G4int i=0; i<fNelements; i++) 99 for(G4int i=0; i<fNelements; i++)
127 { 100 {
128 vfilein >> var; 101 vfilein >> var;
129 fU1.push_back(var*cm); 102 fU1.push_back(var*cm);
130 } 103 }
131 104
132 if (iModel==1) 105 if (iModel==1)
133 { 106 {
134 // read channel dimensions 107 // read channel dimensions
135 vfilein >> fDx; 108 vfilein >> fDx;
136 fDx*=cm; 109 fDx*=cm;
137 // read interpolation step size 110 // read interpolation step size
138 vfilein >> fNpointsx; 111 vfilein >> fNpointsx;
139 112
140 fDy = fDx; 113 fDy = fDx;
141 fNpointsy = 0; 114 fNpointsy = 0;
142 } 115 }
143 else if (iModel==2) 116 else if (iModel==2)
144 { 117 {
145 // read channel dimensions 118 // read channel dimensions
146 vfilein >> fDx >> fDy; 119 vfilein >> fDx >> fDy;
147 fDx*=cm; 120 fDx*=cm;
148 fDy*=cm; 121 fDy*=cm;
149 // read the number of nodes of interpol 122 // read the number of nodes of interpolation
150 vfilein >> fNpointsx >> fNpointsy; 123 vfilein >> fNpointsx >> fNpointsy;
151 } 124 }
152 125
153 //read the height of the potential well, n 126 //read the height of the potential well, necessary only for step length calculation
154 vfilein >> fVmax; 127 vfilein >> fVmax;
155 fVmax*=eV; 128 fVmax*=eV;
156 fVmax2=2.*fVmax; 129 fVmax2=2.*fVmax;
157 130
158 //read the on-zero minimal potential insid 131 //read the on-zero minimal potential inside the crystal,
159 //necessary for angle recalculation for en 132 //necessary for angle recalculation for entrance/exit through
160 //the crystal lateral surface 133 //the crystal lateral surface
161 vfilein >> fVMinCrystal; 134 vfilein >> fVMinCrystal;
162 fVMinCrystal*=eV; 135 fVMinCrystal*=eV;
163 136
164 // to create the class of interpolation fo 137 // to create the class of interpolation for any function
165 fElectricFieldX = 138 fElectricFieldX =
166 new G4ChannelingFastSimInterpolati 139 new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);
167 if(iModel==2) {fElectricFieldY = 140 if(iModel==2) {fElectricFieldY =
168 new G4ChannelingFastSimInterpolati 141 new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);}
169 fElectronDensity = 142 fElectronDensity =
170 new G4ChannelingFastSimInterpolati 143 new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);
171 fMinIonizationEnergy = 144 fMinIonizationEnergy =
172 new G4ChannelingFastSimInterpolati 145 new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel);
173 146
174 // do it for any element of crystal materi 147 // do it for any element of crystal material
175 for(G4int i=0; i<fNelements; i++) 148 for(G4int i=0; i<fNelements; i++)
176 { 149 {
177 fNucleiDensity.push_back( 150 fNucleiDensity.push_back(
178 new G4ChannelingFastSimInterpolati 151 new G4ChannelingFastSimInterpolation(fDx,fDy,fNpointsx,fNpointsy,iModel));
179 } 152 }
180 153
181 if (iModel==1) 154 if (iModel==1)
182 { 155 {
183 G4double ai, bi, ci, di; 156 G4double ai, bi, ci, di;
184 for(G4int i=0; i<fNpointsx; i++) 157 for(G4int i=0; i<fNpointsx; i++)
185 { 158 {
186 //reading the coefficients of cubic sp 159 //reading the coefficients of cubic spline
187 vfilein >> ai >> bi >> ci >> di; 160 vfilein >> ai >> bi >> ci >> di;
188 //setting spline coefficients for elec 161 //setting spline coefficients for electric field
189 unitIF=eV/cm; 162 unitIF=eV/cm;
190 fElectricFieldX->SetCoefficients1D(ai* 163 fElectricFieldX->SetCoefficients1D(ai*unitIF, bi*unitIF,
191 ci* 164 ci*unitIF, di*unitIF, i);
192 165
193 //reading the coefficients of cubic sp 166 //reading the coefficients of cubic spline
194 vfilein >> ai >> bi >> ci >> di; 167 vfilein >> ai >> bi >> ci >> di;
195 //setting spline coefficients for nucl 168 //setting spline coefficients for nuclear density (first element)
196 unitIF=1.; 169 unitIF=1.;
197 fNucleiDensity[0]->SetCoefficients1D(a 170 fNucleiDensity[0]->SetCoefficients1D(ai*unitIF, bi*unitIF,
198 c 171 ci*unitIF, di*unitIF, i);
199 172
200 //reading the coefficients of cubic sp 173 //reading the coefficients of cubic spline
201 vfilein >> ai >> bi >> ci >> di; 174 vfilein >> ai >> bi >> ci >> di;
202 //setting spline coefficients for elec 175 //setting spline coefficients for electron density
203 unitIF=1./cm3; 176 unitIF=1./cm3;
204 fElectronDensity->SetCoefficients1D(ai 177 fElectronDensity->SetCoefficients1D(ai*unitIF, bi*unitIF,
205 ci 178 ci*unitIF, di*unitIF, i);
206 179
207 //reading the coefficients of cubic sp 180 //reading the coefficients of cubic spline
208 vfilein >> ai >> bi >> ci >> di; 181 vfilein >> ai >> bi >> ci >> di;
209 //setting spline coefficients for mini 182 //setting spline coefficients for minimal ionization energy
210 unitIF=eV; 183 unitIF=eV;
211 fMinIonizationEnergy->SetCoefficients1 184 fMinIonizationEnergy->SetCoefficients1D(ai*unitIF, bi*unitIF,
212 185 ci*unitIF, di*unitIF, i);
213 186
214 for(G4int ii=1; ii<fNelements; ii++) 187 for(G4int ii=1; ii<fNelements; ii++)
215 { 188 {
216 //reading the coefficients of cubi 189 //reading the coefficients of cubic spline
217 vfilein >> ai >> bi >> ci >> di; 190 vfilein >> ai >> bi >> ci >> di;
218 //setting spline coefficients for 191 //setting spline coefficients for nuclear density (other elements if any)
219 unitIF=1.; 192 unitIF=1.;
220 fNucleiDensity[ii]->SetCoefficient 193 fNucleiDensity[ii]->SetCoefficients1D(ai*unitIF, bi*unitIF,
221 194 ci*unitIF, di*unitIF, i);
222 } 195 }
223 } 196 }
224 } 197 }
225 else if (iModel==2) 198 else if (iModel==2)
226 { 199 {
227 G4double ai3D, bi3D, ci3D; 200 G4double ai3D, bi3D, ci3D;
228 for(G4int j=0; j<fNpointsy; j++) 201 for(G4int j=0; j<fNpointsy; j++)
229 { 202 {
230 for(G4int i=0; i<fNpointsx+1; i++) 203 for(G4int i=0; i<fNpointsx+1; i++)
231 { 204 {
232 for(G4int k=0; k<2; k++) 205 for(G4int k=0; k<2; k++)
233 { 206 {
234 //reading the coefficients of cubi 207 //reading the coefficients of cubic spline
235 vfilein >> ai3D >> bi3D >> ci3D; 208 vfilein >> ai3D >> bi3D >> ci3D;
236 unitIF=eV; 209 unitIF=eV;
237 //setting spline coefficients for 210 //setting spline coefficients for minimal ionization energy
238 fMinIonizationEnergy->SetCoefficie 211 fMinIonizationEnergy->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
239 212 i, j, k);
240 213
241 //reading the coefficients of cubi 214 //reading the coefficients of cubic spline
242 vfilein >> ai3D >> bi3D >> ci3D; 215 vfilein >> ai3D >> bi3D >> ci3D;
243 //setting spline coefficients for 216 //setting spline coefficients for horizontal electric field
244 unitIF=eV/cm; 217 unitIF=eV/cm;
245 fElectricFieldX->SetCoefficients2D 218 fElectricFieldX->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
246 219 i, j, k);
247 220
248 //reading the coefficients of cubi 221 //reading the coefficients of cubic spline
249 vfilein >> ai3D >> bi3D >> ci3D; 222 vfilein >> ai3D >> bi3D >> ci3D;
250 //setting spline coefficients for 223 //setting spline coefficients for vertical electric field
251 unitIF=eV/cm; 224 unitIF=eV/cm;
252 fElectricFieldY->SetCoefficients2D 225 fElectricFieldY->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
253 226 i, j, k);
254 227
255 //reading the coefficients of cubi 228 //reading the coefficients of cubic spline
256 vfilein >> ai3D >> bi3D >> ci3D; 229 vfilein >> ai3D >> bi3D >> ci3D;
257 //setting spline coefficients for 230 //setting spline coefficients for nuclear density (first element)
258 unitIF=1.; 231 unitIF=1.;
259 fNucleiDensity[0]->SetCoefficients 232 fNucleiDensity[0]->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
260 233 i, j, k);
261 234
262 //reading the coefficients of cubi 235 //reading the coefficients of cubic spline
263 vfilein >> ai3D >> bi3D >> ci3D; 236 vfilein >> ai3D >> bi3D >> ci3D;
264 //setting spline coefficients for 237 //setting spline coefficients for electron density
265 unitIF=1./cm3; 238 unitIF=1./cm3;
266 fElectronDensity->SetCoefficients2 239 fElectronDensity->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF, ci3D*unitIF,
267 240 i, j, k);
268 241
269 for(G4int ii=1; ii<fNelements; ii+ 242 for(G4int ii=1; ii<fNelements; ii++)
270 { 243 {
271 //reading the coefficients of 244 //reading the coefficients of cubic spline
272 vfilein >> ai3D >> bi3D >> ci3 245 vfilein >> ai3D >> bi3D >> ci3D;
273 //setting spline coefficients 246 //setting spline coefficients for nuclear density (other elements if any)
274 unitIF=1.; 247 unitIF=1.;
275 fNucleiDensity[ii]->SetCoeffic 248 fNucleiDensity[ii]->SetCoefficients2D(ai3D*unitIF, bi3D*unitIF,
276 249 ci3D*unitIF,
277 250 i, j, k);
278 } 251 }
279 252
280 } 253 }
281 } 254 }
282 } 255 }
283 } 256 }
284 257
285 vfilein.close(); 258 vfilein.close();
286 259
287 //set special values and coefficients 260 //set special values and coefficients
288 G4double alphahbarc2=std::pow(CLHEP::fine_ 261 G4double alphahbarc2=std::pow(CLHEP::fine_structure_const*CLHEP::hbarc ,2.);
289 fK30=2.*CLHEP::pi*alphahbarc2/CLHEP::elect 262 fK30=2.*CLHEP::pi*alphahbarc2/CLHEP::electron_mass_c2;
290 263
291 for(G4int i=0; i<fNelements; i++) 264 for(G4int i=0; i<fNelements; i++)
292 { 265 {
293 fRF.push_back((std::pow(9*CLHEP::pi*CLHE 266 fRF.push_back((std::pow(9*CLHEP::pi*CLHEP::pi/128/fZ1[i],1/3.))
294 *0.5291772109217*angstrom) 267 *0.5291772109217*angstrom);//Thomas-Fermi screening radius
295 268
296 fTetamax0.push_back(CLHEP::hbarc/(fR0*st 269 fTetamax0.push_back(CLHEP::hbarc/(fR0*std::pow(fAN[i],1./3.)));
297 fTeta10.push_back(CLHEP::hbarc/fRF[i]); 270 fTeta10.push_back(CLHEP::hbarc/fRF[i]);
298 fPu11.push_back(std::pow(fU1[i]/CLHEP::h 271 fPu11.push_back(std::pow(fU1[i]/CLHEP::hbarc,2.));
299 272
300 fK20.push_back(alphahbarc2*4*CLHEP::pi*f 273 fK20.push_back(alphahbarc2*4*CLHEP::pi*fN0[i]*fZ1[i]*fZ1[i]);
301 274
302 fK40.push_back(3.76*std::pow(CLHEP::fine 275 fK40.push_back(3.76*std::pow(CLHEP::fine_structure_const*fZ1[i],2.));
303 276
304 fKD.push_back(fK30*fZ1[i]*fN0[i]); 277 fKD.push_back(fK30*fZ1[i]*fN0[i]);
305 <<
306 fLogPlasmaEdI0.push_back(G4Log((crystal- <<
307 } 278 }
308 279
309 fBB.resize(fNelements); 280 fBB.resize(fNelements);
310 fE1XBbb.resize(fNelements); 281 fE1XBbb.resize(fNelements);
311 fBBDEXP.resize(fNelements); 282 fBBDEXP.resize(fNelements);
312 fPzu11.resize(fNelements); 283 fPzu11.resize(fNelements);
313 fTeta12.resize(fNelements); 284 fTeta12.resize(fNelements);
314 fTetamax2.resize(fNelements); 285 fTetamax2.resize(fNelements);
315 fTetamax12.resize(fNelements); 286 fTetamax12.resize(fNelements);
316 fK2.resize(fNelements); 287 fK2.resize(fNelements);
317 288
318 fChangeStep = CLHEP::pi*std::min(fDx,fDy)/ 289 fChangeStep = CLHEP::pi*std::min(fDx,fDy)/fNsteps;//necessary to define simulation
319 //ste 290 //step = fChannelingStep =
320 // = 291 // = fChannelingStep0*sqrt(pv)
321 } 292 }
322 293
323 //....oooOO0OOooo........oooOO0OOooo........oo 294 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
324 295
325 G4ThreeVector G4ChannelingFastSimCrystalData:: 296 G4ThreeVector G4ChannelingFastSimCrystalData::CoordinatesFromBoxToLattice
326 297 (const G4ThreeVector &pos0)
327 { 298 {
328 G4double x0=pos0.x(),y0=pos0.y(),z0=pos0.z( 299 G4double x0=pos0.x(),y0=pos0.y(),z0=pos0.z();
329 z0+=fHalfDimBoundingBox.z(); 300 z0+=fHalfDimBoundingBox.z();
330 G4double x,y,z; 301 G4double x,y,z;
331 302
332 if (fBent) 303 if (fBent)
333 { 304 {
334 // for bent crystal 305 // for bent crystal
335 G4double rsqrt = std::sqrt(fBendingRsqu 306 G4double rsqrt = std::sqrt(fBendingRsquare -
336 fBending2R*( 307 fBending2R*(x0*fCosMiscutAngle - z0*fSinMiscutAngle) +
337 x0*x0 + z0*z 308 x0*x0 + z0*z0);
338 //transform to co-rotating reference sy 309 //transform to co-rotating reference system connected with "central plane/axis"
339 x = fBendingR - rsqrt; 310 x = fBendingR - rsqrt;
340 y = y0; 311 y = y0;
341 z = fBendingR*std::asin((z0*fCosMiscutA 312 z = fBendingR*std::asin((z0*fCosMiscutAngle + x0*fSinMiscutAngle)/rsqrt);
342 } 313 }
343 else 314 else
344 { 315 {
345 //for straight crystal 316 //for straight crystal
346 x = x0*fCosMiscutAngle - z0*fSinMiscutA 317 x = x0*fCosMiscutAngle - z0*fSinMiscutAngle;
347 y = y0; 318 y = y0;
348 z = x0*fSinMiscutAngle + z0*fCosMiscutA 319 z = x0*fSinMiscutAngle + z0*fCosMiscutAngle;
349 320
350 //for crystalline undulator 321 //for crystalline undulator
351 if(fCU){x-=GetCUx(z);} 322 if(fCU){x-=GetCUx(z);}
352 } 323 }
353 324
354 //calculation of coordinates within a chann 325 //calculation of coordinates within a channel (periodic cell)
355 fNChannelx=std::floor(x/fDx); //remember th 326 fNChannelx=std::floor(x/fDx); //remember the horizontal channel number
356 //to track th 327 //to track the particle
357 x-=fNChannelx*fDx; 328 x-=fNChannelx*fDx;
358 fNChannely=std::floor(y/fDy);//remember the 329 fNChannely=std::floor(y/fDy);//remember the vertical channel number
359 //to track the 330 //to track the particle (=0 for planar case)
360 y-=fNChannely*fDy; 331 y-=fNChannely*fDy;
361 //correction of the longitudinal coordinate 332 //correction of the longitudinal coordinate
362 if (fBent) {fCorrectionZ = fBendingR/(fBend 333 if (fBent) {fCorrectionZ = fBendingR/(fBendingR-fNChannelx*fDx);}
363 334
364 return G4ThreeVector(x,y,z); 335 return G4ThreeVector(x,y,z);
365 } 336 }
366 337
367 //....oooOO0OOooo........oooOO0OOooo........oo 338 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
368 339
369 G4ThreeVector G4ChannelingFastSimCrystalData:: 340 G4ThreeVector G4ChannelingFastSimCrystalData::CoordinatesFromLatticeToBox(
370 const G4ThreeVector &pos) 341 const G4ThreeVector &pos)
371 { 342 {
372 G4double x=pos.x(),y=pos.y(),z=pos.z(); 343 G4double x=pos.x(),y=pos.y(),z=pos.z();
373 344
374 //transform to co-rotating reference system 345 //transform to co-rotating reference system connected with "central plane/axis"
375 x+=fNChannelx*fDx; 346 x+=fNChannelx*fDx;
376 y+=fNChannely*fDy; 347 y+=fNChannely*fDy;
377 348
378 G4double x0,y0,z0; 349 G4double x0,y0,z0;
379 350
380 if (fBent) 351 if (fBent)
381 { 352 {
382 // for bent crystal 353 // for bent crystal
383 G4double rcos = (fBendingR - x)*(1. - s 354 G4double rcos = (fBendingR - x)*(1. - std::cos(z/fBendingR));
384 G4double a = x + rcos; 355 G4double a = x + rcos;
385 G4double b = std::sqrt(x*x + fBending2R 356 G4double b = std::sqrt(x*x + fBending2R*rcos - a*a);
386 357
387 //transform to Box coordinates 358 //transform to Box coordinates
388 x0 = a*fCosMiscutAngle + b*fSinMiscutAn 359 x0 = a*fCosMiscutAngle + b*fSinMiscutAngle;
389 y0 = y; 360 y0 = y;
390 z0 = b*fCosMiscutAngle - a*fSinMiscutAn 361 z0 = b*fCosMiscutAngle - a*fSinMiscutAngle;
391 } 362 }
392 else 363 else
393 { 364 {
394 //for crystalline undulator 365 //for crystalline undulator
395 if(fCU){x+=GetCUx(z);} 366 if(fCU){x+=GetCUx(z);}
396 367
397 //for straight crystal 368 //for straight crystal
398 x0 = x*fCosMiscutAngle + z*fSinMiscutAn 369 x0 = x*fCosMiscutAngle + z*fSinMiscutAngle;
399 y0 = y; 370 y0 = y;
400 z0 =-x*fSinMiscutAngle + z*fCosMiscutAn 371 z0 =-x*fSinMiscutAngle + z*fCosMiscutAngle;
401 } 372 }
402 373
403 return G4ThreeVector(x0,y0,z0-fHalfDimBound 374 return G4ThreeVector(x0,y0,z0-fHalfDimBoundingBox.z());
404 } 375 }
405 376
406 //....oooOO0OOooo........oooOO0OOooo........oo 377 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
407 378
408 G4ThreeVector G4ChannelingFastSimCrystalData:: 379 G4ThreeVector G4ChannelingFastSimCrystalData::ChannelChange(G4double& x,
409 G4double& y, 380 G4double& y,
410 G4double& z) 381 G4double& z)
411 { 382 {
412 383
413 //test of enter in other channel 384 //test of enter in other channel
414 if (x<0) 385 if (x<0)
415 { 386 {
416 fNChannelx-=1; 387 fNChannelx-=1;
417 x+=fDx; //enter in other channel 388 x+=fDx; //enter in other channel
418 //correction of the longitudinal coord 389 //correction of the longitudinal coordinate
419 if (fBent) {fCorrectionZ = fBendingR/( 390 if (fBent) {fCorrectionZ = fBendingR/(fBendingR-fNChannelx*fDx);}
420 } 391 }
421 else if (x>=fDx) 392 else if (x>=fDx)
422 { 393 {
423 fNChannelx+=1; 394 fNChannelx+=1;
424 x-=fDx; //enter in other channel 395 x-=fDx; //enter in other channel
425 //correction of the longitudinal coord 396 //correction of the longitudinal coordinate
426 if (fBent) {fCorrectionZ = fBendingR/( 397 if (fBent) {fCorrectionZ = fBendingR/(fBendingR-fNChannelx*fDx);}
427 } 398 }
428 399
429 //test of enter in other channel 400 //test of enter in other channel
430 if (y<0) 401 if (y<0)
431 { 402 {
432 fNChannely-=1; 403 fNChannely-=1;
433 y+=fDy; //enter in other channel 404 y+=fDy; //enter in other channel
434 } 405 }
435 else if (y>=fDy) 406 else if (y>=fDy)
436 { 407 {
437 fNChannely+=1; 408 fNChannely+=1;
438 y-=fDy; //enter in other channel 409 y-=fDy; //enter in other channel
439 } 410 }
440 411
441 return G4ThreeVector(x,y,z); 412 return G4ThreeVector(x,y,z);
442 } 413 }
443 414
444 //....oooOO0OOooo........oooOO0OOooo........oo 415 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
445 416