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

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

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