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

Last change on this file since 1199 was 1196, checked in by garnier, 15 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.