Geant4 11.1.1
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4KineticTrack Class Reference

#include <G4KineticTrack.hh>

+ Inheritance diagram for G4KineticTrack:

Public Types

enum  CascadeState {
  undefined , outside , going_in , inside ,
  going_out , gone_out , captured , miss_nucleus
}
 

Public Member Functions

 G4KineticTrack ()
 
 G4KineticTrack (const G4KineticTrack &right)
 
 G4KineticTrack (const G4ParticleDefinition *aDefinition, G4double aFormationTime, const G4ThreeVector &aPosition, const G4LorentzVector &a4Momentum)
 
 G4KineticTrack (G4Nucleon *nucleon, const G4ThreeVector &aPosition, const G4LorentzVector &a4Momentum)
 
 ~G4KineticTrack ()
 
G4KineticTrackoperator= (const G4KineticTrack &right)
 
G4bool operator== (const G4KineticTrack &right) const
 
G4bool operator!= (const G4KineticTrack &right) const
 
const G4ParticleDefinitionGetDefinition () const
 
void SetDefinition (const G4ParticleDefinition *aDefinition)
 
G4double GetFormationTime () const
 
void SetFormationTime (G4double aFormationTime)
 
const G4ThreeVectorGetPosition () const
 
void SetPosition (const G4ThreeVector aPosition)
 
const G4LorentzVectorGet4Momentum () const
 
void Set4Momentum (const G4LorentzVector &a4Momentum)
 
void Update4Momentum (G4double aEnergy)
 
void Update4Momentum (const G4ThreeVector &aMomentum)
 
void SetTrackingMomentum (const G4LorentzVector &a4Momentum)
 
void UpdateTrackingMomentum (G4double aEnergy)
 
void UpdateTrackingMomentum (const G4ThreeVector &aMomentum)
 
const G4LorentzVectorGetTrackingMomentum () const
 
G4double SampleResidualLifetime ()
 
void Hit ()
 
void SetNucleon (G4Nucleon *aN)
 
G4bool IsParticipant () const
 
G4KineticTrackVectorDecay ()
 
G4doubleGetActualWidth () const
 
G4double GetActualMass () const
 
G4int GetnChannels () const
 
CascadeState SetState (const CascadeState new_state)
 
CascadeState GetState () const
 
void SetProjectilePotential (const G4double aPotential)
 
G4double GetProjectilePotential () const
 
void SetCreatorModelID (G4int id)
 
G4int GetCreatorModelID () const
 
const G4ParticleDefinitionGetParentResonanceDef () const
 
void SetParentResonanceDef (const G4ParticleDefinition *parent)
 
G4int GetParentResonanceID () const
 
void SetParentResonanceID (const G4int parentID)
 
G4double BrWig (const G4double Gamma, const G4double rmass, const G4double mass) const
 
- Public Member Functions inherited from G4VKineticNucleon
 G4VKineticNucleon ()
 
 G4VKineticNucleon (const G4VKineticNucleon &right)
 
virtual ~G4VKineticNucleon ()
 
const G4VKineticNucleonoperator= (const G4VKineticNucleon &right)
 
G4bool operator== (const G4VKineticNucleon &right) const
 
G4bool operator!= (const G4VKineticNucleon &right) const
 
virtual G4KineticTrackVectorDecay ()
 
virtual const G4LorentzVectorGet4Momentum () const =0
 
virtual const G4ParticleDefinitionGetDefinition () const =0
 
virtual const G4ThreeVectorGetPosition () const =0
 

Detailed Description

Definition at line 56 of file G4KineticTrack.hh.

Member Enumeration Documentation

◆ CascadeState

Enumerator
undefined 
outside 
going_in 
inside 
going_out 
gone_out 
captured 
miss_nucleus 

Definition at line 118 of file G4KineticTrack.hh.

Constructor & Destructor Documentation

◆ G4KineticTrack() [1/4]

G4KineticTrack::G4KineticTrack ( )

Definition at line 67 of file G4KineticTrack.cc.

67 :
68 theDefinition(0),
69 theFormationTime(0),
70 thePosition(0),
71 the4Momentum(0),
72 theFermi3Momentum(0),
73 theTotal4Momentum(0),
74 theNucleon(0),
75 nChannels(0),
76 theActualMass(0),
77 theActualWidth(0),
78 theDaughterMass(0),
79 theDaughterWidth(0),
80 theStateToNucleus(undefined),
81 theProjectilePotential(0),
82 theCreatorModel(-1),
83 theParentResonanceDef(nullptr),
84 theParentResonanceID(0)
85{
86////////////////
87// DEBUG //
88////////////////
89
90/*
91 G4cerr << G4endl << G4endl << G4endl;
92 G4cerr << " G4KineticTrack default constructor invoked! \n";
93 G4cerr << " =========================================== \n" << G4endl;
94*/
95}

Referenced by Decay().

◆ G4KineticTrack() [2/4]

G4KineticTrack::G4KineticTrack ( const G4KineticTrack right)

Definition at line 103 of file G4KineticTrack.cc.

104{
105 theDefinition = right.GetDefinition();
106 theFormationTime = right.GetFormationTime();
107 thePosition = right.GetPosition();
108 the4Momentum = right.GetTrackingMomentum();
109 theFermi3Momentum = right.theFermi3Momentum;
110 theTotal4Momentum = right.theTotal4Momentum;
111 theNucleon=right.theNucleon;
112 nChannels = right.GetnChannels();
113 theActualMass = right.GetActualMass();
114 theActualWidth = new G4double[nChannels];
115 for (G4int i = 0; i < nChannels; i++)
116 {
117 theActualWidth[i] = right.theActualWidth[i];
118 }
119 theDaughterMass = 0;
120 theDaughterWidth = 0;
121 theStateToNucleus = right.theStateToNucleus;
122 theProjectilePotential = right.theProjectilePotential;
123 theCreatorModel = right.GetCreatorModelID();
124 theParentResonanceDef = right.GetParentResonanceDef();
125 theParentResonanceID = right.GetParentResonanceID();
126
127////////////////
128// DEBUG //
129////////////////
130
131/*
132 G4cerr << G4endl << G4endl << G4endl;
133 G4cerr << " G4KineticTrack copy constructor invoked! \n";
134 G4cerr << " ======================================== \n" <<G4endl;
135*/
136}
double G4double
Definition: G4Types.hh:83
int G4int
Definition: G4Types.hh:85
G4double GetFormationTime() const
G4int GetParentResonanceID() const
G4int GetCreatorModelID() const
G4int GetnChannels() const
const G4ThreeVector & GetPosition() const
const G4ParticleDefinition * GetDefinition() const
const G4LorentzVector & GetTrackingMomentum() const
const G4ParticleDefinition * GetParentResonanceDef() const
G4double GetActualMass() const

