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 // Author: Sebastien Incerti 26 // Author: Sebastien Incerti 27 // 31 March 2012 27 // 31 March 2012 28 // on base of G4LivermoreRayleighModel 28 // on base of G4LivermoreRayleighModel 29 // 29 // 30 30 31 #include "G4LivermoreRayleighModel.hh" 31 #include "G4LivermoreRayleighModel.hh" 32 << 33 #include "G4AutoLock.hh" << 34 #include "G4EmParameters.hh" << 35 #include "G4RayleighAngularGenerator.hh" << 36 #include "G4SystemOfUnits.hh" 32 #include "G4SystemOfUnits.hh" >> 33 #include "G4RayleighAngularGenerator.hh" 37 34 38 //....oooOO0OOooo........oooOO0OOooo........oo << 35 using namespace std; 39 << 40 namespace << 41 { << 42 G4Mutex LivermoreRayleighModelMutex = G4MUTEX_ << 43 } << 44 << 45 G4PhysicsFreeVector* G4LivermoreRayleighModel: << 46 G4String G4LivermoreRayleighModel::gDataDirect << 47 36 48 //....oooOO0OOooo........oooOO0OOooo........oo 37 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 49 39 50 G4LivermoreRayleighModel::G4LivermoreRayleighM << 40 G4LivermoreRayleighModel::G4LivermoreRayleighModel() >> 41 :G4VEmModel("LivermoreRayleigh"),isInitialised(false),maxZ(100) 51 { 42 { 52 fParticleChange = nullptr; << 43 fParticleChange = 0; 53 lowEnergyLimit = 10 * CLHEP::eV; << 44 lowEnergyLimit = 10 * eV; >> 45 >> 46 dataCS.resize(maxZ+1,0); 54 47 55 SetAngularDistribution(new G4RayleighAngular 48 SetAngularDistribution(new G4RayleighAngularGenerator()); 56 << 49 57 verboseLevel = 0; << 50 verboseLevel= 0; 58 // Verbosity scale for debugging purposes: 51 // Verbosity scale for debugging purposes: 59 // 0 = nothing << 52 // 0 = nothing 60 // 1 = calculation of cross sections, file o 53 // 1 = calculation of cross sections, file openings... 61 // 2 = entering in methods 54 // 2 = entering in methods 62 55 63 if (verboseLevel > 0) { << 56 if(verboseLevel > 0) >> 57 { 64 G4cout << "G4LivermoreRayleighModel is con 58 G4cout << "G4LivermoreRayleighModel is constructed " << G4endl; 65 } 59 } >> 60 66 } 61 } 67 62 68 //....oooOO0OOooo........oooOO0OOooo........oo 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 69 64 70 G4LivermoreRayleighModel::~G4LivermoreRayleigh 65 G4LivermoreRayleighModel::~G4LivermoreRayleighModel() 71 { << 66 { 72 if (IsMaster()) { << 67 for(G4int i=0; i<=maxZ; ++i) { delete dataCS[i]; } 73 for (G4int i = 0; i <= maxZ; ++i) { << 74 if (nullptr != dataCS[i]) { << 75 delete dataCS[i]; << 76 dataCS[i] = nullptr; << 77 } << 78 } << 79 } << 80 } 68 } 81 69 82 //....oooOO0OOooo........oooOO0OOooo........oo 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 83 71 84 void G4LivermoreRayleighModel::Initialise(cons 72 void G4LivermoreRayleighModel::Initialise(const G4ParticleDefinition* particle, 85 cons << 73 const G4DataVector& cuts) 86 { 74 { 87 if (verboseLevel > 1) { << 75 if (verboseLevel > 1) >> 76 { 88 G4cout << "Calling Initialise() of G4Liver 77 G4cout << "Calling Initialise() of G4LivermoreRayleighModel." << G4endl 89 << "Energy range: " << LowEnergyLim << 78 << "Energy range: " 90 << " GeV" << G4endl; << 79 << LowEnergyLimit() / eV << " eV - " 91 } << 80 << HighEnergyLimit() / GeV << " GeV" 92 << 81 << G4endl; 93 if (IsMaster()) { << 82 } 94 // Initialise element selector << 83 95 InitialiseElementSelectors(particle, cuts) << 84 // Initialise element selector 96 << 85 InitialiseElementSelectors(particle, cuts); 97 // Access to elements << 86 98 const G4ElementTable* elemTable = G4Elemen << 87 // Access to elements 99 std::size_t numElems = (*elemTable).size() << 88 100 for (std::size_t ie = 0; ie < numElems; ++ << 89 char* path = getenv("G4LEDATA"); 101 const G4Element* elem = (*elemTable)[ie] << 90 G4ProductionCutsTable* theCoupleTable = 102 const G4int Z = std::min(maxZ, elem->Get << 91 G4ProductionCutsTable::GetProductionCutsTable(); 103 if (dataCS[Z] == nullptr) { << 92 G4int numOfCouples = theCoupleTable->GetTableSize(); 104 ReadData(Z); << 93 105 } << 94 for(G4int i=0; i<numOfCouples; ++i) 106 } << 95 { 107 } << 96 const G4MaterialCutsCouple* couple = 108 if (isInitialised) { << 97 theCoupleTable->GetMaterialCutsCouple(i); 109 return; << 98 const G4Material* material = couple->GetMaterial(); 110 } << 99 const G4ElementVector* theElementVector = material->GetElementVector(); >> 100 G4int nelm = material->GetNumberOfElements(); >> 101 >> 102 for (G4int j=0; j<nelm; ++j) >> 103 { >> 104 G4int Z = G4lrint((*theElementVector)[j]->GetZ()); >> 105 if(Z < 1) { Z = 1; } >> 106 else if(Z > maxZ) { Z = maxZ; } >> 107 if( (!dataCS[Z]) ) { ReadData(Z, path); } >> 108 } >> 109 } >> 110 // >> 111 >> 112 //fGenerator->PrintGeneratorInformation(); >> 113 >> 114 if(isInitialised) { return; } 111 fParticleChange = GetParticleChangeForGamma( 115 fParticleChange = GetParticleChangeForGamma(); 112 isInitialised = true; 116 isInitialised = true; 113 } << 114 117 115 //....oooOO0OOooo........oooOO0OOooo........oo << 116 << 117 void G4LivermoreRayleighModel::InitialiseLocal << 118 { << 119 SetElementSelectors(masterModel->GetElementS << 120 } 118 } 121 119 122 //....oooOO0OOooo........oooOO0OOooo........oo 120 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 123 121 124 const G4String& G4LivermoreRayleighModel::Find << 122 void G4LivermoreRayleighModel::ReadData(size_t Z, const char* path) 125 { 123 { 126 // no check in this method - environment var << 124 if (verboseLevel > 1) 127 if (gDataDirectory.empty()) { << 125 { 128 auto param = G4EmParameters::Instance(); << 126 G4cout << "Calling ReadData() of G4LivermoreRayleighModel" 129 std::ostringstream ost; << 127 << G4endl; 130 if (param->LivermoreDataDir() == "livermor << 128 } 131 ost << param->GetDirLEDATA() << "/liverm << 129 >> 130 if(dataCS[Z]) { return; } >> 131 >> 132 const char* datadir = path; >> 133 >> 134 if(!datadir) >> 135 { >> 136 datadir = getenv("G4LEDATA"); >> 137 if(!datadir) >> 138 { >> 139 G4Exception("G4LivermoreRayleighModelModel::ReadData()","em0006", >> 140 FatalException, >> 141 "Environment variable G4LEDATA not defined"); >> 142 return; 132 } 143 } 133 else { << 134 ost << param->GetDirLEDATA() << "/epics2 << 135 } << 136 gDataDirectory = ost.str(); << 137 } << 138 return gDataDirectory; << 139 } << 140 << 141 //....oooOO0OOooo........oooOO0OOooo........oo << 142 << 143 void G4LivermoreRayleighModel::ReadData(const << 144 { << 145 if (verboseLevel > 1) { << 146 G4cout << "Calling ReadData() of G4Livermo << 147 } << 148 const G4int Z = std::min(ZZ, maxZ); << 149 << 150 if (nullptr != dataCS[Z]) { << 151 return; << 152 } 144 } 153 145 154 dataCS[Z] = new G4PhysicsFreeVector(); << 146 // 155 << 147 >> 148 dataCS[Z] = new G4LPhysicsFreeVector(); >> 149 >> 150 // Activation of spline interpolation >> 151 //dataCS[Z] ->SetSpline(true); >> 152 156 std::ostringstream ostCS; 153 std::ostringstream ostCS; 157 ostCS << FindDirectoryPath() << "re-cs-" << << 154 ostCS << datadir << "/livermore/rayl/re-cs-" << Z <<".dat"; 158 << 159 std::ifstream finCS(ostCS.str().c_str()); 155 std::ifstream finCS(ostCS.str().c_str()); 160 << 156 161 if (!finCS.is_open()) { << 157 if( !finCS .is_open() ) >> 158 { 162 G4ExceptionDescription ed; 159 G4ExceptionDescription ed; 163 ed << "G4LivermoreRayleighModel data file << 160 ed << "G4LivermoreRayleighModel data file <" << ostCS.str().c_str() 164 << G4endl; << 161 << "> is not opened!" << G4endl; 165 G4Exception("G4LivermoreRayleighModel::Rea << 162 G4Exception("G4LivermoreRayleighModel::ReadData()","em0003",FatalException, 166 "G4LEDATA version should be G4 << 163 ed,"G4LEDATA version should be G4EMLOW6.27 or later."); 167 return; 164 return; 168 } << 165 } 169 else { << 166 else 170 if (verboseLevel > 3) { << 167 { 171 G4cout << "File " << ostCS.str() << " is << 168 if(verboseLevel > 3) { >> 169 G4cout << "File " << ostCS.str() >> 170 << " is opened by G4LivermoreRayleighModel" << G4endl; 172 } 171 } 173 dataCS[Z]->Retrieve(finCS, true); 172 dataCS[Z]->Retrieve(finCS, true); 174 } << 173 } 175 } 174 } 176 175 177 //....oooOO0OOooo........oooOO0OOooo........oo 176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 178 177 179 G4double G4LivermoreRayleighModel::ComputeCros << 178 G4double G4LivermoreRayleighModel::ComputeCrossSectionPerAtom( 180 << 179 const G4ParticleDefinition*, 181 << 180 G4double GammaEnergy, >> 181 G4double Z, G4double, >> 182 G4double, G4double) 182 { 183 { 183 if (verboseLevel > 1) { << 184 if (verboseLevel > 1) 184 G4cout << "G4LivermoreRayleighModel::Compu << 185 { 185 } << 186 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreRayleighModel" 186 << 187 << G4endl; 187 if (GammaEnergy < lowEnergyLimit) { << 188 return 0.0; << 189 } 188 } 190 189 >> 190 if(GammaEnergy < lowEnergyLimit) { return 0.0; } >> 191 191 G4double xs = 0.0; 192 G4double xs = 0.0; >> 193 192 G4int intZ = G4lrint(Z); 194 G4int intZ = G4lrint(Z); 193 if (intZ < 1 || intZ > maxZ) { << 194 return xs; << 195 } << 196 195 197 G4PhysicsFreeVector* pv = dataCS[intZ]; << 196 if(intZ < 1 || intZ > maxZ) { return xs; } 198 197 199 // if element was not initialised << 198 G4LPhysicsFreeVector* pv = dataCS[intZ]; 200 // do initialisation safely for MT mode << 201 if (nullptr == pv) { << 202 InitialiseForElement(nullptr, intZ); << 203 pv = dataCS[intZ]; << 204 if (nullptr == pv) { << 205 return xs; << 206 } << 207 } << 208 199 209 auto n = G4int(pv->GetVectorLength() - 1); << 200 // element was not initialised 210 G4double e = GammaEnergy / MeV; << 201 if(!pv) 211 if (e >= pv->Energy(n)) { << 202 { 212 xs = (*pv)[n] / (e * e); << 203 char* path = getenv("G4LEDATA"); 213 } << 204 ReadData(intZ, path); 214 else if (e >= pv->Energy(0)) { << 205 pv = dataCS[intZ]; 215 xs = pv->Value(e) / (e * e); << 206 if(!pv) { return xs; } 216 } 207 } 217 208 218 if (verboseLevel > 1) { << 209 G4int n = pv->GetVectorLength() - 1; 219 G4cout << "****** DEBUG: tcs value for Z=" << 210 G4double e = GammaEnergy/MeV; 220 G4cout << " cs (Geant4 internal unit)=" < << 211 if(e >= pv->Energy(n)) { 221 G4cout << " -> first E*E*cs value in CS << 212 xs = (*pv)[n]/(e*e); 222 G4cout << " -> last E*E*cs value in CS << 213 } else if(e >= pv->Energy(0)) { 223 G4cout << "******************************* << 214 xs = pv->Value(e)/(e*e); >> 215 } >> 216 >> 217 if(verboseLevel > 0) >> 218 { >> 219 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)=" << e << G4endl; >> 220 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl; >> 221 G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0] << G4endl; >> 222 G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n] << G4endl; >> 223 G4cout << "*********************************************************" << G4endl; 224 } 224 } 225 return xs; 225 return xs; 226 } 226 } 227 227 228 //....oooOO0OOooo........oooOO0OOooo........oo 228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 229 229 230 void G4LivermoreRayleighModel::SampleSecondari << 230 void G4LivermoreRayleighModel::SampleSecondaries( 231 << 231 std::vector<G4DynamicParticle*>*, 232 << 232 const G4MaterialCutsCouple* couple, 233 << 233 const G4DynamicParticle* aDynamicGamma, >> 234 G4double, G4double) 234 { 235 { 235 if (verboseLevel > 1) { 236 if (verboseLevel > 1) { 236 G4cout << "Calling SampleSecondaries() of << 237 G4cout << "Calling SampleSecondaries() of G4LivermoreRayleighModel" >> 238 << G4endl; 237 } 239 } 238 G4double photonEnergy0 = aDynamicGamma->GetK 240 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy(); 239 241 >> 242 // absorption of low-energy gamma >> 243 if (photonEnergy0 <= lowEnergyLimit) >> 244 { >> 245 fParticleChange->ProposeTrackStatus(fStopAndKill); >> 246 fParticleChange->SetProposedKineticEnergy(0.); >> 247 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0); >> 248 return ; >> 249 } >> 250 240 // Select randomly one element in the curren 251 // Select randomly one element in the current material 241 const G4ParticleDefinition* particle = aDyna << 252 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition(); 242 const G4Element* elm = SelectRandomAtom(coup << 253 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0); 243 G4int Z = elm->GetZasInt(); << 254 G4int Z = G4lrint(elm->GetZ()); 244 255 245 // Sample the angle of the scattered photon 256 // Sample the angle of the scattered photon 246 G4ThreeVector photonDirection = GetAngularDi << 257 247 aDynamicGamma, photonEnergy0, Z, couple->G << 258 G4ThreeVector photonDirection = >> 259 GetAngularDistribution()->SampleDirection(aDynamicGamma, >> 260 photonEnergy0, Z, couple->GetMaterial()); 248 fParticleChange->ProposeMomentumDirection(ph 261 fParticleChange->ProposeMomentumDirection(photonDirection); 249 } << 250 << 251 //....oooOO0OOooo........oooOO0OOooo........oo << 252 << 253 void G4LivermoreRayleighModel::InitialiseForEl << 254 { << 255 if (nullptr != dataCS[Z]) { << 256 return; << 257 } << 258 G4AutoLock l(&LivermoreRayleighModelMutex); << 259 if (nullptr == dataCS[Z]) { << 260 ReadData(Z); << 261 } << 262 l.unlock(); << 263 } 262 } 264 263 265 //....oooOO0OOooo........oooOO0OOooo........oo 264 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 266 265