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