[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 | // |
---|
[1340] | 26 | // $Id: G4PAIPhotonModel.cc,v 1.25 2010/10/26 09:16:50 vnivanch Exp $ |
---|
| 27 | // GEANT4 tag $Name: emstand-V09-03-24 $ |
---|
[1055] | 28 | // |
---|
| 29 | // ------------------------------------------------------------------- |
---|
| 30 | // |
---|
| 31 | // GEANT4 Class |
---|
[819] | 32 | // File name: G4PAIPhotonModel.cc |
---|
| 33 | // |
---|
| 34 | // Author: Vladimir.Grichine@cern.ch based on G4PAIModel class |
---|
| 35 | // |
---|
| 36 | // Creation date: 20.05.2004 |
---|
| 37 | // |
---|
| 38 | // Modifications: |
---|
| 39 | // |
---|
| 40 | // 17.08.04 V.Grichine, bug fixed for Tkin<=0 in SampleSecondary |
---|
| 41 | // 16.08.04 V.Grichine, bug fixed in massRatio for DEDX, CrossSection, SampleSecondary |
---|
| 42 | // 11.04.05 Major optimisation of internal interfaces (V.Ivantchenko) |
---|
| 43 | // |
---|
| 44 | |
---|
| 45 | #include "G4Region.hh" |
---|
| 46 | #include "G4PhysicsLogVector.hh" |
---|
| 47 | #include "G4PhysicsFreeVector.hh" |
---|
| 48 | #include "G4PhysicsTable.hh" |
---|
| 49 | #include "G4ProductionCutsTable.hh" |
---|
| 50 | #include "G4MaterialCutsCouple.hh" |
---|
| 51 | #include "G4MaterialTable.hh" |
---|
| 52 | #include "G4SandiaTable.hh" |
---|
| 53 | #include "G4PAIxSection.hh" |
---|
| 54 | |
---|
| 55 | #include "G4PAIPhotonModel.hh" |
---|
| 56 | #include "Randomize.hh" |
---|
| 57 | #include "G4Electron.hh" |
---|
| 58 | #include "G4Positron.hh" |
---|
| 59 | #include "G4Gamma.hh" |
---|
| 60 | #include "G4Poisson.hh" |
---|
| 61 | #include "G4Step.hh" |
---|
| 62 | #include "G4Material.hh" |
---|
| 63 | #include "G4DynamicParticle.hh" |
---|
| 64 | #include "G4ParticleDefinition.hh" |
---|
| 65 | #include "G4ParticleChangeForLoss.hh" |
---|
| 66 | #include "G4GeometryTolerance.hh" |
---|
| 67 | |
---|
| 68 | //////////////////////////////////////////////////////////////////////// |
---|
| 69 | |
---|
| 70 | using namespace std; |
---|
| 71 | |
---|
| 72 | G4PAIPhotonModel::G4PAIPhotonModel(const G4ParticleDefinition* p, const G4String& nam) |
---|
| 73 | : G4VEmModel(nam),G4VEmFluctuationModel(nam), |
---|
| 74 | fLowestKineticEnergy(10.0*keV), |
---|
| 75 | fHighestKineticEnergy(100.*TeV), |
---|
| 76 | fTotBin(200), |
---|
| 77 | fMeanNumber(20), |
---|
| 78 | fParticle(0), |
---|
| 79 | fHighKinEnergy(100.*TeV), |
---|
| 80 | fLowKinEnergy(2.0*MeV), |
---|
| 81 | fTwoln10(2.0*log(10.0)), |
---|
| 82 | fBg2lim(0.0169), |
---|
| 83 | fTaulim(8.4146e-3) |
---|
| 84 | { |
---|
| 85 | fVerbose = 0; |
---|
| 86 | fElectron = G4Electron::Electron(); |
---|
| 87 | fPositron = G4Positron::Positron(); |
---|
| 88 | |
---|
| 89 | fProtonEnergyVector = new G4PhysicsLogVector(fLowestKineticEnergy, |
---|
[1340] | 90 | fHighestKineticEnergy, |
---|
| 91 | fTotBin); |
---|
[819] | 92 | fPAItransferTable = 0; |
---|
| 93 | fPAIphotonTable = 0; |
---|
| 94 | fPAIplasmonTable = 0; |
---|
| 95 | |
---|
| 96 | fPAIdEdxTable = 0; |
---|
| 97 | fSandiaPhotoAbsCof = 0; |
---|
| 98 | fdEdxVector = 0; |
---|
| 99 | |
---|
| 100 | fLambdaVector = 0; |
---|
| 101 | fdNdxCutVector = 0; |
---|
| 102 | fdNdxCutPhotonVector = 0; |
---|
| 103 | fdNdxCutPlasmonVector = 0; |
---|
| 104 | |
---|
[1340] | 105 | fSandiaIntervalNumber = 0; |
---|
| 106 | fMatIndex = 0; |
---|
| 107 | |
---|
| 108 | if(p) { SetParticle(p); } |
---|
| 109 | else { SetParticle(fElectron); } |
---|
| 110 | |
---|
[819] | 111 | isInitialised = false; |
---|
| 112 | } |
---|
| 113 | |
---|
| 114 | //////////////////////////////////////////////////////////////////////////// |
---|
| 115 | |
---|
| 116 | G4PAIPhotonModel::~G4PAIPhotonModel() |
---|
| 117 | { |
---|
| 118 | if(fProtonEnergyVector) delete fProtonEnergyVector; |
---|
| 119 | if(fdEdxVector) delete fdEdxVector ; |
---|
| 120 | if ( fLambdaVector) delete fLambdaVector; |
---|
| 121 | if ( fdNdxCutVector) delete fdNdxCutVector; |
---|
| 122 | |
---|
| 123 | if( fPAItransferTable ) |
---|
| 124 | { |
---|
| 125 | fPAItransferTable->clearAndDestroy(); |
---|
| 126 | delete fPAItransferTable ; |
---|
| 127 | } |
---|
| 128 | if( fPAIphotonTable ) |
---|
| 129 | { |
---|
| 130 | fPAIphotonTable->clearAndDestroy(); |
---|
| 131 | delete fPAIphotonTable ; |
---|
| 132 | } |
---|
| 133 | if( fPAIplasmonTable ) |
---|
| 134 | { |
---|
| 135 | fPAIplasmonTable->clearAndDestroy(); |
---|
| 136 | delete fPAIplasmonTable ; |
---|
| 137 | } |
---|
| 138 | if(fSandiaPhotoAbsCof) |
---|
| 139 | { |
---|
| 140 | for(G4int i=0;i<fSandiaIntervalNumber;i++) |
---|
| 141 | { |
---|
| 142 | delete[] fSandiaPhotoAbsCof[i]; |
---|
| 143 | } |
---|
| 144 | delete[] fSandiaPhotoAbsCof; |
---|
| 145 | } |
---|
| 146 | } |
---|
| 147 | |
---|
| 148 | /////////////////////////////////////////////////////////////////////////////// |
---|
| 149 | |
---|
| 150 | void G4PAIPhotonModel::SetParticle(const G4ParticleDefinition* p) |
---|
| 151 | { |
---|
| 152 | fParticle = p; |
---|
| 153 | fMass = fParticle->GetPDGMass(); |
---|
| 154 | fSpin = fParticle->GetPDGSpin(); |
---|
| 155 | G4double q = fParticle->GetPDGCharge()/eplus; |
---|
| 156 | fChargeSquare = q*q; |
---|
| 157 | fLowKinEnergy *= fMass/proton_mass_c2; |
---|
| 158 | fRatio = electron_mass_c2/fMass; |
---|
| 159 | fQc = fMass/fRatio; |
---|
| 160 | } |
---|
| 161 | |
---|
| 162 | //////////////////////////////////////////////////////////////////////////// |
---|
| 163 | |
---|
| 164 | void G4PAIPhotonModel::Initialise(const G4ParticleDefinition* p, |
---|
| 165 | const G4DataVector&) |
---|
| 166 | { |
---|
[1340] | 167 | // G4cout<<"G4PAIPhotonModel::Initialise for "<<p->GetParticleName()<<G4endl; |
---|
| 168 | if(isInitialised) { return; } |
---|
[819] | 169 | isInitialised = true; |
---|
| 170 | |
---|
| 171 | if(!fParticle) SetParticle(p); |
---|
| 172 | |
---|
[1055] | 173 | fParticleChange = GetParticleChangeForLoss(); |
---|
[819] | 174 | |
---|
| 175 | const G4ProductionCutsTable* theCoupleTable = |
---|
| 176 | G4ProductionCutsTable::GetProductionCutsTable(); |
---|
| 177 | |
---|
| 178 | for(size_t iReg = 0; iReg < fPAIRegionVector.size();++iReg) // region loop |
---|
| 179 | { |
---|
| 180 | const G4Region* curReg = fPAIRegionVector[iReg]; |
---|
| 181 | |
---|
| 182 | vector<G4Material*>::const_iterator matIter = curReg->GetMaterialIterator(); |
---|
| 183 | size_t jMat; |
---|
| 184 | size_t numOfMat = curReg->GetNumberOfMaterials(); |
---|
| 185 | |
---|
| 186 | // for(size_t jMat = 0; jMat < curReg->GetNumberOfMaterials();++jMat){} |
---|
| 187 | const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); |
---|
| 188 | size_t numberOfMat = G4Material::GetNumberOfMaterials(); |
---|
| 189 | |
---|
| 190 | for(jMat = 0 ; jMat < numOfMat; ++jMat) // region material loop |
---|
| 191 | { |
---|
| 192 | const G4MaterialCutsCouple* matCouple = theCoupleTable-> |
---|
| 193 | GetMaterialCutsCouple( *matIter, curReg->GetProductionCuts() ); |
---|
| 194 | fMaterialCutsCoupleVector.push_back(matCouple); |
---|
| 195 | |
---|
| 196 | size_t iMatGlob; |
---|
| 197 | for(iMatGlob = 0 ; iMatGlob < numberOfMat ; iMatGlob++ ) |
---|
| 198 | { |
---|
| 199 | if( *matIter == (*theMaterialTable)[iMatGlob]) break ; |
---|
| 200 | } |
---|
| 201 | fMatIndex = iMatGlob; |
---|
| 202 | |
---|
| 203 | ComputeSandiaPhotoAbsCof(); |
---|
| 204 | BuildPAIonisationTable(); |
---|
| 205 | |
---|
| 206 | fPAIxscBank.push_back(fPAItransferTable); |
---|
| 207 | fPAIphotonBank.push_back(fPAIphotonTable); |
---|
| 208 | fPAIplasmonBank.push_back(fPAIplasmonTable); |
---|
| 209 | fPAIdEdxBank.push_back(fPAIdEdxTable); |
---|
| 210 | fdEdxTable.push_back(fdEdxVector); |
---|
| 211 | |
---|
| 212 | BuildLambdaVector(matCouple); |
---|
| 213 | |
---|
| 214 | fdNdxCutTable.push_back(fdNdxCutVector); |
---|
| 215 | fdNdxCutPhotonTable.push_back(fdNdxCutPhotonVector); |
---|
| 216 | fdNdxCutPlasmonTable.push_back(fdNdxCutPlasmonVector); |
---|
| 217 | fLambdaTable.push_back(fLambdaVector); |
---|
| 218 | |
---|
| 219 | |
---|
| 220 | matIter++; |
---|
| 221 | } |
---|
| 222 | } |
---|
| 223 | } |
---|
| 224 | |
---|
| 225 | ////////////////////////////////////////////////////////////////// |
---|
| 226 | |
---|
[1055] | 227 | void G4PAIPhotonModel::InitialiseMe(const G4ParticleDefinition*) |
---|
| 228 | {} |
---|
| 229 | |
---|
| 230 | ////////////////////////////////////////////////////////////////// |
---|
| 231 | |
---|
[819] | 232 | void G4PAIPhotonModel::ComputeSandiaPhotoAbsCof() |
---|
| 233 | { |
---|
| 234 | G4int i, j, numberOfElements ; |
---|
[1340] | 235 | const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable(); |
---|
[819] | 236 | |
---|
| 237 | G4SandiaTable thisMaterialSandiaTable(fMatIndex) ; |
---|
| 238 | numberOfElements = (*theMaterialTable)[fMatIndex]-> |
---|
| 239 | GetNumberOfElements(); |
---|
| 240 | G4int* thisMaterialZ = new G4int[numberOfElements] ; |
---|
| 241 | |
---|
| 242 | for(i=0;i<numberOfElements;i++) |
---|
| 243 | { |
---|
| 244 | thisMaterialZ[i] = |
---|
| 245 | (G4int)(*theMaterialTable)[fMatIndex]->GetElement(i)->GetZ() ; |
---|
| 246 | } |
---|
| 247 | fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaIntervals |
---|
| 248 | (thisMaterialZ,numberOfElements) ; |
---|
| 249 | |
---|
| 250 | fSandiaIntervalNumber = thisMaterialSandiaTable.SandiaMixing |
---|
| 251 | ( thisMaterialZ , |
---|
| 252 | (*theMaterialTable)[fMatIndex]->GetFractionVector() , |
---|
| 253 | numberOfElements,fSandiaIntervalNumber) ; |
---|
| 254 | |
---|
| 255 | fSandiaPhotoAbsCof = new G4double*[fSandiaIntervalNumber] ; |
---|
| 256 | |
---|
| 257 | for(i=0;i<fSandiaIntervalNumber;i++) fSandiaPhotoAbsCof[i] = new G4double[5] ; |
---|
| 258 | |
---|
| 259 | for( i = 0 ; i < fSandiaIntervalNumber ; i++ ) |
---|
| 260 | { |
---|
| 261 | fSandiaPhotoAbsCof[i][0] = thisMaterialSandiaTable.GetPhotoAbsorpCof(i+1,0) ; |
---|
| 262 | |
---|
| 263 | for( j = 1; j < 5 ; j++ ) |
---|
| 264 | { |
---|
| 265 | fSandiaPhotoAbsCof[i][j] = thisMaterialSandiaTable. |
---|
| 266 | GetPhotoAbsorpCof(i+1,j)* |
---|
| 267 | (*theMaterialTable)[fMatIndex]->GetDensity() ; |
---|
| 268 | } |
---|
| 269 | } |
---|
[1340] | 270 | delete[] thisMaterialZ ; |
---|
[819] | 271 | } |
---|
| 272 | |
---|
| 273 | //////////////////////////////////////////////////////////////////////////// |
---|
| 274 | // |
---|
| 275 | // Build tables for the ionization energy loss |
---|
| 276 | // the tables are built for MATERIALS |
---|
| 277 | // ********* |
---|
| 278 | |
---|
| 279 | void |
---|
| 280 | G4PAIPhotonModel::BuildPAIonisationTable() |
---|
| 281 | { |
---|
| 282 | G4double LowEdgeEnergy , ionloss ; |
---|
| 283 | G4double massRatio, tau, Tmax, Tmin, Tkin, deltaLow, gamma, bg2 ; |
---|
| 284 | /* |
---|
| 285 | if( fPAItransferTable ) |
---|
| 286 | { |
---|
| 287 | fPAItransferTable->clearAndDestroy() ; |
---|
| 288 | delete fPAItransferTable ; |
---|
| 289 | } |
---|
| 290 | */ |
---|
| 291 | fPAItransferTable = new G4PhysicsTable(fTotBin); |
---|
| 292 | /* |
---|
| 293 | if( fPAIratioTable ) |
---|
| 294 | { |
---|
| 295 | fPAIratioTable->clearAndDestroy() ; |
---|
| 296 | delete fPAIratioTable ; |
---|
| 297 | } |
---|
| 298 | */ |
---|
| 299 | fPAIphotonTable = new G4PhysicsTable(fTotBin); |
---|
| 300 | fPAIplasmonTable = new G4PhysicsTable(fTotBin); |
---|
| 301 | /* |
---|
| 302 | if( fPAIdEdxTable ) |
---|
| 303 | { |
---|
| 304 | fPAIdEdxTable->clearAndDestroy() ; |
---|
| 305 | delete fPAIdEdxTable ; |
---|
| 306 | } |
---|
| 307 | */ |
---|
| 308 | fPAIdEdxTable = new G4PhysicsTable(fTotBin); |
---|
| 309 | |
---|
| 310 | // if(fdEdxVector) delete fdEdxVector ; |
---|
| 311 | fdEdxVector = new G4PhysicsLogVector( fLowestKineticEnergy, |
---|
| 312 | fHighestKineticEnergy, |
---|
| 313 | fTotBin ) ; |
---|
| 314 | Tmin = fSandiaPhotoAbsCof[0][0] ; // low energy Sandia interval |
---|
| 315 | deltaLow = 100.*eV; // 0.5*eV ; |
---|
| 316 | |
---|
[1196] | 317 | for (G4int i = 0 ; i <= fTotBin ; i++) //The loop for the kinetic energy |
---|
[819] | 318 | { |
---|
| 319 | LowEdgeEnergy = fProtonEnergyVector->GetLowEdgeEnergy(i) ; |
---|
| 320 | tau = LowEdgeEnergy/proton_mass_c2 ; |
---|
| 321 | // if(tau < 0.01) tau = 0.01 ; |
---|
| 322 | gamma = tau +1. ; |
---|
| 323 | // G4cout<<"gamma = "<<gamma<<endl ; |
---|
| 324 | bg2 = tau*(tau + 2. ) ; |
---|
| 325 | massRatio = electron_mass_c2/proton_mass_c2 ; |
---|
| 326 | Tmax = MaxSecondaryEnergy(fParticle, LowEdgeEnergy); |
---|
| 327 | // G4cout<<"proton Tkin = "<<LowEdgeEnergy/MeV<<" MeV" |
---|
| 328 | // <<" Tmax = "<<Tmax/MeV<<" MeV"<<G4endl; |
---|
| 329 | // Tkin = DeltaCutInKineticEnergyNow ; |
---|
| 330 | |
---|
| 331 | // if ( DeltaCutInKineticEnergyNow > Tmax) // was < |
---|
| 332 | Tkin = Tmax ; |
---|
| 333 | if ( Tkin < Tmin + deltaLow ) // low energy safety |
---|
| 334 | { |
---|
| 335 | Tkin = Tmin + deltaLow ; |
---|
| 336 | } |
---|
| 337 | G4PAIxSection protonPAI( fMatIndex, |
---|
| 338 | Tkin, |
---|
| 339 | bg2, |
---|
| 340 | fSandiaPhotoAbsCof, |
---|
| 341 | fSandiaIntervalNumber ) ; |
---|
| 342 | |
---|
| 343 | |
---|
| 344 | // G4cout<<"ionloss = "<<ionloss*cm/keV<<" keV/cm"<<endl ; |
---|
| 345 | // G4cout<<"n1 = "<<protonPAI.GetIntegralPAIxSection(1)*cm<<" 1/cm"<<endl ; |
---|
| 346 | // G4cout<<"protonPAI.GetSplineSize() = "<< |
---|
| 347 | // protonPAI.GetSplineSize()<<G4endl<<G4endl ; |
---|
| 348 | |
---|
| 349 | G4PhysicsFreeVector* transferVector = new |
---|
| 350 | G4PhysicsFreeVector(protonPAI.GetSplineSize()) ; |
---|
| 351 | G4PhysicsFreeVector* photonVector = new |
---|
| 352 | G4PhysicsFreeVector(protonPAI.GetSplineSize()) ; |
---|
| 353 | G4PhysicsFreeVector* plasmonVector = new |
---|
| 354 | G4PhysicsFreeVector(protonPAI.GetSplineSize()) ; |
---|
| 355 | G4PhysicsFreeVector* dEdxVector = new |
---|
| 356 | G4PhysicsFreeVector(protonPAI.GetSplineSize()) ; |
---|
| 357 | |
---|
| 358 | for( G4int k = 0 ; k < protonPAI.GetSplineSize() ; k++ ) |
---|
| 359 | { |
---|
| 360 | transferVector->PutValue( k , |
---|
| 361 | protonPAI.GetSplineEnergy(k+1), |
---|
| 362 | protonPAI.GetIntegralPAIxSection(k+1) ) ; |
---|
| 363 | photonVector->PutValue( k , |
---|
| 364 | protonPAI.GetSplineEnergy(k+1), |
---|
| 365 | protonPAI.GetIntegralCerenkov(k+1) ) ; |
---|
| 366 | plasmonVector->PutValue( k , |
---|
| 367 | protonPAI.GetSplineEnergy(k+1), |
---|
| 368 | protonPAI.GetIntegralPlasmon(k+1) ) ; |
---|
| 369 | dEdxVector->PutValue( k , |
---|
| 370 | protonPAI.GetSplineEnergy(k+1), |
---|
| 371 | protonPAI.GetIntegralPAIdEdx(k+1) ) ; |
---|
| 372 | } |
---|
| 373 | ionloss = protonPAI.GetMeanEnergyLoss() ; // total <dE/dx> |
---|
| 374 | if ( ionloss <= 0.) ionloss = DBL_MIN ; |
---|
| 375 | fdEdxVector->PutValue(i,ionloss) ; |
---|
| 376 | |
---|
| 377 | fPAItransferTable->insertAt(i,transferVector) ; |
---|
| 378 | fPAIphotonTable->insertAt(i,photonVector) ; |
---|
| 379 | fPAIplasmonTable->insertAt(i,plasmonVector) ; |
---|
| 380 | fPAIdEdxTable->insertAt(i,dEdxVector) ; |
---|
| 381 | |
---|
| 382 | } // end of Tkin loop |
---|
| 383 | // theLossTable->insert(fdEdxVector); |
---|
| 384 | // end of material loop |
---|
| 385 | // G4cout<<"G4PAIonisation::BuildPAIonisationTable() have been called"<<G4endl ; |
---|
| 386 | // G4cout<<"G4PAIonisation::BuildLossTable() have been called"<<G4endl ; |
---|
| 387 | } |
---|
| 388 | |
---|
| 389 | /////////////////////////////////////////////////////////////////////// |
---|
| 390 | // |
---|
| 391 | // Build mean free path tables for the delta ray production process |
---|
| 392 | // tables are built for MATERIALS |
---|
| 393 | // |
---|
| 394 | |
---|
| 395 | void |
---|
| 396 | G4PAIPhotonModel::BuildLambdaVector(const G4MaterialCutsCouple* matCutsCouple) |
---|
| 397 | { |
---|
| 398 | G4int i ; |
---|
| 399 | G4double dNdxCut,dNdxPhotonCut,dNdxPlasmonCut, lambda; |
---|
| 400 | G4double kCarTolerance = G4GeometryTolerance::GetInstance() |
---|
| 401 | ->GetSurfaceTolerance(); |
---|
| 402 | |
---|
| 403 | const G4ProductionCutsTable* theCoupleTable= |
---|
| 404 | G4ProductionCutsTable::GetProductionCutsTable(); |
---|
| 405 | |
---|
| 406 | size_t numOfCouples = theCoupleTable->GetTableSize(); |
---|
| 407 | size_t jMatCC; |
---|
| 408 | |
---|
| 409 | for (jMatCC = 0 ; jMatCC < numOfCouples ; jMatCC++ ) |
---|
| 410 | { |
---|
| 411 | if( matCutsCouple == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break; |
---|
| 412 | } |
---|
| 413 | if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--; |
---|
| 414 | |
---|
| 415 | const vector<G4double>* deltaCutInKineticEnergy = theCoupleTable-> |
---|
| 416 | GetEnergyCutsVector(idxG4ElectronCut); |
---|
| 417 | const vector<G4double>* photonCutInKineticEnergy = theCoupleTable-> |
---|
| 418 | GetEnergyCutsVector(idxG4GammaCut); |
---|
| 419 | |
---|
| 420 | if (fLambdaVector) delete fLambdaVector; |
---|
| 421 | if (fdNdxCutVector) delete fdNdxCutVector; |
---|
| 422 | if (fdNdxCutPhotonVector) delete fdNdxCutPhotonVector; |
---|
| 423 | if (fdNdxCutPlasmonVector) delete fdNdxCutPlasmonVector; |
---|
| 424 | |
---|
| 425 | fLambdaVector = new G4PhysicsLogVector( fLowestKineticEnergy, |
---|
| 426 | fHighestKineticEnergy, |
---|
| 427 | fTotBin ); |
---|
| 428 | fdNdxCutVector = new G4PhysicsLogVector( fLowestKineticEnergy, |
---|
| 429 | fHighestKineticEnergy, |
---|
| 430 | fTotBin ); |
---|
| 431 | fdNdxCutPhotonVector = new G4PhysicsLogVector( fLowestKineticEnergy, |
---|
| 432 | fHighestKineticEnergy, |
---|
| 433 | fTotBin ); |
---|
| 434 | fdNdxCutPlasmonVector = new G4PhysicsLogVector( fLowestKineticEnergy, |
---|
| 435 | fHighestKineticEnergy, |
---|
| 436 | fTotBin ); |
---|
| 437 | |
---|
| 438 | G4double deltaCutInKineticEnergyNow = (*deltaCutInKineticEnergy)[jMatCC]; |
---|
| 439 | G4double photonCutInKineticEnergyNow = (*photonCutInKineticEnergy)[jMatCC]; |
---|
| 440 | |
---|
| 441 | if(fVerbose > 0) |
---|
| 442 | { |
---|
| 443 | G4cout<<"PAIPhotonModel deltaCutInKineticEnergyNow = " |
---|
| 444 | <<deltaCutInKineticEnergyNow/keV<<" keV"<<G4endl; |
---|
| 445 | G4cout<<"PAIPhotonModel photonCutInKineticEnergyNow = " |
---|
| 446 | <<photonCutInKineticEnergyNow/keV<<" keV"<<G4endl; |
---|
| 447 | } |
---|
[1196] | 448 | for ( i = 0 ; i <= fTotBin ; i++ ) |
---|
[819] | 449 | { |
---|
| 450 | dNdxPhotonCut = GetdNdxPhotonCut(i,photonCutInKineticEnergyNow); |
---|
| 451 | dNdxPlasmonCut = GetdNdxPlasmonCut(i,deltaCutInKineticEnergyNow); |
---|
| 452 | |
---|
| 453 | dNdxCut = dNdxPhotonCut + dNdxPlasmonCut; |
---|
| 454 | lambda = dNdxCut <= DBL_MIN ? DBL_MAX: 1.0/dNdxCut; |
---|
| 455 | |
---|
| 456 | if (lambda <= 1000*kCarTolerance) lambda = 1000*kCarTolerance; // Mmm ??? |
---|
| 457 | |
---|
| 458 | fLambdaVector->PutValue(i, lambda); |
---|
| 459 | |
---|
| 460 | fdNdxCutVector->PutValue(i, dNdxCut); |
---|
| 461 | fdNdxCutPhotonVector->PutValue(i, dNdxPhotonCut); |
---|
| 462 | fdNdxCutPlasmonVector->PutValue(i, dNdxPlasmonCut); |
---|
| 463 | } |
---|
| 464 | } |
---|
| 465 | |
---|
| 466 | /////////////////////////////////////////////////////////////////////// |
---|
| 467 | // |
---|
| 468 | // Returns integral PAI cross section for energy transfers >= transferCut |
---|
| 469 | |
---|
| 470 | G4double |
---|
| 471 | G4PAIPhotonModel::GetdNdxCut( G4int iPlace, G4double transferCut) |
---|
| 472 | { |
---|
| 473 | G4int iTransfer; |
---|
| 474 | G4double x1, x2, y1, y2, dNdxCut; |
---|
| 475 | // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl; |
---|
| 476 | // G4cout<<"size = "<<G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) |
---|
| 477 | // <<G4endl; |
---|
| 478 | for( iTransfer = 0 ; |
---|
| 479 | iTransfer < G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ; |
---|
| 480 | iTransfer++) |
---|
| 481 | { |
---|
| 482 | if(transferCut <= (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer)) |
---|
| 483 | { |
---|
| 484 | break ; |
---|
| 485 | } |
---|
| 486 | } |
---|
| 487 | if ( iTransfer >= G4int((*fPAItransferTable)(iPlace)->GetVectorLength()) ) |
---|
| 488 | { |
---|
| 489 | iTransfer = (*fPAItransferTable)(iPlace)->GetVectorLength() - 1 ; |
---|
| 490 | } |
---|
| 491 | y1 = (*(*fPAItransferTable)(iPlace))(iTransfer-1) ; |
---|
| 492 | y2 = (*(*fPAItransferTable)(iPlace))(iTransfer) ; |
---|
| 493 | // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; |
---|
| 494 | x1 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; |
---|
| 495 | x2 = (*fPAItransferTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; |
---|
| 496 | // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; |
---|
| 497 | |
---|
| 498 | if ( y1 == y2 ) dNdxCut = y2 ; |
---|
| 499 | else |
---|
| 500 | { |
---|
| 501 | // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; |
---|
[1055] | 502 | // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; |
---|
| 503 | if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ; |
---|
[819] | 504 | else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; |
---|
| 505 | } |
---|
| 506 | // G4cout<<""<<dNdxCut<<G4endl; |
---|
| 507 | return dNdxCut ; |
---|
| 508 | } |
---|
| 509 | |
---|
| 510 | /////////////////////////////////////////////////////////////////////// |
---|
| 511 | // |
---|
| 512 | // Returns integral PAI cherenkovcross section for energy transfers >= transferCut |
---|
| 513 | |
---|
| 514 | G4double |
---|
| 515 | G4PAIPhotonModel::GetdNdxPhotonCut( G4int iPlace, G4double transferCut) |
---|
| 516 | { |
---|
| 517 | G4int iTransfer; |
---|
| 518 | G4double x1, x2, y1, y2, dNdxCut; |
---|
| 519 | // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl; |
---|
| 520 | // G4cout<<"size = "<<G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) |
---|
| 521 | // <<G4endl; |
---|
| 522 | for( iTransfer = 0 ; |
---|
| 523 | iTransfer < G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) ; |
---|
| 524 | iTransfer++) |
---|
| 525 | { |
---|
| 526 | if(transferCut <= (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer)) |
---|
| 527 | { |
---|
| 528 | break ; |
---|
| 529 | } |
---|
| 530 | } |
---|
| 531 | if ( iTransfer >= G4int((*fPAIphotonTable)(iPlace)->GetVectorLength()) ) |
---|
| 532 | { |
---|
| 533 | iTransfer = (*fPAIphotonTable)(iPlace)->GetVectorLength() - 1 ; |
---|
| 534 | } |
---|
| 535 | y1 = (*(*fPAIphotonTable)(iPlace))(iTransfer-1) ; |
---|
| 536 | y2 = (*(*fPAIphotonTable)(iPlace))(iTransfer) ; |
---|
| 537 | // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; |
---|
| 538 | x1 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; |
---|
| 539 | x2 = (*fPAIphotonTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; |
---|
| 540 | // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; |
---|
| 541 | |
---|
| 542 | if ( y1 == y2 ) dNdxCut = y2 ; |
---|
| 543 | else |
---|
| 544 | { |
---|
| 545 | // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; |
---|
[1055] | 546 | // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; |
---|
| 547 | if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ; |
---|
[819] | 548 | else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; |
---|
| 549 | } |
---|
| 550 | // G4cout<<""<<dNdxPhotonCut<<G4endl; |
---|
| 551 | return dNdxCut ; |
---|
| 552 | } |
---|
| 553 | |
---|
| 554 | /////////////////////////////////////////////////////////////////////// |
---|
| 555 | // |
---|
| 556 | // Returns integral PAI cross section for energy transfers >= transferCut |
---|
| 557 | |
---|
| 558 | G4double |
---|
| 559 | G4PAIPhotonModel::GetdNdxPlasmonCut( G4int iPlace, G4double transferCut) |
---|
| 560 | { |
---|
| 561 | G4int iTransfer; |
---|
| 562 | G4double x1, x2, y1, y2, dNdxCut; |
---|
| 563 | |
---|
| 564 | // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl; |
---|
| 565 | // G4cout<<"size = "<<G4int((*fPAIPlasmonTable)(iPlace)->GetVectorLength()) |
---|
| 566 | // <<G4endl; |
---|
| 567 | for( iTransfer = 0 ; |
---|
| 568 | iTransfer < G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength()) ; |
---|
| 569 | iTransfer++) |
---|
| 570 | { |
---|
| 571 | if(transferCut <= (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer)) |
---|
| 572 | { |
---|
| 573 | break ; |
---|
| 574 | } |
---|
| 575 | } |
---|
| 576 | if ( iTransfer >= G4int((*fPAIplasmonTable)(iPlace)->GetVectorLength()) ) |
---|
| 577 | { |
---|
| 578 | iTransfer = (*fPAIplasmonTable)(iPlace)->GetVectorLength() - 1 ; |
---|
| 579 | } |
---|
| 580 | y1 = (*(*fPAIplasmonTable)(iPlace))(iTransfer-1) ; |
---|
| 581 | y2 = (*(*fPAIplasmonTable)(iPlace))(iTransfer) ; |
---|
| 582 | // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; |
---|
| 583 | x1 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; |
---|
| 584 | x2 = (*fPAIplasmonTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; |
---|
| 585 | // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; |
---|
| 586 | |
---|
| 587 | if ( y1 == y2 ) dNdxCut = y2 ; |
---|
| 588 | else |
---|
| 589 | { |
---|
| 590 | // if ( x1 == x2 ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; |
---|
[1055] | 591 | // if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*G4UniformRand() ; |
---|
| 592 | if ( std::abs(x1-x2) <= eV ) dNdxCut = y1 + (y2 - y1)*0.5 ; |
---|
[819] | 593 | else dNdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; |
---|
| 594 | } |
---|
| 595 | // G4cout<<""<<dNdxPlasmonCut<<G4endl; |
---|
| 596 | return dNdxCut ; |
---|
| 597 | } |
---|
| 598 | |
---|
| 599 | /////////////////////////////////////////////////////////////////////// |
---|
| 600 | // |
---|
| 601 | // Returns integral dEdx for energy transfers >= transferCut |
---|
| 602 | |
---|
| 603 | G4double |
---|
| 604 | G4PAIPhotonModel::GetdEdxCut( G4int iPlace, G4double transferCut) |
---|
| 605 | { |
---|
| 606 | G4int iTransfer; |
---|
| 607 | G4double x1, x2, y1, y2, dEdxCut; |
---|
| 608 | // G4cout<<"iPlace = "<<iPlace<<"; "<<"transferCut = "<<transferCut<<G4endl; |
---|
| 609 | // G4cout<<"size = "<<G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) |
---|
| 610 | // <<G4endl; |
---|
| 611 | for( iTransfer = 0 ; |
---|
| 612 | iTransfer < G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ; |
---|
| 613 | iTransfer++) |
---|
| 614 | { |
---|
| 615 | if(transferCut <= (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer)) |
---|
| 616 | { |
---|
| 617 | break ; |
---|
| 618 | } |
---|
| 619 | } |
---|
| 620 | if ( iTransfer >= G4int((*fPAIdEdxTable)(iPlace)->GetVectorLength()) ) |
---|
| 621 | { |
---|
| 622 | iTransfer = (*fPAIdEdxTable)(iPlace)->GetVectorLength() - 1 ; |
---|
| 623 | } |
---|
| 624 | y1 = (*(*fPAIdEdxTable)(iPlace))(iTransfer-1) ; |
---|
| 625 | y2 = (*(*fPAIdEdxTable)(iPlace))(iTransfer) ; |
---|
| 626 | // G4cout<<"y1 = "<<y1<<"; "<<"y2 = "<<y2<<G4endl; |
---|
| 627 | x1 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1) ; |
---|
| 628 | x2 = (*fPAIdEdxTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; |
---|
| 629 | // G4cout<<"x1 = "<<x1<<"; "<<"x2 = "<<x2<<G4endl; |
---|
| 630 | |
---|
| 631 | if ( y1 == y2 ) dEdxCut = y2 ; |
---|
| 632 | else |
---|
| 633 | { |
---|
| 634 | // if ( x1 == x2 ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; |
---|
[1055] | 635 | // if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*G4UniformRand() ; |
---|
| 636 | if ( std::abs(x1-x2) <= eV ) dEdxCut = y1 + (y2 - y1)*0.5 ; |
---|
[819] | 637 | else dEdxCut = y1 + (transferCut - x1)*(y2 - y1)/(x2 - x1) ; |
---|
| 638 | } |
---|
| 639 | // G4cout<<""<<dEdxCut<<G4endl; |
---|
| 640 | return dEdxCut ; |
---|
| 641 | } |
---|
| 642 | |
---|
| 643 | ////////////////////////////////////////////////////////////////////////////// |
---|
| 644 | |
---|
[1055] | 645 | G4double G4PAIPhotonModel::ComputeDEDXPerVolume(const G4Material*, |
---|
| 646 | const G4ParticleDefinition* p, |
---|
| 647 | G4double kineticEnergy, |
---|
| 648 | G4double cutEnergy) |
---|
[819] | 649 | { |
---|
| 650 | G4int iTkin,iPlace; |
---|
| 651 | size_t jMat; |
---|
[1055] | 652 | |
---|
| 653 | //G4double cut = std::min(MaxSecondaryEnergy(p, kineticEnergy), cutEnergy); |
---|
| 654 | G4double cut = cutEnergy; |
---|
| 655 | |
---|
[819] | 656 | G4double particleMass = p->GetPDGMass(); |
---|
| 657 | G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass; |
---|
| 658 | G4double charge = p->GetPDGCharge()/eplus; |
---|
| 659 | G4double charge2 = charge*charge; |
---|
| 660 | G4double dEdx = 0.; |
---|
[1055] | 661 | const G4MaterialCutsCouple* matCC = CurrentCouple(); |
---|
[819] | 662 | |
---|
| 663 | for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) |
---|
| 664 | { |
---|
| 665 | if( matCC == fMaterialCutsCoupleVector[jMat] ) break; |
---|
| 666 | } |
---|
| 667 | if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--; |
---|
| 668 | |
---|
| 669 | fPAIdEdxTable = fPAIdEdxBank[jMat]; |
---|
| 670 | fdEdxVector = fdEdxTable[jMat]; |
---|
[1196] | 671 | for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++) |
---|
[819] | 672 | { |
---|
| 673 | if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; |
---|
| 674 | } |
---|
| 675 | iPlace = iTkin - 1; |
---|
| 676 | if(iPlace < 0) iPlace = 0; |
---|
[1055] | 677 | dEdx = charge2*( (*fdEdxVector)(iPlace) - GetdEdxCut(iPlace,cut) ) ; |
---|
[819] | 678 | |
---|
| 679 | if( dEdx < 0.) dEdx = 0.; |
---|
| 680 | return dEdx; |
---|
| 681 | } |
---|
| 682 | |
---|
| 683 | ///////////////////////////////////////////////////////////////////////// |
---|
| 684 | |
---|
[1055] | 685 | G4double G4PAIPhotonModel::CrossSectionPerVolume( const G4Material*, |
---|
| 686 | const G4ParticleDefinition* p, |
---|
| 687 | G4double kineticEnergy, |
---|
| 688 | G4double cutEnergy, |
---|
| 689 | G4double maxEnergy ) |
---|
[819] | 690 | { |
---|
| 691 | G4int iTkin,iPlace; |
---|
| 692 | size_t jMat, jMatCC; |
---|
[1055] | 693 | G4double tmax = std::min(MaxSecondaryEnergy(p, kineticEnergy), maxEnergy); |
---|
| 694 | if(cutEnergy >= tmax) return 0.0; |
---|
[819] | 695 | G4double particleMass = p->GetPDGMass(); |
---|
| 696 | G4double scaledTkin = kineticEnergy*proton_mass_c2/particleMass; |
---|
| 697 | G4double charge = p->GetPDGCharge(); |
---|
| 698 | G4double charge2 = charge*charge, cross, cross1, cross2; |
---|
| 699 | G4double photon1, photon2, plasmon1, plasmon2; |
---|
| 700 | |
---|
[1055] | 701 | const G4MaterialCutsCouple* matCC = CurrentCouple(); |
---|
| 702 | |
---|
[819] | 703 | const G4ProductionCutsTable* theCoupleTable= |
---|
| 704 | G4ProductionCutsTable::GetProductionCutsTable(); |
---|
| 705 | |
---|
| 706 | size_t numOfCouples = theCoupleTable->GetTableSize(); |
---|
| 707 | |
---|
| 708 | for (jMatCC = 0 ; jMatCC < numOfCouples ; jMatCC++ ) |
---|
| 709 | { |
---|
| 710 | if( matCC == theCoupleTable->GetMaterialCutsCouple(jMatCC) ) break; |
---|
| 711 | } |
---|
| 712 | if( jMatCC == numOfCouples && jMatCC > 0 ) jMatCC--; |
---|
| 713 | |
---|
| 714 | const vector<G4double>* photonCutInKineticEnergy = theCoupleTable-> |
---|
| 715 | GetEnergyCutsVector(idxG4GammaCut); |
---|
| 716 | |
---|
| 717 | G4double photonCut = (*photonCutInKineticEnergy)[jMatCC] ; |
---|
| 718 | |
---|
| 719 | for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) |
---|
| 720 | { |
---|
| 721 | if( matCC == fMaterialCutsCoupleVector[jMat] ) break; |
---|
| 722 | } |
---|
| 723 | if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--; |
---|
| 724 | |
---|
| 725 | fPAItransferTable = fPAIxscBank[jMat]; |
---|
| 726 | fPAIphotonTable = fPAIphotonBank[jMat]; |
---|
| 727 | fPAIplasmonTable = fPAIplasmonBank[jMat]; |
---|
| 728 | |
---|
[1196] | 729 | for(iTkin = 0 ; iTkin <= fTotBin ; iTkin++) |
---|
[819] | 730 | { |
---|
| 731 | if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; |
---|
| 732 | } |
---|
| 733 | iPlace = iTkin - 1; |
---|
| 734 | if(iPlace < 0) iPlace = 0; |
---|
| 735 | |
---|
| 736 | // G4cout<<"iPlace = "<<iPlace<<"; tmax = " |
---|
| 737 | // <<tmax<<"; cutEnergy = "<<cutEnergy<<G4endl; |
---|
| 738 | photon1 = GetdNdxPhotonCut(iPlace,tmax); |
---|
| 739 | photon2 = GetdNdxPhotonCut(iPlace,photonCut); |
---|
| 740 | |
---|
| 741 | plasmon1 = GetdNdxPlasmonCut(iPlace,tmax); |
---|
| 742 | plasmon2 = GetdNdxPlasmonCut(iPlace,cutEnergy); |
---|
| 743 | |
---|
| 744 | cross1 = photon1 + plasmon1; |
---|
| 745 | // G4cout<<"cross1 = "<<cross1<<G4endl; |
---|
| 746 | cross2 = photon2 + plasmon2; |
---|
| 747 | // G4cout<<"cross2 = "<<cross2<<G4endl; |
---|
| 748 | cross = (cross2 - cross1)*charge2; |
---|
| 749 | // G4cout<<"cross = "<<cross<<G4endl; |
---|
| 750 | |
---|
| 751 | if( cross < 0. ) cross = 0.; |
---|
| 752 | return cross; |
---|
| 753 | } |
---|
| 754 | |
---|
| 755 | /////////////////////////////////////////////////////////////////////////// |
---|
| 756 | // |
---|
| 757 | // It is analog of PostStepDoIt in terms of secondary electron or photon to |
---|
| 758 | // be returned as G4Dynamicparticle*. |
---|
| 759 | // |
---|
| 760 | |
---|
| 761 | void G4PAIPhotonModel::SampleSecondaries(std::vector<G4DynamicParticle*>* vdp, |
---|
| 762 | const G4MaterialCutsCouple* matCC, |
---|
| 763 | const G4DynamicParticle* dp, |
---|
| 764 | G4double tmin, |
---|
| 765 | G4double maxEnergy) |
---|
| 766 | { |
---|
| 767 | size_t jMat; |
---|
| 768 | for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) |
---|
| 769 | { |
---|
| 770 | if( matCC == fMaterialCutsCoupleVector[jMat] ) break; |
---|
| 771 | } |
---|
| 772 | if( jMat == fMaterialCutsCoupleVector.size() && jMat > 0 ) jMat--; |
---|
| 773 | |
---|
| 774 | fPAItransferTable = fPAIxscBank[jMat]; |
---|
| 775 | fPAIphotonTable = fPAIphotonBank[jMat]; |
---|
| 776 | fPAIplasmonTable = fPAIplasmonBank[jMat]; |
---|
| 777 | |
---|
| 778 | fdNdxCutVector = fdNdxCutTable[jMat]; |
---|
| 779 | fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat]; |
---|
| 780 | fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat]; |
---|
| 781 | |
---|
| 782 | G4double tmax = min(MaxSecondaryKinEnergy(dp), maxEnergy); |
---|
| 783 | if( tmin >= tmax && fVerbose > 0) |
---|
| 784 | { |
---|
| 785 | G4cout<<"G4PAIPhotonModel::SampleSecondary: tmin >= tmax "<<G4endl; |
---|
| 786 | } |
---|
| 787 | |
---|
| 788 | G4ThreeVector direction = dp->GetMomentumDirection(); |
---|
| 789 | G4double particleMass = dp->GetMass(); |
---|
| 790 | G4double kineticEnergy = dp->GetKineticEnergy(); |
---|
| 791 | G4double scaledTkin = kineticEnergy*fMass/particleMass; |
---|
| 792 | G4double totalEnergy = kineticEnergy + particleMass; |
---|
| 793 | G4double pSquare = kineticEnergy*(totalEnergy+particleMass); |
---|
| 794 | |
---|
| 795 | G4int iTkin; |
---|
[1196] | 796 | for(iTkin=0;iTkin<=fTotBin;iTkin++) |
---|
[819] | 797 | { |
---|
| 798 | if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; |
---|
| 799 | } |
---|
| 800 | G4int iPlace = iTkin - 1 ; |
---|
| 801 | if(iPlace < 0) iPlace = 0; |
---|
| 802 | |
---|
| 803 | G4double dNdxPhotonCut = (*fdNdxCutPhotonVector)(iPlace) ; |
---|
| 804 | G4double dNdxPlasmonCut = (*fdNdxCutPlasmonVector)(iPlace) ; |
---|
| 805 | G4double dNdxCut = dNdxPhotonCut + dNdxPlasmonCut; |
---|
| 806 | |
---|
| 807 | G4double ratio; |
---|
| 808 | if (dNdxCut > 0.) ratio = dNdxPhotonCut/dNdxCut; |
---|
| 809 | else ratio = 0.; |
---|
| 810 | |
---|
| 811 | if(ratio < G4UniformRand() ) // secondary e- |
---|
| 812 | { |
---|
| 813 | G4double deltaTkin = GetPostStepTransfer(fPAIplasmonTable, fdNdxCutPlasmonVector, |
---|
| 814 | iPlace, scaledTkin); |
---|
| 815 | |
---|
| 816 | // G4cout<<"PAIPhotonModel PlasmonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl ; |
---|
| 817 | |
---|
| 818 | if( deltaTkin <= 0. ) |
---|
| 819 | { |
---|
| 820 | G4cout<<"G4PAIPhotonModel::SampleSecondary e- deltaTkin = "<<deltaTkin<<G4endl; |
---|
| 821 | } |
---|
| 822 | if( deltaTkin <= 0.) return; |
---|
| 823 | |
---|
| 824 | G4double deltaTotalMomentum = sqrt(deltaTkin*(deltaTkin + 2. * electron_mass_c2 )); |
---|
| 825 | G4double totalMomentum = sqrt(pSquare); |
---|
| 826 | G4double costheta = deltaTkin*(totalEnergy + electron_mass_c2) |
---|
| 827 | /(deltaTotalMomentum * totalMomentum); |
---|
| 828 | |
---|
| 829 | if( costheta > 0.99999 ) costheta = 0.99999; |
---|
| 830 | G4double sintheta = 0.0; |
---|
| 831 | G4double sin2 = 1. - costheta*costheta; |
---|
| 832 | if( sin2 > 0.) sintheta = sqrt(sin2); |
---|
| 833 | |
---|
| 834 | // direction of the delta electron |
---|
| 835 | |
---|
| 836 | G4double phi = twopi*G4UniformRand(); |
---|
| 837 | G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta; |
---|
| 838 | |
---|
| 839 | G4ThreeVector deltaDirection(dirx,diry,dirz); |
---|
| 840 | deltaDirection.rotateUz(direction); |
---|
| 841 | |
---|
| 842 | // primary change |
---|
| 843 | |
---|
| 844 | kineticEnergy -= deltaTkin; |
---|
| 845 | G4ThreeVector dir = totalMomentum*direction - deltaTotalMomentum*deltaDirection; |
---|
| 846 | direction = dir.unit(); |
---|
| 847 | fParticleChange->SetProposedMomentumDirection(direction); |
---|
| 848 | |
---|
| 849 | // create G4DynamicParticle object for e- delta ray |
---|
| 850 | |
---|
| 851 | G4DynamicParticle* deltaRay = new G4DynamicParticle; |
---|
| 852 | deltaRay->SetDefinition(G4Electron::Electron()); |
---|
| 853 | deltaRay->SetKineticEnergy( deltaTkin ); |
---|
| 854 | deltaRay->SetMomentumDirection(deltaDirection); |
---|
| 855 | vdp->push_back(deltaRay); |
---|
| 856 | |
---|
| 857 | } |
---|
| 858 | else // secondary 'Cherenkov' photon |
---|
| 859 | { |
---|
| 860 | G4double deltaTkin = GetPostStepTransfer(fPAIphotonTable, fdNdxCutPhotonVector, |
---|
| 861 | iPlace,scaledTkin); |
---|
| 862 | |
---|
| 863 | // G4cout<<"PAIPhotonModel PhotonPostStepTransfer = "<<deltaTkin/keV<<" keV"<<G4endl ; |
---|
| 864 | |
---|
| 865 | if( deltaTkin <= 0. ) |
---|
| 866 | { |
---|
| 867 | G4cout<<"G4PAIPhotonModel::SampleSecondary gamma deltaTkin = "<<deltaTkin<<G4endl; |
---|
| 868 | } |
---|
| 869 | if( deltaTkin <= 0.) return; |
---|
| 870 | |
---|
| 871 | G4double costheta = 0.; // G4UniformRand(); // VG: ??? for start only |
---|
| 872 | G4double sintheta = sqrt((1.+costheta)*(1.-costheta)); |
---|
| 873 | |
---|
| 874 | // direction of the 'Cherenkov' photon |
---|
| 875 | G4double phi = twopi*G4UniformRand(); |
---|
| 876 | G4double dirx = sintheta*cos(phi), diry = sintheta*sin(phi), dirz = costheta; |
---|
| 877 | |
---|
| 878 | G4ThreeVector deltaDirection(dirx,diry,dirz); |
---|
| 879 | deltaDirection.rotateUz(direction); |
---|
| 880 | |
---|
| 881 | // primary change |
---|
| 882 | kineticEnergy -= deltaTkin; |
---|
| 883 | |
---|
| 884 | // create G4DynamicParticle object for photon ray |
---|
| 885 | |
---|
| 886 | G4DynamicParticle* photonRay = new G4DynamicParticle; |
---|
| 887 | photonRay->SetDefinition( G4Gamma::Gamma() ); |
---|
| 888 | photonRay->SetKineticEnergy( deltaTkin ); |
---|
| 889 | photonRay->SetMomentumDirection(deltaDirection); |
---|
| 890 | |
---|
| 891 | vdp->push_back(photonRay); |
---|
| 892 | } |
---|
| 893 | |
---|
| 894 | fParticleChange->SetProposedKineticEnergy(kineticEnergy); |
---|
| 895 | } |
---|
| 896 | |
---|
| 897 | |
---|
| 898 | /////////////////////////////////////////////////////////////////////// |
---|
| 899 | // |
---|
| 900 | // Returns post step PAI energy transfer > cut electron/photon energy according to passed |
---|
| 901 | // scaled kinetic energy of particle |
---|
| 902 | |
---|
| 903 | G4double |
---|
| 904 | G4PAIPhotonModel::GetPostStepTransfer( G4PhysicsTable* pTable, |
---|
| 905 | G4PhysicsLogVector* pVector, |
---|
| 906 | G4int iPlace, G4double scaledTkin ) |
---|
| 907 | { |
---|
| 908 | // G4cout<<"G4PAIPhotonModel::GetPostStepTransfer"<<G4endl ; |
---|
| 909 | |
---|
| 910 | G4int iTkin = iPlace+1, iTransfer; |
---|
| 911 | G4double transfer = 0.0, position, dNdxCut1, dNdxCut2, E1, E2, W1, W2, W ; |
---|
| 912 | |
---|
| 913 | dNdxCut1 = (*pVector)(iPlace) ; |
---|
| 914 | |
---|
| 915 | // G4cout<<"iPlace = "<<iPlace<<endl ; |
---|
| 916 | |
---|
| 917 | if(iTkin == fTotBin) // Fermi plato, try from left |
---|
| 918 | { |
---|
| 919 | position = dNdxCut1*G4UniformRand() ; |
---|
| 920 | |
---|
| 921 | for( iTransfer = 0; |
---|
| 922 | iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ ) |
---|
| 923 | { |
---|
| 924 | if(position >= (*(*pTable)(iPlace))(iTransfer)) break ; |
---|
| 925 | } |
---|
| 926 | transfer = GetEnergyTransfer(pTable,iPlace,position,iTransfer); |
---|
| 927 | } |
---|
| 928 | else |
---|
| 929 | { |
---|
| 930 | dNdxCut2 = (*pVector)(iPlace+1) ; |
---|
| 931 | if(iTkin == 0) // Tkin is too small, trying from right only |
---|
| 932 | { |
---|
| 933 | position = dNdxCut2*G4UniformRand() ; |
---|
| 934 | |
---|
| 935 | for( iTransfer = 0; |
---|
| 936 | iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ ) |
---|
| 937 | { |
---|
| 938 | if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break ; |
---|
| 939 | } |
---|
| 940 | transfer = GetEnergyTransfer(pTable,iPlace+1,position,iTransfer); |
---|
| 941 | } |
---|
| 942 | else // general case: Tkin between two vectors of the material |
---|
| 943 | { |
---|
| 944 | E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ; |
---|
| 945 | E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ; |
---|
| 946 | W = 1.0/(E2 - E1) ; |
---|
| 947 | W1 = (E2 - scaledTkin)*W ; |
---|
| 948 | W2 = (scaledTkin - E1)*W ; |
---|
| 949 | |
---|
| 950 | position = ( dNdxCut1*W1 + dNdxCut2*W2 )*G4UniformRand() ; |
---|
| 951 | |
---|
| 952 | // G4cout<<position<<"\t" ; |
---|
| 953 | |
---|
| 954 | G4int iTrMax1, iTrMax2, iTrMax; |
---|
| 955 | |
---|
| 956 | iTrMax1 = G4int((*pTable)(iPlace)->GetVectorLength()); |
---|
| 957 | iTrMax2 = G4int((*pTable)(iPlace+1)->GetVectorLength()); |
---|
| 958 | |
---|
| 959 | if (iTrMax1 >= iTrMax2) iTrMax = iTrMax2; |
---|
| 960 | else iTrMax = iTrMax1; |
---|
| 961 | |
---|
| 962 | for( iTransfer = 0; iTransfer < iTrMax; iTransfer++ ) |
---|
| 963 | { |
---|
| 964 | if( position >= |
---|
| 965 | ( (*(*pTable)(iPlace))(iTransfer)*W1 + |
---|
| 966 | (*(*pTable)(iPlace+1))(iTransfer)*W2) ) break ; |
---|
| 967 | } |
---|
| 968 | transfer = GetEnergyTransfer(pTable, iPlace, position, iTransfer); |
---|
| 969 | } |
---|
| 970 | } |
---|
| 971 | // G4cout<<"PAIPhotonModel PostStepTransfer = "<<transfer/keV<<" keV"<<G4endl ; |
---|
| 972 | if(transfer < 0.0 ) transfer = 0.0 ; |
---|
| 973 | return transfer ; |
---|
| 974 | } |
---|
| 975 | |
---|
| 976 | /////////////////////////////////////////////////////////////////////// |
---|
| 977 | // |
---|
| 978 | // Returns random PAI energy transfer according to passed |
---|
| 979 | // indexes of particle |
---|
| 980 | |
---|
| 981 | G4double |
---|
| 982 | G4PAIPhotonModel::GetEnergyTransfer( G4PhysicsTable* pTable, G4int iPlace, |
---|
| 983 | G4double position, G4int iTransfer ) |
---|
| 984 | { |
---|
| 985 | G4int iTransferMax; |
---|
| 986 | G4double x1, x2, y1, y2, energyTransfer; |
---|
| 987 | |
---|
| 988 | if(iTransfer == 0) |
---|
| 989 | { |
---|
| 990 | energyTransfer = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer); |
---|
| 991 | } |
---|
| 992 | else |
---|
| 993 | { |
---|
| 994 | iTransferMax = G4int((*pTable)(iPlace)->GetVectorLength()); |
---|
| 995 | |
---|
| 996 | if ( iTransfer >= iTransferMax) iTransfer = iTransferMax - 1; |
---|
| 997 | |
---|
| 998 | y1 = (*(*pTable)(iPlace))(iTransfer-1); |
---|
| 999 | y2 = (*(*fPAItransferTable)(iPlace))(iTransfer); |
---|
| 1000 | |
---|
| 1001 | x1 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer-1); |
---|
| 1002 | x2 = (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer); |
---|
| 1003 | |
---|
| 1004 | if ( x1 == x2 ) energyTransfer = x2; |
---|
| 1005 | else |
---|
| 1006 | { |
---|
| 1007 | if ( y1 == y2 ) energyTransfer = x1 + (x2 - x1)*G4UniformRand(); |
---|
| 1008 | else |
---|
| 1009 | { |
---|
| 1010 | energyTransfer = x1 + (position - y1)*(x2 - x1)/(y2 - y1); |
---|
| 1011 | } |
---|
| 1012 | } |
---|
| 1013 | } |
---|
| 1014 | return energyTransfer; |
---|
| 1015 | } |
---|
| 1016 | |
---|
| 1017 | /////////////////////////////////////////////////////////////////////// |
---|
| 1018 | // |
---|
| 1019 | // Works like AlongStepDoIt method of process family |
---|
| 1020 | |
---|
| 1021 | |
---|
| 1022 | |
---|
| 1023 | |
---|
| 1024 | G4double G4PAIPhotonModel::SampleFluctuations( const G4Material* material, |
---|
| 1025 | const G4DynamicParticle* aParticle, |
---|
| 1026 | G4double&, |
---|
| 1027 | G4double& step, |
---|
| 1028 | G4double&) |
---|
| 1029 | { |
---|
| 1030 | size_t jMat; |
---|
| 1031 | for( jMat = 0 ;jMat < fMaterialCutsCoupleVector.size() ; ++jMat ) |
---|
| 1032 | { |
---|
| 1033 | if( material == fMaterialCutsCoupleVector[jMat]->GetMaterial() ) break; |
---|
| 1034 | } |
---|
| 1035 | if(jMat == fMaterialCutsCoupleVector.size() && jMat > 0) jMat--; |
---|
| 1036 | |
---|
| 1037 | fPAItransferTable = fPAIxscBank[jMat]; |
---|
| 1038 | fPAIphotonTable = fPAIphotonBank[jMat]; |
---|
| 1039 | fPAIplasmonTable = fPAIplasmonBank[jMat]; |
---|
| 1040 | |
---|
| 1041 | fdNdxCutVector = fdNdxCutTable[jMat]; |
---|
| 1042 | fdNdxCutPhotonVector = fdNdxCutPhotonTable[jMat]; |
---|
| 1043 | fdNdxCutPlasmonVector = fdNdxCutPlasmonTable[jMat]; |
---|
| 1044 | |
---|
| 1045 | G4int iTkin, iPlace ; |
---|
| 1046 | |
---|
| 1047 | // G4cout<<"G4PAIPhotonModel::SampleFluctuations"<<G4endl ; |
---|
| 1048 | |
---|
| 1049 | G4double loss, photonLoss, plasmonLoss, charge2 ; |
---|
| 1050 | |
---|
| 1051 | |
---|
| 1052 | G4double Tkin = aParticle->GetKineticEnergy() ; |
---|
| 1053 | G4double MassRatio = proton_mass_c2/aParticle->GetDefinition()->GetPDGMass() ; |
---|
| 1054 | G4double charge = aParticle->GetDefinition()->GetPDGCharge() ; |
---|
| 1055 | charge2 = charge*charge ; |
---|
| 1056 | G4double scaledTkin = Tkin*MassRatio ; |
---|
| 1057 | G4double cof = step*charge2; |
---|
| 1058 | |
---|
[1196] | 1059 | for( iTkin = 0; iTkin <= fTotBin; iTkin++) |
---|
[819] | 1060 | { |
---|
| 1061 | if(scaledTkin < fProtonEnergyVector->GetLowEdgeEnergy(iTkin)) break ; |
---|
| 1062 | } |
---|
| 1063 | iPlace = iTkin - 1 ; |
---|
| 1064 | if( iPlace < 0 ) iPlace = 0; |
---|
| 1065 | |
---|
| 1066 | photonLoss = GetAlongStepTransfer(fPAIphotonTable,fdNdxCutPhotonVector, |
---|
| 1067 | iPlace,scaledTkin,step,cof); |
---|
| 1068 | |
---|
| 1069 | // G4cout<<"PAIPhotonModel AlongStepPhotonLoss = "<<photonLoss/keV<<" keV"<<G4endl ; |
---|
| 1070 | |
---|
| 1071 | plasmonLoss = GetAlongStepTransfer(fPAIplasmonTable,fdNdxCutPlasmonVector, |
---|
| 1072 | iPlace,scaledTkin,step,cof); |
---|
| 1073 | |
---|
| 1074 | // G4cout<<"PAIPhotonModel AlongStepPlasmonLoss = "<<plasmonLoss/keV<<" keV"<<G4endl ; |
---|
| 1075 | |
---|
| 1076 | loss = photonLoss + plasmonLoss; |
---|
| 1077 | |
---|
| 1078 | // G4cout<<"PAIPhotonModel AlongStepLoss = "<<loss/keV<<" keV"<<G4endl ; |
---|
| 1079 | |
---|
| 1080 | return loss; |
---|
| 1081 | } |
---|
| 1082 | |
---|
| 1083 | /////////////////////////////////////////////////////////////////////// |
---|
| 1084 | // |
---|
| 1085 | // Returns along step PAI energy transfer < cut electron/photon energy according to passed |
---|
| 1086 | // scaled kinetic energy of particle and cof = step*charge*charge |
---|
| 1087 | |
---|
| 1088 | G4double |
---|
| 1089 | G4PAIPhotonModel::GetAlongStepTransfer( G4PhysicsTable* pTable, |
---|
| 1090 | G4PhysicsLogVector* pVector, |
---|
| 1091 | G4int iPlace, G4double scaledTkin,G4double step, |
---|
| 1092 | G4double cof ) |
---|
| 1093 | { |
---|
| 1094 | G4int iTkin = iPlace + 1, iTransfer; |
---|
| 1095 | G4double loss = 0., position, E1, E2, W1, W2, W, dNdxCut1, dNdxCut2, meanNumber; |
---|
| 1096 | G4double lambda, stepDelta, stepSum=0. ; |
---|
| 1097 | G4long numOfCollisions=0; |
---|
| 1098 | G4bool numb = true; |
---|
| 1099 | |
---|
| 1100 | dNdxCut1 = (*pVector)(iPlace) ; |
---|
| 1101 | |
---|
| 1102 | // G4cout<<"iPlace = "<<iPlace<<endl ; |
---|
| 1103 | |
---|
| 1104 | if(iTkin == fTotBin) // Fermi plato, try from left |
---|
| 1105 | { |
---|
| 1106 | meanNumber = ((*(*pTable)(iPlace))(0) - dNdxCut1)*cof; |
---|
| 1107 | if(meanNumber < 0.) meanNumber = 0. ; |
---|
| 1108 | // numOfCollisions = RandPoisson::shoot(meanNumber) ; |
---|
| 1109 | if( meanNumber > 0.) lambda = step/meanNumber; |
---|
| 1110 | else lambda = DBL_MAX; |
---|
| 1111 | while(numb) |
---|
| 1112 | { |
---|
| 1113 | stepDelta = CLHEP::RandExponential::shoot(lambda); |
---|
| 1114 | stepSum += stepDelta; |
---|
| 1115 | if(stepSum >= step) break; |
---|
| 1116 | numOfCollisions++; |
---|
| 1117 | } |
---|
| 1118 | |
---|
| 1119 | // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ; |
---|
| 1120 | |
---|
| 1121 | while(numOfCollisions) |
---|
| 1122 | { |
---|
| 1123 | position = dNdxCut1+ |
---|
| 1124 | ((*(*pTable)(iPlace))(0) - dNdxCut1)*G4UniformRand() ; |
---|
| 1125 | |
---|
| 1126 | for( iTransfer = 0; |
---|
| 1127 | iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ ) |
---|
| 1128 | { |
---|
| 1129 | if(position >= (*(*pTable)(iPlace))(iTransfer)) break ; |
---|
| 1130 | } |
---|
| 1131 | loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer); |
---|
| 1132 | numOfCollisions-- ; |
---|
| 1133 | } |
---|
| 1134 | } |
---|
| 1135 | else |
---|
| 1136 | { |
---|
| 1137 | dNdxCut2 = (*pVector)(iPlace+1) ; |
---|
| 1138 | |
---|
| 1139 | if(iTkin == 0) // Tkin is too small, trying from right only |
---|
| 1140 | { |
---|
| 1141 | meanNumber = ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*cof; |
---|
| 1142 | if( meanNumber < 0. ) meanNumber = 0. ; |
---|
| 1143 | // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ; |
---|
| 1144 | if( meanNumber > 0.) lambda = step/meanNumber; |
---|
| 1145 | else lambda = DBL_MAX; |
---|
| 1146 | while(numb) |
---|
| 1147 | { |
---|
| 1148 | stepDelta = CLHEP::RandExponential::shoot(lambda); |
---|
| 1149 | stepSum += stepDelta; |
---|
| 1150 | if(stepSum >= step) break; |
---|
| 1151 | numOfCollisions++; |
---|
| 1152 | } |
---|
| 1153 | |
---|
| 1154 | // G4cout<<"numOfCollisions = "<<numOfCollisions<<G4endl ; |
---|
| 1155 | |
---|
| 1156 | while(numOfCollisions) |
---|
| 1157 | { |
---|
| 1158 | position = dNdxCut2+ |
---|
| 1159 | ((*(*pTable)(iPlace+1))(0) - dNdxCut2)*G4UniformRand(); |
---|
| 1160 | |
---|
| 1161 | for( iTransfer = 0; |
---|
| 1162 | iTransfer < G4int((*pTable)(iPlace+1)->GetVectorLength()); iTransfer++ ) |
---|
| 1163 | { |
---|
| 1164 | if(position >= (*(*pTable)(iPlace+1))(iTransfer)) break ; |
---|
| 1165 | } |
---|
| 1166 | loss += GetEnergyTransfer(pTable,iPlace+1,position,iTransfer); |
---|
| 1167 | numOfCollisions-- ; |
---|
| 1168 | } |
---|
| 1169 | } |
---|
| 1170 | else // general case: Tkin between two vectors of the material |
---|
| 1171 | { |
---|
| 1172 | E1 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin - 1) ; |
---|
| 1173 | E2 = fProtonEnergyVector->GetLowEdgeEnergy(iTkin) ; |
---|
| 1174 | W = 1.0/(E2 - E1) ; |
---|
| 1175 | W1 = (E2 - scaledTkin)*W ; |
---|
| 1176 | W2 = (scaledTkin - E1)*W ; |
---|
| 1177 | |
---|
| 1178 | // G4cout<<"(*(*pTable)(iPlace))(0) = "<< |
---|
| 1179 | // (*(*pTable)(iPlace))(0)<<G4endl ; |
---|
| 1180 | // G4cout<<"(*(*pTable)(iPlace+1))(0) = "<< |
---|
| 1181 | // (*(*pTable)(iPlace+1))(0)<<G4endl ; |
---|
| 1182 | |
---|
| 1183 | meanNumber=( ((*(*pTable)(iPlace))(0)-dNdxCut1)*W1 + |
---|
| 1184 | ((*(*pTable)(iPlace+1))(0)-dNdxCut2)*W2 )*cof; |
---|
| 1185 | if(meanNumber<0.0) meanNumber = 0.0; |
---|
| 1186 | // numOfCollisions = CLHEP::RandPoisson::shoot(meanNumber) ; |
---|
| 1187 | if( meanNumber > 0.) lambda = step/meanNumber; |
---|
| 1188 | else lambda = DBL_MAX; |
---|
| 1189 | while(numb) |
---|
| 1190 | { |
---|
| 1191 | stepDelta = CLHEP::RandExponential::shoot(lambda); |
---|
| 1192 | stepSum += stepDelta; |
---|
| 1193 | if(stepSum >= step) break; |
---|
| 1194 | numOfCollisions++; |
---|
| 1195 | } |
---|
| 1196 | |
---|
| 1197 | // G4cout<<"numOfCollisions = "<<numOfCollisions<<endl ; |
---|
| 1198 | |
---|
| 1199 | while(numOfCollisions) |
---|
| 1200 | { |
---|
| 1201 | position = dNdxCut1*W1 + dNdxCut2*W2 + |
---|
| 1202 | ( ( (*(*pTable)(iPlace ))(0) - dNdxCut1)*W1 + |
---|
| 1203 | |
---|
| 1204 | ( (*(*pTable)(iPlace+1))(0) - dNdxCut2)*W2 )*G4UniformRand(); |
---|
| 1205 | |
---|
| 1206 | // G4cout<<position<<"\t" ; |
---|
| 1207 | |
---|
| 1208 | for( iTransfer = 0; |
---|
| 1209 | iTransfer < G4int((*pTable)(iPlace)->GetVectorLength()); iTransfer++ ) |
---|
| 1210 | { |
---|
| 1211 | if( position >= |
---|
| 1212 | ( (*(*pTable)(iPlace))(iTransfer)*W1 + |
---|
| 1213 | (*(*pTable)(iPlace+1))(iTransfer)*W2) ) |
---|
| 1214 | { |
---|
| 1215 | break ; |
---|
| 1216 | } |
---|
| 1217 | } |
---|
| 1218 | // loss += (*pTable)(iPlace)->GetLowEdgeEnergy(iTransfer) ; |
---|
| 1219 | loss += GetEnergyTransfer(pTable,iPlace,position,iTransfer); |
---|
| 1220 | numOfCollisions-- ; |
---|
| 1221 | } |
---|
| 1222 | } |
---|
| 1223 | } |
---|
| 1224 | |
---|
| 1225 | return loss ; |
---|
| 1226 | |
---|
| 1227 | } |
---|
| 1228 | |
---|
| 1229 | ////////////////////////////////////////////////////////////////////// |
---|
| 1230 | // |
---|
| 1231 | // Returns the statistical estimation of the energy loss distribution variance |
---|
| 1232 | // |
---|
| 1233 | |
---|
| 1234 | |
---|
| 1235 | G4double G4PAIPhotonModel::Dispersion( const G4Material* material, |
---|
| 1236 | const G4DynamicParticle* aParticle, |
---|
| 1237 | G4double& tmax, |
---|
| 1238 | G4double& step ) |
---|
| 1239 | { |
---|
| 1240 | G4double loss, sumLoss=0., sumLoss2=0., sigma2, meanLoss=0.; |
---|
| 1241 | for(G4int i = 0 ; i < fMeanNumber; i++) |
---|
| 1242 | { |
---|
| 1243 | loss = SampleFluctuations(material,aParticle,tmax,step,meanLoss); |
---|
| 1244 | sumLoss += loss; |
---|
| 1245 | sumLoss2 += loss*loss; |
---|
| 1246 | } |
---|
| 1247 | meanLoss = sumLoss/fMeanNumber; |
---|
| 1248 | sigma2 = meanLoss*meanLoss + (sumLoss2-2*sumLoss*meanLoss)/fMeanNumber; |
---|
| 1249 | return sigma2; |
---|
| 1250 | } |
---|
| 1251 | |
---|
[1055] | 1252 | ///////////////////////////////////////////////////////////////////// |
---|
[819] | 1253 | |
---|
[1055] | 1254 | G4double G4PAIPhotonModel::MaxSecondaryEnergy( const G4ParticleDefinition* p, |
---|
| 1255 | G4double kinEnergy) |
---|
| 1256 | { |
---|
| 1257 | G4double tmax = kinEnergy; |
---|
| 1258 | if(p == fElectron) tmax *= 0.5; |
---|
| 1259 | else if(p != fPositron) { |
---|
| 1260 | G4double mass = p->GetPDGMass(); |
---|
| 1261 | G4double ratio= electron_mass_c2/mass; |
---|
| 1262 | G4double gamma= kinEnergy/mass + 1.0; |
---|
| 1263 | tmax = 2.0*electron_mass_c2*(gamma*gamma - 1.) / |
---|
| 1264 | (1. + 2.0*gamma*ratio + ratio*ratio); |
---|
| 1265 | } |
---|
| 1266 | return tmax; |
---|
| 1267 | } |
---|
| 1268 | |
---|
| 1269 | /////////////////////////////////////////////////////////////// |
---|
| 1270 | |
---|
| 1271 | void G4PAIPhotonModel::DefineForRegion(const G4Region* r) |
---|
| 1272 | { |
---|
| 1273 | fPAIRegionVector.push_back(r); |
---|
| 1274 | } |
---|
| 1275 | |
---|
| 1276 | |
---|
[819] | 1277 | // |
---|
| 1278 | // |
---|
| 1279 | ///////////////////////////////////////////////// |
---|
| 1280 | |
---|
| 1281 | |
---|
| 1282 | |
---|
| 1283 | |
---|
| 1284 | |
---|
| 1285 | |
---|