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

Last change on this file since 1347 was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
File size: 28.2 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*/
[1347]495 G4ParticleDefinition* thisDefinition = this->GetDefinition();
496 if(!thisDefinition)
497 {
498 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
499 G4cerr << " track has no particle definition associated."<<G4endl;
500 return 0;
501 }
502 G4DecayTable* theDecayTable = thisDefinition->GetDecayTable();
503 if(!theDecayTable)
504 {
505 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
506 G4cerr << " particle definiton has no decay table associated."<<G4endl;
507 G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
508 return 0;
509 }
[819]510
511 G4int chargeBalance = G4lrint(theDefinition->GetPDGCharge() );
512 G4int baryonBalance = G4lrint(theDefinition->GetBaryonNumber() );
513 G4LorentzVector energyMomentumBalance(Get4Momentum());
514 G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
515 if (theTotalActualWidth !=0)
516 {
517 G4int index;
518 G4double theSumActualWidth = 0.0;
519 G4double* theCumActualWidth = new G4double[nChannels];
520 for (index = nChannels - 1; index >= 0; index--)
521 {
522 theSumActualWidth += theActualWidth[index];
523 theCumActualWidth[index] = theSumActualWidth;
524 // cout << "DECAY Cum. Width " << index << " " << theCumActualWidth[index] << G4endl;
525 }
526 // cout << "DECAY Total Width " << theSumActualWidth << G4endl;
527 // cout << "DECAY Total Width " << theTotalActualWidth << G4endl;
528 G4double r = theTotalActualWidth * G4UniformRand();
[1347]529 G4VDecayChannel* theDecayChannel(0);
[819]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 }
[1347]540
541 delete [] theCumActualWidth;
542
[819]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
688G4double 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
702G4double 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
715G4double 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
729G4double 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
749G4double 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
762G4double 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
776G4double 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
Note: See TracBrowser for help on using the repository browser.