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