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