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 // GEANT 4 class implementation file 30 // 31 // History: first implementation, Maxim K 32 // ------------------------------------------- 33 #include "G4QGSMFragmentation.hh" 34 #include "G4PhysicalConstants.hh" 35 #include "G4SystemOfUnits.hh" 36 #include "Randomize.hh" 37 #include "G4ios.hh" 38 #include "G4FragmentingString.hh" 39 #include "G4DiQuarks.hh" 40 #include "G4Quarks.hh" 41 #include "G4HadronicParameters.hh" 42 #include "G4Pow.hh" 43 44 //#define debug_QGSMfragmentation 45 46 // Class G4QGSMFragmentation 47 //******************************************** 48 49 G4QGSMFragmentation::G4QGSMFragmentation() 50 { 51 SigmaQT = 0.45 * GeV; 52 53 MassCut = 0.35*GeV; 54 55 SetStrangenessSuppression((1.0 - 0.16)/2.) 56 57 // Check if charmed and bottom hadrons are 58 // set the non-zero probabilities for c-cb 59 // else set them to 0.0. If these probabil 60 // hadrons can't/can be created during the 61 // (i.e. not heavy) projectile hadron nucl 62 if ( G4HadronicParameters::Instance()->Ena 63 SetProbCCbar(0.0002); // According to O 64 SetProbBBbar(5.0e-5); // According to O 65 } else { 66 SetProbCCbar(0.0); 67 SetProbBBbar(0.0); 68 } 69 70 SetDiquarkSuppression(0.32); 71 SetDiquarkBreakProbability(0.7); 72 73 SetMinMasses(); 74 75 arho = 0.5; // alpha_rho0 76 aphi = 0.0; // alpha_fi 77 aJPs =-2.2; // alpha_J/Psi 78 aUps =-8.0; // alpha_Y ??? O. Pisk 79 80 aksi =-1.0; 81 alft = 0.5; // 2 * alpha'_R *<Pt^2> 82 83 an = -0.5 ; 84 ala = -0.75; // an - arho/2 + aphi/2 85 alaC = an - arho/2.0 + aJPs/2.0; 86 alaB = an - arho/2.0 + aUps/2.0; 87 aXi = 0.0; // ?? 88 aXiC = 0.0; // ?? 89 aXiB = 0.0; // ?? 90 aXiCC = 0.0; // ?? 91 aXiCB = 0.0; // ?? 92 aXiBB = 0.0; // ?? 93 94 SetFFq2q(); 95 SetFFq2qq(); 96 SetFFqq2q(); 97 SetFFqq2qq(); 98 // d u s c b 99 G4int Index[5][5] = { { 0, 1, 2, 3, 4 } 100 { 1, 5, 6, 7, 8 } 101 { 2, 6, 9, 10, 11 } 102 { 3, 7, 10, 12, 13 } 103 { 4, 8, 11, 13, 14 } 104 for (G4int i = 0; i < 5; i++ ) { 105 for ( G4int j = 0; j < 5; j++ ) { 106 IndexDiQ[i][j] = Index[i][j]; 107 } 108 }; 109 } 110 111 G4QGSMFragmentation::~G4QGSMFragmentation() 112 {} 113 114 //-------------------------------------------- 115 116 G4KineticTrackVector* G4QGSMFragmentation::Fra 117 { 118 119 G4FragmentingString aString(theString); 120 SetMinimalStringMass(&aString); 121 122 #ifdef debug_QGSMfragmentation 123 G4cout<<G4endl<<"QGSM StringFragm: String Ma 124 <<theString.Get4M 125 <<theString.Get4M 126 <<"-------------- 127 G4cout<<"String ends Direct "<<theString.Get 128 <<theString.Get 129 <<theString.Get 130 G4cout<<"Left mom "<<theString.GetLeftParto 131 G4cout<<"Right mom "<<theString.GetRightPart 132 G4cout<<"Check for Fragmentation "<<G4endl; 133 #endif 134 135 // Can no longer modify Parameters for Fragm 136 PastInitPhase=true; 137 138 // Check if string has enough mass to fragme 139 G4KineticTrackVector * LeftVector=NULL; 140 141 if ( !IsItFragmentable(&aString) ) { 142 LeftVector=ProduceOneHadron(&theString); 143 144 #ifdef debug_QGSMfragmentation 145 if ( LeftVector != 0 ) G4cout<<"Non fragm 146 #endif 147 148 if ( LeftVector == nullptr ) LeftVector = 149 return LeftVector; 150 } 151 152 #ifdef debug_QGSMfragmentation 153 G4cout<<"The string will be fragmented. "<<G 154 #endif 155 156 LeftVector = new G4KineticTrackVector; 157 G4KineticTrackVector * RightVector=new G4Kin 158 159 G4ExcitedString *theStringInCMS=CopyExcited( 160 G4LorentzRotation toCms=theStringInCMS->Tran 161 162 G4bool success=false, inner_sucess=true; 163 G4int attempt=0; 164 while ( !success && attempt++ < StringLoopIn 165 { 166 #ifdef debug_QGSMfragmentation 167 G4cout<<"Loop_toFrag "<<theStr 168 <<theStr 169 <<theStr 170 #endif 171 172 G4FragmentingString *currentString=new G4F 173 174 std::for_each(LeftVector->begin(), LeftVec 175 LeftVector->clear(); 176 std::for_each(RightVector->begin(), RightV 177 RightVector->clear(); 178 179 inner_sucess=true; // set false on failur 180 const G4int maxNumberOfLoops = 181 G4int loopCounter = -1; 182 while (! StopFragmenting(currentString) && 183 { // Split current string into hadron + n 184 185 #ifdef debug_QGSMfragm 186 G4cout<<"The string ca 187 #endif 188 G4FragmentingString *newString=0; // us 189 G4KineticTrack * Hadron=Splitup(currentS 190 191 if ( Hadron != 0 ) 192 { 193 #ifdef debug_QGSMfr 194 G4cout<<"Hadron pro 195 #endif 196 // To close the pro 197 if ( currentString->GetDecayDirection 198 LeftVector->push_back(Hadron); 199 else 200 RightVector->push_back(Hadron); 201 202 delete currentString; 203 currentString=newString; 204 205 } else { 206 207 #ifdef debug_QGSMfr 208 G4cout<<"abandon .. 209 #endif 210 211 // Abandon ... start from the beginni 212 if (newString) delete newString; 213 inner_sucess=false; 214 break; 215 } 216 } 217 if ( loopCounter >= maxNumberO 218 inner_sucess=false; 219 } 220 221 // Split current string into 2 final Hadro 222 #ifdef debug_QGSMfragmentation 223 if( inner_sucess ) { 224 G4cout<<"Split remaining str 225 } else { 226 G4cout<<" New attempt to fragment string 227 } 228 #endif 229 // To the close production of 230 if ( inner_sucess && 231 SplitLast(currentString,LeftVector, R 232 { 233 success=true; 234 } 235 delete currentString; 236 } 237 238 delete theStringInCMS; 239 240 if ( ! success ) 241 { 242 std::for_each(LeftVector->begin(), LeftVec 243 LeftVector->clear(); 244 std::for_each(RightVector->begin(), RightV 245 delete RightVector; 246 return LeftVector; 247 } 248 249 // Join Left- and RightVector into LeftVecto 250 while(!RightVector->empty()) /* Loop checki 251 { 252 LeftVector->push_back(RightVector->back( 253 RightVector->erase(RightVector->end()-1) 254 } 255 delete RightVector; 256 257 CalculateHadronTimePosition(theString.Get4Mo 258 259 G4LorentzRotation toObserverFrame(toCms.inve 260 261 for (size_t C1 = 0; C1 < LeftVector->size(); 262 { 263 G4KineticTrack* Hadron = LeftVector->oper 264 G4LorentzVector Momentum = Hadron->Get4Mo 265 Momentum = toObserverFrame*Momentum; 266 Hadron->Set4Momentum(Momentum); 267 G4LorentzVector Coordinate(Hadron->GetPos 268 Momentum = toObserverFrame*Coordinate; 269 Hadron->SetFormationTime(Momentum.e()); 270 G4ThreeVector aPosition(Momentum.vect()); 271 Hadron->SetPosition(theString.GetPosition 272 } 273 return LeftVector; 274 } 275 276 //-------------------------------------------- 277 278 G4bool G4QGSMFragmentation::IsItFragmentable(c 279 { 280 return sqr( MinimalStringMass + MassCut ) < 281 } 282 283 //-------------------------------------------- 284 285 G4bool G4QGSMFragmentation::StopFragmenting(co 286 { 287 SetMinimalStringMass(string); 288 if ( MinimalStringMass < 0.0 ) return 289 290 G4double smass = string->Mass(); 291 G4double x = (string->IsAFourQuarkString()) 292 : 0.66e-6*(smass - MinimalStringMass)*(sma 293 294 G4bool res = true; 295 if(x > 0.0) { 296 res = (x < 200.) ? (G4UniformRand() 297 } 298 return res; 299 } 300 301 //-------------------------------------------- 302 303 G4KineticTrack * G4QGSMFragmentation::Splitup( 304 G4FragmentingStri 305 { 306 #ifdef debug_QGSMfragmentation 307 G4cout<<G4endl; 308 G4cout<<"Start SplitUP (G4VLongitudinal 309 G4cout<<"String partons: " <<string->Ge 310 <<string->Ge 311 <<"Direction " <<string->Ge 312 #endif 313 314 //... random choice of string end to us 315 G4int SideOfDecay = (G4UniformRand() < 316 if (SideOfDecay < 0) 317 { 318 string->SetLeftPartonStable(); 319 } else 320 { 321 string->SetRightPartonStable(); 322 } 323 324 G4ParticleDefinition *newStringEnd; 325 G4ParticleDefinition * HadronDefinition 326 if (string->DecayIsQuark()) 327 { 328 G4double ProbDqADq = GetDiquarkSuppress(); 329 330 G4int NumberOfpossibleBaryons = 2; 331 332 if (string->GetLeftParton()->GetParticleSu 333 if (string->GetRightParton()->GetParticleS 334 335 G4double ActualProb = ProbDqADq ; 336 ActualProb *= (1.0-G4Exp(2.0*(1.0 - 337 338 SetDiquarkSuppression(ActualProb); 339 340 HadronDefinition= QuarkSplitup(strin 341 342 SetDiquarkSuppression(ProbDqADq); 343 } else { 344 HadronDefinition= DiQuarkSplitup(str 345 } 346 347 if ( HadronDefinition == NULL ) return 348 349 #ifdef debug_QGSMfragmentation 350 G4cout<<"The parton "<<string->GetDecay 351 <<" produces hadron "<<HadronDefi 352 <<" and is transformed to "<<newS 353 G4cout<<"The side of the string decay L 354 #endif 355 // create new String from old, ie. keep 356 357 newString=new G4FragmentingString(*stri 358 359 360 #ifdef debug_QGSMfragmentation 361 G4cout<<"An attempt to determine its en 362 #endif 363 G4LorentzVector* HadronMomentum=SplitEa 364 365 delete newString; newString=0; 366 367 G4KineticTrack * Hadron =0; 368 if ( HadronMomentum != 0 ) { 369 370 #ifdef debug_QGSMfragmentation 371 G4cout<<"The attempt was successful 372 #endif 373 G4ThreeVector Pos; 374 Hadron = new G4KineticTrack(HadronDefinit 375 376 newString=new G4FragmentingString(*st 377 378 delete HadronMomentum; 379 } 380 else 381 { 382 #ifdef debug_QGSMfragmentation 383 G4cout<<"The attempt was not successf 384 #endif 385 } 386 387 #ifdef debug_VStringDecay 388 G4cout<<"End SplitUP (G4VLongitudinalSt 389 #endif 390 391 return Hadron; 392 } 393 394 //-------------------------------------------- 395 396 G4ParticleDefinition *G4QGSMFragmentation::DiQ 397 398 { 399 //... can Diquark break or not? 400 if (G4UniformRand() < DiquarkBreakProb ) / 401 { 402 G4int stableQuarkEncoding = decay->GetPD 403 G4int decayQuarkEncoding = (decay->GetPD 404 405 if (G4UniformRand() < 0.5) 406 { 407 G4int Swap = stableQuarkEncoding; 408 stableQuarkEncoding = decayQuarkEncod 409 decayQuarkEncoding = Swap; 410 } 411 412 G4int IsParticle=(decayQuarkEncoding>0) 413 414 G4double StrSup=GetStrangeSuppress(); 415 SetStrangenessSuppression((1.0 - 0.07)/2 416 pDefPair QuarkPair = CreatePartonPair(Is 417 SetStrangenessSuppression(StrSup); 418 419 //... Build new Diquark 420 G4int QuarkEncoding=QuarkPair.second->Ge 421 G4int i10 = std::max(std::abs(QuarkEnco 422 G4int i20 = std::min(std::abs(QuarkEnco 423 G4int spin = (i10 != i20 && G4UniformRan 424 G4int NewDecayEncoding = -1*IsParticle*( 425 426 created = FindParticle(NewDecayEncoding) 427 G4ParticleDefinition * decayQuark=FindPa 428 G4ParticleDefinition * had=hadronizer->B 429 430 DecayQuark = decay->GetPDGEncoding(); 431 NewQuark = NewDecayEncoding; 432 433 return had; 434 435 } else { //... Diquark does not break 436 437 G4int IsParticle=(decay->GetPDGEncoding( 438 G4double StrSup=GetStrangeSuppress(); / 439 SetStrangenessSuppression((1.0 - 0.07)/2 440 pDefPair QuarkPair = CreatePartonPair(Is 441 SetStrangenessSuppression(StrSup); 442 443 created = QuarkPair.second; 444 445 DecayQuark = decay->GetPDGEncoding(); 446 NewQuark = created->GetPDGEncoding(); 447 448 G4ParticleDefinition * had=hadronizer->B 449 return had; 450 } 451 } 452 453 //-------------------------------------------- 454 455 G4LorentzVector * G4QGSMFragmentation::SplitEa 456 457 458 { 459 G4double HadronMass = pHadron->GetPDGMa 460 461 SetMinimalStringMass(NewString); 462 463 if ( MinimalStringMass < 0.0 ) return n 464 465 #ifdef debug_QGSMfragmentation 466 G4cout<<"G4QGSMFragmentation::SplitEand 467 G4cout<<"String 4 mom, String M "<<stri 468 G4cout<<"HadM MinimalStringMassLeft Str 469 <<string->Mass()<<" "<<HadronMass 470 #endif 471 472 if (HadronMass + MinimalStringMass > st 473 { 474 #ifdef debug_QGSMfragmentation 475 G4cout<<"Mass of the string is not su 476 #endif 477 return 0; 478 } // have to start all over! 479 480 // calculate and assign hadron transver 481 G4double StringMT2 = string->MassT2(); 482 G4double StringMT = std::sqrt(StringMT 483 484 G4LorentzVector String4Momentum = strin 485 String4Momentum.setPz(0.); 486 G4ThreeVector StringPt = String4Momentu 487 488 G4ThreeVector HadronPt , RemSysPt; 489 G4double HadronMassT2, ResidualMas 490 491 // Mt distribution is implemented 492 G4double HadronMt, Pt, Pt2, phi; 493 494 //... sample Pt of the hadron 495 G4int attempt=0; 496 do 497 { 498 attempt++; if (attempt > StringLoopIn 499 500 HadronMt = HadronMass - 200.0*G4Log(G 501 Pt2 = sqr(HadronMt)-sqr(HadronMass); 502 phi = 2.*pi*G4UniformRand(); 503 G4ThreeVector SampleQuarkPtw= G4Three 504 HadronPt =SampleQuarkPtw + string->D 505 HadronPt.setZ(0); 506 RemSysPt = StringPt - HadronPt; 507 508 HadronMassT2 = sqr(HadronMass) + Hadr 509 ResidualMassT2=sqr(MinimalStringMass) 510 511 } while (std::sqrt(HadronMassT2) + std: 512 513 //... sample z to define hadron longit 514 //... but first check the available pha 515 516 G4double Pz2 = (sqr(StringMT2 - HadronM 517 4*HadronMassT2 * ResidualMassT2)/4./ 518 519 if ( Pz2 < 0 ) {return 0;} // 520 521 //... then compute allowed z region z_ 522 523 G4double Pz = std::sqrt(Pz2); 524 G4double zMin = (std::sqrt(HadronMassT2 525 G4double zMax = (std::sqrt(HadronMassT2 526 527 if (zMin >= zMax) return 0; // have 528 529 G4double z = GetLightConeZ(zMin, zMax, 530 string->GetDecayParton() 531 HadronPt.x(), HadronPt.y 532 533 //... now compute hadron longitudinal m 534 // longitudinal hadron momentum compone 535 536 HadronPt.setZ( 0.5* string->GetDecayDir 537 (z * string->LightConeDecay() - 538 HadronMassT2/(z * string->LightCone 539 G4double HadronE = 0.5* (z * string->L 540 HadronMassT2/(z * string->LightConeDe 541 542 G4LorentzVector * a4Momentum= new G4Lor 543 544 #ifdef debug_QGSMfragmentation 545 G4cout<<"string->GetDecayDirection() st 546 <<string->GetDecayDirection()<<" 547 G4cout<<"HadronPt,HadronE "<<HadronPt<< 548 G4cout<<"Out of QGSM SplitEandP "<<G4en 549 #endif 550 551 return a4Momentum; 552 } 553 554 //-------------------------------------------- 555 556 G4double G4QGSMFragmentation::GetLightConeZ(G4 557 G4 558 { 559 G4double lambda = 2.0*(sqr(ptx)+sqr(pty))/sq 560 561 #ifdef debug_QGSMfragmentation 562 G4cout<<"GetLightConeZ zmin zmax Parton pHad 563 <<" "/*<< pHadron->GetParticleName() * 564 #endif 565 566 G4double z(0.); 567 G4int DiQold(0), DiQnew(0); 568 G4double d1(-1.0), d2(0.); 569 G4double invD1(0.),invD2(0.), r1(0.),r2(0.), 570 571 G4int absDecayQuarkCode = std::abs( DecayQua 572 G4int absNewQuarkCode = std::abs( NewQuark 573 574 G4int q1(0), q2(0); 575 // q1 = absDecayQuarkCode/1000; q2 = (absDe 576 577 G4int qA(0), qB(0); 578 // qA = absNewQuarkCode/1000; qB = (absNe 579 580 if ( (absDecayQuarkCode < 6) && (absNewQuark 581 d1 = FFq2q[absDecayQuarkCode-1][absNewQuar 582 } 583 584 if ( (absDecayQuarkCode < 6) && (absNewQuark 585 qA = absNewQuarkCode/1000; qB = (absNewQu 586 d1 = FFq2qq[absDecayQuarkCode-1][DiQnew][0] 587 } 588 589 if ( (absDecayQuarkCode > 6) && (absNewQuark 590 q1 = absDecayQuarkCode/1000; q2 = (absDeca 591 d1 = FFqq2q[DiQold][absNewQuarkCode-1][0]; 592 } 593 594 if ( d1 < 0. ) { 595 q1 = absDecayQuarkCode/1000; q2 = (absDeca 596 qA = absNewQuarkCode/1000; qB = (absNewQ 597 d1 = FFqq2qq[DiQold][DiQnew][0]; d2 = FFqq 598 } 599 600 d2 +=lambda; 601 d1+=1.0; d2+=1.0; 602 603 invD1=1./d1; invD2=1./d2; 604 605 const G4int maxNumberOfLoops = 10000; 606 G4int loopCounter = 0; 607 do // Jong's algorithm 608 { 609 r1=G4Pow::GetInstance()->powA(G4UniformRan 610 r2=G4Pow::GetInstance()->powA(G4UniformRan 611 r12=r1+r2; 612 z=r1/r12; 613 } while ( ( (r12 > 1.0) || !((zmin <= z)&&(z 614 ++loopCounter < maxNumberOfLoops ) 615 616 if ( loopCounter >= maxNumberOfLoops ) { 617 z = 0.5*(zmin + zmax); // Just a value be 618 } 619 620 return z; 621 } 622 623 //-------------------------------------------- 624 625 G4bool G4QGSMFragmentation::SplitLast(G4Fragme 626 G4KineticTrackVector * Lef 627 G4KineticTrackVector * Right 628 { 629 //... perform last cluster decay 630 631 G4ThreeVector ClusterVel =string->Get4Mome 632 G4double ResidualMass =string->Mass(); 633 634 #ifdef debug_QGSMfragmentation 635 G4cout<<"Split last----------------------- 636 G4cout<<"StrMass "<<ResidualMass<<" q's " 637 <<string->GetLeftParton()->GetPartic 638 <<string->GetRightParton()->GetParti 639 #endif 640 641 G4int cClusterInterrupt = 0; 642 G4ParticleDefinition *LeftHadron = nullptr 643 G4ParticleDefinition *RightHadron = nullpt 644 const G4int maxNumberOfLoops = 1000; 645 G4int loopCounter = 0; 646 647 G4double LeftHadronMass(0.); G4double Righ 648 do 649 { 650 if (cClusterInterrupt++ >= ClusterLoop 651 LeftHadronMass = -MaxMass; RightHadron 652 653 G4ParticleDefinition * quark = nullptr; 654 string->SetLeftPartonStable(); // to query q 655 656 if (string->DecayIsQuark() && string->Stable 657 { 658 //... there are quarks on cluster ends 659 660 G4int IsParticle=(string->GetLeftParton()- 661 // if we have a quark, we need 662 663 pDefPair QuarkPair = CreatePartonPair(IsPa 664 quark = QuarkPair.second; 665 666 LeftHadron= hadronizer->Build(QuarkPair.fi 667 if ( LeftHadron == NULL ) continue; 668 RightHadron = hadronizer->Build(stri 669 if ( RightHadron == NULL ) continue; 670 } else if( (!string->DecayIsQuark() && stri 671 ( string->DecayIsQuark() && !stri 672 //... there is a Diquark on one of cluster 673 G4int IsParticle; 674 if ( string->StableIsQuark() ) { 675 IsParticle=(string->GetLeftParton()->Get 676 } else { 677 IsParticle=(string->GetLeftParton()->Get 678 } 679 680 //G4double ProbSaS = 1.0 - 2.0 * G 681 //G4double ActualProb = ProbSaS * 1. 682 //SetStrangenessSuppression((1.0-Act 683 684 pDefPair QuarkPair = CreatePartonPai 685 //SetStrangenessSuppression((1.0-Pro 686 quark = QuarkPair.second; 687 LeftHadron=hadronizer->Build(QuarkPa 688 if ( LeftHadron == NULL ) continue; 689 RightHadron = hadronizer->Build(stri 690 if ( RightHadron == NULL ) continue; 691 } else { // Diquark and anti-diquark 692 if (cClusterInterrupt++ >= ClusterLo 693 G4int LeftQuark1= string->GetLeftPar 694 G4int LeftQuark2=(string->GetLeftPar 695 G4int RightQuark1= string->GetRightP 696 G4int RightQuark2=(string->GetRightP 697 if (G4UniformRand()<0.5) { 698 LeftHadron =hadronizer->Build(Fin 699 RightHadron =hadronizer->Build(Fin 700 } else { 701 LeftHadron =hadronizer->Build(Fin 702 RightHadron =hadronizer->Build(Fin 703 } 704 if ( (LeftHadron == NULL) || (RightHadron 705 } 706 LeftHadronMass = LeftHadron->GetPDGMa 707 RightHadronMass = RightHadron->GetPDGM 708 //... repeat procedure, if mass of clu 709 } while ( ( ResidualMass <= LeftHadronMass 710 && ++loopCounter < maxNumberOfLo 711 712 if ( loopCounter >= maxNumberOfLoops ) { 713 return false; 714 } 715 716 //... compute hadron momenta and energies 717 G4LorentzVector LeftMom, RightMom; 718 G4ThreeVector Pos; 719 Sample4Momentum(&LeftMom , LeftHadron->Get 720 &RightMom, RightHadron->Ge 721 LeftMom.boost(ClusterVel); 722 RightMom.boost(ClusterVel); 723 724 #ifdef debug_QGSMfragmentation 725 G4cout<<LeftHadron->GetParticleName()<<" " 726 G4cout<<"Left Hadrom P M "<<LeftMom<<" "< 727 G4cout<<"Right Hadrom P M "<<RightMom<<" " 728 #endif 729 730 LeftVector->push_back(new G4KineticTrack(L 731 RightVector->push_back(new G4KineticTrack( 732 733 return true; 734 } 735 736 //-------------------------------------------- 737 738 void G4QGSMFragmentation::Sample4Momentum(G4Lo 739 G4Lo 740 { 741 #ifdef debug_QGSMfragmentation 742 G4cout<<"Sample4Momentum Last------------- 743 G4cout<<" StrMass "<<InitialMass<<" Mass1 744 G4cout<<" SumMass "<<Mass+AntiMass<<G4end 745 #endif 746 747 G4double r_val = sqr(InitialMass*InitialMa 748 G4double Pabs = (r_val > 0.)? std::sqrt(r_ 749 750 //... sample unit vector 751 G4double pz = 1. - 2.*G4UniformRand(); 752 G4double st = std::sqrt(1. - pz * pz)* 753 G4double phi = 2.*pi*G4UniformRand(); 754 G4double px = st*std::cos(phi); 755 G4double py = st*std::sin(phi); 756 pz *= Pabs; 757 758 Mom->setPx(px); Mom->setPy(py); Mom->setPz 759 Mom->setE(std::sqrt(Pabs*Pabs + Mass*Mass) 760 761 AntiMom->setPx(-px); AntiMom->setPy(-py); 762 AntiMom->setE (std::sqrt(Pabs*Pabs + AntiM 763 } 764 765 //-------------------------------------------- 766 767 void G4QGSMFragmentation::SetFFq2q() // q-> q 768 { 769 for (G4int i=0; i < 5; i++) { 770 FFq2q[i][0][0] = 2.0 ; FFq2q[i][0][1] = -a 771 FFq2q[i][1][0] = 2.0 ; FFq2q[i][1][1] = -a 772 FFq2q[i][2][0] = 1.0 ; FFq2q[i][2][1] = -a 773 FFq2q[i][3][0] = 1.0 ; FFq2q[i][3][1] = -a 774 FFq2q[i][4][0] = 1.0 ; FFq2q[i][4][1] = -a 775 } 776 } 777 778 //-------------------------------------------- 779 780 void G4QGSMFragmentation::SetFFq2qq() // q-> 781 { 782 for (G4int i=0; i < 5; i++) { 783 FFq2qq[i][ 0][0] = 0.0 ; FFq2qq[i][ 0][1] 784 FFq2qq[i][ 1][0] = 0.0 ; FFq2qq[i][ 1][1] 785 FFq2qq[i][ 2][0] = 0.0 ; FFq2qq[i][ 2][1] 786 FFq2qq[i][ 3][0] = 0.0 ; FFq2qq[i][ 3][1] 787 FFq2qq[i][ 4][0] = 0.0 ; FFq2qq[i][ 4][1] 788 FFq2qq[i][ 5][0] = 0.0 ; FFq2qq[i][ 5][1] 789 FFq2qq[i][ 6][0] = 0.0 ; FFq2qq[i][ 6][1] 790 FFq2qq[i][ 7][0] = 0.0 ; FFq2qq[i][ 7][1] 791 FFq2qq[i][ 8][0] = 0.0 ; FFq2qq[i][ 8][1] 792 FFq2qq[i][ 9][0] = 0.0 ; FFq2qq[i][ 9][1] 793 FFq2qq[i][10][0] = 0.0 ; FFq2qq[i][10][1] 794 FFq2qq[i][11][0] = 0.0 ; FFq2qq[i][11][1] 795 FFq2qq[i][12][0] = 0.0 ; FFq2qq[i][12][1] 796 FFq2qq[i][13][0] = 0.0 ; FFq2qq[i][13][1] 797 FFq2qq[i][14][0] = 0.0 ; FFq2qq[i][14][1] 798 } 799 } 800 801 //-------------------------------------------- 802 803 void G4QGSMFragmentation::SetFFqq2q() // q1q2 804 { 805 for (G4int i=0; i < 15; i++) { 806 FFqq2q[i][0][0] = 2.0*(arho - an); FFqq2q[ 807 FFqq2q[i][1][0] = 2.0*(arho - an); FFqq2q[ 808 FFqq2q[i][2][0] = 2.0*(arho - an); FFqq2q[ 809 FFqq2q[i][3][0] = 2.0*(arho - an); FFqq2q[ 810 FFqq2q[i][4][0] = 2.0*(arho - an); FFqq2q[ 811 } 812 } 813 814 //-------------------------------------------- 815 816 void G4QGSMFragmentation::SetFFqq2qq() // q1( 817 { 818 for (G4int i=0; i < 15; i++) { 819 FFqq2qq[i][0][0] = 0. ; FFqq2qq[i][0][1] 820 FFqq2qq[i][1][0] = 0. ; FFqq2qq[i][1][1] 821 FFqq2qq[i][2][0] = 0. ; FFqq2qq[i][2][1] 822 FFqq2qq[i][3][0] = 0. ; FFqq2qq[i][3][1] 823 FFqq2qq[i][4][0] = 0. ; FFqq2qq[i][4][1] 824 } 825 } 826 827