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

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

tag geant4.9.4 beta 1 + modifs locales
File size: 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.