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

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

tag geant4.9.4 beta 1 + modifs locales
File size: 31.3 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.15 2010/06/16 15:34:15 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-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
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)
411// varAngle =2* (1 - std::cos(Theta)) or approximately = Theta*Theta
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 ;
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)   ;
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
503// x = 2*(1 - std::cos(Theta)) ~ Theta^2
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)  ;
514  return ( (a + b)*std::log((x + b)/(x + a))/(a - b)
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        }
688        theta = std::sqrt((*fAngleDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1)) ;
689
690        // G4cout<<iTransfer<<" :  theta = "<<theta<<G4endl ;
691
692        phi = twopi*G4UniformRand() ;
693        dirX = std::sin(theta)*std::cos(phi)  ;
694        dirY = std::sin(theta)*std::sin(phi)  ;
695        dirZ = std::cos(theta)           ;
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          }
763          theta = std::sqrt(((*fAngleDistrTable)(iPlace)->
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() ;
771          dirX = std::sin(theta)*std::cos(phi)  ;
772          dirY = std::sin(theta)*std::sin(phi)  ;
773          dirZ = std::cos(theta)           ;
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.