◆ G4KineticTrack() [3/4]

G4KineticTrack::G4KineticTrack ( const G4ParticleDefinition aDefinition,
G4double  aFormationTime,
const G4ThreeVector aPosition,
const G4LorentzVector a4Momentum 
)

Definition at line 143 of file G4KineticTrack.cc.

146 :
147 theDefinition(aDefinition),
148 theFormationTime(aFormationTime),
149 thePosition(aPosition),
150 the4Momentum(a4Momentum),
151 theFermi3Momentum(0),
152 theTotal4Momentum(a4Momentum),
153 theNucleon(0),
154 theStateToNucleus(undefined),
155 theProjectilePotential(0),
156 theCreatorModel(-1),
157 theParentResonanceDef(nullptr),
158 theParentResonanceID(0)
159{
160 if(G4KaonZero::KaonZero() == theDefinition ||
161 G4AntiKaonZero::AntiKaonZero() == theDefinition)
162 {
163 if(G4UniformRand()<0.5)
164 {
165 theDefinition = G4KaonZeroShort::KaonZeroShort();
166 }
167 else
168 {
169 theDefinition = G4KaonZeroLong::KaonZeroLong();
170 }
171 }
172
173//
174// Get the number of decay channels
175//
176
177 G4DecayTable* theDecayTable = theDefinition->GetDecayTable();
178 if (theDecayTable != 0)
179 {
180 nChannels = theDecayTable->entries();
181
182 }
183 else
184 {
185 nChannels = 0;
186 }
187
188//
189// Get the actual mass value
190//
191
192 theActualMass = GetActualMass();
193
194//
195// Create an array to Store the actual partial widths
196// of the decay channels
197//
198
199 theDaughterMass = 0;
200 theDaughterWidth = 0;
201 theActualWidth = 0;
202 G4bool * theDaughterIsShortLived = 0;
203
204 if(nChannels!=0) theActualWidth = new G4double[nChannels];
205
206 // cout << " ****CONSTR*** ActualMass ******* " << theActualMass << G4endl;
207 G4int index;
208 for (index = nChannels - 1; index >= 0; --index)
209 {
210 G4VDecayChannel* theChannel = theDecayTable->GetDecayChannel(index);
211 G4int nDaughters = theChannel->GetNumberOfDaughters();
212 G4double theMotherWidth;
213 if (nDaughters == 2 || nDaughters == 3)
214 {
215 G4double thePoleMass = theDefinition->GetPDGMass();
216 theMotherWidth = theDefinition->GetPDGWidth();
217 G4double thePoleWidth = theChannel->GetBR()*theMotherWidth;
218 const G4ParticleDefinition* aDaughter;
219 theDaughterMass = new G4double[nDaughters];
220 theDaughterWidth = new G4double[nDaughters];
221 theDaughterIsShortLived = new G4bool[nDaughters];
222 for (G4int n = 0; n < nDaughters; ++n)
223 {
224 aDaughter = theChannel->GetDaughter(n);
225 theDaughterMass[n] = aDaughter->GetPDGMass();
226 theDaughterWidth[n] = aDaughter->GetPDGWidth();
227 theDaughterIsShortLived[n] = aDaughter->IsShortLived();
228 }
229
230//
231// Check whether both the decay products are stable
232//
233
234 G4double theActualMom = 0.0;
235 G4double thePoleMom = 0.0;
236 G4SampleResonance aSampler;
237 if (nDaughters==2)
238 {
239 if ( !theDaughterIsShortLived[0] && !theDaughterIsShortLived[1] )
240 {
241
242 // G4cout << G4endl << "Both the " << nDaughters <<
243 // " decay products are stable!";
244 // cout << " LB: Both decay products STABLE !" << G4endl;
245 // cout << " parent: " << theChannel->GetParentName() << G4endl;
246 // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl;
247 // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl;
248
249 theActualMom = EvaluateCMMomentum(theActualMass,
250 theDaughterMass);
251 thePoleMom = EvaluateCMMomentum(thePoleMass,
252 theDaughterMass);
253 // cout << G4endl;
254 // cout << " LB: ActualMass/DaughterMass " << theActualMass << " " << theDaughterMass << G4endl;
255 // cout << " LB: ActualMom " << theActualMom << G4endl;
256 // cout << " LB: PoleMom " << thePoleMom << G4endl;
257 // cout << G4endl;
258 }
259 else if ( !theDaughterIsShortLived[0] && theDaughterIsShortLived[1] )
260 {
261
262 // G4cout << G4endl << "Only the first of the " << nDaughters <<" decay products is stable!";
263 // cout << " LB: only the first decay product is STABLE !" << G4endl;
264 // cout << " parent: " << theChannel->GetParentName() << G4endl;
265 // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl;
266 // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl;
267
268// global variable definition
269 G4double lowerLimit = aSampler.GetMinimumMass(theChannel->GetDaughter(1));
270 theActualMom = IntegrateCMMomentum(lowerLimit);
271 thePoleMom = IntegrateCMMomentum(lowerLimit, thePoleMass);
272 // cout << " LB Parent Mass = " << G4KineticTrack_Gmass << G4endl;
273 // cout << " LB Actual Mass = " << theActualMass << G4endl;
274 // cout << " LB Daughter1 Mass = " << G4KineticTrack_Gmass1 << G4endl;
275 // cout << " LB Daughter2 Mass = " << G4KineticTrack_Gmass2 << G4endl;
276 // cout << " The Actual Momentum = " << theActualMom << G4endl;
277 // cout << " The Pole Momentum = " << thePoleMom << G4endl;
278 // cout << G4endl;
279
280 }
281 else if ( theDaughterIsShortLived[0] && !theDaughterIsShortLived[1] )
282 {
283
284 // G4cout << G4endl << "Only the second of the " << nDaughters <<
285 // " decay products is stable!";
286 // cout << " LB: only the second decay product is STABLE !" << G4endl;
287 // cout << " parent: " << theChannel->GetParentName() << G4endl;
288 // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl;
289 // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl;
290
291//
292// Swap the content of the theDaughterMass and theDaughterWidth arrays!!!
293//
294
295 G4SwapObj(theDaughterMass, theDaughterMass + 1);
296 G4SwapObj(theDaughterWidth, theDaughterWidth + 1);
297
298// global variable definition
299 G4double lowerLimit = aSampler.GetMinimumMass(theChannel->GetDaughter(0));
300 theActualMom = IntegrateCMMomentum(lowerLimit);
301 thePoleMom = IntegrateCMMomentum(lowerLimit, thePoleMass);
302 // cout << " LB Parent Mass = " << G4KineticTrack_Gmass << G4endl;
303 // cout << " LB Actual Mass = " << theActualMass << G4endl;
304 // cout << " LB Daughter1 Mass = " << G4KineticTrack_Gmass1 << G4endl;
305 // cout << " LB Daughter2 Mass = " << G4KineticTrack_Gmass2 << G4endl;
306 // cout << " The Actual Momentum = " << theActualMom << G4endl;
307 // cout << " The Pole Momentum = " << thePoleMom << G4endl;
308 // cout << G4endl;
309
310 }
311 else if ( theDaughterIsShortLived[0] && theDaughterIsShortLived[1] )
312 {
313
314// G4cout << G4endl << "Both the " << nDaughters <<
315// " decay products are resonances!";
316 // cout << " LB: both decay products are RESONANCES !" << G4endl;
317 // cout << " parent: " << theChannel->GetParentName() << G4endl;
318 // cout << " particle1: " << theChannel->GetDaughterName(0) << G4endl;
319 // cout << " particle2: " << theChannel->GetDaughterName(1) << G4endl;
320
321// global variable definition
322 G4KineticTrack_Gmass = theActualMass;
323 theActualMom = IntegrateCMMomentum2();
324 G4KineticTrack_Gmass = thePoleMass;
325 thePoleMom = IntegrateCMMomentum2();
326 // cout << " LB Parent Mass = " << G4KineticTrack_Gmass << G4endl;
327 // cout << " LB Daughter1 Mass = " << G4KineticTrack_Gmass1 << G4endl;
328 // cout << " LB Daughter2 Mass = " << G4KineticTrack_Gmass2 << G4endl;
329 // cout << " The Actual Momentum = " << theActualMom << G4endl;
330 // cout << " The Pole Momentum = " << thePoleMom << G4endl;
331 // cout << G4endl;
332
333 }
334 }
335 else // (nDaughter==3)
336 {
337
338 G4int nShortLived = 0;
339 if ( theDaughterIsShortLived[0] )
340 {
341 ++nShortLived;
342 }
343 if ( theDaughterIsShortLived[1] )
344 {
345 ++nShortLived;
346 G4SwapObj(theDaughterMass, theDaughterMass + 1);
347 G4SwapObj(theDaughterWidth, theDaughterWidth + 1);
348 }
349 if ( theDaughterIsShortLived[2] )
350 {
351 ++nShortLived;
352 G4SwapObj(theDaughterMass, theDaughterMass + 2);
353 G4SwapObj(theDaughterWidth, theDaughterWidth + 2);
354 }
355 if ( nShortLived == 0 )
356 {
357 theDaughterMass[1]+=theDaughterMass[2];
358 theActualMom = EvaluateCMMomentum(theActualMass,
359 theDaughterMass);
360 thePoleMom = EvaluateCMMomentum(thePoleMass,
361 theDaughterMass);
362 }
363// else if ( nShortLived == 1 )
364 else if ( nShortLived >= 1 )
365 {
366 // need the shortlived particle in slot 1! (very bad style...)
367 G4SwapObj(theDaughterMass, theDaughterMass + 1);
368 G4SwapObj(theDaughterWidth, theDaughterWidth + 1);
369 theDaughterMass[0] += theDaughterMass[2];
370 theActualMom = IntegrateCMMomentum(0.0);
371 thePoleMom = IntegrateCMMomentum(0.0, thePoleMass);
372 }
373// else
374// {
375// throw G4HadronicException(__FILE__, __LINE__, ("can't handle more than one shortlived in 3 particle output channel");
376// }
377
378 }
379
380 //if(nDaughters<3) theChannel->GetAngularMomentum();
381 G4double theMassRatio = thePoleMass / theActualMass;
382 G4double theMomRatio = theActualMom / thePoleMom;
383 // VI 11.06.2015: for l=0 one not need use pow
384 //G4double l=0;
385 //theActualWidth[index] = thePoleWidth * theMassRatio *
386 // std::pow(theMomRatio, (2 * l + 1)) *
387 // (1.2 / (1+ 0.2*std::pow(theMomRatio, (2 * l))));
388 theActualWidth[index] = thePoleWidth * theMassRatio *
389 theMomRatio;
390 delete [] theDaughterMass;
391 theDaughterMass = 0;
392 delete [] theDaughterWidth;
393 theDaughterWidth = 0;
394 delete [] theDaughterIsShortLived;
395 theDaughterIsShortLived = 0;
396 }
397
398 else // nDaughter = 1 ( e.g. K0 decays 50% to Kshort, 50% Klong
399 {
400 theMotherWidth = theDefinition->GetPDGWidth();
401 theActualWidth[index] = theChannel->GetBR()*theMotherWidth;
402 }
403 }
404
405////////////////
406// DEBUG //
407////////////////
408
409// for (G4int y = nChannels - 1; y >= 0; --y)
410// {
411// G4cout << G4endl << theActualWidth[y];
412// }
413// G4cout << G4endl << G4endl << G4endl;
414
415 /*
416 G4cerr << G4endl << G4endl << G4endl;
417 G4cerr << " G4KineticTrack by argument constructor invoked! \n";
418 G4cerr << " =============================================== \n" << G4endl;
419 */
420
421}
bool G4bool
Definition: G4Types.hh:86
#define G4UniformRand()
Definition: Randomize.hh:52
static G4AntiKaonZero * AntiKaonZero()
G4VDecayChannel * GetDecayChannel(G4int index) const
G4int entries() const
static G4KaonZeroLong * KaonZeroLong()
static G4KaonZeroShort * KaonZeroShort()
static G4KaonZero * KaonZero()
Definition: G4KaonZero.cc:103
G4double GetPDGWidth() const
G4DecayTable * GetDecayTable() const
G4double GetMinimumMass(const G4ParticleDefinition *p) const
G4double GetBR() const
G4int GetNumberOfDaughters() const
G4ParticleDefinition * GetDaughter(G4int anIndex)
void G4SwapObj(T *a, T *b)
Definition: templates.hh:112

