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 // $Id: G4LivermorePolarizedComptonModel.cc,v 1.1 2008/10/30 14:16:35 sincerti Exp $ >> 27 // GEANT4 tag $Name: geant4-09-02-patch-01 $ 26 // 28 // 27 // Authors: G.Depaola & F.Longo << 28 // << 29 // History: << 30 // ------- << 31 // << 32 // 05 Apr 2021 J Allison added quantum entan << 33 // If the photons have been "tagged" as "quant << 34 // G4eplusAnnihilation for annihilation into 2 << 35 // here if - and only if - both photons suffer << 36 // predictions from Pryce and Ward, Nature No << 37 // Physical Review 73 (1948) p.440. Experiment << 38 // entanglement in the MeV regime and its appl << 39 // D. Watts, J. Allison et al., Nature Communi << 40 // https://doi.org/10.1038/s41467-021-22907-5. << 41 // << 42 // 02 May 2009 S Incerti as V. Ivanchenko pr << 43 // << 44 // Cleanup initialisation and generation of se << 45 // - apply internal high-ener << 46 // - do not apply low-energy << 47 // - remove GetMeanFreePath m << 48 // - added protection against << 49 // - use G4ElementSelector << 50 29 51 #include "G4LivermorePolarizedComptonModel.hh" 30 #include "G4LivermorePolarizedComptonModel.hh" 52 #include "G4PhysicalConstants.hh" << 53 #include "G4SystemOfUnits.hh" << 54 #include "G4AutoLock.hh" << 55 #include "G4Electron.hh" << 56 #include "G4ParticleChangeForGamma.hh" << 57 #include "G4LossTableManager.hh" << 58 #include "G4VAtomDeexcitation.hh" << 59 #include "G4AtomicShell.hh" << 60 #include "G4Gamma.hh" << 61 #include "G4ShellData.hh" << 62 #include "G4DopplerProfile.hh" << 63 #include "G4Log.hh" << 64 #include "G4Exp.hh" << 65 #include "G4Pow.hh" << 66 #include "G4LogLogInterpolation.hh" << 67 #include "G4PhysicsModelCatalog.hh" << 68 #include "G4EntanglementAuxInfo.hh" << 69 #include "G4eplusAnnihilationEntanglementClipB << 70 31 71 //....oooOO0OOooo........oooOO0OOooo........oo 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 72 33 73 using namespace std; 34 using namespace std; 74 namespace { G4Mutex LivermorePolarizedComptonM << 75 << 76 << 77 G4PhysicsFreeVector* G4LivermorePolarizedCompt << 78 G4ShellData* G4LivermorePolarizedCompton << 79 G4DopplerProfile* G4LivermorePolarizedCompton << 80 G4CompositeEMDataSet* G4LivermorePolarizedComp << 81 35 82 //....oooOO0OOooo........oooOO0OOooo........oo 36 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 83 37 84 G4LivermorePolarizedComptonModel::G4LivermoreP << 38 G4LivermorePolarizedComptonModel::G4LivermorePolarizedComptonModel(const G4ParticleDefinition*, 85 :G4VEmModel(nam),isInitialised(false) << 39 const G4String& nam) 86 { << 40 :G4VEmModel(nam),isInitialised(false) 87 verboseLevel= 1; << 41 { >> 42 lowEnergyLimit = 250 * eV; // SI - Could be 10 eV ? >> 43 highEnergyLimit = 100 * GeV; >> 44 SetLowEnergyLimit(lowEnergyLimit); >> 45 SetHighEnergyLimit(highEnergyLimit); >> 46 >> 47 verboseLevel= 0; 88 // Verbosity scale: 48 // Verbosity scale: 89 // 0 = nothing 49 // 0 = nothing 90 // 1 = warning for energy non-conservation 50 // 1 = warning for energy non-conservation 91 // 2 = details of energy budget 51 // 2 = details of energy budget 92 // 3 = calculation of cross sections, file o 52 // 3 = calculation of cross sections, file openings, sampling of atoms 93 // 4 = entering in methods 53 // 4 = entering in methods 94 54 95 if( verboseLevel>1 ) << 55 G4cout << "Livermore Polarized Compton is constructed " << G4endl 96 G4cout << "Livermore Polarized Compton is << 56 << "Energy range: " 97 << 57 << lowEnergyLimit / keV << " keV - " 98 //Mark this model as "applicable" for atomic << 58 << highEnergyLimit / GeV << " GeV" 99 SetDeexcitationFlag(true); << 59 << G4endl; 100 << 60 101 fParticleChange = nullptr; << 102 fAtomDeexcitation = nullptr; << 103 fEntanglementModelID = G4PhysicsModelCatalog << 104 } 61 } 105 62 106 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 107 64 108 G4LivermorePolarizedComptonModel::~G4Livermore 65 G4LivermorePolarizedComptonModel::~G4LivermorePolarizedComptonModel() 109 { 66 { 110 if(IsMaster()) { << 67 delete meanFreePathTable; 111 delete shellData; << 68 delete crossSectionHandler; 112 shellData = nullptr; << 69 delete scatterFunctionData; 113 delete profileData; << 114 profileData = nullptr; << 115 delete scatterFunctionData; << 116 scatterFunctionData = nullptr; << 117 for(G4int i=0; i<maxZ; ++i) { << 118 if(data[i]) { << 119 delete data[i]; << 120 data[i] = nullptr; << 121 } << 122 } << 123 } << 124 } 70 } 125 71 126 //....oooOO0OOooo........oooOO0OOooo........oo 72 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 127 73 128 void G4LivermorePolarizedComptonModel::Initial 74 void G4LivermorePolarizedComptonModel::Initialise(const G4ParticleDefinition* particle, 129 const G 75 const G4DataVector& cuts) 130 { 76 { 131 if (verboseLevel > 1) << 77 if (verboseLevel > 3) 132 G4cout << "Calling G4LivermorePolarizedCom 78 G4cout << "Calling G4LivermorePolarizedComptonModel::Initialise()" << G4endl; 133 79 134 // Initialise element selector << 80 InitialiseElementSelectors(particle,cuts); 135 if(IsMaster()) { << 136 // Access to elements << 137 const char* path = G4FindDataDir("G4LEDATA << 138 81 139 G4ProductionCutsTable* theCoupleTable = << 82 // Energy limits 140 G4ProductionCutsTable::GetProductionCuts << 83 141 << 84 if (LowEnergyLimit() < lowEnergyLimit) 142 G4int numOfCouples = (G4int)theCoupleTable << 85 { 143 << 86 G4cout << "G4LivermorePolarizedComptonModel: low energy limit increased from " << 144 for(G4int i=0; i<numOfCouples; ++i) { << 87 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit << " eV" << G4endl; 145 const G4Material* material = << 88 SetLowEnergyLimit(lowEnergyLimit); 146 theCoupleTable->GetMaterialCutsCouple( << 147 const G4ElementVector* theElementVector << 148 std::size_t nelm = material->GetNumberOf << 149 << 150 for (std::size_t j=0; j<nelm; ++j) { << 151 G4int Z = G4lrint((*theElementVector)[ << 152 if(Z < 1) { Z = 1; } << 153 else if(Z > maxZ){ Z = maxZ; } << 154 << 155 if( (!data[Z]) ) { ReadData(Z, path); << 156 } << 157 } 89 } 158 << 90 159 // For Doppler broadening << 91 if (HighEnergyLimit() > highEnergyLimit) 160 if(!shellData) { << 92 { 161 shellData = new G4ShellData(); << 93 G4cout << "G4LivermorePolarizedComptonModel: high energy limit decreased from " << 162 shellData->SetOccupancyData(); << 94 HighEnergyLimit()/GeV << " GeV to " << highEnergyLimit << " GeV" << G4endl; 163 G4String file = "/doppler/shell-doppler" << 95 SetHighEnergyLimit(highEnergyLimit); 164 shellData->LoadData(file); << 165 } 96 } 166 if(!profileData) { profileData = new G4Dop << 167 97 168 // Scattering Function << 98 // Reading of data files - all materials are read 169 if(!scatterFunctionData) << 170 { << 171 << 172 G4VDataSetAlgorithm* scatterInterpolation = << 173 G4String scatterFile = "comp/ce-sf-"; << 174 scatterFunctionData = new G4CompositeEMDataS << 175 scatterFunctionData->LoadData(scatterFile); << 176 } << 177 << 178 InitialiseElementSelectors(particle, cuts) << 179 } << 180 << 181 if (verboseLevel > 2) { << 182 G4cout << "Loaded cross section files" << << 183 } << 184 << 185 if( verboseLevel>1 ) { << 186 G4cout << "G4LivermoreComptonModel is init << 187 << "Energy range: " << 188 << LowEnergyLimit() / eV << " eV - " << 189 << HighEnergyLimit() / GeV << " GeV" << 190 << G4endl; << 191 } << 192 // << 193 if(isInitialised) { return; } << 194 99 195 fParticleChange = GetParticleChangeForGamma( << 100 crossSectionHandler = new G4CrossSectionHandler; 196 fAtomDeexcitation = G4LossTableManager::Ins << 101 crossSectionHandler->Clear(); 197 isInitialised = true; << 102 G4String crossSectionFile = "comp/ce-cs-"; 198 } << 103 crossSectionHandler->LoadData(crossSectionFile); 199 104 >> 105 meanFreePathTable = 0; >> 106 meanFreePathTable = crossSectionHandler->BuildMeanFreePathForMaterials(); 200 107 201 void G4LivermorePolarizedComptonModel::Initial << 108 G4VDataSetAlgorithm* scatterInterpolation = new G4LogLogInterpolation; 202 G4VEmModel* masterModel) << 109 G4String scatterFile = "comp/ce-sf-"; 203 { << 110 scatterFunctionData = new G4CompositeEMDataSet(scatterInterpolation, 1., 1.); 204 SetElementSelectors(masterModel->GetElementS << 111 scatterFunctionData->LoadData(scatterFile); 205 } << 206 112 207 //....oooOO0OOooo........oooOO0OOooo........oo << 113 // For Doppler broadening >> 114 shellData.SetOccupancyData(); >> 115 G4String file = "/doppler/shell-doppler"; >> 116 shellData.LoadData(file); 208 117 209 void G4LivermorePolarizedComptonModel::ReadDat << 118 // 210 { << 119 if (verboseLevel > 2) 211 if (verboseLevel > 1) << 120 G4cout << "Loaded cross section files for Livermore Polarized Compton model" << G4endl; 212 { << 121 213 G4cout << "G4LivermorePolarizedComptonMo << 122 G4cout << "Livermore Polarized Compton model is initialized " << G4endl 214 << G4endl; << 123 << "Energy range: " 215 } << 124 << LowEnergyLimit() / keV << " keV - " 216 if(data[Z]) { return; } << 125 << HighEnergyLimit() / GeV << " GeV" 217 const char* datadir = path; << 126 << G4endl; 218 if(!datadir) << 127 219 { << 128 // 220 datadir = G4FindDataDir("G4LEDATA"); << 221 if(!datadir) << 222 { << 223 G4Exception("G4LivermorePolarizedComptonMo << 224 "em0006",FatalException, << 225 "Environment variable G4LEDATA not d << 226 return; << 227 } << 228 } << 229 << 230 data[Z] = new G4PhysicsFreeVector(); << 231 129 232 std::ostringstream ost; << 130 if(isInitialised) return; 233 ost << datadir << "/livermore/comp/ce-cs-" < << 131 234 std::ifstream fin(ost.str().c_str()); << 132 if(pParticleChange) 235 << 133 fParticleChange = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange); 236 if( !fin.is_open()) << 134 else 237 { << 135 fParticleChange = new G4ParticleChangeForGamma(); 238 G4ExceptionDescription ed; << 136 239 ed << "G4LivermorePolarizedComptonModel << 137 isInitialised = true; 240 << "> is not opened!" << G4endl; << 241 G4Exception("G4LivermoreComptonModel::Re << 242 "em0003",FatalException, << 243 ed,"G4LEDATA version should be G4EMLOW8. << 244 return; << 245 } else { << 246 if(verboseLevel > 3) { << 247 G4cout << "File " << ost.str() << 248 << " is opened by G4LivermorePolarizedC << 249 } << 250 data[Z]->Retrieve(fin, true); << 251 data[Z]->ScaleVector(MeV, MeV*barn); << 252 } << 253 fin.close(); << 254 } 138 } 255 139 256 //....oooOO0OOooo........oooOO0OOooo........oo 140 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 257 141 258 G4double G4LivermorePolarizedComptonModel::Com 142 G4double G4LivermorePolarizedComptonModel::ComputeCrossSectionPerAtom( 259 const G 143 const G4ParticleDefinition*, 260 G 144 G4double GammaEnergy, 261 G 145 G4double Z, G4double, 262 G 146 G4double, G4double) 263 { 147 { 264 if (verboseLevel > 3) 148 if (verboseLevel > 3) 265 G4cout << "Calling ComputeCrossSectionPerA 149 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermorePolarizedComptonModel" << G4endl; 266 150 267 G4double cs = 0.0; << 151 G4double cs = crossSectionHandler->FindValue(G4int(Z), GammaEnergy); 268 << 269 if (GammaEnergy < LowEnergyLimit()) << 270 return 0.0; << 271 << 272 G4int intZ = G4lrint(Z); << 273 if(intZ < 1 || intZ > maxZ) { return cs; } << 274 << 275 G4PhysicsFreeVector* pv = data[intZ]; << 276 << 277 // if element was not initialised << 278 // do initialisation safely for MT mode << 279 if(!pv) << 280 { << 281 InitialiseForElement(0, intZ); << 282 pv = data[intZ]; << 283 if(!pv) { return cs; } << 284 } << 285 << 286 G4int n = G4int(pv->GetVectorLength() - 1); << 287 G4double e1 = pv->Energy(0); << 288 G4double e2 = pv->Energy(n); << 289 << 290 if(GammaEnergy <= e1) { cs = GammaEnerg << 291 else if(GammaEnergy <= e2) { cs = pv->Value( << 292 else if(GammaEnergy > e2) { cs = pv->Value( << 293 << 294 return cs; 152 return cs; 295 << 296 } 153 } 297 154 298 //....oooOO0OOooo........oooOO0OOooo........oo 155 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 299 156 300 void G4LivermorePolarizedComptonModel::SampleS 157 void G4LivermorePolarizedComptonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* fvect, 301 const G4MaterialCutsCouple* co 158 const G4MaterialCutsCouple* couple, 302 const G4DynamicParticle* aDyna 159 const G4DynamicParticle* aDynamicGamma, 303 G4double, 160 G4double, 304 G4double) 161 G4double) 305 { 162 { 306 // The scattered gamma energy is sampled acc 163 // The scattered gamma energy is sampled according to Klein - Nishina formula. 307 // The random number techniques of Butcher & 164 // The random number techniques of Butcher & Messel are used (Nuc Phys 20(1960),15). 308 // GEANT4 internal units 165 // GEANT4 internal units 309 // 166 // 310 // Note : Effects due to binding of atomic e 167 // Note : Effects due to binding of atomic electrons are negliged. 311 168 312 if (verboseLevel > 3) 169 if (verboseLevel > 3) 313 G4cout << "Calling SampleSecondaries() of 170 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedComptonModel" << G4endl; 314 171 315 G4double gammaEnergy0 = aDynamicGamma->GetKi 172 G4double gammaEnergy0 = aDynamicGamma->GetKineticEnergy(); 316 << 317 // do nothing below the threshold << 318 // should never get here because the XS is z << 319 if (gammaEnergy0 < LowEnergyLimit()) << 320 return ; << 321 << 322 G4ThreeVector gammaPolarization0 = aDynamicG 173 G4ThreeVector gammaPolarization0 = aDynamicGamma->GetPolarization(); 323 174 324 // Protection: a polarisation parallel to th 175 // Protection: a polarisation parallel to the 325 // direction causes problems; 176 // direction causes problems; 326 // in that case find a random polarization 177 // in that case find a random polarization >> 178 327 G4ThreeVector gammaDirection0 = aDynamicGamm 179 G4ThreeVector gammaDirection0 = aDynamicGamma->GetMomentumDirection(); 328 180 329 // Make sure that the polarization vector is 181 // Make sure that the polarization vector is perpendicular to the 330 // gamma direction. If not 182 // gamma direction. If not >> 183 331 if(!(gammaPolarization0.isOrthogonal(gammaDi 184 if(!(gammaPolarization0.isOrthogonal(gammaDirection0, 1e-6))||(gammaPolarization0.mag()==0)) 332 { // only for testing now 185 { // only for testing now 333 gammaPolarization0 = GetRandomPolarizati 186 gammaPolarization0 = GetRandomPolarization(gammaDirection0); 334 } 187 } 335 else 188 else 336 { 189 { 337 if ( gammaPolarization0.howOrthogonal(ga 190 if ( gammaPolarization0.howOrthogonal(gammaDirection0) != 0) 338 { 191 { 339 gammaPolarization0 = GetPerpendicularPolar 192 gammaPolarization0 = GetPerpendicularPolarization(gammaDirection0, gammaPolarization0); 340 } 193 } 341 } 194 } >> 195 342 // End of Protection 196 // End of Protection 343 197 >> 198 // Within energy limit? >> 199 >> 200 if(gammaEnergy0 <= lowEnergyLimit) >> 201 { >> 202 fParticleChange->ProposeTrackStatus(fStopAndKill); >> 203 fParticleChange->SetProposedKineticEnergy(0.); >> 204 fParticleChange->ProposeLocalEnergyDeposit(gammaEnergy0); >> 205 // SI - IS THE FOLLOWING RETURN NECESSARY ? >> 206 return; >> 207 } >> 208 344 G4double E0_m = gammaEnergy0 / electron_mass 209 G4double E0_m = gammaEnergy0 / electron_mass_c2 ; 345 210 346 // Select randomly one element in the curren 211 // Select randomly one element in the current material 347 //G4int Z = crossSectionHandler->SelectRando << 212 348 const G4ParticleDefinition* particle = aDyn << 213 G4int Z = crossSectionHandler->SelectRandomAtom(couple,gammaEnergy0); 349 const G4Element* elm = SelectRandomAtom(coup << 350 G4int Z = (G4int)elm->GetZ(); << 351 214 352 // Sample the energy and the polarization of 215 // Sample the energy and the polarization of the scattered photon >> 216 353 G4double epsilon, epsilonSq, onecost, sinThe 217 G4double epsilon, epsilonSq, onecost, sinThetaSqr, greject ; 354 218 355 G4double epsilon0Local = 1./(1. + 2*E0_m); << 219 G4double epsilon0 = 1./(1. + 2*E0_m); 356 G4double epsilon0Sq = epsilon0Local*epsilon0 << 220 G4double epsilon0Sq = epsilon0*epsilon0; 357 G4double alpha1 = - G4Log(epsilon0Local); << 221 G4double alpha1 = - std::log(epsilon0); 358 G4double alpha2 = 0.5*(1.- epsilon0Sq); 222 G4double alpha2 = 0.5*(1.- epsilon0Sq); 359 223 360 G4double wlGamma = h_Planck*c_light/gammaEne 224 G4double wlGamma = h_Planck*c_light/gammaEnergy0; 361 G4double gammaEnergy1; 225 G4double gammaEnergy1; 362 G4ThreeVector gammaDirection1; 226 G4ThreeVector gammaDirection1; 363 227 364 do { 228 do { 365 if ( alpha1/(alpha1+alpha2) > G4UniformRan 229 if ( alpha1/(alpha1+alpha2) > G4UniformRand() ) 366 { 230 { 367 epsilon = G4Exp(-alpha1*G4UniformRand()); << 231 epsilon = std::exp(-alpha1*G4UniformRand()); 368 epsilonSq = epsilon*epsilon; 232 epsilonSq = epsilon*epsilon; 369 } 233 } 370 else 234 else 371 { 235 { 372 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4 236 epsilonSq = epsilon0Sq + (1.- epsilon0Sq)*G4UniformRand(); 373 epsilon = std::sqrt(epsilonSq); 237 epsilon = std::sqrt(epsilonSq); 374 } 238 } 375 239 376 onecost = (1.- epsilon)/(epsilon*E0_m); 240 onecost = (1.- epsilon)/(epsilon*E0_m); 377 sinThetaSqr = onecost*(2.-onecost); 241 sinThetaSqr = onecost*(2.-onecost); 378 242 379 // Protection 243 // Protection 380 if (sinThetaSqr > 1.) 244 if (sinThetaSqr > 1.) 381 { 245 { 382 G4cout 246 G4cout 383 << " -- Warning -- G4LivermorePolarizedCom 247 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 384 << "sin(theta)**2 = " 248 << "sin(theta)**2 = " 385 << sinThetaSqr 249 << sinThetaSqr 386 << "; set to 1" 250 << "; set to 1" 387 << G4endl; 251 << G4endl; 388 sinThetaSqr = 1.; 252 sinThetaSqr = 1.; 389 } 253 } 390 if (sinThetaSqr < 0.) 254 if (sinThetaSqr < 0.) 391 { 255 { 392 G4cout 256 G4cout 393 << " -- Warning -- G4LivermorePolarizedCom 257 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 394 << "sin(theta)**2 = " 258 << "sin(theta)**2 = " 395 << sinThetaSqr 259 << sinThetaSqr 396 << "; set to 0" 260 << "; set to 0" 397 << G4endl; 261 << G4endl; 398 sinThetaSqr = 0.; 262 sinThetaSqr = 0.; 399 } 263 } 400 // End protection 264 // End protection 401 265 402 G4double x = std::sqrt(onecost/2.) / (wlG 266 G4double x = std::sqrt(onecost/2.) / (wlGamma/cm);; 403 G4double scatteringFunction = scatterFunct 267 G4double scatteringFunction = scatterFunctionData->FindValue(x,Z-1); 404 greject = (1. - epsilon*sinThetaSqr/(1.+ e 268 greject = (1. - epsilon*sinThetaSqr/(1.+ epsilonSq))*scatteringFunction; 405 269 406 } while(greject < G4UniformRand()*Z); 270 } while(greject < G4UniformRand()*Z); 407 271 408 272 409 // ***************************************** 273 // **************************************************** 410 // Phi determination 274 // Phi determination 411 // ***************************************** 275 // **************************************************** >> 276 412 G4double phi = SetPhi(epsilon,sinThetaSqr); 277 G4double phi = SetPhi(epsilon,sinThetaSqr); 413 278 414 // 279 // 415 // scattered gamma angles. ( Z - axis along 280 // scattered gamma angles. ( Z - axis along the parent gamma) 416 // 281 // >> 282 417 G4double cosTheta = 1. - onecost; 283 G4double cosTheta = 1. - onecost; 418 284 419 // Protection 285 // Protection >> 286 420 if (cosTheta > 1.) 287 if (cosTheta > 1.) 421 { 288 { 422 G4cout 289 G4cout 423 << " -- Warning -- G4LivermorePolarizedCompt 290 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 424 << "cosTheta = " 291 << "cosTheta = " 425 << cosTheta 292 << cosTheta 426 << "; set to 1" 293 << "; set to 1" 427 << G4endl; 294 << G4endl; 428 cosTheta = 1.; 295 cosTheta = 1.; 429 } 296 } 430 if (cosTheta < -1.) 297 if (cosTheta < -1.) 431 { 298 { 432 G4cout 299 G4cout 433 << " -- Warning -- G4LivermorePolarizedCompt 300 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 434 << "cosTheta = " 301 << "cosTheta = " 435 << cosTheta 302 << cosTheta 436 << "; set to -1" 303 << "; set to -1" 437 << G4endl; 304 << G4endl; 438 cosTheta = -1.; 305 cosTheta = -1.; 439 } 306 } 440 // End protection 307 // End protection 441 << 308 >> 309 442 G4double sinTheta = std::sqrt (sinThetaSqr); 310 G4double sinTheta = std::sqrt (sinThetaSqr); 443 311 444 // Protection 312 // Protection 445 if (sinTheta > 1.) 313 if (sinTheta > 1.) 446 { 314 { 447 G4cout 315 G4cout 448 << " -- Warning -- G4LivermorePolarizedCompt 316 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 449 << "sinTheta = " 317 << "sinTheta = " 450 << sinTheta 318 << sinTheta 451 << "; set to 1" 319 << "; set to 1" 452 << G4endl; 320 << G4endl; 453 sinTheta = 1.; 321 sinTheta = 1.; 454 } 322 } 455 if (sinTheta < -1.) 323 if (sinTheta < -1.) 456 { 324 { 457 G4cout 325 G4cout 458 << " -- Warning -- G4LivermorePolarizedCompt 326 << " -- Warning -- G4LivermorePolarizedComptonModel::SampleSecondaries " 459 << "sinTheta = " 327 << "sinTheta = " 460 << sinTheta 328 << sinTheta 461 << "; set to -1" 329 << "; set to -1" 462 << G4endl; 330 << G4endl; 463 sinTheta = -1.; 331 sinTheta = -1.; 464 } 332 } 465 // End protection 333 // End protection 466 334 467 // Check for entanglement and re-sample phi << 468 << 469 const auto* auxInfo << 470 = fParticleChange->GetCurrentTrack()->GetAux << 471 if (auxInfo) { << 472 const auto* entanglementAuxInfo = dynamic_ << 473 if (entanglementAuxInfo) { << 474 auto* clipBoard = dynamic_cast<G4eplusAn << 475 (entanglementAuxInfo->GetEntanglementCli << 476 if (clipBoard) { << 477 // This is an entangled photon from ep << 478 // If this is the first scatter of the << 479 // phi on the clipboard. << 480 // If this is the first scatter of the << 481 // phi of the first scatter of the fir << 482 // theta of the second photon, to samp << 483 if (clipBoard->IsTrack1Measurement()) << 484 // Check we have the relevant track. << 485 // necessary but I want to be sure t << 486 // entangled system are properly pai << 487 // Note: the tracking manager pops t << 488 // will rely on that. (If not, the << 489 // more complicated to ensure we mat << 490 // So our track 1 is clipboard track << 491 if (clipBoard->GetTrackB() == fParti << 492 // This is the first scatter of th << 493 // // Debug << 494 // auto* track1 = fPart << 495 // G4cout << 496 // << "This is the firs << 497 // << "\nTrack: " << tr << 498 // << ", Parent: " << t << 499 // << ", Name: " << cli << 500 // << G4endl; << 501 // // End debug << 502 clipBoard->ResetTrack1Measurement( << 503 // Store cos(theta),phi of first p << 504 clipBoard->SetComptonCosTheta1(cos << 505 clipBoard->SetComptonPhi1(phi); << 506 } << 507 } else if (clipBoard->IsTrack2Measurem << 508 // Check we have the relevant track. << 509 // Remember our track 2 is clipboard << 510 if (clipBoard->GetTrackA() == fParti << 511 // This is the first scatter of th << 512 // // Debug << 513 // auto* track2 = fPart << 514 // G4cout << 515 // << "This is the firs << 516 // << "\nTrack: " << tr << 517 // << ", Parent: " << t << 518 // << ", Name: " << cli << 519 // << G4endl; << 520 // // End debug << 521 clipBoard->ResetTrack2Measurement( << 522 << 523 // Get cos(theta),phi of first pho << 524 const G4double& cosTheta1 = clipBo << 525 const G4double& phi1 = clipBoard-> << 526 // For clarity make aliases for th << 527 const G4double& cosTheta2 = cosThe << 528 G4double& phi2 = phi; << 529 // G4cout << "cosTheta1 << 530 // G4cout << "cosTheta2 << 531 << 532 // Re-sample phi << 533 // Draw the difference of azimutha << 534 // A + B * cos(2*deltaPhi), or rat << 535 // C = A / (A + |B|) and D = B / ( << 536 const G4double sin2Theta1 = 1.-cos << 537 const G4double sin2Theta2 = 1.-cos << 538 << 539 // Pryce and Ward, Nature No 4065 << 540 auto* g4Pow = G4Pow::GetInstance() << 541 const G4double A = << 542 ((g4Pow->powN(1.-cosTheta1,3))+2.) << 543 ((g4Pow->powN(2.-cosTheta1,3)*g4Po << 544 const G4double B = -(sin2Theta1*si << 545 ((g4Pow->powN(2.-cosTheta1,2)*g4Po << 546 << 547 // // Snyder et al, Physical Re << 548 // // (This is an alternative f << 549 // const G4double& k0 = gammaEn << 550 // const G4double k1 = k0/(2.-c << 551 // const G4double k2 = k0/(2.-c << 552 // const G4double gamma1 = k1/k << 553 // const G4double gamma2 = k2/k << 554 // const G4double A1 = gamma1*g << 555 // const G4double B1 = 2.*sin2T << 556 // // That's A1 + B1*sin2(delta << 557 // const G4double A = A1 + 0.5* << 558 // const G4double B = -0.5*B1; << 559 << 560 const G4double maxValue = A + std: << 561 const G4double C = A / maxValue; << 562 const G4double D = B / maxValue; << 563 // G4cout << "A,B,C,D: " << A < << 564 << 565 // Sample delta phi << 566 G4double deltaPhi; << 567 const G4int maxCount = 999999; << 568 G4int iCount = 0; << 569 for (; iCount < maxCount; ++iCount << 570 deltaPhi = twopi * G4UniformRand << 571 if (G4UniformRand() < C + D * co << 572 } << 573 if (iCount >= maxCount ) { << 574 G4cout << "G4LivermorePolarizedC << 575 << "Re-sampled delta phi not fou << 576 << " tries - carrying on anyway. << 577 } << 578 << 579 // Thus, the desired second photon << 580 phi2 = deltaPhi - phi1 + halfpi; << 581 // The minus sign is in above stat << 582 // annihilation photons are in opp << 583 // are measured in the opposite di << 584 // halfpi is added for the followi << 585 // In this function phi is relativ << 586 // SystemOfRefChange below. We kno << 587 // the polarisations of the two an << 588 // to each other, i.e., halfpi dif << 589 // Furthermore, only sin(phi) and << 590 // need to place any range constra << 591 // if (phi2 > pi) { << 592 // phi2 -= twopi; << 593 // } << 594 // if (phi2 < -pi) { << 595 // phi2 += twopi; << 596 // } << 597 } << 598 } << 599 } << 600 } << 601 } << 602 << 603 // End of entanglement << 604 335 605 G4double dirx = sinTheta*std::cos(phi); 336 G4double dirx = sinTheta*std::cos(phi); 606 G4double diry = sinTheta*std::sin(phi); 337 G4double diry = sinTheta*std::sin(phi); 607 G4double dirz = cosTheta ; 338 G4double dirz = cosTheta ; 608 << 339 >> 340 609 // oneCosT , eom 341 // oneCosT , eom 610 342 611 // Doppler broadening - Method based on: 343 // Doppler broadening - Method based on: 612 // Y. Namito, S. Ban and H. Hirayama, 344 // Y. Namito, S. Ban and H. Hirayama, 613 // "Implementation of the Doppler Broadening 345 // "Implementation of the Doppler Broadening of a Compton-Scattered Photon Into the EGS4 Code" 614 // NIM A 349, pp. 489-494, 1994 346 // NIM A 349, pp. 489-494, 1994 615 347 616 // Maximum number of sampling iterations 348 // Maximum number of sampling iterations 617 static G4int maxDopplerIterations = 1000; << 349 >> 350 G4int maxDopplerIterations = 1000; 618 G4double bindingE = 0.; 351 G4double bindingE = 0.; 619 G4double photonEoriginal = epsilon * gammaEn 352 G4double photonEoriginal = epsilon * gammaEnergy0; 620 G4double photonE = -1.; 353 G4double photonE = -1.; 621 G4int iteration = 0; 354 G4int iteration = 0; 622 G4double eMax = gammaEnergy0; 355 G4double eMax = gammaEnergy0; 623 356 624 G4int shellIdx = 0; << 625 << 626 if (verboseLevel > 3) { << 627 G4cout << "Started loop to sample broading << 628 } << 629 << 630 do 357 do 631 { 358 { 632 iteration++; 359 iteration++; 633 // Select shell based on shell occupancy 360 // Select shell based on shell occupancy 634 shellIdx = shellData->SelectRandomShell( << 361 G4int shell = shellData.SelectRandomShell(Z); 635 bindingE = shellData->BindingEnergy(Z,sh << 362 bindingE = shellData.BindingEnergy(Z,shell); 636 363 637 if (verboseLevel > 3) { << 638 G4cout << "Shell ID= " << shellIdx << 639 << " Ebind(keV)= " << bindingE/keV << << 640 } << 641 eMax = gammaEnergy0 - bindingE; 364 eMax = gammaEnergy0 - bindingE; 642 << 365 643 // Randomly sample bound electron moment 366 // Randomly sample bound electron momentum (memento: the data set is in Atomic Units) 644 G4double pSample = profileData->RandomSe << 367 G4double pSample = profileData.RandomSelectMomentum(Z,shell); 645 << 646 if (verboseLevel > 3) { << 647 G4cout << "pSample= " << pSample << G4e << 648 } << 649 // Rescale from atomic units 368 // Rescale from atomic units 650 G4double pDoppler = pSample * fine_struc 369 G4double pDoppler = pSample * fine_structure_const; 651 G4double pDoppler2 = pDoppler * pDoppler 370 G4double pDoppler2 = pDoppler * pDoppler; 652 G4double var2 = 1. + onecost * E0_m; 371 G4double var2 = 1. + onecost * E0_m; 653 G4double var3 = var2*var2 - pDoppler2; 372 G4double var3 = var2*var2 - pDoppler2; 654 G4double var4 = var2 - pDoppler2 * cosTh 373 G4double var4 = var2 - pDoppler2 * cosTheta; 655 G4double var = var4*var4 - var3 + pDoppl 374 G4double var = var4*var4 - var3 + pDoppler2 * var3; 656 if (var > 0.) 375 if (var > 0.) 657 { 376 { 658 G4double varSqrt = std::sqrt(var); 377 G4double varSqrt = std::sqrt(var); 659 G4double scale = gammaEnergy0 / var3; 378 G4double scale = gammaEnergy0 / var3; 660 // Random select either root 379 // Random select either root 661 if (G4UniformRand() < 0.5) photonE = (var4 380 if (G4UniformRand() < 0.5) photonE = (var4 - varSqrt) * scale; 662 else photonE = (var4 + varSqrt) * scale; 381 else photonE = (var4 + varSqrt) * scale; 663 } 382 } 664 else 383 else 665 { 384 { 666 photonE = -1.; 385 photonE = -1.; 667 } 386 } 668 } while ( iteration <= maxDopplerIterations 387 } while ( iteration <= maxDopplerIterations && 669 (photonE < 0. || photonE > eMax || phot 388 (photonE < 0. || photonE > eMax || photonE < eMax*G4UniformRand()) ); 670 << 389 671 // End of recalculation of photon energy wit 390 // End of recalculation of photon energy with Doppler broadening 672 // Revert to original if maximum number of i 391 // Revert to original if maximum number of iterations threshold has been reached 673 if (iteration >= maxDopplerIterations) 392 if (iteration >= maxDopplerIterations) 674 { 393 { 675 photonE = photonEoriginal; 394 photonE = photonEoriginal; 676 bindingE = 0.; 395 bindingE = 0.; 677 } 396 } 678 397 679 gammaEnergy1 = photonE; 398 gammaEnergy1 = photonE; 680 399 681 // 400 // 682 // update G4VParticleChange for the scattere 401 // update G4VParticleChange for the scattered photon 683 // 402 // >> 403 >> 404 // gammaEnergy1 = epsilon*gammaEnergy0; >> 405 >> 406 684 // New polarization 407 // New polarization >> 408 685 G4ThreeVector gammaPolarization1 = SetNewPol 409 G4ThreeVector gammaPolarization1 = SetNewPolarization(epsilon, 686 sinThetaSqr, 410 sinThetaSqr, 687 phi, 411 phi, 688 cosTheta); 412 cosTheta); 689 413 690 // Set new direction 414 // Set new direction 691 G4ThreeVector tmpDirection1( dirx,diry,dirz 415 G4ThreeVector tmpDirection1( dirx,diry,dirz ); 692 gammaDirection1 = tmpDirection1; 416 gammaDirection1 = tmpDirection1; 693 417 694 // Change reference frame. 418 // Change reference frame. >> 419 695 SystemOfRefChange(gammaDirection0,gammaDirec 420 SystemOfRefChange(gammaDirection0,gammaDirection1, 696 gammaPolarization0,gammaPolarization1) 421 gammaPolarization0,gammaPolarization1); 697 422 698 if (gammaEnergy1 > 0.) 423 if (gammaEnergy1 > 0.) 699 { 424 { 700 fParticleChange->SetProposedKineticEnerg 425 fParticleChange->SetProposedKineticEnergy( gammaEnergy1 ) ; 701 fParticleChange->ProposeMomentumDirectio 426 fParticleChange->ProposeMomentumDirection( gammaDirection1 ); 702 fParticleChange->ProposePolarization( ga 427 fParticleChange->ProposePolarization( gammaPolarization1 ); 703 } 428 } 704 else 429 else 705 { 430 { 706 gammaEnergy1 = 0.; << 707 fParticleChange->SetProposedKineticEnerg 431 fParticleChange->SetProposedKineticEnergy(0.) ; 708 fParticleChange->ProposeTrackStatus(fSto 432 fParticleChange->ProposeTrackStatus(fStopAndKill); 709 } 433 } 710 434 711 // 435 // 712 // kinematic of the scattered electron 436 // kinematic of the scattered electron 713 // 437 // >> 438 714 G4double ElecKineEnergy = gammaEnergy0 - gam 439 G4double ElecKineEnergy = gammaEnergy0 - gammaEnergy1 -bindingE; 715 440 716 // SI -protection against negative final ene << 441 // SI - Removed range test 717 // like in G4LivermoreComptonModel.cc << 442 718 if(ElecKineEnergy < 0.0) { << 719 fParticleChange->ProposeLocalEnergyDeposit << 720 return; << 721 } << 722 << 723 G4double ElecMomentum = std::sqrt(ElecKineEn 443 G4double ElecMomentum = std::sqrt(ElecKineEnergy*(ElecKineEnergy+2.*electron_mass_c2)); 724 444 725 G4ThreeVector ElecDirection((gammaEnergy0 * 445 G4ThreeVector ElecDirection((gammaEnergy0 * gammaDirection0 - 726 gammaEnergy1 * gammaDirection1) * ( 446 gammaEnergy1 * gammaDirection1) * (1./ElecMomentum)); 727 447 728 G4DynamicParticle* dp = << 729 new G4DynamicParticle (G4Electron::Electro << 730 fvect->push_back(dp); << 731 << 732 // sample deexcitation << 733 // << 734 if (verboseLevel > 3) { << 735 G4cout << "Started atomic de-excitation " << 736 } << 737 << 738 if(fAtomDeexcitation && iteration < maxDoppl << 739 G4int index = couple->GetIndex(); << 740 if(fAtomDeexcitation->CheckDeexcitationAct << 741 std::size_t nbefore = fvect->size(); << 742 G4AtomicShellEnumerator as = G4AtomicShe << 743 const G4AtomicShell* shell = fAtomDeexci << 744 fAtomDeexcitation->GenerateParticles(fve << 745 std::size_t nafter = fvect->size(); << 746 if(nafter > nbefore) { << 747 for (std::size_t i=nbefore; i<nafter; ++i) { << 748 //Check if there is enough residual energy << 749 if (bindingE >= ((*fvect)[i])->GetKineticE << 750 { << 751 //Ok, this is a valid secondary: << 752 bindingE -= ((*fvect)[i])->GetKineticE << 753 } << 754 else << 755 { << 756 //Invalid secondary: not enough energy << 757 //Keep its energy in the local deposit << 758 delete (*fvect)[i]; << 759 (*fvect)[i]=0; << 760 } << 761 } << 762 } << 763 } << 764 } << 765 //This should never happen << 766 if(bindingE < 0.0) << 767 G4Exception("G4LivermoreComptonModel::Samp << 768 "em2050",FatalException,"Negative local en << 769 << 770 fParticleChange->ProposeLocalEnergyDeposit(b 448 fParticleChange->ProposeLocalEnergyDeposit(bindingE); >> 449 >> 450 G4DynamicParticle* dp = new G4DynamicParticle (G4Electron::Electron(),ElecDirection.unit(),ElecKineEnergy) ; >> 451 fvect->push_back(dp); 771 452 772 } 453 } 773 454 774 //....oooOO0OOooo........oooOO0OOooo........oo 455 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 775 456 776 G4double G4LivermorePolarizedComptonModel::Set 457 G4double G4LivermorePolarizedComptonModel::SetPhi(G4double energyRate, 777 G4double sinSqrTh) 458 G4double sinSqrTh) 778 { 459 { 779 G4double rand1; 460 G4double rand1; 780 G4double rand2; 461 G4double rand2; 781 G4double phiProbability; 462 G4double phiProbability; 782 G4double phi; 463 G4double phi; 783 G4double a, b; 464 G4double a, b; 784 465 785 do 466 do 786 { 467 { 787 rand1 = G4UniformRand(); 468 rand1 = G4UniformRand(); 788 rand2 = G4UniformRand(); 469 rand2 = G4UniformRand(); 789 phiProbability=0.; 470 phiProbability=0.; 790 phi = twopi*rand1; 471 phi = twopi*rand1; 791 472 792 a = 2*sinSqrTh; 473 a = 2*sinSqrTh; 793 b = energyRate + 1/energyRate; 474 b = energyRate + 1/energyRate; 794 475 795 phiProbability = 1 - (a/b)*(std::cos(phi 476 phiProbability = 1 - (a/b)*(std::cos(phi)*std::cos(phi)); >> 477 >> 478 >> 479 796 } 480 } 797 while ( rand2 > phiProbability ); 481 while ( rand2 > phiProbability ); 798 return phi; 482 return phi; 799 } 483 } 800 484 801 485 802 //....oooOO0OOooo........oooOO0OOooo........oo 486 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 803 487 804 G4ThreeVector G4LivermorePolarizedComptonModel 488 G4ThreeVector G4LivermorePolarizedComptonModel::SetPerpendicularVector(G4ThreeVector& a) 805 { 489 { 806 G4double dx = a.x(); 490 G4double dx = a.x(); 807 G4double dy = a.y(); 491 G4double dy = a.y(); 808 G4double dz = a.z(); 492 G4double dz = a.z(); 809 G4double x = dx < 0.0 ? -dx : dx; 493 G4double x = dx < 0.0 ? -dx : dx; 810 G4double y = dy < 0.0 ? -dy : dy; 494 G4double y = dy < 0.0 ? -dy : dy; 811 G4double z = dz < 0.0 ? -dz : dz; 495 G4double z = dz < 0.0 ? -dz : dz; 812 if (x < y) { 496 if (x < y) { 813 return x < z ? G4ThreeVector(-dy,dx,0) : G 497 return x < z ? G4ThreeVector(-dy,dx,0) : G4ThreeVector(0,-dz,dy); 814 }else{ 498 }else{ 815 return y < z ? G4ThreeVector(dz,0,-dx) : G 499 return y < z ? G4ThreeVector(dz,0,-dx) : G4ThreeVector(-dy,dx,0); 816 } 500 } 817 } 501 } 818 502 819 //....oooOO0OOooo........oooOO0OOooo........oo 503 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 820 504 821 G4ThreeVector G4LivermorePolarizedComptonModel 505 G4ThreeVector G4LivermorePolarizedComptonModel::GetRandomPolarization(G4ThreeVector& direction0) 822 { 506 { 823 G4ThreeVector d0 = direction0.unit(); 507 G4ThreeVector d0 = direction0.unit(); 824 G4ThreeVector a1 = SetPerpendicularVector(d0 508 G4ThreeVector a1 = SetPerpendicularVector(d0); //different orthogonal 825 G4ThreeVector a0 = a1.unit(); // unit vector 509 G4ThreeVector a0 = a1.unit(); // unit vector 826 510 827 G4double rand1 = G4UniformRand(); 511 G4double rand1 = G4UniformRand(); 828 512 829 G4double angle = twopi*rand1; // random pola 513 G4double angle = twopi*rand1; // random polar angle 830 G4ThreeVector b0 = d0.cross(a0); // cross pr 514 G4ThreeVector b0 = d0.cross(a0); // cross product 831 515 832 G4ThreeVector c; 516 G4ThreeVector c; 833 517 834 c.setX(std::cos(angle)*(a0.x())+std::sin(ang 518 c.setX(std::cos(angle)*(a0.x())+std::sin(angle)*b0.x()); 835 c.setY(std::cos(angle)*(a0.y())+std::sin(ang 519 c.setY(std::cos(angle)*(a0.y())+std::sin(angle)*b0.y()); 836 c.setZ(std::cos(angle)*(a0.z())+std::sin(ang 520 c.setZ(std::cos(angle)*(a0.z())+std::sin(angle)*b0.z()); 837 521 838 G4ThreeVector c0 = c.unit(); 522 G4ThreeVector c0 = c.unit(); 839 523 840 return c0; 524 return c0; >> 525 841 } 526 } 842 527 843 //....oooOO0OOooo........oooOO0OOooo........oo 528 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 844 529 845 G4ThreeVector G4LivermorePolarizedComptonModel 530 G4ThreeVector G4LivermorePolarizedComptonModel::GetPerpendicularPolarization 846 (const G4ThreeVector& gammaDirection, const G4 531 (const G4ThreeVector& gammaDirection, const G4ThreeVector& gammaPolarization) const 847 { 532 { >> 533 848 // 534 // 849 // The polarization of a photon is always pe 535 // The polarization of a photon is always perpendicular to its momentum direction. 850 // Therefore this function removes those vec 536 // Therefore this function removes those vector component of gammaPolarization, which 851 // points in direction of gammaDirection 537 // points in direction of gammaDirection 852 // 538 // 853 // Mathematically we search the projection o 539 // Mathematically we search the projection of the vector a on the plane E, where n is the 854 // plains normal vector. 540 // plains normal vector. 855 // The basic equation can be found in each g 541 // The basic equation can be found in each geometry book (e.g. Bronstein): 856 // p = a - (a o n)/(n o n)*n 542 // p = a - (a o n)/(n o n)*n 857 543 858 return gammaPolarization - gammaPolarization 544 return gammaPolarization - gammaPolarization.dot(gammaDirection)/gammaDirection.dot(gammaDirection) * gammaDirection; 859 } 545 } 860 546 861 //....oooOO0OOooo........oooOO0OOooo........oo 547 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 862 548 863 G4ThreeVector G4LivermorePolarizedComptonModel 549 G4ThreeVector G4LivermorePolarizedComptonModel::SetNewPolarization(G4double epsilon, 864 G4double sinSqrTh, 550 G4double sinSqrTh, 865 G4double phi, 551 G4double phi, 866 G4double costheta) 552 G4double costheta) 867 { 553 { 868 G4double rand1; 554 G4double rand1; 869 G4double rand2; 555 G4double rand2; 870 G4double cosPhi = std::cos(phi); 556 G4double cosPhi = std::cos(phi); 871 G4double sinPhi = std::sin(phi); 557 G4double sinPhi = std::sin(phi); 872 G4double sinTheta = std::sqrt(sinSqrTh); 558 G4double sinTheta = std::sqrt(sinSqrTh); 873 G4double cosSqrPhi = cosPhi*cosPhi; 559 G4double cosSqrPhi = cosPhi*cosPhi; 874 // G4double cossqrth = 1.-sinSqrTh; 560 // G4double cossqrth = 1.-sinSqrTh; 875 // G4double sinsqrphi = sinPhi*sinPhi; 561 // G4double sinsqrphi = sinPhi*sinPhi; 876 G4double normalisation = std::sqrt(1. - cosS 562 G4double normalisation = std::sqrt(1. - cosSqrPhi*sinSqrTh); 877 563 >> 564 878 // Determination of Theta 565 // Determination of Theta >> 566 >> 567 // ---- MGP ---- Commented out the following 3 lines to avoid compilation >> 568 // warnings (unused variables) >> 569 // G4double thetaProbability; 879 G4double theta; 570 G4double theta; >> 571 // G4double a, b; >> 572 // G4double cosTheta; >> 573 >> 574 /* >> 575 >> 576 depaola method >> 577 >> 578 do >> 579 { >> 580 rand1 = G4UniformRand(); >> 581 rand2 = G4UniformRand(); >> 582 thetaProbability=0.; >> 583 theta = twopi*rand1; >> 584 a = 4*normalisation*normalisation; >> 585 b = (epsilon + 1/epsilon) - 2; >> 586 thetaProbability = (b + a*std::cos(theta)*std::cos(theta))/(a+b); >> 587 cosTheta = std::cos(theta); >> 588 } >> 589 while ( rand2 > thetaProbability ); >> 590 >> 591 G4double cosBeta = cosTheta; >> 592 >> 593 */ >> 594 880 595 881 // Dan Xu method (IEEE TNS, 52, 1160 (2005)) 596 // Dan Xu method (IEEE TNS, 52, 1160 (2005)) >> 597 882 rand1 = G4UniformRand(); 598 rand1 = G4UniformRand(); 883 rand2 = G4UniformRand(); 599 rand2 = G4UniformRand(); 884 600 885 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsi 601 if (rand1<(epsilon+1.0/epsilon-2)/(2.0*(epsilon+1.0/epsilon)-4.0*sinSqrTh*cosSqrPhi)) 886 { 602 { 887 if (rand2<0.5) 603 if (rand2<0.5) 888 theta = pi/2.0; 604 theta = pi/2.0; 889 else 605 else 890 theta = 3.0*pi/2.0; 606 theta = 3.0*pi/2.0; 891 } 607 } 892 else 608 else 893 { 609 { 894 if (rand2<0.5) 610 if (rand2<0.5) 895 theta = 0; 611 theta = 0; 896 else 612 else 897 theta = pi; 613 theta = pi; 898 } 614 } 899 G4double cosBeta = std::cos(theta); 615 G4double cosBeta = std::cos(theta); 900 G4double sinBeta = std::sqrt(1-cosBeta*cosBe 616 G4double sinBeta = std::sqrt(1-cosBeta*cosBeta); 901 617 902 G4ThreeVector gammaPolarization1; 618 G4ThreeVector gammaPolarization1; 903 619 904 G4double xParallel = normalisation*cosBeta; 620 G4double xParallel = normalisation*cosBeta; 905 G4double yParallel = -(sinSqrTh*cosPhi*sinPh 621 G4double yParallel = -(sinSqrTh*cosPhi*sinPhi)*cosBeta/normalisation; 906 G4double zParallel = -(costheta*sinTheta*cos 622 G4double zParallel = -(costheta*sinTheta*cosPhi)*cosBeta/normalisation; 907 G4double xPerpendicular = 0.; 623 G4double xPerpendicular = 0.; 908 G4double yPerpendicular = (costheta)*sinBeta 624 G4double yPerpendicular = (costheta)*sinBeta/normalisation; 909 G4double zPerpendicular = -(sinTheta*sinPhi) 625 G4double zPerpendicular = -(sinTheta*sinPhi)*sinBeta/normalisation; 910 626 911 G4double xTotal = (xParallel + xPerpendicula 627 G4double xTotal = (xParallel + xPerpendicular); 912 G4double yTotal = (yParallel + yPerpendicula 628 G4double yTotal = (yParallel + yPerpendicular); 913 G4double zTotal = (zParallel + zPerpendicula 629 G4double zTotal = (zParallel + zPerpendicular); 914 630 915 gammaPolarization1.setX(xTotal); 631 gammaPolarization1.setX(xTotal); 916 gammaPolarization1.setY(yTotal); 632 gammaPolarization1.setY(yTotal); 917 gammaPolarization1.setZ(zTotal); 633 gammaPolarization1.setZ(zTotal); 918 634 919 return gammaPolarization1; 635 return gammaPolarization1; >> 636 920 } 637 } 921 638 922 //....oooOO0OOooo........oooOO0OOooo........oo 639 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 923 640 924 void G4LivermorePolarizedComptonModel::SystemO 641 void G4LivermorePolarizedComptonModel::SystemOfRefChange(G4ThreeVector& direction0, 925 G4ThreeVector& direction1, 642 G4ThreeVector& direction1, 926 G4ThreeVector& polarization0, 643 G4ThreeVector& polarization0, 927 G4ThreeVector& polarization1) 644 G4ThreeVector& polarization1) 928 { 645 { 929 // direction0 is the original photon directi 646 // direction0 is the original photon direction ---> z 930 // polarization0 is the original photon pola 647 // polarization0 is the original photon polarization ---> x 931 // need to specify y axis in the real refere 648 // need to specify y axis in the real reference frame ---> y 932 G4ThreeVector Axis_Z0 = direction0.unit(); 649 G4ThreeVector Axis_Z0 = direction0.unit(); 933 G4ThreeVector Axis_X0 = polarization0.unit() 650 G4ThreeVector Axis_X0 = polarization0.unit(); 934 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_ 651 G4ThreeVector Axis_Y0 = (Axis_Z0.cross(Axis_X0)).unit(); // to be confirmed; 935 652 936 G4double direction_x = direction1.getX(); 653 G4double direction_x = direction1.getX(); 937 G4double direction_y = direction1.getY(); 654 G4double direction_y = direction1.getY(); 938 G4double direction_z = direction1.getZ(); 655 G4double direction_z = direction1.getZ(); 939 656 940 direction1 = (direction_x*Axis_X0 + directio 657 direction1 = (direction_x*Axis_X0 + direction_y*Axis_Y0 + direction_z*Axis_Z0).unit(); 941 G4double polarization_x = polarization1.getX 658 G4double polarization_x = polarization1.getX(); 942 G4double polarization_y = polarization1.getY 659 G4double polarization_y = polarization1.getY(); 943 G4double polarization_z = polarization1.getZ 660 G4double polarization_z = polarization1.getZ(); 944 661 945 polarization1 = (polarization_x*Axis_X0 + po 662 polarization1 = (polarization_x*Axis_X0 + polarization_y*Axis_Y0 + polarization_z*Axis_Z0).unit(); 946 663 947 } 664 } 948 665 >> 666 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 949 667 950 //....oooOO0OOooo........oooOO0OOooo........oo << 668 G4double G4LivermorePolarizedComptonModel::GetMeanFreePath(const G4Track& track, 951 //....oooOO0OOooo........oooOO0OOooo........oo << 669 G4double, 952 << 670 G4ForceCondition*) 953 void << 671 { 954 G4LivermorePolarizedComptonModel::InitialiseFo << 672 const G4DynamicParticle* photon = track.GetDynamicParticle(); 955 G4int Z) << 673 G4double energy = photon->GetKineticEnergy(); 956 { << 674 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple(); 957 G4AutoLock l(&LivermorePolarizedComptonModel << 675 size_t materialIndex = couple->GetIndex(); 958 if(!data[Z]) { ReadData(Z); } << 676 G4double meanFreePath; 959 l.unlock(); << 677 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex); >> 678 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX; >> 679 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex); >> 680 return meanFreePath; 960 } 681 } >> 682 >> 683 961 684