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
00054
00055
00056
00057 #include "G4LossTableBuilder.hh"
00058 #include "G4SystemOfUnits.hh"
00059 #include "G4PhysicsTable.hh"
00060 #include "G4PhysicsLogVector.hh"
00061 #include "G4PhysicsTableHelper.hh"
00062 #include "G4LPhysicsFreeVector.hh"
00063 #include "G4ProductionCutsTable.hh"
00064 #include "G4MaterialCutsCouple.hh"
00065 #include "G4Material.hh"
00066 #include "G4VEmModel.hh"
00067 #include "G4ParticleDefinition.hh"
00068 #include "G4LossTableManager.hh"
00069
00070
00071
00072 G4LossTableBuilder::G4LossTableBuilder()
00073 {
00074 splineFlag = true;
00075 isInitialized = false;
00076
00077 theDensityFactor = new std::vector<G4double>;
00078 theDensityIdx = new std::vector<G4int>;
00079 theFlag = new std::vector<G4bool>;
00080 }
00081
00082
00083
00084 G4LossTableBuilder::~G4LossTableBuilder()
00085 {
00086 delete theDensityFactor;
00087 delete theDensityIdx;
00088 delete theFlag;
00089 }
00090
00091
00092
00093 void
00094 G4LossTableBuilder::BuildDEDXTable(G4PhysicsTable* dedxTable,
00095 const std::vector<G4PhysicsTable*>& list)
00096 {
00097 size_t n_processes = list.size();
00098
00099 if(1 >= n_processes) { return; }
00100
00101 size_t nCouples = dedxTable->size();
00102 if(0 >= nCouples) { return; }
00103
00104 for (size_t i=0; i<nCouples; ++i) {
00105
00106 G4PhysicsLogVector* pv0 = static_cast<G4PhysicsLogVector*>((*(list[0]))[i]);
00107 if(pv0) {
00108 size_t npoints = pv0->GetVectorLength();
00109 G4PhysicsLogVector* pv = new G4PhysicsLogVector(*pv0);
00110
00111 pv->SetSpline(splineFlag);
00112 for (size_t j=0; j<npoints; ++j) {
00113 G4double dedx = 0.0;
00114 for (size_t k=0; k<n_processes; ++k) {
00115 G4PhysicsVector* pv1 = (*(list[k]))[i];
00116 dedx += (*pv1)[j];
00117 }
00118 pv->PutValue(j, dedx);
00119 }
00120 if(splineFlag) { pv->FillSecondDerivatives(); }
00121 G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv);
00122 }
00123 }
00124 }
00125
00126
00127
00128 void G4LossTableBuilder::BuildRangeTable(const G4PhysicsTable* dedxTable,
00129 G4PhysicsTable* rangeTable,
00130 G4bool isIonisation)
00131
00132 {
00133 size_t nCouples = dedxTable->size();
00134 if(0 >= nCouples) { return; }
00135
00136 size_t n = 100;
00137 G4double del = 1.0/(G4double)n;
00138
00139 for (size_t i=0; i<nCouples; ++i) {
00140 if(isIonisation) {
00141 if( !(*theFlag)[i] ) { continue; }
00142 }
00143 G4PhysicsLogVector* pv = static_cast<G4PhysicsLogVector*>((*dedxTable)[i]);
00144 size_t npoints = pv->GetVectorLength();
00145 size_t bin0 = 0;
00146 G4double elow = pv->Energy(0);
00147 G4double ehigh = pv->Energy(npoints-1);
00148 G4double dedx1 = (*pv)[0];
00149
00150
00151
00152
00153 if(dedx1 == 0.0) {
00154 for (size_t k=1; k<npoints; ++k) {
00155 bin0++;
00156 elow = pv->Energy(k);
00157 dedx1 = (*pv)[k];
00158 if(dedx1 > 0.0) { break; }
00159 }
00160 npoints -= bin0;
00161 }
00162
00163
00164
00165
00166
00167 if(npoints < 2) { npoints = 2; }
00168
00169 delete (*rangeTable)[i];
00170 G4PhysicsLogVector* v;
00171 if(0 == bin0) { v = new G4PhysicsLogVector(*pv); }
00172 else { v = new G4PhysicsLogVector(elow, ehigh, npoints-1); }
00173
00174
00175 if(2 == npoints) {
00176 v->PutValue(0,1000.);
00177 v->PutValue(1,2000.);
00178 G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
00179 return;
00180 }
00181 v->SetSpline(splineFlag);
00182
00183
00184 G4double energy1 = v->Energy(0);
00185 G4double range = 2.*energy1/dedx1;
00186
00187 v->PutValue(0,range);
00188
00189 for (size_t j=1; j<npoints; ++j) {
00190
00191 G4double energy2 = v->Energy(j);
00192 G4double de = (energy2 - energy1) * del;
00193 G4double energy = energy2 + de*0.5;
00194 G4double sum = 0.0;
00195
00196
00197 for (size_t k=0; k<n; ++k) {
00198 energy -= de;
00199 dedx1 = pv->Value(energy);
00200 if(dedx1 > 0.0) { sum += de/dedx1; }
00201 }
00202 range += sum;
00203 v->PutValue(j,range);
00204 energy1 = energy2;
00205 }
00206 if(splineFlag) { v->FillSecondDerivatives(); }
00207 G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v);
00208 }
00209 }
00210
00211
00212
00213 void G4LossTableBuilder::BuildInverseRangeTable(const G4PhysicsTable* rangeTable,
00214 G4PhysicsTable* invRangeTable,
00215 G4bool isIonisation)
00216
00217 {
00218 size_t nCouples = rangeTable->size();
00219 if(0 >= nCouples) { return; }
00220
00221 for (size_t i=0; i<nCouples; ++i) {
00222
00223 if(isIonisation) {
00224 if( !(*theFlag)[i] ) { continue; }
00225 }
00226 G4PhysicsVector* pv = (*rangeTable)[i];
00227 size_t npoints = pv->GetVectorLength();
00228 G4double rlow = (*pv)[0];
00229 G4double rhigh = (*pv)[npoints-1];
00230
00231 delete (*invRangeTable)[i];
00232 G4LPhysicsFreeVector* v = new G4LPhysicsFreeVector(npoints,rlow,rhigh);
00233 v->SetSpline(splineFlag);
00234
00235 for (size_t j=0; j<npoints; ++j) {
00236 G4double e = pv->Energy(j);
00237 G4double r = (*pv)[j];
00238 v->PutValues(j,r,e);
00239 }
00240 if(splineFlag) { v->FillSecondDerivatives(); }
00241
00242 G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v);
00243 }
00244 }
00245
00246
00247
00248 void
00249 G4LossTableBuilder::InitialiseBaseMaterials(G4PhysicsTable* table)
00250 {
00251 size_t nCouples = table->size();
00252 size_t nFlags = theFlag->size();
00253
00254 if(nCouples == nFlags && isInitialized) { return; }
00255
00256 isInitialized = true;
00257
00258
00259
00260
00261
00262 const G4ProductionCutsTable* theCoupleTable=
00263 G4ProductionCutsTable::GetProductionCutsTable();
00264
00265
00266
00267
00268
00269
00270
00271
00272
00273
00274
00275 if(nFlags < nCouples) {
00276 for(size_t i=nFlags; i<nCouples; ++i) { theDensityFactor->push_back(1.0); }
00277 for(size_t i=nFlags; i<nCouples; ++i) { theDensityIdx->push_back(-1); }
00278 for(size_t i=nFlags; i<nCouples; ++i) { theFlag->push_back(true); }
00279 }
00280 for(size_t i=0; i<nCouples; ++i) {
00281
00282
00283
00284 (*theFlag)[i] = table->GetFlag(i);
00285 if ((*theDensityIdx)[i] < 0) {
00286 (*theDensityIdx)[i] = i;
00287 const G4MaterialCutsCouple* couple =
00288 theCoupleTable->GetMaterialCutsCouple(i);
00289 const G4ProductionCuts* pcuts = couple->GetProductionCuts();
00290 const G4Material* mat = couple->GetMaterial();
00291 const G4Material* bmat = mat->GetBaseMaterial();
00292
00293
00294 if(bmat) {
00295 for(size_t j=0; j<nCouples; ++j) {
00296
00297 if(j == i) { continue; }
00298 const G4MaterialCutsCouple* bcouple =
00299 theCoupleTable->GetMaterialCutsCouple(j);
00300
00301 if(bcouple->GetMaterial() == bmat &&
00302 bcouple->GetProductionCuts() == pcuts) {
00303
00304
00305 (*theDensityIdx)[i] = j;
00306 (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity();
00307 (*theFlag)[i] = false;
00308
00309
00310 (*theDensityIdx)[j] = j;
00311 (*theDensityFactor)[j] = 1.0;
00312 (*theFlag)[j] = true;
00313 break;
00314 }
00315 }
00316 }
00317 }
00318 }
00319
00320
00321
00322
00323
00324
00325
00326
00327
00328
00329 }
00330
00331
00332
00333 void G4LossTableBuilder::InitialiseCouples()
00334 {
00335 isInitialized = true;
00336
00337
00338 const G4ProductionCutsTable* theCoupleTable=
00339 G4ProductionCutsTable::GetProductionCutsTable();
00340 size_t nCouples = theCoupleTable->GetTableSize();
00341
00342
00343
00344 theDensityFactor->resize(nCouples, 1.0);
00345 theDensityIdx->resize(nCouples, -1);
00346 theFlag->resize(nCouples, true);
00347
00348 for(size_t i=0; i<nCouples; ++i) {
00349
00350
00351
00352 if ((*theDensityIdx)[i] < 0) {
00353 (*theDensityIdx)[i] = i;
00354 const G4MaterialCutsCouple* couple =
00355 theCoupleTable->GetMaterialCutsCouple(i);
00356 const G4ProductionCuts* pcuts = couple->GetProductionCuts();
00357 const G4Material* mat = couple->GetMaterial();
00358 const G4Material* bmat = mat->GetBaseMaterial();
00359
00360
00361 if(bmat) {
00362 for(size_t j=0; j<nCouples; ++j) {
00363
00364 if(j == i) { continue; }
00365 const G4MaterialCutsCouple* bcouple =
00366 theCoupleTable->GetMaterialCutsCouple(j);
00367
00368 if(bcouple->GetMaterial() == bmat &&
00369 bcouple->GetProductionCuts() == pcuts) {
00370
00371
00372 (*theDensityIdx)[i] = j;
00373 (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity();
00374 (*theFlag)[i] = false;
00375
00376
00377 (*theDensityIdx)[j] = j;
00378 (*theDensityFactor)[j] = 1.0;
00379 (*theFlag)[j] = true;
00380 break;
00381 }
00382 }
00383 }
00384 }
00385 }
00386
00387
00388
00389
00390
00391
00392
00393
00394
00395 }
00396
00397
00398
00399 G4PhysicsTable*
00400 G4LossTableBuilder::BuildTableForModel(G4PhysicsTable* aTable,
00401 G4VEmModel* model,
00402 const G4ParticleDefinition* part,
00403 G4double emin, G4double emax,
00404 G4bool spline)
00405 {
00406
00407 G4PhysicsTable* table = G4PhysicsTableHelper::PreparePhysicsTable(aTable);
00408 if(!table) { return table; }
00409 if(emin >= emax) {
00410 table->clearAndDestroy();
00411 delete table;
00412 table = 0;
00413 return table;
00414 }
00415 InitialiseBaseMaterials(table);
00416
00417 G4int nbins = G4int(std::log10(emax/emin) + 0.5)
00418 *G4LossTableManager::Instance()->GetNumberOfBinsPerDecade();
00419 if(nbins < 3) { nbins = 3; }
00420
00421
00422 const G4ProductionCutsTable* theCoupleTable=
00423 G4ProductionCutsTable::GetProductionCutsTable();
00424 size_t numOfCouples = theCoupleTable->GetTableSize();
00425
00426 G4PhysicsLogVector* aVector = 0;
00427 G4PhysicsLogVector* bVector = 0;
00428
00429 for(size_t i=0; i<numOfCouples; ++i) {
00430
00431
00432
00433 if (GetFlag(i)) {
00434
00435
00436 const G4MaterialCutsCouple* couple =
00437 theCoupleTable->GetMaterialCutsCouple(i);
00438 delete (*table)[i];
00439
00440
00441
00442 const G4Material* mat = couple->GetMaterial();
00443
00444 G4double tmin = std::max(emin,model->MinPrimaryEnergy(mat,part));
00445 if(0.0 >= tmin) { tmin = eV; }
00446 G4int n = nbins + 1;
00447
00448 if(tmin >= emax) {
00449 aVector = 0;
00450 } else if(tmin > emin) {
00451 G4int bin = nbins*G4int(std::log10(emax/tmin) + 0.5);
00452 if(bin < 3) { bin = 3; }
00453 n = bin + 1;
00454 aVector = new G4PhysicsLogVector(tmin, emax, bin);
00455
00456 } else if(!bVector) {
00457 aVector = new G4PhysicsLogVector(emin, emax, nbins);
00458 bVector = aVector;
00459
00460 } else {
00461 aVector = new G4PhysicsLogVector(*bVector);
00462 }
00463
00464 if(aVector) {
00465 aVector->SetSpline(spline);
00466 for(G4int j=0; j<n; ++j) {
00467 aVector->PutValue(j, model->Value(couple, part, aVector->Energy(j)));
00468 }
00469 if(spline) { aVector->FillSecondDerivatives(); }
00470 }
00471 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector);
00472 }
00473 }
00474
00475
00476
00477
00478
00479 return table;
00480 }
00481
00482
00483