#include <G4LossTableBuilder.hh>
Public Member Functions | |
G4LossTableBuilder () | |
virtual | ~G4LossTableBuilder () |
void | BuildDEDXTable (G4PhysicsTable *dedxTable, const std::vector< G4PhysicsTable * > &) |
void | BuildRangeTable (const G4PhysicsTable *dedxTable, G4PhysicsTable *rangeTable, G4bool isIonisation=false) |
void | BuildInverseRangeTable (const G4PhysicsTable *rangeTable, G4PhysicsTable *invRangeTable, G4bool isIonisation=false) |
G4PhysicsTable * | BuildTableForModel (G4PhysicsTable *table, G4VEmModel *model, const G4ParticleDefinition *, G4double emin, G4double emax, G4bool spline) |
void | InitialiseBaseMaterials (G4PhysicsTable *table) |
const std::vector< G4int > * | GetCoupleIndexes () |
const std::vector< G4double > * | GetDensityFactors () |
G4bool | GetFlag (size_t idx) const |
void | SetSplineFlag (G4bool flag) |
void | SetInitialisationFlag (G4bool flag) |
Definition at line 61 of file G4LossTableBuilder.hh.
G4LossTableBuilder::G4LossTableBuilder | ( | ) |
Definition at line 72 of file G4LossTableBuilder.cc.
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 }
G4LossTableBuilder::~G4LossTableBuilder | ( | ) | [virtual] |
void G4LossTableBuilder::BuildDEDXTable | ( | G4PhysicsTable * | dedxTable, | |
const std::vector< G4PhysicsTable * > & | ||||
) |
Definition at line 94 of file G4LossTableBuilder.cc.
References G4PhysicsVector::FillSecondDerivatives(), G4PhysicsVector::PutValue(), G4PhysicsTableHelper::SetPhysicsVector(), and G4PhysicsVector::SetSpline().
00096 { 00097 size_t n_processes = list.size(); 00098 //G4cout << "Nproc= " << n_processes << " Ncoup= " << dedxTable->size() << G4endl; 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 // if ((*theFlag)[i]) { 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 // pv = new G4PhysicsLogVector(elow, ehigh, npoints-1); 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 }
void G4LossTableBuilder::BuildInverseRangeTable | ( | const G4PhysicsTable * | rangeTable, | |
G4PhysicsTable * | invRangeTable, | |||
G4bool | isIonisation = false | |||
) |
Definition at line 213 of file G4LossTableBuilder.cc.
References G4PhysicsVector::Energy(), G4PhysicsVector::FillSecondDerivatives(), G4PhysicsVector::GetVectorLength(), G4LPhysicsFreeVector::PutValues(), G4PhysicsTableHelper::SetPhysicsVector(), and G4PhysicsVector::SetSpline().
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 }
void G4LossTableBuilder::BuildRangeTable | ( | const G4PhysicsTable * | dedxTable, | |
G4PhysicsTable * | rangeTable, | |||
G4bool | isIonisation = false | |||
) |
Definition at line 128 of file G4LossTableBuilder.cc.
References G4PhysicsVector::Energy(), G4PhysicsVector::FillSecondDerivatives(), G4PhysicsVector::GetVectorLength(), CLHEP::detail::n, G4PhysicsVector::PutValue(), G4PhysicsTableHelper::SetPhysicsVector(), G4PhysicsVector::SetSpline(), and G4PhysicsVector::Value().
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 //G4cout << "i= " << i << "npoints= " << npoints << " dedx1= " << dedx1 << G4endl; 00151 00152 // protection for specific cases dedx=0 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 //G4cout<<"New Range vector" << G4endl; 00163 //G4cout<<"nbins= "<<npoints-1<<" elow= "<<elow<<" ehigh= "<<ehigh 00164 // <<" bin0= " << bin0 <<G4endl; 00165 00166 // initialisation of a new vector 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 // dedx is exact zero cannot build range table 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 // assumed dedx proportional to beta 00184 G4double energy1 = v->Energy(0); 00185 G4double range = 2.*energy1/dedx1; 00186 //G4cout << "range0= " << range << G4endl; 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 //G4cout << "j= " << j << " e1= " << energy1 << " e2= " << energy2 00196 // << " n= " << n << G4endl; 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 }
G4PhysicsTable * G4LossTableBuilder::BuildTableForModel | ( | G4PhysicsTable * | table, | |
G4VEmModel * | model, | |||
const G4ParticleDefinition * | , | |||
G4double | emin, | |||
G4double | emax, | |||
G4bool | spline | |||
) |
Definition at line 400 of file G4LossTableBuilder.cc.
References G4PhysicsTable::clearAndDestroy(), G4PhysicsVector::Energy(), G4PhysicsVector::FillSecondDerivatives(), GetFlag(), G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4LossTableManager::GetNumberOfBinsPerDecade(), G4ProductionCutsTable::GetProductionCutsTable(), G4ProductionCutsTable::GetTableSize(), InitialiseBaseMaterials(), G4LossTableManager::Instance(), G4VEmModel::MinPrimaryEnergy(), CLHEP::detail::n, G4PhysicsTableHelper::PreparePhysicsTable(), G4PhysicsVector::PutValue(), G4PhysicsTableHelper::SetPhysicsVector(), G4PhysicsVector::SetSpline(), and G4VEmModel::Value().
00405 { 00406 // check input 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 // Access to materials 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 //G4cout<< "i= " << i << " Flag= " << GetFlag(i) << G4endl; 00432 00433 if (GetFlag(i)) { 00434 00435 // create physics vector and fill it 00436 const G4MaterialCutsCouple* couple = 00437 theCoupleTable->GetMaterialCutsCouple(i); 00438 delete (*table)[i]; 00439 00440 // if start from zero then change the scale 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 G4cout << "G4LossTableBuilder::BuildTableForModel done for " 00476 << part->GetParticleName() << " and "<< model->GetName() 00477 << " " << table << G4endl; 00478 */ 00479 return table; 00480 }
const std::vector< G4int > * G4LossTableBuilder::GetCoupleIndexes | ( | ) | [inline] |
Definition at line 123 of file G4LossTableBuilder.hh.
Referenced by G4VEnergyLossProcess::PreparePhysicsTable(), and G4VEmProcess::PreparePhysicsTable().
00124 { 00125 if(theDensityIdx->size() == 0) { InitialiseCouples(); } 00126 return theDensityIdx; 00127 }
const std::vector< G4double > * G4LossTableBuilder::GetDensityFactors | ( | ) | [inline] |
Definition at line 130 of file G4LossTableBuilder.hh.
Referenced by G4VEnergyLossProcess::PreparePhysicsTable(), and G4VEmProcess::PreparePhysicsTable().
00131 { 00132 if(theDensityIdx->size() == 0) { InitialiseCouples(); } 00133 return theDensityFactor; 00134 }
G4bool G4LossTableBuilder::GetFlag | ( | size_t | idx | ) | const [inline] |
Definition at line 136 of file G4LossTableBuilder.hh.
Referenced by G4VEnergyLossProcess::BuildDEDXTable(), G4VEnergyLossProcess::BuildLambdaTable(), and BuildTableForModel().
void G4LossTableBuilder::InitialiseBaseMaterials | ( | G4PhysicsTable * | table | ) |
Definition at line 249 of file G4LossTableBuilder.cc.
References G4Material::GetBaseMaterial(), G4Material::GetDensity(), G4PhysicsTable::GetFlag(), G4MaterialCutsCouple::GetMaterial(), G4ProductionCutsTable::GetMaterialCutsCouple(), G4MaterialCutsCouple::GetProductionCuts(), and G4ProductionCutsTable::GetProductionCutsTable().
Referenced by BuildTableForModel(), G4VEnergyLossProcess::PreparePhysicsTable(), and G4VEmProcess::PreparePhysicsTable().
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 //G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials Ncouples= " 00259 // << nCouples << " FlagSize= " << nFlags << G4endl; 00260 00261 // variable density check 00262 const G4ProductionCutsTable* theCoupleTable= 00263 G4ProductionCutsTable::GetProductionCutsTable(); 00264 00265 /* 00266 for(size_t i=0; i<nFlags; ++i) { 00267 G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i] 00268 << " tableFlag= " << table->GetFlag(i) << " " 00269 << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName() 00270 << G4endl; 00271 } 00272 */ 00273 00274 // expand vectors 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 // base material is needed only for a couple which is not 00283 // initialised and for which tables will be computed 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 // base material exists - find it and check if it can be reused 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 // based couple exist in the same region 00305 (*theDensityIdx)[i] = j; 00306 (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity(); 00307 (*theFlag)[i] = false; 00308 00309 // ensure that there will no double initialisation 00310 (*theDensityIdx)[j] = j; 00311 (*theDensityFactor)[j] = 1.0; 00312 (*theFlag)[j] = true; 00313 break; 00314 } 00315 } 00316 } 00317 } 00318 } 00319 /* 00320 for(size_t i=0; i<nCouples; ++i) { 00321 G4cout << "CoupleIdx= " << i << " Flag= " << (*theFlag)[i] 00322 << " TableFlag= " << table->GetFlag(i) << " " 00323 << theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetName() 00324 << G4endl; 00325 } 00326 G4cout << "%%%%%% G4LossTableBuilder::InitialiseBaseMaterials end" 00327 << G4endl; 00328 */ 00329 }
void G4LossTableBuilder::SetInitialisationFlag | ( | G4bool | flag | ) | [inline] |
Definition at line 146 of file G4LossTableBuilder.hh.
Referenced by G4LossTableManager::PreparePhysicsTable().
void G4LossTableBuilder::SetSplineFlag | ( | G4bool | flag | ) | [inline] |
Definition at line 141 of file G4LossTableBuilder.hh.
Referenced by G4LossTableManager::SetSplineFlag().