Geant4_10
G4ForwardXrayTR.cc
Go to the documentation of this file.
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4ForwardXrayTR.cc 68037 2013-03-13 14:15:08Z gcosmo $
28 //
29 // G4ForwardXrayTR class -- implementation file
30 //
31 // History:
32 // 1st version 11.09.97 V. Grichine (Vladimir.Grichine@cern.ch )
33 // 2nd version 17.12.97 V. Grichine
34 // 17-09-01, migration of Materials to pure STL (mma)
35 // 10-03-03, migration to "cut per region" (V.Ivanchenko)
36 // 03.06.03, V.Ivanchenko fix compilation warnings
37 
38 #include "G4ForwardXrayTR.hh"
39 
40 #include "globals.hh"
41 #include "G4PhysicalConstants.hh"
42 #include "G4SystemOfUnits.hh"
43 #include "G4Poisson.hh"
44 #include "G4Material.hh"
45 #include "G4PhysicsTable.hh"
46 #include "G4PhysicsVector.hh"
47 #include "G4PhysicsLinearVector.hh"
48 #include "G4PhysicsLogVector.hh"
49 #include "G4ProductionCutsTable.hh"
50 #include "G4GeometryTolerance.hh"
51 
52 // Initialization of local constants
53 
55 
61 
65 
68 
70 
71 
73 //
74 // Constructor for creation of physics tables (angle and energy TR
75 // distributions) for a couple of selected materials.
76 //
77 // Recommended for use in applications with many materials involved,
78 // when only few (usually couple) materials are interested for generation
79 // of TR on the interface between them
80 
81 
83 G4ForwardXrayTR( const G4String& matName1, // G4Material* pMat1,
84  const G4String& matName2, // G4Material* pMat2,
85  const G4String& processName )
86  : G4TransitionRadiation(processName)
87 {
88  fPtrGamma = 0;
91  fSigma1 = fSigma2 = 0.0;
92  fAngleDistrTable = 0;
94  fMatIndex1 = fMatIndex2 = 0;
95 
96  // Proton energy vector initialization
97  //
100  G4int iMat;
101  const G4ProductionCutsTable* theCoupleTable=
103  G4int numOfCouples = theCoupleTable->GetTableSize();
104 
105  G4bool build = true;
106 
107  for(iMat=0;iMat<numOfCouples;iMat++) // check first material name
108  {
109  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(iMat);
110  if( matName1 == couple->GetMaterial()->GetName() )
111  {
112  fMatIndex1 = couple->GetIndex();
113  break;
114  }
115  }
116  if(iMat == numOfCouples)
117  {
118  G4Exception("G4ForwardXrayTR::G4ForwardXrayTR", "ForwardXrayTR01", JustWarning, "Invalid first material name in G4ForwardXrayTR constructor!");
119  build = false;
120  }
121 
122  if(build) {
123  for(iMat=0;iMat<numOfCouples;iMat++) // check second material name
124  {
125  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(iMat);
126  if( matName2 == couple->GetMaterial()->GetName() )
127  {
128  fMatIndex2 = couple->GetIndex();
129  break;
130  }
131  }
132  if(iMat == numOfCouples)
133  {
134  G4Exception("G4ForwardXrayTR::G4ForwardXrayTR", "ForwardXrayTR02", JustWarning, "Invalid second material name in G4ForwardXrayTR constructor!");
135  build = false;
136  }
137  }
138  // G4cout<<"G4ForwardXray constructor is called"<<G4endl;
139  if(build) { BuildXrayTRtables(); }
140 }
141 
143 //
144 // Constructor used by X-ray transition radiation parametrisation models
145 
147 G4ForwardXrayTR( const G4String& processName )
148  : G4TransitionRadiation(processName)
149 {
150  fPtrGamma = 0;
153  fSigma1 = fSigma2 = 0.0;
154  fAngleDistrTable = 0;
155  fEnergyDistrTable = 0;
156  fMatIndex1 = fMatIndex2 = 0;
157 
158  // Proton energy vector initialization
159  //
162 }
163 
164 
166 //
167 // Destructor
168 //
169 
171 {
172  delete fAngleDistrTable;
173  delete fEnergyDistrTable;
174  delete fProtonEnergyVector;
175 }
176 
178  G4double,
180 {
181  *condition = Forced;
182  return DBL_MAX; // so TR doesn't limit mean free path
183 }
184 
186 //
187 // Build physics tables for energy and angular distributions of X-ray TR photon
188 
190 {
191  G4int iMat, jMat, iTkin, iTR, iPlace;
192  const G4ProductionCutsTable* theCoupleTable=
194  G4int numOfCouples = theCoupleTable->GetTableSize();
195 
197 
200 
201 
202  for(iMat=0;iMat<numOfCouples;iMat++) // loop over pairs of different materials
203  {
204  if( iMat != fMatIndex1 && iMat != fMatIndex2 ) continue;
205 
206  for(jMat=0;jMat<numOfCouples;jMat++) // transition iMat -> jMat !!!
207  {
208  if( iMat == jMat || ( jMat != fMatIndex1 && jMat != fMatIndex2 ) )
209  {
210  continue;
211  }
212  else
213  {
214  const G4MaterialCutsCouple* iCouple = theCoupleTable->GetMaterialCutsCouple(iMat);
215  const G4MaterialCutsCouple* jCouple = theCoupleTable->GetMaterialCutsCouple(jMat);
216  const G4Material* mat1 = iCouple->GetMaterial();
217  const G4Material* mat2 = jCouple->GetMaterial();
218 
221 
222  // fGammaTkinCut = fGammaCutInKineticEnergy[jMat]; // TR photon in jMat !
223 
224  fGammaTkinCut = 0.0;
225 
226  if(fGammaTkinCut > fTheMinEnergyTR) // setting of min/max TR energies
227  {
229  }
230  else
231  {
233  }
234  if(fGammaTkinCut > fTheMaxEnergyTR)
235  {
236  fMaxEnergyTR = 2.0*fGammaTkinCut; // usually very low TR rate
237  }
238  else
239  {
241  }
242  for(iTkin=0;iTkin<fTotBin;iTkin++) // Lorentz factor loop
243  {
245  energyVector = new G4PhysicsLogVector( fMinEnergyTR,
246  fMaxEnergyTR,
247  fBinTR );
248 
249  fGamma = 1.0 + (fProtonEnergyVector->
250  GetLowEdgeEnergy(iTkin)/proton_mass_c2);
251 
252  fMaxThetaTR = 10000.0/(fGamma*fGamma);
253 
255  {
257  }
258  else
259  {
261  {
263  }
264  }
265  // G4cout<<G4endl<<"fGamma = "<<fGamma<<" fMaxThetaTR = "<<fMaxThetaTR<<G4endl;
267  angleVector = new G4PhysicsLinearVector( 0.0,
268  fMaxThetaTR,
269  fBinTR );
270  G4double energySum = 0.0;
271  G4double angleSum = 0.0;
272 
273  energyVector->PutValue(fBinTR-1,energySum);
274  angleVector->PutValue(fBinTR-1,angleSum);
275 
276  for(iTR=fBinTR-2;iTR>=0;iTR--)
277  {
278  energySum += fCofTR*EnergySum(energyVector->GetLowEdgeEnergy(iTR),
279  energyVector->GetLowEdgeEnergy(iTR+1));
280 
281  angleSum += fCofTR*AngleSum(angleVector->GetLowEdgeEnergy(iTR),
282  angleVector->GetLowEdgeEnergy(iTR+1));
283 
284  energyVector->PutValue(iTR,energySum);
285  angleVector ->PutValue(iTR,angleSum);
286  }
287  // G4cout<<"sumE = "<<energySum<<"; sumA = "<<angleSum<<G4endl;
288 
289  if(jMat < iMat)
290  {
291  iPlace = fTotBin+iTkin; // (iMat*(numOfMat-1)+jMat)*
292  }
293  else // jMat > iMat right part of matrices (jMat-1) !
294  {
295  iPlace = iTkin; // (iMat*(numOfMat-1)+jMat-1)*fTotBin+
296  }
297  fEnergyDistrTable->insertAt(iPlace,energyVector);
298  fAngleDistrTable->insertAt(iPlace,angleVector);
299  } // iTkin
300  } // jMat != iMat
301  } // jMat
302  } // iMat
303  // G4cout<<"G4ForwardXrayTR::BuildXrayTRtables have been called"<<G4endl;
304 }
305 
307 //
308 // This function returns the spectral and angle density of TR quanta
309 // in X-ray energy region generated forward when a relativistic
310 // charged particle crosses interface between two materials.
311 // The high energy small theta approximation is applied.
312 // (matter1 -> matter2)
313 // varAngle =2* (1 - std::cos(Theta)) or approximately = Theta*Theta
314 //
315 
316 G4double
318  G4double varAngle ) const
319 {
320  G4double formationLength1, formationLength2;
321  formationLength1 = 1.0/
322  (1.0/(fGamma*fGamma)
323  + fSigma1/(energy*energy)
324  + varAngle);
325  formationLength2 = 1.0/
326  (1.0/(fGamma*fGamma)
327  + fSigma2/(energy*energy)
328  + varAngle);
329  return (varAngle/energy)*(formationLength1 - formationLength2)
330  *(formationLength1 - formationLength2);
331 
332 }
333 
334 
336 //
337 // Analytical formula for angular density of X-ray TR photons
338 //
339 
341  G4double varAngle ) const
342 {
343  G4double x, x2, /*a, b,*/ c, d, f, a2, b2, a4, b4;
344  G4double cof1, cof2, cof3;
345  x = 1.0/energy;
346  x2 = x*x;
347  c = 1.0/fSigma1;
348  d = 1.0/fSigma2;
349  f = (varAngle + 1.0/(fGamma*fGamma));
350  a2 = c*f;
351  b2 = d*f;
352  a4 = a2*a2;
353  b4 = b2*b2;
354  //a = std::sqrt(a2);
355  // b = std::sqrt(b2);
356  cof1 = c*c*(0.5/(a2*(x2 +a2)) +0.5*std::log(x2/(x2 +a2))/a4);
357  cof3 = d*d*(0.5/(b2*(x2 +b2)) +0.5*std::log(x2/(x2 +b2))/b4);
358  cof2 = -c*d*(std::log(x2/(x2 +b2))/b2 - std::log(x2/(x2 +a2))/a2)/(a2 - b2) ;
359  return -varAngle*(cof1 + cof2 + cof3);
360 }
361 
363 //
364 // Definite integral of X-ray TR spectral-angle density from energy1
365 // to energy2
366 //
367 
369  G4double energy2,
370  G4double varAngle ) const
371 {
372  return AngleDensity(energy2,varAngle)
373  - AngleDensity(energy1,varAngle);
374 }
375 
377 //
378 // Integral angle distribution of X-ray TR photons based on analytical
379 // formula for angle density
380 //
381 
383  G4double varAngle2 ) const
384 {
385  G4int i;
386  G4double h , sumEven = 0.0 , sumOdd = 0.0;
387  h = 0.5*(varAngle2 - varAngle1)/fSympsonNumber;
388  for(i=1;i<fSympsonNumber;i++)
389  {
390  sumEven += EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle1 + 2*i*h );
392  varAngle1 + (2*i - 1)*h );
393  }
395  varAngle1 + (2*fSympsonNumber - 1)*h );
396 
397  return h*(EnergyInterval(fMinEnergyTR,fMaxEnergyTR,varAngle1)
399  + 4.0*sumOdd + 2.0*sumEven )/3.0;
400 }
401 
403 //
404 // Analytical Expression for spectral density of Xray TR photons
405 // x = 2*(1 - std::cos(Theta)) ~ Theta^2
406 //
407 
409  G4double x ) const
410 {
411  G4double a, b;
412  a = 1.0/(fGamma*fGamma)
413  + fSigma1/(energy*energy) ;
414  b = 1.0/(fGamma*fGamma)
415  + fSigma2/(energy*energy) ;
416  return ( (a + b)*std::log((x + b)/(x + a))/(a - b)
417  + a/(x + a) + b/(x + b) )/energy;
418 
419 }
420 
422 //
423 // The spectral density in some angle interval from varAngle1 to
424 // varAngle2
425 //
426 
428  G4double varAngle1,
429  G4double varAngle2 ) const
430 {
431  return SpectralDensity(energy,varAngle2)
432  - SpectralDensity(energy,varAngle1);
433 }
434 
436 //
437 // Integral spectral distribution of X-ray TR photons based on
438 // analytical formula for spectral density
439 //
440 
442  G4double energy2 ) const
443 {
444  G4int i;
445  G4double h , sumEven = 0.0 , sumOdd = 0.0;
446  h = 0.5*(energy2 - energy1)/fSympsonNumber;
447  for(i=1;i<fSympsonNumber;i++)
448  {
449  sumEven += AngleInterval(energy1 + 2*i*h,0.0,fMaxThetaTR);
450  sumOdd += AngleInterval(energy1 + (2*i - 1)*h,0.0,fMaxThetaTR);
451  }
452  sumOdd += AngleInterval(energy1 + (2*fSympsonNumber - 1)*h,
453  0.0,fMaxThetaTR);
454 
455  return h*( AngleInterval(energy1,0.0,fMaxThetaTR)
456  + AngleInterval(energy2,0.0,fMaxThetaTR)
457  + 4.0*sumOdd + 2.0*sumEven )/3.0;
458 }
459 
461 //
462 // PostStepDoIt function for creation of forward X-ray photons in TR process
463 // on boubndary between two materials with really different plasma energies
464 //
465 
467  const G4Step& aStep)
468 {
469  aParticleChange.Initialize(aTrack);
470  // G4cout<<"call G4ForwardXrayTR::PostStepDoIt"<<G4endl;
471  G4int iMat, jMat, iTkin, iPlace, numOfTR, iTR, iTransfer;
472 
473  G4double energyPos, anglePos, energyTR, theta, phi, dirX, dirY, dirZ;
474  G4double W, W1, W2, E1, E2;
475 
476  G4StepPoint* pPreStepPoint = aStep.GetPreStepPoint();
477  G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint();
479 
480  if (pPostStepPoint->GetStepStatus() != fGeomBoundary)
481  {
482  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
483  }
484  if (aTrack.GetStepLength() <= tol)
485  {
486  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
487  }
488  // Come on boundary, so begin to try TR
489 
490  const G4MaterialCutsCouple* iCouple = pPreStepPoint ->GetPhysicalVolume()->
491  GetLogicalVolume()->GetMaterialCutsCouple();
492  const G4MaterialCutsCouple* jCouple = pPostStepPoint ->GetPhysicalVolume()->
493  GetLogicalVolume()->GetMaterialCutsCouple();
494  const G4Material* iMaterial = iCouple->GetMaterial();
495  const G4Material* jMaterial = jCouple->GetMaterial();
496  iMat = iCouple->GetIndex();
497  jMat = jCouple->GetIndex();
498 
499  // The case of equal or approximate (in terms of plasma energy) materials
500  // No TR photons ?!
501 
502  if ( iMat == jMat
503  || ( (fMatIndex1 >= 0 && fMatIndex1 >= 0)
504  && ( iMat != fMatIndex1 && iMat != fMatIndex2 )
505  && ( jMat != fMatIndex1 && jMat != fMatIndex2 ) )
506 
507  || iMaterial->GetState() == jMaterial->GetState()
508 
509  ||(iMaterial->GetState() == kStateSolid && jMaterial->GetState() == kStateLiquid )
510 
511  ||(iMaterial->GetState() == kStateLiquid && jMaterial->GetState() == kStateSolid ) )
512  {
513  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
514  }
515 
516  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
517  G4double charge = aParticle->GetDefinition()->GetPDGCharge();
518 
519  if(charge == 0.0) // Uncharged particle doesn't Generate TR photons
520  {
521  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
522  }
523  // Now we are ready to Generate TR photons
524 
525  G4double chargeSq = charge*charge;
526  G4double kinEnergy = aParticle->GetKineticEnergy();
527  G4double massRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass();
528  G4double TkinScaled = kinEnergy*massRatio;
529  for(iTkin=0;iTkin<fTotBin;iTkin++)
530  {
531  if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) // <= ?
532  {
533  break;
534  }
535  }
536  if(jMat < iMat)
537  {
538  iPlace = fTotBin + iTkin - 1; // (iMat*(numOfMat - 1) + jMat)*
539  }
540  else
541  {
542  iPlace = iTkin - 1; // (iMat*(numOfMat - 1) + jMat - 1)*fTotBin +
543  }
544  // G4PhysicsVector* energyVector1 = (*fEnergyDistrTable)(iPlace) ;
545  // G4PhysicsVector* energyVector2 = (*fEnergyDistrTable)(iPlace + 1);
546 
547  // G4PhysicsVector* angleVector1 = (*fAngleDistrTable)(iPlace) ;
548  // G4PhysicsVector* angleVector2 = (*fAngleDistrTable)(iPlace + 1) ;
549 
550  G4ParticleMomentum particleDir = aParticle->GetMomentumDirection();
551 
552  if(iTkin == fTotBin) // TR plato, try from left
553  {
554  // G4cout<<iTkin<<" mean TR number = "<<( (*(*fEnergyDistrTable)(iPlace))(0) +
555  // (*(*fAngleDistrTable)(iPlace))(0) )
556  // *chargeSq*0.5<<G4endl;
557 
558  numOfTR = G4Poisson( ( (*(*fEnergyDistrTable)(iPlace))(0) +
559  (*(*fAngleDistrTable)(iPlace))(0) )
560  *chargeSq*0.5 );
561  if(numOfTR == 0)
562  {
563  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
564  }
565  else
566  {
567  // G4cout<<"Number of X-ray TR photons = "<<numOfTR<<G4endl;
568 
570 
571  for(iTR=0;iTR<numOfTR;iTR++)
572  {
573  energyPos = (*(*fEnergyDistrTable)(iPlace))(0)*G4UniformRand();
574  for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
575  {
576  if(energyPos >= (*(*fEnergyDistrTable)(iPlace))(iTransfer)) break;
577  }
578  energyTR = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer);
579 
580  // G4cout<<"energyTR = "<<energyTR/keV<<"keV"<<G4endl;
581 
582  kinEnergy -= energyTR;
583  aParticleChange.ProposeEnergy(kinEnergy);
584 
585  anglePos = (*(*fAngleDistrTable)(iPlace))(0)*G4UniformRand();
586  for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
587  {
588  if(anglePos > (*(*fAngleDistrTable)(iPlace))(iTransfer)) break;
589  }
590  theta = std::sqrt((*fAngleDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1));
591 
592  // G4cout<<iTransfer<<" : theta = "<<theta<<G4endl;
593 
594  phi = twopi*G4UniformRand();
595  dirX = std::sin(theta)*std::cos(phi) ;
596  dirY = std::sin(theta)*std::sin(phi) ;
597  dirZ = std::cos(theta) ;
598  G4ThreeVector directionTR(dirX,dirY,dirZ);
599  directionTR.rotateUz(particleDir);
601  directionTR,
602  energyTR );
603  aParticleChange.AddSecondary(aPhotonTR);
604  }
605  }
606  }
607  else
608  {
609  if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
610  {
611  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
612  }
613  else // general case: Tkin between two vectors of the material
614  {
615  E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1);
617  W = 1.0/(E2 - E1);
618  W1 = (E2 - TkinScaled)*W;
619  W2 = (TkinScaled - E1)*W;
620 
621  // G4cout<<iTkin<<" mean TR number = "<<(((*(*fEnergyDistrTable)(iPlace))(0)+
622  // (*(*fAngleDistrTable)(iPlace))(0))*W1 +
623  // ((*(*fEnergyDistrTable)(iPlace + 1))(0)+
624  // (*(*fAngleDistrTable)(iPlace + 1))(0))*W2)
625  // *chargeSq*0.5<<G4endl;
626 
627  numOfTR = G4Poisson((((*(*fEnergyDistrTable)(iPlace))(0)+
628  (*(*fAngleDistrTable)(iPlace))(0))*W1 +
629  ((*(*fEnergyDistrTable)(iPlace + 1))(0)+
630  (*(*fAngleDistrTable)(iPlace + 1))(0))*W2)
631  *chargeSq*0.5 );
632  if(numOfTR == 0)
633  {
634  return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep);
635  }
636  else
637  {
638  // G4cout<<"Number of X-ray TR photons = "<<numOfTR<<G4endl;
639 
641  for(iTR=0;iTR<numOfTR;iTR++)
642  {
643  energyPos = ((*(*fEnergyDistrTable)(iPlace))(0)*W1+
644  (*(*fEnergyDistrTable)(iPlace + 1))(0)*W2)*G4UniformRand();
645  for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
646  {
647  if(energyPos >= ((*(*fEnergyDistrTable)(iPlace))(iTransfer)*W1+
648  (*(*fEnergyDistrTable)(iPlace + 1))(iTransfer)*W2)) break;
649  }
650  energyTR = ((*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer))*W1+
651  ((*fEnergyDistrTable)(iPlace + 1)->GetLowEdgeEnergy(iTransfer))*W2;
652 
653  // G4cout<<"energyTR = "<<energyTR/keV<<"keV"<<G4endl;
654 
655  kinEnergy -= energyTR;
656  aParticleChange.ProposeEnergy(kinEnergy);
657 
658  anglePos = ((*(*fAngleDistrTable)(iPlace))(0)*W1+
659  (*(*fAngleDistrTable)(iPlace + 1))(0)*W2)*G4UniformRand();
660  for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
661  {
662  if(anglePos > ((*(*fAngleDistrTable)(iPlace))(iTransfer)*W1+
663  (*(*fAngleDistrTable)(iPlace + 1))(iTransfer)*W2)) break;
664  }
665  theta = std::sqrt(((*fAngleDistrTable)(iPlace)->
666  GetLowEdgeEnergy(iTransfer-1))*W1+
667  ((*fAngleDistrTable)(iPlace + 1)->
668  GetLowEdgeEnergy(iTransfer-1))*W2);
669 
670  // G4cout<<iTransfer<<" : theta = "<<theta<<G4endl;
671 
672  phi = twopi*G4UniformRand();
673  dirX = std::sin(theta)*std::cos(phi) ;
674  dirY = std::sin(theta)*std::sin(phi) ;
675  dirZ = std::cos(theta) ;
676  G4ThreeVector directionTR(dirX,dirY,dirZ);
677  directionTR.rotateUz(particleDir);
679  directionTR,
680  energyTR );
681  aParticleChange.AddSecondary(aPhotonTR);
682  }
683  }
684  }
685  }
686  return &aParticleChange;
687 }
688 
690 //
691 // Test function for checking of PostStepDoIt random preparation of TR photon
692 // energy
693 //
694 
695 G4double
697 {
698  G4int iPlace, numOfTR, iTR, iTransfer;
699  G4double energyTR = 0.0; // return this value for no TR photons
700  G4double energyPos ;
701  G4double W1, W2;
702 
703  const G4ProductionCutsTable* theCoupleTable=
705  G4int numOfCouples = theCoupleTable->GetTableSize();
706 
707  // The case of equal or approximate (in terms of plasma energy) materials
708  // No TR photons ?!
709 
710  const G4MaterialCutsCouple* iCouple = theCoupleTable->GetMaterialCutsCouple(iMat);
711  const G4MaterialCutsCouple* jCouple = theCoupleTable->GetMaterialCutsCouple(jMat);
712  const G4Material* iMaterial = iCouple->GetMaterial();
713  const G4Material* jMaterial = jCouple->GetMaterial();
714 
715  if ( iMat == jMat
716 
717  || iMaterial->GetState() == jMaterial->GetState()
718 
719  ||(iMaterial->GetState() == kStateSolid && jMaterial->GetState() == kStateLiquid )
720 
721  ||(iMaterial->GetState() == kStateLiquid && jMaterial->GetState() == kStateSolid ) )
722 
723  {
724  return energyTR;
725  }
726 
727  if(jMat < iMat)
728  {
729  iPlace = (iMat*(numOfCouples - 1) + jMat)*fTotBin + iTkin - 1;
730  }
731  else
732  {
733  iPlace = (iMat*(numOfCouples - 1) + jMat - 1)*fTotBin + iTkin - 1;
734  }
735  G4PhysicsVector* energyVector1 = (*fEnergyDistrTable)(iPlace) ;
736  G4PhysicsVector* energyVector2 = (*fEnergyDistrTable)(iPlace + 1);
737 
738  if(iTkin == fTotBin) // TR plato, try from left
739  {
740  numOfTR = G4Poisson( (*energyVector1)(0) );
741  if(numOfTR == 0)
742  {
743  return energyTR;
744  }
745  else
746  {
747  for(iTR=0;iTR<numOfTR;iTR++)
748  {
749  energyPos = (*energyVector1)(0)*G4UniformRand();
750  for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
751  {
752  if(energyPos >= (*energyVector1)(iTransfer)) break;
753  }
754  energyTR += energyVector1->GetLowEdgeEnergy(iTransfer);
755  }
756  }
757  }
758  else
759  {
760  if(iTkin == 0) // Tkin is too small, neglect of TR photon generation
761  {
762  return energyTR;
763  }
764  else // general case: Tkin between two vectors of the material
765  { // use trivial mean half/half
766  W1 = 0.5;
767  W2 = 0.5;
768  numOfTR = G4Poisson( (*energyVector1)(0)*W1 +
769  (*energyVector2)(0)*W2 );
770  if(numOfTR == 0)
771  {
772  return energyTR;
773  }
774  else
775  {
776  G4cout<<"It is still OK in GetEnergyTR(int,int,int)"<<G4endl;
777  for(iTR=0;iTR<numOfTR;iTR++)
778  {
779  energyPos = ((*energyVector1)(0)*W1+
780  (*energyVector2)(0)*W2)*G4UniformRand();
781  for(iTransfer=0;iTransfer<fBinTR-1;iTransfer++)
782  {
783  if(energyPos >= ((*energyVector1)(iTransfer)*W1+
784  (*energyVector2)(iTransfer)*W2)) break;
785  }
786  energyTR += (energyVector1->GetLowEdgeEnergy(iTransfer))*W1+
787  (energyVector2->GetLowEdgeEnergy(iTransfer))*W2;
788 
789  }
790  }
791  }
792  }
793 
794  return energyTR;
795 }
796 
798 //
799 // Test function for checking of PostStepDoIt random preparation of TR photon
800 // theta angle relative to particle direction
801 //
802 
803 
804 G4double
806 {
807  G4double theta = 0.0;
808 
809  return theta;
810 }
811 
812 
813 
814 // end of G4ForwardXrayTR implementation file
815 //
static G4double fMaxProtonTkin
G4double AngleDensity(G4double energy, G4double varAngle) const
G4double condition(const G4ErrorSymMatrix &m)
G4double AngleInterval(G4double energy, G4double varAngle1, G4double varAngle2) const
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
static G4int fBinTR
Double_t x2[nxs]
Definition: Style.C:19
tuple a
Definition: test.py:11
G4double GetKineticEnergy() const
Float_t d
Definition: plot.C:237
static G4double fTheMaxEnergyTR
G4double GetEnergyTR(G4int iMat, G4int jMat, G4int iTkin) const
const G4DynamicParticle * GetDynamicParticle() const
static G4double fMinProtonTkin
G4ForwardXrayTR(const G4String &matName1, const G4String &matName2, const G4String &processName="XrayTR")
G4double SpectralDensity(G4double energy, G4double x) const
G4StepStatus GetStepStatus() const
const G4String & GetName() const
Definition: G4Material.hh:176
G4double GetSurfaceTolerance() const
G4double AngleSum(G4double varAngle1, G4double varAngle2) const
G4ParticleDefinition * GetDefinition() const
G4double GetThetaTR(G4int iMat, G4int jMat, G4int iTkin) const
tuple x
Definition: test.py:50
G4double GetLowEdgeEnergy(size_t binNumber) const
int G4int
Definition: G4Types.hh:78
G4PhysicsTable * fAngleDistrTable
TFile f
Definition: plotHisto.C:6
static G4double fPlasmaCof
static G4double fTheMinAngle
G4StepPoint * GetPreStepPoint() const
tuple b
Definition: test.py:12
static G4int fSympsonNumber
#define G4UniformRand()
Definition: Randomize.hh:87
G4GLOB_DLL std::ostream G4cout
double energy
Definition: plottest35.C:25
G4VPhysicalVolume * GetPhysicalVolume() const
G4double GetElectronDensity() const
Definition: G4Material.hh:215
G4VParticleChange * PostStepDoIt(const G4Track &aTrack, const G4Step &aStep)
bool G4bool
Definition: G4Types.hh:79
const G4ThreeVector & GetMomentumDirection() const
Hep3Vector & rotateUz(const Hep3Vector &)
Definition: ThreeVector.cc:72
void PutValue(size_t index, G4double theValue)
virtual ~G4ForwardXrayTR()
float proton_mass_c2
Definition: hepunit.py:275
static G4int fTotBin
float electron_mass_c2
Definition: hepunit.py:274
Definition: G4Step.hh:76
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetMeanFreePath(const G4Track &, G4double, G4ForceCondition *condition)
G4PhysicsTable * fEnergyDistrTable
virtual void Initialize(const G4Track &)
static G4ProductionCutsTable * GetProductionCutsTable()
G4double GetPDGMass() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
G4ParticleDefinition * fPtrGamma
G4double SpectralAngleTRdensity(G4double energy, G4double varAngle) const
void SetNumberOfSecondaries(G4int totSecondaries)
G4double EnergySum(G4double energy1, G4double energy2) const
G4StepPoint * GetPostStepPoint() const
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
void AddSecondary(G4Track *aSecondary)
void insertAt(size_t, G4PhysicsVector *)
static G4double fCofTR
static G4double fTheMinEnergyTR
#define G4endl
Definition: G4ios.hh:61
G4State GetState() const
Definition: G4Material.hh:179
static G4double fTheMaxAngle
double G4double
Definition: G4Types.hh:76
tuple c
Definition: test.py:13
G4ForceCondition
G4double GetPDGCharge() const
#define DBL_MAX
Definition: templates.hh:83
G4double EnergyInterval(G4double energy1, G4double energy2, G4double varAngle) const
G4PhysicsLogVector * fProtonEnergyVector
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
const std::vector< G4double > * fGammaCutInKineticEnergy
G4double GetStepLength() const
static G4GeometryTolerance * GetInstance()
const G4Material * GetMaterial() const