Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4KineticTrack Class Reference

#include <G4KineticTrack.hh>

Inheritance diagram for G4KineticTrack:
Collaboration 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)
 
G4int operator== (const G4KineticTrack &right) const
 
G4int 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
 
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)
 
int operator== (const G4VKineticNucleon &right) const
 
int operator!= (const G4VKineticNucleon &right) const
 

Detailed Description

Definition at line 57 of file G4KineticTrack.hh.

Member Enumeration Documentation

Constructor & Destructor Documentation

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 {
84 // DEBUG //
86 
87 /*
88  G4cerr << G4endl << G4endl << G4endl;
89  G4cerr << " G4KineticTrack default constructor invoked! \n";
90  G4cerr << " =========================================== \n" << G4endl;
91 */
92 }
G4KineticTrack::G4KineticTrack ( const G4KineticTrack right)

Definition at line 100 of file G4KineticTrack.cc.

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

Here is the call graph for this function:

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

Definition at line 138 of file G4KineticTrack.cc.

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

Here is the call graph for this function:

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

Definition at line 416 of file G4KineticTrack.cc.

419  : theDefinition(nucleon->GetDefinition()),
420  theFormationTime(0),
421  thePosition(aPosition),
422  the4Momentum(a4Momentum),
423  theFermi3Momentum(nucleon->GetMomentum()),
424  theNucleon(nucleon),
425  nChannels(0),
426  theActualMass(nucleon->GetDefinition()->GetPDGMass()),
427  theActualWidth(0),
428  theDaughterMass(0),
429  theDaughterWidth(0),
430  theStateToNucleus(undefined),
431  theProjectilePotential(0)
432 {
433  theFermi3Momentum.setE(0);
434  Set4Momentum(a4Momentum);
435 }
const G4LorentzVector & GetMomentum() const
Definition: G4Nucleon.hh:71
virtual const G4ParticleDefinition * GetDefinition() const
Definition: G4Nucleon.hh:85
void Set4Momentum(const G4LorentzVector &a4Momentum)
G4double GetPDGMass() const

Here is the call graph for this function:

G4KineticTrack::~G4KineticTrack ( )

Definition at line 438 of file G4KineticTrack.cc.

439 {
440  if (theActualWidth != 0) delete [] theActualWidth;
441  if (theDaughterMass != 0) delete [] theDaughterMass;
442  if (theDaughterWidth != 0) delete [] theDaughterWidth;
443 }

Member Function Documentation

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

Definition at line 379 of file G4KineticTrack.hh.

380 {
381  G4double Norm = CLHEP::twopi;
382  return (Gamma/((mass-rmass)*(mass-rmass)+Gamma*Gamma/4.))/Norm;
383 }
double G4double
Definition: G4Types.hh:76
static constexpr double twopi
Definition: SystemOfUnits.h:55
G4KineticTrackVector * G4KineticTrack::Decay ( )
virtual

Reimplemented from G4VKineticNucleon.

Definition at line 486 of file G4KineticTrack.cc.

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

Here is the call graph for this function:

Here is the caller graph for this function:

const G4LorentzVector & G4KineticTrack::Get4Momentum ( ) const
inlinevirtual

Implements G4VKineticNucleon.

Definition at line 240 of file G4KineticTrack.hh.

241 {
242  return theTotal4Momentum;
243 }

Here is the caller graph for this function:

G4double G4KineticTrack::GetActualMass ( ) const
inline

Definition at line 319 of file G4KineticTrack.hh.

320 {
321  return std::sqrt(std::abs(the4Momentum.mag2()));
322 }
double mag2() const

Here is the call graph for this function:

Here is the caller graph for this function:

G4double * G4KineticTrack::GetActualWidth ( ) const
inline

Definition at line 334 of file G4KineticTrack.hh.

335 {
336  return theActualWidth;
337 }
const G4ParticleDefinition * G4KineticTrack::GetDefinition ( ) const
inlinevirtual

Implements G4VKineticNucleon.

Definition at line 210 of file G4KineticTrack.hh.

211 {
212  return theDefinition;
213 }

Here is the caller graph for this function:

G4double G4KineticTrack::GetFormationTime ( ) const
inline

Definition at line 220 of file G4KineticTrack.hh.

221 {
222  return theFormationTime;
223 }

Here is the caller graph for this function:

G4int G4KineticTrack::GetnChannels ( ) const
inline

Definition at line 324 of file G4KineticTrack.hh.

325 {
326  return nChannels;
327 }

Here is the caller graph for this function:

const G4ThreeVector & G4KineticTrack::GetPosition ( void  ) const
inlinevirtual

Implements G4VKineticNucleon.

