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

Last change on this file since 892 was 819, checked in by garnier, 17 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.