ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4XTRTransparentRegRadModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4XTRTransparentRegRadModel.cc
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 
28 #include <complex>
29 
31 #include "G4PhysicalConstants.hh"
32 #include "Randomize.hh"
33 #include "G4Integrator.hh"
34 #include "G4Gamma.hh"
35 
37 //
38 // Constructor, destructor
39 
41  G4Material* foilMat,G4Material* gasMat,
43  const G4String& processName) :
44  G4VXTRenergyLoss(anEnvelope,foilMat,gasMat,a,b,n,processName)
45 {
46  G4cout<<"Regular transparent X-ray TR radiator EM process is called"<<G4endl;
47 
48  // Build energy and angular integral spectra of X-ray TR photons from
49  // a radiator
50  fExitFlux = true;
51  fAlphaPlate = 10000;
52  fAlphaGas = 1000;
53 
54  // BuildTable();
55 }
56 
58 
60 {
61  ;
62 }
63 
65 //
66 //
67 
69 {
70  G4double result, sum = 0., tmp, cof1, cof2, cofMin, cofPHC,aMa, bMb, sigma;
71  G4int k, kMax, kMin;
72 
73  aMa = GetPlateLinearPhotoAbs(energy);
74  bMb = GetGasLinearPhotoAbs(energy);
75 
76  if(fCompton)
77  {
78  aMa += GetPlateCompton(energy);
79  bMb += GetGasCompton(energy);
80  }
81  aMa *= fPlateThick;
82  bMb *= fGasThick;
83 
84  sigma = aMa + bMb;
85 
86  cofPHC = 4.*pi*hbarc;
87  tmp = (fSigma1 - fSigma2)/cofPHC/energy;
88  cof1 = fPlateThick*tmp;
89  cof2 = fGasThick*tmp;
90 
91  cofMin = energy*(fPlateThick + fGasThick)/fGamma/fGamma;
92  cofMin += (fPlateThick*fSigma1 + fGasThick*fSigma2)/energy;
93  cofMin /= cofPHC;
94 
95  // if (fGamma < 1200) kMin = G4int(cofMin); // 1200 ?
96  // else kMin = 1;
97 
98 
99  kMin = G4int(cofMin);
100  if (cofMin > kMin) kMin++;
101 
102  // tmp = (fPlateThick + fGasThick)*energy*fMaxThetaTR;
103  // tmp /= cofPHC;
104  // kMax = G4int(tmp);
105  // if(kMax < 0) kMax = 0;
106  // kMax += kMin;
107 
108 
109  kMax = kMin + 19; // 5; // 9; // kMin + G4int(tmp);
110 
111  // tmp /= fGamma;
112  // if( G4int(tmp) < kMin ) kMin = G4int(tmp);
113  // G4cout<<"kMin = "<<kMin<<"; kMax = "<<kMax<<G4endl;
114 
115  for( k = kMin; k <= kMax; k++ )
116  {
117  tmp = pi*fPlateThick*(k + cof2)/(fPlateThick + fGasThick);
118  result = (k - cof1)*(k - cof1)*(k + cof2)*(k + cof2);
119 
120  if( k == kMin && kMin == G4int(cofMin) )
121  {
122  sum += 0.5*std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
123  }
124  else
125  {
126  sum += std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
127  }
128  // G4cout<<"k = "<<k<<"; sum = "<<sum<<G4endl;
129  }
130  result = 4.*( cof1 + cof2 )*( cof1 + cof2 )*sum/energy;
131  result *= ( 1. - std::exp(-fPlateNumber*sigma) )/( 1. - std::exp(-sigma) );
132  return result;
133 }
134 
135 
137 //
138 // Approximation for radiator interference factor for the case of
139 // fully Regular radiator. The plate and gas gap thicknesses are fixed .
140 // The mean values of the plate and gas gap thicknesses
141 // are supposed to be about XTR formation zones but much less than
142 // mean absorption length of XTR photons in coresponding material.
143 
144 G4double
146  G4double gamma, G4double varAngle )
147 {
148  /*
149  G4double result, Za, Zb, Ma, Mb, sigma;
150 
151  Za = GetPlateFormationZone(energy,gamma,varAngle);
152  Zb = GetGasFormationZone(energy,gamma,varAngle);
153  Ma = GetPlateLinearPhotoAbs(energy);
154  Mb = GetGasLinearPhotoAbs(energy);
155  sigma = Ma*fPlateThick + Mb*fGasThick;
156 
157  G4complex Ca(1.0+0.5*fPlateThick*Ma/fAlphaPlate,fPlateThick/Za/fAlphaPlate);
158  G4complex Cb(1.0+0.5*fGasThick*Mb/fAlphaGas,fGasThick/Zb/fAlphaGas);
159 
160  G4complex Ha = std::pow(Ca,-fAlphaPlate);
161  G4complex Hb = std::pow(Cb,-fAlphaGas);
162  G4complex H = Ha*Hb;
163  G4complex F1 = (1.0 - Ha)*(1.0 - Hb )/(1.0 - H)
164  * G4double(fPlateNumber) ;
165  G4complex F2 = (1.0-Ha)*(1.0-Ha)*Hb/(1.0-H)/(1.0-H)
166  * (1.0 - std::exp(-0.5*fPlateNumber*sigma)) ;
167  // *(1.0 - std::pow(H,fPlateNumber)) ;
168  G4complex R = (F1 + F2)*OneInterfaceXTRdEdx(energy,gamma,varAngle);
169  // G4complex R = F2*OneInterfaceXTRdEdx(energy,gamma,varAngle);
170  result = 2.0*std::real(R);
171  return result;
172  */
173  // numerically unstable result
174 
175  G4double result, Qa, Qb, Q, aZa, bZb, aMa, bMb, D, sigma;
176 
177  aZa = fPlateThick/GetPlateFormationZone(energy,gamma,varAngle);
178  bZb = fGasThick/GetGasFormationZone(energy,gamma,varAngle);
179  aMa = fPlateThick*GetPlateLinearPhotoAbs(energy);
180  bMb = fGasThick*GetGasLinearPhotoAbs(energy);
181  sigma = aMa*fPlateThick + bMb*fGasThick;
182  Qa = std::exp(-0.5*aMa);
183  Qb = std::exp(-0.5*bMb);
184  Q = Qa*Qb;
185 
186  G4complex Ha( Qa*std::cos(aZa), -Qa*std::sin(aZa) );
187  G4complex Hb( Qb*std::cos(bZb), -Qb*std::sin(bZb) );
188  G4complex H = Ha*Hb;
189  G4complex Hs = conj(H);
190  D = 1.0 /( (1. - Q)*(1. - Q) +
191  4.*Q*std::sin(0.5*(aZa + bZb))*std::sin(0.5*(aZa + bZb)) );
192  G4complex F1 = (1.0 - Ha)*(1.0 - Hb)*(1.0 - Hs)
194  G4complex F2 = (1.0 - Ha)*(1.0 - Ha)*Hb*(1.0 - Hs)*(1.0 - Hs)
195  // * (1.0 - std::pow(H,fPlateNumber)) * D*D;
196  * (1.0 - std::exp(-0.5*fPlateNumber*sigma)) * D*D;
197  G4complex R = (F1 + F2)*OneInterfaceXTRdEdx(energy,gamma,varAngle);
198  result = 2.0*std::real(R);
199  return result;
200 
201 }
202 
203 
204 //
205 //
207 
208 
209 
210 
211 
212 
213 
214