ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4KineticTrack.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4KineticTrack.cc
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 // -----------------------------------------------------------------------------
28 // GEANT 4 class implementation file
29 //
30 // History: first implementation, A. Feliciello, 20th May 1998
31 // -----------------------------------------------------------------------------
32 
33 #include "globals.hh"
34 #include "G4ios.hh"
35 //#include <cmath>
36 
37 #include "Randomize.hh"
38 #include "G4SimpleIntegration.hh"
39 #include "G4ThreeVector.hh"
40 #include "G4LorentzVector.hh"
41 #include "G4KineticTrack.hh"
42 #include "G4KineticTrackVector.hh"
43 #include "G4ParticleDefinition.hh"
44 #include "G4DecayTable.hh"
46 #include "G4DecayProducts.hh"
47 #include "G4LorentzRotation.hh"
48 #include "G4SampleResonance.hh"
49 #include "G4Integrator.hh"
50 #include "G4KaonZero.hh"
51 #include "G4KaonZeroShort.hh"
52 #include "G4KaonZeroLong.hh"
53 #include "G4AntiKaonZero.hh"
54 
55 #include "G4HadTmpUtil.hh"
56 
57 //
58 // Some static clobal for integration
59 //
60 
62 
63 //
64 // Default constructor
65 //
66 
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 }
93 
94 
95 
96 //
97 // Copy constructor
98 //
99 
101 {
102  G4int i;
103  theDefinition = right.GetDefinition();
105  thePosition = right.GetPosition();
109  theNucleon=right.theNucleon;
110  nChannels = right.GetnChannels();
111  theActualMass = right.GetActualMass();
113  for (i = 0; i < nChannels; i++)
114  {
115  theActualWidth[i] = right.theActualWidth[i];
116  }
117  theDaughterMass = 0;
118  theDaughterWidth = 0;
121 
123 // DEBUG //
125 
126 /*
127  G4cerr << G4endl << G4endl << G4endl;
128  G4cerr << " G4KineticTrack copy constructor invoked! \n";
129  G4cerr << " ======================================== \n" <<G4endl;
130 */
131 }
132 
133 
134 //
135 // By argument constructor
136 //
137 
139  G4double aFormationTime,
140  const G4ThreeVector& aPosition,
141  const G4LorentzVector& a4Momentum) :
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 {
154  {
155  if(G4UniformRand()<0.5)
156  {
158  }
159  else
160  {
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 
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 
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,
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 
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
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++;
341  }
342  if ( theDaughterIsShortLived[2] )
343  {
344  nShortLived++;
347  }
348  if ( nShortLived == 0 )
349  {
351  theActualMom = EvaluateCMMomentum(theActualMass,
352  theDaughterMass);
353  thePoleMom = EvaluateCMMomentum(thePoleMass,
355  }
356 // else if ( nShortLived == 1 )
357  else if ( nShortLived >= 1 )
358  {
359  // need the shortlived particle in slot 1! (very bad style...)
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 }
415 
417  const G4ThreeVector& aPosition,
418  const G4LorentzVector& a4Momentum)
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 {
434  Set4Momentum(a4Momentum);
435 }
436 
437 
439 {
440  if (theActualWidth != 0) delete [] theActualWidth;
441  if (theDaughterMass != 0) delete [] theDaughterMass;
442  if (theDaughterWidth != 0) delete [] theDaughterWidth;
443 }
444 
445 
446 
448 {
449  if (this != &right)
450  {
451  theDefinition = right.GetDefinition();
453  the4Momentum = right.the4Momentum;
457  theNucleon=right.theNucleon;
459  if (theActualWidth != 0) delete [] theActualWidth;
460  nChannels = right.GetnChannels();
462  for ( G4int i = 0; i < nChannels; i++)
463  {
464  theActualWidth[i] = right.theActualWidth[i];
465  }
466  }
467  return *this;
468 }
469 
470 
471 
473 {
474  return (this == & right);
475 }
476 
477 
478 
480 {
481  return (this != & right);
482 }
483 
484 
485 
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 definition 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;
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(std::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 }
743 
745 {
746  G4double mass = theActualMass; /* the actual mass value */
747  G4double mass1 = theDaughterMass[0];
748  G4double mass2 = theDaughterMass[1];
749  G4double gamma2 = theDaughterWidth[1];
750 
751  G4double result = (1. / (2 * mass)) *
752  std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
753  ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
754  BrWig(gamma2, mass2, xmass);
755  return result;
756 }
757 
759 {
760  G4double mass = theDefinition->GetPDGMass(); /* the pole mass value */
761  G4double mass1 = theDaughterMass[0];
762  G4double mass2 = theDaughterMass[1];
763  G4double gamma2 = theDaughterWidth[1];
764  G4double result = (1. / (2 * mass)) *
765  std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
766  ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
767  BrWig(gamma2, mass2, xmass);
768  return result;
769 }
770 
772 {
773  const G4double mass = G4KineticTrack_Gmass; /* the actual mass value */
774 // const G4double mass1 = theDaughterMass[0];
775  const G4double mass2 = theDaughterMass[1];
776  const G4double gamma2 = theDaughterWidth[1];
777 
778  const G4double result = (1. / (2 * mass)) *
779  std::sqrt(((mass * mass) - (G4KineticTrack_xmass1 + xmass) * (G4KineticTrack_xmass1 + xmass)) *
780  ((mass * mass) - (G4KineticTrack_xmass1 - xmass) * (G4KineticTrack_xmass1 - xmass))) *
781  BrWig(gamma2, mass2, xmass);
782  return result;
783 }
784 
786 {
788  const G4double mass1 = theDaughterMass[0];
789  const G4double gamma1 = theDaughterWidth[0];
790 // const G4double mass2 = theDaughterMass[1];
791 
792  G4KineticTrack_xmass1 = xmass;
793 
794  const G4double theLowerLimit = 0.0;
795  const G4double theUpperLimit = mass - xmass;
796  const G4int nIterations = 100;
797 
799  G4double result = BrWig(gamma1, mass1, xmass)*
800  integral.Simpson(this, &G4KineticTrack::IntegrandFunction3, theLowerLimit, theUpperLimit, nIterations);
801 
802  return result;
803 }
804 
806 {
807  const G4double theUpperLimit = theActualMass - theDaughterMass[0];
808  const G4int nIterations = 100;
809 
810  if (theLowerLimit>=theUpperLimit) return 0.0;
811 
813  G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction1,
814  theLowerLimit, theUpperLimit, nIterations);
815  return theIntegralOverMass2;
816 }
817 
818 G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit, const G4double poleMass) const
819 {
820  const G4double theUpperLimit = poleMass - theDaughterMass[0];
821  const G4int nIterations = 100;
822 
823  if (theLowerLimit>=theUpperLimit) return 0.0;
824 
826  const G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction2,
827  theLowerLimit, theUpperLimit, nIterations);
828  return theIntegralOverMass2;
829 }
830 
831 
833 {
834  const G4double theLowerLimit = 0.0;
835  const G4double theUpperLimit = theActualMass;
836  const G4int nIterations = 100;
837 
838  if (theLowerLimit>=theUpperLimit) return 0.0;
839 
841  G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction4,
842  theLowerLimit, theUpperLimit, nIterations);
843  return theIntegralOverMass2;
844 }
845 
846 
847 
848 
849 
850 
851 
852