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

Last change on this file since 967 was 819, checked in by garnier, 17 years ago

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