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