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 /// \file XTRTransparentRegRadModel.cc
0027 /// \brief Implementation of the XTRTransparentRegRadModel class
0028
0029 #include "XTRTransparentRegRadModel.hh"
0030
0031 #include "G4Gamma.hh"
0032 #include "G4Integrator.hh"
0033 #include "G4PhysicalConstants.hh"
0034 #include "Randomize.hh"
0035
0036 #include <complex>
0037
0038 using namespace std;
0039
0040 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0041
0042 ////////////////////////////////////////////////////////////////////////////
0043 //
0044 // Constructor, destructor
0045
0046 XTRTransparentRegRadModel::XTRTransparentRegRadModel(G4LogicalVolume* anEnvelope,
0047                                                      G4Material* foilMat, G4Material* gasMat,
0048                                                      G4double a, G4double b, G4int n,
0049                                                      const G4String& processName)
0050   : G4VXTRenergyLoss(anEnvelope, foilMat, gasMat, a, b, n, processName)
0051 {
0052   G4cout << "Regular transparent X-ray TR  radiator EM process is called" << G4endl;
0053
0054   // Build energy and angular integral spectra of X-ray TR photons from
0055   // a radiator
0056   fExitFlux = true;
0057   fAlphaPlate = 10000;
0058   fAlphaGas = 1000;
0059
0060   //  BuildTable();
0061 }
0062
0063 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0064
0065 XTRTransparentRegRadModel::~XTRTransparentRegRadModel()
0066 {
0067   ;
0068 }
0069
0070 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0071
0072 G4double XTRTransparentRegRadModel::SpectralXTRdEdx(G4double energy)
0073 {
0074   G4double result, sum = 0., tmp, cof1, cof2, cofMin, cofPHC, aMa, bMb, sigma;
0075   G4int k, kMax, kMin;
0076
0077   aMa = GetPlateLinearPhotoAbs(energy);
0078   bMb = GetGasLinearPhotoAbs(energy);
0079
0080   // if(fCompton)
0081   {
0082     aMa += GetPlateCompton(energy);
0083     bMb += GetGasCompton(energy);
0084   }
0085   aMa *= fPlateThick;
0086   bMb *= fGasThick;
0087
0088   sigma = aMa + bMb;
0089
0090   cofPHC = 4 * pi * hbarc;
0091   cofPHC *= 200. / 197.;
0092   tmp = (fSigma1 - fSigma2) / cofPHC / energy;
0093   cof1 = fPlateThick * tmp;
0094   cof2 = fGasThick * tmp;
0095
0096   cofMin = energy * (fPlateThick + fGasThick) / fGamma / fGamma;
0097   cofMin += (fPlateThick * fSigma1 + fGasThick * fSigma2) / energy;
0098   cofMin /= cofPHC;
0099
0100   //  if (fGamma < 1200) kMin = G4int(cofMin);  // 1200 ?
0101   // else               kMin = 1;
0102
0103   kMin = G4int(cofMin);
0104   if (cofMin > kMin) kMin++;
0105
0106   // tmp  = (fPlateThick + fGasThick)*energy*fMaxThetaTR;
0107   // tmp /= cofPHC;
0108   // kMax = G4int(tmp);
0109   // if(kMax < 0) kMax = 0;
0110   // kMax += kMin;
0111
0112   kMax = kMin + 9;  // 5; // 9; //   kMin + G4int(tmp);
0113
0114   // tmp /= fGamma;
0115   // if( G4int(tmp) < kMin ) kMin = G4int(tmp);
0116   // G4cout<<"kMin = "<<kMin<<";    kMax = "<<kMax<<G4endl;
0117
0118   for (k = kMin; k <= kMax; k++) {
0119     tmp = pi * fPlateThick * (k + cof2) / (fPlateThick + fGasThick);
0120     result = (k - cof1) * (k - cof1) * (k + cof2) * (k + cof2);
0121
0122     if (k == kMin && kMin == G4int(cofMin)) {
0123       sum += 0.5 * sin(tmp) * sin(tmp) * std::abs(k - cofMin) / result;
0124     }
0125     else {
0126       sum += sin(tmp) * sin(tmp) * std::abs(k - cofMin) / result;
0127     }
0128     //  G4cout<<"k = "<<k<<";    sum = "<<sum<<G4endl;
0129   }
0130   result = 4. * (cof1 + cof2) * (cof1 + cof2) * sum / energy;
0131   result *= (1. - exp(-fPlateNumber * sigma)) / (1. - exp(-sigma));
0132   return result;
0133 }
0134
0135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0136 //
0137 // Approximation for radiator interference factor for the case of
0138 // fully Regular radiator. The plate and gas gap thicknesses are fixed .
0139 // The mean values of the plate and gas gap thicknesses
0140 // are supposed to be about XTR formation zones but much less than
0141 // mean absorption length of XTR photons in coresponding material.
0142
0143 G4double XTRTransparentRegRadModel::GetStackFactor(G4double energy, G4double gamma,
0144                                                    G4double varAngle)
0145 {
0146   /*
0147   G4double result, Za, Zb, Ma, Mb, sigma;
0148
0149   Za = GetPlateFormationZone(energy,gamma,varAngle);
0150   Zb = GetGasFormationZone(energy,gamma,varAngle);
0151   Ma = GetPlateLinearPhotoAbs(energy);
0152   Mb = GetGasLinearPhotoAbs(energy);
0153   sigma = Ma*fPlateThick + Mb*fGasThick;
0154
0155   G4complex Ca(1.0+0.5*fPlateThick*Ma/fAlphaPlate,fPlateThick/Za/fAlphaPlate);
0156   G4complex Cb(1.0+0.5*fGasThick*Mb/fAlphaGas,fGasThick/Zb/fAlphaGas);
0157
0158   G4complex Ha = pow(Ca,-fAlphaPlate);
0159   G4complex Hb = pow(Cb,-fAlphaGas);
0160   G4complex H  = Ha*Hb;
0161   G4complex F1 =   (1.0 - Ha)*(1.0 - Hb )/(1.0 - H)
0162                  * G4double(fPlateNumber) ;
0163   G4complex F2 =   (1.0-Ha)*(1.0-Ha)*Hb/(1.0-H)/(1.0-H)
0164                  * (1.0 - exp(-0.5*fPlateNumber*sigma)) ;
0165   //    *(1.0 - pow(H,fPlateNumber)) ;
0166     G4complex R  = (F1 + F2)*OneInterfaceXTRdEdx(energy,gamma,varAngle);
0167   // G4complex R  = F2*OneInterfaceXTRdEdx(energy,gamma,varAngle);
0168   result       = 2.0*real(R);
0169   return      result;
0170   */
0171   // numerically unstable result
0172
0173   G4double result, Qa, Qb, Q, aZa, bZb, aMa, bMb, D, sigma;
0174
0175   aZa = fPlateThick / GetPlateFormationZone(energy, gamma, varAngle);
0176   bZb = fGasThick / GetGasFormationZone(energy, gamma, varAngle);
0177   aMa = fPlateThick * GetPlateLinearPhotoAbs(energy);
0178   bMb = fGasThick * GetGasLinearPhotoAbs(energy);
0179   sigma = aMa * fPlateThick + bMb * fGasThick;
0180   Qa = exp(-0.5 * aMa);
0181   Qb = exp(-0.5 * bMb);
0182   Q = Qa * Qb;
0183
0184   G4complex Ha(Qa * cos(aZa), -Qa * sin(aZa));
0185   G4complex Hb(Qb * cos(bZb), -Qb * sin(bZb));
0186   G4complex H = Ha * Hb;
0187   G4complex Hs = conj(H);
0188   D = 1.0 / ((1 - Q) * (1 - Q) + 4 * Q * sin(0.5 * (aZa + bZb)) * sin(0.5 * (aZa + bZb)));
0189   G4complex F1 = (1.0 - Ha) * (1.0 - Hb) * (1.0 - Hs) * G4double(fPlateNumber) * D;
0190   G4complex F2 = (1.0 - Ha) * (1.0 - Ha) * Hb * (1.0 - Hs)
0191                  * (1.0 - Hs)
0192                  // * (1.0 - pow(H,fPlateNumber)) * D*D;
0193                  * (1.0 - exp(-0.5 * fPlateNumber * sigma)) * D * D;
0194   G4complex R = (F1 + F2) * OneInterfaceXTRdEdx(energy, gamma, varAngle);
0195   result = 2.0 * real(R);
0196   return result;
0197 }