source: trunk/source/processes/decay/src/G4Decay.cc @ 846

Last change on this file since 846 was 819, checked in by garnier, 16 years ago

import all except CVS
File size: 14.8 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer                                           *
4// *                                                                  *
5// * The  Geant4 software  is  copyright of the Copyright Holders  of *
6// * the Geant4 Collaboration.  It is provided  under  the terms  and *
7// * conditions of the Geant4 Software License,  included in the file *
8// * LICENSE and available at  http://cern.ch/geant4/license .  These *
9// * include a list of copyright holders.                             *
10// *                                                                  *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work  make  any representation or  warranty, express or implied, *
14// * regarding  this  software system or assume any liability for its *
15// * use.  Please see the license in the file  LICENSE  and URL above *
16// * for the full disclaimer and the limitation of liability.         *
17// *                                                                  *
18// * This  code  implementation is the result of  the  scientific and *
19// * technical work of the GEANT4 collaboration.                      *
20// * By using,  copying,  modifying or  distributing the software (or *
21// * any work based  on the software)  you  agree  to acknowledge its *
22// * use  in  resulting  scientific  publications,  and indicate your *
23// * acceptance of all terms of the Geant4 Software license.          *
24// ********************************************************************
25//
26//
27// $Id: G4Decay.cc,v 1.27.2.1 2008/04/17 08:59:24 kurasige Exp $
28// GEANT4 tag $Name: geant4-09-01-patch-02 $
29//
30//
31// --------------------------------------------------------------
32//      GEANT 4 class implementation file
33//
34//      History: first implementation, based on object model of
35//      2nd December 1995, G.Cosmo
36//      7 July 1996 H.Kurashige
37// ------------------------------------------------------------
38//   remove  BuildPhysicsTable()  28 Nov. 1997  H.Kurashige
39//   change  DBL_EPSIRON to DBL_MIN 14 Dec. 1997  H.Kurashige
40//   modified for new ParticleChange 12 Mar. 1998  H.Kurashige
41//   modified for "GoodForTrackingFlag" 19 June 1998  H.Kurashige
42//   rename thePhysicsTable to aPhyscisTable 2 Aug. 1998 H.Kurashige
43//   modified IsApplicable in order to protect the decay from registered
44//   to resonances    12 Dec. 1998   H.Kurashige
45//   remove G4ParticleMomentum  6 Feb. 99 H.Kurashige
46//   modified  IsApplicable to activate G4Decay for resonances  1 Mar. 00 H.Kurashige
47//   Add External Decayer         23 Feb. 2001  H.Kurashige
48//   change LowestBinValue,HighestBinValue and TotBin(200) 9 Feb. 2002
49//
50
51#include "G4Decay.hh"
52#include "G4DynamicParticle.hh"
53#include "G4DecayProducts.hh"
54#include "G4DecayTable.hh"
55#include "G4PhysicsLogVector.hh"
56#include "G4ParticleChangeForDecay.hh"
57#include "G4VExtDecayer.hh"
58
59// constructor
60G4Decay::G4Decay(const G4String& processName)
61                               :G4VRestDiscreteProcess(processName, fDecay),
62                                verboseLevel(1),
63                                HighestValue(20.0),
64                                pExtDecayer(0)
65{
66#ifdef G4VERBOSE
67  if (GetVerboseLevel()>1) {
68    G4cout << "G4Decay  constructor " << "  Name:" << processName << G4endl;
69  }
70#endif
71  pParticleChange = &fParticleChangeForDecay;
72}
73
74G4Decay::~G4Decay()
75{
76  if (pExtDecayer) {
77    delete pExtDecayer;
78  }
79}
80
81G4bool G4Decay::IsApplicable(const G4ParticleDefinition& aParticleType)
82{
83   // check if the particle is stable?
84   if (aParticleType.GetPDGLifeTime() <0.0) {
85     return false;
86   } else if (aParticleType.GetPDGMass() <= 0.0*MeV) {
87     return false;
88   } else {
89     return true; 
90   }
91}
92
93G4double G4Decay::GetMeanLifeTime(const G4Track& aTrack  ,
94                                  G4ForceCondition*)
95{
96   // returns the mean free path in GEANT4 internal units
97   G4double meanlife;
98
99   // get particle
100   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
101   G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
102   G4double aLife = aParticleDef->GetPDGLifeTime();
103
104   // check if the particle is stable?
105   if (aParticleDef->GetPDGStable()) {
106     meanlife = DBL_MAX;
107   
108   } else {
109     meanlife = aLife;
110   }
111
112#ifdef G4VERBOSE
113   if (GetVerboseLevel()>1) {
114     G4cout << "mean life time: "<< meanlife/ns << "[ns]" << G4endl;
115   }
116#endif
117
118   return  meanlife;
119}
120
121G4double G4Decay::GetMeanFreePath(const G4Track& aTrack,G4double, G4ForceCondition*)
122{
123   // get particle
124   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
125   G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
126   G4double aMass = aParticle->GetMass();
127   G4double aLife = aParticleDef->GetPDGLifeTime();
128
129
130    // returns the mean free path in GEANT4 internal units
131   G4double pathlength;
132   G4double aCtau = c_light * aLife;
133
134   // check if the particle is stable?
135   if (aParticleDef->GetPDGStable()) {
136     pathlength = DBL_MAX;
137
138   //check if the particle has very short life time ?
139   } else if (aCtau < DBL_MIN) { 
140     pathlength =  DBL_MIN;
141 
142   } else {
143    //calculate the mean free path
144    // by using normalized kinetic energy (= Ekin/mass)
145     G4double   rKineticEnergy = aParticle->GetKineticEnergy()/aMass; 
146     if ( rKineticEnergy > HighestValue) {
147       // gamma >>  1
148       pathlength = ( rKineticEnergy + 1.0)* aCtau;
149     } else if ( rKineticEnergy < DBL_MIN ) {
150       // too slow particle
151#ifdef G4VERBOSE
152       if (GetVerboseLevel()>1) {
153         G4cout << "G4Decay::GetMeanFreePath()   !!particle stops!!";
154         G4cout << aParticleDef->GetParticleName() << G4endl;
155         G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
156       }
157#endif
158       pathlength = DBL_MIN;
159     } else {
160       // beta <1
161       pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ;
162     }
163   }
164  return  pathlength;
165}
166
167void G4Decay::BuildPhysicsTable(const G4ParticleDefinition&)
168{
169  return;
170}
171
172G4VParticleChange* G4Decay::DecayIt(const G4Track& aTrack, const G4Step& )
173{
174  // The DecayIt() method returns by pointer a particle-change object.
175  // Units are expressed in GEANT4 internal units.
176
177  //   Initialize ParticleChange
178  //     all members of G4VParticleChange are set to equal to
179  //     corresponding member in G4Track
180  fParticleChangeForDecay.Initialize(aTrack);
181
182  // get particle
183  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
184  G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
185
186  // check if  the particle is stable
187  if (aParticleDef->GetPDGStable()) return &fParticleChangeForDecay ;
188 
189
190  //check if thePreAssignedDecayProducts exists
191  const G4DecayProducts* o_products = (aParticle->GetPreAssignedDecayProducts());
192  G4bool isPreAssigned = (o_products != 0);   
193  G4DecayProducts* products = 0;
194
195  // decay table
196  G4DecayTable   *decaytable = aParticleDef->GetDecayTable();
197 
198  // check if external decayer exists
199  G4bool isExtDecayer = (decaytable == 0) && (pExtDecayer !=0);
200
201  // Error due to NO Decay Table
202  if ( (decaytable == 0) && !isExtDecayer &&!isPreAssigned ){
203    if (GetVerboseLevel()>0) {
204      G4cout <<  "G4Decay::DoIt  : decay table not defined  for ";
205      G4cout << aParticle->GetDefinition()->GetParticleName()<< G4endl;
206    }
207    G4Exception( "G4Decay::DecayIt ",
208                 "No Decay Table",JustWarning, 
209                 "Decay table is not defined");
210
211    fParticleChangeForDecay.SetNumberOfSecondaries(0);
212    // Kill the parent particle
213    fParticleChangeForDecay.ProposeTrackStatus( fStopAndKill ) ;
214    fParticleChangeForDecay.ProposeLocalEnergyDeposit(0.0); 
215   
216    ClearNumberOfInteractionLengthLeft();
217    return &fParticleChangeForDecay ;
218  }
219
220  if (isPreAssigned) {
221    // copy decay products
222    products = new G4DecayProducts(*o_products); 
223  } else if ( isExtDecayer ) {
224    // decay according to external decayer
225    products = pExtDecayer->ImportDecayProducts(aTrack);
226  } else {
227    // decay acoording to decay table
228    // choose a decay channel
229    G4VDecayChannel *decaychannel = decaytable->SelectADecayChannel();
230    if (decaychannel == 0 ){
231      // decay channel not found
232      G4Exception("G4Decay::DoIt  : can not determine decay channel ");
233    } else {
234      // execute DecayIt()
235#ifdef G4VERBOSE
236      G4int temp = decaychannel->GetVerboseLevel();
237      if (GetVerboseLevel()>1) {
238        G4cout << "G4Decay::DoIt  : selected decay channel  addr:" << decaychannel <<G4endl;
239        decaychannel->SetVerboseLevel(GetVerboseLevel());
240      }
241#endif
242      products = decaychannel->DecayIt(aParticle->GetMass());
243#ifdef G4VERBOSE
244      if (GetVerboseLevel()>1) {
245        decaychannel->SetVerboseLevel(temp);
246      }
247#endif
248#ifdef G4VERBOSE
249      if (GetVerboseLevel()>2) {
250        if (! products->IsChecked() ) products->DumpInfo();
251      }
252#endif
253    }
254  }
255 
256  // get parent particle information ...................................
257  G4double   ParentEnergy  = aParticle->GetTotalEnergy();
258  G4double   ParentMass    = aParticle->GetMass();
259  if (ParentEnergy < ParentMass) {
260    ParentEnergy = ParentMass;
261    if (GetVerboseLevel()>0) {
262      G4cout << "G4Decay::DoIt  : Total Energy is less than its mass" << G4endl;
263      G4cout << " Particle: " << aParticle->GetDefinition()->GetParticleName();
264      G4cout << " Energy:"    << ParentEnergy/MeV << "[MeV]";
265      G4cout << " Mass:"    << ParentMass/MeV << "[MeV]";
266      G4cout << G4endl;
267    }
268  }
269
270  G4ThreeVector ParentDirection(aParticle->GetMomentumDirection());
271
272  //boost all decay products to laboratory frame
273  G4double energyDeposit = 0.0;
274  G4double finalGlobalTime = aTrack.GetGlobalTime();
275  if (aTrack.GetTrackStatus() == fStopButAlive ){
276    // AtRest case
277    finalGlobalTime += fRemainderLifeTime;
278    energyDeposit += aParticle->GetKineticEnergy();
279    if (isPreAssigned) products->Boost( ParentEnergy, ParentDirection);
280  } else {
281    // PostStep case
282    if (!isExtDecayer) products->Boost( ParentEnergy, ParentDirection);
283  }
284
285   // set polarization for daughter particles
286   DaughterPolarization(aTrack, products);
287
288
289  //add products in fParticleChangeForDecay
290  G4int numberOfSecondaries = products->entries();
291  fParticleChangeForDecay.SetNumberOfSecondaries(numberOfSecondaries);
292#ifdef G4VERBOSE
293  if (GetVerboseLevel()>1) {
294    G4cout << "G4Decay::DoIt  : Decay vertex :";
295    G4cout << " Time: " << finalGlobalTime/ns << "[ns]";
296    G4cout << " X:" << (aTrack.GetPosition()).x() /cm << "[cm]";
297    G4cout << " Y:" << (aTrack.GetPosition()).y() /cm << "[cm]";
298    G4cout << " Z:" << (aTrack.GetPosition()).z() /cm << "[cm]";
299    G4cout << G4endl;
300    G4cout << "G4Decay::DoIt  : decay products in Lab. Frame" << G4endl;
301    products->DumpInfo();
302  }
303#endif
304  G4int index;
305  G4ThreeVector currentPosition;
306  const G4TouchableHandle thand = aTrack.GetTouchableHandle();
307  for (index=0; index < numberOfSecondaries; index++)
308  {
309     // get current position of the track
310     currentPosition = aTrack.GetPosition();
311     // create a new track object
312     G4Track* secondary = new G4Track( products->PopProducts(),
313                                      finalGlobalTime ,
314                                      currentPosition );
315     // switch on good for tracking flag
316     secondary->SetGoodForTrackingFlag();
317     secondary->SetTouchableHandle(thand);
318     // add the secondary track in the List
319     fParticleChangeForDecay.AddSecondary(secondary);
320  }
321  delete products;
322
323  // Kill the parent particle
324  fParticleChangeForDecay.ProposeTrackStatus( fStopAndKill ) ;
325  fParticleChangeForDecay.ProposeLocalEnergyDeposit(energyDeposit); 
326  fParticleChangeForDecay.ProposeGlobalTime( finalGlobalTime );
327  // Clear NumberOfInteractionLengthLeft
328  ClearNumberOfInteractionLengthLeft();
329
330  return &fParticleChangeForDecay ;
331} 
332
333void G4Decay::DaughterPolarization(const G4Track& , G4DecayProducts* )
334{
335}
336
337
338
339void G4Decay::StartTracking(G4Track*)
340{
341  currentInteractionLength = -1.0;
342  ResetNumberOfInteractionLengthLeft();
343 
344  fRemainderLifeTime = -1.0;
345}
346
347void G4Decay::EndTracking()
348{
349  // Clear NumberOfInteractionLengthLeft
350  ClearNumberOfInteractionLengthLeft();
351
352  currentInteractionLength = -1.0;
353}
354
355
356G4double G4Decay::PostStepGetPhysicalInteractionLength(
357                             const G4Track& track,
358                             G4double   previousStepSize,
359                             G4ForceCondition* condition
360                            )
361{
362 
363   // condition is set to "Not Forced"
364  *condition = NotForced;
365
366   // pre-assigned Decay time
367  G4double pTime = track.GetDynamicParticle()->GetPreAssignedDecayProperTime();
368  G4double aLife = track.GetDynamicParticle()->GetDefinition()->GetPDGLifeTime();
369
370  if (pTime < 0.) {
371    // normal case
372    if ( previousStepSize > 0.0){
373      // subtract NumberOfInteractionLengthLeft
374      SubtractNumberOfInteractionLengthLeft(previousStepSize);
375      if(theNumberOfInteractionLengthLeft<0.){
376        theNumberOfInteractionLengthLeft=perMillion;
377      }
378      fRemainderLifeTime = theNumberOfInteractionLengthLeft*aLife;
379    }
380    // get mean free path
381    currentInteractionLength = GetMeanFreePath(track, previousStepSize, condition);
382   
383#ifdef G4VERBOSE
384    if ((currentInteractionLength <=0.0) || (verboseLevel>2)){
385      G4cout << "G4Decay::PostStepGetPhysicalInteractionLength " << G4endl;
386      track.GetDynamicParticle()->DumpInfo();
387      G4cout << " in Material  " << track.GetMaterial()->GetName() <<G4endl;
388      G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" <<G4endl;
389    }
390#endif
391
392    G4double value;
393    if (currentInteractionLength <DBL_MAX) {
394      value = theNumberOfInteractionLengthLeft * currentInteractionLength;
395    } else {
396      value = DBL_MAX;
397    }
398
399    return value;
400
401  } else {
402    //pre-assigned Decay time case
403    // reminder proper time
404    fRemainderLifeTime = pTime - track.GetProperTime();
405    if (fRemainderLifeTime <= 0.0) fRemainderLifeTime = DBL_MIN;
406    G4double  rvalue=0.0; 
407    // use pre-assigned Decay time to determine PIL
408    if (aLife>0.0) {
409      // ordinary particle
410      rvalue = (fRemainderLifeTime/aLife)*GetMeanFreePath(track, previousStepSize
411, condition);
412    } else {
413     // shortlived particle
414      rvalue = c_light * fRemainderLifeTime;
415     // by using normalized kinetic energy (= Ekin/mass)
416     G4double   aMass =  track.GetDynamicParticle()->GetMass();
417     rvalue *= track.GetDynamicParticle()->GetTotalMomentum()/aMass;
418    }
419    return rvalue;
420 
421  }
422}
423
424G4double G4Decay::AtRestGetPhysicalInteractionLength(
425                             const G4Track& track,
426                             G4ForceCondition* condition
427                            )
428{
429     // condition is set to "Not Forced"
430  *condition = NotForced;
431
432  G4double pTime = track.GetDynamicParticle()->GetPreAssignedDecayProperTime();
433  if (pTime >= 0.) {
434    fRemainderLifeTime = pTime - track.GetProperTime();
435    if (fRemainderLifeTime <= 0.0) fRemainderLifeTime = DBL_MIN;
436  } else {
437    fRemainderLifeTime = 
438      theNumberOfInteractionLengthLeft * GetMeanLifeTime(track, condition);
439  }
440  return fRemainderLifeTime;
441}
Note: See TracBrowser for help on using the repository browser.