◆ G4KineticTrack() [4/4]

G4KineticTrack::G4KineticTrack ( G4Nucleon nucleon,
const G4ThreeVector aPosition,
const G4LorentzVector a4Momentum 
)

Definition at line 423 of file G4KineticTrack.cc.

426 : theDefinition(nucleon->GetDefinition()),
427 theFormationTime(0),
428 thePosition(aPosition),
429 the4Momentum(a4Momentum),
430 theFermi3Momentum(nucleon->GetMomentum()),
431 theNucleon(nucleon),
432 nChannels(0),
433 theActualMass(nucleon->GetDefinition()->GetPDGMass()),
434 theActualWidth(0),
435 theDaughterMass(0),
436 theDaughterWidth(0),
437 theStateToNucleus(undefined),
438 theProjectilePotential(0),
439 theCreatorModel(-1),
440 theParentResonanceDef(nullptr),
441 theParentResonanceID(0)
442{
443 theFermi3Momentum.setE(0);
444 Set4Momentum(a4Momentum);
445}
void Set4Momentum(const G4LorentzVector &a4Momentum)
G4bool nucleon(G4int ityp)

◆ ~G4KineticTrack()

G4KineticTrack::~G4KineticTrack ( )

Definition at line 448 of file G4KineticTrack.cc.

