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