source: trunk/source/processes/hadronic/models/util/src/G4KineticTrack.cc@ 1330

Last change on this file since 1330 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
File size: 28.3 KB
RevLine 
[819]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
61static G4double G4KineticTrack_Gmass, G4KineticTrack_xmass1;
62
63//
64// Default constructor
65//
66
67G4KineticTrack::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{
83////////////////
84// DEBUG //
85////////////////
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
100G4KineticTrack::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
122////////////////
123// DEBUG //
124////////////////
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
138G4KineticTrack::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
395////////////////
396// DEBUG //
397////////////////
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
413G4KineticTrack::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
435G4KineticTrack::~G4KineticTrack()
436{
437 if (theActualWidth != 0) delete [] theActualWidth;
438 if (theDaughterMass != 0) delete [] theDaughterMass;
439 if (theDaughterWidth != 0) delete [] theDaughterWidth;
440}
441
442
443
444const G4KineticTrack& G4KineticTrack::operator=(const G4KineticTrack& right)
445{
446 G4int i;
447 if (this != &right)
448 {
449 theDefinition = right.GetDefinition();
450 theFormationTime = right.GetFormationTime();
451 the4Momentum = right.the4Momentum;
452 the4Momentum = right.GetTrackingMomentum();
453 theFermi3Momentum = right.theFermi3Momentum;
454 theTotal4Momentum = right.theTotal4Momentum;
455 theNucleon=right.theNucleon;
456 theStateToNucleus=right.theStateToNucleus;
457 if (theActualWidth != 0) delete [] theActualWidth;
458 nChannels = right.GetnChannels();
459 theActualWidth = new G4double[nChannels];
460 for (i = 0; i < nChannels; i++)
461 {
462 theActualWidth[i] = right.theActualWidth[i];
463 }
464 }
465 return *this;
466}
467
468
469
470G4int G4KineticTrack::operator==(const G4KineticTrack& right) const
471{
472 return (this == & right);
473}
474
475
476
477G4int G4KineticTrack::operator!=(const G4KineticTrack& right) const
478{
479 return (this != & right);
480}
481
482
483
484G4KineticTrackVector* G4KineticTrack::Decay()
485{
486//
487// Select a possible decay channel
488//
[1196]489/*
490 G4int index1;
491 for (index1 = nChannels - 1; index1 >= 0; index1--)
492 G4cout << "DECAY Actual Width IND/ActualW " << index1 << " " << theActualWidth[index1] << G4endl;
493 G4cout << "DECAY Actual Mass " << theActualMass << G4endl;
494*/
[819]495
496 G4int chargeBalance = G4lrint(theDefinition->GetPDGCharge() );
497 G4int baryonBalance = G4lrint(theDefinition->GetBaryonNumber() );
498 G4LorentzVector energyMomentumBalance(Get4Momentum());
499 G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
500 if (theTotalActualWidth !=0)
501 {
502 G4int index;
503 G4double theSumActualWidth = 0.0;
504 G4double* theCumActualWidth = new G4double[nChannels];
505 for (index = nChannels - 1; index >= 0; index--)
506 {
507 theSumActualWidth += theActualWidth[index];
508 theCumActualWidth[index] = theSumActualWidth;
509 // cout << "DECAY Cum. Width " << index << " " << theCumActualWidth[index] << G4endl;
510 }
511 // cout << "DECAY Total Width " << theSumActualWidth << G4endl;
512 // cout << "DECAY Total Width " << theTotalActualWidth << G4endl;
513 G4double r = theTotalActualWidth * G4UniformRand();
514 G4ParticleDefinition* thisDefinition = this->GetDefinition();
515 if(!thisDefinition)
516 {
517 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
518 G4cerr << " track has no particle definition associated."<<G4endl;
519 return 0;
520 }
521 G4DecayTable* theDecayTable = thisDefinition->GetDecayTable();
522 if(!theDecayTable)
523 {
524 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
525 G4cerr << " particle definiton has no decay table associated."<<G4endl;
526 G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
527 return 0;
528 }
529 G4VDecayChannel* theDecayChannel=NULL;
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 if(!theDecayChannel)
542 {
543 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
544 G4cerr << " decay channel has 0x0 channel associated."<<G4endl;
545 G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
546 G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
547 return 0;
548 }
549 G4String theParentName = theDecayChannel->GetParentName();
550 G4double theParentMass = this->GetActualMass();
551 G4double theBR = theActualWidth[index];
552 // cout << "**BR*** DECAYNEW " << theBR << G4endl;
553 G4int theNumberOfDaughters = theDecayChannel->GetNumberOfDaughters();
554 G4String theDaughtersName1 = "";
555 G4String theDaughtersName2 = "";
556 G4String theDaughtersName3 = "";
557
558 G4double masses[3]={0.,0.,0.};
559 G4int shortlivedDaughters[3];
560 G4int numberOfShortliveds(0);
561 G4double SumLongLivedMass(0);
562 for (G4int aD=0; aD < theNumberOfDaughters ; aD++)
563 {
564 G4ParticleDefinition* aDaughter = theDecayChannel->GetDaughter(aD);
565 masses[aD] = aDaughter->GetPDGMass();
566 if ( aDaughter->IsShortLived() )
567 {
568 shortlivedDaughters[numberOfShortliveds]=aD;
569 numberOfShortliveds++;
570 } else {
571 SumLongLivedMass += aDaughter->GetPDGMass();
572 }
573
574 }
575 switch (theNumberOfDaughters)
576 {
577 case 0:
578 break;
579 case 1:
580 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
581 theDaughtersName2 = "";
582 theDaughtersName3 = "";
583 break;
584 case 2:
585 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
586 theDaughtersName2 = theDecayChannel->GetDaughterName(1);
587 theDaughtersName3 = "";
588 if ( numberOfShortliveds == 1)
589 { G4SampleResonance aSampler;
590 G4double massmax=theParentMass - SumLongLivedMass;
591 G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
592 masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
593 } else if ( numberOfShortliveds == 2) {
594 // choose masses one after the other, start with randomly choosen
595 G4int zero= (G4UniformRand() > 0.5) ? 0 : 1;
596 G4int one = 1-zero;
597 G4SampleResonance aSampler;
598 G4double massmax=theParentMass - aSampler.GetMinimumMass(theDecayChannel->GetDaughter(shortlivedDaughters[one]));
599 G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[zero]);
600 masses[shortlivedDaughters[zero]]=aSampler.SampleMass(aDaughter,massmax);
601 massmax=theParentMass - masses[shortlivedDaughters[zero]];
602 aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[one]);
603 masses[shortlivedDaughters[one]]=aSampler.SampleMass(aDaughter,massmax);
604 }
605 break;
606 default:
607 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
608 theDaughtersName2 = theDecayChannel->GetDaughterName(1);
609 theDaughtersName3 = theDecayChannel->GetDaughterName(2);
610 if ( numberOfShortliveds == 1)
611 { G4SampleResonance aSampler;
612 G4double massmax=theParentMass - SumLongLivedMass;
613 G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
614 masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
615 }
616 break;
617 }
618
619//
620// Get the decay products List
621//
622
623 G4GeneralPhaseSpaceDecay thePhaseSpaceDecayChannel(theParentName,
624 theParentMass,
625 theBR,
626 theNumberOfDaughters,
627 theDaughtersName1,
628 theDaughtersName2,
629 theDaughtersName3,
630 masses);
631 G4DecayProducts* theDecayProducts = thePhaseSpaceDecayChannel.DecayIt();
632 if(!theDecayProducts)
633 {
634 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
635 G4cerr << " phase-space decay failed."<<G4endl;
636 G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
637 G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
638 G4cerr << " "<<theNumberOfDaughters<< " Daughter particles: "
639 << theDaughtersName1<<" "<<theDaughtersName2<<" "<<theDaughtersName3<<G4endl;
640 return 0;
641 }
642
643//
644// Create the kinetic track List associated to the decay products
645//
646 G4LorentzRotation toMoving(Get4Momentum().boostVector());
647 G4DynamicParticle* theDynamicParticle;
648 G4double formationTime = 0.0;
649 G4ThreeVector position = this->GetPosition();
650 G4LorentzVector momentum;
651 G4LorentzVector momentumBalanceCMS(0);
652 G4KineticTrackVector* theDecayProductList = new G4KineticTrackVector;
653 G4int dEntries = theDecayProducts->entries();
654 G4ParticleDefinition * aProduct = 0;
655 for (G4int i=dEntries; i > 0; i--)
656 {
657 theDynamicParticle = theDecayProducts->PopProducts();
658 aProduct = theDynamicParticle->GetDefinition();
659 chargeBalance -= G4lrint(aProduct->GetPDGCharge() );
660 baryonBalance -= G4lrint(aProduct->GetBaryonNumber() );
661 momentumBalanceCMS += theDynamicParticle->Get4Momentum();
662 momentum = toMoving*theDynamicParticle->Get4Momentum();
663 energyMomentumBalance -= momentum;
664 theDecayProductList->push_back(new G4KineticTrack (aProduct,
665 formationTime,
666 position,
667 momentum));
668 delete theDynamicParticle;
669 }
670 delete theDecayProducts;
671 delete [] theCumActualWidth;
672 if(getenv("DecayEnergyBalanceCheck"))
673 std::cout << "DEBUGGING energy balance in cms and lab, charge baryon balance : "
674 << momentumBalanceCMS << " "
675 <<energyMomentumBalance << " "
676 <<chargeBalance<<" "
677 <<baryonBalance<<" "
678 <<G4endl;
679 return theDecayProductList;
680 }
681 else
682 {
683 return 0;
684 }
685}
686
687G4double G4KineticTrack::IntegrandFunction1(G4double xmass) const
688{
689 G4double mass = theActualMass; /* the actual mass value */
690 G4double mass1 = theDaughterMass[0];
691 G4double mass2 = theDaughterMass[1];
692 G4double gamma2 = theDaughterWidth[1];
693
694 G4double result = (1. / (2 * mass)) *
695 std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
696 ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
697 BrWig(gamma2, mass2, xmass);
698 return result;
699}
700
701G4double G4KineticTrack::IntegrandFunction2(G4double xmass) const
702{
703 G4double mass = theDefinition->GetPDGMass(); /* the pole mass value */
704 G4double mass1 = theDaughterMass[0];
705 G4double mass2 = theDaughterMass[1];
706 G4double gamma2 = theDaughterWidth[1];
707 G4double result = (1. / (2 * mass)) *
708 std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
709 ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
710 BrWig(gamma2, mass2, xmass);
711 return result;
712}
713
714G4double G4KineticTrack::IntegrandFunction3(G4double xmass) const
715{
716 const G4double mass = G4KineticTrack_Gmass; /* the actual mass value */
717// const G4double mass1 = theDaughterMass[0];
718 const G4double mass2 = theDaughterMass[1];
719 const G4double gamma2 = theDaughterWidth[1];
720
721 const G4double result = (1. / (2 * mass)) *
722 std::sqrt(((mass * mass) - (G4KineticTrack_xmass1 + xmass) * (G4KineticTrack_xmass1 + xmass)) *
723 ((mass * mass) - (G4KineticTrack_xmass1 - xmass) * (G4KineticTrack_xmass1 - xmass))) *
724 BrWig(gamma2, mass2, xmass);
725 return result;
726}
727
728G4double G4KineticTrack::IntegrandFunction4(G4double xmass) const
729{
730 const G4double mass = G4KineticTrack_Gmass;
731 const G4double mass1 = theDaughterMass[0];
732 const G4double gamma1 = theDaughterWidth[0];
733// const G4double mass2 = theDaughterMass[1];
734
735 G4KineticTrack_xmass1 = xmass;
736
737 const G4double theLowerLimit = 0.0;
738 const G4double theUpperLimit = mass - xmass;
739 const G4int nIterations = 100;
740
741 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
742 G4double result = BrWig(gamma1, mass1, xmass)*
743 integral.Simpson(this, &G4KineticTrack::IntegrandFunction3, theLowerLimit, theUpperLimit, nIterations);
744
745 return result;
746}
747
748G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit) const
749{
750 const G4double theUpperLimit = theActualMass - theDaughterMass[0];
751 const G4int nIterations = 100;
752
753 if (theLowerLimit>=theUpperLimit) return 0.0;
754
755 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
756 G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction1,
757 theLowerLimit, theUpperLimit, nIterations);
758 return theIntegralOverMass2;
759}
760
761G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit, const G4double poleMass) const
762{
763 const G4double theUpperLimit = poleMass - theDaughterMass[0];
764 const G4int nIterations = 100;
765
766 if (theLowerLimit>=theUpperLimit) return 0.0;
767
768 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
769 const G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction2,
770 theLowerLimit, theUpperLimit, nIterations);
771 return theIntegralOverMass2;
772}
773
774
775G4double G4KineticTrack::IntegrateCMMomentum2() const
776{
777 const G4double theLowerLimit = 0.0;
778 const G4double theUpperLimit = theActualMass;
779 const G4int nIterations = 100;
780
781 if (theLowerLimit>=theUpperLimit) return 0.0;
782
783 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
784 G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction4,
785 theLowerLimit, theUpperLimit, nIterations);
786 return theIntegralOverMass2;
787}
788
789
790
791
792
793
794
795
Note: See TracBrowser for help on using the repository browser.