// // ******************************************************************** // * License and Disclaimer * // * * // * The Geant4 software is copyright of the Copyright Holders of * // * the Geant4 Collaboration. It is provided under the terms and * // * conditions of the Geant4 Software License, included in the file * // * LICENSE and available at http://cern.ch/geant4/license . These * // * include a list of copyright holders. * // * * // * Neither the authors of this software system, nor their employing * // * institutes,nor the agencies providing financial support for this * // * work make any representation or warranty, express or implied, * // * regarding this software system or assume any liability for its * // * use. Please see the license in the file LICENSE and URL above * // * for the full disclaimer and the limitation of liability. * // * * // * This code implementation is the result of the scientific and * // * technical work of the GEANT4 collaboration. * // * By using, copying, modifying or distributing the software (or * // * any work based on the software) you agree to acknowledge its * // * use in resulting scientific publications, and indicate your * // * acceptance of all terms of the Geant4 Software license. * // ******************************************************************** // // // The lust update: M.V. Kossov, CERN/ITEP(Moscow) 17-June-02 // GEANT4 tag $Name: geant4-09-03 $ // // // G4 Physics class: G4PhotoNuclearCrossSection for gamma+A cross sections // Created: M.V. Kossov, CERN/ITEP(Moscow), 10-OCT-01 // The last update: M.V. Kossov, CERN/ITEP (Moscow) 17-May-02 // //#define debug //#define pdebug //#define debug3 //#define debugn //#define debugs #include "G4PhotoNuclearCrossSection.hh" // Initialization of the statics G4int G4PhotoNuclearCrossSection::lastN=0; // The last N of calculated nucleus G4int G4PhotoNuclearCrossSection::lastZ=0; // The last Z of calculated nucleus G4double G4PhotoNuclearCrossSection::lastSig=0.; // Last value of the Cross Section G4double* G4PhotoNuclearCrossSection::lastGDR=0; // Pointer to the last array of GDR cross sections G4double* G4PhotoNuclearCrossSection::lastHEN=0; // Pointer to the last array of HEn cross sections G4double G4PhotoNuclearCrossSection::lastE=0.; // Last used in the cross section Energy G4double G4PhotoNuclearCrossSection::lastTH=0.; // Last value of the Energy Threshold (A-dependent) G4double G4PhotoNuclearCrossSection::lastSP=0.; // Last value of the ShadowingPomeron (A-dependent) std::vector G4PhotoNuclearCrossSection::GDR; // Vector of pointers to the GDRPhotonuclearCrossSection std::vector G4PhotoNuclearCrossSection::HEN; // Vector of pointers to the HighEnPhotonuclearCrossSect G4PhotoNuclearCrossSection::G4PhotoNuclearCrossSection() { } G4PhotoNuclearCrossSection::~G4PhotoNuclearCrossSection() { std::vector::iterator pos; for(pos=GDR.begin(); posGetNumberOfIsotopes(); G4double cross_section = 0; if (nIso) { G4double psig; G4IsotopeVector* isoVector = anEle->GetIsotopeVector(); G4double* abundVector = anEle->GetRelativeAbundanceVector(); G4double ZZ; G4double AA; for (G4int i = 0; i < nIso; i++) { ZZ = G4double( (*isoVector)[i]->GetZ() ); AA = G4double( (*isoVector)[i]->GetN() ); psig = GetIsoZACrossSection(aPart, ZZ, AA, temperature); cross_section += psig*abundVector[i]; } } else { cross_section = GetIsoZACrossSection(aPart, anEle->GetZ(), anEle->GetN(), temperature); } return cross_section; } G4double G4PhotoNuclearCrossSection::GetIsoZACrossSection(const G4DynamicParticle* aPart, G4double ZZ, G4double AA, G4double /*temperature*/) { static const G4double THmin=2.; // minimum Energy Threshold static const G4double dE=1.; // step for the GDR table static const G4int nL=105; // A#of GDResonance points in E // (each MeV from 2 to 106) static const G4double Emin=THmin+(nL-1)*dE; // minE for the HighE part static const G4double Emax=50000.; // maxE for the HighE part static const G4int nH=224; // A#of HResonance points in lnE static const G4double milE=std::log(Emin); // Low logarithm energy for // the HighE part static const G4double malE=std::log(Emax); // High logarithm energy // (each 2.75 percent) static const G4double dlE=(malE-milE)/(nH-1); // Step in logarithm energy // in the HighE part // //static const G4double shd=1.075-.0023*std::log(2.); // HE PomShadowing(D) static const G4double shd=1.0734; // HE PomShadowing(D) static const G4double shc=0.072; // HE Shadowing constant static const G4double poc=0.0375; // HE Pomeron coefficient static const G4double pos=16.5; // HE Pomeron shift static const G4double reg=.11; // HE Reggeon slope //static const G4double shp=1.075; // HE PomShadowing(P) // Associative memory for acceleration static std::vector colN; // N of calculated nuclei static std::vector colZ; // Z of calculated nuclei static std::vector spA; // shadowing coefficients (A-dependent) static std::vector eTH; // energy threshold (A-dependent) // const G4double Energy = aPart->GetKineticEnergy()/MeV; const G4int targetAtomicNumber = static_cast(AA+.499); //@@ Nat mixture (?!) const G4int targZ = static_cast(ZZ); const G4int targN = targetAtomicNumber-targZ; #ifdef debug G4cout << "G4PhotoNuclearCrossSection::GetCS:N=" << targN << ",Z=" << targZ << ",E=" << Energy << G4endl; #endif if (EnergyGetDefinition()->GetPDGEncoding() == 22 ) { G4double A=targN+targZ; if(targN!=lastN || targZ!=lastZ) // Otherwise the set of parameters is ready { lastN = targN; // The last N of calculated nucleus lastZ = targZ; // The last Z of calculated nucleus G4int n=colN.size(); // Size of Associated Memory G4bool in=false; // The nucleus is in the AssocMem DB if(n) { for(G4int i=0; i=Xh) return YN[N-1];//-| // G4double Xp=0.; // | // G4int j=0; // | // while (X>Xj && j=nLA) k=nLA-1; // Extrapolation from the last bin (U) G4int k1=k-1; G4double xi=LA[k1]; G4double b=(a-xi)/(LA[k]-xi); for(G4int m=0; m1.5) { G4double yi=SL[k1][m]; y[m]=yi+(SL[k][m]-yi)*b; #ifdef debugs if(y[m]<0.) G4cout << "G4PhotNucCS::GetF:y=" << y[m] << ",k=" << k << ",yi=" << yi << ",ya=" << SL[k][m] << ",b=" << b << ",xi=" << xi << ",xa=" << LA[k] << ",a=" << a << G4endl; #endif } else y[m]=0.; } r=1; } if(!h) // High Energy part is not filled { G4int k=0; for(k=1; k=nHA) k=nHA-1; // Extrapolation from the last bin (Pu) G4int k1=k-1; G4double xi=HA[k1]; G4double b=(a-xi)/(HA[k]-xi); for(G4int m=0; m