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

Last change on this file since 1016 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 15.0 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 $
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.