449{
450 if (theActualWidth != 0) delete [] theActualWidth;
451 if (theDaughterMass != 0) delete [] theDaughterMass;
452 if (theDaughterWidth != 0) delete [] theDaughterWidth;
453}

Member Function Documentation

◆ BrWig()

G4double G4KineticTrack::BrWig ( const G4double  Gamma,
const G4double  rmass,
const G4double  mass 
) const
inline

Definition at line 390 of file G4KineticTrack.hh.

391{
392 G4double Norm = CLHEP::twopi;
393 return (Gamma/((mass-rmass)*(mass-rmass)+Gamma*Gamma/4.))/Norm;
394}

◆ Decay()

G4KineticTrackVector * G4KineticTrack::Decay ( )
virtual

Reimplemented from G4VKineticNucleon.

Definition at line 496 of file G4KineticTrack.cc.

497{
498//
499// Select a possible decay channel
500//
501/*
502 G4int index1;
503 for (index1 = nChannels - 1; index1 >= 0; --index1)
504 G4cout << "DECAY Actual Width IND/ActualW " << index1 << " " << theActualWidth[index1] << G4endl;
505 G4cout << "DECAY Actual Mass " << theActualMass << G4endl;
506*/
507 const G4ParticleDefinition* thisDefinition = this->GetDefinition();
508 if(!thisDefinition)
509 {
510 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
511 G4cerr << " track has no particle definition associated."<<G4endl;
512 return 0;
513 }
514 G4DecayTable* theDecayTable = thisDefinition->GetDecayTable();
515 if(!theDecayTable)
516 {
517 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
518 G4cerr << " particle definition has no decay table associated."<<G4endl;
519 G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
520 return 0;
521 }
522
523 G4int chargeBalance = G4lrint(theDefinition->GetPDGCharge() );
524 G4int baryonBalance = G4lrint(theDefinition->GetBaryonNumber() );
525 G4LorentzVector energyMomentumBalance(Get4Momentum());
526 G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
527 if (theTotalActualWidth !=0)
528 {
529
530 //AR-16Aug2016 : Repeat the sampling of the decay channel until is
531 // kinematically above threshold or a max number of attempts is reached
532 G4bool isChannelBelowThreshold = true;
533 const G4int maxNumberOfLoops = 10000;
534 G4int loopCounter = 0;
535
536 G4int chosench;
537 G4String theParentName;
538 G4double theParentMass;
539 G4double theBR;
540 G4int theNumberOfDaughters;
541 G4String theDaughtersName1;
542 G4String theDaughtersName2;
543 G4String theDaughtersName3;
544 G4String theDaughtersName4;
545 G4double masses[4]={0.,0.,0.,0.};
546
547 do {
548
549 G4double theSumActualWidth = 0.0;
550 G4double* theCumActualWidth = new G4double[nChannels]{};
551 for (G4int index = nChannels - 1; index >= 0; --index)
552 {
553 theSumActualWidth += theActualWidth[index];
554 theCumActualWidth[index] = theSumActualWidth;
555 // cout << "DECAY Cum. Width " << index << " " << theCumActualWidth[index] << G4endl;
556 }
557 // cout << "DECAY Total Width " << theSumActualWidth << G4endl;
558 // cout << "DECAY Total Width " << theTotalActualWidth << G4endl;
559 G4double r = theTotalActualWidth * G4UniformRand();
560 G4VDecayChannel* theDecayChannel(0);
561 chosench=-1;
562 for (G4int index = nChannels - 1; index >= 0; --index)
563 {
564 if (r < theCumActualWidth[index])
565 {
566 theDecayChannel = theDecayTable->GetDecayChannel(index);
567 // cout << "DECAY SELECTED CHANNEL" << index << G4endl;
568 chosench=index;
569 break;
570 }
571 }
572
573 delete [] theCumActualWidth;
574
575 if(!theDecayChannel)
576 {
577 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
578 G4cerr << " decay channel has 0x0 channel associated."<<G4endl;
579 G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
580 G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
581 return 0;
582 }
583 theParentName = theDecayChannel->GetParentName();
584 theParentMass = this->GetActualMass();
585 theBR = theActualWidth[chosench];
586 // cout << "**BR*** DECAYNEW " << theBR << G4endl;
587 theNumberOfDaughters = theDecayChannel->GetNumberOfDaughters();
588 theDaughtersName1 = "";
589 theDaughtersName2 = "";
590 theDaughtersName3 = "";
591 theDaughtersName4 = "";
592
593 for (G4int i=0; i < 4; ++i) masses[i]=0.;
594 G4int shortlivedDaughters[4];
595 G4int numberOfShortliveds(0);
596 G4double SumLongLivedMass(0);
597 for (G4int aD=0; aD < theNumberOfDaughters ; ++aD)
598 {
599 const G4ParticleDefinition* aDaughter = theDecayChannel->GetDaughter(aD);
600 masses[aD] = aDaughter->GetPDGMass();
601 if ( aDaughter->IsShortLived() )
602 {
603 shortlivedDaughters[numberOfShortliveds]=aD;
604 ++numberOfShortliveds;
605 } else {
606 SumLongLivedMass += aDaughter->GetPDGMass();
607 }
608
609 }
610 switch (theNumberOfDaughters)
611 {
612 case 0:
613 break;
614 case 1:
615 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
616 theDaughtersName2 = "";
617 theDaughtersName3 = "";
618 theDaughtersName4 = "";
619 break;
620 case 2:
621 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
622 theDaughtersName2 = theDecayChannel->GetDaughterName(1);
623 theDaughtersName3 = "";
624 theDaughtersName4 = "";
625 if ( numberOfShortliveds == 1)
626 { G4SampleResonance aSampler;
627 G4double massmax=theParentMass - SumLongLivedMass;
628 const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
629 masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
630 } else if ( numberOfShortliveds == 2) {
631 // choose masses one after the other, start with randomly choosen
632 G4int zero= (G4UniformRand() > 0.5) ? 0 : 1;
633 G4int one = 1-zero;
634 G4SampleResonance aSampler;
635 G4double massmax=theParentMass - aSampler.GetMinimumMass(theDecayChannel->GetDaughter(shortlivedDaughters[one]));
636 const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[zero]);
637 masses[shortlivedDaughters[zero]]=aSampler.SampleMass(aDaughter,massmax);
638 massmax=theParentMass - masses[shortlivedDaughters[zero]];
639 aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[one]);
640 masses[shortlivedDaughters[one]]=aSampler.SampleMass(aDaughter,massmax);
641 }
642 break;
643 case 3:
644 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
645 theDaughtersName2 = theDecayChannel->GetDaughterName(1);
646 theDaughtersName3 = theDecayChannel->GetDaughterName(2);
647 theDaughtersName4 = "";
648 if ( numberOfShortliveds == 1)
649 { G4SampleResonance aSampler;
650 G4double massmax=theParentMass - SumLongLivedMass;
651 const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
652 masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
653 }
654 break;
655 default:
656 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
657 theDaughtersName2 = theDecayChannel->GetDaughterName(1);
658 theDaughtersName3 = theDecayChannel->GetDaughterName(2);
659 theDaughtersName4 = theDecayChannel->GetDaughterName(3);
660 if ( numberOfShortliveds == 1)
661 { G4SampleResonance aSampler;
662 G4double massmax=theParentMass - SumLongLivedMass;
663 const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
664 masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
665 }
666 if ( theNumberOfDaughters > 4 ) {
668 ed << "More than 4 decay daughters: kept only the first 4" << G4endl;
669 G4Exception( "G4KineticTrack::Decay()", "KINTRK5", JustWarning, ed );
670 }
671 break;
672 }
673
674 //AR-16Aug2016 : Check whether the sum of the masses of the daughters is smaller than the parent mass.
675 // If this is still not the case, but the max number of attempts has been reached,
676 // then the subsequent call thePhaseSpaceDecayChannel.DecayIt() will throw an exception.
677 G4double sumDaughterMasses = 0.0;
678 for (G4int i=0; i < 4; ++i) sumDaughterMasses += masses[i];
679 if ( theParentMass - sumDaughterMasses > 0.0 ) isChannelBelowThreshold = false;
680
681 } while ( isChannelBelowThreshold && ++loopCounter < maxNumberOfLoops ); /* Loop checking, 16.08.2016, A.Ribon */
682
683//
684// Get the decay products List
685//
686
687 G4GeneralPhaseSpaceDecay thePhaseSpaceDecayChannel(theParentName,
688 theParentMass,
689 theBR,
690 theNumberOfDaughters,
691 theDaughtersName1,
692 theDaughtersName2,
693 theDaughtersName3,
694 theDaughtersName4,
695 masses);
696 G4DecayProducts* theDecayProducts = thePhaseSpaceDecayChannel.DecayIt();
697 if(!theDecayProducts)
698 {
700 ed << "Error condition encountered: phase-space decay failed." << G4endl
701 << "\t the decaying particle is: " << thisDefinition->GetParticleName() << G4endl
702 << "\t the channel index is: "<< chosench << " of "<< nChannels << "channels" << G4endl
703 << "\t " << theNumberOfDaughters << " daughter particles: "
704 << theDaughtersName1 << " " << theDaughtersName2 << " " << theDaughtersName3 << " "
705 << theDaughtersName4 << G4endl;
706 G4Exception( "G4KineticTrack::Decay ", "HAD_KINTRACK_001", JustWarning, ed );
707 return 0;
708 }
709
710//
711// Create the kinetic track List associated to the decay products
712//
713// For the decay products of hadronic resonances, we assign as creator model ID
714// the same as their parent
715 G4LorentzRotation toMoving(Get4Momentum().boostVector());
716 G4DynamicParticle* theDynamicParticle;
717 G4double formationTime = 0.0;
719 G4LorentzVector momentum;
720 G4LorentzVector momentumBalanceCMS(0);
721 G4KineticTrackVector* theDecayProductList = new G4KineticTrackVector;
722 G4int dEntries = theDecayProducts->entries();
723 const G4ParticleDefinition * aProduct = 0;
724 // Use the integer round mass in keV to get an unique ID for the parent resonance
725 G4int uniqueID = static_cast< G4int >( round( Get4Momentum().mag() / CLHEP::keV ) );
726 for (G4int i=dEntries; i > 0; --i)
727 {
728 theDynamicParticle = theDecayProducts->PopProducts();
729 aProduct = theDynamicParticle->GetDefinition();
730 chargeBalance -= G4lrint(aProduct->GetPDGCharge() );
731 baryonBalance -= G4lrint(aProduct->GetBaryonNumber() );
732 momentumBalanceCMS += theDynamicParticle->Get4Momentum();
733 momentum = toMoving*theDynamicParticle->Get4Momentum();
734 energyMomentumBalance -= momentum;
735 G4KineticTrack* aDaughter = new G4KineticTrack (aProduct,
736 formationTime,
737 position,
738 momentum);
739 if (aDaughter != nullptr)
740 {
743 aDaughter->SetParentResonanceID(uniqueID);
744 }
745 theDecayProductList->push_back(aDaughter);
746 delete theDynamicParticle;
747 }
748 delete theDecayProducts;
749 if(std::getenv("DecayEnergyBalanceCheck"))
750 std::cout << "DEBUGGING energy balance in cms and lab, charge baryon balance : "
751 << momentumBalanceCMS << " "
752 <<energyMomentumBalance << " "
753 <<chargeBalance<<" "
754 <<baryonBalance<<" "
755 <<G4endl;
756 return theDecayProductList;
757 }
758 else
759 {
760 return 0;
761 }
762}
@ JustWarning
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *description)
Definition: G4Exception.cc:59
std::ostringstream G4ExceptionDescription
Definition: G4Exception.hh:40
G4GLOB_DLL std::ostream G4cerr
#define G4endl
Definition: G4ios.hh:57
G4int entries() const
G4DynamicParticle * PopProducts()
G4ParticleDefinition * GetDefinition() const
G4LorentzVector Get4Momentum() const
void SetCreatorModelID(G4int id)
void SetParentResonanceID(const G4int parentID)
const G4LorentzVector & Get4Momentum() const
void SetParentResonanceDef(const G4ParticleDefinition *parent)
G4double GetPDGCharge() const
const G4String & GetParticleName() const
G4double SampleMass(const G4double poleMass, const G4double gamma, const G4double minMass, const G4double maxMass) const
int G4lrint(double ad)
Definition: templates.hh:134

