| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer << 3 // * DISCLAIMER *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th << 5 // * The following disclaimer summarizes all the specific disclaimers *
6 // * the Geant4 Collaboration. It is provided << 6 // * of contributors to this software. The specific disclaimers,which *
7 // * conditions of the Geant4 Software License << 7 // * govern, are listed with their locations in: *
8 // * LICENSE and available at http://cern.ch/ << 8 // * http://cern.ch/geant4/license *
9 // * include a list of copyright holders. <<
10 // * 9 // * *
11 // * Neither the authors of this software syst 10 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 11 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 12 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 13 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file << 14 // * use. *
16 // * for the full disclaimer and the limitatio <<
17 // * 15 // * *
18 // * This code implementation is the result << 16 // * This code implementation is the intellectual property of the *
19 // * technical work of the GEANT4 collaboratio << 17 // * GEANT4 collaboration. *
20 // * By using, copying, modifying or distri << 18 // * By copying, distributing or modifying the Program (or any work *
21 // * any work based on the software) you ag << 19 // * based on the Program) you indicate your acceptance of this *
22 // * use in resulting scientific publicati << 20 // * statement, and all its terms. *
23 // * acceptance of all terms of the Geant4 Sof <<
24 // ******************************************* 21 // ********************************************************************
25 // 22 //
26 // G4VUserPhysicsList implementation <<
27 // 23 //
28 // Original author: H.Kurashige (Kobe Universi << 24 // $Id: G4VUserPhysicsList.cc,v 1.13.2.1 2001/06/28 19:15:21 gunter Exp $
29 // ------------------------------------------- << 25 // GEANT4 tag $Name: $
>> 26 //
>> 27 //
>> 28 // ------------------------------------------------------------
>> 29 // GEANT 4 class header file
>> 30 //
>> 31 // ------------------------------------------------------------
>> 32 // History
>> 33 // first version 09 Jan. 1998 by H.Kurashige
>> 34 // modified 24 Jan. 1998 by H.Kurashige
>> 35 // modified 06 June 1998 by H.Kurashige
>> 36 // add G4ParticleWithCuts::SetEnergyRange
>> 37 // 18 June 1998 by H.Kurashige
>> 38 // modifeid for short lived particles 27 June 1998 by H.Kurashige
>> 39 // G4BestUnit on output 12 nov. 1998 mma
>> 40 // Added RemoveProcessManager 9 Feb. 1999 by H.Kurashige
>> 41 // Fixed RemoveProcessManager 15 Apr. 1999 by H.Kurashige
>> 42 // Removed ConstructAllParticles() 15 Apr. 1999 by H.Kurashige
>> 43 // modified 08, Mar 2001 by H.Kurashige
>> 44 // ------------------------------------------------------------
30 45
>> 46 #include "globals.hh"
31 #include "G4VUserPhysicsList.hh" 47 #include "G4VUserPhysicsList.hh"
32 << 48 #include "G4ParticleWithCuts.hh"
33 #include "G4Material.hh" <<
34 #include "G4MaterialCutsCouple.hh" <<
35 #include "G4ParticleDefinition.hh" <<
36 #include "G4ParticleTable.hh" <<
37 #include "G4PhysicsListHelper.hh" <<
38 #include "G4ProcessManager.hh" 49 #include "G4ProcessManager.hh"
39 #include "G4ProductionCuts.hh" << 50 #include "G4ParticleTable.hh"
40 #include "G4ProductionCutsTable.hh" << 51 #include "G4ParticleWithCuts.hh"
41 #include "G4Region.hh" << 52 #include "G4Material.hh"
42 #include "G4RegionStore.hh" << 53 #include "G4UserPhysicsListMessenger.hh"
43 #include "G4SystemOfUnits.hh" <<
44 #include "G4UImanager.hh" 54 #include "G4UImanager.hh"
45 #include "G4UnitsTable.hh" 55 #include "G4UnitsTable.hh"
46 #include "G4UserPhysicsListMessenger.hh" <<
47 #include "G4VTrackingManager.hh" <<
48 #include "G4ios.hh" <<
49 #include "globals.hh" <<
50 <<
51 #include <fstream> <<
52 #include <iomanip> <<
53 #include <unordered_set> <<
54 56
55 // This static member is thread local. For eac << 57 #include "G4ios.hh"
56 // size of G4VUPLData instances. << 58 #include "g4std/iomanip"
57 // << 59 #include "g4std/fstream"
58 #define G4MT_theMessenger ((this->subInstanceM <<
59 #define G4MT_thePLHelper ((this->subInstanceMa <<
60 #define fIsPhysicsTableBuilt \ <<
61 ((this->subInstanceManager.offset[this->g4vu <<
62 #define fDisplayThreshold \ <<
63 ((this->subInstanceManager.offset[this->g4vu <<
64 #define theParticleIterator \ <<
65 ((this->subInstanceManager.offset[this->g4vu <<
66 60
67 // This field helps to use the class G4VUPLMan << 61 ////////////////////////////////////////////////////////
68 // << 62 G4VUserPhysicsList::G4VUserPhysicsList()
69 G4VUPLManager G4VUserPhysicsList::subInstanceM << 63 :verboseLevel(1),
>> 64 numberOfMaterial(0),
>> 65 fRetrievePhysicsTable(false),
>> 66 directoryPhysicsTable("."),
>> 67 fIsCheckedForRetrievePhysicsTable(false),
>> 68 fStoredInAscii(true),
>> 69 fIsRestoredCutValues(false)
70 70
71 // ------------------------------------------- <<
72 void G4VUPLData::initialize() <<
73 { 71 {
74 _theParticleIterator = G4ParticleTable::GetP << 72 // default cut value (1.0mm)
75 _theMessenger = nullptr; << 73 defaultCutValue = 1.0*mm;
76 _thePLHelper = G4PhysicsListHelper::GetPhysi <<
77 _fIsPhysicsTableBuilt = false; <<
78 _fDisplayThreshold = 0; <<
79 } <<
80 74
81 // ------------------------------------------- << 75 // set energy range for SetCut calcuration
82 G4VUserPhysicsList::G4VUserPhysicsList() << 76 G4ParticleWithCuts::SetEnergyRange(0.99*keV, 100*TeV);
83 { <<
84 g4vuplInstanceID = subInstanceManager.Create <<
85 // default cut value (1.0mm) <<
86 defaultCutValue = 1.0 * mm; <<
87 77
88 // pointer to the particle table 78 // pointer to the particle table
89 theParticleTable = G4ParticleTable::GetParti 79 theParticleTable = G4ParticleTable::GetParticleTable();
90 << 80 theParticleIterator = theParticleTable->GetIterator();
91 // pointer to the cuts table << 81
92 fCutsTable = G4ProductionCutsTable::GetProdu << 82 // UI Messenger
93 << 83 theMessenger = new G4UserPhysicsListMessenger(this);
94 // set energy range for SetCut calcuration <<
95 fCutsTable->SetEnergyRange(0.99 * keV, 100 * <<
96 <<
97 // UI Messenger <<
98 G4MT_theMessenger = new G4UserPhysicsListMes <<
99 <<
100 // PhysicsListHelper <<
101 G4MT_thePLHelper->SetVerboseLevel(verboseLev <<
102 <<
103 fIsPhysicsTableBuilt = false; <<
104 fDisplayThreshold = 0; <<
105 } <<
106 <<
107 // ------------------------------------------- <<
108 void G4VUserPhysicsList::InitializeWorker() <<
109 { <<
110 // Remember messengers are per-thread, so th <<
111 // worker and due to the presence of "this" <<
112 // G4VUPLData::initialize() <<
113 G4MT_theMessenger = new G4UserPhysicsListMes <<
114 } 84 }
115 85
116 // ------------------------------------------- << 86 ////////////////////////////////////////////////////////
117 void G4VUserPhysicsList::TerminateWorker() << 87 G4VUserPhysicsList::~G4VUserPhysicsList()
118 { 88 {
119 RemoveProcessManager(); << 89 if (theMessenger != 0) {
120 RemoveTrackingManager(); << 90 delete theMessenger;
121 delete G4MT_theMessenger; << 91 theMessenger = 0;
122 G4MT_theMessenger = nullptr; << 92 }
123 } 93 }
124 94
125 // ------------------------------------------- << 95 ////////////////////////////////////////////////////////
126 G4VUserPhysicsList::~G4VUserPhysicsList() << 96 void G4VUserPhysicsList::AddProcessManager(G4ParticleDefinition* newParticle,
127 { << 97 G4ProcessManager* newManager)
128 delete G4MT_theMessenger; << 98 {
129 G4MT_theMessenger = nullptr; << 99 if (newParticle == 0) return;
>> 100 if (newParticle->GetProcessManager() != 0) {
>> 101 #ifdef G4VERBOSE
>> 102 if (verboseLevel >1){
>> 103 G4cout << "G4VUserPhysicsList::AddProcessManager: ";
>> 104 G4cout << newParticle->GetParticleName();
>> 105 G4cout << " already has ProcessManager " << G4endl;
>> 106 }
>> 107 #endif
>> 108 return;
>> 109 }
130 110
131 RemoveProcessManager(); << 111 // create new process manager if newManager == 0
132 RemoveTrackingManager(); << 112 if (newManager == 0){
>> 113 // Add ProcessManager
>> 114 if (newParticle->GetParticleType() == "nucleus") {
>> 115 // Create a copy of the process manager of "GenericIon" in case of "nucleus"
>> 116 G4ParticleDefinition* genericIon =
>> 117 (G4ParticleTable::GetParticleTable())->FindParticle("GenericIon");
133 118
134 // invoke DeleteAllParticle << 119 if (genericIon != 0) {
135 theParticleTable->DeleteAllParticles(); << 120 G4ProcessManager* ionMan = genericIon->GetProcessManager();
136 } << 121 if (ionMan != 0) {
>> 122 newManager = new G4ProcessManager(*ionMan);
>> 123 } else {
>> 124 // no process manager has been registered yet
>> 125 newManager = new G4ProcessManager(newParticle);
>> 126 }
>> 127 } else {
>> 128 // "GenericIon" does not exist
>> 129 newManager = new G4ProcessManager(newParticle);
>> 130 }
137 131
138 // ------------------------------------------- << 132 } else {
139 G4VUserPhysicsList::G4VUserPhysicsList(const G << 133 // create process manager for particles other than "nucleus"
140 : verboseLevel(right.verboseLevel), << 134 newManager = new G4ProcessManager(newParticle);
141 defaultCutValue(right.defaultCutValue), << 135 }
142 isSetDefaultCutValue(right.isSetDefaultCut << 136 }
143 fRetrievePhysicsTable(right.fRetrievePhysi <<
144 fStoredInAscii(right.fStoredInAscii), <<
145 fIsCheckedForRetrievePhysicsTable(right.fI <<
146 fIsRestoredCutValues(right.fIsRestoredCutV <<
147 directoryPhysicsTable(right.directoryPhysi <<
148 fDisableCheckParticleList(right.fDisableCh <<
149 { <<
150 g4vuplInstanceID = subInstanceManager.Create <<
151 // pointer to the particle table <<
152 theParticleTable = G4ParticleTable::GetParti <<
153 theParticleIterator = theParticleTable->GetI <<
154 // pointer to the cuts table <<
155 fCutsTable = G4ProductionCutsTable::GetProdu <<
156 137
157 // UI Messenger << 138 // set particle type
158 G4MT_theMessenger = new G4UserPhysicsListMes << 139 newManager->SetParticleType(newParticle);
159 140
160 // PhysicsListHelper << 141 // add the process manager
161 G4MT_thePLHelper = G4PhysicsListHelper::GetP << 142 newParticle->SetProcessManager(newManager);
162 G4MT_thePLHelper->SetVerboseLevel(verboseLev <<
163 143
164 fIsPhysicsTableBuilt = << 144 #ifdef G4VERBOSE
165 right.GetSubInstanceManager().offset[right << 145 if (verboseLevel >2){
166 fDisplayThreshold = << 146 G4cout << "G4VUserPhysicsList::AddProcessManager: ";
167 right.GetSubInstanceManager().offset[right << 147 G4cout << "adds ProcessManager to ";
168 } << 148 G4cout << newParticle->GetParticleName() << G4endl;
>> 149 newManager->DumpInfo();
>> 150 }
>> 151 #endif
169 152
170 // ------------------------------------------- <<
171 G4VUserPhysicsList& G4VUserPhysicsList::operat <<
172 { <<
173 if (this != &right) { <<
174 verboseLevel = right.verboseLevel; <<
175 defaultCutValue = right.defaultCutValue; <<
176 isSetDefaultCutValue = right.isSetDefaultC <<
177 fRetrievePhysicsTable = right.fRetrievePhy <<
178 fStoredInAscii = right.fStoredInAscii; <<
179 fIsCheckedForRetrievePhysicsTable = right. <<
180 fIsRestoredCutValues = right.fIsRestoredCu <<
181 directoryPhysicsTable = right.directoryPhy <<
182 fIsPhysicsTableBuilt = <<
183 right.GetSubInstanceManager().offset[rig <<
184 fDisplayThreshold = <<
185 right.GetSubInstanceManager().offset[rig <<
186 fDisableCheckParticleList = right.fDisable <<
187 } <<
188 return *this; <<
189 } 153 }
190 154
191 // ------------------------------------------- <<
192 void G4VUserPhysicsList::AddProcessManager(G4P <<
193 { <<
194 if (newParticle == nullptr) return; <<
195 G4Exception("G4VUserPhysicsList::AddProcessM <<
196 "This method is obsolete"); <<
197 } <<
198 155
199 // ------------------------------------------- << 156 ////////////////////////////////////////////////////////
200 void G4VUserPhysicsList::InitializeProcessMana 157 void G4VUserPhysicsList::InitializeProcessManager()
201 { 158 {
202 // Request lock for particle table access. S << 159 // loop over all particles in G4ParticleTable
203 // this critical region. <<
204 #ifdef G4MULTITHREADED <<
205 G4MUTEXLOCK(&G4ParticleTable::particleTableM <<
206 G4ParticleTable::lockCount()++; <<
207 #endif <<
208 G4ParticleDefinition* gion = G4ParticleTable <<
209 <<
210 // loop over all particles in G4ParticleTabl <<
211 theParticleIterator->reset(); 160 theParticleIterator->reset();
212 while ((*theParticleIterator)()) { << 161 while( (*theParticleIterator)() ){
213 G4ParticleDefinition* particle = thePartic 162 G4ParticleDefinition* particle = theParticleIterator->value();
214 G4ProcessManager* pmanager = particle->Get 163 G4ProcessManager* pmanager = particle->GetProcessManager();
215 << 164 if (pmanager==0) {
216 if (pmanager == nullptr) { << 165 // create process manager if the particle has no its one
217 // create process manager if the particl <<
218 pmanager = new G4ProcessManager(particle 166 pmanager = new G4ProcessManager(particle);
219 particle->SetProcessManager(pmanager); 167 particle->SetProcessManager(pmanager);
220 if (particle->GetMasterProcessManager() <<
221 particle->SetMasterProcessManager(pman <<
222 #ifdef G4VERBOSE <<
223 if (verboseLevel > 2) { <<
224 G4cout << "G4VUserPhysicsList::Initial <<
225 "ProcessManager to " <<
226 << particle->GetParticleName() <<
227 } <<
228 #endif <<
229 } <<
230 } <<
231 <<
232 if (gion != nullptr) { <<
233 G4ProcessManager* gionPM = gion->GetProces <<
234 // loop over all particles once again (thi <<
235 theParticleIterator->reset(false); <<
236 while ((*theParticleIterator)()) { <<
237 G4ParticleDefinition* particle = thePart <<
238 if (particle->IsGeneralIon()) { <<
239 particle->SetProcessManager(gionPM); <<
240 #ifdef G4VERBOSE <<
241 if (verboseLevel > 2) { <<
242 G4cout << "G4VUserPhysicsList::Initi <<
243 "ProcessManager to " <<
244 << particle->GetParticleName( <<
245 } <<
246 #endif <<
247 } <<
248 } 168 }
249 } 169 }
250 <<
251 // release lock for particle table access <<
252 #ifdef G4MULTITHREADED <<
253 G4MUTEXUNLOCK(&G4ParticleTable::particleTabl <<
254 #endif <<
255 } 170 }
256 171
257 // ------------------------------------------- << 172 ////////////////////////////////////////////////////////
258 void G4VUserPhysicsList::RemoveProcessManager( 173 void G4VUserPhysicsList::RemoveProcessManager()
259 { 174 {
260 // Request lock for particle table access. S << 175 // loop over all particles in G4ParticleTable
261 // this critical region. <<
262 #ifdef G4MULTITHREADED <<
263 G4MUTEXLOCK(&G4ParticleTable::particleTableM <<
264 G4ParticleTable::lockCount()++; <<
265 #endif <<
266 <<
267 // loop over all particles in G4ParticleTabl <<
268 theParticleIterator->reset(); 176 theParticleIterator->reset();
269 while ((*theParticleIterator)()) { << 177 while( (*theParticleIterator)() ){
270 G4ParticleDefinition* particle = thePartic 178 G4ParticleDefinition* particle = theParticleIterator->value();
271 if (particle->GetInstanceID() < G4Particle << 179 G4ProcessManager* pmanager = particle->GetProcessManager();
272 if (particle->GetParticleSubType() != "g << 180 if (pmanager!=0) delete pmanager;
273 || particle->GetParticleName() == "G << 181 particle->SetProcessManager(0);
274 { << 182 #ifdef G4VERBOSE
275 G4ProcessManager* pmanager = particle- << 183 if (verboseLevel >2){
276 << 184 G4cout << "G4VUserPhysicsList::RemoveProcessManager: ";
277 delete pmanager; << 185 G4cout << "remove ProcessManager from ";
278 #ifdef G4VERBOSE << 186 G4cout << particle->GetParticleName() << G4endl;
279 if (verboseLevel > 2) { <<
280 G4cout << "G4VUserPhysicsList::Remov <<
281 G4cout << "remove ProcessManager fro <<
282 G4cout << particle->GetParticleName( <<
283 } <<
284 #endif <<
285 } <<
286 particle->SetProcessManager(nullptr); <<
287 } 187 }
288 } <<
289 <<
290 // release lock for particle table access <<
291 #ifdef G4MULTITHREADED <<
292 G4MUTEXUNLOCK(&G4ParticleTable::particleTabl <<
293 #endif 188 #endif
>> 189 }
294 } 190 }
295 191
296 // ------------------------------------------- << 192
297 void G4VUserPhysicsList::RemoveTrackingManager << 193 ////////////////////////////////////////////////////////
>> 194 #include "G4Transportation.hh"
>> 195 void G4VUserPhysicsList::AddTransportation()
298 { 196 {
299 // One tracking manager may be registered fo << 197 G4Transportation* theTransportationProcess= new G4Transportation();
300 // to delete every object only once. <<
301 std::unordered_set<G4VTrackingManager*> trac <<
302 198
303 // loop over all particles in G4ParticleTabl << 199 // loop over all particles in G4ParticleTable
304 theParticleIterator->reset(); 200 theParticleIterator->reset();
305 while ((*theParticleIterator)()) { << 201 while( (*theParticleIterator)() ){
306 G4ParticleDefinition* particle = thePartic 202 G4ParticleDefinition* particle = theParticleIterator->value();
307 if (auto* trackingManager = particle->GetT << 203 G4ProcessManager* pmanager = particle->GetProcessManager();
308 #ifdef G4VERBOSE << 204 if (!particle->IsShortLived()) {
309 if (verboseLevel > 2) { << 205 // Add transportation process for all particles other than "shortlived"
310 G4cout << "G4VUserPhysicsList::RemoveT << 206 if ( pmanager == 0) {
311 G4cout << "remove TrackingManager from << 207 // Error !! no process manager
312 G4cout << particle->GetParticleName() << 208 G4Exception("G4VUserPhysicsList::AddTransportation : no process manager!");
>> 209 } else {
>> 210 // add transportation with ordering = ( -1, "first", "first" )
>> 211 pmanager ->AddProcess(theTransportationProcess);
>> 212 pmanager ->SetProcessOrderingToFirst(theTransportationProcess, idxAlongStep);
>> 213 pmanager ->SetProcessOrderingToFirst(theTransportationProcess, idxPostStep);
313 } 214 }
314 #endif << 215 } else {
315 trackingManagers.insert(trackingManager) << 216 // shortlived particle case
316 particle->SetTrackingManager(nullptr); <<
317 } 217 }
318 } 218 }
319 <<
320 for (G4VTrackingManager* tm : trackingManage <<
321 delete tm; <<
322 } <<
323 } 219 }
324 220
325 // ------------------------------------------- <<
326 void G4VUserPhysicsList::SetCuts() <<
327 { <<
328 if (!isSetDefaultCutValue) { <<
329 SetDefaultCutValue(defaultCutValue); <<
330 } <<
331 221
332 #ifdef G4VERBOSE << 222 ////////////////////////////////////////////////////////
333 if (verboseLevel > 1) { << 223 void G4VUserPhysicsList::SetDefaultCutValue(G4double value)
334 G4cout << "G4VUserPhysicsList::SetCuts: << 224 {
335 G4cout << "Cut for gamma: " << GetCutValue << 225 if (value<=0.0) {
336 G4cout << "Cut for e-: " << GetCutValue(" << 226 #ifdef G4VERBOSE
337 G4cout << "Cut for e+: " << GetCutValue(" << 227 if (verboseLevel >0){
338 G4cout << "Cut for proton: " << GetCutVal << 228 G4cout << "G4VUserPhysicsList::SetDefaultCutValue: negative cut values";
>> 229 G4cout << " :" << value/mm << "[mm]" << G4endl;
>> 230 }
>> 231 #endif
>> 232 } else {
>> 233 #ifdef G4VERBOSE
>> 234 if (verboseLevel >1){
>> 235 G4cout << "G4VUserPhysicsList::SetDefaultCutValue:";
>> 236 G4cout << "default cut value is changed to :" ;
>> 237 G4cout << value/mm << "[mm]" << G4endl;
>> 238 }
>> 239 #endif
>> 240 defaultCutValue = value;
>> 241 ResetCuts();
>> 242 }
>> 243 }
>> 244
>> 245 ////////////////////////////////////////////////////////
>> 246 void G4VUserPhysicsList::ResetCuts()
>> 247 {
>> 248 #ifdef G4VERBOSE
>> 249 if (verboseLevel >1) {
>> 250 G4cout << "G4VUserPhysicsList::ResetCuts()" << G4endl;
>> 251 G4cout << " cut values in energy will be calculated later" << G4endl;
339 } 252 }
340 #endif 253 #endif
341 254
342 // dump Cut values if verboseLevel==3 << 255 // Reset cut values for other particles
343 if (verboseLevel > 2) { << 256 // loop over all particles in G4ParticleTable
344 DumpCutValuesTable(); << 257 theParticleIterator->reset();
>> 258 while( (*theParticleIterator)() ){
>> 259 G4ParticleDefinition* particle = theParticleIterator->value();
>> 260 if (!particle->IsShortLived()) particle->ResetCuts();
345 } 261 }
>> 262
>> 263 // inform that cut values are modified to the run manager
>> 264 // i.e state will be changed and SetCuts will be invoked
>> 265 // just before event loop
>> 266
>> 267 // set /control/verbose 0
>> 268 G4int tempVerboseLevel = G4UImanager::GetUIpointer()->GetVerboseLevel();
>> 269 G4UImanager::GetUIpointer()->SetVerboseLevel(0);
>> 270 // issue /run/cutoffModified
>> 271 G4UImanager::GetUIpointer()->ApplyCommand("/run/cutoffModified");
>> 272 // retreive /control/verbose
>> 273 G4UImanager::GetUIpointer()->SetVerboseLevel(tempVerboseLevel);
>> 274
346 } 275 }
347 276
348 // ------------------------------------------- << 277
349 void G4VUserPhysicsList::SetDefaultCutValue(G4 << 278 ////////////////////////////////////////////////////////
350 { << 279 void G4VUserPhysicsList::SetCutValueForOtherThan(G4double cutValue,
351 if (value < 0.0) { << 280 G4ParticleDefinition* first,
352 #ifdef G4VERBOSE << 281 G4ParticleDefinition* second,
353 if (verboseLevel > 0) { << 282 G4ParticleDefinition* third,
354 G4cout << "G4VUserPhysicsList::SetDefaul << 283 G4ParticleDefinition* fourth,
355 << " :" << value / mm << "[mm]" << 284 G4ParticleDefinition* fifth,
>> 285 G4ParticleDefinition* sixth,
>> 286 G4ParticleDefinition* seventh,
>> 287 G4ParticleDefinition* eighth,
>> 288 G4ParticleDefinition* nineth,
>> 289 G4ParticleDefinition* tenth )
>> 290 {
>> 291 // check cut value is positive
>> 292 if (cutValue <= 0.0) {
>> 293 #ifdef G4VERBOSE
>> 294 if (verboseLevel >0){
>> 295 G4cout << "G4VUserPhysicsList::SetCutValueForOtherThan: negative cut values";
>> 296 G4cout << " :" << cutValue/mm << "[mm]" << G4endl;
356 } 297 }
357 #endif 298 #endif
358 return; 299 return;
>> 300 } else {
>> 301 #ifdef G4VERBOSE
>> 302 if (verboseLevel >1) {
>> 303 G4cout << "G4VUserPhysicsList::SetCutValueForOtherThan ";
>> 304 G4cout << " :" << cutValue/mm << "[mm]" << G4endl;
>> 305 }
>> 306 #endif
359 } 307 }
360 308
361 defaultCutValue = value; << 309 // check specified particle types in arguments
362 isSetDefaultCutValue = true; << 310 G4ParticleDefinition* specifiedParticles[10];
>> 311 G4int numberOfSpecifiedParticles = 0;
>> 312 if (first != 0) {
>> 313 specifiedParticles[numberOfSpecifiedParticles] = first;
>> 314 numberOfSpecifiedParticles +=1;
>> 315 }
>> 316 if (second != 0) {
>> 317 specifiedParticles[numberOfSpecifiedParticles] = second;
>> 318 numberOfSpecifiedParticles +=1;
>> 319 }
>> 320 if (third != 0) {
>> 321 specifiedParticles[numberOfSpecifiedParticles] = third;
>> 322 numberOfSpecifiedParticles +=1;
>> 323 }
>> 324 if (fourth != 0) {
>> 325 specifiedParticles[numberOfSpecifiedParticles] = fourth;
>> 326 numberOfSpecifiedParticles +=1;
>> 327 }
>> 328 if (fifth != 0) {
>> 329 specifiedParticles[numberOfSpecifiedParticles] = fifth;
>> 330 numberOfSpecifiedParticles +=1;
>> 331 }
>> 332 if (sixth != 0) {
>> 333 specifiedParticles[numberOfSpecifiedParticles] = sixth;
>> 334 numberOfSpecifiedParticles +=1;
>> 335 }
>> 336 if (seventh != 0) {
>> 337 specifiedParticles[numberOfSpecifiedParticles] = seventh;
>> 338 numberOfSpecifiedParticles +=1;
>> 339 }
>> 340 if (eighth != 0) {
>> 341 specifiedParticles[numberOfSpecifiedParticles] = eighth;
>> 342 numberOfSpecifiedParticles +=1;
>> 343 }
>> 344 if (nineth != 0) {
>> 345 specifiedParticles[numberOfSpecifiedParticles] = nineth;
>> 346 numberOfSpecifiedParticles +=1;
>> 347 }
>> 348 if (tenth != 0) {
>> 349 specifiedParticles[numberOfSpecifiedParticles] = tenth;
>> 350 numberOfSpecifiedParticles +=1;
>> 351 }
363 352
364 // set cut values for gamma at first and for << 353 // Set cut values for other particles
365 SetCutValue(defaultCutValue, "gamma"); << 354 G4bool isSpecified;
366 SetCutValue(defaultCutValue, "e-"); << 355 theParticleIterator->reset();
367 SetCutValue(defaultCutValue, "e+"); << 356 while( (*theParticleIterator)() ){
368 SetCutValue(defaultCutValue, "proton"); << 357 G4ParticleDefinition* particle = theParticleIterator->value();
>> 358 isSpecified = false;
>> 359 for (G4int index = 0; index <numberOfSpecifiedParticles; index++) {
>> 360 isSpecified = isSpecified || (particle == specifiedParticles[index]);
>> 361 }
>> 362 if ( (!isSpecified) && (!particle->IsShortLived()) ){
>> 363 // set cut value
>> 364 SetParticleCuts( cutValue ,particle );
>> 365 // build physics table
>> 366 BuildPhysicsTable(particle);
369 367
370 #ifdef G4VERBOSE << 368 #ifdef G4VERBOSE
371 if (verboseLevel > 1) { << 369 if (verboseLevel >1) G4cout << "Set cuts for " << particle->GetParticleName() << G4endl;
372 G4cout << "G4VUserPhysicsList::SetDefaultC <<
373 << "default cut value is changed to <<
374 } <<
375 #endif 370 #endif
376 } <<
377 <<
378 // ------------------------------------------- <<
379 G4double G4VUserPhysicsList::GetCutValue(const <<
380 { <<
381 std::size_t nReg = (G4RegionStore::GetInstan <<
382 if (nReg == 0) { <<
383 #ifdef G4VERBOSE <<
384 if (verboseLevel > 0) { <<
385 G4cout << "G4VUserPhysicsList::GetCutVal <<
386 << " : No Default Region " << G4e <<
387 } 371 }
388 #endif <<
389 G4Exception("G4VUserPhysicsList::GetCutVal <<
390 return -1. * mm; <<
391 } 372 }
392 G4Region* region = G4RegionStore::GetInstanc <<
393 return region->GetProductionCuts()->GetProdu <<
394 } 373 }
395 374
396 // ------------------------------------------- <<
397 void G4VUserPhysicsList::SetCutValue(G4double <<
398 { <<
399 SetParticleCuts(aCut, name); <<
400 } <<
401 375
402 // ------------------------------------------- << 376 ////////////////////////////////////////////////////////
403 void G4VUserPhysicsList::SetCutValue(G4double << 377 void G4VUserPhysicsList::SetCutValue(G4double aCut, const G4String& name)
404 { 378 {
405 G4Region* region = G4RegionStore::GetInstanc << 379 G4ParticleDefinition* particle = theParticleTable->FindParticle(name);
406 if (region != nullptr) { << 380 #ifdef G4VERBOSE
407 // set cut value << 381 if (particle != 0){
408 SetParticleCuts(aCut, pname, region); << 382 if (verboseLevel >1) {
409 } << 383 G4cout << "G4VUserPhysicsList::SetCutValue :";
410 else { << 384 G4cout << "Set cuts for " << name << G4endl;
411 #ifdef G4VERBOSE << 385 }
412 if (verboseLevel > 0) { << 386 } else {
413 G4cout << "G4VUserPhysicsList::SetCutVal << 387 if (verboseLevel >0) {
414 << " : No Region of " << rname << << 388 G4cout << "G4VUserPhysicsList::SetCutValue :";
>> 389 G4cout << name << " is not found in ParticleTable" << G4endl;
415 } 390 }
416 #endif <<
417 } 391 }
>> 392 #endif
>> 393
>> 394 if (particle != 0){
>> 395 if (!particle->IsShortLived()) {
>> 396 //set cut value
>> 397 SetParticleCuts( aCut ,particle );
>> 398 // build physics table
>> 399 BuildPhysicsTable(particle);
>> 400 }
>> 401 }
418 } 402 }
419 403
420 // ------------------------------------------- << 404
>> 405
>> 406 ////////////////////////////////////////////////////////
421 void G4VUserPhysicsList::SetCutsWithDefault() 407 void G4VUserPhysicsList::SetCutsWithDefault()
422 { 408 {
423 SetDefaultCutValue(defaultCutValue); << 409 // default cut value
424 G4VUserPhysicsList::SetCuts(); << 410 G4double cut = defaultCutValue;
425 } <<
426 411
427 // ------------------------------------------- << 412 #ifdef G4VERBOSE
428 void G4VUserPhysicsList::SetCutsForRegion(G4do << 413 if (verboseLevel >1){
429 { << 414 G4cout << "G4VUserPhysicsList::SetCutsWithDefault:";
430 // set cut values for gamma at first and for << 415 G4cout << "CutLength : " << cut/mm << " (mm)" << G4endl;
431 SetCutValue(aCut, "gamma", rname); << 416 }
432 SetCutValue(aCut, "e-", rname); << 417 #endif
433 SetCutValue(aCut, "e+", rname); << 418
434 SetCutValue(aCut, "proton", rname); << 419 // set cut values for gamma at first and for e- second and next for e+,
435 } << 420 // because some processes for e+/e- need cut values for gamma
>> 421 SetCutValue(cut, "gamma");
>> 422 SetCutValue(cut, "e-");
>> 423 SetCutValue(cut, "e+");
>> 424
>> 425 // set cut values for proton and anti_proton before all other hadrons
>> 426 // because some processes for hadrons need cut values for proton/anti_proton
>> 427 SetCutValue(cut, "proton");
>> 428 SetCutValue(cut, "anti_proton");
>> 429
>> 430 SetCutValueForOthers(cut);
>> 431
>> 432 if (verboseLevel>1) {
>> 433 DumpCutValuesTable();
>> 434 }
>> 435 }
436 436
437 // ------------------------------------------- <<
438 void G4VUserPhysicsList::SetParticleCuts(G4dou <<
439 G4Reg <<
440 { <<
441 SetParticleCuts(cut, particle->GetParticleNa <<
442 } <<
443 437
444 // ------------------------------------------- << 438 ////////////////////////////////////////////////////////
445 void G4VUserPhysicsList::SetParticleCuts(G4dou << 439 void G4VUserPhysicsList::SetCutValueForOthers(G4double cutValue)
446 G4Reg <<
447 { 440 {
448 if (cut < 0.0) { << 441 // check cut value is positive
449 #ifdef G4VERBOSE << 442 if (cutValue <= 0.0) {
450 if (verboseLevel > 0) { << 443 #ifdef G4VERBOSE
451 G4cout << "G4VUserPhysicsList::SetPartic << 444 if (verboseLevel >0){
452 << " :" << cut / mm << "[mm]" << 445 G4cout << "G4VUserPhysicsList::SetCutValueForOthers: negative cut values";
453 << " for " << particleName << G4e << 446 G4cout << " :" << cutValue/mm << "[mm]" << G4endl;
454 } 447 }
455 #endif 448 #endif
456 return; 449 return;
457 } 450 }
458 451
459 G4Region* world_region = << 452 #ifdef G4VERBOSE
460 G4RegionStore::GetInstance()->GetRegion("D << 453 if (verboseLevel >1) {
461 if (region == nullptr) { << 454 G4cout << "G4VUserPhysicsList::SetCutValueForOthers ";
462 std::size_t nReg = G4RegionStore::GetInsta << 455 G4cout << " :" << cutValue/mm << "[mm]" << G4endl;
463 if (nReg == 0) { <<
464 #ifdef G4VERBOSE <<
465 if (verboseLevel > 0) { <<
466 G4cout << "G4VUserPhysicsList::SetPart <<
467 << " : No Default Region " << G <<
468 } <<
469 #endif <<
470 G4Exception("G4VUserPhysicsList::SetPart <<
471 "No Default Region"); <<
472 return; <<
473 } <<
474 region = world_region; <<
475 } <<
476 <<
477 if (!isSetDefaultCutValue) { <<
478 SetDefaultCutValue(defaultCutValue); <<
479 } <<
480 <<
481 G4ProductionCuts* pcuts = region->GetProduct <<
482 if (region != world_region <<
483 && pcuts == G4ProductionCutsTable::GetPr <<
484 { // This region had no unique cuts yet but <<
485 // Need to create a new object before set <<
486 pcuts = new G4ProductionCuts( <<
487 *(G4ProductionCutsTable::GetProductionCu <<
488 region->SetProductionCuts(pcuts); <<
489 } <<
490 pcuts->SetProductionCut(cut, particleName); <<
491 #ifdef G4VERBOSE <<
492 if (verboseLevel > 2) { <<
493 G4cout << "G4VUserPhysicsList::SetParticle <<
494 << " :" << cut / mm << "[mm]" <<
495 << " for " << particleName << G4end <<
496 } 456 }
497 #endif 457 #endif
498 } <<
499 458
500 // ------------------------------------------- << 459 // Sets a cut value to particle types which have not be called SetCuts()
501 void G4VUserPhysicsList::BuildPhysicsTable() <<
502 { <<
503 // Prepare Physics table for all particles <<
504 theParticleIterator->reset(); 460 theParticleIterator->reset();
505 while ((*theParticleIterator)()) { << 461 while( (*theParticleIterator)() ){
506 G4ParticleDefinition* particle = thePartic 462 G4ParticleDefinition* particle = theParticleIterator->value();
507 PreparePhysicsTable(particle); <<
508 } <<
509 463
510 // ask processes to prepare physics table << 464 if (!particle->IsShortLived()) {
511 if (fRetrievePhysicsTable) { << 465 // check if the cut value has already been set
512 fIsRestoredCutValues = fCutsTable->Retriev << 466 if ((particle->GetLengthCuts()<0.0) ||(particle->GetEnergyCuts()==0)) {
513 // check if retrieve Cut Table successfull << 467 // set cut value
514 if (!fIsRestoredCutValues) { << 468 SetParticleCuts( cutValue ,particle );
515 #ifdef G4VERBOSE << 469 // build physics table
516 if (verboseLevel > 0) { << 470 BuildPhysicsTable(particle);
517 G4cout << "G4VUserPhysicsList::BuildPh << 471
518 << " Retrieve Cut Table failed << 472 #ifdef G4VERBOSE
519 } << 473 if (verboseLevel >1)
>> 474 G4cout << "Set cuts for " << particle->GetParticleName() << G4endl;
520 #endif 475 #endif
521 G4Exception("G4VUserPhysicsList::BuildPh <<
522 "Fail to retrieve Production <<
523 } <<
524 else { <<
525 #ifdef G4VERBOSE <<
526 if (verboseLevel > 2) { <<
527 G4cout << "G4VUserPhysicsList::BuildPh <<
528 << " Retrieve Cut Table succes <<
529 } 476 }
530 #endif <<
531 } 477 }
>> 478
532 } 479 }
533 else { << 480 }
534 #ifdef G4VERBOSE <<
535 if (verboseLevel > 2) { <<
536 G4cout << "G4VUserPhysicsList::BuildPhys <<
537 << " does not retrieve Cut Table <<
538 } <<
539 #endif <<
540 } <<
541 481
542 // Sets a value to particle <<
543 // set cut values for gamma at first and for <<
544 G4String particleName; <<
545 G4ParticleDefinition* GammaP = theParticleTa <<
546 if (GammaP != nullptr) BuildPhysicsTable(Gam <<
547 G4ParticleDefinition* EMinusP = theParticleT <<
548 if (EMinusP != nullptr) BuildPhysicsTable(EM <<
549 G4ParticleDefinition* EPlusP = theParticleTa <<
550 if (EPlusP != nullptr) BuildPhysicsTable(EPl <<
551 G4ParticleDefinition* ProtonP = theParticleT <<
552 if (ProtonP != nullptr) BuildPhysicsTable(Pr <<
553 482
554 theParticleIterator->reset(); << 483 ////////////////////////////////////////////////////////
555 while ((*theParticleIterator)()) { << 484 void G4VUserPhysicsList::ReCalcCutValue(const G4String& name)
556 G4ParticleDefinition* particle = thePartic << 485 {
557 if (particle != GammaP && particle != EMin << 486 G4ParticleDefinition* particle = theParticleTable->FindParticle(name);
>> 487 if (particle != 0 ){
>> 488 if (!particle->IsShortLived()) {
>> 489 particle->ReCalcCuts();
558 BuildPhysicsTable(particle); 490 BuildPhysicsTable(particle);
>> 491
>> 492 #ifdef G4VERBOSE
>> 493 if (verboseLevel >1) G4cout << "Recalc cuts for " << name << G4endl;
>> 494 #endif
559 } 495 }
>> 496 }
>> 497
>> 498 #ifdef G4VERBOSE
>> 499 if (( particle ==0 ) && (verboseLevel >0)) {
>> 500 G4cout << name << " is not found in ParticleTable" << G4endl;
560 } 501 }
>> 502 #endif
561 503
562 // Set flag <<
563 fIsPhysicsTableBuilt = true; <<
564 } 504 }
565 505
566 // ------------------------------------------- << 506
567 void G4VUserPhysicsList::BuildPhysicsTable(G4P << 507 ////////////////////////////////////////////////////////
>> 508 void G4VUserPhysicsList::ReCalcCutValueForOthers()
568 { 509 {
569 if (auto* trackingManager = particle->GetTra << 510 #ifdef G4VERBOSE
570 #ifdef G4VERBOSE << 511 if (verboseLevel >1) {
571 if (verboseLevel > 2) { << 512 G4cout << "G4VUserPhysicsList::ReCalcCutValueForOthers ";
572 G4cout << "G4VUserPhysicsList::BuildPhys << 513 }
573 << "Calculate Physics Table for " <<
574 << " via custom TrackingManager" <<
575 } <<
576 #endif 514 #endif
577 trackingManager->BuildPhysicsTable(*partic << 515
578 return; << 516 // Sets a cut value to particle types which have not be called SetCuts()
>> 517 theParticleIterator->reset();
>> 518 while( (*theParticleIterator)() ){
>> 519
>> 520 G4ParticleDefinition* particle = theParticleIterator->value();
>> 521
>> 522 if (!particle->IsShortLived()) {
>> 523 if (particle->GetEnergyCuts()==0) {
>> 524 particle->ReCalcCuts();
>> 525 BuildPhysicsTable(particle);
>> 526
>> 527 #ifdef G4VERBOSE
>> 528 if (verboseLevel >1)
>> 529 G4cout << "ReCalc cuts for " << particle->GetParticleName() << G4endl;
>> 530 #endif
>> 531 }
>> 532 }
>> 533
579 } 534 }
>> 535 }
580 536
581 // Change in order to share physics tables f <<
582 537
583 if (particle->GetMasterProcessManager() == n << 538 ////////////////////////////////////////////////////////
584 #ifdef G4VERBOSE << 539 ////////////////////////////////////////////////////////
585 if (verboseLevel > 0) { << 540 void G4VUserPhysicsList::SetParticleCuts( G4double cut, G4ParticleDefinition* particle)
586 G4cout << "#### G4VUserPhysicsList::Buil << 541 {
587 << particle->GetParticleName() << << 542 if (fRetrievePhysicsTable && !fIsRestoredCutValues) {
>> 543 #ifdef G4VERBOSE
>> 544 if (verboseLevel>2){
>> 545 G4cout << "G4VUserPhysicsList::SetParticleCuts ";
>> 546 G4cout << " Retrieve Cut Values for ";
>> 547 G4cout << particle->GetParticleName() << G4endl;
588 } 548 }
589 #endif 549 #endif
590 return; << 550 RetrieveCutValues(directoryPhysicsTable, fStoredInAscii);
>> 551 fIsRestoredCutValues = true;
591 } 552 }
>> 553
>> 554 if (particle->GetEnergyCuts() == 0) {
>> 555 particle->SetCuts(cut);
>> 556 }
>> 557 }
>> 558
>> 559 ///////////////////////////////////////////////////////////////
>> 560 void G4VUserPhysicsList::BuildPhysicsTable(G4ParticleDefinition* particle)
>> 561 {
592 if (fRetrievePhysicsTable) { 562 if (fRetrievePhysicsTable) {
593 if (!fIsRestoredCutValues) { << 563 #ifdef G4VERBOSE
594 // fail to retrieve cut tables << 564 if (verboseLevel>2){
595 #ifdef G4VERBOSE << 565 G4cout << "G4VUserPhysicsList::BuildPhysicsTable ";
596 if (verboseLevel > 0) { << 566 G4cout << " Retrieve Physics Table for ";
597 G4cout << "G4VUserPhysicsList::BuildPh << 567 G4cout << particle->GetParticleName() << G4endl;
598 << "Physics table can not be re << 568 }
599 } <<
600 #endif <<
601 fRetrievePhysicsTable = false; <<
602 } <<
603 else { <<
604 #ifdef G4VERBOSE <<
605 if (verboseLevel > 2) { <<
606 G4cout << "G4VUserPhysicsList::BuildPh <<
607 << " Retrieve Physics Table for <<
608 } <<
609 #endif 569 #endif
>> 570 if (CheckForRetrievePhysicsTable(directoryPhysicsTable, fStoredInAscii)) {
610 // Retrieve PhysicsTable from files for 571 // Retrieve PhysicsTable from files for proccesses
611 RetrievePhysicsTable(particle, directory 572 RetrievePhysicsTable(particle, directoryPhysicsTable, fStoredInAscii);
612 } <<
613 } <<
614 573
615 #ifdef G4VERBOSE <<
616 if (verboseLevel > 2) { <<
617 G4cout << "G4VUserPhysicsList::BuildPhysic <<
618 << "Calculate Physics Table for " < <<
619 } <<
620 #endif <<
621 // Rebuild the physics tables for every proc <<
622 // if particle is not ShortLived <<
623 if (!particle->IsShortLived()) { <<
624 G4ProcessManager* pManager = particle->Get <<
625 if (pManager == nullptr) { <<
626 #ifdef G4VERBOSE <<
627 if (verboseLevel > 0) { <<
628 G4cout << "G4VUserPhysicsList::BuildPh <<
629 << " : No Process Manager for " <<
630 G4cout << particle->GetParticleName() <<
631 } <<
632 #endif <<
633 G4Exception("G4VUserPhysicsList::BuildPh <<
634 "No process manager"); <<
635 return; 574 return;
636 } << 575 } else {
637 << 576 #ifdef G4VERBOSE
638 // Get processes from master thread; << 577 if (verboseLevel>2){
639 G4ProcessManager* pManagerShadow = particl << 578 G4cout << "CheckForRetrievePhysicsTable failed " << G4endl;
640 <<
641 G4ProcessVector* pVector = pManager->GetPr <<
642 if (pVector == nullptr) { <<
643 #ifdef G4VERBOSE <<
644 if (verboseLevel > 0) { <<
645 G4cout << "G4VUserPhysicsList::BuildPh <<
646 << " : No Process Vector for " <<
647 } 579 }
648 #endif 580 #endif
649 G4Exception("G4VUserPhysicsList::BuildPh <<
650 "No process Vector"); <<
651 return; <<
652 } 581 }
653 #ifdef G4VERBOSE << 582 }
654 if (verboseLevel > 2) { <<
655 G4cout << "G4VUserPhysicsList::BuildPhys <<
656 << G4endl; <<
657 G4cout << " ProcessManager : " << pManag <<
658 << G4endl; <<
659 for (G4int iv1 = 0; iv1 < (G4int)pVector <<
660 G4cout << " " << iv1 << " - " << (*pV <<
661 } <<
662 G4cout << "----------------------------- <<
663 G4ProcessVector* pVectorShadow = pManage <<
664 583
665 for (G4int iv2 = 0; iv2 < (G4int)pVector << 584 #ifdef G4VERBOSE
666 G4cout << " " << iv2 << " - " << (*pV << 585 if (verboseLevel>2){
667 } << 586 G4cout << "G4VUserPhysicsList::BuildPhysicsTable ";
>> 587 G4cout << " for " << particle->GetParticleName() << G4endl;
668 } 588 }
669 #endif 589 #endif
670 for (G4int j = 0; j < (G4int)pVector->size << 590 G4int j;
671 // Andrea July 16th 2013 : migration to << 591 // Rebuild the physics tables for every process for this particle type
672 // Infer if we are in a worker thread or << 592 G4ProcessVector* pVector = (particle->GetProcessManager())->GetProcessList();
673 // Master thread is the one in which the << 593 for ( j=0; j < pVector->entries(); ++j) {
674 // and process manager shadow pointers a << 594 (*pVector)[j]->BuildPhysicsTable(*particle);
675 if (pManagerShadow == pManager) { <<
676 (*pVector)[j]->BuildPhysicsTable(*part <<
677 } <<
678 else { <<
679 (*pVector)[j]->BuildWorkerPhysicsTable <<
680 } <<
681 <<
682 } // End loop on processes vector <<
683 } // End if short-lived <<
684 } <<
685 <<
686 // ------------------------------------------- <<
687 void G4VUserPhysicsList::PreparePhysicsTable(G <<
688 { <<
689 if (auto* trackingManager = particle->GetTra <<
690 trackingManager->PreparePhysicsTable(*part <<
691 return; <<
692 } 595 }
693 << 596 for ( j=0; j < pVector->entries(); ++j) {
694 if ((particle->GetMasterProcessManager()) == << 597 // temporary addition to make the integral schema
695 return; << 598 BuildIntegralPhysicsTable((*pVector)[j], particle);
696 } 599 }
697 // Prepare the physics tables for every proc << 600 }
698 // if particle is not ShortLived <<
699 if (!particle->IsShortLived()) { <<
700 G4ProcessManager* pManager = particle->Get <<
701 if (pManager == nullptr) { <<
702 #ifdef G4VERBOSE <<
703 if (verboseLevel > 0) { <<
704 G4cout << "G4VUserPhysicsList::Prepare <<
705 << ": No Process Manager for " <<
706 G4cout << particle->GetParticleName() <<
707 } <<
708 #endif <<
709 G4Exception("G4VUserPhysicsList::Prepare <<
710 "No process manager"); <<
711 return; <<
712 } <<
713 <<
714 // Get processes from master thread <<
715 G4ProcessManager* pManagerShadow = particl <<
716 601
717 G4ProcessVector* pVector = pManager->GetPr << 602 ///////////////////////////////////////////////////////////////
718 if (pVector == nullptr) { << 603 void G4VUserPhysicsList::BuildIntegralPhysicsTable(G4VProcess* process,
719 #ifdef G4VERBOSE << 604 G4ParticleDefinition* particle)
720 if (verboseLevel > 0) { << 605 {
721 G4cout << "G4VUserPhysicsList::Prepare << 606 //*********************************************************************
722 << ": No Process Vector for " < << 607 // temporary addition to make the integral schema of electromagnetic
723 } <<
724 #endif <<
725 G4Exception("G4VUserPhysicsList::Prepare <<
726 "No process Vector"); <<
727 return; <<
728 } <<
729 for (G4int j = 0; j < (G4int)pVector->size <<
730 // Andrea July 16th 2013 : migration to <<
731 // Infer if we are in a worker thread or <<
732 // Master thread is the one in which the <<
733 // and process manager shadow pointers a <<
734 if (pManagerShadow == pManager) { <<
735 (*pVector)[j]->PreparePhysicsTable(*pa <<
736 } <<
737 else { <<
738 (*pVector)[j]->PrepareWorkerPhysicsTab <<
739 } <<
740 } // End loop on processes vector <<
741 } // End if pn ShortLived <<
742 } <<
743 <<
744 // ------------------------------------------- <<
745 void G4VUserPhysicsList::BuildIntegralPhysicsT <<
746 <<
747 { <<
748 // TODO Should we change this function? <<
749 //****************************************** <<
750 // Temporary addition to make the integral s <<
751 // processes work. 608 // processes work.
>> 609 //
752 610
753 if ((process->GetProcessName() == "Imsc") || << 611 if ( (process->GetProcessName() == "Imsc") ||
754 || (process->GetProcessName() == "IeBrem << 612 (process->GetProcessName() == "IeIoni") ||
755 || (process->GetProcessName() == "IhIoni << 613 (process->GetProcessName() == "IeBrems") ||
756 || (process->GetProcessName() == "IMuBre << 614 (process->GetProcessName() == "Iannihil") ||
757 { << 615 (process->GetProcessName() == "IhIoni") ||
758 #ifdef G4VERBOSE << 616 (process->GetProcessName() == "IMuIoni") ||
759 if (verboseLevel > 2) { << 617 (process->GetProcessName() == "IMuBrems") ||
760 G4cout << "G4VUserPhysicsList::BuildInte << 618 (process->GetProcessName() == "IMuPairProd") ) {
761 << " BuildPhysicsTable is invoked << 619 #ifdef G4VERBOSE
762 << particle->GetParticleName() << << 620 if (verboseLevel>2){
>> 621 G4cout << "G4VUserPhysicsList::BuildIntegralPhysicsTable ";
>> 622 G4cout << " BuildPhysicsTable is invoked for ";
>> 623 G4cout << process->GetProcessName();
>> 624 G4cout << "(" << particle->GetParticleName() << ")" << G4endl;
763 } 625 }
764 #endif 626 #endif
765 process->BuildPhysicsTable(*particle); 627 process->BuildPhysicsTable(*particle);
766 } 628 }
767 } 629 }
768 630
769 // ------------------------------------------- << 631 ///////////////////////////////////////////////////////////////
770 void G4VUserPhysicsList::DumpList() const 632 void G4VUserPhysicsList::DumpList() const
771 { 633 {
772 theParticleIterator->reset(); 634 theParticleIterator->reset();
773 G4int idx = 0; 635 G4int idx = 0;
774 while ((*theParticleIterator)()) { << 636 while( (*theParticleIterator)() ){
775 G4ParticleDefinition* particle = thePartic 637 G4ParticleDefinition* particle = theParticleIterator->value();
776 G4cout << particle->GetParticleName(); 638 G4cout << particle->GetParticleName();
777 if ((idx++ % 4) == 3) { 639 if ((idx++ % 4) == 3) {
778 G4cout << G4endl; 640 G4cout << G4endl;
779 } << 641 } else {
780 else { <<
781 G4cout << ", "; 642 G4cout << ", ";
782 } << 643 }
783 } 644 }
784 G4cout << G4endl; 645 G4cout << G4endl;
785 } 646 }
786 647
787 // ------------------------------------------- << 648 ///////////////////////////////////////////////////////////////
788 void G4VUserPhysicsList::DumpCutValuesTable(G4 << 649 void G4VUserPhysicsList::DumpCutValues(const G4String &particle_name) const
789 { 650 {
790 fDisplayThreshold = flag; << 651 G4ParticleDefinition* particle;
791 } << 652 if ((particle_name == "ALL") || (particle_name == "all")) {
>> 653 theParticleIterator->reset();
>> 654 while( (*theParticleIterator)() ){
>> 655 particle = theParticleIterator->value();
>> 656 DumpCutValues(particle);
>> 657 }
>> 658 } else {
>> 659 particle = theParticleTable->FindParticle(particle_name);
>> 660 if (particle != 0) DumpCutValues(particle);
>> 661 }
>> 662 }
>> 663
>> 664 ///////////////////////////////////////////////////////////////
>> 665 void G4VUserPhysicsList::DumpCutValues( G4ParticleDefinition* particle) const
>> 666 {
>> 667 if (particle == 0) return;
>> 668
>> 669 G4int prec = G4cout.precision(3);
>> 670
>> 671 if (particle->IsShortLived()) {
>> 672 // name field
>> 673 G4cout << " --- " << particle->GetParticleName() << " is a short lived particle ------ " << G4endl;
>> 674 } else {
>> 675 // name field
>> 676 G4cout << " --- " << particle->GetParticleName() << " ------ " << G4endl;
>> 677
>> 678 // cut value in range field
>> 679 G4cout << " - Cut in range = " << G4BestUnit(particle->GetLengthCuts(),"Length") << G4endl;
>> 680
>> 681 // material and energy cut value for the material
>> 682 G4double* theKineticEnergyCuts = particle->GetEnergyCuts();
>> 683
>> 684 if (theKineticEnergyCuts != 0) {
>> 685 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
>> 686 G4cout << " - Material ---------------- Energy Cut ---" << G4endl;
>> 687 for (G4int idx=0; idx<materialTable->entries(); idx++){
>> 688 G4cout << " " << G4std::setw(19) << (*materialTable)[idx]->GetName();
>> 689 G4cout << " : " << G4std::setw(10) << G4BestUnit(theKineticEnergyCuts[idx],"Energy");
>> 690 G4cout << G4endl;
>> 691 }
>> 692
>> 693 } else {
>> 694 G4cout << " - Cuts in energy are not calculated yet --" << G4endl;
>> 695 G4cout << " Enter /run/initialize command to calculate cuts " << G4endl;
>> 696 }
>> 697
>> 698 }
>> 699 G4cout.precision(prec);
>> 700 }
>> 701
>> 702 ///////////////////////////////////////////////////////////////
>> 703 void G4VUserPhysicsList::DumpCutValuesTable() const
>> 704 {
>> 705 // This methods Print out a table of cut value information
>> 706 // for "e-", "gamma", "mu-", "proton" and "neutron"
>> 707 G4int prec = G4cout.precision(3);
>> 708 const G4int size = 5;
>> 709 G4String name[size] = {"gamma", "e-", "mu-", "proton", "neutron"};
>> 710 G4ParticleDefinition* particle[size];
>> 711 G4bool IsOK = true;
>> 712 G4int size_display=2;
>> 713 G4int idx;
>> 714 for (idx=0; idx <size_display; idx++) {
>> 715 particle[idx] = theParticleTable->FindParticle(name[idx]);
>> 716 }
>> 717
>> 718 //line 1 //-- Commented out (M.Asai)
>> 719 //G4cout << "Default cut value in range :";
>> 720 //G4cout << defaultCutValue/mm << "[mm]" << G4endl;
>> 721
>> 722 //line 2
>> 723 G4cout << "============= The cut Energy ==============================" <<G4endl;
>> 724
>> 725 // line 3
>> 726 G4cout << " ";
>> 727 for (idx=0; idx <size_display; idx++) {
>> 728 G4cout << " " << G4std::setw(11) << name[idx] << " ";
>> 729 }
>> 730 G4cout << G4endl;
792 731
793 // ------------------------------------------- << 732 // line 4
794 void G4VUserPhysicsList::DumpCutValuesTableIfR << 733 G4cout << "Cut in range ";
795 { << 734 for (idx=0; idx <size_display; idx++) {
796 if (fDisplayThreshold == 0) return; << 735 if (particle[idx] == 0) {
797 G4ProductionCutsTable::GetProductionCutsTabl << 736 G4cout << " ";
798 fDisplayThreshold = 0; << 737 } else {
>> 738 G4cout << " " << G4std::setw(11) << G4BestUnit(particle[idx]->GetLengthCuts(),"Length");
>> 739 }
>> 740 }
>> 741 G4cout << G4endl;
>> 742
>> 743 // line 5
>> 744 G4cout << "Cut in energy";
>> 745 G4cout << G4endl;
>> 746
>> 747 // line 6 ..
>> 748 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
>> 749 for (G4int J=0; J<materialTable->entries(); J++) {
>> 750 G4cout << " " << G4std::setw(18) << ((*materialTable)[J])->GetName();
>> 751 for (idx=0; idx <size_display; idx++) {
>> 752 if (particle[idx] == 0) {
>> 753 G4cout << " ";
>> 754 } else {
>> 755 if (particle[idx]->GetEnergyCuts() == 0) {
>> 756 G4cout << " ---------- ";
>> 757 IsOK = false;
>> 758 } else {
>> 759 G4cout << " " << G4std::setw(11) << G4BestUnit((particle[idx]->GetEnergyCuts())[J],"Energy");
>> 760 }
>> 761 }
>> 762 }
>> 763 G4cout << G4endl;
>> 764 }
>> 765
>> 766 if (!IsOK) {
>> 767 G4cout << " Cuts in energy have not calculated yet !!" << G4endl;
>> 768 G4cout << " Enter /run/initialize command to calculate cuts " << G4endl;
>> 769 }
>> 770
>> 771 // last line
>> 772 G4cout << "===================================================" << G4endl;
>> 773 G4cout.precision(prec);
799 } 774 }
800 775
801 // ------------------------------------------- << 776
>> 777 ///////////////////////////////////////////////////////////////
>> 778 ///////////////////////////////////////////////////////////////
802 G4bool G4VUserPhysicsList::StorePhysicsTable(c 779 G4bool G4VUserPhysicsList::StorePhysicsTable(const G4String& directory)
803 { 780 {
804 G4bool ascii = fStoredInAscii; << 781 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
805 G4String dir = directory; << 782 numberOfMaterial = matTable->entries();
806 if (dir.empty()) <<
807 dir = directoryPhysicsTable; <<
808 else <<
809 directoryPhysicsTable = dir; <<
810 783
811 // store CutsTable info << 784 G4bool ascii = fStoredInAscii;
812 if (!fCutsTable->StoreCutsTable(dir, ascii)) << 785 G4String dir = directory;
813 G4Exception("G4VUserPhysicsList::StorePhys << 786 if (dir.isNull()) dir = directoryPhysicsTable;
814 "Fail to store Cut Table"); << 787 if (!StoreMaterialInfo(dir, ascii)) return false;
815 return false; << 788 if (!StoreCutValues(dir,ascii)) return false;
816 } << 789 #ifdef G4VERBOSE
817 #ifdef G4VERBOSE << 790 if (verboseLevel>2){
818 if (verboseLevel > 2) { << 791 G4cout << "G4VUserPhysicsList::StorePhysicsTable ";
819 G4cout << "G4VUserPhysicsList::StorePhysic << 792 G4cout << " Store material and cu values successfully" << G4endl;
820 << " Store material and cut values <<
821 } 793 }
822 #endif 794 #endif
823 795
824 G4bool success = true; << 796 G4bool success= true;
825 797
826 // loop over all particles in G4ParticleTabl << 798 // loop over all particles in G4ParticleTable
827 theParticleIterator->reset(); 799 theParticleIterator->reset();
828 while ((*theParticleIterator)()) { << 800 while( (*theParticleIterator)() ){
829 G4ParticleDefinition* particle = thePartic 801 G4ParticleDefinition* particle = theParticleIterator->value();
830 // Store physics tables for every process 802 // Store physics tables for every process for this particle type
831 G4ProcessVector* pVector = (particle->GetP 803 G4ProcessVector* pVector = (particle->GetProcessManager())->GetProcessList();
832 for (G4int j = 0; j < (G4int)pVector->size << 804 G4int j;
833 if (!(*pVector)[j]->StorePhysicsTable(pa << 805 for ( j=0; j < pVector->entries(); ++j) {
834 G4String comment = "Fail to store phys << 806 if (!(*pVector)[j]->StorePhysicsTable(directoryPhysicsTable,ascii)){
835 comment += (*pVector)[j]->GetProcessNa << 807 #ifdef G4VERBOSE
836 comment += "(" + particle->GetParticle << 808 if (verboseLevel>2){
837 G4Exception("G4VUserPhysicsList::Store << 809 G4cout << "G4VUserPhysicsList::StorePhysicsTable ";
838 success = false; << 810 G4cout << " Fail to retrieve for ";
>> 811 G4cout << (*pVector)[j]->GetProcessName();
>> 812 G4cout << "(" << particle->GetParticleName() <<")" << G4endl;
>> 813 }
>> 814 #endif
>> 815 success = false;
839 } 816 }
840 } 817 }
841 // end loop over processes 818 // end loop over processes
842 } 819 }
843 // end loop over particles 820 // end loop over particles
844 return success; 821 return success;
845 } 822 }
846 823
847 // ------------------------------------------- << 824
848 void G4VUserPhysicsList::SetPhysicsTableRetrie << 825 ///////////////////////////////////////////////////////////////
>> 826 G4bool G4VUserPhysicsList::StoreMaterialInfo(const G4String& directory,
>> 827 G4bool ascii)
>> 828 {
>> 829 const G4String fileName = directory + "/" + "material.dat";
>> 830 const G4String key = "MATERIAL";
>> 831 G4std::ofstream fOut;
>> 832
>> 833 // open output file //
>> 834 #ifdef G4USE_STD_NAMESPACE
>> 835 if (!ascii )
>> 836 fOut.open(fileName,G4std::ios::out|G4std::ios::binary);
>> 837 else
>> 838 #endif
>> 839 fOut.open(fileName,G4std::ios::out);
>> 840
>> 841
>> 842 // check if the file has been opened successfully
>> 843 if (!fOut) {
>> 844 #ifdef G4VERBOSE
>> 845 G4cerr << "G4VUserPhysicsList::StoreMaterialInfo ";
>> 846 G4cerr << " Can not open file " << fileName << G4endl;
>> 847 #endif
>> 848 return false;
>> 849 }
>> 850
>> 851 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
>> 852 // number of materials in the table
>> 853 numberOfMaterial = matTable->entries();
>> 854
>> 855 if (ascii) {
>> 856 /////////////// ASCII mode /////////////////
>> 857 // key word
>> 858 fOut << key << G4endl;
>> 859
>> 860 // number of materials in the table
>> 861 fOut << numberOfMaterial << G4endl;
>> 862
>> 863 fOut.setf(G4std::ios::scientific);
>> 864
>> 865 // material name and density
>> 866 for (size_t idx=0; idx<matTable->entries(); ++idx){
>> 867 fOut << G4std::setw(FixedStringLengthForStore) << ((*matTable)[idx])->GetName();
>> 868 fOut << G4std::setw(FixedStringLengthForStore) << ((*matTable)[idx])->GetDensity()/(g/cm3) << G4endl;
>> 869 }
>> 870
>> 871 fOut.unsetf(G4std::ios::scientific);
>> 872
>> 873 } else {
>> 874 /////////////// Binary mode /////////////////
>> 875 char temp[FixedStringLengthForStore];
>> 876 size_t i;
>> 877
>> 878 // key word
>> 879 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
>> 880 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i) temp[i]=key[i];
>> 881 fOut.write(temp, FixedStringLengthForStore);
>> 882
>> 883 // number of materials in the table
>> 884 fOut.write( (char*)(&numberOfMaterial), sizeof (G4int));
>> 885
>> 886 // material name and density
>> 887 for (size_t imat=0; imat<matTable->entries(); ++imat){
>> 888 G4String name = ((*matTable)[imat])->GetName();
>> 889 G4double density = ((*matTable)[imat])->GetDensity();
>> 890 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
>> 891 for (i=0; i<name.length() && i<FixedStringLengthForStore-1; ++i) temp[i]=name[i];
>> 892 fOut.write(temp, FixedStringLengthForStore);
>> 893 fOut.write( (char*)(&density), sizeof (G4double));
>> 894 }
>> 895 }
>> 896
>> 897 fOut.close();
>> 898 return true;
>> 899 }
>> 900
>> 901 ///////////////////////////////////////////////////////////////
>> 902 G4bool G4VUserPhysicsList::StoreCutValues(const G4String& directory,
>> 903 G4bool ascii)
>> 904 {
>> 905 const G4String fileName = directory + "/" + "cut_value.dat";
>> 906 const G4String key = "CUT_VALUE";
>> 907 G4std::ofstream fOut;
>> 908
>> 909 // open output file //
>> 910 #ifdef G4USE_STD_NAMESPACE
>> 911 if (!ascii )
>> 912 fOut.open(fileName,G4std::ios::out|G4std::ios::binary);
>> 913 else
>> 914 #endif
>> 915 fOut.open(fileName,G4std::ios::out);
>> 916
>> 917 // check if the file has been opened successfully
>> 918 if (!fOut) {
>> 919 #ifdef G4VERBOSE
>> 920 G4cerr << "G4VUserPhysicsList::StoreCutValues ";
>> 921 G4cerr << " Can not open file " << fileName << G4endl;
>> 922 #endif
>> 923 return false;
>> 924 }
>> 925
>> 926 char temp[FixedStringLengthForStore];
>> 927 size_t i;
>> 928 if (ascii) {
>> 929 /////////////// ASCII mode /////////////////
>> 930 // key word
>> 931 fOut << key << G4endl;
>> 932
>> 933 fOut.setf(G4std::ios::scientific);
>> 934 // default cut value
>> 935 fOut << "Default " << G4std::setw(20) << defaultCutValue/mm << G4endl;
>> 936 fOut.unsetf(G4std::ios::scientific);
>> 937 } else {
>> 938 /////////////// Binary mode /////////////////
>> 939 // key word
>> 940 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
>> 941 for (i=0; i<key.length() && i<FixedStringLengthForStore-1; ++i) temp[i]=key[i];
>> 942 fOut.write(temp, FixedStringLengthForStore);
>> 943
>> 944 // default cut value
>> 945 fOut.write( (char*)(&defaultCutValue), sizeof (G4double));
>> 946 }
>> 947
>> 948 // loop over all particles in G4ParticleTable
>> 949 theParticleIterator->reset();
>> 950 while( (*theParticleIterator)() ){
>> 951 G4ParticleDefinition* particle = theParticleIterator->value();
>> 952 // Only following particles are concerned
>> 953 if ( (particle->GetParticleName() == "gamma" ) ||
>> 954 (particle->GetParticleName() == "e-" ) ||
>> 955 (particle->GetParticleName() == "e+" ) ||
>> 956 (particle->GetParticleName() == "mu-" ) ||
>> 957 (particle->GetParticleName() == "mu+" ) ||
>> 958 (particle->GetParticleName() == "proton" ) ||
>> 959 (particle->GetParticleName() == "anti_proton" ) ||
>> 960 (particle->GetParticleName() == "neutron" ) ||
>> 961 (particle->GetParticleName() == "anti_neutron" ) ){
>> 962
>> 963
>> 964 // particle name and cut in length
>> 965 if (ascii) {
>> 966 /////////////// ASCII mode /////////////////
>> 967 fOut.setf(G4std::ios::scientific);
>> 968 fOut << G4std::setw(FixedStringLengthForStore) << particle->GetParticleName();
>> 969 fOut << G4std::setw(20) << particle->GetLengthCuts()/mm << G4endl;
>> 970 fOut.unsetf(G4std::ios::scientific);
>> 971 } else {
>> 972 /////////////// Binary mode /////////////////
>> 973 G4String name = particle->GetParticleName();
>> 974 G4double cutLength = particle->GetLengthCuts();
>> 975 for (i=0; i<FixedStringLengthForStore; ++i) temp[i] = '\0';
>> 976 for (i=0; i<name.length() && i<FixedStringLengthForStore-1; ++i) temp[i]=name[i];
>> 977 fOut.write(temp, FixedStringLengthForStore);
>> 978 fOut.write((char*)(&cutLength), sizeof (G4double));
>> 979 }
>> 980
>> 981 G4double* cutArray = particle->GetEnergyCuts();
>> 982 if (cutArray ==0) {
>> 983 #ifdef G4VERBOSE
>> 984 if (verboseLevel>0){
>> 985 G4cout << "G4VUserPhysicsList::StoreCutValues " ;
>> 986 G4cout << ": Energy Cut Values have not be calculated for";
>> 987 G4cout << particle->GetParticleName() << G4endl;
>> 988 }
>> 989 #endif
>> 990 return false;
>> 991 }
>> 992
>> 993 // cut energy for each material
>> 994 if (ascii) {
>> 995 /////////////// ASCII mode /////////////////
>> 996 fOut.setf(G4std::ios::scientific);
>> 997 G4int jj =0;
>> 998 for(size_t idx=0; idx<numberOfMaterial; ++idx, ++jj) {
>> 999 if (jj==4) {
>> 1000 fOut << G4endl;
>> 1001 jj =0;
>> 1002 }
>> 1003 fOut << G4std::setw(20) << cutArray[idx]/keV;
>> 1004 }
>> 1005 fOut << G4endl;
>> 1006 fOut.unsetf(G4std::ios::scientific);
>> 1007 } else {
>> 1008 /////////////// Binary mode /////////////////
>> 1009 fOut.write( (char*)(cutArray), numberOfMaterial*(sizeof (G4double)) );
>> 1010 }
>> 1011 }
>> 1012 }
>> 1013
>> 1014 fOut.close();
>> 1015 return true;
>> 1016 }
>> 1017
>> 1018 ///////////////////////////////////////////////////////////////
>> 1019 void G4VUserPhysicsList::SetPhysicsTableRetrieved(const G4String& directory)
849 { 1020 {
850 fRetrievePhysicsTable = true; 1021 fRetrievePhysicsTable = true;
851 if (!directory.empty()) { << 1022 if(!directory.isNull()) {
852 directoryPhysicsTable = directory; 1023 directoryPhysicsTable = directory;
853 } 1024 }
854 fIsCheckedForRetrievePhysicsTable = false; << 1025 fIsCheckedForRetrievePhysicsTable=false;
855 fIsRestoredCutValues = false; 1026 fIsRestoredCutValues = false;
>> 1027 ResetCuts();
856 } 1028 }
857 1029
858 // ------------------------------------------- << 1030 ///////////////////////////////////////////////////////////////
859 void G4VUserPhysicsList::RetrievePhysicsTable( << 1031 void G4VUserPhysicsList::RetrievePhysicsTable(G4ParticleDefinition* particle,
860 << 1032 const G4String& directory,
>> 1033 G4bool ascii)
861 { 1034 {
862 G4bool success[100]; << 1035 G4int j;
>> 1036 G4bool success[100];
863 // Retrieve physics tables for every process 1037 // Retrieve physics tables for every process for this particle type
864 G4ProcessVector* pVector = (particle->GetPro 1038 G4ProcessVector* pVector = (particle->GetProcessManager())->GetProcessList();
865 for (G4int j = 0; j < (G4int)pVector->size() << 1039 for ( j=0; j < pVector->entries(); ++j) {
866 success[j] = (*pVector)[j]->RetrievePhysic << 1040 success[j] =
>> 1041 (*pVector)[j]->RetrievePhysicsTable(directory,ascii);
867 1042
868 if (!success[j]) { 1043 if (!success[j]) {
869 #ifdef G4VERBOSE << 1044 #ifdef G4VERBOSE
870 if (verboseLevel > 2) { << 1045 if (verboseLevel>2){
871 G4cout << "G4VUserPhysicsList::Retriev << 1046 G4cout << "G4VUserPhysicsList::RetrievePhysicsTable ";
872 << " Fail to retrieve Physics T << 1047 G4cout << " Fail to retrieve for ";
873 << G4endl; << 1048 G4cout << (*pVector)[j]->GetProcessName();
874 G4cout << "Calculate Physics Table for << 1049 G4cout << "(" << particle->GetParticleName() <<")" << G4endl;
875 } 1050 }
876 #endif 1051 #endif
877 (*pVector)[j]->BuildPhysicsTable(*partic 1052 (*pVector)[j]->BuildPhysicsTable(*particle);
878 } 1053 }
879 } 1054 }
880 for (G4int j = 0; j < (G4int)pVector->size() << 1055 for ( j=0; j < pVector->entries(); ++j) {
881 // temporary addition to make the integral 1056 // temporary addition to make the integral schema
882 if (!success[j]) BuildIntegralPhysicsTable << 1057 if (!success[j]) BuildIntegralPhysicsTable((*pVector)[j], particle);
883 } 1058 }
884 } 1059 }
885 1060
886 // ------------------------------------------- << 1061 ///////////////////////////////////////////////////////////////
887 void G4VUserPhysicsList::SetApplyCuts(G4bool v << 1062 G4bool G4VUserPhysicsList::CheckForRetrievePhysicsTable(const G4String& directory,
>> 1063 G4bool ascii)
888 { 1064 {
889 #ifdef G4VERBOSE << 1065 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
890 if (verboseLevel > 2) { << 1066 numberOfMaterial = matTable->entries();
891 G4cout << "G4VUserPhysicsList::SetApplyCut << 1067
>> 1068 if (!fIsCheckedForRetrievePhysicsTable) {
>> 1069 if (CheckMaterialInfo(directory, ascii)) {
>> 1070 if (CheckCutValues(directory,ascii))
>> 1071 fIsCheckedForRetrievePhysicsTable = true;
>> 1072 }
892 } 1073 }
>> 1074 return fIsCheckedForRetrievePhysicsTable;
>> 1075 }
>> 1076
>> 1077 ///////////////////////////////////////////////////////////////
>> 1078 G4bool G4VUserPhysicsList::CheckMaterialInfo(const G4String& directory,
>> 1079 G4bool ascii)
>> 1080 {
>> 1081 const G4String fileName = directory + "/" + "material.dat";
>> 1082 const G4String key = "MATERIAL";
>> 1083 G4std::ifstream fIn;
>> 1084
>> 1085 // open input file //
>> 1086 #ifdef G4USE_STD_NAMESPACE
>> 1087 if (!ascii )
>> 1088 fIn.open(fileName,G4std::ios::in|G4std::ios::binary);
>> 1089 else
>> 1090 #endif
>> 1091 fIn.open(fileName,G4std::ios::in);
>> 1092
>> 1093 // check if the file has been opened successfully
>> 1094 if (!fIn) {
>> 1095 #ifdef G4VERBOSE
>> 1096 G4cerr << "G4VUserPhysicsList::CheckMaterialInfo ";
>> 1097 G4cerr << " Can not open file " << fileName << G4endl;
893 #endif 1098 #endif
894 if (name == "all") { << 1099 return false;
895 theParticleTable->FindParticle("gamma")->S <<
896 theParticleTable->FindParticle("e-")->SetA <<
897 theParticleTable->FindParticle("e+")->SetA <<
898 theParticleTable->FindParticle("proton")-> <<
899 } 1100 }
900 else { << 1101
901 theParticleTable->FindParticle(name)->SetA << 1102 char temp[FixedStringLengthForStore];
>> 1103 size_t i;
>> 1104
>> 1105 // key word
>> 1106 G4String keyword;
>> 1107 if (ascii) {
>> 1108 fIn >> keyword;
>> 1109 } else {
>> 1110 fIn.read(temp, FixedStringLengthForStore);
>> 1111 keyword = (const char*)(temp);
>> 1112 }
>> 1113 if (key!=keyword) {
>> 1114 #ifdef G4VERBOSE
>> 1115 if (verboseLevel>0){
>> 1116 G4cout << "G4VUserPhysicsList::CheckMaterialInfo ";
>> 1117 G4cout << " Can not find key word " << keyword << G4endl;
>> 1118 }
>> 1119 #endif
>> 1120 return false;
902 } 1121 }
903 } <<
904 1122
905 // ------------------------------------------- << 1123 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
906 G4bool G4VUserPhysicsList::GetApplyCuts(const << 1124 numberOfMaterial = matTable->entries();
907 { << 1125 // number of materials in the table
908 return theParticleTable->FindParticle(name)- << 1126 G4int nmat;
909 } << 1127 if (ascii) {
>> 1128 fIn >> nmat;
>> 1129 } else {
>> 1130 fIn.read( (char*)(&nmat), sizeof (G4int));
>> 1131 }
>> 1132 if (nmat!=numberOfMaterial) {
>> 1133 #ifdef G4VERBOSE
>> 1134 if (verboseLevel>0){
>> 1135 G4cout << "G4VUserPhysicsList::CheckMaterialInfo ";
>> 1136 G4cout << "Number of material is inconsistent "<< G4endl;
>> 1137 }
>> 1138 #endif
>> 1139 return false;
>> 1140 }
>> 1141
>> 1142 // list of material
>> 1143 for (G4int idx=0; idx<matTable->entries(); ++idx){
>> 1144 // check eof
>> 1145 if(fIn.eof()) {
>> 1146 #ifdef G4VERBOSE
>> 1147 if (verboseLevel>0){
>> 1148 G4cout << "G4VUserPhysicsList::CheckMaterialInfo ";
>> 1149 G4cout << " encountered End of File" << G4endl;
>> 1150 }
>> 1151 #endif
>> 1152 fIn.close();
>> 1153 return false;
>> 1154 }
>> 1155
>> 1156 // check material name and density
>> 1157 char name[FixedStringLengthForStore];
>> 1158 double density;
>> 1159 if (ascii) {
>> 1160 fIn >> name >> density;
>> 1161 } else {
>> 1162 fIn.read(name, FixedStringLengthForStore);
>> 1163 fIn.read((char*)(&density), sizeof (G4double));
>> 1164 }
>> 1165 if (fIn.fail()) {
>> 1166 #ifdef G4VERBOSE
>> 1167 if (verboseLevel>0){
>> 1168 G4cout << "G4VUserPhysicsList::CheckMaterialInfo ";
>> 1169 G4cout << " Bad data format " << G4endl;
>> 1170 }
>> 1171 #endif
>> 1172 fIn.close();
>> 1173 return false;
>> 1174 }
>> 1175 G4double ratio = abs( (density*g/cm3)/((*matTable)[idx])->GetDensity() );
>> 1176 if ( (name != ((*matTable)[idx])->GetName()) || (0.999>ratio) || (ratio>1.001) ){
>> 1177 #ifdef G4VERBOSE
>> 1178 if (verboseLevel>0){
>> 1179 G4cout << "G4VUserPhysicsList::CheckMaterialInfo ";
>> 1180 G4cout << " Inconsistent material name or density" << G4endl;;
>> 1181 G4cout << G4std::setw(40) << name;
>> 1182 G4cout << G4std::setw(20) << G4std::setiosflags(G4std::ios::scientific) << density << G4endl;
>> 1183 G4cout << G4std::resetiosflags(G4std::ios::scientific);
>> 1184 }
>> 1185 #endif
>> 1186 fIn.close();
>> 1187 return false;
>> 1188 }
>> 1189 }
>> 1190 fIn.close();
>> 1191 return true;
>> 1192 }
>> 1193
>> 1194 ///////////////////////////////////////////////////////////////
>> 1195 G4bool G4VUserPhysicsList::CheckCutValues(const G4String& directory,
>> 1196 G4bool ascii)
>> 1197 {
>> 1198 const G4String fileName = directory + "/" + "cut_value.dat";
>> 1199 const G4String key = "CUT_VALUE";
>> 1200 G4std::ifstream fIn;
>> 1201
>> 1202 // open input file //
>> 1203 #ifdef G4USE_STD_NAMESPACE
>> 1204 if (!ascii )
>> 1205 fIn.open(fileName,G4std::ios::in|G4std::ios::binary);
>> 1206 else
>> 1207 #endif
>> 1208 fIn.open(fileName,G4std::ios::in);
910 1209
911 // ------------------------------------------- << 1210 // check if the file has been opened successfully
912 void G4VUserPhysicsList::CheckParticleList() << 1211 if (!fIn) {
913 { << 1212 #ifdef G4VERBOSE
914 if (!fDisableCheckParticleList) { << 1213 G4cerr << "G4VUserPhysicsList::CheckCutValues ";
915 G4MT_thePLHelper->CheckParticleList(); << 1214 G4cerr << " Can not open file " << fileName << G4endl;
>> 1215 #endif
>> 1216 return false;
916 } 1217 }
917 } <<
918 1218
919 // ------------------------------------------- << 1219 char temp[FixedStringLengthForStore];
920 void G4VUserPhysicsList::AddTransportation() << 1220 size_t i;
921 { <<
922 G4MT_thePLHelper->AddTransportation(); <<
923 } <<
924 1221
925 // ------------------------------------------- << 1222 // key word
926 void G4VUserPhysicsList::UseCoupledTransportat << 1223 G4String keyword;
927 { << 1224 if (ascii) {
928 G4MT_thePLHelper->UseCoupledTransportation(v << 1225 fIn >> keyword;
929 } << 1226 } else {
>> 1227 fIn.read(temp, FixedStringLengthForStore);
>> 1228 keyword = (const char*)(temp);
>> 1229 }
>> 1230 if (key!=keyword) {
>> 1231 #ifdef G4VERBOSE
>> 1232 if (verboseLevel>0){
>> 1233 G4cout << "G4VUserPhysicsList::CheckCutValues ";
>> 1234 G4cout << " Can not find key word " << keyword << G4endl;
>> 1235 }
>> 1236 #endif
>> 1237 }
930 1238
931 // ------------------------------------------- << 1239 // default cut value
932 G4bool G4VUserPhysicsList::RegisterProcess(G4V << 1240 G4double defaultCut;
933 { << 1241 if (ascii) {
934 return G4MT_thePLHelper->RegisterProcess(pro << 1242 fIn >> keyword >> defaultCut;
935 } << 1243 } else {
>> 1244 fIn.read( (char*)(&defaultCut), sizeof (G4double));
>> 1245 }
>> 1246 if (fIn.fail()) {
>> 1247 #ifdef G4VERBOSE
>> 1248 if (verboseLevel>0){
>> 1249 G4cout << "G4VUserPhysicsList::CheckCutValues ";
>> 1250 G4cout << " Bad data format " << G4endl;
>> 1251 }
>> 1252 #endif
>> 1253 fIn.close();
>> 1254 return false;
>> 1255 }
936 1256
937 // ------------------------------------------- << 1257 // check default value
938 G4ParticleTable::G4PTblDicIterator* G4VUserPhy << 1258 G4double ratio = abs(defaultCut / (defaultCutValue/mm) );
939 { << 1259 if ((keyword != "Default")|| (ratio<0.999) ||(ratio>1.001) ){
940 return (subInstanceManager.offset[g4vuplInst << 1260 #ifdef G4VERBOSE
941 } << 1261 if (verboseLevel>0){
>> 1262 G4cout << "G4VUserPhysicsList::CheckCutValues ";
>> 1263 G4cout << " Inconsistent default cut values" << G4endl;;
>> 1264 }
>> 1265 #endif
>> 1266 fIn.close();
>> 1267 return false;
>> 1268 }
942 1269
943 // ------------------------------------------- << 1270 // loop over all particles
944 void G4VUserPhysicsList::SetVerboseLevel(G4int << 1271 while(!fIn.eof()){
945 { << 1272 // read in particle name and cut in length
946 verboseLevel = value; << 1273 char name[FixedStringLengthForStore];
947 // set verboseLevel for G4ProductionCutsTabl << 1274 G4double cutLength;
948 // G4VUserPhysicsList: << 1275 if (ascii) {
949 fCutsTable->SetVerboseLevel(verboseLevel); << 1276 fIn >> name >> cutLength;
>> 1277 } else {
>> 1278 fIn.read(name, FixedStringLengthForStore);
>> 1279 fIn.read((char*)(&cutLength), sizeof (G4double));
>> 1280 }
>> 1281 if (fIn.eof()) break;
>> 1282 if (fIn.fail()) {
>> 1283 #ifdef G4VERBOSE
>> 1284 if (verboseLevel>0){
>> 1285 G4cout << "G4VUserPhysicsList::CheckCutValues ";
>> 1286 G4cout << " Bad data format " << G4endl;
>> 1287 }
>> 1288 #endif
>> 1289 fIn.close();
>> 1290 return false;
>> 1291 }
>> 1292
>> 1293 // Search particle
>> 1294 G4ParticleDefinition* particle = theParticleTable->FindParticle(name);
>> 1295 if (particle==0) {
>> 1296 #ifdef G4VERBOSE
>> 1297 if (verboseLevel>0){
>> 1298 G4cout << "G4VUserPhysicsList::CheckCutValues ";
>> 1299 G4cout << " Particle " << name <<" is not found "<< G4endl;
>> 1300 }
>> 1301 fIn.close();
>> 1302 return false;
>> 1303 #endif
>> 1304 }
>> 1305
>> 1306 // chech cut value in length
>> 1307 ratio = abs(cutLength/ (particle->GetLengthCuts()/mm) );
>> 1308 if ((ratio<0.999) ||(ratio>1.001) ){
>> 1309 #ifdef G4VERBOSE
>> 1310 if (verboseLevel>0){
>> 1311 G4cout << "G4VUserPhysicsList::CheckCutValues ";
>> 1312 G4cout << " Inconsistent cut values for " << name << G4endl;;
>> 1313 }
>> 1314 #endif
>> 1315 fIn.close();
>> 1316 return false;
>> 1317 }
>> 1318
>> 1319 G4double* cutArray = new G4double[numberOfMaterial];
>> 1320 // read in energy cut for all materials
>> 1321 if (ascii) {
>> 1322 /////////////// ASCII mode /////////////////
>> 1323 G4int jj;
>> 1324 for(size_t idx=0; idx<numberOfMaterial; ++idx) {
>> 1325 G4double value;
>> 1326 fIn >> value;
>> 1327 cutArray[idx] = value;
>> 1328 }
>> 1329 } else {
>> 1330 /////////////// Binary mode /////////////////
>> 1331 fIn.read( (char*)(cutArray), numberOfMaterial*(sizeof (G4double)) );
>> 1332 }
>> 1333
>> 1334 }
>> 1335
>> 1336 fIn.close();
>> 1337 return true;
>> 1338 }
>> 1339
>> 1340
>> 1341 ///////////////////////////////////////////////////////////////
>> 1342 G4bool G4VUserPhysicsList::RetrieveCutValues(const G4String& directory,
>> 1343 G4bool ascii )
>> 1344 {
>> 1345 if (!CheckMaterialInfo(directory, ascii)) {
>> 1346 #ifdef G4VERBOSE
>> 1347 G4cout << "G4VUserPhysicsList::RetrieveCutValues ";
>> 1348 G4cout << " Can not retrieve cut values " << G4endl;
>> 1349 #endif
>> 1350 }
>> 1351
>> 1352 // file name
>> 1353 const G4String fileName = directory + "/" + "cut_value.dat";
>> 1354
>> 1355 G4std::ifstream fIn;
>> 1356 // open input file //
>> 1357 #ifdef G4USE_STD_NAMESPACE
>> 1358 if (!ascii )
>> 1359 fIn.open(fileName,G4std::ios::in|G4std::ios::binary);
>> 1360 else
>> 1361 #endif
>> 1362 fIn.open(fileName,G4std::ios::in);
>> 1363
>> 1364 // check if the file has been opened successfully
>> 1365 if (!fIn) {
>> 1366 #ifdef G4VERBOSE
>> 1367 G4cerr << "G4VUserPhysicsList::RetrieveCutValues ";
>> 1368 G4cerr << " Can not open file " << fileName << G4endl;
>> 1369 #endif
>> 1370 return false;
>> 1371 }
950 1372
951 G4MT_thePLHelper->SetVerboseLevel(verboseLev << 1373 char temp[FixedStringLengthForStore];
>> 1374 size_t i;
952 1375
953 #ifdef G4VERBOSE << 1376 // key word
954 if (verboseLevel > 1) { << 1377 const G4String key = "CUT_VALUE";
955 G4cout << "G4VUserPhysicsList::SetVerboseL << 1378 G4String keyword;
956 << " Verbose level is set to " << v << 1379 if (ascii) {
>> 1380 fIn >> keyword;
>> 1381 } else {
>> 1382 fIn.read(temp, FixedStringLengthForStore);
>> 1383 keyword = (const char*)(temp);
>> 1384 }
>> 1385 if (key!=keyword) {
>> 1386 #ifdef G4VERBOSE
>> 1387 if (verboseLevel>0){
>> 1388 G4cout << "G4VUserPhysicsList::RetrieveCutValues ";
>> 1389 G4cout << " Can not find key word " << keyword << G4endl;
>> 1390 }
>> 1391 #endif
>> 1392 }
>> 1393
>> 1394 // default cut value
>> 1395 G4double defaultCut;
>> 1396 if (ascii) {
>> 1397 fIn >> keyword >> defaultCut;
>> 1398 } else {
>> 1399 fIn.read( (char*)(&defaultCut), sizeof (G4double));
>> 1400 }
>> 1401 if (fIn.fail()) {
>> 1402 #ifdef G4VERBOSE
>> 1403 if (verboseLevel>0){
>> 1404 G4cout << "G4VUserPhysicsList::RetrieveCutValues ";
>> 1405 G4cout << " Bad data format " << G4endl;
>> 1406 }
>> 1407 #endif
>> 1408 fIn.close();
>> 1409 return false;
>> 1410 }
>> 1411
>> 1412 // set default value
>> 1413 defaultCutValue = defaultCut;
>> 1414 #ifdef G4VERBOSE
>> 1415 if (verboseLevel>2){
>> 1416 G4cout.setf(G4std::ios::scientific);
>> 1417 G4cout << "Default " << G4std::setw(20) << defaultCutValue/mm << G4endl;
>> 1418 G4cout.unsetf(G4std::ios::scientific);
957 } 1419 }
958 #endif 1420 #endif
>> 1421
>> 1422 // Reset cut values for other particles
>> 1423 // loop over all particles in G4ParticleTable
>> 1424 theParticleIterator->reset();
>> 1425 while( (*theParticleIterator)() ){
>> 1426 G4ParticleDefinition* particle = theParticleIterator->value();
>> 1427 if (!particle->IsShortLived()) particle->ResetCuts();
>> 1428 }
>> 1429
>> 1430 // number of materials in the table
>> 1431 const G4MaterialTable* matTable = G4Material::GetMaterialTable();
>> 1432 numberOfMaterial = matTable->entries();
>> 1433
>> 1434 // loop over all particles
>> 1435 while(!fIn.eof()){
>> 1436 // read in particle name and cut in length
>> 1437 char name[FixedStringLengthForStore];
>> 1438 G4double cutLength;
>> 1439 if (ascii) {
>> 1440 fIn >> name >> cutLength;
>> 1441 } else {
>> 1442 fIn.read(name, FixedStringLengthForStore);
>> 1443 fIn.read((char*)(&cutLength), sizeof (G4double));
>> 1444 }
>> 1445
>> 1446 if (fIn.eof()) break;
>> 1447 if (fIn.fail()) {
>> 1448 #ifdef G4VERBOSE
>> 1449 if (verboseLevel>0){
>> 1450 G4cout << "G4VUserPhysicsList::RetrieveCutValues ";
>> 1451 G4cout << " Bad data format " << G4endl;
>> 1452 }
>> 1453 #endif
>> 1454 fIn.close();
>> 1455 return false;
>> 1456 }
>> 1457
>> 1458 #ifdef G4VERBOSE
>> 1459 if (verboseLevel>2){
>> 1460 G4cout.setf(G4std::ios::scientific);
>> 1461 G4cout << name << G4std::setw(20) << cutLength << G4endl;
>> 1462 G4cout.unsetf(G4std::ios::scientific);
>> 1463 }
>> 1464 #endif
>> 1465
>> 1466 // Search particle
>> 1467 G4ParticleDefinition* particle = theParticleTable->FindParticle(name);
>> 1468 if (particle==0) {
>> 1469 #ifdef G4VERBOSE
>> 1470 if (verboseLevel>0){
>> 1471 G4cout << "G4VUserPhysicsList::RetrieveCutValues ";
>> 1472 G4cout << " Particle " << name <<" is not found "<< G4endl;
>> 1473 }
>> 1474 fIn.close();
>> 1475 return false;
>> 1476 #endif
>> 1477 }
>> 1478
>> 1479 G4double* cutArray = new G4double[numberOfMaterial];
>> 1480
>> 1481 // read in energy cut for all materials
>> 1482 if (ascii) {
>> 1483 /////////////// ASCII mode /////////////////
>> 1484 G4int jj;
>> 1485 for(size_t idx=0; idx<numberOfMaterial; ++idx) {
>> 1486 G4double value;
>> 1487 fIn >> value;
>> 1488 cutArray[idx] = value;
>> 1489 }
>> 1490 } else {
>> 1491 /////////////// Binary mode /////////////////
>> 1492 fIn.read( (char*)(cutArray), numberOfMaterial*(sizeof (G4double)) );
>> 1493 }
>> 1494
>> 1495
>> 1496 // restore cut values
>> 1497 particle->RestoreCuts( cutLength , cutArray) ;
>> 1498 #ifdef G4VERBOSE
>> 1499 if (verboseLevel >2)DumpCutValues(particle);
>> 1500 #endif
>> 1501 }
>> 1502
>> 1503 #ifdef G4VERBOSE
>> 1504 if (verboseLevel >1){
>> 1505 G4cout << "G4VUserPhysicsList::RetrieveCutValues: Cut Values are successfully restored ";
>> 1506 G4cout << "CutLength : " << defaultCutValue/mm << " (mm)" << G4endl;
>> 1507 }
>> 1508 #endif
>> 1509
>> 1510 fIsRestoredCutValues = true;
>> 1511 fIn.close();
>> 1512 return true;
959 } 1513 }
>> 1514
>> 1515
>> 1516
>> 1517
>> 1518
>> 1519
>> 1520
>> 1521
960 1522