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 #include <utility> 27 28 #include "G4QGSParticipants.hh" 29 #include "globals.hh" 30 #include "G4SystemOfUnits.hh" 31 #include "G4LorentzVector.hh" 32 #include "G4V3DNucleus.hh" 33 #include "G4ParticleTable.hh" 34 #include "G4IonTable.hh" 35 #include "G4PhysicalConstants.hh" 36 37 #include "G4Exp.hh" 38 #include "G4Log.hh" 39 #include "G4Pow.hh" 40 41 //#define debugQGSParticipants 42 //#define debugPutOnMassShell 43 44 // Class G4QGSParticipants 45 46 // Promoting model parameters from local varia 47 G4ThreadLocal G4int G4QGSParticipants_NPart = 48 49 G4QGSParticipants::G4QGSParticipants() 50 : theDiffExcitaton(), ModelMode(SOFT), nCutM 51 ThresholdParameter(0.0*GeV), QGSMThreshold 52 theNucleonRadius(1.5*fermi), theCurrentVel 53 theProjectileSplitable(nullptr), Regge(nul 54 InteractionMode(ALL), alpha(-0.5), beta(2. 55 NumberOfInvolvedNucleonsOfTarget(0), Numbe 56 ProjectileResidual4Momentum(G4LorentzVecto 57 ProjectileResidualCharge(0), ProjectileRes 58 TargetResidual4Momentum(G4LorentzVector()) 59 TargetResidualCharge(0), TargetResidualExc 60 CofNuclearDestruction(0.0), R2ofNuclearDes 61 ExcitationEnergyPerWoundedNucleon(0.0), Do 62 Pt2ofNuclearDestruction(0.0), MaxPt2ofNucl 63 { 64 for (size_t i=0; i < 250; ++i) { 65 TheInvolvedNucleonsOfTarget[i] = nullptr; 66 TheInvolvedNucleonsOfProjectile[i] = nullp 67 } 68 // Parameters setting 69 SetCofNuclearDestruction( 1.0 ); 70 SetR2ofNuclearDestruction( 1.5*fermi*fermi ) 71 SetDofNuclearDestruction( 0.3 ); 72 SetPt2ofNuclearDestruction( 0.075*GeV*GeV ); 73 SetMaxPt2ofNuclearDestruction( 1.0*GeV*GeV ) 74 SetExcitationEnergyPerWoundedNucleon( 40.0*M 75 76 sigmaPt=0.25*sqr(GeV); 77 } 78 79 G4QGSParticipants::G4QGSParticipants(const G4Q 80 : G4VParticipants(), ModelMode(right.ModelMo 81 ThresholdParameter(right.ThresholdParamete 82 QGSMThreshold(right.QGSMThreshold), 83 theNucleonRadius(right.theNucleonRadius), 84 theCurrentVelocity(right.theCurrentVelocit 85 theProjectileSplitable(right.theProjectile 86 Regge(right.Regge), InteractionMode(right. 87 alpha(right.alpha), beta(right.beta), sigm 88 NumberOfInvolvedNucleonsOfTarget(right.Num 89 NumberOfInvolvedNucleonsOfProjectile(right 90 ProjectileResidual4Momentum(right.Projecti 91 ProjectileResidualMassNumber(right.Project 92 ProjectileResidualCharge(right.ProjectileR 93 ProjectileResidualExcitationEnergy(right.P 94 TargetResidual4Momentum(right.TargetResidu 95 TargetResidualMassNumber(right.TargetResid 96 TargetResidualCharge(right.TargetResidualC 97 TargetResidualExcitationEnergy(right.Targe 98 CofNuclearDestruction(right.CofNuclearDest 99 R2ofNuclearDestruction(right.R2ofNuclearDe 100 ExcitationEnergyPerWoundedNucleon(right.Ex 101 DofNuclearDestruction(right.DofNuclearDest 102 Pt2ofNuclearDestruction(right.Pt2ofNuclear 103 MaxPt2ofNuclearDestruction(right.MaxPt2ofN 104 { 105 for (size_t i=0; i < 250; ++i) { 106 TheInvolvedNucleonsOfTarget[i] = right.The 107 TheInvolvedNucleonsOfProjectile[i] = right 108 } 109 } 110 111 G4QGSParticipants::~G4QGSParticipants() {} 112 113 void G4QGSParticipants::BuildInteractions(cons 114 { 115 theProjectile = thePrimary; 116 117 Regge = new G4Reggeons(theProjectile.GetDefi 118 119 SetProjectileNucleus( 0 ); 120 121 NumberOfInvolvedNucleonsOfProjectile= 0; 122 G4LorentzVector tmp( 0.0, 0.0, 0.0, 0.0 ); 123 ProjectileResidualMassNumber = 0; 124 ProjectileResidualCharge = 0; 125 ProjectileResidualExcitationEnergy = 0.0; 126 ProjectileResidual4Momentum = tmp; 127 128 NumberOfInvolvedNucleonsOfTarget= 0; 129 TargetResidualMassNumber = theNucleus- 130 TargetResidualCharge = theNucleus- 131 TargetResidualExcitationEnergy = 0.0; 132 133 theNucleus->StartLoop(); 134 G4Nucleon* NuclearNucleon; 135 while ( ( NuclearNucleon = theNucleus->GetNe 136 tmp+=NuclearNucleon->Get4Momentum(); 137 } 138 TargetResidual4Momentum = tmp; 139 140 if ( std::abs( theProjectile.GetDefinition() 141 // Projectile is a hadron : meson or baryo 142 ProjectileResidualMassNumber = std::abs( t 143 ProjectileResidualCharge = G4int( theProje 144 ProjectileResidualExcitationEnergy = 0.0; 145 ProjectileResidual4Momentum.setVect( thePr 146 ProjectileResidual4Momentum.setE( theProje 147 } 148 149 #ifdef debugQGSParticipants 150 G4cout <<G4endl<< "G4QGSParticipants::Buil 151 << "thePrimary " << thePrimary.GetD 152 <<ProjectileResidual4Momentum<<G4en 153 G4cout << "Target A and Z " << theNucleus 154 << TargetResidual4Momentum<<G4endl; 155 #endif 156 157 G4bool Success( true ); 158 159 const G4int maxNumberOfLoops = 1000; 160 G4int loopCounter = 0; 161 do 162 { 163 const G4int maxNumberOfInternalLoops = 100 164 G4int internalLoopCounter = 0; 165 do 166 { 167 if(std::abs(theProjectile.GetDefinition( 168 { 169 SelectInteractions(theProjectile); // 170 } 171 else 172 { 173 GetList( theProjectile ); // Get list 174 } 175 176 if ( theInteractions.size() == 0 ) retur 177 178 StoreInvolvedNucleon(); // Store p 179 180 ReggeonCascade(); // Make re 181 182 Success = PutOnMassShell(); // Ascribe 183 184 if(!Success) PrepareInitialState( thePri 185 186 } while( (!Success) && ++internalLoopCount 187 188 if ( internalLoopCounter >= maxNumberOfInt 189 Success = false; 190 } 191 192 if ( Success ) { 193 #ifdef debugQGSParticipants 194 G4cout<<G4endl<<"PerformDiffractiveCol 195 #endif 196 197 PerformDiffractiveCollisions(); 198 199 #ifdef debugQGSParticipants 200 G4cout<<G4endl<<"SplitHadrons();" <<G4 201 #endif 202 203 SplitHadrons(); 204 205 if( theProjectileSplitable && theProject 206 #ifdef debugQGSParticipants 207 G4cout<<"Perform non-Diffractive Col 208 #endif 209 Success = DeterminePartonMomenta(); 210 } 211 212 if(!Success) PrepareInitialState( thePri 213 } 214 } while( (!Success) && ++loopCounter < maxNu 215 216 if ( loopCounter >= maxNumberOfLoops ) { 217 Success = false; 218 #ifdef debugQGSParticipants 219 G4cout<<"NOT Successful ======" <<G4endl 220 #endif 221 } 222 223 if ( Success ) { 224 CreateStrings(); // To create strings 225 226 GetResiduals(); // To calculate residual 227 228 #ifdef debugQGSParticipants 229 G4cout<<"All O.K. ======" <<G4endl; 230 #endif 231 } 232 233 // clean-up, if necessary 234 #ifdef debugQGSParticipants 235 G4cout<<"Clearing "<<G4endl; 236 G4cout<<"theInteractions.size() "<<theInte 237 #endif 238 239 if( Regge ) delete Regge; 240 241 std::for_each( theInteractions.begin(), theI 242 theInteractions.clear(); 243 244 // Erasing of target involved nucleons. 245 #ifdef debugQGSParticipants 246 G4cout<<"Erasing of involved target nucleo 247 #endif 248 249 if ( NumberOfInvolvedNucleonsOfTarget != 0 ) 250 for ( G4int i = 0; i < NumberOfInvolvedNu 251 G4VSplitableHadron* aNucleon = TheInvolv 252 if ( (aNucleon != 0 ) && (aNucleon->GetS 253 } 254 } 255 256 // Erasing of projectile involved nucleons. 257 if ( NumberOfInvolvedNucleonsOfProjectile != 258 for ( G4int i = 0; i < NumberOfInvolvedNu 259 G4VSplitableHadron* aNucleon = TheInvol 260 if ( aNucleon ) delete aNucleon; 261 } 262 } 263 264 #ifdef debugQGSParticipants 265 G4cout<<"Delition of target nucleons from 266 <<G4endl<<G4endl; 267 #endif 268 std::for_each(theTargets.begin(), theTargets 269 theTargets.clear(); 270 271 if ( theProjectileSplitable ) { 272 delete theProjectileSplitable; 273 theProjectileSplitable = 0; 274 } 275 } 276 277 //============================================ 278 void G4QGSParticipants::PrepareInitialState( c 279 { 280 // Clearing of the arrays 281 // Erasing of the projectile 282 G4InteractionContent* anIniteraction = theIn 283 G4VSplitableHadron* pProjectile = anIniterac 284 if( pProjectile ) delete pProjectile; 285 286 std::for_each(theInteractions.begin(), theIn 287 theInteractions.clear(); 288 289 // Erasing of the envolved nucleons and targ 290 theNucleus->StartLoop(); 291 G4Nucleon* aNucleon; 292 while ( ( aNucleon = theNucleus->GetNextNucl 293 { 294 if ( aNucleon->AreYouHit() ) { 295 G4VSplitableHadron* splaNucleon = aNucle 296 if ( (splaNucleon != 0) && (splaNucleon- 297 aNucleon->Hit(nullptr); 298 NumberOfInvolvedNucleonsOfTarget--; 299 } 300 } 301 302 // Erasing of nuclear nucleons participated 303 std::for_each(theTargets.begin(), theTargets 304 theTargets.clear(); 305 306 // Preparation to a new attempt 307 theProjectile = thePrimary; 308 309 theNucleus->Init(theNucleus->GetMassNumber() 310 theNucleus->SortNucleonsIncZ(); 311 DoLorentzBoost(-theCurrentVelocity); // Lor 312 313 if (theNucleus->GetMassNumber() == 1) 314 { 315 G4ThreeVector aPos = G4ThreeVector(0.,0.,0 316 theNucleus->StartLoop(); 317 G4Nucleon* tNucleon=theNucleus->GetNextNuc 318 tNucleon->SetPosition(aPos); 319 } 320 321 G4LorentzVector Tmp( 0.0, 0.0, 0.0, 0.0 ); 322 NumberOfInvolvedNucleonsOfTarget= 0; 323 TargetResidualMassNumber = theNucleus- 324 TargetResidualCharge = theNucleus- 325 TargetResidualExcitationEnergy = 0.0; 326 327 G4Nucleon* NuclearNucleon; 328 while ( ( NuclearNucleon = theNucleus->GetNe 329 {Tmp+=NuclearNucleon->Get4Moment 330 331 TargetResidual4Momentum = Tmp; 332 } 333 334 //============================================ 335 void G4QGSParticipants::GetList( const G4React 336 #ifdef debugQGSParticipants 337 G4cout<<G4endl<<"G4QGSParticipants::GetLis 338 #endif 339 340 // Direction: True - Proj, False - Target 341 theProjectileSplitable = new G4QGSMSplitable 342 theProjectileSplitable->SetStatus(1); 343 344 G4LorentzVector aPrimaryMomentum(thePrimary. 345 G4LorentzVector aNucleonMomentum(0.,0.,0., 9 346 347 G4double SS=(aPrimaryMomentum + aNucleonMome 348 349 Regge->SetS(SS); 350 351 //-------------------------------------- 352 theNucleus->StartLoop(); 353 G4Nucleon * tNucleon = theNucleus->GetNextNu 354 355 if ( ! tNucleon ) { 356 #ifdef debugQGSParticipants 357 G4cout << "QGSM - BAD situation: pNucleo 358 #endif 359 return; 360 } 361 362 G4double theNucleusOuterR = theNucleus->GetO 363 364 if (theNucleus->GetMassNumber() == 1) 365 { 366 G4ThreeVector aPos = G4ThreeVector(0.,0.,0 367 tNucleon->SetPosition(aPos); 368 theNucleusOuterR = 0.; 369 } 370 371 // Determination of participating nucleons o 372 373 std::for_each(theInteractions.begin(), theIn 374 theInteractions.clear(); 375 376 G4int MaxPower=thePrimary.GetMomentum().mag( 377 378 const G4int maxNumberOfLoops = 1000; 379 380 G4int NumberOfTries = 0; 381 G4double Scale = 1.0; 382 383 G4int loopCounter = -1; 384 while( (theInteractions.size() == 0) && ++lo 385 { 386 InteractionMode = ALL; // Mode = ALL, WIT 387 388 // choose random impact parameter of a col 389 std::pair<G4double, G4double> theImpactPar 390 391 NumberOfTries++; 392 if( NumberOfTries == 100*(NumberOfTries/10 393 394 theImpactParameter = theNucleus->ChooseImp 395 G4double impactX = theImpactParameter.firs 396 G4double impactY = theImpactParameter.seco 397 398 #ifdef debugQGSParticipants 399 G4cout<<"InteractionMode "<<InteractionM 400 G4cout<<"Impact parameter (fm ) "<<std:: 401 G4int nucleonCount = -1; 402 #endif 403 404 // loop over nucleons to find collisions 405 theNucleus->StartLoop(); 406 G4QGSParticipants_NPart = 0; 407 408 G4double Power=MaxPower; 409 410 while( (tNucleon = theNucleus->GetNextNucl 411 { 412 if(Power <= 0.) break; 413 414 G4LorentzVector nucleonMomentum=tNucleon 415 416 G4double Distance2 = sqr(impactX - tNucl 417 sqr(impactY - tNucl 418 419 G4double Pint(0.); // 420 G4double Pprd(0.), Ptrd(0.), Pdd(0.); // 421 G4double Pnd (0.), Pnvr(0.); // 422 G4int NcutPomerons(0); // 423 424 Regge->GetProbabilities(std::sqrt(Distan 425 Pint, Pprd, Ptrd, Pdd, Pnd, Pnvr); 426 #ifdef debugQGSParticipants 427 nucleonCount++; 428 G4cout<<"Nucleon & its impact paramete 429 G4cout<<"Probability of interaction: 430 G4cout<<"Probability of PrD, TrD, DD: "<< 431 G4cout<<"Probability of NonDiff, QuarkExc.: 432 #endif 433 434 if (Pint > G4UniformRand()) 435 { // An inte 436 437 G4double rndNumber = G4UniformRand(); 438 G4int InteractionType(0); 439 440 if((InteractionMode==ALL)||(Interactio 441 { 442 if( rndNumber < Pprd ) { 443 else if( rndNumber < Pprd+Ptrd) { 444 else if( rndNumber < Pprd+Ptrd+Pdd) { 445 else if( rndNumber < Pprd+Ptrd+Pdd+Pnd ) { 446 NcutPomerons = Regge->ncPom 447 else {Intera 448 } 449 else // InteractionMode == NON_DIFF 450 { 451 InteractionMode = NON_DIFF; 452 if( rndNumber < Ptrd ) {Interact 453 else if( rndNumber < Ptrd + Pnd) {Interact 454 } 455 456 if( (InteractionType == NonD) && (Ncut 457 458 G4QGSParticipants_NPart ++; 459 G4QGSMSplitableHadron* aTargetSPB = ne 460 tNucleon->Hit(aTargetSPB); 461 462 #ifdef debugQGSParticipants 463 G4cout<<"An interaction is happend." 464 G4cout<<"Target nucleon - "<<nucleon 465 <<tNucleon->GetDefinition()->G 466 G4cout<<"Interaction type:"<<Interac 467 <<" (0 -PrD, 1 - TrD, 2 - DD, 468 G4cout<<"New Inter. mode:"<<Interac 469 <<" (0 -ALL, 1 - WITHOUT_R, 2 470 if( InteractionType == NonD ) 471 G4cout<<"Number of cutted pomerons 472 #endif 473 474 if((InteractionType == PrD) || (Intera 475 (InteractionType == Qexc)) 476 { // 477 #ifdef debugQGSParticipants 478 G4cout<<"Diffractive-like interact 479 #endif 480 481 G4InteractionContent * aInteraction 482 theProjectileSplitable->SetStatus(1* 483 484 aInteraction->SetTarget(aTargetSPB); 485 aInteraction->SetTargetNucleon(tNucl 486 aTargetSPB->SetCollisionCount(0); 487 aTargetSPB->SetStatus(1); 488 489 aInteraction->SetNumberOfDiffractive 490 aInteraction->SetNumberOfSoftCollisi 491 aInteraction->SetStatus(InteractionT 492 theInteractions.push_back(aInteracti 493 } 494 else 495 { // non 496 #ifdef debugQGSParticipants 497 G4cout<<"Non-diffractive interacti 498 #endif 499 500 G4int nCuts; 501 502 G4int Vncut=0; 503 for(nCuts = 0; nCuts < NcutPomerons; nCuts 504 { 505 if( G4UniformRand() < Power/MaxPower ){V 506 } 507 nCuts=Vncut; 508 509 if( nCuts == 0 ) {delete aTargetSPB; tNucl 510 511 #ifdef debugQGSParticipants 512 G4cout<<"Number of cuts in the int 513 #endif 514 515 aTargetSPB->IncrementCollisionCount(nCuts) 516 aTargetSPB->SetStatus(0); 517 theTargets.push_back(aTargetSPB); 518 519 theProjectileSplitable->IncrementCollision 520 theProjectileSplitable->SetStatus(0* 521 522 G4InteractionContent * aInteraction = new 523 aInteraction->SetTarget(aTargetSPB); 524 aInteraction->SetTargetNucleon(tNucl 525 aInteraction->SetNumberOfSoftCollisions(nC 526 aInteraction->SetStatus(InteractionT 527 theInteractions.push_back(aInteraction); 528 } 529 } // End of if (Pint > G4UniformRand( 530 } // End of while( (tNucleon = theNucl 531 532 #ifdef debugQGSParticipants 533 G4cout << G4endl<<"Number of wounded nuc 534 #endif 535 536 } // End of while( (theInteractions.size() 537 538 if ( loopCounter >= maxNumberOfLoops ) { 539 #ifdef debugQGSParticipants 540 G4cout <<"BAD situation: forced loop exi 541 #endif 542 // Perhaps there is something to set here. 543 // Decrease impact parameter ?? 544 // Select collisions with only diffraction 545 // Selecy only non-diffractive interaction 546 } 547 //------------------------------------------ 548 std::vector<G4InteractionContent*>::iterator 549 550 if( theInteractions.size() != 0) 551 { 552 if( InteractionMode == ALL ) // It can be 553 { // Only the 554 i = theInteractions.end()-1; 555 556 while ( theInteractions.size() != 1 ) 557 { 558 G4InteractionContent* anInteraction = 559 G4Nucleon * pNucleon = anInteraction-> 560 delete anInteraction->GetTarget(); 561 delete *i; 562 i=theInteractions.erase(i); 563 i--; 564 } 565 } 566 else 567 { // All quark 568 i = theInteractions.begin(); 569 while ( i != theInteractions.end() ) 570 { 571 G4InteractionContent* anInteraction = 572 573 if( anInteraction->GetStatus() == Qexc 574 { 575 G4Nucleon* aTargetNucleon = a 576 aTargetNucleon->Hit(nullptr); 577 578 delete anInteraction->GetTarget(); 579 delete *i; 580 i=theInteractions.erase(i); 581 } 582 else 583 { 584 i++; 585 } 586 } 587 } 588 } 589 } 590 591 //============================================ 592 void G4QGSParticipants::StoreInvolvedNucleon() 593 { //To store nucleons involved in the interact 594 595 NumberOfInvolvedNucleonsOfTarget = 0; 596 597 theNucleus->StartLoop(); 598 599 G4Nucleon* aNucleon; 600 while ( ( aNucleon = theNucleus->GetNextNucl 601 if ( aNucleon->AreYouHit() ) { 602 TheInvolvedNucleonsOfTarget[NumberOfInvo 603 NumberOfInvolvedNucleonsOfTarget++; 604 } 605 } 606 607 #ifdef debugQGSParticipants 608 G4cout << G4endl<<"G4QGSParticipants::Stor 609 <<"Stored # of wounded nucleons of 610 << NumberOfInvolvedNucleonsOfTarget 611 #endif 612 return; 613 } 614 615 //============================================ 616 617 void G4QGSParticipants::ReggeonCascade() 618 { // Implementation of the reggeon theory insp 619 #ifdef debugQGSParticipants 620 G4cout << G4endl<<"Reggeon cascading ..... 621 G4cout<<"C of nucl. desctruction "<<GetCof 622 <<" R2 "<<GetR2ofNuclearDestruction( 623 #endif 624 625 G4int InitNINt = NumberOfInvolvedNucleonsOfT 626 627 // Reggeon cascading in target nucleus 628 for ( G4int InvTN = 0; InvTN < InitNINt; Inv 629 G4Nucleon* aTargetNucleon = TheInvolvedNuc 630 631 G4double CreationTime = aTargetNucleon->Ge 632 633 G4double XofWoundedNucleon = aTargetNucleo 634 G4double YofWoundedNucleon = aTargetNucleo 635 636 G4V3DNucleus* theTargetNucleus = theNucleu 637 theTargetNucleus->StartLoop(); 638 639 #ifdef debugQGSParticipants 640 G4int TrgNuc=0; 641 #endif 642 G4Nucleon* Neighbour(0); 643 while ( ( Neighbour = theTargetNucleus->Ge 644 #ifdef debugQGSParticipants 645 TrgNuc++; 646 #endif 647 if ( ! Neighbour->AreYouHit() ) { 648 G4double impact2 = sqr( XofWoundedNucl 649 sqr( YofWoundedNucl 650 651 if ( G4UniformRand() < GetCofNuclearDe 652 G4Exp( -impact2 653 ) { 654 // The neighbour nucleon is involved 655 #ifdef debugQGSParticipants 656 G4cout<<"Target nucleon involved i 657 <<Neighbour->GetDefinition() 658 #endif 659 TheInvolvedNucleonsOfTarget[ NumberO 660 NumberOfInvolvedNucleonsOfTarget++; 661 662 G4QGSMSplitableHadron* targetSplitab 663 664 Neighbour->Hit( targetSplitable ); 665 targetSplitable->SetTimeOfCreation( 666 targetSplitable->SetStatus( 2 ); 667 targetSplitable->SetCollisionCount(0 668 669 G4InteractionContent * anInteraction 670 anInteraction->SetTarget(targetSplit 671 anInteraction->SetTargetNucleon(Neig 672 673 anInteraction->SetNumberOfDiffractiv 674 anInteraction->SetNumberOfSoftCollis 675 anInteraction->SetStatus(3); 676 theInteractions.push_back(anInteract 677 } 678 } 679 } 680 } 681 682 #ifdef debugQGSParticipants 683 G4cout <<"Number of new involved nucleons 684 #endif 685 return; 686 } 687 688 //============================================ 689 690 G4bool G4QGSParticipants::PutOnMassShell() { 691 692 G4bool isProjectileNucleus = false; 693 if ( GetProjectileNucleus() ) { 694 isProjectileNucleus = true; 695 } 696 697 #ifdef debugPutOnMassShell 698 G4cout <<G4endl<< "PutOnMassShell start .. 699 if ( isProjectileNucleus ) {G4cout << "Put 700 #endif 701 702 G4LorentzVector Pprojectile( theProjectile.G 703 if ( Pprojectile.z() < 0.0 ) { 704 return false; 705 } 706 707 G4bool isOk = true; 708 709 G4LorentzVector Ptarget( 0.0, 0.0, 0.0, 0.0 710 G4LorentzVector PtargetResidual( 0.0, 0.0, 0 711 G4double SumMasses = 0.0; 712 G4V3DNucleus* theTargetNucleus = GetTargetNu 713 G4double TargetResidualMass = 0.0; 714 715 #ifdef debugPutOnMassShell 716 G4cout << "Target : "; 717 #endif 718 719 isOk = ComputeNucleusProperties( theTargetNu 720 TargetResid 721 TargetResid 722 723 if ( ! isOk ) return false; 724 725 G4double Mprojectile = 0.0; 726 G4double M2projectile = 0.0; 727 G4LorentzVector Pproj( 0.0, 0.0, 0.0, 0.0 ); 728 G4LorentzVector PprojResidual( 0.0, 0.0, 0.0 729 G4V3DNucleus* thePrNucleus = GetProjectileNu 730 G4double PrResidualMass = 0.0; 731 732 if ( ! isProjectileNucleus ) { // hadron-nu 733 Mprojectile = Pprojectile.mag(); 734 M2projectile = Pprojectile.mag2(); 735 SumMasses += Mprojectile + 20.0*MeV; 736 } else { // nucleus-nucleus or antinucleus- 737 738 #ifdef debugPutOnMassShell 739 G4cout << "Projectile : "; 740 #endif 741 742 isOk = ComputeNucleusProperties( thePrNucl 743 Projectil 744 Projectil 745 if ( ! isOk ) return false; 746 } 747 748 G4LorentzVector Psum = Pprojectile + Ptarget 749 G4double SqrtS = Psum.mag(); 750 G4double S = Psum.mag2(); 751 752 #ifdef debugPutOnMassShell 753 G4cout << "Pproj "<<Pprojectile<<G4endl; 754 G4cout << "Ptarg "<<Ptarget<<G4endl; 755 G4cout << "Psum " << Psum/GeV << " GeV" << 756 << "SumMasses, PrResidualMass and T 757 << PrResidualMass/GeV << " " << Tar 758 G4cout << "Ptar res. "<<PtargetResidual<<G 759 #endif 760 761 if ( SqrtS < SumMasses ) { 762 return false; // It is impossible to simu 763 } 764 765 // Try to consider also the excitation energ 766 // possible, with the available energy; othe 767 768 G4double savedSumMasses = SumMasses; 769 if ( isProjectileNucleus ) { 770 SumMasses -= std::sqrt( sqr( PrResidualMas 771 SumMasses += std::sqrt( sqr( PrResidualMas 772 + PprojResidual.pe 773 } 774 SumMasses -= std::sqrt( sqr( TargetResidualM 775 SumMasses += std::sqrt( sqr( TargetResidualM 776 + PtargetResidual.pe 777 778 if ( SqrtS < SumMasses ) { 779 SumMasses = savedSumMasses; 780 if ( isProjectileNucleus ) { 781 ProjectileResidualExcitationEnergy = 0.0 782 } 783 TargetResidualExcitationEnergy = 0.0; 784 } 785 786 TargetResidualMass += TargetResidualExcitati 787 if ( isProjectileNucleus ) { 788 PrResidualMass += ProjectileResidualExcita 789 } 790 791 #ifdef debugPutOnMassShell 792 if ( isProjectileNucleus ) { 793 G4cout << "PrResidualMass ProjResidualEx 794 << ProjectileResidualExcitationEnergy < 795 } 796 G4cout << "TargetResidualMass TargetResidu 797 << TargetResidualExcitationEnergy < 798 << "Sum masses " << SumMasses/GeV < 799 #endif 800 801 // Sampling of nucleons what can transfer to 802 if ( isProjectileNucleus && thePrNucleus-> 803 isOk = GenerateDeltaIsobar( SqrtS, NumberO 804 TheInvolvedNuc 805 } 806 if ( theTargetNucleus->GetMassNumber() != 1 807 isOk = isOk && 808 GenerateDeltaIsobar( SqrtS, NumberO 809 TheInvolvedNuc 810 } 811 if ( ! isOk ) return false; 812 813 // Now we know that it is kinematically poss 814 // of the involved nucleons (or correspondin 815 // We have to sample the kinematical variabl 816 // of the final state. The sampled kinematic 817 // Notice that the sampling of the transvers 818 // Fermi motion. 819 820 // If target is nucleon - return ? 821 822 G4LorentzRotation toCms( -1*Psum.boostVector 823 G4LorentzVector Ptmp = toCms*Pprojectile; 824 if ( Ptmp.pz() <= 0.0 ) { // "String" movin 825 return false; 826 } 827 828 G4LorentzRotation toLab( toCms.inverse() ); 829 830 G4double YprojectileNucleus = 0.0; 831 if ( isProjectileNucleus ) { 832 Ptmp = toCms*Pproj; 833 YprojectileNucleus = Ptmp.rapidity(); 834 } 835 Ptmp = toCms*Ptarget; 836 G4double YtargetNucleus = Ptmp.rapidity(); 837 838 // Ascribing of the involved nucleons Pt and 839 G4double DcorP = 0.0; 840 if ( isProjectileNucleus ) { 841 DcorP = GetDofNuclearDestruction() / thePr 842 } 843 G4double DcorT = GetDofNuclearDestruct 844 G4double AveragePt2 = GetPt2ofNuclearDestru 845 G4double maxPtSquare = GetMaxPt2ofNuclearDes 846 847 #ifdef debugPutOnMassShell 848 if ( isProjectileNucleus ) { 849 G4cout << "Y projectileNucleus " << Ypro 850 } 851 G4cout << "Y targetNucleus " << Ytarge 852 << "Dcor " << GetDofNuclearDestruct 853 << " DcorP DcorT " << DcorP << " " 854 #endif 855 856 G4double M2proj = M2projectile; // Initiali 857 G4double WplusProjectile = 0.0; 858 G4double M2target = 0.0; 859 G4double WminusTarget = 0.0; 860 G4int NumberOfTries = 0; 861 G4double ScaleFactor = 1.0; 862 G4bool OuterSuccess = true; 863 864 const G4int maxNumberOfLoops = 1000; 865 G4int loopCounter = 0; 866 do { 867 G4double sqrtM2proj = 0.0, sqrtM2target = 868 OuterSuccess = true; 869 const G4int maxNumberOfTries = 1000; 870 do { 871 NumberOfTries++; 872 if ( NumberOfTries == 100*(NumberOfTries 873 // After many tries, it is convenient 874 // AveragePt2, so that the sampled mom 875 // involved nucleons (or corresponding delta 876 // it is more likely to satisfy the mo 877 ScaleFactor /= 2.0; 878 DcorP *= ScaleFactor; 879 DcorT *= ScaleFactor; 880 AveragePt2 *= ScaleFactor; 881 } 882 if ( isProjectileNucleus ) { 883 // Sampling of kinematical properties 884 isOk = SamplingNucleonKinematics( Aver 885 theP 886 PrRe 887 Numb 888 TheI 889 } 890 // Sampling of kinematical properties of 891 isOk = isOk && 892 SamplingNucleonKinematics( Averag 893 theTar 894 Target 895 Number 896 TheInv 897 898 if ( M2proj < 0.0 ) { 899 if( M2proj < -0.000001 ) { 900 G4ExceptionDescription ed; 901 ed << "Projectile " << theProjectile.GetDe 902 << " Target (Z,A)=(" << theTargetNucle 903 << ") M2proj=" << M2proj << " -> 904 G4Exception( "G4QGSParticipants::PutOnMass 905 "HAD_QGSPARTICIPANTS_002", JustWarn 906 } 907 M2proj = 0.0; 908 } 909 sqrtM2proj = std::sqrt( M2proj ); 910 if ( M2target < 0.0 ) { 911 G4ExceptionDescription ed; 912 ed << "Projectile " << theProjectile.G 913 << " Target (Z,A)=(" << theTargetN 914 << ") M2target=" << M2target << " 915 G4Exception( "G4QGSParticipants::PutOn 916 "HAD_QGSPARTICIPANTS_003", 917 M2target = 0.0; 918 }; 919 sqrtM2target = std::sqrt( M2target ); 920 921 #ifdef debugPutOnMassShell 922 G4cout << "SqrtS, Mp+Mt, Mp, Mt " << S 923 << ( sqrtM2proj + sqrtM2target 924 << sqrtM2proj/GeV << " " << sqr 925 #endif 926 927 if ( ! isOk ) return false; 928 } while ( ( SqrtS < ( sqrtM2proj + sqrtM2t 929 ++NumberOfTries < maxNumberOfTri 930 if ( NumberOfTries >= maxNumberOfTries ) { 931 return false; 932 } 933 if ( isProjectileNucleus ) { 934 isOk = CheckKinematics( S, SqrtS, M2proj 935 NumberOfInvolved 936 TheInvolvedNucle 937 WminusTarget, Wp 938 } 939 isOk = isOk && 940 CheckKinematics( S, SqrtS, M2proj, 941 NumberOfInvolvedNu 942 WminusTarget, Wplu 943 if ( ! isOk ) return false; 944 } while ( ( ! OuterSuccess ) && 945 ++loopCounter < maxNumberOfLoops ) 946 if ( loopCounter >= maxNumberOfLoops ) { 947 return false; 948 } 949 950 // Now the sampling is completed, and we can 951 // whole system. This is done first in the c 952 // to the lab frame. The transverse momentum 953 // the recoil of each hadron (nucleon or del 954 // to conserve (by construction) the transve 955 956 if ( ! isProjectileNucleus ) { // hadron-nu 957 958 G4double Pzprojectile = WplusProjectile/2. 959 G4double Eprojectile = WplusProjectile/2. 960 Pprojectile.setPz( Pzprojectile ); 961 Pprojectile.setE( Eprojectile ); 962 963 #ifdef debugPutOnMassShell 964 G4cout << "Proj after in CMS " << Pproje 965 #endif 966 967 Pprojectile.transform( toLab ); 968 theProjectile.SetMomentum( Pprojectile.vec 969 theProjectile.SetTotalEnergy( Pprojectile. 970 971 if ( theProjectileSplitable ) theProjectil 972 973 #ifdef debugPutOnMassShell 974 G4cout << "Final proj. mom in Lab. " <<t 975 <<t 976 #endif 977 978 } else { // nucleus-nucleus or antinucleus- 979 980 isOk = FinalizeKinematics( WplusProjectile 981 ProjectileResid 982 TheInvolvedNucl 983 984 #ifdef debugPutOnMassShell 985 G4cout << "Projectile Residual4Momentum 986 #endif 987 988 if ( ! isOk ) return false; 989 990 ProjectileResidual4Momentum.transform( toL 991 992 #ifdef debugPutOnMassShell 993 G4cout << "Projectile Residual4Momentum 994 #endif 995 996 } 997 998 isOk = FinalizeKinematics( WminusTarget, fal 999 TargetResidualMas 1000 TheInvolvedNucle 1001 1002 #ifdef debugPutOnMassShell 1003 G4cout << "Target Residual4Momentum in CM 1004 #endif 1005 1006 if ( ! isOk ) return false; 1007 1008 TargetResidual4Momentum.transform( toLab ); 1009 1010 #ifdef debugPutOnMassShell 1011 G4cout << "Target Residual4Momentum in La 1012 #endif 1013 1014 return true; 1015 1016 } 1017 1018 //=========================================== 1019 1020 G4ThreeVector G4QGSParticipants::GaussianPt( 1021 // @@ this method is used in FTFModel as w 1022 1023 G4double Pt2( 0.0 ), Pt(0.0); 1024 if ( AveragePt2 > 0.0 ) { 1025 G4double x = maxPtSquare/AveragePt2; 1026 Pt2 = (x < 200) ? 1027 -AveragePt2 * G4Log( 1.0 + G4UniformRan 1028 : -AveragePt2 * G4Log( 1.0 - G4UniformR 1029 Pt = std::sqrt( Pt2 ); 1030 } 1031 G4double phi = G4UniformRand() * twopi; 1032 1033 return G4ThreeVector( Pt*std::cos(phi), Pt* 1034 } 1035 //=========================================== 1036 1037 G4bool G4QGSParticipants:: 1038 ComputeNucleusProperties( G4V3DNucleus* nucle 1039 G4LorentzVector& nu 1040 G4LorentzVector& re 1041 G4double& sumMasses 1042 G4double& residualE 1043 G4double& residualM 1044 G4int& residualMass 1045 G4int& residualChar 1046 1047 // This method, which is called only by Put 1048 // - either the target nucleus (which is n 1049 // of hadronic interaction (hadron-nucle 1050 // - or the projectile nucleus or antinucl 1051 // or antinucleus-nucleus interaction. 1052 // This method assumes that the all the par 1053 // the action of this method consists in mo 1054 // first one. The return value is "false" o 1055 // is null. 1056 1057 if ( ! nucleus ) return false; 1058 1059 G4double ExcitationEPerWoundedNucleon = Get 1060 1061 // Loop over the nucleons of the nucleus. 1062 // The nucleons that have been involved in 1063 // Reggeon Cascading) will be candidate to 1064 // All the remaining nucleons will be the n 1065 // The variable sumMasses is the amount of 1066 // 1. transverse mass of each involved 1067 // 2. 20.0*MeV separation energy for ea 1068 // 3. transverse mass of the residual n 1069 // In this first evaluation of sumMasses, t 1070 // (residualExcitationEnergy, estimated by 1071 // nucleon) is not taken into account. 1072 G4Nucleon* aNucleon = 0; 1073 nucleus->StartLoop(); 1074 while ( ( aNucleon = nucleus->GetNextNucleo 1075 nucleusMomentum += aNucleon->Get4Momentum 1076 if ( aNucleon->AreYouHit() ) { // Involv 1077 // Consider in sumMasses the nominal, i 1078 // (not the current masses, which could 1079 1080 sumMasses += std::sqrt( sqr( aNucleon-> 1081 + aNucleon->Ge 1082 sumMasses += 20.0*MeV; // Separation e 1083 1084 //residualExcitationEnergy += Excitatio 1085 residualExcitationEnergy += -Excitation 1086 residualMassNumber--; 1087 // The absolute value below is needed o 1088 residualCharge -= std::abs( G4int( aNuc 1089 } else { // Spectator nucleons 1090 residualMomentum += aNucleon->Get4Momen 1091 } 1092 } 1093 #ifdef debugPutOnMassShell 1094 G4cout << "ExcitationEnergyPerWoundedNucl 1095 << "\t Residual Charge, MassNumber 1096 << G4endl << "\t Initial Momentum 1097 << G4endl << "\t Residual Momentum 1098 #endif 1099 1100 residualMomentum.setPz( 0.0 ); 1101 residualMomentum.setE( 0.0 ); 1102 if ( residualMassNumber == 0 ) { 1103 residualMass = 0.0; 1104 residualExcitationEnergy = 0.0; 1105 } else { 1106 residualMass = G4ParticleTable::GetPartic 1107 GetIonMass( residualChar 1108 if ( residualMassNumber == 1 ) { 1109 residualExcitationEnergy = 0.0; 1110 } 1111 } 1112 sumMasses += std::sqrt( sqr( residualMass ) 1113 return true; 1114 } 1115 1116 1117 //=========================================== 1118 1119 G4bool G4QGSParticipants:: 1120 GenerateDeltaIsobar( const G4double sqrtS, 1121 const G4int numberOfInvo 1122 G4Nucleon* involvedNucle 1123 G4double& sumMasses ) { 1124 1125 // This method, which is called only by Put 1126 // re-interpret some of the involved nucleo 1127 // - either by replacing a proton (2212) wi 1128 // - or by replacing a neutron (2112) with 1129 // The on-shell mass of these delta-isobars 1130 // the corresponding nucleon on-shell mass. 1131 // the max number of deltas compatible with 1132 // The delta-isobars are considered with th 1133 // corresponding nucleons. 1134 // This method assumes that all the paramet 1135 // the action of this method consists in mo 1136 // sumMasses. The return value is "false" o 1137 // have unphysical values. 1138 1139 if ( sqrtS < 0.0 || numberOfInvolvedNucle 1140 1141 //const G4double ProbDeltaIsobar = 0.05; / 1142 //const G4double ProbDeltaIsobar = 0.25; / 1143 const G4double probDeltaIsobar = 0.10; // 1144 1145 G4int maxNumberOfDeltas = G4int( (sqrtS - s 1146 G4int numberOfDeltas = 0; 1147 1148 for ( G4int i = 0; i < numberOfInvolvedNucl 1149 //G4cout << "i maxNumberOfDeltas probDelt 1150 // << " " << probDeltaIsobar << G4e 1151 if ( G4UniformRand() < probDeltaIsobar & 1152 numberOfDeltas++; 1153 if ( ! involvedNucleons[i] ) continue; 1154 G4VSplitableHadron* splitableHadron = i 1155 G4double massNuc = std::sqrt( sqr( spli 1156 + splitab 1157 //AR The absolute value below is needed 1158 G4int pdgCode = std::abs( splitableHadr 1159 const G4ParticleDefinition* old_def = s 1160 G4int newPdgCode = pdgCode/10; newPdgCo 1161 if ( splitableHadron->GetDefinition()-> 1162 const G4ParticleDefinition* ptr = 1163 G4ParticleTable::GetParticleTable()-> 1164 splitableHadron->SetDefinition( ptr ); 1165 G4double massDelta = std::sqrt( sqr( sp 1166 + split 1167 //G4cout << i << " " << sqrtS/GeV << " 1168 // << " " << massNuc << G4endl; 1169 if ( sqrtS < sumMasses + massDelta - ma 1170 splitableHadron->SetDefinition( old_d 1171 break; 1172 } else { // Change is accepted 1173 sumMasses += ( massDelta - massNuc ); 1174 } 1175 } 1176 } 1177 //G4cout << "maxNumberOfDeltas numberOfDelt 1178 // << numberOfDeltas << G4endl; 1179 return true; 1180 } 1181 1182 1183 //=========================================== 1184 1185 G4bool G4QGSParticipants:: 1186 SamplingNucleonKinematics( G4double averagePt 1187 const G4double max 1188 G4double dCor, 1189 G4V3DNucleus* nucl 1190 const G4LorentzVec 1191 const G4double res 1192 const G4int residu 1193 const G4int number 1194 G4Nucleon* involve 1195 G4double& mass2 ) 1196 1197 // This method, which is called only by Put 1198 // - either the target nucleons: this for 1199 // (hadron-nucleus, nucleus-nucleus, ant 1200 // - or the projectile nucleons or antinuc 1201 // nucleus-nucleus or antinucleus-nucleu 1202 // This method assumes that all the paramet 1203 // the action of this method consists in ch 1204 // whose pointers are in the vector involve 1205 // variable mass2. 1206 1207 if ( ! nucleus ) return false; 1208 1209 if ( residualMassNumber == 0 && numberOfI 1210 dCor = 0.0; 1211 averagePt2 = 0.0; 1212 } 1213 1214 G4bool success = true; 1215 1216 G4double SumMasses = residualMass; 1217 for ( G4int i = 0; i < numberOfInvolvedNucl 1218 G4Nucleon* aNucleon = involvedNucleons[i] 1219 if ( ! aNucleon ) continue; 1220 SumMasses += aNucleon->GetSplitableHadron 1221 } 1222 1223 const G4int maxNumberOfLoops = 1000; 1224 G4int loopCounter = 0; 1225 do { 1226 1227 success = true; 1228 G4ThreeVector ptSum( 0.0, 0.0, 0.0 ); 1229 G4double xSum = 0.0; 1230 1231 for ( G4int i = 0; i < numberOfInvolvedNu 1232 G4Nucleon* aNucleon = involvedNucleons[ 1233 if ( ! aNucleon ) continue; 1234 G4ThreeVector tmpPt = GaussianPt( avera 1235 ptSum += tmpPt; 1236 G4ThreeVector tmpX = GaussianPt( dCor*d 1237 G4double x = tmpX.x() + 1238 aNucleon->GetSplitableHadr 1239 if ( x < 0.0 || x > 1.0 ) { 1240 success = false; 1241 break; 1242 } 1243 xSum += x; 1244 //AR The energy is in the lab (instead 1245 G4LorentzVector tmp( tmpPt.x(), tmpPt.y 1246 aNucleon->SetMomentum( tmp ); 1247 } 1248 1249 if ( xSum < 0.0 || xSum > 1.0 ) success 1250 1251 if ( ! success ) continue; 1252 1253 G4double deltaPx = ( ptSum.x() - pResidua 1254 G4double deltaPy = ( ptSum.y() - pResidua 1255 G4double delta = 0.0; 1256 if ( residualMassNumber == 0 ) { 1257 delta = ( xSum - 1.0 ) / numberOfInvolv 1258 } else { 1259 delta = 0.0; 1260 } 1261 1262 xSum = 1.0; 1263 mass2 = 0.0; 1264 for ( G4int i = 0; i < numberOfInvolvedNu 1265 G4Nucleon* aNucleon = involvedNucleons[ 1266 if ( ! aNucleon ) continue; 1267 G4double x = aNucleon->Get4Momentum().p 1268 xSum -= x; 1269 if ( residualMassNumber == 0 ) { 1270 if ( x <= 0.0 || x > 1.0 ) { 1271 success = false; 1272 break; 1273 } 1274 } else { 1275 if ( x <= 0.0 || x > 1.0 || xSum 1276 success = false; 1277 break; 1278 } 1279 } 1280 G4double px = aNucleon->Get4Momentum(). 1281 G4double py = aNucleon->Get4Momentum(). 1282 mass2 += ( sqr( aNucleon->GetSplitableH 1283 + sqr( px ) + sqr( py ) ) 1284 G4LorentzVector tmp( px, py, x, aNucleo 1285 aNucleon->SetMomentum( tmp ); 1286 } 1287 1288 if ( success && residualMassNumber != 0 1289 mass2 += ( sqr( residualMass ) + pResid 1290 } 1291 1292 #ifdef debugPutOnMassShell 1293 G4cout << "success " << success << G4en 1294 #endif 1295 1296 } while ( ( ! success ) && 1297 ++loopCounter < maxNumberOfLoops 1298 if ( loopCounter >= maxNumberOfLoops ) { 1299 success = false; 1300 } 1301 1302 return success; 1303 } 1304 1305 1306 //=========================================== 1307 1308 G4bool G4QGSParticipants:: 1309 CheckKinematics( const G4double sValue, 1310 const G4double sqrtS, 1311 const G4double projectileMas 1312 const G4double targetMass2, 1313 const G4double nucleusY, 1314 const G4bool isProjectileNuc 1315 const G4int numberOfInvolved 1316 G4Nucleon* involvedNucleons[ 1317 G4double& targetWminus, 1318 G4double& projectileWplus, 1319 G4bool& success ) { 1320 1321 // This method, which is called only by Put 1322 // kinematics is acceptable or not. 1323 // This method assumes that all the paramet 1324 // notice that the input boolean parameter 1325 // only in the case of nucleus or antinucle 1326 // The action of this method consists in co 1327 // and setting the parameter success to fal 1328 // be rejeted. 1329 1330 G4double decayMomentum2 = sqr( sValue ) + s 1331 - 2.0*sValue*proj 1332 - 2.0*projectileM 1333 targetWminus = ( sValue - projectileMass2 + 1334 / 2.0 / sqrtS; 1335 projectileWplus = sqrtS - targetMass2/targe 1336 G4double projectilePz = projectileWplus/2.0 1337 G4double projectileE = projectileWplus/2.0 1338 G4double projectileY(1.0e5); 1339 if (projectileE - projectilePz > 0.) { 1340 projectileY = 0.5 * G4Log( (proje 1341 (proje 1342 } 1343 G4double targetPz = -targetWminus/2.0 + tar 1344 G4double targetE = targetWminus/2.0 + tar 1345 G4double targetY = 0.5 * G4Log( (targetE + 1346 1347 #ifdef debugPutOnMassShell 1348 G4cout << "decayMomentum2 " << decayMomen 1349 << "\t targetWminus projectileWplu 1350 << "\t projectileY targetY " << pr 1351 #endif 1352 1353 for ( G4int i = 0; i < numberOfInvolvedNucl 1354 G4Nucleon* aNucleon = involvedNucleons[i] 1355 if ( ! aNucleon ) continue; 1356 G4LorentzVector tmp = aNucleon->Get4Momen 1357 G4double mt2 = sqr( tmp.x() ) + sqr( tmp. 1358 sqr( aNucleon->GetSplitabl 1359 G4double x = tmp.z(); 1360 G4double pz = -targetWminus*x/2.0 + mt2/( 1361 G4double e = targetWminus*x/2.0 + mt2/( 1362 if ( isProjectileNucleus ) { 1363 pz = projectileWplus*x/2.0 - mt2/(2.0*p 1364 e = projectileWplus*x/2.0 + mt2/(2.0*p 1365 } 1366 G4double nucleonY = 0.5 * G4Log( (e + pz) 1367 1368 #ifdef debugPutOnMassShell 1369 G4cout << "i nY pY nY-AY AY " << i << " 1370 #endif 1371 1372 if ( std::abs( nucleonY - nucleusY ) > 2 1373 ( isProjectileNucleus && targetY > 1374 ( ! isProjectileNucleus && project 1375 success = false; 1376 break; 1377 } 1378 } 1379 return true; 1380 } 1381 1382 1383 //=========================================== 1384 1385 G4bool G4QGSParticipants:: 1386 FinalizeKinematics( const G4double w, 1387 const G4bool isProjectile 1388 const G4LorentzRotation& 1389 const G4double residualMa 1390 const G4int residualMassN 1391 const G4int numberOfInvol 1392 G4Nucleon* involvedNucleo 1393 G4LorentzVector& residual4Momen 1394 1395 // This method, which is called only by Put 1396 // this method is called when we are sure t 1397 // acceptable. 1398 // This method assumes that all the paramet 1399 // notice that the input boolean parameter 1400 // only in the case of nucleus or antinucle 1401 // because the sign of pz (in the center-of 1402 // with respect to the case of a normal had 1403 // The action of this method consists in mo 1404 // (in the lab frame) and computing the res 1405 // frame). 1406 1407 G4ThreeVector residual3Momentum( 0.0, 0.0, 1408 1409 for ( G4int i = 0; i < numberOfInvolvedNucl 1410 G4Nucleon* aNucleon = involvedNucleons[i] 1411 if ( ! aNucleon ) continue; 1412 G4LorentzVector tmp = aNucleon->Get4Momen 1413 residual3Momentum -= tmp.vect(); 1414 G4double mt2 = sqr( tmp.x() ) + sqr( tmp. 1415 sqr( aNucleon->GetSplitabl 1416 G4double x = tmp.z(); 1417 G4double pz = -w * x / 2.0 + mt2 / ( 2. 1418 G4double e = w * x / 2.0 + mt2 / ( 2. 1419 // Reverse the sign of pz in the case of 1420 if ( isProjectileNucleus ) pz *= -1.0; 1421 tmp.setPz( pz ); 1422 tmp.setE( e ); 1423 tmp.transform( boostFromCmsToLab ); 1424 aNucleon->SetMomentum( tmp ); 1425 G4VSplitableHadron* splitableHadron = aNu 1426 splitableHadron->Set4Momentum( tmp ); 1427 #ifdef debugPutOnMassShell 1428 G4cout << "Target involved nucleon No, 1429 << i<<" "<<aNucleon->GetDefinitio 1430 #endif 1431 } 1432 1433 G4double residualMt2 = sqr( residualMass ) 1434 + sqr( residual3Moment 1435 1436 #ifdef debugPutOnMassShell 1437 G4cout <<G4endl<< "w residual3Momentum.z( 1438 #endif 1439 1440 G4double residualPz = 0.0; 1441 G4double residualE = 0.0; 1442 if ( residualMassNumber != 0 ) { 1443 residualPz = -w * residual3Momentum.z() / 1444 residualMt2 / ( 2.0 * w * r 1445 residualE = w * residual3Momentum.z() / 1446 residualMt2 / ( 2.0 * w * r 1447 // Reverse the sign of residualPz in the 1448 if ( isProjectileNucleus ) residualPz *= 1449 } 1450 1451 residual4Momentum.setPx( residual3Momentum. 1452 residual4Momentum.setPy( residual3Momentum. 1453 residual4Momentum.setPz( residualPz ); 1454 residual4Momentum.setE( residualE ); 1455 1456 return true; 1457 } 1458 1459 //=========================================== 1460 void G4QGSParticipants::PerformDiffractiveCol 1461 { 1462 #ifdef debugQGSParticipants 1463 G4cout<<G4endl<<"PerformDiffractiveCollisi 1464 <<"theInteractions.size() "<<theInte 1465 #endif 1466 1467 unsigned int i; 1468 for (i = 0; i < theInteractions.size(); i++ 1469 { 1470 G4InteractionContent* anIniteraction = th 1471 #ifdef debugQGSParticipants 1472 G4cout<<"Interaction # and its status " 1473 <<i<<" "<<theInteractions[i]->Get 1474 #endif 1475 1476 G4int InterStatus = theInteractions[i]->G 1477 if ( (InterStatus == PrD) || (InterStatus 1478 { // Selection of diffractive interactio 1479 #ifdef debugQGSParticipants 1480 G4cout<<"Simulation of diffractive in 1481 #endif 1482 1483 G4VSplitableHadron* aTarget = anInitera 1484 1485 #ifdef debugQGSParticipants 1486 G4cout<<"The proj. before inter " 1487 <<theProjectileSplitable->Get4M 1488 <<theProjectileSplitable->Get4M 1489 G4cout<<"The targ. before inter " <<a 1490 <<aTarget->Get4Momentum().mag() 1491 #endif 1492 1493 if ( InterStatus == PrD ) 1494 theSingleDiffExcitation.ExcitePartici 1495 1496 if ( InterStatus == TrD ) 1497 theSingleDiffExcitation.ExcitePartici 1498 1499 if ( InterStatus == DD ) 1500 theDiffExcitaton.ExciteParticipants(t 1501 1502 #ifdef debugQGSParticipants 1503 G4cout<<"The proj. after inter " <<t 1504 <<theProjectileSplitable->Get4M 1505 G4cout<<"The targ. after inter " <<aTarget 1506 <<aTarget->Get4Momentum().mag() 1507 #endif 1508 } 1509 1510 if ( InterStatus == Qexc ) 1511 { // Quark exchange process 1512 #ifdef debugQGSParticipants 1513 G4cout<<"Simulation of interaction wi 1514 #endif 1515 G4VSplitableHadron* aTarget = anInitera 1516 1517 #ifdef debugQGSParticipants 1518 G4cout<<"The proj. before inter " <<t 1519 <<theProjectileSplitable->Get4M 1520 G4cout<<"The targ. before inter "<<aT 1521 <<aTarget->Get4Momentum().mag() 1522 #endif 1523 1524 theQuarkExchange.ExciteParticipants(the 1525 1526 #ifdef debugQGSParticipants 1527 G4cout<<"The proj. after inter " <<t 1528 <<theProjectileSplitable->Get4M 1529 G4cout<<"The targ. after inter " <<aTarget 1530 <<aTarget->Get4Momentum().mag() 1531 #endif 1532 } 1533 } 1534 } 1535 1536 //=========================================== 1537 G4bool G4QGSParticipants::DeterminePartonMome 1538 { 1539 if ( ! theProjectileSplitable ) return fals 1540 1541 const G4double aHugeValue = 1.0e+10; 1542 1543 #ifdef debugQGSParticipants 1544 G4cout<<G4endl<<"DeterminePartonMomenta() 1545 G4cout<<"theProjectile status (0 -nondiff 1546 #endif 1547 1548 if (theProjectileSplitable->GetStatus() != 1549 1550 G4LorentzVector Projectile4Momentum = theP 1551 G4LorentzVector Psum = Projectile4Momentum; 1552 1553 G4double VqM_pr(0.), VaqM_pr(0.), VqM_tr(35 1554 if (std::abs(theProjectile.GetDefinition()- 1555 1556 #ifdef debugQGSParticipants 1557 G4cout<<"Projectile 4 momentum "<<Psum<<G 1558 <<"Target nucleon momenta at start" 1559 G4int NuclNo=0; 1560 #endif 1561 1562 std::vector<G4VSplitableHadron*>::iterator 1563 1564 for (i = theTargets.begin(); i != theTarget 1565 { 1566 Psum += (*i)->Get4Momentum(); 1567 #ifdef debugQGSParticipants 1568 G4cout<<"Nusleus nucleon # and its 4Mom 1569 NuclNo++; 1570 #endif 1571 } 1572 1573 G4LorentzRotation toCms( -1*Psum.boostVecto 1574 1575 G4LorentzVector Ptmp = toCms*Projectile4Mom 1576 1577 toCms.rotateZ( -1*Ptmp.phi() ); 1578 toCms.rotateY( -1*Ptmp.theta() ); 1579 G4LorentzRotation toLab(toCms.inverse()); 1580 Projectile4Momentum.transform( toCms ); 1581 // Ptarget.transform( toCms ); 1582 1583 #ifdef debugQGSParticipants 1584 G4cout<<G4endl<<"In CMS---------------"<< 1585 G4cout<<"Projectile 4 Mom "<<Projectile4M 1586 NuclNo=0; 1587 #endif 1588 1589 G4LorentzVector Target4Momentum(0.,0.,0.,0. 1590 for(i = theTargets.begin(); i != theTargets 1591 { 1592 G4LorentzVector tmp= (*i)->Get4Momentum() 1593 (*i)->Set4Momentum( tmp ); 1594 #ifdef debugQGSParticipants 1595 G4cout<<"Target nucleon # and 4Mom "<<" 1596 NuclNo++; 1597 #endif 1598 Target4Momentum += tmp; 1599 } 1600 1601 G4double S = Psum.mag2(); 1602 G4double SqrtS = std::sqrt(S); 1603 1604 #ifdef debugQGSParticipants 1605 G4cout<<"Sum of target nucleons 4 momentu 1606 G4cout<<"Target nucleons mom: px, py, z_1 1607 NuclNo=0; 1608 #endif 1609 1610 //G4double PplusProjectile = Projectile4Mom 1611 G4double PminusTarget = Target4Momentum. 1612 1613 for(i = theTargets.begin(); i != theTargets 1614 { 1615 G4LorentzVector tmp = (*i)->Get4Momentum( 1616 1617 //AR-19Jan2017 : the following line is ca 1618 // is built in optimized mo 1619 // To fix it, I get the mas 1620 // mass from the Lorentz ve 1621 // square root. If the mass 1622 // exception is thrown, and 1623 //G4double Mass = tmp.mag(); 1624 G4double Mass2 = tmp.mag2(); 1625 G4double Mass = 0.0; 1626 if ( Mass2 < 0.0 ) { 1627 G4ExceptionDescription ed; 1628 ed << "Projectile " << theProjectile.Ge 1629 << "Â 4-momentum " << Psum << G4end 1630 ed << "LorentzVector tmp " << tmp << " 1631 G4Exception( "G4QGSParticipants::Determ 1632 "HAD_QGSPARTICIPANTS_001", 1633 } else { 1634 Mass = std::sqrt( Mass2 ); 1635 } 1636 1637 tmp.setPz(tmp.minus()/PminusTarget); tm 1638 (*i)->Set4Momentum(tmp); 1639 #ifdef debugQGSParticipants 1640 G4cout<<"Target nucleons # and mom: "<< 1641 NuclNo++; 1642 #endif 1643 } 1644 1645 //+++++++++++++++++++++++++++++++++++++++++ 1646 1647 G4double SigPt = sigmaPt; 1648 G4Parton* aParton(0); 1649 G4ThreeVector aPtVector(0.,0.,0.); 1650 G4LorentzVector tmp(0.,0.,0.,0.); 1651 1652 G4double Mt(0.); 1653 G4double ProjSumMt(0.), ProjSumMt2perX(0.); 1654 G4double TargSumMt(0.), TargSumMt2perX(0.); 1655 1656 1657 G4double aBeta = beta; // Member of the c 1658 1659 const G4ParticleDefinition* theProjectileDe 1660 if (theProjectileDefinition == G4PionMinus: 1661 if (theProjectileDefinition == G4Gamma::Gam 1662 if (theProjectileDefinition == G4PionPlus:: 1663 if (theProjectileDefinition == G4PionZero:: 1664 if (theProjectileDefinition == G4KaonPlus:: 1665 if (theProjectileDefinition == G4KaonMinus: 1666 1667 G4double Xmin = 0.; 1668 1669 G4bool Success = true; G4int attempt = 0; 1670 const G4int maxNumberOfAttempts = 1000; 1671 do 1672 { 1673 attempt++; if( attempt == 100*(attempt/ 1674 1675 ProjSumMt=0.; ProjSumMt2perX=0.; 1676 TargSumMt=0.; TargSumMt2perX=0.; 1677 1678 Success = true; 1679 G4int nSeaPair = theProjectileSplitable-> 1680 #ifdef debugQGSParticipants 1681 G4cout<<"attempt ---------------------- 1682 G4cout<<"nSeaPair of proj "<<nSeaPair<< 1683 #endif 1684 1685 G4double SumPx = 0.; 1686 G4double SumPy = 0.; 1687 G4double SumZ = 0.; 1688 G4int NumberOfUnsampledSeaQuarks = 2*nSea 1689 1690 G4double Qmass=0.; 1691 for (G4int aSeaPair = 0; aSeaPair < nSeaP 1692 { 1693 aParton = theProjectileSplitable->GetNe 1694 #ifdef debugQGSParticipants 1695 G4cout<<"Sea quarks: "<<aSeaPair<<" " 1696 #endif 1697 aPtVector = GaussianPt(SigPt, aHugeValu 1698 tmp.setPx(aPtVector.x()); tmp.setPy(aPt 1699 SumPx += aPtVector.x(); SumPy += aPtV 1700 Mt = std::sqrt(aPtVector.mag2()+sqr(Qma 1701 ProjSumMt += Mt; 1702 1703 // Sampling of Z fraction 1704 tmp.setPz(SampleX(Xmin, NumberOfUnsampl 1705 SumZ += tmp.z(); 1706 1707 NumberOfUnsampledSeaQuarks--; 1708 ProjSumMt2perX +=sqr(Mt)/tmp.pz(); 1709 tmp.setE(sqr(Mt)); 1710 aParton->Set4Momentum(tmp); 1711 1712 aParton = theProjectileSplitable->GetNe 1713 #ifdef debugQGSParticipants 1714 G4cout<<" "<<aParton->GetDefinition() 1715 G4cout<<" "<<tmp<<" "<<S 1716 #endif 1717 aPtVector = GaussianPt(SigPt, aHugeValu 1718 tmp.setPx(aPtVector.x()); tmp.setPy(aPt 1719 SumPx += aPtVector.x(); SumPy += aPtV 1720 Mt = std::sqrt(aPtVector.mag2()+sqr(Qma 1721 ProjSumMt += Mt; 1722 1723 // Sampling of Z fraction 1724 tmp.setPz(SampleX(Xmin, NumberOfUnsampl 1725 SumZ += tmp.z(); 1726 1727 NumberOfUnsampledSeaQuarks--; 1728 ProjSumMt2perX +=sqr(Mt)/tmp.pz(); 1729 tmp.setE(sqr(Mt)); 1730 aParton->Set4Momentum(tmp); 1731 #ifdef debugQGSParticipants 1732 G4cout<<" "<<tmp<<" "<<S 1733 #endif 1734 } 1735 1736 // For valence quark 1737 aParton = theProjectileSplitable->GetNext 1738 #ifdef debugQGSParticipants 1739 G4cout<<"Val quark of Pr"<<" "<<aParton 1740 #endif 1741 aPtVector = GaussianPt(SigPt, aHugeValue) 1742 tmp.setPx(aPtVector.x()); tmp.setPy(aPtVe 1743 SumPx += aPtVector.x(); SumPy += aPtVec 1744 Mt = std::sqrt(aPtVector.mag2()+sqr(VqM_p 1745 ProjSumMt += Mt; 1746 1747 // Sampling of Z fraction 1748 tmp.setPz(SampleX(Xmin, NumberOfUnsampled 1749 SumZ += tmp.z(); 1750 1751 ProjSumMt2perX +=sqr(Mt)/tmp.pz(); 1752 tmp.setE(sqr(Mt)); 1753 aParton->Set4Momentum(tmp); 1754 1755 // For valence di-quark 1756 aParton = theProjectileSplitable->GetNext 1757 #ifdef debugQGSParticipants 1758 G4cout<<" "<<aParton->GetDefinition()-> 1759 G4cout<<" "<<tmp<<" "<<Sum 1760 #endif 1761 tmp.setPx(-SumPx); tmp.setPy(-SumPy); 1762 //Uzhi 2019 Mt = std::sqrt(aPtVector.mag 1763 Mt = std::sqrt( sqr(SumPx) + sqr(SumPy) + 1764 ProjSumMt += Mt; 1765 tmp.setPz(1.-SumZ); 1766 1767 ProjSumMt2perX +=sqr(Mt)/tmp.pz(); // QQ 1768 tmp.setE(sqr(Mt)); 1769 aParton->Set4Momentum(tmp); 1770 #ifdef debugQGSParticipants 1771 G4cout<<" "<<tmp<<" "<<Sum 1772 NuclNo=0; 1773 #endif 1774 1775 // End of work with the projectile 1776 1777 // Work with target nucleons 1778 1779 for(i = theTargets.begin(); i != theTarge 1780 { 1781 nSeaPair = (*i)->GetSoftCollisionCount( 1782 #ifdef debugQGSParticipants 1783 G4cout<<"nSeaPair of target N "<<nSea 1784 <<"Target nucleon 4Mom "<<(*i)- 1785 #endif 1786 1787 SumPx = (*i)->Get4Momentum().px() * (-1 1788 SumPy = (*i)->Get4Momentum().py() * (-1 1789 SumZ = 0.; 1790 1791 NumberOfUnsampledSeaQuarks = 2*nSeaPair 1792 1793 Qmass=0; 1794 for (G4int aSeaPair = 0; aSeaPair < nSe 1795 { 1796 aParton = (*i)->GetNextParton(); // 1797 #ifdef debugQGSParticipants 1798 G4cout<<"Sea quarks: "<<aSeaPair<<" 1799 #endif 1800 aPtVector = GaussianPt(SigPt, aHugeVa 1801 tmp.setPx(aPtVector.x()); tmp.setPy(a 1802 SumPx += aPtVector.x(); SumPy += aP 1803 Mt=std::sqrt(aPtVector.mag2()+sqr(Qma 1804 TargSumMt += Mt; 1805 1806 // Sampling of Z fraction 1807 tmp.setPz(SampleX(Xmin, NumberOfUnsam 1808 SumZ += tmp.z(); 1809 tmp.setPz((*i)->Get4Momentum().pz()*t 1810 NumberOfUnsampledSeaQuarks--; 1811 TargSumMt2perX +=sqr(Mt)/tmp.pz(); 1812 tmp.setE(sqr(Mt)); 1813 aParton->Set4Momentum(tmp); 1814 1815 aParton = (*i)->GetNextAntiParton(); 1816 #ifdef debugQGSParticipants 1817 G4cout<<" "<<aParton->GetDefinition 1818 G4cout<<" "<<tmp<<" "< 1819 #endif 1820 aPtVector = GaussianPt(SigPt, aHugeVa 1821 tmp.setPx(aPtVector.x()); tmp.setPy(a 1822 SumPx += aPtVector.x(); SumPy += aP 1823 Mt=std::sqrt(aPtVector.mag2()+sqr(Qma 1824 TargSumMt += Mt; 1825 1826 // Sampling of Z fraction 1827 tmp.setPz(SampleX(Xmin, NumberOfUnsam 1828 SumZ += tmp.z(); 1829 tmp.setPz((*i)->Get4Momentum().pz()*t 1830 NumberOfUnsampledSeaQuarks--; 1831 TargSumMt2perX +=sqr(Mt)/tmp.pz(); 1832 tmp.setE(sqr(Mt)); 1833 aParton->Set4Momentum(tmp); 1834 #ifdef debugQGSParticipants 1835 G4cout<<" "<<tmp<<" "< 1836 #endif 1837 } 1838 1839 // Valence quark 1840 aParton = (*i)->GetNextParton(); // f 1841 #ifdef debugQGSParticipants 1842 G4cout<<"Val quark of Tr"<<" "<<aPart 1843 #endif 1844 aPtVector = GaussianPt(SigPt, aHugeValu 1845 tmp.setPx(aPtVector.x()); tmp.setPy(aPt 1846 SumPx += aPtVector.x(); SumPy += aPtV 1847 Mt=std::sqrt(aPtVector.mag2()+sqr(VqM_t 1848 TargSumMt += Mt; 1849 1850 // Sampling of Z fraction 1851 tmp.setPz(SampleX(Xmin, NumberOfUnsampl 1852 SumZ += tmp.z(); 1853 tmp.setPz((*i)->Get4Momentum().pz()*tmp 1854 TargSumMt2perX +=sqr(Mt)/tmp.pz(); 1855 tmp.setE(sqr(Mt)); 1856 aParton->Set4Momentum(tmp); 1857 1858 // Valence di-quark 1859 aParton = (*i)->GetNextAntiParton(); 1860 #ifdef debugQGSParticipants 1861 G4cout<<" "<<aParton->GetDefinition() 1862 G4cout<<" "<<tmp<<" "<<S 1863 #endif 1864 tmp.setPx(-SumPx); tmp 1865 //Uzhi 2019 Mt=std::sqrt(aPtVector.mag 1866 Mt=std::sqrt( sqr(SumPx) + sqr(SumPy) + 1867 TargSumMt += Mt; 1868 1869 tmp.setPz((*i)->Get4Momentum().pz()*(1. 1870 TargSumMt2perX +=sqr(Mt)/tmp.pz(); 1871 tmp.setE(sqr(Mt)); 1872 aParton->Set4Momentum(tmp); 1873 #ifdef debugQGSParticipants 1874 G4cout<<" "<<tmp<<" "<<1 1875 #endif 1876 1877 } // End of for(i = theTargets.begin(); 1878 1879 if( ProjSumMt + TargSumMt > Sqr 1880 Success = false; continue;} 1881 if( std::sqrt(ProjSumMt2perX) + std::sqrt 1882 Success = false; continue;} 1883 1884 } while( (!Success) && 1885 attempt < maxNumberOfAttempts ); 1886 1887 if ( attempt >= maxNumberOfAttempts ) { 1888 return false; 1889 } 1890 1891 //+++++++++++++++++++++++++++++++++++++++++ 1892 1893 G4double DecayMomentum2 = sqr(S) + sqr(Proj 1894 - 2.0*S*ProjSumMt2perX - 2.0*S 1895 1896 G4double targetWminus=( S - ProjSumMt2perX 1897 G4double projectileWplus = SqrtS - TargSumM 1898 1899 G4LorentzVector Tmp(0.,0.,0.,0.); 1900 G4double z(0.); 1901 1902 G4int nSeaPair = theProjectileSplitable->Ge 1903 #ifdef debugQGSParticipants 1904 G4cout<<"Backward transformation ======== 1905 G4cout<<"nSeaPair of proj "<<nSeaPair<<G4 1906 #endif 1907 1908 for (G4int aSeaPair = 0; aSeaPair < nSeaPai 1909 { 1910 aParton = theProjectileSplitable->GetNext 1911 #ifdef debugQGSParticipants 1912 G4cout<<"Sea quarks: "<<aSeaPair<<" "<< 1913 #endif 1914 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 1915 1916 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e() 1917 Tmp.setE( projectileWplus*z/2.0 + Tmp.e() 1918 Tmp.transform( toLab ); 1919 1920 aParton->Set4Momentum(Tmp); 1921 1922 aParton = theProjectileSplitable->GetNext 1923 #ifdef debugQGSParticipants 1924 G4cout<<" "<<aParton->GetDefinition()-> 1925 G4cout<<" "<<Tmp<<" "<<Tmp 1926 #endif 1927 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 1928 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e() 1929 Tmp.setE( projectileWplus*z/2.0 + Tmp.e() 1930 Tmp.transform( toLab ); 1931 1932 aParton->Set4Momentum(Tmp); 1933 #ifdef debugQGSParticipants 1934 G4cout<<" "<<Tmp<<" "<<Tmp 1935 #endif 1936 } 1937 1938 // For valence quark 1939 aParton = theProjectileSplitable->GetNextPa 1940 #ifdef debugQGSParticipants 1941 G4cout<<"Val quark of Pr"<<" "<<aParton-> 1942 #endif 1943 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 1944 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/( 1945 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/( 1946 Tmp.transform( toLab ); 1947 1948 aParton->Set4Momentum(Tmp); 1949 1950 // For valence di-quark 1951 aParton = theProjectileSplitable->GetNextAn 1952 #ifdef debugQGSParticipants 1953 G4cout<<" "<<aParton->GetDefinition()->Ge 1954 G4cout<<" "<<Tmp<<" "<<Tmp.m 1955 #endif 1956 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 1957 Tmp.setPz(projectileWplus*z/2.0 - Tmp.e()/( 1958 Tmp.setE( projectileWplus*z/2.0 + Tmp.e()/( 1959 Tmp.transform( toLab ); 1960 1961 aParton->Set4Momentum(Tmp); 1962 1963 #ifdef debugQGSParticipants 1964 G4cout<<" "<<Tmp<<" "<<Tmp.m 1965 NuclNo=0; 1966 #endif 1967 1968 // End of work with the projectile 1969 1970 // Work with target nucleons 1971 for(i = theTargets.begin(); i != theTargets 1972 { 1973 nSeaPair = (*i)->GetSoftCollisionCount()- 1974 #ifdef debugQGSParticipants 1975 G4cout<<"nSeaPair of target and N# "<<n 1976 NuclNo++; 1977 #endif 1978 for (G4int aSeaPair = 0; aSeaPair < nSeaP 1979 { 1980 aParton = (*i)->GetNextParton(); // f 1981 #ifdef debugQGSParticipants 1982 G4cout<<"Sea quarks: "<<aSeaPair<<" " 1983 #endif 1984 Tmp =aParton->Get4Momentum(); z=Tmp.z() 1985 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e() 1986 Tmp.setE( targetWminus*z/2.0 + Tmp.e() 1987 Tmp.transform( toLab ); 1988 1989 aParton->Set4Momentum(Tmp); 1990 1991 aParton = (*i)->GetNextAntiParton(); 1992 #ifdef debugQGSParticipants 1993 G4cout<<" "<<aParton->GetDefinition() 1994 G4cout<<" "<<Tmp<<" "<<T 1995 #endif 1996 Tmp =aParton->Get4Momentum(); z=Tmp.z() 1997 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e() 1998 Tmp.setE( targetWminus*z/2.0 + Tmp.e() 1999 Tmp.transform( toLab ); 2000 2001 aParton->Set4Momentum(Tmp); 2002 #ifdef debugQGSParticipants 2003 G4cout<<" "<<Tmp<<" "<<T 2004 #endif 2005 } 2006 2007 // Valence quark 2008 2009 aParton = (*i)->GetNextParton(); // for 2010 #ifdef debugQGSParticipants 2011 G4cout<<"Val quark of Tr"<<" "<<aParton 2012 #endif 2013 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 2014 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/( 2015 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/( 2016 Tmp.transform( toLab ); 2017 2018 aParton->Set4Momentum(Tmp); 2019 2020 // Valence di-quark 2021 aParton = (*i)->GetNextAntiParton(); // 2022 #ifdef debugQGSParticipants 2023 G4cout<<" "<<aParton->GetDefinition()-> 2024 G4cout<<" "<<Tmp<<" "<<Tmp 2025 #endif 2026 Tmp =aParton->Get4Momentum(); z=Tmp.z(); 2027 Tmp.setPz(-targetWminus*z/2.0 + Tmp.e()/( 2028 Tmp.setE( targetWminus*z/2.0 + Tmp.e()/( 2029 Tmp.transform( toLab ); 2030 2031 aParton->Set4Momentum(Tmp); 2032 #ifdef debugQGSParticipants 2033 G4cout<<" "<<Tmp<<" "<<Tmp 2034 NuclNo++; 2035 #endif 2036 } // End of for(i = theTargets.begin(); i 2037 2038 return true; 2039 } 2040 2041 //=========================================== 2042 G4double G4QGSParticipants:: 2043 SampleX(G4double, G4int nSea, G4int, G4double 2044 { 2045 G4double Oalfa = 1./(alpha + 1.); 2046 G4double Obeta = 1./(aBeta + (alpha + 1.)*n 2047 2048 G4double Ksi1, Ksi2, r1, r2, r12; 2049 const G4int maxNumberOfLoops = 1000; 2050 G4int loopCounter = 0; 2051 do 2052 { 2053 Ksi1 = G4UniformRand(); r1 = G4Pow::GetIn 2054 Ksi2 = G4UniformRand(); r2 = G4Pow::GetIn 2055 r12=r1+r2; 2056 } while( ( r12 > 1.) && 2057 ++loopCounter < maxNumberOfLoops ) 2058 if ( loopCounter >= maxNumberOfLoops ) { 2059 return 0.5; // Just an acceptable value, 2060 } 2061 2062 G4double result = r1/r12; 2063 return result; 2064 } 2065 2066 //=========================================== 2067 void G4QGSParticipants::CreateStrings() 2068 { 2069 2070 #ifdef debugQGSParticipants 2071 G4cout<<"CreateStrings() ................ 2072 #endif 2073 2074 if ( ! theProjectileSplitable ) { 2075 #ifdef debugQGSParticipants 2076 G4cout<<"BAD situation: theProjectileSp 2077 #endif 2078 return; 2079 } 2080 2081 #ifdef debugQGSParticipants 2082 G4cout<<"theProjectileSplitable->GetStatu 2083 G4LorentzVector str4Mom; 2084 #endif 2085 2086 if( theProjectileSplitable->GetStatus() == 2087 { 2088 G4ThreeVector Position = theProjectileSpl 2089 2090 G4PartonPair * aPair = new G4PartonPair(t 2091 t 2092 G4PartonPair::DIFFRACTIVE 2093 #ifdef debugQGSParticipants 2094 G4cout << "Pr. Diffr. String: Qs 4mom X 2095 G4cout << " " << aPair->Ge 2096 << aPair->GetParton1()->Get4Momentum 2097 << aPair->GetParton1()->GetX() 2098 G4cout << " " << aPair->Ge 2099 << aPair->GetParton2()->Get4Momentum 2100 << aPair->GetParton2()->GetX() 2101 str4Mom += aPair->GetParton1()->Get4Mom 2102 str4Mom += aPair->GetParton2()->Get4Mom 2103 #endif 2104 2105 thePartonPairs.push_back(aPair); 2106 } 2107 2108 G4int N_EnvTarg = NumberOfInvolvedNucleonsO 2109 2110 for ( G4int i = 0; i < N_EnvTarg; i++ ) { 2111 G4Nucleon* aTargetNucleon = TheInvolvedNu 2112 2113 G4VSplitableHadron* HitTargetNucleon = aT 2114 2115 #ifdef debugQGSParticipants 2116 G4cout<<"Involved Nucleon # and its sta 2117 #endif 2118 if( HitTargetNucleon->GetStatus() >= 1) 2119 { 2120 G4ThreeVector Position = HitTargetN 2121 2122 G4PartonPair * aPair = new G4PartonPair 2123 2124 G4PartonPair::DIFFRACTI 2125 #ifdef debugQGSParticipants 2126 G4cout << "Tr. Diffr. String: Qs 4mom 2127 G4cout << "Diffr. String " << aPair->GetPar 2128 << aPair->GetParton1()->Get4Moment 2129 << aPair->GetParton1()->GetX() 2130 G4cout << " " << aPair->GetPar 2131 << aPair->GetParton2()->Get4Moment 2132 << aPair->GetParton2()->GetX() 2133 2134 str4Mom += aPair->GetParton1()->Get4Momentu 2135 str4Mom += aPair->GetParton2()->Get4Momentu 2136 #endif 2137 2138 thePartonPairs.push_back(aPair); 2139 } 2140 } 2141 2142 //----------------------------------------- 2143 #ifdef debugQGSParticipants 2144 G4int IntNo=0; 2145 G4cout<<"Strings created in soft interact 2146 #endif 2147 std::vector<G4InteractionContent*>::iterato 2148 i = theInteractions.begin(); 2149 while ( i != theInteractions.end() ) /* Lo 2150 { 2151 G4InteractionContent* anIniteraction = *i 2152 G4PartonPair * aPair = NULL; 2153 2154 #ifdef debugQGSParticipants 2155 G4cout<<"An interaction # and soft coll 2156 <<anIniteraction->GetNumberOfSoft 2157 IntNo++; 2158 #endif 2159 if (anIniteraction->GetNumberOfSoftCollis 2160 { 2161 G4VSplitableHadron* pProjectile = anIni 2162 G4VSplitableHadron* pTarget = anIni 2163 2164 for (G4int j = 0; j < anIniteraction->G 2165 { 2166 aPair = new G4PartonPair(pTarget->Get 2167 G4PartonPair::SOFT, G4PartonPair::TA 2168 #ifdef debugQGSParticipants 2169 G4cout << "SoftPair " << aPair->GetParton 2170 << aPair->GetParton1()->Get4Momentum() < 2171 << aPair->GetParton1()->Get4Momentum().m 2172 G4cout << " " << aPair->GetParton 2173 << aPair->GetParton2()->Get4Momentum() < 2174 <<aPair->GetParton2()->Get4Momentum().m 2175 str4Mom += aPair->GetParton1()->Get4Momen 2176 str4Mom += aPair->GetParton2()->Get4Momen 2177 #endif 2178 2179 thePartonPairs.push_back(aPair); 2180 2181 aPair = new G4PartonPair(pProjectile->GetNe 2182 G4PartonPair::SOFT, G4PartonPair::PR 2183 #ifdef debugQGSParticipants 2184 G4cout << "SoftPair " << aPair->GetParton 2185 << aPair->GetParton1()->Get4Momentum() < 2186 << aPair->GetParton1()->Get4Moment 2187 G4cout << " " << aPair->GetParton 2188 << aPair->GetParton2()->Get4Momentum() < 2189 << aPair->GetParton2()->Get4Momentum().m 2190 #endif 2191 #ifdef debugQGSParticipants 2192 str4Mom += aPair->GetParton1()->Get4Momen 2193 str4Mom += aPair->GetParton2()->Get4Momen 2194 #endif 2195 2196 thePartonPairs.push_back(aPair); 2197 } 2198 2199 delete *i; 2200 i=theInteractions.erase(i); // i now 2201 } 2202 else 2203 { 2204 i++; 2205 } 2206 } // End of while ( i != theInteractions.e 2207 #ifdef debugQGSParticipants 2208 G4cout << "Sum of strings 4 momenta " << 2209 #endif 2210 } 2211 2212 //=========================================== 2213 2214 void G4QGSParticipants::GetResiduals() { 2215 // This method is needed for the correct ap 2216 2217 #ifdef debugQGSParticipants 2218 G4cout << "GetResiduals(): GetProjectileN 2219 #endif 2220 2221 #ifdef debugQGSParticipants 2222 G4cout << "NumberOfInvolvedNucleonsOfTarg 2223 #endif 2224 2225 G4double DeltaExcitationE = TargetResidualE 2226 G4LorentzVector DeltaPResidualNucleus = Tar 2227 G4d 2228 2229 for ( G4int i = 0; i < NumberOfInvolvedNucl 2230 G4Nucleon* aNucleon = TheInvolvedNucleons 2231 2232 #ifdef debugQGSParticipants 2233 G4VSplitableHadron* targetSplitable = a 2234 G4cout << i << " Hit? " << aNucleon->Ar 2235 if ( targetSplitable ) G4cout << i << " 2236 #endif 2237 2238 G4LorentzVector tmp = -DeltaPResidualNucl 2239 aNucleon->SetMomentum( tmp ); 2240 aNucleon->SetBindingEnergy( DeltaExcitati 2241 } 2242 2243 //------------------------------------- 2244 if( TargetResidualMassNumber != 0 ) 2245 { 2246 G4ThreeVector bstToCM =TargetResidual4Mom 2247 2248 G4V3DNucleus* theTargetNucleus = GetTarge 2249 G4LorentzVector residualMomentum(0.,0.,0. 2250 G4Nucleon* aNucleon = 0; 2251 theTargetNucleus->StartLoop(); 2252 while ( ( aNucleon = theTargetNucleus->Ge 2253 if ( !aNucleon->AreYouHit() ) { 2254 G4LorentzVector tmp=aNucleon->Get4Mom 2255 aNucleon->SetMomentum(tmp); 2256 residualMomentum +=tmp; 2257 } 2258 } 2259 2260 residualMomentum/=TargetResidualMassNumbe 2261 2262 G4double Mass = TargetResidual4Momentum.m 2263 G4double SumMasses=0.; 2264 2265 aNucleon = 0; 2266 theTargetNucleus->StartLoop(); 2267 while ( ( aNucleon = theTargetNucleus->Ge 2268 if ( !aNucleon->AreYouHit() ) { 2269 G4LorentzVector tmp=aNucleon->Get4Mom 2270 G4double E=std::sqrt(tmp.vect().mag2( 2271 sqr(aNucleon->Ge 2272 tmp.setE(E); aNucleon->SetMomentum(t 2273 SumMasses+=E; 2274 } 2275 } 2276 2277 G4double Chigh=Mass/SumMasses; G4double C 2278 const G4int maxNumberOfLoops = 1000; 2279 G4int loopCounter = 0; 2280 do 2281 { 2282 C=(Chigh+Clow)/2.; 2283 2284 SumMasses=0.; 2285 aNucleon = 0; 2286 theTargetNucleus->StartLoop(); 2287 while ( ( aNucleon = theTargetNucleus-> 2288 if ( !aNucleon->AreYouHit() ) { 2289 G4LorentzVector tmp=aNucleon->Get4M 2290 G4double E=std::sqrt(tmp.vect().mag 2291 sqr(aNucleon-> 2292 SumMasses+=E; 2293 } 2294 } 2295 2296 if(SumMasses > Mass) {Chigh=C;} 2297 else {Clow =C;} 2298 2299 } while( (Chigh-Clow > 0.01) && 2300 ++loopCounter < maxNumberOfLoops 2301 if ( loopCounter >= maxNumberOfLoops ) { 2302 #ifdef debugQGSParticipants 2303 G4cout <<"BAD situation: forced loop 2304 #endif 2305 // Perhaps there is something to set he 2306 } else { 2307 aNucleon = 0; 2308 theTargetNucleus->StartLoop(); 2309 while ( ( aNucleon = theTargetNucleus-> 2310 if ( !aNucleon->AreYouHit() ) { 2311 G4LorentzVector tmp=aNucleon->Get4M 2312 G4double E=std::sqrt(tmp.vect().mag 2313 sqr(aNucleon-> 2314 tmp.setE(E); tmp.boost(-bstToCM); 2315 aNucleon->SetMomentum(tmp); 2316 } 2317 } 2318 } 2319 2320 } // End of if( TargetResidualMassNumber ! 2321 //------------------------------------- 2322 2323 #ifdef debugQGSParticipants 2324 G4cout << "End GetResiduals ------------- 2325 #endif 2326 2327 } 2328 2329 2330 //=========================================== 2331 void G4QGSParticipants::PerformSoftCollisions 2332 { 2333 std::vector<G4InteractionContent*>::iterato 2334 G4LorentzVector str4Mom; 2335 i = theInteractions.begin(); 2336 while ( i != theInteractions.end() ) /* Lo 2337 { 2338 G4InteractionContent* anIniteraction = *i 2339 G4PartonPair * aPair = NULL; 2340 if (anIniteraction->GetNumberOfSoftCollis 2341 { 2342 G4VSplitableHadron* pProjectile = anIni 2343 G4VSplitableHadron* pTarget = anIni 2344 for (G4int j = 0; j < anIniteraction->G 2345 { 2346 aPair = new G4PartonPair(pTarget->Get 2347 G4PartonPair::SOFT, G4PartonPair::TA 2348 #ifdef debugQGSParticipants 2349 G4cout << "SoftPair " << aPair->GetParton 2350 << aPair->GetParton1()->Get4Momentum() < 2351 << aPair->GetParton1()->GetX() << " " << 2352 G4cout << " " << aPair->GetParton 2353 << aPair->GetParton2()->Get4Momentum() < 2354 << aPair->GetParton2()->GetX() << " " << 2355 #endif 2356 #ifdef debugQGSParticipants 2357 str4Mom += aPair->GetParton1()->Get4Momen 2358 str4Mom += aPair->GetParton2()->Get4Momen 2359 #endif 2360 thePartonPairs.push_back(aPair); 2361 aPair = new G4PartonPair(pProjectile->GetNe 2362 G4PartonPair::SOFT, G4PartonPair::PR 2363 #ifdef debugQGSParticipants 2364 G4cout << "SoftPair " << aPair->GetParton 2365 << aPair->GetParton1()->Get4Momentum() < 2366 << aPair->GetParton1()->GetX() << " " << 2367 G4cout << " " << aPair->GetParton 2368 << aPair->GetParton2()->Get4Momentum() < 2369 << aPair->GetParton2()->GetX() << " " << 2370 #endif 2371 #ifdef debugQGSParticipants 2372 str4Mom += aPair->GetParton1()->Get4Momen 2373 str4Mom += aPair->GetParton2()->Get4Momen 2374 #endif 2375 thePartonPairs.push_back(aPair); 2376 } 2377 delete *i; 2378 i=theInteractions.erase(i); // i now 2379 } else { 2380 i++; 2381 } 2382 } 2383 #ifdef debugQGSParticipants 2384 G4cout << " string 4 mom " << str4Mom << 2385 #endif 2386 } 2387 2388 2389 //******************************************* 2390 G4VSplitableHadron* G4QGSParticipants::Select 2391 { 2392 // Check reaction threshold - goes to Chec 2393 2394 theProjectileSplitable = new G4QGSMSplitabl 2395 theProjectileSplitable->SetStatus(1); 2396 2397 G4LorentzVector aPrimaryMomentum(thePrimary 2398 G4LorentzVector aTargetNMomentum(0.,0.,0.,9 2399 2400 if ((!(aPrimaryMomentum.e()>-1)) && (!(aPri 2401 { 2402 throw G4HadronicException(__FILE__, __LIN 2403 "G4GammaParticipants::SelectInter 2404 } 2405 G4double S = (aPrimaryMomentum + aTargetNMo 2406 G4double ThresholdMass = thePrimary.GetMass 2407 ModelMode = SOFT; 2408 2409 #ifdef debugQGSParticipants 2410 G4cout << "Gamma projectile - SelectInter 2411 G4cout<<"Energy and Nucleus "<<thePrimary 2412 G4cout << "SqrtS ThresholdMass ModelMode 2413 #endif 2414 2415 if (sqr(ThresholdMass + ThresholdParameter) 2416 { 2417 ModelMode = DIFFRACTIVE; 2418 } 2419 if (sqr(ThresholdMass + QGSMThreshold) > S) 2420 { 2421 ModelMode = DIFFRACTIVE; 2422 } 2423 2424 // first find the collisions HPW 2425 std::for_each(theInteractions.begin(), theI 2426 theInteractions.clear(); 2427 2428 G4int theCurrent = G4int(theNucleus->GetMas 2429 G4int NucleonNo=0; 2430 2431 theNucleus->StartLoop(); 2432 G4Nucleon * pNucleon = 0; 2433 2434 while( (pNucleon = theNucleus->GetNextNucle 2435 { if(NucleonNo == theCurrent) break; Nucleo 2436 2437 if ( pNucleon ) { 2438 G4QGSMSplitableHadron* aTarget = new G4QG 2439 2440 pNucleon->Hit(aTarget); 2441 2442 if ( (0.06 > G4UniformRand() &&(ModelMode 2443 { 2444 G4InteractionContent * aInteraction = n 2445 theProjectileSplitable->SetStatus(1*the 2446 2447 aInteraction->SetTarget(aTarget); 2448 aInteraction->SetTargetNucleon(pNucleon 2449 aTarget->SetCollisionCount(0); 2450 aTarget->SetStatus(1); 2451 2452 aInteraction->SetNumberOfDiffractiveCol 2453 aInteraction->SetNumberOfSoftCollisions 2454 aInteraction->SetStatus(1); 2455 2456 theInteractions.push_back(aInteraction) 2457 } 2458 else 2459 { 2460 // nondiffractive soft interaction occu 2461 aTarget->IncrementCollisionCount(1); 2462 aTarget->SetStatus(0); 2463 theTargets.push_back(aTarget); 2464 2465 theProjectileSplitable->IncrementCollis 2466 theProjectileSplitable->SetStatus(0*the 2467 2468 G4InteractionContent * aInteraction = n 2469 aInteraction->SetTarget(aTarget); 2470 aInteraction->SetTargetNucleon(pNucleon 2471 aInteraction->SetNumberOfSoftCollisions 2472 aInteraction->SetStatus(3); 2473 theInteractions.push_back(aInteraction) 2474 } 2475 } 2476 return theProjectileSplitable; 2477 } 2478