00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033 #include "globals.hh"
00034 #include "G4ios.hh"
00035 #include <cmath>
00036
00037 #include "Randomize.hh"
00038 #include "G4SimpleIntegration.hh"
00039 #include "G4ThreeVector.hh"
00040 #include "G4LorentzVector.hh"
00041 #include "G4KineticTrack.hh"
00042 #include "G4KineticTrackVector.hh"
00043 #include "G4ParticleDefinition.hh"
00044 #include "G4DecayTable.hh"
00045 #include "G4GeneralPhaseSpaceDecay.hh"
00046 #include "G4DecayProducts.hh"
00047 #include "G4LorentzRotation.hh"
00048 #include "G4SampleResonance.hh"
00049 #include "G4Integrator.hh"
00050 #include "G4KaonZero.hh"
00051 #include "G4KaonZeroShort.hh"
00052 #include "G4KaonZeroLong.hh"
00053 #include "G4AntiKaonZero.hh"
00054
00055 #include "G4HadTmpUtil.hh"
00056
00057
00058
00059
00060
00061 static G4double G4KineticTrack_Gmass, G4KineticTrack_xmass1;
00062
00063
00064
00065
00066
00067 G4KineticTrack::G4KineticTrack() :
00068 theDefinition(0),
00069 theFormationTime(0),
00070 thePosition(0),
00071 the4Momentum(0),
00072 theFermi3Momentum(0),
00073 theTotal4Momentum(0),
00074 theNucleon(0),
00075 nChannels(0),
00076 theActualMass(0),
00077 theActualWidth(0),
00078 theDaughterMass(0),
00079 theDaughterWidth(0),
00080 theStateToNucleus(undefined),
00081 theProjectilePotential(0)
00082 {
00084
00086
00087
00088
00089
00090
00091
00092 }
00093
00094
00095
00096
00097
00098
00099
00100 G4KineticTrack::G4KineticTrack(const G4KineticTrack &right) : G4VKineticNucleon()
00101 {
00102 G4int i;
00103 theDefinition = right.GetDefinition();
00104 theFormationTime = right.GetFormationTime();
00105 thePosition = right.GetPosition();
00106 the4Momentum = right.GetTrackingMomentum();
00107 theFermi3Momentum = right.theFermi3Momentum;
00108 theTotal4Momentum = right.theTotal4Momentum;
00109 theNucleon=right.theNucleon;
00110 nChannels = right.GetnChannels();
00111 theActualMass = right.GetActualMass();
00112 theActualWidth = new G4double[nChannels];
00113 for (i = 0; i < nChannels; i++)
00114 {
00115 theActualWidth[i] = right.theActualWidth[i];
00116 }
00117 theDaughterMass = 0;
00118 theDaughterWidth = 0;
00119 theStateToNucleus=right.theStateToNucleus;
00120 theProjectilePotential=right.theProjectilePotential;
00121
00123
00125
00126
00127
00128
00129
00130
00131 }
00132
00133
00134
00135
00136
00137
00138 G4KineticTrack::G4KineticTrack(G4ParticleDefinition* aDefinition,
00139 G4double aFormationTime,
00140 G4ThreeVector aPosition,
00141 G4LorentzVector& a4Momentum) :
00142 theDefinition(aDefinition),
00143 theFormationTime(aFormationTime),
00144 thePosition(aPosition),
00145 the4Momentum(a4Momentum),
00146 theFermi3Momentum(0),
00147 theTotal4Momentum(a4Momentum),
00148 theNucleon(0),
00149 theStateToNucleus(undefined),
00150 theProjectilePotential(0)
00151 {
00152 if(G4KaonZero::KaonZero() == theDefinition ||
00153 G4AntiKaonZero::AntiKaonZero() == theDefinition)
00154 {
00155 if(G4UniformRand()<0.5)
00156 {
00157 theDefinition = G4KaonZeroShort::KaonZeroShort();
00158 }
00159 else
00160 {
00161 theDefinition = G4KaonZeroLong::KaonZeroLong();
00162 }
00163 }
00164
00165
00166
00167
00168
00169 G4DecayTable* theDecayTable = theDefinition->GetDecayTable();
00170 if (theDecayTable != 0)
00171 {
00172 nChannels = theDecayTable->entries();
00173
00174 }
00175 else
00176 {
00177 nChannels = 0;
00178 }
00179
00180
00181
00182
00183
00184 theActualMass = GetActualMass();
00185
00186
00187
00188
00189
00190
00191 theDaughterMass = 0;
00192 theDaughterWidth = 0;
00193 theActualWidth = 0;
00194 G4bool * theDaughterIsShortLived = 0;
00195
00196 if(nChannels!=0) theActualWidth = new G4double[nChannels];
00197
00198
00199 G4int index;
00200 for (index = nChannels - 1; index >= 0; index--)
00201 {
00202 G4VDecayChannel* theChannel = theDecayTable->GetDecayChannel(index);
00203 G4int nDaughters = theChannel->GetNumberOfDaughters();
00204 G4double theMotherWidth;
00205 if (nDaughters == 2 || nDaughters == 3)
00206 {
00207 G4double thePoleMass = theDefinition->GetPDGMass();
00208 theMotherWidth = theDefinition->GetPDGWidth();
00209 G4double thePoleWidth = theChannel->GetBR()*theMotherWidth;
00210 G4ParticleDefinition* aDaughter;
00211 theDaughterMass = new G4double[nDaughters];
00212 theDaughterWidth = new G4double[nDaughters];
00213 theDaughterIsShortLived = new G4bool[nDaughters];
00214 G4int n;
00215 for (n = 0; n < nDaughters; n++)
00216 {
00217 aDaughter = theChannel->GetDaughter(n);
00218 theDaughterMass[n] = aDaughter->GetPDGMass();
00219 theDaughterWidth[n] = aDaughter->GetPDGWidth();
00220 theDaughterIsShortLived[n] = aDaughter->IsShortLived();
00221 }
00222
00223
00224
00225
00226
00227 G4double theActualMom = 0.0;
00228 G4double thePoleMom = 0.0;
00229 G4SampleResonance aSampler;
00230 if (nDaughters==2)
00231 {
00232 if ( !theDaughterIsShortLived[0] && !theDaughterIsShortLived[1] )
00233 {
00234
00235
00236
00237
00238
00239
00240
00241
00242 theActualMom = EvaluateCMMomentum(theActualMass,
00243 theDaughterMass);
00244 thePoleMom = EvaluateCMMomentum(thePoleMass,
00245 theDaughterMass);
00246
00247
00248
00249
00250
00251 }
00252 else if ( !theDaughterIsShortLived[0] && theDaughterIsShortLived[1] )
00253 {
00254
00255
00256
00257
00258
00259
00260
00261
00262 G4double lowerLimit = aSampler.GetMinimumMass(theChannel->GetDaughter(1));
00263 theActualMom = IntegrateCMMomentum(lowerLimit);
00264 thePoleMom = IntegrateCMMomentum(lowerLimit, thePoleMass);
00265
00266
00267
00268
00269
00270
00271
00272
00273 }
00274 else if ( theDaughterIsShortLived[0] && !theDaughterIsShortLived[1] )
00275 {
00276
00277
00278
00279
00280
00281
00282
00283
00284
00285
00286
00287
00288 G4SwapObj(theDaughterMass, theDaughterMass + 1);
00289 G4SwapObj(theDaughterWidth, theDaughterWidth + 1);
00290
00291
00292 G4double lowerLimit = aSampler.GetMinimumMass(theChannel->GetDaughter(0));
00293 theActualMom = IntegrateCMMomentum(lowerLimit);
00294 thePoleMom = IntegrateCMMomentum(lowerLimit, thePoleMass);
00295
00296
00297
00298
00299
00300
00301
00302
00303 }
00304 else if ( theDaughterIsShortLived[0] && theDaughterIsShortLived[1] )
00305 {
00306
00307
00308
00309
00310
00311
00312
00313
00314
00315 G4KineticTrack_Gmass = theActualMass;
00316 theActualMom = IntegrateCMMomentum2();
00317 G4KineticTrack_Gmass = thePoleMass;
00318 thePoleMom = IntegrateCMMomentum2();
00319
00320
00321
00322
00323
00324
00325
00326 }
00327 }
00328 else
00329 {
00330
00331 int nShortLived = 0;
00332 if ( theDaughterIsShortLived[0] )
00333 {
00334 nShortLived++;
00335 }
00336 if ( theDaughterIsShortLived[1] )
00337 {
00338 nShortLived++;
00339 G4SwapObj(theDaughterMass, theDaughterMass + 1);
00340 G4SwapObj(theDaughterWidth, theDaughterWidth + 1);
00341 }
00342 if ( theDaughterIsShortLived[2] )
00343 {
00344 nShortLived++;
00345 G4SwapObj(theDaughterMass, theDaughterMass + 2);
00346 G4SwapObj(theDaughterWidth, theDaughterWidth + 2);
00347 }
00348 if ( nShortLived == 0 )
00349 {
00350 theDaughterMass[1]+=theDaughterMass[2];
00351 theActualMom = EvaluateCMMomentum(theActualMass,
00352 theDaughterMass);
00353 thePoleMom = EvaluateCMMomentum(thePoleMass,
00354 theDaughterMass);
00355 }
00356
00357 else if ( nShortLived >= 1 )
00358 {
00359
00360 G4SwapObj(theDaughterMass, theDaughterMass + 1);
00361 G4SwapObj(theDaughterWidth, theDaughterWidth + 1);
00362 theDaughterMass[0] += theDaughterMass[2];
00363 theActualMom = IntegrateCMMomentum(0.0);
00364 thePoleMom = IntegrateCMMomentum(0.0, thePoleMass);
00365 }
00366
00367
00368
00369
00370
00371 }
00372
00373 G4double l=0;
00374
00375 G4double theMassRatio = thePoleMass / theActualMass;
00376 G4double theMomRatio = theActualMom / thePoleMom;
00377 theActualWidth[index] = thePoleWidth * theMassRatio *
00378 std::pow(theMomRatio, (2 * l + 1)) *
00379 (1.2 / (1+ 0.2*std::pow(theMomRatio, (2 * l))));
00380 delete [] theDaughterMass;
00381 theDaughterMass = 0;
00382 delete [] theDaughterWidth;
00383 theDaughterWidth = 0;
00384 delete [] theDaughterIsShortLived;
00385 theDaughterIsShortLived = 0;
00386 }
00387
00388 else
00389 {
00390 theMotherWidth = theDefinition->GetPDGWidth();
00391 theActualWidth[index] = theChannel->GetBR()*theMotherWidth;
00392 }
00393 }
00394
00396
00398
00399
00400
00401
00402
00403
00404
00405
00406
00407
00408
00409
00410
00411 }
00412
00413 G4KineticTrack::G4KineticTrack(G4Nucleon * nucleon,
00414 G4ThreeVector aPosition,
00415 G4LorentzVector& a4Momentum)
00416 : theDefinition(nucleon->GetDefinition()),
00417 theFormationTime(0),
00418 thePosition(aPosition),
00419 the4Momentum(a4Momentum),
00420 theFermi3Momentum(nucleon->GetMomentum()),
00421 theNucleon(nucleon),
00422 nChannels(0),
00423 theActualMass(nucleon->GetDefinition()->GetPDGMass()),
00424 theActualWidth(0),
00425 theDaughterMass(0),
00426 theDaughterWidth(0),
00427 theStateToNucleus(undefined),
00428 theProjectilePotential(0)
00429 {
00430 theFermi3Momentum.setE(0);
00431 Set4Momentum(a4Momentum);
00432 }
00433
00434
00435 G4KineticTrack::~G4KineticTrack()
00436 {
00437 if (theActualWidth != 0) delete [] theActualWidth;
00438 if (theDaughterMass != 0) delete [] theDaughterMass;
00439 if (theDaughterWidth != 0) delete [] theDaughterWidth;
00440 }
00441
00442
00443
00444 G4KineticTrack& G4KineticTrack::operator=(const G4KineticTrack& right)
00445 {
00446 if (this != &right)
00447 {
00448 theDefinition = right.GetDefinition();
00449 theFormationTime = right.GetFormationTime();
00450 the4Momentum = right.the4Momentum;
00451 the4Momentum = right.GetTrackingMomentum();
00452 theFermi3Momentum = right.theFermi3Momentum;
00453 theTotal4Momentum = right.theTotal4Momentum;
00454 theNucleon=right.theNucleon;
00455 theStateToNucleus=right.theStateToNucleus;
00456 if (theActualWidth != 0) delete [] theActualWidth;
00457 nChannels = right.GetnChannels();
00458 theActualWidth = new G4double[nChannels];
00459 for ( G4int i = 0; i < nChannels; i++)
00460 {
00461 theActualWidth[i] = right.theActualWidth[i];
00462 }
00463 }
00464 return *this;
00465 }
00466
00467
00468
00469 G4int G4KineticTrack::operator==(const G4KineticTrack& right) const
00470 {
00471 return (this == & right);
00472 }
00473
00474
00475
00476 G4int G4KineticTrack::operator!=(const G4KineticTrack& right) const
00477 {
00478 return (this != & right);
00479 }
00480
00481
00482
00483 G4KineticTrackVector* G4KineticTrack::Decay()
00484 {
00485
00486
00487
00488
00489
00490
00491
00492
00493
00494 G4ParticleDefinition* thisDefinition = this->GetDefinition();
00495 if(!thisDefinition)
00496 {
00497 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
00498 G4cerr << " track has no particle definition associated."<<G4endl;
00499 return 0;
00500 }
00501 G4DecayTable* theDecayTable = thisDefinition->GetDecayTable();
00502 if(!theDecayTable)
00503 {
00504 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
00505 G4cerr << " particle definiton has no decay table associated."<<G4endl;
00506 G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
00507 return 0;
00508 }
00509
00510 G4int chargeBalance = G4lrint(theDefinition->GetPDGCharge() );
00511 G4int baryonBalance = G4lrint(theDefinition->GetBaryonNumber() );
00512 G4LorentzVector energyMomentumBalance(Get4Momentum());
00513 G4double theTotalActualWidth = this->EvaluateTotalActualWidth();
00514 if (theTotalActualWidth !=0)
00515 {
00516 G4int index;
00517 G4double theSumActualWidth = 0.0;
00518 G4double* theCumActualWidth = new G4double[nChannels];
00519 for (index = nChannels - 1; index >= 0; index--)
00520 {
00521 theSumActualWidth += theActualWidth[index];
00522 theCumActualWidth[index] = theSumActualWidth;
00523
00524 }
00525
00526
00527 G4double r = theTotalActualWidth * G4UniformRand();
00528 G4VDecayChannel* theDecayChannel(0);
00529 G4int chosench=-1;
00530 for (index = nChannels - 1; index >= 0; index--)
00531 {
00532 if (r < theCumActualWidth[index])
00533 {
00534 theDecayChannel = theDecayTable->GetDecayChannel(index);
00535
00536 chosench=index;
00537 break;
00538 }
00539 }
00540
00541 delete [] theCumActualWidth;
00542
00543 if(!theDecayChannel)
00544 {
00545 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
00546 G4cerr << " decay channel has 0x0 channel associated."<<G4endl;
00547 G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
00548 G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
00549 return 0;
00550 }
00551 G4String theParentName = theDecayChannel->GetParentName();
00552 G4double theParentMass = this->GetActualMass();
00553 G4double theBR = theActualWidth[index];
00554
00555 G4int theNumberOfDaughters = theDecayChannel->GetNumberOfDaughters();
00556 G4String theDaughtersName1 = "";
00557 G4String theDaughtersName2 = "";
00558 G4String theDaughtersName3 = "";
00559
00560 G4double masses[3]={0.,0.,0.};
00561 G4int shortlivedDaughters[3];
00562 G4int numberOfShortliveds(0);
00563 G4double SumLongLivedMass(0);
00564 for (G4int aD=0; aD < theNumberOfDaughters ; aD++)
00565 {
00566 G4ParticleDefinition* aDaughter = theDecayChannel->GetDaughter(aD);
00567 masses[aD] = aDaughter->GetPDGMass();
00568 if ( aDaughter->IsShortLived() )
00569 {
00570 shortlivedDaughters[numberOfShortliveds]=aD;
00571 numberOfShortliveds++;
00572 } else {
00573 SumLongLivedMass += aDaughter->GetPDGMass();
00574 }
00575
00576 }
00577 switch (theNumberOfDaughters)
00578 {
00579 case 0:
00580 break;
00581 case 1:
00582 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
00583 theDaughtersName2 = "";
00584 theDaughtersName3 = "";
00585 break;
00586 case 2:
00587 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
00588 theDaughtersName2 = theDecayChannel->GetDaughterName(1);
00589 theDaughtersName3 = "";
00590 if ( numberOfShortliveds == 1)
00591 { G4SampleResonance aSampler;
00592 G4double massmax=theParentMass - SumLongLivedMass;
00593 G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
00594 masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
00595 } else if ( numberOfShortliveds == 2) {
00596
00597 G4int zero= (G4UniformRand() > 0.5) ? 0 : 1;
00598 G4int one = 1-zero;
00599 G4SampleResonance aSampler;
00600 G4double massmax=theParentMass - aSampler.GetMinimumMass(theDecayChannel->GetDaughter(shortlivedDaughters[one]));
00601 G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[zero]);
00602 masses[shortlivedDaughters[zero]]=aSampler.SampleMass(aDaughter,massmax);
00603 massmax=theParentMass - masses[shortlivedDaughters[zero]];
00604 aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[one]);
00605 masses[shortlivedDaughters[one]]=aSampler.SampleMass(aDaughter,massmax);
00606 }
00607 break;
00608 default:
00609 theDaughtersName1 = theDecayChannel->GetDaughterName(0);
00610 theDaughtersName2 = theDecayChannel->GetDaughterName(1);
00611 theDaughtersName3 = theDecayChannel->GetDaughterName(2);
00612 if ( numberOfShortliveds == 1)
00613 { G4SampleResonance aSampler;
00614 G4double massmax=theParentMass - SumLongLivedMass;
00615 G4ParticleDefinition * aDaughter=theDecayChannel->GetDaughter(shortlivedDaughters[0]);
00616 masses[shortlivedDaughters[0]]= aSampler.SampleMass(aDaughter,massmax);
00617 }
00618 break;
00619 }
00620
00621
00622
00623
00624
00625 G4GeneralPhaseSpaceDecay thePhaseSpaceDecayChannel(theParentName,
00626 theParentMass,
00627 theBR,
00628 theNumberOfDaughters,
00629 theDaughtersName1,
00630 theDaughtersName2,
00631 theDaughtersName3,
00632 masses);
00633 G4DecayProducts* theDecayProducts = thePhaseSpaceDecayChannel.DecayIt();
00634 if(!theDecayProducts)
00635 {
00636 G4cerr << "Error condition encountered in G4KineticTrack::Decay()"<<G4endl;
00637 G4cerr << " phase-space decay failed."<<G4endl;
00638 G4cerr << " particle was "<<thisDefinition->GetParticleName()<<G4endl;
00639 G4cerr << " channel index "<< chosench << "of "<<nChannels<<"channels"<<G4endl;
00640 G4cerr << " "<<theNumberOfDaughters<< " Daughter particles: "
00641 << theDaughtersName1<<" "<<theDaughtersName2<<" "<<theDaughtersName3<<G4endl;
00642 return 0;
00643 }
00644
00645
00646
00647
00648 G4LorentzRotation toMoving(Get4Momentum().boostVector());
00649 G4DynamicParticle* theDynamicParticle;
00650 G4double formationTime = 0.0;
00651 G4ThreeVector position = this->GetPosition();
00652 G4LorentzVector momentum;
00653 G4LorentzVector momentumBalanceCMS(0);
00654 G4KineticTrackVector* theDecayProductList = new G4KineticTrackVector;
00655 G4int dEntries = theDecayProducts->entries();
00656 G4ParticleDefinition * aProduct = 0;
00657 for (G4int i=dEntries; i > 0; i--)
00658 {
00659 theDynamicParticle = theDecayProducts->PopProducts();
00660 aProduct = theDynamicParticle->GetDefinition();
00661 chargeBalance -= G4lrint(aProduct->GetPDGCharge() );
00662 baryonBalance -= G4lrint(aProduct->GetBaryonNumber() );
00663 momentumBalanceCMS += theDynamicParticle->Get4Momentum();
00664 momentum = toMoving*theDynamicParticle->Get4Momentum();
00665 energyMomentumBalance -= momentum;
00666 theDecayProductList->push_back(new G4KineticTrack (aProduct,
00667 formationTime,
00668 position,
00669 momentum));
00670 delete theDynamicParticle;
00671 }
00672 delete theDecayProducts;
00673 if(getenv("DecayEnergyBalanceCheck"))
00674 std::cout << "DEBUGGING energy balance in cms and lab, charge baryon balance : "
00675 << momentumBalanceCMS << " "
00676 <<energyMomentumBalance << " "
00677 <<chargeBalance<<" "
00678 <<baryonBalance<<" "
00679 <<G4endl;
00680 return theDecayProductList;
00681 }
00682 else
00683 {
00684 return 0;
00685 }
00686 }
00687
00688 G4double G4KineticTrack::IntegrandFunction1(G4double xmass) const
00689 {
00690 G4double mass = theActualMass;
00691 G4double mass1 = theDaughterMass[0];
00692 G4double mass2 = theDaughterMass[1];
00693 G4double gamma2 = theDaughterWidth[1];
00694
00695 G4double result = (1. / (2 * mass)) *
00696 std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
00697 ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
00698 BrWig(gamma2, mass2, xmass);
00699 return result;
00700 }
00701
00702 G4double G4KineticTrack::IntegrandFunction2(G4double xmass) const
00703 {
00704 G4double mass = theDefinition->GetPDGMass();
00705 G4double mass1 = theDaughterMass[0];
00706 G4double mass2 = theDaughterMass[1];
00707 G4double gamma2 = theDaughterWidth[1];
00708 G4double result = (1. / (2 * mass)) *
00709 std::sqrt(std::max((((mass * mass) - (mass1 + xmass) * (mass1 + xmass)) *
00710 ((mass * mass) - (mass1 - xmass) * (mass1 - xmass))),0.0)) *
00711 BrWig(gamma2, mass2, xmass);
00712 return result;
00713 }
00714
00715 G4double G4KineticTrack::IntegrandFunction3(G4double xmass) const
00716 {
00717 const G4double mass = G4KineticTrack_Gmass;
00718
00719 const G4double mass2 = theDaughterMass[1];
00720 const G4double gamma2 = theDaughterWidth[1];
00721
00722 const G4double result = (1. / (2 * mass)) *
00723 std::sqrt(((mass * mass) - (G4KineticTrack_xmass1 + xmass) * (G4KineticTrack_xmass1 + xmass)) *
00724 ((mass * mass) - (G4KineticTrack_xmass1 - xmass) * (G4KineticTrack_xmass1 - xmass))) *
00725 BrWig(gamma2, mass2, xmass);
00726 return result;
00727 }
00728
00729 G4double G4KineticTrack::IntegrandFunction4(G4double xmass) const
00730 {
00731 const G4double mass = G4KineticTrack_Gmass;
00732 const G4double mass1 = theDaughterMass[0];
00733 const G4double gamma1 = theDaughterWidth[0];
00734
00735
00736 G4KineticTrack_xmass1 = xmass;
00737
00738 const G4double theLowerLimit = 0.0;
00739 const G4double theUpperLimit = mass - xmass;
00740 const G4int nIterations = 100;
00741
00742 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
00743 G4double result = BrWig(gamma1, mass1, xmass)*
00744 integral.Simpson(this, &G4KineticTrack::IntegrandFunction3, theLowerLimit, theUpperLimit, nIterations);
00745
00746 return result;
00747 }
00748
00749 G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit) const
00750 {
00751 const G4double theUpperLimit = theActualMass - theDaughterMass[0];
00752 const G4int nIterations = 100;
00753
00754 if (theLowerLimit>=theUpperLimit) return 0.0;
00755
00756 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
00757 G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction1,
00758 theLowerLimit, theUpperLimit, nIterations);
00759 return theIntegralOverMass2;
00760 }
00761
00762 G4double G4KineticTrack::IntegrateCMMomentum(const G4double theLowerLimit, const G4double poleMass) const
00763 {
00764 const G4double theUpperLimit = poleMass - theDaughterMass[0];
00765 const G4int nIterations = 100;
00766
00767 if (theLowerLimit>=theUpperLimit) return 0.0;
00768
00769 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
00770 const G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction2,
00771 theLowerLimit, theUpperLimit, nIterations);
00772 return theIntegralOverMass2;
00773 }
00774
00775
00776 G4double G4KineticTrack::IntegrateCMMomentum2() const
00777 {
00778 const G4double theLowerLimit = 0.0;
00779 const G4double theUpperLimit = theActualMass;
00780 const G4int nIterations = 100;
00781
00782 if (theLowerLimit>=theUpperLimit) return 0.0;
00783
00784 G4Integrator<const G4KineticTrack, G4double(G4KineticTrack::*)(G4double) const> integral;
00785 G4double theIntegralOverMass2 = integral.Simpson(this, &G4KineticTrack::IntegrandFunction4,
00786 theLowerLimit, theUpperLimit, nIterations);
00787 return theIntegralOverMass2;
00788 }
00789
00790
00791
00792
00793
00794
00795
00796