G4XTRTransparentRegRadModel.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 
00028 #include <complex>
00029 
00030 #include "G4XTRTransparentRegRadModel.hh"
00031 #include "G4PhysicalConstants.hh"
00032 #include "Randomize.hh"
00033 #include "G4Integrator.hh"
00034 #include "G4Gamma.hh"
00035 
00037 //
00038 // Constructor, destructor
00039 
00040 G4XTRTransparentRegRadModel::G4XTRTransparentRegRadModel(G4LogicalVolume *anEnvelope,
00041                                          G4Material* foilMat,G4Material* gasMat, 
00042                                          G4double a, G4double b, G4int n,
00043                                          const G4String& processName) :
00044   G4VXTRenergyLoss(anEnvelope,foilMat,gasMat,a,b,n,processName)
00045 {
00046   G4cout<<"Regular transparent X-ray TR  radiator EM process is called"<<G4endl;
00047 
00048   // Build energy and angular integral spectra of X-ray TR photons from
00049   // a radiator
00050   fExitFlux   = true;
00051   fAlphaPlate = 10000;
00052   fAlphaGas   = 1000;
00053 
00054   //  BuildTable();
00055 }
00056 
00058 
00059 G4XTRTransparentRegRadModel::~G4XTRTransparentRegRadModel()
00060 {
00061   ;
00062 }
00063 
00065 //
00066 //
00067 
00068 G4double G4XTRTransparentRegRadModel::SpectralXTRdEdx(G4double energy)
00069 {
00070   G4double result, sum = 0., tmp, cof1, cof2, cofMin, cofPHC,aMa, bMb, sigma;
00071   G4int k, kMax, kMin;
00072 
00073   aMa = GetPlateLinearPhotoAbs(energy);
00074   bMb = GetGasLinearPhotoAbs(energy);
00075 
00076   if(fCompton)
00077   {
00078     aMa += GetPlateCompton(energy);
00079     bMb += GetGasCompton(energy);
00080   }
00081   aMa *= fPlateThick;
00082   bMb *= fGasThick;
00083 
00084   sigma = aMa + bMb;
00085    
00086   cofPHC  = 4*pi*hbarc;
00087   tmp     = (fSigma1 - fSigma2)/cofPHC/energy;  
00088   cof1    = fPlateThick*tmp;
00089   cof2    = fGasThick*tmp;
00090 
00091   cofMin  =  energy*(fPlateThick + fGasThick)/fGamma/fGamma;
00092   cofMin += (fPlateThick*fSigma1 + fGasThick*fSigma2)/energy;
00093   cofMin /= cofPHC;
00094 
00095   //  if (fGamma < 1200) kMin = G4int(cofMin);  // 1200 ?
00096   // else               kMin = 1;
00097 
00098 
00099   kMin = G4int(cofMin);
00100   if (cofMin > kMin) kMin++;
00101 
00102   // tmp  = (fPlateThick + fGasThick)*energy*fMaxThetaTR;
00103   // tmp /= cofPHC;
00104   // kMax = G4int(tmp);
00105   // if(kMax < 0) kMax = 0;
00106   // kMax += kMin;
00107   
00108 
00109   kMax = kMin + 19; // 5; // 9; //   kMin + G4int(tmp);
00110 
00111   // tmp /= fGamma;
00112   // if( G4int(tmp) < kMin ) kMin = G4int(tmp);
00113   // G4cout<<"kMin = "<<kMin<<";    kMax = "<<kMax<<G4endl;
00114 
00115   for( k = kMin; k <= kMax; k++ )
00116   {
00117     tmp    = pi*fPlateThick*(k + cof2)/(fPlateThick + fGasThick);
00118     result = (k - cof1)*(k - cof1)*(k + cof2)*(k + cof2);
00119 
00120     if( k == kMin && kMin == G4int(cofMin) )
00121     {
00122       sum   += 0.5*std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
00123     }
00124     else
00125     {
00126       sum   += std::sin(tmp)*std::sin(tmp)*std::abs(k-cofMin)/result;
00127     }
00128     //  G4cout<<"k = "<<k<<";    sum = "<<sum<<G4endl;    
00129   }
00130   result = 4.*( cof1 + cof2 )*( cof1 + cof2 )*sum/energy;
00131   result *= ( 1. - std::exp(-fPlateNumber*sigma) )/( 1. - std::exp(-sigma) );  
00132   return result;
00133 }
00134 
00135 
00137 //
00138 // Approximation for radiator interference factor for the case of
00139 // fully Regular radiator. The plate and gas gap thicknesses are fixed .
00140 // The mean values of the plate and gas gap thicknesses 
00141 // are supposed to be about XTR formation zones but much less than 
00142 // mean absorption length of XTR photons in coresponding material.
00143 
00144 G4double 
00145 G4XTRTransparentRegRadModel::GetStackFactor( G4double energy, 
00146                                          G4double gamma, G4double varAngle )
00147 {
00148   /*
00149   G4double result, Za, Zb, Ma, Mb, sigma;
00150   
00151   Za = GetPlateFormationZone(energy,gamma,varAngle);
00152   Zb = GetGasFormationZone(energy,gamma,varAngle);
00153   Ma = GetPlateLinearPhotoAbs(energy);
00154   Mb = GetGasLinearPhotoAbs(energy);
00155   sigma = Ma*fPlateThick + Mb*fGasThick;
00156 
00157   G4complex Ca(1.0+0.5*fPlateThick*Ma/fAlphaPlate,fPlateThick/Za/fAlphaPlate); 
00158   G4complex Cb(1.0+0.5*fGasThick*Mb/fAlphaGas,fGasThick/Zb/fAlphaGas); 
00159 
00160   G4complex Ha = std::pow(Ca,-fAlphaPlate);  
00161   G4complex Hb = std::pow(Cb,-fAlphaGas);
00162   G4complex H  = Ha*Hb;
00163   G4complex F1 =   (1.0 - Ha)*(1.0 - Hb )/(1.0 - H)
00164                  * G4double(fPlateNumber) ;
00165   G4complex F2 =   (1.0-Ha)*(1.0-Ha)*Hb/(1.0-H)/(1.0-H)
00166                  * (1.0 - std::exp(-0.5*fPlateNumber*sigma)) ;
00167   //    *(1.0 - std::pow(H,fPlateNumber)) ;
00168     G4complex R  = (F1 + F2)*OneInterfaceXTRdEdx(energy,gamma,varAngle);
00169   // G4complex R  = F2*OneInterfaceXTRdEdx(energy,gamma,varAngle);
00170   result       = 2.0*std::real(R);  
00171   return      result;
00172   */
00173    // numerically unstable result
00174 
00175   G4double result, Qa, Qb, Q, aZa, bZb, aMa, bMb, D, sigma; 
00176  
00177   aZa   = fPlateThick/GetPlateFormationZone(energy,gamma,varAngle);
00178   bZb   = fGasThick/GetGasFormationZone(energy,gamma,varAngle);
00179   aMa   = fPlateThick*GetPlateLinearPhotoAbs(energy);
00180   bMb   = fGasThick*GetGasLinearPhotoAbs(energy);
00181   sigma = aMa*fPlateThick + bMb*fGasThick;
00182   Qa    = std::exp(-0.5*aMa);
00183   Qb    = std::exp(-0.5*bMb);
00184   Q     = Qa*Qb;
00185 
00186   G4complex Ha( Qa*std::cos(aZa), -Qa*std::sin(aZa)   );  
00187   G4complex Hb( Qb*std::cos(bZb), -Qb*std::sin(bZb)    );
00188   G4complex H  = Ha*Hb;
00189   G4complex Hs = conj(H);
00190   D            = 1.0 /( (1 - Q)*(1 - Q) + 
00191                   4*Q*std::sin(0.5*(aZa + bZb))*std::sin(0.5*(aZa + bZb)) );
00192   G4complex F1 = (1.0 - Ha)*(1.0 - Hb)*(1.0 - Hs)
00193                  * G4double(fPlateNumber)*D;
00194   G4complex F2 = (1.0 - Ha)*(1.0 - Ha)*Hb*(1.0 - Hs)*(1.0 - Hs)
00195                    // * (1.0 - std::pow(H,fPlateNumber)) * D*D;
00196                  * (1.0 - std::exp(-0.5*fPlateNumber*sigma)) * D*D;
00197   G4complex R  = (F1 + F2)*OneInterfaceXTRdEdx(energy,gamma,varAngle);
00198   result       = 2.0*std::real(R); 
00199   return      result;
00200   
00201 }
00202 
00203 
00204 //
00205 //
00207 
00208 
00209 
00210 
00211 
00212 
00213 
00214 

Generated on Mon May 27 17:50:28 2013 for Geant4 by  doxygen 1.4.7