Referenced by G4LMsdGenerator::ApplyYourself(), G4QMDCollision::CalKinematicsOfBinaryCollisions(), G4DecayKineticTracks::Decay(), G4NeutrinoNucleusModel::FinalBarion(), G4NeutrinoNucleusModel::FinalMeson(), and G4BCDecay::GetFinalState().

◆ Get4Momentum()

const G4LorentzVector & G4KineticTrack::Get4Momentum ( ) const
inlinevirtual

◆ GetActualMass()

G4double G4KineticTrack::GetActualMass ( ) const
inline

◆ GetActualWidth()

G4double * G4KineticTrack::GetActualWidth ( ) const
inline

Definition at line 345 of file G4KineticTrack.hh.

346{
347 return theActualWidth;
348}

◆ GetCreatorModelID()

G4int G4KineticTrack::GetCreatorModelID ( ) const
inline

Definition at line 443 of file G4KineticTrack.hh.

444{
445 return theCreatorModel;
446}

Referenced by Decay(), G4DecayKineticTracks::Decay(), G4KineticTrack(), operator=(), and G4DecayStrongResonances::Propagate().

◆ GetDefinition()

const G4ParticleDefinition * G4KineticTrack::GetDefinition ( ) const
inlinevirtual

Implements G4VKineticNucleon.

