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 // 26 // >> 27 // $Id: G4eLowEnergyLoss.cc 107396 2017-11-10 08:28:08Z gcosmo $ >> 28 // GEANT4 tag $Name: geant4-09-01-ref-00 $ 27 // 29 // 28 // ------------------------------------------- 30 // ----------------------------------------------------------- 29 // GEANT 4 class implementation file 31 // GEANT 4 class implementation file 30 // 32 // 31 // History: based on object model of 33 // History: based on object model of 32 // 2nd December 1995, G.Cosmo 34 // 2nd December 1995, G.Cosmo 33 // ---------- G4eLowEnergyLoss physics pr 35 // ---------- G4eLowEnergyLoss physics process ----------- 34 // by Laszlo Urban, 20 March 19 36 // by Laszlo Urban, 20 March 1997 35 // ******************************************* 37 // ************************************************************** 36 // It calculates the energy loss of e+/e-. 38 // It calculates the energy loss of e+/e-. 37 // ------------------------------------------- 39 // -------------------------------------------------------------- 38 // 40 // 39 // 08-05-97: small changes by L.Urban 41 // 08-05-97: small changes by L.Urban 40 // 27-05-98: several bugs and inconsistencies 42 // 27-05-98: several bugs and inconsistencies are corrected, 41 // new table (the inverse of the ran 43 // new table (the inverse of the range table) added , 42 // AlongStepDoit uses now this new t 44 // AlongStepDoit uses now this new table. L.Urban 43 // 08-09-98: cleanup 45 // 08-09-98: cleanup 44 // 24-09-98: rndmStepFlag false by default (no 46 // 24-09-98: rndmStepFlag false by default (no randomization of the step) 45 // 14-10-98: messenger file added. 47 // 14-10-98: messenger file added. 46 // 16-10-98: public method SetStepFunction() 48 // 16-10-98: public method SetStepFunction() 47 // 20-01-99: important correction in AlongStep 49 // 20-01-99: important correction in AlongStepDoIt , L.Urban 48 // 10/02/00 modifications , new e.m. structur 50 // 10/02/00 modifications , new e.m. structure, L.Urban 49 // 11/04/00: Bug fix in dE/dx fluctuation simu 51 // 11/04/00: Bug fix in dE/dx fluctuation simulation, Veronique Lefebure 50 // 19-09-00 change of fluctuation sampling V. 52 // 19-09-00 change of fluctuation sampling V.Ivanchenko 51 // 20/09/00 update fluctuations V.Ivanchenko 53 // 20/09/00 update fluctuations V.Ivanchenko 52 // 18/10/01 add fluorescence AlongStepDoIt V. 54 // 18/10/01 add fluorescence AlongStepDoIt V.Ivanchenko 53 // 18/10/01 Revision to improve code quality 55 // 18/10/01 Revision to improve code quality and consistency with design, MGP 54 // 19/10/01 update according to new design, V 56 // 19/10/01 update according to new design, V.Ivanchenko 55 // 24/10/01 MGP - Protection against negative 57 // 24/10/01 MGP - Protection against negative energy loss in AlongStepDoIt 56 // 26/10/01 VI Clean up access to deexcitatio 58 // 26/10/01 VI Clean up access to deexcitation 57 // 23/11/01 VI Move static member-functions f 59 // 23/11/01 VI Move static member-functions from header to source 58 // 28/05/02 VI Remove flag fStopAndKill 60 // 28/05/02 VI Remove flag fStopAndKill 59 // 03/06/02 MGP - Restore fStopAndKill 61 // 03/06/02 MGP - Restore fStopAndKill 60 // 28/10/02 VI Optimal binning for dE/dx 62 // 28/10/02 VI Optimal binning for dE/dx 61 // 21/01/03 VI cut per region 63 // 21/01/03 VI cut per region 62 // 01/06/04 VI check if stopped particle has 64 // 01/06/04 VI check if stopped particle has AtRest processes 63 // 65 // 64 // ------------------------------------------- 66 // -------------------------------------------------------------- 65 67 66 #include "G4eLowEnergyLoss.hh" 68 #include "G4eLowEnergyLoss.hh" 67 #include "G4SystemOfUnits.hh" 69 #include "G4SystemOfUnits.hh" 68 #include "" 70 #include "" 69 #include "G4Poisson.hh" 71 #include "G4Poisson.hh" 70 #include "G4ProductionCutsTable.hh" 72 #include "G4ProductionCutsTable.hh" 71 73 72 // 74 // 73 75 74 // Initialisation of static data members 76 // Initialisation of static data members 75 // ------------------------------------- 77 // ------------------------------------- 76 // Contributing processes : ion.loss + soft br 78 // Contributing processes : ion.loss + soft brems->NbOfProcesses is initialized 77 // to 2 . YOU DO NOT HAVE TO CHANGE this varia 79 // to 2 . YOU DO NOT HAVE TO CHANGE this variable for a 'normal' run. 78 // 80 // 79 // You have to change NbOfProcesses if you inv 81 // You have to change NbOfProcesses if you invent a new process contributing 80 // to the continuous energy loss. 82 // to the continuous energy loss. 81 // The NbOfProcesses data member can be change 83 // The NbOfProcesses data member can be changed using the (public static) 82 // functions Get/Set/Plus/MinusNbOfProcesses ( 84 // functions Get/Set/Plus/MinusNbOfProcesses (see G4eLowEnergyLoss.hh) 83 85 84 G4int G4eLowEnergyLoss::NbOfProcess 86 G4int G4eLowEnergyLoss::NbOfProcesses = 2; 85 87 86 G4int G4eLowEnergyLoss::CounterOfEl 88 G4int G4eLowEnergyLoss::CounterOfElectronProcess = 0; 87 G4int G4eLowEnergyLoss::CounterOfPo 89 G4int G4eLowEnergyLoss::CounterOfPositronProcess = 0; 88 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfE 90 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfElectronProcess = 89 new 91 new G4PhysicsTable*[10]; 90 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfP 92 G4PhysicsTable** G4eLowEnergyLoss::RecorderOfPositronProcess = 91 new 93 new G4PhysicsTable*[10]; 92 94 93 95 94 G4PhysicsTable* G4eLowEnergyLoss::theDEDXElec 96 G4PhysicsTable* G4eLowEnergyLoss::theDEDXElectronTable = 0; 95 G4PhysicsTable* G4eLowEnergyLoss::theDEDXPosi 97 G4PhysicsTable* G4eLowEnergyLoss::theDEDXPositronTable = 0; 96 G4PhysicsTable* G4eLowEnergyLoss::theRangeEle 98 G4PhysicsTable* G4eLowEnergyLoss::theRangeElectronTable = 0; 97 G4PhysicsTable* G4eLowEnergyLoss::theRangePos 99 G4PhysicsTable* G4eLowEnergyLoss::theRangePositronTable = 0; 98 G4PhysicsTable* G4eLowEnergyLoss::theInverseR 100 G4PhysicsTable* G4eLowEnergyLoss::theInverseRangeElectronTable = 0; 99 G4PhysicsTable* G4eLowEnergyLoss::theInverseR 101 G4PhysicsTable* G4eLowEnergyLoss::theInverseRangePositronTable = 0; 100 G4PhysicsTable* G4eLowEnergyLoss::theLabTimeE 102 G4PhysicsTable* G4eLowEnergyLoss::theLabTimeElectronTable = 0; 101 G4PhysicsTable* G4eLowEnergyLoss::theLabTimeP 103 G4PhysicsTable* G4eLowEnergyLoss::theLabTimePositronTable = 0; 102 G4PhysicsTable* G4eLowEnergyLoss::theProperTi 104 G4PhysicsTable* G4eLowEnergyLoss::theProperTimeElectronTable = 0; 103 G4PhysicsTable* G4eLowEnergyLoss::theProperTi 105 G4PhysicsTable* G4eLowEnergyLoss::theProperTimePositronTable = 0; 104 106 105 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCo 107 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCoeffATable = 0; 106 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCo 108 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCoeffBTable = 0; 107 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCo 109 G4PhysicsTable* G4eLowEnergyLoss::theeRangeCoeffCTable = 0; 108 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCo 110 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCoeffATable = 0; 109 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCo 111 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCoeffBTable = 0; 110 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCo 112 G4PhysicsTable* G4eLowEnergyLoss::thepRangeCoeffCTable = 0; 111 113 112 G4double G4eLowEnergyLoss::LowerBoundE 114 G4double G4eLowEnergyLoss::LowerBoundEloss = 10.*eV ; 113 G4double G4eLowEnergyLoss::UpperBoundE 115 G4double G4eLowEnergyLoss::UpperBoundEloss = 100.*GeV ; 114 G4int G4eLowEnergyLoss::NbinEloss = 116 G4int G4eLowEnergyLoss::NbinEloss = 360 ; 115 G4double G4eLowEnergyLoss::RTable ; 117 G4double G4eLowEnergyLoss::RTable ; 116 G4double G4eLowEnergyLoss::LOGRTable ; 118 G4double G4eLowEnergyLoss::LOGRTable ; 117 119 118 120 119 G4EnergyLossMessenger* G4eLowEnergyLoss::eLoss 121 G4EnergyLossMessenger* G4eLowEnergyLoss::eLossMessenger = 0; 120 122 121 // 123 // 122 124 123 // constructor and destructor 125 // constructor and destructor 124 126 125 G4eLowEnergyLoss::G4eLowEnergyLoss(const G4Str 127 G4eLowEnergyLoss::G4eLowEnergyLoss(const G4String& processName) 126 : G4RDVeLowEnergyLoss (processName), 128 : G4RDVeLowEnergyLoss (processName), 127 theLossTable(0), 129 theLossTable(0), 128 MinKineticEnergy(1.*eV), 130 MinKineticEnergy(1.*eV), 129 Charge(-1.), 131 Charge(-1.), 130 lastCharge(0.), 132 lastCharge(0.), 131 theDEDXTable(0), 133 theDEDXTable(0), 132 CounterOfProcess(0), 134 CounterOfProcess(0), 133 RecorderOfProcess(0), 135 RecorderOfProcess(0), 134 fdEdx(0), 136 fdEdx(0), 135 fRangeNow(0), 137 fRangeNow(0), 136 linLossLimit(0.05), 138 linLossLimit(0.05), 137 theFluo(false) 139 theFluo(false) 138 { 140 { 139 141 140 //create (only once) EnergyLoss messenger 142 //create (only once) EnergyLoss messenger 141 if(!eLossMessenger) eLossMessenger = new G4En 143 if(!eLossMessenger) eLossMessenger = new G4EnergyLossMessenger(); 142 } 144 } 143 145 144 // 146 // 145 147 146 G4eLowEnergyLoss::~G4eLowEnergyLoss() 148 G4eLowEnergyLoss::~G4eLowEnergyLoss() 147 { 149 { 148 if (theLossTable) 150 if (theLossTable) 149 { 151 { 150 theLossTable->clearAndDestroy(); 152 theLossTable->clearAndDestroy(); 151 delete theLossTable; 153 delete theLossTable; 152 } 154 } 153 } 155 } 154 156 155 void G4eLowEnergyLoss::SetNbOfProcesses(G4int 157 void G4eLowEnergyLoss::SetNbOfProcesses(G4int nb) 156 { 158 { 157 NbOfProcesses=nb; 159 NbOfProcesses=nb; 158 } 160 } 159 161 160 void G4eLowEnergyLoss::PlusNbOfProcesses() 162 void G4eLowEnergyLoss::PlusNbOfProcesses() 161 { 163 { 162 NbOfProcesses++; 164 NbOfProcesses++; 163 } 165 } 164 166 165 void G4eLowEnergyLoss::MinusNbOfProcesses() 167 void G4eLowEnergyLoss::MinusNbOfProcesses() 166 { 168 { 167 NbOfProcesses--; 169 NbOfProcesses--; 168 } 170 } 169 171 170 G4int G4eLowEnergyLoss::GetNbOfProcesses() 172 G4int G4eLowEnergyLoss::GetNbOfProcesses() 171 { 173 { 172 return NbOfProcesses; 174 return NbOfProcesses; 173 } 175 } 174 176 175 void G4eLowEnergyLoss::SetLowerBoundEloss(G4do 177 void G4eLowEnergyLoss::SetLowerBoundEloss(G4double val) 176 { 178 { 177 LowerBoundEloss=val; 179 LowerBoundEloss=val; 178 } 180 } 179 181 180 void G4eLowEnergyLoss::SetUpperBoundEloss(G4do 182 void G4eLowEnergyLoss::SetUpperBoundEloss(G4double val) 181 { 183 { 182 UpperBoundEloss=val; 184 UpperBoundEloss=val; 183 } 185 } 184 186 185 void G4eLowEnergyLoss::SetNbinEloss(G4int nb) 187 void G4eLowEnergyLoss::SetNbinEloss(G4int nb) 186 { 188 { 187 NbinEloss=nb; 189 NbinEloss=nb; 188 } 190 } 189 191 190 G4double G4eLowEnergyLoss::GetLowerBoundEloss( 192 G4double G4eLowEnergyLoss::GetLowerBoundEloss() 191 { 193 { 192 return LowerBoundEloss; 194 return LowerBoundEloss; 193 } 195 } 194 196 195 G4double G4eLowEnergyLoss::GetUpperBoundEloss( 197 G4double G4eLowEnergyLoss::GetUpperBoundEloss() 196 { 198 { 197 return UpperBoundEloss; 199 return UpperBoundEloss; 198 } 200 } 199 201 200 G4int G4eLowEnergyLoss::GetNbinEloss() 202 G4int G4eLowEnergyLoss::GetNbinEloss() 201 { 203 { 202 return NbinEloss; 204 return NbinEloss; 203 } 205 } 204 // 206 // 205 207 206 void G4eLowEnergyLoss::BuildDEDXTable( 208 void G4eLowEnergyLoss::BuildDEDXTable( 207 const G4ParticleDefin 209 const G4ParticleDefinition& aParticleType) 208 { 210 { 209 ParticleMass = aParticleType.GetPDGMass(); 211 ParticleMass = aParticleType.GetPDGMass(); 210 Charge = aParticleType.GetPDGCharge()/eplus; 212 Charge = aParticleType.GetPDGCharge()/eplus; 211 213 212 // calculate data members LOGRTable,RTable 214 // calculate data members LOGRTable,RTable first 213 215 214 G4double lrate = std::log(UpperBoundEloss/Lo 216 G4double lrate = std::log(UpperBoundEloss/LowerBoundEloss); 215 LOGRTable=lrate/NbinEloss; 217 LOGRTable=lrate/NbinEloss; 216 RTable =std::exp(LOGRTable); 218 RTable =std::exp(LOGRTable); 217 // Build energy loss table as a sum of the e 219 // Build energy loss table as a sum of the energy loss due to the 218 // different processes. 220 // different processes. 219 // 221 // 220 222 221 const G4ProductionCutsTable* theCoupleTable= 223 const G4ProductionCutsTable* theCoupleTable= 222 G4ProductionCutsTable::GetProductionCu 224 G4ProductionCutsTable::GetProductionCutsTable(); 223 size_t numOfCouples = theCoupleTable->GetTab 225 size_t numOfCouples = theCoupleTable->GetTableSize(); 224 226 225 // create table for the total energy loss 227 // create table for the total energy loss 226 228 227 if (&aParticleType==G4Electron::Electron()) 229 if (&aParticleType==G4Electron::Electron()) 228 { 230 { 229 RecorderOfProcess=RecorderOfElectronProc 231 RecorderOfProcess=RecorderOfElectronProcess; 230 CounterOfProcess=CounterOfElectronProces 232 CounterOfProcess=CounterOfElectronProcess; 231 if (CounterOfProcess == NbOfProcesses) 233 if (CounterOfProcess == NbOfProcesses) 232 { 234 { 233 if (theDEDXElectronTable) 235 if (theDEDXElectronTable) 234 { 236 { 235 theDEDXElectronTable->clearAndDes 237 theDEDXElectronTable->clearAndDestroy(); 236 delete theDEDXElectronTable; 238 delete theDEDXElectronTable; 237 } 239 } 238 theDEDXElectronTable = new G4PhysicsT 240 theDEDXElectronTable = new G4PhysicsTable(numOfCouples); 239 theDEDXTable = theDEDXElectronTable; 241 theDEDXTable = theDEDXElectronTable; 240 } 242 } 241 } 243 } 242 if (&aParticleType==G4Positron::Positron()) 244 if (&aParticleType==G4Positron::Positron()) 243 { 245 { 244 RecorderOfProcess=RecorderOfPositronProce 246 RecorderOfProcess=RecorderOfPositronProcess; 245 CounterOfProcess=CounterOfPositronProcess 247 CounterOfProcess=CounterOfPositronProcess; 246 if (CounterOfProcess == NbOfProcesses) 248 if (CounterOfProcess == NbOfProcesses) 247 { 249 { 248 if (theDEDXPositronTable) 250 if (theDEDXPositronTable) 249 { 251 { 250 theDEDXPositronTable->clearAndDest 252 theDEDXPositronTable->clearAndDestroy(); 251 delete theDEDXPositronTable; 253 delete theDEDXPositronTable; 252 } 254 } 253 theDEDXPositronTable = new G4PhysicsTa 255 theDEDXPositronTable = new G4PhysicsTable(numOfCouples); 254 theDEDXTable = theDEDXPositronTable; 256 theDEDXTable = theDEDXPositronTable; 255 } 257 } 256 } 258 } 257 259 258 if (CounterOfProcess == NbOfProcesses) 260 if (CounterOfProcess == NbOfProcesses) 259 { 261 { 260 // fill the tables 262 // fill the tables 261 // loop for materials 263 // loop for materials 262 G4double LowEdgeEnergy , Value; 264 G4double LowEdgeEnergy , Value; 263 G4bool isOutRange; 265 G4bool isOutRange; 264 G4PhysicsTable* pointer; 266 G4PhysicsTable* pointer; 265 267 266 for (size_t J=0; J<numOfCouples; J++) 268 for (size_t J=0; J<numOfCouples; J++) 267 { 269 { 268 // create physics vector and fill it 270 // create physics vector and fill it 269 271 270 G4PhysicsLogVector* aVector = new G4P 272 G4PhysicsLogVector* aVector = new G4PhysicsLogVector( 271 LowerBoundEloss, UpperBoun 273 LowerBoundEloss, UpperBoundEloss, NbinEloss); 272 274 273 // loop for the kinetic energy 275 // loop for the kinetic energy 274 for (G4int i=0; i<NbinEloss; i++) 276 for (G4int i=0; i<NbinEloss; i++) 275 { 277 { 276 LowEdgeEnergy = aVector->GetLowE 278 LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ; 277 //here comes the sum of the diff 279 //here comes the sum of the different tables created by the 278 //processes (ionisation,bremsstr 280 //processes (ionisation,bremsstrahlung,etc...) 279 Value = 0.; 281 Value = 0.; 280 for (G4int process=0; process < 282 for (G4int process=0; process < NbOfProcesses; process++) 281 { 283 { 282 pointer= RecorderOfProcess[ 284 pointer= RecorderOfProcess[process]; 283 Value += (*pointer)[J]->Get 285 Value += (*pointer)[J]->GetValue(LowEdgeEnergy,isOutRange); 284 } 286 } 285 287 286 aVector->PutValue(i,Value) ; 288 aVector->PutValue(i,Value) ; 287 } 289 } 288 290 289 theDEDXTable->insert(aVector) ; 291 theDEDXTable->insert(aVector) ; 290 292 291 } 293 } 292 294 293 295 294 //reset counter to zero 296 //reset counter to zero 295 if (&aParticleType==G4Electron::Electron( 297 if (&aParticleType==G4Electron::Electron()) CounterOfElectronProcess=0; 296 if (&aParticleType==G4Positron::Positron( 298 if (&aParticleType==G4Positron::Positron()) CounterOfPositronProcess=0; 297 299 298 ParticleMass = aParticleType.GetPDGMass() 300 ParticleMass = aParticleType.GetPDGMass(); 299 301 300 if (&aParticleType==G4Electron::Electron( 302 if (&aParticleType==G4Electron::Electron()) 301 { 303 { 302 // Build range table 304 // Build range table 303 theRangeElectronTable = BuildRangeTable 305 theRangeElectronTable = BuildRangeTable(theDEDXElectronTable, 304 306 theRangeElectronTable, 305 LowerBoundEloss, 307 LowerBoundEloss,UpperBoundEloss,NbinEloss); 306 308 307 // Build lab/proper time tables 309 // Build lab/proper time tables 308 theLabTimeElectronTable = BuildLabTimeT 310 theLabTimeElectronTable = BuildLabTimeTable(theDEDXElectronTable, 309 theLabTimeElectronTab 311 theLabTimeElectronTable, 310 LowerBoundEloss,Upper 312 LowerBoundEloss,UpperBoundEloss,NbinEloss); 311 theProperTimeElectronTable = BuildPrope 313 theProperTimeElectronTable = BuildProperTimeTable(theDEDXElectronTable, 312 theProperTimeElect 314 theProperTimeElectronTable, 313 LowerBoundEloss,Up 315 LowerBoundEloss,UpperBoundEloss,NbinEloss); 314 316 315 // Build coeff tables for the energy lo 317 // Build coeff tables for the energy loss calculation 316 theeRangeCoeffATable = BuildRangeCoeffA 318 theeRangeCoeffATable = BuildRangeCoeffATable(theRangeElectronTable, 317 theeRangeCoeffATa 319 theeRangeCoeffATable, 318 LowerBoundEloss,U 320 LowerBoundEloss,UpperBoundEloss,NbinEloss); 319 321 320 theeRangeCoeffBTable = BuildRangeCoeffB 322 theeRangeCoeffBTable = BuildRangeCoeffBTable(theRangeElectronTable, 321 theeRangeCoeffBTa 323 theeRangeCoeffBTable, 322 LowerBoundEloss,U 324 LowerBoundEloss,UpperBoundEloss,NbinEloss); 323 325 324 theeRangeCoeffCTable = BuildRangeCoeffC 326 theeRangeCoeffCTable = BuildRangeCoeffCTable(theRangeElectronTable, 325 theeRangeCoeffCTa 327 theeRangeCoeffCTable, 326 LowerBoundEloss,U 328 LowerBoundEloss,UpperBoundEloss,NbinEloss); 327 329 328 // invert the range table 330 // invert the range table 329 theInverseRangeElectronTable = BuildInv 331 theInverseRangeElectronTable = BuildInverseRangeTable(theRangeElectronTable, 330 theeRangeCoeffAT 332 theeRangeCoeffATable, 331 theeRangeCoeffBT 333 theeRangeCoeffBTable, 332 theeRangeCoeffCT 334 theeRangeCoeffCTable, 333 theInverseRangeE 335 theInverseRangeElectronTable, 334 LowerBoundEloss, 336 LowerBoundEloss,UpperBoundEloss,NbinEloss); 335 } 337 } 336 if (&aParticleType==G4Positron::Positron( 338 if (&aParticleType==G4Positron::Positron()) 337 { 339 { 338 // Build range table 340 // Build range table 339 theRangePositronTable = BuildRangeTable 341 theRangePositronTable = BuildRangeTable(theDEDXPositronTable, 340 342 theRangePositronTable, 341 LowerBoundEloss, 343 LowerBoundEloss,UpperBoundEloss,NbinEloss); 342 344 343 345 344 // Build lab/proper time tables 346 // Build lab/proper time tables 345 theLabTimePositronTable = BuildLabTimeT 347 theLabTimePositronTable = BuildLabTimeTable(theDEDXPositronTable, 346 theLabTimePositronTab 348 theLabTimePositronTable, 347 LowerBoundEloss,Upper 349 LowerBoundEloss,UpperBoundEloss,NbinEloss); 348 theProperTimePositronTable = BuildPrope 350 theProperTimePositronTable = BuildProperTimeTable(theDEDXPositronTable, 349 theProperTimePosit 351 theProperTimePositronTable, 350 LowerBoundEloss,Up 352 LowerBoundEloss,UpperBoundEloss,NbinEloss); 351 353 352 // Build coeff tables for the energy lo 354 // Build coeff tables for the energy loss calculation 353 thepRangeCoeffATable = BuildRangeCoeffA 355 thepRangeCoeffATable = BuildRangeCoeffATable(theRangePositronTable, 354 thepRangeCoeffATa 356 thepRangeCoeffATable, 355 LowerBoundEloss,U 357 LowerBoundEloss,UpperBoundEloss,NbinEloss); 356 358 357 thepRangeCoeffBTable = BuildRangeCoeffB 359 thepRangeCoeffBTable = BuildRangeCoeffBTable(theRangePositronTable, 358 thepRangeCoeffBTa 360 thepRangeCoeffBTable, 359 LowerBoundEloss,U 361 LowerBoundEloss,UpperBoundEloss,NbinEloss); 360 362 361 thepRangeCoeffCTable = BuildRangeCoeffC 363 thepRangeCoeffCTable = BuildRangeCoeffCTable(theRangePositronTable, 362 thepRangeCoeffCTa 364 thepRangeCoeffCTable, 363 LowerBoundEloss,U 365 LowerBoundEloss,UpperBoundEloss,NbinEloss); 364 366 365 // invert the range table 367 // invert the range table 366 theInverseRangePositronTable = BuildInv 368 theInverseRangePositronTable = BuildInverseRangeTable(theRangePositronTable, 367 thepRangeCoeffAT 369 thepRangeCoeffATable, 368 thepRangeCoeffBT 370 thepRangeCoeffBTable, 369 thepRangeCoeffCT 371 thepRangeCoeffCTable, 370 theInverseRangeP 372 theInverseRangePositronTable, 371 LowerBoundEloss, 373 LowerBoundEloss,UpperBoundEloss,NbinEloss); 372 } 374 } 373 375 374 if(verboseLevel > 1) { 376 if(verboseLevel > 1) { 375 G4cout << (*theDEDXElectronTable) << G4 377 G4cout << (*theDEDXElectronTable) << G4endl; 376 } 378 } 377 379 378 380 379 // make the energy loss and the range tab 381 // make the energy loss and the range table available 380 G4EnergyLossTables::Register(&aParticleTy 382 G4EnergyLossTables::Register(&aParticleType, 381 (&aParticleType==G4Electron::Electron() 383 (&aParticleType==G4Electron::Electron())? 382 theDEDXElectronTable: theDEDXPositronTa 384 theDEDXElectronTable: theDEDXPositronTable, 383 (&aParticleType==G4Electron::Electron() 385 (&aParticleType==G4Electron::Electron())? 384 theRangeElectronTable: theRangePositron 386 theRangeElectronTable: theRangePositronTable, 385 (&aParticleType==G4Electron::Electron() 387 (&aParticleType==G4Electron::Electron())? 386 theInverseRangeElectronTable: theInvers 388 theInverseRangeElectronTable: theInverseRangePositronTable, 387 (&aParticleType==G4Electron::Electron() 389 (&aParticleType==G4Electron::Electron())? 388 theLabTimeElectronTable: theLabTimePosi 390 theLabTimeElectronTable: theLabTimePositronTable, 389 (&aParticleType==G4Electron::Electron() 391 (&aParticleType==G4Electron::Electron())? 390 theProperTimeElectronTable: theProperTi 392 theProperTimeElectronTable: theProperTimePositronTable, 391 LowerBoundEloss, UpperBoundEloss, 1.,Nb 393 LowerBoundEloss, UpperBoundEloss, 1.,NbinEloss); 392 } 394 } 393 } 395 } 394 396 395 // 397 // 396 398 397 G4VParticleChange* G4eLowEnergyLoss::AlongStep 399 G4VParticleChange* G4eLowEnergyLoss::AlongStepDoIt( const G4Track& trackData, 398 400 const G4Step& stepData) 399 { 401 { 400 // compute the energy loss after a Step 402 // compute the energy loss after a Step 401 403 402 static const G4double faclow = 1.5 ; 404 static const G4double faclow = 1.5 ; 403 405 404 // get particle and material pointers from t 406 // get particle and material pointers from trackData 405 const G4DynamicParticle* aParticle = trackDa 407 const G4DynamicParticle* aParticle = trackData.GetDynamicParticle(); 406 G4double E = aParticle->GetKineticEnerg 408 G4double E = aParticle->GetKineticEnergy(); 407 409 408 const G4MaterialCutsCouple* couple = trackDa 410 const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple(); 409 411 410 G4double Step = stepData.GetStepLength(); 412 G4double Step = stepData.GetStepLength(); 411 413 412 aParticleChange.Initialize(trackData); 414 aParticleChange.Initialize(trackData); 413 //fParticleChange.Initialize(trackData); 415 //fParticleChange.Initialize(trackData); 414 416 415 G4double MeanLoss, finalT; 417 G4double MeanLoss, finalT; 416 418 417 if (E < MinKineticEnergy) finalT = 0.; 419 if (E < MinKineticEnergy) finalT = 0.; 418 420 419 else if ( E< faclow*LowerBoundEloss) 421 else if ( E< faclow*LowerBoundEloss) 420 { 422 { 421 if (Step >= fRangeNow) finalT = 0.; 423 if (Step >= fRangeNow) finalT = 0.; 422 // else finalT = E*(1.-Step/fRangeNow) ; 424 // else finalT = E*(1.-Step/fRangeNow) ; 423 else finalT = E*(1.-std::sqrt(Step/fRangeN 425 else finalT = E*(1.-std::sqrt(Step/fRangeNow)) ; 424 } 426 } 425 427 426 else if (E>=UpperBoundEloss) finalT = E - St 428 else if (E>=UpperBoundEloss) finalT = E - Step*fdEdx; 427 429 428 else if (Step >= fRangeNow) finalT = 0.; 430 else if (Step >= fRangeNow) finalT = 0.; 429 431 430 else 432 else 431 { 433 { 432 if(Step/fRangeNow < linLossLimit) finalT = 434 if(Step/fRangeNow < linLossLimit) finalT = E-Step*fdEdx ; 433 else 435 else 434 { 436 { 435 if (Charge<0.) finalT = G4EnergyLossTabl 437 if (Charge<0.) finalT = G4EnergyLossTables::GetPreciseEnergyFromRange 436 (G4Electron::Elec 438 (G4Electron::Electron(),fRangeNow-Step,couple); 437 else finalT = G4EnergyLossTabl 439 else finalT = G4EnergyLossTables::GetPreciseEnergyFromRange 438 (G4Positron::Posi 440 (G4Positron::Positron(),fRangeNow-Step,couple); 439 } 441 } 440 } 442 } 441 443 442 if(finalT < MinKineticEnergy) finalT = 0. ; 444 if(finalT < MinKineticEnergy) finalT = 0. ; 443 445 444 MeanLoss = E-finalT ; 446 MeanLoss = E-finalT ; 445 447 446 //now the loss with fluctuation 448 //now the loss with fluctuation 447 if ((EnlossFlucFlag) && (finalT > 0.) && (fi 449 if ((EnlossFlucFlag) && (finalT > 0.) && (finalT < E)&&(E > LowerBoundEloss)) 448 { 450 { 449 finalT = E-GetLossWithFluct(aParticle,coup 451 finalT = E-GetLossWithFluct(aParticle,couple,MeanLoss,Step); 450 if (finalT < 0.) finalT = 0.; 452 if (finalT < 0.) finalT = 0.; 451 } 453 } 452 454 453 // kill the particle if the kinetic energy < 455 // kill the particle if the kinetic energy <= 0 454 if (finalT <= 0. ) 456 if (finalT <= 0. ) 455 { 457 { 456 finalT = 0.; 458 finalT = 0.; 457 if(Charge > 0.0) aParticleChange.ProposeTr 459 if(Charge > 0.0) aParticleChange.ProposeTrackStatus(fStopButAlive); 458 else aParticleChange.ProposeTr 460 else aParticleChange.ProposeTrackStatus(fStopAndKill); 459 } 461 } 460 462 461 G4double edep = E - finalT; 463 G4double edep = E - finalT; 462 464 463 aParticleChange.ProposeEnergy(finalT); 465 aParticleChange.ProposeEnergy(finalT); 464 466 465 // Deexcitation of ionised atoms 467 // Deexcitation of ionised atoms 466 std::vector<G4DynamicParticle*>* deexcitatio 468 std::vector<G4DynamicParticle*>* deexcitationProducts = 0; 467 if (theFluo) deexcitationProducts = Deexcite 469 if (theFluo) deexcitationProducts = DeexciteAtom(couple,E,edep); 468 470 469 size_t nSecondaries = 0; 471 size_t nSecondaries = 0; 470 if (deexcitationProducts != 0) nSecondaries 472 if (deexcitationProducts != 0) nSecondaries = deexcitationProducts->size(); 471 aParticleChange.SetNumberOfSecondaries(nSeco 473 aParticleChange.SetNumberOfSecondaries(nSecondaries); 472 474 473 if (nSecondaries > 0) { 475 if (nSecondaries > 0) { 474 476 475 const G4StepPoint* preStep = stepData.GetP 477 const G4StepPoint* preStep = stepData.GetPreStepPoint(); 476 const G4StepPoint* postStep = stepData.Get 478 const G4StepPoint* postStep = stepData.GetPostStepPoint(); 477 G4ThreeVector r = preStep->GetPosition(); 479 G4ThreeVector r = preStep->GetPosition(); 478 G4ThreeVector deltaR = postStep->GetPositi 480 G4ThreeVector deltaR = postStep->GetPosition(); 479 deltaR -= r; 481 deltaR -= r; 480 G4double t = preStep->GetGlobalTime(); 482 G4double t = preStep->GetGlobalTime(); 481 G4double deltaT = postStep->GetGlobalTime( 483 G4double deltaT = postStep->GetGlobalTime(); 482 deltaT -= t; 484 deltaT -= t; 483 G4double time, q; 485 G4double time, q; 484 G4ThreeVector position; 486 G4ThreeVector position; 485 487 486 for (size_t i=0; i<nSecondaries; i++) { 488 for (size_t i=0; i<nSecondaries; i++) { 487 489 488 G4DynamicParticle* part = (*deexcitation 490 G4DynamicParticle* part = (*deexcitationProducts)[i]; 489 if (part != 0) { 491 if (part != 0) { 490 G4double eSecondary = part->GetKinetic 492 G4double eSecondary = part->GetKineticEnergy(); 491 edep -= eSecondary; 493 edep -= eSecondary; 492 if (edep > 0.) 494 if (edep > 0.) 493 { 495 { 494 q = G4UniformRand(); 496 q = G4UniformRand(); 495 time = deltaT*q + t; 497 time = deltaT*q + t; 496 position = deltaR*q; 498 position = deltaR*q; 497 position += r; 499 position += r; 498 G4Track* newTrack = new G4Track(part, ti 500 G4Track* newTrack = new G4Track(part, time, position); 499 aParticleChange.AddSecondary(newTrack); 501 aParticleChange.AddSecondary(newTrack); 500 } 502 } 501 else 503 else 502 { 504 { 503 edep += eSecondary; 505 edep += eSecondary; 504 delete part; 506 delete part; 505 part = 0; 507 part = 0; 506 } 508 } 507 } 509 } 508 } 510 } 509 } 511 } 510 delete deexcitationProducts; 512 delete deexcitationProducts; 511 513 512 aParticleChange.ProposeLocalEnergyDeposit(ed 514 aParticleChange.ProposeLocalEnergyDeposit(edep); 513 515 514 return &aParticleChange; 516 return &aParticleChange; 515 } 517 } 516 518 517 // 519 // 518 520 519 521