| 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 // Class Description 26 // Class Description
27 // Manager of NetronHP 27 // Manager of NetronHP
28 // << 28 //
29 // 121031 First implementation done by T. Koi 29 // 121031 First implementation done by T. Koi (SLAC/PPA)
30 // P. Arce, June-2014 Conversion neutron_hp to 30 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
31 // V. Ivanchenko, July-2023 Basic revision of <<
32 // 31 //
33 #include "G4ParticleHPManager.hh" 32 #include "G4ParticleHPManager.hh"
34 <<
35 #include "G4Exception.hh" <<
36 #include "G4HadronicException.hh" <<
37 #include "G4ParticleHPMessenger.hh" <<
38 #include "G4ParticleDefinition.hh" <<
39 #include "G4HadronicParameters.hh" <<
40 #include "G4ParticleHPThreadLocalManager.hh" 33 #include "G4ParticleHPThreadLocalManager.hh"
41 #include "G4SystemOfUnits.hh" << 34 #include "G4ParticleHPMessenger.hh"
42 << 35 #include "G4HadronicException.hh"
43 #include <zlib.h> <<
44 #include <fstream> <<
45 36
46 G4ParticleHPManager* G4ParticleHPManager::inst << 37 //G4ThreadLocal G4ParticleHPManager* G4ParticleHPManager::instance = NULL;
>> 38 G4ParticleHPManager* G4ParticleHPManager::instance = G4ParticleHPManager::GetInstance();
47 39
48 G4ParticleHPManager::G4ParticleHPManager() 40 G4ParticleHPManager::G4ParticleHPManager()
49 : theMinEnergyDBRC(0.1 * CLHEP::eV), << 41 :/*RWB(NULL),*/
50 theMaxEnergyDBRC(210. * CLHEP::eV), << 42 verboseLevel(1)
51 theMaxEnergyDoppler(30. * CLHEP::keV) << 43 ,USE_ONLY_PHOTONEVAPORATION(false)
52 { << 44 ,SKIP_MISSING_ISOTOPES(false)
53 messenger = new G4ParticleHPMessenger(this); << 45 ,NEGLECT_DOPPLER(false)
54 verboseLevel = G4HadronicParameters::Instanc << 46 ,DO_NOT_ADJUST_FINAL_STATE(false)
55 const char* ss = G4FindDataDir("NeutronHPNam << 47 ,PRODUCE_FISSION_FRAGMENTS(false)
56 if (nullptr != ss) { CHECK_HP_NAMES = true; << 48 ,theElasticCrossSections(NULL)
57 ss = G4FindDataDir("G4PHP_DO_NOT_CHECK_DIFF_ << 49 ,theCaptureCrossSections(NULL)
58 if (nullptr != ss) { PHP_CHECK = false; } << 50 //,theInelasticCrossSections(NULL)
59 ss = G4FindDataDir("G4PHP_MULTIPLICITY_METHO << 51 ,theFissionCrossSections(NULL)
60 if (nullptr != ss && "BetweenInts" == G4Stri << 52 ,theElasticFSs(NULL)
61 ss = G4FindDataDir("G4ParticleHPDebug"); << 53 //,theInelasticFSs(NULL)
62 if (nullptr != ss) { DEBUG = true; } << 54 ,theCaptureFSs(NULL)
63 << 55 ,theFissionFSs(NULL)
64 // identify and check data path once - it sh << 56 ,theTSCoherentCrossSections(NULL)
65 const char* nch = G4FindDataDir("G4NEUTRONHP << 57 ,theTSIncoherentCrossSections(NULL)
66 if (nullptr == nch) { << 58 ,theTSInelasticCrossSections(NULL)
67 G4Exception("G4ParticleHPManager::G4Partic << 59 ,theTSCoherentFinalStates(NULL)
68 FatalException, "G4NEUTRONXSDA << 60 ,theTSIncoherentFinalStates(NULL)
69 } else { << 61 ,theTSInelasticFinalStates(NULL)
70 fDataPath[0] = G4String(nch); << 62 {
71 } << 63 messenger = new G4ParticleHPMessenger( this );
72 // path may be defined by two environment va << 64 if ( getenv( "G4NEUTRONHP_DO_NOT_ADJUST_FINAL_STATE" ) || getenv("G4PHP_DO_NOT_ADJUST_FINAL_STATE") ) DO_NOT_ADJUST_FINAL_STATE = true;
73 // it is not mandatory to access PHP data - << 65 if ( getenv( "G4NEUTRONHP_USE_ONLY_PHOTONEVAPORATION" ) ) USE_ONLY_PHOTONEVAPORATION = true;
74 const char* ttp = G4FindDataDir("G4PARTICLEH << 66 if ( getenv( "G4NEUTRONHP_NEGLECT_DOPPLER" ) || getenv("G4PHP_NEGLECT_DOPPLER") ) NEGLECT_DOPPLER = true;
75 G4String tendl = (nullptr == ttp) ? G4String << 67 if ( getenv( "G4NEUTRONHP_SKIP_MISSING_ISOTOPES" ) ) SKIP_MISSING_ISOTOPES = true;
76 const char* ssp = G4FindDataDir("G4PROTONHPD << 68 if ( getenv( "G4NEUTRONHP_PRODUCE_FISSION_FRAGMENTS" ) ) PRODUCE_FISSION_FRAGMENTS = true;
77 fDataPath[1] = (nullptr == ssp) ? tendl + "/ <<
78 <<
79 ssp = G4FindDataDir("G4DEUTERONHPDATA"); <<
80 fDataPath[2] = (nullptr == ssp) ? tendl + "/ <<
81 <<
82 ssp = G4FindDataDir("G4TRITONHPDATA"); <<
83 fDataPath[3] = (nullptr == ssp) ? tendl + "/ <<
84 <<
85 ssp = G4FindDataDir("G4HE3HPDATA"); <<
86 fDataPath[4] = (nullptr == ssp) ? tendl + "/ <<
87 <<
88 ssp = G4FindDataDir("G4ALPHAHPDATA"); <<
89 fDataPath[5] = (nullptr == ssp) ? tendl + "/ <<
90 } 69 }
91 <<
92 G4ParticleHPManager::~G4ParticleHPManager() 70 G4ParticleHPManager::~G4ParticleHPManager()
93 { 71 {
94 delete messenger; << 72 delete messenger;
95 } 73 }
96 <<
97 G4ParticleHPManager* G4ParticleHPManager::GetI <<
98 { <<
99 static G4ParticleHPManager manager; <<
100 if (instance == nullptr) { <<
101 instance = &manager; <<
102 } <<
103 return instance; <<
104 } <<
105 <<
106 G4int G4ParticleHPManager::GetPHPIndex(const G <<
107 G4int pdg = part->GetPDGEncoding(); <<
108 G4int idx; <<
109 if (pdg == 2112) { idx = 0; } <<
110 else if (pdg == 2212) { idx = 1; } <<
111 else if (pdg == 1000010020) { idx = 2; } <<
112 else if (pdg == 1000010030) { idx = 3; } <<
113 else if (pdg == 1000020030) { idx = 4; } <<
114 else if (pdg == 1000020040) { idx = 5; } <<
115 else { <<
116 idx = 0; <<
117 G4ExceptionDescription ed; <<
118 ed << "Particle " << part->GetParticleName <<
119 << " cannot be handled by the ParticleH <<
120 G4Exception("G4ParticleHPManager::G4Partic <<
121 FatalException, ed, ""); <<
122 } <<
123 return idx; <<
124 } <<
125 <<
126 const G4String& <<
127 G4ParticleHPManager::GetParticleHPPath(const G <<
128 return fDataPath[GetPHPIndex(part)]; <<
129 } <<
130 <<
131 void G4ParticleHPManager::OpenReactionWhiteBoa 74 void G4ParticleHPManager::OpenReactionWhiteBoard()
132 { 75 {
133 G4ParticleHPThreadLocalManager::GetInstance( << 76 // if ( RWB != NULL ) {
>> 77 // G4cout << "Warning: G4ParticleHPReactionWhiteBoard is tried doubly opening" << G4endl;
>> 78 // RWB = new G4ParticleHPReactionWhiteBoard();
>> 79 // }
>> 80 //
>> 81 // RWB = new G4ParticleHPReactionWhiteBoard();
>> 82 G4ParticleHPThreadLocalManager::GetInstance()->OpenReactionWhiteBoard();
134 } 83 }
135 <<
136 G4ParticleHPReactionWhiteBoard* G4ParticleHPMa 84 G4ParticleHPReactionWhiteBoard* G4ParticleHPManager::GetReactionWhiteBoard()
137 { 85 {
138 return G4ParticleHPThreadLocalManager::GetIn << 86 // if ( RWB == NULL ) {
>> 87 // G4cout << "Warning: try to access G4ParticleHPReactionWhiteBoard before opening" << G4endl;
>> 88 // RWB = new G4ParticleHPReactionWhiteBoard();
>> 89 // }
>> 90 // return RWB;
>> 91 return G4ParticleHPThreadLocalManager::GetInstance()->GetReactionWhiteBoard();
139 } 92 }
140 <<
141 void G4ParticleHPManager::CloseReactionWhiteBo 93 void G4ParticleHPManager::CloseReactionWhiteBoard()
142 { 94 {
143 G4ParticleHPThreadLocalManager::GetInstance( << 95 G4ParticleHPThreadLocalManager::GetInstance()->CloseReactionWhiteBoard();
144 } 96 }
145 97
146 void G4ParticleHPManager::GetDataStream(const << 98 #include "zlib.h"
>> 99 #include <fstream>
>> 100 void G4ParticleHPManager::GetDataStream( G4String filename , std::istringstream& iss )
147 { 101 {
148 G4String* data = nullptr; << 102 G4String* data=NULL;
149 G4String compfilename(filename); << 103 G4String compfilename(filename);
150 compfilename += ".z"; << 104 compfilename += ".z";
151 auto in = new std::ifstream(compfilename, st << 105 std::ifstream* in = new std::ifstream ( compfilename , std::ios::binary | std::ios::ate );
152 if (in->good()) { << 106 if ( in->good() )
153 // Use the compressed file << 107 {
154 std::streamoff file_size = in->tellg(); << 108 // Use the compressed file
155 in->seekg(0, std::ios::beg); << 109 G4int file_size = in->tellg();
156 auto compdata = new Bytef[file_size]; << 110 in->seekg( 0 , std::ios::beg );
157 << 111 Bytef* compdata = new Bytef[ file_size ];
158 while (*in) { // Loop checking, 11.05.201 << 112
159 in->read((char*)compdata, file_size); << 113 while ( *in ) { // Loop checking, 11.05.2015, T. Koi
160 } << 114 in->read( (char*)compdata , file_size );
161 <<
162 auto complen = (uLongf)(file_size * 4); <<
163 auto uncompdata = new Bytef[complen]; <<
164 <<
165 while (Z_OK != uncompress(uncompdata, &com <<
166 { // Loop checking, 11.05.2015, T. Koi <<
167 delete[] uncompdata; <<
168 complen *= 2; <<
169 uncompdata = new Bytef[complen]; <<
170 } <<
171 delete[] compdata; <<
172 // Now "complen" has uncomplessed size <<
173 data = new G4String((char*)uncompdata, (G4 <<
174 delete[] uncompdata; <<
175 } <<
176 else { <<
177 // Use regular text file <<
178 std::ifstream thefData(filename, std::ios: <<
179 if (thefData.good()) { <<
180 std::streamoff file_size = thefData.tell <<
181 thefData.seekg(0, std::ios::beg); <<
182 auto filedata = new char[file_size]; <<
183 while (thefData) { // Loop checking, 11 <<
184 thefData.read(filedata, file_size); <<
185 } 115 }
186 thefData.close(); <<
187 data = new G4String(filedata, file_size) <<
188 delete[] filedata; <<
189 } <<
190 else { <<
191 // found no data file <<
192 // set error bit to the stream <<
193 iss.setstate(std::ios::badbit); <<
194 } <<
195 } <<
196 if (data != nullptr) { <<
197 iss.str(*data); <<
198 G4String id; <<
199 iss >> id; <<
200 if (id == "G4NDL") { <<
201 // Register information of file <<
202 G4String source; <<
203 iss >> source; <<
204 register_data_file(filename, source); <<
205 } <<
206 else { <<
207 iss.seekg(0, std::ios::beg); <<
208 } <<
209 } <<
210 in->close(); <<
211 delete in; <<
212 delete data; <<
213 } <<
214 116
215 void G4ParticleHPManager::GetDataStream2(const << 117 uLongf complen = (uLongf) ( file_size*4 );
216 { << 118 Bytef* uncompdata = new Bytef[complen];
217 // Checking existance of data file <<
218 119
219 G4String compfilename(filename); << 120 while ( Z_OK != uncompress ( uncompdata , &complen , compdata , file_size ) ) { // Loop checking, 11.05.2015, T. Koi
220 compfilename += ".z"; << 121 //G4cout << "Too small, retry 2 times bigger size." << G4endl;
221 auto in = new std::ifstream(compfilename, st << 122 delete[] uncompdata;
222 if (in->good()) { << 123 complen *= 2;
223 // Compressed file is exist << 124 uncompdata = new Bytef[complen];
224 in->close(); << 125 }
225 } << 126 delete [] compdata;
226 else { << 127 // Now "complen" has uncomplessed size
227 std::ifstream thefData(filename, std::ios: << 128 data = new G4String ( (char*)uncompdata , (G4long)complen );
228 if (thefData.good()) { << 129 delete [] uncompdata;
229 // Regular text file is exist << 130 } else {
230 thefData.close(); << 131 // Use regular text file
231 } << 132 std::ifstream thefData( filename , std::ios::in | std::ios::ate );
232 else { << 133 if ( thefData.good() ) {
233 // found no data file << 134 G4int file_size = thefData.tellg();
234 // set error bit to the stream << 135 thefData.seekg( 0 , std::ios::beg );
235 iss.setstate(std::ios::badbit); << 136 char* filedata = new char[ file_size ];
236 } << 137 while ( thefData ) { // Loop checking, 11.05.2015, T. Koi
237 } << 138 thefData.read( filedata , file_size );
238 delete in; << 139 }
>> 140 thefData.close();
>> 141 data = new G4String ( filedata , file_size );
>> 142 delete [] filedata;
>> 143 } else {
>> 144 // found no data file
>> 145 // set error bit to the stream
>> 146 iss.setstate( std::ios::badbit );
>> 147 }
>> 148 }
>> 149 if ( data != NULL ) {
>> 150 iss.str(*data);
>> 151 G4String id;
>> 152 iss >> id;
>> 153 if ( id == "G4NDL" ) {
>> 154 //Register information of file
>> 155 G4String source;
>> 156 iss >> source;
>> 157 register_data_file(filename,source);
>> 158 } else {
>> 159 iss.seekg( 0 , std::ios::beg );
>> 160 }
>> 161 }
>> 162 //G4cout << iss.rdbuf()->in_avail() << G4endl;
>> 163 in->close(); delete in;
>> 164 delete data;
>> 165 }
>> 166 // Checking existance of data file
>> 167 void G4ParticleHPManager::GetDataStream2( G4String filename , std::istringstream& iss )
>> 168 {
>> 169 G4String compfilename(filename);
>> 170 compfilename += ".z";
>> 171 std::ifstream* in = new std::ifstream ( compfilename , std::ios::binary | std::ios::ate );
>> 172 if ( in->good() )
>> 173 {
>> 174 // Compressed file is exist
>> 175 in->close();
>> 176 } else {
>> 177 std::ifstream thefData( filename , std::ios::in | std::ios::ate );
>> 178 if ( thefData.good() ) {
>> 179 // Regular text file is exist
>> 180 thefData.close();
>> 181 } else {
>> 182 // found no data file
>> 183 // set error bit to the stream
>> 184 iss.setstate( std::ios::badbit );
>> 185 }
>> 186 }
>> 187 delete in;
239 } 188 }
240 189
241 void G4ParticleHPManager::SetVerboseLevel(G4in << 190 void G4ParticleHPManager::SetVerboseLevel( G4int newValue )
242 { 191 {
243 G4cout << "You are setting a new verbose lev << 192 G4cout << "You are setting a new verbose level for Particle HP package." << G4endl;
244 G4cout << "the new value will be used in who << 193 G4cout << "the new value will be used in whole of the Particle HP package, i.e., models and cross sections for Capture, Elastic, Fission and Inelastic interaction." << G4endl;
245 "cross sections for Capture, Elast << 194 verboseLevel = newValue;
246 << G4endl; <<
247 verboseLevel = newValue; <<
248 } 195 }
249 196
250 void G4ParticleHPManager::register_data_file(c << 197 void G4ParticleHPManager::register_data_file(G4String filename, G4String source)
251 { 198 {
252 mDataEvaluation.insert(std::pair<G4String, G << 199 mDataEvaluation.insert( std::pair < G4String , G4String > ( filename , source ) );
253 } 200 }
254 201
255 void G4ParticleHPManager::DumpDataSource() con << 202 void G4ParticleHPManager::DumpDataSource()
256 { 203 {
257 G4cout << "Data source of this Partile HP ca << 204
258 for (const auto& it : mDataEvaluation) { << 205 G4cout << "Data source of this Partile HP calculation are " << G4endl;
259 G4cout << it.first << " " << it.second << << 206 for ( std::map< G4String , G4String >::iterator
260 } << 207 it = mDataEvaluation.begin() ; it != mDataEvaluation.end() ; it++ ) {
261 G4cout << G4endl; << 208 G4cout << it->first << " " << it->second << G4endl;
>> 209 }
>> 210 G4cout << G4endl;
262 } 211 }
263 212
264 void G4ParticleHPManager::DumpSetting() << 213 G4PhysicsTable* G4ParticleHPManager::GetInelasticCrossSections(const G4ParticleDefinition* particle ){
265 { << 214 if ( theInelasticCrossSections.end() != theInelasticCrossSections.find( particle ) )
266 if(isPrinted) { return; } << 215 return theInelasticCrossSections.find( particle )->second;
267 G4cout << G4endl << 216 else
268 << "================================= << 217 return NULL;
269 << "====== ParticleHP Physics P << 218 }
270 << "================================= << 219
271 << " Use only photo-evaporation << 220 void G4ParticleHPManager::RegisterInelasticCrossSections( const G4ParticleDefinition* particle, G4PhysicsTable* val ){
272 << " Skip missing isotopes << 221 theInelasticCrossSections.insert( std::pair<const G4ParticleDefinition* , G4PhysicsTable* >( particle , val ) );
273 << " Neglect Doppler << 222 }
274 << " Do not adjust final state << 223
275 << " Produce fission fragments << 224 std::vector<G4ParticleHPChannelList*>* G4ParticleHPManager::GetInelasticFinalStates(const G4ParticleDefinition* particle) {
276 << " Use WendtFissionModel << 225 if ( theInelasticFSs.end() != theInelasticFSs.find( particle ) )
277 << " Use NRESP71Model << 226 return theInelasticFSs.find( particle )->second;
278 << " Use DBRC << 227 else
279 << " PHP use Poisson << 228 return NULL;
280 << " PHP check << 229 }
281 << " CHECK HP NAMES << 230
282 << " Enable DEBUG << 231 void G4ParticleHPManager::RegisterInelasticFinalStates( const G4ParticleDefinition* particle , std::vector<G4ParticleHPChannelList*>* val ) {
283 << " Use probability tables from << 232 theInelasticFSs.insert ( std::pair<const G4ParticleDefinition*,std::vector<G4ParticleHPChannelList*>*>( particle , val ) );
284 << "================================= <<
285 isPrinted = true; <<
286 } 233 }
287 234