Geant4 Cross Reference |
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 26 27 #include "G4RegularXTRadiator.hh" 27 #include "G4RegularXTRadiator.hh" 28 28 29 #include "G4Gamma.hh" 29 #include "G4Gamma.hh" 30 #include "G4PhysicalConstants.hh" 30 #include "G4PhysicalConstants.hh" 31 31 32 ////////////////////////////////////////////// 32 //////////////////////////////////////////////////////////////////////////// 33 // Constructor, destructor 33 // Constructor, destructor 34 G4RegularXTRadiator::G4RegularXTRadiator(G4Log 34 G4RegularXTRadiator::G4RegularXTRadiator(G4LogicalVolume* anEnvelope, 35 G4Mat 35 G4Material* foilMat, 36 G4Mat 36 G4Material* gasMat, G4double a, 37 G4dou 37 G4double b, G4int n, 38 const 38 const G4String& processName) 39 : G4VXTRenergyLoss(anEnvelope, foilMat, gasM 39 : G4VXTRenergyLoss(anEnvelope, foilMat, gasMat, a, b, n, processName) 40 { 40 { 41 G4cout << "Regular X-ray TR radiator EM proc 41 G4cout << "Regular X-ray TR radiator EM process is called" << G4endl; 42 42 43 // Build energy and angular integral spectra 43 // Build energy and angular integral spectra of X-ray TR photons from 44 // a radiator 44 // a radiator 45 45 46 fAlphaPlate = 10000; 46 fAlphaPlate = 10000; 47 fAlphaGas = 1000; 47 fAlphaGas = 1000; 48 G4cout << "fAlphaPlate = " << fAlphaPlate << 48 G4cout << "fAlphaPlate = " << fAlphaPlate << " ; fAlphaGas = " << fAlphaGas 49 << G4endl; 49 << G4endl; 50 } 50 } 51 51 52 ////////////////////////////////////////////// 52 /////////////////////////////////////////////////////////////////////////// 53 G4RegularXTRadiator::~G4RegularXTRadiator() = 53 G4RegularXTRadiator::~G4RegularXTRadiator() = default; 54 54 55 void G4RegularXTRadiator::ProcessDescription(s 55 void G4RegularXTRadiator::ProcessDescription(std::ostream& out) const 56 { 56 { 57 out << "Simulation of X-ray transition radia 57 out << "Simulation of X-ray transition radiation generated by\n" 58 "relativistic charged particles cross 58 "relativistic charged particles crossing the interface between\n" 59 "two materials. Thicknesses of plates 59 "two materials. Thicknesses of plates and gaps are fixed.\n"; 60 } 60 } 61 61 62 ////////////////////////////////////////////// 62 /////////////////////////////////////////////////////////////////////////// 63 G4double G4RegularXTRadiator::SpectralXTRdEdx( 63 G4double G4RegularXTRadiator::SpectralXTRdEdx(G4double energy) 64 { 64 { 65 G4double result, sum = 0., tmp, cof1, cof2, 65 G4double result, sum = 0., tmp, cof1, cof2, cofMin, cofPHC, theta2, theta2k; 66 G4double aMa, bMb, sigma, dump; 66 G4double aMa, bMb, sigma, dump; 67 G4int k, kMax, kMin; 67 G4int k, kMax, kMin; 68 68 69 aMa = fPlateThick * GetPlateLinearPhotoAbs 69 aMa = fPlateThick * GetPlateLinearPhotoAbs(energy); 70 bMb = fGasThick * GetGasLinearPhotoAbs(ene 70 bMb = fGasThick * GetGasLinearPhotoAbs(energy); 71 sigma = 0.5 * (aMa + bMb); 71 sigma = 0.5 * (aMa + bMb); 72 dump = std::exp(-fPlateNumber * sigma); 72 dump = std::exp(-fPlateNumber * sigma); 73 if(verboseLevel > 2) 73 if(verboseLevel > 2) 74 G4cout << " dump = " << dump << G4endl; 74 G4cout << " dump = " << dump << G4endl; 75 cofPHC = 4 * pi * hbarc; 75 cofPHC = 4 * pi * hbarc; 76 tmp = (fSigma1 - fSigma2) / cofPHC / ener 76 tmp = (fSigma1 - fSigma2) / cofPHC / energy; 77 cof1 = fPlateThick * tmp; 77 cof1 = fPlateThick * tmp; 78 cof2 = fGasThick * tmp; 78 cof2 = fGasThick * tmp; 79 79 80 cofMin = energy * (fPlateThick + fGasThick) 80 cofMin = energy * (fPlateThick + fGasThick) / fGamma / fGamma; 81 cofMin += (fPlateThick * fSigma1 + fGasThick 81 cofMin += (fPlateThick * fSigma1 + fGasThick * fSigma2) / energy; 82 cofMin /= cofPHC; 82 cofMin /= cofPHC; 83 83 84 theta2 = cofPHC / (energy * (fPlateThick + f 84 theta2 = cofPHC / (energy * (fPlateThick + fGasThick)); 85 85 86 kMin = G4int(cofMin); 86 kMin = G4int(cofMin); 87 if(cofMin > kMin) 87 if(cofMin > kMin) 88 kMin++; 88 kMin++; 89 89 90 kMax = kMin + 49; 90 kMax = kMin + 49; 91 91 92 if(verboseLevel > 2) 92 if(verboseLevel > 2) 93 { 93 { 94 G4cout << cof1 << " " << cof2 << " 94 G4cout << cof1 << " " << cof2 << " " << cofMin << G4endl; 95 G4cout << "kMin = " << kMin << "; kMax 95 G4cout << "kMin = " << kMin << "; kMax = " << kMax << G4endl; 96 } 96 } 97 for(k = kMin; k <= kMax; ++k) 97 for(k = kMin; k <= kMax; ++k) 98 { 98 { 99 tmp = pi * fPlateThick * (k + cof2) / ( 99 tmp = pi * fPlateThick * (k + cof2) / (fPlateThick + fGasThick); 100 result = (k - cof1) * (k - cof1) * (k + co 100 result = (k - cof1) * (k - cof1) * (k + cof2) * (k + cof2); 101 if(k == kMin && kMin == G4int(cofMin)) 101 if(k == kMin && kMin == G4int(cofMin)) 102 { 102 { 103 sum += 103 sum += 104 0.5 * std::sin(tmp) * std::sin(tmp) * 104 0.5 * std::sin(tmp) * std::sin(tmp) * std::abs(k - cofMin) / result; 105 } 105 } 106 else 106 else 107 { 107 { 108 sum += std::sin(tmp) * std::sin(tmp) * s 108 sum += std::sin(tmp) * std::sin(tmp) * std::abs(k - cofMin) / result; 109 } 109 } 110 theta2k = std::sqrt(theta2 * std::abs(k - 110 theta2k = std::sqrt(theta2 * std::abs(k - cofMin)); 111 111 112 if(verboseLevel > 2) 112 if(verboseLevel > 2) 113 { 113 { 114 G4cout << k << " " << theta2k << " 114 G4cout << k << " " << theta2k << " " 115 << std::sin(tmp) * std::sin(tmp) 115 << std::sin(tmp) * std::sin(tmp) * std::abs(k - cofMin) / result 116 << " " << sum << G4endl; 116 << " " << sum << G4endl; 117 } 117 } 118 } 118 } 119 result = 2 * (cof1 + cof2) * (cof1 + cof2) * 119 result = 2 * (cof1 + cof2) * (cof1 + cof2) * sum / energy; 120 result *= (1 - dump + 2 * dump * fPlateNumbe 120 result *= (1 - dump + 2 * dump * fPlateNumber); 121 121 122 return result; 122 return result; 123 } 123 } 124 124 125 ////////////////////////////////////////////// 125 /////////////////////////////////////////////////////////////////////////// 126 // Approximation for radiator interference fac 126 // Approximation for radiator interference factor for the case of 127 // fully Regular radiator. The plate and gas g 127 // fully Regular radiator. The plate and gas gap thicknesses are fixed. 128 // The mean values of the plate and gas gap th 128 // The mean values of the plate and gas gap thicknesses 129 // are supposed to be about XTR formation zone 129 // are supposed to be about XTR formation zones but much less than 130 // mean absorption length of XTR photons in co 130 // mean absorption length of XTR photons in corresponding material. 131 131 132 G4double G4RegularXTRadiator::GetStackFactor(G 132 G4double G4RegularXTRadiator::GetStackFactor(G4double energy, G4double gamma, 133 G 133 G4double varAngle) 134 { 134 { 135 // some gamma (10000/1000) like algorithm 135 // some gamma (10000/1000) like algorithm 136 136 137 G4double result, Za, Zb, Ma, Mb; 137 G4double result, Za, Zb, Ma, Mb; 138 138 139 Za = GetPlateFormationZone(energy, gamma, va 139 Za = GetPlateFormationZone(energy, gamma, varAngle); 140 Zb = GetGasFormationZone(energy, gamma, varA 140 Zb = GetGasFormationZone(energy, gamma, varAngle); 141 141 142 Ma = GetPlateLinearPhotoAbs(energy); 142 Ma = GetPlateLinearPhotoAbs(energy); 143 Mb = GetGasLinearPhotoAbs(energy); 143 Mb = GetGasLinearPhotoAbs(energy); 144 144 145 G4complex Ca(1.0 + 0.5 * fPlateThick * Ma / 145 G4complex Ca(1.0 + 0.5 * fPlateThick * Ma / fAlphaPlate, 146 fPlateThick / Za / fAlphaPlate) 146 fPlateThick / Za / fAlphaPlate); 147 G4complex Cb(1.0 + 0.5 * fGasThick * Mb / fA 147 G4complex Cb(1.0 + 0.5 * fGasThick * Mb / fAlphaGas, 148 fGasThick / Zb / fAlphaGas); 148 fGasThick / Zb / fAlphaGas); 149 149 150 G4complex Ha = std::pow(Ca, -fAlphaPlate); 150 G4complex Ha = std::pow(Ca, -fAlphaPlate); 151 G4complex Hb = std::pow(Cb, -fAlphaGas); 151 G4complex Hb = std::pow(Cb, -fAlphaGas); 152 G4complex H = Ha * Hb; 152 G4complex H = Ha * Hb; 153 153 154 G4complex F1 = (1.0 - Ha) * (1.0 - Hb) / (1. 154 G4complex F1 = (1.0 - Ha) * (1.0 - Hb) / (1.0 - H) * G4double(fPlateNumber); 155 155 156 G4complex F2 = (1.0 - Ha) * (1.0 - Ha) * Hb 156 G4complex F2 = (1.0 - Ha) * (1.0 - Ha) * Hb / (1.0 - H) / (1.0 - H) * 157 (1.0 - std::pow(H, fPlateNumb 157 (1.0 - std::pow(H, fPlateNumber)); 158 158 159 G4complex R = (F1 + F2) * OneInterfaceXTRdEd 159 G4complex R = (F1 + F2) * OneInterfaceXTRdEdx(energy, gamma, varAngle); 160 160 161 result = 2.0 * std::real(R); 161 result = 2.0 * std::real(R); 162 162 163 return result; 163 return result; 164 } 164 } 165 165