Geant4 Cross Reference |
1 // 1 // 2 // ******************************************* 2 // ******************************************************************** 3 // * License and Disclaimer 3 // * License and Disclaimer * 4 // * 4 // * * 5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. 9 // * include a list of copyright holders. * 10 // * 10 // * * 11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. * 17 // * 17 // * * 18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************* 24 // ******************************************************************** 25 // 25 // >> 26 // $Id: G4LivermorePolarizedRayleighModel.cc 93810 2015-11-02 11:27:56Z gcosmo $ 26 // 27 // 27 // Author: Sebastien Incerti 28 // Author: Sebastien Incerti 28 // 30 October 2008 29 // 30 October 2008 29 // on base of G4LowEnergyPolarizedRayl 30 // on base of G4LowEnergyPolarizedRayleigh developed by R. Capra 30 // 31 // 31 // History: 32 // History: 32 // -------- 33 // -------- 33 // 02 May 2009 S Incerti as V. Ivanchenko pr 34 // 02 May 2009 S Incerti as V. Ivanchenko proposed in G4LivermoreRayleighModel.cc 34 // 35 // 35 // Cleanup initialisation and generation of se 36 // Cleanup initialisation and generation of secondaries: 36 // - apply internal high-ener 37 // - apply internal high-energy limit only in constructor 37 // - do not apply low-energy 38 // - do not apply low-energy limit (default is 0) 38 // - remove GetMeanFreePath m 39 // - remove GetMeanFreePath method and table 39 // - remove initialisation of 40 // - remove initialisation of element selector 40 // - use G4ElementSelector 41 // - use G4ElementSelector 41 42 42 #include "G4LivermorePolarizedRayleighModel.hh 43 #include "G4LivermorePolarizedRayleighModel.hh" 43 #include "G4PhysicalConstants.hh" 44 #include "G4PhysicalConstants.hh" 44 #include "G4SystemOfUnits.hh" 45 #include "G4SystemOfUnits.hh" 45 #include "G4LogLogInterpolation.hh" 46 #include "G4LogLogInterpolation.hh" 46 #include "G4CompositeEMDataSet.hh" 47 #include "G4CompositeEMDataSet.hh" 47 #include "G4AutoLock.hh" << 48 48 49 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 50 50 51 using namespace std; 51 using namespace std; 52 namespace { G4Mutex LivermorePolarizedRayleigh << 53 52 54 //....oooOO0OOooo........oooOO0OOooo........oo 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 55 54 56 G4PhysicsFreeVector* G4LivermorePolarizedRayle << 55 G4int G4LivermorePolarizedRayleighModel::maxZ = 100; 57 G4PhysicsFreeVector* G4LivermorePolarizedRayle << 56 G4LPhysicsFreeVector* G4LivermorePolarizedRayleighModel::dataCS[] = {0}; >> 57 G4VEMDataSet* G4LivermorePolarizedRayleighModel::formFactorData = 0; 58 58 59 G4LivermorePolarizedRayleighModel::G4Livermore 59 G4LivermorePolarizedRayleighModel::G4LivermorePolarizedRayleighModel(const G4ParticleDefinition*, 60 const G4String& nam) 60 const G4String& nam) 61 :G4VEmModel(nam),fParticleChange(nullptr),is << 61 :G4VEmModel(nam),fParticleChange(0),isInitialised(false) 62 { 62 { 63 lowEnergyLimit = 250 * CLHEP::eV; << 63 fParticleChange =0; >> 64 lowEnergyLimit = 250 * eV; >> 65 //SetLowEnergyLimit(lowEnergyLimit); >> 66 //SetHighEnergyLimit(highEnergyLimit); 64 // 67 // 65 verboseLevel= 0; 68 verboseLevel= 0; 66 // Verbosity scale: 69 // Verbosity scale: 67 // 0 = nothing 70 // 0 = nothing 68 // 1 = warning for energy non-conservation 71 // 1 = warning for energy non-conservation 69 // 2 = details of energy budget 72 // 2 = details of energy budget 70 // 3 = calculation of cross sections, file o 73 // 3 = calculation of cross sections, file openings, sampling of atoms 71 // 4 = entering in methods 74 // 4 = entering in methods 72 75 73 if(verboseLevel > 0) { 76 if(verboseLevel > 0) { 74 G4cout << "Livermore Polarized Rayleigh is 77 G4cout << "Livermore Polarized Rayleigh is constructed " << G4endl 75 << "Energy range: " << LowEnergyLim << 78 << "Energy range: " 76 << HighEnergyLimit() / CLHEP::GeV << " Ge << 79 << LowEnergyLimit() / eV << " eV - " 77 << G4endl; << 80 << HighEnergyLimit() / GeV << " GeV" >> 81 << G4endl; 78 } 82 } 79 } 83 } 80 84 81 //....oooOO0OOooo........oooOO0OOooo........oo 85 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 82 86 83 G4LivermorePolarizedRayleighModel::~G4Livermor 87 G4LivermorePolarizedRayleighModel::~G4LivermorePolarizedRayleighModel() 84 { 88 { 85 if(IsMaster()) { << 89 if(IsMaster()) { 86 for(G4int i=0; i<maxZ; ++i) { << 90 for(G4int i=0; i<maxZ; ++i) { 87 if(dataCS[i]) { << 91 if(dataCS[i]) { 88 delete dataCS[i]; << 92 delete dataCS[i]; 89 dataCS[i] = nullptr; << 93 dataCS[i] = 0; 90 delete formFactorData[i]; << 94 } 91 formFactorData[i] = nullptr; << 95 } 92 } << 96 delete formFactorData; 93 } << 97 formFactorData = 0; 94 } << 98 >> 99 } 95 } 100 } 96 101 97 //....oooOO0OOooo........oooOO0OOooo........oo 102 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 98 103 99 void G4LivermorePolarizedRayleighModel::Initia 104 void G4LivermorePolarizedRayleighModel::Initialise(const G4ParticleDefinition* particle, 100 const G 105 const G4DataVector& cuts) 101 { 106 { 102 // Rayleigh process: Th << 107 // Rayleigh process: The Quantum Theory of Radiation 103 // W. << 108 // W. Heitler, Oxford at the Clarendon Press, Oxford (1954) 104 // Scattering function: A << 109 // Scattering function: A simple model of photon transport 105 // D. << 110 // D.E. Cullen, Nucl. Instr. Meth. in Phys. Res. B 101 (1995) 499-510 106 // Polarization of the outcoming photon: Be << 111 // Polarization of the outcoming photon: Beam test of a prototype detector array for the PoGO astronomical hard X-ray/soft gamma-ray polarimeter 107 // X- << 112 // T. Mizuno et al., Nucl. Instr. Meth. in Phys. Res. A 540 (2005) 158-168 108 // T. << 113 109 if (verboseLevel > 3) 114 if (verboseLevel > 3) 110 G4cout << "Calling G4LivermorePolarizedRay 115 G4cout << "Calling G4LivermorePolarizedRayleighModel::Initialise()" << G4endl; 111 116 >> 117 112 if(IsMaster()) { 118 if(IsMaster()) { 113 119 >> 120 // Form Factor >> 121 >> 122 G4VDataSetAlgorithm* ffInterpolation = new G4LogLogInterpolation; >> 123 G4String formFactorFile = "rayl/re-ff-"; >> 124 formFactorData = new G4CompositeEMDataSet(ffInterpolation,1.,1.); >> 125 formFactorData->LoadData(formFactorFile); >> 126 114 // Initialise element selector 127 // Initialise element selector 115 InitialiseElementSelectors(particle, cuts) 128 InitialiseElementSelectors(particle, cuts); 116 129 117 // Access to elements 130 // Access to elements 118 const char* path = G4FindDataDir("G4LEDATA << 131 char* path = getenv("G4LEDATA"); 119 auto elmTable = G4Element::GetElementTable << 132 G4ProductionCutsTable* theCoupleTable = 120 for (auto const & elm : *elmTable) { << 133 G4ProductionCutsTable::GetProductionCutsTable(); 121 G4int Z = std::min(elm->GetZasInt(), max << 134 G4int numOfCouples = theCoupleTable->GetTableSize(); 122 if( nullptr == dataCS[Z] ) { ReadData(Z, << 135 123 } << 136 for(G4int i=0; i<numOfCouples; ++i) >> 137 { >> 138 const G4MaterialCutsCouple* couple = >> 139 theCoupleTable->GetMaterialCutsCouple(i); >> 140 const G4Material* material = couple->GetMaterial(); >> 141 const G4ElementVector* theElementVector = material->GetElementVector(); >> 142 G4int nelm = material->GetNumberOfElements(); >> 143 >> 144 for (G4int j=0; j<nelm; ++j) >> 145 { >> 146 G4int Z = G4lrint((*theElementVector)[j]->GetZ()); >> 147 if(Z < 1) { Z = 1; } >> 148 else if(Z > maxZ) { Z = maxZ; } >> 149 if( (!dataCS[Z]) ) { ReadData(Z, path); } >> 150 } >> 151 } 124 } 152 } 125 153 126 if(isInitialised) { return; } 154 if(isInitialised) { return; } 127 fParticleChange = GetParticleChangeForGamma( 155 fParticleChange = GetParticleChangeForGamma(); 128 isInitialised = true; 156 isInitialised = true; >> 157 129 } 158 } 130 159 131 160 132 //....oooOO0OOooo........oooOO0OOooo........oo 161 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 133 162 134 void G4LivermorePolarizedRayleighModel::Initia << 163 void G4LivermorePolarizedRayleighModel::InitialiseLocal(const G4ParticleDefinition*, 135 const G4ParticleDefinition*, G << 164 G4VEmModel* masterModel) 136 { 165 { 137 SetElementSelectors(masterModel->GetElementS 166 SetElementSelectors(masterModel->GetElementSelectors()); 138 } 167 } 139 168 140 //....oooOO0OOooo........oooOO0OOooo........oo 169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 141 170 142 void G4LivermorePolarizedRayleighModel::ReadDa << 171 void G4LivermorePolarizedRayleighModel::ReadData(size_t Z, const char* path) 143 { 172 { 144 if (verboseLevel > 1) { << 173 if (verboseLevel > 1) 145 G4cout << "Calling ReadData() of G4Livermo << 174 { 146 << G4endl; << 175 G4cout << "Calling ReadData() of G4LivermoreRayleighModel" 147 } << 176 << G4endl; >> 177 } 148 178 149 if(nullptr != dataCS[Z]) { return; } << 179 if(dataCS[Z]) { return; } 150 180 151 const char* datadir = path; 181 const char* datadir = path; 152 182 153 if(nullptr == datadir) << 183 if(!datadir) 154 { 184 { 155 datadir = G4FindDataDir("G4LEDATA"); << 185 datadir = getenv("G4LEDATA"); 156 if(nullptr == datadir) << 186 if(!datadir) 157 { 187 { 158 G4Exception("G4LivermoreRayleighModelModel 188 G4Exception("G4LivermoreRayleighModelModel::ReadData()","em0006", 159 FatalException, 189 FatalException, 160 "Environment variable G4LEDATA not d 190 "Environment variable G4LEDATA not defined"); 161 return; 191 return; 162 } 192 } 163 } 193 } 164 dataCS[Z] = new G4PhysicsFreeVector(); << 194 165 formFactorData[Z] = new G4PhysicsFreeVector( << 195 // >> 196 >> 197 dataCS[Z] = new G4LPhysicsFreeVector(); >> 198 >> 199 // Activation of spline interpolation >> 200 //dataCS[Z] ->SetSpline(true); 166 201 167 std::ostringstream ostCS; 202 std::ostringstream ostCS; 168 ostCS << datadir << "/livermore/rayl/re-cs-" 203 ostCS << datadir << "/livermore/rayl/re-cs-" << Z <<".dat"; 169 std::ifstream finCS(ostCS.str().c_str()); 204 std::ifstream finCS(ostCS.str().c_str()); 170 205 171 if( !finCS .is_open() ) { << 206 if( !finCS .is_open() ) 172 G4ExceptionDescription ed; << 207 { 173 ed << "G4LivermorePolarizedRayleighModel d << 208 G4ExceptionDescription ed; 174 << "> is not opened!" << G4endl; << 209 ed << "G4LivermorePolarizedRayleighModel data file <" << ostCS.str().c_str() 175 G4Exception("G4LivermorePolarizedRayleighM << 210 << "> is not opened!" << G4endl; 176 FatalException, << 211 G4Exception("G4LivermorePolarizedRayleighModel::ReadData()","em0003",FatalException, 177 ed,"G4LEDATA version should be G4EMLOW8.0 << 212 ed,"G4LEDATA version should be G4EMLOW6.27 or later."); 178 return; << 213 return; 179 } else { << 214 } 180 if(verboseLevel > 3) { << 215 else 181 G4cout << "File " << ostCS.str() << 216 { 182 << " is opened by G4LivermoreRayleighMo << 217 if(verboseLevel > 3) { 183 } << 218 G4cout << "File " << ostCS.str() 184 dataCS[Z]->Retrieve(finCS, true); << 219 << " is opened by G4LivermoreRayleighModel" << G4endl; 185 } << 220 } 186 << 221 dataCS[Z]->Retrieve(finCS, true); 187 std::ostringstream ostFF; << 222 } 188 ostFF << datadir << "/livermore/rayl/re-ff-" << 223 } 189 std::ifstream finFF(ostFF.str().c_str()); << 190 << 191 if( !finFF.is_open() ) { << 192 G4ExceptionDescription ed; << 193 ed << "G4LivermorePolarizedRayleighModel d << 194 << "> is not opened!" << G4endl; << 195 G4Exception("G4LivermorePolarizedRayleighM << 196 FatalException, << 197 ed,"G4LEDATA version should be G4EMLOW8.0 << 198 return; << 199 } else { << 200 if(verboseLevel > 3) { << 201 G4cout << "File " << ostFF.str() << 202 << " is opened by G4LivermoreRa << 203 } << 204 formFactorData[Z]->Retrieve(finFF, true); << 205 } << 206 } << 207 224 208 //....oooOO0OOooo........oooOO0OOooo........oo << 225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 209 226 210 G4double G4LivermorePolarizedRayleighModel::Co << 227 G4double G4LivermorePolarizedRayleighModel::ComputeCrossSectionPerAtom( 211 const 228 const G4ParticleDefinition*, 212 G4double GammaEnergy, 229 G4double GammaEnergy, 213 G4double Z, G4double, 230 G4double Z, G4double, 214 G4double, G4double) 231 G4double, G4double) 215 { << 232 { 216 if (verboseLevel > 1) { << 233 if (verboseLevel > 1) 217 G4cout << "G4LivermoreRayleighModel::Compu << 234 { 218 << G4endl; << 235 G4cout << "G4LivermoreRayleighModel::ComputeCrossSectionPerAtom()" 219 } << 236 << G4endl; >> 237 } 220 238 221 if(GammaEnergy < lowEnergyLimit) { return 0. << 239 if(GammaEnergy < lowEnergyLimit) { return 0.0; } >> 240 >> 241 G4double xs = 0.0; 222 242 223 G4double xs = 0.0; << 243 G4int intZ = G4lrint(Z); 224 244 225 G4int intZ = G4lrint(Z); << 245 if(intZ < 1 || intZ > maxZ) { return xs; } 226 if(intZ < 1 || intZ > maxZ) { return xs; } << 227 246 228 G4PhysicsFreeVector* pv = dataCS[intZ]; << 247 G4LPhysicsFreeVector* pv = dataCS[intZ]; 229 248 230 // if element was not initialised << 249 // if element was not initialised 231 // do initialisation safely for MT mode << 250 // do initialisation safely for MT mode 232 if(nullptr == pv) { << 251 if(!pv) { 233 InitialiseForElement(0, intZ); << 252 InitialiseForElement(0, intZ); 234 pv = dataCS[intZ]; << 253 pv = dataCS[intZ]; 235 if(nullptr == pv) { return xs; } << 254 if(!pv) { return xs; } 236 } << 255 } >> 256 >> 257 G4int n = pv->GetVectorLength() - 1; >> 258 G4double e = GammaEnergy/MeV; >> 259 if(e >= pv->Energy(n)) { >> 260 xs = (*pv)[n]/(e*e); >> 261 } else if(e >= pv->Energy(0)) { >> 262 xs = pv->Value(e)/(e*e); >> 263 } 237 264 238 G4int n = G4int(pv->GetVectorLength() - 1); << 265 /* if(verboseLevel > 0) 239 G4double e = GammaEnergy/MeV; << 266 { 240 if(e >= pv->Energy(n)) { << 267 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)=" 241 xs = (*pv)[n]/(e*e); << 268 << e << G4endl; 242 } else if(e >= pv->Energy(0)) { << 269 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl; 243 xs = pv->Value(e)/(e*e); << 270 G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0] 244 } << 271 << G4endl; 245 return xs; << 272 G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n] 246 } << 273 << G4endl; >> 274 G4cout << "*********************************************************" >> 275 << G4endl; >> 276 } >> 277 */ >> 278 >> 279 return xs; >> 280 } 247 281 248 //....oooOO0OOooo........oooOO0OOooo........oo 282 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 249 283 250 void G4LivermorePolarizedRayleighModel::Sample << 284 void G4LivermorePolarizedRayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/, 251 std::vector<G4DynamicParticle* << 285 const G4MaterialCutsCouple* couple, 252 const G4MaterialCutsCouple* couple, << 286 const G4DynamicParticle* aDynamicGamma, 253 const G4DynamicParticle* aDynamicGamma, << 287 G4double, 254 G4double, G4double) << 288 G4double) 255 { 289 { 256 if (verboseLevel > 3) 290 if (verboseLevel > 3) 257 G4cout << "Calling SampleSecondaries() of 291 G4cout << "Calling SampleSecondaries() of G4LivermorePolarizedRayleighModel" << G4endl; 258 292 259 G4double photonEnergy0 = aDynamicGamma->GetK 293 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy(); 260 294 261 if (photonEnergy0 <= lowEnergyLimit) 295 if (photonEnergy0 <= lowEnergyLimit) 262 { 296 { 263 fParticleChange->ProposeTrackStatus(fStopA << 297 fParticleChange->ProposeTrackStatus(fStopAndKill); 264 fParticleChange->SetProposedKineticEnergy( << 298 fParticleChange->SetProposedKineticEnergy(0.); 265 fParticleChange->ProposeLocalEnergyDeposit << 299 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0); 266 return; << 300 return ; 267 } 301 } 268 302 269 G4ParticleMomentum photonDirection0 = aDynam 303 G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection(); 270 304 271 // Select randomly one element in the curren 305 // Select randomly one element in the current material >> 306 // G4int Z = crossSectionHandler->SelectRandomAtom(couple,photonEnergy0); 272 const G4ParticleDefinition* particle = aDyn 307 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition(); 273 const G4Element* elm = SelectRandomAtom(coup 308 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0); 274 G4int Z = elm->GetZasInt(); << 309 G4int Z = (G4int)elm->GetZ(); 275 310 276 G4double outcomingPhotonCosTheta = GenerateC 311 G4double outcomingPhotonCosTheta = GenerateCosTheta(photonEnergy0, Z); 277 G4double outcomingPhotonPhi = GeneratePhi(ou 312 G4double outcomingPhotonPhi = GeneratePhi(outcomingPhotonCosTheta); 278 G4double beta = GeneratePolarizationAngle(); << 313 G4double beta=GeneratePolarizationAngle(); 279 314 280 // incomingPhoton reference frame: 315 // incomingPhoton reference frame: 281 // z = versor parallel to the incomingPhoton 316 // z = versor parallel to the incomingPhotonDirection 282 // x = versor parallel to the incomingPhoton 317 // x = versor parallel to the incomingPhotonPolarization 283 // y = defined as z^x 318 // y = defined as z^x 284 319 285 // outgoingPhoton reference frame: 320 // outgoingPhoton reference frame: 286 // z' = versor parallel to the outgoingPhoto 321 // z' = versor parallel to the outgoingPhotonDirection 287 // x' = defined as x-x*z'z' normalized 322 // x' = defined as x-x*z'z' normalized 288 // y' = defined as z'^x' << 323 // y' = defined as z'^x' >> 324 289 G4ThreeVector z(aDynamicGamma->GetMomentumDi 325 G4ThreeVector z(aDynamicGamma->GetMomentumDirection().unit()); 290 G4ThreeVector x(GetPhotonPolarization(*aDyna 326 G4ThreeVector x(GetPhotonPolarization(*aDynamicGamma)); 291 G4ThreeVector y(z.cross(x)); 327 G4ThreeVector y(z.cross(x)); 292 328 293 // z' = std::cos(phi)*std::sin(theta) << 329 // z' = std::cos(phi)*std::sin(theta) x + std::sin(phi)*std::sin(theta) y + std::cos(theta) z 294 // x + std::sin(phi)*std::sin(theta) y + std << 295 G4double xDir; 330 G4double xDir; 296 G4double yDir; 331 G4double yDir; 297 G4double zDir; 332 G4double zDir; 298 zDir=outcomingPhotonCosTheta; 333 zDir=outcomingPhotonCosTheta; 299 xDir=std::sqrt(1-outcomingPhotonCosTheta*out 334 xDir=std::sqrt(1-outcomingPhotonCosTheta*outcomingPhotonCosTheta); 300 yDir=xDir; 335 yDir=xDir; 301 xDir*=std::cos(outcomingPhotonPhi); 336 xDir*=std::cos(outcomingPhotonPhi); 302 yDir*=std::sin(outcomingPhotonPhi); 337 yDir*=std::sin(outcomingPhotonPhi); 303 338 304 G4ThreeVector zPrime((xDir*x + yDir*y + zDir 339 G4ThreeVector zPrime((xDir*x + yDir*y + zDir*z).unit()); 305 G4ThreeVector xPrime(x.perpPart(zPrime).unit 340 G4ThreeVector xPrime(x.perpPart(zPrime).unit()); 306 G4ThreeVector yPrime(zPrime.cross(xPrime)); 341 G4ThreeVector yPrime(zPrime.cross(xPrime)); 307 342 308 // outgoingPhotonPolarization is directed as << 343 // outgoingPhotonPolarization is directed as x' std::cos(beta) + y' std::sin(beta) 309 // x' std::cos(beta) + y' std::sin(beta) << 310 G4ThreeVector outcomingPhotonPolarization(xP 344 G4ThreeVector outcomingPhotonPolarization(xPrime*std::cos(beta) + yPrime*std::sin(beta)); 311 345 312 fParticleChange->ProposeMomentumDirection(zP 346 fParticleChange->ProposeMomentumDirection(zPrime); 313 fParticleChange->ProposePolarization(outcomi 347 fParticleChange->ProposePolarization(outcomingPhotonPolarization); 314 fParticleChange->SetProposedKineticEnergy(ph 348 fParticleChange->SetProposedKineticEnergy(photonEnergy0); >> 349 315 } 350 } 316 351 317 //....oooOO0OOooo........oooOO0OOooo........oo 352 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 318 353 319 G4double G4LivermorePolarizedRayleighModel::Ge << 354 G4double G4LivermorePolarizedRayleighModel::GenerateCosTheta(G4double incomingPhotonEnergy, G4int zAtom) const 320 { 355 { 321 // d sigma 356 // d sigma k0 322 // --------- = r0^2 * pi * F^2(x, Z) * ( 2 357 // --------- = r0^2 * pi * F^2(x, Z) * ( 2 - sin^2 theta) * std::sin (theta), x = ---- std::sin(theta/2) 323 // d theta 358 // d theta hc 324 359 325 // d sigma 360 // d sigma k0 1 - y 326 // --------- = r0^2 * pi * F^2(x, Z) * ( 1 + 361 // --------- = r0^2 * pi * F^2(x, Z) * ( 1 + y^2), x = ---- std::sqrt ( ------- ), y = std::cos(theta) 327 // d y 362 // d y hc 2 328 363 329 // Z 364 // Z 330 // F(x, Z) ~ -------- 365 // F(x, Z) ~ -------- 331 // a + bx 366 // a + bx 332 // 367 // 333 // The time to exit from the outer loop grow 368 // The time to exit from the outer loop grows as ~ k0 334 // On pcgeant2 the time is ~ 1 s for k0 ~ 1 369 // On pcgeant2 the time is ~ 1 s for k0 ~ 1 MeV on the oxygen element. A 100 GeV 335 // event will take ~ 10 hours. 370 // event will take ~ 10 hours. 336 // 371 // 337 // On the avarage the inner loop does 1.5 it 372 // On the avarage the inner loop does 1.5 iterations before exiting 338 const G4double xxfact = CLHEP::cm/(CLHEP::h_ << 373 339 const G4double xFactor = incomingPhotonEnerg << 374 const G4double xFactor = (incomingPhotonEnergy*cm)/(h_Planck*c_light); >> 375 //const G4VEMDataSet * formFactorData = GetScatterFunctionData(); 340 376 341 G4double cosTheta; 377 G4double cosTheta; 342 G4double fCosTheta; 378 G4double fCosTheta; 343 G4double x; 379 G4double x; 344 G4double fValue; 380 G4double fValue; 345 381 346 if (incomingPhotonEnergy > 5.*CLHEP::MeV) << 382 if (incomingPhotonEnergy > 5.*MeV) 347 { 383 { 348 cosTheta = 1.; 384 cosTheta = 1.; 349 } 385 } 350 else 386 else 351 { 387 { 352 do 388 do 353 { 389 { 354 do 390 do 355 { 391 { 356 cosTheta = 2.*G4UniformRand()-1.; 392 cosTheta = 2.*G4UniformRand()-1.; 357 fCosTheta = (1.+cosTheta*cosTheta)/2.; 393 fCosTheta = (1.+cosTheta*cosTheta)/2.; 358 } 394 } 359 while (fCosTheta < G4UniformRand()); 395 while (fCosTheta < G4UniformRand()); 360 396 361 x = xFactor*std::sqrt((1.-cosTheta)/2.); 397 x = xFactor*std::sqrt((1.-cosTheta)/2.); 362 398 363 if (x > 1.e+005) 399 if (x > 1.e+005) 364 fValue = formFactorData[Z]->Value(x); << 400 fValue = formFactorData->FindValue(x, zAtom-1); 365 else 401 else 366 fValue = formFactorData[Z]->Value(0.); << 402 fValue = formFactorData->FindValue(0., zAtom-1); 367 403 368 fValue /= Z; << 404 fValue/=zAtom; 369 fValue *= fValue; << 405 fValue*=fValue; 370 } 406 } 371 while(fValue < G4UniformRand()); 407 while(fValue < G4UniformRand()); 372 } 408 } 373 409 374 return cosTheta; 410 return cosTheta; 375 } 411 } 376 412 377 //....oooOO0OOooo........oooOO0OOooo........oo 413 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 378 414 379 G4double G4LivermorePolarizedRayleighModel::Ge 415 G4double G4LivermorePolarizedRayleighModel::GeneratePhi(G4double cosTheta) const 380 { 416 { 381 // d sigma 417 // d sigma 382 // --------- = alpha * ( 1 - sin^2 (theta) * 418 // --------- = alpha * ( 1 - sin^2 (theta) * cos^2 (phi) ) 383 // d phi 419 // d phi 384 420 385 // On the average the loop takes no more tha 421 // On the average the loop takes no more than 2 iterations before exiting 386 422 387 G4double phi; 423 G4double phi; 388 G4double cosPhi; 424 G4double cosPhi; 389 G4double phiProbability; 425 G4double phiProbability; 390 G4double sin2Theta; 426 G4double sin2Theta; 391 427 392 sin2Theta=1.-cosTheta*cosTheta; 428 sin2Theta=1.-cosTheta*cosTheta; 393 429 394 do 430 do 395 { 431 { 396 phi = CLHEP::twopi * G4UniformRand(); << 432 phi = twopi * G4UniformRand(); 397 cosPhi = std::cos(phi); 433 cosPhi = std::cos(phi); 398 phiProbability= 1. - sin2Theta*cosPhi*co 434 phiProbability= 1. - sin2Theta*cosPhi*cosPhi; 399 } 435 } 400 while (phiProbability < G4UniformRand()); 436 while (phiProbability < G4UniformRand()); 401 437 402 return phi; 438 return phi; 403 } 439 } 404 440 405 //....oooOO0OOooo........oooOO0OOooo........oo 441 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 406 442 407 G4double G4LivermorePolarizedRayleighModel::Ge 443 G4double G4LivermorePolarizedRayleighModel::GeneratePolarizationAngle(void) const 408 { 444 { 409 // Rayleigh polarization is always on the x' 445 // Rayleigh polarization is always on the x' direction >> 446 410 return 0; 447 return 0; 411 } 448 } 412 449 413 //....oooOO0OOooo........oooOO0OOooo........oo 450 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 414 451 415 G4ThreeVector G4LivermorePolarizedRayleighMode 452 G4ThreeVector G4LivermorePolarizedRayleighModel::GetPhotonPolarization(const G4DynamicParticle& photon) 416 { 453 { 417 // From G4VLowEnergyDiscretePhotonProcess.cc << 454 >> 455 // SI - From G4VLowEnergyDiscretePhotonProcess.cc >> 456 418 G4ThreeVector photonMomentumDirection; 457 G4ThreeVector photonMomentumDirection; 419 G4ThreeVector photonPolarization; 458 G4ThreeVector photonPolarization; 420 459 421 photonPolarization = photon.GetPolarization( 460 photonPolarization = photon.GetPolarization(); 422 photonMomentumDirection = photon.GetMomentum 461 photonMomentumDirection = photon.GetMomentumDirection(); 423 462 424 if ((!photonPolarization.isOrthogonal(photon 463 if ((!photonPolarization.isOrthogonal(photonMomentumDirection, 1e-6)) || photonPolarization.mag()==0.) 425 { 464 { 426 // if |photonPolarization|==0. or |photo 465 // if |photonPolarization|==0. or |photonPolarization * photonDirection0| > 1e-6 * |photonPolarization ^ photonDirection0| 427 // then polarization is choosen randomly << 466 // then polarization is choosen randomly. >> 467 428 G4ThreeVector e1(photonMomentumDirection 468 G4ThreeVector e1(photonMomentumDirection.orthogonal().unit()); 429 G4ThreeVector e2(photonMomentumDirection 469 G4ThreeVector e2(photonMomentumDirection.cross(e1).unit()); 430 470 431 G4double angle(G4UniformRand() * CLHEP:: << 471 G4double angle(G4UniformRand() * twopi); 432 472 433 e1*=std::cos(angle); 473 e1*=std::cos(angle); 434 e2*=std::sin(angle); 474 e2*=std::sin(angle); 435 475 436 photonPolarization=e1+e2; 476 photonPolarization=e1+e2; 437 } 477 } 438 else if (photonPolarization.howOrthogonal(ph 478 else if (photonPolarization.howOrthogonal(photonMomentumDirection) != 0.) 439 { 479 { 440 // if |photonPolarization * photonDirect 480 // if |photonPolarization * photonDirection0| != 0. 441 // then polarization is made orthonormal << 481 // then polarization is made orthonormal; >> 482 442 photonPolarization=photonPolarization.pe 483 photonPolarization=photonPolarization.perpPart(photonMomentumDirection); 443 } 484 } 444 485 445 return photonPolarization.unit(); 486 return photonPolarization.unit(); 446 } 487 } 447 488 448 //....oooOO0OOooo........oooOO0OOooo........oo 489 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 449 490 450 void G4LivermorePolarizedRayleighModel::Initia << 491 #include "G4AutoLock.hh" 451 const G4ParticleDefinition*, G << 492 namespace { G4Mutex LivermorePolarizedRayleighModelMutex = G4MUTEX_INITIALIZER; } 452 { << 493 453 G4AutoLock l(&LivermorePolarizedRayleighMode << 494 void G4LivermorePolarizedRayleighModel::InitialiseForElement(const G4ParticleDefinition*, 454 if(nullptr == dataCS[Z]) { ReadData(Z); } << 495 G4int Z) 455 l.unlock(); << 496 { 456 } << 497 G4AutoLock l(&LivermorePolarizedRayleighModelMutex); >> 498 // G4cout << "G4LivermoreRayleighModel::InitialiseForElement Z= " >> 499 // << Z << G4endl; >> 500 if(!dataCS[Z]) { ReadData(Z); } >> 501 l.unlock(); >> 502 } 457 503