| 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 // ------------------------------------------- 27 // -------------------------------------------------------------------
28 // 28 //
29 // GEANT4 Class 29 // GEANT4 Class
30 // File name: G4PAIPhotModel.cc 30 // File name: G4PAIPhotModel.cc
31 // 31 //
32 // Author: Vladimir.Grichine@cern.ch on base o 32 // Author: Vladimir.Grichine@cern.ch on base of G4PAIModel MT interface
33 // 33 //
34 // Creation date: 07.10.2013 34 // Creation date: 07.10.2013
35 // 35 //
36 // Modifications: 36 // Modifications:
37 // 37 //
38 // 38 //
39 39
40 #include "G4PAIPhotModel.hh" 40 #include "G4PAIPhotModel.hh"
41 41
42 #include "G4SystemOfUnits.hh" 42 #include "G4SystemOfUnits.hh"
43 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh"
44 #include "G4Region.hh" 44 #include "G4Region.hh"
>> 45 #include "G4PhysicsLogVector.hh"
>> 46 #include "G4PhysicsFreeVector.hh"
>> 47 #include "G4PhysicsTable.hh"
45 #include "G4ProductionCutsTable.hh" 48 #include "G4ProductionCutsTable.hh"
46 #include "G4MaterialCutsCouple.hh" 49 #include "G4MaterialCutsCouple.hh"
47 #include "G4MaterialTable.hh" 50 #include "G4MaterialTable.hh"
>> 51 #include "G4SandiaTable.hh"
>> 52 #include "G4OrderedTable.hh"
48 #include "G4RegionStore.hh" 53 #include "G4RegionStore.hh"
49 54
50 #include "Randomize.hh" 55 #include "Randomize.hh"
51 #include "G4Electron.hh" 56 #include "G4Electron.hh"
52 #include "G4Positron.hh" 57 #include "G4Positron.hh"
53 #include "G4Gamma.hh" 58 #include "G4Gamma.hh"
54 #include "G4Poisson.hh" 59 #include "G4Poisson.hh"
55 #include "G4Step.hh" 60 #include "G4Step.hh"
56 #include "G4Material.hh" 61 #include "G4Material.hh"
57 #include "G4DynamicParticle.hh" 62 #include "G4DynamicParticle.hh"
58 #include "G4ParticleDefinition.hh" 63 #include "G4ParticleDefinition.hh"
59 #include "G4ParticleChangeForLoss.hh" 64 #include "G4ParticleChangeForLoss.hh"
60 #include "G4PAIPhotData.hh" 65 #include "G4PAIPhotData.hh"
61 #include "G4DeltaAngle.hh" 66 #include "G4DeltaAngle.hh"
62 67
63 ////////////////////////////////////////////// 68 ////////////////////////////////////////////////////////////////////////
64 69
65 using namespace std; 70 using namespace std;
66 71
67 G4PAIPhotModel::G4PAIPhotModel(const G4Particl 72 G4PAIPhotModel::G4PAIPhotModel(const G4ParticleDefinition* p, const G4String& nam)
68 : G4VEmModel(nam),G4VEmFluctuationModel(nam) 73 : G4VEmModel(nam),G4VEmFluctuationModel(nam),
69 fVerbose(0), 74 fVerbose(0),
70 fModelData(nullptr), 75 fModelData(nullptr),
71 fParticle(nullptr) 76 fParticle(nullptr)
72 { 77 {
73 fElectron = G4Electron::Electron(); 78 fElectron = G4Electron::Electron();
74 fPositron = G4Positron::Positron(); 79 fPositron = G4Positron::Positron();
75 80
76 fParticleChange = nullptr; 81 fParticleChange = nullptr;
77 82
78 if(p) { SetParticle(p); } 83 if(p) { SetParticle(p); }
79 else { SetParticle(fElectron); } 84 else { SetParticle(fElectron); }
80 85
81 // default generator 86 // default generator
82 SetAngularDistribution(new G4DeltaAngle()); 87 SetAngularDistribution(new G4DeltaAngle());
83 fLowestTcut = 12.5*CLHEP::eV; 88 fLowestTcut = 12.5*CLHEP::eV;
84 } 89 }
85 90
86 ////////////////////////////////////////////// 91 ////////////////////////////////////////////////////////////////////////////
87 92
88 G4PAIPhotModel::~G4PAIPhotModel() 93 G4PAIPhotModel::~G4PAIPhotModel()
89 { 94 {
>> 95 //G4cout << "G4PAIPhotModel::~G4PAIPhotModel() " << this << G4endl;
90 if(IsMaster()) { delete fModelData; fModelDa 96 if(IsMaster()) { delete fModelData; fModelData = nullptr; }
91 } 97 }
92 98
93 ////////////////////////////////////////////// 99 ////////////////////////////////////////////////////////////////////////////
94 100
95 void G4PAIPhotModel::Initialise(const G4Partic 101 void G4PAIPhotModel::Initialise(const G4ParticleDefinition* p,
96 const G4DataVector& cuts) 102 const G4DataVector& cuts)
97 { 103 {
98 if(fVerbose > 1) << 104 if(fVerbose > 0)
99 { 105 {
100 G4cout<<"G4PAIPhotModel::Initialise for "< 106 G4cout<<"G4PAIPhotModel::Initialise for "<<p->GetParticleName()<<G4endl;
101 } 107 }
102 SetParticle(p); 108 SetParticle(p);
103 fParticleChange = GetParticleChangeForLoss() 109 fParticleChange = GetParticleChangeForLoss();
104 110
105 if( IsMaster() ) 111 if( IsMaster() )
106 { 112 {
>> 113
>> 114 InitialiseElementSelectors(p, cuts);
>> 115
107 delete fModelData; 116 delete fModelData;
108 fMaterialCutsCoupleVector.clear(); 117 fMaterialCutsCoupleVector.clear();
109 118
110 G4double tmin = LowEnergyLimit()*fRatio; 119 G4double tmin = LowEnergyLimit()*fRatio;
111 G4double tmax = HighEnergyLimit()*fRatio; 120 G4double tmax = HighEnergyLimit()*fRatio;
112 fModelData = new G4PAIPhotData(tmin, tmax, 121 fModelData = new G4PAIPhotData(tmin, tmax, fVerbose);
113 122
114 // Prepare initialization 123 // Prepare initialization
115 const G4MaterialTable* theMaterialTable = 124 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
116 size_t numOfMat = G4Material::GetNumberO 125 size_t numOfMat = G4Material::GetNumberOfMaterials();
117 size_t numRegions = fPAIRegionVector.size( 126 size_t numRegions = fPAIRegionVector.size();
118 127
119 // protect for unit tests 128 // protect for unit tests
120 if(0 == numRegions) { 129 if(0 == numRegions) {
121 G4Exception("G4PAIModel::Initialise()"," 130 G4Exception("G4PAIModel::Initialise()","em0106",JustWarning,
122 "no G4Regions are registered 131 "no G4Regions are registered for the PAI model - World is used");
123 fPAIRegionVector.push_back(G4RegionStore 132 fPAIRegionVector.push_back(G4RegionStore::GetInstance()
124 ->GetRegion("DefaultRegionForTheWorld 133 ->GetRegion("DefaultRegionForTheWorld", false));
125 numRegions = 1; 134 numRegions = 1;
126 } 135 }
127 136
128 for( size_t iReg = 0; iReg < numRegions; + 137 for( size_t iReg = 0; iReg < numRegions; ++iReg )
129 { 138 {
130 const G4Region* curReg = fPAIRegionVecto 139 const G4Region* curReg = fPAIRegionVector[iReg];
131 G4Region* reg = const_cast<G4Region*>(cu 140 G4Region* reg = const_cast<G4Region*>(curReg);
132 141
133 for(size_t jMat = 0; jMat < numOfMat; ++ 142 for(size_t jMat = 0; jMat < numOfMat; ++jMat)
134 { 143 {
135 G4Material* mat = (*theMaterialTable)[jMat]; 144 G4Material* mat = (*theMaterialTable)[jMat];
136 const G4MaterialCutsCouple* cutCouple = reg- 145 const G4MaterialCutsCouple* cutCouple = reg->FindCouple(mat);
137 if(nullptr != cutCouple) << 146 //G4cout << "Couple <" << fCutCouple << G4endl;
>> 147 if(cutCouple)
138 { 148 {
139 if(fVerbose > 1) << 149 if(fVerbose>0)
140 { 150 {
141 G4cout << "Reg <" <<curReg->GetName() << 151 G4cout << "Reg <" <<curReg->GetName() << "> mat <"
142 << mat->GetName() << "> fCouple= " 152 << mat->GetName() << "> fCouple= "
143 << cutCouple << ", idx= " << cutCouple-> 153 << cutCouple << ", idx= " << cutCouple->GetIndex()
144 <<" " << p->GetParticleName() 154 <<" " << p->GetParticleName()
145 <<", cuts.size() = " << cuts.size() << G 155 <<", cuts.size() = " << cuts.size() << G4endl;
146 } 156 }
147 // check if this couple is not already ini 157 // check if this couple is not already initialized
148 158
149 size_t n = fMaterialCutsCoupleVector.size( 159 size_t n = fMaterialCutsCoupleVector.size();
150 160
151 G4bool isnew = true; 161 G4bool isnew = true;
152 if( 0 < n ) 162 if( 0 < n )
153 { 163 {
154 for(size_t i=0; i<fMaterialCutsCoupleVec 164 for(size_t i=0; i<fMaterialCutsCoupleVector.size(); ++i)
155 { 165 {
156 if(cutCouple == fMaterialCutsCoupleVec 166 if(cutCouple == fMaterialCutsCoupleVector[i]) {
157 isnew = false; 167 isnew = false;
158 break; 168 break;
159 } 169 }
160 } 170 }
161 } 171 }
162 // initialise data banks 172 // initialise data banks
163 if(isnew) { 173 if(isnew) {
164 fMaterialCutsCoupleVector.push_back(cutC 174 fMaterialCutsCoupleVector.push_back(cutCouple);
165 G4double deltaCutInKinEnergy = cuts[cutC 175 G4double deltaCutInKinEnergy = cuts[cutCouple->GetIndex()];
166 fModelData->Initialise(cutCouple, deltaC 176 fModelData->Initialise(cutCouple, deltaCutInKinEnergy, this);
167 } 177 }
168 } 178 }
169 } 179 }
170 } 180 }
171 InitialiseElementSelectors(p, cuts); <<
172 } 181 }
173 } 182 }
174 183
175 ////////////////////////////////////////////// 184 /////////////////////////////////////////////////////////////////////////
176 185
177 void G4PAIPhotModel::InitialiseLocal(const G4P 186 void G4PAIPhotModel::InitialiseLocal(const G4ParticleDefinition*,
178 G4VEmModel* masterModel) 187 G4VEmModel* masterModel)
179 { 188 {
180 fModelData = static_cast<G4PAIPhotModel*>(ma 189 fModelData = static_cast<G4PAIPhotModel*>(masterModel)->GetPAIPhotData();
181 fMaterialCutsCoupleVector = static_cast<G4PA 190 fMaterialCutsCoupleVector = static_cast<G4PAIPhotModel*>(masterModel)->GetVectorOfCouples();
182 SetElementSelectors( masterModel->GetElement 191 SetElementSelectors( masterModel->GetElementSelectors() );
183 } 192 }
184 193
185 ////////////////////////////////////////////// 194 //////////////////////////////////////////////////////////////////////////////
186 195
187 G4double G4PAIPhotModel::MinEnergyCut(const G4 196 G4double G4PAIPhotModel::MinEnergyCut(const G4ParticleDefinition*,
188 const G4MaterialCutsCouple*) 197 const G4MaterialCutsCouple*)
189 { 198 {
190 return fLowestTcut; 199 return fLowestTcut;
191 } 200 }
192 201
193 ////////////////////////////////////////////// 202 //////////////////////////////////////////////////////////////////////////////
194 203
195 G4double G4PAIPhotModel::ComputeDEDXPerVolume( 204 G4double G4PAIPhotModel::ComputeDEDXPerVolume(const G4Material*,
196 const G4ParticleDefinition* p, 205 const G4ParticleDefinition* p,
197 G4double kineticEnergy, 206 G4double kineticEnergy,
198 G4double cutEnergy) 207 G4double cutEnergy)
199 { 208 {
200 G4int coupleIndex = FindCoupleIndex(CurrentC 209 G4int coupleIndex = FindCoupleIndex(CurrentCouple());
201 if(0 > coupleIndex) { return 0.0; } 210 if(0 > coupleIndex) { return 0.0; }
202 211
203 G4double cut = std::min(MaxSecondaryEnergy(p 212 G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy);
>> 213
204 G4double scaledTkin = kineticEnergy*fRatio; 214 G4double scaledTkin = kineticEnergy*fRatio;
205 G4double dedx = fChargeSquare*fModelData->DE << 215
206 return dedx; << 216 return fChargeSquare*fModelData->DEDXPerVolume(coupleIndex, scaledTkin, cut);
207 } 217 }
208 218
209 ////////////////////////////////////////////// 219 /////////////////////////////////////////////////////////////////////////
210 220
211 G4double G4PAIPhotModel::CrossSectionPerVolume 221 G4double G4PAIPhotModel::CrossSectionPerVolume( const G4Material*,
212 const G4ParticleDefinition* p, 222 const G4ParticleDefinition* p,
213 G4double kineticEnergy, 223 G4double kineticEnergy,
214 G4double cutEnergy, 224 G4double cutEnergy,
215 G4double maxEnergy ) 225 G4double maxEnergy )
216 { 226 {
217 G4int coupleIndex = FindCoupleIndex(CurrentC 227 G4int coupleIndex = FindCoupleIndex(CurrentCouple());
218 if(0 > coupleIndex) { return 0.0; } << 228
>> 229 if(0 > coupleIndex) return 0.0;
219 230
220 G4double tmax = std::min(MaxSecondaryEnergy( 231 G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy);
221 if(tmax <= cutEnergy) { return 0.0; } << 232
>> 233 if(tmax <= cutEnergy) return 0.0;
222 234
223 G4double scaledTkin = kineticEnergy*fRatio; 235 G4double scaledTkin = kineticEnergy*fRatio;
224 G4double xs = fChargeSquare*fModelData->Cros << 236 G4double xsc = fChargeSquare*fModelData->CrossSectionPerVolume(coupleIndex,
225 scal << 237 scaledTkin,
226 return xs; << 238 cutEnergy, tmax);
>> 239 return xsc;
227 } 240 }
228 241
229 ////////////////////////////////////////////// 242 ///////////////////////////////////////////////////////////////////////////
230 // 243 //
231 // It is analog of PostStepDoIt in terms of se 244 // It is analog of PostStepDoIt in terms of secondary electron.
232 // 245 //
233 246
234 void G4PAIPhotModel::SampleSecondaries(std::ve 247 void G4PAIPhotModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp,
235 const G4MaterialCutsCouple* matCC, 248 const G4MaterialCutsCouple* matCC,
236 const G4DynamicParticle* dp, 249 const G4DynamicParticle* dp,
237 G4double tmin, 250 G4double tmin,
238 G4double maxEnergy) 251 G4double maxEnergy)
239 { 252 {
240 G4int coupleIndex = FindCoupleIndex(matCC); 253 G4int coupleIndex = FindCoupleIndex(matCC);
241 if(0 > coupleIndex) { return; } 254 if(0 > coupleIndex) { return; }
242 255
243 SetParticle(dp->GetDefinition()); 256 SetParticle(dp->GetDefinition());
244 257
245 G4double kineticEnergy = dp->GetKineticEnerg 258 G4double kineticEnergy = dp->GetKineticEnergy();
246 259
247 G4double tmax = MaxSecondaryEnergy(fParticle 260 G4double tmax = MaxSecondaryEnergy(fParticle, kineticEnergy);
248 261
249 if( maxEnergy < tmax) tmax = maxEnergy; 262 if( maxEnergy < tmax) tmax = maxEnergy;
250 if( tmin >= tmax) return; 263 if( tmin >= tmax) return;
251 264
252 G4ThreeVector direction = dp->GetMomentumDir 265 G4ThreeVector direction = dp->GetMomentumDirection();
253 G4double scaledTkin = kineticEnergy*fRat 266 G4double scaledTkin = kineticEnergy*fRatio;
254 G4double totalEnergy = kineticEnergy + fM 267 G4double totalEnergy = kineticEnergy + fMass;
255 G4double totalMomentum = sqrt(kineticEnergy 268 G4double totalMomentum = sqrt(kineticEnergy*(totalEnergy + fMass));
256 G4double plRatio = fModelData->GetPla 269 G4double plRatio = fModelData->GetPlasmonRatio(coupleIndex, scaledTkin);
257 270
258 if( G4UniformRand() <= plRatio ) 271 if( G4UniformRand() <= plRatio )
259 { 272 {
260 G4double deltaTkin = fModelData->SamplePos 273 G4double deltaTkin = fModelData->SamplePostStepPlasmonTransfer(coupleIndex, scaledTkin);
261 274
262 // G4cout<<"G4PAIPhotModel::SampleSecondar 275 // G4cout<<"G4PAIPhotModel::SampleSecondaries; dp "<<dp->GetParticleDefinition()->GetParticleName()
263 // <<"; Tkin = "<<kineticEnergy/keV<<" keV 276 // <<"; Tkin = "<<kineticEnergy/keV<<" keV; transfer = "<<deltaTkin/keV<<" keV "<<G4endl;
264 277
265 if( deltaTkin <= 0. && fVerbose > 0) 278 if( deltaTkin <= 0. && fVerbose > 0)
266 { 279 {
267 G4cout<<"G4PAIPhotModel::SampleSecondary << 280 G4cout<<"G4PAIPhotModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl;
268 } 281 }
269 if( deltaTkin <= 0.) { return; } << 282 if( deltaTkin <= 0.) return;
270 283
271 if( deltaTkin > tmax) { deltaTkin = tmax; << 284 if( deltaTkin > tmax) deltaTkin = tmax;
272 285
273 const G4Element* anElement = SelectTargetA << 286 const G4Element* anElement = SelectRandomAtom(matCC,fParticle,kineticEnergy);
274 << 287 G4int Z = G4lrint(anElement->GetZ());
275 G4int Z = anElement->GetZasInt(); <<
276 288
277 auto deltaRay = new G4DynamicParticle(fEle << 289 G4DynamicParticle* deltaRay = new G4DynamicParticle(fElectron,
278 GetAngularDistribution()->SampleDirectio 290 GetAngularDistribution()->SampleDirection(dp, deltaTkin,
279 Z, matCC->GetMaterial()), 291 Z, matCC->GetMaterial()),
280 deltaTkin); 292 deltaTkin);
281 293
282 // primary change 294 // primary change
283 295
284 kineticEnergy -= deltaTkin; 296 kineticEnergy -= deltaTkin;
285 297
286 if( kineticEnergy <= 0. ) // kill primary 298 if( kineticEnergy <= 0. ) // kill primary as local? energy deposition
287 { 299 {
288 fParticleChange->SetProposedKineticEnerg 300 fParticleChange->SetProposedKineticEnergy(0.0);
289 fParticleChange->ProposeLocalEnergyDepos 301 fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy+deltaTkin);
>> 302 // fParticleChange->ProposeTrackStatus(fStopAndKill);
290 return; 303 return;
291 } 304 }
292 else 305 else
293 { 306 {
294 G4ThreeVector dir = totalMomentum*direct 307 G4ThreeVector dir = totalMomentum*direction - deltaRay->GetMomentum();
295 direction = dir.unit(); 308 direction = dir.unit();
296 fParticleChange->SetProposedKineticEnerg 309 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
297 fParticleChange->SetProposedMomentumDire 310 fParticleChange->SetProposedMomentumDirection(direction);
298 vdp->push_back(deltaRay); 311 vdp->push_back(deltaRay);
299 } 312 }
300 } 313 }
301 else // secondary X-ray CR photon 314 else // secondary X-ray CR photon
302 { 315 {
303 G4double deltaTkin = fModelData->Sampl 316 G4double deltaTkin = fModelData->SamplePostStepPhotonTransfer(coupleIndex, scaledTkin);
304 317
305 // G4cout<<"PAIPhotonModel PhotonPostStep 318 // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl;
306 319
307 if( deltaTkin <= 0. ) 320 if( deltaTkin <= 0. )
308 { 321 {
309 G4cout<<"G4PAIPhotonModel::SampleSeconda 322 G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl;
310 } 323 }
311 if( deltaTkin <= 0.) return; 324 if( deltaTkin <= 0.) return;
312 325
313 if( deltaTkin >= kineticEnergy ) // stop p 326 if( deltaTkin >= kineticEnergy ) // stop primary
314 { 327 {
315 deltaTkin = kineticEnergy; 328 deltaTkin = kineticEnergy;
316 kineticEnergy = 0.0; 329 kineticEnergy = 0.0;
317 } 330 }
318 G4double costheta = 0.; // G4UniformRand() 331 G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only
319 G4double sintheta = sqrt((1.+costheta)*(1. 332 G4double sintheta = sqrt((1.+costheta)*(1.-costheta));
320 333
321 // direction of the 'Cherenkov' photon 334 // direction of the 'Cherenkov' photon
322 G4double phi = twopi*G4UniformRand(); 335 G4double phi = twopi*G4UniformRand();
323 G4double dirx = sintheta*cos(phi), diry = 336 G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta;
324 337
325 G4ThreeVector deltaDirection(dirx,diry,dir 338 G4ThreeVector deltaDirection(dirx,diry,dirz);
326 deltaDirection.rotateUz(direction); 339 deltaDirection.rotateUz(direction);
327 340
328 if( kineticEnergy > 0.) // primary change 341 if( kineticEnergy > 0.) // primary change
329 { 342 {
330 kineticEnergy -= deltaTkin; 343 kineticEnergy -= deltaTkin;
331 fParticleChange->SetProposedKineticEnerg 344 fParticleChange->SetProposedKineticEnergy(kineticEnergy);
332 } 345 }
333 else // stop primary, but pass X-ray CR 346 else // stop primary, but pass X-ray CR
334 { 347 {
335 // fParticleChange->ProposeLocalEnergyDe 348 // fParticleChange->ProposeLocalEnergyDeposit(deltaTkin);
336 fParticleChange->SetProposedKineticEnerg 349 fParticleChange->SetProposedKineticEnergy(0.0);
337 } 350 }
338 // create G4DynamicParticle object for pho 351 // create G4DynamicParticle object for photon ray
339 352
340 auto photonRay = new G4DynamicParticle; << 353 G4DynamicParticle* photonRay = new G4DynamicParticle;
341 photonRay->SetDefinition( G4Gamma::Gamma() 354 photonRay->SetDefinition( G4Gamma::Gamma() );
342 photonRay->SetKineticEnergy( deltaTkin ); 355 photonRay->SetKineticEnergy( deltaTkin );
343 photonRay->SetMomentumDirection(deltaDirec 356 photonRay->SetMomentumDirection(deltaDirection);
344 357
345 vdp->push_back(photonRay); 358 vdp->push_back(photonRay);
346 } 359 }
347 return; 360 return;
348 } 361 }
349 362
350 ////////////////////////////////////////////// 363 ///////////////////////////////////////////////////////////////////////
351 364
352 G4double G4PAIPhotModel::SampleFluctuations( << 365 G4double G4PAIPhotModel::SampleFluctuations( const G4MaterialCutsCouple* matCC,
353 const G4MaterialCutsC << 366 const G4DynamicParticle* aParticle,
354 const G4DynamicPartic << 367 G4double, G4double step,
355 const G4double, const << 368 G4double eloss)
356 const G4double step, <<
357 { 369 {
358 // return 0.; 370 // return 0.;
359 G4int coupleIndex = FindCoupleIndex(matCC); 371 G4int coupleIndex = FindCoupleIndex(matCC);
360 if(0 > coupleIndex) { return eloss; } 372 if(0 > coupleIndex) { return eloss; }
361 373
362 SetParticle(aParticle->GetDefinition()); 374 SetParticle(aParticle->GetDefinition());
363 375
>> 376
364 // G4cout << "G4PAIPhotModel::SampleFluctuat 377 // G4cout << "G4PAIPhotModel::SampleFluctuations step(mm)= "<< step/mm
365 // << " Eloss(keV)= " << eloss/keV << " in 378 // << " Eloss(keV)= " << eloss/keV << " in "
366 // << matCC->GetMaterial()->GetName() << G4e 379 // << matCC->GetMaterial()->GetName() << G4endl;
>> 380
367 381
368 G4double Tkin = aParticle->GetKineticE 382 G4double Tkin = aParticle->GetKineticEnergy();
369 G4double scaledTkin = Tkin*fRatio; 383 G4double scaledTkin = Tkin*fRatio;
370 384
371 G4double loss = fModelData->SampleAlongStepP 385 G4double loss = fModelData->SampleAlongStepPhotonTransfer(coupleIndex, Tkin,
372 scaledTkin, << 386 scaledTkin,
373 step*fChargeSquare); << 387 step*fChargeSquare);
374 loss += fModelData->SampleAlongStepPlasmonTr << 388 loss += fModelData->SampleAlongStepPlasmonTransfer(coupleIndex, Tkin,
375 << 389 scaledTkin,
>> 390 step*fChargeSquare);
>> 391
376 392
377 // G4cout<<" PAIPhotModel::SampleFluctuatio 393 // G4cout<<" PAIPhotModel::SampleFluctuations loss = "<<loss/keV<<" keV, on step = "
378 // <<step/mm<<" mm"<<G4endl; 394 // <<step/mm<<" mm"<<G4endl;
379 return loss; 395 return loss;
380 396
381 } 397 }
382 398
383 ////////////////////////////////////////////// 399 //////////////////////////////////////////////////////////////////////
384 // 400 //
385 // Returns the statistical estimation of the e 401 // Returns the statistical estimation of the energy loss distribution variance
386 // 402 //
387 403
388 404
389 G4double G4PAIPhotModel::Dispersion(const G4Ma << 405 G4double G4PAIPhotModel::Dispersion( const G4Material* material,
390 const G4Dy << 406 const G4DynamicParticle* aParticle,
391 const G4double tcut, << 407 G4double tmax,
392 const G4double tmax, << 408 G4double step )
393 const G4double step) <<
394 { 409 {
395 G4double particleMass = aParticle->GetMass( 410 G4double particleMass = aParticle->GetMass();
396 G4double electronDensity = material->GetElec 411 G4double electronDensity = material->GetElectronDensity();
397 G4double kineticEnergy = aParticle->GetKinet 412 G4double kineticEnergy = aParticle->GetKineticEnergy();
398 G4double q = aParticle->GetCharge()/eplus; 413 G4double q = aParticle->GetCharge()/eplus;
399 G4double etot = kineticEnergy + particleMass 414 G4double etot = kineticEnergy + particleMass;
400 G4double beta2 = kineticEnergy*(kineticEnerg 415 G4double beta2 = kineticEnergy*(kineticEnergy + 2.0*particleMass)/(etot*etot);
401 G4double siga = (tmax/beta2 - 0.5*tcut) * t << 416 G4double siga = (1.0/beta2 - 0.5) * twopi_mc2_rcl2 * tmax * step
402 * electronDensity * q * q; 417 * electronDensity * q * q;
403 418
404 return siga; 419 return siga;
>> 420 /*
>> 421 G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.;
>> 422 for(G4int i = 0; i < fMeanNumber; i++)
>> 423 {
>> 424 loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss);
>> 425 sumLoss += loss;
>> 426 sumLoss2 += loss*loss;
>> 427 }
>> 428 meanLoss = sumLoss/fMeanNumber;
>> 429 sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber;
>> 430 return sigma2;
>> 431 */
405 } 432 }
406 433
407 ////////////////////////////////////////////// 434 /////////////////////////////////////////////////////////////////////
408 435
409 G4double G4PAIPhotModel::MaxSecondaryEnergy( c 436 G4double G4PAIPhotModel::MaxSecondaryEnergy( const G4ParticleDefinition* p,
410 G4double kinEnergy) 437 G4double kinEnergy)
411 { 438 {
412 SetParticle(p); 439 SetParticle(p);
413 G4double tmax = kinEnergy; 440 G4double tmax = kinEnergy;
414 if(p == fElectron) { tmax *= 0.5; } 441 if(p == fElectron) { tmax *= 0.5; }
415 else if(p != fPositron) { 442 else if(p != fPositron) {
416 G4double ratio= electron_mass_c2/fMass; 443 G4double ratio= electron_mass_c2/fMass;
417 G4double gamma= kinEnergy/fMass + 1.0; 444 G4double gamma= kinEnergy/fMass + 1.0;
418 tmax = 2.0*electron_mass_c2*(gamma*gamma - 445 tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) /
419 (1. + 2.0*gamma*ratio + rati 446 (1. + 2.0*gamma*ratio + ratio*ratio);
420 } 447 }
421 return tmax; 448 return tmax;
422 } 449 }
423 450
424 ////////////////////////////////////////////// 451 ///////////////////////////////////////////////////////////////
425 452
426 void G4PAIPhotModel::DefineForRegion(const G4R 453 void G4PAIPhotModel::DefineForRegion(const G4Region* r)
427 { 454 {
428 fPAIRegionVector.push_back(r); 455 fPAIRegionVector.push_back(r);
429 } 456 }
430 457
431 // 458 //
432 // 459 //
433 ////////////////////////////////////////////// 460 /////////////////////////////////////////////////
434 461