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 91231 2015-06-26 10:40:45Z 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" >> 44 #include "G4ParticleTable.hh" >> 45 #include "G4VelocityTable.hh" 34 #include "G4VAuxiliaryTrackInformation.hh" 46 #include "G4VAuxiliaryTrackInformation.hh" >> 47 #include "G4PhysicsModelCatalog.hh" 35 48 36 #include <iostream> 49 #include <iostream> 37 #include <iomanip> 50 #include <iomanip> 38 51 39 // ------------------------------------------- << 52 G4ThreadLocal G4Allocator<G4Track> *aTrackAllocator = 0; 40 G4Allocator<G4Track>*& aTrackAllocator() << 53 41 { << 54 G4ThreadLocal G4VelocityTable* G4Track::velTable=0; 42 G4ThreadLocalStatic G4Allocator<G4Track>* _i << 43 return _instance; << 44 } << 45 55 46 // ------------------------------------------- << 56 /////////////////////////////////////////////////////////// 47 G4Track::G4Track(G4DynamicParticle* apValueDyn 57 G4Track::G4Track(G4DynamicParticle* apValueDynamicParticle, 48 G4double aValueTime, 58 G4double aValueTime, 49 const G4ThreeVector& aValuePo 59 const G4ThreeVector& aValuePosition) 50 : fPosition(aValuePosition) << 60 /////////////////////////////////////////////////////////// 51 , fGlobalTime(aValueTime) << 61 : fCurrentStepNumber(0), fPosition(aValuePosition), 52 , fVelocity(c_light) << 62 fGlobalTime(aValueTime), fLocalTime(0.), 53 { << 63 fTrackLength(0.), 54 fpDynamicParticle = (apValueDynamicParticle) << 64 fParentID(0), fTrackID(0), 55 ? apValueDynamicParticle : << 65 fVelocity(c_light), >> 66 fpDynamicParticle(apValueDynamicParticle), >> 67 fTrackStatus(fAlive), >> 68 fBelowThreshold(false), fGoodForTracking(false), >> 69 fStepLength(0.0), fWeight(1.0), >> 70 fpStep(0), >> 71 fVtxKineticEnergy(0.0), >> 72 fpLVAtVertex(0), fpCreatorProcess(0), >> 73 fCreatorModelIndex(-1), >> 74 fpUserInformation(0), >> 75 prev_mat(0), groupvel(0), >> 76 prev_velocity(0.0), prev_momentum(0.0), >> 77 is_OpticalPhoton(false), >> 78 useGivenVelocity(false), >> 79 fpAuxiliaryTrackInformationMap(0) >> 80 { >> 81 static G4ThreadLocal G4bool isFirstTime = true; >> 82 static G4ThreadLocal G4ParticleDefinition* fOpticalPhoton =0; >> 83 if ( isFirstTime ) { >> 84 isFirstTime = false; >> 85 // set fOpticalPhoton >> 86 fOpticalPhoton = G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton"); >> 87 } 56 // check if the particle type is Optical Pho 88 // check if the particle type is Optical Photon 57 is_OpticalPhoton = << 89 is_OpticalPhoton = (fpDynamicParticle->GetDefinition() == fOpticalPhoton); 58 (fpDynamicParticle->GetDefinition()->GetPD << 90 >> 91 if (velTable ==0 ) velTable = G4VelocityTable::GetVelocityTable(); >> 92 >> 93 fVelocity = CalculateVelocity(); >> 94 59 } 95 } 60 96 61 // ------------------------------------------- << 97 ////////////////// 62 G4Track::G4Track() 98 G4Track::G4Track() 63 : fVelocity(c_light) << 99 ////////////////// 64 , fpDynamicParticle(new G4DynamicParticle()) << 100 : fCurrentStepNumber(0), 65 {} << 101 fGlobalTime(0), fLocalTime(0.), 66 << 102 fTrackLength(0.), 67 // ------------------------------------------- << 103 fParentID(0), fTrackID(0), 68 G4Track::G4Track(const G4Track& right,G4bool c << 104 fVelocity(c_light), 69 : fVelocity(c_light) << 105 fpDynamicParticle(0), >> 106 fTrackStatus(fAlive), >> 107 fBelowThreshold(false), fGoodForTracking(false), >> 108 fStepLength(0.0), fWeight(1.0), >> 109 fpStep(0), >> 110 fVtxKineticEnergy(0.0), >> 111 fpLVAtVertex(0), fpCreatorProcess(0), >> 112 fCreatorModelIndex(-1), >> 113 fpUserInformation(0), >> 114 prev_mat(0), groupvel(0), >> 115 prev_velocity(0.0), prev_momentum(0.0), >> 116 is_OpticalPhoton(false), >> 117 useGivenVelocity(false), >> 118 fpAuxiliaryTrackInformationMap(0) 70 { 119 { 71 fCopyTouchables = copyTouchables; << 72 *this = right; << 73 fCopyTouchables = true; << 74 } 120 } 75 121 76 // ------------------------------------------- << 122 ////////////////// 77 G4Track::~G4Track() << 123 G4Track::G4Track(const G4Track& right) >> 124 ////////////////// >> 125 : fCurrentStepNumber(0), >> 126 fGlobalTime(0), fLocalTime(0.), >> 127 fTrackLength(0.), >> 128 fParentID(0), fTrackID(0), >> 129 fVelocity(c_light), >> 130 fpDynamicParticle(0), >> 131 fTrackStatus(fAlive), >> 132 fBelowThreshold(false), fGoodForTracking(false), >> 133 fStepLength(0.0), fWeight(1.0), >> 134 fpStep(0), >> 135 fVtxKineticEnergy(0.0), >> 136 fpLVAtVertex(0), fpCreatorProcess(0), >> 137 fCreatorModelIndex(-1), >> 138 fpUserInformation(0), >> 139 prev_mat(0), groupvel(0), >> 140 prev_velocity(0.0), prev_momentum(0.0), >> 141 is_OpticalPhoton(false), >> 142 useGivenVelocity(false), >> 143 fpAuxiliaryTrackInformationMap(0) 78 { 144 { 79 delete fpDynamicParticle; << 145 *this = right; 80 delete fpUserInformation; << 81 ClearAuxiliaryTrackInformation(); << 82 } 146 } 83 147 84 // ------------------------------------------- << 148 /////////////////// 85 G4Track& G4Track::operator=(const G4Track& rig << 149 G4Track::~G4Track() >> 150 /////////////////// 86 { 151 { 87 if(this != &right) << 152 delete fpDynamicParticle; 88 { << 153 delete fpUserInformation; 89 fPosition = right.fPosition; << 154 ClearAuxiliaryTrackInformation(); 90 fGlobalTime = right.fGlobalTime; << 155 } 91 fLocalTime = right.fLocalTime; << 156 92 fTrackLength = right.fTrackLength; << 157 ////////////////// 93 fWeight = right.fWeight; << 158 G4Track & G4Track::operator=(const G4Track &right) 94 fStepLength = right.fStepLength; << 159 ////////////////// 95 << 160 { 96 // additional fields required for geometri << 161 if (this != &right) { 97 if(fCopyTouchables) { << 162 fPosition = right.fPosition; 98 fpTouchable = right.fpTouchable; << 163 fGlobalTime = right.fGlobalTime; 99 fpNextTouchable = right.fpNextTouchable; << 164 fLocalTime = right.fLocalTime; 100 fpOriginTouchable = right.fpOriginToucha << 165 fTrackLength = right.fTrackLength; 101 } << 166 fWeight = right.fWeight; 102 << 167 fStepLength = right.fStepLength; 103 // Track ID (and Parent ID) is not copied << 168 104 fTrackID = 0; << 169 // Track ID (and Parent ID) is not copied and set to zero for new track 105 fParentID = 0; << 170 fTrackID = 0; 106 << 171 fParentID =0; 107 // CurrentStepNumber is set to be 0 << 172 108 fCurrentStepNumber = 0; << 173 // CurrentStepNumber is set to be 0 109 << 174 fCurrentStepNumber = 0; 110 // Creator model ID << 175 111 fCreatorModelID = right.fCreatorModelID; << 176 // velocity information 112 << 177 fVelocity = right.fVelocity; 113 // Parent resonance << 178 114 fParentResonanceDef = right.fParentResonan << 179 // dynamic particle information 115 fParentResonanceID = right.fParentResonan << 180 fpDynamicParticle = new G4DynamicParticle(*(right.fpDynamicParticle)); >> 181 >> 182 // track status and flags for tracking >> 183 fTrackStatus = right.fTrackStatus; >> 184 fBelowThreshold = right.fBelowThreshold; >> 185 fGoodForTracking = right.fGoodForTracking; 116 186 117 // velocity information << 187 // Step information (Step Length, Step Number, pointer to the Step,) 118 fVelocity = right.fVelocity; << 188 // are not copied 119 << 189 fpStep=0; 120 // dynamic particle information << 190 121 delete fpDynamicParticle; << 191 // vertex information 122 fpDynamicParticle = new G4DynamicParticle( << 192 fVtxPosition = right.fVtxPosition; 123 << 193 fpLVAtVertex = right.fpLVAtVertex; 124 // track status and flags for tracking << 194 fVtxKineticEnergy = right.fVtxKineticEnergy; 125 fTrackStatus = right.fTrackStatus; << 195 fVtxMomentumDirection = right.fVtxMomentumDirection; 126 fBelowThreshold = right.fBelowThreshold; << 196 127 fGoodForTracking = right.fGoodForTracking; << 197 // CreatorProcess and UserInformation are not copied 128 << 198 fpCreatorProcess = 0; 129 // Step information (Step Length, Step Num << 199 fpUserInformation = 0; 130 // are not copied << 200 131 fpStep = nullptr; << 201 prev_mat = right.prev_mat; 132 << 202 groupvel = right.groupvel; 133 // vertex information << 203 prev_velocity = right.prev_velocity; 134 fVtxPosition = right.fVtxPosition << 204 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 205 149 is_OpticalPhoton = right.is_OpticalPhoton; << 206 is_OpticalPhoton = right.is_OpticalPhoton; 150 useGivenVelocity = right.useGivenVelocity; << 207 useGivenVelocity = right.useGivenVelocity; 151 208 152 ClearAuxiliaryTrackInformation(); << 209 fpAuxiliaryTrackInformationMap = 0; 153 } 210 } 154 return *this; 211 return *this; 155 } 212 } 156 213 157 // ------------------------------------------- << 214 /////////////////// 158 void G4Track::CopyTrackInfo(const G4Track& rig << 215 void G4Track::CopyTrackInfo(const G4Track& right) >> 216 ////////////////// 159 { 217 { 160 fCopyTouchables = copyTouchables; << 161 *this = right; 218 *this = right; 162 fCopyTouchables = true; << 163 } 219 } 164 220 165 // ------------------------------------------- << 221 /////////////////// >> 222 G4double G4Track::CalculateVelocity() const >> 223 /////////////////// >> 224 { >> 225 if (useGivenVelocity) return fVelocity; >> 226 >> 227 G4double velocity = c_light ; >> 228 >> 229 G4double mass = fpDynamicParticle->GetMass(); >> 230 >> 231 // special case for photons >> 232 if ( is_OpticalPhoton ) return CalculateVelocityForOpticalPhoton(); >> 233 >> 234 // particles other than optical photon >> 235 if (mass<DBL_MIN) { >> 236 // Zero Mass >> 237 velocity = c_light; >> 238 } else { >> 239 G4double T = (fpDynamicParticle->GetKineticEnergy())/mass; >> 240 if (T > GetMaxTOfVelocityTable()) { >> 241 velocity = c_light; >> 242 } else if (T<DBL_MIN) { >> 243 velocity =0.; >> 244 } else if (T<GetMinTOfVelocityTable()) { >> 245 velocity = c_light*std::sqrt(T*(T+2.))/(T+1.0); >> 246 } else { >> 247 velocity = velTable->Value(T); >> 248 } >> 249 >> 250 } >> 251 return velocity ; >> 252 } >> 253 >> 254 /////////////////// 166 G4double G4Track::CalculateVelocityForOpticalP 255 G4double G4Track::CalculateVelocityForOpticalPhoton() const >> 256 /////////////////// 167 { 257 { 168 G4double velocity = c_light; << 258 >> 259 G4double velocity = c_light ; >> 260 169 261 170 G4Material* mat = nullptr; << 262 G4Material* mat=0; 171 G4bool update_groupvel = false; 263 G4bool update_groupvel = false; 172 if(fpStep != nullptr) << 264 if ( fpStep !=0 ){ 173 { << 265 mat= this->GetMaterial(); // Fix for repeated volumes 174 mat = this->GetMaterial(); // Fix for r << 266 }else{ 175 } << 267 if (fpTouchable!=0){ 176 else << 268 mat=fpTouchable->GetVolume()->GetLogicalVolume()->GetMaterial(); 177 { << 178 if(fpTouchable) << 179 { << 180 mat = fpTouchable->GetVolume()->GetLogic << 181 } 269 } 182 } 270 } 183 // check if previous step is in the same vol 271 // check if previous step is in the same volume 184 // and get new GROUPVELOCITY table if neces << 272 // and get new GROUPVELOCITY table if necessary 185 if((mat != nullptr) && ((mat != prev_mat) || << 273 if ((mat != 0) && ((mat != prev_mat)||(groupvel==0))) { 186 { << 274 groupvel = 0; 187 groupvel = nullptr; << 275 if(mat->GetMaterialPropertiesTable() != 0) 188 if(mat->GetMaterialPropertiesTable() != nu << 276 groupvel = mat->GetMaterialPropertiesTable()->GetProperty("GROUPVEL"); 189 groupvel = mat->GetMaterialPropertiesTab << 190 update_groupvel = true; 277 update_groupvel = true; 191 } 278 } 192 prev_mat = mat; 279 prev_mat = mat; 193 << 280 194 if(groupvel != nullptr) << 281 if (groupvel != 0 ) { 195 { << 196 // light velocity = c/(rindex+d(rindex)/d( 282 // light velocity = c/(rindex+d(rindex)/d(log(E_phot))) 197 // values stored in GROUPVEL material prop 283 // values stored in GROUPVEL material properties vector 198 velocity = prev_velocity; << 284 velocity = prev_velocity; 199 << 285 200 // check if momentum is same as in the pre 286 // check if momentum is same as in the previous step 201 // and calculate group velocity if necess << 287 // and calculate group velocity if necessary 202 G4double current_momentum = fpDynamicParti 288 G4double current_momentum = fpDynamicParticle->GetTotalMomentum(); 203 if(update_groupvel || (current_momentum != << 289 if( update_groupvel || (current_momentum != prev_momentum) ) { 204 { << 290 velocity = 205 velocity = groupvel->Value(current_ << 291 groupvel->Value(current_momentum); 206 prev_velocity = velocity; 292 prev_velocity = velocity; 207 prev_momentum = current_momentum; 293 prev_momentum = current_momentum; 208 } 294 } 209 } << 295 } >> 296 >> 297 return velocity ; >> 298 } 210 299 211 return velocity; << 300 /////////////////// >> 301 void G4Track::SetVelocityTableProperties(G4double t_max, G4double t_min, G4int nbin) >> 302 /////////////////// >> 303 { >> 304 G4VelocityTable::SetVelocityTableProperties(t_max, t_min, nbin); >> 305 velTable = G4VelocityTable::GetVelocityTable(); 212 } 306 } 213 307 214 // ------------------------------------------- << 308 /////////////////// 215 void G4Track::SetAuxiliaryTrackInformation(G4i << 309 G4double G4Track::GetMaxTOfVelocityTable() 216 G4VAuxiliaryTrackInformation* in << 310 /////////////////// >> 311 { return G4VelocityTable::GetMaxTOfVelocityTable(); } >> 312 >> 313 /////////////////// >> 314 G4double G4Track::GetMinTOfVelocityTable() >> 315 /////////////////// >> 316 { return G4VelocityTable::GetMinTOfVelocityTable(); } >> 317 >> 318 /////////////////// >> 319 G4int G4Track::GetNbinOfVelocityTable() >> 320 /////////////////// >> 321 { return G4VelocityTable::GetNbinOfVelocityTable(); } >> 322 >> 323 /////////////////// >> 324 void G4Track::SetAuxiliaryTrackInformation(G4int idx, G4VAuxiliaryTrackInformation* info) const >> 325 /////////////////// 217 { 326 { 218 if(fpAuxiliaryTrackInformationMap == nullptr << 327 if(!fpAuxiliaryTrackInformationMap) 219 { << 328 { fpAuxiliaryTrackInformationMap = new std::map<G4int,G4VAuxiliaryTrackInformation*>; } 220 fpAuxiliaryTrackInformationMap = << 329 if(idx<0 || idx>=G4PhysicsModelCatalog::Entries()) 221 new std::map<G4int, G4VAuxiliaryTrackInf << 222 } << 223 if(G4PhysicsModelCatalog::GetModelIndex(id) << 224 { 330 { 225 G4ExceptionDescription ED; 331 G4ExceptionDescription ED; 226 ED << "Process/model ID <" << id << "> is << 332 ED << "Process/model index <" << idx << "> is invalid."; 227 G4Exception("G4VAuxiliaryTrackInformation: << 333 G4Exception("G4VAuxiliaryTrackInformation::G4VAuxiliaryTrackInformation()","TRACK0982", 228 "TRACK0982", FatalException, E << 334 FatalException, ED); 229 } 335 } 230 (*fpAuxiliaryTrackInformationMap)[id] = info << 336 (*fpAuxiliaryTrackInformationMap)[idx] = info; 231 } 337 } 232 338 233 // ------------------------------------------- << 339 /////////////////// 234 G4VAuxiliaryTrackInformation* << 340 G4VAuxiliaryTrackInformation* G4Track::GetAuxiliaryTrackInformation(G4int idx) const 235 G4Track::GetAuxiliaryTrackInformation(G4int id << 341 /////////////////// 236 { 342 { 237 if(fpAuxiliaryTrackInformationMap == nullptr << 343 if(!fpAuxiliaryTrackInformationMap) return 0; 238 return nullptr; << 344 std::map<G4int,G4VAuxiliaryTrackInformation*>::iterator itr 239 << 345 = fpAuxiliaryTrackInformationMap->find(idx); 240 auto itr = fpAuxiliaryTrackInformationMap->f << 346 if(itr==fpAuxiliaryTrackInformationMap->end()) return 0; 241 if(itr == fpAuxiliaryTrackInformationMap->ce << 347 else return (*itr).second; 242 return nullptr; << 243 return (*itr).second; << 244 } 348 } 245 349 246 // ------------------------------------------- << 350 /////////////////// 247 void G4Track::RemoveAuxiliaryTrackInformation( << 351 void G4Track::RemoveAuxiliaryTrackInformation(G4int idx) >> 352 /////////////////// 248 { 353 { 249 if(fpAuxiliaryTrackInformationMap != nullptr << 354 if(fpAuxiliaryTrackInformationMap && idx>=0 && idx<G4PhysicsModelCatalog::Entries()) 250 fpAuxiliaryTrackInformationMap->find(id) << 355 fpAuxiliaryTrackInformationMap->erase(idx); 251 { << 252 fpAuxiliaryTrackInformationMap->erase(id); << 253 } << 254 } 356 } 255 357 256 // ------------------------------------------- << 358 /////////////////// 257 void G4Track::RemoveAuxiliaryTrackInformation( 359 void G4Track::RemoveAuxiliaryTrackInformation(G4String& name) >> 360 /////////////////// 258 { 361 { 259 if(fpAuxiliaryTrackInformationMap != nullptr << 362 if(fpAuxiliaryTrackInformationMap) 260 { 363 { 261 G4int id = G4PhysicsModelCatalog::GetModel << 364 G4int idx = G4PhysicsModelCatalog::GetIndex(name); 262 RemoveAuxiliaryTrackInformation(id); << 365 RemoveAuxiliaryTrackInformation(idx); 263 } 366 } 264 } 367 } 265 368 266 // ------------------------------------------- << 369 /////////////////// 267 void G4Track::ClearAuxiliaryTrackInformation() 370 void G4Track::ClearAuxiliaryTrackInformation() >> 371 /////////////////// 268 { 372 { 269 if(fpAuxiliaryTrackInformationMap == nullptr << 373 if(!fpAuxiliaryTrackInformationMap) return; 270 return; << 374 for(std::map<G4int,G4VAuxiliaryTrackInformation*>::iterator itr=fpAuxiliaryTrackInformationMap->begin(); 271 for(const auto& itr : *fpAuxiliaryTrackInfor << 375 itr!=fpAuxiliaryTrackInformationMap->end(); itr++) 272 { << 376 { delete (*itr).second; } 273 delete itr.second; << 274 } << 275 delete fpAuxiliaryTrackInformationMap; 377 delete fpAuxiliaryTrackInformationMap; 276 fpAuxiliaryTrackInformationMap = nullptr; << 378 fpAuxiliaryTrackInformationMap = 0; 277 } 379 } >> 380 >> 381 278 382