- Timestamp:
- Dec 22, 2010, 3:52:27 PM (14 years ago)
- Location:
- trunk/source/processes/hadronic/models/radioactive_decay/src
- Files:
-
- 2 edited
Legend:
- Unmodified
- Added
- Removed
-
trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecay.cc
r1340 r1347 121 121 // Constructor 122 122 // 123 G4RadioactiveDecay::G4RadioactiveDecay 124 (const G4String& processName) 125 :G4VRestDiscreteProcess(processName, fDecay), HighestBinValue(10.0), 126 LowestBinValue(1.0e-3), TotBin(200), verboseLevel(0) 123 G4RadioactiveDecay::G4RadioactiveDecay (const G4String& processName) 124 :G4VRestDiscreteProcess(processName, fDecay), HighestBinValue(10.0), 125 LowestBinValue(1.0e-3), TotBin(200), verboseLevel(0) 127 126 { 128 127 #ifdef G4VERBOSE … … 136 135 aPhysicsTable = 0; 137 136 pParticleChange = &fParticleChangeForRadDecay; 138 137 139 138 // 140 139 // Now register the Isotopetable with G4IonTable. … … 156 155 NSourceBin = 1; 157 156 SBin[0] = 0.* s; 158 SBin[1] = 1 e10* s;157 SBin[1] = 1.* s; 159 158 SProfile[0] = 1.; 160 SProfile[1] = 1.;159 SProfile[1] = 0.; 161 160 NDecayBin = 1; 162 DBin[0] = 9.9e9* s ;163 DBin[1] = 1 e10* s;161 DBin[0] = 0. * s ; 162 DBin[1] = 1. * s; 164 163 DProfile[0] = 1.; 165 164 DProfile[1] = 0.; 165 decayWindows[0] = 0; 166 G4RadioactivityTable* rTable = new G4RadioactivityTable() ; 167 theRadioactivityTables.push_back(rTable); 166 168 NSplit = 1; 167 169 AnalogueMC = true ; … … 187 189 delete aPhysicsTable; 188 190 } 189 // delete theIsotopeTable;190 delete theRadioactiveDecaymessenger;191 // delete theIsotopeTable; 192 delete theRadioactiveDecaymessenger; 191 193 } 192 194 … … 197 199 // 198 200 G4bool G4RadioactiveDecay::IsApplicable( const G4ParticleDefinition & 199 201 aParticle) 200 202 { 201 203 // … … 224 226 // 225 227 G4bool G4RadioactiveDecay::IsLoaded(const G4ParticleDefinition & 226 aParticle)228 aParticle) 227 229 { 228 230 // … … 232 234 // 233 235 return std::binary_search(LoadedNuclei.begin(), 234 235 236 LoadedNuclei.end(), 237 aParticle.GetParticleName()); 236 238 } 237 239 //////////////////////////////////////////////////////////////////////////////// … … 242 244 void G4RadioactiveDecay::SelectAVolume(const G4String aVolume) 243 245 { 244 246 245 247 G4LogicalVolumeStore *theLogicalVolumes; 246 248 G4LogicalVolume *volume; … … 331 333 G4cout << " RDM removed from all volumes" << G4endl; 332 334 #endif 333 335 334 336 } 335 337 … … 340 342 // 341 343 G4bool G4RadioactiveDecay::IsRateTableReady(const G4ParticleDefinition & 342 aParticle)344 aParticle) 343 345 { 344 346 // … … 351 353 { 352 354 if (theDecayRateTableVector[i].GetIonName() == aParticleName) 353 355 return true ; 354 356 } 355 357 return false; … … 363 365 // 364 366 void G4RadioactiveDecay::GetDecayRateTable(const G4ParticleDefinition & 365 aParticle)367 aParticle) 366 368 { 367 369 … … 376 378 } 377 379 #ifdef G4VERBOSE 378 379 380 381 382 380 if (GetVerboseLevel()>0) 381 { 382 G4cout <<"The DecayRate Table for " 383 << aParticleName << " is selected." << G4endl; 384 } 383 385 #endif 384 386 } … … 388 390 // GetTaoTime 389 391 // 390 // to perform the conv ilution of the sourcetimeprofile function with the392 // to perform the convolution of the source time profile function with the 391 393 // decay constants in the decay chain. 392 394 // … … 408 410 } 409 411 taotime += SProfile[nbin] * (1-std::exp(-(t-SBin[nbin])/tao)); 412 if (taotime < 0.) { 413 G4cout <<" Tao time: " <<taotime << " reset to zero!"<<G4endl; 414 taotime = 0.; 415 } 416 410 417 #ifdef G4VERBOSE 411 418 if (GetVerboseLevel()>1) 412 419 {G4cout <<" Tao time: " <<taotime <<G4endl;} 413 420 #endif 414 return 415 } 416 421 return taotime ; 422 } 423 417 424 //////////////////////////////////////////////////////////////////////////////// 418 425 // … … 433 440 #ifdef G4VERBOSE 434 441 if (GetVerboseLevel()>1) 435 {G4cout <<" Decay time: " <<decaytime <<"[ns]" <<G4endl;}442 {G4cout <<" Decay time: " <<decaytime/s <<"[s]" <<G4endl;} 436 443 #endif 437 444 return decaytime; … … 447 454 G4int G4RadioactiveDecay::GetDecayTimeBin(const G4double aDecayTime) 448 455 { 449 for (G4int i = 0; i < NDecayBin; i++) 450 { 451 if ( aDecayTime > DBin[i]) return i+1; 452 } 453 return 1; 456 G4int i = 0; 457 while ( aDecayTime > DBin[i] ) i++; 458 return i; 454 459 } 455 460 //////////////////////////////////////////////////////////////////////////////// … … 472 477 G4ParticleDefinition* theParticleDef = theParticle->GetDefinition(); 473 478 G4double theLife = theParticleDef->GetPDGLifeTime(); 474 479 475 480 #ifdef G4VERBOSE 476 481 if (GetVerboseLevel()>2) … … 485 490 else if (theLife < 0.0) {meanlife = DBL_MAX;} 486 491 else {meanlife = theLife;} 487 // set the meanlife to zero for excited istopes which is not in the RDM database492 // set the meanlife to zero for excited istopes which is not in the RDM database 488 493 if (((const G4Ions*)(theParticleDef))->GetExcitationEnergy() > 0. && meanlife == DBL_MAX) {meanlife = 0.;} 489 494 } 490 495 #ifdef G4VERBOSE 491 496 if (GetVerboseLevel()>1) 492 {G4cout <<"mean life time: " <<meanlife <<"[ns]" <<G4endl;}493 #endif 494 497 {G4cout <<"mean life time: " <<meanlife/s <<"[s]" <<G4endl;} 498 #endif 499 495 500 return meanlife; 496 501 } … … 507 512 // constants 508 513 G4bool isOutRange ; 509 514 510 515 // get particle 511 516 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle(); 512 517 513 518 // returns the mean free path in GEANT4 internal units 514 519 G4double pathlength; … … 516 521 G4double aCtau = c_light * aParticleDef->GetPDGLifeTime(); 517 522 G4double aMass = aParticle->GetMass(); 518 523 519 524 #ifdef G4VERBOSE 520 525 if (GetVerboseLevel()>2) { … … 525 530 } 526 531 #endif 527 532 528 533 // check if the particle is stable? 529 534 if (aParticleDef->GetPDGStable()) { 530 535 pathlength = DBL_MAX; 531 536 532 537 } else if (aCtau < 0.0) { 533 538 pathlength = DBL_MAX; 534 539 535 540 //check if the particle has very short life time ? 536 541 } else if (aCtau < DBL_MIN) { 537 542 pathlength = DBL_MIN; 538 543 539 544 //check if zero mass 540 545 } else if (aMass < DBL_MIN) { … … 545 550 } 546 551 #endif 547 548 549 550 551 552 553 554 555 556 557 558 559 560 561 562 #ifdef G4VERBOSE 563 564 565 566 567 568 #endif 569 570 571 572 573 574 575 #ifdef G4VERBOSE 576 577 578 579 #endif 580 552 } else { 553 //calculate the mean free path 554 // by using normalized kinetic energy (= Ekin/mass) 555 G4double rKineticEnergy = aParticle->GetKineticEnergy()/aMass; 556 if ( rKineticEnergy > HighestBinValue) { 557 // beta >> 1 558 pathlength = ( rKineticEnergy + 1.0)* aCtau; 559 } else if ( rKineticEnergy > LowestBinValue) { 560 // check if aPhysicsTable exists 561 if (aPhysicsTable == 0) BuildPhysicsTable(*aParticleDef); 562 // beta is in the range valid for PhysicsTable 563 pathlength = aCtau * 564 ((*aPhysicsTable)(0))-> GetValue(rKineticEnergy,isOutRange); 565 } else if ( rKineticEnergy < DBL_MIN ) { 566 // too slow particle 567 #ifdef G4VERBOSE 568 if (GetVerboseLevel()>2) { 569 G4cout << "G4Decay::GetMeanFreePath() !!particle stops!!"; 570 G4cout << aParticleDef->GetParticleName() << G4endl; 571 G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]"; 572 } 573 #endif 574 pathlength = DBL_MIN; 575 } else { 576 // beta << 1 577 pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ; 578 } 579 } 580 #ifdef G4VERBOSE 581 if (GetVerboseLevel()>1) { 582 G4cout << "mean free path: "<< pathlength/m << "[m]" << G4endl; 583 } 584 #endif 585 return pathlength; 581 586 } 582 587 //////////////////////////////////////////////////////////////////////////////// … … 604 609 G4int i; 605 610 for ( i = 0 ; i < TotBin ; i++ ) { 606 607 608 611 gammainv = 1.0/(aVector->GetLowEdgeEnergy(i) + 1.0); 612 beta = std::sqrt((1.0 - gammainv)*(1.0 +gammainv)); 613 aVector->PutValue(i, beta/gammainv); 609 614 } 610 615 aPhysicsTable->insert(aVector); … … 636 641 G4cout << "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files." << G4endl; 637 642 throw G4HadronicException(__FILE__, __LINE__, 638 "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files.");643 "Please setenv G4RADIOACTIVEDATA to point to the radioactive decay data files."); 639 644 } 640 645 G4String dirName = getenv("G4RADIOACTIVEDATA"); … … 647 652 G4String file = os.str(); 648 653 649 654 650 655 std::ifstream DecaySchemeFile(file); 651 656 G4bool found(false); … … 664 669 modeSumBR[i] = 0.0; 665 670 } 666 671 667 672 G4bool complete(false); 668 673 char inputChars[80]={' '}; … … 718 723 a/=1000.; 719 724 c/=1000.; 720 725 721 726 // cout<< "The decay mode is [LoadTable] "<<theDecayMode<<G4endl; 722 727 723 728 switch (theDecayMode) { 724 729 case IT: … … 750 755 e0 = c*MeV/0.511; 751 756 n = aBetaFermiFunction->GetFFN(e0); 752 757 753 758 // now to work out the histogram and initialise the random generator 754 759 G4int npti = 100; … … 765 770 // Special treatment for K-40 (problem #1068) (flei,06/05/2010) 766 771 if (theParentNucleus.GetParticleName() == "K40[0.0]") f *= 767 (std::pow((g*g-1),3)+std::pow((ee-g),6)+7*(g*g-1)*std::pow((ee-g),2)*(g*g-1+std::pow((ee-g),2)));772 (std::pow((g*g-1),3)+std::pow((ee-g),6)+7*(g*g-1)*std::pow((ee-g),2)*(g*g-1+std::pow((ee-g),2))); 768 773 pdf[ptn] = f*aBetaFermiFunction->GetFF(e); 769 774 } … … 798 803 if (e0 > 0.) { 799 804 G4BetaFermiFunction* aBetaFermiFunction = new G4BetaFermiFunction (A, -(Z-1)); 800 805 801 806 n = aBetaFermiFunction->GetFFN(e0); 802 807 803 808 // now to work out the histogram and initialise the random generator 804 809 G4int npti = 100; … … 824 829 theDecayTable->Insert(aBetaPlusChannel); 825 830 modeSumBR[2] += b; 826 831 827 832 delete[] pdf; 828 833 delete aBetaFermiFunction; … … 909 914 G4Exception("G4RadioactiveDecay::LoadDecayTable()", "601", 910 915 FatalException, "Error in decay mode selection"); 911 916 912 917 } 913 918 } … … 972 977 // 973 978 void G4RadioactiveDecay::SetDecayRate(G4int theZ, G4int theA, G4double theE, 974 975 979 G4int theG, std::vector<G4double> theRates, 980 std::vector<G4double> theTaos) 976 981 { 977 982 //fill the decay rate vector … … 993 998 // 2) Add the coefficiencies to the decay rate table vector 994 999 // 995 1000 996 1001 // 997 1002 // Create and initialise variables used in the method. … … 999 1004 1000 1005 theDecayRateVector.clear(); 1001 1006 1002 1007 G4int nGeneration = 0; 1003 1008 std::vector<G4double> rates; 1004 1009 std::vector<G4double> taos; 1005 1010 1006 1011 // start rate is -1. 1012 // Eq.4.26 of the Technical Note 1007 1013 rates.push_back(-1.); 1008 1014 // … … 1012 1018 G4double E = ((const G4Ions*)(&theParentNucleus))->GetExcitationEnergy(); 1013 1019 G4double tao = theParentNucleus.GetPDGLifeTime(); 1020 if (tao < 0.) tao = 1e-30; 1014 1021 taos.push_back(tao); 1015 1022 G4int nEntry = 0; 1016 1023 1017 1024 //fill the decay rate with the intial isotope data 1018 1025 SetDecayRate(Z,A,E,nGeneration,rates,taos); … … 1034 1041 G4AlphaDecayChannel *theAlphaChannel = 0; 1035 1042 G4RadioactiveDecayMode theDecayMode; 1036 // G4NuclearLevelManager levelManager;1037 //const G4NuclearLevel* level;1038 1043 G4double theBR = 0.0; 1039 1044 G4int AP = 0; … … 1056 1061 // 1057 1062 theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable()); 1058 1063 1059 1064 while (!stable) { 1060 1065 nGeneration++; … … 1077 1082 aParentNucleus->SetDecayTable(LoadDecayTable(*aParentNucleus)); 1078 1083 } 1084 1085 G4DecayTable *theDecayTable = new G4DecayTable(); 1079 1086 aTempDecayTable = aParentNucleus->GetDecayTable(); 1087 for (i=0; i< 7; i++) brs[i] = 0.0; 1088 1080 1089 // 1081 1090 // 1082 1091 // Go through the decay table and to combine the same decay channels 1083 1092 // 1084 for (i=0; i< 7; i++) brs[i] = 0.0;1085 1086 G4DecayTable *theDecayTable = new G4DecayTable();1087 1088 1093 for (i=0; i<aTempDecayTable->entries(); i++) { 1089 1094 theChannel = aTempDecayTable->GetDecayChannel(i); … … 1099 1104 1100 1105 if (std::abs(daughterExcitation - level->Energy()) < levelTolerance) { 1101 1102 // Level hafe life is in ns and the gate as 1 micros by default 1103 if ( theDecayMode == 0 && level->HalfLife()*ns >= halflifethreshold ){ 1104 // some further though may needed here 1106 // Level half-life is in ns and the threshold is set to 1 micros by default, user can set it via the UI command 1107 if (level->HalfLife()*ns >= halflifethreshold ){ 1108 // save the metastable nucleus 1105 1109 theDecayTable->Insert(theChannel); 1106 1110 } … … 1120 1124 brs[3] = brs[4] =brs[5] = 0.0; 1121 1125 for (i= 0; i<7; i++){ 1122 if (brs[i] > 0 ) {1126 if (brs[i] > 0.) { 1123 1127 switch ( i ) { 1124 1128 case 0: … … 1127 1131 // Decay mode is isomeric transition. 1128 1132 // 1129 1133 1130 1134 theITChannel = new G4ITDecayChannel 1131 1135 (0, (const G4Ions*) aParentNucleus, brs[0]); 1132 1136 1133 1137 theDecayTable->Insert(theITChannel); 1134 1138 break; 1135 1139 1136 1140 case 1: 1137 1141 // … … 1142 1146 brs[1], 0.*MeV, 0.*MeV, 1, false, 0); 1143 1147 theDecayTable->Insert(theBetaMinusChannel); 1144 1148 1145 1149 break; 1146 1150 1147 1151 case 2: 1148 1152 // … … 1154 1158 theDecayTable->Insert(theBetaPlusChannel); 1155 1159 break; 1156 1160 1157 1161 case 6: 1158 1162 // … … 1164 1168 theDecayTable->Insert(theAlphaChannel); 1165 1169 break; 1166 1170 1167 1171 default: 1168 1172 break; … … 1170 1174 } 1171 1175 } 1172 1173 1176 // 1174 1177 // loop over all braches in theDecayTable … … 1179 1182 theBR = theChannel->GetBR(); 1180 1183 theDaughterNucleus = theNuclearDecayChannel->GetDaughterNucleus(); 1181 1182 // 1183 // now test if the daughterNucleus is a valid one 1184 // 1185 if (IsApplicable(*theDaughterNucleus) && theBR 1186 && aParentNucleus != theDaughterNucleus ) { 1184 // first check if the decay of the original nucleus is an IT channel, if true create a new groud-level nucleus 1185 if (theNuclearDecayChannel->GetDecayMode() == IT && nGeneration == 1 ) { 1186 A = ((const G4Ions*)(theDaughterNucleus))->GetAtomicMass(); 1187 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber(); 1188 theDaughterNucleus=theIonTable->GetIon(Z,A,0.); 1189 } 1190 if (IsApplicable(*theDaughterNucleus) && theBR && aParentNucleus != theDaughterNucleus) { 1187 1191 // need to make sure daugher has decaytable 1188 1192 if (!IsLoaded(*theDaughterNucleus)){ … … 1194 1198 Z = ((const G4Ions*)(theDaughterNucleus))->GetAtomicNumber(); 1195 1199 E = ((const G4Ions*)(theDaughterNucleus))->GetExcitationEnergy(); 1196 1200 1197 1201 TaoPlus = theDaughterNucleus->GetPDGLifeTime(); 1198 1202 // cout << TaoPlus <<G4endl; 1199 if (TaoPlus > 0.) { 1200 // first set the taos, one simply need to add to the parent ones 1201 taos.clear(); 1202 taos = TP; 1203 taos.push_back(TaoPlus); 1204 // now calculate the coefficiencies 1205 // 1206 // they are in two parts, first the les than n ones 1207 rates.clear(); 1208 size_t k; 1209 for (k = 0; k < RP.size(); k++){ 1203 if (TaoPlus <= 0.) TaoPlus = 1e-30; 1204 1205 // first set the taos, one simply need to add to the parent ones 1206 taos.clear(); 1207 taos = TP; 1208 taos.push_back(TaoPlus); 1209 // now calculate the coefficiencies 1210 // 1211 // they are in two parts, first the less than n ones 1212 // Eq 4.24 of the TN 1213 rates.clear(); 1214 size_t k; 1215 for (k = 0; k < RP.size(); k++){ 1216 if ((TP[k]-TaoPlus) == 0.) { 1217 theRate = 1e30; 1218 } else { 1210 1219 theRate = TP[k]/(TP[k]-TaoPlus) * theBR * RP[k]; 1211 rates.push_back(theRate);1212 1220 } 1213 // 1214 // the sencond part: the n:n coefficiency 1215 theRate = 0.; 1216 for (k = 0; k < RP.size(); k++){ 1217 theRate -=TaoPlus/(TP[k]-TaoPlus) * theBR * RP[k]; 1221 rates.push_back(theRate); 1222 } 1223 // 1224 // the sencond part: the n:n coefficiency 1225 // Eq 4.25 of the TN. Note Yn+1 is zero apart from Y1 which is -1 as treated at line 1013 1226 // 1227 theRate = 0.; 1228 G4double aRate; 1229 for (k = 0; k < RP.size(); k++){ 1230 if ((TP[k]-TaoPlus) == 0.) { 1231 aRate = 1e30; 1232 } else { 1233 aRate = TaoPlus/(TP[k]-TaoPlus) * theBR * RP[k]; 1218 1234 } 1219 rates.push_back(theRate); 1220 SetDecayRate (Z,A,E,nGeneration,rates,taos); 1221 theDecayRateVector.push_back(theDecayRate); 1222 nEntry++; 1223 } 1235 theRate -= aRate; 1236 } 1237 rates.push_back(theRate); 1238 SetDecayRate (Z,A,E,nGeneration,rates,taos); 1239 theDecayRateVector.push_back(theDecayRate); 1240 nEntry++; 1224 1241 } 1225 1242 } … … 1233 1250 if (nS == nT) stable = true; 1234 1251 } 1235 1252 1236 1253 //end of while loop 1237 1254 // the calculation completed here 1238 1239 1255 1256 1240 1257 // fill the first part of the decay rate table 1241 1258 // which is the name of the original particle (isotope) … … 1246 1263 // now fill the decay table with the newly completed decay rate vector 1247 1264 // 1248 1265 1249 1266 theDecayRateTable.SetItsRates(theDecayRateVector); 1250 1267 1251 1268 // 1252 1269 // finally add the decayratetable to the tablevector … … 1254 1271 theDecayRateTableVector.push_back(theDecayRateTable); 1255 1272 } 1256 1273 1257 1274 //////////////////////////////////////////////////////////////////////////////// 1258 1275 // … … 1267 1284 std::ifstream infile ( filename, std::ios::in ); 1268 1285 if ( !infile ) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException, "Unable to open source data file" ); 1269 1270 floatbin, flux;1286 1287 G4double bin, flux; 1271 1288 NSourceBin = -1; 1272 1289 while (infile >> bin >> flux ) { … … 1277 1294 } 1278 1295 SetAnalogueMonteCarlo(0); 1296 infile.close(); 1297 1279 1298 #ifdef G4VERBOSE 1280 1299 if (GetVerboseLevel()>1) … … 1294 1313 std::ifstream infile ( filename, std::ios::in); 1295 1314 if ( !infile ) G4Exception(__FILE__, G4inttostring(__LINE__), FatalException, "Unable to open bias data file" ); 1296 1297 float bin, flux; 1315 1316 G4double bin, flux; 1317 G4int dWindows = 0; 1318 G4int i ; 1319 1320 theRadioactivityTables.clear(); 1321 // for (i = 0; i<100; i++) decayWindows[i] = -1; 1322 1298 1323 NDecayBin = -1; 1299 1324 while (infile >> bin >> flux ) { … … 1302 1327 DBin[NDecayBin] = bin * s; 1303 1328 DProfile[NDecayBin] = flux; 1329 if (flux > 0.) { 1330 decayWindows[NDecayBin] = dWindows; 1331 dWindows++; 1332 G4RadioactivityTable *rTable = new G4RadioactivityTable() ; 1333 theRadioactivityTables.push_back(rTable); 1334 } 1304 1335 } 1305 G4int i ;1306 1336 for ( i = 1; i<= NDecayBin; i++) DProfile[i] += DProfile[i-1]; 1307 1337 for ( i = 0; i<= NDecayBin; i++) DProfile[i] /= DProfile[NDecayBin]; 1338 // converted to accumulated probabilities 1339 // 1308 1340 SetAnalogueMonteCarlo(0); 1341 infile.close(); 1342 1309 1343 #ifdef G4VERBOSE 1310 1344 if (GetVerboseLevel()>1) … … 1315 1349 //////////////////////////////////////////////////////////////////////////////// 1316 1350 // 1317 //1318 1351 // DecayIt 1319 1352 // 1320 G4VParticleChange* G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step& )1321 { 1322 // 1353 G4VParticleChange* 1354 G4RadioactiveDecay::DecayIt(const G4Track& theTrack, const G4Step&) 1355 { 1323 1356 // Initialize the G4ParticleChange object. Get the particle details and the 1324 1357 // decay table. 1325 // 1358 1326 1359 fParticleChangeForRadDecay.Initialize(theTrack); 1327 1360 const G4DynamicParticle* theParticle = theTrack.GetDynamicParticle(); … … 1329 1362 1330 1363 // First check whether RDM applies to the current logical volume 1331 // 1332 if(!std::binary_search(ValidVolumes.begin(), 1333 ValidVolumes.end(), 1334 theTrack.GetVolume()->GetLogicalVolume()->GetName())) 1335 { 1336 #ifdef G4VERBOSE 1337 if (GetVerboseLevel()>0) 1338 { 1339 G4cout <<"G4RadioactiveDecay::DecayIt : " 1340 << theTrack.GetVolume()->GetLogicalVolume()->GetName() 1341 << " is not selected for the RDM"<< G4endl; 1342 G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl; 1343 G4cout << " The Valid volumes are " << G4endl; 1344 for (size_t i = 0; i< ValidVolumes.size(); i++) 1345 G4cout << ValidVolumes[i] << G4endl; 1346 } 1347 #endif 1348 fParticleChangeForRadDecay.SetNumberOfSecondaries(0); 1349 // 1350 // 1351 // Kill the parent particle. 1352 // 1353 fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ; 1354 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0); 1355 ClearNumberOfInteractionLengthLeft(); 1356 return &fParticleChangeForRadDecay; 1364 1365 if (!std::binary_search(ValidVolumes.begin(), ValidVolumes.end(), 1366 theTrack.GetVolume()->GetLogicalVolume()->GetName())) { 1367 #ifdef G4VERBOSE 1368 if (GetVerboseLevel()>0) { 1369 G4cout <<"G4RadioactiveDecay::DecayIt : " 1370 << theTrack.GetVolume()->GetLogicalVolume()->GetName() 1371 << " is not selected for the RDM"<< G4endl; 1372 G4cout << " There are " << ValidVolumes.size() << " volumes" << G4endl; 1373 G4cout << " The Valid volumes are " << G4endl; 1374 for (size_t i = 0; i< ValidVolumes.size(); i++) G4cout << ValidVolumes[i] << G4endl; 1357 1375 } 1358 1376 #endif 1377 fParticleChangeForRadDecay.SetNumberOfSecondaries(0); 1378 1379 // Kill the parent particle. 1380 1381 fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ; 1382 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0); 1383 ClearNumberOfInteractionLengthLeft(); 1384 return &fParticleChangeForRadDecay; 1385 } 1386 1359 1387 // now check is the particle is valid for RDM 1360 // 1361 if (!(IsApplicable(*theParticleDef))) 1362 { 1363 // 1364 // The particle is not a Ion or outside the nucleuslimits for decay 1365 // 1366 #ifdef G4VERBOSE 1367 if (GetVerboseLevel()>0) 1368 { 1369 G4cerr <<"G4RadioactiveDecay::DecayIt : " 1370 <<theParticleDef->GetParticleName() 1371 << " is not a valid nucleus for the RDM"<< G4endl; 1372 } 1373 #endif 1374 fParticleChangeForRadDecay.SetNumberOfSecondaries(0); 1375 // 1376 // 1377 // Kill the parent particle. 1378 // 1379 fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ; 1380 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0); 1381 ClearNumberOfInteractionLengthLeft(); 1382 return &fParticleChangeForRadDecay; 1388 1389 if (!(IsApplicable(*theParticleDef))) { 1390 // 1391 // The particle is not a Ion or outside the nucleuslimits for decay 1392 // 1393 #ifdef G4VERBOSE 1394 if (GetVerboseLevel()>0) { 1395 G4cerr <<"G4RadioactiveDecay::DecayIt : " 1396 <<theParticleDef->GetParticleName() 1397 << " is not a valid nucleus for the RDM"<< G4endl; 1383 1398 } 1384 1399 #endif 1400 fParticleChangeForRadDecay.SetNumberOfSecondaries(0); 1401 1402 // 1403 // Kill the parent particle. 1404 // 1405 fParticleChangeForRadDecay.ProposeTrackStatus( fStopAndKill ) ; 1406 fParticleChangeForRadDecay.ProposeLocalEnergyDeposit(0.0); 1407 ClearNumberOfInteractionLengthLeft(); 1408 return &fParticleChangeForRadDecay; 1409 } 1410 1385 1411 if (!IsLoaded(*theParticleDef)) 1386 1412 { … … 1388 1414 } 1389 1415 G4DecayTable *theDecayTable = theParticleDef->GetDecayTable(); 1390 1416 1391 1417 if (theDecayTable == 0 || theDecayTable->entries() == 0 ) 1392 1418 { … … 1423 1449 G4ThreeVector currentPosition; 1424 1450 currentPosition = theTrack.GetPosition(); 1425 1451 1426 1452 // check whether use Analogue or VR implementation 1427 1453 // … … 1454 1480 G4double ParentEnergy = theParticle->GetTotalEnergy(); 1455 1481 G4ThreeVector ParentDirection(theParticle->GetMomentumDirection()); 1456 1482 1457 1483 if (theTrack.GetTrackStatus() == fStopButAlive ) 1458 1484 { … … 1501 1527 (products->PopProducts(), finalGlobalTime, currentPosition); 1502 1528 secondary->SetGoodForTrackingFlag(); 1503 1529 secondary->SetTouchableHandle(theTrack.GetTouchableHandle()); 1504 1530 fParticleChangeForRadDecay.AddSecondary(secondary); 1505 1531 } … … 1521 1547 if (!IsRateTableReady(*theParticleDef)) { 1522 1548 // if the decayrates are not ready, calculate them and 1523 1549 // add to the rate table vector 1524 1550 AddDecayRateTable(*theParticleDef); 1525 1551 } … … 1538 1564 G4double taotime; 1539 1565 G4double decayRate; 1540 1566 1541 1567 size_t i; 1542 1568 size_t j; … … 1557 1583 for (G4int n = 0; n < NSplit; n++) 1558 1584 { 1559 /*1560 //1561 1585 // Get the decay time following the decay probability function 1562 1586 // suppllied by user 1563 //1587 1564 1588 G4double theDecayTime = GetDecayTime(); 1565 1566 1589 G4int nbin = GetDecayTimeBin(theDecayTime); 1567 1590 1568 // c laculate the first part of the weight function1591 // calculate the first part of the weight function 1569 1592 1570 G4double weight1 =1./DProfile[nbin-1] 1571 *(DBin[nbin]-DBin[nbin-1]) 1572 /NSplit; 1573 if (nbin > 1) { 1574 weight1 = 1./(DProfile[nbin]-DProfile[nbin-2]) 1575 *(DBin[nbin]-DBin[nbin-1]) 1576 /NSplit;} 1593 G4double weight1 = 1.; 1594 if (nbin == 1) { 1595 weight1 = 1./DProfile[nbin-1] 1596 *(DBin[nbin]-DBin[nbin-1])/NSplit; 1597 } else if (nbin > 1) { 1598 weight1 = 1./(DProfile[nbin]-DProfile[nbin-2]) 1599 *(DBin[nbin]-DBin[nbin-1])/NSplit; 1600 } 1601 1577 1602 // it should be calculated in seconds 1578 1603 weight1 /= s ; 1579 */ 1580 // 1604 1581 1605 // loop over all the possible secondaries of the nucleus 1582 1606 // the first one is itself. 1583 // 1584 for ( 1607 1608 for (i = 0; i<theDecayRateVector.size(); i++){ 1585 1609 PZ = theDecayRateVector[i].GetZ(); 1586 1610 PA = theDecayRateVector[i].GetA(); … … 1589 1613 PR = theDecayRateVector[i].GetDecayRateC(); 1590 1614 1591 //1592 // Get the decay time following the decay probability function1593 // suppllied by user1594 //1595 G4double theDecayTime = GetDecayTime();1596 1597 G4int nbin = GetDecayTimeBin(theDecayTime);1598 1599 // claculate the first part of the weight function1600 1601 G4double weight1 =1./DProfile[nbin-1]1602 *(DBin[nbin]-DBin[nbin-1])1603 /NSplit;1604 if (nbin > 1) {1605 weight1 = 1./(DProfile[nbin]-DProfile[nbin-2])1606 *(DBin[nbin]-DBin[nbin-1])1607 /NSplit;}1608 // it should be calculated in seconds1609 weight1 /= s ;1610 1611 1615 // a temprary products buffer and its contents is transfered to 1612 1616 // the products at the end of the loop 1613 // 1617 1614 1618 G4DecayProducts *tempprods = 0; 1615 1619 1616 1620 // calculate the decay rate of the isotope 1617 // one need to fold the the source bias function with the decaytime 1618 // 1621 // decayRate is the radioactivity of isotope (PZ,PA,PE) at the 1622 // time 'theDecayTime' 1623 // it will be used to calculate the statistical weight of the 1624 // decay products of this isotope 1625 1626 // G4cout <<"PA= "<< PA << " PZ= " << PZ << " PE= "<< PE <<G4endl; 1619 1627 decayRate = 0.; 1620 1628 for ( j = 0; j < PT.size(); j++){ 1621 1629 taotime = GetTaoTime(theDecayTime,PT[j]); 1622 1630 decayRate -= PR[j] * taotime; 1631 // Eq.4.23 of of the TN 1632 // note the negative here is required as the rate in the eqation is defined to be negative, 1633 // i.e. decay away, but we need pasitive value here. 1634 1635 // G4cout << j << "\t"<< PT[j]/s <<"\t"<<PR[j]<< "\t" << decayRate << G4endl ; 1623 1636 } 1624 1625 // decayRatehe radioactivity of isotope (PZ,PA,PE) at the 1626 // time 'theDecayTime' 1627 // it will be used to calculate the statistical weight of the 1628 // decay products of this isotope 1629 1630 1631 // 1637 1632 1638 // now calculate the statistical weight 1633 // 1634 1635 G4double weight = weight1*decayRate; 1639 1640 // one need to fold the the source bias function with the decaytime 1641 // also need to include the track weight! (F.Lei, 28/10/10) 1642 G4double weight = weight1*decayRate*theTrack.GetWeight(); 1643 1644 // add the isotope to the radioactivity tables 1645 // G4cout <<theDecayTime/s <<"\t"<<nbin<<G4endl; 1646 //G4cout << theTrack.GetWeight() <<"\t"<<weight1<<"\t"<<decayRate<< G4endl; 1647 theRadioactivityTables[decayWindows[nbin-1]]->AddIsotope(PZ,PA,PE,weight); 1648 1636 1649 // decay the isotope 1637 1650 theIonTable = (G4IonTable *)(G4ParticleTable::GetParticleTable()->GetIonTable()); 1638 1651 parentNucleus = theIonTable->GetIon(PZ,PA,PE); 1639 1652 1640 1653 // decide whther to apply branching ratio bias or not 1641 1654 // … … 1644 1657 ndecaych = G4int(theDecayTable->entries()*G4UniformRand()); 1645 1658 G4VDecayChannel *theDecayChannel = theDecayTable->GetDecayChannel(ndecaych); 1646 if (theDecayChannel == 0) 1647 { 1648 // Decay channel not found. 1649 #ifdef G4VERBOSE 1650 if (GetVerboseLevel()>0) 1651 { 1652 G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel"; 1653 G4cerr <<G4endl; 1654 theDecayTable ->DumpInfo(); 1655 } 1656 #endif 1659 if (theDecayChannel == 0) { 1660 // Decay channel not found. 1661 #ifdef G4VERBOSE 1662 if (GetVerboseLevel()>0) { 1663 G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel"; 1664 G4cerr <<G4endl; 1665 theDecayTable ->DumpInfo(); 1657 1666 } 1658 else 1659 { 1660 // A decay channel has been identified, so execute the DecayIt. 1661 G4double tempmass = parentNucleus->GetPDGMass(); 1662 tempprods = theDecayChannel->DecayIt(tempmass); 1663 weight *= (theDecayChannel->GetBR())*(theDecayTable->entries()); 1664 } 1665 } 1666 else { 1667 #endif 1668 } else { 1669 // A decay channel has been identified, so execute the DecayIt. 1670 G4double tempmass = parentNucleus->GetPDGMass(); 1671 tempprods = theDecayChannel->DecayIt(tempmass); 1672 weight *= (theDecayChannel->GetBR())*(theDecayTable->entries()); 1673 } 1674 } else { 1667 1675 tempprods = DoDecay(*parentNucleus); 1668 1676 } 1669 // 1677 1670 1678 // save the secondaries for buffers 1671 // 1679 1672 1680 numberOfSecondaries = tempprods->entries(); 1673 1681 currentTime = finalGlobalTime + theDecayTime; 1674 for (index=0; index < numberOfSecondaries; index++) 1675 { 1676 asecondaryparticle = tempprods->PopProducts(); 1677 if (asecondaryparticle->GetDefinition()->GetBaryonNumber() < 5){ 1678 pw.push_back(weight); 1679 ptime.push_back(currentTime); 1680 secondaryparticles.push_back(asecondaryparticle); 1681 } 1682 for (index=0; index < numberOfSecondaries; index++) { 1683 asecondaryparticle = tempprods->PopProducts(); 1684 if (asecondaryparticle->GetDefinition()->GetBaryonNumber() < 5){ 1685 pw.push_back(weight); 1686 ptime.push_back(currentTime); 1687 secondaryparticles.push_back(asecondaryparticle); 1682 1688 } 1683 // 1689 } 1690 1684 1691 delete tempprods; 1685 1692 1686 1693 //end of i loop 1687 1694 } 1688 1695 1689 1696 // end of n loop 1690 } 1697 } 1698 1691 1699 // now deal with the secondaries in the two stl containers 1692 1700 // and submmit them back to the tracking manager … … 1699 1707 secondaryparticles[index], ptime[index], currentPosition); 1700 1708 secondary->SetGoodForTrackingFlag(); 1701 1709 secondary->SetTouchableHandle(theTrack.GetTouchableHandle()); 1702 1710 secondary->SetWeight(pw[index]); 1703 1711 fParticleChangeForRadDecay.AddSecondary(secondary); 1704 1712 } 1705 1713 // … … 1708 1716 //theTrack.SetTrackStatus(fStopButAlive); 1709 1717 //energyDeposit += theParticle->GetKineticEnergy(); 1710 1718 1711 1719 } 1712 1720 1713 1721 // 1714 1722 // Kill the parent particle. … … 1722 1730 // 1723 1731 ClearNumberOfInteractionLengthLeft(); 1724 1732 1725 1733 return &fParticleChangeForRadDecay ; 1726 1734 } 1727 1735 } 1728 1736 1729 //////////////////////////////////////////////////////////////////////////////// 1730 // 1737 /////////////////////////////////////////////////////////////////// 1731 1738 // 1732 1739 // DoDecay 1733 1740 // 1734 G4DecayProducts* G4RadioactiveDecay::DoDecay( G4ParticleDefinition& theParticleDef )1735 { 1736 G4DecayProducts *products = 0; 1737 //1738 // 1741 G4DecayProducts* 1742 G4RadioactiveDecay::DoDecay( G4ParticleDefinition& theParticleDef ) 1743 { 1744 G4DecayProducts* products = 0; 1745 1739 1746 // follow the decaytable and generate the secondaries... 1740 // 1741 // 1742 #ifdef G4VERBOSE 1743 if (GetVerboseLevel()>0) 1744 { 1745 G4cout<<"Begin of DoDecay..."<<G4endl; 1747 1748 #ifdef G4VERBOSE 1749 if (GetVerboseLevel()>0) G4cout<<"Begin of DoDecay..."<<G4endl; 1750 #endif 1751 1752 G4DecayTable* theDecayTable = theParticleDef.GetDecayTable(); 1753 1754 // Choose a decay channel. 1755 1756 #ifdef G4VERBOSE 1757 if (GetVerboseLevel()>0) G4cout <<"Selecte a channel..."<<G4endl; 1758 #endif 1759 1760 G4VDecayChannel* theDecayChannel = theDecayTable->SelectADecayChannel(); 1761 if (theDecayChannel == 0) { 1762 // Decay channel not found. 1763 1764 G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel"; 1765 G4cerr <<G4endl; 1766 theDecayTable ->DumpInfo(); 1767 } else { 1768 1769 // A decay channel has been identified, so execute the DecayIt. 1770 1771 #ifdef G4VERBOSE 1772 if (GetVerboseLevel()>1) { 1773 G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel addr:"; 1774 G4cerr <<theDecayChannel <<G4endl; 1746 1775 } 1747 1776 #endif 1748 G4DecayTable *theDecayTable = theParticleDef.GetDecayTable(); 1749 // 1750 // Choose a decay channel. 1751 // 1752 #ifdef G4VERBOSE 1753 if (GetVerboseLevel()>0) 1754 { 1755 G4cout <<"Selecte a channel..."<<G4endl; 1756 } 1757 #endif 1758 G4VDecayChannel *theDecayChannel = theDecayTable->SelectADecayChannel(); 1759 if (theDecayChannel == 0) 1760 { 1761 // Decay channel not found. 1762 // 1763 G4cerr <<"G4RadioactiveDecay::DoIt : can not determine decay channel"; 1764 G4cerr <<G4endl; 1765 theDecayTable ->DumpInfo(); 1766 } 1767 else 1768 { 1769 // 1770 // A decay channel has been identified, so execute the DecayIt. 1771 // 1772 #ifdef G4VERBOSE 1773 if (GetVerboseLevel()>1) 1774 { 1775 G4cerr <<"G4RadioactiveDecay::DoIt : selected decay channel addr:"; 1776 G4cerr <<theDecayChannel <<G4endl; 1777 } 1778 #endif 1779 1780 G4double tempmass = theParticleDef.GetPDGMass(); 1781 // 1782 1783 products = theDecayChannel->DecayIt(tempmass); 1784 1785 } 1777 1778 G4double tempmass = theParticleDef.GetPDGMass(); 1779 products = theDecayChannel->DecayIt(tempmass); 1780 } 1781 1786 1782 return products; 1787 1788 } 1789 1790 1791 1792 1793 1794 1795 1796 1797 1783 } -
trunk/source/processes/hadronic/models/radioactive_decay/src/G4RadioactiveDecaymessenger.cc
r1340 r1347 123 123 icmCmd->SetParameterName("applyICM",true); 124 124 icmCmd->SetDefaultValue(true); 125 icmCmd->AvailableForStates(G4State_PreInit);125 //icmCmd->AvailableForStates(G4State_PreInit); 126 126 // 127 127 // Command contols whether ARM will be applied or not … … 131 131 armCmd->SetParameterName("applyARM",true); 132 132 armCmd->SetDefaultValue(true); 133 armCmd->AvailableForStates(G4State_PreInit);133 //armCmd->AvailableForStates(G4State_PreInit); 134 134 // 135 135 // Command to set the h-l thresold for isomer production … … 140 140 // hlthCmd->SetRange("hlThreshold>0."); 141 141 hlthCmd->SetUnitCategory("Time"); 142 hlthCmd->AvailableForStates(G4State_PreInit);142 // hlthCmd->AvailableForStates(G4State_PreInit); 143 143 // 144 144 // Command to define the incident particle source time profile. … … 209 209 SetNucleusLimits(nucleuslimitsCmd->GetNewNucleusLimitsValue(newValues));} 210 210 else if (command==analoguemcCmd) { 211 G4int vl; 212 const char* t = newValues; 213 std::istringstream is(t); 214 is >> vl; 215 theRadioactiveDecayContainer->SetAnalogueMonteCarlo(vl!=0);} 211 theRadioactiveDecayContainer->SetAnalogueMonteCarlo(analoguemcCmd->GetNewBoolValue(newValues));} 216 212 else if (command==fbetaCmd) { 217 G4int vl; 218 const char* t = newValues; 219 std::istringstream is(t); 220 is >> vl; 221 theRadioactiveDecayContainer->SetFBeta(vl!=0);} 213 theRadioactiveDecayContainer->SetFBeta(fbetaCmd->GetNewBoolValue(newValues));} 222 214 else if (command==avolumeCmd) {theRadioactiveDecayContainer-> 223 215 SelectAVolume(newValues);} … … 229 221 DeselectAllVolumes();} 230 222 else if (command==brbiasCmd) { 231 G4int vl; 232 const char* t = newValues; 233 std::istringstream is(t); 234 is >> vl; 235 theRadioactiveDecayContainer->SetBRBias(vl!=0);} 223 theRadioactiveDecayContainer->SetBRBias(brbiasCmd->GetNewBoolValue(newValues));} 236 224 else if (command==sourcetimeprofileCmd) {theRadioactiveDecayContainer-> 237 225 SetSourceTimeProfile(newValues);}
Note: See TracChangeset
for help on using the changeset viewer.