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