source: trunk/source/processes/electromagnetic/xrays/src/G4ForwardXrayTR.cc@ 1201

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 31.2 KB
Line 
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//
27// $Id: G4ForwardXrayTR.cc,v 1.14 2007/05/11 14:23:04 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-03-cand-01 $
29//
30// G4ForwardXrayTR class -- implementation file
31
32// GEANT 4 class implementation file --- Copyright CERN 1995
33// CERN Geneva Switzerland
34
35// For information related to this code, please, contact
36// CERN, CN Division, ASD Group
37// History:
38// 1st version 11.09.97 V. Grichine (Vladimir.Grichine@cern.ch )
39// 2nd version 17.12.97 V. Grichine
40// 17-09-01, migration of Materials to pure STL (mma)
41// 10-03-03, migration to "cut per region" (V.Ivanchenko)
42// 03.06.03, V.Ivanchenko fix compilation warnings
43
44#include "G4ForwardXrayTR.hh"
45
46#include "globals.hh"
47#include "G4Poisson.hh"
48#include "G4Material.hh"
49#include "G4PhysicsTable.hh"
50#include "G4PhysicsVector.hh"
51#include "G4PhysicsLinearVector.hh"
52#include "G4PhysicsLogVector.hh"
53#include "G4ProductionCutsTable.hh"
54#include "G4GeometryTolerance.hh"
55
56// Table initialization
57
58// G4PhysicsTable* G4ForwardXrayTR::fAngleDistrTable = NULL ;
59// G4PhysicsTable* G4ForwardXrayTR::fEnergyDistrTable = NULL ;
60
61
62// Initialization of local constants
63
64G4int G4ForwardXrayTR::fSympsonNumber = 100 ;
65
66G4double G4ForwardXrayTR::fTheMinEnergyTR = 1.0*keV ;
67G4double G4ForwardXrayTR::fTheMaxEnergyTR = 100.0*keV ;
68G4double G4ForwardXrayTR::fTheMaxAngle = 1.0e-3 ;
69G4double G4ForwardXrayTR::fTheMinAngle = 5.0e-6 ;
70G4int G4ForwardXrayTR::fBinTR = 50 ;
71
72G4double G4ForwardXrayTR::fMinProtonTkin = 100.0*GeV ;
73G4double G4ForwardXrayTR::fMaxProtonTkin = 100.0*TeV ;
74G4int G4ForwardXrayTR::fTotBin = 50 ;
75// Proton energy vector initialization
76
77G4PhysicsLogVector* G4ForwardXrayTR::
78fProtonEnergyVector = new G4PhysicsLogVector(fMinProtonTkin,
79 fMaxProtonTkin,
80 fTotBin ) ;
81
82G4double G4ForwardXrayTR::fPlasmaCof = 4.0*pi*fine_structure_const*
83 hbarc*hbarc*hbarc/electron_mass_c2 ;
84
85G4double G4ForwardXrayTR::fCofTR = fine_structure_const/pi ;
86
87using namespace std;
88
89/* ************************************************************************
90
91
92///////////////////////////////////////////////////////////////////////
93//
94// Constructor for preparation tables with angle and energy TR distributions
95// in all materials involved in test program. Lorentz factors correspond to
96// kinetic energies of protons between 100*GeV and 100*TeV, ~ 10^2-10^5
97//
98// Recommended only for use in applications with
99// few light materials involved !!!!!!!!!!!!!!
100
101G4ForwardXrayTR::G4ForwardXrayTR()
102 : G4TransitionRadiation("XrayTR")
103{
104 G4int iMat, jMat, iTkin, iTR, iPlace ;
105 static
106 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
107 G4int numOfMat = G4Material::GetNumberOfMaterials();
108 fGammaCutInKineticEnergy = new G4double[numOfMat] ;
109 fGammaCutInKineticEnergy = fPtrGamma->GetEnergyCuts() ;
110 fMatIndex1 = -1 ;
111 fMatIndex2 = -1 ;
112 fAngleDistrTable = new G4PhysicsTable(numOfMat*(numOfMat - 1)*fTotBin) ;
113 fEnergyDistrTable = new G4PhysicsTable(numOfMat*(numOfMat - 1)*fTotBin) ;
114
115 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(fMinProtonTkin,
116 fMaxProtonTkin,
117 fTotBin ) ;
118
119 for(iMat=0;iMat<numOfMat;iMat++) // loop over pairs of different materials
120 {
121 for(jMat=0;jMat<numOfMat;jMat++) // transition iMat -> jMat !!!
122 {
123 if(iMat == jMat) continue ; // no TR !!
124 else
125 {
126 const G4Material* mat1 = (*theMaterialTable)[iMat] ;
127 const G4Material* mat2 = (*theMaterialTable)[jMat] ;
128
129 fSigma1 = fPlasmaCof*(mat1->GetElectronDensity()) ;
130 fSigma2 = fPlasmaCof*(mat2->GetElectronDensity()) ;
131
132// fGammaTkinCut = fGammaCutInKineticEnergy[jMat] ; // TR photon in jMat !
133 fGammaTkinCut = 0.0 ;
134
135 if(fGammaTkinCut > fTheMinEnergyTR) // setting of min/max TR energies
136 {
137 fMinEnergyTR = fGammaTkinCut ;
138 }
139 else
140 {
141 fMinEnergyTR = fTheMinEnergyTR ;
142 }
143 if(fGammaTkinCut > fTheMaxEnergyTR)
144 {
145 fMaxEnergyTR = 2.0*fGammaTkinCut ; // usually very low TR rate
146 }
147 else
148 {
149 fMaxEnergyTR = fTheMaxEnergyTR ;
150 }
151 for(iTkin=0;iTkin<fTotBin;iTkin++) // Lorentz factor loop
152 {
153 G4PhysicsLogVector*
154 energyVector = new G4PhysicsLogVector(fMinEnergyTR,
155 fMaxEnergyTR,
156 fBinTR ) ;
157 G4PhysicsLinearVector*
158 angleVector = new G4PhysicsLinearVector( 0.0,
159 fMaxThetaTR,
160 fBinTR ) ;
161 G4double energySum = 0.0 ;
162 G4double angleSum = 0.0 ;
163 fGamma = 1.0 + (aVector->GetLowEdgeEnergy(iTkin)/proton_mass_c2) ;
164 fMaxThetaTR = 10000.0/(fGamma*fGamma) ;
165 if(fMaxThetaTR > fTheMaxAngle)
166 {
167 fMaxThetaTR = fTheMaxAngle ;
168 }
169 else
170 {
171 if(fMaxThetaTR < fTheMinAngle)
172 {
173 fMaxThetaTR = fTheMinAngle ;
174 }
175 }
176 energyVector->PutValue(fBinTR-1,energySum) ;
177 angleVector->PutValue(fBinTR-1,angleSum) ;
178
179 for(iTR=fBinTR-2;iTR>=0;iTR--)
180 {
181 energySum += fCofTR*EnergySum(energyVector->GetLowEdgeEnergy(iTR),
182 energyVector->GetLowEdgeEnergy(iTR+1)) ;
183
184 angleSum += fCofTR*AngleSum(angleVector->GetLowEdgeEnergy(iTR),
185 angleVector->GetLowEdgeEnergy(iTR+1)) ;
186 energyVector->PutValue(iTR,energySum) ;
187 angleVector->PutValue(iTR,angleSum) ;
188 }
189 if(jMat < iMat)
190 {
191 iPlace = (iMat*(numOfMat-1)+jMat)*fTotBin+iTkin ;
192 }
193 else // jMat > iMat right part of matrices (jMat-1) !
194 {
195 iPlace = (iMat*(numOfMat-1)+jMat-1)*fTotBin+iTkin ;
196 }
197 fEnergyDistrTable->insertAt(iPlace,energyVector) ;
198 fAngleDistrTable->insertAt(iPlace,angleVector) ;
199 } // iTkin
200 } // jMat != iMat
201 } // jMat
202 } // iMat
203}
204
205
206**************************************************************** */
207
208
209//////////////////////////////////////////////////////////////////////
210//
211// Constructor for creation of physics tables (angle and energy TR
212// distributions) for a couple of selected materials.
213//
214// Recommended for use in applications with many materials involved,
215// when only few (usually couple) materials are interested for generation
216// of TR on the interface between them
217
218
219G4ForwardXrayTR::
220G4ForwardXrayTR( const G4String& matName1, // G4Material* pMat1,
221 const G4String& matName2, // G4Material* pMat2,
222 const G4String& processName )
223 : G4TransitionRadiation(processName)
224{
225 // fMatIndex1 = pMat1->GetIndex() ;
226 // fMatIndex2 = pMat2->GetIndex() ;
227 G4int iMat;
228 const G4ProductionCutsTable* theCoupleTable=
229 G4ProductionCutsTable::GetProductionCutsTable();
230 G4int numOfCouples = theCoupleTable->GetTableSize();
231
232 for(iMat=0;iMat<numOfCouples;iMat++) // check first material name
233 {
234 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(iMat);
235 if( matName1 == couple->GetMaterial()->GetName() )
236 {
237 fMatIndex1 = couple->GetIndex() ;
238 break ;
239 }
240 }
241 if(iMat == numOfCouples)
242 {
243 G4Exception("Invalid first material name in G4ForwardXrayTR constructor") ;
244 }
245
246 for(iMat=0;iMat<numOfCouples;iMat++) // check second material name
247 {
248 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(iMat);
249 if( matName2 == couple->GetMaterial()->GetName() )
250 {
251 fMatIndex2 = couple->GetIndex() ;
252 break ;
253 }
254 }
255 if(iMat == numOfCouples)
256 {
257 G4Exception("Invalid second material name in G4ForwardXrayTR constructor") ;
258 }
259 // G4cout<<"G4ForwardXray constructor is called"<<G4endl ;
260 BuildXrayTRtables() ;
261}
262
263/////////////////////////////////////////////////////////////////////////
264//
265// Constructor used by X-ray transition radiation parametrisation models
266
267G4ForwardXrayTR::
268G4ForwardXrayTR( const G4String& processName )
269 : G4TransitionRadiation(processName)
270{
271 ;
272}
273
274
275//////////////////////////////////////////////////////////////////////
276//
277// Destructor
278//
279
280G4ForwardXrayTR::~G4ForwardXrayTR()
281{
282 ;
283}
284
285//////////////////////////////////////////////////////////////////////////////
286//
287// Build physics tables for energy and angular distributions of X-ray TR photon
288
289void G4ForwardXrayTR::BuildXrayTRtables()
290{
291 G4int iMat, jMat, iTkin, iTR, iPlace ;
292 const G4ProductionCutsTable* theCoupleTable=
293 G4ProductionCutsTable::GetProductionCutsTable();
294 G4int numOfCouples = theCoupleTable->GetTableSize();
295
296 fGammaCutInKineticEnergy = theCoupleTable->GetEnergyCutsVector(idxG4GammaCut);
297
298 fAngleDistrTable = new G4PhysicsTable(2*fTotBin) ;
299 fEnergyDistrTable = new G4PhysicsTable(2*fTotBin) ;
300
301
302 for(iMat=0;iMat<numOfCouples;iMat++) // loop over pairs of different materials
303 {
304 if( iMat != fMatIndex1 && iMat != fMatIndex2 ) continue ;
305
306 for(jMat=0;jMat<numOfCouples;jMat++) // transition iMat -> jMat !!!
307 {
308 if( iMat == jMat || ( jMat != fMatIndex1 && jMat != fMatIndex2 ) )
309 {
310 continue ;
311 }
312 else
313 {
314 const G4MaterialCutsCouple* iCouple = theCoupleTable->GetMaterialCutsCouple(iMat);
315 const G4MaterialCutsCouple* jCouple = theCoupleTable->GetMaterialCutsCouple(jMat);
316 const G4Material* mat1 = iCouple->GetMaterial() ;
317 const G4Material* mat2 = jCouple->GetMaterial() ;
318
319 fSigma1 = fPlasmaCof*(mat1->GetElectronDensity()) ;
320 fSigma2 = fPlasmaCof*(mat2->GetElectronDensity()) ;
321
322 // fGammaTkinCut = fGammaCutInKineticEnergy[jMat] ; // TR photon in jMat !
323
324 fGammaTkinCut = 0.0 ;
325
326 if(fGammaTkinCut > fTheMinEnergyTR) // setting of min/max TR energies
327 {
328 fMinEnergyTR = fGammaTkinCut ;
329 }
330 else
331 {
332 fMinEnergyTR = fTheMinEnergyTR ;
333 }
334 if(fGammaTkinCut > fTheMaxEnergyTR)
335 {
336 fMaxEnergyTR = 2.0*fGammaTkinCut ; // usually very low TR rate
337 }
338 else
339 {
340 fMaxEnergyTR = fTheMaxEnergyTR ;
341 }
342 for(iTkin=0;iTkin<fTotBin;iTkin++) // Lorentz factor loop
343 {
344 G4PhysicsLogVector*
345 energyVector = new G4PhysicsLogVector( fMinEnergyTR,
346 fMaxEnergyTR,
347 fBinTR ) ;
348
349 fGamma = 1.0 + (fProtonEnergyVector->
350 GetLowEdgeEnergy(iTkin)/proton_mass_c2) ;
351
352 fMaxThetaTR = 10000.0/(fGamma*fGamma) ;
353
354 if(fMaxThetaTR > fTheMaxAngle)
355 {
356 fMaxThetaTR = fTheMaxAngle ;
357 }
358 else
359 {
360 if(fMaxThetaTR < fTheMinAngle)
361 {
362 fMaxThetaTR = fTheMinAngle ;
363 }
364 }
365 // G4cout<<G4endl<<"fGamma = "<<fGamma<<" fMaxThetaTR = "<<fMaxThetaTR<<G4endl ;
366 G4PhysicsLinearVector*
367 angleVector = new G4PhysicsLinearVector( 0.0,
368 fMaxThetaTR,
369 fBinTR ) ;
370 G4double energySum = 0.0 ;
371 G4double angleSum = 0.0 ;
372
373 energyVector->PutValue(fBinTR-1,energySum) ;
374 angleVector->PutValue(fBinTR-1,angleSum) ;
375
376 for(iTR=fBinTR-2;iTR>=0;iTR--)
377 {
378 energySum += fCofTR*EnergySum(energyVector->GetLowEdgeEnergy(iTR),
379 energyVector->GetLowEdgeEnergy(iTR+1)) ;
380
381 angleSum += fCofTR*AngleSum(angleVector->GetLowEdgeEnergy(iTR),
382 angleVector->GetLowEdgeEnergy(iTR+1)) ;
383
384 energyVector->PutValue(iTR,energySum) ;
385 angleVector ->PutValue(iTR,angleSum) ;
386 }
387 // G4cout<<"sumE = "<<energySum<<" ; sumA = "<<angleSum<<G4endl ;
388
389 if(jMat < iMat)
390 {
391 iPlace = fTotBin+iTkin ; // (iMat*(numOfMat-1)+jMat)*
392 }
393 else // jMat > iMat right part of matrices (jMat-1) !
394 {
395 iPlace = iTkin ; // (iMat*(numOfMat-1)+jMat-1)*fTotBin+
396 }
397 fEnergyDistrTable->insertAt(iPlace,energyVector) ;
398 fAngleDistrTable->insertAt(iPlace,angleVector) ;
399 } // iTkin
400 } // jMat != iMat
401 } // jMat
402 } // iMat
403 // G4cout<<"G4ForwardXrayTR::BuildXrayTRtables have been called"<<G4endl ;
404}
405
406///////////////////////////////////////////////////////////////////////
407//
408// This function returns the spectral and angle density of TR quanta
409// in X-ray energy region generated forward when a relativistic
410// charged particle crosses interface between two materials.
411// The high energy small theta approximation is applied.
412// (matter1 -> matter2)
413// varAngle =2* (1 - cos(Theta)) or approximately = Theta*Theta
414//
415
416G4double
417G4ForwardXrayTR::SpectralAngleTRdensity( G4double energy,
418 G4double varAngle ) const
419{
420 G4double formationLength1, formationLength2 ;
421 formationLength1 = 1.0/
422 (1.0/(fGamma*fGamma)
423 + fSigma1/(energy*energy)
424 + varAngle) ;
425 formationLength2 = 1.0/
426 (1.0/(fGamma*fGamma)
427 + fSigma2/(energy*energy)
428 + varAngle) ;
429 return (varAngle/energy)*(formationLength1 - formationLength2)
430 *(formationLength1 - formationLength2) ;
431
432}
433
434
435//////////////////////////////////////////////////////////////////
436//
437// Analytical formula for angular density of X-ray TR photons
438//
439
440G4double G4ForwardXrayTR::AngleDensity( G4double energy,
441 G4double varAngle ) const
442{
443 G4double x, x2, a, b, c, d, f, a2, b2, a4, b4 ;
444 G4double cof1, cof2, cof3 ;
445 x = 1.0/energy ;
446 x2 = x*x ;
447 c = 1.0/fSigma1 ;
448 d = 1.0/fSigma2 ;
449 f = (varAngle + 1.0/(fGamma*fGamma)) ;
450 a2 = c*f ;
451 b2 = d*f ;
452 a4 = a2*a2 ;
453 b4 = b2*b2 ;
454 a = sqrt(a2) ;
455 b = sqrt(b2) ;
456 cof1 = c*c*(0.5/(a2*(x2 +a2)) +0.5*log(x2/(x2 +a2))/a4) ;
457 cof3 = d*d*(0.5/(b2*(x2 +b2)) +0.5*log(x2/(x2 +b2))/b4) ;
458 cof2 = -c*d*(log(x2/(x2 +b2))/b2 - log(x2/(x2 +a2))/a2)/(a2 - b2) ;
459 return -varAngle*(cof1 + cof2 + cof3) ;
460}
461
462/////////////////////////////////////////////////////////////////////
463//
464// Definite integral of X-ray TR spectral-angle density from energy1
465// to energy2
466//
467
468G4double G4ForwardXrayTR::EnergyInterval( G4double energy1,
469 G4double energy2,
470 G4double varAngle ) const
471{
472 return AngleDensity(energy2,varAngle)
473 - AngleDensity(energy1,varAngle) ;
474}
475
476//////////////////////////////////////////////////////////////////////
477//
478// Integral angle distribution of X-ray TR photons based on analytical
479// formula for angle density
480//
481
482G4double G4ForwardXrayTR::AngleSum( G4double varAngle1,
483 G4double varAngle2 ) const
484{
485 G4int i ;
486 G4double h , sumEven = 0.0 , sumOdd = 0.0 ;
487 h = 0.5*(varAngle2 - varAngle1)/fSympsonNumber ;
488 for(i=1;i<fSympsonNumber;i++)
489 {
490 sumEven += EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle1 + 2*i*h ) ;
491 sumOdd += EnergyInterval(fMinEnergyTR,fMaxEnergyTR,
492 varAngle1 + (2*i - 1)*h ) ;
493 }
494 sumOdd += EnergyInterval(fMinEnergyTR,fMaxEnergyTR,
495 varAngle1 + (2*fSympsonNumber - 1)*h ) ;
496
497 return h*(EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle1)
498 + EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle2)
499 + 4.0*sumOdd + 2.0*sumEven )/3.0 ;
500}
501
502/////////////////////////////////////////////////////////////////////
503//
504// Analytical Expression for spectral density of Xray TR photons
505// x = 2*(1 - cos(Theta)) ~ Theta^2
506//
507
508G4double G4ForwardXrayTR::SpectralDensity( G4double energy,
509 G4double x ) const
510{
511 G4double a, b ;
512 a = 1.0/(fGamma*fGamma)
513 + fSigma1/(energy*energy) ;
514 b = 1.0/(fGamma*fGamma)
515 + fSigma2/(energy*energy) ;
516 return ( (a + b)*log((x + b)/(x + a))/(a - b)
517 + a/(x + a) + b/(x + b) )/energy ;
518
519}
520
521////////////////////////////////////////////////////////////////////
522//
523// The spectral density in some angle interval from varAngle1 to
524// varAngle2
525//
526
527G4double G4ForwardXrayTR::AngleInterval( G4double energy,
528 G4double varAngle1,
529 G4double varAngle2 ) const
530{
531 return SpectralDensity(energy,varAngle2)
532 - SpectralDensity(energy,varAngle1) ;
533}
534
535////////////////////////////////////////////////////////////////////
536//
537// Integral spectral distribution of X-ray TR photons based on
538// analytical formula for spectral density
539//
540
541G4double G4ForwardXrayTR::EnergySum( G4double energy1,
542 G4double energy2 ) const
543{
544 G4int i ;
545 G4double h , sumEven = 0.0 , sumOdd = 0.0 ;
546 h = 0.5*(energy2 - energy1)/fSympsonNumber ;
547 for(i=1;i<fSympsonNumber;i++)
548 {
549 sumEven += AngleInterval(energy1 + 2*i*h,0.0,fMaxThetaTR);
550 sumOdd += AngleInterval(energy1 + (2*i - 1)*h,0.0,fMaxThetaTR) ;
551 }
552 sumOdd += AngleInterval(energy1 + (2*fSympsonNumber - 1)*h,
553 0.0,fMaxThetaTR) ;
554
555 return h*( AngleInterval(energy1,0.0,fMaxThetaTR)
556 + AngleInterval(energy2,0.0,fMaxThetaTR)
557 + 4.0*sumOdd + 2.0*sumEven )/3.0 ;
558}
559
560/////////////////////////////////////////////////////////////////////////
561//
562// PostStepDoIt function for creation of forward X-ray photons in TR process
563// on boubndary between two materials with really different plasma energies
564//
565
566G4VParticleChange* G4ForwardXrayTR::PostStepDoIt(const G4Track& aTrack,
567 const G4Step& aStep)
568{
569 aParticleChange.Initialize(aTrack);
570 // G4cout<<"call G4ForwardXrayTR::PostStepDoIt"<<G4endl ;
571 G4int iMat, jMat, iTkin, iPlace, numOfTR, iTR, iTransfer ;
572
573 G4double energyPos, anglePos, energyTR, theta, phi, dirX, dirY, dirZ ;
574 G4double W, W1, W2, E1, E2 ;
575
576 G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
577 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
578 G4double tol=0.5*G4GeometryTolerance::GetInstance()->GetSurfaceTolerance();
579
580 if (pPostStepPoint->GetStepStatus() != fGeomBoundary)
581 {
582 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
583 }
584 if (aTrack.GetStepLength() <= tol)
585 {
586 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
587 }
588 // Come on boundary, so begin to try TR
589
590 const G4MaterialCutsCouple* iCouple = pPreStepPoint ->GetPhysicalVolume()->
591 GetLogicalVolume()->GetMaterialCutsCouple();
592 const G4MaterialCutsCouple* jCouple = pPostStepPoint ->GetPhysicalVolume()->
593 GetLogicalVolume()->GetMaterialCutsCouple();
594 const G4Material* iMaterial = iCouple->GetMaterial();
595 const G4Material* jMaterial = jCouple->GetMaterial();
596 iMat = iCouple->GetIndex();
597 jMat = jCouple->GetIndex();
598
599 // The case of equal or approximate (in terms of plasma energy) materials
600 // No TR photons ?!
601
602 if ( iMat == jMat
603 || ( (fMatIndex1 >= 0 && fMatIndex1 >= 0)
604 && ( iMat != fMatIndex1 && iMat != fMatIndex2 )
605 && ( jMat != fMatIndex1 && jMat != fMatIndex2 ) )
606
607 || iMaterial->GetState() == jMaterial->GetState()
608
609 ||(iMaterial->GetState() == kStateSolid && jMaterial->GetState() == kStateLiquid )
610
611 ||(iMaterial->GetState() == kStateLiquid && jMaterial->GetState() == kStateSolid ) )
612 {
613 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep) ;
614 }
615
616 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
617 G4double charge = aParticle->GetDefinition()->GetPDGCharge();
618
619 if(charge == 0.0) // Uncharged particle doesn't Generate TR photons
620 {
621 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
622 }
623 // Now we are ready to Generate TR photons
624
625 G4double chargeSq = charge*charge ;
626 G4double kinEnergy = aParticle->GetKineticEnergy() ;
627 G4double massRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass() ;
628 G4double TkinScaled = kinEnergy*massRatio ;
629 for(iTkin=0;iTkin<fTotBin;iTkin++)
630 {
631 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) // <= ?
632 {
633 break ;
634 }
635 }
636 if(jMat < iMat)
637 {
638 iPlace = fTotBin + iTkin - 1 ; // (iMat*(numOfMat - 1) + jMat)*
639 }
640 else
641 {
642 iPlace = iTkin - 1 ; // (iMat*(numOfMat - 1) + jMat - 1)*fTotBin +
643 }
644 // G4PhysicsVector* energyVector1 = (*fEnergyDistrTable)(iPlace) ;
645 // G4PhysicsVector* energyVector2 = (*fEnergyDistrTable)(iPlace + 1) ;
646
647 // G4PhysicsVector* angleVector1 = (*fAngleDistrTable)(iPlace) ;
648 // G4PhysicsVector* angleVector2 = (*fAngleDistrTable)(iPlace + 1) ;
649
650 G4ParticleMomentum particleDir = aParticle->GetMomentumDirection() ;
651
652 if(iTkin == fTotBin) // TR plato, try from left
653 {
654 // G4cout<<iTkin<<" mean TR number = "<<( (*(*fEnergyDistrTable)(iPlace))(0) +
655 // (*(*fAngleDistrTable)(iPlace))(0) )
656 // *chargeSq*0.5<<G4endl ;
657
658 numOfTR = G4Poisson( ( (*(*fEnergyDistrTable)(iPlace))(0) +
659 (*(*fAngleDistrTable)(iPlace))(0) )
660 *chargeSq*0.5 ) ;
661 if(numOfTR == 0)
662 {
663 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
664 }
665 else
666 {
667 // G4cout<<"Number of X-ray TR photons = "<<numOfTR<<G4endl ;
668
669 aParticleChange.SetNumberOfSecondaries(numOfTR);
670
671 for(iTR=0;iTR<numOfTR;iTR++)
672 {
673 energyPos = (*(*fEnergyDistrTable)(iPlace))(0)*G4UniformRand() ;
674 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
675 {
676 if(energyPos >= (*(*fEnergyDistrTable)(iPlace))(iTransfer)) break ;
677 }
678 energyTR = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
679
680 // G4cout<<"energyTR = "<<energyTR/keV<<"keV"<<G4endl ;
681
682 kinEnergy -= energyTR ;
683 aParticleChange.ProposeEnergy(kinEnergy);
684
685 anglePos = (*(*fAngleDistrTable)(iPlace))(0)*G4UniformRand() ;
686 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
687 {
688 if(anglePos > (*(*fAngleDistrTable)(iPlace))(iTransfer)) break ;
689 }
690 theta = sqrt((*fAngleDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1)) ;
691
692 // G4cout<<iTransfer<<" : theta = "<<theta<<G4endl ;
693
694 phi = twopi*G4UniformRand() ;
695 dirX = sin(theta)*cos(phi) ;
696 dirY = sin(theta)*sin(phi) ;
697 dirZ = cos(theta) ;
698 G4ThreeVector directionTR(dirX,dirY,dirZ) ;
699 directionTR.rotateUz(particleDir) ;
700 G4DynamicParticle* aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(),
701 directionTR,
702 energyTR ) ;
703 aParticleChange.AddSecondary(aPhotonTR) ;
704 }
705 }
706 }
707 else
708 {
709 if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
710 {
711 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
712 }
713 else // general case: Tkin between two vectors of the material
714 {
715 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ;
716 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ;
717 W = 1.0/(E2 - E1) ;
718 W1 = (E2 - TkinScaled)*W ;
719 W2 = (TkinScaled - E1)*W ;
720
721 // G4cout<<iTkin<<" mean TR number = "<<(((*(*fEnergyDistrTable)(iPlace))(0)+
722 // (*(*fAngleDistrTable)(iPlace))(0))*W1 +
723 // ((*(*fEnergyDistrTable)(iPlace + 1))(0)+
724 // (*(*fAngleDistrTable)(iPlace + 1))(0))*W2)
725 // *chargeSq*0.5<<G4endl ;
726
727 numOfTR = G4Poisson((((*(*fEnergyDistrTable)(iPlace))(0)+
728 (*(*fAngleDistrTable)(iPlace))(0))*W1 +
729 ((*(*fEnergyDistrTable)(iPlace + 1))(0)+
730 (*(*fAngleDistrTable)(iPlace + 1))(0))*W2)
731 *chargeSq*0.5 ) ;
732 if(numOfTR == 0)
733 {
734 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
735 }
736 else
737 {
738 // G4cout<<"Number of X-ray TR photons = "<<numOfTR<<G4endl ;
739
740 aParticleChange.SetNumberOfSecondaries(numOfTR);
741 for(iTR=0;iTR<numOfTR;iTR++)
742 {
743 energyPos = ((*(*fEnergyDistrTable)(iPlace))(0)*W1+
744 (*(*fEnergyDistrTable)(iPlace + 1))(0)*W2)*G4UniformRand() ;
745 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
746 {
747 if(energyPos >= ((*(*fEnergyDistrTable)(iPlace))(iTransfer)*W1+
748 (*(*fEnergyDistrTable)(iPlace + 1))(iTransfer)*W2)) break ;
749 }
750 energyTR = ((*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer))*W1+
751 ((*fEnergyDistrTable)(iPlace + 1)->GetLowEdgeEnergy(iTransfer))*W2 ;
752
753 // G4cout<<"energyTR = "<<energyTR/keV<<"keV"<<G4endl ;
754
755 kinEnergy -= energyTR ;
756 aParticleChange.ProposeEnergy(kinEnergy);
757
758 anglePos = ((*(*fAngleDistrTable)(iPlace))(0)*W1+
759 (*(*fAngleDistrTable)(iPlace + 1))(0)*W2)*G4UniformRand() ;
760 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
761 {
762 if(anglePos > ((*(*fAngleDistrTable)(iPlace))(iTransfer)*W1+
763 (*(*fAngleDistrTable)(iPlace + 1))(iTransfer)*W2)) break ;
764 }
765 theta = sqrt(((*fAngleDistrTable)(iPlace)->
766 GetLowEdgeEnergy(iTransfer-1))*W1+
767 ((*fAngleDistrTable)(iPlace + 1)->
768 GetLowEdgeEnergy(iTransfer-1))*W2) ;
769
770 // G4cout<<iTransfer<<" : theta = "<<theta<<G4endl ;
771
772 phi = twopi*G4UniformRand() ;
773 dirX = sin(theta)*cos(phi) ;
774 dirY = sin(theta)*sin(phi) ;
775 dirZ = cos(theta) ;
776 G4ThreeVector directionTR(dirX,dirY,dirZ) ;
777 directionTR.rotateUz(particleDir) ;
778 G4DynamicParticle* aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(),
779 directionTR,
780 energyTR ) ;
781 aParticleChange.AddSecondary(aPhotonTR) ;
782 }
783 }
784 }
785 }
786 return &aParticleChange ;
787}
788
789////////////////////////////////////////////////////////////////////////////
790//
791// Test function for checking of PostStepDoIt random preparation of TR photon
792// energy
793//
794
795G4double
796G4ForwardXrayTR::GetEnergyTR(G4int iMat, G4int jMat, G4int iTkin) const
797{
798 G4int iPlace, numOfTR, iTR, iTransfer ;
799 G4double energyTR = 0.0 ; // return this value for no TR photons
800 G4double energyPos ;
801 G4double W1, W2;
802
803 const G4ProductionCutsTable* theCoupleTable=
804 G4ProductionCutsTable::GetProductionCutsTable();
805 G4int numOfCouples = theCoupleTable->GetTableSize();
806
807 // The case of equal or approximate (in terms of plasma energy) materials
808 // No TR photons ?!
809
810 const G4MaterialCutsCouple* iCouple = theCoupleTable->GetMaterialCutsCouple(iMat);
811 const G4MaterialCutsCouple* jCouple = theCoupleTable->GetMaterialCutsCouple(jMat);
812 const G4Material* iMaterial = iCouple->GetMaterial();
813 const G4Material* jMaterial = jCouple->GetMaterial();
814
815 if ( iMat == jMat
816
817 || iMaterial->GetState() == jMaterial->GetState()
818
819 ||(iMaterial->GetState() == kStateSolid && jMaterial->GetState() == kStateLiquid )
820
821 ||(iMaterial->GetState() == kStateLiquid && jMaterial->GetState() == kStateSolid ) )
822
823 {
824 return energyTR ;
825 }
826
827 if(jMat < iMat)
828 {
829 iPlace = (iMat*(numOfCouples - 1) + jMat)*fTotBin + iTkin - 1 ;
830 }
831 else
832 {
833 iPlace = (iMat*(numOfCouples - 1) + jMat - 1)*fTotBin + iTkin - 1 ;
834 }
835 G4PhysicsVector* energyVector1 = (*fEnergyDistrTable)(iPlace) ;
836 G4PhysicsVector* energyVector2 = (*fEnergyDistrTable)(iPlace + 1) ;
837
838 if(iTkin == fTotBin) // TR plato, try from left
839 {
840 numOfTR = G4Poisson( (*energyVector1)(0) ) ;
841 if(numOfTR == 0)
842 {
843 return energyTR ;
844 }
845 else
846 {
847 for(iTR=0;iTR<numOfTR;iTR++)
848 {
849 energyPos = (*energyVector1)(0)*G4UniformRand() ;
850 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
851 {
852 if(energyPos >= (*energyVector1)(iTransfer)) break ;
853 }
854 energyTR += energyVector1->GetLowEdgeEnergy(iTransfer) ;
855 }
856 }
857 }
858 else
859 {
860 if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
861 {
862 return energyTR ;
863 }
864 else // general case: Tkin between two vectors of the material
865 { // use trivial mean half/half
866 W1 = 0.5 ;
867 W2 = 0.5 ;
868 numOfTR = G4Poisson( (*energyVector1)(0)*W1 +
869 (*energyVector2)(0)*W2 ) ;
870 if(numOfTR == 0)
871 {
872 return energyTR ;
873 }
874 else
875 {
876 G4cout<<"It is still OK in GetEnergyTR(int,int,int)"<<G4endl;
877 for(iTR=0;iTR<numOfTR;iTR++)
878 {
879 energyPos = ((*energyVector1)(0)*W1+
880 (*energyVector2)(0)*W2)*G4UniformRand() ;
881 for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
882 {
883 if(energyPos >= ((*energyVector1)(iTransfer)*W1+
884 (*energyVector2)(iTransfer)*W2)) break ;
885 }
886 energyTR += (energyVector1->GetLowEdgeEnergy(iTransfer))*W1+
887 (energyVector2->GetLowEdgeEnergy(iTransfer))*W2 ;
888
889 }
890 }
891 }
892 }
893
894 return energyTR ;
895}
896
897////////////////////////////////////////////////////////////////////////////
898//
899// Test function for checking of PostStepDoIt random preparation of TR photon
900// theta angle relative to particle direction
901//
902
903
904G4double
905G4ForwardXrayTR::GetThetaTR(G4int, G4int, G4int) const
906{
907 G4double theta = 0.0 ;
908
909 return theta ;
910}
911
912
913
914// end of G4ForwardXrayTR implementation file
915//
916///////////////////////////////////////////////////////////////////////////
Note: See TracBrowser for help on using the repository browser.