[1316] | 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 | // $Id: G4Penelope08RayleighModel.cc,v 1.1 2010/03/17 14:18:50 pandola Exp $ |
---|
| 27 | // GEANT4 tag $Name: geant4-09-04-beta-cand-01 $ |
---|
| 28 | // |
---|
| 29 | // Author: Luciano Pandola |
---|
| 30 | // |
---|
| 31 | // History: |
---|
| 32 | // -------- |
---|
| 33 | // 03 Dec 2009 L Pandola First implementation |
---|
| 34 | // |
---|
| 35 | |
---|
| 36 | #include "G4Penelope08RayleighModel.hh" |
---|
| 37 | #include "G4PenelopeSamplingData.hh" |
---|
| 38 | #include "G4ParticleDefinition.hh" |
---|
| 39 | #include "G4MaterialCutsCouple.hh" |
---|
| 40 | #include "G4ProductionCutsTable.hh" |
---|
| 41 | #include "G4DynamicParticle.hh" |
---|
| 42 | #include "G4PhysicsTable.hh" |
---|
| 43 | #include "G4ElementTable.hh" |
---|
| 44 | #include "G4Element.hh" |
---|
| 45 | #include "G4PenelopeIntegrator.hh" |
---|
| 46 | #include "G4PhysicsFreeVector.hh" |
---|
| 47 | |
---|
| 48 | |
---|
| 49 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 50 | |
---|
| 51 | |
---|
| 52 | G4Penelope08RayleighModel::G4Penelope08RayleighModel(const G4ParticleDefinition*, |
---|
| 53 | const G4String& nam) |
---|
| 54 | :G4VEmModel(nam),isInitialised(false),logAtomicCrossSection(0), |
---|
| 55 | atomicFormFactor(0),logFormFactorTable(0),pMaxTable(0),samplingTable(0) |
---|
| 56 | { |
---|
| 57 | fIntrinsicLowEnergyLimit = 100.0*eV; |
---|
| 58 | fIntrinsicHighEnergyLimit = 100.0*GeV; |
---|
| 59 | // SetLowEnergyLimit(fIntrinsicLowEnergyLimit); |
---|
| 60 | SetHighEnergyLimit(fIntrinsicHighEnergyLimit); |
---|
| 61 | // |
---|
| 62 | verboseLevel= 0; |
---|
| 63 | // Verbosity scale: |
---|
| 64 | // 0 = nothing |
---|
| 65 | // 1 = warning for energy non-conservation |
---|
| 66 | // 2 = details of energy budget |
---|
| 67 | // 3 = calculation of cross sections, file openings, sampling of atoms |
---|
| 68 | // 4 = entering in methods |
---|
| 69 | |
---|
| 70 | //build the energy grid. It is the same for all materials |
---|
| 71 | G4double logenergy = log(fIntrinsicLowEnergyLimit/2.); |
---|
| 72 | G4double logmaxenergy = log(1.5*fIntrinsicHighEnergyLimit); |
---|
| 73 | //finer grid below 160 keV |
---|
| 74 | G4double logtransitionenergy = log(160*keV); |
---|
| 75 | G4double logfactor1 = log(10.)/250.; |
---|
| 76 | G4double logfactor2 = logfactor1*10; |
---|
| 77 | logEnergyGridPMax.push_back(logenergy); |
---|
| 78 | do{ |
---|
| 79 | if (logenergy < logtransitionenergy) |
---|
| 80 | logenergy += logfactor1; |
---|
| 81 | else |
---|
| 82 | logenergy += logfactor2; |
---|
| 83 | logEnergyGridPMax.push_back(logenergy); |
---|
| 84 | }while (logenergy < logmaxenergy); |
---|
| 85 | } |
---|
| 86 | |
---|
| 87 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 88 | |
---|
| 89 | G4Penelope08RayleighModel::~G4Penelope08RayleighModel() |
---|
| 90 | { |
---|
| 91 | std::map <const G4int,G4PhysicsFreeVector*>::iterator i; |
---|
| 92 | if (logAtomicCrossSection) |
---|
| 93 | { |
---|
| 94 | for (i=logAtomicCrossSection->begin();i != logAtomicCrossSection->end();i++) |
---|
| 95 | if (i->second) delete i->second; |
---|
| 96 | delete logAtomicCrossSection; |
---|
| 97 | } |
---|
| 98 | |
---|
| 99 | if (atomicFormFactor) |
---|
| 100 | { |
---|
| 101 | for (i=atomicFormFactor->begin();i != atomicFormFactor->end();i++) |
---|
| 102 | if (i->second) delete i->second; |
---|
| 103 | delete atomicFormFactor; |
---|
| 104 | } |
---|
| 105 | |
---|
| 106 | ClearTables(); |
---|
| 107 | } |
---|
| 108 | |
---|
| 109 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 110 | void G4Penelope08RayleighModel::ClearTables() |
---|
| 111 | { |
---|
| 112 | std::map <const G4Material*,G4PhysicsFreeVector*>::iterator i; |
---|
| 113 | |
---|
| 114 | if (logFormFactorTable) |
---|
| 115 | { |
---|
| 116 | for (i=logFormFactorTable->begin(); i != logFormFactorTable->end(); i++) |
---|
| 117 | if (i->second) delete i->second; |
---|
| 118 | delete logFormFactorTable; |
---|
| 119 | logFormFactorTable = 0; //zero explicitely |
---|
| 120 | } |
---|
| 121 | |
---|
| 122 | if (pMaxTable) |
---|
| 123 | { |
---|
| 124 | for (i=pMaxTable->begin(); i != pMaxTable->end(); i++) |
---|
| 125 | if (i->second) delete i->second; |
---|
| 126 | delete pMaxTable; |
---|
| 127 | pMaxTable = 0; //zero explicitely |
---|
| 128 | } |
---|
| 129 | |
---|
| 130 | std::map<const G4Material*,G4PenelopeSamplingData*>::iterator ii; |
---|
| 131 | if (samplingTable) |
---|
| 132 | { |
---|
| 133 | for (ii=samplingTable->begin(); ii != samplingTable->end(); ii++) |
---|
| 134 | if (ii->second) delete ii->second; |
---|
| 135 | delete samplingTable; |
---|
| 136 | samplingTable = 0; //zero explicitely |
---|
| 137 | } |
---|
| 138 | |
---|
| 139 | return; |
---|
| 140 | } |
---|
| 141 | |
---|
| 142 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 143 | |
---|
| 144 | void G4Penelope08RayleighModel::Initialise(const G4ParticleDefinition* , |
---|
| 145 | const G4DataVector& ) |
---|
| 146 | { |
---|
| 147 | if (verboseLevel > 3) |
---|
| 148 | G4cout << "Calling G4Penelope08RayleighModel::Initialise()" << G4endl; |
---|
| 149 | |
---|
| 150 | //clear tables depending on materials, not the atomic ones |
---|
| 151 | ClearTables(); |
---|
| 152 | |
---|
| 153 | //create new tables |
---|
| 154 | // |
---|
| 155 | // logAtomicCrossSection and atomicFormFactor are created only once, |
---|
| 156 | // since they are never cleared |
---|
| 157 | if (!logAtomicCrossSection) |
---|
| 158 | logAtomicCrossSection = new std::map<const G4int,G4PhysicsFreeVector*>; |
---|
| 159 | if (!atomicFormFactor) |
---|
| 160 | atomicFormFactor = new std::map<const G4int,G4PhysicsFreeVector*>; |
---|
| 161 | |
---|
| 162 | if (!logFormFactorTable) |
---|
| 163 | logFormFactorTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; |
---|
| 164 | if (!pMaxTable) |
---|
| 165 | pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; |
---|
| 166 | if (!samplingTable) |
---|
| 167 | samplingTable = new std::map<const G4Material*,G4PenelopeSamplingData*>; |
---|
| 168 | |
---|
| 169 | |
---|
| 170 | if (verboseLevel > 0) { |
---|
| 171 | G4cout << "Penelope08 Rayleigh model is initialized " << G4endl |
---|
| 172 | << "Energy range: " |
---|
| 173 | << LowEnergyLimit() / keV << " keV - " |
---|
| 174 | << HighEnergyLimit() / GeV << " GeV" |
---|
| 175 | << G4endl; |
---|
| 176 | } |
---|
| 177 | |
---|
| 178 | if(isInitialised) return; |
---|
| 179 | fParticleChange = GetParticleChangeForGamma(); |
---|
| 180 | isInitialised = true; |
---|
| 181 | } |
---|
| 182 | |
---|
| 183 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 184 | |
---|
| 185 | G4double G4Penelope08RayleighModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*, |
---|
| 186 | G4double energy, |
---|
| 187 | G4double Z, |
---|
| 188 | G4double, |
---|
| 189 | G4double, |
---|
| 190 | G4double) |
---|
| 191 | { |
---|
| 192 | // Cross section of Rayleigh scattering in Penelope2008 is calculated by the EPDL97 |
---|
| 193 | // tabulation, Cuellen et al. (1997), with non-relativistic form factors from Hubbel |
---|
| 194 | // et al. J. Phys. Chem. Ref. Data 4 (1975) 471; Erratum ibid. 6 (1977) 615. |
---|
| 195 | |
---|
| 196 | if (verboseLevel > 3) |
---|
| 197 | G4cout << "Calling CrossSectionPerAtom() of G4Penelope08RayleighModel" << G4endl; |
---|
| 198 | |
---|
| 199 | G4int iZ = (G4int) Z; |
---|
| 200 | |
---|
| 201 | //read data files |
---|
| 202 | if (!logAtomicCrossSection->count(iZ)) |
---|
| 203 | ReadDataFile(iZ); |
---|
| 204 | //now it should be ok |
---|
| 205 | if (!logAtomicCrossSection->count(iZ)) |
---|
| 206 | { |
---|
| 207 | G4cout << "Problem in G4Penelope08RayleighModel::ComputeCrossSectionPerAtom" |
---|
| 208 | << G4endl; |
---|
| 209 | G4Exception(); |
---|
| 210 | } |
---|
| 211 | |
---|
| 212 | G4double cross = 0; |
---|
| 213 | |
---|
| 214 | G4PhysicsFreeVector* atom = logAtomicCrossSection->find(iZ)->second; |
---|
| 215 | if (!atom) |
---|
| 216 | { |
---|
| 217 | G4cout << "Problem in G4Penelope08RayleighModel::ComputeCrossSectionPerAtom" |
---|
| 218 | << G4endl; |
---|
| 219 | G4Exception(); |
---|
| 220 | } |
---|
| 221 | G4double logene = log(energy); |
---|
| 222 | G4double logXS = atom->Value(logene); |
---|
| 223 | cross = exp(logXS); |
---|
| 224 | |
---|
| 225 | if (verboseLevel > 2) |
---|
| 226 | G4cout << "Compton cross Section at " << energy/keV << " keV for Z=" << Z << |
---|
| 227 | " = " << cross/barn << " barn" << G4endl; |
---|
| 228 | return cross; |
---|
| 229 | } |
---|
| 230 | |
---|
| 231 | |
---|
| 232 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 233 | void G4Penelope08RayleighModel::BuildFormFactorTable(const G4Material* material) |
---|
| 234 | { |
---|
| 235 | |
---|
| 236 | /* |
---|
| 237 | 1) get composition and equivalent molecular density |
---|
| 238 | */ |
---|
| 239 | |
---|
| 240 | G4int nElements = material->GetNumberOfElements(); |
---|
| 241 | const G4ElementVector* elementVector = material->GetElementVector(); |
---|
| 242 | const G4double* fractionVector = material->GetFractionVector(); |
---|
| 243 | |
---|
| 244 | std::vector<G4double> *StechiometricFactors = new std::vector<G4double>; |
---|
| 245 | for (G4int i=0;i<nElements;i++) |
---|
| 246 | { |
---|
| 247 | G4double fraction = fractionVector[i]; |
---|
| 248 | G4double atomicWeigth = (*elementVector)[i]->GetA()/(g/mole); |
---|
| 249 | StechiometricFactors->push_back(fraction/atomicWeigth); |
---|
| 250 | } |
---|
| 251 | //Find max |
---|
| 252 | G4double MaxStechiometricFactor = 0.; |
---|
| 253 | for (G4int i=0;i<nElements;i++) |
---|
| 254 | { |
---|
| 255 | if ((*StechiometricFactors)[i] > MaxStechiometricFactor) |
---|
| 256 | MaxStechiometricFactor = (*StechiometricFactors)[i]; |
---|
| 257 | } |
---|
| 258 | if (MaxStechiometricFactor<1e-16) |
---|
| 259 | { |
---|
| 260 | G4cout << "G4Penelope08RayleighModel::BuildFormFactorTable" << G4endl; |
---|
| 261 | G4cout << "Problem with the mass composition of " << material->GetName() << G4endl; |
---|
| 262 | G4Exception(); |
---|
| 263 | } |
---|
| 264 | //Normalize |
---|
| 265 | for (G4int i=0;i<nElements;i++) |
---|
| 266 | (*StechiometricFactors)[i] /= MaxStechiometricFactor; |
---|
| 267 | |
---|
| 268 | // Equivalent atoms per molecule |
---|
| 269 | G4double atomsPerMolecule = 0; |
---|
| 270 | for (G4int i=0;i<nElements;i++) |
---|
| 271 | atomsPerMolecule += (*StechiometricFactors)[i]; |
---|
| 272 | |
---|
| 273 | /* |
---|
| 274 | CREATE THE FORM FACTOR TABLE |
---|
| 275 | */ |
---|
| 276 | G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector(logQSquareGrid.size()); |
---|
| 277 | theFFVec->SetSpline(true); |
---|
| 278 | |
---|
| 279 | for (size_t k=0;k<logQSquareGrid.size();k++) |
---|
| 280 | { |
---|
| 281 | G4double ff2 = 0; //squared form factor |
---|
| 282 | for (G4int i=0;i<nElements;i++) |
---|
| 283 | { |
---|
| 284 | G4int iZ = (G4int) (*elementVector)[i]->GetZ(); |
---|
| 285 | G4PhysicsFreeVector* theAtomVec = atomicFormFactor->find(iZ)->second; |
---|
| 286 | G4double f = (*theAtomVec)[k]; //the q-grid is always the same |
---|
| 287 | ff2 += f*f*(*StechiometricFactors)[i]; |
---|
| 288 | } |
---|
| 289 | if (ff2) |
---|
| 290 | theFFVec->PutValue(k,logQSquareGrid[k],log(ff2)); //NOTICE: THIS IS log(Q^2) vs. log(F^2) |
---|
| 291 | } |
---|
| 292 | logFormFactorTable->insert(std::make_pair(material,theFFVec)); |
---|
| 293 | |
---|
| 294 | delete StechiometricFactors; |
---|
| 295 | |
---|
| 296 | return; |
---|
| 297 | } |
---|
| 298 | |
---|
| 299 | |
---|
| 300 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 301 | |
---|
| 302 | void G4Penelope08RayleighModel::SampleSecondaries(std::vector<G4DynamicParticle*>* , |
---|
| 303 | const G4MaterialCutsCouple* couple, |
---|
| 304 | const G4DynamicParticle* aDynamicGamma, |
---|
| 305 | G4double, |
---|
| 306 | G4double) |
---|
| 307 | { |
---|
| 308 | // Sampling of the Rayleigh final state (namely, scattering angle of the photon) |
---|
| 309 | // from the Penelope2008 model. The scattering angle is sampled from the atomic |
---|
| 310 | // cross section dOmega/d(cosTheta) from Born ("Atomic Phyisics", 1969), disregarding |
---|
| 311 | // anomalous scattering effects. The Form Factor F(Q) function which appears in the |
---|
| 312 | // analytical cross section is retrieved via the method GetFSquared(); atomic data |
---|
| 313 | // are tabulated for F(Q). Form factor for compounds is calculated according to |
---|
| 314 | // the additivity rule. The sampling from the F(Q) is made via a Rational Inverse |
---|
| 315 | // Transform with Aliasing (RITA) algorithm; RITA parameters are calculated once |
---|
| 316 | // for each material and managed by G4PenelopeSamplingData objects. |
---|
| 317 | // The sampling algorithm (rejection method) has efficiency 67% at low energy, and |
---|
| 318 | // increases with energy. For E=100 keV the efficiency is 100% and 86% for |
---|
| 319 | // hydrogen and uranium, respectively. |
---|
| 320 | |
---|
| 321 | if (verboseLevel > 3) |
---|
| 322 | G4cout << "Calling SamplingSecondaries() of G4Penelope08RayleighModel" << G4endl; |
---|
| 323 | |
---|
| 324 | G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy(); |
---|
| 325 | |
---|
| 326 | if (photonEnergy0 <= fIntrinsicLowEnergyLimit) |
---|
| 327 | { |
---|
| 328 | fParticleChange->ProposeTrackStatus(fStopAndKill); |
---|
| 329 | fParticleChange->SetProposedKineticEnergy(0.); |
---|
| 330 | fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0); |
---|
| 331 | return ; |
---|
| 332 | } |
---|
| 333 | |
---|
| 334 | G4ParticleMomentum photonDirection0 = aDynamicGamma->GetMomentumDirection(); |
---|
| 335 | |
---|
| 336 | const G4Material* theMat = couple->GetMaterial(); |
---|
| 337 | |
---|
| 338 | //1) Verify if tables are ready |
---|
| 339 | if (!pMaxTable || !samplingTable) |
---|
| 340 | { |
---|
| 341 | G4cout << " G4Penelope08RayleighModel::SampleSecondaries" << G4endl; |
---|
| 342 | G4cout << "Looks like the model is _not_ properly initialized" << G4endl; |
---|
| 343 | G4Exception(); |
---|
| 344 | } |
---|
| 345 | |
---|
| 346 | //2) retrieve or build the sampling table |
---|
| 347 | if (!(samplingTable->count(theMat))) |
---|
| 348 | InitializeSamplingAlgorithm(theMat); |
---|
| 349 | G4PenelopeSamplingData* theDataTable = samplingTable->find(theMat)->second; |
---|
| 350 | |
---|
| 351 | //3) retrieve or build the pMax data |
---|
| 352 | if (!pMaxTable->count(theMat)) |
---|
| 353 | GetPMaxTable(theMat); |
---|
| 354 | G4PhysicsFreeVector* thePMax = pMaxTable->find(theMat)->second; |
---|
| 355 | |
---|
| 356 | G4double cosTheta = 1.0; |
---|
| 357 | |
---|
| 358 | //OK, ready to go! |
---|
| 359 | G4double qmax = 2.0*photonEnergy0/electron_mass_c2; //this is non-dimensional now |
---|
| 360 | |
---|
| 361 | if (qmax < 1e-10) //very low momentum transfer |
---|
| 362 | { |
---|
| 363 | G4bool loopAgain=false; |
---|
| 364 | do |
---|
| 365 | { |
---|
| 366 | loopAgain = false; |
---|
| 367 | cosTheta = 1.0-2.0*G4UniformRand(); |
---|
| 368 | G4double G = 0.5*(1+cosTheta*cosTheta); |
---|
| 369 | if (G4UniformRand()>G) |
---|
| 370 | loopAgain = true; |
---|
| 371 | }while(loopAgain); |
---|
| 372 | } |
---|
| 373 | else //larger momentum transfer |
---|
| 374 | { |
---|
| 375 | size_t nData = theDataTable->GetNumberOfStoredPoints(); |
---|
| 376 | G4double LastQ2inTheTable = theDataTable->GetX(nData-1); |
---|
| 377 | G4double q2max = std::min(qmax*qmax,LastQ2inTheTable); |
---|
| 378 | |
---|
| 379 | G4bool loopAgain = false; |
---|
| 380 | G4double MaxPValue = thePMax->Value(photonEnergy0); |
---|
| 381 | G4double xx=0; |
---|
| 382 | |
---|
| 383 | //Sampling by rejection method. The rejection function is |
---|
| 384 | //G = 0.5*(1+cos^2(theta)) |
---|
| 385 | // |
---|
| 386 | do{ |
---|
| 387 | loopAgain = false; |
---|
| 388 | G4double RandomMax = G4UniformRand()*MaxPValue; |
---|
| 389 | xx = theDataTable->SampleValue(RandomMax); |
---|
| 390 | //xx is a random value of q^2 in (0,q2max),sampled according to |
---|
| 391 | //F(Q^2) via the RITA algorithm |
---|
| 392 | if (xx > q2max) |
---|
| 393 | loopAgain = true; |
---|
| 394 | cosTheta = 1.0-2.0*xx/q2max; |
---|
| 395 | G4double G = 0.5*(1+cosTheta*cosTheta); |
---|
| 396 | if (G4UniformRand()>G) |
---|
| 397 | loopAgain = true; |
---|
| 398 | }while(loopAgain); |
---|
| 399 | } |
---|
| 400 | |
---|
| 401 | G4double sinTheta = std::sqrt(1-cosTheta*cosTheta); |
---|
| 402 | |
---|
| 403 | // Scattered photon angles. ( Z - axis along the parent photon) |
---|
| 404 | G4double phi = twopi * G4UniformRand() ; |
---|
| 405 | G4double dirX = sinTheta*std::cos(phi); |
---|
| 406 | G4double dirY = sinTheta*std::sin(phi); |
---|
| 407 | G4double dirZ = cosTheta; |
---|
| 408 | |
---|
| 409 | // Update G4VParticleChange for the scattered photon |
---|
| 410 | G4ThreeVector photonDirection1(dirX, dirY, dirZ); |
---|
| 411 | |
---|
| 412 | photonDirection1.rotateUz(photonDirection0); |
---|
| 413 | fParticleChange->ProposeMomentumDirection(photonDirection1) ; |
---|
| 414 | fParticleChange->SetProposedKineticEnergy(photonEnergy0) ; |
---|
| 415 | |
---|
| 416 | |
---|
| 417 | return; |
---|
| 418 | } |
---|
| 419 | |
---|
| 420 | |
---|
| 421 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 422 | |
---|
| 423 | void G4Penelope08RayleighModel::ReadDataFile(const G4int Z) |
---|
| 424 | { |
---|
| 425 | if (verboseLevel > 2) |
---|
| 426 | { |
---|
| 427 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
| 428 | G4cout << "Going to read Rayleigh data files for Z=" << Z << G4endl; |
---|
| 429 | } |
---|
| 430 | |
---|
| 431 | char* path = getenv("G4LEDATA"); |
---|
| 432 | if (!path) |
---|
| 433 | { |
---|
| 434 | G4String excep = "G4Penelope08RayleighModel - G4LEDATA environment variable not set!"; |
---|
| 435 | G4Exception(excep); |
---|
| 436 | } |
---|
| 437 | |
---|
| 438 | /* |
---|
| 439 | Read first the cross section file |
---|
| 440 | */ |
---|
| 441 | std::ostringstream ost; |
---|
| 442 | if (Z>9) |
---|
| 443 | ost << path << "/penelope/rayleigh/pdgra" << Z << ".p08"; |
---|
| 444 | else |
---|
| 445 | ost << path << "/penelope/rayleigh/pdgra0" << Z << ".p08"; |
---|
| 446 | std::ifstream file(ost.str().c_str()); |
---|
| 447 | if (!file.is_open()) |
---|
| 448 | { |
---|
| 449 | G4String excep = "G4Penelope08RayleighModel - data file " + G4String(ost.str()) + " not found!"; |
---|
| 450 | G4Exception(excep); |
---|
| 451 | } |
---|
| 452 | G4int readZ =0; |
---|
| 453 | size_t nPoints= 0; |
---|
| 454 | file >> readZ >> nPoints; |
---|
| 455 | //check the right file is opened. |
---|
| 456 | if (readZ != Z || nPoints <= 0) |
---|
| 457 | { |
---|
| 458 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
| 459 | G4cout << "Corrupted data file for Z=" << Z << G4endl; |
---|
| 460 | G4Exception(); |
---|
| 461 | } |
---|
| 462 | G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector((size_t)nPoints); |
---|
| 463 | G4double ene=0,f1=0,f2=0,xs=0; |
---|
| 464 | for (size_t i=0;i<nPoints;i++) |
---|
| 465 | { |
---|
| 466 | file >> ene >> f1 >> f2 >> xs; |
---|
| 467 | //dimensional quantities |
---|
| 468 | ene *= eV; |
---|
| 469 | xs *= cm2; |
---|
| 470 | theVec->PutValue(i,log(ene),log(xs)); |
---|
| 471 | if (file.eof() && i != (nPoints-1)) //file ended too early |
---|
| 472 | { |
---|
| 473 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
| 474 | G4cout << "Corrupted data file for Z=" << Z << G4endl; |
---|
| 475 | G4cout << "Found less than " << nPoints << "entries " <<G4endl; |
---|
| 476 | G4Exception(); |
---|
| 477 | } |
---|
| 478 | } |
---|
| 479 | if (!logAtomicCrossSection) |
---|
| 480 | { |
---|
| 481 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
| 482 | G4cout << "Problem with allocation of logAtomicCrossSection data table " << G4endl; |
---|
| 483 | G4Exception(); |
---|
| 484 | } |
---|
| 485 | logAtomicCrossSection->insert(std::make_pair(Z,theVec)); |
---|
| 486 | file.close(); |
---|
| 487 | |
---|
| 488 | /* |
---|
| 489 | Then read the form factor file |
---|
| 490 | */ |
---|
| 491 | std::ostringstream ost2; |
---|
| 492 | if (Z>9) |
---|
| 493 | ost2 << path << "/penelope/rayleigh/pdaff" << Z << ".p08"; |
---|
| 494 | else |
---|
| 495 | ost2 << path << "/penelope/rayleigh/pdaff0" << Z << ".p08"; |
---|
| 496 | file.open(ost2.str().c_str()); |
---|
| 497 | if (!file.is_open()) |
---|
| 498 | { |
---|
| 499 | G4String excep = "G4Penelope08RayleighModel - data file " + G4String(ost2.str()) + " not found!"; |
---|
| 500 | G4Exception(excep); |
---|
| 501 | } |
---|
| 502 | file >> readZ >> nPoints; |
---|
| 503 | //check the right file is opened. |
---|
| 504 | if (readZ != Z || nPoints <= 0) |
---|
| 505 | { |
---|
| 506 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
| 507 | G4cout << "Corrupted data file for Z=" << Z << G4endl; |
---|
| 508 | G4Exception(); |
---|
| 509 | } |
---|
| 510 | G4PhysicsFreeVector* theFFVec = new G4PhysicsFreeVector((size_t)nPoints); |
---|
| 511 | G4double q=0,ff=0,incoh=0; |
---|
| 512 | G4bool fillQGrid = false; |
---|
| 513 | //fill this vector only the first time. |
---|
| 514 | if (!logQSquareGrid.size()) |
---|
| 515 | fillQGrid = true; |
---|
| 516 | for (size_t i=0;i<nPoints;i++) |
---|
| 517 | { |
---|
| 518 | file >> q >> ff >> incoh; |
---|
| 519 | //q and ff are dimensionless (q is in units of (m_e*c) |
---|
| 520 | theFFVec->PutValue(i,q,ff); |
---|
| 521 | if (fillQGrid) |
---|
| 522 | { |
---|
| 523 | logQSquareGrid.push_back(2.0*log(q)); |
---|
| 524 | } |
---|
| 525 | if (file.eof() && i != (nPoints-1)) //file ended too early |
---|
| 526 | { |
---|
| 527 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
| 528 | G4cout << "Corrupted data file for Z=" << Z << G4endl; |
---|
| 529 | G4cout << "Found less than " << nPoints << "entries " <<G4endl; |
---|
| 530 | G4Exception(); |
---|
| 531 | } |
---|
| 532 | } |
---|
| 533 | if (!atomicFormFactor) |
---|
| 534 | { |
---|
| 535 | G4cout << "G4Penelope08RayleighModel::ReadDataFile()" << G4endl; |
---|
| 536 | G4cout << "Problem with allocation of atomicFormFactor data table " << G4endl; |
---|
| 537 | G4Exception(); |
---|
| 538 | } |
---|
| 539 | atomicFormFactor->insert(std::make_pair(Z,theFFVec)); |
---|
| 540 | file.close(); |
---|
| 541 | return; |
---|
| 542 | } |
---|
| 543 | |
---|
| 544 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 545 | |
---|
| 546 | G4double G4Penelope08RayleighModel::GetFSquared(const G4Material* mat, const G4double QSquared) |
---|
| 547 | { |
---|
| 548 | G4double f2 = 0; |
---|
| 549 | G4double logQSquared = log(QSquared); |
---|
| 550 | //last value of the table |
---|
| 551 | G4double maxlogQ2 = logQSquareGrid[logQSquareGrid.size()-1]; |
---|
| 552 | //If the table has not been built for the material, do it! |
---|
| 553 | if (!logFormFactorTable->count(mat)) |
---|
| 554 | BuildFormFactorTable(mat); |
---|
| 555 | |
---|
| 556 | //now it should be all right |
---|
| 557 | G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second; |
---|
| 558 | |
---|
| 559 | if (!theVec) |
---|
| 560 | { |
---|
| 561 | G4cout << " G4Penelope08RayleighModel::GetFSquared()" << G4endl; |
---|
| 562 | G4cout << "Unable to retrieve F squared table for " << mat->GetName(); |
---|
| 563 | G4Exception(); |
---|
| 564 | } |
---|
| 565 | if (logQSquared < -20) // Q < 1e-9 |
---|
| 566 | { |
---|
| 567 | G4double logf2 = (*theVec)[0]; //first value of the table |
---|
| 568 | f2 = exp(logf2); |
---|
| 569 | } |
---|
| 570 | else if (logQSquared > maxlogQ2) |
---|
| 571 | f2 =0; |
---|
| 572 | else |
---|
| 573 | { |
---|
| 574 | //log(Q^2) vs. log(F^2) |
---|
| 575 | G4double logf2 = theVec->Value(logQSquared); |
---|
| 576 | f2 = exp(logf2); |
---|
| 577 | |
---|
| 578 | } |
---|
| 579 | if (verboseLevel > 3) |
---|
| 580 | { |
---|
| 581 | G4cout << "G4Penelope08RayleighModel::GetFSquared() in " << mat->GetName() << G4endl; |
---|
| 582 | G4cout << "Q^2 = " << QSquared << " (units of 1/(m_e*c); F^2 = " << f2 << G4endl; |
---|
| 583 | } |
---|
| 584 | return f2; |
---|
| 585 | } |
---|
| 586 | |
---|
| 587 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 588 | |
---|
| 589 | void G4Penelope08RayleighModel::InitializeSamplingAlgorithm(const G4Material* mat) |
---|
| 590 | { |
---|
| 591 | G4double q2min = 0; |
---|
| 592 | G4double q2max = 0; |
---|
| 593 | const size_t np = 150; //hard-coded in Penelope |
---|
| 594 | for (size_t i=1;i<logQSquareGrid.size();i++) |
---|
| 595 | { |
---|
| 596 | G4double Q2 = exp(logQSquareGrid[i]); |
---|
| 597 | if (GetFSquared(mat,Q2) > 1e-35) |
---|
| 598 | { |
---|
| 599 | q2max = exp(logQSquareGrid[i-1]); |
---|
| 600 | } |
---|
| 601 | } |
---|
| 602 | |
---|
| 603 | size_t nReducedPoints = np/4; |
---|
| 604 | |
---|
| 605 | //check for errors |
---|
| 606 | if (np < 16) |
---|
| 607 | { |
---|
| 608 | G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl; |
---|
| 609 | G4cout << "Too few points to initialize the sampling algorithm" << G4endl; |
---|
| 610 | G4Exception(); |
---|
| 611 | } |
---|
| 612 | if (q2min > (q2max-1e-10)) |
---|
| 613 | { |
---|
| 614 | G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl; |
---|
| 615 | G4cout << "Too narrow grid to initialize the sampling algorithm" << G4endl; |
---|
| 616 | G4Exception(); |
---|
| 617 | } |
---|
| 618 | |
---|
| 619 | //This is subroutine RITAI0 of Penelope |
---|
| 620 | //Create an object of type G4Penelope08RayleighSamplingData --> store in a map::Material* |
---|
| 621 | |
---|
| 622 | //temporary vectors --> Then everything is stored in G4PenelopeSamplingData |
---|
| 623 | G4DataVector* x = new G4DataVector(); |
---|
| 624 | |
---|
| 625 | /******************************************************************************* |
---|
| 626 | Start with a grid of NUNIF points uniformly spaced in the interval q2min,q2max |
---|
| 627 | ********************************************************************************/ |
---|
| 628 | size_t NUNIF = std::min(std::max(((size_t)8),nReducedPoints),np/2); |
---|
| 629 | const G4int nip = 51; //hard-coded in Penelope |
---|
| 630 | |
---|
| 631 | G4double dx = (q2max-q2min)/((G4double) NUNIF-1); |
---|
| 632 | x->push_back(q2min); |
---|
| 633 | for (size_t i=1;i<NUNIF-1;i++) |
---|
| 634 | { |
---|
| 635 | G4double app = q2min + i*dx; |
---|
| 636 | x->push_back(app); //increase |
---|
| 637 | } |
---|
| 638 | x->push_back(q2max); |
---|
| 639 | |
---|
| 640 | if (verboseLevel> 3) |
---|
| 641 | G4cout << "Vector x has " << x->size() << " points, while NUNIF = " << NUNIF << G4endl; |
---|
| 642 | |
---|
| 643 | //Allocate temporary storage vectors |
---|
| 644 | G4DataVector* area = new G4DataVector(); |
---|
| 645 | G4DataVector* a = new G4DataVector(); |
---|
| 646 | G4DataVector* b = new G4DataVector(); |
---|
| 647 | G4DataVector* c = new G4DataVector(); |
---|
| 648 | G4DataVector* err = new G4DataVector(); |
---|
| 649 | |
---|
| 650 | for (size_t i=0;i<NUNIF-1;i++) //build all points but the last |
---|
| 651 | { |
---|
| 652 | //Temporary vectors for this loop |
---|
| 653 | G4DataVector* pdfi = new G4DataVector(); |
---|
| 654 | G4DataVector* pdfih = new G4DataVector(); |
---|
| 655 | G4DataVector* sumi = new G4DataVector(); |
---|
| 656 | |
---|
| 657 | G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1)); |
---|
| 658 | G4double pdfmax = 0; |
---|
| 659 | for (G4int k=0;k<nip;k++) |
---|
| 660 | { |
---|
| 661 | G4double xik = (*x)[i]+k*dxi; |
---|
| 662 | G4double pdfk = std::max(GetFSquared(mat,xik),0.); |
---|
| 663 | pdfi->push_back(pdfk); |
---|
| 664 | pdfmax = std::max(pdfmax,pdfk); |
---|
| 665 | if (k < (nip-1)) |
---|
| 666 | { |
---|
| 667 | G4double xih = xik + 0.5*dxi; |
---|
| 668 | G4double pdfIK = std::max(GetFSquared(mat,xih),0.); |
---|
| 669 | pdfih->push_back(pdfIK); |
---|
| 670 | pdfmax = std::max(pdfmax,pdfIK); |
---|
| 671 | } |
---|
| 672 | } |
---|
| 673 | |
---|
| 674 | //Simpson's integration |
---|
| 675 | G4double cons = dxi*0.5*(1./3.); |
---|
| 676 | sumi->push_back(0.); |
---|
| 677 | for (G4int k=1;k<nip;k++) |
---|
| 678 | { |
---|
| 679 | G4double previous = (*sumi)[k-1]; |
---|
| 680 | G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]); |
---|
| 681 | sumi->push_back(next); |
---|
| 682 | } |
---|
| 683 | |
---|
| 684 | G4double lastIntegral = (*sumi)[sumi->size()-1]; |
---|
| 685 | area->push_back(lastIntegral); |
---|
| 686 | //Normalize cumulative function |
---|
| 687 | G4double factor = 1.0/lastIntegral; |
---|
| 688 | for (size_t k=0;k<sumi->size();k++) |
---|
| 689 | (*sumi)[k] *= factor; |
---|
| 690 | |
---|
| 691 | //When the PDF vanishes at one of the interval end points, its value is modified |
---|
| 692 | if ((*pdfi)[0] < 1e-35) |
---|
| 693 | (*pdfi)[0] = 1e-5*pdfmax; |
---|
| 694 | if ((*pdfi)[pdfi->size()-1] < 1e-35) |
---|
| 695 | (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; |
---|
| 696 | |
---|
| 697 | G4double pli = (*pdfi)[0]*factor; |
---|
| 698 | G4double pui = (*pdfi)[pdfi->size()-1]*factor; |
---|
| 699 | G4double B_temp = 1.0-1.0/(pli*pui*dx*dx); |
---|
| 700 | G4double A_temp = (1.0/(pli*dx))-1.0-B_temp; |
---|
| 701 | G4double C_temp = 1.0+A_temp+B_temp; |
---|
| 702 | if (C_temp < 1e-35) |
---|
| 703 | { |
---|
| 704 | a->push_back(0.); |
---|
| 705 | b->push_back(0.); |
---|
| 706 | c->push_back(1.); |
---|
| 707 | } |
---|
| 708 | else |
---|
| 709 | { |
---|
| 710 | a->push_back(A_temp); |
---|
| 711 | b->push_back(B_temp); |
---|
| 712 | c->push_back(C_temp); |
---|
| 713 | } |
---|
| 714 | |
---|
| 715 | //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation |
---|
| 716 | //and the true pdf, extended over the interval (X(I),X(I+1)) |
---|
| 717 | G4int icase = 1; //loop code |
---|
| 718 | G4bool reLoop = false; |
---|
| 719 | err->push_back(0.); |
---|
| 720 | do |
---|
| 721 | { |
---|
| 722 | reLoop = false; |
---|
| 723 | (*err)[i] = 0.; //zero variable |
---|
| 724 | for (G4int k=0;k<nip;k++) |
---|
| 725 | { |
---|
| 726 | G4double rr = (*sumi)[k]; |
---|
| 727 | G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/ |
---|
| 728 | ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i])); |
---|
| 729 | if (k == 0 || k == nip-1) |
---|
| 730 | (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]); |
---|
| 731 | else |
---|
| 732 | (*err)[i] += std::fabs(pap-(*pdfi)[k]); |
---|
| 733 | } |
---|
| 734 | (*err)[i] *= dxi; |
---|
| 735 | |
---|
| 736 | //If err(I) is too large, the pdf is approximated by a uniform distribution |
---|
| 737 | if ((*err)[i] > 0.1*(*area)[i] && icase == 1) |
---|
| 738 | { |
---|
| 739 | (*b)[i] = 0; |
---|
| 740 | (*a)[i] = 0; |
---|
| 741 | (*c)[i] = 1.; |
---|
| 742 | icase = 2; |
---|
| 743 | reLoop = true; |
---|
| 744 | } |
---|
| 745 | }while(reLoop); |
---|
| 746 | |
---|
| 747 | if (pdfi) delete pdfi; |
---|
| 748 | if (pdfih) delete pdfih; |
---|
| 749 | if (sumi) delete sumi; |
---|
| 750 | } //end of first loop over i |
---|
| 751 | |
---|
| 752 | //Now assign last point |
---|
| 753 | (*x)[x->size()-1] = q2max; |
---|
| 754 | a->push_back(0.); |
---|
| 755 | b->push_back(0.); |
---|
| 756 | c->push_back(0.); |
---|
| 757 | err->push_back(0.); |
---|
| 758 | area->push_back(0.); |
---|
| 759 | |
---|
| 760 | if (x->size() != NUNIF || a->size() != NUNIF || |
---|
| 761 | err->size() != NUNIF || area->size() != NUNIF) |
---|
| 762 | { |
---|
| 763 | G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl; |
---|
| 764 | G4cout << "Problem in building the Table for Sampling: array dimensions do not match " << G4endl; |
---|
| 765 | G4Exception(); |
---|
| 766 | } |
---|
| 767 | |
---|
| 768 | /******************************************************************************* |
---|
| 769 | New grid points are added by halving the sub-intervals with the largest absolute error |
---|
| 770 | This is done up to np=150 points in the grid |
---|
| 771 | ********************************************************************************/ |
---|
| 772 | do |
---|
| 773 | { |
---|
| 774 | G4double maxError = 0.0; |
---|
| 775 | size_t iErrMax = 0; |
---|
| 776 | for (size_t i=0;i<err->size()-2;i++) |
---|
| 777 | { |
---|
| 778 | //maxError is the lagest of the interval errors err[i] |
---|
| 779 | if ((*err)[i] > maxError) |
---|
| 780 | { |
---|
| 781 | maxError = (*err)[i]; |
---|
| 782 | iErrMax = i; |
---|
| 783 | } |
---|
| 784 | } |
---|
| 785 | |
---|
| 786 | //OK, now I have to insert one new point in the position iErrMax |
---|
| 787 | G4double newx = 0.5*((*x)[iErrMax]+(*x)[iErrMax+1]); |
---|
| 788 | |
---|
| 789 | x->insert(x->begin()+iErrMax+1,newx); |
---|
| 790 | //Add place-holders in the other vectors |
---|
| 791 | area->insert(area->begin()+iErrMax+1,0.); |
---|
| 792 | a->insert(a->begin()+iErrMax+1,0.); |
---|
| 793 | b->insert(b->begin()+iErrMax+1,0.); |
---|
| 794 | c->insert(c->begin()+iErrMax+1,0.); |
---|
| 795 | err->insert(err->begin()+iErrMax+1,0.); |
---|
| 796 | |
---|
| 797 | //Now calculate the other parameters |
---|
| 798 | for (size_t i=iErrMax;i<=iErrMax+1;i++) |
---|
| 799 | { |
---|
| 800 | //Temporary vectors for this loop |
---|
| 801 | G4DataVector* pdfi = new G4DataVector(); |
---|
| 802 | G4DataVector* pdfih = new G4DataVector(); |
---|
| 803 | G4DataVector* sumi = new G4DataVector(); |
---|
| 804 | |
---|
| 805 | G4double dx = (*x)[i+1]-(*x)[i]; |
---|
| 806 | G4double dxi = ((*x)[i+1]-(*x)[i])/(G4double (nip-1)); |
---|
| 807 | G4double pdfmax = 0; |
---|
| 808 | for (G4int k=0;k<nip;k++) |
---|
| 809 | { |
---|
| 810 | G4double xik = (*x)[i]+k*dxi; |
---|
| 811 | G4double pdfk = std::max(GetFSquared(mat,xik),0.); |
---|
| 812 | pdfi->push_back(pdfk); |
---|
| 813 | pdfmax = std::max(pdfmax,pdfk); |
---|
| 814 | if (k < (nip-1)) |
---|
| 815 | { |
---|
| 816 | G4double xih = xik + 0.5*dxi; |
---|
| 817 | G4double pdfIK = std::max(GetFSquared(mat,xih),0.); |
---|
| 818 | pdfih->push_back(pdfIK); |
---|
| 819 | pdfmax = std::max(pdfmax,pdfIK); |
---|
| 820 | } |
---|
| 821 | } |
---|
| 822 | |
---|
| 823 | //Simpson's integration |
---|
| 824 | G4double cons = dxi*0.5*(1./3.); |
---|
| 825 | sumi->push_back(0.); |
---|
| 826 | for (G4int k=1;k<nip;k++) |
---|
| 827 | { |
---|
| 828 | G4double previous = (*sumi)[k-1]; |
---|
| 829 | G4double next = previous + cons*((*pdfi)[k-1]+4.0*(*pdfih)[k-1]+(*pdfi)[k]); |
---|
| 830 | sumi->push_back(next); |
---|
| 831 | } |
---|
| 832 | G4double lastIntegral = (*sumi)[sumi->size()-1]; |
---|
| 833 | (*area)[i] = lastIntegral; |
---|
| 834 | |
---|
| 835 | //Normalize cumulative function |
---|
| 836 | G4double factor = 1.0/lastIntegral; |
---|
| 837 | for (size_t k=0;k<sumi->size();k++) |
---|
| 838 | (*sumi)[k] *= factor; |
---|
| 839 | |
---|
| 840 | //When the PDF vanishes at one of the interval end points, its value is modified |
---|
| 841 | if ((*pdfi)[0] < 1e-35) |
---|
| 842 | (*pdfi)[0] = 1e-5*pdfmax; |
---|
| 843 | if ((*pdfi)[pdfi->size()-1] < 1e-35) |
---|
| 844 | (*pdfi)[pdfi->size()-1] = 1e-5*pdfmax; |
---|
| 845 | |
---|
| 846 | G4double pli = (*pdfi)[0]*factor; |
---|
| 847 | G4double pui = (*pdfi)[pdfi->size()-1]*factor; |
---|
| 848 | G4double B_temp = 1.0-1.0/(pli*pui*dx*dx); |
---|
| 849 | G4double A_temp = (1.0/(pli*dx))-1.0-B_temp; |
---|
| 850 | G4double C_temp = 1.0+A_temp+B_temp; |
---|
| 851 | if (C_temp < 1e-35) |
---|
| 852 | { |
---|
| 853 | (*a)[i]= 0.; |
---|
| 854 | (*b)[i] = 0.; |
---|
| 855 | (*c)[i] = 0; |
---|
| 856 | } |
---|
| 857 | else |
---|
| 858 | { |
---|
| 859 | (*a)[i]= A_temp; |
---|
| 860 | (*b)[i] = B_temp; |
---|
| 861 | (*c)[i] = C_temp; |
---|
| 862 | } |
---|
| 863 | //OK, now get ERR(I), the integral of the absolute difference between the rational interpolation |
---|
| 864 | //and the true pdf, extended over the interval (X(I),X(I+1)) |
---|
| 865 | G4int icase = 1; //loop code |
---|
| 866 | G4bool reLoop = false; |
---|
| 867 | do |
---|
| 868 | { |
---|
| 869 | reLoop = false; |
---|
| 870 | (*err)[i] = 0.; //zero variable |
---|
| 871 | for (G4int k=0;k<nip;k++) |
---|
| 872 | { |
---|
| 873 | G4double rr = (*sumi)[k]; |
---|
| 874 | G4double pap = (*area)[i]*(1.0+((*a)[i]+(*b)[i]*rr)*rr)*(1.0+((*a)[i]+(*b)[i]*rr)*rr)/ |
---|
| 875 | ((1.0-(*b)[i]*rr*rr)*(*c)[i]*((*x)[i+1]-(*x)[i])); |
---|
| 876 | if (k == 0 || k == nip-1) |
---|
| 877 | (*err)[i] += 0.5*std::fabs(pap-(*pdfi)[k]); |
---|
| 878 | else |
---|
| 879 | (*err)[i] += std::fabs(pap-(*pdfi)[k]); |
---|
| 880 | } |
---|
| 881 | (*err)[i] *= dxi; |
---|
| 882 | |
---|
| 883 | //If err(I) is too large, the pdf is approximated by a uniform distribution |
---|
| 884 | if ((*err)[i] > 0.1*(*area)[i] && icase == 1) |
---|
| 885 | { |
---|
| 886 | (*b)[i] = 0; |
---|
| 887 | (*a)[i] = 0; |
---|
| 888 | (*c)[i] = 1.; |
---|
| 889 | icase = 2; |
---|
| 890 | reLoop = true; |
---|
| 891 | } |
---|
| 892 | }while(reLoop); |
---|
| 893 | if (pdfi) delete pdfi; |
---|
| 894 | if (pdfih) delete pdfih; |
---|
| 895 | if (sumi) delete sumi; |
---|
| 896 | } |
---|
| 897 | }while(x->size() < np); |
---|
| 898 | |
---|
| 899 | if (x->size() != np || a->size() != np || |
---|
| 900 | err->size() != np || area->size() != np) |
---|
| 901 | { |
---|
| 902 | G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl; |
---|
| 903 | G4cout << "Problem in building the extended Table for Sampling: array dimensions do not match " << G4endl; |
---|
| 904 | G4Exception(); |
---|
| 905 | } |
---|
| 906 | |
---|
| 907 | /******************************************************************************* |
---|
| 908 | Renormalization |
---|
| 909 | ********************************************************************************/ |
---|
| 910 | G4double ws = 0; |
---|
| 911 | for (size_t i=0;i<np-1;i++) |
---|
| 912 | ws += (*area)[i]; |
---|
| 913 | ws = 1.0/ws; |
---|
| 914 | G4double errMax = 0; |
---|
| 915 | for (size_t i=0;i<np-1;i++) |
---|
| 916 | { |
---|
| 917 | (*area)[i] *= ws; |
---|
| 918 | (*err)[i] *= ws; |
---|
| 919 | errMax = std::max(errMax,(*err)[i]); |
---|
| 920 | } |
---|
| 921 | |
---|
| 922 | //Vector with the normalized cumulative distribution |
---|
| 923 | G4DataVector* PAC = new G4DataVector(); |
---|
| 924 | PAC->push_back(0.); |
---|
| 925 | for (size_t i=0;i<np-1;i++) |
---|
| 926 | { |
---|
| 927 | G4double previous = (*PAC)[i]; |
---|
| 928 | PAC->push_back(previous+(*area)[i]); |
---|
| 929 | } |
---|
| 930 | (*PAC)[PAC->size()-1] = 1.; |
---|
| 931 | |
---|
| 932 | /******************************************************************************* |
---|
| 933 | Pre-calculated limits for the initial binary search for subsequent sampling |
---|
| 934 | ********************************************************************************/ |
---|
| 935 | |
---|
| 936 | //G4DataVector* ITTL = new G4DataVector(); |
---|
| 937 | std::vector<size_t> *ITTL = new std::vector<size_t>; |
---|
| 938 | std::vector<size_t> *ITTU = new std::vector<size_t>; |
---|
| 939 | |
---|
| 940 | //Just create place-holders |
---|
| 941 | for (size_t i=0;i<np;i++) |
---|
| 942 | { |
---|
| 943 | ITTL->push_back(0); |
---|
| 944 | ITTU->push_back(0); |
---|
| 945 | } |
---|
| 946 | |
---|
| 947 | G4double bin = 1.0/(np-1); |
---|
| 948 | (*ITTL)[0]=0; |
---|
| 949 | for (size_t i=1;i<(np-1);i++) |
---|
| 950 | { |
---|
| 951 | G4double ptst = i*bin; |
---|
| 952 | G4bool found = false; |
---|
| 953 | for (size_t j=(*ITTL)[i-1];j<np && !found;j++) |
---|
| 954 | { |
---|
| 955 | if ((*PAC)[j] > ptst) |
---|
| 956 | { |
---|
| 957 | (*ITTL)[i] = j-1; |
---|
| 958 | (*ITTU)[i-1] = j; |
---|
| 959 | found = true; |
---|
| 960 | } |
---|
| 961 | } |
---|
| 962 | } |
---|
| 963 | (*ITTU)[ITTU->size()-2] = ITTU->size()-1; |
---|
| 964 | (*ITTU)[ITTU->size()-1] = ITTU->size()-1; |
---|
| 965 | (*ITTL)[ITTL->size()-1] = ITTU->size()-2; |
---|
| 966 | |
---|
| 967 | if (ITTU->size() != np || ITTU->size() != np) |
---|
| 968 | { |
---|
| 969 | G4cout << "G4Penelope08RayleighModel::InitializeSamplingAlgorithm" << G4endl; |
---|
| 970 | G4cout << "Problem in building the Limit Tables for Sampling: array dimensions do not match " << G4endl; |
---|
| 971 | G4Exception(); |
---|
| 972 | } |
---|
| 973 | |
---|
| 974 | |
---|
| 975 | /******************************************************************************** |
---|
| 976 | Copy tables |
---|
| 977 | ********************************************************************************/ |
---|
| 978 | G4PenelopeSamplingData* theTable = new G4PenelopeSamplingData(np); |
---|
| 979 | for (size_t i=0;i<np;i++) |
---|
| 980 | { |
---|
| 981 | theTable->AddPoint((*x)[i],(*PAC)[i],(*a)[i],(*b)[i],(*ITTL)[i],(*ITTU)[i]); |
---|
| 982 | } |
---|
| 983 | |
---|
| 984 | if (verboseLevel > 2) |
---|
| 985 | { |
---|
| 986 | G4cout << "*************************************************************************" << G4endl; |
---|
| 987 | G4cout << "Sampling table for Penelope Rayleigh scattering in " << mat->GetName() << G4endl; |
---|
| 988 | theTable->DumpTable(); |
---|
| 989 | } |
---|
| 990 | samplingTable->insert(std::make_pair(mat,theTable)); |
---|
| 991 | |
---|
| 992 | |
---|
| 993 | //Clean up temporary vectors |
---|
| 994 | if (x) delete x; |
---|
| 995 | if (a) delete a; |
---|
| 996 | if (b) delete b; |
---|
| 997 | if (c) delete c; |
---|
| 998 | if (err) delete err; |
---|
| 999 | if (area) delete area; |
---|
| 1000 | if (PAC) delete PAC; |
---|
| 1001 | if (ITTL) delete ITTL; |
---|
| 1002 | if (ITTU) delete ITTU; |
---|
| 1003 | |
---|
| 1004 | //DONE! |
---|
| 1005 | return; |
---|
| 1006 | |
---|
| 1007 | } |
---|
| 1008 | |
---|
| 1009 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 1010 | |
---|
| 1011 | void G4Penelope08RayleighModel::GetPMaxTable(const G4Material* mat) |
---|
| 1012 | { |
---|
| 1013 | if (!pMaxTable) |
---|
| 1014 | { |
---|
| 1015 | G4cout << "G4Penelope08RayleighModel::BuildPMaxTable" << G4endl; |
---|
| 1016 | G4cout << "Going to instanziate the pMaxTable !" << G4endl; |
---|
| 1017 | G4cout << "That should _not_ be here! " << G4endl; |
---|
| 1018 | pMaxTable = new std::map<const G4Material*,G4PhysicsFreeVector*>; |
---|
| 1019 | } |
---|
| 1020 | //check if the table is already there |
---|
| 1021 | if (pMaxTable->count(mat)) |
---|
| 1022 | return; |
---|
| 1023 | |
---|
| 1024 | //otherwise build it |
---|
| 1025 | if (!samplingTable) |
---|
| 1026 | { |
---|
| 1027 | G4cout << "G4Penelope08RayleighModel::BuildPMaxTable" << G4endl; |
---|
| 1028 | G4cout << "WARNING: samplingTable is not properly instantiated" << G4endl; |
---|
| 1029 | G4Exception(); |
---|
| 1030 | } |
---|
| 1031 | |
---|
| 1032 | if (!samplingTable->count(mat)) |
---|
| 1033 | InitializeSamplingAlgorithm(mat); |
---|
| 1034 | |
---|
| 1035 | G4PenelopeSamplingData *theTable = samplingTable->find(mat)->second; |
---|
| 1036 | size_t tablePoints = theTable->GetNumberOfStoredPoints(); |
---|
| 1037 | |
---|
| 1038 | size_t nOfEnergyPoints = logEnergyGridPMax.size(); |
---|
| 1039 | G4PhysicsFreeVector* theVec = new G4PhysicsFreeVector(nOfEnergyPoints); |
---|
| 1040 | |
---|
| 1041 | const size_t nip = 51; //hard-coded in Penelope |
---|
| 1042 | |
---|
| 1043 | for (size_t ie=0;ie<logEnergyGridPMax.size();ie++) |
---|
| 1044 | { |
---|
| 1045 | G4double energy = exp(logEnergyGridPMax[ie]); |
---|
| 1046 | G4double Qm = 2.0*energy/electron_mass_c2; //this is non-dimensional now |
---|
| 1047 | G4double Qm2 = Qm*Qm; |
---|
| 1048 | G4double firstQ2 = theTable->GetX(0); |
---|
| 1049 | G4double lastQ2 = theTable->GetX(tablePoints-1); |
---|
| 1050 | G4double thePMax = 0; |
---|
| 1051 | |
---|
| 1052 | if (Qm2 > firstQ2) |
---|
| 1053 | { |
---|
| 1054 | if (Qm2 < lastQ2) |
---|
| 1055 | { |
---|
| 1056 | //bisection to look for the index of Qm |
---|
| 1057 | size_t lowerBound = 0; |
---|
| 1058 | size_t upperBound = tablePoints-1; |
---|
| 1059 | while (lowerBound <= upperBound) |
---|
| 1060 | { |
---|
| 1061 | size_t midBin = (lowerBound + upperBound)/2; |
---|
| 1062 | if( Qm2 < theTable->GetX(midBin)) |
---|
| 1063 | { upperBound = midBin-1; } |
---|
| 1064 | else |
---|
| 1065 | { lowerBound = midBin+1; } |
---|
| 1066 | } |
---|
| 1067 | //upperBound is the output (but also lowerBounf --> should be the same!) |
---|
| 1068 | G4double Q1 = theTable->GetX(upperBound); |
---|
| 1069 | G4double Q2 = Qm2; |
---|
| 1070 | G4double DQ = (Q2-Q1)/((G4double)(nip-1)); |
---|
| 1071 | G4double theA = theTable->GetA(upperBound); |
---|
| 1072 | G4double theB = theTable->GetB(upperBound); |
---|
| 1073 | G4double thePAC = theTable->GetPAC(upperBound); |
---|
| 1074 | G4DataVector* fun = new G4DataVector(); |
---|
| 1075 | for (size_t k=0;k<nip;k++) |
---|
| 1076 | { |
---|
| 1077 | G4double qi = Q1 + k*DQ; |
---|
| 1078 | G4double tau = (qi-Q1)/ |
---|
| 1079 | (theTable->GetX(upperBound+1)-Q1); |
---|
| 1080 | G4double con1 = 2.0*theB*tau; |
---|
| 1081 | G4double ci = 1.0+theA+theB; |
---|
| 1082 | G4double con2 = ci-theA*tau; |
---|
| 1083 | G4double etap = 0; |
---|
| 1084 | if (std::fabs(con1) > 1.0e-16*std::fabs(con2)) |
---|
| 1085 | etap = con2*(1.0-std::sqrt(1.0-2.0*tau*con1/(con2*con2)))/con1; |
---|
| 1086 | else |
---|
| 1087 | etap = tau/con2; |
---|
| 1088 | G4double theFun = (theTable->GetPAC(upperBound+1)-thePAC)* |
---|
| 1089 | (1.0+(theA+theB*etap)*etap)*(1.0+(theA+theB*etap)*etap)/ |
---|
| 1090 | ((1.0-theB*etap*etap)*ci*(theTable->GetX(upperBound+1)-Q1)); |
---|
| 1091 | fun->push_back(theFun); |
---|
| 1092 | } |
---|
| 1093 | //Now intergrate numerically the fun Cavalieri-Simpson's method |
---|
| 1094 | G4DataVector* sum = new G4DataVector; |
---|
| 1095 | G4double CONS = DQ*(1./12.); |
---|
| 1096 | G4double HCONS = 0.5*CONS; |
---|
| 1097 | sum->push_back(0.); |
---|
| 1098 | G4double secondPoint = (*sum)[0] + |
---|
| 1099 | (5.0*(*fun)[0]+8.0*(*fun)[1]-(*fun)[2])*CONS; |
---|
| 1100 | sum->push_back(secondPoint); |
---|
| 1101 | for (size_t hh=2;hh<nip-1;hh++) |
---|
| 1102 | { |
---|
| 1103 | G4double previous = (*sum)[hh-1]; |
---|
| 1104 | G4double next = previous+(13.0*((*fun)[hh-1]+(*fun)[hh])- |
---|
| 1105 | (*fun)[hh+1]-(*fun)[hh-2])*HCONS; |
---|
| 1106 | sum->push_back(next); |
---|
| 1107 | } |
---|
| 1108 | G4double last = (*sum)[nip-2]+(5.0*(*fun)[nip-1]+8.0*(*fun)[nip-2]- |
---|
| 1109 | (*fun)[nip-3])*CONS; |
---|
| 1110 | sum->push_back(last); |
---|
| 1111 | thePMax = thePAC + (*sum)[sum->size()-1]; //last point |
---|
| 1112 | if (fun) delete fun; |
---|
| 1113 | if (sum) delete sum; |
---|
| 1114 | } |
---|
| 1115 | else |
---|
| 1116 | { |
---|
| 1117 | thePMax = 1.0; |
---|
| 1118 | } |
---|
| 1119 | } |
---|
| 1120 | else |
---|
| 1121 | { |
---|
| 1122 | thePMax = theTable->GetPAC(0); |
---|
| 1123 | } |
---|
| 1124 | |
---|
| 1125 | //Write number in the table |
---|
| 1126 | theVec->PutValue(ie,energy,thePMax); |
---|
| 1127 | } |
---|
| 1128 | |
---|
| 1129 | pMaxTable->insert(std::make_pair(mat,theVec)); |
---|
| 1130 | return; |
---|
| 1131 | |
---|
| 1132 | } |
---|
| 1133 | |
---|
| 1134 | //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... |
---|
| 1135 | |
---|
| 1136 | void G4Penelope08RayleighModel::DumpFormFactorTable(const G4Material* mat) |
---|
| 1137 | { |
---|
| 1138 | G4cout << "*****************************************************************" << G4endl; |
---|
| 1139 | G4cout << "G4Penelope08RayleighModel: Form Factor Table for " << mat->GetName() << G4endl; |
---|
| 1140 | //try to use the same format as Penelope-Fortran, namely Q (/m_e*c) and F |
---|
| 1141 | G4cout << "Q/(m_e*c) F(Q) " << G4endl; |
---|
| 1142 | G4cout << "*****************************************************************" << G4endl; |
---|
| 1143 | if (!logFormFactorTable->count(mat)) |
---|
| 1144 | BuildFormFactorTable(mat); |
---|
| 1145 | |
---|
| 1146 | G4PhysicsFreeVector* theVec = logFormFactorTable->find(mat)->second; |
---|
| 1147 | for (size_t i=0;i<theVec->GetVectorLength();i++) |
---|
| 1148 | { |
---|
| 1149 | G4double logQ2 = theVec->GetLowEdgeEnergy(i); |
---|
| 1150 | G4double Q = exp(0.5*logQ2); |
---|
| 1151 | G4double logF2 = (*theVec)[i]; |
---|
| 1152 | G4double F = exp(0.5*logF2); |
---|
| 1153 | G4cout << Q << " " << F << G4endl; |
---|
| 1154 | } |
---|
| 1155 | //DONE |
---|
| 1156 | return; |
---|
| 1157 | } |
---|