G4Decay.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 //
00027 // $Id$
00028 //
00029 // 
00030 // --------------------------------------------------------------
00031 //      GEANT 4 class implementation file
00032 //
00033 //      History: first implementation, based on object model of
00034 //      2nd December 1995, G.Cosmo
00035 //      7 July 1996 H.Kurashige
00036 // ------------------------------------------------------------
00037 //   remove  BuildPhysicsTable()  28 Nov. 1997  H.Kurashige
00038 //   change  DBL_EPSIRON to DBL_MIN 14 Dec. 1997  H.Kurashige
00039 //   modified for new ParticleChange 12 Mar. 1998  H.Kurashige
00040 //   modified for "GoodForTrackingFlag" 19 June 1998  H.Kurashige
00041 //   rename thePhysicsTable to aPhyscisTable 2 Aug. 1998 H.Kurashige
00042 //   modified IsApplicable in order to protect the decay from registered 
00043 //   to resonances    12 Dec. 1998   H.Kurashige 
00044 //   remove G4ParticleMomentum  6 Feb. 99 H.Kurashige
00045 //   modified  IsApplicable to activate G4Decay for resonances  1 Mar. 00 H.Kurashige 
00046 //   Add External Decayer         23 Feb. 2001  H.Kurashige
00047 //   change LowestBinValue,HighestBinValue and TotBin(200) 9 Feb. 2002
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 // constructor
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   // set Process Sub Type
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    // check if the particle is stable?
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    // returns the mean free path in GEANT4 internal units
00105    G4double meanlife;
00106 
00107    // get particle 
00108    const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00109    const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
00110    G4double aLife = aParticleDef->GetPDGLifeTime();
00111 
00112    // check if the particle is stable?
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    // get particle 
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     // returns the mean free path in GEANT4 internal units
00139    G4double pathlength;
00140    G4double aCtau = c_light * aLife;
00141 
00142    // check if the particle is stable?
00143    if (aParticleDef->GetPDGStable()) {
00144      pathlength = DBL_MAX;
00145 
00146    //check if the particle has very short life time ?
00147    } else if (aCtau < DBL_MIN) { 
00148      pathlength =  DBL_MIN;
00149  
00150    } else {
00151     //calculate the mean free path 
00152     // by using normalized kinetic energy (= Ekin/mass)
00153      G4double   rKineticEnergy = aParticle->GetKineticEnergy()/aMass; 
00154      if ( rKineticEnergy > HighestValue) {
00155        // gamma >>  1
00156        pathlength = ( rKineticEnergy + 1.0)* aCtau;
00157      } else if ( rKineticEnergy < DBL_MIN ) {
00158        // too slow particle
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        // beta <1 
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   // The DecayIt() method returns by pointer a particle-change object.
00183   // Units are expressed in GEANT4 internal units.
00184 
00185   //   Initialize ParticleChange
00186   //     all members of G4VParticleChange are set to equal to 
00187   //     corresponding member in G4Track
00188   fParticleChangeForDecay.Initialize(aTrack);
00189 
00190   // get particle 
00191   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
00192   const G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
00193 
00194   // check if  the particle is stable
00195   if (aParticleDef->GetPDGStable()) return &fParticleChangeForDecay ;
00196  
00197 
00198   //check if thePreAssignedDecayProducts exists
00199   const G4DecayProducts* o_products = (aParticle->GetPreAssignedDecayProducts());
00200   G4bool isPreAssigned = (o_products != 0);   
00201   G4DecayProducts* products = 0;
00202 
00203   // decay table
00204   G4DecayTable   *decaytable = aParticleDef->GetDecayTable();
00205  
00206   // check if external decayer exists
00207   G4bool isExtDecayer = (decaytable == 0) && (pExtDecayer !=0);
00208 
00209   // Error due to NO Decay Table 
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     // Kill the parent particle
00221     fParticleChangeForDecay.ProposeTrackStatus( fStopAndKill ) ;
00222     fParticleChangeForDecay.ProposeLocalEnergyDeposit(0.0); 
00223     
00224     ClearNumberOfInteractionLengthLeft();
00225     return &fParticleChangeForDecay ;
00226   }
00227 
00228   if (isPreAssigned) {
00229     // copy decay products 
00230     products = new G4DecayProducts(*o_products); 
00231   } else if ( isExtDecayer ) {
00232     // decay according to external decayer
00233     products = pExtDecayer->ImportDecayProducts(aTrack);
00234   } else {
00235     // decay acoording to decay table
00236     // choose a decay channel
00237     G4VDecayChannel *decaychannel = decaytable->SelectADecayChannel();
00238     if (decaychannel == 0 ){
00239       // decay channel not found
00240       G4Exception("G4Decay::DoIt", "DECAY003", FatalException,
00241                   " can not determine decay channel ");
00242     } else {
00243       // execute DecayIt() 
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   // get parent particle information ...................................
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   //boost all decay products to laboratory frame
00285   G4double energyDeposit = 0.0;
00286   G4double finalGlobalTime = aTrack.GetGlobalTime();
00287   G4double finalLocalTime = aTrack.GetLocalTime();
00288   if (aTrack.GetTrackStatus() == fStopButAlive ){
00289     // AtRest case
00290     finalGlobalTime += fRemainderLifeTime;
00291     finalLocalTime += fRemainderLifeTime;
00292     energyDeposit += aParticle->GetKineticEnergy();
00293     if (isPreAssigned) products->Boost( ParentEnergy, ParentDirection);
00294   } else {
00295     // PostStep case
00296     if (!isExtDecayer) products->Boost( ParentEnergy, ParentDirection);
00297   }
00298 
00299    // set polarization for daughter particles
00300    DaughterPolarization(aTrack, products);
00301 
00302 
00303   //add products in fParticleChangeForDecay
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      // get current position of the track
00324      currentPosition = aTrack.GetPosition();
00325      // create a new track object
00326      G4Track* secondary = new G4Track( products->PopProducts(),
00327                                       finalGlobalTime ,
00328                                       currentPosition );
00329      // switch on good for tracking flag
00330      secondary->SetGoodForTrackingFlag();
00331      secondary->SetTouchableHandle(thand);
00332      // add the secondary track in the List
00333      fParticleChangeForDecay.AddSecondary(secondary);
00334   }
00335   delete products;
00336 
00337   // Kill the parent particle
00338   fParticleChangeForDecay.ProposeTrackStatus( fStopAndKill ) ;
00339   fParticleChangeForDecay.ProposeLocalEnergyDeposit(energyDeposit); 
00340   fParticleChangeForDecay.ProposeLocalTime( finalLocalTime );
00341 
00342   // Clear NumberOfInteractionLengthLeft
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   // Clear NumberOfInteractionLengthLeft
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    // condition is set to "Not Forced"
00379   *condition = NotForced;
00380 
00381    // pre-assigned Decay time
00382   G4double pTime = track.GetDynamicParticle()->GetPreAssignedDecayProperTime();
00383   G4double aLife = track.GetDynamicParticle()->GetDefinition()->GetPDGLifeTime();
00384 
00385   if (pTime < 0.) {
00386     // normal case 
00387     if ( previousStepSize > 0.0){
00388       // subtract NumberOfInteractionLengthLeft 
00389       SubtractNumberOfInteractionLengthLeft(previousStepSize);
00390       if(theNumberOfInteractionLengthLeft<0.){
00391         theNumberOfInteractionLengthLeft=perMillion;
00392       }
00393       fRemainderLifeTime = theNumberOfInteractionLengthLeft*aLife;
00394     }
00395     // get mean free path
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     //pre-assigned Decay time case
00418     // reminder proper time
00419     fRemainderLifeTime = pTime - track.GetProperTime();
00420     if (fRemainderLifeTime <= 0.0) fRemainderLifeTime = DBL_MIN;
00421     
00422     G4double  rvalue=0.0; 
00423     // use pre-assigned Decay time to determine PIL
00424     if (aLife>0.0) {
00425       // ordinary particle
00426       rvalue = (fRemainderLifeTime/aLife)*GetMeanFreePath(track, previousStepSize, condition);
00427     } else {
00428      // shortlived particle
00429       rvalue = c_light * fRemainderLifeTime;
00430      // by using normalized kinetic energy (= Ekin/mass)
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      // condition is set to "Not Forced"
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   // set Process Sub Type
00463   if ( pExtDecayer !=0 ) {
00464     SetProcessSubType(static_cast<int>(DECAY_External));
00465   }
00466 }

Generated on Mon May 27 17:47:59 2013 for Geant4 by  doxygen 1.4.7