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