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 file 29 // 30 // History: first implementation 31 // HPW, 10DEC 98, the decay part original 32 // in his FTF-test-program. 33 // 34 // M.Kelsey, 28 Jul 2011 -- Replace loop 35 // with new utility class, simplify cleanup 36 // 37 // A.Ribon, 27 Oct 2021 -- Extended the m 38 // to deal with projectile hyper 39 // 40 // ------------------------------------------- 41 42 #include <algorithm> 43 #include <vector> 44 45 #include "G4GeneratorPrecompoundInterface.hh" 46 #include "G4PhysicalConstants.hh" 47 #include "G4SystemOfUnits.hh" 48 #include "G4DynamicParticleVector.hh" 49 #include "G4KineticTrackVector.hh" 50 #include "G4Proton.hh" 51 #include "G4Neutron.hh" 52 #include "G4Lambda.hh" 53 54 #include "G4Deuteron.hh" 55 #include "G4Triton.hh" 56 #include "G4He3.hh" 57 #include "G4Alpha.hh" 58 59 #include "G4V3DNucleus.hh" 60 #include "G4Nucleon.hh" 61 62 #include "G4AntiProton.hh" 63 #include "G4AntiNeutron.hh" 64 #include "G4AntiLambda.hh" 65 #include "G4AntiDeuteron.hh" 66 #include "G4AntiTriton.hh" 67 #include "G4AntiHe3.hh" 68 #include "G4AntiAlpha.hh" 69 70 #include "G4HyperTriton.hh" 71 #include "G4HyperH4.hh" 72 #include "G4HyperAlpha.hh" 73 #include "G4HyperHe5.hh" 74 #include "G4DoubleHyperH4.hh" 75 #include "G4DoubleHyperDoubleNeutron.hh" 76 77 #include "G4AntiHyperTriton.hh" 78 #include "G4AntiHyperH4.hh" 79 #include "G4AntiHyperAlpha.hh" 80 #include "G4AntiHyperHe5.hh" 81 #include "G4AntiDoubleHyperH4.hh" 82 #include "G4AntiDoubleHyperDoubleNeutron.hh" 83 84 #include "G4FragmentVector.hh" 85 #include "G4ReactionProduct.hh" 86 #include "G4ReactionProductVector.hh" 87 #include "G4PreCompoundModel.hh" 88 #include "G4ExcitationHandler.hh" 89 #include "G4DecayKineticTracks.hh" 90 #include "G4HadronicInteractionRegistry.hh" 91 92 #include "G4PhysicsModelCatalog.hh" 93 #include "G4HyperNucleiProperties.hh" 94 //-------------------------------------------- 95 #include "Randomize.hh" 96 #include "G4Log.hh" 97 98 //#define debugPrecoInt 99 100 G4GeneratorPrecompoundInterface::G4GeneratorPr 101 : CaptureThreshold(70*MeV), DeltaM(5.0*MeV), 102 { 103 proton = G4Proton::Proton(); 104 neutron = G4Neutron::Neutron(); 105 lambda = G4Lambda::Lambda(); 106 107 deuteron=G4Deuteron::Deuteron(); 108 triton =G4Triton::Triton(); 109 He3 =G4He3::He3(); 110 He4 =G4Alpha::Alpha(); 111 112 ANTIproton=G4AntiProton::AntiProton(); 113 ANTIneutron=G4AntiNeutron::AntiNeutron(); 114 115 ANTIdeuteron=G4AntiDeuteron::AntiDeuteron() 116 ANTItriton =G4AntiTriton::AntiTriton(); 117 ANTIHe3 =G4AntiHe3::AntiHe3(); 118 ANTIHe4 =G4AntiAlpha::AntiAlpha(); 119 120 if(preModel) { SetDeExcitation(preModel); } 121 else { 122 G4HadronicInteraction* hadi = 123 G4HadronicInteractionRegistry::Ins 124 G4VPreCompoundModel* pre = static_cast<G 125 if(!pre) { pre = new G4PreCompoundModel( 126 SetDeExcitation(pre); 127 } 128 129 secID = G4PhysicsModelCatalog::GetModelID(" 130 } 131 132 G4GeneratorPrecompoundInterface::~G4GeneratorP 133 { 134 } 135 136 G4ReactionProductVector* G4GeneratorPrecompoun 137 Propagate(G4KineticTrackVector* theSecondaries 138 { 139 #ifdef debugPrecoInt 140 G4cout<<G4endl<<"G4GeneratorPrecompoundI 141 G4cout<<"Target A and Z "<<theNucleus->G 142 G4cout<<"Directly produced particles num 143 #endif 144 145 G4ReactionProductVector * theTotalResult = 146 147 // decay the strong resonances 148 G4DecayKineticTracks decay(theSecondaries); 149 #ifdef debugPrecoInt 150 G4cout<<"Final stable particles number " 151 #endif 152 153 // prepare the fragment (it is assumed that 154 G4int anA=theNucleus->GetMassNumber(); 155 G4int aZ=theNucleus->GetCharge(); 156 // G4double TargetNucleusMass = G4NucleiPrope 157 158 G4int numberOfEx = 0; 159 G4int numberOfCh = 0; 160 G4int numberOfHoles = 0; 161 162 G4double R = theNucleus->GetNuclearRadius() 163 164 G4LorentzVector captured4Momentum(0.,0.,0., 165 G4LorentzVector Residual4Momentum(0.,0.,0., 166 G4LorentzVector Secondary4Momentum(0.,0.,0. 167 168 // loop over secondaries 169 G4KineticTrackVector::iterator iter; 170 for(iter=theSecondaries->begin(); iter !=th 171 { 172 const G4ParticleDefinition* part = (*ite 173 G4double e = (*iter)->Get4Momentum().e() 174 G4double mass = (*iter)->Get4Momentum(). 175 G4ThreeVector mom = (*iter)->Get4Momentu 176 if((part != proton && part != neutron) | 177 ((*iter)->GetPosition().mag() > R) 178 G4ReactionProduct * theNew = new G4Re 179 theNew->SetMomentum(mom); 180 theNew->SetTotalEnergy(e); 181 theNew->SetCreatorModelID((*iter)->GetCreat 182 theNew->SetParentResonanceDef((*iter)->GetP 183 theNew->SetParentResonanceID((*iter)->GetPa 184 theTotalResult->push_back(theNew); 185 Secondary4Momentum += (*iter)->Get4Mo 186 #ifdef debugPrecoInt 187 G4cout<<"Secondary 4Mom "<<part->G 188 <<(*iter)->Get4Momentum().ma 189 #endif 190 } else { 191 if( e-mass > -CaptureThreshold*G4Log( 192 G4ReactionProduct * theNew = new G 193 theNew->SetMomentum(mom); 194 theNew->SetTotalEnergy(e); 195 theNew->SetCreatorModelID((*iter)->GetCr 196 theNew->SetParentResonanceDef((*iter)->G 197 theNew->SetParentResonanceID((*iter)->Ge 198 theTotalResult->push_back(theNew); 199 Secondary4Momentum += (*iter)->Get 200 #ifdef debugPrecoInt 201 G4cout<<"Secondary 4Mom "<<part 202 <<(*iter)->Get4Momentum() 203 #endif 204 } else { 205 // within the nucleus, neutron or 206 // now calculate A, Z of the frag 207 ++anA; 208 ++numberOfEx; 209 G4int Z = G4int(part->GetPDGCharge 210 aZ += Z; 211 numberOfCh += Z; 212 captured4Momentum += (*iter)->Get4 213 #ifdef debugPrecoInt 214 G4cout<<"Captured 4Mom "<<part 215 #endif 216 } 217 } 218 delete (*iter); 219 } 220 delete theSecondaries; 221 222 // loop over wounded nucleus 223 G4Nucleon * theCurrentNucleon = 224 theNucleus->StartLoop() ? theNucleus- 225 while(theCurrentNucleon) /* Loop checkin 226 { 227 if(theCurrentNucleon->AreYouHit()) { 228 ++numberOfHoles; 229 ++numberOfEx; 230 --anA; 231 aZ -= G4int(theCurrentNucleon->GetDef 232 233 Residual4Momentum -= theCurrentNucleo 234 } 235 theCurrentNucleon = theNucleus->GetNextN 236 } 237 238 #ifdef debugPrecoInt 239 G4cout<<G4endl; 240 G4cout<<"Secondary 4Mom "<<Secondary4Mom 241 G4cout<<"Captured 4Mom "<<captured4Mome 242 G4cout<<"Sec + Captured "<<Secondary4Mom 243 G4cout<<"Residual4Mom "<<Residual4Mome 244 G4cout<<"Sum 4 momenta " 245 <<Secondary4Momentum + captured4Mo 246 #endif 247 248 // Check that we use QGS model; loop over w 249 G4bool QGSM(false); 250 theCurrentNucleon = theNucleus->StartLoop() 251 while(theCurrentNucleon) /* Loop checking 252 { 253 if(theCurrentNucleon->AreYouHit()) 254 { 255 if(theCurrentNucleon->Get4Momentum().ma 256 theCurrentNucleon->GetDefinition()-> 257 } 258 theCurrentNucleon = theNucleus->GetNextN 259 } 260 261 #ifdef debugPrecoInt 262 if(!QGSM){ 263 G4cout<<G4endl; 264 G4cout<<"Residual A and Z "<<anA<<" "< 265 G4cout<<"Residual 4Mom "<<Residual4Mo 266 if(numberOfEx == 0) 267 {G4cout<<"Residual 4Mom = 0 means tha 268 } 269 #endif 270 271 if(anA == 0) return theTotalResult; 272 273 G4LorentzVector exciton4Momentum(0.,0.,0.,0 274 if(anA >= aZ) 275 { 276 if(!QGSM) 277 { // FTF model was used 278 G4double fMass = G4NucleiProperties::Ge 279 280 // G4LorentzVector exciton4Momentum = Res 281 exciton4Momentum = Residual4Momentum + c 282 //exciton4Momentum.setE(std::sqrt(exciton4Mome 283 G4double ActualMass = exciton4Momentum.m 284 if(ActualMass <= fMass ) { 285 exciton4Momentum.setE(std::sqrt(excito 286 } 287 288 #ifdef debugPrecoInt 289 G4double exEnergy = 0.0; 290 if(ActualMass <= fMass ) {exEnergy = 291 else {exEnergy = 292 G4cout<<"Ground state residual Mass " 293 #endif 294 } 295 else 296 { // QGS model was used 297 G4double InitialTargetMass = 298 G4NucleiProperties::GetNuclearMa 299 300 exciton4Momentum = 301 GetPrimaryProjectile()->Get4Mome 302 -Secondary4Momentum; 303 304 G4double fMass = G4NucleiProperties::Get 305 G4double ActualMass = exciton4Momentum.ma 306 307 #ifdef debugPrecoInt 308 G4cout<<G4endl; 309 G4cout<<"Residual A and Z "<<anA<<" "< 310 G4cout<<"Residual4Momentum "<<exciton4 311 G4cout<<"ResidualMass, GroundStateMass 312 <<ActualMass - fMass<<G4endl; 313 #endif 314 315 if(ActualMass - fMass < 0.) 316 { 317 G4double ResE = std::sqrt(exciton4Moment 318 exciton4Momentum.setE(ResE); 319 #ifdef debugPrecoInt 320 G4cout<<"ActualMass - fMass < 0. "<<A 321 #endif 322 } 323 } 324 325 // Need to de-excite the remnant nucleus o 326 G4Fragment anInitialState(anA, aZ, exciton 327 anInitialState.SetNumberOfParticles(number 328 anInitialState.SetNumberOfCharged(numberOf 329 anInitialState.SetNumberOfHoles(numberOfHo 330 anInitialState.SetCreatorModelID(secID); 331 332 G4ReactionProductVector * aPrecoResult = 333 theDeExcitation->DeExcite(anInitialState 334 // fill pre-compound part into the result, 335 #ifdef debugPrecoInt 336 G4cout<<"Target fragment number "<<aPrecoR 337 #endif 338 for(unsigned int ll=0; ll<aPrecoResult->si 339 { 340 theTotalResult->push_back(aPrecoResult-> 341 #ifdef debugPrecoInt 342 G4cout<<"Fragment "<<ll<<" " 343 <<aPrecoResult->operator[](ll)->GetDefin 344 <<aPrecoResult->operator[](ll)->GetMomen 345 <<aPrecoResult->operator[](ll)->GetTotal 346 <<aPrecoResult->operator[](ll)->GetDefin 347 #endif 348 } 349 delete aPrecoResult; 350 } 351 352 return theTotalResult; 353 } 354 355 G4HadFinalState* G4GeneratorPrecompoundInterfa 356 ApplyYourself(const G4HadProjectile &, G4Nucle 357 { 358 G4cout << "G4GeneratorPrecompoundInterface: 359 << G4endl; 360 G4cout << "This class is only a mediator be 361 G4cout << "Please remove from your physics 362 throw G4HadronicException(__FILE__, __LINE_ 363 return new G4HadFinalState; 364 } 365 void G4GeneratorPrecompoundInterface::Propagat 366 { 367 outFile << "G4GeneratorPrecompoundInterface 368 << "energy model through the wounded 369 << "Low energy protons and neutron pr 370 << "the high energy generator and wit 371 << "nucleus and the captured particle 372 << "fragment is passed to the Geant4 373 << "Nuclear de-excitation:\n"; 374 // preco 375 376 } 377 378 379 G4ReactionProductVector* G4GeneratorPrecompoun 380 PropagateNuclNucl(G4KineticTrackVector* theSec 381 G4V3DNucleus* theProjectileNucleus) 382 { 383 #ifdef debugPrecoInt 384 G4cout<<G4endl<<"G4GeneratorPrecompoundInte 385 G4cout<<"Projectile A and Z (and numberOfLa 386 <<theProjectileNucleus->GetCharge()<<" (" 387 <<theProjectileNucleus->GetNumberOfLambdas( 388 G4cout<<"Target A and Z "<<theNucleus-> 389 <<theNucleus->GetCharge()<<" (" 390 <<theNucleus->GetNumberOfLambdas()<<")"<<G4 391 G4cout<<"Directly produced particles number 392 G4cout<<"Projectile 4Mom and mass "<<GetPri 393 <<GetPri 394 #endif 395 396 // prepare the target residual (assumed to 397 G4int anA=theNucleus->GetMassNumber(); 398 G4int aZ=theNucleus->GetCharge(); 399 //G4int aL=theNucleus->GetNumberOfLambdas() 400 G4int numberOfEx = 0; 401 G4int numberOfCh = 0; 402 G4int numberOfHoles = 0; 403 G4double exEnergy = 0.0; 404 G4double R = theNucleus->GetNuclearRadius() 405 G4LorentzVector Target4Momentum(0.,0.,0.,0. 406 407 // loop over the wounded target nucleus 408 G4Nucleon * theCurrentNucleon = 409 theNucleus->StartLoop() ? theNucleus- 410 while(theCurrentNucleon) /* Loop checking 411 { 412 if(theCurrentNucleon->AreYouHit()) { 413 ++numberOfHoles; 414 ++numberOfEx; 415 --anA; 416 aZ -= G4int(theCurrentNucleon->GetDef 417 eplus + 0.1); 418 exEnergy += theCurrentNucleon->GetBin 419 Target4Momentum -=theCurrentNucleon-> 420 } 421 theCurrentNucleon = theNucleus->GetNextN 422 } 423 424 #ifdef debugPrecoInt 425 G4cout<<"Residual Target A Z (numberOfLambd 426 <<") "<<exEnergy<<" "<<Target4Momentu 427 #endif 428 429 // prepare the projectile residual - which 430 431 G4bool ProjectileIsAntiNucleus= 432 GetPrimaryProjectile()->GetDefinition 433 434 G4ThreeVector bst = GetPrimaryProjectile()- 435 436 G4int anAb=theProjectileNucleus->GetMassNum 437 G4int aZb=theProjectileNucleus->GetCharge() 438 G4int aLb=theProjectileNucleus->GetNumberOf 439 G4int numberOfExB = 0; 440 G4int numberOfChB = 0; 441 G4int numberOfHolesB = 0; 442 G4double exEnergyB = 0.0; 443 G4double Rb = theProjectileNucleus->GetNucl 444 G4LorentzVector Projectile4Momentum(0.,0.,0 445 446 // loop over the wounded projectile nucleus 447 theCurrentNucleon = 448 theProjectileNucleus->StartLoop() ? t 449 while(theCurrentNucleon) /* Loop checkin 450 { 451 if(theCurrentNucleon->AreYouHit()) { 452 ++numberOfHolesB; 453 ++numberOfExB; 454 --anAb; 455 if(!ProjectileIsAntiNucleus) { 456 aZb -= G4int(theCurrentNucleon->Ge 457 eplus + 0.1); 458 if (theCurrentNucleon->GetParticleType() 459 } else { 460 aZb += G4int(theCurrentNucleon->Ge 461 eplus - 0.1); 462 if (theCurrentNucleon->GetParticleType() 463 } 464 exEnergyB += theCurrentNucleon->GetBi 465 Projectile4Momentum -=theCurrentNucle 466 } 467 theCurrentNucleon = theProjectileNucleus 468 } 469 470 G4bool ExistTargetRemnant = G4double (nu 471 0.3* G4double (numberOfHoles + anA); 472 G4bool ExistProjectileRemnant= G4double (nu 473 0.3*G4double (numberOfHolesB + anAb); 474 475 #ifdef debugPrecoInt 476 G4cout<<"Projectile residual A Z (numberOfL 477 <<") "<<exEnergyB<<" "<<Projectile4Momentum 478 G4cout<<" ExistTargetRemnant ExistProjectil 479 <<ExistTargetRemnant<<" "<< ExistProj 480 #endif 481 //----------------------------------------- 482 // decay the strong resonances 483 G4ReactionProductVector * theTotalResult = 484 G4DecayKineticTracks decay(theSecondaries); 485 486 MakeCoalescence(theSecondaries); 487 488 #ifdef debugPrecoInt 489 G4cout<<"Secondary stable particles numb 490 #endif 491 492 #ifdef debugPrecoInt 493 G4LorentzVector secondary4Momemtum(0,0,0,0) 494 G4int SecondrNum(0); 495 #endif 496 497 // Loop over secondaries. 498 // We are assuming that only protons and ne 499 // and only antiprotons and antineutrons - 500 // not instead lambdas (or hyperons more ge 501 // (or anti-hyperons more generally) - for 502 // to be eventually reviewed later on, in p 503 // anti-hypernuclei are introduced, instead 504 // anti-hypernuclei which currently exist i 505 G4KineticTrackVector::iterator iter; 506 for(iter=theSecondaries->begin(); iter !=th 507 { 508 const G4ParticleDefinition* part = (*ite 509 G4LorentzVector aTrack4Momentum=(*iter)- 510 511 if( part != proton && part != neutron && 512 (part != ANTIproton && Projectile 513 (part != ANTIneutron && Projectile 514 { 515 G4ReactionProduct * theNew = new G4Re 516 theNew->SetMomentum(aTrack4Momentum.v 517 theNew->SetTotalEnergy(aTrack4Momentu 518 theNew->SetCreatorModelID((*iter)->GetCreat 519 theNew->SetParentResonanceDef((*iter)->GetP 520 theNew->SetParentResonanceID((*iter)->GetPa 521 theTotalResult->push_back(theNew); 522 #ifdef debugPrecoInt 523 SecondrNum++; 524 secondary4Momemtum += (*iter)->Get4Mo 525 G4cout<<"Secondary "<<SecondrNum<<" 526 <<theNew->GetDefinition()->GetP 527 <<theNew->GetMomentum()<<" "<<t 528 529 #endif 530 delete (*iter); 531 continue; 532 } 533 534 G4bool CanBeCapturedByTarget = false; 535 if( part == proton || part == neutron) 536 { 537 CanBeCapturedByTarget = ExistTargetRe 538 (-CaptureThreshold*G4Log( G4Uni 539 (aTrack4Momentum + Target4 540 aTrack4Momentum.mag() - Target4 541 ((*iter)->GetPosition().mag 542 } 543 // --------------------------- 544 G4LorentzVector Position((*iter)->GetPos 545 Position.boost(bst); 546 547 G4bool CanBeCapturedByProjectile = false 548 549 if( !ProjectileIsAntiNucleus && 550 ( part == proton || part == neutro 551 { 552 CanBeCapturedByProjectile = ExistProj 553 (-CaptureThreshold*G4Log( G4Uni 554 (aTrack4Momentum + Project 555 aTrack4Momentum.mag() - Project 556 (Position.vect().mag() < Rb 557 } 558 559 if( ProjectileIsAntiNucleus && 560 ( part == ANTIproton || part == AN 561 { 562 CanBeCapturedByProjectile = ExistProj 563 (-CaptureThreshold*G4Log( G4Uni 564 (aTrack4Momentum + Project 565 aTrack4Momentum.mag() - Project 566 (Position.vect().mag() < Rb 567 } 568 569 if(CanBeCapturedByTarget && CanBeCapture 570 { 571 if(G4UniformRand() < 0.5) 572 { CanBeCapturedByTarget = true; CanBe 573 else 574 { CanBeCapturedByTarget = false; CanB 575 } 576 577 if(CanBeCapturedByTarget) 578 { 579 // within the target nucleus, neutron 580 // now calculate A, Z of the fragmen 581 // number of exciton states 582 #ifdef debugPrecoInt 583 G4cout<<"Track is CapturedByTarget "< 584 <<aTrack4Momentum<<" "<<aTrack4 585 #endif 586 ++anA; 587 ++numberOfEx; 588 G4int Z = G4int(part->GetPDGCharge()/ 589 aZ += Z; 590 numberOfCh += Z; 591 Target4Momentum +=aTrack4Momentum; 592 delete (*iter); 593 } else if(CanBeCapturedByProjectile) 594 { 595 // within the projectile nucleus, neu 596 // now calculate A, Z of the fragmen 597 // number of exciton states 598 #ifdef debugPrecoInt 599 G4cout<<"Track is CapturedByProjectil 600 <<aTrack4Momentum<<" "<<aTrack4 601 #endif 602 ++anAb; 603 ++numberOfExB; 604 G4int Z = G4int(part->GetPDGCharge()/ 605 if( ProjectileIsAntiNucleus ) Z=-Z; 606 aZb += Z; 607 numberOfChB += Z; 608 Projectile4Momentum +=aTrack4Momentum 609 delete (*iter); 610 } else 611 { // the track is not captured 612 G4ReactionProduct * theNew = new G4Re 613 theNew->SetMomentum(aTrack4Momentum.v 614 theNew->SetTotalEnergy(aTrack4Momentu 615 theNew->SetCreatorModelID((*iter)->GetCreat 616 theNew->SetParentResonanceDef((*iter)->GetP 617 theNew->SetParentResonanceID((*iter)->GetPa 618 theTotalResult->push_back(theNew); 619 620 #ifdef debugPrecoInt 621 SecondrNum++; 622 secondary4Momemtum += (*iter)->Get4Mo 623 /* 624 G4cout<<"Secondary "<<SecondrNum<<" 625 <<theNew->GetDefinition()->GetP 626 <<secondary4Momemtum<<G4endl; 627 */ 628 #endif 629 delete (*iter); 630 continue; 631 } 632 } 633 delete theSecondaries; 634 //----------------------------------------- 635 636 #ifdef debugPrecoInt 637 G4cout<<"Final target residual A Z (numberO 638 <<") "<<exEnergy<<" "<<Target4Momentum<< 639 #endif 640 641 if(0!=anA ) 642 { 643 // We assume that the target residual is 644 G4double fMass = G4NucleiProperties::Ge 645 646 if((anA == theNucleus->GetMassNumber()) 647 {Target4Momentum.setE(fMass);} 648 649 G4double RemnMass=Target4Momentum.mag(); 650 651 if(RemnMass < fMass) 652 { 653 RemnMass=fMass + exEnergy; 654 Target4Momentum.setE(std::sqrt(Target 655 RemnMass*RemnMass)); 656 } else 657 { exEnergy=RemnMass-fMass;} 658 659 if( exEnergy < 0.) exEnergy=0.; 660 661 // Need to de-excite the remnant nucleus 662 G4Fragment anInitialState(anA, aZ, Targe 663 anInitialState.SetNumberOfParticles(numb 664 anInitialState.SetNumberOfCharged(number 665 anInitialState.SetNumberOfHoles(numberOf 666 anInitialState.SetCreatorModelID(secID); 667 668 G4ReactionProductVector * aPrecoResult = 669 theDeExcitation->DeExcite(anInitia 670 671 #ifdef debugPrecoInt 672 G4cout<<"Target fragment number "<<aP 673 #endif 674 675 // fill pre-compound part into the resul 676 for(unsigned int ll=0; ll<aPrecoResult-> 677 { 678 theTotalResult->push_back(aPrecoResul 679 #ifdef debugPrecoInt 680 G4cout<<"Target fragment "<<ll<<" 681 <<aPrecoResult->operator[](l 682 <<aPrecoResult->operator[](l 683 <<aPrecoResult->operator[](l 684 <<aPrecoResult->operator[](l 685 #endif 686 } 687 delete aPrecoResult; 688 } 689 690 //----------------------------------------- 691 if((anAb == theProjectileNucleus->GetMassNu 692 {Projectile4Momentum = GetPrimaryProjectile 693 694 #ifdef debugPrecoInt 695 G4cout<<"Final projectile residual A Z (num 696 <<aLb<<") "<<exEnergyB<<" "<<Projectile4Mom 697 <<Projectile4Momen 698 #endif 699 700 if(0!=anAb) 701 { 702 // The projectile residual can be a hype 703 G4double fMass = 0.0; 704 if ( aLb > 0 ) { 705 fMass = G4HyperNucleiProperties::GetNuclearM 706 } else { 707 fMass = G4NucleiProperties::GetNuclearMass(a 708 } 709 G4double RemnMass=Projectile4Momentum.ma 710 711 if(RemnMass < fMass) 712 { 713 RemnMass=fMass + exEnergyB; 714 Projectile4Momentum.setE(std::sqrt(Pr 715 RemnMass*RemnMass)); 716 } else 717 { exEnergyB=RemnMass-fMass;} 718 719 if( exEnergyB < 0.) exEnergyB=0.; 720 721 G4ThreeVector bstToCM =Projectile4Moment 722 Projectile4Momentum.boost(bstToCM); 723 724 // Need to de-excite the remnant nucleus 725 G4Fragment anInitialState(anAb, aZb, aLb 726 anInitialState.SetNumberOfParticles(numb 727 anInitialState.SetNumberOfCharged(number 728 anInitialState.SetNumberOfHoles(numberOf 729 anInitialState.SetCreatorModelID(secID); 730 731 G4ReactionProductVector * aPrecoResult = 732 theDeExcitation->DeExcite(anInitia 733 734 #ifdef debugPrecoInt 735 G4cout<<"Projectile fragment number " 736 #endif 737 738 // fill pre-compound part into the resul 739 for(unsigned int ll=0; ll<aPrecoResult-> 740 { 741 G4LorentzVector tmp=G4LorentzVector(a 742 a 743 tmp.boost(-bstToCM); // Transformatio 744 aPrecoResult->operator[](ll)->SetMome 745 aPrecoResult->operator[](ll)->SetTota 746 747 if(ProjectileIsAntiNucleus) 748 { 749 const G4ParticleDefinition * aFrag 750 const G4ParticleDefinition * LastF 751 if (aFragment == proton) {Las 752 else if(aFragment == neutron) {Las 753 else if(aFragment == lambda) {Las 754 else if(aFragment == deuteron){Las 755 else if(aFragment == triton) {Las 756 else if(aFragment == He3) {Las 757 else if(aFragment == He4) {Las 758 else {} 759 760 if (aLb > 0) { // Anti-hypernucle 761 if (aFragment == G4HyperTriton: 762 LastFragment=G4AntiHyperTriton::Definition 763 } else if (aFragment == G4HyperH4::Def 764 LastFragment=G4AntiHyperH4::Definition(); 765 } else if (aFragment == G4HyperAlpha:: 766 LastFragment=G4AntiHyperAlpha::Definition( 767 } else if (aFragment == G4HyperHe5::De 768 LastFragment=G4AntiHyperHe5::Definition(); 769 } else if (aFragment == G4DoubleHyperH 770 LastFragment=G4AntiDoubleHyperH4::Definiti 771 } else if (aFragment == G4DoubleHyperD 772 LastFragment=G4AntiDoubleHyperDoubleNeutro 773 } 774 } 775 776 aPrecoResult->operator[](ll)->SetD 777 } 778 779 #ifdef debugPrecoInt 780 G4cout<<"Projectile fragment "<<ll 781 <<aPrecoResult->operator[](l 782 <<aPrecoResult->operator[](l 783 <<aPrecoResult->operator[](l 784 <<aPrecoResult->operator[](l 785 #endif 786 787 theTotalResult->push_back(aPrecoResul 788 } 789 790 delete aPrecoResult; 791 } 792 793 return theTotalResult; 794 } 795 796 797 void G4GeneratorPrecompoundInterface::MakeCoal 798 // This method replaces pairs of proton-neut 799 // antiproton-antineutron - in the case of a 800 // momentum, with, respectively, deuterons a 801 // Note that in the case of hypernuclei or a 802 // are not considered for coalescence becaus 803 // are assumed not to exist. 804 805 if (!tracks) return; 806 807 G4double MassCut = deuteron->GetPDGMass() + 808 809 for ( std::size_t i = 0; i < tracks->size(); 810 811 G4KineticTrack* trackP = (*tracks)[i]; 812 if ( ! trackP ) continue; 813 if (trackP->GetDefinition() != proton) con 814 815 G4LorentzVector Prot4Mom = trackP->Get4Mom 816 G4LorentzVector ProtSPposition = G4Lorentz 817 818 for ( std::size_t j = 0; j < tracks->size( 819 820 G4KineticTrack* trackN = (*tracks)[j]; 821 if (! trackN ) continue; 822 if (trackN->GetDefinition() != neutron) 823 824 G4LorentzVector Neut4Mom = trackN->Get4M 825 G4LorentzVector NeutSPposition = G4Loren 826 G4double EffMass = (Prot4Mom + Neut4Mom) 827 828 if ( EffMass <= MassCut ) { // && (EffD 829 G4KineticTrack* aDeuteron = 830 new G4KineticTrack( deuteron, 831 (trackP->GetForm 832 (trackP->GetPosi 833 ( Prot4Mom 834 aDeuteron->SetCreatorModelID(secID); 835 tracks->push_back(aDeuteron); 836 delete trackP; delete trackN; 837 (*tracks)[i] = nullptr; (*tracks)[j] = 838 break; 839 } 840 } 841 } 842 843 // Find and remove null pointers created by 844 for ( G4int jj = (G4int)tracks->size()-1; jj 845 if ( ! (*tracks)[jj] ) tracks->erase(track 846 } 847 } 848