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 ]

Diff markup

Differences between /examples/extended/electromagnetic/TestEm10/src/XTRTransparentRegRadModel.cc (Version 11.3.0) and /examples/extended/electromagnetic/TestEm10/src/XTRTransparentRegRadModel.cc (Version 10.4.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
                                                   >>  26 // $Id: XTRTransparentRegRadModel.cc 94932 2015-12-18 09:21:29Z gcosmo $
 26 //                                                 27 //
 27 /// \file electromagnetic/TestEm10/src/XTRTran     28 /// \file electromagnetic/TestEm10/src/XTRTransparentRegRadModel.cc
 28 /// \brief Implementation of the XTRTransparen     29 /// \brief Implementation of the XTRTransparentRegRadModel class
 29 //                                                 30 //
 30 //                                                 31 //
 31                                                    32 
 32 #include "XTRTransparentRegRadModel.hh"        <<  33 #include <complex>
 33                                                    34 
 34 #include "G4Gamma.hh"                          <<  35 #include "XTRTransparentRegRadModel.hh"
                                                   >>  36 #include "Randomize.hh"
 35 #include "G4Integrator.hh"                         37 #include "G4Integrator.hh"
                                                   >>  38 #include "G4Gamma.hh"
 36 #include "G4PhysicalConstants.hh"                  39 #include "G4PhysicalConstants.hh"
 37 #include "Randomize.hh"                        << 
 38                                                << 
 39 #include <complex>                             << 
 40                                                    40 
 41 using namespace std;                               41 using namespace std;
 42                                                    42 
 43 //....oooOO0OOooo........oooOO0OOooo........oo << 
 44                                                << 
 45 //////////////////////////////////////////////     43 ////////////////////////////////////////////////////////////////////////////
 46 //                                                 44 //
 47 // Constructor, destructor                         45 // Constructor, destructor
 48                                                    46 
 49 XTRTransparentRegRadModel::XTRTransparentRegRa <<  47 XTRTransparentRegRadModel::XTRTransparentRegRadModel(G4LogicalVolume *anEnvelope,
 50                                                <<  48                                          G4Material* foilMat,G4Material* gasMat, 
 51                                                <<  49                                          G4double a, G4double b, G4int n,
 52                                                <<  50                                          const G4String& processName) :
 53   : G4VXTRenergyLoss(anEnvelope, foilMat, gasM <<  51   G4VXTRenergyLoss(anEnvelope,foilMat,gasMat,a,b,n,processName)
 54 {                                                  52 {
 55   G4cout << "Regular transparent X-ray TR  rad <<  53   G4cout<<"Regular transparent X-ray TR  radiator EM process is called"<<G4endl;
 56                                                    54 
 57   // Build energy and angular integral spectra     55   // Build energy and angular integral spectra of X-ray TR photons from
 58   // a radiator                                    56   // a radiator
 59   fExitFlux = true;                            <<  57   fExitFlux   = true;
 60   fAlphaPlate = 10000;                             58   fAlphaPlate = 10000;
 61   fAlphaGas = 1000;                            <<  59   fAlphaGas   = 1000;
 62                                                    60 
 63   //  BuildTable();                                61   //  BuildTable();
 64 }                                                  62 }
 65                                                    63 
 66 //....oooOO0OOooo........oooOO0OOooo........oo <<  64 ///////////////////////////////////////////////////////////////////////////
 67                                                    65 
 68 XTRTransparentRegRadModel::~XTRTransparentRegR     66 XTRTransparentRegRadModel::~XTRTransparentRegRadModel()
 69 {                                                  67 {
 70   ;                                                68   ;
 71 }                                                  69 }
 72                                                    70 
 73 //....oooOO0OOooo........oooOO0OOooo........oo <<  71 ///////////////////////////////////////////////////////////////////////////
                                                   >>  72 //
                                                   >>  73 //
 74                                                    74 
 75 G4double XTRTransparentRegRadModel::SpectralXT     75 G4double XTRTransparentRegRadModel::SpectralXTRdEdx(G4double energy)
 76 {                                                  76 {
 77   G4double result, sum = 0., tmp, cof1, cof2,  <<  77   G4double result, sum = 0., tmp, cof1, cof2, cofMin, cofPHC,aMa, bMb, sigma;
 78   G4int k, kMax, kMin;                             78   G4int k, kMax, kMin;
 79                                                    79 
 80   aMa = GetPlateLinearPhotoAbs(energy);            80   aMa = GetPlateLinearPhotoAbs(energy);
 81   bMb = GetGasLinearPhotoAbs(energy);              81   bMb = GetGasLinearPhotoAbs(energy);
 82                                                    82 
 83   // if(fCompton)                                  83   // if(fCompton)
 84   {                                                84   {
 85     aMa += GetPlateCompton(energy);                85     aMa += GetPlateCompton(energy);
 86     bMb += GetGasCompton(energy);                  86     bMb += GetGasCompton(energy);
 87   }                                                87   }
 88   aMa *= fPlateThick;                              88   aMa *= fPlateThick;
 89   bMb *= fGasThick;                                89   bMb *= fGasThick;
 90                                                    90 
 91   sigma = aMa + bMb;                               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;
 92                                                    98 
 93   cofPHC = 4 * pi * hbarc;                     <<  99   cofMin  =  energy*(fPlateThick + fGasThick)/fGamma/fGamma;
 94   cofPHC *= 200. / 197.;                       << 100   cofMin += (fPlateThick*fSigma1 + fGasThick*fSigma2)/energy;
 95   tmp = (fSigma1 - fSigma2) / cofPHC / energy; << 
 96   cof1 = fPlateThick * tmp;                    << 
 97   cof2 = fGasThick * tmp;                      << 
 98                                                << 
 99   cofMin = energy * (fPlateThick + fGasThick)  << 
100   cofMin += (fPlateThick * fSigma1 + fGasThick << 
101   cofMin /= cofPHC;                               101   cofMin /= cofPHC;
102                                                   102 
103   //  if (fGamma < 1200) kMin = G4int(cofMin);    103   //  if (fGamma < 1200) kMin = G4int(cofMin);  // 1200 ?
104   // else               kMin = 1;                 104   // else               kMin = 1;
105                                                   105 
                                                   >> 106 
106   kMin = G4int(cofMin);                           107   kMin = G4int(cofMin);
107   if (cofMin > kMin) kMin++;                      108   if (cofMin > kMin) kMin++;
108                                                   109 
109   // tmp  = (fPlateThick + fGasThick)*energy*f    110   // tmp  = (fPlateThick + fGasThick)*energy*fMaxThetaTR;
110   // tmp /= cofPHC;                               111   // tmp /= cofPHC;
111   // kMax = G4int(tmp);                           112   // kMax = G4int(tmp);
112   // if(kMax < 0) kMax = 0;                       113   // if(kMax < 0) kMax = 0;
113   // kMax += kMin;                                114   // kMax += kMin;
                                                   >> 115   
114                                                   116 
115   kMax = kMin + 9;  // 5; // 9; //   kMin + G4 << 117   kMax = kMin + 9; // 5; // 9; //   kMin + G4int(tmp);
116                                                   118 
117   // tmp /= fGamma;                               119   // tmp /= fGamma;
118   // if( G4int(tmp) < kMin ) kMin = G4int(tmp)    120   // if( G4int(tmp) < kMin ) kMin = G4int(tmp);
119   // G4cout<<"kMin = "<<kMin<<";    kMax = "<<    121   // G4cout<<"kMin = "<<kMin<<";    kMax = "<<kMax<<G4endl;
120                                                   122 
121   for (k = kMin; k <= kMax; k++) {             << 123   for( k = kMin; k <= kMax; k++ )
122     tmp = pi * fPlateThick * (k + cof2) / (fPl << 124   {
123     result = (k - cof1) * (k - cof1) * (k + co << 125     tmp    = pi*fPlateThick*(k + cof2)/(fPlateThick + fGasThick);
                                                   >> 126     result = (k - cof1)*(k - cof1)*(k + cof2)*(k + cof2);
124                                                   127 
125     if (k == kMin && kMin == G4int(cofMin)) {  << 128     if( k == kMin && kMin == G4int(cofMin) )
126       sum += 0.5 * sin(tmp) * sin(tmp) * std:: << 129     {
                                                   >> 130       sum   += 0.5*sin(tmp)*sin(tmp)*std::abs(k-cofMin)/result;
127     }                                             131     }
128     else {                                     << 132     else
129       sum += sin(tmp) * sin(tmp) * std::abs(k  << 133     {
                                                   >> 134       sum   += sin(tmp)*sin(tmp)*std::abs(k-cofMin)/result;
130     }                                             135     }
131     //  G4cout<<"k = "<<k<<";    sum = "<<sum< << 136     //  G4cout<<"k = "<<k<<";    sum = "<<sum<<G4endl;    
132   }                                               137   }
133   result = 4. * (cof1 + cof2) * (cof1 + cof2)  << 138   result = 4.*( cof1 + cof2 )*( cof1 + cof2 )*sum/energy;
134   result *= (1. - exp(-fPlateNumber * sigma))  << 139   result *= ( 1. - exp(-fPlateNumber*sigma) )/( 1. - exp(-sigma) );  
135   return result;                                  140   return result;
136 }                                                 141 }
137                                                   142 
138 //....oooOO0OOooo........oooOO0OOooo........oo << 143 
                                                   >> 144 ///////////////////////////////////////////////////////////////////////////
139 //                                                145 //
140 // Approximation for radiator interference fac    146 // Approximation for radiator interference factor for the case of
141 // fully Regular radiator. The plate and gas g    147 // fully Regular radiator. The plate and gas gap thicknesses are fixed .
142 // The mean values of the plate and gas gap th << 148 // The mean values of the plate and gas gap thicknesses 
143 // are supposed to be about XTR formation zone << 149 // are supposed to be about XTR formation zones but much less than 
144 // mean absorption length of XTR photons in co    150 // mean absorption length of XTR photons in coresponding material.
145                                                   151 
146 G4double XTRTransparentRegRadModel::GetStackFa << 152 G4double 
147                                                << 153 XTRTransparentRegRadModel::GetStackFactor( G4double energy, 
                                                   >> 154                                          G4double gamma, G4double varAngle )
148 {                                                 155 {
149   /*                                              156   /*
150   G4double result, Za, Zb, Ma, Mb, sigma;         157   G4double result, Za, Zb, Ma, Mb, sigma;
151                                                << 158   
152   Za = GetPlateFormationZone(energy,gamma,varA    159   Za = GetPlateFormationZone(energy,gamma,varAngle);
153   Zb = GetGasFormationZone(energy,gamma,varAng    160   Zb = GetGasFormationZone(energy,gamma,varAngle);
154   Ma = GetPlateLinearPhotoAbs(energy);            161   Ma = GetPlateLinearPhotoAbs(energy);
155   Mb = GetGasLinearPhotoAbs(energy);              162   Mb = GetGasLinearPhotoAbs(energy);
156   sigma = Ma*fPlateThick + Mb*fGasThick;          163   sigma = Ma*fPlateThick + Mb*fGasThick;
157                                                   164 
158   G4complex Ca(1.0+0.5*fPlateThick*Ma/fAlphaPl << 165   G4complex Ca(1.0+0.5*fPlateThick*Ma/fAlphaPlate,fPlateThick/Za/fAlphaPlate); 
159   G4complex Cb(1.0+0.5*fGasThick*Mb/fAlphaGas, << 166   G4complex Cb(1.0+0.5*fGasThick*Mb/fAlphaGas,fGasThick/Zb/fAlphaGas); 
160                                                   167 
161   G4complex Ha = pow(Ca,-fAlphaPlate);         << 168   G4complex Ha = pow(Ca,-fAlphaPlate);  
162   G4complex Hb = pow(Cb,-fAlphaGas);              169   G4complex Hb = pow(Cb,-fAlphaGas);
163   G4complex H  = Ha*Hb;                           170   G4complex H  = Ha*Hb;
164   G4complex F1 =   (1.0 - Ha)*(1.0 - Hb )/(1.0    171   G4complex F1 =   (1.0 - Ha)*(1.0 - Hb )/(1.0 - H)
165                  * G4double(fPlateNumber) ;       172                  * G4double(fPlateNumber) ;
166   G4complex F2 =   (1.0-Ha)*(1.0-Ha)*Hb/(1.0-H    173   G4complex F2 =   (1.0-Ha)*(1.0-Ha)*Hb/(1.0-H)/(1.0-H)
167                  * (1.0 - exp(-0.5*fPlateNumbe    174                  * (1.0 - exp(-0.5*fPlateNumber*sigma)) ;
168   //    *(1.0 - pow(H,fPlateNumber)) ;            175   //    *(1.0 - pow(H,fPlateNumber)) ;
169     G4complex R  = (F1 + F2)*OneInterfaceXTRdE    176     G4complex R  = (F1 + F2)*OneInterfaceXTRdEdx(energy,gamma,varAngle);
170   // G4complex R  = F2*OneInterfaceXTRdEdx(ene    177   // G4complex R  = F2*OneInterfaceXTRdEdx(energy,gamma,varAngle);
171   result       = 2.0*real(R);                  << 178   result       = 2.0*real(R);  
172   return      result;                             179   return      result;
173   */                                              180   */
174   // numerically unstable result               << 181    // numerically unstable result
175                                                   182 
176   G4double result, Qa, Qb, Q, aZa, bZb, aMa, b << 183   G4double result, Qa, Qb, Q, aZa, bZb, aMa, bMb, D, sigma; 
                                                   >> 184  
                                                   >> 185   aZa   = fPlateThick/GetPlateFormationZone(energy,gamma,varAngle);
                                                   >> 186   bZb   = fGasThick/GetGasFormationZone(energy,gamma,varAngle);
                                                   >> 187   aMa   = fPlateThick*GetPlateLinearPhotoAbs(energy);
                                                   >> 188   bMb   = fGasThick*GetGasLinearPhotoAbs(energy);
                                                   >> 189   sigma = aMa*fPlateThick + bMb*fGasThick;
                                                   >> 190   Qa    = exp(-0.5*aMa);
                                                   >> 191   Qb    = exp(-0.5*bMb);
                                                   >> 192   Q     = Qa*Qb;
177                                                   193 
178   aZa = fPlateThick / GetPlateFormationZone(en << 194   G4complex Ha( Qa*cos(aZa), -Qa*sin(aZa)   );  
179   bZb = fGasThick / GetGasFormationZone(energy << 195   G4complex Hb( Qb*cos(bZb), -Qb*sin(bZb)    );
180   aMa = fPlateThick * GetPlateLinearPhotoAbs(e << 196   G4complex H  = Ha*Hb;
181   bMb = fGasThick * GetGasLinearPhotoAbs(energ << 
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);                         197   G4complex Hs = conj(H);
191   D = 1.0 / ((1 - Q) * (1 - Q) + 4 * Q * sin(0 << 198   D            = 1.0 /( (1 - Q)*(1 - Q) + 
192   G4complex F1 = (1.0 - Ha) * (1.0 - Hb) * (1. << 199                   4*Q*sin(0.5*(aZa + bZb))*sin(0.5*(aZa + bZb)) );
193   G4complex F2 = (1.0 - Ha) * (1.0 - Ha) * Hb  << 200   G4complex F1 = (1.0 - Ha)*(1.0 - Hb)*(1.0 - Hs)
194                  * (1.0 - Hs)                  << 201                  * G4double(fPlateNumber)*D;
195                  // * (1.0 - pow(H,fPlateNumbe << 202   G4complex F2 = (1.0 - Ha)*(1.0 - Ha)*Hb*(1.0 - Hs)*(1.0 - Hs)
196                  * (1.0 - exp(-0.5 * fPlateNum << 203                    // * (1.0 - pow(H,fPlateNumber)) * D*D;
197   G4complex R = (F1 + F2) * OneInterfaceXTRdEd << 204                  * (1.0 - exp(-0.5*fPlateNumber*sigma)) * D*D;
198   result = 2.0 * real(R);                      << 205   G4complex R  = (F1 + F2)*OneInterfaceXTRdEdx(energy,gamma,varAngle);
199   return result;                               << 206   result       = 2.0*real(R); 
                                                   >> 207   return      result;
                                                   >> 208   
200 }                                                 209 }
                                                   >> 210 
                                                   >> 211 
                                                   >> 212 //
                                                   >> 213 //
                                                   >> 214 ////////////////////////////////////////////////////////////////////////////
                                                   >> 215 
                                                   >> 216 
                                                   >> 217 
                                                   >> 218 
                                                   >> 219 
                                                   >> 220 
                                                   >> 221 
                                                   >> 222 
201                                                   223