Definition at line 230 of file G4KineticTrack.hh.

231 {
232  return thePosition;
233 }

Here is the caller graph for this function:

G4double G4KineticTrack::GetProjectilePotential ( ) const
inline

Definition at line 421 of file G4KineticTrack.hh.

422 {
423  return theProjectilePotential;
424 }

Here is the caller graph for this function:

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

Definition at line 402 of file G4KineticTrack.hh.

403 {
404  return theStateToNucleus;
405 }

Here is the caller graph for this function:

const G4LorentzVector & G4KineticTrack::GetTrackingMomentum ( ) const
inline

Definition at line 245 of file G4KineticTrack.hh.

246 {
247  return the4Momentum;
248 }

Here is the caller graph for this function:

void G4KineticTrack::Hit ( )
inline

Definition at line 386 of file G4KineticTrack.hh.

387 {
388  if(theNucleon)
389  {
390  theNucleon->Hit(1);
391  }
392 }
void Hit(G4VSplitableHadron *aHit)
Definition: G4Nucleon.hh:90

Here is the call graph for this function:

G4bool G4KineticTrack::IsParticipant ( ) const
inline

Definition at line 395 of file G4KineticTrack.hh.

396 {
397  if(!theNucleon) return true;
398  return theNucleon->AreYouHit();
399 }
G4bool AreYouHit() const
Definition: G4Nucleon.hh:97

Here is the call graph for this function:

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

Definition at line 479 of file G4KineticTrack.cc.

480 {
481  return (this != & right);
482 }
G4KineticTrack & G4KineticTrack::operator= ( const G4KineticTrack right)

Definition at line 447 of file G4KineticTrack.cc.

448 {
449  if (this != &right)
450  {
451  theDefinition = right.GetDefinition();
452  theFormationTime = right.GetFormationTime();
453  the4Momentum = right.the4Momentum;
454  the4Momentum = right.GetTrackingMomentum();
455  theFermi3Momentum = right.theFermi3Momentum;
456  theTotal4Momentum = right.theTotal4Momentum;
457  theNucleon=right.theNucleon;
458  theStateToNucleus=right.theStateToNucleus;
459  if (theActualWidth != 0) delete [] theActualWidth;
460  nChannels = right.GetnChannels();
461  theActualWidth = new G4double[nChannels];
462  for ( G4int i = 0; i < nChannels; i++)
463  {
464  theActualWidth[i] = right.theActualWidth[i];
465  }
466  }
467  return *this;
468 }
G4int GetnChannels() const
int G4int
Definition: G4Types.hh:78
G4double GetFormationTime() const
const G4LorentzVector & GetTrackingMomentum() const
double G4double
Definition: G4Types.hh:76
const G4ParticleDefinition * GetDefinition() const

Here is the call graph for this function:

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

Definition at line 472 of file G4KineticTrack.cc.

473 {
474  return (this == & right);
475 }
G4double G4KineticTrack::SampleResidualLifetime ( )
inline

Definition at line 357 of file G4KineticTrack.hh.

358 {
359  G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
360  G4double tau = CLHEP::hbar_Planck * (-1.0 / theTotalActualWidth);
361  G4double theResidualLifetime = tau * G4Log(G4UniformRand());
362  return theResidualLifetime*the4Momentum.gamma();
363 }
#define G4UniformRand()
Definition: Randomize.hh:97
G4double G4Log(G4double x)
Definition: G4Log.hh:230
double G4double
Definition: G4Types.hh:76
static constexpr double hbar_Planck

Here is the call graph for this function:

Here is the caller graph for this function:

void G4KineticTrack::Set4Momentum ( const G4LorentzVector a4Momentum)
inline

Definition at line 250 of file G4KineticTrack.hh.

251 {
252 // set the4Momentum and update theTotal4Momentum
253 
254  theTotal4Momentum=a4Momentum;
255  the4Momentum = theTotal4Momentum;
256  theFermi3Momentum=G4LorentzVector(0);
257 }
CLHEP::HepLorentzVector G4LorentzVector

Here is the caller graph for this function:

void G4KineticTrack::SetDefinition ( const G4ParticleDefinition aDefinition)
inline

Definition at line 215 of file G4KineticTrack.hh.

216 {
217  theDefinition = aDefinition;
218 }
void G4KineticTrack::SetFormationTime ( G4double  aFormationTime)
inline

Definition at line 225 of file G4KineticTrack.hh.

226 {
227  theFormationTime = aFormationTime;
228 }

Here is the caller graph for this function:

void G4KineticTrack::SetNucleon ( G4Nucleon aN)
inline

Definition at line 106 of file G4KineticTrack.hh.

