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

Last change on this file since 1347 was 1337, checked in by garnier, 15 years ago

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