00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #include "G4LivermoreRayleighModel.hh"
00032 #include "G4SystemOfUnits.hh"
00033 #include "G4RayleighAngularGenerator.hh"
00034
00035 using namespace std;
00036
00037
00038
00039
00040 G4LivermoreRayleighModel::G4LivermoreRayleighModel()
00041 :G4VEmModel("LivermoreRayleigh"),isInitialised(false),maxZ(100)
00042 {
00043 fParticleChange = 0;
00044 lowEnergyLimit = 10 * eV;
00045
00046 dataCS.resize(maxZ+1,0);
00047
00048 SetAngularDistribution(new G4RayleighAngularGenerator());
00049
00050 verboseLevel= 0;
00051
00052
00053
00054
00055
00056 if(verboseLevel > 0)
00057 {
00058 G4cout << "G4LivermoreRayleighModel is constructed " << G4endl;
00059 }
00060
00061 }
00062
00063
00064
00065 G4LivermoreRayleighModel::~G4LivermoreRayleighModel()
00066 {
00067 for(G4int i=0; i<=maxZ; ++i) { delete dataCS[i]; }
00068 }
00069
00070
00071
00072 void G4LivermoreRayleighModel::Initialise(const G4ParticleDefinition* particle,
00073 const G4DataVector& cuts)
00074 {
00075 if (verboseLevel > 1)
00076 {
00077 G4cout << "Calling Initialise() of G4LivermoreRayleighModel." << G4endl
00078 << "Energy range: "
00079 << LowEnergyLimit() / eV << " eV - "
00080 << HighEnergyLimit() / GeV << " GeV"
00081 << G4endl;
00082 }
00083
00084
00085 InitialiseElementSelectors(particle, cuts);
00086
00087
00088
00089 char* path = getenv("G4LEDATA");
00090 G4ProductionCutsTable* theCoupleTable =
00091 G4ProductionCutsTable::GetProductionCutsTable();
00092 G4int numOfCouples = theCoupleTable->GetTableSize();
00093
00094 for(G4int i=0; i<numOfCouples; ++i)
00095 {
00096 const G4MaterialCutsCouple* couple =
00097 theCoupleTable->GetMaterialCutsCouple(i);
00098 const G4Material* material = couple->GetMaterial();
00099 const G4ElementVector* theElementVector = material->GetElementVector();
00100 G4int nelm = material->GetNumberOfElements();
00101
00102 for (G4int j=0; j<nelm; ++j)
00103 {
00104 G4int Z = G4lrint((*theElementVector)[j]->GetZ());
00105 if(Z < 1) { Z = 1; }
00106 else if(Z > maxZ) { Z = maxZ; }
00107 if( (!dataCS[Z]) ) { ReadData(Z, path); }
00108 }
00109 }
00110
00111
00112
00113
00114 if(isInitialised) { return; }
00115 fParticleChange = GetParticleChangeForGamma();
00116 isInitialised = true;
00117
00118 }
00119
00120
00121
00122 void G4LivermoreRayleighModel::ReadData(size_t Z, const char* path)
00123 {
00124 if (verboseLevel > 1)
00125 {
00126 G4cout << "Calling ReadData() of G4LivermoreRayleighModel"
00127 << G4endl;
00128 }
00129
00130 if(dataCS[Z]) { return; }
00131
00132 const char* datadir = path;
00133
00134 if(!datadir)
00135 {
00136 datadir = getenv("G4LEDATA");
00137 if(!datadir)
00138 {
00139 G4Exception("G4LivermoreRayleighModelModel::ReadData()","em0006",
00140 FatalException,
00141 "Environment variable G4LEDATA not defined");
00142 return;
00143 }
00144 }
00145
00146
00147
00148 dataCS[Z] = new G4LPhysicsFreeVector();
00149
00150
00151
00152
00153 std::ostringstream ostCS;
00154 ostCS << datadir << "/livermore/rayl/re-cs-" << Z <<".dat";
00155 std::ifstream finCS(ostCS.str().c_str());
00156
00157 if( !finCS .is_open() )
00158 {
00159 G4ExceptionDescription ed;
00160 ed << "G4LivermoreRayleighModel data file <" << ostCS.str().c_str()
00161 << "> is not opened!" << G4endl;
00162 G4Exception("G4LivermoreRayleighModel::ReadData()","em0003",FatalException,
00163 ed,"G4LEDATA version should be G4EMLOW6.27 or later.");
00164 return;
00165 }
00166 else
00167 {
00168 if(verboseLevel > 3) {
00169 G4cout << "File " << ostCS.str()
00170 << " is opened by G4LivermoreRayleighModel" << G4endl;
00171 }
00172 dataCS[Z]->Retrieve(finCS, true);
00173 }
00174 }
00175
00176
00177
00178 G4double G4LivermoreRayleighModel::ComputeCrossSectionPerAtom(
00179 const G4ParticleDefinition*,
00180 G4double GammaEnergy,
00181 G4double Z, G4double,
00182 G4double, G4double)
00183 {
00184 if (verboseLevel > 1)
00185 {
00186 G4cout << "Calling ComputeCrossSectionPerAtom() of G4LivermoreRayleighModel"
00187 << G4endl;
00188 }
00189
00190 if(GammaEnergy < lowEnergyLimit) { return 0.0; }
00191
00192 G4double xs = 0.0;
00193
00194 G4int intZ = G4lrint(Z);
00195
00196 if(intZ < 1 || intZ > maxZ) { return xs; }
00197
00198 G4LPhysicsFreeVector* pv = dataCS[intZ];
00199
00200
00201 if(!pv)
00202 {
00203 char* path = getenv("G4LEDATA");
00204 ReadData(intZ, path);
00205 pv = dataCS[intZ];
00206 if(!pv) { return xs; }
00207 }
00208
00209 G4int n = pv->GetVectorLength() - 1;
00210 G4double e = GammaEnergy/MeV;
00211 if(e >= pv->Energy(n)) {
00212 xs = (*pv)[n]/(e*e);
00213 } else if(e >= pv->Energy(0)) {
00214 xs = pv->Value(e)/(e*e);
00215 }
00216
00217 if(verboseLevel > 0)
00218 {
00219 G4cout << "****** DEBUG: tcs value for Z=" << Z << " at energy (MeV)=" << e << G4endl;
00220 G4cout << " cs (Geant4 internal unit)=" << xs << G4endl;
00221 G4cout << " -> first E*E*cs value in CS data file (iu) =" << (*pv)[0] << G4endl;
00222 G4cout << " -> last E*E*cs value in CS data file (iu) =" << (*pv)[n] << G4endl;
00223 G4cout << "*********************************************************" << G4endl;
00224 }
00225 return xs;
00226 }
00227
00228
00229
00230 void G4LivermoreRayleighModel::SampleSecondaries(
00231 std::vector<G4DynamicParticle*>*,
00232 const G4MaterialCutsCouple* couple,
00233 const G4DynamicParticle* aDynamicGamma,
00234 G4double, G4double)
00235 {
00236 if (verboseLevel > 1) {
00237 G4cout << "Calling SampleSecondaries() of G4LivermoreRayleighModel"
00238 << G4endl;
00239 }
00240 G4double photonEnergy0 = aDynamicGamma->GetKineticEnergy();
00241
00242
00243 if (photonEnergy0 <= lowEnergyLimit)
00244 {
00245 fParticleChange->ProposeTrackStatus(fStopAndKill);
00246 fParticleChange->SetProposedKineticEnergy(0.);
00247 fParticleChange->ProposeLocalEnergyDeposit(photonEnergy0);
00248 return ;
00249 }
00250
00251
00252 const G4ParticleDefinition* particle = aDynamicGamma->GetDefinition();
00253 const G4Element* elm = SelectRandomAtom(couple,particle,photonEnergy0);
00254 G4int Z = G4lrint(elm->GetZ());
00255
00256
00257
00258 G4ThreeVector photonDirection =
00259 GetAngularDistribution()->SampleDirection(aDynamicGamma,
00260 photonEnergy0, Z, couple->GetMaterial());
00261 fParticleChange->ProposeMomentumDirection(photonDirection);
00262 }
00263
00264