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 #include "G4Decay.hh"
00051
00052 #include "G4PhysicalConstants.hh"
00053 #include "G4SystemOfUnits.hh"
00054 #include "G4DynamicParticle.hh"
00055 #include "G4DecayProducts.hh"
00056 #include "G4DecayTable.hh"
00057 #include "G4VDecayChannel.hh"
00058 #include "G4PhysicsLogVector.hh"
00059 #include "G4ParticleChangeForDecay.hh"
00060 #include "G4VExtDecayer.hh"
00061
00062
00063 G4Decay::G4Decay(const G4String& processName)
00064 :G4VRestDiscreteProcess(processName, fDecay),
00065 verboseLevel(1),
00066 HighestValue(20.0),
00067 fRemainderLifeTime(-1.0),
00068 pExtDecayer(0)
00069 {
00070
00071 SetProcessSubType(static_cast<int>(DECAY));
00072
00073 #ifdef G4VERBOSE
00074 if (GetVerboseLevel()>1) {
00075 G4cout << "G4Decay constructor " << " Name:" << processName << G4endl;
00076 }
00077 #endif
00078
00079 pParticleChange = &fParticleChangeForDecay;
00080 }
00081
00082 G4Decay::~G4Decay()
00083 {
00084 if (pExtDecayer) {
00085 delete pExtDecayer;
00086 }
00087 }
00088
00089 G4bool G4Decay::IsApplicable(const G4ParticleDefinition& aParticleType)
00090 {
00091
00092 if (aParticleType.GetPDGLifeTime() <0.0) {
00093 return false;
00094 } else if (aParticleType.GetPDGMass() <= 0.0*MeV) {
00095 return false;
00096 } else {
00097 return true;
00098 }
00099 }
00100
00101 G4double G4Decay::GetMeanLifeTime(const G4Track& aTrack ,
00102 G4ForceCondition*)
00103 {
00104
00105 G4double meanlife;
00106
00107
00108 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00109 const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
00110 G4double aLife = aParticleDef->GetPDGLifeTime();
00111
00112
00113 if (aParticleDef->GetPDGStable()) {
00114 meanlife = DBL_MAX;
00115
00116 } else {
00117 meanlife = aLife;
00118 }
00119
00120 #ifdef G4VERBOSE
00121 if (GetVerboseLevel()>1) {
00122 G4cout << "mean life time: "<< meanlife/ns << "[ns]" << G4endl;
00123 }
00124 #endif
00125
00126 return meanlife;
00127 }
00128
00129 G4double G4Decay::GetMeanFreePath(const G4Track& aTrack,G4double, G4ForceCondition*)
00130 {
00131
00132 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00133 const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
00134 G4double aMass = aParticle->GetMass();
00135 G4double aLife = aParticleDef->GetPDGLifeTime();
00136
00137
00138
00139 G4double pathlength;
00140 G4double aCtau = c_light * aLife;
00141
00142
00143 if (aParticleDef->GetPDGStable()) {
00144 pathlength = DBL_MAX;
00145
00146
00147 } else if (aCtau < DBL_MIN) {
00148 pathlength = DBL_MIN;
00149
00150 } else {
00151
00152
00153 G4double rKineticEnergy = aParticle->GetKineticEnergy()/aMass;
00154 if ( rKineticEnergy > HighestValue) {
00155
00156 pathlength = ( rKineticEnergy + 1.0)* aCtau;
00157 } else if ( rKineticEnergy < DBL_MIN ) {
00158
00159 #ifdef G4VERBOSE
00160 if (GetVerboseLevel()>1) {
00161 G4cout << "G4Decay::GetMeanFreePath() !!particle stops!!";
00162 G4cout << aParticleDef->GetParticleName() << G4endl;
00163 G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
00164 }
00165 #endif
00166 pathlength = DBL_MIN;
00167 } else {
00168
00169 pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ;
00170 }
00171 }
00172 return pathlength;
00173 }
00174
00175 void G4Decay::BuildPhysicsTable(const G4ParticleDefinition&)
00176 {
00177 return;
00178 }
00179
00180 G4VParticleChange* G4Decay::DecayIt(const G4Track& aTrack, const G4Step& )
00181 {
00182
00183
00184
00185
00186
00187
00188 fParticleChangeForDecay.Initialize(aTrack);
00189
00190
00191 const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00192 const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
00193
00194
00195 if (aParticleDef->GetPDGStable()) return &fParticleChangeForDecay ;
00196
00197
00198
00199 const G4DecayProducts* o_products = (aParticle->GetPreAssignedDecayProducts());
00200 G4bool isPreAssigned = (o_products != 0);
00201 G4DecayProducts* products = 0;
00202
00203
00204 G4DecayTable *decaytable = aParticleDef->GetDecayTable();
00205
00206
00207 G4bool isExtDecayer = (decaytable == 0) && (pExtDecayer !=0);
00208
00209
00210 if ( (decaytable == 0) && !isExtDecayer &&!isPreAssigned ){
00211 if (GetVerboseLevel()>0) {
00212 G4cout << "G4Decay::DoIt : decay table not defined for ";
00213 G4cout << aParticle->GetDefinition()->GetParticleName()<< G4endl;
00214 }
00215 G4Exception( "G4Decay::DecayIt ",
00216 "DECAY101",JustWarning,
00217 "Decay table is not defined");
00218
00219 fParticleChangeForDecay.SetNumberOfSecondaries(0);
00220
00221 fParticleChangeForDecay.ProposeTrackStatus( fStopAndKill ) ;
00222 fParticleChangeForDecay.ProposeLocalEnergyDeposit(0.0);
00223
00224 ClearNumberOfInteractionLengthLeft();
00225 return &fParticleChangeForDecay ;
00226 }
00227
00228 if (isPreAssigned) {
00229
00230 products = new G4DecayProducts(*o_products);
00231 } else if ( isExtDecayer ) {
00232
00233 products = pExtDecayer->ImportDecayProducts(aTrack);
00234 } else {
00235
00236
00237 G4VDecayChannel *decaychannel = decaytable->SelectADecayChannel();
00238 if (decaychannel == 0 ){
00239
00240 G4Exception("G4Decay::DoIt", "DECAY003", FatalException,
00241 " can not determine decay channel ");
00242 } else {
00243
00244 #ifdef G4VERBOSE
00245 G4int temp = decaychannel->GetVerboseLevel();
00246 if (GetVerboseLevel()>1) {
00247 G4cout << "G4Decay::DoIt : selected decay channel addr:" << decaychannel <<G4endl;
00248 decaychannel->SetVerboseLevel(GetVerboseLevel());
00249 }
00250 #endif
00251 products = decaychannel->DecayIt(aParticle->GetMass());
00252 #ifdef G4VERBOSE
00253 if (GetVerboseLevel()>1) {
00254 decaychannel->SetVerboseLevel(temp);
00255 }
00256 #endif
00257 #ifdef G4VERBOSE
00258 if (GetVerboseLevel()>2) {
00259 if (! products->IsChecked() ) products->DumpInfo();
00260 }
00261 #endif
00262 }
00263 }
00264
00265
00266 G4double ParentEnergy = aParticle->GetTotalEnergy();
00267 G4double ParentMass = aParticle->GetMass();
00268 if (ParentEnergy < ParentMass) {
00269 if (GetVerboseLevel()>0) {
00270 G4cout << "G4Decay::DoIt : Total Energy is less than its mass" << G4endl;
00271 G4cout << " Particle: " << aParticle->GetDefinition()->GetParticleName();
00272 G4cout << " Energy:" << ParentEnergy/MeV << "[MeV]";
00273 G4cout << " Mass:" << ParentMass/MeV << "[MeV]";
00274 G4cout << G4endl;
00275 }
00276 G4Exception( "G4Decay::DecayIt ",
00277 "DECAY102",JustWarning,
00278 "Total Energy is less than its mass");
00279 ParentEnergy = ParentMass;
00280 }
00281
00282 G4ThreeVector ParentDirection(aParticle->GetMomentumDirection());
00283
00284
00285 G4double energyDeposit = 0.0;
00286 G4double finalGlobalTime = aTrack.GetGlobalTime();
00287 G4double finalLocalTime = aTrack.GetLocalTime();
00288 if (aTrack.GetTrackStatus() == fStopButAlive ){
00289
00290 finalGlobalTime += fRemainderLifeTime;
00291 finalLocalTime += fRemainderLifeTime;
00292 energyDeposit += aParticle->GetKineticEnergy();
00293 if (isPreAssigned) products->Boost( ParentEnergy, ParentDirection);
00294 } else {
00295
00296 if (!isExtDecayer) products->Boost( ParentEnergy, ParentDirection);
00297 }
00298
00299
00300 DaughterPolarization(aTrack, products);
00301
00302
00303
00304 G4int numberOfSecondaries = products->entries();
00305 fParticleChangeForDecay.SetNumberOfSecondaries(numberOfSecondaries);
00306 #ifdef G4VERBOSE
00307 if (GetVerboseLevel()>1) {
00308 G4cout << "G4Decay::DoIt : Decay vertex :";
00309 G4cout << " Time: " << finalGlobalTime/ns << "[ns]";
00310 G4cout << " X:" << (aTrack.GetPosition()).x() /cm << "[cm]";
00311 G4cout << " Y:" << (aTrack.GetPosition()).y() /cm << "[cm]";
00312 G4cout << " Z:" << (aTrack.GetPosition()).z() /cm << "[cm]";
00313 G4cout << G4endl;
00314 G4cout << "G4Decay::DoIt : decay products in Lab. Frame" << G4endl;
00315 products->DumpInfo();
00316 }
00317 #endif
00318 G4int index;
00319 G4ThreeVector currentPosition;
00320 const G4TouchableHandle thand = aTrack.GetTouchableHandle();
00321 for (index=0; index < numberOfSecondaries; index++)
00322 {
00323
00324 currentPosition = aTrack.GetPosition();
00325
00326 G4Track* secondary = new G4Track( products->PopProducts(),
00327 finalGlobalTime ,
00328 currentPosition );
00329
00330 secondary->SetGoodForTrackingFlag();
00331 secondary->SetTouchableHandle(thand);
00332
00333 fParticleChangeForDecay.AddSecondary(secondary);
00334 }
00335 delete products;
00336
00337
00338 fParticleChangeForDecay.ProposeTrackStatus( fStopAndKill ) ;
00339 fParticleChangeForDecay.ProposeLocalEnergyDeposit(energyDeposit);
00340 fParticleChangeForDecay.ProposeLocalTime( finalLocalTime );
00341
00342
00343 ClearNumberOfInteractionLengthLeft();
00344
00345 return &fParticleChangeForDecay ;
00346 }
00347
00348 void G4Decay::DaughterPolarization(const G4Track& , G4DecayProducts* )
00349 {
00350 }
00351
00352
00353
00354 void G4Decay::StartTracking(G4Track*)
00355 {
00356 currentInteractionLength = -1.0;
00357 ResetNumberOfInteractionLengthLeft();
00358
00359 fRemainderLifeTime = -1.0;
00360 }
00361
00362 void G4Decay::EndTracking()
00363 {
00364
00365 ClearNumberOfInteractionLengthLeft();
00366
00367 currentInteractionLength = -1.0;
00368 }
00369
00370
00371 G4double G4Decay::PostStepGetPhysicalInteractionLength(
00372 const G4Track& track,
00373 G4double previousStepSize,
00374 G4ForceCondition* condition
00375 )
00376 {
00377
00378
00379 *condition = NotForced;
00380
00381
00382 G4double pTime = track.GetDynamicParticle()->GetPreAssignedDecayProperTime();
00383 G4double aLife = track.GetDynamicParticle()->GetDefinition()->GetPDGLifeTime();
00384
00385 if (pTime < 0.) {
00386
00387 if ( previousStepSize > 0.0){
00388
00389 SubtractNumberOfInteractionLengthLeft(previousStepSize);
00390 if(theNumberOfInteractionLengthLeft<0.){
00391 theNumberOfInteractionLengthLeft=perMillion;
00392 }
00393 fRemainderLifeTime = theNumberOfInteractionLengthLeft*aLife;
00394 }
00395
00396 currentInteractionLength = GetMeanFreePath(track, previousStepSize, condition);
00397
00398 #ifdef G4VERBOSE
00399 if ((currentInteractionLength <=0.0) || (verboseLevel>2)){
00400 G4cout << "G4Decay::PostStepGetPhysicalInteractionLength " << G4endl;
00401 track.GetDynamicParticle()->DumpInfo();
00402 G4cout << " in Material " << track.GetMaterial()->GetName() <<G4endl;
00403 G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" <<G4endl;
00404 }
00405 #endif
00406
00407 G4double value;
00408 if (currentInteractionLength <DBL_MAX) {
00409 value = theNumberOfInteractionLengthLeft * currentInteractionLength;
00410 } else {
00411 value = DBL_MAX;
00412 }
00413
00414 return value;
00415
00416 } else {
00417
00418
00419 fRemainderLifeTime = pTime - track.GetProperTime();
00420 if (fRemainderLifeTime <= 0.0) fRemainderLifeTime = DBL_MIN;
00421
00422 G4double rvalue=0.0;
00423
00424 if (aLife>0.0) {
00425
00426 rvalue = (fRemainderLifeTime/aLife)*GetMeanFreePath(track, previousStepSize, condition);
00427 } else {
00428
00429 rvalue = c_light * fRemainderLifeTime;
00430
00431 G4double aMass = track.GetDynamicParticle()->GetMass();
00432 rvalue *= track.GetDynamicParticle()->GetTotalMomentum()/aMass;
00433 }
00434 return rvalue;
00435 }
00436 }
00437
00438 G4double G4Decay::AtRestGetPhysicalInteractionLength(
00439 const G4Track& track,
00440 G4ForceCondition* condition
00441 )
00442 {
00443
00444 *condition = NotForced;
00445
00446 G4double pTime = track.GetDynamicParticle()->GetPreAssignedDecayProperTime();
00447 if (pTime >= 0.) {
00448 fRemainderLifeTime = pTime - track.GetProperTime();
00449 if (fRemainderLifeTime <= 0.0) fRemainderLifeTime = DBL_MIN;
00450 } else {
00451 fRemainderLifeTime =
00452 theNumberOfInteractionLengthLeft * GetMeanLifeTime(track, condition);
00453 }
00454 return fRemainderLifeTime;
00455 }
00456
00457
00458 void G4Decay::SetExtDecayer(G4VExtDecayer* val)
00459 {
00460 pExtDecayer = val;
00461
00462
00463 if ( pExtDecayer !=0 ) {
00464 SetProcessSubType(static_cast<int>(DECAY_External));
00465 }
00466 }