Geant4  10.02.p01
G4KineticTrack.cc
Go to the documentation of this file.
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 {
433  theFermi3Momentum.setE(0);
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 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  G4int index;
520  G4double theSumActualWidth = 0.0;
521  G4double* theCumActualWidth = new G4double[nChannels];
522  for (index = nChannels - 1; index >= 0; index--)
523  {
524  theSumActualWidth += theActualWidth[index];
525  theCumActualWidth[index] = theSumActualWidth;
526  // cout << "DECAY Cum. Width " << index << " " << theCumActualWidth[index] << G4endl;
527  }
528  // cout << "DECAY Total Width " << theSumActualWidth << G4endl;
529  // cout << "DECAY Total Width " << theTotalActualWidth << G4endl;
530  G4double r = theTotalActualWidth * G4UniformRand();
531  G4VDecayChannel* theDecayChannel(0);
532  G4int chosench=-1;
533  for (index = nChannels - 1; index >= 0; index--)
534  {
535  if (r < theCumActualWidth[index])
536  {
537  theDecayChannel = theDecayTable->GetDecayChannel(index);
538  // cout << "DECAY SELECTED CHANNEL" << index << G4endl;
539  chosench=index;
540  break;
541  }
542  }
543 
544  delete [] theCumActualWidth;
545 
546  if(!theDecayChannel)
547  {
548  G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
549  G4cerr << " decay channel has 0x0 channel associated."<<G4endl;
550  G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
551  G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
552  return 0;
553  }
554  G4String theParentName = theDecayChannel->GetParentName();
555  G4double theParentMass = this->GetActualMass();
556  G4double theBR = theActualWidth[index];
557  // cout << "**BR*** DECAYNEW " << theBR << G4endl;
558  G4int theNumberOfDaughters = theDecayChannel->GetNumberOfDaughters();
559  G4String theDaughtersName1 = "";
560  G4String theDaughtersName2 = "";
561  G4String theDaughtersName3 = "";
562  G4String theDaughtersName4 = "";
563 
564  G4double masses[4]={0.,0.,0.,0.};
565  G4int shortlivedDaughters[4];
566  G4int numberOfShortliveds(0);
567  G4double SumLongLivedMass(0);
568  for (G4int aD=0; aD < theNumberOfDaughters ; aD++)
569  {
570  const G4ParticleDefinition* aDaughter = theDecayChannel->GetDaughter(aD);
571  masses[aD] = aDaughter->GetPDGMass();
572  if ( aDaughter->IsShortLived() )
573  {
574  shortlivedDaughters[numberOfShortliveds]=aD;
575  numberOfShortliveds++;
576  } else {
577  SumLongLivedMass += aDaughter->GetPDGMass();
578  }
579 
580  }
581  switch (theNumberOfDaughters)
582  {
583  case 0:
584  break;
585  case 1:
586  theDaughtersName1 = theDecayChannel->GetDaughterName(0);
587  theDaughtersName2 = "";
588  theDaughtersName3 = "";
589  theDaughtersName4 = "";
590  break;
591  case 2:
592  theDaughtersName1 = theDecayChannel->GetDaughterName(0);
593  theDaughtersName2 = theDecayChannel->GetDaughterName(1);
594  theDaughtersName3 = "";
595  theDaughtersName4 = "";
596  if ( numberOfShortliveds == 1)
597  { G4SampleResonance aSampler;
598  G4double massmax=theParentMass - SumLongLivedMass;
599  const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
600  masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
601  } else if ( numberOfShortliveds == 2) {
602  // choose masses one after the other, start with randomly choosen
603  G4int zero= (G4UniformRand() > 0.5) ? 0 : 1;
604  G4int one = 1-zero;
605  G4SampleResonance aSampler;
606  G4double massmax=theParentMass - aSampler.GetMinimumMass(theDecayChannel->GetDaughter(shortlivedDaughters[one]));
607  const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[zero]);
608  masses[shortlivedDaughters[zero]]=aSampler.SampleMass(aDaughter,massmax);
609  massmax=theParentMass - masses[shortlivedDaughters[zero]];
610  aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[one]);
611  masses[shortlivedDaughters[one]]=aSampler.SampleMass(aDaughter,massmax);
612  }
613  break;
614  case 3:
615  theDaughtersName1 = theDecayChannel->GetDaughterName(0);
616  theDaughtersName2 = theDecayChannel->GetDaughterName(1);
617  theDaughtersName3 = theDecayChannel->GetDaughterName(2);
618  theDaughtersName4 = "";
619  if ( numberOfShortliveds == 1)
620  { G4SampleResonance aSampler;
621  G4double massmax=theParentMass - SumLongLivedMass;
622  const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
623  masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
624  }
625  break;
626  default:
627  theDaughtersName1 = theDecayChannel->GetDaughterName(0);
628  theDaughtersName2 = theDecayChannel->GetDaughterName(1);
629  theDaughtersName3 = theDecayChannel->GetDaughterName(2);
630  theDaughtersName4 = theDecayChannel->GetDaughterName(3);
631  if ( numberOfShortliveds == 1)
632  { G4SampleResonance aSampler;
633  G4double massmax=theParentMass - SumLongLivedMass;
634  const G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
635  masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
636  }
637  if ( theNumberOfDaughters > 4 ) {
639  ed << "More than 4 decay daughters: kept only the first 4" << G4endl;
640  G4Exception( "G4KineticTrack::Decay()", "KINTRK5", JustWarning, ed );
641  }
642  break;
643  }
644 
645 //
646 // Get the decay products List
647 //
648 
649  G4GeneralPhaseSpaceDecay thePhaseSpaceDecayChannel(theParentName,
650  theParentMass,
651  theBR,
652  theNumberOfDaughters,
653  theDaughtersName1,
654  theDaughtersName2,
655  theDaughtersName3,
656  theDaughtersName4,
657  masses);
658  G4DecayProducts* theDecayProducts = thePhaseSpaceDecayChannel.DecayIt();
659  if(!theDecayProducts)
660  {
662  ed << "Error condition encountered: phase-space decay failed." << G4endl
663  << "\t the decaying particle is: " << thisDefinition->GetParticleName() << G4endl
664  << "\t the channel index is: "<< chosench << " of "<< nChannels << "channels" << G4endl
665  << "\t " << theNumberOfDaughters << " daughter particles: "
666  << theDaughtersName1 << " " << theDaughtersName2 << " " << theDaughtersName3 << " "
667  << theDaughtersName4 << G4endl;
668  G4Exception( "G4KineticTrack::Decay ", "HAD_KINTRACK_001", JustWarning, ed );
669  return 0;
670  }
671 
672 //
673 // Create the kinetic track List associated to the decay products
674 //
675  G4LorentzRotation toMoving(Get4Momentum().boostVector());
676  G4DynamicParticle* theDynamicParticle;
677  G4double formationTime = 0.0;
679  G4LorentzVector momentum;
680  G4LorentzVector momentumBalanceCMS(0);
681  G4KineticTrackVector* theDecayProductList = new G4KineticTrackVector;
682  G4int dEntries = theDecayProducts->entries();
683  const G4ParticleDefinition * aProduct = 0;
684  for (G4int i=dEntries; i > 0; i--)
685  {
686  theDynamicParticle = theDecayProducts->PopProducts();
687  aProduct = theDynamicParticle->GetDefinition();
688  chargeBalance -= G4lrint(aProduct->GetPDGCharge() );
689  baryonBalance -= G4lrint(aProduct->GetBaryonNumber() );
690  momentumBalanceCMS += theDynamicParticle->Get4Momentum();
691  momentum = toMoving*theDynamicParticle->Get4Momentum();
692  energyMomentumBalance -= momentum;
693  theDecayProductList->push_back(new G4KineticTrack (aProduct,
694  formationTime,
695  position,
696  momentum));
697  delete theDynamicParticle;
698  }
699  delete theDecayProducts;
700  if(getenv("DecayEnergyBalanceCheck"))
701  std::cout << "DEBUGGING energy balance in cms and lab, charge baryon balance : "
702  << momentumBalanceCMS << " "
703  <<energyMomentumBalance << " "
704  <<chargeBalance<<" "
705  <<baryonBalance<<" "
706  <<G4endl;
707  return theDecayProductList;
708  }
709  else
710  {
711  return 0;
712  }
713 }
714 
716 {
717  G4double mass = theActualMass; /* the actual mass value */
718  G4double mass1 = theDaughterMass[0];
719  G4double mass2 = theDaughterMass[1];
720  G4double gamma2 = theDaughterWidth[1];
721 
722  G4double result = (1. / (2 * mass)) *
723  std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
724  ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
725  BrWig(gamma2, mass2, xmass);
726  return result;
727 }
728 
730 {
731  G4double mass = theDefinition->GetPDGMass(); /* the pole mass value */
732  G4double mass1 = theDaughterMass[0];
733  G4double mass2 = theDaughterMass[1];
734  G4double gamma2 = theDaughterWidth[1];
735  G4double result = (1. / (2 * mass)) *
736  std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
737  ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
738  BrWig(gamma2, mass2, xmass);
739  return result;
740 }
741 
743 {
744  const G4double mass = G4KineticTrack_Gmass; /* the actual mass value */
745 // const G4double mass1 = theDaughterMass[0];
746  const G4double mass2 = theDaughterMass[1];
747  const G4double gamma2 = theDaughterWidth[1];
748 
749  const G4double result = (1. / (2 * mass)) *
750  std::sqrt(((mass * mass) - (G4KineticTrack_xmass1 + xmass) * (G4KineticTrack_xmass1 + xmass)) *
751  ((mass * mass) - (G4KineticTrack_xmass1 - xmass) * (G4KineticTrack_xmass1 - xmass))) *
752  BrWig(gamma2, mass2, xmass);
753  return result;
754 }
755 
757 {
758  const G4double mass = G4KineticTrack_Gmass;
759  const G4double mass1 = theDaughterMass[0];
760  const G4double gamma1 = theDaughterWidth[0];
761 // const G4double mass2 = theDaughterMass[1];
762 
763  G4KineticTrack_xmass1 = xmass;
764 
765  const G4double theLowerLimit = 0.0;
766  const G4double theUpperLimit = mass - xmass;
767  const G4int nIterations = 100;
768 
770  G4double result = BrWig(gamma1, mass1, xmass)*
771  integral.Simpson(this, &G4KineticTrack::IntegrandFunction3, theLowerLimit, theUpperLimit, nIterations);
772 
773  return result;
774 }
775 
777 {
778  const G4double theUpperLimit = theActualMass - theDaughterMass[0];
779  const G4int nIterations = 100;
780 
781  if (theLowerLimit>=theUpperLimit) return 0.0;
782 
784  G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction1,
785  theLowerLimit, theUpperLimit, nIterations);
786  return theIntegralOverMass2;
787 }
788 
789 G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit, const G4double poleMass) const
790 {
791  const G4double theUpperLimit = poleMass - theDaughterMass[0];
792  const G4int nIterations = 100;
793 
794  if (theLowerLimit>=theUpperLimit) return 0.0;
795 
797  const G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction2,
798  theLowerLimit, theUpperLimit, nIterations);
799  return theIntegralOverMass2;
800 }
801 
802 
804 {
805  const G4double theLowerLimit = 0.0;
806  const G4double theUpperLimit = theActualMass;
807  const G4int nIterations = 100;
808 
809  if (theLowerLimit>=theUpperLimit) return 0.0;
810 
812  G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction4,
813  theLowerLimit, theUpperLimit, nIterations);
814  return theIntegralOverMass2;
815 }
816 
817 
818 
819 
820 
821 
822 
823 
static G4KaonZero * KaonZero()
Definition: G4KaonZero.cc:104
G4double SampleMass(const G4double poleMass, const G4double gamma, const G4double minMass, const G4double maxMass) const
G4double GetBR() const
G4double IntegrateCMMomentum2() const
void G4SwapObj(T *a, T *b)
Definition: templates.hh:129
G4int operator==(const G4KineticTrack &right) const
G4double theFormationTime
G4int GetnChannels() const
std::ostringstream G4ExceptionDescription
Definition: globals.hh:76
G4LorentzVector theTotal4Momentum
CLHEP::Hep3Vector G4ThreeVector
CLHEP::HepLorentzRotation G4LorentzRotation
const G4String & GetParentName() const
G4int GetNumberOfDaughters() const
G4double * theDaughterWidth
const G4ThreeVector & GetPosition() const
G4KineticTrackVector * Decay()
G4ParticleDefinition * GetDaughter(G4int anIndex)
static G4ThreadLocal G4double G4KineticTrack_Gmass
static G4KaonZeroLong * KaonZeroLong()
G4VDecayChannel * GetDecayChannel(G4int index) const
G4double GetActualMass() const
G4ParticleDefinition * GetDefinition() const
G4ThreeVector thePosition
#define G4ThreadLocal
Definition: tls.hh:89
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
G4double IntegrateCMMomentum(const G4double lowerLimit) const
G4double IntegrandFunction3(G4double xmass) const
G4double * theDaughterMass
G4LorentzVector theFermi3Momentum
G4int entries() const
G4DecayTable * GetDecayTable() const
G4double IntegrandFunction2(G4double xmass) const
G4int operator!=(const G4KineticTrack &right) const
G4double GetPDGWidth() const
#define position
Definition: xmlparse.cc:622
#define G4UniformRand()
Definition: Randomize.hh:97
G4double IntegrandFunction4(G4double xmass) const
static G4KaonZeroShort * KaonZeroShort()
G4bool nucleon(G4int ityp)
G4double GetFormationTime() const
bool G4bool
Definition: G4Types.hh:79
G4LorentzVector the4Momentum
G4double EvaluateTotalActualWidth()
G4KineticTrack & operator=(const G4KineticTrack &right)
static G4ThreadLocal G4double G4KineticTrack_xmass1
const G4int n
void Set4Momentum(const G4LorentzVector &a4Momentum)
G4LorentzVector Get4Momentum() const
const G4String & GetDaughterName(G4int anIndex) const
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double GetMinimumMass(const G4ParticleDefinition *p) const
virtual G4DecayProducts * DecayIt(G4double mass=0.0)
G4double EvaluateCMMomentum(const G4double mass, const G4double *m_ij) const
static G4AntiKaonZero * AntiKaonZero()
G4double GetPDGMass() const
int G4lrint(double ad)
Definition: templates.hh:163
T max(const T t1, const T t2)
brief Return the largest of the two arguments
const G4LorentzVector & GetTrackingMomentum() const
G4DynamicParticle * PopProducts()
G4double * theActualWidth
G4double IntegrandFunction1(G4double xmass) const
#define G4endl
Definition: G4ios.hh:61
G4double Simpson(T &typeT, F f, G4double a, G4double b, G4int n)
G4int entries() const
G4Nucleon * theNucleon
G4double theActualMass
double G4double
Definition: G4Types.hh:76
const G4LorentzVector & Get4Momentum() const
G4double GetPDGCharge() const
G4double theProjectilePotential
const G4ParticleDefinition * GetDefinition() const
G4double BrWig(const G4double Gamma, const G4double rmass, const G4double mass) const
G4GLOB_DLL std::ostream G4cerr
CLHEP::HepLorentzVector G4LorentzVector
const G4ParticleDefinition * theDefinition
CascadeState theStateToNucleus