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