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 #include "G4EmBiasingManager.hh"
00050 #include "G4SystemOfUnits.hh"
00051 #include "G4MaterialCutsCouple.hh"
00052 #include "G4ProductionCutsTable.hh"
00053 #include "G4ProductionCuts.hh"
00054 #include "G4Region.hh"
00055 #include "G4RegionStore.hh"
00056 #include "G4Track.hh"
00057 #include "G4Electron.hh"
00058 #include "G4VEmModel.hh"
00059 #include "G4LossTableManager.hh"
00060 #include "G4ParticleChangeForLoss.hh"
00061 #include "G4ParticleChangeForGamma.hh"
00062
00063
00064
00065 G4EmBiasingManager::G4EmBiasingManager()
00066 : nForcedRegions(0),nSecBiasedRegions(0),eIonisation(0),
00067 currentStepLimit(0.0),startTracking(true)
00068 {
00069 fSafetyMin = 1.e-6*mm;
00070 theElectron = G4Electron::Electron();
00071 }
00072
00073
00074
00075 G4EmBiasingManager::~G4EmBiasingManager()
00076 {}
00077
00078
00079
00080 void G4EmBiasingManager::Initialise(const G4ParticleDefinition& part,
00081 const G4String& procName, G4int verbose)
00082 {
00083
00084
00085
00086 const G4ProductionCutsTable* theCoupleTable=
00087 G4ProductionCutsTable::GetProductionCutsTable();
00088 size_t numOfCouples = theCoupleTable->GetTableSize();
00089
00090 if(0 < nForcedRegions) { idxForcedCouple.resize(numOfCouples, -1); }
00091 if(0 < nSecBiasedRegions) { idxSecBiasedCouple.resize(numOfCouples, -1); }
00092
00093
00094 for (size_t j=0; j<numOfCouples; ++j) {
00095 const G4MaterialCutsCouple* couple =
00096 theCoupleTable->GetMaterialCutsCouple(j);
00097 const G4ProductionCuts* pcuts = couple->GetProductionCuts();
00098 if(0 < nForcedRegions) {
00099 for(G4int i=0; i<nForcedRegions; ++i) {
00100 if(forcedRegions[i]) {
00101 if(pcuts == forcedRegions[i]->GetProductionCuts()) {
00102 idxForcedCouple[j] = i;
00103 break;
00104 }
00105 }
00106 }
00107 }
00108 if(0 < nSecBiasedRegions) {
00109 for(G4int i=0; i<nSecBiasedRegions; ++i) {
00110 if(secBiasedRegions[i]) {
00111 if(pcuts == secBiasedRegions[i]->GetProductionCuts()) {
00112 idxSecBiasedCouple[j] = i;
00113 break;
00114 }
00115 }
00116 }
00117 }
00118 }
00119 if (nForcedRegions > 0 && 0 < verbose) {
00120 G4cout << " Forced Interaction is activated for "
00121 << part.GetParticleName() << " and "
00122 << procName
00123 << " inside G4Regions: " << G4endl;
00124 for (G4int i=0; i<nForcedRegions; ++i) {
00125 const G4Region* r = forcedRegions[i];
00126 if(r) { G4cout << " " << r->GetName() << G4endl; }
00127 }
00128 }
00129 if (nSecBiasedRegions > 0 && 0 < verbose) {
00130 G4cout << " Secondary biasing is activated for "
00131 << part.GetParticleName() << " and "
00132 << procName
00133 << " inside G4Regions: " << G4endl;
00134 for (G4int i=0; i<nSecBiasedRegions; ++i) {
00135 const G4Region* r = secBiasedRegions[i];
00136 if(r) {
00137 G4cout << " " << r->GetName()
00138 << " BiasingWeight= " << secBiasedWeight[i] << G4endl;
00139 }
00140 }
00141 }
00142 }
00143
00144
00145
00146 void G4EmBiasingManager::ActivateForcedInteraction(G4double val,
00147 const G4String& rname)
00148 {
00149 G4RegionStore* regionStore = G4RegionStore::GetInstance();
00150 G4String name = rname;
00151 if(name == "" || name == "world" || name == "World") {
00152 name = "DefaultRegionForTheWorld";
00153 }
00154 const G4Region* reg = regionStore->GetRegion(name, false);
00155 if(!reg) {
00156 G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
00157 << " G4Region <"
00158 << rname << "> is unknown" << G4endl;
00159 return;
00160 }
00161
00162
00163 if (0 < nForcedRegions) {
00164 for (G4int i=0; i<nForcedRegions; ++i) {
00165 if (reg == forcedRegions[i]) {
00166 lengthForRegion[i] = val;
00167 return;
00168 }
00169 }
00170 }
00171 if(val < 0.0) {
00172 G4cout << "### G4EmBiasingManager::ForcedInteraction WARNING: "
00173 << val << " < 0.0, so no activation for the G4Region <"
00174 << rname << ">" << G4endl;
00175 return;
00176 }
00177
00178
00179 forcedRegions.push_back(reg);
00180 lengthForRegion.push_back(val);
00181 ++nForcedRegions;
00182 }
00183
00184
00185
00186 void
00187 G4EmBiasingManager::ActivateSecondaryBiasing(const G4String& rname,
00188 G4double factor,
00189 G4double energyLimit)
00190 {
00191
00192
00193
00194 G4RegionStore* regionStore = G4RegionStore::GetInstance();
00195 G4String name = rname;
00196 if(name == "" || name == "world" || name == "World") {
00197 name = "DefaultRegionForTheWorld";
00198 }
00199 const G4Region* reg = regionStore->GetRegion(name, false);
00200 if(!reg) {
00201 G4cout << "### G4EmBiasingManager::ActivateBremsstrahlungSplitting WARNING: "
00202 << " G4Region <"
00203 << rname << "> is unknown" << G4endl;
00204 return;
00205 }
00206
00207
00208 G4int nsplit = 0;
00209 G4double w = factor;
00210
00211
00212 if(factor >= 1.0) {
00213 nsplit = G4lrint(factor);
00214 w = 1.0/G4double(nsplit);
00215
00216
00217 } else if(0.0 < factor) {
00218 nsplit = 1;
00219 w = 1.0/factor;
00220 }
00221
00222
00223 if (0 < nSecBiasedRegions) {
00224 for (G4int i=0; i<nSecBiasedRegions; ++i) {
00225 if (reg == secBiasedRegions[i]) {
00226 secBiasedWeight[i] = w;
00227 nBremSplitting[i] = nsplit;
00228 secBiasedEnegryLimit[i] = energyLimit;
00229 return;
00230 }
00231 }
00232 }
00233
00234
00235
00236
00237
00238
00239
00240 secBiasedRegions.push_back(reg);
00241 secBiasedWeight.push_back(w);
00242 nBremSplitting.push_back(nsplit);
00243 secBiasedEnegryLimit.push_back(energyLimit);
00244 ++nSecBiasedRegions;
00245
00246 }
00247
00248
00249
00250 G4double G4EmBiasingManager::GetStepLimit(G4int coupleIdx,
00251 G4double previousStep)
00252 {
00253 if(startTracking) {
00254 startTracking = false;
00255 G4int i = idxForcedCouple[coupleIdx];
00256 if(i < 0) {
00257 currentStepLimit = DBL_MAX;
00258 } else {
00259 currentStepLimit = lengthForRegion[i];
00260 if(currentStepLimit > 0.0) { currentStepLimit *= G4UniformRand(); }
00261 }
00262 } else {
00263 currentStepLimit -= previousStep;
00264 }
00265 if(currentStepLimit < 0.0) { currentStepLimit = 0.0; }
00266 return currentStepLimit;
00267 }
00268
00269
00270
00271 G4double
00272 G4EmBiasingManager::ApplySecondaryBiasing(
00273 std::vector<G4DynamicParticle*>& vd,
00274 const G4Track& track,
00275 G4VEmModel* currentModel,
00276 G4ParticleChangeForLoss* pPartChange,
00277 G4double& eloss,
00278 G4int coupleIdx,
00279 G4double tcut,
00280 G4double safety)
00281 {
00282 G4int index = idxSecBiasedCouple[coupleIdx];
00283 G4double weight = 1.0;
00284 if(0 <= index) {
00285 size_t n = vd.size();
00286
00287
00288
00289
00290 if(0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {
00291
00292 G4int nsplit = nBremSplitting[index];
00293
00294
00295 if(0 == nsplit) {
00296 if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }
00297
00298
00299 } if(1 == nsplit) {
00300 weight = ApplyRussianRoulette(vd, index);
00301
00302
00303 } else {
00304 G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
00305 G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();
00306
00307 weight = ApplySplitting(vd, track, currentModel, index, tcut);
00308
00309 pPartChange->SetProposedKineticEnergy(tmpEnergy);
00310 pPartChange->ProposeMomentumDirection(tmpMomDir);
00311 }
00312 }
00313 }
00314 return weight;
00315 }
00316
00317
00318
00319 G4double
00320 G4EmBiasingManager::ApplySecondaryBiasing(
00321 std::vector<G4DynamicParticle*>& vd,
00322 const G4Track& track,
00323 G4VEmModel* currentModel,
00324 G4ParticleChangeForGamma* pPartChange,
00325 G4double& eloss,
00326 G4int coupleIdx,
00327 G4double tcut,
00328 G4double safety)
00329 {
00330 G4int index = idxSecBiasedCouple[coupleIdx];
00331 G4double weight = 1.0;
00332 if(0 <= index) {
00333 size_t n = vd.size();
00334
00335
00336
00337
00338 if(0 < n && vd[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {
00339
00340 G4int nsplit = nBremSplitting[index];
00341
00342
00343 if(0 == nsplit) {
00344 if(safety > fSafetyMin) { ApplyRangeCut(vd, track, eloss, safety); }
00345
00346
00347 } if(1 == nsplit) {
00348 weight = ApplyRussianRoulette(vd, index);
00349
00350
00351 } else {
00352 G4double tmpEnergy = pPartChange->GetProposedKineticEnergy();
00353 G4ThreeVector tmpMomDir = pPartChange->GetProposedMomentumDirection();
00354
00355 weight = ApplySplitting(vd, track, currentModel, index, tcut);
00356
00357 pPartChange->SetProposedKineticEnergy(tmpEnergy);
00358 pPartChange->ProposeMomentumDirection(tmpMomDir);
00359 }
00360 }
00361 }
00362 return weight;
00363 }
00364
00365
00366
00367 G4double
00368 G4EmBiasingManager::ApplySecondaryBiasing(std::vector<G4Track*>& track,
00369 G4int coupleIdx)
00370 {
00371 G4int index = idxSecBiasedCouple[coupleIdx];
00372 G4double weight = 1.0;
00373 if(0 <= index) {
00374 size_t n = track.size();
00375
00376
00377
00378
00379 if(0 < n && track[0]->GetKineticEnergy() < secBiasedEnegryLimit[index]) {
00380
00381 G4int nsplit = nBremSplitting[index];
00382
00383
00384 if(1 == nsplit) {
00385 weight = secBiasedWeight[index];
00386 for(size_t k=0; k<n; ++k) {
00387 if(G4UniformRand()*weight > 1.0) {
00388 const G4Track* t = track[k];
00389 delete t;
00390 track[k] = 0;
00391 }
00392 }
00393 }
00394 }
00395 }
00396 return weight;
00397 }
00398
00399
00400
00401 void
00402 G4EmBiasingManager::ApplyRangeCut(std::vector<G4DynamicParticle*>& vd,
00403 const G4Track& track,
00404 G4double& eloss, G4double safety)
00405 {
00406 size_t n = vd.size();
00407 if(!eIonisation) {
00408 eIonisation = G4LossTableManager::Instance()->GetEnergyLossProcess(theElectron);
00409 }
00410 if(eIonisation) {
00411 for(size_t k=0; k<n; ++k) {
00412 const G4DynamicParticle* dp = vd[k];
00413 if(dp->GetDefinition() == theElectron) {
00414 G4double e = dp->GetKineticEnergy();
00415 if(eIonisation->GetRangeForLoss(e, track.GetMaterialCutsCouple()) < safety) {
00416 eloss += e;
00417 delete dp;
00418 vd[k] = 0;
00419 }
00420 }
00421 }
00422 }
00423 }
00424
00425
00426
00427 G4double
00428 G4EmBiasingManager::ApplySplitting(std::vector<G4DynamicParticle*>& vd,
00429 const G4Track& track,
00430 G4VEmModel* currentModel,
00431 G4int index,
00432 G4double tcut)
00433 {
00434
00435
00436 G4double weight = 1.0;
00437 size_t n = vd.size();
00438 G4double w = secBiasedWeight[index];
00439
00440 if(1 != n || 1.0 <= w) { return weight; }
00441
00442 G4double trackWeight = track.GetWeight();
00443 const G4DynamicParticle* dynParticle = track.GetDynamicParticle();
00444
00445 G4int nsplit = nBremSplitting[index];
00446
00447
00448 if(1 < nsplit && trackWeight>w) {
00449
00450 weight = w;
00451
00452 if(nsplit > (G4int)tmpSecondaries.size()) {
00453 tmpSecondaries.reserve(nsplit);
00454 }
00455 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
00456 for(G4int k=1; k<nsplit; ++k) {
00457 tmpSecondaries.clear();
00458 currentModel->SampleSecondaries(&tmpSecondaries, couple, dynParticle, tcut);
00459 for (size_t kk=0; kk<tmpSecondaries.size(); ++kk) {
00460 vd.push_back(tmpSecondaries[kk]);
00461 }
00462 }
00463 }
00464 return weight;
00465 }
00466
00467