source: trunk/source/processes/electromagnetic/xrays/src/G4VXTRenergyLoss.cc @ 1350

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

tag geant4.9.4 beta 1 + modifs locales
File size: 45.8 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: G4VXTRenergyLoss.cc,v 1.45 2010/06/16 15:34:15 gcosmo Exp $
28// GEANT4 tag $Name: geant4-09-04-beta-01 $
29//
30// History:
31// 2001-2002 R&D by V.Grichine
32// 19.06.03 V. Grichine, modifications in BuildTable for the integration
33//                       in respect of angle: range is increased, accuracy is
34//                       improved
35// 28.07.05, P.Gumplinger add G4ProcessType to constructor
36// 28.09.07, V.Ivanchenko general cleanup without change of algorithms
37//
38
39#include "G4Timer.hh"
40
41#include "G4VXTRenergyLoss.hh"
42#include "G4Poisson.hh"
43#include "G4MaterialTable.hh"
44#include "G4VDiscreteProcess.hh"
45#include "G4VParticleChange.hh"
46#include "G4VSolid.hh"
47
48#include "G4RotationMatrix.hh"
49#include "G4ThreeVector.hh"
50#include "G4AffineTransform.hh"
51#include "G4SandiaTable.hh"
52
53#include "G4PhysicsVector.hh"
54#include "G4PhysicsFreeVector.hh"
55#include "G4PhysicsLinearVector.hh"
56
57////////////////////////////////////////////////////////////////////////////
58//
59// Constructor, destructor
60
61G4VXTRenergyLoss::G4VXTRenergyLoss(G4LogicalVolume *anEnvelope,
62                                 G4Material* foilMat,G4Material* gasMat,
63                                 G4double a, G4double b,
64                                 G4int n,const G4String& processName,
65                                 G4ProcessType type) :
66  G4VDiscreteProcess(processName, type),
67  fGammaCutInKineticEnergy(0),
68  fGammaTkinCut(0),
69  fAngleDistrTable(0),
70  fEnergyDistrTable(0),
71  fPlatePhotoAbsCof(0),
72  fGasPhotoAbsCof(0),
73  fAngleForEnergyTable(0)
74{
75  verboseLevel = 1;
76
77  // Initialization of local constants
78  fTheMinEnergyTR = 1.0*keV;
79  fTheMaxEnergyTR = 100.0*keV;
80  fTheMaxAngle    = 1.0e-3;
81  fTheMinAngle    = 5.0e-6;
82  fBinTR          = 50;
83
84  fMinProtonTkin  = 100.0*GeV;
85  fMaxProtonTkin  = 100.0*TeV;
86  fTotBin         = 50;
87
88  // Proton energy vector initialization
89
90  fProtonEnergyVector = new G4PhysicsLogVector(fMinProtonTkin,
91                                               fMaxProtonTkin,
92                                               fTotBin  );
93
94  fXTREnergyVector = new G4PhysicsLogVector(fTheMinEnergyTR,
95                                            fTheMaxEnergyTR,
96                                            fBinTR  );
97
98  fPlasmaCof = 4.0*pi*fine_structure_const*hbarc*hbarc*hbarc/electron_mass_c2;
99
100  fCofTR     = fine_structure_const/pi;
101
102  fEnvelope  = anEnvelope ;
103
104  fPlateNumber = n ;
105  if(verboseLevel > 0)
106    G4cout<<"### G4VXTRenergyLoss: the number of TR radiator plates = "
107          <<fPlateNumber<<G4endl ;
108  if(fPlateNumber == 0)
109  {
110    G4Exception("G4VXTRenergyLoss: No plates in X-ray TR radiator") ;
111  }
112  // default is XTR dEdx, not flux after radiator
113
114  fExitFlux      = false;
115  fAngleRadDistr = false;
116  fCompton       = false;
117
118  fLambda = DBL_MAX;
119  // Mean thicknesses of plates and gas gaps
120
121  fPlateThick = a ;
122  fGasThick   = b ;
123  fTotalDist  = fPlateNumber*(fPlateThick+fGasThick) ;
124  if(verboseLevel > 0)
125    G4cout<<"total radiator thickness = "<<fTotalDist/cm<<" cm"<<G4endl ;
126
127  // index of plate material
128  fMatIndex1 = foilMat->GetIndex()  ;
129  if(verboseLevel > 0)
130    G4cout<<"plate material = "<<foilMat->GetName()<<G4endl ;
131
132  // index of gas material
133  fMatIndex2 = gasMat->GetIndex()  ;
134  if(verboseLevel > 0)
135    G4cout<<"gas material = "<<gasMat->GetName()<<G4endl ;
136
137  // plasma energy squared for plate material
138
139  fSigma1 = fPlasmaCof*foilMat->GetElectronDensity()  ;
140  //  fSigma1 = (20.9*eV)*(20.9*eV) ;
141  if(verboseLevel > 0)
142    G4cout<<"plate plasma energy = "<<std::sqrt(fSigma1)/eV<<" eV"<<G4endl ;
143
144  // plasma energy squared for gas material
145
146  fSigma2 = fPlasmaCof*gasMat->GetElectronDensity()  ;
147  if(verboseLevel > 0)
148    G4cout<<"gas plasma energy = "<<std::sqrt(fSigma2)/eV<<" eV"<<G4endl ;
149
150  // Compute cofs for preparation of linear photo absorption
151
152  ComputePlatePhotoAbsCof() ;
153  ComputeGasPhotoAbsCof() ;
154
155  pParticleChange = &fParticleChange;
156
157}
158
159///////////////////////////////////////////////////////////////////////////
160
161G4VXTRenergyLoss::~G4VXTRenergyLoss()
162{
163  if(fEnvelope) delete fEnvelope;
164}
165
166///////////////////////////////////////////////////////////////////////////////
167//
168// Returns condition for application of the model depending on particle type
169
170
171G4bool G4VXTRenergyLoss::IsApplicable(const G4ParticleDefinition& particle)
172{
173  return  ( particle.GetPDGCharge() != 0.0 ) ;
174}
175
176/////////////////////////////////////////////////////////////////////////////////
177//
178// Calculate step size for XTR process inside raaditor
179
180G4double G4VXTRenergyLoss::GetMeanFreePath(const G4Track& aTrack,
181                                           G4double, // previousStepSize,
182                                           G4ForceCondition* condition)
183{
184  G4int iTkin, iPlace;
185  G4double lambda, sigma, kinEnergy, mass, gamma;
186  G4double charge, chargeSq, massRatio, TkinScaled;
187  G4double E1,E2,W,W1,W2;
188
189  *condition = NotForced;
190 
191  if( aTrack.GetVolume()->GetLogicalVolume() != fEnvelope ) lambda = DBL_MAX;
192  else
193  {
194    const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
195    kinEnergy = aParticle->GetKineticEnergy();
196    mass      = aParticle->GetDefinition()->GetPDGMass();
197    gamma     = 1.0 + kinEnergy/mass;
198    if(verboseLevel > 1)
199    {
200      G4cout<<" gamma = "<<gamma<<";   fGamma = "<<fGamma<<G4endl;
201    }
202
203    if ( std::fabs( gamma - fGamma ) < 0.05*gamma ) lambda = fLambda;
204    else
205    {
206      charge = aParticle->GetDefinition()->GetPDGCharge();
207      chargeSq  = charge*charge;
208      massRatio = proton_mass_c2/mass;
209      TkinScaled = kinEnergy*massRatio;
210
211      for(iTkin = 0; iTkin < fTotBin; iTkin++)
212      {
213        if( TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin))  break ;   
214      }
215      iPlace = iTkin - 1 ;
216
217      if(iTkin == 0) lambda = DBL_MAX; // Tkin is too small, neglect of TR photon generation
218      else          // general case: Tkin between two vectors of the material
219      {
220        if(iTkin == fTotBin) 
221        {
222          sigma = (*(*fEnergyDistrTable)(iPlace))(0)*chargeSq;
223        }
224        else
225        {
226          E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ; 
227          E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin)     ;
228           W = 1.0/(E2 - E1) ;
229          W1 = (E2 - TkinScaled)*W ;
230          W2 = (TkinScaled - E1)*W ;
231          sigma = ( (*(*fEnergyDistrTable)(iPlace  ))(0)*W1 +
232                (*(*fEnergyDistrTable)(iPlace+1))(0)*W2   )*chargeSq;
233     
234        }
235        if (sigma < DBL_MIN) lambda = DBL_MAX;
236        else                 lambda = 1./sigma; 
237        fLambda = lambda;
238        fGamma  = gamma;   
239        if(verboseLevel > 1)
240        {
241          G4cout<<" lambda = "<<lambda/mm<<" mm"<<G4endl;
242        }
243      }
244    }
245  } 
246  return lambda;
247}
248
249//////////////////////////////////////////////////////////////////////////
250//
251// Interface for build table from physics list
252
253void G4VXTRenergyLoss::BuildPhysicsTable(const G4ParticleDefinition& pd)
254{
255  if(pd.GetPDGCharge()  == 0.) 
256  {
257    G4Exception("G4VXTRenergyLoss::BuildPhysicsTable", "Notification", JustWarning,
258                 "XTR initialisation for neutral particle ?!" );   
259  }
260  BuildTable();
261  if (fAngleRadDistr) 
262  {
263    if(verboseLevel > 0)
264      G4cout<<"Build angle distribution according the transparent regular radiator"
265            <<G4endl;
266    BuildAngleTable();
267  }
268}
269
270
271//////////////////////////////////////////////////////////////////////////
272//
273// Build integral energy distribution of XTR photons
274
275void G4VXTRenergyLoss::BuildTable()
276{
277  G4int iTkin, iTR, iPlace;
278  G4double radiatorCof = 1.0;           // for tuning of XTR yield
279
280  fEnergyDistrTable = new G4PhysicsTable(fTotBin);
281
282  fGammaTkinCut = 0.0;
283 
284  // setting of min/max TR energies
285 
286  if(fGammaTkinCut > fTheMinEnergyTR)  fMinEnergyTR = fGammaTkinCut ;
287  else                                 fMinEnergyTR = fTheMinEnergyTR ;
288       
289  if(fGammaTkinCut > fTheMaxEnergyTR) fMaxEnergyTR = 2.0*fGammaTkinCut ; 
290  else                                fMaxEnergyTR = fTheMaxEnergyTR ;
291
292  G4cout.precision(4) ;
293  G4Timer timer ;
294  timer.Start() ;
295
296  if(verboseLevel > 0) {
297    G4cout<<G4endl;
298    G4cout<<"Lorentz Factor"<<"\t"<<"XTR photon number"<<G4endl;
299    G4cout<<G4endl;
300  }
301  for( iTkin = 0 ; iTkin < fTotBin ; iTkin++ )      // Lorentz factor loop
302  {
303    G4PhysicsLogVector* energyVector = new G4PhysicsLogVector( fMinEnergyTR,
304                                                               fMaxEnergyTR,
305                                                               fBinTR  ) ;
306
307    fGamma = 1.0 + (fProtonEnergyVector->
308                    GetLowEdgeEnergy(iTkin)/proton_mass_c2) ;
309
310    fMaxThetaTR = 25.0/(fGamma*fGamma) ;  // theta^2
311
312    fTheMinAngle = 1.0e-3 ; // was 5.e-6, e-6 !!!, e-5, e-4
313 
314    if( fMaxThetaTR > fTheMaxAngle )    fMaxThetaTR = fTheMaxAngle; 
315    else
316      {
317        if( fMaxThetaTR < fTheMinAngle )  fMaxThetaTR = fTheMinAngle;
318      }
319    G4PhysicsLinearVector* angleVector = new G4PhysicsLinearVector(0.0,
320                                                                   fMaxThetaTR,
321                                                                   fBinTR  );
322
323    G4double energySum = 0.0;
324    G4double angleSum  = 0.0;
325
326    G4Integrator<G4VXTRenergyLoss,G4double(G4VXTRenergyLoss::*)(G4double)> integral;
327
328    energyVector->PutValue(fBinTR-1,energySum);
329    angleVector->PutValue(fBinTR-1,angleSum);
330
331    for( iTR = fBinTR - 2 ; iTR >= 0 ; iTR-- )
332      {
333        energySum += radiatorCof*fCofTR*integral.Legendre10(
334                     this,&G4VXTRenergyLoss::SpectralXTRdEdx,
335                     energyVector->GetLowEdgeEnergy(iTR),
336                     energyVector->GetLowEdgeEnergy(iTR+1) ); 
337
338        //    angleSum  += fCofTR*integral.Legendre96(
339        //       this,&G4VXTRenergyLoss::AngleXTRdEdx,
340        //       angleVector->GetLowEdgeEnergy(iTR),
341        //       angleVector->GetLowEdgeEnergy(iTR+1) );
342
343        energyVector->PutValue(iTR,energySum/fTotalDist);
344        //  angleVector ->PutValue(iTR,angleSum);
345      }
346    if(verboseLevel > 0)
347      {
348        G4cout
349          // <<iTkin<<"\t"
350          //   <<"fGamma = "
351          <<fGamma<<"\t"  //  <<"  fMaxThetaTR = "<<fMaxThetaTR
352          //  <<"sumN = "
353          <<energySum      // <<" ; sumA = "<<angleSum
354          <<G4endl;
355      }
356    iPlace = iTkin;
357    fEnergyDistrTable->insertAt(iPlace,energyVector);
358     //  fAngleDistrTable->insertAt(iPlace,angleVector);
359  }     
360  timer.Stop();
361  G4cout.precision(6);
362  if(verboseLevel > 0) {
363    G4cout<<G4endl;
364    G4cout<<"total time for build X-ray TR energy loss tables = "
365          <<timer.GetUserElapsed()<<" s"<<G4endl;
366  }
367  fGamma = 0.;
368  return ;
369}
370
371//////////////////////////////////////////////////////////////////////////
372//
373//
374
375void G4VXTRenergyLoss::BuildEnergyTable()
376{
377}
378
379////////////////////////////////////////////////////////////////////////
380//
381// Build XTR angular distribution at given energy based on the model
382// of transparent regular radiator
383
384void G4VXTRenergyLoss::BuildAngleTable()
385{
386  G4int iTkin, iTR;
387  G4double  energy;
388
389  fGammaTkinCut = 0.0;
390                             
391  // setting of min/max TR energies
392 
393  if(fGammaTkinCut > fTheMinEnergyTR)  fMinEnergyTR = fGammaTkinCut;
394  else                                 fMinEnergyTR = fTheMinEnergyTR;
395       
396  if(fGammaTkinCut > fTheMaxEnergyTR) fMaxEnergyTR = 2.0*fGammaTkinCut; 
397  else                                fMaxEnergyTR = fTheMaxEnergyTR;
398
399  G4cout.precision(4);
400  G4Timer timer;
401  timer.Start();
402  if(verboseLevel > 0) {
403    G4cout<<G4endl;
404    G4cout<<"Lorentz Factor"<<"\t"<<"XTR photon number"<<G4endl;
405    G4cout<<G4endl;
406  }
407  for( iTkin = 0 ; iTkin < fTotBin ; iTkin++ )      // Lorentz factor loop
408  {
409   
410    fGamma = 1.0 + (fProtonEnergyVector->
411                            GetLowEdgeEnergy(iTkin)/proton_mass_c2) ;
412
413    fMaxThetaTR = 25.0/(fGamma*fGamma) ;  // theta^2
414
415    fTheMinAngle = 1.0e-3 ; // was 5.e-6, e-6 !!!, e-5, e-4
416 
417    if( fMaxThetaTR > fTheMaxAngle )    fMaxThetaTR = fTheMaxAngle; 
418    else
419    {
420       if( fMaxThetaTR < fTheMinAngle )  fMaxThetaTR = fTheMinAngle;
421    }
422
423    fAngleForEnergyTable = new G4PhysicsTable(fBinTR);
424
425    for( iTR = 0; iTR < fBinTR; iTR++ )
426    {
427      // energy = fMinEnergyTR*(iTR+1);
428
429      energy = fXTREnergyVector->GetLowEdgeEnergy(iTR);
430
431      G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(fBinTR);
432
433      angleVector = GetAngleVector(energy,fBinTR);
434      // G4cout<<G4endl;
435
436      fAngleForEnergyTable->insertAt(iTR,angleVector) ;
437    }
438   
439    fAngleBank.push_back(fAngleForEnergyTable); 
440  }     
441  timer.Stop();
442  G4cout.precision(6);
443  if(verboseLevel > 0) {
444    G4cout<<G4endl;
445    G4cout<<"total time for build XTR angle for given energy tables = "
446          <<timer.GetUserElapsed()<<" s"<<G4endl;
447  }
448  fGamma = 0.;
449 
450  return;
451} 
452
453/////////////////////////////////////////////////////////////////////////
454//
455// Vector of angles and angle integral distributions
456
457G4PhysicsFreeVector* G4VXTRenergyLoss::GetAngleVector(G4double energy, G4int n)
458{
459  G4double theta=0., result, tmp=0., cof1, cof2, cofMin, cofPHC, angleSum  = 0.;
460  G4int iTheta, k, kMax, kMin;
461
462  G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(n);
463 
464  cofPHC  = 4*pi*hbarc;
465  tmp     = (fSigma1 - fSigma2)/cofPHC/energy;
466  cof1    = fPlateThick*tmp;
467  cof2    = fGasThick*tmp;
468
469  cofMin  =  energy*(fPlateThick + fGasThick)/fGamma/fGamma;
470  cofMin += (fPlateThick*fSigma1 + fGasThick*fSigma2)/energy;
471  cofMin /= cofPHC;
472
473  kMin = G4int(cofMin);
474  if (cofMin > kMin) kMin++;
475
476  kMax = kMin + fBinTR -1;
477  if(verboseLevel > 2)
478  {
479    G4cout<<"n-1 = "<<n-1<<"; theta = "
480          <<std::sqrt(fMaxThetaTR)*fGamma<<"; tmp = "
481          <<0. 
482          <<";    angleSum = "<<angleSum<<G4endl; 
483  }
484  angleVector->PutValue(n-1,fMaxThetaTR, angleSum);
485
486  for( iTheta = n - 2 ; iTheta >= 1 ; iTheta-- )
487  {
488
489    k = iTheta- 1 + kMin;
490
491    tmp    = pi*fPlateThick*(k + cof2)/(fPlateThick + fGasThick);
492
493    result = (k - cof1)*(k - cof1)*(k + cof2)*(k + cof2);
494
495    tmp = std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
496
497    if( k == kMin && kMin == G4int(cofMin) )
498    {
499      angleSum   += 0.5*tmp; // 0.5*std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
500    }
501    else
502    {
503      angleSum   += tmp; // std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
504    }
505    theta = std::abs(k-cofMin)*cofPHC/energy/(fPlateThick + fGasThick);
506    if(verboseLevel > 2)
507    {
508      G4cout<<"iTheta = "<<iTheta<<"; k = "<<k<<"; theta = "
509            <<std::sqrt(theta)*fGamma<<"; tmp = "
510            <<tmp // std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result
511            <<";    angleSum = "<<angleSum<<G4endl;
512    }
513    angleVector->PutValue( iTheta, theta, angleSum );       
514  }
515  if (theta > 0.)
516  {
517    angleSum += 0.5*tmp;
518    theta = 0.;
519  }
520  if(verboseLevel > 2)
521  {
522    G4cout<<"iTheta = "<<iTheta<<"; theta = "
523          <<std::sqrt(theta)*fGamma<<"; tmp = "
524          <<tmp
525          <<";    angleSum = "<<angleSum<<G4endl;
526  }
527  angleVector->PutValue( iTheta, theta, angleSum );
528
529  return angleVector;
530}
531
532////////////////////////////////////////////////////////////////////////
533//
534// Build XTR angular distribution based on the model of transparent regular radiator
535
536void G4VXTRenergyLoss::BuildGlobalAngleTable()
537{
538  G4int iTkin, iTR, iPlace;
539  G4double radiatorCof = 1.0;           // for tuning of XTR yield
540  G4double angleSum;
541  fAngleDistrTable = new G4PhysicsTable(fTotBin);
542
543  fGammaTkinCut = 0.0;
544 
545  // setting of min/max TR energies
546 
547  if(fGammaTkinCut > fTheMinEnergyTR)  fMinEnergyTR = fGammaTkinCut ;
548  else                                 fMinEnergyTR = fTheMinEnergyTR ;
549       
550  if(fGammaTkinCut > fTheMaxEnergyTR) fMaxEnergyTR = 2.0*fGammaTkinCut ; 
551  else                                fMaxEnergyTR = fTheMaxEnergyTR ;
552
553  G4cout.precision(4) ;
554  G4Timer timer ;
555  timer.Start() ;
556  if(verboseLevel > 0) {
557    G4cout<<G4endl;
558    G4cout<<"Lorentz Factor"<<"\t"<<"XTR photon number"<<G4endl;
559    G4cout<<G4endl;
560  }
561  for( iTkin = 0 ; iTkin < fTotBin ; iTkin++ )      // Lorentz factor loop
562  {
563   
564    fGamma = 1.0 + (fProtonEnergyVector->
565                            GetLowEdgeEnergy(iTkin)/proton_mass_c2) ;
566
567    fMaxThetaTR = 25.0/(fGamma*fGamma) ;  // theta^2
568
569    fTheMinAngle = 1.0e-3 ; // was 5.e-6, e-6 !!!, e-5, e-4
570 
571    if( fMaxThetaTR > fTheMaxAngle )    fMaxThetaTR = fTheMaxAngle; 
572    else
573    {
574       if( fMaxThetaTR < fTheMinAngle )  fMaxThetaTR = fTheMinAngle;
575    }
576    G4PhysicsLinearVector* angleVector = new G4PhysicsLinearVector(0.0,
577                                                               fMaxThetaTR,
578                                                               fBinTR      );
579
580    angleSum  = 0.0;
581
582    G4Integrator<G4VXTRenergyLoss,G4double(G4VXTRenergyLoss::*)(G4double)> integral;
583
584   
585    angleVector->PutValue(fBinTR-1,angleSum);
586
587    for( iTR = fBinTR - 2 ; iTR >= 0 ; iTR-- )
588    {
589
590      angleSum  += radiatorCof*fCofTR*integral.Legendre96(
591                   this,&G4VXTRenergyLoss::AngleXTRdEdx,
592                   angleVector->GetLowEdgeEnergy(iTR),
593                   angleVector->GetLowEdgeEnergy(iTR+1) );
594
595      angleVector ->PutValue(iTR,angleSum);
596    }
597    if(verboseLevel > 1) {
598      G4cout
599        // <<iTkin<<"\t"
600        //   <<"fGamma = "
601        <<fGamma<<"\t"  //  <<"  fMaxThetaTR = "<<fMaxThetaTR
602        //  <<"sumN = "<<energySum      // <<" ; sumA = "
603        <<angleSum
604        <<G4endl;
605    }
606    iPlace = iTkin;
607    fAngleDistrTable->insertAt(iPlace,angleVector);
608  }     
609  timer.Stop();
610  G4cout.precision(6);
611  if(verboseLevel > 0) {
612    G4cout<<G4endl;
613    G4cout<<"total time for build X-ray TR angle tables = "
614          <<timer.GetUserElapsed()<<" s"<<G4endl;
615  }
616  fGamma = 0.;
617 
618  return;
619} 
620
621
622//////////////////////////////////////////////////////////////////////////////
623//
624// The main function which is responsible for the treatment of a particle passage
625// trough G4Envelope with discrete generation of G4Gamma
626
627G4VParticleChange* G4VXTRenergyLoss::PostStepDoIt( const G4Track& aTrack, 
628                                                  const G4Step&  aStep   )
629{
630  G4int iTkin, iPlace;
631  G4double energyTR, theta,theta2, phi, dirX, dirY, dirZ;
632 
633
634  fParticleChange.Initialize(aTrack);
635
636  if(verboseLevel > 1)
637  {
638    G4cout<<"Start of G4VXTRenergyLoss::PostStepDoIt "<<G4endl ;
639    G4cout<<"name of current material =  "
640          <<aTrack.GetVolume()->GetLogicalVolume()->GetMaterial()->GetName()<<G4endl ;
641  }
642  if( aTrack.GetVolume()->GetLogicalVolume() != fEnvelope ) 
643  {
644    if(verboseLevel > 0)
645    {
646      G4cout<<"Go out from G4VXTRenergyLoss::PostStepDoIt: wrong volume "<<G4endl;
647    }
648    return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
649  }
650  else
651  {
652    G4StepPoint* pPostStepPoint        = aStep.GetPostStepPoint();
653    const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
654   
655    // Now we are ready to Generate one TR photon
656
657    G4double kinEnergy = aParticle->GetKineticEnergy() ;
658    G4double mass      = aParticle->GetDefinition()->GetPDGMass() ;
659    G4double gamma     = 1.0 + kinEnergy/mass ;
660
661    if(verboseLevel > 1 )
662    {
663      G4cout<<"gamma = "<<gamma<<G4endl ;
664    }
665    G4double         massRatio   = proton_mass_c2/mass ;
666    G4double          TkinScaled = kinEnergy*massRatio ;
667    G4ThreeVector      position  = pPostStepPoint->GetPosition();
668    G4ParticleMomentum direction = aParticle->GetMomentumDirection();
669    G4double           startTime = pPostStepPoint->GetGlobalTime();
670
671    for( iTkin = 0; iTkin < fTotBin; iTkin++ )
672    {
673      if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin))  break;   
674    }
675    iPlace = iTkin - 1;
676
677    if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
678    {
679      if( verboseLevel > 0)
680      {
681        G4cout<<"Go out from G4VXTRenergyLoss::PostStepDoIt:iTkin = "<<iTkin<<G4endl;
682      }
683      return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
684    } 
685    else          // general case: Tkin between two vectors of the material
686    {
687      fParticleChange.SetNumberOfSecondaries(1);
688
689      energyTR = GetXTRrandomEnergy(TkinScaled,iTkin);
690
691      if( verboseLevel > 1)
692      {
693        G4cout<<"energyTR = "<<energyTR/keV<<" keV"<<G4endl;
694      }
695      if (fAngleRadDistr)
696      {
697        // theta = std::fabs(G4RandGauss::shoot(0.0,pi/gamma));
698        theta2 = GetRandomAngle(energyTR,iTkin);
699        if(theta2 > 0.) theta = std::sqrt(theta2);
700        else            theta = theta2;
701      }
702      else theta = std::fabs(G4RandGauss::shoot(0.0,pi/gamma));
703
704      if( theta >= 0.1 ) theta = 0.1;
705
706      // G4cout<<" : theta = "<<theta<<endl ;
707
708      phi = twopi*G4UniformRand();
709
710      dirX = std::sin(theta)*std::cos(phi);
711      dirY = std::sin(theta)*std::sin(phi);
712      dirZ = std::cos(theta);
713
714      G4ThreeVector directionTR(dirX,dirY,dirZ);
715      directionTR.rotateUz(direction);
716      directionTR.unit();
717
718      G4DynamicParticle* aPhotonTR = new G4DynamicParticle(G4Gamma::Gamma(),
719                                                           directionTR, energyTR);
720
721      // A XTR photon is set on the particle track inside the radiator
722      // and is moved to the G4Envelope surface for standard X-ray TR models
723      // only. The case of fExitFlux=true
724
725      if( fExitFlux )
726      {
727        const G4RotationMatrix* rotM = pPostStepPoint->GetTouchable()->GetRotation();
728        G4ThreeVector transl = pPostStepPoint->GetTouchable()->GetTranslation();
729        G4AffineTransform transform = G4AffineTransform(rotM,transl);
730        transform.Invert();
731        G4ThreeVector localP = transform.TransformPoint(position);
732        G4ThreeVector localV = transform.TransformAxis(directionTR);
733
734        G4double distance = fEnvelope->GetSolid()->DistanceToOut(localP, localV);
735        if(verboseLevel > 1)
736        {
737          G4cout<<"distance to exit = "<<distance/mm<<" mm"<<G4endl;
738        }
739        position         += distance*directionTR;
740        startTime        += distance/c_light;
741      }
742      G4Track* aSecondaryTrack = new G4Track( aPhotonTR, 
743                                                startTime, position );
744      aSecondaryTrack->SetTouchableHandle(
745                         aStep.GetPostStepPoint()->GetTouchableHandle());
746      aSecondaryTrack->SetParentID( aTrack.GetTrackID() );
747
748      fParticleChange.AddSecondary(aSecondaryTrack);
749      fParticleChange.ProposeEnergy(kinEnergy);     
750    }
751  }
752  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
753}
754
755///////////////////////////////////////////////////////////////////////
756//
757// This function returns the spectral and angle density of TR quanta
758// in X-ray energy region generated forward when a relativistic
759// charged particle crosses interface between two materials.
760// The high energy small theta approximation is applied.
761// (matter1 -> matter2, or 2->1)
762// varAngle =2* (1 - std::cos(theta)) or approximately = theta*theta
763//
764
765G4complex G4VXTRenergyLoss::OneInterfaceXTRdEdx( G4double energy,
766                                           G4double gamma,
767                                           G4double varAngle ) 
768{
769  G4complex Z1    = GetPlateComplexFZ(energy,gamma,varAngle) ;
770  G4complex Z2    = GetGasComplexFZ(energy,gamma,varAngle) ;
771
772  G4complex zOut  = (Z1 - Z2)*(Z1 - Z2)
773                    * (varAngle*energy/hbarc/hbarc) ; 
774  return    zOut  ;
775
776}
777
778
779//////////////////////////////////////////////////////////////////////////////
780//
781// For photon energy distribution tables. Integrate first over angle
782//
783
784G4double G4VXTRenergyLoss::SpectralAngleXTRdEdx(G4double varAngle)
785{
786  G4double result =  GetStackFactor(fEnergy,fGamma,varAngle);
787  if(result < 0.0) result = 0.0;
788  return result;
789}
790
791/////////////////////////////////////////////////////////////////////////
792//
793// For second integration over energy
794 
795G4double G4VXTRenergyLoss::SpectralXTRdEdx(G4double energy)
796{
797  G4int i, iMax = 8;
798  G4double result = 0.0;
799
800  G4double lim[8] = { 0.0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0 };
801
802  for( i = 0; i < iMax; i++ ) lim[i] *= fMaxThetaTR;
803
804  G4Integrator<G4VXTRenergyLoss,G4double(G4VXTRenergyLoss::*)(G4double)> integral;
805
806  fEnergy = energy;
807
808  for( i = 0; i < iMax-1; i++ )
809  {
810    result += integral.Legendre96(this,&G4VXTRenergyLoss::SpectralAngleXTRdEdx,
811                           lim[i],lim[i+1]);
812    // result += integral.Legendre10(this,&G4VXTRenergyLoss::SpectralAngleXTRdEdx,
813    //                     lim[i],lim[i+1]);
814  }
815
816  return result;
817} 
818 
819//////////////////////////////////////////////////////////////////////////
820//
821// for photon angle distribution tables
822//
823
824G4double G4VXTRenergyLoss::AngleSpectralXTRdEdx(G4double energy)
825{
826  G4double result =  GetStackFactor(energy,fGamma,fVarAngle);
827  if(result < 0) result = 0.0;
828  return result;
829} 
830
831///////////////////////////////////////////////////////////////////////////
832//
833// The XTR angular distribution based on transparent regular radiator
834
835G4double G4VXTRenergyLoss::AngleXTRdEdx(G4double varAngle) 
836{
837  // G4cout<<"angle2 = "<<varAngle<<"; fGamma = "<<fGamma<<G4endl;
838 
839  G4double result;
840  G4double sum = 0., tmp1, tmp2, tmp=0., cof1, cof2, cofMin, cofPHC, energy1, energy2;
841  G4int k, kMax, kMin, i;
842
843  cofPHC  = twopi*hbarc;
844
845  cof1    = (fPlateThick + fGasThick)*(1./fGamma/fGamma + varAngle);
846  cof2    = fPlateThick*fSigma1 + fGasThick*fSigma2;
847
848  // G4cout<<"cof1 = "<<cof1<<"; cof2 = "<<cof2<<"; cofPHC = "<<cofPHC<<G4endl;
849
850  cofMin  =  std::sqrt(cof1*cof2); 
851  cofMin /= cofPHC;
852
853  kMin = G4int(cofMin);
854  if (cofMin > kMin) kMin++;
855
856  kMax = kMin + 9; //  9; // kMin + G4int(tmp);
857
858  // G4cout<<"cofMin = "<<cofMin<<"; kMin = "<<kMin<<"; kMax = "<<kMax<<G4endl;
859
860  for( k = kMin; k <= kMax; k++ )
861  {
862    tmp1 = cofPHC*k;
863    tmp2 = std::sqrt(tmp1*tmp1-cof1*cof2);
864    energy1 = (tmp1+tmp2)/cof1;
865    energy2 = (tmp1-tmp2)/cof1;
866
867    for(i = 0; i < 2; i++)
868    {
869      if( i == 0 )
870      {
871        if (energy1 > fTheMaxEnergyTR || energy1 < fTheMinEnergyTR) continue;
872        tmp1 = ( energy1*energy1*(1./fGamma/fGamma + varAngle) + fSigma1 )
873          * fPlateThick/(4*hbarc*energy1);
874        tmp2 = std::sin(tmp1);
875        tmp  = energy1*tmp2*tmp2;
876        tmp2 = fPlateThick/(4*tmp1);
877        tmp1 = hbarc*energy1/( energy1*energy1*(1./fGamma/fGamma + varAngle) + fSigma2 );
878        tmp *= (tmp1-tmp2)*(tmp1-tmp2);
879        tmp1 = cof1/(4*hbarc) - cof2/(4*hbarc*energy1*energy1);
880        tmp2 = std::abs(tmp1);
881        if(tmp2 > 0.) tmp /= tmp2;
882        else continue;
883      }
884      else
885      {
886        if (energy2 > fTheMaxEnergyTR || energy2 < fTheMinEnergyTR) continue;
887        tmp1 = ( energy2*energy2*(1./fGamma/fGamma + varAngle) + fSigma1 )
888          * fPlateThick/(4*hbarc*energy2);
889        tmp2 = std::sin(tmp1);
890        tmp  = energy2*tmp2*tmp2;
891        tmp2 = fPlateThick/(4*tmp1);
892        tmp1 = hbarc*energy2/( energy2*energy2*(1./fGamma/fGamma + varAngle) + fSigma2 );
893        tmp *= (tmp1-tmp2)*(tmp1-tmp2);
894        tmp1 = cof1/(4*hbarc) - cof2/(4*hbarc*energy2*energy2);
895        tmp2 = std::abs(tmp1);
896        if(tmp2 > 0.) tmp /= tmp2;
897        else continue;
898      }
899      sum += tmp;
900    }
901    // G4cout<<"k = "<<k<<"; energy1 = "<<energy1/keV<<" keV; energy2 = "<<energy2/keV
902    //  <<" keV; tmp = "<<tmp<<"; sum = "<<sum<<G4endl;
903  }
904  result = 4.*pi*fPlateNumber*sum*varAngle;
905  result /= hbarc*hbarc;
906
907  // old code based on general numeric integration
908  // fVarAngle = varAngle;
909  // G4Integrator<G4VXTRenergyLoss,G4double(G4VXTRenergyLoss::*)(G4double)> integral;
910  // result = integral.Legendre10(this,&G4VXTRenergyLoss::AngleSpectralXTRdEdx,
911  //                         fMinEnergyTR,fMaxEnergyTR);
912  return result;
913}
914
915
916//////////////////////////////////////////////////////////////////////
917//////////////////////////////////////////////////////////////////////
918//////////////////////////////////////////////////////////////////////
919//
920// Calculates formation zone for plates. Omega is energy !!!
921
922G4double G4VXTRenergyLoss::GetPlateFormationZone( G4double omega ,
923                                                G4double gamma ,
924                                                G4double varAngle    ) 
925{
926  G4double cof, lambda ;
927  lambda = 1.0/gamma/gamma + varAngle + fSigma1/omega/omega ;
928  cof = 2.0*hbarc/omega/lambda ;
929  return cof ;
930}
931
932//////////////////////////////////////////////////////////////////////
933//
934// Calculates complex formation zone for plates. Omega is energy !!!
935
936G4complex G4VXTRenergyLoss::GetPlateComplexFZ( G4double omega ,
937                                             G4double gamma ,
938                                             G4double varAngle    ) 
939{
940  G4double cof, length,delta, real_v, image_v ;
941
942  length = 0.5*GetPlateFormationZone(omega,gamma,varAngle) ;
943  delta  = length*GetPlateLinearPhotoAbs(omega) ;
944  cof    = 1.0/(1.0 + delta*delta) ;
945
946  real_v  = length*cof ;
947  image_v = real_v*delta ;
948
949  G4complex zone(real_v,image_v); 
950  return zone ;
951}
952
953////////////////////////////////////////////////////////////////////////
954//
955// Computes matrix of Sandia photo absorption cross section coefficients for
956// plate material
957
958void G4VXTRenergyLoss::ComputePlatePhotoAbsCof() 
959{
960  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
961  const G4Material* mat = (*theMaterialTable)[fMatIndex1];
962  fPlatePhotoAbsCof = mat->GetSandiaTable();
963
964  return;
965}
966
967
968
969//////////////////////////////////////////////////////////////////////
970//
971// Returns the value of linear photo absorption coefficient (in reciprocal
972// length) for plate for given energy of X-ray photon omega
973
974G4double G4VXTRenergyLoss::GetPlateLinearPhotoAbs(G4double omega) 
975{
976  //  G4int i ;
977  G4double omega2, omega3, omega4 ; 
978
979  omega2 = omega*omega ;
980  omega3 = omega2*omega ;
981  omega4 = omega2*omega2 ;
982
983  G4double* SandiaCof = fPlatePhotoAbsCof->GetSandiaCofForMaterial(omega);
984  G4double cross = SandiaCof[0]/omega  + SandiaCof[1]/omega2 +
985                   SandiaCof[2]/omega3 + SandiaCof[3]/omega4;
986  return cross;
987}
988
989
990//////////////////////////////////////////////////////////////////////
991//
992// Calculates formation zone for gas. Omega is energy !!!
993
994G4double G4VXTRenergyLoss::GetGasFormationZone( G4double omega ,
995                                              G4double gamma ,
996                                              G4double varAngle   ) 
997{
998  G4double cof, lambda ;
999  lambda = 1.0/gamma/gamma + varAngle + fSigma2/omega/omega ;
1000  cof = 2.0*hbarc/omega/lambda ;
1001  return cof ;
1002}
1003
1004
1005//////////////////////////////////////////////////////////////////////
1006//
1007// Calculates complex formation zone for gas gaps. Omega is energy !!!
1008
1009G4complex G4VXTRenergyLoss::GetGasComplexFZ( G4double omega ,
1010                                           G4double gamma ,
1011                                           G4double varAngle    ) 
1012{
1013  G4double cof, length,delta, real_v, image_v ;
1014
1015  length = 0.5*GetGasFormationZone(omega,gamma,varAngle) ;
1016  delta  = length*GetGasLinearPhotoAbs(omega) ;
1017  cof    = 1.0/(1.0 + delta*delta) ;
1018
1019  real_v   = length*cof ;
1020  image_v  = real_v*delta ;
1021
1022  G4complex zone(real_v,image_v); 
1023  return zone ;
1024}
1025
1026
1027
1028////////////////////////////////////////////////////////////////////////
1029//
1030// Computes matrix of Sandia photo absorption cross section coefficients for
1031// gas material
1032
1033void G4VXTRenergyLoss::ComputeGasPhotoAbsCof() 
1034{
1035  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1036  const G4Material* mat = (*theMaterialTable)[fMatIndex2];
1037  fGasPhotoAbsCof = mat->GetSandiaTable();
1038  return;
1039}
1040
1041//////////////////////////////////////////////////////////////////////
1042//
1043// Returns the value of linear photo absorption coefficient (in reciprocal
1044// length) for gas
1045
1046G4double G4VXTRenergyLoss::GetGasLinearPhotoAbs(G4double omega) 
1047{
1048  G4double omega2, omega3, omega4 ; 
1049
1050  omega2 = omega*omega ;
1051  omega3 = omega2*omega ;
1052  omega4 = omega2*omega2 ;
1053
1054  G4double* SandiaCof = fGasPhotoAbsCof->GetSandiaCofForMaterial(omega);
1055  G4double cross = SandiaCof[0]/omega  + SandiaCof[1]/omega2 +
1056                   SandiaCof[2]/omega3 + SandiaCof[3]/omega4;
1057  return cross;
1058
1059}
1060
1061//////////////////////////////////////////////////////////////////////
1062//
1063// Calculates the product of linear cof by formation zone for plate.
1064// Omega is energy !!!
1065
1066G4double G4VXTRenergyLoss::GetPlateZmuProduct( G4double omega ,
1067                                             G4double gamma ,
1068                                             G4double varAngle   ) 
1069{
1070  return GetPlateFormationZone(omega,gamma,varAngle)
1071    * GetPlateLinearPhotoAbs(omega) ;
1072}
1073//////////////////////////////////////////////////////////////////////
1074//
1075// Calculates the product of linear cof by formation zone for plate.
1076// G4cout and output in file in some energy range.
1077
1078void G4VXTRenergyLoss::GetPlateZmuProduct() 
1079{
1080  std::ofstream outPlate("plateZmu.dat", std::ios::out ) ;
1081  outPlate.setf( std::ios::scientific, std::ios::floatfield );
1082
1083  G4int i ;
1084  G4double omega, varAngle, gamma ;
1085  gamma = 10000. ;
1086  varAngle = 1/gamma/gamma ;
1087  if(verboseLevel > 0)
1088    G4cout<<"energy, keV"<<"\t"<<"Zmu for plate"<<G4endl ;
1089  for(i=0;i<100;i++)
1090  {
1091    omega = (1.0 + i)*keV ;
1092    if(verboseLevel > 1)
1093      G4cout<<omega/keV<<"\t"<<GetPlateZmuProduct(omega,gamma,varAngle)<<"\t";
1094    if(verboseLevel > 0)
1095      outPlate<<omega/keV<<"\t\t"<<GetPlateZmuProduct(omega,gamma,varAngle)<<G4endl ;
1096  }
1097  return  ;
1098}
1099
1100//////////////////////////////////////////////////////////////////////
1101//
1102// Calculates the product of linear cof by formation zone for gas.
1103// Omega is energy !!!
1104
1105G4double G4VXTRenergyLoss::GetGasZmuProduct( G4double omega ,
1106                                             G4double gamma ,
1107                                             G4double varAngle   ) 
1108{
1109  return GetGasFormationZone(omega,gamma,varAngle)*GetGasLinearPhotoAbs(omega) ;
1110}
1111//////////////////////////////////////////////////////////////////////
1112//
1113// Calculates the product of linear cof byformation zone for gas.
1114// G4cout and output in file in some energy range.
1115
1116void G4VXTRenergyLoss::GetGasZmuProduct() 
1117{
1118  std::ofstream outGas("gasZmu.dat", std::ios::out ) ;
1119  outGas.setf( std::ios::scientific, std::ios::floatfield );
1120  G4int i ;
1121  G4double omega, varAngle, gamma ;
1122  gamma = 10000. ;
1123  varAngle = 1/gamma/gamma ;
1124  if(verboseLevel > 0)
1125    G4cout<<"energy, keV"<<"\t"<<"Zmu for gas"<<G4endl ;
1126  for(i=0;i<100;i++)
1127  {
1128    omega = (1.0 + i)*keV ;
1129    if(verboseLevel > 1)
1130      G4cout<<omega/keV<<"\t"<<GetGasZmuProduct(omega,gamma,varAngle)<<"\t" ;
1131    if(verboseLevel > 0)
1132      outGas<<omega/keV<<"\t\t"<<GetGasZmuProduct(omega,gamma,varAngle)<<G4endl ;
1133  }
1134  return  ;
1135}
1136
1137////////////////////////////////////////////////////////////////////////
1138//
1139// Computes Compton cross section for plate material in 1/mm
1140
1141G4double G4VXTRenergyLoss::GetPlateCompton(G4double omega) 
1142{
1143  G4int i, numberOfElements;
1144  G4double xSection = 0., nowZ, sumZ = 0.;
1145
1146  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1147  numberOfElements = (*theMaterialTable)[fMatIndex1]->GetNumberOfElements() ;
1148
1149  for( i = 0; i < numberOfElements; i++ )
1150  {
1151    nowZ      = (*theMaterialTable)[fMatIndex1]->GetElement(i)->GetZ();
1152    sumZ     += nowZ;
1153    xSection += GetComptonPerAtom(omega,nowZ); // *nowZ;
1154  }
1155  xSection /= sumZ;
1156  xSection *= (*theMaterialTable)[fMatIndex1]->GetElectronDensity();
1157  return xSection;
1158}
1159
1160
1161////////////////////////////////////////////////////////////////////////
1162//
1163// Computes Compton cross section for gas material in 1/mm
1164
1165G4double G4VXTRenergyLoss::GetGasCompton(G4double omega) 
1166{
1167  G4int i, numberOfElements;
1168  G4double xSection = 0., nowZ, sumZ = 0.;
1169
1170  const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
1171  numberOfElements = (*theMaterialTable)[fMatIndex2]->GetNumberOfElements() ;
1172
1173  for( i = 0; i < numberOfElements; i++ )
1174  {
1175    nowZ      = (*theMaterialTable)[fMatIndex2]->GetElement(i)->GetZ();
1176    sumZ     += nowZ;
1177    xSection += GetComptonPerAtom(omega,nowZ); // *nowZ;
1178  }
1179  xSection /= sumZ;
1180  xSection *= (*theMaterialTable)[fMatIndex2]->GetElectronDensity();
1181  return xSection;
1182}
1183
1184////////////////////////////////////////////////////////////////////////
1185//
1186// Computes Compton cross section per atom with Z electrons for gamma with
1187// the energy GammaEnergy
1188
1189G4double G4VXTRenergyLoss::GetComptonPerAtom(G4double GammaEnergy, G4double Z) 
1190{
1191  G4double CrossSection = 0.0 ;
1192  if ( Z < 0.9999 )                 return CrossSection;
1193  if ( GammaEnergy < 0.1*keV      ) return CrossSection;
1194  if ( GammaEnergy > (100.*GeV/Z) ) return CrossSection;
1195
1196  static const G4double a = 20.0 , b = 230.0 , c = 440.0;
1197
1198  static const G4double
1199  d1= 2.7965e-1*barn, d2=-1.8300e-1*barn, d3= 6.7527   *barn, d4=-1.9798e+1*barn,
1200  e1= 1.9756e-5*barn, e2=-1.0205e-2*barn, e3=-7.3913e-2*barn, e4= 2.7079e-2*barn,
1201  f1=-3.9178e-7*barn, f2= 6.8241e-5*barn, f3= 6.0480e-5*barn, f4= 3.0274e-4*barn;
1202
1203  G4double p1Z = Z*(d1 + e1*Z + f1*Z*Z), p2Z = Z*(d2 + e2*Z + f2*Z*Z),
1204           p3Z = Z*(d3 + e3*Z + f3*Z*Z), p4Z = Z*(d4 + e4*Z + f4*Z*Z);
1205
1206  G4double T0  = 15.0*keV;
1207  if (Z < 1.5) T0 = 40.0*keV;
1208
1209  G4double X   = std::max(GammaEnergy, T0) / electron_mass_c2;
1210  CrossSection = p1Z*std::log(1.+2.*X)/X
1211               + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
1212
1213  //  modification for low energy. (special case for Hydrogen)
1214
1215  if (GammaEnergy < T0) 
1216  {
1217    G4double dT0 = 1.*keV;
1218    X = (T0+dT0) / electron_mass_c2 ;
1219    G4double sigma = p1Z*std::log(1.+2*X)/X
1220                    + (p2Z + p3Z*X + p4Z*X*X)/(1. + a*X + b*X*X + c*X*X*X);
1221    G4double   c1 = -T0*(sigma-CrossSection)/(CrossSection*dT0);
1222    G4double   c2 = 0.150;
1223    if (Z > 1.5) c2 = 0.375-0.0556*std::log(Z);
1224    G4double    y = std::log(GammaEnergy/T0);
1225    CrossSection *= std::exp(-y*(c1+c2*y));
1226  }
1227  //  G4cout << "e= " << GammaEnergy << " Z= " << Z << " cross= " << CrossSection << G4endl;
1228  return CrossSection; 
1229}
1230
1231
1232///////////////////////////////////////////////////////////////////////
1233//
1234// This function returns the spectral and angle density of TR quanta
1235// in X-ray energy region generated forward when a relativistic
1236// charged particle crosses interface between two materials.
1237// The high energy small theta approximation is applied.
1238// (matter1 -> matter2, or 2->1)
1239// varAngle =2* (1 - std::cos(theta)) or approximately = theta*theta
1240//
1241
1242G4double
1243G4VXTRenergyLoss::OneBoundaryXTRNdensity( G4double energy,G4double gamma,
1244                                         G4double varAngle ) const
1245{
1246  G4double  formationLength1, formationLength2 ;
1247  formationLength1 = 1.0/
1248  (1.0/(gamma*gamma)
1249  + fSigma1/(energy*energy)
1250  + varAngle) ;
1251  formationLength2 = 1.0/
1252  (1.0/(gamma*gamma)
1253  + fSigma2/(energy*energy)
1254  + varAngle) ;
1255  return (varAngle/energy)*(formationLength1 - formationLength2)
1256              *(formationLength1 - formationLength2)  ;
1257
1258}
1259
1260G4double G4VXTRenergyLoss::GetStackFactor( G4double energy, G4double gamma,
1261                                                     G4double varAngle )
1262{
1263  // return stack factor corresponding to one interface
1264
1265  return std::real( OneInterfaceXTRdEdx(energy,gamma,varAngle) );
1266}
1267
1268//////////////////////////////////////////////////////////////////////////////
1269//
1270// For photon energy distribution tables. Integrate first over angle
1271//
1272
1273G4double G4VXTRenergyLoss::XTRNSpectralAngleDensity(G4double varAngle)
1274{
1275  return OneBoundaryXTRNdensity(fEnergy,fGamma,varAngle)*
1276         GetStackFactor(fEnergy,fGamma,varAngle)             ;
1277}
1278
1279/////////////////////////////////////////////////////////////////////////
1280//
1281// For second integration over energy
1282 
1283G4double G4VXTRenergyLoss::XTRNSpectralDensity(G4double energy)
1284{
1285  fEnergy = energy ;
1286  G4Integrator<G4VXTRenergyLoss,G4double(G4VXTRenergyLoss::*)(G4double)> integral ;
1287  return integral.Legendre96(this,&G4VXTRenergyLoss::XTRNSpectralAngleDensity,
1288                             0.0,0.2*fMaxThetaTR) +
1289         integral.Legendre10(this,&G4VXTRenergyLoss::XTRNSpectralAngleDensity,
1290                             0.2*fMaxThetaTR,fMaxThetaTR) ;
1291} 
1292 
1293//////////////////////////////////////////////////////////////////////////
1294//
1295// for photon angle distribution tables
1296//
1297
1298G4double G4VXTRenergyLoss::XTRNAngleSpectralDensity(G4double energy)
1299{
1300  return OneBoundaryXTRNdensity(energy,fGamma,fVarAngle)*
1301         GetStackFactor(energy,fGamma,fVarAngle)             ;
1302} 
1303
1304///////////////////////////////////////////////////////////////////////////
1305//
1306//
1307
1308G4double G4VXTRenergyLoss::XTRNAngleDensity(G4double varAngle) 
1309{
1310  fVarAngle = varAngle ;
1311  G4Integrator<G4VXTRenergyLoss,G4double(G4VXTRenergyLoss::*)(G4double)> integral ;
1312  return integral.Legendre96(this,&G4VXTRenergyLoss::XTRNAngleSpectralDensity,
1313                             fMinEnergyTR,fMaxEnergyTR) ;
1314}
1315
1316//////////////////////////////////////////////////////////////////////////////
1317//
1318// Check number of photons for a range of Lorentz factors from both energy
1319// and angular tables
1320
1321void G4VXTRenergyLoss::GetNumberOfPhotons()
1322{
1323  G4int iTkin ;
1324  G4double gamma, numberE ;
1325
1326  std::ofstream outEn("numberE.dat", std::ios::out ) ;
1327  outEn.setf( std::ios::scientific, std::ios::floatfield );
1328
1329  std::ofstream outAng("numberAng.dat", std::ios::out ) ;
1330  outAng.setf( std::ios::scientific, std::ios::floatfield );
1331
1332  for(iTkin=0;iTkin<fTotBin;iTkin++)      // Lorentz factor loop
1333  {
1334     gamma = 1.0 + (fProtonEnergyVector->
1335                            GetLowEdgeEnergy(iTkin)/proton_mass_c2) ;
1336     numberE = (*(*fEnergyDistrTable)(iTkin))(0) ;
1337     //  numberA = (*(*fAngleDistrTable)(iTkin))(0) ;
1338     if(verboseLevel > 1)
1339       G4cout<<gamma<<"\t\t"<<numberE<<"\t"    //  <<numberA
1340             <<G4endl ; 
1341     if(verboseLevel > 0)
1342       outEn<<gamma<<"\t\t"<<numberE<<G4endl ; 
1343  }
1344  return ;
1345} 
1346
1347/////////////////////////////////////////////////////////////////////////
1348//
1349// Returns randon energy of a X-ray TR photon for given scaled kinetic energy
1350// of a charged particle
1351
1352G4double G4VXTRenergyLoss::GetXTRrandomEnergy( G4double scaledTkin, G4int iTkin )
1353{
1354  G4int iTransfer, iPlace  ;
1355  G4double transfer = 0.0, position, E1, E2, W1, W2, W ;
1356
1357  iPlace = iTkin - 1 ;
1358
1359  //  G4cout<<"iPlace = "<<iPlace<<endl ;
1360
1361  if(iTkin == fTotBin) // relativistic plato, try from left
1362  {
1363      position = (*(*fEnergyDistrTable)(iPlace))(0)*G4UniformRand() ;
1364
1365      for(iTransfer=0;;iTransfer++)
1366      {
1367        if(position >= (*(*fEnergyDistrTable)(iPlace))(iTransfer)) break ;
1368      }
1369      transfer = GetXTRenergy(iPlace,position,iTransfer);
1370  }
1371  else
1372  {
1373    E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ; 
1374    E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin)     ;
1375    W  = 1.0/(E2 - E1) ;
1376    W1 = (E2 - scaledTkin)*W ;
1377    W2 = (scaledTkin - E1)*W ;
1378
1379    position =( (*(*fEnergyDistrTable)(iPlace))(0)*W1 + 
1380                    (*(*fEnergyDistrTable)(iPlace+1))(0)*W2 )*G4UniformRand() ;
1381
1382        // G4cout<<position<<"\t" ;
1383
1384    for(iTransfer=0;;iTransfer++)
1385    {
1386          if( position >=
1387          ( (*(*fEnergyDistrTable)(iPlace))(iTransfer)*W1 + 
1388            (*(*fEnergyDistrTable)(iPlace+1))(iTransfer)*W2) ) break ;
1389    }
1390    transfer = GetXTRenergy(iPlace,position,iTransfer);
1391   
1392  } 
1393  //  G4cout<<"XTR transfer = "<<transfer/keV<<" keV"<<endl ;
1394  if(transfer < 0.0 ) transfer = 0.0 ;
1395  return transfer ;
1396}
1397
1398////////////////////////////////////////////////////////////////////////
1399//
1400// Returns approximate position of X-ray photon energy during random sampling
1401// over integral energy distribution
1402
1403G4double G4VXTRenergyLoss::GetXTRenergy( G4int    iPlace, 
1404                                       G4double position, 
1405                                       G4int    iTransfer )
1406{
1407  G4double x1, x2, y1, y2, result ;
1408
1409  if(iTransfer == 0)
1410  {
1411    result = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
1412  } 
1413  else
1414  {
1415    y1 = (*(*fEnergyDistrTable)(iPlace))(iTransfer-1) ;
1416    y2 = (*(*fEnergyDistrTable)(iPlace))(iTransfer) ;
1417
1418    x1 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
1419    x2 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
1420
1421    if ( x1 == x2 )    result = x2 ;
1422    else
1423    {
1424      if ( y1 == y2  ) result = x1 + (x2 - x1)*G4UniformRand() ;
1425      else
1426      {
1427        result = x1 + (position - y1)*(x2 - x1)/(y2 - y1) ;
1428      }
1429    }
1430  }
1431  return result ;
1432}
1433
1434/////////////////////////////////////////////////////////////////////////
1435//
1436//  Get XTR photon angle at given energy and Tkin
1437
1438G4double G4VXTRenergyLoss::GetRandomAngle( G4double energyXTR, G4int iTkin )
1439{
1440  G4int iTR, iAngle;
1441  G4double position, angle;
1442
1443  if (iTkin == fTotBin) iTkin--;
1444
1445  fAngleForEnergyTable = fAngleBank[iTkin];
1446
1447  for( iTR = 0; iTR < fBinTR; iTR++ )
1448  {
1449    if( energyXTR < fXTREnergyVector->GetLowEdgeEnergy(iTR) )  break;   
1450  }
1451  if (iTR == fBinTR) iTR--;
1452     
1453  position = (*(*fAngleForEnergyTable)(iTR))(0)*G4UniformRand() ;
1454
1455  for(iAngle = 0;;iAngle++)
1456  {
1457    if(position >= (*(*fAngleForEnergyTable)(iTR))(iAngle)) break ;
1458  }
1459  angle = GetAngleXTR(iTR,position,iAngle);
1460  return angle;
1461}
1462
1463////////////////////////////////////////////////////////////////////////
1464//
1465// Returns approximate position of X-ray photon angle at given energy during random sampling
1466// over integral energy distribution
1467
1468G4double G4VXTRenergyLoss::GetAngleXTR( G4int    iPlace, 
1469                                       G4double position, 
1470                                       G4int    iTransfer )
1471{
1472  G4double x1, x2, y1, y2, result ;
1473
1474  if(iTransfer == 0)
1475  {
1476    result = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
1477  } 
1478  else
1479  {
1480    y1 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer-1) ;
1481    y2 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer) ;
1482
1483    x1 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ;
1484    x2 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ;
1485
1486    if ( x1 == x2 )    result = x2 ;
1487    else
1488    {
1489      if ( y1 == y2  ) result = x1 + (x2 - x1)*G4UniformRand() ;
1490      else
1491      {
1492        result = x1 + (position - y1)*(x2 - x1)/(y2 - y1) ;
1493      }
1494    }
1495  }
1496  return result ;
1497}
1498
1499
1500//
1501//
1502///////////////////////////////////////////////////////////////////////
1503
Note: See TracBrowser for help on using the repository browser.