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 // $Id: G4ContinuousGainOfEnergy.cc,v 1.5 2010-11-11 11:51:56 ldesorgh Exp $ >> 27 // GEANT4 tag $Name: not supported by cvs2svn $ 26 // 28 // 27 << 28 #include "G4ContinuousGainOfEnergy.hh" 29 #include "G4ContinuousGainOfEnergy.hh" 29 << 30 #include "G4EmCorrections.hh" << 31 #include "G4LossTableManager.hh" << 32 #include "G4Material.hh" << 33 #include "G4MaterialCutsCouple.hh" << 34 #include "G4ParticleChange.hh" << 35 #include "G4ParticleDefinition.hh" << 36 #include "G4PhysicalConstants.hh" << 37 #include "G4Step.hh" 30 #include "G4Step.hh" 38 #include "G4SystemOfUnits.hh" << 31 #include "G4ParticleDefinition.hh" 39 #include "G4VEmFluctuationModel.hh" << 40 #include "G4VEmModel.hh" 32 #include "G4VEmModel.hh" 41 #include "G4VEnergyLossProcess.hh" << 33 #include "G4VEmFluctuationModel.hh" 42 #include "G4VParticleChange.hh" 34 #include "G4VParticleChange.hh" >> 35 #include "G4UnitsTable.hh" >> 36 #include "G4AdjointCSManager.hh" >> 37 #include "G4LossTableManager.hh" 43 38 44 ////////////////////////////////////////////// 39 /////////////////////////////////////////////////////// 45 G4ContinuousGainOfEnergy::G4ContinuousGainOfEn << 40 // 46 << 41 G4ContinuousGainOfEnergy::G4ContinuousGainOfEnergy(const G4String& name, 47 : G4VContinuousProcess(name, type) << 42 G4ProcessType type): G4VContinuousProcess(name, type) 48 {} << 43 { 49 44 50 ////////////////////////////////////////////// << 45 51 G4ContinuousGainOfEnergy::~G4ContinuousGainOfE << 46 linLossLimit=0.05; >> 47 lossFluctuationArePossible =true; >> 48 lossFluctuationFlag=true; >> 49 is_integral = false; >> 50 >> 51 //Will be properly set in SetDirectParticle() >> 52 IsIon=false; >> 53 massRatio =1.; >> 54 chargeSqRatio=1.; >> 55 preStepChargeSqRatio=1.; >> 56 >> 57 //Some initialization >> 58 currentCoupleIndex=9999999; >> 59 currentCutInRange=0.; >> 60 currentMaterialIndex=9999999; >> 61 currentTcut=0.; >> 62 preStepKinEnergy=0.; >> 63 preStepRange=0.; >> 64 preStepScaledKinEnergy=0.; >> 65 >> 66 currentCouple=0; >> 67 } 52 68 53 ////////////////////////////////////////////// 69 /////////////////////////////////////////////////////// 54 void G4ContinuousGainOfEnergy::ProcessDescript << 70 // >> 71 G4ContinuousGainOfEnergy::~G4ContinuousGainOfEnergy() 55 { 72 { 56 out << "Continuous process acting on adjoint << 73 57 "continuous gain of energy of charged << 74 } 58 "tracked back.\n"; << 75 /////////////////////////////////////////////////////// >> 76 // >> 77 >> 78 void G4ContinuousGainOfEnergy::PreparePhysicsTable( >> 79 const G4ParticleDefinition& ) >> 80 {//theDirectEnergyLossProcess->PreparePhysicsTable(part); >> 81 >> 82 ; 59 } 83 } 60 84 61 ////////////////////////////////////////////// 85 /////////////////////////////////////////////////////// 62 void G4ContinuousGainOfEnergy::SetDirectPartic << 86 // 63 { << 87 64 fDirectPartDef = p; << 88 void G4ContinuousGainOfEnergy::BuildPhysicsTable(const G4ParticleDefinition&) 65 if(fDirectPartDef->GetParticleType() == "nuc << 89 {//theDirectEnergyLossProcess->BuildPhysicsTable(part); 66 { << 90 ; 67 fIsIon = true; << 91 } 68 fMassRatio = proton_mass_c2 / fDirectPartD << 92 69 } << 93 /////////////////////////////////////////////////////// >> 94 // >> 95 void G4ContinuousGainOfEnergy::SetDirectParticle(G4ParticleDefinition* p) >> 96 {theDirectPartDef=p; >> 97 if (theDirectPartDef->GetParticleType()== "nucleus") { >> 98 IsIon=true; >> 99 massRatio = proton_mass_c2/theDirectPartDef->GetPDGMass(); >> 100 G4double q=theDirectPartDef->GetPDGCharge(); >> 101 chargeSqRatio=q*q; >> 102 >> 103 >> 104 } >> 105 70 } 106 } 71 107 72 ////////////////////////////////////////////// 108 /////////////////////////////////////////////////////// >> 109 // >> 110 // 73 G4VParticleChange* G4ContinuousGainOfEnergy::A 111 G4VParticleChange* G4ContinuousGainOfEnergy::AlongStepDoIt(const G4Track& track, 74 << 112 const G4Step& step) 75 { 113 { 76 // Caution in this method the step length sh << 114 77 // A problem is that this is computed by the << 115 //Caution in this method the step length should be the true step length 78 // not know the energy at the end of the adj << 116 // A problem is that this is compute by the multiple scattering that does not know the energy at the end of the adjoint step. This energy is used during the 79 // during the forward sim. Nothing we can re << 117 //Forward sim. Nothing we can really do against that at this time. This is inherent to the MS method 80 // time. This is inherent to the MS method << 118 // 81 << 119 >> 120 >> 121 82 aParticleChange.Initialize(track); 122 aParticleChange.Initialize(track); 83 << 123 84 // Get the actual (true) Step length 124 // Get the actual (true) Step length >> 125 //---------------------------------- 85 G4double length = step.GetStepLength(); 126 G4double length = step.GetStepLength(); 86 G4double degain = 0.0; << 127 G4double degain = 0.0; 87 << 128 >> 129 >> 130 88 // Compute this for weight change after cont 131 // Compute this for weight change after continuous energy loss 89 G4double DEDX_before = << 132 //------------------------------------------------------------- 90 fDirectEnergyLossProcess->GetDEDX(fPreStep << 133 G4double DEDX_before = theDirectEnergyLossProcess->GetDEDX(preStepKinEnergy, currentCouple); 91 << 134 92 // For the fluctuation we generate a new dyn << 135 93 // = preEnergy+egain and then compute the fl << 136 94 // case. << 137 // For the fluctuation we generate a new dynamic particle with energy =preEnergy+egain >> 138 // and then compute the fluctuation given in the direct case. >> 139 //----------------------------------------------------------------------- 95 G4DynamicParticle* dynParticle = new G4Dynam 140 G4DynamicParticle* dynParticle = new G4DynamicParticle(); 96 *dynParticle = *(track.Get << 141 *dynParticle = *(track.GetDynamicParticle()); 97 dynParticle->SetDefinition(fDirectPartDef); << 142 dynParticle->SetDefinition(theDirectPartDef); 98 G4double Tkin = dynParticle->GetKineticEnerg << 143 G4double Tkin = dynParticle->GetKineticEnergy(); 99 << 144 100 G4double dlength = length; << 145 101 if(Tkin != fPreStepKinEnergy && fIsIon) << 146 size_t n=1; 102 { << 147 if (is_integral ) n=10; 103 G4double chargeSqRatio = fCurrentModel->Ge << 148 n=1; 104 fDirectPartDef, fCurrentMaterial, Tkin); << 149 G4double dlength= length/n; 105 fDirectEnergyLossProcess->SetDynamicMassCh << 150 for (size_t i=0;i<n;i++) { 106 } << 151 if (Tkin != preStepKinEnergy && IsIon) { 107 << 152 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin); 108 G4double r = fDirectEnergyLossProcess->GetRa << 153 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio); 109 if(dlength <= fLinLossLimit * r) << 154 110 { << 155 } 111 degain = DEDX_before * dlength; << 156 112 } << 157 G4double r = theDirectEnergyLossProcess->GetRange(Tkin, currentCouple); 113 else << 158 if( dlength <= linLossLimit * r ) { 114 { << 159 degain = DEDX_before*dlength; 115 G4double x = r + dlength; << 160 } 116 G4double E = fDirectEnergyLossProcess->Get << 161 else { 117 if(fIsIon) << 162 G4double x = r + dlength; 118 { << 163 //degain = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple) - theDirectEnergyLossProcess->GetKineticEnergy(r,currentCouple); 119 G4double chargeSqRatio = fCurrentModel-> << 164 G4double E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple); 120 fDirectPartDef, fCurrentMaterial, E); << 165 if (IsIon){ 121 fDirectEnergyLossProcess->SetDynamicMass << 166 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E); 122 G4double x1 = fDirectEnergyLossProcess-> << 167 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio); 123 << 168 G4double x1= theDirectEnergyLossProcess->GetRange(E, currentCouple); 124 G4int ii = 0; << 169 while (std::abs(x-x1)>0.01*x) { 125 constexpr G4int iimax = 100; << 170 E = theDirectEnergyLossProcess->GetKineticEnergy(x,currentCouple); 126 while(std::abs(x - x1) > 0.01 * x) << 171 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,E); 127 { << 172 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio); 128 E = fDirectEnergyLossProcess->GetKinet << 173 x1= theDirectEnergyLossProcess->GetRange(E, currentCouple); 129 chargeSqRatio = fCurrentModel->GetChar << 174 130 fDirectPartDef, fCurrentMaterial, E) << 175 } 131 fDirectEnergyLossProcess->SetDynamicMa << 176 } 132 << 177 133 x1 = fDirectEnergyLossProcess->GetRang << 178 degain=E-Tkin; 134 ++ii; << 179 135 if(ii >= iimax) << 180 136 { << 181 137 break; << 182 } 138 } << 183 //G4cout<<degain<<G4endl; 139 } << 184 G4double tmax = currentModel->MaxSecondaryKinEnergy(dynParticle); 140 } << 185 tmax = std::min(tmax,currentTcut); 141 << 186 142 degain = E - Tkin; << 187 143 } << 188 dynParticle->SetKineticEnergy(Tkin+degain); 144 G4double tmax = fCurrentModel->MaxSecondaryK << 189 145 fCurrentTcut = std::min(fCurrentTcut, tmax); << 190 // Corrections, which cannot be tabulated for ions 146 << 191 //---------------------------------------- 147 dynParticle->SetKineticEnergy(Tkin + degain) << 192 G4double esecdep=0;//not used in most models 148 << 193 currentModel->CorrectionsAlongStep(currentCouple, dynParticle, degain,esecdep, dlength); 149 // Corrections, which cannot be tabulated fo << 194 150 fCurrentModel->CorrectionsAlongStep(fCurrent << 195 // Sample fluctuations 151 << 196 //------------------- 152 // Sample fluctuations << 197 153 G4double deltaE = 0.; << 198 154 if(fLossFluctuationFlag) << 199 G4double deltaE =0.; 155 { << 200 if (lossFluctuationFlag ) { 156 deltaE = fCurrentModel->GetModelOfFluctuat << 201 deltaE = currentModel->GetModelOfFluctuations()-> 157 fCurrentCouple, dynParticle, fCurrentTcu << 202 SampleFluctuations(currentMaterial,dynParticle,tmax,dlength,degain)-degain; 158 - degain; << 203 } 159 } << 204 160 << 205 G4double egain=degain+deltaE; 161 G4double egain = degain + deltaE; << 206 if (egain <=0) egain=degain; 162 if(egain <= 0.) << 207 Tkin+=egain; 163 egain = degain; << 208 dynParticle->SetKineticEnergy(Tkin); 164 Tkin += egain; << 209 } 165 dynParticle->SetKineticEnergy(Tkin); << 210 >> 211 >> 212 166 213 >> 214 167 delete dynParticle; 215 delete dynParticle; 168 << 216 169 if(fIsIon) << 217 if (IsIon){ 170 { << 218 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,Tkin); 171 G4double chargeSqRatio = fCurrentModel->Ge << 219 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatio); 172 fDirectPartDef, fCurrentMaterial, Tkin); << 220 173 fDirectEnergyLossProcess->SetDynamicMassCh << 174 } 221 } 175 << 222 176 G4double DEDX_after = fDirectEnergyLossProce << 223 G4double DEDX_after = theDirectEnergyLossProcess->GetDEDX(Tkin, currentCouple); 177 G4double weight_correction = DEDX_after / DE << 224 178 << 225 >> 226 G4double weight_correction=DEDX_after/DEDX_before; >> 227 >> 228 179 aParticleChange.ProposeEnergy(Tkin); 229 aParticleChange.ProposeEnergy(Tkin); 180 << 230 181 // Caution!!! It is important to select the << 231 //we still need to register in the particleChange the modification of the weight of the particle 182 // as the current weight and not the weight << 232 G4double new_weight=weight_correction*track.GetWeight(); 183 // the track is changed after having applied << 184 << 185 G4double new_weight = << 186 weight_correction * step.GetPostStepPoint( << 187 aParticleChange.SetParentWeightByProcess(fal 233 aParticleChange.SetParentWeightByProcess(false); 188 aParticleChange.ProposeParentWeight(new_weig 234 aParticleChange.ProposeParentWeight(new_weight); >> 235 189 236 190 return &aParticleChange; 237 return &aParticleChange; 191 } << 192 238 >> 239 } 193 ////////////////////////////////////////////// 240 /////////////////////////////////////////////////////// >> 241 // 194 void G4ContinuousGainOfEnergy::SetLossFluctuat 242 void G4ContinuousGainOfEnergy::SetLossFluctuations(G4bool val) 195 { 243 { 196 if(val && !fLossFluctuationArePossible) << 244 if(val && !lossFluctuationArePossible) return; 197 return; << 245 lossFluctuationFlag = val; 198 fLossFluctuationFlag = val; << 199 } 246 } 200 << 201 ////////////////////////////////////////////// 247 /////////////////////////////////////////////////////// 202 G4double G4ContinuousGainOfEnergy::GetContinuo << 248 // 203 << 204 << 205 { << 206 DefineMaterial(track.GetMaterialCutsCouple() << 207 << 208 fPreStepKinEnergy = track.GetKineticEnergy() << 209 fCurrentModel = fDirectEnergyLossProcess << 210 track.GetKineticEnergy() * fMassRatio, fCu << 211 G4double emax_model = fCurrentMode << 212 G4double preStepChargeSqRatio = 0.; << 213 if(fIsIon) << 214 { << 215 G4double chargeSqRatio = fCurrentModel->Ge << 216 fDirectPartDef, fCurrentMaterial, fPreSt << 217 preStepChargeSqRatio = chargeSqRatio; << 218 fDirectEnergyLossProcess->SetDynamicMassCh << 219 << 220 } << 221 << 222 G4double maxE = 1.1 * fPreStepKinEnergy; << 223 << 224 if(fPreStepKinEnergy < fCurrentTcut) << 225 maxE = std::min(fCurrentTcut, maxE); << 226 << 227 maxE = std::min(emax_model * 1.001, maxE); << 228 << 229 G4double preStepRange = << 230 fDirectEnergyLossProcess->GetRange(fPreSte << 231 << 232 if(fIsIon) << 233 { << 234 G4double chargeSqRatioAtEmax = fCurrentMod << 235 fDirectPartDef, fCurrentMaterial, maxE); << 236 fDirectEnergyLossProcess->SetDynamicMassCh << 237 << 238 } << 239 249 240 G4double r1 = fDirectEnergyLossProcess->GetR << 241 250 242 if(fIsIon) << 243 fDirectEnergyLossProcess->SetDynamicMassCh << 244 << 245 251 246 return std::max(r1 - preStepRange, 0.001 * m << 252 G4double G4ContinuousGainOfEnergy::GetContinuousStepLimit(const G4Track& track, >> 253 G4double , G4double , G4double& ) >> 254 { >> 255 G4double x = DBL_MAX; >> 256 x=.1*mm; >> 257 >> 258 >> 259 DefineMaterial(track.GetMaterialCutsCouple()); >> 260 >> 261 preStepKinEnergy = track.GetKineticEnergy(); >> 262 preStepScaledKinEnergy = track.GetKineticEnergy()*massRatio; >> 263 currentModel = theDirectEnergyLossProcess->SelectModelForMaterial(preStepScaledKinEnergy,currentCoupleIndex); >> 264 G4double emax_model=currentModel->HighEnergyLimit(); >> 265 if (IsIon) { >> 266 chargeSqRatio = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,preStepKinEnergy); >> 267 preStepChargeSqRatio = chargeSqRatio; >> 268 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio); >> 269 } >> 270 >> 271 >> 272 G4double maxE =1.1*preStepKinEnergy; >> 273 /*if (preStepKinEnergy< 0.05*MeV) maxE =2.*preStepKinEnergy; >> 274 else if (preStepKinEnergy< 0.1*MeV) maxE =1.5*preStepKinEnergy; >> 275 else if (preStepKinEnergy< 0.5*MeV) maxE =1.25*preStepKinEnergy;*/ >> 276 >> 277 if (preStepKinEnergy < currentTcut) maxE = std::min(currentTcut,maxE); >> 278 >> 279 maxE=std::min(emax_model*1.001,maxE); >> 280 >> 281 preStepRange = theDirectEnergyLossProcess->GetRange(preStepKinEnergy, currentCouple); >> 282 >> 283 if (IsIon) { >> 284 G4double chargeSqRatioAtEmax = currentModel->GetChargeSquareRatio(theDirectPartDef,currentMaterial,maxE); >> 285 theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,chargeSqRatioAtEmax); >> 286 } >> 287 >> 288 G4double r1 = theDirectEnergyLossProcess->GetRange(maxE, currentCouple); >> 289 >> 290 if (IsIon) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,preStepChargeSqRatio); >> 291 >> 292 >> 293 >> 294 x=r1-preStepRange; >> 295 x=std::max(r1-preStepRange,0.001*mm); >> 296 >> 297 return x; >> 298 >> 299 247 } 300 } 248 << 301 #include "G4EmCorrections.hh" 249 ////////////////////////////////////////////// 302 /////////////////////////////////////////////////////// 250 void G4ContinuousGainOfEnergy::SetDynamicMassC << 303 // 251 << 304 252 { << 305 void G4ContinuousGainOfEnergy::SetDynamicMassCharge(const G4Track& ,G4double energy) 253 G4double ChargeSqRatio = << 306 { 254 G4LossTableManager::Instance()->EmCorrecti << 307 255 fDirectPartDef, fCurrentMaterial, energy << 308 G4double ChargeSqRatio= G4LossTableManager::Instance()->EmCorrections()->EffectiveChargeSquareRatio(theDirectPartDef,currentMaterial,energy); 256 if(fDirectEnergyLossProcess) << 309 if (theDirectEnergyLossProcess) theDirectEnergyLossProcess->SetDynamicMassCharge(massRatio,ChargeSqRatio); 257 fDirectEnergyLossProcess->SetDynamicMassCh << 258 } 310 } 259 311