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 #include "G4VEmModel.hh"
00053 #include "G4LossTableManager.hh"
00054 #include "G4ProductionCutsTable.hh"
00055 #include "G4ParticleChangeForLoss.hh"
00056 #include "G4ParticleChangeForGamma.hh"
00057
00058
00059
00060
00061 G4VEmModel::G4VEmModel(const G4String& nam):
00062 flucModel(0),anglModel(0), name(nam), lowLimit(0.1*CLHEP::keV),
00063 highLimit(100.0*CLHEP::TeV),eMinActive(0.0),eMaxActive(DBL_MAX),
00064 polarAngleLimit(CLHEP::pi),secondaryThreshold(DBL_MAX),
00065 theLPMflag(false),flagDeexcitation(false),flagForceBuildTable(false),
00066 pParticleChange(0),xSectionTable(0),theDensityFactor(0),theDensityIdx(0),
00067 fCurrentCouple(0),fCurrentElement(0),
00068 nsec(5)
00069 {
00070 xsec.resize(nsec);
00071 nSelectors = 0;
00072 G4LossTableManager::Instance()->Register(this);
00073 }
00074
00075
00076
00077 G4VEmModel::~G4VEmModel()
00078 {
00079 G4LossTableManager::Instance()->DeRegister(this);
00080 G4int n = elmSelectors.size();
00081 if(n > 0) {
00082 for(G4int i=0; i<n; ++i) {
00083 delete elmSelectors[i];
00084 }
00085 }
00086 delete anglModel;
00087 if(xSectionTable) {
00088 xSectionTable->clearAndDestroy();
00089 delete xSectionTable;
00090 }
00091 }
00092
00093
00094
00095 G4ParticleChangeForLoss* G4VEmModel::GetParticleChangeForLoss()
00096 {
00097 G4ParticleChangeForLoss* p = 0;
00098 if (pParticleChange) {
00099 p = static_cast<G4ParticleChangeForLoss*>(pParticleChange);
00100 } else {
00101 p = new G4ParticleChangeForLoss();
00102 pParticleChange = p;
00103 }
00104 return p;
00105 }
00106
00107
00108
00109 G4ParticleChangeForGamma* G4VEmModel::GetParticleChangeForGamma()
00110 {
00111 G4ParticleChangeForGamma* p = 0;
00112 if (pParticleChange) {
00113 p = static_cast<G4ParticleChangeForGamma*>(pParticleChange);
00114 } else {
00115 p = new G4ParticleChangeForGamma();
00116 pParticleChange = p;
00117 }
00118 return p;
00119 }
00120
00121
00122
00123 void G4VEmModel::InitialiseElementSelectors(const G4ParticleDefinition* p,
00124 const G4DataVector& cuts)
00125 {
00126
00127 G4LossTableManager* man = G4LossTableManager::Instance();
00128
00129
00130
00131
00132
00133 G4bool spline = false;
00134
00135
00136
00137
00138
00139
00140 G4int nbins = G4int(man->GetNumberOfBinsPerDecade()
00141 * std::log10(highLimit/lowLimit) / 6.0);
00142 if(nbins < 5) { nbins = 5; }
00143
00144 G4ProductionCutsTable* theCoupleTable=
00145 G4ProductionCutsTable::GetProductionCutsTable();
00146 G4int numOfCouples = theCoupleTable->GetTableSize();
00147
00148
00149 if(numOfCouples > nSelectors) {
00150 elmSelectors.resize(numOfCouples,0);
00151 nSelectors = numOfCouples;
00152 }
00153
00154
00155 for(G4int i=0; i<numOfCouples; ++i) {
00156 fCurrentCouple = theCoupleTable->GetMaterialCutsCouple(i);
00157 const G4Material* material = fCurrentCouple->GetMaterial();
00158 G4int idx = fCurrentCouple->GetIndex();
00159
00160
00161 G4bool create = true;
00162 if(elmSelectors[i]) {
00163 if(material == elmSelectors[i]->GetMaterial()) { create = false; }
00164 else { delete elmSelectors[i]; }
00165 }
00166 if(create) {
00167 elmSelectors[i] = new G4EmElementSelector(this,material,nbins,
00168 lowLimit,highLimit,spline);
00169 }
00170 elmSelectors[i]->Initialise(p, cuts[idx]);
00171
00172 }
00173 }
00174
00175
00176
00177 G4double G4VEmModel::ComputeDEDXPerVolume(const G4Material*,
00178 const G4ParticleDefinition*,
00179 G4double,G4double)
00180 {
00181 return 0.0;
00182 }
00183
00184
00185
00186 G4double G4VEmModel::CrossSectionPerVolume(const G4Material* material,
00187 const G4ParticleDefinition* p,
00188 G4double ekin,
00189 G4double emin,
00190 G4double emax)
00191 {
00192 SetupForMaterial(p, material, ekin);
00193 G4double cross = 0.0;
00194 const G4ElementVector* theElementVector = material->GetElementVector();
00195 const G4double* theAtomNumDensityVector = material->GetVecNbOfAtomsPerVolume();
00196 G4int nelm = material->GetNumberOfElements();
00197 if(nelm > nsec) {
00198 xsec.resize(nelm);
00199 nsec = nelm;
00200 }
00201 for (G4int i=0; i<nelm; i++) {
00202 cross += theAtomNumDensityVector[i]*
00203 ComputeCrossSectionPerAtom(p,(*theElementVector)[i],ekin,emin,emax);
00204 xsec[i] = cross;
00205 }
00206 return cross;
00207 }
00208
00209
00210
00211 void G4VEmModel::StartTracking(G4Track*)
00212 {}
00213
00214
00215
00216 const G4Element* G4VEmModel::SelectRandomAtom(const G4Material* material,
00217 const G4ParticleDefinition* pd,
00218 G4double kinEnergy,
00219 G4double tcut,
00220 G4double tmax)
00221 {
00222 const G4ElementVector* theElementVector = material->GetElementVector();
00223 G4int n = material->GetNumberOfElements() - 1;
00224 fCurrentElement = (*theElementVector)[n];
00225 if (n > 0) {
00226 G4double x = G4UniformRand()*
00227 G4VEmModel::CrossSectionPerVolume(material,pd,kinEnergy,tcut,tmax);
00228 for(G4int i=0; i<n; ++i) {
00229 if (x <= xsec[i]) {
00230 fCurrentElement = (*theElementVector)[i];
00231 break;
00232 }
00233 }
00234 }
00235 return fCurrentElement;
00236 }
00237
00238
00239
00240 G4double G4VEmModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition*,
00241 G4double, G4double, G4double,
00242 G4double, G4double)
00243 {
00244 return 0.0;
00245 }
00246
00247
00248
00249 void G4VEmModel::DefineForRegion(const G4Region*)
00250 {}
00251
00252
00253
00254 G4double G4VEmModel::ChargeSquareRatio(const G4Track& track)
00255 {
00256 return GetChargeSquareRatio(track.GetParticleDefinition(),
00257 track.GetMaterial(), track.GetKineticEnergy());
00258 }
00259
00260
00261
00262 G4double G4VEmModel::GetChargeSquareRatio(const G4ParticleDefinition* p,
00263 const G4Material*, G4double)
00264 {
00265 G4double q = p->GetPDGCharge()/CLHEP::eplus;
00266 return q*q;
00267 }
00268
00269
00270
00271 G4double G4VEmModel::GetParticleCharge(const G4ParticleDefinition* p,
00272 const G4Material*, G4double)
00273 {
00274 return p->GetPDGCharge();
00275 }
00276
00277
00278
00279 void G4VEmModel::CorrectionsAlongStep(const G4MaterialCutsCouple*,
00280 const G4DynamicParticle*,
00281 G4double&,G4double&,G4double)
00282 {}
00283
00284
00285
00286 G4double G4VEmModel::Value(const G4MaterialCutsCouple* couple,
00287 const G4ParticleDefinition* p, G4double e)
00288 {
00289 fCurrentCouple = couple;
00290 return e*e*CrossSectionPerVolume(couple->GetMaterial(),p,e,0.0,DBL_MAX);
00291 }
00292
00293
00294
00295 G4double G4VEmModel::MinPrimaryEnergy(const G4Material*,
00296 const G4ParticleDefinition*)
00297 {
00298 return 0.0;
00299 }
00300
00301
00302
00303 G4double G4VEmModel::MaxSecondaryEnergy(const G4ParticleDefinition*,
00304 G4double kineticEnergy)
00305 {
00306 return kineticEnergy;
00307 }
00308
00309
00310
00311 void G4VEmModel::SetupForMaterial(const G4ParticleDefinition*,
00312 const G4Material*, G4double)
00313 {}
00314
00315
00316
00317 void
00318 G4VEmModel::SetParticleChange(G4VParticleChange* p, G4VEmFluctuationModel* f)
00319 {
00320 if(p && pParticleChange != p) { pParticleChange = p; }
00321 flucModel = f;
00322 }
00323
00324
00325
00326 void G4VEmModel::SetCrossSectionTable(G4PhysicsTable* p)
00327 {
00328 if(p != xSectionTable) {
00329 if(xSectionTable) {
00330 xSectionTable->clearAndDestroy();
00331 delete xSectionTable;
00332 }
00333 xSectionTable = p;
00334 }
00335 }
00336
00337