- Timestamp:
- Apr 6, 2009, 12:30:29 PM (15 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/cross_sections/src/G4CrossSectionDataStore.cc
r819 r962 24 24 // ******************************************************************** 25 25 // 26 // $Id: G4CrossSectionDataStore.cc,v 1.16 2009/01/24 11:54:47 vnivanch Exp $ 27 // GEANT4 tag $Name: geant4-09-02-ref-02 $ 28 // 29 // ------------------------------------------------------------------- 30 // 31 // GEANT4 Class file 32 // 33 // 34 // File name: G4CrossSectionDataStore 35 // 36 // 37 // Modifications: 38 // 23.01.2009 V.Ivanchenko add destruction of data sets 39 // 26 40 // 27 41 … … 31 45 #include "G4HadTmpUtil.hh" 32 46 #include "Randomize.hh" 33 47 #include "G4Nucleus.hh" 48 49 G4CrossSectionDataStore::G4CrossSectionDataStore() : 50 NDataSetList(0), verboseLevel(0) 51 {} 52 53 G4CrossSectionDataStore::~G4CrossSectionDataStore() 54 {} 34 55 35 56 G4double … … 140 161 141 162 142 std::pair<G4double, G4double> 143 G4CrossSectionDataStore::SelectRandomIsotope(const G4DynamicParticle* particle, 144 const G4Material* aMaterial) 163 G4Element* G4CrossSectionDataStore::SampleZandA(const G4DynamicParticle* particle, 164 const G4Material* aMaterial, 165 G4Nucleus& target) 166 { 167 G4double aTemp = aMaterial->GetTemperature(); 168 const G4int nElements = aMaterial->GetNumberOfElements(); 169 const G4ElementVector* theElementVector = aMaterial->GetElementVector(); 170 G4Element* anElement = (*theElementVector)[0]; 171 G4VCrossSectionDataSet* inCharge; 172 G4int i; 173 174 // compounds 175 if(1 < nElements) { 176 G4double* xsec = new G4double [nElements]; 177 const G4double* theAtomsPerVolumeVector = aMaterial->GetVecNbOfAtomsPerVolume(); 178 G4double cross = 0.0; 179 for(i=0; i<nElements; i++) { 180 anElement= (*theElementVector)[i]; 181 inCharge = whichDataSetInCharge(particle, anElement); 182 cross += theAtomsPerVolumeVector[i]* 183 inCharge->GetCrossSection(particle, anElement, aTemp); 184 xsec[i] = cross; 185 } 186 cross *= G4UniformRand(); 187 188 for(i=0; i<nElements; i++) { 189 if( cross <= xsec[i] ) { 190 anElement = (*theElementVector)[i]; 191 break; 192 } 193 } 194 delete [] xsec; 195 } 196 197 // element have been selected 198 inCharge = whichDataSetInCharge(particle, anElement); 199 G4double ZZ = anElement->GetZ(); 200 G4double AA; 201 202 // Collect abundance weighted cross sections and A values for each isotope 203 // in each element 204 205 const G4int nIsoPerElement = anElement->GetNumberOfIsotopes(); 206 207 // user-defined isotope abundances 208 if (0 < nIsoPerElement) { 209 G4IsotopeVector* isoVector = anElement->GetIsotopeVector(); 210 AA = G4double((*isoVector)[0]->GetN()); 211 if(1 < nIsoPerElement) { 212 213 G4double* xsec = new G4double [nIsoPerElement]; 214 G4double iso_xs = 0.0; 215 G4double cross = 0.0; 216 217 G4double* abundVector = anElement->GetRelativeAbundanceVector(); 218 G4bool elementXS = false; 219 for (i = 0; i<nIsoPerElement; i++) { 220 if (inCharge->IsZAApplicable(particle, ZZ, G4double((*isoVector)[i]->GetN()))) { 221 iso_xs = inCharge->GetIsoCrossSection(particle, (*isoVector)[i], aTemp); 222 } else if (elementXS == false) { 223 iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp); 224 elementXS = true; 225 } 226 227 cross += abundVector[i]*iso_xs; 228 xsec[i] = cross; 229 } 230 cross *= G4UniformRand(); 231 for (i = 0; i<nIsoPerElement; i++) { 232 if(cross <= xsec[i]) { 233 AA = G4double((*isoVector)[i]->GetN()); 234 break; 235 } 236 } 237 delete [] xsec; 238 } 239 // natural abundances 240 } else { 241 242 G4StableIsotopes theDefaultIsotopes; 243 G4int Z = G4int(ZZ + 0.5); 244 const G4int nIso = theDefaultIsotopes.GetNumberOfIsotopes(Z); 245 G4int index = theDefaultIsotopes.GetFirstIsotope(Z); 246 AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index)); 247 248 if(1 < nIso) { 249 250 G4double* xsec = new G4double [nIso]; 251 G4double iso_xs = 0.0; 252 G4double cross = 0.0; 253 G4bool elementXS= false; 254 255 for (i = 0; i<nIso; i++) { 256 AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index+i)); 257 if (inCharge->IsZAApplicable(particle, ZZ, AA )) { 258 iso_xs = inCharge->GetIsoZACrossSection(particle, ZZ, AA, aTemp); 259 } else if (elementXS == false) { 260 iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp); 261 elementXS = true; 262 } 263 cross += theDefaultIsotopes.GetAbundance(index+i)*iso_xs; 264 xsec[i] = cross; 265 } 266 cross *= G4UniformRand(); 267 for (i = 0; i<nIso; i++) { 268 if(cross <= xsec[i]) { 269 AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index+i)); 270 break; 271 } 272 } 273 delete [] xsec; 274 } 275 } 276 //G4cout << "XS: " << particle->GetDefinition()->GetParticleName() 277 // << " e(GeV)= " << particle->GetKineticEnergy()/GeV 278 // << " in " << aMaterial->GetName() 279 // << " ZZ= " << ZZ << " AA= " << AA << " " << anElement->GetName() << G4endl; 280 281 target.SetParameters(AA, ZZ); 282 return anElement; 283 } 284 285 /* 286 G4Element* G4CrossSectionDataStore::SampleZandA(const G4DynamicParticle* particle, 287 const G4Material* aMaterial, 288 G4Nucleus& target) 145 289 { 146 290 static G4StableIsotopes theDefaultIsotopes; // natural abundances and … … 229 373 // cross section is zero 230 374 375 anElement = (*theElementVector)[0]; 376 G4double ZZ = anElement->GetZ(); 377 G4double AA = AvaluesPerElement[0][0]; 378 if (crossSectionTotal != 0.) { 379 G4double random = G4UniformRand(); 380 for(G4int i=0; i < nElements; i++) { 381 if(i!=0) runningSum[i] += runningSum[i-1]; 382 if(random <= runningSum[i]/crossSectionTotal) { 383 anElement = (*theElementVector)[i]; 384 ZZ = anElement->GetZ(); 385 386 // Compare random number to running sum over isotope xc to choose A 387 388 G4int nIso = awicsPerElement[i].size(); 389 G4double* running = new G4double[nIso]; 390 for (G4int j=0; j < nIso; j++) { 391 running[j] = awicsPerElement[i][j]; 392 if(j!=0) running[j] += running[j-1]; 393 } 394 395 G4double trial = G4UniformRand(); 396 for (G4int j=0; j < nIso; j++) { 397 AA = AvaluesPerElement[i][j]; 398 if (trial <= running[j]/running[nIso-1]) break; 399 } 400 delete [] running; 401 break; 402 } 403 } 404 } 405 406 //G4cout << "XS: " << particle->GetDefinition()->GetParticleName() 407 // << " e(GeV)= " << particle->GetKineticEnergy()/GeV 408 // << " in " << aMaterial->GetName() 409 // << " ZZ= " << ZZ << " AA= " << AA << " " << anElement->GetName() << G4endl; 410 411 target.SetParameters(AA, ZZ); 412 return anElement; 413 } 414 415 std::pair<G4double, G4double> 416 G4CrossSectionDataStore::SelectRandomIsotope(const G4DynamicParticle* particle, 417 const G4Material* aMaterial) 418 { 419 static G4StableIsotopes theDefaultIsotopes; // natural abundances and 420 // stable isotopes 421 G4double aTemp = aMaterial->GetTemperature(); 422 G4int nElements = aMaterial->GetNumberOfElements(); 423 const G4ElementVector* theElementVector = aMaterial->GetElementVector(); 424 const G4double* theAtomsPerVolumeVector = aMaterial->GetVecNbOfAtomsPerVolume(); 425 std::vector<std::vector<G4double> > awicsPerElement; 426 std::vector<std::vector<G4double> > AvaluesPerElement; 427 G4Element* anElement; 428 429 // Collect abundance weighted cross sections and A values for each isotope 430 // in each element 431 432 for (G4int i = 0; i < nElements; i++) { 433 anElement = (*theElementVector)[i]; 434 G4int nIsoPerElement = anElement->GetNumberOfIsotopes(); 435 std::vector<G4double> isoholder; 436 std::vector<G4double> aholder; 437 G4double iso_xs = DBL_MIN; 438 439 if (nIsoPerElement) { // user-defined isotope abundances 440 G4IsotopeVector* isoVector = anElement->GetIsotopeVector(); 441 G4double* abundVector = anElement->GetRelativeAbundanceVector(); 442 G4VCrossSectionDataSet* inCharge = whichDataSetInCharge(particle, anElement); 443 G4bool elementXS = false; 444 for (G4int j = 0; j < nIsoPerElement; j++) { 445 if (inCharge->IsZAApplicable(particle, (*isoVector)[j]->GetZ(), 446 (*isoVector)[j]->GetN() ) ) { 447 iso_xs = inCharge->GetIsoCrossSection(particle, (*isoVector)[j], aTemp); 448 } else if (elementXS == false) { 449 iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp); 450 elementXS = true; 451 } 452 453 isoholder.push_back(abundVector[j]*iso_xs); 454 aholder.push_back(G4double((*isoVector)[j]->GetN())); 455 } 456 457 } else { // natural abundances 458 G4int ZZ = G4lrint(anElement->GetZ()); 459 nIsoPerElement = theDefaultIsotopes.GetNumberOfIsotopes(ZZ); 460 G4int index = theDefaultIsotopes.GetFirstIsotope(ZZ); 461 G4double AA; 462 G4double abundance; 463 464 G4VCrossSectionDataSet* inCharge = whichDataSetInCharge(particle, anElement); 465 G4bool elementXS = false; 466 467 for (G4int j = 0; j < nIsoPerElement; j++) { 468 AA = G4double(theDefaultIsotopes.GetIsotopeNucleonCount(index+j)); 469 aholder.push_back(AA); 470 if (inCharge->IsZAApplicable(particle, G4double(ZZ), AA) ) { 471 iso_xs = inCharge->GetIsoZACrossSection(particle, G4double(ZZ), AA, aTemp); 472 } else if (elementXS == false) { 473 iso_xs = inCharge->GetCrossSection(particle, anElement, aTemp); 474 elementXS = true; 475 } 476 abundance = theDefaultIsotopes.GetAbundance(index+j)/100.0; 477 isoholder.push_back(abundance*iso_xs); 478 } 479 } 480 481 awicsPerElement.push_back(isoholder); 482 AvaluesPerElement.push_back(aholder); 483 } 484 485 // Calculate running sums for isotope selection 486 487 G4double crossSectionTotal = 0; 488 G4double xSectionPerElement; 489 std::vector<G4double> runningSum; 490 491 for (G4int i=0; i < nElements; i++) { 492 xSectionPerElement = 0; 493 for (G4int j=0; j < G4int(awicsPerElement[i].size()); j++) 494 xSectionPerElement += awicsPerElement[i][j]; 495 runningSum.push_back(theAtomsPerVolumeVector[i]*xSectionPerElement); 496 crossSectionTotal += runningSum[i]; 497 } 498 499 // Compare random number to running sum over element xc to choose Z 500 501 // Initialize Z and A to first element and first isotope in case 502 // cross section is zero 503 231 504 G4double ZZ = (*theElementVector)[0]->GetZ(); 232 505 G4double AA = AvaluesPerElement[0][0]; … … 259 532 return std::pair<G4double, G4double>(ZZ, AA); 260 533 } 261 534 */ 262 535 263 536 void 264 537 G4CrossSectionDataStore::AddDataSet(G4VCrossSectionDataSet* aDataSet) 265 538 { 266 if (NDataSetList == NDataSetMax) { 267 G4cout << "WARNING: G4CrossSectionDataStore::AddDataSet: "<<G4endl; 268 G4cout << " reached maximum number of data sets"; 269 G4cout << " data set not added !!!!!!!!!!!!!!!!"; 270 return; 271 } 272 DataSetList[NDataSetList] = aDataSet; 273 NDataSetList++; 274 } 275 539 DataSetList.push_back(aDataSet); 540 NDataSetList++; 541 } 276 542 277 543 void … … 279 545 BuildPhysicsTable(const G4ParticleDefinition& aParticleType) 280 546 { 281 if (NDataSetList == 0) 282 { 283 G4Exception("G4CrossSectionDataStore", "007", FatalException, 284 "BuildPhysicsTable: no data sets registered"); 285 return; 286 } 287 for (G4int i = NDataSetList-1; i >= 0; i--) { 547 if(NDataSetList > 0) { 548 for (G4int i=0; i<NDataSetList; i++) { 288 549 DataSetList[i]->BuildPhysicsTable(aParticleType); 289 } 550 } 551 } 290 552 } 291 553 … … 295 557 DumpPhysicsTable(const G4ParticleDefinition& aParticleType) 296 558 { 297 if (NDataSetList == 0) { 298 G4cout << "WARNING - G4CrossSectionDataStore::DumpPhysicsTable: no data sets registered"<<G4endl; 299 return; 300 } 301 for (G4int i = NDataSetList-1; i >= 0; i--) { 302 DataSetList[i]->DumpPhysicsTable(aParticleType); 303 } 304 } 559 if (NDataSetList == 0) { 560 G4cout << "WARNING - G4CrossSectionDataStore::DumpPhysicsTable: " 561 << " no data sets registered"<<G4endl; 562 return; 563 } 564 for (G4int i = NDataSetList-1; i >= 0; i--) { 565 DataSetList[i]->DumpPhysicsTable(aParticleType); 566 } 567 }
Note: See TracChangeset
for help on using the changeset viewer.