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