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
00050
00051
00052
00053
00054
00055
00056
00057
00058
00059
00060
00061
00062
00063
00064
00065
00066
00067
00068
00069 #include "G4VMultipleScattering.hh"
00070 #include "G4PhysicalConstants.hh"
00071 #include "G4SystemOfUnits.hh"
00072 #include "G4LossTableManager.hh"
00073 #include "G4MaterialCutsCouple.hh"
00074 #include "G4Step.hh"
00075 #include "G4ParticleDefinition.hh"
00076 #include "G4VEmFluctuationModel.hh"
00077 #include "G4UnitsTable.hh"
00078 #include "G4ProductionCutsTable.hh"
00079 #include "G4Electron.hh"
00080 #include "G4GenericIon.hh"
00081 #include "G4TransportationManager.hh"
00082 #include "G4SafetyHelper.hh"
00083
00084
00085
00086 G4VMultipleScattering::G4VMultipleScattering(const G4String& name, G4ProcessType):
00087 G4VContinuousDiscreteProcess("msc", fElectromagnetic),
00088 numberOfModels(0),
00089 firstParticle(0),
00090 currParticle(0),
00091 stepLimit(fUseSafety),
00092 skin(1.0),
00093 facrange(0.04),
00094 facgeom(2.5),
00095 latDisplasment(true),
00096 isIon(false)
00097 {
00098 SetVerboseLevel(1);
00099 SetProcessSubType(fMultipleScattering);
00100 if("ionmsc" == name) { firstParticle = G4GenericIon::GenericIon(); }
00101
00102 geomMin = 1.e-6*CLHEP::mm;
00103 lowestKinEnergy = 1*eV;
00104
00105
00106 polarAngleLimit = 0.0;
00107
00108 physStepLimit = gPathLength = tPathLength = 0.0;
00109 fIonisation = 0;
00110
00111 pParticleChange = &fParticleChange;
00112 safetyHelper = 0;
00113 fPositionChanged = false;
00114 isActive = false;
00115
00116 modelManager = new G4EmModelManager();
00117 emManager = G4LossTableManager::Instance();
00118 emManager->Register(this);
00119
00120 warn = 0;
00121 }
00122
00123
00124
00125 G4VMultipleScattering::~G4VMultipleScattering()
00126 {
00127 if(1 < verboseLevel) {
00128 G4cout << "G4VMultipleScattering destruct " << GetProcessName()
00129 << G4endl;
00130 }
00131 delete modelManager;
00132 emManager->DeRegister(this);
00133 }
00134
00135
00136
00137 void G4VMultipleScattering::AddEmModel(G4int order, G4VEmModel* p,
00138 const G4Region* region)
00139 {
00140 G4VEmFluctuationModel* fm = 0;
00141 modelManager->AddEmModel(order, p, fm, region);
00142 if(p) { p->SetParticleChange(pParticleChange); }
00143 }
00144
00145
00146
00147 void G4VMultipleScattering::SetModel(G4VMscModel* p, G4int index)
00148 {
00149 ++warn;
00150 if(warn < 10) {
00151 G4cout << "### G4VMultipleScattering::SetModel is obsolete method "
00152 << "and will be removed for the next release." << G4endl;
00153 G4cout << " Please, use SetEmModel" << G4endl;
00154 }
00155 G4int n = mscModels.size();
00156 if(index >= n) { for(G4int i=n; i<=index; ++i) {mscModels.push_back(0);} }
00157 mscModels[index] = p;
00158 }
00159
00160
00161
00162 G4VMscModel* G4VMultipleScattering::Model(G4int index)
00163 {
00164 ++warn;
00165 if(warn < 10) {
00166 G4cout << "### G4VMultipleScattering::Model is obsolete method "
00167 << "and will be removed for the next release." << G4endl;
00168 G4cout << " Please, use EmModel" << G4endl;
00169 }
00170 G4VMscModel* p = 0;
00171 if(index >= 0 && index < G4int(mscModels.size())) { p = mscModels[index]; }
00172 return p;
00173 }
00174
00175
00176
00177 void G4VMultipleScattering::SetEmModel(G4VMscModel* p, G4int index)
00178 {
00179 G4int n = mscModels.size();
00180 if(index >= n) { for(G4int i=n; i<=index; ++i) { mscModels.push_back(0); } }
00181 mscModels[index] = p;
00182 }
00183
00184
00185
00186 G4VMscModel* G4VMultipleScattering::EmModel(G4int index)
00187 {
00188 G4VMscModel* p = 0;
00189 if(index >= 0 && index < G4int(mscModels.size())) { p = mscModels[index]; }
00190 return p;
00191 }
00192
00193
00194
00195 G4VEmModel*
00196 G4VMultipleScattering::GetModelByIndex(G4int idx, G4bool ver) const
00197 {
00198 return modelManager->GetModel(idx, ver);
00199 }
00200
00201
00202
00203 void
00204 G4VMultipleScattering::PreparePhysicsTable(const G4ParticleDefinition& part)
00205 {
00206 if(!firstParticle) { firstParticle = ∂ }
00207 if(part.GetParticleType() == "nucleus") {
00208 SetStepLimitType(fMinimal);
00209 SetLateralDisplasmentFlag(false);
00210 SetRangeFactor(0.2);
00211 if(&part == G4GenericIon::GenericIon()) { firstParticle = ∂ }
00212 isIon = true;
00213 }
00214
00215 emManager->PreparePhysicsTable(&part, this);
00216 currParticle = 0;
00217
00218 if(1 < verboseLevel) {
00219 G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for "
00220 << GetProcessName()
00221 << " and particle " << part.GetParticleName()
00222 << " local particle " << firstParticle->GetParticleName()
00223 << " isIon= " << isIon
00224 << G4endl;
00225 }
00226
00227 if(firstParticle == &part) {
00228
00229 InitialiseProcess(firstParticle);
00230
00231
00232 numberOfModels = modelManager->NumberOfModels();
00233 for(G4int i=0; i<numberOfModels; ++i) {
00234 G4VMscModel* msc = static_cast<G4VMscModel*>(modelManager->GetModel(i));
00235 msc->SetIonisation(0, firstParticle);
00236 if(0 == i) { currentModel = msc; }
00237 if(isIon) {
00238 msc->SetStepLimitType(fMinimal);
00239 msc->SetLateralDisplasmentFlag(false);
00240 msc->SetRangeFactor(0.2);
00241 } else {
00242 msc->SetStepLimitType(StepLimitType());
00243 msc->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
00244 msc->SetSkin(Skin());
00245 msc->SetRangeFactor(RangeFactor());
00246 msc->SetGeomFactor(GeomFactor());
00247 }
00248 msc->SetPolarAngleLimit(polarAngleLimit);
00249 G4double emax =
00250 std::min(msc->HighEnergyLimit(),emManager->MaxKinEnergy());
00251 msc->SetHighEnergyLimit(emax);
00252 }
00253
00254 modelManager->Initialise(firstParticle, G4Electron::Electron(),
00255 10.0, verboseLevel);
00256
00257 if(!safetyHelper) {
00258 safetyHelper = G4TransportationManager::GetTransportationManager()
00259 ->GetSafetyHelper();
00260 safetyHelper->InitialiseHelper();
00261 }
00262 }
00263 }
00264
00265
00266
00267 void G4VMultipleScattering::BuildPhysicsTable(const G4ParticleDefinition& part)
00268 {
00269 G4String num = part.GetParticleName();
00270 if(1 < verboseLevel) {
00271 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for "
00272 << GetProcessName()
00273 << " and particle " << num
00274 << G4endl;
00275 }
00276
00277 emManager->BuildPhysicsTable(firstParticle);
00278
00279
00280 if(1 < verboseLevel ||
00281 (0 < verboseLevel && (num == "e-" ||
00282 num == "e+" || num == "mu+" ||
00283 num == "mu-" || num == "proton"||
00284 num == "pi+" || num == "pi-" ||
00285 num == "kaon+" || num == "kaon-" ||
00286 num == "alpha" || num == "anti_proton" ||
00287 num == "GenericIon")))
00288 {
00289 G4cout << G4endl << GetProcessName()
00290 << ": for " << num
00291 << " SubType= " << GetProcessSubType()
00292 << G4endl;
00293 PrintInfo();
00294 modelManager->DumpModelList(verboseLevel);
00295 }
00296
00297 if(1 < verboseLevel) {
00298 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() done for "
00299 << GetProcessName()
00300 << " and particle " << num
00301 << G4endl;
00302 }
00303 }
00304
00305
00306
00307 void G4VMultipleScattering::PrintInfoDefinition()
00308 {
00309 if (0 < verboseLevel) {
00310 G4cout << G4endl << GetProcessName()
00311 << ": for " << firstParticle->GetParticleName()
00312 << " SubType= " << GetProcessSubType()
00313 << G4endl;
00314 PrintInfo();
00315 modelManager->DumpModelList(verboseLevel);
00316 }
00317 }
00318
00319
00320
00321 void G4VMultipleScattering::StartTracking(G4Track* track)
00322 {
00323 G4VEnergyLossProcess* eloss = 0;
00324 if(track->GetParticleDefinition() != currParticle) {
00325 currParticle = track->GetParticleDefinition();
00326 fIonisation = emManager->GetEnergyLossProcess(currParticle);
00327 eloss = fIonisation;
00328 }
00329
00330 if(1 == numberOfModels) {
00331 currentModel->StartTracking(track);
00332 if(eloss) { currentModel->SetIonisation(fIonisation, currParticle); }
00333
00334
00335 } else {
00336 for(G4int i=0; i<numberOfModels; ++i) {
00337 G4VMscModel* msc = static_cast<G4VMscModel*>(modelManager->GetModel(i));
00338 msc->StartTracking(track);
00339 if(eloss) { msc->SetIonisation(fIonisation, currParticle); }
00340 }
00341 }
00342 }
00343
00344
00345
00346 G4double G4VMultipleScattering::AlongStepGetPhysicalInteractionLength(
00347 const G4Track& track,
00348 G4double,
00349 G4double currentMinimalStep,
00350 G4double&,
00351 G4GPILSelection* selection)
00352 {
00353
00354 *selection = NotCandidateForSelection;
00355 physStepLimit = gPathLength = tPathLength = currentMinimalStep;
00356
00357 G4double ekin = track.GetKineticEnergy();
00358
00359 if(isIon) {
00360 ekin *= proton_mass_c2/track.GetParticleDefinition()->GetPDGMass();
00361 }
00362
00363
00364 if(1 < numberOfModels) {
00365 currentModel = static_cast<G4VMscModel*>(
00366 SelectModel(ekin,track.GetMaterialCutsCouple()->GetIndex()));
00367 }
00368
00369
00370 if(currentModel->IsActive(ekin) && gPathLength >= geomMin
00371 && ekin >= lowestKinEnergy) {
00372 isActive = true;
00373 tPathLength = currentModel->ComputeTruePathLengthLimit(track, gPathLength);
00374 if (tPathLength < physStepLimit) {
00375 *selection = CandidateForSelection;
00376 }
00377 } else { isActive = false; }
00378
00379
00380
00381
00382
00383
00384
00385
00386 return gPathLength;
00387 }
00388
00389
00390
00391 G4double
00392 G4VMultipleScattering::PostStepGetPhysicalInteractionLength(
00393 const G4Track&, G4double, G4ForceCondition* condition)
00394 {
00395 *condition = Forced;
00396
00397 return DBL_MAX;
00398 }
00399
00400
00401
00402 G4VParticleChange*
00403 G4VMultipleScattering::AlongStepDoIt(const G4Track& track, const G4Step& step)
00404 {
00405 fParticleChange.ProposeMomentumDirection(step.GetPostStepPoint()->GetMomentumDirection());
00406 fNewPosition = step.GetPostStepPoint()->GetPosition();
00407 fParticleChange.ProposePosition(fNewPosition);
00408 fPositionChanged = false;
00409
00410 G4double geomLength = step.GetStepLength();
00411
00412
00413 if(!isActive) {
00414 tPathLength = geomLength;
00415
00416
00417 } else {
00418 G4double range =
00419 currentModel->GetRange(currParticle,track.GetKineticEnergy(),
00420 track.GetMaterialCutsCouple());
00421
00422 G4double trueLength = currentModel->ComputeTrueStepLength(geomLength);
00423
00424
00425
00426
00427
00428
00429
00430
00431
00432
00433
00434
00435
00436 if (trueLength <= physStepLimit) {
00437 tPathLength = trueLength;
00438 } else {
00439 tPathLength = physStepLimit - 0.5*geomMin;
00440 }
00441
00442
00443 if(tPathLength + geomMin < range && tPathLength > geomMin) {
00444
00445 G4double preSafety = step.GetPreStepPoint()->GetSafety();
00446 G4double postSafety= preSafety - geomLength;
00447 G4bool safetyRecomputed = false;
00448 if( postSafety < geomMin ) {
00449 safetyRecomputed = true;
00450 postSafety = currentModel->ComputeSafety(fNewPosition,0.0);
00451 }
00452 G4ThreeVector displacement =
00453 currentModel->SampleScattering(step.GetPostStepPoint()->GetMomentumDirection(),postSafety);
00454
00455 G4double r2 = displacement.mag2();
00456
00457
00458
00459
00460 if(r2 > 0.0) {
00461
00462 fPositionChanged = true;
00463 G4double fac = 1.0;
00464
00465
00466 if(r2 > postSafety*postSafety) {
00467 if(!safetyRecomputed) {
00468 postSafety = currentModel->ComputeSafety(fNewPosition, 0.0);
00469 }
00470
00471 if(r2 > postSafety*postSafety) { fac = 0.99*postSafety/std::sqrt(r2); }
00472 }
00473
00474 fNewPosition += fac*displacement;
00475
00476 }
00477 }
00478 }
00479 fParticleChange.ProposeTrueStepLength(tPathLength);
00480
00481 return &fParticleChange;
00482 }
00483
00484
00485
00486 G4VParticleChange*
00487 G4VMultipleScattering::PostStepDoIt(const G4Track& track, const G4Step&)
00488 {
00489 fParticleChange.Initialize(track);
00490
00491 if(fPositionChanged) {
00492 safetyHelper->ReLocateWithinVolume(fNewPosition);
00493 fParticleChange.ProposePosition(fNewPosition);
00494 }
00495
00496 return &fParticleChange;
00497 }
00498
00499
00500
00501 G4double G4VMultipleScattering::GetContinuousStepLimit(
00502 const G4Track& track,
00503 G4double previousStepSize,
00504 G4double currentMinimalStep,
00505 G4double& currentSafety)
00506 {
00507 G4GPILSelection selection = NotCandidateForSelection;
00508 G4double x = AlongStepGetPhysicalInteractionLength(track,previousStepSize,
00509 currentMinimalStep,
00510 currentSafety, &selection);
00511 return x;
00512 }
00513
00514
00515
00516 G4double G4VMultipleScattering::ContinuousStepLimit(
00517 const G4Track& track,
00518 G4double previousStepSize,
00519 G4double currentMinimalStep,
00520 G4double& currentSafety)
00521 {
00522 return GetContinuousStepLimit(track,previousStepSize,currentMinimalStep,
00523 currentSafety);
00524 }
00525
00526
00527
00528 G4double G4VMultipleScattering::GetMeanFreePath(
00529 const G4Track&, G4double, G4ForceCondition* condition)
00530 {
00531 *condition = Forced;
00532 return DBL_MAX;
00533 }
00534
00535
00536
00537 G4bool
00538 G4VMultipleScattering::StorePhysicsTable(const G4ParticleDefinition* part,
00539 const G4String& directory,
00540 G4bool ascii)
00541 {
00542 G4bool yes = true;
00543 if(part != firstParticle) { return yes; }
00544 G4int nmod = modelManager->NumberOfModels();
00545 const G4String ss[4] = {"1","2","3","4"};
00546 for(G4int i=0; i<nmod; ++i) {
00547 G4VEmModel* msc = modelManager->GetModel(i);
00548 yes = true;
00549 G4PhysicsTable* table = msc->GetCrossSectionTable();
00550 if (table) {
00551 G4int j = std::min(i,3);
00552 G4String name =
00553 GetPhysicsTableFileName(part,directory,"LambdaMod"+ss[j],ascii);
00554 yes = table->StorePhysicsTable(name,ascii);
00555
00556 if ( yes ) {
00557 if ( verboseLevel>0 ) {
00558 G4cout << "Physics table are stored for " << part->GetParticleName()
00559 << " and process " << GetProcessName()
00560 << " with a name <" << name << "> " << G4endl;
00561 }
00562 } else {
00563 G4cout << "Fail to store Physics Table for " << part->GetParticleName()
00564 << " and process " << GetProcessName()
00565 << " in the directory <" << directory
00566 << "> " << G4endl;
00567 }
00568 }
00569 }
00570 return yes;
00571 }
00572
00573
00574
00575 G4bool
00576 G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition*,
00577 const G4String&,
00578 G4bool)
00579 {
00580 return true;
00581 }
00582
00583
00584
00585 void G4VMultipleScattering::SetIonisation(G4VEnergyLossProcess* p)
00586 {
00587 for(G4int i=0; i<numberOfModels; ++i) {
00588 G4VMscModel* msc = static_cast<G4VMscModel*>(modelManager->GetModel(i));
00589 msc->SetIonisation(p, firstParticle);
00590 }
00591 }
00592
00593
00594