Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
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"
45 #include "G4GeneralPhaseSpaceDecay.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 
61 static G4double G4KineticTrack_Gmass, G4KineticTrack_xmass1;
62 
63 //
64 // Default constructor
65 //
66 
67 G4KineticTrack::G4KineticTrack() :
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 
100 G4KineticTrack::G4KineticTrack(const G4KineticTrack &right) : G4VKineticNucleon()
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 }
132 
133 
134 //
135 // By argument constructor
136 //
137 
138 G4KineticTrack::G4KineticTrack(G4ParticleDefinition* aDefinition,
139  G4double aFormationTime,
140  G4ThreeVector aPosition,
141  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 {
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  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  int 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  G4double l=0;
374  //if(nDaughters<3) theChannel->GetAngularMomentum();
375  G4double theMassRatio = thePoleMass / theActualMass;
376  G4double theMomRatio = theActualMom / thePoleMom;
377  theActualWidth[index] = thePoleWidth * theMassRatio *
378  std::pow(theMomRatio, (2 * l + 1)) *
379  (1.2 / (1+ 0.2*std::pow(theMomRatio, (2 * l))));
380  delete [] theDaughterMass;
381  theDaughterMass = 0;
382  delete [] theDaughterWidth;
383  theDaughterWidth = 0;
384  delete [] theDaughterIsShortLived;
385  theDaughterIsShortLived = 0;
386  }
387 
388  else // nDaughter = 1 ( e.g. K0 decays 50% to Kshort, 50% Klong
389  {
390  theMotherWidth = theDefinition->GetPDGWidth();
391  theActualWidth[index] = theChannel->GetBR()*theMotherWidth;
392  }
393  }
394 
396 // DEBUG //
398 
399 // for (G4int y = nChannels - 1; y >= 0; y--)
400 // {
401 // G4cout << G4endl << theActualWidth[y];
402 // }
403 // G4cout << G4endl << G4endl << G4endl;
404 
405  /*
406  G4cerr << G4endl << G4endl << G4endl;
407  G4cerr << " G4KineticTrack by argument constructor invoked! \n";
408  G4cerr << " =============================================== \n" << G4endl;
409  */
410 
411 }
412 
413 G4KineticTrack::G4KineticTrack(G4Nucleon * nucleon,
414  G4ThreeVector aPosition,
415  G4LorentzVector& a4Momentum)
416  : theDefinition(nucleon->GetDefinition()),
417  theFormationTime(0),
418  thePosition(aPosition),
419  the4Momentum(a4Momentum),
420  theFermi3Momentum(nucleon->GetMomentum()),
421  theNucleon(nucleon),
422  nChannels(0),
423  theActualMass(nucleon->GetDefinition()->GetPDGMass()),
424  theActualWidth(0),
425  theDaughterMass(0),
426  theDaughterWidth(0),
427  theStateToNucleus(undefined),
428  theProjectilePotential(0)
429 {
430  theFermi3Momentum.setE(0);
431  Set4Momentum(a4Momentum);
432 }
433 
434 
435 G4KineticTrack::~G4KineticTrack()
436 {
437  if (theActualWidth != 0) delete [] theActualWidth;
438  if (theDaughterMass != 0) delete [] theDaughterMass;
439  if (theDaughterWidth != 0) delete [] theDaughterWidth;
440 }
441 
442 
443 
444 G4KineticTrack& G4KineticTrack::operator=(const G4KineticTrack& right)
445 {
446  if (this != &right)
447  {
448  theDefinition = right.GetDefinition();
449  theFormationTime = right.GetFormationTime();
450  the4Momentum = right.the4Momentum;
451  the4Momentum = right.GetTrackingMomentum();
452  theFermi3Momentum = right.theFermi3Momentum;
453  theTotal4Momentum = right.theTotal4Momentum;
454  theNucleon=right.theNucleon;
455  theStateToNucleus=right.theStateToNucleus;
456  if (theActualWidth != 0) delete [] theActualWidth;
457  nChannels = right.GetnChannels();
458  theActualWidth = new G4double[nChannels];
459  for ( G4int i = 0; i < nChannels; i++)
460  {
461  theActualWidth[i] = right.theActualWidth[i];
462  }
463  }
464  return *this;
465 }
466 
467 
468 
469 G4int G4KineticTrack::operator==(const G4KineticTrack& right) const
470 {
471  return (this == & right);
472 }
473 
474 
475 
476 G4int G4KineticTrack::operator!=(const G4KineticTrack& right) const
477 {
478  return (this != & right);
479 }
480 
481 
482 
483 G4KineticTrackVector* G4KineticTrack::Decay()
484 {
485 //
486 // Select a possible decay channel
487 //
488 /*
489  G4int index1;
490  for (index1 = nChannels - 1; index1 >= 0; index1--)
491  G4cout << "DECAY Actual Width IND/ActualW " << index1 << " " << theActualWidth[index1] << G4endl;
492  G4cout << "DECAY Actual Mass " << theActualMass << G4endl;
493 */
494  G4ParticleDefinition* thisDefinition = this->GetDefinition();
495  if(!thisDefinition)
496  {
497  G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
498  G4cerr << " track has no particle definition associated."<<G4endl;
499  return 0;
500  }
501  G4DecayTable* theDecayTable = thisDefinition->GetDecayTable();
502  if(!theDecayTable)
503  {
504  G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
505  G4cerr << " particle definiton has no decay table associated."<<G4endl;
506  G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
507  return 0;
508  }
509 
510  G4int chargeBalance = G4lrint(theDefinition->GetPDGCharge() );
511  G4int baryonBalance = G4lrint(theDefinition->GetBaryonNumber() );
512  G4LorentzVector energyMomentumBalance(Get4Momentum());
513  G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
514  if (theTotalActualWidth !=0)
515  {
516  G4int index;
517  G4double theSumActualWidth = 0.0;
518  G4double* theCumActualWidth = new G4double[nChannels];
519  for (index = nChannels - 1; index >= 0; index--)
520  {
521  theSumActualWidth += theActualWidth[index];
522  theCumActualWidth[index] = theSumActualWidth;
523  // cout << "DECAY Cum. Width " << index << " " << theCumActualWidth[index] << G4endl;
524  }
525  // cout << "DECAY Total Width " << theSumActualWidth << G4endl;
526  // cout << "DECAY Total Width " << theTotalActualWidth << G4endl;
527  G4double r = theTotalActualWidth * G4UniformRand();
528  G4VDecayChannel* theDecayChannel(0);
529  G4int chosench=-1;
530  for (index = nChannels - 1; index >= 0; index--)
531  {
532  if (r < theCumActualWidth[index])
533  {
534  theDecayChannel = theDecayTable->GetDecayChannel(index);
535  // cout << "DECAY SELECTED CHANNEL" << index << G4endl;
536  chosench=index;
537  break;
538  }
539  }
540 
541  delete [] theCumActualWidth;
542 
543  if(!theDecayChannel)
544  {
545  G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
546  G4cerr << " decay channel has 0x0 channel associated."<<G4endl;
547  G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
548  G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
549  return 0;
550  }
551  G4String theParentName = theDecayChannel->GetParentName();
552  G4double theParentMass = this->GetActualMass();
553  G4double theBR = theActualWidth[index];
554  // cout << "**BR*** DECAYNEW " << theBR << G4endl;
555  G4int theNumberOfDaughters = theDecayChannel->GetNumberOfDaughters();
556  G4String theDaughtersName1 = "";
557  G4String theDaughtersName2 = "";
558  G4String theDaughtersName3 = "";
559 
560  G4double masses[3]={0.,0.,0.};
561  G4int shortlivedDaughters[3];
562  G4int numberOfShortliveds(0);
563  G4double SumLongLivedMass(0);
564  for (G4int aD=0; aD < theNumberOfDaughters ; aD++)
565  {
566  G4ParticleDefinition* aDaughter = theDecayChannel->GetDaughter(aD);
567  masses[aD] = aDaughter->GetPDGMass();
568  if ( aDaughter->IsShortLived() )
569  {
570  shortlivedDaughters[numberOfShortliveds]=aD;
571  numberOfShortliveds++;
572  } else {
573  SumLongLivedMass += aDaughter->GetPDGMass();
574  }
575 
576  }
577  switch (theNumberOfDaughters)
578  {
579  case 0:
580  break;
581  case 1:
582  theDaughtersName1 = theDecayChannel->GetDaughterName(0);
583  theDaughtersName2 = "";
584  theDaughtersName3 = "";
585  break;
586  case 2:
587  theDaughtersName1 = theDecayChannel->GetDaughterName(0);
588  theDaughtersName2 = theDecayChannel->GetDaughterName(1);
589  theDaughtersName3 = "";
590  if ( numberOfShortliveds == 1)
591  { G4SampleResonance aSampler;
592  G4double massmax=theParentMass - SumLongLivedMass;
593  G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
594  masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
595  } else if ( numberOfShortliveds == 2) {
596  // choose masses one after the other, start with randomly choosen
597  G4int zero= (G4UniformRand() > 0.5) ? 0 : 1;
598  G4int one = 1-zero;
599  G4SampleResonance aSampler;
600  G4double massmax=theParentMass - aSampler.GetMinimumMass(theDecayChannel->GetDaughter(shortlivedDaughters[one]));
601  G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[zero]);
602  masses[shortlivedDaughters[zero]]=aSampler.SampleMass(aDaughter,massmax);
603  massmax=theParentMass - masses[shortlivedDaughters[zero]];
604  aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[one]);
605  masses[shortlivedDaughters[one]]=aSampler.SampleMass(aDaughter,massmax);
606  }
607  break;
608  default:
609  theDaughtersName1 = theDecayChannel->GetDaughterName(0);
610  theDaughtersName2 = theDecayChannel->GetDaughterName(1);
611  theDaughtersName3 = theDecayChannel->GetDaughterName(2);
612  if ( numberOfShortliveds == 1)
613  { G4SampleResonance aSampler;
614  G4double massmax=theParentMass - SumLongLivedMass;
615  G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
616  masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
617  }
618  break;
619  }
620 
621 //
622 // Get the decay products List
623 //
624 
625  G4GeneralPhaseSpaceDecay thePhaseSpaceDecayChannel(theParentName,
626  theParentMass,
627  theBR,
628  theNumberOfDaughters,
629  theDaughtersName1,
630  theDaughtersName2,
631  theDaughtersName3,
632  masses);
633  G4DecayProducts* theDecayProducts = thePhaseSpaceDecayChannel.DecayIt();
634  if(!theDecayProducts)
635  {
636  G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
637  G4cerr << " phase-space decay failed."<<G4endl;
638  G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
639  G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
640  G4cerr << " "<<theNumberOfDaughters<< " Daughter particles: "
641  << theDaughtersName1<<" "<<theDaughtersName2<<" "<<theDaughtersName3<<G4endl;
642  return 0;
643  }
644 
645 //
646 // Create the kinetic track List associated to the decay products
647 //
648  G4LorentzRotation toMoving(Get4Momentum().boostVector());
649  G4DynamicParticle* theDynamicParticle;
650  G4double formationTime = 0.0;
651  G4ThreeVector position = this->GetPosition();
652  G4LorentzVector momentum;
653  G4LorentzVector momentumBalanceCMS(0);
654  G4KineticTrackVector* theDecayProductList = new G4KineticTrackVector;
655  G4int dEntries = theDecayProducts->entries();
656  G4ParticleDefinition * aProduct = 0;
657  for (G4int i=dEntries; i > 0; i--)
658  {
659  theDynamicParticle = theDecayProducts->PopProducts();
660  aProduct = theDynamicParticle->GetDefinition();
661  chargeBalance -= G4lrint(aProduct->GetPDGCharge() );
662  baryonBalance -= G4lrint(aProduct->GetBaryonNumber() );
663  momentumBalanceCMS += theDynamicParticle->Get4Momentum();
664  momentum = toMoving*theDynamicParticle->Get4Momentum();
665  energyMomentumBalance -= momentum;
666  theDecayProductList->push_back(new G4KineticTrack (aProduct,
667  formationTime,
668  position,
669  momentum));
670  delete theDynamicParticle;
671  }
672  delete theDecayProducts;
673  if(getenv("DecayEnergyBalanceCheck"))
674  std::cout << "DEBUGGING energy balance in cms and lab, charge baryon balance : "
675  << momentumBalanceCMS << " "
676  <<energyMomentumBalance << " "
677  <<chargeBalance<<" "
678  <<baryonBalance<<" "
679  <<G4endl;
680  return theDecayProductList;
681  }
682  else
683  {
684  return 0;
685  }
686 }
687 
688 G4double G4KineticTrack::IntegrandFunction1(G4double xmass) const
689 {
690  G4double mass = theActualMass; /* the actual mass value */
691  G4double mass1 = theDaughterMass[0];
692  G4double mass2 = theDaughterMass[1];
693  G4double gamma2 = theDaughterWidth[1];
694 
695  G4double result = (1. / (2 * mass)) *
696  std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
697  ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
698  BrWig(gamma2, mass2, xmass);
699  return result;
700 }
701 
702 G4double G4KineticTrack::IntegrandFunction2(G4double xmass) const
703 {
704  G4double mass = theDefinition->GetPDGMass(); /* the pole mass value */
705  G4double mass1 = theDaughterMass[0];
706  G4double mass2 = theDaughterMass[1];
707  G4double gamma2 = theDaughterWidth[1];
708  G4double result = (1. / (2 * mass)) *
709  std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
710  ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
711  BrWig(gamma2, mass2, xmass);
712  return result;
713 }
714 
715 G4double G4KineticTrack::IntegrandFunction3(G4double xmass) const
716 {
717  const G4double mass = G4KineticTrack_Gmass; /* the actual mass value */
718 // const G4double mass1 = theDaughterMass[0];
719  const G4double mass2 = theDaughterMass[1];
720  const G4double gamma2 = theDaughterWidth[1];
721 
722  const G4double result = (1. / (2 * mass)) *
723  std::sqrt(((mass * mass) - (G4KineticTrack_xmass1 + xmass) * (G4KineticTrack_xmass1 + xmass)) *
724  ((mass * mass) - (G4KineticTrack_xmass1 - xmass) * (G4KineticTrack_xmass1 - xmass))) *
725  BrWig(gamma2, mass2, xmass);
726  return result;
727 }
728 
729 G4double G4KineticTrack::IntegrandFunction4(G4double xmass) const
730 {
731  const G4double mass = G4KineticTrack_Gmass;
732  const G4double mass1 = theDaughterMass[0];
733  const G4double gamma1 = theDaughterWidth[0];
734 // const G4double mass2 = theDaughterMass[1];
735 
736  G4KineticTrack_xmass1 = xmass;
737 
738  const G4double theLowerLimit = 0.0;
739  const G4double theUpperLimit = mass - xmass;
740  const G4int nIterations = 100;
741 
742  G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
743  G4double result = BrWig(gamma1, mass1, xmass)*
744  integral.Simpson(this, &G4KineticTrack::IntegrandFunction3, theLowerLimit, theUpperLimit, nIterations);
745 
746  return result;
747 }
748 
749 G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit) const
750 {
751  const G4double theUpperLimit = theActualMass - theDaughterMass[0];
752  const G4int nIterations = 100;
753 
754  if (theLowerLimit>=theUpperLimit) return 0.0;
755 
756  G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
757  G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction1,
758  theLowerLimit, theUpperLimit, nIterations);
759  return theIntegralOverMass2;
760 }
761 
762 G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit, const G4double poleMass) const
763 {
764  const G4double theUpperLimit = poleMass - theDaughterMass[0];
765  const G4int nIterations = 100;
766 
767  if (theLowerLimit>=theUpperLimit) return 0.0;
768 
769  G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
770  const G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction2,
771  theLowerLimit, theUpperLimit, nIterations);
772  return theIntegralOverMass2;
773 }
774 
775 
776 G4double G4KineticTrack::IntegrateCMMomentum2() const
777 {
778  const G4double theLowerLimit = 0.0;
779  const G4double theUpperLimit = theActualMass;
780  const G4int nIterations = 100;
781 
782  if (theLowerLimit>=theUpperLimit) return 0.0;
783 
784  G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
785  G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction4,
786  theLowerLimit, theUpperLimit, nIterations);
787  return theIntegralOverMass2;
788 }
789 
790 
791 
792 
793 
794 
795 
796 
Generated on Sat May 25 2013 14:34:08 for Geant4 by   doxygen 1.8.4