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