Geant4_10
G4XTRRegularRadModel.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 
28 #include <complex>
29 
30 #include "G4XTRRegularRadModel.hh"
31 #include "G4PhysicalConstants.hh"
32 #include "Randomize.hh"
33 
34 #include "G4Gamma.hh"
35 using namespace std;
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<<" XTR Regular discrete radiator model is called"<<G4endl ;
48 
49  fExitFlux = true;
50 
51  // Build energy and angular integral spectra of X-ray TR photons from
52  // a radiator
53 
54  // BuildTable() ;
55 }
56 
58 
60 {
61  ;
62 }
63 
65 //
66 //
67 
69 {
70  G4double result, sum = 0., tmp, cof1, cof2, cofMin, cofPHC, theta2, theta2k;
71  G4double aMa, bMb ,sigma, dump;
72  G4int k, kMax, kMin;
73 
74  aMa = fPlateThick*GetPlateLinearPhotoAbs(energy);
75  bMb = fGasThick*GetGasLinearPhotoAbs(energy);
76  sigma = 0.5*(aMa + bMb);
77  dump = std::exp(-fPlateNumber*sigma);
78  if(verboseLevel > 2) G4cout<<" dump = "<<dump<<G4endl;
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 = 2*( 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 *= dump*( -1 + dump + 2*fPlateNumber );
141  /*
142  fEnergy = energy;
143  // G4Integrator<G4VXTRenergyLoss,G4double(G4VXTRenergyLoss::*)(G4double)> integral;
144  G4Integrator<G4TransparentRegXTRadiator,G4double(G4VXTRenergyLoss::*)(G4double)> integral;
145 
146  tmp = integral.Legendre96(this,&G4VXTRenergyLoss::SpectralAngleXTRdEdx,
147  0.0,0.3*fMaxThetaTR) +
148  integral.Legendre96(this,&G4VXTRenergyLoss::SpectralAngleXTRdEdx,
149  0.3*fMaxThetaTR,0.6*fMaxThetaTR) +
150  integral.Legendre96(this,&G4VXTRenergyLoss::SpectralAngleXTRdEdx,
151  0.6*fMaxThetaTR,fMaxThetaTR) ;
152  result += tmp;
153  */
154  return result;
155 }
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  G4double result, Qa, Qb, Q, aZa, bZb, aMa, bMb, I2 ;
172 
173  aZa = fPlateThick/GetPlateFormationZone(energy,gamma,varAngle) ;
174  bZb = fGasThick/GetGasFormationZone(energy,gamma,varAngle) ;
175 
176  aMa = fPlateThick*GetPlateLinearPhotoAbs(energy) ;
177  bMb = fGasThick*GetGasLinearPhotoAbs(energy) ;
178 
179  Qa = std::exp(-aMa) ;
180  Qb = std::exp(-bMb) ;
181  Q = Qa*Qb ;
182 
183  // G4complex Ca(1.0+0.5*fPlateThick*Ma,fPlateThick/Za) ;
184  // G4complex Cb(1.0+0.5*fGasThick*Mb,fGasThick/Zb) ;
185 
186  G4complex Ha( std::exp(-0.5*aMa)*std::cos(aZa),
187  -std::exp(-0.5*aMa)*std::sin(aZa) ) ;
188 
189  G4complex Hb( std::exp(-0.5*bMb)*std::cos(bZb),
190  -std::exp(-0.5*bMb)*std::sin(bZb) ) ;
191 
192  G4complex H = Ha*Hb ;
193 
194  G4complex Hs = std::conj(H) ;
195 
196  // G4complex F1 = ( 0.5*(1+Qa)*(1+H) - Ha - Qa*Hb )/(1-H) ;
197 
198  G4complex F2 = (1.0-Ha)*(Qa-Ha)*Hb*(1.0-Hs)*(Q-Hs) ;
199 
200  F2 *= std::pow(Q,G4double(fPlateNumber)) - std::pow(H,fPlateNumber) ;
201 
202  result = ( 1 - std::pow(Q,G4double(fPlateNumber)) )/( 1 - Q ) ;
203 
204  result *= (1 - Qa)*(1 + Qa - 2*std::sqrt(Qa)*std::cos(aZa)) ;
205 
206  result /= (1 - std::sqrt(Q))*(1 - std::sqrt(Q)) +
207  4*std::sqrt(Q)*std::sin(0.5*(aZa+bZb))*std::sin(0.5*(aZa+bZb)) ;
208 
209  I2 = 1.; // 2.0*std::real(F2) ;
210 
211  I2 /= (1 - std::sqrt(Q))*(1 - std::sqrt(Q)) +
212  4*std::sqrt(Q)*std::sin(0.5*(aZa+bZb))*std::sin(0.5*(aZa+bZb)) ;
213 
214  I2 /= Q*( (std::sqrt(Q)-std::cos(aZa+bZb))*(std::sqrt(Q)-std::cos(aZa+bZb)) +
215  std::sin(aZa+bZb)*std::sin(aZa+bZb) ) ;
216 
217  G4complex stack = 2.*I2*F2;
218  stack += result;
219  stack *= OneInterfaceXTRdEdx(energy,gamma,varAngle);
220 
221  // result += I2 ;
222  result = std::real(stack);
223 
224  return result ;
225 }
226 
227 
228 //
229 //
231 
232 
233 
234 
235 
236 
237 
238 
G4double GetStackFactor(G4double energy, G4double gamma, G4double varAngle)
Float_t tmp
Definition: plot.C:37
G4double GetGasFormationZone(G4double, G4double, G4double)
G4int verboseLevel
Definition: G4VProcess.hh:368
tuple a
Definition: test.py:11
G4double GetPlateLinearPhotoAbs(G4double)
G4double GetGasLinearPhotoAbs(G4double)
G4double G4NeutronHPJENDLHEData::G4double result
int G4int
Definition: G4Types.hh:78
G4double GetPlateFormationZone(G4double, G4double, G4double)
G4XTRRegularRadModel(G4LogicalVolume *anEnvelope, G4Material *, G4Material *, G4double, G4double, G4int, const G4String &processName="XTRegularModel")
tuple b
Definition: test.py:12
Char_t n[5]
std::complex< G4double > G4complex
Definition: G4Types.hh:81
G4GLOB_DLL std::ostream G4cout
double energy
Definition: plottest35.C:25
G4double SpectralXTRdEdx(G4double energy)
G4complex OneInterfaceXTRdEdx(G4double energy, G4double gamma, G4double varAngle)
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76