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 // $Id: G4LowEnergyPolarizedCompton.cc,v 1.25 2008/05/02 19:23:38 pia Exp $ >> 28 // GEANT4 tag $Name: $ 27 // 29 // 28 // ------------------------------------------- 30 // ------------------------------------------------------------ 29 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file 30 // CERN Geneva Switzerland 32 // CERN Geneva Switzerland 31 // 33 // 32 34 33 // --------- G4LowEnergyPolarizedCompton class 35 // --------- G4LowEnergyPolarizedCompton class ----- 34 // 36 // 35 // by G.Depaola & F.Longo (21 may 20 37 // by G.Depaola & F.Longo (21 may 2001) 36 // 38 // 37 // 21 May 2001 - MGP Modified to inherit 39 // 21 May 2001 - MGP Modified to inherit from G4VDiscreteProcess 38 // Applies same algorit 40 // Applies same algorithm as LowEnergyCompton 39 // if the incoming phot 41 // if the incoming photon is not polarised 40 // Temporary protection 42 // Temporary protection to avoid crash in the case 41 // of polarisation || i 43 // of polarisation || incident photon direction 42 // 44 // 43 // 17 October 2001 - F.Longo - Revised accordi 45 // 17 October 2001 - F.Longo - Revised according to a design iteration 44 // 46 // 45 // 21 February 2002 - F.Longo Revisions with A 47 // 21 February 2002 - F.Longo Revisions with A.Zoglauer and G.Depaola 46 // - better descrip 48 // - better description of parallelism 47 // - system of ref 49 // - system of ref change method improved 48 // 22 January 2003 - V.Ivanchenko Cut per reg 50 // 22 January 2003 - V.Ivanchenko Cut per region 49 // 24 April 2003 - V.Ivanchenko Cut per reg 51 // 24 April 2003 - V.Ivanchenko Cut per region mfpt 50 // 52 // 51 // 53 // 52 // ******************************************* 54 // ************************************************************ 53 // 55 // 54 // Corrections by Rui Curado da Silva (2000) 56 // Corrections by Rui Curado da Silva (2000) 55 // New Implementation by G.Depaola & F.Longo 57 // New Implementation by G.Depaola & F.Longo 56 // 58 // 57 // - sampling of phi 59 // - sampling of phi 58 // - polarization of scattered photon 60 // - polarization of scattered photon 59 // 61 // 60 // ------------------------------------------- 62 // -------------------------------------------------------------- 61 63 62 #include "G4LowEnergyPolarizedCompton.hh" 64 #include "G4LowEnergyPolarizedCompton.hh" 63 #include "Randomize.hh" 65 #include "Randomize.hh" 64 #include "G4PhysicalConstants.hh" << 65 #include "G4SystemOfUnits.hh" << 66 #include "G4ParticleDefinition.hh" 66 #include "G4ParticleDefinition.hh" 67 #include "G4Track.hh" 67 #include "G4Track.hh" 68 #include "G4Step.hh" 68 #include "G4Step.hh" 69 #include "G4ForceCondition.hh" 69 #include "G4ForceCondition.hh" 70 #include "G4Gamma.hh" 70 #include "G4Gamma.hh" 71 #include "G4Electron.hh" 71 #include "G4Electron.hh" 72 #include "G4DynamicParticle.hh" 72 #include "G4DynamicParticle.hh" 73 #include "G4VParticleChange.hh" 73 #include "G4VParticleChange.hh" 74 #include "G4ThreeVector.hh" 74 #include "G4ThreeVector.hh" 75 #include "G4RDVCrossSectionHandler.hh" 75 #include "G4RDVCrossSectionHandler.hh" 76 #include "G4RDCrossSectionHandler.hh" 76 #include "G4RDCrossSectionHandler.hh" 77 #include "G4RDVEMDataSet.hh" 77 #include "G4RDVEMDataSet.hh" 78 #include "G4RDCompositeEMDataSet.hh" 78 #include "G4RDCompositeEMDataSet.hh" 79 #include "G4RDVDataSetAlgorithm.hh" 79 #include "G4RDVDataSetAlgorithm.hh" 80 #include "G4RDLogLogInterpolation.hh" 80 #include "G4RDLogLogInterpolation.hh" 81 #include "G4RDVRangeTest.hh" 81 #include "G4RDVRangeTest.hh" 82 #include "G4RDRangeTest.hh" 82 #include "G4RDRangeTest.hh" 83 #include "G4MaterialCutsCouple.hh" 83 #include "G4MaterialCutsCouple.hh" 84 84 85 // constructor 85 // constructor 86 86 87 G4LowEnergyPolarizedCompton::G4LowEnergyPolari 87 G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton(const G4String& processName) 88 : G4VDiscreteProcess(processName), 88 : G4VDiscreteProcess(processName), 89 lowEnergyLimit (250*eV), // i 89 lowEnergyLimit (250*eV), // initialization 90 highEnergyLimit(100*GeV), 90 highEnergyLimit(100*GeV), 91 intrinsicLowEnergyLimit(10*eV), 91 intrinsicLowEnergyLimit(10*eV), 92 intrinsicHighEnergyLimit(100*GeV) 92 intrinsicHighEnergyLimit(100*GeV) 93 { 93 { 94 if (lowEnergyLimit < intrinsicLowEnergyLimit 94 if (lowEnergyLimit < intrinsicLowEnergyLimit || 95 highEnergyLimit > intrinsicHighEnergyLim 95 highEnergyLimit > intrinsicHighEnergyLimit) 96 { 96 { 97 G4Exception("G4LowEnergyPolarizedCompton 97 G4Exception("G4LowEnergyPolarizedCompton::G4LowEnergyPolarizedCompton()", 98 "OutOfRange", FatalException 98 "OutOfRange", FatalException, 99 "Energy outside intrinsic pr 99 "Energy outside intrinsic process validity range!"); 100 } 100 } 101 101 102 crossSectionHandler = new G4RDCrossSectionHa 102 crossSectionHandler = new G4RDCrossSectionHandler; 103 103 104 104 105 G4RDVDataSetAlgorithm* scatterInterpolation 105 G4RDVDataSetAlgorithm* scatterInterpolation = new G4RDLogLogInterpolation; 106 G4String scatterFile = "comp/ce-sf-"; 106 G4String scatterFile = "comp/ce-sf-"; 107 scatterFunctionData = new G4RDCompositeEMDat 107 scatterFunctionData = new G4RDCompositeEMDataSet(scatterInterpolation,1.,1.); 108 scatterFunctionData->LoadData(scatterFile); 108 scatterFunctionData->LoadData(scatterFile); 109 109 110 meanFreePathTable = 0; 110 meanFreePathTable = 0; 111 111 112 rangeTest = new G4RDRangeTest; 112 rangeTest = new G4RDRangeTest; 113 113 114 // For Doppler broadening 114 // For Doppler broadening 115 shellData.SetOccupancyData(); 115 shellData.SetOccupancyData(); 116 116 117 117 118 if (verboseLevel > 0) 118 if (verboseLevel > 0) 119 { 119 { 120 G4cout << GetProcessName() << " is crea 120 G4cout << GetProcessName() << " is created " << G4endl 121 << "Energy range: " 121 << "Energy range: " 122 << lowEnergyLimit / keV << " keV - " 122 << lowEnergyLimit / keV << " keV - " 123 << highEnergyLimit / GeV << " GeV" 123 << highEnergyLimit / GeV << " GeV" 124 << G4endl; 124 << G4endl; 125 } 125 } 126 } 126 } 127 127 128 // destructor 128 // destructor 129 129 130 G4LowEnergyPolarizedCompton::~G4LowEnergyPolar 130 G4LowEnergyPolarizedCompton::~G4LowEnergyPolarizedCompton() 131 { 131 { 132 delete meanFreePathTable; 132 delete meanFreePathTable; 133 delete crossSectionHandler; 133 delete crossSectionHandler; 134 delete scatterFunctionData; 134 delete scatterFunctionData; 135 delete rangeTest; 135 delete rangeTest; 136 } 136 } 137 137 138 138 139 void G4LowEnergyPolarizedCompton::BuildPhysics 139 void G4LowEnergyPolarizedCompton::BuildPhysicsTable(const G4ParticleDefinition& ) 140 { 140 { 141 141 142 crossSectionHandler->Clear(); 142 crossSectionHandler->Clear(); 143 G4String crossSectionFile = "comp/ce-cs-"; 143 G4String crossSectionFile = "comp/ce-cs-"; 144 crossSectionHandler->LoadData(crossSectionFi 144 crossSectionHandler->LoadData(crossSectionFile); 145 delete meanFreePathTable; 145 delete meanFreePathTable; 146 meanFreePathTable = crossSectionHandler->Bui 146 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials(); 147 147 148 // For Doppler broadening 148 // For Doppler broadening 149 G4String file = "/doppler/shell-doppler"; 149 G4String file = "/doppler/shell-doppler"; 150 shellData.LoadData(file); 150 shellData.LoadData(file); 151 151 152 } 152 } 153 153 154 G4VParticleChange* G4LowEnergyPolarizedCompton 154 G4VParticleChange* G4LowEnergyPolarizedCompton::PostStepDoIt(const G4Track& aTrack, 155 const G4Step& aStep) 155 const G4Step& aStep) 156 { 156 { 157 // The scattered gamma energy is sampled acc 157 // The scattered gamma energy is sampled according to Klein - Nishina formula. 158 // The random number techniques of Butcher & 158 // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15). 159 // GEANT4 internal units 159 // GEANT4 internal units 160 // 160 // 161 // Note : Effects due to binding of atomic e 161 // Note : Effects due to binding of atomic electrons are negliged. 162 162 163 aParticleChange.Initialize(aTrack); 163 aParticleChange.Initialize(aTrack); 164 164 165 // Dynamic particle quantities 165 // Dynamic particle quantities 166 const G4DynamicParticle* incidentPhoton = aT 166 const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle(); 167 G4double gammaEnergy0 = incidentPhoton->GetK 167 G4double gammaEnergy0 = incidentPhoton->GetKineticEnergy(); 168 G4ThreeVector gammaPolarization0 = incidentP 168 G4ThreeVector gammaPolarization0 = incidentPhoton->GetPolarization(); 169 169 170 // gammaPolarization0 = gammaPolarization0. 170 // gammaPolarization0 = gammaPolarization0.unit(); // 171 171 172 // Protection: a polarisation parallel to th 172 // Protection: a polarisation parallel to the 173 // direction causes problems; 173 // direction causes problems; 174 // in that case find a random polarization 174 // in that case find a random polarization 175 175 176 G4ThreeVector gammaDirection0 = incidentPhot 176 G4ThreeVector gammaDirection0 = incidentPhoton->GetMomentumDirection(); 177 // ---- MGP ---- Next two lines commented ou 177 // ---- MGP ---- Next two lines commented out to remove compilation warnings 178 // G4double scalarproduct = gammaPolarizatio 178 // G4double scalarproduct = gammaPolarization0.dot(gammaDirection0); 179 // G4double angle = gammaPolarization0.angle 179 // G4double angle = gammaPolarization0.angle(gammaDirection0); 180 180 181 // Make sure that the polarization vector is 181 // Make sure that the polarization vector is perpendicular to the 182 // gamma direction. If not 182 // gamma direction. If not 183 183 184 if(!(gammaPolarization0.isOrthogonal(gammaDi 184 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0)) 185 { // only for testing now 185 { // only for testing now 186 gammaPolarization0 = GetRandomPolarizati 186 gammaPolarization0 = GetRandomPolarization(gammaDirection0); 187 } 187 } 188 else 188 else 189 { 189 { 190 if ( gammaPolarization0.howOrthogonal(ga 190 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0) 191 { 191 { 192 gammaPolarization0 = GetPerpendicularPolar 192 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0); 193 } 193 } 194 } 194 } 195 195 196 // End of Protection 196 // End of Protection 197 197 198 // Within energy limit? 198 // Within energy limit? 199 199 200 if(gammaEnergy0 <= lowEnergyLimit) 200 if(gammaEnergy0 <= lowEnergyLimit) 201 { 201 { 202 aParticleChange.ProposeTrackStatus(fStop 202 aParticleChange.ProposeTrackStatus(fStopAndKill); 203 aParticleChange.ProposeEnergy(0.); 203 aParticleChange.ProposeEnergy(0.); 204 aParticleChange.ProposeLocalEnergyDeposi 204 aParticleChange.ProposeLocalEnergyDeposit(gammaEnergy0); 205 return G4VDiscreteProcess::PostStepDoIt( 205 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep); 206 } 206 } 207 207 208 G4double E0_m = gammaEnergy0 / electron_mass 208 G4double E0_m = gammaEnergy0 / electron_mass_c2 ; 209 209 210 // Select randomly one element in the curren 210 // Select randomly one element in the current material 211 211 212 const G4MaterialCutsCouple* couple = aTrack. 212 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple(); 213 G4int Z = crossSectionHandler->SelectRandomA 213 G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0); 214 214 215 // Sample the energy and the polarization of 215 // Sample the energy and the polarization of the scattered photon 216 216 217 G4double epsilon, epsilonSq, onecost, sinThe 217 G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ; 218 218 219 G4double epsilon0 = 1./(1. + 2*E0_m); 219 G4double epsilon0 = 1./(1. + 2*E0_m); 220 G4double epsilon0Sq = epsilon0*epsilon0; 220 G4double epsilon0Sq = epsilon0*epsilon0; 221 G4double alpha1 = - std::log(epsilon0); 221 G4double alpha1 = - std::log(epsilon0); 222 G4double alpha2 = 0.5*(1.- epsilon0Sq); 222 G4double alpha2 = 0.5*(1.- epsilon0Sq); 223 223 224 G4double wlGamma = h_Planck*c_light/gammaEne 224 G4double wlGamma = h_Planck*c_light/gammaEnergy0; 225 G4double gammaEnergy1; 225 G4double gammaEnergy1; 226 G4ThreeVector gammaDirection1; 226 G4ThreeVector gammaDirection1; 227 227 228 do { 228 do { 229 if ( alpha1/(alpha1+alpha2) > G4UniformRan 229 if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) 230 { 230 { 231 epsilon = std::exp(-alpha1*G4UniformRand() 231 epsilon = std::exp(-alpha1*G4UniformRand()); 232 epsilonSq = epsilon*epsilon; 232 epsilonSq = epsilon*epsilon; 233 } 233 } 234 else 234 else 235 { 235 { 236 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4 236 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand(); 237 epsilon = std::sqrt(epsilonSq); 237 epsilon = std::sqrt(epsilonSq); 238 } 238 } 239 239 240 onecost = (1.- epsilon)/(epsilon*E0_m); 240 onecost = (1.- epsilon)/(epsilon*E0_m); 241 sinThetaSqr = onecost*(2.-onecost); 241 sinThetaSqr = onecost*(2.-onecost); 242 242 243 // Protection 243 // Protection 244 if (sinThetaSqr > 1.) 244 if (sinThetaSqr > 1.) 245 { 245 { 246 if (verboseLevel>0) G4cout 246 if (verboseLevel>0) G4cout 247 << " -- Warning -- G4LowEnergyPolarizedCom 247 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt " 248 << "sin(theta)**2 = " 248 << "sin(theta)**2 = " 249 << sinThetaSqr 249 << sinThetaSqr 250 << "; set to 1" 250 << "; set to 1" 251 << G4endl; 251 << G4endl; 252 sinThetaSqr = 1.; 252 sinThetaSqr = 1.; 253 } 253 } 254 if (sinThetaSqr < 0.) 254 if (sinThetaSqr < 0.) 255 { 255 { 256 if (verboseLevel>0) G4cout 256 if (verboseLevel>0) G4cout 257 << " -- Warning -- G4LowEnergyPolarizedCom 257 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt " 258 << "sin(theta)**2 = " 258 << "sin(theta)**2 = " 259 << sinThetaSqr 259 << sinThetaSqr 260 << "; set to 0" 260 << "; set to 0" 261 << G4endl; 261 << G4endl; 262 sinThetaSqr = 0.; 262 sinThetaSqr = 0.; 263 } 263 } 264 // End protection 264 // End protection 265 265 266 G4double x = std::sqrt(onecost/2.) / (wlG 266 G4double x = std::sqrt(onecost/2.) / (wlGamma/cm);; 267 G4double scatteringFunction = scatterFunct 267 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1); 268 greject = (1. - epsilon*sinThetaSqr/(1.+ e 268 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction; 269 269 270 } while(greject < G4UniformRand()*Z); 270 } while(greject < G4UniformRand()*Z); 271 271 272 272 273 // ***************************************** 273 // **************************************************** 274 // Phi determination 274 // Phi determination 275 // ***************************************** 275 // **************************************************** 276 276 277 G4double phi = SetPhi(epsilon,sinThetaSqr); 277 G4double phi = SetPhi(epsilon,sinThetaSqr); 278 278 279 // 279 // 280 // scattered gamma angles. ( Z - axis along 280 // scattered gamma angles. ( Z - axis along the parent gamma) 281 // 281 // 282 282 283 G4double cosTheta = 1. - onecost; 283 G4double cosTheta = 1. - onecost; 284 284 285 // Protection 285 // Protection 286 286 287 if (cosTheta > 1.) 287 if (cosTheta > 1.) 288 { 288 { 289 if (verboseLevel>0) G4cout 289 if (verboseLevel>0) G4cout 290 << " -- Warning -- G4LowEnergyPolarizedCompt 290 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt " 291 << "cosTheta = " 291 << "cosTheta = " 292 << cosTheta 292 << cosTheta 293 << "; set to 1" 293 << "; set to 1" 294 << G4endl; 294 << G4endl; 295 cosTheta = 1.; 295 cosTheta = 1.; 296 } 296 } 297 if (cosTheta < -1.) 297 if (cosTheta < -1.) 298 { 298 { 299 if (verboseLevel>0) G4cout 299 if (verboseLevel>0) G4cout 300 << " -- Warning -- G4LowEnergyPolarizedCompt 300 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt " 301 << "cosTheta = " 301 << "cosTheta = " 302 << cosTheta 302 << cosTheta 303 << "; set to -1" 303 << "; set to -1" 304 << G4endl; 304 << G4endl; 305 cosTheta = -1.; 305 cosTheta = -1.; 306 } 306 } 307 // End protection 307 // End protection 308 308 309 309 310 G4double sinTheta = std::sqrt (sinThetaSqr); 310 G4double sinTheta = std::sqrt (sinThetaSqr); 311 311 312 // Protection 312 // Protection 313 if (sinTheta > 1.) 313 if (sinTheta > 1.) 314 { 314 { 315 if (verboseLevel>0) G4cout 315 if (verboseLevel>0) G4cout 316 << " -- Warning -- G4LowEnergyPolarizedCompt 316 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt " 317 << "sinTheta = " 317 << "sinTheta = " 318 << sinTheta 318 << sinTheta 319 << "; set to 1" 319 << "; set to 1" 320 << G4endl; 320 << G4endl; 321 sinTheta = 1.; 321 sinTheta = 1.; 322 } 322 } 323 if (sinTheta < -1.) 323 if (sinTheta < -1.) 324 { 324 { 325 if (verboseLevel>0) G4cout 325 if (verboseLevel>0) G4cout 326 << " -- Warning -- G4LowEnergyPolarizedCompt 326 << " -- Warning -- G4LowEnergyPolarizedCompton::PostStepDoIt " 327 << "sinTheta = " 327 << "sinTheta = " 328 << sinTheta 328 << sinTheta 329 << "; set to -1" 329 << "; set to -1" 330 << G4endl; 330 << G4endl; 331 sinTheta = -1.; 331 sinTheta = -1.; 332 } 332 } 333 // End protection 333 // End protection 334 334 335 335 336 G4double dirx = sinTheta*std::cos(phi); 336 G4double dirx = sinTheta*std::cos(phi); 337 G4double diry = sinTheta*std::sin(phi); 337 G4double diry = sinTheta*std::sin(phi); 338 G4double dirz = cosTheta ; 338 G4double dirz = cosTheta ; 339 339 340 340 341 // oneCosT , eom 341 // oneCosT , eom 342 342 343 343 344 344 345 // Doppler broadening - Method based on: 345 // Doppler broadening - Method based on: 346 // Y. Namito, S. Ban and H. Hirayama, 346 // Y. Namito, S. Ban and H. Hirayama, 347 // "Implementation of the Doppler Broadening 347 // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code" 348 // NIM A 349, pp. 489-494, 1994 348 // NIM A 349, pp. 489-494, 1994 349 349 350 // Maximum number of sampling iterations 350 // Maximum number of sampling iterations 351 351 352 G4int maxDopplerIterations = 1000; 352 G4int maxDopplerIterations = 1000; 353 G4double bindingE = 0.; 353 G4double bindingE = 0.; 354 G4double photonEoriginal = epsilon * gammaEn 354 G4double photonEoriginal = epsilon * gammaEnergy0; 355 G4double photonE = -1.; 355 G4double photonE = -1.; 356 G4int iteration = 0; 356 G4int iteration = 0; 357 G4double eMax = gammaEnergy0; 357 G4double eMax = gammaEnergy0; 358 358 359 do 359 do 360 { 360 { 361 iteration++; 361 iteration++; 362 // Select shell based on shell occupancy 362 // Select shell based on shell occupancy 363 G4int shell = shellData.SelectRandomShel 363 G4int shell = shellData.SelectRandomShell(Z); 364 bindingE = shellData.BindingEnergy(Z,she 364 bindingE = shellData.BindingEnergy(Z,shell); 365 365 366 eMax = gammaEnergy0 - bindingE; 366 eMax = gammaEnergy0 - bindingE; 367 367 368 // Randomly sample bound electron moment 368 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units) 369 G4double pSample = profileData.RandomSel 369 G4double pSample = profileData.RandomSelectMomentum(Z,shell); 370 // Rescale from atomic units 370 // Rescale from atomic units 371 G4double pDoppler = pSample * fine_struc 371 G4double pDoppler = pSample * fine_structure_const; 372 G4double pDoppler2 = pDoppler * pDoppler 372 G4double pDoppler2 = pDoppler * pDoppler; 373 G4double var2 = 1. + onecost * E0_m; 373 G4double var2 = 1. + onecost * E0_m; 374 G4double var3 = var2*var2 - pDoppler2; 374 G4double var3 = var2*var2 - pDoppler2; 375 G4double var4 = var2 - pDoppler2 * cosTh 375 G4double var4 = var2 - pDoppler2 * cosTheta; 376 G4double var = var4*var4 - var3 + pDoppl 376 G4double var = var4*var4 - var3 + pDoppler2 * var3; 377 if (var > 0.) 377 if (var > 0.) 378 { 378 { 379 G4double varSqrt = std::sqrt(var); 379 G4double varSqrt = std::sqrt(var); 380 G4double scale = gammaEnergy0 / var3; 380 G4double scale = gammaEnergy0 / var3; 381 // Random select either root 381 // Random select either root 382 if (G4UniformRand() < 0.5) photonE = (var4 382 if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale; 383 else photonE = (var4 + varSqrt) * scale; 383 else photonE = (var4 + varSqrt) * scale; 384 } 384 } 385 else 385 else 386 { 386 { 387 photonE = -1.; 387 photonE = -1.; 388 } 388 } 389 } while ( iteration <= maxDopplerIterations 389 } while ( iteration <= maxDopplerIterations && 390 (photonE < 0. || photonE > eMax || phot 390 (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) ); 391 391 392 // End of recalculation of photon energy wit 392 // End of recalculation of photon energy with Doppler broadening 393 // Revert to original if maximum number of i 393 // Revert to original if maximum number of iterations threshold has been reached 394 if (iteration >= maxDopplerIterations) 394 if (iteration >= maxDopplerIterations) 395 { 395 { 396 photonE = photonEoriginal; 396 photonE = photonEoriginal; 397 bindingE = 0.; 397 bindingE = 0.; 398 } 398 } 399 399 400 gammaEnergy1 = photonE; 400 gammaEnergy1 = photonE; 401 401 402 // G4cout << "--> PHOTONENERGY1 = " << photo 402 // G4cout << "--> PHOTONENERGY1 = " << photonE/keV << G4endl; 403 403 404 404 405 /// Doppler Broadeing 405 /// Doppler Broadeing 406 406 407 407 408 408 409 409 410 // 410 // 411 // update G4VParticleChange for the scattere 411 // update G4VParticleChange for the scattered photon 412 // 412 // 413 413 414 // gammaEnergy1 = epsilon*gammaEnergy0; 414 // gammaEnergy1 = epsilon*gammaEnergy0; 415 415 416 416 417 // New polarization 417 // New polarization 418 418 419 G4ThreeVector gammaPolarization1 = SetNewPol 419 G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon, 420 sinThetaSqr, 420 sinThetaSqr, 421 phi, 421 phi, 422 cosTheta); 422 cosTheta); 423 423 424 // Set new direction 424 // Set new direction 425 G4ThreeVector tmpDirection1( dirx,diry,dirz 425 G4ThreeVector tmpDirection1( dirx,diry,dirz ); 426 gammaDirection1 = tmpDirection1; 426 gammaDirection1 = tmpDirection1; 427 427 428 // Change reference frame. 428 // Change reference frame. 429 429 430 SystemOfRefChange(gammaDirection0,gammaDirec 430 SystemOfRefChange(gammaDirection0,gammaDirection1, 431 gammaPolarization0,gammaPolarization1) 431 gammaPolarization0,gammaPolarization1); 432 432 433 if (gammaEnergy1 > 0.) 433 if (gammaEnergy1 > 0.) 434 { 434 { 435 aParticleChange.ProposeEnergy( gammaEner 435 aParticleChange.ProposeEnergy( gammaEnergy1 ) ; 436 aParticleChange.ProposeMomentumDirection 436 aParticleChange.ProposeMomentumDirection( gammaDirection1 ); 437 aParticleChange.ProposePolarization( gam 437 aParticleChange.ProposePolarization( gammaPolarization1 ); 438 } 438 } 439 else 439 else 440 { 440 { 441 aParticleChange.ProposeEnergy(0.) ; 441 aParticleChange.ProposeEnergy(0.) ; 442 aParticleChange.ProposeTrackStatus(fStop 442 aParticleChange.ProposeTrackStatus(fStopAndKill); 443 } 443 } 444 444 445 // 445 // 446 // kinematic of the scattered electron 446 // kinematic of the scattered electron 447 // 447 // 448 448 449 G4double ElecKineEnergy = gammaEnergy0 - gam 449 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE; 450 450 451 451 452 // Generate the electron only if with large 452 // Generate the electron only if with large enough range w.r.t. cuts and safety 453 453 454 G4double safety = aStep.GetPostStepPoint()-> 454 G4double safety = aStep.GetPostStepPoint()->GetSafety(); 455 455 456 456 457 if (rangeTest->Escape(G4Electron::Electron() 457 if (rangeTest->Escape(G4Electron::Electron(),couple,ElecKineEnergy,safety)) 458 { 458 { 459 G4double ElecMomentum = std::sqrt(ElecKi 459 G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2)); 460 G4ThreeVector ElecDirection((gammaEnergy 460 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 - 461 gammaEnergy1 * gammaDirection1) * ( 461 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum)); 462 G4DynamicParticle* electron = new G4Dyna 462 G4DynamicParticle* electron = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ; 463 aParticleChange.SetNumberOfSecondaries(1 463 aParticleChange.SetNumberOfSecondaries(1); 464 aParticleChange.AddSecondary(electron); 464 aParticleChange.AddSecondary(electron); 465 // aParticleChange.ProposeLocalEner 465 // aParticleChange.ProposeLocalEnergyDeposit(0.); 466 aParticleChange.ProposeLocalEnergyDeposi 466 aParticleChange.ProposeLocalEnergyDeposit(bindingE); 467 } 467 } 468 else 468 else 469 { 469 { 470 aParticleChange.SetNumberOfSecondaries(0 470 aParticleChange.SetNumberOfSecondaries(0); 471 aParticleChange.ProposeLocalEnergyDeposi 471 aParticleChange.ProposeLocalEnergyDeposit(ElecKineEnergy+bindingE); 472 } 472 } 473 473 474 return G4VDiscreteProcess::PostStepDoIt( aTr 474 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep); 475 475 476 } 476 } 477 477 478 478 479 G4double G4LowEnergyPolarizedCompton::SetPhi(G 479 G4double G4LowEnergyPolarizedCompton::SetPhi(G4double energyRate, 480 G4double sinSqrTh) 480 G4double sinSqrTh) 481 { 481 { 482 G4double rand1; 482 G4double rand1; 483 G4double rand2; 483 G4double rand2; 484 G4double phiProbability; 484 G4double phiProbability; 485 G4double phi; 485 G4double phi; 486 G4double a, b; 486 G4double a, b; 487 487 488 do 488 do 489 { 489 { 490 rand1 = G4UniformRand(); 490 rand1 = G4UniformRand(); 491 rand2 = G4UniformRand(); 491 rand2 = G4UniformRand(); 492 phiProbability=0.; 492 phiProbability=0.; 493 phi = twopi*rand1; 493 phi = twopi*rand1; 494 494 495 a = 2*sinSqrTh; 495 a = 2*sinSqrTh; 496 b = energyRate + 1/energyRate; 496 b = energyRate + 1/energyRate; 497 497 498 phiProbability = 1 - (a/b)*(std::cos(phi 498 phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi)); 499 499 500 500 501 501 502 } 502 } 503 while ( rand2 > phiProbability ); 503 while ( rand2 > phiProbability ); 504 return phi; 504 return phi; 505 } 505 } 506 506 507 507 508 G4ThreeVector G4LowEnergyPolarizedCompton::Set 508 G4ThreeVector G4LowEnergyPolarizedCompton::SetPerpendicularVector(G4ThreeVector& a) 509 { 509 { 510 G4double dx = a.x(); 510 G4double dx = a.x(); 511 G4double dy = a.y(); 511 G4double dy = a.y(); 512 G4double dz = a.z(); 512 G4double dz = a.z(); 513 G4double x = dx < 0.0 ? -dx : dx; 513 G4double x = dx < 0.0 ? -dx : dx; 514 G4double y = dy < 0.0 ? -dy : dy; 514 G4double y = dy < 0.0 ? -dy : dy; 515 G4double z = dz < 0.0 ? -dz : dz; 515 G4double z = dz < 0.0 ? -dz : dz; 516 if (x < y) { 516 if (x < y) { 517 return x < z ? G4ThreeVector(-dy,dx,0) : G 517 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy); 518 }else{ 518 }else{ 519 return y < z ? G4ThreeVector(dz,0,-dx) : G 519 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0); 520 } 520 } 521 } 521 } 522 522 523 G4ThreeVector G4LowEnergyPolarizedCompton::Get 523 G4ThreeVector G4LowEnergyPolarizedCompton::GetRandomPolarization(G4ThreeVector& direction0) 524 { 524 { 525 G4ThreeVector d0 = direction0.unit(); 525 G4ThreeVector d0 = direction0.unit(); 526 G4ThreeVector a1 = SetPerpendicularVector(d0 526 G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal 527 G4ThreeVector a0 = a1.unit(); // unit vector 527 G4ThreeVector a0 = a1.unit(); // unit vector 528 528 529 G4double rand1 = G4UniformRand(); 529 G4double rand1 = G4UniformRand(); 530 530 531 G4double angle = twopi*rand1; // random pola 531 G4double angle = twopi*rand1; // random polar angle 532 G4ThreeVector b0 = d0.cross(a0); // cross pr 532 G4ThreeVector b0 = d0.cross(a0); // cross product 533 533 534 G4ThreeVector c; 534 G4ThreeVector c; 535 535 536 c.setX(std::cos(angle)*(a0.x())+std::sin(ang 536 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x()); 537 c.setY(std::cos(angle)*(a0.y())+std::sin(ang 537 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y()); 538 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang 538 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z()); 539 539 540 G4ThreeVector c0 = c.unit(); 540 G4ThreeVector c0 = c.unit(); 541 541 542 return c0; 542 return c0; 543 543 544 } 544 } 545 545 546 546 547 G4ThreeVector G4LowEnergyPolarizedCompton::Get 547 G4ThreeVector G4LowEnergyPolarizedCompton::GetPerpendicularPolarization 548 (const G4ThreeVector& gammaDirection, const G4 548 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const 549 { 549 { 550 550 551 // 551 // 552 // The polarization of a photon is always pe 552 // The polarization of a photon is always perpendicular to its momentum direction. 553 // Therefore this function removes those vec 553 // Therefore this function removes those vector component of gammaPolarization, which 554 // points in direction of gammaDirection 554 // points in direction of gammaDirection 555 // 555 // 556 // Mathematically we search the projection o 556 // Mathematically we search the projection of the vector a on the plane E, where n is the 557 // plains normal vector. 557 // plains normal vector. 558 // The basic equation can be found in each g 558 // The basic equation can be found in each geometry book (e.g. Bronstein): 559 // p = a - (a o n)/(n o n)*n 559 // p = a - (a o n)/(n o n)*n 560 560 561 return gammaPolarization - gammaPolarization 561 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection; 562 } 562 } 563 563 564 564 565 G4ThreeVector G4LowEnergyPolarizedCompton::Set 565 G4ThreeVector G4LowEnergyPolarizedCompton::SetNewPolarization(G4double epsilon, 566 G4double sinSqrTh, 566 G4double sinSqrTh, 567 G4double phi, 567 G4double phi, 568 G4double costheta) 568 G4double costheta) 569 { 569 { 570 G4double rand1; 570 G4double rand1; 571 G4double rand2; 571 G4double rand2; 572 G4double cosPhi = std::cos(phi); 572 G4double cosPhi = std::cos(phi); 573 G4double sinPhi = std::sin(phi); 573 G4double sinPhi = std::sin(phi); 574 G4double sinTheta = std::sqrt(sinSqrTh); 574 G4double sinTheta = std::sqrt(sinSqrTh); 575 G4double cosSqrPhi = cosPhi*cosPhi; 575 G4double cosSqrPhi = cosPhi*cosPhi; 576 // G4double cossqrth = 1.-sinSqrTh; 576 // G4double cossqrth = 1.-sinSqrTh; 577 // G4double sinsqrphi = sinPhi*sinPhi; 577 // G4double sinsqrphi = sinPhi*sinPhi; 578 G4double normalisation = std::sqrt(1. - cosS 578 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh); 579 579 580 580 581 // Determination of Theta 581 // Determination of Theta 582 582 583 // ---- MGP ---- Commented out the following 583 // ---- MGP ---- Commented out the following 3 lines to avoid compilation 584 // warnings (unused variables) 584 // warnings (unused variables) 585 // G4double thetaProbability; 585 // G4double thetaProbability; 586 G4double theta; 586 G4double theta; 587 // G4double a, b; 587 // G4double a, b; 588 // G4double cosTheta; 588 // G4double cosTheta; 589 589 590 /* 590 /* 591 591 592 depaola method 592 depaola method 593 593 594 do 594 do 595 { 595 { 596 rand1 = G4UniformRand(); 596 rand1 = G4UniformRand(); 597 rand2 = G4UniformRand(); 597 rand2 = G4UniformRand(); 598 thetaProbability=0.; 598 thetaProbability=0.; 599 theta = twopi*rand1; 599 theta = twopi*rand1; 600 a = 4*normalisation*normalisation; 600 a = 4*normalisation*normalisation; 601 b = (epsilon + 1/epsilon) - 2; 601 b = (epsilon + 1/epsilon) - 2; 602 thetaProbability = (b + a*std::cos(theta 602 thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b); 603 cosTheta = std::cos(theta); 603 cosTheta = std::cos(theta); 604 } 604 } 605 while ( rand2 > thetaProbability ); 605 while ( rand2 > thetaProbability ); 606 606 607 G4double cosBeta = cosTheta; 607 G4double cosBeta = cosTheta; 608 608 609 */ 609 */ 610 610 611 611 612 // Dan Xu method (IEEE TNS, 52, 1160 (2005)) 612 // Dan Xu method (IEEE TNS, 52, 1160 (2005)) 613 613 614 rand1 = G4UniformRand(); 614 rand1 = G4UniformRand(); 615 rand2 = G4UniformRand(); 615 rand2 = G4UniformRand(); 616 616 617 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsi 617 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi)) 618 { 618 { 619 if (rand2<0.5) 619 if (rand2<0.5) 620 theta = pi/2.0; 620 theta = pi/2.0; 621 else 621 else 622 theta = 3.0*pi/2.0; 622 theta = 3.0*pi/2.0; 623 } 623 } 624 else 624 else 625 { 625 { 626 if (rand2<0.5) 626 if (rand2<0.5) 627 theta = 0; 627 theta = 0; 628 else 628 else 629 theta = pi; 629 theta = pi; 630 } 630 } 631 G4double cosBeta = std::cos(theta); 631 G4double cosBeta = std::cos(theta); 632 G4double sinBeta = std::sqrt(1-cosBeta*cosBe 632 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta); 633 633 634 G4ThreeVector gammaPolarization1; 634 G4ThreeVector gammaPolarization1; 635 635 636 G4double xParallel = normalisation*cosBeta; 636 G4double xParallel = normalisation*cosBeta; 637 G4double yParallel = -(sinSqrTh*cosPhi*sinPh 637 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation; 638 G4double zParallel = -(costheta*sinTheta*cos 638 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation; 639 G4double xPerpendicular = 0.; 639 G4double xPerpendicular = 0.; 640 G4double yPerpendicular = (costheta)*sinBeta 640 G4double yPerpendicular = (costheta)*sinBeta/normalisation; 641 G4double zPerpendicular = -(sinTheta*sinPhi) 641 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation; 642 642 643 G4double xTotal = (xParallel + xPerpendicula 643 G4double xTotal = (xParallel + xPerpendicular); 644 G4double yTotal = (yParallel + yPerpendicula 644 G4double yTotal = (yParallel + yPerpendicular); 645 G4double zTotal = (zParallel + zPerpendicula 645 G4double zTotal = (zParallel + zPerpendicular); 646 646 647 gammaPolarization1.setX(xTotal); 647 gammaPolarization1.setX(xTotal); 648 gammaPolarization1.setY(yTotal); 648 gammaPolarization1.setY(yTotal); 649 gammaPolarization1.setZ(zTotal); 649 gammaPolarization1.setZ(zTotal); 650 650 651 return gammaPolarization1; 651 return gammaPolarization1; 652 652 653 } 653 } 654 654 655 655 656 void G4LowEnergyPolarizedCompton::SystemOfRefC 656 void G4LowEnergyPolarizedCompton::SystemOfRefChange(G4ThreeVector& direction0, 657 G4ThreeVector& direction1, 657 G4ThreeVector& direction1, 658 G4ThreeVector& polarization0, 658 G4ThreeVector& polarization0, 659 G4ThreeVector& polarization1) 659 G4ThreeVector& polarization1) 660 { 660 { 661 // direction0 is the original photon directi 661 // direction0 is the original photon direction ---> z 662 // polarization0 is the original photon pola 662 // polarization0 is the original photon polarization ---> x 663 // need to specify y axis in the real refere 663 // need to specify y axis in the real reference frame ---> y 664 G4ThreeVector Axis_Z0 = direction0.unit(); 664 G4ThreeVector Axis_Z0 = direction0.unit(); 665 G4ThreeVector Axis_X0 = polarization0.unit() 665 G4ThreeVector Axis_X0 = polarization0.unit(); 666 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_ 666 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed; 667 667 668 G4double direction_x = direction1.getX(); 668 G4double direction_x = direction1.getX(); 669 G4double direction_y = direction1.getY(); 669 G4double direction_y = direction1.getY(); 670 G4double direction_z = direction1.getZ(); 670 G4double direction_z = direction1.getZ(); 671 671 672 direction1 = (direction_x*Axis_X0 + directio 672 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit(); 673 G4double polarization_x = polarization1.getX 673 G4double polarization_x = polarization1.getX(); 674 G4double polarization_y = polarization1.getY 674 G4double polarization_y = polarization1.getY(); 675 G4double polarization_z = polarization1.getZ 675 G4double polarization_z = polarization1.getZ(); 676 676 677 polarization1 = (polarization_x*Axis_X0 + po 677 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit(); 678 678 679 } 679 } 680 680 681 681 682 G4bool G4LowEnergyPolarizedCompton::IsApplicab 682 G4bool G4LowEnergyPolarizedCompton::IsApplicable(const G4ParticleDefinition& particle) 683 { 683 { 684 return ( &particle == G4Gamma::Gamma() ); 684 return ( &particle == G4Gamma::Gamma() ); 685 } 685 } 686 686 687 687 688 G4double G4LowEnergyPolarizedCompton::GetMeanF 688 G4double G4LowEnergyPolarizedCompton::GetMeanFreePath(const G4Track& track, 689 G4double, 689 G4double, 690 G4ForceCondition*) 690 G4ForceCondition*) 691 { 691 { 692 const G4DynamicParticle* photon = track.GetD 692 const G4DynamicParticle* photon = track.GetDynamicParticle(); 693 G4double energy = photon->GetKineticEnergy() 693 G4double energy = photon->GetKineticEnergy(); 694 const G4MaterialCutsCouple* couple = track.G 694 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple(); 695 size_t materialIndex = couple->GetIndex(); 695 size_t materialIndex = couple->GetIndex(); 696 G4double meanFreePath; 696 G4double meanFreePath; 697 if (energy > highEnergyLimit) meanFreePath = 697 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex); 698 else if (energy < lowEnergyLimit) meanFreePa 698 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX; 699 else meanFreePath = meanFreePathTable->FindV 699 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex); 700 return meanFreePath; 700 return meanFreePath; 701 } 701 } 702 702 703 703 704 704 705 705 706 706 707 707 708 708 709 709 710 710 711 711 712 712 713 713 714 714 715 715 716 716