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: Zhuxin Li@CENBG 26 // Author: Zhuxin Li@CENBG 27 // 11 March 2020 27 // 11 March 2020 28 // on the base of G4LivermoreGammaConv << 28 // on the base of G4LivermoreGammaConversionModel 29 // derives from G4BetheHeitler5DModel << 29 // derives from G4BetheHeitler5DModel 30 // ------------------------------------------- 30 // ------------------------------------------------------------------- 31 31 32 #include "G4LivermoreGammaConversion5DModel.hh 32 #include "G4LivermoreGammaConversion5DModel.hh" 33 << 34 #include "G4AutoLock.hh" << 35 #include "G4Electron.hh" 33 #include "G4Electron.hh" >> 34 #include "G4Positron.hh" 36 #include "G4EmParameters.hh" 35 #include "G4EmParameters.hh" 37 #include "G4Exp.hh" << 38 #include "G4ParticleChangeForGamma.hh" 36 #include "G4ParticleChangeForGamma.hh" 39 #include "G4PhysicalConstants.hh" << 40 #include "G4PhysicsFreeVector.hh" 37 #include "G4PhysicsFreeVector.hh" >> 38 #include "G4PhysicsLogVector.hh" >> 39 #include "G4ProductionCutsTable.hh" >> 40 #include "G4PhysicalConstants.hh" 41 #include "G4SystemOfUnits.hh" 41 #include "G4SystemOfUnits.hh" >> 42 #include "G4Exp.hh" >> 43 #include "G4AutoLock.hh" 42 44 43 namespace << 45 namespace { G4Mutex LivermoreGammaConversion5DModelMutex = G4MUTEX_INITIALIZER; } 44 { << 45 G4Mutex LivermoreGammaConversion5DModelMutex = << 46 } << 47 46 48 //....oooOO0OOooo........oooOO0OOooo........oo 47 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 49 //....oooOO0OOooo........oooOO0OOooo........oo 48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 50 49 >> 50 const G4int G4LivermoreGammaConversion5DModel::maxZ; >> 51 G4double G4LivermoreGammaConversion5DModel::lowEnergyLimit = 2.*CLHEP::electron_mass_c2; 51 G4PhysicsFreeVector* G4LivermoreGammaConversio 52 G4PhysicsFreeVector* G4LivermoreGammaConversion5DModel::data[] = {nullptr}; 52 G4String G4LivermoreGammaConversion5DModel::gD << 53 53 54 G4LivermoreGammaConversion5DModel::G4Livermore << 54 G4LivermoreGammaConversion5DModel::G4LivermoreGammaConversion5DModel(const G4ParticleDefinition* p, 55 << 55 const G4String& nam) 56 : G4BetheHeitler5DModel(p, nam) << 56 : G4BetheHeitler5DModel(p, nam), fParticleChange(nullptr) 57 { 57 { 58 fParticleChange = nullptr; << 59 lowEnergyLimit = 2. * CLHEP::electron_mass_c << 60 verboseLevel = 0; 58 verboseLevel = 0; 61 // Verbosity scale for debugging purposes: 59 // Verbosity scale for debugging purposes: 62 // 0 = nothing << 60 // 0 = nothing 63 // 1 = calculation of cross sections, file o 61 // 1 = calculation of cross sections, file openings... 64 // 2 = entering in methods 62 // 2 = entering in methods 65 if (verboseLevel > 0) { << 63 if(verboseLevel > 0) >> 64 { 66 G4cout << "G4LivermoreGammaConversion5DMod 65 G4cout << "G4LivermoreGammaConversion5DModel is constructed " << G4endl; 67 } 66 } 68 } 67 } 69 68 70 //....oooOO0OOooo........oooOO0OOooo........oo 69 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 71 70 72 G4LivermoreGammaConversion5DModel::~G4Livermor 71 G4LivermoreGammaConversion5DModel::~G4LivermoreGammaConversion5DModel() 73 { 72 { 74 if (IsMaster()) { << 73 if(IsMaster()) { 75 for (G4int i = 0; i < maxZ; ++i) { << 74 for(G4int i=0; i<maxZ; ++i) { 76 if (data[i]) { << 75 if(data[i]) { 77 delete data[i]; << 76 delete data[i]; 78 data[i] = nullptr; << 77 data[i] = nullptr; 79 } << 80 } << 81 } << 82 } << 83 << 84 //....oooOO0OOooo........oooOO0OOooo........oo << 85 << 86 void G4LivermoreGammaConversion5DModel::Initia << 87 << 88 { << 89 G4BetheHeitler5DModel::Initialise(particle, << 90 << 91 if (verboseLevel > 1) { << 92 G4cout << "Calling Initialise() of G4Liver << 93 << "Energy range: " << LowEnergyLim << 94 << " GeV isMater: " << IsMaster() < << 95 } << 96 << 97 if (IsMaster()) { << 98 // Initialise element selector << 99 InitialiseElementSelectors(particle, cuts) << 100 << 101 // Access to elements << 102 const G4ElementTable* elemTable = G4Elemen << 103 std::size_t numElems = (*elemTable).size() << 104 for (std::size_t ie = 0; ie < numElems; ++ << 105 const G4Element* elem = (*elemTable)[ie] << 106 const G4int Z = std::min(maxZ, elem->Get << 107 if (data[Z] == nullptr) { << 108 ReadData(Z); << 109 } 78 } 110 } 79 } 111 } 80 } 112 << 113 if (isInitialised) { << 114 return; << 115 } << 116 fParticleChange = GetParticleChangeForGamma( << 117 isInitialised = true; << 118 } 81 } 119 82 120 //....oooOO0OOooo........oooOO0OOooo........oo 83 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 121 84 122 const G4String& G4LivermoreGammaConversion5DMo << 85 void 123 { << 86 G4LivermoreGammaConversion5DModel::Initialise( const G4ParticleDefinition* particle, 124 // no check in this method - environment var << 87 const G4DataVector& cuts) 125 if (gDataDirectory.empty()) { << 88 { 126 auto param = G4EmParameters::Instance(); << 89 G4BetheHeitler5DModel::Initialise(particle, cuts); 127 std::ostringstream ost; << 90 128 if (param->LivermoreDataDir() == "livermor << 91 if (verboseLevel > 1) 129 ost << param->GetDirLEDATA() << "/liverm << 92 { 130 useSpline = true; << 93 G4cout << "Calling Initialise() of G4LivermoreGammaConversion5DModel." >> 94 << G4endl >> 95 << "Energy range: " >> 96 << LowEnergyLimit() / MeV << " MeV - " >> 97 << HighEnergyLimit() / GeV << " GeV isMater: " << IsMaster() >> 98 << G4endl; >> 99 } >> 100 >> 101 if(!fParticleChange) { >> 102 fParticleChange = GetParticleChangeForGamma(); >> 103 } >> 104 >> 105 if(IsMaster()) >> 106 { >> 107 // Initialise element selector >> 108 InitialiseElementSelectors(particle, cuts); >> 109 // Access to elements >> 110 char* path = std::getenv("G4LEDATA"); >> 111 G4ProductionCutsTable* theCoupleTable = >> 112 G4ProductionCutsTable::GetProductionCutsTable(); >> 113 G4int numOfCouples = theCoupleTable->GetTableSize(); >> 114 for(G4int i=0; i<numOfCouples; ++i) >> 115 { >> 116 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i); >> 117 SetCurrentCouple(couple); >> 118 const G4Material* mat = couple->GetMaterial(); >> 119 const G4ElementVector* theElementVector = mat->GetElementVector(); >> 120 G4int nelm = mat->GetNumberOfElements(); >> 121 for (G4int j=0; j<nelm; ++j) >> 122 { >> 123 G4int Z = std::max(1, std::min((*theElementVector)[j]->GetZasInt(), maxZ)); >> 124 if(!data[Z]) { ReadData(Z, path); } >> 125 } >> 126 } >> 127 } >> 128 } >> 129 >> 130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 131 >> 132 void G4LivermoreGammaConversion5DModel::ReadData(size_t Z, const char* path) >> 133 { >> 134 if (verboseLevel > 1) >> 135 { >> 136 G4cout << "Calling ReadData() of G4LivermoreGammaConversion5DModel" >> 137 << G4endl; >> 138 } >> 139 >> 140 if(data[Z]) { return; } >> 141 const char* datadir = path; >> 142 if(!datadir) >> 143 { >> 144 datadir = std::getenv("G4LEDATA"); >> 145 if(!datadir) >> 146 { >> 147 G4Exception("G4LivermoreGammaConversion5DModel::ReadData()", >> 148 "em0006",FatalException, >> 149 "Environment variable G4LEDATA not defined"); >> 150 return; 131 } 151 } 132 else { << 133 ost << param->GetDirLEDATA() << "/epics2 << 134 } << 135 gDataDirectory = ost.str(); << 136 } 152 } 137 return gDataDirectory; << 153 data[Z] = new G4PhysicsFreeVector(); 138 } << 139 << 140 //....oooOO0OOooo........oooOO0OOooo........oo << 141 << 142 void G4LivermoreGammaConversion5DModel::ReadDa << 143 { << 144 if (verboseLevel > 1) { << 145 G4cout << "Calling ReadData() of G4Livermo << 146 } << 147 << 148 if (data[Z]) { << 149 return; << 150 } << 151 << 152 std::ostringstream ost; 154 std::ostringstream ost; 153 ost << FindDirectoryPath() << "pp-cs-" << Z << 155 ost << datadir << "/epics2017/pair/pp-cs-" << Z <<".dat"; 154 << 155 data[Z] = new G4PhysicsFreeVector(useSpline) << 156 << 157 std::ifstream fin(ost.str().c_str()); 156 std::ifstream fin(ost.str().c_str()); 158 << 157 159 if (!fin.is_open()) { << 158 if( !fin.is_open()) >> 159 { 160 G4ExceptionDescription ed; 160 G4ExceptionDescription ed; 161 ed << "G4LivermoreGammaConversion5DModel d << 161 ed << "G4LivermoreGammaConversion5DModel data file <" << ost.str().c_str() 162 << G4endl; << 162 << "> is not opened!" << G4endl; 163 G4Exception("G4LivermoreGammaConversion5DM << 163 G4Exception("G4LivermoreGammaConversion5DModel::ReadData()", 164 "G4LEDATA version should be G4 << 164 "em0003",FatalException, >> 165 ed,"G4LEDATA version should be G4EMLOW6.27 or later."); 165 return; 166 return; 166 } << 167 } 167 else { << 168 else 168 if (verboseLevel > 1) { << 169 { 169 G4cout << "File " << ost.str() << " is o << 170 if(verboseLevel > 1) { G4cout << "File " << ost.str() 170 } << 171 << " is opened by G4LivermoreGammaConversion5DModel" << G4endl;} 171 data[Z]->Retrieve(fin, true); << 172 data[Z]->Retrieve(fin, true); 172 } << 173 } 173 // Activation of spline interpolation << 174 if (useSpline) data[Z]->FillSecondDerivative << 175 } 174 } 176 175 177 //....oooOO0OOooo........oooOO0OOooo........oo 176 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 178 177 179 G4double << 178 G4double 180 G4LivermoreGammaConversion5DModel::ComputeCros << 179 G4LivermoreGammaConversion5DModel::ComputeCrossSectionPerAtom( 181 << 180 const G4ParticleDefinition* particle, G4double GammaEnergy, G4double Z, 182 << 181 G4double, G4double, G4double) 183 { << 182 { 184 if (verboseLevel > 1) { << 183 if (verboseLevel > 1) 185 G4cout << "G4LivermoreGammaConversion5DMod << 184 { >> 185 G4cout << "G4LivermoreGammaConversion5DModel::ComputeCrossSectionPerAtom() Z= " >> 186 << Z << G4endl; 186 } 187 } 187 G4double xs = 0.0; 188 G4double xs = 0.0; 188 if (GammaEnergy < lowEnergyLimit) { << 189 if (GammaEnergy < lowEnergyLimit) { return xs; } 189 return xs; << 190 190 } << 191 << 192 G4int intZ = std::max(1, std::min(G4lrint(Z) 191 G4int intZ = std::max(1, std::min(G4lrint(Z), maxZ)); 193 G4PhysicsFreeVector* pv = data[intZ]; 192 G4PhysicsFreeVector* pv = data[intZ]; 194 // if element was not initialised 193 // if element was not initialised 195 // do initialisation safely for MT mode 194 // do initialisation safely for MT mode 196 if (!pv) { << 195 if(!pv) >> 196 { 197 InitialiseForElement(particle, intZ); 197 InitialiseForElement(particle, intZ); 198 pv = data[intZ]; 198 pv = data[intZ]; 199 if (!pv) { << 199 if(!pv) { return xs; } 200 return xs; << 201 } << 202 } 200 } 203 // x-section is taken from the table 201 // x-section is taken from the table 204 xs = pv->Value(GammaEnergy); << 202 xs = pv->Value(GammaEnergy); 205 if (verboseLevel > 0) { << 203 if(verboseLevel > 0) 206 G4cout << "*** Gamma conversion xs for Z=" << 204 { 207 << " cs=" << xs / millibarn << " m << 205 G4cout << "*** Gamma conversion xs for Z=" << Z << " at energy E(MeV)=" >> 206 << GammaEnergy/MeV << " cs=" << xs/millibarn << " mb" << G4endl; 208 } 207 } 209 return xs; 208 return xs; 210 } 209 } 211 210 212 //....oooOO0OOooo........oooOO0OOooo........oo 211 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 213 212 214 void G4LivermoreGammaConversion5DModel::Initia << 213 void G4LivermoreGammaConversion5DModel::InitialiseForElement( >> 214 const G4ParticleDefinition*, >> 215 G4int Z) 215 { 216 { 216 G4AutoLock l(&LivermoreGammaConversion5DMode 217 G4AutoLock l(&LivermoreGammaConversion5DModelMutex); 217 if (!data[Z]) { << 218 if(!data[Z]) { ReadData(Z); } 218 ReadData(Z); << 219 } << 220 l.unlock(); 219 l.unlock(); 221 } 220 } 222 221 223 //....oooOO0OOooo........oooOO0OOooo........oo 222 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 224 223