source: trunk/source/processes/electromagnetic/lowenergy/src/G4PenelopeCompton.cc@ 830

Last change on this file since 830 was 819, checked in by garnier, 17 years ago

import all except CVS
File size: 23.3 KB
RevLine 
[819]1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26// $Id: G4PenelopeCompton.cc,v 1.26 2006/06/29 19:40:41 gunter Exp $
27// GEANT4 tag $Name: geant4-09-01-patch-02 $
28//
29// Author: Luciano Pandola
30//
31// History:
32// --------
33// 12 Feb 2003 MG Pia const argument in SelectRandomAtomForCompton
34// Migration to "cuts per region"
35// 14 Feb 2003 MG Pia Corrected compilation errors and warnings
36// from SUN
37// Modified some variables to lowercase initial
38// 10 Mar 2003 V.Ivanchenko Remove CutPerMaterial warning
39// 13 Mar 2003 L.Pandola Code "cleaned"
40// 20 Mar 2003 L.Pandola ReadData() changed (performance improved)
41// 26 Mar 2003 L.Pandola Added fluorescence
42// 24 May 2003 MGP Removed memory leak
43// 09 Mar 2004 L.Pandola Bug fixed in the generation of final state
44// (bug report # 585)
45// 17 Mar 2004 L.Pandola Removed unnecessary calls to std::pow(a,b)
46// 18 Mar 2004 L.Pandola Use of std::map (code review)
47//
48// -------------------------------------------------------------------
49
50#include "G4PenelopeCompton.hh"
51#include "Randomize.hh"
52#include "G4ParticleDefinition.hh"
53#include "G4Track.hh"
54#include "G4Step.hh"
55#include "G4ForceCondition.hh"
56#include "G4Gamma.hh"
57#include "G4Electron.hh"
58#include "G4DynamicParticle.hh"
59#include "G4VParticleChange.hh"
60#include "G4ThreeVector.hh"
61#include "G4EnergyLossTables.hh"
62#include "G4VCrossSectionHandler.hh"
63#include "G4CrossSectionHandler.hh"
64#include "G4VEMDataSet.hh"
65#include "G4EMDataSet.hh"
66#include "G4CompositeEMDataSet.hh"
67#include "G4VDataSetAlgorithm.hh"
68#include "G4LogLogInterpolation.hh"
69#include "G4VRangeTest.hh"
70#include "G4RangeTest.hh"
71#include "G4ProductionCutsTable.hh"
72#include "G4AtomicTransitionManager.hh"
73#include "G4AtomicShell.hh"
74#include "G4AtomicDeexcitation.hh"
75#include "G4PenelopeIntegrator.hh"
76#include "G4MaterialCutsCouple.hh"
77
78
79G4PenelopeCompton::G4PenelopeCompton(const G4String& processName)
80 : G4VDiscreteProcess(processName),
81 lowEnergyLimit(250*eV),
82 highEnergyLimit(100*GeV),
83 intrinsicLowEnergyLimit(10*eV),
84 intrinsicHighEnergyLimit(100*GeV),
85 energyForIntegration(0.0),
86 ZForIntegration(1),
87 nBins(200),
88 cutForLowEnergySecondaryPhotons(250.0*eV)
89{
90 if (lowEnergyLimit < intrinsicLowEnergyLimit ||
91 highEnergyLimit > intrinsicHighEnergyLimit)
92 {
93 G4Exception("G4PenelopeCompton::G4PenelopeCompton - energy outside intrinsic process validity range");
94 }
95
96 meanFreePathTable = 0;
97 ionizationEnergy = new std::map<G4int,G4DataVector*>;
98 hartreeFunction = new std::map<G4int,G4DataVector*>;
99 occupationNumber = new std::map<G4int,G4DataVector*>;
100
101 rangeTest = new G4RangeTest;
102
103 ReadData(); //Read data from file
104
105 if (verboseLevel > 0)
106 {
107 G4cout << GetProcessName() << " is created " << G4endl
108 << "Energy range: "
109 << lowEnergyLimit / keV << " keV - "
110 << highEnergyLimit / GeV << " GeV"
111 << G4endl;
112 }
113}
114
115G4PenelopeCompton::~G4PenelopeCompton()
116{
117 delete meanFreePathTable;
118 delete rangeTest;
119
120 for (size_t i1=0;i1<matCrossSections->size();i1++)
121 {
122 delete (*matCrossSections)[i1];
123 }
124
125 delete matCrossSections;
126
127 for (G4int Z=1;Z<100;Z++)
128 {
129 if (ionizationEnergy->count(Z)) delete (ionizationEnergy->find(Z)->second);
130 if (hartreeFunction->count(Z)) delete (hartreeFunction->find(Z)->second);
131 if (occupationNumber->count(Z)) delete (occupationNumber->find(Z)->second);
132 }
133 delete ionizationEnergy;
134 delete hartreeFunction;
135 delete occupationNumber;
136}
137
138void G4PenelopeCompton::BuildPhysicsTable(const G4ParticleDefinition& )
139{
140 G4DataVector energyVector;
141 G4double dBin = std::log10(highEnergyLimit/lowEnergyLimit)/nBins;
142 G4int i;
143 for (i=0;i<nBins;i++)
144 {
145 energyVector.push_back(std::pow(10.,std::log10(lowEnergyLimit)+i*dBin));
146 }
147
148 const G4MaterialTable* materialTable = G4Material::GetMaterialTable();
149 G4int nMaterials = G4Material::GetNumberOfMaterials();
150 G4VDataSetAlgorithm* algo = new G4LogLogInterpolation();
151
152 size_t nOfBins = energyVector.size();
153 size_t bin=0;
154
155 G4DataVector* energies;
156 G4DataVector* data;
157
158 matCrossSections = new std::vector<G4VEMDataSet*>;
159
160
161 G4int m;
162 for (m=0; m<nMaterials; m++)
163 {
164 const G4Material* material= (*materialTable)[m];
165 G4int nElements = material->GetNumberOfElements();
166 const G4ElementVector* elementVector = material->GetElementVector();
167 const G4double* nAtomsPerVolume = material->GetAtomicNumDensityVector();
168
169 G4VEMDataSet* setForMat = new G4CompositeEMDataSet(algo,1.,1.);
170
171 for (i=0; i<nElements; i++) {
172
173 G4int Z = (G4int) (*elementVector)[i]->GetZ();
174 G4double density = nAtomsPerVolume[i];
175 G4double cross=0.0;
176 energies = new G4DataVector;
177 data = new G4DataVector;
178
179
180 for (bin=0; bin<nOfBins; bin++)
181 {
182 G4double e = energyVector[bin];
183 energies->push_back(e);
184 cross = density * CrossSection(e,Z);
185 data->push_back(cross);
186 }
187
188 G4VEMDataSet* elSet = new G4EMDataSet(i,energies,data,algo,1.,1.);
189 setForMat->AddComponent(elSet);
190 }
191
192 matCrossSections->push_back(setForMat);
193 }
194
195
196 //Build the mean free path table!
197 G4double matCS = 0.0;
198 G4VEMDataSet* matCrossSet = new G4CompositeEMDataSet(algo,1.,1.);
199 G4VEMDataSet* materialSet = new G4CompositeEMDataSet(algo,1.,1.);
200
201
202 for (m=0; m<nMaterials; m++)
203 {
204 energies = new G4DataVector;
205 data = new G4DataVector;
206 const G4Material* material= (*materialTable)[m];
207 material= (*materialTable)[m];
208 for (bin=0; bin<nOfBins; bin++)
209 {
210 G4double energy = energyVector[bin];
211 energies->push_back(energy);
212 matCrossSet = (*matCrossSections)[m];
213 matCS = 0.0;
214 G4int nElm = matCrossSet->NumberOfComponents();
215 for(G4int j=0; j<nElm; j++) {
216 matCS += matCrossSet->GetComponent(j)->FindValue(energy);
217 }
218 if (matCS > 0.)
219 {
220 data->push_back(1./matCS);
221 }
222 else
223 {
224 data->push_back(DBL_MAX);
225 }
226 }
227 G4VEMDataSet* dataSet = new G4EMDataSet(m,energies,data,algo,1.,1.);
228 materialSet->AddComponent(dataSet);
229 }
230 meanFreePathTable = materialSet;
231}
232
233G4VParticleChange* G4PenelopeCompton::PostStepDoIt(const G4Track& aTrack,
234 const G4Step& aStep)
235{
236 //Penelope model
237
238 aParticleChange.Initialize(aTrack);
239
240 // Dynamic particle quantities
241 const G4DynamicParticle* incidentPhoton = aTrack.GetDynamicParticle();
242 G4double photonEnergy0 = incidentPhoton->GetKineticEnergy();
243
244 if (photonEnergy0 <= lowEnergyLimit)
245 {
246 aParticleChange.ProposeTrackStatus(fStopAndKill);
247 aParticleChange.ProposeEnergy(0.);
248 aParticleChange.ProposeLocalEnergyDeposit(photonEnergy0);
249 return G4VDiscreteProcess::PostStepDoIt(aTrack,aStep);
250 }
251
252 G4ParticleMomentum photonDirection0 = incidentPhoton->GetMomentumDirection();
253
254 const G4MaterialCutsCouple* couple = aTrack.GetMaterialCutsCouple();
255 const G4Material* material = couple->GetMaterial();
256 G4int Z = SelectRandomAtomForCompton(material,photonEnergy0);
257 const G4int nmax = 64;
258 G4double rn[nmax],pac[nmax];
259
260 G4double ki,ki1,ki2,ki3,taumin,a1,a2;
261 G4double tau,TST;
262 G4double S=0.0;
263 G4double epsilon,cosTheta;
264 G4double harFunc = 0.0;
265 G4int occupNb= 0;
266 G4double ionEnergy=0.0;
267 G4int nosc = occupationNumber->find(Z)->second->size();
268 G4int iosc = nosc;
269 ki = photonEnergy0/electron_mass_c2;
270 ki2 = 2*ki+1.0;
271 ki3 = ki*ki;
272 ki1 = ki3-ki2-1.0;
273 taumin = 1.0/ki2;
274 a1 = std::log(ki2);
275 a2 = a1+2.0*ki*(1.0+ki)/(ki2*ki2);
276 if (photonEnergy0 > 5*MeV)
277 {
278 do{
279 do{
280 if ((a2*G4UniformRand()) < a1)
281 {
282 tau = std::pow(taumin,G4UniformRand());
283 }
284 else
285 {
286 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
287 }
288 //rejection function
289 TST = (1+tau*(ki1+tau*(ki2+tau*ki3)))/(ki3*tau*(1.0+tau*tau));
290 }while (G4UniformRand()> TST);
291 epsilon=tau;
292 cosTheta = 1.0 - (1.0-tau)/(ki*tau);
293 //Target shell electrons
294 TST = Z*G4UniformRand();
295 iosc = nosc;
296 S=0.0;
297 for (G4int j=0;j<nosc;j++)
298 {
299 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[j];
300 S = S + occupNb;
301 if (S > TST) iosc = j;
302 if (S > TST) break;
303 }
304 ionEnergy = (*(ionizationEnergy->find(Z)->second))[iosc];
305 }while((epsilon*photonEnergy0-photonEnergy0+ionEnergy) >0);
306 }
307
308 else //photonEnergy0<5 MeV
309 {
310 //Incoherent scattering function for theta=PI
311 G4double s0=0.0;
312 G4double pzomc=0.0,rni=0.0;
313 G4double aux=0.0;
314 for (G4int i=0;i<nosc;i++){
315 ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
316 if (photonEnergy0 > ionEnergy)
317 {
318 G4double aux = photonEnergy0*(photonEnergy0-ionEnergy)*2.0;
319 harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
320 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
321 pzomc = harFunc*(aux-electron_mass_c2*ionEnergy)/
322 (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
323 if (pzomc > 0)
324 {
325 rni = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*(std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
326 }
327 else
328 {
329 rni = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*(std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
330 }
331 s0 = s0 + occupNb*rni;
332 }
333 }
334
335 //Sampling tau
336 G4double cdt1;
337 do
338 {
339 if ((G4UniformRand()*a2) < a1)
340 {
341 tau = std::pow(taumin,G4UniformRand());
342 }
343 else
344 {
345 tau = std::sqrt(1.0+G4UniformRand()*(taumin*taumin-1.0));
346 }
347 cdt1 = (1.0-tau)/(ki*tau);
348 S=0.0;
349 //Incoherent scattering function
350 for (G4int i=0;i<nosc;i++){
351 ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
352 if (photonEnergy0 > ionEnergy) //sum only on excitable levels
353 {
354 aux = photonEnergy0*(photonEnergy0-ionEnergy)*cdt1;
355 harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
356 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
357 pzomc = harFunc*(aux-electron_mass_c2*ionEnergy)/
358 (electron_mass_c2*std::sqrt(2.0*aux+ionEnergy*ionEnergy));
359 if (pzomc > 0)
360 {
361 rn[i] = 1.0-0.5*std::exp(0.5-(std::sqrt(0.5)+std::sqrt(2.0)*pzomc)*(std::sqrt(0.5)+std::sqrt(2.0)*pzomc));
362 }
363 else
364 {
365 rn[i] = 0.5*std::exp(0.5-(std::sqrt(0.5)-std::sqrt(2.0)*pzomc)*(std::sqrt(0.5)-std::sqrt(2.0)*pzomc));
366 }
367 S = S + occupNb*rn[i];
368 pac[i] = S;
369 }
370 else
371 {
372 pac[i] = S-(1e-06);
373 }
374 }
375 //Rejection function
376 TST = S*(1.0+tau*(ki1+tau*(ki2+tau*ki3)))/(ki3*tau*(1.0+tau*tau));
377 }while ((G4UniformRand()*s0) > TST);
378
379 //Target electron shell
380 cosTheta = 1.0 - cdt1;
381 G4double fpzmax=0.0,fpz=0.0;
382 G4double A=0.0;
383 do
384 {
385 do
386 {
387 TST =S*G4UniformRand();
388 iosc=nosc;
389 for (G4int i=0;i<nosc;i++){
390 if (pac[i]>TST) iosc = i;
391 if (pac[i]>TST) break;
392 }
393 A = G4UniformRand()*rn[iosc];
394 harFunc = (*(hartreeFunction->find(Z)->second))[iosc]/fine_structure_const;
395 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[iosc];
396 if (A < 0.5) {
397 pzomc = (std::sqrt(0.5)-std::sqrt(0.5-std::log(2.0*A)))/
398 (std::sqrt(2.0)*harFunc);
399 }
400 else
401 {
402 pzomc = (std::sqrt(0.5-std::log(2.0-2.0*A))-std::sqrt(0.5))/
403 (std::sqrt(2.0)*harFunc);
404 }
405 } while (pzomc < -1);
406 // F(EP) rejection
407 G4double XQC = 1.0+tau*(tau-2.0*cosTheta);
408 G4double AF = std::sqrt(XQC)*(1.0+tau*(tau-cosTheta)/XQC);
409 if (AF > 0) {
410 fpzmax = 1.0+AF*0.2;
411 }
412 else
413 {
414 fpzmax = 1.0-AF*0.2;
415 }
416 fpz = 1.0+AF*std::max(std::min(pzomc,0.2),-0.2);
417 }while ((fpzmax*G4UniformRand())>fpz);
418
419 //Energy of the scattered photon
420 G4double T = pzomc*pzomc;
421 G4double b1 = 1.0-T*tau*tau;
422 G4double b2 = 1.0-T*tau*cosTheta;
423 if (pzomc > 0.0)
424 {
425 epsilon = (tau/b1)*(b2+std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
426 }
427 else
428 {
429 epsilon = (tau/b1)*(b2-std::sqrt(std::abs(b2*b2-b1*(1.0-T))));
430 }
431 }
432
433
434 G4double sinTheta = std::sqrt(1-cosTheta*cosTheta);
435 G4double phi = twopi * G4UniformRand() ;
436 G4double dirx = sinTheta * std::cos(phi);
437 G4double diry = sinTheta * std::sin(phi);
438 G4double dirz = cosTheta ;
439
440 // Update G4VParticleChange for the scattered photon
441
442 G4ThreeVector photonDirection1(dirx,diry,dirz);
443 photonDirection1.rotateUz(photonDirection0);
444 aParticleChange.ProposeMomentumDirection(photonDirection1) ;
445 G4double photonEnergy1 = epsilon * photonEnergy0;
446
447 if (photonEnergy1 > 0.)
448 {
449 aParticleChange.ProposeEnergy(photonEnergy1) ;
450 }
451 else
452 {
453 aParticleChange.ProposeEnergy(0.) ;
454 aParticleChange.ProposeTrackStatus(fStopAndKill);
455 }
456
457
458 // Kinematics of the scattered electron
459
460
461 G4double diffEnergy = photonEnergy0*(1-epsilon);
462 ionEnergy = (*(ionizationEnergy->find(Z)->second))[iosc];
463 //G4double eKineticEnergy = diffEnergy - ionEnergy;
464 G4double Q2 = photonEnergy0*photonEnergy0+photonEnergy1*(photonEnergy1-2.0*photonEnergy0*cosTheta);
465 G4double cosThetaE; //scattering angle for the electron
466 if (Q2 > 1.0e-12)
467 {
468 cosThetaE = (photonEnergy0-photonEnergy1*cosTheta)/std::sqrt(Q2);
469 }
470 else
471 {
472 cosThetaE = 1.0;
473 }
474 G4double sinThetaE = std::sqrt(1-cosThetaE*cosThetaE);
475
476
477
478 const G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
479 const G4AtomicShell* shell = transitionManager->Shell(Z,iosc);
480 G4double bindingEnergy = shell->BindingEnergy();
481 G4int shellId = shell->ShellId();
482 //G4cout << bindingEnergy/keV << " " << ionEnergy/keV << " keV" << G4endl;
483 ionEnergy = std::max(bindingEnergy,ionEnergy); //protection against energy non-conservation
484 G4double eKineticEnergy = diffEnergy - ionEnergy;
485
486 size_t nTotPhotons=0;
487 G4int nPhotons=0;
488
489 const G4ProductionCutsTable* theCoupleTable=
490 G4ProductionCutsTable::GetProductionCutsTable();
491 size_t indx = couple->GetIndex();
492 G4double cutg = (*(theCoupleTable->GetEnergyCutsVector(0)))[indx];
493 cutg = std::min(cutForLowEnergySecondaryPhotons,cutg);
494
495 G4double cute = (*(theCoupleTable->GetEnergyCutsVector(1)))[indx];
496 cute = std::min(cutForLowEnergySecondaryPhotons,cute);
497
498 std::vector<G4DynamicParticle*>* photonVector=0;
499 G4DynamicParticle* aPhoton;
500 G4AtomicDeexcitation deexcitationManager;
501
502 if (Z>5 && (ionEnergy > cutg || ionEnergy > cute))
503 {
504 photonVector = deexcitationManager.GenerateParticles(Z,shellId);
505 nTotPhotons = photonVector->size();
506 for (size_t k=0;k<nTotPhotons;k++){
507 aPhoton = (*photonVector)[k];
508 if (aPhoton)
509 {
510 G4double itsCut = cutg;
511 if (aPhoton->GetDefinition() == G4Electron::Electron()) itsCut = cute;
512 G4double itsEnergy = aPhoton->GetKineticEnergy();
513 if (itsEnergy > itsCut && itsEnergy <= ionEnergy)
514 {
515 nPhotons++;
516 ionEnergy -= itsEnergy;
517 }
518 else
519 {
520 delete aPhoton;
521 (*photonVector)[k]=0;
522 }
523 }
524 }
525 }
526 G4double energyDeposit =ionEnergy; //il deposito locale e' quello che rimane
527 G4int nbOfSecondaries=nPhotons;
528
529 // Generate the electron only if with large enough range w.r.t. cuts and safety
530 G4double safety = aStep.GetPostStepPoint()->GetSafety();
531 G4DynamicParticle* electron = 0;
532 if (rangeTest->Escape(G4Electron::Electron(),couple,eKineticEnergy,safety))
533 {
534 G4double xEl = sinThetaE * std::cos(phi+pi);
535 G4double yEl = sinThetaE * std::sin(phi+pi);
536 G4double zEl = cosThetaE;
537 G4ThreeVector eDirection(xEl,yEl,zEl); //electron direction
538 eDirection.rotateUz(photonDirection0);
539 electron = new G4DynamicParticle (G4Electron::Electron(),
540 eDirection,eKineticEnergy) ;
541 nbOfSecondaries++;
542 }
543 else
544 {
545
546 energyDeposit += eKineticEnergy;
547 }
548
549 aParticleChange.SetNumberOfSecondaries(nbOfSecondaries);
550 if (electron) aParticleChange.AddSecondary(electron);
551 for (size_t ll=0;ll<nTotPhotons;ll++)
552 {
553 aPhoton = (*photonVector)[ll];
554 if (aPhoton) aParticleChange.AddSecondary(aPhoton);
555 }
556 delete photonVector;
557 if (energyDeposit < 0)
558 {
559 G4cout << "WARNING-"
560 << "G4PenelopeCompton::PostStepDoIt - Negative energy deposit"
561 << G4endl;
562 energyDeposit=0;
563 }
564 aParticleChange.ProposeLocalEnergyDeposit(energyDeposit);
565
566
567 return G4VDiscreteProcess::PostStepDoIt( aTrack, aStep);
568}
569
570G4bool G4PenelopeCompton::IsApplicable(const G4ParticleDefinition& particle)
571{
572 return ( &particle == G4Gamma::Gamma() );
573}
574
575G4double G4PenelopeCompton::GetMeanFreePath(const G4Track& track,
576 G4double, // previousStepSize
577 G4ForceCondition*)
578{
579 const G4DynamicParticle* photon = track.GetDynamicParticle();
580 G4double energy = photon->GetKineticEnergy();
581 G4Material* material = track.GetMaterial();
582 size_t materialIndex = material->GetIndex();
583
584 G4double meanFreePath;
585 if (energy > highEnergyLimit) meanFreePath = meanFreePathTable->FindValue(highEnergyLimit,materialIndex);
586 else if (energy < lowEnergyLimit) meanFreePath = DBL_MAX;
587 else meanFreePath = meanFreePathTable->FindValue(energy,materialIndex);
588 return meanFreePath;
589}
590
591
592void G4PenelopeCompton::ReadData()
593{
594 char* path = getenv("G4LEDATA");
595 if (!path)
596 {
597 G4String excep = "G4PenelopeCompton - G4LEDATA environment variable not set!";
598 G4Exception(excep);
599 }
600 G4String pathString(path);
601 G4String pathFile = pathString + "/penelope/compton-pen.dat";
602 std::ifstream file(pathFile);
603 std::filebuf* lsdp = file.rdbuf();
604
605 if (!(lsdp->is_open()))
606 {
607 G4String excep = "G4PenelopeCompton - data file " + pathFile + " not found!";
608 G4Exception(excep);
609 }
610
611 G4int k1,test,test1;
612 G4double a1,a2;
613 G4int Z=1,nLevels=0;
614 G4DataVector* f;
615 G4DataVector* u;
616 G4DataVector* j;
617
618 do{
619 f = new G4DataVector;
620 u = new G4DataVector;
621 j = new G4DataVector;
622 file >> Z >> nLevels;
623 for (G4int h=0;h<nLevels;h++){
624 file >> k1 >> a1 >> a2;
625 f->push_back((G4double) k1);
626 u->push_back(a1);
627 j->push_back(a2);
628 }
629 ionizationEnergy->insert(std::make_pair(Z,u));
630 hartreeFunction->insert(std::make_pair(Z,j));
631 occupationNumber->insert(std::make_pair(Z,f));
632 file >> test >> test1; //-1 -1 close the data for each Z
633 if (test > 0) {
634 G4String excep = "G4PenelopeCompton - data file corrupted!";
635 G4Exception(excep);
636 }
637 }while (test != -2); //the very last Z is closed with -2 instead of -1
638}
639
640G4double G4PenelopeCompton::CrossSection(G4double energy,G4int Z)
641{
642 G4double cs=0.0;
643 energyForIntegration=energy;
644 ZForIntegration = Z;
645 if (energy< 5*MeV)
646 {
647 G4PenelopeIntegrator<G4PenelopeCompton,G4double (G4PenelopeCompton::*)(G4double)> theIntegrator;
648 cs = theIntegrator.Calculate(this,&G4PenelopeCompton::DifferentialCrossSection,-1.0,1.0,1e-05);
649 }
650 else
651 {
652 G4double ki=energy/electron_mass_c2;
653 G4double ki3=ki*ki;
654 G4double ki2=1.0+2*ki;
655 G4double ki1=ki3-ki2-1.0;
656 G4double t0=1.0/(ki2);
657 G4double csl = 0.5*ki3*t0*t0+ki2*t0+ki1*std::log(t0)-(1.0/t0);
658 G4int nosc = occupationNumber->find(Z)->second->size();
659 for (G4int i=0;i<nosc;i++)
660 {
661 G4double ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
662 G4double tau=(energy-ionEnergy)/energy;
663 if (tau > t0)
664 {
665 G4double csu = 0.5*ki3*tau*tau+ki2*tau+ki1*std::log(tau)-(1.0/tau);
666 G4int f = (G4int) (*(occupationNumber->find(Z)->second))[i];
667 cs = cs + f*(csu-csl);
668 }
669 }
670 cs=pi*classic_electr_radius*classic_electr_radius*cs/(ki*ki3);
671 }
672 return cs;
673}
674
675
676G4double G4PenelopeCompton::DifferentialCrossSection(G4double cosTheta)
677{
678 const G4double k2 = std::sqrt(2.0);
679 const G4double k1 = std::sqrt(0.5);
680 const G4double k12 = 0.5;
681 G4double cdt1 = 1.0-cosTheta;
682 G4double energy = energyForIntegration;
683 G4int Z = ZForIntegration;
684 G4double ionEnergy=0.0,Pzimax=0.0,XKN=0.0;
685 G4double diffCS=0.0;
686 G4double x=0.0,siap=0.0;
687 G4double harFunc=0.0;
688 G4int occupNb;
689 //energy of Compton line;
690 G4double EOEC = 1.0+(energy/electron_mass_c2)*cdt1;
691 G4double ECOE = 1.0/EOEC;
692 //Incoherent scattering function (analytical profile)
693 G4double sia = 0.0;
694 G4int nosc = occupationNumber->find(Z)->second->size();
695 for (G4int i=0;i<nosc;i++){
696 ionEnergy = (*(ionizationEnergy->find(Z)->second))[i];
697 //Sum only of those shells for which E>Eion
698 if (energy > ionEnergy)
699 {
700 G4double aux = energy * (energy-ionEnergy)*cdt1;
701 Pzimax = (aux - electron_mass_c2*ionEnergy)/(electron_mass_c2*std::sqrt(2*aux+ionEnergy*ionEnergy));
702 harFunc = (*(hartreeFunction->find(Z)->second))[i]/fine_structure_const;
703 occupNb = (G4int) (*(occupationNumber->find(Z)->second))[i];
704 x = harFunc*Pzimax;
705 if (x > 0)
706 {
707 siap = 1.0-0.5*std::exp(k12-(k1+k2*x)*(k1+k2*x));
708 }
709 else
710 {
711 siap = 0.5*std::exp(k12-(k1-k2*x)*(k1-k2*x));
712 }
713 sia = sia + occupNb*siap; //sum of all contributions;
714 }
715 }
716 XKN = EOEC+ECOE-1+cosTheta*cosTheta;
717 diffCS = pi*classic_electr_radius*classic_electr_radius*ECOE*ECOE*XKN*sia;
718 return diffCS;
719}
720
721G4int G4PenelopeCompton::SelectRandomAtomForCompton(const G4Material* material,G4double energy) const
722{
723 G4int nElements = material->GetNumberOfElements();
724 //Special case: the material consists of one element
725 if (nElements == 1)
726 {
727 G4int Z = (G4int) material->GetZ();
728 return Z;
729 }
730
731 //Composite material
732 const G4ElementVector* elementVector = material->GetElementVector();
733 size_t materialIndex = material->GetIndex();
734
735 G4VEMDataSet* materialSet = (*matCrossSections)[materialIndex];
736 G4double materialCrossSection0 = 0.0;
737 G4DataVector cross;
738 cross.clear();
739 G4int i;
740 for (i=0;i<nElements;i++)
741 {
742 G4double cr = (materialSet->GetComponent(i))->FindValue(energy);
743 materialCrossSection0 += cr;
744 cross.push_back(materialCrossSection0); //cumulative cross section
745 }
746
747 G4double random = G4UniformRand()*materialCrossSection0;
748 for (i=0;i<nElements;i++)
749 {
750 if (random <= cross[i]) return (G4int) (*elementVector)[i]->GetZ();
751 }
752 //It should never get here
753 return 0;
754}
755
Note: See TracBrowser for help on using the repository browser.