Geant4 Cross Reference |
1 // 2 // ******************************************************************** 3 // * License and Disclaimer * 4 // * * 5 // * The Geant4 software is copyright of the Copyright Holders of * 6 // * the Geant4 Collaboration. It is provided under the terms and * 7 // * conditions of the Geant4 Software License, included in the file * 8 // * LICENSE and available at http://cern.ch/geant4/license . These * 9 // * include a list of copyright holders. * 10 // * * 11 // * Neither the authors of this software system, nor their employing * 12 // * institutes,nor the agencies providing financial support for this * 13 // * work make any representation or warranty, express or implied, * 14 // * regarding this software system or assume any liability for its * 15 // * use. Please see the license in the file LICENSE and URL above * 16 // * for the full disclaimer and the limitation of liability. * 17 // * * 18 // * This code implementation is the result of the scientific and * 19 // * technical work of the GEANT4 collaboration. * 20 // * By using, copying, modifying or distributing the software (or * 21 // * any work based on the software) you agree to acknowledge its * 22 // * use in resulting scientific publications, and indicate your * 23 // * acceptance of all terms of the Geant4 Software license. * 24 // ******************************************************************** 25 // 26 // ------------------------------------------------------------------- 27 // 28 // GEANT4 Class file 29 // 30 // 31 // File name: G4LossTableBuilder 32 // 33 // Author: Vladimir Ivanchenko 34 // 35 // Creation date: 03.01.2002 36 // 37 // Modifications: 38 // 39 // 23-01-03 V.Ivanchenko Cut per region 40 // 21-07-04 V.Ivanchenko Fix problem of range for dedx=0 41 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivanchenko) 42 // 07-12-04 Fix of BuildDEDX table (V.Ivanchenko) 43 // 27-03-06 Add bool options isIonisation (V.Ivanchenko) 44 // 16-01-07 Fill new (not old) DEDX table (V.Ivanchenko) 45 // 12-02-07 Use G4LPhysicsFreeVector for the inverse range table (V.Ivanchenko) 46 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko) 47 // 48 // Class Description: 49 // 50 // ------------------------------------------------------------------- 51 // 52 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 53 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 54 55 #include "G4LossTableBuilder.hh" 56 #include "G4SystemOfUnits.hh" 57 #include "G4PhysicsTable.hh" 58 #include "G4PhysicsLogVector.hh" 59 #include "G4PhysicsTableHelper.hh" 60 #include "G4PhysicsFreeVector.hh" 61 #include "G4ProductionCutsTable.hh" 62 #include "G4MaterialCutsCouple.hh" 63 #include "G4Material.hh" 64 #include "G4VEmModel.hh" 65 #include "G4ParticleDefinition.hh" 66 #include "G4LossTableManager.hh" 67 #include "G4EmParameters.hh" 68 69 G4bool G4LossTableBuilder::baseMatFlag = false; 70 std::vector<G4double>* G4LossTableBuilder::theDensityFactor = nullptr; 71 std::vector<G4int>* G4LossTableBuilder::theDensityIdx = nullptr; 72 std::vector<G4bool>* G4LossTableBuilder::theFlag = nullptr; 73 74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 75 76 G4LossTableBuilder::G4LossTableBuilder(G4bool master) 77 : isInitializer(master) 78 { 79 theParameters = G4EmParameters::Instance(); 80 if (nullptr == theFlag) { 81 theDensityFactor = new std::vector<G4double>; 82 theDensityIdx = new std::vector<G4int>; 83 theFlag = new std::vector<G4bool>; 84 } 85 } 86 87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 88 89 G4LossTableBuilder::~G4LossTableBuilder() 90 { 91 if (isInitializer) { 92 delete theDensityFactor; 93 delete theDensityIdx; 94 delete theFlag; 95 theDensityFactor = nullptr; 96 theDensityIdx = nullptr; 97 theFlag = nullptr; 98 } 99 } 100 101 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 102 103 const std::vector<G4int>* G4LossTableBuilder::GetCoupleIndexes() 104 { 105 return theDensityIdx; 106 } 107 108 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 109 110 const std::vector<G4double>* G4LossTableBuilder::GetDensityFactors() 111 { 112 return theDensityFactor; 113 } 114 115 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 116 117 G4bool G4LossTableBuilder::GetFlag(std::size_t idx) 118 { 119 return (idx < theFlag->size()) ? (*theFlag)[idx] : false; 120 } 121 122 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 123 124 G4bool G4LossTableBuilder::GetBaseMaterialFlag() 125 { 126 return baseMatFlag; 127 } 128 129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 130 131 void 132 G4LossTableBuilder::BuildDEDXTable(G4PhysicsTable* dedxTable, 133 const std::vector<G4PhysicsTable*>& list) 134 { 135 InitialiseBaseMaterials(dedxTable); 136 std::size_t n_processes = list.size(); 137 if(1 >= n_processes) { return; } 138 139 std::size_t nCouples = dedxTable->size(); 140 //G4cout << "Nproc= " << n_processes << " nCouples=" << nCouples << " Nv= " 141 // << dedxTable->size() << G4endl; 142 if(0 >= nCouples) { return; } 143 144 for (std::size_t i=0; i<nCouples; ++i) { 145 auto pv0 = static_cast<G4PhysicsLogVector*>((*(list[0]))[i]); 146 //if (0 == i) G4cout << i << ". pv0=" << pv0 << " t:" << list[0] << G4endl; 147 if(pv0 == nullptr) { continue; } 148 std::size_t npoints = pv0->GetVectorLength(); 149 auto pv = new G4PhysicsLogVector(*pv0); 150 for (std::size_t j=0; j<npoints; ++j) { 151 G4double dedx = 0.0; 152 for (std::size_t k=0; k<n_processes; ++k) { 153 const G4PhysicsVector* pv1 = (*(list[k]))[i]; 154 //if (0 == i) G4cout << " " << k << ". pv1=" << pv1 << " t:" << list[k] << G4endl; 155 dedx += (*pv1)[j]; 156 } 157 pv->PutValue(j, dedx); 158 } 159 if(splineFlag) { pv->FillSecondDerivatives(); } 160 G4PhysicsTableHelper::SetPhysicsVector(dedxTable, i, pv); 161 } 162 //G4cout << "### G4LossTableBuilder::BuildDEDXTable " << G4endl; 163 //G4cout << *dedxTable << G4endl; 164 } 165 166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 167 168 void G4LossTableBuilder::BuildRangeTable(const G4PhysicsTable* dedxTable, 169 G4PhysicsTable* rangeTable) 170 // Build range table from the energy loss table 171 { 172 //G4cout << "### G4LossTableBuilder::BuildRangeTable: DEDX table" << G4endl; 173 //G4cout << *const_cast<G4PhysicsTable*>(dedxTable) << G4endl; 174 const std::size_t nCouples = dedxTable->size(); 175 if(0 >= nCouples) { return; } 176 177 const std::size_t n = 100; 178 const G4double del = 1.0/(G4double)n; 179 180 for (std::size_t i=0; i<nCouples; ++i) { 181 auto pv = static_cast<G4PhysicsLogVector*>((*dedxTable)[i]); 182 if((pv == nullptr) || (isBaseMatActive && !(*theFlag)[i])) { continue; } 183 std::size_t npoints = pv->GetVectorLength(); 184 std::size_t bin0 = 0; 185 G4double elow = pv->Energy(0); 186 G4double ehigh = pv->Energy(npoints-1); 187 G4double dedx1 = (*pv)[0]; 188 189 // protection for specific cases dedx=0 190 if(dedx1 == 0.0) { 191 for (std::size_t k=1; k<npoints; ++k) { 192 ++bin0; 193 elow = pv->Energy(k); 194 dedx1 = (*pv)[k]; 195 if(dedx1 > 0.0) { break; } 196 } 197 npoints -= bin0; 198 } 199 200 // initialisation of a new vector 201 if(npoints < 3) { npoints = 3; } 202 203 delete (*rangeTable)[i]; 204 G4PhysicsLogVector* v; 205 if(0 == bin0) { v = new G4PhysicsLogVector(*pv); } 206 else { v = new G4PhysicsLogVector(elow, ehigh, npoints-1, splineFlag); } 207 208 // assumed dedx proportional to beta 209 G4double energy1 = v->Energy(0); 210 G4double range = 2.*energy1/dedx1; 211 /* 212 G4cout << "New Range vector Npoints=" << v->GetVectorLength() 213 << " coupleIdx=" << i << " spline=" << v->GetSpline() 214 << " Elow=" << v->GetMinEnergy() <<" Ehigh=" << v->GetMinEnergy() 215 << " DEDX(Elow)=" << dedx1 << " R(Elow)=" << range << G4endl; 216 */ 217 v->PutValue(0,range); 218 219 for (std::size_t j=1; j<npoints; ++j) { 220 221 G4double energy2 = v->Energy(j); 222 G4double de = (energy2 - energy1) * del; 223 G4double energy = energy2 + de*0.5; 224 G4double sum = 0.0; 225 std::size_t idx = j - 1; 226 for (std::size_t k=0; k<n; ++k) { 227 energy -= de; 228 dedx1 = pv->Value(energy, idx); 229 if(dedx1 > 0.0) { sum += de/dedx1; } 230 } 231 range += sum; 232 /* 233 if(energy < 10.) 234 G4cout << "j= " << j << " e1= " << energy1 << " e2= " << energy2 235 << " n= " << n << " range=" << range<< G4endl; 236 */ 237 v->PutValue(j,range); 238 energy1 = energy2; 239 } 240 if(splineFlag) { v->FillSecondDerivatives(); } 241 G4PhysicsTableHelper::SetPhysicsVector(rangeTable, i, v); 242 } 243 //G4cout << "### Range table" << G4endl; 244 //G4cout << *rangeTable << G4endl; 245 } 246 247 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 248 249 void 250 G4LossTableBuilder::BuildInverseRangeTable(const G4PhysicsTable* rangeTable, 251 G4PhysicsTable* invRangeTable) 252 // Build inverse range table from the energy loss table 253 { 254 std::size_t nCouples = rangeTable->size(); 255 if(0 >= nCouples) { return; } 256 257 for (std::size_t i=0; i<nCouples; ++i) { 258 G4PhysicsVector* pv = (*rangeTable)[i]; 259 if((pv == nullptr) || (isBaseMatActive && !(*theFlag)[i])) { continue; } 260 std::size_t npoints = pv->GetVectorLength(); 261 262 delete (*invRangeTable)[i]; 263 auto v = new G4PhysicsFreeVector(npoints, splineFlag); 264 265 for (std::size_t j=0; j<npoints; ++j) { 266 G4double e = pv->Energy(j); 267 G4double r = (*pv)[j]; 268 v->PutValues(j,r,e); 269 } 270 if (splineFlag) { v->FillSecondDerivatives(); } 271 v->EnableLogBinSearch(theParameters->NumberForFreeVector()); 272 273 G4PhysicsTableHelper::SetPhysicsVector(invRangeTable, i, v); 274 } 275 //G4cout << "### Inverse range table" << G4endl; 276 //G4cout << *invRangeTable << G4endl; 277 } 278 279 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 280 281 void G4LossTableBuilder::InitialiseBaseMaterials(const G4PhysicsTable* table) 282 { 283 if(!isInitializer) { return; } 284 const G4ProductionCutsTable* theCoupleTable= 285 G4ProductionCutsTable::GetProductionCutsTable(); 286 std::size_t nCouples = theCoupleTable->GetTableSize(); 287 std::size_t nFlags = theFlag->size(); 288 /* 289 G4cout << "### InitialiseBaseMaterials: nCouples=" << nCouples 290 << " nFlags=" << nFlags << " isInit:" << isInitialized 291 << " baseMat:" << baseMatFlag << G4endl; 292 */ 293 // define base material flag 294 if(isBaseMatActive && !baseMatFlag) { 295 for(G4int i=0; i<(G4int)nCouples; ++i) { 296 if(nullptr != theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial()->GetBaseMaterial()) { 297 baseMatFlag = true; 298 isInitialized = false; 299 break; 300 } 301 } 302 } 303 304 if(nFlags != nCouples) { isInitialized = false; } 305 if(isInitialized) { return; } 306 307 // reserve memory 308 theFlag->resize(nCouples, true); 309 theDensityFactor->resize(nCouples,1.0); 310 theDensityIdx->resize(nCouples, 0); 311 312 // define default flag and index of used material cut couple 313 for (G4int i=0; i<(G4int)nCouples; ++i) { 314 (*theFlag)[i] = (nullptr == table) ? true : table->GetFlag(i); 315 (*theDensityIdx)[i] = i; 316 } 317 isInitialized = true; 318 if (!baseMatFlag) { return; } 319 320 // use base materials 321 for (G4int i=0; i<(G4int)nCouples; ++i) { 322 // base material is needed only for a couple which is not 323 // initialised and for which tables will be computed 324 auto couple = theCoupleTable->GetMaterialCutsCouple(i); 325 auto pcuts = couple->GetProductionCuts(); 326 auto mat = couple->GetMaterial(); 327 auto bmat = mat->GetBaseMaterial(); 328 329 // base material exists - find it and check if it can be reused 330 if(nullptr != bmat) { 331 for(G4int j=0; j<(G4int)nCouples; ++j) { 332 if(j == i) { continue; } 333 auto bcouple = theCoupleTable->GetMaterialCutsCouple(j); 334 335 if(bcouple->GetMaterial() == bmat && 336 bcouple->GetProductionCuts() == pcuts) { 337 338 // based couple exist in the same region 339 (*theDensityFactor)[i] = mat->GetDensity()/bmat->GetDensity(); 340 (*theDensityIdx)[i] = j; 341 (*theFlag)[i] = false; 342 343 // ensure that there will no double initialisation 344 (*theDensityFactor)[j] = 1.0; 345 (*theDensityIdx)[j] = j; 346 (*theFlag)[j] = true; 347 break; 348 } 349 } 350 } 351 } 352 } 353 354 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 355 356 G4PhysicsTable* 357 G4LossTableBuilder::BuildTableForModel(G4PhysicsTable* aTable, 358 G4VEmModel* model, 359 const G4ParticleDefinition* part, 360 G4double emin, G4double emax, 361 G4bool spline) 362 { 363 // check input 364 G4PhysicsTable* table = G4PhysicsTableHelper::PreparePhysicsTable(aTable); 365 if (nullptr == table) { return table; } 366 if (aTable != nullptr && aTable != table) { 367 aTable->clearAndDestroy(); 368 delete aTable; 369 } 370 371 InitialiseBaseMaterials(table); 372 G4int nbins = theParameters->NumberOfBinsPerDecade(); 373 374 // Access to materials 375 const G4ProductionCutsTable* theCoupleTable= 376 G4ProductionCutsTable::GetProductionCutsTable(); 377 std::size_t numOfCouples = theCoupleTable->GetTableSize(); 378 /* 379 G4cout << " G4LossTableBuilder::BuildTableForModel Ncouple=" << numOfCouples 380 << " isMaster=" << isInitializer << " model:" << model->GetName() 381 << " " << part->GetParticleName() << G4endl; 382 */ 383 G4PhysicsLogVector* aVector = nullptr; 384 385 for(G4int i=0; i<(G4int)numOfCouples; ++i) { 386 //G4cout << i << ". " << (*theFlag)[i] << " " << table->GetFlag(i) << G4endl; 387 if (table->GetFlag(i)) { 388 389 // create physics vector and fill it 390 auto couple = theCoupleTable->GetMaterialCutsCouple(i); 391 delete (*table)[i]; 392 393 // if start from zero then change the scale 394 const G4Material* mat = couple->GetMaterial(); 395 396 G4double tmin = std::max(emin, model->MinPrimaryEnergy(mat,part)); 397 if(0.0 >= tmin) { tmin = CLHEP::eV; } 398 G4int n = nbins; 399 400 if(tmin >= emax) { 401 aVector = nullptr; 402 } else { 403 n *= G4lrint(std::log10(emax/tmin)); 404 n = std::max(n, 3); 405 aVector = new G4PhysicsLogVector(tmin, emax, n, spline); 406 } 407 408 if(nullptr != aVector) { 409 //G4cout << part->GetParticleName() << " in " << mat->GetName() 410 // << " emin= " << tmin << " emax=" << emax << " n=" << n << G4endl; 411 for(G4int j=0; j<=n; ++j) { 412 G4double e = aVector->Energy(j); 413 G4double y = model->Value(couple, part, e); 414 //G4cout << " " << j << ") E=" << e << " y=" << y << G4endl; 415 aVector->PutValue(j, y); 416 } 417 if(spline) { aVector->FillSecondDerivatives(); } 418 } 419 G4PhysicsTableHelper::SetPhysicsVector(table, i, aVector); 420 } 421 } 422 /* 423 G4cout << "G4LossTableBuilder::BuildTableForModel done for " 424 << part->GetParticleName() << " and "<< model->GetName() 425 << " " << table << G4endl; 426 */ 427 //G4cout << *table << G4endl; 428 return table; 429 } 430 431 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 432