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