Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // 27 // ------------------------------------------- 28 // GEANT 4 class implementation file 29 // 30 // History: first implementation, A. Feli 31 // ------------------------------------------- 32 33 #include "globals.hh" 34 #include "G4ios.hh" 35 //#include <cmath> 36 37 #include "Randomize.hh" 38 #include "G4SimpleIntegration.hh" 39 #include "G4ThreeVector.hh" 40 #include "G4LorentzVector.hh" 41 #include "G4KineticTrack.hh" 42 #include "G4KineticTrackVector.hh" 43 #include "G4ParticleDefinition.hh" 44 #include "G4DecayTable.hh" 45 #include "G4GeneralPhaseSpaceDecay.hh" 46 #include "G4DecayProducts.hh" 47 #include "G4LorentzRotation.hh" 48 #include "G4SampleResonance.hh" 49 #include "G4Integrator.hh" 50 #include "G4KaonZero.hh" 51 #include "G4KaonZeroShort.hh" 52 #include "G4KaonZeroLong.hh" 53 #include "G4AntiKaonZero.hh" 54 55 #include "G4HadTmpUtil.hh" 56 57 // 58 // Some static clobal for integration 59 // 60 61 static G4ThreadLocal G4double G4KineticTrack_ 62 63 // 64 // Default constructor 65 // 66 67 G4KineticTrack::G4KineticTrack() : 68 theDefinition(0), 69 theFormationTime(0), 70 thePosition(0), 71 the4Momentum(0), 72 theFermi3Momentum(0), 73 theTotal4Momentum(0), 74 theNucleon(0), 75 nChannels(0), 76 theActualMass(0), 77 theActualWidth(0), 78 theDaughterMass(0), 79 theDaughterWidth(0), 80 theStateToNucleus(undefined), 81 theProjectilePotential(0), 82 theCreatorModel(-1), 83 theParentResonanceDef(nullptr) 84 theParentResonanceID(0) 85 { 86 //////////////// 87 // DEBUG // 88 //////////////// 89 90 /* 91 G4cerr << G4endl << G4endl << G4endl; 92 G4cerr << " G4KineticTrack default construc 93 G4cerr << " =============================== 94 */ 95 } 96 97 98 99 // 100 // Copy constructor 101 // 102 103 G4KineticTrack::G4KineticTrack(const G4Kinetic 104 { 105 theDefinition = right.GetDefinition(); 106 theFormationTime = right.GetFormationTime(); 107 thePosition = right.GetPosition(); 108 the4Momentum = right.GetTrackingMomentum(); 109 theFermi3Momentum = right.theFermi3Momentum; 110 theTotal4Momentum = right.theTotal4Momentum; 111 theNucleon=right.theNucleon; 112 nChannels = right.GetnChannels(); 113 theActualMass = right.GetActualMass(); 114 theActualWidth = new G4double[nChannels]; 115 for (G4int i = 0; i < nChannels; i++) 116 { 117 theActualWidth[i] = right.theActualWidth[i 118 } 119 theDaughterMass = 0; 120 theDaughterWidth = 0; 121 theStateToNucleus = right.theStateToNucleus; 122 theProjectilePotential = right.theProjectileP 123 theCreatorModel = right.GetCreatorModelID(); 124 theParentResonanceDef = right.GetParentResona 125 theParentResonanceID = right.GetParentResonan 126 127 //////////////// 128 // DEBUG // 129 //////////////// 130 131 /* 132 G4cerr << G4endl << G4endl << G4endl; 133 G4cerr << " G4KineticTrack copy constructor 134 G4cerr << " =============================== 135 */ 136 } 137 138 139 // 140 // By argument constructor 141 // 142 143 G4KineticTrack::G4KineticTrack(const G4Particl 144 G4double aForma 145 const G4ThreeVe 146 const G4Lorentz 147 theDefinition(aDefinition), 148 theFormationTime(aFormationTime), 149 thePosition(aPosition), 150 the4Momentum(a4Momentum), 151 theFermi3Momentum(0), 152 theTotal4Momentum(a4Momentum), 153 theNucleon(0), 154 theStateToNucleus(undefined), 155 theProjectilePotential(0), 156 theCreatorModel(-1), 157 theParentResonanceDef(nullptr) 158 theParentResonanceID(0) 159 { 160 if(G4KaonZero::KaonZero() == theDefinition | 161 G4AntiKaonZero::AntiKaonZero() == theDefin 162 { 163 if(G4UniformRand()<0.5) 164 { 165 theDefinition = G4KaonZeroShort::KaonZer 166 } 167 else 168 { 169 theDefinition = G4KaonZeroLong::KaonZero 170 } 171 } 172 173 // 174 // Get the number of decay channels 175 // 176 177 G4DecayTable* theDecayTable = theDefinition-> 178 if (theDecayTable != nullptr) 179 { 180 nChannels = theDecayTable->entries(); 181 182 } 183 else 184 { 185 nChannels = 0; 186 } 187 188 // 189 // Get the actual mass value 190 // 191 192 theActualMass = GetActualMass(); 193 194 // 195 // Create an array to Store the actual parti 196 // of the decay channels 197 // 198 199 theDaughterMass = 0; 200 theDaughterWidth = 0; 201 theActualWidth = 0; 202 G4bool * theDaughterIsShortLived = nullptr; 203 204 if(nChannels!=0) theActualWidth = new G4doub 205 206 // cout << " ****CONSTR*** ActualMass ***** 207 G4int index; 208 for (index = nChannels - 1; index >= 0; --in 209 { 210 G4VDecayChannel* theChannel = theDecayTa 211 G4int nDaughters = theChannel->GetNumber 212 G4double theMotherWidth; 213 if (nDaughters == 2 || nDaughters == 3) 214 { 215 G4double thePoleMass = theDefinitio 216 theMotherWidth = theDefinition->GetP 217 G4double thePoleWidth = theChannel-> 218 const G4ParticleDefinition* aDaughte 219 theDaughterMass = new G4double[nDaug 220 theDaughterWidth = new G4double[nDau 221 theDaughterIsShortLived = new G4bool[nDaug 222 for (G4int n = 0; n < nDaughters; ++ 223 { 224 aDaughter = theChannel->GetDaugh 225 theDaughterMass[n] = aDaughter-> 226 theDaughterWidth[n] = aDaughter- 227 theDaughterIsShortLived[n] = aDaughter 228 } 229 230 // 231 // Check whether both the decay prod 232 // 233 234 G4double theActualMom = 0.0; 235 G4double thePoleMom = 0.0; 236 G4SampleResonance aSampler; 237 if (nDaughters==2) 238 { 239 if ( !theDaughterIsShortLived[0] && !t 240 { 241 242 // G4cout << G4endl << "Bot 243 // " decay 244 // cout << " LB: Both de 245 // cout << " parent: 246 // cout << " particle1: 247 // cout << " particle2: 248 249 theActualMom = EvaluateCMMomentum(theAct 250 theDaughterMass); 251 thePoleMom = EvaluateCMMomentum(thePoleM 252 theDaughterMass); 253 // cout << G4endl; 254 // cout << " LB: ActualMass/Daugh 255 // cout << " LB: ActualMom " << t 256 // cout << " LB: PoleMom " << t 257 // cout << G4endl; 258 } 259 else if ( !theDaughterIsShortLived[0] 260 { 261 262 // G4cout << G4endl << "Onl 263 // cout << " LB: only the first 264 // cout << " parent: " << th 265 // cout << " particle1: " << th 266 // cout << " particle2: " << th 267 268 // global variable definition 269 G4double lowerLimit = aSampler.GetMinimu 270 theActualMom = IntegrateCMMomentum(lower 271 thePoleMom = IntegrateCMMomentum(lowerLi 272 // cout << " LB Parent Mass = " < 273 // cout << " LB Actual Mass = " < 274 // cout << " LB Daughter1 Mass = 275 // cout << " LB Daughter2 Mass = 276 // cout << " The Actual Momentum 277 // cout << " The Pole Momentum 278 // cout << G4endl; 279 280 } 281 else if ( theDaughterIsShortLived[0] & 282 { 283 284 // G4cout << G4endl << "Onl 285 // " decay 286 // cout << " LB: only th 287 // cout << " parent: " << th 288 // cout << " particle1: " << th 289 // cout << " particle2: " << th 290 291 // 292 // Swap the content of the theDa 293 // 294 295 G4SwapObj(theDaughterMass, theDaughterMa 296 G4SwapObj(theDaughterWidth, theDaughterW 297 298 // global variable definition 299 G4double lowerLimit = aSampler.GetMinimu 300 theActualMom = IntegrateCMMomentum(lower 301 thePoleMom = IntegrateCMMomentum(lowerLi 302 // cout << " LB Parent Mass = " < 303 // cout << " LB Actual Mass = " < 304 // cout << " LB Daughter1 Mass = 305 // cout << " LB Daughter2 Mass = 306 // cout << " The Actual Momentum 307 // cout << " The Pole Momentum 308 // cout << G4endl; 309 310 } 311 else if ( theDaughterIsShortLived[0] & 312 { 313 314 // G4cout << G4endl << "Both the 315 // " decay produc 316 // cout << " LB: both decay prod 317 // cout << " parent: " << th 318 // cout << " particle1: " << th 319 // cout << " particle2: " << th 320 321 // global variable definition 322 G4KineticTrack_Gmass = theActualMass; 323 theActualMom = IntegrateCMMomentum2(); 324 G4KineticTrack_Gmass = thePoleMass; 325 thePoleMom = IntegrateCMMomentum2(); 326 // cout << " LB Parent Mass = " < 327 // cout << " LB Daughter1 Mass = 328 // cout << " LB Daughter2 Mass = 329 // cout << " The Actual Mom 330 // cout << " The Pole Momen 331 // cout << G4endl; 332 333 } 334 } 335 else // (nDaughter==3) 336 { 337 338 G4int nShortLived = 0; 339 if ( theDaughterIsShortLived[0] ) 340 { 341 ++nShortLived; 342 } 343 if ( theDaughterIsShortLived[1] ) 344 { 345 ++nShortLived; 346 G4SwapObj(theDaughterMass, theDaughterMa 347 G4SwapObj(theDaughterWidth, theDaughterW 348 } 349 if ( theDaughterIsShortLived[2] ) 350 { 351 ++nShortLived; 352 G4SwapObj(theDaughterMass, theDaughterMa 353 G4SwapObj(theDaughterWidth, theDaughterW 354 } 355 if ( nShortLived == 0 ) 356 { 357 theDaughterMass[1]+=theDaughterMass[2]; 358 theActualMom = EvaluateCMMomentum(theAct 359 theDaughterMass); 360 thePoleMom = EvaluateCMMomentum(thePoleM 361 theDaughterMass); 362 } 363 // else if ( nShortLived == 1 ) 364 else if ( nShortLived >= 1 ) 365 { 366 // need the shortlived particle in slot 367 G4SwapObj(theDaughterMass, theDaughterMa 368 G4SwapObj(theDaughterWidth, theDaughterW 369 theDaughterMass[0] += theDaughterMass[2] 370 theActualMom = IntegrateCMMomentum(0.0); 371 thePoleMom = IntegrateCMMomentum(0.0, th 372 } 373 // else 374 // { 375 // throw G4HadronicException(__FILE__, __ 376 // } 377 378 } 379 380 //if(nDaughters<3) theChannel->GetAngularM 381 G4double theMassRatio = thePoleMass / theA 382 G4double theMomRatio = theActualMom 383 // VI 11.06.2015: for l=0 one not need use 384 //G4double l=0; 385 //theActualWidth[index] = thePoleWid 386 // std::pow(t 387 // (1.2 / (1+ 388 theActualWidth[index] = thePoleWidth 389 theMomRatio; 390 delete [] theDaughterMass; 391 theDaughterMass = nullptr; 392 delete [] theDaughterWidth; 393 theDaughterWidth = nullptr; 394 delete [] theDaughterIsShortLived; 395 theDaughterIsShortLived = nullptr; 396 } 397 398 else // nDaughter = 1 ( e.g. K0 decays 399 { 400 theMotherWidth = theDefinition->GetPDGWidt 401 theActualWidth[index] = theChannel->GetBR( 402 } 403 } 404 405 //////////////// 406 // DEBUG // 407 //////////////// 408 409 // for (G4int y = nChannels - 1; y >= 0; --y) 410 // { 411 // G4cout << G4endl << theActualWidth[y]; 412 // } 413 // G4cout << G4endl << G4endl << G4endl; 414 415 /* 416 G4cerr << G4endl << G4endl << G4endl; 417 G4cerr << " G4KineticTrack by argument cons 418 G4cerr << " =============================== 419 */ 420 421 } 422 423 G4KineticTrack::G4KineticTrack(G4Nucleon * nuc 424 const G4ThreeVector& aPosition, 425 const G4Lorentz 426 : theDefinition(nucleon->GetDefinition() 427 theFormationTime(0), 428 thePosition(aPosition), 429 the4Momentum(a4Momentum), 430 theFermi3Momentum(nucleon->GetMomentum()), 431 theNucleon(nucleon), 432 nChannels(0), 433 theActualMass(nucleon->GetDefinition()->GetP 434 theActualWidth(0), 435 theDaughterMass(0), 436 theDaughterWidth(0), 437 theStateToNucleus(undefined), 438 theProjectilePotential(0), 439 theCreatorModel(-1), 440 theParentResonanceDef(nullptr), 441 theParentResonanceID(0) 442 { 443 theFermi3Momentum.setE(0); 444 Set4Momentum(a4Momentum); 445 } 446 447 448 G4KineticTrack::~G4KineticTrack() 449 { 450 if (theActualWidth != 0) delete [] theActualW 451 if (theDaughterMass != 0) delete [] theDaught 452 if (theDaughterWidth != 0) delete [] theDaugh 453 } 454 455 456 457 G4KineticTrack& G4KineticTrack::operator=(cons 458 { 459 if (this != &right) 460 { 461 theDefinition = right.GetDefinition(); 462 theFormationTime = right.GetFormationTime 463 the4Momentum = right.the4Momentum; 464 the4Momentum = right.GetTrackingMomentum( 465 theFermi3Momentum = right.theFermi3Moment 466 theTotal4Momentum = right.theTotal4Moment 467 theNucleon = right.theNucleon; 468 theStateToNucleus = right.theStateToNucle 469 delete [] theActualWidth; 470 nChannels = right.GetnChannels(); 471 theActualWidth = new G4double[nChannels]; 472 for (G4int i = 0; i < nChannels; ++i) the 473 theCreatorModel = right.GetCreatorModelID 474 theParentResonanceDef = right.GetParentRe 475 theParentResonanceID = right.GetParentRes 476 } 477 return *this; 478 } 479 480 481 482 G4bool G4KineticTrack::operator==(const G4Kine 483 { 484 return (this == & right); 485 } 486 487 488 489 G4bool G4KineticTrack::operator!=(const G4Kine 490 { 491 return (this != & right); 492 } 493 494 495 496 G4KineticTrackVector* G4KineticTrack::Decay() 497 { 498 // 499 // Select a possible decay channel 500 // 501 /* 502 G4int index1; 503 for (index1 = nChannels - 1; index1 >= 0; 504 G4cout << "DECAY Actual Width IND/Actual 505 G4cout << "DECAY Actual Mass " << theAct 506 */ 507 const G4ParticleDefinition* thisDefinition = 508 if(!thisDefinition) 509 { 510 G4cerr << "Error condition encountered in 511 G4cerr << " track has no particle definit 512 return nullptr; 513 } 514 G4DecayTable* theDecayTable = thisDefinition 515 if(!theDecayTable) 516 { 517 G4cerr << "Error condition encountered in 518 G4cerr << " particle definition has no de 519 G4cerr << " particle was "<<thisDefinitio 520 return nullptr; 521 } 522 523 G4int chargeBalance = G4lrint(theDefinition-> 524 G4int baryonBalance = G4lrint(theDefinition-> 525 G4LorentzVector energyMomentumBalance(Get4Mom 526 G4double theTotalActualWidth = this->Evaluate 527 if (theTotalActualWidth !=0) 528 { 529 530 //AR-16Aug2016 : Repeat the sampling of t 531 // kinematically above thre 532 G4bool isChannelBelowThreshold = true; 533 const G4int maxNumberOfLoops = 10000; 534 G4int loopCounter = 0; 535 536 G4int chosench; 537 G4String theParentName; 538 G4double theParentMass; 539 G4double theBR; 540 G4int theNumberOfDaughters; 541 G4String theDaughtersName1; 542 G4String theDaughtersName2; 543 G4String theDaughtersName3; 544 G4String theDaughtersName4; 545 G4double masses[4]={0.,0.,0.,0.}; 546 547 do { 548 549 G4double theSumActualWidth = 0.0; 550 G4double* theCumActualWidth = new G4dou 551 for (G4int index = nChannels - 1; index 552 { 553 theSumActualWidth += theActualWidth 554 theCumActualWidth[index] = theSumAc 555 // cout << "DECAY Cum. Width " << index 556 } 557 // cout << "DECAY Total Width " << the 558 // cout << "DECAY Total Width " << the 559 G4double r = theTotalActualWidth * G4Un 560 G4VDecayChannel* theDecayChannel(nullpt 561 chosench=-1; 562 for (G4int index = nChannels - 1; index 563 { 564 if (r < theCumActualWidth[index]) 565 { 566 theDecayChannel = theDecayTable 567 // cout << "DECAY SELECTED CHANN 568 chosench=index; 569 break; 570 } 571 } 572 573 delete [] theCumActualWidth; 574 575 if(theDecayChannel == nullptr) 576 { 577 G4cerr << "Error condition encountere 578 G4cerr << " decay channel has 0x0 ch 579 G4cerr << " particle was "<<thisDefi 580 G4cerr << " channel index "<< chosen 581 return nullptr; 582 } 583 theParentName = theDecayChannel->GetPar 584 theParentMass = this->GetActualMass(); 585 theBR = theActualWidth[chosench]; 586 // cout << "**BR*** DECAYNEW " << 587 theNumberOfDaughters = theDecayChannel- 588 theDaughtersName1 = ""; 589 theDaughtersName2 = ""; 590 theDaughtersName3 = ""; 591 theDaughtersName4 = ""; 592 593 for (G4int i=0; i < 4; ++i) masses[i]=0 594 G4int shortlivedDaughters[4]; 595 G4int numberOfShortliveds(0); 596 G4double SumLongLivedMass(0); 597 for (G4int aD=0; aD < theNumberOfDaught 598 { 599 const G4ParticleDefinition* aDaughte 600 masses[aD] = aDaughter->GetPDGMass() 601 if ( aDaughter->IsShortLived() ) 602 { 603 shortlivedDaughters[numberOfShortlived 604 ++numberOfShortliveds; 605 } else { 606 SumLongLivedMass += aDaughter->GetPDGM 607 } 608 609 } 610 switch (theNumberOfDaughters) 611 { 612 case 0: 613 break; 614 case 1: 615 theDaughtersName1 = theDecayChan 616 theDaughtersName2 = ""; 617 theDaughtersName3 = ""; 618 theDaughtersName4 = ""; 619 break; 620 case 2: 621 theDaughtersName1 = theDecayChan 622 theDaughtersName2 = theDecayChan 623 theDaughtersName3 = ""; 624 theDaughtersName4 = ""; 625 if ( numberOfShortliveds == 1) 626 { G4SampleResonance aSampler; 627 G4double massmax=theParentMa 628 const G4ParticleDefinition * aDaughter=t 629 masses[shortlivedDaughters[0]]= aS 630 } else if ( numberOfShortliveds == 2) 631 // choose masses one after the oth 632 G4int zero= (G4UniformRand() > 0.5 633 G4int one = 1-zero; 634 G4SampleResonance aSampler; 635 G4double massmax=theParentMass - aSample 636 const G4ParticleDefinition * aDaughter=t 637 masses[shortlivedDaughters[zero]]=aSampl 638 massmax=theParentMass - masses[shortlive 639 aDaughter=theDecayChannel->GetDaughter(s 640 masses[shortlivedDaughters[one]]=aSample 641 } 642 break; 643 case 3: 644 theDaughtersName1 = theDecayChan 645 theDaughtersName2 = theDecayChan 646 theDaughtersName3 = theDecayChan 647 theDaughtersName4 = ""; 648 if ( numberOfShortliveds == 1) 649 { G4SampleResonance aSampler; 650 G4double massmax=theParentMa 651 const G4ParticleDefinition * aDaught 652 masses[shortlivedDaughters[0]]= aS 653 } 654 break; 655 default: 656 theDaughtersName1 = theDecayChan 657 theDaughtersName2 = theDecayChan 658 theDaughtersName3 = theDecayChan 659 theDaughtersName4 = theDecayChan 660 if ( numberOfShortliveds == 1) 661 { G4SampleResonance aSampler; 662 G4double massmax=theParentMa 663 const G4ParticleDefinition * aDaught 664 masses[shortlivedDaughters[0]]= aS 665 } 666 if ( theNumberOfDaughters > 4 ) 667 G4ExceptionDescription ed; 668 ed << "More than 4 decay daugh 669 G4Exception( "G4KineticTrack:: 670 } 671 break; 672 } 673 674 //AR-16Aug2016 : Check whether the s 675 // If this is still no 676 // then the subsequent 677 G4double sumDaughterMasses = 0.0; 678 for (G4int i=0; i < 4; ++i) sumDaugh 679 if ( theParentMass - sumDaughterMass 680 681 } while ( isChannelBelowThreshold && ++lo 682 683 // 684 // Get the decay products List 685 // 686 687 G4GeneralPhaseSpaceDecay thePhaseSpaceDec 688 689 690 691 692 th 693 th 694 th 695 masses); 696 G4DecayProducts* theDecayProducts = thePh 697 if(theDecayProducts == nullptr) 698 { 699 G4ExceptionDescription ed; 700 ed << "Error condition encountered: pha 701 << "\t the decaying particle is: " < 702 << "\t the channel index is: "<< cho 703 << "\t " << theNumberOfDaughters << 704 << theDaughtersName1 << " " << theDa 705 << theDaughtersName4 << G4endl; 706 G4Exception( "G4KineticTrack::Decay ", 707 return nullptr; 708 } 709 710 // 711 // Create the kinetic track List associat 712 // 713 // For the decay products of hadronic res 714 // the same as their parent 715 G4LorentzRotation toMoving(Get4Momentum() 716 G4DynamicParticle* theDynamicParticle; 717 G4double formationTime = 0.0; 718 G4ThreeVector position = this->GetPositio 719 G4LorentzVector momentum; 720 G4LorentzVector momentumBalanceCMS(0); 721 G4KineticTrackVector* theDecayProductList 722 G4int dEntries = theDecayProducts->entrie 723 const G4ParticleDefinition * aProduct = n 724 // Use the integer round mass in keV to g 725 G4int uniqueID = static_cast< G4int >( ro 726 for (G4int i=dEntries; i > 0; --i) 727 { 728 theDynamicParticle = theDecayProducts->PopP 729 aProduct = theDynamicParticle->GetDef 730 chargeBalance -= G4lrint(aProduct->Ge 731 baryonBalance -= G4lrint(aProduct->Ge 732 momentumBalanceCMS += theDynamicParticle->G 733 momentum = toMoving*theDynamicParticl 734 energyMomentumBalance -= momentum; 735 G4KineticTrack* aDaughter = new G4Kin 736 737 738 739 if (aDaughter != nullptr) 740 { 741 aDaughter->SetCreatorModelID(GetC 742 aDaughter->SetParentResonanceDef( 743 aDaughter->SetParentResonanceID(u 744 } 745 theDecayProductList->push_back(aDaugh 746 delete theDynamicParticle; 747 } 748 delete theDecayProducts; 749 if(std::getenv("DecayEnergyBalanceCheck") 750 std::cout << "DEBUGGING energy balance 751 << momentumBalanceCMS << " " 752 <<energyMomentumBalance << " " 753 <<chargeBalance<<" " 754 <<baryonBalance<<" " 755 <<G4endl; 756 return theDecayProductList; 757 } 758 else 759 { 760 return nullptr; 761 } 762 } 763 764 G4double G4KineticTrack::IntegrandFunction1(G4 765 { 766 G4double mass = theActualMass; /* the actu 767 G4double mass1 = theDaughterMass[0]; 768 G4double mass2 = theDaughterMass[1]; 769 G4double gamma2 = theDaughterWidth[1]; 770 771 G4double result = (1. / (2 * mass)) * 772 std::sqrt(std::max((((mass * mass) - (mass 773 ((mass * mass) - (mass1 - xmass) * (mas 774 BrWig(gamma2, mass2, xmass); 775 return result; 776 } 777 778 G4double G4KineticTrack::IntegrandFunction2(G4 779 { 780 G4double mass = theDefinition->GetPDGMass(); 781 G4double mass1 = theDaughterMass[0]; 782 G4double mass2 = theDaughterMass[1]; 783 G4double gamma2 = theDaughterWidth[1]; 784 G4double result = (1. / (2 * mass)) * 785 std::sqrt(std::max((((mass * mass) - (mass 786 ((mass * mass) - (mass1 - xmass) * (mas 787 BrWig(gamma2, mass2, xmass); 788 return result; 789 } 790 791 G4double G4KineticTrack::IntegrandFunction3(G4 792 { 793 const G4double mass = G4KineticTrack_Gmass; 794 // const G4double mass1 = theDaughterMass[0]; 795 const G4double mass2 = theDaughterMass[1]; 796 const G4double gamma2 = theDaughterWidth[1]; 797 798 const G4double result = (1. / (2 * mass)) * 799 std::sqrt(((mass * mass) - (G4KineticTrack 800 ((mass * mass) - (G4KineticTrack_xmass1 - x 801 BrWig(gamma2, mass2, xmass); 802 return result; 803 } 804 805 G4double G4KineticTrack::IntegrandFunction4(G4 806 { 807 const G4double mass = G4KineticTrack_Gmass; 808 const G4double mass1 = theDaughterMass[0]; 809 const G4double gamma1 = theDaughterWidth[0]; 810 // const G4double mass2 = theDaughterMass[1]; 811 812 G4KineticTrack_xmass1 = xmass; 813 814 const G4double theLowerLimit = 0.0; 815 const G4double theUpperLimit = mass - xmass; 816 const G4int nIterations = 100; 817 818 G4Integrator<const G4KineticTrack, G4double( 819 G4double result = BrWig(gamma1, mass1, xmass 820 integral.Simpson(this, &G4KineticTrack::In 821 822 return result; 823 } 824 825 G4double G4KineticTrack::IntegrateCMMomentum(c 826 { 827 const G4double theUpperLimit = theActualMass 828 const G4int nIterations = 100; 829 830 if (theLowerLimit>=theUpperLimit) return 0.0 831 832 G4Integrator<const G4KineticTrack, G4double( 833 G4double theIntegralOverMass2 = integral.Sim 834 theLowerLimit, theUpperLimit, n 835 return theIntegralOverMass2; 836 } 837 838 G4double G4KineticTrack::IntegrateCMMomentum(c 839 { 840 const G4double theUpperLimit = poleMass - th 841 const G4int nIterations = 100; 842 843 if (theLowerLimit>=theUpperLimit) return 0.0 844 845 G4Integrator<const G4KineticTrack, G4double( 846 const G4double theIntegralOverMass2 = integr 847 theLowerLimit, theUpperLimit, 848 return theIntegralOverMass2; 849 } 850 851 852 G4double G4KineticTrack::IntegrateCMMomentum2( 853 { 854 const G4double theLowerLimit = 0.0; 855 const G4double theUpperLimit = theActualMass 856 const G4int nIterations = 100; 857 858 if (theLowerLimit>=theUpperLimit) return 0.0 859 860 G4Integrator<const G4KineticTrack, G4double( 861 G4double theIntegralOverMass2 = integral.Sim 862 theLowerLimit, theUpperLimit, n 863 return theIntegralOverMass2; 864 } 865