Definition at line 221 of file G4KineticTrack.hh.

222{
223 return theDefinition;
224}

Referenced by G4CollisionManager::AddCollision(), G4KineticTrackVector::BoostBeam(), G4CollisionComposite::CrossSection(), G4CollisionNN::CrossSection(), G4XAnnihilationChannel::CrossSection(), G4XAqmTotal::CrossSection(), G4XnpElasticLowE::CrossSection(), G4XnpTotalLowE::CrossSection(), G4XPDGElastic::CrossSection(), G4XPDGTotal::CrossSection(), G4XResonance::CrossSection(), Decay(), G4DecayKineticTracks::Decay(), G4VXResonance::DegeneracyFactor(), G4VXResonance::DetailedBalance(), G4VElasticCollision::FinalState(), G4VScatteringCollision::FinalState(), G4Absorber::FindAbsorbers(), G4VCrossSectionSource::FindKeyParticle(), G4VCrossSectionSource::FindLightParticle(), G4Absorber::FindProducts(), G4QGSMFragmentation::FragmentString(), G4ExcitedStringDecay::FragmentStrings(), G4KineticTrack(), G4BCDecay::GetCollisions(), G4MesonAbsorption::GetFinalState(), G4ParticleTypeConverter::GetGenericType(), G4ConcreteMesonBaryonToResonance::GetOutgoingParticle(), G4Scatterer::GetTimeToInteraction(), G4CollisionMesonBaryonElastic::IsInCharge(), G4CollisionNNElastic::IsInCharge(), G4CollisionnpElastic::IsInCharge(), G4ConcreteNNTwoBodyResonance::IsInCharge(), G4GeneralNNCollision::IsInCharge(), G4VXResonance::IsospinCorrection(), G4GeneratorPrecompoundInterface::MakeCoalescence(), operator=(), G4CollisionManager::Print(), G4CollisionInitialState::Print(), G4IntraNucleiCascader::processSecondary(), G4DecayStrongResonances::Propagate(), G4CascadeInterface::Propagate(), G4IntraNucleiCascader::releaseSecondary(), G4Scatterer::Scatter(), G4RKPropagation::Transport(), and G4Absorber::WillBeAbsorbed().

