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 // History: 26 // History: 27 // 2001-2002 R&D by V.Grichine 27 // 2001-2002 R&D by V.Grichine 28 // 19.06.03 V. Grichine, modifications in Buil 28 // 19.06.03 V. Grichine, modifications in BuildTable for the integration 29 // in respect of angle: 29 // in respect of angle: range is increased, accuracy is 30 // improved 30 // improved 31 // 28.07.05, P.Gumplinger add G4ProcessType to 31 // 28.07.05, P.Gumplinger add G4ProcessType to constructor 32 // 28.09.07, V.Ivanchenko general cleanup with 32 // 28.09.07, V.Ivanchenko general cleanup without change of algorithms 33 // 33 // 34 34 35 #include "G4VXTRenergyLoss.hh" 35 #include "G4VXTRenergyLoss.hh" 36 36 37 #include "G4AffineTransform.hh" 37 #include "G4AffineTransform.hh" 38 #include "G4DynamicParticle.hh" 38 #include "G4DynamicParticle.hh" 39 #include "G4EmProcessSubType.hh" 39 #include "G4EmProcessSubType.hh" 40 #include "G4Integrator.hh" 40 #include "G4Integrator.hh" 41 #include "G4MaterialTable.hh" 41 #include "G4MaterialTable.hh" 42 #include "G4ParticleMomentum.hh" 42 #include "G4ParticleMomentum.hh" 43 #include "G4PhysicalConstants.hh" 43 #include "G4PhysicalConstants.hh" 44 #include "G4PhysicsFreeVector.hh" 44 #include "G4PhysicsFreeVector.hh" 45 #include "G4PhysicsLinearVector.hh" 45 #include "G4PhysicsLinearVector.hh" 46 #include "G4PhysicsLogVector.hh" 46 #include "G4PhysicsLogVector.hh" 47 #include "G4RotationMatrix.hh" 47 #include "G4RotationMatrix.hh" 48 #include "G4SandiaTable.hh" 48 #include "G4SandiaTable.hh" 49 #include "G4SystemOfUnits.hh" 49 #include "G4SystemOfUnits.hh" 50 #include "G4ThreeVector.hh" 50 #include "G4ThreeVector.hh" 51 #include "G4Timer.hh" 51 #include "G4Timer.hh" 52 #include "G4VDiscreteProcess.hh" 52 #include "G4VDiscreteProcess.hh" 53 #include "G4VParticleChange.hh" 53 #include "G4VParticleChange.hh" 54 #include "G4VSolid.hh" 54 #include "G4VSolid.hh" 55 #include "G4PhysicsModelCatalog.hh" 55 #include "G4PhysicsModelCatalog.hh" 56 56 57 ////////////////////////////////////////////// 57 //////////////////////////////////////////////////////////////////////////// 58 // Constructor, destructor 58 // Constructor, destructor 59 G4VXTRenergyLoss::G4VXTRenergyLoss(G4LogicalVo 59 G4VXTRenergyLoss::G4VXTRenergyLoss(G4LogicalVolume* anEnvelope, 60 G4Material* 60 G4Material* foilMat, G4Material* gasMat, 61 G4double a, 61 G4double a, G4double b, G4int n, 62 const G4Str 62 const G4String& processName, 63 G4ProcessTy 63 G4ProcessType type) 64 : G4VDiscreteProcess(processName, type) 64 : G4VDiscreteProcess(processName, type) 65 , fGammaCutInKineticEnergy(nullptr) 65 , fGammaCutInKineticEnergy(nullptr) 66 , fAngleDistrTable(nullptr) 66 , fAngleDistrTable(nullptr) 67 , fEnergyDistrTable(nullptr) 67 , fEnergyDistrTable(nullptr) 68 , fAngleForEnergyTable(nullptr) 68 , fAngleForEnergyTable(nullptr) 69 , fPlatePhotoAbsCof(nullptr) 69 , fPlatePhotoAbsCof(nullptr) 70 , fGasPhotoAbsCof(nullptr) 70 , fGasPhotoAbsCof(nullptr) 71 , fGammaTkinCut(0.0) 71 , fGammaTkinCut(0.0) 72 { 72 { 73 verboseLevel = 1; 73 verboseLevel = 1; 74 secID = G4PhysicsModelCatalog::GetModelID("m 74 secID = G4PhysicsModelCatalog::GetModelID("model_XTRenergyLoss"); 75 SetProcessSubType(fTransitionRadiation); 75 SetProcessSubType(fTransitionRadiation); 76 76 77 fPtrGamma = nullptr; 77 fPtrGamma = nullptr; 78 fMinEnergyTR = fMaxEnergyTR = fMaxThetaTR = 78 fMinEnergyTR = fMaxEnergyTR = fMaxThetaTR = fGamma = fEnergy = 0.0; 79 fVarAngle = fLambda = fTotalDist = fPlateThi 79 fVarAngle = fLambda = fTotalDist = fPlateThick = fGasThick = 0.0; 80 fAlphaPlate = 100.; 80 fAlphaPlate = 100.; 81 fAlphaGas = 40.; 81 fAlphaGas = 40.; 82 82 83 fTheMinEnergyTR = CLHEP::keV * 1.; // 1.; / 83 fTheMinEnergyTR = CLHEP::keV * 1.; // 1.; // 84 fTheMaxEnergyTR = CLHEP::keV * 100.; // 40.; 84 fTheMaxEnergyTR = CLHEP::keV * 100.; // 40.; // 85 85 86 fTheMinAngle = 1.e-8; // 86 fTheMinAngle = 1.e-8; // 87 fTheMaxAngle = 4.e-4; 87 fTheMaxAngle = 4.e-4; 88 88 89 fTotBin = 50; // number of bins in log sca 89 fTotBin = 50; // number of bins in log scale 90 fBinTR = 100; // number of bins in TR vec << 90 fBinTR = 100; // number of bins in TR vectors 91 fKrange = 229; << 91 92 // min/max angle2 in log-vectors 92 // min/max angle2 in log-vectors 93 93 94 fMinThetaTR = 3.0e-9; 94 fMinThetaTR = 3.0e-9; 95 fMaxThetaTR = 1.0e-4; 95 fMaxThetaTR = 1.0e-4; 96 96 97 97 98 // Proton energy vector initialization 98 // Proton energy vector initialization 99 fProtonEnergyVector = 99 fProtonEnergyVector = 100 new G4PhysicsLogVector(fMinProtonTkin, fMa 100 new G4PhysicsLogVector(fMinProtonTkin, fMaxProtonTkin, fTotBin); 101 101 102 fXTREnergyVector = 102 fXTREnergyVector = 103 new G4PhysicsLogVector(fTheMinEnergyTR, fT 103 new G4PhysicsLogVector(fTheMinEnergyTR, fTheMaxEnergyTR, fBinTR); 104 104 105 fEnvelope = anEnvelope; 105 fEnvelope = anEnvelope; 106 106 107 fPlateNumber = n; 107 fPlateNumber = n; 108 if(verboseLevel > 0) 108 if(verboseLevel > 0) 109 G4cout << "### G4VXTRenergyLoss: the numbe 109 G4cout << "### G4VXTRenergyLoss: the number of TR radiator plates = " 110 << fPlateNumber << G4endl; 110 << fPlateNumber << G4endl; 111 if(fPlateNumber == 0) 111 if(fPlateNumber == 0) 112 { 112 { 113 G4Exception("G4VXTRenergyLoss::G4VXTRenerg 113 G4Exception("G4VXTRenergyLoss::G4VXTRenergyLoss()", "VXTRELoss01", 114 FatalException, "No plates in 114 FatalException, "No plates in X-ray TR radiator"); 115 } 115 } 116 // default is XTR dEdx, not flux after radia 116 // default is XTR dEdx, not flux after radiator 117 fExitFlux = false; 117 fExitFlux = false; 118 // default angle distribution according nume 118 // default angle distribution according numerical integration 119 fFastAngle = false; // no angle accordin 119 fFastAngle = false; // no angle according sum of delta-functions by default 120 fAngleRadDistr = true; 120 fAngleRadDistr = true; 121 fCompton = false; 121 fCompton = false; 122 122 123 fLambda = DBL_MAX; 123 fLambda = DBL_MAX; 124 124 125 // Mean thicknesses of plates and gas gaps 125 // Mean thicknesses of plates and gas gaps 126 fPlateThick = a; 126 fPlateThick = a; 127 fGasThick = b; 127 fGasThick = b; 128 fTotalDist = fPlateNumber * (fPlateThick + 128 fTotalDist = fPlateNumber * (fPlateThick + fGasThick); 129 if(verboseLevel > 0) 129 if(verboseLevel > 0) 130 G4cout << "total radiator thickness = " << 130 G4cout << "total radiator thickness = " << fTotalDist / cm << " cm" 131 << G4endl; 131 << G4endl; 132 132 133 // index of plate material 133 // index of plate material 134 fMatIndex1 = (G4int)foilMat->GetIndex(); << 134 fMatIndex1 = foilMat->GetIndex(); 135 if(verboseLevel > 0) 135 if(verboseLevel > 0) 136 G4cout << "plate material = " << foilMat-> 136 G4cout << "plate material = " << foilMat->GetName() << G4endl; 137 137 138 // index of gas material 138 // index of gas material 139 fMatIndex2 = (G4int)gasMat->GetIndex(); << 139 fMatIndex2 = gasMat->GetIndex(); 140 if(verboseLevel > 0) 140 if(verboseLevel > 0) 141 G4cout << "gas material = " << gasMat->Get 141 G4cout << "gas material = " << gasMat->GetName() << G4endl; 142 142 143 // plasma energy squared for plate material 143 // plasma energy squared for plate material 144 fSigma1 = fPlasmaCof * foilMat->GetElectronD 144 fSigma1 = fPlasmaCof * foilMat->GetElectronDensity(); 145 if(verboseLevel > 0) 145 if(verboseLevel > 0) 146 G4cout << "plate plasma energy = " << std: 146 G4cout << "plate plasma energy = " << std::sqrt(fSigma1) / eV << " eV" 147 << G4endl; 147 << G4endl; 148 148 149 // plasma energy squared for gas material 149 // plasma energy squared for gas material 150 fSigma2 = fPlasmaCof * gasMat->GetElectronDe 150 fSigma2 = fPlasmaCof * gasMat->GetElectronDensity(); 151 if(verboseLevel > 0) 151 if(verboseLevel > 0) 152 G4cout << "gas plasma energy = " << std::s 152 G4cout << "gas plasma energy = " << std::sqrt(fSigma2) / eV << " eV" 153 << G4endl; 153 << G4endl; 154 154 155 // Compute cofs for preparation of linear ph 155 // Compute cofs for preparation of linear photo absorption 156 ComputePlatePhotoAbsCof(); 156 ComputePlatePhotoAbsCof(); 157 ComputeGasPhotoAbsCof(); 157 ComputeGasPhotoAbsCof(); 158 158 159 pParticleChange = &fParticleChange; 159 pParticleChange = &fParticleChange; 160 } 160 } 161 161 162 ////////////////////////////////////////////// 162 /////////////////////////////////////////////////////////////////////////// 163 G4VXTRenergyLoss::~G4VXTRenergyLoss() 163 G4VXTRenergyLoss::~G4VXTRenergyLoss() 164 { 164 { 165 delete fProtonEnergyVector; 165 delete fProtonEnergyVector; 166 delete fXTREnergyVector; 166 delete fXTREnergyVector; 167 if(fEnergyDistrTable) 167 if(fEnergyDistrTable) 168 { 168 { 169 fEnergyDistrTable->clearAndDestroy(); 169 fEnergyDistrTable->clearAndDestroy(); 170 delete fEnergyDistrTable; 170 delete fEnergyDistrTable; 171 } 171 } 172 if(fAngleRadDistr) 172 if(fAngleRadDistr) 173 { 173 { 174 fAngleDistrTable->clearAndDestroy(); 174 fAngleDistrTable->clearAndDestroy(); 175 delete fAngleDistrTable; 175 delete fAngleDistrTable; 176 } 176 } 177 if(fAngleForEnergyTable) 177 if(fAngleForEnergyTable) 178 { 178 { 179 fAngleForEnergyTable->clearAndDestroy(); 179 fAngleForEnergyTable->clearAndDestroy(); 180 delete fAngleForEnergyTable; 180 delete fAngleForEnergyTable; 181 } 181 } 182 } 182 } 183 183 184 void G4VXTRenergyLoss::ProcessDescription(std: 184 void G4VXTRenergyLoss::ProcessDescription(std::ostream& out) const 185 { 185 { 186 out << "Base class for 'fast' parameterisati 186 out << "Base class for 'fast' parameterisation model describing X-ray " 187 "transition\n" 187 "transition\n" 188 "radiation. Angular distribution is v 188 "radiation. Angular distribution is very rough.\n"; 189 } 189 } 190 190 191 ////////////////////////////////////////////// 191 /////////////////////////////////////////////////////////////////////////////// 192 // Returns condition for application of the mo 192 // Returns condition for application of the model depending on particle type 193 G4bool G4VXTRenergyLoss::IsApplicable(const G4 193 G4bool G4VXTRenergyLoss::IsApplicable(const G4ParticleDefinition& particle) 194 { 194 { 195 return (particle.GetPDGCharge() != 0.0); 195 return (particle.GetPDGCharge() != 0.0); 196 } 196 } 197 197 198 ////////////////////////////////////////////// 198 ///////////////////////////////////////////////////////////////////////////////// 199 // Calculate step size for XTR process inside 199 // Calculate step size for XTR process inside raaditor 200 G4double G4VXTRenergyLoss::GetMeanFreePath(con 200 G4double G4VXTRenergyLoss::GetMeanFreePath(const G4Track& aTrack, G4double, 201 G4F 201 G4ForceCondition* condition) 202 { 202 { 203 G4int iTkin, iPlace; 203 G4int iTkin, iPlace; 204 G4double lambda, sigma, kinEnergy, mass, gam 204 G4double lambda, sigma, kinEnergy, mass, gamma; 205 G4double charge, chargeSq, massRatio, TkinSc 205 G4double charge, chargeSq, massRatio, TkinScaled; 206 G4double E1, E2, W, W1, W2; 206 G4double E1, E2, W, W1, W2; 207 207 208 *condition = NotForced; 208 *condition = NotForced; 209 209 210 if(aTrack.GetVolume()->GetLogicalVolume() != 210 if(aTrack.GetVolume()->GetLogicalVolume() != fEnvelope) 211 lambda = DBL_MAX; 211 lambda = DBL_MAX; 212 else 212 else 213 { 213 { 214 const G4DynamicParticle* aParticle = aTrac 214 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 215 kinEnergy = aPart 215 kinEnergy = aParticle->GetKineticEnergy(); 216 mass = aParticle->GetDefinition()->GetPDG 216 mass = aParticle->GetDefinition()->GetPDGMass(); 217 gamma = 1.0 + kinEnergy / mass; 217 gamma = 1.0 + kinEnergy / mass; 218 if(verboseLevel > 1) 218 if(verboseLevel > 1) 219 { 219 { 220 G4cout << " gamma = " << gamma << "; f 220 G4cout << " gamma = " << gamma << "; fGamma = " << fGamma << G4endl; 221 } 221 } 222 222 223 if(std::fabs(gamma - fGamma) < 0.05 * gamm 223 if(std::fabs(gamma - fGamma) < 0.05 * gamma) 224 lambda = fLambda; 224 lambda = fLambda; 225 else 225 else 226 { 226 { 227 charge = aParticle->GetDefinition()- 227 charge = aParticle->GetDefinition()->GetPDGCharge(); 228 chargeSq = charge * charge; 228 chargeSq = charge * charge; 229 massRatio = proton_mass_c2 / mass; 229 massRatio = proton_mass_c2 / mass; 230 TkinScaled = kinEnergy * massRatio; 230 TkinScaled = kinEnergy * massRatio; 231 231 232 for(iTkin = 0; iTkin < fTotBin; ++iTkin) 232 for(iTkin = 0; iTkin < fTotBin; ++iTkin) 233 { 233 { 234 if(TkinScaled < fProtonEnergyVector->G 234 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) 235 break; 235 break; 236 } 236 } 237 iPlace = iTkin - 1; 237 iPlace = iTkin - 1; 238 238 239 if(iTkin == 0) 239 if(iTkin == 0) 240 lambda = DBL_MAX; // Tkin is too smal 240 lambda = DBL_MAX; // Tkin is too small, neglect of TR photon generation 241 else // general case: Tkin between two 241 else // general case: Tkin between two vectors of the material 242 { 242 { 243 if(iTkin == fTotBin) 243 if(iTkin == fTotBin) 244 { 244 { 245 sigma = (*(*fEnergyDistrTable)(iPlac 245 sigma = (*(*fEnergyDistrTable)(iPlace))(0) * chargeSq; 246 } 246 } 247 else 247 else 248 { 248 { 249 E1 = fProtonEnergyVector->GetLowE 249 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1); 250 E2 = fProtonEnergyVector->GetLowE 250 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin); 251 W = 1.0 / (E2 - E1); 251 W = 1.0 / (E2 - E1); 252 W1 = (E2 - TkinScaled) * W; 252 W1 = (E2 - TkinScaled) * W; 253 W2 = (TkinScaled - E1) * W; 253 W2 = (TkinScaled - E1) * W; 254 sigma = ((*(*fEnergyDistrTable)(iPla 254 sigma = ((*(*fEnergyDistrTable)(iPlace))(0) * W1 + 255 (*(*fEnergyDistrTable)(iPla 255 (*(*fEnergyDistrTable)(iPlace + 1))(0) * W2) * 256 chargeSq; 256 chargeSq; 257 } 257 } 258 if(sigma < DBL_MIN) 258 if(sigma < DBL_MIN) 259 lambda = DBL_MAX; 259 lambda = DBL_MAX; 260 else 260 else 261 lambda = 1. / sigma; 261 lambda = 1. / sigma; 262 fLambda = lambda; 262 fLambda = lambda; 263 fGamma = gamma; 263 fGamma = gamma; 264 if(verboseLevel > 1) 264 if(verboseLevel > 1) 265 { 265 { 266 G4cout << " lambda = " << lambda / m 266 G4cout << " lambda = " << lambda / mm << " mm" << G4endl; 267 } 267 } 268 } 268 } 269 } 269 } 270 } 270 } 271 return lambda; 271 return lambda; 272 } 272 } 273 273 274 ////////////////////////////////////////////// 274 ////////////////////////////////////////////////////////////////////////// 275 // Interface for build table from physics list 275 // Interface for build table from physics list 276 void G4VXTRenergyLoss::BuildPhysicsTable(const 276 void G4VXTRenergyLoss::BuildPhysicsTable(const G4ParticleDefinition& pd) 277 { 277 { 278 if(pd.GetPDGCharge() == 0.) 278 if(pd.GetPDGCharge() == 0.) 279 { 279 { 280 G4Exception("G4VXTRenergyLoss::BuildPhysic 280 G4Exception("G4VXTRenergyLoss::BuildPhysicsTable", "Notification", 281 JustWarning, "XTR initialisati 281 JustWarning, "XTR initialisation for neutral particle ?!"); 282 } 282 } 283 BuildEnergyTable(); 283 BuildEnergyTable(); 284 284 285 if(fAngleRadDistr) 285 if(fAngleRadDistr) 286 { 286 { 287 if(verboseLevel > 0) 287 if(verboseLevel > 0) 288 { 288 { 289 G4cout 289 G4cout 290 << "Build angle for energy distributio 290 << "Build angle for energy distribution according the current radiator" 291 << G4endl; 291 << G4endl; 292 } 292 } 293 BuildAngleForEnergyBank(); 293 BuildAngleForEnergyBank(); 294 } 294 } 295 } 295 } 296 296 297 ////////////////////////////////////////////// 297 ////////////////////////////////////////////////////////////////////////// 298 // Build integral energy distribution of XTR p 298 // Build integral energy distribution of XTR photons 299 void G4VXTRenergyLoss::BuildEnergyTable() 299 void G4VXTRenergyLoss::BuildEnergyTable() 300 { 300 { 301 G4int iTkin, iTR, iPlace; 301 G4int iTkin, iTR, iPlace; 302 G4double radiatorCof = 1.0; // for tuning o 302 G4double radiatorCof = 1.0; // for tuning of XTR yield 303 G4double energySum = 0.0; 303 G4double energySum = 0.0; 304 304 305 fEnergyDistrTable = new G4PhysicsTable(fTotB 305 fEnergyDistrTable = new G4PhysicsTable(fTotBin); 306 if(fAngleRadDistr) 306 if(fAngleRadDistr) 307 fAngleDistrTable = new G4PhysicsTable(fTot 307 fAngleDistrTable = new G4PhysicsTable(fTotBin); 308 308 309 fGammaTkinCut = 0.0; 309 fGammaTkinCut = 0.0; 310 310 311 // setting of min/max TR energies 311 // setting of min/max TR energies 312 if(fGammaTkinCut > fTheMinEnergyTR) 312 if(fGammaTkinCut > fTheMinEnergyTR) 313 fMinEnergyTR = fGammaTkinCut; 313 fMinEnergyTR = fGammaTkinCut; 314 else 314 else 315 fMinEnergyTR = fTheMinEnergyTR; 315 fMinEnergyTR = fTheMinEnergyTR; 316 316 317 if(fGammaTkinCut > fTheMaxEnergyTR) 317 if(fGammaTkinCut > fTheMaxEnergyTR) 318 fMaxEnergyTR = 2.0 * fGammaTkinCut; 318 fMaxEnergyTR = 2.0 * fGammaTkinCut; 319 else 319 else 320 fMaxEnergyTR = fTheMaxEnergyTR; 320 fMaxEnergyTR = fTheMaxEnergyTR; 321 321 322 G4Integrator<G4VXTRenergyLoss, G4double (G4V 322 G4Integrator<G4VXTRenergyLoss, G4double (G4VXTRenergyLoss::*)(G4double)> 323 integral; 323 integral; 324 324 325 G4cout.precision(4); 325 G4cout.precision(4); 326 G4Timer timer; 326 G4Timer timer; 327 timer.Start(); 327 timer.Start(); 328 328 329 if(verboseLevel > 0) 329 if(verboseLevel > 0) 330 { 330 { 331 G4cout << G4endl; 331 G4cout << G4endl; 332 G4cout << "Lorentz Factor" 332 G4cout << "Lorentz Factor" 333 << "\t" 333 << "\t" 334 << "XTR photon number" << G4endl; 334 << "XTR photon number" << G4endl; 335 G4cout << G4endl; 335 G4cout << G4endl; 336 } 336 } 337 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // 337 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // Lorentz factor loop 338 { 338 { 339 auto energyVector = << 339 G4PhysicsLogVector* energyVector = 340 new G4PhysicsLogVector(fMinEnergyTR, fMa 340 new G4PhysicsLogVector(fMinEnergyTR, fMaxEnergyTR, fBinTR); 341 341 342 fGamma = 342 fGamma = 343 1.0 + (fProtonEnergyVector->GetLowEdgeEn 343 1.0 + (fProtonEnergyVector->GetLowEdgeEnergy(iTkin) / proton_mass_c2); 344 344 345 // if(fMaxThetaTR > fTheMaxAngle) fMax << 345 fMaxThetaTR = 25. * 2500.0 / (fGamma * fGamma); // theta^2 346 // else if(fMaxThetaTR < fTheMinAngle) << 346 >> 347 if(fMaxThetaTR > fTheMaxAngle) >> 348 fMaxThetaTR = fTheMaxAngle; >> 349 else if(fMaxThetaTR < fTheMinAngle) >> 350 fMaxThetaTR = fTheMinAngle; 347 351 348 energySum = 0.0; 352 energySum = 0.0; 349 353 350 energyVector->PutValue(fBinTR - 1, energyS 354 energyVector->PutValue(fBinTR - 1, energySum); 351 355 352 for(iTR = fBinTR - 2; iTR >= 0; --iTR) 356 for(iTR = fBinTR - 2; iTR >= 0; --iTR) 353 { 357 { 354 // Legendre96 or Legendre10 358 // Legendre96 or Legendre10 355 359 356 energySum += radiatorCof * fCofTR * 360 energySum += radiatorCof * fCofTR * 357 361 358 // integral.Legendre10(this, &G4VXTRenergyLo 362 // integral.Legendre10(this, &G4VXTRenergyLoss::SpectralXTRdEdx, 359 363 360 integral.Legendre96(this, & 364 integral.Legendre96(this, &G4VXTRenergyLoss::SpectralXTRdEdx, 361 365 362 energyV 366 energyVector->GetLowEdgeEnergy(iTR), 363 energyV 367 energyVector->GetLowEdgeEnergy(iTR + 1)); 364 368 365 energyVector->PutValue(iTR, energySum / 369 energyVector->PutValue(iTR, energySum / fTotalDist); 366 } 370 } 367 iPlace = iTkin; 371 iPlace = iTkin; 368 fEnergyDistrTable->insertAt(iPlace, energy 372 fEnergyDistrTable->insertAt(iPlace, energyVector); 369 373 370 if(verboseLevel > 0) 374 if(verboseLevel > 0) 371 { 375 { 372 G4cout << fGamma << "\t" << energySum << 376 G4cout << fGamma << "\t" << energySum << G4endl; 373 } 377 } 374 } 378 } 375 timer.Stop(); 379 timer.Stop(); 376 G4cout.precision(6); 380 G4cout.precision(6); 377 if(verboseLevel > 0) 381 if(verboseLevel > 0) 378 { 382 { 379 G4cout << G4endl; 383 G4cout << G4endl; 380 G4cout << "total time for build X-ray TR e 384 G4cout << "total time for build X-ray TR energy loss tables = " 381 << timer.GetUserElapsed() << " s" < 385 << timer.GetUserElapsed() << " s" << G4endl; 382 } 386 } 383 fGamma = 0.; 387 fGamma = 0.; 384 return; 388 return; 385 } 389 } 386 390 387 ////////////////////////////////////////////// 391 ////////////////////////////////////////////////////////////////////////// 388 // Bank of angle distributions for given energ 392 // Bank of angle distributions for given energies (slow!) 389 393 390 void G4VXTRenergyLoss::BuildAngleForEnergyBank 394 void G4VXTRenergyLoss::BuildAngleForEnergyBank() 391 { 395 { 392 396 393 if( ( this->GetProcessName() == "TranspRegXT 397 if( ( this->GetProcessName() == "TranspRegXTRadiator" || 394 this->GetProcessName() == "TranspRegXT 398 this->GetProcessName() == "TranspRegXTRmodel" || 395 this->GetProcessName() == "RegularXTRa 399 this->GetProcessName() == "RegularXTRadiator" || 396 this->GetProcessName() == "RegularXTRmodel" 400 this->GetProcessName() == "RegularXTRmodel" ) && fFastAngle ) // ffastAngle=true! 397 { 401 { 398 BuildAngleTable(); // by sum of delta-func 402 BuildAngleTable(); // by sum of delta-functions 399 return; 403 return; 400 } 404 } 401 G4int i, iTkin, iTR; 405 G4int i, iTkin, iTR; 402 G4double angleSum = 0.0; 406 G4double angleSum = 0.0; 403 407 404 fGammaTkinCut = 0.0; 408 fGammaTkinCut = 0.0; 405 409 406 // setting of min/max TR energies 410 // setting of min/max TR energies 407 if(fGammaTkinCut > fTheMinEnergyTR) 411 if(fGammaTkinCut > fTheMinEnergyTR) 408 fMinEnergyTR = fGammaTkinCut; 412 fMinEnergyTR = fGammaTkinCut; 409 else 413 else 410 fMinEnergyTR = fTheMinEnergyTR; 414 fMinEnergyTR = fTheMinEnergyTR; 411 415 412 if(fGammaTkinCut > fTheMaxEnergyTR) 416 if(fGammaTkinCut > fTheMaxEnergyTR) 413 fMaxEnergyTR = 2.0 * fGammaTkinCut; 417 fMaxEnergyTR = 2.0 * fGammaTkinCut; 414 else 418 else 415 fMaxEnergyTR = fTheMaxEnergyTR; 419 fMaxEnergyTR = fTheMaxEnergyTR; 416 420 417 auto energyVector = << 421 G4PhysicsLogVector* energyVector = 418 new G4PhysicsLogVector(fMinEnergyTR, fMaxE 422 new G4PhysicsLogVector(fMinEnergyTR, fMaxEnergyTR, fBinTR); 419 423 420 G4Integrator<G4VXTRenergyLoss, G4double (G4V 424 G4Integrator<G4VXTRenergyLoss, G4double (G4VXTRenergyLoss::*)(G4double)> 421 integral; 425 integral; 422 426 423 G4cout.precision(4); 427 G4cout.precision(4); 424 G4Timer timer; 428 G4Timer timer; 425 timer.Start(); 429 timer.Start(); 426 430 427 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // 431 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // Lorentz factor loop 428 { 432 { 429 fGamma = 433 fGamma = 430 1.0 + (fProtonEnergyVector->GetLowEdgeEn 434 1.0 + (fProtonEnergyVector->GetLowEdgeEnergy(iTkin) / proton_mass_c2); 431 435 432 if(fMaxThetaTR > fTheMaxAngle) 436 if(fMaxThetaTR > fTheMaxAngle) 433 fMaxThetaTR = fTheMaxAngle; 437 fMaxThetaTR = fTheMaxAngle; 434 else if(fMaxThetaTR < fTheMinAngle) 438 else if(fMaxThetaTR < fTheMinAngle) 435 fMaxThetaTR = fTheMinAngle; 439 fMaxThetaTR = fTheMinAngle; 436 440 437 fAngleForEnergyTable = new G4PhysicsTable( 441 fAngleForEnergyTable = new G4PhysicsTable(fBinTR); 438 442 439 for(iTR = 0; iTR < fBinTR; ++iTR) 443 for(iTR = 0; iTR < fBinTR; ++iTR) 440 { 444 { 441 angleSum = 0.0; 445 angleSum = 0.0; 442 fEnergy = energyVector->GetLowEdgeEnerg 446 fEnergy = energyVector->GetLowEdgeEnergy(iTR); 443 447 444 // log-vector to increase number of thin 448 // log-vector to increase number of thin bins for small angles 445 auto angleVector = new G4PhysicsLogVecto << 449 G4PhysicsLogVector* angleVector = new G4PhysicsLogVector(fMinThetaTR, fMaxThetaTR, fBinTR); 446 450 447 451 448 452 449 angleVector->PutValue(fBinTR - 1, angleS 453 angleVector->PutValue(fBinTR - 1, angleSum); 450 454 451 for(i = fBinTR - 2; i >= 0; --i) 455 for(i = fBinTR - 2; i >= 0; --i) 452 { 456 { 453 // Legendre96 or Legendre10 457 // Legendre96 or Legendre10 454 458 455 angleSum += 459 angleSum += 456 integral.Legendre10(this, &G4VXTRene 460 integral.Legendre10(this, &G4VXTRenergyLoss::SpectralAngleXTRdEdx, 457 angleVector->Get 461 angleVector->GetLowEdgeEnergy(i), 458 angleVector->Get 462 angleVector->GetLowEdgeEnergy(i + 1)); 459 463 460 angleVector->PutValue(i, angleSum); 464 angleVector->PutValue(i, angleSum); 461 } 465 } 462 fAngleForEnergyTable->insertAt(iTR, angl 466 fAngleForEnergyTable->insertAt(iTR, angleVector); 463 } 467 } 464 fAngleBank.push_back(fAngleForEnergyTable) 468 fAngleBank.push_back(fAngleForEnergyTable); 465 } 469 } 466 timer.Stop(); 470 timer.Stop(); 467 G4cout.precision(6); 471 G4cout.precision(6); 468 if(verboseLevel > 0) 472 if(verboseLevel > 0) 469 { 473 { 470 G4cout << G4endl; 474 G4cout << G4endl; 471 G4cout << "total time for build X-ray TR a 475 G4cout << "total time for build X-ray TR angle for energy loss tables = " 472 << timer.GetUserElapsed() << " s" < 476 << timer.GetUserElapsed() << " s" << G4endl; 473 } 477 } 474 fGamma = 0.; 478 fGamma = 0.; 475 delete energyVector; 479 delete energyVector; 476 } 480 } 477 481 478 ////////////////////////////////////////////// 482 //////////////////////////////////////////////////////////////////////// 479 // Build XTR angular distribution at given ene 483 // Build XTR angular distribution at given energy based on the model 480 // of transparent regular radiator 484 // of transparent regular radiator 481 void G4VXTRenergyLoss::BuildAngleTable() 485 void G4VXTRenergyLoss::BuildAngleTable() 482 { 486 { 483 G4int iTkin, iTR; 487 G4int iTkin, iTR; 484 G4double energy; 488 G4double energy; 485 489 486 fGammaTkinCut = 0.0; 490 fGammaTkinCut = 0.0; 487 491 488 // setting of min/max TR energies 492 // setting of min/max TR energies 489 if(fGammaTkinCut > fTheMinEnergyTR) 493 if(fGammaTkinCut > fTheMinEnergyTR) 490 fMinEnergyTR = fGammaTkinCut; 494 fMinEnergyTR = fGammaTkinCut; 491 else 495 else 492 fMinEnergyTR = fTheMinEnergyTR; 496 fMinEnergyTR = fTheMinEnergyTR; 493 497 494 if(fGammaTkinCut > fTheMaxEnergyTR) 498 if(fGammaTkinCut > fTheMaxEnergyTR) 495 fMaxEnergyTR = 2.0 * fGammaTkinCut; 499 fMaxEnergyTR = 2.0 * fGammaTkinCut; 496 else 500 else 497 fMaxEnergyTR = fTheMaxEnergyTR; 501 fMaxEnergyTR = fTheMaxEnergyTR; 498 502 499 G4cout.precision(4); 503 G4cout.precision(4); 500 G4Timer timer; 504 G4Timer timer; 501 timer.Start(); 505 timer.Start(); 502 if(verboseLevel > 0) 506 if(verboseLevel > 0) 503 { 507 { 504 G4cout << G4endl << "Lorentz Factor" << "\ 508 G4cout << G4endl << "Lorentz Factor" << "\t" 505 << "XTR photon number" << G4endl << 509 << "XTR photon number" << G4endl << G4endl; 506 } 510 } 507 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // 511 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // Lorentz factor loop 508 { 512 { 509 fGamma = 513 fGamma = 510 1.0 + (fProtonEnergyVector->GetLowEdgeEn 514 1.0 + (fProtonEnergyVector->GetLowEdgeEnergy(iTkin) / proton_mass_c2); 511 515 512 // fMaxThetaTR = 25. * 2500.0 / (fGamma * 516 // fMaxThetaTR = 25. * 2500.0 / (fGamma * fGamma); // theta^2 513 517 514 if(fMaxThetaTR > fTheMaxAngle) 518 if(fMaxThetaTR > fTheMaxAngle) 515 fMaxThetaTR = fTheMaxAngle; 519 fMaxThetaTR = fTheMaxAngle; 516 else 520 else 517 { 521 { 518 if(fMaxThetaTR < fTheMinAngle) 522 if(fMaxThetaTR < fTheMinAngle) 519 fMaxThetaTR = fTheMinAngle; 523 fMaxThetaTR = fTheMinAngle; 520 } 524 } 521 525 522 fAngleForEnergyTable = new G4PhysicsTable( 526 fAngleForEnergyTable = new G4PhysicsTable(fBinTR); 523 527 524 for(iTR = 0; iTR < fBinTR; ++iTR) 528 for(iTR = 0; iTR < fBinTR; ++iTR) 525 { 529 { 526 energy = fXTREnergyVector->GetLowEdgeEne 530 energy = fXTREnergyVector->GetLowEdgeEnergy(iTR); 527 531 528 G4PhysicsFreeVector* angleVector = GetAn 532 G4PhysicsFreeVector* angleVector = GetAngleVector(energy, fBinTR); 529 533 530 fAngleForEnergyTable->insertAt(iTR, angl 534 fAngleForEnergyTable->insertAt(iTR, angleVector); 531 } 535 } 532 fAngleBank.push_back(fAngleForEnergyTable) 536 fAngleBank.push_back(fAngleForEnergyTable); 533 } 537 } 534 timer.Stop(); 538 timer.Stop(); 535 G4cout.precision(6); 539 G4cout.precision(6); 536 if(verboseLevel > 0) 540 if(verboseLevel > 0) 537 { 541 { 538 G4cout << G4endl; 542 G4cout << G4endl; 539 G4cout << "total time for build XTR angle 543 G4cout << "total time for build XTR angle for given energy tables = " 540 << timer.GetUserElapsed() << " s" < 544 << timer.GetUserElapsed() << " s" << G4endl; 541 } 545 } 542 fGamma = 0.; 546 fGamma = 0.; 543 547 544 return; 548 return; 545 } 549 } 546 550 547 ////////////////////////////////////////////// 551 ///////////////////////////////////////////////////////////////////////// 548 // Vector of angles and angle integral distrib 552 // Vector of angles and angle integral distributions 549 G4PhysicsFreeVector* G4VXTRenergyLoss::GetAngl 553 G4PhysicsFreeVector* G4VXTRenergyLoss::GetAngleVector(G4double energy, G4int n) 550 { 554 { 551 G4double theta = 0., result, tmp = 0., cof1, 555 G4double theta = 0., result, tmp = 0., cof1, cof2, cofMin, cofPHC, 552 angleSum = 0.; 556 angleSum = 0.; 553 G4int iTheta, k, kMin; 557 G4int iTheta, k, kMin; 554 558 555 auto angleVector = new G4PhysicsFreeVector(n << 559 G4PhysicsFreeVector* angleVector = new G4PhysicsFreeVector(n); 556 560 557 cofPHC = 4. * pi * hbarc; 561 cofPHC = 4. * pi * hbarc; 558 tmp = (fSigma1 - fSigma2) / cofPHC / ener 562 tmp = (fSigma1 - fSigma2) / cofPHC / energy; 559 cof1 = fPlateThick * tmp; 563 cof1 = fPlateThick * tmp; 560 cof2 = fGasThick * tmp; 564 cof2 = fGasThick * tmp; 561 565 562 cofMin = energy * (fPlateThick + fGasThick) 566 cofMin = energy * (fPlateThick + fGasThick) / fGamma / fGamma; 563 cofMin += (fPlateThick * fSigma1 + fGasThick 567 cofMin += (fPlateThick * fSigma1 + fGasThick * fSigma2) / energy; 564 cofMin /= cofPHC; 568 cofMin /= cofPHC; 565 569 566 kMin = G4int(cofMin); 570 kMin = G4int(cofMin); 567 if(cofMin > kMin) 571 if(cofMin > kMin) 568 kMin++; 572 kMin++; 569 573 570 if(verboseLevel > 2) 574 if(verboseLevel > 2) 571 { 575 { 572 G4cout << "n-1 = " << n - 1 576 G4cout << "n-1 = " << n - 1 573 << "; theta = " << std::sqrt(fMaxTh 577 << "; theta = " << std::sqrt(fMaxThetaTR) * fGamma 574 << "; tmp = " << 0. << "; angleS 578 << "; tmp = " << 0. << "; angleSum = " << angleSum << G4endl; 575 } 579 } 576 580 577 for(iTheta = n - 1; iTheta >= 1; --iTheta) 581 for(iTheta = n - 1; iTheta >= 1; --iTheta) 578 { 582 { 579 k = iTheta - 1 + kMin; 583 k = iTheta - 1 + kMin; 580 tmp = pi * fPlateThick * (k + cof2) / ( 584 tmp = pi * fPlateThick * (k + cof2) / (fPlateThick + fGasThick); 581 result = (k - cof1) * (k - cof1) * (k + co 585 result = (k - cof1) * (k - cof1) * (k + cof2) * (k + cof2); 582 tmp = std::sin(tmp) * std::sin(tmp) * s 586 tmp = std::sin(tmp) * std::sin(tmp) * std::abs(k - cofMin) / result; 583 587 584 if(k == kMin && kMin == G4int(cofMin)) 588 if(k == kMin && kMin == G4int(cofMin)) 585 { 589 { 586 // angleSum += 0.5 * tmp; 590 // angleSum += 0.5 * tmp; 587 angleSum += tmp; // ATLAS TB 591 angleSum += tmp; // ATLAS TB 588 } 592 } 589 else if(iTheta == n - 1) 593 else if(iTheta == n - 1) 590 ; 594 ; 591 else 595 else 592 { 596 { 593 angleSum += tmp; 597 angleSum += tmp; 594 } 598 } 595 theta = std::abs(k - cofMin) * cofPHC / en 599 theta = std::abs(k - cofMin) * cofPHC / energy / (fPlateThick + fGasThick); 596 600 597 if(verboseLevel > 2) 601 if(verboseLevel > 2) 598 { 602 { 599 G4cout << "iTheta = " << iTheta << "; k 603 G4cout << "iTheta = " << iTheta << "; k = " << k 600 << "; theta = " << std::sqrt(thet 604 << "; theta = " << std::sqrt(theta) * fGamma << "; tmp = " << tmp 601 << "; angleSum = " << angleSum 605 << "; angleSum = " << angleSum << G4endl; 602 } 606 } 603 angleVector->PutValue(iTheta, theta, angle 607 angleVector->PutValue(iTheta, theta, angleSum); 604 } 608 } 605 if(theta > 0.) 609 if(theta > 0.) 606 { 610 { 607 // angleSum += 0.5 * tmp; 611 // angleSum += 0.5 * tmp; 608 angleSum += 0.; // ATLAS TB 612 angleSum += 0.; // ATLAS TB 609 theta = 0.; 613 theta = 0.; 610 } 614 } 611 if(verboseLevel > 2) 615 if(verboseLevel > 2) 612 { 616 { 613 G4cout << "iTheta = " << iTheta << "; thet 617 G4cout << "iTheta = " << iTheta << "; theta = " << std::sqrt(theta) * fGamma 614 << "; tmp = " << tmp << "; angle 618 << "; tmp = " << tmp << "; angleSum = " << angleSum << G4endl; 615 } 619 } 616 angleVector->PutValue(iTheta, theta, angleSu 620 angleVector->PutValue(iTheta, theta, angleSum); 617 621 618 return angleVector; 622 return angleVector; 619 } 623 } 620 624 621 ////////////////////////////////////////////// 625 //////////////////////////////////////////////////////////////////////// 622 // Build XTR angular distribution based on the 626 // Build XTR angular distribution based on the model of transparent regular 623 // radiator 627 // radiator 624 void G4VXTRenergyLoss::BuildGlobalAngleTable() 628 void G4VXTRenergyLoss::BuildGlobalAngleTable() 625 { 629 { 626 G4int iTkin, iTR, iPlace; 630 G4int iTkin, iTR, iPlace; 627 G4double radiatorCof = 1.0; // for tuning o 631 G4double radiatorCof = 1.0; // for tuning of XTR yield 628 G4double angleSum; 632 G4double angleSum; 629 fAngleDistrTable = new G4PhysicsTable(fTotBi 633 fAngleDistrTable = new G4PhysicsTable(fTotBin); 630 634 631 fGammaTkinCut = 0.0; 635 fGammaTkinCut = 0.0; 632 636 633 // setting of min/max TR energies 637 // setting of min/max TR energies 634 if(fGammaTkinCut > fTheMinEnergyTR) 638 if(fGammaTkinCut > fTheMinEnergyTR) 635 fMinEnergyTR = fGammaTkinCut; 639 fMinEnergyTR = fGammaTkinCut; 636 else 640 else 637 fMinEnergyTR = fTheMinEnergyTR; 641 fMinEnergyTR = fTheMinEnergyTR; 638 642 639 if(fGammaTkinCut > fTheMaxEnergyTR) 643 if(fGammaTkinCut > fTheMaxEnergyTR) 640 fMaxEnergyTR = 2.0 * fGammaTkinCut; 644 fMaxEnergyTR = 2.0 * fGammaTkinCut; 641 else 645 else 642 fMaxEnergyTR = fTheMaxEnergyTR; 646 fMaxEnergyTR = fTheMaxEnergyTR; 643 647 644 G4cout.precision(4); 648 G4cout.precision(4); 645 G4Timer timer; 649 G4Timer timer; 646 timer.Start(); 650 timer.Start(); 647 if(verboseLevel > 0) 651 if(verboseLevel > 0) 648 { 652 { 649 G4cout << G4endl; 653 G4cout << G4endl; 650 G4cout << "Lorentz Factor" 654 G4cout << "Lorentz Factor" 651 << "\t" 655 << "\t" 652 << "XTR photon number" << G4endl; 656 << "XTR photon number" << G4endl; 653 G4cout << G4endl; 657 G4cout << G4endl; 654 } 658 } 655 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // 659 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // Lorentz factor loop 656 { 660 { 657 fGamma = 661 fGamma = 658 1.0 + (fProtonEnergyVector->GetLowEdgeEn 662 1.0 + (fProtonEnergyVector->GetLowEdgeEnergy(iTkin) / proton_mass_c2); 659 663 660 // fMaxThetaTR = 25.0 / (fGamma * fGamma); 664 // fMaxThetaTR = 25.0 / (fGamma * fGamma); // theta^2 661 // fMaxThetaTR = 1.e-4; // theta^2 665 // fMaxThetaTR = 1.e-4; // theta^2 662 666 663 if(fMaxThetaTR > fTheMaxAngle) 667 if(fMaxThetaTR > fTheMaxAngle) 664 fMaxThetaTR = fTheMaxAngle; 668 fMaxThetaTR = fTheMaxAngle; 665 else 669 else 666 { 670 { 667 if(fMaxThetaTR < fTheMinAngle) 671 if(fMaxThetaTR < fTheMinAngle) 668 fMaxThetaTR = fTheMinAngle; 672 fMaxThetaTR = fTheMinAngle; 669 } 673 } 670 auto angleVector = << 674 G4PhysicsLinearVector* angleVector = 671 // G4PhysicsLogVector* angleVector = 675 // G4PhysicsLogVector* angleVector = 672 new G4PhysicsLinearVector(0.0, fMaxTheta 676 new G4PhysicsLinearVector(0.0, fMaxThetaTR, fBinTR); 673 // new G4PhysicsLogVector(1.e-8, fMaxThet 677 // new G4PhysicsLogVector(1.e-8, fMaxThetaTR, fBinTR); 674 678 675 angleSum = 0.0; 679 angleSum = 0.0; 676 680 677 G4Integrator<G4VXTRenergyLoss, G4double (G 681 G4Integrator<G4VXTRenergyLoss, G4double (G4VXTRenergyLoss::*)(G4double)> 678 integral; 682 integral; 679 683 680 angleVector->PutValue(fBinTR - 1, angleSum 684 angleVector->PutValue(fBinTR - 1, angleSum); 681 685 682 for(iTR = fBinTR - 2; iTR >= 0; --iTR) 686 for(iTR = fBinTR - 2; iTR >= 0; --iTR) 683 { 687 { 684 angleSum += radiatorCof * fCofTR * 688 angleSum += radiatorCof * fCofTR * 685 integral.Legendre96(this, &G 689 integral.Legendre96(this, &G4VXTRenergyLoss::AngleXTRdEdx, 686 angleVec 690 angleVector->GetLowEdgeEnergy(iTR), 687 angleVec 691 angleVector->GetLowEdgeEnergy(iTR + 1)); 688 692 689 angleVector->PutValue(iTR, angleSum); 693 angleVector->PutValue(iTR, angleSum); 690 } 694 } 691 if(verboseLevel > 1) 695 if(verboseLevel > 1) 692 { 696 { 693 G4cout << fGamma << "\t" << angleSum << 697 G4cout << fGamma << "\t" << angleSum << G4endl; 694 } 698 } 695 iPlace = iTkin; 699 iPlace = iTkin; 696 fAngleDistrTable->insertAt(iPlace, angleVe 700 fAngleDistrTable->insertAt(iPlace, angleVector); 697 } 701 } 698 timer.Stop(); 702 timer.Stop(); 699 G4cout.precision(6); 703 G4cout.precision(6); 700 if(verboseLevel > 0) 704 if(verboseLevel > 0) 701 { 705 { 702 G4cout << G4endl; 706 G4cout << G4endl; 703 G4cout << "total time for build X-ray TR a 707 G4cout << "total time for build X-ray TR angle tables = " 704 << timer.GetUserElapsed() << " s" < 708 << timer.GetUserElapsed() << " s" << G4endl; 705 } 709 } 706 fGamma = 0.; 710 fGamma = 0.; 707 711 708 return; 712 return; 709 } 713 } 710 714 711 ////////////////////////////////////////////// 715 ////////////////////////////////////////////////////////////////////////////// 712 // The main function which is responsible for 716 // The main function which is responsible for the treatment of a particle 713 // passage through G4Envelope with discrete ge 717 // passage through G4Envelope with discrete generation of G4Gamma 714 G4VParticleChange* G4VXTRenergyLoss::PostStepD 718 G4VParticleChange* G4VXTRenergyLoss::PostStepDoIt(const G4Track& aTrack, 715 719 const G4Step& aStep) 716 { 720 { 717 G4int iTkin; 721 G4int iTkin; 718 G4double energyTR, theta, theta2, phi, dirX, 722 G4double energyTR, theta, theta2, phi, dirX, dirY, dirZ; 719 723 720 fParticleChange.Initialize(aTrack); 724 fParticleChange.Initialize(aTrack); 721 725 722 if(verboseLevel > 1) 726 if(verboseLevel > 1) 723 { 727 { 724 G4cout << "Start of G4VXTRenergyLoss::Post 728 G4cout << "Start of G4VXTRenergyLoss::PostStepDoIt " << G4endl; 725 G4cout << "name of current material = " 729 G4cout << "name of current material = " 726 << aTrack.GetVolume()->GetLogicalVo 730 << aTrack.GetVolume()->GetLogicalVolume()->GetMaterial()->GetName() 727 << G4endl; 731 << G4endl; 728 } 732 } 729 if(aTrack.GetVolume()->GetLogicalVolume() != 733 if(aTrack.GetVolume()->GetLogicalVolume() != fEnvelope) 730 { 734 { 731 if(verboseLevel > 0) 735 if(verboseLevel > 0) 732 { 736 { 733 G4cout << "Go out from G4VXTRenergyLoss: 737 G4cout << "Go out from G4VXTRenergyLoss::PostStepDoIt: wrong volume " 734 << G4endl; 738 << G4endl; 735 } 739 } 736 return G4VDiscreteProcess::PostStepDoIt(aT 740 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 737 } 741 } 738 else 742 else 739 { 743 { 740 G4StepPoint* pPostStepPoint = aStep 744 G4StepPoint* pPostStepPoint = aStep.GetPostStepPoint(); 741 const G4DynamicParticle* aParticle = aTrac 745 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 742 746 743 // Now we are ready to Generate one TR pho 747 // Now we are ready to Generate one TR photon 744 G4double kinEnergy = aParticle->GetKinetic 748 G4double kinEnergy = aParticle->GetKineticEnergy(); 745 G4double mass = aParticle->GetDefinit 749 G4double mass = aParticle->GetDefinition()->GetPDGMass(); 746 G4double gamma = 1.0 + kinEnergy / mas 750 G4double gamma = 1.0 + kinEnergy / mass; 747 751 748 if(verboseLevel > 1) 752 if(verboseLevel > 1) 749 { 753 { 750 G4cout << "gamma = " << gamma << G4endl; 754 G4cout << "gamma = " << gamma << G4endl; 751 } 755 } 752 G4double massRatio = proton_mass 756 G4double massRatio = proton_mass_c2 / mass; 753 G4double TkinScaled = kinEnergy * 757 G4double TkinScaled = kinEnergy * massRatio; 754 G4ThreeVector position = pPostStepPo 758 G4ThreeVector position = pPostStepPoint->GetPosition(); 755 G4ParticleMomentum direction = aParticle-> 759 G4ParticleMomentum direction = aParticle->GetMomentumDirection(); 756 G4double startTime = pPostStepPo 760 G4double startTime = pPostStepPoint->GetGlobalTime(); 757 761 758 for(iTkin = 0; iTkin < fTotBin; ++iTkin) 762 for(iTkin = 0; iTkin < fTotBin; ++iTkin) 759 { 763 { 760 if(TkinScaled < fProtonEnergyVector->Get 764 if(TkinScaled < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) 761 break; 765 break; 762 } 766 } 763 767 764 if(iTkin == 0) // Tkin is too small, negl 768 if(iTkin == 0) // Tkin is too small, neglect of TR photon generation 765 { 769 { 766 if(verboseLevel > 0) 770 if(verboseLevel > 0) 767 { 771 { 768 G4cout << "Go out from G4VXTRenergyLos 772 G4cout << "Go out from G4VXTRenergyLoss::PostStepDoIt:iTkin = " << iTkin 769 << G4endl; 773 << G4endl; 770 } 774 } 771 return G4VDiscreteProcess::PostStepDoIt( 775 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 772 } 776 } 773 else // general case: Tkin between two ve 777 else // general case: Tkin between two vectors of the material 774 { 778 { 775 fParticleChange.SetNumberOfSecondaries(1 779 fParticleChange.SetNumberOfSecondaries(1); 776 780 777 energyTR = GetXTRrandomEnergy(TkinScaled 781 energyTR = GetXTRrandomEnergy(TkinScaled, iTkin); 778 782 779 if(verboseLevel > 1) 783 if(verboseLevel > 1) 780 { 784 { 781 G4cout << "energyTR = " << energyTR / 785 G4cout << "energyTR = " << energyTR / keV << " keV" << G4endl; 782 } 786 } 783 if(fAngleRadDistr) 787 if(fAngleRadDistr) 784 { 788 { 785 theta2 = GetRandomAngle(energyTR, iTki 789 theta2 = GetRandomAngle(energyTR, iTkin); 786 if(theta2 > 0.) 790 if(theta2 > 0.) 787 theta = std::sqrt(theta2); 791 theta = std::sqrt(theta2); 788 else 792 else 789 theta = 0.; 793 theta = 0.; 790 } 794 } 791 else 795 else 792 theta = std::fabs(G4RandGauss::shoot(0 796 theta = std::fabs(G4RandGauss::shoot(0.0, pi / gamma)); 793 797 794 if(theta >= 0.1) 798 if(theta >= 0.1) 795 theta = 0.1; 799 theta = 0.1; 796 800 797 phi = twopi * G4UniformRand(); 801 phi = twopi * G4UniformRand(); 798 802 799 dirX = std::sin(theta) * std::cos(phi); 803 dirX = std::sin(theta) * std::cos(phi); 800 dirY = std::sin(theta) * std::sin(phi); 804 dirY = std::sin(theta) * std::sin(phi); 801 dirZ = std::cos(theta); 805 dirZ = std::cos(theta); 802 806 803 G4ThreeVector directionTR(dirX, dirY, di 807 G4ThreeVector directionTR(dirX, dirY, dirZ); 804 directionTR.rotateUz(direction); 808 directionTR.rotateUz(direction); 805 directionTR.unit(); 809 directionTR.unit(); 806 810 807 auto aPhotonTR = << 811 G4DynamicParticle* aPhotonTR = 808 new G4DynamicParticle(G4Gamma::Gamma() 812 new G4DynamicParticle(G4Gamma::Gamma(), directionTR, energyTR); 809 813 810 // A XTR photon is set on the particle t 814 // A XTR photon is set on the particle track inside the radiator 811 // and is moved to the G4Envelope surfac 815 // and is moved to the G4Envelope surface for standard X-ray TR models 812 // only. The case of fExitFlux=true 816 // only. The case of fExitFlux=true 813 817 814 if(fExitFlux) 818 if(fExitFlux) 815 { 819 { 816 const G4RotationMatrix* rotM = 820 const G4RotationMatrix* rotM = 817 pPostStepPoint->GetTouchable()->GetR 821 pPostStepPoint->GetTouchable()->GetRotation(); 818 G4ThreeVector transl = pPostStepPoint- 822 G4ThreeVector transl = pPostStepPoint->GetTouchable()->GetTranslation(); 819 G4AffineTransform transform = G4Affine 823 G4AffineTransform transform = G4AffineTransform(rotM, transl); 820 transform.Invert(); 824 transform.Invert(); 821 G4ThreeVector localP = transform.Trans 825 G4ThreeVector localP = transform.TransformPoint(position); 822 G4ThreeVector localV = transform.Trans 826 G4ThreeVector localV = transform.TransformAxis(directionTR); 823 827 824 G4double distance = 828 G4double distance = 825 fEnvelope->GetSolid()->DistanceToOut 829 fEnvelope->GetSolid()->DistanceToOut(localP, localV); 826 if(verboseLevel > 1) 830 if(verboseLevel > 1) 827 { 831 { 828 G4cout << "distance to exit = " << d 832 G4cout << "distance to exit = " << distance / mm << " mm" << G4endl; 829 } 833 } 830 position += distance * directionTR; 834 position += distance * directionTR; 831 startTime += distance / c_light; 835 startTime += distance / c_light; 832 } 836 } 833 G4Track* aSecondaryTrack = new G4Track(a 837 G4Track* aSecondaryTrack = new G4Track(aPhotonTR, startTime, position); 834 aSecondaryTrack->SetTouchableHandle( 838 aSecondaryTrack->SetTouchableHandle( 835 aStep.GetPostStepPoint()->GetTouchable 839 aStep.GetPostStepPoint()->GetTouchableHandle()); 836 aSecondaryTrack->SetParentID(aTrack.GetT 840 aSecondaryTrack->SetParentID(aTrack.GetTrackID()); 837 841 838 fParticleChange.AddSecondary(aSecondaryT 842 fParticleChange.AddSecondary(aSecondaryTrack); 839 fParticleChange.ProposeEnergy(kinEnergy) 843 fParticleChange.ProposeEnergy(kinEnergy); 840 } 844 } 841 } 845 } 842 return G4VDiscreteProcess::PostStepDoIt(aTra 846 return G4VDiscreteProcess::PostStepDoIt(aTrack, aStep); 843 } 847 } 844 848 845 ////////////////////////////////////////////// 849 /////////////////////////////////////////////////////////////////////// 846 // This function returns the spectral and angl 850 // This function returns the spectral and angle density of TR quanta 847 // in X-ray energy region generated forward wh 851 // in X-ray energy region generated forward when a relativistic 848 // charged particle crosses interface between 852 // charged particle crosses interface between two materials. 849 // The high energy small theta approximation i 853 // The high energy small theta approximation is applied. 850 // (matter1 -> matter2, or 2->1) 854 // (matter1 -> matter2, or 2->1) 851 // varAngle =2* (1 - std::cos(theta)) or appro 855 // varAngle =2* (1 - std::cos(theta)) or approximately = theta*theta 852 G4complex G4VXTRenergyLoss::OneInterfaceXTRdEd 856 G4complex G4VXTRenergyLoss::OneInterfaceXTRdEdx(G4double energy, G4double gamma, 853 857 G4double varAngle) 854 { 858 { 855 G4complex Z1 = GetPlateComplexFZ(energy, gam 859 G4complex Z1 = GetPlateComplexFZ(energy, gamma, varAngle); 856 G4complex Z2 = GetGasComplexFZ(energy, gamma 860 G4complex Z2 = GetGasComplexFZ(energy, gamma, varAngle); 857 861 858 G4complex zOut = (Z1 - Z2) * (Z1 - Z2) * (va 862 G4complex zOut = (Z1 - Z2) * (Z1 - Z2) * (varAngle * energy / hbarc / hbarc); 859 return zOut; 863 return zOut; 860 } 864 } 861 865 862 ////////////////////////////////////////////// 866 ////////////////////////////////////////////////////////////////////////////// 863 // For photon energy distribution tables. Inte 867 // For photon energy distribution tables. Integrate first over angle 864 G4double G4VXTRenergyLoss::SpectralAngleXTRdEd 868 G4double G4VXTRenergyLoss::SpectralAngleXTRdEdx(G4double varAngle) 865 { 869 { 866 G4double result = GetStackFactor(fEnergy, fG 870 G4double result = GetStackFactor(fEnergy, fGamma, varAngle); 867 if(result < 0.0) 871 if(result < 0.0) 868 result = 0.0; 872 result = 0.0; 869 return result; 873 return result; 870 } 874 } 871 875 872 ////////////////////////////////////////////// 876 ///////////////////////////////////////////////////////////////////////// 873 // For second integration over energy 877 // For second integration over energy 874 G4double G4VXTRenergyLoss::SpectralXTRdEdx(G4d 878 G4double G4VXTRenergyLoss::SpectralXTRdEdx(G4double energy) 875 { 879 { 876 G4int i; 880 G4int i; 877 static constexpr G4int iMax = 8; 881 static constexpr G4int iMax = 8; 878 G4double angleSum = 0.0; 882 G4double angleSum = 0.0; 879 883 880 G4double lim[iMax] = { 0.0, 0.01, 0.02, 0.05 884 G4double lim[iMax] = { 0.0, 0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0 }; 881 885 882 for(i = 0; i < iMax; ++i) 886 for(i = 0; i < iMax; ++i) 883 lim[i] *= fMaxThetaTR; 887 lim[i] *= fMaxThetaTR; 884 888 885 G4Integrator<G4VXTRenergyLoss, G4double (G4V 889 G4Integrator<G4VXTRenergyLoss, G4double (G4VXTRenergyLoss::*)(G4double)> 886 integral; 890 integral; 887 891 888 fEnergy = energy; 892 fEnergy = energy; 889 { 893 { 890 for(i = 0; i < iMax - 1; ++i) 894 for(i = 0; i < iMax - 1; ++i) 891 { 895 { 892 angleSum += integral.Legendre96( 896 angleSum += integral.Legendre96( 893 this, &G4VXTRenergyLoss::SpectralAngle 897 this, &G4VXTRenergyLoss::SpectralAngleXTRdEdx, lim[i], lim[i + 1]); 894 } 898 } 895 } 899 } 896 return angleSum; 900 return angleSum; 897 } 901 } 898 902 899 ////////////////////////////////////////////// 903 ////////////////////////////////////////////////////////////////////////// 900 // for photon angle distribution tables 904 // for photon angle distribution tables 901 G4double G4VXTRenergyLoss::AngleSpectralXTRdEd 905 G4double G4VXTRenergyLoss::AngleSpectralXTRdEdx(G4double energy) 902 { 906 { 903 G4double result = GetStackFactor(energy, fGa 907 G4double result = GetStackFactor(energy, fGamma, fVarAngle); 904 if(result < 0) 908 if(result < 0) 905 result = 0.0; 909 result = 0.0; 906 return result; 910 return result; 907 } 911 } 908 912 909 ////////////////////////////////////////////// 913 /////////////////////////////////////////////////////////////////////////// 910 // The XTR angular distribution based on trans 914 // The XTR angular distribution based on transparent regular radiator 911 G4double G4VXTRenergyLoss::AngleXTRdEdx(G4doub 915 G4double G4VXTRenergyLoss::AngleXTRdEdx(G4double varAngle) 912 { 916 { 913 G4double result; 917 G4double result; 914 G4double sum = 0., tmp1, tmp2, tmp = 0., cof 918 G4double sum = 0., tmp1, tmp2, tmp = 0., cof1, cof2, cofMin, cofPHC, energy1, 915 energy2; 919 energy2; 916 G4int k, kMax, kMin, i; 920 G4int k, kMax, kMin, i; 917 921 918 cofPHC = twopi * hbarc; 922 cofPHC = twopi * hbarc; 919 923 920 cof1 = (fPlateThick + fGasThick) * (1. / fGa 924 cof1 = (fPlateThick + fGasThick) * (1. / fGamma / fGamma + varAngle); 921 cof2 = fPlateThick * fSigma1 + fGasThick * f 925 cof2 = fPlateThick * fSigma1 + fGasThick * fSigma2; 922 926 923 cofMin = std::sqrt(cof1 * cof2); 927 cofMin = std::sqrt(cof1 * cof2); 924 cofMin /= cofPHC; 928 cofMin /= cofPHC; 925 929 926 kMin = G4int(cofMin); 930 kMin = G4int(cofMin); 927 if(cofMin > kMin) 931 if(cofMin > kMin) 928 kMin++; 932 kMin++; 929 933 930 kMax = kMin + 9; 934 kMax = kMin + 9; 931 935 932 for(k = kMin; k <= kMax; ++k) 936 for(k = kMin; k <= kMax; ++k) 933 { 937 { 934 tmp1 = cofPHC * k; 938 tmp1 = cofPHC * k; 935 tmp2 = std::sqrt(tmp1 * tmp1 - cof1 * c 939 tmp2 = std::sqrt(tmp1 * tmp1 - cof1 * cof2); 936 energy1 = (tmp1 + tmp2) / cof1; 940 energy1 = (tmp1 + tmp2) / cof1; 937 energy2 = (tmp1 - tmp2) / cof1; 941 energy2 = (tmp1 - tmp2) / cof1; 938 942 939 for(i = 0; i < 2; ++i) 943 for(i = 0; i < 2; ++i) 940 { 944 { 941 if(i == 0) 945 if(i == 0) 942 { 946 { 943 if(energy1 > fTheMaxEnergyTR || energy 947 if(energy1 > fTheMaxEnergyTR || energy1 < fTheMinEnergyTR) 944 continue; 948 continue; 945 949 946 tmp1 = 950 tmp1 = 947 (energy1 * energy1 * (1. / fGamma / 951 (energy1 * energy1 * (1. / fGamma / fGamma + varAngle) + fSigma1) * 948 fPlateThick / (4 * hbarc * energy1); 952 fPlateThick / (4 * hbarc * energy1); 949 tmp2 = std::sin(tmp1); 953 tmp2 = std::sin(tmp1); 950 tmp = energy1 * tmp2 * tmp2; 954 tmp = energy1 * tmp2 * tmp2; 951 tmp2 = fPlateThick / (4. * tmp1); 955 tmp2 = fPlateThick / (4. * tmp1); 952 tmp1 = 956 tmp1 = 953 hbarc * energy1 / 957 hbarc * energy1 / 954 (energy1 * energy1 * (1. / fGamma / 958 (energy1 * energy1 * (1. / fGamma / fGamma + varAngle) + fSigma2); 955 tmp *= (tmp1 - tmp2) * (tmp1 - tmp2); 959 tmp *= (tmp1 - tmp2) * (tmp1 - tmp2); 956 tmp1 = cof1 / (4. * hbarc) - cof2 / (4 960 tmp1 = cof1 / (4. * hbarc) - cof2 / (4. * hbarc * energy1 * energy1); 957 tmp2 = std::abs(tmp1); 961 tmp2 = std::abs(tmp1); 958 962 959 if(tmp2 > 0.) 963 if(tmp2 > 0.) 960 tmp /= tmp2; 964 tmp /= tmp2; 961 else 965 else 962 continue; 966 continue; 963 } 967 } 964 else 968 else 965 { 969 { 966 if(energy2 > fTheMaxEnergyTR || energy 970 if(energy2 > fTheMaxEnergyTR || energy2 < fTheMinEnergyTR) 967 continue; 971 continue; 968 972 969 tmp1 = 973 tmp1 = 970 (energy2 * energy2 * (1. / fGamma / 974 (energy2 * energy2 * (1. / fGamma / fGamma + varAngle) + fSigma1) * 971 fPlateThick / (4. * hbarc * energy2) 975 fPlateThick / (4. * hbarc * energy2); 972 tmp2 = std::sin(tmp1); 976 tmp2 = std::sin(tmp1); 973 tmp = energy2 * tmp2 * tmp2; 977 tmp = energy2 * tmp2 * tmp2; 974 tmp2 = fPlateThick / (4. * tmp1); 978 tmp2 = fPlateThick / (4. * tmp1); 975 tmp1 = 979 tmp1 = 976 hbarc * energy2 / 980 hbarc * energy2 / 977 (energy2 * energy2 * (1. / fGamma / 981 (energy2 * energy2 * (1. / fGamma / fGamma + varAngle) + fSigma2); 978 tmp *= (tmp1 - tmp2) * (tmp1 - tmp2); 982 tmp *= (tmp1 - tmp2) * (tmp1 - tmp2); 979 tmp1 = cof1 / (4. * hbarc) - cof2 / (4 983 tmp1 = cof1 / (4. * hbarc) - cof2 / (4. * hbarc * energy2 * energy2); 980 tmp2 = std::abs(tmp1); 984 tmp2 = std::abs(tmp1); 981 985 982 if(tmp2 > 0.) 986 if(tmp2 > 0.) 983 tmp /= tmp2; 987 tmp /= tmp2; 984 else 988 else 985 continue; 989 continue; 986 } 990 } 987 sum += tmp; 991 sum += tmp; 988 } 992 } 989 } 993 } 990 result = 4. * pi * fPlateNumber * sum * varA 994 result = 4. * pi * fPlateNumber * sum * varAngle; 991 result /= hbarc * hbarc; 995 result /= hbarc * hbarc; 992 996 993 return result; 997 return result; 994 } 998 } 995 999 996 ////////////////////////////////////////////// 1000 ////////////////////////////////////////////////////////////////////// 997 // Calculates formation zone for plates. Omega 1001 // Calculates formation zone for plates. Omega is energy !!! 998 G4double G4VXTRenergyLoss::GetPlateFormationZo 1002 G4double G4VXTRenergyLoss::GetPlateFormationZone(G4double omega, G4double gamma, 999 1003 G4double varAngle) 1000 { 1004 { 1001 G4double cof, lambda; 1005 G4double cof, lambda; 1002 lambda = 1.0 / gamma / gamma + varAngle + f 1006 lambda = 1.0 / gamma / gamma + varAngle + fSigma1 / omega / omega; 1003 cof = 2.0 * hbarc / omega / lambda; 1007 cof = 2.0 * hbarc / omega / lambda; 1004 return cof; 1008 return cof; 1005 } 1009 } 1006 1010 1007 ///////////////////////////////////////////// 1011 ////////////////////////////////////////////////////////////////////// 1008 // Calculates complex formation zone for plat 1012 // Calculates complex formation zone for plates. Omega is energy !!! 1009 G4complex G4VXTRenergyLoss::GetPlateComplexFZ 1013 G4complex G4VXTRenergyLoss::GetPlateComplexFZ(G4double omega, G4double gamma, 1010 1014 G4double varAngle) 1011 { 1015 { 1012 G4double cof, length, delta, real_v, image_ 1016 G4double cof, length, delta, real_v, image_v; 1013 1017 1014 length = 0.5 * GetPlateFormationZone(omega, 1018 length = 0.5 * GetPlateFormationZone(omega, gamma, varAngle); 1015 delta = length * GetPlateLinearPhotoAbs(om 1019 delta = length * GetPlateLinearPhotoAbs(omega); 1016 cof = 1.0 / (1.0 + delta * delta); 1020 cof = 1.0 / (1.0 + delta * delta); 1017 1021 1018 real_v = length * cof; 1022 real_v = length * cof; 1019 image_v = real_v * delta; 1023 image_v = real_v * delta; 1020 1024 1021 G4complex zone(real_v, image_v); 1025 G4complex zone(real_v, image_v); 1022 return zone; 1026 return zone; 1023 } 1027 } 1024 1028 1025 ///////////////////////////////////////////// 1029 //////////////////////////////////////////////////////////////////////// 1026 // Computes matrix of Sandia photo absorption 1030 // Computes matrix of Sandia photo absorption cross section coefficients for 1027 // plate material 1031 // plate material 1028 void G4VXTRenergyLoss::ComputePlatePhotoAbsCo 1032 void G4VXTRenergyLoss::ComputePlatePhotoAbsCof() 1029 { 1033 { 1030 const G4MaterialTable* theMaterialTable = G 1034 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 1031 const G4Material* mat = ( 1035 const G4Material* mat = (*theMaterialTable)[fMatIndex1]; 1032 fPlatePhotoAbsCof = m 1036 fPlatePhotoAbsCof = mat->GetSandiaTable(); 1033 1037 1034 return; 1038 return; 1035 } 1039 } 1036 1040 1037 ///////////////////////////////////////////// 1041 ////////////////////////////////////////////////////////////////////// 1038 // Returns the value of linear photo absorpti 1042 // Returns the value of linear photo absorption coefficient (in reciprocal 1039 // length) for plate for given energy of X-ra 1043 // length) for plate for given energy of X-ray photon omega 1040 G4double G4VXTRenergyLoss::GetPlateLinearPhot 1044 G4double G4VXTRenergyLoss::GetPlateLinearPhotoAbs(G4double omega) 1041 { 1045 { 1042 G4double omega2, omega3, omega4; 1046 G4double omega2, omega3, omega4; 1043 1047 1044 omega2 = omega * omega; 1048 omega2 = omega * omega; 1045 omega3 = omega2 * omega; 1049 omega3 = omega2 * omega; 1046 omega4 = omega2 * omega2; 1050 omega4 = omega2 * omega2; 1047 1051 1048 const G4double* SandiaCof = fPlatePhotoAbsC 1052 const G4double* SandiaCof = fPlatePhotoAbsCof->GetSandiaCofForMaterial(omega); 1049 G4double cross = SandiaCof[0] / 1053 G4double cross = SandiaCof[0] / omega + SandiaCof[1] / omega2 + 1050 SandiaCof[2] / omega3 + Sa 1054 SandiaCof[2] / omega3 + SandiaCof[3] / omega4; 1051 return cross; 1055 return cross; 1052 } 1056 } 1053 1057 1054 ///////////////////////////////////////////// 1058 ////////////////////////////////////////////////////////////////////// 1055 // Calculates formation zone for gas. Omega i 1059 // Calculates formation zone for gas. Omega is energy !!! 1056 G4double G4VXTRenergyLoss::GetGasFormationZon 1060 G4double G4VXTRenergyLoss::GetGasFormationZone(G4double omega, G4double gamma, 1057 1061 G4double varAngle) 1058 { 1062 { 1059 G4double cof, lambda; 1063 G4double cof, lambda; 1060 lambda = 1.0 / gamma / gamma + varAngle + f 1064 lambda = 1.0 / gamma / gamma + varAngle + fSigma2 / omega / omega; 1061 cof = 2.0 * hbarc / omega / lambda; 1065 cof = 2.0 * hbarc / omega / lambda; 1062 return cof; 1066 return cof; 1063 } 1067 } 1064 1068 1065 ///////////////////////////////////////////// 1069 ////////////////////////////////////////////////////////////////////// 1066 // Calculates complex formation zone for gas 1070 // Calculates complex formation zone for gas gaps. Omega is energy !!! 1067 G4complex G4VXTRenergyLoss::GetGasComplexFZ(G 1071 G4complex G4VXTRenergyLoss::GetGasComplexFZ(G4double omega, G4double gamma, 1068 G 1072 G4double varAngle) 1069 { 1073 { 1070 G4double cof, length, delta, real_v, image_ 1074 G4double cof, length, delta, real_v, image_v; 1071 1075 1072 length = 0.5 * GetGasFormationZone(omega, g 1076 length = 0.5 * GetGasFormationZone(omega, gamma, varAngle); 1073 delta = length * GetGasLinearPhotoAbs(omeg 1077 delta = length * GetGasLinearPhotoAbs(omega); 1074 cof = 1.0 / (1.0 + delta * delta); 1078 cof = 1.0 / (1.0 + delta * delta); 1075 1079 1076 real_v = length * cof; 1080 real_v = length * cof; 1077 image_v = real_v * delta; 1081 image_v = real_v * delta; 1078 1082 1079 G4complex zone(real_v, image_v); 1083 G4complex zone(real_v, image_v); 1080 return zone; 1084 return zone; 1081 } 1085 } 1082 1086 1083 ///////////////////////////////////////////// 1087 //////////////////////////////////////////////////////////////////////// 1084 // Computes matrix of Sandia photo absorption 1088 // Computes matrix of Sandia photo absorption cross section coefficients for 1085 // gas material 1089 // gas material 1086 void G4VXTRenergyLoss::ComputeGasPhotoAbsCof( 1090 void G4VXTRenergyLoss::ComputeGasPhotoAbsCof() 1087 { 1091 { 1088 const G4MaterialTable* theMaterialTable = G 1092 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 1089 const G4Material* mat = ( 1093 const G4Material* mat = (*theMaterialTable)[fMatIndex2]; 1090 fGasPhotoAbsCof = m 1094 fGasPhotoAbsCof = mat->GetSandiaTable(); 1091 return; 1095 return; 1092 } 1096 } 1093 1097 1094 ///////////////////////////////////////////// 1098 ////////////////////////////////////////////////////////////////////// 1095 // Returns the value of linear photo absorpti 1099 // Returns the value of linear photo absorption coefficient (in reciprocal 1096 // length) for gas 1100 // length) for gas 1097 G4double G4VXTRenergyLoss::GetGasLinearPhotoA 1101 G4double G4VXTRenergyLoss::GetGasLinearPhotoAbs(G4double omega) 1098 { 1102 { 1099 G4double omega2, omega3, omega4; 1103 G4double omega2, omega3, omega4; 1100 1104 1101 omega2 = omega * omega; 1105 omega2 = omega * omega; 1102 omega3 = omega2 * omega; 1106 omega3 = omega2 * omega; 1103 omega4 = omega2 * omega2; 1107 omega4 = omega2 * omega2; 1104 1108 1105 const G4double* SandiaCof = fGasPhotoAbsCof 1109 const G4double* SandiaCof = fGasPhotoAbsCof->GetSandiaCofForMaterial(omega); 1106 G4double cross = SandiaCof[0] / 1110 G4double cross = SandiaCof[0] / omega + SandiaCof[1] / omega2 + 1107 SandiaCof[2] / omega3 + Sa 1111 SandiaCof[2] / omega3 + SandiaCof[3] / omega4; 1108 return cross; 1112 return cross; 1109 } 1113 } 1110 1114 1111 ///////////////////////////////////////////// 1115 ////////////////////////////////////////////////////////////////////// 1112 // Calculates the product of linear cof by fo 1116 // Calculates the product of linear cof by formation zone for plate. 1113 // Omega is energy !!! 1117 // Omega is energy !!! 1114 G4double G4VXTRenergyLoss::GetPlateZmuProduct 1118 G4double G4VXTRenergyLoss::GetPlateZmuProduct(G4double omega, G4double gamma, 1115 1119 G4double varAngle) 1116 { 1120 { 1117 return GetPlateFormationZone(omega, gamma, 1121 return GetPlateFormationZone(omega, gamma, varAngle) * 1118 GetPlateLinearPhotoAbs(omega); 1122 GetPlateLinearPhotoAbs(omega); 1119 } 1123 } 1120 ///////////////////////////////////////////// 1124 ////////////////////////////////////////////////////////////////////// 1121 // Calculates the product of linear cof by fo 1125 // Calculates the product of linear cof by formation zone for plate. 1122 // G4cout and output in file in some energy r 1126 // G4cout and output in file in some energy range. 1123 void G4VXTRenergyLoss::GetPlateZmuProduct() 1127 void G4VXTRenergyLoss::GetPlateZmuProduct() 1124 { 1128 { 1125 std::ofstream outPlate("plateZmu.dat", std: 1129 std::ofstream outPlate("plateZmu.dat", std::ios::out); 1126 outPlate.setf(std::ios::scientific, std::io 1130 outPlate.setf(std::ios::scientific, std::ios::floatfield); 1127 1131 1128 G4int i; 1132 G4int i; 1129 G4double omega, varAngle, gamma; 1133 G4double omega, varAngle, gamma; 1130 gamma = 10000.; 1134 gamma = 10000.; 1131 varAngle = 1 / gamma / gamma; 1135 varAngle = 1 / gamma / gamma; 1132 if(verboseLevel > 0) 1136 if(verboseLevel > 0) 1133 G4cout << "energy, keV" << "\t" << "Zmu f 1137 G4cout << "energy, keV" << "\t" << "Zmu for plate" << G4endl; 1134 for(i = 0; i < 100; ++i) 1138 for(i = 0; i < 100; ++i) 1135 { 1139 { 1136 omega = (1.0 + i) * keV; 1140 omega = (1.0 + i) * keV; 1137 if(verboseLevel > 1) 1141 if(verboseLevel > 1) 1138 G4cout << omega / keV << "\t" 1142 G4cout << omega / keV << "\t" 1139 << GetPlateZmuProduct(omega, gam 1143 << GetPlateZmuProduct(omega, gamma, varAngle) << "\t"; 1140 if(verboseLevel > 0) 1144 if(verboseLevel > 0) 1141 outPlate << omega / keV << "\t\t" 1145 outPlate << omega / keV << "\t\t" 1142 << GetPlateZmuProduct(omega, g 1146 << GetPlateZmuProduct(omega, gamma, varAngle) << G4endl; 1143 } 1147 } 1144 return; 1148 return; 1145 } 1149 } 1146 1150 1147 ///////////////////////////////////////////// 1151 ////////////////////////////////////////////////////////////////////// 1148 // Calculates the product of linear cof by fo 1152 // Calculates the product of linear cof by formation zone for gas. 1149 // Omega is energy !!! 1153 // Omega is energy !!! 1150 G4double G4VXTRenergyLoss::GetGasZmuProduct(G 1154 G4double G4VXTRenergyLoss::GetGasZmuProduct(G4double omega, G4double gamma, 1151 G 1155 G4double varAngle) 1152 { 1156 { 1153 return GetGasFormationZone(omega, gamma, va 1157 return GetGasFormationZone(omega, gamma, varAngle) * 1154 GetGasLinearPhotoAbs(omega); 1158 GetGasLinearPhotoAbs(omega); 1155 } 1159 } 1156 1160 1157 ///////////////////////////////////////////// 1161 ////////////////////////////////////////////////////////////////////// 1158 // Calculates the product of linear cof by fo 1162 // Calculates the product of linear cof by formation zone for gas. 1159 // G4cout and output in file in some energy r 1163 // G4cout and output in file in some energy range. 1160 void G4VXTRenergyLoss::GetGasZmuProduct() 1164 void G4VXTRenergyLoss::GetGasZmuProduct() 1161 { 1165 { 1162 std::ofstream outGas("gasZmu.dat", std::ios 1166 std::ofstream outGas("gasZmu.dat", std::ios::out); 1163 outGas.setf(std::ios::scientific, std::ios: 1167 outGas.setf(std::ios::scientific, std::ios::floatfield); 1164 G4int i; 1168 G4int i; 1165 G4double omega, varAngle, gamma; 1169 G4double omega, varAngle, gamma; 1166 gamma = 10000.; 1170 gamma = 10000.; 1167 varAngle = 1 / gamma / gamma; 1171 varAngle = 1 / gamma / gamma; 1168 if(verboseLevel > 0) 1172 if(verboseLevel > 0) 1169 G4cout << "energy, keV" << "\t" << "Zmu f 1173 G4cout << "energy, keV" << "\t" << "Zmu for gas" << G4endl; 1170 for(i = 0; i < 100; ++i) 1174 for(i = 0; i < 100; ++i) 1171 { 1175 { 1172 omega = (1.0 + i) * keV; 1176 omega = (1.0 + i) * keV; 1173 if(verboseLevel > 1) 1177 if(verboseLevel > 1) 1174 G4cout << omega / keV << "\t" << GetGas 1178 G4cout << omega / keV << "\t" << GetGasZmuProduct(omega, gamma, varAngle) 1175 << "\t"; 1179 << "\t"; 1176 if(verboseLevel > 0) 1180 if(verboseLevel > 0) 1177 outGas << omega / keV << "\t\t" 1181 outGas << omega / keV << "\t\t" 1178 << GetGasZmuProduct(omega, gamma 1182 << GetGasZmuProduct(omega, gamma, varAngle) << G4endl; 1179 } 1183 } 1180 return; 1184 return; 1181 } 1185 } 1182 1186 1183 ///////////////////////////////////////////// 1187 //////////////////////////////////////////////////////////////////////// 1184 // Computes Compton cross section for plate m 1188 // Computes Compton cross section for plate material in 1/mm 1185 G4double G4VXTRenergyLoss::GetPlateCompton(G4 1189 G4double G4VXTRenergyLoss::GetPlateCompton(G4double omega) 1186 { 1190 { 1187 G4int i, numberOfElements; 1191 G4int i, numberOfElements; 1188 G4double xSection = 0., nowZ, sumZ = 0.; 1192 G4double xSection = 0., nowZ, sumZ = 0.; 1189 1193 1190 const G4MaterialTable* theMaterialTable = G 1194 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 1191 numberOfElements = (G4int)(*theMaterialTabl << 1195 numberOfElements = (*theMaterialTable)[fMatIndex1]->GetNumberOfElements(); 1192 1196 1193 for(i = 0; i < numberOfElements; ++i) 1197 for(i = 0; i < numberOfElements; ++i) 1194 { 1198 { 1195 nowZ = (*theMaterialTable)[fMatIndex1]->G 1199 nowZ = (*theMaterialTable)[fMatIndex1]->GetElement(i)->GetZ(); 1196 sumZ += nowZ; 1200 sumZ += nowZ; 1197 xSection += GetComptonPerAtom(omega, nowZ 1201 xSection += GetComptonPerAtom(omega, nowZ); 1198 } 1202 } 1199 xSection /= sumZ; 1203 xSection /= sumZ; 1200 xSection *= (*theMaterialTable)[fMatIndex1] 1204 xSection *= (*theMaterialTable)[fMatIndex1]->GetElectronDensity(); 1201 return xSection; 1205 return xSection; 1202 } 1206 } 1203 1207 1204 ///////////////////////////////////////////// 1208 //////////////////////////////////////////////////////////////////////// 1205 // Computes Compton cross section for gas mat 1209 // Computes Compton cross section for gas material in 1/mm 1206 G4double G4VXTRenergyLoss::GetGasCompton(G4do 1210 G4double G4VXTRenergyLoss::GetGasCompton(G4double omega) 1207 { 1211 { 1208 G4double xSection = 0., sumZ = 0.; << 1212 G4int i, numberOfElements; >> 1213 G4double xSection = 0., nowZ, sumZ = 0.; 1209 1214 1210 const G4MaterialTable* theMaterialTable = G 1215 const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); 1211 G4int numberOfElements = (G4int)(*theMateri << 1216 numberOfElements = (*theMaterialTable)[fMatIndex2]->GetNumberOfElements(); 1212 1217 1213 for (G4int i = 0; i < numberOfElements; ++i << 1218 for(i = 0; i < numberOfElements; ++i) 1214 { 1219 { 1215 G4double nowZ = (*theMaterialTable)[fMatI << 1220 nowZ = (*theMaterialTable)[fMatIndex2]->GetElement(i)->GetZ(); 1216 sumZ += nowZ; 1221 sumZ += nowZ; 1217 xSection += GetComptonPerAtom(omega, nowZ 1222 xSection += GetComptonPerAtom(omega, nowZ); 1218 } 1223 } 1219 if (sumZ > 0.0) { xSection /= sumZ; } << 1224 xSection /= sumZ; 1220 xSection *= (*theMaterialTable)[fMatIndex2] 1225 xSection *= (*theMaterialTable)[fMatIndex2]->GetElectronDensity(); 1221 return xSection; 1226 return xSection; 1222 } 1227 } 1223 1228 1224 ///////////////////////////////////////////// 1229 //////////////////////////////////////////////////////////////////////// 1225 // Computes Compton cross section per atom wi 1230 // Computes Compton cross section per atom with Z electrons for gamma with 1226 // the energy GammaEnergy 1231 // the energy GammaEnergy 1227 G4double G4VXTRenergyLoss::GetComptonPerAtom( 1232 G4double G4VXTRenergyLoss::GetComptonPerAtom(G4double GammaEnergy, G4double Z) 1228 { 1233 { 1229 G4double CrossSection = 0.0; 1234 G4double CrossSection = 0.0; 1230 if(Z < 0.9999) 1235 if(Z < 0.9999) 1231 return CrossSection; 1236 return CrossSection; 1232 if(GammaEnergy < 0.1 * keV) 1237 if(GammaEnergy < 0.1 * keV) 1233 return CrossSection; 1238 return CrossSection; 1234 if(GammaEnergy > (100. * GeV / Z)) 1239 if(GammaEnergy > (100. * GeV / Z)) 1235 return CrossSection; 1240 return CrossSection; 1236 1241 1237 static constexpr G4double a = 20.0; 1242 static constexpr G4double a = 20.0; 1238 static constexpr G4double b = 230.0; 1243 static constexpr G4double b = 230.0; 1239 static constexpr G4double c = 440.0; 1244 static constexpr G4double c = 440.0; 1240 1245 1241 static constexpr G4double d1 = 2.7965e-1 * 1246 static constexpr G4double d1 = 2.7965e-1 * barn, d2 = -1.8300e-1 * barn, 1242 d3 = 6.7527 * bar 1247 d3 = 6.7527 * barn, d4 = -1.9798e+1 * barn, 1243 e1 = 1.9756e-5 * 1248 e1 = 1.9756e-5 * barn, e2 = -1.0205e-2 * barn, 1244 e3 = -7.3913e-2 * 1249 e3 = -7.3913e-2 * barn, e4 = 2.7079e-2 * barn, 1245 f1 = -3.9178e-7 * 1250 f1 = -3.9178e-7 * barn, f2 = 6.8241e-5 * barn, 1246 f3 = 6.0480e-5 * 1251 f3 = 6.0480e-5 * barn, f4 = 3.0274e-4 * barn; 1247 1252 1248 G4double p1Z = Z * (d1 + e1 * Z + f1 * Z * 1253 G4double p1Z = Z * (d1 + e1 * Z + f1 * Z * Z); 1249 G4double p2Z = Z * (d2 + e2 * Z + f2 * Z * 1254 G4double p2Z = Z * (d2 + e2 * Z + f2 * Z * Z); 1250 G4double p3Z = Z * (d3 + e3 * Z + f3 * Z * 1255 G4double p3Z = Z * (d3 + e3 * Z + f3 * Z * Z); 1251 G4double p4Z = Z * (d4 + e4 * Z + f4 * Z * 1256 G4double p4Z = Z * (d4 + e4 * Z + f4 * Z * Z); 1252 1257 1253 G4double T0 = 15.0 * keV; 1258 G4double T0 = 15.0 * keV; 1254 if(Z < 1.5) 1259 if(Z < 1.5) 1255 T0 = 40.0 * keV; 1260 T0 = 40.0 * keV; 1256 1261 1257 G4double X = std::max(GammaEnergy, T0) / el 1262 G4double X = std::max(GammaEnergy, T0) / electron_mass_c2; 1258 CrossSection = 1263 CrossSection = 1259 p1Z * std::log(1. + 2. * X) / X + 1264 p1Z * std::log(1. + 2. * X) / X + 1260 (p2Z + p3Z * X + p4Z * X * X) / (1. + a * 1265 (p2Z + p3Z * X + p4Z * X * X) / (1. + a * X + b * X * X + c * X * X * X); 1261 1266 1262 // modification for low energy. (special c 1267 // modification for low energy. (special case for Hydrogen) 1263 if(GammaEnergy < T0) 1268 if(GammaEnergy < T0) 1264 { 1269 { 1265 G4double dT0 = 1. * keV; 1270 G4double dT0 = 1. * keV; 1266 X = (T0 + dT0) / electron_mass 1271 X = (T0 + dT0) / electron_mass_c2; 1267 G4double sigma = 1272 G4double sigma = 1268 p1Z * std::log(1. + 2. * X) / X + 1273 p1Z * std::log(1. + 2. * X) / X + 1269 (p2Z + p3Z * X + p4Z * X * X) / (1. + a 1274 (p2Z + p3Z * X + p4Z * X * X) / (1. + a * X + b * X * X + c * X * X * X); 1270 G4double c1 = -T0 * (sigma - CrossSection 1275 G4double c1 = -T0 * (sigma - CrossSection) / (CrossSection * dT0); 1271 G4double c2 = 0.150; 1276 G4double c2 = 0.150; 1272 if(Z > 1.5) 1277 if(Z > 1.5) 1273 c2 = 0.375 - 0.0556 * std::log(Z); 1278 c2 = 0.375 - 0.0556 * std::log(Z); 1274 G4double y = std::log(GammaEnergy / T0); 1279 G4double y = std::log(GammaEnergy / T0); 1275 CrossSection *= std::exp(-y * (c1 + c2 * 1280 CrossSection *= std::exp(-y * (c1 + c2 * y)); 1276 } 1281 } 1277 return CrossSection; 1282 return CrossSection; 1278 } 1283 } 1279 1284 1280 ///////////////////////////////////////////// 1285 /////////////////////////////////////////////////////////////////////// 1281 // This function returns the spectral and ang 1286 // This function returns the spectral and angle density of TR quanta 1282 // in X-ray energy region generated forward w 1287 // in X-ray energy region generated forward when a relativistic 1283 // charged particle crosses interface between 1288 // charged particle crosses interface between two materials. 1284 // The high energy small theta approximation 1289 // The high energy small theta approximation is applied. 1285 // (matter1 -> matter2, or 2->1) 1290 // (matter1 -> matter2, or 2->1) 1286 // varAngle =2* (1 - std::cos(theta)) or appr 1291 // varAngle =2* (1 - std::cos(theta)) or approximately = theta*theta 1287 G4double G4VXTRenergyLoss::OneBoundaryXTRNden 1292 G4double G4VXTRenergyLoss::OneBoundaryXTRNdensity(G4double energy, 1288 1293 G4double gamma, 1289 1294 G4double varAngle) const 1290 { 1295 { 1291 G4double formationLength1, formationLength2 1296 G4double formationLength1, formationLength2; 1292 formationLength1 = 1297 formationLength1 = 1293 1.0 / (1.0 / (gamma * gamma) + fSigma1 / 1298 1.0 / (1.0 / (gamma * gamma) + fSigma1 / (energy * energy) + varAngle); 1294 formationLength2 = 1299 formationLength2 = 1295 1.0 / (1.0 / (gamma * gamma) + fSigma2 / 1300 1.0 / (1.0 / (gamma * gamma) + fSigma2 / (energy * energy) + varAngle); 1296 return (varAngle / energy) * (formationLeng 1301 return (varAngle / energy) * (formationLength1 - formationLength2) * 1297 (formationLength1 - formationLength2 1302 (formationLength1 - formationLength2); 1298 } 1303 } 1299 1304 1300 G4double G4VXTRenergyLoss::GetStackFactor(G4d 1305 G4double G4VXTRenergyLoss::GetStackFactor(G4double energy, G4double gamma, 1301 G4d 1306 G4double varAngle) 1302 { 1307 { 1303 // return stack factor corresponding to one 1308 // return stack factor corresponding to one interface 1304 return std::real(OneInterfaceXTRdEdx(energy 1309 return std::real(OneInterfaceXTRdEdx(energy, gamma, varAngle)); 1305 } 1310 } 1306 1311 1307 ///////////////////////////////////////////// 1312 ////////////////////////////////////////////////////////////////////////////// 1308 // For photon energy distribution tables. Int 1313 // For photon energy distribution tables. Integrate first over angle 1309 G4double G4VXTRenergyLoss::XTRNSpectralAngleD 1314 G4double G4VXTRenergyLoss::XTRNSpectralAngleDensity(G4double varAngle) 1310 { 1315 { 1311 return OneBoundaryXTRNdensity(fEnergy, fGam 1316 return OneBoundaryXTRNdensity(fEnergy, fGamma, varAngle) * 1312 GetStackFactor(fEnergy, fGamma, varA 1317 GetStackFactor(fEnergy, fGamma, varAngle); 1313 } 1318 } 1314 1319 1315 ///////////////////////////////////////////// 1320 ///////////////////////////////////////////////////////////////////////// 1316 // For second integration over energy 1321 // For second integration over energy 1317 G4double G4VXTRenergyLoss::XTRNSpectralDensit 1322 G4double G4VXTRenergyLoss::XTRNSpectralDensity(G4double energy) 1318 { 1323 { 1319 fEnergy = energy; 1324 fEnergy = energy; 1320 G4Integrator<G4VXTRenergyLoss, G4double (G4 1325 G4Integrator<G4VXTRenergyLoss, G4double (G4VXTRenergyLoss::*)(G4double)> 1321 integral; 1326 integral; 1322 return integral.Legendre96(this, &G4VXTRene 1327 return integral.Legendre96(this, &G4VXTRenergyLoss::XTRNSpectralAngleDensity, 1323 0.0, 0.2 * fMaxT 1328 0.0, 0.2 * fMaxThetaTR) + 1324 integral.Legendre10(this, &G4VXTRene 1329 integral.Legendre10(this, &G4VXTRenergyLoss::XTRNSpectralAngleDensity, 1325 0.2 * fMaxThetaT 1330 0.2 * fMaxThetaTR, fMaxThetaTR); 1326 } 1331 } 1327 1332 1328 ///////////////////////////////////////////// 1333 ////////////////////////////////////////////////////////////////////////// 1329 // for photon angle distribution tables 1334 // for photon angle distribution tables 1330 G4double G4VXTRenergyLoss::XTRNAngleSpectralD 1335 G4double G4VXTRenergyLoss::XTRNAngleSpectralDensity(G4double energy) 1331 { 1336 { 1332 return OneBoundaryXTRNdensity(energy, fGamm 1337 return OneBoundaryXTRNdensity(energy, fGamma, fVarAngle) * 1333 GetStackFactor(energy, fGamma, fVarA 1338 GetStackFactor(energy, fGamma, fVarAngle); 1334 } 1339 } 1335 1340 1336 ///////////////////////////////////////////// 1341 /////////////////////////////////////////////////////////////////////////// 1337 G4double G4VXTRenergyLoss::XTRNAngleDensity(G 1342 G4double G4VXTRenergyLoss::XTRNAngleDensity(G4double varAngle) 1338 { 1343 { 1339 fVarAngle = varAngle; 1344 fVarAngle = varAngle; 1340 G4Integrator<G4VXTRenergyLoss, G4double (G4 1345 G4Integrator<G4VXTRenergyLoss, G4double (G4VXTRenergyLoss::*)(G4double)> 1341 integral; 1346 integral; 1342 return integral.Legendre96(this, &G4VXTRene 1347 return integral.Legendre96(this, &G4VXTRenergyLoss::XTRNAngleSpectralDensity, 1343 fMinEnergyTR, fM 1348 fMinEnergyTR, fMaxEnergyTR); 1344 } 1349 } 1345 1350 1346 ///////////////////////////////////////////// 1351 ////////////////////////////////////////////////////////////////////////////// 1347 // Check number of photons for a range of Lor 1352 // Check number of photons for a range of Lorentz factors from both energy 1348 // and angular tables 1353 // and angular tables 1349 void G4VXTRenergyLoss::GetNumberOfPhotons() 1354 void G4VXTRenergyLoss::GetNumberOfPhotons() 1350 { 1355 { 1351 G4int iTkin; 1356 G4int iTkin; 1352 G4double gamma, numberE; 1357 G4double gamma, numberE; 1353 1358 1354 std::ofstream outEn("numberE.dat", std::ios 1359 std::ofstream outEn("numberE.dat", std::ios::out); 1355 outEn.setf(std::ios::scientific, std::ios:: 1360 outEn.setf(std::ios::scientific, std::ios::floatfield); 1356 1361 1357 std::ofstream outAng("numberAng.dat", std:: 1362 std::ofstream outAng("numberAng.dat", std::ios::out); 1358 outAng.setf(std::ios::scientific, std::ios: 1363 outAng.setf(std::ios::scientific, std::ios::floatfield); 1359 1364 1360 for(iTkin = 0; iTkin < fTotBin; ++iTkin) / 1365 for(iTkin = 0; iTkin < fTotBin; ++iTkin) // Lorentz factor loop 1361 { 1366 { 1362 gamma = 1367 gamma = 1363 1.0 + (fProtonEnergyVector->GetLowEdgeE 1368 1.0 + (fProtonEnergyVector->GetLowEdgeEnergy(iTkin) / proton_mass_c2); 1364 numberE = (*(*fEnergyDistrTable)(iTkin))( 1369 numberE = (*(*fEnergyDistrTable)(iTkin))(0); 1365 if(verboseLevel > 1) 1370 if(verboseLevel > 1) 1366 G4cout << gamma << "\t\t" << numberE << 1371 G4cout << gamma << "\t\t" << numberE << "\t" << G4endl; 1367 if(verboseLevel > 0) 1372 if(verboseLevel > 0) 1368 outEn << gamma << "\t\t" << numberE << 1373 outEn << gamma << "\t\t" << numberE << G4endl; 1369 } 1374 } 1370 return; 1375 return; 1371 } 1376 } 1372 1377 1373 ///////////////////////////////////////////// 1378 ///////////////////////////////////////////////////////////////////////// 1374 // Returns random energy of a X-ray TR photon 1379 // Returns random energy of a X-ray TR photon for given scaled kinetic energy 1375 // of a charged particle 1380 // of a charged particle 1376 G4double G4VXTRenergyLoss::GetXTRrandomEnergy 1381 G4double G4VXTRenergyLoss::GetXTRrandomEnergy(G4double scaledTkin, G4int iTkin) 1377 { 1382 { 1378 G4int iTransfer, iPlace; 1383 G4int iTransfer, iPlace; 1379 G4double transfer = 0.0, position, E1, E2, 1384 G4double transfer = 0.0, position, E1, E2, W1, W2, W; 1380 1385 1381 iPlace = iTkin - 1; 1386 iPlace = iTkin - 1; 1382 1387 1383 if(iTkin == fTotBin) // relativistic plato 1388 if(iTkin == fTotBin) // relativistic plato, try from left 1384 { 1389 { 1385 position = (*(*fEnergyDistrTable)(iPlace) 1390 position = (*(*fEnergyDistrTable)(iPlace))(0) * G4UniformRand(); 1386 1391 1387 for(iTransfer = 0;; ++iTransfer) 1392 for(iTransfer = 0;; ++iTransfer) 1388 { 1393 { 1389 if(position >= (*(*fEnergyDistrTable)(i 1394 if(position >= (*(*fEnergyDistrTable)(iPlace))(iTransfer)) 1390 break; 1395 break; 1391 } 1396 } 1392 transfer = GetXTRenergy(iPlace, position, 1397 transfer = GetXTRenergy(iPlace, position, iTransfer); 1393 } 1398 } 1394 else 1399 else 1395 { 1400 { 1396 E1 = fProtonEnergyVector->GetLowEdgeEnerg 1401 E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1); 1397 E2 = fProtonEnergyVector->GetLowEdgeEnerg 1402 E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin); 1398 W = 1.0 / (E2 - E1); 1403 W = 1.0 / (E2 - E1); 1399 W1 = (E2 - scaledTkin) * W; 1404 W1 = (E2 - scaledTkin) * W; 1400 W2 = (scaledTkin - E1) * W; 1405 W2 = (scaledTkin - E1) * W; 1401 1406 1402 position = ((*(*fEnergyDistrTable)(iPlace 1407 position = ((*(*fEnergyDistrTable)(iPlace))(0) * W1 + 1403 (*(*fEnergyDistrTable)(iPlace 1408 (*(*fEnergyDistrTable)(iPlace + 1))(0) * W2) * 1404 G4UniformRand(); 1409 G4UniformRand(); 1405 1410 1406 for(iTransfer = 0;; ++iTransfer) 1411 for(iTransfer = 0;; ++iTransfer) 1407 { 1412 { 1408 if(position >= ((*(*fEnergyDistrTable)( 1413 if(position >= ((*(*fEnergyDistrTable)(iPlace))(iTransfer) *W1 + 1409 (*(*fEnergyDistrTable)( 1414 (*(*fEnergyDistrTable)(iPlace + 1))(iTransfer) *W2)) 1410 break; 1415 break; 1411 } 1416 } 1412 transfer = GetXTRenergy(iPlace, position, 1417 transfer = GetXTRenergy(iPlace, position, iTransfer); 1413 } 1418 } 1414 if(transfer < 0.0) 1419 if(transfer < 0.0) 1415 transfer = 0.0; 1420 transfer = 0.0; 1416 return transfer; 1421 return transfer; 1417 } 1422 } 1418 1423 1419 ///////////////////////////////////////////// 1424 //////////////////////////////////////////////////////////////////////// 1420 // Returns approximate position of X-ray phot 1425 // Returns approximate position of X-ray photon energy during random sampling 1421 // over integral energy distribution 1426 // over integral energy distribution 1422 G4double G4VXTRenergyLoss::GetXTRenergy(G4int 1427 G4double G4VXTRenergyLoss::GetXTRenergy(G4int iPlace, G4double, G4int iTransfer) 1423 { 1428 { 1424 G4double x1, x2, y1, y2, result; 1429 G4double x1, x2, y1, y2, result; 1425 1430 1426 if(iTransfer == 0) 1431 if(iTransfer == 0) 1427 { 1432 { 1428 result = (*fEnergyDistrTable)(iPlace)->Ge 1433 result = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer); 1429 } 1434 } 1430 else 1435 else 1431 { 1436 { 1432 y1 = (*(*fEnergyDistrTable)(iPlace))(iTra 1437 y1 = (*(*fEnergyDistrTable)(iPlace))(iTransfer - 1); 1433 y2 = (*(*fEnergyDistrTable)(iPlace))(iTra 1438 y2 = (*(*fEnergyDistrTable)(iPlace))(iTransfer); 1434 1439 1435 x1 = (*fEnergyDistrTable)(iPlace)->GetLow 1440 x1 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer - 1); 1436 x2 = (*fEnergyDistrTable)(iPlace)->GetLow 1441 x2 = (*fEnergyDistrTable)(iPlace)->GetLowEdgeEnergy(iTransfer); 1437 1442 1438 if(x1 == x2) 1443 if(x1 == x2) 1439 result = x2; 1444 result = x2; 1440 else 1445 else 1441 { 1446 { 1442 if(y1 == y2) 1447 if(y1 == y2) 1443 result = x1 + (x2 - x1) * G4UniformRa 1448 result = x1 + (x2 - x1) * G4UniformRand(); 1444 else 1449 else 1445 { 1450 { 1446 result = x1 + (x2 - x1) * G4UniformRa 1451 result = x1 + (x2 - x1) * G4UniformRand(); 1447 } 1452 } 1448 } 1453 } 1449 } 1454 } 1450 return result; 1455 return result; 1451 } 1456 } 1452 1457 1453 ///////////////////////////////////////////// 1458 ///////////////////////////////////////////////////////////////////////// 1454 // Get XTR photon angle at given energy and 1459 // Get XTR photon angle at given energy and Tkin 1455 1460 1456 G4double G4VXTRenergyLoss::GetRandomAngle(G4d 1461 G4double G4VXTRenergyLoss::GetRandomAngle(G4double energyXTR, G4int iTkin) 1457 { 1462 { 1458 G4int iTR, iAngle; 1463 G4int iTR, iAngle; 1459 G4double position, angle; 1464 G4double position, angle; 1460 1465 1461 if(iTkin == fTotBin) 1466 if(iTkin == fTotBin) 1462 --iTkin; 1467 --iTkin; 1463 1468 1464 fAngleForEnergyTable = fAngleBank[iTkin]; 1469 fAngleForEnergyTable = fAngleBank[iTkin]; 1465 1470 1466 for(iTR = 0; iTR < fBinTR; ++iTR) 1471 for(iTR = 0; iTR < fBinTR; ++iTR) 1467 { 1472 { 1468 if(energyXTR < fXTREnergyVector->GetLowEd 1473 if(energyXTR < fXTREnergyVector->GetLowEdgeEnergy(iTR)) 1469 break; 1474 break; 1470 } 1475 } 1471 if(iTR == fBinTR) 1476 if(iTR == fBinTR) 1472 --iTR; 1477 --iTR; 1473 1478 1474 position = (*(*fAngleForEnergyTable)(iTR))( 1479 position = (*(*fAngleForEnergyTable)(iTR))(0) * G4UniformRand(); 1475 // position = (*(*fAngleForEnergyTable)(iTR 1480 // position = (*(*fAngleForEnergyTable)(iTR))(1) * G4UniformRand(); // ATLAS TB 1476 1481 1477 for(iAngle = 0;; ++iAngle) 1482 for(iAngle = 0;; ++iAngle) 1478 // for(iAngle = 1;; ++iAngle) // ATLAS TB 1483 // for(iAngle = 1;; ++iAngle) // ATLAS TB 1479 { 1484 { 1480 if(position >= (*(*fAngleForEnergyTable)( 1485 if(position >= (*(*fAngleForEnergyTable)(iTR))(iAngle)) 1481 break; 1486 break; 1482 } 1487 } 1483 angle = GetAngleXTR(iTR, position, iAngle); 1488 angle = GetAngleXTR(iTR, position, iAngle); 1484 return angle; 1489 return angle; 1485 } 1490 } 1486 1491 1487 ///////////////////////////////////////////// 1492 //////////////////////////////////////////////////////////////////////// 1488 // Returns approximate position of X-ray phot 1493 // Returns approximate position of X-ray photon angle at given energy during 1489 // random sampling over integral energy distr 1494 // random sampling over integral energy distribution 1490 1495 1491 G4double G4VXTRenergyLoss::GetAngleXTR(G4int 1496 G4double G4VXTRenergyLoss::GetAngleXTR(G4int iPlace, G4double position, 1492 G4int 1497 G4int iTransfer) 1493 { 1498 { 1494 G4double x1, x2, y1, y2, result; 1499 G4double x1, x2, y1, y2, result; 1495 1500 1496 if( iTransfer == 0 ) 1501 if( iTransfer == 0 ) 1497 // if( iTransfer == 1 ) // ATLAS TB 1502 // if( iTransfer == 1 ) // ATLAS TB 1498 { 1503 { 1499 result = (*fAngleForEnergyTable)(iPlace)- 1504 result = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer); 1500 } 1505 } 1501 else 1506 else 1502 { 1507 { 1503 y1 = (*(*fAngleForEnergyTable)(iPlace))(i 1508 y1 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer - 1); 1504 y2 = (*(*fAngleForEnergyTable)(iPlace))(i 1509 y2 = (*(*fAngleForEnergyTable)(iPlace))(iTransfer); 1505 1510 1506 x1 = (*fAngleForEnergyTable)(iPlace)->Get 1511 x1 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer - 1); 1507 x2 = (*fAngleForEnergyTable)(iPlace)->Get 1512 x2 = (*fAngleForEnergyTable)(iPlace)->GetLowEdgeEnergy(iTransfer); 1508 1513 1509 if(x1 == x2) result = x2; 1514 if(x1 == x2) result = x2; 1510 else 1515 else 1511 { 1516 { 1512 if( y1 == y2 ) result = x1 + (x2 - x1) 1517 if( y1 == y2 ) result = x1 + (x2 - x1) * G4UniformRand(); 1513 else 1518 else 1514 { 1519 { 1515 result = x1 + (position - y1) * (x2 - 1520 result = x1 + (position - y1) * (x2 - x1) / (y2 - y1); 1516 // result = x1 + 0.1*(position - y1) 1521 // result = x1 + 0.1*(position - y1) * (x2 - x1) / (y2 - y1); // ATLAS TB 1517 // result = x1 + 0.05*(position - y1) 1522 // result = x1 + 0.05*(position - y1) * (x2 - x1) / (y2 - y1); // ATLAS TB 1518 } 1523 } 1519 } 1524 } 1520 } 1525 } 1521 return result; 1526 return result; 1522 } 1527 } 1523 1528