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

Last change on this file since 961 was 961, checked in by garnier, 15 years ago

update processes
File size: 15.1 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.30 2008/09/19 03:19:53 kurasige Exp $
28// GEANT4 tag $Name: geant4-09-02-ref-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  // set Process Sub Type
67  SetProcessSubType(static_cast<int>(DECAY));
68
69#ifdef G4VERBOSE
70  if (GetVerboseLevel()>1) {
71    G4cout << "G4Decay  constructor " << "  Name:" << processName << G4endl;
72  }
73#endif
74
75  pParticleChange = &fParticleChangeForDecay;
76}
77
78G4Decay::~G4Decay()
79{
80  if (pExtDecayer) {
81    delete pExtDecayer;
82  }
83}
84
85G4bool G4Decay::IsApplicable(const G4ParticleDefinition& aParticleType)
86{
87   // check if the particle is stable?
88   if (aParticleType.GetPDGLifeTime() <0.0) {
89     return false;
90   } else if (aParticleType.GetPDGMass() <= 0.0*MeV) {
91     return false;
92   } else {
93     return true; 
94   }
95}
96
97G4double G4Decay::GetMeanLifeTime(const G4Track& aTrack  ,
98                                  G4ForceCondition*)
99{
100   // returns the mean free path in GEANT4 internal units
101   G4double meanlife;
102
103   // get particle
104   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
105   G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
106   G4double aLife = aParticleDef->GetPDGLifeTime();
107
108   // check if the particle is stable?
109   if (aParticleDef->GetPDGStable()) {
110     meanlife = DBL_MAX;
111   
112   } else {
113     meanlife = aLife;
114   }
115
116#ifdef G4VERBOSE
117   if (GetVerboseLevel()>1) {
118     G4cout << "mean life time: "<< meanlife/ns << "[ns]" << G4endl;
119   }
120#endif
121
122   return  meanlife;
123}
124
125G4double G4Decay::GetMeanFreePath(const G4Track& aTrack,G4double, G4ForceCondition*)
126{
127   // get particle
128   const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
129   G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
130   G4double aMass = aParticle->GetMass();
131   G4double aLife = aParticleDef->GetPDGLifeTime();
132
133
134    // returns the mean free path in GEANT4 internal units
135   G4double pathlength;
136   G4double aCtau = c_light * aLife;
137
138   // check if the particle is stable?
139   if (aParticleDef->GetPDGStable()) {
140     pathlength = DBL_MAX;
141
142   //check if the particle has very short life time ?
143   } else if (aCtau < DBL_MIN) { 
144     pathlength =  DBL_MIN;
145 
146   } else {
147    //calculate the mean free path
148    // by using normalized kinetic energy (= Ekin/mass)
149     G4double   rKineticEnergy = aParticle->GetKineticEnergy()/aMass; 
150     if ( rKineticEnergy > HighestValue) {
151       // gamma >>  1
152       pathlength = ( rKineticEnergy + 1.0)* aCtau;
153     } else if ( rKineticEnergy < DBL_MIN ) {
154       // too slow particle
155#ifdef G4VERBOSE
156       if (GetVerboseLevel()>1) {
157         G4cout << "G4Decay::GetMeanFreePath()   !!particle stops!!";
158         G4cout << aParticleDef->GetParticleName() << G4endl;
159         G4cout << "KineticEnergy:" << aParticle->GetKineticEnergy()/GeV <<"[GeV]";
160       }
161#endif
162       pathlength = DBL_MIN;
163     } else {
164       // beta <1
165       pathlength = (aParticle->GetTotalMomentum())/aMass*aCtau ;
166     }
167   }
168  return  pathlength;
169}
170
171void G4Decay::BuildPhysicsTable(const G4ParticleDefinition&)
172{
173  return;
174}
175
176G4VParticleChange* G4Decay::DecayIt(const G4Track& aTrack, const G4Step& )
177{
178  // The DecayIt() method returns by pointer a particle-change object.
179  // Units are expressed in GEANT4 internal units.
180
181  //   Initialize ParticleChange
182  //     all members of G4VParticleChange are set to equal to
183  //     corresponding member in G4Track
184  fParticleChangeForDecay.Initialize(aTrack);
185
186  // get particle
187  const G4DynamicParticle* aParticle = aTrack.GetDynamicParticle();
188  G4ParticleDefinition* aParticleDef = aParticle->GetDefinition();
189
190  // check if  the particle is stable
191  if (aParticleDef->GetPDGStable()) return &fParticleChangeForDecay ;
192 
193
194  //check if thePreAssignedDecayProducts exists
195  const G4DecayProducts* o_products = (aParticle->GetPreAssignedDecayProducts());
196  G4bool isPreAssigned = (o_products != 0);   
197  G4DecayProducts* products = 0;
198
199  // decay table
200  G4DecayTable   *decaytable = aParticleDef->GetDecayTable();
201 
202  // check if external decayer exists
203  G4bool isExtDecayer = (decaytable == 0) && (pExtDecayer !=0);
204
205  // Error due to NO Decay Table
206  if ( (decaytable == 0) && !isExtDecayer &&!isPreAssigned ){
207    if (GetVerboseLevel()>0) {
208      G4cout <<  "G4Decay::DoIt  : decay table not defined  for ";
209      G4cout << aParticle->GetDefinition()->GetParticleName()<< G4endl;
210    }
211    G4Exception( "G4Decay::DecayIt ",
212                 "No Decay Table",JustWarning, 
213                 "Decay table is not defined");
214
215    fParticleChangeForDecay.SetNumberOfSecondaries(0);
216    // Kill the parent particle
217    fParticleChangeForDecay.ProposeTrackStatus( fStopAndKill ) ;
218    fParticleChangeForDecay.ProposeLocalEnergyDeposit(0.0); 
219   
220    ClearNumberOfInteractionLengthLeft();
221    return &fParticleChangeForDecay ;
222  }
223
224  if (isPreAssigned) {
225    // copy decay products
226    products = new G4DecayProducts(*o_products); 
227  } else if ( isExtDecayer ) {
228    // decay according to external decayer
229    products = pExtDecayer->ImportDecayProducts(aTrack);
230  } else {
231    // decay acoording to decay table
232    // choose a decay channel
233    G4VDecayChannel *decaychannel = decaytable->SelectADecayChannel();
234    if (decaychannel == 0 ){
235      // decay channel not found
236      G4Exception("G4Decay::DoIt  : can not determine decay channel ");
237    } else {
238      // execute DecayIt()
239#ifdef G4VERBOSE
240      G4int temp = decaychannel->GetVerboseLevel();
241      if (GetVerboseLevel()>1) {
242        G4cout << "G4Decay::DoIt  : selected decay channel  addr:" << decaychannel <<G4endl;
243        decaychannel->SetVerboseLevel(GetVerboseLevel());
244      }
245#endif
246      products = decaychannel->DecayIt(aParticle->GetMass());
247#ifdef G4VERBOSE
248      if (GetVerboseLevel()>1) {
249        decaychannel->SetVerboseLevel(temp);
250      }
251#endif
252#ifdef G4VERBOSE
253      if (GetVerboseLevel()>2) {
254        if (! products->IsChecked() ) products->DumpInfo();
255      }
256#endif
257    }
258  }
259 
260  // get parent particle information ...................................
261  G4double   ParentEnergy  = aParticle->GetTotalEnergy();
262  G4double   ParentMass    = aParticle->GetMass();
263  if (ParentEnergy < ParentMass) {
264    ParentEnergy = ParentMass;
265    if (GetVerboseLevel()>0) {
266      G4cout << "G4Decay::DoIt  : Total Energy is less than its mass" << G4endl;
267      G4cout << " Particle: " << aParticle->GetDefinition()->GetParticleName();
268      G4cout << " Energy:"    << ParentEnergy/MeV << "[MeV]";
269      G4cout << " Mass:"    << ParentMass/MeV << "[MeV]";
270      G4cout << G4endl;
271    }
272  }
273
274  G4ThreeVector ParentDirection(aParticle->GetMomentumDirection());
275
276  //boost all decay products to laboratory frame
277  G4double energyDeposit = 0.0;
278  G4double finalGlobalTime = aTrack.GetGlobalTime();
279  if (aTrack.GetTrackStatus() == fStopButAlive ){
280    // AtRest case
281    finalGlobalTime += fRemainderLifeTime;
282    energyDeposit += aParticle->GetKineticEnergy();
283    if (isPreAssigned) products->Boost( ParentEnergy, ParentDirection);
284  } else {
285    // PostStep case
286    if (!isExtDecayer) products->Boost( ParentEnergy, ParentDirection);
287  }
288
289   // set polarization for daughter particles
290   DaughterPolarization(aTrack, products);
291
292
293  //add products in fParticleChangeForDecay
294  G4int numberOfSecondaries = products->entries();
295  fParticleChangeForDecay.SetNumberOfSecondaries(numberOfSecondaries);
296#ifdef G4VERBOSE
297  if (GetVerboseLevel()>1) {
298    G4cout << "G4Decay::DoIt  : Decay vertex :";
299    G4cout << " Time: " << finalGlobalTime/ns << "[ns]";
300    G4cout << " X:" << (aTrack.GetPosition()).x() /cm << "[cm]";
301    G4cout << " Y:" << (aTrack.GetPosition()).y() /cm << "[cm]";
302    G4cout << " Z:" << (aTrack.GetPosition()).z() /cm << "[cm]";
303    G4cout << G4endl;
304    G4cout << "G4Decay::DoIt  : decay products in Lab. Frame" << G4endl;
305    products->DumpInfo();
306  }
307#endif
308  G4int index;
309  G4ThreeVector currentPosition;
310  const G4TouchableHandle thand = aTrack.GetTouchableHandle();
311  for (index=0; index < numberOfSecondaries; index++)
312  {
313     // get current position of the track
314     currentPosition = aTrack.GetPosition();
315     // create a new track object
316     G4Track* secondary = new G4Track( products->PopProducts(),
317                                      finalGlobalTime ,
318                                      currentPosition );
319     // switch on good for tracking flag
320     secondary->SetGoodForTrackingFlag();
321     secondary->SetTouchableHandle(thand);
322     // add the secondary track in the List
323     fParticleChangeForDecay.AddSecondary(secondary);
324  }
325  delete products;
326
327  // Kill the parent particle
328  fParticleChangeForDecay.ProposeTrackStatus( fStopAndKill ) ;
329  fParticleChangeForDecay.ProposeLocalEnergyDeposit(energyDeposit); 
330  fParticleChangeForDecay.ProposeGlobalTime( finalGlobalTime );
331  // Clear NumberOfInteractionLengthLeft
332  ClearNumberOfInteractionLengthLeft();
333
334  return &fParticleChangeForDecay ;
335} 
336
337void G4Decay::DaughterPolarization(const G4Track& , G4DecayProducts* )
338{
339}
340
341
342
343void G4Decay::StartTracking(G4Track*)
344{
345  currentInteractionLength = -1.0;
346  ResetNumberOfInteractionLengthLeft();
347 
348  fRemainderLifeTime = -1.0;
349}
350
351void G4Decay::EndTracking()
352{
353  // Clear NumberOfInteractionLengthLeft
354  ClearNumberOfInteractionLengthLeft();
355
356  currentInteractionLength = -1.0;
357}
358
359
360G4double G4Decay::PostStepGetPhysicalInteractionLength(
361                             const G4Track& track,
362                             G4double   previousStepSize,
363                             G4ForceCondition* condition
364                            )
365{
366 
367   // condition is set to "Not Forced"
368  *condition = NotForced;
369
370   // pre-assigned Decay time
371  G4double pTime = track.GetDynamicParticle()->GetPreAssignedDecayProperTime();
372  G4double aLife = track.GetDynamicParticle()->GetDefinition()->GetPDGLifeTime();
373
374  if (pTime < 0.) {
375    // normal case
376    if ( previousStepSize > 0.0){
377      // subtract NumberOfInteractionLengthLeft
378      SubtractNumberOfInteractionLengthLeft(previousStepSize);
379      if(theNumberOfInteractionLengthLeft<0.){
380        theNumberOfInteractionLengthLeft=perMillion;
381      }
382      fRemainderLifeTime = theNumberOfInteractionLengthLeft*aLife;
383    }
384    // get mean free path
385    currentInteractionLength = GetMeanFreePath(track, previousStepSize, condition);
386   
387#ifdef G4VERBOSE
388    if ((currentInteractionLength <=0.0) || (verboseLevel>2)){
389      G4cout << "G4Decay::PostStepGetPhysicalInteractionLength " << G4endl;
390      track.GetDynamicParticle()->DumpInfo();
391      G4cout << " in Material  " << track.GetMaterial()->GetName() <<G4endl;
392      G4cout << "MeanFreePath = " << currentInteractionLength/cm << "[cm]" <<G4endl;
393    }
394#endif
395
396    G4double value;
397    if (currentInteractionLength <DBL_MAX) {
398      value = theNumberOfInteractionLengthLeft * currentInteractionLength;
399    } else {
400      value = DBL_MAX;
401    }
402
403    return value;
404
405  } else {
406    //pre-assigned Decay time case
407    // reminder proper time
408    fRemainderLifeTime = pTime - track.GetProperTime();
409    if (fRemainderLifeTime <= 0.0) fRemainderLifeTime = DBL_MIN;
410   
411    G4double  rvalue=0.0; 
412    // use pre-assigned Decay time to determine PIL
413    if (aLife>0.0) {
414      // ordinary particle
415      rvalue = (fRemainderLifeTime/aLife)*GetMeanFreePath(track, previousStepSize, condition);
416    } else {
417     // shortlived particle
418      rvalue = c_light * fRemainderLifeTime;
419     // by using normalized kinetic energy (= Ekin/mass)
420     G4double   aMass =  track.GetDynamicParticle()->GetMass();
421     rvalue *= track.GetDynamicParticle()->GetTotalMomentum()/aMass;
422    }
423    return rvalue;
424  }
425}
426
427G4double G4Decay::AtRestGetPhysicalInteractionLength(
428                             const G4Track& track,
429                             G4ForceCondition* condition
430                            )
431{
432     // condition is set to "Not Forced"
433  *condition = NotForced;
434
435  G4double pTime = track.GetDynamicParticle()->GetPreAssignedDecayProperTime();
436  if (pTime >= 0.) {
437    fRemainderLifeTime = pTime - track.GetProperTime();
438    if (fRemainderLifeTime <= 0.0) fRemainderLifeTime = DBL_MIN;
439  } else {
440    fRemainderLifeTime = 
441      theNumberOfInteractionLengthLeft * GetMeanLifeTime(track, condition);
442  }
443  return fRemainderLifeTime;
444}
445
446
447void G4Decay::SetExtDecayer(G4VExtDecayer* val)
448{
449  pExtDecayer = val;
450
451  // set Process Sub Type
452  if ( pExtDecayer !=0 ) {
453    SetProcessSubType(static_cast<int>(DECAY_External));
454  }
455}
Note: See TracBrowser for help on using the repository browser.