◆ GetFormationTime()

◆ GetnChannels()

G4int G4KineticTrack::GetnChannels ( ) const
inline

Definition at line 335 of file G4KineticTrack.hh.

336{
337 return nChannels;
338}

Referenced by G4KineticTrack(), and operator=().

◆ GetParentResonanceDef()

const G4ParticleDefinition * G4KineticTrack::GetParentResonanceDef ( ) const
inline

Definition at line 449 of file G4KineticTrack.hh.

450{
451 return theParentResonanceDef;
452}

Referenced by G4KineticTrack(), operator=(), and G4DecayStrongResonances::Propagate().

◆ GetParentResonanceID()

G4int G4KineticTrack::GetParentResonanceID ( ) const
inline

Definition at line 461 of file G4KineticTrack.hh.

462{
463 return theParentResonanceID;
464}

Referenced by G4KineticTrack(), operator=(), and G4DecayStrongResonances::Propagate().

◆ GetPosition()

◆ GetProjectilePotential()

G4double G4KineticTrack::GetProjectilePotential ( ) const
inline

Definition at line 432 of file G4KineticTrack.hh.

433{
434 return theProjectilePotential;
435}

Referenced by G4RKPropagation::Transport().

◆ GetState()

G4KineticTrack::CascadeState G4KineticTrack::GetState ( ) const
inline

Definition at line 413 of file G4KineticTrack.hh.

414{
415 return theStateToNucleus;
416}

Referenced by G4RKPropagation::Transport().

◆ GetTrackingMomentum()

const G4LorentzVector & G4KineticTrack::GetTrackingMomentum ( ) const
inline

◆ Hit()

void G4KineticTrack::Hit ( )
inline

Definition at line 397 of file G4KineticTrack.hh.

398{
399 if(theNucleon)
400 {
401 theNucleon->Hit(1);
402 }
403}
void Hit(G4VSplitableHadron *aHit)
Definition: G4Nucleon.hh:91

◆ IsParticipant()

G4bool G4KineticTrack::IsParticipant ( ) const
inline

Definition at line 406 of file G4KineticTrack.hh.

407{
408 if(!theNucleon) return true;
409 return theNucleon->AreYouHit();
410}
G4bool AreYouHit() const
Definition: G4Nucleon.hh:98

◆ operator!=()

G4bool G4KineticTrack::operator!= ( const G4KineticTrack right) const

Definition at line 489 of file G4KineticTrack.cc.

490{
491 return (this != & right);
492}

◆ operator=()

G4KineticTrack & G4KineticTrack::operator= ( const G4KineticTrack right)

Definition at line 457 of file G4KineticTrack.cc.

458{
459 if (this != &right)
460 {
461 theDefinition = right.GetDefinition();
462 theFormationTime = right.GetFormationTime();
463 the4Momentum = right.the4Momentum;
464 the4Momentum = right.GetTrackingMomentum();
465 theFermi3Momentum = right.theFermi3Momentum;
466 theTotal4Momentum = right.theTotal4Momentum;
467 theNucleon = right.theNucleon;
468 theStateToNucleus = right.theStateToNucleus;
469 if (theActualWidth != 0) delete [] theActualWidth;
470 nChannels = right.GetnChannels();
471 theActualWidth = new G4double[nChannels];
472 for (G4int i = 0; i < nChannels; ++i) theActualWidth[i] = right.theActualWidth[i];
473 theCreatorModel = right.GetCreatorModelID();
474 theParentResonanceDef = right.GetParentResonanceDef();
475 theParentResonanceID = right.GetParentResonanceID();
476 }
477 return *this;
478}

◆ operator==()

G4bool G4KineticTrack::operator== ( const G4KineticTrack right) const

Definition at line 482 of file G4KineticTrack.cc.

483{
484 return (this == & right);
485}

◆ SampleResidualLifetime()

G4double G4KineticTrack::SampleResidualLifetime ( )
inline

Definition at line 368 of file G4KineticTrack.hh.

369{
370 G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
371 G4double tau = CLHEP::hbar_Planck * (-1.0 / theTotalActualWidth);
372 G4double theResidualLifetime = tau * G4Log(G4UniformRand());
373 return theResidualLifetime*the4Momentum.gamma();
374}
G4double G4Log(G4double x)
Definition: G4Log.hh:227

Referenced by G4BCDecay::GetCollisions().

◆ Set4Momentum()

void G4KineticTrack::Set4Momentum ( const G4LorentzVector a4Momentum)
inline

Definition at line 261 of file G4KineticTrack.hh.

262{
263// set the4Momentum and update theTotal4Momentum
264
265 theTotal4Momentum=a4Momentum;
266 the4Momentum = theTotal4Momentum;
267 theFermi3Momentum=G4LorentzVector(0);
268}
CLHEP::HepLorentzVector G4LorentzVector

