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 // G4MicroElecSurface.cc, 27 // G4MicroElecSurface.cc, 28 // 2020/05/20 P. Caron, C. Ingui << 28 // 2020/05/20 P. Caron, C. Inguimbert are with ONERA [b] 29 // Q. Gibaru is with << 29 // Q. Gibaru is with CEA [a], ONERA [b] and CNES [c] 30 // M. Raine and D. La << 30 // M. Raine and D. Lambert are with CEA [a] 31 // 31 // 32 // A part of this work has been funded by the 32 // A part of this work has been funded by the French space agency(CNES[c]) 33 // [a] CEA, DAM, DIF - 91297 ARPAJON, France 33 // [a] CEA, DAM, DIF - 91297 ARPAJON, France 34 // [b] ONERA - DPHY, 2 avenue E.Belin, 31055 T 34 // [b] ONERA - DPHY, 2 avenue E.Belin, 31055 Toulouse, France 35 // [c] CNES, 18 av.E.Belin, 31401 Toulouse CED 35 // [c] CNES, 18 av.E.Belin, 31401 Toulouse CEDEX, France 36 // 36 // 37 // Based on the following publications 37 // Based on the following publications 38 // 38 // 39 // - Q.Gibaru, C.Inguimbert, P.Caron, M.Raine, << 39 // - Q.Gibaru, C.Inguimbert, P.Caron, M.Raine, D.Lambert, J.Puech, 40 // Geant4 physics processes for microdosimet << 40 // Geant4 physics processes for microdosimetry and secondary electron emission simulation : 41 // Extension of MicroElec to very low energi << 41 // Extension of MicroElec to very low energies and new materials 42 // NIM B, 2020, in review. << 42 // NIM B, 2020, in review. 43 // 43 // 44 // - Modele de transport d'electrons a basse e << 44 // 45 // applications spatiales (OSMOSEE, GEANT4), << 45 // - Modèle de transport d'électrons à basse énergie (10 eV- 2 keV) pour >> 46 // applications spatiales (OSMOSEE, GEANT4), PhD dissertation, 2017. >> 47 // 46 // 48 // 47 ////////////////////////////////////////////// 49 //////////////////////////////////////////////////////////////////////// 48 50 49 #include "G4MicroElecSurface.hh" 51 #include "G4MicroElecSurface.hh" >> 52 50 #include "G4ios.hh" 53 #include "G4ios.hh" 51 #include "G4PhysicalConstants.hh" 54 #include "G4PhysicalConstants.hh" 52 #include "G4EmProcessSubType.hh" 55 #include "G4EmProcessSubType.hh" 53 #include "G4GeometryTolerance.hh" 56 #include "G4GeometryTolerance.hh" >> 57 54 #include "G4SystemOfUnits.hh" 58 #include "G4SystemOfUnits.hh" 55 59 56 60 57 //....oooOO0OOooo........oooOO0OOooo........oo 61 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 58 62 59 G4MicroElecSurface::G4MicroElecSurface(const G 63 G4MicroElecSurface::G4MicroElecSurface(const G4String& processName,G4ProcessType type) 60 : G4VDiscreteProcess(processName, type), 64 : G4VDiscreteProcess(processName, type), 61 oldMomentum(0.,0.,0.), previousMomentum(0. 65 oldMomentum(0.,0.,0.), previousMomentum(0.,0.,0.), 62 theGlobalNormal(0.,0.,0.), theFacetNormal( 66 theGlobalNormal(0.,0.,0.), theFacetNormal(0.,0.,0.) 63 { 67 { 64 if ( verboseLevel > 0) 68 if ( verboseLevel > 0) 65 { << 69 { 66 G4cout << GetProcessName() << " is created << 70 G4cout << GetProcessName() << " is created " << G4endl; 67 } << 71 } 68 72 69 isInitialised=false; 73 isInitialised=false; 70 SetProcessSubType(fSurfaceReflection); << 74 SetProcessSubType(25); 71 75 72 theStatus = UndefinedSurf; 76 theStatus = UndefinedSurf; 73 material1 = nullptr; 77 material1 = nullptr; 74 material2 = nullptr; 78 material2 = nullptr; 75 79 76 kCarTolerance = G4GeometryTolerance::GetInst 80 kCarTolerance = G4GeometryTolerance::GetInstance()->GetSurfaceTolerance(); 77 theParticleMomentum = 0.; 81 theParticleMomentum = 0.; 78 82 79 flag_franchissement_surface = false; 83 flag_franchissement_surface = false; 80 flag_normal = false; 84 flag_normal = false; 81 flag_reflexion = false; 85 flag_reflexion = false; 82 teleportToDo = teleportDone = false; 86 teleportToDo = teleportDone = false; 83 87 84 ekint = thetat = thetaft = energyThreshold = 88 ekint = thetat = thetaft = energyThreshold = crossingProbability = 0.0; 85 } 89 } 86 90 87 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 88 92 89 G4MicroElecSurface::~G4MicroElecSurface() 93 G4MicroElecSurface::~G4MicroElecSurface() 90 {} 94 {} 91 95 92 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 93 97 94 G4bool 98 G4bool 95 G4MicroElecSurface::IsApplicable(const G4Parti 99 G4MicroElecSurface::IsApplicable(const G4ParticleDefinition& aParticleType) 96 { 100 { 97 return ( aParticleType.GetPDGEncoding() == 1 101 return ( aParticleType.GetPDGEncoding() == 11 ); 98 } 102 } 99 103 100 //....oooOO0OOooo........oooOO0OOooo........oo 104 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 101 105 102 void G4MicroElecSurface::Initialise() << 103 { << 104 if (isInitialised) { return; } << 105 << 106 G4ProductionCutsTable* theCoupleTable = << 107 G4ProductionCutsTable::GetProductionCutsTa << 108 G4int numOfCouples = (G4int)theCoupleTable-> << 109 G4cout << numOfCouples << G4endl; << 110 << 111 for (G4int i = 0; i < numOfCouples; ++i) << 112 { << 113 const G4Material* material = << 114 theCoupleTable->GetMaterialCutsCouple(i) << 115 << 116 G4cout << this->GetProcessName() << ", Mat << 117 << " / " << numOfCouples << " : " < << 118 if (material->GetName() == "Vacuum") << 119 { << 120 tableWF[material->GetName()] = 0; contin << 121 } << 122 G4String mat = material->GetName(); << 123 G4MicroElecMaterialStructure str = G4Micro << 124 tableWF[mat] = str.GetWorkFunction(); << 125 } << 126 << 127 isInitialised = true; << 128 } << 129 << 130 //....oooOO0OOooo........oooOO0OOooo........oo << 131 << 132 void G4MicroElecSurface::BuildPhysicsTable(con 106 void G4MicroElecSurface::BuildPhysicsTable(const G4ParticleDefinition&) 133 { 107 { 134 if (isInitialised) { return; } 108 if (isInitialised) { return; } 135 109 136 G4ProductionCutsTable* theCoupleTable = 110 G4ProductionCutsTable* theCoupleTable = 137 G4ProductionCutsTable::GetProductionCutsTa 111 G4ProductionCutsTable::GetProductionCutsTable(); 138 G4int numOfCouples = (G4int)theCoupleTable-> << 112 G4int numOfCouples = theCoupleTable->GetTableSize(); 139 G4cout << "G4MicroElecSurface::Initialise: N 113 G4cout << "G4MicroElecSurface::Initialise: Ncouples= " 140 << numOfCouples << G4endl; 114 << numOfCouples << G4endl; 141 115 142 for (G4int i = 0; i < numOfCouples; ++i) << 116 for (G4int i = 0; i < numOfCouples; ++i) { 143 { << 144 const G4Material* material = 117 const G4Material* material = 145 theCoupleTable->GetMaterialCutsCouple(i) 118 theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial(); 146 119 147 G4cout << "G4Surface, Material " << i + 1 << 120 G4cout << "G4Surface, Material " << i + 1 << " / " << numOfCouples << " : " << material->GetName() << G4endl; 148 << numOfCouples << " : " << materia << 121 if (material->GetName() == "Vacuum") { tableWF[material->GetName()] = 0; continue; } 149 if (material->GetName() == "Vacuum") << 150 { << 151 tableWF[material->GetName()] = 0; << 152 continue; << 153 } << 154 G4String mat = material->GetName(); 122 G4String mat = material->GetName(); 155 G4MicroElecMaterialStructure str = G4Micro 123 G4MicroElecMaterialStructure str = G4MicroElecMaterialStructure(mat); 156 tableWF[mat] = str.GetWorkFunction(); 124 tableWF[mat] = str.GetWorkFunction(); 157 } 125 } 158 isInitialised = true; 126 isInitialised = true; 159 } 127 } 160 128 161 //....oooOO0OOooo........oooOO0OOooo........oo 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 162 130 163 G4VParticleChange* G4MicroElecSurface::PostSte 131 G4VParticleChange* G4MicroElecSurface::PostStepDoIt(const G4Track& aTrack, const G4Step& aStep) 164 { 132 { 165 << 166 if (!isInitialised) { Initialise(); } << 167 theStatus = UndefinedSurf; 133 theStatus = UndefinedSurf; 168 << 134 169 // Definition of the parameters for the part << 135 //Definition of the parameters for the particle 170 aParticleChange.Initialize(aTrack); 136 aParticleChange.Initialize(aTrack); 171 aParticleChange.ProposeVelocity(aTrack.GetVe 137 aParticleChange.ProposeVelocity(aTrack.GetVelocity()); 172 << 138 173 G4StepPoint* pPreStepPoint = aStep.GetPreSt 139 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint(); 174 G4StepPoint* pPostStepPoint = aStep.GetPostS 140 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint(); 175 << 141 176 material1 = pPreStepPoint -> GetMaterial(); 142 material1 = pPreStepPoint -> GetMaterial(); 177 material2 = pPostStepPoint -> GetMaterial(); 143 material2 = pPostStepPoint -> GetMaterial(); 178 << 144 179 theStatus = UndefinedSurf; << 180 << 181 const G4DynamicParticle* aParticle = aTrack. 145 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 182 << 146 183 theParticleMomentum = aParticle->GetTotalMom 147 theParticleMomentum = aParticle->GetTotalMomentum(); 184 previousMomentum = oldMomentum; 148 previousMomentum = oldMomentum; 185 oldMomentum = aParticle->GetMomentumDirectio 149 oldMomentum = aParticle->GetMomentumDirection(); 186 << 150 187 << 151 //First case: not a boundary 188 // Fisrt case: not a boundary << 152 if (pPostStepPoint->GetStepStatus() != fGeomBoundary || 189 if (pPostStepPoint->GetStepStatus() != fGeom << 153 pPostStepPoint->GetPhysicalVolume() == pPreStepPoint->GetPhysicalVolume()) 190 || pPostStepPoint->GetPhysicalVolume()-> << 154 { 191 { << 155 theStatus = NotAtBoundarySurf; 192 theStatus = NotAtBoundarySurf; << 156 flag_franchissement_surface = false; 193 flag_franchissement_surface = false; << 157 flag_reflexion = false; 194 flag_reflexion = false; << 158 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 195 return G4VDiscreteProcess::PostStepDoIt(aT << 159 196 } << 160 } 197 theStatus = UndefinedSurf; 161 theStatus = UndefinedSurf; 198 << 162 199 // Third case: same material << 163 //Third case: same material 200 if (material1 == material2) 164 if (material1 == material2) 201 { << 165 { 202 theStatus = SameMaterialSurf; << 166 theStatus = SameMaterialSurf; 203 return G4VDiscreteProcess::PostStepDoIt(aT << 167 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 204 } << 168 205 if (verboseLevel > 3) << 169 } 206 { << 170 if (verboseLevel > 0) 207 G4cout << G4endl << " Electron at Boundary << 171 { 208 G4VPhysicalVolume* thePrePV = pPreStepPoin << 172 G4cout << G4endl << " Electron at Boundary! " << G4endl; 209 G4VPhysicalVolume* thePostPV = pPostStepPo << 173 G4VPhysicalVolume* thePrePV = pPreStepPoint->GetPhysicalVolume(); 210 if (thePrePV) G4cout << " thePrePV: " << << 174 G4VPhysicalVolume* thePostPV = pPostStepPoint->GetPhysicalVolume(); 211 if (thePostPV) G4cout << " thePostPV: " << << 175 if (thePrePV) G4cout << " thePrePV: " << thePrePV->GetName() << G4endl; 212 G4cout << " Old Momentum Direction: " << o << 176 if (thePostPV) G4cout << " thePostPV: " << thePostPV->GetName() << G4endl; 213 } << 177 G4cout << " Old Momentum Direction: " << oldMomentum << G4endl; 214 << 178 } 215 // Definition of the parameters for the surf << 179 >> 180 //Definition of the parameters for the surface 216 G4ThreeVector theGlobalPoint = pPostStepPoin 181 G4ThreeVector theGlobalPoint = pPostStepPoint->GetPosition(); 217 << 182 218 G4Navigator* theNavigator = G4Transportation << 183 G4Navigator* theNavigator = 219 GetTransportationManager()->GetNavigatorFo << 184 G4TransportationManager::GetTransportationManager()-> 220 << 185 GetNavigatorForTracking(); >> 186 221 G4bool valid; 187 G4bool valid; >> 188 222 theGlobalNormal = theNavigator->GetGlobalExi 189 theGlobalNormal = theNavigator->GetGlobalExitNormal(theGlobalPoint, &valid); 223 << 190 // G4cout << "Global exit normal = " << theGlobalNormal << " valid = " << valid << G4endl; >> 191 224 if (valid) 192 if (valid) 225 { << 226 theGlobalNormal = -theGlobalNormal; << 227 } << 228 else << 229 { << 230 G4ExceptionDescription ed; << 231 ed << " G4MicroElecSurface/PostStepDoIt(): << 232 << " The Navigator reports that it retu << 233 << G4endl; << 234 G4Exception("G4MuElecSurf::PostStepDoIt", << 235 EventMustBeAborted, ed, << 236 "Invalid Surface Normal - Geom << 237 } << 238 << 239 // Exception: the particle is not in the rig << 240 if (oldMomentum * theGlobalNormal > 0.0) << 241 { << 242 theGlobalNormal = -theGlobalNormal; << 243 } << 244 << 245 if (aTrack.GetStepLength()<=kCarTolerance/2 << 246 { << 247 if (flag_reflexion == true) << 248 { << 249 flag_reflexion = false; << 250 return G4VDiscreteProcess::PostStepDoIt( << 251 } << 252 << 253 theStatus = StepTooSmallSurf; << 254 << 255 G4double energyThreshold_surface = 0.0*eV; << 256 << 257 WorkFunctionTable::iterator postStepWF; << 258 postStepWF = tableWF.find(pPostStepPoint-> << 259 WorkFunctionTable::iterator preStepWF; << 260 preStepWF = tableWF.find(pPreStepPoint->Ge << 261 << 262 if (postStepWF == tableWF.end()) << 263 { 193 { 264 G4String str = "Material "; << 194 theGlobalNormal = -theGlobalNormal; 265 str += pPostStepPoint->GetMaterial()->Ge << 266 G4Exception("G4Surface::G4Surface", "em0 << 267 return nullptr; << 268 } 195 } 269 else if (preStepWF == tableWF.end()) << 196 else 270 { 197 { 271 G4String str = "Material "; << 198 G4ExceptionDescription ed; 272 str += pPreStepPoint->GetMaterial()->Get << 199 ed << " G4MicroElecSurface/PostStepDoIt(): " 273 G4Exception("G4Surface::G4Surface", "em0 << 200 << " The Navigator reports that it returned an invalid normal.\n" 274 return nullptr; << 201 << "PV: " << pPreStepPoint->GetPhysicalVolume()->GetName() >> 202 << " TrackID= " << aTrack.GetTrackID() >> 203 << " Ekin(MeV)= " << aTrack.GetKineticEnergy() >> 204 << " position: " << theGlobalPoint >> 205 << " direction: " << oldMomentum >> 206 << G4endl; >> 207 G4Exception("G4MuElecSurf::PostStepDoIt", "OpBoun01", >> 208 FatalException, ed, >> 209 "Invalid Surface Normal - Geometry must return valid surface normal"); >> 210 return 0; 275 } 211 } 276 else << 212 >> 213 //Exception: the particle is not in the right direction >> 214 if (oldMomentum * theGlobalNormal > 0.0) 277 { 215 { 278 G4double thresholdNew_surface = postStep << 216 theGlobalNormal = -theGlobalNormal; 279 G4double thresholdOld_surface = preStepW << 280 energyThreshold_surface = thresholdNew_s << 281 } 217 } >> 218 >> 219 //Second case: step too small >> 220 //Corrections bug rotation + réflexion >> 221 if (aTrack.GetStepLength()<=kCarTolerance) >> 222 { >> 223 theStatus = StepTooSmallSurf; >> 224 >> 225 WorkFunctionTable::iterator postStepWF; >> 226 postStepWF = tableWF.find(pPostStepPoint->GetMaterial()->GetName()); >> 227 WorkFunctionTable::iterator preStepWF; >> 228 preStepWF = tableWF.find(pPreStepPoint->GetMaterial()->GetName()); >> 229 >> 230 if (postStepWF == tableWF.end()) { >> 231 G4String str = "Material "; >> 232 str += pPostStepPoint->GetMaterial()->GetName() + " not found!"; >> 233 G4Exception("G4Surface::G4Surface", "em0002", FatalException, str); >> 234 return 0; >> 235 } >> 236 >> 237 else if (preStepWF == tableWF.end()) { >> 238 G4String str = "Material "; >> 239 str += pPreStepPoint->GetMaterial()->GetName() + " not found!"; >> 240 G4Exception("G4Surface::G4Surface", "em0002", FatalException, str); >> 241 return 0; >> 242 } >> 243 >> 244 if (pPreStepPoint->GetMaterial() != pPostStepPoint->GetMaterial()) { >> 245 >> 246 flag_franchissement_surface = false; 282 247 283 if (flag_franchissement_surface == true) << 248 if (flag_reflexion == true && flag_normal == true) { 284 { << 249 aParticleChange.ProposeMomentumDirection(-Reflexion(aStep.GetPostStepPoint())); 285 aParticleChange.ProposeEnergy(aStep.GetP << 250 flag_reflexion = false; 286 flag_franchissement_surface = false; << 251 flag_normal = false; 287 } << 252 } 288 if (flag_reflexion == true && flag_normal << 253 } 289 { << 254 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 290 aParticleChange.ProposeMomentumDirection << 291 flag_reflexion = false; << 292 flag_normal = false; << 293 } 255 } 294 return G4VDiscreteProcess::PostStepDoIt(aT << 295 } << 296 << 297 flag_normal = (theGlobalNormal == G4ThreeVec << 298 || theGlobalNormal == G4ThreeVec << 299 256 >> 257 flag_normal = (theGlobalNormal.x() == 0.0 && theGlobalNormal.y() == 0.0); >> 258 300 G4LogicalSurface* Surface = nullptr; 259 G4LogicalSurface* Surface = nullptr; 301 << 260 302 Surface = G4LogicalBorderSurface::GetSurface 261 Surface = G4LogicalBorderSurface::GetSurface 303 (pPreStepPoint ->GetPhysical << 262 (pPreStepPoint ->GetPhysicalVolume(), 304 pPostStepPoint->GetPhysical << 263 pPostStepPoint->GetPhysicalVolume()); 305 << 264 306 if (Surface == nullptr) 265 if (Surface == nullptr) 307 { << 308 G4bool enteredDaughter=(pPostStepPoint->Ge << 309 == pPreStepPoint->Get << 310 if(enteredDaughter) << 311 { 266 { 312 Surface = G4LogicalSkinSurface::GetSurfa << 267 G4bool enteredDaughter=(pPostStepPoint->GetPhysicalVolume() 313 (pPostStepPoint->GetPhysicalVo << 268 ->GetMotherLogical() == 314 << 269 pPreStepPoint->GetPhysicalVolume() 315 if(Surface == nullptr) << 270 ->GetLogicalVolume()); 316 Surface = G4LogicalSkinSurface::GetSurfa << 271 if(enteredDaughter) 317 (pPreStepPoint->GetPhysicalVol << 272 { 318 } << 273 Surface = G4LogicalSkinSurface::GetSurface 319 else << 274 (pPostStepPoint->GetPhysicalVolume()-> 320 { << 275 GetLogicalVolume()); 321 Surface = G4LogicalSkinSurface::GetSurfa << 276 322 (pPreStepPoint->GetPhysicalVol << 277 if(Surface == nullptr) 323 << 278 Surface = G4LogicalSkinSurface::GetSurface 324 if(Surface == nullptr) << 279 (pPreStepPoint->GetPhysicalVolume()-> 325 Surface = G4LogicalSkinSurface::GetSurfa << 280 GetLogicalVolume()); 326 (pPostStepPoint->GetPhysicalVo << 281 } >> 282 else >> 283 { >> 284 Surface = G4LogicalSkinSurface::GetSurface >> 285 (pPreStepPoint->GetPhysicalVolume()-> >> 286 GetLogicalVolume()); >> 287 >> 288 if(Surface == nullptr) >> 289 Surface = G4LogicalSkinSurface::GetSurface >> 290 (pPostStepPoint->GetPhysicalVolume()-> >> 291 GetLogicalVolume()); >> 292 } 327 } 293 } 328 } << 294 329 << 295 330 // Definition of the parameters for the surf << 331 G4VPhysicalVolume* thePrePV = pPreStepPoint- 296 G4VPhysicalVolume* thePrePV = pPreStepPoint->GetPhysicalVolume(); 332 G4VPhysicalVolume* thePostPV = pPostStepPoin 297 G4VPhysicalVolume* thePostPV = pPostStepPoint->GetPhysicalVolume(); 333 << 298 334 energyThreshold = 0.0*eV; << 299 if (thePostPV) 335 G4double energyDelta = 0; << 336 << 337 if ((thePrePV)&&(thePostPV)) << 338 { << 339 WorkFunctionTable::iterator postStepWF; << 340 postStepWF = tableWF.find(thePostPV->GetLo << 341 WorkFunctionTable::iterator preStepWF; << 342 preStepWF = tableWF.find(thePrePV->GetLogi << 343 << 344 if (postStepWF == tableWF.end()) << 345 { << 346 G4String str = "Material "; << 347 str += thePostPV->GetLogicalVolume()->Ge << 348 G4Exception("G4Surface::G4Surface", "em0 << 349 return nullptr; << 350 } << 351 else if (preStepWF == tableWF.end()) << 352 { << 353 G4String str = "Material "; << 354 str += thePrePV->GetLogicalVolume()->Get << 355 G4Exception("G4Surface::G4Surface", "em0 << 356 return nullptr; << 357 } << 358 else << 359 { 300 { 360 G4double thresholdNew = postStepWF->seco << 301 WorkFunctionTable::iterator postStepWF; 361 G4double thresholdOld = preStepWF->secon << 302 postStepWF = tableWF.find(thePostPV->GetLogicalVolume()->GetMaterial()->GetName()); 362 << 303 WorkFunctionTable::iterator preStepWF; 363 energyThreshold = thresholdNew - thresho << 304 preStepWF = tableWF.find(thePrePV->GetLogicalVolume()->GetMaterial()->GetName()); 364 energyDelta = thresholdOld- thresholdNew << 305 >> 306 if (postStepWF == tableWF.end()) { >> 307 G4String str = "Material "; >> 308 str += thePostPV->GetLogicalVolume()->GetMaterial()->GetName() + " not found!"; >> 309 G4Exception("G4Surface::G4Surface", "em0002", FatalException, str); >> 310 return 0; >> 311 } >> 312 >> 313 else if (preStepWF == tableWF.end()) { >> 314 G4String str = "Material "; >> 315 str += thePrePV->GetLogicalVolume()->GetMaterial()->GetName() + " not found!"; >> 316 G4Exception("G4Surface::G4Surface", "em0002", FatalException, str); >> 317 return 0; >> 318 } >> 319 >> 320 else >> 321 { >> 322 G4double thresholdNew = postStepWF->second; >> 323 G4double thresholdOld = preStepWF->second; >> 324 energyThreshold = thresholdNew - thresholdOld; >> 325 } 365 } 326 } 366 } << 327 367 << 328 ekint = pPreStepPoint->GetKineticEnergy(); 368 ekint=aStep.GetPostStepPoint()->GetKineticEn << 369 thetat= GetIncidentAngle(); //angle d'incide 329 thetat= GetIncidentAngle(); //angle d'incidence 370 G4double ekinNormalt=ekint*std::cos(thetat)* 330 G4double ekinNormalt=ekint*std::cos(thetat)*std::cos(thetat); 371 << 331 372 thetaft=std::asin(std::sqrt(ekint/(ekint+ene << 332 G4double atet = std::sqrt(ekint/(ekint+energyThreshold))*std::sin(thetat); 373 if(std::sqrt(ekint/(ekint+energyThreshold))* << 333 thetaft = (atet > 1.0) ? pi*0.5 : std::asin(atet);//Angle de réfraction 374 { << 334 375 thetaft=std::asin(1.0); << 376 } << 377 << 378 G4double aleat=G4UniformRand(); 335 G4double aleat=G4UniformRand(); 379 << 336 380 G4double waveVectort=std::sqrt(2*9.1093826E- << 337 const G4double waveVectort=std::sqrt(2*9.1093826E-31*1.602176487E-19)/(6.6260755E-34/(2.0*pi)); 381 << 338 382 // Parameter for an exponential barrier of p << 339 //Parameter for an exponential barrier of potential (Thèse P68) 383 G4double at=0.5E-10; << 340 const G4double at=0.5E-10; 384 << 341 >> 342 //G4double modif already declared in .hh 385 crossingProbability=0; 343 crossingProbability=0; 386 << 344 387 G4double kft=waveVectort*std::sqrt(ekint+ene 345 G4double kft=waveVectort*std::sqrt(ekint+energyThreshold)*std::cos(thetaft); 388 G4double kit=waveVectort*std::sqrt(ekinNorma 346 G4double kit=waveVectort*std::sqrt(ekinNormalt); 389 << 347 390 crossingProbability=1-(std::pow(std::sinh(pi << 348 G4double yy = std::sinh(pi*at*(kit-kft))/std::sinh(pi*at*(kit+kft)); 391 << 349 crossingProbability = 1 - yy*yy; 392 // First case: the electron crosses the surf << 350 393 if((aleat<=crossingProbability)&&(ekint> ene << 351 //First case: the electron crosses the surface 394 { << 352 if((aleat<=crossingProbability)&&(ekint>std::abs(energyThreshold))) 395 if (aStep.GetPreStepPoint()->GetMaterial() << 396 != aStep.GetPostStepPoint()->GetMaterial( << 397 { 353 { 398 aParticleChange.ProposeEnergy(ekint - en << 354 if (pPreStepPoint->GetMaterial() != pPostStepPoint->GetMaterial()) { 399 flag_franchissement_surface = true; << 355 flag_franchissement_surface = true; 400 } << 356 } 401 << 357 402 thetaft=std::abs(thetaft-thetat); << 358 thetaft=std::abs(thetaft-thetat); 403 << 359 404 G4ThreeVector zVerst = aStep.GetPostStepPo << 360 G4ThreeVector zVerst = aStep.GetPostStepPoint()->GetMomentumDirection(); 405 G4ThreeVector xVerst = zVerst.orthogonal() << 361 G4ThreeVector xVerst = zVerst.orthogonal(); 406 G4ThreeVector yVerst = zVerst.cross(xVerst << 362 G4ThreeVector yVerst = zVerst.cross(xVerst); 407 << 363 408 G4double xDirt = std::sqrt(1. - std::cos(t << 364 G4double cost = std::cos(thetaft); 409 G4double yDirt = xDirt; << 365 G4double xDirt = std::sqrt(1. - cost*cost); 410 << 366 G4double yDirt = xDirt; 411 G4ThreeVector zPrimeVerst=((xDirt*xVerst + << 367 >> 368 G4ThreeVector zPrimeVerst = xDirt*xVerst + yDirt*yVerst + cost*zVerst; 412 369 413 aParticleChange.ProposeMomentumDirection(z << 370 aParticleChange.ProposeMomentumDirection(zPrimeVerst.unit()); 414 } << 415 else if ((aleat > crossingProbability) && (e << 416 { << 417 flag_reflexion = true; << 418 if (flag_normal) << 419 { << 420 aParticleChange.ProposeMomentumDirection << 421 } 371 } 422 else << 372 else if ((aleat > crossingProbability) && (ekint>std::abs(energyThreshold))) 423 { 373 { 424 aParticleChange.ProposeMomentumDirection << 374 flag_reflexion = true; 425 } << 375 if (flag_normal) { aParticleChange.ProposeMomentumDirection(-oldMomentum.unit()); } 426 } << 376 else { aParticleChange.ProposeMomentumDirection(Reflexion(aStep.GetPostStepPoint())); } 427 else << 377 428 { << 429 if (flag_normal) << 430 { << 431 aParticleChange.ProposeMomentumDirection << 432 } << 433 else << 434 { << 435 aParticleChange.ProposeMomentumDirection << 436 } 378 } >> 379 else { >> 380 >> 381 if (flag_normal) { aParticleChange.ProposeMomentumDirection(-oldMomentum.unit()); } >> 382 else { aParticleChange.ProposeMomentumDirection(Reflexion(aStep.GetPostStepPoint())); } 437 flag_reflexion = true; 383 flag_reflexion = true; 438 } 384 } >> 385 439 return G4VDiscreteProcess::PostStepDoIt(aTra 386 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 440 } 387 } 441 388 442 //....oooOO0OOooo........oooOO0OOooo........oo 389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 443 390 444 G4double G4MicroElecSurface::GetMeanFreePath(c 391 G4double G4MicroElecSurface::GetMeanFreePath(const G4Track&, G4double, 445 G << 392 G4ForceCondition* condition) 446 { 393 { 447 *condition = Forced; 394 *condition = Forced; 448 return DBL_MAX; 395 return DBL_MAX; 449 } 396 } 450 397 451 //....oooOO0OOooo........oooOO0OOooo........oo 398 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 452 399 453 G4double G4MicroElecSurface::GetIncidentAngle( 400 G4double G4MicroElecSurface::GetIncidentAngle() 454 { 401 { 455 theFacetNormal=theGlobalNormal; 402 theFacetNormal=theGlobalNormal; 456 403 457 G4double PdotN = oldMomentum * theFacetNorma 404 G4double PdotN = oldMomentum * theFacetNormal; 458 G4double magP= oldMomentum.mag(); 405 G4double magP= oldMomentum.mag(); 459 G4double magN= theFacetNormal.mag(); 406 G4double magN= theFacetNormal.mag(); 460 407 461 G4double incidentangle = pi - std::acos(Pdot 408 G4double incidentangle = pi - std::acos(PdotN/(magP*magN)); 462 409 463 return incidentangle; 410 return incidentangle; 464 } 411 } 465 412 466 //....oooOO0OOooo........oooOO0OOooo........oo 413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 467 414 468 G4ThreeVector G4MicroElecSurface::Reflexion(co 415 G4ThreeVector G4MicroElecSurface::Reflexion(const G4StepPoint* PostStepPoint) 469 { 416 { 470 // Normale << 417 //Normale 471 G4double Nx = theGlobalNormal.x(); 418 G4double Nx = theGlobalNormal.x(); 472 G4double Ny = theGlobalNormal.y(); 419 G4double Ny = theGlobalNormal.y(); 473 G4double Nz = theGlobalNormal.z(); 420 G4double Nz = theGlobalNormal.z(); 474 421 475 // PostStepPoint << 422 //PostStepPoint 476 G4double PSx = PostStepPoint->GetPosition(). 423 G4double PSx = PostStepPoint->GetPosition().x(); 477 G4double PSy = PostStepPoint->GetPosition(). 424 G4double PSy = PostStepPoint->GetPosition().y(); 478 G4double PSz = PostStepPoint->GetPosition(). 425 G4double PSz = PostStepPoint->GetPosition().z(); 479 426 480 // P(alpha,beta,gamma) - PostStep avec trans << 427 //P(alpha,beta,gamma) - PostStep avec translation momentum 481 G4double alpha = PSx + oldMomentum.x(); 428 G4double alpha = PSx + oldMomentum.x(); 482 G4double beta = PSy + oldMomentum.y(); 429 G4double beta = PSy + oldMomentum.y(); 483 G4double gamma = PSz + oldMomentum.z(); 430 G4double gamma = PSz + oldMomentum.z(); 484 G4double r = theGlobalNormal.mag(); 431 G4double r = theGlobalNormal.mag(); 485 G4double x, y, z, d, A, B, PM2x, PM2y, PM2z; 432 G4double x, y, z, d, A, B, PM2x, PM2y, PM2z; 486 d = -(Nx*PSx + Ny*PSy + Nz*PSz); 433 d = -(Nx*PSx + Ny*PSy + Nz*PSz); 487 434 488 if (Ny == 0 && Nx == 0) << 435 if (Ny == 0 && Nx == 0) { 489 { << 490 gamma = -gamma; 436 gamma = -gamma; 491 } << 437 } 492 else << 438 493 { << 439 else { 494 if (Ny == 0) << 440 if (Ny == 0) { 495 { << 441 496 A = (Nz*Nz*alpha) + (Nx*Nx*PSx) + (Nx*Nz 442 A = (Nz*Nz*alpha) + (Nx*Nx*PSx) + (Nx*Nz*(PSz - gamma)); 497 B = r*r; 443 B = r*r; 498 444 499 // M(x,y,z) - Projection de P sur la sur << 445 //M(x,y,z) - Projection de P sur la surface 500 x = A / B; 446 x = A / B; 501 y = beta; 447 y = beta; 502 z = (x - alpha)*(Nz / Nx) + gamma; 448 z = (x - alpha)*(Nz / Nx) + gamma; 503 } 449 } 504 else << 450 505 { << 451 else { >> 452 506 A = (r*r) / Ny; 453 A = (r*r) / Ny; 507 B = (beta / Ny)*(Nx*Nx + Nz*Nz) - (Nx*al 454 B = (beta / Ny)*(Nx*Nx + Nz*Nz) - (Nx*alpha + Nz*gamma + d); 508 455 509 //M(x,y,z) - Projection de P sur la surf 456 //M(x,y,z) - Projection de P sur la surface 510 y = B / A; 457 y = B / A; 511 x = (y - beta)*(Nx / Ny) + alpha; 458 x = (y - beta)*(Nx / Ny) + alpha; 512 z = (y - beta)*(Nz / Ny) + gamma; 459 z = (y - beta)*(Nz / Ny) + gamma; 513 } 460 } 514 461 515 // Vecteur 2*PM << 462 //Vecteur 2*PM 516 PM2x = 2 * (x - alpha); PM2y = 2 * (y - b 463 PM2x = 2 * (x - alpha); PM2y = 2 * (y - beta); PM2z = 2 * (z - gamma); 517 464 518 // Nouveau point P << 465 //Nouveau point P 519 alpha += PM2x; beta += PM2y; gamma += PM2z 466 alpha += PM2x; beta += PM2y; gamma += PM2z; 520 467 521 } 468 } 522 469 523 G4ThreeVector newMomentum = G4ThreeVector(al << 470 G4ThreeVector newMomentum = G4ThreeVector(alpha-PSx,beta-PSy,gamma-PSz); >> 471 524 return newMomentum.unit(); 472 return newMomentum.unit(); 525 } 473 } 526 474 527 //....oooOO0OOooo........oooOO0OOooo........oo 475 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 528 476 529 G4MicroElecSurfaceStatus G4MicroElecSurface::G 477 G4MicroElecSurfaceStatus G4MicroElecSurface::GetStatus() const 530 { 478 { 531 return theStatus; 479 return theStatus; 532 } 480 } 533 481 534 //....oooOO0OOooo........oooOO0OOooo........oo 482 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 535 483 536 484 537 void G4MicroElecSurface::SetFlagFranchissement 485 void G4MicroElecSurface::SetFlagFranchissement() 538 { 486 { 539 flag_franchissement_surface = false; 487 flag_franchissement_surface = false; 540 } 488 } 541 489 542 //....oooOO0OOooo........oooOO0OOooo........oo 490 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 543 491 544 492