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