| 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
27 #include "G4Channeling.hh"
28
29 #include "Randomize.hh"
30
31 #include "G4ChannelingTrackData.hh"
32 #include "G4TouchableHistory.hh"
33 #include "G4SystemOfUnits.hh"
34 #include "G4LambdacPlus.hh"
35 #include "G4PhysicsModelCatalog.hh"
36
37 G4Channeling::G4Channeling():
38 G4VDiscreteProcess("channeling"),
39 fChannelingID(G4PhysicsModelCatalog::GetModelI
40 fTimeStepMin(0.),
41 fTimeStepMax(0.),
42 fTransverseVariationMax(2.E-2 * CLHEP::angstro
43 k010(G4ThreeVector(0.,1.,0.)),
44 fSpin(G4ThreeVector(0.,0.,0.))
45 {}
46
47 //....oooOO0OOooo........oooOO0OOooo........oo
48
49 G4Channeling::~G4Channeling(){;}
50
51 //....oooOO0OOooo........oooOO0OOooo........oo
52
53 G4ChannelingTrackData* G4Channeling::GetTrackD
54 G4ChannelingTrackData* trackdata =
55 (G4ChannelingTrackData*)(aTrack.GetAux
56 if(trackdata == nullptr){
57 trackdata = new G4ChannelingTrackData(
58 aTrack.SetAuxiliaryTrackInformation(fC
59 }
60 return trackdata;
61 }
62
63 //....oooOO0OOooo........oooOO0OOooo........oo
64
65 void G4Channeling::GetEF(const G4Track& aTrack
66 G4ThreeVector& pos,
67 G4ThreeVector& out){
68 out = G4ThreeVector((GetMatData(aTrack)->G
69 (GetMatData(aTrack)->G
70 0.);
71 }
72
73 //....oooOO0OOooo........oooOO0OOooo........oo
74
75 void G4Channeling::PosToLattice(G4StepPoint* s
76 G4TouchableHistory* theTouchable = (G4Touc
77
78 pos -= theTouchable->GetTranslation();
79 pos = ((*theTouchable->GetRotation()).inve
80 }
81
82 //....oooOO0OOooo........oooOO0OOooo........oo
83
84 G4bool G4Channeling::UpdateParameters(const G4
85
86 G4LogicalCrystalVolume* aLCV = (G4LogicalC
87
88 G4StepPoint* postStepPoint = aTrack.GetSte
89 G4StepPoint* preStepPoint = aTrack.GetStep
90
91 G4ThreeVector posPost = postStepPoint->Get
92 aLCV->RotateToLattice(posPost);
93 G4ThreeVector posPre = preStepPoint->GetPo
94 aLCV->RotateToLattice(posPre);
95
96 G4double integrationLimit = std::fabs(posP
97
98 if(integrationLimit > 0.){
99 //------------------------------------
100 // Check if the crystal is bent
101 //------------------------------------
102 G4bool isBent = GetMatData(aTrack)->Is
103
104 //------------------------------------
105 // Get the momentum in the world refer
106 // frame and rotate to the solid refer
107 //------------------------------------
108 G4TouchableHistory* theTouchable = (G4
109 G4ThreeVector momWorld = aTrack.GetSte
110 G4ThreeVector mom = (*theTouchable->Ge
111
112 //------------------------------------
113 // Get the momentum in the solid refer
114 // frame and rotate to the crystal ref
115 //------------------------------------
116 aLCV->RotateToLattice(mom);
117
118 //------------------------------------
119 // Get the momentum in the crystal ref
120 // frame and rotate to the reference f
121 // solidal to the bent planes
122 //------------------------------------
123 if(isBent){
124 PosToLattice(preStepPoint,posPre);
125 G4ThreeVector axis010 = (*theTouch
126 mom.rotate(axis010,-posPre.z()/Get
127 }
128
129 //------------------------------------
130 // Take the position stored in the tra
131 // If the particle enters the crystal,
132 // the position in the channel is rand
133 // generated using a uniform distribut
134 //------------------------------------
135 G4ThreeVector pos;
136 if(GetTrackData(aTrack)->GetPosCh().x(
137 G4double posX = G4UniformRand() *
138 G4double posY = G4UniformRand() *
139 pos = G4ThreeVector(posX,posY,0.);
140 }
141 else{
142 pos = GetTrackData(aTrack)->GetPos
143 }
144
145 G4double step=0., stepTot=0.;
146 G4double nud =0., eld =0.;
147 G4double efx =0., efy =0.;
148 G4double nud_temp =0., eld_temp =0.
149
150 G4double beta = aTrack.GetVelocity()/C
151 G4double Z = GetParticleDefinition(aTr
152
153 const G4double oneSixth = 1./6.;
154 G4ThreeVector posk1,posk2,posk3,posk4,
155 G4ThreeVector momk1,momk2,momk3,momk4,
156 G4ThreeVector pos_temp, efxy;
157
158 do{
159 //--------------------------------
160 // Limit the variable step length
161 // integration via the selected al
162 // and update variables for the in
163 //--------------------------------
164
165 UpdateIntegrationStep(aTrack,mom,s
166 if(step + stepTot > integrationLim
167 step = integrationLimit - step
168 }
169
170 //--------------------------------
171 // Function integration algorithm
172 // 4th Order Runge-Kutta
173 //--------------------------------
174
175 GetEF(aTrack,pos,efxy);
176 posk1 = step / mom.z() * mom;
177 momk1 = step / beta * Z * efxy;
178 if(isBent) momk1.setX(momk1.x() -
179
180 GetEF(aTrack,pos_temp = pos + posk
181 posk2 = step / mom.z() * (mom + mo
182 momk2 = step / beta * Z * efxy;
183 if(isBent) momk2.setX(momk2.x() -
184
185 GetEF(aTrack,pos_temp = pos + posk
186 posk3 = step / mom.z() * (mom + mo
187 momk3 = step / beta * Z * efxy;
188 if(isBent) momk3.setX(momk3.x() -
189
190 GetEF(aTrack,pos_temp = pos + posk
191 posk4 = step / mom.z() * (mom + mo
192 momk4 = step / beta * Z * efxy;
193 if(isBent) momk4.setX(momk4.x() -
194
195 pos = pos + oneSixth * (posk1 + 2.
196 mom = mom + oneSixth * (momk1 + 2.
197
198 //--------------------------------
199 // Function integration algorithm
200 // 2th Order Velocity-Verlet
201 //--------------------------------
202
203 /*
204 GetEF(aTrack,pos,efxy);
205 posk1 = pos + (step * 0.5 / mom.z(
206 //momk1 = mom + step * 0.5 / betaZ
207 momk1 = mom;
208 if(isBent) momk1.setX(momk1.x() -
209
210 GetEF(aTrack,posk1,efxy);
211 pos = pos + (step / momk1.z()) * m
212 //mom = mom + step / betaZ * efxy;
213 mom = mom;
214 if(isBent) mom.setX(mom.x() - step
215 */
216
217 //--------------------------------
218 // Update the total step and the e
219 // and nuclei density experienced
220 // the particle during its motion
221 //--------------------------------
222
223 stepTot += step;
224
225 nud_temp = GetMatData(aTrack)->Get
226 eld_temp = GetMatData(aTrack)->Get
227
228 if(nud_temp < 0.) {nud_temp = 0.;}
229 if(eld_temp < 0.) {eld_temp = 0.;}
230
231 nud += (step * nud_temp);
232 eld += (step * eld_temp);
233
234 efx += (step * GetMatData(aTrack)-
235 efy += (step * GetMatData(aTrack)-
236
237 } while(stepTot<integrationLimit);
238
239 nud /= stepTot;
240 eld /= stepTot;
241
242 if(nud < 1.E-10) {nud = 1.E-10;}
243 if(eld < 1.E-10) {eld = 1.E-10;}
244
245 GetTrackData(aTrack)->SetNuD(nud);
246 GetTrackData(aTrack)->SetElD(eld);
247
248 GetTrackData(aTrack)->SetEFX(efx);
249 GetTrackData(aTrack)->SetEFY(efy);
250
251 GetTrackData(aTrack)->SetMomCh(mom);
252 GetTrackData(aTrack)->SetPosCh(pos);
253 return true;
254 }
255 else{
256 return false;
257 }
258
259 return false;
260 }
261
262 //....oooOO0OOooo........oooOO0OOooo........oo
263
264 G4bool G4Channeling::
265 UpdateIntegrationStep(const G4Track& aTrack,
266 G4ThreeVector& mom,
267 G4double& step){
268
269 if(mom.x() != 0.0 || mom.y() != 0.0){
270 double xy2 = mom.x() * mom.x() + mom.y
271
272 if(xy2!=0.){
273 step = std::fabs(fTransverseVariat
274 if(step < fTimeStepMin) step = fTi
275 else{
276 fTimeStepMax = std::sqrt( fTra
277 / std::fab
278
279 if(step > fTimeStepMax) step =
280 }
281 }
282 else{
283 step = fTimeStepMin;
284 }
285
286 return true;
287 }
288 else{
289 step = fTimeStepMin;
290 }
291 return false;
292 }
293
294 //....oooOO0OOooo........oooOO0OOooo........oo
295
296 G4double G4Channeling::
297 GetMeanFreePath(const G4Track& aTrack,
298 G4double, // previousStepSize
299 G4ForceCondition* condition){
300
301 //----------------------------------------
302 // the condition is forced to check if
303 // the volume has a lattice at each step.
304 // if it hasn't, return DBL_MAX
305 //----------------------------------------
306
307 *condition = Forced;
308
309 G4LogicalVolume* aLV = aTrack.GetVolume()-
310 G4LogicalVolume* aNLV = aTrack.GetNextVolu
311
312 if(G4LogicalCrystalVolume::IsLattice(aLV)
313 G4LogicalCrystalVolume::IsLattice(aNLV)
314 G4double osc_per = GetOscillationPerio
315 fTimeStepMin = osc_per * 2.E-4;
316 return osc_per * 0.01;
317 }
318 else{
319 GetTrackData(aTrack)->Reset();
320 return DBL_MAX;
321 }
322 }
323
324 //....oooOO0OOooo........oooOO0OOooo........oo
325
326 G4VParticleChange* G4Channeling::
327 PostStepDoIt(const G4Track& aTrack,
328 const G4Step&){
329
330 //----------------------------------------
331 // check if the volume has a lattice
332 // and if the particle is in channeling.
333 // If it is so, the particle is forced
334 // to follow the channeling plane
335 // direction. If the particle has
336 // dechanneled or exited the crystal,
337 // the outgoing angle is evaluated
338 //----------------------------------------
339
340 aParticleChange.Initialize(aTrack);
341 G4LogicalVolume* aLV = aTrack.GetVolume()-
342 G4LogicalVolume* aNLV = aTrack.GetNextVolu
343
344
345 if(G4LogicalCrystalVolume::IsLattice(aLV)
346 G4LogicalCrystalVolume::IsLattice(aNLV)
347
348 G4bool bModifiedTraj = UpdateParameter
349
350 if(bModifiedTraj==true){
351 //--------------------------------
352 // Get the momentum in the referen
353 // solidal to the bent planes and
354 // to the reference frame
355 //--------------------------------
356 G4LogicalCrystalVolume* aLCV = (G4
357 G4ThreeVector momCh = GetTrackData
358
359 G4StepPoint* postStepPoint = aTrac
360 G4TouchableHistory* theTouchable =
361
362 if(GetMatData(aTrack)->IsBent()){
363 G4ThreeVector posPost = postSt
364 PosToLattice(postStepPoint,pos
365 G4ThreeVector axis010 = (*theT
366 momCh.rotate(axis010,posPost.z
367 }
368
369 //--------------------------------
370 // Get the momentum in the crystal
371 // frame and rotate to the solid r
372 //--------------------------------
373
374 aLCV->RotateToSolid(momCh);
375
376 //--------------------------------
377 // Get the momentum in the solid r
378 // frame and rotate to the world r
379 //--------------------------------
380 G4ThreeVector mom = ((*theTouchabl
381
382 aParticleChange.ProposeMomentumDir
383 aParticleChange.ProposePolarizatio
384 }
385 }
386 else{
387 // if the volume has no lattice it res
388 GetTrackData(aTrack)->Reset();
389 }
390
391 return &aParticleChange;
392 }
393
394 //....oooOO0OOooo........oooOO0OOooo........oo
395