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