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 29 // ------------------------------------------- 30 // GEANT 4 class implementation file 31 // 32 // ---------------- G4FTFModel ---------- 33 // by Gunter Folger, May 1998. 34 // class implementing the excitation in 35 // 36 // Vladimir Uzhinsky, November 37 // simulation of nucleus-nucleus interac 38 // ------------------------------------------- 39 40 #include <utility> 41 42 #include "G4FTFModel.hh" 43 #include "G4ios.hh" 44 #include "G4PhysicalConstants.hh" 45 #include "G4SystemOfUnits.hh" 46 #include "G4FTFParameters.hh" 47 #include "G4FTFParticipants.hh" 48 #include "G4DiffractiveSplitableHadron.hh" 49 #include "G4InteractionContent.hh" 50 #include "G4LorentzRotation.hh" 51 #include "G4ParticleDefinition.hh" 52 #include "G4ParticleTable.hh" 53 #include "G4IonTable.hh" 54 #include "G4KineticTrack.hh" 55 #include "G4HyperNucleiProperties.hh" 56 #include "G4Exp.hh" 57 #include "G4Log.hh" 58 59 //============================================ 60 61 //#define debugFTFmodel 62 //#define debugReggeonCascade 63 //#define debugPutOnMassShell 64 //#define debugAdjust 65 //#define debugBuildString 66 67 68 //============================================ 69 70 G4FTFModel::G4FTFModel( const G4String& modelN 71 G4VPartonStringModel( modelName ), 72 theExcitation( new G4DiffractiveExcitation() 73 theElastic( new G4ElasticHNScattering() ), 74 theAnnihilation( new G4FTFAnnihilation() ) 75 { 76 // ---> JVY theParameters = 0; 77 theParameters = new G4FTFParameters(); 78 // 79 NumberOfInvolvedNucleonsOfTarget = 0; 80 NumberOfInvolvedNucleonsOfProjectile= 0; 81 for ( G4int i = 0; i < 250; ++i ) { 82 TheInvolvedNucleonsOfTarget[i] = 0; 83 TheInvolvedNucleonsOfProjectile[i] = 0; 84 } 85 86 //LowEnergyLimit = 2000.0*MeV; 87 LowEnergyLimit = 1000.0*MeV; 88 89 HighEnergyInter = true; 90 91 G4LorentzVector tmp( 0.0, 0.0, 0.0, 0.0 ); 92 ProjectileResidual4Momentum = tmp; 93 ProjectileResidualMassNumber = 0; 94 ProjectileResidualCharge = 0; 95 ProjectileResidualLambdaNumber = 0; 96 ProjectileResidualExcitationEnergy = 0.0; 97 98 TargetResidual4Momentum = tmp; 99 TargetResidualMassNumber = 0; 100 TargetResidualCharge = 0; 101 TargetResidualExcitationEnergy = 0.0; 102 103 Bimpact = -1.0; 104 BinInterval = false; 105 Bmin = 0.0; 106 Bmax = 0.0; 107 NumberOfProjectileSpectatorNucleons = 0; 108 NumberOfTargetSpectatorNucleons = 0; 109 NumberOfNNcollisions = 0; 110 111 SetEnergyMomentumCheckLevels( 2.0*perCent, 1 112 } 113 114 115 //============================================ 116 117 struct DeleteVSplitableHadron { void operator( 118 119 120 //============================================ 121 122 G4FTFModel::~G4FTFModel() { 123 // Because FTF model can be called for vari 124 // 125 // ---> NOTE (JVY): This statement below is 126 // theParameters must be erased at the end 127 // Thus the delete is also in G4FTFModel::G 128 // ---> JVY 129 // 130 if ( theParameters != 0 ) delete theParam 131 if ( theExcitation != 0 ) delete theExcit 132 if ( theElastic != 0 ) delete theElast 133 if ( theAnnihilation != 0 ) delete theAnnih 134 135 // Erasing of strings created at annihilati 136 if ( theAdditionalString.size() != 0 ) { 137 std::for_each( theAdditionalString.begin( 138 DeleteVSplitableHadron() ) 139 } 140 theAdditionalString.clear(); 141 142 // Erasing of target involved nucleons. 143 if ( NumberOfInvolvedNucleonsOfTarget != 0 144 for ( G4int i = 0; i < NumberOfInvolvedNu 145 G4VSplitableHadron* aNucleon = TheInvol 146 if ( aNucleon ) delete aNucleon; 147 } 148 } 149 150 // Erasing of projectile involved nucleons. 151 if ( NumberOfInvolvedNucleonsOfProjectile ! 152 for ( G4int i = 0; i < NumberOfInvolvedNu 153 G4VSplitableHadron* aNucleon = TheInvol 154 if ( aNucleon ) delete aNucleon; 155 } 156 } 157 } 158 159 160 //============================================ 161 162 void G4FTFModel::Init( const G4Nucleus& aNucle 163 164 theProjectile = aProjectile; 165 166 G4double PlabPerParticle( 0.0 ); // Laborat 167 168 #ifdef debugFTFmodel 169 G4cout << "FTF init Proj Name " << theProjec 170 << "FTF init Proj Mass " << theProjec 171 << " " << theProjectile.GetMomentum() 172 << "FTF init Proj B Q " << theProjec 173 << " " << (G4int) theProjectile.GetDe 174 << "FTF init Target A Z " << aNucleus 175 << " " << aNucleus.GetZ_asInt() << G4 176 #endif 177 178 theParticipants.Clean(); 179 180 theParticipants.SetProjectileNucleus( 0 ); 181 182 G4LorentzVector tmp( 0.0, 0.0, 0.0, 0.0 ); 183 ProjectileResidualMassNumber = 0; 184 ProjectileResidualCharge = 0; 185 ProjectileResidualLambdaNumber = 0; 186 ProjectileResidualExcitationEnergy = 0.0; 187 ProjectileResidual4Momentum = tmp; 188 189 TargetResidualMassNumber = aNucleus.Ge 190 TargetResidualCharge = aNucleus.Ge 191 TargetResidualExcitationEnergy = 0.0; 192 TargetResidual4Momentum = tmp; 193 G4double TargetResidualMass = G4ParticleTabl 194 ->GetIonMass( 195 196 TargetResidual4Momentum.setE( TargetResidual 197 198 if ( std::abs( theProjectile.GetDefinition() 199 // Projectile is a hadron : meson or baryo 200 ProjectileResidualMassNumber = std::abs( t 201 ProjectileResidualCharge = G4int( theProje 202 PlabPerParticle = theProjectile.GetMomentu 203 ProjectileResidualExcitationEnergy = 0.0; 204 //G4double ProjectileResidualMass = thePro 205 ProjectileResidual4Momentum.setVect( thePr 206 ProjectileResidual4Momentum.setE( theProje 207 if ( PlabPerParticle < LowEnergyLimit ) { 208 HighEnergyInter = false; 209 } else { 210 HighEnergyInter = true; 211 } 212 } else { 213 if ( theProjectile.GetDefinition()->GetBar 214 // Projectile is a nucleus 215 ProjectileResidualMassNumber = theProjec 216 ProjectileResidualCharge = G4int( thePro 217 ProjectileResidualLambdaNumber = theProj 218 PlabPerParticle = theProjectile.GetMomen 219 if ( PlabPerParticle < LowEnergyLimit ) 220 HighEnergyInter = false; 221 } else { 222 HighEnergyInter = true; 223 } 224 theParticipants.InitProjectileNucleus( P 225 ProjectileResidualLambdaNumber 226 } else if ( theProjectile.GetDefinition()- 227 // Projectile is an anti-nucleus 228 ProjectileResidualMassNumber = std::abs( 229 ProjectileResidualCharge = std::abs( G4i 230 ProjectileResidualLambdaNumber = theProj 231 PlabPerParticle = theProjectile.GetMomen 232 if ( PlabPerParticle < LowEnergyLimit ) 233 HighEnergyInter = false; 234 } else { 235 HighEnergyInter = true; 236 } 237 theParticipants.InitProjectileNucleus( P 238 P 239 theParticipants.GetProjectileNucleus()-> 240 G4Nucleon* aNucleon; 241 while ( ( aNucleon = theParticipants.Get 242 if ( aNucleon->GetDefinition() == G4Pr 243 aNucleon->SetParticleType( G4AntiPro 244 } else if ( aNucleon->GetDefinition() 245 aNucleon->SetParticleType( G4AntiNeu 246 } else if ( aNucleon->GetDefinition() 247 aNucleon->SetParticleType( G4AntiLambda::D 248 } 249 } 250 } 251 252 G4ThreeVector BoostVector = theProjectile. 253 theParticipants.GetProjectileNucleus()->Do 254 theParticipants.GetProjectileNucleus()->Do 255 ProjectileResidualExcitationEnergy = 0.0; 256 //G4double ProjectileResidualMass = thePro 257 ProjectileResidual4Momentum.setVect( thePr 258 ProjectileResidual4Momentum.setE( theProje 259 } 260 261 // Init target nucleus (assumed to be never 262 theParticipants.Init( aNucleus.GetA_asInt(), 263 264 NumberOfProjectileSpectatorNucleons = std::a 265 NumberOfTargetSpectatorNucleons = aNucleus.G 266 NumberOfNNcollisions = 0; 267 268 // reset/recalculate everything for the new 269 theParameters->InitForInteraction( theProjec 270 aNucleus. 271 272 if ( theAdditionalString.size() != 0 ) { 273 std::for_each( theAdditionalString.begin() 274 DeleteVSplitableHadron() ); 275 } 276 theAdditionalString.clear(); 277 278 #ifdef debugFTFmodel 279 G4cout << "FTF end of Init" << G4endl << G4e 280 #endif 281 282 // In the case of Hydrogen target, for non-i 283 // do NOT simulate quasi-elastic (by forcing 284 // elastic scatering in theParameters - whic 285 // This is necessary because in this case qu 286 // with only one nucleon would be identical 287 // and the latter is already included in the 288 // (i.e. G4HadronElasticProcess). 289 if ( std::abs( theProjectile.GetDefinition() 290 aNucleus.GetA_asInt() < 2 ) theParamete 291 292 if ( SampleBinInterval() ) theParticipants.S 293 } 294 295 296 //============================================ 297 298 G4ExcitedStringVector* G4FTFModel::GetStrings( 299 300 #ifdef debugFTFmodel 301 G4cout << "G4FTFModel::GetStrings() " << G4e 302 #endif 303 304 G4ExcitedStringVector* theStrings = new G4Ex 305 theParticipants.GetList( theProjectile, theP 306 307 SetImpactParameter( theParticipants.GetImpac 308 309 StoreInvolvedNucleon(); 310 311 G4bool Success( true ); 312 313 if ( HighEnergyInter ) { 314 ReggeonCascade(); 315 316 #ifdef debugFTFmodel 317 G4cout << "FTF PutOnMassShell " << G4endl; 318 #endif 319 320 Success = PutOnMassShell(); 321 322 #ifdef debugFTFmodel 323 G4cout << "FTF PutOnMassShell Success? " 324 #endif 325 326 } 327 328 #ifdef debugFTFmodel 329 G4cout << "FTF ExciteParticipants " << G4end 330 #endif 331 332 if ( Success ) Success = ExciteParticipants( 333 334 #ifdef debugFTFmodel 335 G4cout << "FTF ExciteParticipants Success? " 336 #endif 337 338 if ( Success ) { 339 340 #ifdef debugFTFmodel 341 G4cout << "FTF BuildStrings "; 342 #endif 343 344 BuildStrings( theStrings ); 345 346 #ifdef debugFTFmodel 347 G4cout << "FTF BuildStrings " << theString 348 << "FTF GetResiduals of Nuclei " << 349 #endif 350 351 GetResiduals(); 352 353 /* 354 if ( theParameters != 0 ) { 355 delete theParameters; 356 theParameters = 0; 357 } 358 */ 359 } else if ( ! GetProjectileNucleus() ) { 360 // Erase the hadron projectile 361 std::vector< G4VSplitableHadron* > primari 362 theParticipants.StartLoop(); 363 while ( theParticipants.Next() ) { /* Loo 364 const G4InteractionContent& interaction 365 // Do not allow for duplicates 366 if ( primaries.end() == 367 std::find( primaries.begin(), prima 368 primaries.push_back( interaction.GetPr 369 } 370 } 371 std::for_each( primaries.begin(), primarie 372 primaries.clear(); 373 } 374 375 // Cleaning of the memory 376 G4VSplitableHadron* aNucleon = 0; 377 378 // Erase the projectile nucleons 379 for ( G4int i = 0; i < NumberOfInvolvedNucle 380 aNucleon = TheInvolvedNucleonsOfProjectile 381 if ( aNucleon ) delete aNucleon; 382 } 383 NumberOfInvolvedNucleonsOfProjectile = 0; 384 385 // Erase the target nucleons 386 for ( G4int i = 0; i < NumberOfInvolvedNucle 387 aNucleon = TheInvolvedNucleonsOfTarget[i]- 388 if ( aNucleon ) delete aNucleon; 389 } 390 NumberOfInvolvedNucleonsOfTarget = 0; 391 392 #ifdef debugFTFmodel 393 G4cout << "End of FTF. Go to fragmentation" 394 << "To continue - enter 1, to stop - 395 #endif 396 397 theParticipants.Clean(); 398 399 return theStrings; 400 } 401 402 403 //============================================ 404 405 void G4FTFModel::StoreInvolvedNucleon() { 406 //To store nucleons involved in the interact 407 408 NumberOfInvolvedNucleonsOfTarget = 0; 409 410 G4V3DNucleus* theTargetNucleus = GetTargetNu 411 theTargetNucleus->StartLoop(); 412 413 G4Nucleon* aNucleon; 414 while ( ( aNucleon = theTargetNucleus->GetNe 415 if ( aNucleon->AreYouHit() ) { 416 TheInvolvedNucleonsOfTarget[NumberOfInvo 417 NumberOfInvolvedNucleonsOfTarget++; 418 } 419 } 420 421 #ifdef debugFTFmodel 422 G4cout << "G4FTFModel::StoreInvolvedNucleon 423 G4cout << "NumberOfInvolvedNucleonsOfTarget 424 << G4endl << G4endl; 425 #endif 426 427 428 if ( ! GetProjectileNucleus() ) return; // 429 430 // The projectile is a nucleus or an anti-nu 431 432 NumberOfInvolvedNucleonsOfProjectile = 0; 433 434 G4V3DNucleus* theProjectileNucleus = GetProj 435 theProjectileNucleus->StartLoop(); 436 437 G4Nucleon* aProjectileNucleon; 438 while ( ( aProjectileNucleon = theProjectile 439 if ( aProjectileNucleon->AreYouHit() ) { 440 // Projectile nucleon was involved in th 441 TheInvolvedNucleonsOfProjectile[NumberOf 442 NumberOfInvolvedNucleonsOfProjectile++; 443 } 444 } 445 446 #ifdef debugFTFmodel 447 G4cout << "NumberOfInvolvedNucleonsOfProject 448 << G4endl << G4endl; 449 #endif 450 return; 451 } 452 453 454 //============================================ 455 456 void G4FTFModel::ReggeonCascade() { 457 // Implementation of the reggeon theory insp 458 459 #ifdef debugReggeonCascade 460 G4cout << "G4FTFModel::ReggeonCascade ------ 461 << "theProjectile.GetTotalMomentum() 462 << "theProjectile.GetTotalEnergy() " 463 << "ExcitationE/WN " << theParameters 464 #endif 465 466 G4int InitNINt = NumberOfInvolvedNucleonsOfT 467 468 // Reggeon cascading in target nucleus 469 for ( G4int InvTN = 0; InvTN < InitNINt; Inv 470 G4Nucleon* aTargetNucleon = TheInvolvedNuc 471 472 G4double CreationTime = aTargetNucleon->Ge 473 474 G4double XofWoundedNucleon = aTargetNucleo 475 G4double YofWoundedNucleon = aTargetNucleo 476 477 G4V3DNucleus* theTargetNucleus = GetTarget 478 theTargetNucleus->StartLoop(); 479 480 G4Nucleon* Neighbour(0); 481 while ( ( Neighbour = theTargetNucleus->Ge 482 if ( ! Neighbour->AreYouHit() ) { 483 G4double impact2 = sqr( XofWoundedNucl 484 sqr( YofWoundedNucl 485 486 if ( G4UniformRand() < theParameters-> 487 G4Exp( -impact2 488 ) { 489 // The neighbour nucleon is involved 490 TheInvolvedNucleonsOfTarget[ NumberO 491 NumberOfInvolvedNucleonsOfTarget++; 492 493 G4VSplitableHadron* targetSplitable; 494 targetSplitable = new G4DiffractiveS 495 496 Neighbour->Hit( targetSplitable ); 497 targetSplitable->SetTimeOfCreation( 498 targetSplitable->SetStatus( 3 ); // 499 } 500 } 501 } 502 } 503 504 #ifdef debugReggeonCascade 505 G4cout << "Final NumberOfInvolvedNucleonsOfT 506 << NumberOfInvolvedNucleonsOfTarget < 507 #endif 508 509 if ( ! GetProjectileNucleus() ) return; 510 511 // Nucleus-Nucleus Interaction : Destruction 512 G4int InitNINp = NumberOfInvolvedNucleonsOfP 513 514 //for ( G4int InvPN = 0; InvPN < NumberOfInv 515 for ( G4int InvPN = 0; InvPN < InitNINp; Inv 516 G4Nucleon* aProjectileNucleon = TheInvolve 517 518 G4double CreationTime = aProjectileNucleon 519 520 G4double XofWoundedNucleon = aProjectileNu 521 G4double YofWoundedNucleon = aProjectileNu 522 523 G4V3DNucleus* theProjectileNucleus = GetPr 524 theProjectileNucleus->StartLoop(); 525 526 G4Nucleon* Neighbour( 0 ); 527 while ( ( Neighbour = theProjectileNucleus 528 if ( ! Neighbour->AreYouHit() ) { 529 G4double impact2= sqr( XofWoundedNucle 530 sqr( YofWoundedNucle 531 532 if ( G4UniformRand() < theParameters-> 533 G4Exp( -impact2 534 ) { 535 // The neighbour nucleon is involved 536 TheInvolvedNucleonsOfProjectile[ Num 537 NumberOfInvolvedNucleonsOfProjectile 538 539 G4VSplitableHadron* projectileSplita 540 projectileSplitable = new G4Diffract 541 542 Neighbour->Hit( projectileSplitable 543 projectileSplitable->SetTimeOfCreati 544 projectileSplitable->SetStatus( 3 ); 545 } 546 } 547 } 548 } 549 550 #ifdef debugReggeonCascade 551 G4cout << "NumberOfInvolvedNucleonsOfProject 552 << NumberOfInvolvedNucleonsOfProjecti 553 #endif 554 } 555 556 557 //============================================ 558 559 G4bool G4FTFModel::PutOnMassShell() { 560 561 G4bool isProjectileNucleus = false; 562 if ( GetProjectileNucleus() ) isProjectileNu 563 564 #ifdef debugPutOnMassShell 565 G4cout << "PutOnMassShell start " << G4endl; 566 if ( isProjectileNucleus ) { 567 G4cout << "PutOnMassShell for Nucleus_Nucl 568 } 569 #endif 570 571 G4LorentzVector Pprojectile( theProjectile.G 572 if ( Pprojectile.z() < 0.0 ) return false; 573 574 G4bool isOk = true; 575 576 G4LorentzVector Ptarget( 0.0, 0.0, 0.0, 0.0 577 G4LorentzVector PtargetResidual( 0.0, 0.0, 0 578 G4double SumMasses = 0.0; 579 G4V3DNucleus* theTargetNucleus = GetTargetNu 580 G4double TargetResidualMass = 0.0; 581 582 #ifdef debugPutOnMassShell 583 G4cout << "Target : "; 584 #endif 585 isOk = ComputeNucleusProperties( theTargetNu 586 TargetResid 587 TargetResid 588 if ( ! isOk ) return false; 589 590 G4double Mprojectile = 0.0; 591 G4double M2projectile = 0.0; 592 G4LorentzVector Pproj( 0.0, 0.0, 0.0, 0.0 ); 593 G4LorentzVector PprojResidual( 0.0, 0.0, 0.0 594 G4V3DNucleus* thePrNucleus = GetProjectileNu 595 G4double PrResidualMass = 0.0; 596 597 if ( ! isProjectileNucleus ) { // hadron-nu 598 Mprojectile = Pprojectile.mag(); 599 M2projectile = Pprojectile.mag2(); 600 SumMasses += Mprojectile + 20.0*MeV; 601 } else { // nucleus-nucleus or antinucleus- 602 #ifdef debugPutOnMassShell 603 G4cout << "Projectile : "; 604 #endif 605 isOk = ComputeNucleusProperties( thePrNucl 606 Projectil 607 Projectil 608 if ( ! isOk ) return false; 609 } 610 611 G4LorentzVector Psum = Pprojectile + Ptarget 612 G4double SqrtS = Psum.mag(); 613 G4double S = Psum.mag2(); 614 615 #ifdef debugPutOnMassShell 616 G4cout << "Psum " << Psum/GeV << " GeV" << G 617 << "SumMasses, PrResidualMass and Tar 618 << PrResidualMass/GeV << " " << Targe 619 #endif 620 621 if ( SqrtS < SumMasses ) return false; // I 622 623 // Try to consider also the excitation energ 624 // possible, with the available energy; othe 625 G4double savedSumMasses = SumMasses; 626 if ( isProjectileNucleus ) { 627 SumMasses -= std::sqrt( sqr( PrResidualMas 628 SumMasses += std::sqrt( sqr( PrResidualMas 629 + PprojResidual.pe 630 } 631 SumMasses -= std::sqrt( sqr( TargetResidualM 632 SumMasses += std::sqrt( sqr( TargetResidualM 633 + PtargetResidual.pe 634 635 if ( SqrtS < SumMasses ) { 636 SumMasses = savedSumMasses; 637 if ( isProjectileNucleus ) ProjectileResid 638 TargetResidualExcitationEnergy = 0.0; 639 } 640 641 TargetResidualMass += TargetResidualExcitati 642 if ( isProjectileNucleus ) PrResidualMass += 643 644 #ifdef debugPutOnMassShell 645 if ( isProjectileNucleus ) { 646 G4cout << "PrResidualMass ProjResidualExci 647 << ProjectileResidualExcitationEnergy << 648 } 649 G4cout << "TargetResidualMass TargetResidual 650 << TargetResidualExcitationEnergy << 651 << "Sum masses " << SumMasses/GeV << 652 #endif 653 654 // Sampling of nucleons what can transfer to 655 if ( isProjectileNucleus && thePrNucleus-> 656 isOk = GenerateDeltaIsobar( SqrtS, Numbe 657 TheInvolvedN 658 } 659 if ( theTargetNucleus->GetMassNumber() != 1 660 isOk = isOk && GenerateDeltaIsobar( Sqrt 661 TheI 662 } 663 if ( ! isOk ) return false; 664 665 // Now we know that it is kinematically poss 666 // of the involved nucleons (or correspondin 667 // We have to sample the kinematical variabl 668 // of the final state. The sampled kinematic 669 // Notice that the sampling of the transvers 670 // Fermi motion. 671 672 G4LorentzRotation toCms( -1*Psum.boostVector 673 G4LorentzVector Ptmp = toCms*Pprojectile; 674 if ( Ptmp.pz() <= 0.0 ) return false; // "S 675 676 G4LorentzRotation toLab( toCms.inverse() ); 677 678 G4double YprojectileNucleus = 0.0; 679 if ( isProjectileNucleus ) { 680 Ptmp = toCms*Pproj; 681 YprojectileNucleus = Ptmp.rapidity(); 682 } 683 Ptmp = toCms*Ptarget; 684 G4double YtargetNucleus = Ptmp.rapidity(); 685 686 // Ascribing of the involved nucleons Pt and 687 G4double DcorP = 0.0; 688 if ( isProjectileNucleus ) DcorP = theParame 689 G4double DcorT = theParameters->GetDof 690 G4double AveragePt2 = theParameters->GetPt2 691 G4double maxPtSquare = theParameters->GetMax 692 693 #ifdef debugPutOnMassShell 694 if ( isProjectileNucleus ) { 695 G4cout << "Y projectileNucleus " << Yproje 696 } 697 G4cout << "Y targetNucleus " << YtargetN 698 << "Dcor " << theParameters->GetDofNu 699 << " DcorP DcorT " << DcorP << " " << 700 #endif 701 702 G4double M2proj = M2projectile; // Initiali 703 G4double WplusProjectile = 0.0; 704 G4double M2target = 0.0; 705 G4double WminusTarget = 0.0; 706 G4int NumberOfTries = 0; 707 G4double ScaleFactor = 2.0; 708 G4bool OuterSuccess = true; 709 710 const G4int maxNumberOfLoops = 1000; 711 G4int loopCounter = 0; 712 do { // while ( ! OuterSuccess ) 713 OuterSuccess = true; 714 const G4int maxNumberOfInnerLoops = 10000; 715 do { // while ( SqrtS < Mprojectile + std 716 NumberOfTries++; 717 if ( NumberOfTries == 100*(NumberOfTries 718 // After many tries, it is convenient 719 // AveragePt2, so that the sampled mom 720 // involved nucleons (or corresponding delta 721 // it is more likely to satisfy the mo 722 ScaleFactor /= 2.0; 723 DcorP *= ScaleFactor; 724 DcorT *= ScaleFactor; 725 AveragePt2 *= ScaleFactor; 726 } 727 if ( isProjectileNucleus ) { 728 // Sampling of kinematical properties 729 isOk = SamplingNucleonKinematics( Aver 730 theP 731 PrRe 732 Numb 733 TheI 734 } 735 // Sampling of kinematical properties of 736 isOk = isOk && SamplingNucleonKinemati 737 738 739 740 741 #ifdef debugPutOnMassShell 742 G4cout << "SqrtS, Mp+Mt, Mp, Mt " << Sqr 743 << ( std::sqrt( M2proj ) + std::s 744 << std::sqrt( M2proj )/GeV << " " 745 #endif 746 if ( ! isOk ) return false; 747 } while ( ( SqrtS < std::sqrt( M2proj ) + 748 NumberOfTries < maxNumberOfInner 749 if ( NumberOfTries >= maxNumberOfInnerLoop 750 #ifdef debugPutOnMassShell 751 G4cout << "BAD situation: forced exit of 752 #endif 753 return false; 754 } 755 if ( isProjectileNucleus ) { 756 isOk = CheckKinematics( S, SqrtS, M2proj 757 NumberOfInvolved 758 TheInvolvedNucle 759 WminusTarget, Wp 760 } 761 isOk = isOk && CheckKinematics( S, SqrtS 762 NumberOf 763 WminusTa 764 if ( ! isOk ) return false; 765 } while ( ( ! OuterSuccess ) && 766 ++loopCounter < maxNumberOfLoops ) 767 if ( loopCounter >= maxNumberOfLoops ) { 768 #ifdef debugPutOnMassShell 769 G4cout << "BAD situation: forced exit of t 770 #endif 771 return false; 772 } 773 774 // Now the sampling is completed, and we can 775 // whole system. This is done first in the c 776 // to the lab frame. The transverse momentum 777 // the recoil of each hadron (nucleon or del 778 // to conserve (by construction) the transve 779 780 if ( ! isProjectileNucleus ) { // hadron-nu 781 782 G4double Pzprojectile = WplusProjectile/2. 783 G4double Eprojectile = WplusProjectile/2. 784 Pprojectile.setPz( Pzprojectile ); 785 Pprojectile.setE( Eprojectile ); 786 787 #ifdef debugPutOnMassShell 788 G4cout << "Proj after in CMS " << Pproject 789 #endif 790 791 Pprojectile.transform( toLab ); 792 theProjectile.SetMomentum( Pprojectile.vec 793 theProjectile.SetTotalEnergy( Pprojectile. 794 795 theParticipants.StartLoop(); 796 theParticipants.Next(); 797 G4VSplitableHadron* primary = theParticipa 798 primary->Set4Momentum( Pprojectile ); 799 800 #ifdef debugPutOnMassShell 801 G4cout << "Final proj. mom in Lab. " << pr 802 #endif 803 804 } else { // nucleus-nucleus or antinucleus- 805 806 isOk = FinalizeKinematics( WplusProjectile 807 ProjectileResid 808 TheInvolvedNucl 809 810 #ifdef debugPutOnMassShell 811 G4cout << "Projectile Residual4Momentum in 812 #endif 813 814 if ( ! isOk ) return false; 815 816 ProjectileResidual4Momentum.transform( toL 817 818 #ifdef debugPutOnMassShell 819 G4cout << "Projectile Residual4Momentum in 820 #endif 821 822 } 823 824 isOk = FinalizeKinematics( WminusTarget, fal 825 TargetResidualMas 826 TheInvolvedNucleo 827 828 #ifdef debugPutOnMassShell 829 G4cout << "Target Residual4Momentum in CMS " 830 #endif 831 832 if ( ! isOk ) return false; 833 834 TargetResidual4Momentum.transform( toLab ); 835 836 #ifdef debugPutOnMassShell 837 G4cout << "Target Residual4Momentum in Lab " 838 #endif 839 840 return true; 841 842 } 843 844 845 //============================================ 846 847 G4bool G4FTFModel::ExciteParticipants() { 848 849 #ifdef debugBuildString 850 G4cout << "G4FTFModel::ExciteParticipants() 851 #endif 852 853 G4bool Success( false ); 854 G4int MaxNumOfInelCollisions = G4int( thePar 855 if ( MaxNumOfInelCollisions > 0 ) { // Pla 856 G4double ProbMaxNumber = theParameters->Ge 857 if ( G4UniformRand() < ProbMaxNumber ) Max 858 } else { 859 // Plab < Pbound, normal application of FT 860 MaxNumOfInelCollisions = 1; 861 } 862 863 #ifdef debugBuildString 864 G4cout << "MaxNumOfInelCollisions per hadron 865 #endif 866 867 G4int CurrentInteraction( 0 ); 868 theParticipants.StartLoop(); 869 870 G4bool InnerSuccess( true ); 871 while ( theParticipants.Next() ) { /* Loop 872 CurrentInteraction++; 873 const G4InteractionContent& collision = th 874 G4VSplitableHadron* projectile = collision 875 G4Nucleon* ProjectileNucleon = collision.G 876 G4VSplitableHadron* target = collision.Get 877 G4Nucleon* TargetNucleon = collision.GetTa 878 879 #ifdef debugBuildString 880 G4cout << G4endl << "Interaction # Status 881 << collision.GetStatus() << G4endl 882 << target << G4endl << "projectile- 883 << projectile->GetStatus() << " " < 884 << "projectile->GetSoftC target->G 885 << " " << target->GetSoftCollisionC 886 #endif 887 888 if ( collision.GetStatus() ) { 889 if ( G4UniformRand() < theParameters->Ge 890 // Elastic scattering 891 892 #ifdef debugBuildString 893 G4cout << "Elastic scattering" << G4en 894 #endif 895 896 if ( ! HighEnergyInter ) { 897 G4bool Annihilation = false; 898 G4bool Result = AdjustNucleons( proj 899 Targ 900 if ( ! Result ) continue; 901 } 902 InnerSuccess = theElastic->ElasticSca 903 } else if ( G4UniformRand() > theParamet 904 // Inelastic scattering 905 906 #ifdef debugBuildString 907 G4cout << "Inelastic interaction" << G 908 << "MaxNumOfInelCollisions per 909 #endif 910 911 if ( ! HighEnergyInter ) { 912 G4bool Annihilation = false; 913 G4bool Result = AdjustNucleons( proj 914 Targ 915 if ( ! Result ) continue; 916 } 917 if ( G4UniformRand() < 918 ( 1.0 - target->GetSoftCollisionC 919 ( 1.0 - projectile->GetSoftCollis 920 //if ( ! HighEnergyInter ) { 921 // G4bool Annihilation = false; 922 // G4bool Result = AdjustNucleons( 923 // 924 // if ( ! Result ) continue; 925 //} 926 if ( theExcitation->ExciteParticipan 927 InnerSuccess = true; 928 NumberOfNNcollisions++; 929 #ifdef debugBuildString 930 G4cout << "FTF excitation Successf 931 // G4cout << "After pro " << proj 932 // << projectile->Get4Momen 933 // << "After tar " << targ 934 // << target->Get4Momentum( 935 #endif 936 } else { 937 InnerSuccess = theElastic->Elastic 938 #ifdef debugBuildString 939 G4cout << "FTF excitation Non Inne 940 << InnerSuccess << G4endl; 941 #endif 942 } 943 } else { // The inelastic interactiti 944 #ifdef debugBuildString 945 G4cout << "Elastic scat. at rejectio 946 #endif 947 //if ( ! HighEnergyInter ) { 948 // G4bool Annihilation = false; 949 // G4bool Result = AdjustNucleons( 950 // 951 // if ( ! Result) continue; 952 //} 953 InnerSuccess = theElastic->ElasticSc 954 } 955 } else { // Annihilation 956 957 #ifdef debugBuildString 958 G4cout << "Annihilation" << G4endl; 959 #endif 960 961 // At last, annihilation 962 if ( ! HighEnergyInter ) { 963 G4bool Annihilation = true; 964 G4bool Result = AdjustNucleons( proj 965 Targ 966 if ( ! Result ) continue; 967 } 968 969 G4VSplitableHadron* AdditionalString = 970 if ( theAnnihilation->Annihilate( proj 971 InnerSuccess = true; 972 #ifdef debugBuildString 973 G4cout << "Annihilation successfull. 974 << AdditionalString << G4endl 975 //G4cout << "After pro " << project 976 //G4cout << "After tar " << target- 977 #endif 978 979 if ( AdditionalString != 0 ) theAddi 980 981 NumberOfNNcollisions++; 982 983 // Skipping possible interactions of 984 while ( theParticipants.Next() ) { 985 G4InteractionContent& acollision = 986 G4VSplitableHadron* NextProjectile 987 G4VSplitableHadron* NextTargetNucl 988 if ( projectile == NextProjectileN 989 acollision.SetStatus( 0 ); 990 } 991 } 992 993 // Continue the interactions 994 theParticipants.StartLoop(); 995 for ( G4int i = 0; i < CurrentIntera 996 997 /* 998 if ( target->GetStatus() == 4 ) { 999 // Skipping possible interactions 1000 while ( theParticipants.Next() ) 1001 G4InteractionContent& acollisio 1002 G4VSplitableHadron* NextProject 1003 G4VSplitableHadron* NextTargetN 1004 if ( target == NextTargetNucleo 1005 } 1006 } 1007 theParticipants.StartLoop(); 1008 for ( G4int I = 0; I < CurrentInter 1009 */ 1010 1011 } 1012 } 1013 } 1014 1015 if( InnerSuccess ) Success = true; 1016 1017 #ifdef debugBuildString 1018 G4cout << "----------------------------- 1019 << "projectile->GetStatus target-> 1020 << " " << target->GetStatus() << G 1021 << projectile->GetSoftCollisionCou 1022 << G4endl << "ExciteParticipants() 1023 #endif 1024 1025 } // end of while ( theParticipants.Next() 1026 1027 return Success; 1028 } 1029 1030 1031 //=========================================== 1032 1033 G4bool G4FTFModel::AdjustNucleons( G4VSplitab 1034 G4Nucleon* 1035 G4VSplitab 1036 G4Nucleon* 1037 G4bool Ann 1038 1039 #ifdef debugAdjust 1040 G4cout << "AdjustNucleons ----------------- 1041 << "Proj is nucleus? " << GetProject 1042 << "Proj 4mom " << SelectedAntiBaryo 1043 << "Targ 4mom " << SelectedTargetNuc 1044 << "Pr ResidualMassNumber Pr Residua 1045 << ProjectileResidualMassNumber << " 1046 << ProjectileResidualExcitationEnerg 1047 << "Tr ResidualMassNumber Tr Residua 1048 << TargetResidualMassNumber << " " < 1049 << TargetResidualExcitationEnergy << 1050 << "Collis. pr tr " << SelectedAntiB 1051 << SelectedTargetNucleon->GetSoftCol 1052 #endif 1053 1054 if ( SelectedAntiBaryon->GetSoftCollisionCo 1055 SelectedTargetNucleon->GetSoftCollisio 1056 return true; // Selected hadrons were adj 1057 } 1058 1059 G4int interactionCase = 0; 1060 if ( ( ! GetProjectileNucleus() && 1061 SelectedAntiBaryon->GetSoftCollis 1062 SelectedTargetNucleon->GetSoftCol 1063 || 1064 ( SelectedAntiBaryon->GetSoftCollis 1065 SelectedTargetNucleon->GetSoftCol 1066 // The case of hadron-nucleus interaction 1067 // the case when projectile nuclear nucle 1068 // a collision, but target nucleon did no 1069 interactionCase = 1; 1070 #ifdef debugAdjust 1071 G4cout << "case 1, hA prcol=0 trcol=0, AA 1072 #endif 1073 if ( TargetResidualMassNumber < 1 ) { 1074 return false; 1075 } 1076 if ( SelectedAntiBaryon->Get4Momentum().r 1077 return false; 1078 } 1079 if ( TargetResidualMassNumber == 1 ) { 1080 TargetResidualMassNumber = 0; 1081 TargetResidualCharge = 0; 1082 TargetResidualExcitationEnergy = 0.0; 1083 SelectedTargetNucleon->Set4Momentum( Ta 1084 TargetResidual4Momentum = G4LorentzVect 1085 return true; 1086 } 1087 } else if ( SelectedAntiBaryon->GetSoftColl 1088 SelectedTargetNucleon->GetSoftC 1089 // It is assumed that in the case there i 1090 interactionCase = 2; 1091 #ifdef debugAdjust 1092 G4cout << "case 2, prcol=0 trcol#0" << G 1093 #endif 1094 if ( ProjectileResidualMassNumber < 1 ) { 1095 return false; 1096 } 1097 if ( ProjectileResidual4Momentum.rapidity 1098 SelectedTargetNucleon->Get4Momentum( 1099 return false; 1100 } 1101 if ( ProjectileResidualMassNumber == 1 ) 1102 ProjectileResidualMassNumber = 0; 1103 ProjectileResidualCharge = 0; 1104 ProjectileResidualExcitationEnergy = 0. 1105 SelectedAntiBaryon->Set4Momentum( Proje 1106 ProjectileResidual4Momentum = G4Lorentz 1107 return true; 1108 } 1109 } else { // It has to be a nucleus-nucleus 1110 interactionCase = 3; 1111 #ifdef debugAdjust 1112 G4cout << "case 3, prcol=0 trcol=0" << G 1113 #endif 1114 if ( ! GetProjectileNucleus() ) { 1115 return false; 1116 } 1117 #ifdef debugAdjust 1118 G4cout << "Proj res Init " << ProjectileR 1119 << "Targ res Init " << TargetResid 1120 << "ProjectileResidualMassNumber P 1121 << ProjectileResidualMassNumber << 1122 << " (" << ProjectileResidualLambd 1123 << "TargetResidualMassNumber Targe 1124 << " " << TargetResidualCharge << 1125 #endif 1126 } 1127 1128 CommonVariables common; 1129 G4int returnCode = AdjustNucleonsAlgorithm_ 1130 1131 1132 G4bool returnResult = false; 1133 if ( returnCode == 0 ) { 1134 returnResult = true; // Successfully end 1135 } else if ( returnCode == 1 ) { 1136 // The part before sampling has been succ 1137 returnResult = AdjustNucleonsAlgorithm_Sa 1138 if ( returnResult ) { // The sampling ha 1139 AdjustNucleonsAlgorithm_afterSampling( 1140 1141 } 1142 } 1143 1144 return returnResult; 1145 } 1146 1147 //------------------------------------------- 1148 1149 G4int G4FTFModel::AdjustNucleonsAlgorithm_bef 1150 1151 1152 1153 1154 1155 1156 // First of the three utility methods used 1157 // This method returns an integer code - in 1158 // "0" : successfully ended and nothing e 1159 // "1" : successfully completed, but the 1160 // "99" : unsuccessfully ended, nothing el 1161 G4int returnCode = 99; 1162 1163 G4double ExcitationEnergyPerWoundedNucleon 1164 1165 // some checks and initializations 1166 if ( interactionCase == 1 ) { 1167 common.Psum = SelectedAntiBaryon->Get4Mom 1168 #ifdef debugAdjust 1169 G4cout << "Targ res Init " << TargetResid 1170 #endif 1171 common.Pprojectile = SelectedAntiBaryon-> 1172 } else if ( interactionCase == 2 ) { 1173 common.Psum = ProjectileResidual4Momentum 1174 common.Pprojectile = ProjectileResidual4M 1175 } else if ( interactionCase == 3 ) { 1176 common.Psum = ProjectileResidual4Momentum 1177 common.Pprojectile = ProjectileResidual4M 1178 } 1179 1180 // transform momenta to cms and then rotate 1181 common.toCms = G4LorentzRotation( -1*common 1182 common.Ptmp = common.toCms * common.Pprojec 1183 common.toCms.rotateZ( -1*common.Ptmp.phi() 1184 common.toCms.rotateY( -1*common.Ptmp.theta( 1185 common.Pprojectile.transform( common.toCms 1186 common.toLab = common.toCms.inverse(); 1187 common.SqrtS = common.Psum.mag(); 1188 common.S = sqr( common.SqrtS ); 1189 1190 // get properties of the target residual an 1191 G4bool Stopping = false; 1192 if ( interactionCase == 1 ) { 1193 common.TResidualMassNumber = TargetResidu 1194 common.TResidualCharge = TargetResidual 1195 - G4int( TargetN 1196 common.TResidualExcitationEnergy = Targ 1197 - Exci 1198 if ( common.TResidualMassNumber <= 1 ) { 1199 common.TResidualExcitationEnergy = 0.0; 1200 } 1201 if ( common.TResidualMassNumber != 0 ) { 1202 common.TResidualMass = G4ParticleTable: 1203 ->GetIonMass( co 1204 } 1205 common.TNucleonMass = TargetNucleon->GetD 1206 common.SumMasses = SelectedAntiBaryon-> 1207 + common.TResidualMass 1208 #ifdef debugAdjust 1209 G4cout << "Annihilation " << Annihilation 1210 #endif 1211 } else if ( interactionCase == 2 ) { 1212 common.Ptarget = common.toCms * SelectedT 1213 common.TResidualMassNumber = ProjectileRe 1214 common.TResidualCharge = ProjectileResi 1215 - std::abs( G4in 1216 common.TResidualExcitationEnergy = Proj 1217 - Exci 1218 if ( common.TResidualMassNumber <= 1 ) { 1219 common.TResidualExcitationEnergy = 0.0; 1220 } 1221 if ( common.TResidualMassNumber != 0 ) { 1222 common.TResidualMass = G4ParticleTable: 1223 ->GetIonMass( co 1224 } 1225 common.TNucleonMass = ProjectileNucleon-> 1226 common.SumMasses = SelectedTargetNucleo 1227 + common.TResidualMass 1228 #ifdef debugAdjust 1229 G4cout << "SelectedTN.mag() PNMass + PRes 1230 << SelectedTargetNucleon->Get4Mome 1231 << common.TNucleonMass << " " << c 1232 #endif 1233 } else if ( interactionCase == 3 ) { 1234 common.Ptarget = common.toCms * TargetRes 1235 common.PResidualMassNumber = ProjectileRe 1236 common.PResidualCharge = ProjectileResi 1237 - std::abs( G4in 1238 common.PResidualLambdaNumber = Projectile 1239 if ( ProjectileNucleon->GetDefinition() = 1240 ProjectileNucleon->GetDefinition() == G4An 1241 --common.PResidualLambdaNumber; 1242 } 1243 common.PResidualExcitationEnergy = Proj 1244 - Exci 1245 if ( common.PResidualMassNumber <= 1 ) { 1246 common.PResidualExcitationEnergy = 0.0; 1247 } 1248 if ( common.PResidualMassNumber != 0 ) { 1249 if ( common.PResidualMassNumber == 1 ) 1250 if ( std::abs( common.PResidualCharge 1251 common.PResidualMass = G4Proton::De 1252 } else if ( common.PResidualLambdaNum 1253 common.PResidualMass = G4Lambda::De 1254 } else { 1255 common.PResidualMass = G4Neutron::D 1256 } 1257 } else { 1258 if ( common.PResidualLambdaNumber > 0 1259 if ( common.PResidualMassNumber == 1260 common.PResidualMass = G4Lambda:: 1261 if ( std::abs( common.PResidualCharge ) 1262 common.PResidualMass += G4Proton::Def 1263 } else if ( common.PResidualLambdaNumbe 1264 common.PResidualMass += G4Neutron::De 1265 } else { 1266 common.PResidualMass += G4Lambda::Def 1267 } 1268 } else { 1269 common.PResidualMass = G4HyperNucleiPro 1270 std::abs( common.PRes 1271 common.PResidualLambdaN 1272 } 1273 } else { 1274 common.PResidualMass = G4ParticleTable::G 1275 GetIonMass( std::a 1276 } 1277 } 1278 } 1279 common.PNucleonMass = ProjectileNucleon-> 1280 common.TResidualMassNumber = TargetResidu 1281 common.TResidualCharge = TargetResidual 1282 - G4int( TargetN 1283 common.TResidualExcitationEnergy = Targ 1284 - Exci 1285 if ( common.TResidualMassNumber <= 1 ) { 1286 common.TResidualExcitationEnergy = 0.0; 1287 } 1288 if ( common.TResidualMassNumber != 0 ) { 1289 common.TResidualMass = G4ParticleTable: 1290 GetIonMass( comm 1291 } 1292 common.TNucleonMass = TargetNucleon->GetD 1293 common.SumMasses = common.PNucleonMass 1294 + common.TResidualMass 1295 #ifdef debugAdjust 1296 G4cout << "PNucleonMass PResidualMass TNu 1297 << " " << common.PResidualMass << 1298 << common.TResidualMass << G4endl 1299 << "PResidualExcitationEnergy " << 1300 << "TResidualExcitationEnergy " << 1301 #endif 1302 } // End-if on interactionCase 1303 1304 if ( ! Annihilation ) { 1305 if ( common.SqrtS < common.SumMasses ) { 1306 #ifdef debugAdjust 1307 G4cout << "SqrtS < SumMasses " << commo 1308 #endif 1309 return returnCode; // Unsuccessfully e 1310 } 1311 if ( interactionCase == 1 || interactio 1312 if ( common.SqrtS < common.SumMasses + 1313 #ifdef debugAdjust 1314 G4cout << "TResidualExcitationEnergy 1315 #endif 1316 common.TResidualExcitationEnergy = co 1317 #ifdef debugAdjust 1318 G4cout << "TResidualExcitationEnergy 1319 #endif 1320 Stopping = true; 1321 return returnCode; // unsuccessfully 1322 } 1323 } else if ( interactionCase == 3 ) { 1324 #ifdef debugAdjust 1325 G4cout << "SqrtS < SumMasses + PResidua 1326 << common.SqrtS << " " << common 1327 << G4endl; 1328 #endif 1329 if ( common.SqrtS < common.SumMasses 1330 + common.TResidualE 1331 Stopping = true; 1332 if ( common.PResidualExcitationEnergy 1333 common.TResidualExcitationEnergy = 1334 } else if ( common.TResidualExcitatio 1335 common.PResidualExcitationEnergy = 1336 } else { 1337 G4double Fraction = ( common.SqrtS 1338 / ( common.PResidualExcitationEn 1339 common.PResidualExcitationEnergy *= 1340 common.TResidualExcitationEnergy *= 1341 } 1342 } 1343 } 1344 } // End-if on ! Annihilation 1345 1346 if ( Annihilation ) { 1347 #ifdef debugAdjust 1348 G4cout << "SqrtS < SumMasses - TNucleonMa 1349 << common.SumMasses - common.TNucl 1350 #endif 1351 if ( common.SqrtS < common.SumMasses - co 1352 return returnCode; // unsuccessfully e 1353 } 1354 #ifdef debugAdjust 1355 G4cout << "SqrtS < SumMasses " << common. 1356 #endif 1357 if ( common.SqrtS < common.SumMasses ) { 1358 if ( interactionCase == 2 || interact 1359 common.TResidualExcitationEnergy = 0. 1360 } 1361 common.TNucleonMass = common.SqrtS - 1362 - common.TResidua 1363 #ifdef debugAdjust 1364 G4cout << "TNucleonMass " << common.TNu 1365 #endif 1366 common.SumMasses = common.SqrtS - commo 1367 Stopping = true; 1368 #ifdef debugAdjust 1369 G4cout << "SqrtS < SumMasses " << commo 1370 #endif 1371 } 1372 if ( interactionCase == 1 || interactio 1373 if ( common.SqrtS < common.SumMasses + 1374 common.TResidualExcitationEnergy = co 1375 Stopping = true; 1376 } 1377 } else if ( interactionCase == 3 ) { 1378 if ( common.SqrtS < common.SumMasses 1379 + common.TResidualE 1380 Stopping = true; 1381 if ( common.PResidualExcitationEnergy 1382 common.TResidualExcitationEnergy = 1383 } else if ( common.TResidualExcitatio 1384 common.PResidualExcitationEnergy = 1385 } else { 1386 G4double Fraction = ( common.SqrtS 1387 ( common.PResidualExcitationEnerg 1388 common.PResidualExcitationEnergy *= 1389 common.TResidualExcitationEnergy *= 1390 } 1391 } 1392 } 1393 } // End-if on Annihilation 1394 1395 #ifdef debugAdjust 1396 G4cout << "Stopping " << Stopping << G4endl 1397 #endif 1398 1399 if ( Stopping ) { 1400 // All 3-momenta of particles = 0 1401 common.Ptmp.setPx( 0.0 ); common.Ptmp.set 1402 // New projectile 1403 if ( interactionCase == 1 ) { 1404 common.Ptmp.setE( SelectedAntiBaryon->G 1405 } else if ( interactionCase == 2 ) { 1406 common.Ptmp.setE( common.TNucleonMass ) 1407 } else if ( interactionCase == 3 ) { 1408 common.Ptmp.setE( common.PNucleonMass ) 1409 } 1410 #ifdef debugAdjust 1411 G4cout << "Proj stop " << common.Ptmp << 1412 #endif 1413 common.Pprojectile = common.Ptmp; 1414 common.Pprojectile.transform( common.toLa 1415 //---AR-Jul2019 : To avoid unphysical pro 1416 // original momentum of th 1417 G4LorentzVector saveSelectedAntiBaryon4Mo 1418 saveSelectedAntiBaryon4Momentum.transform 1419 //--- 1420 SelectedAntiBaryon->Set4Momentum( common. 1421 // New target nucleon 1422 if ( interactionCase == 1 || interactio 1423 common.Ptmp.setE( common.TNucleonMass ) 1424 } else if ( interactionCase == 2 ) { 1425 common.Ptmp.setE( SelectedTargetNucleon 1426 } 1427 #ifdef debugAdjust 1428 G4cout << "Targ stop " << common.Ptmp << 1429 #endif 1430 common.Ptarget = common.Ptmp; 1431 common.Ptarget.transform( common.toLab ); 1432 //---AR-Jul2019 : To avoid unphysical tar 1433 // momentum of the target 1434 G4LorentzVector saveSelectedTargetNucleon 1435 saveSelectedTargetNucleon4Momentum.transf 1436 //--- 1437 SelectedTargetNucleon->Set4Momentum( comm 1438 // New target residual 1439 if ( interactionCase == 1 || interactio 1440 common.Ptmp.setPx( 0.0 ); common.Ptmp.s 1441 TargetResidualMassNumber = common 1442 TargetResidualCharge = common 1443 TargetResidualExcitationEnergy = common 1444 //---AR-Jul2019 : To avoid unphysical t 1445 // original momentum of 1446 // This is a rough and s 1447 //common.Ptmp.setE( common.TResidualMas 1448 common.Ptmp.setPx( -saveSelectedTargetN 1449 common.Ptmp.setPy( -saveSelectedTargetN 1450 common.Ptmp.setPz( -saveSelectedTargetN 1451 common.Ptmp.setE( std::sqrt( sqr( commo 1452 //--- 1453 #ifdef debugAdjust 1454 G4cout << "Targ Resi stop " << common.P 1455 #endif 1456 common.Ptmp.transform( common.toLab ); 1457 TargetResidual4Momentum = common.Ptmp; 1458 } 1459 // New projectile residual 1460 if ( interactionCase == 2 || interactio 1461 common.Ptmp.setPx( 0.0 ); common.Ptmp.s 1462 if ( interactionCase == 2 ) { 1463 ProjectileResidualMassNumber = 1464 ProjectileResidualCharge = 1465 ProjectileResidualLambdaNumber = 0; // 1466 ProjectileResidualExcitationEnergy = 1467 common.Ptmp.setE( common.TResidualMas 1468 } else { 1469 ProjectileResidualMassNumber = 1470 ProjectileResidualCharge = 1471 ProjectileResidualLambdaNumber = common 1472 ProjectileResidualExcitationEnergy = 1473 //---AR-Jul2019 : To avoid unphysical 1474 // saved original mome 1475 // This is a rough and 1476 //common.Ptmp.setE( common.PResidualM 1477 common.Ptmp.setPx( -saveSelectedAntiB 1478 common.Ptmp.setPy( -saveSelectedAntiB 1479 common.Ptmp.setPz( -saveSelectedAntiB 1480 common.Ptmp.setE( std::sqrt( sqr( com 1481 //--- 1482 } 1483 #ifdef debugAdjust 1484 G4cout << "Proj Resi stop " << common.P 1485 #endif 1486 common.Ptmp.transform( common.toLab ); 1487 ProjectileResidual4Momentum = common.Pt 1488 } 1489 return returnCode = 0; // successfully e 1490 } // End-if on Stopping 1491 1492 // Initializations before sampling 1493 if ( interactionCase == 1 ) { 1494 common.Mprojectile = common.Pprojectile. 1495 common.M2projectile = common.Pprojectile. 1496 common.TResidual4Momentum = common.toCms 1497 common.YtargetNucleus = common.TResidual4 1498 common.TResidualMass += common.TResidualE 1499 } else if ( interactionCase == 2 ) { 1500 common.Mtarget = common.Ptarget.mag(); 1501 common.M2target = common.Ptarget.mag2(); 1502 common.TResidual4Momentum = common.toCms 1503 common.YprojectileNucleus = common.TResid 1504 common.TResidualMass += common.TResidualE 1505 } else if ( interactionCase == 3 ) { 1506 common.PResidual4Momentum = common.toCms 1507 common.YprojectileNucleus = common.PResid 1508 common.TResidual4Momentum = common.toCms* 1509 common.YtargetNucleus = common.TResidual4 1510 common.PResidualMass += common.PResidualE 1511 common.TResidualMass += common.TResidualE 1512 } 1513 #ifdef debugAdjust 1514 G4cout << "YprojectileNucleus " << common.Y 1515 #endif 1516 1517 return returnCode = 1; // successfully com 1518 } 1519 1520 1521 //------------------------------------------- 1522 1523 G4bool G4FTFModel::AdjustNucleonsAlgorithm_Sa 1524 1525 // Second of the three utility methods used 1526 // This method returns "false" if it fails 1527 1528 // Ascribing of the involved nucleons Pt an 1529 G4double Dcor = theParameters->GetDofNuclea 1530 G4double DcorP = 0.0, DcorT = 0.0; 1531 if ( ProjectileResidualMassNumber != 0 ) Dc 1532 if ( TargetResidualMassNumber != 0 ) Dc 1533 G4double AveragePt2 = theParameters->GetPt2 1534 G4double maxPtSquare = theParameters->GetMa 1535 1536 G4double ScaleFactor = 1.0; 1537 G4bool OuterSuccess = true; 1538 const G4int maxNumberOfLoops = 1000; 1539 const G4int maxNumberOfTries = 10000; 1540 G4int loopCounter = 0; 1541 G4int NumberOfTries = 0; 1542 do { // Outmost do while loop 1543 OuterSuccess = true; 1544 G4bool loopCondition = false; 1545 do { // Intermediate do while loop 1546 if ( NumberOfTries == 100*(NumberOfTrie 1547 // At large number of tries it would 1548 ScaleFactor /= 2.0; 1549 DcorP *= ScaleFactor; 1550 DcorT *= ScaleFactor; 1551 AveragePt2 *= ScaleFactor; 1552 #ifdef debugAdjust 1553 //G4cout << "NumberOfTries ScaleFacto 1554 #endif 1555 } 1556 1557 // Some kinematics 1558 if ( interactionCase == 1 ) { 1559 } else if ( interactionCase == 2 ) { 1560 #ifdef debugAdjust 1561 G4cout << "ProjectileResidualMassNumb 1562 #endif 1563 if ( ProjectileResidualMassNumber > 1 1564 common.PtNucleon = GaussianPt( Aver 1565 } else { 1566 common.PtNucleon = G4ThreeVector( 0 1567 } 1568 common.PtResidual = - common.PtNucleo 1569 common.Mprojectile = std::sqrt( sqr 1570 + std::sqrt( sqr 1571 #ifdef debugAdjust 1572 G4cout << "SqrtS < Mtarget + Mproject 1573 << " " << common.Mprojectile < 1574 #endif 1575 common.M2projectile = sqr( common.Mpr 1576 if ( common.SqrtS < common.Mtarget + 1577 OuterSuccess = false; 1578 loopCondition = true; 1579 continue; 1580 } 1581 } else if ( interactionCase == 3 ) { 1582 if ( ProjectileResidualMassNumber > 1 1583 common.PtNucleonP = GaussianPt( Ave 1584 } else { 1585 common.PtNucleonP = G4ThreeVector( 1586 } 1587 common.PtResidualP = - common.PtNucle 1588 if ( TargetResidualMassNumber > 1 ) { 1589 common.PtNucleonT = GaussianPt( Ave 1590 } else { 1591 common.PtNucleonT = G4ThreeVector( 1592 } 1593 common.PtResidualT = - common.PtNucle 1594 common.Mprojectile = std::sqrt( sqr 1595 + std::sqrt( sqr 1596 common.M2projectile = sqr( common.Mpr 1597 common.Mtarget = std::sqrt( sqr( co 1598 + std::sqrt( sqr( co 1599 common.M2target = sqr( common.Mtarget 1600 if ( common.SqrtS < common.Mprojectil 1601 OuterSuccess = false; 1602 loopCondition = true; 1603 continue; 1604 } 1605 } // End-if on interactionCase 1606 1607 G4int numberOfTimesExecuteInnerLoop = 1 1608 if ( interactionCase == 3 ) numberOfTim 1609 for ( G4int iExecute = 0; iExecute < nu 1610 1611 G4bool InnerSuccess = true; 1612 G4bool isTargetToBeHandled = ( intera 1613 ( inte 1614 G4bool condition = false; 1615 if ( isTargetToBeHandled ) { 1616 condition = ( TargetResidualMassNum 1617 } else { // Projectile to be handled 1618 condition = ( ProjectileResidualMas 1619 } 1620 if ( condition ) { 1621 const G4int maxNumberOfInnerLoops = 1622 G4int innerLoopCounter = 0; 1623 do { // Inner do while loop 1624 InnerSuccess = true; 1625 if ( isTargetToBeHandled ) { 1626 G4double Xcenter = 0.0; 1627 if ( interactionCase == 1 ) { 1628 common.PtNucleon = GaussianPt 1629 common.PtResidual = - common. 1630 common.Mtarget = std::sqrt( 1631 + std::sqrt( 1632 if ( common.SqrtS < common.Mp 1633 InnerSuccess = false; 1634 continue; 1635 } 1636 Xcenter = std::sqrt( sqr( com 1637 / common.Mtarget; 1638 } else { 1639 Xcenter = std::sqrt( sqr( com 1640 / common.Mtarget; 1641 } 1642 G4ThreeVector tmpX = GaussianPt 1643 common.XminusNucleon = Xcenter 1644 if ( common.XminusNucleon <= 0. 1645 InnerSuccess = false; 1646 continue; 1647 } 1648 common.XminusResidual = 1.0 - c 1649 } else { // Projectile to be han 1650 G4ThreeVector tmpX = GaussianPt 1651 G4double Xcenter = 0.0; 1652 if ( interactionCase == 2 ) { 1653 Xcenter = std::sqrt( sqr( com 1654 / common.Mprojectil 1655 } else { 1656 Xcenter = std::sqrt( sqr( com 1657 / common.Mprojectil 1658 } 1659 common.XplusNucleon = Xcenter + 1660 if ( common.XplusNucleon <= 0.0 1661 InnerSuccess = false; 1662 continue; 1663 } 1664 common.XplusResidual = 1.0 - co 1665 } // End-if on isTargetToBeHandl 1666 } while ( ( ! InnerSuccess ) && 1667 ++innerLoopCounter < ma 1668 if ( innerLoopCounter >= maxNumberO 1669 #ifdef debugAdjust 1670 G4cout << "BAD situation: forced 1671 #endif 1672 return false; 1673 } 1674 } else { // condition is false 1675 if ( isTargetToBeHandled ) { 1676 common.XminusNucleon = 1.0; 1677 common.XminusResidual = 1.0; // 1678 } else { // Projectile to be handl 1679 common.XplusNucleon = 1.0; 1680 common.XplusResidual = 1.0; // 1681 } 1682 } // End-if on condition 1683 1684 } // End of for loop on iExecute 1685 1686 if ( interactionCase == 1 ) { 1687 common.M2target = ( sqr( common.TN 1688 / common.XminusN 1689 + ( sqr( common.TR 1690 / common.XminusR 1691 loopCondition = ( common.SqrtS < comm 1692 } else if ( interactionCase == 2 ) { 1693 #ifdef debugAdjust 1694 G4cout << "TNucleonMass PtNucleon Xpl 1695 << common.PtNucleon << " " << 1696 << "TResidualMass PtResidual X 1697 << common.PtResidual << " " < 1698 #endif 1699 common.M2projectile = ( sqr( commo 1700 / common.Xpl 1701 + ( sqr( commo 1702 / common.Xpl 1703 #ifdef debugAdjust 1704 G4cout << "SqrtS < Mtarget + std::sqr 1705 << common.Mtarget << " " << st 1706 << common.Mtarget + std::sqrt( 1707 #endif 1708 loopCondition = ( common.SqrtS < comm 1709 } else if ( interactionCase == 3 ) { 1710 #ifdef debugAdjust 1711 G4cout << "PtNucleonP " << common.PtN 1712 << "XplusNucleon XplusResidual 1713 << " " << common.XplusResidual 1714 << "PtNucleonT " << common.PtN 1715 << "XminusNucleon XminusResidu 1716 << " " << common.XminusResidua 1717 #endif 1718 common.M2projectile = ( sqr( common 1719 / common.Xplu 1720 + ( sqr( common 1721 / common.Xplu 1722 common.M2target = ( sqr( common.TN 1723 / common.XminusN 1724 + ( sqr( common.TR 1725 / common.XminusR 1726 loopCondition = ( common.SqrtS < ( 1727 + std::sqrt( common.M2target ) ) 1728 } // End-if on interactionCase 1729 1730 } while ( loopCondition && 1731 ++NumberOfTries < maxNumberOfTr 1732 if ( NumberOfTries >= maxNumberOfTries ) 1733 #ifdef debugAdjust 1734 G4cout << "BAD situation: forced exit o 1735 #endif 1736 return false; 1737 } 1738 1739 // kinematics 1740 G4double Yprojectile = 0.0, YprojectileNu 1741 G4double DecayMomentum2 = sqr( common.S ) 1742 - 2.0 * ( commo 1743 + com 1744 if ( interactionCase == 1 ) { 1745 common.WminusTarget = ( common.S - co 1746 + std::sqrt( De 1747 common.WplusProjectile = common.SqrtS - 1748 common.Pzprojectile = common.WplusPro 1749 - common.M2projec 1750 common.Eprojectile = common.WplusPro 1751 + common.M2projec 1752 Yprojectile = 0.5 * G4Log( ( common. 1753 / ( common. 1754 #ifdef debugAdjust 1755 G4cout << "DecayMomentum2 " << DecayMom 1756 << "WminusTarget WplusProjectile 1757 << " " << common.WplusProjectile 1758 << "Yprojectile " << Yprojectile 1759 #endif 1760 common.Mt2targetNucleon = sqr( common.T 1761 common.PztargetNucleon = - common.Wminu 1762 + common.Mt2ta 1763 / ( 2.0 * co 1764 common.EtargetNucleon = common.Wminus 1765 + common.Mt2tar 1766 / ( 2.0 * com 1767 YtargetNucleon = 0.5 * G4Log( ( commo 1768 / ( commo 1769 #ifdef debugAdjust 1770 G4cout << "YtN Ytr YtN-Ytr " << " " << 1771 << " " << YtargetNucleon - commo 1772 << "YtN Ypr YtN-Ypr " << " " << 1773 << " " << YtargetNucleon - Yproj 1774 #endif 1775 if ( std::abs( YtargetNucleon - common. 1776 Yprojectile < YtargetNucleon ) { 1777 OuterSuccess = false; 1778 continue; 1779 } 1780 } else if ( interactionCase == 2 ) { 1781 common.WplusProjectile = ( common.S + 1782 + std::sqrt( 1783 common.WminusTarget = common.SqrtS - co 1784 common.Pztarget = - common.WminusTarget 1785 common.Etarget = common.WminusTarget 1786 Ytarget = 0.5 * G4Log( ( common.Etarg 1787 / ( common.Etarg 1788 #ifdef debugAdjust 1789 G4cout << "DecayMomentum2 " << DecayMom 1790 << "WminusTarget WplusProjectile 1791 << " " << common.WplusProjectile 1792 << "Ytarget " << Ytarget << G4en 1793 #endif 1794 common.Mt2projectileNucleon = sqr( comm 1795 common.PzprojectileNucleon = common.W 1796 - common.M 1797 / ( 2.0 1798 common.EprojectileNucleon = common.W 1799 + common.M 1800 / ( 2.0 1801 YprojectileNucleon = 0.5 * G4Log( ( c 1802 / ( c 1803 #ifdef debugAdjust 1804 G4cout << "YpN Ypr YpN-Ypr " << " " << 1805 << " " << YprojectileNucleon - c 1806 << "YpN Ytr YpN-Ytr " << " " << 1807 << " " << YprojectileNucleon - Y 1808 #endif 1809 if ( std::abs( YprojectileNucleon - com 1810 Ytarget > YprojectileNucleon ) { 1811 OuterSuccess = false; 1812 continue; 1813 } 1814 } else if ( interactionCase == 3 ) { 1815 common.WplusProjectile = ( common.S + 1816 + std::sqrt( 1817 common.WminusTarget = common.SqrtS - co 1818 common.Mt2projectileNucleon = sqr( comm 1819 common.PzprojectileNucleon = common.W 1820 - common.M 1821 / ( 2.0 1822 common.EprojectileNucleon = common.W 1823 + common.M 1824 / ( 2.0 1825 YprojectileNucleon = 0.5 * G4Log( ( c 1826 / ( c 1827 common.Mt2targetNucleon = sqr( common.T 1828 common.PztargetNucleon = - common.Wminu 1829 + common.Mt2ta 1830 / ( 2.0 * co 1831 common.EtargetNucleon = common.Wminu 1832 + common.Mt2ta 1833 / ( 2.0 * co 1834 YtargetNucleon = 0.5 * G4Log( ( commo 1835 / ( commo 1836 if ( std::abs( YtargetNucleon - common. 1837 std::abs( YprojectileNucleon - com 1838 YprojectileNucleon < YtargetNucleo 1839 OuterSuccess = false; 1840 continue; 1841 } 1842 } // End-if on interactionCase 1843 1844 } while ( ( ! OuterSuccess ) && 1845 ++loopCounter < maxNumberOfLoops 1846 if ( loopCounter >= maxNumberOfLoops ) { 1847 #ifdef debugAdjust 1848 G4cout << "BAD situation: forced exit of 1849 #endif 1850 return false; 1851 } 1852 1853 return true; 1854 } 1855 1856 //------------------------------------------- 1857 1858 void G4FTFModel::AdjustNucleonsAlgorithm_afte 1859 1860 1861 1862 // Third of the three utility methods used 1863 // and transform back. 1864 1865 // New projectile 1866 if ( interactionCase == 1 ) { 1867 common.Pprojectile.setPz( common.Pzprojec 1868 common.Pprojectile.setE( common.Eprojecti 1869 } else if ( interactionCase == 2 ) { 1870 common.Pprojectile.setPx( common.PtNucleo 1871 common.Pprojectile.setPy( common.PtNucleo 1872 common.Pprojectile.setPz( common.Pzprojec 1873 common.Pprojectile.setE( common.Eprojecti 1874 } else if ( interactionCase == 3 ) { 1875 common.Pprojectile.setPx( common.PtNucleo 1876 common.Pprojectile.setPy( common.PtNucleo 1877 common.Pprojectile.setPz( common.Pzprojec 1878 common.Pprojectile.setE( common.Eprojecti 1879 } 1880 #ifdef debugAdjust 1881 G4cout << "Proj after in CMS " << common.Pp 1882 #endif 1883 common.Pprojectile.transform( common.toLab 1884 SelectedAntiBaryon->Set4Momentum( common.Pp 1885 #ifdef debugAdjust 1886 G4cout << "Proj after in Lab " << common.Pp 1887 #endif 1888 1889 // New target nucleon 1890 if ( interactionCase == 1 ) { 1891 common.Ptarget.setPx( common.PtNucleon.x( 1892 common.Ptarget.setPy( common.PtNucleon.y( 1893 common.Ptarget.setPz( common.PztargetNucl 1894 common.Ptarget.setE( common.EtargetNucleo 1895 } else if ( interactionCase == 2 ) { 1896 common.Ptarget.setPz( common.Pztarget ); 1897 common.Ptarget.setE( common.Etarget ); 1898 } else if ( interactionCase == 3 ) { 1899 common.Ptarget.setPx( common.PtNucleonT.x 1900 common.Ptarget.setPy( common.PtNucleonT.y 1901 common.Ptarget.setPz( common.PztargetNucl 1902 common.Ptarget.setE( common.EtargetNucleo 1903 } 1904 #ifdef debugAdjust 1905 G4cout << "Targ after in CMS " << common.Pt 1906 #endif 1907 common.Ptarget.transform( common.toLab ); 1908 SelectedTargetNucleon->Set4Momentum( common 1909 #ifdef debugAdjust 1910 G4cout << "Targ after in Lab " << common.Pt 1911 #endif 1912 1913 // New target residual 1914 if ( interactionCase == 1 || interactionC 1915 TargetResidualMassNumber = common.T 1916 TargetResidualCharge = common.T 1917 TargetResidualExcitationEnergy = common.T 1918 #ifdef debugAdjust 1919 G4cout << "TargetResidualMassNumber Targe 1920 << TargetResidualMassNumber << " " 1921 << TargetResidualExcitationEnergy 1922 #endif 1923 if ( TargetResidualMassNumber != 0 ) { 1924 G4double Mt2 = 0.0; 1925 if ( interactionCase == 1 ) { 1926 Mt2 = sqr( common.TResidualMass ) + c 1927 TargetResidual4Momentum.setPx( common 1928 TargetResidual4Momentum.setPy( common 1929 } else { // interactionCase == 3 1930 Mt2 = sqr( common.TResidualMass ) + c 1931 TargetResidual4Momentum.setPx( common 1932 TargetResidual4Momentum.setPy( common 1933 } 1934 G4double Pz = - common.WminusTarget * c 1935 + Mt2 / ( 2.0 * common.Wm 1936 G4double E = common.WminusTarget * c 1937 + Mt2 / ( 2.0 * common.Wm 1938 TargetResidual4Momentum.setPz( Pz ); 1939 TargetResidual4Momentum.setE( E ) ; 1940 TargetResidual4Momentum.transform( comm 1941 } else { 1942 TargetResidual4Momentum = G4LorentzVect 1943 } 1944 #ifdef debugAdjust 1945 G4cout << "Tr N R " << common.Ptarget << 1946 #endif 1947 } 1948 1949 // New projectile residual 1950 if ( interactionCase == 2 || interactionC 1951 if ( interactionCase == 2 ) { 1952 ProjectileResidualMassNumber = co 1953 ProjectileResidualCharge = co 1954 ProjectileResidualExcitationEnergy = co 1955 ProjectileResidualLambdaNumber = co 1956 } else { // interactionCase == 3 1957 ProjectileResidualMassNumber = co 1958 ProjectileResidualCharge = co 1959 ProjectileResidualExcitationEnergy = co 1960 ProjectileResidualLambdaNumber = co 1961 } 1962 #ifdef debugAdjust 1963 G4cout << "ProjectileResidualMassNumber P 1964 << ProjectileResidualMassNumber << 1965 << ProjectileResidualLambdaNumber 1966 << ProjectileResidualExcitationEne 1967 #endif 1968 if ( ProjectileResidualMassNumber != 0 ) 1969 G4double Mt2 = 0.0; 1970 if ( interactionCase == 2 ) { 1971 Mt2 = sqr( common.TResidualMass ) + c 1972 ProjectileResidual4Momentum.setPx( co 1973 ProjectileResidual4Momentum.setPy( co 1974 } else { // interactionCase == 3 1975 Mt2 = sqr( common.PResidualMass ) + c 1976 ProjectileResidual4Momentum.setPx( co 1977 ProjectileResidual4Momentum.setPy( co 1978 } 1979 G4double Pz = common.WplusProjectile 1980 - Mt2 / ( 2.0 * common.Wp 1981 G4double E = common.WplusProjectile 1982 + Mt2 / ( 2.0 * common.Wp 1983 ProjectileResidual4Momentum.setPz( Pz ) 1984 ProjectileResidual4Momentum.setE( E ); 1985 ProjectileResidual4Momentum.transform( 1986 } else { 1987 ProjectileResidual4Momentum = G4Lorentz 1988 } 1989 #ifdef debugAdjust 1990 G4cout << "Pr N R " << common.Pprojectile 1991 << " " << ProjectileResidual 1992 #endif 1993 } 1994 1995 } 1996 1997 1998 //=========================================== 1999 2000 void G4FTFModel::BuildStrings( G4ExcitedStrin 2001 // Loop over all collisions; find all prima 2002 // (targets may be duplicate in the List (t 2003 2004 G4ExcitedString* FirstString( 0 ); // I 2005 G4ExcitedString* SecondString( 0 ); // t 2006 2007 if ( ! GetProjectileNucleus() ) { 2008 2009 std::vector< G4VSplitableHadron* > primar 2010 theParticipants.StartLoop(); 2011 while ( theParticipants.Next() ) { /* Lo 2012 const G4InteractionContent& interaction 2013 // do not allow for duplicates ... 2014 if ( interaction.GetStatus() ) { 2015 if ( primaries.end() == std::find( pr 2016 in 2017 primaries.push_back( interaction.Ge 2018 } 2019 } 2020 } 2021 2022 #ifdef debugBuildString 2023 G4cout << "G4FTFModel::BuildStrings()" << 2024 << "Number of projectile strings " 2025 #endif 2026 2027 for ( unsigned int ahadron = 0; ahadron < 2028 G4bool isProjectile( true ); 2029 //G4cout << "primaries[ ahadron ] " << 2030 //if ( primaries[ ahadron ]->GetStatus( 2031 FirstString = 0; SecondString = 0; 2032 if ( primaries[ahadron]->GetStatus() == 2033 theExcitation->CreateStrings( primarie 2034 FirstStr 2035 NumberOfProjectileSpectatorNucleons--; 2036 } else if ( primaries[ahadron]->GetStat 2037 && primaries[ahadron]->GetSoft 2038 theExcitation->CreateStrings( primarie 2039 FirstStr 2040 NumberOfProjectileSpectatorNucleons--; 2041 } else if ( primaries[ahadron]->GetStat 2042 && primaries[ahadron]->GetSoft 2043 G4LorentzVector ParticleMomentum=prima 2044 G4KineticTrack* aTrack = new G4Kinetic 2045 2046 2047 2048 FirstString = new G4ExcitedString( aTr 2049 } else if (primaries[ahadron]->GetStatu 2050 G4LorentzVector ParticleMomentum=prima 2051 G4KineticTrack* aTrack = new G4Kinetic 2052 2053 2054 2055 FirstString = new G4ExcitedString( aTr 2056 NumberOfProjectileSpectatorNucleons--; 2057 } else { 2058 G4cout << "Something wrong in FTF Mod 2059 } 2060 2061 if ( FirstString != 0 ) strings->push_ 2062 if ( SecondString != 0 ) strings->push_ 2063 2064 #ifdef debugBuildString 2065 G4cout << "FirstString & SecondString? 2066 if ( FirstString->IsExcited() ) { 2067 G4cout << "Quarks on the FirstString 2068 << " " << FirstString->GetLeft 2069 } else { 2070 G4cout << "Kinetic track is stored" < 2071 } 2072 #endif 2073 2074 } 2075 2076 #ifdef debugBuildString 2077 if ( FirstString->IsExcited() ) { 2078 G4cout << "Check 1 string " << strings- 2079 << " " << strings->operator[](0) 2080 } 2081 #endif 2082 2083 std::for_each( primaries.begin(), primari 2084 primaries.clear(); 2085 2086 } else { // Projectile is a nucleus 2087 2088 #ifdef debugBuildString 2089 G4cout << "Building of projectile-like st 2090 #endif 2091 2092 G4bool isProjectile = true; 2093 for ( G4int ahadron = 0; ahadron < Number 2094 2095 #ifdef debugBuildString 2096 G4cout << "Nucleon #, status, intCount 2097 << TheInvolvedNucleonsOfProjecti 2098 << " " << TheInvolvedNucleonsOfP 2099 ->GetSoftCollisionCoun 2100 #endif 2101 2102 G4VSplitableHadron* aProjectile = 2103 TheInvolvedNucleonsOfProjectile[ ah 2104 2105 #ifdef debugBuildString 2106 G4cout << G4endl << "ahadron aProjectil 2107 << " " << aProjectile->GetStatus 2108 #endif 2109 2110 FirstString = 0; SecondString = 0; 2111 if ( aProjectile->GetStatus() == 0 ) { 2112 2113 #ifdef debugBuildString 2114 G4cout << "Case1 aProjectile->GetStat 2115 #endif 2116 2117 theExcitation->CreateStrings( 2118 TheInvolvedNucleon 2119 isProjectile, Firs 2120 NumberOfProjectileSpectatorNucleons-- 2121 } else if ( aProjectile->GetStatus() == 2122 // Nucleon took part in diffractive i 2123 2124 #ifdef debugBuildString 2125 G4cout << "Case2 aProjectile->GetStat 2126 #endif 2127 2128 theExcitation->CreateStrings( 2129 TheInvolvedNucleon 2130 isProjectile, Firs 2131 NumberOfProjectileSpectatorNucleons-- 2132 } else if ( aProjectile->GetStatus() == 2133 HighEnergyInter ) { 2134 // Nucleon was considered as a parici 2135 // but the interaction was skipped du 2136 // It is now considered as an involve 2137 2138 #ifdef debugBuildString 2139 G4cout << "Case3 aProjectile->GetStat 2140 #endif 2141 2142 G4LorentzVector ParticleMomentum = aP 2143 G4KineticTrack* aTrack = new G4Kineti 2144 2145 2146 2147 FirstString = new G4ExcitedString( aT 2148 2149 #ifdef debugBuildString 2150 G4cout << " Strings are built for nuc 2151 << " the interaction was skipp 2152 #endif 2153 2154 } else if ( aProjectile->GetStatus() == 2155 // Nucleon which was involved in the 2156 2157 #ifdef debugBuildString 2158 G4cout << "Case4 aProjectile->GetStat 2159 #endif 2160 2161 G4LorentzVector ParticleMomentum = aP 2162 G4KineticTrack* aTrack = new G4Kineti 2163 2164 2165 2166 FirstString = new G4ExcitedString( aT 2167 2168 #ifdef debugBuildString 2169 G4cout << " Strings are build for inv 2170 #endif 2171 2172 if ( aProjectile->GetStatus() == 2 ) 2173 } else { 2174 2175 #ifdef debugBuildString 2176 G4cout << "Case5 " << G4endl; 2177 #endif 2178 2179 //TheInvolvedNucleonsOfProjectile[ ah 2180 //G4cout << TheInvolvedNucleonsOfProj 2181 2182 #ifdef debugBuildString 2183 G4cout << " No string" << G4endl; 2184 #endif 2185 2186 } 2187 2188 if ( FirstString != 0 ) strings->push_ 2189 if ( SecondString != 0 ) strings->push_ 2190 } 2191 } 2192 2193 #ifdef debugBuildString 2194 G4cout << "Building of target-like strings" 2195 #endif 2196 2197 G4bool isProjectile = false; 2198 for ( G4int ahadron = 0; ahadron < NumberOf 2199 G4VSplitableHadron* aNucleon = TheInvolve 2200 2201 #ifdef debugBuildString 2202 G4cout << "Nucleon #, status, intCount " 2203 << aNucleon->GetStatus() << " " << 2204 #endif 2205 2206 FirstString = 0 ; SecondString = 0; 2207 2208 if ( aNucleon->GetStatus() == 0 ) { // A 2209 theExcitation->CreateStrings( aNucleon, 2210 FirstStri 2211 NumberOfTargetSpectatorNucleons--; 2212 2213 #ifdef debugBuildString 2214 G4cout << " 1 case A string is build" < 2215 #endif 2216 2217 } else if ( aNucleon->GetStatus() == 1 & 2218 // A nucleon took part in diffractive i 2219 theExcitation->CreateStrings( aNucleon, 2220 FirstStri 2221 2222 #ifdef debugBuildString 2223 G4cout << " 2 case A string is build, n 2224 #endif 2225 2226 NumberOfTargetSpectatorNucleons--; 2227 2228 } else if ( aNucleon->GetStatus() == 1 & 2229 HighEnergyInter ) { 2230 // A nucleon was considered as a partic 2231 // its interactions were skipped. It wi 2232 // at high energies. 2233 2234 G4LorentzVector ParticleMomentum = aNuc 2235 G4KineticTrack* aTrack = new G4KineticT 2236 2237 2238 2239 2240 FirstString = new G4ExcitedString( aTra 2241 2242 #ifdef debugBuildString 2243 G4cout << "3 case A string is build" << 2244 #endif 2245 2246 } else if ( aNucleon->GetStatus() == 1 & 2247 ! HighEnergyInter ) { 2248 // A nucleon was considered as a partic 2249 // its interactions were skipped. It wi 2250 // at low energies energies. 2251 aNucleon->SetStatus( 5 ); // 4->5 2252 // ??? delete aNucleon; 2253 2254 #ifdef debugBuildString 2255 G4cout << "4 case A string is not build 2256 #endif 2257 2258 } else if ( aNucleon->GetStatus() == 2 | 2259 aNucleon->GetStatus() == 3 ) 2260 G4LorentzVector ParticleMomentum = aNuc 2261 G4KineticTrack* aTrack = new G4KineticT 2262 2263 2264 2265 FirstString = new G4ExcitedString( aTra 2266 2267 #ifdef debugBuildString 2268 G4cout << "5 case A string is build" << 2269 #endif 2270 2271 if ( aNucleon->GetStatus() == 2 ) Numbe 2272 2273 } else { 2274 2275 #ifdef debugBuildString 2276 G4cout << "6 case No string" << G4endl; 2277 #endif 2278 2279 } 2280 2281 if ( FirstString != 0 ) strings->push_ba 2282 if ( SecondString != 0 ) strings->push_ba 2283 2284 } 2285 2286 #ifdef debugBuildString 2287 G4cout << G4endl << "theAdditionalString.si 2288 << G4endl << G4endl; 2289 #endif 2290 2291 isProjectile = true; 2292 if ( theAdditionalString.size() != 0 ) { 2293 for ( unsigned int ahadron = 0; ahadron 2294 //if ( theAdditionalString[ ahadron ]-> 2295 FirstString = 0; SecondString = 0; 2296 theExcitation->CreateStrings( theAdditi 2297 FirstStri 2298 if ( FirstString != 0 ) strings->push_ 2299 if ( SecondString != 0 ) strings->push_ 2300 } 2301 } 2302 2303 //for ( unsigned int ahadron = 0; ahadron < 2304 // G4cout << ahadron << " " << strings->op 2305 // << " " << strings->operator[]( a 2306 //} 2307 //G4cout << "------------------------" << G 2308 2309 return; 2310 } 2311 2312 2313 //=========================================== 2314 2315 void G4FTFModel::GetResiduals() { 2316 // This method is needed for the correct ap 2317 2318 #ifdef debugFTFmodel 2319 G4cout << "GetResiduals(): HighEnergyInter? 2320 << HighEnergyInter << " " << GetProj 2321 #endif 2322 2323 if ( HighEnergyInter ) { 2324 2325 #ifdef debugFTFmodel 2326 G4cout << "NumberOfInvolvedNucleonsOfTarg 2327 #endif 2328 2329 G4double DeltaExcitationE = TargetResidua 2330 G4double( Num 2331 G4LorentzVector DeltaPResidualNucleus = T 2332 G 2333 2334 for ( G4int i = 0; i < NumberOfInvolvedNu 2335 G4Nucleon* aNucleon = TheInvolvedNucleo 2336 2337 #ifdef debugFTFmodel 2338 G4VSplitableHadron* targetSplitable = a 2339 G4cout << i << " Hit? " << aNucleon->Ar 2340 if ( targetSplitable ) G4cout << i << " 2341 #endif 2342 2343 G4LorentzVector tmp = -DeltaPResidualNu 2344 aNucleon->SetMomentum( tmp ); 2345 aNucleon->SetBindingEnergy( DeltaExcita 2346 } 2347 2348 if ( TargetResidualMassNumber != 0 ) { 2349 G4ThreeVector bstToCM = TargetResidual4 2350 2351 G4V3DNucleus* theTargetNucleus = GetTar 2352 G4LorentzVector residualMomentum( 0.0, 2353 G4Nucleon* aNucleon = 0; 2354 theTargetNucleus->StartLoop(); 2355 while ( ( aNucleon = theTargetNucleus-> 2356 if ( ! aNucleon->AreYouHit() ) { 2357 G4LorentzVector tmp = aNucleon->Get 2358 aNucleon->SetMomentum( tmp ); 2359 residualMomentum += tmp; 2360 } 2361 } 2362 2363 residualMomentum /= TargetResidualMassN 2364 2365 G4double Mass = TargetResidual4Momentum 2366 G4double SumMasses = 0.0; 2367 2368 aNucleon = 0; 2369 theTargetNucleus->StartLoop(); 2370 while ( ( aNucleon = theTargetNucleus-> 2371 if ( ! aNucleon->AreYouHit() ) { 2372 G4LorentzVector tmp = aNucleon->Get 2373 G4double E = std::sqrt( tmp.vect(). 2374 sqr( aNucle 2375 tmp.setE( E ); aNucleon->SetMoment 2376 SumMasses += E; 2377 } 2378 } 2379 2380 G4double Chigh = Mass / SumMasses; G4do 2381 const G4int maxNumberOfLoops = 1000; 2382 G4int loopCounter = 0; 2383 do { 2384 C = ( Chigh + Clow ) / 2.0; 2385 SumMasses = 0.0; 2386 aNucleon = 0; 2387 theTargetNucleus->StartLoop(); 2388 while ( ( aNucleon = theTargetNucleus 2389 if ( ! aNucleon->AreYouHit() ) { 2390 G4LorentzVector tmp = aNucleon->G 2391 G4double E = std::sqrt( tmp.vect( 2392 sqr( aNuc 2393 SumMasses += E; 2394 } 2395 } 2396 2397 if ( SumMasses > Mass ) Chigh = C; 2398 else Clow = C; 2399 2400 } while ( Chigh - Clow > 0.01 && 2401 ++loopCounter < maxNumberOfLo 2402 if ( loopCounter >= maxNumberOfLoops ) 2403 #ifdef debugFTFmodel 2404 G4cout << "BAD situation: forced exit 2405 << "\t return immediately from 2406 #endif 2407 return; 2408 } 2409 2410 aNucleon = 0; 2411 theTargetNucleus->StartLoop(); 2412 while ( ( aNucleon = theTargetNucleus-> 2413 if ( !aNucleon->AreYouHit() ) { 2414 G4LorentzVector tmp = aNucleon->Get 2415 G4double E = std::sqrt( tmp.vect(). 2416 sqr( aNucle 2417 tmp.setE( E ); tmp.boost( -bstToCM 2418 aNucleon->SetMomentum( tmp ); 2419 } 2420 } 2421 } 2422 2423 if ( ! GetProjectileNucleus() ) return; 2424 2425 #ifdef debugFTFmodel 2426 G4cout << "NumberOfInvolvedNucleonsOfProj 2427 << G4endl << "ProjectileResidualEx 2428 << ProjectileResidualExcitationEne 2429 #endif 2430 2431 DeltaExcitationE = ProjectileResidualExci 2432 G4double( NumberOfInvo 2433 DeltaPResidualNucleus = ProjectileResidua 2434 G4double( NumberO 2435 2436 for ( G4int i = 0; i < NumberOfInvolvedNu 2437 G4Nucleon* aNucleon = TheInvolvedNucleo 2438 2439 #ifdef debugFTFmodel 2440 G4VSplitableHadron* projSplitable = aNu 2441 G4cout << i << " Hit? " << aNucleon->Ar 2442 if ( projSplitable ) G4cout << i << "St 2443 #endif 2444 2445 G4LorentzVector tmp = -DeltaPResidualNu 2446 aNucleon->SetMomentum( tmp ); 2447 aNucleon->SetBindingEnergy( DeltaExcita 2448 } 2449 2450 if ( ProjectileResidualMassNumber != 0 ) 2451 G4ThreeVector bstToCM = ProjectileResid 2452 2453 G4V3DNucleus* theProjectileNucleus = Ge 2454 G4LorentzVector residualMomentum( 0.0, 2455 G4Nucleon* aNucleon = 0; 2456 theProjectileNucleus->StartLoop(); 2457 while ( ( aNucleon = theProjectileNucle 2458 if ( ! aNucleon->AreYouHit() ) { 2459 G4LorentzVector tmp = aNucleon->Get 2460 aNucleon->SetMomentum( tmp ); 2461 residualMomentum += tmp; 2462 } 2463 } 2464 2465 residualMomentum /= ProjectileResidualM 2466 2467 G4double Mass = ProjectileResidual4Mome 2468 G4double SumMasses= 0.0; 2469 2470 aNucleon = 0; 2471 theProjectileNucleus->StartLoop(); 2472 while ( ( aNucleon = theProjectileNucle 2473 if ( ! aNucleon->AreYouHit() ) { 2474 G4LorentzVector tmp = aNucleon->Get 2475 G4double E=std::sqrt( tmp.vect().ma 2476 sqr(aNucleon- 2477 tmp.setE( E ); aNucleon->SetMoment 2478 SumMasses += E; 2479 } 2480 } 2481 2482 G4double Chigh = Mass / SumMasses; G4do 2483 const G4int maxNumberOfLoops = 1000; 2484 G4int loopCounter = 0; 2485 do { 2486 C = ( Chigh + Clow ) / 2.0; 2487 2488 SumMasses = 0.0; 2489 aNucleon = 0; 2490 theProjectileNucleus->StartLoop(); 2491 while ( ( aNucleon = theProjectileNuc 2492 if ( ! aNucleon->AreYouHit() ) { 2493 G4LorentzVector tmp = aNucleon->G 2494 G4double E = std::sqrt( tmp.vect( 2495 sqr( aNuc 2496 SumMasses += E; 2497 } 2498 } 2499 2500 if ( SumMasses > Mass) Chigh = C; 2501 else Clow = C; 2502 2503 } while ( Chigh - Clow > 0.01 && 2504 ++loopCounter < maxNumberOfLo 2505 if ( loopCounter >= maxNumberOfLoops ) 2506 #ifdef debugFTFmodel 2507 G4cout << "BAD situation: forced exit 2508 << "\t return immediately from 2509 #endif 2510 return; 2511 } 2512 2513 aNucleon = 0; 2514 theProjectileNucleus->StartLoop(); 2515 while ( ( aNucleon = theProjectileNucle 2516 if ( ! aNucleon->AreYouHit() ) { 2517 G4LorentzVector tmp = aNucleon->Get 2518 G4double E = std::sqrt( tmp.vect(). 2519 sqr( aNucle 2520 tmp.setE( E ); tmp.boost( -bstToCM 2521 aNucleon->SetMomentum( tmp ); 2522 } 2523 } 2524 } // End of if ( ProjectileResidualMass 2525 2526 #ifdef debugFTFmodel 2527 G4cout << "End projectile" << G4endl; 2528 #endif 2529 2530 } else { // Related to the condition: if ( 2531 2532 #ifdef debugFTFmodel 2533 G4cout << "Low energy interaction: Target 2534 << "Tr ResidualMassNumber Tr Resid 2535 << TargetResidualMassNumber << " " 2536 << TargetResidualExcitationEnergy 2537 #endif 2538 2539 G4int NumberOfTargetParticipant( 0 ); 2540 for ( G4int i = 0; i < NumberOfInvolvedNu 2541 G4Nucleon* aNucleon = TheInvolvedNucleo 2542 G4VSplitableHadron* targetSplitable = a 2543 if ( targetSplitable->GetSoftCollisionC 2544 } 2545 2546 G4double DeltaExcitationE( 0.0 ); 2547 G4LorentzVector DeltaPResidualNucleus = G 2548 2549 if ( NumberOfTargetParticipant != 0 ) { 2550 DeltaExcitationE = TargetResidualExcita 2551 DeltaPResidualNucleus = TargetResidual4 2552 } 2553 2554 for ( G4int i = 0; i < NumberOfInvolvedNu 2555 G4Nucleon* aNucleon = TheInvolvedNucleo 2556 G4VSplitableHadron* targetSplitable = a 2557 if ( targetSplitable->GetSoftCollisionC 2558 G4LorentzVector tmp = -DeltaPResidual 2559 aNucleon->SetMomentum( tmp ); 2560 aNucleon->SetBindingEnergy( DeltaExci 2561 } else { 2562 delete targetSplitable; 2563 targetSplitable = 0; 2564 aNucleon->Hit( targetSplitable ); 2565 aNucleon->SetBindingEnergy( 0.0 ); 2566 } 2567 } 2568 2569 #ifdef debugFTFmodel 2570 G4cout << "NumberOfTargetParticipant " << 2571 << "TargetResidual4Momentum " << 2572 #endif 2573 2574 if ( ! GetProjectileNucleus() ) return; 2575 2576 #ifdef debugFTFmodel 2577 G4cout << "Low energy interaction: Projec 2578 << "Pr ResidualMassNumber Pr Resid 2579 << ProjectileResidualMassNumber << 2580 << ProjectileResidualExcitationEne 2581 #endif 2582 2583 G4int NumberOfProjectileParticipant( 0 ); 2584 for ( G4int i = 0; i < NumberOfInvolvedNu 2585 G4Nucleon* aNucleon = TheInvolvedNucleo 2586 G4VSplitableHadron* projectileSplitable 2587 if ( projectileSplitable->GetSoftCollis 2588 } 2589 2590 #ifdef debugFTFmodel 2591 G4cout << "NumberOfProjectileParticipant" 2592 #endif 2593 2594 DeltaExcitationE = 0.0; 2595 DeltaPResidualNucleus = G4LorentzVector( 2596 2597 if ( NumberOfProjectileParticipant != 0 ) 2598 DeltaExcitationE = ProjectileResidualEx 2599 DeltaPResidualNucleus = ProjectileResid 2600 } 2601 //G4cout << "DeltaExcitationE DeltaPResid 2602 // << " " << DeltaPResidualNucleus 2603 for ( G4int i = 0; i < NumberOfInvolvedNu 2604 G4Nucleon* aNucleon = TheInvolvedNucleo 2605 G4VSplitableHadron* projectileSplitable 2606 if ( projectileSplitable->GetSoftCollis 2607 G4LorentzVector tmp = -DeltaPResidual 2608 aNucleon->SetMomentum( tmp ); 2609 aNucleon->SetBindingEnergy( DeltaExci 2610 } else { 2611 delete projectileSplitable; 2612 projectileSplitable = 0; 2613 aNucleon->Hit( projectileSplitable ); 2614 aNucleon->SetBindingEnergy( 0.0 ); 2615 } 2616 } 2617 2618 #ifdef debugFTFmodel 2619 G4cout << "NumberOfProjectileParticipant 2620 << "ProjectileResidual4Momentum " 2621 #endif 2622 2623 } // End of the condition: if ( HighEnergy 2624 2625 #ifdef debugFTFmodel 2626 G4cout << "End GetResiduals --------------- 2627 #endif 2628 2629 } 2630 2631 2632 //=========================================== 2633 2634 G4ThreeVector G4FTFModel::GaussianPt( G4doubl 2635 2636 G4double Pt2( 0.0 ), Pt( 0.0 ); 2637 2638 if (AveragePt2 > 0.0) { 2639 const G4double ymax = maxPtSquare/Average 2640 if ( ymax < 200. ) { 2641 Pt2 = -AveragePt2 * G4Log( 1.0 + G4Unif 2642 } else { 2643 Pt2 = -AveragePt2 * G4Log( 1.0 - G4Unif 2644 } 2645 Pt = std::sqrt( Pt2 ); 2646 } 2647 2648 G4double phi = G4UniformRand() * twopi; 2649 2650 return G4ThreeVector( Pt*std::cos(phi), Pt* 2651 } 2652 2653 //=========================================== 2654 2655 G4bool G4FTFModel:: 2656 ComputeNucleusProperties( G4V3DNucleus* nucle 2657 G4LorentzVector& nu 2658 G4LorentzVector& re 2659 G4double& sumMasses 2660 G4double& residualE 2661 G4double& residualM 2662 G4int& residualMass 2663 G4int& residualChar 2664 2665 // This method, which is called only by Put 2666 // - either the target nucleus (which is n 2667 // of hadronic interaction (hadron-nucle 2668 // - or the projectile nucleus or antinucl 2669 // or antinucleus-nucleus interaction. 2670 // This method assumes that the all the par 2671 // the action of this method consists in mo 2672 // first one. The return value is "false" o 2673 // is null. 2674 2675 if ( ! nucleus ) return false; 2676 2677 G4double ExcitationEnergyPerWoundedNucleon 2678 theParameters->GetExcitationEnergyPerWoun 2679 2680 // Loop over the nucleons of the nucleus. 2681 // The nucleons that have been involved in 2682 // Reggeon Cascading) will be candidate to 2683 // All the remaining nucleons will be the n 2684 // The variable sumMasses is the amount of 2685 // 1. transverse mass of each involved 2686 // 2. 20.0*MeV separation energy for ea 2687 // 3. transverse mass of the residual n 2688 // In this first evaluation of sumMasses, t 2689 // (residualExcitationEnergy, estimated by 2690 // nucleon) is not taken into account. 2691 G4int residualNumberOfLambdas = 0; // Proj 2692 G4Nucleon* aNucleon = 0; 2693 nucleus->StartLoop(); 2694 while ( ( aNucleon = nucleus->GetNextNucleo 2695 nucleusMomentum += aNucleon->Get4Momentum 2696 if ( aNucleon->AreYouHit() ) { // Involv 2697 // Consider in sumMasses the nominal, i 2698 // (not the current masses, which could 2699 sumMasses += std::sqrt( sqr( aNucleon-> 2700 + aNucleon->Ge 2701 sumMasses += 20.0*MeV; // Separation e 2702 2703 //residualExcitationEnergy += Excitatio 2704 residualExcitationEnergy += -Excitation 2705 2706 residualMassNumber--; 2707 // The absolute value below is needed o 2708 residualCharge -= std::abs( G4int( aNuc 2709 } else { // Spectator nucleons 2710 residualMomentum += aNucleon->Get4Momen 2711 if ( aNucleon->GetDefinition() == G4Lam 2712 aNucleon->GetDefinition() == G4AntiLambd 2713 ++residualNumberOfLambdas; 2714 } 2715 } 2716 } 2717 #ifdef debugPutOnMassShell 2718 G4cout << "ExcitationEnergyPerWoundedNucleo 2719 << "\t Residual Charge, MassNumber ( 2720 << residualMassNumber << " (" << residualN 2721 << G4endl << "\t Initial Momentum " 2722 << G4endl << "\t Residual Momentum 2723 #endif 2724 residualMomentum.setPz( 0.0 ); 2725 residualMomentum.setE( 0.0 ); 2726 if ( residualMassNumber == 0 ) { 2727 residualMass = 0.0; 2728 residualExcitationEnergy = 0.0; 2729 } else { 2730 if ( residualMassNumber == 1 ) { 2731 if ( std::abs( residualCharge ) == 1 ) 2732 residualMass = G4Proton::Definition() 2733 } else if ( residualNumberOfLambdas == 2734 residualMass = G4Lambda::Definition() 2735 } else { 2736 residualMass = G4Neutron::Definition( 2737 } 2738 residualExcitationEnergy = 0.0; 2739 } else { 2740 if ( residualNumberOfLambdas > 0 ) { 2741 if ( residualMassNumber == 2 ) { 2742 residualMass = G4Lambda::Definition()->Ge 2743 if ( std::abs( residualCharge ) == 2744 residualMass += G4Proton::Definit 2745 } else if ( residualNumberOfLambdas == 1 2746 residualMass += G4Neutron::Definition() 2747 } else { 2748 residualMass += G4Lambda::Definition()- 2749 } 2750 } else { 2751 residualMass = G4HyperNucleiProperties::G 2752 residualNumberOfLambdas ) 2753 } 2754 } else { 2755 residualMass = G4ParticleTable::GetPa 2756 GetIonMass( std::abs( residu 2757 } 2758 } 2759 residualMass += residualExcitationEnergy; 2760 } 2761 sumMasses += std::sqrt( sqr( residualMass ) 2762 return true; 2763 } 2764 2765 2766 //=========================================== 2767 2768 G4bool G4FTFModel:: 2769 GenerateDeltaIsobar( const G4double sqrtS, 2770 const G4int numberOfInvo 2771 G4Nucleon* involvedNucle 2772 G4double& sumMasses ) { 2773 2774 // This method, which is called only by Put 2775 // re-interpret some of the involved nucleo 2776 // - either by replacing a proton (2212) wi 2777 // - or by replacing a neutron (2112) with 2778 // The on-shell mass of these delta-isobars 2779 // the corresponding nucleon on-shell mass. 2780 // the max number of deltas compatible with 2781 // The delta-isobars are considered with th 2782 // corresponding nucleons. 2783 // This method assumes that all the paramet 2784 // the action of this method consists in mo 2785 // sumMasses. The return value is "false" o 2786 // have unphysical values. 2787 2788 if ( sqrtS < 0.0 || numberOfInvolvedNucle 2789 2790 const G4double probDeltaIsobar = 0.05; 2791 2792 G4int maxNumberOfDeltas = G4int( (sqrtS - s 2793 G4int numberOfDeltas = 0; 2794 2795 for ( G4int i = 0; i < numberOfInvolvedNucl 2796 2797 if ( G4UniformRand() < probDeltaIsobar & 2798 numberOfDeltas++; 2799 if ( ! involvedNucleons[i] ) continue; 2800 // Skip any eventual lambda (that can b 2801 if ( involvedNucleons[i]->GetDefinition 2802 involvedNucleons[i]->GetDefinition() == 2803 G4VSplitableHadron* splitableHadron = i 2804 G4double massNuc = std::sqrt( sqr( spli 2805 + splitab 2806 // The absolute value below is needed i 2807 G4int pdgCode = std::abs( splitableHadr 2808 const G4ParticleDefinition* old_def = s 2809 G4int newPdgCode = pdgCode/10; newPdgCo 2810 if ( splitableHadron->GetDefinition()-> 2811 const G4ParticleDefinition* ptr = 2812 G4ParticleTable::GetParticleTable()-> 2813 splitableHadron->SetDefinition( ptr ); 2814 G4double massDelta = std::sqrt( sqr( sp 2815 + split 2816 //G4cout << i << " " << sqrtS/GeV << " 2817 // << " " << massNuc << G4endl; 2818 if ( sqrtS < sumMasses + massDelta - ma 2819 splitableHadron->SetDefinition( old_d 2820 break; 2821 } else { // Change is accepted 2822 sumMasses += ( massDelta - massNuc ); 2823 } 2824 } 2825 } 2826 2827 return true; 2828 } 2829 2830 2831 //=========================================== 2832 2833 G4bool G4FTFModel:: 2834 SamplingNucleonKinematics( G4double averagePt 2835 const G4double max 2836 G4double dCor, 2837 G4V3DNucleus* nucl 2838 const G4LorentzVec 2839 const G4double res 2840 const G4int residu 2841 const G4int number 2842 G4Nucleon* involve 2843 G4double& mass2 ) 2844 2845 // This method, which is called only by Put 2846 // - either the target nucleons: this for 2847 // (hadron-nucleus, nucleus-nucleus, ant 2848 // - or the projectile nucleons or antinuc 2849 // nucleus-nucleus or antinucleus-nucleu 2850 // This method assumes that all the paramet 2851 // the action of this method consists in ch 2852 // whose pointers are in the vector involve 2853 // variable mass2. 2854 #ifdef debugPutOnMassShell 2855 G4cout << "G4FTFModel::SamplingNucleonKinem 2856 G4cout << " averagePt2= " << averagePt2 << 2857 << " dCor= " << dCor << " resMass(GeV)= " 2858 << " resMassN= " << residualMassNumber 2859 << " nNuc= " << numberOfInvolvedNucl 2860 << " lv= " << pResidual << G4endl; 2861 #endif 2862 2863 if ( ! nucleus || numberOfInvolvedNucleon 2864 2865 if ( residualMassNumber == 0 && numberOfI 2866 dCor = 0.0; 2867 averagePt2 = 0.0; 2868 } 2869 2870 G4bool success = true; 2871 2872 G4double SumMasses = residualMass; 2873 G4double invN = 1.0 / (G4double)numberOfInv 2874 2875 // to avoid problems due to precision lost 2876 const G4double eps = 1.e-10; 2877 const G4int maxNumberOfLoops = 1000; 2878 G4int loopCounter = 0; 2879 do { 2880 2881 success = true; 2882 2883 // Sampling of nucleon Pt 2884 G4ThreeVector ptSum( 0.0, 0.0, 0.0 ); 2885 if( averagePt2 > 0.0 ) { 2886 for ( G4int i = 0; i < numberOfInvolved 2887 G4Nucleon* aNucleon = involvedNucleons[i]; 2888 if ( ! aNucleon ) continue; 2889 G4ThreeVector tmpPt = GaussianPt( averagePt 2890 ptSum += tmpPt; 2891 G4LorentzVector tmp( tmpPt.x(), tmpPt.y(), 2892 aNucleon->SetMomentum( tmp ); 2893 } 2894 } 2895 2896 G4double deltaPx = ( ptSum.x() - pResidua 2897 G4double deltaPy = ( ptSum.y() - pResidua 2898 2899 SumMasses = residualMass; 2900 for ( G4int i = 0; i < numberOfInvolvedNu 2901 G4Nucleon* aNucleon = involvedNucleons[ 2902 if ( ! aNucleon ) continue; 2903 G4double px = aNucleon->Get4Momentum(). 2904 G4double py = aNucleon->Get4Momentum(). 2905 G4double MtN = std::sqrt( sqr( aNucleon 2906 + sqr( px ) + 2907 SumMasses += MtN; 2908 G4LorentzVector tmp( px, py, 0.0, MtN); 2909 aNucleon->SetMomentum( tmp ); 2910 } 2911 2912 // Sampling X of nucleon 2913 G4double xSum = 0.0; 2914 2915 for ( G4int i = 0; i < numberOfInvolvedNu 2916 G4Nucleon* aNucleon = involvedNucleons[ 2917 if ( ! aNucleon ) continue; 2918 2919 G4double x = 0.0; 2920 if( 0.0 != dCor ) { 2921 G4ThreeVector tmpX = GaussianPt( dCor*dCor, 2922 x = tmpX.x(); 2923 } 2924 x += aNucleon->Get4Momentum().e()/SumMa 2925 if ( x < -eps || x > 1.0 + eps ) { 2926 success = false; 2927 break; 2928 } 2929 x = std::min(1.0, std::max(x, 0.0)); 2930 xSum += x; 2931 // The energy is in the lab (instead of 2932 2933 G4LorentzVector tmp( aNucleon->Get4Mome 2934 aNucleon->Get4Mome 2935 x, aNucleon->Get4M 2936 aNucleon->SetMomentum( tmp ); 2937 } 2938 2939 if ( xSum < -eps || xSum > 1.0 + eps ) su 2940 if ( ! success ) continue; 2941 2942 G4double delta = ( residualMassNumber == 2943 2944 xSum = 1.0; 2945 mass2 = 0.0; 2946 for ( G4int i = 0; i < numberOfInvolvedNu 2947 G4Nucleon* aNucleon = involvedNucleons[ 2948 if ( ! aNucleon ) continue; 2949 G4double x = aNucleon->Get4Momentum().p 2950 xSum -= x; 2951 2952 if ( residualMassNumber == 0 ) { 2953 if ( x <= -eps || x > 1.0 + eps ) { 2954 success = false; 2955 break; 2956 } 2957 } else { 2958 if ( x <= -eps || x > 1.0 + eps || xS 2959 success = false; 2960 break; 2961 } 2962 } 2963 x = std::min( 1.0, std::max(x, eps) ); 2964 2965 mass2 += sqr( aNucleon->Get4Momentum(). 2966 2967 G4LorentzVector tmp( aNucleon->Get4Mome 2968 x, aNucleon->Get4M 2969 aNucleon->SetMomentum( tmp ); 2970 } 2971 if ( ! success ) continue; 2972 xSum = std::min( 1.0, std::max(xSum, eps) 2973 2974 if ( residualMassNumber > 0 ) mass2 += ( 2975 2976 #ifdef debugPutOnMassShell 2977 G4cout << "success: " << success << " Mt( 2978 << std::sqrt( mass2 )/GeV << G4endl; 2979 #endif 2980 2981 } while ( ( ! success ) && 2982 ++loopCounter < maxNumberOfLoops 2983 return ( loopCounter < maxNumberOfLoops ); 2984 } 2985 2986 2987 //=========================================== 2988 2989 G4bool G4FTFModel:: 2990 CheckKinematics( const G4double sValue, 2991 const G4double sqrtS, 2992 const G4double projectileMas 2993 const G4double targetMass2, 2994 const G4double nucleusY, 2995 const G4bool isProjectileNuc 2996 const G4int numberOfInvolved 2997 G4Nucleon* involvedNucleons[ 2998 G4double& targetWminus, 2999 G4double& projectileWplus, 3000 G4bool& success ) { 3001 3002 // This method, which is called only by Put 3003 // kinematics is acceptable or not. 3004 // This method assumes that all the paramet 3005 // notice that the input boolean parameter 3006 // only in the case of nucleus or antinucle 3007 // The action of this method consists in co 3008 // and setting the parameter success to fal 3009 // be rejeted. 3010 3011 G4double decayMomentum2 = sqr( sValue ) + s 3012 - 2.0*( sValue*( 3013 + project 3014 targetWminus = ( sValue - projectileMass2 + 3015 / 2.0 / sqrtS; 3016 projectileWplus = sqrtS - targetMass2/targe 3017 G4double projectilePz = projectileWplus/2.0 3018 G4double projectileE = projectileWplus/2.0 3019 G4double projectileY = 0.5 * G4Log( (proje 3020 (proje 3021 G4double targetPz = -targetWminus/2.0 + tar 3022 G4double targetE = targetWminus/2.0 + tar 3023 G4double targetY = 0.5 * G4Log( (targetE + 3024 3025 #ifdef debugPutOnMassShell 3026 G4cout << "decayMomentum2 " << decayMomentu 3027 << "\t targetWminus projectileWplus 3028 << "\t projectileY targetY " << proj 3029 if ( isProjectileNucleus ) { 3030 G4cout << "Order# of Wounded nucleon i, n 3031 } else { 3032 G4cout << "Order# of Wounded nucleon i, n 3033 } 3034 G4cout << G4endl; 3035 #endif 3036 3037 for ( G4int i = 0; i < numberOfInvolvedNucl 3038 G4Nucleon* aNucleon = involvedNucleons[i] 3039 if ( ! aNucleon ) continue; 3040 G4LorentzVector tmp = aNucleon->Get4Momen 3041 G4double mt2 = sqr( tmp.x() ) + sqr( tmp. 3042 sqr( aNucleon->GetSplitabl 3043 G4double x = tmp.z(); 3044 G4double pz = -targetWminus*x/2.0 + mt2/( 3045 G4double e = targetWminus*x/2.0 + mt2/( 3046 if ( isProjectileNucleus ) { 3047 pz = projectileWplus*x/2.0 - mt2/(2.0*p 3048 e = projectileWplus*x/2.0 + mt2/(2.0*p 3049 } 3050 G4double nucleonY = 0.5 * G4Log( (e + pz) 3051 3052 #ifdef debugPutOnMassShell 3053 if( isProjectileNucleus ) { 3054 G4cout << " " << i << " " << nucleonY < 3055 } else { 3056 G4cout << " " << i << " " << nucleonY < 3057 } 3058 G4cout << G4endl; 3059 #endif 3060 3061 if ( std::abs( nucleonY - nucleusY ) > 2 3062 ( isProjectileNucleus && targetY > 3063 ( ! isProjectileNucleus && project 3064 success = false; 3065 break; 3066 } 3067 } 3068 return true; 3069 } 3070 3071 3072 //=========================================== 3073 3074 G4bool G4FTFModel:: 3075 FinalizeKinematics( const G4double w, 3076 const G4bool isProjectile 3077 const G4LorentzRotation& 3078 const G4double residualMa 3079 const G4int residualMassN 3080 const G4int numberOfInvol 3081 G4Nucleon* involvedNucleo 3082 G4LorentzVector& residual4Momen 3083 3084 // This method, which is called only by Put 3085 // this method is called when we are sure t 3086 // acceptable. 3087 // This method assumes that all the paramet 3088 // notice that the input boolean parameter 3089 // only in the case of nucleus or antinucle 3090 // because the sign of pz (in the center-of 3091 // with respect to the case of a normal had 3092 // The action of this method consists in mo 3093 // (in the lab frame) and computing the res 3094 // frame). 3095 3096 G4ThreeVector residual3Momentum( 0.0, 0.0, 3097 3098 for ( G4int i = 0; i < numberOfInvolvedNucl 3099 G4Nucleon* aNucleon = involvedNucleons[i] 3100 if ( ! aNucleon ) continue; 3101 G4LorentzVector tmp = aNucleon->Get4Momen 3102 residual3Momentum -= tmp.vect(); 3103 G4double mt2 = sqr( tmp.x() ) + sqr( tmp. 3104 sqr( aNucleon->GetSplitabl 3105 G4double x = tmp.z(); 3106 G4double pz = -w * x / 2.0 + mt2 / ( 2. 3107 G4double e = w * x / 2.0 + mt2 / ( 2. 3108 // Reverse the sign of pz in the case of 3109 if ( isProjectileNucleus ) pz *= -1.0; 3110 tmp.setPz( pz ); 3111 tmp.setE( e ); 3112 tmp.transform( boostFromCmsToLab ); 3113 aNucleon->SetMomentum( tmp ); 3114 G4VSplitableHadron* splitableHadron = aNu 3115 splitableHadron->Set4Momentum( tmp ); 3116 } 3117 3118 G4double residualMt2 = sqr( residualMass ) 3119 + sqr( residual3Moment 3120 3121 #ifdef debugPutOnMassShell 3122 if ( isProjectileNucleus ) { 3123 G4cout << "Wminus Proj and residual3Momen 3124 } else { 3125 G4cout << "Wplus Targ and residual3Momen 3126 } 3127 #endif 3128 3129 G4double residualPz = 0.0; 3130 G4double residualE = 0.0; 3131 if ( residualMassNumber != 0 ) { 3132 residualPz = -w * residual3Momentum.z() / 3133 residualMt2 / ( 2.0 * w * r 3134 residualE = w * residual3Momentum.z() / 3135 residualMt2 / ( 2.0 * w * r 3136 // Reverse the sign of residualPz in the 3137 if ( isProjectileNucleus ) residualPz *= 3138 } 3139 3140 residual4Momentum.setPx( residual3Momentum. 3141 residual4Momentum.setPy( residual3Momentum. 3142 residual4Momentum.setPz( residualPz ); 3143 residual4Momentum.setE( residualE ); 3144 3145 return true; 3146 } 3147 3148 3149 //=========================================== 3150 3151 void G4FTFModel::ModelDescription( std::ostre 3152 desc << " FTF (Fritiof) Mod 3153 << "The FTF model is based on the well 3154 << "model (B. Andersson et al., Nucl. 3155 << "(1987)). Its first program impleme 3156 << "by B. Nilsson-Almquist and E. Sten 3157 << "Comm. 43, 387 (1987)). The Fritiof 3158 << "that all hadron-hadron interaction 3159 << "reactions, h_1+h_2->h_1'+h_2' wher 3160 << "are excited states of the hadrons 3161 << "mass spectra. The excited hadrons 3162 << "QCD-strings, and the corresponding 3163 << "fragmentation model is applied for 3164 << "their decays. 3165 << " The Fritiof model assumes that 3166 << "a hadron-nucleus interaction a str 3167 << "from the projectile can interact w 3168 << "nuclear nucleons and becomes into 3169 << "states. The probability of multipl 3170 << "calculated in the Glauber approxim 3171 << "of secondary particles was neglect 3172 << "to these, the original Fritiof mod 3173 << "cribe a nuclear destruction and sl 3174 << " In order to overcome the diffic 3175 << "the model by the reggeon theory in 3176 << "nuclear desctruction (Kh. Abdel-Wa 3177 << "nsky, Phys. Atom. Nucl. 60, 828 (1 3178 << "(1997)). Momenta of the nucleons e 3179 << "leus in the reggeon cascading are 3180 << "to a Fermi motion algorithm presen 3181 << "Collaboration (M.I. Adamovich et a 3182 << "A358, 337 (1997)). 3183 << " New features were also added to 3184 << "implemented in Geant4: a simulatio 3185 << "ron-nucleon scatterings, a simulat 3186 << "reactions like NN>NN* in hadron-nu 3187 << "a separate simulation of single di 3188 << " diffractive events. These allowed 3189 << "model parameter tuning a wide set 3190 << "data. 3191 } 3192 3193