| 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 // G4Track class implementation <<
27 // 26 //
28 // Author: Katsuya Amako, KEK - 1996 << 27 // $Id: G4Track.cc 77080 2013-11-21 10:33:55Z gcosmo $
29 // Revisions: Hisaya Kurashige, 1998-2011 << 28 //
30 // ------------------------------------------- << 29 //
>> 30 //---------------------------------------------------------------
>> 31 //
>> 32 // G4Track.cc
>> 33 //
>> 34 //---------------------------------------------------------------
>> 35 // Add copy constructor Hisaya Feb. 07 01
>> 36 // Fix GetVelocity Hisaya Feb. 17 01
>> 37 // Modification for G4TouchableHandle 22 Oct. 2001 R.Chytracek//
>> 38 // Fix GetVelocity (bug report #741) Horton-Smith Apr 14 2005
>> 39 // Remove massless check in GetVelocity 02 Apr. 09 H.Kurashige
>> 40 // Use G4VelocityTable 17 AUg. 2011 H.Kurashige
31 41
32 #include "G4Track.hh" 42 #include "G4Track.hh"
33 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh"
34 #include "G4VAuxiliaryTrackInformation.hh" << 44 #include "G4ParticleTable.hh"
>> 45 #include "G4VelocityTable.hh"
35 46
36 #include <iostream> 47 #include <iostream>
37 #include <iomanip> 48 #include <iomanip>
38 49
39 // ------------------------------------------- << 50 G4ThreadLocal G4Allocator<G4Track> *aTrackAllocator = 0;
40 G4Allocator<G4Track>*& aTrackAllocator() << 51
41 { << 52 G4ThreadLocal G4VelocityTable* G4Track::velTable=0;
42 G4ThreadLocalStatic G4Allocator<G4Track>* _i <<
43 return _instance; <<
44 } <<
45 53
46 // ------------------------------------------- << 54 ///////////////////////////////////////////////////////////
47 G4Track::G4Track(G4DynamicParticle* apValueDyn 55 G4Track::G4Track(G4DynamicParticle* apValueDynamicParticle,
48 G4double aValueTime, 56 G4double aValueTime,
49 const G4ThreeVector& aValuePo 57 const G4ThreeVector& aValuePosition)
50 : fPosition(aValuePosition) << 58 ///////////////////////////////////////////////////////////
51 , fGlobalTime(aValueTime) << 59 : fCurrentStepNumber(0), fPosition(aValuePosition),
52 , fVelocity(c_light) << 60 fGlobalTime(aValueTime), fLocalTime(0.),
53 { << 61 fTrackLength(0.),
54 fpDynamicParticle = (apValueDynamicParticle) << 62 fParentID(0), fTrackID(0),
55 ? apValueDynamicParticle : << 63 fVelocity(c_light),
>> 64 fpDynamicParticle(apValueDynamicParticle),
>> 65 fTrackStatus(fAlive),
>> 66 fBelowThreshold(false), fGoodForTracking(false),
>> 67 fStepLength(0.0), fWeight(1.0),
>> 68 fpStep(0),
>> 69 fVtxKineticEnergy(0.0),
>> 70 fpLVAtVertex(0), fpCreatorProcess(0),
>> 71 fCreatorModelIndex(-1),
>> 72 fpUserInformation(0),
>> 73 prev_mat(0), groupvel(0),
>> 74 prev_velocity(0.0), prev_momentum(0.0),
>> 75 is_OpticalPhoton(false),
>> 76 useGivenVelocity(false)
>> 77 {
>> 78 static G4ThreadLocal G4bool isFirstTime = true;
>> 79 static G4ThreadLocal G4ParticleDefinition* fOpticalPhoton =0;
>> 80 if ( isFirstTime ) {
>> 81 isFirstTime = false;
>> 82 // set fOpticalPhoton
>> 83 fOpticalPhoton = G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
>> 84 }
56 // check if the particle type is Optical Pho 85 // check if the particle type is Optical Photon
57 is_OpticalPhoton = << 86 is_OpticalPhoton = (fpDynamicParticle->GetDefinition() == fOpticalPhoton);
58 (fpDynamicParticle->GetDefinition()->GetPD << 87
>> 88 if (velTable ==0 ) velTable = G4VelocityTable::GetVelocityTable();
>> 89
>> 90 fVelocity = CalculateVelocity();
>> 91
59 } 92 }
60 93
61 // ------------------------------------------- << 94 //////////////////
62 G4Track::G4Track() 95 G4Track::G4Track()
63 : fVelocity(c_light) << 96 //////////////////
64 , fpDynamicParticle(new G4DynamicParticle()) << 97 : fCurrentStepNumber(0),
65 {} << 98 fGlobalTime(0), fLocalTime(0.),
66 << 99 fTrackLength(0.),
67 // ------------------------------------------- << 100 fParentID(0), fTrackID(0),
68 G4Track::G4Track(const G4Track& right,G4bool c << 101 fVelocity(c_light),
69 : fVelocity(c_light) << 102 fpDynamicParticle(0),
>> 103 fTrackStatus(fAlive),
>> 104 fBelowThreshold(false), fGoodForTracking(false),
>> 105 fStepLength(0.0), fWeight(1.0),
>> 106 fpStep(0),
>> 107 fVtxKineticEnergy(0.0),
>> 108 fpLVAtVertex(0), fpCreatorProcess(0),
>> 109 fCreatorModelIndex(-1),
>> 110 fpUserInformation(0),
>> 111 prev_mat(0), groupvel(0),
>> 112 prev_velocity(0.0), prev_momentum(0.0),
>> 113 is_OpticalPhoton(false),
>> 114 useGivenVelocity(false)
>> 115 {
>> 116 }
>> 117
>> 118 //////////////////
>> 119 G4Track::G4Track(const G4Track& right)
>> 120 //////////////////
>> 121 : fCurrentStepNumber(0),
>> 122 fGlobalTime(0), fLocalTime(0.),
>> 123 fTrackLength(0.),
>> 124 fParentID(0), fTrackID(0),
>> 125 fVelocity(c_light),
>> 126 fpDynamicParticle(0),
>> 127 fTrackStatus(fAlive),
>> 128 fBelowThreshold(false), fGoodForTracking(false),
>> 129 fStepLength(0.0), fWeight(1.0),
>> 130 fpStep(0),
>> 131 fVtxKineticEnergy(0.0),
>> 132 fpLVAtVertex(0), fpCreatorProcess(0),
>> 133 fCreatorModelIndex(-1),
>> 134 fpUserInformation(0),
>> 135 prev_mat(0), groupvel(0),
>> 136 prev_velocity(0.0), prev_momentum(0.0),
>> 137 is_OpticalPhoton(false),
>> 138 useGivenVelocity(false)
70 { 139 {
71 fCopyTouchables = copyTouchables; <<
72 *this = right; 140 *this = right;
73 fCopyTouchables = true; <<
74 } 141 }
75 142
76 // ------------------------------------------- << 143 ///////////////////
77 G4Track::~G4Track() 144 G4Track::~G4Track()
>> 145 ///////////////////
78 { 146 {
79 delete fpDynamicParticle; << 147 delete fpDynamicParticle;
80 delete fpUserInformation; << 148 delete fpUserInformation;
81 ClearAuxiliaryTrackInformation(); << 149 }
82 } <<
83 <<
84 // ------------------------------------------- <<
85 G4Track& G4Track::operator=(const G4Track& rig <<
86 { <<
87 if(this != &right) <<
88 { <<
89 fPosition = right.fPosition; <<
90 fGlobalTime = right.fGlobalTime; <<
91 fLocalTime = right.fLocalTime; <<
92 fTrackLength = right.fTrackLength; <<
93 fWeight = right.fWeight; <<
94 fStepLength = right.fStepLength; <<
95 <<
96 // additional fields required for geometri <<
97 if(fCopyTouchables) { <<
98 fpTouchable = right.fpTouchable; <<
99 fpNextTouchable = right.fpNextTouchable; <<
100 fpOriginTouchable = right.fpOriginToucha <<
101 } <<
102 150
103 // Track ID (and Parent ID) is not copied << 151 //////////////////
104 fTrackID = 0; << 152 G4Track & G4Track::operator=(const G4Track &right)
105 fParentID = 0; << 153 //////////////////
106 << 154 {
107 // CurrentStepNumber is set to be 0 << 155 if (this != &right) {
108 fCurrentStepNumber = 0; << 156 fPosition = right.fPosition;
109 << 157 fGlobalTime = right.fGlobalTime;
110 // Creator model ID << 158 fLocalTime = right.fLocalTime;
111 fCreatorModelID = right.fCreatorModelID; << 159 fTrackLength = right.fTrackLength;
112 << 160 fWeight = right.fWeight;
113 // Parent resonance << 161 fStepLength = right.fStepLength;
114 fParentResonanceDef = right.fParentResonan << 162
115 fParentResonanceID = right.fParentResonan << 163 // Track ID (and Parent ID) is not copied and set to zero for new track
>> 164 fTrackID = 0;
>> 165 fParentID =0;
>> 166
>> 167 // CurrentStepNumber is set to be 0
>> 168 fCurrentStepNumber = 0;
>> 169
>> 170 // velocity information
>> 171 fVelocity = right.fVelocity;
>> 172
>> 173 // dynamic particle information
>> 174 fpDynamicParticle = new G4DynamicParticle(*(right.fpDynamicParticle));
>> 175
>> 176 // track status and flags for tracking
>> 177 fTrackStatus = right.fTrackStatus;
>> 178 fBelowThreshold = right.fBelowThreshold;
>> 179 fGoodForTracking = right.fGoodForTracking;
116 180
117 // velocity information << 181 // Step information (Step Length, Step Number, pointer to the Step,)
118 fVelocity = right.fVelocity; << 182 // are not copied
119 << 183 fpStep=0;
120 // dynamic particle information << 184
121 delete fpDynamicParticle; << 185 // vertex information
122 fpDynamicParticle = new G4DynamicParticle( << 186 fVtxPosition = right.fVtxPosition;
123 << 187 fpLVAtVertex = right.fpLVAtVertex;
124 // track status and flags for tracking << 188 fVtxKineticEnergy = right.fVtxKineticEnergy;
125 fTrackStatus = right.fTrackStatus; << 189 fVtxMomentumDirection = right.fVtxMomentumDirection;
126 fBelowThreshold = right.fBelowThreshold; << 190
127 fGoodForTracking = right.fGoodForTracking; << 191 // CreatorProcess and UserInformation are not copied
128 << 192 fpCreatorProcess = 0;
129 // Step information (Step Length, Step Num << 193 fpUserInformation = 0;
130 // are not copied << 194
131 fpStep = nullptr; << 195 prev_mat = right.prev_mat;
132 << 196 groupvel = right.groupvel;
133 // vertex information << 197 prev_velocity = right.prev_velocity;
134 fVtxPosition = right.fVtxPosition << 198 prev_momentum = right.prev_momentum;
135 fpLVAtVertex = right.fpLVAtVertex <<
136 fVtxKineticEnergy = right.fVtxKineticE <<
137 fVtxMomentumDirection = right.fVtxMomentum <<
138 <<
139 // CreatorProcess and UserInformation are <<
140 fpCreatorProcess = nullptr; <<
141 delete fpUserInformation; <<
142 fpUserInformation = nullptr; <<
143 <<
144 prev_mat = right.prev_mat; <<
145 groupvel = right.groupvel; <<
146 prev_velocity = right.prev_velocity; <<
147 prev_momentum = right.prev_momentum; <<
148 <<
149 is_OpticalPhoton = right.is_OpticalPhoton; <<
150 useGivenVelocity = right.useGivenVelocity; <<
151 199
152 ClearAuxiliaryTrackInformation(); << 200 is_OpticalPhoton = right.is_OpticalPhoton;
>> 201 useGivenVelocity = right.useGivenVelocity;
153 } 202 }
154 return *this; 203 return *this;
155 } 204 }
156 205
157 // ------------------------------------------- << 206 ///////////////////
158 void G4Track::CopyTrackInfo(const G4Track& rig << 207 void G4Track::CopyTrackInfo(const G4Track& right)
>> 208 //////////////////
159 { 209 {
160 fCopyTouchables = copyTouchables; <<
161 *this = right; 210 *this = right;
162 fCopyTouchables = true; <<
163 } 211 }
164 212
165 // ------------------------------------------- << 213 ///////////////////
>> 214 G4double G4Track::CalculateVelocity() const
>> 215 ///////////////////
>> 216 {
>> 217 if (useGivenVelocity) return fVelocity;
>> 218
>> 219 G4double velocity = c_light ;
>> 220
>> 221 G4double mass = fpDynamicParticle->GetMass();
>> 222
>> 223 // special case for photons
>> 224 if ( is_OpticalPhoton ) return CalculateVelocityForOpticalPhoton();
>> 225
>> 226 // particles other than optical photon
>> 227 if (mass<DBL_MIN) {
>> 228 // Zero Mass
>> 229 velocity = c_light;
>> 230 } else {
>> 231 G4double T = (fpDynamicParticle->GetKineticEnergy())/mass;
>> 232 if (T > GetMaxTOfVelocityTable()) {
>> 233 velocity = c_light;
>> 234 } else if (T<DBL_MIN) {
>> 235 velocity =0.;
>> 236 } else if (T<GetMinTOfVelocityTable()) {
>> 237 velocity = c_light*std::sqrt(T*(T+2.))/(T+1.0);
>> 238 } else {
>> 239 velocity = velTable->Value(T);
>> 240 }
>> 241
>> 242 }
>> 243 return velocity ;
>> 244 }
>> 245
>> 246 ///////////////////
166 G4double G4Track::CalculateVelocityForOpticalP 247 G4double G4Track::CalculateVelocityForOpticalPhoton() const
>> 248 ///////////////////
167 { 249 {
168 G4double velocity = c_light; << 250
>> 251 G4double velocity = c_light ;
>> 252
169 253
170 G4Material* mat = nullptr; << 254 G4Material* mat=0;
171 G4bool update_groupvel = false; 255 G4bool update_groupvel = false;
172 if(fpStep != nullptr) << 256 if ( fpStep !=0 ){
173 { << 257 mat= this->GetMaterial(); // Fix for repeated volumes
174 mat = this->GetMaterial(); // Fix for r << 258 }else{
175 } << 259 if (fpTouchable!=0){
176 else << 260 mat=fpTouchable->GetVolume()->GetLogicalVolume()->GetMaterial();
177 { <<
178 if(fpTouchable) <<
179 { <<
180 mat = fpTouchable->GetVolume()->GetLogic <<
181 } 261 }
182 } 262 }
183 // check if previous step is in the same vol 263 // check if previous step is in the same volume
184 // and get new GROUPVELOCITY table if neces << 264 // and get new GROUPVELOCITY table if necessary
185 if((mat != nullptr) && ((mat != prev_mat) || << 265 if ((mat != 0) && ((mat != prev_mat)||(groupvel==0))) {
186 { << 266 groupvel = 0;
187 groupvel = nullptr; << 267 if(mat->GetMaterialPropertiesTable() != 0)
188 if(mat->GetMaterialPropertiesTable() != nu << 268 groupvel = mat->GetMaterialPropertiesTable()->GetProperty("GROUPVEL");
189 groupvel = mat->GetMaterialPropertiesTab <<
190 update_groupvel = true; 269 update_groupvel = true;
191 } 270 }
192 prev_mat = mat; 271 prev_mat = mat;
193 << 272
194 if(groupvel != nullptr) << 273 if (groupvel != 0 ) {
195 { <<
196 // light velocity = c/(rindex+d(rindex)/d( 274 // light velocity = c/(rindex+d(rindex)/d(log(E_phot)))
197 // values stored in GROUPVEL material prop 275 // values stored in GROUPVEL material properties vector
198 velocity = prev_velocity; << 276 velocity = prev_velocity;
199 << 277
200 // check if momentum is same as in the pre 278 // check if momentum is same as in the previous step
201 // and calculate group velocity if necess << 279 // and calculate group velocity if necessary
202 G4double current_momentum = fpDynamicParti 280 G4double current_momentum = fpDynamicParticle->GetTotalMomentum();
203 if(update_groupvel || (current_momentum != << 281 if( update_groupvel || (current_momentum != prev_momentum) ) {
204 { << 282 velocity =
205 velocity = groupvel->Value(current_ << 283 groupvel->Value(current_momentum);
206 prev_velocity = velocity; 284 prev_velocity = velocity;
207 prev_momentum = current_momentum; 285 prev_momentum = current_momentum;
208 } 286 }
209 } << 287 }
210 << 288
211 return velocity; << 289 return velocity ;
212 } <<
213 <<
214 // ------------------------------------------- <<
215 void G4Track::SetAuxiliaryTrackInformation(G4i <<
216 G4VAuxiliaryTrackInformation* in <<
217 { <<
218 if(fpAuxiliaryTrackInformationMap == nullptr <<
219 { <<
220 fpAuxiliaryTrackInformationMap = <<
221 new std::map<G4int, G4VAuxiliaryTrackInf <<
222 } <<
223 if(G4PhysicsModelCatalog::GetModelIndex(id) <<
224 { <<
225 G4ExceptionDescription ED; <<
226 ED << "Process/model ID <" << id << "> is <<
227 G4Exception("G4VAuxiliaryTrackInformation: <<
228 "TRACK0982", FatalException, E <<
229 } <<
230 (*fpAuxiliaryTrackInformationMap)[id] = info <<
231 } 290 }
232 291
233 // ------------------------------------------- << 292 ///////////////////
234 G4VAuxiliaryTrackInformation* << 293 void G4Track::SetVelocityTableProperties(G4double t_max, G4double t_min, G4int nbin)
235 G4Track::GetAuxiliaryTrackInformation(G4int id << 294 ///////////////////
236 { 295 {
237 if(fpAuxiliaryTrackInformationMap == nullptr << 296 G4VelocityTable::SetVelocityTableProperties(t_max, t_min, nbin);
238 return nullptr; << 297 velTable = G4VelocityTable::GetVelocityTable();
239 <<
240 auto itr = fpAuxiliaryTrackInformationMap->f <<
241 if(itr == fpAuxiliaryTrackInformationMap->ce <<
242 return nullptr; <<
243 return (*itr).second; <<
244 } 298 }
245 299
246 // ------------------------------------------- << 300 ///////////////////
247 void G4Track::RemoveAuxiliaryTrackInformation( << 301 G4double G4Track::GetMaxTOfVelocityTable()
248 { << 302 ///////////////////
249 if(fpAuxiliaryTrackInformationMap != nullptr << 303 { return G4VelocityTable::GetMaxTOfVelocityTable(); }
250 fpAuxiliaryTrackInformationMap->find(id) <<
251 { <<
252 fpAuxiliaryTrackInformationMap->erase(id); <<
253 } <<
254 } <<
255 304
256 // ------------------------------------------- << 305 ///////////////////
257 void G4Track::RemoveAuxiliaryTrackInformation( << 306 G4double G4Track::GetMinTOfVelocityTable()
258 { << 307 ///////////////////
259 if(fpAuxiliaryTrackInformationMap != nullptr << 308 { return G4VelocityTable::GetMinTOfVelocityTable(); }
260 { <<
261 G4int id = G4PhysicsModelCatalog::GetModel <<
262 RemoveAuxiliaryTrackInformation(id); <<
263 } <<
264 } <<
265 309
266 // ------------------------------------------- << 310 ///////////////////
267 void G4Track::ClearAuxiliaryTrackInformation() << 311 G4int G4Track::GetNbinOfVelocityTable()
268 { << 312 ///////////////////
269 if(fpAuxiliaryTrackInformationMap == nullptr << 313 { return G4VelocityTable::GetNbinOfVelocityTable(); }
270 return; <<
271 for(const auto& itr : *fpAuxiliaryTrackInfor <<
272 { <<
273 delete itr.second; <<
274 } <<
275 delete fpAuxiliaryTrackInformationMap; <<
276 fpAuxiliaryTrackInformationMap = nullptr; <<
277 } <<
278 314