Referenced by G4KineticTrackVector::Boost(), G4KineticTrackVector::BoostBeam(), G4CollisionNN::CrossSection(), G4VElasticCollision::FinalState(), G4QGSMFragmentation::FragmentString(), G4KineticTrack(), G4ExcitedString::LorentzRotate(), G4ExcitedString::TransformToCenterOfMass(), and Update4Momentum().

◆ SetCreatorModelID()

void G4KineticTrack::SetCreatorModelID ( G4int  id)
inline

Definition at line 438 of file G4KineticTrack.hh.

439{
440 theCreatorModel = id;
441}

Referenced by Decay(), G4DecayKineticTracks::Decay(), G4GeneratorPrecompoundInterface::MakeCoalescence(), and G4QuasiElasticChannel::Scatter().

◆ SetDefinition()

void G4KineticTrack::SetDefinition ( const G4ParticleDefinition aDefinition)
inline

Definition at line 226 of file G4KineticTrack.hh.

227{
228 theDefinition = aDefinition;
229}

◆ SetFormationTime()

void G4KineticTrack::SetFormationTime ( G4double  aFormationTime)
inline

Definition at line 236 of file G4KineticTrack.hh.

237{
238 theFormationTime = aFormationTime;
239}

Referenced by G4QGSMFragmentation::FragmentString().

◆ SetNucleon()

void G4KineticTrack::SetNucleon ( G4Nucleon aN)
inline

Definition at line 105 of file G4KineticTrack.hh.

105{theNucleon = aN;}

◆ SetParentResonanceDef()

void G4KineticTrack::SetParentResonanceDef ( const G4ParticleDefinition parent)
inline

Definition at line 455 of file G4KineticTrack.hh.

456{
457 theParentResonanceDef = parentDef;
458}

Referenced by Decay(), and G4DecayKineticTracks::Decay().

◆ SetParentResonanceID()

void G4KineticTrack::SetParentResonanceID ( const G4int  parentID)
inline

Definition at line 467 of file G4KineticTrack.hh.

468{
469 theParentResonanceID = parentID;
470}

Referenced by Decay(), and G4DecayKineticTracks::Decay().

◆ SetPosition()

void G4KineticTrack::SetPosition ( const G4ThreeVector  aPosition)
inline

Definition at line 246 of file G4KineticTrack.hh.

247{
248 thePosition = aPosition;
249}

Referenced by G4QGSMFragmentation::FragmentString(), G4ExcitedStringDecay::FragmentStrings(), and G4KineticTrackVector::Shift().

◆ SetProjectilePotential()

void G4KineticTrack::SetProjectilePotential ( const G4double  aPotential)
inline

Definition at line 427 of file G4KineticTrack.hh.

428{
429 theProjectilePotential = aPotential;
430}

◆ SetState()

G4KineticTrack::CascadeState G4KineticTrack::SetState ( const CascadeState  new_state)
inline

Definition at line 419 of file G4KineticTrack.hh.

420{
421 CascadeState old_state=theStateToNucleus;
422 theStateToNucleus=new_state;
423 return old_state;
424}

Referenced by G4BinaryCascade::ApplyYourself(), and G4RKPropagation::Transport().

◆ SetTrackingMomentum()

void G4KineticTrack::SetTrackingMomentum ( const G4LorentzVector a4Momentum)
inline

Definition at line 294 of file G4KineticTrack.hh.

295{
296// set the4Momentum and update theTotal4Momentum, keep the mass of aMomentum
297
298 the4Momentum = aMomentum;
299 theTotal4Momentum=the4Momentum+theFermi3Momentum;
300// keep mass of aMomentum for the total momentum
301 G4double mass2 = aMomentum.mag2();
302 G4double p2=theTotal4Momentum.vect().mag2();
303 theTotal4Momentum.setE(std::sqrt(mass2+p2));
304}
double mag2() const
Hep3Vector vect() const

Referenced by G4RKPropagation::Transport(), and UpdateTrackingMomentum().

◆ Update4Momentum() [1/2]

void G4KineticTrack::Update4Momentum ( const G4ThreeVector aMomentum)
inline

Definition at line 286 of file G4KineticTrack.hh.

287{
288// update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
289// updates theTotal4Momentum as well.
290 G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
291 Set4Momentum(G4LorentzVector(aMomentum, newE));
292}

◆ Update4Momentum() [2/2]

void G4KineticTrack::Update4Momentum ( G4double  aEnergy)
inline

Definition at line 270 of file G4KineticTrack.hh.

271{
272// update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
273// updates theTotal4Momentum as well.
274 G4double newP(0);
275 G4double mass2=theTotal4Momentum.mag2();
276 if ( sqr(aEnergy) > mass2 )
277 {
278 newP = std::sqrt(sqr(aEnergy) - mass2 );
279 } else
280 {
281 aEnergy=std::sqrt(mass2);
282 }
283 Set4Momentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
284}
Hep3Vector unit() const
T sqr(const T &x)
Definition: templates.hh:128

◆ UpdateTrackingMomentum() [1/2]

void G4KineticTrack::UpdateTrackingMomentum ( const G4ThreeVector aMomentum)
inline

Definition at line 322 of file G4KineticTrack.hh.

323{
324// update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
325// updates theTotal4Momentum as well.
326 G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
327 SetTrackingMomentum(G4LorentzVector(aMomentum, newE));
328}
void SetTrackingMomentum(const G4LorentzVector &a4Momentum)

◆ UpdateTrackingMomentum() [2/2]

void G4KineticTrack::UpdateTrackingMomentum ( G4double  aEnergy)
inline

Definition at line 306 of file G4KineticTrack.hh.

307{
308// update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
309// updates theTotal4Momentum as well.
310 G4double newP(0);
311 G4double mass2=theTotal4Momentum.mag2();
312 if ( sqr(aEnergy) > mass2 )
313 {
314 newP = std::sqrt(sqr(aEnergy) - mass2 );
315 } else
316 {
317 aEnergy=std::sqrt(mass2);
318 }
319 SetTrackingMomentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
320}

The documentation for this class was generated from the following files: