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
00032
00033
00034
00035
00036
00037
00038
00039
00040
00041
00042
00043
00044
00045
00046
00047
00048
00049
00050
00051
00052
00053 #include "G4ContinuumGammaTransition.hh"
00054 #include "G4VLevelDensityParameter.hh"
00055 #include "G4ConstantLevelDensityParameter.hh"
00056 #include "G4RandGeneralTmp.hh"
00057 #include "G4PhysicalConstants.hh"
00058 #include "G4SystemOfUnits.hh"
00059 #include "Randomize.hh"
00060 #include "G4Pow.hh"
00061
00062
00063
00064
00065
00066 G4ContinuumGammaTransition::G4ContinuumGammaTransition(
00067 const G4NuclearLevelManager* levelManager,
00068 G4int Z, G4int A,
00069 G4double excitation,
00070 G4int verbose):
00071 _nucleusA(A), _nucleusZ(Z), _excitation(excitation), _levelManager(levelManager)
00072 {
00073 G4double eTolerance = 0.;
00074 G4int lastButOne = _levelManager->NumberOfLevels() - 2;
00075 if (lastButOne >= 0)
00076 {
00077 eTolerance = (_levelManager->MaxLevelEnergy() -
00078 _levelManager->GetLevel(lastButOne)->Energy());
00079 if (eTolerance < 0.) eTolerance = 0.;
00080 }
00081
00082
00083 _verbose = verbose;
00084 _eGamma = 0.;
00085 _gammaCreationTime = 0.;
00086
00087 _maxLevelE = _levelManager->MaxLevelEnergy() + eTolerance;
00088 _minLevelE = _levelManager->MinLevelEnergy();
00089
00090
00091 _eMin = 0.001 * MeV;
00092
00093 G4double energyGDR = (40.3 / G4Pow::GetInstance()->powZ(_nucleusA,0.2) ) * MeV;
00094
00095 G4double widthGDR = 0.30 * energyGDR;
00096
00097 G4double factor = 5;
00098 _eMax = energyGDR + factor * widthGDR;
00099 if (_eMax > excitation) _eMax = _excitation;
00100
00101 }
00102
00103
00104
00105
00106
00107 G4ContinuumGammaTransition::~G4ContinuumGammaTransition()
00108 {}
00109
00110 void G4ContinuumGammaTransition::SelectGamma()
00111 {
00112
00113 _eGamma = 0.;
00114
00115 G4int nBins = 200;
00116 G4double sampleArray[200];
00117 G4int i;
00118 for (i=0; i<nBins; i++)
00119 {
00120 G4double e = _eMin + ( (_eMax - _eMin) / nBins) * i;
00121 sampleArray[i] = E1Pdf(e);
00122
00123 if(_verbose > 10)
00124 G4cout << "*---* G4ContinuumTransition: e = " << e
00125 << " pdf = " << sampleArray[i] << G4endl;
00126 }
00127 G4RandGeneralTmp randGeneral(sampleArray, nBins);
00128 G4double random = randGeneral.shoot();
00129
00130 _eGamma = _eMin + (_eMax - _eMin) * random;
00131
00132 G4double finalExcitation = _excitation - _eGamma;
00133
00134 if(_verbose > 10) {
00135 G4cout << "*---*---* G4ContinuumTransition: eGamma = " << _eGamma
00136 << " finalExcitation = " << finalExcitation
00137 << " random = " << random << G4endl;
00138 }
00139
00140 if(finalExcitation < _minLevelE/2.)
00141 {
00142 _eGamma = _excitation;
00143 finalExcitation = 0.;
00144 }
00145
00146 if (finalExcitation < _maxLevelE && finalExcitation > 0.)
00147 {
00148 G4double levelE = _levelManager->NearestLevel(finalExcitation)->Energy();
00149 G4double diff = finalExcitation - levelE;
00150 _eGamma = _eGamma + diff;
00151 }
00152
00153 _gammaCreationTime = GammaTime();
00154
00155 if(_verbose > 10) {
00156 G4cout << "*---*---* G4ContinuumTransition: _gammaCreationTime = "
00157 << _gammaCreationTime/second << G4endl;
00158 }
00159 return;
00160 }
00161
00162 G4double G4ContinuumGammaTransition::GetGammaEnergy()
00163 {
00164 return _eGamma;
00165 }
00166
00167 G4double G4ContinuumGammaTransition::GetGammaCreationTime()
00168 {
00169 return _gammaCreationTime;
00170 }
00171
00172
00173 void G4ContinuumGammaTransition::SetEnergyFrom(G4double energy)
00174 {
00175 if (energy > 0.) _excitation = energy;
00176 }
00177
00178
00179 G4double G4ContinuumGammaTransition::E1Pdf(G4double e)
00180 {
00181 G4double theProb = 0.0;
00182 G4double U = std::max(0.0, _excitation - e);
00183
00184 if(e < 0.0 || _excitation < 0.0) { return theProb; }
00185
00186 G4ConstantLevelDensityParameter ldPar;
00187 G4double aLevelDensityParam =
00188 ldPar.LevelDensityParameter(_nucleusA,_nucleusZ,_excitation);
00189
00190
00191
00192 G4double coeff = std::exp(2.0*(std::sqrt(aLevelDensityParam*U)
00193 - std::sqrt(aLevelDensityParam*_excitation)));
00194
00195
00196
00197
00198
00199
00200
00201
00202
00203
00204
00205 G4double sigma0 = 2.5 * _nucleusA;
00206
00207 G4double Egdp = (40.3 /G4Pow::GetInstance()->powZ(_nucleusA,0.2) )*MeV;
00208 G4double GammaR = 0.30 * Egdp;
00209
00210 const G4double normC = 1.0 / (pi * hbarc)*(pi * hbarc);
00211
00212 G4double numerator = sigma0 * e*e * GammaR*GammaR;
00213 G4double denominator = (e*e - Egdp*Egdp)* (e*e - Egdp*Egdp) + GammaR*GammaR*e*e;
00214
00215
00216 G4double sigmaAbs = numerator/denominator ;
00217
00218 if(_verbose > 20) {
00219 G4cout << ".. " << Egdp << " .. " << GammaR
00220 << " .. " << normC << " .. " << sigmaAbs
00221 << " .. " << e*e << " .. " << coeff
00222 << G4endl;
00223 }
00224
00225
00226 theProb = sigmaAbs * e*e * coeff;
00227
00228 return theProb;
00229 }
00230
00231
00232 G4double G4ContinuumGammaTransition::GammaTime()
00233 {
00234
00235 G4double GammaR = 0.30 * (40.3 /G4Pow::GetInstance()->powZ(_nucleusA,0.2) )*MeV;
00236 G4double tau = hbar_Planck/GammaR;
00237 G4double creationTime = -tau*std::log(G4UniformRand());
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250
00251
00252
00253
00254
00255 return creationTime;
00256 }
00257
00258
00259
00260
00261
00262
00263
00264
00265
00266
00267
00268