106 {theNucleon = aN;}
void G4KineticTrack::SetPosition ( const G4ThreeVector  aPosition)
inline

Definition at line 235 of file G4KineticTrack.hh.

236 {
237  thePosition = aPosition;
238 }

Here is the caller graph for this function:

void G4KineticTrack::SetProjectilePotential ( const G4double  aPotential)
inline

Definition at line 416 of file G4KineticTrack.hh.

417 {
418  theProjectilePotential = aPotential;
419 }
G4KineticTrack::CascadeState G4KineticTrack::SetState ( const CascadeState  new_state)
inline

Definition at line 408 of file G4KineticTrack.hh.

409 {
410  CascadeState old_state=theStateToNucleus;
411  theStateToNucleus=new_state;
412  return old_state;
413 }

Here is the caller graph for this function:

void G4KineticTrack::SetTrackingMomentum ( const G4LorentzVector a4Momentum)
inline

Definition at line 283 of file G4KineticTrack.hh.

284 {
285 // set the4Momentum and update theTotal4Momentum, keep the mass of aMomentum
286 
287  the4Momentum = aMomentum;
288  theTotal4Momentum=the4Momentum+theFermi3Momentum;
289 // keep mass of aMomentum for the total momentum
290  G4double mass2 = aMomentum.mag2();
291  G4double p2=theTotal4Momentum.vect().mag2();
292  theTotal4Momentum.setE(std::sqrt(mass2+p2));
293 }
Hep3Vector vect() const
double mag2() const
double mag2() const
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4KineticTrack::Update4Momentum ( G4double  aEnergy)
inline

Definition at line 259 of file G4KineticTrack.hh.

260 {
261 // update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
262 // updates theTotal4Momentum as well.
263  G4double newP(0);
264  G4double mass2=theTotal4Momentum.mag2();
265  if ( sqr(aEnergy) > mass2 )
266  {
267  newP = std::sqrt(sqr(aEnergy) - mass2 );
268  } else
269  {
270  aEnergy=std::sqrt(mass2);
271  }
272  Set4Momentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
273 }
Hep3Vector vect() const
void Set4Momentum(const G4LorentzVector &a4Momentum)
Hep3Vector unit() const
double mag2() const
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
CLHEP::HepLorentzVector G4LorentzVector

Here is the call graph for this function:

void G4KineticTrack::Update4Momentum ( const G4ThreeVector aMomentum)
inline

Definition at line 275 of file G4KineticTrack.hh.

276 {
277 // update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
278 // updates theTotal4Momentum as well.
279  G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
280  Set4Momentum(G4LorentzVector(aMomentum, newE));
281 }
void Set4Momentum(const G4LorentzVector &a4Momentum)
double mag2() const
double mag2() const
double G4double
Definition: G4Types.hh:76
CLHEP::HepLorentzVector G4LorentzVector

Here is the call graph for this function:

void G4KineticTrack::UpdateTrackingMomentum ( G4double  aEnergy)
inline

Definition at line 295 of file G4KineticTrack.hh.

296 {
297 // update the4Momentum with aEnergy at constant mass (the4Momentum.mag()
298 // updates theTotal4Momentum as well.
299  G4double newP(0);
300  G4double mass2=theTotal4Momentum.mag2();
301  if ( sqr(aEnergy) > mass2 )
302  {
303  newP = std::sqrt(sqr(aEnergy) - mass2 );
304  } else
305  {
306  aEnergy=std::sqrt(mass2);
307  }
308  SetTrackingMomentum(G4LorentzVector(newP*the4Momentum.vect().unit(), aEnergy));
309 }
Hep3Vector vect() const
Hep3Vector unit() const
double mag2() const
void SetTrackingMomentum(const G4LorentzVector &a4Momentum)
T sqr(const T &x)
Definition: templates.hh:145
double G4double
Definition: G4Types.hh:76
CLHEP::HepLorentzVector G4LorentzVector

Here is the call graph for this function:

void G4KineticTrack::UpdateTrackingMomentum ( const G4ThreeVector aMomentum)
inline

Definition at line 311 of file G4KineticTrack.hh.

312 {
313 // update the4Momentum with aMomentum at constant mass (the4Momentum.mag()
314 // updates theTotal4Momentum as well.
315  G4double newE=std::sqrt(theTotal4Momentum.mag2() + aMomentum.mag2());
316  SetTrackingMomentum(G4LorentzVector(aMomentum, newE));
317 }
double mag2() const
double mag2() const
void SetTrackingMomentum(const G4LorentzVector &a4Momentum)
double G4double
Definition: G4Types.hh:76
CLHEP::HepLorentzVector G4LorentzVector

Here is the call graph for this function:


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