TestEm10/src/XTRTransparentRegRadModel.cc

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 /// \file electromagnetic/TestEm10/src/XTRTransparentRegRadModel.cc
0028 /// \brief Implementation of the XTRTransparentRegRadModel class
0029 //
0030 //
0031
0032 #include "XTRTransparentRegRadModel.hh"
0033
0034 #include "G4Gamma.hh"
0035 #include "G4Integrator.hh"
0036 #include "G4PhysicalConstants.hh"
0037 #include "Randomize.hh"
0038
0039 #include <complex>
0040
0041 using namespace std;
0042
0043 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0044
0045 ////////////////////////////////////////////////////////////////////////////
0046 //
0047 // Constructor, destructor
0048
0049 XTRTransparentRegRadModel::XTRTransparentRegRadModel(G4LogicalVolume* anEnvelope,
0050                                                      G4Material* foilMat, G4Material* gasMat,
0051                                                      G4double a, G4double b, G4int n,
0052                                                      const G4String& processName)
0053   : G4VXTRenergyLoss(anEnvelope, foilMat, gasMat, a, b, n, processName)
0054 {
0055   G4cout << "Regular transparent X-ray TR  radiator EM process is called" << G4endl;
0056
0057   // Build energy and angular integral spectra of X-ray TR photons from
0058   // a radiator
0059   fExitFlux = true;
0060   fAlphaPlate = 10000;
0061   fAlphaGas = 1000;
0062
0063   //  BuildTable();
0064 }
0065
0066 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0067
0068 XTRTransparentRegRadModel::~XTRTransparentRegRadModel()
0069 {
0070   ;
0071 }
0072
0073 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0074
0075 G4double XTRTransparentRegRadModel::SpectralXTRdEdx(G4double energy)
0076 {
0077   G4double result, sum = 0., tmp, cof1, cof2, cofMin, cofPHC, aMa, bMb, sigma;
0078   G4int k, kMax, kMin;
0079
0080   aMa = GetPlateLinearPhotoAbs(energy);
0081   bMb = GetGasLinearPhotoAbs(energy);
0082
0083   // if(fCompton)
0084   {
0085     aMa += GetPlateCompton(energy);
0086     bMb += GetGasCompton(energy);
0087   }
0088   aMa *= fPlateThick;
0089   bMb *= fGasThick;
0090
0091   sigma = aMa + bMb;
0092
0093   cofPHC = 4 * pi * hbarc;
0094   cofPHC *= 200. / 197.;
0095   tmp = (fSigma1 - fSigma2) / cofPHC / energy;
0096   cof1 = fPlateThick * tmp;
0097   cof2 = fGasThick * tmp;
0098
0099   cofMin = energy * (fPlateThick + fGasThick) / fGamma / fGamma;
0100   cofMin += (fPlateThick * fSigma1 + fGasThick * fSigma2) / energy;
0101   cofMin /= cofPHC;
0102
0103   //  if (fGamma < 1200) kMin = G4int(cofMin);  // 1200 ?
0104   // else               kMin = 1;
0105
0106   kMin = G4int(cofMin);
0107   if (cofMin > kMin) kMin++;
0108
0109   // tmp  = (fPlateThick + fGasThick)*energy*fMaxThetaTR;
0110   // tmp /= cofPHC;
0111   // kMax = G4int(tmp);
0112   // if(kMax < 0) kMax = 0;
0113   // kMax += kMin;
0114
0115   kMax = kMin + 9;  // 5; // 9; //   kMin + G4int(tmp);
0116
0117   // tmp /= fGamma;
0118   // if( G4int(tmp) < kMin ) kMin = G4int(tmp);
0119   // G4cout<<"kMin = "<<kMin<<";    kMax = "<<kMax<<G4endl;
0120
0121   for (k = kMin; k <= kMax; k++) {
0122     tmp = pi * fPlateThick * (k + cof2) / (fPlateThick + fGasThick);
0123     result = (k - cof1) * (k - cof1) * (k + cof2) * (k + cof2);
0124
0125     if (k == kMin && kMin == G4int(cofMin)) {
0126       sum += 0.5 * sin(tmp) * sin(tmp) * std::abs(k - cofMin) / result;
0127     }
0128     else {
0129       sum += sin(tmp) * sin(tmp) * std::abs(k - cofMin) / result;
0130     }
0131     //  G4cout<<"k = "<<k<<";    sum = "<<sum<<G4endl;
0132   }
0133   result = 4. * (cof1 + cof2) * (cof1 + cof2) * sum / energy;
0134   result *= (1. - exp(-fPlateNumber * sigma)) / (1. - exp(-sigma));
0135   return result;
0136 }
0137
0138 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0139 //
0140 // Approximation for radiator interference factor for the case of
0141 // fully Regular radiator. The plate and gas gap thicknesses are fixed .
0142 // The mean values of the plate and gas gap thicknesses
0143 // are supposed to be about XTR formation zones but much less than
0144 // mean absorption length of XTR photons in coresponding material.
0145
0146 G4double XTRTransparentRegRadModel::GetStackFactor(G4double energy, G4double gamma,
0147                                                    G4double varAngle)
0148 {
0149   /*
0150   G4double result, Za, Zb, Ma, Mb, sigma;
0151
0152   Za = GetPlateFormationZone(energy,gamma,varAngle);
0153   Zb = GetGasFormationZone(energy,gamma,varAngle);
0154   Ma = GetPlateLinearPhotoAbs(energy);
0155   Mb = GetGasLinearPhotoAbs(energy);
0156   sigma = Ma*fPlateThick + Mb*fGasThick;
0157
0158   G4complex Ca(1.0+0.5*fPlateThick*Ma/fAlphaPlate,fPlateThick/Za/fAlphaPlate);
0159   G4complex Cb(1.0+0.5*fGasThick*Mb/fAlphaGas,fGasThick/Zb/fAlphaGas);
0160
0161   G4complex Ha = pow(Ca,-fAlphaPlate);
0162   G4complex Hb = pow(Cb,-fAlphaGas);
0163   G4complex H  = Ha*Hb;
0164   G4complex F1 =   (1.0 - Ha)*(1.0 - Hb )/(1.0 - H)
0165                  * G4double(fPlateNumber) ;
0166   G4complex F2 =   (1.0-Ha)*(1.0-Ha)*Hb/(1.0-H)/(1.0-H)
0167                  * (1.0 - exp(-0.5*fPlateNumber*sigma)) ;
0168   //    *(1.0 - pow(H,fPlateNumber)) ;
0169     G4complex R  = (F1 + F2)*OneInterfaceXTRdEdx(energy,gamma,varAngle);
0170   // G4complex R  = F2*OneInterfaceXTRdEdx(energy,gamma,varAngle);
0171   result       = 2.0*real(R);
0172   return      result;
0173   */
0174   // numerically unstable result
0175
0176   G4double result, Qa, Qb, Q, aZa, bZb, aMa, bMb, D, sigma;
0177
0178   aZa = fPlateThick / GetPlateFormationZone(energy, gamma, varAngle);
0179   bZb = fGasThick / GetGasFormationZone(energy, gamma, varAngle);
0180   aMa = fPlateThick * GetPlateLinearPhotoAbs(energy);
0181   bMb = fGasThick * GetGasLinearPhotoAbs(energy);
0182   sigma = aMa * fPlateThick + bMb * fGasThick;
0183   Qa = exp(-0.5 * aMa);
0184   Qb = exp(-0.5 * bMb);
0185   Q = Qa * Qb;
0186
0187   G4complex Ha(Qa * cos(aZa), -Qa * sin(aZa));
0188   G4complex Hb(Qb * cos(bZb), -Qb * sin(bZb));
0189   G4complex H = Ha * Hb;
0190   G4complex Hs = conj(H);
0191   D = 1.0 / ((1 - Q) * (1 - Q) + 4 * Q * sin(0.5 * (aZa + bZb)) * sin(0.5 * (aZa + bZb)));
0192   G4complex F1 = (1.0 - Ha) * (1.0 - Hb) * (1.0 - Hs) * G4double(fPlateNumber) * D;
0193   G4complex F2 = (1.0 - Ha) * (1.0 - Ha) * Hb * (1.0 - Hs)
0194                  * (1.0 - Hs)
0195                  // * (1.0 - pow(H,fPlateNumber)) * D*D;
0196                  * (1.0 - exp(-0.5 * fPlateNumber * sigma)) * D * D;
0197   G4complex R = (F1 + F2) * OneInterfaceXTRdEdx(energy, gamma, varAngle);
0198   result = 2.0 * real(R);
0199   return result;
0200 }