Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/electromagnetic/TestEm10/src/XTRTransparentRegRadModel.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

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