| 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 // G4VelocityTable class implementation 26 // G4VelocityTable class implementation
27 // 27 //
28 // Author: Hisaya Kurashige, 17 August 2011 28 // Author: Hisaya Kurashige, 17 August 2011
29 // ------------------------------------------- 29 // --------------------------------------------------------------------
30 30
31 #include "G4VelocityTable.hh" 31 #include "G4VelocityTable.hh"
32 #include "G4PhysicalConstants.hh" 32 #include "G4PhysicalConstants.hh"
33 #include "G4StateManager.hh" 33 #include "G4StateManager.hh"
34 #include "G4ApplicationState.hh" 34 #include "G4ApplicationState.hh"
35 #include "G4Log.hh" 35 #include "G4Log.hh"
36 #include "G4Exp.hh" 36 #include "G4Exp.hh"
37 37
38 #include "G4ios.hh" 38 #include "G4ios.hh"
39 39
40 G4ThreadLocal G4VelocityTable* G4VelocityTable 40 G4ThreadLocal G4VelocityTable* G4VelocityTable::theInstance = nullptr;
41 41
42 // ------------------------------------------- 42 // --------------------------------------------------------------------
43 G4VelocityTable::G4VelocityTable() 43 G4VelocityTable::G4VelocityTable()
44 { 44 {
45 PrepareVelocityTable(); 45 PrepareVelocityTable();
46 } 46 }
47 47
48 // ------------------------------------------- 48 // --------------------------------------------------------------------
49 G4VelocityTable::~G4VelocityTable() 49 G4VelocityTable::~G4VelocityTable()
50 { 50 {
51 dataVector.clear(); 51 dataVector.clear();
52 binVector.clear(); 52 binVector.clear();
53 } 53 }
54 54
55 // ------------------------------------------- 55 // --------------------------------------------------------------------
56 void G4VelocityTable::PrepareVelocityTable() 56 void G4VelocityTable::PrepareVelocityTable()
57 { 57 {
58 // const G4double g4log10 = std::log(10.); 58 // const G4double g4log10 = std::log(10.);
59 59
60 dataVector.clear(); 60 dataVector.clear();
61 binVector.clear(); 61 binVector.clear();
62 dBin = G4Log(maxT / minT) / NbinT; 62 dBin = G4Log(maxT / minT) / NbinT;
63 baseBin = G4Log(minT) / dBin; 63 baseBin = G4Log(minT) / dBin;
64 64
65 numberOfNodes = NbinT + 1; 65 numberOfNodes = NbinT + 1;
66 dataVector.reserve(numberOfNodes); 66 dataVector.reserve(numberOfNodes);
67 binVector.reserve(numberOfNodes); 67 binVector.reserve(numberOfNodes);
68 68
69 binVector.push_back(minT); 69 binVector.push_back(minT);
70 dataVector.push_back(0.0); 70 dataVector.push_back(0.0);
71 71
72 for(std::size_t i = 1; i < numberOfNodes - 1 72 for(std::size_t i = 1; i < numberOfNodes - 1; ++i)
73 { 73 {
74 binVector.push_back(G4Exp((baseBin + i) * 74 binVector.push_back(G4Exp((baseBin + i) * dBin));
75 dataVector.push_back(0.0); 75 dataVector.push_back(0.0);
76 } 76 }
77 binVector.push_back(maxT); 77 binVector.push_back(maxT);
78 dataVector.push_back(0.0); 78 dataVector.push_back(0.0);
79 79
80 edgeMin = binVector[0]; 80 edgeMin = binVector[0];
81 edgeMax = binVector[numberOfNodes - 1]; 81 edgeMax = binVector[numberOfNodes - 1];
82 82
83 for(G4int i = 0; i <= NbinT; ++i) 83 for(G4int i = 0; i <= NbinT; ++i)
84 { 84 {
85 G4double T = binVector[i]; 85 G4double T = binVector[i];
86 dataVector[i] = c_light * std::sqrt(T * (T 86 dataVector[i] = c_light * std::sqrt(T * (T + 2.)) / (T + 1.0);
87 } 87 }
88 88
89 return; 89 return;
90 } 90 }
91 91
92 // ------------------------------------------- 92 // --------------------------------------------------------------------
93 std::size_t G4VelocityTable::FindBinLocation(G 93 std::size_t G4VelocityTable::FindBinLocation(G4double theEnergy) const
94 { 94 {
95 // For G4PhysicsLogVector, FindBinLocation i 95 // For G4PhysicsLogVector, FindBinLocation is implemented using
96 // a simple arithmetic calculation. 96 // a simple arithmetic calculation.
97 // 97 //
98 // Because this is a virtual function, it is 98 // Because this is a virtual function, it is accessed through a
99 // pointer to the G4PhysicsVector object for 99 // pointer to the G4PhysicsVector object for most usages. In this
100 // case, 'inline' will not be invoked. Howev 100 // case, 'inline' will not be invoked. However, there is a possibility
101 // that the user access to the G4PhysicsLogV 101 // that the user access to the G4PhysicsLogVector object directly and
102 // not through pointers or references. In th 102 // not through pointers or references. In this case, the 'inline' will
103 // be invoked. (See R.B.Murray, "C++ Strateg 103 // be invoked. (See R.B.Murray, "C++ Strategies and Tactics", Chap.6.6)
104 104
105 // const G4double g4log10 = G4Log(10.); 105 // const G4double g4log10 = G4Log(10.);
106 return std::size_t(G4Log(theEnergy) / dBin - 106 return std::size_t(G4Log(theEnergy) / dBin - baseBin);
107 } 107 }
108 108
109 // ------------------------------------------- 109 // --------------------------------------------------------------------
110 G4double G4VelocityTable::Value(G4double theEn 110 G4double G4VelocityTable::Value(G4double theEnergy)
111 { 111 {
112 // Use cache for speed up - check if the val 112 // Use cache for speed up - check if the value 'theEnergy' is same as the
113 // last call. If it is same, then use the la 113 // last call. If it is same, then use the last bin location. Also the
114 // value 'theEnergy' lies between the last e 114 // value 'theEnergy' lies between the last energy and low edge of of the
115 // bin of last call, then the last bin locat 115 // bin of last call, then the last bin location is used
116 116
117 if(theEnergy == lastEnergy) 117 if(theEnergy == lastEnergy)
118 { 118 {
119 } 119 }
120 else if(theEnergy < lastEnergy && theEnergy 120 else if(theEnergy < lastEnergy && theEnergy >= binVector[lastBin])
121 { 121 {
122 lastEnergy = theEnergy; 122 lastEnergy = theEnergy;
123 lastValue = Interpolation(); 123 lastValue = Interpolation();
124 } 124 }
125 else if(theEnergy <= edgeMin) 125 else if(theEnergy <= edgeMin)
126 { 126 {
127 lastBin = 0; 127 lastBin = 0;
128 lastEnergy = edgeMin; 128 lastEnergy = edgeMin;
129 lastValue = dataVector[0]; 129 lastValue = dataVector[0];
130 } 130 }
131 else if(theEnergy >= edgeMax) 131 else if(theEnergy >= edgeMax)
132 { 132 {
133 lastBin = numberOfNodes - 1; 133 lastBin = numberOfNodes - 1;
134 lastEnergy = edgeMax; 134 lastEnergy = edgeMax;
135 lastValue = dataVector[lastBin]; 135 lastValue = dataVector[lastBin];
136 } 136 }
137 else 137 else
138 { 138 {
139 lastBin = (std::size_t)(G4Log(theEnergy) / 139 lastBin = (std::size_t)(G4Log(theEnergy) / dBin - baseBin);
140 if(lastBin == numberOfNodes) 140 if(lastBin == numberOfNodes)
141 { 141 {
142 --lastBin; 142 --lastBin;
143 } // VI: fix possible precision lost 143 } // VI: fix possible precision lost
144 lastEnergy = theEnergy; 144 lastEnergy = theEnergy;
145 lastValue = Interpolation(); 145 lastValue = Interpolation();
146 } 146 }
147 return lastValue; 147 return lastValue;
148 } 148 }
149 149
150 // ------------------------------------------- 150 // --------------------------------------------------------------------
151 G4VelocityTable* G4VelocityTable::GetVelocityT 151 G4VelocityTable* G4VelocityTable::GetVelocityTable()
152 { 152 {
153 if(theInstance == nullptr) 153 if(theInstance == nullptr)
154 { 154 {
155 static G4ThreadLocalSingleton<G4VelocityTa 155 static G4ThreadLocalSingleton<G4VelocityTable> inst;
156 theInstance = inst.Instance(); 156 theInstance = inst.Instance();
157 } 157 }
158 return theInstance; 158 return theInstance;
159 } 159 }
160 160
161 // ------------------------------------------- 161 // --------------------------------------------------------------------
162 void G4VelocityTable:: 162 void G4VelocityTable::
163 SetVelocityTableProperties(G4double t_max, G4d 163 SetVelocityTableProperties(G4double t_max, G4double t_min, G4int nbin)
164 { 164 {
165 if(theInstance == nullptr) 165 if(theInstance == nullptr)
166 { 166 {
167 GetVelocityTable(); 167 GetVelocityTable();
168 } 168 }
169 169
170 G4StateManager* stateManager = G4StateMan 170 G4StateManager* stateManager = G4StateManager::GetStateManager();
171 G4ApplicationState currentState = stateManag 171 G4ApplicationState currentState = stateManager->GetCurrentState();
172 172
173 // check if state is outside event loop 173 // check if state is outside event loop
174 if(!(currentState == G4State_Idle || current 174 if(!(currentState == G4State_Idle || currentState == G4State_PreInit))
175 { 175 {
176 G4Exception("G4VelocityTable::SetVelocityT 176 G4Exception("G4VelocityTable::SetVelocityTableProperties()", "Track101",
177 JustWarning, 177 JustWarning,
178 "Can modify only in PreInit or 178 "Can modify only in PreInit or Idle state : Method ignored.");
179 return; 179 return;
180 } 180 }
181 181
182 if(nbin > 100) 182 if(nbin > 100)
183 theInstance->NbinT = nbin; 183 theInstance->NbinT = nbin;
184 if((t_min < t_max) && (t_min > 0.)) 184 if((t_min < t_max) && (t_min > 0.))
185 { 185 {
186 theInstance->minT = t_min; 186 theInstance->minT = t_min;
187 theInstance->maxT = t_max; 187 theInstance->maxT = t_max;
188 } 188 }
189 theInstance->PrepareVelocityTable(); 189 theInstance->PrepareVelocityTable();
190 } 190 }
191 191
192 // ------------------------------------------- 192 // --------------------------------------------------------------------
193 G4double G4VelocityTable::GetMaxTOfVelocityTab 193 G4double G4VelocityTable::GetMaxTOfVelocityTable()
194 { 194 {
195 return GetVelocityTable()->maxT; 195 return GetVelocityTable()->maxT;
196 } 196 }
197 197
198 // ------------------------------------------- 198 // --------------------------------------------------------------------
199 G4double G4VelocityTable::GetMinTOfVelocityTab 199 G4double G4VelocityTable::GetMinTOfVelocityTable()
200 { 200 {
201 return GetVelocityTable()->minT; 201 return GetVelocityTable()->minT;
202 } 202 }
203 203
204 // ------------------------------------------- 204 // --------------------------------------------------------------------
205 G4int G4VelocityTable::GetNbinOfVelocityTable( 205 G4int G4VelocityTable::GetNbinOfVelocityTable()
206 { 206 {
207 return GetVelocityTable()->NbinT; 207 return GetVelocityTable()->NbinT;
208 } 208 }
209 209